Transition of BH feeding from the quiescent regime into star-forming cold disk regime

TL;DR

This study uses hydrodynamical simulations to explore how supermassive black holes transition from low-efficiency accretion to cold disk feeding, revealing conditions that trigger star formation and explain observed luminosity discrepancies.

Contribution

It provides a detailed simulation-based analysis of the transition from radiatively inefficient to cold disk accretion in SMBHs, linking accretion rates to star formation and luminosity.

Findings

Low Bondi rates suppress BH feeding and luminosity.

Higher Bondi rates lead to cold disk formation and increased luminosity.

Cold disks become gravitationally unstable, promoting star formation.

Abstract

We study the properties of rotating accretion flows onto supermassive black holes (SMBHs) using axisymmetric two-dimensional hydrodynamical simulations with radiative cooling and BH feedback. The simulations resolve the accretion dynamics of gas outside from the BH influence radius through an inner accretion disk. For lower Bondi accretion rates in units of the Eddington rate ($\dot{M}_{\rm B}\ll 10^{-3}~\dot{M}_{\rm Edd}$), the BH feeding is suppressed due to turbulent motion by several orders of magnitudes from the Bondi rate. Thus, the radiative luminosity results in as low as $\sim 10^{-10}-10^{-7}~L_{\rm Edd}$, where $L_{\rm Edd}$ is the Eddington luminosity. For higher rates of $\dot{M}_{\rm B}> 10^{-3}~\dot{M}_{\rm Edd}$, the optically-thin accreting gas cools via free-free emission and forms a geometrically-thin disk, which feeds the BH efficiently and increases the radiative luminosity to $> 10^{-3}~L_{\rm Edd}$. The transitional behavior of accreting BHs in galactic nuclei from radiatively inefficient phases to cold disk accretion naturally explains (1) the reason for the offset between the observed luminosities and theoretical predictions for nearby quiescent SMBHs, and (2) the conditions to fuel gas into the nuclear SMBH. In addition, the cold disk formed in galactic nuclei tends to be gravitationally unstable and leads to star formation when the Bondi rate is as high as $ \dot{M}_{\rm B} > 10^{-2}~M_\odot~{\rm yr}^{-1}$. This is a plausible explanation of the correlation observed between star formation rates and BH feeding rates in Seyfert galaxies.

Equations

$$R_{\mathrm{B}}\equiv \frac{G{M}_{\bullet }}{c_s^2},$$

$$R_{\mathrm{C}}\equiv \frac{j^2}{G{M}_{\bullet }},$$

$$R_{\mathrm{Sch}}\equiv \frac{2G{M}_{\bullet }}{c^2}\ll R_{\mathrm{C}},R_{\mathrm{B}},$$

$$\frac{d\rho}{dt}+\rho\nabla\cdot\mbox{\boldmath$v$}=0,$$

$$\rho\frac{d\mbox{\boldmath$v$}}{dt}=-\nabla p-\rho\nabla\Phi+\nabla\cdot\mbox{\boldmath$\sigma$}+f_{\rm rad}\mbox{\boldmath$e_{r}$},$$

$$\rho\frac{de}{dt}=-p\nabla\cdot\mbox{\boldmath$v$}+(\mbox{\boldmath$\sigma$}\cdot\nabla)\mbox{\boldmath$v$}-\Lambda+\Gamma,$$

$$\sigma_{ij}=\rho\nu\left[\left(\frac{\partial v_{j}}{\partial x_{i}}+\frac{\partial v_{i}}{\partial x_{j}}\right)-\frac{2}{3}(\nabla\cdot\mbox{\boldmath$v$})\delta_{ij}\right],$$

$$\nu =\alpha \frac{c_{\mathrm{s}}^2}{{\Omega }_{\mathrm{K}}},$$

$$\alpha ={\alpha }_0\left\{\exp \left[-{\left(\frac{{\rho }_{\mathrm{cr}}}{\rho }\right)}^2\right]+\max \left(0,-\frac{\partial \ln j}{\partial \ln r}\right)\right\},$$

$$\Gamma -\Lambda =n^2(S_{\mathrm{br}}+S_{\mathrm{Comp}}+S_{\mathrm{ph}}),$$

$$S_{\mathrm{Comp}}=\frac{k_{\mathrm{B}}(T_{\mathrm{C}}-T)}{m_{\mathrm{e}}{c}^2}\frac{{\sigma }_{\mathrm{es}}{L}_{\mathrm{bol}}}{\pi n{r}^2},$$

$$f_{\mathrm{rad}}=\frac{\rho {\kappa }_{\mathrm{es}}}{c}\frac{L_{\mathrm{bol}}}{4\pi r^2}+\frac{\Gamma }{c}.$$

$$\log ϵ=\{\begin{array}{ll}-1.0-(0.0162/\dot{m}{)}^4& \mathrm{for}0.023\le \dot{m},\text{vspace}1mm\\ {\sum }_n{a}_n(\log \dot{m}{)}^n& \mathrm{for}10^{-4}<\dot{m}<0.023,\text{vspace}1mm\\ {\sum }_n{b}_n(\log \dot{m}{)}^n& \mathrm{for}10^{-8}<\dot{m}\le 10^{-4},\end{array}$$

$$\rho ={\rho }_{\infty }{\left[1+(\gamma -1)\frac{G{M}_{\bullet }}{c_{\infty }^{2}r}-\frac{(\gamma -1)}{2}\frac{j_{\infty }^2}{c_{\infty }^2{\varpi }^2}\right]}^{1/(\gamma -1)}$$

$${\rho }_{\mathrm{cr}}={\rho }_{\infty }{\left(1+\frac{\gamma -1}{2\beta }\right)}^{1/(\gamma -1)},$$

$$\beta =\frac{R_{\mathrm{C}}}{R_{\mathrm{B}}}=\frac{j_{\infty }^2{c}_{\infty }^2}{G^2{M}_{\bullet }^2}.$$

$${\dot{M}}_{\mathrm{B}}$$

$$\simeq 8.5\times 10^{-3}{\rho }_{-22}{M}_8^2{T}_7^{-3/2}{\mathrm{M}}_{\odot }{\mathrm{yr}}^{-1},$$

$${\dot{m}}_{\mathrm{B}}\equiv \frac{{\dot{M}}_{\mathrm{B}}}{{\dot{M}}_{\mathrm{Edd}}}\simeq 3.7\times 10^{-3}{\rho }_{-22}{M}_8{T}_7^{-3/2},$$

$${\dot{M}}_{\mathrm{in}}$$

$${\dot{M}}_{\mathrm{out}}$$

$$\rho$$

$$=\frac{e^{3/2}}{6\alpha }{\rho }_{\infty }{\left(\frac{c_{\infty }}{c_{\mathrm{s}}}\right)}^3{\left(\frac{r}{R_{\mathrm{B}}}\right)}^{-3},$$

$$\frac{t_{\mathrm{th},\mathrm{vis}}}{t_{\mathrm{cool}}}\simeq 1.48T_7^{1/2}{\left(\frac{\alpha }{0.01}\right)}^{-1}\left(\frac{{\dot{m}}_{\mathrm{B}}}{10^{-3}}\right){\left(\frac{\beta }{0.1}\right)}^{1/2},$$

$$\frac{{\dot{M}}_{\mathrm{B}}}{{\dot{M}}_{\mathrm{Edd}}}\text{ga}{\dot{m}}_{\mathrm{cool}}\equiv 6.8\times 10^{-4}{T}_7^{-1/2}\left(\frac{\alpha }{0.01}\right){\left(\frac{\beta }{0.1}\right)}^{-1/2}.$$

$$Q\equiv \frac{c_{\mathrm{s}}\Omega }{\pi G\Sigma }.$$

$$Q$$

$$\simeq 2.1T_{\min ,4}^{3/2}{M}_8^{-1}\left(\frac{\alpha }{0.01}\right){\left(\frac{{\dot{m}}_{\mathrm{B}}}{10^{-3}}\right)}^{-1},$$

$${\alpha }_{\mathrm{eff}}={\alpha }_0+{\alpha }_{\max }\exp (-Q^A/B),$$

$$\frac{{\dot{M}}_{\mathrm{B}}}{{\dot{M}}_{\mathrm{Edd}}}\text{ga}5.6\times 10^{-3}{T}_{\min ,4}^{3/2}{M}_8^{-1},$$

$${\dot{M}}_{\mathrm{B}}\text{ga}1.3\times 10^{-2}{T}_{\min ,4}^{3/2}{\mathrm{M}}_{\odot }{\mathrm{yr}}^{-1},$$

$$\Sigma \simeq \frac{c_{\mathrm{s}}\Omega }{\pi GQ}\simeq 2.4\times 10^4\mathrm{g}\mathrm{c}{\mathrm{m}}^{-2}{T}_{\min ,4}^{1/2}{M}_9^{-1}{Q}^{-1}{\hat{r}}_3^{-3/2},$$

Full text

1 Introduction

Supermassive black holes (SMBHs) are almost ubiquitously harbored at the centers of massive nearby galaxies. The existence of SMBHs is consistent with the number and energetics of high-redshift bright quasars (QSOs), which are associated with efficient gas accretion onto SMBHs (Soltan, 1982; Yu & Tremaine, 2002). Through BH feeding and feedback processes, SMBHs are believed to coevolve with their host galaxies over the cosmic time (e.g., Kormendy & Ho, 2013, and references therein). In the local universe, however only a few percent of SMBHs are observed as luminous active galactic nuclei (AGN). A majority of them are nearly quiescent and known as low luminosity AGN (Ho, 2008, 2009), or they may be undetectably faint.

While low-luminosity AGN are energetically unimpressive, the study of such AGN phenomena brings us important insights and intriguing questions regarding the physics of BH accretion. One of the puzzles is that nearby, silent SMBHs exhibit levels of activities much lower than that expected from gas supplying rates onto the galactic nuclei (e.g., Pellegrini, 2005; Ho, 2008, 2009). This luminosity deficit problem clearly suggests the lack of our understanding on how accreting gas finally reaches to the central SMBHs. Simultaneously, we need to address how such quiescent SMBHs turn into active phases, and how their host galaxies control the cycle of the two phases.

The basic physics of the accretion process was first studied for a spherically symmetric flow onto a BH with a mass of $M_{\bullet}$ without radiative effects (Bondi & Hoyle, 1944; Bondi, 1952). This approach remains a good approximation only if rotation is so minor as to allow quasi-spherical accretion all the way down to the BH from the so-called Bondi radius defined by

[TABLE]

within which the negative gravitational energy dominates the thermal energy of the gas with a sound speed of $c_{s}$, where $G$ is the gravitational constant. However, it is hard to conceive of cases where the gas could have such a tiny amount of angular momentum. The characteristic radius where the centrifugal force balances the gravity of the BH is given by

[TABLE]

where $j$ is the angular momentum per unit mass. This radius may or may not be larger than the Bondi radius, but is always larger than the Schwarzschild radius,

[TABLE]

where $c$ is the speed of light. This means that the inflowing gas must always fall to a rotating (thin or thick) disk before reaching the BH, and subsequently its mass transport inwards is determined by the physical processes that transport the angular momentum outwards (Lynden-Bell & Pringle, 1974; Pringle, 1981).

At sufficiently low accretion rates, the infalling matter is not dense enough to cool. Thus, the flow is essentially adiabatic, and then there would be no net accretion in the absence of viscosity. However, magneto-hydrodynamical (MHD) simulations have shown that even small entrained magnetic fields will induce the magneto-rotational instability (MRI) in a sufficiently ionized disk (Balbus & Hawley, 1991; Matsumoto & Tajima, 1995; Stone et al., 1996; Balbus & Hawley, 1998; Stone & Pringle, 2001; Hawley et al., 2001; Machida et al., 2001; McKinney & Gammie, 2004; Ohsuga et al., 2009; Bai, 2011; Narayan et al., 2012; Suzuki & Inutsuka, 2014). MRI-driven turbulence causes an effective shear viscosity which permits further inflows and is often described with the $\alpha$-viscosity prescription (e.g., Shakura & Sunyaev, 1973). In this very low accretion rate domain, several papers have shown that rotating accretion flows become convectively unstable and have been termed convection dominated accretion flows (CDAFs; Igumenshchev & Abramowicz, 1999; Stone et al., 1999; Narayan et al., 2000; Quataert & Gruzinov, 2000; Igumenshchev & Abramowicz, 2000).

Although other solutions of hot accretion flows have been proposed by many authors for different boundary conditions (e.g., Ichimaru, 1977; Narayan & Yi, 1994, 1995; Blandford & Begelman, 1999, 2004), accreting matter captured by the BH gravity through the Bondi radius preferentially settles down to a CDAF within the centrifugal radius (Inayoshi et al. 2018, Paper I). In the CDAF solution, convective flows form a geometrically-thick disk, and the inflow rate decreases toward the center as $\dot{M}_{\rm in}\propto r$ within the disk. Physically, the reason for this is that the flow is not strongly bound to the central BH in the absence of cooling, so that at each radius gravitational energy release is sufficient to reverse the flow of some material and consequently to reduce the net accretion rate. As a result, the net accretion rate onto the central BH is extremely small compared to the rate measured at the Bondi radius. As shown in a schematic overview in Fig. 1, the solution describes the situation in the majority of the locally observed BHs such as the ones in the Milky Way, M31 up to the giant one in M87 (Inayoshi et al., 2018), illustrated by the pink band in Fig. 1.

Now, let us consider higher density environments which could be occasioned either by cooling flows within galaxies (e.g., Ciotti & Ostriker, 2001; Dekel & Birnboim, 2006; Gaspari et al., 2013, 2015; Gan et al., 2018) or galaxy mergers (e.g., Sanders et al., 1988; Juneau et al., 2009; Hopkins & Quataert, 2010). At some critical density, the rotating disk cools and collapses to a thin disc and then, as we shall see, in most cases the disks are locally unstable to the gravitational instabilities (Toomre, 1964; Gammie, 2001). Then spiral modes develop and provide the angular momentum and mass transport (e.g., Thompson et al., 2005; Rice et al., 2005; Lodato & Natarajan, 2006; Hopkins & Quataert, 2010; Kuiper et al., 2011; Meyer et al., 2017; Gan et al., 2018). There is both observational and theoretical evidence for massive star formation to occur within the gravitationally unstable disks (Goodman, 2003; Levin & Beloborodov, 2003; Levin, 2007; Nayakshin et al., 2007).

In this paper, we aim to pinpoint just where and when accretion onto BHs can change drastically from almost being unobservably faint to emitting more than a billion times the solar luminosity. For this purpose, we perform axisymmetric two-dimensional hydrodynamic simulations including a variety of physical processes, e.g., radiative cooling and feedback associated with BH accretion. We investigate the gas dynamics of accretion in the cooling transition and elucidate the characteristics of the flow in the gravitationally unstable domain. It occurs at the level where the accretion rate from the Bondi radius is $\dot{m}_{\rm cool}\sim O(10^{-3})$ measured in units of the Eddington accretion rate (see the yellow shaded region in Fig. 1). There is another flow boundary at still higher inflow rates of $\dot{m}_{\rm fb}\sim O(10^{-2})$, where the radiative feedback from the central BH, i.e., gas heating at larger scales around $R_{\rm B}$, causes the flow to reverse and become intermittent. This critical accretion rate has also been investigated by many authors (Ciotti & Ostriker, 2001; Proga, 2007; Ciotti et al., 2009; Milosavljević et al., 2009; Park & Ricotti, 2011; Novak et al., 2011; Choi et al., 2014; Inayoshi et al., 2016) that basically focused on the feedback effect on the gas dynamics at $r\la R_{\rm B}$ or even larger galactic scales $\sim O({\rm kpc})$. At even higher accretion rates of $\dot{m}\ga 1$, emergent photons are efficiently trapped within optically thick accretion flows (e.g., Ohsuga et al., 2005; Jiang et al., 2014; Sa̧dowski et al., 2015), and thus radiative feedback hardly affects gas supplied from the Bondi scales when the rate is as high as $\dot{m}\sim O(10^{3})$ (Inayoshi et al., 2016; Takeo et al., 2018, 2019). We do not reach those levels of accretion in this work, but will study the feedback self-consistently coupled with the dynamics of the accretion disk in future work.

We here focus on the accretion dynamics at larger scales covering the Bondi radii, unlike previous work that has investigated the accretion flows at smaller scales assuming a compact torus in hydrostatic equilibrium as the initial state. This basic spirit of our simulations allows us to consider plausible initial and boundary conditions, some of which can be directly measured in the nuclear regions surrounding quiescent SMBHs. As shown in Fig. 1, we are also able to test our theoretical model compared to the energy output (e.g., radiative luminosity) associated with actual BH feeding at the centers of nearby galaxies.

The rest of this paper is organized as follows. In §2, we describe the methodology of our numerical simulations. In §3, we show our simulation results and explain their physical properties. In §4.1, we discuss the conditions for the transition from the radiatively inefficient, hot, torus-like accretion domain to radiatively efficient, cold, thin disk accretion domain. In §4.2, we discuss the possibility of gravitational instability of the disk and the properties of marginally unstable disk structure. In §4.3, we describe the relation between the Bondi accretion rates and radiative luminosities and explain the reason for the luminosity deficit of SMBHs in the local Universe. Moreover, in §4.4 we argue the AGN-starburst connection, i.e., the relation between BH feeding rates and star formation rates in the circumnuclear regions. Finally, we summarize our main conclusions in §5.

2 Methodology

2.1 Basic equations

We solve the axisymmetric two-dimensional hydrodynamical equations using the open source code PLUTO (Mignone et al., 2007) with modifications described in Kuiper et al. (2010, 2011). The basic equations are following: the equation of continuity,

[TABLE]

and the equation of motion,

[TABLE]

where $\rho$ is the density, $v$ is the velocity, and $p$ is the gas pressure, the gravitational potential is set to $\Phi=-GM_{\bullet}/r$, $r$ is the distance from the central BH, $\sigma$ is the stress tensor due to viscosity, and $f_{\rm rad}$ is the outward net radiation force in the radial direction. The time derivative is the Lagrangian derivative, given by $d/dt\equiv\partial/\partial t$ + v$$\cdot\nabla.

We solve the energy equation of

[TABLE]

where $e$ is the internal energy per mass. The equation of state of the ideal gas is assumed as $p=(\gamma-1)\rho e$, where $\gamma=5/3$. The first and second terms on the right-hand-side present compressional heating (or expansion cooling) and viscous heating. The last two terms are radiative cooling and heating (in units of ${\rm erg~{}s^{-1}~{}cm^{-3}}$).

The viscous stress tensor is given by

[TABLE]

where $\nu$ is the shear viscosity. Note that the bulk viscosity is neglected here. The shear viscosity is calculated with the $\alpha$-prescription (Shakura & Sunyaev, 1973),

[TABLE]

where $\Omega_{\rm K}\equiv(GM_{\bullet}/r^{3})^{1/2}$. We calculate the viscous parameter by mimicking some properties of the MRI as

[TABLE]

where the strength of viscosity is set to $\alpha_{0}=0.01$ according to MHD simulations of the global disk structure (e.g., Zhu & Stone, 2018; Takasao et al., 2018), and $\rho_{\rm cr}$ is a threshold of the density above which viscosity turns on. We adopt a maximum value of the density at an initial condition as the threshold value (see §2.2). Under this model, the viscous process is active primarily in an accretion disk ($r\la R_{\rm C}$), where the rotational velocity has a significant fraction of the Keplerian velocity. On the other hand, angular momentum transported from the disk would be accumulated outside it, where no viscous processes operate. Such rotating flows with negative gradients of the specific angular momentum outward, i.e., $\partial j/\partial r<0$, are unstable and become turbulent, leading to angular momentum transport in three-dimensional (3D) simulations (Chandrasekhar, 1961). Since our simulations do not capture this 3D effect, instead we add the second term in Eq. (9) so that a steady state of the accretion flow exists. We note that the treatment of the rotational instability does not affect our results (see Appendix in Inayoshi et al. 2018).

In our simulations, radiative processes are taken into account in what follows. The radiative heating and cooling terms are expressed as

[TABLE]

where each term corresponds to the rate associated with bremsstrahlung cooling ($S_{\rm br}$), Compton heating/cooling ($S_{\rm Comp}$), and the sum of photoionization heating, line and recombination continuum cooling ($S_{\rm ph}$). The Compton heating/cooling rate is given by

[TABLE]

where $L_{\rm bol}$ is the AGN bolometric luminosity, $T_{\rm C}$ is the Compton temperature, $n$ is the number density of electron, $k_{\rm B}$ is the Boltzmann constant, $m_{\rm e}$ is the electron mass, $c$ is the speed of light, and $\sigma_{\rm es}$ is the cross section of electron scattering. We adopt $T_{\rm C}=10^{8}~{}{\rm K}$ as our fiducial value for observed low-luminosity AGN (Xie, Yuan & Ho, 2017). Note that Compton temperature for low-luminosity AGN ($L_{\rm X}/L_{\rm Edd}<10^{-3}$) is higher than that estimated in previous works using composite spectra of high luminosity AGNs; $T_{\rm C}\simeq 2\times 10^{7}~{}{\rm K}$ (Sazonov et al., 2004). The functional forms of $S_{\rm br}$ and $S_{\rm ph}$ are given in Sazonov et al. (2005). In this paper, we consider situations where accretion flows are optically thin. Therefore, the outward radiation force term in equation of motion can be expressed by

[TABLE]

We estimate the bolometric AGN luminosity as $L_{\rm bol}=\epsilon\dot{M}_{\bullet}c^{2}$, where $\dot{M}_{\bullet}$ is the gas accretion rate through the inner-most grid and $\epsilon$ is the radiative efficiency.

[TABLE]

where $\dot{m}\equiv\dot{M}_{\bullet}/\dot{M}_{\rm Edd}$. The fitted values are $a_{0}=-0.807$, $a_{1}=0.27$, $a_{n}=0$ ($n\geq 2$), $b_{0}=-1.749$, $b_{1}=-0.267$, $b_{2}=-0.07492$ and $b_{n}=0$ ($n\geq 3$). In Fig. 2, we show the radiative efficiency model we adopt (blue solid curve), which consists of the following two models. The radiative efficiency at a higher rate of $\dot{m}>10^{-4}$ is obtained by a semi-analytical model (green dashed curve; Xie & Yuan, 2012). For lower rates of $\dot{m}\leq 10^{-4}$, the fitted values are based on the result of general relativity and radiation transfer simulations (orange squares; see Table 1 in Ryan et al. 2017), where gas accretion onto a Kerr BH with a spin of $a_{\bullet}=0.5$ has been studied. Our model for the radiative efficiency is significantly higher than that for an ADAF-like model (e.g., Ciotti et al., 2009, red dotted), which is often used in previous work, but does not take into account energy transfer by Coulomb collision between ions and electrons properly.

2.2 Boundary and initial conditions

To compute the basic equations, we employ spherical coordinates in a computational domain of $r_{\rm in}\leq r\leq r_{\rm out}$ and $\epsilon\leq\theta\leq\pi-\epsilon$, where $\epsilon$ is set to $0.01$ to avoid numerical singularity at poles. In the radial direction, we set up logarithmically-spaced grids, the number of which is $N_{r}=512$. In the polar direction, we use static refinement in order to resolve the scale height of a cold disk, setting two types of uniformly spaced grids; $\Delta\theta=3(\pi-2\theta_{0})/N_{\theta}$ at $\theta_{0}\leq\theta\leq\pi-\theta_{0}$ (near the mid-plane) and $\Delta\theta=3(\theta_{0}-\epsilon)/N_{\theta}$ otherwise (near the poles). We set $N_{\theta}=384~{}(=128\times 3)$ and $\theta_{0}=1.3$ so that the disk scale height is sufficiently resolved. As our fiducial case, we set $r_{\rm in}=1.7\times 10^{-2}~{}R_{\rm B}$ and $r_{\rm max}=33~{}R_{\rm B}$. The grid structure of our simulations is illustrated in Fig. 3.

As initial conditions, we adopt a rotational equilibrium distribution ($v_{r}=v_{\theta}=0$) with a constant specific angular momentum of $j_{\infty}$:

[TABLE]

(Fishbone & Moncrief, 1976; Papaloizou & Pringle, 1984), where $\varpi=r\sin\theta$ is the cylindrical radius, $\rho_{\infty}$ is the ambient gas density, and $c_{\infty}$ is the sound speed of the gas. The first and second terms on the right-hand side present density enhancement via gravity of the central BH inside the Bondi radius. The third term expresses the centrifugal force for a given $j_{\infty}$, which leads to a maximum value of the density at $r=R_{\rm C}[=j^{2}_{\infty}/(GM_{\bullet})]$ and $\theta=\pi/2$,

[TABLE]

where $\beta$ is defined by the ratio of the centrifugal radius and Bondi radius,

[TABLE]

Note that the corresponding specific angular momentum is given by $j_{\infty}=\sqrt{\beta}R_{\rm B}c_{\infty}$. Without viscosity, the density never exceeds this value because of the centrifugal barrier. In other words, a high-density region with $\rho>\rho_{\rm cr}$ must be formed by angular momentum transport due to viscosity. In Fig. 3, the density distribution of the equilibrium torus for $\beta=0.1$ is shown.

Although we employ static mesh refinement near the equatorial region, a disk scale height cannot be resolved well if the gas temperature drops below $\sim 0.01~{}T_{\infty}$, where $T_{\infty}$ is the temperature of the ambient gas. In order to avoid this issue, we set a minimum temperature below which radiative cooling turns off, i.e., $\Lambda(T\leq T_{\rm min})=0$. We adopt $T_{\rm min}=0.1~{}T_{\infty}$ as our fiducial value and study two cases with $T_{\rm min}=0.2~{}T_{\infty}$ and $0.05~{}T_{\infty}$. As discussed in §4.2, the Toomre parameter depends on the choice of $T_{\rm min}$ somewhat sensitively. We will conduct an extrapolation of our results for lower values of $T_{\rm min}\sim 10^{4}~{}{\rm K}~{}(\simeq 10^{-3}~{}T_{\infty})$ and construct the global structure of the accretion disk.

We adopt the outflow boundary condition at the innermost/outermost grid zones from Stone & Norman (1992), where zero gradients cross the boundary are imposed on physical quantities in order to avoid spurious reflection of wave energy at the boundary. In addition, we impose $v_{r}\leq 0$ at the inner boundary (i.e., inflowing gas from ghost cells is prohibited). At the poles, the reflective condition is imposed on the circumferential component of the velocity $v_{\theta}$. After the transition to a cold disk, we assume an isothermal Keplerian disk in the inner ghost grids at $\theta=\theta_{\rm cold}$, where $T(r_{\rm in},\theta_{\rm cold})\leq 1.2~{}T_{\rm min}$111 When the gas density is not high enough to permit cooling, the temperature is higher than the initial value of $T_{\infty}\gg T_{\rm min}$.; in practice, the rotational velocity is set to $v_{\phi}(r<r_{\rm in},\theta_{\rm cold})=v_{\rm Kep}~{}[1-(5/2)(H/r)^{2}]^{1/2}$.

2.3 Basic quantities and parameter choices

Before we discuss the simulation results, we introduce a dimensionless physical quantity which characterizes BH accretion systems. As a reference of the accretion rate, we define the Bondi accretion rate for adiabatic and fully ionized gas (i.e., $\gamma=5/3$ and $\mu=0.62$)222Note that both of $\gamma$ and $\mu$ are solved self-consistently in our simulations, but fixed to evaluate the Bondi rate.,

[TABLE]

where $q(\gamma)=1/4$ for $\gamma=5/3$, $\rho_{-22}=\rho_{\infty}/(10^{-22}~{}{\rm g~{}{\rm cm}^{-3}})$, $M_{8}=M_{\bullet}/(10^{8}~{}{\rm M}_{\odot})$ and $T_{7}=T_{\infty}/(10^{7}~{}{\rm K})$. The accretion rate normalized by the Eddington rate is given by

[TABLE]

where $\dot{M}_{\rm Edd}(\equiv 10~{}L_{\rm Edd}/c^{2})=2.3~{}M_{8}~{}{\rm M}_{\odot}~{}{\rm yr}^{-1}$ is the Eddington accretion rate with a $10\%$ radiative efficiency.

In this work, we set the BH mass to $M_{\bullet}=10^{8}~{}{\rm M}_{\odot}$, the temperature of the ambient gas to $T_{\infty}=10^{7}~{}{\rm K}$ ($c_{\infty}=470~{}{\rm km~{}s}^{-1}$) at the initial state and $\beta=0.1$ and $\alpha_{0}=0.01$ as our fiducial case. We study the effects of radiative cooling and heating, varying the gas density $\rho_{\infty}$ and thus $\dot{m}_{\rm B}$.

3 Results

Fig. 4 shows the time evolution of gas accretion rates through the inner-most grid ($r=r_{\rm min}$) for four different values of $\rho_{\infty}=10^{-24}$ (red), $10^{-23}$ (green), $3\times 10^{-23}$ (magenta) and $10^{-22}~{}{\rm g~{}cm}^{-3}$ (blue), each of which corresponds to $\dot{m}_{\rm B}\simeq 4\times 10^{-5}$, $4\times 10^{-4}$, $1.2\times 10^{-3}$ and $4\times 10^{-3}$. The accretion rates are shown in units of ${\rm M}_{\odot}~{}{\rm yr}^{-1}$ (top) and are normalized by the Bondi accretion rate with $\gamma=5/3$ for each case (bottom). The simulation time is normalized by the dynamical timescale at the Bondi radius, $t_{\rm dyn}(\equiv R_{\rm B}/c_{\infty})\simeq 4.1\times 10^{3}~{}{\rm yr}~{}M_{8}T_{7}^{-3/2}$. For the lowest density ($\rho_{-22}=0.01$ and $\dot{m}_{\rm B}\simeq 4\times 10^{-5}$), the accretion rate increases and saturates at $t\simeq 2~{}t_{\rm dyn}$. The saturated value is as small as $\simeq 0.02~{}\dot{M}_{\rm B}$, which means that the BH accretion rate is significantly reduced from the Bondi rate. As long as $\dot{m}_{\rm B}\la 10^{-3}$, where the accreting gas is adiabatic, the accretion rate is simply proportional to the ambient gas density. Therefore, the normalized accretion rate behaves in a self-similar way.

On the other hand, for the highest density ($\rho_{-22}=1$ and $\dot{m}_{\rm B}=4\times 10^{-3}$), the self-similar behavior of the accretion rates is no longer valid. The gas accretion rate further increases even after a few dynamical timescales and ultimately reaches the Bondi rate for isothermal gas ($\gamma=1$), which is $e^{3/2}(\simeq 4.5)$ times higher than the Bondi rate for $\gamma=5/3$ (note that for the case with $\dot{m}_{\rm B}=1.2\times 10^{-3}$, a steady state is realized at $t\ga 230~{}t_{\rm dyn}$). This result clearly shows that radiative cooling boosts the normalized accretion rate by three orders of magnitude. As a result, the radiative luminosity from the nuclear regions is also enhanced by many orders of magnitude, $L/L_{\rm Edd}\simeq 6.3\times 10^{-3}$ (see discussion in §4.3).

Fig. 5 presents the radial structure of the angle-integrated mass inflow (dashed) and net accretion (solid) rates for three different values of $\rho_{\infty}$ ($4\times 10^{-5}\leq\dot{m}_{\rm B}\leq 4\times 10^{-3}$). The inflow and outflow rate are calculated as

[TABLE]

and the net accretion rate is defined as $\dot{M}_{\rm net}\equiv-\dot{M}_{\rm in}-\dot{M}_{\rm out}$. Note that $\dot{M}_{\rm in}<0$ and $\dot{M}_{\rm out}>0$. For the two cases with lower densities ($\dot{m}_{\rm B}\simeq 4\times 10^{-5}$ and $4\times 10^{-4}$), the inflow rates decrease towards the center, following $\dot{M}_{\rm in}\propto r^{p}$ (red and green dashed curves). The power-law index is consistent with that in CDAFs ($p\simeq 1$). This fact indicates that convective motion driven by viscous heating suppresses the gas inflow (e.g., Narayan et al. 2000; Quataert & Gruzinov 2000, see also Inayoshi et al. 2018 for details). The net accretion rates are almost constant at $\sim 0.02~{}\dot{M}_{\rm B}$ within the Bondi radius. For the highest-density case ($\dot{m}_{\rm B}\simeq 4\times 10^{-3}$), the inflow rate becomes comparable to the net accretion rate within the centrifugal radius. As shown in Fig 4, the net accretion rate, normalized by the Bondi rate for $\gamma=5/3$, is two orders of magnitude higher than those in the lower-density cases. The radial dependence of the inflow rate becomes flatter ($p\sim 0.3$) because convective motion ceases due to efficient radiative cooling.

In Fig. 6, we present the distribution of the gas density (left) and temperature (right) at two different elapsed times for the run with $\dot{m}_{\rm B}=4\times 10^{-3}$. The velocity vectors are over plotted in the temperature distribution. At the early stage of $t\la 10~{}t_{\rm dyn}$ (left panel), the accretion flow is highly turbulent (but subsonic) and forms a geometrically-thick torus inside the Bondi radius. Since the accreting gas is still adiabatic, the structure of the accretion flow is similar to that of a CDAF. We note that 3D simulations of non-axisymmetric accretion flows performed by Igumenshchev et al. (2000) suggest that the results are essentially similar to those obtained in 2D simulations, and the turbulent eddies are nearly axisymmetric and transport angular momentum inwards. This allows us to justify that our 2D simulations can capture the essential physics. At the late stage of $t\ga 10~{}t_{\rm dyn}$ (right panel), the accreting gas begins to cool, collapse towards the equator and forms a geometrically-thin disk. The disk increases its mass via gas supply from the Bondi radius, and feeds the central BH with the aid of viscous processes transporting the angular momentum outwards. Since radiative cooling is inefficient above the disk, the cold disk is embedded in hot gas ($T\gg T_{\infty}$). As shown in Fig 4, the accretion flow settles into a steady state at $t\ga 100~{}t_{\rm dyn}$, where the BH accretion rate is as high as $\simeq 0.02~{}\dot{M}_{\rm Edd}$. In the steady state, the radiative luminosity from the BH ($r<r_{\rm in}$) is on the order of $\sim 0.01~{}L_{\rm Edd}$. The combination of heating and outward momentum caused by radiation produces weak outflows at $r\la R_{\rm B}$. A small fraction of the gas can escape from the Bondi radius with a velocity of $v_{r}\simeq 2~{}c_{\infty}$, but the mass outflow does not reduce the BH accretion rate significantly.

In Fig. 7, we show the density distribution and velocity vectors at three elapsed times around $t\simeq 120~{}t_{\rm dyn}$ over a short duration with $\Delta t\la 4~{}t_{\rm dyn}$ for the model with $\dot{m}_{\rm B}\simeq 4\times 10^{-3}$. Since the outflow velocity is comparable to the sound speed at the Bondi radius, a large fraction of the outflowing gas to the polar regions cannot escape outside the Bondi radius. Instead, the gas is compressed by the ambient gas and forms clumps around $\sim R_{\rm B}$. The clumps are captured by the gravity of the BH and accrete back onto the disk again. Such a burst of gas accretion increases the energy output from the disk (viscous dissipation) and BH (radiation). As a result, mass outflows are launched to the bipolar directions episodically.

Fig. 8 shows the radial profiles of the gas density (top) and temperature (bottom) at the mid-plane for $\dot{m}_{\rm B}\simeq 4\times 10^{-3}$ at different elapsed times, which correspond to filled circles shown in Fig. 4. At the early stage (red curves, epoch 1), the gas density and temperature follow $\rho\propto r^{-1/2}$ and $T\propto r^{-1}$ within the centrifugal radius, respectively. The radial structures are consistent with those of a CDAF solution (e.g., Quataert & Gruzinov, 2000). At the late stages, radiative cooling induces gas collapse onto the mid-plane. This transition begins to occur around the centrifugal radius (see green curves, epoch 2), where the cooling condition is satisfied first in the flow ($t_{\rm cool}<t_{\rm th,vis}$; see §4.1). After the transition, the gas density in the disk dramatically increases to $\simeq 10^{8}~{}\rho_{\infty}$ at $r\simeq r_{\rm in}$, and the disk outer edge moves outwards by angular momentum transport (see epoch $3-6$). The gas temperature drops to the threshold value we set ($T_{\rm min}=0.1~{}T_{\infty}$) via optically-thin cooling. Note that this transition occurs in an unstable way because the cooling timescale becomes shorter with decreasing temperature as $t_{\rm cool}\propto T^{s}$, where $s\simeq 0.5-1.5$ at $T\la 10^{7}~{}{\rm K}$ due to free-free emission, metal lines and recombination. In the steady state, the density profile at the mid-plane follows

[TABLE]

where $\dot{M}_{\rm d}$ is the accretion rate through the disk, comparable to the Bondi rate for isothermal gas, and $c_{\rm s}$ is the sound speed of the gas in the disk with a temperature of $T\simeq T_{\rm min}$. The solid line in the top panel shows Eq. (21) at $r\leq R_{\rm C}(=0.1~{}R_{\rm B})$. The temperature exterior to the cold disk, on the other hand, becomes hotter due to Compton heating ($T_{\rm C}=10^{8}~{}{\rm K}$) as the BH feeding rate increases.

In Fig. 9, we present the angular profiles of the density (top) and temperature (bottom) as a function of the polar angle $\theta$, at radial positions of $r_{\rm in}\leq r\leq R_{\rm B}$. Inside the Bondi radius, the density increases toward the center and the mid-plane ($\theta=\pi/2$). Within the cold disk ($0.4\la\theta/\pi\la 0.6$), the temperature is constant at $T=T_{\rm min}(=0.1~{}T_{\infty})$. Above the disk, the gas is heated up with decreasing radii. Note that the cold disk region is located within the region where grids in the polar direction are refined ($\theta_{0}\leq\theta\leq\pi-\theta_{0}$).

4 Physical interpretation

4.1 The formation of a cold accretion disk

For lower densities at the Bondi radius ($\dot{M}_{\rm B}/\dot{M}_{\rm Edd}\ll 10^{-3}$), the accreting gas onto a BH forms a geometrically-thick torus. Namely, the gas density, temperature and inflow rate are consistent with those of a CDAF solution: $\rho\propto r^{-1/2}$, $T\propto r^{-1}$ and $\dot{M}\propto r$ (e.g., Quataert & Gruzinov, 2000). Inside the torus structure, viscous energy dissipation heats the gas and drives convective motion on a thermal timescale of $t_{\rm th,vis}\simeq 1/[\gamma(\gamma-1)\alpha\Omega]^{-1}\propto r^{3/2}$. With the CDAF solution, the bremsstrahlung cooling rate per volume is $Q^{-}_{\rm br}\propto\rho^{2}T^{1/2}\propto r^{-3/2}$ and the cooling timescale is $t_{\rm cool}\propto\rho T/Q^{-}_{\rm br}\propto r^{0}$. Fig. 10 shows the ratio of $t_{\rm th,vis}/t_{\rm cool}$ as a function of $\dot{m}_{\rm B}$ for different values of $\beta(=R_{\rm C}/R_{\rm B})$ and $T_{\infty}$, adopting the radial density and temperature profile of a CDAF solution inside the Bondi radius (Inayoshi et al., 2018). The ratio of the two timescales is estimated as

[TABLE]

where the ratio is evaluated at the outer edge of the torus ($r\simeq 2R_{\rm C}$). We note that this scaling relation agrees with the results shown in Fig. 10, as long as the angular momentum is sufficiently small ($\beta\la 0.3$). From Eq. (22), we obtain the critical accretion rate above which bremsstrahlung cooling plays an important role in the accretion flow as

[TABLE]

In Fig. 10, we also present the results of numerical simulations where a cold accretion disk forms (filled squares) and those where the accretion flow is adiabatic (open squares). For a wide range of the angular momentum ($0.01\leq\beta\leq 0.3$) and ambient gas temperature ($5\times 10^{6}\leq T_{\infty}/{\rm K}\leq 10^{7}$), the transition from adiabatic accretion flows to cold disk accretion occurs at $4\times 10^{-4}\la\dot{m}_{\rm cool}\la 2\times 10^{-3}$. We note that Eq. (23) is not valid for lower temperatures where metal-line cooling dominates bremsstrahlung cooling significantly ($T_{\infty}\ll 5\times 10^{6}~{}{\rm K}$).

It is worth comparing our results to those obtained in previous work which investigated the radiative cooling effects in hot and optically-thin accretion flows. Early studies by Yuan (2001, 2003) have discussed the properties of luminous hot accretion flows (LHAFs) in the inner region ($r\la 10^{2}~{}R_{\rm Sch}$), where the sum of viscous heating and compressional heating balance with cooling via Coulomb collisions between ions and electrons. For higher rates of $\dot{M}\ga\dot{M}_{\rm LHAF}\simeq 1.5~{}\alpha^{0.7}\dot{M}_{\rm Edd}$333 The critical rate would depend on the choice of the outer edge of the disk. For $R_{\rm out}=10^{4}~{}R_{\rm Sch}$, Yuan (2001) obtained $\dot{M}_{\rm LHAF}/\dot{M}_{\rm Edd}\simeq 0.25$, $5\times 10^{-2}$ and $5\times 10^{-4}$ for $\alpha=0.1$, $10^{-2}$ and $10^{-3}$, respectively., hot accretion flows no longer exist because radiative cooling efficiently carries the energy away. We note that the critical value of $\dot{M}_{\rm LHAF}$ is two orders of magnitude higher than the critical rate given by Eq. (23) assuming $\alpha=0.01$. Another qualitative difference is that accreting gas from the Bondi radius begins to cool at the disk outer radius, while in LHAFs the inner region of the disk at $<10^{2}~{}R_{\rm Sch}$ collapses via cooling. Wu et al. (2016) have also studied the cooling transition of LHAFs at $r\la 10^{2}~{}R_{\rm Sch}$, performing hydrodynamical simulations. They found a critical accretion rate of $\dot{M}_{\rm LHAF}\sim 3\alpha\dot{M}_{\rm Edd}$, above which cold and dense clumpy and/or filamentary structures are formed within hot gas. We note that in their simulations, gas supply from larger scales is not considered unlike our simulations. The different treatment of the outer boundary conditions would affect the cooling conditions and significantly lower the critical accretion rate given in Eq. (23).

In this paper, we do not discuss the regime of higher accretion rate where various feedback processes would be important, e.g., radiation, winds and jets (e.g., Sazonov et al. 2005; Ciotti & Ostriker 2001; Proga 2007; Yuan et al. 2009; Ciotti et al. 2009; Booth & Schaye 2009; Milosavljević et al. 2009; Yuan & Li 2011; Novak et al. 2011; Choi et al. 2012, 2014; Costa et al. 2014; Xie et al. 2017; Yuan et al. 2018; Yoon et al. 2018; Bu & Yang 2019, and see Fabian 2012 references therein). We define the critical BH accretion rate $\dot{m}_{\rm fb}$, above which a steady state no longer exists due to feedback but gas accretion occurs episodically. The exact value of $\dot{m}_{\rm fb}$ is still uncertain, depending on which feedback processes dominate. Instead of discussing each process, we here adopt a critical radiation luminosity of $L_{\rm bol}/L_{\rm Edd}\simeq 0.02$ as a reference (e.g., McClintock & Remillard, 2006). Since the BH accretion rate through a cold (isothermal) disk is $\sim e^{3/2}(\simeq 4.5)$ times higher than the Bondi accretion rate estimated with $\gamma=5/3$, the critical value of the Bondi rate is $\dot{m}_{\rm fb,B}\equiv\dot{m}_{\rm fb}/e^{3/2}\simeq 6.7\times 10^{-3}$, where $L_{\rm bol}(\dot{m}_{\rm fb})=0.02~{}L_{\rm Edd}$ (see Fig. 12 and discussion in §4.3). We leave the study of the strong feedback domain for future work.

4.2 Toomre instability and an extension of the disk solution towards the BH

We now focus on the intermediate regime of $\dot{m}_{\rm cool}\la\dot{m}_{\rm B}\la\dot{m}_{\rm fb,B}$, where a dense and geometrically-thin disk forms due to radiative cooling, but feedback does not affect the BH accretion. The cold disk may then fragment into clumps and form stars by a spiral-mode gravitational instability, which is characterized by the Toomre parameter (Toomre, 1964):

[TABLE]

Assuming a steady state, the gas surface density is estimated as $\Sigma\simeq\dot{M}_{\rm d}/(3\pi\nu)$, where $\dot{M}_{\rm d}$ is the gas accretion rate through the disk. Thus, we obtain

[TABLE]

where $T_{\rm min,4}=T_{\rm min}/(10^{4}~{}{\rm K})$. Fig. 11 shows the radial profile of the Toomre parameter for $T_{\rm min}=0.1~{}T_{\infty}~{}(=10^{6}~{}{\rm K})$ (red solid curve). In the cold disk region, the Toomre’parameter is almost constant as a function of radius and has a minimum value of $Q\simeq 10^{3}$ at the centrifugal radius ($R_{\rm C}=0.1~{}R_{\rm B}$). This result agrees with the analytical estimate shown by the star symbol (see Eq. 25) within a factor of two. Dashed curves present the cases for $T_{\rm min}=10^{6}~{}{\rm K}$ (red), $2\times 10^{6}~{}{\rm K}$ (green) and $5\times 10^{5}~{}{\rm K}$ (blue), but with an inner-most radius of $r_{\rm in}=3.4\times 10^{-2}~{}R_{\rm B}$, which is twice larger than our fiducial case (solid). The Toomre values at $r=R_{\rm C}$ decreases for lower $T_{\rm min}$, following $Q\propto T_{\rm min}^{3/2}$. Therefore, we would expect that once radiative cooling leads to the formation of a cold disk and the temperature decreases below $10^{4}~{}{\rm K}$, the disk would become unstable because of its self-gravity.

When the Toomre parameter decreases to $\sim O(1)$, the accretion disk becomes gravitationally unstable (Toomre, 1964) and induces non-axisymmetric structure like spiral arms transporting angular momentum outward efficiently (Gammie, 2001; Goodman, 2003; Rice et al., 2005; Vorobyov & Basu, 2006; Kratter et al., 2008; Machida et al., 2010; Kuiper et al., 2011; Zhu et al., 2012; Takahashi et al., 2016). The strength of gravitational torque is characterized by an effective viscosity of $\alpha_{\rm eff}\sim O(1)$. To mimic this effect, we parameterize the effective viscosity as a function of $Q$ in the following:

[TABLE]

where $\alpha_{0}\simeq 0.01$ is the viscous parameter induced by the MRI and $\alpha_{\rm max}\sim 0.1-1$ is the maximum viscous strength. Hydrodynamical simulations of a self-gravitating disk suggest that the disk tends to fragment when the strength of gravitational torque exceeds a critical threshold value of $\alpha_{\rm max}$. The threshold value depends on mass accretion from an infalling envelope, realistic radiative cooling, and radiative trapping of energy inside clumps (Zhu et al., 2012). In the context of massive star formation, the critical value for fragmentation is suggested to be $\alpha_{\rm max}\sim 0.1-1$ (e.g., Kratter et al., 2008)444Numerical simulations for isolated disks without infalling flow from larger scales suggest lower critical values of $\alpha_{\rm max}\simeq 0.06$ (Rice et al., 2005).. Based on those results, we adopt $\alpha_{\rm max}=1$ as a reference value. We also note that the choice of the two free parameters of $(A,B)$ hardly affects the strength of gravitational torque (see Appendix A) and thus they are set to $(A,B)=(10,1)$. Using Eqs. (25) and (26), we obtain the critical accretion rate above which the disk becomes gravitationally unstable ($Q\sim 1$ and $\alpha_{\rm eff}\ga\alpha_{0}$) as

[TABLE]

or

[TABLE]

where $\alpha_{\rm eff}\geq 0.03$ is adopted as the threshold value. The minimum temperature is as high as $T_{\rm min}\simeq 3000~{}{\rm K}$, which is determined by H- continuum radiation. Thus, once radiative cooling leads to the formation of a dense disk in the Bondi radius of a BH with $M_{\bullet}\ga 10^{8}~{}{\rm M}_{\odot}$, it would turn gravitationally unstable immediately. For less massive BHs with $M_{\bullet}\sim 10^{7}~{}{\rm M}_{\odot}$, accretion rates as high as $\dot{m}_{\rm fb,B}\sim 10^{-2}$ are required to trigger the onset of gravitational instability.

As a specific example, we briefly discuss the disk structure of the nuclear BH in M87 galaxy ($M_{\bullet}=6.15\times 10^{9}~{}{\rm M}_{\odot}$ and $\dot{m}_{\rm B}=1.34\times 10^{-3}$; Gebhardt et al. 2011; Russell et al. 2015), which satisfies Eq. (27). Therefore, the nuclear accretion disk within the Bondi radius becomes unstable, leading to efficient fragmentation and star formation (e.g., Goodman, 2003; Tan & Blackman, 2005)555 The nuclear disk in M87 is no longer stable because $\psi<10^{-2}$, where $Q\la 1$ and $\alpha_{\rm eff}\sim\alpha_{\rm max}$ (see Appendix A).. Since star formation will consume the disk mass and the newly formed stars also will heat the disk, the remaining gas can be stabilized, result in $Q\sim 1$. In the marginally unstable disk, the surface density and disk mass are given by

[TABLE]

and

[TABLE]

where $\beta(=R_{\rm C}/R_{\rm B})=0.1$ is adopted. A fraction $\epsilon_{\star}$ of the disk mass would be consumed by star formation. Assuming the star formation efficiency $\epsilon_{\star}=0.01$ inferred from the Kennicut-Schmidt law (Schmidt, 1959; Kennicutt, 1998; Krumholz & McKee, 2005), the total stellar mass is $M_{\star}\sim 10^{6}~{}{\rm M}_{\odot}$ and the expected number of supernovae (SNe) is $N_{\rm SN}\simeq M_{\star}/m_{\star}\sim 10^{4}$, where $m_{\star}\sim 100~{}{\rm M}_{\odot}$ is the stellar mass required to generate an SN for a Kroupa initial mass function (Kroupa, 2001). Since the total SN energy $N_{SN}E_{\rm SN}\simeq 10^{55}E_{51}$ erg is injected over the lifetime of massive stars $t_{\rm life}\sim 20$ Myr, thus the kinetic luminosity is $L_{\rm kin}\simeq 2\times 10^{40}~{}{\rm erg~{}s}^{-1}$. Let us assume that super-bubbles expand as pressure-driven snowplows with no radiative cooling in their interior. Then, the criterion for “breakout” from the cold disk is expressed as $L_{\rm kin}\gg\rho_{\rm a}v_{\rm t}^{3}A$, where $\rho_{\rm a}$ is the ambient gas density, $v_{\rm t}$ is the turbulent velocity in the disk, and $A~{}(\sim H^{2})$ is the total surface area of the bubbles (Mac Low & McCray, 1988; Koo & McKee, 1992; Kim et al., 2017). In this case for M87 ($\rho_{\rm a}\sim 10^{-16}~{}{\rm g~{}cm}^{-3}$, $v_{\rm t}\sim c_{\rm s}\sim 10~{}{\rm km~{}s}^{-1}$ and $H\sim 2\times 10^{17}$ cm at $r\simeq R_{\rm C}$), the ratio of the SN kinetic luminosity to that of turbulent gas is estimated as $L_{\rm kin}/(\rho_{\rm a}v_{\rm t}^{3}A)\sim 5000\gg 1$. This implies that SN feedback would likely blow the gas away from the cold disk, making star formation episodes brief.

In a giant elliptical galaxy harboring an SMBH with a mass of $M_{\bullet}\ga 10^{9}~{}{\rm M}_{\odot}$ fed at $\dot{m}_{\rm B}\ga\dot{m}_{\rm cool}$, star formation activities in the nuclear region within the Bondi radius would be ubiquitous (e.g., Tan & Blackman, 2005; Gan et al., 2018). In fact, this hypothesis would explain the following observational features of the M87 system: the existence of a cold disk in the nuclear region and a low level of star formation. The cold disk, which is traced by H$\alpha$ emission, has a size of $\sim 20-40~{}{\rm pc}$ (e.g., Harms et al., 1994; Walsh et al., 2013), suggesting that the formation of the cold disk is triggered by radiative cooling inside the Bondi radius ($\sim 120~{}{\rm pc}$). The star formation rate for M87 is estimated as ${\rm SFR}\la 5\times 10^{-2}~{}{\rm M}_{\odot}~{}{\rm yr}^{-1}$ based on the $70\mu m$ infrared luminosity (Temi et al., 2009; Calzetti et al., 2010). The SFR would be considered as the upper limit because the observed aperture is considerably larger than the Bondi scales, and there could be a level of contamination from synchrotron emission associated with the jet for M87. This rate is consistent with our order-of-magnitude estimate of ${\rm SFR}\la M_{\star}/t_{\rm life}\simeq 5\times 10^{-2}~{}{\rm M}_{\odot}~{}{\rm yr}^{-1}$, implying that a large fraction of the gas supplied from the Bondi radius is consumed by star formation, namely, ${\rm SFR}/\dot{M}_{\rm B}\la 0.26$. Presumably, a significant fraction of other giant ellipticals had experienced short episodes of star formation even after the galaxies have terminated major episodes of star formation at $z\la 1$ (e.g., Thomas et al., 2005; Pérez-González et al., 2008). At lower redshifts, the number fraction of elliptical galaxies having ongoing nuclear star bursts, the so-called blue ellipticals, is as small as $\simeq 3~{}\%$ (Tojeiro et al., 2013). In addition, most nuclear disks in the centers of nearby ellipticals seem stable against star formation, i.e., $Q>1$ (Boizelle et al., 2017). Thus, the quiescent phases are expected to last the order of billions of years.

4.3 Radiative luminosity vs. Bondi accretion rate

We here estimate the radiation luminosity ($L_{\rm bol}$) produced from accretion flows and BH feeding rate ($\dot{M}_{\bullet}$) for a given Bondi accretion rate at a range of $10^{-6}\la\dot{M}_{\rm B}/\dot{M}_{\rm Edd}\la 10^{-2}$. The results are summarized in Fig. 12.

Below the cooling threshold ($\dot{m}_{\rm B}<\dot{m}_{\rm cool}\sim 10^{-3}$), the accretion flows are adiabatic and the gas inflow rate decreases towards the center due to convection as shown in Fig. 5. In Paper I (Inayoshi et al., 2018), we estimated the transition radius within which the inflow rate becomes constant because energy transport by thermal conduction would dominate that by convection. The inferred accretion rate onto the central BH is

[TABLE]

where $f_{\rm c}$ is the conductivity suppression factor because thermal conduction in the directions perpendicular to magnetic fields can be suppressed. The value of the suppression factor has been discussed by various theoretical arguments and estimated as $f_{\rm c}\sim 0.1$ (e.g., Narayan & Medvedev, 2001; Maron et al., 2004). For such low BH feeding rates, the radiative efficiency decreases with the accretion rate at $10^{-9}\la\dot{M}_{\bullet}/\dot{M}_{\rm Edd}\la 10^{-5}$ and has been estimated with radiation magneto-hydrodynamical simulations (Ryan et al., 2017, $L\propto\dot{M}_{\bullet}^{1.7}$, see also orange squares in Fig. 2). Therefore, adopting the radiation efficiency, the radiation luminosity model in this domain can be estimated as

[TABLE]

which is shown by the red shaded region in Fig. 12. The width of the region reflects the uncertainties of the conductivity suppression factor $0.03\la f_{\rm c}\la 0.3$. We note that if the BH feeding rate is equal to the Bondi accretion rate ($\dot{M}_{\bullet}=\dot{M}_{\rm B}$; black dashed), the radiation luminosities are much higher than that given by Eq. (32).

For higher Bondi rates ($\dot{m}_{\rm B}\ga\dot{m}_{\rm cool}$), the accretion flow cools and forms a dense accretion disk where the accretion rate is almost constant as a function of radius (see Fig. 5). Because the equation of state of the gas changes from adiabatic to isothermal one, the actual inflow rate becomes comparable to the Bondi rate for isothermal gas and thus we obtain $\dot{M}_{\bullet}\simeq e^{3/2}\dot{M}_{\rm B}$. Therefore, the radiation luminosity after the transition is estimated as

[TABLE]

which would be valid for $\dot{m}_{\rm cool}\la\dot{m}_{\rm B}\la\dot{m}_{\rm fb,B}$, i.e., unless $L_{\rm bol}\ga 0.02~{}L_{\rm Edd}$ where radiative feedback suppresses gas inflows from larger radii. This level of luminosities ($L_{\rm bol}/L_{\rm Edd}\simeq 10^{-3}-10^{-2}$) is consistent with the observational results for nearby Seyfert galaxies (Koss et al., 2017).

For comparison, we present the observational results of Sgr A⋆ and BHs in M31 and M87 (blue asterisks in Fig. 12; see also Inayoshi et al. 2018). Those nearby SMBHs are known as quiescently accreting ones at rates of $\dot{M}_{\rm B}/\dot{M}_{\rm Edd}\sim 10^{-5}-10^{-3}$. The detailed gas properties in the nuclei are studied nicely by Chandra X-ray observations resolving the angular sizes of their Bondi radii (e.g., Baganoff et al., 2003; Garcia et al., 2010; Russell et al., 2015, and references therein). Our theoretical estimate agrees well with the observational results for Sgr A⋆ and M31, and further explains the observed luminosity of M87, which is several times higher than the adiabatic ones (red dashed) because of radiative cooling.

4.4 Gas accretion from the circumnuclear disk to the Bondi scale

So far, we have focused on the accretion physics around and inside the Bondi radius of the central BH, where a small centrifugal radius is set so that $R_{\rm C}<R_{\rm B}$. We here briefly discuss gas supply to the Bondi scale from a circumnuclear disk with a size of $\sim 10-100$ pc (CND; Thompson et al., 2005; Ballantyne, 2008)666This implicitly assumes a larger centrifugal radius. However, the gas inflow velocity could be a significant fraction of the Keplerian value in a CND because of the Toomre instability and/or turbulence associated with star formation activities (Wada & Norman, 2002).. Since CNDs would be massive reservoirs of molecular gas outside but very close to the Bondi radius (e.g., Hicks et al. 2009 and Izumi et al. 2016 for warm and cold molecular gas, respectively), we naturally expect active star formation and/or gas feeding to the nuclear SMBH. Interestingly, observations imply that the BH accretion rate would correlate with the SFR in the CNDs of Seyfert galaxies (e.g., Imanishi et al., 2011; Diamond-Stanic & Rieke, 2012; Esquej et al., 2014; Izumi et al., 2016). This could suggest that star formation activities in CNDs (outside the Bondi scales) control the amount of gas fueling (e.g., Shlosman et al., 1989; Goodman, 2003; Thompson et al., 2005; Tan & Blackman, 2005; Vollmer et al., 2008; Kawakatu & Wada, 2008; Chamani et al., 2017; Gan et al., 2018).

Numerical simulations studying the structure of a CND (Wada & Norman, 2002; Wada et al., 2009) have found that SNe associated with active star formation lead to turbulence in a CND and thus the turbulent viscosity efficiently transports angular momentum outward and mass inward. The turbulent velocity associated with energy injection by SNe is estimated as777 This equilibrium state would be realized unless super-bubbles created by SNe break out the disk scale height (cf. discussion in §4.2). When the break-out timescale is longer than the sound crossing timescale in the disk, the disk would be highly turbulent with typical velocity of $v_{\rm t}$.

[TABLE]

where $\eta$ is the efficiency per mass which denotes the fraction of the energy from SNe and $\tau_{\star}$ is the star formation timescale. Presumably, the strength of an effective viscosity is as high as $\alpha_{\rm eff}\sim O(1)$ and thus the gas accretion rate through the disk is estimated as

[TABLE]

Assuming that the surface density follows a power law in the CND, $\Sigma=\Sigma_{0}(r/r_{\rm out})^{-\gamma}$, the CND mass is $M_{\rm d}=2\pi\Sigma_{0}r_{\rm out}^{2}[1-(r_{\rm in}/r_{\rm out})^{2-\gamma}]/(2-\gamma)$ and thus we obtain

[TABLE]

where $r_{\rm in(out)}$ is the inner (outer) radius of the CND and $x\equiv r_{\rm in}/r_{\rm out}$. The gas accretion rate at the inner radius $r_{\rm in}$ is given by

[TABLE]

where ${\rm SFR}\equiv M_{\rm d}/\tau_{\star}$, $\eta=10^{-3}~{}{\rm M}_{\odot}^{-1}$, $E_{\rm SN}=10^{51}~{}{\rm erg}$ and $\gamma=3/2$. Adopting $r_{\rm in}\simeq R_{\rm B}$, the gas accretion rate through the Bondi radius is given by

[TABLE]

Therefore, star formation activities in the nuclear region in a CND at a rate of ${\rm SFR}\simeq 0.1~{}{\rm M}_{\odot}~{}{\rm yr}^{-1}$ provide an inflow rate at $\sim 10^{-2}~{}{\rm M}_{\odot}~{}{\rm yr}^{-1}$ within the Bondi radius of an SMBH. For the typical mass of BHs hosted in Seyfert galaxies/AGN ($M_{\bullet}\sim 10^{7}-10^{8}~{}{\rm M}_{\odot}$), the Bondi accretion rate is high enough to form a cold, dense accretion disk, leading to a high BH feeding rate of $\dot{M}_{\bullet}\simeq\dot{M}_{\rm B}$ and a high luminosity of $L_{\rm bol}/L_{\rm Edd}\simeq 10^{-3}-10^{-2}$ (see also Koss et al., 2017). The result is consistent with the fact that the BH feeding rate inferred from the nuclear luminosity tightly correlates with the SFR in the CND ($<100~{}{\rm pc}$), i.e., $\dot{M}_{\bullet}/{\rm SFR}\sim O(0.1)$ (e.g., Diamond-Stanic & Rieke, 2012; Esquej et al., 2014; Izumi et al., 2016). On the other hand, for higher BH masses ($M_{\bullet}\gg 10^{8}~{}{\rm M}_{\odot}$), the Bondi rate would not be high enough to trigger the cooling transition.

4.5 Caveats

4.5.1 magnetic fields

In this work, we do not take into account MHD effects explicitly, but adopt the $\alpha$-viscosity prescription to treat angular momentum transport, mimicking some properties of MRI turbulence. The properties of radiatively inefficient MHD accretion flows are affected qualitatively and quantitatively by the initial configuration and boundary conditions for magnetic fields (e.g., Stone & Pringle, 2001; Igumenshchev et al., 2003; McKinney et al., 2012; Narayan et al., 2012; Yuan et al., 2012). Assuming that a toroidal magnetic field or multiple poloidal loops is initially set or injected with accreting matter from the outer boundary, the accretion flow is dominated by turbulent motion driven by MRI and heating associated with magnetic reconnection. Through the dissipation process, the strength of the magnetic field is saturated at some levels where the ratio of the gas pressure to the magnetic pressure is $\beta\equiv P_{\rm gas}/P_{\rm m}\simeq 10^{2}-10^{3}$. This type of accretion flow is qualitatively similar to the results obtained from hydrodynamical simulations adopting $\alpha$-viscosity (e.g., Inayoshi et al., 2018).

On the other hand, when a magnetic field with poloidal topology is dragged inward by accreting matter and accumulated at the vicinity of the BH, the accretion flow behaves in a very different way. In fact, the flow is less turbulent but is likely to produce outflows and/or jets associated with the amplified strong magnetic field. Yuan et al. (2012) studied the detailed properties of this type of MHD flow and concluded that the flow is more related to the so-called adiabatic inflow-outflow solution (Blandford & Begelman, 2004), where matter can accrete at very low rates due to strong outflows which carry energy and angular momentum away.

It is also worthy investigating this issue for gas accretion from outside the Bondi radius where the geometry and binding energy of the injected mass are different from those studied in previous work. In this case, a more realistic situation is that magnetic fields at larger scales fluctuate with time due to gas turbulent motion, and then are dragged inward by accreting matter. We leave this to future work (but note that Igumenshchev & Narayan (2002) briefly discussed the issue for spherical accretion).

4.5.2 mechanical feedback due to winds

Mechanical feedback due to AGN winds is expected to play an important role on the evolution of galaxies when the BH feeding rate is sufficiently low (e.g., Murray et al., 2005; Ostriker et al., 2010; Gaspari et al., 2012; Choi et al., 2015). Assuming a phenomenological model (Ostriker et al., 2010), one can estimate an upper bound on the momentum input by winds as

[TABLE]

where $\dot{M}_{\rm w}$ is the wind mass-loss rate, $v_{\rm w}$ is the wind velocity, $\dot{M}_{\rm inf}$ is the mass inflow rate at the inner grid of $r=r_{\rm in}$, $\eta$ is the ratio of the wind mass-loss rate to the BH feeding rate ($\dot{M}_{\bullet}\equiv\dot{M}_{\rm inf}-\dot{M}_{\rm w}$), and $\dot{M}_{\bullet}$ is assumed to be 10% of the inflow rate, i.e., $\eta=9$. The outflowing material pushes the accreting gas bound by the BH gravity and could form an momentum-driven expanding shell. The strength of gravitational force exerting the shell with a mass of $M_{\rm sh}$ at a radius of $r$ is estimated as

[TABLE]

where mass accretion to the shell from larger scales is neglected.

Let us compare the strength of momentum input and gravitational force on a mass shell for the two different accretion regimes. For adiabatic cases with $\dot{m}_{\rm B}<10^{-3}$, where the mass inflow rate is reduced from the Bondi rate due to turbulent motion as $\dot{M}_{\rm inf}/\dot{M}_{\rm B}\la 10^{-2}$, the strength of mechanical feedback is estimated as $\dot{P}_{\rm w}\la 3\times 10^{-3}~{}\dot{m}_{\rm B}L_{\rm Edd}/c$. Using the density profile of a CDAF near the mid-plane (Eq. 21 in Inayoshi et al. 2018), one can estimate the gravitational force on the shell as $|F_{\rm grav}|\sim 6\times 10^{-2}~{}\dot{m}_{\rm B}(L_{\rm Edd}/c)(r/R_{\rm C})^{1/2}$ and thus $\dot{P}_{\rm w}/|F_{\rm grav}|\la 0.05$. Next, for a higher rate of $\dot{m}_{\rm B}\sim 10^{-3}$, where the inflowing rate through a cold disk is comparable to the Bondi rate without suppression of accretion, the level of momentum input is boosted up to $\dot{P}_{\rm w}\sim 3\times 10^{-1}~{}\dot{m}_{\rm B}L_{\rm Edd}/c$, but is not strong enough to blow the dense accretion disk away (i.e., $\dot{P}_{\rm w}/|F_{\rm grav}|\ll 1$). Therefore, even in this case where a generous fraction of the inflowing gas is ejected in a high velocity, mechanical feedback does not affect the gas dynamics both in CDAF and cold disk accretion cases.

On the other hand, mechanical feedback due to winds would be likely to blow the gas away near the poles where the gas density is sufficiently low and weak outflows are produced even without explicitly assuming AGN winds (see Fig. 6). As a natural outcome, biconical outflows with opening angles of $\sim 30^{\circ}~{}(70^{\circ})$ would be produced from geometrically thick (thin) disk in lower (higher) accretion-rate regimes. To explore the nature and origin of outflows both for CDAF and cold disk accretion solutions is left for future investigations (but note that Proga et al. 2000 discussed wind production from the nuclear disk at an accretion rate of $\sim\dot{M}_{\rm Edd}$).

5 Summary and conclusion

We study the properties of rotating accretion flows onto SMBHs using axisymmetric two-dimensional radiation hydrodynamical simulations with radiative cooling and BH feedback. The simulations resolve the accretion dynamics of gas outside the BH influence radius through an inner accretion disk. For lower Bondi accretion rates at $\dot{M}_{\rm B}\ll 10^{-3}~{}\dot{M}_{\rm Edd}$, the BH feeding is suppressed by several orders of magnitude from the Bondi rate due to turbulent motion with outflows to the Bondi radius nearly balancing inflows. Thus, the radiative luminosity results in as low as $\sim 10^{-10}-10^{-7}~{}L_{\rm Edd}$. For higher rates of $\dot{M}_{\rm B}\ga 10^{-3}~{}\dot{M}_{\rm Edd}$, the optically-thin accreting gas cools via free-free emission and forms a geometrically-thin disk, which feeds the BH efficiently and increases the radiative luminosity to $\ga 10^{-3}~{}L_{\rm Edd}$.

The transitional behavior of accreting BHs in galactic nuclei from the radiative inefficient phases to cold disk accretion naturally explains (1) the reason for the offset between the observed luminosities and the simplest theoretical predictions for nearby quiescent SMBHs, and (2) the conditions to fuel gas into the nuclear SMBH. In addition, the cold disks formed in galactic nuclei tend to be gravitationally unstable and lead to star formation when the Bondi rate is as high as $\dot{M}_{\rm B}\ga 10^{-2}~{}{\rm M}_{\odot}~{}{\rm yr}^{-1}$. This plausibly explains the correlation between star formation rates and BH feeding rates in Seyfert galaxies.

Acknowledgements

We thank Feng Yuan, Luis Ho, Takuma Izumi, Defu Bu and Zhaoming Gan for useful discussions. This work is partially supported by the National Key R&D Program of China (2016YFA0400702), and the National Science Foundation of China (11721303), the Simons Foundation through the Simons Society of Fellows (Kohei Inayoshi), and supported by Program for Establishing a Consortium for the Development of Human Resources in Science and Technology, Japan Science and Technology Agency and Japan Society for the Promotion of Science (JSPS) KAKENHI (18K13584; Kohei Ichikawa). RK acknowledges financial support via the Emmy Noether Research Group on Accretion Flows and Feedback in Realistic Models of Massive Star Formation funded by the German Research Foundation (DFG) under grant no. KU 2849/3-1 and KU 2849/3-2. Numerical computations were carried out on Cray XC50 at the Center for Computational Astrophysics of the National Astronomical Observatory of Japan.

Appendix A Effective viscosity by the global instability in a gravitationally unstabel disk

After the transition to a geometrically-thin, cold and dense accretion disk, the disk is likely to be unstable to its self-gravity, where the Toomre parameter becomes close to unity (i.e., $Q\la 1$). In order to capture the Toomre instability, we adopt a toy model for the effective viscosity in a gravitationally unstable disk as

[TABLE]

In spite of uncertainties of those parameters, Takahashi et al. (2013) found that for a wide range of $(A,B)$ this model for the effective viscosity can reproduce the results of three-dimensional hydrodynamical simulations of a circumstellar disk under continuous gas supply from larger scales (Machida et al., 2010).

In Fig. 13, we show the effective viscosity (top) and Toomre parameter (bottom) as a function of

[TABLE]

for different combinations of $(A,B)$. For $\psi\ga 1$, the Toomre parameter is significantly above unity and thus $\alpha_{\rm eff}\simeq\alpha_{0}(=0.01)$. With $\psi$ decreasing ($10^{-2}<\psi<1$), the Toomre parameter is almost constant at $Q\simeq 1$ and the effective viscosity increases. At the intermediate region, gravitational torques could transport mass and angular momentum so that the disk structure results in a marginally stable state. For $\psi\ll 10^{-2}$, the effective viscosity is saturated at $\alpha_{\rm eff}=1$ and the Toomre parameter becomes substantially lower than unity. Since the disk no longer keeps stable accretion, therefore the disk would be likely to fragment into clumps. We set the two free parameters to $(A,B)=(10,1)$ in §4.2.