Multimessenger Probes of the Supermassive Black Hole Binary Population: The Role of Pulsar Timing Arrays

TL;DR

This work addresses how pulsar timing array measurements of the nanohertz gravitational-wave background can inform the demographics and dynamical evolution of the SMBHB population beyond electromagnetic constraints. Using a six-parameter phenom population model and BayesianInference on simulated PTA data (A4cast and SimGWB) with normalizing-flow emulators, the authors quantify information gain and prior dominance for SMBHB parameters. They find that PTA data meaningfully constrains the total lifetime $\tau_f$, the hardening rate $\nu_{\text{inner}}$, and the GSMF characteristic mass $m_{\psi,0}$, but remains largely prior-dominated for $M_{\rm BH}$--$M_{\rm Bulge}$ relation parameters ($\mu$, $\varepsilon_\mu$) and the GSMF normalization $\psi_0$ in realistic data regimes. The results highlight the selective power of PTAs in probing SMBHB evolution and emphasize the need for multimessenger approaches to fully characterize the SMBHB population and its astrophysical drivers.

Abstract

By inferring the gravitational wave background (GWB) from a population of supermassive black hole binaries (SMBHBs), pulsar timing arrays (PTAs) enable the study of massive black holes. In many ways, PTAs manifest the promise of a multimessenger approach to astronomy: they can constrain SMBHB population characteristics that are otherwise difficult to constrain using electromagnetic observations, such as hardening mechanisms at sub-parsec separations. In this work, we quantify this multimessenger promise using Bayesian inference of many realizations of simulated PTA data, while adopting a model for the SMBHBs that has been successfully applied to the 15-year data set of the North American Nanohertz Observatory for Gravitational Waves (NANOGrav). Our analyses of 200 realistic, simulated NANOGrav data sets show that there is a greater than 50\% chance of reducing the prior uncertainty in the SMBHB hardening rate by more than 50\%, and in the SMBHB evolutionary lifetime by 25--75\%. Additionally, there is an 88\% chance that PTA data can reduce the prior uncertainty in the characteristic mass variable of the galaxy stellar mass function (GSMF) by 25--50\%. For $M_{\text{BH}}$--$M_{\text{Bulge}}$ parameters (in a model without redshift evolution) and the overall normalization parameter of the GSMF, PTA data can provide only marginal information gain beyond the constraints from electromagnetic observations. Our work delineates the domains over which electromagnetic and gravitational-wave data constrain the demographics and dynamics of the supermassive black-hole binary population, offering a clearer picture of the impact of population multi-messenger astrophysics probes with PTAs.

Figures (6)

• Figure 1: A set of contours and histograms showing the multivariate distribution of the GWB-induced PSD of the timing residuals that is used for simulations in this work. The solid blue contours highlight the distribution, while the gray lines indicate the 200 samples drawn from it, which are used to inject into the simulated data sets A4cast and simGWB. As shown in the figure, the 200 drawn samples cover most of the distribution.
• Figure 2: A set of plots depicting the Bayesian inference posterior distributions for $\bm{\theta}$ obtained by analyzing all realizations of A4cast (orange squares) and SimGWB (blue circles). The error bars range over the 68% HDI, whereas the circles and squares show the MAP values. Each panel of the figure belongs to one parameter in $\bm{\theta}$. Attached to each panel, the prior probability distribution and its 68% HDI interval (dashed horizontal lines) are shown. The actual values of each parameter used in all simulations are shown with a dashed line across the main panels. The main takeaway messages of the plots are (i) the prior-dominant posterior distributions for the $M_{\text{BH}}$--$M_{\text{Bulge}}$ parameters regardless of the data set and (ii) informative posterior distributions for the parameters that cannot be easily constrained using electromagnetic observations (i.e., $\tau_f$ and $\nu_{\rm{inner}}$).
• Figure 3: A series of plots showing the relationship between the black hole mass function $\Phi$ and the mass of black holes. The label of each subplot indicates the phenom model parameter that varies over its prior distribution samples, while the rest of the parameters are fixed to their values in \ref{['astro-chosen']}. Each thin band in the plot corresponds to a unique draw from the relevant prior distribution. An arbitrary redshift of 1 is chosen for the plots. Note that $\tau_f$ and $\nu_{\rm{inner}}$ do not affect the black hole mass function; hence, they are not relevant for this plot. Compared to the other two parameters, changes in $\psi_0$ and $m_{\psi,0}$ over their prior distribution correspond to a larger change in the black hole mass function both at low and high black hole mass values.
• Figure 4: Heatmap plots depicting the reduction in uncertainty (68% HDI) of the posterior relative to the prior (left panel) and the posterior-to-prior probability ratio for the small region around the actual values (values of \ref{['astro-chosen']}$\pm$1%) (right panel). The numbers within each cell indicate the percentage of total realizations that fall within a range of the estimated quantity in each panel. The two data sets are color-coded: different shades of blue (lower triangular region of each cell) correspond to simGWB, and different shades of orange (upper triangular region of each cell) correspond to A4Cast. For instance, for $\tau_f$, 50.5% of the realizations see a reduction in uncertainty between 25 to 50% if the inferences are done on simGWB, while this number is 52% for A4Cast data. Overall, the figure quantifies the extent of information gain after incorporating the PTA data.
• Figure 5: A heatmap plot depicting the range of random forest regression scores for each parameter at each frequency bin. The different shades of blue provide a visual clue for the score: the darker the blue, the larger the score. Completely independent decision trees are constructed for each frequency bin; thus, scores are not comparable across frequency bins. The figure indicates that $\mu$ and $\varepsilon_{\mu}$ have negligible influence on the changes of in the amplitude of the GWB-induced timing residual PSD, in contrast to the strong impact of $m_{\psi,0}$, $\nu_{\rm{inner}}$, and $\tau_f$.