Expectations for the first supermassive black-hole binary resolved by PTAs I: Model efficacy

TL;DR

This study evaluates how effectively three PTA SMBHB detection strategies can identify and characterize an individual binary amid a Gaussian GW background using IPTA-like simulations. It compares a full deterministic continuous-wave model that includes Earth- and pulsar-term effects (CW_FS), an Earth-term-only template (CW_ET), and a frequency-resolved Spike Pixel anisotropy approach (SP). Across 40 simulated datasets, CW_FS consistently yields higher Bayes factors and more precise sky localization and frequency estimates, while SP can offer early localization advantages in some low-frequency scenarios; CW_ET struggles at low frequencies due to pulsar-term interference. The results underscore the importance of incorporating the full signal structure for early and accurate detection, while anisotropy-based methods provide useful cross-checks and rapid localization in certain regimes. The work also discusses limitations of current Bayesian tools (e.g., Savage-Dickey Bayes factors) at higher S/N and outlines future directions, including transdimensional anisotropy models and incorporating HD correlations in full signal searches.

Abstract

One of the most promising targets for Pulsar Timing Arrays (PTAs) is identifying an individual supermassive black hole binary (SMBHB) out of the population of binaries theorized to produce a gravitational wave background (GWB). In this work, we emulate an evolving PTA dataset, complete with an increasing number of pulsars and timing baseline, into which we inject a single binary on top of a Gaussian GWB signal. We vary the binary's source parameters, including sky position and frequency, and create an ensemble of simulated datasets with which we assess current Bayesian binary search techniques. We apply two waveform-based template models and a frequency-resolved anisotropy search to these simulations to understand how they compare in their detection and characterization abilities. We find that a template-based search including the full gravitational-wave signal structure (i.e., both Earth and pulsar effects of an incident GW) returns the highest Bayes Factors (BF), exceeding our estimator's capabilities by (S/N)~9-19, and has the most robust parameter estimation. Our anisotropy model attains a realization-median BF>10 at 7<(S/N)<15. Interestingly, despite being a deterministic model, the Earth-term template struggles to detect and characterize low-frequency binaries (5 nHz). These binaries require higher (S/N)~16-19 to reach the same BF threshold. This is likely due to neglected confusion effects between the pulsar and Earth terms. By contrast, the frequency-resolved anisotropy model shows promise for parameter estimation despite treating a binary's GW signal as excess directional GW power without phase modeling. Sky location and frequency parameter constraints returned by the anisotropy model are only surpassed by the Earth term template model at (S/N)~12-13. Milestones for a first detection using the full-signal GW model are included in a companion paper (Petrov et al. 2025).

Figures (6)

• Figure 1: Sky map of all pulsar positions in our simulated array, along with the two sky locations into which we inject a CW signal. The pulsars are colored according to the time slice in which they are added to the dataset, with the earliest-added (and consequently longest-timed) pulsars represented by darker orange stars, while pulsars added to the array more recently are represented by progressively lighter orange stars. The injected sky locations for Simulations A and B are shown with a yellow triangle and square, respectively. By the 20-year dataset, the array includes all 116 pulsars.
• Figure 2: We show the median and full range of the Savage-Dickey Bayes factor for each binary search technique with respect to the dataset timespan. The corresponding $\mathrm{(S/N)_{\Lambda}}$ is plotted on the top x-axis. The full signal model $\mathrm{CW_{FS}}$ (orange circles) recovers higher BFs in comparison to the coarse SP anisotropy model (blue triangles) and the $\mathrm{CW_{ET}}$ (green $\times$s). When there is insufficient MCMC sampling coverage of the CW amplitude's posterior tail, an upper limit is reported for that realization. We do not report BFs if half the realizations require such upper limits.
• Figure 3: Since the Savage-Dickey Bayes factor becomes unreliable as the binary signal becomes stronger in the data, we turn to the ratio of log-likelihoods for each model to a HD GWB. Each log-likelihood was calculated using the injected values for each model respectively. Plotted are the medians over 100 noise realizations of the datasets, which include the ten used in the Bayesian analyses. Here we can observe the marked advantage that the full signal template model (orange dots) has over other models. We also show the true model (pink +)---which we did not apply to the dataset---to demonstrate how even the $\mathrm{CW_{FS}}$ model is impaired by misspecification by relying on a CURN model over HD for the GWB. The $\mathrm{CW_{ET}}$ model displays a higher likelihood ratio than SP, contrary to the measured BFs in \ref{['fig:sdbf_acrossmodels']}. In the background, a histogram displays the number of pulsars at each time slice of the dataset.
• Figure 4: We average over the four binary injections and all respective noise realizations to show the overall behavior of each model in constraining the GW source's sky location. The grey solid, dashed, and dotted horizontal lines denote 100%, 50%, and 10% of the sky, respectively. We also show a dashed black curve representing a $\mathrm{(S/N)}^{-2}$ relationship expected for sky localization SV2010. The colorful curves and markers represent the median 90% credible region across injections and realizations, while the shaded regions enclose the median range of values. The $\mathrm{CW_{FS}}$ model returns the strongest constraints on the sky location and does so earlier than the other models. SP is the second to do so, but is surpassed by the $\mathrm{CW_{ET}}$ model at the 15-year time slice. All models narrow the source's location to $\leq$ 10% of the sky by the final time slice, when $\mathrm{(S/N)_{\Lambda}}\sim22$.
• Figure 5: Realization-median 68% credible intervals as a function of time slice and $\mathrm{(S/N)_\Lambda}$ for all three models. In each panel, the true injected value is marked by a horizontal dashed black line. Of note are the biases incurred by the models in GWB parameter estimation in the 20 nHz injections and the lag in constraints returned by the $\mathrm{CW_{ET}}$ model for the 5 nHz injections. In the 20 nHz injections, SP retains a slight edge over $\mathrm{CW_{ET}}$, but both models are largely comparable. The SP model's frequency resolution is limited by the discrete binning, which shifts slightly with increasing array timespans. This creates the sawtooth pattern seen in the SP frequency constraints, wherein the CW is identified, but in each subsequent dataset the frequency bin location shifts slightly. Frequency bins are marked by the grey dotted lines.