Super-Eddington accretion in high-redshift black holes and the emergence of jetted AGN

In this work we study the co-evolution of central black holes (BHs) and host galaxies by utilizing an advanced iteration of the DELPHI semi-analytical model of galaxy formation and evolution. Based on dark matter halo merger trees spanning the redshift range from $z=20$ to $z=4$, it now incorporates essential components such as gas heating and cooling, cold and hot BH accretion, jet and radiative AGN feedback. We show how different BH growth models impact quasar and galaxy observables at $z \geq 5$, providing predictions that will help discriminate between super-Eddington and Eddington-limited accretion models: despite being both consistent with observed properties of SMBHs and their host galaxies at $z \sim 5-7$, they become very clearly distinguishable at higher redshift and in the intermediate mass regime. We find that the super-Eddington model, unlike the Eddington-limited scenario, predicts a gap in the BH mass function corresponding to the intermediate-mass range $10^4 \mathrm{M_\odot}<M_\mathrm{bh}<10^6 \mathrm{M_\odot}$. Additionally, it predicts black holes up to two orders of magnitude more massive for the same stellar mass at $z=9$. The resulting velocity dispersion - BH mass relation at $z \geq 5$ is consistent with local measurements, suggesting that its slope and normalisation are independent of redshift. Depending on the Eddington ratio, we also model the emergence of AGN jets, predicting their duty cycle across as a function of BH mass and their potential impact on the observed number density distribution of high-redshift AGN in the hard X-ray band.


INTRODUCTION
The conspicuous presence of supermassive black holes (SMBHs) at z > 6 − 7 (see for instance Mortlock et al. 2011;Wu et al. 2015;Bañados et al. 2018;Matsuoka et al. 2018;Wang et al. 2021;Larson et al. 2023) implies that their seeds must have grown in mass by up to 7−8 orders of magnitudes in just a few hundreds million years.This can be explained thanks to massive BH seeds formation mechanisms or very efficient early growth model (see Volonteri et al. 2021, for a review), and possibly a combination of the two.From this perspective, Bromm & Loeb (2003) first proposed the possibility of forming direct-collapse black hole seeds with masses ⋆ piana@ntnu.edu.twM bh ∼ 10 4 − 10 5 M⊙, but their specific environmental requirements are supposed to make them rare objects, and so far there is no observational evidence of their existence.For this reasons, models of super-Eddington BH growth via radiatively inefficient slim accretion disks are more and more often being regarded as the solution to this conundrum (e.g.Pezzulli et al. 2016).Simulations have shown that the feedback from mild super-Eddington accretion rates interferes with the accretion flow itself, making the growth process discontinuous and overall rather inefficient (Regan et al. 2019), unless the spin of the black hole remains low, reducing the feedback efficiency (Lupi et al. 2023).If on the other hand the black hole can enter a regime of hyper-Eddington accretion with Ṁbh > ∼ 500 ṀE, then a prolonged isothermal steady accretion flow can form, favouring a much faster growth (Inayoshi et al. 2016;Sugimura et al. 2017;Takeo et al. 2018).Still, observed SMBHs with M bh > ∼ 10 7 M⊙ do not seem to show any indication of undergoing such strong accretion episodes (e.g.Trakhtenbrot et al. 2017), suggesting that these might be characteristic of earlier phases of BH growth.If so, it is reasonable to expect the huge amount of feedback energy released during a hyper-Eddington accretion episode to leave some imprints in the statistical properties of the young host galaxies.Indeed, the observed correlations between the central black hole mass and the mass, luminosity and velocity dispersion of the galactic bulge (Ferrarese & Merritt 2000;Gebhardt et al. 2000;Marconi & Hunt 2003;Gültekin et al. 2009;Graham 2014) hint to a co-evolution between the BH and its host galaxy (see Shankar 2009;Kormendy & Ho 2013, for a review), and the widespread detection of gas outflows in galaxy with active BHs (e.g.King & Pounds 2015) allow us to think that the mechanism responsible for such co-evolution is the coupling between supernovae (SN) and BH feedback energies and the interstellar medium (ISM).In particular, SN feedback is effective in removing gas from the central regions of the galaxy, especially during the initial phases of galaxy evolution, when the mass is low and the potential well shallow, hence hampering black hole growth (Bower et al. 2017;Lupi et al. 2019).BH feedback, on the other hand, becomes important at higher masses, and can usually take two different forms (see Morganti 2017;Cielo et al. 2018, for comprehensive reviews): the radiative (quasar ) mode is usually associated to high-luminosity AGN with high accretion rates, and powered by radiation emitted from the accretion disk (Shakura & Sunyaev 1973), with photons that are able to transfer their momentum to the IGM (intergalactic medium) particles.These AGN are usually associated to fast outflows with high velocities (> 500 km/s), and their bolometric luminosity correlates with the outflow mass and size (Fiore et al. 2017), suggesting that the AGN radiative feedback might indeed be the main driver of outflows.The jet (sometimes called radio) mode, instead, is supposed to extract energy from the magnetic field of rotating black holes (Blandford & Znajek 1977), and is considered to be dominant in low-power AGN, though it has also been occasionally observed in luminous quasars.To be more precise, it is thought that for Eddington ratios λE = Ṁbh / ṀE < ∼ 0.01 (where Ṁbh is the black hole accretion rate and ṀE is the Eddington accretion rate) the black hole accretion disk is geometrically thick and optically thin, allowing the formation of a weak but steady jet.At higher accretion rates, for 0.01 < ∼ λE < ∼ 1, the disk enters an unstable phase and the gas inflow towards the black hole becomes time-dependent.If λE > ∼ 1 then we have a strong optically-thick radiation-pressure-driven wind which again thickens the disk into a torus, producing jets which will be more collimated for wider tori (for a review about relativistic jet formation in AGN see Blandford et al. 2019).The effect of jets, which are observed in ≈ 10% of all AGN, goes well beyond the boundaries of the galaxies: their radio cavities filled with plasma outgrow the host halo, preventing the intergalactic medium (IGM) in the cluster from cooling down and pushing it outwards, hence disrupting gas accre-tion onto the galaxy.The importance of this effect on the IGM is still unclear though, as it can be degenerate with that of stellar feedback.In addition, the amount of energy deposited into the ISM by the jets can vary a lot, depending on their collimation and on the clumpiness of the medium (Wylezalek & Morganti 2018;Wagner et al. 2012;Mukherjee et al. 2016).
In this work, we use a refined version of the DELPHI cosmological semi-analytic model (Dayal et al. 2019;Piana et al. 2021Piana et al. , 2022) ) -benchmarked against galaxy and AGN UV luminosity functions, stellar mass density and UV luminosity density at z > 4 -to probe the observable properties of the high-redshift AGN population in different scenarios (see Trebitsch et al. 2023, for an alternative version of DEL-PHI).In particular, in addition to the previous version of DELPHI, we model the mass budgets of the cold and hot gas phases, which in turn fuel cold and hot black hole accretion (e.g.Raimundo et al. 2017;Storchi-Bergmann & Schnorr-Müller 2019), AGN radiative and jet-mode feedback, and the galactic gas cycle of outflows and re-accretion onto the galaxy.The goal is to show how we can observationally distinguish between Eddington-limited and super-Eddington accretion models, helping in defining what is the most typical BH growth path.One of our objectives is to predict the emergence of jets at z > 4, and to assess their impact on the observational properties of AGN at high-redshift.More specifically, we will determine the jet duty cycles as a function of black hole mass, and how the the AGN number densities evolve over the redshift and luminosity space.Hard X-rays from distant AGN are less affected by line-of-sight effects, and also at these energy bands the contribution from the host galaxies is generally insignificant.We will therefore compare our results with reference to those by keV X-ray surveys (Ueda et al. 2003;Barger et al. 2005;Ueda et al. 2014;Miyaji et al. 2015;Aird et al. 2015).

MODEL
Galaxies evolve as their host dark matter halos grow.Our semi-analytic model is built on the merger trees of 550 halos with final masses M h = 10 8 − 10 13.5 M⊙, whose merger history and mass evolution are followed from z = 20 to z = 4 in time steps of 20 Myr (though we will test some of our results also on a merger tree with a time step of 10 Myr; see also §2.6).The merger tree algorithm follows the one described in Parkinson et al. (2008).Each of the 550 final halos is assigned a number density consistent with the Sheth-Tormen halo mass function (HMF, Sheth & Tormen 1999) at z = 4, and all of its progenitors along the merger tree are then assigned the same number density as the final halo, so to reproduce the correct HMF at each redshift.The merger tree has a mass resolution of 10 8 M⊙, which constrains the mass with which new halos form along the merger tree.The newlyformed halos are immediately seeded with gas mass, proportionally to the cosmic Ω b /Ωm ratio.The starting leaves of the merger trees can also be seeded with black holes: these seeds are assumed to be 10 .Schematic view of a dark matter halo merger tree branch: progenitors at z + ∆z of a halo at z each bring their own contribution of dark matter, (hot and cold) gas, stellar and black hole mass.In addition, smoothly accreted dark matter and gas from the intergalactic space and from the gas reservoir around the halo are accreted according to equations 1 and 2.
and the Lyman-Werner (LW) background impinging on the halo is 30J21 (Dayal et al. 2017), where J21 is the LW background expressed in units of 10 −21 erg s −1 Hz −1 cm −2 sr −1 (Sugimura et al. 2014).Starting leaves at z > 13 not meeting these criteria are instead seeded with a 150 M⊙ stellar black hole, resulting from the collapse of Pop-III stars in primordial mini-halos (SBH; for a review about BH seeds formation channels see for example Latif & Ferrara 2016;Inayoshi et al. 2020;Volonteri et al. 2021).It is worth to point out that the evolving BH population is fully dominated by descendants of SBHs, since the number density of DCBH seeds is 2-3 orders of magnitudes lower (Dayal et al. 2019).In this section we describe our treatment of the evolution of the dark matter and baryonic mass components along the merger trees, from time step to time step.

Overview
Since gas in galaxies represents the fuel for both star formation activity and black hole growth, it is then essential to follow accurately the evolution of the different gas phases to describe galaxy evolution.In order to do so, at each time step we model the mass changes of the following key components: dark matter halo, hot gas, cold gas, black hole, stellar mass, and the halo gas reservoir, in addition to the gas accreted from the IGM.Schematic plots illustrating the relationships between the components are presented in Figures 1 and 2.
Beside the mass contributions from all of its progenitor, each halo at a time step z will accrete an unresolved amount of dark matter from the intergalactic space according to where the sum runs over all the j progenitors at z + ∆z.Together with the unresolved dark matter, the galaxy accretes from the IGM a gas mass proportional to the cosmic baryonic fraction.The progenitors of a halo will bring in also their content of gas, stars and black hole mass.We can then write where the index j runs over all the progenitors of the halo at the previous time step, and τs is the time step employed in our model, which in our fiducial case is 20 Myr.In the later subsections we will omit the dependence on z from the equations, for simplicity, and we will assume the different terms all refers to the same time step z, unless otherwise specified.

The gas phases
The differential change of the cold gas mass within a time step is modelled by: In this equation, Ṁ cold acc 1 is the cold gas accretion rate onto the galaxy, Ṁcool is the gas cooling rate, Ṁheated is the gas heating rate because of AGN feedback, Ṁ * is the star formation rate, Ṁ cold bh is the cold gas mass accretion rate of the black hole, Ṁ ej * and Ṁ ej bh represent the gas ejection rates from the galaxy by SN and AGN feedback respectively.
Conversely, the change in the amount of hot gas mass is given by Similarly to what we did for the cold gas, Ṁ hot acc is the hot gas accretion rate onto the galaxy and Ṁ hot bh is the hot gas black hole accretion rate.From equations (3) and (4) notice that while we are assuming that stars are formed only from cold ISM (interstellar medium), black holes can accrete both from the hot and cold ISM.In addition, we assume that the effect of SN feedback is solely to drive gas outflows from the host galaxy, while black hole feedback has the potential of  1).See text for the definition of all the terms included here (equations 3 and 4).The hot gas and gas reservoir components represent additions to the model presented in Dayal et al. (2019) and Piana et al. (2021).
both driving gas outflows and heating up part of the cold ISM, depending on whether the jet is active or not.
The terms M cold acc and hot M hot acc represent the cold and hot gas masses accreted from the IGM and from the gas reservoir formed around the halo.Simulations have shown that galaxies will accrete preferentially hot or cold gas depending on their halo mass.If the halo mass M h is lower than a critical value M crit h ∼ 10 12 M⊙, there is no virial shock surrounding the halo and the gas is accreted from the IGM onto the galaxy in cold mode (Kereš et al. 2005;Dekel & Birnboim 2006;Ocvirk et al. 2008).Above the same critical halo mass, on the other hand, pressure builds up and a virial shock develops, leading to quasi-spherical hot gas accretion mode.In this case though, accretion of cold gas filaments travelling inwards from large scales is still possible (for a review see Dayal et al. 2019).In our model, for simplicity, we introduce a parameter f cold = 0.6 that represents the fraction of gas that is accreted cold.We checked that if we implement a linear halo mass dependence of f cold that follows the trend recovered by numerical simulations (see Figure 4 of Ocvirk et al. 2008), the change in our results is negligible. in our model, the total gas accretion rate onto the galaxy is defined as where the term Ṁret represents the return rate of coronal gas mass Mres that falls back onto the galaxy from the reservoir, and reads Ṁret = αretMres/τ dyn , with τ dyn = Rvir/Vvir, being the halo dynamical timescale an αret a free parameter.The gas reservoir Mres is formed by the gas ejected by SN and AGN feedback.We then define and At each time step z a part of the total hot gas mass in the galaxy is then allowed to cool down at a rate Ṁcool .It is generally assumed that the cooled-down mass corresponds to the gas mass enclosed by the cooling radius, within which the cooling timescale is shorter than the free-fall timescale.We need then to impose a density profile for the hot gas, which we assume to follow an isothermal distribution where Rvir is the virial radius of the dark matter halo, defined according to Barkana & Loeb (2001) as Here ∆c = 18π 2 + 82(Ωm(z) − 1) − 39(Ωm(z) − 1) 2 .At each time step we compute the cooling timescale at which the hot diffuse gas is able to cool, defined as where µ = 0.59 for a fully ionised primordial gas, mp is the proton mass, kB the Boltzmann constant, Tvir the virial temperature of the halo and Λ(T, Z) is the cooling function as computed by Sutherland & Dopita (1993).In our case we assume all of our galaxies have metallicity Z = 0.05Z⊙.Once we define the virial velocity of the halo as we can compute the corresponding virial temperature The portion of gas that can fuel black hole or stellar growth corresponds to the amount of gas that has enough time to cool and to fall to the centre of the potential well, and that should be contained within both the cooling radius and the free-fall radius.If we assume that the cooling timescale at the cooling radius is similar to the halo dynamical timescale τ dyn = Rvir/Vvir, we can derive the cooling radius by equating the two quantities, obtaining Similarly, to estimate the free-fall radius, we can compare the dynamical timescale with the free-fall timescale, defined as where ρ is the total density and f hot = ρ hot /ρ the mass fraction of the hot gas component.Note that here we are assuming that the hot gas density profile tracks that of the total halo mass.By imposing τ ff = τ dyn we define the corresponding free-fall radius The accretion radius, within which all the hot gas cools down and is funnelled towards the disk, corresponds then to and is evaluated at each time step.We can then write the instantaneous gas cooling rate

Star formation
Star formation occurs in molecular clouds, so it arises from the cold phase of the gas in the galaxy.At each time step, after implementing the gas cooling mechanism, a fraction f eff * of the cold gas mass forms new stars, and the feedback of the star-forming activity contributes to photo-evaporate part of the remaining cold gas out of the host galaxy.The effective star formation efficiency f eff * is defined as the minimum value between the star formation efficiency whose corresponding SN-II feedback is enough to expel the rest of the gas from the host galaxy and an upper threshold value f * = 0.02.We define ESN = f w * E51ν Ṁ * τs = f w * vs 2 Ṁ * τs as the energy produced by supernovae each time step.Here, E51 = 10 51 erg is the energy imparted onto the ISM by each SN-II explosion and ν = [134M⊙] −1 is the number of SNII per stellar mass formed by a Salpeter initial mass function between 0.1 and 100 M⊙; f w * is the coupling factor between the SN energy and the gas and vs is computed to be 611 km/s.We can then write the star formation rate at time step z as The energy required to unbind the cold gas remaining after the star formation burst is where the escape velocity can be written in terms of the halo rotational velocity as ve = √ 2 Vvir.Equating ESN and Eej we obtain the feedback-limited star formation efficiency The effective star formation efficiency then reads while the gas mass ejected by SN feedback in one time step is

Black hole accretion
The central black hole grows through mergers and accretion of both hot and cold gas.Because of its limited physical size, its accretion rate strongly depends on the environment and on the gas supply in the centre of the galaxy.In particular, Bower et al. (2017) made the point that SN feedback in low-mass halos are very effective in driving cold gas outflows away from the galactic centre, hence starving the central black hole.In such cases, the cold gas bubbles heated by the SN energy are hotter than the gas in the external regions of the galaxy, and, being buoyant, will travel outwards.In high-mass halos, on the other hand, galaxies have a higher virial temperature, and these same bubbles are not buoyant anymore.The gas remains then in the central region of the galaxies, providing the fuel for black hole growth.They estimated the threshold halo mass above which SN feedback is not effective anymore to be where ∆z = with what is shown also in numerical simulations (Rosas-Guevara et al. 2016;Lupi et al. 2019).Since what marks the end of the stunted black hole accretion regime, physically speaking, is the rise of the halo virial temperature, and given that the two values are found to be similar, we assume M th h ≡ M crit h .We treat this critical halo mass as a free parameter of the model, and tune it to be and we write the mass accreted by the black hole at each time step as where the inflow of sparse hot gas towards the central black hole is allowed at all times.
In our model, as in many others, hot gas accretion onto the central black hole takes the form of Bondi-Hoyle-Lyttleton accretion mechanism.In this case we have that the hot gas mass accretion rate of the black hole is where ρ bh is the density of the gas surrounding the black hole and cs is the sound speed, here approximated by the halo virial velocity Vvir (Croton et al. 2016).To find ρ bh we use the so-called maximal cooling flow model (Nulsen & Fabian 2000), and we equate the sound travel time across a shell of diameter twice the Bondi radius to the local cooling time, obtaining Episodes of cold gas accretion onto the black hole instead are assumed to be triggered by major mergers, characterised by a halo mass ratio and to last until the cold gas mass fraction mc = of the new host has decreased below a fraction fc of its value at the moment of the merger.For the duration of the accretion episode, the black hole is allowed to accrete a fixed fraction of the total cold gas mass present in the galaxy after the star formation burst has taken place where we defined

AGN feedback
In our model we assume that gas outflows are launched by radiative feedback during cold gas accretion while jets contribute to heat up the cold gas in the galaxy.Hence, at each time step, the ejected gas mass can be written as where f w qso and f w jet are coupling constants that describe how much of the quasar and jet energy couples to the gas, and for simplicity are tuned to be equal.To compute the total quasar energy Eqso = Lqsoτs emitted in a time step, we employ the solution as computed from simulations of relativistic slim accretion disk of Sadowski (2009) and fitted by Madau et al. (2014), which takes into account the spin of the black hole and is applicable also in the case of super-Eddington accretion.In this formulation the Eddington rate is defined as ṀE = 16 LE/c 2 , with LE being the Eddington luminosity.Given our accretion rate Ṁbh , we can then compute the bolometric luminosity Lqso emitted by the black hole as (34) Here a = 0.5 is the dimensionless spin parameter, taken to be equal for all black holes.This fit is shown to yield acceptable residuals with respect to the numerical results within the ranges 0 < a < 0.998.Physically, this means that the radiative efficiency will be lower for black holes characterised by higher (super-Eddington) accretion rates and higher spins.
Finally, we assume that the jet mode is turned on .Mass evolution for different galaxy components for 1phase model (no gas heating and cooling, solid lines), compared with the fid model (dashed lines), here used as a reference.In each panel we follow the growth of the main branch of the merger tree of halos of different mass.We show the total gas mass content (M gas tot ) at each time step, the cumulative black hole and stellar mass (M bh and M * ), the cumulative gas mass ejected by AGN and stellar feedback (M ej bh and M ej * ).The vertical grey lines indicate all the time steps at which a major merger occurred.
only when λE ⩽ 0.01 or λE ⩾ 1, and its power is computed according to the Blandford-Znajek power defined in Tchekhovskoy et al. (2011) where ϕ is the dimensionless magnetic flux and f(a) = a 2 (1 + √ 1 − a 2 ) −2 .In this work for each black hole seed we randomly draw ϕ from a uniform distribution of values between 1 and 50.Every time there is a black hole merger, the resulting ϕ will be that of the black hole with higher mass.Part of the jet energy goes into heating up cold gas present in the galaxy, according to the equation where f h jet represents the fraction of jet energy that goes into heating up the gas.
The gas masses ejected in the form of outflows by SN and AGN feedback are accumulated into the gas reservoir (Mres), whose change within a single time step is then computed as In each panel we follow the growth of the main branch of the merger tree of the same halos as in Fig. 3.We show the cold and hot gas mass content (M cold and M hot ), and the cumulative heated (M heated ) and cooled (M cool ) gas masses.Notice that the last two quantities and the hot gas mass content are defined only for the fid case.The vertical grey lines indicate all the time steps at which a major merger occurred.

Model setup and parameters
We summarise the free parameters and the adopted values in Table 1.To compare the effects, several different versions of the model are considered, and their differences are outlined in Table 2.In the fiducial (fid) model, we allow both hot and cold black hole accretion, with no explicit Eddington limit.The dark matter merger trees have a time resolution of 20 Myr.In the single-phase (1phase) model, all the gas is assumed to be cold, and no heating mechanisms are considered.Hence, only cold accretion is allowed and all of the feedback energy goes into driving gas outflows.In the Eddington-limited (EDDlim) model, black hole accretion cannot exceed the Eddington accretion rate.Finally, to explore the time resolution effect of the dark matter merger tree, in the 10myr model we consider a merger tree with a finer time step, 10 Myr. −2 Eddington ratio λ E -fid cold gas fraction m c -fid Eddington ratio λ E -EDDlim cold gas fraction m c -EDDlim Figure 5. Evolution of the Eddington ratio λ E = Ṁbh / ṀE (black lines) and of the cold gas fraction mc = M cold /(M h + M * + M cold + M hot + M bh ), (red lines) for the same halos as in Figure 3, showing the results for both the fid (solid lines) and EDDlim (dashed lines) scenarios.The vertical grey lines indicate all time steps at which a major merger occurred.

The importance of including gas heating and cooling
In order to assess the importance of incorporating into DEL-PHI cooling and heating processes, we first compare the time evolution of the main branches of selected halos in the fid and 1phase models in Figure 3. Starting from the final z = 4 halo, the main branch is built by selecting at each time step the most massive progenitor.With grey vertical lines, we also show the moments in which the halo undergoes major mergers, which according to our model trigger BH accretion episodes of cold gas (or prolong an ongoing accretion episode) and hence correspond to the start of faster BH growth phases.First of all, we notice that differences between the models arise only when the halo mass crosses the M crit h .In fact in the fid model only these higher-mass halos accrete hot gas, below this threshold all of the gas is cold (or assumed to cool down immediately).However, past M crit h the presence of hot gas in the fid model limits the growth of the host galaxy, and as a consequence the 1phase model can produce more stars and bigger black holes.In turn this corresponds to more gas mass ejected from the galaxy and lower gas mass available for the subsequent time step.Eventually, both black holes and stellar masses turn out to be ≈ 0.2 − 0.3 dex higher.Notice that in the 1phase scenario Bondi accretion starts as soon as the black hole is seeded, while in the fid model, where it requires the presence of hot gas, Bondi accretion occurs only for M h > M crit h .

The impact of super-Eddington accretion
In Section 3.1 we showed that distinguishing between hot and cold gas is important to properly describe the gas cycle of the galaxies.Here, we want to see how the same gas cycle can be affected by different growth models.We do this by comparing the fid model, in which BHs are allowed to grow at super-Eddington rate, and the EDDlim model, where the BH accretion rate is capped at the Eddington limit.In Figure 4 we plot the evolution of the cold and hot gas masses (M cold and M hot ) at each step, together with the amount of cold gas mass that gets heated up (by AGN feedback, M heated ) and the amount of hot gas that cools down (M cool ).In both models, as soon as the halo outgrows the critical mass (M h > M crit h ), the hot gas mass dominates the total gas budget, as the galaxy switches from cold to hotand-cold accretion mode from the IGM.A fraction of the hot gas cools down, and this fraction is lower for higher-mass halos, since from eq. 18 we can derive that the gas cooling rate follows Ṁcool ∝ M hot r cool Vvir/Rvir 2 ∝ M h 1/2 .It is also interesting to notice that while the amount of gas in outflows increases with halo mass, the total amount of heated gas mass (in fid ), result of the jet phase of the AGN life cycle, does not seem to depend on the mass.This suggests that jet feedback plays a minor role for the host galaxy evolution, and is important only during the first quasar cycle and for halos with M h ∼ 10 11.5 − 10 12.5 M⊙, where it is comparable to the amount of cold gas present in the galaxy.In the ED-Dlim model, the AGN enters its jetted phase only at a later point: this delay allows the galaxy to accumulate up to 50% more cold gas mass, explaining the relatively higher stellar and black hole masses already found at early times in Fig. 3.At the same time, once the jetted phase kicks in, given the higher BH mass and higher accretion rate, the AGN is able to heat up up to 10-15 times more gas than in the fid model, with an impact on the later evolution of the galaxy.As a consequence, EDDlim galaxies will end up containing more hot gas mass than fid galaxies, though later growth tends to smooth out all differences.
In Figure 5 we plot the redshift evolution of the Eddington ratio and the cold gas fraction for the same halos as in Figure 3 for both the fid and EDDlim models.Again, the vertical lines indicate the occurrence of major mergers.In the first scenario, shortly after crossing the M crit h threshold, all halos go through a first jetted phase with accretion rates up to λE ∼ 10 3 .Subsequent major mergers trigger weaker BH accretion episodes, with λE ∼ 0.1 − 1, consistently with observation showing that the most massive black holes at z > 4 mostly grow at sub-Eddington rates (Trakhtenbrot et al. 2017).Since the heating feedback occurs only during the jet phase, for λE ⩽ 0.01 or λE ⩾ 1, the final heated gas mass budget will be dominated by that first hyper-Eddington accretion episode at M h > ∼ M crit h , which potentially occurs in all halos.This explains why the total heated gas mass in Figure 4 (red lines) appears to be the more or less independent on the final halo mass.In the EDDlim scenario, where we do impose an Eddington cap on the accretion rate, the merger-triggered cold accretion episodes alternate with more quiet periods characterised by λE < 0.01, but by z = 4 we recover values of Eddington ratio and cold gas fraction very similar to what we find for the fid model.This is indicative of the self-regulating action of active black holes: the feedback from bigger black holes impacts the ISM more strongly, slowing down subsequent BH growth.Hence, given different accretion models, if we evolve the galaxy-BH system for long enough we expect differences to be smoothed out due to the action of AGN feedback.
As we see in Figure 6, in both scenarios our galaxies are divided between those hosting an active black hole accreting both hot and cold gas mass with log(λE) > ∼ − 1, and those hosting quieter black holes accreting only from the hot gas phase with log(λE) < ∼ − 1.5.The separation between these two populations decreases at higher galaxy masses.In fact, the average λE of active black holes shows a decreasing trend with mass, as the cold gas supply to the black hole is limited by feedback; at the same time, inactive black holes accrete at Bondi rate Ṁ hot bh ∝ M bh M h 2/3 and their average Eddington rate increases with mass.In addition, in the fid we see a third smaller population of hyper-Eddington black holes with log(λE) > ∼ 2.5 at M * ∼ 10 9−9.5 M⊙.One of the consequences of this picture, is that in the fid (i.e.super-Eddington) model we find the highest values of (see for instance Ghodla & Eldridge 2023, for a similar conclusion).In this phase, black holes are growing by up to a couple of orders of magnitudes from M bh ∼ 10 4 M⊙ to M bh ∼ 10 6 −10 6.5 M⊙ within a single time step (i.e.20 Myr), suggesting that the lifetime of intermediate-mass black holes (IMBH) in this scenario is very short, explaining why they still seem to elude detections.This appears clear from Figure 7, in which we show the evolution of the black hole mass function for both the fid and the EDDlim case.In the latter, the black hole mass functions is totally depopulated in the mass range 10 4 M⊙ < ∼ M bh < ∼ 10 6 M⊙.This result should be taken as an indication of the IMBH number density to be very low, rather than actually zero, and exactly how low depends on how smooth the passage from the stunted BH accretion rates in halos with M h < M crit h to the hyper-Eddington regime when M h > M crit h .In addition, it is easy to imagine that environmental factors, such as the ionization state of the surrounding IGM, might have an impact on the value of M crit h , which would then not be universal as it is assumed here.It is clear though that in this scenario we still expect the number density of IMBH to be much lower than that of SMBH already at z = 6.In the EDDlim scenario the IMBH number densities at z = 6 are of the order of 10 −4 Mpc −3 dex −1 , marking a clear difference with respect to fid.This difference becomes negligible at the high-mass end of the BHMF though, for M bh > ∼ 10 6 − 10 7 M⊙.Future missions targeting the IMBH mass range at high redshift like LISA will be able to put more constraints on the black hole accretion models by observing also the intermediatemass range of the BHMF.When comparing our SMBH mass function with that from other works, we see that we overpredict the observational constraint set at z = 6 by (Willott et al. 2010) by almost one order of magnitude.At the same time we fall in between the two fits proposed by Volonteri & Reines (2016), who derive the BHMF by fitting the observed local M bh − M * for different galaxy samples and applying the fit to the observed galaxy mass function.In addition we reproduce the z = 6 results from the GALFORM semi-analytic model (Griffin et al. 2020) but with the difference that we can reproduce better the high-mass end of the BHMF for M bh > 10 8 M⊙.It is worth to notice that the CAT model (Trinca et al. 2022, light green points) also reproduces a gap in the IMBH mass range, but for a different physical reason: in this case the populations from light and heavy BH seeds occupy two separated regions of the plot.

Mass correlation between galaxy and black hole
From Figures 5 and 7 we can start to see that both individual and statistical observables derived from the fid and EDDlim models tend to converge towards lower redshift (z ⩽ 6) and higher black hole masses.If we then want to observationally distinguish between the two models we cannot rely on z < 6 observations of individual sources, but we need to be able to probe a wider population at a higher redshift.With this in mind we plot, in Figure 8, the redshift evolution of the M * − M bh relation for both the fid (cyan dots) and EDDlim models (grey dots).While at z = 5 the difference between the two models for M bh > ∼ 10 6 M⊙ is negligible, due to the self-regulating BH action we have already seen in Figure 5, the two high-mass sequences are offset by almost 1 dex at z = 7, and by almost 2 dex at z = 9.At these epochs, galaxies with M * ∼ 10 9.5 M⊙ in the fid and EDDlim models respectively host black holes of masses M bh ∼ 10 6.5 M⊙ and M bh ∼ 10 8 M⊙.This is indicative of how black hole growth at earlier times proceeds much faster in the fid case, but as time goes on the black holes of the EDDlim model catch up.A faster early growth corresponds to a greater impact of AGN feedback, resulting in a lower cold gas mass (see Figure 3), and hence in a slower subsequent growth.This provides a clear way to distinguish between the two different black hole growth models at z > ∼ 7. JWST surveys have already started Figure 8. Galaxy mass-black hole mass relation for all the galaxies in our model at z = 5, 6 and 10.We compare the results from the fid (cyan dots) and EDDlim (grey dots) models.In the z=5 panel we also show z ∼ 0 observational data from (Graham 2014, G14, blue and red diamonds; direct black hole measurements for blue and red galaxies respectively) and from (Reines & Volonteri 2015, R15, stars) for quiescent elliptical galaxies (red), S/S0 with classical bulges (salmon) and pseudobulges (purple), broad-line low-luminosity AGN (violet), and black holes with reverberation mapping measurements (pink).We also show the z < ∼ 2.5 correlation derived in (Suh et al. 2020, S20) with its 1-σ scatter (grey shaded area) and the one shown in (Kormendy & Ho 2013, KH13, dashed black line) for local BH and bulge masses.The horizontal lines in the middle panel correspond to the estimated mass of the highest-z (z ∼ 7) SMBHs (Mortlock et al. 2011;Wang et al. 2021;Bañados et al. 2018, respectively M11, W21, B18).Black diamonds represent ALMA C-II-detected quasars from (Decarli et al. 2018, D18), red triangles are taken from (Izumi et al. 2019, I19), green triangles from (Maiolino et al. 2023a, M23) and the blue circles from (Kocevski et al. 2023, K23).In the third panel, together with our z = 10 results, we show the quasars reported in (Larson et al. 2023, L23, green star) and in (Maiolino et al. 2023b, M23, fuchsia pentagon) at z ∼ 9 and z ∼ 11.
to populate the M * − M bh (e.g.Larson et al. 2023;Maiolino et al. 2023b), and once the sample of AGN at z > ∼ 7 reaches a size for sufficient statistical significance in the analysis, we will be able to more firmly establish whether the typical black hole growth model allows for super-Eddington accretion or not.Looking at the morphology of the relation at M h > M crit h , we notice that both of our models show a bend in the slope of the relation at M * ∼ 10 9.5 M⊙ and z ∼ 5, due to the decreasing trend of the Eddington ratio as we move to higher stellar masses (see Figure 6).This is consistent with observational measurements in the local Universe taken from (Reines & Volonteri 2015, star symbols of different shades of red), who collected data of both quiescent elliptical galaxies and low-luminosity AGN, with the latter falling by more than one order of magnitude below the (tighter) correlations found for massive ellipticals, and characterised by a steeper slope.A mass/morphology-dependent slope is noticeable also if we look at the results from (Graham 2014, red and blue diamonds), who divided their sample into early-type and late-type galaxies.Despite the big scatter of the inferred relationship, the different slope at high and low masses seems to be a consistent feature, indicating that black holes are not able to accrete as efficiently in higher-mass galaxies.We also show the observational results from (Suh et al. 2020, light grey shaded area) who study the same relationship for a sample of low-luminosity broad-line AGN, finding no significant evolution for 0 < z < 2.5.Similar findings have been shown in other works (Salviander et al. 2015;Sun et al. 2015).In recent years there have been studies supporting the idea of an evolving normalization of the M bh − M * correlation, with an increasing black hole-tostellar mass ratio with increasing redshift (see for instance Decarli et al. 2010;Merloni et al. 2010;Trakhtenbrot et al. 2015;Caplar et al. 2018).If our z = 5 results overpredict the normalization of the relationship with respect to local measurements (left panel), they seem to agree better with higher redshift measurements (middle panel Decarli et al. 2018;Izumi et al. 2019;Maiolino et al. 2023a;Kocevski et al. 2023).However, these comparisons are to be taken carefully, since possible selection bias in the search for black hole hosts at higher redshift might skew the results towards overmassive black holes (Lauer et al. 2007), and the completeness level of high-redshift samples is not clear.In our case, the normalization of the relationship evolves (slowly) with M crit h (z): black holes of the same mass reside in progressively more massive galaxies towards lower redshift.
Despite the fact that the M * −M bh relationship can generally be used as a SMBH mass predictor in the local Universe, there are instances in which it seems to break down.For example, for red nuggets galaxies, characterised by a relatively small radius with respect to galaxies of the same mass, the M * − M bh relation seems to underestimate the black hole mass.The central velocity dispersion σc appears to better predict SMBH masses for these galaxies (e.g.Matt et al. 2023).In addition, the σc − M bh relation exhibits a smaller scatter in the local Universe (for instance Kormendy & Ho 2013), suggesting that σc is a direct tracer of the host dark matter halo mass Ferrarese (2002); Zahid et al. (2018), which in turn ultimately regulates the amount of gas that fuels both black hole and star formation activities.This motivates a closer look at the evolution of this relationship.Here, we do not directly model the value of the central velocity dispersion σc in our galaxies, so in order to estimate it we rely on the M * −σc fits provided by observational studies.Many suggest that the M * − σc relation has a mild but nonnegligible redshift evolution, at least for z > 1. Cannarozzo et al. (2020), based on an extended sample of early-type galaxies (ETGs) at 0 < z < 2.5, successfully fit the M * − σc correlation with two different models: the first one, Mevo, with an evolving slope, is parameterised as with α2 = 0.18 and β2(z) = 0.48 log(1 + z) + 2.21.We take our M * −M bh (Figure 8), we compute the velocity dispersion from our stellar masses with both fits (red dots for eq.38 and orange for eq. 3. 3), and we plot then the σc − M bh relation for our fiducial model in Figure 9.However, it is important to specify that since the fitted sample is formed by ETGs with M * > ∼ 3×10 10 M⊙, our predictions would be meaningful only for the most massive galaxies, assuming these to be the direct progenitors of ETGs at z < ∼ 2.5.We also show the theoretical predictions from the Illustris and BlueTides numerical simulations (Sijacki et al. 2015;Huang et al. 2018) and the GALFORM semi-analytical model (Malbon et al. 2007).While our results match very well with the prediction from BlueTides at z = 9, we are not able to reproduce the lower-mass BH population of Illustris, shown for z = 4.The same true for the z = 6 GALFORM results (Malbon et al. 2007), that lie slightly below ours, as they predict lower-mass black holes residing in galaxies with the same σc at higher redshift.However, our black holes grow bigger, and can reproduce well the high-z quasar population found at z > 5 (see Figure 8).For comparison we also plot the dynamically-measured black hole catalogues from the local Universe collected in Gültekin et al. (2009), Kormendy & Ho (2013) and Graham & Scott (2013).We find that at z = 5 − 9 our results of the σc − M bh relation from the Mevo model follow the same trend set by local observations, and they seem to indicate that the σc −M bh correlation is indeed independent on redshift.Given the different redshift of the sample from which the assumed M * − σ scaling, no strong conclusion can be drawn from a comparison with our results.The theoretical derivation of the σc − M bh relation from King (2003) lies extremely close to these local observational results, and is also independent on redshift.We also plot the low-mass samples presented in Barth et al. (2005) and Xiao et al. (2011) down to M bh ∼ 10 5 M⊙, also observed at z = 0. Our fid model does not produce these intermediatemass black holes, due to the fact that we consider a universal definition of M crit h .The mass range corresponding to M bh ∼ 10 5 −10 6 M⊙ could in fact be populated by black holes hosted by halos for which M crit h is lower than the universal value assumed in this work.In any case, this should not affect the conclusions drawn for SMBH.

The jet duty cycle
Given our jet model, we are able to compute the cumulative jet duty cycle (i.e.how much of their total lifetime BHs spend in the jetted phase) as a function of final black hole mass.This is defined here as the intrinsic fraction of lifetime a central black hole spends in the jetted phase since the seed is born.Observationally speaking, the duty cycle is often estimated by computing the ratio of the observed jetted sources over the total number of AGN.In Figure 10 we show the average τjet/τ bh in each BH mass bin of 0.5 dex for all the SMBH in the main branches of our merger trees in the fid model, together with its 1-σ deviation (blue points, line and shaded area).In order to assess if the outcome depends on the time resolution of the model, we compare it with the results from the 10myr model (red points, line and shaded area), built exactly like the fid model except for the fact that in this case we use a merger tree built with a time step of 10 Myr instead of 20 Myr.Despite the jet fractional lifetimes of each individual black hole can change in the two models, we notice that their statistical distributions show a remarkably precise overlap.In particular, the average jet duty cycle declines from τjet/τ bh ∼ 1 from SMBH with final mass M bh ∼ 10 6 M⊙ to τjet/τ bh ∼ 0.2 for M bh ∼ 10 11 M⊙.This is consistent with the picture we have drawn in the previous sections, according to which black holes will go through a jetted super-Eddington growth phase early in their lifetime, and then their Eddington ratio will slowly decrease as they grow in mass.

AGN number densities at z > 5
Observations of absolute AGN number densities, when compared to results from theoretical models, have to rely on obscuration models, which might be different for different observational samples and observed spectral bands.Also for this reason, high-redshift AGN surveys observe in the X-ray band, which is less subject to contamination from the host galaxy and less affected by obscuration effects than other bands.However, it is not clear how much of the total observed X-ray luminosity is coming from the central source and how much of it is contributed by the jets.In fact, if at low redshift relativistic electrons in jets mostly cool off emitting synchrotron radiation in the radio band, at high redshift the higher CMB energy density UCMB ∝ (1 + z) 4 means that these same electrons cool off through inverse-Compton scattering against CMB photons.However, we do not expect the particle density and magnetic field to be homogeneous across the jet lobes, and the high-energy regions of the lobes (hot spots) will still be dominated by synchrotron emissions.This might affect the number count of observed high-redshift jetted radio and X-ray sources (in misaligned jets the presence of hot spots might result in ad- ditional detection of radio sources) and the total X-ray luminosity (Worrall 2009;Ghisellini et al. 2013Ghisellini et al. , 2014;;Fabian et al. 2014).Given these uncertainties, we model the hard X-ray luminosity as a combination of the X-ray emission coming from black hole accretion plus a contribution from the jets L X,[2−10 keV] = Lqso/K X,[2−10keV] + αXLjet.Here, K X,[2−10keV] is the hard X-ray bolometric correction derived in Duras et al. (2020) for a wide range of luminosities and redshift.At each time step we then select all of the AGN with LX > 10 41 erg/s and we divide them up into 1-dex-wide bins of LX and add up their number densities to compute the XLF.
In Figure 11 we show our intrinsic XLF in both cases αX = 0 (orange line) and αX = 1 (black line) at z = 5 − 9.These are upper limits to the observed XLF, since we are not taking into account torus or dust obscuration models (Ueda et al. 2014;Merloni et al. 2014).Since the fraction of obscured AGN is higher at lower luminosities, our results overpredict the observed XLF especially at lower and intermediate luminosities, with a gap of ∼ 1.5 dex at L X,[2−10 keV] ∼ 10 43 − 10 44 erg/s, while they are in much better agreement at high luminosities.Keeping into account that obscuration models generally produce an obscured fraction of ∼ 90% at L X,[2−10 keV] ∼ 10 43 erg/s and of 10 − 15% at L X,[2−10 keV] > ∼ 10 45 erg/s (Ueda et al. 2014), our results still lie above the expected XLF by ∼ 0.5 mag.At z = 7, we show also the results from several other theoretical models and simulations, re-adapting a Figure taken from Amarantidis et al. (2019), where the shaded areas are defined by the modelled intrinsic and corrected XLF.In this case, our upper limits fall well within the range of predictions from the other models.Nevertheless, we generate a flatter XLF, due to the fact that in our fid model SMBH grow very fast early epochs, thanks to an early hyper-Eddington accretion phase, while most of other models implement an Eddington-limited accretion mechanism, making black hole growth more regular throughout time.In addition, in the αX = 0 model we see a dip in the XLF at L X,[2−10 keV] ∼ 10 42 − 10 43 erg/s, corresponding to a lack of intermediate-mass black holes (see Figure 7).The dip does not show up if we allow the jet power to contribute to the total X-ray luminosity (αX = 1), as in this case the black holes accreting at λE < 0.01 are also forming jets, which give an additional contribution to the X-ray luminosity, pushing the corresponding black holes towards higher luminosity bins.
In Figure 12 we apply the luminosity-dependent obscuration model from Ueda et al. (2014) by computing the unobscured fraction where ψ = 0.43 [1 + min(z, 2)] 0.48 − 0.24 (log LX − 43.75).This function saturates at 0.008 at the low-luminosity end and at 0.73 at the high-luminosity end.We then plot both the intrinsic and corrected number densities our fid model predicts at z = 5−10 for different luminosity bins and in the case αX = 0. We also show the results from Aird et al. (2015) and Ueda et al. (2014) -circles and stars -for the redshift range z = 0 − 5.The first work uses surveys with Chandra, ASCA and ROSAT to build AGN samples detected in the soft and hard X-ray bands to separately model the evolution of the the absorbed and unabsorbed AGN X-ray luminosity function to derive the evolution of the total AGN space density (dash-dotted lines).Ueda et al. (2014) uses multiple surveys (made with Swift/BAT, MAXI, ASCA, XMM-Newton, Chandra and ROSAT ) to make a population synthesis model, and shows their results for the Compton-thin AGN XLF in Figure 12.For the lowest and highest luminosity bins (10 42 erg/s < L X,[2−10 keV] < 10 43 erg/s and 10 45 erg/s < L X,[2−10 keV] < 10 46 erg/s) our z ∼ 5 results are in good agreement with the observations.For the intermediate luminosity bins, our number densities seem up to 1 − 1.5 orders of magnitude times higher than those inferred from observations.Parameters tuning can alleviate the discrepancy, but we find it hard to reproduce at the same time the number densities of intermediate and high-luminosity AGN, as we end up overestimating the former or underestimating the latter.Some cautions may be added about this point: first of all, the modelled results from Ueda et al. (2014) shown here are valid for the (intrinsic) Compton-thin population, while our intrinsic (i.e.uncorrected) X-ray luminosity function takes into account also the contribution from Compton-thick AGN.If their fraction does not exceed 42 43 44 45  (2015), taking them into account would not significantly reduce the discrepancy.However the estimates of their number density at such high redshift are still very uncertain.Secondly, the approximations used in the model are expected to affect our results.Our black hole accretion model is aimed at studying the link between the unresolved large-scale gas dynamics of the host galaxy and the time-averaged BH growth rate, but by not considering the microphysics of the BH accretion mechanism we are potentially overlooking shortterm variabilities which have an impact on our results for the AGN luminosity functions.Recent X-ray studies showed that variable sources could oscillate in luminosity by up to a factor of 10 or more, showing up in some surveys and not in others (Wolf et al. 2020;Kammoun et al. 2023), and this is likely to occur at high redshift as well.Also for this reason, it is not fully clear how complete and statistically significant the AGN samples of X-ray surveys at high redshift are.One solution could be to try to put together samples from different surveys, as done in both Ueda et al. (2014) and Aird et al. (2015).Yet, some part of the most luminous AGN population might still be going undetected, especially at high redshift, both in optical and X-ray surveys (see for instance Onken et al. 2023;Wolf et al. 2024).This might entail changes to the incompleteness corrections that have been applied to previous samples, and, consequently, changes to the observed bright end of the quasar luminosity function.To all these considerations, we should add that in the future it is worth it to test also different AGN obscuration models.Some works suggest that the AGN obscured fraction might not depend on the AGN bolometric luminosity, as considered in Ueda et al. (2014); Merloni et al. (2014), but rather on the Eddington ratio (e.g.Ricci et al. 2017).
The luminosity-dependent X-ray bolometric correction as well might carry some dependence on the Eddington ratio (see Figure 9 of Duras et al. 2020).Finally, we point out that we expect the XLF to fall off towards lower luminosities, as the stellar mass function flattens out at lower stellar masses, while the occupation fraction of X-ray luminous AGN decreases.We predict the XLF to peak at lower luminosities towards higher redshift, and specifically at LX ∼ 10 43.5 erg/s at z = 5 and LX ∼ 10 42.5 erg/s at z = 9, consistently with the evolution of M crit h .In this plot the jet contribution to the X-ray luminosity is not included (i.e.α X = 0, see text for more details).As our computation extends only to z ∼ 5, observational results for the models computed for lower z in (Aird et al. 2015, dash-dotted lines, modelled total XLF) and (Ueda et al. 2014, stars, modelled Compton-thin AGN XLF only) are shown for reference.Notice that for 10 45 erg/s < L X < 10 46 erg/s our results are shown in light green while the observation-based data in dark green for clarity.The discrepancy between our results and the observation results with assumed absorption model could be related to the uncertainties of obscuration details, as discussed in the text.

SUMMARY
In this work, building on the cosmological semianalytic model (Dayal et al. 2017;Piana et al. 2021Piana et al. , 2022)), we explore the cosmic growth of supermassive black holes, and how the evolution of their host galaxies is affected by the black hole accretion and feedback histories.We include gas cooling and heating mechanisms, hot and cold BH accretion, radiative (quasar) and jet (radio) BH feedback modes, and gas re-accretion onto the galaxy, adding them to what remains the key assumption of the model: the critical halo mass threshold M crit h below which BH growth is hindered by SN feedback (Rosas-Guevara et al. 2016;Bower et al. 2017;Lupi et al. 2019).Above this scale, the BH can accrete both from the hot and cold gas phases, respectively through con-tinuous Bondi accretion and through merger-induced accretion episodes.We aim at showing how we can observationally distinguish between different typical BH growth models and at the same time at studying the emergence of jetted AGN at z > 5.The main results and implications are summarised below.
(i) Given our assumptions, the jet-mode feedback is subdominant with respect to the radiative feedback, impacting ≈ 1% − 10 % of the total gas mass affected by AGN feedback at z > 4 (Figure 3), but it plays an important role for halos with M h ∼ 10 12 M⊙.
(ii) By comparing our fiducial super-Eddington (fid) model to the Eddington-limited (EDDlim) model, we find the the former predicts a much lower number density of black holes in the range 10 4 − 10 6 M⊙ at z > 5.This is due to the fact that BH in that mass range are extremely short-lived, resulting from a very fast growth with λE ∼ 10 3 (Figures 5 and 7).
(iii) In both models, the average Eddington ratio tends to decrease as the host galaxy mass increases and the cold gas fraction decreases (Figure 6).This means that active black holes grow faster in higher-z and lower-mass galaxies, and for this reason we observe a bend in the modelled M * − M bh (at z = 5).At the same time, the normalization of the relationship is slightly higher than the one inferred in the local Universe, consistently with what expected from the redshift evolution of M crit h .(iv) The fid model predicts BH up to two order of magnitudes more massive than the EDDlim model for the same galaxy mass at z = 9, with BH masses up to M bh ∼ 10 8.5 M⊙ (Figure 8).Moving to lower redshift, at z = 5, the two models are distinguishable only in terms of the predicted IMBH (M bh ∼ 10 4−6 M⊙) number densities, which are negligible in the super-Eddington scenario.These are clear observational predictions that will allow us to discriminate between a super-Eddington and an Eddington-limited typical BH accretion model once the on-going JWST surveys and future missions like LISA will populate the high-z BH-host galaxy correlation planes and BH mass function.
(v) Our results for the M bh − σ relationship depend on the assumed scaling between the stellar mass and the central velocity dispersion.If we apply the Mevo fit derived in Cannarozzo et al. (2020) for galaxies at 0 < ∼ z < ∼ 2.5, we find that our z ∼ 5 − 9 results are independent on redshift, and consistent with local observations (Figure 9).Given that the fitted sample was formed by ETGs, we can suppose that it should be valid for our most massive galaxies if we assume these to be the direct progenitors of the low-z ETGs.
(vi) Given our AGN jet model, we show that by z ∼ 5 SMBHs with M bh ∼ 10 6 M⊙ have spent close to 95% of their lifetime in jetted mode, while the most massive black holes only 20%, independently on the time resolution of our merger tree (Figure 10).
(vii) We have found that the current model that we have adopted either overpredicts the AGN number densities for AGNs with 10 43 erg/s < L X,[2−10 keV] < 10 45 erg/s at z = 5 (Figure 12) or underpredicts the number densities of highluminosity AGN.Nevertheless, if we compare our XLF with the results at z ∼ 7 from other theoretical models and numerical simulations, we find a relative good agreement (Figure 11).In our f id model, in which lower-mass black holes can go through strong hyper-Eddington accretion phases, we observe a generally flatter XLF than in models implementing Eddington-limited accretion rates. .Finally, it is worth touching upon a few other caveats.We are not taking into account any metallicity evolution, which directly impacts the gas cooling rates in the galaxy, and if it is true that low-mass galaxies at high redshift have been consistently observed with low metallicity, this has not been the case for high-redshift AGN (Decarli et al. 2018;Venemans et al. 2019).In particular, higher metallicities means softer spectra, which in turn means that higher stellar masses are required to reproduce the same stellar UV luminosity function.When considering black hole growth, though, AGN luminosity depends only on the amount of gas accreted, not on its metallicity.In this sense, black hole activity is self-regulated: higher accretion rates from higher cooling rates will eventually lead to more feedback and more ejected gas mass, and hence less fuel and lower accretion rates in the subsequent time step.Therefore, we do not expect the metallicity to significantly affect our resulting black hole masses.We will further study and discuss this in a future work, where we plan to implement a consistent metallicity treatment and to extend the model to z < 4. Secondly, we are neglecting stellar mass losses due to winds, which would return part of the stellar mass to the gas phase and inject more energy into the ISM.However, the energy injection is dominated by high-mass stars, but is still subdominant with respect to SN feedback.At the same time, the gas mass returned to the ISM through winds is mostly contributed by lower-mass stars, which do not evolve much during the first Gyr of life of the Universe.In addition, our transitional halo mass M crit h is unlikely to be a universal value, and should depend on the environment in which galaxies sit, and on the reionization state of the Universe.When taken into account, we would expect some BH to start growing in lower-mass galaxies and some in higher-mass galaxies, hence the IMBH gap visible in the BH mass function in the fid case would partially fill up.In addition, Ocvirk et al. (2008) showed that at high redshift the fraction of gas f cold accreted in cold phase by the galaxy is not a step function of the halo mass, as we assume here, rather a smooth decreasing function.However, when we take this effect into account we notice no significant change in our results.Indeed, the self-regulating action of BH growth and feedback in high-mass halos, as shown in Figure 5, smooths out the difference between the two implementations.We plan to address these issues in future works, as we expand the model down to z = 0.
Figure3.Mass evolution for different galaxy components for 1phase model (no gas heating and cooling, solid lines), compared with the fid model (dashed lines), here used as a reference.In each panel we follow the growth of the main branch of the merger tree of halos of different mass.We show the total gas mass content (M gas tot ) at each time step, the cumulative black hole and stellar mass (M bh and M * ), the cumulative gas mass ejected by AGN and stellar feedback (M ej bh and M ej * ).The vertical grey lines indicate all the time steps at which a major merger occurred.

Figure 4 .
Figure4.Mass evolution for different galaxy components for EDDlim model (solid lines), compared with the fid model (dashed lines), here used as a reference.In each panel we follow the growth of the main branch of the merger tree of the same halos as in Fig.3.We show the cold and hot gas mass content (M cold and M hot ), and the cumulative heated (M heated ) and cooled (M cool ) gas masses.Notice that the last two quantities and the hot gas mass content are defined only for the fid case.The vertical grey lines indicate all the time steps at which a major merger occurred.

Figure 6 .
Figure6.Galaxy mass-Eddington ratio relation at different redshift.We compare the results from the fid (cyan dots) and EDDlim (grey dots) models.Notice that even for the highest values of the intrinsic Eddington ratio Ṁbh / ṀE , the observed ratio will reach at most values of Lqso/L E ∼ 4 (eq.33).

Figure 7 .
Figure7.Total black hole mass function for the fid and the ED-Dlim scenarios (solid and dashed blue lines respectively) at z = 6.We compare our results with those from(Volonteri & Reines 2016), who transpose the observed galaxy mass function to the black hole mass function by applying different fits of the M bh −M * relationship, derived from different samples: local galaxies fromMarconi & Hunt (2003) andHäring & Rix (2004) for the vanilla fit (van, dashed-dotted cyan line), low-luminosity AGN in lowermass galaxies for the low-mass fit (lm, long-dashed red line).We also show the observational results from(Willott et al. 2010, orange shaded area), as well as results from the GALFORM(Griffin et al. 2020, dashed dark green line) and CAT(Trinca et al. 2022, grey diamonds and black pentagons respectively to indicate descendants of light and heavy BH seeds, born respectively with M light ∼ 100 M ⊙ and M heavy = 10 5 M ⊙ ) semi-analytic models.

Figure 9 .
Figure9.Central velocity dispersion-back hole mass relation for all the galaxies in our fiducial model at different redshift.σc is computed by applying the two fits inCannarozzo et al. (2020), one assuming an evolving slop (Mevo) and the other one assuming a redshift-independent slope Mconst.We also plot the local dynamically-measured black hole catalogues collected in(Gültekin et al. 2009,  Gu09, cyan and green stars respectively for ellipticals and spirals) and from (Graham 2014, Gr14, cyan and green diamonds for ellipticals and spirals), which partially overlap.The (local) low-mass black hole measurements are form(Barth et al. 2005, B09)  and(Xiao et al.  2011, X11).The maroon shaded area represents the z = 4 black hole population from the Illustris simulationSijacki et al. (2015) while the grey one is the z = 9 prediction from the BlueTides simulationsHuang et al. (2018).The light green solid line is the z = 6 result from the Galform model(Malbon et al. 2007), while the dotted line is results from the theoretical argument exposed in(King 2003,  K03).Finally, as a benchmark, we plot the (local) fit from (Kormendy & Ho 2013, KH13, black line).As we have seen in Figure7, our model does not produce intermediate mass black holes at these redshift.At lower redshift, the gap would be populated by black holes accreting at Bondi rate.

Figure 10 .
Figure10.Fraction of lifetime BHs spend in the jetted phase as a function of black hole mass for M bh > 10 6 M ⊙ at z = 4.We compare the results of our fiducial (blue line) to those of the 10myr (red line) scenario, in which the time step is τs = 10 Myr, as opposed to τs = 20 Myr for the fiducial case.The line are the averages of the scatter points computed in bins of 0.5 dex in black hole mass, and the shaded areas represent the 1-σ deviations from the average.

Figure 11 .
Figure11.Intrinsic AGN X-ray luminosity function for the α X = 0 and α X = 1 cases.Our results can be considered as upper limits of the observed XLF, as we know that a fraction of AGN are obscured by the central torus, which absorbs part of the light emitted by the black holeUeda et al. (2014);Merloni et al. (2014).In the left panel we show our z = 5 results together with observations of the hard XLF at z = 4 − 5 taken fromRichards et al. (2006),Fiore et al. (2012) andUeda et al. (2014), at z = 4 − 5.2 fromFontanot et al. (2007), at z = 3.8 − 5.2 fromGlikman et al. (2011), and at z = 3.5 − 5 fromAird et al. (2015).The hatched regions define the region of space forbidden by our model.In the middle panel we compare our z = 7 results with those from other theoretical models and simulations at 7 < z < 8, taken from Figure3ofAmarantidis et al. (2019).In the right panel we show our results at z = 9.The shaded areas are defined by the corrected and the intrinsic XLF.

Figure 12 .
Figure12.Evolution of the total AGN intrinsic (solid lines) and obscuration-corrected (dashed lines) number densities for the fid model, as a function of hard X-ray luminosity and redshift.In this plot the jet contribution to the X-ray luminosity is not included (i.e.α X = 0, see text for more details).As our computation extends only to z ∼ 5, observational results for the models computed for lower z in(Aird et al. 2015, dash-dotted lines, modelled total XLF) and(Ueda et al. 2014, stars, modelled   Compton-thin AGN XLF only) are shown for reference.Notice that for 10 45 erg/s < L X < 10 46 erg/s our results are shown in light green while the observation-based data in dark green for clarity.The discrepancy between our results and the observation results with assumed absorption model could be related to the uncertainties of obscuration details, as discussed in the text.

Table 1 .
Model parameters, default values and defining equation.

Table 2 .
List of models.