Rapid growth of black holes accompanied with hot or warm outflows exposed to anisotropic super-Eddington radiation

Rapid growth of black holes accompanied with hot or warm outflows exposed to anisotropic... Abstract We perform two-dimensional radiation hydrodynamical simulations of accretion flows on to a black hole (BH) with a mass of 103 ≤ MBH/ M⊙ ≲ 106 in order to study rapid growth of BHs in the early Universe. For spherically symmetric flows, hyper-Eddington accretion from outside the Bondi radius can occur unimpeded by radiation feedback when MBH ≳ 104 M⊙(n∞/105 cm − 3) − 1(T∞/104 K)3/2, where the density and temperature of ambient gas are initially set to n∞ = 105 cm−3 and T∞ = 104 K. Here, we study accretion flows exposed to anisotropic radiation from a nuclear accretion disc with a luminosity higher than the Eddington value (LEdd) due to collimation towards the bipolar directions. We find that, unlike the spherically symmetric case, even less massive BHs with MBH < 104 M⊙ can be fed at high accretion rates of ≳ LEdd/c2 through the equatorial region, while ionized regions expand towards the poles producing hot outflows with T ∼ 105 K. For more massive BHs with MBH ≳ 5 × 105 M⊙, intense inflows of neutral gas through the equator totally cover the central radiating region due to the non-radial gas motions. Because of efficient recombination by hydrogen, the entire flow settles in neutral and warm gas with T ≃ 8000 K. The BH is fed at a rate of ∼5 × 104LEdd/c2 (a half of the inflow rate from the Bondi radius). Moreover, radiation momentum absorbed by neutral hydrogen produces warm outflows towards the bipolar directions at ∼ 10 per cent of the BH feeding rate and with a velocity several times higher than the escaping value. black hole physics, quasars: supermassive black holes, cosmology: theory 1 INTRODUCTION Observations of bright quasars at high redshifts (z ≳ 6) have revealed that supermassive black holes (SMBHs) with MBH ≳ 109 M⊙ are formed within ≲ 1 Gyr after the beginning of the Universe (e.g. Fan et al. 2004; Mortlock et al. 2011; Wu et al. 2015). Although SMBHs play crucial roles in the co-evolution with their host galaxies through feedback processes (e.g. Silk & Rees 1998; King 2003; Murray, Quataert & Thompson 2005; Kormendy & Ho 2013), their formation and growth processes are still open questions. One possible formation scenario of high-z SMBHs is that stellar mass black holes (BHs) formed by collapse of first generation stars (Population III, hereafter PopIII) with ∼100 M⊙ grow via gas accretion (e.g. Madau & Rees 2001; Haiman & Loeb 2001; Volonteri, Haardt & Madau 2003). In the standard accretion physics, BH growth would be suppressed by radiation emitted from a nuclear accretion disc because the radiation luminosity increases with the accretion rate on to the BH. If the radiative efficiency of the disc is η ≃ 0.1 (Shakura & Sunyaev 1973; Soltan 1982; Yu & Tremaine 2002), the radiation luminosity exceeds the Eddington luminosity (LEdd) at the accretion rate of ≳ LEdd/(ηc2), above which the radiation force dominates the gravity of the central BH. Assuming this Eddington-limited accretion rate, it takes ≳ 1 Gyr for stellar-mass seed BHs to grow up to ∼109 M⊙. This fact requires more rapid accretion mechanisms to form SMBHs by z ∼ 6, when the corresponding cosmic age is ∼0.8 Gyr. As alternative possibilities, formation of massive seed BHs with ≳ 103-105 M⊙ has been proposed: direct-collapse of supermassive stars in protogalaxies (e.g. Inayoshi, Omukai & Tasker 2014; Regan, Johansson & Haehnelt 2014; Visbal, Haiman & Bryan 2014; Becerra et al. 2015; Inayoshi & Tanaka 2015; Chon et al. 2016; Latif, Schleicher & Hartwig 2016) or runaway collisions of stars in a dense cluster (e.g. Devecchi & Volonteri 2009; Katz, Sijacki & Haehnelt 2015; Yajima & Khochfar 2016; Stone, Küpper & Ostriker 2017; Sakurai et al. 2017). However, since the BH growth time even from such massive seeds would be shorten only by a factor of 2 in the case of the Eddington limited accretion, a high duty cycle ( ∼ 100 per cent) is still required to explain the existence of SMBHs with ∼109 M⊙ at high redshifts. The possibility of rapid growth of seed BHs has been discussed (e.g. Volonteri & Rees 2005; Tanaka & Haiman 2009; Alexander & Natarajan 2014; Madau, Haardt & Dotti 2014; Pacucci & Ferrara 2015; Pezzulli, Valiante & Schneider 2016; Valiante et al. 2016). By means of radiation hydrodynamical simulations of accretion flows on to a BH, Ohsuga et al. (2005) concluded that super-Eddington accretion can be realized as long as sufficient gas is supplied at the vicinity of the central BH (see also Ohsuga et al. 2009; Ohsuga & Mineshige 2011; Jiang, Stone & Davis 2014; Sa̧dowski et al. 2015; Takahashi et al. 2016). This is because diffusive photons are efficiently trapped within optically thick accretion flows and the radiation force on to the gas is alleviated (e.g. Begelman 1979). However, it is unclear how sufficient amount of gas is supplied into the nuclear region. In fact, radiation heating and ionization effectively suppress gas accretion from larger scales (Ciotti & Ostriker 2001; Alvarez, Wise & Abel 2009; Ciotti, Ostriker & Proga 2009; Milosavljević et al. 2009; Park & Ricotti 2011, 2012). Inayoshi, Haiman & Ostriker (2016) (hereafter IHO16) have found a self-consistent solution of steady hyper-Eddington accretion in spherically symmetric systems with radiation hydrodynamical simulations. This accretion flow has a high accretion rate of ≳ 5000LEdd/c2, unimpeded by negative feedback due to radiation force and heating simultaneously. The solution is composed of a radiation-dominated core, where photon trapping due to electron scattering works efficiently, and an accreting envelope that follows a Bondi profile with 8000 K. When the emergent radiation luminosity is limited to ≲ LEdd because of photon trapping, radiation from the central region does not affect the gas dynamics at larger scales. The hyper-Eddington accretion is then realized when a BH is embedded in a dense gas cloud, for which the following condition is satisfied,   \begin{equation} M_{\rm BH} \gtrsim 10^4 \,{\rm M}_{\odot }\left(\frac{n_\infty }{10^5 \,{\rm cm}^{-3}}\right)^{-1}\left(\frac{T_\infty }{10^4\, {\rm K}}\right)^{3/2}, \end{equation} (1)where n∞ and T∞ are density and temperature of the ambient gas, respectively. Hereafter, we define MBH, x = (MBH/10x M⊙), n∞, 5 = (n∞/105 cm−3), and T∞, 4 = (T∞/104 K). Note that equation (1) is identical to the condition where the Bondi radius is larger than the size of the ionized region. Even if a high luminosity with >LEdd emerges from the centre, the above conditions are not affected as long as the luminosity is as low as (10 − 30) × LEdd (Sakurai, Inayoshi & Haiman 2016). To satisfy the condition of equation (1), gas-rich regions with n ≳ 107 cm − 3 are required for stellar-mass PopIII BHs with masses of ≲ 100 M⊙. During the assembly of a protogalaxy, some PopIII BHs would fall into a gaseous circumnuclear disc, where the gas accretion rate is likely to be ≳ 5000LEdd/c2 (Ryu et al. 2016, see also Lupi et al. 2016). Radiation flux from an accretion disc around the nuclear BH is more likely to be anisotropic, i.e. preferentially collimated around the rotation axis of the disc (Ohsuga et al. 2009; Ohsuga & Mineshige 2011; Jiang et al. 2014; Sa̧dowski et al. 2015; Takahashi et al. 2016). On the other hand, the disc-like accretion could reduce the feedback effect because intense radiation is significantly attenuated by the dense accreting material in the disc (Ohsuga & Mineshige 2007; Sugimura et al. 2016). In this paper, we investigate how intermediate massive BHs with 103 ≤ MBH ≤ 5 × 105 M⊙, embedded in an initially uniform gas with n∞ = 105 cm−3 and T∞ = 104 K, can grow via accretion when the inflow gas is exposed to the super-Eddington, anisotropic radiation. To address this issue, we perform two-dimensional radiation hydrodynamical simulations, including multifrequency radiation transfer and non-equilibrium primordial chemistry. We confirm that hyper-Eddington accretion from the Bondi radius is likely to occur when radiation is more anisotropic and the BH mass is higher. Furthermore, we demonstrate that for MBH ≳ 5 × 105 M⊙, neutral gas inflows through the equator totally cover the central radiating region due to the non-radial gas motions, and thus the emergent radiation in all directions is blocked. Because of efficient recombination by hydrogen, the entire flow results in neutral and warm gas. In this case, the BH feeding rate is as high as the Bondi accretion rate, and also warm outflows are produced by absorption of radiation momentum. The rest of this paper is organized as follows. In Section 2, we describe the methodology for radiation hydrodynamical simulations. In Section 3, we show our simulation results and discuss the effects of anisotropic radiation. Finally, we summarize the main conclusions of this paper and discuss caveats of our simulations in Section 4. 2 SIMULATION METHOD We perform two-dimensional hydrodynamical simulations including one-dimensional radiation transfer and chemical reaction networks. Our purpose is to study necessary conditions for hyper-Eddington accretion led by sufficient supply of gas from large scales. Gas accretion begins from the Bondi radius, within which the gravity of the central BH dominates the gas-pressure gradient force, given by Bondi (1952) as   \begin{equation} R_{\rm B}\equiv \frac{GM_{\rm BH}}{c_{\infty }^2} \simeq 1.97\times 10^{18} M_{\rm BH,4} T_{\infty ,4}^{-1} \ {\rm cm}, \end{equation} (2)where $$c_{\infty }=\sqrt{\gamma \mathcal {R} T_{\infty } /\bar{\mu }}$$ is the sound speed, γ is the specific heat ratio, $$\mathcal {R}$$ is the gas constant, and $$\bar{\mu }$$ is the mean molecular weight. As a reference of the accretion rate, we define the Bondi accretion rate for γ = 1 as   \begin{equation} \dot{M}_{\rm B} \equiv {\pi} e^{3/2} \rho _\infty \frac{G^2 M_{\rm BH}^2}{c_{\infty }^3}, \end{equation} (3)and the Eddington accretion rate as $$\dot{M}_{\rm Edd} \equiv L_{\rm Edd}/c^2$$, where LEdd = 4πcGMBH/κes and κes is the electron scattering opacity. Note that the Bondi radius and rate are calculated for γ = 1, $$\bar{\mu }=1.23$$, and T∞ = 104 K as a reference value through this paper, and that these are self-consistently calculated in our simulations. 2.1 Basic equations We consider a gas sphere exposed to intense radiation from the accretion flow on to the central BH, solving two-dimensional hydrodynamical equations assuming axisymmetric flows. We here use the hydrodynamical simulation code taken from Takahashi & Ohsuga (2013). The mass flux for the ideal fluid is computed using the Harten–Lax–vanLeer Riemann solver (Harten, Lax & van Leer 1983), and the second-order accuracy in space and time is ensured (van Leer 1977). We adopt the spherical coordinates of (r, θ, ϕ), defining the polar axis (θ = 0 and π) as directions perpendicular to the disc plane. The basic equations of hydrodynamics are the following: the equation of continuity   \begin{equation} \frac{\partial \rho }{\partial t} + \nabla \cdot (\rho {\boldsymbol v}) = 0, \end{equation} (4)and the equations of motion,   \begin{eqnarray} \frac{\partial (\rho v_r)}{\partial t} + \nabla \cdot (\rho v_r {\boldsymbol v}) = -\frac{\partial p}{\partial r} + \rho \left( \frac{v_{\theta }^2}{r} + \frac{v_{\phi }^2}{r} \right) - \rho \frac{\partial \psi}{\partial r}+ f_{\rm rad}, \nonumber\\ \end{eqnarray} (5)  \begin{equation} \frac{\partial \left(\rho rv_{\theta }\right)}{\partial t} + \nabla \cdot (\rho r v_{\theta } {\boldsymbol v}) = -\frac{\partial p}{\partial \theta } + \rho v_{\phi }^2\cot {\theta }, \end{equation} (6)  \begin{equation} \frac{\partial \left(\rho rv_{\phi }\sin {\theta }\right)}{\partial t} + \nabla \cdot \left(\rho r v_{\phi }\sin {\theta } {\boldsymbol v}\right) = 0, \end{equation} (7)where ρ is the gas density, $${\boldsymbol v}=(v_r,v_{\theta }, v_{\phi })$$ is the velocity, p is the gas pressure, and frad is the radiation force. We consider the gravity of the central BH (r = 0) and neglect the gas self-gravity. Since the general relativistic effect is negligible, the gravitational potential is given by ψ = −GMBH/r. We also solve the energy equation   \begin{equation} \frac{\partial e}{\partial t} + \nabla \cdot [(e+p){\boldsymbol v}] = -\frac{GM_{\rm BH}\rho }{r^2}v_r -\Lambda + \Gamma , \end{equation} (8)where $$e=e_{\rm int}+\rho |{\boldsymbol v}|^2/2$$ is the gas total energy density, eint is the gas internal energy density, Λ is the net cooling rate per volume, and Γ is the heating rate due to radiation from the central region. We assume the equation of state of ideal gas as p = (γ − 1)eint for γ = 5/3. We consider radiative cooling by bound–bound transitions of H, He, He+ atoms and free–free emission (Glover & Jappsen 2007). To estimate their rates, we solve chemical reaction networks including six species of H, H+, He, He+, He++, and e−. The abundance of He nuclei relative to H nuclei is set to 8.33 × 10−2. Here, photoionization, collisional ionization, and radiative recombination are considered (Abel et al. 1997; Glover & Jappsen 2007). Instead of treating photoionization by diffusive recombination photons, we adopt the on-the-spot approximation, where the case A recombination rate is replaced by the case B rate. To update the chemical abundances stably, we adopt a semi-implicit method (Anninos et al. 1997), setting the time-steps shorter than chemical time-scales for all the cells defined by $$t_{\rm chem} \equiv (x_{\rm e}+0.001x_{\rm H})/\dot{x}_{\rm e}$$, where xe and xH are the abundance of electrons and neutral hydrogens, respectively (Whalen & Norman 2006, 2008). Note that the cooling/heating part of the energy equation (equation 8) is updated with an implicit method in order to solve it stably and save computation time. We solve the multifrequency, steady radiative transfer equation because the light crossing time is much shorter than the hydrodynamical time-scale,   \begin{equation} \frac{1}{r^2}\frac{{\rm d}}{{\rm d}r}(r^2F_{\nu }) = -\rho \kappa _{\nu }cE_{\nu }, \end{equation} (9)where Fν is the radiation flux, Eν is the radiation energy density, and κν is the absorption opacity. The frequency range is set to hνmin(=13.6eV) ≤ hν ≤ hνmax(=10keV), where h is the Planck constant. Note that the radial component of the radiation flux is calculated (see Section 4 for more details). Since the ionized gas is optically thin to electron scattering and bound-free transitions, we assume Fν ≈ cEν on the right-hand side of equation (9). Most ionizing photons are absorbed by neutral hydrogen at the edge of the ionized region. Non-radial components of diffusive photons due to recombination are neglected in our simulation (see also Section 4). The ionization rate coefficients kph of H, He, He+, and photoionization heating rate Γ(=ΓH) are calculated following the photon-conserving manner (Whalen & Norman 2006). Here, ΓH is the heating rate due to H ionization. We consider radiation force due to electron scattering and bound-free absorption of H atoms as   \begin{equation} f_{\rm rad} = \frac{nx_{\rm e}}{c}\int _{\nu _{\rm min}}^{\nu _{\rm max}}\sigma _{\rm es}F_{\nu }{\rm d}\nu + \frac{\Gamma _{\rm H}}{c}. \end{equation} (10) In order to integrate equation (9), we simply set a radiation source with a single power-law spectrum at the centre, Lν = L0(ν/νmin)−α at νmin ≤ ν ≤ νmax and Lν = 0 at ν ≤ νmin and νmax ≤ ν, where α = 1.5 is adopted. The normalization factor L0 is set by the total luminosity ($$L=\int _{\nu _{\rm min}}^{\nu _{\rm max}} L_{\nu } {\rm d}\nu$$). We set a model for radiation luminosity emitted by an accretion disc, which is unresolved, at the inner boundary,   \begin{equation} \frac{\displaystyle L}{\displaystyle L_{\rm Edd}}= \begin{array}{ll}2\left[1+\ln {(\dot{m}/20)}\right] & {\rm for}\ \dot{m}\ge 20, \\ \dot{m}/10 & {\rm for}\ \dot{m}<20, \end{array} \end{equation} (11)(Watarai et al. 2000), where $$\dot{m} \equiv \dot{M}/\dot{M}_{\rm Edd}$$. Furthermore, we assume anisotropic radiation fields as   \begin{equation} F_{\nu }(r=r_{\rm min}, \theta ) = \frac{\displaystyle (N+1)L_{\nu }}{\displaystyle 4{\pi} r_{\rm min}^2}\cos ^N{\theta }, \end{equation} (12)where the radiation flux Fν(rmin) is normalized so that $$L_{\nu }(r_{\rm min})=\int F_{\nu }(r_{\rm min}) r_{\rm min}^2 {\rm d}\Omega$$ and rmin is the radius of the inner most grid in the computation domain (see Section 2.2). Unlike Sugimura et al. (2016), shadow regions are not assumed in our simulations (see Section 4 for more details). Here, we explore the effects of radiation anisotropy; the isotropic case N = 0 (Model A), the anisotropic cases N = 2 (Model B), and N = 4 (Model C) for MBH = 103 M⊙. We also study cases for higher BH masses with MBH = 5 × 104 M⊙ (Model D) and 5 × 105 M⊙ (Model E) for N = 4. Our models are summarized in Table 1. Table 1. Model  MBH( M⊙)  N  $$\dot{M}/\dot{M}_{\rm Edd}$$  L/LEdd  A  103  0  ∼0.9  ∼0.09  B  103  2  ≈2.7  ≈0.27  C  103  4  ≈44  ≈3.5  D  5 × 104  4  ∼2 × 103  ∼11  E  5 × 105  4  ∼5 × 104  ∼18  Model  MBH( M⊙)  N  $$\dot{M}/\dot{M}_{\rm Edd}$$  L/LEdd  A  103  0  ∼0.9  ∼0.09  B  103  2  ≈2.7  ≈0.27  C  103  4  ≈44  ≈3.5  D  5 × 104  4  ∼2 × 103  ∼11  E  5 × 105  4  ∼5 × 104  ∼18  Column: (1) model ID, (2) BH mass, (3) radiation anisotropy, namely Frad(θ) ∝ cos Nθ, (4) accretion rate, and (5) corresponding luminosity. In columns 4 and 5, we show the time-averaged values in Model A, the values at the end of simulations in Models B–D, and the time-averaged values over the last one cycle of episodic accretion in Model E. View Large 2.2 Initial and boundary conditions We set a computational domain of rmin ≤ r ≤ rmax and 0 ≤ θ ≤ π. Here, we adopt rmin = 0.04RB, rmax = 30RB for Models A, B, and C, and rmin = 0.007RB, rmax = 6RB for Models D and E. In Models A–C, the computational domain is located relatively outward because the size of the ionized region tends to larger than the Bondi radius. We set logarithmically spaced grids in the radial direction to simulate the flow within the Bondi radius in high resolution. We set uniformly spaced grids in the polar direction. The number of the grid points is set to (Nr, Nθ, Nϕ) = (100, 120, 1). Note that our simulations do not consider accretion flows within rmin. Instead, we assume properties of radiation emitted from the central region (see Section 2.1) and discuss gas accretion outside from the Bondi radii. As our initial conditions, we set a neutral uniform and static ($${\boldsymbol v}= 0$$) gas cloud with the density n∞ = 105 cm−3 and temperature T∞ = 104 K. The BH masses are assumed to be constant throughout our simulations. We impose the absorption inner boundary conditions for the gas density, gas pressure, and velocity to be damped smoothly (e.g. Kato, Mineshige & Shibata 2004), and the free outer boundary conditions for three components of the velocity and the specific entropy. We also fix the same gas density at r = rmax as the initial value for grids with an inflow velocity i.e. vr(r = rmax) < 0, otherwise the free boundary conditions are imposed for the density. The reflection symmetry with respect to the polar axis is imposed for non-radial components of the velocity. 3 RESULTS 3.1 Effects of anisotropic radiation Fig. 1 presents the time evolution of gas accretion rates on to a BH with MBH = 103 M⊙ for Model A (N = 0), B (N = 2), and C (N = 4). For the isotropic radiation, the accretion occurs episodically and the time-averaged rate ($${\sim } 0.9 \dot{M}_{\rm Edd}$$) is below the Eddington rate. For the anisotropic radiation, the accretion is less episodic and their rates tend to increase with time because of continuous accretion of neutral gas through the equatorial plane as shown later. The time-averaged accretion rate becomes higher with N (i.e. more anisotropic radiation). In Model C, the accretion rate is as high as $${\sim } 44 \dot{M}_{\rm Edd}$$ at the end of the simulation, where the radiation luminosity is higher than the Eddington value, namely ∼3.5LEdd. We continue the simulation until the ionization front reaches the outer boundary at t ≃ 1.2 × 105 yr, which is much longer than $$t_{\rm dyn} \equiv {\pi} \sqrt{R_{\rm B}^3/(8GM_{\rm BH})} \simeq 8.4 \times 10^3 {M_{\rm BH,3}} {\rm yr}$$. The accretion rate and the radiative luminosity at the end of each simulation are summarized in Table 1. Figure 1. View largeDownload slide Time evolution of gas accretion rates on to a 103 M⊙ BH in cases with isotropic radiation (Model A) and anisotropic radiation, Frad ∝ cos 2θ (Model B) and Frad ∝ cos 4θ (Model C). In Model A, the accretion occurs episodically and the time-averaged rate is below $$\dot{M}_{\rm Edd}$$. The accretion rate in Model C exceeds the Eddington rate ($$\dot{M}>10\dot{M}_{\rm Edd}$$). Figure 1. View largeDownload slide Time evolution of gas accretion rates on to a 103 M⊙ BH in cases with isotropic radiation (Model A) and anisotropic radiation, Frad ∝ cos 2θ (Model B) and Frad ∝ cos 4θ (Model C). In Model A, the accretion occurs episodically and the time-averaged rate is below $$\dot{M}_{\rm Edd}$$. The accretion rate in Model C exceeds the Eddington rate ($$\dot{M}>10\dot{M}_{\rm Edd}$$). In Fig. 2, we show two-dimensional distribution of the gas density and temperature at t = 1.1 × 105 yr in Models A and C, respectively. In the isotropic case (Model A), radiation from the central region heats up and ionizes the surrounding gas equally in all directions. When the ionized region expands and the temperature outside the Bondi radius increases, gas supply from large radii is suppressed since $$\dot{M}_{\rm B} \propto T_\infty ^{-3/2}$$. In this quiescent phase, thermal gas pressure is initially balanced inside and outside the ionized region. However, since the ionized gas within a new sonic point can accrete on to the centre, the gas pressure in the ionized region is gradually reduced due to gas depletion. As a result, the gas at the ionization front is accelerated inward by the positive pressure gradient (∂p/∂r > 0), leading to episodic accretion (see Fig. 1) as shown in previous studies (e.g. Ciotti & Ostriker 2001; Park & Ricotti 2011, 2012, IHO16). Note that radiation heating is a dominant process for suppressing the gas accretion, since the radiation luminosity is below the Eddington value (i.e. L < LEdd). Figure 2. View largeDownload slide Two-dimensional distribution of gas density (left-hand panel) and temperature (right-hand panel) for Model A (upper) and Model C (bottom), respectively, at t = 1.06 × 105 yr (≈13tdyn). In Model C (anisotropic radiation), the ionized region does not exist around the equatorial plane, from which a higher accretion rate is allowed. Figure 2. View largeDownload slide Two-dimensional distribution of gas density (left-hand panel) and temperature (right-hand panel) for Model A (upper) and Model C (bottom), respectively, at t = 1.06 × 105 yr (≈13tdyn). In Model C (anisotropic radiation), the ionized region does not exist around the equatorial plane, from which a higher accretion rate is allowed. In the anisotropic case (Model C), the shape of the ionized region is no longer isotropic. Although the ionization front expands towards the polar directions, the radiation flux in the equatorial direction (θ = π/2) is not intense enough to form an ionized region. Since the radiation luminosity viewed from the polar directions is significantly higher than the Eddington luminosity, strong outflows are launched by the radiation force on to electrons. The outflows collide with the ambient gas at r ≳ RB and form strong shocks with a high temperature of ≳ 5 × 104 K. The gas near the equator remains neutral and the temperature is kept at T ≃ 8000 K due to efficient atomic hydrogen cooling. Therefore, the gas can accrete on to the central BH through the neutral region increasing the density, unimpeded by radiation feedback. The radial profiles at θ = π/2 of the density and inflow velocity approach those for a Bondi solution with T ∼ 8000 K (ρ ∝ r−3/2 and vr ∝ r−1/2 at r ≲ RB). The accretion rate through the equatorial plane is estimated as $$\dot{M}_{\rm HI} \simeq \dot{M}_{\rm B}(\Omega _{\rm HI} / 4{\pi} )$$, where ΩHI is the solid angle of the neutral region. Since a half angle of the neutral region is ΘHI ≡ π/2 − θ ≈ 5° and $$\dot{M}_{\rm B} \approx 670 \dot{M}_{\rm Edd}$$, the inflow rate of the neutral gas is estimated as $$\dot{M}_{\rm HI}\simeq 58 \dot{M}_{\rm Edd}$$. This value agrees to the simulation result of $${\simeq } 44 \dot{M}_{\rm Edd}$$. Unlike our set-ups, Sugimura et al. (2016) assumed a shadow region with a half angle of ΘShad = 10°–45°, within which the radiation flux is completely shut off and accretion of neutral gas is allowed. In this case, they also found that gas accretion rates are given by $$\dot{M} \simeq \dot{M}_{\rm B}(\Omega _{\rm Shad} / 4{\pi} )$$, which is consistent with our results. 3.2 Transition to wholly neutral accretion Next, we examine a case with the highest BH mass of 5 × 105 M⊙ (Model E) but with the same radiation anisotropy (N = 4) as in Model C. Fig. 3 shows two-dimensional distribution of the gas density (left-hand panels) and temperature (right-hand panels) in the early stage at t ≃ 0.14tdyn (top panels) and the late stage at t ≃ 0.49tdyn (bottom panels), respectively. The time evolution of the accretion flow in Model E consists of the following three phases. Initially, the ionized region expands to r ∼ RB preferentially towards the polar directions, and gas accretion is allowed through the neutral region (top panels in Fig. 3). Unlike in Model C with lower MBH (i.e. lower $$\dot{m}$$), the opening angle of the ionized region in Model E gradually shrinks with time, allowing penetration of neutral gas towards the centre through the equator. The dense inflow totally covers the central region, where ionizing photons from the sink cell are completely absorbed by neutral hydrogen. As a result, the radiation flux in all directions (even in the polar directions) is blocked, and the ionized region is confined within r ∼ 0.01RB due to radiative recombination (bottom panels in Fig. 3). The radiation momentum absorbed by neutral gas produces warm outflows towards the bipolar directions. Figure 3. View largeDownload slide Two-dimensional distribution of gas density (left-hand panel) and temperature (right-hand panel) in Model E (MBH = 5 × 105 M⊙ and N = 4) in the early stage at t = 5.9 × 105 yr (≃0.14tdyn) (top), and in the late stage (after the transition) at t = 2.1 × 106 yr (≃0.49tdyn) (bottom), respectively. Figure 3. View largeDownload slide Two-dimensional distribution of gas density (left-hand panel) and temperature (right-hand panel) in Model E (MBH = 5 × 105 M⊙ and N = 4) in the early stage at t = 5.9 × 105 yr (≃0.14tdyn) (top), and in the late stage (after the transition) at t = 2.1 × 106 yr (≃0.49tdyn) (bottom), respectively. In order to understand how and why this transition of the accretion episode occurs, we examine structure of the gas flow in the central region (rmin ≤ r ≲ 0.04RB). Fig. 4 shows the density distribution and the velocity vectors of the flows. Each panel shows the time sequence from the left- (before the transition) to the right-hand side (after the transition). Fig. 5 also presents profiles of the electron fraction, the ram pressure ($$\rho v_{\theta }^2$$) and the gas pressure (pgas) as functions of the polar angle at r = 0.009RB. Before the transition shown in Fig. 4 (panel i1), geometrically thin and dense shock layers form along directions of θ ≈ 55° and 70°, located near the boundary between ionized and neutral gas (see Fig. 5). The inflow from the neutral region crosses the shock with a supersonic velocity, namely |vθ| > cs, and is decelerated (cs is the local sound speed). The shocked gas is ionized and accelerated outward in radial direction by strong radiation force. As shown in Fig. 5, ram pressure of the neutral gas inflow is comparable to gas pressure of the ionized region before the transition. As a result, the neutral gas inflow pushes the shock layer from the equator to the polars. Panel (i2) in Fig. 4 shows that dense clumps form at the shock front by compression caused by continuous accretion of neutral gas and strong radiation force from the centre. Although some clumps are ejected by radiation force, such clump formation occurs more frequently in the late stage. Eventually, the central radiating region is totally covered by dense clumps of neutral gas. Figure 4. View largeDownload slide Two-dimensional distribution of gas density and velocity structure in the central region (rmin ≤ r ≲ 0.04RB) for Model E (MBH = 5 × 105 M⊙). Each panel corresponds to the time at (i1) t = 1.7 × 105 yr (≃0.04tdyn), (i2) t = 5.9 × 105 yr (≃0.14tdyn), and (i3) t = 2.1 × 106 yr (≃0.49tdyn). Figure 4. View largeDownload slide Two-dimensional distribution of gas density and velocity structure in the central region (rmin ≤ r ≲ 0.04RB) for Model E (MBH = 5 × 105 M⊙). Each panel corresponds to the time at (i1) t = 1.7 × 105 yr (≃0.04tdyn), (i2) t = 5.9 × 105 yr (≃0.14tdyn), and (i3) t = 2.1 × 106 yr (≃0.49tdyn). Figure 5. View largeDownload slide Profiles of ram pressure $$\rho v_{\theta }^2$$ (blue), gas pressure pgas (red), and the electron fraction (green) in the Northern hemisphere from the polar (θ = 0°) to the equator (θ = 90°) at r = 0.009RB in Model E. The elapsed time is the same as that of panel (i1) in Fig. 4. Figure 5. View largeDownload slide Profiles of ram pressure $$\rho v_{\theta }^2$$ (blue), gas pressure pgas (red), and the electron fraction (green) in the Northern hemisphere from the polar (θ = 0°) to the equator (θ = 90°) at r = 0.009RB in Model E. The elapsed time is the same as that of panel (i1) in Fig. 4. Figs 6–8 present radial profiles of the gas density, temperature, and radial velocity along three different directions: θ = 0° (polar), 16°, and 90° (equator) in the early stage at t = 0.14tdyn (blue) and in the late stage at t = 0.49tdyn (red). The black dashed curve represents the isothermal Bondi profile with T = 104K. In the early stage before the transition, radiation force and heating affect properties of the accreting material. Inside the ionized region, the gas is depleted and a density cavity is created (θ = 0° and 16°). Colliding outflows due to the radiation force on electrons (dashed curves in Fig. 8) produce strong shocks with a high temperature of ∼106 K. On the other hand, neutral gas accretes through the equator, keeping the temperature at T ≃ 8000 K, as shown in Section 3.1 (Model C). After the transition, the ionized gas becomes neutral due to radiative recombination and the density cavity is filled. The gas temperature settles down to ∼104 K in all directions because of radiative cooling by atomic hydrogen. Along the direction of θ = 16°, the ionized region still exists but its size is too small (∼0.01RB) to affect gas dynamics at the Bondi radius. As a result of absorption of anisotropic radiation by neutral hydrogen, the radiation momentum produces warm outflows with T ≃ 8000 K towards narrow conical regions along the polars (θ ∼ 0°–20°). Figure 6. View largeDownload slide Radial profiles of the gas density in Model E along three different directions of θ = 0° (polar), θ = 16°, and θ = 90° (equator) in the early stage at t = 0.14tdyn (blue) and in the late stage at t = 0.49tdyn (red). The black curve represents an isothermal Bondi profile with T = 104K. Figure 6. View largeDownload slide Radial profiles of the gas density in Model E along three different directions of θ = 0° (polar), θ = 16°, and θ = 90° (equator) in the early stage at t = 0.14tdyn (blue) and in the late stage at t = 0.49tdyn (red). The black curve represents an isothermal Bondi profile with T = 104K. Figure 7. View largeDownload slide Same as Fig. 6 but for the gas temperature. Figure 7. View largeDownload slide Same as Fig. 6 but for the gas temperature. Figure 8. View largeDownload slide Same as Fig. 6 but for the radial velocity. Solid and dashed curves present the inflow and outflow velocity, respectively. Figure 8. View largeDownload slide Same as Fig. 6 but for the radial velocity. Solid and dashed curves present the inflow and outflow velocity, respectively. Fig. 9 shows the time evolution of the inflow (green long-dashed) and outflow (blue short-dashed) rate measured at the Bondi radius, and the BH feeding rate (red solid). The normalized BH feeding rate of $$\dot{M}/\dot{M}_{\rm Edd}$$ continuously increases with the inflow rate from the Bondi radius. An open circle indicates the epoch of the transition to the wholly neutral accretion phase. After the transition, the feeding rate becomes highly episodic because strong radiation force and high ram pressure are tightly balanced at the innermost region. At the end of the simulation, the central BH is fed through the equator at the time-averaged rate of ∼5 × 104LEdd/c2, which corresponds to ∼ 50 per cent of the inflow rate from the Bondi radius. On the other hand, the accreting material is accelerated outward into the bipolar directions (θ ≲ 20°) due to absorption of the radiation momentum with ≫LEdd/c. The outflow rate ends up $${\simeq } 5\times 10^3 \dot{M}_{\rm Edd}$$, which is ∼ 10 per cent of the BH feeding rate. Since the Mach number of the outflowing matter is as high as four at the Bondi radius (see Fig. 8), the typical velocity is estimated as ≲ 40 km s − 1. Sugimura et al. (2016) reported that hot and ionized outflows with T ≃ 7 × 104 K are launched from a rapid accreting system exposed to strong anisotropic radiation, which is consistent with our cases where the transition does not occur as in Model C. On the other hand, our simulation in Model E suggests that warm and neutral outflows with T ≃ 8000K are likely to be produced because of the transition to the wholly neutral phase. Figure 9. View largeDownload slide Time evolution of inflow (green long-dashed) and outflow (blue short-dashed) rates at the Bondi radius and the BH feeding rate (red solid) in Model E (MBH = 5 × 105 M⊙). Figure 9. View largeDownload slide Time evolution of inflow (green long-dashed) and outflow (blue short-dashed) rates at the Bondi radius and the BH feeding rate (red solid) in Model E (MBH = 5 × 105 M⊙). In Fig. 10, we show the BH feeding rates for two cases with MBH = 5 × 104 M⊙ (blue, Model D) and 5 × 105 M⊙ (red, Model E) for comparison. The normalized accretion rate in Model D is saturated at a constant value unlike in Model E. In fact, we do not find a transition to a wholly neutral phase within the entire simulation time. Even without the transition, neutral gas inflows have non-radial motions and dense clumps are formed at the boundary between the ionized and neutral region as in Model E. However, ram pressure of the infalling clumps is too weak to cover the central radiation source. Therefore, this result implies the critical BH mass for the transition is located around MBH ∼ 5 × 105 M⊙ for the initial conditions we adopted (n∞ = 105 cm−3 and T∞ = 104 K). Figure 10. View largeDownload slide Time evolution of the accretion rates for two cases with MBH = 5 × 104 M⊙ (Model D) and MBH = 5 × 105 M⊙ (Model E). An open circle presents the epoch where the ionization region collapses by radiative recombination. The transition occurs for MBH ≳ 5 × 105 M⊙. Figure 10. View largeDownload slide Time evolution of the accretion rates for two cases with MBH = 5 × 104 M⊙ (Model D) and MBH = 5 × 105 M⊙ (Model E). An open circle presents the epoch where the ionization region collapses by radiative recombination. The transition occurs for MBH ≳ 5 × 105 M⊙. Transition to the wholly neutral phase is more likely to occur in more massive BH cases. The physical reason is that the size of the ionized region relative to the Bondi radius becomes smaller with the BH mass, following $$R_{\rm HII}/R_{\rm B}\propto M_{\rm BH}^{-3/2}$$, as found by IHO16. When the radiation flux is anisotropic as in our two-dimensional simulations, neutral gas accretes throughout the equatorial region with a significant fraction of solid angle (e.g. ΩHI/4π ≃ 0.7 for Model E) even if we do not assume shadow structures artificially. The size of the neutral region becomes larger with the BH mass, and eventually the transition occurs for the highest BH mass (Model E). 4 SUMMARY AND DISCUSSION In this paper, we study rapid growth of BHs with a mass range of 103 ≤ MBH ≤ 5 × 105 M⊙ via gas accretion exposed to anisotropic super-Eddington radiation. Radiation flux from an accretion disc around the nuclear BH is likely to be collimated towards the bipolar directions. To address the effects of the anisotropic radiation on the accretion flow, we perform two-dimensional radiation hydrodynamical simulations, including multifrequency radiation transfer and non-equilibrium primordial chemistry. We here consider a BH embedded in an initially uniform gas cloud with the density of n∞ = 105 cm−3 and temperature of T∞ = 104 K. We have found that hyper-Eddington accretion from the Bondi radius can occur when the radiation is more anisotropic and the BH mass is higher. Moreover, when the BH mass is as large as MBH ≳ 5 × 105 M⊙, ionized regions around the BH collapse due to super-sonic non-radial motions of neutral gas inflows and efficient radiative recombination. Due to absorption of the radiation momentum, the neutral gas is accelerated outward to narrow conical regions along the polars, producing warm outflows with T ≃ 8000 K. We adopt a simple model for anisotropic radiation from the central accretion disc, namely, Frad(θ) ∝ cos Nθ. As a result of the anisotropic radiation, rapid accretion is allowed from the equatorial plane. This assumption would be reasonable at the vicinity of the central BH because a super-Eddington accretion disc is optically and geometrically thick, so that the radiation flux is mildly collimated to the polars (e.g. Ohsuga et al. 2005; Jiang et al. 2014; Sa̧dowski et al. 2015). However, although the radiation propagates outward, such anisotropic radiation fields could be smoothed out because of radiative diffusion. In order to explore the anisotropy of the emergent radiation, therefore radiation-hydrodynamical simulations with a wide dynamical range are required. Recently, Sugimura et al. (2016) have performed two-dimensional simulations of accreting BHs under anisotropic radiation and found that anisotropic radiation allows gas accretion through the equatorial plane. Although our calculations have been pursued independently from their work, our simulation set-ups are similar, mainly except two things. First, they assumed the existence of a shadow region caused by attenuation due to a geometrically thick disc and/or disc winds well inside their computation domains. Secondly, our simulations resolve even smaller scales, namely ∼ 1 per cent of the Bondi radius. We can stress that our set-up is more preferable to follow non-radial motions of dense and neutral inflows, which lead to a transition to a wholly neutral accretion phase. Conversely, Sugimura et al. (2016) did not find such transitions, probably because (1) the BH mass adopted in their simulations is lower than that in our Model E, and (2) the radius of the inner most grid in their simulations is larger than that in our simulations. We discuss the subsequent evolution of the accretion flows in Model E (MBH = 105.7 M⊙). The gas profile approaches an isothermal Bondi solution with T ∼ 104 K, where $$\dot{M}_{\rm B}/\dot{M}_{\rm Edd}\simeq 3 \times 10^5 n_{\infty .5}M_{\rm BH,5.7}$$ if the ram pressure of the neutral inflow dominates the radiation force. This condition is written as $$L/L_{\rm Edd} \lesssim 139 M_{\rm BH,5.7}^{3/2}n_{\infty ,5} (r_{\star }/7\times 10^{17}\,{\rm cm})^{-1/2}$$ (Sakurai et al. 2016), where r⋆ is the location of the photosphere, from which the radiation emerges actually. Note that the photosphere is located inside the innermost grid (r⋆ < rmin = 7 × 1017 cm). Assuming that the BH is fed at the Bondi rate, the radiation luminosity towards the polar directions is estimated as L/LEdd ≈ 106 (see equation 12). Therefore, we expect that the overall behaviours of the accretion flows in Model E would be maintained. We set the ambient gas density to n∞ = 105 cm−3. If PopIII stars emit ionizing photons, the gas inside the host dark matter halo with ∼105–6 M⊙ (the so-called mini-halo, z ∼ 20) is evaporated by the irradiation, reducing the density much below that we assume (e.g. Kitayama et al. 2004; Johnson & Bromm 2007). However, when the halo grows to be as massive as ≳ 107 M⊙ with a virial temperature of ≳ 104 K, (the so-called atomic-cooling halo, z ∼ 10), the gas re-collapses on to the halo, and a density profile with a core and an envelope following ρ ∝ r−2 is built up due to Lyman α cooling. Alternatively, when remnant BHs are located at the centre of an atomic-cooling halo, the gas density would be high enough to realize the conditions we assume. The massive BHs we consider in Models D and E would be formed by, for example, direct-collapse of supermassive stars in atomic-cooling haloes. Since a supermassive star has a bloated cold envelope with ≃5000 K, its radiation feedback would be too weak to photoevaporate a surrounding medium (e.g. Chon, Hosokawa & Yoshida 2018). The supermassive star finally collapses to an equally massive BH (∼104–5 M⊙) without mass-loss and explosions. Thus, the gas density profile (ρ ∝ r−2) in the halo would not be broken before the formation of the BH seed (see also IHO16). As shown in Fig. 4, non-radial motions of neutral inflows near the inner boundary are essential for leading to the transition to a wholly neutral phase. However, our simulations neglect non-radial components in radiation flux $$F_{\rm rad}^{\theta }$$, which can be produced by radiative recombination in the ionized region.1 The radiation flux due to recombination is estimated as $$F_{\rm rad}^{\theta } \sim 4{\pi} \eta r_{\rm min}$$, where $$\eta \sim \langle h\nu _{\rm rec} \rangle \alpha _{\rm B}n_{\rm e}^2 \sim 5\times 10^{-17}(n_{\rm e}/10^4\, {\rm cm}^{-3})^2\,{\rm erg\,cm^{-3}\,s^{-1}}$$ is the emissivity, 〈hνrec〉 is the average energy of recombination photons, αB is the case B recombination rate coefficient, and ne is the number density of electrons. Thus, the non-radial radiation flux is estimated as $$|F_{\rm rad}^{\theta }/F_{\rm rad}^r| \sim 10^{-7}$$ at the innermost region, and the radiation force in non-radial directions is negligible. Throughout this paper, we neglect formation of molecular hydrogen (H2), which is the main coolant in primordial gas below ∼104 K. We briefly mention that H2 cooling is unlikely to affect the properties of warm/neutral outflows. In a warm and dense primordial gas cloud (n ≳ 105 cm − 3 and T ∼ 8000 K), the H2 fraction relative to atoms is determined by the balance between formation through electron-catalysed reactions and collisional dissociation with gas particles (Omukai 2001). The resultant fraction is estimated as ∼10−8 as long as the temperature is kept at >4000 K, below which the collisional dissociation becomes inefficient (Inayoshi et al. 2014). As shown in Fig. 7, the temperature cools down to ≃1000–2000 K due to adiabatic expansion at r > 2 × 1019 cm in the outflowing matter (θ = 0°). In this lower temperature region, far-ultraviolet radiation in the Lyman-Werner (LW) bands (hνLW = 10.2–13.6 eV) from the central region plays an important role to determine the H2 fraction. We can estimate the LW flux in the optically thin limit as   \begin{eqnarray} F_{\rm LW,thin}\, = \, f_{\rm LW}(N+1)\frac{2\left[1+\ln (\dot{m}/20)\right]\,L_{\rm Edd}}{4{\pi} r^2\langle \nu _{\rm LW}\rangle }\nonumber \\ \simeq \, 1.2\times 10^{-11} \,{\rm erg \,s^{-1}\, cm^{-2} \,Hz^{-1}} \nonumber\\ \times \left(\frac{f_{\rm LW}}{0.3}\right) \left(\frac{N+1}{5}\right) \!\left(\frac{M_{\rm BH}}{5\times 10^5 \,{\rm M}_{\odot }}\right)^{-1} \!\left(\frac{r}{R_{\rm B}}\right)^{-2}, \,\,\,\nonumber\\ \end{eqnarray} (13)where fLW is the ratio of the LW flux to the bolometric value and 〈νLW〉 = 12eV/h( = 2.9 × 1015Hz). Since the H2 column density is $$N_{\rm H_2}\sim 10^{17} {\rm cm}^{-2}(n_\infty /10^5\, {\rm cm}^{-3})(r/R_{\rm B})$$ at most, the emergent LW flux can be reduced by the H2 self-shielding effect. The shielding factor is simply estimated as $$f_{\rm shield}\simeq (N_{\rm H_2}/10^{14}\,{\rm cm}^{-2})^{-3/4}$$ for a higher column density (Draine & Bertoldi 1996, see also Wolcott-Green, Haiman & Bryan 2011), and thus FLW = fshieldFLW, thin ≃ 7 × 10−14 erg s−1 cm−2 Hz−1, which is higher than the critical LW flux,2FLW, crit ∼ 10−19–10−16 erg s−1 cm−2 Hz−1 (e.g. Sugimura, Omukai & Inoue 2014; Inayoshi & Tanaka 2015; Latif, Schleicher & Hartwig 2016; Wolcott-Green, Haiman & Bryan 2017) required to dissociate H2 and suppress H2 cooling in the warm and neutral gas even at r ∼ 1020 cm. We assume the Ly α cooling (2p → 1s state) in the optically thin limit. In reality, Ly α photons are frequently scattered in neutral gas before escaping from the accretion system. In such a cloud where Ly α photons are trapped, two-photon emission (2s → 1s state) and H− free–bound emission (H + e− → H− + γ) dominate as radiative cooling processes (e.g. Omukai 2001; Schleicher, Spaans & Glover 2010). As shown by the previous studies, thermal evolution of a collapsing gas does not change significantly from that assuming the optically thin Ly α cooling (see also Spaans & Silk 2006; Ge & Wise 2017). Furthermore, additional cooling processes, e.g. fine-structure lines of C ii and O i, and thermal emission by dust grains would be important in a metal-enriched gas cloud. The metallicity effects on hyper-Eddington accretion are left for future investigations. Acknowledgements We thank Zoltán Haiman and Kazuyuki Sugimura for useful discussions. This work is partially supported by the Simons Foundation through the Simons Society of Fellows (KI) and by Japan Society for the Promotion of Science Grant-in-Aid for Scientific Research (A) (17H01102 KO), for Scientific Research (C) (15K05036 KO; 16K05309 KO and 17K05383 SM)) and for Young Scientists (17K14260 HRT). This research was also supported by the Ministry of Education, Culture, Sports, Science and Technology as ‘Priority Issue on Post-K computer’ (Elucidation of the Fundamental Laws and Evolution of the Universe) and the Joint Institute for Computational Fundamental Science. Numerical computations were carried out on Cray XC30 at Center for Computational Astrophysics, National Astronomical Observatory of Japan. Footnotes 1 The non-radial components in radiation flux is unlikely to be produced by electron scattering because the ionized gas is optically thin. In fact, before the transition in Model E, the optical depth is as low as 0.16 at r < RHII along a line of θ = 10°, which is inside the ionized region. 2 The critical LW flux depends on the radiation spectrum, namely on the ratio of the flux at 2 eV to the LW flux (e.g. Sugimura, Omukai & Inoue 2014). This is because photons with ∼2 eV dissociate H− ions, which is an intermediate product in the H2 formation channels (H + e− → H− + γ; H− + H → H2 + e−) and contribute to H2 dissociation indirectly (Omukai 2001). REFERENCES Abel T., Anninos P., Zhang Y., Norman M. L., 1997, New Astron. , 2, 181 https://doi.org/10.1016/S1384-1076(97)00010-9 CrossRef Search ADS   Alexander T., Natarajan P., 2014, Science , 345, 1330 https://doi.org/10.1126/science.1251053 CrossRef Search ADS PubMed  Alvarez M. A., Wise J. H., Abel T., 2009, ApJ , 701, L133 CrossRef Search ADS   Anninos P., Zhang Y., Abel T., Norman M. L., 1997, New Astron. , 2, 209 https://doi.org/10.1016/S1384-1076(97)00009-2 CrossRef Search ADS   Becerra F., Greif T. H., Springel V., Hernquist L. E., 2015, MNRAS , 446, 2380 CrossRef Search ADS   Begelman M. C., 1979, MNRAS , 187, 237 CrossRef Search ADS   Bondi H., 1952, MNRAS , 112, 195 CrossRef Search ADS   Chon S., Hirano S., Hosokawa T., Yoshida N., 2016, ApJ , 832, 134 CrossRef Search ADS   Chon S., Hosokawa T., Yoshida N., 2018, MNRAS , in press Ciotti L., Ostriker J. P., 2001, ApJ , 551, 131 CrossRef Search ADS   Ciotti L., Ostriker J. P., Proga D., 2009, ApJ , 699, 89 CrossRef Search ADS   Devecchi B., Volonteri M., 2009, ApJ , 694, 302 CrossRef Search ADS   Draine B. T., Bertoldi F., 1996, ApJ , 468, 269 CrossRef Search ADS   Fan X. et al.  , 2004, AJ , 128, 515 CrossRef Search ADS   Ge Q., Wise J., 2017, MNRAS , 472, 2773 Glover S. C. O., Jappsen A.-K., 2007, ApJ , 666, 1 CrossRef Search ADS   Haiman Z., Loeb A., 2001, ApJ , 552, 459 CrossRef Search ADS   Harten A., Lax P. D., van Leer B., 1983, SIAM Rev. , 25, 35 https://doi.org/10.1137/1025002 CrossRef Search ADS   Inayoshi K., Tanaka T. L., 2015, MNRAS , 450, 4350 CrossRef Search ADS   Inayoshi K., Omukai K., Tasker E., 2014, MNRAS , 445, L109 CrossRef Search ADS   Inayoshi K., Haiman Z., Ostriker J. P., 2016, MNRAS , 459, 3738 (IHO16) CrossRef Search ADS   Jiang Y.-F., Stone J. M., Davis S. W., 2014, ApJ , 796, 106 CrossRef Search ADS   Johnson J. L., Bromm V., 2007, MNRAS , 374, 1557 CrossRef Search ADS   Kato Y., Mineshige S., Shibata K., 2004, ApJ , 605, 307 CrossRef Search ADS   Katz H., Sijacki D., Haehnelt M. G., 2015, MNRAS , 451, 2352 CrossRef Search ADS   King A., 2003, ApJ , 596, L27 CrossRef Search ADS   Kitayama T., Yoshida N., Susa H., Umemura M., 2004, ApJ , 613, 631 CrossRef Search ADS   Kormendy J., Ho L. C., 2013, ARA&A , 51, 511 CrossRef Search ADS   Latif M. A., Schleicher D. R. G., Hartwig T., 2016, MNRAS , 458, 233 CrossRef Search ADS   Lupi A., Haardt F., Dotti M., Fiacconi D., Mayer L., Madau P., 2016, MNRAS , 456, 2993 CrossRef Search ADS   Madau P., Rees M. J., 2001, ApJ , 551, L27 CrossRef Search ADS   Madau P., Haardt F., Dotti M., 2014, ApJ , 784, L38 CrossRef Search ADS   Milosavljević M., Bromm V., Couch S. M., Oh S. P., 2009, ApJ , 698, 766 CrossRef Search ADS   Mortlock D. J. et al.  , 2011, Nature , 474, 616 CrossRef Search ADS PubMed  Murray N., Quataert E., Thompson T. A., 2005, ApJ , 618, 569 CrossRef Search ADS   Ohsuga K., Mineshige S., 2007, ApJ , 670, 1283 CrossRef Search ADS   Ohsuga K., Mineshige S., 2011, ApJ , 736, 2 CrossRef Search ADS   Ohsuga K., Mori M., Nakamoto T., Mineshige S., 2005, ApJ , 628, 368 CrossRef Search ADS   Ohsuga K., Mineshige S., Mori M., Kato Y., 2009, PASJ , 61, L7 CrossRef Search ADS   Omukai K., 2001, ApJ , 546, 635 CrossRef Search ADS   Pacucci F., Ferrara A., 2015, MNRAS , 448, 104 CrossRef Search ADS   Park K., Ricotti M., 2011, ApJ , 739, 2 CrossRef Search ADS   Park K., Ricotti M., 2012, ApJ , 747, 9 CrossRef Search ADS   Pezzulli E., Valiante R., Schneider R., 2016, MNRAS , 458, 3047 CrossRef Search ADS   Regan J. A., Johansson P. H., Haehnelt M. G., 2014, MNRAS , 439, 1160 CrossRef Search ADS   Ryu T., Tanaka T. L., Perna R., Haiman Z., 2016, MNRAS , 460, 4122 CrossRef Search ADS   Sakurai Y., Inayoshi K., Haiman Z., 2016, MNRAS , 461, 4496 CrossRef Search ADS   Sakurai Y., Yoshida N., Fujii M. S., Hirano S., 2017, MNRAS , 472, 1677 Schleicher D. R. G., Spaans M., Glover S. C. O., 2010, ApJ , 712, L69 CrossRef Search ADS   Shakura N. I., Sunyaev R. A., 1973, A&A , 24, 337 Silk J., Rees M. J., 1998, A&A , 331, L1 Soltan A., 1982, MNRAS , 200, 115 CrossRef Search ADS   Spaans M., Silk J., 2006, ApJ , 652, 902 CrossRef Search ADS   Stone N. C., Küpper A. H. W., Ostriker J. P., 2017, MNRAS , 467, 4180 Sugimura K., Omukai K., Inoue A. K., 2014, MNRAS , 445, 544 CrossRef Search ADS   Sugimura K., Hosokawa T., Yajima H., Omukai K., 2016, MNRAS , 469, 62 CrossRef Search ADS   Sa̧dowski A., Narayan R., Tchekhovskoy A., Abarca D., Zhu Y., McKinney J. C., 2015, MNRAS , 447, 49 CrossRef Search ADS   Takahashi H. R., Ohsuga K., 2013, ApJ , 772, 127 CrossRef Search ADS   Takahashi H. R., Ohsuga K., Kawashima T., Sekiguchi Y., 2016, ApJ , 826, 23 CrossRef Search ADS   Tanaka T., Haiman Z., 2009, ApJ , 696, 1798 CrossRef Search ADS   Valiante R., Schneider R., Volonteri M., Omukai K., 2016, MNRAS , 457, 3356 CrossRef Search ADS   van Leer B., 1977, J. Comput. Phys. , 23, 263 https://doi.org/10.1016/0021-9991(77)90094-8 CrossRef Search ADS   Visbal E., Haiman Z., Bryan G. L., 2014, MNRAS , 445, 1056 CrossRef Search ADS   Volonteri M., Rees M. J., 2005, ApJ , 633, 624 CrossRef Search ADS   Volonteri M., Haardt F., Madau P., 2003, ApJ , 582, 559 CrossRef Search ADS   Watarai K.-y., Fukue J., Takeuchi M., Mineshige S., 2000, PASJ , 52, 133 CrossRef Search ADS   Whalen D., Norman M. L., 2006, ApJS , 162, 281 CrossRef Search ADS   Whalen D. J., Norman M. L., 2008, ApJ , 672, 287 CrossRef Search ADS   Wolcott-Green J., Haiman Z., Bryan G. L., 2011, MNRAS , 418, 838 CrossRef Search ADS   Wolcott-Green J., Haiman Z., Bryan G. L., 2017, MNRAS , 469, 3329 Wu X.-B. et al.  , 2015, Nature , 518, 512 CrossRef Search ADS PubMed  Yajima H., Khochfar S., 2016, MNRAS , 457, 2423 CrossRef Search ADS   Yu Q., Tremaine S., 2002, MNRAS , 335, 965 CrossRef Search ADS   © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Monthly Notices of the Royal Astronomical Society Oxford University Press

Rapid growth of black holes accompanied with hot or warm outflows exposed to anisotropic super-Eddington radiation

Loading next page...
 
/lp/ou_press/rapid-growth-of-black-holes-accompanied-with-hot-or-warm-outflows-muaGyMzucO
Publisher
The Royal Astronomical Society
Copyright
© 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society
ISSN
0035-8711
eISSN
1365-2966
D.O.I.
10.1093/mnras/sty264
Publisher site
See Article on Publisher Site

Abstract

Abstract We perform two-dimensional radiation hydrodynamical simulations of accretion flows on to a black hole (BH) with a mass of 103 ≤ MBH/ M⊙ ≲ 106 in order to study rapid growth of BHs in the early Universe. For spherically symmetric flows, hyper-Eddington accretion from outside the Bondi radius can occur unimpeded by radiation feedback when MBH ≳ 104 M⊙(n∞/105 cm − 3) − 1(T∞/104 K)3/2, where the density and temperature of ambient gas are initially set to n∞ = 105 cm−3 and T∞ = 104 K. Here, we study accretion flows exposed to anisotropic radiation from a nuclear accretion disc with a luminosity higher than the Eddington value (LEdd) due to collimation towards the bipolar directions. We find that, unlike the spherically symmetric case, even less massive BHs with MBH < 104 M⊙ can be fed at high accretion rates of ≳ LEdd/c2 through the equatorial region, while ionized regions expand towards the poles producing hot outflows with T ∼ 105 K. For more massive BHs with MBH ≳ 5 × 105 M⊙, intense inflows of neutral gas through the equator totally cover the central radiating region due to the non-radial gas motions. Because of efficient recombination by hydrogen, the entire flow settles in neutral and warm gas with T ≃ 8000 K. The BH is fed at a rate of ∼5 × 104LEdd/c2 (a half of the inflow rate from the Bondi radius). Moreover, radiation momentum absorbed by neutral hydrogen produces warm outflows towards the bipolar directions at ∼ 10 per cent of the BH feeding rate and with a velocity several times higher than the escaping value. black hole physics, quasars: supermassive black holes, cosmology: theory 1 INTRODUCTION Observations of bright quasars at high redshifts (z ≳ 6) have revealed that supermassive black holes (SMBHs) with MBH ≳ 109 M⊙ are formed within ≲ 1 Gyr after the beginning of the Universe (e.g. Fan et al. 2004; Mortlock et al. 2011; Wu et al. 2015). Although SMBHs play crucial roles in the co-evolution with their host galaxies through feedback processes (e.g. Silk & Rees 1998; King 2003; Murray, Quataert & Thompson 2005; Kormendy & Ho 2013), their formation and growth processes are still open questions. One possible formation scenario of high-z SMBHs is that stellar mass black holes (BHs) formed by collapse of first generation stars (Population III, hereafter PopIII) with ∼100 M⊙ grow via gas accretion (e.g. Madau & Rees 2001; Haiman & Loeb 2001; Volonteri, Haardt & Madau 2003). In the standard accretion physics, BH growth would be suppressed by radiation emitted from a nuclear accretion disc because the radiation luminosity increases with the accretion rate on to the BH. If the radiative efficiency of the disc is η ≃ 0.1 (Shakura & Sunyaev 1973; Soltan 1982; Yu & Tremaine 2002), the radiation luminosity exceeds the Eddington luminosity (LEdd) at the accretion rate of ≳ LEdd/(ηc2), above which the radiation force dominates the gravity of the central BH. Assuming this Eddington-limited accretion rate, it takes ≳ 1 Gyr for stellar-mass seed BHs to grow up to ∼109 M⊙. This fact requires more rapid accretion mechanisms to form SMBHs by z ∼ 6, when the corresponding cosmic age is ∼0.8 Gyr. As alternative possibilities, formation of massive seed BHs with ≳ 103-105 M⊙ has been proposed: direct-collapse of supermassive stars in protogalaxies (e.g. Inayoshi, Omukai & Tasker 2014; Regan, Johansson & Haehnelt 2014; Visbal, Haiman & Bryan 2014; Becerra et al. 2015; Inayoshi & Tanaka 2015; Chon et al. 2016; Latif, Schleicher & Hartwig 2016) or runaway collisions of stars in a dense cluster (e.g. Devecchi & Volonteri 2009; Katz, Sijacki & Haehnelt 2015; Yajima & Khochfar 2016; Stone, Küpper & Ostriker 2017; Sakurai et al. 2017). However, since the BH growth time even from such massive seeds would be shorten only by a factor of 2 in the case of the Eddington limited accretion, a high duty cycle ( ∼ 100 per cent) is still required to explain the existence of SMBHs with ∼109 M⊙ at high redshifts. The possibility of rapid growth of seed BHs has been discussed (e.g. Volonteri & Rees 2005; Tanaka & Haiman 2009; Alexander & Natarajan 2014; Madau, Haardt & Dotti 2014; Pacucci & Ferrara 2015; Pezzulli, Valiante & Schneider 2016; Valiante et al. 2016). By means of radiation hydrodynamical simulations of accretion flows on to a BH, Ohsuga et al. (2005) concluded that super-Eddington accretion can be realized as long as sufficient gas is supplied at the vicinity of the central BH (see also Ohsuga et al. 2009; Ohsuga & Mineshige 2011; Jiang, Stone & Davis 2014; Sa̧dowski et al. 2015; Takahashi et al. 2016). This is because diffusive photons are efficiently trapped within optically thick accretion flows and the radiation force on to the gas is alleviated (e.g. Begelman 1979). However, it is unclear how sufficient amount of gas is supplied into the nuclear region. In fact, radiation heating and ionization effectively suppress gas accretion from larger scales (Ciotti & Ostriker 2001; Alvarez, Wise & Abel 2009; Ciotti, Ostriker & Proga 2009; Milosavljević et al. 2009; Park & Ricotti 2011, 2012). Inayoshi, Haiman & Ostriker (2016) (hereafter IHO16) have found a self-consistent solution of steady hyper-Eddington accretion in spherically symmetric systems with radiation hydrodynamical simulations. This accretion flow has a high accretion rate of ≳ 5000LEdd/c2, unimpeded by negative feedback due to radiation force and heating simultaneously. The solution is composed of a radiation-dominated core, where photon trapping due to electron scattering works efficiently, and an accreting envelope that follows a Bondi profile with 8000 K. When the emergent radiation luminosity is limited to ≲ LEdd because of photon trapping, radiation from the central region does not affect the gas dynamics at larger scales. The hyper-Eddington accretion is then realized when a BH is embedded in a dense gas cloud, for which the following condition is satisfied,   \begin{equation} M_{\rm BH} \gtrsim 10^4 \,{\rm M}_{\odot }\left(\frac{n_\infty }{10^5 \,{\rm cm}^{-3}}\right)^{-1}\left(\frac{T_\infty }{10^4\, {\rm K}}\right)^{3/2}, \end{equation} (1)where n∞ and T∞ are density and temperature of the ambient gas, respectively. Hereafter, we define MBH, x = (MBH/10x M⊙), n∞, 5 = (n∞/105 cm−3), and T∞, 4 = (T∞/104 K). Note that equation (1) is identical to the condition where the Bondi radius is larger than the size of the ionized region. Even if a high luminosity with >LEdd emerges from the centre, the above conditions are not affected as long as the luminosity is as low as (10 − 30) × LEdd (Sakurai, Inayoshi & Haiman 2016). To satisfy the condition of equation (1), gas-rich regions with n ≳ 107 cm − 3 are required for stellar-mass PopIII BHs with masses of ≲ 100 M⊙. During the assembly of a protogalaxy, some PopIII BHs would fall into a gaseous circumnuclear disc, where the gas accretion rate is likely to be ≳ 5000LEdd/c2 (Ryu et al. 2016, see also Lupi et al. 2016). Radiation flux from an accretion disc around the nuclear BH is more likely to be anisotropic, i.e. preferentially collimated around the rotation axis of the disc (Ohsuga et al. 2009; Ohsuga & Mineshige 2011; Jiang et al. 2014; Sa̧dowski et al. 2015; Takahashi et al. 2016). On the other hand, the disc-like accretion could reduce the feedback effect because intense radiation is significantly attenuated by the dense accreting material in the disc (Ohsuga & Mineshige 2007; Sugimura et al. 2016). In this paper, we investigate how intermediate massive BHs with 103 ≤ MBH ≤ 5 × 105 M⊙, embedded in an initially uniform gas with n∞ = 105 cm−3 and T∞ = 104 K, can grow via accretion when the inflow gas is exposed to the super-Eddington, anisotropic radiation. To address this issue, we perform two-dimensional radiation hydrodynamical simulations, including multifrequency radiation transfer and non-equilibrium primordial chemistry. We confirm that hyper-Eddington accretion from the Bondi radius is likely to occur when radiation is more anisotropic and the BH mass is higher. Furthermore, we demonstrate that for MBH ≳ 5 × 105 M⊙, neutral gas inflows through the equator totally cover the central radiating region due to the non-radial gas motions, and thus the emergent radiation in all directions is blocked. Because of efficient recombination by hydrogen, the entire flow results in neutral and warm gas. In this case, the BH feeding rate is as high as the Bondi accretion rate, and also warm outflows are produced by absorption of radiation momentum. The rest of this paper is organized as follows. In Section 2, we describe the methodology for radiation hydrodynamical simulations. In Section 3, we show our simulation results and discuss the effects of anisotropic radiation. Finally, we summarize the main conclusions of this paper and discuss caveats of our simulations in Section 4. 2 SIMULATION METHOD We perform two-dimensional hydrodynamical simulations including one-dimensional radiation transfer and chemical reaction networks. Our purpose is to study necessary conditions for hyper-Eddington accretion led by sufficient supply of gas from large scales. Gas accretion begins from the Bondi radius, within which the gravity of the central BH dominates the gas-pressure gradient force, given by Bondi (1952) as   \begin{equation} R_{\rm B}\equiv \frac{GM_{\rm BH}}{c_{\infty }^2} \simeq 1.97\times 10^{18} M_{\rm BH,4} T_{\infty ,4}^{-1} \ {\rm cm}, \end{equation} (2)where $$c_{\infty }=\sqrt{\gamma \mathcal {R} T_{\infty } /\bar{\mu }}$$ is the sound speed, γ is the specific heat ratio, $$\mathcal {R}$$ is the gas constant, and $$\bar{\mu }$$ is the mean molecular weight. As a reference of the accretion rate, we define the Bondi accretion rate for γ = 1 as   \begin{equation} \dot{M}_{\rm B} \equiv {\pi} e^{3/2} \rho _\infty \frac{G^2 M_{\rm BH}^2}{c_{\infty }^3}, \end{equation} (3)and the Eddington accretion rate as $$\dot{M}_{\rm Edd} \equiv L_{\rm Edd}/c^2$$, where LEdd = 4πcGMBH/κes and κes is the electron scattering opacity. Note that the Bondi radius and rate are calculated for γ = 1, $$\bar{\mu }=1.23$$, and T∞ = 104 K as a reference value through this paper, and that these are self-consistently calculated in our simulations. 2.1 Basic equations We consider a gas sphere exposed to intense radiation from the accretion flow on to the central BH, solving two-dimensional hydrodynamical equations assuming axisymmetric flows. We here use the hydrodynamical simulation code taken from Takahashi & Ohsuga (2013). The mass flux for the ideal fluid is computed using the Harten–Lax–vanLeer Riemann solver (Harten, Lax & van Leer 1983), and the second-order accuracy in space and time is ensured (van Leer 1977). We adopt the spherical coordinates of (r, θ, ϕ), defining the polar axis (θ = 0 and π) as directions perpendicular to the disc plane. The basic equations of hydrodynamics are the following: the equation of continuity   \begin{equation} \frac{\partial \rho }{\partial t} + \nabla \cdot (\rho {\boldsymbol v}) = 0, \end{equation} (4)and the equations of motion,   \begin{eqnarray} \frac{\partial (\rho v_r)}{\partial t} + \nabla \cdot (\rho v_r {\boldsymbol v}) = -\frac{\partial p}{\partial r} + \rho \left( \frac{v_{\theta }^2}{r} + \frac{v_{\phi }^2}{r} \right) - \rho \frac{\partial \psi}{\partial r}+ f_{\rm rad}, \nonumber\\ \end{eqnarray} (5)  \begin{equation} \frac{\partial \left(\rho rv_{\theta }\right)}{\partial t} + \nabla \cdot (\rho r v_{\theta } {\boldsymbol v}) = -\frac{\partial p}{\partial \theta } + \rho v_{\phi }^2\cot {\theta }, \end{equation} (6)  \begin{equation} \frac{\partial \left(\rho rv_{\phi }\sin {\theta }\right)}{\partial t} + \nabla \cdot \left(\rho r v_{\phi }\sin {\theta } {\boldsymbol v}\right) = 0, \end{equation} (7)where ρ is the gas density, $${\boldsymbol v}=(v_r,v_{\theta }, v_{\phi })$$ is the velocity, p is the gas pressure, and frad is the radiation force. We consider the gravity of the central BH (r = 0) and neglect the gas self-gravity. Since the general relativistic effect is negligible, the gravitational potential is given by ψ = −GMBH/r. We also solve the energy equation   \begin{equation} \frac{\partial e}{\partial t} + \nabla \cdot [(e+p){\boldsymbol v}] = -\frac{GM_{\rm BH}\rho }{r^2}v_r -\Lambda + \Gamma , \end{equation} (8)where $$e=e_{\rm int}+\rho |{\boldsymbol v}|^2/2$$ is the gas total energy density, eint is the gas internal energy density, Λ is the net cooling rate per volume, and Γ is the heating rate due to radiation from the central region. We assume the equation of state of ideal gas as p = (γ − 1)eint for γ = 5/3. We consider radiative cooling by bound–bound transitions of H, He, He+ atoms and free–free emission (Glover & Jappsen 2007). To estimate their rates, we solve chemical reaction networks including six species of H, H+, He, He+, He++, and e−. The abundance of He nuclei relative to H nuclei is set to 8.33 × 10−2. Here, photoionization, collisional ionization, and radiative recombination are considered (Abel et al. 1997; Glover & Jappsen 2007). Instead of treating photoionization by diffusive recombination photons, we adopt the on-the-spot approximation, where the case A recombination rate is replaced by the case B rate. To update the chemical abundances stably, we adopt a semi-implicit method (Anninos et al. 1997), setting the time-steps shorter than chemical time-scales for all the cells defined by $$t_{\rm chem} \equiv (x_{\rm e}+0.001x_{\rm H})/\dot{x}_{\rm e}$$, where xe and xH are the abundance of electrons and neutral hydrogens, respectively (Whalen & Norman 2006, 2008). Note that the cooling/heating part of the energy equation (equation 8) is updated with an implicit method in order to solve it stably and save computation time. We solve the multifrequency, steady radiative transfer equation because the light crossing time is much shorter than the hydrodynamical time-scale,   \begin{equation} \frac{1}{r^2}\frac{{\rm d}}{{\rm d}r}(r^2F_{\nu }) = -\rho \kappa _{\nu }cE_{\nu }, \end{equation} (9)where Fν is the radiation flux, Eν is the radiation energy density, and κν is the absorption opacity. The frequency range is set to hνmin(=13.6eV) ≤ hν ≤ hνmax(=10keV), where h is the Planck constant. Note that the radial component of the radiation flux is calculated (see Section 4 for more details). Since the ionized gas is optically thin to electron scattering and bound-free transitions, we assume Fν ≈ cEν on the right-hand side of equation (9). Most ionizing photons are absorbed by neutral hydrogen at the edge of the ionized region. Non-radial components of diffusive photons due to recombination are neglected in our simulation (see also Section 4). The ionization rate coefficients kph of H, He, He+, and photoionization heating rate Γ(=ΓH) are calculated following the photon-conserving manner (Whalen & Norman 2006). Here, ΓH is the heating rate due to H ionization. We consider radiation force due to electron scattering and bound-free absorption of H atoms as   \begin{equation} f_{\rm rad} = \frac{nx_{\rm e}}{c}\int _{\nu _{\rm min}}^{\nu _{\rm max}}\sigma _{\rm es}F_{\nu }{\rm d}\nu + \frac{\Gamma _{\rm H}}{c}. \end{equation} (10) In order to integrate equation (9), we simply set a radiation source with a single power-law spectrum at the centre, Lν = L0(ν/νmin)−α at νmin ≤ ν ≤ νmax and Lν = 0 at ν ≤ νmin and νmax ≤ ν, where α = 1.5 is adopted. The normalization factor L0 is set by the total luminosity ($$L=\int _{\nu _{\rm min}}^{\nu _{\rm max}} L_{\nu } {\rm d}\nu$$). We set a model for radiation luminosity emitted by an accretion disc, which is unresolved, at the inner boundary,   \begin{equation} \frac{\displaystyle L}{\displaystyle L_{\rm Edd}}= \begin{array}{ll}2\left[1+\ln {(\dot{m}/20)}\right] & {\rm for}\ \dot{m}\ge 20, \\ \dot{m}/10 & {\rm for}\ \dot{m}<20, \end{array} \end{equation} (11)(Watarai et al. 2000), where $$\dot{m} \equiv \dot{M}/\dot{M}_{\rm Edd}$$. Furthermore, we assume anisotropic radiation fields as   \begin{equation} F_{\nu }(r=r_{\rm min}, \theta ) = \frac{\displaystyle (N+1)L_{\nu }}{\displaystyle 4{\pi} r_{\rm min}^2}\cos ^N{\theta }, \end{equation} (12)where the radiation flux Fν(rmin) is normalized so that $$L_{\nu }(r_{\rm min})=\int F_{\nu }(r_{\rm min}) r_{\rm min}^2 {\rm d}\Omega$$ and rmin is the radius of the inner most grid in the computation domain (see Section 2.2). Unlike Sugimura et al. (2016), shadow regions are not assumed in our simulations (see Section 4 for more details). Here, we explore the effects of radiation anisotropy; the isotropic case N = 0 (Model A), the anisotropic cases N = 2 (Model B), and N = 4 (Model C) for MBH = 103 M⊙. We also study cases for higher BH masses with MBH = 5 × 104 M⊙ (Model D) and 5 × 105 M⊙ (Model E) for N = 4. Our models are summarized in Table 1. Table 1. Model  MBH( M⊙)  N  $$\dot{M}/\dot{M}_{\rm Edd}$$  L/LEdd  A  103  0  ∼0.9  ∼0.09  B  103  2  ≈2.7  ≈0.27  C  103  4  ≈44  ≈3.5  D  5 × 104  4  ∼2 × 103  ∼11  E  5 × 105  4  ∼5 × 104  ∼18  Model  MBH( M⊙)  N  $$\dot{M}/\dot{M}_{\rm Edd}$$  L/LEdd  A  103  0  ∼0.9  ∼0.09  B  103  2  ≈2.7  ≈0.27  C  103  4  ≈44  ≈3.5  D  5 × 104  4  ∼2 × 103  ∼11  E  5 × 105  4  ∼5 × 104  ∼18  Column: (1) model ID, (2) BH mass, (3) radiation anisotropy, namely Frad(θ) ∝ cos Nθ, (4) accretion rate, and (5) corresponding luminosity. In columns 4 and 5, we show the time-averaged values in Model A, the values at the end of simulations in Models B–D, and the time-averaged values over the last one cycle of episodic accretion in Model E. View Large 2.2 Initial and boundary conditions We set a computational domain of rmin ≤ r ≤ rmax and 0 ≤ θ ≤ π. Here, we adopt rmin = 0.04RB, rmax = 30RB for Models A, B, and C, and rmin = 0.007RB, rmax = 6RB for Models D and E. In Models A–C, the computational domain is located relatively outward because the size of the ionized region tends to larger than the Bondi radius. We set logarithmically spaced grids in the radial direction to simulate the flow within the Bondi radius in high resolution. We set uniformly spaced grids in the polar direction. The number of the grid points is set to (Nr, Nθ, Nϕ) = (100, 120, 1). Note that our simulations do not consider accretion flows within rmin. Instead, we assume properties of radiation emitted from the central region (see Section 2.1) and discuss gas accretion outside from the Bondi radii. As our initial conditions, we set a neutral uniform and static ($${\boldsymbol v}= 0$$) gas cloud with the density n∞ = 105 cm−3 and temperature T∞ = 104 K. The BH masses are assumed to be constant throughout our simulations. We impose the absorption inner boundary conditions for the gas density, gas pressure, and velocity to be damped smoothly (e.g. Kato, Mineshige & Shibata 2004), and the free outer boundary conditions for three components of the velocity and the specific entropy. We also fix the same gas density at r = rmax as the initial value for grids with an inflow velocity i.e. vr(r = rmax) < 0, otherwise the free boundary conditions are imposed for the density. The reflection symmetry with respect to the polar axis is imposed for non-radial components of the velocity. 3 RESULTS 3.1 Effects of anisotropic radiation Fig. 1 presents the time evolution of gas accretion rates on to a BH with MBH = 103 M⊙ for Model A (N = 0), B (N = 2), and C (N = 4). For the isotropic radiation, the accretion occurs episodically and the time-averaged rate ($${\sim } 0.9 \dot{M}_{\rm Edd}$$) is below the Eddington rate. For the anisotropic radiation, the accretion is less episodic and their rates tend to increase with time because of continuous accretion of neutral gas through the equatorial plane as shown later. The time-averaged accretion rate becomes higher with N (i.e. more anisotropic radiation). In Model C, the accretion rate is as high as $${\sim } 44 \dot{M}_{\rm Edd}$$ at the end of the simulation, where the radiation luminosity is higher than the Eddington value, namely ∼3.5LEdd. We continue the simulation until the ionization front reaches the outer boundary at t ≃ 1.2 × 105 yr, which is much longer than $$t_{\rm dyn} \equiv {\pi} \sqrt{R_{\rm B}^3/(8GM_{\rm BH})} \simeq 8.4 \times 10^3 {M_{\rm BH,3}} {\rm yr}$$. The accretion rate and the radiative luminosity at the end of each simulation are summarized in Table 1. Figure 1. View largeDownload slide Time evolution of gas accretion rates on to a 103 M⊙ BH in cases with isotropic radiation (Model A) and anisotropic radiation, Frad ∝ cos 2θ (Model B) and Frad ∝ cos 4θ (Model C). In Model A, the accretion occurs episodically and the time-averaged rate is below $$\dot{M}_{\rm Edd}$$. The accretion rate in Model C exceeds the Eddington rate ($$\dot{M}>10\dot{M}_{\rm Edd}$$). Figure 1. View largeDownload slide Time evolution of gas accretion rates on to a 103 M⊙ BH in cases with isotropic radiation (Model A) and anisotropic radiation, Frad ∝ cos 2θ (Model B) and Frad ∝ cos 4θ (Model C). In Model A, the accretion occurs episodically and the time-averaged rate is below $$\dot{M}_{\rm Edd}$$. The accretion rate in Model C exceeds the Eddington rate ($$\dot{M}>10\dot{M}_{\rm Edd}$$). In Fig. 2, we show two-dimensional distribution of the gas density and temperature at t = 1.1 × 105 yr in Models A and C, respectively. In the isotropic case (Model A), radiation from the central region heats up and ionizes the surrounding gas equally in all directions. When the ionized region expands and the temperature outside the Bondi radius increases, gas supply from large radii is suppressed since $$\dot{M}_{\rm B} \propto T_\infty ^{-3/2}$$. In this quiescent phase, thermal gas pressure is initially balanced inside and outside the ionized region. However, since the ionized gas within a new sonic point can accrete on to the centre, the gas pressure in the ionized region is gradually reduced due to gas depletion. As a result, the gas at the ionization front is accelerated inward by the positive pressure gradient (∂p/∂r > 0), leading to episodic accretion (see Fig. 1) as shown in previous studies (e.g. Ciotti & Ostriker 2001; Park & Ricotti 2011, 2012, IHO16). Note that radiation heating is a dominant process for suppressing the gas accretion, since the radiation luminosity is below the Eddington value (i.e. L < LEdd). Figure 2. View largeDownload slide Two-dimensional distribution of gas density (left-hand panel) and temperature (right-hand panel) for Model A (upper) and Model C (bottom), respectively, at t = 1.06 × 105 yr (≈13tdyn). In Model C (anisotropic radiation), the ionized region does not exist around the equatorial plane, from which a higher accretion rate is allowed. Figure 2. View largeDownload slide Two-dimensional distribution of gas density (left-hand panel) and temperature (right-hand panel) for Model A (upper) and Model C (bottom), respectively, at t = 1.06 × 105 yr (≈13tdyn). In Model C (anisotropic radiation), the ionized region does not exist around the equatorial plane, from which a higher accretion rate is allowed. In the anisotropic case (Model C), the shape of the ionized region is no longer isotropic. Although the ionization front expands towards the polar directions, the radiation flux in the equatorial direction (θ = π/2) is not intense enough to form an ionized region. Since the radiation luminosity viewed from the polar directions is significantly higher than the Eddington luminosity, strong outflows are launched by the radiation force on to electrons. The outflows collide with the ambient gas at r ≳ RB and form strong shocks with a high temperature of ≳ 5 × 104 K. The gas near the equator remains neutral and the temperature is kept at T ≃ 8000 K due to efficient atomic hydrogen cooling. Therefore, the gas can accrete on to the central BH through the neutral region increasing the density, unimpeded by radiation feedback. The radial profiles at θ = π/2 of the density and inflow velocity approach those for a Bondi solution with T ∼ 8000 K (ρ ∝ r−3/2 and vr ∝ r−1/2 at r ≲ RB). The accretion rate through the equatorial plane is estimated as $$\dot{M}_{\rm HI} \simeq \dot{M}_{\rm B}(\Omega _{\rm HI} / 4{\pi} )$$, where ΩHI is the solid angle of the neutral region. Since a half angle of the neutral region is ΘHI ≡ π/2 − θ ≈ 5° and $$\dot{M}_{\rm B} \approx 670 \dot{M}_{\rm Edd}$$, the inflow rate of the neutral gas is estimated as $$\dot{M}_{\rm HI}\simeq 58 \dot{M}_{\rm Edd}$$. This value agrees to the simulation result of $${\simeq } 44 \dot{M}_{\rm Edd}$$. Unlike our set-ups, Sugimura et al. (2016) assumed a shadow region with a half angle of ΘShad = 10°–45°, within which the radiation flux is completely shut off and accretion of neutral gas is allowed. In this case, they also found that gas accretion rates are given by $$\dot{M} \simeq \dot{M}_{\rm B}(\Omega _{\rm Shad} / 4{\pi} )$$, which is consistent with our results. 3.2 Transition to wholly neutral accretion Next, we examine a case with the highest BH mass of 5 × 105 M⊙ (Model E) but with the same radiation anisotropy (N = 4) as in Model C. Fig. 3 shows two-dimensional distribution of the gas density (left-hand panels) and temperature (right-hand panels) in the early stage at t ≃ 0.14tdyn (top panels) and the late stage at t ≃ 0.49tdyn (bottom panels), respectively. The time evolution of the accretion flow in Model E consists of the following three phases. Initially, the ionized region expands to r ∼ RB preferentially towards the polar directions, and gas accretion is allowed through the neutral region (top panels in Fig. 3). Unlike in Model C with lower MBH (i.e. lower $$\dot{m}$$), the opening angle of the ionized region in Model E gradually shrinks with time, allowing penetration of neutral gas towards the centre through the equator. The dense inflow totally covers the central region, where ionizing photons from the sink cell are completely absorbed by neutral hydrogen. As a result, the radiation flux in all directions (even in the polar directions) is blocked, and the ionized region is confined within r ∼ 0.01RB due to radiative recombination (bottom panels in Fig. 3). The radiation momentum absorbed by neutral gas produces warm outflows towards the bipolar directions. Figure 3. View largeDownload slide Two-dimensional distribution of gas density (left-hand panel) and temperature (right-hand panel) in Model E (MBH = 5 × 105 M⊙ and N = 4) in the early stage at t = 5.9 × 105 yr (≃0.14tdyn) (top), and in the late stage (after the transition) at t = 2.1 × 106 yr (≃0.49tdyn) (bottom), respectively. Figure 3. View largeDownload slide Two-dimensional distribution of gas density (left-hand panel) and temperature (right-hand panel) in Model E (MBH = 5 × 105 M⊙ and N = 4) in the early stage at t = 5.9 × 105 yr (≃0.14tdyn) (top), and in the late stage (after the transition) at t = 2.1 × 106 yr (≃0.49tdyn) (bottom), respectively. In order to understand how and why this transition of the accretion episode occurs, we examine structure of the gas flow in the central region (rmin ≤ r ≲ 0.04RB). Fig. 4 shows the density distribution and the velocity vectors of the flows. Each panel shows the time sequence from the left- (before the transition) to the right-hand side (after the transition). Fig. 5 also presents profiles of the electron fraction, the ram pressure ($$\rho v_{\theta }^2$$) and the gas pressure (pgas) as functions of the polar angle at r = 0.009RB. Before the transition shown in Fig. 4 (panel i1), geometrically thin and dense shock layers form along directions of θ ≈ 55° and 70°, located near the boundary between ionized and neutral gas (see Fig. 5). The inflow from the neutral region crosses the shock with a supersonic velocity, namely |vθ| > cs, and is decelerated (cs is the local sound speed). The shocked gas is ionized and accelerated outward in radial direction by strong radiation force. As shown in Fig. 5, ram pressure of the neutral gas inflow is comparable to gas pressure of the ionized region before the transition. As a result, the neutral gas inflow pushes the shock layer from the equator to the polars. Panel (i2) in Fig. 4 shows that dense clumps form at the shock front by compression caused by continuous accretion of neutral gas and strong radiation force from the centre. Although some clumps are ejected by radiation force, such clump formation occurs more frequently in the late stage. Eventually, the central radiating region is totally covered by dense clumps of neutral gas. Figure 4. View largeDownload slide Two-dimensional distribution of gas density and velocity structure in the central region (rmin ≤ r ≲ 0.04RB) for Model E (MBH = 5 × 105 M⊙). Each panel corresponds to the time at (i1) t = 1.7 × 105 yr (≃0.04tdyn), (i2) t = 5.9 × 105 yr (≃0.14tdyn), and (i3) t = 2.1 × 106 yr (≃0.49tdyn). Figure 4. View largeDownload slide Two-dimensional distribution of gas density and velocity structure in the central region (rmin ≤ r ≲ 0.04RB) for Model E (MBH = 5 × 105 M⊙). Each panel corresponds to the time at (i1) t = 1.7 × 105 yr (≃0.04tdyn), (i2) t = 5.9 × 105 yr (≃0.14tdyn), and (i3) t = 2.1 × 106 yr (≃0.49tdyn). Figure 5. View largeDownload slide Profiles of ram pressure $$\rho v_{\theta }^2$$ (blue), gas pressure pgas (red), and the electron fraction (green) in the Northern hemisphere from the polar (θ = 0°) to the equator (θ = 90°) at r = 0.009RB in Model E. The elapsed time is the same as that of panel (i1) in Fig. 4. Figure 5. View largeDownload slide Profiles of ram pressure $$\rho v_{\theta }^2$$ (blue), gas pressure pgas (red), and the electron fraction (green) in the Northern hemisphere from the polar (θ = 0°) to the equator (θ = 90°) at r = 0.009RB in Model E. The elapsed time is the same as that of panel (i1) in Fig. 4. Figs 6–8 present radial profiles of the gas density, temperature, and radial velocity along three different directions: θ = 0° (polar), 16°, and 90° (equator) in the early stage at t = 0.14tdyn (blue) and in the late stage at t = 0.49tdyn (red). The black dashed curve represents the isothermal Bondi profile with T = 104K. In the early stage before the transition, radiation force and heating affect properties of the accreting material. Inside the ionized region, the gas is depleted and a density cavity is created (θ = 0° and 16°). Colliding outflows due to the radiation force on electrons (dashed curves in Fig. 8) produce strong shocks with a high temperature of ∼106 K. On the other hand, neutral gas accretes through the equator, keeping the temperature at T ≃ 8000 K, as shown in Section 3.1 (Model C). After the transition, the ionized gas becomes neutral due to radiative recombination and the density cavity is filled. The gas temperature settles down to ∼104 K in all directions because of radiative cooling by atomic hydrogen. Along the direction of θ = 16°, the ionized region still exists but its size is too small (∼0.01RB) to affect gas dynamics at the Bondi radius. As a result of absorption of anisotropic radiation by neutral hydrogen, the radiation momentum produces warm outflows with T ≃ 8000 K towards narrow conical regions along the polars (θ ∼ 0°–20°). Figure 6. View largeDownload slide Radial profiles of the gas density in Model E along three different directions of θ = 0° (polar), θ = 16°, and θ = 90° (equator) in the early stage at t = 0.14tdyn (blue) and in the late stage at t = 0.49tdyn (red). The black curve represents an isothermal Bondi profile with T = 104K. Figure 6. View largeDownload slide Radial profiles of the gas density in Model E along three different directions of θ = 0° (polar), θ = 16°, and θ = 90° (equator) in the early stage at t = 0.14tdyn (blue) and in the late stage at t = 0.49tdyn (red). The black curve represents an isothermal Bondi profile with T = 104K. Figure 7. View largeDownload slide Same as Fig. 6 but for the gas temperature. Figure 7. View largeDownload slide Same as Fig. 6 but for the gas temperature. Figure 8. View largeDownload slide Same as Fig. 6 but for the radial velocity. Solid and dashed curves present the inflow and outflow velocity, respectively. Figure 8. View largeDownload slide Same as Fig. 6 but for the radial velocity. Solid and dashed curves present the inflow and outflow velocity, respectively. Fig. 9 shows the time evolution of the inflow (green long-dashed) and outflow (blue short-dashed) rate measured at the Bondi radius, and the BH feeding rate (red solid). The normalized BH feeding rate of $$\dot{M}/\dot{M}_{\rm Edd}$$ continuously increases with the inflow rate from the Bondi radius. An open circle indicates the epoch of the transition to the wholly neutral accretion phase. After the transition, the feeding rate becomes highly episodic because strong radiation force and high ram pressure are tightly balanced at the innermost region. At the end of the simulation, the central BH is fed through the equator at the time-averaged rate of ∼5 × 104LEdd/c2, which corresponds to ∼ 50 per cent of the inflow rate from the Bondi radius. On the other hand, the accreting material is accelerated outward into the bipolar directions (θ ≲ 20°) due to absorption of the radiation momentum with ≫LEdd/c. The outflow rate ends up $${\simeq } 5\times 10^3 \dot{M}_{\rm Edd}$$, which is ∼ 10 per cent of the BH feeding rate. Since the Mach number of the outflowing matter is as high as four at the Bondi radius (see Fig. 8), the typical velocity is estimated as ≲ 40 km s − 1. Sugimura et al. (2016) reported that hot and ionized outflows with T ≃ 7 × 104 K are launched from a rapid accreting system exposed to strong anisotropic radiation, which is consistent with our cases where the transition does not occur as in Model C. On the other hand, our simulation in Model E suggests that warm and neutral outflows with T ≃ 8000K are likely to be produced because of the transition to the wholly neutral phase. Figure 9. View largeDownload slide Time evolution of inflow (green long-dashed) and outflow (blue short-dashed) rates at the Bondi radius and the BH feeding rate (red solid) in Model E (MBH = 5 × 105 M⊙). Figure 9. View largeDownload slide Time evolution of inflow (green long-dashed) and outflow (blue short-dashed) rates at the Bondi radius and the BH feeding rate (red solid) in Model E (MBH = 5 × 105 M⊙). In Fig. 10, we show the BH feeding rates for two cases with MBH = 5 × 104 M⊙ (blue, Model D) and 5 × 105 M⊙ (red, Model E) for comparison. The normalized accretion rate in Model D is saturated at a constant value unlike in Model E. In fact, we do not find a transition to a wholly neutral phase within the entire simulation time. Even without the transition, neutral gas inflows have non-radial motions and dense clumps are formed at the boundary between the ionized and neutral region as in Model E. However, ram pressure of the infalling clumps is too weak to cover the central radiation source. Therefore, this result implies the critical BH mass for the transition is located around MBH ∼ 5 × 105 M⊙ for the initial conditions we adopted (n∞ = 105 cm−3 and T∞ = 104 K). Figure 10. View largeDownload slide Time evolution of the accretion rates for two cases with MBH = 5 × 104 M⊙ (Model D) and MBH = 5 × 105 M⊙ (Model E). An open circle presents the epoch where the ionization region collapses by radiative recombination. The transition occurs for MBH ≳ 5 × 105 M⊙. Figure 10. View largeDownload slide Time evolution of the accretion rates for two cases with MBH = 5 × 104 M⊙ (Model D) and MBH = 5 × 105 M⊙ (Model E). An open circle presents the epoch where the ionization region collapses by radiative recombination. The transition occurs for MBH ≳ 5 × 105 M⊙. Transition to the wholly neutral phase is more likely to occur in more massive BH cases. The physical reason is that the size of the ionized region relative to the Bondi radius becomes smaller with the BH mass, following $$R_{\rm HII}/R_{\rm B}\propto M_{\rm BH}^{-3/2}$$, as found by IHO16. When the radiation flux is anisotropic as in our two-dimensional simulations, neutral gas accretes throughout the equatorial region with a significant fraction of solid angle (e.g. ΩHI/4π ≃ 0.7 for Model E) even if we do not assume shadow structures artificially. The size of the neutral region becomes larger with the BH mass, and eventually the transition occurs for the highest BH mass (Model E). 4 SUMMARY AND DISCUSSION In this paper, we study rapid growth of BHs with a mass range of 103 ≤ MBH ≤ 5 × 105 M⊙ via gas accretion exposed to anisotropic super-Eddington radiation. Radiation flux from an accretion disc around the nuclear BH is likely to be collimated towards the bipolar directions. To address the effects of the anisotropic radiation on the accretion flow, we perform two-dimensional radiation hydrodynamical simulations, including multifrequency radiation transfer and non-equilibrium primordial chemistry. We here consider a BH embedded in an initially uniform gas cloud with the density of n∞ = 105 cm−3 and temperature of T∞ = 104 K. We have found that hyper-Eddington accretion from the Bondi radius can occur when the radiation is more anisotropic and the BH mass is higher. Moreover, when the BH mass is as large as MBH ≳ 5 × 105 M⊙, ionized regions around the BH collapse due to super-sonic non-radial motions of neutral gas inflows and efficient radiative recombination. Due to absorption of the radiation momentum, the neutral gas is accelerated outward to narrow conical regions along the polars, producing warm outflows with T ≃ 8000 K. We adopt a simple model for anisotropic radiation from the central accretion disc, namely, Frad(θ) ∝ cos Nθ. As a result of the anisotropic radiation, rapid accretion is allowed from the equatorial plane. This assumption would be reasonable at the vicinity of the central BH because a super-Eddington accretion disc is optically and geometrically thick, so that the radiation flux is mildly collimated to the polars (e.g. Ohsuga et al. 2005; Jiang et al. 2014; Sa̧dowski et al. 2015). However, although the radiation propagates outward, such anisotropic radiation fields could be smoothed out because of radiative diffusion. In order to explore the anisotropy of the emergent radiation, therefore radiation-hydrodynamical simulations with a wide dynamical range are required. Recently, Sugimura et al. (2016) have performed two-dimensional simulations of accreting BHs under anisotropic radiation and found that anisotropic radiation allows gas accretion through the equatorial plane. Although our calculations have been pursued independently from their work, our simulation set-ups are similar, mainly except two things. First, they assumed the existence of a shadow region caused by attenuation due to a geometrically thick disc and/or disc winds well inside their computation domains. Secondly, our simulations resolve even smaller scales, namely ∼ 1 per cent of the Bondi radius. We can stress that our set-up is more preferable to follow non-radial motions of dense and neutral inflows, which lead to a transition to a wholly neutral accretion phase. Conversely, Sugimura et al. (2016) did not find such transitions, probably because (1) the BH mass adopted in their simulations is lower than that in our Model E, and (2) the radius of the inner most grid in their simulations is larger than that in our simulations. We discuss the subsequent evolution of the accretion flows in Model E (MBH = 105.7 M⊙). The gas profile approaches an isothermal Bondi solution with T ∼ 104 K, where $$\dot{M}_{\rm B}/\dot{M}_{\rm Edd}\simeq 3 \times 10^5 n_{\infty .5}M_{\rm BH,5.7}$$ if the ram pressure of the neutral inflow dominates the radiation force. This condition is written as $$L/L_{\rm Edd} \lesssim 139 M_{\rm BH,5.7}^{3/2}n_{\infty ,5} (r_{\star }/7\times 10^{17}\,{\rm cm})^{-1/2}$$ (Sakurai et al. 2016), where r⋆ is the location of the photosphere, from which the radiation emerges actually. Note that the photosphere is located inside the innermost grid (r⋆ < rmin = 7 × 1017 cm). Assuming that the BH is fed at the Bondi rate, the radiation luminosity towards the polar directions is estimated as L/LEdd ≈ 106 (see equation 12). Therefore, we expect that the overall behaviours of the accretion flows in Model E would be maintained. We set the ambient gas density to n∞ = 105 cm−3. If PopIII stars emit ionizing photons, the gas inside the host dark matter halo with ∼105–6 M⊙ (the so-called mini-halo, z ∼ 20) is evaporated by the irradiation, reducing the density much below that we assume (e.g. Kitayama et al. 2004; Johnson & Bromm 2007). However, when the halo grows to be as massive as ≳ 107 M⊙ with a virial temperature of ≳ 104 K, (the so-called atomic-cooling halo, z ∼ 10), the gas re-collapses on to the halo, and a density profile with a core and an envelope following ρ ∝ r−2 is built up due to Lyman α cooling. Alternatively, when remnant BHs are located at the centre of an atomic-cooling halo, the gas density would be high enough to realize the conditions we assume. The massive BHs we consider in Models D and E would be formed by, for example, direct-collapse of supermassive stars in atomic-cooling haloes. Since a supermassive star has a bloated cold envelope with ≃5000 K, its radiation feedback would be too weak to photoevaporate a surrounding medium (e.g. Chon, Hosokawa & Yoshida 2018). The supermassive star finally collapses to an equally massive BH (∼104–5 M⊙) without mass-loss and explosions. Thus, the gas density profile (ρ ∝ r−2) in the halo would not be broken before the formation of the BH seed (see also IHO16). As shown in Fig. 4, non-radial motions of neutral inflows near the inner boundary are essential for leading to the transition to a wholly neutral phase. However, our simulations neglect non-radial components in radiation flux $$F_{\rm rad}^{\theta }$$, which can be produced by radiative recombination in the ionized region.1 The radiation flux due to recombination is estimated as $$F_{\rm rad}^{\theta } \sim 4{\pi} \eta r_{\rm min}$$, where $$\eta \sim \langle h\nu _{\rm rec} \rangle \alpha _{\rm B}n_{\rm e}^2 \sim 5\times 10^{-17}(n_{\rm e}/10^4\, {\rm cm}^{-3})^2\,{\rm erg\,cm^{-3}\,s^{-1}}$$ is the emissivity, 〈hνrec〉 is the average energy of recombination photons, αB is the case B recombination rate coefficient, and ne is the number density of electrons. Thus, the non-radial radiation flux is estimated as $$|F_{\rm rad}^{\theta }/F_{\rm rad}^r| \sim 10^{-7}$$ at the innermost region, and the radiation force in non-radial directions is negligible. Throughout this paper, we neglect formation of molecular hydrogen (H2), which is the main coolant in primordial gas below ∼104 K. We briefly mention that H2 cooling is unlikely to affect the properties of warm/neutral outflows. In a warm and dense primordial gas cloud (n ≳ 105 cm − 3 and T ∼ 8000 K), the H2 fraction relative to atoms is determined by the balance between formation through electron-catalysed reactions and collisional dissociation with gas particles (Omukai 2001). The resultant fraction is estimated as ∼10−8 as long as the temperature is kept at >4000 K, below which the collisional dissociation becomes inefficient (Inayoshi et al. 2014). As shown in Fig. 7, the temperature cools down to ≃1000–2000 K due to adiabatic expansion at r > 2 × 1019 cm in the outflowing matter (θ = 0°). In this lower temperature region, far-ultraviolet radiation in the Lyman-Werner (LW) bands (hνLW = 10.2–13.6 eV) from the central region plays an important role to determine the H2 fraction. We can estimate the LW flux in the optically thin limit as   \begin{eqnarray} F_{\rm LW,thin}\, = \, f_{\rm LW}(N+1)\frac{2\left[1+\ln (\dot{m}/20)\right]\,L_{\rm Edd}}{4{\pi} r^2\langle \nu _{\rm LW}\rangle }\nonumber \\ \simeq \, 1.2\times 10^{-11} \,{\rm erg \,s^{-1}\, cm^{-2} \,Hz^{-1}} \nonumber\\ \times \left(\frac{f_{\rm LW}}{0.3}\right) \left(\frac{N+1}{5}\right) \!\left(\frac{M_{\rm BH}}{5\times 10^5 \,{\rm M}_{\odot }}\right)^{-1} \!\left(\frac{r}{R_{\rm B}}\right)^{-2}, \,\,\,\nonumber\\ \end{eqnarray} (13)where fLW is the ratio of the LW flux to the bolometric value and 〈νLW〉 = 12eV/h( = 2.9 × 1015Hz). Since the H2 column density is $$N_{\rm H_2}\sim 10^{17} {\rm cm}^{-2}(n_\infty /10^5\, {\rm cm}^{-3})(r/R_{\rm B})$$ at most, the emergent LW flux can be reduced by the H2 self-shielding effect. The shielding factor is simply estimated as $$f_{\rm shield}\simeq (N_{\rm H_2}/10^{14}\,{\rm cm}^{-2})^{-3/4}$$ for a higher column density (Draine & Bertoldi 1996, see also Wolcott-Green, Haiman & Bryan 2011), and thus FLW = fshieldFLW, thin ≃ 7 × 10−14 erg s−1 cm−2 Hz−1, which is higher than the critical LW flux,2FLW, crit ∼ 10−19–10−16 erg s−1 cm−2 Hz−1 (e.g. Sugimura, Omukai & Inoue 2014; Inayoshi & Tanaka 2015; Latif, Schleicher & Hartwig 2016; Wolcott-Green, Haiman & Bryan 2017) required to dissociate H2 and suppress H2 cooling in the warm and neutral gas even at r ∼ 1020 cm. We assume the Ly α cooling (2p → 1s state) in the optically thin limit. In reality, Ly α photons are frequently scattered in neutral gas before escaping from the accretion system. In such a cloud where Ly α photons are trapped, two-photon emission (2s → 1s state) and H− free–bound emission (H + e− → H− + γ) dominate as radiative cooling processes (e.g. Omukai 2001; Schleicher, Spaans & Glover 2010). As shown by the previous studies, thermal evolution of a collapsing gas does not change significantly from that assuming the optically thin Ly α cooling (see also Spaans & Silk 2006; Ge & Wise 2017). Furthermore, additional cooling processes, e.g. fine-structure lines of C ii and O i, and thermal emission by dust grains would be important in a metal-enriched gas cloud. The metallicity effects on hyper-Eddington accretion are left for future investigations. Acknowledgements We thank Zoltán Haiman and Kazuyuki Sugimura for useful discussions. This work is partially supported by the Simons Foundation through the Simons Society of Fellows (KI) and by Japan Society for the Promotion of Science Grant-in-Aid for Scientific Research (A) (17H01102 KO), for Scientific Research (C) (15K05036 KO; 16K05309 KO and 17K05383 SM)) and for Young Scientists (17K14260 HRT). This research was also supported by the Ministry of Education, Culture, Sports, Science and Technology as ‘Priority Issue on Post-K computer’ (Elucidation of the Fundamental Laws and Evolution of the Universe) and the Joint Institute for Computational Fundamental Science. Numerical computations were carried out on Cray XC30 at Center for Computational Astrophysics, National Astronomical Observatory of Japan. Footnotes 1 The non-radial components in radiation flux is unlikely to be produced by electron scattering because the ionized gas is optically thin. In fact, before the transition in Model E, the optical depth is as low as 0.16 at r < RHII along a line of θ = 10°, which is inside the ionized region. 2 The critical LW flux depends on the radiation spectrum, namely on the ratio of the flux at 2 eV to the LW flux (e.g. Sugimura, Omukai & Inoue 2014). This is because photons with ∼2 eV dissociate H− ions, which is an intermediate product in the H2 formation channels (H + e− → H− + γ; H− + H → H2 + e−) and contribute to H2 dissociation indirectly (Omukai 2001). REFERENCES Abel T., Anninos P., Zhang Y., Norman M. L., 1997, New Astron. , 2, 181 https://doi.org/10.1016/S1384-1076(97)00010-9 CrossRef Search ADS   Alexander T., Natarajan P., 2014, Science , 345, 1330 https://doi.org/10.1126/science.1251053 CrossRef Search ADS PubMed  Alvarez M. A., Wise J. H., Abel T., 2009, ApJ , 701, L133 CrossRef Search ADS   Anninos P., Zhang Y., Abel T., Norman M. L., 1997, New Astron. , 2, 209 https://doi.org/10.1016/S1384-1076(97)00009-2 CrossRef Search ADS   Becerra F., Greif T. H., Springel V., Hernquist L. E., 2015, MNRAS , 446, 2380 CrossRef Search ADS   Begelman M. C., 1979, MNRAS , 187, 237 CrossRef Search ADS   Bondi H., 1952, MNRAS , 112, 195 CrossRef Search ADS   Chon S., Hirano S., Hosokawa T., Yoshida N., 2016, ApJ , 832, 134 CrossRef Search ADS   Chon S., Hosokawa T., Yoshida N., 2018, MNRAS , in press Ciotti L., Ostriker J. P., 2001, ApJ , 551, 131 CrossRef Search ADS   Ciotti L., Ostriker J. P., Proga D., 2009, ApJ , 699, 89 CrossRef Search ADS   Devecchi B., Volonteri M., 2009, ApJ , 694, 302 CrossRef Search ADS   Draine B. T., Bertoldi F., 1996, ApJ , 468, 269 CrossRef Search ADS   Fan X. et al.  , 2004, AJ , 128, 515 CrossRef Search ADS   Ge Q., Wise J., 2017, MNRAS , 472, 2773 Glover S. C. O., Jappsen A.-K., 2007, ApJ , 666, 1 CrossRef Search ADS   Haiman Z., Loeb A., 2001, ApJ , 552, 459 CrossRef Search ADS   Harten A., Lax P. D., van Leer B., 1983, SIAM Rev. , 25, 35 https://doi.org/10.1137/1025002 CrossRef Search ADS   Inayoshi K., Tanaka T. L., 2015, MNRAS , 450, 4350 CrossRef Search ADS   Inayoshi K., Omukai K., Tasker E., 2014, MNRAS , 445, L109 CrossRef Search ADS   Inayoshi K., Haiman Z., Ostriker J. P., 2016, MNRAS , 459, 3738 (IHO16) CrossRef Search ADS   Jiang Y.-F., Stone J. M., Davis S. W., 2014, ApJ , 796, 106 CrossRef Search ADS   Johnson J. L., Bromm V., 2007, MNRAS , 374, 1557 CrossRef Search ADS   Kato Y., Mineshige S., Shibata K., 2004, ApJ , 605, 307 CrossRef Search ADS   Katz H., Sijacki D., Haehnelt M. G., 2015, MNRAS , 451, 2352 CrossRef Search ADS   King A., 2003, ApJ , 596, L27 CrossRef Search ADS   Kitayama T., Yoshida N., Susa H., Umemura M., 2004, ApJ , 613, 631 CrossRef Search ADS   Kormendy J., Ho L. C., 2013, ARA&A , 51, 511 CrossRef Search ADS   Latif M. A., Schleicher D. R. G., Hartwig T., 2016, MNRAS , 458, 233 CrossRef Search ADS   Lupi A., Haardt F., Dotti M., Fiacconi D., Mayer L., Madau P., 2016, MNRAS , 456, 2993 CrossRef Search ADS   Madau P., Rees M. J., 2001, ApJ , 551, L27 CrossRef Search ADS   Madau P., Haardt F., Dotti M., 2014, ApJ , 784, L38 CrossRef Search ADS   Milosavljević M., Bromm V., Couch S. M., Oh S. P., 2009, ApJ , 698, 766 CrossRef Search ADS   Mortlock D. J. et al.  , 2011, Nature , 474, 616 CrossRef Search ADS PubMed  Murray N., Quataert E., Thompson T. A., 2005, ApJ , 618, 569 CrossRef Search ADS   Ohsuga K., Mineshige S., 2007, ApJ , 670, 1283 CrossRef Search ADS   Ohsuga K., Mineshige S., 2011, ApJ , 736, 2 CrossRef Search ADS   Ohsuga K., Mori M., Nakamoto T., Mineshige S., 2005, ApJ , 628, 368 CrossRef Search ADS   Ohsuga K., Mineshige S., Mori M., Kato Y., 2009, PASJ , 61, L7 CrossRef Search ADS   Omukai K., 2001, ApJ , 546, 635 CrossRef Search ADS   Pacucci F., Ferrara A., 2015, MNRAS , 448, 104 CrossRef Search ADS   Park K., Ricotti M., 2011, ApJ , 739, 2 CrossRef Search ADS   Park K., Ricotti M., 2012, ApJ , 747, 9 CrossRef Search ADS   Pezzulli E., Valiante R., Schneider R., 2016, MNRAS , 458, 3047 CrossRef Search ADS   Regan J. A., Johansson P. H., Haehnelt M. G., 2014, MNRAS , 439, 1160 CrossRef Search ADS   Ryu T., Tanaka T. L., Perna R., Haiman Z., 2016, MNRAS , 460, 4122 CrossRef Search ADS   Sakurai Y., Inayoshi K., Haiman Z., 2016, MNRAS , 461, 4496 CrossRef Search ADS   Sakurai Y., Yoshida N., Fujii M. S., Hirano S., 2017, MNRAS , 472, 1677 Schleicher D. R. G., Spaans M., Glover S. C. O., 2010, ApJ , 712, L69 CrossRef Search ADS   Shakura N. I., Sunyaev R. A., 1973, A&A , 24, 337 Silk J., Rees M. J., 1998, A&A , 331, L1 Soltan A., 1982, MNRAS , 200, 115 CrossRef Search ADS   Spaans M., Silk J., 2006, ApJ , 652, 902 CrossRef Search ADS   Stone N. C., Küpper A. H. W., Ostriker J. P., 2017, MNRAS , 467, 4180 Sugimura K., Omukai K., Inoue A. K., 2014, MNRAS , 445, 544 CrossRef Search ADS   Sugimura K., Hosokawa T., Yajima H., Omukai K., 2016, MNRAS , 469, 62 CrossRef Search ADS   Sa̧dowski A., Narayan R., Tchekhovskoy A., Abarca D., Zhu Y., McKinney J. C., 2015, MNRAS , 447, 49 CrossRef Search ADS   Takahashi H. R., Ohsuga K., 2013, ApJ , 772, 127 CrossRef Search ADS   Takahashi H. R., Ohsuga K., Kawashima T., Sekiguchi Y., 2016, ApJ , 826, 23 CrossRef Search ADS   Tanaka T., Haiman Z., 2009, ApJ , 696, 1798 CrossRef Search ADS   Valiante R., Schneider R., Volonteri M., Omukai K., 2016, MNRAS , 457, 3356 CrossRef Search ADS   van Leer B., 1977, J. Comput. Phys. , 23, 263 https://doi.org/10.1016/0021-9991(77)90094-8 CrossRef Search ADS   Visbal E., Haiman Z., Bryan G. L., 2014, MNRAS , 445, 1056 CrossRef Search ADS   Volonteri M., Rees M. J., 2005, ApJ , 633, 624 CrossRef Search ADS   Volonteri M., Haardt F., Madau P., 2003, ApJ , 582, 559 CrossRef Search ADS   Watarai K.-y., Fukue J., Takeuchi M., Mineshige S., 2000, PASJ , 52, 133 CrossRef Search ADS   Whalen D., Norman M. L., 2006, ApJS , 162, 281 CrossRef Search ADS   Whalen D. J., Norman M. L., 2008, ApJ , 672, 287 CrossRef Search ADS   Wolcott-Green J., Haiman Z., Bryan G. L., 2011, MNRAS , 418, 838 CrossRef Search ADS   Wolcott-Green J., Haiman Z., Bryan G. L., 2017, MNRAS , 469, 3329 Wu X.-B. et al.  , 2015, Nature , 518, 512 CrossRef Search ADS PubMed  Yajima H., Khochfar S., 2016, MNRAS , 457, 2423 CrossRef Search ADS   Yu Q., Tremaine S., 2002, MNRAS , 335, 965 CrossRef Search ADS   © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society

Journal

Monthly Notices of the Royal Astronomical SocietyOxford University Press

Published: May 1, 2018

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 12 million articles from more than
10,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Unlimited reading

Read as many articles as you need. Full articles with original layout, charts and figures. Read online, from anywhere.

Stay up to date

Keep up with your field with Personalized Recommendations and Follow Journals to get automatic updates.

Organize your research

It’s easy to organize your research with our built-in tools.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

Monthly Plan

  • Read unlimited articles
  • Personalized recommendations
  • No expiration
  • Print 20 pages per month
  • 20% off on PDF purchases
  • Organize your research
  • Get updates on your journals and topic searches

$49/month

Start Free Trial

14-day Free Trial

Best Deal — 39% off

Annual Plan

  • All the features of the Professional Plan, but for 39% off!
  • Billed annually
  • No expiration
  • For the normal price of 10 articles elsewhere, you get one full year of unlimited access to articles.

$588

$360/year

billed annually
Start Free Trial

14-day Free Trial