Hyper-Eddington growth of IMBHs

Hyper-Eddington accretion flows onto massive black holes


We study very-high rate, spherically symmetric accretion flows onto massive black holes (BH; ) embedded in dense metal-poor clouds, performing one-dimensional radiation hydrodynamical simulations. We find solutions from outside the Bondi radius at hyper-Eddington rates, unimpeded by radiation feedback when , where and are the density and temperature of ambient gas. Accretion rates in this regime are steady, and larger than , where is the Eddington luminosity. At lower Bondi rates, the accretion is episodic due to radiative feedback and the average rate is below the Eddington rate. In the hyper-Eddington case, the solution consists of a radiation-dominated central core, where photon trapping due to electron scattering is important, and an accreting envelope which follows a Bondi profile with . When the emergent luminosity is limited to because of photon trapping, radiation from the central region does not affect the gas dynamics at larger scales. We apply our result to the rapid formation of massive BHs in protogalaxies with a virial temperature of . Once a seed BH forms at the center of the galaxy, it can grow to a maximum via gas accretion independent of the initial BH mass. Finally, we discuss possible observational signatures of rapidly accreting BHs with/without allowance for dust. We suggest that these systems could explain Ly emitters without X-rays and nearby luminous infrared sources with hot dust emission, respectively.

black hole physics, cosmology: theory, quasars: supermassive black holes

1 Introduction

Can very high accretion rates onto black holes (BHs) occur stably? The existence of bright quasars (QSOs) at high-redshift () provides a challenging puzzle about the origin of supermassive black holes (SMBHs) with masses (Fan, 2006; Willott et al., 2010; Mortlock et al., 2011; Venemans et al., 2013; Bañados et al., 2014; Wu et al., 2015). To form such massive objects within the first billion years after the Big Bang, rapid growth of seed BHs in the early Universe is required.

The initial seeds may originate via several scenarios (e.g. Volonteri, 2012; Haiman, 2013, references therein). One possibility is massive BHs formed through collapse of Population III (Pop III) stars with typical mass of (Yoshida, Omukai & Hernquist, 2008; Hosokawa et al., 2011; Stacy, Greif & Bromm, 2012; Hirano et al., 2014; Hosokawa et al., 2015), which grow up to SMBHs via continuous gas accretion or frequent major mergers of their host galaxies (Madau & Rees, 2001; Haiman & Loeb, 2001; Volonteri, Haardt & Madau, 2003; Li et al., 2007; Alexander & Natarajan, 2014; Kulier et al., 2015). A second possibility is the so-called direct collapse scenario, in which a supermassive star collapses by general relativistic instability and turns to a massive seed BH with mass of (Loeb & Rasio, 1994; Oh & Haiman, 2002; Bromm & Loeb, 2003; Lodato & Natarajan, 2006; Begelman, Volonteri & Rees, 2006; Regan, Johansson & Haehnelt, 2014; Inayoshi, Omukai & Tasker, 2014; Becerra et al., 2015; Latif, Schleicher & Hartwig, 2016). In this scenario, the constraint on BH growth time is alleviated because the initial BH mass is larger, but the duty cycle needs to be (Tanaka, 2014), which is difficult to achieve. As an alternative, a massive BH could form via runaway stellar collisions in a globular cluster (Sanders, 1970; Portegies Zwart et al., 2004; Freitag, Gürkan & Rasio, 2006) and in protogalaxies (Omukai, Schneider & Haiman, 2008; Devecchi & Volonteri, 2009; Katz, Sijacki & Haehnelt, 2015; Yajima & Khochfar, 2016). In any of these cases, seed BHs with still need to grow by subsequent accretion and mergers up to SMBHs with by (e.g. Tanaka & Haiman, 2009). The Soltan-Paczynski argument demands that most BH growth at occurs in the accretion regime (Soltan, 1982; Yu & Tremaine, 2002).

The basic physics of the accretion process was first studied for a spherically symmetric flow without radiative feedback (Bondi, 1952). Solutions with radiative cooling and heating have been studied by many authors (e.g. Shapiro, 1973; Park, 1990a, b; Nobili, Turolla & Zampieri, 1991; Ciotti, Ostriker & Proga, 2009, 2010). With only two-body radiative processes included, the solutions are characterized by the normalized luminosity and accretion rate,


where is the Eddington luminosity, is the electron scattering opacity, and (Chang & Ostriker, 1985)1. In Fig. 1, we show schematically the characteristic behavior of the luminosity and BH growth rate as a function of the Bondi accretion rate, respectively. For low , the accreting gas is optically thin to electron scattering, and radiation processes (cooling, heating and pressure) are inefficient (green region in Fig. 1). Therefore, the accretion flow results in an approximately adiabatic Bondi solution, where the accretion occurs within the Bondi radius , where is the sound speed of the ambient gas. The luminosity from the accreting gas is as low as (Shapiro, 1973; Park, 1990b).

On the other hand, for (blue region in Fig. 1), the hydrodynamical reaction to radiation processes has a major role in determining the properties of the accretion flow. Heating by hot Compton radiation with from the gas near the BH strongly suppresses the accretion from the Bondi radius. As a result, steady and self-consistent solutions exist only for and (Ostriker et al., 1976; Park, 1990a, b; Nobili, Turolla & Zampieri, 1991). However, these solutions have been found to be unstable, resulting in a highly time-dependent flow (e.g. Cowie, Ostriker & Stark, 1978; Zampieri, Miller & Turolla, 1996; Novak, Ostriker & Ciotti, 2012; Gan et al., 2014). We here examine a yet larger regime and explore the existence of self-consistent solutions with very high , where the Compton heating does not suppress the accretion, which we hereafter refer to as the hyper-Eddington regime.

The accreting gas in fact forms a tiny disk around the central BH, even with small angular momentum. Then, the gravitational energy is released via radiation through dissipative processes in the accretion disk more efficiently than for spherically symmetric flows. The typical efficiency is for a thin disk around a non-spining BH (Shakura & Sunyaev, 1973). Several works have studied accretion flows by adding the radiation from a central disk by hand to a Bondi-like solution, and found that the accretion flow does not approach a steady state but becomes episodic due to ionization and heating (Ciotti & Ostriker, 2001; Milosavljević et al., 2009; Milosavljević, Couch & Bromm, 2009; Novak, Ostriker & Ciotti, 2011; Park & Ricotti, 2011, 2012; Novak, Ostriker & Ciotti, 2012). In two-dimensional radiation hydrodynamical (RHD) simulations, these authors found that the time-averaged accretion rate onto a massive BH () is limited to the Eddington accretion rate (). Moreover, gas accretion from cosmological scales is strongly reduced by radiative heating, because of the shallow gravitational potential of dark matter (DM) halos hosting the BHs (Alvarez, Wise & Abel, 2009). Again, radiative heating can shut down gas supply from large scales and potentially quench the BH growth in this intermediate accretion rate domain (cf. blue region in Fig. 1).

In addition, the radiation force can affect the properties of accreting gas near the BH. For , the accretion luminosity released near the Schwarzschild radius () is estimated as , which naively exceeds the Eddington luminosity. In fact, however, the luminosity in a spherical inflow does not exceed because the photons are trapped and advected inward with the falling gas within a characteristic “trapping” radius before escaping via radiative diffusion (Begelman, 1978). Under steady-state conditions, the trapping radius is given by


For high , the luminosity released above () should thus be smaller than (Begelman, 1979). In the non-spherical case, RHD simulations (some of which include magnetic fields) suggest that (1) the gas can accrete with through a geometrically thick disk, where photon trapping occurs, and (2) the radiation flux is highly anisotropic and the luminosity can at least moderately exceed the Eddington luminosity in the polar direction (e.g. Ohsuga et al., 2005, 2009; Yang et al., 2014; Fragile, Olejar & Anninos, 2014; McKinney et al., 2014; Sa̧dowski et al., 2015; Takahashi & Ohsuga, 2015; McKinney, Dai & Avara, 2015). A recent three-dimensional simulation (Jiang, Stone & Davis, 2014) found that high- accretion is possible though magnetic-field advection of radiation in the disk makes photon trapping less efficient. However, all of the above multi-dimensional RHD simulations have been limited to the region near the central BH ().

In this paper, we investigate another domain and consider solutions with very high accretion rates (), examining self-consistent solutions of hyper-Eddington accretion from large scales () onto massive BHs (). We perform one-dimensional hydrodynamical simulations which include multi-frequency radiation transfer and non-equilibrium chemistry. We first run simulations of the outer region which resolves the Bondi radius () and find that hyper-Eddington accretion from the Bondi radius is realized without the negative effects of radiative feedback when (red region in Fig. 1). Consistent with other authors, we find that for , gas accretion leads to episodic inflow due to radiative feedback and the average accretion rate is limited to . We study several cases for different parameters (e.g. and ) and clarify the physics of the transition from the episodic state to a steady state. We find that this transition is well-explained by the comparison of the Bondi radius and the size of the ionized region (H region). Second, for the hyper-Eddington solutions, further simulations of the inner region, which resolve the trapping radius (), are conducted. We confirm that when the emergent luminosity is limited to because of photon trapping, radiation from the inner region does not affect the gas dynamics in the outer region. Moreover, we examine whether our solution would connect with a hyper-Eddington accretion disk solution smoothly at small radii ().

We briefly highlight the primary new results in our paper compared to previous works. For the outer-region simulations, our treatment is based on that of 2D-RHD simulations (Milosavljević, Couch & Bromm, 2009; Park & Ricotti, 2011, 2012), where the radiative efficiency is assumed to be constant () and the average accretion rate is . We update the model of the radiative efficiency including the photon trapping effect for , and find a new pathway to form the accretion solutions with a hyper-Eddington rate. Recently, Pacucci, Volonteri & Ferrara (2015) also studied solutions with high accretion-rates with an 1D-RHD simulation using a gray-approximation around the Bondi radius, and noted the possibility of super-Eddington accretion. Our new points are that combination of photon trapping and radiative recombination induces the transition to a steady state regime for and that the final steady state of the accreting gas is essentially an isothermal Bondi solution. Moreover, we study the gas properties in the inner region, which is not resolved by Pacucci, Volonteri & Ferrara (2015), by solving the radiation transfer equations which include the photon trapping effect self-consistently. With this addition, we show that our solutions may plausibly be connected to hyper-Eddington disk solutions at the central region near the BH.

The rest of this paper is organized as follows. In §2, we describe the methodology of our RHD simulation. In §3, we show the simulation results and discuss the necessary conditions required for a steady, hyper-Eddington accretion flow to exist. In §4, we present a discussion of the growth of seed BHs in the early Universe, and two other possible observational consequences and tests of our new solutions. In §5 and §6, we discuss caveats of our simulations and summarize the main conclusions of this paper.

2 Simulation method

2.1 Our strategy

We here use the hydrodynamical simulation code (ZEUS, Stone & Norman 1992) including multi-frequency radiation transfer, photoionization and heating, and primordial chemical network. We consider the simplest case, spherically symmetric accretion (e.g. Park & Ricotti, 2011). As a reference of the accretion rate, we define the Bondi rate for isothermal gas, i.e. the specific heat ratio ,


where is the ambient gas density. The sound speed is given by , where is the gas temperature, and is the mean molecular weight. Throughout this paper, the reference Bondi radius and Bondi accretion rate are estimated by using K, and . But both of and are solved for self-consistently in our simulations.

Figure 1: A schematic overview of BH accretion solutions. The characteristic behavior of the luminosity and BH growth rate are shown as a function of the Bondi accretion rate. In the lowest regime (green), the accretion flow results in an approximately adiabatic Bondi solution, where and . In the intermediate regime (blue), the accretion flow becomes unstable and highly time-dependent because of radiative feedback, where and on time average. In highest regime (red), a new self-consistent solution of steady hyper-Eddington accretion is found in this paper. In this new regime, the luminosity is limited to because of photon trapping, and a very-high accretion rate is realized (). The region between the intermittent and hyper-Eddington accretion is still uncertain because the hydrodynamical instability (e.g. RT instability) might decrease the critical accretion rate.

We need to consider gas dynamics over a wide range of spacial scales. Fig. 2 shows a schematic picture of the flow. In our case, the typical values of the Bondi and the trapping radius are


where , and . The required dynamical range of our simulation box is at least orders of magnitude in radius, which is computationally prohibited. Thus, we separate our simulation region into two regions of an outer region () and an inner region (), and first perform simulations of the outer region. For the cases with , we run several simulations by setting the gas properties (i.e. the density, thermal energy density, and velocity) at the inner-boundary of the outer-region simulation to the outer-boundary conditions of the inner-region simulation, following each domain appropriate time resolution. Finally, we obtain self-consistent solutions of the accretion flow onto a BH with hyper-Eddington accretion rates () by combining the results of the inner and outer regions. Note that since solutions with do not approach a steady state due to radiative feedback in the outer region, we neither conduct the inner-region simulations nor obtain a fully self-consistent solution in this unstable regime.

Figure 2: A schematic picture of a spherically symmetric accretion flow onto a massive BH at a hyper-Eddington accretion rate (). There are three characteristic scales: the Bondi radius (), photosphere (), and the trapping radius (). The dashed curve marks the boundary between the two regions simulated separately: the outer region () and the inner region ().

2.2 Basic equations

The basic equations of hydrodynamics we solve are the following: the equation of continuity


the equation of motion


where is the gas density, is the radial velocity (inflow; ), is the gas pressure, the gravitational potential with a general relativistic correction is set to (Paczyńsky & Wiita, 1980), and is the outward net radiation force in the radial direction.

We solve the energy equation including radiative cooling and heating,


where is the specific energy (erg g). The equation of state of the ideal gas is assumed as , where . The first term of the right-hand side is the compressional heating term. The last two terms are radiative cooling and heating, whose rates are and in units of erg s cm. The cooling rate is estimated as


where each term corresponds to the cooling rate associated with H, He, He atoms, free-free emission, and chemical reactions (see below). For the outer-region simulation, we assume optically thin cooling rates of H atoms (Ly, ), He atoms ( state) and He ions, and free-free transitions (Glover & Jappsen, 2007). In the inner region, since the gas is opaque to Ly photons, we solve the level population of H atoms ( and state) including the Ly trapping effect and estimate the cooling rate of two-photon emission (Omukai, 2001). In addition to the H transitions, free-bound emission of H () contributes as a cooling process. Thus, for the inner-region simulation, . We show the details of our treatment of the Ly trapping, continuum radiation cooling, and opacity in the Appendix A.

We estimate the cooling rate by solving a chemical reaction network of metal-free gas, which is composed of seven species (H, H, e, H, He, He, and He). Since the reactions relevant to H occur faster than the gas dynamical timescale, the H fraction is assumed to be in equilibrium (see Appendix A3). The chemical reactions include photoionization, collisional ionization, radiative recombination and collisional recombination of H, He and He. Instead of considering photoionization by diffuse photons, we adopt the on-the-spot approximation where the case A radiative recombination rate coefficient is replaced by that for case B. To provide a stable, positive definite and first-order accurate solution of the chemical network, we use a method based on a semi-implicit formulation (Anninos et al., 1997). The order of the updating is H, H, He, He, He and e (Whalen & Norman, 2006). For the inner-region simulation, inside the photosphere where all reactions are balanced, the chemical abundances are determined by solving the Saha equations instead of the non-equilibrium reaction network.

To ensure the accuracy of solutions of the hydrodynamical equations coupled with radiative cooling/heating and primordial chemistry, the time step must be shorter than the Courant time (the Courant number is set to 0.1), cooling/heating time and chemical time . The cooling/heating time and chemical reaction time are given by


where and are the electron and neutral fraction (Whalen & Norman, 2006, 2008). We set the time step to the lowest value among these timescales and integrate the energy equation (8) by an implicit method as well.

To solve the above basic equations, we employ spherical coordinates with a logarithmically-spaced grid in the radial direction: the position of -th grid is given by for , where is the radius of the inner boundary, is the size of the inner-most grid-cell, is the size ratio between consecutive grids, is the number of grids, and the radius of the outer boundary is given by . In our simulations, the size and layout of the coordinate grid is characterized by the four parameters (, , , and ). The number of the grid cells is for the outer (inner) region simulations. The value of is chosen to resolve the ionization front and the Bondi radius. Our fiducial value is . We check the convergence of our results, changing from to and .

2.3 Radiation transfer

To estimate the radiative heating rate , ionization rate , and radiation pressure force, we need to solve for radiation quantities in the accreting gas. We here assume a steady state radiation field because the light crossing time (, where is the optical depth) is much shorter than the time step of hydrodynamics outside the trapping radius. Thus, we have set the position of the innermost grid to in the inner-region simulation.  In what follows, we describe how to solve the radiation transfer (RT) equations in the two regions, respectively.

Outer region

For the outer-region simulation, we follow the same treatment as Park & Ricotti (2011, 2012). The RT equation is simply given by


where is the radiative flux, the radiation energy density, the emissivity, and the absorption opacity. Since the gas is transparent to photons with eV, we solve the RT equation for ionizing photons ( eV). Inside the ionized region, the gas is optically thin even to ionizing photons and thus . Therefore, the radiation transfer equation is reduced to a simpler form:


where the frequency range is eV keV.

We consider ionization and radiative heating due to bound-free absorption of H, He and He. The ionization and heating rates are estimated so that the number of photons emitted along any line of sight will always equal the number of photo-ionizations in that direction over any time interval (Whalen & Norman, 2006; Whalen et al., 2008) as


where is the mean intensity which conserves the number of ionizing photons, ( H, He, He) is the bound-free cross section and is the ionization threshold, is the photoelectron’s energy (). The acceleration due to the momentum transfer from radiation is


In the outer-region simulation, for simplicity, we set a radiation source with a single power-law spectrum


and ( eV) as the inner boundary condition. This assumption somewhat overestimates the effect of radiative feedback because this spectrum would be harder than that expected from self-consistent solutions (see below). We here neglect the secondary ionization and heating due to high-energy electrons with eV produced by the primary ionization. Since this assumption is not valid for the case with a power-law spectrum, we underestimate the effect of ionizing radiation. However, radiation from the inner region would have a thermal spectrum with , whose peak energy is eV, in the hyper-Eddington cases we are interested in (see §3.5). We here set the power-law index to .

The normalization of the luminosity at the inner boundary is estimated from the bolometric luminosity given by , where is the radiative efficiency. We here consider two models for : (1) standard model ( constant) and (2) trapping model. The standard model corresponds to radiation from a thin disk around the BH (Shakura & Sunyaev, 1973). For comparison with previous studies (e.g. Park & Ricotti, 2011, 2012), we consider several cases of . As the trapping model, we adopt a simple prescription for the efficiency given by


where for and for . Since the trapping effect becomes important for , we set a transition point around . For high , the luminosity has a maximum value of () because of the photon trapping effect. We note that the choice of the maximum luminosity is motivated by Begelman (1979), who argued that the emergent luminosity is limited to even in a hyper-Eddington accretion phase. As we discuss in §3, physical processes to provide our main conclusion are not radiation pressure but photo-ionization heating. Therefore, this choice does not change our main results significantly.

Inner region

In the inner region, a photosphere forms because of the high density and opacity. We define the position of the photosphere as


where and are the optical depth due to absorption (H free-bound and H free-free transition) and scattering (H Rayleigh and electron scattering) processes.

As the initial conditions, we set uniform profiles of the density, temperature, velocity and neutral H fraction (), whose values are given by those at the inner-boundary of the outer-region simulation. If we run simulations under these conditions, the density increases at the center before a sufficiently large H region forms. Then, photoionization and recombination are tightly balanced at the boundary of the small H region embedded by dense neutral gas. As a result, the chemical reaction time becomes much shorter than the gas dynamical time, and thus we need a long computational time to reach the final steady state. To avoid this problem, we adopt a value of the equilibrium opacity (Eqs. 71 and 72) instead of the non-equilibrium opacity at the early stage of the simulation (note that the chemical reaction network is solved outside the photosphere; see §2.2). Because the equilibrium opacity is higher than the non-equilibrium opacity, a photosphere forms in the H region before the density increases. Since the chemical abundance inside the photosphere is given by the Saha equations, this approach saves significant computation time. After one dynamical time at the outer-boundary ( s), we use the non-equilibrium opacity self-consistently and obtain the final steady state at s. We show the details about the treatment of opacities in Appendix A3.

Outside the photosphere (), we can use the same RT equation as in the outer region (Eq. 13) for ionizing photons with eV. For simplicity, we assume that the gas is optically thin against low-energy photons with eV except Ly photons ( eV). We here neglect radiative cooling via other Lyman series and Balmer series photons since their total cooling rate is comparable or smaller than that of two-photon emission (Omukai, 2001; Schleicher, Spaans & Glover, 2010; Johnson et al., 2012).

Inside the photosphere, radiation pressure dominates over the gas pressure and photon trapping is important. To capture these effects, we consider the following steady RT equations including the first-order correction in velocity (Mihalas & Klein, 1982; Lowrie, Morel & Hittinger, 1999; Jiang, Stone & Davis, 2014)


where all radiation moments are measured in the rest frame for an observer at infinity. The last two terms of the right-hand side in Eq. (20), , are the second-order corrections in . These terms were included by Lowrie, Morel & Hittinger (1999) to ensure the correct thermal equilibrium state in moving fluids when scattering dominates the opacity.

The Lorentz transformation between radiation moments in the rest frame () and those in the fluid frame () is given by , and . Inside the photosphere, the radiation is thermalized and isotropic. Thus, , where is the Planck function, and the radiative flux in the fluid frame is approximated in the diffusion limit as . Thus, the differences in and between the rest and fluid frame is negligible because . However, the difference in the radiative flux is not negligible inside the photosphere () because (Mihalas & Mihalas, 1984).

In the optically thick regime, we can approximate even in the rest frame. Thus, . Then, the RT equations in the rest frame are


We note that Eq. (23) is identical to the Lorentz transformation of the radiative flux. The first and second term corresponds to the diffusion and advection term, respectively. After the frequency integration, the ratio of these terms is


which means that the radiative flux in the rest frame has a positive (negative) value at ().

In order to solve the RT equations both inside and outside the photosphere, we need boundary conditions. The photospheric luminosity and the effective temperature are eigen-values to be determined by the requirement that integration inward from the photosphere reaches acceptable central conditions. The photospheric luminosity is determined by physical conditions inside the photosphere (e.g. disk formation and/or shocks). Begelman (1979) estimated the maximum value of the emergent luminosity from a spherical accretion flow with a high as , based on analytical arguments. Referring to this value, we produced several simulations for , and studied whether the accretion flows approach a steady state under these boundary conditions. Then, we determine the physically-correct value of such that our solution smoothly connects with a small central accretion disk, which would form within (see §3.5). Once setting and , the bolometric radiative flux is given by . At the photosphere, , where is the Stefan-Boltzmann constant. Assuming that the radiation has a thermal spectrum with , we obtain . Using this value, we solve the RT equation inside and outside the photosphere. Finally, for the gas at , we solve the temperature profile using


which is given by the frequency integration of Eq. (23), i.e. , where is the Rosseland mean opacity.

3 Results

3.1 Ordinary stellar-mass BH ()

In this section, we discuss the case of an ordinary stellar-remnant BH with . This corresponds to an accreting BH which forms through the gravitational collapse of a massive star. As initial conditions of the outer region, we adopt a neutral uniform gas with , and for comparison with previous works. The corresponding Bondi and Eddington accretion rates are and (), respectively.

Figure 3: Time evolution of the gas accretion rate onto a massive stellar-remnant BH ( and ; in the outer region). Red solid and blue dashed curves correspond to the cases of the standard model with and the trapping model with , respectively. Open circles mark the six epochs at which we show radial profiles in Fig. 4. The corresponding Bondi accretion rate is . Even in the trapping model, where the luminosity does not exceed , the time averaged accretion rate is limited to . This result demonstrates that gas accretion is quenched not by radiation pressure but by radiative heating.

Fig. 3 shows the time evolution of the accretion rate (in Eddington units) for the standard model (red) with , and the trapping model (blue) with . For both cases, the accretion is episodic. During the brief burst phases, the rates are larger than the Eddington accretion rate, whereas the long-term average rate is for both cases. The peak values for the trapping model are several times larger because the maximum luminosity in the trapping model is limited to and thus radiation pressure is less efficient than in the standard model. However, the average accretion rate is still even in the trapping model. This fact demonstrates that rapid accretion onto the BH is strongly quenched primarily not by radiation pressure but by radiative heating.

The physical reason why the accretion history becomes burst-like can be understood as follows. Fig. 4 shows radial profiles of the (a) number density and (b) temperature during one cycle in the trapping model. The corresponding epochs are indicated by open circles in Fig. 3. Once the gas accretes onto the BH, the emergent radiation starts to ionize and heat the gas around the BH. The hot H region expands until photoionization and recombination are balanced ( cm). The H region becomes larger than the original Bondi radius, cm (Fig. 4b). In the burst phase (curves 1–2 in Fig. 4), the BH’s gravity accelerates the inflow within a new sonic point in the H region ( cm for ), while the gas pressure (and partially the radiation pressure force) pushes the gas outward at cm and the accretion rate is reduced. As a result, in the quiescent phase (curves 3–5 in Fig. 4), a density cavity forms within the H region, where the outward and inward gas pressure forces are balanced. However, as the depletion of the hot gas proceeds and the outward gas pressure force decreases, a density bump forms inside the cavity at cm (see the sharp feature in the curve 5 in Fig. 4a). Then, the gas is accelerated by the inward gas pressure (i.e. ) due to this bump. Finally, the density bump falls into the central BH and produces burst-like accretion again (curve 6 in Fig. 4a).

A similar episodic behavior of the accretion rate for a constant radiative efficiency has been already described in previous works (e.g. Ciotti & Ostriker, 2001; Milosavljević et al., 2009; Milosavljević, Couch & Bromm, 2009; Park & Ricotti, 2011, 2012). Although the radiation pressure force just after the burst-like accretion is stronger for the constant efficiency, the episodic behavior is due to radiative heating and leads to a similar behavior in the trapping model.

Figure 4: Radial profiles of the (a) number density and (b) temperature of the accretion flow over one cycle for and (in the outer region). The curves show the profiles at (1) yr (red solid), (2) yr (green long-dashed), (3) yr (blue short-dashed), (4) yr (magenta dotted), (5) yr (light-blue long-dashed-dotted) and (6) yr (orange short-dashed-dotted) The corresponding epochs are indicated by open circles in Fig. 3. The horizontal lines are initial conditions ( and ) and the initial Bondi radius is cm. In this case where the size of the H region () is larger than , the gas dynamics is affected by radiation feedback (mainly heating), resulting in episodic accretion.

3.2 Higher BH mass cases

Next, we move onto the cases of more massive BHs. Fig. 5 presents the accretion history for different BH masses (). We adopt the radiative efficiency of the trapping model (Eq. 18). For the lowest BH mass, the accretion occurs episodically as in the case of ordinary stellar-remnant massive BH as in §3.1. With the BH mass increasing, the accretion rate becomes less time-dependent and the average accretion rate is , which corresponds to the case of a radiative efficiency many papers assume motivated by a thin disk (e.g. Shakura & Sunyaev, 1973). For the larger BH mass, we find a big jump of the accretion rate from to at yr for . After the transition for each case, the accretion rate approaches a constant value, which is identical to the Bondi accretion rate.

Figure 5: Accretion rate history for a higher-mass BH (in the outer region). The curves show the cases for (green long dashed), (blue short dashed), (magenta dotted), (red solid), and (orange dot-dashed). The radiative efficiency is assumed to be (the trapping model). The density of the ambient gas is . Open circles mark the six epochs at which we show the radial profiles in Fig. 6. For lower BH mass (), the average accretion rate is limited to , which is similar to that shown in Fig. 3. For higher BH mass (), a big jump of the accretion rate occurs and the accretion rate approaches a constant value, where ().

In Fig. 6, we present radial profiles of the density, temperature, and inflow rate at different times for the case with the jump of the accretion rate ( and ). In the early phase, the accretion history has several quiescent and burst phases as in the lower BH mass case (see §3.1). In this case, however, the dense shell which developed at the edge of the H region pushes the hot ionized gas inward at yr. This results in the big jump of the accretion rate (phase 3–5 in Figs. 5 and 6). Note that the inward acceleration by the positive pressure gradient force is subdominant compared to that by the BH’s gravity, in contrast to the burst-like accretion with a low accretion rate (see §3.1). Fig. 6(b) shows that the size of the H region shrinks and it disappears. Finally, the temperature profile is nearly isothermal with . The final profile of the inflow rate is almost constant (Fig. 6c), which is a steady and approximately isothermal Bondi solution with .

Figure 6: Time evolution of radial profiles of the (a) number density, (b) temperature and (c) inflow rate () for and (in the outer region). The curves show the profiles at (1) yr (red solid), (2) yr (green long-dashed), (3) yr (blue short-dashed), (4) yr (magenta dotted), (5) yr (orange dashed-dotted), and (6) yr (black solid). The initial Bondi radius is cm. Note that in panel (c) for clarify, we do not show positive values (phase 2; green) at cm. In this case where the H region is confined within the Bondi radius (), the big jump of the accretion rate occurs. Finally, the accretion flow becomes a steady and isothermal Bondi solution with .

We note the final fates with , where the accretion rates are almost constant at the end of the simulations. For these marginal cases between the burst-like accretion () and the hyper-Eddington accretion (), a dense shell develops just outside the H region but radiative heating still suppresses the strong inflow. We do not find a transition to a steady hyper-Eddington phase at least within the simulation time, which is times longer than the dynamical time at . However, the dense shell might fall into the BH at some point in a longer-term simulation if the gravity slightly exceeds the outward force of gas pressure. If in the case of density inversion flows, there is a non-linear density fluctuation due to the Rayleigh-Taylor instability (see §5.1), then accretion might increase even for the marginal cases of .

Figure 7: Summary of the results for different values of the BH mass and number density of the ambient gas . Each symbol indicates whether the final result is hyper-Eddington accretion (circle), or constant or episodic accretion at the rate of (square). The dashed line marks the boundary between the two accretion modes: . The region above the boundary indicates the conditions required to realize steady hyper-Eddington accretion.

3.3 Parameter dependence

To discuss the conditions required to realize a steady hyper-Eddington accretion flow, we conduct four additional simulations for different BH mass () and number density of the ambient gas ( and ).

The results for are similar to the cases shown in §3.2. For a larger BH mass or higher density, the accretion rate jumps to the corresponding Bondi rate. On the other hand, for smaller BH mass or lower density, the accretion rate is limited to . Fig. 7 summarizes our results. Each symbol indicates whether the final result is steady hyper-Eddington accretion (circle) or sub-Eddington accretion (; square). The dashed line marks the boundary between the two accretion modes: .

3.4 Analytic argument

We here give a simple analytic argument for the conditions required for hyper-Eddington accretion. As we explained in §3.2, the relation between the size of the H region and the Bondi radius is important to determine whether the transition to the hyper-Eddington accretion phase occurs. For the lower BH mass case, the ionizing front propagates outside the Bondi radius and never shrinks (; see Fig. 4 b). The radiation heating and pressure in this case can affect the gas dynamics at the Bondi radius. Thus, the accretion is suppressed and (blue curve in Fig. 5). For the higher BH mass, the H region is always confined within the Bondi radius (; Fig. 6 b) and shrinks dramatically at yr. As a result, the accretion from the Bondi radius cannot halt due to radiation feedback, and the accretion flow becomes steady with (see Fig. 5).

The size of an H region in an uniform-density medium with is estimated as


where is the ionizing photon number flux (in units s) and is the H radiative recombination coefficient (case B). For the trapping model (), the maximum value of is , where is the average energy of ionizing photons. We obtain


where is the temperature inside the H region and we set eV. This value is larger by a factor of than the actual value because Eq. (26) neglects the fact that the density profile has a steep slope (; ) within . We set . Thus, the condition for the transition to the hyper-Eddington accretion () is written as




For , the critical accretion ratio is , which agrees well with the numerical simulation (see Figs. 5 and 6). In the following discussions (§4 and §5), we set the critical rate required for hyper-Eddington accretion to , which is given by the dashed line in Fig. 7.

3.5 Inner-region simulations

For the case of the hyper-Eddington accretion flow, we further conduct simulations of the inner region, resolving the trapping radius. Here, we study whether gas and radiation in the inner region affect or modify the hyper-Eddington solution of the outer region. As we mentioned in §2.3.2, we have run several simulations for . Since the choice of is still arbitrary, we attempt to determine the physically correct value of so that our solution smoothly connects with a small accretion disk well inside .

First, we show the results of the inner-region simulation for the hyper-Eddington solution shown in Fig. 6 ( and ) with . Fig. 8 presents radial profiles of the (a) temperature and (b) inflow velocity at the end of the simulations ( s) after a steady state is reached. For , the temperature profile is almost constant () outside the photosphere. The size of the photosphere is cm and the effective temperature is (). Inside the photosphere, the gas temperature reaches a few at the inner boundary ( cm). The inflow velocity is close to free-fall () outside the photosphere, where , while the radiation pressure force becomes important inside the photosphere, the inflow velocity decreases and the density is slightly piled-up within cm.

Figure 8: The profile of the (a) gas temperature and (b) inflow velocity at the end of the simulation for and (in the inner region). The curves correspond to the cases with the photospheric luminosity of (red solid), (green long-dashed), (blue short-dashed), (magenta dotted), (light-blue long-dashed-dotted), and (orange short-dashed-dotted). Open circles indicate the location of the photosphere for each case. For , the comoving luminosity exceeds the Eddington luminosity within (see Fig. 9b), and the velocity begins to decelerate due to radiation pressure and would approach a settling solution (). A solution with could connect with a hyper-Eddington accretion disk solution smoothly at small radii.
Figure 9: The profile of the radiative luminosity at the end of the simulation for , (in the inner region). In panel (a), where , red curve shows the luminosity in the rest (observed) frame; for (solid) and (dotted). Green short-dashed and blue long-dashed curve present the luminosity in the fluid’s comoving frame and the advection luminosity, respectively. The open circle indicates the position of the photosphere . The horizontal line is the Eddington luminosity. The rest-frame luminosity (red curve) decreases inward from because of photon trapping, and becomes negative within the trapping radius . In other words, most of the energy released at is trapped within the inflow and falls into the BH, while the leaked energy can reach outside the photosphere () by radiative diffusion (green curve). Panel (b) shows the radial profile of the comoving luminosity with different values of .

Fig. 9 (a) shows the profile of the radiative luminosity at the final state ( s) for the case with . Each curve corresponds to the luminosity in the rest (observed) frame (red), in the fluid (comoving) frame (green), and the advection luminosity (blue). Note that is given by the diffusion approximation, and as shown in Eq. (