Semi-empirical Framework of Supermassive Black Hole Evolution: Highlighting a possible tension between Demographics and Gravitational Wave Background

TL;DR

The paper tackles the evolution of the supermassive black hole population by combining a continuity equation for mass growth via accretion with a Smoluchowski coagulation term for binary mergers, all constrained by local BH demographics, clustering, and nano-Hz GW data from pulsar timing arrays. It formulates a six-parameter Bayesian framework tied to observable inputs: the AGN luminosity function, the Eddington ratio distribution, and a merger-rate model anchored to galaxy pair counts, then derives the resulting BH mass function, clustering bias, and GW background. The results show accretion-dominated growth across most of cosmic history, while mergers contribute modestly to the high-mass end and are insufficient to account for the full PTA GW background, implying additional GW sources or revised demographics may be needed. This data-driven approach provides robust, model-independent constraints and highlights tensions between BH demographics and the GW background, with clear implications for future surveys and GW observatories.

Abstract

The evolution of the supermassive Black Hole (BH) population across cosmic times remains a central unresolved issue in modern astrophysics, due to the many noticeable uncertainties in the involved physical processes that span a huge range of spatial, temporal and energy scales. Here we tackle the problem via a semi-empirical approach with minimal assumptions and data-driven inputs. This is based on a continuity plus Smoluchowski equation framework that allows to unitarily describe the two primary modes of BH growth: gas accretion and binary mergers. Key quantities related to the latter processes are incorporated through educated parameterizations, and then constrained in a Bayesian setup from joint observational estimates of the local BH mass function, of the large-scale BH clustering, and of the nano-Hz stochastic gravitational wave (GW) background measured from Pulsar Timimg Array (PTA) experiments. We find that the BH accretion-related parameters are strongly dependent on the local BH mass function determination: higher normalizations and flatter high-mass slopes in the latter imply lower radiative efficiencies and mean Eddington ratios with a stronger redshift evolution. Additionally, the binary BH merger rate is estimated to be a fraction $\lesssim 10^{-1}$ of the galaxy merger rate derived from galaxy pairs counts by \texttt{JWST}, and constrained not to exceed the latter at $\gtrsim 2σ$. Relatedly, we highlight hints of a possible tension between current constraints on BH demographics and the interpretation of the nano-Hz GW background as predominantly caused by binary BH mergers. Specifically, we bound the latter's contribution to $\lesssim 30-50\%$ at $\sim 3σ$, suggesting that additional astrophysical/cosmological sources are needed to explain the residual part of the signal measured by PTA experiments.

Figures (14)

• Figure 1: Schematics showing the basic methodology underlying our semi-empirical framework (see Section \ref{['sec|theory']} for details).
• Figure 2: The AGN bolometric luminosity functions derived from multi-wavelength data by Shen2020 (we use the fit labeled 'A' in their Table 4; see also Section \ref{['sec|accretion']}) as a function of redshift (color-coded), which is the basic input in our semi-empirical framework.
• Figure 3: MCMC posterior distributions for the parameters exploited in our semi-empirical framework: radiative efficiency of the thin disk phase $\epsilon_{\rm thin}$; mean Eddington ratio $\langle\log \lambda_0 \rangle$ and its redshift evolution parameter $\xi$; mass-dependence parameter $\eta$ in the supermassive BH merger rate and rescaling parameter $f_{\bullet\bullet\rightarrow \bullet}$ of the BH to galaxy merger rate; scatter $\sigma_{\log M_{\rm H}}$ in the abundance matching relationship between halo and BH mass. Coloured contours/lines refer to the analysis with different determinations of the local BH mass function: orange is for that obtained combining the local SMF with the $M_\bullet-M_\star$ relation, and magenta adds a narrow prior around $f_{\bullet\bullet\rightarrow\bullet}\approx 1$; blue is for that obtained combining the VDF with the observed $M_\bullet-\sigma$ relation; green is for that obtained combining the VDF with the debiased $M_\bullet-\sigma$ relation. The marginalized distributions are in arbitrary units (normalized to unity at their maximum value). Colored crosses mark bestfit values. In the uppermost panel, the gray shaded areas highlight non-physical values of the thin-disk radiative efficiency $\epsilon_{\rm thin}$ which are excluded by hard priors.
• Figure 4: The supermassive BH mass function. Colored symbols refer to the local BH mass function determined from SMF and observed $M_\bullet-M_\star$ relation (green squares), from VDF and observed $M_\bullet-\sigma$ relation (blue circles), and from VDF and debiased $M_\bullet-\sigma$ relation (red rhomboids). These determinations bracket the overall systematics uncertainty on the estimates of the local BH mass function. For comparison, the dotted line is the recent determination based on the $M_\bullet-M_\star$ scaling relation by Liepold2024; the gray shaded area is the determination by Shen2020, obtained convolving the AGN luminosity function with the AGN luminosity vs. Eddington ratio relationship from SDSS Shen2009Nobuta2012 and further correcting for obscured AGNs and inactive BHs. Grey symbols display recent estimates of the active BH mass function at $z\sim 4-6$ including BAL quasars and little red dots by He2024 (reverse triangles), Lai2024 (triangles), Taylor2025 (pentagons), Matthee2024 (hexagons), Kokorev2024 (plus signs), and Wu2022 (crosses). Solid lines with shaded areas illustrate the outcomes (median and $2\sigma$ credible intervals) for the local BH mass function reconstructed from our semi-empirical framework and MCMC analysis, with the same color-code as in previous Figure, while dashed lines illustrate the corresponding active BH mass functions at $z\sim 4.5$.
• Figure 5: Large-scale bias as a function of BH mass. Data are from X-ray AGNs at $z\sim 0.04$ by Powell2018Allevato2021 (circles), from X-ray (squares) and optical (crosses) AGNs at $z\sim 0.25$ by Krumpe2015Krumpe2023. Colored solid lines with shaded areas illustrate the outcomes (median and $2\sigma$ credible intervals) from our semi-empirical framework model and MCMC analysis, with the same color code as in previous Figures. The inset shows the BH bias as a function of redshift, averaged over the supermassive BH mass function. Colored lines are the outcomes from our semi-empirical framework, while symbols report observational estimates from optical (plus signs; Shen2009White2012Ross2020), X-ray (diamonds; Allevato2011Allevato2014Allevato2016) and radio-selected (inverse triangles; Magliocchetti2017Mazumder2022Hale2025) samples. Note that the redshift dependence has not been fitted upon, and constitutes just an a-posteriori validation of our semi-empirical framework.