The Black Hole Mass Function Across Cosmic Times II. Heavy Seeds and (Super)Massive Black Holes

This is the second paper in a series aimed at modeling the black hole (BH) mass function, from the stellar to the (super)massive regime. In the present work we focus on (super)massive BHs and provide an ab-initio computation of their mass function across cosmic times. We consider two main mechanisms to grow the central BH, that are expected to cooperate in the high-redshift star-forming progenitors of local massive galaxies. The first is the gaseous dynamical friction process, that can cause the migration toward the nuclear regions of stellar-mass BHs originated during the intense bursts of star formation in the gas-rich host progenitor galaxy, and the buildup of a central heavy BH seed $M_\bullet\sim 10^{3-5}\, M_\odot$ within short timescales $\lesssim$ some $10^7$ yr. The second mechanism is the standard Eddington-type gas disk accretion onto the heavy BH seed, through which the central BH can become (super)massive $M_\bullet\sim 10^{6-10}\, M_\odot$ within the typical star-formation duration $\lesssim 1$ Gyr of the host. We validate our semi-empirical approach by reproducing the observed redshift-dependent bolometric AGN luminosity functions and Eddington ratio distributions, and the relationship between the star-formation and the bolometric luminosity of the accreting central BH. We then derive the relic (super)massive BH mass function at different redshifts via a generalized continuity equation approach, and compare it with present observational estimates. Finally, we reconstruct the overall BH mass function from the stellar to the (super)massive regime, over more than ten orders of magnitudes in BH mass.


INTRODUCTION
The formation of (super)massive black holes (BHs) with masses M • ∼ 10 6−10 M and their role in galaxy evolution constitute crucial yet longstanding problems in modern astrophysics and cosmology. These monsters are thought to have grown mainly by gaseous accretion onto a disk surrounding the BH (e.g., Lynden-Bell 1969;Shakura & Sunyaev 1973) that energizes the spectacular broadband emission of active galactic nuclei (AGNs), and leaves a BH relic ubiquitously found at the center of massive galaxies in the local Universe (e.g., Kormendy & Ho 2013; also textbooks by Mo et al. 2010, Cimatti et al. 2020. This paradigm has recently received an astonishing confirmation by the EHT collaboration (2019,2022) via the imaging of the BH shadow caused by gravitational light bending and photon capture at the event horizon of M87 and Sgr A .
Accreting supermassive BHs can have a profound impact on the evolution of the host galaxies (see review by Alexander & Hickox 2012), as testified by the observed tight relationships between the relic BH masses and the physical properties of the hosts, most noticeably the stellar mass or velocity dispersion of the bulge component (e.g., Magorrian 1998;Ferrarese & Merritt 2000;Gebhardt et al. 2000;Tremaine et al. 2002;Kormendy & Ho 2013;McConnell & Ma 2013;Reines & Volonteri 2015;Sahu et al. 2019;Shankar et al. 2016Shankar et al. , 2020aZhu et al. 2021). These suggest that (apart from short-time stochastic fluctuations) the BH and the bulge stellar mass must have co-evolved over comparable timescales, possibly determined by the energy feedback from the BH itself on the gas/dust content of the host (see Tinsley 1980;Silk & Rees 1998;Fabian 1999;King 2005;Lapi et al. 2006Lapi et al. , 2014Lapi et al. , 2018; for a review, see King & Pounds 2015). In fact, targeted X-ray observations in the high-redshift star-forming progenitors of local massive galaxies have started to reveal the early growth of a dust-enshrouded (super)massive BH in their nuclear regions (e.g., Mullaney et al. 2012;Page et al. 2012;Delvecchio et al. 2015;Rodighiero et al. 2015Rodighiero et al. , 2019Fiore et al. 2017;Stanley et al. 2015Stanley et al. , 2017Massardi et al. 2018;Combes et al. 2019;D'Amato et al. 2020), before it attains a high enough mass and power to manifest as a bright AGN and to eventually reduce/quench star formation and partly evacuate gas and dust from the host (e.g., Granato et al. 2001Granato et al. , 2004Lapi et al. 2014Lapi et al. , 2018. Another, albeit more indirect, indication of coevolution for the bulk of the BH and the host stellar mass comes from the similarity between the activity timescales of central BH to the transition timescale of (green valley) galaxies from the blue cloud to the red sequence (see Wang et al. 2017;Lin et al. 2021Lin et al. , 2022 this is true apart from rejuvenations at late cosmic times, see Martin-Navarro et al. 2021).
However, two recent pieces of evidence may suggest that standard disk accretion is not the only process at work in growing a BH to the (super)massive regime. The first is the discovery of an increasing number of active BHs with masses M • 10 9 M at very high redshifts z 7 (e.g., Mortlock et al. 2011;Wu et al. 2015;Venemans et al. 2017;Reed et al. 2019;Banados et al. 2018Banados et al. , 2021Wang et al. 2019Wang et al. , 2021, when the age of the Universe was shorter than 0.8 Gyr. The second is the robust measurements of extreme BH masses M • 10 9−10 M at the center of early-type galaxies with stellar mass M 10 11 M (e.g., McConnell et al. 2011;Ferre-Mateu et al. 2015;Thomas et al. 2016;Mehrgan et al. 2019;Dullo et al. 2021), that have formed most of their old stellar component during a star-formation episode lasting some 10 8 yr at z 1, as demonstrated by astro-archeological measurements of their stellar ages and α-enhanced metal content (e.g., Thomas et al. 2005Thomas et al. , 2010Gallazzi et al. 2006Gallazzi et al. , 2014Johansson et al. 2012;Maiolino & Mannucci 2019;Morishita et al. 2019;Saracco et al. 2020). These observations concur to raise the issue of how billion solar mass BHs may have grown in less than a Gyr. In fact, this is somewhat challenging if standard disk accretion starts from a light seed ∼ 10 2 M of stellar origin and proceeds with the typical Eddington ratios λ 1 as estimated out to z ∼ 6 in active BHs (see Vestergaard & Osmer 2009;Nobuta et al. 2012;Kelly & Shen 2013;Dai et al. 2014;Kim & Im 2019;Duras et al. 2020;Ananna et al. 2022), that would require an overall time 0.8/λ Gyr to attain ∼ 10 9 M . Solutions may invoke mechanisms able to rapidly produce heavy BH seeds 10 3−5 M , so reducing the time required to attain the final masses by standard disk accretion (see Natarajan 2104, Mayer & Bonoli 2019and Volonteri et al. 2021 for exhaustive reviews). Viable possibilities comprise: direct collapse of gas clouds within a (proto)galaxy, possibly induced by galaxy mergers or enhanced matter inflow along cosmic filaments (e.g., Lodato & Natarajan 2007;Mayer et al. 2010Mayer et al. , 2015Di Matteo et al. 2012; merging of stars inside globular or nuclear star clusters (e.g., Portegies Zwart et al. 2004;Devecchi et al. 2012;Latif & Ferrara 2016); migration of stellar BHs towards the nuclear galaxy regions via dynamical friction against the dense gas-rich environment in strongly star-forming progenitors of local massive galaxies (e.g., Boco et al. 2020Boco et al. , 2021. 1 Such a complex picture for the overall (super)massive BH growth may in principle be probed via one of the most fundamental quantities for demographic studies of the BH population, namely the BH mass function, that expresses the number density of BHs per comoving volume and unit BH mass as a function of redshift. For (super)massive BHs, where most of the mass is accumulated through gas disk accretion, this is usually estimated (but still subject to systematics) from the AGN luminosity functions via Soltan (1982)-type or continuity equation arguments (e.g., Small & Blandford 1992;Haehnelt et al. 1998;Salucci et al. 1999;Yu & Tremaine 2002;Yu & Lu 2004Merloni & Heinz 2008;Cao 2010;Kelly & Merloni 2012;Aversa et al. 2015;Shankar et al. 2004Shankar et al. , 2009Shankar et al. , 2013, or from local galaxy mass/luminosity/velocity dispersion functions and scaling relation among these properties and the BH mass (e.g., Vika et al. 2009;Li et al. 2011;Mutlu-Pakdil et al. 2016;Shankar et al. 2016Shankar et al. , 2020a. In a future perspective, a precise assessment of the relic BH mass function is also important to work out detailed predictions for the gravitational wave emission expected from mergers of (super)massive BHs, that will constitute the primary targets of the upcoming Laser Interferometer Space Antenna mission (e.g., Sesana et al. 2016;Ricarte & 1 When referring to the dynamical friction mechanism, the term 'seeds' is used in a broader sense with respect to the classic meaning in the literature. A seed is usually referred to as the first compact object on which subsequent disk accretion occurs, eventually leading to the formation of a supermassive BH. The heavy seeds formed with the dynamical friction mechanism are by-products of multiple mergers of already-existing stellar-mass BHs (that in turn could be referred as light seeds) forming across a wide redshift range.
Natarajan 2018; for a review, Barausse & Lapi 2021 and references therein) and of ongoing and future Pulsar-Timing Array experiments (e.g., Antoniadis et al. 2022). Thus a theoretical grasp on the (super)massive BH mass function across cosmic times is of crucial importance. This is the second paper in a series aimed at modeling the BH mass function, from the stellar to the intermediate and (super)massive regime. In Sicilia et al. (2022; hereafter paper I) we have focused on the stellar mass BH relic mass function, while in the present work we provide an ab-initio computation of the redshift-dependent mass function for (super)massive BHs. We consider two mechanisms to grow the central BH, that likely cooperate in the high-redshift star-forming progenitors of local massive galaxies. The first one is the gaseous dynamical friction introduced by Boco et al. (2020), which can cause the migration of stellar-mass BHs originated during the intense bursts of star formation in the gas-rich central regions of the host progenitor galaxy, and the buildup of heavy BH seeds 10 5 M within short timescales some 10 7 yr. The second mechanism is the standard Eddington-type gas disk accretion onto the heavy seed, through which the central BH can become (super)massive within the typical star-formation timescales 1 Gyr of the host galaxy.
Our approach is semi-empirical, requires minimal modeling and a few educated assumptions, and is original in at least three respects: (i) we start from the galaxy SFR functions and derive BH-related statistics by jointly modeling the evolution of the central BH mass and the stellar mass of the host; (ii) we explicitly compute (and do not assume a priori) the heavy seed mass function by exploiting the distribution of stellar mass BHs originated from star-formation (Sicilia et al. 2022) and their migration rates due to dynamical friction ); (iii) we determine the detailed shape of the BH growth curve during disk accretion (in particular, we set the Eddington ratio) by requiring that the final BH and host stellar mass satisfy a Magorrian-like relationship, and that the star-formation timescale of the galaxy host is set by the main-sequence relation. We validate our approach by reproducing the observed redshiftdependent bolometric AGN luminosity functions and Eddington ratio distributions, and the relationship between the star-formation of the host galaxy and the bolometric luminosity of the accreting central BH. We then derive the relic (super)massive BH mass function at different redshifts via a generalized continuity equation approach, and compare it with present observational estimates. At the same time, we provide a robust theoretical basis for a physically-motivated heavy seed distribution at high redshifts. Finally, we put together the results from paper I and the present work to reconstruct the overall BH mass function from the stellar to the intermediate to the (super)massive regime, over more than ten orders of magnitudes in BH mass.
The plan of the paper is straightforward: in Section 2 we describe our semi-empirical framework, in Section 3 we present and discuss our results, and in Section 4 we summarize our main findings and outline future perspectives. In the Appendix we recall the basics of the gaseous dynamical friction mechanism (see Appendix A) and of the continuity equation technique (see Appendix B) exploited in the computations of the main text. Throughout this work, we adopt the standard flat ΛCDM cosmology (Planck Collaboration 2020) with rounded parameter values: matter density Ω M = 0.3, dark energy density Ω Λ = 0.7, baryon density Ω b = 0.05, Hubble constant H 0 = 100 h km s −1 Mpc −1 with h = 0.7, and mass variance σ 8 = 0.8 on a scale of 8 h −1 Mpc. A Kroupa (2001) initial mass function (IMF) in the star mass range m ∼ 0.1 − 150 M is adopted.

THEORETICAL BACKGROUND
We aim to derive the redshift-dependent (super)massive relic BH mass function dN/d log M • dV , i.e., the number density of massive BHs per unit comoving volume V and BH mass M • . We rely on two main mechanisms to grow the central BH mass, which are likely to cooperate in the gas-rich star-forming progenitors of local massive galaxies (hosting massive relic BHs): gas disk accretion and stellar BH migration via gaseous dynamical friction. Both processes will require the joint modeling of the stellar and BH mass growth history in a galaxy of given SFR and redshift.

Stellar mass growth
As to the stellar mass growth, we assume a simple two-stage star formation historẏ where ψ is the SFR, R ≈ 0.45 is the IMF-dependent gas fraction restituted to the ISM during stelllar evolution (with the quoted value applying for a Kroupa IMF), and Θ H [·] is the Heaviside step function. Basically, this reflects a constant SFR, that is then abruptly quenched by the radiative/kinetic power associated to SN explosions/stellar winds and/or the central BH activity at around the age τ b . This temporal evolution renders to a good approximation the behavior expected from state-of-the-art in-situ galaxy evolution models (e.g., Pantoni et al. 2019;Lapi et al. 2020), and is also indicated by SED-modeling studies of high-redshift dusty star-forming galaxies (e.g., Papovich et al. 2011;Smit et al. 2012;Moustakas et al. 2013;Steinhardt et al. 2014;Cassará et al. 2016;Citro et al. 2016;Gonzalez Delgado et al. 2017;Carnall et al. 2019;Williams et al. 2021;Pantoni et al. 2021), and by the observed fraction of IR-detected host galaxies in X-ray (e.g., Mullaney et al. 2012;Page et al. 2012;Rosario et al. 2012;Azadi et al. 2015;Stanley et al. 2015;Carraro et al. 2020) and in IR or optically selected AGNs (e.g., Mor et al. 2012;Wang et al. 2013;Willott et al. 2015;Stanley et al. 2017;Dai et al. 2018;Bianchini et al. 2019;Nguyen et al. 2020;Wang et al. 2021). Correspondingly, the stellar mass increases as and hereafter we will indicate for convenience M ,relic ≡ M (τ b ) = (1 − R) ψ τ b . Note, however, that this is the stellar mass just before the quenching at τ b , not the relic stellar mass at z ≈ 0; in fact, at late cosmic times this may be further increased by dry mergers, especially in very massive galaxies (e.g., Rodriguez-Gomez et al. 2016;Buitrago et al. 2017;Lapi et al. 2018). We can estimate the value of the star-formation duration τ b by requiring that, just before the quenching, the SFR ψ and the stellar mass M (τ b ) = M ,relic satisfy the redshift-dependent main sequence relationship ψ MS (M , z) (see Daddi et al. 2007;Rodighiero et al. 2011Rodighiero et al. , 2015Sargent et al. 2012;Speagle et al. 2014;Whitaker et al. 2014;Schreiber et al. 2015;Caputi et al. 2017;Bisigello et al. 2018;Boogaard et al. 2018;Leja et al. 2022;Rinaldi et al. 2022;Popesso et al. 2022); in other words, the condition sets the timescale τ b (ψ, z) for any galaxy with SFR ψ and redshift z. We adopt as our reference the main-sequence determination by Speagle et al. (2014) log ψ MS (M , z) M yr −1 ≈ (−6.51 + 0.11 t z ) + (0.84 − 0.026 t z ) log where t z is the age of the Universe at redshift z in units of Gyr. We will show in Sect. 3.2 the effect of adopting a different main sequence prescription.
2.2. BH growth due to dynamical friction Boco et al. (2020Boco et al. ( , 2021 have pointed out that the central BH can grow, especially in the early stages, by a continuous rain of stellar mass BHs that are funnelled toward the nuclear region via dynamical friction against the gas-rich background of high-redshift star-forming galaxies. The related growth rate is computed following Boco et al. (2020; see also Appendix A), to which we refer the interested reader for details. For consistency, in the present work we initialize the computation basing on the stellar BH mass function and light seed distribution derived in paper I, along the following lines.
First of all, we extract from the stellar and binary evolutionary code SEVN (see Spera et al. 2019) the so-called stellar term, i.e. the number of BHs originated per unit of stellar mass formed M SFR and BH mass m • : This includes three different contributions from isolated stars evolving into BHs ( → •), from stars that are originally in binary systems but end up as an isolated BH because one of the companions has been ejected or destroyed or cannibalized ( → •), and from stars in binary systems that evolve into binary BHs ( → ••). All these terms are strongly dependent on metallicity Z, which affects the efficiency of the various processes involved in stellar and binary evolution, like mass loss rates, mass transfers, core-collapse physics, etc. (see paper I for details).
We then derive the birthrate of stellar BHs with mass m • at time τ in an individual galaxy with SFR ψ at redshift z from the expression The integrand is the product of the stellar term from the previous Eq. (5) and of the metallicity distribution dp/d log Z.
For the latter we adopt a lognormal shape centered around the fundamental metallicity relation log Z FMR (ψ, M ) by Mannucci et al. (2011; for a review, see Maiolino & Mannucci 2019), with a dispersion of ∆ log Z FMR ≈ 0.15 dex, and M (τ ) given by Eq.
(2). We then compute the migration rate per unit stellar BH mass due to gaseous dynamical friction at time τ inside a galaxy with SFR ψ at redshift z here dp/dr and dp/dv r,θ are the probability distributions of initial radii and velocities and τ DF (m • , r, v r , v θ ) is the dynamical friction timescale against the gaseous background for a compact remnant of mass m • . All these quantities are recalled in Appendix A and detailed in Boco et al. (2020); in the latter paper the reader can find an account of how the dynamical friction timescale depends on such quantities and on the parameters ruling the gas distribution. Finally, the growth rate of the central BH due to the dynamical friction mechanism is jusṫ where the step function θ H (·) specifies that the mechanism is no longer active after τ b since the gaseous medium is expected to have been at least partly evacuated from the nuclear regions due to feedback processes.

BH growth due to gas accretion
In parallel, the central BH can grow due to standard gas disk accretion. We adopt a BH accretion rate curve with shape (e.g., Yu & Lu 2004Shen 2009;Li 2012;Lapi et al. 2014;Aversa et al. 2015) This describes a growth due to disk accretion in two stages, separated at the galaxy age τ b where star formation is quenched. The rationale behind the above expression is the following: at early epochs (ages τ τ b ) when there is plenty of material to accrete onto the BH in the nuclear galaxy regions, a demand-limited, Eddington-type BH accretion rate over a characteristic e−folding timescale τ ef is assumed. At late times (ages τ τ b ) the BH mass and radiative/kinetic power may be so large as to quench the star formation and partly evacuate gas from the host; however, if residual gas mass is still present in the central regions, it can be accreted in a supply-driven fashion, thus originating the exponentially declining part of the accretion curve with a characteristic timescale τ d . The IR-detected fraction of X-ray selected AGNs (seee Mullaney et al. 2012;Page et al. 2012;Rosario et al. 2012;Azadi et al. 2015;Stanley et al. 2015;Carraro et al. 2020) suggests τ d ≈ 2 τ ef , as shown by Lapi et al. (2014) and adopted by Aversa et al. (2015) and Mancuso et al. (2016bMancuso et al. ( , 2017. Eventually, we consider the accretion to stop for ages τ τ b + ζ τ d with ζ ≈ 3 (our results are anyway weakly affected by the value of this latter quantity); this is reasonable since at that point the accretion rate becomes so small withṀ •,acc τ ef /M • 10 −2 as to enter in an ADAF (i.e., advection-dominated accretion flow) regime, where the mass growth can be safely neglected with respect to that accumulated during the slim/thin disk accretion. We will show in Sect. 3.2 the effect of adopting a different, scale-free declining portion (e.g., Shen 2009) of the BH growth curve, that also avoids the inclusion of the quantity ζ.
Provided that L = Ṁ • c 2 is the accretion luminosity and λ ≡ L/L Edd the (luminous) Eddington ratio in terms of the Eddington luminosity 2 L Edd , the e−folding time of the early growth reads where t Edd = M • c 2 /L Edd ≈ 0.45 Gyr is the Eddington timescale and is the radiative efficiency. As to the latter, it is worth considering that in the early stages the demand-limited accretion may be prone to the development of a slim accretion disk (e.g., Abramowicz et al. 1988), while at late-times the supply-limited accretion tends to originate a classic thin-disk accretion (e.g., Shakura & Sunyaev 1973). To describe both conditions we express the radiative efficiency via the prescription by Aversa et al. (2015) valid for both thin and slim disks (see also Mineshige et al. 2000;Watarai et al. 2001;Li 2012;Madau et al. 2014): here thin is the efficiency during the thin-disk phase, which may range from ≈ 0.057 for a nonrotating BH to ≈ 0.32 for a maximally rotating Kerr BH (see Thorne 1974). We will adopt thin ≈ 0.15 as our fiducial value (see Davis & Laor 2011;Raimundo et al. 2012;Trakhtenbrot et al. 2017;Shankar et al. 2020b), but will show in Sect. 3.2 the effect of adopting a larger efficiency thin ≈ 0.3. Note that in principle thin may even depend on the galactic age since in the early stages the accretion is likely chaotic and so the spin of the BH should stay rather small, while in the late stages a coherent accretion is expected to set in and the spin can rapidly increase to maximal values (see Lapi et al. 2014); however, we neglect such spin/efficiency evolution in the present framework.

Overall BH growth
The overall growth of the central BH mass due to both dynamical friction and gaseous accretion writeṡ Given Eqs. (8) and (9), the previous equation can be formally integrated to yield the overall central BH mass growth where the value on the last line corresponds to the final, relic BH mass The dynamical friction process dominates in the initial growth stage for τ τ ef ; we will show that it provides, by inducing the migration of stellar BHs originated from star formation, heavy seeds of order 10 3−5 M within some 10 7 yr, before standard Eddington-type accretion takes over as the dominant mechanism for BH growth. Remarkably, our modelling above, at variance with other approaches in the literature, does not require assumptions regarding the seed BH from which to start gas accretion: the light seeds are provided by star formation and stellar evolution, and the heavy seeds by the gaseous dynamical friction mechanism, in a consistent way.
As a consequence, for any galaxy with SFR ψ at redshift z, the evolution of the BH mass is completely specified by assigning the Eddington factor λ of the early growth stage, which determines the radiative efficiency via Eq. (11) and hence the e−folding timescale τ ef . We empirically determine the Eddington ratio (see Shankar et al. 2020b) by requiring that the relic BH and stellar mass just after the quenching at τ b satisfy a Magorrian-like relation M •,Mag (M , z), with a possible redshift dependence. In other words, from the condition one can determine λ(ψ, z) for any galaxy with SFR ψ and redshift z. We rely on the debiased determination of the Magorrian relationship by Shankar et al. (2016Shankar et al. ( , 2020a log , but the latest studies suggest a mild evolution with η ≈ 0.2, that we take as our fiducial value (our results are anyway weakly affected by this choice). We will show in Sect. 3.2 the effect of adopting a different Magorrian relationship.
We stress that at least two low-redshift processes, that can in principle affect the BH mass function, have not been considered in our framework: (i) relic supermassive BHs can be reactivated by accretion of gas funnelled toward the central regions by galaxy mergers or internal disk instabilities (e.g., Di Matteo et al. 2005;Capelo et al. 2015), that can trigger spectacular radio-mode activity in terms of relativistic jets; (ii) relic supermassive BHs can coalesce following, with some delay, a galaxy merger. In fact, the impact of these processes on the supermassive BH mass function is still somewhat debated: an important role of galaxy mergers in reproducing the massive end of the mass function has been claimed in semi-analytic models (e.g., Marulli et al. 2008;Bonoli et al. 2009), while other semi-empirical and numerical approaches have instead pointed out a much more limited relevance of mergers on the BH mass function (e.g., Aversa et al. 2015;Steinborn et al. 2018;McAlpine et al. 2020). The detailed treatment of galaxy and BH mergers is beyond the main scope of the present paper, and is deferred to future work.

BH growth rate function
Toward a statistical description, we start from the SFR function dN/d log ψ dV , i.e., the number density of galaxies with given SFR ψ per unit comoving cosmological volume V at redshift z. For this we adopt the determination by Boco et al. (2021, their Fig. 1; for an analytic Schechter fit see Eq. 2 and Table 1 in Mancuso et al. 2016a) derived from an educated combination of the dust-corrected UV (e.g., Oesch et al. 2018;Bouwens et al. 2021), IR (e.g., Gruppioni et al. 2020;Zavala et al. 2021), and radio (e.g. Novak 2017; Ocran 2020) luminosity functions, appropriately converted into SFR (see Kennicutt & Evans 2012) using our assumed Kroupa (2001) IMF.
We first compute the central BH growth rate function whereṀ • (τ |ψ, z) is provided by Eq. (13) and dτ /d logṀ • is the related time spent by the BH in a logarithmic bin of given growth rate; the summation allows for multiple solutions τ i of the equationṀ • (τ |ψ, z) =Ṁ • , that are typically two for the growth curve assumed in this work. Note that the SFR dependence inṀ • (τ |ψ, z) is twofold: on the one hand it is related to the growth rate of heavy seeds by migration of stellar-mass BHs, whose birthrate ultimately depends on star formation; on the other hand, such a dependence is encoded in the Eddington ratio λ(ψ, z) derived after Eq. (14) and in the radiative efficiency given by Eq. (11). To allow for some scatter induced by the Magorrian-like relationship, one can write the first equality follows trivially from the properties of the Dirac δ D [·] function, while in the second we have substituted a log-logistic distribution with dispersionσ logṀ• ( √ 3/π) σ logṀ• in terms of the standard log-normal dispersion σ logṀ• . The reason for using a log-logistic distribution in place of the standard log-normal one is that having heavier tails it tends to maintain intrinsic power-law distributions at the high-mass end, as indicated by the data relating to the AGN luminosity functions and BH mass function (see discussion by Ren & Trenti 2021). Agreement with the latter statistics requires to adopt σ logṀ• ≈ 0.3 − 0.4 dex, in line with the scatter of the Magorrian.

AGN luminosity functions, Eddington ratios and mean SFRs
The broadband emission of AGNs is energized by the gas accretion onto the (super)massive BHs; thus a relevant statistics to validate our semi-empirical approach is the redshift-dependent bolometric AGNs luminosity function. This may be computed analogously to Eq. (16) as where the times τ i are now determined from the condition L AGN = Ṁ •,acc c 2 /(1 − ), with the gas accretion curvė M •,acc (τ |ψ, z) specified by Eq. (9). Notice that the Eddington ratio λ and the radiative efficiency here are not free parameters but are self-consistenly computed, for any SFR ψ and redshift z, by Eqs. (14) and (11). We will compare our results with the bolometric luminosity function determination by Shen et al. (2020), reconstructed from a large compilation of rest-frame B-band/UV (e.g., Hopkins et al. 2007;Giallongo et al. 2012;Manti et al. 2017;Kulkarni et al. 2018), soft/hard X-ray (e.g., Fiore et al. 2012;Ueda et al. 2014;Aird et al. 2015aAird et al. , 2015bMiyaji et al. 2015), and IR data (e.g., Assef et al. 2011;Lacy et al. 2015) collected in the past decades (see Shen et al. 2020 for details concerning bolometric and obscuration corrections). Notice that the integrand in Eq. (18) constitutes the number density of galaxies d 2 N/d log ψ d log L AGN dV per comoving volume in bins of SFR and AGN luminosity. Thus it may be exploited to build up the so-called coevolution plane SFR vs. L AGN , and the mean relationship between these two quantities. Finally, the previous expressions can also be adapted to derive the Eddington-ratio distribution by simply substituting L AGN with λ = L AGN /L Edd .

Relic BH mass function
To derive the relic (super)massive BH mass function dN/d log M • dV we exploit a generalized version of the continuity equation (see Yu & Lu 2004Aversa et al. 2015), whose derivation is recalled in Appendix B; this is basically a technique to relate the BH growth functions (or AGN luminosity functions) to the BH mass functions. The outcome reads Here the quantity dN/d logṀ • dV is the growth rate function from Eq. (16), all the integrand is computed at the maximum accretion rate for a given relic BH massṀ , and the various quantities implicitly entering there (e.g., λ, τ ef , τ d ) must be referred to a relic BH mass M • and redshift z. We stress that in our framework, at variance with many previous approaches based on continuity equation, the input AGN luminosity functions are not just taken from observations, but are derived from the galaxy statistics via Eq. (18). The related relic (super)massive BH mass density can be computed as where typically the integral is taken over BH masses M • 10 6 M . Fig. 1 summarizes in an illustrative way all the steps followed to compute the (super)massive BH mass function and described in this Section.

RESULTS AND DISCUSSION
In this Section we will show results of our empirical model concerning the growth of the central BH mass, AGN luminosity functions and Eddington ratio distribution, relationship of the AGN luminosity with the host SFR, and BH mass function. We will highlight the role played by the gaseous dynamical friction process in providing a physical mechanism to originate heavy seeds, so allowing the growth of the central BH to the supermassive regime at moderate Eddington ratios within the typical star-formation timescale of the host. We will also discuss the dependence of our basic results on various assumptions.

Basic results
To start with, in the top row of Fig. 2 we illustrate the time evolution of the central BH mass (left panels) and BH growth rate (right panels) in a prototypical star-forming galaxy with SFR ψ ∼ 300 M yr −1 at reference redshifts z ≈ 2 (top and middle rows) and z ≈ 6 (bottom row). In the top row the final BH mass is assumed to satisfy the average Magorrian relationship, while in the middle and bottom rows it is taken as a 3σ upper outlier with respect to the Magorrian; these latter instances are representative of extremely massive BHs, that are possibly sampled because of observational biases (especially at high redshifts). The overall growth is illustrated as black solid lines, and the corresponding Eddington ratio is reported in the first entry of the legend, while the contribution from migration of stellar BHs via gaseous dynamical friction is shown by the blue dot-dashed lines. It is seen that in all these cases the evolution of the total BH mass at small galactic ages is dominated by the growth due to migration of stellar BHs via gaseous dynamical friction; such a process can effectively build up a heavy central BH seed of mass M • ∼ 10 3−5 within 10 8 yr. Thereafter Eddington-type gas disk accretion takes over and can grow the central BH to the (super)massive regime M • 10 8−9 M . Remarkably, the overall effect of the early growth by dynamical friction is twofold. First, it allows the central BH to attain the final mass within a rather short timescale of some 10 8 yr; this can contribute to alleviate, or even to solve, the high-redshift quasar problem, i.e. the buildup of billion-solar-mass BHs in quasar hosts at z 6, when the age of the universe 1 Gyr constitutes a demanding constraint. Second, such a growth can be obtained with reasonable values of the Eddington ratios λ ∼ 0.3, that are in sound agreement with the observational determinations (see below); even in the extreme instance of an upper 3σ outlier of the Magorrian at z ≈ 6 (bottom panels), the growth can be achieved with sub-Eddington conditions λ 1.
In addition, in Fig. 2 we also illustrate what happens in the absence of the dynamical friction process, hence enforcing a BH growth by pure disk accretion. In particular, the orange dashed lines depict the evolution of a central BH with the same final mass as the solid lines but starting from a stellar mass seed ≈ 10 2 M ; such a case is seen to imply an appreciably higher Eddington ratio (reported in the last entry of the legend). In other words, growing the BH from light seeds of stellar origin to the supermassive regime would require a time 0.8/λ Gyr. Thus especially at high redshifts and/or for upper outliers of the Magorrian relationship (i.e., BHs with billion solar masses), the growth of the central BH should proceed at appreciably high values of λ, and possibly in super-Eddington conditions (as in the bottom panels). Though this instance can be partially justified theoretically (e.g., Li 2012;Madau et al. 2014) and there are hints of a few cases at z 6 (e.g., Fujimoto et al. 2022), it struggles somewhat against the bulk of present observational estimates at z 6 (see references below and Fig. 5). On the other hand, the red dotted lines refer to the evolution of a central BH growing to the same final mass and with the same Eddington ratio λ as the solid lines; such a case is seen to imply that the initial seeds must be 10 4 M . Therefore a specific mechanism, alternative to dynamical friction, must be in any case envisaged to obtain such massive seeds (e.g., Volonteri et al. 2021; see also Sect. 1).
In Fig. 3 we illustrate the growth rate function of the central BH at different redshifts z ∼ 1 − 8 (color-coded). As it can be seen from Eq. (16), its shape as a function ofṀ • is determined by a combination of galaxy statistics (i.e., the SFR functions) and the time spent by the central BH in a given bin of growth rate. The latter is in turn determined by the shape of the BH growth rate as a function of galactic age plotted in Fig. 2: at early times the BH rate grows almost linearly due to dynamical friction, at intermediate times it raises almost exponentially over the timescale τ ef due to disk accretion, and at late times it diminishes exponentially over the timescale τ d . The redshift evolution mirrors that of galaxy statistics, with the knee of the function first increasing toward largerṀ • out to z ≈ 2 and then receding at higher redshifts. The turnover of the function at z 1 at low accretion rates reflects the progressive inefficiency of the dynamical friction process (in turn mirroring the decreased efficiency of star-formation and stellar mass BH generation) toward late cosmic times.
In Fig. 4 we show the bolometric AGN luminosity functions at different redshifts z ≈ 1, 2, 4 and 6, computed from Eq. (18). The results from our approach are compared with the observational estimates collected by Shen et al. (2020; see full list of references therein) from selections in the rest-frame B-band/UV (e.g., Hopkins et al. 2007;Giallongo et al. 2012;Manti et al. 2017;Kulkarni et al. 2018), soft/hard X-ray (e.g., Fiore et al. 2012;Ueda et al. 2014;Aird et al. 2015aAird et al. , 2015bMiyaji et al. 2015), and IR (e.g., Assef et al. 2011;Lacy et al. 2015), converted using appropriate bolometric corrections (see Table 1 and Section 3 in Shen et al. 2020). The agreement is pretty good, both in terms of shape and redshift evolution. It is worth mentioning that the number density for AGNs with bright luminosities (especially toward high redshifts) may be overestimated in the data due to the uncertainties in the bolometric corrections. Note that we do not attempt a comparison with the observed AGN luminosity functions at z 1 since our framework does not include BH reactivations from late-time mergers and disk-instabilities (see Sect. 2); the latter are known to be a fundamental ingredient in determining the low-z AGN luminosity functions, especially at the faint end, though reproducing these observables has demonstrated to be a highly non-trivial task even for detailed models incorporating the aforementioned processes (see Griffin et al. 2019;Izquierdo-Villalba et al. 2020).
In Fig. 5 we illustrate the Eddington ratio distribution dN/dV d log λ and the average Eddington ratio log λ with its dispersion as a function of redshift. This is compared with observational estimates from different samples (see Duras et al. 2020;Kim & Im 2019;Vignali et al. 2018;Dai et al. 2014;Nobuta et al. 2012;Vestergaard & Osmer 2009). In our fiducial framework, the average Eddington ratio slowly increases from values λ ≈ 0.1 at z ≈ 1 to values λ ∼ 0.6 at z 4. The Eddington ratio distribution is quite broad, with a 1σ dispersion of 0.4 dex almost independent of the redshift. The outcome from our approach is in good agreement with the observational estimates, although the latter, being mainly based on single-epoch estimators, are still subject to considerable uncertainties, especially toward high redshift. Note that recently in the literature a lot of attention has been paid to observational estimates of the Eddington ratio distribution as a function of host galactic properties, most noticeably stellar mass and specific SFR (e.g., Bongiorno et al. 2016;Georgakakis et al. 2017;Aird et al. 2018Aird et al. , 2022Yang et al. 2019;Ananna et al. 2022;Carraro et al. 2020Carraro et al. , 2022; however, the estimates are still subject to considerable uncertainties, especially at z 0.5 and for massive galaxies. The comparison with, and the interpretation of such distributions is beyond the scope of the present paper, and we defer it to a future work. In the above Figs. 4 and 5 we also illustrate (dashed lines) the expected luminosity functions and average Eddington ratio when the dynamical friction mechanism is switched off and light BH seeds ≈ 10 2 M are assumed (a value taken as representative for the most massive seeds of stellar origin). The results on the luminosity functions are almost indistinguishable from our fiducial case, since by construction our approach imposes that, with or without dynamical friction, the final BH masses must adhere to the same Magorrian relationship and are obtained within the same timescales set by the main sequence; in turn, this implies that the peak AGN luminosities are very close to each other. However, without dynamical friction, this is at the cost of increasing somewhat the average Eddington ratio, because the growth starts from a lighter seed. Albeit in a statistical sense these higher values of λ are still within the large dispersion of the observational data, the problem may be exacerbated for the very massive BHs M • 10 9 M , and especially so at high z 6 whose formation would require λ ∼ a few (see Fig. 2 and related discussion above).
In Fig. 6 we illustrate the coevolution plane at a reference redshift z ≈ 2; this represents the number density of objects in the SFR ψ vs. AGN bolometric luminosity L AGN diagram (grey-scale color-coded); the average relationship and its 1 − 2σ scatter (solid line and shaded areas) are computed from such a distribution, taking into account the typical SFR detection threshold of present observations, around ψ ≈ 150 M yr −1 . The distribution of objects in the coevolution plane is again determined mainly by the number density of galaxies with a given value of the SFR, implying that galaxies with higher SFRs are rarer, and by the time a galaxy spends in different AGN luminosity bins. The average SFR and its scatter, computed taking into account the typical SFR detection threshold mentioned above, stays roughly constant with AGN luminosity out to L AGN ≈ 10 46 erg s −1 and then slowly increases. Such a rise occurs just because, statistically, to achieve a higher AGN luminosity, the BH must reside in a more massive galaxy with a higher initial SFR. For comparison, in Fig. 6 we report various observational determinations (see Page et al. 2012;Netzer et al. 2015;Stanley et al. 2015Stanley et al. , 2017Fan et al. 2016;Bianchini et al. 2019;Rodighiero et al. 2019) concerning different primary AGN selections in the optical, X-ray, IR or mixed (color-coded); detections are highlighted with full symbols and stacked data with open symbols. Our findings are remarkably consistent with observations, with the detections being distributed around the average relationship within its scatter, and the stacked measurements settling at the margin of the expected 2σ dispersion.
In Fig. 7 we illustrate the (super)massive relic BH mass function as derived from the continuity equation Eq. (19) at different redshifts z ≈ 0 − 8 (color-coded). The redshift evolution is quite strong down to z ≈ 2, with the knee (characteristic BH mass) strongly increasing from M • 10 7 M at z 8 up to M • 10 9 M for z 2; the evolution slows down considerably, especially at the high mass end, for z 2, such that essentially below z ≈ 1 the mass function undergoes only a minor evolution.
In Fig. 8 we show the related BH mass density computed after Eq. (20). It increases quite steeply from ρ • 10 3 M Mpc −3 at z 6 up to some ρ • 10 5 M Mpc −3 at z 1. The local BH mass density amounts to ρ • ≈ 6 × 10 5 M Mpc −3 , in sound agreement with the available observational determinations (see Shankar et al. 2004, 2009: Hopkins et al. 2007Marconi et al. 2004;Graham et al. 2007;Yu & Lu 2008). Fig. 8 also displays the contribution to the mass density from different BH mass ranges, to highlight that at z 6 and for M • 10 9 M , more massive BHs tend to accumulate their mass faster, displaying a kind of downsizing behavior.
In Fig. 9 we present the local BH mass function, and compare it with theoretical and observational estimates. In particular, the green shaded area refers to the uncertainty region in the current estimates of the BH mass function (see Shankar et al. 2016Shankar et al. , 2020a, obtained when combining the local stellar mass/velocity dispersion functions with various literature relationships linking BH mass to stellar mass/velocity dispersion of the host. We also report for comparison the classic estimates by Marconi et al. (2004; see also Shankar et al. 2009) via a simplified continuity equation approach, and by Vika et al. (2009) via an object-by-object analysis of the BH mass-host luminosity relationship. Our mass function is in agreement with most determinations for BH masses M • some 10 8 M . At the high-mass end it lies well within the Shankar et al. (2020a) uncertainty region, but it declines substantially slower with respect to the classic estimates by Marconi et al. (2004) and Vika et al. (2009).
We stress that to obtain a BH mass function with a steep behavior at the high mass end is a non-trivial task. Specifically, in our framework we determine λ from the empirical Magorrian relation and main sequence timescale, obtaining values λ < 1 that are in good agreement with the observed Eddington ratios; we also predict AGN luminosity functions closely matching the data. However, when inserted into the continuity equation these low λ values originate a rather flat BH mass function at the high mass end since large BH masses correspond to moderate peak AGN luminosities (approximately L AGN ∝ λ M • holds) falling in the rather flat portion of the luminosity function. Even the slightly higher λ values we obtain when switching off dynamical friction (see dashed line in Fig. 5) are not sufficient to appreciably steepen the BH mass function, which features a high mass end similar to our fiducial case. Contrariwise, in other literature studies (e.g., Aversa et al. 2015;Shen et al. 2020) a steep behavior of the mass function is enforced by starting from the observed AGN luminosity functions (not self-consistently predicting them, as in this work) and by assuming values λ 1 designed on purpose. For example, the redshift-dependent parameterization λ(z) ≈ 4 {1 − 0.5 × erfc[(z − 2)/3]} proposed by Aversa et al. (2015) works quite well in producing a steep BH mass function, but at the price of assuming λ values somewhat in tension with the observed average Eddington ratios (see dotted lines in Figs. 5 and 9). Insisting on such high λ values in a self-consistent approach while maintaining a good prediction of the AGN luminosity functions is still possible, but requires BH growth timescales 100 Myr, much shorter than derived via the main sequence prescription.

Robustness of results against main assumptions
In Fig. 10 we highlight the dependence of our results concerning the AGN luminosity function, redshift evolution of the average Eddington ratio, and local (super)massive BH mass function on various assumptions/relationships used in our reference framework.
First, we vary the main sequence relationship, switching from Speagle et al. (2014) to the recent determination by Popesso et al. (2022). In analogy with Eq. (4), this can be rendered as With respect to the almost linear relation by Speagle et al. (2014), the above is characterized by a steepening toward the lower stellar masses and a progressive flattening toward higher stellar masses, that have some relevance in galaxy formation since they may be interpreted as the effects of stellar feedback and mass quenching, respectively (e.g., Lapi et al. 2018;Daddi et al. 2022). Second, we vary the shape of the declining portion of the accretion rate curve in Eq. (9). In particular, we switch from an exponential to a scale-free, powerlaw shape M •,acc (τ ) =Ṁ •,acc ( Here ω > 1 rules the steepness of the decline, and we set ω ≈ 2.5 as in Shen (2009). Correspondingly, the overall BH growth at late times (cf. Eq. 12) follows 22) and the relic mass for τ . Such a powerlaw behavior is often adopted in empirical BH evolution models and generically ascribed to a residual accretion related to viscosity in a thin accretion disk (e.g., Yu & Lu 2008;Shen 2009); it has also been claimed to be consistent with a few numerical simulations present in the literature (see discussion by Habouzit et al. 2022).
Third, we vary the adopted Magorrian relationship (cf. Eq. 15) from the debiased determination by Shankar et al. ( , 2020a based on dynamical BH masses to that by Reines & Volonteri (2015) based on single-epoch virial estimators for locally active BHs (calibrated on a subsample of reverberation-mapped AGNs): where for consistency we retain the same redshift dependence adopted in Eq. (15). Finally, we vary the radiative efficiency thin of the thin disk regime (see Eq. 11) from our fiducial value 0.15 to 0.3; the latter is close to the limit applying for maximally spinning BHs. In fact, some theoretical works (e.g., Volonteri et al. 2013;Sesana et al. 2014;Griffin et al. 2019;Izquierdo-Villalba et al. 2020) have pointed out that the population of high-z BHs might be maximally spinning, so it is interesting to check the effect of this assumption especially on the AGN luminosity function at high redshift z 6. Fig. 10 shows that the most critical assumptions are, not surprisingly, the adopted main sequence and Magorrian relationships, that clearly affect the timescale of BH growth and the final BH masses, hence the resulting AGN luminosity function and BH mass function. As for the Popesso main sequence, it causes both a reduced number density of galaxies with high SFR, and a smaller stellar mass at a given SFR. This yields smaller BH masses hence a lower and steeper BH mass function. At the same time, with the Popesso main sequence shorter timescales are available for BH growth, implying minor variations in the Eddington ratio and correspondingly lower luminosities. As for the Magorrian relation by Reines & Volonteri (2015), it is flatter than our reference case and tends to yield lower BH masses for stellar masses M a few 10 10 M , and viceversa. Overall, this naturally originates an AGN luminosity function and a local BH mass function pumped at the faint end and depressed at the bright one, while the change in the average Eddington ratio is minor. Adopting a power-law shape of the declining portion in the BH accretion rate curve affects somewhat the AGN luminosity functions, while the impact on the Eddington ratio distribution and on the BH mass function is limited. Finally, we also highlight that adopting a high value thin ≈ 0.3 of the thin disk radiative efficiency implies higher Eddington ratio λ. This is seen by combining Eqs. (10) and (11) given that τ ef stays put since it is determined for a final BH mass by the Magorrian relation and the main sequence timescale. In the end this originates higher AGN luminosity functions, which better agree with observational estimates for z 6; this is particularly interesting since, as mentioned above, higher efficiencies associated to quickly spinning BHs are mostly expected toward such high redshifts.

The Overall BH mass function
In Fig. 11 we illustrate the overall black hole mass function, from the stellar to the (super)massive regime, over more than ten orders of magnitude in BH mass. The stellar regime for M • 10 2 M is taken from paper I and strictly associated to the star formation process in galaxies. In our framework, the intermediate mass regime M • ∼ 10 2−5 is mainly associated to the formation of heavy BH seeds by migration of stellar BHs via gaseous dynamical friction at the center of star-forming galaxies; the migrating stellar mass BHs are a very tiny fraction of the total number, so that the number density of these intermediate BHs is substantially lower than the stellar one. Finally, the (super)massive regime M • ∼ 10 6−10 M is mainly populated by the BHs that have grown to large masses (from heavy seeds) via Eddington-type gas disk accretion. Most of such massive BHs are active at high redshifts z 6, so that the BH mass function in the intermediate and (super)massive regime is continuously connected. On the other hand, moving toward lower redshifts the mass function in the (super)massive range increases because relic BHs grown by disk accretion accumulate, while the number of intermediate mass BHs diminishes since the dynamical friction process becomes less efficient and the overall production of stellar mass BHs also lowers (following the progressive decline in the amount of star-formation within galaxies). This is at the origin of the discontinuity (or gap; see Trinca et al. 2022 andSpinoso et al. 2022 for a similar behavior) between the intermediate and (super)massive mass function around M • 10 6 M ; it is pleasing that this transition occurs at around the typical value usually considered to separate intermediate from supermassive BHs.
For reference, in Fig. 11 we have also illustrated as colored boxes the mass and density ranges expected from other classic seed formation channels (taken from Volonteri et al. 2021, see their Fig. 1): remnants of the first massive pop-III stars (red box), direct collapse of primordial gas clouds (green box), and runaway stellar or BH mergers in compact primeval star clusters (yellow box). These distributions mainly originates in (proto)galaxies at z 10, and are then progressively eroded (but not substantially refurnished) at lower redshifts, when the seeds merge together or accrete gas and become more massive BHs (e.g., Mayer & Bonoli 2019;Volonteri et al. 2021;Trinca et al. 2022;Spinoso et al. 2022). This is at variance with our framework, where heavy seeds are continuously produced across cosmic times by the migration and merging of stellar-mass BHs associated to star formation in galaxies. In view of the above, if present, such classic seed formation channels are expected to enhance somewhat the BH mass function in the range M • ∼ 10 2−5 M especially at redshifts z 8. At later cosmic times, classic formation channels will feature a substantially eroded distribution in the intermediate mass range, so that their impact on our BH mass function should be minor. However, in a future work it would be interesting to perform a detailed investigation of the cooperative action of all these seed formation mechanisms across cosmic history.

SUMMARY AND OUTLOOK
In this work we have provided an ab-initio computation of the (super)massive BH mass function across cosmic times (see Fig. 1). To this purpose, we have started from the redshift-dependent galaxy statistics (constituted by the SFR functions) and have modeled the joint evolution of the central BH mass and the stellar mass of the host (see Fig.  2). We have considered two mechanisms to grow the central BH, that are expected to cooperate in the high-redshift star-forming progenitors of local massive galaxies. One is the gaseous dynamical friction envisaged by Boco et al. (2020), that can cause the migration of stellar-mass BHs originated during the intense bursts of star formation toward the gas-rich central regions of the host progenitor galaxies; this leads to the buildup of an heavy BH seeds 10 5 M within short timescales a few 10 7 yr. The second mechanism is the standard Eddington-type gas disk accretion onto the heavy seed, through which the central BH can become (super)massive within the typical star-formation timescales 1 Gyr of the host galaxy, as set by the galaxy main sequence. We have self-consistently combined these mechanisms to compute the overall growth rate functions of the central (super)massive BHs (see Fig. 3).
We have validated our approach by consistently reproducing the observed redshift-dependent bolometric AGN luminosity functions (Fig. 4), the observed Eddington ratio distributions (Fig. 5), and the observed relationship between the star-formation of the host galaxy and the bolometric luminosity of the accreting central BH (Fig. 6). We have then derived the relic (super)massive BH mass function (Fig. 7) and BH mass density (Fig. 8) via a generalized continuity equation approach, finding a pleasing agreement with the most recent observational estimates at z ≈ 0 (Fig. 9). All in all, we have found that the present (super)massive BH mass density amounts to ρ • ≈ 6 × 10 5 M Mpc −3 , in accord with available estimates.
We have stressed that in the absence of the dynamical friction process, statistical observables like the AGN luminosity functions are not substantially affected, since most of the BH mass is accumulated in the gas disk accretion phase. However, to attain BH masses 10 9 M within the typical star-formation duration 1 Gyr of the host galaxy without such dynamical friction process is challenging, especially at high redshifts z 6 or for overmassive BHs that are upper outliers of the average Magorrian relationship. In such a case, the BH growth must proceed at appreciably high Eddington ratios λ 1 and/or starting from heavy BH seeds M • ∼ 10 3−5 M . The first instance can be partially justified theoretically but struggles somewhat against present observational estimates, the second would require a specific mechanism, alternative to gaseous dynamical friction, designed to obtain such massive seeds.
Finally, putting together the results from paper I and the present work, we have reconstructed the overall BH mass function from the stellar to the (super)massive regime over more than ten orders of magnitude in BH mass. At the same time, we have provided a robust theoretical basis for a physically-motivated heavy seed distribution as a function of redshift. At variance with classic seed production channels, in our framework the heavy seed distribution is time-dependent: heavy seeds are continuously produced by the merging of light seeds originated from star formation via the gaseous dynamical friction mechanism; but they also grow via standard Eddington-type accretion, and soon leave the intermediate mass regime to become (super)massive. It would be extremely interesting to implement such a time-dependent seed distribution in analytic and numerical models of BH formation and evolution.
In a future perspective, our semi-empirical approach could be exploited to populate a N −body simulation, in order to build up a mock catalog encapsulating the three-dimensional spatial distribution and clustering of heavy seeds and of (super)massive BHs within their galactic hosts (e.g., Allevato et al. 2021). Another development could be a more detailed comparison of the properties of (super)massive BHs and host star forming galaxies, for example in terms of Eddington distributions as a function of BH environment and host galaxy properties (SFR, stellar mass, nuclear obscuration, etc.; e.g., Aird et al. 2018;Ananna et al. 2022). Moreover, we plan to work out predictions for upcoming or future observations via space instruments like JWST and Athena. Specifically, young BHs lying at the center or wandering in the nuclear regions of dusty starforming hosts may be detectable, even in the early stages of growth, via their X-ray and/or strongly extincted UV emissions; the latter could constitute a probe for the existence and abundance of intermediate mass BHs and could provide a characterization of their main growth mechanisms. Finally, we aim to exploit the BH mass function derived here to estimate the rate of (super)massive BH mergers. Although their effect on the overall mass function is expected to be mild and confined at the very massive end and late cosmic times, these events can constitute powerful sources of gravitational waves (e.g., Barausse & Lapi 2021). Thus we will provide detailed forecasts for their detectability by the Laser Interferometer Space Antenna mission and by present and future Pulsar-Timing Array experiments.

A. MIGRATION OF STELLAR BHs VIA GASEOUS DYNAMICAL FRICTION
In this Appendix we recall some details of the mechanism envisaged by Boco et al. (2020) to grow heavy seeds via migration of stellar mass BHs due to gaseous dynamical friction.
In the local Universe, supermassive BHs are hosted at the center of massive spheroidal galaxies. Thus their heavy seeds must have formed in the progenitors of such systems at intermediate/high redshifts, which are known to be dusty star-forming galaxies. These objects, detected and investigated mainly in the far-IR/(sub)mm band by groundbased interferometers like ALMA, feature large SFRs ψ 10 2 − 10 3 M yr −1 and huge molecular gas reservoirs M gas 10 10 − 10 11 M concentrated in a compact region of a few kpc (see Scoville et al. 2014Scoville et al. , 2016Ikarashi et al. 2015;Simpson et al. 2015;Barro et al. 2016;Spilker et al. 2016;Tadaki et al. 2017aTadaki et al. , 2017bTadaki et al. , 2018Lang et al. 2018;Talia et al. 2018Talia et al. , 2021Smail et al. 2021). These conditions are prompt for the efficient sinking of many compact objects toward the nuclear regions via gaseous dynamical friction (e.g., Ostriker 1999;Sanchez-Salcedo & Brandenburg 2001;Escala et al. 2004;Tanaka & Haiman 2009;Tagawa et al. 2016;Boco et al. 2020).
Specifically, Boco et al. (2020) have run a series of dynamical simulations and derived a fitting formula for the corresponding migration timescale of stellar mass BHs: here M gas is the total gas mass, R e is the half mass radius of the gas distribution, m • is the mass of the migrating compact object, and j are the initial specific energy and angular momentum of the compact object, r c ( ) is the circular radius that the compact object would have if it were on a circular orbit with energy , and j c ( ) is the angular momentum associated to that orbit. The precise values of the exponents (a, b, c, β, γ) and of the normalization factor N depend on the specific shape of the gas density profile. In the present work we adopt the fiducial setup of Boco et al. (2020), namely, a 3D Sersic gas density profile ρ(r) ∝ r −α e −k (r/Re) 1/n with n = 1.5, α = 1−1.188/2n+0.22/4n 2 ∼ 0.6 and half-mass radius R e ∼ 1 kpc. Then the values for the parameters in Eq. (A1) read a ≈ −0.95, b ≈ 0.45, c ≈ −1.2, β ≈ 1.5, γ ≈ 2.5 and N ≈ 3.4 × 10 8 yr. The effect of different setups on the dynamical friction timescale is discussed in Boco et al. (2020). Given the high SFR ongoing in the progenitors of local spheroidal galaxies, a lot of stars and compact remnants are formed in a short timescale within the nuclear regions. We assume that stars are initially distributed in space as the gas density profile ρ, so that the probability distribution for a star to be born at distance r from the galactic center is dp/dr ∝ r 2 ρ(r). After 10 7 yr massive stars (m 7−8 M ) undergo a supernova explosion possibly leaving a stellar mass BH. We assume that the latter follow the same velocity distribution of the progenitor stars, which is in turn related to that of the star-forming molecular gas cloud. In particular, we assume a Gaussian distributions of radial and tangential velocities: dp/dv r,θ ∝ e −v 2 r,θ /2σ 2 with dispersion σ(r) found by solving the isotropic Jeans equation: σ 2 (r) ∝ ρ −1 (r) ∞ r dr ρ(r ) r −2 r 0 dr r 2 ρ(r ). From these distributions the initial positions and velocities of stellar BHs, their initial energy and angular momentum can be easily extracted.
For a galaxy with spatially-averaged SFR ψ, we compute the associated relic stellar mass M from the well-established galaxy main sequence relationships (e.g., Speagle et al. 2014) and then estimate the initial gas mass M gas , entering the dynamical friction timescale, via the redshift-dependent M gas − M relation from abundance matching techniques (see Moster et al. 2013Moster et al. , 2018Aversa et al. 2015;Shi et al. 2017;Behroozi et al. 2019). Finally, the dynamical friction timescale τ (m • , r, v R , v θ ) can be computed from Eq. (A1), and the convolution of the stellar BH birthrate with the aforementioned distributions of initial position and velocity yields the migration rates according to Eq. (7).

B. CONTINUITY EQUATION
In this Appendix we provide some details on how to solve the continuity equation in the integral formulation, along the lines envisaged by Yu & Lu (2004 and Aversa et al. (2015). The continuity equation links the BH mass function N (M • , t) ≡ dN/dM • dV and the BH growth rate function N (Ṁ • , t) ≡ dN/dṀ • dV according to where M • (Ṁ • , t) is the minimum BH mass that has accreted atṀ • . After differentiating both sides by logṀ • one obtains Rearranging the expression and integrating over cosmic time, one finally can write the solution of the continuity equation in the form where all the integrand is calculated at the maximum growth rateṀ • (M • , t) for a given BH mass. This is a very general solution of the continuity equation, that holds even when the parameters defining the growth curve, e.g. the e−folding time, depend on M • ,Ṁ • and cosmic time t. As a simple application, consider the very special case when the growth of the BH occurs by gas accretion in an Eddington-limited regime (with constant e−folding time τ ef independent of accretion rate and cosmic time) up to a time τ b , so thatṀ In such an instance one has that the maximum BH accretion rate attained by a BH with final mass M • iṡ M • (M • ) = M • /τ ef , hence ∂ logṀ• log M • = 1. Moreover, i dτ i /d logṀ • = τ ef ln(10). All in all, the solution writes which is the classic expression derived by Marconi et al. (2004) and Shankar et al. (2004).
We thank the referee for a competent and constructive report. We acknowledge S. Bonoli Figure 1. Schematics showing the main steps to compute the (super)massive relic BH mass function. The starting point is the stellar term from paper I, representing the number of BHs originated per unit of star formed mass, and includes contributions from the evolution of isolated or binary stars into isolated or binary BHs (light seeds; see Eq. 5). This is coupled to the metallicity distribution (extracted from the fundamental metallicity relation) and with the timescale for gaseous dynamical friction by Boco et al. (2020; see Eq. A1) to derive the growth rate of the central BH by migration of stellar remnants (see Eqs. 7 and 8). In parallel, galaxy statistics provided by the SFR functions are coupled with model growth curves of the stellar and BH mass (see Eqs. 2 and 13); the latter includes the growth by dynamical friction migration and by gaseous Eddington-type accretion. Crucial parameters of these growth curves, like the star formation duration and the Eddington factor, are derived by requiring consistency with the main sequence of star-forming galaxies (see Eqs. 3 and 4) and with the local Magorrian relationship (see Eqs. 14 and 15). The main outcome of this procedure are BH growth rate functions (see Eq. 16), and byproducts are Eddington ratio functions, AGN luminosity functions (Eq. 18) and the coevolution plane SFR vs. LAGN. Finally, a generalized continuity equation approach allows to convert the growth rate functions into the (super)massive BH mass function (Eq. 19). Figure 2. Time evolution of the central BH mass (left panels) and BH growth rate (right panels) in a star-forming galaxy with SFR ψ ∼ 300 M yr −1 at z ≈ 2 (top and middle rows) and z ≈ 6 (bottom row); in the top row the final BH mass is on the average Magorrian relationship, in the middle and bottom rows it is a 3σ upper outlier of the Magorrian relationship. The overall growth of the central (super)massive BH is illustrated by black solid line (and the corresponding Eddington ratio λ is indicated in the first entry of the legend) while the contribution from migration of stellar BHs via gaseous dynamical friction is shown by the blue dot-dashed line. The red dotted line represents the evolution of a central BH growing by pure disk accretion (i.e., without dynamical friction) with the same final mass and with the same λ as the solid line, implying that the initial seed must be 10 4 M . Finally, the orange dashed line shows the evolution of a central BH growing by pure disk accretion (i.e., without dynamical friction) with the same final mass as the solid line from a stellar mass seed ≈ 100 M , implying an appreciably higher Eddington ratios (indicated in the last entry of the legend). . The Eddington ratio distribution and average Eddington ratio as a function of redshift z. The intensity of the black and white background illustrates the Eddington ratio distribution, while the black solid line is the average relationship expected from our approach (dark and light grey shades represents the 1σ and 2σ dispersion). In addition, dashed line is the average Eddington ratio when the gaseous dynamical friction mechanism is switched off, and the dotted line is the average Eddington ratio adopted on an empirical basis by Aversa et al. (2015).   Rodighiero et al. (2019;pentagons). Symbol colors refer to observational selection in the optical (blue), X-ray (orange), IR (green), or mixed (magenta); moreover, filled symbols refer to detections, while empty symbols refer to stacking estimates.  . The (super)massive BH mass function at z ≈ 0. Solid line illustrates the outcome of our framework, while dotted line is the mass function originated when coupling the observed AGN luminosity functions with the average Eddington ratio adopted by Aversa et al. (2015). Observational estimates are from Marconi et al. (2004; blue circles), Vika et al. (2009;red diamonds), and Shankar et al. (2009Shankar et al. ( , 2016Shankar et al. ( , 2020; green shaded area); the latter reflects the overall uncertainty region when determining the BH mass function from the local stellar mass/velocity dispersion functions combined with various relationships of these observables with the BH mass. Figure 10. Dependence of our results concerning AGN luminosity function at z ∼ 2 (top left) and z ∼ 6 (top right), average Eddington ratio as a function of redshift (middle), and local supermassive BH mass function (bottom) to various assumptions/relationships employed in this work. In all panels solid lines refer to our fiducial assumptions, dashed lines to our results when the main sequence relation by Popesso et al. (2022) is used in place of Speagle et al. (2014), dot-dashed line to our results when the BH accretion rate curves is characterized by a powerlaw decline instead of an exponential one, dotted lines to our results when the Magorrian relation by Reines & Volonteri (2015) for AGNs is employed in place of the one by Shankar et al. (2016Shankar et al. ( , 2020a, and triple-dot-dashed line to our results when the thin disk efficiency thin ≈ 0.3 is adopted instead of our fiducial value thin ≈ 0.15. See Sect. 3.2 for details. Figure 11. The overall BH mass function from our semi-empirical framework, from the stellar to the intermediate to the (super)massive regime, at different redshifts z ≈ 0 (cyan), 1 (orange), 2 (green), 4 (red), 6 (purple), and 8 (brown). The colored boxes illustrate the mass and density ranges from other seed formation channels (see Volonteri et al. 2021): remnants of the first massive pop-III stars (red box), direct collapse of primordial gas clouds (green box), and runaway stellar or BH mergers in compact primeval star clusters (yellow box).