Supermassive black holes in cosmological simulations I: M_BH − M_* relation and black hole mass function

The past decade has seen significant progress in understanding galaxy formation and evolution using large-scale cosmological simulations. While these simulations produce galaxies in overall good agreement with observations, they employ different sub-grid models for galaxies and supermassive black holes (BHs). We investigate the impact of the sub-grid models on the BH mass properties of the Illustris, TNG100, TNG300, Horizon-AGN, EAGLE, and SIMBA simulations, focusing on the M_BH − M_* relation and the BH mass function. All simulations predict tight M_BH − M_* relations, and struggle to produce BHs of M_BH ≤ 10^7.5 M☉ in galaxies of M_* ∼ 10^10.5-10^11.5 M☉. While the time evolution of the mean M_BH − M_* relation is mild (ΔM_BH ≤ 1 dex for 0 ≤ z ≤ 5) for all the simulations, its linearity (shape) and normalization varies from simulation to simulation. The strength of SN feedback has a large impact on the linearity and time evolution for M_* ≤ 10^10.5 M☉. We find that the low-mass end is a good discriminant of the simulation models, and highlights the need for new observational constraints. At the high-mass end, strong AGN feedback can suppress the time evolution of the relation normalization. Compared with observations of the local Universe, we find an excess of BHs with M_BH ≥ 10⁹ M☉ in most of the simulations. The BH mass function is dominated by efficiently accreting BHs (log₁₀ f_Edd ≥ −2) at high redshifts, and transitions progressively from the high-mass to the low-mass end to be governed by inactive BHs. The transition time and the contribution of active BHs are different among the simulations, and can be used to evaluate models against observations.