A Joint Optimal Search for Gravitational Waves from Resolved and Unresolved Supermassive Binary Black Holes with Pulsar Timing Arrays

Abstract: We introduce, from first principles, a joint model of the gravitational wave background (GWB) and brightest supermassive black hole binary (SMBHB) sources that may be individually resolvable in Pulsar Timing Array (PTA) searches for gravitational waves. We propose the characteristic number of SMBHB sources, $N_{\rm c}$, as a detection statistic for the astrophysical origin of the GWB. We then demonstrate how the brightest SMBHBs assist in resolving $N_{\rm c}$. Applying our method to the simulated NANOGrav 15-year data, which replicates all aspects of real data's known noise, observations, and the inferred GWB power spectrum, we demonstrate direct astrophysical limits on the strain amplitude of individually resolvable SMBHBs. We find that 21 of 114 SMBHB candidates from active galactic nuclei observations are in tension with the NANOGrav's observations. In contrast, only one candidate is in tension with the NANOGrav data based on the upper limits reported in the original analysis. Constraining the Poisson-specific characteristic number of SMBHBs, $N_{\rm c}$, at ${\rm yr}<sup>{-1}$, we outline implications for the population properties of SMBHBs. Based on our new model applied to the simulated NANOGrav data, we calculate the probability of detecting GWs from isolated SMBHB in the 15-year data to be 2\% at the ${\rm SNR}=5$ level. Our projection towards the expected NANOGrav 20-year data suggests an increase to 5\%. With this, we estimate the probability of finding an outlier with an SNR of 2 in the NANOGrav 20-year data to be $40\%$.

• The paper presents a unified Bayesian framework that jointly analyzes stochastic and resolvable gravitational wave signals in pulsar timing data.
• It employs astrophysical priors, including characteristic source numbers and strain amplitudes, to enhance parameter estimation and challenge Gaussian assumptions.
• The results yield stringent upper limits on continuous wave strains and improve detection prospects for individual SMBHB signals in future PTA observations.

A Joint Hierarchical Model for SMBHB Gravitational Wave Searches with Pulsar Timing Arrays

Introduction and Motivation

This work presents a comprehensive hierarchical statistical framework to search for nanohertz gravitational waves (GWs) generated by supermassive black hole binaries (SMBHBs) using pulsar timing arrays (PTAs). Rather than segmenting PTA data analysis into purely stochastic GW background (GWB) searches and independent resolved continuous wave (CW) searches, this model jointly describes both the unresolved stochastic background and potentially individually resolvable sources within a unified Bayesian framework. Such an approach leverages astrophysical priors directly informed by population synthesis, population hyperparameters, and the physics of SMBHB inspirals, addressing critical limitations in current PTA analyses which separate the inference of stochastic backgrounds and bright resolvable binaries.

Astrophysical Model for the Complete GW Signal

The formalism starts with detailed characterization of the contributions to the PTA GW signal from both unresolved sources forming the GWB and the brightest, possibly resolvable, SMBHBs. The stochastic GWB is understood as a sum over a cosmological population of SMBHBs, whose individual strains add incoherently, with the total strain described as a stochastic process. Deviations from Gaussianity emerge when the number of contributing binaries per frequency bin is low—an inevitable consequence at the highest frequencies accessible to PTAs due to the decreasing temporal window per bin and the decreasing number density of sources.

The model is effectively parameterized by two key astrophysical hyperparameters:

• $N_{\rm c}$, the characteristic number of sources contributing to the background in a reference frequency bin,
• $h_{\rm c}$, the characteristic strain amplitude at the reference frequency.

The probability density function (PDF) of the net characteristic strain, and the PDF for the brightest source’s strain, are constructed analytically using the theory of Poisson processes and source intensity functions.

Figure 1: An illustration of the PDF of the total characteristic strain squared for SMBHBs, showing the distribution for the total, the brightest, and all-but-the-brightest sources as a function of frequency.

A pivotal aspect is that, at high $N_{\rm c}$, the total strain distribution converges to Gaussian, whereas for small $N_{\rm c}$, finite population effects induce pronounced spectral non-Gaussianity and strong outliers—the onset of which signals the potential resolvability of individual binaries. The model utilizes analytic and numerical results for the SMBHB luminosity function, connecting physical parameters (e.g., mass, redshift, chirp mass distributions) to the observed strain PDFs.

Hierarchical Data Analysis Approach

The PTA likelihood is constructed as usual for ensemble pulsar timing residuals, but—in contrast to standard approaches—the prior on the set of frequency-binned GWB and CW amplitudes is replaced with the joint astrophysical prior conditioned on $(N_{\rm c}, h_{\rm c})$. This enables consistent propagation of population-level information into the inference regarding both the background and resolvable sources, and vice versa. The framework offers direct calculation of evidence for or against the discrete Poissonian nature of the GWB as opposed to an idealized Gaussian model.

Figure 2: Characteristic number of SMBHB sources, $N_{\rm c}$, as a detection statistic for the astrophysical origin of the GWB. Posterior distributions illustrate transitions between non-Gaussian (finite-source) and Gaussian (large-source) regimes.

In synthetic studies, the methodology efficiently discriminates between cases dominated by a stochastic background (large $N_{\rm c}$, Gaussian regime) and those with apparent non-Gaussian power spectral outliers due to finite-source effects (small $N_{\rm c}$). The analysis thus operationalizes the detection of the breakdown of the Gaussian-background hypothesis, providing a robust test for the astrophysical origin of the GWB.

Parameter Estimation and Practical Sensitivity

Applied to datasets emulating the NANOGrav 15-year data, the joint model enables joint posterior estimation on $N_{\rm c}$, $h_{\rm c}$, and, if present, the amplitude and frequency of the brightest CW component. The analysis finds that, given current sensitivity, the NANOGrav data remain consistent with a purely stochastic background, with no compelling evidence for a resolved individual SMBHB. However, the framework sharply constrains the allowed strain amplitudes for potential continuous-wave sources, as well as the population-level parameters themselves.

Figure 3: The estimation of parameters of the joint model for GWB and CW sources as a function of the simulated CW brightness $h_{\rm c}$0.

Figure 4: Marginalized posterior for three parameters of the joint hierarchical model: characteristic number of SMBHBs $h_{\rm c}$1, characteristic strain $h_{\rm c}$2, and frequency of the fitted brightest source $h_{\rm c}$3.

A major outcome is the computation of astrophysically motivated upper limits on the strain $h_{\rm c}$4 of the brightest expected CW, using hierarchical priors. These limits are systematically stronger than generic, non-informative analyses, substantially increasing the statistical tension for a subset of AGN-selected SMBHB candidates purported in EM surveys.

Figure 5: Marginalized posteriors and posterior-predictive samples for the characteristic strain of unresolved GWB sources (left) and resolved sources (right) per GW frequency bin. At high frequency, the signal is dominated by the brightest sources, but sensitivity degrades due to rising noise.

Figure 6: Constraints on the CW strain amplitude $h_{\rm c}$5 for the brightest sources across the sky, showing improved (lowered) upper limits with the joint model relative to previous NANOGrav analyses. More AGN-selected SMBHB candidates are disfavored by the new limits.

Quantitatively, the computed probability of a detectable ($h_{\rm c}$6) individually resolved CW with current data is calculated to be 2%, increasing to 5% with an anticipated 20-year dataset, reflecting the underlying SMBHB population statistics and PTA sensitivity evolution.

Figure 7: Posterior-predictive signal-to-noise ratio $h_{\rm c}$7 and detection probabilities for CWs in the NANOGrav data. Maximum $h_{\rm c}$8 and RMS statistics across frequency bins are also shown.

Comprehensive tabulation of probabilities for CW signals exceeding SNR thresholds—taking posterior uncertainties in both astrophysics and instrument characteristics into account—enables direct computation of false-alarm and true-alarm rates for future PTA campaigns.

Astrophysical Interpretation and Comparison

Posterior inference on $h_{\rm c}$9 is recast in terms of physical properties of the SMBHB population: the peak chirp mass $N_{\rm c}$0 and the mass density of SMBHs $N_{\rm c}$1. This allows direct confrontation with population synthesis predictions, previously derived GWB amplitude priors, and empirical galactic scaling relations. The results of the hierarchical analysis are in good agreement with population-level constraints derived from EPTA data and with local black hole demographics, but with the prospect of breaking degeneracies between abundance and mass scale through future improved parameter estimation.

Figure 8: Posterior on the SMBHB mass density $N_{\rm c}$2 and the mass scale $N_{\rm c}$3, compared to astrophysical priors and population synthesis models. The thick brown line shows EPTA-based astrophysical priors.

The framework is readily extensible to incorporate effects such as non-circular orbits, environmental coupling, and time-varying population statistics—a direction of substantial importance for ongoing and next-generation PTA datasets.

Conclusion

This work establishes a rigorous, physically motivated joint hierarchical approach for PTA analysis of nanohertz GWs from SMBHBs, integrating both the stochastic GWB and resolved CW searches within a consistent Bayesian model. The method leverages detailed population-level priors, robustly tests the Gaussian-background assumption, and produces stringent, astrophysically interpretable upper limits on the strains of individual sources as well as background properties. The approach is both more powerful and more precise than segmented analyses, and offers a pathway for more definitive identification of the astrophysical GW background and its constituents with forthcoming PTA datasets.

(2606.18241)