The multi-messenger view of pulsar timing array black holes with the Horizon-AGN simulation

¹ Université Paris Cité, CNRS, Astroparticule et Cosmologie, F-75013 Paris, France
² Institut d’Astrophysique de Paris, UMR 7095, CNRS and Sorbonne Université, 98 bis boulevard Arago, 75014 Paris, France
³ IRAP, Université de Toulouse, CNRS, UPS, CNES, 9 Avenue du Colonel Roche, BP 44346, 31028 Toulouse Cedex 4, France
⁴ Institute for Astronomy, University of Edinburgh, Royal Observatory, Edinburgh EH9 3HJ, UK

^★ Corresponding author: This email address is being protected from spambots. You need JavaScript enabled to view it.

Received: 16 May 2025
Accepted: 3 November 2025

Abstract

We used the HORIZON-AGN cosmological simulation to study the properties of supermassive black hole binaries (MBHBs) with the largest contribution to the gravitational wave background (GWB) signal expected for the pulsar timing array (PTA) band. We developed a pipeline to generate realistic populations of MBHBs, which enabled us to estimate both the characteristic strain and GWB time series observable by PTA experiments. We identified potential continuous wave (CW) candidates standing above the background noise, using toy PTA sensitivities representing the current EPTA and future SKA. We estimated the probability of detecting at least one CW with a signal-to-noise ratio of S/N > 3 to be 4% (20%) for EPTA (SKA)-like sensitivities, assuming a ten-year baseline. We found that the GWB is dominated by hundreds to thousands of binaries at redshifts in the range 0.05 − 1, with chirp masses of 10^8.5 − 10^9.5 M_⊙, primarily hosted in quiescent massive galaxies residing in halos of mass ∼10¹³ M_⊙. CW candidates have larger masses and lower redshifts, and they tend to be found in even more massive halos, typical of galaxy groups and clusters. The majority of these systems would appear as active galactic nuclei (AGNs), rather than quasars, because of their low Eddington ratios. Nevertheless, CW candidates with f_Edd > 10⁻³ can still outshine their hosts, particularly in radio and X-ray bands, suggesting that they could serve as the most promising route for identification. Our findings imply that optical and near-infrared (NIR) searches based on light curve variability are challenging and biased toward more luminous systems. Finally, we highlight important caveats in the common method used to compare PTA observations with theoretical models. We find that GWB spectral inferences used by PTAs could be biased toward shallower slopes and higher amplitudes at f = 1/yr, thereby reducing the apparent tension between astrophysical expectations and PTA observations.