Probing Picohertz Gravitational Waves with Pulsars

With periods much longer than the duration of current pulsar timing surveys, gravitational waves in the picohertz (pHz) regime are not detectable in the typical analysis framework for pulsar timing data. However, signatures of these low-frequency signals persist in the slow variation of pulsar timing parameters. In this work, we present the results of the first Bayesian search for continuous pHz gravitational waves using the drift of two sensitive pulsar timing parameters – time derivative of pulsar binary orbital period $\dot{P}_b$ and second order time derivative of pulsar spin period $\ddot{P}$. We apply our new technique to a dataset with more than double the number of pulsars as previous searches in this frequency band, achieving an order-of-magnitude sensitivity improvement. No continuous wave signal is detected in current data; however, we show that future observations by the Square Kilometre Array will provide significantly improved sensitivity and the opportunity to observe continuous pHz signals, including the early stages of supermassive black hole mergers. We explore the detection prospects for this signal by extending existing population models into the pHz regime, finding that future observations will probe phenomenologically-interesting parameter space. Our new Bayesian technique and leading sensitivity in this frequency domain paves the way for new discoveries in both black hole astrophysics and the search for new physics in the early universe.

💡 Research Summary

The paper introduces a novel Bayesian search for continuous gravitational waves (CWs) in the picohertz (pHz) band (10⁻¹⁸–10⁻¹⁶ Hz) using secular drifts in two pulsar timing parameters: the binary orbital period derivative (˙P_b) and the second derivative of the spin period (¨P). Traditional pulsar timing arrays (PTAs) are insensitive to such ultra‑low frequencies because the lowest detectable frequency is set by the inverse of the observation span; pHz waves have periods of thousands of years, so only a tiny segment of the waveform is sampled. Nevertheless, a passing pHz GW induces a slowly varying redshift, an effective velocity, and consequently a constant acceleration (a_GW) and jerk (j_GW) at the Earth–pulsar system. These manifest as additional contributions to ˙P_b and ¨P that can be modeled analytically.

The authors develop an analytic model linking the GW strain amplitude h₀, frequency f_GW, sky location (α, δ), inclination i, polarization angle ϕ, and initial orbital phase Ψ₀ to the observable accelerations and jerks. They construct a Gaussian likelihood assuming uncorrelated pulsar uncertainties, and adopt uniform priors for angular parameters while using a log‑uniform prior for h₀ (switching to a linear‑uniform prior when the data are non‑informative). A Markov Chain Monte Carlo (MCMC) sampler explores the seven‑dimensional posterior, and detection significance is quantified via Bayes factors (BF) computed with the Savage‑Dickey density ratio. A detection is claimed only if the BF exceeds the null distribution threshold and the marginalized posterior for h₀ excludes zero.

To assess detection prospects, the authors extend an existing supermassive black‑hole binary (SMBHB) population synthesis model (originally calibrated to the NANOGrav 15‑yr stochastic background) down to pHz frequencies. The model incorporates dual‑AGN counts, a mass‑ratio distribution, dynamical friction, and stellar hardening, but deliberately omits gas‑driven hardening, which remains uncertain. By anchoring the pHz SMBHB population to the nHz background, they generate realistic CW strain distributions for the pHz band.

The data set comprises 30 millisecond pulsars (MSPs) with measured ˙P_b (an increase from the 14 used in prior work) and 29 MSPs with ¨P measurements. Galactic accelerations and Shklovskii effects are subtracted using the MWPotential2014 model via the galpy package. Pulsars whose distance uncertainties are smaller than a GW wavelength (~100 pc) are retained, reducing the sample to roughly half for the ¨P analysis. Observation baselines range from 4 to 22 years with µs‑level timing residuals.

Applying the Bayesian pipeline to the current data yields no significant CW detection; the Bayes factor falls below the detection threshold, and the 95 % upper limit on h₀ reaches ≈10⁻²⁰, an order‑of‑magnitude improvement over previous pHz searches. This demonstrates that the drift‑parameter method substantially enhances sensitivity in a frequency regime previously inaccessible to PTAs.

Looking forward, the authors project that the Square Kilometre Array (SKA) will dramatically increase the number of precisely timed MSPs and extend observation baselines, potentially improving strain sensitivity by one to two orders of magnitude. Such gains would enable the detection of individual pHz CW sources, notably SMBHBs with total masses 10⁸–10¹⁰ M⊙ at separations ≲1 pc, whose GW emission peaks in the pHz band. Moreover, the method opens a window onto exotic sources such as cosmic‑string networks, first‑order phase transitions, and primordial inflationary or pre‑heating signals that could populate this ultra‑low‑frequency band.

In summary, the paper delivers the first Bayesian framework that leverages secular drifts in pulsar timing parameters to search for continuous pHz gravitational waves, achieves a tenfold sensitivity gain over earlier attempts, and outlines a clear path—via next‑generation radio telescopes—to probe a previously hidden segment of the gravitational‑wave spectrum, with profound implications for supermassive black‑hole astrophysics and early‑Universe physics.