The multimessenger view of Pulsar Timing Array black holes with the Horizon-AGN simulation

Problem Statement
The origin and evolution of supermassive black holes (MBHs) remain uncertain, particularly regarding their mergers in the Local Universe. While Pulsar Timing Arrays (PTAs) have recently detected a gravitational wave background (GWB) signal consistent with a stochastic superposition of massive black hole binaries (MBHBs), there is a persistent tension between the observed amplitude and theoretical predictions, which often favor lower amplitudes. Furthermore, the interpretation of the GWB spectrum is complicated by the finite duration of datasets and the potential presence of individually resolvable continuous wave (CW) sources standing above the background. Current methods for comparing theoretical models with PTA data often rely on simplified assumptions that may introduce biases. This paper aims to address these issues by utilizing the Horizon-AGN cosmological simulation to study the properties of MBHBs contributing to the PTA band, generate realistic populations for GWB and CW predictions, and critically evaluate the methodologies used to infer GWB spectral properties from timing residuals.

Methodology
The study utilizes the Horizon-AGN hydrodynamical cosmological simulation, which models a $142 \, \text{cMpc}^3$ volume with a standard $\Lambda$CDM cosmology. The authors constructed a catalog of MBH mergers by tracking MBH IDs through the simulation and identifying numerical mergers. To account for the physical delay between numerical merging (at $\sim 4 \, \text{kpc}$ separation) and actual coalescence, they applied a post-processing pipeline incorporating:

1. Dynamical Friction: Calculated based on stellar dynamical friction within an isothermal sphere approximation.
2. Binary Hardening: Modeled via stellar hardening, viscous torques in circumbinary disks, and gravitational wave (GW) emission, following established prescriptions (e.g., Sesana & Khan 2015; Dotti et al. 2015).
3. Eccentricity Evolution: Tracked via loss-cone scattering, viscous drag, and GW emission, assigning initial eccentricities based on previous numerical results.

To estimate the GWB, the authors employed two distinct approaches:

• Gaussian Ensemble (GE) Approach: A deterministic calculation of the characteristic strain spectrum assuming a continuous distribution of sources, often used in literature.
• Population Approach: A stochastic method generating 2,000 "Universe realizations" by drawing inspiralling binaries from a Poisson distribution based on the merger catalog. This accounts for cosmic variance and the discrete nature of sources.

The authors simulated PTA timing residuals for these populations using two methods:

1. Gaussian Ensemble: Converting characteristic strain to timing residuals power spectral density (PSD) via an averaged pulsar response.
2. Population Approach: Summing the individual Earth-term GW strain contributions of all binaries in a realization, explicitly modeling the finite observation duration ($T_{\text{obs}}$) and spectral leakage effects via sinc functions.

They compared these theoretical predictions against "toy" PTA datasets representing the current EPTA (25 pulsars, 10.33 years) and future SKA (80 pulsars, 10.33 years) sensitivities. They also identified potential CW candidates (individual binaries with Signal-to-Noise Ratio $> 3$) and analyzed their electromagnetic properties (optical, NIR, radio, X-ray) by modeling AGN emission and host galaxy properties.

Key Contributions and Results

• GWB Spectral Inference and Bias: The study highlights a significant methodological caveat in comparing theoretical models to PTA data. The standard Gaussian Ensemble approach, which assumes infinite observation duration, fails to account for spectral leakage caused by finite $T_{\text{obs}}$. When the population approach is used with realistic PTA analysis (including timing model marginalization), the inferred GWB spectrum is biased toward shallower slopes and higher amplitudes at $f = 1/\text{yr}$ compared to the true underlying spectrum. This bias, driven by spectral leakage from low-frequency sources, may partially explain the tension between theoretical predictions and observed PTA amplitudes.
• Stochasticity and Source Counts: At low frequencies ($\sim 3 \, \text{nHz}$), the GWB is dominated by hundreds to thousands of binaries, confirming the validity of the stochastic background approximation. However, at higher frequencies ($\sim 10\text{--}50 \, \text{nHz}$), the number of contributing sources drops significantly (to a few tens), leading to a loss of stochasticity and increased cosmic variance. Eccentric populations maintain stochasticity to higher frequencies than circular ones due to emission at higher orbital harmonics.
• Source Properties (Background vs. CW):
  • Background: The GWB is dominated by binaries with chirp masses $M_c \sim 10^{8.5} - 10^{9.5} \, M_\odot$ at redshifts $z \in 0.05 - 1$. These are primarily hosted in massive galaxies ($M_\star > 10^{11} \, M_\odot$) within halos of mass $\sim 10^{13} \, M_\odot$.
  • CW Candidates: Individual resolvable sources (CW candidates) tend to be more massive and located at lower redshifts ($z < 0.4$) compared to the general background. They are predominantly hosted in the central galaxies of massive groups and clusters ($M_{\text{halo}} > 10^{13} \, M_\odot$).
• Electromagnetic Counterparts:
  • Nature: The majority of CW candidates appear as Active Galactic Nuclei (AGN) rather than quasars, characterized by low Eddington ratios ($f_{\text{Edd}} < 10^{-2}$).
  • Detectability: In optical and near-infrared (NIR) bands, most CW candidates are outshined by their host galaxies, making light-curve variability searches challenging and biased toward rare, luminous systems.
  • Best Wavelengths: Radio and X-ray bands offer the best prospects for identification. For candidates with $f_{\text{Edd}} > 10^{-3}$, the AGN emission can outshine the host galaxy in these bands. Even for low-accretion systems, optimistic estimates suggest up to 80% (radio) and 30% (X-ray) could be detectable above the host.
• Detection Probabilities: Assuming a 10-year baseline, the probability of detecting at least one CW candidate with $S/N > 3$ is estimated at 4% for EPTA-like sensitivities and 20% for SKA-like sensitivities.

Significance
The paper claims that its primary significance lies in providing a unified, multi-messenger framework for understanding the PTA signal using a state-of-the-art cosmological simulation. Crucially, it challenges the standard methodology for comparing theoretical GWB models with observational data. By demonstrating that finite observation duration and spectral leakage can artificially inflate the inferred GWB amplitude and flatten the spectral slope, the authors suggest that the apparent tension between astrophysical models and PTA observations may be partially an artifact of analysis techniques rather than a failure of the underlying astrophysics. Additionally, the study provides specific guidance for electromagnetic follow-up, suggesting that searches should prioritize radio and X-ray bands and focus on massive galaxy groups and clusters, rather than relying solely on optical variability in quasar-like systems.