ACCRETION ONTO INTERMEDIATE-MASS BLACK HOLES IN DENSE PROTOGALACTIC CLOUDS

Known stellar-mass and intermediate-mass black hole candidates accrete through mass transfer from a stellar companion without exception. Rapidly accreting massive black hole candidates, such as in Seyfert galaxies and quasars, appear to be situated close to dusty, molecular, and highly inhomogeneous gas reservoirs, where cold, dense clumps are in free fall or are embedded in a hot, ionized, and tenuous confining medium. In neither case does the accretion flow resemble the laminar flow of the Bondi–Hoyle–Littleton solution. There is a distinct possibility, however, that massive black hole progenitor "seeds," with masses M_bh ≳ 100 M_☉ (Madau & Rees 2001), did pass through an early phase of quasi-radial accretion in protogalaxies, prior to significant dust enrichment and the emergence of a pervasive cold phase in the interstellar medium (e.g., Volonteri & Rees 2005; Alvarez et al. 2006, 2008; Johnson & Bromm 2007; Pelupessy et al. 2007; Di Matteo et al. 2008; Greif et al. 2008). This was the period immediately following the onset of atomic cooling, the epoch when protogalaxies experienced rapid growth due to cosmic infall and yet had potential wells too shallow to retain all of their newly synthesized metals.

Since the accretion rate in quasi-radial, Bondi-like accretion is proportional to gas density, Bondi-like accretion presents an appealing formation process for massive black holes: given the presence of sufficiently dense gas in a protogalactic nucleus, Eddington-limited accretion and mass exponentiation on a Salpeter timescale would proceed; stellar-size seed black holes would reach masses ≳10⁹ M_☉ inferred in high-redshift quasars in the cosmic time available to them (e.g., Haiman & Loeb 2001; Volonteri et al. 2003; Tanaka & Haiman 2008). A potential flaw in this reasoning was noticed early, namely, that in the event of radiatively efficient accretion, radiative heating and radiation pressure may reduce the accretion rate (Shvartsman 1971; Buff & McCray 1974; Hatchett et al. 1976; Ostriker et al. 1976; Cowie et al. 1978; Stellingwerf 1982; Begelman 1985; Ricotti et al. 2008; Alvarez et al. 2008). The possibility of a transition to radiative-feedback-driven temporal intermittency, which has been evident in one-dimensional simulations of Compton-heated (Ciotti & Ostriker 2001) as well as Compton- and photoionization-heated and radiation pressure-driven (Sazonov et al. 2005; Ciotti & Ostriker 2007) accretion onto massive black holes in elliptical galaxies, has further complicated the assessment of quasi-radial accretion's viability as a massive black hole formation process. If the accretion is indeed intermittent and episodic, how long is the duty cycle, how large is the average accretion rate, and does it permit evolution of seed black holes into quasars?  Does episodic accretion proceed in bursts, where episodes of Eddington-limited rapid accretion alternate with near-complete radiative quenching and outflow, as one-dimensional models appear to suggest, or does the multidimensional nature of realistic irradiated accretion flows imply a damping of the fluctuations and convergence to a quasi-steady state characterized by a reduced accretion rate?

Here, we present the first results from our program to investigate radiatively efficient accretion onto black holes embedded in a realistic galactic gaseous medium with the help of high-resolution hydrodynamic simulations. We do not attempt to simulate fluid flow and radiation emission mechanisms near the event horizon of the black hole; instead, we simply fix the shape of the spectral energy distribution of the radiation produced near the event horizon and simulate the irradiated fluid flow on radial scales spanning the sonic radius, where the photoheated fluid falls supersonically toward the black hole, and a much larger radius in the gas cloud that is entirely unaffected by the black hole's radiative output, where boundary conditions for the accretion flow are set. This work is organized as follows. In Section 2, we concisely describe our numerical algorithm. A more detailed description and a suite of tests will be provided in a longer follow-up paper (S. M. Couch et al. 2009, in preparation). In Section 3, we describe our results with particular focus on the intermittency of the accretion flow. In Section 4, we compare our results to analytical models (Milosavljević et al. 2009, and references therein) and discuss the outlook for massive black hole formation through quasi-radial accretion.

Here we describe numerical methodologies that we have implemented in the parallel adaptive-mesh-refinement (AMR) code FLASH (Fryxell et al. 2000), version 2.5, to facilitate our simulations of low-angular momentum accretion flow with a central ionizing source. The hydrodynamic solver used was the piecewise parabolic method (PPM), which is a higher order Godunov method. Contact steepening that sharpens contact discontinuities in the PPM module was disabled after we noticed that it led to spurious flow near stationary ionization fronts. By default, FLASH does not conserve angular momentum in two-dimensional cylindrical simulations. To simulate fluid motion in the azimuthal direction, we introduce the specific angular momentum l ≡ Rv_ϕ as an additional conserved "mass scalar" quantity that is advected with the fluid. Equivalence with three-dimensional cylindrically symmetric equations of fluid dynamics is achieved by adding the centrifugal force to the total acceleration, ${\bf a}={\bf a}_{\rm grav} + {\bf a}_{\rm rad} + l^2\hat{\bf R}/R^3$.

2.1. The Central Hole and the Central Ionizing Source

The simulations were carried out in cylindrical symmetry on a rectangular (R, z) grid in two spatial dimensions, where R is the distance from the axis of symmetry. In what follows, r will denote the distance from the origin, r² = R² + z². A small spherical region r < r_hole was excluded from the hydrodynamical update; we refer to this region as the "hole." A unidirectional "outflow" boundary condition, allowing only outflow from the simulated spatial domain and inflow into the central hole, was imposed at r = r_hole. The gravitational acceleration was set to that of a Newtonian point mass located at the origin of the grid ${\bf a}_{\rm grav}=-G M_{\rm hole}\hat{\bf r}/ r^2$. The mass of the material flowing into the hole was added to the point mass M_hole. Self-gravity of the fluid was ignored because, typically, the mass of the photoionized gas around the black hole was less than 1% of the mass of the black hole.

An isotropic source of ionizing radiation was placed at the origin of the grid. The total luminosity of the ionizing source between energies hν_min = 13.6eV and hν_max = 100keV was set to be proportional to the accretion rate into the hole,

$$\begin{equation} \int _{\nu _{\rm min}}^{\nu _{\rm max}} L_\nu d\nu = \epsilon \dot{M}_{\rm hole} c^2, \end{equation} \tag{ 1 }$$

where is a radiative efficiency which was uniformly set to 0.1. The radiative spectral energy distribution was set to the power law L_ν ∝ ν^−3/2 and was calculated in 40 logarithmically spaced energy bins.

2.2. Composition, Chemistry, and Thermodynamics

2.3. Line-of-Sight Radiative Transfer and Radiation Pressure

Here, we discuss our implementation of multifrequency photon-conserving line-of-sight radiative transfer and the resulting radiation pressure force. We calculate the line-of-sight optical depth τ^(k)_ν from the central source to the faces of any rectangular zone along N_ray = 4 separate rays, k = 1, ..., N_ray; the rays are equally spaced in angle as seen from the central source. We also calculate the optical depth to absorption within the zone along each of the rays $\hat{\tau }_\nu ^{(k)}$. The effective flux at the zone center is calculated via (see, e.g., Mellema et al. 2006, and references therein)

$$\begin{equation} F_\nu = \frac{L_\nu }{4\pi r^2} \frac{\sum _{k=1}^{N_{\rm rays}} \exp \big(-\tau _\nu ^{(k)} \big) \big[1-\exp \big(-\hat{\tau }_\nu ^{(k)} \big)\big]}{\sum _{k=1}^{N_{\rm rays}} \hat{\tau }_\nu ^{(k)}}. \end{equation} \tag{ 2 }$$

The total photon absorption rates per chemical species i within the zone are then calculated via Γ_i,ν = σ_i,νn_iVF_ν/hν, where V is the zone volume. In the limit in which the zone is very optically thick, $\hat{\tau }_\nu ^{(k)}\gg 1$, we obtain

$$\begin{equation} \sum _i \Gamma _{i,\nu } \approx (h\nu)^{-1} F_\nu \big\langle \exp \big(-\tau _\nu ^{(k)}\big)\big\rangle \Omega, \end{equation} \tag{ 3 }$$

where Ω is the angular size of the zone as seen from the central source, and the extinction external to the zone is averaged over rays. The approximation in Equation (3) becomes exact in the limit N_rays → ∞, as it has to be to conserve photons.

We calculate radiation pressure associated with photoionization and Thomson scattering. The radiation pressure force is added to the gravitational force of the central point mass. We ignore line resonance radiation altogether and also ignore any diffusion of ionizing radiation (see, e.g., Ritzerveld 2005).

2.4. Mesh Refinement and Initial Conditions

The central accretion rate, shown in Figure 1, oscillates between values close to, and occasionally mildly exceeding the Eddington limit, $\dot{M}_{\rm max}\sim \dot{M}_{\rm Edd}\sim 2\times 10^{-6}\;M_\odot \hbox{ yr}^{-1}$, and a rate lower by an order of magnitude, $\dot{M}_{\rm min}\sim 2\times 10^{-7} \;M_\odot \hbox{ yr}^{-1}$. The accretion is approximately periodic with a mild trend toward an increasing period separating consecutive peak episodes; the period varies between 250yr and 350yr. Accretion rate falloff following a maximum is roughly exponential, as may be expected when photoheating and photoionization radiation pressure in the ionized region surrounding the black hole drive down the accretion rate (see Section 3 in Milosavljević et al. 2009). The average accretion rate and luminosity are $\langle \dot{M}\rangle = 4.6\times 10^{19} \hbox{ g s}^{-1}$ and 〈L〉 = 4.2 × 10³⁹ergs⁻¹, which says that on average, the black hole accretes at 32% of the Eddington limit. The average accretion rate is still only ∼0.2% of the isothermal "Bondi" accretion rate $\dot{M}_{\rm Bondi} =e^{3/2}\pi (GM_{\rm bh})^2 n m_p/c_{\rm s}^{3}(\infty)$, calculated ignoring radiative feedback, for an ambient sound speed of c_s(∞) = 14kms⁻¹.

Zoom In Zoom Out Reset image size

During a central accretion maximum, electron scattering and photoionization radiation pressure drive an outflow in the ionized gas within the H ii region that has neutral fractions χ_H ∼ 10⁻⁴to10⁻⁵ (Figure 2, lower panels). This leads to rarefaction and exponential drop in central accretion. Meanwhile, as radiation pressure subsides, gas near the edge of the H ii region accelerates inward. This acceleration is driven by a gas pressure imbalance near the edge; the imbalance was inherited from the preceding accretion maximum when an outward radiation pressure force balanced an inward gas pressure gradient force. The outflow intersects with the inflow, and the inflowing gas ultimately arrives at the edge of the central hole and gives rise to a new accretion maximum. The longest timescale in the cycle is the inward acceleration and infall time, which is here a factor of ∼3 shorter than the radial sound crossing time of the H ii region with typical temperature T_H II ∼ 4 × 10⁴K.

Zoom In Zoom Out Reset image size

To investigate the relative importance of radiation pressure compared to heating, we continued the simulation shown in Figures 1–3 with the radiation pressure force artificially set to zero. The heating-only simulation exhibited a significantly higher, super-Eddington mean accretion rate and a rapid, monthly variability reflecting an episodic heating to the Compton temperature and convective transport at r ∼ r_hole. The longer period (∼300yr) behavior observed in the simulation with the radiation pressure was absent in the heating-only simulation. The very high accretion rate observed in the heating-only simulation is a potential artifact of r_hole being much larger than the sonic radius in this simulation, but if we spatially resolved the sonic radius in the high-accretion-rate regime, the Courant-limited time step would have become prohibitively short.

Zoom In Zoom Out Reset image size

Figure 3(a) shows that the radius of the H ii region surrounding the black hole reaches an asymptotic value of r_ion ∼ 6 ×  10¹⁶cm. This only somewhat smaller, by 46%, than the crude estimate obtained for hydrogen-only cloud in Milosavljević et al. (2009). The region r < r_ion is not entirely free of neutral gas; during accretion minima, as seen in the two left panels of Figure 2, dense clumps form in the gas that has been expelled by radiation pressure and is now falling back toward the black hole, and the regions in the dense clumps' shadows recombine. The range of radii that is shaded in Figure 3(a) contains a multiphase medium. The densest neutral clumps that form near the edges of the H ii region are six times denser than the ambient gas, but this density is not sufficient to resist photoevaporation and dispersal during accretion maxima. Figure 3(b) shows that the flow into the central hole is typically supersonic with sonic radius located at r_s ∼ 2 × 10¹⁴cm, but during accretion maxima, shock heating in the infalling gas drives subsonic flow to the brink of the excised hole at r_hole = 10¹⁴cm.

We have simulated accretion from a uniform, high-density (10⁷cm⁻³)¹ metal-free protogalactic cloud onto an intermediate-mass black hole. At radii from the black hole smaller than the size of the excised region, 10¹⁴cm, the accretion was assumed to be radiatively efficient. We find that a compact H ii region forms around the black hole. The accretion is intermittent due to alternating radiation pressure-driven expulsion and external pressure-driven fallback. The average accretion rate is ≲1% of the Bondi accretion rate calculated ignoring the radiation's influence, and ${\sim } \frac{1}{3}$ of the Eddington-limited rate.

The numerical results are consistent with the predictions of our toy model for episodic quasi-radial accretion (Milosavljević et al. 2009), where we suggested that photoionization heating and photoionization radiation pressure within the compact H ii region are the dominant accretion-suppression mechanisms. We also suggested that Lyα line resonance radiation pressure, which we do not account for in the simulations presented here, should further reduce the accretion rate. The intermittency is qualitatively similar to that seen in one-dimensional simulations of the accretion of hot irradiated gas onto black holes in elliptical galaxies by Ciotti & Ostriker (2007), where Compton heating and photoionization heating were both observed to be significant drivers of the accretion intermittency; we emphasize, however, that the gas densities, temperatures, and gravitational potential depths in elliptical galaxies are very different than the ones considered here.

The simulation in which the accreting fluid was endowed with small angular momentum corresponding to maximum circularization radius equal to r_hole exhibited similar intermittency and only a slightly reduced average accretion rate, consistent with Krumholz et al. (2005, and references therein). In contrast with the nonrotating simulations that did not contain axial outflows, narrow polar outflows form that resemble some of those found by Proga (2007) and Proga et al. (2008) in simulations of accretion flows irradiated by a quasar. The accretion intermittency and inflow–outflow behavior that we observe is qualitatively similar to that seen in Proga et al., but additional tests are needed to ascertain whether our axial outflows are physical.

The overall goal is to combine the highly resolved simulations presented here with large-scale cosmological simulations of how the first galaxies formed at the end of the cosmic dark ages (e.g., Greif et al. 2008). Those simulations can resolve the central gas flows to within ∼0.01pc, where we can implement subgrid prescriptions based on the calculations done here. An important challenge will be to derive prescriptions for the feedback from black hole accretion that are applicable to a wide range of protogalactic conditions. We will report on these efforts elsewhere (S. M. Couch et al. 2009, in preparation).

We acknowledge valuable discussions with Daniel Proga. The software used in this work was in part developed by the DOE-supported ASC/Alliance Center for Astrophysical Thermonuclear Flashes at the University of Chicago. The authors acknowledge the Texas Advanced Computing Center (TACC) at the University of Texas at Austin for providing high-performance computing resources that have contributed to this research. V.B. and M.M. acknowledge support from NSF grant AST-0708795.