<figure><img src="image_1.png" /></figure>

HAL
open science

# The mass assembly of high-redshift black holes

Olmo Piana, Pratika Dayal, Marta Volonteri, Tirthankar Roy Choudhury

► To cite this version:

Olmo Piana, Pratika Dayal, Marta Volonteri, Tirthankar Roy Choudhury. The mass assembly of high-redshift black holes. Monthly Notices of the Royal Astronomical Society, 2021, 500 (2), pp.2146-2158. 10.1093/mnras/staa3363. hal-02973143

HAL Id: hal-02973143

https://hal.science/hal-02973143

Submitted on 28 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The mass assembly of high-redshift black holes

Olmo Piana, Pratika Dayal, Marta Volonteri, Tirthankar Roy Choudhury

► To cite this version:

Olmo Piana, Pratika Dayal, Marta Volonteri, Tirthankar Roy Choudhury. The mass assembly of high-redshift black holes. Monthly Notices of the Royal Astronomical Society, 2020, 500 (2), pp.2146-2158. 10.1093/mnras/staa3363. hal-02973143

HAL Id: hal-02973143

https://hal.science/hal-02973143v1

Submitted on 28 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# The mass assembly of high-redshift black holes

Olmo Piana<sup>1,★</sup> Pratika Dayal<sup>1,†</sup> Marta Volonteri<sup>2</sup> and Tirthankar Roy Choudhury<sup>3</sup>

Accepted 2020 October 22. Received 2020 October 12; in original form 2020 April 2

## ABSTRACT

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_{\text{crit}}^h$) 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_{\text{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_{\text{h}|z=4} \sim 10^{11-12} M_{\odot}$ haloes, accretion in secondary branches plays a significant role in assembling the black hole mass in higher mass haloes ($M_{\text{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.

**Key words:** galaxies: evolution – galaxies: haloes – galaxies: high-redshift – quasars: supermassive black holes – galaxies: statistics.

## 1 INTRODUCTION

High-redshift observations have discovered that about 200 quasars, powered by accretion on to supermassive black holes (SMBHs), 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} \sim 1.2 \times 10^{10} M_{\odot}$ at a redshift of $z = 6.3$ (Wu et al. 2015) and $M_{bh} \sim 7.8 \times 10^8 M_{\odot}$ at $z = 7.54$ (Bañados et al. 2018). The discoveries of such early SMBHs have reignited the – still largely unanswered – question of how they could have grown so large so fast (for a recent review see Inayoshi, Visbal & Haiman 2020). Growing a black hole of $10^9 M_{\odot}$ from a seed of about $100 M_{\odot}$ (at a reasonable radiative efficiency of 10 per cent) would require near-Eddington accretion for almost 800 Myr – this corresponds to the age of the universe at $z \sim 6.8$ (Haiman & Loeb 2001; Tanaka 2014). The problem is exacerbated by the fact that high-velocity recoil kicks (of up to 500–1000 km s$^{-1}$) can significantly reduce the efficiency of black hole growth through mergers (Fitchett 1983; Favata, Hughes & Holz 2004; Haiman 2004; Volonteri & Rees 2006) by expelling them from the shallow potential wells of low-mass haloes. 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 disc model (Haiman 2004; Yoo & Miralda-Escudé 2004; Volonteri & Rees 2005; Lupi et al. 2016; Pezzulli, Valiante & Schneider 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} \sim 10^4-10^5 M_{\odot}$ (Haehnelt & Rees 1993; Loeb & Rasio 1994; Eisenstein & Loeb 1995; Bromm & Loeb 2003; Begelman, Volonteri & Rees 2006; Habouzit et al. 2016; Luo et al. 2018). Alternatively, Boco, Lapi & Danese (2020) recently proposed that heavy black hole seeds of masses $10^4-10^6 M_{\odot}$ 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. Bower et al. 2006; Croton et al. 2006; Volonteri, Lodato & Natarajan 2008; Barausse 2012; Hirschmann et al. 2012; Sesana et al. 2014; Lacey et al. 2016; Mutch et al. 2016; Valiante et al. 2016; Marshall et al. 2020) and hydrodynamic simulations (e.g. Dubois et al. 2014; Vogelsberger et al. 2014; Khandai et al. 2015; Schaye et al. 2015; Sijacki et al. 2015; Volonteri et al. 2016; Di Matteo et al. 2017; Latif, Volonteri & Wise 2018; Habouzit et al. 2019; Huang et al. 2020). 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, Devriendt & Slyz 2019), and to follow the (sub-)kiloparsec-scale gas inflows from the cosmic web that drive black hole growth at later evolutionary stages

<footnote n="1"><sup>1</sup>Kapteyn Astronomical Institute, University of Groningen, PO Box 800, NL-9700 AV Groningen, the Netherlands</footnote>

<footnote n="2"><sup>2</sup>Institut d'Astrophysique de Paris, Sorbonne Universites, UPMC Univ Paris 6 et CNRS, UMR 7095, 98 bis bd Arago, F-75014 Paris, France</footnote>

<footnote n="3"><sup>3</sup>National Centre for Radio Astrophysics, Tata Institute of Fundamental Research, Pune 411007, India</footnote>

<footnote>* E-mail: p.o.piana@rug.nl</footnote>
<figure><img src="image_2.png" /></figure>

# The mass assembly of high-redshift black holes

Olmo Piana<sup>1,★</sup> Pratika Dayal<sup>1,†</sup> Marta Volonteri<sup>2</sup> and Tirthankar Roy Choudhury<sup>3</sup>

Accepted 2020 October 22. Received 2020 October 12; in original form 2020 April 2

## ABSTRACT

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_{\text{crit}}^h$) 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_{\text{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_{\text{h}|z=4} \sim 10^{11-12} M_{\odot}$ haloes, accretion in secondary branches plays a significant role in assembling the black hole mass in higher mass haloes ($M_{\text{h}|z=4} \gtrsim 10^{12} M_{\odot}$); (ii) while the Eddington ratio increases with decreasing redshift for low-mass ($M_{\text{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.

**Key words:** galaxies: evolution – galaxies: haloes – galaxies: high-redshift – quasars: supermassive black holes – galaxies: statistics.

## 1 INTRODUCTION

High-redshift observations have discovered that about 200 quasars, powered by accretion on to supermassive black holes (SMBHs), 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_{\text{bh}} \sim 1.2 \times 10^{10} M_{\odot}$ at a redshift of $z = 6.3$ (Wu et al. 2015) and $M_{\text{bh}} \sim 7.8 \times 10^8 M_{\odot}$ at $z = 7.54$ (Bañados et al. 2018). The discoveries of such early SMBHs have reignited the – still largely unanswered – question of how they could have grown so large so fast (for a recent review see Inayoshi, Visbal & Haiman 2020). Growing a black hole of $10^9 M_{\odot}$ from a seed of about $100 M_{\odot}$ (at a reasonable radiative efficiency of 10 per cent) would require near-Eddington accretion for almost 800 Myr – this corresponds to the age of the universe at $z \sim 6.8$ (Haiman & Loeb 2001; Tanaka 2014). The problem is exacerbated by the fact that high-velocity recoil kicks (of up to 500–1000 km s$^{-1}$) can significantly reduce the efficiency of black hole growth through mergers (Fitchett 1983; Favata, Hughes & Holz 2004; Haiman 2004; Volonteri & Rees 2006) by expelling them from the shallow potential wells of low-mass haloes. 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 disc model (Haiman 2004; Yoo & Miralda-Escudé 2004; Volonteri & Rees 2005; Lupi et al. 2016; Pezzulli, Valiante & Schneider 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_{\text{bh}} \sim 10^4-10^5 M_{\odot}$ (Haehnelt & Rees 1993; Loeb & Rasio 1994; Eisenstein & Loeb 1995; Bromm & Loeb 2003; Begelman, Volonteri & Rees 2006; Habouzit et al. 2016; Luo et al. 2018). Alternatively, Boco, Lapi & Danese (2020) recently proposed that heavy black hole seeds of masses $10^4-10^6 M_{\odot}$ 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. Bower et al. 2006; Croton et al. 2006; Volonteri, Lodato & Natarajan 2008; Barausse 2012; Hirschmann et al. 2012; Sesana et al. 2014; Lacey et al. 2016; Mutch et al. 2016; Valiante et al. 2016; Marshall et al. 2020) and hydrodynamic simulations (e.g. Dubois et al. 2014; Vogelsberger et al. 2014; Khandai et al. 2015; Schaye et al. 2015; Sijacki et al. 2015; Volonteri et al. 2016; Di Matteo et al. 2017; Latif, Volonteri & Wise 2018; Habouzit et al. 2019; Huang et al. 2020). 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, Devriendt & Slyz 2019), and to follow the (sub-)kiloparsec-scale gas inflows from the cosmic web that drive black hole growth at later evolutionary stages

<footnote n="1"><sup>1</sup>Kapteyn Astronomical Institute, University of Groningen, PO Box 800, NL-9700 AV Groningen, the Netherlands</footnote>

<footnote n="2"><sup>2</sup>Institut d'Astrophysique de Paris, Sorbonne Universites, UPMC Univ Paris 6 et CNRS, UMR 7095, 98 bis bd Arago, F-75014 Paris, France</footnote>

<footnote n="3"><sup>3</sup>National Centre for Radio Astrophysics, Tata Institute of Fundamental Research, Pune 411007, India</footnote>

<footnote>* E-mail: p.o.piana@rug.nl</footnote>
(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 low-redshift black hole-host galaxy data. Episodes of fast black hole growth are triggered by galaxy mergers, disc 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 time-scale (Tremmel et al. 2018) 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 SAM Delphi that has been shown to reproduce all key observables for galaxies and active galactic nuclei (AGN) at $z \gtrsim 5$ [including the AGN and stellar ultraviolet (UV) luminosity function, the stellar mass function, the stellar mass density, and the black hole mass function (BHMF); 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, Volonteri & Dubois 2017) finding that black hole growth is suppressed in low-mass haloes, due to the fact that supernova (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 cannot accrete below a critical halo mass threshold $M_{\text{h}}^{\text{crit}}$. So far, Delphi is the only SAM explicitly implementing such a mechanism to regulate black hole growth.

The cosmological parameters used in this work are $\Omega_m$, $\Omega_\Lambda$, $\Omega_b$, $h$, $n_s$, $\sigma_8 = 0.3089$, 0.6911, 0.049, 0.67, 0.96, 0.81, respectively (Planck Collaboration XIII 2016). We quote all quantities in comoving units, unless stated otherwise, and express all magnitudes in the standard AB system (Oke & Gunn 1983).

We start by describing the theoretical model in Section 2. We then discuss the assembly histories of black holes and their host haloes in Section 3 and present the resulting key properties of black holes (occupation fraction, duty cycle, and the mass function) in Section 4 before concluding in Section 5.

## 2 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. 2014, 2019), which jointly tracks the assembly of the dark matter, baryonic and black hole components of high-redshift ($z \gtrsim 4$) galaxies. In brief, starting at $z = 4$, we build analytic merger trees for 550 haloes, equally separated in log space in the mass range $M_h = 10^8-10^{13.5} M_\odot$, 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 \sim 5-20$.

The very first progenitors of each tree are seeded with an initial gas mass proportional to the halo mass such that $M_{\text{gi}} = (\Omega_b/\Omega_m)M_h$. Such first haloes irradiated by a Lyman-Werner (LW) background of intensity ($J_{\text{LW}}$) larger than a critical value $J_{\text{crit}} = \alpha J_{21}$ (where $J_{21} = 10^{-21}$ erg s$^{-1}$ Hz$^{-1}$ cm$^{-2}$ sr$^{-1}$ and $\alpha$ is a free parameter) are assigned a DCBH seed with a mass between $10^3$ and $10^4 M_\odot$. First haloes (at $z > 13$) not fulfilling this criterion are instead assigned a stellar black hole (SBH) seed with a mass of $150 M_\odot$. Starting from these first progenitors and going forward in time, at each redshift step the haloes 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 (SNeII) 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) over 0.1–100 $M_\odot$, has a metallicity $Z = 0.05$ $Z_\odot$ and an age of 2 Myr. A fraction of the gas mass $M_{*}^{\text{gf}}(z)$ left in the halo after the star formation and SNeII feedback can be accreted on to the black hole. The black hole accretion rate depends on the mass and the redshift of the host halo through a critical halo mass $M_{\text{h}}^{\text{crit}}(z) = 10^{11.25}[\Omega_m(1+z)^3 + \Omega_\lambda]^{-0.125}$, below which the black hole is stuck in a stunted accretion mode. Once the halo outgrows $M_{\text{h}}^{\text{crit}}$, at each time step its black hole is allowed to accrete either a fixed fraction $f_{\text{bh}}^{\text{ac}}$ 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_{\text{h}}^{\text{crit}}$ 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; Anglés-Alcázar et al. 2017; Bower et al. 2017; Habouzit 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_{*} \approx 10^{10} M_\odot$, and SN feedback is not effective in keeping the gas away from the central regions anymore. In our work, the transitional host mass $M_{\text{h}}^{\text{crit}}(z)$ has been tuned to simultaneously reproduce the number densities of low-luminosity AGN at $z \ge 5$ and the BHMF at $z = 6$ (see Dayal et al. 2020).

The mass accreted by the black hole at any redshift step can then be written as

$$M_{\text{bh}}^{\text{ac}}(z) = \min[f_{\text{Edd}}M_{\text{Edd}}(z), (1 - \epsilon_r)f_{\text{bh}}^{\text{ac}}M_{*}^{\text{gf}}(z)], \quad (1)$$

where $M_{\text{Edd}}(z) = \dot{M}_{\text{Edd}} \Delta t$ is the mass accreted by the black hole in a time step assuming the Eddington growth rate, $f_{\text{bh}}^{\text{ac}} = 5.5 \times 10^{-4}$ represents the maximum fraction of the total gas mass left in the host galaxy that can be accreted on to the black hole, $\epsilon_r = 0.1$ is the radiative efficiency of the black hole, and $f_{\text{Edd}}$ is defined as

$$f_{\text{Edd}} = \begin{cases} 7.5 \times 10^{-5} & M_h(z) < M_h^{\text{crit}}(z), \\ 1 & M_h(z) \ge M_h^{\text{crit}}(z). \end{cases} \quad (2)$$

A fixed fraction $f_{\text{bh}}^w = 3 \times 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 haloes 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 UV heating background, reionization can photoevaporate the gas mass from low-mass haloes (Dayal et al. 2019). We implement this feedback in its
(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 low-redshift black hole-host galaxy data. Episodes of fast black hole growth are triggered by galaxy mergers, disc 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 time-scale (Tremmel et al. 2018) 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 SAM Delphi that has been shown to reproduce all key observables for galaxies and active galactic nuclei (AGN) at $z \gtrsim 5$ [including the AGN and stellar ultraviolet (UV) luminosity function, the stellar mass function, the stellar mass density, and the black hole mass function (BHMF); 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, Volonteri & Dubois 2017) finding that black hole growth is suppressed in low-mass haloes, due to the fact that supernova (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 cannot accrete below a critical halo mass threshold $M_{\text{crit}}^h$. So far, Delphi is the only SAM explicitly implementing such a mechanism to regulate black hole growth.

The cosmological parameters used in this work are $\Omega_m$, $\Omega_\Lambda$, $\Omega_b$, $h$, $n_s$, $\sigma_8 = 0.3089$, 0.6911, 0.049, 0.67, 0.96, 0.81, respectively (Planck Collaboration XIII 2016). We quote all quantities in comoving units, unless stated otherwise, and express all magnitudes in the standard AB system (Oke & Gunn 1983).

We start by describing the theoretical model in Section 2. We then discuss the assembly histories of black holes and their host haloes in Section 3 and present the resulting key properties of black holes (occupation fraction, duty cycle, and the mass function) in Section 4 before concluding in Section 5.

## 2 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. 2014, 2019), which jointly tracks the assembly of the dark matter, baryonic and black hole components of high-redshift ($z \gtrsim 4$) galaxies. In brief, starting at $z = 4$, we build analytic merger trees for 550 haloes, equally separated in log space in the mass range $M_h = 10^8-10^{13.5} M_\odot$, 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 \sim 5-20$.

The very first progenitors of each tree are seeded with an initial gas mass proportional to the halo mass such that $M_{\text{gi}} = (\Omega_b/\Omega_m)M_h$. Such first haloes irradiated by a Lyman-Werner (LW) background of intensity ($J_{\text{LW}}$) larger than a critical value $J_{\text{crit}} = \alpha J_{21}$ (where $J_{21} = 10^{-21}$ erg s$^{-1}$ Hz$^{-1}$ cm$^{-2}$ sr$^{-1}$ and $\alpha$ is a free parameter) are assigned a DCBH seed with a mass between $10^3$ and $10^4 M_\odot$. First haloes (at $z > 13$) not fulfilling this criterion are instead assigned a stellar black hole (SBH) seed with a mass of $150 M_\odot$. Starting from these first progenitors and going forward in time, at each redshift step the haloes 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 (SNeII) 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) over 0.1–100 $M_\odot$, has a metallicity $Z = 0.05$ $Z_\odot$ and an age of 2 Myr. A fraction of the gas mass $M_{*}^{\text{gf}}(z)$ left in the halo after the star formation and SNeII feedback can be accreted on to the black hole. The black hole accretion rate depends on the mass and the redshift of the host halo through a critical halo mass $M_{\text{crit}}^h(z) = 10^{11.25}[\Omega_m(1+z)^3 + \Omega_\Lambda]^{-0.125}$, below which the black hole is stuck in a stunted accretion mode. Once the halo outgrows $M_{\text{crit}}^h$, at each time step its black hole is allowed to accrete either a fixed fraction $f_{\text{bh}}^{\text{ac}}$ 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_{\text{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; Anglés-Alcázar et al. 2017; Bower et al. 2017; Habouzit 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_* \approx 10^{10} M_\odot$, and SN feedback is not effective in keeping the gas away from the central regions anymore. In our work, the transitional host mass $M_{\text{crit}}^h(z)$ has been tuned to simultaneously reproduce the number densities of low-luminosity AGN at $z \ge 5$ and the BHMF at $z = 6$ (see Dayal et al. 2020).

The mass accreted by the black hole at any redshift step can then be written as

$$M_{\text{bh}}^{\text{ac}}(z) = \min[f_{\text{Edd}}M_{\text{Edd}}(z), (1 - \epsilon_r)f_{\text{bh}}^{\text{ac}}M_{*}^{\text{gf}}(z)], \quad (1)$$

where $M_{\text{Edd}}(z) = \dot{M}_{\text{Edd}} \Delta t$ is the mass accreted by the black hole in a time step assuming the Eddington growth rate, $f_{\text{bh}}^{\text{ac}} = 5.5 \times 10^{-4}$ represents the maximum fraction of the total gas mass left in the host galaxy that can be accreted on to the black hole, $\epsilon_r = 0.1$ is the radiative efficiency of the black hole, and $f_{\text{Edd}}$ is defined as

$$f_{\text{Edd}} = \begin{cases} 7.5 \times 10^{-5} & M_h(z) < M_h^{\text{crit}}(z), \\ 1 & M_h(z) \ge M_h^{\text{crit}}(z). \end{cases} \quad (2)$$

A fixed fraction $f_{\text{bh}}^w = 3 \times 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 haloes 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 UV heating background, reionization can photoevaporate the gas mass from low-mass haloes (Dayal et al. 2019). We implement this feedback in its
<table><caption>**Table 1.** Here we list the free parameters of the model in column 1, we describe them in column 2, and we show their values in the *ins1* and the *tdf4* scenarios in columns 3 and 4, respectively.</caption><thead><tr><th>Parameter</th><th>Description</th><th>ins1</th><th>tdf4</th></tr></thead><tbody><tr><td>ε<sub>r</sub></td><td>Radiative efficiency of black hole accretion</td><td>0.1</td><td>0.1</td></tr><tr><td>f<sub>*</sub></td><td>Star formation efficiency threshold</td><td>0.02</td><td>0.02</td></tr><tr><td>f<sub>*</sub><sup>w</sup></td><td>Fraction of SN energy that couples to the gas</td><td>0.1</td><td>0.1</td></tr><tr><td>f<sub>bh</sub><sup>w</sup></td><td>Fraction of AGN energy that couples to the gas</td><td>0.003</td><td>0.003</td></tr><tr><td>f<sub>bh</sub><sup>ac</sup></td><td>Fraction of available gas mass that black holes can accrete</td><td>5.5 × 10<sup>-4</sup></td><td>5.5 × 10<sup>-4</sup></td></tr><tr><td>f<sub>Edd</sub>(M<sub>h</sub> &lt; M<sub>h</sub><sup>crit</sup>)</td><td>Black hole accretion rate in fraction of Eddington</td><td>7.5 × 10<sup>-5</sup></td><td>7.5 × 10<sup>-5</sup></td></tr><tr><td>f<sub>Edd</sub>(M<sub>h</sub> &gt; M<sub>h</sub><sup>crit</sup>)</td><td>Black hole accretion rate in fraction of Eddington</td><td>1</td><td>1</td></tr><tr><td>α</td><td>LW background threshold for DCBH formation (in units of J<sub>21</sub>)</td><td>30</td><td>300</td></tr><tr><td>Reionization feedback</td><td>–</td><td>No</td><td>Yes</td></tr><tr><td>Delayed mergers</td><td>Dynamical friction acting to delay the merging of the baryonic components</td><td>No</td><td>Yes</td></tr></tbody></table>

maximal form by suppressing the gas content of all haloes with a virial velocity $V_{vir}(z) < 40$ km s$^{-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 haloes, no UV feedback and DCBH seeds using a value of $\alpha = 30$; (ii) the delayed model (tdf4) includes the prescription for the dynamical friction delay, a higher value of $\alpha = 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 summarized in Table 1.

We now briefly discuss the interplay between the baryonic processes implemented in our model and black hole growth. Gas accretion on to 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_\ast^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_h^{crit})$ (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_h^{crit}$ 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_h^{crit}$, $f_{bh}^{ac}$, and $f_{bh}^{w}$) by looking at how changing them would affect the AGN luminosity function and the BHMF. An increase in $M_h^{crit}$ leads to a decrease of the amplitude of the BHMF at $M_{bh} \gtrsim 10^4 M_\odot$, to a lower normalization of the AGN luminosity function and to a slightly higher number density of very luminous AGN. In fact a higher $M_h^{crit}$ 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_{bh}^{ac}$ directly intervenes in equation (1): a higher $f_{bh}^{ac}$ – or a lower coupling factor $f_{bh}^{w}$ – 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_{bh}^{ac}$ or in the coupling factor $f_{bh}^{w}$ would mean less gas mass available for the subsequent time step, compensating the change in $f_{bh}^{ac}$ and $f_{bh}^{w}$. In this regime black hole growth is self-regulated. These parameters are in turn tuned to reproduce simultaneously both the AGN UV luminosity function and the BHMF.

# 3 THE MASS ASSEMBLY OF EARLY BLACK HOLES AND THEIR HOST HALOES

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 Section 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.

## 3.1 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 overshines 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 submillimetre 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 UV indicators (e.g. McLure & Jarvis 2002; Vestergaard 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 tdf4 scenarios, colour 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 haloes ($M_h \lesssim 10^9 M_\odot$) do not accrete at all, since SNeII feedback ejects all of the gas mass from the host galaxies. The double tail at $M_h \sim 10^{9-11} M_\odot$ is due to the presence of two types of black hole seeds: the upper part is composed of (mixed) black holes resulting from the mergers of DCBHs and SBHs, while the lower part is composed of SBHs. The black holes in this halo mass range show *stunted growth*, given that they accrete at a rate $\dot{M}_{bh} = 7.5 \times 10^{-5} \dot{M}_{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 (Eddington-limited phase); this phase ends when the gas-limited condition $M_*^{gf}(z) f_{bh}^{ac} (1 - \epsilon_r) < M_{Edd}(z)$ is met. This happens for the highest mass haloes with $M_h \gtrsim 10^{12} M_\odot$ at $z = 5$. Our black hole mass–halo mass relation is therefore described by a triple power law at $z \sim 5$: flat at $M_h \lesssim M_h^{crit}$, almost vertical at $M_h^{crit} \lesssim M_h \lesssim 10^{11.6} M_\odot$, and $\log(M_{bh}/M_\odot) \propto \log(M_h/M_\odot)^{1.4}$ at $M_h \gtrsim 10^{11.6} M_\odot$. On the other hand, notice that black holes at $z \gtrsim 7$ never enter
<table><caption>**Table 1.** Here we list the free parameters of the model in column 1, we describe them in column 2, and we show their values in the *ins1* and the *tdf4* scenarios in columns 3 and 4, respectively.</caption><thead><tr><th>Parameter</th><th>Description</th><th>ins1</th><th>tdf4</th></tr></thead><tbody><tr><td>ε<sub>r</sub></td><td>Radiative efficiency of black hole accretion</td><td>0.1</td><td>0.1</td></tr><tr><td>f<sub>*</sub></td><td>Star formation efficiency threshold</td><td>0.02</td><td>0.02</td></tr><tr><td>f<sup>w</sup><sub>*</sub></td><td>Fraction of SN energy that couples to the gas</td><td>0.1</td><td>0.1</td></tr><tr><td>f<sup>w</sup><sub>bh</sub></td><td>Fraction of AGN energy that couples to the gas</td><td>0.003</td><td>0.003</td></tr><tr><td>f<sup>ac</sup><sub>bh</sub></td><td>Fraction of available gas mass that black holes can accrete</td><td>5.5 × 10<sup>-4</sup></td><td>5.5 × 10<sup>-4</sup></td></tr><tr><td>f<sub>Edd</sub>(M<sub>h</sub> &lt; M<sub>h</sub><sup>crit</sup>)</td><td>Black hole accretion rate in fraction of Eddington</td><td>7.5 × 10<sup>-5</sup></td><td>7.5 × 10<sup>-5</sup></td></tr><tr><td>f<sub>Edd</sub>(M<sub>h</sub> &gt; M<sub>h</sub><sup>crit</sup>)</td><td>Black hole accretion rate in fraction of Eddington</td><td>1</td><td>1</td></tr><tr><td>α</td><td>LW background threshold for DCBH formation (in units of J<sub>21</sub>)</td><td>30</td><td>300</td></tr><tr><td>Reionization feedback</td><td></td><td>No</td><td>Yes</td></tr><tr><td>Delayed mergers</td><td>Dynamical friction acting to delay the merging of the baryonic components</td><td>No</td><td>Yes</td></tr></tbody></table>

maximal form by suppressing the gas content of all haloes with a virial velocity $V_{vir}(z) < 40$ km s$^{-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 haloes, no UV feedback and DCBH seeds using a value of $\alpha = 30$; (ii) the delayed model (tdf4) includes the prescription for the dynamical friction delay, a higher value of $\alpha = 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 summarized in Table 1.

We now briefly discuss the interplay between the baryonic processes implemented in our model and black hole growth. Gas accretion on to 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_h^{crit})$ (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_h^{crit}$ 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_h^{crit}$, $f_{bh}^{ac}$, and $f_{bh}^w$) by looking at how changing them would affect the AGN luminosity function and the BHMF. An increase in $M_h^{crit}$ leads to a decrease of the amplitude of the BHMF at $M_{bh} \gtrsim 10^4 M_\odot$, to a lower normalization of the AGN luminosity function and to a slightly higher number density of very luminous AGN. In fact a higher $M_h^{crit}$ 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_{bh}^{ac}$ directly intervenes in equation (1): a higher $f_{bh}^{ac}$ – or a lower coupling factor $f_{bh}^w$ – 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_{bh}^{ac}$ or in the coupling factor $f_{bh}^w$ would mean less gas mass available for the subsequent time step, compensating the change in $f_{bh}^{ac}$ and $f_{bh}^w$. In this regime black hole growth is self-regulated. These parameters are in turn tuned to reproduce simultaneously both the AGN UV luminosity function and the BHMF.

# 3 THE MASS ASSEMBLY OF EARLY BLACK HOLES AND THEIR HOST HALOES

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 Section 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.

## 3.1 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 overshines 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 submillimetre 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 UV indicators (e.g. McLure & Jarvis 2002; Vestergaard 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 tdf4 scenarios, colour 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 haloes ($M_h \lesssim 10^9 M_\odot$) do not accrete at all, since SNeII feedback ejects all of the gas mass from the host galaxies. The double tail at $M_h \sim 10^{9-11} M_\odot$ is due to the presence of two types of black hole seeds: the upper part is composed of (mixed) black holes resulting from the mergers of DCBHs and SBHs, while the lower part is composed of SBHs. The black holes in this halo mass range show *stunted growth*, given that they accrete at a rate $\dot{M}_{bh} = 7.5 \times 10^{-5} \dot{M}_{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 (Eddington-limited phase); this phase ends when the gas-limited condition $M_*^{gf}(z) f_{bh}^{ac} (1 - \epsilon_r) < M_{Edd}(z)$ is met. This happens for the highest mass haloes with $M_h \gtrsim 10^{12} M_\odot$ at $z = 5$. Our black hole mass–halo mass relation is therefore described by a triple power law at $z \sim 5$: flat at $M_h \lesssim M_h^{crit}$, almost vertical at $M_h^{crit} \lesssim M_h \lesssim 10^{11.6} M_\odot$, and $\log(M_{bh}/M_\odot) \propto \log(M_h/M_\odot)^{1.4}$ at $M_h \gtrsim 10^{11.6} M_\odot$. On the other hand, notice that black holes at $z \gtrsim 7$ never enter
<figure><img src="image_4.png" /></figure>

**Figure 1.** The halo mass–black hole mass relation at $z=5, 7$, and $9$, as marked. The upper and lower panels show results for the *ins1* and *tdf4* cases, respectively. The different symbols in each panels show the black hole accretion channels as marked in the top bar (green filled circles represent the gas-limited accretion channel, blue filled triangles are for stunted accretion mode, and orange empty squares are for Eddington-limited accretion); the cyan filled diamonds represent DCBHs. In the $z=5$ panel, the solid, dashed, and dotted lines show the (average) observational results inferred for the local universe (Ferrarese 2002), for $0.1 \lesssim z \lesssim 0.3$ (Bandara, Crampton & Simard 2009), and theoretical predictions at $z=5$ (Croton 2009), respectively.

the last accretion phase. Finally, the slight increase of $M_{\mathrm{h}}^{\mathrm{crit}}$ with decreasing redshift naturally results in more massive haloes (and hence black holes) accreting at the Eddington rate with decreasing redshift.

We then look at results from the *tdf4* 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_{\mathrm{bh}}-M_{\mathrm{h}}$ relation in this scenario follows a similar slope as the *ins1* model for $M_{\mathrm{h}} \gtrsim 10^{11} M_{\odot}$ 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 *tdf4* model is lower by a factor of $\approx 50$. In addition, this model has a larger scatter (by about a factor of 2) 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 anticorrelate with the width of the scatter in the black hole mass–host mass relation.

We find that for haloes with $M_{\mathrm{h}} \gtrsim 10^{11} M_{\odot}$, for both models, while the slope of our $M_{\mathrm{bh}}-M_{\mathrm{h}}$ relation is consistent with a number of other (observational and theoretical) works, its normalization 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 the *James Webb Space Telescope (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_{\mathrm{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
<figure><img src="image_3.png" /></figure>

**Figure 1.** The halo mass–black hole mass relation at $z=5, 7$, and $9$, as marked. The upper and lower panels show results for the *ins1* and *tdf4* cases, respectively. The different symbols in each panels show the black hole accretion channels as marked in the top bar (green filled circles represent the gas-limited accretion channel, blue filled triangles are for stunted accretion mode, and orange empty squares are for Eddington-limited accretion); the cyan filled diamonds represent DCBHs. In the $z=5$ panel, the solid, dashed, and dotted lines show the (average) observational results inferred for the local universe (Ferrarese 2002), for $0.1 \lesssim z \lesssim 0.3$ (Bandara, Crampton & Simard 2009), and theoretical predictions at $z=5$ (Croton 2009), respectively.

the last accretion phase. Finally, the slight increase of $M_{\mathrm{h}}^{\mathrm{crit}}$ with decreasing redshift naturally results in more massive haloes (and hence black holes) accreting at the Eddington rate with decreasing redshift.

We then look at results from the *tdf4* 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_{\mathrm{bh}}-M_{\mathrm{h}}$ relation in this scenario follows a similar slope as the *ins1* model for $M_{\mathrm{h}} \gtrsim 10^{11} M_{\odot}$ 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 *tdf4* model is lower by a factor of $\approx 50$. In addition, this model has a larger scatter (by about a factor of 2) 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 anticorrelate with the width of the scatter in the black hole mass–host mass relation.

We find that for haloes with $M_{\mathrm{h}} \gtrsim 10^{11} M_{\odot}$, for both models, while the slope of our $M_{\mathrm{bh}}-M_{\mathrm{h}}$ relation is consistent with a number of other (observational and theoretical) works, its normalization 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 the *James Webb Space Telescope (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_{\mathrm{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
<figure><img src="image_5.png" /><figcaption>**Figure 2.** The stellar mass–black hole mass relation at $z = 5, 7$, and 9, as marked. The upper and lower panels show results for the *ins1* and *tdf4* cases, respectively. The different symbols in each panels show the black hole accretion channels as marked in the top bar (green filled circles represent the gas-limited accretion channel, blue filled triangles are for stunted accretion mode, and orange empty squares are for Eddington-limited accretion); the cyan filled diamonds represent DCBHs. In the $z = 5$ panel, we plot the observational results for local moderate-luminosity AGN (solid line) and high-mass ellipticals (long-dashed line) inferred in Reines & Volonteri (2015), the results of the SIMBA simulation at $z = 5$ from Thomas et al. (2019) (short-dashed line), and those observed at $z = 0–2$ from Suh et al. (2020) (grey shaded areas).</figcaption></figure>

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; Decarli et al. 2010; Merloni et al. 2010; Caplar, Lilly & Trakhtenbrot 2018), while other authors found no significant difference in the mass ratios of high- and low-redshift AGN, both with observational (Shields et al. 2003; Shankar, Bernardi & Haiman 2009; Salviander, Shields & Bonning 2015; Shen et al. 2015) and theoretical (Thomas et al. 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_{*} \sim 10^{8.5-9.5} M_{\odot}$, 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_{*} \gtrsim 10^{10} M_{\odot}$ sit ~1 dex higher with respect to the same low-redshift measurements and the relation for moderate-luminosity 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 reionization 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
<figure><img src="image_4.png" /><figcaption>**Figure 2.** The stellar mass–black hole mass relation at $z = 5, 7$, and 9, as marked. The upper and lower panels show results for the *ins1* and *tdf4* cases, respectively. The different symbols in each panels show the black hole accretion channels as marked in the top bar (green filled circles represent the gas-limited accretion channel, blue filled triangles are for stunted accretion mode, and orange empty squares are for Eddington-limited accretion); the cyan filled diamonds represent DCBHs. In the $z = 5$ panel, we plot the observational results for local moderate-luminosity AGN (solid line) and high-mass ellipticals (long-dashed line) inferred in Reines & Volonteri (2015), the results of the SIMBA simulation at $z = 5$ from Thomas et al. (2019) (short-dashed line), and those observed at $z = 0–2$ from Suh et al. (2020) (grey shaded areas).</figcaption></figure>

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; Decarli et al. 2010; Merloni et al. 2010; Caplar, Lilly & Trakhtenbrot 2018), while other authors found no significant difference in the mass ratios of high- and low-redshift AGN, both with observational (Shields et al. 2003; Shankar, Bernardi & Haiman 2009; Salviander, Shields & Bonning 2015; Shen et al. 2015) and theoretical (Thomas et al. 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_{*} \sim 10^{8.5-9.5} M_{\odot}$, 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_{*} \gtrsim 10^{10} M_{\odot}$ sit ~1 dex higher with respect to the same low-redshift measurements and the relation for moderate-luminosity 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 reionization 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
<figure><img src="image_6.png" /><figcaption>**Figure 3.** The redshift evolution of the (logarithm of the) average Eddington ratio $\bar{\lambda}_{\text{Edd}}$ for the different black hole mass bins marked. Solid and dot-dashed lines show results for the *ins1* and *tdf4* scenarios, respectively.</figcaption></figure>

massive galaxies above $10^{9.5} \text{ M}_\odot$ 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} \text{ M}_\odot$ 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 $\lambda_{\text{Edd}} = M_{\text{bh}}^{\text{ac}}/M_{\text{Edd}}$, which is indicative of the black hole growth efficiency. By construction, we have that $\lambda_{\text{Edd}} = 7.5 \times 10^{-5}$ for halo masses $M_h < M_h^{\text{crit}} \sim 10^{11.5} \text{ M}_\odot$, while $\lambda_{\text{Edd}}$ is generally between 0.1 and 1 for black holes in higher mass haloes. The results for the average Eddington ratios ($\bar{\lambda}_{\text{Edd}}$) as a function of redshift and for different black hole mass bins are shown in Fig. 3. Low-mass black holes ($M_{\text{bh}} \sim 10^{3-4} \text{ M}_\odot$) reside in progressively more massive haloes with decreasing $z$, eventually starting to accrete at the Eddington rate (also see Fig. 1), leading to an increases of $\bar{\lambda}_{\text{Edd}}$ from $10^{-3}$ at $z=9$ to 0.1 by $z=4$. Moving on, approximately 90 per cent of the black holes in the mass range $10^4-10^5 \text{ M}_\odot$ are in the Eddington-limited regime at $z=4$, with $\bar{\lambda}_{\text{Edd}} \approx 0.9$, while the rest of the black holes sits in low-mass haloes and accretes at $\bar{\lambda}_{\text{Edd}} = 7.5 \times 10^{-5}$. Black holes with masses $10^{5-6} \text{ M}_\odot$ generally accrete at $\bar{\lambda}_{\text{Edd}} = 1$ at all redshifts, as also seen from Fig. 1. As black holes grow above $10^6 \text{ M}_\odot$ and switch to the gas-limited accretion phase, the average Eddington ratio drops from $\bar{\lambda}_{\text{Edd}} = 1$ to progressively lower values $\sim 0.1$ (see e.g. Merloni 2004). This is because as $M_{\text{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 ratio distribution at high redshift, which extends to super-Eddington regimes (see for instance Kelly & Shen 2013). In this sense, the results shown in this plot are meant to portrait a qualitative more than quantitative trend.

Compared to the *ins1* model, in the *tdf4* scenario, the largest difference is seen for the lowest mass bins ($M_{\text{bh}} = 10^{3-5} \text{ M}_\odot$). We also see a slight increase in $\bar{\lambda}$ for $M_{\text{bh}} \gtrsim 10^6 \text{ M}_\odot$. This is because a

delay in the baryonic mass assembly in the *tdf4* 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_{\text{Edd}}(z) < M_{\text{bh}}^{\text{gf}}(z) f_{\text{bh}}^{\text{ac}}(1-\epsilon_r)$ is met for a longer amount of time. Additional delays in black holes mergers that are not taken into account here, such as the black hole binary inspiralling phase and the final parsec time-scale, might further increase the average Eddington ratio.

## 3.2 Joint assembly of black holes and their host haloes

We now study the relative contribution of mergers and accretion to the joint assembly of black holes and their host haloes at $z > 4$. For both black holes and haloes, 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 on to the black hole and halo, respectively. Here we explore results for 21 galaxies over three halo mass bins (at $z=4$) each: low-mass haloes ($M_h \sim 10^{9.9-10.1} \text{ M}_\odot$), intermediate-mass haloes ($M_h \sim 10^{11.65-11.85} \text{ M}_\odot$), and high-mass haloes ($M_h \sim 10^{13.3-13.5} \text{ M}_\odot$). We start by discussing results from the *ins1* scenario. Naturally, both the *ins1* and *tdf4* scenarios show the same results for halo assembly; the dynamical friction delay only affects the baryonic component.

Low-mass haloes (panel a of Fig. 4) start assembling their dark matter mass at $z \sim 14$. Most of the halo growth is driven by accretion down to $z \sim 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 \sim 13-4$. The accretion mode is subdominant given the negligible accretion rates ($\sim 7.5 \times 10^{-5} M_{\text{Edd}}$), and contributes less than 10 per cent of the total black hole mass by $z \sim 4$.

Intermediate-mass haloes (panel b of Fig. 4) naturally start assembling earlier, at $z \sim 17.5$. Accretion dominates the mass assembly until $z \sim 7$, below which point mergers start contributing equally. While mergers dominate the black hole mass build-up down to $z \sim 7$, thereafter Eddington-limited accretion takes over, dominating the black hole mass budget by $\sim 2$ orders of magnitude by $z=4$. By then, more than 95 per cent 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 haloes, 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 haloes is dominated by accretion until $z \sim 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) in the range $z \sim 17.5-9$; accretion dominates instead the later phases, contributing about 90 per cent of the final mass at $z \sim 4$. It is interesting also to notice that for high-mass haloes the variance in the merged components is generally lower, and significant only at $z \le 7$, when stochastic major mergers take place. Further, more than one progenitor is able to overcome $M_h^{\text{crit}}$, and on average only about 20 (75 per cent) of the mass in assembled in the major branch at $z \sim 10(4)$ as shown in Fig. 5.

In the *tdf4* 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 haloes).

To summarize, while black holes residing in low-mass haloes can grow only through (rare) mergers, those in intermediate- and high-mass mass haloes predominantly grow through mergers only in the
<figure><img src="image_5.png" /><figcaption>**Figure 3.** The redshift evolution of the (logarithm of the) average Eddington ratio $\bar{\lambda}_{\text{Edd}}$ for the different black hole mass bins marked. Solid and dot-dashed lines show results for the *ins1* and *tdf4* scenarios, respectively.</figcaption></figure>

massive galaxies above $10^{9.5} \text{ M}_\odot$ 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} \text{ M}_\odot$ 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 $\lambda_{\text{Edd}} = M_{\text{bh}}^{\text{ac}}/M_{\text{Edd}}$, which is indicative of the black hole growth efficiency. By construction, we have that $\lambda_{\text{Edd}} = 7.5 \times 10^{-5}$ for halo masses $M_h < M_h^{\text{crit}} \sim 10^{11.5} \text{ M}_\odot$, while $\lambda_{\text{Edd}}$ is generally between 0.1 and 1 for black holes in higher mass haloes. The results for the average Eddington ratios ($\bar{\lambda}_{\text{Edd}}$) as a function of redshift and for different black hole mass bins are shown in Fig. 3. Low-mass black holes ($M_{\text{bh}} \sim 10^{3-4} \text{ M}_\odot$) reside in progressively more massive haloes with decreasing $z$, eventually starting to accrete at the Eddington rate (also see Fig. 1), leading to an increases of $\bar{\lambda}_{\text{Edd}}$ from $10^{-3}$ at $z=9$ to 0.1 by $z=4$. Moving on, approximately 90 per cent of the black holes in the mass range $10^4-10^5 \text{ M}_\odot$ are in the Eddington-limited regime at $z=4$, with $\bar{\lambda}_{\text{Edd}} \approx 0.9$, while the rest of the black holes sits in low-mass haloes and accretes at $\bar{\lambda}_{\text{Edd}} = 7.5 \times 10^{-5}$. Black holes with masses $10^{5-6} \text{ M}_\odot$ generally accrete at $\bar{\lambda}_{\text{Edd}} = 1$ at all redshifts, as also seen from Fig. 1. As black holes grow above $10^6 \text{ M}_\odot$ and switch to the gas-limited accretion phase, the average Eddington ratio drops from $\bar{\lambda}_{\text{Edd}} = 1$ to progressively lower values $\sim 0.1$ (see e.g. Merloni 2004). This is because as $M_{\text{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 ratio 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 *tdf4* scenario, the largest difference is seen for the lowest mass bins ($M_{\text{bh}} = 10^{3-5} \text{ M}_\odot$). We also see a slight increase in $\bar{\lambda}$ for $M_{\text{bh}} \gtrsim 10^6 \text{ M}_\odot$. This is because a

delay in the baryonic mass assembly in the *tdf4* 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_{\text{Edd}}(z) < M_{\text{bh}}^{\text{gf}}(z) f_{\text{bh}}^{\text{ac}} (1 - \epsilon_r)$ is met for a longer amount of time. Additional delays in black holes mergers that are not taken into account here, such as the black hole binary inspiralling phase and the final parsec time-scale, might further increase the average Eddington ratio.

## 3.2 Joint assembly of black holes and their host haloes

We now study the relative contribution of mergers and accretion to the joint assembly of black holes and their host haloes at $z > 4$. For both black holes and haloes, 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 on to the black hole and halo, respectively. Here we explore results for 21 galaxies over three halo mass bins (at $z=4$) each: low-mass haloes ($M_h \sim 10^{9.9-10.1} \text{ M}_\odot$), intermediate-mass haloes ($M_h \sim 10^{11.65-11.85} \text{ M}_\odot$), and high-mass haloes ($M_h \sim 10^{13.3-13.5} \text{ M}_\odot$). We start by discussing results from the *ins1* scenario. Naturally, both the *ins1* and *tdf4* scenarios show the same results for halo assembly; the dynamical friction delay only affects the baryonic component.

Low-mass haloes (panel a of Fig. 4) start assembling their dark matter mass at $z \sim 14$. Most of the halo growth is driven by accretion down to $z \sim 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 \sim 13-4$. The accretion mode is subdominant given the negligible accretion rates ($\sim 7.5 \times 10^{-5} M_{\text{Edd}}$), and contributes less than 10 per cent of the total black hole mass by $z \sim 4$.

Intermediate-mass haloes (panel b of Fig. 4) naturally start assembling earlier, at $z \sim 17.5$. Accretion dominates the mass assembly until $z \sim 7$, below which point mergers start contributing equally. While mergers dominate the black hole mass build-up down to $z \sim 7$, thereafter Eddington-limited accretion takes over, dominating the black hole mass budget by $\sim 2$ orders of magnitude by $z=4$. By then, more than 95 per cent 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 haloes, 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 haloes is dominated by accretion until $z \sim 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) in the range $z \sim 17.5-9$; accretion dominates instead the later phases, contributing about 90 per cent of the final mass at $z \sim 4$. It is interesting also to notice that for high-mass haloes the variance in the merged components is generally lower, and significant only at $z \le 7$, when stochastic major mergers take place. Further, more than one progenitor is able to overcome $M_h^{\text{crit}}$, and on average only about 20 (75 per cent) of the mass in assembled in the major branch at $z \sim 10(4)$ as shown in Fig. 5.

In the *tdf4* 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 haloes).

To summarize, while black holes residing in low-mass haloes can grow only through (rare) mergers, those in intermediate- and high-mass mass haloes predominantly grow through mergers only in the
<figure><img src="image_7.png" /><figcaption>**Figure 4.** The cumulative mass assembly of the black hole and the host dark matter halo as a result of mergers ($M_{bh}^{mer}$ and $M_{h}^{mer}$, represented respectively by the orange and the cyan lines) and accretion ($M_{bh}^{acc}$ and $M_{h}^{acc}$, red and blue lines) as a function of redshift. We show results for $z = 4$ halo masses in the range $M_{h} \sim 10^{9.9-10.1} M_{\odot}$ (upper panel), $M_{h} \sim 10^{11.65-11.85} M_{\odot}$ (middle panel), and $M_{h} \sim 10^{13.3-13.5} M_{\odot}$ (lower panel), as marked, averaging over 21 haloes in each panel. The solid and dashed lines in each panel show results for the *ins1* and *tdf4* cases, respectively; shaded regions show the $1\sigma$ variance for the *ins1* scenario. As expected, while the halo assembly is identical in both models, black hole assembly is delayed in the *tdf4* model.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** The redshift evolution of the ratio of the cumulative mass accreted by black holes in the major branch with respect to that accreted in the whole merger tree of a galaxy (averaged over 21 haloes) for $z = 4$ halo masses of $M_{h} \sim 10^{11.65-11.85} M_{\odot}$ (red lines) and $M_{h} \sim 10^{13.3-13.5} M_{\odot}$ (blue lines), as marked. The solid and dashed lines show results for the *ins1* and *tdf4* cases, respectively; shaded regions show the $1\sigma$ variance for the *ins1* scenario.</figcaption></figure>

early phases, before switching to an accretion-dominated growth phase. As expected, mergers are less important in the *tdf4* scenario, contributing $\sim 10^{-3.5}$ per cent ($10^{-2}$ per cent) of the total mass in intermediate (high) mass haloes.

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 haloes (when present) are starved throughout the redshift range $z = 14–4$, driving the mass ratio to a constant decrease. Intermediate-mass haloes show a flat evolution at $z > 13$, when the black hole seed formation efficiency is still 100 per cent (i.e. every newly born halo has a black hole). At $z \sim 13–6$, the ratio drops as new haloes form without black hole seeds. At $z \lesssim 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} \sim 10^{-7.4}$ at $z = 6$ to $M_{bh}/M_{h} \sim 10^{-4.5}$ at $z = 4$. While the evolution of the highest mass haloes is qualitatively similar to intermediate-mass haloes, their black holes already enter the Eddington-limited growth channel at $z \approx 8$. Sustaining this growth for a longer time, they show a mass ratio that increases from $M_{bh}/M_{h} \sim 10^{-6.8}$ at $z = 8$ to $M_{bh}/M_{h} \sim 10^{-3.6}$ at $z = 4$. Again, intermediate-mass haloes show a larger variance with respect to high-mass haloes as a result of the diverse merging and accretion histories of their black holes.

# 4 KEY PROPERTIES OF THE BLACK HOLE POPULATION IN EARLY GALAXIES

We now explore the key properties of black holes in early galaxies, specifically focusing on their occupation fraction and duty cycles in this section.

## 4.1 The black hole occupation fraction

We define the black hole occupation fraction ($f_{bh}$) as the fraction of haloes hosting a *central* black hole. In the case of merging systems,
<figure><img src="image_6.png" /><figcaption>**Figure 4.** The cumulative mass assembly of the black hole and the host dark matter halo as a result of mergers ($M_{bh}^{mer}$ and $M_{h}^{mer}$, represented respectively by the orange and the cyan lines) and accretion ($M_{bh}^{acc}$ and $M_{h}^{acc}$, red and blue lines) as a function of redshift. We show results for $z = 4$ halo masses in the range $M_{h} \sim 10^{9.9-10.1} M_{\odot}$ (upper panel), $M_{h} \sim 10^{11.65-11.85} M_{\odot}$ (middle panel), and $M_{h} \sim 10^{13.3-13.5} M_{\odot}$ (lower panel), as marked, averaging over 21 haloes in each panel. The solid and dashed lines in each panel show results for the *ins1* and *tdf4* cases, respectively; shaded regions show the $1\sigma$ variance for the *ins1* scenario. As expected, while the halo assembly is identical in both models, black hole assembly is delayed in the *tdf4* model.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** The redshift evolution of the ratio of the cumulative mass accreted by black holes in the major branch with respect to that accreted in the whole merger tree of a galaxy (averaged over 21 haloes) for $z = 4$ halo masses of $M_{h} \sim 10^{11.65-11.85} M_{\odot}$ (red lines) and $M_{h} \sim 10^{13.3-13.5} M_{\odot}$ (blue lines), as marked. The solid and dashed lines show results for the *ins1* and *tdf4* cases, respectively; shaded regions show the $1\sigma$ variance for the *ins1* scenario.</figcaption></figure>

early phases, before switching to an accretion-dominated growth phase. As expected, mergers are less important in the *tdf4* scenario, contributing $\sim 10^{-3.5}$ per cent ($10^{-2}$ per cent) of the total mass in intermediate (high) mass haloes.

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 haloes (when present) are starved throughout the redshift range $z = 14–4$, driving the mass ratio to a constant decrease. Intermediate-mass haloes show a flat evolution at $z > 13$, when the black hole seed formation efficiency is still 100 per cent (i.e. every newly born halo has a black hole). At $z \sim 13–6$, the ratio drops as new haloes form without black hole seeds. At $z \lesssim 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} \sim 10^{-7.4}$ at $z = 6$ to $M_{bh}/M_{h} \sim 10^{-4.5}$ at $z = 4$. While the evolution of the highest mass haloes is qualitatively similar to intermediate-mass haloes, their black holes already enter the Eddington-limited growth channel at $z \approx 8$. Sustaining this growth for a longer time, they show a mass ratio that increases from $M_{bh}/M_{h} \sim 10^{-6.8}$ at $z = 8$ to $M_{bh}/M_{h} \sim 10^{-3.6}$ at $z = 4$. Again, intermediate-mass haloes show a larger variance with respect to high-mass haloes as a result of the diverse merging and accretion histories of their black holes.

# 4 KEY PROPERTIES OF THE BLACK HOLE POPULATION IN EARLY GALAXIES

We now explore the key properties of black holes in early galaxies, specifically focusing on their occupation fraction and duty cycles in this section.

## 4.1 The black hole occupation fraction

We define the black hole occupation fraction ($f_{bh}$) as the fraction of haloes hosting a *central* black hole. In the case of merging systems,
<figure><img src="image_9.png" /><figcaption>**Figure 6.** As a function of $z$, we show the ratio of the cumulative black hole mass to the cumulative halo mass (averaged over 21 haloes) for $z = 4$ halo masses of $M_h \sim 10^{9.6-10.1} M_\odot$ (solid dark blue line), $M_h \sim 10^{11.65-11.85} M_\odot$ (dashed red line), and $M_h \sim 10^{13.3-13.5} M_\odot$ (dotted light blue line), as marked. The upper and lower panels show results for the *ins1* and *tdf4* cases, respectively; shaded regions show the associated $1\sigma$ variance.</figcaption></figure>

this refers to the fraction of central galaxies hosting a black hole. We note that this specification is important only in the *tdf4* scenario; in the *ins1* case galaxies merge at the same time as their haloes. In our model, haloes are populated with SBH seeds down to $z = 13$. Lower $z$ haloes contain a black hole if they are either seeded with a DCBH or if they gain one through mergers. The latter (and 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 per cent of haloes contain a central black hole) for a minimum halo mass of $M_h \sim 10^{8.4} M_\odot$ at $z \sim 12$, which increases to $M_h \sim 10^{9.5} M_\odot$ by $z = 5$. Reaching $f_{bh} = 1$ naturally requires a higher halo mass of $M_h \sim 10^{8.9} (10^{10.5}) M_\odot$ at $z \sim 12 (5)$. As seen from the same plot, there is no sensible difference in the *ins1* and *tdf4* 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 *tdf4* 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 deconstruct $f_{bh}$ as a function of halo mass for different black hole mass bins, and show the results in Fig. 8. First, 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

<figure><img src="image_10.png" /><figcaption>**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 *tdf4* (dashed lines) scenarios, as marked.</figcaption></figure>

haloes. In terms of the most massive black holes, while by $z \simeq 9$ the most massive haloes ($M_h \sim 10^{11.5} M_\odot$) host black holes as massive as $\sim 10^6 M_\odot$, this increases to $M_{bh} \sim 10^{9.5} M_\odot$ by $z = 5$ for haloes with $M_h \sim 10^{13.4} M_\odot$. 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$ haloes with $M_h \sim 10^{11-12} M_\odot$ host black holes with $M_{bh} \sim 10^4-10^6 M_\odot$, at $z = 5$ the same halo masses host black holes spanning four orders of magnitude in mass, between $10^{3.5}$ and $10^{7.5} M_\odot$. This is mainly due to the black holes residing in these haloes entering the Eddington-limited 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_h^{crit}$, are amplified by the exponential growth of black holes during the Eddington-limited growth phase.

In the *tdf4* case, haloes 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 \gtrsim 9$ in intermediate- and high-mass haloes. As accretion overtakes the mass build-up at lower $z$, results from the *ins1* and *tdf4* models come into closer agreement.

## 4.2 Time-weighed duty cycle

We now explore the black hole duty cycle defined as the fractional lifetime a black hole spent in the stunted, Eddington-limited, and gas-limited accretion modes. This is done for all the $z = 4$ haloes that host a black hole, accounting for all of the progenitors along the merger tree, and the results are shown in Fig. 9. The time-weighed duty cycle for each accretion mode is then defined as $\tau_i = t_i/t_{tot}$, where $t_i$ indicates the time spent in each of the three accretion channels (*i*) and $t_{tot}$ 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 haloes, where they are subject to stunted accretion. In particular, black hole progenitors of haloes with $M_h|_{z=4} < M_h^{crit}$ spend all of their lifetime in a stunted accretion mode or – if they are in SN-feedback-dominated haloes (which are more common among the progenitors of high-mass
<figure><img src="image_8.png" /><figcaption>**Figure 6.** As a function of $z$, we show the ratio of the cumulative black hole mass to the cumulative halo mass (averaged over 21 haloes) for $z = 4$ halo masses of $M_h \sim 10^{9.6-10.1} M_\odot$ (solid dark blue line), $M_h \sim 10^{11.65-11.85} M_\odot$ (dashed red line), and $M_h \sim 10^{13.3-13.5} M_\odot$ (dotted light blue line), as marked. The upper and lower panels show results for the *ins1* and *tdf4* cases, respectively; shaded regions show the associated $1\sigma$ variance.</figcaption></figure>

this refers to the fraction of central galaxies hosting a black hole. We note that this specification is important only in the *tdf4* scenario; in the *ins1* case galaxies merge at the same time as their haloes. In our model, haloes are populated with SBH seeds down to $z = 13$. Lower $z$ haloes contain a black hole if they are either seeded with a DCBH or if they gain one through mergers. The latter (and 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 per cent of haloes contain a central black hole) for a minimum halo mass of $M_h \sim 10^{8.4} M_\odot$ at $z \sim 12$, which increases to $M_h \sim 10^{9.5} M_\odot$ by $z = 5$. Reaching $f_{bh} = 1$ naturally requires a higher halo mass of $M_h \sim 10^{8.9} (10^{10.5}) M_\odot$ at $z \sim 12 (5)$. As seen from the same plot, there is no sensible difference in the *ins1* and *tdf4* 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 *tdf4* 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 deconstruct $f_{bh}$ as a function of halo mass for different black hole mass bins, and show the results in Fig. 8. First, 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

<figure><img src="image_9.png" /><figcaption>**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 *tdf4* (dashed lines) scenarios, as marked.</figcaption></figure>

haloes. In terms of the most massive black holes, while by $z \simeq 9$ the most massive haloes ($M_h \sim 10^{11.5} M_\odot$) host black holes as massive as $\sim 10^6 M_\odot$, this increases to $M_{bh} \sim 10^{9.5} M_\odot$ by $z = 5$ for haloes with $M_h \sim 10^{13.4} M_\odot$. 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$ haloes with $M_h \sim 10^{11-12} M_\odot$ host black holes with $M_{bh} \sim 10^4-10^6 M_\odot$, at $z = 5$ the same halo masses host black holes spanning four orders of magnitude in mass, between $10^{3.5}$ and $10^{7.5} M_\odot$. This is mainly due to the black holes residing in these haloes entering the Eddington-limited 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_h^{crit}$, are amplified by the exponential growth of black holes during the Eddington-limited growth phase.

In the *tdf4* case, haloes 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 \gtrsim 9$ in intermediate- and high-mass haloes. As accretion overtakes the mass build-up at lower $z$, results from the *ins1* and *tdf4* models come into closer agreement.

## 4.2 Time-weighed duty cycle

We now explore the black hole duty cycle defined as the fractional lifetime a black hole spent in the stunted, Eddington-limited, and gas-limited accretion modes. This is done for all the $z = 4$ haloes that host a black hole, accounting for all of the progenitors along the merger tree, and the results are shown in Fig. 9. The time-weighed duty cycle for each accretion mode is then defined as $\tau_i = t_i/t_{tot}$, where $t_i$ indicates the time spent in each of the three accretion channels (*i*) and $t_{tot}$ 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 haloes, where they are subject to stunted accretion. In particular, black hole progenitors of haloes with $M_h|_{z=4} < M_h^{crit}$ spend all of their lifetime in a stunted accretion mode or – if they are in SN-feedback-dominated haloes (which are more common among the progenitors of high-mass
<figure><img src="image_11.png" /><figcaption>**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.</figcaption></figure>

haloes) – not accreting at all. For $M_h|_{z=4} \gtrsim 10^{11.5}$ M$_\odot$ haloes, 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 haloes in their merger trees. The black holes in intermediate-mass haloes spend at most 5 per cent and 1 per cent of their lifetime accreting in the Eddington- and gas-limited modes, respectively; these drop to 0.1 per cent and 1 per cent for the highest mass haloes.

The most important difference upon implementing dynamical friction delay and UV feedback is that the black hole lifetime spent in non-accreting haloes increases by at least a factor of 3 for the highest mass haloes (and up to an order of magnitude for the lowest mass haloes) due to UV suppression of their gas mass.

Finally, the variations (up to one order of magnitude) between the duty cycles of haloes of very similar final masses are driven by the assembly histories of both black holes and their host haloes, which determine both the halo and black hole masses, as well as the amount of gas mass available for accretion on to the black hole.

## 4.3 Black hole predictions at $z \ge 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-hand 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-hand panel). We use our intrinsic and obscured XLFs 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 realizations 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_h^{crit}$ 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$ BHMF. 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} \lesssim 10^7$ M$_\odot$ follow a power law with their number densities being virtually constant at $z \sim 5-7$. On the other hand, for $M_{bh} \gtrsim 10^7$ M$_\odot$, the number density of high-mass black holes shows an exponential cut-off, especially at $z \sim 5-6$, in addition to showing a rapid evolution between $z \sim 6$ and 7. In particular, the mass bin $10^6 \lesssim M_{bh} \lesssim 10^7$ M$_\odot$ experiences the greatest variation of number
<figure><img src="image_10.png" /><figcaption>**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.</figcaption></figure>

haloes) – not accreting at all. For $M_h|_{z=4} \gtrsim 10^{11.5}$ M$_\odot$ haloes, 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 haloes in their merger trees. The black holes in intermediate-mass haloes spend at most 5 per cent and 1 per cent of their lifetime accreting in the Eddington- and gas-limited modes, respectively; these drop to 0.1 per cent and 1 per cent for the highest mass haloes.

The most important difference upon implementing dynamical friction delay and UV feedback is that the black hole lifetime spent in non-accreting haloes increases by at least a factor of 3 for the highest mass haloes (and up to an order of magnitude for the lowest mass haloes) due to UV suppression of their gas mass.

Finally, the variations (up to one order of magnitude) between the duty cycles of haloes of very similar final masses are driven by the assembly histories of both black holes and their host haloes, which determine both the halo and black hole masses, as well as the amount of gas mass available for accretion on to the black hole.

## 4.3 Black hole predictions at $z \ge 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-hand 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-hand panel). We use our intrinsic and obscured XLFs 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 realizations 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_h^{crit}$ 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$ BHMF. 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} \lesssim 10^7$ M$_\odot$ follow a power law with their number densities being virtually constant at $z \sim 5-7$. On the other hand, for $M_{bh} \gtrsim 10^7$ M$_\odot$, the number density of high-mass black holes shows an exponential cut-off, especially at $z \sim 5-6$, in addition to showing a rapid evolution between $z \sim 6$ and 7. In particular, the mass bin $10^6 \lesssim M_{bh} \lesssim 10^7$ M$_\odot$ experiences the greatest variation of number
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Logarithm of the fraction ($\tau_i$) of the cumulative lifetime of the black hole spent in the stunted (blue line), Eddington-limited (cyan line), and gas-limited (red line) accretion modes as a function of the $z = 4$ halo mass ($M_h|_{z=4}$). The orange line shows black holes in SNeII-feedback-dominated galaxies that do not accrete at all.</figcaption></figure>

density in the shortest time-scale, increasing by almost three orders of magnitude between $z = 7$ and 6 (i.e. in less than 300 Myr). 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} \sim 10^{4.5}$ and $\sim 10^7 M_{\odot}$ have very similar number densities. The exponential growth of black holes in the Eddington-limited accretion regime (i.e. all black holes with $M_{bh} \gtrsim 10^4 M_{\odot}$) makes them move from one black hole mass bin to the next 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_{bh} < 10^7 M_{\odot}$ at $z = 5$ is a direct effect of the critical halo mass $M_{\text{crit}}$ imposed in our model.

As noted before, we reiterate that the black hole accretion channel strongly depends on its mass: while low-mass black holes with $M_{bh} \sim 10^4 M_{\odot}$ accrete at very low ($7.5 \times 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 (Volonteri et al. 2016) and the AGN synthesis model by Merloni & Heinz (2008), although these models seem to hint to a characteristic knee mass lower by $\approx 0.5$ dex. Finally at $z = 5-6$, our results are within the range bracketed by the

GALFORM (Griffin et al. 2019) SAM (higher than our BHMF by more than 1 dex at $z = 6$ for $10^6 \lesssim M_{bh} \lesssim 10^8 M_{\odot}$) and the *Dragons* SAM (Qin et al. 2017), which predicts a lower number density (by about 0.5 dex for $10^6 \lesssim M_{bh} \lesssim 10^8 M_{\odot}$ 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 BHMF at $z = 5$, but at the same time finding a slight mismatch between our XLF and the one inferred at $z = 5$ by Ueda et al. (2014) at a luminosity of $L_X \sim 10^{44}$ erg s$^{-1}$, might point to a tension between these results. Nevertheless, the different AGN samples, selection criteria, and corrections used in Ueda et al. (2014) and Willott et al. (2010) to compute, respectively, the XLF and the BHMF do not allow us to draw any solid conclusion about this.

# 5 CONCLUSIONS AND DISCUSSION

In this work, we have used the *Delphi* SAM to study the joint assembly of high-redshift ($z > 4$) black holes and their host haloes.
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Logarithm of the fraction ($\tau_i$) of the cumulative lifetime of the black hole spent in the stunted (blue line), Eddington-limited (cyan line), and gas-limited (red line) accretion modes as a function of the $z = 4$ halo mass ($M_h|_{z=4}$). The orange line shows black holes in SNeII-feedback-dominated galaxies that do not accrete at all.</figcaption></figure>

density in the shortest time-scale, increasing by almost three orders of magnitude between $z = 7$ and 6 (i.e. in less than 300 Myr). 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} \sim 10^{4.5}$ and $\sim 10^7 M_{\odot}$ have very similar number densities. The exponential growth of black holes in the Eddington-limited accretion regime (i.e. all black holes with $M_{bh} \gtrsim 10^4 M_{\odot}$) makes them move from one black hole mass bin to the next 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_{bh} < 10^7 M_{\odot}$ at $z = 5$ is a direct effect of the critical halo mass $M_{\text{crit}}$ imposed in our model.

As noted before, we reiterate that the black hole accretion channel strongly depends on its mass: while low-mass black holes with $M_{bh} \sim 10^4 M_{\odot}$ accrete at very low ($7.5 \times 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 (Volonteri et al. 2016) and the AGN synthesis model by Merloni & Heinz (2008), although these models seem to hint to a characteristic knee mass lower by $\approx 0.5$ dex. Finally at $z = 5-6$, our results are within the range bracketed by the

GALFORM (Griffin et al. 2019) SAM (higher than our BHMF by more than 1 dex at $z = 6$ for $10^6 \lesssim M_{bh} \lesssim 10^8 M_{\odot}$) and the *Dragons* SAM (Qin et al. 2017), which predicts a lower number density (by about 0.5 dex for $10^6 \lesssim M_{bh} \lesssim 10^8 M_{\odot}$ 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 BHMF at $z = 5$, but at the same time finding a slight mismatch between our XLF and the one inferred at $z = 5$ by Ueda et al. (2014) at a luminosity of $L_X \sim 10^{44} \text{ erg s}^{-1}$, might point to a tension between these results. Nevertheless, the different AGN samples, selection criteria, and corrections used in Ueda et al. (2014) and Willott et al. (2010) to compute, respectively, the XLF and the BHMF do not allow us to draw any solid conclusion about this.

## 5 CONCLUSIONS AND DISCUSSION

In this work, we have used the *Delphi* SAM to study the joint assembly of high-redshift ($z > 4$) black holes and their host haloes.
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Evolution of the AGN X-ray luminosity function (XLF) at $z = 5, 7$, and 8. The black shaded areas represent our results, for which the upper and lower limit are set, respectively, by the intrinsic and the obscured XLF, where the obscured fraction as a function of luminosity is taken from Ueda et al. (2014). The red error bars represent the $1\sigma$ uncertainties on our lower limit. The other coloured shaded areas represent results from various models taken from fig. 3 of Amarantidis et al. (2019), their upper and lower limits corresponding to the minimum and maximum obscuration effect. The orange solid line shows the observational results from Ueda et al. (2014).</figcaption></figure>

We implement a physically motivated black hole growth mechanism that includes stunted growth in (gas-poor) low-mass haloes while high-mass haloes above the critical halo mass $M_{\text{crit}}(z)$ can accrete at the Eddington rate. After tuning the parameters of the model for the fiducial *ins1* scenario (see Table 1), we also consider a more realistic second case (*tdf4*) in which the mergers of the baryonic components are delayed by a dynamical friction time-scale, a maximal UV feedback scenario is implemented and the LW background threshold for DCBH seed formation is higher. We summarize here the key results.

(i) First of all, as a consequence of our seeding mechanism and accretion model, in haloes with $M_h > M_{\text{crit}}$ the black hole mass is assembled mostly via mergers until $z \sim 8$. Notably, this remains true even when the delay in mergers due to dynamical friction is added to the model. In fact, at such high redshift the mass ratios between haloes are on average closer to unity, and the merging time-scales are shorter. At $z < 8$ accretion takes over.

(ii) By $z = 4$ we find that in haloes with $M_h \approx 10^{11.75}$ ($10^{13.4}$) $M_{\odot}$ up to 95 per cent (75 per cent) of the final black hole mass is accreted within the major branch of the merger tree, while the rest comes from mergers with black holes hosted in secondary branches. While bigger haloes have more progenitors growing above $M_{\text{crit}}$, in intermediate-mass haloes only the major branch of the tree is able to fuel Eddington-limited accretion.

(iii) Finally, the average Eddington ratio $\bar{\lambda}_{\text{Edd}}$ seems to have a different evolution for black holes of different masses: for $M_{\text{bh}} \gtrsim 10^5 M_{\odot}$ we have $\bar{\lambda}_{\text{Edd}} \approx 1$ at $z > 6$, while it decreases down to $\approx 0.1$ at lower redshift. Low-mass black holes with $M_{\text{bh}} < 10^5 M_{\odot}$, instead,

show an opposite trend, increasing from $\bar{\lambda}_{\text{Edd}} \ll 1$ at high redshift to $\bar{\lambda}_{\text{Edd}} \approx 0.1-1$ at $z = 4$.

We tested our black hole seeding mechanism, finding that varying the initial DCBH seed mass our results remain effectively unchanged. This is in agreement with previous 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 disc (Dotti et al. 2010). We also point out that the black hole merging time-scales 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 LW background necessary to form DCBH seeds is still very poorly
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Evolution of the AGN X-ray luminosity function (XLF) at $z = 5, 7$, and 8. The black shaded areas represent our results, for which the upper and lower limit are set, respectively, by the intrinsic and the obscured XLF, where the obscured fraction as a function of luminosity is taken from Ueda et al. (2014). The red error bars represent the $1\sigma$ uncertainties on our lower limit. The other coloured shaded areas represent results from various models taken from fig. 3 of Amarantidis et al. (2019), their upper and lower limits corresponding to the minimum and maximum obscuration effect. The orange solid line shows the observational results from Ueda et al. (2014).</figcaption></figure>

We implement a physically motivated black hole growth mechanism that includes stunted growth in (gas-poor) low-mass haloes while high-mass haloes above the critical halo mass $M_{\text{crit}}^h(z)$ can accrete at the Eddington rate. After tuning the parameters of the model for the fiducial *ins1* scenario (see Table 1), we also consider a more realistic second case (*tdf4*) in which the mergers of the baryonic components are delayed by a dynamical friction time-scale, a maximal UV feedback scenario is implemented and the LW background threshold for DCBH seed formation is higher. We summarize here the key results.

(i) First of all, as a consequence of our seeding mechanism and accretion model, in haloes with $M_h > M_{\text{crit}}^h$ the black hole mass is assembled mostly via mergers until $z \sim 8$. Notably, this remains true even when the delay in mergers due to dynamical friction is added to the model. In fact, at such high redshift the mass ratios between haloes are on average closer to unity, and the merging time-scales are shorter. At $z < 8$ accretion takes over.

(ii) By $z = 4$ we find that in haloes with $M_h \approx 10^{11.75} (10^{13.4}) M_{\odot}$ up to 95 per cent (75 per cent) of the final black hole mass is accreted within the major branch of the merger tree, while the rest comes from mergers with black holes hosted in secondary branches. While bigger haloes have more progenitors growing above $M_{\text{crit}}^h$, in intermediate-mass haloes only the major branch of the tree is able to fuel Eddington-limited accretion.

(iii) Finally, the average Eddington ratio $\bar{\lambda}_{\text{Edd}}$ seems to have a different evolution for black holes of different masses: for $M_{\text{bh}} \gtrsim 10^5 M_{\odot}$ we have $\bar{\lambda}_{\text{Edd}} \approx 1$ at $z > 6$, while it decreases down to $\approx 0.1$ at lower redshift. Low-mass black holes with $M_{\text{bh}} < 10^5 M_{\odot}$, instead,

show an opposite trend, increasing from $\bar{\lambda}_{\text{Edd}} \ll 1$ at high redshift to $\bar{\lambda}_{\text{Edd}} \approx 0.1-1$ at $z = 4$.

We tested our black hole seeding mechanism, finding that varying the initial DCBH seed mass our results remain effectively unchanged. This is in agreement with previous 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 disc (Dotti et al. 2010). We also point out that the black hole merging time-scales 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 LW background necessary to form DCBH seeds is still very poorly
<figure><img src="image_14.png" /></figure>

**Figure 11.** Redshift evolution of the black hole mass function (BHMF) for the *ins1* scenario (solid cyan line). This is deconstructed to show the mass function for black holes accreting in the stunted mode (dashed blue line), Eddington-limited mode (long-dashed red line), and the gas-limited modes (green dotted line). We also show results from other theoretical and observational works at *z* = 5 from the EAGLE simulations (Rosas-Guevara et al. 2016, black solid line), the DRAGONS model (Qin et al. 2017), the Horizon-AGN simulation (Volonteri et al. 2016, green shaded area), and from Merloni & Heinz (2008, purple shaded area); and at *z* = 6 we show theoretical results from the GALFORM model (Griffin et al. 2019, green dashed line) and from Volonteri & Reines (2016, black dashed line) and observational results (Willott et al. 2010, orange shaded area).

constrained in the literature, and values of $\alpha$ 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_{\odot}$ 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.

## ACKNOWLEDGEMENTS

We thank the anonymous referee for helping us to improve the clarity of the paper. OP and PD acknowledge support from the European Commission's and University of Groningen's CO-FUND Rosalind Franklin program. PD also acknowledges support from the European Research Council's starting grant ERC StG-717001 ('DELPHI') and from the NWO grant 016.VIDI.189.162 ('ODIN'). TRC acknowledges the support of the Department of Atomic Energy, Government of India, under project no. 12-R&D-TFR-5.02-0700.

## DATA AVAILABILITY

The data underlying this paper, which have been produced by our model, are available on reasonable request to the corresponding author.

## REFERENCES

* Amarantidis S. et al., 2019, MNRAS, 485, 2694

* Anglés-Alcázar D., Faucher-Giguère C.-A., Quataert E., Hopkins P. F., Feldmann R., Torrey P., Wetzel A., Kereš D., 2017, MNRAS, 472, L109

* Bandara K., Crampton D., Simard L., 2009, ApJ, 704, 1135

* Barausse E., 2012, MNRAS, 423, 2533

* Bañados E. et al., 2018, Nature, 553, 473

* Beckmann R. S., Devriendt J., Slyz A., 2019, MNRAS, 483, 3488

* Begelman M. C., Volonteri M., Rees M. J., 2006, MNRAS, 370, 289

* Blecha L. et al., 2016, MNRAS, 456, 961

* Boco L., Lapi A., Danese L., 2020, ApJ, 891, 94

* Bower R. G., Benson A. J., Malbon R., Helly J. C., Frenk C. S., Baugh C. M., Cole S., Lacey C. G., 2006, MNRAS, 370, 645

* Bower R. G., Schaye J., Frenk C. S., Theuns T., Schaller M., Crain R. A., McAlpine S., 2017, MNRAS, 465, 32

* Bromm V., Loeb A., 2003, ApJ, 596, 34

* Caplar N., Lilly S. J., Trakhtenbrot B., 2018, ApJ, 867, 148

* Cattaneo A., 2002, MNRAS, 333, 353

* Costa T., Sijacki D., Trenti M., Haehnelt M. G., 2014, MNRAS, 439, 2146

* Croton D. J., 2009, MNRAS, 394, 1109

* Croton D. J. et al., 2006, MNRAS, 365, 11

* Dayal P., Ferrara A., Dunlop J. S., Pacucci F., 2014, MNRAS, 445, 2545

* Dayal P., Rossi E. M., Shiralilou B., Piana O., Choudhury T. R., Volonteri M., 2019, MNRAS, 486, 2336

* Dayal P. et al., 2020, MNRAS, 495, 3065

* Decarli R., Falomo R., Treves A., Labita M., Kotilainen J. K., Scarpa R., 2010, MNRAS, 402, 2453

* Decarli R. et al., 2018, ApJ, 854, 97

* Di Matteo T., Croft R. A. C., Feng Y., Waters D., Wilkins S., 2017, MNRAS, 467, 4243
<figure><img src="image_13.png" /></figure>

**Figure 11.** Redshift evolution of the black hole mass function (BHMF) for the *ins1* scenario (solid cyan line). This is deconstructed to show the mass function for black holes accreting in the stunted mode (dashed blue line), Eddington-limited mode (long-dashed red line), and the gas-limited modes (green dotted line). We also show results from other theoretical and observational works at *z* = 5 from the EAGLE simulations (Rosas-Guevara et al. 2016, black solid line), the DRAGONS model (Qin et al. 2017), the Horizon-AGN simulation (Volonteri et al. 2016, green shaded area), and from Merloni & Heinz (2008, purple shaded area); and at *z* = 6 we show theoretical results from the GALFORM model (Griffin et al. 2019, green dashed line) and from Volonteri & Reines (2016, black dashed line) and observational results (Willott et al. 2010, orange shaded area).

constrained in the literature, and values of $\alpha$ 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_{\odot}$ 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.

## ACKNOWLEDGEMENTS

We thank the anonymous referee for helping us to improve the clarity of the paper. OP and PD acknowledge support from the European Commission's and University of Groningen's CO-FUND Rosalind Franklin program. PD also acknowledges support from the European Research Council's starting grant ERC StG-717001 ('DELPHI') and from the NWO grant 016.VIDI.189.162 ('ODIN'). TRC acknowledges the support of the Department of Atomic Energy, Government of India, under project no. 12-R&D-TFR-5.02-0700.

## DATA AVAILABILITY

The data underlying this paper, which have been produced by our model, are available on reasonable request to the corresponding author.

## REFERENCES

* Amarantidis S. et al., 2019, MNRAS, 485, 2694

* Anglés-Alcázar D., Faucher-Giguère C.-A., Quataert E., Hopkins P. F., Feldmann R., Torrey P., Wetzel A., Kereš D., 2017, MNRAS, 472, L109

* Bandara K., Crampton D., Simard L., 2009, ApJ, 704, 1135

* Barausse E., 2012, MNRAS, 423, 2533

* Bañados E. et al., 2018, Nature, 553, 473

* Beckmann R. S., Devriendt J., Slyz A., 2019, MNRAS, 483, 3488

* Begelman M. C., Volonteri M., Rees M. J., 2006, MNRAS, 370, 289

* Blecha L. et al., 2016, MNRAS, 456, 961

* Boco L., Lapi A., Danese L., 2020, ApJ, 891, 94

* Bower R. G., Benson A. J., Malbon R., Helly J. C., Frenk C. S., Baugh C. M., Cole S., Lacey C. G., 2006, MNRAS, 370, 645

* Bower R. G., Schaye J., Frenk C. S., Theuns T., Schaller M., Crain R. A., McAlpine S., 2017, MNRAS, 465, 32

* Bromm V., Loeb A., 2003, ApJ, 596, 34

* Caplar N., Lilly S. J., Trakhtenbrot B., 2018, ApJ, 867, 148

* Cattaneo A., 2002, MNRAS, 333, 353

* Costa T., Sijacki D., Trenti M., Haehnelt M. G., 2014, MNRAS, 439, 2146

* Croton D. J., 2009, MNRAS, 394, 1109

* Croton D. J. et al., 2006, MNRAS, 365, 11

* Dayal P., Ferrara A., Dunlop J. S., Pacucci F., 2014, MNRAS, 445, 2545

* Dayal P., Rossi E. M., Shiralilou B., Piana O., Choudhury T. R., Volonteri M., 2019, MNRAS, 486, 2336

* Dayal P. et al., 2020, MNRAS, 495, 3065

* Decarli R., Falomo R., Treves A., Labita M., Kotilainen J. K., Scarpa R., 2010, MNRAS, 402, 2453

* Decarli R. et al., 2018, ApJ, 854, 97

* Di Matteo T., Croft R. A. C., Feng Y., Waters D., Wilkins S., 2017, MNRAS, 467, 4243
Dotti M., Volonteri M., Perego A., Colpi M., Ruszkowski M., Haardt F., 2010, MNRAS, 402, 682

Dubois Y., Volonteri M., Silk J., Devriendt J., Slyz A., Teyssier R., 2015, MNRAS, 452, 1502

Dubois Y. et al., 2014, MNRAS, 444, 1453

Duras F. et al., 2020, A&A, 636, A73

Eisenstein D. J., Loeb A., 1995, ApJ, 443, 11

Fan X. et al., 2001, AJ, 122, 2833

Favata M., Hughes S. A., Holz D. E., 2004, ApJ, 607, L5

Ferrarese L., 2002, ApJ, 578, 90

Fitchett M. J., 1983, MNRAS, 203, 1049

Griffin A. J., Lacey C. G., Gonzalez-Perez V., Lagos C. d. P., Baugh C. M., Fanidakis N., 2019, MNRAS, 487, 198

Habouzit M., Volonteri M., Dubois Y., 2017, MNRAS, 468, 3935

Habouzit M., Volonteri M., Latif M., Dubois Y., Peirani S., 2016, MNRAS, 463, 529

Habouzit M., Volonteri M., Somerville R. S., Dubois Y., Peirani S., Pichon C., Devriendt J., 2019, MNRAS, 489, 1206

Haehnelt M. G., Rees M. J., 1993, MNRAS, 263, 168

Haiman Z., 2004, ApJ, 613, 36

Haiman Z., Loeb A., 2001, ApJ, 552, 459

Hirschmann M., Somerville R. S., Naab T., Burkert A., 2012, MNRAS, 426, 237

Huang K.-W., Ni Y., Feng Y., Di Matteo T., 2020, MNRAS, 496, 1

Inayoshi K., Visbal E., Haiman Z., 2020, ARA&A, 58, 27

Kelly B. C., Shen Y., 2013, ApJ, 764, 45

Khandai N., Di Matteo T., Croft R., Wilkins S., Feng Y., Tucker E., DeGraff C., Liu M.-S., 2015, MNRAS, 450, 1349

Lacey C., Cole S., 1993, MNRAS, 262, 627

Lacey C. G. et al., 2016, MNRAS, 462, 3854

Latif M. A., Volonteri M., Wise J. H., 2018, MNRAS, 476, 5016

Lauer T. R., Tremaine S., Richstone D., Faber S. M., 2007, ApJ, 670, 249

Loeb A., Rasio F. A., 1994, ApJ, 432, 52

Luo Y., Ardaneh K., Shlosman I., Nagamine K., Wise J. H., Begelman M. C., 2018, MNRAS, 476, 3523

Lupi A., Haardt F., Dotti M., Fiacconi D., Mayer L., Madau P., 2016, MNRAS, 456, 2993

Lupi A., Volonteri M., Decarli R., Bovino S., Silk J., Bergeron J., 2019, MNRAS, 488, 4004

Marshall M. A., Mutch S. J., Qin Y., Poole G. B., Wyithe J. S. B., 2020, MNRAS, 494, 2747

Matsuoka Y. et al., 2018, PASJ, 70, S35

Mayer L., 2019, in Latif M., Schleicher D., eds, *Formation of the First Black Holes*. World Scientific Press, Singapore, p. 195

Mazzucchelli C. et al., 2017, ApJ, 849, 91

McLure R. J., Jarvis M. J., 2002, MNRAS, 337, 109

Merloni A., 2004, MNRAS, 353, 1035

Merloni A., Heinz S., 2008, MNRAS, 388, 1011

Merloni A. et al., 2010, ApJ, 708, 137

Mortlock D. J. et al., 2011, Nature, 474, 616

Mutch S. J., Geil P. M., Poole G. B., Angel P. W., Duffy A. R., Mesinger A., Wyithe J. S. B., 2016, MNRAS, 462, 250

Oke J. B., Gunn J. E., 1983, ApJ, 266, 713

Peng C. Y., 2007, ApJ, 671, 1098

Peng C. Y., Impey C. D., Rix H.-W., Kochanek C. S., Keeton C. R., Falco E. E., Lehár J., McLeod B. A., 2006, ApJ, 649, 616

Pezzulli E., Valiante R., Schneider R., 2016, MNRAS, 458, 3047

Pfister H., Volonteri M., Dubois Y., Dotti M., Colpi M., 2019, MNRAS, 486, 101

Planck Collaboration XIII, 2016, A&A, 594, A13

Qin Y. et al., 2017, MNRAS, 472, 2009

Reines A. E., Volonteri M., 2015, ApJ, 813, 82

Ricarte A., Natarajan P., 2018, MNRAS, 481, 3278

Richardson M. L. A., Scannapieco E., Devriendt J., Slyz A., Thacker R. J., Dubois Y., Wurster J., Silk J., 2016, ApJ, 825, 83

Rosas-Guevara Y., Bower R. G., Schaye J., McAlpine S., Dalla Vecchia C., Frenk C. S., Schaller M., Theuns T., 2016, MNRAS, 462, 190

Salpeter E. E., 1955, ApJ, 121, 161

Salviander S., Shields G. A., Bonning E. W., 2015, ApJ, 799, 173

Schaye J. et al., 2015, MNRAS, 446, 521

Sesana A., Barausse E., Dotti M., Rossi E. M., 2014, ApJ, 794, 104

Shankar F., Bernardi M., Haiman Z., 2009, ApJ, 694, 867

Shao Y. et al., 2017, ApJ, 845, 138

Shen Y. et al., 2015, ApJ, 805, 96

Sheth R. K., Tormen G., 2002, MNRAS, 329, 61

Shields G. A., Gebhardt K., Salviander S., Wills B. J., Xie B., Brotherton M. S., Yuan J., Dietrich M., 2003, ApJ, 583, 124

Sijacki D., Vogelsberger M., Genel S., Springel V., Torrey P., Snyder G. F., Nelson D., Hernquist L., 2015, MNRAS, 452, 575

Suh H., Civano F., Trakhtenbrot B., Shankar F., Hasinger G., Sanders D. B., Allevato V., 2020, ApJ, 889, 32

Tanaka T. L., 2014, Classical Quantum Gravity, 31, 244005

Thomas N., Davé R., Anglés-Alcázar D., Jarvis M., 2019, MNRAS, 487, 5764

Tremmel M., Governato F., Volonteri M., Quinn T., Pontzen A., 2018, MNRAS, 475, 4967

Ueda Y., Akiyama M., Hasinger G., Miyaji T., Watson M. G., 2014, ApJ, 786, 104

Valiante R., Schneider R., Volonteri M., Omukai K., 2016, MNRAS, 457, 3356

Venemans B. P., Walter F., Zschaechner L., Decarli R., De Rosa G., Findlay J. R., McMahon R. G., Sutherland W. J., 2016, ApJ, 816, 37

Vestergaard M., 2002, ApJ, 571, 733

Vogelsberger M. et al., 2014, MNRAS, 444, 1518

Volonteri M., Dubois Y., Pichon C., Devriendt J., 2016, MNRAS, 460, 2979

Volonteri M., Lodato G., Natarajan P., 2008, MNRAS, 383, 1079

Volonteri M., Rees M. J., 2005, ApJ, 633, 624

Volonteri M., Rees M. J., 2006, ApJ, 650, 669

Volonteri M., Reines A. E., 2016, ApJ, 820, L6

Volonteri M., Stark D. P., 2011, MNRAS, 417, 2085

Willott C. J. et al., 2010, AJ, 140, 546

Wise J. H., 2019, in Latif M., Schleicher D., eds, *Formation of the First Black Holes*. World Scientific Press, Singapore, p. 177

Wu X.-B. et al., 2015, Nature, 518, 512

Yoo J., Miralda-Escudé J., 2004, ApJ, 614, L25

This paper has been typeset from a TeX/LA TeX file prepared by the author.
Dotti M., Volonteri M., Perego A., Colpi M., Ruszkowski M., Haardt F., 2010, MNRAS, 402, 682

Dubois Y., Volonteri M., Silk J., Devriendt J., Slyz A., Teyssier R., 2015, MNRAS, 452, 1502

Dubois Y. et al., 2014, MNRAS, 444, 1453

Duras F. et al., 2020, A&A, 636, A73

Eisenstein D. J., Loeb A., 1995, ApJ, 443, 11

Fan X. et al., 2001, AJ, 122, 2833

Favata M., Hughes S. A., Holz D. E., 2004, ApJ, 607, L5

Ferrarese L., 2002, ApJ, 578, 90

Fitchett M. J., 1983, MNRAS, 203, 1049

Griffin A. J., Lacey C. G., Gonzalez-Perez V., Lagos C. d. P., Baugh C. M., Fanidakis N., 2019, MNRAS, 487, 198

Habouzit M., Volonteri M., Dubois Y., 2017, MNRAS, 468, 3935

Habouzit M., Volonteri M., Latif M., Dubois Y., Peirani S., 2016, MNRAS, 463, 529

Habouzit M., Volonteri M., Somerville R. S., Dubois Y., Peirani S., Pichon C., Devriendt J., 2019, MNRAS, 489, 1206

Haehnelt M. G., Rees M. J., 1993, MNRAS, 263, 168

Haiman Z., 2004, ApJ, 613, 36

Haiman Z., Loeb A., 2001, ApJ, 552, 459

Hirschmann M., Somerville R. S., Naab T., Burkert A., 2012, MNRAS, 426, 237

Huang K.-W., Ni Y., Feng Y., Di Matteo T., 2020, MNRAS, 496, 1

Inayoshi K., Visbal E., Haiman Z., 2020, ARA&A, 58, 27

Kelly B. C., Shen Y., 2013, ApJ, 764, 45

Khandai N., Di Matteo T., Croft R., Wilkins S., Feng Y., Tucker E., DeGraff C., Liu M.-S., 2015, MNRAS, 450, 1349

Lacey C., Cole S., 1993, MNRAS, 262, 627

Lacey C. G. et al., 2016, MNRAS, 462, 3854

Latif M. A., Volonteri M., Wise J. H., 2018, MNRAS, 476, 5016

Lauer T. R., Tremaine S., Richstone D., Faber S. M., 2007, ApJ, 670, 249

Loeb A., Rasio F. A., 1994, ApJ, 432, 52

Luo Y., Ardaneh K., Shlosman I., Nagamine K., Wise J. H., Begelman M. C., 2018, MNRAS, 476, 3523

Lupi A., Haardt F., Dotti M., Fiacconi D., Mayer L., Madau P., 2016, MNRAS, 456, 2993

Lupi A., Volonteri M., Decarli R., Bovino S., Silk J., Bergeron J., 2019, MNRAS, 488, 4004

Marshall M. A., Mutch S. J., Qin Y., Poole G. B., Wyithe J. S. B., 2020, MNRAS, 494, 2747

Matsuoka Y. et al., 2018, PASJ, 70, S35

Mayer L., 2019, in Latif M., Schleicher D., eds, *Formation of the First Black Holes*. World Scientific Press, Singapore, p. 195

Mazzucchelli C. et al., 2017, ApJ, 849, 91

McLure R. J., Jarvis M. J., 2002, MNRAS, 337, 109

Merloni A., 2004, MNRAS, 353, 1035

Merloni A., Heinz S., 2008, MNRAS, 388, 1011

Merloni A. et al., 2010, ApJ, 708, 137

Mortlock D. J. et al., 2011, Nature, 474, 616

Mutch S. J., Geil P. M., Poole G. B., Angel P. W., Duffy A. R., Mesinger A., Wyithe J. S. B., 2016, MNRAS, 462, 250

Oke J. B., Gunn J. E., 1983, ApJ, 266, 713

Peng C. Y., 2007, ApJ, 671, 1098

Peng C. Y., Impey C. D., Rix H.-W., Kochanek C. S., Keeton C. R., Falco E. E., Lehár J., McLeod B. A., 2006, ApJ, 649, 616

Pezzulli E., Valiante R., Schneider R., 2016, MNRAS, 458, 3047

Pfister H., Volonteri M., Dubois Y., Dotti M., Colpi M., 2019, MNRAS, 486, 101

Planck Collaboration XIII, 2016, A&A, 594, A13

Qin Y. et al., 2017, MNRAS, 472, 2009

Reines A. E., Volonteri M., 2015, ApJ, 813, 82

Ricarte A., Natarajan P., 2018, MNRAS, 481, 3278

Richardson M. L. A., Scannapieco E., Devriendt J., Slyz A., Thacker R. J., Dubois Y., Wurster J., Silk J., 2016, ApJ, 825, 83

Rosas-Guevara Y., Bower R. G., Schaye J., McAlpine S., Dalla Vecchia C., Frenk C. S., Schaller M., Theuns T., 2016, MNRAS, 462, 190

Salpeter E. E., 1955, ApJ, 121, 161

Salviander S., Shields G. A., Bonning E. W., 2015, ApJ, 799, 173

Schaye J. et al., 2015, MNRAS, 446, 521

Sesana A., Barausse E., Dotti M., Rossi E. M., 2014, ApJ, 794, 104

Shankar F., Bernardi M., Haiman Z., 2009, ApJ, 694, 867

Shao Y. et al., 2017, ApJ, 845, 138

Shen Y. et al., 2015, ApJ, 805, 96

Sheth R. K., Tormen G., 2002, MNRAS, 329, 61

Shields G. A., Gebhardt K., Salviander S., Wills B. J., Xie B., Brotherton M. S., Yuan J., Dietrich M., 2003, ApJ, 583, 124

Sijacki D., Vogelsberger M., Genel S., Springel V., Torrey P., Snyder G. F., Nelson D., Hernquist L., 2015, MNRAS, 452, 575

Suh H., Civano F., Trakhtenbrot B., Shankar F., Hasinger G., Sanders D. B., Allevato V., 2020, ApJ, 889, 32

Tanaka T. L., 2014, Classical Quantum Gravity, 31, 244005

Thomas N., Davé R., Anglés-Alcázar D., Jarvis M., 2019, MNRAS, 487, 5764

Tremmel M., Governato F., Volonteri M., Quinn T., Pontzen A., 2018, MNRAS, 475, 4967

Ueda Y., Akiyama M., Hasinger G., Miyaji T., Watson M. G., 2014, ApJ, 786, 104

Valiante R., Schneider R., Volonteri M., Omukai K., 2016, MNRAS, 457, 3356

Venemans B. P., Walter F., Zschaechner L., Decarli R., De Rosa G., Findlay J. R., McMahon R. G., Sutherland W. J., 2016, ApJ, 816, 37

Vestergaard M., 2002, ApJ, 571, 733

Vogelsberger M. et al., 2014, MNRAS, 444, 1518

Volonteri M., Dubois Y., Pichon C., Devriendt J., 2016, MNRAS, 460, 2979

Volonteri M., Lodato G., Natarajan P., 2008, MNRAS, 383, 1079

Volonteri M., Rees M. J., 2005, ApJ, 633, 624

Volonteri M., Rees M. J., 2006, ApJ, 650, 669

Volonteri M., Reines A. E., 2016, ApJ, 820, L6

Volonteri M., Stark D. P., 2011, MNRAS, 417, 2085

Willott C. J. et al., 2010, AJ, 140, 546

Wise J. H., 2019, in Latif M., Schleicher D., eds, *Formation of the First Black Holes*. World Scientific Press, Singapore, p. 177

Wu X.-B. et al., 2015, Nature, 518, 512

Yoo J., Miralda-Escudé J., 2004, ApJ, 614, L25

This paper has been typeset from a TeX/LATEX file prepared by the author.