Probing Picohertz Gravitational Waves with Pulsars

TL;DR

This paper introduces a Bayesian method to detect ultra-low-frequency gravitational waves in pulsar timing data, achieving improved sensitivity and exploring future detection prospects with upcoming telescopes like the SKA.

Contribution

The authors develop a novel Bayesian approach to search for picohertz gravitational waves using pulsar timing parameters, significantly enhancing sensitivity over previous methods.

Findings

No signals detected in current data.

Future SKA observations will improve detection prospects.

Extended population models suggest promising parameter space coverage.

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.

Tables

Table 1. Table 1: Prior distributions for model parameters.

Parameter	Prior Distribution	Domain
$\alpha$	uniform
$\delta$	$\frac{1}{2}\cos \left(b\right)$
$i$	$\frac{1}{2}\cos \left(i\right)$
$\varphi$	uniform
${\Psi }_0$	uniform
${\log }_{10}{f}_{\mathrm{GW}}$	uniform
${\log }_{10}{h}_0$	uniform

Table 2. Table 2: Complete table of pulsar data including Galactic coordinates, distance, and fractional acceleration in orbital period Δ ​ P ˙ b / P b \Delta\dot{P}_{b}/P_{b} . This quantity includes all contributions that should sum up to a GW a_{\text{GW}} as in Eq.  3 . Uncertainties are reported as 1 standard deviation. The contribution from Galactic acceleration is calculated using galpy  [ Bovy2015 ] , and all distances are measured using parallax. The distance and timing parameter data can be found in  [ Moran:2023 , DeRoccoDror2024 ] . The pulsars satisfying the distance measurement criterion discussed in Sec.  IV are in bold face.

Pulsar	$l$ (deg)	$b$ (deg)	d (kpc)	$\Delta {\dot{P}}_b/P_b\left({10}^{-18}{s}^{-1}\right)$
B1259-63	304.184	-0.992	2.6(4)	130(70)
B1534+12	19.848	48.341	0.94(7)	0.2(1)
B1913+16	49.968	2.122	4(1)	-1.26(8)
J0348+0432	183.337	-36.774	1.6(9)	-2(6)
J0437-4715	253.394	-41.963	0.16(1)	0.05(6)
J0613-0200	210.413	-9.305	1.11(5)	-0.1(2)
J0737-3039A/B	245.236	-4.505	0.74(6)	-0.5(40)
J0740+6620	149.730	29.599	1.2(2)	-0.04(60)
J0751+1807	202.730	21.086	1.4(4)	-0.3(10)
J1012+5307	160.347	50.858	0.85(1)	0.2(2)
J1017-7156	291.558	-12.553	1(2)	0.5(6)
J1022+1001	231.795	51.101	0.85(6)	-0.08(11)
J1125-6014	292.504	0.894	1.5(8)	-0.02(60)
J1455-3330	330.722	22.562	0.76(6)	0.6(3)
J1537+1155	19.850	48.340	1.2(2)	-0.09(40)
J1600-3053	344.090	16.451	1.39(4)	0.18(8)
J1603-7202	316.630	-14.496	0.5(3)	0.29(8)
J1614-2230	352.636	20.192	0.69(5)	0.3(2)
J1640+2224	41.051	38.271	1.0(2)	0.4(2)
J1713+0747	28.751	25.223	1.31(4)	-0.03(2)
J1738+0333	27.721	17.742	1.7(1)	0.03(30)
J1741+1351	37.885	21.641	1.8(4)	23(6)
J1909-3744	359.731	-19.596	1.14(1)	0.05(3)
J1933-6211	334.431	-28.632	1.2(4)	-0.3(1)
J2043+1711	61.919	-15.313	1.6(2)	107(60)
J2129-5721	338.005	-43.570	0.6(4)	2.5(2)
J2145-0750	47.777	-42.084	0.62(2)	0.11(5)
J2222-0137	62.018	-46.075	0.2681(1)	0.00(5)
J2234+0611	72.991	-43.006	0.95(4)	-0.4(9)
J2339-0533	81.349	-62.476	1.4(3)	8(5)

Table 3. Table 3: Pulsar timing parameters including Galactic coordinates, distance, fractional spin-down acceleration P ¨ obs / P \ddot{P}_{\text{obs}}/P as in Eq.  4 , variance induced by red noise σ RN \sigma_{\text{RN}} , and references. The uncertainties are reported as 1 standard deviation. For pulsars without a parallax measurement, the electron density model YMW 17  [ yao2017new ] is used to estimate the distance. A relative uncertainty of 40 % 40\% is assumed in such cases. References for both the pulsar distances and the timing parameter data are included.

Pulsar	$l$ (deg)	$b$ (deg)	$d$ (kpc)	${\ddot{P}}_{\text{obs}}/P$ (${10}^{-30}{s}^{-2}$)	${\sigma }_{\text{RN}}$ (${10}^{-30}{s}^{-2}$)	Ref.
J0030+0451	113.141	-57.611	0.3(1)[Arzoumanian_2018]	$-4\left(4\right)$	0.54	[Desvignes2016]
J0034-0534	111.492	-68.069	1.3(5)[yao2017new]	$0\left(20\right)$	0	[Desvignes2016]
J0218+4232	139.508	-17.527	3(1)[Verbiest_2014]	$-2\left(5\right)$	0.14	[Desvignes2016]
J0437-4715	253.394	-41.963	0.16(1)[Moran:2023]	$-1\left(1\right)$	0	[Reardon2016]
J0610-2100	227.747	-18.184	3(1)[yao2017new]	$0\left(50\right)$	0	[Desvignes2016]
J0613-0200	210.413	-9.305	1.11(5)[Moran:2023]	$0.6\left(6\right)$	0.13	[Desvignes2016]
J0621+1002	200.570	-2.013	0.4(2)[yao2017new]	$-70\left(30\right)$	1.28	[Desvignes2016]
J0711-6830	279.531	-23.280	0.11(4)[yao2017new]	$1\left(1\right)$	0	[Reardon2016]
J0751+1807	202.730	21.086	1.4(4)[Moran:2023]	$0\left(2\right)$	0	[Desvignes2016]
J0900-3144	256.162	9.486	0.8(4)[Desvignes2016]	$-10\left(20\right)$	0	[Desvignes2016]
J1012+5307	160.347	50.858	0.85(1)[Moran:2023]	$0.4\left(7\right)$	0.01	[Desvignes2016]
J1022+1001	231.795	51.101	0.85(6)[Moran:2023]	$-2\left(1\right)$	0.01	[Desvignes2016]
J1045-4509	280.851	12.254	0.3(2)[Bobakov_2019]	$-2\left(7\right)$	0.05	[Reardon2016]
J1455-3330	330.722	22.562	0.76(6)[Moran:2023]	$6\left(20\right)$	0.17	[Desvignes2016]
J1600-3053	344.090	16.451	1.39(4)[Moran:2023]	$4\left(5\right)$	0.06	[Reardon2016]
J1603-7202	316.630	-14.496	0.5(3)[Moran:2023]	$1\left(4\right)$	0.02	[Reardon2016]
J1640+2224	41.051	38.271	1.0(2)[Moran:2023]	$-0.9\left(9\right)$	0	[Desvignes2016]
J1643-1224	5.669	21.218	0.8(2)[Desvignes2016]	$-2\left(2\right)$	0.01	[Desvignes2016]
J1713+0747	28.751	25.223	1.31(4)[Moran:2023]	$-0.5\left(5\right)$	0.38	[Desvignes2016]
J1721-2457	8.267	6.751	1.4(6)[yao2017new]	$-30\left(70\right)$	0.01	[Desvignes2016]
J1730-2304	3.137	6.023	0.44(5)[Shamohammadi:2024uvg]	$0\left(2\right)$	0	[Reardon2016]
J1732-5049	340.029	-9.454	1.9(7)[yao2017new]	$20\left(20\right)$	0.03	[Reardon2016]
J1738+0333	27.721	17.742	1.7(1)[yao2017new]	$-30\left(90\right)$	0	[Desvignes2016]
J1744-1134	14.794	9.180	0.38(1)[Shamohammadi:2024uvg]	$0.8\left(8\right)$	0.02	[Desvignes2016]
J1751-2857	0.646	-1.124	1.1(4)[yao2017new]	$-10\left(50\right)$	0	[Desvignes2016]
J1801-1417	14.546	4.162	2(1)[Shamohammadi:2024uvg]	$-30\left(100\right)$	0.02	[Desvignes2016]
J1802-2124	8.382	0.611	0.6(4)[Desvignes2016]	$10\left(60\right)$	0.01	[Desvignes2016]
J1804-2717	3.508	-2.736	0.8(3)[yao2017new]	$-40\left(40\right)$	0	[Desvignes2016]
J1843-1113	22.055	-3.397	5(3)[gaia_2023]	$-7\left(20\right)$	0.05	[Desvignes2016]
J1853+1303	44.875	5.367	0.9(4)[Desvignes2016]	$-30\left(20\right)$	0	[Desvignes2016]
B1855+09	42.290	3.060	0.9(2)[Verbiest_2012]	$1\left(2\right)$	0.03	[Desvignes2016]
J1909-3744	359.731	-19.596	1.14(1)[Moran:2023]	$0.6\left(9\right)$	0.02	[Desvignes2016]
J1910+1256	46.564	-1.795	2.5(5)[K_t_kc__2022]	$30\left(20\right)$	0	[Desvignes2016]
J1911+1347	25.137	-9.579	2.0(3)[K_t_kc__2022]	$14\left(8\right)$	0	[Desvignes2016]
J1911-1114	47.518	1.809	1.1(4)[yao2017new]	$20\left(50\right)$	0	[Desvignes2016]
J1918-0642	30.027	-9.123	1.2(3)[Shamohammadi:2024uvg]	$0\left(8\right)$	2.46	[Desvignes2016]
B1953+29	65.839	0.443	6(3)[yao2017new]	$-20\left(50\right)$	0	[Desvignes2016]
J2010-1323	29.446	-23.540	3(2)[Shamohammadi:2024uvg]	$20\left(20\right)$	0	[Desvignes2016]
J2019+2425	64.746	-6.624	1.2(5)[yao2017new]	$-500\left(900\right)$	0	[Desvignes2016]
J2033+1734	60.857	-13.154	1.7(7)[yao2017new]	$40\left(100\right)$	0	[Desvignes2016]
J2124-3358	10.925	-45.438	0.48(4)[Shamohammadi:2024uvg]	$0\left(3\right)$	0.02	[Reardon2016]
J2129-5721	338.005	-43.570	0.6(4)[Moran:2023]	$-1\left(2\right)$	0	[Reardon2016]
J2145-0750	47.777	-42.084	0.62(2)[Moran:2023]	$-2\left(1\right)$	0.28	[Desvignes2016]
J2229+2643	87.693	-26.284	1.8(7)[yao2017new]	$-20\left(20\right)$	0	[Desvignes2016]
J2317+1439	91.361	-42.360	2.0(4)[Arzoumanian_2018]	$-1\left(3\right)$	0	[Desvignes2016]
J2322+2057	96.515	-37.310	0.8(2)[Shamohammadi:2024uvg]	$30\left(70\right)$	0	[Desvignes2016]

Equations

$$v_{\mathrm{GW}}=\sum _{A=+,\times }{F}_A(\hat{\Omega })\left[h_A(t,0)-h_A(t-d_p,d_p)\right],$$

$$F_A(\hat{\Omega })=\frac{\widehat{d_p^i}\widehat{d_p^j}{\hat{e}}_{ij}^A(\hat{\Omega })}{2(1+\hat{\Omega }\cdot \widehat{d_p})},$$

$$\frac{{\dot{P}}_{b,\mathrm{obs}}}{P_b}=\frac{{\dot{P}}_{b,\text{int}}}{P_b}-a_{\text{shk}}-a_{\text{MW}}-a_{\text{GW}},$$

$$\frac{{\ddot{P}}_{\mathrm{obs}}}{P}=j_{\text{GW}}.$$

$$\mathcal{L}(h_0,f_{\mathrm{GW}},\Theta |{\vec{y}}_{\mathrm{GW}})=\prod _{p=1}^N\frac{1}{\sqrt{2\pi }{\sigma }_p}\exp \left[-\frac{(y_{\mathrm{GW},p}-{\bar{y}}_{\mathrm{GW},p}(h_0,f_{\mathrm{GW}},\Theta ){)}^2}{2{\sigma }_p^2}\right],$$

$$P(\theta |D)=\frac{P(D|\theta )P(\theta )}{P(D)},$$

$$P(\theta )=\pi (\alpha )\pi (\delta )\pi (i)\pi (\varphi )\pi ({\Psi }_0)\pi (f_{\text{GW}})\pi (h_0).$$

$$\frac{d^3{\Phi }_{\mathrm{SMBHB}}}{d{M}_{\mathrm{BH},1}dzdq}=R\frac{d^2{\Phi }_{\mathrm{Q}}}{d{M}_{\mathrm{BH}}dz}p(q),$$

$$v_{GW}=\sum _{A=+,\times }{F}_A(\hat{\Omega })(h_A(t,0)-h_A(t-d_p,\overrightarrow{d_p})=\frac{\widehat{d_p^i}\widehat{d_p^j}}{2(1+\hat{\Omega }\cdot \widehat{d_p})}[e_{ij}^{+}(\hat{\Omega })(h_{+}(t,0)-h_{+}(t-d_p,\overrightarrow{d_p}))+e_{ij}^{\times }(\hat{\Omega })(h_{\times }(t,0)-h_{\times }(t-d_p,\overrightarrow{d_p}))],$$

$$h_0=\frac{2(G\mathcal{M}{)}^{5/3}}{d_L}{\omega }_0^{2/3}$$

$$\{\begin{array}{l}{h}_{+}=h_0\cos (2\varphi )(1+{\cos }^2(i))\sin (2\Psi )\\ +2h_0\sin (2\varphi )\cos (i)\cos (2\Psi )\\ h_{\times }=h_0\sin (2\varphi )(1+{\cos }^2(i))\sin (2\Psi )\\ -2h_0\cos (2\varphi )\cos (i)\cos (2\Psi )\end{array}$$

$$v_{GW}=\frac{\widehat{d_p^i}\widehat{d_p^j}}{2(1+\hat{\Omega }\cdot \widehat{d_p})}[e_{ij}^{+}(\hat{\Omega })(h_0\cos (2\varphi )(1+{\cos }^2(i))\sin (2{\Psi }_0+2{\omega }_{0}t)+2h_0\sin (2\varphi )\cos (i)\cos (2{\Psi }_0+2{\omega }_{0}t)-h_0\cos (2\varphi )(1+{\cos }^2(i))\sin (2{\Psi }_0+2{\omega }_0(t-\hat{\Omega }\cdot \overrightarrow{d_p}))-2h_0\sin (2\varphi )\cos (i)\cos (2{\Psi }_0+2{\omega }_0(t-\hat{\Omega }\cdot \overrightarrow{d_p})))+e_{ij}^{\times }(\hat{\Omega })(h_0\sin (2\varphi )(1+{\cos }^2(i))\sin (2{\Psi }_0+2{\omega }_{0}t)-2h_0\cos (2\varphi )\cos (i)\cos (2{\Psi }_0+2{\omega }_{0}t)-h_0\sin (2\varphi )(1+{\cos }^2(i))\sin (2{\Psi }_0+2{\omega }_0(t-\hat{\Omega }\cdot \overrightarrow{d_p}))+2h_0\cos (2\varphi )\cos (i)\cos (2{\Psi }_0+2{\omega }_0(t-\hat{\Omega }\cdot \overrightarrow{d_p})))]$$

$$\{\begin{array}{l}\hat{\alpha }=-\sin (\alpha )\hat{i}+\cos (\alpha )\hat{j}\\ \hat{\delta }=\cos (\delta )\cos (\alpha )\hat{i}+\cos (\delta )\sin (\alpha )\hat{j}-\sin (\delta )\hat{k}\end{array}$$

$$\{\begin{array}{l}{e}_{+}=\hat{\alpha }\hat{\alpha }-\hat{\delta }\hat{\delta }\\ e_{\times }=\hat{\alpha }\hat{\delta }+\hat{\delta }\hat{\alpha }\end{array}$$

$$a_{GW}=\frac{\widehat{d_p^i}\widehat{d_p^j}}{1+\hat{\Omega }\cdot \widehat{d_p}}{\omega }_0{h}_0[e_{ij}^{+}(\hat{\Omega })[\cos (2\varphi )(1+{\cos }^2(i))(\cos (2{\Psi }_0)-\cos (2{\Psi }_0-2{\omega }_0(1+\hat{\Omega }\cdot \widehat{d_p})d_p))+2\sin (2\varphi )\cos (i)(-\sin (2{\Psi }_0)+\sin (2{\Psi }_0-2{\omega }_0(1+\hat{\Omega }\cdot \widehat{d_p})d_p))]+e_{ij}^{\times }(\hat{\Omega })[\sin (2\varphi )(1+{\cos }^2(i))(\cos (2{\Psi }_0)-\cos (2{\Psi }_0-2{\omega }_0(1+\hat{\Omega }\cdot \widehat{d_p})d_p))-2\cos (2\varphi )\cos (i)(-\sin (2{\Psi }_0)+\sin (2{\Psi }_0-2{\omega }_0(1+\hat{\Omega }\cdot \widehat{d_p})d_p))]]$$

$$j_{GW}=-2\frac{\widehat{d_p^i}\widehat{d_p^j}}{1+\hat{\Omega }\cdot \widehat{d_p}}{\omega }_0^2{h}_0[e_{ij}^{+}(\hat{\Omega })[\cos (2\varphi )(1+{\cos }^2(i))(\sin (2{\Psi }_0)-\sin (2{\Psi }_0-2{\omega }_0(1+\hat{\Omega }\cdot \widehat{d_p})d_p))+2\sin (2\varphi )\cos (i)(\cos (2{\Psi }_0)-\cos (2{\Psi }_0-2{\omega }_0(1+\hat{\Omega }\cdot \widehat{d_p})d_p))]+e_{ij}^{\times }(\hat{\Omega })[\sin (2\varphi )(1+{\cos }^2(i))(\sin (2{\Psi }_0)-\sin (2{\Psi }_0-2{\omega }_0(1+\hat{\Omega }\cdot \widehat{d_p})d_p))-2\cos (2\varphi )\cos (i)(\cos (2{\Psi }_0)-\cos (2{\Psi }_0-2{\omega }_0(1+\hat{\Omega }\cdot \widehat{d_p})d_p))]]$$

$$\{\begin{array}{l}\widehat{d_p}=\sin ({\delta }_p)\cos ({\alpha }_p)\hat{i}+\sin ({\delta }_p)\sin ({\alpha }_p)\hat{j}+\cos ({\delta }_p)\hat{k}\\ \hat{\Omega }=-\sin (\delta )\cos (\alpha )\hat{i}+\sin (\delta )\sin (\alpha )\hat{j}+\cos (\delta )\hat{k}\end{array}$$

Full text

I Introduction

Gravitational waves (GW) have opened an entirely new observational window on the Universe. [maggiore2008gravitational, afzal2023nanograv, 2016gravitational] Over the past decade, multiple detection techniques have dramatically advanced the field [sathyaprakash2021gravitational]. Ground-based interferometers such LIGO, Virgo, and KAGRA operate in the Hz to kHz range, successfully detecting hundreds of stellar-mass compact binary coalescences [Abadie_2010]. Space-based interferometry missions, notably LISA [Amaro-Seoane2017], are designed to explore the mHz regime. Pulsar timing arrays (PTAs) [Antoniadis2022, Chen2021, Goncharov2021, nanograv] operate within the nHz regime and have recently uncovered evidence for the stochastic gravitational wave background (GWB) likely produced by supermassive black hole binaries (SMBHB), and potentially other cosmological sources. At the lowest frequencies, between $10^{-18}$ and $10^{-16}$ Hz, the imprints of primordial GWs can be inferred from polarization measurements in the cosmic microwave background (CMB) [Planck2018Inflation].

This leaves a significant gap in GW frequency band between $10^{-16}$ Hz and 1 nHz, called the picohertz (pHz) regime, where exciting physics may be hiding. First of all, here lies a population of pHz frequency SMBHBs – progenitors of the nHz SMBHB population found in the PTA band [Begelman1980]. All these binaries emit continuous GWs (CW), and a GWB arises from the superposition of all these individual CW sources. In the pHz regime, environmental effects likely dominate over GW emission in the binary evolution, which should influence the SMBHB demographics. Detecting GWs in this pHz frequency regime can teach us about the SMBHB population, and in turn give us insight into the physics driving SMBHB evolution.

A range of beyond-the-Standard-Model and early-Universe processes could also contribute to the GW signals in this band. These include networks of cosmic strings [Ghoshal2023, Chang2020], first-order phase transitions associated with symmetry breaking events in the early Universe [Freese2023, Weir2018], turbulent processes during the quantum chromodynamics epoch [Neronov2021, Brandenburg2021, Moore2021], and primordial GWs generated during inflation [Liddle2000, Kamionkowski2016] and subsequent preheating phases [Boyle2005]. Detection of pHz GWs would thus open a new observational frontier for testing new physics.

Traditional PTAs are sensitive to nHz GWs. To probe the pHz frequency regime, we use millisecond pulsars (MSPs) in a novel way: instead of analyzing the timing residuals, we search for drifts in MSP parameters induced by ultra-low frequency GWs. In this pHz frequency range, GWs have periods of hundreds to thousands of years. For observation times of only a few decades, we capture just a small segment of the gravitational waveform. These pHz GWs induce distortions in the spacetime around the pulsars that manifest as secular changes in the observed pulsar parameters, e.g. $\dot{P}_{b},~\ddot{P}$. $\dot{P}_{b}$ is the time derivative of the orbital period of a pulsar with a binary stellar companion, while $\ddot{P}$ is the second-order time derivative of the pulsar spin period.

The extensive data accumulated by PTA campaigns, combined with increasingly sophisticated pulsar noise models [Goncharov2024, Larsen2024, Agazie2023], provide an unprecedented opportunity to explore the pHz frequency regime. The concept that pHz GWs can induce measurable perturbations in pulsar timing parameters was first proposed by Bertotti1983, and later expanded upon by DeRocco2023, DeRoccoDror2024 in a modern theoretical framework. Here we present a new Bayesian framework that uses MSP timing parameters to search for CWs in the pHz frequency band. We develop an analytic model that describes the CW-induced drifts in the observed pulsar parameters $\dot{P}_{b}$ and $\ddot{P}$.

For our study, we extend an existing GWB model for SMBHBs [Casey-Clyde2022] to pHz frequencies. The population is matched in a self-consistent way to the evidence for a GWB in the NANOGrav 15-year data [nanograv]. From this extended model, we derive the characteristic strain of the background and generate synthetic realizations of pHz CW sources. Our analysis improves the upper limits on CW strain amplitudes by up to an order of magnitude relative to DeRoccoDror2024, and achieves a new leading constraint. Furthermore, through detailed simulations, we explore the CW scaling laws in the pHz regime and forecast strain amplitude sensitivities accessible by next-generation telescopes such as the Square Kilometre Array (SKA). We show that by timing new pulsars and improving timing precision with SKA, the technique presented here could potentially detect individual pHz sources, opening a new avenue to probe fundamental physics not accessible by any other means.

The paper is organized as follows. In Section II, we introduce the theoretical basis for the pHz GW detection method and the Bayesian framework. Section III describes the construction and extension of the SMBHB population model into the sub-nHz regime. Section IV presents the datasets and processing techniques employed in this analysis. In Section V, we summarize the results of our search. Section VI discusses the implications of our findings and prospects for future work.

II Methods

II.1 Analytical Model

In traditional PTA analyses, the lowest detectable GW frequency is determined by the inverse of the observation time, with $f_{\mathrm{GW}}^{\mathrm{min}}\geq T_{\mathrm{obs}}^{-1}$. Potential GW signals at frequencies below $f_{\mathrm{GW}}^{\mathrm{min}}$ are absorbed into the timing model noise during the fitting procedure. Nevertheless, these pHz GWs still induce systematic shifts in the measured pulsar parameters.

Specifically, at pHz frequencies, a passing GW induces an additional effective velocity between the pulsar and the Solar System Barycenter (SSB). The induced velocity is given by

[TABLE]

where the subscript $A=+,\times$ labels the GW polarizations (plus and cross, respectively), and $\mathbb{d}_{p}$ denotes the vector pointing from the SSB to pulsar $a$. The functions $h_{+,\times}(t,\mathbb{x})$ represent the GW strain at position $\mathbb{x}$ and time $t$. The antenna beam pattern, $F_{A}(\hat{\Omega})$, is defined as [Romano_2024]

[TABLE]

where $\hat{e}_{ij}^{A}$ is the polarization tensor of the GW, and $\hat{\Omega}$ is the unit vector of GW propagation. We work in the transverse traceless gauge and $i,j$ are the spatial indices.

Here we focus on the CW search using $\dot{P}_{b}$ and $\ddot{P}$. At sufficiently low GW frequencies, the additional velocity term induced by the GW also creates an apparent acceleration. This manifests as a constant offset in the pulsar timing parameters over the observation time. See Fig. 1 for a pictorial illustration.

To leading order, the observed derivative of the binary orbital period can be written as [Bell1996]

[TABLE]

where the terms on the right hand side correspond to, in order: $\dot{P}_{b,\text{int}}/{P_{b}}$, the intrinsic period derivative due to GW emission from the binary; $a_{\text{shk}}$, the kinematic Shklovskii effect due to the motion of the pulsar proportional to the relative motion of the pulsar perpendicular to the line of sight; $a_{\text{MW}}$, the acceleration due to the Milky Way gravitational potential; and $a_{\text{GW}}$, the contribution from pHz frequency GWs, which constitutes the signal we seek. In particular, $a_{\text{skh}}=\mu^{2}d_{p}$, where $\mu$ is the proper motion of the pulsar.

Similarly, the observed second derivative of the spin period is related to the jerk induced by GWs to leading order through

[TABLE]

In principle, there are kinetic and galactic contributions to $\ddot{P}_{\mathrm{obs}}$ as well, but the non-relativistic nature of the galaxy makes them subdominant to $j_{\mathrm{GW}}$ by several orders of magnitude [DeRoccoDror2024]. Intrinsic contributions due to magnetic dipole braking are of order $(\dot{P}/P)^{2}$, which typically evaluate to $\sim 10^{-35}s^{-2}$, much below the observed uncertainty of $\ddot{P}_{\mathrm{obs}}/P$ [LorimerKramer2012]. The higher order parameters such as $\dddot{P}$ can be derived in a similar fashion, and could be used for GW search in the future analysis. Explicit expressions for $a_{\mathrm{GW}}$ and $j_{\mathrm{GW}}$ induced by a CW are derived in Appendix A.

II.2 Bayesian Analysis

For this study, we consider GWs emitted by SMBHBs in a circular orbit. However, the method is easily generalizable to other CW sources. We have used 7 parameters characterizing the CW signal: $h_{0}$, the GW strain amplitude; $f_{\text{GW}}$, the GW frequency; and $\mathbb{\Theta}$, which consist of $\{\alpha,\delta,i,\phi,\Psi_{0}\}$. In order $\mathbb{\Theta}$ include the right ascension and declination $\alpha,\delta$, the orientation of the GW source $i,\phi$, and the orbital phase of the SMBH binary $\Psi_{0}$. Due to the relatively high dimensionality of the problem, we adopt a Bayesian approach.

We construct a likelihood making the standard assumption that pulsar uncertainties are uncorrelated [LorimerKramer2012]. The likelihood is given by

[TABLE]

where $\vec{y}_{\mathrm{GW}}$ the collection of accelerations or jerks $y_{\mathrm{GW},p}=\{a_{\mathrm{GW},p},j_{\mathrm{GW},p}\}$ [calculated using Eq. 3 or 4] with corresponding uncertainty $\sigma_{p}$ for each pulsar $p$. $N$ is the number of pulsars in the data.

Bayes’ theorem states that the posterior distribution of parameters $\theta$ given data $D$ is:

[TABLE]

where: $P(\theta|D)$ is the posterior distribution, $P(D|\theta)$ is the likelihood function, $P(\theta)$ is the prior distribution, and the denominator is the evidence.

In this study, we adopt uniform priors for all parameters, as listed in Table 1. However, we note that informative priors could apply if, for example, MMA observations are available. Particular care should be given to the strain amplitude $h_{0}$. By the nature of a GW search where the signal is weak, $h_{0}$ resides in a linear space spanning multiple orders of magnitude. When looking for a GW signal, we use the logarithmic uniform prior so that smaller values of $h_{0}$ are sampled sufficiently. In the case where the data are not sufficiently informative so that the likelihood function goes flat for small values of $h_{0}$, we use the linear-uniform prior, which corresponds to a prior in the log space where small values of $h_{0}$ are exponentially suppressed. It is worth noting that if the data are informative enough, the prior should not influence the posterior distribution after convergence. However, when the data are too noisy to support a detection, the linear flat prior gives a conservative estimate of the upper limit of $h_{0}$ in the sense that false positives tend to be avoided. The full prior is the product of the individual priors:

[TABLE]

We run Markov Chain Monte Carlo (MCMC) to determine the posterior distribution and interpret the posterior for either detection of a signal or setting an upper bound using existing data. Our method is naturally poised for using more informative priors based on multi-messenger astronomy observations, the details of which are beyond the scope of this paper. Such informative priors will help break the degeneracies and pick out the signal more accurately.

To characterize a detection, we employ a combined analysis based on the Bayes factor (BF) and the highest posterior density region. Specifically, the marginalized posterior distribution for the strain amplitude must be clearly distinct from the null hypothesis to declare a detection. We compute the BFs using the Savage-Dickey density ratio method. To set a detection threshold, we simulate the null distribution of the BF for the data under analysis.

If the marginalized posterior for the strain amplitude were to support a CW detection, we could infer the other relevant source parameters by analyzing the high density regions of the posterior distributions. Since the current data do not support a CW detection, we defer detailed parameter estimation to future analyses.

III Supermassive Black Hole Binary Population Model

In order to understand the detection prospects of our method, we extend the SMBHB population synthesis model from [Casey-Clyde2022] to the pHz regime. This model takes the GWB amplitude measurement as an input, and uses it to constrain SMBHB population characteristics, including the minimum black hole mass contributing to the background, the number density of SMBHBs per unit volume, and the volume that encloses the GWB. Here we use this population synthesis code to anchor the pHz SMBHB population to the nHz population, creating a self-consistent model of our GW sources.

In this model, we use dual active galactic nuclei (AGN) as tracers of SMBHBs: we assume a fixed ratio between dual AGN density $\Phi_{\mathrm{Q}}$ and SMBHB density $\Phi_{\mathrm{SMBHB}}$, depending on a mass ratio distribution $p(q)$:

[TABLE]

where $M_{\mathrm{BH,1}}$ is the mass of the primary, $M_{\mathrm{BH}}$ is the total mass of the binary, $z$ is the redshift, and $q$ is the mass ratio. The model parameters [Casey-Clyde2022] including the constant of proportionality $R$ are then extracted by fitting the integrated GWB characteristic strain to the evidence of GWB in the NANOGrav 15-year data.

Towards the pHz regime, several effects that shrink the SMBHB orbits are likely to compete with or even dominate over the GW emission of the binary system. The dual–AGN that contain the SMBHs which will eventually comprise the binary initially evolves under dynamical friction [BinneyTremaine2008]. As the SMBHs approach a separation of $\sim 100$ pc, the evolution of the dual–AGN becomes dominated by three-body interactions with individual stars, also known as “stellar hardening”. At separations of $\sim 1$ pc stellar hardening becomes inefficient as the SMBHs clears their vicinity of nearby stars. At this stage interactions between the SMBHB and nearby gas may become important for further shrinking the SMBH separation, though the details of this process remain an open question. Hard binaries of SMBHs start to form at this separation. For a SMBHB of mass $10^{8}\sim 10^{10}\mathrm{{M}}_{{\odot}}$ and separation $\sim 1~\text{pc}$, the emitted GW frequencies are $\sim 10-100~\text{pHz}$.

In our model, we include both dynamical friction and stellar hardening effects. We do not include the gas hardening effect and leave it for future investigation. Our method of pHz GW detection would be an important step to probe the various models of SMBHB evolution, manifested as a turnover that deviates from the power-law spectrum of GWB in the PTA frequency band. The GWB from our population model is orders of magnitude below the CW upper limits, so the noise due to GWB can be reasonably ignore in this work.

IV Dataset

For the $\dot{P}_{b}$ analysis, we use the updated data from [Moran:2023], which includes 29 pulsars. For our purposes, we independently estimate the Galactic accelerations of the pulsars using the MWPotential2014 model as implemented in the galpy package [Bovy2015].

The original data set from [Chakrabarti_2021] used in DeRoccoDror2024 contains 14 pulsars, 13 of which are also included in our updated sample. After cross-checking with other recent data releases, we adopt the more updated measurements from [Moran:2023] for all overlapping pulsars. In total, we obtain 30 pulsars with updated timing measurements, which is double the number used in the analysis of DeRoccoDror2024. The observed $\dot{P}_{b}$’s range in order of magnitude from $10^{-21}$ to $10^{-8}$, with the majority being around $10^{-14}\sim 10^{-12}$. We subtract non-GW contributions from the observed $\dot{P}_{b}$ and report the GW induced $\Delta\dot{P}_{b}/P_{b}$ data immediately used for analysis. The timing baselines of the pulsars range from 4 years to 22 years, with $\mu s$ level RMS timing residuals. For a full list of the pulsars, see Table 2.

Pulsar distance measurements are essential for reducing degeneracies in this analysis. If the resulting dataset is large enough, we require $f\Delta d_{p}\leq 1$; namely, we constrain the pulsar distance up to a precision less than a radiation wavelength, which helps mitigate the degeneracies associated with $2\pi$ periodicity of trig functions. For the pHz regime, this typically corresponds to a wavelength of $\sim 100$ pc. Half of the pulsars in the $\dot{P}_{b}$ data meet this criterion and are used in the analysis, as we use $f_{\mathrm{GW}}=100$ pHz to define the cutoff. When we perform our MCMC analysis, we use the pulsar distance measurements as a prior, and sample in the restricted error range in the analysis.

Additionally, to ensure that the Taylor expansions used in the analytic modeling remain valid, we cut off the analysis at $(4T_{\mathrm{max}})^{-1}$, consistent with [DeRoccoDror2024], where $T_{\text{max}}$ is the longest observation time in the data.

We account for the Shklovskii contribution $a_{\text{kin}}$ using the analytic formula in Section II. The uncertainty in $d$ dominates the uncertainty in $a_{\text{kin}}$, allowing us to treat $\mu$ as effectively fixed with a small uncertainty. When sampling $d$, we will propagate the proper motion uncertainty and add it in quadrature to the total uncertainty in $a_{\text{kin}}$. We verify that our cut on $\dot{P}_{b}$ data has resulted in a better upper limit by simulating null distributions, finding that including all pulsars shifts the null Bayes factor distribution towards larger values, which in turn biases the posterior toward higher strain amplitudes. The values of $\mu$ used in our analysis are taken from [Moran:2023].

For the $\ddot{P}$ analysis, we use the same data as in [DeRoccoDror2024], given by [Desvignes2016] and [Reardon2016]. A full list can be found in Table 3.The pulsars have timing baselines ranging from 7 to 18 years, with a typical RMS timing residual of a few $\mu s$. We use the parallax measurement of the pulsar distance whenever it is available. Otherwise, we estimate using dispersion measure (DM) and a model of Galactic free electron density model. A total of 13 pulsar distances in the data are estimated using DM and the Galactic free electron density model YMW 17 [yao2017new]. Since they lack a well-defined uncertainty in distance, we conservatively assume a $40\%$ relative uncertainty on these distances. Only one fourth of the pulsar distances in the $\ddot{P}$ data satisfy $f\Delta d_{p}\leq 1$, and dropping the majority of the data will make the set prone to outliers and lose information from a larger data set. Thus, we retain all pulsars in the $\ddot{P}$ data for analysis.

We also include an estimate of the red noise following the procedure described in [DeRoccoDror2024]. The red noise is treated as a stochastic background given by a power-law spectrum, and we estimate the ensemble variance of the two timing parameters of interest due to the red noise. For the current data, red noise contributions are negligible in the $\dot{P}_{b}$ measurements but relevant for $\ddot{P}$. We incorporate the red noise contribution into the total uncertainty budget for $\ddot{P}$. We emphasize that as timing precision improves, red noise may no longer remain subdominant to other sources of uncertainty in $\dot{P}_{b}$ and will need to be treated carefully in future analyses. The full data can be found in Appendix C.

V Results

We search for CW signals using $\dot{P}_{b}$ and $\ddot{P}$ independently. For the $\dot{P}_{b}$ data, a total of 15 pulsars satisfy our criteria in Sec. IV, which we use for GW searches and to place upper limit on the strain amplitude. To determine a detection threshold, we calculate the null distributions for both data sets through simulations with $1,000$ samples each. The null distributions exhibit small tails extending beyond $\text{BF}=1$, with $100\%$ of the samples falling below $\text{BF}=3$. We choose the conservative $\text{BF}=3$ as the detection threshold. (See Fig. 3)

We then search for CW signals in the real data, running up to $3,000,000$ samples in MCMC and discarding the first $20\%$ as burn-in. For the search we use log uniform prior on $h_{0}$. No CW signals are found in either data set. Fig. 4 and Fig. 5 show the posterior distributions of the relevant parameters for $\dot{P}_{b}$ search and $\ddot{P}$ search respectively. We recover the priors from the analysis except for a constrained strain amplitude. We report $BF=0.6913\pm 0.0058$ for the $\dot{P}_{b}$ search and $BF=0.4728\pm 0.0042$ for the $\ddot{P}$ search. Based on the null distribution and the posterior distributions, we do not detect a signal.

V.1 Upper limits

In the absence of detection, we set upper limits on the CW strain amplitude as a function of frequency by taking the $95\%$ cutoff of the marginalized posterior distribution of $h_{0}$, as shown in Fig. 2. For setting the upper limit we adopt a linear uniform prior on $h_{0}$. Specifically, for each frequency bin we apply a narrow Gaussian prior centered on that frequency with a standard deviation of 0.1 in the log space of frequency. The prior on $h_{0}$ is chosen to be uniform in linear space, which guarantees that the posterior is agnostic to the lower cut-off of the $h_{0}$ prior. If a log-uniform prior were used, a smaller lower cut-off in $h_{0}$ prior would always result in a smaller $h_{0}$ upper limit.

With our new analytic model, Bayesian framework, and the updated $\dot{P}_{b}$ data, we improve the $h_{0}$ upper limits by 9.53 times at a reference frequency of $f_{\text{GW}}=100~\text{pHz}$, an order of magnitude over previous results [DeRocco2023]. For the current data, the $\ddot{P}$ analysis achieves tighter strain limits than $\dot{P}_{b}$ for $f_{\text{GW}}\gtrsim 1.5~\text{pHz}$. Nevertheless, a combined analysis across multiple timing parameters remains valuable: when the sensitivity reaches levels comparable to expected CW signals, joint analysis of $\dot{P}_{b}$, $\ddot{P}$, and potentially other parameters could break degeneracies in source properties including CW frequency and sky locations, providing stronger constraints on strain amplitude.

On top of the fully marginalized upper limit, we also calculate the sky location dependent upper limit, presented in Fig. 7. This is done by fixing $\alpha,\,\delta$ in the Bayesian analysis for the corresponding sky location. In case of no detection, we expect the sensitivity distribution to correlate with the pulsar distribution by virtue of the antenna pattern function defined in Eq. 2. Our result conforms to this expectation. In the calculation presented here, we use a reference GW frequency of $100~\text{pHz}$. We note that the detection sensitivity can be different by up to an order of magnitude in different regions of sky. The upper limit in strain is as low as $h_{0}\simeq 6.7\times 10^{-13}$ for particular sky locations, one seventh of the marginalized result over all sky at the same frequency.

V.2 Scaling laws

The scaling of the upper limit $h_{0}$ with GW frequency $f_{\text{GW}}$ follows $h_{0}\propto f_{\text{GW}}^{-1}$ for $\dot{P}_{b}$, and $h_{0}\propto f_{\text{GW}}^{-2}$ for $\ddot{P}$ (see Fig. 2). These scaling relations are consistent with our analytic model, where $a_{\text{\text{GW}}}\propto h_{0}f_{\text{GW}}$ and $j_{\text{GW}}\propto h_{0}f_{\text{GW}}^{2}$.

At very low frequencies, the scaling of sensitivity changes. For $f_{\text{GW}}\lesssim 10~\text{pHz}$, the GW wavelength exceeds $d_{p}^{-1}$ for the most sensitive pulsars. In this regime, the earth term and pulsar term in Eq. 1 become “in phase,” which results in a weakened $h_{0}$ upper limit for extremely low frequency.

V.3 SMBHB synthesis model

The integrated GWB spectrum predicted from our SMBHB population model deviates significantly from a simple power-law extrapolation in the sub-nHz regime. In particular, a turnover in the spectrum appears around $60~\text{pHz}$, corresponding to the strongest GWB signal across the entire frequency range. A simple power-law assumes that GW emission dominates SMBHB orbital evolution, while our model incorporates environmental effects such as stellar hardening and dynamical friction. At pHz frequencies, these environmental processes dominate, shortening the binary residence timescale ($f_{\mathrm{orb}}/\dot{f}_{\mathrm{orb}}$) compared to purely GW-driven evolution (Fig. 6), where $f_{\mathrm{orb}}$ is the binary orbital frequency. The shorter residence time leads to a suppression of GWB power at the lowest frequencies, manifested as a reduced characteristic strain in Fig. 2.

We generate $1,000$ realizations of the SMBHB population to estimate the statistical uncertainty, represented by the quantile shading in Fig. 2. The distribution of SMBHBs across frequency bins for different masses follows a similar pattern (Fig. 8), with a turnover occurring near $50~\text{pHz}$, matching the feature in the GWB spectrum. The overall population decreases with increasing SMBHB mass, consistent with the expected mass function.

We extract the loudest individual CW sources in each frequency bin. We find no strong scaling between the strain amplitude of the loudest sources and frequency. To extract these sources, we find the largest strain amplitude from all the bins that contains a SMBHB, for each frequency. The loudest sources typically exhibit strain amplitudes around $1\sim 4\times 10^{-16}$. The majority of individual sources lie several orders of magnitude below this level.

V.4 Future Detection Prospects

To forecast detection prospects, we conduct simulations to determine the scaling relations between detection sensitivity and (1) data size and (2) data uncertainty levels. To generate the simulated pulsar data, we draw from the extrapolated sky location distribution of known pulsars inspired by method described in [Xin2021]. Both conservative and more optimistic estimates for SKA1 and 2 are included here. The SKA LOW survey will concentrate on regions within $\pm 5^{\circ}$ from the galactic plane, while the SKA MID survey will concentrate on regions within $\pm 10^{\circ}$ of the galactic plane [Braun:2015zta]. Upcoming experiments will dramatically expand the PTA by increasing the number of precisely timed MSPs to the order of several hundred, representing more than an order-of-magnitude increase in size relative to current data [Xin2021]. Distance measurements will also improve significantly with SKA astrometry [Smits:2011], reducing uncertainties to levels that satisfy our distance-precision criterion for essentially all pulsars. According to [Xin2021], SKA1 is expected to produce precise measurements for 675 pulsars, with an additional 120 achieved by SKA2. Because SKA1 consists of southern hemisphere telescopes, we apply a cut-off at $30^{\circ}$ in the declination angle of the pulsars. We then draw pulsar locations from our empirical sky location distribution extracted from all the pulsars in IPTA data [Antoniadis2022].

We verify that, when probing strain amplitudes above potential CW sources, uncertainties in pulsar distances alone do not significantly impact the upper limits on $h_{0}$. Each simulation assumes a null CW signal. We find the detection sensitivity scales as $\propto\sigma_{p}/\sqrt{N}$ when $N$ is asymptotically large (above $\sim 100$), where $\sigma_{p}$ denotes the uncertainty in data used for inference and $N$ is the number of pulsars (Fig. 9). When $N$ is below 50, boundary effects kick in and the scaling deviates from $\propto 1/\sqrt{N}$—the exponent on $N$ in the power law scaling is slightly smaller than $-1/2$.

For a conservative forecast, we assume each pulsar retains the typical uncertainty level seen in the current data. We further assume that $80\%$ of the $\sim 800$ precisely timed pulsars are in binary systems [Lorimer2008], each with a measurable $\dot{P}_{b}$, and that all have measurable $\ddot{P}$. In addition, we take the observation time for all the new pulsars to be 10 years, which is below the longest observation time in the current data. Under these assumptions, we compare the forecast sensitivities derived from $\dot{P}_{b}$ and $\ddot{P}$ analyses with the strain amplitudes of the loudest individual sources from our SMBHB population model. As shown in Fig. 2, SKA-based data will enable an additional order-of-magnitude improvement over current sensitivities, solely due to the increase in data size. Our simulations also show that with more pulsars, the detection sensitivity will become more isotropic (Fig. 10). SKA data is expected to reduce the variance between the most and least sensitive sky locations by $\sim 40\%$ for $\ddot{P}$ analysis. However, the strain amplitudes of individual SMBHB CW sources drawn from our GWB synthesis model are still $\sim 3$ orders of magnitude below this forecast sensitivity.

V.5 Injection of a CW

To illustrate that our method is not only able to set upper limits, but capable of detecting an individual CW source, we draw the loudest SMBHB from our population and simulate its GW imprints on the timing parameters. We then retrieve the CW signal from the simulated data using our method. We focus on $\ddot{P}$ in this section and predict the required precision of $\ddot{P}$ measurement for a realistic detection.

In our simulations, we choose NGC 6240 as the host galaxy of the SMBHB. NGC 6240 is a nearby ultraluminous infrared galaxy and prominent candidate for dual AGN, which in turn strongly indicates the presence of a SMBHB [sobolenko2022ngc6240]. We assume the CW from the SMBHB has the strongest strain amplitude predicted by our population model and a redshifted frequency of $300~\text{pHz}$.

We use BF as an indicator of detection. In Fig. 11, we show how the BF scales with the number of pulsars. The blue line assumes the current uncertainty level in the data, the black dashed line assumes an universal improvement in the uncertainty of measured $\ddot{P}$ by 2 orders of magnitude in all pulsars, and the red dotted line is the detection threshold set by the null distribution. It is evident that the increasing number of pulsars alone by SKA will not facilitate a detection—more precise measurements of pulsar timing parameters are needed.

To forecast detection prospects with improving data uncertainty, we carry out simulations under two scenarios, as shown in Fig. 12. The first is when the parameter measurements for all the pulsars improve at the same rate, plotted as a black solid line, while the second is considering a scenario where 10 pulsars have improving data uncertainty and the remaining pulsars have their precision fixed at the typical uncertainty in current data. Our results suggest a uniform improvement of two orders of magnitude on the existing sensitivity ($\delta(\ddot{P}/P)\sim\times 10^{-30}s^{-2}$) is necessary, or, alternatively a three order of magnitude improvement for 10 out of the 795 pulsars. So there is a clear trade-off between precision required and number of pulsars that can achieve this precision—the larger the number of pulsars achieving the precision, the lower the precision requirement. Though it may seem as if a uniform improvement of two orders of magnitude would be easier to achieve, it may actually prove the case that an improvement on only 10 pulsars is more realistic. Measuring $\ddot{P}$ requires extremely stable pulsars and long timing baselines, hence it is likely that even with SKA data, only a select few pulsars will be able to provide precise measurements of their $\ddot{P}$.

In Appendix B, we show the detailed search results of the two cases discussed above, along with the results of a search using 795 pulsars with typical uncertainty of current $\ddot{P}$ data, in which we detect no signal (Fig. 13). Fig. 14 assumes that 10 of these pulsars achieve 3 orders of magnitude improvement in $\ddot{P}$ uncertainty, and is able to recover $h_{0}$ fairly well. This search also produces posteriors with peaks near the true values for other parameters, such as $f_{\mathrm{GW}}$ and SMBHB sky location. The high degrees of degeneracies between the various parameters lead to multiple peaks in the posteriors. We verify that even smaller uncertainties in the data eventually make the posteriors converge to the injected parameter values. Fig. 15 assumes that all of 795 pulsars achieve 3 orders of magnitude improvement in $\ddot{P}$ uncertainty, and reaches better convergence to the true values of the injected CW parameters.

In summary, we find that to detect a realistic CW source, a $\sim 2$ order of magnitude improvement in the uncertainty of $\ddot{P}$ is required for all pulsars that SKA will observe. Interestingly, however, a trade-off exists between the number of well-measured $\ddot{P}$s and the required precision; in particular, we show only 10 out of 795 pulsars in SKA would need to improve by $\sim 3$ orders of magnitude in order to achieve a detection. In which direction this trade-off will play out when SKA begins taking data remains an open question.

VI Discussion and Conclusions

As difficult as the required improvement in timing parameter measurement for a future detection may seem, multiple factors work in our favor. Current studies indicate the timing precision of SKA could be up to 2 orders of magnitude better than today [liu2014pulsar]. On top of that, it has been empirically shown that a longer observation time reduces the standard deviations in timing parameters [andrews2020nanograv], and simulations suggest a scaling of $\sim T_{\mathrm{obs}}^{-1/2}$ [liu2019correlated]. The typical timing precision of current operational PTAs is $\sim 1~\mu s$, and the current uncertainty of the period normalized $\dot{P}_{b}/P_{b}$, $\ddot{P}/P$ sits at $10^{-20}s^{-1}$, $10^{-30}s^{-2}$. Optimistically the uncertainties in the timing parameters scale linearly with the timing precision. With 2 orders or magnitude improvement in the timing precision, we would achieve $\sigma(\dot{P}_{b}/P_{b})\sim 10^{-22}s^{-1}$  /  $\sigma(\ddot{P}/P)\sim 10^{-32}s^{-2}$. The observation times of current PTAs range from a few years to two decades. Longer timespan will only improve the uncertainty in timing parameters, though the improvement due to longer timespan will be incremental compared to the improvement due to better timing precision by the scaling relations. Referring to Fig 12, it is likely that some potential CW sources may fall within our sensitivity range with SKA data.

However, there are a few caveats. Currently, the dominant sources of uncertainty in both $\dot{P}_{b}$ and $\ddot{P}$ are due to the measured PTA timing parameters, which are expected to improve with longer time baselines and next-generation telescopes such as SKA. As timing precision improves, however, red noise and uncertainties in galactic acceleration will not be negligible. These effects present their own challenges and motivate the exploration of higher-order timing parameters with better-controlled uncertainties. More details of recent developments can be found for red noise mitigation in [Goncharov_2020] and for Galactic acceleration modeling in [ibata2023chartinggalacticaccelerationfield]. In particular, since red noise is covariant with a GWB signal, it will be interesting to include a GWB when extracting timing parameters, which has not been done. We also acknowledge that only a very small fraction of pulsars will have precise measurements of higher-order timing parameters, due to the stringent requirement of pulsar stability and long timing baseline. However, we have shown in Sec. V.5 that a small number of precise measurements could potentially facilitate a successful detection.

Targeted searches may also help improve detection prospects [Agarwal:2025cag]. In our framework, targeted searches are achieved by adopting informative priors on CW parameters. These priors are often informed by multi-messenger astrophysical observations. For example, recent work has shown that dual AGNs can serve as tracers of SMBHBs, where the sky location of candidate dual AGNs can be determined through optical, radio, and X-ray observations, while redshifts are measured spectroscopically [Casey-Clyde2022]. Moreover, the orientation of SMBHB candidates can be constrained through observations of jet alignment [Dotti2010], broad-line region kinematics [Afanasiev2019], variability and polarization of quasar light [Hutsemekers2005], and host galaxy morphology [Zakamska2006].

The framework presented here can also be used to search for alternative deterministic GW signals by adapting the CW model parameters. This opens opportunities to search for a variety of exotic physics, including ultralight dark matter and dark matter substructure [Brito2015, NANOGrav2023NewPhysics].

The exact effect of gas hardening on binary formation and evolution is still an open research question. For our GWB synthesis model, the energy loss of the binary by gas hardening is likely to be comparable with GW emission and stellar hardening at a transitional binary separation of $\sim 1$ pc [Koss2023UGC4211, Izquierdo‐Villalba2021]. Below this separation, the energy loss is dominated by GW emission; above this separation, energy loss is dominated by stellar hardening. Inclusion of gas hardening in the model may increase the turnover frequency in the GWB spectrum, since it will further reduce the residence time of a binary and thus GWs emitted at a certain frequency. Towards the very beginning of binary formation, gas hardening is expected to be subdominant to stellar hardening and dynamical friction, so the low-frequency asymptotic behavior of our GWB spectrum should remain unaltered.

Looking forward, there are several promising extensions to pursue. First, a dedicated analysis targeting the pHz GWB could provide independent constraints on the properties of sources such as SMBHBs, topological defects, or early-Universe phase transitions. Second, a key technical challenge in CW searches is the presence of a background signal, which can act as a source of contamination. Future work should incorporate the treatment of an unresolved GWB as a correlated noise component when searching for individual CW sources, a step necessary for maintaining sensitivity in a background-dominated regime. Third, the joint analysis of multiple timing parameters, such as $\dot{P}_{b}$, $\ddot{P}$, and higher-order derivatives, offers a pathway to break degeneracies between CW source properties and enhance detection significance. Finally, extending the framework to accommodate alternative GW morphologies, including non-monochromatic signals and bursts with memory, would substantially broaden the range of astrophysical and cosmological phenomena accessible to pHz GW searches. Each of these future developments is essential for maximizing the scientific return of next-generation PTA data and for exploiting the full potential of this method.

In summary, we have carried out the most comprehensive search to date for pHz gravitational waves using secular pulsar parameters, achieving an order-of-magnitude improvement in strain limits over earlier studies. By introducing a Bayesian detection framework, doubling the pulsar sample, and connecting the analysis directly to population synthesis models, we establish this approach as a robust tool for probing the ultra-low frequency regime. Looking ahead, the dramatic increase in pulsar discoveries and timing precision with SKA will transform these methods from upper-limit exercises into genuine discovery opportunities. This opens a realistic path to testing models of supermassive black hole binary evolution, while simultaneously constraining exotic physics such as topological defects and early-Universe phase transitions.

Acknowledgements.

We thank Jeffrey S. Hazboun for the discussion of SKA pulsar population forecasting. The research of J.D. is supported in part by the U.S. Department of Energy grant number DE-SC0025569. W.D. was supported by NSF grant PHY-2210361 and the Maryland Center for Fundamental Physics.

Appendix A Derivation of analytic model

GWs are ripples of spacetime, and thus will induce an apparent velocity between the pulsars and the solar system barycenter. This, and its derivatives, is the main quantity used to detect a signal of GW in this method discussed. Here we derive the observable shifts in timing parameters due to the CWs emitted by a SMBHB, which is assumed to be the main source of CW in this paper. The general idea can be adapted to model other sources.

We start from the well known expression

[TABLE]

where $\hat{\Omega}$ is the unit vector of GW propagation direction, while $\hat{d_{p}}$ is the vector pointing to the pulsar under consideration. We assume the CWs are emitted by the fiducial circular orbit SMBHBs

[TABLE]

where $\mathscr{M}$ is the chirp mass, $d_{L}$ is the luminosity distance.

[TABLE]

where $\phi$, $i$ defines the orientation of the GW source, and $\Psi=\Psi_{0}+\omega_{0}t$, with $\omega_{0}$ the orbital frequency of the binary gravitational wave source.

Expand by virtue of equation 2, we have

[TABLE]

The direction of the source can be parametrized by the right ascension $\alpha$ and declination $\delta$. In the Cartesian coordinates,

[TABLE]

and

[TABLE]

In the analysis, we need to use the acceleration and jerk $a_{GW}$, $j_{GW}$, which we derive here. Since we are working in the pHz frequency range with typical periods much longer than the observation time, we can safely assume that $\omega_{0}t<<1$ in our derivations, and thus expand $\cos(2\omega_{0}t)$ as 1, and $\sin(2\omega_{0}t)$ as $2\omega_{0}t$. Now we have

[TABLE]

[TABLE]

The unit vectors $\hat{d_{p}}$ and $\hat{\Omega}$ can be parametrized by the equatorial ascension($\alpha_{p},\alpha$) and declination ($\delta_{p},\delta$) of the pulsars and source in the Cartesian coordinates, which read

[TABLE]

One can then write down explicit expressions for $a_{\text{GW}},j_{\text{GW}}$ in terms of the parameters above, for which the numerical values are given in the data.

Appendix B CW detection

Detailed search results using simulated data with CW signal injected, as described in Sec. V.5.

Appendix C Pulsar data

Table 3 includes the data used for $\ddot{P}$ search, and Table 2 includes the data used for $\dot{P}_{b}$ search. By the criterion described in Sec. IV, the 15 pulsars in bold font in Table 2 are used in setting the upper limits.