Greene 2012; Kormendy & Ho 2013). These objects were likely already in place early on in the history of our Universe. Evidence for the presence of these massive objects in the early Universe (z > 5) include observations of extremely bright quasars powered by 10 9-10 M BHs when the Universe was about 1 Gyr (Mortlock et al. 2011;Bañados et al. 2018), but also relativistic jets and centrally driven winds at various redshifts (down to z = 0; Förster Schreiber et al. 2014;Genzel et al. 2014;Cheung et al. 2016). Accreting BHs, On the population of BHs 1941 known as active galactic nuclei (AGNs) and observed in a broad range of redshifts, can release large amounts of energy in their host galaxies, via AGN feedback. One of the major outstanding issues in modern astrophysics is how BHs form and evolve from high redshift until today, and how they interact with their host galaxies.

Over the last decade, we have been able to numerically address galaxy formation in a cosmological context. Large-scale cosmological simulations with a box length of ∼ 100 cMpc and a spatial resolution of ∼ 1 kpc have successfully demonstrated that it was possible to achieve reasonable agreements with current observational constraints, in terms of e.g. galaxy clustering, stellar mass content of galaxies/haloes, galaxy star formation rates, morphologies, sizes, colour bimodality (e.g. Di Matteo et al. 2008;Ocvirk, Pichon & Teyssier 2008;Dubois et al. 2014a;Genel et al. 2014;Hirschmann et al. 2014;Vogelsberger et al. 2014a, b;Schaye et al. 2015;Davé, Thompson & Hopkins 2016;Pillepich et al. 2018a;Davé et al. 2019).

Some discrepancies among galaxy properties have been found as well, and they help us to improve the simulation sub-grid physics models. For example, the galaxy sizes in Illustris were too large (Nelson et al. 2015;Pillepich et al. 2018a) compared to observational constraints. The galaxy mass function of Illustris showed an excess of low-mass galaxies with M 10 10 M (Pillepich et al. 2018a). Horizon-AGN also produced an excess of galaxies with M 10 10 M at z = 0, before the knee of the galaxy mass function. The gas fractions in the Illustris galaxies were elevated compared to observations, while the halo gas fractions were too low. This was a consequence of the hot bubble mode of the AGN feedback of Illustris, whose injection of thermal energy was displaced from the galaxies, and impacting more the host haloes. As a result, the stellar-to-halo mass ratios of Illustris were on the upper limit of observational constraints, and too high for massive haloes of M h a few 10 12 M (Vogelsberger et al. 2014b). Most cosmological simulations (Horizon-AGN, EAGLE, MUFASA, Illustris) present lower fractions of quenched galaxies (i.e. galaxies with suppressed star formation rates) for M 10 10 M compared to SDSS constraints (e.g. Donnari et al. 2019;Hahn et al. 2019). Many papers address comparisons between these simulations and observations, but they often use different observational constraints and definitions of the studied quantities; drawing global conclusions from these analyses is difficult.

The resolution of large-scale cosmological simulations is not sufficient to resolve processes across the entire dynamical range needed, from BH accretion discs to large-scale filaments. Processes related to BH formation, growth, and feedback, as well as any other baryonic processes taking place at the galactic scale are modelled with sub-grid physics. Simulations all use different sub-grid models, e.g. different location and mass for BH seeding, different models to compute the accretion on to BHs, different efficiencies and models to inject energy released by AGN, some models assuming that AGN feedback channels explicitly depend on BH mass, some others assuming a uniform feedback. Therefore, it is crucial to understand whether the BH populations produced by these simulations are in good agreement with observational constraints, and how the discrepancies could help us to improve sub-grid models.

We can generally divide the population of BHs in three categories. The first one is the population of high-redshift quasars, i.e. the most extremely bright and massive BHs that we observed at high redshift, in optical and near-IR surveys. These monsters probe the most extreme end of the BH mass spectrum, and are a challenge for our understanding of BH formation/accretion: they would require an almost continuous growth at the Eddington limit for the entire life of these objects. On the other end of the BH mass spectrum, we find the tiny BHs of M BH 10 4 M found in dwarf galaxies of M 10 9 M (Reines & Comastri 2016;Greene, Strader & Ho 2019), and theoretically expected to be the relics of BH seeds formed at high redshift. These two first regimes are unfortunately not covered by the large-scale cosmological simulations studied here. Dwarf galaxies are indeed barely resolved and BH formation subgrid models are too simplistic. These simulations are also limited by their volumes to form and evolve the rare quasars (but see Di Matteo et al. 2017, for a larger simulation of 400 h -1 Mpc run down to z = 7). Between these two BH mass ranges, cosmological simulations fully cover the regime of all the other BHs with mass M BH = 10 5 -10 10 M observed at z = 4 -0. When accreting matter, these BHs power the AGN that we see e.g. in X-ray surveys (e.g. XMM-Newton, Chandra, NuSTAR) with X-ray luminosities of L x = 10 42 -10 46 erg s -1 (Ueda et al. 2014;Aird et al. 2015;Buchner et al. 2015;Miyaji et al. 2015, and references therein). This is the category of BHs that we investigate in this series of papers.

With time, we have accumulated evidence for correlations between the mass of BHs and the properties of their host galaxies, such as the total stellar mass, bulge mass, velocity dispersion, Sérsic index, and infrared luminosity (Marconi & Hunt 2003;Häring & Rix 2004;Merloni 2004;Shankar et al. 2004;Gültekin et al. 2009;Shankar, Weinberg & Miralda-Escudé 2009b;Kelly & Merloni 2012;Kormendy & Ho 2013;McConnell & Ma 2013;Läsker et al. 2014;Graham & Scott 2015;Reines & Volonteri 2015). These scaling relations are mostly observed in the local Universe since it is extremely challenging to measure these quantities at higher redshifts. For this reason, we only study the M BH -M relation in this paper (M being easier to measure), even if BH mass scales more tightly with the galaxy velocity dispersion locally (Saglia et al. 2016;Shankar et al. 2016;de Nicola, Marconi & Longo 2019;Marsden et al. 2020). The scaling relations potentially imply the co-evolution of BHs and their host galaxies. However, several works have also emphasized the possibility that scaling relations may not require co-evolution processes according to the central-limit theorem (Hirschmann et al. 2010;Jahnke & Macciò 2011). In these studies, scaling relations can be explained by the hierarchical assembly of BH and bulge growth by galaxy mergers. A large scatter in the high-redshift Universe will be reduced with time. Assuming that BHs only grow by mergers, Hirschmann et al. (2010) show that the scatter of the M BH -M bulge relation decreases with an increasing number of mergers, and so with time.

The time evolution of such scaling relations is difficult to investigate in high-redshift observations. Studies are often based on broad-line AGN targets, where BH masses are estimated from the virial method. The underlying assumption is that broad-line AGN behave the same way as any other galaxies, and obey the same scaling relations, which may not be the case (Reines & Volonteri 2015). Most observational studies (including sometimes fainter AGN as well) have found mild or null evolution of the scaling relations (Shields et al. 2003;Jahnke et al. 2009;Salviander & Shields 2013;Schramm & Silverman 2013;Shen et al. 2015). For instance, Salviander & Shields (2013) find no time evolution of the M BHσ relation of quasars for 0.1 < z < 1.2 and for 10 7.5 < M BH /M < 10 9 ; the relation evolves for lower and more massive BHs, but this can be attributed to observational uncertainties, as well as the large intrinsic scatter in the relation. Similar results are obtained in Cisternas et al. (2011) for 32 broad-line AGNs with 10 7.2 < M BH /M < 10 8.7 in the redshift range 0.3 < z < 0.9, and in Schramm & Silverman (2013) with 18 X-ray selected broad-line AGNs of the CDFS in the range 0.5 < z < 1.2. Jahnke et al. (2009) find no evolution of the M BH -M relation while comparing 1 < z < 2 objects with the local M BH -M bulge relation.

The evolution of the scaling relation is often studied by comparing the high-redshift M BH -M to the same relation scaled from the local bulge mass relation. Going from the M BH -M plane to the M BH -M bulge plane implies a redistribution of the stellar mass from the galactic disc to form a bulge, by mergers or secular processes. The studies above find that a large fraction of the high-redshift galaxies have, indeed, a disc component. The mild evolution found in observations tends to indicate that no addition of stellar mass would be required to build the bulges, and the M BH -M bulge relation, from the high-redshift galaxies. The presence of discs in the highredshift galaxies also means that their bulges can be considered as undermassive, which implies that the time evolution of the M BH -M bulge relation may be stronger than the evolution of the M BH -M relation. It also favours the idea that the assembly of BHs takes place before the assembly of their host galaxy bulges. More recently, Sun et al. (2015) (70 Herschel-detected broad-line AGN) found no evidence for redshift evolution of the M BH -M relation since z = 2. Mullaney et al. (2012) reached the same conclusion based on an X-ray stacking study (z 2), as well as Suh et al. (2020) with a sample of 100 X-ray selected moderately luminous AGN at z 2.5, combined with brighter AGN from the literature. However, recently Ding et al. (2020) used 32 X-ray selected broad-line AGNs within 1.2 < z < 1.7 and showed that the ratio M BH /M and M BH /M bulge were larger than in the local Universe (∼ 0.43 dex), meaning that the relations are slightly offset towards more massive BHs or lower stellar mass/bulge mass compared to the relation at z = 0. Finally, McLure et al. (2006) find that the ratio M BH /Mfoot_1 of radioloud AGN galaxies increases with redshift. Since their sample is based on the most massive early-type galaxies in the redshift range 0 < z < 2, we cannot rule out that the results are biased. Nevertheless, it may also indicate that BH growth in these massive galaxies is completed first, and then galaxies catch up. In this paper, we will see that there is no consensus on the time evolution of the M BH -M mean relation in cosmological simulations.

Several observational studies have pointed out that the M BH -M relation may not be linear over the full range of stellar masses. The shape of the relation could depend e.g. on BH mass, stellar mass, galaxy types, star-forming activity. The correlation between the quiescence of the host galaxies and the shape of the M BH -M relation has been discussed in e.g. Graham & Scott (2015). Different slopes of the relation have been found for samples of earlytype/core-Sersic and late-type/Sersic galaxies in observations (Davis, Graham & Cameron 2018;Sahu, Graham & Davis 2019), but also for different types of galaxies and activity of the BHs (i.e. AGN or quiescent BHs; Reines & Volonteri 2015). Recently, a non-linear shape of the M BH -M relation was also found in cosmological simulations (e.g. Bower et al. 2017;Habouzit, Volonteri & Dubois 2017;Anglés-Alcázar et al. 2017c) and semi-analytical galaxy formation models (e.g. Fontanot, Monaco & Shankar 2015). In the following, we compare the shape and normalization of the M BH -M relation in several cosmological simulations to understand how this relation is impacted by the simulations subgrid physics.

In addition to scaling relations between BH mass and galaxy properties, the BH mass function describes the evolution of the mass distribution of BHs (Kelly & Merloni 2012, for a review), and is crucial to understand BH mass growth and the effects of BH accretion and feedback. The comparison with observational constraints can be tricky as these constraints often rely on some physical modelling/assumptions. The empirical BH mass functions have been used to constrain physical models of BH evolution. However, due to the uncertainties, it is actually challenging to rule out a significant number of models. Estimating the BH mass function suffers from incompleteness and mass estimations. BH mass measurements are only possible for a small subset of local BHs. Therefore, these analyses rely on scaling relations (with bulge mass, bulge luminosity, width of the broad emission lines) to estimate BH mass for a large number of BHs, and not direct measurements of BH mass. In this case, the local BH mass function is estimated by combining one of the scaling relations with the local number density of galaxies (as a function of their velocity dispersion or bulge luminosity). This method is employed, for example, in Merloni & Heinz (2008), Shankar et al. (2009). Several works have also investigated the evolution of the BH mass function with redshift, by combining the local BH mass function to the AGN luminosity function (e.g. Shankar et al. 2009;Cao 2010).

Empirical BH mass functions (Merloni & Heinz 2008;Shankar et al. 2009;Cao 2010) and theoretical models based on semianalytical models (Shen 2009;Volonteri & Begelman 2010;Fanidakis et al. 2011) or hydrodynamical simulations (Di Matteo et al. 2008;Hopkins et al. 2008) agree well at z = 0, but show large discrepancies at larger redshifts (e.g. z = 2; Kelly & Merloni 2012). These differences particularly emerge at the low-mass end of the BH mass function for M BH 10 8.5 M . Some observational constraints predict a turnover for lower mass BHs (Merloni & Heinz 2008) when some others do not (Shankar et al. 2009;Cao 2010). The contribution of active and inactive BHs to the BH mass function is also crucial to understand the build-up of the BH population with time. In this paper, we compare the BH mass functions of the different simulations, and analyse how they evolve with time, and the relative contribution of active and inactive BHs.

Over the last decade, we have pushed the development of largescale cosmological simulations to understand galaxy formation in a cosmological context. These simulations were mainly meant to reproduce galaxy properties measured in observations, and not so much the BH properties. Observational constraints become more and more accurate and numerous, so today is an excellent time to review the BH populations in cosmological simulations and how they compare to observations. In this first paper of a series, we aim at drawing conclusions on the BH populations from the six large-scale cosmological simulations: Illustris, TNG100, TNG300, Horizon-AGN, EAGLE, and SIMBA. We compare these simulations to one another in a coherent way, i.e. applying the same post-processing analyses to all the simulations. We also compare their results with observational constraints, when possible. In this paper, we do not make close apple-to-apple comparisons of observations with simulated mock observations; our comparisons with observations are mostly to guide the eye of the readers. We point out the different properties that are successfully reproduced and the properties that are not. Our main goal is to understand how the different sub-grid models of these simulations can lead to differences in the BH populations, and to evaluate the range of possible values produced by the simulations for a given quantity. This will help us to understand what we are missing in our modelling of galaxy and BH formation/evolution in large-scale cosmological simulations, but also what we need to focus on and improve in the future. In this first paper, we address the mass properties of the simulated BH populations, i.e. the M BH -M scaling relation and the BH mass function.

In Section 2, we describe the six Illustris, TNG100, TNG300, Horizon-AGN, EAGLE, and SIMBA large-scale simulations. We present the M BH -M diagram in Section 3, and analyse the scaling relations of the simulations in Section 4. In Section 5, we focus on the Illustris and TNG100 simulations, which share the same initial conditions and are, therefore, ideal to understand how the different sub-grid models create different features in the M BH -M relation. Finally, in Section 6 we present the BH mass function for all the simulations. We discuss our results and conclude in Sections 7 and 8.

In this section, we introduce the six large-scale cosmological simulations. Only BH-related subgrid models and the main aspects of SN feedback modelling are described. The analysis of the BH population of the Illustris simulation has been presented in Sijacki et al. (2015), of the TNG simulations in Weinberger et al. (2017Weinberger et al. ( , 2018)), Habouzit et al. (2019), and references therein. The analysis of the Horizon-AGN BH population has been conducted in Volonteri et al. (2016), and the analysis of EAGLE in Schaye et al. (2015), Rosas-Guevara et al. (2016), McAlpine et al. (2017, 2018). Finally, the BH population of SIMBA has been studied in Davé et al. (2019), Thomas et al. (2019). In the following, we apply the same cut in the total stellar mass of galaxies, i.e. we only consider galaxies with M 5 × 10 8 M , which is a common minimum galaxy stellar mass resolved in all simulations.

The Illustris simulation consists of a volume of (106 cMpc) 3 , and was run with the moving Voronoi mesh code AREPO (Springel 2010). The gravitational softening for DM particles is 1.4 ckpc (comoving kpc), and the collisionless baryonic particle softening is 1.4 ckpc for z 1, and 0.7 pkpc otherwise.

The BH seeding of the simulation is based on dark matter halo mass. Every halo reaching a mass of M h,thres = 7.10 × 10 10 M is seeded with a BH of initial mass M seed = 1.42 × 10 5 M . The accretion on to BHs follows the Bondi-Hoyle-Lyttleton formalism and is capped at the Eddington limit:

with cs the kernel weighted ambient sound speed, r the radiative efficiency, and α the boost factor. The boost factor is introduced to compensate for the fact that the ISM multiphase structure is not resolved (Springel et al. 2005;Booth & Schaye 2009). As simulations tend to underestimate the density around the BHs, the boost factor allows us to increase the accretion rates on to the BHs. Illustris employs a two-mode AGN feedback, both with release of thermal energy. AGN are able to deposit thermal energy in their surroundings with a coupling efficiency of 0.05 in the highaccretion mode, i.e. for BHs with high Eddington ratios of f Edd = ṀBH / ṀEdd > 0.05. For BHs with f Edd < 0.05, thermal energy is released in hot bubbles within a radius of ∼100 kpc around the BHs and couple to the gas with an efficiency of 0.35.

Illustris employs a kinetic wind model for the SN feedback, which is characterized by a mass loading factor and an initial wind velocity. The initial wind velocity can be written as

with k wind = 3.7 a dimensionless efficiency factor, and σ DM the 1D dark matter velocity dispersion.

The mass loading factor is defined by

with the available SN II energy per formed stellar mass e wind = 1.89 × 10 -foot_2 × 10 51 = 1.09 × 1.73 × 10 -2 × E SN II , with 1.73 × 10 -2 the number of SN II per stellar mass formed, and E SN II = 10 51 erg M -1 the available energy available for each SN II .

The simulations IllustrisTNG (hereafter TNG100 and TNG300; Marinacci et al. 2018;Naiman et al. 2018;Nelson et al. 2018;Pillepich et al. 2018b;Springel et al. 2018) have a volume of (111 cMpc) 3 and (302 cMpc) 3 , respectively. TNG100 shares its initial conditions with the previous Illustris simulation. The softening of the stellar particles and of DM particles is the same ( DM = 1.48, 2.96 ckpc for TNG100 and TNG300, respectively) up to z = 1, and fixed at their z = 1 proper values for later redshifts ( DM,z=0 = 0.74, 1.48 kpc for TNG100 and TNG300, respectively). The minimum softening of the gas is gas,min = 0.19, 0.37 ckpc for TNG100 and TNG300. Dark matter haloes with a mass exceeding M h,thres = 7.38 × 10 10 M are seeded in their centre with BHs. The initial BH mass is set to M seed = 1.18 × 10 6 M , one order of magnitude higher than in the previous Illustris simulation.

The accretion on to BHs also follows the Bondi-Hoyle-Lyttleton model, but no boost factor is added. The TNG100 and TNG300 are magnetohydrodynamical simulations; therefore the kernel weighted ambient sound speed is now written as cs = (c 2 s,therm + (B 2 /4πρ)) 1/2 , and includes a term for the magnetic signal propagation speed. The addition of the magnetic fields can impact the relationship between the BHs and the properties of their host galaxies. In TNG, the normalization of the M BH -M mean relation is higher with magnetic fields (Pillepich et al. 2018a). Illustris and TNG have two different technical implementations of the Bondi model: Illustris only uses the parent gas cell of the BHs to compute the accretion rates on to the BHs, whereas TNG100 employs a kernel-weighted accretion rate over about 256 neighbouring cells. As a consequence, the accretion rates on to the TNG BHs correlate with the gas properties of the central region of the galaxies, while the rates in the Illustris depend on the gas properties at the location of the BHs and could be more stochastic.

TNG includes a two mode AGN feedback: injection of thermal energy in the surroundings of the BHs accreting at high rates, and directional injection of momentum, with each event oriented in a random direction 2 for low accretion rate BHs (Weinberger et al. 2017(Weinberger et al. , 2018)). The transition between modes is set by the threshold:

Only the TNG and SIMBA simulations include a dependence on BH mass for the transition between AGN feedback modes.

The stellar and gas properties in the TNG simulations have been shown to be in generally better agreement with many observations than the Illustris simulation (Pillepich et al. 2018a, b). In particular, the gas fraction was too high in galaxies and too small in the circumgalactic medium in massive haloes. These discrepancies can be attributed to the low-accretion AGN feedback mode of Illustris, for which thermal energy was injected in bubbles displaced from the centre of galaxies (Sijacki et al. 2015). Galaxies were also too large in Illustris by a factor of a few for galaxies of M 10 10.7 M , while the mass-size relation in TNG100 is now in better agreement (Genel et al. 2018).

In the following, we describe the modelling of SN feedback in TNG (both with the velocity of the winds and the wind mass loading factor), its evolution with redshift and halo mass, and how it compares to the previous Illustris SN feedback model. Our description below only discusses the strength of the SN feedback at time of injection of the winds, and do not mean that the final impact of the winds will follow the same trends. 3 The feedback of the SNe will have an important impact on the evolution of the BH populations, as demonstrated later in this paper.

The velocity of the winds in TNG100 is written as

with k wind = 7.4 (k wind = 3.7, for Illustris), σ DM the 1D dark matter velocity dispersion, H 0 , H(z) the Hubble constants, and v wind,min = 350 km s -1 the wind velocity floor imposed in TNG.

Compared to the Illustris model, TNG100 has a new dependence on the Hubble constant, as well as the addition of a velocity floor (Pillepich et al. 2018a). This overall parametrization implies that at injection, the TNG winds are faster (higher wind velocity), for all halo masses and at all redshifts. The addition of v w, min prevents low-mass galaxies to have a unphysically low velocity of the winds, and makes the feedback in these galaxies more impactful than in the Illustris galaxies. In conclusion, SN feedback in TNG is globally more impactful in low-mass galaxies at all times.

In TNG, a given fraction τ wind of the energy is released as thermal energy, and the mass loading factor is defined by

with [1.4-5.6] an efficiency factor which depends on the metallicity of the gas cell (1.4 for gas cell metallicity of 0.02): The metallicity dependence makes the SN wind mass loading factor (and the available wind energy) smaller in metal-enriched environments. The overall effect is a more effective SN feedback at injection in low-mass galaxies in TNG than in Illustris. The mass loading factor at injection is smaller in TNG than in Illustris at z = 0 for all halo masses. The mass loading factor always increases with decreasing redshift in Illustris, and is constant (for the smallest haloes of M h 10 10 M ) or decreasing (for more massive haloes) in TNG. Consequently, the wind mass loading factor of TNG is higher than in Illustris only at high redshift (z 2 -3, depending on halo mass) for haloes of M h 10 11 M .

We discussed here the strength of the feedback at injection. The final effect seen in the TNG simulation is an overall stronger SN feedback in low-mass galaxies. This is shown by the lower normalization of the galaxy stellar mass function in TNG compared to Illustris (Fig. 10) in the low-mass regime (see also Pillepich et al. 2018a). In this paper, we will analyse the effect of the modelling of SN feedback on the growth of BHs.

Horizon-AGN (Dubois et al. 2014a(Dubois et al. , 2016) ) has a volume of (142 cMpc) 3 , and was run with the adaptive mesh refinement (AMR) code RAMSES (Teyssier 2002). The seeding of Horizon-AGN is different from the precedent simulations, which all use a fixed threshold in the dark matter halo mass. In Horizon-AGN, BHs form in dense gas cells (i.e. in cells for which the gas density exceeds the threshold for star formation, here n 0 = 0.1 H cm -3 ) with stellar velocity dispersion greater than 100 km s -1 . BH seeds are not formed closer than 50 ckpc of an existing BH. Finally, BHs are prevented from forming after z = 1.5. At that time, all the progenitors of the M 10 10 M galaxies at z = 0 are formed and seeded with BHs (Volonteri et al. 2016).

The accretion rate on to BHs is computed following the Bondi-Hoyle-Lyttleton formalism with a boost factor α = (ρ/ρ 0 ) 2 when the density ρ is higher than the resolution-dependent threshold ρ 0 and α = 1 otherwise (Booth & Schaye 2009).

Horizon-AGN includes a two mode AGN feedback. In the highaccretion mode (f Edd > 0.01), thermal energy is isotropically released within a sphere of radius a few resolution elements. The energy deposition rate is ĖAGN = f r ṀBH c 2 = 0.015 ṀBH c 2 . In the lowaccretion mode, energy is injected into a bipolar outflow with a velocity of 10 4 km s -1 , to mimic the formation of a jet. The energy rate in this mode is ĖAGN = 0.1 ṀBH c 2 . The technical details of BH formation, growth, and feedback modelling of Horizon-AGN can be found in Dubois et al. (2012).

Horizon-AGN employs a kinetic SN feedback, including momentum, mechanical energy, and metals from type II, Type Ia SNe, and stellar winds (details in Kaviraj et al. 2017). The feedback is modelled as kinetic release of energy on time-scale below 50 Myr, and a thermal energy after 50 Myr. The feedback is also pulsed: energy is accumulated until sufficient to propagate the blast wave to at least two cells. The energy released depends on the SSP modelled asssumed and the metallicity of the gas, and is about e SN ∼ 10 49 erg M -1 .

The EAGLE simulation has a volume of (100 cMpc) 3 (Crain et al. 2015;Schaye et al. 2015;McAlpine et al. 2016), and was run with the code ANARCHY (Dalla Vecchia et al., in preparation;Schaller et al. 2015), which is based on the Smoothed Particle Hydrodynamics (SPH) method and is a modified version of the code GADGET3 (Springel 2005). The DM mass resolution is M DM,res = 9.7 × 10 6 M , and the baryonic mass resolution M baryon,res = 1.81 × 10 6 M . The gravitational interaction between particles is determined by a Plummer potential with a softening length of 2.66 ckpc and limited to a maximum physical length of 0.70 pkpc All DM haloes more massive than M DM = 1.48 × 10 10 M are seeded with a BH of M seed = 1.48 × 10 5 M .

Accretion rates on to BHs are computed using a modified Bondi-Hoyle-Lyttleton model (Rosas-Guevara et al. 2015, 2016):

with c s the sound speed of the surrounding gas, V the SPH-average circular speed of the gas around BHs, C visc is a free parameter referring to the viscosity of the sub-grid accretion disc. In the default Bondi model, i.e. in the absence of angular momentum of the gas, the gas within the Bondi radius (= GM BH /c 2 s ) directly falls on to the BHs. In the presence of angular momentum, the gas will instead form an accretion disc. Taking into account the angular momentum of the infalling material reduces BH accretion rates in small galaxies compared to the default Bondi model. No additional boost factor is present in the accretion rate model.

The simulation employs a single-mode AGN feedback (Booth & Schaye 2009). A fraction of the accreted gas is stochastically injected as thermal energy with a net efficiency of r f = 0.1 × 0.15 = 0.015. To prevent AGN feedback from becoming inefficient due to numerical loses and overcooling, the injection of thermal energy only occurs when the BH has accreted a sufficient amount of mass whose equivalent energy can raise the temperature of a gas particle by T = 10 8.5 K. Feedback from SNe is injected stochastically as thermal energy (Dalla Vecchia & Schaye 2012). SN energy of 10 51 erg M -1 is released 30 Myr after the stellar particle is born. The available energy injected depends on the local gas metallicity and density, and therefore, the feedback is stronger at higher redshift. The metallicity dependence is physically motivated: the efficiency is reduced for metallicities Z > 0.1Z where metal-line cooling becomes important. The efficiency increases with higher densities to compensate for numerical overcooling (the galaxies would be too compact without the density dependence; Crain et al. 2015).

The simulation SIMBA (Davé et al. 2019) was run with the code GIZMO (Hopkins 2015(Hopkins , 2017)), with a box side length of 147 cMpc. The DM mass resolution is M DM,res = 9.6 × 10 7 M , and the baryonic mass resolution M baryon,DM = 1.8 × 10 7 M . The softening is adaptive for all particle types with minimum value of = 0.74 pkpc.

A friends-of-friends algorithm identifies galaxies on the fly, and seed them with a BH of M seed = 1.4 × 10 4 M when their mass exceed the limit M 10 9.5 M . This relatively high stellar mass threshold is motivated by higher resolution simulations showing that stellar feedback strongly suppresses black hole growth in lower mass galaxies (Habouzit et al. 2017;Anglés-Alcázar et al. 2017c).

While all the simulations presented above use a single model to compute the accretion on to the BHs, SIMBA uses a two mode model: torque-limited accretion model for cold gas (T < 10 5 K) and the Bondi-Hoyle-Lyttleton model for hot gas (T > 10 5 K). The accretion rate is written as

where r = 0.1, and ṀTorque is the gas inflow rate driven by gravitational instabilities from the scale of the galaxy to the accretion disc of the BH, i.e. here within R 0 = 2h -1 kpc (Hopkins & Quataert 2011;Anglés-Alcázar et al. 2015, 2017a). The gas inflow rate ṀTorque is defined as

with T a normalization parameter, f d the mass fraction (gas + stellar content) of the disc, f gas the gas fraction of the disc, M enc (R 0 ) the mass of the gas and stellar content. A normalization parameter α = 0.1 is included in the Bondi mode to represent the efficiency of gas transport from the accretion disc down to the black hole (for consistency with the gravitational torque rate).

The simulation employs a modified version of the kinetic AGN feedback model of Anglés-Alcázar et al. (2017a). In the high accretion rate mode (f Edd 0.2, radiative mode), high mass loading outflows are employed, and the velocity of the outflows increases with log 10 M BH . In this mode the outflow velocities are typically < 1000 km s -1 . In the low accretion rate mode of AGN feedback (f Edd < 0.2, jet mode), lower mass loading but faster outflows are used. In this mode the velocity of the jets increases for decreasing f Edd . Full speed jets with v jet ∼ 8000 km s -1 are only achieved at f Edd < 0.02. Only BHs of M BH 10 7.5 M are allowed to transition into this jet mode. Finally, X-ray feedback (following Choi et al. 2012) is included for galaxies with low gas fraction (f gas < 0.2) and velocity jets of v jet 7000 km s -1 (Davé et al. 2019). Finally, the feedback kinetic efficiencies = 0.03 (radiative feedback mode) and = 0.3 (jet), listed in Table 1, depend on the outflow velocity. The numbers quoted here are representative of BHs with M BH = 10 9 M .

Stellar feedback in SIMBA is modelled via hydrodynamically decoupled winds as in Mufasa (Davé et al. 2016). The mass loading factor depends on the stellar mass of the galaxy following scalings from the FIRE simulations (Anglés-Alcázar et al. 2017b). Similar to TNG, the mass loading factor is flat in galaxies with M 3 × 10 8 M and is reduced at high redshift to avoid excessive feedback in these galaxies, and allows them to grow when poorly resolved. The typical wind velocity are substantially lower than TNG100 (∼1.6 × v circ ), and 30 per cent of the winds are ejected hot. SIMBA also includes Type Ia SNe and AGB wind heating (Davé et al. 2019).

In the following, we generally refer to the mass of the BHs as the mass of the most massive BHs in the galaxies. For Illustris and TNG, BH mass refers to the sum of the BHs in the galaxies, but since BHs are automatically merged when their galaxies merge it does not make any difference (in practice only one or two galaxies per snapshots host several BHs).

In most of the simulations BHs (not in Horizon-AGN) are re-positioned every time-step to the local gravitational potential minimum of the galaxies or within smaller regions, which favours the accretion of gas on to the BHs. For example in EAGLE, BHs are moved to a particle with a lower potential but only within their kernel, and only if the relative velocity of the particle is less than 0.25 c s (Schaye et al. 2015), and if its distance is within three gravitational softening lengths. Finally, BHs merge if they are close enough to each other. Some simulations have additional criteria, such as Horizon-AGN, EAGLE, SIMBA, and TNG. For instance, in the EAGLE simulation BHs are merged if their separation distance is smaller than both the SPH smoothing kernel of the BH and three gravitational softening lengths and if the relative velocities of the BHs at separations of the SPH kernel are less than the circular velocity at this distance (Salcido et al. 2016). In TNG, BHs merge if within the accretion/feedback region of another BH.

Sub-grid models usually involve some efficiency parameters that can be tuned in order to produce a population of objects, here BHs, Table 1. Parameters and models of BH and galaxy formation/evolution in the simulation Illustris, TNG100, TNG300, Horizon-AGN, EAGLE, and SIMBA. We include the quantities related to the volume and resolution of the simulations, the seeding prescriptions (i.e. minimum halo mass seeded, minimum cell density, velocity dispersion), and the BH mass of the seeds, the parameters of the BH accretion models (models and boost factors), the SN feedback models (models, efficiencies, and energy released per core collapse SN), and finally the parameters related to AGN feedback (number of modes, models, efficiencies, transition between modes). Seeding prescriptions Häring & Rix (2004) as a solid yellow line. We also show the two M BH -M relations of Reines & Volonteri (2015) in dashed orange lines: the line on the top for elliptical and spiral/S0 galaxies with classical bulges, and the bottom line for the broad-line AGN. The full sample of Reines & Volonteri (2015) is shown as well with coloured dots. We also show the sample of spiral galaxies of Davis et al. (2018) with dynamical BH mass measurements. Finally, we also reproduce here the sample of AGN of Baron & Ménard (2019), which includes obscured AGN, and their scaling M BH -M scaling relation as a dotted pink line. We draw here a cartoon black region including a large fraction of the observations that we use in the following.

matching a given observational constraint. These models include the modelling of the accretion on to the BHs and of AGN feedback. The parameters related to seeding, e.g. the initial mass of the BH particles, are also often chosen to produce a BH population matching one of the z = 0 empirical scaling relations. We show some of the empirical relations commonly used in the literature in Fig. 1, i.e. the relation M BH -M bulge of Häring & Rix (2004), Kormendy & Ho (2013), McConnell & Ma (2013), and the relation M BH -M of Reines & Volonteri (2015), Baron & Ménard (2019).

All the simulations were calibrated to one of the empirical M BH -M bulge relations, with some variations in the computation of M bulge . The Illustris simulation was calibrated to the M BH -M bulge relation of Kormendy & Ho (2013), assuming that the total stellar mass within the stellar half-mass radius of the simulated galaxies was a proxy for M bulge . Horizon-AGN was calibrated on the local scaling M BH -M bulge relation of Häring & Rix (2004) (Dubois et al. 2012) (using M bulge for the simulated galaxies as well), to determine what fraction of rest-mass accreted energy should be released as AGN feedback. The simulation EAGLE was calibrated to the McConnell & Ma (2013) relation between the mass of BHs and the bulge mass of their host galaxies, assuming that the stellar mass of the observed galaxies was dominated by their bulge. The simulation was therefore calibrated assuming that M bulge = M , which may be a correct for the highest mass galaxies. The SIMBA simulation was calibrated on the amplitude of the M BH -M bulge relation of Kormendy & Ho (2013) (assuming that M bulge = M ) at log 10 M BH /M ∼ 8-9 through the efficiency parameter of the accretion model/AGN model (see Anglés-Alcázar et al. 2017a, for more details). All simulations based on Bondi accretion (Illustris, TNG, Horizon-AGN, EAGLE) match the normalization of the M BH -M bulge relation by tuning a feedback efficiency parameter, while SIMBA (based on gravitational torque accretion) calibrates the accretion efficiency to match the M BH -M bulge relation. This difference is primarily driven by the BH mass dependence of the accretion parametrization (Anglés-Alcázar et al. 2015, 2017a).

More generally, the simulations are calibrated with observational constraints of galaxies, either by directly using the fits of these constraints (EAGLE), or by simple agreement with the constraints (Illustris, IllustrisTNG). IllustrisTNG was not directly calibrated by eye with any new observational constraints, but instead the sub-grid models of the previous Illustris were adapted to better match the Illustris galaxy properties to observations used for the calibration of the Illustris simulation. The EAGLE simulation was calibrated to fits of observational data representing the galaxy mass function at z = 0, galaxy sizes at a function of galaxy mass at z = 0.1, and the local M BH -M relation (Crain et al. 2015;Schaye et al. 2015). The Illustris simulation was calibrated based on cosmic star formation rate density, galaxy mass function at z = 0, stellar to halo mass ratios at z = 0, as well the gas metallicity mass relation at z = 0, and the local M BH -M relation (Vogelsberger et al. 2013;Torrey et al. 2014). Similarly, IllustrisTNG was calibrated based on the cosmic star formation rate density, the galaxy mass function, the stellar-tohalo ratios, and the M BH -M as well, by comparison to the results of the Illustris simulation. Calibrations to the gas fraction at z = 0, and to the galaxy sizes as a function of galaxy masses at z = 0 were added. The SIMBA simulation was calibrated on the galaxy stellar mass function and its evolution with time. No calibration on galaxy sizes or galaxy/halo gas content was made. Finally, the Horizon-AGN simulation was simply calibrated with the M BH -M relation to finalize the AGN feedback model, and the rest of the subgrid physics (star formation, SN feedback) was just the result of the underlying model (no calibration on the galaxy stellar mass function for instance). Reines & Volonteri (2015) for the local Universe (z ∼ 0) in black dots in the last panel on the right (uncertainties on BH masses are ∼ 0.5 and ∼ 0.3 dex for the stellar masses). To guide the eye, we define and duplicate in each panel a cartoon region of the diagram including most of the observations. Two main discrepancies emerge from this figure: some simulations do not produce the most massive BHs observed in galaxies with M ∼ 10 11 M and the broad-line AGN of M BH ∼ 10 6-7 M observed in galaxies of M ∼ 10 10.5-11 M .

We compute in post-processing the luminosity of the BHs with the model of Churazov et al. (2005), i.e. explicitly distinguishing radiatively efficient and radiatively inefficient AGN. The bolometric luminosity of radiatively efficient BHs, i.e. with an Eddington ratio of f Edd > 0.1, is defined as

BHs with small Eddington ratio of f Edd 0.1 are considered radiatively inefficient and their bolometric luminosities are computed as

The distinction between radiatively efficient and inefficient AGN is often not made when computing the luminosity of the AGN in simulations. Therefore, our post-processing computation of the luminosities differs from the intrinsic luminosity that goes into AGN feedback in the simulations. We use the radiative efficiency parameter which were used to derive the accretion rate self-consistently in the simulations, i.e. r = 0.2 for Illustris, TNG100, and TNG300, and r = 0.1 for Horizon-AGN, EAGLE, and SIMBA.

We focus our investigation on the scaling relation between BH mass and the total stellar mass of their host galaxies, a quantity that can be measured beyond the local Universe. Empirical scaling relations have been derived with other galaxy properties (e.g. stellar velocity dispersion, bulge luminosity, bulge mass), and we discuss these quantities in Appendix A. Given the differences found in the M BHσ relation of Illustris and TNG100 and observations (Li et al. 2020), we prefer to not investigate this relation here. In this section, we present several versions of the M BH -M diagrams of the simulations. While we do not intend to broadly discuss accretion properties of the BHs and the correlation with star-forming activity of their host galaxies, we present here some insights into these correlations.

Before investigating the M BH -M diagram in cosmological simulations, we show in Fig. 1 several observational samples. The sample of Reines & Volonteri (2015) includes a small sample of dwarf galaxies (shown in pink), 262 broad-line AGN (light blue), reverberationmapped AGN (in purple), and 79 galaxies with dynamical BH mass measurements (dark blue). The BHs present in RGG118 (Baldassare et al. 2015) and Pox 52 (Barth et al. 2004) are shown in green.

We add the sample of Davis et al. (2018) with dynamical mass measurements in dark blue. Finally, we show the recent observational sample of Baron & Ménard (2019), which include obscured AGN. In the following, we compare the simulations to the sample of Reines & Volonteri (2015), which covers the same M BH -M regions as most other observations. In Fig. 1 we add several empirical scaling relations that are commonly used in the literature, both M BH -M bulge (Häring & Rix 2004;Kormendy & Ho 2013;McConnell & Ma 2013) and M BH -M relations (Reines & Volonteri 2015;Baron & Ménard 2019). The empirical scaling relations and even the samples themselves can be biased, in many different ways. That is why we do not aim at an apple-to-apple comparison with observations, which would require us to take into account all the selection effects and biases (e.g. galaxy selection, BH mass measurements). Moreover, most of the samples are biased towards luminous unobscured AGN (but see Baron & Ménard 2019), or massive BHs whose sphere of influence can be resolved for mass dynamical measurements. For example, selection bias in dynamically measured BH sample, whose host galaxies often appear to have higher velocity dispersion σ , could artificially enhance BH masses (Bernardi et al. 2007). Since the empirical M BHσ relations are employed to estimate BH masses of AGN, this could impact more broadly the observational samples by a factor of a few (Shankar et al. 2016)

We present the M BH -M diagrams of the simulations at z = 0 in Fig. 2 in hexabins colour-coded by the number of BHs in each bin. For comparison, we show the observational sample of Reines & Volonteri (2015) in the right-hand panel of the figure. Their AGN sample from SDSS goes out to z 0.055, corresponding to a distance of ∼ 230 Mpc (similar to the simulations length) but does not cover the entire sky. To guide the eye, we report these observations with a black region in all the panels. This black region is not rigorously defined, but is drawn by eye to include all the observations except a few isolated data points. The region purposely looks like a cartoon to show that we do not make any analysis of the observations here.

The simulated populations of BHs/galaxies form a tight correlation between BH mass and the stellar mass of the host galaxies in the diagrams, and in most cases, the scatter of the observed BH population is not fully reproduced in the simulations. Here, we list some of the differences between simulated and observed samples of BHs.

(i) The broad-line AGN located in the region enclosed within M BH = 10 6-8 M and M = 10 10.5-11 M are not very well reproduced by the simulations, in particular by the TNG100 and the Horizon-AGN simulations.

(ii) The BHs of M BH = 10 7 M in galaxies of M 10 11 M are not produced by any of the simulations, with the exception of SIMBA and TNG300.

(iii) Similarly, some of the simulations, such as TNG100, Horizon-AGN, and especially EAGLE, also have a hard time producing the most massive BHs that we observe in galaxies with total stellar mass of a few 10 10 M and 10 11 M . These BHs are the most massive BHs at fixed stellar mass, and their number in observations is low. Therefore, their assembly in the different simulations could be limited by the small simulated volumes (∼ (100 cMpc) 3 ).

(iv) All the simulations predict BHs of M BH ∼ 10 10 M in galaxies of M ∼ 10 12 M at z = 0 (coloured dots located beyond the left-hand side of the black line region, which is constrained by this observational sample). The mass of BHs in massive galaxies is often determined by dynamical measurements, which is only possible for close systems. In observations, these rare massive galaxies are often found on the centre of clusters.

(v) All the simulations (with the exception of TNG) predict BHs of M BH 10 6 M in galaxies of M = 10 10-11 M at z = 0, although with different number densities. SIMBA produces quite a lot of these BHs compared to the other simulations. This region of the diagram is only slightly populated in the observed plane. In the observation panel on the left of Fig. 2, these BHs appear below the black line region. This BH mass range is not probed by the TNG model that employs a higher seeding mass.

(vi) Only very few BHs have been observed with M BH > 10 8 M in galaxies of M < 10 10 M ; only one BH is present in the sample of Reines & Volonteri (2015) (black dot above the black shape in this galaxy mass range). These BHs are rare in simulations. TNG300 is the simulation forming the highest number of these BHs, due to its larger volume. Additionally, we also note that the observations shown here only include a small number of BHs with M BH > 10 6.5 M in galaxies of M < 10 10 M (BHs above the black lines), while most of the simulations predict a non-negligible number of them. EAGLE is the only simulation in good agreement in this stellar and BH mass range.

Here, we compared the simulated BH populations with the observational sample of Reines & Volonteri (2015), which was made to only include dynamical measurements or estimates based on single-epoch virial masses for broad-line AGN. As shown in Fig. 1, this sample includes massive BHs in massive inactive galaxies, which are elliptical and spiral/S0 galaxies with classical bulges, and lower-mass BHs found in spiral and pseudo-bulge galaxies. Recently, Baron & Ménard (2019) has discussed the possibility of using only narrow emission lines to estimate BH masses (while broad lines are generally used), therefore allowing BH mass determination for obscured AGN (type II) in addition to the non-obscured AGN (type I) used in Reines & Volonteri (2015). The sample of Baron & Ménard (2019) covered the region between the lower mass broadline AGN of Reines & Volonteri (2015) and its dynamical mass measurement BHs, i.e. the region defined by M BH = 10 7 -10 9 M BHs in M 10 10 M galaxies. The sample of Baron & Ménard (2019) confirms that the simulations are generally missing some of the most massive BHs at fixed stellar mass below M = 10 11 M , and confirms that compared to the observations the simulations do not produce enough BHs in the bottom right-hand side of the M BH -M diagrams, i.e. BHs of M BH = 10 7 -10 8.5 M in galaxies of M 10 11 M . Reines & Volonteri (2015) establish two distinct scaling relations for massive BHs in quiescent elliptical galaxies and lower mass BHs observed as AGN, showing that these populations may be distinct populations of BHs. They also discuss the possibility that they could also be sub-populations of the global BH-galaxy population, and that e.g. quiescent low-mass BHs could overlap with the AGN but would simply not be detectable. The work of Baron & Ménard (2019) reinforces this idea, and shows that AGN (obscured AGN) in massive galaxies can overlap with quiescent BHs in quiescent elliptical galaxies. Whether these populations are part of the same population will become clearer in time as we increase our ability to accurately measure BH masses in the local Universe.

In Fig. 3, we show the M BH -M diagrams colour-coded by BH bolometric luminosity. While we do not intend to analyse in detail the accretion properties of the AGN in this paper (this is the focus of the second paper of our series), the accretion rates, and luminosities are important to understand the different processes involved in the evolution of the BH populations. Below, we describe some specific aspects of the BH sub-grid physics that can be identified in Fig. 3.

The high seeding mass of the TNG simulations can be easily seen in Fig. 3: no BHs form with masses lower than 10 6 M . As indicated in the last bottom left-hand panel of Fig. 2, lower-mass BHs have been observed in the local Universe. These BHs of M BH ∼ 10 4-5 M are detected as AGN in dwarf galaxies of M ∼ 10 9 M (Reines, Greene & Geha 2013;Baldassare et al. 2015;Reines & Volonteri 2015;Mezcua et al. 2016); probably only the most massive/luminous BHs are detected. The presence of BHs in dwarf galaxies is also not covered in the simulation SIMBA, which starts forming BHs in galaxies with M 10 9.5 M . While modelling BH formation in low-mass galaxies is crucial to understand BH formation in the highredshift Universe, as well as to understand the current populations of BHs in local dwarf galaxies, the regime of low-mass galaxies is barely resolved in such large-scale simulations of 100 cMpc on a side (but see Habouzit et al. 2017, and references therein for BH formation in low-mass galaxies from high redshift to low redshift). All the other simulations employ lower seeding masses than the TNG simulations. In Horizon-AGN fewer BHs of M BH 10 6 M are present at z = 0 than at higher redshift. This is because BH formation stops at z = 1.5 by design in this simulation.

In simulations, the initial mass of BHs was shown to affect the lowmass end of the M BH -M diagram, and the overall normalization of the scaling relation (e.g. Bower et al. 2017). For higher stellar masses, the relations for different seeding masses converge to the same relation (whose normalization depends on the BH accretion efficiency parameter) as a result of self-regulation. The seed mass is important for the simulations using the Bondi accretion model (e.g. Illustris, TNG, Horizon-AGN, EAGLE), for which there is a degeneracy between the seeding mass and the boost factor to produce the same normalization of the M BH -M mean relation. In simulations using a torque accretion model (SIMBA), the accretion rate on to the BHs does not strongly depend on their masses ( ṀBH ∝ M et al. 2015). We point here again that the seeding in large-scale cosmological simulations is often not physically motivated (e.g. fixed BH mass, threshold in galaxy, or halo mass to form a BH), but instead parameters of the seeding models are chosen to reproduce the local scaling relation.

Getting the BHs on the main scaling relation depends both on the seeding BH mass and BH accretion rate. The parametrizations of these models are degenerate. Therefore, different choices allow the BH population to get on to the same empirical scaling relation. For example, the TNG model seeds BHs with a higher initial BH mass than the Illustris model (by one order of magnitude) but also does not include any boost factor in the accretion model (while Illustris does include a boost factor of 100).

The colour code of Fig. 3 provides information on the accretion rates on to the BHs. At z ∼ 2, it is clear that the accretion properties are very different among the different simulations. For example, BHs of M BH ∼ 10 6 M embedded in galaxies of M ∼ 10 10 M have high bolometric luminosities of L bol 10 44 erg s -1 (limit commonly employed to define an AGN) in TNG, Horizon-AGN, and SIMBA, while a large fraction of these BHs in Illustris and EAGLE have much lower luminosities. Interestingly, we see that the initial growth of BHs in SIMBA is quite efficient, and it allows them to get rapidly on the track of the main scaling relation even if only galaxies of M 10 9.5 M get seeded. The accretion model of SIMBA based on gravitational torques is more efficient in the regime of these low-mass BHs than the Bondi accretion (which scales as M 2 BH ) (Anglés- Alcázar et al. 2015Alcázar et al. , 2017a)). In SIMBA, the torque model and the Bondi model always co-exist, but in practice the torque model dominates at early times in gas-rich galaxies: cold gas dominates the gas reservoir in the centre of galaxies, and there is only little hot gas to be accreted through the Bondi model (Angles-Alcazar, in preparation). We colour-code the simulated BHs with the specific star formation rate (sSFR) of their host galaxies. We set sSFR = 10 -13 yr -1 for galaxies with lower sSFR. The limit sSFR = 10 -11 yr -1 is often use to define star-forming galaxies (with higher sSFR) or quiescent galaxies (with lower sSFR). In the figure, galaxies with blue colours are forming stars, and yellow to red galaxies are quiescent. We show in the last right-hand panel the observed star-forming and quiescent galaxies of Terrazas et al. (2017) for references, colour-coded by their sSFR as well. To guide the eye, we again reproduce the black region representing the observational sample of Reines & Volonteri (2015) in each panel. A similar figure (Fig. B1) shows the behaviour at higher redshifts.

In Fig. 3, we can identify a clear cut at fixed BH mass in the bolometric luminosity of BHs in the TNG100 and TNG300 simulations (at M BH ∼ 10 8 M ). Above this characteristic BH mass the accretion rate on to the BHs and consequently their luminosity is strongly reduced. In TNG, most of the M BH ∼ 10 8 M BHs with low accretion rates transition from the thermal high accretion feedback mode to the kinetic wind low accretion feedback mode, which is the mode responsible for efficiently quenching massive galaxies (Weinberger et al. 2018;Habouzit et al. 2019).

In EAGLE, we note the presence of faint AGN or inactive BHs with L bol 10 38 erg s -1 and M BH 10 7 M in galaxies of M ∼ 10 10 M (mostly visible at z = 4, 2). This is due to AGN feedback. BHs in the EAGLE simulation are first regulated by the efficient SN feedback in low-mass galaxies of M 10 10 M , and then have a period of efficient growth, before being regulated by AGN feedback in galaxies of M 10 10 M (Bower et al. 2017;McAlpine et al. 2017McAlpine et al. , 2018)).

In the other simulations, the different AGN feedback models do not produce any strong signature in the M BH -M diagrams. We note a small signature in SIMBA at z = 0, i.e. a weak horizontal line of demarcation at M BH 10 7.5 M , where BHs seem to have lower luminosities on average. In SIMBA, only BHs of M BH 10 7.5 M and with accretion rates such that f Edd 0.2 are allowed to enter the efficient kinetic mode of AGN feedback. Once they enter this mode, BHs regulate themselves, as well as the star formation in their host galaxies (Thomas et al. 2019) similar to TNG.

In Fig. 4, we reproduce again the M BH -M diagram at z = 0 for the different simulations, but this time we colour-code the BHs by the specific star formation rate (sSFR) of their host galaxies. In the literature, we often use the convenient limit sSFR = 10 -11 yr -1 to separate star-forming galaxies from galaxies with low star formation rate or quenched galaxies (i.e. with low SFR sustained in time). In Fig. 4, star-forming galaxies appear in bluish colours and quiescent galaxies in reddish ones. On the last panel, we show the sample of Terrazas et al. (2017) foot_4 that only include BHs with dynamical mass measurements, and galaxies with SFR estimated from far-infrared fluxes (from IRAS).

We find a good qualitative agreement between the observations and the simulations. In Terrazas et al. (2017), observed galaxies of M = 10 10 -10 12 M hosting the most massive BHs (M BH 10 8.5 M ) have lower sSFR (sSFR 10 -11 yr -1 ). Similarly, in Reines & Volonteri (2015) these galaxies are ellipticals, and galaxies with classical bulges (mostly early-type galaxies). In the simulations, these galaxies also tend to have low specific star formation rates, and to be quenched. We note that some of the simulations have a large fraction of quenched galaxies with sSFR 10 -13 yr -1 (Illustris, TNG, SIMBA), while the other simulations (Horizon-AGN, EAGLE) have quenched galaxies with a broader range of sSFR (10 -13 sSFR 10 -11 yr -1 ). While all simulations rely on AGN feedback to solve the overcooling process in massive galaxies and produce a galaxy stellar mass function in good agreement with observations, they all use different feedback modellings. Consequently and as seen here, these modellings produce different behaviours and level of quenching in the simulations.

In the sample of Terrazas et al. (2017), most of the galaxies with lower-mass BHs (M BH 10 8 M ) are forming stars more efficiently and have higher sSFR of sSFR > 10 -11 yr -1 . This is true for galaxies with M 10 10.5 M , and we note that, in this sample, galaxies with stellar masses of M ∼ 10 10 M have lower sSFR. In the sample of Reines & Volonteri (2015), these galaxies are also forming stars as they are mostly spiral galaxies and pseudo-bulges galaxies (found in late-type galaxies).

As we discussed earlier in this paper, the region of the M BH -M diagram corresponding to the low-mass BHs in the observational samples of Reines & Volonteri (2015) and Terrazas et al. (2017) does not overlap very well with the low-mass BHs formed in the simulations. In the observations BHs of M BH = 10 6-8 M are found in galaxies of M = 10 10.5-11.5 M , while in the simulations these BHs are generally found in less massive galaxies. If we compare the sSFR of these observed galaxies with BHs of M BH = 10 6-8 M with the sSFR of lower stellar mass galaxies with the same-mass BHs, we find a good agreement. Moreover, we note that some of the samemass BHs in even lower mass galaxies also have lower sSFR in the simulations, as found in the observations of Terrazas et al. (2017).

To summarize, we do find a good qualitative agreement between the sSFR properties of the simulated and observed galaxies hosting the most massive BHs, and we find the same trend for lower mass BH host galaxies (higher sSFR for BHs in more massive galaxies) but the observed and simulated galaxies populate different regions of the M BH -M diagram.

We do not aim at comparing in detail the samples of Terrazas et al. (2017) and Reines & Volonteri (2015), but we note some correlations in the following. The quenched galaxies of Terrazas et al. (2017) overlap with the inactive galaxies of Reines & Volonteri (2015) (elliptical and spiral/S0 galaxies), whose BHs are also inactive and have dynamical mass measurement. BHs with lower masses in Reines & Volonteri (2015) are mostly AGN, whose masses are estimated from broad-line emission, and BHs found in spiral and pseudo-bulge galaxies; they cover the same region as the starforming galaxies of Terrazas et al. (2017). From the comparison of these two samples, we note a correlation between quenched galaxies and inactive massive BHs, and similarly between lower mass BHs detected as AGN and star-forming galaxies. Such correlations need to be taken with a grain of salt, as distinct regions of the M BH -M diagram may be populated by different selections and methods of mass measurements/estimates. Mass estimates for the AGN can be based on empirical relations derived from specific galaxy populations (e.g. star-forming galaxies). While BHs with mass dynamical measurements can also be biased to galaxies with e.g. higher velocity dispersion (Bernardi et al. 2007;Shankar et al. 2016, and references therein).

We now describe what we find in the simulations at z = 0 (in the last row of Fig. 3). In the TNG and EAGLE simulations, the massive BHs in massive galaxies (M BH 10 8.5 M , M 10 10.5 ) tend to have low accretion rates. In Illustris, Horizon-AGN, we find both non-active and active BHs among the massive BHs. There is also a significant high fraction of massive BHs with high accretion rates in SIMBA. Comparing Figs 3 and4, we find a clear correlation between quenched galaxies and non-active massive BHs in Illustris, TNG, as in some of the observations mentioned above. The correlation is less obvious in Horizon-AGN and EAGLE. The massive galaxies at z = 0 in SIMBA can both feed the massive BHs and have low sSFR. The correlation stands for lower mass BHs in the TNG simulations, in Horizon-AGN, in SIMBA, i.e. the host galaxies of lower mass BHs form stars more efficiently. The picture seems harder to establish for the low-mass galaxies in the Illustris and EAGLE simulations. Here we only look at the instantaneous accretion rates on to the BHs and luminosities, but the AGN luminosities can vary substantially within short amounts of time (e.g. fig. 1 of Rosas-Guevara et al. 2019). Correlations between star formation rates and average AGN luminosities over a few hundreds Myr could be stronger. In a forthcoming paper of our series, we will study the correlation between BH activity and star formation of the host galaxies in a more quantitative way. In the following section, we discuss in more details the time evolution, normalization, scatter, and shape of the M BH -M relation in the simulations.

In Fig. 4, all the simulations seem to favour a linear M BH -M relation for galaxies of M 10 11 M , a mass range dominated by quiescent galaxies (in red). However, for lower mass galaxies which are dominated by star-forming galaxies, a linear relation is not supported by all the simulations. Observationally, this correlation between the quiescence of the host galaxies and the shape of the M BH -M relation has been discussed (e.g. Graham & Scott 2015).

Different slopes of the scaling relation have been found for samples of early-type/core-Sersic and late-type/Sersic galaxies in observations (Davis et al. 2018;Sahu et al. 2019). In the next section, we investigate the time evolution, scatter, shape, and normalization of the M BH -M relation.

We now turn to study the evolution of the M BH -M relation and its scatter with time. Our goal is to identify the different features seen in the simulations; the physical interpretations of these features is given in the next section. When possible, we show observational constraints to guide the eye, but we do not aim at building apple-toapple comparisons between the observations and simulations.

We show in the top panels of Fig. 5 the median of the M BH -M relation for several redshifts. We obtain almost identical relations for the mean relations. For reference, we show a grey shaded area in each panel enclosing the local M BH -M bulge empirical relations of Kormendy & Ho (2013), McConnell & Ma (2013), Häring & Rix (2004), which were used to calibrate the simulations (with different assumptions for M bulge ). In Fig. 5 we show the total stellar mass of the simulated galaxies M , which should be comparable to the bulge mass for the highest mass galaxies. In the bottom panels, we show the 15th-85th percentiles of the M BH distributions in stellar mass bins. Please note that in this section, we use the same terminology scaling relation for the mean/median of the M BH -M population in the simulations and for the empirical scaling relations (which are not simply defined by the mean/median of the observational samples).

We find that the scaling relations of all the simulations evolve with time. However, the evolution of the overall normalization is small, particularly in the redshift range that can currently be probed by observations 0 z 2 (red, yellow, and brown lines in Fig. 5). In this redshift range, the time evolution is smaller than one order of magnitude in BH mass. The strongest evolution is found for the Illustris simulation for these redshifts and for M 10 11 M . When looking at the time evolution of the overall normalization, two trends emerge. Half of the simulations have a lower overall median M BH -M relation with decreasing redshift; this is the case for the Illustris, Horizon-AGN, and EAGLE simulations. In the EAGLE simulation, the relation only evolves with time for intermediate-mass galaxies. The other TNG and SIMBA simulations have scaling relations with higher overall normalization.

The time evolution of the scaling relation is likely the result of several physical processes affecting the growth of both the BHs and their host galaxies. Among these processes, the ability of BHs to accrete gas, the feedback of the AGN and of SNe, the growth of the host galaxies, can play a crucial role. There are different hypotheses to explain the two trends in the M BH -M diagram. A higher overall normalization at higher redshift (Illustris, Horizon-AGN, EAGLE) can be the signature of a more rapid growth of the BHs at higher redshifts. Theoretically, we expect galaxies to be gasrich and more compact at higher redshift, which favours efficient BH accretion. While at lower redshift (e.g. z 1), BH growth could be less efficient. The second hypothesis is that BH growth has a constant efficiency with time, but that galaxies have a relative growth faster at low redshift than at high redshift compared to their BHs. It would lead to a shift of the M BH -M relation towards more massive galaxies at lower redshifts. An increase of the overall scaling 2004) (lower normalization of the relation). Bottom panels: 15-85th percentile of the distributions. While in general the median of the M BH -M relation does not evolve much with time, we do note some differences between the simulations. These differences can also be appreciated in the previous figure. Some simulations increase their scatter toward lower mass BHs with time (at fixed stellar mass), i.e. the median relation moves down in the M BH -M diagram for lower redshift. This is true for Illustris, Horizon-AGN, and EAGLE. However, the TNG and the SIMBA simulations evolve upward in the diagram.

relations with time (TNG, SIMBA) could be explained by a relative more efficient growth of BHs at low redshifts with respect to their galaxies, or by a smaller growth of the galaxies with respect to their BHs. We investigate how these hypotheses could drive the two evolution patterns found here in more detail in the following section of the paper.

For reference, most of the observational studies have concluded that the time evolution of the scaling relation between z ∼ 2 and z ∼ 0 was consistent with no evolution (Shields et al. 2003;Jahnke et al. 2009;Cisternas et al. 2011;Salviander & Shields 2013;Schramm & Silverman 2013;Sun et al. 2015;Suh et al. 2020). A few studies have shown that the ratio M BH /M could be larger at higher redshifts (e.g. McLure et al. 2006;Ding et al. 2020). For example, Ding et al. (2020) find that the small evolution in their sample (32 X-ray selected broad-line AGNs, 1.2 < z < 1.7) goes into this direction, i.e. higher BH masses or lower stellar masses at higher redshift.

Finally, studying the time evolution of the scaling relation for the most massive galaxies of M 10 11 M is problematic with the current large-scale simulations of side length ∼ 100 cMpc, as they suffer from low number statistics. The large volume of TNG300 captures the evolution of a higher number of massive galaxies, especially at higher redshift. Indeed, there are only 270 galaxies in TNG100 and 4378 in TNG300 with 10 log 10 M /M 11 at z = 4. Only two galaxies with log 10 M /M 11 in TNG100, and 66 in TNG300. We find that there is a rapid evolution of the massive end (log 10 M /M 10.5) of the relation for 3 z 5 in TNG300, a regime that is not accurately covered by the other simulations. 5 At z 3, we note the absence of an evolution of the relation with time. The BHs embedded in these massive galaxies are regulated by the efficient kinetic feedback. Therefore, their growth is likely driven by mergers only (and not by gas accretion); the growth 5 SIMBA (147 cMpc box length) also produces a high number of massive galaxies, e.g. 30 galaxies of log 10 M /M 11 at z = 4.

of their host galaxies is likely also mostly driven by mergers. From the central-limit theorem we expect the M BH -M relation to be mostly independent of redshift for the most massive galaxies, and the distributions to have a smaller scatter. This is true for the TNG simulations (Weinberger et al. 2017(Weinberger et al. , 2018)), and can also be seen in the other simulations with an efficient AGN feedback even if lacking statistics for massive galaxies. In SIMBA, we find an evolution of the relation for these massive galaxies of M 10 11 M : BH growth (through the Bondi accretion channel) exceeds the relative galaxy stellar mass growth by mergers (Cui et al., in preparation).

We now turn to study the evolution of the scatter of the M BH -M relation with both with stellar mass and redshift. We define the scatter as the 15th-85th percentiles of the BH mass distributions in stellar mass bins. We provide the stellar mass bins (bin size is 0.4 dex), redshifts, mean and median BH mass in the bins, and the 15th-85th percentiles in Tables C1 andC2. We split the tables in two: simulations with an increasing M BH -M overall normalization with decreasing redshift (Illustris, Horizon-AGN, EAGLE), and those with a decreasing overall normalization (TNG100, SIMBA). We find some differences between all of these simulations. For most of the simulations, the scatter of the scaling relation is below 1 dex in log 10 M BH /M , except for EAGLE which has the largest scatter (e.g. > 1 dex for M = 10 10 , 10 10.5 , 10 11 M for z 3). Horizon-AGN has the smallest scatter, below half a dex in log 10 M BH /M (for all redshifts and M bins).

For galaxies of M 10 11 M , the scatter in Illustris and TNG100 generally increases with decreasing redshift to z = 0. The scatter decreases in EAGLE with decreasing redshift, slightly oscillates in SIMBA for different redshifts, and does not evolve with redshift in Horizon-AGN. For all the simulations, the evolution of the scatter with redshift is smaller than 1 dex in BH mass. We define the time variation of the scatter by the difference of the 15th-85th percentiles between two given redshifts, i.e. σ = |σ iσ j |. Horizon-AGN has the smallest scatter evolution with σ 0.1 dex, while EAGLE has a variation of almost one order of magnitude in BH mass (i.e. σ = 0.84 dex), the largest among all the simulations. According to the central-limit theorem (i.e. assuming that growth is driven only by mergers) any distribution with a large scatter at high redshift would have a smaller scatter with time. The mild evolution of the scatter that we find here for M 10 11 M suggests that BH growth by gas accretion, in addition to mergers, play a role in the build-up of the scaling relation.

Comparing the time evolution of the scatter in simulations with observations is challenging. The scatter in observations is often defined by the 1σ confidence interval of the empirical scaling relations (about half a dex in log 10 M BH /M ), which relies more on the uncertainties of the mass estimates 6 rather than the intrinsic scatter of the BH mass distribution in stellar mass bins, as derived here for the simulations. Given that it is harder to estimate BH masses accurately at higher redshift, the scatter in observations is expected to increase with increasing redshifts. Recently, the intrinsic scatter of the M BH -M relation was found to be similar at z ∼ 1.5 and z = 0 (e.g. Ding et al. 2020).

Variations of the scatter as a function of the stellar mass (for M 10 11 M ) are small in the simulations, on average. The amplitude of those would probably be lower than the various uncertainties in estimating BH masses and host stellar masses in observations.

Interestingly, the simulations produce different evolutions of the scatter with M . At fixed redshift, the scatter decreases with larger M in SIMBA, and slightly increases for Illustris and Horizon-AGN. There is a large variation of the scatter in TNG100, TNG300, and EAGLE: the scatter is smaller in low-mass galaxies of M 10 9.5-10 M (depending on the redshift), is larger for galaxies of M ∼ 10 10 M , and decreases for more massive galaxies. The dependence with stellar masses in these simulations correlates with the efficiency of BH growth (see the next section with the shape of the M BH -M relation). In the low-mass galaxies, BH growth is not efficient and consequently a small scatter is found. However, when BHs start growing efficiently in galaxies of M ∼ 10 10 M the scatter is more pronounced (e.g. see McAlpine et al. 2018, for EAGLE). Finally, when the growth of the BHs is regulated by AGN feedback in more massive galaxies the scatter decreases. With the larger volume of TNG300 we find that the scatter for even more massive galaxies of M 10 11 M in the redshift range 0 z 3 (for which we have more statistics) is even smaller.

For reference, we compute the scatter (15th-85th percentiles) of the two observational samples of Reines & Volonteri (2015), Baron & Ménard (2019) for M 10 11 M and z ∼ 0. Our simple method does not reflect the scatter found in observations, and we do not add any correction (e.g. completeness of the sample, volume, low statistics). 6 Uncertainties on BH mass estimates are of about ∼ 0.5 dex (Vestergaard & Peterson 2006;Reines & Volonteri 2015). The uncertainties depend on the method used to measure/estimate BH masses. Uncertainties are smaller for masers, reverberation mapping, dynamical measurements, and larger for single epoch measurements.

The scatter found in observations also probably does not represent accurately the intrinsic scatter of the entire BH population in the Universe. This is even more true at the BH low-mass end due to very low number of detections. We discuss this further in the discussion section. We find that the percentile 15th-85th varies in the range σ = 0.9-1.9 in the stellar mass range M = 10 9.5 -10 11 M for the sample of Reines & Volonteri (2015), and within σ = 1.2, 1.6 for Baron & Ménard (2019). The scatter found in the simulations is smaller than for these two observational samples with on average more than a dex in log 10 M BH of scatter. The scatter increases with increasing M for Reines & Volonteri (2015), from σ = 1.1 (M = 10 9.5 M ) to σ = 1.9 (M = 10 11 M ), but the scatter decreases from σ = 1.6 (M = 10 9.5 M ) to σ = 1.2 (M = 10 11 M ) for Baron & Ménard (2019).

Some essential physical processes are not consistently modelled in large-scale simulations but could modulate the accretion, growth, and feedback of BHs, possibly impacting the scaling relation and its scatter. Here we only discuss the impact of BH spin, and mention more processes in the discussion. BH spin is closely tied to both accretion and the energy that can be released by AGN feedback. BHs with M BH < 10 8 M build their mass by coherent accretion of gas and a few mergers: they have high spins (Dubois, Volonteri & Silk 2014b). More massive BHs of M BH > 10 8 M experience more mergers and less coherent accretion: they have more moderate spins (see also Izquierdo-Villalba et al. 2020). BHs with high spins will release more specific energy than non-rotating BHs, which alters as well the amount of gas that is accreted by the BHs, and consequently, their spins. Adding the spin evolution in simulations increases the scatter, especially in the massive end of the scaling relation, where the relation is very tight. In the TNG model, the scatter is mostly increased for stellar masses of M ∼ 10 10 M and higher (Bustamante & Springel 2019).

The M BH -M relations presented in Fig. 5 have different shapes, and we investigate qualitatively the presence of a change of slope in the relations for the six simulations.

The median M BH -M relation of Illustris is linear at all redshifts and for all stellar masses. The linear relationship between BH mass and galaxy total stellar mass indicates that, on average, the BHs and their host galaxies grow at similar rates. When looking at the population of BHs statistically, there is no galaxy mass regime in Illustris at which BH growth is strongly regulated with respect to the growth of the galaxies and vice-versa. We also cannot identify any change of slope of the scaling relation of Horizon-AGN for z 2. However, we do see a transition around M ∼ 10 10 M at lower redshift. This is a signature of BH seeding: the formation of BHs stops at z = 1.5 in Horizon-AGN, and there are no more newly formed BHs to bring down the scaling relation at lower redshift.

All the other simulations have a linear evolution in massive galaxies (M 10 10.5 M ), but present a change of slope in the regime of the lower mass galaxies (M 10 10.5 M ). This change of slope reflects a transition between two phases: a first phase in which BHs have a hard time growing in low-mass galaxies, and a second phase for which BH growth is more efficient in more massive galaxies. At the transition between these two regimes, we can define a characteristic galaxy stellar mass, from which BHs can start growing efficiently. In TNG100, we note a change of slope at z 2 only, while in TNG300 the change of slope is present at all redshifts. The transition between inefficient BH growth and efficient BH growth changes with redshift in the range M ∼ 10 9.5 -10 10 M : the transition takes place in lower mass galaxies with time. The inefficient BH growth phase is longer at earlier times than at later times in TNG. The EAGLE simulation also presents different phases of BH-galaxy growth. There is a first phase with a lower relative growth of BHs compared to their host galaxies in low-mass galaxies of M 10 9.5-10 M . Galaxies of M ∼ 10 10 M experience a phase of efficient BH growth. BHs hosted by more massive galaxies of M 10 10-10.5 M have again a less efficient growth than their galaxies. As in TNG, the characteristic galaxy mass between the first phase and the second phase depends on redshift. However, in EAGLE the inefficient BH growth phase is shorter at earlier times than at later times. We do not discuss the shape in SIMBA, since BH seeding takes place in galaxies of M = 10 9.5 M , i.e. where we could expect a change of slope. However, we note that the M BH -M mean relation transitions to a steeper slope at M ∼ 10 10.8 M for z 1.

The regime of galaxies with M 10 10 M is important to understand the build-up of the BH population and the M BH -M scaling relation. As shown below, this regime is subject to resolution effects for the range of resolutions covered in the simulations presented here. In the following sections, we investigate qualitatively the physical processes responsible for the shape of the M BH -M relation and its evolution with time. Before going into details, we mention here that the M BH -M relation can show some imprints of the strength and modelling of stellar feedback. In Habouzit et al. (2017), we discussed the build-up of the BH population in the M BH -M diagram in largescale cosmological hydrodynamical simulations with a physical model of BH formation, and highlighted the ability of SN feedback to stunt BH growth in low-mass galaxies of M 10 9 M (see also Dubois et al. 2015;Anglés-Alcázar et al. 2017c;Bower et al. 2017;McAlpine et al. 2018;Dekel, Lapiner & Dubois 2019). The role of SN feedback was also pointed out in the semi-analytical galaxy formation model MORGANA (Fontanot et al. 2015). It has been shown in several simulations that below the characteristic mass, SNdriven winds are fast enough to overcome the escape velocity of the gravitational potential of the galaxies (Dubois et al. 2015), therefore depleting the centre of galaxies from their cold gas reservoir (which could have fed the BHs). BHs are only able to grow when their host galaxies are massive enough to overcome the effect of SN feedback.

In all these simulations BHs are re-positioned to the potential minimum (of the whole galaxies or within smaller regions) every time-step and thus do not move around in their galaxies, except in Horizon-AGN (see also Habouzit et al. 2017). Doing so favours the accretion on to BHs. Off-centre BHs could also stunt BH growth, as shown in the cosmological simulation NewHorizon (Dekel et al. 2019), and the zoom-in MARVEL-ous dwarf simulations (Bellovary et al. 2019) and the zoom-in FIRE simulations (Anglés-Alcázar et al. 2017c;C ¸atmabacak et al. 2020). Off-centre (or wandering) BHs have been found in recent high-resolution radio observations (Reines et al. 2020).

TNG100 and TNG300 use the same galaxy and BH models but have different resolutions (baryonic resolution of 1.4 × 10 6 M and 1.1 × 10foot_6 M , and spatial resolution of 0.74 and 1.48 kpc, for TNG100 and TNG300), and volumes (100 3 , 300 3 cMpc 3 ). The parameters of the sub-grid models do not depend on mass or spatial resolution.

The better resolution of TNG100 allows us to more accurately resolve the BH surrounding gas density, and thus their gas accretion (Weinberger et al. 2018). As a result, the BHs and galaxies tend to have slightly higher masses in TNG100 at fixed halo mass. We find that the shape of the scaling relation at the low-mass end (log 10 M /M ∼ 9) is different. While we note the presence of a change of slope at all redshifts in TNG300, this is not the case for TNG100 at z = 1-0. The better resolution of TNG100 leads to more efficient growth of the BHs embedded in galaxies of M 10 10 M . Indeed, the simulation is more resolved and the surroundings of BHs have higher gas densities, which leads to increased accretion rates on to the BHs. The normalization of the scaling relation of TNG100 should, therefore, be higher, particularly in the regime where gas accretion is more important than BH-BH mergers for BH growth (M 10 10.5 M ). The resolution also affects the strength of SN winds (Pillepich et al. 2018a); the overall impact of the SN feedback is weaker in the more resolved simulation TNG100. The consequence of the more efficient SN feedback in low-mass galaxies of the TNG model compared to the Illustris model is not as strong in the more resolved simulation TNG100. No other simulation project produced lower resolution volume, thus we cannot assess their level of resolution convergence in a similar manner as the TNG model.

Illustris and TNG100 share the same initial conditions, and therefore one can follow the evolution of individual matched galaxies in the two simulations, and understand how the galaxy/BH sub-grid models affect the build-up of the scaling relation with time. We use the catalogue of Lovell et al. (2014Lovell et al. ( , 2018) that provides the IDs of the matched galaxies. We use the Illustris and TNG100 examples to interpret the evolution of the Horizon-AGN, EAGLE, and SIMBA M BH -M relations.

We compare the final BH mass of all the Illustris and TNG100 matched galaxies in Fig. 6. Our sample includes about 18 000 matched galaxies at z = 0. We find that TNG100 BHs are more massive than their matched BHs in Illustris in low-mass galaxies, i.e. for BHs of M BH 10 8 M . In TNG100, M BH, TNG ∼ 10 8 M is the characteristic mass (yellow dashed line in Fig. 6) for which many BHs transition from the thermal high accretion mode to the kinetic low accretion mode of AGN feedback. The kinetic mode is responsible for efficiently quenching star formation in the BH host galaxies, but also for stunting BH gas accretion. Consequently, it leads to a change of slope 7 in Fig. 6 at this characteristic BH mass of M BH, TNG ∼ 10 8 M . BHs with larger masses have a hard time growing in TNG100, but BHs in the matched Illustris galaxies are still able to do so because of a less effective low accretion mode of AGN feedback. The ability of the Illustris massive BHs to grow when the TNG100 BHs are stunted compensates the lower initial mass of the Illustris BHs. As a result, we find that the most massive matched BHs in Illustris and TNG100 with M BH 10 8.5 M have similar masses in Fig. 6. 10 9 M ). For reference, the BH seeding mass is M BH ∼ 10 5 M in Illustris, and M BH ∼ 10 6 M in TNG100. The black line indicates hypothetic identical mass BH pairs. In low-mass galaxies (i.e. with M BH 10foot_7 M ) BHs are more massive in TNG100. However, when the TNG100 BHs transition to the kinetic AGN feedback mode (upper side of the yellow dashed line), their growth as well as their host galaxies growth is strongly regulated, while the same BHs in Illustris are still able to grow. Consequently, the Illustris BHs catch up their TNG counterparts, and in the most massive matched galaxies TNG100 and Illustris BHs eventually end up with similar masses.

We investigate the time evolution of three different samples of ∼50 matched galaxies each. Given the differences found in Fig. 6, we choose these samples to include galaxies with three key stellar masses at z = 0: (i) Low-mass galaxies: M , TNG, z=0 ∼ 10foot_8 M , corresponding to BHs with M BH,TNG ∼ 10 6-7 M at z = 0, (ii) Intermediate mass galaxies: M , TNG, z=0 ∼ 10 10 M , corresponding to BHs of M BH,TNG ∼ 10 7-8 M at z = 0, (iii) Massive galaxies: M , TNG, z=0 ∼ 10 11 M , corresponding to BHs of M BH,TNG ∼ 10 8-9 M at z = 0.

Because the M BH -M relation is tight in TNG100, the corresponding TNG BH mass ranges that we quote above are within one order of magnitude in BH mass. The M BH -M relation in the Illustris simulation has a larger scatter; the BH mass ranges covered by these galaxy mass samples are larger, between one and two orders of magnitude in BH mass. We first select galaxies that fulfilled our final stellar mass criterion in TNG100, and then look for their matched galaxies in Illustris. We then follow back in time all these galaxies using the galaxy merger trees of the two simulations (Rodriguez-Gomez et al. 2015), which are based on stellar particles and gas cells. We show the time evolution of both their stellar masses (left-hand panels) and their BH masses (right-hand panels) in Fig. 7.

In the following, we start by analysing the stellar mass content of galaxies in the three samples. We find that the relatively lowmass galaxies shown in the top left-hand panel (i.e. M , TNG, z=0 ∼ 10 9 M ) tend to be slightly more massive in Illustris than in TNG100. The TNG100 SN feedback is stronger in low-mass galaxies at all redshifts (addition of velocity floor) and at lower redshifts (new scaling of the wind velocity with the Hubble constant) than in Illustris (see our section 2.2). For the intermediate galaxies (M , TNG, z=0 ∼ 10 10 M , middle left-hand panel), the overall stellar mass content of the matched Illustris and TNG100 galaxies is similar. This is again the case for the evolution of the most massive galaxies (M , TNG, z=0 ∼ 10 11 M , bottom left-hand panel), but we note that some matched galaxies have very different evolutionary paths in this mass range.

While most of the matched galaxies evolve their mass in a similar way, we find that their BHs have different growth histories (righthand side panels in Fig. 7). We set a null BH mass for galaxies that have not been seeded with a BH yet, which translates into the vertical lines in all panels of Fig. 7 for TNG BHs. In the following, we investigate how the Illustris and TNG100 sub-grid models affect the different phases of BH evolution (see the items below) with the right-hand panels of Fig. 7, and what the consequences on the BHs/galaxies co-evolution are with Fig. 8. In Fig. 8 we show the time evolution of the matched BHs and their host galaxies in the M BH -M diagram, dividing in the same three stellar mass bins.

(i) Initial mass of BHs: the initial mass of the TNG100 BHs is higher (by a factor of 10) than the seeding mass of the matched Illustris BHs. As a result, recently formed low-mass galaxies have more massive BHs in TNG100 than in Illustris, as shown in the top right-hand panel of Fig. 7.

(ii) First phase of BH growth: the first phase of BH growth is delayed in TNG100 and takes place earlier in Illustris (middle and bottom right-hand panels of Fig. 7). The stronger feedback from SNe in TNG100 is responsible for delaying the first episode of gas accretion, and not the differences in the seeding or accretion models. Indeed, both Illustris and TNG100 are seeded in dark matter haloes of same mass and thus at about the same cosmic time; the Illustris BHs are not seeded before the TNG100 BHs. At birth, the accretion rates on to the Illustris and TNG100 BHs, given by the Bondi model ṀBH ∝ αM 2 BH , should be similar 8 : the seeding mass (M seed ∼ 10 6 , 10 5 M for TNG100 and Illustris, respectively) and the boost factor parameters (α = 1, 100) compensate one another. The stronger SN feedback is responsible for the change of slope in the M BH -M diagram of TNG100 and the presence of a characteristic mass in the stellar mass range M = 10 9 -10 10 M . With the left-hand panels of Fig. 8, we find that the TNG100 characteristic galaxy stellar mass is moving toward less massive stellar mass with time, and is absent for z 1. Indeed, the first quiescent phase of BH growth is very short in the top panel of Fig. 8 (orange to red coloured lines) for galaxies of M ∼ 10 8.5 -10 9 M , and longer for the galaxies of same mass but at high redshift (bottom left-hand panel, z = 5 -3, green to yellow coloured lines). We find here that the modelling of SN feedback in TNG100 impacts the growth of BHs more substantially at high redshift. The stronger SN feedback in TNG100 is also responsible for the slightly lower normalization of the M BH -M mean relation in TNG100 compared to Illustris, at the low-mass end (M 10 10 M ) at high redshift (see Fig. 5). At high redshift (z ∼ 5) and because of the less effective SN feedback, the Illustris BHs can grow compared to the TNG matched BHs. Therefore, the Illustris BHs embedded in galaxies of M ∼ 10 9 M (e.g. bottom right-hand panel of Fig. 8) are on average slightly more massive, even if the seeding mass of Illustris is smaller than in TNG100.

(iii) Later phases of BH growth: Once the TNG100 BHs start growing, they grow more efficiently than the Illustris BHs, as Figure 7. Comparison of the stellar mass and BH mass in TNG100 (y-axes) and Illustris (x-axes) matched galaxies. The black lines indicate hypothetic identical host stellar and BH masses. Matched galaxies have very similar time evolution of their stellar mass (left-hand panels). The average evolution of the galaxy stellar mass is similar in the two simulations, with the exception that the low-mass galaxies have lower stellar mass in TNG100 (top panels, M , TNG, z=0 ∼ 10 9 M ). BHs in these galaxies are more massive in TNG100 because of the larger seeding mass (right top panel). The initial growth of BHs takes place later in TNG100 than in Illustris. However, when they start growing, they grow more efficiently than in Illustris (middle right-hand panel). However, when the TNG BHs transition to the kinetic AGN feedback mode (bottom right-hand panel), their growth as well as their host galaxies growth is strongly regulated, while the same BHs in Illustris are still able to grow. Consequently, in the most massive matched galaxies (M 10 10.5 M ) BHs have similar masses.

demonstrated in the middle and bottom right-hand panels of Fig. 7.

The difference can also be appreciated in Fig. 8 (e.g. middle panels). While most of the TNG100 BHs grow efficiently, there is a flattening in the M BH -M diagram of Illustris: BH growth is there less efficient than their host galaxy growth, on average. In addition to the potential different gas content in the surroundings of the BHs (see Fig. 10, top panels), there are two important differences in the accretion schemes of the two simulations. Illustris uses the gas content of the BH parent cell to compute the accretion rates, while TNG100 uses a kernel. The accretion rates on to the TNG100 BHs could also be boosted by the presence of the magnetic fields in the simulation (see fig. 8 of Pillepich et al. 2018a, where a higher normalization of M BH -M h is found with magnetic field at z = 0, a higher (or smaller) gas fraction, and smaller (or larger) galaxy size for lower-mass galaxies of a few 10 10 M (or for more massive galaxies)).

All the TNG100 BHs grow efficiently, which leads to little diversity of growth histories in the M BH -M plane in Fig. 8. Instead, the Illustris BHs alternate between episodes of high accretion rates and more quiescent periods with lower accretion rates. BH growth is more diverse in Illustris, at all galaxy mass. Consequently, the scatter of the Illustris M BH -M relation is almost independent of stellar mass (see Fig. 8), and also, on average, larger than in TNG100, except for the phase of efficient TNG100 BH growth that we are discussing here.

The scatter of TNG100 depends on the stellar mass: smaller scatter in the first inefficient BH growth in low-mass galaxies, smaller in the phase where BH growth is regulated by AGN feedback, and larger in the phase of efficient BH growth.

(iv) Towards BH quiescence: The impact of the strong kinetic low accretion mode AGN feedback in TNG100 can be seen in the bottom right-hand panel of Fig. 7. While the growth of the TNG100 BHs is completely stunted because of AGN feedback for BHs of M BH ∼ 10 8 M (horizontal lines in the panel), the Illustris BHs are less regulated by the Illustris AGN feedback (injection of thermal energy in bubbles displaced from the galaxies) and are, therefore, still able to grow. As noticed in Fig. 6, the BHs of these massive galaxies reach similar masses by z = 0.

In Appendix D, we analyse in detail eight pairs of matched Illustris and TNG100 galaxies, with final galaxy mass of M = 10 10 M and M = 10 11 M at z = 0. The phases of BH growth and feedback described above are illustrated for individual galaxies. The properties Reines & Volonteri (2015) (z ∼ 0). The initial growth of the TNG BHs is delayed because of the stronger SN feedback in low-mass galaxies. However, when they start growing, they grow very efficiently. The early growth histories of the Illustris BHs are more diverse and their growth is on average less efficient. There is no strong decrease of BH growth in the massive galaxies of M ∼ 10 11 M in Illustris, whereas their growth is stunted in TNG because of the efficient kinetic AGN feedback.

of BHs and their galaxies are followed in time, and used to understand the build-up of the M BH -M relation with time in Illustris and TNG100, as described in the next section. We provide in Fig. 9 an illustration of the different BH growth histories that two TNG100 and Illustris matched BHs can have.

In Fig. 11, we present a summary of the features discussed in this section with the scaling relation of the simulations at z = 3 and z = 0.

In Illustris, the overall normalization of the M BH -M mean relation decreases with time from z = 3 to z = 0, at all galaxy stellar masses. This can be due to two effects: a less efficient growth of BHs with time and/or a more efficient growth of the host galaxies than their BHs with time. In this section, we use the conclusions from Figs 7, 8 and 10 which shows the gas content of the Illustris and TNG galaxies, the mass of the BHs as a function of their halo mass, and the stellar mass functions of the simulations.

We find that the average mass of gas in the half-mass radius of the galaxies M gas, 1/2 decreases with time. This could favour the first option above (less efficient BH growth at later times), however one has to consider the relative growth efficiency of BHs and galaxies. We note here that the average gas content decreases with a smaller amplitude in Illustris than in TNG100: statistically the Illustris galaxies have more gas than the TNG100 galaxies (for any galaxy stellar mass). Having a decreasing average gas mass does not ensure that the gas content decreases in time for a given galaxy. In fact, when analysing the evolution of individual galaxies (Appendix B, Fig. D1) we find that the gas mass of the galaxies does not evolve significantly with time for z 1 for many galaxies of M ,z=0 ∼ 10 10 M . However, the inner gas content of these galaxies diminishes with time with respect to their total stellar mass and their SFR is constant in the range z 1: galaxies still form stars with the same efficiency and grow their total stellar content more efficiently with respect to their inner gas content. Consequently, galaxy growth becomes more efficient than BH feeding (which scales with the inner gas content, see Fig. D1) with time.

The M BH -M h, DM diagram in Fig. 10 shows very little evolution between the average BH mass and the dark matter mass of their host haloes, indicating that the efficiency of BH growth does not change much with time, on average. The more efficient growth at low redshifts of the galaxies compared to their BHs is also corroborated by the evolution of the galaxy stellar mass function, shown in Fig. 10 (bottom panels). In Illustris, we find that the galaxy mass function keeps evolving in the redshift range z = 1 -0. The mean M BH -M is shifted toward more massive galaxies with time.

TNG100 shows the opposite trend: the overall normalization of the mean M BH -M relation increases with time in the redshift range 4 z 0 and for M 10 10.5 M . Our results for the larger simulation TNG300 indicates that the evolution at high redshift z 3 could extend to larger stellar mass galaxies.

A more efficient growth of BHs at low redshift, or a slower galaxy growth compared to BH growth in galaxies of M 10 10.5 M , could explain the trend identified in TNG100. For TNG100, there is no evolution of the galaxy stellar mass function for M 10 10.5 M and z 1. Moreover, we find a time evolution of the M BH -M h, DM mean relation, suggesting that the time evolution in the M BH -M diagram is driven by changes in the efficiency of BH growth, and that these changes do not affect the host galaxies as much as the BHs. As explained previously, the stronger modelling of the SN feedback is responsible for stunting BH growth at high redshift in TNG100 (see Fig. 11). We find that the impact of SN feedback on BH growth in low-mass galaxies is weaker at low redshift z = 1 -0 than at higher redshifts. As a result, the overall normalization of the mean M BH -M relation increases with time.

In more massive galaxies of M 10 10.5 M and low redshift, the efficient kinetic AGN feedback takes over and stops BH growth via gas accretion. It is responsible for the absence of time evolution in massive galaxies at low redshift. Top panels: Mean mass of the gas contained in the half-mass radius of the galaxies; the entire population at z = 0 is shown in black. The gas content in the inner region of the galaxies diminishes with time. Middle panels: Mean BH mass as a function of the dark matter mass of their host haloes. We show in black all the haloes with BHs at z = 0. At fixed dark matter halo mass, we find a larger evolution in TNG100, and almost no evolution in Illustris. Bottom panels: Galaxy stellar mass functions. Grey symbols show the observational constraints of Baldry, Glazebrook & Driver (2008) (z ∼ 0). The low-mass end of the galaxy mass function is overestimated in Illustris (M 10 10.5 M ); a better agreement is found for TNG100.

The overall normalization of the M BH -M relation decreases with time in Horizon-AGN (Fig. 5 for all stellar masses, except at z = 0 in low-mass galaxies with M 10 10 M whose normalization is higher than at higher redshift.

As in Illustris and TNG100, the amount of gas in the Horizon-AGN galaxies diminishes with time (Beckmann et al. 2017). As shown in the previous section with the TNG100 simulation counterexample, this does not explain on its own the decrease of the normalization with time. One has to consider both the relative growth of the BHs and the growth of their galaxies. To build a complete picture, we analyse the stellar mass function of Horizon-AGN with redshift and compare it to the observational constraints of Baldry et al. (2008) at z = 0 (not shown). 9The stellar mass function of Horizon-AGN cleary overpredicts the observational constraints at z = 0. The simulation produces a higher number of galaxies below the knee of the function (i.e. M 10 11 M ) at z = 0 and at higher redshifts (z 3). It also overpredicts the number of more massive galaxies of M > 10 11 M at z = 0. The stellar mass function evolves towards more massive galaxies from z = 4 to z = 0. We also note that the stellar mass function at z = 0 for lower-mass galaxies of M 10 10.3 M is slightly smaller than at higher redshifts.

We now compare the time evolution of the stellar mass function for given stellar mass ranges, to the evolution of BH mass function for the corresponding BH mass ranges. From Fig. 3 we know that in Horizon-AGN galaxies of M ∼ 10 10.5 -10 11 M host BHs of M BH ∼ 10 8 M at z = 0, and BHs of M BH ∼ 10 7.5 -10 8.5 M at z = 2. Here we anticipate on our following section where we discuss the BH mass function and its evolution with time. The corresponding BH mass function (shown in Fig. 14) does not evolve in the redshift range z = 2-0 for these BHs of M BH = 10 7.5 -10 8 M . Statistically, the galaxies in this stellar mass range grow more than their BHs. This effect can also be seen by comparing Figs 3 and 4 at z = 0. We find a larger fraction of quite faint AGN with L bol 10 42 M in galaxies of M ∼ 10 10.5 -10 11 M (Fig. 3, bottom panel), while a large fraction of these galaxies still form stars with sSFR 10 -10 yr -1 .

Lower mass galaxies with M ∼ 10 10 M host BHs of M BH ∼ 10 6.5 -10 7.5 M at z = 0-2, and BHs of M BH ∼ 10 7 -10 7.5 M at z = 4. The BH mass function for M BH ∼ 10 7 M decreases from z = 1 to z = 0. There are less low-mass BHs of 10 7 M in the simulation at z = 0 than at z = 1. This is due to the non-replenishment of low-mass BHs in the simulation after z = 1.5, when BH formation ceases.

To conclude, we find that the overall decrease of the M BH -M normalization in Horizon-AGN is due to a more important growth of the galaxies than the BHs (suffering from reduced amount of gas in their surroundings), as shown by the galaxy stellar mass function and the BH mass function. As in Illustris, the M BH -M relation of Horizon-AGN is shifted towards more massive galaxies with time. The higher M BH -M normalization for M 10 10 M at z < 1 is explained by the suppression of BH sink particle formation at z = 1.5 in the simulation (see Fig. 11). At fixed galaxy stellar mass, the tail of the BH distribution present at high redshift and composed of newly formed BHs does not exist anymore at z = 0, and leads to a higher mean BH mass value.

The overall normalization of the EAGLE simulation decreases with time for galaxies with M 10 10.7 M , and does not evolve with time for more massive galaxies. The relation stops evolving at z 1 for all galaxy stellar masses, and is unchanged to z = 0. The evolution with time of the M BH -M mean relation has been analysed in detail in Rosas-Guevara et al. (2016), McAlpine et al. (2018), Bower et al. (2017), we briefly summarize their findings.

As in the TNG simulations, the non-linear shape of the M BH -M relation in low-mass galaxies arises from the low efficiency of BH growth, because of the strength of SN feedback (see Fig. 11). When a weaker SN feedback is employed, the relation becomes linear (Crain et al. 2015, their Fig. 10). An even stronger feedback than the fiducial EAGLE model stunts BH growth even more and the simulation does not form any massive BHs of M BH > 10 7 M . As shown in our Fig. 10 with TNG100, the relation between BH mass and dark matter halo mass changes as well with time, demonstrating that the change in the M BH -M normalization is mainly due to a variation in the efficiency of BH growth and not galaxy growth (Rosas-Guevara et al. 2016, their Fig. 6). The change of overall normalization is due to the prevention of BH growth because of SN feedback.

In EAGLE, the characteristic mass at which BHs start growing efficiently, i.e. M ∼ 10 10 M , depends on redshift (Fig. 5). At high redshift, BHs start growing in less massive galaxies. This is interesting because while the TNG simulations have a similar shape because of SN feedback, the evolution with redshift of the characteristic mass follows the opposite trend. McAlpine et al. (2018) show that the transition to the BH efficient growth phase in EAGLE is set by the development of the a hot halo, which traps the SNdriven winds. The transition takes place at fixed virial temperature of T vir ∼ 10 5.6 K, independently of redshift (their fig. 5). The mass of the galaxies reaching this fixed temperature changes with redshift, which cause the evolution of the characteristic mass with redshift.

We also note here that at some level the transition to the rapid BH growth phase also depends on the modelling of BH accretion. In EAGLE, a modified version of the Bondi accretion model is employed, and takes into account the angular momentum of the gas to be accreted by the BHs. Instead of directly falling on to the BHs, the gas first settles into an accretion disc. When a higher viscosity is assumed (i.e. lower C visc parameter in equation 12), the phase of rapid BH growth takes place in less massive galaxies. The phase starts in more massive galaxies for a smaller viscosity parameter (Crain et al. 2015). This shows that the characteristic mass at which BHs start growing efficiently can also depend on the modelling of accretion rate on to the BHs.

The mean M BH -M relation does not evolve at z 2, but its normalization increases with decreasing redshift for z = 2-0 and all galaxy stellar masses (Thomas et al. 2019). The shape of the relation at high redshift resembles the shape identified in TNG. However, since in SIMBA the seeding takes place in galaxies of M 10 9.5 M , the non-linear shape may be due to the seeding rather than SN feedback (see Fig. 11). We note here also that the SN feedback in SIMBA is less efficient (lower wind velocities) than in TNG.

While all the simulations studied here use a version of the Bondi accretion model, SIMBA uses a two-mode model, i.e. a combination of the torque model (for the cold gas) and the Bondi model (for the hot gas). In practice, the torque model dominates at early times (i.e. for most BHs, except the most massive). In that case, the accretion rates on to the BHs is mostly a function of the gas and stellar content in the inner galactic disc, i.e. within a sphere of radius 2 h -1 ckpc around the BHs, rather than the BH mass (only proportional to M 1/6 BH in the torque model, against M 2 BH in the Bondi model). The growth of BHs being tied to the growth of the inner galactic discs, there is no evolution with redshift.

More massive BHs also accrete hot gas via the Bondi model. These BHs have, on average, lower accretion rates and Eddington ratios. For these BHs, the impact of AGN feedback increases as f Edd decreases (AGN wind velocity increases with decreasing f Edd ) and contributes to increase in the amount of hot gas in the galaxies. A Merloni & Heinz (2008), Shankar et al. (2009), andCao (2010). We reproduce the constraint of Merloni & Heinz (2008) for z = 0.1 in all panels (in black). The overall normalization of the BH mass functions increases with time. However, there is no consensus among these constraints on the low-mass end of the BH mass function, or the details of its time evolution. runaway process takes place: increase of BH mass, more important AGN feedback, enhancement of the hot environment, and in return increase of BH mass again via the Bondi accretion channel. In the second paper of this series, we studied in detail the f Edd distributions of all the simulations, but we give here a few key features. At z = 2 a large fraction of the BHs in the mass range M BH = 10 8 -10 9 M enter the regime of low accretion rates of log 10 f Edd -2, while the distribution of lower mass BHs still peaks at log 10 f Edd -1. By z = 0 the majority of the BHs (independently of their mass) have low Eddington ratios (log 10 f Edd -2). Therefore, between z = 2 and z = 0 more and more BHs enter a regime in which they quench their host galaxies, but also increase the hot environment of their close surrounding (gas ejected in jets is heated to the virial temperature of the haloes, assuming that they thermalize their energy into the surrounding hot gas) and favour their growth through the Bondi channel. The less efficient growth of the galaxies and the additional growth channel of BH growth both contribute to increase the overall normalization of the mean M BH -M relation in the redshift range z = 2-0.

In SIMBA, we also note a redshift evolution of the massive end of the M BH -M relation, for M 10 11 M . The Illustris and Horizon-AGN also present an evolution with redshift; there is an overall decrease of the M BH -M relation with time. As we explained already, in these simulations the galaxies keep growing in mass while BH growth is slower, the mean relations actually move towards more massive galaxies. The more efficient AGN feedback of the TNG and EAGLE simulations suppresses both BH gas accretion and galaxy growth, and the growth of both systems is dominated by mergers in this regime, preventing any time evolution of the M BH -M relation. In SIMBA, there is an overall increase of the relation with time. This is due to a more important BH growth via gas accretion with respect to galaxy growth in this regime of massive galaxies (Cui et al., in preparation).

In this section, we present the observational constraints of Merloni & Heinz (2008), Cao (2010), andShankar et al. (2009). BH mass functions are often derived from a given relation between BH mass and a given property of their host galaxies in the local Universe, which can be their total stellar masses, or more likely their stellar velocity dispersion or the luminosity of their bulges. In this case, deriving estimations of the BH mass function at higher redshifts means assuming an evolution of the given scaling relation with time (which is not well constrained at the moment). Other methods combine the estimate of the BH mass function at z = 0 and either an AGN light curve (e.g. Shankar et al. 2009;Cao 2010) or the BH mass function of AGN (e.g. Merloni & Heinz 2008). For both, the underlying goal is to estimate how much the BHs accrete as a function of BH mass and time (i.e. the average growth rate). In the first case, the methods need to state/assume for how long BHs are efficiently accreting, i.e. AGN lifetime/duty cycle, radiative efficiency, which can be a function of BH mass, redshift, BH populations, obscured/nonobscured AGN, etc. Different light curves have been employed. For example, Cao (2010) use a light curve with a power-law decay mimicking the decrease of BH gas accretion due to AGN feedback. In the second case, the methods rely again on scaling relations to derive the distribution of BH accretion rates empirically (instead of assuming a given AGN light curve). The review of Kelly & Merloni (2012) provides a full description of all the methods that have been used in the literature and the pros and cons for each of them.

The three constraints that we show in Fig. 12 do not reach a consensus on either the normalization of the BH mass function, the low-mass end, and do present some differences in their time evolutions. We summarize this below.

(i) Normalization of the BH mass function. All the constraints have an increased overall normalization with time (i.e. decreasing redshift). The normalization of the functions increases at all BH masses in Cao (2010) and Shankar et al. (2009), while the normalization of the low-mass end (M BH 10 7 M ) of the Merloni & Heinz (2008) constraints does not evolve with time. The shapes of the BH mass functions are also different, and this, to some degree, depends on the Eddington ratio distribution employed in the models (Kelly & Merloni 2012;Shankar, Weinberg & Miralda-Escudé 2013).

(ii) Low-mass end of the BH mass function. The constraint of Merloni & Heinz (2008) presents a turn over at the low-mass end for z 1. Instead, the BH mass functions of Cao (2010) and Shankar et al. (2009) more or less reach a plateau at the low-mass end. In these models, these differences are the by-product of the assumptions made when imposing the boundary conditions at the low-mass end (Shankar et al. 2009;Cao 2010).

(iii) Time evolution. While all the overall constraints evolve upward with time, we note some differences in the time evolutions.

The BH mass function of Shankar et al. (2009) evolves gradually in intervals of z = 1, at all BH masses. The massive end of the BH mass function evolves slightly more rapidly in Shankar et al. (2009), and even more so in Cao (2010). In Merloni & Heinz (2008), the BH mass function does not evolve much at high redshift, between z = 4 and z = 3 independently of BH mass. However, a stronger evolution is observed for z 3 (i.e. more or less at the peak of cosmic AGN activity) for M BH 10 7.5 M . The peak of the BH mass functions evolves towards higher mass BHs in the constraints of Merloni & Heinz (2008) and Cao (2010); the identification of a peak in Shankar et al. (2009) is harder. In Merloni & Heinz (2008), the low-mass end does not evolve in the redshift range z = 4-1, while the high-mass end does. The evolution can, therefore, be seen as antihierarchical: the massive end of the BH mass function assembles first, before the full low-mass BH population does (at z ∼ 0). The build-up of the BH mass growth is different in the constraints of Shankar et al. (2009) and Cao (2010). The BH mass functions, instead, build up uniformly for all BH masses up to z ∼ 2. At z ∼ 2 while the massive end is already built and does not evolve much anymore, the low-mass end keeps building up with time up to z = 0.

Given these discrepancies between the observational constraints, one needs to be careful when addressing the differences with largescale cosmological simulations. Moreover, to derive their BH mass functions, the studies described here assume a radiative efficiency, which is different (lower) from what is assumed in the large-scale simulations. For example, Merloni & Heinz (2008) assume a radiative r = 0.07. Since r is largely degenerate with other parameters in the simulations and used more or less as a free parameter to match a given scaling relation, it will not affect strongly our qualitative comparison in the following.

In Fig. 13, we show the simulation BH mass functions. We select all galaxies above the resolution limit M 5 × 10 8 M . We compare the simulation results to the observational constraints presented above: Merloni & Heinz (2008), Cao (2010) (solid grey lines), and Shankar et al. (2009) (dashed grey lines).

The global trends of the observations are reproduced, i.e. that the normalization of the BH mass function increases with time, and less massive BHs are more abundant than their most massive counterpart. While these global trends are in agreement with observations, we do find some noticeable differences with observations. In the following, we start by presenting the results for z > 0 and then we address the comparison between the constraints and the simulations at z = 0. We also discuss the time evolution of the BH mass function obtained in the simulations.

In Merloni & Heinz (2008), the BH mass function peaks at about M BH ∼ 10 7-8 M . We use this peak as a reference and call the low-mass end the part of the function below that peak. We find that the Illustris, TNG100, and Horizon-AGN simulations, overpredict the number of lower mass BHs (M BH 10 7-8 M ). TNG300 overpredicts the number of BHs only for the less massive BHs of M BH ∼ 10 6.5 M . The EAGLE and SIMBA simulations have a lower normalization of the low-mass end. The number of BHs produced is lower than the constraints of Merloni & Heinz (2008) at z 3. EAGLE has a good agreement with the constraints of Merloni & Heinz (2008) at lower redshifts, while we find a good agreement between the simulation SIMBA and the constraints of Shankar et al. (2009) and Cao (2010) at all redshift.

For more massive BHs of M BH 10 8 M , Horizon-AGN and EAGLE provide a good match to the constraints of Merloni & Heinz (2008). Illustris, TNG, and SIMBA seem to form too many massive BHs when compared to Merloni & Heinz (2008), but in good agreement with Cao (2010) and Shankar et al. (2009) for z 1. Interestingly, the simulation SIMBA produces a very good agreement with the BH mass function derived in Cao (2010) (for z 1) and Shankar et al. (2009) for the entire BH mass range and all redshifts. Illustris and TNG produce a slightly higher number of massive BHs with respect to the constraints presented here (at z = 1, 2 for Illustris, and at z 2 for TNG).

Finally, all the simulations except EAGLE overpredict the number of massive BHs with 10 9 -10 10 M compared to the constraints of Merloni & Heinz (2008), for all redshifts.

By z = 0, we find that the overall normalizations of the BH mass functions in the simulations are in good agreement with the constraints of Merloni & Heinz (2008). At higher redshift, Merloni & Heinz (2008) predict a BH mass distribution that peaks at M BH 10 7-8 M , with fewer lower mass BHs. At z = 0 the distribution flattens at the BH low-mass end (presence of more low-mass BHs of M BH 10 7 M at z ∼ 0 than at z 1). We see the same behaviour in Illustris. Most of the other simulations (TNG100, HorizonAGN, and SIMBA) predict a drop of the BH mass function for the less massive BHs, with the exception of the EAGLE simulation, whose shape is very similar to the shape found in Shankar et al. (2009). The low-mass end of Horizon-AGN is the lowest at z = 0; the drop increases between z = 1 and z = 0. The number of BHs with M BH ∼ 10 6 M is reduced by one order of magnitude. This is due to the modelling of BH formation in the simulation; BHs cannot form after z < 1.5, and therefore there are no newly formed BHs to fill the gap at the low-mass end for z < 1.

At z = 0, the constraint of Merloni & Heinz (2008) predicts that one BH of M BH ∼ 10 10 M should be found in a volume of (100 cMpc) 3 (grey line in the bottom panels of Fig. 13). Except for EAGLE, all the simulations produce more than a single BH at this mass. This can also clearly be seen in the bottom panels of Fig. 2, where the number of these massive BHs (located at the top right-hand side of the black line region) is higher in the simulations than in observations.

In Fig. 14, we show the same BH mass functions for the same redshift but on a single panel for each simulation, where it is easier to see the build-up of the BH mass function with time. The black lines indicate the observational constraint of Merloni & Heinz (2008) for z = 0.1. In all the simulations, the normalization of the BH mass function evolves over one order of magnitude. In most of them, the normalization increases for all BH masses with time: quite strongly between z = 4 to z = 1, and then slows down to z = 0. Some of the M BH ∼ 10 7-9 M populations are almost already in place at z = 1.

The homogeneous increase of the normalization in the simulations is more similar to the constraints of Shankar et al. (2009), Cao (2010) Figure 13. Redshift evolution of the BH mass function. Grey shaded areas represent the observational constraints of Merloni & Heinz (2008) for the same redshifts as for the simulations, with the exception of z = 0 for which we show the observations at z = 0.1. We also show the observational constraints of Cao (2010) (solid grey lines) and Shankar et al. (2009) (dashed grey lines). The global trend seen in the observational constraints is reproduced, but some discrepancies are noticeable. At the low-mass end, some simulations produce too many M BH ∼ 10 7 M BHs at some redshift (Illustris, TNG, Horizon-AGN), or not enough (EAGLE), compared to the constraints of Merloni & Heinz (2008). All the simulations (except EAGLE) overproduce BHs of M BH ∼ 10 10 M at z = 0. (Merloni & Heinz 2008). We note for some of the simulations an excess of the most massive BHs. While some simulations do not show a strong change of their normalization between z = 2 and z = 1, this is not true for all of them. than Merloni & Heinz (2008). However, at z 2-1 the massive end of their BH mass functions (M BH 10 9 M ) is already built and do not evolve anymore (antihierarchical build-up), only the lowmass end still assembles. In the simulations, the massive end of the BH mass function (M BH 10 9 M ) keeps evolving after z ∼ 2-1. However, it seems that the simulations also present an excess of the most massive BHs compared to the observational constraints. If the growth of these massive BHs was more regulated, it could lead to a more antihierarchical build-up of the BH mass function.

To understand the build-up of the BH mass function, we divide in Fig. 15 the function by efficiently accreting BHs (also called active BHs and defined as log 10 f Edd -2, dotted lines) and BHs with lower Eddington ratios (log 10 f Edd < -2, dashed lines). The global BH mass functions are shown with grey lines. At early times (z 4) the BH mass function is driven by BHs with high accretion rates in all the simulations. In the EAGLE simulation, the massive end of the function is already dominated by inactive BHs by that time. For the other simulations, the contribution of inactive BHs becomes predominant around z = 2 for the massive BHs, and then gradually covers the lower mass end of the BH mass function. The effect is sometimes called downsizing or antihierarchical growth (while using a different definition than in the previous section): at high redshift the inactive BHs are massive, and at lower redshift inactive BHs encompass more lower mass BHs. In other words, massive BHs at z = 0 are less active than the same mass BHs at earlier times (e.g. Merloni 2004;Greene & Ho 2007;Vestergaard et al. 2008;Vestergaard & Osmer 2009;Hirschmann et al. 2012). At z = 0, we still note some differences among the simulations. While in Illustris, Horizon-AGN, and EAGLE, the contribution from active BHs is low (at all BH masses), which means that the BH mass function is the result of BHs with log 10 f Edd < -2, the contribution of inactive and active BHs are more or less the same for the low-mass end of TNG100. In SIMBA, the massive end is dominated by inactive BHs, but the active BHs still completely take over for the low-mass of the BH mass function (for M BH 10 8.5 M ). The large fraction of active BHs, i.e. BHs with luminosity of L bol 10 42 erg s -1 , for this mass range at z = 0 in SIMBA can be seen in Fig. 3. While comparing the results for Illustris and TNG, we see a very strong decrease in the number of active BHs (dotted lines) at the characteristic mass of M BH ∼ 10 8 M in TNG compared to Illustris. This is a signature of the kinetic AGN feedback mode, which is able to suppress BH accretion in these massive BHs significantly.

The picture that emerges from our analysis is the following. In simulations, BHs grow primarily as active BHs rapidly accreting (i.e. log 10 f Edd -2). Indeed, at high redshift we see that the contribution of non-active BHs (i.e. log 10 f Edd < -2) is modest. The massive end of the BH mass function is built first with active BHs, and is rapidly dominated by non-active BHs as soon as z = 2. The above picture is more or less the same for all the simulations. Still, we have noted some differences among the simulations (as stated above), which indicate that the simulations can have different growth histories for their BHs. The final contribution of active BHs and inactive BHs to the BH mass function at z = 0 is also different for the different simulations. To stress the differences at z = 0, we show the active BH mass function of the simulations in Fig. 16.

A few papers investigate the empirical active BH mass function (Marconi et al. 2004;Greene & Ho 2007;Vestergaard et al. 2008;Vestergaard & Osmer 2009;Kelly et al. 2010;Schulze & Wisotzki 2010), by estimation from broad-line AGN. At z = 0, it appears that the number of broad-line AGN is much lower than the empirical estimates of the full BH mass function, at all masses. At higher redshift (e.g. z = 2), AGN still seem to represent a small fraction of the BHs, except for the most massive BHs. Kelly et al. (2010) find that the number of broad-line AGN with M BH 10 9 M (z = 2) could be high and within the scatter of all the empirical estimates of the BH Figure 16. BH mass function at $z=0$ of the rapidly growing BHs, defined as log 10 f Edd -2. We add the constraints of Schulze & Wisotzki (2010) assuming the same definition log 10 f Edd -2. mass function. In Fig. 16, we also show the observational constraints of Schulze & Wisotzki (2010) which use the same definition for active BHs as the one that we employ for the simulations, i.e. log 10 f Edd -2. 10 The active BH mass function in the constraints is reconstructed: it takes into account the fact that BHs in the observational sample of Schulze & Wisotzki (2010) are flux selected (only AGN above the flux limit are detected) and not BH mass selected. These selection effects primarily impact the low-mass end of the active BH mass function; the uncorrected data of Schulze & Wisotzki (2010) present a decrease of the function towards the low-mass end (which is also present in Greene & Ho 2007, for the same reason 11 ). We find that all the simulations (except EAGLE) have a higher normalization of their low-mass end of the active BH mass function (M BH 10 8 M ). As previously said, we note the complete absence of active BHs among the M BH 10 8 M BHs in some of the simulations (EAGLE, TNG). Horizon-AGN, SIMBA, and Illustris have higher normalization at the high-mass end than the constraints.

While all the simulations follow, by design, one of the empirical M BH -M scaling relations, they struggle to produce the diversity of BHs 10 The definition of active BHs is important, and impacts the normalization but also the shape of the active BH mass function, both in simulations and observations. However, our main conclusions on the evolution and z = 0 contributions of active and inactive BHs to the BH mass function are the same with a different definition such as log 10 f Edd -1. 11 Possible explanations for the discrepancies among the estimates of the active BH mass function can be found in Schulze & Wisotzki (2010). The authors particularly discuss the discrepancies of several orders of magnitude ( = 10 -8 -10 -3 cMpc -3 dex -1 ) for BHs of M BH 10 7.5 M with respect to the studies of Greene & Ho (2007), Greene & Ho (2009), Vestergaard et al. (2008), Vestergaard & Osmer (2009). observed in galaxies in the local Universe (e.g. Reines & Volonteri 2015;Baron & Ménard 2019).

7.1.1 BHs in low-mass galaxies of M ∼ 10 9 M searches for low-mass BHs in dwarf galaxies have led to the discovery of BHs (mostly AGN) of M BH 10 6 M in galaxies of M ∼ 10 9 M (Greene 2012;Reines et al. 2013;Baldassare et al. 2015). The seeding mass in the TNG model is M BH ∼ 10 6 M , and SIMBA starts seeding galaxies of M 10 9.5 M (Fig. 2). By their design, these simulations do not cover the regime of low-mass BHs (M BH 10 6 M ) in low-mass galaxies. The predictions from all the simulations on the low-mass galaxy regime are uncertain due to limited resolution, and seeding prescriptions. While massive BHs of M BH 10 7 M have not been found yet in observations of low-mass galaxies, the TNG, Horizon-AGN and SIMBA (a few only) predict some. While these BHs are inactive in SIMBA, some of them can be qualified as AGN in TNG and Horizon-AGN (Fig. 3). If existing in the local Universe, these massive AGN should have been easier to identify in low-mass galaxies than the lower mass AGN that have been already detected.

7.1.2 BHs in intermediate mass galaxies of M ∼ 10 10-11 M Some BHs of M BH 10 8 M in galaxies of M ∼ 10 10 M have been observed, but are often lacking in the simulations. In EAGLE, all these BHs were found in very compact satellite galaxies that experienced a combination of early stellar assembly and tidal stripping of their stellar component (Barber et al. 2016). In general, the simulations do not produce the lowest mass and highest mass BHs at fixed galaxy stellar mass found in observations in the range M ∼ 10 10-11.5 M . the low-mass AGN (broadline AGN in Reines & Volonteri 2015) of M BH ∼ 10 6 -10 7.5 M observed in galaxies of M ∼ 10 10 -10 11.5 M are not formed in many simulations.

We demonstrated with the BH mass function that there is probably an excess of the most massive BHs (M BH ∼ 10 9 -10 10 M ) in simulations, compared to the constraints of Merloni & Heinz (2008). The most massive of these simulated BHs (M BH ∼ 10 10 M ) are embedded in massive galaxies of M ∼ 10 12 M , a region not covered by Reines & Volonteri (2015). The luminosity function of such observed galaxies with M ∼ 10 12 M at z ∼ 0 is within the range ∼ 5 × 10 -6 -10 -5.5 cMpc -3 dex -1 (Baldry et al. 2008;Bernardi et al. 2013;Muzzin et al. 2013;D'Souza, Vegetti & Kauffmann 2015). All the simulations tend to have extended and tight M BH -M relation after M 10 11.5 M . BHs of M BH 10 10.5 have been found in BCGs (brightest cluster galaxies) with BH dynamical mass measurements (Kormendy & Ho 2013;McConnell & Ma 2013), and up to M BH 10 12 M for BH mass estimates (e.g. mass inferred from the Fundamental Plane, BCG K-band luminosities, Xray, and radio luminosities). The most massive BHs reside in BCGs with bulge mass up to M bulge < 10 12.5 M . A non-negligible fraction of the BHs residing in BGGs and BCGs are overmassive with respect to the bulge stellar mass or the stellar velocity dispersion of their host galaxies (Bogdán et al. 2018;Phipps et al. 2019), compared to the scaling relation of McConnell & Ma (2013) for example. In the comparison of Bogdán et al. (2018) with the Horizon-AGN simulation, a very good agreement is found in the plane M BH -M 500 (M 500 the mass enclosed in a sphere of mean mass density 500 times the Universe critical density). While a relatively good agreement is also found for the M BH -M bulge relation (tighter relation in Horizon-AGN). This tends to show that cosmological simulations produce tighter relations even for the most massive systems when we compare to actual samples of BHs in BCGs/BGGs. Interestingly, a few simulated galaxies were identified in Horizon-AGN with M bulge 10 12.5 M , i.e. more massive than the observed ones (Bogdán et al. 2018;Phipps et al. 2019), and overmassive with respect to their M 500 . These simulated galaxies could be too massive while providing a good agreement for the mass of their BHs. Further comparisons with simulations are needed to understand this regime of very massive systems, and whether the growth of their BHs is triggered in a different way (enhanced number of galaxy mergers, BH mergers, dry mergers, but also low-angular momentum and cold gas flows) than in more normal galaxies.

Two scaling relations are favoured in the local sample of Reines & Volonteri (2015): one for inactive BHs in elliptical and bulge galaxies with dynamical BH mass measurements, and one with a lower normalization for active BHs (i.e. broad-line AGN with masses estimated via reverberation mapping or single-epoch virial theorem). While the two subsamples seem distinct at first glance, Reines & Volonteri (2015) discuss how the quiescent BH population also likely overlaps with the lower AGN relation, but for the most part, they are not detected. The galaxies with pseudo-bulges and dynamical BH masses also overlap with the broad-line AGN. The large sample of Baron & Ménard (2019) with BH mass estimates from narrowline emissions includes a larger number of more massive BHs (M BH ∼ 10 7.5 -10 9 M ) in galaxies of M ∼ 10 10 M . It shows that AGN also populate the region between the two scaling relations for inactive and active BHs of Reines & Volonteri (2015), i.e. BHs of M BH ∼ 10 7 -10 8.5 M in galaxies of M ∼ 10 11 -10 11.8 M . This suggests again that the two previous samples of inactive and active BHs may just be subgroups of a more global population of BHs. The sample of Baron & Ménard (2019) also reinforces our finding that the population of BHs in simulations may not successfully reproduce the diversity of observed BHs (smaller scatter in simulations than observations at fixed M ). Given the low number of objects in the observational samples, incompleteness, selection biases, samples based on AGN (which are rare objects), the uncertainties on BH mass estimates, the scatter of the M BH -M relations probably does not reflect accurately the intrinsic scatter of the BH population in the Universe. In the simulations we consider all the BHs, AGN, or inactive, and we do not take into account the volume of the samples and simulations that we compare. Finally, here the galaxy total stellar mass is used for both simulations and observations, but these computations may not be the same and affect our comparison. For all these reasons, our comparison with observations need to be taken with a grain of salt. As observational samples become larger, we must include volume-weighted comparisons with simulations, and not only compare which regions of the M BH -M observational plane are produced by simulations but also whether the corresponding galaxies have the same types/morphologies as in the observations. Here, we started this process by comparing the star forming properties of the observed and simulated galaxies in the M BH -M plane.

The lower scatter of the M BH -M relation in the simulations may be due to the lack of stochasticity in the galaxy/BH subgrid models. They all use free parameters for the efficiency of physical processes (e.g. radiative efficiency, BH accretion, SN feedback, AGN feedback). These parameters are often set the same for all galaxies, independently of e.g. their mass, evolution, environment, redshift. Instead, these parameters could be more physically motivated, and/or randomly drawn from distributions rather than being fixed values. The latter solution would allow us to investigate the likely larger scatter of the scaling relations.

Instead of seeding all galaxies/haloes reaching a given mass, the seeding could be based on the local properties of the gas such as density (Dubois et al. 2014a) and gas metallicity (Bellovary et al. 2011;Habouzit et al. 2017;Tremmel et al. 2017). The least massive supermassive BH observed in the local Universe has a mass of M BH ∼ 5 × 10 4 M and is located in a M ∼ 2.5 × 10 9 M galaxy (Baldassare et al. 2015). While large-scale simulations of ∼ 100 cMpc do not have the resolution to implement completely physical BH formation mechanisms, they could use seeding mass allowing for the formation of the low-mass BHs observed locally (i.e. M BH = 10 4 M ). Decreasing the initial mass of BHs is also challenging numerically. To follow the dynamics of lower mass BHs and their growth, the seed mass of the BH particles must be sufficiently higher than the mass of the other particles. This allows us to avoid scattering of the BHs with the other particles, i.e. stellar, dark matter, and gas particles (Tremmel et al. 2015;Pfister et al. 2019). The lower the BH seed mass, the higher the resolution of the simulation must be. To account for the stochasticity of BH formation and unresolved early growth in primordial haloes, a distribution of initial BH masses could be employed.

In simulations, BH growth is self-regulated by AGN feedback. A lower AGN feedback efficiency leads to more massive BHs for the same galaxy properties and a higher feedback efficiency to less massive BHs (Di Matteo, Springel & Hernquist 2005;Springel et al. 2005). The radiative efficiency of BH accretion, which is one of the factors determining the strength of AGN feedback, is usually set to a fixed value for all BHs (10 per cent or 20 per cent). The radiative efficiency is closely related to the spin of the BHs, since the spin sets the marginally stable orbit beyond which matter will fall on to the BHs without losing further energy. The BH accretion disc extends closer or further to the BH depending on its spin. Assuming a distribution of radiative efficiencies rather than a fixed value would naturally increase the M BH -M scatter. BH spins are not followed consistently in the simulations studied here. Doing so could both allow us to self-consistently determine the radiative efficiency for every BH, and also increase the scatter especially for the massive end of M 10 10 M (Dubois et al. 2014b;Bustamante & Springel 2019). Consequently, this would affect AGN feedback, since the amount of energy that is released from AGN feedback is higher for highly spinning BHs, and lower for non-spinning BHs. The efficiency of AGN feedback also depends on a variety of other factors that the large-scale cosmological simulations do not resolve, such as the phase structure of the gas into which the energy is deposited. The multiphase ISM is not modelled in most of the simulations, and this can significantly impact both the accretion on to the BHs and the coupling efficiency of AGN feedback, and therefore also the scatter in the M BH -M plane. Better simulation resolutions could help capturing greater starburst and accretion rates on to the BHs, potentially increasing the scatter as well. The modelling of BH accretion itself can affect the scatter, both by the dependence on BH mass (strong dependence in the Bondi model, and weak dependence in the torque model employed in SIMBA), but also by the size of the region considered to compute the accretion rates (parent cell of the BHs, or kernel).

The fact that simulations capture a mean M BH -M relation in agreement with empirical scaling relations but not the full scatter at fixed stellar mass (if this is confirmed by new larger observational samples) may not seem like a major problem at first glance, but it may be. The observational samples of Reines & Volonteri (2015), Baron & Ménard (2019) demonstrate that there is a large number of BHs with M BH 10 8 M in galaxies of M = 10 10.5-11.7 M in the local Universe. These galaxies would be considered as having undermassive BHs in most of the simulations compared to their mean relation. We need to understand whether these observed galaxies represent simply the tail of the BH mass distribution at fixed stellar mass, or if they represent the mean behaviour of the BH population. In the latter case, BHs being less massive than what is found in the simulations would mean that they have experienced less growth and/or mergers and that they have probably injected less energy via AGN feedback in their host galaxies through their history. Would we find the same galaxy properties for these galaxies of M = 10 10. 5-11.7 M in the simulations if they would host a BH of M BH 10 6.5-8 M instead of M BH > 10 8 M ? For example, TNG galaxies quench through the kinetic low-accretion mode of AGN feedback, taking place in galaxies with BHs of M BH 10 8 M . The fraction of quenched massive galaxies may change if a large fraction of these galaxies would host less massive BHs, or we would need to model quenching mechanism differently.

A linear M BH -M relation favours a coherent growth of the BHs and their host galaxies, while a non-linear relation shows that in some galaxy mass regimes the BHs or their host galaxies grow more efficiently than the other. In this paper, we have shown that the BHs formed in simulations with a strong SN feedback, such as EAGLE and TNG (at high redshift), have a reduced growth just after their birth until their host galaxies reach a characteristic mass of 10 9.5 -10 10 M ; after this mass BHs are free to grow efficiently and to catch up the growth of their host galaxies. Only very few AGN candidates have been found in observations at high redshift z 4-5 (Willott et al. 2010;Cowie, Barger & Hasinger 2012;Fiore et al. 2012;Treister et al. 2013;Giallongo et al. 2015;Weigel et al. 2015;Cowie et al. 2020). The lack of AGN at high redshift could be due to high level of obscuration, or to reduced growth of BHs at early times (and so lowermass BHs), in agreement with the shape of the M BH -M relation in TNG and EAGLE. In the upcoming years, combining X-ray missions (i.e. Athena, LynX, AXIS) to JWST will help us to constraint the existence or absence of AGN at high redshift and the properties of their host galaxies (e.g. stellar masses). We discussed the impact of SN feedback on the low-mass end of the M BH -M relation, but it is interesting to note that the simulations studied here also have different implementations of the quasar mode of the AGN feedback. It may, to some extend, impact the normalization and shape of the low-mass end of the relation. The non-linear relationship in low-mass galaxies is a key point to be explored further in observations in the local Universe (Fig. 1). Interestingly, we note that while BH growth is stunted in the low-mass galaxies at all redshifts in EAGLE, this is not the case in TNG100, for which the reduced growth is only observed for z 2 and not in the local Universe (Fig. 5). The low-mass regime is also a key point to explore in cosmological simulations with higher resolution (Habouzit et al. 2017). Higher resolution simulations also suggest that resolving the multiphase ISM is key to capturing the suppression of BH growth by stellar feedback independently of the accretion model (Anglés-Alcázar et al. 2017c, 2020).

As demonstrated in this paper, the shape and the scatter of the M BH -M relation and the BH mass function are all affected by the simulation resolution. With a higher resolution, TNG100 resolves better the gas in the BH surroundings and in the whole galaxies, which explains the higher normalization of the low-mass end of the BH mass function in TNG100 (i.e. the excess of BHs with M BH 10 8.5 M ) compared to TNG300. The lower resolution of TNG300 makes the feedback from SNe more effective, and regulates BH growth more efficiently in low-mass galaxies of M 10 10 M at all redshifts. The main goal of this paper was to carefully compare the BH populations of the last decade cosmological simulations and understand how the subgrid physics affect them, and not to make perfect apple-to-apple comparisons with observations (which is particularly challenging for the scaling relation and BH mass function). We stress here that we have not attempted to understand whether all the observational constraints used in this paper were actually consistent with one another. This is a global concern. Complementary observational constraints often use samples that do not overlap, can be biased in different ways, and assume relations established with other empirical constraints. For example, assuming a scaling relation is in most of the cases needed to determine the BH mass function. The scaling relation may be different from the one we found for a given simulation, but we still compare the simulated BH mass function to the one derived from a different scaling relation. A global effort to address more consistent comparisons is needed.

We have compared six large-scale cosmological hydrodynamical simulations of 100 cMpc length on a side: Illustris, TNG100, TNG300, Horizon-AGN, EAGLE, and SIMBA. We have studied the mass properties of the supermassive BHs of these simulations, more specifically their relations with the stellar mass of their host galaxies, and their mass functions. Our main findings are summarized as follows.

(i) All the simulations produce a M BH -M relation in general agreement with the observations. This is partially by design since they all follow a given empirical scaling relation (e.g. Kormendy & Ho 2013) to calibrate their subgrid models.

(ii) Most of the simulations produce a tight M BH -M relation, and a smaller diversity of BH masses at fixed stellar masses than in the observational samples of Reines & Volonteri (2015) (Fig. 2) and of Baron & Ménard (2019). The most striking difference from observations is that some simulations tend to miss lower mass BHs (M BH 10 7.5 M ) at fixed stellar mass for M 10 11.5 M . (iii) Simulated massive galaxies in the M BH -M diagrams (Fig. 4) have sSFR qualifying them as quiescent galaxies (with sSFR 10 11 yr -1 ), in agreement with observations of Terrazas et al. (2017) and the quiescent elliptical galaxies of Reines & Volonteri (2015). Galaxies hosting less massive BHs covering the region M ∼ 10 10-11.5 M and M BH ∼ 10 6-8 M (i.e. covering or close to the observed AGN region of Reines & Volonteri 2015, that consists of a high fraction of spiral galaxies) are generally star-forming galaxies (with sSFR 10 11 yr -1 ). For most of the simulations, these galaxies also have luminosities typical of AGN (i.e. L bol ∼ 10 43 erg s -1 , Fig. 2).

(iv) The evolution of the median/mean M BH -M relation of the simulations is small in the redshift range 0 z 5 ( 1 dex of log 10 M BH /M , Fig. 5). The normalization of the relation decreases with time in Illustris, Horizon-AGN, and EAGLE, and increases in TNG and SIMBA.

(v) We analysed matched Illustris and TNG100 galaxies and find that the evolution with time of the M BH -M mean relation depends on the subgrid models of galaxy formation and BH physics. The overall decrease of the normalization in Illustris can be explained by a more efficient growth of the galaxies with respect to their BHs (particularly for z 2). For TNG100, the higher normalization with time for M 10 10 M is due to reduced BH growth by strong SN feedback at high redshift in low-mass galaxies.

(vi) The overall decrease of the M BH -M normalization with time in Horizon-AGN is due to a more efficient relative growth of galaxies compared to their BHs with time, as in Illustris. In EAGLE, BH growth is initially reduced by SN feedback (as in TNG) and the rapid growth phase starts in haloes of fixed virial temperature (i.e. in less massive galaxies at higher redshift), which explains the overall normalization increase with time (McAlpine et al. 2018). The normalization of the relation in SIMBA only increases for z 2: the low accretion AGN feedback mode quenches galaxies, but also increases their hot environment, which in return favours an additional growth channel for the BHs (Bondi accretion for the hot gas). The growth of BHs becomes more efficient than for their host galaxies.

(vii) We find that the change of slope of the M BH -M relation arises from the modelling of SN feedback (e.g. its strength, its dependence with time). In TNG and EAGLE the modelling of SN feedback is strong enough to stunt BH growth in low-mass galaxies of M 10 9.5 -10 10 M , creating non-linear relations (Fig. 5). Horizon-AGN and Illustris have weaker SN feedback, leading to linear M BH-M relations. Simulations do not agree on the linearity or nonlinearity of the low-mass end of the relation, which is therefore a key point to be explored in observations in the local Universe and at high redshift.

(viii) The time variation of the M BH -M scatter does not exceed one dex in BH mass for 10 9 M /M 10 11 in all the simulations. In general, the scatter stays the same in Horizon-AGN, increases in Illustris and SIMBA, and decreases in EAGLE and TNG100 with time.

(ix) The scatter of the M BH -M relation has different dependence on M in different simulations (for M 10 11 M ). In general the scatter is constant in Horizon-AGN, slightly increases in Illustris, and decreases in SIMBA. For TNG and EAGLE, the scatter is small when SN feedback efficiently regulates BH growth, larger when BHs grow efficiently, and decreases when they are regulated by AGN feedback.

(x) Different simulations show different behaviours at the highmass end of the M BH -M relation (M 10 11 M ). In TNG300 (largest simulated volume and more statistics in this regime), the relation does not evolve with time, because BH and galaxy growth are driven by mergers in this regime (Pillepich et al. 2018b;Weinberger et al. 2018), and gas accretion is stunted by AGN feedback. The scatter is smaller in this regime and decreases with time, in agreement with predictions from the central-limit theorem. However in SIMBA, BH growth exceeds galaxy growth by mergers in this regime; gas accretion plays a more important role (Cui et al., in preparation). Consequently, the massive end of the relation evolves with time, and its scatter does not decrease.

We analysed the BH mass functions of the simulations, and summarize below our main results.

(i) All the simulations have a BH mass function with an increasing overall normalization with time (Figs 13 and 14), following the hierarchical build-up of their host galaxies.

(ii) Most of the simulations, except the EAGLE simulation, overpredict the number of the most massive BHs (M BH 10 9 M ) at z = 0 compared to the constraint of Merloni & Heinz (2008). These BHs are embedded in the most massive galaxies of M ∼ 10 12 M (Fig. 2). The TNG100 simulation also predicts an excess of lower mass BHs with M BH = 10 7 -10 9 M (probably due to the high seeding mass).

(iii) We find that the most massive BHs are less active at z = 0 than the same mass BHs at earlier times in all the simulations, and therefore all the simulations show evidence of an antihierarchical growth.

(iv) We identify some differences in the contribution of rapidly growing BHs (active BHs with log 10 f Edd -2) and inactive BHs to the BH mass function among the simulations. At high redshift (z = 4), the BH mass function is dominated by active BHs for most of the simulations, except EAGLE with already a large contribution from inactive BHs. Inactive BHs dominate the massive end of the mass function at z = 2, and the entire mass function by z = 0. However in SIMBA, active BHs still dominate the BH mass function for M BH 10 8.5 M .

(v) At z = 0, all the simulations (except EAGLE) overpredict the number of active BHs at the low-mass end (M BH 10 8 M ) compared to the constraints of Schulze & Wisotzki (2010). At the high-mass end, half of the simulations do not have any active BHs of M BH 10 8 M , the other half have higher normalization of more than an order of magnitude.

While the scaling relation M BH -M is the most convenient because stellar mass is the easiest quantity to measure both in observations and simulations, it is not the tightest relation that has been found in observations. The relation with bulge mass has a smaller scatter (Magorrian et al. 1998;Häring & Rix 2004;Gültekin et al. 2009), and the relation with stellar velocity dispersion of the host galaxies an even tighter scatter (Ferrarese & Merritt 2000;Gebhardt et al. 2000;Kormendy & Ho 2013;McConnell & Ma 2013;Woo et al. 2013;Saglia et al. 2016;Shankar et al. 2016;de Nicola et al. 2019). Measuring the mass of bulges in observations is challenging and difficult to do in a systematic way for all galaxy morphologies. Bulge + disc decomposition relies on several aspects/assumptions (each study having its own) and can suffer from lack of information on the inclination, surface brightness limits, spatial resolution, and signal-to-noise ratio (see discussions in Simard et al. 2011;Dubois et al. 2012;Volonteri et al. 2016). It has been shown that the M BH σ plane could be more fundamental than those with galaxy stellar luminosity or mass, or effective radius (for a review Marsden et al. 2020). At relatively low redshift, there is no observational evidence for a redshift evolution of the M BHσ relation, as demonstrated by Shen et al. (2015) using about 90 SDSS quasars at z < 1 and by indirect estimates (Shankar, Bernardi & Haiman 2009a;Zhang, Lu & Yu 2012).

Regarding cosmological simulations, Li et al. (2020) recently studied the relation M BHσ in Illustris and TNG100. The σ values were computed using rest-frame SDSS-r band luminosityweighted stellar line-of-sight velocity dispersion measured within a given projected radius, as it is the case for observations. They demonstrated some discrepancies between the simulations and the observations, where BHs in Illustris and TNG100 were either too massive at fixed σ , or conversely σ were too small with respect to the BHs (Li et al. 2020). The same M BHσ scaling relation was also studied for Illustris (Sijacki et al. 2015), Horizon-AGN (Volonteri et al. 2016), EAGLE (van Son et al. 2019) and SIMBA (Thomas et al. 2019), but without an observationally friendly frame. Interestingly, SIMBA shows a broader scatter in M BHσ than for the other scaling relations (Thomas et al. 2019).

In Fig. B1, we show the simulation M BH -M diagrams colour coded by the specific star formation rate of the BH host galaxies. , 0 for all the simulations. We colour-code the simulated BHs with the specific star formation rate (sSFR) of their host galaxies. We set sSFR = 10 -13 yr -1 for galaxies with lower sSFR. The limit sSFR = 10 -11 yr -1 is often used to define star-forming galaxies (with higher sSFR) or quiescent galaxies (with lower sSFR). In the figure, galaxies with blue colours are forming stars, and yellow to red galaxies are quiescent. We show in the last right-hand panel the observed star-forming and quiescent galaxies of Terrazas et al. (2017) for references, colour-coded by their sSFR as well. To guide the eye, we again reproduce the black region representing the observational sample of Reines & Volonteri (2015) at z = 0 in each panel.

We present the scatter of the M BH -M relations for the simulation Illustris, Horizon-AGN, and EAGLE in Table C1, for different bins Table C1. Time evolution of the scaling relation M BH -M . We provide the stellar mass bins (log 10 M /M ), the redshift in the range z = 3 -0, the corresponding mean ( log 10 M BH /M ) and median (log 10 MBH /M ) of BH mass in the given stellar mass bin, and the scatter (15th-85th percentiles) of the mean scaling relation. For reference, we find that the 15th-85th percentiles of the observations of Reines & Volonteri (2015) (at z ∼ 0) is comprised in the range σ = 0.9-1.9 for M = 10 9.5 -10 11 M . For the observational sample of Baron & Ménard (2019), we find σ = 1.2-1.6 in the same M range. The simulations produce a smaller scatter than found in the observations.

M redshift M BH MBH σ Illustris 9.5 z = 3 6 . 7 6 6 . 7 7 0 . 5 4 9.5 z = 2 6 . 6 2 6 . 6 2 0 . 5 9 9.5 z = 1 6 . 3 6 6 . 3 6 0 . 7 4 9.5 z = 0 5 . 9 0 5 . 8 6 0 . 8 8 10 z = 3 7 . 3 6 7 . 3 5 0 . 6 0 10.5 z = 3 7 . 9 0 7 . 9 2 0 . 4 5

. Time evolution of four pairs of matched Illustris (dashed lines) and TNG100 (solid lines) galaxies with M gal, TNG, z=0 = 10 10 M . First columns: BH mass and host galaxy stellar mass evolution (top), and the ratio of these quantities M BH /M gal (bottom). Second columns: ratios M gas, 1/2 /M gal of the gas mass in the half-mass stellar radius to the stellar mass in twice the half-mass radius. Third columns: galaxy sizes (R 1/2 ), and star formation rates (SFR). Bottom panel: M BH -M gal diagram. In these galaxies evolving to M gal = 10 10 M at z = 0, we note that a large fraction of their BHs are more massive in TNG100 than in Illustris. This result holds even for the BHs embedded in galaxies with similar stellar mass evolution; the right-hand panels show that these galaxies have however different evolution of their sizes and SFR history with time. The TNG100 galaxies are seeded with massive BHs of ∼ 10 6 M . The start of the efficient BH accretion phase is delayed. After the low accretion phase, the TNG100 BHs grow more efficiently than their galaxies, as shown by the increasing M BH /M gal ratios with time. In Illustris, we note a large fraction of systems with decreasing M BH /M gal ratios with time, reflecting the more efficient growth of the galaxies compared to their BHs. In this mass range, most of the matched BHs in the two simulations are more massive in TNG100 than in Illustris.

parent cell as in Illustris), but also to the addition of the magnetic fields, which boosts the accretion in TNG100. We stress here that at fixed environment (galaxy properties, gas content, etc) and without considering the above factors, the accretion rates on to the TNG100 BHs will be more important in TNG100 as soon as the mass of the TNG100 BHs is one order of magnitude larger than their matched Illustris BHs (due to the choice of the parameters in the Bondi modelling).

The largest differences in the evolution of the BH mass (up to one order of magnitude in M BH ) is found for the stellar mass regime M , TNG(z=0) = 10 10 M . The evolutionary paths of BHs in matched galaxies that reach M , TNG, z=0 ∼ 10 11 M at z = 0 is more diverse than for galaxies with M ,TNG, z=0 = 10 10 M . We illustrate this with Fig. D2. While the BHs tend to be always more massive when reaching z = 0 in TNG100, we find some periods of time for which the matched Illustris BHs are more massive.