The mass assembly of high-redshift black holes

We use the Delphi semi-analytic model to study the mass assembly and properties of high-redshift ($z>4$) black holes over a wide mass range, $10^3<M_{bh}/M_\odot<10^{10}$. Our black hole growth implementation includes a critical halo mass ($M_{h}^{crit}$) below which the black hole is starved and above which it is allowed to grow either at the Eddington limit or proportionally to the gas content of the galaxy. As a consequence, after an initial growth phase dominated by black hole mergers down to $z \sim 7 (9)$, supermassive black holes in $z=4$ halo masses of $M_h|_{z=4} \sim 10^{11.75} ~ (10^{13.4}) M_\odot$ mainly grow by gas accretion from the interstellar medium. In particular, we find that: (i) while most of the accretion occurs in the major branch for $M_h|_{z=4} \sim 10^{11-12} M_\odot$ halos, accretion in secondary branches plays a significant role in assembling the black hole mass in higher-mass halos ($M_h|_{z=4} \gtrsim 10^{12} M_\odot$); (ii) while the Eddington ratio increases with decreasing redshift for low-mass ($M_{bh}<10^5 M_\odot$) black holes, it shows the opposite trend for larger masses. In addition, since the accretion rate depends on the gas mass present in the host halo, the duty cycle of the Eddington-limited accretion phase -- which can last up to $\approx 650$ Myr -- is crucially linked to the joint assembly history of the black hole and its host halo.


INTRODUCTION
High-redshift observations have discovered that about 200 quasars, powered by accretion onto supermassive black holes (SMBH), were already in place in the first billion years of the Universe (e.g. Fan et al. 2001;Mortlock et al. 2011;Mazzucchelli et al. 2017;Bañados et al. 2018;Matsuoka et al. 2018). The masses of these black holes can be as high as M bh ∼ 1.2 × 10 10 M at a redshift of z = 6.3 (Wu et al. 2015) and M bh ∼ 7.8 × 10 8 M at z = 7.54 (Bañados et al. 2018). The discoveries of such early SMBHs have re-ignited the -still largely unanswered -question of how they could have grown so large so fast (for a recent review see Inayoshi et al. 2019). Growing a black hole of 10 9 M from a seed of about 100M (at a reasonable radiative efficiency of 10%) would require near-Eddington accretion for almost 800 Myr -this corresponds to the age of the universe at z ∼ 6.8 (Haiman & Loeb 2001;Tanaka 2014). The prob-piana@astro.rug.nl lem is exacerbated by the fact that high-velocity recoil kicks (of up to 500 − 1000kms −1 ) can significantly reduce the efficiency of black hole growth through mergers (Fitchett 1983;Favata et al. 2004;Haiman 2004; by expelling them from the shallow potential wells of low-mass halos. While some studies have tried to alleviate this problem by allowing black holes to undergo short episodes of super-Eddington accretion within the context of a slim accretion disk model (Haiman 2004;Yoo & Miralda-Escudé 2004;Volonteri & Rees 2005;Lupi et al. 2016;Pezzulli et al. 2016;Mayer 2019), others have suggested that pristine gas clouds could collapse isothermally and without fragmentation to yield direct-collapse black hole (DCBH) seeds with mass M bh ∼ 10 4 − 10 5 M (Haehnelt & Rees 1993;Loeb & Rasio 1994;Eisenstein & Loeb 1995;Bromm & Loeb 2003;Begelman et al. 2006;Habouzit et al. 2016;Luo et al. 2018). Alternatively, Boco et al. (2020) recently proposed that heavy black hole seeds of masses 10 4 −10 6 M could form via merging of compact stellar remnants in young galaxies with dense gas environments.
A huge amount of theoretical effort has been dedicated to modelling black hole evolution at high redshift using both semi-analytic models (SAMs; e.g. Volonteri et al. 2008;Croton et al. 2006;Bower et al. 2006;Hirschmann et al. 2012;Barausse 2012;Sesana et al. 2014;Mutch et al. 2016;Lacey et al. 2016;Valiante et al. 2016;Marshall et al. 2019) and hydrodynamic simulations (e.g. Schaye et al. 2015;Khandai et al. 2015;Di Matteo et al. 2017;Dubois et al. 2014;Vogelsberger et al. 2014;Sijacki et al. 2015;Latif et al. 2018;Habouzit et al. 2019;Huang et al. 2019). While SAMs offer speed and flexibility in exploring a larger parameter space and uncovering the properties of the black hole population as a whole, simulations offer the advantage of tracking the dynamics of the different galactic components. Nevertheless, zoom-in simulations are required both to reach the parsec-scale resolution needed to accurately model black hole accretion during the early growth phase (Lupi et al. 2019), when the internal structure of the host gas cloud plays a major role (Beckmann et al. 2019), and to follow the (sub-) kiloparsec-scale gas inflows from the cosmic web which drive black hole growth at later evolutionary stages (Costa et al. 2014;Richardson et al. 2016); these techniques, though, do not allow a statistical study of the black hole population as a whole.
In general these works are tuned to reproduce the lowredshift black hole-host galaxy data. Episodes of fast black hole growth are triggered by galaxy mergers, disk instabilities or bursts of star formation, with the black hole accretion rate often assumed to be equal to the Bondi-Hoyle rate and capped at one or two times the Eddington limit. Since the efficiency of black hole growth through mergers, which are much more common at high redshift, is hampered by the black hole merging timescale (?) and by recoil kicks (Haiman 2004), most works agree that SMBHs at z > 6 mainly grow by gas accretion (see Wise 2019, for a review). Despite this progress in understanding the evolution of the high-redshift black hole population, many questions still remain open.
In this paper we want to study how the black hole physical properties and growth histories vary as a function of host halo mass across the entire population. To this end, we use the semi-analytic model Delphi, that has been shown to reproduce all key observables for galaxies and active galactic nuclei (AGN) at z > ∼ 5 (including the AGN and stellar ultra-violet luminosity function, the stellar mass function, the stellar mass density and the black hole mass function, Dayal et al. 2019) in addition to reionization observables (the electron scattering optical depth and emissivity constraints, Dayal et al. 2020). The key strengths of this model lie in that: (i) it is seeded with two types of black hole seeds (stellar and direct collapse); (ii) the black hole accretion rate is linked to the host halo mass; and (iii) it uses a minimal set of free parameters to model star formation, black hole accretion and their associated feedback. In particular, we make use of the results from numerical simulations (Dubois et al. 2015;Bower et al. 2017;Habouzit et al. 2017) finding that black hole growth is suppressed in low-mass halos, due to the fact that SN feedback in low-mass galaxies hinders the accumulation of gas mass around the central regions. With this picture in mind, we assume that black holes can not accrete below a critical halo mass threshold M crit h . So far, Delphi is the only SAM explicitly implementing such a mechanism to regulate BH growth.
We start by describing the theoretical model in Sec. 2. We then discuss the assembly histories of black holes and their host halos in Sec. 3 and present the resulting key properties of black holes (occupation fraction, duty cycle and the mass function) in Sec. 4 before concluding in Sec. 5.

THEORETICAL MODEL
In this work, we use the semi-analytic code Delphi (Dark matter and the emergence of galaxies in the epoch of reionization; Dayal et al. 2014Dayal et al. , 2019, which jointly tracks the assembly of the dark matter, baryonic and black hole components of high-redshift (z > ∼ 4) galaxies. In brief, starting at z = 4, we build analytic merger trees for 550 halos, equally separated in log space in the mass range M h = 10 8 − 10 13.5 M , up to z = 20 in time steps of 20 Myr. Each z = 4 halo is assigned a number density according to the Sheth-Tormen halo mass function (HMF, Sheth & Tormen 2002) at z = 4. We then assign the same number density also to all the progenitors of the halo; we have confirmed that the resulting HMFs are in accord with the Sheth-Tormen HMF at all redshifts for z ∼ 5 − 20.
The very first progenitors of each tree are seeded with an initial gas mass proportional to the halo mass such that Mgi = (Ω b /Ωm)M h . Such first halos irradiated by a Lyman-Werner (LW) background of intensity (JLW ) larger than a critical value Jcrit = αJ21 (where J21 = 10 −21 ergs −1 Hz −1 cm −2 sr −1 and α is a free parameter) are assigned a DCBH seed with a mass between 10 3 − 10 4 M . First halos (at z > 13) not fulfilling this criteria are instead assigned a stellar black hole (SBH) seed with a mass of 150 M . Starting from these first progenitors and going forward in time, at each redshift step the halos and their baryonic components can grow both through the smooth accretion of dark matter and gas from the intergalactic medium (IGM) and through mergers; similarly, black holes grow through accretion from the interstellar medium (ISM) and mergers. At a given time, the star formation efficiency is defined as the minimum between the type-II Supernovae (SNII) energy required to unbind the rest of the gas and an upper threshold f * . In the interest of simplicity, each newly formed stellar population follows a Salpeter initial mass function (Salpeter 1955)  the black hole is stuck in a stunted accretion mode. Once the halo outgrows M crit h , at each time step its black hole is allowed to accrete either a fixed fraction f ac bh of the gas mass present in the galaxy or a mass corresponding to what it would accrete if it were growing at the Eddington rate, whichever lower. The physical reasoning behind M crit h is the following: while gas inflows towards the central black hole are disrupted by SN feedback in low-mass galaxies (see for instance Dubois et al. 2015;Habouzit et al. 2017;Bower et al. 2017;Anglés-Alcázar et al. 2017), Lupi et al. (2019) have shown that the black hole mass growth speeds up at later times, once the host galaxy reaches a stellar mass of M * ≈ 10 10 M , and SN feedback is not effective in keeping the gas away from the central regions anymore. In our work, the transitional host mass M crit h (z) has been tuned to simultaneously reproduce the number densities of low-luminosity AGN at z 5 and the black hole mass function at z = 6 (see Dayal et al. 2020).
The mass accreted by the black hole at any redshift step can then be written as where M Edd (z) =Ṁ Edd × ∆t is the mass accreted by the black hole in a time step assuming the Eddington growth rate, f ac bh = 5.5 × 10 −4 represents the maximum fraction of the total gas mass left in the host galaxy that can be accreted onto the black hole, r = 0.1 is the radiative efficiency of the black hole and f Edd is defined as A fixed fraction f w bh = 3 × 10 −3 of the total energy emitted by the accreting black hole is allowed to couple to the gas content of the host halo.
While we assume the mergers between dark matter halos to occur instantaneously, we explore two prescriptions for the mergers of their baryonic components: instantaneous, and merging after a delay induced by dynamical friction (Lacey & Cole 1993). Finally, we include the impact of reionization feedback. Creating an ultra-violet (UV) heating background, reionization can photo-evaporate the gas mass from low-mass halos (Dayal et al. 2019). We implement this feedback in its maximal form by suppressing the gas content of all halos with a virial velocity Vvir(z) < 40 kms −1 .
To test the robustness of our approximations, we present results for two cases: (i) the instantaneous model (ins1) corresponds to galaxies and black holes merging instantaneously with their host halos, no UV feedback and DCBH seeds using a value of α = 30; (ii) the delayed model (tdf 4) includes the prescription for the dynamical friction delay, a higher value of α = 300 for DCBH seed formation (which effectively reduces the number of DCBH seeds by a factor of ≈ 50) and the impact of the reionization feedback. These models and their free parameter values used are summarised in In Table 1.
We now briefly discuss the interplay between the baryonic processes implemented in our model and black hole growth. Gas accretion onto the host galaxy (either from the IGM or through galaxy mergers), star formation, SN feedback and AGN feedback all contribute in regulating the amount of gas mass in the galaxies, thereby having a direct effect on black hole growth. The extent of their impact changes across the different growth phases of the black hole and its host. Increasing the star formation rate threshold f * or the coupling factor f w * means that SN feedback are more effective in ejecting gas, especially in low-mass galaxies. These parameters are tuned to match the galaxy UV luminosity function, the stellar mass function, the stellar mass density and the star formation rate density. As the host halo grows and the potential well deepens, though, SN feedback becomes less effective. The value of f Edd (M h < M crit h ) (the fraction of the Eddington rate at which the black hole is allowed to accrete), is tuned to match the faint end of the AGN UV luminosity function. Once M h > M crit h and the Eddington-limited accretion phase has started, the impact of AGN feedback grows, along with the black hole mass, with an exponential trend. We tested the values of the parameters strictly related to black hole growth (M crit h , f ac bh and f w bh ) by looking at how changing them would affect the AGN luminosity function and the black hole mass function. An increase in M crit h leads to a decrease of the amplitude of the black hole mass function at M bh > ∼ 10 4 M , to a lower normalisation of the AGN luminosity function and to a slightly higher number density of very luminous AGN. In fact a higher M crit h means a higher gas mass at the onset of the Eddington-limited accretion phase, which can then be sustained for a longer period of time. The parameter f ac bh directly intervenes in Eq. 1: a higher f ac bh -or a lower coupling factor f w bh -translates into a later switch of the growth mode to sub-Eddington accretion rates, yielding a longer duty cycle. Once the black hole has exited the Eddington-limited phase, an increase in f ac bh or in the coupling factor f w bh would mean less gas mass available for the subsequent time step, compensating the change in f ac bh and f w bh . In this regime black hole growth is self-regulated. These parameters are in turn tuned to reproduce simultaneously both the AGN UV LF and the black hole mass function.

THE MASS ASSEMBLY OF EARLY BLACK HOLES AND THEIR HOST HALOS
In this section, we start by discussing the evolution of the black hole accretion rates (expressed in terms of the Eddington fraction) as a function of the halo mass at different redshifts in Sec. 3.1. We then discuss the merger-and accretion-driven assembly of the black hole and dark matter components of early galaxies in Section 3.2.

Black hole accretion rates
The relation between black holes and their host galaxies remains poorly constrained at high redshifts. Given that the only confirmed black holes at high redshifts are those powering luminous quasars, where the light from the accreting black hole over-shines that from the host galaxy, the stellar mass/bulge mass and the stellar velocity of the host cannot be measured. The best estimates of the dynamical mass for quasars are then obtained through measurements of the cold molecular gas properties in sub-millimetre observations (e.g. Venemans et al. 2016;Shao et al. 2017;Decarli et al. 2018, and references therein). Black hole masses instead are generally inferred using rest-frame ultra-violet indicators (e.g. Vestergaard 2002;McLure & Jarvis 2002).
In Fig. 1 we show the black hole mass-halo mass relation at z = 5, 7 and 9 for both the ins1 and the tdf 4 scenarios, color-coded by the black hole accretion channel. We start by discussing the results of the ins1 model, which yields the upper limit of the black hole mass. Independently of redshift, black holes residing in low-mass halos (M h < ∼ 10 9 M ) do not accrete at all, since SNII feedback ejects all of the gas mass from the host galaxies. The double tail at M h ∼ 10 9−11 M is due to the presence of two types of black hole seeds: the upper part is comprised of (mixed) black holes resulting from the mergers of DCBHs and SBHs, while the lower part is comprised of SBHs. The black holes in this halo mass range show stunted growth, given that they accrete at a rateṀ bh = 7.5 × 10 −5Ṁ Edd . Interestingly, since DCBH seeds are assigned solely on the basis of halo bias (see discussion in Dayal et al. 2019), they cover almost the entire halo mass range at all the redshifts studied. Once the critical halo mass threshold is reached, black holes can start accreting at the Eddington rate (Eddingtonlimited phase); this phase ends when the gas-limited condi- This happens for the highest-mass halos with M h 10 12 M at z = 5. Our black hole mass-halo mass relation is therefore described by a triple power-law at z ∼ 5: On the other hand, notice that black holes at z > ∼ 7 never enter the last accretion phase. Finally, the slight increase of M crit h with decreasing redshift naturally results in more massive halos (and hence black holes) accreting at the Eddington rate with decreasing redshift.
We then look at results from the tdf 4 scenario that includes the impact of the UV background and a dynamical friction-induced time delay in mergers. As seen from the same figure, the M bh − M h relation in this scenario follows a similar slope as the ins1 model for M h > ∼ 10 11 M at all the redshifts considered. However, as a result of the higher LW threshold used, the total number of DCBH black hole seeds in the tdf 4 model is lower by a factor of ≈ 50. In addition, this model has a larger scatter (by about a factor of two) in the black hole masses at a given halo mass compared to the ins1 model. This is in accord with previous works (e.g. Peng 2007) that find the number of merging events to anti-correlate with the width of the scatter in the black hole mass-host mass relation.
We find that for halos with M h > ∼ 10 11 M , for both models, while the slope of our M bh − M h relation is consistent with a number of other (observational and theoretical) works, its normalisation is bracketed by the observational/theoretical ranges. This is quite encouraging given the different proxies used to estimate the halo masses in low-z observational works, ranging from bulge properties (Ferrarese 2002) to gravitational lens modelling (Bandara et al. 2009), as well as theoretical model results extrapolating local relations to z = 5 (Croton 2009).
From an observational perspective the stellar mass is easier to measure than the halo mass, especially since telescopes like JWST will be able to estimate stellar masses also at high redshifts. Therefore, as a complement to Fig. 1, we show in Fig. 2 the M * − M bh relationship as predicted from our model. A debate has been going on in the literature regarding the evolution of this relation (see the introduction of Suh et al. 2020, for a review on the current state of the discussion). In brief, some works inferred a positive evolution with redshift of the M bh − M * ratio, suggesting that at earlier times black hole growth was more efficient than that of the host galaxies (see for instance Peng et al. 2006;Merloni et al. 2010;Decarli et al. 2010;Caplar et al. 2018  2019) studies. To complicate things further, selection biases at high redshift might affect these measurements by favouring the detection of overmassive active black holes (Lauer et al. 2007). From Fig. 2 it appears that the M bh − M * relation of the low-luminosity subset of our AGN sample, corresponding to M * ∼ 10 8.5−9.5 M , does not evolve significantly with redshift, and is compatible with the results from Suh et al. (2020). On the other hand, galaxies with masses M * > ∼ 10 10 M sit ∼ 1 dex higher with respect to the same low-redshift measurements and the relation for moderateluminosity local AGN shown in Reines & Volonteri (2015), pointing towards a positive evolution of the high-mass end of the M bh − M * relationship. This dual behaviour of the relationship is also found in other studies (e.g. Merloni et al. 2010), but we point out that any comparison with observations at z < 4 and our results cannot yield any solid conclusion, as our model does not extend to lower redshift, and any extrapolation of higher-redshift results would not take into account the different physical processes intervening in more recent epochs, such as Helium reionisation and galaxy quenching processes. For the same mass range, our relation also seems to be in very good agreement with the relation inferred for AGN in old high-mass elliptical galaxies (Reines & Volonteri 2015). In conclusion, our results point to black holes hosted by massive galaxies above 10 9.5 M to grow efficiently and be more massive than black holes in galaxies with the same stellar mass at lower redshift, while black holes hosted in low-mass galaxies below 10 9.5 M are less massive and follow a steeper relation (Volonteri & Stark 2011). The different behaviour at low and high stellar mass is simply a reflection of the behaviour at low and high halo mass.
With this picture in mind, we can now look at the evolution of the Eddington ratio, defined as λ Edd =Ṁ ac bh /Ṁ Edd , which is indicative of the black hole growth efficiency. By construction, we have that λ Edd = 7.5×10 −5 for halo masses M h < M crit h ∼ 10 11.5 M , while λ Edd is generally between 0.1 and 1 for black holes in higher-mass halos. The results for the average Eddington ratios (λ Edd ) as a function of redshift and for different black hole mass bins are shown in Fig. 3. Low-mass black holes (M bh ∼ 10 3−4 M ) reside in progressively more massive halos with decreasing z, eventually starting to accrete at the Eddington rate (also see Fig. 1), leading to an increases of λ Edd from 10 −3 at z = 9 to 0.1 by z = 4. Moving on, approximately 90% of the black holes in the mass range 10 4 − 10 5 M are in the Eddingtonlimited regime at z = 4, with λ Edd ≈ 0.9, while the rest of the black holes sits in low-mass halos and accretes at λ Edd = 7.5 × 10 −5 . Black holes with masses 10 5−6 M generally accrete at λ Edd = 1 at all redshifts, as also seen from Fig. 1. As black holes grow above 10 6 M and switch to the gas-limited accretion phase, the average Eddington ratio drops from λ Edd = 1 to progressively lower values ∼ 0.1 (see e.g. Merloni 2004). This is because as M bh increases, so does the Eddington mass, with the consequence that more massive black holes switch to the gas-limited phase at higher redshifts. Note that by construction our black hole accretion rate is capped at the Eddington rate, so we do expect a tension between our results and the observed Eddington rato distribution at high redshift, which extends to super-Eddington regimes (see for instance Kelly & Shen 2013). In this sense, the results showed in this plot are meant to portrait a qualitative more than quantitative trend. Compared to the ins1 model, in the tdf 4 scenario, the largest difference is seen for the lowest mass bins (M bh = 10 3−5 M ). We also see a slight increase in λ for M bh 10 6 M . This is because a delay in the baryonic mass assembly in the tdf 4 scenario causes (i) the black holes to enter the Eddington-limited phase with a lower mass, and (ii) to stay in this phase longer, since the condition M Edd (z) < M gf * (z)f ac bh (1 − r ) is met for a longer amount of time. Additional delays in black holes mergers which are not taken into account here, such as the black hole binary inspiraling phase and the final parsec timescale, might further increase the average Eddington ratio.

Joint assembly of black holes and their host halos
We now study the relative contribution of mergers and accretion to the joint assembly of black holes and their host halos at z > 4. For both black holes and halos, the merged mass is that brought in by all the progenitors of the previous z-step; the accreted mass is the ISM and IGM mass that is accreted onto the black hole and halo, respectively. Here we explore results for 21 galaxies over three halo mass bins (at z = 4) each: low-mass halos (M h ∼ 10 9.9−10.1 M ), intermediate mass halos (M h ∼ 10 11.65−11.85 M ) and high-mass halos (M h ∼ 10 13.3−13.5 M ). We start by discussing results from the ins1 scenario. Naturally, both the ins1 and tdf 4 scenarios show the same results for halo assembly; the dynamical friction delay only affects the baryonic component.
Low-mass halos (panel a of Fig. 4) start assembling their dark matter mass at z ∼ 14. Most of the halo growth is driven by accretion down to z ∼ 4, at which point mergers start contributing to the halo mass assembly in equal measure. On the other hand, their black hole assembly is dominated by stochastic mergers across the whole redshift range z ∼ 13 − 4. The accretion mode is sub-dominant given the negligible accretion rates (∼ 7.5 × 10 −5 M edd ), and contributes less than 10% of the total black hole mass by z ∼ 4.
Intermediate-mass halos (panel b of the same figure) naturally start assembling earlier, at z ∼ 17.5. Accretion dominates the mass assembly until z ∼ 7, below which point mergers start contributing equally. While mergers dominate the black hole mass build-up down to z ∼ 7, thereafter Eddington-limited accretion takes over, dominating the black hole mass budget by ∼ 2 orders of magnitude by z = 4. By then, more than 95% of this black hole accretion has taken place in the major branch (see Fig. 5), wherein black holes accrete in the Eddington-limited phase. The high variance here is due to the fact that different halos, with their different assembly histories, overcome the critical halo mass for Eddington-limited accretion at different redshifts.
The dark matter mass assembly of high-mass halos is dominated by accretion until z ∼ 6, below which mergers take over as the main growth mechanism (see also Cattaneo 2002). As for black holes, mergers dominate (and increase the mass by about three orders of magnitude) between z ∼ 17.5 − 9; accretion dominates instead the later phases, contributing about 90% of the final mass at z ∼ 4. It is interesting also to notice that for high-mass halos the variance in the merged components is generally lower, and significant only at z 7, when stochastic major mergers take place. Further, more than one progenitor is able to overcome M crit h , and on average only about 20 (75%) of the mass in assembled in the major branch at z ∼ 10 (4) as shown in Fig. 5.
In the tdf 4 case the black hole mass assembly proceeds slower, due to fewer mergers (because of dynamical friction) that bring in lower gas masses (because of the photoevaporation of gas from low-mass halos).
To summarise, while black holes residing in low-mass halos can grow only through (rare) mergers, those in intermediate-and high-mass mass halos predominantly grow   through mergers only in the early phases, before switching to an accretion-dominated growth phase. As expected, mergers are less important in the tdf 4 scenario, contributing ∼ 10 −3.5 % (10 −2 %) of the total mass in intermediate (high) mass halos. Finally, in Fig. 6 we plot the redshift evolution of the ratio of the cumulative black hole mass to the cumulative halo mass for the same halo mass bins as in Fig. 4. As shown, black holes in low-mass halos (when present) are starved throughout the redshift range z = 14 − 4, driving the mass ratio to a constant decrease. Intermediate-mass halos show a flat evolution at z > 13, when the black hole seed formation efficiency is still 100% (i.e. every newly-born halo has a black hole). At z ∼ 13 − 6, the ratio drops as new halos form without black hole seeds. At z < ∼ 6, the ratio shows a steep rise as the black holes residing in the biggest progenitors enter the Eddington-limited accretion phase. This results in a rise from M bh /M h ∼ 10 −7.4 at z = 6 to M bh /M h ∼ 10 −4.5 at z = 4. While the evolution of the highest mass halos is qualitatively similar to intermediate mass halos, their black holes already enter the Eddington-limited growth channel at z ≈ 8. Sustaining this growth for a longer time, they show a mass ratio that increases from M bh /M h ∼ 10 −6.8 at z = 8 to M bh /M h ∼ 10 −3.6 at z = 4. Again, intermediate mass halos show a larger variance with respect to high-mass halos as a result of the diverse merging and accretion histories of their black holes.

KEY PROPERTIES OF THE BLACK HOLE POPULATION IN EARLY GALAXIES
We now explore the key properties of BHs in early galaxies, specifically focusing on their occupation fraction and duty cycles in this section.

The black hole occupation fraction
We define the black hole occupation fraction (f bh ) as the fraction of halos hosting a central black hole. In the case of merging systems, this refers to the fraction of central galaxies hosting a black hole. We note that this specification is important only in the tdf 4 scenario; in the ins1 case galaxies merge at the same time as their halos. In our model, halos are populated with SBH seeds down to z = 13. Lower-z halos contain a black hole if they are either seeded with a DCBH or if they gain one through mergers. The latter (and  Figure 7. The redshift evolution of the minimum halo mass for which the central black hole occupation fraction f bh = 1 (thick lines) and f bh = 0.1 (thin lines) for both the ins1 (solid lines) and the tdf 4 (dashed lines) scenarios, as marked.
dominant) mechanism naturally implies that the minimum halo mass for which f bh = 1 increases with decreasing redshift. Indeed, as shown in Fig. 7, we find that f bh = 0.1 (i.e. 10% of halos contain a central black hole) for a minimum halo mass of M h ∼ 10 8.4 M at z ∼ 12, which increases to M h ∼ 10 9.5 M by z = 5. Reaching f bh = 1 naturally requires a higher halo mass of M h ∼ 10 8.9 (10 10.5 )M at z ∼ 12 (5). As seen from the same plot, there is no sensible difference in the ins1 and tdf 4 results for f bh = 0.1. However, the minimum halo mass for f bh = 1 is about 0.3 − 0.4 dex higher in the tdf 4 model as compared to ins1 model. This is because a merging galaxy and its black hole reach the central galaxy at a later time step as compared to the instantaneous (ins1) scenario.
We then now deconstruct f bh as a function of halo mass for different black hole mass bins, and show the results in Fig. 8. Firstly, as expected, independent of the redshift and model used, there is a clear trend of black holes of increasing mass being hosted in increasingly massive halos. In terms of the most massive black holes, while by z 9 the most massive halos (M h ∼ 10 11.5 M ) host black holes as massive as ∼ 10 6 M , this increases to M bh ∼ 10 9.5 M by z = 5 for halos with M h ∼ 10 13.4 M . At each redshift the overlap between the different occupation fraction curves is indicative of the intrinsic scatter of the M bh −M h relation shown in Fig. 1. Indeed, while at z = 9 halos with M h ∼ 10 11−12 M host black holes with M bh ∼ 10 4 −10 6 M , at z = 5 the same halo masses host black holes spanning four orders of magnitude in mass, between 10 3.5 − 10 7.5 M . This is mainly due to the black holes residing in these halos entering the Eddingtonlimited accretion regime. Further, deviations from the average M bh − M h relation, which depend on the black hole mass at the time the halo reaches M crit h , are amplified by 2 logM bh < 3 3 logM bh < 4 4 logM bh < 5 5 logM bh < 6 6 logM bh < 7 7 logM bh < 8 8 logM bh < 9 9 logM bh < 10 all f bh Figure 8. The central black hole occupation fraction (f bh ) as a function of halo mass for the different black hole mass bins noted. The occupation fractions are shown using lines for the ins1 case and as shaded areas for the tdf4 case. The solid thin cyan line shows the cumulative black hole fraction for the ins1 case.
the exponential growth of black holes during the Eddingtonlimited growth phase.
In the tdf 4 case, halos of a given mass range typically host black holes of slightly lower masses compared to ins1. This is due to a combination of delayed mergers and (to a lower extent) lower gas masses, given that mergers are the dominant mode of black hole growth at z > ∼ 9 in intermediate-and high-mass halos. As accretion overtakes the mass build-up at lower z, results from the ins1 and tdf 4 models come into closer agreement.

Time-weighed duty cycle
We now explore the black hole duty cycle defined as the fractional lifetime a black hole spent in the stunted, Eddingtonlimited and gas-limited accretion modes. This is done for all the z = 4 halos that host a black hole, accounting for all of the progenitors along the merger tree. The timeweighed duty cycle for each accretion mode is then defined as τi = ti/ttot, where ti indicates the time spent in each of the three accretion channels (i) and ttot is the total lifetime of the black hole.
As expected in a hierarchical growth scenario, most of the lifetime of the progenitors of the final black hole is spent in small halos, where they are subject to stunted accretion. In particular, black hole progenitors of halos with M h |z=4 < M crit h spend all of their lifetime in a stunted accretion mode or -if they are in SN-feedback dominated halos (which are more common among the progenitors of highmass halos) -not accreting at all. For M h |z=4 10 11.5 M halos, the time spent in both the Eddington and gas-limited regimes decreases with increasing halo mass, again due to an increasing fraction of low-mass halos in their merger trees. The black holes in intermediate-mass halos spend at most 5% and 1% of their lifetime accreting in the Eddington-and gas-limited modes, respectively; these drop to 0.1% and 1% for the highest mass halos.
The most important difference upon implementing dynamical friction delay and UV feedback is that the black hole lifetime spent in non-accreting halos increases by at least a factor of 3 for the highest-mass halos (and up to an order of magnitude for the lowest-mass halos) due to UV suppression of their gas mass.
Finally, the variations (up to one order of magnitude) between the duty cycles of halos of very similar final masses is driven by the assembly histories of both black holes and their host halos, which determine both the halo and black hole masses, as well as the amount of gas mass available for accretion onto the black hole.

Black hole predictions at z 5
We conclude our analysis by proposing a couple of testable predictions of our model. In Fig. 10 we show the resulting X-ray luminosity function (XLF) at z = 5, 7 and 8. Since the model provides us with the bolometric luminosity information of each AGN, we build the luminosity function just by binning our AGN sample. We adopt the Duras et al. (2020) X-ray bolometric correction, while the obscured fraction as a function of luminosity is taken from Ueda et al. (2014). We compare our results to the z = 5 observational XLF from Ueda et al. (2014) (left panel), and we also compare the 7 < z < 8 results of other SAMs and simulations taken from Fig. 3 of Amarantidis et al. (2019) to our XLF at z = 7 (central panel) and z = 8 (right panel). We use our intrinsic and obscured X-ray luminosity functions as upper and lower limit on the XLF, since for the other results considered here the upper and lower limits correspond to the estimations for the maximum and minimum obscuration effects. We point out that our lower limit represents the median obscured XLF obtained from 1000 realisations of the obscuration mechanism: for each realisation, in each luminosity bin, we randomly draw a subsample of the AGN populating that bin corresponding to the obscured fraction and manually set their AGN luminosity to zero. The red error bars represent the 16th and 84th percentiles of the resulting set of obscured XLFs. We can easily notice in the z=7 and z=8 panels that introducing the threshold halo mass M crit h for black hole growth delays the build-up of our XLF with respect to that of the other models, but this delay disappears at z=5 when comparing it to the observational results.
Together with the luminosity function, the mass function is the other main predictable statistical observable, though as of now it is still very loosely constrained at high redshift (Willott et al. 2010). In Fig. 11 we show our z = 5−7 black hole mass function. Here we focus on the results from the ins1 model and deconstruct it to show the contributions from all three black hole accretion channels. The general behaviour follows that expected from the hierarchical structure formation model: black holes with M bh < ∼ 10 7 M follow a power law with their number densities being virtually constant at z ∼ 5−7. On the other hand, for M bh 10 7 M , the number density of high-mass black holes show an exponential cut-off, specially at z ∼ 5 − 6, in addition to showing a rapid evolution between z ∼ 6−7. In particular, the mass bin 10 6 M M bh 10 7 M experiences the greatest variation of number density in the shortest timescale, increasing by almost three orders of magnitude between z = 7 and z = 6 (i.e. in less than 300Myr). As a consequence of the different evolution in number densities of the low-and high-mass ends of the BHMF, at z = 5 black holes with M bh ∼ 10 4.5 M and M bh ∼ 10 7 M have very similar number densities.  one in a very short time span, while the halo mass grows much slower. Therefore, the shallow slope of the BHMF in the intermediate mass range 10 4 M < M bh < 10 7 M at z = 5 is a direct effect of the critical halo mass M crit h imposed in our model.
As noted before, we reiterate that the black hole accretion channel strongly depends on its mass: while lowmass black holes with M bh ∼ 10 4 M accrete at very low (7.5 × 10 −5 ) Eddington rates, intermediate-mass black holes can accrete at the Eddington rate while the most massive black holes show a variable Eddington ratio.
In the same figure, we also compare our BHMF to other theoretical and observational works. At z = 6, our results are in qualitative agreement with the observational BHMF built (using spectroscopically-selected quasars) by Willott et al. (2010), as well as the theoretical results of Volonteri & Reines (2016), based on extrapolating the observed z = 0 black hole-host galaxy correlations to higher redshifts. At z = 5, our results are in reasonable agreement with those from the EAGLE simulation (Rosas-Guevara et al. 2016). We are also in qualitative agreement with the results from the Horizon-AGN project  and the AGN synthesis model by Merloni & Heinz (2008), although these models seem to hint to a characteristic knee mass lower by ≈ 0.5 dex. Finally at z = 5 − 6, our results are within the range bracketed by the GALFORM (Griffin et al. 2019) semi-analytic model (higher than our BHMF by more than 1 dex at z = 6 for 10 6 M < ∼ M bh < ∼ 10 8 M ) and the Dragons semi-analytic model (Qin et al. 2017), which predicts a lower number density (by about 0.5 dex for 10 6 M < ∼ M bh < ∼ 10 8 M at z = 5). This mismatch is indicative of the different black hole seeding and growth mechanisms implementations. Griffin et al. (2019) assume that black hole growth proceeds proportionally to the star formation rate in the host galaxy, and that accretion episodes can occur also during galaxy mergers able to drive gas towards the central regions. Given the absence of any halo mass threshold for the growth of black holes, they start assembling at earlier times with respect to our model. We then expect their BHMF to evolve faster than ours at higher redshift. In Qin et al. (2017), instead, black holes are allowed to accrete a (variable) fraction of the hot gas present in the galaxy, while cold gas can only be accreted only during mergers; there is no such distinction between cold and hot gas accretion in our model, which might explain why we produce more massive black holes.
The fact that we are able to reproduce very well the observed black hole mass function at z = 5, but at the same time finding a slight mismatch between our X-ray luminosity function and the one inferred at z = 5 by Ueda et al. (2014) at a luminosity of LX ∼ 10 44 erg s −1 , might point to a tension between these results. Nevertheless, the differ-works (e.g. Ricarte & Natarajan 2018) that conclude that a hybrid scenario with a bimodal black hole seed mass distribution (such as ours) yields basically the same results as a light-seed one, given the low number density of DCBHs compared to SBHs.
To conclude, the reader should keep in mind a few caveats. First, we are not accounting for gravitational recoil kicks in mergers between massive black holes, which could either offset the central (merged) black hole or eject it altogether, hampering its subsequent growth (Blecha et al. 2016). Nevertheless, spin alignment between the merging black holes might reduce the effectiveness of this mechanism, especially if the merging system is embedded in a circumbinary gaseous disk (Dotti et al. 2010). We also point out that the black hole merging timescales are still probably severely underestimated: not only we do not model the black hole binary inspiralling phase after the galaxies have merged, but we are also not accounting for the final parsec problem, which requires stellar scattering processes for allowing the black holes to coalesce together. In addition, we do not consider that in this low-mass regime the stellar and gaseous components can actually drive the dynamics of the black holes, introducing stochasticity in their orbits and possibly affecting the black hole merger rate, as shown by Pfister et al. (2019). This last mechanism might widely stretch or even prevent black holes from merging at all. We should also mention that the intensity of the Lyman-Werner background necessary to form DCBH seeds is still very poorly constrained in the literature, and values of α higher (lower) than the ones we have chosen in Table 1 would result in a lower (higher) number density of DCBH seeds forming at 8 < z < 13. Finally, we manage to form black holes as massive as 10 10 M at z = 5 by allowing them to accrete a fraction of the total gas mass present in the host galaxy, rather than only the gas confined in the central regions. For these reasons the black hole masses that our model ins1 yields should be considered as an upper limit.
We will dedicate a future work to study how the black hole growth portrayed here affects the evolution of the stellar mass of the host galaxy.