BLACK HOLE FORMATION IN PRIMORDIAL GALAXIES: CHEMICAL AND RADIATIVE CONDITIONS

Supermassive black holes with masses M>10⁸–10⁹ M_☉ are known to have existed at very early times in the history of the universe (Jiang et al. 2009), prompting one to ask how black holes of such a size could have formed so quickly. They may have started life as stellar-mass objects, the compact remnants of the first generation of stars, and they may have grown to their present size by accretion. However, forming supermassive black holes in this way presents a number of difficulties. The first stars are thought to have had masses on the order of 100 M_☉ or more, making them strong sources of ionizing radiation (Abel et al. 2002; Bromm & Larson 2004; Glover 2005; Yoshida et al. 2008). The dark matter halos in which these stars formed were small, with masses of only 10⁶ M_☉, and models of the effects of the ionizing radiation produced by a primordial star within one of these small halos show that it readily photo-evaporates the bulk of the gas from the halo, leaving little gas available to be accreted by the black hole (Johnson & Bromm 2007; Milosavljevic et al. 2009; Alvarez et al. 2009). Although accretion is likely more efficient in the 5–6σ peaks that harbor the observed supermassive black holes, it is still constrained by the Eddington limit. As discussed by Shapiro (2005), seed masses of at least 10⁵ M_☉ are required for a magnetohydrodynamic disk or a standard thin disk. Similarly, Schleicher et al. (2010) found that high Eddington ratios and unusually low mass–energy conversion efficiencies are needed for stellar progenitors.

Previous studies have therefore examined the possibility of forming these large seed black holes by direct gravitational collapse (Eisenstein & Loeb 1995; Koushiappas et al. 2004; Begelman et al. 2006; Spaans & Silk 2006). Gas falling into a halo with a mass M>5 × 10⁷[(1 + z)/10]^−3/2 M_☉ (corresponding to a virial temperature T_vir>10⁴ K) will be shock-heated up to T_vir, and so will become hot enough to cool via the electronic emission lines of atomic hydrogen. The fate of the gas in one of these so-called "atomic cooling" halos will then depend on its subsequent thermal evolution. If the gas can cool efficiently as it collapses, significantly lowering both its temperature and its Jeans mass, then it is very difficult to prevent it from fragmenting, and hence forming stars rather than a massive black hole (Bromm & Loeb 2003; Omukai et al. 2008; Clark et al. 2009). On the other hand, if cooling remains inefficient and the effective equation of state for the gas is stiffer than isothermal, then fragmentation is suppressed (Li et al. 2005; Lodato & Natarajan 2006) and the formation of a massive black hole is a plausible outcome. In gas that has been enriched with heavy elements and dust from the first generation of stars, it appears that cooling and fragmentation cannot be avoided (Omukai et al. 2008). On the other hand, in primordial gas the only effective low-temperature coolant is molecular hydrogen (H₂), and if enough of this can be destroyed by photodissociation, then cooling can be suppressed (Begelman et al. 2006; Dijkstra et al. 2008; Regan & Haehnelt 2009; Shang et al. 2010). The soft-UV background needs to be strong enough to suppress H₂ formation until densities of ∼10⁵ cm⁻³ are reached, following which collisional dissociation will suppress the H₂ abundance for gas temperatures greater than a few 1000 K. The required background field depends on the chemical model, the spectral shape, and the potential presence of additional electrons created in shock fronts. It is on the order of J₂₁ = 10³–10⁵ (Omukai 2001; Bromm & Loeb 2003; Shang et al. 2010). These values are much larger than estimates of the mean strength of the background, which range from J₂₁ ∼ 0.1 (Johnson et al. 2008) to J₂₁ ∼ 40 (Dijkstra et al. 2008). However, recent work has shown that the radiation background is highly inhomogeneous, owing to the strong clustering of the first generation of star-forming galaxies (Dijkstra et al. 2008; Ahn et al. 2009). In regions close to these galaxies, the ultraviolet background can be very large, and Dijkstra et al. (2008) have recently demonstrated that enough massive primordial galaxies with a local value of J₂₁>10³ may exist. In addition, early reionization in local patches of the universe may suppress star formation in some primordial halos in a similar fashion (Johnson 2009). They may thus remain pristine, but grow in mass until they can collapse.

In the absence of H₂ cooling, emission from the Lyman series lines of atomic hydrogen, especially Lyα, plays a central role in regulating the temperature of the collapsing gas. Hydrodynamic models of the collapsing gas typically assume that the Lyman series lines remain effectively optically thin, allowing the optically thin form of the cooling rate to be used. These models find that the temperature of the collapsing gas decreases with increasing density, albeit very slowly, typically reaching T ∼ 6000–7000 K at densities on the order of n ∼ 10⁹ cm^-3. However, the large hydrogen column densities present in these protogalaxies produce very large optical depths in the Lyman series lines, and so the assumption that the optically thin form of the cooling rate can safely be used is highly questionable. Spaans & Silk (2006) examined the effects of Lyα photon trapping within an analytic framework for the evolution of the gas and showed that in the absence of additional cooling mechanisms, it gives rise to a temperature evolution described by

$$\begin{equation} \gamma _{\rm eff} =1+\frac{d\log T}{d\log \rho }\sim 1-\frac{0.5-7/18\cdot Bn^{7/18}}{\log C n^{1/2}+Bn^{7/18}}, \end{equation} \tag{ 1 }$$

where ρ is the mass density, T the temperature, n the number density, B ∼ 0.5–0.1 cm^7/6, and C ∼ 10⁻³⁶M_h/(10⁷ M_☉) cm^3/2, with M_h the halo mass. For halo masses in the range 10⁷–10⁹ M_☉, this yields γ_eff − 1 ∼ 0.01–0.5, indicating an effective equation of state that is stiffer than isothermal, and hence a temperature that increases with increasing density. This result suggests that the neglect of the effects of photon trapping is not only a poor approximation, but also may lead to qualitatively incorrect results, as it is crucial for fragmentation whether the effective equation of state is harder than isothermal.

In this paper, we re-examine the issue of Lyα photon trapping, using a considerably more detailed treatment than in Spaans & Silk (2006). In particular, we address the issue of whether there are other cooling mechanisms that can compensate for Lyα if the latter is strongly suppressed.

To follow the thermal and chemical evolution of the collapsing gas, we make use of a one-zone treatment in which the form of the density evolution is prescribed in advance. One-zone models are widely used for studying the chemistry and thermodynamics of primordial or low metallicity gas (e.g., Omukai 2001; Omukai et al. 2005; Cazaux & Spaans 2009; Glover & Savin 2009; Schleicher et al. 2009), as they enable one to model the chemistry in great detail and to include a large number of different cooling processes, without the computational efficiency concerns inherent to a three-dimensional treatment.

We extend the one-zone model developed by Glover & Savin (2009) and Schleicher et al. (2009) to include effects that become important at higher temperatures and in the presence of Lyα trapping. Near T ∼ 10⁴ K, H⁻ formation cooling may contribute to the overall cooling rate (Omukai 2001). We assume that a typical electron undergoing radiative attachment to form H⁻ has an energy on the order of k_BT, where k_B is Boltzmann's constant, and hence write the H⁻ formation cooling rate as

$$\begin{equation} \Lambda _{\rm H^-}=k_{\rm H^-} n_{\rm H} n_{\rm e} k_{\rm B} T, \end{equation} \tag{ 2 }$$

where $k_{\rm H^-}$ is the H⁻ formation rate coefficient, n_H is the atomic hydrogen number density, and n_e is the electron number density. To model the effects of Lyα trapping, we include different level populations as separate species in the code. In our model, we consider energy levels up to n = 5. For the first excited state, we distinguish between the 2s and the 2p states, as only decays from the 2p state to the ground state will produce Lyα photons; two-photon decays from the 2s state will produce continuum photons that will not be trapped. The other states are considered as averages over the angular momentum quantum numbers. Radiative decay rates A_ij, collisional excitation and de-excitation rates C_ij, collisional ionization and three-body recombination rates, and collisional transition rates between the 2s and 2p state are all adopted from Omukai (2001). We define the transition rate R_ij from level i to level j as

$$\begin{eqnarray} R_{ij}&=&A_{ij}\beta _{{\rm {\scriptsize esc}},ij}(1+Q_{ij})+C_{ij},\quad i>j, \end{eqnarray}\vspace*{-12pt} \tag{ 3 }$$

$$\begin{eqnarray} R_{ij}&=&\frac{g_j}{g_i}A_{ji}\beta _{{\rm {\scriptsize esc}},ji}Q_{ji}+C_{ij}, \quad i<j, \end{eqnarray} \tag{ 4 }$$

where $\beta _{{\rm {\scriptsize esc}},ij}$ is the escape probability for the i → j transition, g_i = 2i² denotes the statistical weight of level i and $Q_{ij}=c^2 J_{{\rm {\scriptsize cont}},ij}/(2h\nu _{ij}^3)$, where ν_ij is the frequency of the transition and $J_{{\rm {\scriptsize cont}},ij}$ the average intensity of the background radiation field at this frequency. The spectral shape of the background radiation field is taken to be that of a blackbody with a temperature of 10⁵ K (Glover & Savin 2009), and the normalization is set by fixing the value of J₂₁. We assume that UV photons more energetic than 13.6 eV are absorbed by atomic hydrogen in the intergalactic medium, and we do not consider the effects of an extragalactic X-ray background. For radiative or collisional transitions to the first excited state, we assume that the reaction products will be distributed according to the statistical weights of the 2s and 2p states. However, our results are not sensitive to this assumption, as collisional transitions between these two states occur rapidly in the conditions of interest, and so the ratio between the 2s and 2p level populations is always very close to equilibrium.

For the escape probability $\beta _{{\rm {\scriptsize esc}},ij}$, we adopt the expression

$$\begin{equation} \beta _{{\rm {\scriptsize esc}},ij}=\frac{1-\mathrm{exp}(-\tau _{ij})}{\tau _{ij}}\mathrm{exp}(-\beta t_{{\rm {\scriptsize ph}}} / t_{{\rm {\scriptsize coll}}}), \end{equation} \tag{ 5 }$$

where τ_ij is the optical depth at line center of the i → j transition, $t_{{\rm {\scriptsize coll}}}$ the collapse time of the gas, and $t_{{\rm {\scriptsize ph}}}$ is the photon diffusion time, i.e., the time required for a photon to diffuse out of the optically thick gas, and β a geometrical factor. Even small deviations from spherical symmetry lead to escape along a preferred direction and β = 3 (Dijkstra et al. 2006; Spaans & Silk 2006). In our model, the collapse time is the free-fall time, corrected by a factor that takes into account the effective equation of state (see Omukai et al. 2005; Schleicher et al. 2009). For most line transitions, $t_{{\rm {\scriptsize ph}}} \ll t_{{\rm {\scriptsize coll}}}$, such that the exponential factor is negligible. For these lines, we follow Omukai (2000) and assume that the dominant contribution to the optical depth comes from material within one local Jeans length, and that the density, temperature, etc., do not vary significantly on this scale. This is justified provided that the collapse is close to isothermal, and that the initial mass of the collapsing gas is comparable to the Jeans mass, as in this case the density and velocity profiles of the gas will come to resemble the Larson–Penston similarity solution (Larson 1969; Penston 1969), which has just these properties. For Lyα photons and other direct transitions to the ground state, we compute the optical depth τ and the photon diffusion timescale following Spaans & Silk (2006), in order to take into account the significant line broadening that occurs during the diffusion process. For the Lyα line, this yields a timescale $t_{{\rm {\scriptsize ph}}}=L (a\tau _{21})^{1/3}/c$ with the natural-to-thermal line width a and an optical depth τ₂₁ = 1.04 × 10⁻¹³N_HT^−0.5₄, with L the distance to the edge of the halo, and N_H the column density of atomic hydrogen. For the Lyα line, the natural-to-thermal line width is taken as a = 4.7 × 10⁻⁴T^−1/2₄ (Spaans & Silk 2006), while we correct this value for the increased lifetime in case of other line transitions to the ground state. To calculate the column density, we assume that the density profile of the gas scales with radial distance r as r^−2.2, as indicated by previous numerical simulations (e.g., Wise & Abel 2007). The cooling functions for the hydrogen lines are then evaluated based on the escape probability, the level populations and the strength of the background radiation.

The temperature evolution that we find with our approach for several different values of J₂₁ is plotted as a function of density in Figure 1. Even for J₂₁ = 1, the temperature remains much higher during the collapse than in the radiation-free case, and never drops significantly below 1000 K.

Zoom In Zoom Out Reset image size

As previously noted, much higher values of J₂₁ can be obtained locally, in the presence of a luminous neighbor within a distance of 10 kpc or less (Dijkstra et al. 2008). As we increase J₂₁, we find systematically higher temperatures at each density; and for J₂₁ = 100, we find that the effects of H₂ collisional dissociation start to become important. In this case, H₂ cooling is marginally effective at densities n < 10⁸ cm^-3, but at higher densities, collisional dissociation destroys the H₂, causing a sharp rise in the temperature. Such a radiation field may be present in a fraction of 10⁻⁵ − 10⁻³ of all atomic cooling halos (Dijkstra et al. 2008). For J₂₁ = 10³ and J₂₁ = 10⁴, we find, in common with previous studies, that H₂ cooling never becomes important, and the gas temperature remains high. Nevertheless, in contrast to the predictions of Spaans & Silk (2006), we find that the temperature does decline with increasing density, although it never drops below 5000 K over the range of densities examined here. With γ_eff = 0.97–0.98, we find that the effective equation of state of the gas is just slightly softer than isothermal. We note, though, that this estimate provides just a lower limit for the value of γ_eff, which is strictly defined as dlog T/dlog ρ at constant entropy.

In order to understand why the gas is still able to cool, we explore the case with J₂₁ = 1000 in more detail in Figures 2 and 3. In Figure 2, we show the evolution of the fractional abundances of free electrons and H₂ molecules, as well as the fractional level populations of the 2s, 2p, and 3 states of atomic hydrogen (denoted henceforth as H(2s), H(2p), and H(3)). In Figure 3, we show the contributions that the main cooling processes make to the total cooling rate. The Lyα optical depth depends on the size of the protogalaxy, and we consider two examples, one with a virial radius of 500 pc (corresponding to a halo mass of ∼10⁷ M_☉), and a second with a virial radius of 3 kpc (corresponding to a halo mass of 2 × 10⁹ M_☉). In each case, the abundance of the 2p state is initially very small, but increases roughly linearly with increasing density, owing to the increasing importance of Lyα photon trapping. However, collisional excitation to the H(3) state eventually comes to dominate over radiative decay from the 2p state, following which the abundance of the 2p state remains roughly constant.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The abundances of the H(2s), H(2p), and H(3) states are clearly larger in the case with the larger virial radius, as Lyα photons are trapped more effectively, directly boosting the population of the 2p state. The population of the 2s state is then increased due to collisional coupling. The same holds true for the H(3) state. However, the electron and H₂ abundances are not affected. Despite these changes in the hydrogen level populations, the total cooling rate does not change significantly. Cooling due to the 2p–1s line decreases, but this is balanced by additional cooling from the 2s–1s and 3–2 transitions, due to the increased populations of these states. At densities higher than 10⁸ cm⁻³, collisional dissociation of H₂ becomes the dominant cooling process (even after correcting for H₂ formation heating) and cools the gas to about 5000 K. H⁻ formation cooling is found to contribute significantly in a broad range of densities.

Although we find that the effective equation of state of the gas may be softer than isothermal, it is nevertheless true that γ_eff ∼ 1 throughout the collapse. Moreover, the temperature evolution of the gas is surprisingly similar to that obtained in previous three-dimensional simulations that assume optically thin Lyα cooling. The outcome of the initial collapse may thus be similar to that found in previous simulations, i.e., the formation of a single massive bound object at the center of the halo (see, e.g., Wise et al. 2008; Begelman & Shlosman 2009; Regan & Haehnelt 2009; Shang et al. 2010). However, it is also highly likely that gas that falls in later, with higher angular momentum, will begin to buildup an accretion disk surrounding this object. The high gas temperatures imply a rapid accretion flow onto such a disk, which may quickly become gravitationally unstable. If the gas in the disk is also able to cool effectively, which our results suggest will be the case even if Lyα cooling is completely discounted, then we would expect it to fragment and form stars (Levin 2007). The latter may be drawn to the center by dynamical friction and merge with the central black hole (Devecchi & Volonteri 2009). Owing to the high temperatures, we would expect the stars formed in the disk to accrete more rapidly than standard Population III stars, by a factor of 10 or more, and hence their final masses may be very large, M ∼ 100–500 M_☉. However, recent works also indicate the possibility of quasi-stable cold self-gravitating accretion disks that do not fragment, which could feed the central object during the further evolution (Levine et al. 2008; Begelman & Shlosman 2009).

The research leading to these results has received funding from the European Community's Seventh Framework Programme (/FP7/2007-2013/) under grant agreement No 229517. D.R.G.S. and S.C.O.G. are grateful for subsidies from the Landesstiftung Baden-Württemberg via their program International Collaboration II. We thank the referee Nick Gnedin for helpful comments that improved our manuscript.