The Black Hole Mass Function Across Cosmic Times. I. Stellar Black Holes and Light Seed Distribution

Alex Sicilia; Andrea Lapi; Lumen Boco; Mario Spera; Ugo N. Di Carlo; Michela Mapelli; Francesco Shankar; David M. Alexander; Alessandro Bressan; Luigi Danese

doi:10.3847/1538-4357/ac34fb

1. Introduction

The formation and evolution of black holes (BHs) in the universe is one of the major issues to be addressed by the modern research in astrophysics and cosmology. In the mass range m_• ∼ 5–150 M_⊙, BHs are originated from the final, often dramatic stages in the evolution of massive stars (possibly hosted in binary systems). These compact remnants can produce luminous X-ray binaries (e.g., Mapelli et al. 2010; Farr et al. 2011; Inoue et al. 2016), can constitute powerful sources of gravitational waves for ground-based detectors like the current LIGO/Virgo facility (e.g., Belczynski et al. 2010; Dominik et al. 2015; Spera & Mapelli 2017; Boco et al. 2019; Spera et al. 2019; Abbott et al. 2021a, 2021b), can possibly energize short gamma-ray bursts and associated kilonovas (e.g., Abbott et al. 2020, 2021c; Ackley et al. 2020; Gompertz et al. 2020), can inject strong energy inputs in the primeval universe (e.g., Mirabel et al. 2011; Justham & Schawinski 2012; Artale et al. 2015; Madau & Fragos 2017; Lehmer et al. 2021), and can provide light seeds for the subsequent growth of more massive BHs (e.g., Madau et al. 2014; Volonteri et al. 2015; Lupi et al. 2016; Pacucci et al. 2017; Boco et al. 2020; Das et al. 2021). At the other end, in the range M_• ∼ 10⁶–10¹⁰ M_⊙, supermassive BHs grow mainly by gaseous accretion that energizes the spectacular broadband emission of active galactic nuclei (AGNs). Such an activity can have a profound impact on galaxy evolution (e.g., Alexander & Hickox 2012; Lapi et al. 2014, 2018), as testified by the strict relationships between the relic BH masses and the physical properties of the hosts (e.g., Kormendy & Ho 2013; Shankar et al. 2016, 2020; Zhu et al. 2021). The intermediate-mass range m_• ∼ 10³–10⁶ M_⊙ is the most uncertain. So far, only tentative evidence of these systems has been identified (see Paynter et al. 2021). However, the chase is open in view of their astrophysical relevance. Most noticeably, they can provide heavy seeds for quick (super)massive BHs growth (e.g., Mayer & Bonoli 2019; Boco et al. 2020), as it seems required by the puzzling observations of an increasing numbers of giant monsters M_• ≳ 10⁹ M_⊙ when the age of the universe was shorter than ≲0.8 Gyr (e.g., Mortlock et al. 2011; Venemans et al. 2017; Banados et al. 2018). Moreover, such intermediate-mass BHs will constitute important targets for space-based gravitational wave detectors like LISA and DECIGO (see eLisa Consortium 2013; Kawamura et al. 2021; also Barausse & Lapi 2021; Boco et al. 2021b).

One of the most fundamental quantities for demographic studies of the BH population is constituted by the relic mass function, namely the number density of BHs per comoving volume and unit BH mass, as a function of redshift. In the supermassive regime, where most of the BH mass is accumulated through gas accretion, this is usually determined from the AGN luminosity function via Soltan (1982)-type or continuity equation arguments (e.g., Small & Blandford 1992; Salucci et al. 1999; Shankar et al. 2004, 2009, 2013, 2020; Yu & Lu 2004; Merloni & Heinz 2008; Kelly & Merloni 2012; Aversa et al. 2015), or from galaxy statistics and scaling relations among galactic and BH properties (e.g., Vika et al. 2009; Li et al. 2011; Mutlu-Pakdil et al. 2016). At the other end, the stellar BH mass function is largely unknown, and the most promising messenger to constrain it is constituted by the gravitational wave emission from coalescing binary BH systems (e.g., Kovetz et al. 2017; Perna et al. 2019; Ding et al. 2020; Tang et al. 2020); a first determination of the primary mass distribution for merging binary BHs has been established by the LIGO/Virgo collaboration (see Abbott et al. 2021a), although it depends somewhat on some assumptions (see Baxter et al. 2021). Therefore a theoretical grasp on the stellar BH mass function at different redshifts is of crucial importance.

In the present work, we provide an ab initio computation of the stellar BH relic mass function across cosmic times, by coupling the state-of-the-art stellar and binary evolutionary code SEVN to redshift-dependent galaxy statistics and empirical scaling relations involving metallicity, star formation rate and stellar mass. Note that here we mainly consider the standard, and likely dominant production channel of stellar-mass BHs constituted by isolated single/binary star evolution, and defer to a future paper the treatment of other formation mechanisms like dynamical effects in dense star clusters (e.g., Miller & Hamilton 2002; Antonini et al. 2019), AGN disks (e.g., McKernan et al. 2012; Yang et al. 2019), hierarchical triples (e.g., Kimpson et al. 2016; Fragione et al. 2019) that are thought to produce corrections in the high-mass tail of the mass function, toward the intermediate-mass BH regime.

The plan of the paper is straightforward. In Section 2, we introduce the theoretical background underlying our computation of the stellar BH relic mass function; specifically, we highlight the role of quantities related to galaxy formation (Section 2.1) and to stellar evolution (Section 2.2), and show how to include the effects of binary BH mergers (Section 2.3). In Section 3, we present our results concerning the stellar BH cosmic birthrate, the redshift-dependent stellar BH relic mass function and relic mass density; we also compare the primary mass distribution for merging BH binaries with the recent estimates from gravitational wave observations by LIGO/Virgo. In Section 4, we critically discuss how our results depend on the adopted binary stellar evolution code. Moreover, we estimate the effects of binary BH formation in dense environments like (young) star clusters. We also highlight the relevance of our computations in providing a light seed distributions for BH growth at high redshift. Finally, in Section 5, we summarize our findings and outlook future developments.

Throughout this work, we adopt the standard flat ΛCDM cosmology (Planck Collaboration et al. 2020) with rounded parameter values: matter density Ω_M = 0.3, dark energy density Ω_Λ = 0.7, baryon density Ω_b = 0.05, Hubble constant H₀ = 100 h km s⁻¹ Mpc⁻¹ with h = 0.7. A Kroupa (2001) initial mass function (IMF) in the star mass range m_⋆ ∼ 0.1–150 M_⊙ and a value Z_⊙ ≈ 0.015 for the solar metallicity are adopted (see Caffau et al. 2011).

2. Theoretical Background

We aim at deriving the stellar BH relic mass function, i.e., the number density of stellar BHs per unit comoving volume V and remnant mass m_• accumulated down to redshift z:

$\begin{eqnarray}&&\displaystyle \frac{{{dN}}_{\bullet }}{{dVd}\mathrm{log}\,{m}_{\bullet }}({m}_{\bullet }| \gt z)={\int }_{z}^{\infty }{dz}^{\prime} \,\displaystyle \frac{{{dt}}_{z^{\prime} }}{{dz}^{\prime} }\,\displaystyle \frac{d{\dot{N}}_{\bullet }}{{dV}\,d\mathrm{log}\,{m}_{\bullet }}({m}_{\bullet }| z^{\prime} ),\end{eqnarray} \tag{ 1 }$

and the related stellar BH relic mass density:

$\begin{eqnarray}&&{\rho }_{\bullet }(z)=\int d\mathrm{log}\,{m}_{\bullet }{m}_{\bullet }\displaystyle \frac{{{dN}}_{\bullet }}{{dVd}\mathrm{log}\,{m}_{\bullet }}({m}_{\bullet }| z).\end{eqnarray} \tag{ 2 }$

In Equation (1), the quantity t_z is the cosmic time corresponding to redshift z and $d{\dot{N}}_{\bullet }/{dV}\,d\mathrm{log}\,{m}_{\bullet }$ is the stellar BH cosmic birthrate, i.e., the production rate of BH of given mass per unit comoving volume. In turn, the latter can be expressed as

$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{d{\dot{N}}_{\bullet }}{{dVd}\mathrm{log}\,{m}_{\bullet }}({m}_{\bullet }| z)=\displaystyle \int d\mathrm{log}Z\\ \quad \times \,\displaystyle \frac{{{dN}}_{\bullet }}{{{dM}}_{\mathrm{SFR}}d\mathrm{log}\,{m}_{\bullet }}({m}_{\bullet }| Z)\displaystyle \frac{d{\dot{M}}_{\mathrm{SFR}}}{{dVd}\mathrm{log}Z}(Z| z),\end{array}\end{eqnarray} \tag{ 3 }$

where Z is the metallicity and M_SFR is the star-formed mass. The first and second factors in the integrand are usually referred to as the stellar and the galactic term, respectively (see Chruslinska & Nelemans 2019; Boco et al. 2021a); we will now describe each of them in some detail, starting with the latter. The main steps toward the computation of the stellar BH relic mass function are sketched in Figure 1.

Figure 1. Refer to the following caption and surrounding text. — **Figure 1.** Schematics showing the main steps to compute the stellar BH relic mass function (Equation (1)). This is obtained by integration over redshift of the BH cosmic birthrate (Equation (3)), which is in turn determined via the convolution of the galactic and the stellar terms. The galactic term (Equation (4)) is computed by convolving the galaxy statistics (based on the SFR function) with the metallicity distribution of a galaxy at given SFR and stellar mass (assumed as a log-normal shape around the fundamental metallicity relation; see Equation (6)) and with the stellar-mass distribution of a galaxy at given SFR and redshift (built from the redshift-dependent main sequence of star-forming galaxies; see Equation (5)). The stellar term (Equation (7)) is computed from the stellar and binary evolutionary code `SEVN` by summing up, for a given IMF, the probabilities per unit star-formed mass that (i) a single star evolves into a BH remnant (single stellar evolution); (ii) a binary stellar system evolves into a single BH remnant or two BHs no longer bounded (failed binaries); (iii) a binary stellar system evolves into a binary BH (binaries) that may eventually coealesce into a single BH via emission of gravitational waves.
Download figure:
Standard image High-resolution image

2.1. The Galactic Term

The galactic term in Equation (3) represents the cosmic SFR density per unit cosmic volume and metallicity; in other words, it constitutes the classic "Madau" plot (see Madau & Dickinson 2014 and references therein) sliced in metallicity bins, and as such is mainly related to galaxy formation and evolution. It may be estimated in various ways, starting from different galaxy statistics and empirical scaling laws. We compute it as

$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{d{\dot{M}}_{\mathrm{SFR}}}{{dV}\,d\mathrm{log}Z}(Z| z)=\displaystyle \int d\mathrm{log}\psi \,\psi \,\displaystyle \frac{{dN}}{{dV}\,d\mathrm{log}\psi }(\psi ,z)\\ \quad \times \displaystyle \int d\mathrm{log}{M}_{\star }\,\displaystyle \frac{{dp}}{d\mathrm{log}{M}_{\star }}({M}_{\star }| \psi ,z)\displaystyle \frac{{dp}}{d\mathrm{log}Z}(Z| {M}_{\star },\psi ),\end{array}\end{eqnarray} \tag{ 4 }$

where ψ is the SFR and M_⋆ the stellar mass of a galaxy. Three ingredients enter into the above expression. The first is constituted by the redshift-dependent SFR functions ${dN}/{dV}\,d\mathrm{log}\psi$ , expressing the number of galaxies per unit cosmological volume and SFR bin. For these, we adopt the determination by Boco et al. (2021a, their Figure 1; for an analytic Schechter fit see Equation (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 et al. 2017; Ocran 2021) luminosity functions, appropriately converted in SFR (see Kennicutt & Evans 2012).

The second ingredient is the probability distribution of stellar mass at given SFR and redshift:

$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{dp}}{d\mathrm{log}{M}_{\star }}({M}_{\star }| \psi ,z)\\ \quad \propto \,\left\{\begin{array}{ll}{M}_{\star } & {M}_{\star }\lt {M}_{\star ,\mathrm{MS}}(\psi ,z)\\ {M}_{\star ,\mathrm{MS}}\,\exp \left\{-\tfrac{{\left[\mathrm{log}{M}_{\star }-\mathrm{log}{M}_{\star ,\mathrm{MS}}\left(\psi ,z\right)\right]}^{2}}{2\,{\sigma }_{\mathrm{log}{M}_{\star }}^{2}}\right\} & {M}_{\star }\geqslant {M}_{\star ,\mathrm{MS}}(\psi ,z)\end{array}\right.,\end{array}\end{eqnarray} \tag{ 5 }$

where M_⋆,MS(ψ, z) is the observed redshift-dependent galaxy main sequence with log-normal scatter ${\sigma }_{\mathrm{log}{M}_{\star }}\approx 0.2$ dex (we adopt the determination by Speagle et al. 2014 for an anlytic fit, see their Equation (28)). The main sequence is a relationship between SFR and stellar mass followed by the majority of star-forming galaxies, apart from some outliers located above the average SFR at given stellar mass (see Daddi et al. 2007; Rodighiero et al. 2011, 2015; Sargent 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). The expression in Equation (5) holds for an approximately constant SFR history, which is indicated both by in situ galaxy formation scenarios (see Mancuso et al. 2016b; Pantoni et al. 2019; Lapi et al. 2020) and by observations of ETG progenitors (that have on overage slowly rising star formation history with typical duration of ≲1 Gyr; see Papovich et al. 2011; Smit et al. 2012; Moustakas et al. 2013; Steinhardt et al. 2014; Cassará et al. 2016; Citro et al. 2016) and late-type galaxies (that have on the average slowly declining star formation history over a long timescale of several gigayears; e.g., see Chiappini et al. 1997; Courteau et al. 2014; Pezzulli & Fraternali 2016; Grisoni et al. 2017). In this vein, off-main-sequence objects can be simply viewed as galaxies caught in an early evolutionary stage that are still accumulating their stellar mass (which grows almost linearly with time for a constant SFR), and are thus found to be preferentially located above the main sequence or, better, to the left of it. As time goes by and the stellar mass increases, the galaxy moves toward the average main-sequence relationship, around which it will spend most of its lifetime before being quenched due to gas exhaustion or feedback processes.

The third ingredient is the probability distribution of metallicity at given stellar mass and SFR:

$\begin{eqnarray}&&\displaystyle \frac{{dp}}{d\mathrm{log}Z}(Z| {M}_{\star },\psi )\propto \exp \left\{-\displaystyle \frac{{[\mathrm{log}Z-\mathrm{log}{Z}_{\mathrm{FMR}}({M}_{\star },\psi )]}^{2}}{2\,{\sigma }_{\mathrm{log}Z}^{2}}\right\},\end{eqnarray} \tag{ 6 }$

which is assumed to be log-normal with mean $\mathrm{log}{Z}_{\mathrm{FMR}}({M}_{\star },\psi )$ and scatter ${\sigma }_{\mathrm{log}Z}\approx 0.15$ dex provided by the observed fundamental metallicity relation (we adopt the determination by Mannucci et al. 2011; for an analytic fit, see their Equation (2)). This is a relationship among metallicity, stellar mass and SFR of a galaxy, which is found to be closely independent of redshift at least out to z ≲ 4 (see Mannucci et al. 2010, 2011; Hunt et al. 2016; Curti et al. 2020; Sanders et al. 2021). In the aforementioned in situ galaxy formation scenarios (see Pantoni et al. 2019; Lapi et al. 2020; Boco et al. 2021a), the fundamental metallicity relation naturally stems from an early, rapid increase of the metallicity with galactic age, which for a roughly constant star formation history, just amounts to the ratio M_⋆/ψ; a late-time saturation of the metallicity to a mass-dependent value which is determined by the balance among metal dilution, astration, removal by feedback, and partial restitution by stellar evolution processes and galactic fountains (if present).

In Figure 2, we illustrate the galactic term given by Equation (4). Specifically, we have color-coded the cosmic SFR density $d{\dot{M}}_{\mathrm{SFR}}/{dV}\,d\mathrm{log}Z$ per unit comoving volume as a function of redshift z and metallicity $\mathrm{log}Z$ . Most of the cosmic star formation occurs around redshift z ∼ 2–3 at rather high metallicity, close to solar values. The metallicity at which most star formation takes place decreases with redshift, but very mildly out to z ∼ 4–5. However, we note that at any redshift there is a long, pronounced tail of star formation occurring at low metallicities, down to 0.01 Z_⊙. We discuss possible variations of the prescriptions entering the galactic term in Section 4.1.

Figure 2. Refer to the following caption and surrounding text. — **Figure 2.** The galactic term $d{\dot{M}}_{\star }/{dV}\,d\mathrm{log}Z$ of Equation (4) expressing the amount of SFR (color-coded) per unit comoving volume as a function of redshift z (on x-axis) and of metallicity Z (on y-axis). The black horizontal line highlights the solar value Z_⊙.
Download figure:
Standard image High-resolution image

**Figure 2.** The galactic term $d{\dot{M}}_{\star }/{dV}\,d\mathrm{log}Z$ of Equation (4) expressing the amount of SFR (color-coded) per unit comoving volume as a function of redshift z (on x-axis) and of metallicity Z (on y-axis). The black horizontal line highlights the solar value Z_⊙.
Download figure:
Standard image High-resolution image

2.2. The Stellar Term

The stellar term in Equation (3) represents the number of BHs formed per unit of star-formed mass and remnant mass; as such, it is related to the evolution of isolated or binary stars, and must be evaluated via detailed simulations of stellar astrophysics. The stellar term can be split into various contributions:

$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{{dN}}_{\bullet }}{{{dM}}_{\mathrm{SFR}}d\mathrm{log}\,{m}_{\bullet }}({m}_{\bullet }| Z)=\displaystyle \frac{{{dN}}_{\star \to \bullet }}{{{dM}}_{\mathrm{SFR}}d\mathrm{log}\,{m}_{\bullet }}({m}_{\bullet }| Z)\\ \quad +\,\displaystyle \frac{{{dN}}_{\star \star \to \bullet }}{{{dM}}_{\mathrm{SFR}}d\mathrm{log}\,{m}_{\bullet }}({m}_{\bullet }| Z)+\displaystyle \sum _{i=1,2}\displaystyle \frac{{{dN}}_{\star \star \to \bullet \bullet }}{{{dM}}_{\mathrm{SFR}}d\mathrm{log}\,{m}_{\bullet ,i}}({m}_{\bullet }| Z).\end{array}\end{eqnarray} \tag{ 7 }$

The first comes from the evolution of isolated, massive stars that evolve into BHs at the end of their life (hereafter referred to as "single stellar evolution"); the second comes from stars that are originally in binary systems but end up as an isolated BH because one of the companion has been ejected or destroyed or cannibalized (hereafter "failed BH binaries"); the third comes from stars in binary systems that evolve into a binary BH with primary mass m_•,1 and secondary mass m_•,2 (hereafter "binaries"). All these terms are strongly dependent on metallicity Z, since this quantity affects the efficiency of the various processes involved in stellar and binary evolution, like mass loss rates, mass transfers, core-collapse physics, etc. Note that other crucial ingredients implicitly entering the above terms are the IMF ϕ(m_⋆), i.e., the distribution of star masses m_⋆ per unit mass formed in stars, and the binary fraction f_⋆⋆, i.e., the mass fraction of stars originally born in binary systems; both are rather uncertain quantities. As a reference, we adopt a universal (i.e., equal for every galaxies at any cosmic time) Kroupa (2001) IMF and a binary mass fraction f_⋆⋆ ≈ 0.5 constant with the star mass m_⋆. We discuss the impact of these choices on our results in Section 4.1.

To compute the stellar term, we exploit the outcomes of the SEVN stellar and binary evolution code, which directly provides each of the above contributions (see Spera et al. 2019 for details). However, to add some more grasp on the physics involved, we can provide a handy approximation for the first term referring to isolated stars:

$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{{dN}}_{\star \to \bullet }}{{{dM}}_{\mathrm{SFR}}d\mathrm{log}\,{m}_{\bullet }}({m}_{\bullet }| Z)=(1-{f}_{\star \star })\\ \times \ \int {{dm}}_{\star }\phi ({m}_{\star })\displaystyle \frac{{{dp}}_{\star \to \bullet }}{d\mathrm{log}\,{m}_{\bullet }}({m}_{\bullet }| {m}_{\star },Z),\end{array}\end{eqnarray} \tag{ 8 }$

where f_⋆⋆ is the binary fraction, ϕ(m_⋆) is the IMF, and

$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{{dp}}_{\star \to \bullet }}{d\mathrm{log}\,{m}_{\bullet }}({m}_{\bullet }| {m}_{\star },Z)\\ \propto \,\exp \left\{-\displaystyle \frac{{[\mathrm{log}\,{m}_{\bullet }-\mathrm{log}\,{m}_{\star \to \bullet }({m}_{\star },Z)]}^{2}}{2\,{\sigma }_{\mathrm{log}\,{m}_{\bullet }}^{2}}\right\}\end{array}\end{eqnarray} \tag{ 9 }$

is an approximately log-normal distribution centered around the metallicity-dependent relationship m_⋆→•(m_⋆, Z) between the remnant mass and the initial zero-age main-sequence star mass, with dispersion ${\sigma }_{\mathrm{log}\,{m}_{\bullet }}\approx 0.15$ dex (see Spera et al. 2015; Spera & Mapelli 2017; Boco et al. 2019). Unfortunately, for the other contributions of the stellar term related to binary evolution, an analogous analytic approximation is not viable, so one must trust the outputs of the SEVN code. We discuss the impact of adopting different prescriptions and codes for the computation of the stellar term in Section 4.1.

In Figure 3, we illustrate the stellar term of Equation (7), split in its various contributions. Specifically, we have color-coded the distributions ${{dN}}_{\bullet }/{{dM}}_{\mathrm{SFR}}\,d\mathrm{log}\,{m}_{\bullet }$ per unit star mass formed M_SFR and BH remnant mass m_•. The top left panel refers to isolated stars evolving in BHs. It is seen that a roughly constant number of remnants with masses m_• ∼ 5–30 M_⊙ is produced per logarithmic BH mass bin at any metallicity. Then there is a peak around m_• ∼ 30–60 M_⊙, with the larger values applying to more metal-poor conditions, where stellar winds are not powerful enough to substantially erode the stellar envelope before the final collapse. Finally, the distribution rapidly falls off for larger m_• ≳ 50–60 M_⊙ even at low metallicity,¹⁰ due to the presence of pair-instability and pulsational pair-instability supernovae (see Woosley et al. 2002; Belczynski et al. 2016; Spera & Mapelli 2017; Woosley 2017).

Figure 3. Refer to the following caption and surrounding text. — **Figure 3.** The stellar term ${{dN}}_{\bullet }/{{dM}}_{\mathrm{SFR}}\,d\mathrm{log}\,{m}_{\bullet }$ of Equation (7) expressing the number of BHs per unit star-formed mass (color-coded) as a function of BH mass m_• (on the x-axis) and of metallicity Z (on the y-axis). Different panels refer to isolated stars evolving into single BH (top left); binary stars failing to form a compact binary and instead originating a single BH (top right); binary stars evolving in binary BHs (bottom left); and a summation of these contributions (bottom right).
Download figure:
Standard image High-resolution image

**Figure 3.** The stellar term ${{dN}}_{\bullet }/{{dM}}_{\mathrm{SFR}}\,d\mathrm{log}\,{m}_{\bullet }$ of Equation (7) expressing the number of BHs per unit star-formed mass (color-coded) as a function of BH mass m_• (on the x-axis) and of metallicity Z (on the y-axis). Different panels refer to isolated stars evolving into single BH (top left); binary stars failing to form a compact binary and instead originating a single BH (top right); binary stars evolving in binary BHs (bottom left); and a summation of these contributions (bottom right).
Download figure:
Standard image High-resolution image

The top right panel refers to binary stellar systems failing to form a compact binary, and evolving instead into isolated BHs; this may happen because one of the progenitor stars has been ejected far away or destroyed or cannibalized during binary stellar evolution. The distribution of failed BH binaries differs substantially from single stellar evolution, being skewed toward more massive BHs, and with an appreciable number of remnants of mass m_• ∼ 50–160 M_⊙ produced especially at low metallicities. Such massive BH remnants are mainly formed when two (nondegenerate) companion stars merge during a common-envelope phase, possibly leaving a big BH remnant. Finally, note that, at high metallicity, a nonnegligible fraction of BHs in this channel has formed after the low-mass companion star had been ejected far away or destroyed by stellar winds and/or supernova explosions (see Spera et al. 2019 for details). The bottom left panel refers to binary stellar systems evolving into binary BHs; note that the distribution of primary and secondary BHs in the final configuration have been summed over. Although the overall number of binary BHs is substantially lower than the single BHs originated from the other two channels, most of them have a very similar mass spectrum; these remnants have formed from binary stars that underwent minor mass transfer episodes. However, at low metallicity, stellar winds are reduced and hence the mass exchanged or lost during binary evolution may be significantly larger, implying a more extended tail toward masses m_• ≲ 100 M_⊙ with respect to single stellar evolution. Finally, the bottom right panel illustrates the sum of all the previous formation channels.

2.3. Binary BH Mergers

Tight BH binaries may be able to progressively lose their energy via gravitational wave emission and to merge in a single, more massive BH. The cosmic merging rate of binary BHs can be computed as

$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{d{\dot{N}}_{\bullet \bullet \to \bullet }}{{dV}\,d\mathrm{log}\,{m}_{\bullet ,i}}({m}_{\bullet }| z)=\displaystyle \int {{dt}}_{{\rm{d}}}\displaystyle \int d\mathrm{log}Z\,\displaystyle \frac{{{dN}}_{\star \star \to \bullet \bullet \to \bullet }}{{{dM}}_{\mathrm{SFR}}d\mathrm{log}\,{m}_{\bullet ,i}}({m}_{\bullet }| Z)\\ \times \displaystyle \frac{{dp}}{{{dt}}_{{\rm{d}}}}({t}_{{\rm{d}}}| Z)\displaystyle \frac{d{\dot{M}}_{\mathrm{SFR}}}{{dVd}\mathrm{log}Z}(Z| {z}_{t-{t}_{{\rm{d}}}})\,\,\,\,\,\,\,\,\,i=1,2\end{array}\end{eqnarray} \tag{ 10 }$

where z_t is the redshift corresponding to a cosmic time t, and t_d is the time delay between the formation of the progenitor binary and the merger of the compact binary. The integrand in the expression above is the product of three terms. The rightmost one is the galactic term that must be computed at a time t − t_d. The middle one represents the probability distribution of delay times, possibly dependent on metallicity. The leftmost one is the stellar term that represents the number of binary stellar systems first evolving in a compact binary and then being able to merge within the Hubble time, per unit of formed star mass and compact remnant mass (primary m_•,1 and secondary m_•,2); this is a fraction of the third term on the right-hand side of Equation (7). The merger rate is then given by convolution of such a specific stellar term with the galactic term, weighted by the time delay distribution.

As a consequence of binary BH mergers, the stellar BH relic mass function may be somewhat distorted, since a number of low-mass BHs are shifted to larger masses. This amounts to applying a merging correction to the original cosmic birthrate $d\dot{N}/{dV}\,{{dm}}_{\bullet }$ , which can be computed as follows:¹¹

$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{d{\dot{N}}_{\bullet ,\mathrm{mergcorr}}}{{{dVdm}}_{\bullet }}({m}_{\bullet }| z)=\displaystyle \frac{d{\dot{N}}_{\bullet }}{{{dVdm}}_{\bullet }}({m}_{\bullet }| z)+\displaystyle \frac{1}{2}\,{\displaystyle \int }_{{m}_{\bullet }/2}^{{m}_{\bullet }}{{dm}}_{\bullet ,1}\,\\ \quad \times \,\displaystyle \frac{d{\dot{N}}_{\bullet \bullet \to \bullet }}{{dV}\,{{dm}}_{\bullet ,1}}({m}_{\bullet ,1}| z)\displaystyle \frac{\tfrac{d{\dot{N}}_{\bullet \bullet \to \bullet }}{{dV}\,{{dm}}_{\bullet ,2}}({m}_{\bullet }-{m}_{\bullet ,1}| z)}{{\displaystyle \int }_{0}^{{m}_{\bullet ,1}}{{dm}}_{\bullet ,2}\,\tfrac{d{\dot{N}}_{\bullet \bullet \to \bullet }}{{dV}\,{{dm}}_{\bullet ,2}}({m}_{\bullet ,2}| z)}\\ \\ \quad +\,\displaystyle \frac{1}{2}{\displaystyle \int }_{0}^{{m}_{\bullet }/2}{{dm}}_{\bullet ,2}\,\displaystyle \frac{d{\dot{N}}_{\bullet \bullet \to \bullet }}{{dV}\,{{dm}}_{\bullet ,2}}({m}_{\bullet ,2}| z)\\ \quad \times \,\displaystyle \frac{\tfrac{d{\dot{N}}_{\bullet \bullet \to \bullet }}{{dV}\,{{dm}}_{\bullet ,1}}({m}_{\bullet }-{m}_{\bullet ,2}| z)}{{\displaystyle \int }_{{m}_{\bullet ,2}}^{\infty }{{dm}}_{\bullet ,1}\,\tfrac{d{\dot{N}}_{\bullet \bullet \to \bullet }}{{dV}\,{{dm}}_{\bullet ,1}}({m}_{\bullet ,1}| z)}-\displaystyle \sum _{i=1,2}\,\displaystyle \frac{d{\dot{N}}_{\bullet \bullet \to \bullet }}{{dV}\,{{dm}}_{\bullet ,i}}({m}_{\bullet }| z);\end{array}\end{eqnarray} \tag{ 11 }$

in the positive terms on the right-hand side, the normalizations of the rates have been chosen so as to ensure mass conservation and self-consistency of the integrated merger rates.

We caveat that additional distorsions of the relic BH mass function may be induced by dynamical effects and/or hierarchical mergers in dense environments (see Section 4.2 for a preliminary discussion), by mass growth due to accretion in X-ray binaries or in AGN disks, etc.; these processes are expected to be relevant especially at the high-mass end, toward the intermediate-mass BH regime.

3. Results

In Figure 4, we illustrate the stellar BH cosmic birthrate of Equation (3) as a function of redshift z and remnant mass m_•. Specifically, we have color-coded the birthrate $d{\dot{N}}_{\bullet }/{dV}\,d\mathrm{log}\,{m}_{\bullet }$ per unit comoving volume V and BH mass m_•. At z ∼ 1–4, most of the BHs are formed with a rather flat distribution in $\mathrm{log}\,{m}_{\bullet }$ for m_• ∼ 5–50 M_⊙, though there is a nonnegligible tail up to m_• ≲ 160 M_⊙ due to stellar mergers (see previous section for details). Moving toward high redshift z ≳ 4 the mass distribution becomes skewed toward masses m_• ∼ 30–50 M_⊙ since these are preferentially produced in the low metallicity environments, while the tail for m_• ≳ 50 M_⊙ tends to be reduced because of the decrease in the number density of star-forming galaxies with appreciable SFR.

Figure 4. Refer to the following caption and surrounding text. — **Figure 4.** The stellar BH cosmic birthrate $d\dot{N}/{dV}\,d\mathrm{log}\,{m}_{\bullet }$ (color-coded) as a function of redshift z (on the x-axis) and BH mass m_• (on the y-axis).
Download figure:
Standard image High-resolution image

**Figure 4.** The stellar BH cosmic birthrate $d\dot{N}/{dV}\,d\mathrm{log}\,{m}_{\bullet }$ (color-coded) as a function of redshift z (on the x-axis) and BH mass m_• (on the y-axis).
Download figure:
Standard image High-resolution image

In Figure 5, we illustrate the stellar BH relic mass function ${{dN}}_{\bullet }/{dV}\,d\mathrm{log}\,{m}_{\bullet }$ of Equation (1) as a function of the remnant mass m_• for different redshifts z ∼ 0–10. At given redshift, the mass function features a roughly constant behavior for m_• ∼ 5–50 M_⊙, followed by a quite steep decline for m_• ≳ 50 M_⊙. Noticeably, there are bumps at around m_• ∼ 20 M_⊙, m_• ∼ 30–50 M_⊙, and m_• ∼ 120 M_⊙, which are more pronounced at high redshift (where metallicity is smaller) and progressively washed out toward the local Universe. Such features reflect regions of higher BH numbers that are evident in the stellar term of Figure 3 (in the form of almost vertical darker strips). Specifically, the first two originate from the detailed shape of the metallicity-dependent BH mass spectrum from single stellar evolution, while the third one mainly depends on binary evolution effects.

Figure 5. Refer to the following caption and surrounding text. — **Figure 5.** The stellar BH relic mass function ${dN}/{dV}\,d\mathrm{log}\,{m}_{\bullet }$ as a function of the BH mass m_• at different redshifts z ∼ 0 (cyan), z ∼ 1 (orange), z ∼ 2 (green), z ∼ 4 (red), z ∼ 6 (violet), z ∼ 8 (brown), and z ∼ 10 (pink).
Download figure:
Standard image High-resolution image

**Figure 5.** The stellar BH relic mass function ${dN}/{dV}\,d\mathrm{log}\,{m}_{\bullet }$ as a function of the BH mass m_• at different redshifts z ∼ 0 (cyan), z ∼ 1 (orange), z ∼ 2 (green), z ∼ 4 (red), z ∼ 6 (violet), z ∼ 8 (brown), and z ∼ 10 (pink).
Download figure:
Standard image High-resolution image

The mass function increases for decreasing redshift, quite rapidly down to z ∼ 2–3 and then more mildly toward z ∼ 0. Smoothing out minor features, the mass function can be analytically rendered in the range m_• ∼ 5–160 M_⊙ via a Schechter+Gaussian shape (Saunders et al. 1990):

$\begin{eqnarray}&&\begin{array}{c}\displaystyle \frac{{dN}}{{dV}\,d\mathrm{log}\,{m}_{\bullet }}\simeq { \mathcal N }{\left(\displaystyle \frac{{m}_{\bullet }}{{{ \mathcal M }}_{\bullet }}\right)}^{1-\alpha }\,{e}^{-{m}_{\bullet }/{{ \mathcal M }}_{\bullet }}+{{ \mathcal N }}_{{\rm{G}}}\\ \quad \times \,\displaystyle \frac{1}{\sqrt{2\pi {\sigma }_{G}^{2}}}\,{e}^{-{(\mathrm{log}\,{m}_{\bullet }-\mathrm{log}{{ \mathcal M }}_{\bullet ,{\rm{G}}})}^{2}/2\,{\sigma }_{{\rm{G}}}^{2}};\end{array}\end{eqnarray} \tag{ 12 }$

the Schechter function has normalization ${ \mathcal N }$ , low-mass power-law index α and characteristic mass ${{ \mathcal M }}_{\bullet }$ , while the Gaussian function (log-normal) has normalization ${{ \mathcal N }}_{{\rm{G}}}$ , mean ${{ \mathcal M }}_{\bullet ,G}$ and dispersion σ_G. The optimal values of these parameters have been determined via a Levenberg–Marquardt least-squares fit to the numerical mass function for m_• ∼ 5–160 M_⊙; the results for representative redshifts in the range z ∼ 0–10 are reported in Table 1.

Table 1. Fits to Stellar BH Relic Mass Function via Equation (12)

	Field, f_field = 1						Field+Cluster, f_field = 0.6
z	$\mathrm{log}{ \mathcal N }$	$\mathrm{log}{{ \mathcal M }}_{\bullet }$	α	$\mathrm{log}{{ \mathcal N }}_{{\rm{G}}}$	$\mathrm{log}{{ \mathcal M }}_{\bullet ,{\rm{G}}}$	σ_G	$\mathrm{log}{ \mathcal N }$	$\mathrm{log}{{ \mathcal M }}_{\bullet }$	α	$\mathrm{log}{{ \mathcal N }}_{{\rm{G}}}$	$\mathrm{log}{{ \mathcal M }}_{\bullet ,{\rm{G}}}$	σ_G
	(Mpc⁻³)	(M_⊙)		(Mpc⁻³)	(M_⊙)		(Mpc⁻³)	(M_⊙)		(Mpc⁻³)	(M_⊙)
0	5.623	0.607	−3.781	2.413	2.021	0.052	6.078	0.704	−2.717	3.496	1.808	0.1846
1	5.429	0.609	−3.859	2.309	2.023	0.051	5.887	0.709	−2.785	3.304	1.843	0.173
2	5.107	0.612	−3.914	2.064	2.024	0.051	5.592	0.713	−2.823	3.008	1.866	0.165
4	4.344	0.634	−3.902	1.419	2.037	0.049	4.796	0.747	−2.782	2.101	1.952	0.132
6	3.614	0.659	−3.866	0.806	2.054	0.045	4.112	0.785	−2.718	1.359	2.012	0.107
8	2.894	0.676	−3.868	0.197	2.066	0.043	3.457	0.816	−2.660	0.685	2.046	0.091
10	2.305	0.680	−3.884	-0.344	2.072	0.042	2.897	0.831	−2.623	0.113	2.059	0.0841

Note. Fits valid in the range m_• ∼ 5–160 M_⊙. Typical relative uncertainties on the parameters are ≲10%, and typical values of the reduced ${\chi }_{r}^{2}\lesssim 1$ are found. The results for f_field = 0.6 are based on the simulations by Di Carlo et al. (2019) and Ding et al. (2020) for young open star clusters.

Download table as: ASCII Typeset image

In Figure 6, we highlight the contribution of the different stellar evolution channels to the stellar BH relic mass function at z ∼ 0 and at z ∼ 10. In the range m_• ∼ 5–50 M_⊙, the single stellar evolution and the failed BH binaries channels are very similar and dominate over BH binaries. For m_• ≳ 50 M_⊙, the single stellar evolution contribution sharply dies (due to the mass gap from pair-instability SNe, see Section 2.2) and the binary BH channel abruptly decreases (due to mass loss in common-envelope phase, see Section 2.2), while the contribution from failed BH binaries dominates largely. Such a behavior in the relative contributions is basically independent of redshift.

Figure 6. Refer to the following caption and surrounding text. — **Figure 6.** The stellar BH relic mass function ${dN}/{dV}\,d\mathrm{log}\,{m}_{\bullet }$ as a function of the BH mass m_• at redshift z ∼ 0 (blue) and z ∼ 10 (green), highlighting the contribution to the total (solid lines) by single BHs formed from isolated stars (dotted–dashed lines), failed BH binaries (dashed lines), and BH binaries (dotted lines).
Download figure:
Standard image High-resolution image

In Figure 7, we illustrate the stellar BH relic mass density ρ_•(z) of Equation (2) as a function of redshift z. The mass density increases quite steeply toward smaller redshifts, changing from values ≲10⁵ M_⊙ Mpc⁻³ at z ∼ 10 to 10⁷ M_⊙ Mpc⁻³ around z ∼ 2–3 and then saturating for lower redshifts up to 5 × 10⁷ M_⊙ Mpc⁻³ at z ∼ 0. The evolution with redshift plainly follows that of the mass density ${\rho }_{\star }=(1-R){\int }_{0}^{z}{dz}^{\prime} \int d\mathrm{log}\psi \,\psi \,({dN}/{dV}/d\mathrm{log}\psi )$ in stars within galaxies, that has been computed by integrating the SFR functions and taking into account the recycling fraction R ≈ 0.45; this corresponds to a local energy density Ω_⋆ ≈ 5 × 10⁻³, consistent with classic estimates (e.g., Fukugita & Peebles 2004). On the gross average, we find a ratio of the mass densities ρ_•/ρ_⋆ ≈ 10%; we have checked that this mainly reflects the fraction of stars in the Kroupa IMF with mass m_⋆ ≳ 20 M_⊙ originating BH remnants, appropriately lowered by mass loss in stellar winds and by binary evolution effects. The local density in stellar-mass BHs is substantially larger than that in supermassive BHs with masses M_• ∼ 10⁶–10¹⁰ M_⊙, which amounts to ∼1–2 × 10⁵ M_⊙ Mpc⁻³ (see Shankar et al. 2020; see also Kormendy & Ho 2013; Aversa et al. 2015), implying that most of the BH mass in the local Universe is of stellar origin. The ratio of the local mass densities in supermassive to stellar-mass BHs is around ∼2 × 10⁻³. Finally, we estimate a local BH energy density Ω_• ≈ 4 × 10⁻⁴, corresponding to a BH-to-baryon ratio Ω_•/Ω_b ≲ 10⁻², i.e., the total mass in stellar BHs amounts to ≲1% of the local baryonic matter.

Figure 7. Refer to the following caption and surrounding text. — **Figure 7.** The relic density ρ_• in stellar-mass BHs as a function of redshift z (solid blue line); contribution from different BH mass ranges are also displayed: m_•/M_⊙ ≲ 20 (dotted), 20 ≲ m_•/M_⊙ ≲ 50 (dashed) and m_•/M_⊙ ≳ 50 (dotted–dashed). For reference, the stellar-mass density in galaxies is also illustrated (orange line). Moreover, the observational estimate of the mass density in supermassive BHs with M_• ∼ 10⁶–10¹⁰ M_⊙ at z ∼ 0 is shown (Shankar et al. 2020; red shaded area in the bottom left corner).
Download figure:
Standard image High-resolution image

For the sake of completeness, in Figure 7, we illustrate the contributions to the BH mass density from different mass ranges: m_• ≲ 20 M_⊙, m_• ∼ 20–50 M_⊙, and m_• ≳ 50 M_⊙. As expected from the shape of the mass function, most of the mass density is contributed by BHs with mass in the interval m_• ∼ 20–50 M_⊙. We also note that the mass density of BHs with masses m_• ≲ 20 M_⊙ tends to decline more steeply than the total one toward high redshift, because the production of such small-mass BHs is disfavored in the lower-metallicity environment of high-z galaxies; the opposite holds for the mass density of BHs with larger masses m_• ≳ 50 M_⊙.

In Figure 8, we evaluate the impact of binary BH mergers on the relic mass function. In the top left panel, we have color-coded the stellar term ${\sum }_{i=\mathrm{1,2}}\,{{dN}}_{\star \star \to \bullet \bullet \to \bullet }/{{dM}}_{\mathrm{SFR}}\,d\mathrm{log}\,{m}_{\bullet ,i}$ representing the number of binary stellar systems that first evolve in a BH binary and then are able to merge within the Hubble time, as a function of metallicity Z and BH mass m_•. The remnants contributing to this term are actually a subsample of those in the binary stellar term of Figure 3 (bottom left panel). It is evident that the number of merging binary BHs is a small fraction of the total, while its distribution lacks the high-mass tail. This is because, to form a tight compact binary, the progenitor stars must have undergone a substantial phase of common envelope, and an ensuing mass loss/envelope ejection (if not, the two progenitors merge prematurely and a single BH remnant is left). The end product is a binary that can be sufficiently tight to merge within the Hubble time, but typically made of BHs that cannot be more massive than m_• ≲ 40–50 M_⊙ (see Giacobbo & Mapelli 2018; Spera et al. 2019).

Figure 8. Refer to the following caption and surrounding text. — **Figure 8.** Impact of binary BH mergers on the stellar BH relic mass function. Top left panel: stellar term ${\sum }_{i=\mathrm{1,2}}\,{{dN}}_{\star \star \to \bullet \bullet \to \bullet }/{{dM}}_{\mathrm{SFR}}\,d\mathrm{log}\,{m}_{\bullet ,i}$ representing the number of binary stellar systems that first evolve in a BH binary and then are able to merge within the Hubble time (color-coded), as function of metallicity Z and BH mass m_•. Top right panel: cosmic merging rate ${\sum }_{i=\mathrm{1,2}}\,d{\dot{N}}_{\bullet \bullet \to \bullet }/{dV}\,d\mathrm{log}\,{m}_{\bullet ,i}$ (color-coded; see Equation (10)) as a function of BH mass m_• and redshift z. Bottom panel: stellar BH relic mass function ${{dN}}_{\bullet }/{dV}\,d\mathrm{log}\,{m}_{\bullet }$ as a function of BH mass m_• at redshift z ∼ 0, without (red lines) and with (green line) correction for binary BH merging effects; solid lines refer to the total mass function and dashed line to the contribution from binary BHs.
Download figure:
Standard image High-resolution image

**Figure 8.** Impact of binary BH mergers on the stellar BH relic mass function. Top left panel: stellar term ${\sum }_{i=\mathrm{1,2}}\,{{dN}}_{\star \star \to \bullet \bullet \to \bullet }/{{dM}}_{\mathrm{SFR}}\,d\mathrm{log}\,{m}_{\bullet ,i}$ representing the number of binary stellar systems that first evolve in a BH binary and then are able to merge within the Hubble time (color-coded), as function of metallicity Z and BH mass m_•. Top right panel: cosmic merging rate ${\sum }_{i=\mathrm{1,2}}\,d{\dot{N}}_{\bullet \bullet \to \bullet }/{dV}\,d\mathrm{log}\,{m}_{\bullet ,i}$ (color-coded; see Equation (10)) as a function of BH mass m_• and redshift z. Bottom panel: stellar BH relic mass function ${{dN}}_{\bullet }/{dV}\,d\mathrm{log}\,{m}_{\bullet }$ as a function of BH mass m_• at redshift z ∼ 0, without (red lines) and with (green line) correction for binary BH merging effects; solid lines refer to the total mass function and dashed line to the contribution from binary BHs.
Download figure:
Standard image High-resolution image

In the top right panel, we have color-coded the cosmic merging rate ${\sum }_{i=\mathrm{1,2}}\,d{\dot{N}}_{\bullet \bullet \to \bullet }/{dV}\,d\mathrm{log}\,{m}_{\bullet ,i}$ constituting the destruction term in Equation (10), as a function of BH mass m_• and redshift z. It is seen that most binary BH mergers occur for redshift z ∼ 1–5 in the mass ranges m_• ∼ 15–40 M_⊙ and m_• ∼ 5–8 M_⊙. In the bottom panel, we show the stellar BH relic mass function at z ∼ 0 with and without the correction for binary BH mergers; the solid lines are the total mass function, while the dashed lines are the contribution from BH binaries. As expected, the main changes are at the high-mass end, where an appreciable number of BHs in binaries are shifted toward larger masses by the merging process; however, the net effect on the total mass function is minor.

The contribution to the stellar BH mass function from merging binary BHs can be probed via gravitational wave observations. Recently, the LIGO/Virgo collaborations (Abbott et al. 2021a; see also Baxter et al. 2021) have estimated the primary mass distribution for BH binaries that coalesce around z ≈ 0; given the quite small cosmological volume probed by the current gravitational wave detectors, this approximately corresponds to the merger rate integrated in a narrow redshift interval Δz ≲ 0.5. The expectation from this work is illustrated as a red solid line in Figure 9, and compared with the estimates by Abbott et al. (2021a) (blue shaded area) for a power-law + peak model and by Baxter et al. (2021) (orange shaded area) for their astrophysically motivated model (for both determinations the 68% and 95% credible intervals are represented with dark and light shades). Our result is in remarkable agreement with these estimates up to 40 M_⊙. However, the observed primary mass distribution declines gently for m_• ≳ 40 M_⊙ out to m_• ∼ 80–100 M_⊙, while our model dies off, since stellar evolution effects hinder the presence of very massive BHs in coalescing binaries. This occurs mainly for two reasons (see also Section 3): the mass gap m_• ∼ 50–120 M_⊙ for the production of BH due to pair-instability and pulsational pair-instability supernovae; and the substantial mass loss during the common-envelope phase needed to produce a hardened compact binary that can merge within reasonable timescales. We also stress that such a sharp decline is not dependent on the specific galactic prescriptions nor scatter in the adopted relations, but it is instead common to any approach including the production and possible merging of only isolated BH binaries. A viable solution is explored in the next section.

Figure 9. Refer to the following caption and surrounding text. — **Figure 9.** The BH mass function of merging BH binaries as a function of primary BH mass at z ≈ 0. The outcome of this work is illustrated by the red solid line. Estimates from the analysis of gravitational wave observations by Abbott et al. (2021a) and Baxter et al. (2021) are reported as blue and orange shaded areas, respectively (for both determinations the 68% and 95% credible intervals are represented with dark and light shades).
Download figure:
Standard image High-resolution image

4. Discussion

In this section, we discuss three interesting issues that could potentially affect our results: the impact of the adopted modeling prescriptions, the dynamical formation channel of BH binaries in dense environment such as star clusters, and the implication of our work in the formation of light (super)massive BH seeds at high redshift.

4.1. Impact of Modeling Prescriptions

We warn the reader that the stellar BH relic mass function derived in the present work is somewhat dependent on the prescriptions entering the galactic and the stellar terms discussed in Sections 2.1 and 2.2. It is beyond the scope of this paper to address the issue in detail by exploring the variety of modeling recipes, numerical approaches, and associated parameter space present in the literature. However, to provide the reader with a grasp of the related impact on our main results, we discuss the main source of uncertainties and show their impact on the local BH mass function.

As to the galactic term, the main source of uncertainty is constituted by the adopted main sequence and fundamental metallicity relationships. To estimate the related degree of uncertainty, we recomputed the galactic term by using in turn the main-sequence relation by Boogaard et al. (2018) in place of our reference by Speagle et al. (2014), and the fundamental metallicity relation by Hunt et al. (2016) in place of our reference by Mannucci et al. (2011). The results are illustrated in Figure 10, and it is seen that the related changes on the stellar BH relic mass function are minor.

Figure 10. Refer to the following caption and surrounding text. — **Figure 10.** Impact on the BH mass function of different prescriptions for the galactic term (see Section 4.1). Top panels: the galactic term of Equation (4) when assuming our reference Speagle et al. (2014) main sequence and Mannucci et al. (2011) fundamental metallicity relation (left), when changing the fundamental metallicity relation to that by Hunt et al. (2016; middle), and when changing the main-sequence relation to that by Boogaard et al. (2018; right). Bottom panel: the relic stellar BH mass function at redshift z ∼ 0 for the three cases just described above.
Download figure:
Standard image High-resolution image

As to the stellar term, the main source of uncertainty is constituted by the numerical treatment in the SEVN code of stellar winds, of pair-instability and pulsational pair-instability supernovae, of supernova explosions and natal kicks, of mass transfers and common-envelope phase in binary systems: all these physical processes can have profound impact on the formation and evolution of isolated or binary stellar BHs (see reviews by Postnov & Yungelson 2006; Ivanova et al. 2013; Mapelli 2020; and references therein). To estimate the related degree of uncertainty, we recomputed the stellar term for binaries/failed binaries and the associated relic BH mass function by running a different binary stellar evolution code. We rely on COSMIC (Breivik et al. 2020) since it constitutes a state-of-the-art, community-developed software with extensive online documentation on the Github repository (see https://cosmic-popsynth.github.io/) and it upgrades the classic and widespread code BSE (Hurley et al. 2002). As for the parameters regulating mass transfer, common envelope, tidal evolution, SN explosion, and compact-object birth kicks, we keep the same choices adopted for SEVN, which are described in Spera et al. (2019). The results are presented in Figure 11. The top left and middle panels compare the (total) stellar term from the two codes. It is apparent that the main difference involves the distribution of remnants with mass m_• ≳50–60 M_⊙ that are associated with compact binaries and especially with failed BH binaries (i.e., binary stars have coalesced before evolving into two distinct binary BHs; see also Section 2.2). Specifically, COSMIC tends to produce fewer failed BH binaries than SEVN; this is related to the treatment of the common-envelope stage, which is indeed one of the most uncertain processes in the evolution of binary stars. The bottom panel illustrates the effects on the relic stellar BH mass function at z ∼ 0. All in all, the differences are minor out to m_• ∼ 60 M_⊙. At larger masses, the mass function built from the COSMIC stellar term decreases rapidly, to imply a paucity of remnants with masses m_• ≳ 100 M_⊙; contrariwise, the mass function built from the SEVN stellar term is decreasing mildly and actually features a bump around m_• ∼ 150 M_⊙, before a cutoff. We stress that the major differences occur in a range where the mass function ${dN}/{{dm}}_{\bullet }$ is already decreasing much faster than ${m}_{\bullet }^{-1}$ .

Figure 11. Refer to the following caption and surrounding text. — **Figure 11.** Impact on the BH mass function of different codes and prescriptions for the computation of the stellar term (see Section 4.1). Top panels: the (total) stellar term of Equation (7) from the `SEVN` (left) and the `COSMIC` (middle) binary stellar evolution codes assuming a common-envelope parameter α_CE = 5 (see Section 4.1); results with `COSMIC` (right) for α_CE = 1 are also shown. Bottom panel: the relic stellar BH mass function at redshift z ∼ 0 from the `SEVN` (solid blue line) and the `COSMIC` (dotted blue line) with α_CE = 5, and from `COSMIC` (dashed blue line) with α_CE = 1; the contribution from binaries is also highlighted (green lines).
Download figure:
Standard image High-resolution image

Among the many parameters entering binary evolution codes, a major role is played by the common-envelope parameter α_CE, that quantifies the energy available to unbind the envelope (see Ivanova et al. 2013). As in Spera et al. (2019), we have set α_CE = 5 throughout this paper and in the above comparison between SEVN and COSMIC. To have a grasp on the impact of such a choice, we have also run COSMIC with α_CE = 1, and show in Figure 11 the resulting stellar term (top right panel) and BH mass function at z = 0 (bottom panel, dashed lines). The impact of α_CE is rather limited, and mainly affects the high-mass tail of the distribution, which is dominated by binary stellar evolution.

Another crucial parameter entering the stellar term is the binary mass fraction f_⋆⋆, for which we have assumed a constant value 0.5. However, this is a rather uncertain quantity, that may well depend on the star mass m_⋆, on the properties of the host galaxy and/or of the local environment, and even on redshift. In fact, observational constraints (e.g., Raghavan et al. 2010; Sana et al. 2012; Li et al. 2013; Sota et al. 2014; Dunstall et al. 2015; Luo et al. 2021) suggest values in the range f_⋆⋆ ≈ 0.3–0.7, with a possible increase toward more massive stars. To bracket the effects of different binary fractions on the stellar term and on the relic BH mass function, in Figure 12, we compare the results for our reference f_⋆⋆ = 0.5 to the extreme cases f_⋆⋆ = 0 and f_⋆⋆ = 1. Plainly adopting f_⋆⋆ = 0 (i.e., no stars born in binaries) cuts the tail of the stellar term and of the BH mass function for m_• ≳ 50–60 M_⊙, that are mostly produced by binary stellar evolution effects; meanwhile, for smaller masses, the mass function is increased by the relatively more numerous single stars. Contrariwise, the other extreme case f_⋆⋆ = 1 (all stars born in binaries), tends to enhance the high-mass tail of the mass function and reduce somewhat the number density of BHs with masses m_• ≲ 50–60 M_⊙. All in all, it is seen that the differences on the mass function are minor, and especially so for f_⋆⋆ ≳ 0.5.

Figure 12. Refer to the following caption and surrounding text. — **Figure 12.** Impact of the binary fraction f_⋆⋆ on our results (see Section 4.1). Top panels: the (total) stellar term of Equation (7) for f_⋆⋆ = 0 (left; no stars born in binaries), 0.5 (middle; our fiducial case), and 1 (right; all stars born in binaries). Bottom panel: the relic stellar BH mass function at redshift z ∼ 0 for f_⋆⋆ = 0 (dotted line), 0.5 (solid), and 1 (dashed).
Download figure:
Standard image High-resolution image

A final caveat about the IMF is in order. Throughout the paper, we have self-consistently adopted as a reference the Kroupa (2001) IMF, mainly because it constitutes a standard both in the galaxy formation and the stellar evolution communities (together with the Chabrier (2003) IMF, which anyway would yield almost indistinguishable results) and because we had prompt availability of stellar evolution simulations based on that. Other classic choices like the Salpeter (1955), the Kennicut (1983), and the Scalo (1986) IMF differ somewhat for the relative amount of star formation occurring below and above 1 M_⊙; such a difference is exacerbated in bottom-heavy IMFs like the one proposed to apply in massive ellipticals (e.g., van Dokkum & Conroy 2010) or top-heavy IMFs like the one proposed to apply in the early Universe (e.g., Larson 1998) or in starburst galaxies (e.g., Lacey et al. 2010). It is important to stress that the IMF enters nontrivially both in the galactic and in the stellar term. In the galactic term, an IMF is needed to convert the observed galaxy luminosity functions into a statistics based on an intrinsic physical quantity such as the SFR or the stellar mass; moreover, the determination itself of SFR, stellar masses (hence of the main sequence of star-forming galaxies), star formation history, and metallicity via broadband SED fitting rely on the assumption of a specific IMF. Typically, at given SFR, more top-heavy IMFs are proportionally richer in massive short-lived stars, and tend to yield a larger restframe UV luminosity/ionizing power, a faster chemical enrichment, and a smaller stellar mass locked in long-lived stars. In the stellar term, the IMF is plainly an essential ingredient in determining the relative proportion of compact remnants of different masses (see Equation (8)). More top-heavy IMFs tend to produce more massive stars per unit SFR, thus in principle, more massive BH remnants (and a smaller relative amount of neutron stars), although this is not so straightforward due to mass loss and mass transfers during binary evolution. From the above considerations, one could be led to speculate that precision determinations of the BH mass function in an extended mass range would constitute a probe for the IMF; however, this plan, at least presently, struggles against the many and large uncertainties associated with the treatment of binary evolution effects, and against the unexplored degeneracies with the bunch of poorly constrained parameters of galaxy and stellar evolution.

4.2. Dynamical Channel

In Figure 9, we highlighted a possible mismatch at m_• ≳ 40 M_⊙ between the primary mass distribution for merging BH binaries from isolated binary evolution and the estimates from gravitational wave observations. A viable solution could be that such large primary masses are produced in binary systems formed within the dense environment of young stellar clusters, open clusters, globular clusters, or nuclear star clusters (e.g., Rodriguez et al. 2015, 2021; Antonini & Rasio 2016; Di Carlo et al. 2019, 2020; Kumamoto et al. 2019; Arca Sedda et al. 2020; Banerjee 2021; Mapelli et al. 2021). The central density of a star cluster can be so high that the orbits of binary stars are continuously perturbed by dynamical encounters with other members. Massive BHs m_• ≳ 40 M_⊙ in the pair-instability mass gap can then be originated by hardening of BH binaries via dynamical exchanges in three-body encounters, and via the merging of massive progenitor stars; in addition, runaway collisions (i.e., a fast sequence of mergers; e.g., Portegies Zwart et al. 2004; Giersz et al. 2015; Mapelli 2016) in the densest cores of clusters with low metallicity can even produce intermediate-mass BHs with m_• ≳ some 10² M_⊙.

To have a grasp on these effects, we proceed as follows. First, we construct the stellar term of Equation (7) from the simulations by Di Carlo et al. (2020), which include dynamical effects in young star clusters. With respect to isolated conditions, an appreciable number of merging binaries with a primary mass m_• ≳ 40 M_⊙ is originated via dynamical exchanges, especially at low metallicities. In more detail, BHs with mass m_• ∼ 40–65 M_⊙ can be formed (even in the field) from the evolution of single stars or stars in loose binaries that retain a fraction of their envelope up to their final collapse. BHs with mass m_• ≳ 65 M_⊙ can be originated via collisions of massive stars (eased in star clusters because dynamical encounters can induce fast merging before mass transfer episodes peel-off the primary star). If produced in the field, BHs from both these channels will remain isolated or locked in loose binaries unable to merge; contrariwise, in a star cluster they can acquire a new companion via dynamical exchanges and merge by dynamical hardening and gravitational wave emission. Note that, in contrast, repeated mergers of BHs are suppressed in young star clusters, because of their low escape velocity.

Second, we assume that a fraction f_field of the star formation occurs in the field and the complementary fraction 1 − f_field occurs in young star clusters (actually most of the stars are formed in young star clusters, but only a fraction of these may be subject to dynamical effects before exiting from the cluster or before the star cluster itself dissolves). Observations (see Goddard et al. 2010; Johnson et al. 2016; Chandar et al. 2017; Adamo et al. 2020) and cluster formation models (Kruijssen 2012; Elmegreen 2018; Pfeffer et al. 2018; El-Badry et al. 2019; Grudic et al. 2021) indicate that such a fraction is highly uncertain and possibly dependent on properties like the SFR spatial density and redshift; in this exploratory computation, we let the fraction f_field vary from 0.2 to 1 (which corresponds to isolated binaries only), and we split the galactic term of Equation (4) accordingly. Finally, we combine the stellar and galactic term so derived to compute the merger rate of Equation (10) and the expected primary mass distribution.

The outcome is illustrated in the top panel of Figure 13 for different values of f_field in the range 0.2–1. As expected, increasing the fraction of SFR in star clusters (i.e., decreasing f_field) produces a progressively more extended tail toward high primary masses, to the point that values f_field ≲ 0.8 actually can reconcile the theoretical prediction with the observational estimates. For completeness, in the bottom panel of Figure 13, we show how the dynamical evolution channel affects the relic BH mass function at z ∼ 0. The marked difference with respect to the model with only isolated binaries is the absence of the drop at around m_• ∼ 60 M_⊙ and of the abrupt cutoff for m_• ∼ 150 M_⊙. Instead, the mass function declines smoothly for m_• ≳ 60 M_⊙. The analytical fits in terms of Equation (12) for the field + cluster mass function with f_field = 0.6 at representative redshifts are reported in Table 1.

Figure 13. Refer to the following caption and surrounding text. — **Figure 13.** Top panel: impact of the dynamical formation channel on the BH primary mass distribution of merging BH binaries z ≈ 0, computed from the simulations for young open star clusters by Di Carlo et al. (2020, see Section 4.2 for details). Lines refer to different fractions of star formation occurring in the field f_field ≈ 1 (solid; only isolated binaries), 0.8 (dashed), 0.6 (dotted–dashed), 0.4 (dotted–dotted–dashed), and 0.2 (dotted). Estimates from the analysis of gravitational wave observations by Abbott et al. (2021a) and Baxter et al. (2021) are reported as blue and orange shaded areas, respectively (for both determination the 68% and 95% credible intervals are represented with dark and light shades). Bottom panel: impact of the dynamical formation channel on the relic BH mass function at z ∼ 0. Dotted lines refer to the contribution from star clusters and solid lines to the total (field + star clusters) BH mass function. Color code refers to different fractions of star formation occurring in the field f_field ≈ 1 (red), 0.8 (blue), 0.6 (green), 0.4 (purple), and 0.2 (yellow).
Download figure:
Standard image High-resolution image

Two caveats are in order here. First, in the simulations by Di Carlo et al. (2019, 2020) the treatment of stellar mergers is based on simplified assumptions: no mass loss and chemical mixing during the merger, instantaneous recovery of hydrostatic equilibrium after the merger, and rejuvenation of the merger product according to the simple Hurley et al. (2002) prescriptions. Plainly, these details of the process are quite uncertain, and dedicated hydrodynamical simulations are required to have a better understanding of the final outcome. Second, in estimating the impact of the dynamical formation channel, we have considered only young star clusters for the sake of simplicity (and because the prompt availability of in-house dynamical simulations, which are extremely time demanding to run from scratch). Globular clusters and nuclear star clusters could also be effective environments to build up massive binary BHs, since hierarchical mergers are more efficient in very rich and compact stellar systems (e.g., Miller & Hamilton 2002; Antonini et al. 2019; Mapelli et al. 2021). Hence, including models for globular and nuclear star clusters could allow to reproduce the observed primary BH mass function with an even larger value of f_field.

4.3. BH Seeds at High Redshift

We stress that the stellar BH mass function at high redshift z ≳ 6 derived in the present paper actually provides a light BH seed distribution in primordial galaxies, as originated by stellar and binary evolution processes. This is complementary to the seed distributions expected from other classic formation channels. The most relevant are basically three (see review by Volonteri 2010 and references therein): BHs formed by Population III stars in metal-free environments like high-redshift minihalos at z ≳ 20 (e.g., Madau & Rees 2001), runaway stellar collisions in metal-poor nuclear star clusters (e.g., Devecchi et al. 2012), and direct collapse scenarios (e.g., Lodato & Natarajan 2006).

The light seed distributions from these models (extracted from Figure 4 of Volonteri 2010) are compared in Figure 14 with the mass function at z ∼ 6–10 from stellar BHs presented in this work. Given the appreciable contribution of the latter for masses m_• ≲ 150 M_⊙ at redshift z ≲ 10, it will be extremely relevant to include it in numerical simulations, semi-analytic and semiempirical models of BH formation and evolution at high redshift.

Figure 14. Refer to the following caption and surrounding text. — **Figure 14.** The light seed distribution at redshift z ∼ 10 (solid), 8 (dotted–dashed), 6 (dashed), and 0 (dotted) from the stellar BH mass function of this work is compared with that expected from Population III remnants (orange line), runaway stellar collisions in nuclear star clusters (green line), and direct collapse scenarios (blue line); see discussion in Section 4.3.
Download figure:
Standard image High-resolution image

5. Summary

In this work, we have provided an ab initio computation of the relic stellar black hole (BH) mass function across cosmic times. To this purpose, we have exploited the state-of-the-art stellar and binary evolutionary code SEVN, and have coupled its outputs with redshift-dependent galaxy statistics and empirical scaling relations involving galaxy metallicity, star formation rate, and stellar mass. Our main findings are summarized below.

1.
The relic mass function ${dN}/{dVd}\mathrm{log}\,{m}_{\bullet }$ as a function of the BH mass m_• features a rather flat shape up to m_• ≈ 50 M_⊙ and then a log-normal decline for larger masses, while its normalization increases with decreasing redshift, quite rapidly down to z ∼ 2–3 and then more mildly toward z ∼ 0, see Figure 5 and Table 1. The local stellar BH mass function could be eventually probed via microlensing observations (see Paczynski 1986; for a review, Mao 2012).
2.
The local relic BH mass function for m_• ≲ 50 M_⊙ is comparably contributed from isolated stars evolving into BH and from binary stellar systems ending up in single BH (failed BH binaries), while binary BHs are subdominant; for higher masses m_• ≳ 50 M_⊙ the single stellar evolution and binary BH channels abruptly decrease, while the failed BH binaries dominates largely. See Figure 6.
3.
The local stellar BH relic mass density amounts to ρ_• ≈ 5 × 10⁷ M_⊙ Mpc⁻³, exceeding that in supermassive BHs by more than two orders of magnitude; this translates into an energy density parameter Ω_• ≈ 5 × 10⁻⁴, implying that the total mass in stellar BHs amounts to ≲1% of the local baryonic matter. See Figure 7.
4.
The stellar BH relic mass function can be distorted by binary BH mergers, which can redistribute remnants from the low to the high-mass range. However, such a reshaping is found to have a minor effect on the mass function at the high-mass end. See Figure 8.
5.
The distribution of primary masses for merging BH binaries is found to be in remarkable agreement with the recent estimates from gravitational wave observations by LIGO/Virgo out to m_• ≲ 40 M_⊙. For larger masses, the observed distribution declines gently while the theoretical one dies off due to a twofold reason: the mass gap from pair-instability and pulsational pair-instability supernovae; and the substantial mass loss during the common-envelope phase needed to produce a tight BH binary that can merge within the Hubble time. We have proposed as a viable solution to consider the dynamical formation of merging BH binaries in dense environment like (young) star clusters. See Figures 9–13 and Table 1.
6.
We have discussed the impact on the mass function of adopting different physical prescriptions entering the galactic and the stellar term. As for the galactic term, minor differences arise when changing the adopted main sequence of star-forming galaxies and the fundamental metallicity relationship. As for the stellar term, the main differences concern the high-mass end of the mass function, and are due to the numerical treatment of binary stellar evolution effects, which are still considerably uncertain even from a theoretical point of view. See Figures 11 and 10.
7.
The BH mass function derived here can provide a firm theoretical basis for a physically motivated light seed distribution at high redshift; as expected, for masses m_• ≲ 150 M_⊙ and redshifts z ≲ 10 it overcomes complementary light seed formation channels, like Population III stars, stellar mergers in nuclear star clusters and direct collapse scenarios that are aimed to produce larger mass seeds. It will be worth implementing the light seed distribution from the present work in semi-analytic and numerical models of (super)massive BH formation and evolution, as it can be quite relevant for predictions concerning future gravitational wave observations via Einstein Telescope and LISA. See Figure 14.

In a future perspective, it would be interesting to exploit the galactic and stellar evolution prescriptions adopted here to populate an N-body simulation via subhalo clustering abundance matching technique (e.g., Ronconi et al. 2020) in order to derive a mock catalog encapsulating the three-dimensional spatial distribution of stellar-mass BHs and of their galactic hosts. Finally, we stress that our work can constitute a starting point to investigate the origin of heavy seeds and the growth of (super)massive BHs in high-redshift star-forming galaxies, that we will pursue in forthcoming papers.

We thank the referee for a competent and constructive report. We acknowledge C. Baccigalupi and M. Massardi for interesting discussions and critical reading. This work has been supported by the EU H2020-MSCA-ITN-2019 Project 860744 "BiD4BESt: Big Data applications for black hole Evolution STudies". A.L. acknowledges funding from the PRIN MIUR 2017 prot. 20173ML3WW, "Opening the ALMA Window on the Cosmic Evolution of Gas, Stars and Supermassive Black Holes". U.N.D.C. and M.M. acknowledge financial support from the European Research Council for the ERC Consolidator grant DEMOBLACK, under contract no. 770017.