---
abstract: 'Discovery of high-redshift ($z > 6$) supermassive black holes (BHs) may indicate that the rapid (or super-Eddington) gas accretion has aided their quick growth. Here, we study such rapid accretion of the primordial gas on to intermediate-mass ($10^2 - 10^5~M_\odot$) BHs under anisotropic radiation feedback. We perform two-dimensional radiation hydrodynamics simulations that solve the flow structure across the Bondi radius, from far outside of the Bondi radius down to a central part which is larger than a circum-BH accretion disc. The radiation from the unresolved circum-BH disc is analytically modeled considering self-shadowing effect. We show that the flow settles into a steady state, where the flow structure consists of two distinct parts: (1) bipolar ionized outflowing regions, where the gas is pushed outward by thermal gas pressure and super-Eddington radiation pressure, and (2) an equatorial neutral inflowing region, where the gas falls toward the central BH without affected by radiation feedback. The resulting accretion rate is much higher than that in the case of isotropic radiation, far exceeding the Eddington-limited rate to reach a value slightly lower than the Bondi one. The opening angle of the equatorial inflowing region is determined by the luminosity and directional dependence of the central radiation. We find that photoevaporation from its surfaces set the critical opening angle of about ten degrees below which the accretion to the BH is quenched. We suggest that the shadowing effect allows even stellar-remnant BHs to grow rapidly enough to become high-redshift supermassive BHs.'
author:
- |
    Kazuyuki Sugimura,$^1$[^1] Takashi Hosokawa,$^{2,3,4}$ Hidenobu Yajima$^{1,5}$ and Kazuyuki Omukai$^{1,3}$\
    $^1$Astronomical Institute, Tohoku University, Aoba, Sendai 980-8578, Japan\
    $^2$Department of Physics, Kyoto University, Sakyo, Kyoto 606-8502, Japan\
    $^3$Kavli Institute for Theoretical Physics, University of California, Santa Barbara, California 93106, USA\
    $^4$Department of Physics and Research Center for the Early Universe, the University of Tokyo, Bunkyo, Tokyo 113-0033, Japan\
    $^5$Frontier Research Institute for Interdisciplinary Sciences, Tohoku University, Aoba, Sendai 980-8578, Japan
title: Rapid Black Hole Growth under Anisotropic Radiation Feedback
---

-1cm

quasars: supermassive black holes-cosmology: theory.

Introduction {#sec:intro}
============

Discovery of high-$z$ ($z\gtrsim6$) quasars suggests that supermassive black holes (SMBHs) already exist when the age the Universe is less than $1{\,\mathrm{Gyr}}$ [see, e.g., @Fan:2001aa; @Willott:2010aa; @Mortlock:2011aa; @Venemans:2013aa; @Wu:2015aa]. This poses a question about the formation mechanism of SMBHs in such a short interval. Among the scenarios for the SMBH seed formation [see, e.g., @Volonteri:2012ab; @Haiman:2013aa for a review], including the dense stellar cluster scenario [see, e.g., @Omukai:2008aa; @Devecchi:2009aa; @Katz:2015ab; @Yajima:2016aa and reference therein], following two are the most studied: the direct collapse BH (DCBH) and the population III (Pop III) remnant BH scenarios.

In the former scenario, supermassive stars of $\sim 10^5\,M_\odot$ collapse to form seed BHs with approximately the same mass. Specifically, supermassive stars are envisaged to form in exceptional environments in the high-$z$ Universe, for example, in atomic-cooling halos where the ${\mathrm{H_2}}$ cooling is totally suppressed by very strong far ultraviolet (FUV) irradiation [e.g., @Sugimura:2014aa]. While the seed BHs in this case are rather massive with $\sim
10^5\,M_\odot$, their number density might be too small to explain all the observed high-$z$ SMBHs due to the stringent necessary conditions [@Dijkstra:2008aa; @Dijkstra:2014aa; @Agarwal:2012aa; @Sugimura:2014aa; @Sugimura:2016aa; @Inayoshi:2015ab; @Chon:2016aa].

In the latter scenario, the remnant BHs of Pop III stars [@Yoshida:2008aa; @Hosokawa:2011aa; @Hosokawa:2016aa] are thought as SMBH seeds [@Alvarez:2009aa; @Jeon:2012aa]. Contrary to the DCBH scenario, they are abundant but the problem is whether they can actually grow to the SMBHs from smaller initial mass of $\lesssim 10^3\,M_\odot$ [@Susa:2014aa; @Hirano:2015aa] within the available time. Although BHs can acquire the mass by collisions with other BHs [@Tanikawa:2011aa], the BH collisions often result in ejection of the merged BHs from the host halo due to the recoil of gravitational wave emission [e.g., @Baker:2006aa; @Koppitz:2007aa]. Thus, the feasibility of this scenario relies on whether the rapid accretion on to seed BHs is possible or not [@Madau:2014aa; @Alexander:2014aa; @Volonteri:2015aa].

Recently, a number of authors have studied the BH accretion under radiation feedback [e.g., @Milosavljevic:2009aa; @Milosavljevic:2009ab; @Park:2011aa; @Park:2012aa; @Park:2013aa]. They solve the gas dynamics over the scale of the Bondi radius, where the accretion rate on to the circum-BH disc is physically determined. Although the central circum-BH disc is not spatially resolved, subgrid models that provide analytic prescriptions of its emissivity have been used. They have shown that the accretion rate is significantly reduced to $\lesssim1\%$ of that without radiation feedback (i.e., the Bondi rate) in case with modest BH mass and ambient density (e.g., $10^2\,M_\odot$ and $10^5{\,\mathrm{cm^{-3}}}$). Only in case with very high BH mass and/or ambient density (e.g., $10^4\,M_\odot$ and $10^5{\,\mathrm{cm^{-3}}}$), the accretion rate reaches to the Bondi value because of inefficient radiation feedback, as recently shown by [@Inayoshi:2016ac] [see also @Li:2011aa; @Pacucci:2015ab; @Park:2016ab for other mechanisms of efficient accretion]. However, all those calculations assume isotropic radiation (in either one- or two-dimensional simulations), whereas in reality the radiation from the BH accretion disc should be anisotropic. The flow structure will be significantly altered in such anisotropic radiation field. Although the BH accretion under anisotropic radiation has been studied in the context of active galactic nuclei (AGN) with the BH mass $\ga 10^6\,M_\odot$ [@Proga:2007aa; @Kurosawa:2009aa; @Novak:2011aa; @Barai:2012aa], the nature of accretion on to stellar-mass BHs would be quite different.

The anisotropic BH irradiation has been examined with different models of the BH accretion discs, including the “standard disc” for moderate accretion rates [@Shakura:1973aa], and “slim disc” for the higher rates [@Abramowicz:1988aa]. In particular, recent multi-dimensional simulations have investigated inner structure of the slim disc within roughly a hundred Schwarzschild radii, showing that the accretion rates can indeed exceed the Eddington-limited rate [e.g., @Ohsuga:2005aa; @Jiang:2014aa; @McKinney:2014aa; @Fragile:2014aa; @Takahashi:2015aa; @Sc-adowski:2016aa]. These studies show that the high-energy photons are predominantly emitted in polar directions from the inner part of the disc. However, the outer structure of the disc, which is not solved in the above simulations, should also modify the anisotropic radiation field. For instance, disc winds such as the line-driven AGN winds launched from the outer region will absorb a part of photons coming from the inner region [e.g., @Proga:2000aa; @Proga:2004aa; @Nomura:2016aa]. Since numerical simulations solving the whole structure of the disc are still infeasible, it is very uncertain how much anisotropy the BH accretion discs actually create.

In this paper, we will investigate accretion of the primordial gas on to BHs under the anisotropic radiation feedback from the central circum-BH accretion discs, considering the shadowing effect by the outer part of the discs. We perform a set of proof-of-concept two-dimensional (2D) radiation hydrodynamics (RHD) simulations, assuming that BHs are initially embedded in homogeneous and static media. We do not attempt to simulate the realistic directional dependence of BH irradiation in consideration of its high uncertainties; instead, we model it in a simple fashion to study how the anisotropy of radiation changes the nature of accretion flows. As confirmed later by our results, the shadowing effect dramatically enhances the accretion rate. This mechanism might give a new pathway from the remnant BHs of Pop III stars to SMBHs within a limited timescale of $\la 1{\,\mathrm{Gyr}}$ after the Big Bang.

The paper is organized as follows. In Sec. \[sec:sph\_acc\], we briefly review the basics of spherical gas accretion on to a BH. In Sec. \[sec:num-method\], we describe the numerical method and cases considered. In Sec. \[sec:result\], we present the main results of our simulations. The conclusions and discussions are given in Sec. \[sec:conclusion\].

Basics {#sec:sph_acc}
======

For later reference, we first briefly summarize the basics of spherical gas accretion on to a central BH under radiation feedback. We consider a system where a BH is embedded in a static and homogeneous medium. We take the BH mass $M_{\mathrm{BH}}=10^3\,M_\odot$, ambient density $n_\infty=10^5\,{\mathrm{cm^{-3}}}$ and ambient temperature $T_{\mathrm{HI}}=10^4\,{\mathrm{K}}$ as a fiducial parameter set.

If we ignore the effect of feedback, the mass accretion will proceed at the Bondi rate in this case, $$\begin{aligned}
\dot{M}_{\mathrm{B}}
&=\frac{4\pi\lambda_{\mathrm{B}}\rho_\infty G^2 M_{\mathrm{BH}}^2}{c_{\mathrm{s,HI}}^3}{\nonumber\\}&=1.7\times 10^{-3}
\,\left(\frac{n_\infty}{10^5\,{\mathrm{cm^{-3}}}}\right){\nonumber\\}&\qquad \times \left(\frac{M_{\mathrm{BH}}}{10^3\,M_\odot}\right)^{2}
 \left(\frac{T_{\mathrm{HI}}}{10^4\,{\mathrm{K}}}\right)^{-3/2}
 {\mathrm{M_\odot\, yr^{-1}}}
 \,,\label{eq:1}\end{aligned}$$ where we take $\lambda_{\mathrm{B}}=(1/4)\left[2/(5-3\gamma)\right]^{(5-3\gamma)/2(\gamma-1)}
= 1.12$ assuming the gas is isothermal (the polytropic index $\gamma=1$). For a neutral primordial gas with helium-to-hydrogen ratio in the number of nuclei $y_{\mathrm{He}}=0.0972$, the mean molecular weight $\mu=(1+4y_{\mathrm{He}})/(1+y_{\mathrm{He}})=1.3$, the mass density of the medium $\rho_\infty=n_\infty(1+4y_{\mathrm{He}})m_{\mathrm{p}}=2.3\times10^{-19}{\,\mathrm{g\,cm^{-3}}}$ with $m_{\mathrm{p}}$ the proton mass and the (isothermal) sound speed $c_{\mathrm{s,HI}}=(k_{\mathrm{B}}T_{\mathrm{HI}}/\mu m_{\mathrm{p}})^{1/2}
=8.1\,\left(T_{\mathrm{HI}}/10^4\,{\mathrm{K}}\right)^{1/2} \,{\mathrm{km\, s^{-1}}}$. The Bondi radius, defined as $$\begin{aligned}
 r_{\mathrm{B}}
&=\frac{GM_{\mathrm{BH}}}{c_{\mathrm{s,HI}}^2}{\nonumber\\}&=1.4\times 10^{4} 
 \left(\frac{M_{\mathrm{BH}}}{10^3\,M_\odot}\right)
 \left(\frac{T_{\mathrm{HI}}}{10^4\,{\mathrm{K}}}\right)^{-1}
 {\mathrm{AU}}\,,
\label{eq:2}\end{aligned}$$ demarcates the inner region where the gravitational energy dominates the thermal energy and the outer region where the thermal energy dominates. Correspondingly, the gas is approximately in free fall inside $r_{\mathrm{B}}$, whereas the pressure equilibrium is almost achieved outside.

The Eddington luminosity $L_{\mathrm{E}}$ is the critical luminosity above which the outward radiation force via the Thomson scattering exceeds the inward gravitational pull of the BH in fully ionized hydrogen gas, $$\begin{aligned}
L_{\mathrm{E}}=\frac{4\pi G M_{\mathrm{BH}} c m_{\mathrm{p}}}{\sigma_{\mathrm{T}}}
&=3.3\times 10^7
\left(\frac{M_{\mathrm{BH}}}{10^3M_\odot}\right) L_\odot\,,
\label{eq:12}\end{aligned}$$ where $\sigma_{{\mathrm{T}}}$ is the Thomson scattering cross section. Note that the Eddington luminosity does not always provide physical limit because gas pressure is not considered in the above argument. In addition, the radiation force becomes less effective in a partially ionized gas.

The (efficiency-independent) Eddington-limited accretion rate is defined as $$\begin{aligned}
 \dot{M}_{\mathrm{E}}
&=\frac{L_{\mathrm{E}}}{c^2}
=2.2\times 10^{-6} \left(\frac{M_{\mathrm{BH}}}{10^3M_\odot}\right)\,M_\odot\, {\mathrm{yr^{-1}}}\,,
\label{eq:8}\end{aligned}$$ and the condition for the luminosity to be sub-critical can be rewritten as $\dot{M}<\dot{M}_{\mathrm{E}}/\eta$ with the radiative efficiency $\eta$. The radiative efficiency $\eta\approx 0.1$ for a standard accretion disc is widely used in the previous works. [see, e.g., @Milosavljevic:2009aa; @Park:2011aa; @Park:2012aa]. Note that in some literatures the efficiency-dependent Eddington-limited accretion rate, $\dot{M}_{\mathrm{E}}/\eta$ in our definition, is used instead. For large $M_{\mathrm{BH}}$ and/or $n_\infty$, the Bondi rate $\dot{M}_{\mathrm{B}}$ can be much larger than the Eddington rate $\dot{M}_{\mathrm{E}}$, e.g., $\dot{M}_{\mathrm{B}}/\dot{M}_{\mathrm{E}}=7.8\times 10^{2}
\,\left(n_\infty/10^5\,{\mathrm{cm^{-3}}}\right)
\left(M_{\mathrm{BH}}/10^3\,M_\odot\right)$, because $\dot{M}_{\mathrm{B}}$ is proportional to $M_{\mathrm{BH}}^2\,n_\infty$ while $\dot{M}_{\mathrm{E}}$ to $M_{\mathrm{BH}}$.

High energy photons emitted by the BH accretion disk create a surrounding [H[ii]{} ]{}bubble. With the power-law spectrum $L_\nu \propto
\nu^{-1.5}$, which is often postulated in the literature, the photoionized gas is heated up to $T_{\mathrm{HII}}\sim 7\times10^4{\,\mathrm{K}}$ owing to helium ionization heating. The high thermal pressure of the [H[ii]{} ]{}bubble, together with the outward radiation pressure, can significantly reduce the accretion rate [@Milosavljevic:2009aa; @Park:2011aa; @Park:2012aa]. The size of the [H[ii]{} ]{}bubble is estimated by the Strömgren radius, $$\begin{aligned}
  r_{\mathrm{HII}} &=
6.8\times 10^{4}
\left(\frac{T_{\mathrm{HII}}}{7\times10^4{\,\mathrm{K}}}\right)^{\frac{1}{3}}{\nonumber\\}&\qquad\times\left(\frac{L}{3.3\times10^7L_\odot}\right)^{\frac{1}{3}}
\left(\frac{n_{\mathrm{HII}}}{10^5\,{\mathrm{cm^{-3}}}}\right)^{-\frac{2}{3}} 
 {\,\mathrm{AU}}\,,
\label{eq:13}\end{aligned}$$ which is obtained by equating the ionizing photon emissivity $\dot{N}_{\mathrm{ion}}=\int_{\nu_{\mathrm{T}}}^\infty d\nu L_\nu/h\nu$ ($=L/3h\nu_{\mathrm{T}}$ for the above spectrum with $L_\nu \propto
\nu^{-1.5}$) with the recombination rate within the [H[ii]{} ]{}bubble $\alpha_{\mathrm{B}}(4\pi/3)r_{\mathrm{HII}}^3 n_{\mathrm{HII}}^2$, where $n_{\mathrm{HII}}$ is the number density of hydrogen nuclei inside the bubble, $h\nu_{\mathrm{T}}=13.6{\,\mathrm{eV}}$ the hydrogen ionization energy, and $\alpha_{\mathrm{B}}=4.6\times10^{-14}{\,\mathrm{cm^3\,s^{-1}}}$ the case B hydrogen recombination coefficient at $7\times10^4{\,\mathrm{K}}$ [@Ferland:1992aa]. In the above, we take $n_{\mathrm{HII}}=n_\infty\,(=10^5{\,\mathrm{cm^{-3}}})$ and $L=L_{\mathrm{E}}$ as reference values.

Suppose that an ionizing source is suddenly turned on at the centre. The ionization front first propagates up to $r_{\mathrm{HII}}$ with $n_{\mathrm{HII}}$ kept almost constant. Then the bubble expands until pressure equilibrium with the surrounding medium is reached with $n_{\mathrm{HII}}=(c_{\mathrm{s,HI}}/c_{\mathrm{s,HII}})^2\,n_\infty \sim 0.07\,n_\infty$, where $c_{\mathrm{s,HII}} = (2T_{\mathrm{HII}}/T_{\mathrm{HI}})^{1/2}c_{\mathrm{s,HI}}$ is the sound speed of the ionized gas, with the factor of $2$ accounting for the increase of the particle number by ionization. Note that we have neglected the effect of helium in estimating $r_{\mathrm{HII}}$ and $n_{\mathrm{HII}}$, because it modifies them only slightly. For the flow from the [H[ii]{} ]{}bubble in pressure equilibrium with the surrounding neutral medium, the Bondi radius and rate are given by $r_{\mathrm{B,HII}}\sim
0.07\,r_{\mathrm{B}}$ and $\dot{M}_{\mathrm{B,HII}} \sim 1\times 10^{-3}
\dot{M}_{\mathrm{B}}$, respectively. This clearly shows that the photo-ionization feedback can considerably suppress the accretion.

In order for the [H[ii]{} ]{}bubble to be trapped around the BH, however, the condition $r_{\mathrm{HII}}>r_{\mathrm{B}}$ must be satisfied [@Inayoshi:2016ac]. Otherwise, the gas originally in between $r_{\mathrm{HII}}$ and $r_{\mathrm{B}}$ would accumulate around the periphery of the [H[ii]{} ]{}bubble. This leads to the enhancement of density $n_{\mathrm{HII}}$ and thus to shrinkage of the bubble. In the end, the [H[ii]{} ]{}bubble disappears and radiation feedback no longer affects the accretion. At that time, the accretion rate returns to the original Bondi value for the neutral gas $\dot{M}_{\mathrm{B}}$, instead of that for the [H[ii]{} ]{}bubble $\dot{M}_{\mathrm{B,HII}}$. For this to happen, a system with a massive BH and/or dense ambient medium, namely $(M_{\mathrm{BH}}/10^4
M_\odot)(n_\infty/10^5\,{\mathrm{cm^{-3}}}) \gtrsim 1$ [@Inayoshi:2016ac], is required when the BH luminosity is close to the Eddington value as $L\approx L_{\mathrm{E}}$.

Recall that the above argument is based on the assumption of the spherical symmetry, which should be modified in realistic situations with anisotropic BH irradiation. We expect that the flow through shadowed equatorial regions, if exist, enhance the accretion rate. In what follows, we will see what kind of the flow structure appears for such cases by using numerical simulations.

NUMERICAL METHOD {#sec:num-method}
================

We study accretion of primordial gas on to BHs under anisotropic radiation by performing a series of 2D RHD simulations (Sec. \[sec:rhd\_sim\]). Specifically, we solve the dynamics of the flow around the Bondi radius (see Fig. \[fig:acc\_whole\]b), where the accretion rate on to the BH and disc system is determined. We mask the inner circum-BH accretion disc (see Fig. \[fig:acc\_whole\]a) by the central sink region and inject ionizing photons at the inner boundary $R_{\rm in}$ according to a simple parametric sub-grid model that represents various directional dependences of BH irradiation (Sec. \[sec:radiation\]).

Two-dimensional radiation hydrodynamics simulations {#sec:rhd_sim}
---------------------------------------------------

We use a modified version of the public multi-dimensional magneto-hydrodynamics code [Pluto 3.0]{} [@Mignone:2007aa], which has been applied to studies on the present-day high-mass star formation [e.g., @Kuiper:2010aa; @Kuiper:2010ab; @Kuiper:2011aa; @Kuiper:2013aa] and Pop III star formation [@Hosokawa:2016aa].

Here, we have tuned the code used for the Pop III star formation [@Hosokawa:2016aa] to fit our study of the BH accretion. As in @Kuiper:2010aa, we adopt a 2D polar coordinate system assuming the axial symmetry. We calculate only the gravity of the central BH and neglect the gas self-gravity, as in the previous studies [e.g., @Park:2011aa; @Milosavljevic:2009ab; @Inayoshi:2016ac but also see Li 2011]. We assume that the outer edge of the accretion disc, i.e., the centrifugal radius $r_{\rm cen}$, is much smaller than the sink radius $R_{\mathrm{in}}$. We thus ignore the angular momentum of the flow in the computational domain. Other modifications we have added are summarized as follows.

### Chemical and thermal processes {#sec:chemistry}

To solve the chemical and thermal processes, we use the same methods developed in [@Hosokawa:2016aa] with several modifications. Unlike in [@Hosokawa:2016aa], we omit ${\mathrm{H}}_2$ chemistry assuming that ${\mathrm{H}}_2$ is completely photo-dissociated by the central FUV irradiation.[^2] We have added the ${\mathrm{He}}$ chemistry, since hard UV photons from BH accretion discs create a large helium photoionized region embedded in an [H[ii]{} ]{}region.

In summary, we solve the chemical network with six species: ${\mathrm{H}}$, ${\mathrm{H^+}}$, ${\mathrm{e}}$, ${\mathrm{He}}$, ${\mathrm{He^+}}$, and ${\mathrm{He^{2+}}}$, which consists of the following chemical processes: photoionization of ${\mathrm{H}}$, ${\mathrm{He}}$ and ${\mathrm{He^+}}$; collisional ionization of ${\mathrm{H}}$, ${\mathrm{He}}$ and ${\mathrm{He^+}}$; recombination of ${\mathrm{H^+}}$, ${\mathrm{He^+}}$ and ${\mathrm{He^{2+}}}$. Accordingly we consider the following thermal processes: photoionization heating of ${\mathrm{H}}$, ${\mathrm{He}}$, and ${\mathrm{He^+}}$; recombination cooling of ${\mathrm{H^+}}$, ${\mathrm{He^+}}$, and ${\mathrm{He^{2+}}}$; excitation cooling of ${\mathrm{H}}$, ${\mathrm{He}}$, and ${\mathrm{He^+}}$; collisional ionization cooling of ${\mathrm{H}}$, ${\mathrm{He}}$, and ${\mathrm{He^+}}$; free-free cooling of ${\mathrm{H}}$, ${\mathrm{He}}$, and ${\mathrm{He^{+}}}$; Compton cooling by cosmic microwave background (CMB) photons. Complete lists of our adopted chemical and thermal processes are available in Appendix \[sec:chem\_detail\].

We turn off the cooling when the temperature falls below $10^4{\,\mathrm{K}}$, as in the previous 2D simulations [e.g., @Park:2011aa]. We neglect secondary ionization and heating caused by X-ray photoionization [@Shull:1979ab; @Shull:1985aa; @Ricotti:2002aa]. We have confirmed with test calculations that these processes hardly affect the gas dynamics though the ionization degree is only slightly enhanced just outside the [H[ii]{} ]{}bubble.

### Transfer of ionizing photons {#sec:rad_transfer}

As in @Hosokawa:2016aa, we only solve the transfer of ionizing photons directly coming from the central accretion disc. Diffuse recombination photons are considered by way of the on-the-spot approximation. The radiation transfer is successively solved with the chemistry from the innermost cell, where photons are injected according to the sub-grid radiation model (see Sec. \[sec:radiation\] below). We do not consider the absorption between the radiation source and the inner boundary, which is currently masked by the sink cell.

Regarding the transport of ionizing photons, we have made the following major updates. First, we solve the frequency-dependent transfer with 128 logarithmically-spaced frequency bins between $13.6{\,\mathrm{eV}}$ and $1{\,\mathrm{keV}}$, to consider the photoionization of ${\mathrm{H}}$, ${\mathrm{He}}$ and ${\mathrm{He^+}}$ with different threshold energies. Second, we consider the radiation pressure via Thomson scattering and photoionization. As will be seen in Section \[sec:result\], the radiation pressure becomes important when the luminosity exceeds the Eddington limit.

We simply assume that photons with energy below $13.6{\,\mathrm{eV}}$ freely escape from the system [@Park:2011aa]. Although the radiation pressure of accumulated Ly$\alpha$ photons would affect the gas dynamics in spherically symmetric systems, it is probably not the case in realistic systems with channels for Ly$\alpha$ photons to escape [e.g., @McKee:2008aa; @Milosavljevic:2009aa].

Subgrid model for the irradiation by BH {#sec:radiation}
---------------------------------------

In our model of anisotropic BH irradiation, we assume that ionizing photons are emitted from the inner hot part of a circum-BH accretion disc but a portion of them are absorbed (or scattered) by outer structures (see Fig. \[fig:acc\_whole\]a). We inject ionizing photons at the inner boundary depending on the inflow rate into the sink cell, according to the model described here. We first describe the structure of the BH accretion disc in Section \[sec:acc\_disc\], which motivates our subgrid model. Then we give the expressions for the luminosity and directional dependence in Sections \[sec:efficiency\] and \[sec:direction\_dep\], respectively.

### BH accretion disc with shadowing effect {#sec:acc_disc}

Fig. \[fig:acc\_whole\](a) shows the expected inner structure including a BH accretion disc that motivates our subgrid model. Below we explain the inner and outer parts of the structure presented in Fig. \[fig:acc\_whole\](a) in this order. We also describe resulting directional dependences of the BH irradiation.

In the inner part, we see that the circum-BH disc consists of the two different types of accretion discs. One is the innermost geometrically thick (aspect ratio $\sim 1$) slim disc appearing inside the photon trapping radius $r_{\mathrm{tr}} \equiv (\dot{M}/\dot{M}_{\mathrm{E}})\,r_{\mathrm{Sch}}$, where the cooling via radial advection balances with the viscous heating [e.g., @Begelman:1978aa; @Abramowicz:1988aa]. The other is the geometrically thin (aspect ratio $\ll 1$) standard accretion disc appearing outside $r_{\mathrm{tr}}$, where the radiative loss from the disc surfaces is the main cooling process [e.g., @Shakura:1973aa]. When $\dot{M}/\dot{M}_{\mathrm{E}} < 1$, the slim disc disappears and the standard accretion disc extends all the way to the inner disk edge. We model the luminosity based on this consideration in Section \[sec:efficiency\]. Since the surface temperature of the disc increases with decreasing the radius $r$, ionizing photons mostly come from the hot innermost part.

In the outer part, a disc wind might be launched from the disc surface photo-heated by the high-energy photons from the inner region [see @Proga:2000aa; @Proga:2004aa; @Nomura:2013aa; @Nomura:2016aa for line-driven disc wind of AGNs]. In addition, around the outer edge of the disc, the vertically falling flow due to the centrifugal barrier might collide with the one coming from the opposite side of the equatorial plane and form a shocked region. Fig. \[fig:acc\_whole\](a) depicts these structures, both of which can absorb (or scatter) the ionizing photons coming from the inner part, forming a shadowed region behind them. We assume that the shadowing effect considered here is caused by the outer disc structures, and thus the appearance of a slim disc is not essential in forming the shadowed region. The outer structure of the disc is highly uncertain and probably varies depending on the BH mass, accretion rate, angular momentum, metallicity of inflowing gas, etc.. In Section \[sec:direction\_dep\], we model the shadowing effect with a simple parametric fashion.

### Luminosity {#sec:efficiency}

In our simulation, we determine the luminosity of the BH radiation $L$ depending on $\dot{M}$ evaluated at the inner boundary at each time step. To model the luminosity, we adopt the fitting formula [@Watarai:2000aa], $$\begin{aligned}
 L = 
\begin{cases}
{\displaystyle}2\, L_{\mathrm{E}}\,\left[1+\ln\left(\frac{\dot{m}}{20}\right)\right] & \dot{m} > 20\\
{\displaystyle}0.1\, L_{\mathrm{E}}\, \dot{m}  & \dot{m} < 20
\end{cases}\,,
\label{eq:3}\end{aligned}$$ where $\dot{m}\equiv \dot{M}/\dot{M}_{\mathrm{E}}$. [@Watarai:2000aa] obtained this formula by fitting the dependence of the luminosity on the accretion rate in the 1D stationary disc model, taking into account the (dis-)appearance of the slim disc depending on $\dot{M}$ . When $\dot{M}$ is low ($\dot{m} < 20$), the radiative efficiency is fixed at 10%, which agrees with that of the standard disc. For rapid accretion with $\dot{m} > 20$, the second term $2 L_{\mathrm{E}}\, \ln(\dot{m}/20)$ represents the luminosity from the innermost slim disc, where the photon advection reduces the radiative efficiency. Note that the luminosity $L$ increases logarithmically with $\dot{M}$ and can even exceed $L_{\mathrm{E}}$ because a large fraction of the emitted photons escape from the disc surfaces in vertical directions [see e.g., @Abramowicz:1988aa; @Watarai:2000aa; @Ohsuga:2005aa; @Jiang:2014aa; @Sc-adowski:2016aa]. The first term $2 L_{\mathrm{E}}$ corresponds to the luminosity from the outer standard disc in $r > r_{\mathrm{tr}}$ (see Fig. \[fig:acc\_whole\]a), given approximately by the energy generation rate due to the gravitational energy released by $r_{\mathrm{tr}}$, $G\dot{M}/r_{\mathrm{tr}} \sim L_{\mathrm{E}}$ [e.g., @Begelman:1978aa; @Kato:1998aa].

The spectrum of the BH radiation is simply assumed to be the power-law with $L_\nu\propto \nu^{-1.5}$ for $h\nu>13.6{\,\mathrm{eV}}$, where $L=\int_{h\nu>13.6{\mathrm{eV}}}L_\nu d\nu$, as often assumed in the literature [e.g., @Park:2011aa; @Park:2012aa; @Milosavljevic:2009ab]. [@Park:2011aa] have shown that the qualitative properties of accretion do not depend on the spectral shape.

### Directional dependence {#sec:direction_dep}

We inject ionizing photons at the inner boundary with the directional dependence described below. Specifically, we multiply the anisotropy factor $\mathcal{F}(\theta)$ normalized as $\int
\mathcal{F}(\theta)d\Omega=4\pi$ with an isotropic radiation flux $L/4\pi R_{\mathrm{in}}^2$ at the inner boundary $R_{\mathrm{in}}$. With this definition, $\mathcal{F}(\theta) = 1$ represents the isotropic radiation (Fig. \[fig:dirdep\]). We use the latitudinal angle $\theta$ defined as the angle measured from the equatorial plane for our convenience.

Motivated by the expected disc structure described in Section \[sec:acc\_disc\] (also see Fig. \[fig:acc\_whole\]a), we model $\mathcal{F}(\theta)$ as $$\begin{aligned}
 \mathcal{F}(\theta) &= C\, f_{\mathrm{disc}}(\theta)\, f_{\mathrm{shadow}}(\theta)\,,
\label{eq:5}\end{aligned}$$ where $C$ is the normalization factor. In this expression, the inner anisotropy factor $f_{\mathrm{disc}}$ that represents the directional dependences of the radiation emitted from the inner part of the disc is multiplied by the outer one $f_{\mathrm{shadow}}$ to take into account the outer shadowing effect.

For the inner anisotropy factor $f_{\mathrm{disc}}$, we simply assume $$\begin{aligned}
f_{\mathrm{disc}}(\theta)\propto\sin\theta\,,
\label{eq:11} \end{aligned}$$ which corresponds to radiation from an infinitely thin disc (recall that we define $\theta$ as the angle from the equatorial plane). Although numerical simulations suggest somewhat steeper $\theta$-dependence especially in the polar directions [e.g., @Ohsuga:2005aa; @Sc-adowski:2016aa], such deviations cause little effects on our results because the mass accretion predominantly occurs through the infalling region near the equatorial plane. For the disc radiation without the outer shadowing effect (i.e., $f_{\mathrm{shadow}}=1$), the normalized anisotropy factor is $\mathcal{F}(\theta) = 2\sin \theta$ (Fig. \[fig:dirdep\]).

We model the outer anisotropy factor $f_{\mathrm{shadow}}$ as $$\begin{aligned}
 f_{\mathrm{shadow}}(\theta)
 =  
\begin{cases}
{\displaystyle}\exp\left[-\left(\frac{\theta-\tilde{\theta}_{\mathrm{shadow}}}{\delta\theta}\right)^2\right]
&0 < \theta < \tilde{\theta}_{\mathrm{shadow}}\\[0.4cm]
1  &  \tilde{\theta}_{\mathrm{shadow}} < \theta < 90^\circ
\end{cases}
\label{eq:10}\end{aligned}$$ where $\tilde{\theta}_{\mathrm{shadow}}=\theta_{\mathrm{shadow}}+2\,\delta\theta$, $\theta_{\mathrm{shadow}}$ is the opening angle of the shadow, and $\delta\theta$ the thickness of the transition region. Here, we assume $f_{\mathrm{shadow}}$ is symmetric about the equatorial plane. We adopt the finite transition region setting $\delta\theta = 6^\circ$ to avoid artificial ionization structure that appears with $\delta\theta\rightarrow 0$. Our conclusions are independent of the arbitrary choice of a small value for $\delta\theta$. We show $\mathcal{F}(\theta)$ for the disc radiation with the outer shadowing effect with $\theta_{\mathrm{shadow}}=45^\circ$ and $22.5^\circ$ in Fig. \[fig:dirdep\]. With the expression given by equation , the outer anisotropy factor begins to decrease even for $\theta > \theta_{\mathrm{shadow}}$, and takes a value of $\sim 0.01$ at $\theta = \theta_{\mathrm{shadow}}$. Although we fix the shadowing profile $f_{\mathrm{shadow}}(\theta)$ during each simulation run for simplicity, it probably depends on accretion rates in reality.[^3] In view of large uncertainties in the shadowing effect, we perform a number of simulations varying $\theta_{\mathrm{shadow}}$ as a free parameter (see Sec. \[sec:parameters\]).

Cases considered {#sec:parameters}
----------------

[lccccccccccc]{} run& $M_{\mathrm{BH}}\,[M_\odot]$ & $n_\infty\,[{\mathrm{cm^{-3}}}]$ & $\theta_{\mathrm{shadow}}^{a}$& $N_r\times N_\theta$ & $R_{\mathrm{in}}\,[{\mathrm{AU}}]$ & $R_{\mathrm{out}}\,[{\mathrm{AU}}]$ & $t_{\mathrm{end}}\,[{\mathrm{yr}}]$\
Di & $10^3$ & $10^5$ & [**isotropic**]{}$^{b}$ & $512\times144$ & $3\times10^2$ & $6\times10^5$ & $5\times10^5$\
Ddn & $10^3$ & $10^5$ & [**disc**]{}$^{c}$ & $512\times144$ & $3\times10^2$ & $6\times10^5$ & $5\times10^5$\
Dds$^{d}$ & $10^3$ & $10^5$ & [$45^\circ$]{} & $512\times144$ & $2\times10^3$ & $3\times10^6$ & $2\times10^6$\
\
s075 & $10^3$ & $10^5$ & [$33.75^\circ$]{} & $256\times72$ & $2\times10^3$ & $3\times10^6$ & $2\times10^6$\
s050 & $10^3$ & $10^5$ & [$22.5^\circ$]{} & $256\times72$ & $2\times10^3$ & $3\times10^6$ & $2\times10^6$\
s025 & $10^3$ & $10^5$ & [$11.25^\circ$]{} & $256\times72$ & $2\times10^3$ & $3\times10^6$ & $2\times10^6$\
\
M1e2 & [$10^2$]{} & $10^5$ & $45^\circ$ & $256\times72$ & $2\times10^2$ & $1.5\times10^6$ & $2\times10^6$\
M1e4 & [$10^4$]{} & $10^5$ & $45^\circ$ & $256\times72$ & $2\times10^4$ & $2\times10^7$ & $2\times10^7$\
M1e5 & [$10^5$]{} & $10^5$ & $45^\circ$ & $256\times72$ & $2\times10^5$ & $1\times10^8$ & $5\times10^7$\
\
n1e3 & $10^3$ & [$10^3$]{} & $45^\circ$ & $256\times72$ & $2\times10^3$ & $1\times10^7$ & $5\times10^7$\
n1e4 & $10^3$ & [$10^4$]{} & $45^\circ$ & $256\times72$ & $2\times10^3$ & $6\times10^6$ & $2\times10^7$\
n1e6 & $10^3$ & [$10^6$]{} & $45^\circ$ & $256\times72$ & $2\times10^3$ & $2\times10^6$ & $2\times10^6$\

\

NOTES.$^{a}$Disc radiation with shadowing effect is assumed except for Di and Ddn runs; $^{b}$isotropic radiation; $^{c}$disc radiation without shadowing effect; $^{d}$Dds run is also called s100, M1e3 and n1e5 runs.

We perform a set of simulations to see how the directional dependence of BH irradiation affects the nature of accretion. Table \[tab:model\] summarizes model parameters and numerical settings adopted for the cases examined. In all the cases, we initially set a static and homogeneous neutral medium with the number density $n_\infty$ and the temperature $T_{\mathrm{HI}}=10^4{\,\mathrm{K}}$ around a central BH. The BH mass $M_{\mathrm{BH}}$ is fixed constant during the calculation for simplicity.

In Section \[sec:structure\], we perform three high-resolution simulations, called “D-series” (for “Directional”), with different types of the directional dependence of the BH irradiation. For “Di run” (“i” for “isotropic”), we assume the isotropic irradiation, i.e., $f_{\mathrm{disc}} = f_{\mathrm{shadow}}=1$ in equation . The anisotropic disc radiation without the outer shadowing effect, i.e., $f_{\mathrm{shadow}}=1$, is assumed for “Ddn run” (“dn” for “disc no-shadow”), and both $\theta$-dependences of $f_{\mathrm{disc}}$ and $f_{\mathrm{shadow}}$ are allowed for “Dds run” (“ds” for “disc shadow”). Below we take $\theta_{\mathrm{shadow}} = 45^\circ$ as the fiducial value for the shadow opening angle. For the other parameters, we take $M_{\mathrm{BH}}=10^3M_\odot$ and $n_\infty=10^{5}{\,\mathrm{cm^{-3}}}$. Note that, for this set of $M_{\mathrm{BH}}$ and $n_\infty$, previous studies with isotropic BH irradiation have shown that the accretion rate is significantly reduced by radiation feedback ([@Milosavljevic:2009ab; @Park:2012aa; @Inayoshi:2016ac]; see also Sec. \[sec:sph\_acc\]).

In Sec. \[sec:dependence\], we study how the BH accretion changes with different shadow size $\theta_{\mathrm{shadow}}$, BH mass $M_{\mathrm{BH}}$, and ambient density $n_\infty$. First, to see the $\theta_{\mathrm{shadow}}$-dependence, we perform three simulations of “s-series” (for “shadow”) with different values of $\theta_{\mathrm{shadow}}$ (Sec. \[sec:sdep\]). Specifically, we take $\theta_{\mathrm{shadow}}=11.25^\circ$, $22.5^\circ$, $33.75^\circ$ and $45^\circ$. Second, we study the $M_{\mathrm{BH}}$-dependence with the “M-series”, where we take $M_{\mathrm{BH}}=10^2$, $10^3$, $10^4$ and $10^5\,M_\odot$ (Sec. \[sec:Mdep\]). Finally, the $n_\infty$-dependence is examined with the “n-series”, where we take different values of $n_\infty=10^3$, $10^4$, $10^5$ and $10^6{\,\mathrm{cm^{-3}}}$ (Sec. \[sec:ndep\]). In the above simulations, we take the fiducial values of $M_{\mathrm{BH}}=10^3M_\odot$, $n_\infty=10^{5}{\,\mathrm{cm^{-3}}}$ and $\theta_{\mathrm{shadow}} =45^\circ$ unless otherwise stated. We discuss the parameters relevant to the growth of the remnant BHs of Pop III stars in Section \[sec:conclusion\].

For each case, the inner and outer boundaries $R_{\mathrm{in}}$ and $R_{\mathrm{out}}$ are determined in the following way. We choose small enough $R_{\mathrm{in}}$ to correctly evaluate $\dot{M}$. To be more specific, $R_{\mathrm{in}}$ is taken to be much smaller than the Bondi radius for a neutral (ionized) gas when the dominant component of the accreting gas is neutral (ionized). We choose large enough $R_{\mathrm{out}}$ to keep an [H[ii]{} ]{}bubble within a simulation region. We only allow the flow going out of the computational domain at the inner boundary at $R_{\mathrm{in}}$, where $\dot{M}$ is evaluated. Across the outer boundary at $R_{\mathrm{out}}$, however, both the inflow and outflow are allowed. In the angular direction, the computational domain is $0<\theta<90^\circ$ under the assumption of the equatorial symmetry.

The grid numbers are taken to be $N_r\times N_\theta=512\times144$ and $256\times 72$ for the high- and medium-resolution simulations, respectively (see Table \[tab:model\]). In order to simultaneously resolve the Bondi and Strömgren radii, which are different typically by $3-4$ orders of magnitude, we increase the radial cell size $\Delta
r$ with the fixed size ratio $\Delta r_i/\Delta r_{i-1}$ $(>1)$. We set $\Delta r_1 = 0.1 R_{\mathrm{in}}$ at the inner boundary. The grids in the angular direction are homogeneously distributed over $0<\theta<90^\circ$, and thus the grid size is $\Delta\theta =
90^\circ/N_\theta$.

We have tested the convergence of the numerical results by varying the grid numbers or inner boundary radius (see Appendix \[sec:res\_check\]). We follow the evolution over the duration $t_{\mathrm{end}}$, until the accretion reaches a steady state in Di and Ddn runs, or until $\dot{M}$ reaches almost constant in the other runs.

Results {#sec:result}
=======

Structures of flows {#sec:structure}
-------------------

. \[tab:d-model\]

  run    subgrid radiation type   $\theta_{\mathrm{inflow}}(r_{\mathrm{B}})$$^{a}$   $\dot{M}/\dot{M}_{\mathrm{B}}$$^{d}$
  ----- ------------------------ -------------------------------------------------- --------------------------------------
  Di           isotropic                             ...$^{b}$                                   $0.17\%^{e}$
  Ddn             disc                               ...$^{c}$                                   $0.13\%^{e}$
  Dds       disc $+$ shadow                          $40^\circ$                                   $59\%^{f}$

  : Summary of the results in Sec. \[sec:structure\]

\

NOTES. $^{a}$opening angle of equatorial neutral inflow region at $r_{\mathrm{B}}$ (see text); $^{b}$no equatorial neutral region; $^{c}$equatorial neutral region does not reach $R_{\mathrm{in}}$; $^{d}$accretion rate normalized by Bondi one; $^{e}$averaged between $t=4\times 10^5{\,\mathrm{yr}}$ and $5\times 10^5{\,\mathrm{yr}}$; $^{f}$evaluated at the end of simulation.

In this section, we perform the simulations of “D-series”, in order to see how the flow structure changes with different directional dependences of the radiation fields. The basic results for these cases are summarized in Table \[tab:d-model\].

### Case with isotropic radiation {#sec:spherical_sym}

We first describe “Di run”, for which we assume the isotropic BH irradiation with $M_{\mathrm{BH}}=10^3M_\odot$ and $n_\infty=10^5{\,\mathrm{cm^{-3}}}$. The behaviour of the accretion flow is qualitatively the same as those obtained in the earlier works [@Milosavljevic:2009ab; @Park:2011aa; @Park:2012aa], and so we refer the readers to the above literature for full details.

Fig. \[fig:mdot\](a) shows the time evolution of the accretion rate and luminosity in our 2D simulation, along with the result of 1D calculation with the same parameter set. As explained by [@Park:2011aa], the accretion rate oscillates by repeating the following three phases: (a) high thermal pressure of the hot [H[ii]{} ]{}bubble suppresses the gas inflow and forms a dense shell of the swept up neutral gas; (b) the accretion rate, and hence the luminosity, decreases because the density of the [H[ii]{} ]{}bubble decreases due to the bubble expansion and/or gas inflow into the sink, leading to the contraction of the [H[ii]{} ]{}bubble with the dense shell; (c) an accretion burst caused by the collapse of the shell dramatically increases the luminosity and revive the large [H[ii]{} ]{}bubble again. The interval time between bursts roughly corresponds to the sound crossing time across the [H[ii]{} ]{}bubble. The accretion history in the 2D simulation is identical to the 1D result in an early stage ($t\lesssim 5\times10^4{\,\mathrm{yr}}$), but deviates from it later on because the spherical symmetry breaks down due to the growth of numerical perturbations by the instability of expanding ionization front [e.g., @Garcia-Segura:1996aa; @Whalen:2008ac; @Whalen:2008aa; @Park:2014ab]. Although the qualitative features are similar in the 1D and 2D cases, the accretion variability is slightly weaker in the latter case. The peaks of accretion burst in different directions are smoothed out because they are not exactly synchronized in the 2D case.

The average accretion rate between $t=4\times 10^5{\,\mathrm{yr}}$ and $5\times
10^5{\,\mathrm{yr}}$ is only $\dot{M}=1.7\times10^{-3}\dot{M}_{\mathrm{B}}$. Such a low rate is consistent with the Bondi rate in the ionized medium, $\dot{M}_{\mathrm{B,HII}} \sim 1 \times 10^{-3} \dot{M}_{\mathrm{B}}$ (see Sec. \[sec:sph\_acc\]). Our result is in good agreement with the previous ones by @Milosavljevic:2009ab and @Park:2012aa, who provided $\dot{M}\sim 2\times10^{-3}\dot{M}_{\mathrm{B}}$ and $\dot{M}\sim
10^{-2}\dot{M}_{\mathrm{B}}$, respectively. The differences of a factor of a few might come from differences in the adopted chemistry, because the accretion rate is sensitive to the thermal structure within the [H[ii]{} ]{}bubble [@Park:2011aa; @Park:2012aa].

Fig. \[fig:snap\_Di\](a) and (b) show the structures of accretion flows before and after an accretion burst, respectively. As explained above, the [H[ii]{} ]{}bubble shrinks before the burst (Fig. \[fig:snap\_Di\]a) and expands again due to the enhanced luminosity after the burst (Fig. \[fig:snap\_Di\]b). Whereas the Bondi radius $r_{\mathrm{B}}$ for the ambient neutral gas is illustrated in the figure, $r_{\mathrm{B,HII}}$ for the ionized medium is, although resolved in our simulations, too small to be shown. The velocity field of the ionized gas does not exhibit a systematic inflow but subsonic turbulent structure since the gas pressure dominates the gravity outside the Bondi radius. These snapshots are also very similar to those shown in the previous works [@Milosavljevic:2009ab; @Park:2011aa; @Park:2012aa].

Note that there appears a thin finger-like structure of the neutral gas along the $z$ axis in Fig. \[fig:snap\_Di\](b). Since both initial condition and BH irradiation are spherically symmetric, the flow patterns should also be spherically symmetric at least in a statistical sense. Thus, this is an artifact of our 2D simulation, presenting its limitation. Any flows toward the $z$ axis inevitably collide each other on the axis due to the assumed axisymmetry, creating a high density neutral gas column that shadows the cells behind it. We expect this artifact vanishes in future 3D simulations.

### Case with disc radiation without shadowing effect {#sec:disc_rad}

Next, we consider “Ddn run”, for which we assume the disc radiation without shadowing effect. Specifically, we adopt the anisotropy factor $\mathcal{F}=2\sin\theta$ by taking $f_{\mathrm{shadow}}=1$ in equation .

Fig. \[fig:mdot\](b) presents that, as in the case with isotropic radiation, the accretion rate and luminosity initially show strong oscillatory behaviours. However, the oscillation settles down in $3\times10^5{\,\mathrm{Myr}}$, after which only weak variability remains. We consider that the initial strong oscillation occurs due to the artificial initial condition of a static homogeneous medium. The accretion variability gradually ceases, after which the flow structure reaches a quasi-steady state. In this case, burst accretions coming from different directions cannot be synchronized due to the aspherical shape of the [H[ii]{} ]{}bubble created by the anisotropic BH irradiation, resulting in the less variable accretion rate. The mean accretion rate between $t=4\times 10^5{\,\mathrm{yr}}$ and $5\times 10^5{\,\mathrm{yr}}$, when the oscillation has already abated, is $\dot{M}=1.3\times10^{-3}\dot{M}_{\mathrm{B}}$. This low rate is nearly the same as that in the isotropic irradiation case, and also well approximated by the Bondi rate from the ionized medium (see Sec. \[sec:sph\_acc\]). The inflows through the equatorial neutral region have little contribution to the accretion rate, as will be seen below.

![Same as Fig. \[fig:snap\_Di\] but for the case with disc radiation without shadowing effect. The structure of accretion flow at the end of the simulation is shown on the larger and smaller scales in panels (a) and (b), respectively. The white dashed square in panel (a) represents the region plotted in (b). In addition to the Bondi radius for the ambient neutral gas ($r_{\mathrm{B}}$, dashed black) shown in both panels (a) and (b), that for the ionized medium ($r_{\mathrm{B,HII}}$, dashed white) is shown in (b).[]{data-label="fig:snap_Dd"}](figure/snap_Ddn.eps){width="7cm"}

Fig. \[fig:snap\_Dd\] shows the structure of the accretion flow at the end of the simulation. The whole [H[ii]{} ]{}bubble is shown in the upper panel, while the central region on the scale of $r_{\mathrm{B}}$ is enlarged in the lower panel. In Fig. \[fig:snap\_Dd\](a), we see that the pressure equilibrium is approximately realized throughout the simulation region. The [H[ii]{} ]{}bubble is squeezed in the equatorial directions because of the anisotropic irradiation. The gas within the [H[ii]{} ]{}bubble moves upward until colliding with the ambient neutral medium. A high-density region near the $z$ axis just outside the [H[ii]{} ]{}bubble is again likely to be an artifact, as seen in Sec. \[sec:spherical\_sym\]. Hereafter, we will ignore this kind of features since it hardly affects our conclusion. In Fig. \[fig:snap\_Dd\](b), the gas is ionized in most of the region except the equatorial thin neutral layer, which extends inward across the Bondi radius but does not reach the sink. Most of the outflow within the [H[ii]{} ]{}region is launched from this equatorial neutral layer. In this figure, we do not clearly see the inflow, which is actually limited to a very central part, because the Bondi radius for the ionized medium $r_{\mathrm{B,HII}}$ is currently much smaller than the size of the plotted area.

With the current spatial resolution, the thickness of the equatorial neutral layer is limited by the angular cell size of $\Delta \theta =
0.6^\circ$. The ionizing photon flux injected into the cells closest to the equatorial plane ($0 \leq \theta \leq \Delta\theta$) is reduced to $\sim 0.01$ of the angle-averaged value for the assumed anisotropy. This flux is, however, still large enough to make the [H[ii]{} ]{}region extend beyond the sink radius. In fact, the Strömgren radius with $n=
n_\infty =10^5{\,\mathrm{cm^{-3}}}$ and $L = 10^{-3}L_{\mathrm{E}}$ gives $r_{\mathrm{HII}}
\simeq 7\times10^3{\,\mathrm{AU}}$ (equation \[eq:13\]), which is much larger than the sink radius $R_{\mathrm{in}} = 3 \times 10^2{\,\mathrm{AU}}$. The horizontal extension of the bubble is even larger than this because the actual bubble density is smaller than $n_\infty$. If we could perform a simulation with much higher resolution, a thinner equatorial neutral region would reach the sink as the less ionizing photon flux is injected for the smaller $\theta$. Such a very thin neutral region, however, is not expected to affect the overall accretion, because the mass that can be carried through such a very thin region is severely limited and is further reduced by the mass loss into the [H[ii]{} ]{}bubble, as will be shown later (Sec. \[sec:shadow\_rad\] - \[sec:mass\_loss\_inner\]). Moreover, diffuse recombination photons processed within the bubble, which are not considered in the current simulations, would eliminate such a very thin neutral region [see, e.g., @Hollenbach:1994aa; @Tanaka:2013aa].

### Case with disc radiation with shadowing effect {#sec:shadow_rad}

Finally, we describe “Dds run”, in which the inner disc radiation is modified by the outer shadowing effect, as $\mathcal{F}=C\,f_{\mathrm{disc}}\,f_{\mathrm{shadow}}a$ (equation \[eq:5\]). We adopt $\theta_{\mathrm{shadow}}=45^\circ$ for the outer anisotropy factor $f_{\mathrm{shadow}}$ given by equation .

Fig. \[fig:mdot\](c) shows that the accretion rate converges towards an constant value $\dot{M} \simeq 0.59~\dot{M}_{\mathrm{B}}$ in $1{\,\mathrm{Myr}}$, which is much higher than the value $\dot{M} \lesssim
2\times10^{-3}~\dot{M}_{\mathrm{B}}$ obtained in the former two cases with isotropic radiation and disc radiation without shadowing effect. This accretion rate is also “super-critical” and 400 times larger than the Eddington-limited rate $\dot{M}_{\mathrm{E}}$. By the end of the simulation, the luminosity also converges to $L \simeq 8L_{\mathrm{E}}$, i.e., a super-Eddington luminosity realized by the high accretion rate.

Fig. \[fig:snap\_Ds\] shows the structure of the accretion flow at the end of the simulation. The whole [H[ii]{} ]{}bubble is shown within the large plotted area of Fig. \[fig:snap\_Ds\](a), while the central regions over $\sim 10^6$ AU and $\sim 10^5$ AU scales are enlarged in Figs. \[fig:snap\_Ds\](b) and \[fig:snap\_Ds\](c). Owing to the shadowing effect, we see in Fig. \[fig:snap\_Ds\](a) that the large horizontal neutral region cuts into the central part with the [H[ii]{} ]{}bubbles bound to the bipolar regions. Similarly to the case with disc radiation without shadowing effect (Sec. \[sec:disc\_rad\]), the gas within the bipolar [H[ii]{} ]{}bubbles flows outward and collides with the ambient neutral medium. The pressure equilibrium is approximately achieved throughout the simulation region, although the thermal pressure slightly decreases before the collision because the total pressure including the ram pressure is balanced. The size of the bipolar [H[ii]{} ]{}bubbles is much larger than in the former cases owing to the much higher luminosity $L \simeq 8L_{\mathrm{E}}$. Note that the injected radiation is super-Eddington only in the polar region with $\theta \gtrsim 50^\circ$.

Fig. \[fig:snap\_Ds\](b) shows that the gas inside the equatorial neutral region is almost at rest in the pressure equilibrium, while that on the surfaces is photoevaporated to join the ionized outflow. In Fig. \[fig:snap\_Ds\](c), where the structure over the scale of the Bondi radius is presented, the gas flows into the central sink through the equatorial neutral region. The density and pressure increase with decreasing $r$ for $r \lesssim r_{\mathrm{B}}$, as expected for the Bondi flow. As seen in Fig. \[fig:snap\_Ds\](b), the neutral gas is photoevaporated into the [H[ii]{} ]{}regions, where the acceleration by the radiation pressure is stronger than the gravitational pull owing to the super-Eddington fluxes in the polar directions. The accretion proceeds only through the solid angle covered by the equatorial neutral region. Note that, in our simulations, we neglect possible photoevaporation outflow coming out from the sink. We discuss it later in Sec. \[sec:mass\_loss\_inner\]. In the next section, we investigate the structure of the flow in more detail.

### Analysis of flow structure in case with shadowing effect {#sec:anl_modeling}

In this section, we develop an analytical model and compare it with our result to examine the inflow-outflow structure presented in Fig. \[fig:snap\_Ds\]. The overall structure of our model is schematically depicted in Fig. \[fig:anl\_model\] and can be summarized as follows: in the equatorial neutral region where ionizing photons cannot penetrate, the gas inflows in a Bondi accretion fashion; in the bipolar [H[ii]{} ]{}regions where ionizing photons heat up the gas via photoionization, the outflows are launched due to the thermal and radiation pressure; through their boundaries, the photoevaporating gas is lost from the neutral region and supplied into the [H[ii]{} ]{}regions.

We begin with considering the density profiles of the inflow and outflow. The radial density profile in the equatorial neutral region $n_{\mathrm{inflow}}(r)$ is well approximated by that of the Bondi solution, which we further simplify as $$\begin{aligned}
  n_{\mathrm{inflow}}(r) &=
\begin{cases}
{\displaystyle}n_\infty\left(\frac{r}{r_{\mathrm{B}}}\right)^{-3/2}&r < r_{\mathrm{B}}\\[0.4cm]
n_\infty&r > r_{\mathrm{B}}
\end{cases}\label{eq:9} \,.
 \end{aligned}$$ This expression slightly over- and underestimates the obtained density profile at $r\ll r_{\mathrm{B}}$ and $r\sim r_{\mathrm{B}}$. As for the bipolar ionized outflows, the density profile $n_{\mathrm{outflow}}(r)$ can be estimated by assuming the pressure equilibrium at the conical boundaries between the neutral and ionized gas, $$\begin{aligned}
n_{\mathrm{outflow}}(r)&=n_{\mathrm{inflow}}(r)\left(\frac{T_{\mathrm{HI}}}{2T_{\mathrm{HII}}}\right)\,,
\label{eq:7}
 \end{aligned}$$ where $T_{\mathrm{HI}} \simeq 10^4{\,\mathrm{K}}$ and $T_{\mathrm{HII}} \simeq
7\times10^4{\,\mathrm{K}}$ are the temperature in the neutral and ionized regions, respectively. Although not very precise, this simple expression captures the qualitative features of the bipolar ionized outflows.

![The opening angle of the equatorial neutral inflow region $\theta_{\mathrm{inflow}}$ as a function of $r$ (solid red), along with the angle-dependence of the radius of the [[H[ii]{} ]{}]{} region $r=r_{\mathrm{HII}}(\theta)$ (dashed orange; equation \[eq:6\]). The dotted part of the red line marks a few innermost cells that are artificially ionized as we neglect the absorption within the sink (also see the text). []{data-label="fig:th_in_Ds"}](figure/th_in_Ds.eps){width="8.5cm"}

Fig. \[fig:th\_in\_Ds\] shows the opening angle of the equatorial neutral region $\theta_{\mathrm{inflow}}$ as a function of $r$. In practice, we define the neutral region as the region where the ionization degree of hydrogen is less than $50\%$. Although $\theta_{\mathrm{inflow}}$ decreases as $r$ decreases, the neutral inflow region reaches the inner boundary at $r = R_{\mathrm{in}}$ with a finite angle, unlike in the case with disc radiation without shadowing effect. The sharp drop of $\theta_{\mathrm{inflow}}$ around $R_{\mathrm{in}} = 2000$ AU is caused by photoionization of a few innermost cells due to our ignorance of the consumption of ionizing photons within the sink, although this does not affect our conclusion. The opening angle at the Bondi radius, $\theta_{\mathrm{inflow}}(r_{\mathrm{B}})=40^\circ$, is similar to the assumed shadow opening angle $\theta_{\mathrm{shadow}} =45^\circ$ (equation \[eq:10\]).

In order to estimate $\theta_{\mathrm{inflow}}$, we calculate the radius of the [[H[ii]{} ]{}]{} region $r_{\mathrm{HII}}$ in each direction $\theta$ with the modeled density profiles. The supply rate of ionizing photons per unit solid angle is given by $\dot{N}_{\mathrm{ion}}\,\mathcal{F}(\theta)/4\pi$, where $\dot{N}_{\mathrm{ion}}$ is the total ionizing photon emissivity of the central accretion disc, given by $\dot{N}_{\mathrm{ion}}=L/3h\nu_{\mathrm{T}}$ for the assumed spectral shape of $L_\nu\propto \nu^{-1.5}$. Equating this supply rate with the recombination rate, we have $$\frac{L}{3h \nu_{\mathrm{T}}}\frac{\mathcal{F}(\theta)}{4 \pi}
= \int_{R_{\mathrm{in}}}^{r_{\mathrm{HII}}}(\theta)
\alpha_{\mathrm{B}}\,n_{\mathrm{outflow}}^2\,r^2\,{\mathrm{d}}r ,
\label{eq:uve}$$ where $n_{\mathrm{outflow}}$ is given by equation . Performing the integration in equation , we finally get $$\begin{aligned}
 r_{\mathrm{HII}}(\theta) 
=
\begin{cases}
{\displaystyle}R_{\mathrm{in}}  \exp\left[
\frac{L\,\mathcal{F}(\theta)\,T_{\mathrm{HII}}^2}{3\pi h\nu_{\mathrm{T}}\, \alpha_{\mathrm{B}}\, n_\infty^2\, r_{\mathrm{B}}^3\, T_{\mathrm{HI}}^2}
\right]&r_{\mathrm{HII}} < r_{\mathrm{B}}\\[0.4cm]
{\displaystyle}\left[
\frac{L\,\mathcal{F}(\theta)\,T_{\mathrm{HII}}^2}
{\pi h\nu_{\mathrm{T}}\,\alpha_{\mathrm{B}}\,n_\infty^2\,T_{\mathrm{HI}}^2} 
- 3r_{\mathrm{B}}^3 \ln\left(\frac{r_{\mathrm{B}}}{R_{\mathrm{in}}}\right)
+ r_{\mathrm{B}}^3
\right]^{1/3} \hspace{-2cm}\\[0.4cm]
&r_{\mathrm{HII}} > r_{\mathrm{B}}
 \end{cases}\,.
\label{eq:6}\end{aligned}$$ Here, we show the relation $r=r_{\mathrm{HII}}(\theta)$, or equivalently $\theta=r_{\mathrm{HII}}^{-1}(r)$, in Fig. \[fig:th\_in\_Ds\] with $L =
8.2\,L_{\mathrm{E}}$ (Fig. \[fig:mdot\]). We see that equation  qualitatively reproduces $\theta_{\mathrm{inflow}} (r)$ obtained in the simulation with a small deviation of a few degrees at each $r$. Such a deviation mainly comes from approximate modelling of $n_{\mathrm{outflow}}$ in equation . For example, in an outer part of the [H[ii]{} ]{}bubble where helium is not doubly ionized, the gas is no longer heated up to $7\times10^4{\,\mathrm{K}}$ by the ${\mathrm{He^+}}$ photoionization, resulting in the higher density than that estimated by equation  with $T_{\mathrm{HII}} = 7 \times 10^4{\,\mathrm{K}}$. Nonetheless, our simple modelling with equation  describes the numerical results well.

![The radial dependence of the equatorial inflow rate $\dot{M}_{\mathrm{inflow}}(r)$ (solid blue) and bipolar outflow rate $\dot{M}_{\mathrm{outflow}}(r)$ (solid red), defined in equation . The blue arrow marks the estimated value of the inflow rate at the inner boundary $\dot{M}_{\mathrm{inflow}}(R_{\mathrm{in}})$ (equation \[eq:17\]). The red dashed line represents the analytical estimate of the outflow rate $\dot{M}_{\mathrm{outflow}}(r)$ (equation \[eq:14\]). The values of the Bondi radius $r_{\mathrm{B}}$ and accretion rate $\dot{M}_{\mathrm{B}}$ are marked by the black arrows. []{data-label="fig:mdoio_Ds"}](figure/mdotio_Ds.eps){width="8.5cm"}

Finally, we investigate the flow rates in the neutral and ionized regions. Fig. \[fig:mdoio\_Ds\] shows the equatorial inflow rate $\dot{M}_{\mathrm{inflow}}$ and bipolar outflow rate $\dot{M}_{\mathrm{outflow}}$ through a spherical surface with radius $r$, $$\begin{aligned}
 \dot{M}_{\mathrm{inflow}}(r) &=  - 4\pi r^2\int_0^{\theta_{\mathrm{inflow}}}  \rho\,v_r\,\cos \theta d\theta\,,{\nonumber\\}\dot{M}_{\mathrm{outflow}}(r) &=  4\pi r^2\int_{\theta_{\mathrm{inflow}}}^{\pi/2} \rho\,v_r\, \cos \theta d\theta\,,
\label{eq:16}\end{aligned}$$ where $v_r$ is the outward velocity. We have multiplied a factor of two to take into account the equatorial symmetry. The net accretion rate $\dot{M} \equiv \dot{M}_{\mathrm{inflow}}-\dot{M}_{\mathrm{outflow}}$ is almost constant with $r$, consistent with a quasi-steady flow structure. The value of $\dot{M}$ is equal to the inflow rate $\dot{M}_{\mathrm{inflow}}$ at the inner boundary $R_{\mathrm{in}}$, where $\dot{M}_{\mathrm{outflow}} = 0$ is imposed as the boundary condition.

We now estimate $\dot{M}_{\mathrm{inflow}}(R_{\mathrm{in}})$ ($=\dot{M}$) from the Bondi-like accretion through a solid angle corresponding to the opening angle $\theta_{\mathrm{inflow}}(r_{\mathrm{B}})$. With the solid angle $\Delta\Omega_{\mathrm{inflow}}(r_{\mathrm{B}})=4\pi
\,\sin\theta_{\mathrm{inflow}}(r_{\mathrm{B}})$, we obtain $$\begin{aligned}
 \dot{M}_{\mathrm{inflow}}(R_{\mathrm{in}})=
\frac{\Delta\Omega_{\mathrm{inflow}}(r_{\mathrm{B}})}{4\pi}\dot{M}_{\mathrm{B}}\,.
\label{eq:17}\end{aligned}$$ Fig. \[fig:mdoio\_Ds\] demonstrates that equation  well reproduces $\dot{M}_{\mathrm{inflow}}(R_{\mathrm{in}})$ measured in the simulation,[^4] although slightly overestimated owing to the photoevaporation mass loss. Note that the inflow rate $\dot{M}_{\mathrm{inflow}}$ increases with $r$ and exceeds the Bondi rate $\dot{M}_{\mathrm{B}}$ at $\sim 3000$ AU. This is because the circulation flows are generated to compensate the photoevaporation mass loss from the equatorial neutral region (see below) and contribute to the accretion rate together with the Bondi-like inflow.

The outflow rate can be modeled as the photoevaporation mass loss. The mass-loss flux from the surfaces of neutral region is estimated as $f_{\mathrm{outflow}}\,\rho_{\mathrm{b,HII}}\, c_{\mathrm{s,HII}}$, where $f_{\mathrm{outflow}}$ is an $O(1)$ correction factor, $\rho_{\mathrm{b,HII}}$ the density at the bottom of the ionized layer and $c_{\mathrm{s,HII}}$ the sound velocity for ionized gas [e.g., @Hollenbach:1994aa; @Tanaka:2013aa]. The outflow rate $\dot{M}_{\mathrm{outflow}}$ through a given radius $r$ is obtained by integrating the mass-loss fluxes between $R_{\mathrm{in}}$ and $r$ because the outflow is in a quasi-steady state. In reality, additional mass loss may happen even inside $R_{\mathrm{in}}$, as we will discuss in Sec. \[sec:mass\_loss\_inner\]. With $\rho_{\mathrm{b,HII}} \simeq m_{\mathrm{p}}
n_{\mathrm{outflow}}$ (equation \[eq:7\]) and $c_{\mathrm{s,HII}}=(2T_{\mathrm{HII}}/T_{\mathrm{HI}})^{1/2}c_{\mathrm{s,HI}}$, we obtain $$\begin{aligned}
 \dot{M}_{\mathrm{outflow}}(r)
\simeq 4\pi\int_{R_{\mathrm{in}}}^r &f_{\mathrm{outflow}}\,m_{\mathrm{p}}\,c_{\mathrm{s,HII}}\,n_{\mathrm{outflow}} r' {\mathrm{d}}r'{\nonumber\\}=
2\pi\,f_{\mathrm{outflow}}\,& m_{\mathrm{p}}\, n_\infty\,c_{\mathrm{s,HI}}
\left(\frac{2T_{\mathrm{HII}}}{T_{\mathrm{HI}}}\right)^{1/2}{\nonumber\\}\times&
\begin{cases}
4 r_{\mathrm{B}}^{3/2}  \left(r^{1/2}-R_{\mathrm{in}}^{1/2}\right)
& r < r_{\mathrm{B}} \\
r^2 +  3r_{\mathrm{B}}^2 -4r_{\mathrm{B}}^{3/2}R_{\mathrm{in}}^{1/2}
&
r > r_{\mathrm{B}}
\end{cases}\,,
\label{eq:14}\end{aligned}$$ where a factor of two is multiplied in the first equality to take into account both top and bottom surfaces. We see in Fig. \[fig:mdoio\_Ds\] that modeled $\dot{M}_{\mathrm{outflow}}$ with the best-fit value of $f_{\mathrm{outflow}}=0.7$ reproduces the simulation result with remarkable agreement.

### Possible mass loss from neutral inflow inside the sink {#sec:mass_loss_inner}

As mentioned above, we neglect the possible mass loss from the innermost part of the flow masked by the sink. Since the size of the accretion disc is supposed to be much smaller than the sink radius $R_{\mathrm{in}}$ (see Fig. \[fig:acc\_whole\]), we neglect the centrifugal effect and assume the similar flow structure extends inward. As an upper limit for the mass-loss rate, we evaluate the integral with the same integrand as equation  but for the different range of $0< r <
R_{\mathrm{in}}$, and obtain $$\begin{aligned}
 \dot{M}_{\mathrm{loss}}(<R_{\mathrm{in}})
\lesssim 0.1 \left[\frac{R_{\mathrm{in}}}{(r_{\mathrm{B}}/7)}\right]^{1/2} \dot{M}_{\mathrm{B}}\,.
\label{eq:4}\end{aligned}$$ Here, we take $R_{\mathrm{in}} \approx r_{\mathrm{B}}/7$ of our simulation setup (see Table \[tab:model\]) as a reference value.

The mass supply rate to the accretion disc can be conservatively estimated by $\dot{M} - \dot{M}_{\mathrm{loss}}(<R_{\mathrm{in}})$, meaning that $\dot{M}$ measured in the simulation slightly overestimates the true value. This would be alleviated by taking a smaller value for $R_{\mathrm{in}}$. We have performed a test run with the smaller sink radius (see Appendix \[sec:res\_check\]), but found no remarkable differences of the accretion rate. Note that the effect of the angular momentum becomes important on the smaller scale. It is not allowed to take an arbitrary small sink radius without considering the effect of the angular momentum.

If the accretion disc is spatially resolved, we expect further mass loss happens due to, e.g., the disc winds [e.g., @Blandford:1999aa; @Zahra-Zeraatgari:2016aa; @Begelman:2016aa] and/or jets from a close vicinity of the BH [e.g., @Ohsuga:2005aa; @Jiang:2014aa; @Yuan:2015aa; @Sc-adowski:2016aa]. Outflows from the sink caused by such phenomena may change the outer gas dynamics on the scale of the Bondi radius. This should be studied in future work.

Parameter dependence {#sec:dependence}
--------------------

Here, we study how the flow structure changes with variation of the simulation parameters: the shadow opening angle $\theta_{\mathrm{shadow}}$ (in Sec. \[sec:sdep\]), BH mass $M_{\mathrm{BH}}$ (in Sec. \[sec:Mdep\]), and ambient density $n_\infty$ (in Sec. \[sec:ndep\]). In Sec. \[sec:comp\_prev\], we compare our results with previous 1D calculations.

### Dependence on shadow size {#sec:sdep}

  run           $\theta_{\mathrm{shadow}}$$^{a}$   $\theta_{\mathrm{inflow}}(r_{\mathrm{B}})$$^{b}$   $\dot{M}/\dot{M}_{\mathrm{B}}$$^{c}$
  ------------ ---------------------------------- -------------------------------------------------- --------------------------------------
  s100 (Dds)               $45^\circ$                                 $40^\circ$                                     $59\%$
  s075                   $37.75^\circ$                                $29^\circ$                                     $42\%$
  s050                     $25^\circ$                                 $19^\circ$                                     $25\%$
  s025                   $11.25^\circ$                                $9^\circ$                                     $6.5\%$

  : Summary of the $\theta_{\mathrm{shadow}}$ dependence.[]{data-label="tab:s-model"}

\

NOTES. $^{a}$shadow opening angle of our subgrid model (equation \[eq:10\]); $^{b}$opening angle of equatorial neutral inflow region at $r_{\mathrm{B}}$; $^{c}$accretion rate normalized by Bondi one.

Considering uncertainties in the anisotropic shadowing effect (Sec. \[sec:direction\_dep\]), we study the cases with different shadow opening angles $\theta_{\mathrm{shadow}}$ by reducing it from $45^\circ$ in “Dds run” (here we also call it “s100 run”) to $33.75^\circ$ (“s075 run”), $22.5^\circ$ (“s050 run”) and $11.25^\circ$ (“s025 run”). We call this series of runs as “s-series”. We take $M_{\mathrm{BH}}=10^3M_\odot$ and $n_\infty=10^5{\,\mathrm{cm^{-3}}}$ for the s-series.

Our main findings are as follows: in all the runs of the s-series, the overall flow structures are similar and the accretion rates $\dot{M}$ are much higher than in the cases without the shadow (i.e., Di and Ddn runs). The obtained accretion rates and opening angles of the equatorial neutral region at Bondi radius $\theta_{\mathrm{inflow}}(r_{\mathrm{B}})$ are summarized in Table \[tab:s-model\]. The values of $\theta_{\mathrm{inflow}}(r_{\mathrm{B}})$ agree well with the prediction by equation  with only small offsets $\lesssim 3^\circ$. Note also that $\theta_{\mathrm{inflow}}(r_{\mathrm{B}}) \simeq \theta_{\mathrm{shadow}}$ despite the gradual transition between the shadowed and non-shadowed regions modeled as in equation .

The equatorial inflow rates $\dot{M}_{\mathrm{inflow}} (r)$ (see equation \[eq:16\]) are shown in Fig. \[fig:mdotio\_sdep\](a). The values of $\dot{M}_{\mathrm{inflow}}$ at $R_{\mathrm{in}}$ agree well with the rates estimated by the Bondi flow through the solid angle of $4\pi\sin\theta_{\mathrm{inflow}}(r_{\mathrm{B}})$ (equation \[eq:17\]; arrows in Fig. \[fig:mdotio\_sdep\] a), but with slight downward offset due to the photoevaporation mass loss. Fig. \[fig:mdotio\_sdep\](b) shows the outflow rates in polar directions $\dot{M}_{\mathrm{outflow}}(r)$ (again, see equation \[eq:16\]). As seen in Sec. \[sec:shadow\_rad\], the estimate by equation  with $f_{\mathrm{outflow}}=0.7$ gives a good fit to the numerical results. Small differences of $\sim$ a few $\times$ 10 % among them are comparable to the intrinsic fluctuations of the outflow rates present in the quasi-steady states.

![The net accretion rate $\dot{M}$ normalized by the Bondi rate $\dot{M}_{\mathrm{B}}$ as a function of the opening angle of the horizontal neutral layer at the Bondi radius $\theta_{\mathrm{inflow}}(r_{\mathrm{B}})$. The crosses show the runs with shadow opening angles $\theta_{\mathrm{shadow}}=45^\circ$, $33.75^\circ$, $22.5^\circ$ and $11.25^\circ$. The solid line represents the relation given by equation  assuming $\dot{M}_{\mathrm{loss}}=0.07\dot{M}_{\mathrm{B}}$. []{data-label="fig:mdot_th_indep"}](figure/mdot_th_indep.eps){width="8.5cm"}

The net accretion rates $\dot{M}$ are plotted as crosses in Fig. \[fig:mdot\_th\_indep\] against $\theta_{\mathrm{inflow}}(r_{\mathrm{B}})$. Since equation  slightly overestimates $\dot{M}_{\mathrm{inflow}}(R_{\mathrm{in}})$ due to the photoevaporation mass loss, we modify equation  assuming a constant mass loss rate $\dot{M}_{\mathrm{loss}}$ in all cases, as $$\begin{aligned}
 \dot{M} = \dot{M}_{\mathrm{inflow}}(R_{\mathrm{in}}) =
\frac{\Delta\Omega_{\mathrm{inflow}}(r_{\mathrm{B}})}{4\pi}\dot{M}_{\mathrm{B}}
- \dot{M}_{\mathrm{loss}} \,.
\label{eq:18}\end{aligned}$$ We find that $\dot{M}_{\mathrm{loss}}=0.07\,\dot{M}_{\mathrm{B}}$ gives the best fit to the simulated results with errors less than 2% of $\dot{M}_{\mathrm{B}}$. This good agreement also supports the above assumption of constant $\dot{M}_{\mathrm{loss}}$. The value of $\dot{M}_{\mathrm{loss}}$ is similar to but smaller than $\dot{M}_{\mathrm{outflow}}(r_{\mathrm{B}})$ (equation \[eq:14\]) partly due to the contribution from the circulation flows, as mentioned in Sec. \[sec:anl\_modeling\]. Moreover, by setting $\dot{M} = 0$ in equation , we get the critical opening angle $\theta_{\mathrm{cr}} \simeq 4^\circ$, below which the equatorial neutral flow disappears by photoevaporation. This value will be raised up to $\theta_{\mathrm{cr}} \simeq 10^\circ$ if we include the mass loss inside the sink, which is currently ignored (see Sec. \[sec:mass\_loss\_inner\]).

### Dependence on BH mass {#sec:Mdep}

  run           $M_{\mathrm{BH}}\,[M_\odot]$   $\theta_{\mathrm{inflow}}(r_{\mathrm{B}})$   $\dot{M}/\dot{M}_{\mathrm{B}}$
  ------------ ------------------------------ -------------------------------------------- --------------------------------
  M1e2                     $10^2$                              $38^\circ$                               $55\%$
  M1e3 (Dds)               $10^3$                              $40^\circ$                               $59\%$
  M1e4                     $10^4$                              $43^\circ$                               $61\%$
  M1e5                     $10^5$                              $48^\circ$                               $67\%$

  : Summary of the $M_{\mathrm{BH}}$ dependence.[]{data-label="tab:M-model"}

\

Next, we study the dependence on the BH mass $M_{\mathrm{BH}}$ by performing a set of simulations termed “M-series”, where $M_{\mathrm{BH}}=10^2\,M_\odot$ (“M1e2 run”), $10^3\,M_\odot$ (“M1e3 run” identical to “Dds run”), $10^4\,M_\odot$ (“M1e4 run” ) and $10^5\,M_\odot$ (“M1e5 run”). The other parameters are set to $n_\infty=10^5{\,\mathrm{cm^{-3}}}$ and $\theta_{\mathrm{shadow}}=45^\circ$.

We find that the accretion proceeds roughly at the Bondi rate in a quasi-steady fashion for all the runs. Flow properties at the end of calculation are summarized in Table \[tab:M-model\]. We see that, for all the runs, the neutral region spans the opening angle $\theta_{\mathrm{inflow}}(r_{\mathrm{B}})\simeq \theta_{\mathrm{shadow}} (=45^\circ)$ at the Bondi radius, and that the Bondi-like accretion proceeds through this solid angle with the rate $\dot{M} \simeq 0.5-0.7\,
\dot{M}_{\mathrm{B}}$.

Note, however, that the opening angle and thus the accretion rate increase gradually with the BH mass. This dependence can be understood as follows. Recall that the luminosity $L$ is approximately proportional to the Eddington value or the mass $M_{\mathrm{BH}}$ in the super-Eddington regime (equations \[eq:12\] and \[eq:3\]). The radius of the [H[ii]{} ]{}region thus varies as $r_{\mathrm{HII}}(\theta) \propto M_{\mathrm{BH}}^{2/3}$ (equation \[eq:13\]), while the Bondi radius follows $r_{\mathrm{B}} \propto
M_{\mathrm{BH}}$ (equation \[eq:2\]). Since $r_{\mathrm{HII}}(\theta)$ is an increasing function of $\theta$ (equation \[eq:6\]; see also Fig. \[fig:th\_in\_Ds\]), this means that $\theta_{\mathrm{inflow}}(r_{\mathrm{B}})$, obtained by solving $r_{\mathrm{B}}=r_{\mathrm{HII}}(\theta)$ with respect to $\theta$, increases with $M_{\mathrm{BH}}$. Equation  indeed explains the variation of the opening angle $\theta_{\mathrm{inflow}}(r_{\mathrm{B}})$ in the numerical results within the error of $4^\circ$. Similarly, the accretion rates $\dot{M}/\dot{M}_{\mathrm{B}}$ estimated by equation  reproduce the results with errors $\lesssim10\%$.

### Dependence on ambient density {#sec:ndep}

  run           $n_\infty\,[{\mathrm{cm^{-3}}}]$   $\theta_{\mathrm{inflow}}(r_{\mathrm{B}})$   $\dot{M}/\dot{M}_{\mathrm{B}}$
  ------------ ---------------------------------- -------------------------------------------- --------------------------------
  n1e3                       $10^3$                                $36^\circ$                               $59\%$
  n1e4                       $10^4$                                $38^\circ$                               $54\%$
  n1e5 (Dds)                 $10^5$                                $40^\circ$                               $59\%$
  n1e6                       $10^6$                                $44^\circ$                               $71\%$

  : Summary of the $n_\infty$ dependence.[]{data-label="tab:n-model"}

\

Motivated by a wide variety of the environment in the vicinity of BHs, we finally investigate the cases with different ambient densities, termed “n-series”, where $n_\infty$ is $10^3{\,\mathrm{cm^{-3}}}$ (“n1e3 run”), $10^4{\,\mathrm{cm^{-3}}}$ (“n1e4 run”), $10^5{\,\mathrm{cm^{-3}}}$ (“n1e5 run” identical to “Dds run”) and $10^6{\,\mathrm{cm^{-3}}}$ (“n1e6 run”). The other parameters are set to $M_{\mathrm{BH}}=10^3\,M_\odot$ and $\theta_{\mathrm{shadow}}=45^\circ$.

The flow characteristics are similar regardless of $n_\infty$ with the neutral-region opening angle $\theta_{\mathrm{inflow}}(r_{\mathrm{B}})\simeq
\theta_{\mathrm{shadow}} (=45^\circ)$ and the mass accretion rates comparable to the Bondi rates, $\dot{M}/\dot{M}_{\mathrm{B}} \simeq 0.6-0.7$, in all the cases (Table \[tab:n-model\]). The increasing trend of $\theta_{\mathrm{inflow}}(r_{\mathrm{B}})$ with $n_\infty$ can be understood again as in Sec. \[sec:Mdep\]. Now the luminosity $L$ is almost independent of $n_\infty$ in the super-Eddington regime (equation \[eq:3\]), which leads to $r_{\mathrm{HII}}(\theta) \propto n_\infty^{-2/3}$ (equation \[eq:6\]) while $r_{\mathrm{B}}$ is independent of $n_\infty$ (equation \[eq:2\]). Since $r_{\mathrm{HII}}(\theta)$ is an increasing function of $\theta$, it follows that $\theta_{\mathrm{inflow}}(r_{\mathrm{B}})$ increases with $n_\infty$. Again, the analytic estimates of $\theta_{\mathrm{inflow}}(r_{\mathrm{B}})$ with equation  and $\dot{M}/\dot{M}_{\mathrm{B}}$ with equation  agree well with the numerical results, with errors less than $4^\circ$ and $10\%$, respectively.

Comparison with previous works {#sec:comp_prev}
------------------------------

Recently, @Inayoshi:2016ac and @Sakurai:2016aa have investigated necessary conditions to overcome the radiative feedback to achieve the Bondi-like accretion using 1D calculations under the spherical symmetry. These authors simulated accretion on to the BH in the same setting but with different prescriptions on the BH irradiation, i.e., whether or not the BH luminosity is capped by $L_{\mathrm{E}}$. They concluded that the efficient Bondi-like accretion appears when the following condition is satisfied: $M_{\mathrm{BH}}\,n_\infty \gtrsim 10^9\,
M_\odot\,{\mathrm{cm^{-3}}}$ (see Sec. \[sec:sph\_acc\]).

However, our 2D simulations suggest that the above criterion needs to be modified. We find the efficient accretion at about the Bondi rate is always possible, as long as the shadow opening angle has a certain size ($\theta_{\mathrm{shadow}} \gtrsim O(10)^\circ$), because the inflows from equatorial shadowed regions are allowed in 2D simulations. We here emphasize that neither $M_{\mathrm{BH}}$ nor $n_\infty$ appears in this condition, and that the efficient accretion is possible even with $M_{\mathrm{BH}}\,n_\infty \ll 10^9\, M_\odot\,{\mathrm{cm^{-3}}}$. The accretion rate through the shadowed direction is set by the Bondi rate, which depends on $M_{\mathrm{BH}}$ and $n_\infty$. For the accretion rate to largely exceeds the Eddington rate, $\dot{M}_{\mathrm{B}}/\dot{M}_{\mathrm{E}}=
\left(M_{\mathrm{BH}}\,n_\infty/10^5\,M_\odot\,{\mathrm{cm^{-3}}}\right) \gg 1$ should be met in addition to the above condition for the shadow size. Note, however, this condition for $M_{\mathrm{BH}}\,n_\infty$ is much easier to be satisfied than the condition for the 1D calculations, i.e., $M_{\mathrm{BH}}\,n_\infty \gtrsim 10^9\, M_\odot\,{\mathrm{cm^{-3}}}$.

Our 2D simulations show that the flow structure qualitatively differs from that in 1D even in the cases with $M_{\mathrm{BH}}\,n_\infty \gtrsim
10^9\, M_\odot\,{\mathrm{cm^{-3}}}$. For example, with $M_{\mathrm{BH}}=10^5\,M_\odot$ and $n_\infty=10^5{\,\mathrm{cm^{-3}}}$ (M1e5 run), the large bipolar [H[ii]{} ]{}bubbles persist in a steady state, whereas the spherical Bondi-like accretion quenches the [H[ii]{} ]{}bubble in the 1D test run. This is partly because the enhanced ionizing radiation in the polar directions due to the assumed directional dependence increases the size of the [H[ii]{} ]{}bubbles.

As seen above, $M_{\mathrm{BH}}\,n_\infty$ is not the key parameter to demarcate the regimes for efficient/inefficient accretion in 2D simulations, unlike in the 1D cases. The efficient accretion is possible if only the shadow size is sufficiently large, irrespective of $M_{\mathrm{BH}}$ or $n_\infty$.

conclusions and discussion {#sec:conclusion}
==========================

We have studied the black hole (BH) accretion of the primordial gas under anisotropic irradiation by the circum-BH accretion disc. Using two-dimensional radiation hydrodynamics simulations, we have solved the dynamics of the accretion flow spatially resolving both Bondi radius and the size of the [H[ii]{} ]{}region, which can differ by 3-4 orders of magnitude. We do not resolve the central accretion disc which emits anisotropic radiation, but inject ionizing photons at the inner boundary of the computational domain according to the subgrid prescription. To see how the anisotropy of the BH irradiation affects the flow structure, we first perform simulations with the three different types of the directional dependence: isotropic radiation, anisotropic radiation from the disc with and without the shadowing effect. For the case with the anisotropic shadowing effect, we have also studied the dependence of the flow structure on the shadow opening angle $\theta_{\mathrm{shadow}}$, BH mass $M_{\mathrm{BH}}$, and ambient density $n_\infty$.

With the isotropic irradiation, the accretion rate varies periodically as a result of recurrent formation and collapse of a hot and low-density [H[ii]{} ]{}bubble around the BH. The time-averaged accretion rate is only 0.2% of the original Bondi rate $\dot{M}_{\mathrm{B}}$ for the neutral medium and roughly given by the Bondi rate for the ionized medium [e.g., @Park:2011aa; @Park:2012aa; @Milosavljevic:2009ab]. Even with the anisotropy of the BH irradiation, the accretion rate is still similar to that in the isotropic case unless the shadowing effect is included. The flow structure in this case, however, is qualitatively different from the isotropic case: the large periodic variation, which has been reported in previous studies, disappears.

Unlike in the former two cases, the accretion rate becomes much higher in the case with the shadowing effect. For example, in the case with $M_{\mathrm{BH}}=10^3\,M_\odot$, $n_\infty=10^5{\,\mathrm{cm^{-3}}}$ and $\theta_{\mathrm{shadow}}=45^\circ$, the accretion rate reaches as high as 60% of the Bondi rate $\dot{M}_{\mathrm{B}}$ and is “super-critical” with 400 times larger than the Eddington-limited rate $\dot{M}_{\mathrm{E}}$. The flow structure in the steady state consists of the equatorial Bondi-like neutral inflow and bipolar ionized outflow. Since the radiation is confined to the polar directions, the rapid accretion proceeds in spite of the BH luminosity eight times larger than the Eddington value.

We have investigated such steady flow structure with the analytical models. The opening angle of the equatorial neutral layer $\theta_{\mathrm{inflow}}$ is derived from the balance between the supply and consumption rates of ionizing photons in each direction (equation \[eq:6\]). In turn, the accretion rate $\dot{M}$ is modeled assuming a Bondi-like flow through this equatorial layer also considering the photoevaporation mass loss from its surfaces (equation \[eq:18\]). We have also found that, in order for the equatorial Bondi-like inflow to be maintained, $\theta_{\mathrm{inflow}}$ at the Bondi radius must be above a critical value, which is $\simeq
4^\circ$ for $M_{\mathrm{BH}}=10^3\,M_\odot$ and $n_\infty=10^5{\,\mathrm{cm^{-3}}}$. This value is raised up to $\sim 10^\circ$ if we account for the mass loss inside the sink. The parameter dependence of the flow structure found in our simulations is well reproduced by this analytical model.

Our results highlight the importance of the directional dependence of BH irradiation, especially in the equatorial directions, in determining $\dot{M}$. However, our current knowledge about the actual anisotropy is very limited. Although not exactly the system of our interest, line-driven disc winds around a supermassive BH (SMBH) with $M_{\mathrm{BH}}
\sim 10^8\,M_\odot$ is shown to create anisotropic radiation fields by blocking high-energy photons [e.g., @Proga:2000aa; @Proga:2004aa; @Nomura:2016aa]. @Proga:2000aa suggest that the opening angle of the resulting shadow is $\theta_{\mathrm{shadow}} \simeq 12^\circ$, which is larger than the critical angle. In reality, the accretion disc may undergo the precession with variable angular momentum of accreting gas and change the orientation of the shadowed region in time, resulting in the destruction of the pre-existing neutral inflowing region. In any case, it is clearly awaited to study the structure of the inner part and resulting anisotropy of the BH irradiation. Our current study is complementary to such future works, because our results provide outer boundary conditions for them.

We have found that the required condition for the rapid accretion is substantially relaxed from that obtained for the isotropic irradiation. For these cases, the accretion rate is reduced to $\lesssim
0.01\,\dot{M}_{\mathrm{B}}$ by the radiative feedback unless the condition $(M_{\mathrm{BH}}/10^4M_\odot)(n_\infty/10^5{\,\mathrm{cm^{-3}}})\gtrsim 1$ [e.g., @Inayoshi:2016ac] is satisfied. This condition requires the ambient density $n_\infty$ as high as $10^6{\,\mathrm{cm^{-3}}}$ even for the most massive Pop III remnants with $M_{\mathrm{BH}}\sim 10^3M_\odot$ in [@Hirano:2015aa]. Such high ambient density seems difficult to achieve because the typical central density of the first galaxies at $z\sim 15$ is estimated as $10^5{\,\mathrm{cm^{-3}}}$ [e.g., @Oh:2002aa; @Volonteri:2005aa], although it is theoretically possible if the BH resides at the very centre of a halo with the density profile $\rho \propto r^{-2}$ [e.g., @Wise:2007aa; @Inayoshi:2016ac]. With the shadowing effect, BHs in a central part of the first galaxies can grow much more quickly.

A long-term evolution of such fast mass growth would be as follows. Suppose that a seed BH with mass $M_{\mathrm{BH}}=10^3M_\odot$ is embedded in an ambient medium with $n_\infty= 10^5{\,\mathrm{cm^{-3}}}$. Using the shadow opening angle of $\theta_{\mathrm{shadow}}=12^\circ$ suggested by @Proga:2000aa for line-driven SMBH winds, we obtain the accretion rate of $\dot{M} \sim 0.1\,\dot{M}_{\mathrm{B}} \sim
2\times10^{-4}(M_{\mathrm{BH}}/10^3M_\odot)^2\, M_\odot{\,\mathrm{yr^{-1}}}$ with equations  and . Integrating this expression, we obtain the growth history of the BH mass as $$M_{\mathrm{BH}}(t)\sim  \frac{10^3\,M_\odot}{1-\left[(t-t_0)/5{\,\mathrm{Myr}}\right]}\,,$$ where $t_0$ is the initial time of the accretion. At a face value, the BH mass diverges within a short timescale of $5{\,\mathrm{Myr}}$. In reality, however, the BH mass growth via accretion should be limited by changes in the environmental conditions, such as exhaustion of the ambient gas by accretion. If the remnant BHs of Pop III stars with $\sim
10^{2-3}\,M_\odot$ grow immediately in a few Myr time-scale by accretion to $\sim 10^{5-6} M_\odot$, they subsequently evolve in the same way as direct collapse BHs and can eventually grow to $\sim 10^9M_\odot$ SMBHs via gas accretion and/or mergers by $z \sim 7$ [see, e.g., @Tanaka:2009aa]. With the shadowing effect, Pop III remnants can be seeds for high-$z$ SMBHs. It does not mean, however, that all the Pop III remnants experience such rapid growth. For example, if they stay in a low-density region with $n<10 {\,\mathrm{cm^{-3}}}$, as suggested by [@Alvarez:2009aa], they hardly grow in mass even at the Bondi accretion rate.

Although we have assumed the weak rotation of the ambient gas, with which a BH accretion disk should be much smaller than the Bondi radius, accreting gas may have higher amount of the angular momentum in general. Regarding the SMBH accretion, [@Li:2013aa] has shown that the accretion rate is considerably reduced by the rotation [see also @Proga:2003ab; @Proga:2003aa]. If the flow predominantly comes from the equatorial plane, the amount of the angular momentum carried, and hence the impact of the rotational support, would be increased. To investigate this effect, mechanisms for the angular momentum transfer should also be considered (see below).

Our simulations also neglect the gas self-gravity. The torque caused by the self-gravity can play an important role in the angular momentum transport [see, e.g., @Shlosman:1989aa; @Shlosman:2016aa]. In addition, the inward force of the self-gravity can enhance the gas accretion on to BHs [e.g., @Li:2011aa], on the scales larger than both of the following two: (1) the radius where the enclosed gas mass equals to the BH mass, $r_{\mathrm{eq}}=(3\,M_{\mathrm{BH}}/4\pi\rho)^{1/3}\sim
10^5\, (M_{\mathrm{BH}}/10^3\,
M_\odot)^{1/3}(n/10^5\,{\,\mathrm{cm^{-3}}})^{-1/3}\,{\mathrm{AU}}$, and (2) the Jeans length $\lambda_{\mathrm{J}}=\sqrt{\pi}c_{\mathrm{s}}/\sqrt{\rho G}\sim
10^6\,(n/10^5\,{\,\mathrm{cm^{-3}}})^{-1/2}\,{\mathrm{AU}}$. In our cases, however, the Bondi radii are smaller than those scales and thus the flow structures on the scale of Bondi radius are hardly affected by this effect, although the size of [H[ii]{} ]{}bubbles can exceed them in some cases.

To be realistic, the following improvements are needed. First of all, the gas dynamics should be followed in 3D, in particular, to see the effect of the gravitational torque. Next, diffuse recombination photons can modify the structure of neutral layer. Furthermore, solving the gas dynamics beyond our computational domain, i.e, in the outer molecular and photodissociation regions [e.g., @Ricotti:2001aa], will be needed to see the large-scale flow structure. Finally, it is crucial to perform numerical simulations for the inner part dedicated to resolving the generation of anisotropic radiation fields, which are inevitably coupled to our simulations through the outer boundary conditions for them. We have clearly shown that the interplay of multi-scale processes is essential in understanding the BH accretion.

Acknowledgements {#acknowledgements .unnumbered}
================

The authors would like to thank Kazumi Kashiyama, Rohta Takahashi, Sanemichi Takahashi and Kenji Toma for fruitful discussions. The numerical simulations were performed on the Cray XC30 at CfCA of the National Astronomical Observatory of Japan, as well as on the computer cluster, [Draco]{}, at Frontier Research Institute for Interdisciplinary Sciences of Tohoku University. This work is supported in part by MEXT/JSPS KAKENHI Grant Number 15J03873 (KS), 25800102, 15H00776 and 16H05996 (TH), 15H06022 (HY) and 25287040 (KO).

natexlab\#1[\#1]{}

, T., [Anninos]{}, P., [Zhang]{}, Y., & [Norman]{}, M. L. 1997, New Astron., 2, 181

, M. A., [Czerny]{}, B., [Lasota]{}, J. P., & [Szuszkiewicz]{}, E. 1988, , 332, 646

, B., [Khochfar]{}, S., [Johnson]{}, J. L., [et al.]{} 2012, , 425, 2854

, S. M. V., & [Pequignot]{}, D. 1973, , 25, 137

, T., & [Natarajan]{}, P. 2014, Science, 345, 1330

, M. A., [Wise]{}, J. H., & [Abel]{}, T. 2009, , 701, L133

, P., [Zhang]{}, Y., [Abel]{}, T., & [Norman]{}, M. L. 1997, New Astron., 2, 209

, J. G., [Centrella]{}, J., [Choi]{}, D.-I., [et al.]{} 2006, , 653, L93

, P., [Proga]{}, D., & [Nagamine]{}, K. 2012, , 424, 728

, M. C. 1978, , 184, 53

, M. C., & [Volonteri]{}, M. 2016, arXiv:1609.07137v1

, J. H. 1981, , 197, 553

, R. D., & [Begelman]{}, M. C. 1999, , 303, L1

, I., [Burgess]{}, A., [Fursa]{}, D. V., & [Tully]{}, J. A. 2000, , 146, 481

, R. 1992, , 78, 341

, S., [Hirano]{}, S., [Hosokawa]{}, T., & [Yoshida]{}, N. 2016, , 832, 134

, R. E. S., & [Harrington]{}, J. P. 1989, , 239, 869

, B., & [Volonteri]{}, M. 2009, , 694, 302

, M., [Ferrara]{}, A., & [Mesinger]{}, A. 2014, , 442, 2036

, M., [Haiman]{}, Z., [Mesinger]{}, A., & [Wyithe]{}, J. S. B. 2008, , 391, 1961

, B. T. 2011, [Physics of the Interstellar and Intergalactic Medium]{} (Princeton University Press)

, X., [Narayanan]{}, V. K., [Lupton]{}, R. H., [et al.]{} 2001, , 122, 2833

, G. J., [Peterson]{}, B. M., [Horne]{}, K., [Welsh]{}, W. F., & [Nahar]{}, S. N. 1992, , 387, 95

, P. C., [Olejar]{}, A., & [Anninos]{}, P. 2014, , 796, 22

, G., & [Franco]{}, J. 1996, , 469, 171

, S. C. O., & [Jappsen]{}, A.-K. 2007, , 666, 1

, Z. 2013, in Astrophysics and Space Science Library, Vol. 396, The First Galaxies, ed. T. [Wiklind]{}, B. [Mobasher]{}, & V. [Bromm]{}, 293

, S., [Hosokawa]{}, T., [Yoshida]{}, N., [Omukai]{}, K., & [Yorke]{}, H. W. 2015, , 448, 568

, D., [Johnstone]{}, D., [Lizano]{}, S., & [Shu]{}, F. 1994, , 428, 654

, T., [Hirano]{}, S., [Kuiper]{}, R., [et al.]{} 2016, , 824, 119

, T., [Omukai]{}, K., [Yoshida]{}, N., & [Yorke]{}, H. W. 2011, Science, 334, 1250

, D. G., & [Storey]{}, P. J. 1998, , 297, 1073

, K., [Haiman]{}, Z., & [Ostriker]{}, J. P. 2016, , 459, 3738

, K., & [Tanaka]{}, T. L. 2015, , 450, 4350

, R. K., [Langer]{}, W. D., & [Evans]{}, K. 1987, [Elementary processes in Hydrogen-Helium plasmas - Cross sections and reaction rate coefficients]{} (Springer)

, M., [Pawlik]{}, A. H., [Greif]{}, T. H., [et al.]{} 2012, , 754, 34

, Y.-F., [Stone]{}, J. M., & [Davis]{}, S. W. 2014, , 796, 106

, S., [Fukue]{}, J., & [Mineshige]{}, S. 1998, [Black-hole accretion disks]{} (Kyoto University Press)

, H., [Sijacki]{}, D., & [Haehnelt]{}, M. G. 2015, , 451, 2352

, M., [Lane]{}, N. F., [Dalgarno]{}, A., & [Dixson]{}, R. G. 1993, , 405, 801

, M., [Pollney]{}, D., [Reisswig]{}, C., [et al.]{} 2007, Physical Review Letters, 99, 041102

, R., [Klahr]{}, H., [Beuther]{}, H., & [Henning]{}, T. 2010, , 722, 1556

—. 2011, , 732, 20

, R., [Klahr]{}, H., [Dullemond]{}, C., [Kley]{}, W., & [Henning]{}, T. 2010, , 511, A81

, R., & [Klessen]{}, R. S. 2013, , 555, A7

, R., & [Proga]{}, D. 2009, , 693, 1929

, J., [Ostriker]{}, J., & [Sunyaev]{}, R. 2013, , 767, 105

, Y. 2011, ArXiv e-prints, arXiv:1109.3442

, P., [Haardt]{}, F., & [Dotti]{}, M. 2014, , 784, L38

, C. F., & [Tan]{}, J. C. 2008, , 681, 771

, J. C., [Tchekhovskoy]{}, A., [Sadowski]{}, A., & [Narayan]{}, R. 2014, , 441, 3177

, A., [Bodo]{}, G., [Massaglia]{}, S., [et al.]{} 2007, , 170, 228

, M., [Bromm]{}, V., [Couch]{}, S. M., & [Oh]{}, S. P. 2009, , 698, 766

, M., [Couch]{}, S. M., & [Bromm]{}, V. 2009, , 696, L146

, D. J., [Warren]{}, S. J., [Venemans]{}, B. P., [et al.]{} 2011, , 474, 616

, M., [Ohsuga]{}, K., [Takahashi]{}, H. R., [Wada]{}, K., & [Yoshida]{}, T. 2016, , 68, 16

, M., [Ohsuga]{}, K., [Wada]{}, K., [Susa]{}, H., & [Misawa]{}, T. 2013, , 65, arXiv:1212.3075

, G. S., [Ostriker]{}, J. P., & [Ciotti]{}, L. 2011, , 737, 26

, S. P., & [Haiman]{}, Z. 2002, , 569, 558

, K., [Mori]{}, M., [Nakamoto]{}, T., & [Mineshige]{}, S. 2005, , 628, 368

, K., [Schneider]{}, R., & [Haiman]{}, Z. 2008, , 686, 801

, D. E. 1989, [Astrophysics of gaseous nebulae and active galactic nuclei]{} (University Science Books)

, F., & [Ferrara]{}, A. 2015, , 448, 104 , F., [Salpeter]{}, E. E., & [Stahler]{}, S. W. 1983, , 271, 632

, K., & [Ricotti]{}, M. 2011, , 739, 2

—. 2012, , 747, 9

—. 2013, , 767, 163

, K., [Ricotti]{}, M., [Di Matteo]{}, T., & [Reynolds]{}, C. S. 2014, , 437, 2856

, K., [Ricotti]{}, M., [Natarajan]{}, P., [Bogdanovi[ć]{}]{}, T., & [Wise]{}, J. H. 2016, , 818, 184

, G., [Fender]{}, R. P., [Begelman]{}, M. C., [et al.]{} 2012, , 422, 11

, D. 2007, , 661, 693

, D., & [Begelman]{}, M. C. 2003, , 582, 69

—. 2003, , 592, 767

, D., & [Kallman]{}, T. R. 2004, , 616, 688

, D., [Stone]{}, J. M., & [Kallman]{}, T. R. 2000, , 543, 686

, M., [Gnedin]{}, N. Y., & [Shull]{}, J. M. 2001, , 560, 580

—. 2002, , 575, 33

, Y., [Inayoshi]{}, K., & [Haiman]{}, Z. 2016, , 461, 4496

, A., & [Narayan]{}, R. 2016, , 456, 3929

, N. I., & [Sunyaev]{}, R. A. 1973, , 24, 337

, P. R., & [Kang]{}, H. 1987, , 318, 32

, I., [Choi]{}, J.-H., [Begelman]{}, M. C., & [Nagamine]{}, K. 2016, , 456, 500

, I., [Frank]{}, J., & [Begelman]{}, M. C. 1989, , 338, 45

, J. M. 1979, , 234, 761

, J. M., & [van Steenberg]{}, M. E. 1985, , 298, 268

, K., [Coppola]{}, C. M., [Omukai]{}, K., [Galli]{}, D., & [Palla]{}, F. 2016, , 456, 270

, K., [Omukai]{}, K., & [Inoue]{}, A. K. 2014, , 445, 544

, H., [Hasegawa]{}, K., & [Tominaga]{}, N. 2014, , 792, 32

, H. R., & [Ohsuga]{}, K. 2015, , 67, 60

, K. E. I., [Nakamoto]{}, T., & [Omukai]{}, K. 2013, , 773, 155

, T., & [Haiman]{}, Z. 2009, , 696, 1798

, A., & [Umemura]{}, M. 2011, , 728, L31

, B. P., [Findlay]{}, J. R., [Sutherland]{}, W. J., [et al.]{} 2013, , 779, 24

, M. 2012, Science, 337, 544

, M., & [Rees]{}, M. J. 2005, , 633, 624

, M., [Silk]{}, J., & [Dubus]{}, G. 2015, , 804, 148

, K.-y., [Fukue]{}, J., [Takeuchi]{}, M., & [Mineshige]{}, S. 2000, , 52, 133

, D., & [Norman]{}, M. L. 2008, , 673, 664

, D. J., & [Norman]{}, M. L. 2008, , 672, 287

, C. J., [Delorme]{}, P., [Reyl[é]{}]{}, C., [et al.]{} 2010, , 139, 906

, J. H., & [Abel]{}, T. 2007, , 665, 899

, X.-B., [Wang]{}, F., [Fan]{}, X., [et al.]{} 2015, , 518, 512

, H., & [Khochfar]{}, S. 2016, , 457, 2423

, M., [Sadeghpour]{}, H. R., & [Dalgarno]{}, A. 1998, , 496, 1044

, N., [Omukai]{}, K., & [Hernquist]{}, L. 2008, Science, 321, 669

, F., [Gan]{}, Z., [Narayan]{}, R., [et al.]{} 2015, , 804, 101

, F., [Abbassi]{}, S., & [Mosallanezhad]{}, A. 2016, , 823, 92

, B., [Dalgarno]{}, A., [Kimura]{}, M., & [Lane]{}, N. F. 1989, , 40, 2340

details of chemical and thermal modelling {#sec:chem_detail}
=========================================

Reaction rates {#sec:reaction_rates}
--------------

  No.       Reaction                                                                        Rate coeff. ${\mathrm{[cm^3\,s^{-1}]}}$                                                                                          Ref.
  --------- ------------------------------------------------------------------------------- ------------------------------------------------------------------------------------------------------------------------------- ------
  1         ${\mathrm{H}}+{\mathrm{e}} \rightarrow {\mathrm{H^+}} + 2{\mathrm{e}}$          $k_1=$                                                                                                                            1
  2         ${\mathrm{He}}+{\mathrm{e}} \rightarrow {\mathrm{He^+}} + 2{\mathrm{e}}$        $k_2=$                                                                                                                            1
  3         ${\mathrm{He^+}}+{\mathrm{e}} \rightarrow {\mathrm{He^{2+}}} + 2{\mathrm{e}}$   $k_3=$                                                                                                                            2
  $4^{a}$   ${\mathrm{H^+}}+{\mathrm{e}} \rightarrow {\mathrm{H}} + h\nu$                   $k_4=2.753\times10^{-14}  (T_{\mathrm{K}}/315614)^{-3/2}  (1+(T_{\mathrm{K}}/115188)^{-0.407})^{-2.242}$                          3
  $5^{b}$   ${\mathrm{He^+}}+{\mathrm{e}} \rightarrow {\mathrm{He}} + h\nu$                 $k_5=k_{\mathrm{5rr}}+k_{\mathrm{5di}}$                                                                                         
                                                                                            $k_{\mathrm{5rr}}=$                                                                                                               4
                                                                                            $k_{\mathrm{5di}}=1.9\times10^{-3}  T_{\mathrm{K}}^{-3/2} \exp[-473421/T_{\mathrm{K}}]  (1 + 0.3\exp[-94684/T_{\mathrm{K}}])$     5
  $6^{c}$   ${\mathrm{He^{2+}}}+{\mathrm{e}} \rightarrow {\mathrm{He^+}} + h\nu$            $k_6=5.08\times10^{-13} (T_{\mathrm{K}}/40000)^{-0.8163-0.0208\log_{10}(T_{\mathrm{K}}/40000)}$                                   6
  7         ${\mathrm{He^+}}+{\mathrm{H}} \rightarrow {\mathrm{He}} + {\mathrm{H^+}}$       $k_7=1.25\times10^{-15}(T_{\mathrm{K}}/300)^{0.25}$                                                                               7
  8         ${\mathrm{H^+}}+{\mathrm{He}} \rightarrow {\mathrm{H}} + {\mathrm{He^+}}$       $k_8=$                                                                                                                            8
  9         $2{\mathrm{H}} \rightarrow {\mathrm{H^+}} + {\mathrm{H}} + {\mathrm{e}}$        $k_9=1.7\times10^{-4}k_1$                                                                                                         9

\

NOTES. The $T_{\mathrm{K}}$ and $T_{\mathrm{eV}}$ are the gas temperature in units of K and eV, respectively; $^{a}$Case B; $^{b}$radiative [Case B; singlet; our fit to @Hummer:1998aa] and dielectric [@Aldrovandi:1973aa] recombination; $^{c}$Case B [@Draine:2011aa with typo about the charge dependence corrected].\
REFERENCES. (1) [@Janev:1987aa]; (2) [@Abel:1997aa from Aladdin database 1989]; (3) [@Ferland:1992aa]; (4) [@Hummer:1998aa]; (5) [@Aldrovandi:1973aa]; (6) [@Draine:2011aa]; (7) [@Zygelman:1989aa]; (8) [@Kimura:1993aa]; (9) [@Palla:1983aa].

In Table \[tab:reaction\_rates\], we summarize the chemical reactions considered in this work, which are adopted following [@Glover:2007aa], [@Abel:1997aa] and [@Anninos:1997aa]. We adopt the Case B recombination rates for the recombination of ${\mathrm{H^+}}$, ${\mathrm{He^+}}$ and ${\mathrm{He^{2+}}}$. We neglect the ${\mathrm{He^+}}$ recombination through the quasi-stable triplet state $2^3S$ of ${\mathrm{He}}$, assuming that ${\mathrm{He}}$ in that state is easily photoionized by the BH irradiation [see, e.g., @Clegg:1989aa].

Cross sections {#sec:cross_sections}
--------------

  No.   Reaction                                                               Cross section ${\mathrm{[cm^2]}}$                                                  Ref.
  ----- ---------------------------------------------------------------------- ----------------------------------- --------------------------------------------- ------
  1     ${\mathrm{H}}+h\nu \rightarrow {\mathrm{H^+}} + {\mathrm{e}}$          $\sigma_{\nu,1}=$                   $h\nu_{\mathrm{T,1}}= 13.60{\,\mathrm{eV}}$     1
  2     ${\mathrm{He}}+h\nu \rightarrow {\mathrm{He^+}} + {\mathrm{e}}$        $\sigma_{\nu,2}=$                   $h\nu_{\mathrm{T,2}}= 24.58{\,\mathrm{eV}}$     2
  3     ${\mathrm{He^+}}+h\nu \rightarrow {\mathrm{He^{2+}}} + {\mathrm{e}}$   $\sigma_{\nu,3}=$                   $h\nu_{\mathrm{T,3}}= 54.40{\,\mathrm{eV}}$     1

\

REFERENCES. (1) [@Osterbrock:1989aa]; (2) [@Yan:1998aa].

In Table \[tab:cross\_sections\], we summarize the cross sections considered in this work.

Heating and cooling rates {#sec:heat_cool}
-------------------------

[lllc]{} No. & Process & Rate \[${\mathrm{erg\, cm^{-3}\, s^{-1}}}$\] & Ref.\
\
1 & ${\mathrm{H}}$ photoionization & $\Gamma_1$ (see text)&\
2 & ${\mathrm{He}}$ photoionization & $\Gamma_2$ (see text)&\
3 & ${\mathrm{He^+}}$ photoionization & $\Gamma_3$ (see text)&\
\
$1^{a}$ & ${\mathrm{H^+}}$ recombination & $\Lambda_1=$ & 1\
$2^{b}$ & ${\mathrm{He^+}}$ recombination & & 2,3\
$3^{c}$ & ${\mathrm{He^{2+}}}$ recombination & $\Lambda_3= 1.38\times10^{-16}T_{\mathrm{K}} \left(0.684-0.0416\,\ln(T_{\mathrm{K}}/40000)\right)  
k_6\,n({\mathrm{e}})\,n({\mathrm{He^{2+}}})$ & 4\
4 & ${\mathrm{H}}$ excitation & $\Lambda_4=7.50\times10^{-19}  \left(1 + (T_{\mathrm{K}}/100000)^{1/2}\right)^{-1}  \exp[-118348/T_{\mathrm{K}}]
\,n({\mathrm{e}})\,n({\mathrm{H}})$ & 5\
$5^{d}$ & ${\mathrm{He}}$ excitation & $\Lambda_5=1.1\times10^{-19}  T_{\mathrm{K}}^{0.082}  \exp[-230000/T_{\mathrm{K}}]
\,n({\mathrm{e}})\,n({\mathrm{He}})$ & 6\
6 & ${\mathrm{He^+}}$ excitation & $\Lambda_6=5.54\times10^{-17} T_{\mathrm{K}}^{-0.397}  \left(1 + (T_{\mathrm{K}}/100000)^{1/2}\right)^{-1}  \exp[-473638/T_{\mathrm{K}}]
\,n({\mathrm{e}})\,n({\mathrm{He^{+}}})$ & 5\
7 & ${\mathrm{H}}$ ionization & $\Lambda_7=2.18\times10^{-11}k_1\,n({\mathrm{e}})\,n({\mathrm{H}})$ & 7\
$8^{c}$ & ${\mathrm{He}}$ ionization & $\Lambda_8=3.94\times10^{-11}k_2\,n({\mathrm{e}})\,n({\mathrm{He}})$ & 7\
9 & ${\mathrm{He^+}}$ ionization & $\Lambda_9=8.72\times10^{-11}k_3\,n({\mathrm{e}})\,n({\mathrm{He^{+}}})$ & 7\
10 & Free-free & & 8\
$11^{e}$ & Compton & $\Lambda_{11}=1.017\times10^{-37} T_{\mathrm{CMB}}^4  (T_{\mathrm{K}} - T_{\mathrm{CMB}}) \,n({\mathrm{e}})$ & 5\

\

NOTES. $^{a}$Case B, our fit to [@Ferland:1992aa]; $^{b}$radiative [Case B; singlet; our fit to @Hummer:1998aa] and dielectric [@Black:1981aa] recombination cooling; $^{c}$Case B [@Draine:2011aa with typo about the charge dependence corrected]; $^{d}$singlet; $^{e}$$T_{\mathrm{CMB}}=2.73(1 + z)$ with $z=15$.\
REFERENCES. (1) [@Ferland:1992aa]; (2) [@Hummer:1998aa]; (3) [@Black:1981aa] (4) [@Draine:2011aa]; (5) [@Cen:1992aa]; (6) [@Bray:2000aa]; (7) [@Anninos:1997aa]; (8) [@Shapiro:1987aa].

In Table \[tab:cool\_rates\], we summarize the heating and cooling processes considered in this work. Here, $n(X)$ is the number density of species $X$ in units of ${\mathrm{cm^{-3}}}$. We calculate the photoionization heating rates as $\Gamma_i=\int (4\pi j_\nu/h\nu)
n(X_i)\sigma_{\nu,i} (h\nu - h\nu_{{\mathrm{T}},i}){\mathrm{d}}\nu$ (see Table \[tab:cross\_sections\]), with $X_1={\mathrm{H}}$, $X_2= {\mathrm{He}}$ and $X_3={\mathrm{He^{+}}}$.

Resolution check {#sec:res_check}
================

![Same as Fig. \[fig:mdot\] but for the runs in App. \[sec:res\_check\]. The physical parameters are the same as Dds run but the resolution is different in each run. See the text for details.[]{data-label="fig:mdot_rdep"}](figure/mdot_rdep.eps){width="8.5cm"}

To check the resolution dependence of our results, we here see how the evolution of $\dot{M}$ is affected by numerical settings, namely the number of grids and sink size $R_{\mathrm{in}}$. Taking the same physical parameters as Dds run, we perform additional simulations with different resolutions, as shown in Fig. \[fig:mdot\_rdep\]. Here, we take $N_r
\times N_\theta = 512\times 144$; $N_r \times N_\theta = 256\times 72$; $N_r \times N_\theta = 128\times 36$; $N_r \times N_\theta = 256\times
72$ with $R_{\mathrm{in}}$ halved and doubled from the fiducial value. Note that our main results are obtained with the high- and medium-resolution simulations with $N_r \times N_\theta = 512\times 144$ and $256\times
72$, respectively.

The dependence on the number of grids is checked by comparing the results with $N_r \times N_\theta = 512\times 144$, $256\times 72$, and $128\times 36$ (Fig. \[fig:mdot\_rdep\]). The strong variability of $\dot{M}$ for $t \lesssim 10^6{\,\mathrm{yr}}$ seen with the highest-resolution is smoothed out with the lower resolutions. However, the values of $\dot{M}$ at the end of the simulations are almost the same in all three cases. This confirms that the conclusion of this paper does not depend on the number of grids.

Fig. \[fig:mdot\_rdep\] also shows the evolution of $\dot{M}$ for the cases with the different sink sizes. The differences of accretion rates are less than 10 % for $t > 1.5 \times 10^6$ years in all three cases. We also find that the values of $\dot{M}$ at the end of the simulations decrease only by 4% by halving the sink size from the fiducial value. Such a trend is consistent with the estimated mass-loss rate from the region between the halved and fiducial inner boundaries (see equation \[eq:14\]), although it is also within the numerical error. This ensures that the dependence of our results on the sink size is weak.

[^1]: E-mail: sugimura@astr.tohoku.ac.jp

[^2]: For a case with $M_{\mathrm{BH}}=10^3\,M_\odot$ and $L=L_{\mathrm{E}}$ with the spectrum $L_\nu \propto \nu^{-1.5}$, the specific FUV intensity at $r_{\mathrm{B}}$ is $J_{21}\sim 10^9$ (in units of $10^{-21}{\,\mathrm{erg\,s^{-1}\,Hz^{-1}\,sr^{-1}\,cm^{-2}}}$), while the critical intensity for totally suppressing ${\mathrm{H_2}}$ formation in atomic cooling halos is $J_{21,cr}\sim 10^3$ [see, e.g., @Sugimura:2014aa].

[^3]: Observations of Galactic stellar BHs support that disc winds and associated shadowed regions only exist in the high/soft state and disappear in the low/hard state [e.g., see @Ponti:2012aa].

[^4]: With equation , we are now able to explain why the accretion rate temporarily exceeds the Bondi rate in the early stage of the simulation in Fig. \[fig:mdot\](c). In the beginning, the hot bipolar [H[ii]{} ]{}bubbles compresses the equatorial neutral layer. Consequently, the corresponding Bondi-like accretion rate increases and becomes larger than the original Bondi rate $\dot{M}_{\mathrm{B}}$ even after multiplying the fraction of the neutral solid angle $\Delta\Omega_{\mathrm{inflow}}(r_{\mathrm{B}})/(4\pi)$.
