The Mass Function of Supermassive Black Holes in the Direct-collapse Scenario

A key challenge to the theory of the formation of supermassive black holes (SMBHs) in the early universe is the observation of very massive (M ≈ 10⁹ ${M}_{\odot }$) and luminous (L ≳ 10¹³ ${L}_{\odot }$) quasars already in place by z ∼ 7, when the universe was just ∼800 Myr old (e.g., Fan et al. 2006; Mortlock et al. 2011; Wu et al. 2015; Bañados et al. 2018; see also the review by Woods et al. 2018). Objects that accumulate at least a billion ${M}_{\odot }$ in less than a billion years after the big bang put a strain on the normal ideas of Eddington-limited growth of black hole seeds that originate from Population III stellar remnants. Starting from a seed mass M₀, Eddington-limited growth leads to a mass $M(t)={M}_{0}\exp [{\epsilon }^{-1}(1-\epsilon )t/{t}_{{\rm{E}}}]$, where ${t}_{{\rm{E}}}\approx 450\,\mathrm{Myr}$ and (≈0.1) is a radiative efficiency factor. Population III stars are thought to have masses ≲40 ${M}_{\odot }$ (Hosokawa et al. 2011), and their remnants would be less massive, so that there is apparently not enough time available to reach M ∼ 10⁹ ${M}_{\odot }$. These constraints show that a combination of both more-massive initial seeds and a super-Eddington growth rate may be necessary to account for the observed SMBHs at z ∼ 7.

One promising pathway is that of direct-collapse black holes (DCBHs; Bromm & Loeb 2003). The idea is that Lyman–Werner (LW) photons (having energies 11.2–13.6 eV) from the first Population III stars can propagate far from their sources and dissociate H₂ in other primordial gas clouds. Without H₂ cooling these gas clouds equilibrate to temperatures T ∼ 8000 K set by atomic cooling, which means that the Jeans mass is ∼10⁵ ${M}_{\odot }$ at a number density n = 10⁴ cm⁻³, as opposed to ∼10³ ${M}_{\odot }$ in a normal Population III star formation environment or ∼1 ${M}_{\odot }$ in present-day star formation. Due to their large masses, these collapsing cloud fragments may be able to collapse directly into black holes, after a brief period as a supermassive star (Bromm & Loeb 2003) or quasi-star (Begelman et al. 2006, 2008), if the infalling matter can overcome the angular momentum barrier and disruptive effects of radiative feedback. An interesting joint solution to these barriers is proposed by Sakurai et al. (2016) based on the episodic accretion scenario of Vorobyov et al. (2013) that is powered by gravitational instability in a circumstellar disk. In this model the episodic accretion results in a lower surface temperature of a supermassive star, thereby also reducing the effect of radiative feedback that can limit mass accumulation in the case of normal Population III star formation (Hosokawa et al. 2011). The DCBH model has been extensively developed in the context of galaxy-formation models, resulting in a scenario where the formation of atomic cooling halos is seeded by the first stars, and the subsequent DCBH produce LW radiation that triggers the formation of other atomic cooling halos and DCBH in a kind of chain reaction process (Yue et al. 2014). A rapid period of growth of atomic cooling halos, and therefore DCBH formation, ensues, with the growth rate at any time related to the instantaneous number of DCBHs. The rapidly growing phase of DCBH creation is also a period of possible rapid mass growth through super-Eddington accretion (Inayoshi & Haiman 2016; Pacucci et al. 2017). The whole process comes to a rapid halt, however, when the gas in the atomic cooling halos is photo-evaporated by the ambient radiation field. According to Yue et al. (2014) the DCBH era lasts from z ≈ 20 to z ≈ 13, or a time period T ≈ 150 Myr, after which DCBH formation is completely suppressed. Here, we adopt the picture emerging from their semi-analytic model; however, we note that numerical simulations (e.g., Agarwal et al. 2012; Chon et al. 2016; Habouzit et al. 2016) have not reached a consensus on the DCBH formation rate or the termination time of their formation.

In this Letter, we seek a simple model of the growth of DCBHs in the early universe that captures just the essential features of the scenario described above, in order to reach an analytic understanding of the mass and luminosity functions of observable quasars that form through the DCBH scenario.

Yue et al. (2014) estimated that the rapid formation of DCBHs occurs between z ≈ 20 and z ≈ 13, after which it ends abruptly. The semi-analytic model of Dijkstra et al. (2014) is broadly consistent with this and shows that the number density of DCBHs ${n}_{\mathrm{DCBH}}$ grows rapidly during a similar interval.

Figure 1 shows the results of Dijkstra et al. (2014) plotted against time after the big bang, for two possible values of the critical flux of LW photons ${J}_{\mathrm{crit}}$ (written in units of 10⁻²¹ erg s⁻¹ cm⁻² Hz⁻¹ sr⁻¹) that is needed to create the atomic cooling halos, which are DCBH progenitors. The calculation of ${J}_{\mathrm{crit}}$ (e.g., Shang et al. 2010) includes the additional effect of near-infrared radiation with energy above 0.76 eV that can also inhibit the formation of H₂ by dissociating the H⁻ ion, which is an intermediary in the H₂ formation chain. Although significant uncertainties exist in the actual number of DCBHs created, due to uncertainties in the value of ${J}_{\mathrm{crit}}$ that depends on the source density and spectra, escape fraction from the host halos, etc., the slope of growth, $\lambda (t)\equiv d{\rm{ln}}{n}_{{\rm{DCBH}}}/{dt}$, is similar in all of their modeled cases. For the two cases shown here, we find an average value λ measured between data points from z = 20 to z = 12, which corresponds roughly to the period of rapid DCBH formation. For their canonical model ${J}_{\mathrm{crit}}=300$ the best-fit line yields λ = 27.7 Gyr⁻¹, and for ${J}_{\mathrm{crit}}=100$ the best fit is λ = 17.7 Gyr⁻¹. Here we use the canonical model and adopt λ = 28.0 Gyr⁻¹, which is slightly steeper than the best fit. This in the interest of rounding off, and also because there is evidence that the DCBH growth era ended just prior to this data point, at z ≈ 13 (Yue et al. 2014); this would tend to drive the slope to a slightly greater value.

Zoom In Zoom Out Reset image size

Each DCBH can be modeled as growing in mass at an exponential rate, but the starting times of the accretion process will be spread throughout the DCBH formation era. However, the super-Eddington growth will cease for all DCBHs at about the same time, so that there will be a distribution of accretion times among the population of DCBH. This is a key part of our model as developed in Section 3.

The growth of an individual DCBH is thought to proceed by default at an Eddington-limited rate, but periods of super-Eddington growth are also possible (Pacucci et al. 2017). The Eddington luminosity is

$$\begin{eqnarray}&&{L}_{{\rm{E}}}=\displaystyle \frac{4\pi {{cGm}}_{{\rm{p}}}M}{{\sigma }_{T}}.\end{eqnarray} \tag{ 1 }$$

where M is the black hole mass, ${m}_{{\rm{p}}}$ is the proton mass, and ${\sigma }_{T}={(8\pi /3)({e}^{2}/{m}_{{\rm{e}}}{c}^{2})}^{2}$ is the Thomson cross section in which ${m}_{{\rm{e}}}$ is the electron mass. At this luminosity the radiation pressure can balance the gravitational pressure. The accretion of mass to very small radii comparable to the Schwarzschild radius will release a significant portion of the rest mass energy, hence the luminosity is normally estimated as ${L}_{\mathrm{acc}}=\epsilon {\dot{M}}_{\mathrm{acc}}{c}^{2}$, where ${\dot{M}}_{\mathrm{acc}}$ is the mass accretion rate and is the radiative efficiency, typically set to 0.1. Because the accretor will gain rest mass at the rate $\dot{M}=(1-\epsilon ){\dot{M}}_{\mathrm{acc}}$, we equate ${L}_{\mathrm{acc}}$ with ${L}_{{\rm{E}}}$ to find that

$$\begin{eqnarray}&&\displaystyle \frac{{dM}}{{dt}}={\gamma }_{0}\,M\Rightarrow M(t)={M}_{0}\exp ({\gamma }_{0}t),\end{eqnarray} \tag{ 2 }$$

where ${\gamma }_{0}=(1-\epsilon )/(\epsilon {t}_{{\rm{E}}})$ and ${t}_{{\rm{E}}}=2{e}^{4}/(3{{Gm}}_{{\rm{p}}}{m}_{{\rm{e}}}^{2}{c}^{3})=450\,\mathrm{Myr}$ is the Eddington time. Here we follow Pacucci et al. (2017) in accounting for the idea that accretion (especially of the super-Eddington variety) may be episodic, by identifying the duty cycle ${ \mathcal D }$ as the fraction of time spent accreting, and the Eddington ratio ${f}_{\mathrm{Edd}}$ that is =1 for Eddington-limited accretion but can be <1 for sub-Eddington accretion and >1 for super-Eddington accretion. We use a generalized accretion rate γ = χγ₀, where $\chi ={ \mathcal D }{f}_{\mathrm{Edd}}$ is a correction factor to account for the fact that the accretion rate could be super-Eddington for some periods of time. The quantitatively relevant parameter is $\chi ={ \mathcal D }{f}_{\mathrm{Edd}}$. Pacucci et al. (2017) found that objects with M ≳ 10⁵ ${M}_{\odot }$ can have high-efficiency accretion, $0.5\leqslant { \mathcal D }\leqslant 1$ and $1\leqslant {f}_{\mathrm{Edd}}\leqslant 100$, but objects with M ≲ 10⁵ M_⊙ have low-efficiency accretion, $0\leqslant { \mathcal D }\leqslant 0.5$ and $0\leqslant {f}_{\mathrm{Edd}}\leqslant 1$. The simplest assumption is that χ = 1 for Eddington-limited growth, but our model allows for the putative super-Eddington growth in the DCBH formation era.

We assume that the distribution of initial black hole masses is lognormal; i.e., the differential number density per logarithmic mass bin is distributed normally:

$$\begin{eqnarray}&&\displaystyle \frac{{dn}}{d\mathrm{ln}{M}_{0}}=\displaystyle \frac{1}{\sqrt{2\pi }{\sigma }_{0}}\exp \left(-\displaystyle \frac{{(\mathrm{ln}{M}_{0}-{\mu }_{0})}^{2}}{2{\sigma }_{0}^{2}}\right).\end{eqnarray} \tag{ 3 }$$

Here μ₀ and σ₀ are the mean and standard deviation of the distribution of $\mathrm{ln}$ M₀, respectively. A lognormal distribution for the birth mass function of DCBH seeds is consistent with the results of Ferrara et al. (2014) for intermediate masses ($4.75\lt \mathrm{log}(M/{M}_{\odot })\lt 6.25$), and we fit those results with μ₀ = 11.7 (corresponding to a peak at $\mathrm{log}M/{M}_{\odot }=5.1$) and σ₀ = 1.0.

Because the growth law implies that

$$\begin{eqnarray}&&\mathrm{ln}M(t)=\mathrm{ln}{M}_{0}+\gamma \,t,\end{eqnarray} \tag{ 4 }$$

we can write the mass function at a later time as

$$\begin{eqnarray}&&\displaystyle \frac{{dn}}{d\mathrm{ln}M}=\displaystyle \frac{1}{\sqrt{2\pi }{\sigma }_{0}}\exp \left(-\displaystyle \frac{{(\mathrm{ln}M-{\mu }_{0}-\gamma t)}^{2}}{2{\sigma }_{0}^{2}}\right).\end{eqnarray} \tag{ 5 }$$

The accretion time t may not be a fixed constant that applies to all objects, therefore we can integrate over a function f(t) (which has units of inverse time) that describes the distribution of accretion times. In this case the final observed mass function f($\mathrm{ln}$ M) ≡ dn/d $\mathrm{ln}$ M is

$$\begin{eqnarray}&&{\int }_{0}^{T}\displaystyle \frac{1}{\sqrt{2\pi }{\sigma }_{0}}\exp \left(-\displaystyle \frac{{(\mathrm{ln}M-{\mu }_{0}-\gamma t^{\prime} )}^{2}}{2{\sigma }_{0}^{2}}\right)f(t^{\prime} )\,{dt}^{\prime} .\end{eqnarray} \tag{ 6 }$$

Here f(t') is a normalized distribution of accretion times t' and T is the maximum possible accretion time. The function f(t') is determined by considering the creation rate of black holes in the DCBH scenario. The number density n of black holes grows in a type of chain reaction (Dijkstra et al. 2014; Yue et al. 2014) with the instantaneous creation rate dn/dt = λ(t) n. The simplest case, where λ(t) = λ, has a constant value leads to pure exponential growth. If this growth continues from a time t = 0 when the first DCBH is created until a time T when the creation of all DCBHs is terminated, then each black hole that was created at time t has an accretion lifetime t' = T − t in the range [0, T]. The normalized distribution of accretion lifetimes t' is then

$$\begin{eqnarray}&&f(t^{\prime} )=\displaystyle \frac{\lambda \,\exp (-\lambda \,t^{\prime} )}{[1-\exp (-\lambda \,T)]}.\end{eqnarray} \tag{ 7 }$$

Using the indefinite integral identity

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \int \exp [-({{ax}}^{2}+{bx}+c)]{dx}\\ \quad =\displaystyle \frac{1}{2}\sqrt{\displaystyle \frac{\pi }{a}}\exp \left(\displaystyle \frac{{b}^{2}-4{ac}}{4a}\right)\mathrm{erf}\left(\sqrt{a}\left[x+\displaystyle \frac{b}{2a}\right]\right),\end{array}\end{eqnarray} \tag{ 8 }$$

valid for a > 0, we evaluate the integral in Equation (6) using Equation (7) and obtain a full expression

$$\begin{eqnarray}\begin{array}{rcl}f(\mathrm{ln}M) & = & \displaystyle \frac{\alpha \,\exp \,(\alpha {\mu }_{0}+{\alpha }^{2}{\sigma }_{0}^{2}/2)}{2[1-\exp \,(-\alpha \beta )]}\,{M}^{-\alpha }\\ & & \times \left[\mathrm{erf}\left(\displaystyle \frac{1}{\sqrt{2}}\left(\alpha {\sigma }_{0}-\displaystyle \frac{\mathrm{ln}M-{\mu }_{0}-\beta }{{\sigma }_{0}}\right)\right)\right.\\ & & -\left.\mathrm{erf}\left(\displaystyle \frac{1}{\sqrt{2}}\left(\alpha {\sigma }_{0}-\displaystyle \frac{\mathrm{ln}M-{\mu }_{0}}{{\sigma }_{0}}\right)\right)\right].\end{array}\end{eqnarray} \tag{ 9 }$$

Here α ≡ λ/γ, the dimensionless ratio of the growth rate of DCBH formation to the growth rate of the mass of individual DCBHs, and β ≡ γ T, the dimensionless number of DCBH growth times within the DCBH formation era.

In the limit $T\to \infty$, the function becomes

$$\begin{eqnarray}\begin{array}{rcl}f(\mathrm{ln}M) & = & \displaystyle \frac{\alpha }{2}\exp \,(\alpha {\mu }_{0}+{\alpha }^{2}{\sigma }_{0}^{2}/2)\,{M}^{-\alpha }\\ & & \times \mathrm{erfc}\left(\displaystyle \frac{1}{\sqrt{2}}\left(\alpha {\sigma }_{0}-\displaystyle \frac{\mathrm{ln}M-{\mu }_{0}}{{\sigma }_{0}}\right)\right),\end{array}\end{eqnarray} \tag{ 10 }$$

which is the modified lognormal power-law (MLP) function (Basu et al. 2015). Equation (9) represents a tapered version of the MLP, with the break in the power law occurring at $\mathrm{log}M\approx ({\mu }_{0}+\beta )/2.3$, meaning that the peak of the original lognormal is shifted in $\mathrm{ln}$ M by an amount β = γT. Henceforth, we refer to Equation (9) as the tapered power-law (TPL) function.

Figure 2 shows the TPL function using parameter values obtained from models of the DCBH growth era. We pick μ₀ = 11.7 and σ₀ = 1.0 based on the model of Ferrara et al. (2014). From Dijkstra et al. (2014; see Figure 1) we adopt λ = 28.0 Gyr⁻¹ for the era of rapid DCBH formation using their canonical model. The length of the DCBH growth era is T = 0.15 Gyr (Yue et al. 2014). For accretion growth during this period we expect that super-Eddington growth can occur (Dijkstra et al. 2014) but with a wide range of possible values. We pick a series of values χ = [1, 2, 3] covering Eddington-limited growth and two values of super-Eddington growth. As γ = χγ₀ = 20χ Gyr⁻¹, this leads to α = [1.4, 0.7, 0.47] and β = [3, 6, 9] for our adopted values of λ and T. Figure 2 shows that the super-Eddington growth models allow for the development of a mass function that has both a visually evident power-law profile as well as a notable break in the power law at high mass. This break is a marker of the end of the DCBH growth era, as both the creation of new DCBH as well as their super-Eddington growth ceases after the time interval T.

Zoom In Zoom Out Reset image size

Once the DCBH growth era has ended at z ≈ 13, the population of DCBHs may continue to undergo Eddington-limited accretion, and the luminosity function can be estimated using Equation (1). Over time, the mass function f($\mathrm{ln}$ M) will retain its shape but move to the right because $\mathrm{ln}$ M at the end of the DCBH era will shift by an amount γΔt_z, where Δt_z is the time interval between the end of the DCBH era (z ≈ 13) and an observable redshift z. However, the duty cycle ${ \mathcal D }$, and therefore $\chi ={ \mathcal D }{f}_{\mathrm{Edd}}$, may be ≪1 in this era, rendering mass growth to small fractional levels. A random sampling of ${ \mathcal D }$ and ${f}_{\mathrm{Edd}}$ for individual object growth after z ≈ 13 shows that the overall distribution maintains its shape and moves to the right in $\mathrm{log}M$. We also note that the mass growth of SMBHs may be quenched above ∼10¹⁰ ${M}_{\odot }$ (Inayoshi & Haiman 2016; Ichikawa & Inayoshi 2017), in agreement with results of current quasar surveys (Ghisellini et al. 2010; Trakhtenbrot 2014).

Assuming that observed quasars are undergoing Eddington-limited accretion, we use Equation (1) to transform the mass function into a luminosity function. We expect that the mass of the quasars are not growing substantially during this time, for reasons discussed above, and we are really most interested in fitting the shape of the function, which should remain much the same for a variety of redshifts in the post-DCBH-growth era. Figure 3 shows the inferred quasar luminosity function (QLF) $\phi (L)\propto {dn}/d\mathrm{log}L$ for a suitable pair of values for α and β, overlaid on z = 3 quasar bolometric luminosities compiled by Hopkins et al. (2007). Here we are only interested in fitting the shape of the luminosity function and not the absolute number of sources. The normalization can be scaled to fit the observed number at any redshift.

Zoom In Zoom Out Reset image size

We hold (μ₀, σ₀) fixed at their model-inspired values (11.7, 1.0) as they are most important in determining the unobserved low end of the luminosity function. We effectively fit observations with two parameters (α, β). This is in contrast to the usual practice of fitting the observed QLF with a double power law (e.g., Hopkins et al. 2007; Masters et al. 2012; Schindler et al. 2019) that requires three parameters: the two power-law indices and a joining point.

In our model, the values of α and β that fit the QLF are not merely mathematical parameters. Instead, they reveal the history of the putative DCBH growth era. For the adopted DCBH number growth rate λ = 28.0 Gyr⁻¹ and a duration T = 0.15 Gyr, and individual masses growing at a rate γ = χγ₀, the two fitted parameters are related to the super-Eddington factor χ by

$$\begin{eqnarray}&&\alpha =1.4{\chi }^{-1},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&\beta =3\chi .\end{eqnarray} \tag{ 12 }$$

We find an excellent fit to the QLF with [α, β] = [0.5, 8.4]. Both the values of α and β imply a super-Eddington factor χ = 2.8, revealing the self-consistency of our model. In principle the QLF could have been fit with any α and β that could individually imply very different values of χ. In that case our underlying model would be inconsistent, or at least need to explore values of λ and T that were quite different than those implied by current models of the DCBH growth era.

To elaborate on the above point, our model could in principle also be applied to SMBH formation from alternate scenarios such as Population III remnants (Madau & Rees 2001; Whalen & Fryer 2012) or mergers of primordial Population III stars (Boekholt et al. 2018; Reinoso et al. 2018). It could apply as long as the black hole production could be described as growing exponentially at some rate λ and for a finite time T, during which the individual masses grew at an Eddington-limited or super-Eddington rate.

We have presented an analytic model that captures some essential features of the DCBH growth scenario and uses them to derive an analytic mass function and by implication a luminosity function. A double power-law function has been commonly used in the literature to mathematically fit the QLF. Here, we instead use a physically motivated formula based on the scenario of the DCBH growth era that has been developed by many researchers. It is not a double power law at high mass and luminosity, but rather a TPL. We believe that the rapid fall off in the QLF at high luminosity is better modeled as a tapered part of a power law rather than as a second power law. The break point of the power law identifies the end of the era of DCBH creation.

We have fit an observed QLF with a power-law index α = 0.5 and the break-point-related parameter β = 8.4. These are consistent with a period of rapid mass growth of DCBH with super-Eddington factor χ = 2.8, for a time period T = 150 Myr during which the growth rate of the number density ${n}_{\mathrm{DCBH}}$ was λ = 28.0 Gyr⁻¹. In principle, the best-fit values to QLF data can be used to constrain such theoretical models of DCBH growth.

Our model has two key components. Initially, high mass ∼10⁵ ${M}_{\odot }$ seeds grow rapidly in number during a limited time period in the early universe, because DCBH formation leads to the emission of LW photons that seed the formation of other DCBH. These objects also live within gas-rich halos and undergo super-Eddington mass accretion. Then at some time both the formation of DCBH, as well as the super-Eddington accretion of the existing DCBH, comes to a rapid halt due to the photoevaporation of the host halos. What remains is a TPL distribution of masses and therefore also of luminosity if the observed quasars are undergoing subsequent Eddington-limited accretion. Future modeling can relax some of these assumptions, for example the formation of DCBH may continue long enough to outlive the period of rapid (super-Eddington) mass growth, especially if driven by mechanisms other than the LW flux (Wise et al. 2019), and the super-Eddington accretion may not apply to all objects (Pacucci et al. 2017; Latif et al. 2018).

We thank the referee for constructive comments. S.B. was supported by a Discovery Grant from NSERC.