Finding Supermassive Black Hole Binary Mergers in Pulsar Timing Array Data

TL;DR

This paper addresses detecting gravitational-wave memory from supermassive black hole binary mergers in pulsar timing array data by introducing a full inspiral–merger–ringdown (IMR) waveform with null memory generated from a numerical-relativity surrogate (NRHybSur3dq8_CCE). The authors develop a Bayesian analysis pipeline that projects the signal into PTA observables via the Earth term, models realistic pulsar noise and a stochastic gravitational-wave background, and uses parallel-tempered MCMC to sample multimodal posteriors. They demonstrate precise recovery of chirp mass and luminosity distance (to within ~0.1%) across simulated mergers and show strong memory-detection evidence (ln B > 10) in favorable scenarios, while also quantifying limitations of the traditional memory-burst model. The work highlights the importance of physically consistent merger waveforms for parameter estimation, notes potential extensions to spinning and eccentric SMBHBs, and discusses prospects for electromagnetic counterparts and the broader scientific impact of opening a new window on SMBHB populations with PTAs.

Abstract

Galaxy observations suggest there is about one merger of supermassive black hole binaries (SMBHB) throughout the observable universe in a year. Here, we introduce the methodology to search for gravitational waves from these events with Pulsar Timing Arrays (PTAs). Modelling the inspiral, the merger, the ringdown, and the gravitational wave memory components of the signal in simulated data, we demonstrate a proof of principle for detection and parameter estimation. We study a few representative SMBHB mergers with chirp masses spanning $10^{8} - 10^{10}~M_\odot$ at distances from a few Mpc to 100~Mpc to asses their detectability in PTA observations. Assuming the fixed binary inclination angle of $90^{\circ}$ corresponding to the maximum displacement memory signal, these signals appear distinct for a PTA with 25 pulsars timed for 13 years with 100 ns precision. We demonstrate the capabilities of PTAs to constrain chirp masses and distances of detected merging binaries, as well as to place limits. The sky position uncertainties of the order of $1^{\circ}$, which we find in this optimistic example, could potentially enable electromagnetic follow-up and multi-messenger observations of SMBHB mergers. Finally, we show that the measurement uncertainties on the parameters of simulated merging binaries depend weakly on the presence of the gravitational wave background with Hellings-Downs correlations in our simulated data.

Figures (12)

• Figure 1: Gravitational wave strain and timing residuals from a merger of a non-spinning supermassive black hole binary (SMBHB) with parameters $\mathcal{M}_c = 10^{10} M_{\odot}$, $q=1$, and $D_\text{L} = 1000$ Mpc. (a) GW strain waveform $h(t)$ for three different combinations of radiative modes: $(l,m) = (2\pm2),(2\pm1),(3\pm3)$ without memory (red), the null memory-only $(2,0)$ mode (green), and the complete signal including both oscillatory and memory contributions (blue). (b) The corresponding response induced by the GW source located at $(\text{ra},\text{dec}) = (0^{\circ},0^{\circ})$ on a pulsar at $(\text{ra},\text{dec}) = (258.4564^{\circ},7.7937^{\circ})$. Colours match those in (a) and the merger occurs at MJD 58000. (c) Post-fit residuals after subtraction of the best-fit pulsar spin frequency and its derivatives. The full signal is shown in blue, the linear and quadratic spin-down model in orange, and the resulting post-fit residuals in green.
• Figure 2: The figure illustrates the posterior distributions from our simulation recovery studies with 25 pulsars. Simulation studies were conducted with 2 types of datasets - wn+rn+mem and wn+rn+gwb+mem and the contours are shown in blue and orange respectively. Three subfigures show a different SMBHB merger. (a) A merger with $\mathcal{M}_c = 10^{8} M_{\odot}$ and $D_\text{L} = 3$ Mpc, (b) A merger with $\mathcal{M}_c = 10^{9} M_{\odot}$ and $D_\text{L} = 20$ Mpc and, (c) A merger with $\mathcal{M}_c = 10^{10} M_{\odot}$ and $D_\text{L} = 100$ Mpc. Maximum-a-posteriori (MAP) values are indicated by dotted lines in the respective contour colors (blue and orange), and the true simulated values are shown as red dotted lines.
• Figure 3: Sky maps of the posterior distributions of simulated supermassive black hole binary mergers, shown in equatorial coordinates (RA, Dec in degrees). The color scale represents the relative posterior probability density derived from our HEALPix projections of the samples, while the yellow cross marks the true (simulated) sky location of the source. Each square projection shows a zoomed region around the source, with the scale bar indicating the angular width on the sky. The three subfigures correspond to mergers with increasing chirp mass and luminosity distance: (a) $\mathcal{M}_c = 10^8 M_\odot$, $D_L = 2$ Mpc; (b) $\mathcal{M}_c = 10^9 M_\odot$, $D_L = 20$ Mpc; (c) $\mathcal{M}_c = 10^{10} M_\odot$, $D_L = 100$ Mpc. As the chirp mass increases, the posterior uncertainty in sky localization decreases, demonstrating improved measurement precision for more massive mergers.
• Figure 4: Timing residuals from a pulsar located at (ra,dec) = ($258.4564^{\circ}$, $7.7937^{\circ}$) induced by a gravitational wave from a SMBHB merger with $\mathcal{M}_c = 10^{10} M_\odot$, $D_L = 100$ Mpc and $q=1$ located at ($0^\circ$ ,$0^\circ$). (a) Prefit timing residuals for memory burst model (orange) and the SMBHB merger model (blue). Both models are adjusted to have the same final memory offset. The SMBHB model shows a gradual pre-merger rise at angle $\alpha$, whereas the burst model shows an abrupt change at angle $\beta$. (b) Post-fit residuals after subtraction of a linear spindown model. The observable residual bump for the SMBHB model (blue) is suppressed compared to the burst model (orange), due to pre-merger curvature being partially absorbed by the fit. For the equal-offset case, the post-fit RMS difference between the models is $0.399\,\mu\mathrm{s}$. (c) Difference in post-fit timing residuals between the burst model and SMBHB model for varying mass ratios. The discrepancy is largest for the equal-mass case ($q=1$) and decreases for more unequal-mass binaries. The RMS difference drops by $\sim 0.1\,\mu\mathrm{s}$ when the mass ratio increases from $q=1$ to $q=7$.
• Figure 5: Lower limits on luminosity distance as a function of chirp mass in the signal-free simulated data. The filled contours show the marginalized posterior density in the $\mathcal{M}$–$D_L$ plane using the full SMBHB merger waveform model. The red curve represents the exclusion boundary below which a source would have been detectable with $95\%$ confidence. The orange curve shows the corresponding lower limit derived from the $95\%$ upper limit on memory strain amplitude using the memory burst model, converted to distance via Equation \ref{['eq:strain_mass_dist']}.