Probing Picohertz Gravitational Waves with Pulsars

TL;DR

This work tackles the gap in gravitational-wave searches at picohertz frequencies by exploiting secular pulsar timing parameters, using a Bayesian framework to search for continuous pHz GWs from supermassive black hole binaries. It derives analytic expressions for GW-induced drifts in $\dot{P}_b$ and $\ddot{P}$, extends SMBHB population models to the pHz regime, and applies the method to an expanded pulsar dataset to set world-leading upper limits on the CW strain $h_0$. Although no detection emerges yet, the study demonstrates an order-of-magnitude improvement over previous limits and outlines clear SKA-era prospects, including the potential to observe early SMBHB evolution and to probe new physics in the early universe. The results establish a robust framework for pHz GW searches that can leverage future pulsar discoveries, precise astrometry, and targeted priors to turn upper limits into potential detections and population-level constraints.

Abstract

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.

Figures (15)

• Figure 1: Pulsar timing parameters capture the metric perturbations due to pHz GWs. LHS: the pulsars (red) and the Earth (blue) are both affected by the passing GW. RHS: short observation baselines ($T_{\mathrm{obs}}$) allow for treating the perturbations as stationary in time. Pulsars at different locations are probing the CWs emitted by the source at different times. For a 10 pHz CW, the period is $\sim3,000$ years, while the observation time is $\sim10$ years. The metric perturbations manifest as measurable redshift $v_{\mathrm{GW}}$ and the associated derivatives $a_{\mathrm{GW}}$ (acceleration) and $j_{\mathrm{GW}}$ (jerk). The explicit expressions of these quantities in terms of the CW waveform are given in Appendix \ref{['sec:analytic derivation']}.
• Figure 2: CW strain amplitude upper limits and forecast detection sensitivities of pulsar timing parameters $\dot{P}_b$ and $\ddot{P}$. The fiducial CW sources are modeled as circular orbit SMBHBs. The current upper limit in this work is an order of magnitude better than best constraints from previous work in this frequency range DeRoccoDror2024. Forecast sensitivities with SKA data provide another order of magnitude improvement by increase of data size alone, without accounting for the improvement in data uncertainty. Current results are comparable to PTA upper limit on CWs (purple dashed line) for certain frequencies. Expected individual CW sources (denoted by the stars) are at least 3 orders of magnitude in strain amplitude below the current constraint according to our SMBHB population model. The orange solid line is the GWB spectrum drawn from our population model by matching self-consistently with GWB at $f_{\mathrm{GW}}=1~\mathrm{yr}^{-1}$ detected by NANOGrav. The shaded regions centered on this curve denotes the uncertainty quantiles. The black dashed line represents the power-law GWB spectrum matched to NANOGrav GWB at $f_{\mathrm{GW}}=1~\mathrm{yr}^{-1}$.
• Figure 3: Null distributions of the $\dot{P}_b$ and $\ddot{P}$ data used to place current upper bounds on CW strain amplitude. Both have a small tail surpassing BF=1. In each case, a null CW signal is injected to the simulated data, and $1,000$ realizations of BFs are extracted.
• Figure 4: $\dot{P}_b$ search on current data with uniform prior of $f_{\text{GW}}$. Upper limit on $h_0$ can be extracted from the posterior distribution. The priors (blue dashed line) are also shown. The contours are set at $50\%, 80\%, 95\%$. All the angles and phase are in radians, and $f_{\text{GW}}$ is in $\log_{10} \text{Hz}$. The following figures follow the same convention.
• Figure 5: $\ddot{P}$ search on current data with uniform prior of $f_{\text{GW}}$. Upper limit on $h_0$ can be extracted from the posterior distribution. The priors (blue dashed line) are also shown.