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

TL;DR

This study uses 2D radiation hydrodynamical simulations to explore how anisotropic super-Eddington radiation influences black hole growth and outflows, revealing conditions under which rapid accretion occurs and the formation of warm, outflowing gas.

Contribution

It demonstrates that anisotropic radiation allows less massive black holes to grow rapidly and produces warm outflows, contrasting with previous spherically symmetric models.

Findings

Less massive BHs can accrete at high rates due to anisotropic radiation.

Hot outflows with temperatures around 10^5 K are produced in bipolar directions.

Neutral, warm gas with ~8000 K covers the central region in massive BHs.

Abstract

We perform two-dimensional radiation hydrodynamical simulations of accretion flows onto a black hole (BH) with a mass of $10^3\leq M_{\rm BH}/M_{\odot} \lesssim 10^6$ in order to study rapid growth of BHs in the early Universe. For spherically symmetric flows, hyper-Eddington accretion onto the BH from outside the Bondi radius can occur unimpeded by radiation feedback only when the BH mass is higher than $\simeq 10^4~M_{\odot}(n_\infty/10^5~{\rm cm}^{-3})^{-1}(T_\infty/10^4~{\rm K})^{3/2}$, where $n_\infty$ and $T_\infty$ are the density and temperature of ambient gas. Here, we study the properties of accretion flows exposed to anisotropic radiation from a nuclear accretion disk with a luminosity higher than the Eddington value ($L_{\rm Edd}$) due to collimation toward the bipolar directions. We find that, unlike the spherically symmetric case, even less massive BHs with $M_{\rm BH} < 10^4~M_{\odot}$ can be fed by surrounding gas at high accretion rates of $\gtrsim L_{\rm Edd}/c^2$ through the equatorial plane, while ionized regions expand to the polar directions producing hot outflows with $T\sim 10^5$K. For more massive BHs with $M_{\rm BH}\gtrsim 5\times 10^5~M_{\odot}$, neutral gas through the equatorial plane totally covers 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 with $T \simeq 8000~{\rm K}$ . The central BH is fed through the equator at the averaged rate of $\sim 5\times 10^4~L_{\rm Edd}/c^2$, which corresponds to $\sim 50~\%$ of the inflow rate from the Bondi radius. Moreover, radiation momentum absorbed by neutral hydrogen produces warm outflows toward the bipolar directions at $\sim 30~\%$ of the BH feeding rate and with a typical velocity of $\simeq 50~{\rm km~s}^{-1}$.

Figures

Tables

Model	$M_{\mathrm{BH}}\left(M_{\odot }\right)$	$N$	$\dot{M}/{\dot{M}}_{\mathrm{Edd}}$	$L/L_{\mathrm{Edd}}$
A	${10}^3$	0	$\sim 0.9$	$\sim 0.09$
B	${10}^3$	2	$\approx 2.7$	$\approx 0.27$
C	${10}^3$	4	$\approx 44$	$\approx 3.5$
D	$5\times {10}^4$	4	$\sim 1\times {10}^3$	$\sim 10$
E	$5\times {10}^5$	4	$\sim 5\times {10}^4$	$\sim 18$

Equations

$$M_{\mathrm{BH}}\text{ga}10^4{\mathrm{M}}_{\odot }{\left(\frac{n_{\infty }}{10^5{\mathrm{cm}}^{-3}}\right)}^{-1}{\left(\frac{T_{\infty }}{10^4\mathrm{K}}\right)}^{3/2},$$

$$R_{\mathrm{B}}\equiv \frac{G{M}_{\mathrm{BH}}}{c_{\infty }^2}\simeq 1.97\times 10^{18}{M}_{\mathrm{BH},4}{T}_{\infty ,4}^{-1}\mathrm{cm},$$

$${\dot{M}}_{\mathrm{B}}\equiv \pi e^{3/2}{\rho }_{\infty }\frac{G^2{M}_{\mathrm{BH}}^2}{c_{\infty }^3},$$

$$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \mathbf{v})=0,$$

$$\frac{\partial \left(\rho v_r\right)}{\partial t}+\nabla \cdot (\rho v_r\mathbf{v})=-\frac{\partial p}{\partial r}+\rho \left(\frac{v_{\theta }^2}{r}+\frac{v_{\varphi }^2}{r}\right)-\rho \frac{\partial \psi }{\partial r}+f_{\mathrm{rad}},$$

$$\frac{\partial \left(\rho r{v}_{\theta }\right)}{\partial t}+\nabla \cdot (\rho r{v}_{\theta }\mathbf{v})=-\frac{\partial p}{\partial \theta }+\rho v_{\varphi }^2\cot \theta ,$$

$$\frac{\partial \left(\rho r{v}_{\varphi }\sin \theta \right)}{\partial t}+\nabla \cdot \left(\rho r{v}_{\varphi }\sin \theta \mathbf{v}\right)=0,$$

$$\frac{\partial e}{\partial t}+\nabla \cdot [(e+p)\mathbf{v}]=-\frac{G{M}_{\mathrm{BH}}\rho }{r^2}{v}_r-\Lambda +\Gamma ,$$

$$\frac{1}{r^2}\frac{\mathrm{d}}{\mathrm{d}r}(r^2{F}_{\nu })=-\rho {\kappa }_{\nu }c{E}_{\nu },$$

$$f_{\mathrm{rad}}=\frac{n{x}_{\mathrm{e}}}{c}{\int }_{{\nu }_{\min }}^{{\nu }_{\max }}{\sigma }_{\mathrm{es}}{F}_{\nu }\mathrm{d}\nu +\frac{{\Gamma }_{\mathrm{H}}}{c}.$$

$$\frac{\displaystyle L}{\displaystyle L_{\rm Edd}}=\left\{\begin{array}[]{ll}2\left[1+\ln{(\dot{m}/20)}\right]&{\rm for}\ \dot{m}\geq 20,\\ \dot{m}/10&{\rm for}\ \dot{m}<20,\end{array}\right.$$

$$F_{\nu }(r=r_{\min },\theta )=\frac{(N+1)L_{\nu }}{4\pi r_{\min }^2}{\cos }^N\theta ,$$

$$F_{\mathrm{LW},\mathrm{thin}}$$

$$\simeq 1.2\times 10^{-11}\mathrm{erg}{\mathrm{s}}^{-1}\mathrm{c}{\mathrm{m}}^{-2}\mathrm{H}{\mathrm{z}}^{-1}$$

$$\times \left(\frac{f_{\mathrm{LW}}}{0.3}\right)\left(\frac{N+1}{5}\right){\left(\frac{M_{\mathrm{BH}}}{5\times 10^5{\mathrm{M}}_{\odot }}\right)}^{-1}{\left(\frac{r}{R_{\mathrm{B}}}\right)}^{-2},$$

Full text

1 Introduction

Observations of bright quasars at high redshifts $(z\ga 6)$ have revealed that supermassive black holes (SMBHs) with $M_{\rm BH}\ga 10^{9}~{}{\rm M}_{\odot}$ are formed within $\la 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 $\sim 100~{}{\rm M}_{\odot}$ 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 disk because the radiation luminosity increases with the accretion rate onto the BH. If the radiative efficiency of the disk is $\eta\simeq 0.1$ (Shakura & Sunyaev, 1973; Soltan, 1982; Yu & Tremaine, 2002), the radiation luminosity exceeds the Eddington luminosity ($L_{\rm Edd}$) at the accretion rate of $\ga L_{\rm Edd}/(\eta c^{2})$, above which the radiation force dominates the gravity of the central BH. Assuming this Eddington-limited accretion rate, it takes $\ga 1$ Gyr for stellar-mass seed BHs to grow up to $\sim 10^{9}~{}{\rm M}_{\odot}$. This fact requires more rapid accretion mechanisms to form SMBHs by $z\sim 6$, when the corresponding cosmic age is $\sim 0.8$ Gyr.

As alternative possibilities, formation of massive seed BHs with $\ga 10^{3}-10^{5}~{}{\rm M}_{\odot}$ has been proposed: direct-collapse of supermassive stars in protogalaxies (e.g., Regan, Johansson & Haehnelt, 2014; Inayoshi, Omukai & Tasker, 2014; Becerra et al., 2015; Inayoshi & Tanaka, 2015; Visbal, Haiman & Bryan, 2014; Latif, Schleicher & Hartwig, 2016; Chon et al., 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 two in the case of the Eddington limited accretion, a high duty cycle ($\sim 100\%$) is still required to explain the existence of SMBHs with $\sim 10^{9}~{}{\rm M}_{\odot}$ 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; Valiante et al., 2016; Pezzulli, Valiante & Schneider, 2016). By means of radiation hydrodynamical simulations of accretion flows onto 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 onto 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; Milosavljević et al., 2009; Alvarez, Wise & Abel, 2009; Ciotti, Ostriker & Proga, 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 $\ga 5000~{}L_{\rm Edd}/c^{2}$, 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 which follows a Bondi profile with $8000~{}{\rm K}$. When the emergent radiation luminosity is limited to $\la L_{\rm Edd}$ 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,

[TABLE]

where $n_{\infty}$ and $T_{\infty}$ are density and temperature of the ambient gas, respectively. Hereafter, we define $M_{{\rm BH},x}=(M_{\rm BH}/10^{x}~{}{\rm M}_{\odot})$, $n_{\infty,5}=(n_{\infty}/10^{5}~{}{\rm cm}^{-3})$, and $T_{\infty,4}=(T_{\infty}/10^{4}~{}{\rm K})$. Note that Eq.(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 $>L_{\rm Edd}$ emerges from the center, the above conditions are not affected as long as the luminosity is as low as $(10-30)\times L_{\rm Edd}$ (Sakurai, Inayoshi & Haiman, 2016). To satisfy the condition of Eq. (1), gas-rich regions with $n\ga 10^{7}~{}{\rm cm}^{-3}$ are required for stellar-mass PopIII BHs with masses of $\la 100~{}{\rm M}_{\odot}$. During the assembly of a protogalaxy, some PopIII BHs would fall into a gaseous circumnuclear disk, where the gas accretion rate is likely to be $\ga 5000~{}L_{\rm Edd}/c^{2}$ (Ryu et al. 2016, see also Lupi et al. 2016).

Radiation flux from an accretion disk around the nuclear BH is more likely to be anisotropic, i.e., preferentially collimated around the rotation axis of the disk (Ohsuga et al., 2009; Ohsuga & Mineshige, 2011; Jiang, Stone & Davis, 2014; Sa̧dowski et al., 2015; Takahashi et al., 2016). On the other hand, the disk-like accretion could reduce the feedback effect because intense radiation is significantly attenuated by the dense accreting material in the disk (Ohsuga & Mineshige, 2007; Sugimura et al., 2016). In this paper, we investigate how intermediate massive BHs with $10^{3}\leq M_{\rm BH}\leq 5\times 10^{5}~{}{\rm M}_{\odot}$, embedded in an initially uniform gas with $n_{\infty}=10^{5}~{}{\rm cm}^{-3}$ and $T_{\infty}=10^{4}~{}{\rm 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 $M_{\rm BH}\ga 5\times 10^{5}~{}{\rm M}_{\odot}$, 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)

[TABLE]

where $c_{\infty}=\sqrt{\gamma\mathcal{R}T_{\infty}/\bar{\mu}}$ is the sound speed, $\gamma$ 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 $\gamma=1$ as

[TABLE]

and the Eddington accretion rate as $\dot{M}_{\rm Edd}\equiv L_{\rm Edd}/c^{2}$, where $L_{\rm Edd}=4\pi cGM_{\rm BH}/\kappa_{\rm es}$ and $\kappa_{\rm es}$ is the electron scattering opacity. Note that the Bondi radius and rate are calculated for $\gamma=1$, $\bar{\mu}=1.23$ and $T_{\infty}=10^{4}{\rm 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 onto the central BH. We adopt the spherical coordinates of $(r,\theta,\phi)$, defining the polar axis ($\theta=0$ and $\pi$) as directions perpendicular to the disk plane.

We solve two-dimensional hydrodynamical equations assuming axisymmetric flows; the equation of continuity

[TABLE]

and the equations of motion,

[TABLE]

[TABLE]

[TABLE]

where $\rho$ is the gas density, ${\boldsymbol{v}}=(v_{r},v_{\theta},v_{\phi})$ is the velocity, $p$ is the gas pressure, and $f_{\rm rad}$ 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 $\psi=-GM_{\rm BH}/r$.

We also solve the energy equation

[TABLE]

where $e=e_{\rm int}+\rho|{\boldsymbol{v}}|^{2}/2$ is the gas total energy density, $e_{\rm int}$ is the gas internal energy density, $\Lambda$ is the net cooling rate per volume and $\Gamma$ is the heating rate due to radiation from the central region. We assume the equation of state of ideal gas as $p=(\gamma-1)e_{\rm int}$ for $\gamma=5/3$.

We consider radiative cooling by bound-bound transitions of ${\rm 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 ${\rm e^{-}}$. The abundance of He nuclei relative to H nuclei is set to $8.33\times 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 time steps shorter than chemical timescales defined by $t_{\rm chem}\equiv(x_{\rm e}+0.001x_{\rm H})/\dot{x}_{\rm e}$, where $x_{\rm e}$ and $x_{\rm H}$ are the abundance of electrons and neutral hydrogens, respectively (Whalen & Norman, 2006, 2008).

We solve the multi-frequency, steady radiative transfer equation because the light crossing time is much shorter than the hydrodynamical timescale,

[TABLE]

where $F_{\nu}$ is the radiation flux, $E_{\nu}$ is the radiation energy density, and $\kappa_{\nu}$ is the absorption opacity. Note that the radial component of the radiation flux is calculated (see §4 for more details). Since the ionized gas is optically thin to electron scattering and bound-free transitions, we assume $F_{\nu}\approx cE_{\nu}$ on the right-hand-side of Eq. (9). The frequency range is set to $h\nu_{\rm min}(=13.6~{}{\rm eV})\leq h\nu\leq h\nu_{\rm max}(=10~{}{\rm keV})$, where $h$ is the Planck constant. The ionization rate coefficients $k_{\rm ph}$ of ${\rm H,He,He^{+}}$ and photoionization heating rate $\Gamma(=\Gamma_{\rm H})$ are calculated following the photon-conserving manner (Whalen & Norman, 2006). Here $\Gamma_{\rm 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

[TABLE]

In order to integrate Eq. (9), we simply set a radiation source with a single power-law spectrum at the center, $L_{\nu}=L_{0}(\nu/\nu_{\rm min})^{-\alpha}$ at $\nu_{\rm min}\leq\nu\leq\nu_{\rm max}$ and $L_{\nu}=0$ at $\nu\leq\nu_{\rm min}$ and $\nu_{\rm max}\leq\nu$, where $\alpha=1.5$ is adopted. The normalization factor $L_{0}$ 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 disk, which is unresolved, at the inner boundary,

[TABLE]

(Watarai et al., 2000), where $\dot{m}\equiv\dot{M}/\dot{M}_{\rm Edd}$. Furthermore, we assume anisotropic radiation fields as

[TABLE]

where the radiation flux $F_{\nu}(r_{\rm min})$ is normalized so that $L_{\nu}(r_{\rm min})=\int F_{\nu}(r_{\rm min})r_{\rm min}^{2}{\rm d}\Omega$ and $r_{\rm min}$ is the radius of the inner-most grid in the computation domain (see §2.2). Unlike Sugimura et al. (2016), shadow regions are not assumed in our simulations (see §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 $M_{\rm BH}=10^{3}~{}{\rm M}_{\odot}$. We also study cases for higher BH masses with $M_{\rm BH}=5\times 10^{4}~{}{\rm M}_{\odot}$ (Model D) and $5\times 10^{5}~{}{\rm M}_{\odot}$ (Model E) for $N=4$. Our models are summarized in Table 1.

2.2 Initial and boundary conditions

We set a computational domain of $r_{\rm min}\leq r\leq r_{\rm max}$ and $0\leq\theta\leq\pi$. Here we adopt $r_{\rm min}=0.035~{}R_{\rm B}$, $r_{\rm max}=30~{}R_{\rm B}$ for Models A, B, and C, and $r_{\rm min}=0.01~{}R_{\rm B}$, $r_{\rm max}=6~{}R_{\rm B}$ 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. Thus, we set logarithmically-spaced grids in the radial direction for Models A-C to simulate the flow within the Bondi radius in high resolution. We also employ power-law-spaced grids in Models D and E to resolve the structure of outflows in the outer region ($r>R_{\rm B}$), as well as inflows in the inner region ($r<R_{\rm B}$). We set uniformly-spaced grids in the polar direction. The number of the grid points is set to $(N_{r},N_{\theta},N_{\phi})=(100,120,1)$. Note that our simulations do not consider accretion flows within $r_{\rm min}$. Instead, we assume properties of radiation emitted from the central region (see §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_{\infty}=10^{5}~{}{\rm cm}^{-3}$ and temperature $T_{\infty}=10^{4}~{}{\rm 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=r_{\rm max}$ as the initial value for grids with a inflow velocity i.e., $v_{r}(r=r_{\rm max})<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

Figure 1 presents the time evolution of gas accretion rates onto a BH with $M_{\rm BH}=10^{3}~{}{\rm M}_{\odot}$ 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 $\sim 3.5~{}L_{\rm Edd}$. We continue the simulation until the ionization front reaches the outer boundary at $t\simeq 1.2\times 10^{5}$ 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.

In Figure 2, we show two-dimensional distribution of the gas density and temperature at $t=1.1\times 10^{5}$ 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}$, leading to episodic accretion (see Figure 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.

In the anisotropic case (Model C), the shape of the ionized region is no longer isotropic. While the ionization front expands toward the polar directions, the radiation flux in the equatorial direction ($\theta=\pi/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 onto electrons. The outflows collide with the ambient gas at $r\ga R_{\rm B}$ and form strong shocks with a high temperature of $\ga 5\times 10^{4}~{}{\rm K}$. The gas near the equator remains neutral and the temperature is kept at $T\simeq 8000~{}{\rm K}$ due to efficient atomic hydrogen cooling. Therefore, the gas can accrete onto the central BH through the neutral region increasing the density, unimpeded by radiation feedback. The radial profiles at $\theta=\pi/2$ of the density and inflow velocity approach those for a Bondi solution with $T\sim 8000~{}{\rm K}$ ($\rho\propto r^{-3/2}$ and $v_{\rm r}\propto r^{-1/2}$ at $r\la R_{\rm B}$).

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 $\Omega_{\rm HI}$ is the solid angle of the neutral region. Since a half angle of the neutral region is $\Theta_{\rm HI}\equiv\pi/2-\theta\approx 5^{\circ}$ 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 setups, Sugimura et al. (2016) assumed a shadow region with a half angle of $\Theta_{\rm Shad}=10^{\circ}-45^{\circ}$, 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\times 10^{5}~{}{\rm M}_{\odot}$ (Model E) but with the same radiation anisotropy ($N=4$) as in Model C. Figure 3 shows two-dimensional distribution of the gas density (left panels) and temperature (right panels) in the early stage at $t\simeq 0.15~{}t_{\rm dyn}$ (top panels) and the late stage at $t\simeq 0.49~{}t_{\rm dyn}$ (bottom panels), respectively. Initially, the ionized region expands to $r\sim R_{\rm B}$ preferentially toward the polar directions, and gas accretion is allowed through the neutral region (top panels in Figure 3). Unlike in Model C with lower $M_{\rm BH}$, however, neutral and dense gas inflows through the equator totally cover the central region. Since ionizing photons from the center are absorbed by the surrounding neutral gas, the radiation flux in all directions (even in the polar directions) is blocked. As a result, the ionized region is confined within $r\sim 0.04R_{\rm B}$ due to radiative recombination (bottom panels in Figure 3).

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 ($r_{\rm min}\leq r\la 0.05R_{\rm B}$). Figure 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 (after the transition). Figure 5 also presents profiles of the ionization degree, the ram pressure ($\rho v_{\theta}^{2}$) and the gas pressure ($p_{\rm gas}$) as functions of the polar angle at $r=0.02R_{\rm B}$. Before the transition shown in Figure 4 (panel i1), a geometrically thin and dense shock layer forms along the direction of $\theta\approx 45^{\circ}$, where the boundary between ionized and neutral gas is located (see Figure 5). The inflow from the neutral region crosses the shock with a supersonic velocity, namely $|v_{\theta}|>c_{\rm s}$ and is decelerated ($c_{\rm s}$ is the local sound speed). The shocked gas is ionized and accelerated outward in radial direction by strong radiation force. As shown in Figure 5, ram pressure of the neutral gas inflow slightly dominates 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 Figure 4 shows that a dense clump forms at the shock front by compression caused by continuous accretion of neutral gas and strong radiation force from the center. 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.

Figures 6, 7 and 8 present radial profiles of the gas density, temperature, and radial velocity along three different directions: $\theta=0^{\circ}$ (polar), $17^{\circ}$, and $90^{\circ}$ (equator) in the early stage at $t=0.15~{}t_{\rm dyn}$ (blue) and in the late stage at $t=0.49~{}t_{\rm dyn}$ (red). The black dashed curve represents the isothermal Bondi profile with $T=10^{4}{\rm K}$. 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 ($\theta=0^{\circ}$ and $17^{\circ}$). Near the edge of the ionized region, colliding outflows due to the radiation force on electrons (dashed curves in Figure 8) produce strong shocks with a high temperature of $\sim 10^{6}~{}{\rm K}$. On the other hand, neutral gas accretes through the equator keeping the temperature at $T\simeq 8000~{}{\rm K}$, as shown in §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 $\sim 10^{4}$ K in all directions because of radiative cooling by atomic hydrogen. Along the direction of $\theta=17^{\circ}$, the ionized region still exists but its size is too small ($\sim 0.04R_{\rm B}$) 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\simeq 8000$ K toward narrow conical regions along the polars ($\theta\sim 0^{\circ}-20^{\circ}$).

Figure 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 $\sim 5\times 10^{4}~{}L_{\rm Edd}/c^{2}$, which corresponds to $\sim 50~{}\%$ of the inflow rate from the Bondi radius. On the other hand, the accreting material is accelerated outward into the bipolar directions ($\theta\la 20^{\circ}$) due to absorption of the radiation momentum with $\gg L_{\rm Edd}/c$. The outflow rate ends up $\simeq 1.5\times 10^{4}~{}\dot{M}_{\rm Edd}$, which is $\sim 30~{}\%$ of the BH feeding rate. Since the Mach number of the outflowing matter is as high as five at the Bondi radius (see Figure 8), the typical velocity is estimated as $\simeq 50~{}{\rm km~{}s}^{-1}$. Sugimura et al. (2016) reported that hot and ionized outflows with $T\simeq 7\times 10^{4}~{}{\rm 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\simeq 8000~{}{\rm K}$ are likely to be produced because of the transition to the wholly neutral phase.

In Figure 10, we show the BH feeding rates for two cases with $M_{\rm BH}=5\times 10^{4}~{}{\rm M}_{\odot}$ (blue, Model D) and $5\times 10^{5}~{}{\rm M}_{\odot}$ (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, radiation force from the center dominates ram-pressure of the infalling clumps, and thus such clumps do not cover the central radiation source. Therefore, this result implies the critical BH mass for the transition is located around $M_{\rm BH}\sim 5\times 10^{5}~{}{\rm M}_{\odot}$ for the initial conditions we adopted ($n_{\infty}=10^{5}~{}{\rm cm}^{-3}$ and $T_{\infty}=10^{4}~{}{\rm K}$).

4 Summary and discussion

In this paper, we study rapid growth of BHs with a mass range of $10^{3}\leq M_{\rm BH}\leq 5\times 10^{5}~{}{\rm M}_{\odot}$ via gas accretion exposed to anisotropic super-Eddington radiation. Radiation flux from an accretion disk around the nuclear BH is likely to be collimated toward 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_{\infty}=10^{5}~{}{\rm cm}^{-3}$ and temperature of $T_{\infty}=10^{4}~{}{\rm 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 $M_{\rm BH}\ga 5\times 10^{5}~{}{\rm M}_{\odot}$, 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\simeq 8000~{}{\rm K}$.

We adopt a simple model for anisotropic radiation from the central accretion disk, namely, $F_{\rm rad}(\theta)\propto\cos^{N}\theta$. 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 disk is optically and geometrically thick, so that the radiation flux is mildly collimated to the polars (e.g., Ohsuga et al., 2005; Jiang, Stone & Davis, 2014; Sa̧dowski et al., 2015). However, while 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 setups are similar, mainly except two things. First, they assumed the existence of a shadow region caused by attenuation due to a geometrically thick disk and/or disk winds well inside their computation domains. Second, our simulations resolve even smaller scales, namely $\sim 1\%$ of the Bondi radius. We can stress that our setup 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 ($M_{\rm BH}=10^{5.7}~{}{\rm M}_{\odot}$). The gas profile approaches an isothermal Bondi solution with $T\sim 10^{4}$ 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}\la 116~{}M_{\rm BH,5.7}^{3/2}n_{\infty,5}(r_{\star}/10^{18}~{}{\rm cm})^{-1/2}$ (Sakurai, Inayoshi & Haiman, 2016), where $r_{\star}$ is the location of the photosphere, from which the radiation emerges actually. Note that the photosphere is located inside the innermost grid ($r_{\star}<r_{\rm min}=10^{18}~{}{\rm cm})$. Assuming that the BH is fed at the Bondi rate, the radiation luminosity toward the polar directions is estimated as $L/L_{\rm Edd}\approx 106$ (see Eq.12). Therefore, we expect that the overall behaviors of the accretion flows in Model E would be maintained.

As shown in Figure 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 region111 The non-radial components in radiation flux is unlikely to be produced by electron scattering because the ionized gas is optically thin even near the innermost region.. 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 10^{-18}(n_{\rm e}/10^{3}~{}{\rm cm}^{-3})^{2}~{}{\rm erg\ cm^{-3}s^{-1}}$ is the emissivity, $\langle h\nu_{\rm rec}\rangle$ is the average energy of recombination photons, $\alpha_{\rm B}$ is the case B recombination rate coefficient, and $n_{\rm e}$ is the numb density of electrons. Thus the non-radial radiation flux is estimated as $|F_{\rm rad}^{\theta}/F_{\rm rad}^{r}|\sim 10^{-8}$ 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 $\sim 10^{4}~{}{\rm 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\ga 10^{5}~{}{\rm cm}^{-3}$ and $T\sim 8000~{}{\rm K}$), the H2 fraction relative to atoms is determined by the balance between formation through electron-catalyzed reactions and collisional dissociation with gas particles (Omukai, 2001). The resultant fraction is estimated as $\sim 10^{-8}$ as long as the temperature is kept at $>4000~{}{\rm K}$, below which the collisional dissociation becomes inefficient (Inayoshi, Omukai & Tasker, 2014). As shown in Figure 7, the temperature cools down to $\simeq 1000-2000~{}{\rm K}$ due to adiabatic expansion at $r>2\times 10^{19}$ cm in the outflowing matter ($\theta=0^{\circ}$). In this lower-temperature region, far-ultraviolet radiation in the Lyman-Werner (LW) bands ($h\nu_{\rm 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

[TABLE]

where $f_{\rm LW}$ is the ratio of the LW flux to the bolometric value and $\langle\nu_{\rm LW}\rangle=12~{}{\rm eV}/h~{}(=2.9\times 10^{15}~{}{\rm Hz})$. 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 $F_{\rm LW}=f_{\rm shield}F_{\rm LW,thin}\simeq 7\times 10^{-14}~{}{\rm erg~{}s^{-1}~{}cm^{-2}~{}Hz^{-1}}$, which is higher than the critical LW flux222The 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 $\sim 2$ eV dissociate H- ions, which is an intermediate product in the H2 formation channels (${\rm H}+{\rm e}^{-}\rightarrow{\rm H}^{-}+\gamma$; ${\rm H}^{-}+{\rm H}\rightarrow{\rm H}_{2}+{\rm e}^{-}$) and contribute to H2 dissociation indirectly (Omukai, 2001)., $F_{\rm LW,crit}\sim 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\sim 10^{20}$ cm.

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 JSPS Grant-in-Aid for Scientific Research (C) (15K05036 K.O.) and for Young Scientists (17K14260 H.R.T). This research was also supported by MEXT as ”Priority Issue on Post-K computer” (Elucidation of the Fundamental Laws and Evolution of the Universe) and JICFuS. Numerical computations were carried out on Cray XC30 at Center for Computational Astrophysics, National Astronomical Observatory of Japan.