The effect of angular momentum and radiation on accretion

Stunted accretion growth of black holes by combined effect of the flow angular momentum and radiation feedback

Abstract

Accretion on to seed black holes (BHs) is believed to play a crucial role in formation of supermassive BHs observed at high-redshift (). Here, we investigate the combined effect of gas angular momentum and radiation feedback on the accretion flow, by performing 2D axially symmetric radiation hydrodynamics simulations that solve the flow structure across the Bondi radius and the outer part of the accretion disc simultaneously. The accreting gas with finite angular momentum forms a rotationally-supported disc inside the Bondi radius, where the accretion proceeds by the angular momentum transport due to assumed -type viscosity. We find that the interplay of radiation and angular momentum significantly suppresses accretion even if the radiative feedback is weakened in an equatorial shadowing region. The accretion rate is times the Bondi value, where is the viscosity parameter. By developing an analytical model, we show that such a great reduction of the accretion rate persists unless the angular momentum is so small that the corresponding centrifugal radius is times the Bondi radius. We argue that BHs are hard to grow quickly via rapid mass accretion considering the angular momentum barrier presented in this paper.

keywords:
quasars: supermassive black holes-cosmology: theory.

1 Introduction

Observations of supermassive black holes (SMBHs) with mass at redshift , or after the big bang, severely constrain their formation mechanism (e.g., Fan et al., 2001; Willott et al., 2010; Mortlock et al., 2011; Venemans et al., 2013; Wu et al., 2015; see also Gallerani et al., 2017 for review). Heavy seed BHs have been invoked in several scenarios to reconcile the short available time with their hugeness (see, e.g., Volonteri, 2012; Haiman, 2013, for a review), which includes (1) Pop III remnant BHs with mass of (e.g., Yoshida et al., 2008; Hosokawa et al., 2011, 2016; Susa et al., 2014; Hirano et al., 2015; Stacy et al., 2016); (2) direct collapse BHs with formed via the collapse of supermassive stars (e.g., Omukai, 2001; Bromm & Loeb, 2003; Hosokawa et al., 2012; Sugimura et al., 2014, 2016; Inayoshi et al., 2014; Chon et al., 2016; Umeda et al., 2016); and (3) massive BHs with formed as a consequence of stellar mergers in dense clusters (e.g., Omukai et al., 2008; Devecchi & Volonteri, 2009; Katz et al., 2015; Tagawa et al., 2015; Yajima & Khochfar, 2016; Sakurai et al., 2017).

In all the cases, the seeds have to increase their mass further by several orders of magnitude. Although Pop III remnants with relatively small mass have been claimed to become SMBHs by by very rapid (super-Eddington) gas accretion (Madau et al., 2014; Alexander & Natarajan, 2014; Volonteri et al., 2015), whether this actually occurs or not is still very uncertain. This is because our knowledge on the realistic accretion process is still limited.

Bondi (1952) analytically investigated the most simplistic case of steady, spherical accretion of the polytropic gas from a homogeneous and static medium (see also, e.g., Shapiro & Teukolsky, 1983). For a BH with mass surrounded by a gas with density and sound speed , the so-called Bondi accretion rate is

(1)

Here, is the non-dimensional factor depending on the polytropic index ( for the isothermal, i.e., , case) and for a neutral primordial gas with temperature . In the second expression of Eq. (1), we use the relation between the number density of hydrogen atoms and for the primordial gas.

Although the Bondi rate provides a rough estimate of an accretion rate onto a BH, various effects reduce the rate in realistic situations. For instance, radiation feedback from circum-BH discs can significantly disturb the accretion flow so that the rate is decreased down to (e.g., Milosavljević et al., 2009; Park & Ricotti, 2011, 2012; Park et al., 2017). Recent studies have proposed pathways to alleviate the feedback, in such cases as where the radiation is highly obscured around the equatorial plane by disc winds (Sugimura et al., 2017; Takeo et al., 2018) or where the radiation is trapped in a dense (super-Eddington) accretion flow (Inayoshi et al., 2016; Sakurai et al., 2016). However, these studies have considered the accretion flow without angular momentum, assuming that the size of the circum-BH disc, if any, is much smaller than the Bondi radius. It is most likely that such an assumption is not always reasonable, but we do not know the condition on angular momentum that such a simplified treatment is justified.

Accretion to active galactic nuclei (AGNs) has been investigated in the last few decades (e.g., Ciotti & Ostriker, 2001; Wada & Norman, 2002; Kawakatu & Wada, 2008; Kurosawa & Proga, 2009; Novak et al., 2011; Kawaguchi & Mori, 2010; Barai et al., 2012; Yuan et al., 2012). In such context, effect of gas angular momentum has partly been investigated. Even without radiation feedback, moderate gas angular momentum significantly suppresses the accretion in low-luminosity AGN (Proga & Begelman, 2003b, a; Cuadra et al., 2006; Li et al., 2013; Gaspari et al., 2015; Inayoshi et al., 2018). This implies the effect of angular momentum, as well as radiation feedback, needs to be considered in studying the accretion on to the seed BHs.

In this work, to see the combined effect of the angular momentum and the radiation feedback on seed BH accretion, we perform a set of 2D axisymmetric radiation hydrodynamics simulations, considering both finite gas angular momentum and radiation from the circum-BH discs. We follow formation of rotationally-supported discs and subsequent viscous accretion through them. We find that angular momentum of the gas, in cooperation with radiation feedback, suppresses the accretion on to seed BHs. To understand its mechanism, we also develop an analytical model for accretion from a rotating medium.

The paper is organized as follows. In Sec. 2, we describe the numerical method and the parameter sets explored. In Sec. 3, we present the results of our simulations. In Sec. 4, we develop an analytical model, which allow us to obtain the condition for suppression of accretion. In Sec. 5, we discuss the possible growth history of Pop III remnant BHs based on our findings, as well as the caveats of our simulations. Finally, we summarize in Sec. 6.

2 Method

We perform axisymmetric 2D radiation hydrodynamics (RHD) simulations, by using a modified version of a public multidimensional magnetohydrodynamics (MHD) code PLUTO 4.1 (Mignone et al., 2007), which is mostly the same as used in our previous work (Sugimura et al., 2017; see also Kuiper et al., 2010a, b, 2011; Kuiper & Klessen, 2013; Hosokawa et al., 2016). Major update here is implementation of physics related to the rotation of gas. In the rest of this paper, we use both spherical and cylindrical coordinates interchangeably, although the spherical one was actually used in the calculation.

2.1 Modelling accretion on to a BH under radiation feedback

Figure 1: Directional dependence of ionizing flux. The radial extent represents the strength of the flux compared to the isotropic case.

We briefly describe the properties of the code that are common to the previous work, where we followed accretion on to BHs under radiation feedback (see Sugimura et al., 2017, for details).

We put a BH at the center of computational domain and treat it as a sink. Through the sink surface, the gas is allowed to flow in, while the radiation is emitted according to a semi-analytical model described below. We solve the coupled equations of hydrodynamics, radial multi-frequency radiation transport and primordial gas chemistry. In the current version of the code, helium ionization is considered, hydrogen molecules are assumed to be completely destroyed, and the self-gravity of gas is ignored.

We consider ionizing photons emitted from the circum-BH discs with a semi-analytical prescription. Using the accretion rate evaluated at the sink surface, the luminosity is given by (Watarai et al., 2000)

(2)

with . Here, is the Eddington luminosity,

(3)

and the (efficiency-independent) Eddington accretion rate,

(4)

Note that the decrease of radiative efficiency at is caused by photon trapping in the slim discs (Abramowicz et al., 1988). We assume that the spectral energy distribution is given by the power law with the frequency minimum at (e.g., Park & Ricotti, 2011).

Radiation from the sink is supposed to be anisotropic, partly because the photons from hot inner part of the circm-BH disc is obscured by some outer materials. For example, Proga et al. (2000) suggested that line-driven winds from a SMBH accretion disc make a shadowing region with opening angle from the equatorial plane . However, especially in the case with smaller BH mass or lower metallicity, the anisotropy is highly uncertain due to the lack of knowledge on obscuring materials, which are presumably (failed) winds or coronae above the disc1 (see, e.g., Hollenbach et al., 1994; Begelman et al., 1983; Woods et al., 1996; Proga et al., 2000; Wada, 2012; Suzuki & Inutsuka, 2014; Nomura et al., 2016). We thus do not attempt to realistically model the anisotropy, but instead, consider the two limiting cases: isotropic radiation and anisotropic radiation with a shadowing region (Fig. 1; Equation 9 of Sugimura et al., 2017). We expect the reality lies somewhere between the two.

2.2 Modelling gas rotation

We assume that the surrounding gas has the initial profile of specific angular momentum

(5)

with the Bondi radius

(6)

That is, the gas has constant throughout the computational domain, except near the rotation axis with , where the angular velocity is constant instead. Below, we interchangeably use and the centrifugal radius , where the centrifugal force and the BH gravity balances, to indicate the angular momentum of the accreted gas, as is related to as

(7)

Angular momentum must be transported for a gas with finite to reach the vicinity of the BH. We assume the -type viscous stress (Shakura & Sunyaev, 1973; see also Igumenshchev & Abramowicz, 1999; Stone et al., 1999; Li et al., 2013), which is possibly due to turbulence driven by the magnetorotational instability (MRI; Balbus & Hawley, 1998). The gravitational torque is insignificant in the cases studied here, because the disc is highly gravitationally stable (see Appendix E). Note that while the value of corresponding to the actual MRI turbulence is yet to be fully understood, Bai & Stone (2013), for example, have reported that is in the case with a weak vertical magnetic field but can be even larger than unity in the strong field case.

To assure that the viscosity works only inside the disc (e.g., Stone et al., 1999), we introduce a confinement factor and use the following expression for viscosity:2

(8)

with isothermal sound speed , adiabatic index and Keplerian angular velocity . Here, we identify the disc region based on the degree of rotational support against the BH gravity: we set for and for , and linearly interpolate between them, with the threshold value and the transition width in most of our calculations. We check effect of changing in Sec. A.2. Note that we have seen in test runs that if the viscosity is not limited to the disc, the disc will disappear and the flow becomes spherical due to angular momentum transport in the entire computational domain.

2.3 Cases considered

radiation
0.1 0.01 anisotropic
0.1 0.1 anisotropic
0.3 0.01 anisotropic
0.1 0.01 isotropic
0.1 0.1 isotropic
0.3 0.01 isotropic
0.1 0.01 no
0.1 0.1 no
0.3 0.01 no
Non-rotating case
0 anisotropic
0 isotropic
0 no

NOTES.—We set , and in all runs.

Anisotropic radiation with shadow.

Results from analytical estimate and the simulations in Sugimura et al. (2017).

Accretion rate at the end of simulation (averaged near the end of simulation if oscillating).

The accretion rate can also be normalized by using the relation (See Eqs. 1 and 4).

Table 1: Summary of runs.

Our runs are summarized in Table 1. In all runs, we set and fix it constant during the calculations. The surrounding medium is assumed to be the neutral primordial gas with density and temperature . In the fiducial run, we assume anisotropic radiation field with and . We also perform runs with isotropic radiation field, and in some cases change to or to to see the parameter dependence.

We start from the homogeneous and quasi-static () initial condition. We turn the radiation off for the first to allow the flow to settle in a steady state with a rotationally-supported disc. Note that this period is longer than either the dynamical time at , , or the viscous time at , . We then turn on the radiation and follow the evolution until in the anisotropic radiation runs and in the isotropic radiation runs.

We impose the reflection symmetry with respect to the equatorial plane, as well as axisymmetry around the rotation axis, and thus extends from to . In the -direction, our computational domain ranges from to . Note that is smaller than either (see above) or the Bondi radius for H ii gas with , while is larger than the size of H ii region. We use logarithmic grids in the direction and homogeneous ones in the directions, with . We discuss the dependence of our results on the numerical configuration in Appendix D.

At the outer boundary, we let the flow go out from the computational domain but not come into it. At the inner boundary, we impose the same boundary condition in most cases, but when the gas is judged to belong to the Keplerian accretion disc (specifically, when and , where quantities with subscripts “in” are evaluated at ), we determine the physical quantities in the ghost cells according to the radial profile of the isothermal Keplerian disc: , , , and for the other variables (see Appendix A.2). We set a temperature floor at for simplicity (e.g., Sugimura et al., 2017), as well as minimum density and maximum velocity for numerical reasons.

Figure 2: Time evolution of the (a) accretion rate and (b) luminosity in the runs with and . Both fiducial anisotropic radiation case (blue) and isotropic radiation case (orange) are shown. In these runs, the radiation is turned off for the initial years. The bottom thin-dashed lines represent the averaged values for . The results for the non-rotating case (Sugimura et al., 2017) are also shown by arrows for comparison.

3 Results

Results of our simulations are summarized in Table 1, along with our previous cases for the non-rotating gas (Sugimura et al., 2017), for comparison. Below, we first describe the fiducial anisotropic radiation run in Sec. 3.1. We then present the isotropic radiation run with the same parameter set in Sec. 3.2 and finally see the cases with different parameter sets in Sec. 3.3.

3.1 The fiducial run

Figure 3: The gas distribution on the scales of (a) and (b) just before turning on the radiation in the fiducial anisotropic radiation run. In each panel, the four quadrants (clockwise from top left) represent number density , temperature , neutral fraction of hydrogen and specific angular momentum shown by the corresponding centrifugal radius . The arrows represent the velocity vector , shown only when . The contours in the bottom left panel represent 0.5 (white), 0.6 (pink), 0.7 (orange), 0.8 (red) and 0.9 (dark red). The Bondi radii for neutral and ionized gases are shown as dashed black and white circles, respectively.

We first describe the fiducial run, where the parameter set is given as follows: , , , and . The radiation field is assumed to be anisotropic with shadow (see Fig. 1). Such wide obscuration can be regarded as an extreme of strong anisotropy. Starting from the homogeneous initial condition, we follow the evolution of flow with the radiation turned off until . We then turn on the radiation and follow the evolution until .

Figs. 2(a) and (b) show the time evolution of the accretion rate and the luminosity , respectively, along with the result for the case with isotropic radiation, which will be described in the next section. Fig. 2(a) shows that during the former period without radiation the accretion rate once converges to the constant value , only less than of the Bondi rate (). This remarkable reduction in is totally attributable to the effect of angular momentum.

Fig. 3 shows the gas distribution just before turning on the radiation. A flared disc formed inside can be seen. In the disc, the gravity is balanced with the centrifugal force inside and with the pressure gradient outside , so that the dynamical equilibrium () is approximately maintained throughout. In the bipolar regions, a gas with low angular momentum directly flows into the sink without hitting on the disc.

Before turning on the radiation, the gas is almost isothermal at due to the efficient Ly cooling throughout the computational domain (see Fig. 3). As an experiment, we have rerun the isothermal simulation setting and have confirmed that the result does not change. Therefore, in the rest of this paper, we adopt the isothermal equation of state with in cases without radiation, to save computational costs.

After turning on the radiation, the accretion rate decreases further, as seen in Fig. 2(a). It reaches the smaller constant value at , which is about 0.1 of the value before turning on the radiation and even less than of the Bondi rate. As shown in Sugimura et al. (2017), the accretion rate with this parameter set but without rotation is very high with . Therefore, this result demonstrates that the accretion rate is largely reduced by the interplay of angular momentum and radiation. The amount of the reduction will be analytically understood in Sec. 4. In Fig. 2(b), we see that the luminosity behaves in the same way as following Eq. (2). The luminosity is generally sub-Eddington, with at .

Figure 4: Same as Fig. 3 but at the end of the simulation, when the gas in the polar regions are photoionized by the anisotropic radiation. Here, the gas distribution is plotted on the scales of (a) , (b) and (c) .
Figure 5: Same as Fig. 3 but (a, b) before and (c, d) after an accretion burst in the isotropic radiation run.

Fig. 4 shows the flow structure at . The accretion occurs through the neutral disc remaining inside the shadow. The boundary between the ionized and neutral regions is determined by the shadow angle, from the equatorial plane. Outside the Bondi radius for the ionized gas, , material on the surface of the neutral gas is photoevaporated and flows out in the vertical direction, as the sound speed, , is larger than the escape velocity, . Conversely, inside , where , the photoionized gas falls back to the disc again. The outflows in the polar directions blow away the low- gases initially locating near the poles. Recall that accretion of such material boosted before the radiation is turned on.

In Fig. 4(c), a low density region appears near the poles as a result of the centrifugal barrier.3 However, structures along the -axis are partly artifacts of assumed axisymmetry (e.g., Sugimura et al., 2017). In addition, jets or outflows launched near the BH, if included, would largely change such features. Since these structures hardly affect anyway because of the small solid angle, we do not attempt to study them in more detail below.

Although the ionization boundary looks similar in shape to that in the non-rotating case, the flow structure is significantly altered by the angular momentum. Whereas the gas in the shadowed region falls freely in the non-rotating case, it has to pass through the accretion disc slowly in a viscous timescale otherwise. In the non-rotating case with large , Takeo et al. (2018) found that, for a massive enough () BH surrounded by a non-rotating medium, the radiation is obscured by the gas which is pushed by the ram pressure and intruded into the polar regions, thereby mitigating the radiation feedback. In the rotating case, however, the ram pressure of the flow is so weak that we do not see such a phenomenon. We have confirmed this by performing an additional simulation with enhanced to .

3.2 Isotropic radiation run

Next, we describe the isotropic radiation run. The set-up is the same as in the previous run, except that the radiation is isotropic.

Figs. 2(a) and (b) show the time evolution of and , respectively. Fig. 2(a) shows that after the radiation is turned on, oscillates violently repeating burst and quiescent phases. This behaviour is similar to what is found in the non-rotating case (e.g., Sugimura et al., 2017), but the averaged accretion rate in the last years, , is about ten times smaller. Here, we see again the reduction of by the interplay of angular momentum and radiation. Again, Fig. 2(b) shows that oscillates in the same way as following Eq. (2). Its absolute value is generally small and even at the burst phases.

Figs. 5 (a, b) and (c, d) show the gas distribution before and after an accretion burst, respectively. The H ii region contracts in a quiescent phase, and expands in a burst phase (e.g., Park & Ricotti, 2011). The neutral disc initially formed (Fig. 3) is completely ionized and vanishes soon after the radiation is turned on. The H ii regions is elongated along the -axis (Figs. 5a and c), because the density decreases near the axis.

3.3 Cases with different sets of the parameters

Figure 6: Same as Fig. 2(a) but for different sets of ). The grey lines are the results for the case with , which have been already shown in Fig. 2(a).

In order to investigate the dependence of the accretion flow on and , we perform a set of simulations assuming different values for these two parameters. While we have assumed and in the previous runs, we here change to or to . For each set of the parameters, we perform runs with both anisotropic and isotropic radiation fields.

Fig. 6 shows the evolution of in those runs (see also Table 1). We can see the parameter dependence of for each prescription of the radiation. First, before turning on the radiation, the value to which converges increases with increasing but is almost independent of . Next, in the anisotropic radiation case, at the end of the simulation increases with and slightly decreases with . The accretion rate in this case, and hence its parameter dependence, will be analytically understood in the next section. Finally, in the isotropic radiation cases, the average value of increases with but hardly depends on .

4 Analytical formulation

In the previous section, we found that the angular momentum can largely suppress the accretion in the case with anisotropic radiation, where the accretion rate would be high without angular momentum. In order to understand such suppression, we here develop an analytical model describing accretion through a neutral disc connected to a medium. This model is motivated by the fact that the accretion occurs through the neutral disc remaining inside the shadow in the anisotropic radiation runs (Fig. 4). Note that this model is applicable to neither the case with isotropic radiation where the disc is completely photoionized (Fig. 5), nor that without radiation where polar inflows of low angular momentum gas contribute to the accretion (Fig. 3). In the latter case, however, such polar inflows might be prevented by jets or outflows launched near the BH.

After developing the model, we derive the critical angular momentum of the medium needed to suppress the accretion, as the model predicts that the accretion rate goes back to the Bondi rate in the case with low angular momentum. Using the critical value, we can roughly estimate the impact of the angular momentum on the accretion rate form the property of the medium (e.g., Sec. 5.2).

Below, we first develop the analytical model in Sec. 4.1. Then, in Sec. 4.2, we interpret the simulation results obtained in Sec. 3 with the model. Furthermore, we derive the condition for accretion suppression in Sec. 4.3.

4.1 Analytical model of accretion through a disc connected to a medium

Here, to understand the mechanism of the suppression of accretion by the angular momentum, as well as to analytically estimate the extent of the suppression, we develop an analytical model for the accretion through a neutral disc connected to a rotating medium with constant specific angular momentum . In this model, we suppose that the accretion is suppressed because the gas stagnates at the centrifugal radius and accumulates between the centrifugal and the Bondi radii, and thus the pressure within the Bondi radius is enhanced.

We develop the model by connecting the following two types of solutions at the centrifugal radius . Outside , where the dynamical equilibrium is hold between the gravity, pressure gradient and centrifugal force, the accretion is associated with the inward gas supply that compensates the depletion to the inner disc at . Inside , where the Keplerian disc is formed, the accretion is caused by the viscous angular momentum loss, which transports the inward mass flux supplied at . Below, we assume that the gas is isothermal.

We start by obtaining the surface density for the outer dynamically equilibrium distribution that is connected to the medium with (Papaloizou & Pringle, 1984). Here, we assume that the gas has everywhere. Then, the equation for the force balance can be analytically solved, as shown in Appendix A.1. Inside , where the scale height is smaller than the radius and the distribution is disc-like, is given by (Appendix A.1)

(9)

Next, we obtain the relation between the accretion rate and for the inner viscous Keplerian disc (e.g., Shakura & Sunyaev, 1973; Kato et al., 1998). In Appendix A.2, assuming that the disc is thin and adopting the -type viscosity, , we obtain

(10)

In a steady accretion disc becomes radially constant by adjusting .

We then connect the two solutions at by substituting Eq. (9) into Eq. (10). We obtain

(11)

where we have used at . Here, if is so small that Eq. (11) yields , the assumption of the outer dynamically-equilibrium distribution breaks down, because the gas cannot be supplied to the disc at a rate exceeding . Thus, imposing that does not exceed , we modify Eq. (11) as

(12)

which gives an analytical estimate for the accretion rate from a rotating medium.

Figure 7: The parameter dependence of for the accretion from an isothermal medium with . We normalize by . The horizontal axis represents , while the colors correspond to (orange) and 0.01 (blue). The lines are the analytical estimates given by Eq. (12), whereas the filled dots are the numerical values obtained in Appendix B. We overplot the results of anisotropic radiation runs in Sec. 3 with open stars.

In Fig. 7, the parameter dependence of the analytical accretion rate is shown, along with the numerical results obtained in the equivalent settings (see Appendix B). In the figure, we normalize by . Note that how much is suppressed with respect to depends only on the combination of and (not on either , or ), because, with Eqs. (1) and (11), can be rewritten as

(13)

where we have used . The agreement of the analytical and numerical results is remarkable considering the simplicity of the model. In the case that the accretion is suppressed by the angular momentum, i.e., , the dependence is that is proportional to (Eq. 11). As for the dependence, rapidly increases with decreasing for , mainly because of the exponential increase of under the assumption of the dynamical equilibrium (see Eq. 9). When reaches , however, the dynamical equilibrium is broken and the increase of stops. For lower , is equal to , with little effect of the angular momentum. The disagreement between the numerical and analytical results can be partly attributed to the considerable thickness of disc at (see Fig. 10c), which is not consistent with our assumption of thin inner disc. As the thickness at increases with and the aspect ratio reaches unity when , our model is reliable only when is sufficiently smaller than .

4.2 Interpreting the simulation results with the analytical model

In this section, we understand the simulation results in Sec. 3 with the analytical model developed in the previous section, focusing mainly on the case with anisotropic radiation. We begin with overplotting the simulation results (star symbols) in Fig. 7, where the parameter dependence of in the model (lines), as well as that for the isothermal accretion (dot symbols), is shown. For all the three cases with examined in this paper, the simulation results are well reproduced by the model. The downward shifts of by at most a factor of three compared to the isothermal case (dot symbols) are partly attributable to the photoevaporation mass-loss from the surface of the disc. Note that we have also confirmed that the equatorial gas profile in the anisotropic radiation run with and is consistent with that of the model (Appendix C). The agreement both in the accretion rate and the gas profile suggests that the accretion is suppressed by the same mechanism as in the model, i.e., it is due to the stagnation of the gas at the centrifugal radius and the consequent pressure enhancement within the Bondi radius.

Below, we make some remarks on the simulation results in the case without radiation, as well as in the case with isotropic radiation. In the former case, the accretion rates are larger than in the model, because, in addition to the accretion through the disc, the polar low- gas directly flows into the sink and contributes to the accretion rate (see Fig. 3). Recall that in the anisotropic radiation runs, polar photoevaporative outflows prevent such inflows. Conversely, in the latter case, the accretion rates are smaller than in the model, because the accretion disc is totally photoevaporated (see Fig. 5). In the anisotropic radiation runs, however, the disc survives without being photoevaporated inside the shadow.

4.3 Condition for accretion suppression by the angular momentum

Figure 8: The critical centrifugal radius normalized by as a function of . The accretion is significantly suppressed by the angular momentum if .

Next, we obtain the critical angular momentum above which the accretion is significantly suppressed.

As mentioned in Sec. 4.1, we suppose that the effect of angular momentum becomes negligible when Eq. (13) yields . Thus, we define the critical centrifugal radius by the condition . In Fig. 8, we plot as a function of . The former slowly increases with the latter, with at . For a wide range of with , the suppression of accretion becomes important as exceeds , with decreased to when (see Fig. 7). Recall that the case with is beyond the scope of our analytical model (Sec. 4.1).

In our simulation runs in Sec. 3, we adopt and and thus they corresponds to the points far above the critical curve in Fig. 8. Therefore, suppression of the accretion by the angular momentum is expected according to our condition. In fact, we have observed significant reduction of the accretion rate in those runs, consistent with the expectation.

5 Discussion

5.1 Growth of Pop III remnant BHs

Rapid growth of BHs is critically important for the formation of SMBHs in the early universe, especially if their growth starts from “light seeds”, i.e., the Pop III remnant BHs. Reflecting the expected diversity of the stellar mass of the Pop III stars (e.g., Hirano et al., 2014, 2015), their remnant BHs would have a variety of masses, . The subsequent growth of such BHs should depend on the environments surrounding them (Fig. 9). In the standard bottom-up structure formation paradigm in the CDM universe, Pop III stars normally form in mini-halos with masses of at (e.g, Yoshida et al., 2006). In addition, Pop III stars also form in atomic-cooling halos with , if prior star formation is hampered by Lyman-Werner radiation and/or baryonic streaming motions (e.g., Wise & Abel, 2007b; Visbal et al., 2014; Tanaka & Li, 2014; Hirano et al., 2017; Schauer et al., 2017). In the following, we discuss the possibility of rapid growth of Pop III remnant BHs, separately considering their different formation sites of mini-halos and atomic-cooling halos. We argue that the BH growth via rapid mass accretion is not easily realized for both cases. In particular, the angular momentum effect studied in the previous sections is a critical obstacle for the growth of BHs formed in atomic-cooling halos.

Figure 9: Two evolutionary paths of a Pop III remnant BH and its environment. A Pop III remnant BH can form either in a minihalo (Sec. 5.1.1) or in an atomic-cooling (AC) halo (Sec. 5.1.2). See text for details.

Pop III remnants formed in mini-halos

Let us first consider the growth of Pop III remnant BHs in a mini-halo with . Massive Pop III stars with produce intense ionizing radiation before collapsing into BHs. This ionizing radiation feedback evacuates a large fraction of gas from the mini-halo, because of its shallow gravitational potential(Kitayama et al., 2004). Rapid mass accretion on to the remnant BH is thus not expected, at least when the BH is still harbored in the same mini-halo (Johnson & Bromm, 2007).

Through the assembly of DM halos, a significant fraction of Pop III remnant BHs fall into an atomic-cooling halo with at . In such a massive and gas-rich halo, the BHs could grow via accretion if they could sink to the galactic disc due to dynamical friction. We assume that a fraction of the gas in the halo forms stars, i.e., , and those stars exert friction on the remnant BHs. Here, we adopt (Planck Collaboration et al., 2016). If the density distribution of the stars is approximated by a singular isothermal sphere, the dynamical friction timescale for a BH4 is estimated as (Binney & Tremaine, 1987)

(14)

where the crossing time is and the Coulomb logarithm is set to . In the above expression, we have defined and . We assume that the BH is located at the outskirts of the galactic disc, i.e., (Mo et al., 1998), where is the dimensionless spin parameter and is the virial radius of the halo (Barkana & Loeb, 2001), with . Using this relation, the ratio of the dynamical friction timescale to the Hubble timescale is estimated as

(15)

with . This suggests that dynamical friction is inefficient for BHs with . As a result, Pop III remnant BHs coming originally from mini-halos will just continue to wander in the outskirts of the galactic disc. Such BHs will hardly grow via accretion, simply because the density of the surrounding gas is low (, also see §5.1.2 below). In this case, the rapid growth of BHs is unlikely to occur regardless of whether the angular momentum of the gas affects the BH feeding.

Pop III remnants formed in atomic-cooling halos

Pop III stars can be also formed in a galactic disc in an atomic-cooling halo with (or ). Such a halo can hold the gas against the stellar feedback (Kitayama et al., 2004). The remnant BHs are initially embedded in the gas-rich disc.

Let us suppose that the BH is embedded in an isothermal, exponential disc with a gas temperature . The gas density at the mid-plane within the disc radius () is estimated as (Oh & Haiman, 2002)

(16)

Assuming that a remnant BH has a peculiar velocity comparable to the circular velocity of the gas disc,

(17)

which is higher than the sound speed of the gas . Thus, the Bondi-Hoyle-Lyttleton (BHL) accretion rate is reduced by a factor of from the value in Eq. (1),

(18)

Note that Eq. (18) is valid for , and the BHL rate weakly depends on and . The typical BH growth timescale is . Thus, we obtain , which means that Pop III remnant BHs with formed inside an atomic-cooling halo may undergo rapid accretion. However, the angular momentum of the accreting gas here comes into play to prevent the BH growth.

According to cosmological simulations of the first galaxy, the gas flow in an atomic-cooling halo is turbulent in general (e.g., Wise & Abel, 2007a). We consider the same density fluctuation of the turbulent medium as in our Galaxy, , where is the characteristic spacial length of the fluctuation (Armstrong et al., 1995; Draine, 2011). Accreting gas with such density fluctuation brings a net angular momentum, which is estimated as (Ipser & Price, 1977; Ioka et al., 2017; Matsumoto et al., 2018). Thus, the ratio of the centrifugal radius to the Bondi radius is given by

(19)

This value is comparable to the critical one of obtained in Sec. 4.3, above which the accretion rate is suppressed by a factor of . As a result, even if Pop III remnant BHs is embedded in the gas disc, the BH growth timescale becomes even longer than the age of the universe at ,

(20)

where . Therefore, we conclude that Pop III remnant BHs with are unlikely to experience rapid growth via accretion.

Possible pathways for rapid growth of Pop III remnant BHs

A number of previous studies (e.g., Madau et al., 2014; Alexander & Natarajan, 2014; Volonteri et al., 2015; Tagawa et al., 2015; Ryu et al., 2016) have investigated the possibility of rapid (super-Eddington, i.e., and thus ) gas accretion of Pop III remnant BHs to explain the existence of high-z SMBHs. Since our speculation seems different from those works, we discuss what kind of uncertainties lead to the discrepancies, comparing to our arguments in §5.1.1 and 5.1.2.

Tagawa et al. (2015) and Ryu et al. (2016) considered the evolution of Pop III remnant BHs in an atomic-cooing halo, performing N-body simulations. In fact, they concluded that some remnant BHs can fall into the central region much faster than we estimated. This is mainly because the gas density is assumed to be higher than or to follow an isothermal singular profile (). Because of the higher gas densities (at smaller scales), dynamical friction allows the BHs to quickly sink into the galactic center. When the remnant BH reaches the central region with , the Bondi accretion rate on to it becomes high enough () to realize hyper-Eddington accretion without impeded by radiation feedback (Inayoshi et al., 2016; Sakurai et al., 2016; Sugimura et al., 2017; Takeo et al., 2018). This process may quickly form intermediate massive BHs with at the center of the protogalaxy.

Obviously, a key uncertainty is the density structure of the gas containing the BHs. Tagawa et al. (2015) and Ryu et al. (2016) assume that the density profile does not change during the orbital evolution of BHs. In reality, however, the star formation will easily occur in the high density regions, so that the original density structure could be modified by feedback effects such as supernova explosions (e.g., Dubois et al., 2015; Yajima et al., 2017). In future work, we will study the evolution of the BH orbital motions and ambient density structure, self-consistently incorporating the feedback effects. Such treatment will also allow us to accurately estimate how much angular momentum is brought by the accreting gas, and to assess whether super-Eddington accretion is possible circumventing the angular momentum barrier presented in this paper.

5.2 Caveats

We have made a number of simplifications and approximations in this work, and now discuss their significances.

First, we have examined only the two limiting cases of the anisotropy of radiation field. In the case of a non-rotating medium, Sugimura et al. (2017) have found that there is a critical shadowing angle ( from the equatorial plane) above which the efficient accretion is realized by the neutral Bondi-like inflows through the equatorial layer that exceeds the photoevaporative mass-loss from the surfaces. In the rotating case, however, the above condition needs to be modified, because the photoevaporative mass-loss from the surfaces of the rotationally-supported disc can be significant. We have seen that the disc is completely ionized by the isotropic radiation, while it is not photoevaporated inside the shadow in the case of the anisotropic radiation. We will study the dependence of on the anisotropy more in the future.

Second, the actual anisotropy created inside the sink is highly uncertain, although strongly depends on it. In the literature, generation of (failed) winds or coronae above the disc have been investigated (e.g., Hollenbach et al., 1994; Begelman et al., 1983; Woods et al., 1996; Proga et al., 2000; Wada, 2012; Suzuki & Inutsuka, 2014; Nomura et al., 2016). We expect that the materials associated with such structures obscure the outward radiation and create its anisotropy. In a future work, we plan to study such process, considering its dependence on , , the metallicity of gas, etc..

Third, we have adopted the -type viscosity to mimic the angular momentum transport via the turbulence driven by the MRI. Although our results depend on the value of , as well as where the viscosity works, i.e., the confinement factor in Eq. (8), it is computationally too expensive to perform 3D magnetohydrodynamics simulations of the same problem. The unstable non-axisymmetric modes can also affect the flow in the 3D simulations (Papaloizou & Pringle, 1984). In the cases studied here, the Toomre parameter is above unity and the disc is gravitationally stable. In a case with different parameter set, however, the gravitational instability can play a role in transporting the angular momentum depending on and (see Appendix E).

Fourth, there is an uncertainty regarding the accumulation of the angular momentum outside the outer edge of the disc. Although it does not affect if the medium has moderate angular momentum, as studied in Sec. 3, it can largely reduce in the case with small angular momentum (; see Appendix B). While the accumulation occurs in a highly symmetric system, it becomes irrelevant if the direction of the (net) angular momentum of the accreted gas varies before the accumulation proceeds, e.g., in the case of accretion from a turbulent medium. To resolve the uncertainty, we need to consider realistic environments, such as the ones based on galactic-scale simulations. In addition, the disc is known to be Rayleigh unstable when the angular momentum decreases outward, so that the accumulated angular momentum would be transported in 3D simulations (see, e.g., Inayoshi et al., 2018, and reference therein).

Furthermore, although we assume that the dominant cooling process in the neutral gas is the Ly cooling, the free-bound cooling becomes dominant and cools the gas to when (Omukai, 2001). In addition, the temperature might drop even to if molecules somehow form in spite of the UV irradiation from the BH neighbourhood. Consideration of these processes could lead to modification of .

6 Summary

We have investigated the combined effect of gas angular momentum and radiation feedback on seed BH accretion, by preforming a suit of 2D axisymmetric simulations considering both finite gas angular momentum and radiation from the circum-BH disc. The BH is located at the center of a rotating medium, whose centrifugal radius is typically a tenth of the Bondi radius. We follow the formation of the rotationally-supported disc, through which the accretion proceeds by the angular momentum transport due to the assumed -type viscosity.

We have found that the accretion is strongly suppressed by the gas angular momentum. The accretion rate is reduced by one order of magnitude even without radiation feedback, but it becomes even smaller with the feedback. In particular, the accretion rate in the case with anisotropic radiation field, which would be in the same order as the Bondi rate without gas rotation (Sugimura et al., 2017), is reduced by a factor of . Our results clearly indicate the importance of the interplay of the angular momentum and radiation feedback.

We have also developed an analytical model that describes accretion through a neutral disc connected to a medium. This model is capable of reproducing the accretion rate obtained in the simulations with anisotropic radiation. Furthermore, the model suggests the presence of the critical angular momentum above which the accretion is significantly suppressed. The corresponding critical centrifugal radius is times the Bondi radius for a wide range of , suggesting that even such small angular momentum is enough to reduce the accretion rate. This provides a useful estimate for the impact of the angular momentum on the accretion rate, for example, in future cosmological simulations of SMBH formation.

Finally, we have discussed the implications of our findings on the growth of Pop III remnant BHs. In the literature, those BHs are claimed to grow to high- SMBHs by very rapid (super-Eddington) accretion. However, the angular momentum effect studied in this paper can be a crucial obstacle for such rapid mass accretion. Whilst the condition for the formation of direct-collapse BHs is known to be difficult to achieve (e.g., Sugimura et al., 2014), it should be equally challenging for the Pop III remnant BHs to rapidly grow via the super-Eddington accretion. Clearly, further studies are needed to unveil the nature.

Acknowledgements

The authors would like to thank Ken Ohsuga, 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), 17H04827 (HY) and 17H01102 and 17H06360 (KO) and by the Simons Foundation through the Simons Society of Fellows (KI).

Appendix A Inner and outer solutions of analytical model

Our model consists of an outer dynamically equilibrium distribution (Sec. A.1) and an inner viscous Keplerian disc (Sec. A.2), which are connected at the centrifugal radius. Below, we describe the outer and inner solutions in this order.

a.1 Outer dynamically equilibrium distribution

Here, we describe an outer dynamically equilibrium distribution connected to a homogeneous medium (Papaloizou & Pringle, 1984). We assume that the gas is isothermal and that the angular momentum everywhere.

We start by describing the equation for the dynamical equilibrium between the gravity, pressure gradient and centrifugal force,

(21)

where and are the unit vectors in the and directions, respectively. Using the isothermal equation of state, , Eq. (21) can be rewritten as

(22)

Then, imposing the outer boundary condition of constant density except near the pole, i.e., as with finite , we obtain

(23)

where we have used and in the second equality. The density rapidly decreases towards the pole, i.e., as , due to the assumed constant angular momentum.

The density profile corresponds to that of a thin disc at . Near the equatorial plane, where , we can approximate Eq. (23) as

(24)

with the equatorial density,

(25)

and the scale height,

(26)

At , where the aspect ratio is small, the majority of gas is confined to a layer with . In such layer, the approximate form of Eq. (24) is valid and the integration in the -direction yields the expression for surface density,

(27)

which can be rewritten as Eq. (9) with Eq. (24). Note that has its maximum at , as the pressure gradient is balanced with the gravity at but with the centrifugal force at .

a.2 Inner viscous Keplerian disc

Next, we describe a thin Keplerian disc extended towards the BH (e.g., Shakura & Sunyaev, 1973; Kato et al., 1998). Through the disc, the steady accretion is driven by the angular momentum transport via the -type viscosity. Again, we assume the isothermal gas.

The gas distribution is governed by the following equations. Since the disc is vertically hydrostatic, the -dependence of density is the same as Eq. (24) and thus the relation between and is given by Eq. (27) (though is different from Eq. 25). Radially, the gravity is balanced with the centrifugal force with angular velocity

(28)

The vertically-integrated mass conservation equation can be written as

(29)

with the accretion rate . Finally, the vertically-integrated angular momentum conservation equation is given by

(30)

with the net angular momentum flux . Note that we have used Eqs. (28) and (29) to obtain Eq. (30). Finally, we adopt the -type viscosity

(31)

which is the same as Eq. (8) with .

The structure of the disc is uniquely determined once , , and are given. Here, can be determined from the inner boundary condition: we impose the torque free condition, i.e., the second viscous stress term in Eq. (30) is set to zero, at the inner boundary where the first advection term, which is proportional to , is also small.5 Thus, neglecting in the right-hand side of Eq. (30), we finally obtain (Eq. 10). For given , with Eq. (10), as well as Eqs. (27),  (