High-redshift SMBHs can grow from stellar-mass seeds via chaotic accretion

We test the evolution of black hole mass and spin, as well as the final SMBH mass distributions, under various accretion alignment conditions, using a numerical model. The model is initialised with a seed black hole of mass $M_{0}=10{\,\rm M_{\odot}}$, with a random spin magnitude and direction. A series of accretion events are then used to grow the SMBH. Each event is limited in mass by the self-gravity condition (eq. 3); we also checked that imposing an additional maximum mass of an individual event, to simulate the possible mass limit of gas clouds feeding the SMBH, has no effect on our results. The initial direction of the disc’s spin is chosen completely randomly in one model; in two others, we add an additional parameter $p_{\rm aligned}=0.2;0.4$, and randomly select this fraction of accretion events to have spins aligned in the positive-z direction. The code then checks the condition for hole-disc alignment (eq. 4) and calculates the new direction of the hole’s angular momentum by considering the vector sum of $\vec{J_{\rm d}}$ and $\vec{J_{\rm h}}$. The mass and angular momentum of the black hole are updated by adding the appropriate fractions of the disc mass (eq. 3) and angular momentum (eq. 5). The new angular momentum is used to calculate the spin parameter and the radiative efficiency to be used for the next accretion event. The duration of a single episode is calculated as the time for the AGN to consume the disc mass at the given luminosity:

We set the luminosity to be equal to the Eddington luminosity of the black hole. While the luminosity is likely to vary a lot during SMBH growth, we are mainly interested in the early Universe and in extreme black holes which are likely to have grown in dense gas environments where prolonged feeding at the Eddington (or higher) rate was common. If the black hole is fed at a mildly super-Eddington rate, the accretion disc sheds material at radii larger than ISCO and the black hole does not grow faster than at the Eddington rate (Shakura & Sunyaev, 1973; Pringle, 1981).

The code runs, updating the SMBH mass and spin and advancing the time, until either $700$ Myr have passed or the SMBH mass exceeds $10^{11}{\,\rm M_{\odot}}$. We implicitly assume that the SMBH mass and spin do not change in between accretion events, so the time tracked in the model is actually the “accretion time”, defined over sufficiently long timescales as

where $t$ is actual time elapsed. We run 150 simulations in total, 50 each with $p_{\rm aligned}=0$, $0.2$ and $0.4$. The only differences between the simulations in each set are the random numbers used to define the initial SMBH spin magnitude and direction and the directions of individual accretion episodes.

The time coordinate in our models is the accretion time rather than total time of the black hole’s existence. So the fact that $t_{\rm acc}\sim 520$ Myr is required for production of a Pōniuā’ena -mass SMBH means that the duty cycle of its growth has to be almost $100\%$. This may be unrealistic combined with the necessity of having a low rate of alignment of accretion directions; in fact, one might expect accretion directions to be highly correlated if accretion is essentially continuous. Nevertheless, there are several ways of mitigating this tension.

The masses of individual accretion episodes in our model are limited by the self-gravity condition (Pringle, 1981; King & Pringle, 2007). The validity of this condition rests on two assumptions: that the disc is geometrically thin and optically thick, i.e. that the standard thin-disc equations are applicable to describing its structure, and that cooling is efficient at radii $R\lower 2.15277pt\hbox{$\;\buildrel>\over{\sim}\;$}10^{3}r_{\rm s}$, where $r_{\rm s}=2GM_{\rm BH}/c^{2}$ is the Schwarzschild radius of the SMBH. Thin discs are the standard solution describing an accretion flow with mass flow rate $0.01\dot{M}_{\rm Edd}\lower 2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}\dot{M}\lower 2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}2\dot{M}_{\rm Edd}$ (Best & Heckman, 2012; Sądowski et al., 2013; Inayoshi et al., 2016). The importance of self-gravity in the outer regions of these discs has been recognised by numerous authors (e.g., Shlosman et al., 1990; Gammie, 2001; Rice et al., 2005). In addition, cooling is typically efficient in these regions and leads to star formation (Shlosman & Begelman, 1989; Collin & Zahn, 1999; Goodman, 2003); even if the forming stars heat the disc up and prevent further fragmentation, this does not shut off the growth of these stars (Nayakshin, 2006). The stellar rings around Sgr A^∗​, the central SMBH of the Milky Way, provide strong evidence for the validity of this scenario (Paumard et al., 2006; Bonnell & Rice, 2008; Bartko et al., 2009). Cooling in primordial accretion discs is probably less efficient due to low metallicity, but this is unlikely to stabilize them significantly against self-gravity (Tanaka & Omukai, 2014).

The gas reservoir surrounding the SMBH is likely to be turbulent, perhaps composed of streams and orbiting clumps. In this case, the streams of gas approaching the SMBH may easily have different directions, even if the reservoir as a whole has some angular momentum. The radius of the feeding reservoir immediately available to the SMBH is $\sim R_{\rm infl}$, the radius of influence of the black hole’s gravity (Dehnen & King, 2013), which is much higher than the maximum radius of a non-self-gravitating accretion disc. As a result, even a small difference in the impact parameter of two streams falling toward the SMBH will lead to significant differences in the angle of the forming disc (Nayakshin et al., 2012). This suggests that even essentially continuous accretion can be almost entirely chaotic.

Streams falling in from various directions also continuously change the orientation of the accretion disc. This happens on approximately the dynamical timescale of the reservoir,

where we scale the reservoir size to 100 pc and the velocity dispersion to 200 km s^-1. This should have little effect on the initial alignment between the disc and the hole, because only the warped region of the disc is affected by this. So long as $M_{\rm BH}\lower 2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}10^{8}{\,\rm M_{\odot}}$, the warp radius is smaller than the self-gravity radius and the initial alignment between the black hole and the disc is not affected by the outer regions of the disc.

Some galaxies have disky structures surrounding their central SMBHs on scales of order $0.1-1$ pc, e.g. NGC 4258 (Herrnstein et al., 2005), IC 2560 (Yamauchi et al., 2012) and NGC 1068 (Imanishi et al., 2020). Although these are almost certainly strongly self-gravitating systems (e.g., Kumar, 1999; Lodato & Bertin, 2003; Martin, 2008; Tilak et al., 2008), we can consider, for the sake of argument, what happens if such systems are prevented from fragmenting. Over time, the accretion disc aligns with the black hole in a prograde fashion. However, this this happens on the viscous timescale at the radius where the disc angular momentum starts to dominate over that of the hole, which is (cf. King et al., 2005)

where $M_{8}$ is the black hole mass scaled to $10^{8}\,{\,\rm M_{\odot}}$ and assuming Eddington-rate accretion. The ratio between the two timescales is

This ratio is much lower than unity except for the smallest black holes, and even then only if the feeding reservoir is as large as 100 pc. If the black hole is initially fed inside a cloud with radius $R_{\rm cl}\sim 24$ pc, the dynamical timescale is shorter than the alignment timescale even for a $10\,{\,\rm M_{\odot}}$ seed black hole. As a result, new gas streams effectively randomize the direction of the accretion disc much faster than the disc can cause the SMBH to align with it.

The reservoir is depleted by accretion and small-scale momentum-driven feedback, but outflows are unlikely to affect dense gas outside of the sphere of influence for as long as the SMBH mass is lower than the critical $M_{\sigma}$ mass (King, 2010). The depleted reservoir is continuously refilled by gas falling in from larger distances throughout the host halo. Assuming that the halo is gas-rich, as appropriate for high redshift, the infall rate may be a significant fraction of the dynamical rate (cf. Dehnen & King, 2013)

where we scale the gas fraction $f_{\rm g}$ to the cosmological value $0.16$ and the velocity dispersion $\sigma$ to $200$ km s^-1. This rate is much higher than the Eddington accretion rate even for a SMBH with the mass equivalent to Pōniuā’ena ​. Considering that the most massive black holes should reside in haloes with the highest velocity dispersion, refilling the reservoir surrounding the black hole is even easier.

The seed mass we chose, $M_{0}=10{\,\rm M_{\odot}}$, is also rather conservative. It is comparable to that expected for supernova remnants of present-day stars. Early, Population III, stars are expected to be more massive on average than present-day ones (Hirano et al., 2014; Susa et al., 2014; Stacy et al., 2016); as a result, their explosion remnants might be more massive as well. The final mass of the black hole at any given time is approximately proportional to the seed mass, so a more massive seed makes it slightly easier to grow an extremely massive SMBH quickly. Furthermore, the effective seed mass might be higher if the black hole forms in a cluster and merges with a few other similar black holes while its mass is still small. However, the constraint on SMBH spin, and hence on the maximum level of alignment, hardly changes if we assume higher seed masses - this can be seen in Figure 1.

Our model is insensitive to the precise redshift of the formation of the black hole seed. In principle, even a seed forming as late as at $z\sim 18$, corresponding to $t=180$ Myr after the Big Bang, can give rise to Pōniuā’ena ​, provided that it maintains a duty cycle $\delta=1$, so that $t_{\rm acc}=t$. If the seed forms at $t=100$ Myr, a duty cycle of $520/600\mathrel{\raise 1.29167pt\hbox{$\sim$}\mkern-14.0mu\lower 1.72218pt\hbox{$-$}}0.87$ is enough.

The possibilities described in sections 5.1 - 5.3 allow for some variation of model parameters, but the observational constraints are still tight. They also provide a stringent test for our model. In particular, the model would be falsified if an SMBH would be detected at a time $t_{\rm obs}$ with a mass greater than

This corresponds to $M=10^{9}{\,\rm M_{\odot}}$ at $t=545$ Myr after the Big Bang, or, equivalently, $z\mathrel{\raise 1.29167pt\hbox{$\sim$}\mkern-14.0mu\lower 1.72218pt\hbox{$-$}}9.1$. The growth timescale of $27.6$ Myr corresponds to the purely chaotic SMBH growth model, and we use the large seed mass $M_{0}=100{\,\rm M_{\odot}}$ for this estimate. Quasars powered by such SMBHs - if they exist - should be discovered by upcoming Euclid (Roche et al., 2012) and Athena (Nandra et al., 2013; Aird et al., 2013) missions. Discoveries of slightly less extreme SMBHs, including that of Pōniuā’ena ​, provide constraints on parameter combinations encompassing the seed black hole mass, duty cycle, formation time and alignment fraction of accretion episodes.

Current observational estimates of the space density of extreme high-redshift quasars with $L>10^{46}$ erg s^-1 are $n_{\rm QSO}\lower 2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}10^{-7}$ Mpc^-3 (Kulkarni et al., 2019; Wang et al., 2019; Shen et al., 2020). This luminosity corresponds to the Eddington limit for a $7.7\times 10^{7}{\,\rm M_{\odot}}$ SMBH. Quasars with bolometric luminosity $L\sim 2\times 10^{47}$ erg s^-1, corresponding to the luminosity of Pōniuā’ena (Yang et al., 2020), are even rarer: $n_{\rm Pon}\lower 2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}10^{-8}$ Mpc^-3. Given that we predict an almost 100% duty cycle for these objects, the space density of SMBHs should be not much higher than that of quasars. These values allow us to estimate the rarity of conditions that give rise to such extreme black holes so quickly.

The rate of formation of the first stars is estimated at $\dot{M}_{*}\mathrel{\raise 1.29167pt\hbox{$\sim$}\mkern-14.0mu\lower 1.72218pt\hbox{$-$}}2-3\times 10^{-5}{\,\rm M_{\odot}}$ yr^-1 Mpc^-3 at $z=20$ (Johnson et al., 2008); this drops down to $\dot{M}_{*}\mathrel{\raise 1.29167pt\hbox{$\sim$}\mkern-14.0mu\lower 1.72218pt\hbox{$-$}}10^{-6}{\,\rm M_{\odot}}$ yr^-1 Mpc^-3 at $z=30$. These Population III stars were probably more massive, on average, than stars forming today, with an average mass of $\sim 10-100{\,\rm M_{\odot}}$ (Hirano et al., 2014; Susa et al., 2014; Stacy et al., 2016). Approximately half of the stars may end up forming black holes at the end of their lives (Heger & Woosley, 2002; Hirano et al., 2014). This leads to the formation rate of seed black holes of $\dot{N}_{\rm seed}\mathrel{\raise 1.29167pt\hbox{$\sim$}\mkern-14.0mu\lower 1.72218pt\hbox{$-$}}5\times 10^{-9}-5\times 10^{-8}$ yr^-1 Mpc^-3 at $z=30$ and $\dot{N}_{\rm seed}\mathrel{\raise 1.29167pt\hbox{$\sim$}\mkern-14.0mu\lower 1.72218pt\hbox{$-$}}10^{-7}-1.5\times 10^{-6}$ yr^-1 Mpc^-3 at $z=20$. Assuming that the seeds of the most extreme quasars must have been produced before $z=18$ (corresponding to $t=180$ Myr, leaving $520$ Myr for Pōniuā’ena to grow), this gives $t_{\rm form}\sim 80$ Myr for seed formation. Over this period, even the most conservative estimate provides a seed BH density $n_{\rm seed}\sim 0.4$ Mpc^-3, and it is possible that this density was at least an order of magnitude higher. These values are consistent with estimates based on numerical modelling of Population III star formation (Madau & Rees, 2001; Volonteri, 2010). Comparing this value with the quasar space density, we see that if only one in $4\times 10^{7}$ of initial seeds encounters favourable growth conditions, the population of extreme quasars observed at $z\sim 6-7$ is reproduced.

This result suggests an explanation for the issue of ‘stellar black hole starvation’, which is expected at high redshifts (e.g. Pacucci et al., 2015, 2017; Smith et al., 2018). It happens due to two reasons: many nascent stellar black holes can be ejected from their parent star-forming halos (Whalen & Fryer, 2012), and those that remain cannot accrete material efficiently until the fossil HII region that surrounded the progenitor star cools down (Johnson et al., 2007). However, a small fraction of black holes can avoid these issues if they are born within a larger protogalactic halo with escape velocity $v_{\rm esc}>100$ km s^-1. A natal kick with $v<v_{\rm esc}$ allows the black hole to escape the fossil HII region very quickly. Small-scale feedback can enhance gas mixing (Dehnen & King, 2013) and suppress fragmentation (Alvarez et al., 2009) around these black holes, promoting accretion at near-Eddington rates. The expected rarity of these objects can also explain why they are not reproduced in numerical and semi-analytical simulations (Smith et al., 2018; Griffin et al., 2019).

We expect that the key parameters determining the growth rate of the SMBH vary smoothly over their allowed range; furthermore, there are no discontinuities or sharp changes in the dependence of the final SMBH mass on these parameters. Due to this, we predict that the mass function of black holes, and of quasars, should be approximately continuous in mass/luminosity, consistent with available observational constraints. However, the slope of the mass function may be very shallow, i.e. the density of less massive black holes might not be much higher than that of the most massive ones. The total X-ray radiation produced by high-redshift quasars suggests that luminous accretion contributed only $<1000{\,\rm M_{\odot}}$ Mpc^-3 to the SMBH mass density by redshift $z=6$ (Treister et al., 2013). The observed quasars already account for $\rho\mathrel{\raise 1.29167pt\hbox{$\sim$}\mkern-14.0mu\lower 1.72218pt\hbox{$-$}}\langle M_{\rm BH}\rangle n_{\rm QSO}\mathrel{\raise 1.29167pt\hbox{$\sim$}\mkern-14.0mu\lower 1.72218pt\hbox{$-$}}100{\,\rm M_{\odot}}$ Mpc^-3, i.e. $>10\%$ of this limit. This means that the smaller black holes cannot have grown very much at these early times. In addition, the lower duty cycle of these smaller black holes means that the quasar density is significantly lower than the SMBH density. These two factors might explain why the smaller SMBHs at high redshift have not been detected yet.

Our results provide two important and potentially testable predictions regarding the properties of extreme high-redshift SMBHs: their spins must be low and they should reside in dense gas environments.

Typical spins of black holes growing through purely chaotic accretion settle to $a<0.3$ in $t_{\rm acc}<350$ Myr. At this time, the black hole mass is only $M\sim 10^{6}{\,\rm M_{\odot}}$. Similarly, if accretion is slightly aligned with $p_{\rm aligned}=0.2$, black hole spin is $a<0.5$ when its mass is $10^{6}{\,\rm M_{\odot}}$ or higher. If accretion is aligned more, the black hole does not reach a mass of $10^{8}{\,\rm M_{\odot}}$ over $600$ Myr of accretion. So the most massive black holes, those with $M>10^{8}{\,\rm M_{\odot}}$, observed when the Universe was $\sim 700$ Myr old or younger, should all have spins lower than $a=0.5$. Among more moderate black holes, with $10^{6}{\,\rm M_{\odot}}<M<10^{8}{\,\rm M_{\odot}}$, the spin distribution should be wide, with both low and high spin values possible. Black holes with masses $M<10^{6}msun$ should generally spin rapidly, although low spin values should be possible in this population as well. This contrasts with available observational constraints on SMBHs in the Local Universe, where black holes with masses $M\lower 2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}3\times 10^{7}{\,\rm M_{\odot}}$ seem to have spin values close to maximal, and even more massive black holes generally have spins $a>0.4$ (Vasudevan et al., 2016; Reynolds, 2019, 2020).

There are a few ways of reconciling our prediction with the observationally-estimated spins. First of all, the observed high-redshift quasars represent an extreme population (see Section 5.5) and their spins need not be representative of the high-redshift population as a whole. We predict that as SMBHs with more modest masses are discovered at high redshift, and as their spin estimates are obtained, the overall population will come to resemble the locally-observed one more closely, even if not exactly. Secondly, the global mass growth of supermassive black holes is expected to be dominated by mergers since $z\sim 0.5$, i.e. for the past 3-4 Gyr (Bassini et al., 2019). Mergers tend to flatten the spin distribution (Berti & Volonteri, 2008), so even if the extreme quasars have spins representative of the early SMBH population as a whole, the local population will have much more broadly distributed spin values. There is tentative observational evidence, based on the evolution AGN jet power, that the typical spins have been increasing at least since $z=1$ (Martínez-Sansigre & Rawlings, 2011). Finally, observational biases lead to easier identification of AGN powered by rapidly spinning SMBHs, and the currently observed population is consistent with a non-negligible number of slowly-spinning black holes (Vasudevan et al., 2016). The situation should become significantly clearer once the so-far elusive low mass SMBH population at high redshifts will be discovered, with upcoming observatories such as Lynx (Gaskin et al., 2019) and Roman (Green et al., 2012). In addition, future missions, including Athena (Barret & Cappi, 2019) and AXIS (Mushotzky et al., 2019), will be able to measure SMBH spins with unprecedented precision even at high redshift, providing better constraints on the evolution of this quantity.

As spin determines radiative efficiency, information about quasar luminosity might also be used to put constraints on black hole spin. The average radiative efficiency over the whole growth period is $\langle\epsilon\rangle\mathrel{\raise 1.29167pt\hbox{$\sim$}\mkern-14.0mu\lower 1.72218pt\hbox{$-$}}0.058$ in the purely chaotic models and increases to $\langle\epsilon\rangle\mathrel{\raise 1.29167pt\hbox{$\sim$}\mkern-14.0mu\lower 1.72218pt\hbox{$-$}}0.095$ in the $p_{\rm aligned}=0.4$ models. In order to power the observed luminosity of Pōniuā’ena ​, $L_{\rm Pon}\mathrel{\raise 1.29167pt\hbox{$\sim$}\mkern-14.0mu\lower 1.72218pt\hbox{$-$}}1.9\times 10^{47}$ erg s^-1, the SMBH should be accreting material at a rate $\dot{M}_{\rm Pon}\mathrel{\raise 1.29167pt\hbox{$\sim$}\mkern-14.0mu\lower 1.72218pt\hbox{$-$}}3.34/\epsilon\,{\,\rm M_{\odot}}$ yr^-1. This translates to $57.64\,{\,\rm M_{\odot}}$ yr^-1 and $35.2\,{\,\rm M_{\odot}}$ yr^-1 for the two alignment cases, respectively. Assuming this represents the rate of spherical infall of fully ionised material at the outer edge of the accretion disc, $r_{\rm out}\mathrel{\raise 1.29167pt\hbox{$\sim$}\mkern-14.0mu\lower 1.72218pt\hbox{$-$}}0.01$ pc, the expected column density is

where $\mu\mathrel{\raise 1.29167pt\hbox{$\sim$}\mkern-14.0mu\lower 1.72218pt\hbox{$-$}}0.5$ is the mean atomic weight per particle and $m_{\rm p}$ is the proton mass. The two cases of low- and high-spin SMBH give expected hydrogen column density around the Compton-thick limit, $N_{\rm H}=1.5\times 10^{24}$ cm^-2. Future observations might be able to measure the X-ray spectrum of these AGN and determine their Compton thickness, providing constraints on matter infall rates and, hence, the radiative efficiency of accretion.

Our models are consistent with the general expectation of numerical simulations that the most massive SMBHs should be found in dense environments (Overzier et al., 2009; Costa et al., 2014; Habouzit et al., 2019). Some observational campaigns confirm this (e.g., Kim et al., 2009; Husband et al., 2013), while others don’t (Morselli et al., 2014; Mazzucchelli et al., 2017). The discrepancy in results may arise due to various selection effects, low-number statistics and stochasticity (Buchner et al., 2019). In addition, the host galaxies of high-redshift quasars tend to have very high star formation rates (SFR; Trakhtenbrot et al., 2017; Decarli et al., 2018), indicating high gas content. The reason for this correlation in our model is the requirement of a massive gas reservoir that can feed the SMBH while withstanding AGN feedback; in some cases, this feedback may be positive and induce a part of the observed high SFR. In contrast, other SMBH origin models require the seeds to form in the largest dark matter overdensities.

Some recent works have investigated the evolution of SMBH spin and mass from first principles, both at high redshift or during the whole age of the Universe.

Sesana et al. (2014) used a semi-analytical model of galaxy evolution to track the mass and spin evolution of black holes from high-redshift seeds to the present day. They use several prescriptions linking the large-scale dynamics of the galaxy to that of the accretion flow on to the SMBH. In all of them, the flow inherits some memory of the angular momentum of the whole system, leading to partially aligned accretion. Predictably, then, the calculated SMBH spin is generally high, especially at high redshift, when gas accretion is the dominant mode of mass growth. The authors did not set out to reproduce the observed extreme quasars, but the few examples of individual SMBH mass and spin evolution show very little growth in the first Gyr after the Big bang. At later times, spin tends to decline somewhat, coming into agreement with available low-redshift observational constraints.

Griffin et al. (2019) use the Galform semi-analytical galaxy evolution model to track the mass and spin evolution of black holes. They use stellar black hole seeds and include mergers and chaotic cold gas accretion as modes of evolution of SMBH spin. They recover the present-day SMBH mass function and find that SMBHs should generally have rather low spins $a<0.5$, but they do not reproduce the extreme high-redshift quasars. This is most likely because of limited sample size which only allows them to account for SMBH space densities $>10^{-7}$ Mpc^-3 at $z>6$.

Zhang & Lu (2019) used a model composed of two SMBH growth phases: an initial super-Eddington coherent phase followed by a number of Eddington-limited chaotic episodes. During the latter, the SMBH spin decreases, but the equilibrium values are typically higher than in our models, mainly due to the somewhat less stringent constraint on maximum disc mass, which leads to more common alignment of the disc and the black hole at the early stages of evolution. This led the authors to conclude, in (Zhang et al., 2020), that a period of super-Eddington growth is necessary to produce the observed high-redshift quasars. In that paper, they added an additional constraint that the currently-observed Eddington ratio of the high-redshift quasars is the maximum allowed in the chaotic accretion periods of SMBH growth, which further constrained the amount of mass obtained during chaotic accretion episodes. Given that individual episodes of SMBH growth are short compared to the total growth time and that the AGN luminosity can vary significantly during each episode, we do not think such a constraint is warranted.