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

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

Kohei Inayoshi, Kohei Ichikawa, Jeremiah P. Ostriker, and Rolf Kuiper

Kavli Institute for Astronomy and Astrophysics, Peking University, Beijing 100871, China,
Frontier Research Institute for Interdisciplinary Sciences, Tohoku University, Sendai 980-8578, Japan
Department of Astronomy, Columbia University, 550 W. 120th Street, New York, NY 10027, USA,
Institute of Astronomy and Astrophysics, University of Tübingen, Auf der Morgenstelle 10, D-72076 Tübingen, Germany
E-mail: (KI)

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 (), the BH feeding is suppressed due to turbulent motion by several orders of magnitude from the Bondi rate with outflows to the Bondi radius nearly balancing inflows. Thus, the radiative luminosity results in as low as , where is the Eddington luminosity. For higher rates of , 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 . 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 . This is a plausible explanation of the correlation observed between star formation rates and BH feeding rates in Seyfert galaxies.

galaxies: nuclei – galaxies: Seyfert – quasars: supermassive black holes

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 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


within which the negative gravitational energy dominates the thermal energy of the gas with a sound speed of , where 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


where 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,


where 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 -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 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.

Figure 1: A schematic overview of BH accretion solutions. The characteristic bolometric luminosities are shown as a function of the gas inflow rates (i.e., the Bondi accretion rates). Both values are normalized by the Eddington luminosity and accretion rate. For lower rates at , the accreting gas behaves adiabatically and forms a geometrically-thick disk. As a result, matter is convectively unstable and the BH accretion rate is well below the Bondi rate (), and the radiative luminosity is as small as (Inayoshi et al., 2018). The pink band encloses the CDAF solutions which explain nearby silent SMBHs. For higher , the gas cools and forms a geometrically-thin disk. The BH accretion rate increases by several orders of magnitude and reaches the Bondi rate (), which causes a sharp rise of the radiative luminosity. The cold accretion disk is likely to be unstable by its self-gravity (), fragment into small clumps and triggers star bursts.

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 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 , where the radiative feedback from the central BH, i.e., gas heating at larger scales around , 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 or even larger galactic scales . At even higher accretion rates of , 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 (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,


and the equation of motion,


where is the density, is the velocity, and is the gas pressure, the gravitational potential is set to , is the distance from the central BH, is the stress tensor due to viscosity, and is the outward net radiation force in the radial direction. The time derivative is the Lagrangian derivative, given by + .

We solve the energy equation of


where is the internal energy per mass. The equation of state of the ideal gas is assumed as , where . 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 ).

The viscous stress tensor is given by


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


where . We calculate the viscous parameter by mimicking some properties of the MRI as


where the strength of viscosity is set to according to MHD simulations of the global disk structure (e.g., Zhu & Stone, 2018; Takasao et al., 2018), and 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 (), 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., , 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


where each term corresponds to the rate associated with bremsstrahlung cooling (), Compton heating/cooling (), and the sum of photoionization heating, line and recombination continuum cooling (). The Compton heating/cooling rate is given by


where is the AGN bolometric luminosity, is the Compton temperature, is the number density of electron, is the Boltzmann constant, is the electron mass, is the speed of light, and is the cross section of electron scattering. We adopt as our fiducial value for observed low-luminosity AGN (Xie, Yuan & Ho, 2017). Note that Compton temperature for low-luminosity AGN () is higher than that estimated in previous works using composite spectra of high luminosity AGNs; (Sazonov et al., 2004). The functional forms of and 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


We estimate the bolometric AGN luminosity as , where is the gas accretion rate through the inner-most grid and is the radiative efficiency.


where . The fitted values are , , (), , , and (). 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 is obtained by a semi-analytical model (green dashed curve; Xie & Yuan, 2012). For lower rates of , 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 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 and , where is set to to avoid numerical singularity at poles. In the radial direction, we set up logarithmically-spaced grids, the number of which is . 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; at (near the mid-plane) and otherwise (near the poles). We set and so that the disk scale height is sufficiently resolved. As our fiducial case, we set and . The grid structure of our simulations is illustrated in Fig. 3.

As initial conditions, we adopt a rotational equilibrium distribution () with a constant specific angular momentum of :


(Fishbone & Moncrief, 1976; Papaloizou & Pringle, 1984), where is the cylindrical radius, is the ambient gas density, and 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 , which leads to a maximum value of the density at and ,


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


Note that the corresponding specific angular momentum is given by . Without viscosity, the density never exceeds this value because of the centrifugal barrier. In other words, a high-density region with must be formed by angular momentum transport due to viscosity. In Fig. 3, the density distribution of the equilibrium torus for is shown.

Figure 2: Radiative efficiency models for a hot accretion disk as a function of the BH feeding rate in units of . Dashed (green) curve and symbols (orange square) show the efficiency obtained from a semi-analytical model (Xie & Yuan, 2012), and by general relativity and radiation transfer simulations (Ryan et al., 2017). Solid (blue) curve presents the values we adopt in this work (linear interpolation of the log values at ). Dotted (red) curve is an ADAF-like model (e.g., Ciotti et al., 2009), which is often used in previous work.

Although we employ static mesh refinement near the equatorial region, a disk scale height cannot be resolved well if the gas temperature drops below , where 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., . We adopt as our fiducial value and study two cases with and . As discussed in §4.2, the Toomre parameter depends on the choice of somewhat sensitively. We will conduct an extrapolation of our results for lower values of 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 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 . After the transition to a cold disk, we assume an isothermal Keplerian disk in the inner ghost grids at , where 111 When the gas density is not high enough to permit cooling, the temperature is higher than the initial value of .; in practice, the rotational velocity is set to .

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., and )222Note that both of and are solved self-consistently in our simulations, but fixed to evaluate the Bondi rate.,


where for , , and . The accretion rate normalized by the Eddington rate is given by


where is the Eddington accretion rate with a radiative efficiency.

In this work, we set the BH mass to , the temperature of the ambient gas to () at the initial state and and as our fiducial case. We study the effects of radiative cooling and heating, varying the gas density and thus .

Figure 3: Density distribution of the initial equilibrium torus given by Eq. (14) for . The grid structure for our simulations is shown. In the polar direction, we use static refinement at (near the mid-plane) in order to resolve the scale height of a cold disk (see §3).

3 Results

Fig. 4 shows the time evolution of gas accretion rates through the inner-most grid () for four different values of (red), (green), (magenta) and (blue), each of which corresponds to , , and . The accretion rates are shown in units of (top) and are normalized by the Bondi accretion rate with for each case (bottom). The simulation time is normalized by the dynamical timescale at the Bondi radius, . For the lowest density ( and ), the accretion rate increases and saturates at . The saturated value is as small as , which means that the BH accretion rate is significantly reduced from the Bondi rate. As long as , 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 ( and ), 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 (), which is times higher than the Bondi rate for (note that for the case with , a steady state is realized at ). 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, (see discussion in §4.3).

Figure 4: Time evolution of the accretion rates for four different values of (red), (green), (magenta) and (blue), which correspond to , , and . The horizontal lines show the Bondi accretion rate for adiabatic gas with (short dashed) and for isothermal gas with (long dashed). The accretion rates are shown in units of (top) and are normalized by the Bondi rates with (bottom). While for lower densities (), the accretion rates are significantly reduced from the Bondi rate, for higher densities (), the accretion rate continues to increase and ultimately reaches the Bondi rate for isothermal gas, which is times higher than that for adiabatic gas (note that for the case with , a steady state is realized at .). Open circles mark the epochs at which the radial and angular structure of the density and temperature are shown in Figs. 8 and 9.
Figure 5: Radial structures of the inflow rates (dashed) and net accretion rates (solid) for (red), (green), and (blue). The accretion rates are normalized by the Bondi accretion rates for adiabatic gas with . The elapsed times correspond to those when the accretion flows are in steady states (). For lower-density cases, the inflow rates decrease toward the center, following (black dotted line). The net accretion rates become a constant value of within the Bondi radius. For the highest-density case, the net accretion rate is almost comparable to the inflow rate and depends on the radius weakly (; at ).
Figure 6: Distribution of the gas density (left) and temperature (right) for the case with the highest density, . At the early stage (left panel), the accretion flow is highly turbulent and forms a geometrically-thick disk inside the Bondi radius. In the late stage (right panel), the accreting gas begins to cool, collapses towards the equator and forms a geometrically-thin disk.
Figure 7: Density distribution and velocity vectors at three elapsed times around over a short duration with for the model with . A significant fraction of outflowing material to the polar regions cannot escape and forms clumps around due to compression with the ambient gas (left). Those clumps are captured by the gravity of the BH and accrete onto the disk (middle). Outflows in the bipolar directions are produced by energy release from the disk and BH.
Figure 8: Radial profile of the gas density (top) and temperature (bottom) for the case with the highest density, . The curves show the profiles at (black dashed), (red), (green), (blue), (magenta), (cyan), and (orange).
Figure 9: Angular profile of the gas density (top) and temperature (bottom) for the case with the highest density, in the final time step (). Different curves show the profiles at different radial positions of . The cold disk is located within two vertical lines (orange) where grids in the polar direction are refined.

Fig. 5 presents the radial structure of the angle-integrated mass inflow (dashed) and net accretion (solid) rates for three different values of (). The inflow and outflow rate are calculated as


and the net accretion rate is defined as . Note that and . For the two cases with lower densities ( and ), the inflow rates decrease towards the center, following (red and green dashed curves). The power-law index is consistent with that in CDAFs (). 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 within the Bondi radius. For the highest-density case (), 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 , is two orders of magnitude higher than those in the lower-density cases. The radial dependence of the inflow rate becomes flatter () 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 . The velocity vectors are over plotted in the temperature distribution. At the early stage of (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 (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 (). As shown in Fig 4, the accretion flow settles into a steady state at , where the BH accretion rate is as high as . In the steady state, the radiative luminosity from the BH () is on the order of . The combination of heating and outward momentum caused by radiation produces weak outflows at . A small fraction of the gas can escape from the Bondi radius with a velocity of , 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 over a short duration with for the model with . 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 . 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 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 and 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 (; see §4.1). After the transition, the gas density in the disk dramatically increases to at , and the disk outer edge moves outwards by angular momentum transport (see epoch ). The gas temperature drops to the threshold value we set () via optically-thin cooling. Note that this transition occurs in an unstable way because the cooling timescale becomes shorter with decreasing temperature as , where at due to free-free emission, metal lines and recombination. In the steady state, the density profile at the mid-plane follows


where is the accretion rate through the disk, comparable to the Bondi rate for isothermal gas, and is the sound speed of the gas in the disk with a temperature of . The solid line in the top panel shows Eq. (21) at . The temperature exterior to the cold disk, on the other hand, becomes hotter due to Compton heating () 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 , at radial positions of . Inside the Bondi radius, the density increases toward the center and the mid-plane (). Within the cold disk (), the temperature is constant at . 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 ().

4 Physical interpretation

4.1 The formation of a cold accretion disk

For lower densities at the Bondi radius (), 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: , and (e.g., Quataert & Gruzinov, 2000). Inside the torus structure, viscous energy dissipation heats the gas and drives convective motion on a thermal timescale of . With the CDAF solution, the bremsstrahlung cooling rate per volume is and the cooling timescale is . Fig. 10 shows the ratio of as a function of for different values of and , 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


where the ratio is evaluated at the outer edge of the torus (). We note that this scaling relation agrees with the results shown in Fig. 10, as long as the angular momentum is sufficiently small (). From Eq. (22), we obtain the critical accretion rate above which bremsstrahlung cooling plays an important role in the accretion flow as


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 () and ambient gas temperature (), the transition from adiabatic accretion flows to cold disk accretion occurs at . We note that Eq. (23) is not valid for lower temperatures where metal-line cooling dominates bremsstrahlung cooling significantly ().

Figure 10: The ratio of the viscous heating timescale to the radiative cooling timescale as a function of , for different values of the angular momentum () and ambient gas temperature (). The shaded region at presents the transition region from adiabatic flows to cold disk accretion. For and , we show our numerical results where adiabatic accretion is maintained (open) and where cold accretion disk forms (filled).

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 (), where the sum of viscous heating and compressional heating balance with cooling via Coulomb collisions between ions and electrons. For higher rates of 333 The critical rate would depend on the choice of the outer edge of the disk. For , Yuan (2001) obtained , and for , and , respectively., hot accretion flows no longer exist because radiative cooling efficiently carries the energy away. We note that the critical value of is two orders of magnitude higher than the critical rate given by Eq. (23) assuming . 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 collapses via cooling. Wu et al. (2016) have also studied the cooling transition of LHAFs at , performing hydrodynamical simulations. They found a critical accretion rate of , 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 , above which a steady state no longer exists due to feedback but gas accretion occurs episodically. The exact value of is still uncertain, depending on which feedback processes dominate. Instead of discussing each process, we here adopt a critical radiation luminosity of as a reference (e.g., McClintock & Remillard, 2006). Since the BH accretion rate through a cold (isothermal) disk is times higher than the Bondi accretion rate estimated with , the critical value of the Bondi rate is , where (see Fig. 12 and discussion in §4.3). We leave the study of the strong feedback domain for future work.

Figure 11: Radial profiles of the Toomre parameter for at the final time step () with an inner-most radius of (fiducial case, solid) and (dashed). Each curve presents the profile for a minimum temperature of (green), (red) and (blue), respectively. The minimum values are located at . Star symbols indicate those shown in Eq. (25).

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

We now focus on the intermediate regime of , 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):


Assuming a steady state, the gas surface density is estimated as , where is the gas accretion rate through the disk. Thus, we obtain


where . Fig. 11 shows the radial profile of the Toomre parameter for (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 at the centrifugal radius (). 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 (red), (green) and (blue), but with an inner-most radius of , which is twice larger than our fiducial case (solid). The Toomre values at decreases for lower , following . Therefore, we would expect that once radiative cooling leads to the formation of a cold disk and the temperature decreases below , the disk would become unstable because of its self-gravity.

When the Toomre parameter decreases to , 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 . To mimic this effect, we parameterize the effective viscosity as a function of in the following:


where is the viscous parameter induced by the MRI and 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 . 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 (e.g., Kratter et al., 2008)444Numerical simulations for isolated disks without infalling flow from larger scales suggest lower critical values of (Rice et al., 2005).. Based on those results, we adopt as a reference value. We also note that the choice of the two free parameters of hardly affects the strength of gravitational torque (see Appendix A) and thus they are set to . Using Eqs. (25) and (26), we obtain the critical accretion rate above which the disk becomes gravitationally unstable ( and ) as




where is adopted as the threshold value. The minimum temperature is as high as , 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 , it would turn gravitationally unstable immediately. For less massive BHs with , accretion rates as high as 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 ( and ; 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 , where and (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 . In the marginally unstable disk, the surface density and disk mass are given by




where is adopted. A fraction of the disk mass would be consumed by star formation. Assuming the star formation efficiency inferred from the Kennicut-Schmidt law (Schmidt, 1959; Kennicutt, 1998; Krumholz & McKee, 2005), the total stellar mass is and the expected number of supernovae (SNe) is , where is the stellar mass required to generate an SN for a Kroupa initial mass function (Kroupa, 2001). Since the total SN energy erg is injected over the lifetime of massive stars Myr, thus the kinetic luminosity is . 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 , where is the ambient gas density, is the turbulent velocity in the disk, and is the total surface area of the bubbles (Mac Low & McCray, 1988; Koo & McKee, 1992; Kim et al., 2017). In this case for M87 (, and cm at ), the ratio of the SN kinetic luminosity to that of turbulent gas is estimated as . 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 fed at , 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 emission, has a size of (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 (). The star formation rate for M87 is estimated as based on the 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 , implying that a large fraction of the gas supplied from the Bondi radius is consumed by star formation, namely, . 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 (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 (Tojeiro et al., 2013). In addition, most nuclear disks in the centers of nearby ellipticals seem stable against star formation, i.e., (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 () produced from accretion flows and BH feeding rate () for a given Bondi accretion rate at a range of . The results are summarized in Fig. 12.

Below the cooling threshold (), 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


where 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 (e.g., Narayan & Medvedev, 2001; Maron et al., 2004). For such low BH feeding rates, the radiative efficiency decreases with the accretion rate at and has been estimated with radiation magneto-hydrodynamical simulations (Ryan et al., 2017, , see also orange squares in Fig. 2). Therefore, adopting the radiation efficiency, the radiation luminosity model in this domain can be estimated as


which is shown by the red shaded region in Fig. 12. The width of the region reflects the uncertainties of the conductivity suppression factor . We note that if the BH feeding rate is equal to the Bondi accretion rate (