Low-density, radiatively inefficient rotating-accretion flow on to a black hole

Low-density, radiatively inefficient rotating-accretion flow on to a black hole Abstract We study low-density axisymmetric accretion flows on to black holes (BHs) with two-dimensional hydrodynamical simulations, adopting the α-viscosity prescription. When the gas angular momentum is low enough to form a rotationally supported disc within the Bondi radius (RB), we find a global steady accretion solution. The solution consists of a rotational equilibrium distribution around r ∼ RB, where the density follows ρ ∝ (1 + RB/r)3/2, surrounding a geometrically thick and optically thin accretion disc at the centrifugal radius RC(<RB), where thermal energy generated by viscosity is transported via convection. Physical properties of the inner solution agree with those expected in convection-dominated accretion flows (ρ ∝ r−1/2). In the inner solution, the gas inflow rate decreases towards the centre due to convection ($$\dot{M}\propto r$$), and the net accretion rate (including both inflows and outflows) is strongly suppressed by several orders of magnitude from the Bondi accretion rate $$\dot{M}_{\rm B}$$. The net accretion rate depends on the viscous strength, following $$\dot{M}/\dot{M}_{\rm B}\propto (\alpha /0.01)^{0.6}$$. This solution holds for low accretion rates of $$\dot{M}_{\rm B}/\dot{M}_{\rm Edd}\lesssim 10^{-3}$$ having minimal radiation cooling, where $$\dot{M}_{\rm Edd}$$ is the Eddington accretion rate. In a hot plasma at the bottom (r < 10−3 RB), thermal conduction would dominate the convective energy flux. Since suppression of the accretion by convection ceases, the final BH feeding rate is found to be $$\dot{M}/\dot{M}_{\rm B}\sim 10^{-3}$$–10−2. This rate is as low as $$\dot{M}/\dot{M}_{\rm Edd}\sim 10^{-7}$$–10−6 inferred for SgrA* and the nuclear BHs in M31 and M87, and can explain their low luminosities, without invoking any feedback mechanism. accretion, accretion discs, black hole physics, X-rays: galaxies 1 INTRODUCTION The basic physics of gas accretion on to an astrophysical object was first studied for a spherically symmetric flow (Bondi 1952). When rotation is not included, gas accretion on to a black hole (BH) with a mass of M• begins from a characteristic radius, where the negative gravitational energy becomes greater than the thermal energy of the gas. The so-called Bondi radius is given by   \begin{equation} R_{\rm B}\equiv \frac{GM_\bullet }{c_{\rm s}^2}, \end{equation} (1)where cs is the sound speed of the gas. One can safely assume roughly free fall within this radius at a rate of $$\dot{M}_{\rm B}$$, so-called the Bondi rate. If radiation emitted from the inner region of the flow exceeds several percent of the Eddington luminosity (LEdd ≡ 4πcGM•/κes, where κes is the electron scattering opacity, then this feedback can add enough energy and momentum at the Bondi radius to reverse the flow (Ostriker et al. 1976; Ciotti & Ostriker 2001; Proga 2007; Ciotti, Ostriker & Proga 2009; Milosavljević et al. 2009; Park & Ricotti 2011; Inayoshi, Haiman & Ostriker 2016). However, the problem is rather simple for low accretion rates (Shapiro 1973; Park 1990). Now let us add angular momentum per unit mass j to the flow. This effect introduces a new physical scale of   \begin{equation} R_{\rm C}\equiv \frac{j^2}{GM_\bullet }, \end{equation} (2)which is the centrifugal radius where the centrifugal force balances the gravity of the BH. If this radius is larger than the Schwarzschild radius   \begin{equation} R_{\rm Sch}\equiv \frac{GM_\bullet }{c^2}, \end{equation} (3)the accreting gas will simply settle into a rotational equilibrium distribution with no net inflow (Fishbone & Moncrief 1976; Papaloizou & Pringle 1984; Li, Ostriker & Sunyaev 2013). However, if viscosity exists, then angular momentum transport is possible and an accretion disc forms at the bottom of the distribution. The effective shear viscosity driven by magneto rotational instability (MRI) is often described with the standard α-prescription (e.g. Shakura & Sunyaev 1973) as   \begin{equation} \nu =\alpha c_{\rm s}H, \end{equation} (4)where H is the scale height of the inner flow and α presents the strength of viscosity estimated as α ∼ O(10−2) by magnetohydrodynamical (MHD) simulations (Balbus & Hawley 1991; Matsumoto & Tajima 1995; Stone et al. 1996; Balbus & Hawley 1998; Sano et al. 2004). There have been a number of analytical solutions and numerical studies for rotating flows with viscous angular momentum transport (e.g. Shakura & Sunyaev 1973; see reviews by Pringle 1981; Kato, Fukue & Mineshige 2008, references therein). Among those, we here focus on accretion flows which cannot lose internal energy via radiative cooling because of very low gas density. Such radiatively inefficient accretion flows are quite interesting since many observed BHs accrete at rates of only a small fraction of the Bondi accretion rate and their radiation luminosity is as low as ∼10−1 to 10−9 LEdd (Ho 2008, 2009). Sagittarius A* (Sgr A*) is inferred to be accreting at a rate of 10−3 to $$10^{-2}\ \dot{M}_{\rm B}$$ (e.g. Yuan, Quataert & Narayan 2003; Quataert 2004), where $$\dot{M}_{\rm B}\simeq 10^{-5}\ {\rm M}_{\odot }\, {\rm yr}^{-1}$$ is measured from the temperature and density near the Bondi radius with X-ray observations (Baganoff et al. 2003). Because of such a low accretion rate, the bolometric luminosity of Sgr A* (M• ≃ 4 × 106 M⊙; Ghez et al. 2003) is as small as Lbol ∼ 1036 erg s−1 ∼ 2 × 10−9 LEdd. The second example is a BH at the centre of the giant elliptical galaxy M87. The gas accretion rate at the vicinity of the BH is estimated as ≲ 9.2 × 10− 4 M⊙ yr− 1 (Kuo et al. 2014), which is lower than $${\sim } 10^{-2}\ \dot{M}_{\rm B}$$ (Russell et al. 2015). Since the BH mass is estimated as $$M_\bullet = 6.6^{+0.4}_{-0.4}\times 10^9\ {\rm M}_{\odot }$$ (Gebhardt et al. 2011) and $$M_\bullet = 3.5^{0.9}_{-0.7}\times 10^9\ {\rm M}_{\odot }$$ (Walsh et al. 2013), the bolometric luminosity of Lbol ∼ 2 × 1041 erg s−1 is ∼3 × 10−7 LEdd. The third example is a BH at the centre of the Andromeda Galaxy (M31). The estimated BH mass is $$M_\bullet \simeq 1.4^{+0.7}_{-0.3}\times 10^8\ {\rm M}_{\odot }$$ (Bender et al. 2005). The Bondi accretion rate and the X-ray luminosity are estimated as $$\dot{M}_{\rm B}\simeq 7\times 10^{-5}\ {\rm M}_{\odot }\, {\rm yr}^{-1}$$ and LX ≃ 2 × 1036 erg s−1 ≃ 10−10 LEdd, respectively (Garcia et al. 2010). The corresponding bolometric luminosity is inferred as ≃10−9 LEdd by assuming the bolometric correction factor to be ≃10 (Hopkins, Richards & Hernquist 2007). In very low-accretion-rate flows, the gas never cools via radiation and the accretion disc becomes hot and thick. This thick-disc solution has been found by Ichimaru (1977) and studied in the subsequent works by Narayan & Yi (1994, 1995). This solution is well known in a simplified self-similar version as advection-dominated accretion flows (ADAFs). Their physical properties can be understood by considering a self-similar solution, which suggests that low-density accretion flow is hot and geometrically extended in the polar direction and the thermal energy is advected with the flow on to the BH. Moreover, as a related solution, the so-called adiabatic inflow–outflow solution has been proposed by Blandford & Begelman (1999, 2004), where matter can accrete at very low rates due to strong outflows which carry energy and angular momentum away. As pointed out by Narayan & Yi (1994, 1995), ADAF solutions tend to be convectively unstable because the gas entropy increases towards the centre due to thermal energy release via viscosity. Narayan, Igumenshchev & Abramowicz (2000) and Quataert & Gruzinov (2000) conducted stability analyses for rotating adiabatic gas against convection motions, and found that the density profile in a marginally stable state follows ρ ∝ r−1/2, whose slope is rather flatter than ρ ∝ r−3/2 expected from ADAFs. This solution is known as convection-dominated accretion flows (CDAF). Numerical simulations of radiatively inefficient accretion flows have also supported that CDAF solutions can exist in a weak viscosity and no cooling regime (Igumenshchev & Abramowicz 1999, 2000; Stone, Pringle & Begelman 1999; Igumenshchev, Narayan & Abramowicz 2003). Many previous works with multidimensional MHD simulations have studied the properties of the accretion flow at the vicinity of the central BH (r ≲ a few × 102 RSch) starting from torus-like initial reservoirs of gas (Hawley, Balbus & Stone 2001; Machida, Matsumoto & Mineshige 2001; Stone & Pringle 2001; Igumenshchev, Narayan & Abramowicz 2003; McKinney & Gammie 2004; Ohsuga et al. 2009; Narayan et al. 2012; see also Yuan, Wu & Bu 2012, who studied the gas dynamics over a wider range of spacial scales). However, because of the limitation of the computational domain, it remains unclear how and whether the accretion flow, which begins from large radii (r ∼ RB), can lose the angular momentum and reach smaller radii (r ≪ RB), where the gas is tightly bound to the central BH. To explicitly connect large and small scales in a self-consistent solution, we study dynamics of gas accretion on to a BH over a wide range of spatial scales, and are able to connect the inner region, where a disc forms, with the region well outside the Bondi radius, where the accretion originates in the first place. Li et al. (2013) studied a similar problem for one specific parameter set, i.e. RC/RB = 0.02 and α ∼ 10−3, with two-dimensional hydrodynamical simulations assuming equatorial symmetry. This work followed gas accretion from outside the Bondi radius to within the centrifugal radius. Proga & Begelman (2003) have conducted two-dimensional MHD simulations for this problem with a more limited dynamic range (r < 103 RSch), and found that gas accretion is allowed by viscosity driven by the MRI. However, several crucial questions remain unanswered. How do physical parameters, e.g. RC/RB and α, determine the actual inflow rate on to a BH? What do physical properties of the accretion flow inside the centrifugal radius look like? We perform two-dimensional viscous-hydrodynamical simulations with a large dynamic range including angular momentum transport with the α-viscosity prescription, and address the above questions for radiatively inefficient rotating-accretion flows on to a BH. Here, we focus on cases where the gas angular momentum is low enough to form a rotationally supported disc/torus within the Bondi radius, i.e.   \begin{equation} R_{\rm Sch}\ll R_{\rm C}<R_{\rm B}. \end{equation} (5)In this case, we find a global steady accretion solution which consists of three parts. In the outer region (r ∼ RB), a rotational equilibrium distribution is developed where the density distribution follows ρ ∝ (1 + RB/r)3/2 and the gas motion is very subsonic. The solution deviates only slightly from the rotational equilibrium flow with no radial motion (Fishbone & Moncrief 1976). Interior to this (r ≲ 2 RC), a geometrically thick accretion torus is formed, where thermal energy generated by viscosity is transported via strong convection motions. Physical properties of the inner solution agree with those expected in a CDAF solution, e.g. ρ ∝ r−1/2 (Narayan et al. 2000; Quataert & Gruzinov 2000). In the inner CDAF solution, the gas inflow rate decreases towards the centre due to convection ($$\dot{M}\propto r$$), and the net accretion rate (including both inflows and outflows) is strongly suppressed by several orders of magnitude from the Bondi accretion rate. We also study the effect of the viscous strength on the results and find that the net accretion rate is approximated as $$\dot{M}\sim (\alpha /0.01)^{0.6}(r/R_{\rm B})\dot{M}_{\rm B}$$. Finally, in a hot plasma at the bottom of this solution (r ≲ 10− 3 RB), thermal conduction would dominate the convective energy flux and the flow would resemble an optically thin accretion disc. Since suppression of the accretion by convection ceases, the final BH feeding rate is determined as $$\dot{M}\sim 10^{-3}-10^{-2}\ \dot{M}_{\rm B}$$. This rate can explain why Sgr A* and the BHs in M31 and M87 accrete at such low accretion rates, without invoking any feedback mechanism. The rest of this paper is organized as follows. In Section 2, we describe the methodology of our numerical simulations. In Section 3, we show our simulation results and explain their physical properties. We compare the results to those in CDAF solutions and quantify the effects of energy and angular momentum transport by convection (Section 3.1), and study the effect of the choice of viscous parameters on the results (Section 3.2). In Section 4, we give an analytical argument extending the accretion solution further inwards and discuss a global solution of radiatively inefficient rotating-accretion flows. In Section 5, we summarize our main conclusions and discuss the importance for observations of very low-luminosity BHs. 2 METHODOLOGY 2.1 Basic quantities We introduce basic dimensionless physical quantities which characterize BH accretion systems. If feedback by radiation and/or momentum from the central BH is negligible, gas accretion begins from the Bondi radius, RB. The standard accretion rate from the Bondi radius is given by   \begin{eqnarray} \dot{M}_{\rm B}&=& 4\pi q(\gamma ) \rho _\infty \frac{G^2M_\bullet ^2}{c_{\infty }^3}, \nonumber \\ &\simeq & 8.8\times 10^{-5}\ \rho _{-22} M_6^2 c_7^{-3}\ {\rm M}_{\odot }\, {\rm yr}^{-1}, \end{eqnarray} (6)where q(γ) = 1/4 for γ = 5/3, $$\rho _{-22}=\rho _\infty /(10^{-22}\ {\rm g\ {\rm cm}^{-3}})$$, M6 = M•/(106 M⊙) and c7 = c∞/(107 cm s−1) and the accretion rate normalized by the Eddington rate is given by   \begin{eqnarray} \dot{m}_{\rm B}\equiv \frac{\dot{M}_{\rm B}}{\dot{M}_{\rm Edd}} =3.8\times 10^{-3}\ \rho _{-22} M_6 c_7^{-3}, \end{eqnarray} (7)where $$\dot{M}_{\rm Edd}(\equiv 10\ L_{\rm Edd}/c^2) = 2.3\times 10^{-2}\ M_6\ {\rm M}_{\odot }\, {\rm yr}^{-1}$$ is the Eddington accretion rate with a 10 per cent radiative efficiency. As shown in Section 1, there are three characteristic physical scales: the Bondi radius RB, the centrifugal radius RC and the Schwarzschild radius RSch. From these three radii, we can define two ratios as   \begin{equation} \frac{R_{\rm Sch}}{R_{\rm B}}=\frac{2c_\infty ^2}{c^2}, \end{equation} (8)and   \begin{equation} \frac{R_{\rm C}}{R_{\rm B}}=\frac{j^2c_\infty ^2}{G^2 M_\bullet ^2}. \end{equation} (9)Assuming a constant specific angular momentum of $$j=\sqrt{\beta }R_{\rm B}c_\infty$$, the ratio is written as RC/RB = β. Since we discuss cases where RSch ≪ RC < RB, we assume β < 1. In axisymmetric two-dimensional hydrodynamical simulations, we need angular momentum transport via viscous processes in order to allow gas accretion on to a BH at the centre. We here model the effects of viscosity with the standard α-prescription proposed by Shakura & Sunyaev (1973) instead of solving the time evolution of MHD equations. According to MHD simulations, effective viscosity driven by the MRI provides angular momentum transport and the strength is estimated as α ∼ O(10−2) (Balbus & Hawley 1991; Matsumoto & Tajima 1995; Stone et al. 1996; Balbus & Hawley 1998; Sano et al. 2004). Throughout this paper, we focus on adiabatic gas without radiative cooling in order to keep numerical results scale-free. The assumptions are justified for very low accretion-rate system, i.e. $$\dot{M}_{\rm B}\ll \dot{M}_{\rm Edd}$$ (see discussion in Section 4). In this limit, we have only three important non-dimensional parameters of RSch/RB, RC/RB( = β) and α to define a BH accretion system. 2.2 Basic equations We solve the axisymmetric two-dimensional hydrodynamical equations using an open code pluto (Mignone et al. 2007; Kuiper et al. 2010, 2011). The basic equations are following: the equation of continuity,   \begin{equation} \frac{{\rm d}\rho }{{\rm d}t}+\rho \nabla \cdot \boldsymbol v=0, \end{equation} (10)and the equation of motion,   \begin{equation} \rho \frac{{\rm d}\boldsymbol v}{{\rm d}t}=-\nabla p -\rho \nabla \Phi + \nabla \cdot \boldsymbol \sigma , \end{equation} (11)where ρ is the density, $$\boldsymbol v$$ is the velocity, and p is the gas pressure, the gravitational potential is set to Φ = −GM•/r, and $$\boldsymbol \sigma$$ is the stress tensor due to viscosity. The time derivative is the Lagrangian derivative, given by d/dt ≡ ∂/∂t + $$\boldsymbol v$$· ∇. We solve the energy equation of   \begin{equation} \rho \frac{{\rm d}e}{{\rm d}t}=-p\nabla \cdot \boldsymbol v + (\boldsymbol \sigma \cdot \nabla ) \boldsymbol v, \end{equation} (12)where e is the internal energy per mass. The equation of state of the ideal gas is assumed as p = (γ − 1)ρe, where γ = 1.6.1 The two terms on the right-hand side present compressional heating (or expansion cooling) and viscous heating. The viscous stress tensor is given by   \begin{equation} \sigma _{ij}=\rho \nu \left[ \left(\frac{\mathrm{\partial} v_j}{\mathrm{\partial} x_i} + \frac{\mathrm{\partial} v_i}{\mathrm{\partial} x_j}\right) -\frac{2}{3} (\nabla \cdot \boldsymbol v)\delta _{ij} \right], \end{equation} (13)where ν is the shear viscosity. Note that the bulk viscosity is neglected here. The shear viscosity is calculated with the α-prescription (Shakura & Sunyaev 1973),   \begin{equation} \nu = \alpha \frac{c_{\rm s}^2}{\Omega _{\rm K}}, \end{equation} (14)where ΩK = (GM•/r3)1/2. We calculate the viscous parameter by mimicking some properties of the MRI as   \begin{equation} \alpha = \alpha _0 \exp {\left[ -\left(\frac{\rho _{\rm cr}}{\rho }\right)^2\right]}, \end{equation} (15)where ρ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 Section 2.3). Under this model, the viscous process is active almost only within an accretion disc (r ≲ RC), where the rotational velocity has a significant fraction of the Keplerarian value. On the other hand, since the viscosity becomes zero outside the disc, angular momentum transported from the disc tends to be accumulated outside it as shown by the ‘bump’ in the angular momentum distribution in Fig. A4, where no viscous processes operate. Exterior to the peak in the angular momentum distribution where the specific angular momentum has a negative gradient outwards, i.e. ∂j/∂r < 0, so-called Rayleigh's criterion. In reality, such rotating flows are unstable and become turbulent (e.g. Chandrasekhar 1961), leading to angular momentum transport in three-dimensional simulations. However, because of limitations of our two-dimensional simulations, this could not occur and angular momentum flowing out from the central regions would accumulate outside of the region where ρ/ρcr ∼ 1. Thus, we modify equation (15) by adding the second term as   \begin{equation} \alpha = \alpha _0 \left\lbrace \exp {\left[ -\left(\frac{\rho _{\rm cr}}{\rho }\right)^2\right]} + {\rm max}\left(0, -\frac{\mathrm{\partial} \ln j}{\mathrm{\partial} \ln r}\right)\right\rbrace , \end{equation} (16)In a physical sense, the second term is necessary so that a steady state of the accretion flow exists. Note that the treatment of the rotational instability does not affect our results (see Appendix). 2.3 Boundary and initial conditions We set a computational domain of rmin ≤ r ≤ rmax and ε ≤ θ ≤ π − ε, where ε is set to 0.01 radian to avoid numerical singularity at poles. We set logarithmically spaced grids in the radial direction and uniformly spaced grids in the polar direction. The number of the grid points is set to (Nr, Nθ) = (512, 256). We have checked convergence for our simulation results, changing the number of the grids (see Appendix). As our fiducial case, we set rmin = 7.7 × 10−3 RB and rmax = 54 RB. Moreover, we calculate two cases with rmin = 3.0 × 10−3 RB and 1.3 × 10−2 RB (see Sections 3.1.2 and 3.2). Note that the relative size of each cell is given by $$\Delta r/r=(r_{\rm max}/r_{\rm min})^{1/N_r}$$. We adopt the outflow boundary condition at the innermost/outermost grid (e.g. 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. We also impose vr ≤ 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θ. As initial conditions, we adopt a rotational equilibrium distribution (vr = vθ = 0) with a constant specific angular momentum of j∞,   \begin{equation} \frac{\rho }{\rho _\infty }=\left[1+(\gamma -1)\frac{GM_\bullet }{c_\infty ^2 r} - \frac{(\gamma -1)}{2}\frac{j_\infty ^2}{c_{\rm \infty }^2\varpi ^2} \right] ^{\frac{1}{\gamma -1}} \end{equation} (17)(Papaloizou & Pringle 1984; see also Fishbone & Moncrief 1976), where ϖ = r sin θ is the cylindrical radius. Hereafter, we refer to this as the ‘Fishbone–Moncrief’ solution. 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∞, which leads to a maximum value of the density at r = RC and θ = 0   \begin{equation} \rho _{\rm cr}=\rho _\infty \left(1+\frac{\gamma -1}{2\beta }\right)^{1/(\gamma -1)}, \end{equation} (18)where β = RC/RB as defined below equation (9). Without viscosity, the density never exceeds this value because of the centrifugal barrier. In other words, a high-density region with ρ > ρcr must be formed by angular momentum transport due to viscosity. In order to study the dependence of our numerical results on the choice of the initial conditions, we also conduct a simulation which starts from the Bondi accretion solution with a constant angular momentum j∞. Note that the numerical result with the Bondi profiles as initial conditions approaches that adopting a rotational equilibrium distribution given by equation (17) (see Appendix). 3 RESULTS In this section, we show results of our two-dimensional simulations for accretion flows. We set all physical quantities as M• = 106 M⊙, $$\rho _\infty = 10^{-22}\ {\rm g\ {\rm cm}^{-3}}$$, and c∞ = 100 km s−1 ($$\dot{m}_{\rm B}\simeq 3.8\times 10^{-3}$$). As already mentioned, in adiabatic cases, our results do not depend on $$\dot{m}_{\rm B}$$ but on the two dimensionless parameters of β = RC/RB and α. In Section 3.1, we first describe overall properties of accretion flows with different initial angular momentum β for our fiducial case with α = 0.01. In Section 3.2, the dependence of our results on the choice of viscous parameter is discussed. 3.1 Fiducial case with α = 0.01 Fig. 1 shows the time evolution of gas accretion rates on to a BH (i.e. a sink cell at r = rmin) for different initial values of angular momentum: RC/RB = 0.1 (red), 0.2 (green), and 0.3 (blue). The viscous parameter is set to α = 0.01. We normalize the simulation time by the dynamical time-scale at the Bondi radius, $$t_{\rm dyn}(\equiv R_{\rm B}/c_\infty ) \simeq 1.3\times 10^{11}\ M_6 c_7^{-3}\ {\rm s}$$. For all the cases, the accretion rates increase with fluctuations at t ≲ 4 tdyn and saturate at almost constant values by t ≃ 5 tdyn. The behaviour of the accretion rates in all the cases is similar. In fact, the time-averaged values over 6 ≤ t/tdyn ≤ 16 are $$\dot{M}/\dot{M}_{\rm B}\simeq 7\times 10^{-3}$$, which does not depend on the initial angular momentum. The results show that gas rotation reduces the inflow rates by two orders of magnitude from the Bondi rate. Note that the suppression factor of the accretion rate is close to the viscous parameter we assume (α = 0.01). Of course in the limit of α = 0, there is no net accretion. Figure 1. View largeDownload slide Time evolution of the net accretion rate (in units of the Bondi rate) on to a BH, i.e. a sink cell for different values of angular momentum: RC/RB = 0.1 (red), 0.2 (green), and 0.3 (blue). The viscous parameter is set to α = 0.01. The time-averaged accretion rates (t > 6 tdyn) are $$\dot{M}\simeq 7\times 10^{-3}\ \dot{M}_{\rm B}$$. Figure 1. View largeDownload slide Time evolution of the net accretion rate (in units of the Bondi rate) on to a BH, i.e. a sink cell for different values of angular momentum: RC/RB = 0.1 (red), 0.2 (green), and 0.3 (blue). The viscous parameter is set to α = 0.01. The time-averaged accretion rates (t > 6 tdyn) are $$\dot{M}\simeq 7\times 10^{-3}\ \dot{M}_{\rm B}$$. Fig. 2 shows the radial structure of the angle-integrated mass inflow (red) and outflow (green) rates for the case with RC/RB = 0.1 and α = 0.01 at t ≃ 13 tdyn. Those rates are defined by   \begin{eqnarray} \dot{M}_{\rm in}(r)=2\pi r^2\int \rho \cdot {\rm min}(v_r,0)\sin \theta\, {\rm d}\theta , \end{eqnarray} (19)  \begin{eqnarray} \dot{M}_{\rm out}(r)=2\pi r^2\int \rho \cdot {\rm max}(v_r,0)\sin \theta\, {\rm d}\theta , \end{eqnarray} (20)where $$\dot{M}_{\rm in}< 0$$ and $$\dot{M}_{\rm out}> 0$$ (e.g. Stone et al. 1999). Note that $$\dot{M}_{\rm out}$$ in equation (20) does not necessarily mean the outflow rates of gas escaping to the infinity, but takes account into gas flows with vr > 0. We also show the net accretion rate defined by $$-\dot{M}\equiv -\dot{M}_{\rm in}-\dot{M}_{\rm out}$$ (blue). Inside the Bondi radius (r ≲ RB), both inflows and outflows exist and those rates, which decrease towards the centre following ∝ r, are tightly balanced. Within the centrifugal radius (r ≲ RC = 0.1 RB), the net accretion rate becomes a constant value ($${\simeq } 7\times 10^{-3}\ \dot{M}_{\rm B}$$). We note that the outflow rate dominates the inflow rate outside the Bondi radius. In Appendix, we study the time evolution of the outflowing gas and make sure that this component does not affect stationarity of the accretion flow (see in Fig. A2). Figure 2. View largeDownload slide Radial structure of the angle-integrated mass inflow and outflow rates for the case with RC/RB = 0.1 and α = 0.01 at t ≃ 13 tdyn. Each curve presents the inflow rate ($$-\dot{M}_{\rm in}$$ in equation 19; red), the outflow rate ($$\dot{M}_{\rm out}$$ in equation 20; green), and the net accretion rate ($$-\dot{M}\equiv -\dot{M}_{\rm in}-\dot{M}_{\rm out}$$; blue), respectively. Inside the Bondi radius (r ≲ RB), the inflow and outflow rate decrease towards the centre following ∝ r and are tightly balanced. Within the centrifugal radius, where a geometrically thick accretion disc forms, the net accretion rate becomes nearly constant. Figure 2. View largeDownload slide Radial structure of the angle-integrated mass inflow and outflow rates for the case with RC/RB = 0.1 and α = 0.01 at t ≃ 13 tdyn. Each curve presents the inflow rate ($$-\dot{M}_{\rm in}$$ in equation 19; red), the outflow rate ($$\dot{M}_{\rm out}$$ in equation 20; green), and the net accretion rate ($$-\dot{M}\equiv -\dot{M}_{\rm in}-\dot{M}_{\rm out}$$; blue), respectively. Inside the Bondi radius (r ≲ RB), the inflow and outflow rate decrease towards the centre following ∝ r and are tightly balanced. Within the centrifugal radius, where a geometrically thick accretion disc forms, the net accretion rate becomes nearly constant. In Fig. 3, we present the distribution of the gas density (top) and temperature (bottom) for two cases with RC/RB = 0.3 (left) and 0.1 (right). We also plot isodensity contours and the velocity vectors in the top and bottom panel, respectively. The elapsed time is set to t = 13 tdyn as shown in Fig. 2. For two cases, the density distribution around the Bondi radius is approximated by a rotational equilibrium distribution with a constant specific angular momentum, i.e. the Fishbone–Moncrief solution (see equation 17 and Fig. 5 below). Red dashed curve presents the boundary where ρ = 0 for the initial distribution. At the vicinity of the centrifugal radius (r ≲ 2 RC) and inside the zero-density boundary for the initial conditions (red dashed), the density distribution deviates from the equilibrium solution since the angular momentum begins to be transported outwards and the profile approaches j ∝ r1/2 (see Figs 6 and A4 below). The gas temperature increases towards the centre by compressional heating due to the gravity of the BH and energy dissipation due to viscosity. Gas motions in non-azimuthal directions are subsonic and form circulations. We note that some fraction of the gas is outflowing from inside the Bondi radius (thick black curve). In the polar directions, weakly collimated hot outflows are launched (see also green curve in Fig. 2 at r ≳ 2 RB). The outflow speed increases with the angular momentum of the system. Figure 3. View largeDownload slide Distribution of the gas density (top) and temperature (bottom) for different values of angular momentum with RC/RB = 0.3 (left) and 0.1 (right). Isodensity contours and velocity vectors are shown in the top and bottom panel, respectively. The elapsed time is set to t = 13 tdyn as shown in Fig. 2. In the top panel, the zero-density boundary for the initial distribution is shown by red dashed curve, inside which the density distribution deviates from the equilibrium solution. The thick black curve in the bottom panel shows the location of the Bondi radius. Figure 3. View largeDownload slide Distribution of the gas density (top) and temperature (bottom) for different values of angular momentum with RC/RB = 0.3 (left) and 0.1 (right). Isodensity contours and velocity vectors are shown in the top and bottom panel, respectively. The elapsed time is set to t = 13 tdyn as shown in Fig. 2. In the top panel, the zero-density boundary for the initial distribution is shown by red dashed curve, inside which the density distribution deviates from the equilibrium solution. The thick black curve in the bottom panel shows the location of the Bondi radius. Fig. 4 shows the density distribution at smaller scale for RC/RB = 0.3 (left) and 0.1 (right), respectively. For the two cases, a dense torus structure is formed around the central BH, where the gas motions are very subsonic even inside the Bondi radius and highly convective (see Section 3.1.1). From the torus surface with ρ ≃ 2 ρ∞, outflows are launched towards the poles. Except in the polar directions, the gas is both inflowing and outflowing, and the two rates are almost balanced. We emphasize that mirror-symmetry of the accretion flow cross the equatorial plane (θ = π/2) is broken. Although a previous study by Li et al. (2013) imposed equatorial mirror-symmetry and found that a large fraction of the gas is outflowing through the equator coherently, the equatorial outflow is an artefact of the imposed symmetry (see also Roberts et al. 2017). Figure 4. View largeDownload slide Distribution of the gas density in the inner region for different values of angular momentum with RC/RB = 0.3 (left) and 0.1 (right). Isodensity contours and velocity vectors are shown, and the elapsed time is set to t = 13 tdyn as shown in Fig. 2. Figure 4. View largeDownload slide Distribution of the gas density in the inner region for different values of angular momentum with RC/RB = 0.3 (left) and 0.1 (right). Isodensity contours and velocity vectors are shown, and the elapsed time is set to t = 13 tdyn as shown in Fig. 2. Angular momentum affects the gas density near the central hole. Namely, the absolute value is higher with lower angular momentum. This is because concentration of the density is suppressed by convective motions inside the centrifugal radius. As we will discuss in Section 3.1.1, the density profile becomes flatter (ρ ∝ r−1/2) inside the convective envelope (r ≲ 2 RC) rather than ρ ∝ r−3/2. 3.1.1 Comparison to self-similar solutions In order to understand properties of the accretion flows, we consider time-averaged profiles of physical quantities and compare to those of self-similar solutions for radiatively inefficient accretion flows. In the following, we show time-averaged values over 6 ≤ t/tdyn ≤ 16. Fig. 5 shows radial profiles of the gas density along the equator (θ = π/2) for three cases with RC/RB = 0.1 (red), 0.2 (green), and 0.3 (blue). The profiles consist of two components. In the outer region (r > 2 RC), the density profiles are similar to the equilibrium solution of ρ/ρ∞ = (1 + RB/r)3/2 (black dotted). On the other hand, the density follows ρ ∝ r−1/2 in the inner region (r < 2 RC) and is well approximated by   \begin{equation} \rho \simeq \rho _\infty \frac{(1+2\beta )^{3/2}}{2\beta }\left(\frac{r}{R_{\rm B}}\right)^{-1/2}, \end{equation} (21)(see dashed lines in Fig. 5). The slope of the density profile in this region agrees with that for a CDAF solution (ρ ∝ r−1/2; Narayan et al. 2000; Quataert & Gruzinov 2000) rather than that for an ADAF solution (ρ ∝ r−3/2; Narayan & Yi 1994, 1995). Figure 5. View largeDownload slide Radial profiles of the gas density along the equator (θ = π/2) for the cases with RC/RB = 0.1 (red) 0.2 (green) and 0.3 (blue). The density profiles are time-averaged over 6 ≤ t/tdyn ≤ 16, and consist of two components. In the outer region, they are approximated by an equilibrium solution of rotating gas, ‘Fishbone–Moncrief solution’, following ρ/ρ∞ = (1 + RB/r)3/2 (dotted curve). In the inner region, the density profiles approach ρ ∝ r−1/2 (dashed lines; equation 21), whose slope is the same as in CDAF solutions. The transition occurs at r ≃ 2 RC for each case (open circles). Figure 5. View largeDownload slide Radial profiles of the gas density along the equator (θ = π/2) for the cases with RC/RB = 0.1 (red) 0.2 (green) and 0.3 (blue). The density profiles are time-averaged over 6 ≤ t/tdyn ≤ 16, and consist of two components. In the outer region, they are approximated by an equilibrium solution of rotating gas, ‘Fishbone–Moncrief solution’, following ρ/ρ∞ = (1 + RB/r)3/2 (dotted curve). In the inner region, the density profiles approach ρ ∝ r−1/2 (dashed lines; equation 21), whose slope is the same as in CDAF solutions. The transition occurs at r ≃ 2 RC for each case (open circles). Fig. 6 shows radial profiles along the equator of (a) the sound speed, (b) the rotational velocity and (c) the Bernoulli number, which is defined by   \begin{equation} Be \equiv \frac{v^2}{2}+\frac{c_{\rm s}^2}{\gamma -1}-\frac{GM_\bullet }{r}. \end{equation} (22)Those profiles inside the centrifugal radius (r < 2 RC) are summarized as $$c_{\rm s}^2\propto r^{-1}$$, vϕ/vK ≃ 1, and $$Be/v_{\rm K}^2\simeq 0$$, where $$v_{\rm K}=\sqrt{GM_\bullet /r}$$ is the Keplerian velocity. Since gas accretion begins from the Bondi radius, the Bernoulli number is close to zero unless energy dissipation via viscosity is significant. Within the centrifugal radius (r ≲ RC), the rotational velocity dominates and is approximated as ∼vK. Therefore, the sound speed follows as $$c_{\rm s}^2\propto r^{-1}$$ as shown in panel (a). We note that the profiles of the sound speed do not depend on the angular momentum of the flow. In fact, the profiles can be fit with   \begin{equation} c_{\rm s}^2=c_{\infty }^2 \left[ 1+f(\gamma )\frac{R_{\rm B}}{r} \right], \end{equation} (23)where we adopt a functional form of f(γ) as   \begin{equation} f(\gamma )\equiv \frac{\gamma -1}{\gamma +1-(\gamma -1)/2} \end{equation} (24)Moreover, the rotational velocity inside the centrifugal radius is approximated as   \begin{equation} \frac{v_\phi }{v_{\rm K}} = g(\gamma ) \equiv \sqrt{\frac{2-(\gamma -1)}{\gamma +1-(\gamma -1)/2}}. \end{equation} (25)For γ = 1.6, f(γ) = 0.26 and g(γ) = 0.78. These functions of f and g have been calculated by Quataert & Gruzinov (2000) for marginally stable rotating-accretion flows against convection motions, i.e. CDAF solutions. Figure 6. View largeDownload slide Radial profile of the time-averaged (a) sound speed, (b) rotational velocity, and (c) Bernoulli number along the equator (θ = π/2) for three different cases with RC/RB = 0.1 (red), 0.2 (green), and 0.3 (blue). The dotted curves in the top two panels (equations 23 and 25) are those of the CDAF solutions, respectively. Figure 6. View largeDownload slide Radial profile of the time-averaged (a) sound speed, (b) rotational velocity, and (c) Bernoulli number along the equator (θ = π/2) for three different cases with RC/RB = 0.1 (red), 0.2 (green), and 0.3 (blue). The dotted curves in the top two panels (equations 23 and 25) are those of the CDAF solutions, respectively. Fig. 7 shows profiles of the density as a function of the polar angle θ, at radial positions of 0.02 ≤ r/RB ≤ 1.0. In the outer region (r ≳ 2 RC), the density profile follows an equilibrium distribution and explained by the Fishbone–Moncrief solution (short-dashed). In the inner region (r < 2 RC = 0.2 RB), the density concentration to the mid-plane is higher. The angular profiles in the inner region can be explained by long-dashed curve, ρ(θ) ∝ (sin θ)2[1/(γ − 1) − 1/2][ = (sin θ)2.33 for γ = 1.6] (Quataert & Gruzinov 2000). This result also indicates that the profiles of the accretion flow within the centrifugal radius are self-similar, and those properties approach those in CDAFs. Fig. 8 also presents angular profiles of the Bernoulli number. In the inner region (r < 2 RC), the value of Be is smaller than $${\lesssim } 0.2\ v_{\rm Kep}^2$$ except in the vicinity of the poles. We note that the value increases rapidly towards the poles, where weak outflows are launched. Figure 7. View largeDownload slide Angular profiles of the time-averaged density at radial positions of 0.02 ≤ r/RB ≤ 1.0. The results in the inner region (r < 2 RC) and outer region (r ≥ 2 RC) are shown by solid and dotted curves, respectively. Analytical angular profiles for CDAF solutions ρ(θ) ∝ (sin θ)[2/(γ − 1) − 1] (long-dashed) and for the Fishbone–Moncrief solution at r = RB (short-dashed) are shown for comparison to the numerical results. The profiles are normalized by the maximum value for each case. Figure 7. View largeDownload slide Angular profiles of the time-averaged density at radial positions of 0.02 ≤ r/RB ≤ 1.0. The results in the inner region (r < 2 RC) and outer region (r ≥ 2 RC) are shown by solid and dotted curves, respectively. Analytical angular profiles for CDAF solutions ρ(θ) ∝ (sin θ)[2/(γ − 1) − 1] (long-dashed) and for the Fishbone–Moncrief solution at r = RB (short-dashed) are shown for comparison to the numerical results. The profiles are normalized by the maximum value for each case. Figure 8. View largeDownload slide Angular profiles of the time-averaged Bernoulli number at radial positions of 0.02 ≤ r/RB ≤ 1.0. In the inner region (r < 2 RC), the value of Be is smaller than $${\sim } 0.2\ v_{\rm Kep}^2$$. Note that the value increases rapidly towards the poles, where weak outflows are launched. Figure 8. View largeDownload slide Angular profiles of the time-averaged Bernoulli number at radial positions of 0.02 ≤ r/RB ≤ 1.0. In the inner region (r < 2 RC), the value of Be is smaller than $${\sim } 0.2\ v_{\rm Kep}^2$$. Note that the value increases rapidly towards the poles, where weak outflows are launched. 3.1.2 Angular momentum and energy transport via convection motions Angular momentum transport by convection has been discussed with analytical methods and numerical simulations by previous works (e.g. Igumenshchev & Abramowicz 2000; Igumenshchev, Abramowicz & Narayan 2000; Narayan et al. 2000; Quataert & Gruzinov 2000; Igumenshchev et al. 2003). To analyse the effect, we calculate the r–ϕ component of Reynolds stress, $$\tau _{r \phi }\equiv -\rho \langle v^{\prime }_r v^{\prime }_\phi \rangle$$, where $$v^{\prime }_i\equiv v_i-\langle v_i \rangle$$ and 〈 · 〉 means the time-averaged value over 6 ≤ t/tdyn ≤ 16. In practice, the mass-weighted Reynolds stress integrated over the polar angle is calculated as   \begin{equation} \tau _{r\phi }(r) = -\int ^{\pi}_0 \langle \rho v^{\prime }_r v^{\prime }_\phi \rangle \sin \theta\, {\rm d} \theta . \end{equation} (26)In Fig. 9, we show the radial profile of the Reynolds stress normalized by $$\rho v_{\rm K}^2$$ for three cases of RC/RB = 0.1 (red), 0.2 (green), and 0.3 (blue). The Reynolds stress does not depend on the initial angular momentum, and the value is approximated as τrϕ ≃ Ar−3/2, where A is positive. Note that the r–ϕ component of the viscous stress has a negative value (i.e. σrϕr3/2 ≃ const. < 0). Therefore, the positive sign of A means that convective motions transport angular momentum inwards, while the standard α-viscosity transports it outwards. Igumenshchev et al. (2000) have conducted three-dimensional simulations relaxing the assumption of axisymmetry and found that the results are essentially similar to those obtained in two-dimensional simulations. In fact, the convective eddies are nearly axisymmetric and transport angular momentum inwards. This allows us to justify that our two-dimensional simulations can capture the important physics. Figure 9. View largeDownload slide Radial profiles of the r − ϕ component of the Reynolds stress τrϕ for different three cases with RC/RB = 0.1 (red), 0.2 (green), and 0.3 (blue). Solid and dashed curves show positive and negative values, respectively. The values are normalized by $$\rho v_{\rm K}^2$$. Since τrϕ ∝ Ar−3/2 is positive, convection motions transport angular momentum inwards within the centrifugal radius. Figure 9. View largeDownload slide Radial profiles of the r − ϕ component of the Reynolds stress τrϕ for different three cases with RC/RB = 0.1 (red), 0.2 (green), and 0.3 (blue). Solid and dashed curves show positive and negative values, respectively. The values are normalized by $$\rho v_{\rm K}^2$$. Since τrϕ ∝ Ar−3/2 is positive, convection motions transport angular momentum inwards within the centrifugal radius. Fig. 10 shows radial profiles of three kinds of torques for the case with RC/RB = 0.1; the Reynolds stress (red), the viscosity (green) and the advection (blue). Those three terms are calculated as   \begin{equation} T_{\rm Rey} \equiv - \frac{\mathrm{\partial} }{\mathrm{\partial} r } \left\langle 2\pi r^2 \Sigma v^{\prime }_r v^{\prime }_\phi \right\rangle , \end{equation} (27)  \begin{equation} T_{\rm vis} \equiv \frac{\mathrm{\partial} }{\mathrm{\partial} r} \left\langle 2\pi r^3 \nu \Sigma \frac{\mathrm{\partial} (v_\phi /r)}{\mathrm{\partial} r}\right\rangle , \end{equation} (28)  \begin{equation} T_{\rm adv} \equiv -\frac{\mathrm{\partial} }{\mathrm{\partial} r} \left[ \left\langle 2\pi r^2 \Sigma v_r \right\rangle \cdot \left\langle v_\phi \right\rangle \right], \end{equation} (29)where the surface density is given by $$\Sigma \approx \int _0^{\pi} \rho r\sin ^2 \theta {\rm d}\theta$$. This approximation is valid when the density concentrates near the equator, i.e. ρ(θ) ∝ (sin θ)2.33. The sign of the torque suggests that convection motions transport angular momentum inwards (red solid) rather than outwards (red dotted), as expected from Fig. 9. Within the centrifugal radius (r ≲ RC), the angular momentum is transported inwards via advection (Tadv > 0) and outwards via viscosity (Tvis < 0). Both of them are tightly balanced, and the residual agrees with the Reynolds torque within a factor of 2. Figure 10. View largeDownload slide Radial profiles of the time-averaged torque (over 6 ≤ t/tdyn ≤ 16) for the case with RC/RB = 0.1. Each curve shows the Reynolds torque (red), the viscous torque (green), and the advection torque (blue). The Reynolds component due to convection transports the angular momentum inwards (solid) and outwards (dotted), respectively. Within the centrifugal radius (r ≲ RC), the angular momentum transports via advection (inwards) and viscosity (outwards) are almost balanced. Figure 10. View largeDownload slide Radial profiles of the time-averaged torque (over 6 ≤ t/tdyn ≤ 16) for the case with RC/RB = 0.1. Each curve shows the Reynolds torque (red), the viscous torque (green), and the advection torque (blue). The Reynolds component due to convection transports the angular momentum inwards (solid) and outwards (dotted), respectively. Within the centrifugal radius (r ≲ RC), the angular momentum transports via advection (inwards) and viscosity (outwards) are almost balanced. Convective motions can transport thermal energy as well. In order to evaluate the convective energy flux Fconv, we consider the time-averaged energy equation as2  \begin{equation} \frac{1}{r^2}\frac{\mathrm{\partial} }{\mathrm{\partial} r}(r^2F_{\rm conv}) = \left\langle \sigma _{r\phi }r\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left(\frac{v_\phi }{r}\right) \right\rangle - \left\langle \rho Tv_r\frac{{\rm d}s}{{\rm d}r} \right\rangle . \end{equation} (30)The two terms on the right-hand side present time-averaged values of the energy generation rate by viscous dissipation ($$Q^+_{\rm vis}$$) and the energy advection rate ($$Q^-_{\rm adv}$$), respectively. Contribution due to the convective energy flux can be approximated by the difference between $$\langle Q^+_{\rm vis}\rangle$$ and $$\langle Q^-_{\rm adv}\rangle$$. Integrating the energy equation over angles except the polar regions, we obtain   \begin{equation} \frac{{\rm d}L_{\rm conv}}{{\rm d}r}(r)=2\pi r^2 \int _{D_\theta } (\langle Q^+_{\rm vis}\rangle -\langle Q^-_{\rm adv}\rangle )\sin \theta\, {\rm d}\theta , \end{equation} (31)where the convective luminosity is defined as $$L_{\rm conv}=2\pi r^2 \int _{D_\theta }F_{\rm conv} \sin \theta\, {\rm d}\theta$$ and Dθ = [π/6, 5π/6]. Fig. 11 shows the radial profiles of dLconv/d ln r for different values of RC/RB (dashed). Since the derivative is positive at 0.02 ≲ r/RB ≲ 1, thermal energy is transported outwards by convection. Those values are almost constant, namely $${\rm d}L_{\rm conv}/{\rm d}\ln r\simeq (0.02{\rm -}0.04)\times \dot{M}_{\rm B}c_\infty ^2$$. We note that the energy advection dominates at the innermost region (r < 0.02 RB). This is due to our inner boundary conditions, where we do not consider energy output from the sink cell at the centre. In order to study the effect of the boundary condition, we conduct a simulation for RC/RB = 0.1 with a smaller rmin, which is set to 3 × 10−3 RB (red solid). This result clearly shows that a smaller innermost grid does not change the absolute value of dLconv/d ln r but makes the convection-dominated region larger. We integrate the value from 10−2 RB to RB and calculate the convective luminosity (black solid). Then, the convective luminosity is written as $$L_{\rm conv}\sim \ \eta _{\rm conv}\dot{M}_{\rm B}c_\infty ^2$$, where the efficiency is approximated as ηconv ≃ 0.2. Figure 11. View largeDownload slide Radial profiles of the convective luminosity, dLconv/d ln r (dashed) for different three cases with RC/RB = 0.1 (red), 0.2 (green), and 0.3 (blue). In order to study the effect of the boundary condition, we conduct a simulation for RC/RB = 0.1 with a smaller rmin = 3 × 10−3RB (red solid). The integrated convective luminosity is shown by black solid, presenting $$L_{\rm conv}\simeq 0.2\ \dot{M}_{\rm B}c_\infty ^2$$. Figure 11. View largeDownload slide Radial profiles of the convective luminosity, dLconv/d ln r (dashed) for different three cases with RC/RB = 0.1 (red), 0.2 (green), and 0.3 (blue). In order to study the effect of the boundary condition, we conduct a simulation for RC/RB = 0.1 with a smaller rmin = 3 × 10−3RB (red solid). The integrated convective luminosity is shown by black solid, presenting $$L_{\rm conv}\simeq 0.2\ \dot{M}_{\rm B}c_\infty ^2$$. 3.2 Effects of viscous parameter Next, we discuss the dependence of our results on the viscous parameter α, conducting several simulations with the same parameters as shown in Section 3.1 except varying α. Fig. 12 shows the time evolution of the accretion rate on to a BH at the centre for different values of α (0.003 ≤ α ≤ 0.1). To keep the numerical simulations stable, the viscous parameter is assumed to be 0.01 at t ≤ 2 tdyn, increase (or decrease) linearly proportional to the time until t = 4 tdyn, and keep a constant value what we consider at t > 4 tdyn. The accretion rates increase with the viscous parameter because angular momentum transport becomes more efficient. For the lowest value of α, inflows with lower angular momentum through the polar region is not negligible because inflows driven by viscosity through the torus become less efficient. For the highest value of α, the fluctuation in the accretion rate is suppressed by strong viscosity, which is consistent with that reported by previous numerical simulations by Igumenshchev & Abramowicz (1999, 2000). Figure 12. View largeDownload slide Same as Fig. 1, but for different viscous parameters: α = 0.003 (magenta), 0.01 (red), 0.03 (green), and 0.1 (blue). The angular momentum is set so that RC/RB = 0.1. The accretion rate increases with α because angular momentum transport is more efficient. Figure 12. View largeDownload slide Same as Fig. 1, but for different viscous parameters: α = 0.003 (magenta), 0.01 (red), 0.03 (green), and 0.1 (blue). The angular momentum is set so that RC/RB = 0.1. The accretion rate increases with α because angular momentum transport is more efficient. In Fig. 13, we show radial profiles of the gas density for three different values of the viscous parameter: α = 0.01 (red), 0.03 (green), and 0.1 (blue). The profiles are approximated by an equilibrium solution following ρ/ρ∞ = (1 + RB/r)3/2 (dotted curve) in the outer region and ρ ∝ r−1/2, which is consistent with that in CDAF solutions in the inner region (r < 0.2 RB). For all the cases, the overall behaviour is similar, i.e. the density profile is not dependent on the choice of viscous parameter. However, for the highest value of α( = 0.1), the density slope becomes as steep as −3/2 at the innermost region (r < 0.02 RB). In this case, convective motions are suppressed by strong viscosity and thus the energy generated by viscous heating is not transported by convection but by advection with inflows, i.e. ADAFs (Narayan & Yi 1994, 1995). Figure 13. View largeDownload slide Same as Fig. 5 but for different viscous parameters: α = 0.01 (red), 0.03 (green), and 0.1 (blue). The angular momentum is set to RC/RB = 0.1. The profiles are approximated by an equilibrium solution of ρ/ρ∞ = (1 + RB/r)3/2 (black dotted) in the outer region and by ρ ∝ r−1/2 in the inner region (r < 0.2 RB), which is consistent with that in CDAF solutions, For the highest α, the density profile approaches ρ ∝ r−3/2 in the innermost region (see the text for explanation). Figure 13. View largeDownload slide Same as Fig. 5 but for different viscous parameters: α = 0.01 (red), 0.03 (green), and 0.1 (blue). The angular momentum is set to RC/RB = 0.1. The profiles are approximated by an equilibrium solution of ρ/ρ∞ = (1 + RB/r)3/2 (black dotted) in the outer region and by ρ ∝ r−1/2 in the inner region (r < 0.2 RB), which is consistent with that in CDAF solutions, For the highest α, the density profile approaches ρ ∝ r−3/2 in the innermost region (see the text for explanation). Radial profiles of other physical quantities, e.g. cs, vϕ, and Be, do not change significantly from those with α = 0.01. However, for the highest value of α( = 0.1), the normalized rotational velocity decreases at r < 0.03 RB to the centre to vϕ/vK ≃ 0.3. Such a small ratio is found in ADAF solutions (Narayan & Yi 1995). Fig. 14 shows angular profiles of the gas density at r = 0.02 RB for different viscous parameters: α = 0.01 (red), 0.03 (green), and 0.1 (blue). For lower values of α( < 0.1), the angular profiles of the density follow that of a CDAF solution (black dashed). On the other hand, for the highest α( = 0.1), where the radial density profile is as steep as −3/2, the angular profile is also different from the CDAF one. Such non-equatorial symmetric profiles produced by outflows towards one of the poles have been reported in previous numerical simulations by Igumenshchev & Abramowicz (2000) (their Model G, where α = 0.1). Figure 14. View largeDownload slide Same as Fig. 7 at r = 0.02 RB but for different viscous parameters: α = 0.01 (red), 0.03 (green), and 0.1 (blue). The angular momentum is set to RC/RB = 0.1. Dashed curve presents an angular profile expected in CDAF solutions: ρ(θ) ∝ (sin θ)[2/(γ − 1) − 1]. Figure 14. View largeDownload slide Same as Fig. 7 at r = 0.02 RB but for different viscous parameters: α = 0.01 (red), 0.03 (green), and 0.1 (blue). The angular momentum is set to RC/RB = 0.1. Dashed curve presents an angular profile expected in CDAF solutions: ρ(θ) ∝ (sin θ)[2/(γ − 1) − 1]. Fig. 15 presents the dependence of the net accretion rate on the viscous parameter α. Each symbol shows the total net accretion rates (blue cross) and the rate through the equatorial region of π/6 ≤ θ ≤ 5π/6 (red square). The latter can clarify the effect of angular momentum transport via viscosity in the torus. Those results for α ≤ 0.2 can be fit by $$\dot{M}/\dot{M}_{\rm B}\simeq 0.12\ \alpha ^\delta$$, where δ = 0.62 ± 0.24 (dashed). We note that the relation is qualitatively different from that considered in previous works; self-similar solutions or spherical accretion solutions without treating convection ($$\dot{M}\simeq \alpha \dot{M}_{\rm B}$$; Narayan & Fabian 2011). Therefore, this relation between $$\dot{M}$$ and α is an essential result from the self-consistent solution, connecting the Fishbone–Moncrief quasi-static solution (large scales) to existing CDAF solutions (small scales), and show the importance of the global solution because without our simulates, it was unclear how physical quantities on larger scales determine the properties of the accretion flow on small scales. Figure 15. View largeDownload slide Dependence of the net gas inflow rate on the viscous parameter α. The net rate increases with the viscous parameter because angular momentum is transported more efficiently. The total net accretion rates (blue cross) and the inflow rate through the equatorial region (red square) are shown. The results for α ≤ 0.2 can be fit by $$\dot{M}/\dot{M}_{\rm B}\simeq 0.12\ \alpha ^{0.62}$$ (dashed). Figure 15. View largeDownload slide Dependence of the net gas inflow rate on the viscous parameter α. The net rate increases with the viscous parameter because angular momentum is transported more efficiently. The total net accretion rates (blue cross) and the inflow rate through the equatorial region (red square) are shown. The results for α ≤ 0.2 can be fit by $$\dot{M}/\dot{M}_{\rm B}\simeq 0.12\ \alpha ^{0.62}$$ (dashed). 4 GLOBAL SOLUTIONS OF RADIATIVELY INEFFICIENT ACCRETION FLOWS We now briefly discuss the global solution of radiatively inefficient rotating-accretion flows, extending our results down to the central BH. We also give an analytical expression for the net BH feeding rate in the global solution, which can be compared to observations of low-luminosity BHs such as Sgr A* and the BH in M87. 4.1 BH feeding rate and energy loss via radiation As shown in Fig. 2, the mass inflow rate decreases towards the centre as $$\dot{M}\simeq \rho |v_r| r^2 \propto r$$. This is because the density profile of a CDAF solution follows ρ ∝ r−1/2 and the radial velocity follows |vr| ∼ ν/r ∝ r−1/2 in the α-viscosity model, respectively. Then, the net accretion rate (including both inflows and outflows) is independent of radius and at a much lower rate than the Bondi value. Because of limitation of computational time, we do not extend our computational domain down to the central BH (r ∼ RSch). Instead, we perform two additional simulations with different locations of the innermost grid of rmin = 3 × 10−3 RB (red) and 1.3 × 10−2 RB (blue). Fig. 16 shows radial profiles of the angle-integrated inflow rate (dashed) and the net accretion rate (solid). As shown clearly, the net accretion rate decreases with rmin. Therefore, the time-averaged value of the net accretion rate can be approximated as   \begin{equation} \dot{M}\simeq \left(\frac{\alpha }{0.01}\right)^\delta \left(\frac{r_{\rm min}}{R_{\rm B}}\right) \dot{M}_{\rm B}, \end{equation} (32)where δ ≃ 0.62. Assuming that the net accretion rate is constant at r ≲ 2 RC, we can estimate the net radial velocity as $$\overline{v_r} \equiv -\dot{M}/4\pi \rho r^2$$, or for β( = RC/RB) ≪ 1   \begin{equation} \overline{v_r} \simeq - \frac{\beta }{2} \left(\frac{\alpha }{0.01}\right)^\delta \left(\frac{r_{\rm min}}{R_{\rm B}}\right) \left(\frac{r}{R_{\rm B}}\right)^{-3/2} \ c_\infty . \end{equation} (33)Using the net inflow velocity, the dynamical time-scale is estimated as $$t_{\rm dyn}=r/\overline{v_r}\propto r^{5/2}$$. Inside the torus, heating via viscous dissipation dominates and the heating time-scale is given by tvis ≃ 1/[γ(γ − 1)αΩ]−1 ∝ r3/2. The ratio of the two time-scales is given by   \begin{equation} \frac{t_{\rm vis}}{t_{\rm dyn}}\sim 0.6 \left(\frac{\alpha }{0.01}\right)^{\delta -1} \left(\frac{r_{\rm min}}{10^{-2}\ R_{\rm B}}\right) \left(\frac{r}{R_{\rm C}}\right)^{-1}. \end{equation} (34)Thus, we find tvis ≲ tdyn except at r ≲ 6(β/0.1) rmin. Figure 16. View largeDownload slide Radial structure of the angle-integrated mass inflow and outflow rates for different sizes of the innermost grid rmin. Physical parameters are set to RC/RB = 0.1 and α = 0.01, and the elapsed time is t ≃ 13 tdyn. Solid and dashed curves present the inflow rate and the net accretion rate (including both inflows and outflows), respectively. Since the density follows ρ ∝ r−1/2 and the radial velocity follows |vr| ∼ ν/r ∝ r−1/2, the inflow rate is proportional to the radius ($$\dot{M}\propto r$$) within the centrifugal radius. The net accretion rate decreases with rmin and can be approximated by equation (32). The dotted line shows $$\dot{M}=(r_{\rm min}/R_{\rm B})\dot{M}_{\rm B}$$. Figure 16. View largeDownload slide Radial structure of the angle-integrated mass inflow and outflow rates for different sizes of the innermost grid rmin. Physical parameters are set to RC/RB = 0.1 and α = 0.01, and the elapsed time is t ≃ 13 tdyn. Solid and dashed curves present the inflow rate and the net accretion rate (including both inflows and outflows), respectively. Since the density follows ρ ∝ r−1/2 and the radial velocity follows |vr| ∼ ν/r ∝ r−1/2, the inflow rate is proportional to the radius ($$\dot{M}\propto r$$) within the centrifugal radius. The net accretion rate decreases with rmin and can be approximated by equation (32). The dotted line shows $$\dot{M}=(r_{\rm min}/R_{\rm B})\dot{M}_{\rm B}$$. Next, in order to evaluate the effect of radiative cooling, we compare the heating time-scale to the cooling time-scale. Since ρ ∝ r−1/2 and T ∝ r−1 inside the centrifugal radius, the bremsstrahlung cooling rate per volume is $$Q^-_{\rm br} \propto \rho ^2 T^{1/2} \propto r^{-3/2}$$ and the cooling time-scale is $$t_{\rm cool}\propto \rho T/Q^-_{\rm br} \propto r^0$$. Since tvis ≲ tdyn at r ≃ RC, thus we estimate the ratio of tvis to tcool as   \begin{equation} \frac{t_{\rm vis}}{t_{\rm cool}}\sim 0.25 \left(\frac{\alpha }{0.01}\right)^{-1} \left(\frac{r}{R_{\rm C}}\right)^{3/2} \left(\frac{T_\infty }{10^7\ {\rm K}}\right)^{1/2} \left(\frac{\dot{m}_{\rm B}}{10^{-3}}\right), \end{equation} (35)where γ = 1.6 and β = 0.1 are assumed. Bremsstrahlung cooling does not play an important role in the accretion flow as long as $$\dot{m}_{\rm B}<4\times 10^{-3}$$, which is consistent with the results shown in Li et al. (2013). Since $$Q^-_{\rm br} \propto r^{-3/2}$$, the CDAF is sandwiched between an inner region where gravitational energy would be released close to the BH and an outer radiating region (Narayan et al. 2000; Quataert & Gruzinov 2000; Ball, Narayan & Quataert 2001). 4.2 Thermal conductivity In a hot accretion flow, thermal conductivity potentially transports energy outwards and could affect the gas dynamics (e.g. Johnson & Quataert 2007; Sharma, Quataert & Stone 2008; Shcherbakov & Baganoff 2010). In the classic picture, the conductive energy flux is estimated as   \begin{equation} \boldsymbol F_{\rm cond}=-\kappa \boldsymbol \nabla T, \end{equation} (36)where κ is the conduction coefficient given by Spitzer (1962) as   \begin{equation} \kappa = 5.0\times 10^{-7} \left(\frac{\ln \Lambda _{\rm c}}{37}\right)^{-1}T^{5/2} f_{\rm c}, \end{equation} (37)in units of erg s−1 cm−1 K−1, and fc is the conductivity suppression factor because thermal conduction in the perpendicular directions to magnetic fields can be suppressed. The value of the suppression factor has been discussed by various theoretical arguments and estimated as fc ∼ 0.1 (e.g. Narayan & Medvedev 2001; Maron, Chandran & Blackman 2004). However, the precise value of the suppression factor is highly uncertain and depends on the nature of magnetic turbulence. It is known that there are two instabilities of magnetized plasma: the magneto-thermal instability (MTI; Balbus 2000) and the heat-flux-driven buoyancy instability (HBI; Quataert 2008). The MTI occurs when the temperature decreases with height and the HBI does in the opposite case. McCourt et al. (2011) have performed long-term and global MHD simulations in order to follow the evolution of both instabilities into the non-linear regime. They have found the following two results: (1) the MTI can drive strong turbulence leading to the isotropic configuration of magnetic fields (fc ∼ 0.3) and produce a large convective energy flux of $${\gtrsim}0.01 \rho c_{\rm s}^3$$ and (2) the HBI reorients the magnetic field and suppresses the conductive heat flux through the plasma. Since the value of the suppression factor is uncertain, we adopt fc = 0.1 as a fiducial one. Using our accretion solution inside the centrifugal radius, where the temperature follows T ≃ T∞f(γ)RB/r, the conductive luminosity (Lcond ≡ 4πr2Fcond) is estimated as   \begin{equation} L_{\rm cond} \simeq 6.4\times 10^{29}\ M_6 c_7^{-2} \left(\frac{f_{\rm c}}{0.1}\right) \left(\frac{r}{R_{\rm B}}\right)^{-5/2} \ {\rm erg\ s^{-1}}. \end{equation} (38)On the other hand, the convective luminosity is roughly given by   \begin{eqnarray} L_{\rm conv} &\simeq & \eta _{\rm conv} \left(\frac{\alpha }{0.01}\right)^\delta \dot{M}_{\rm B} c_\infty ^2, \nonumber \\ &\simeq & 1.1\times 10^{35}\ \rho _{-22}M_6^2 c_7^{-1} \left(\frac{\eta _{\rm conv}}{0.2}\right) \left(\frac{\alpha }{0.01}\right)^\delta \ {\rm erg\ s^{-1}}. \end{eqnarray} (39)Thus, the ratio of the two luminosities is   \begin{eqnarray} \frac{L_{\rm cond}}{L_{\rm conv}} &\simeq & 4.5\times 10^{-8} \left(\frac{f_{\rm c}}{0.1}\right) \left(\frac{\eta _{\rm conv}}{0.2}\right)^{-1} \left(\frac{\alpha }{0.01}\right)^{-\delta } \nonumber \\ &&\ \times \left(\frac{\dot{m}_{\rm B}}{10^{-3}}\right)^{-1} \left(\frac{T_\infty }{10^7\ {\rm K}}\right)^{-2} \left(\frac{r}{R_{\rm B}}\right)^{-5/2}. \end{eqnarray} (40)Therefore, we can estimate a characteristic radius where Lconv ≃ Lcond as   \begin{eqnarray} \frac{R_{\rm tr}}{R_{\rm B}} & \simeq & 1.5\times 10^{-3} \left(\frac{f_{\rm c}}{\eta _{\rm conv}}\right)^{2/5} \left(\frac{\alpha }{0.01}\right)^{-2\delta /5} \nonumber \\ && \ \times \left(\frac{\dot{m}_{\rm B}}{10^{-3}}\right)^{-2/5} \left(\frac{T_\infty }{10^7\ {\rm K}}\right)^{-4/5}. \end{eqnarray} (41)Inside this radius, thermal conduction would dominate the energy transport over convection, and the temperature profile would be flatter rather than the adiabatic scaling law (T ∝ r−1). The ratio of cs/vK would decline and the flow would resemble an optically thin viscous disc. Note that thermal conduction would hardly affect the dynamics of accretion flows at the vicinity of the BH (r < 100 RSch), where the conductive energy flux is limited to a few per cent of the saturated value of $${\sim } \rho c_{\rm s}^3$$ (Foucart et al. 2016, 2017). This fraction is consistent with that obtained in our solution at r ≳ Rtr, namely $$L_{\rm cond}/(4\pi r^2 \rho c_{\rm s}^3)\sim 1.5\times 10^{-2}(r/R_{\rm tr})^{-5/2}$$. To explore the nature of thermal conduction on the convective motions is left for future investigations. 4.3 Connection to the innermost region As discussed in Section 4.2, some physical processes, e.g. thermal conduction, can transport the energy outwards instead of convection. Once additional energy transport operates, if any, the temperature profile is no longer that for an adiabatic flow. We characterize the transition radius as Rtr = ξRB (ξ ≪ 1). The critical radius where Lconv ≃ Lcond corresponds to ξ = ξcond (≃10−3). Since the temperature follows T ∝ r−p (0 ≲ p < 1) within Rtr, the accretion disc becomes thinner and thus the density follows ρ ∝ r−3 + 3p/2, assuming that the accretion rate is constant through the disc as $$\dot{M}\simeq \xi (\alpha /0.01)^\delta \dot{M}_{\rm B}$$ (see equation 32). Inside the transition radius, the optical depth to electron scattering increases and is estimated as   \begin{eqnarray} \tau _{\rm es}(r)&=& \int _r^{R_{\rm tr}} \rho \kappa _{\rm es}\,{\rm d}r, \nonumber \\ &\simeq & \frac{40\xi ^{1/2}}{4-3p} \frac{\dot{m}_{\rm B}}{\beta } \left(\frac{c_\infty }{c}\right) \left(\frac{r}{R_{\rm tr}}\right)^{-2+3p/2}, \end{eqnarray} (42)where β ≪ 1 is assumed. We evaluate the optical depth at the innermost stable circular orbit (ISCO; r ∼ 3 RSch) for p = 0 as   \begin{eqnarray} \tau _{\rm es}^{\rm ISCO} \simeq 0.022 \left(\frac{\beta }{0.1}\right)^{-1} \left(\frac{\xi }{\xi _{\rm cond}}\right)^{5/2} \left(\frac{\dot{m}_{\rm B}}{10^{-3}}\right) \left(\frac{T_\infty }{10^7\ {\rm K}}\right)^{3/2}. \end{eqnarray} (43)Since the accretion flow at the ISCO is optically thin even in an isothermal case (p = 0), the gas is more likely to be optically thin for 0 < p < 1. In this paper, our numerical simulations do not treat thermal conductivity because the effect is subdominant in the computation domain of r ≥ rmin( > Rtr). We plan to numerically study the transition to the innermost solution in future work. 4.4 Radiation luminosities of BHs with very low accretion rates Finally, we estimate radiation luminosities produced from a low-density, radiatively inefficient accretion flow. Although our simulations do not treat the inner region (r ≲ Rtr), where energy transport by thermal conduction would dominate that by convection, we can infer the accretion rate on to the central BH from equations (32) and (41),   \begin{eqnarray} \frac{\dot{M}}{\dot{M}_{\rm Edd}}\simeq 1.5\times 10^{-6} \left(\frac{\alpha }{0.01}\right)^{3\delta /5} \left(\frac{T_\infty }{10^7\ {\rm K}}\right)^{-4/5} \left(\frac{\dot{m}_{\rm B}}{10^{-3}}\right)^{3/5}, \end{eqnarray} (44)where fc/ηconv = 1 is set. In order to estimate the radiation luminosity, we adopt a radiative-efficiency model for different values of $$\dot{M}/\dot{M}_{\rm Edd}$$, provided by axisymmetric numerical simulations of accretion flows around a supermassive BH with M• = 108 M⊙ on a small scale (r ≤ 100 RSch), including general relativistic MHD and frequency-dependent radiation transport (Ryan et al. 2017) (see also Ohsuga et al. 2009; Sa̧dowski et al. 2017). We note that their simulation results including radiation processes and electron thermodynamics, are no longer scale-free, i.e. the radiative efficiency we adopt probably depend on the choice of the BH mass. Nevertheless, it is worth demonstrating the radiation luminosity produced from a self-consistent, radiatively inefficient accretion flow, assuming that the radiative efficiency hardly depends on the BH mass. Fig. 17 shows the radiation luminosities estimated from our results (red square) as a function of the Bondi rate normalized by the Eddington luminosity. The radiation luminosity produced in the vicinity of the BH is extremely low because (1) the radiative efficiency decreases with the accretion rate in such low-accretion-rate regime ($$L\propto \dot{M}^{1.7}$$; see fig. 1 in Ryan et al. 2017) and (2) the net accretion rate at r ≲ RC is strongly suppressed by convective motions from the Bondi accretion rate $$\dot{M}_{\rm B}$$. We can approximate the radiation luminosity for $$\dot{M}/\dot{M}_{\rm Edd}\lesssim 10^{-3}$$ as   \begin{equation} \frac{L_{\rm bol}}{L_{\rm Edd}} \simeq 3\times 10^{-5} \left(\frac{\alpha }{0.01}\right)^\delta \frac{\dot{M}_{\rm B}}{\dot{M}_{\rm Edd}}, \end{equation} (45)or   \begin{equation} L_{\rm bol} \simeq 3\times 10^{-6} \dot{M}_{\rm B}c^2 \left(\frac{\alpha }{0.01}\right)^\delta . \end{equation} (46)We also present observational results of SgrA* and BHs in M31 and M87 in (blue symbols in Fig. 17). Those supermassive BHs in the nearby Universe are known as accreting BHs at low rates of $$\dot{M}_{\rm B}/\dot{M}_{\rm Edd}\sim 10^{-5}{\rm -}10^{-3}$$. Moreover, gas properties in the nuclei are studied since the angular sizes of their Bondi radii ( ≳ 1 arcsec) can be resolved well by Chandra observations (e.g. Baganoff et al. 2003; Garcia et al. 2010; Russell et al. 2015, and references therein).3 Our theoretical estimate agrees well with the observational results for SgrA* and M31. Although the observed luminosity of M87 is several times higher than our estimate, this would be because suppression of the accretion by convective motions in the CDAF regime becomes less efficient for $$\dot{M}_{\rm B}/\dot{M}_{\rm Edd}\gtrsim 10^{-3}$$, and/or because there is a level of contamination from the dominant emission due to the jet. Figure 17. View largeDownload slide Bolometric radiation luminosity produced by radiatively inefficient accretion flows (red square). The luminosities are calculated from our results given by equations (32) and (41), combining with a radiative efficiency provided by MHD simulations on a small scale (r ≤ 100 RSch), including general relativity and radiation transfer (Ryan et al. 2017). We also present observational results for SgrA* and BHs in M31 and M87 (blue asterisk). Figure 17. View largeDownload slide Bolometric radiation luminosity produced by radiatively inefficient accretion flows (red square). The luminosities are calculated from our results given by equations (32) and (41), combining with a radiative efficiency provided by MHD simulations on a small scale (r ≤ 100 RSch), including general relativity and radiation transfer (Ryan et al. 2017). We also present observational results for SgrA* and BHs in M31 and M87 (blue asterisk). Isolated BHs in our Galaxy can also accrete interstellar medium (ISM) at very low rates. Assuming typical density of the ISM and molecular clouds (n ∼ 1–100 cm−3) and the BH has a peculiar velocity of V ∼ 40–100 km s−1, the Bondi–Hoyle–Lyttleton accretion rate is as low as $$\dot{M}_{\rm B}/\dot{M}_{\rm Edd}\sim 10^{-8}{\rm -}10^{-5}$$ for M• ∼ 10 M⊙. Some of those isolated BHs in our Galaxy could be observed as high-energy sources (Fujita et al. 1998; Armitage & Natarajan 1999; Agol & Kamionkowski 2002; Ioka et al. 2017; Matsumoto, Teraki & Ioka 2018). However, since those previous studies assumed self-similar ADAF solutions, it is worth revisiting their arguments considering our self-consistent accretion solutions. Note that Perna et al. (2003) have discussed the same problem for neutron stars, and pointed out that the reduction of accretion rates from the Bondi rates can explain why isolated neutron stars accreting the ISM are rarely observed in X-rays. 5 SUMMARY We study radiatively inefficient rotating-accretion flows on to a BH with two-dimensional hydrodynamical simulations. We consider axisymmetric accretion flows and allow angular momentum transport adopting the α viscosity prescription. When the gas angular momentum is low enough to form a rotationally supported disc/torus within the Bondi radius (RB), we find a global steady accretion solution. The solution consists of three phases: (1) a rotational equilibrium distribution at r ∼ RB, where the density distribution follows ρ ∝ (1 + RB/r)3/2 and the gas is very subsonic; (2) a geometrically thick accretion torus at the centrifugal radius RC( < RB), where thermal energy generated by viscosity is transported via strong convection motions expected in CDAF; and (3) an estimated inner optically thin disc controlled by conduction and viscosity. Physical properties of the accretion flow in the intermediate region (2) agree with those in CDAFs (e.g. ρ ∝ r−1/2). In the CDAF solution, the gas accretion rate decreases towards the centre due to convection ($$\dot{M}\propto r$$), and the net accretion rate is strongly suppressed by several orders of magnitude from the Bondi accretion rate. We find that the net accretion rate depends on the viscous strength and the size of the innermost radius of the computational domain (rmin), and can be approximated as $$\dot{M}\sim (\alpha /0.01)^{0.6}(r_{\rm min}/R_{\rm B})\dot{M}_{\rm B}$$. This solution holds for low accretion rates of $$\dot{M}_{\rm B}/\dot{M}_{\rm Edd}\lesssim 10^{-3}$$ having minimal radiation cooling. In a hot plasma at the bottom of this solution (r < 10− 3 RB ≲ rmin), thermal conduction would dominate the convective energy flux. Since suppression of the accretion by convection ceases, the final BH feeding rate is found to be $$\dot{M}/\dot{M}_{\rm B} \sim 10^{-2}{\rm -}10^{-3}$$. This rate is as low as $$\dot{M}/\dot{M}_{\rm Edd} \sim 10^{-7}{\rm -}10^{-6}$$ inferred for SgrA* and the nuclear BHs in M31 and M87, and can explain the low accretion rates and low luminosities in these sources, without invoking any feedback mechanism. Finally, we briefly mention several effects which we do not treat in our simulations. First, stellar winds and X-ray emission from massive stars around the Bondi radius are neglected. These effects heat the gas to T ∼ keV and could drive outflows as in the Galactic centre around SgrA* (e.g. Najarro et al. 1997; Baganoff et al. 2003; Quataert 2004). Secondly, non-axisymmetric perturbations would lead to instability and allow additional angular momentum transport in the accretion flow (e.g. Papaloizou & Pringle 1984; Chandrasekhar 1961). In order to include these effects, we further need to conduct three-dimensional hydrodynamical simulations of accretion flows at large scales (e.g. Gaspari, Ruszkowski & Oh 2013). Thirdly, thermal conductivity in the inner region (r < Rtr ∼ 10−3 RB) plays an important role as an efficient process carrying the energy outwards instead of convection. In fact, the location of the transition radius would determine the net accretion rate on to the nuclear disc and the radiative luminosity produced in the vicinity of the BH as shown in Fig. 17. We need to investigate the detailed properties of radiatively inefficient accretion flows, in order to make a prediction for future observations with the Event Horizon Telescope.4 Acknowledgements We thank James Stone, Charles Gammie, Ramesh Narayan, Eliot Quataert, Lorenzo Sironi, Daniel Wang, Kengo Tomida, Kazumi Kashiyama, Kohei Ichikawa, Takashi Hosokawa, and Kazuyuki Sugimura for useful discussions. This work is partially supported by the Simons Foundation through the Simons Society of Fellows (KI), Simons Fellowship in Theoretical Physics (ZH), and NASA grant NNX15AB19G (ZH). 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. Numerical computations were carried out on Cray XC30 at the Center for Computational Astrophysics of the National Astronomical Observatory of Japan. Footnotes 1 For a spherically symmetric flow, the sonic point is estimated as Rsonic = [(5 − 3γ)/4]RB. Since for γ = 5/3, the sonic point is located at the origin, the inner boundary conditions would affect gas dynamics in the computational domain. Thus, we set a slightly smaller value of γ = 1.6, for which the sonic point can be resolved (Rsonic ≃ 0.05 RB). 2 The heating term due to convective viscosity and other higher order terms are neglected. Our purpose is not showing precise formulae but giving an order-of-magnitude estimate for the convective energy flux. 3 The bolometric luminosity of the M31 BH is estimated from the X-ray luminosity, assuming a bolometric correction factor of fbol = 10 for such a low-luminosity BH (Hopkins et al. 2007). 4 http://eventhorizontelescope.org REFERENCES Agol E., Kamionkowski M., 2002, MNRAS , 334, 553 CrossRef Search ADS   Armitage P. J., Natarajan P., 1999, ApJ , 523, L7 CrossRef Search ADS   Baganoff F. K. et al.  , 2003, ApJ , 591, 891 CrossRef Search ADS   Balbus S. A., 2000, ApJ , 534, 420 CrossRef Search ADS   Balbus S. A., Hawley J. F., 1991, ApJ , 376, 214 CrossRef Search ADS   Balbus S. A., Hawley J. F., 1998, Rev. Mod. Phys. , 70, 1 CrossRef Search ADS   Ball G. H., Narayan R., Quataert E., 2001, ApJ , 552, 221 CrossRef Search ADS   Bender R. et al.  , 2005, ApJ , 631, 280 CrossRef Search ADS   Blandford R. D., Begelman M. C., 1999, MNRAS , 303, L1 CrossRef Search ADS   Blandford R. D., Begelman M. C., 2004, MNRAS , 349, 68 CrossRef Search ADS   Bondi H., 1952, MNRAS , 112, 195 CrossRef Search ADS   Chandrasekhar S., 1961, Hydrodynamic and Hydromagnetic Stability , International Series of Monographs on Physics. Clarendon, Oxford Ciotti L., Ostriker J. P., 2001, ApJ , 551, 131 CrossRef Search ADS   Ciotti L., Ostriker J. P., Proga D., 2009, ApJ , 699, 89 CrossRef Search ADS   Fishbone L. G., Moncrief V., 1976, ApJ , 207, 962 CrossRef Search ADS   Foucart F., Chandra M., Gammie C. F., Quataert E., 2016, MNRAS , 456, 1332 CrossRef Search ADS   Foucart F., Chandra M., Gammie C. F., Quataert E., Tchekhovskoy A., 2017, MNRAS , 470, 2240 CrossRef Search ADS   Fujita Y., Inoue S., Nakamura T., Manmoto T., Nakamura K. E., 1998, ApJ , 495, L85 CrossRef Search ADS   Garcia M. R. et al.  , 2010, ApJ , 710, 755 CrossRef Search ADS   Gaspari M., Ruszkowski M., Oh S. P., 2013, MNRAS , 432, 3401 CrossRef Search ADS   Gebhardt K., Adams J., Richstone D., Lauer T. R., Faber S. M., Gültekin K., Murphy J., Tremaine S., 2011, ApJ , 729, 119 CrossRef Search ADS   Ghez A. M. et al.  , 2003, ApJ , 586, L127 CrossRef Search ADS   Hawley J. F., Balbus S. A., Stone J. M., 2001, ApJ , 554, L49 CrossRef Search ADS   Ho L. C., 2008, ARA&A , 46, 475 CrossRef Search ADS   Ho L. C., 2009, ApJ , 699, 626 CrossRef Search ADS   Hopkins P. F., Richards G. T., Hernquist L., 2007, ApJ , 654, 731 CrossRef Search ADS   Ichimaru S., 1977, ApJ , 214, 840 CrossRef Search ADS   Igumenshchev I. V., Abramowicz M. A., 1999, MNRAS , 303, 309 CrossRef Search ADS   Igumenshchev I. V., Abramowicz M. A., 2000, ApJS , 130, 463 CrossRef Search ADS   Igumenshchev I. V., Abramowicz M. A., Narayan R., 2000, ApJ , 537, L27 CrossRef Search ADS   Igumenshchev I. V., Narayan R., Abramowicz M. A., 2003, ApJ , 592, 1042 CrossRef Search ADS   Inayoshi K., Haiman Z., Ostriker J. P., 2016, MNRAS , 459, 3738 CrossRef Search ADS   Ioka K., Matsumoto T., Teraki Y., Kashiyama K., Murase K., 2017, MNRAS , 470, 3332 CrossRef Search ADS   Johnson B. M., Quataert E., 2007, ApJ , 660, 1273 CrossRef Search ADS   Kato S., Fukue J., Mineshige S., 2008, Black-Hole Accretion Disks – Towards a New Paradigm  –. Kyoto Univ. Press, Kyoto, Japan Kuiper R., Klahr H., Beuther H., Henning T., 2010, ApJ , 722, 1556 CrossRef Search ADS   Kuiper R., Klahr H., Beuther H., Henning T., 2011, ApJ , 732, 20 CrossRef Search ADS   Kuo C. Y. et al.  , 2014, ApJ , 783, L33 CrossRef Search ADS   Li J., Ostriker J., Sunyaev R., 2013, ApJ , 767, 105 CrossRef Search ADS   Machida M., Matsumoto R., Mineshige S., 2001, PASJ , 53, L1 CrossRef Search ADS   Maron J., Chandran B. D., Blackman E., 2004, Phys. Rev. Lett. , 92, 045001 CrossRef Search ADS PubMed  Matsumoto R., Tajima T., 1995, ApJ , 445, 767 CrossRef Search ADS   Matsumoto T., Teraki Y., Ioka K., 2018, MNRAS , 475, 1251 CrossRef Search ADS   McCourt M., Parrish I. J., Sharma P., Quataert E., 2011, MNRAS , 413, 1295 CrossRef Search ADS   McKinney J. C., Gammie C. F., 2004, ApJ , 611, 977 CrossRef Search ADS   Mignone A., Bodo G., Massaglia S., Matsakos T., Tesileanu O., Zanni C., Ferrari A., 2007, ApJS , 170, 228 CrossRef Search ADS   Milosavljević M., Bromm V., Couch S. M., Oh S. P., 2009, ApJ , 698, 766 CrossRef Search ADS   Najarro F., Krabbe A., Genzel R., Lutz D., Kudritzki R. P., Hillier D. J., 1997, A&A , 325, 700 Narayan R., Fabian A. C., 2011, MNRAS , 415, 3721 CrossRef Search ADS   Narayan R., Medvedev M. V., 2001, ApJ , 562, L129 CrossRef Search ADS   Narayan R., Yi I., 1994, ApJ , 428, L13 CrossRef Search ADS   Narayan R., Yi I., 1995, ApJ , 444, 231 CrossRef Search ADS   Narayan R., Igumenshchev I. V., Abramowicz M. A., 2000, ApJ , 539, 798 CrossRef Search ADS   Narayan R., Sa̧dowski A., Penna R. F., Kulkarni A. K., 2012, MNRAS , 426, 3241 CrossRef Search ADS   Ohsuga K., Mineshige S., Mori M., Kato Y., 2009, PASJ , 61, L7 CrossRef Search ADS   Ostriker J. P., Weaver R., Yahil A., McCray R., 1976, ApJ , 208, L61 CrossRef Search ADS   Papaloizou J. C. B., Pringle J. E., 1984, MNRAS , 208, 721 CrossRef Search ADS   Park M.-G., 1990, ApJ , 354, 64 CrossRef Search ADS   Park K., Ricotti M., 2011, ApJ , 739, 2 CrossRef Search ADS   Perna R., Narayan R., Rybicki G., Stella L., Treves A., 2003, ApJ , 594, 936 CrossRef Search ADS   Pringle J. E., 1981, ARA&A , 19, 137 CrossRef Search ADS   Proga D., 2007, ApJ , 661, 693 CrossRef Search ADS   Proga D., Begelman M. C., 2003, ApJ , 592, 767 CrossRef Search ADS   Quataert E., 2004, ApJ , 613, 322 CrossRef Search ADS   Quataert E., 2008, ApJ , 673, 758 CrossRef Search ADS   Quataert E., Gruzinov A., 2000, ApJ , 539, 809 CrossRef Search ADS   Roberts S. R., Jiang Y.-F., Wang Q. D., Ostriker J. P., 2017, MNRAS , 466, 1477 CrossRef Search ADS   Russell H. R., Fabian A. C., McNamara B. R., Broderick A. E., 2015, MNRAS , 451, 588 CrossRef Search ADS   Ryan B. R., Ressler S. M., Dolence J. C., Tchekhovskoy A., Gammie C., Quataert E., 2017, ApJ , 844, L24 CrossRef Search ADS   Sa̧dowski A., Wielgus M., Narayan R., Abarca D., McKinney J. C., Chael A., 2017, MNRAS , 466, 705 CrossRef Search ADS   Sano T., Inutsuka S.-I., Turner N. J., Stone J. M., 2004, ApJ , 605, 321 CrossRef Search ADS   Shakura N. I., Sunyaev R. A., 1973, A&A , 24, 337 Shapiro S. L., 1973, ApJ , 180, 531 CrossRef Search ADS   Sharma P., Quataert E., Stone J. M., 2008, MNRAS , 389, 1815 CrossRef Search ADS   Shcherbakov R. V., Baganoff F. K., 2010, ApJ , 716, 504 CrossRef Search ADS   Spitzer L., 1962, Physics of Fully Ionized Gases . Interscience, New York Stone J. M., Norman M. L., 1992, ApJS , 80, 753 CrossRef Search ADS   Stone J. M., Pringle J. E., 2001, MNRAS , 322, 461 CrossRef Search ADS   Stone J. M., Hawley J. F., Gammie C. F., Balbus S. A., 1996, ApJ , 463, 656 CrossRef Search ADS   Stone J. M., Pringle J. E., Begelman M. C., 1999, MNRAS , 310, 1002 CrossRef Search ADS   Walsh J. L., Barth A. J., Ho L. C., Sarzi M., 2013, ApJ , 770, 86 CrossRef Search ADS   Yuan F., Quataert E., Narayan R., 2003, ApJ , 598, 301 CrossRef Search ADS   Yuan F., Wu M., Bu D., 2012, ApJ , 761, 129 CrossRef Search ADS   APPENDIX: NUMERICAL TESTS In this appendix, we discuss the effects of numerical resolutions, outer boundary conditions, initial conditions, prescriptions of viscosity on our results. As a reference, we define our fiducial case with RC/RB = 0.1 and α = 0.01, where (1) the number of grid cells is Nr × Nθ = 512 × 256, (2) the initial condition is given by a rotating equilibrium density configuration by equation (17), and (3) the additional viscosity which would be driven by rotational instability is assumed (see equation 16). Fig. A1 shows the time evolution of the net accretion rates on to the BH (i.e. a sink cell at the centre) with different numerical resolutions: Nr × Nθ = 512 × 256 (fiducial, red), 256 × 256 (low resolution, green), and 1024 × 512 (high resolution, blue). For the three cases, the time-averaged accretion rate is $$\dot{M}/\dot{M}_{\rm B}=7.1\times 10^{-3}$$ (red), 7.9 × 10−3 (green), and 6.3 × 10−3 (blue), respectively. The estimated errors are at most ∼ 20 per cent. Figure A1. View largeDownload slide Time evolution of the net accretion on to a sink cell (in units of the Bondi rate) for RC/RB = 0.1 and α = 0.01. Simulation results with three different resolutions are shown: Nr × Nθ = 512 × 256 (fiducial, red), 256 × 256 (low resolution, green), and 1024 × 512 (high resolution, blue). Figure A1. View largeDownload slide Time evolution of the net accretion on to a sink cell (in units of the Bondi rate) for RC/RB = 0.1 and α = 0.01. Simulation results with three different resolutions are shown: Nr × Nθ = 512 × 256 (fiducial, red), 256 × 256 (low resolution, green), and 1024 × 512 (high resolution, blue). As seen in Figs 2 and 3, a significant fraction of the gas is outflowing towards the polar directions and the rate dominates the inflow rate outside the Bondi radius. In order to ensure stationarity of the accretion system, we study the time evolution of the outflow rate. Fig. A2 shows radial profiles of the outflow rates (i.e. vr > 0, see equation 20) for different elapsed times (6 ≤ t/tdyn ≤ 23). This component propagates outwards, accumulating the ambient gas which has a uniform distribution and no accretion initially at r > RB. However, the outflow rates inside the Bondi radius seems steady, and the rates even outside the Bondi radius have nearly converged (r ≲ 7 RB). Therefore, our long-term simulations allow us to ensure the stationarity of the accretion flow in the inner region. Figure A2. View largeDownload slide Time evolution of radial profiles of the outflow rates for RC/RB = 0.1 and α = 0.01. Different curves show the time evolution: t/tdyn = 6 (1, red), 10 (2, green), 14 (3, blue), 18 (4, magenta), and 23 (5, light blue). Although the outflow propagates, accumulating the ambient gas, the outflow rate inside the Bondi radius is steady. This result ensures the stationarity of the accretion flow in the inner region. Figure A2. View largeDownload slide Time evolution of radial profiles of the outflow rates for RC/RB = 0.1 and α = 0.01. Different curves show the time evolution: t/tdyn = 6 (1, red), 10 (2, green), 14 (3, blue), 18 (4, magenta), and 23 (5, light blue). Although the outflow propagates, accumulating the ambient gas, the outflow rate inside the Bondi radius is steady. This result ensures the stationarity of the accretion flow in the inner region. Fig. A3 shows the effect of initial conditions on the time evolution of the accretion rate. Red and green curve present the fiducial case and the result with a Bondi accretion solution added a constant angular momentum j∞ as the initial conditions. In the latter case, the accretion rate decreases rapidly by several orders of magnitude in the early stage due to the centrifugal force, begins to increase gradually at t ≳ 2 tdyn and approaches that in the fiducial case at t ≳ 6 tdyn. Therefore, this result clearly shows that the net accretion rate does not depend on the choice of the initial conditions. Figure A3. View largeDownload slide Time evolution of the net accretion on to a sink cell (in units of the Bondi rate) for RC/RB = 0.1 and α = 0.01. Red curve presents the fiducial case as shown in Fig. 1. Green one is the result where the simulation starts from the Bondi solution with a constant angular momentum j∞. Blue one is the result where the second term in equation (16), which provides additional viscosity led by rotational instability, is neglected. Figure A3. View largeDownload slide Time evolution of the net accretion on to a sink cell (in units of the Bondi rate) for RC/RB = 0.1 and α = 0.01. Red curve presents the fiducial case as shown in Fig. 1. Green one is the result where the simulation starts from the Bondi solution with a constant angular momentum j∞. Blue one is the result where the second term in equation (16), which provides additional viscosity led by rotational instability, is neglected. Finally, we study the effect of the additional viscosity given by the second term in equation (16), which would be led by rotational instability. Blue curve in Fig. A3 shows the accretion rate without the additional term, and almost identical to that in the fiducial case. Fig. A4 shows radial profiles of the specific angular momentum normalized by $$j_\infty (=\sqrt{\beta }R_{\rm B}c_\infty )$$ for RC/RB = 0.1 and α = 0.01. Solid (red) curve presents the result at t = 13 tdyn in our fiducial case. The profile has a Keplerian-like profile expected from a CDAF solution inside the centrifugal radius, $$j=g(\gamma )\sqrt{GM_\bullet r}$$ (dotted black), and a constant value outside the torus (r ≳ 2 RC). Two dashed curves show the cases without additional viscosity in the second term in equation (16) at two different elapsed time of t = 13 tdyn (green) and t = 16 tdyn (blue). Compared to the fiducial case, the profile at the same elapsed time (green curve) have an excess of the angular momentum from the boundary value of j∞ outside the torus. Since viscosity does not work outside the torus, the angular momentum is accumulated and the bump structure grows with time (see blue curve). As discussed in Section 2, the regions with negative gradients of the angular momentum (i.e. dj/dr < 0) are unstable and would produce turbulence, which allows additional angular momentum transport. In order to capture the physics of the instability, we need three-dimensional hydrodynamical simulations. However, we have confirmed that the treatment of the rotational instability does not affect our results. Figure A4. View largeDownload slide Radial profiles of the specific angular momentum normalized by $$j_\infty (=\sqrt{\beta }R_{\rm B}c_\infty )$$ for RC/RB = 0.1 and α = 0.01. Solid (red) curve presents the profile at t = 13 tdyn in our fiducial case, where the second term in equation (16) is included to provide additional viscosity led by rotational instability. Dashed curves show the results at t = 13 tdyn (green) and 16 tdyn (blue), respectively, in cases where the additional viscosity is neglected. Dotted (black) curve is that expected in a CDAF solution, $$j=g(\gamma )\sqrt{GM_\bullet r}$$. Figure A4. View largeDownload slide Radial profiles of the specific angular momentum normalized by $$j_\infty (=\sqrt{\beta }R_{\rm B}c_\infty )$$ for RC/RB = 0.1 and α = 0.01. Solid (red) curve presents the profile at t = 13 tdyn in our fiducial case, where the second term in equation (16) is included to provide additional viscosity led by rotational instability. Dashed curves show the results at t = 13 tdyn (green) and 16 tdyn (blue), respectively, in cases where the additional viscosity is neglected. Dotted (black) curve is that expected in a CDAF solution, $$j=g(\gamma )\sqrt{GM_\bullet r}$$. © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Monthly Notices of the Royal Astronomical Society Oxford University Press

Low-density, radiatively inefficient rotating-accretion flow on to a black hole

Loading next page...
 
/lp/ou_press/low-density-radiatively-inefficient-rotating-accretion-flow-on-to-a-uWeDbiqX30
Publisher
Oxford University Press
Copyright
© 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society
ISSN
0035-8711
eISSN
1365-2966
D.O.I.
10.1093/mnras/sty276
Publisher site
See Article on Publisher Site

Abstract

Abstract We study low-density axisymmetric accretion flows on to black holes (BHs) with two-dimensional hydrodynamical simulations, adopting the α-viscosity prescription. When the gas angular momentum is low enough to form a rotationally supported disc within the Bondi radius (RB), we find a global steady accretion solution. The solution consists of a rotational equilibrium distribution around r ∼ RB, where the density follows ρ ∝ (1 + RB/r)3/2, surrounding a geometrically thick and optically thin accretion disc at the centrifugal radius RC(<RB), where thermal energy generated by viscosity is transported via convection. Physical properties of the inner solution agree with those expected in convection-dominated accretion flows (ρ ∝ r−1/2). In the inner solution, the gas inflow rate decreases towards the centre due to convection ($$\dot{M}\propto r$$), and the net accretion rate (including both inflows and outflows) is strongly suppressed by several orders of magnitude from the Bondi accretion rate $$\dot{M}_{\rm B}$$. The net accretion rate depends on the viscous strength, following $$\dot{M}/\dot{M}_{\rm B}\propto (\alpha /0.01)^{0.6}$$. This solution holds for low accretion rates of $$\dot{M}_{\rm B}/\dot{M}_{\rm Edd}\lesssim 10^{-3}$$ having minimal radiation cooling, where $$\dot{M}_{\rm Edd}$$ is the Eddington accretion rate. In a hot plasma at the bottom (r < 10−3 RB), thermal conduction would dominate the convective energy flux. Since suppression of the accretion by convection ceases, the final BH feeding rate is found to be $$\dot{M}/\dot{M}_{\rm B}\sim 10^{-3}$$–10−2. This rate is as low as $$\dot{M}/\dot{M}_{\rm Edd}\sim 10^{-7}$$–10−6 inferred for SgrA* and the nuclear BHs in M31 and M87, and can explain their low luminosities, without invoking any feedback mechanism. accretion, accretion discs, black hole physics, X-rays: galaxies 1 INTRODUCTION The basic physics of gas accretion on to an astrophysical object was first studied for a spherically symmetric flow (Bondi 1952). When rotation is not included, gas accretion on to a black hole (BH) with a mass of M• begins from a characteristic radius, where the negative gravitational energy becomes greater than the thermal energy of the gas. The so-called Bondi radius is given by   \begin{equation} R_{\rm B}\equiv \frac{GM_\bullet }{c_{\rm s}^2}, \end{equation} (1)where cs is the sound speed of the gas. One can safely assume roughly free fall within this radius at a rate of $$\dot{M}_{\rm B}$$, so-called the Bondi rate. If radiation emitted from the inner region of the flow exceeds several percent of the Eddington luminosity (LEdd ≡ 4πcGM•/κes, where κes is the electron scattering opacity, then this feedback can add enough energy and momentum at the Bondi radius to reverse the flow (Ostriker et al. 1976; Ciotti & Ostriker 2001; Proga 2007; Ciotti, Ostriker & Proga 2009; Milosavljević et al. 2009; Park & Ricotti 2011; Inayoshi, Haiman & Ostriker 2016). However, the problem is rather simple for low accretion rates (Shapiro 1973; Park 1990). Now let us add angular momentum per unit mass j to the flow. This effect introduces a new physical scale of   \begin{equation} R_{\rm C}\equiv \frac{j^2}{GM_\bullet }, \end{equation} (2)which is the centrifugal radius where the centrifugal force balances the gravity of the BH. If this radius is larger than the Schwarzschild radius   \begin{equation} R_{\rm Sch}\equiv \frac{GM_\bullet }{c^2}, \end{equation} (3)the accreting gas will simply settle into a rotational equilibrium distribution with no net inflow (Fishbone & Moncrief 1976; Papaloizou & Pringle 1984; Li, Ostriker & Sunyaev 2013). However, if viscosity exists, then angular momentum transport is possible and an accretion disc forms at the bottom of the distribution. The effective shear viscosity driven by magneto rotational instability (MRI) is often described with the standard α-prescription (e.g. Shakura & Sunyaev 1973) as   \begin{equation} \nu =\alpha c_{\rm s}H, \end{equation} (4)where H is the scale height of the inner flow and α presents the strength of viscosity estimated as α ∼ O(10−2) by magnetohydrodynamical (MHD) simulations (Balbus & Hawley 1991; Matsumoto & Tajima 1995; Stone et al. 1996; Balbus & Hawley 1998; Sano et al. 2004). There have been a number of analytical solutions and numerical studies for rotating flows with viscous angular momentum transport (e.g. Shakura & Sunyaev 1973; see reviews by Pringle 1981; Kato, Fukue & Mineshige 2008, references therein). Among those, we here focus on accretion flows which cannot lose internal energy via radiative cooling because of very low gas density. Such radiatively inefficient accretion flows are quite interesting since many observed BHs accrete at rates of only a small fraction of the Bondi accretion rate and their radiation luminosity is as low as ∼10−1 to 10−9 LEdd (Ho 2008, 2009). Sagittarius A* (Sgr A*) is inferred to be accreting at a rate of 10−3 to $$10^{-2}\ \dot{M}_{\rm B}$$ (e.g. Yuan, Quataert & Narayan 2003; Quataert 2004), where $$\dot{M}_{\rm B}\simeq 10^{-5}\ {\rm M}_{\odot }\, {\rm yr}^{-1}$$ is measured from the temperature and density near the Bondi radius with X-ray observations (Baganoff et al. 2003). Because of such a low accretion rate, the bolometric luminosity of Sgr A* (M• ≃ 4 × 106 M⊙; Ghez et al. 2003) is as small as Lbol ∼ 1036 erg s−1 ∼ 2 × 10−9 LEdd. The second example is a BH at the centre of the giant elliptical galaxy M87. The gas accretion rate at the vicinity of the BH is estimated as ≲ 9.2 × 10− 4 M⊙ yr− 1 (Kuo et al. 2014), which is lower than $${\sim } 10^{-2}\ \dot{M}_{\rm B}$$ (Russell et al. 2015). Since the BH mass is estimated as $$M_\bullet = 6.6^{+0.4}_{-0.4}\times 10^9\ {\rm M}_{\odot }$$ (Gebhardt et al. 2011) and $$M_\bullet = 3.5^{0.9}_{-0.7}\times 10^9\ {\rm M}_{\odot }$$ (Walsh et al. 2013), the bolometric luminosity of Lbol ∼ 2 × 1041 erg s−1 is ∼3 × 10−7 LEdd. The third example is a BH at the centre of the Andromeda Galaxy (M31). The estimated BH mass is $$M_\bullet \simeq 1.4^{+0.7}_{-0.3}\times 10^8\ {\rm M}_{\odot }$$ (Bender et al. 2005). The Bondi accretion rate and the X-ray luminosity are estimated as $$\dot{M}_{\rm B}\simeq 7\times 10^{-5}\ {\rm M}_{\odot }\, {\rm yr}^{-1}$$ and LX ≃ 2 × 1036 erg s−1 ≃ 10−10 LEdd, respectively (Garcia et al. 2010). The corresponding bolometric luminosity is inferred as ≃10−9 LEdd by assuming the bolometric correction factor to be ≃10 (Hopkins, Richards & Hernquist 2007). In very low-accretion-rate flows, the gas never cools via radiation and the accretion disc becomes hot and thick. This thick-disc solution has been found by Ichimaru (1977) and studied in the subsequent works by Narayan & Yi (1994, 1995). This solution is well known in a simplified self-similar version as advection-dominated accretion flows (ADAFs). Their physical properties can be understood by considering a self-similar solution, which suggests that low-density accretion flow is hot and geometrically extended in the polar direction and the thermal energy is advected with the flow on to the BH. Moreover, as a related solution, the so-called adiabatic inflow–outflow solution has been proposed by Blandford & Begelman (1999, 2004), where matter can accrete at very low rates due to strong outflows which carry energy and angular momentum away. As pointed out by Narayan & Yi (1994, 1995), ADAF solutions tend to be convectively unstable because the gas entropy increases towards the centre due to thermal energy release via viscosity. Narayan, Igumenshchev & Abramowicz (2000) and Quataert & Gruzinov (2000) conducted stability analyses for rotating adiabatic gas against convection motions, and found that the density profile in a marginally stable state follows ρ ∝ r−1/2, whose slope is rather flatter than ρ ∝ r−3/2 expected from ADAFs. This solution is known as convection-dominated accretion flows (CDAF). Numerical simulations of radiatively inefficient accretion flows have also supported that CDAF solutions can exist in a weak viscosity and no cooling regime (Igumenshchev & Abramowicz 1999, 2000; Stone, Pringle & Begelman 1999; Igumenshchev, Narayan & Abramowicz 2003). Many previous works with multidimensional MHD simulations have studied the properties of the accretion flow at the vicinity of the central BH (r ≲ a few × 102 RSch) starting from torus-like initial reservoirs of gas (Hawley, Balbus & Stone 2001; Machida, Matsumoto & Mineshige 2001; Stone & Pringle 2001; Igumenshchev, Narayan & Abramowicz 2003; McKinney & Gammie 2004; Ohsuga et al. 2009; Narayan et al. 2012; see also Yuan, Wu & Bu 2012, who studied the gas dynamics over a wider range of spacial scales). However, because of the limitation of the computational domain, it remains unclear how and whether the accretion flow, which begins from large radii (r ∼ RB), can lose the angular momentum and reach smaller radii (r ≪ RB), where the gas is tightly bound to the central BH. To explicitly connect large and small scales in a self-consistent solution, we study dynamics of gas accretion on to a BH over a wide range of spatial scales, and are able to connect the inner region, where a disc forms, with the region well outside the Bondi radius, where the accretion originates in the first place. Li et al. (2013) studied a similar problem for one specific parameter set, i.e. RC/RB = 0.02 and α ∼ 10−3, with two-dimensional hydrodynamical simulations assuming equatorial symmetry. This work followed gas accretion from outside the Bondi radius to within the centrifugal radius. Proga & Begelman (2003) have conducted two-dimensional MHD simulations for this problem with a more limited dynamic range (r < 103 RSch), and found that gas accretion is allowed by viscosity driven by the MRI. However, several crucial questions remain unanswered. How do physical parameters, e.g. RC/RB and α, determine the actual inflow rate on to a BH? What do physical properties of the accretion flow inside the centrifugal radius look like? We perform two-dimensional viscous-hydrodynamical simulations with a large dynamic range including angular momentum transport with the α-viscosity prescription, and address the above questions for radiatively inefficient rotating-accretion flows on to a BH. Here, we focus on cases where the gas angular momentum is low enough to form a rotationally supported disc/torus within the Bondi radius, i.e.   \begin{equation} R_{\rm Sch}\ll R_{\rm C}<R_{\rm B}. \end{equation} (5)In this case, we find a global steady accretion solution which consists of three parts. In the outer region (r ∼ RB), a rotational equilibrium distribution is developed where the density distribution follows ρ ∝ (1 + RB/r)3/2 and the gas motion is very subsonic. The solution deviates only slightly from the rotational equilibrium flow with no radial motion (Fishbone & Moncrief 1976). Interior to this (r ≲ 2 RC), a geometrically thick accretion torus is formed, where thermal energy generated by viscosity is transported via strong convection motions. Physical properties of the inner solution agree with those expected in a CDAF solution, e.g. ρ ∝ r−1/2 (Narayan et al. 2000; Quataert & Gruzinov 2000). In the inner CDAF solution, the gas inflow rate decreases towards the centre due to convection ($$\dot{M}\propto r$$), and the net accretion rate (including both inflows and outflows) is strongly suppressed by several orders of magnitude from the Bondi accretion rate. We also study the effect of the viscous strength on the results and find that the net accretion rate is approximated as $$\dot{M}\sim (\alpha /0.01)^{0.6}(r/R_{\rm B})\dot{M}_{\rm B}$$. Finally, in a hot plasma at the bottom of this solution (r ≲ 10− 3 RB), thermal conduction would dominate the convective energy flux and the flow would resemble an optically thin accretion disc. Since suppression of the accretion by convection ceases, the final BH feeding rate is determined as $$\dot{M}\sim 10^{-3}-10^{-2}\ \dot{M}_{\rm B}$$. This rate can explain why Sgr A* and the BHs in M31 and M87 accrete at such low accretion rates, without invoking any feedback mechanism. The rest of this paper is organized as follows. In Section 2, we describe the methodology of our numerical simulations. In Section 3, we show our simulation results and explain their physical properties. We compare the results to those in CDAF solutions and quantify the effects of energy and angular momentum transport by convection (Section 3.1), and study the effect of the choice of viscous parameters on the results (Section 3.2). In Section 4, we give an analytical argument extending the accretion solution further inwards and discuss a global solution of radiatively inefficient rotating-accretion flows. In Section 5, we summarize our main conclusions and discuss the importance for observations of very low-luminosity BHs. 2 METHODOLOGY 2.1 Basic quantities We introduce basic dimensionless physical quantities which characterize BH accretion systems. If feedback by radiation and/or momentum from the central BH is negligible, gas accretion begins from the Bondi radius, RB. The standard accretion rate from the Bondi radius is given by   \begin{eqnarray} \dot{M}_{\rm B}&=& 4\pi q(\gamma ) \rho _\infty \frac{G^2M_\bullet ^2}{c_{\infty }^3}, \nonumber \\ &\simeq & 8.8\times 10^{-5}\ \rho _{-22} M_6^2 c_7^{-3}\ {\rm M}_{\odot }\, {\rm yr}^{-1}, \end{eqnarray} (6)where q(γ) = 1/4 for γ = 5/3, $$\rho _{-22}=\rho _\infty /(10^{-22}\ {\rm g\ {\rm cm}^{-3}})$$, M6 = M•/(106 M⊙) and c7 = c∞/(107 cm s−1) and the accretion rate normalized by the Eddington rate is given by   \begin{eqnarray} \dot{m}_{\rm B}\equiv \frac{\dot{M}_{\rm B}}{\dot{M}_{\rm Edd}} =3.8\times 10^{-3}\ \rho _{-22} M_6 c_7^{-3}, \end{eqnarray} (7)where $$\dot{M}_{\rm Edd}(\equiv 10\ L_{\rm Edd}/c^2) = 2.3\times 10^{-2}\ M_6\ {\rm M}_{\odot }\, {\rm yr}^{-1}$$ is the Eddington accretion rate with a 10 per cent radiative efficiency. As shown in Section 1, there are three characteristic physical scales: the Bondi radius RB, the centrifugal radius RC and the Schwarzschild radius RSch. From these three radii, we can define two ratios as   \begin{equation} \frac{R_{\rm Sch}}{R_{\rm B}}=\frac{2c_\infty ^2}{c^2}, \end{equation} (8)and   \begin{equation} \frac{R_{\rm C}}{R_{\rm B}}=\frac{j^2c_\infty ^2}{G^2 M_\bullet ^2}. \end{equation} (9)Assuming a constant specific angular momentum of $$j=\sqrt{\beta }R_{\rm B}c_\infty$$, the ratio is written as RC/RB = β. Since we discuss cases where RSch ≪ RC < RB, we assume β < 1. In axisymmetric two-dimensional hydrodynamical simulations, we need angular momentum transport via viscous processes in order to allow gas accretion on to a BH at the centre. We here model the effects of viscosity with the standard α-prescription proposed by Shakura & Sunyaev (1973) instead of solving the time evolution of MHD equations. According to MHD simulations, effective viscosity driven by the MRI provides angular momentum transport and the strength is estimated as α ∼ O(10−2) (Balbus & Hawley 1991; Matsumoto & Tajima 1995; Stone et al. 1996; Balbus & Hawley 1998; Sano et al. 2004). Throughout this paper, we focus on adiabatic gas without radiative cooling in order to keep numerical results scale-free. The assumptions are justified for very low accretion-rate system, i.e. $$\dot{M}_{\rm B}\ll \dot{M}_{\rm Edd}$$ (see discussion in Section 4). In this limit, we have only three important non-dimensional parameters of RSch/RB, RC/RB( = β) and α to define a BH accretion system. 2.2 Basic equations We solve the axisymmetric two-dimensional hydrodynamical equations using an open code pluto (Mignone et al. 2007; Kuiper et al. 2010, 2011). The basic equations are following: the equation of continuity,   \begin{equation} \frac{{\rm d}\rho }{{\rm d}t}+\rho \nabla \cdot \boldsymbol v=0, \end{equation} (10)and the equation of motion,   \begin{equation} \rho \frac{{\rm d}\boldsymbol v}{{\rm d}t}=-\nabla p -\rho \nabla \Phi + \nabla \cdot \boldsymbol \sigma , \end{equation} (11)where ρ is the density, $$\boldsymbol v$$ is the velocity, and p is the gas pressure, the gravitational potential is set to Φ = −GM•/r, and $$\boldsymbol \sigma$$ is the stress tensor due to viscosity. The time derivative is the Lagrangian derivative, given by d/dt ≡ ∂/∂t + $$\boldsymbol v$$· ∇. We solve the energy equation of   \begin{equation} \rho \frac{{\rm d}e}{{\rm d}t}=-p\nabla \cdot \boldsymbol v + (\boldsymbol \sigma \cdot \nabla ) \boldsymbol v, \end{equation} (12)where e is the internal energy per mass. The equation of state of the ideal gas is assumed as p = (γ − 1)ρe, where γ = 1.6.1 The two terms on the right-hand side present compressional heating (or expansion cooling) and viscous heating. The viscous stress tensor is given by   \begin{equation} \sigma _{ij}=\rho \nu \left[ \left(\frac{\mathrm{\partial} v_j}{\mathrm{\partial} x_i} + \frac{\mathrm{\partial} v_i}{\mathrm{\partial} x_j}\right) -\frac{2}{3} (\nabla \cdot \boldsymbol v)\delta _{ij} \right], \end{equation} (13)where ν is the shear viscosity. Note that the bulk viscosity is neglected here. The shear viscosity is calculated with the α-prescription (Shakura & Sunyaev 1973),   \begin{equation} \nu = \alpha \frac{c_{\rm s}^2}{\Omega _{\rm K}}, \end{equation} (14)where ΩK = (GM•/r3)1/2. We calculate the viscous parameter by mimicking some properties of the MRI as   \begin{equation} \alpha = \alpha _0 \exp {\left[ -\left(\frac{\rho _{\rm cr}}{\rho }\right)^2\right]}, \end{equation} (15)where ρ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 Section 2.3). Under this model, the viscous process is active almost only within an accretion disc (r ≲ RC), where the rotational velocity has a significant fraction of the Keplerarian value. On the other hand, since the viscosity becomes zero outside the disc, angular momentum transported from the disc tends to be accumulated outside it as shown by the ‘bump’ in the angular momentum distribution in Fig. A4, where no viscous processes operate. Exterior to the peak in the angular momentum distribution where the specific angular momentum has a negative gradient outwards, i.e. ∂j/∂r < 0, so-called Rayleigh's criterion. In reality, such rotating flows are unstable and become turbulent (e.g. Chandrasekhar 1961), leading to angular momentum transport in three-dimensional simulations. However, because of limitations of our two-dimensional simulations, this could not occur and angular momentum flowing out from the central regions would accumulate outside of the region where ρ/ρcr ∼ 1. Thus, we modify equation (15) by adding the second term as   \begin{equation} \alpha = \alpha _0 \left\lbrace \exp {\left[ -\left(\frac{\rho _{\rm cr}}{\rho }\right)^2\right]} + {\rm max}\left(0, -\frac{\mathrm{\partial} \ln j}{\mathrm{\partial} \ln r}\right)\right\rbrace , \end{equation} (16)In a physical sense, the second term is necessary so that a steady state of the accretion flow exists. Note that the treatment of the rotational instability does not affect our results (see Appendix). 2.3 Boundary and initial conditions We set a computational domain of rmin ≤ r ≤ rmax and ε ≤ θ ≤ π − ε, where ε is set to 0.01 radian to avoid numerical singularity at poles. We set logarithmically spaced grids in the radial direction and uniformly spaced grids in the polar direction. The number of the grid points is set to (Nr, Nθ) = (512, 256). We have checked convergence for our simulation results, changing the number of the grids (see Appendix). As our fiducial case, we set rmin = 7.7 × 10−3 RB and rmax = 54 RB. Moreover, we calculate two cases with rmin = 3.0 × 10−3 RB and 1.3 × 10−2 RB (see Sections 3.1.2 and 3.2). Note that the relative size of each cell is given by $$\Delta r/r=(r_{\rm max}/r_{\rm min})^{1/N_r}$$. We adopt the outflow boundary condition at the innermost/outermost grid (e.g. 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. We also impose vr ≤ 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θ. As initial conditions, we adopt a rotational equilibrium distribution (vr = vθ = 0) with a constant specific angular momentum of j∞,   \begin{equation} \frac{\rho }{\rho _\infty }=\left[1+(\gamma -1)\frac{GM_\bullet }{c_\infty ^2 r} - \frac{(\gamma -1)}{2}\frac{j_\infty ^2}{c_{\rm \infty }^2\varpi ^2} \right] ^{\frac{1}{\gamma -1}} \end{equation} (17)(Papaloizou & Pringle 1984; see also Fishbone & Moncrief 1976), where ϖ = r sin θ is the cylindrical radius. Hereafter, we refer to this as the ‘Fishbone–Moncrief’ solution. 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∞, which leads to a maximum value of the density at r = RC and θ = 0   \begin{equation} \rho _{\rm cr}=\rho _\infty \left(1+\frac{\gamma -1}{2\beta }\right)^{1/(\gamma -1)}, \end{equation} (18)where β = RC/RB as defined below equation (9). Without viscosity, the density never exceeds this value because of the centrifugal barrier. In other words, a high-density region with ρ > ρcr must be formed by angular momentum transport due to viscosity. In order to study the dependence of our numerical results on the choice of the initial conditions, we also conduct a simulation which starts from the Bondi accretion solution with a constant angular momentum j∞. Note that the numerical result with the Bondi profiles as initial conditions approaches that adopting a rotational equilibrium distribution given by equation (17) (see Appendix). 3 RESULTS In this section, we show results of our two-dimensional simulations for accretion flows. We set all physical quantities as M• = 106 M⊙, $$\rho _\infty = 10^{-22}\ {\rm g\ {\rm cm}^{-3}}$$, and c∞ = 100 km s−1 ($$\dot{m}_{\rm B}\simeq 3.8\times 10^{-3}$$). As already mentioned, in adiabatic cases, our results do not depend on $$\dot{m}_{\rm B}$$ but on the two dimensionless parameters of β = RC/RB and α. In Section 3.1, we first describe overall properties of accretion flows with different initial angular momentum β for our fiducial case with α = 0.01. In Section 3.2, the dependence of our results on the choice of viscous parameter is discussed. 3.1 Fiducial case with α = 0.01 Fig. 1 shows the time evolution of gas accretion rates on to a BH (i.e. a sink cell at r = rmin) for different initial values of angular momentum: RC/RB = 0.1 (red), 0.2 (green), and 0.3 (blue). The viscous parameter is set to α = 0.01. We normalize the simulation time by the dynamical time-scale at the Bondi radius, $$t_{\rm dyn}(\equiv R_{\rm B}/c_\infty ) \simeq 1.3\times 10^{11}\ M_6 c_7^{-3}\ {\rm s}$$. For all the cases, the accretion rates increase with fluctuations at t ≲ 4 tdyn and saturate at almost constant values by t ≃ 5 tdyn. The behaviour of the accretion rates in all the cases is similar. In fact, the time-averaged values over 6 ≤ t/tdyn ≤ 16 are $$\dot{M}/\dot{M}_{\rm B}\simeq 7\times 10^{-3}$$, which does not depend on the initial angular momentum. The results show that gas rotation reduces the inflow rates by two orders of magnitude from the Bondi rate. Note that the suppression factor of the accretion rate is close to the viscous parameter we assume (α = 0.01). Of course in the limit of α = 0, there is no net accretion. Figure 1. View largeDownload slide Time evolution of the net accretion rate (in units of the Bondi rate) on to a BH, i.e. a sink cell for different values of angular momentum: RC/RB = 0.1 (red), 0.2 (green), and 0.3 (blue). The viscous parameter is set to α = 0.01. The time-averaged accretion rates (t > 6 tdyn) are $$\dot{M}\simeq 7\times 10^{-3}\ \dot{M}_{\rm B}$$. Figure 1. View largeDownload slide Time evolution of the net accretion rate (in units of the Bondi rate) on to a BH, i.e. a sink cell for different values of angular momentum: RC/RB = 0.1 (red), 0.2 (green), and 0.3 (blue). The viscous parameter is set to α = 0.01. The time-averaged accretion rates (t > 6 tdyn) are $$\dot{M}\simeq 7\times 10^{-3}\ \dot{M}_{\rm B}$$. Fig. 2 shows the radial structure of the angle-integrated mass inflow (red) and outflow (green) rates for the case with RC/RB = 0.1 and α = 0.01 at t ≃ 13 tdyn. Those rates are defined by   \begin{eqnarray} \dot{M}_{\rm in}(r)=2\pi r^2\int \rho \cdot {\rm min}(v_r,0)\sin \theta\, {\rm d}\theta , \end{eqnarray} (19)  \begin{eqnarray} \dot{M}_{\rm out}(r)=2\pi r^2\int \rho \cdot {\rm max}(v_r,0)\sin \theta\, {\rm d}\theta , \end{eqnarray} (20)where $$\dot{M}_{\rm in}< 0$$ and $$\dot{M}_{\rm out}> 0$$ (e.g. Stone et al. 1999). Note that $$\dot{M}_{\rm out}$$ in equation (20) does not necessarily mean the outflow rates of gas escaping to the infinity, but takes account into gas flows with vr > 0. We also show the net accretion rate defined by $$-\dot{M}\equiv -\dot{M}_{\rm in}-\dot{M}_{\rm out}$$ (blue). Inside the Bondi radius (r ≲ RB), both inflows and outflows exist and those rates, which decrease towards the centre following ∝ r, are tightly balanced. Within the centrifugal radius (r ≲ RC = 0.1 RB), the net accretion rate becomes a constant value ($${\simeq } 7\times 10^{-3}\ \dot{M}_{\rm B}$$). We note that the outflow rate dominates the inflow rate outside the Bondi radius. In Appendix, we study the time evolution of the outflowing gas and make sure that this component does not affect stationarity of the accretion flow (see in Fig. A2). Figure 2. View largeDownload slide Radial structure of the angle-integrated mass inflow and outflow rates for the case with RC/RB = 0.1 and α = 0.01 at t ≃ 13 tdyn. Each curve presents the inflow rate ($$-\dot{M}_{\rm in}$$ in equation 19; red), the outflow rate ($$\dot{M}_{\rm out}$$ in equation 20; green), and the net accretion rate ($$-\dot{M}\equiv -\dot{M}_{\rm in}-\dot{M}_{\rm out}$$; blue), respectively. Inside the Bondi radius (r ≲ RB), the inflow and outflow rate decrease towards the centre following ∝ r and are tightly balanced. Within the centrifugal radius, where a geometrically thick accretion disc forms, the net accretion rate becomes nearly constant. Figure 2. View largeDownload slide Radial structure of the angle-integrated mass inflow and outflow rates for the case with RC/RB = 0.1 and α = 0.01 at t ≃ 13 tdyn. Each curve presents the inflow rate ($$-\dot{M}_{\rm in}$$ in equation 19; red), the outflow rate ($$\dot{M}_{\rm out}$$ in equation 20; green), and the net accretion rate ($$-\dot{M}\equiv -\dot{M}_{\rm in}-\dot{M}_{\rm out}$$; blue), respectively. Inside the Bondi radius (r ≲ RB), the inflow and outflow rate decrease towards the centre following ∝ r and are tightly balanced. Within the centrifugal radius, where a geometrically thick accretion disc forms, the net accretion rate becomes nearly constant. In Fig. 3, we present the distribution of the gas density (top) and temperature (bottom) for two cases with RC/RB = 0.3 (left) and 0.1 (right). We also plot isodensity contours and the velocity vectors in the top and bottom panel, respectively. The elapsed time is set to t = 13 tdyn as shown in Fig. 2. For two cases, the density distribution around the Bondi radius is approximated by a rotational equilibrium distribution with a constant specific angular momentum, i.e. the Fishbone–Moncrief solution (see equation 17 and Fig. 5 below). Red dashed curve presents the boundary where ρ = 0 for the initial distribution. At the vicinity of the centrifugal radius (r ≲ 2 RC) and inside the zero-density boundary for the initial conditions (red dashed), the density distribution deviates from the equilibrium solution since the angular momentum begins to be transported outwards and the profile approaches j ∝ r1/2 (see Figs 6 and A4 below). The gas temperature increases towards the centre by compressional heating due to the gravity of the BH and energy dissipation due to viscosity. Gas motions in non-azimuthal directions are subsonic and form circulations. We note that some fraction of the gas is outflowing from inside the Bondi radius (thick black curve). In the polar directions, weakly collimated hot outflows are launched (see also green curve in Fig. 2 at r ≳ 2 RB). The outflow speed increases with the angular momentum of the system. Figure 3. View largeDownload slide Distribution of the gas density (top) and temperature (bottom) for different values of angular momentum with RC/RB = 0.3 (left) and 0.1 (right). Isodensity contours and velocity vectors are shown in the top and bottom panel, respectively. The elapsed time is set to t = 13 tdyn as shown in Fig. 2. In the top panel, the zero-density boundary for the initial distribution is shown by red dashed curve, inside which the density distribution deviates from the equilibrium solution. The thick black curve in the bottom panel shows the location of the Bondi radius. Figure 3. View largeDownload slide Distribution of the gas density (top) and temperature (bottom) for different values of angular momentum with RC/RB = 0.3 (left) and 0.1 (right). Isodensity contours and velocity vectors are shown in the top and bottom panel, respectively. The elapsed time is set to t = 13 tdyn as shown in Fig. 2. In the top panel, the zero-density boundary for the initial distribution is shown by red dashed curve, inside which the density distribution deviates from the equilibrium solution. The thick black curve in the bottom panel shows the location of the Bondi radius. Fig. 4 shows the density distribution at smaller scale for RC/RB = 0.3 (left) and 0.1 (right), respectively. For the two cases, a dense torus structure is formed around the central BH, where the gas motions are very subsonic even inside the Bondi radius and highly convective (see Section 3.1.1). From the torus surface with ρ ≃ 2 ρ∞, outflows are launched towards the poles. Except in the polar directions, the gas is both inflowing and outflowing, and the two rates are almost balanced. We emphasize that mirror-symmetry of the accretion flow cross the equatorial plane (θ = π/2) is broken. Although a previous study by Li et al. (2013) imposed equatorial mirror-symmetry and found that a large fraction of the gas is outflowing through the equator coherently, the equatorial outflow is an artefact of the imposed symmetry (see also Roberts et al. 2017). Figure 4. View largeDownload slide Distribution of the gas density in the inner region for different values of angular momentum with RC/RB = 0.3 (left) and 0.1 (right). Isodensity contours and velocity vectors are shown, and the elapsed time is set to t = 13 tdyn as shown in Fig. 2. Figure 4. View largeDownload slide Distribution of the gas density in the inner region for different values of angular momentum with RC/RB = 0.3 (left) and 0.1 (right). Isodensity contours and velocity vectors are shown, and the elapsed time is set to t = 13 tdyn as shown in Fig. 2. Angular momentum affects the gas density near the central hole. Namely, the absolute value is higher with lower angular momentum. This is because concentration of the density is suppressed by convective motions inside the centrifugal radius. As we will discuss in Section 3.1.1, the density profile becomes flatter (ρ ∝ r−1/2) inside the convective envelope (r ≲ 2 RC) rather than ρ ∝ r−3/2. 3.1.1 Comparison to self-similar solutions In order to understand properties of the accretion flows, we consider time-averaged profiles of physical quantities and compare to those of self-similar solutions for radiatively inefficient accretion flows. In the following, we show time-averaged values over 6 ≤ t/tdyn ≤ 16. Fig. 5 shows radial profiles of the gas density along the equator (θ = π/2) for three cases with RC/RB = 0.1 (red), 0.2 (green), and 0.3 (blue). The profiles consist of two components. In the outer region (r > 2 RC), the density profiles are similar to the equilibrium solution of ρ/ρ∞ = (1 + RB/r)3/2 (black dotted). On the other hand, the density follows ρ ∝ r−1/2 in the inner region (r < 2 RC) and is well approximated by   \begin{equation} \rho \simeq \rho _\infty \frac{(1+2\beta )^{3/2}}{2\beta }\left(\frac{r}{R_{\rm B}}\right)^{-1/2}, \end{equation} (21)(see dashed lines in Fig. 5). The slope of the density profile in this region agrees with that for a CDAF solution (ρ ∝ r−1/2; Narayan et al. 2000; Quataert & Gruzinov 2000) rather than that for an ADAF solution (ρ ∝ r−3/2; Narayan & Yi 1994, 1995). Figure 5. View largeDownload slide Radial profiles of the gas density along the equator (θ = π/2) for the cases with RC/RB = 0.1 (red) 0.2 (green) and 0.3 (blue). The density profiles are time-averaged over 6 ≤ t/tdyn ≤ 16, and consist of two components. In the outer region, they are approximated by an equilibrium solution of rotating gas, ‘Fishbone–Moncrief solution’, following ρ/ρ∞ = (1 + RB/r)3/2 (dotted curve). In the inner region, the density profiles approach ρ ∝ r−1/2 (dashed lines; equation 21), whose slope is the same as in CDAF solutions. The transition occurs at r ≃ 2 RC for each case (open circles). Figure 5. View largeDownload slide Radial profiles of the gas density along the equator (θ = π/2) for the cases with RC/RB = 0.1 (red) 0.2 (green) and 0.3 (blue). The density profiles are time-averaged over 6 ≤ t/tdyn ≤ 16, and consist of two components. In the outer region, they are approximated by an equilibrium solution of rotating gas, ‘Fishbone–Moncrief solution’, following ρ/ρ∞ = (1 + RB/r)3/2 (dotted curve). In the inner region, the density profiles approach ρ ∝ r−1/2 (dashed lines; equation 21), whose slope is the same as in CDAF solutions. The transition occurs at r ≃ 2 RC for each case (open circles). Fig. 6 shows radial profiles along the equator of (a) the sound speed, (b) the rotational velocity and (c) the Bernoulli number, which is defined by   \begin{equation} Be \equiv \frac{v^2}{2}+\frac{c_{\rm s}^2}{\gamma -1}-\frac{GM_\bullet }{r}. \end{equation} (22)Those profiles inside the centrifugal radius (r < 2 RC) are summarized as $$c_{\rm s}^2\propto r^{-1}$$, vϕ/vK ≃ 1, and $$Be/v_{\rm K}^2\simeq 0$$, where $$v_{\rm K}=\sqrt{GM_\bullet /r}$$ is the Keplerian velocity. Since gas accretion begins from the Bondi radius, the Bernoulli number is close to zero unless energy dissipation via viscosity is significant. Within the centrifugal radius (r ≲ RC), the rotational velocity dominates and is approximated as ∼vK. Therefore, the sound speed follows as $$c_{\rm s}^2\propto r^{-1}$$ as shown in panel (a). We note that the profiles of the sound speed do not depend on the angular momentum of the flow. In fact, the profiles can be fit with   \begin{equation} c_{\rm s}^2=c_{\infty }^2 \left[ 1+f(\gamma )\frac{R_{\rm B}}{r} \right], \end{equation} (23)where we adopt a functional form of f(γ) as   \begin{equation} f(\gamma )\equiv \frac{\gamma -1}{\gamma +1-(\gamma -1)/2} \end{equation} (24)Moreover, the rotational velocity inside the centrifugal radius is approximated as   \begin{equation} \frac{v_\phi }{v_{\rm K}} = g(\gamma ) \equiv \sqrt{\frac{2-(\gamma -1)}{\gamma +1-(\gamma -1)/2}}. \end{equation} (25)For γ = 1.6, f(γ) = 0.26 and g(γ) = 0.78. These functions of f and g have been calculated by Quataert & Gruzinov (2000) for marginally stable rotating-accretion flows against convection motions, i.e. CDAF solutions. Figure 6. View largeDownload slide Radial profile of the time-averaged (a) sound speed, (b) rotational velocity, and (c) Bernoulli number along the equator (θ = π/2) for three different cases with RC/RB = 0.1 (red), 0.2 (green), and 0.3 (blue). The dotted curves in the top two panels (equations 23 and 25) are those of the CDAF solutions, respectively. Figure 6. View largeDownload slide Radial profile of the time-averaged (a) sound speed, (b) rotational velocity, and (c) Bernoulli number along the equator (θ = π/2) for three different cases with RC/RB = 0.1 (red), 0.2 (green), and 0.3 (blue). The dotted curves in the top two panels (equations 23 and 25) are those of the CDAF solutions, respectively. Fig. 7 shows profiles of the density as a function of the polar angle θ, at radial positions of 0.02 ≤ r/RB ≤ 1.0. In the outer region (r ≳ 2 RC), the density profile follows an equilibrium distribution and explained by the Fishbone–Moncrief solution (short-dashed). In the inner region (r < 2 RC = 0.2 RB), the density concentration to the mid-plane is higher. The angular profiles in the inner region can be explained by long-dashed curve, ρ(θ) ∝ (sin θ)2[1/(γ − 1) − 1/2][ = (sin θ)2.33 for γ = 1.6] (Quataert & Gruzinov 2000). This result also indicates that the profiles of the accretion flow within the centrifugal radius are self-similar, and those properties approach those in CDAFs. Fig. 8 also presents angular profiles of the Bernoulli number. In the inner region (r < 2 RC), the value of Be is smaller than $${\lesssim } 0.2\ v_{\rm Kep}^2$$ except in the vicinity of the poles. We note that the value increases rapidly towards the poles, where weak outflows are launched. Figure 7. View largeDownload slide Angular profiles of the time-averaged density at radial positions of 0.02 ≤ r/RB ≤ 1.0. The results in the inner region (r < 2 RC) and outer region (r ≥ 2 RC) are shown by solid and dotted curves, respectively. Analytical angular profiles for CDAF solutions ρ(θ) ∝ (sin θ)[2/(γ − 1) − 1] (long-dashed) and for the Fishbone–Moncrief solution at r = RB (short-dashed) are shown for comparison to the numerical results. The profiles are normalized by the maximum value for each case. Figure 7. View largeDownload slide Angular profiles of the time-averaged density at radial positions of 0.02 ≤ r/RB ≤ 1.0. The results in the inner region (r < 2 RC) and outer region (r ≥ 2 RC) are shown by solid and dotted curves, respectively. Analytical angular profiles for CDAF solutions ρ(θ) ∝ (sin θ)[2/(γ − 1) − 1] (long-dashed) and for the Fishbone–Moncrief solution at r = RB (short-dashed) are shown for comparison to the numerical results. The profiles are normalized by the maximum value for each case. Figure 8. View largeDownload slide Angular profiles of the time-averaged Bernoulli number at radial positions of 0.02 ≤ r/RB ≤ 1.0. In the inner region (r < 2 RC), the value of Be is smaller than $${\sim } 0.2\ v_{\rm Kep}^2$$. Note that the value increases rapidly towards the poles, where weak outflows are launched. Figure 8. View largeDownload slide Angular profiles of the time-averaged Bernoulli number at radial positions of 0.02 ≤ r/RB ≤ 1.0. In the inner region (r < 2 RC), the value of Be is smaller than $${\sim } 0.2\ v_{\rm Kep}^2$$. Note that the value increases rapidly towards the poles, where weak outflows are launched. 3.1.2 Angular momentum and energy transport via convection motions Angular momentum transport by convection has been discussed with analytical methods and numerical simulations by previous works (e.g. Igumenshchev & Abramowicz 2000; Igumenshchev, Abramowicz & Narayan 2000; Narayan et al. 2000; Quataert & Gruzinov 2000; Igumenshchev et al. 2003). To analyse the effect, we calculate the r–ϕ component of Reynolds stress, $$\tau _{r \phi }\equiv -\rho \langle v^{\prime }_r v^{\prime }_\phi \rangle$$, where $$v^{\prime }_i\equiv v_i-\langle v_i \rangle$$ and 〈 · 〉 means the time-averaged value over 6 ≤ t/tdyn ≤ 16. In practice, the mass-weighted Reynolds stress integrated over the polar angle is calculated as   \begin{equation} \tau _{r\phi }(r) = -\int ^{\pi}_0 \langle \rho v^{\prime }_r v^{\prime }_\phi \rangle \sin \theta\, {\rm d} \theta . \end{equation} (26)In Fig. 9, we show the radial profile of the Reynolds stress normalized by $$\rho v_{\rm K}^2$$ for three cases of RC/RB = 0.1 (red), 0.2 (green), and 0.3 (blue). The Reynolds stress does not depend on the initial angular momentum, and the value is approximated as τrϕ ≃ Ar−3/2, where A is positive. Note that the r–ϕ component of the viscous stress has a negative value (i.e. σrϕr3/2 ≃ const. < 0). Therefore, the positive sign of A means that convective motions transport angular momentum inwards, while the standard α-viscosity transports it outwards. Igumenshchev et al. (2000) have conducted three-dimensional simulations relaxing the assumption of axisymmetry and found that the results are essentially similar to those obtained in two-dimensional simulations. In fact, the convective eddies are nearly axisymmetric and transport angular momentum inwards. This allows us to justify that our two-dimensional simulations can capture the important physics. Figure 9. View largeDownload slide Radial profiles of the r − ϕ component of the Reynolds stress τrϕ for different three cases with RC/RB = 0.1 (red), 0.2 (green), and 0.3 (blue). Solid and dashed curves show positive and negative values, respectively. The values are normalized by $$\rho v_{\rm K}^2$$. Since τrϕ ∝ Ar−3/2 is positive, convection motions transport angular momentum inwards within the centrifugal radius. Figure 9. View largeDownload slide Radial profiles of the r − ϕ component of the Reynolds stress τrϕ for different three cases with RC/RB = 0.1 (red), 0.2 (green), and 0.3 (blue). Solid and dashed curves show positive and negative values, respectively. The values are normalized by $$\rho v_{\rm K}^2$$. Since τrϕ ∝ Ar−3/2 is positive, convection motions transport angular momentum inwards within the centrifugal radius. Fig. 10 shows radial profiles of three kinds of torques for the case with RC/RB = 0.1; the Reynolds stress (red), the viscosity (green) and the advection (blue). Those three terms are calculated as   \begin{equation} T_{\rm Rey} \equiv - \frac{\mathrm{\partial} }{\mathrm{\partial} r } \left\langle 2\pi r^2 \Sigma v^{\prime }_r v^{\prime }_\phi \right\rangle , \end{equation} (27)  \begin{equation} T_{\rm vis} \equiv \frac{\mathrm{\partial} }{\mathrm{\partial} r} \left\langle 2\pi r^3 \nu \Sigma \frac{\mathrm{\partial} (v_\phi /r)}{\mathrm{\partial} r}\right\rangle , \end{equation} (28)  \begin{equation} T_{\rm adv} \equiv -\frac{\mathrm{\partial} }{\mathrm{\partial} r} \left[ \left\langle 2\pi r^2 \Sigma v_r \right\rangle \cdot \left\langle v_\phi \right\rangle \right], \end{equation} (29)where the surface density is given by $$\Sigma \approx \int _0^{\pi} \rho r\sin ^2 \theta {\rm d}\theta$$. This approximation is valid when the density concentrates near the equator, i.e. ρ(θ) ∝ (sin θ)2.33. The sign of the torque suggests that convection motions transport angular momentum inwards (red solid) rather than outwards (red dotted), as expected from Fig. 9. Within the centrifugal radius (r ≲ RC), the angular momentum is transported inwards via advection (Tadv > 0) and outwards via viscosity (Tvis < 0). Both of them are tightly balanced, and the residual agrees with the Reynolds torque within a factor of 2. Figure 10. View largeDownload slide Radial profiles of the time-averaged torque (over 6 ≤ t/tdyn ≤ 16) for the case with RC/RB = 0.1. Each curve shows the Reynolds torque (red), the viscous torque (green), and the advection torque (blue). The Reynolds component due to convection transports the angular momentum inwards (solid) and outwards (dotted), respectively. Within the centrifugal radius (r ≲ RC), the angular momentum transports via advection (inwards) and viscosity (outwards) are almost balanced. Figure 10. View largeDownload slide Radial profiles of the time-averaged torque (over 6 ≤ t/tdyn ≤ 16) for the case with RC/RB = 0.1. Each curve shows the Reynolds torque (red), the viscous torque (green), and the advection torque (blue). The Reynolds component due to convection transports the angular momentum inwards (solid) and outwards (dotted), respectively. Within the centrifugal radius (r ≲ RC), the angular momentum transports via advection (inwards) and viscosity (outwards) are almost balanced. Convective motions can transport thermal energy as well. In order to evaluate the convective energy flux Fconv, we consider the time-averaged energy equation as2  \begin{equation} \frac{1}{r^2}\frac{\mathrm{\partial} }{\mathrm{\partial} r}(r^2F_{\rm conv}) = \left\langle \sigma _{r\phi }r\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left(\frac{v_\phi }{r}\right) \right\rangle - \left\langle \rho Tv_r\frac{{\rm d}s}{{\rm d}r} \right\rangle . \end{equation} (30)The two terms on the right-hand side present time-averaged values of the energy generation rate by viscous dissipation ($$Q^+_{\rm vis}$$) and the energy advection rate ($$Q^-_{\rm adv}$$), respectively. Contribution due to the convective energy flux can be approximated by the difference between $$\langle Q^+_{\rm vis}\rangle$$ and $$\langle Q^-_{\rm adv}\rangle$$. Integrating the energy equation over angles except the polar regions, we obtain   \begin{equation} \frac{{\rm d}L_{\rm conv}}{{\rm d}r}(r)=2\pi r^2 \int _{D_\theta } (\langle Q^+_{\rm vis}\rangle -\langle Q^-_{\rm adv}\rangle )\sin \theta\, {\rm d}\theta , \end{equation} (31)where the convective luminosity is defined as $$L_{\rm conv}=2\pi r^2 \int _{D_\theta }F_{\rm conv} \sin \theta\, {\rm d}\theta$$ and Dθ = [π/6, 5π/6]. Fig. 11 shows the radial profiles of dLconv/d ln r for different values of RC/RB (dashed). Since the derivative is positive at 0.02 ≲ r/RB ≲ 1, thermal energy is transported outwards by convection. Those values are almost constant, namely $${\rm d}L_{\rm conv}/{\rm d}\ln r\simeq (0.02{\rm -}0.04)\times \dot{M}_{\rm B}c_\infty ^2$$. We note that the energy advection dominates at the innermost region (r < 0.02 RB). This is due to our inner boundary conditions, where we do not consider energy output from the sink cell at the centre. In order to study the effect of the boundary condition, we conduct a simulation for RC/RB = 0.1 with a smaller rmin, which is set to 3 × 10−3 RB (red solid). This result clearly shows that a smaller innermost grid does not change the absolute value of dLconv/d ln r but makes the convection-dominated region larger. We integrate the value from 10−2 RB to RB and calculate the convective luminosity (black solid). Then, the convective luminosity is written as $$L_{\rm conv}\sim \ \eta _{\rm conv}\dot{M}_{\rm B}c_\infty ^2$$, where the efficiency is approximated as ηconv ≃ 0.2. Figure 11. View largeDownload slide Radial profiles of the convective luminosity, dLconv/d ln r (dashed) for different three cases with RC/RB = 0.1 (red), 0.2 (green), and 0.3 (blue). In order to study the effect of the boundary condition, we conduct a simulation for RC/RB = 0.1 with a smaller rmin = 3 × 10−3RB (red solid). The integrated convective luminosity is shown by black solid, presenting $$L_{\rm conv}\simeq 0.2\ \dot{M}_{\rm B}c_\infty ^2$$. Figure 11. View largeDownload slide Radial profiles of the convective luminosity, dLconv/d ln r (dashed) for different three cases with RC/RB = 0.1 (red), 0.2 (green), and 0.3 (blue). In order to study the effect of the boundary condition, we conduct a simulation for RC/RB = 0.1 with a smaller rmin = 3 × 10−3RB (red solid). The integrated convective luminosity is shown by black solid, presenting $$L_{\rm conv}\simeq 0.2\ \dot{M}_{\rm B}c_\infty ^2$$. 3.2 Effects of viscous parameter Next, we discuss the dependence of our results on the viscous parameter α, conducting several simulations with the same parameters as shown in Section 3.1 except varying α. Fig. 12 shows the time evolution of the accretion rate on to a BH at the centre for different values of α (0.003 ≤ α ≤ 0.1). To keep the numerical simulations stable, the viscous parameter is assumed to be 0.01 at t ≤ 2 tdyn, increase (or decrease) linearly proportional to the time until t = 4 tdyn, and keep a constant value what we consider at t > 4 tdyn. The accretion rates increase with the viscous parameter because angular momentum transport becomes more efficient. For the lowest value of α, inflows with lower angular momentum through the polar region is not negligible because inflows driven by viscosity through the torus become less efficient. For the highest value of α, the fluctuation in the accretion rate is suppressed by strong viscosity, which is consistent with that reported by previous numerical simulations by Igumenshchev & Abramowicz (1999, 2000). Figure 12. View largeDownload slide Same as Fig. 1, but for different viscous parameters: α = 0.003 (magenta), 0.01 (red), 0.03 (green), and 0.1 (blue). The angular momentum is set so that RC/RB = 0.1. The accretion rate increases with α because angular momentum transport is more efficient. Figure 12. View largeDownload slide Same as Fig. 1, but for different viscous parameters: α = 0.003 (magenta), 0.01 (red), 0.03 (green), and 0.1 (blue). The angular momentum is set so that RC/RB = 0.1. The accretion rate increases with α because angular momentum transport is more efficient. In Fig. 13, we show radial profiles of the gas density for three different values of the viscous parameter: α = 0.01 (red), 0.03 (green), and 0.1 (blue). The profiles are approximated by an equilibrium solution following ρ/ρ∞ = (1 + RB/r)3/2 (dotted curve) in the outer region and ρ ∝ r−1/2, which is consistent with that in CDAF solutions in the inner region (r < 0.2 RB). For all the cases, the overall behaviour is similar, i.e. the density profile is not dependent on the choice of viscous parameter. However, for the highest value of α( = 0.1), the density slope becomes as steep as −3/2 at the innermost region (r < 0.02 RB). In this case, convective motions are suppressed by strong viscosity and thus the energy generated by viscous heating is not transported by convection but by advection with inflows, i.e. ADAFs (Narayan & Yi 1994, 1995). Figure 13. View largeDownload slide Same as Fig. 5 but for different viscous parameters: α = 0.01 (red), 0.03 (green), and 0.1 (blue). The angular momentum is set to RC/RB = 0.1. The profiles are approximated by an equilibrium solution of ρ/ρ∞ = (1 + RB/r)3/2 (black dotted) in the outer region and by ρ ∝ r−1/2 in the inner region (r < 0.2 RB), which is consistent with that in CDAF solutions, For the highest α, the density profile approaches ρ ∝ r−3/2 in the innermost region (see the text for explanation). Figure 13. View largeDownload slide Same as Fig. 5 but for different viscous parameters: α = 0.01 (red), 0.03 (green), and 0.1 (blue). The angular momentum is set to RC/RB = 0.1. The profiles are approximated by an equilibrium solution of ρ/ρ∞ = (1 + RB/r)3/2 (black dotted) in the outer region and by ρ ∝ r−1/2 in the inner region (r < 0.2 RB), which is consistent with that in CDAF solutions, For the highest α, the density profile approaches ρ ∝ r−3/2 in the innermost region (see the text for explanation). Radial profiles of other physical quantities, e.g. cs, vϕ, and Be, do not change significantly from those with α = 0.01. However, for the highest value of α( = 0.1), the normalized rotational velocity decreases at r < 0.03 RB to the centre to vϕ/vK ≃ 0.3. Such a small ratio is found in ADAF solutions (Narayan & Yi 1995). Fig. 14 shows angular profiles of the gas density at r = 0.02 RB for different viscous parameters: α = 0.01 (red), 0.03 (green), and 0.1 (blue). For lower values of α( < 0.1), the angular profiles of the density follow that of a CDAF solution (black dashed). On the other hand, for the highest α( = 0.1), where the radial density profile is as steep as −3/2, the angular profile is also different from the CDAF one. Such non-equatorial symmetric profiles produced by outflows towards one of the poles have been reported in previous numerical simulations by Igumenshchev & Abramowicz (2000) (their Model G, where α = 0.1). Figure 14. View largeDownload slide Same as Fig. 7 at r = 0.02 RB but for different viscous parameters: α = 0.01 (red), 0.03 (green), and 0.1 (blue). The angular momentum is set to RC/RB = 0.1. Dashed curve presents an angular profile expected in CDAF solutions: ρ(θ) ∝ (sin θ)[2/(γ − 1) − 1]. Figure 14. View largeDownload slide Same as Fig. 7 at r = 0.02 RB but for different viscous parameters: α = 0.01 (red), 0.03 (green), and 0.1 (blue). The angular momentum is set to RC/RB = 0.1. Dashed curve presents an angular profile expected in CDAF solutions: ρ(θ) ∝ (sin θ)[2/(γ − 1) − 1]. Fig. 15 presents the dependence of the net accretion rate on the viscous parameter α. Each symbol shows the total net accretion rates (blue cross) and the rate through the equatorial region of π/6 ≤ θ ≤ 5π/6 (red square). The latter can clarify the effect of angular momentum transport via viscosity in the torus. Those results for α ≤ 0.2 can be fit by $$\dot{M}/\dot{M}_{\rm B}\simeq 0.12\ \alpha ^\delta$$, where δ = 0.62 ± 0.24 (dashed). We note that the relation is qualitatively different from that considered in previous works; self-similar solutions or spherical accretion solutions without treating convection ($$\dot{M}\simeq \alpha \dot{M}_{\rm B}$$; Narayan & Fabian 2011). Therefore, this relation between $$\dot{M}$$ and α is an essential result from the self-consistent solution, connecting the Fishbone–Moncrief quasi-static solution (large scales) to existing CDAF solutions (small scales), and show the importance of the global solution because without our simulates, it was unclear how physical quantities on larger scales determine the properties of the accretion flow on small scales. Figure 15. View largeDownload slide Dependence of the net gas inflow rate on the viscous parameter α. The net rate increases with the viscous parameter because angular momentum is transported more efficiently. The total net accretion rates (blue cross) and the inflow rate through the equatorial region (red square) are shown. The results for α ≤ 0.2 can be fit by $$\dot{M}/\dot{M}_{\rm B}\simeq 0.12\ \alpha ^{0.62}$$ (dashed). Figure 15. View largeDownload slide Dependence of the net gas inflow rate on the viscous parameter α. The net rate increases with the viscous parameter because angular momentum is transported more efficiently. The total net accretion rates (blue cross) and the inflow rate through the equatorial region (red square) are shown. The results for α ≤ 0.2 can be fit by $$\dot{M}/\dot{M}_{\rm B}\simeq 0.12\ \alpha ^{0.62}$$ (dashed). 4 GLOBAL SOLUTIONS OF RADIATIVELY INEFFICIENT ACCRETION FLOWS We now briefly discuss the global solution of radiatively inefficient rotating-accretion flows, extending our results down to the central BH. We also give an analytical expression for the net BH feeding rate in the global solution, which can be compared to observations of low-luminosity BHs such as Sgr A* and the BH in M87. 4.1 BH feeding rate and energy loss via radiation As shown in Fig. 2, the mass inflow rate decreases towards the centre as $$\dot{M}\simeq \rho |v_r| r^2 \propto r$$. This is because the density profile of a CDAF solution follows ρ ∝ r−1/2 and the radial velocity follows |vr| ∼ ν/r ∝ r−1/2 in the α-viscosity model, respectively. Then, the net accretion rate (including both inflows and outflows) is independent of radius and at a much lower rate than the Bondi value. Because of limitation of computational time, we do not extend our computational domain down to the central BH (r ∼ RSch). Instead, we perform two additional simulations with different locations of the innermost grid of rmin = 3 × 10−3 RB (red) and 1.3 × 10−2 RB (blue). Fig. 16 shows radial profiles of the angle-integrated inflow rate (dashed) and the net accretion rate (solid). As shown clearly, the net accretion rate decreases with rmin. Therefore, the time-averaged value of the net accretion rate can be approximated as   \begin{equation} \dot{M}\simeq \left(\frac{\alpha }{0.01}\right)^\delta \left(\frac{r_{\rm min}}{R_{\rm B}}\right) \dot{M}_{\rm B}, \end{equation} (32)where δ ≃ 0.62. Assuming that the net accretion rate is constant at r ≲ 2 RC, we can estimate the net radial velocity as $$\overline{v_r} \equiv -\dot{M}/4\pi \rho r^2$$, or for β( = RC/RB) ≪ 1   \begin{equation} \overline{v_r} \simeq - \frac{\beta }{2} \left(\frac{\alpha }{0.01}\right)^\delta \left(\frac{r_{\rm min}}{R_{\rm B}}\right) \left(\frac{r}{R_{\rm B}}\right)^{-3/2} \ c_\infty . \end{equation} (33)Using the net inflow velocity, the dynamical time-scale is estimated as $$t_{\rm dyn}=r/\overline{v_r}\propto r^{5/2}$$. Inside the torus, heating via viscous dissipation dominates and the heating time-scale is given by tvis ≃ 1/[γ(γ − 1)αΩ]−1 ∝ r3/2. The ratio of the two time-scales is given by   \begin{equation} \frac{t_{\rm vis}}{t_{\rm dyn}}\sim 0.6 \left(\frac{\alpha }{0.01}\right)^{\delta -1} \left(\frac{r_{\rm min}}{10^{-2}\ R_{\rm B}}\right) \left(\frac{r}{R_{\rm C}}\right)^{-1}. \end{equation} (34)Thus, we find tvis ≲ tdyn except at r ≲ 6(β/0.1) rmin. Figure 16. View largeDownload slide Radial structure of the angle-integrated mass inflow and outflow rates for different sizes of the innermost grid rmin. Physical parameters are set to RC/RB = 0.1 and α = 0.01, and the elapsed time is t ≃ 13 tdyn. Solid and dashed curves present the inflow rate and the net accretion rate (including both inflows and outflows), respectively. Since the density follows ρ ∝ r−1/2 and the radial velocity follows |vr| ∼ ν/r ∝ r−1/2, the inflow rate is proportional to the radius ($$\dot{M}\propto r$$) within the centrifugal radius. The net accretion rate decreases with rmin and can be approximated by equation (32). The dotted line shows $$\dot{M}=(r_{\rm min}/R_{\rm B})\dot{M}_{\rm B}$$. Figure 16. View largeDownload slide Radial structure of the angle-integrated mass inflow and outflow rates for different sizes of the innermost grid rmin. Physical parameters are set to RC/RB = 0.1 and α = 0.01, and the elapsed time is t ≃ 13 tdyn. Solid and dashed curves present the inflow rate and the net accretion rate (including both inflows and outflows), respectively. Since the density follows ρ ∝ r−1/2 and the radial velocity follows |vr| ∼ ν/r ∝ r−1/2, the inflow rate is proportional to the radius ($$\dot{M}\propto r$$) within the centrifugal radius. The net accretion rate decreases with rmin and can be approximated by equation (32). The dotted line shows $$\dot{M}=(r_{\rm min}/R_{\rm B})\dot{M}_{\rm B}$$. Next, in order to evaluate the effect of radiative cooling, we compare the heating time-scale to the cooling time-scale. Since ρ ∝ r−1/2 and T ∝ r−1 inside the centrifugal radius, the bremsstrahlung cooling rate per volume is $$Q^-_{\rm br} \propto \rho ^2 T^{1/2} \propto r^{-3/2}$$ and the cooling time-scale is $$t_{\rm cool}\propto \rho T/Q^-_{\rm br} \propto r^0$$. Since tvis ≲ tdyn at r ≃ RC, thus we estimate the ratio of tvis to tcool as   \begin{equation} \frac{t_{\rm vis}}{t_{\rm cool}}\sim 0.25 \left(\frac{\alpha }{0.01}\right)^{-1} \left(\frac{r}{R_{\rm C}}\right)^{3/2} \left(\frac{T_\infty }{10^7\ {\rm K}}\right)^{1/2} \left(\frac{\dot{m}_{\rm B}}{10^{-3}}\right), \end{equation} (35)where γ = 1.6 and β = 0.1 are assumed. Bremsstrahlung cooling does not play an important role in the accretion flow as long as $$\dot{m}_{\rm B}<4\times 10^{-3}$$, which is consistent with the results shown in Li et al. (2013). Since $$Q^-_{\rm br} \propto r^{-3/2}$$, the CDAF is sandwiched between an inner region where gravitational energy would be released close to the BH and an outer radiating region (Narayan et al. 2000; Quataert & Gruzinov 2000; Ball, Narayan & Quataert 2001). 4.2 Thermal conductivity In a hot accretion flow, thermal conductivity potentially transports energy outwards and could affect the gas dynamics (e.g. Johnson & Quataert 2007; Sharma, Quataert & Stone 2008; Shcherbakov & Baganoff 2010). In the classic picture, the conductive energy flux is estimated as   \begin{equation} \boldsymbol F_{\rm cond}=-\kappa \boldsymbol \nabla T, \end{equation} (36)where κ is the conduction coefficient given by Spitzer (1962) as   \begin{equation} \kappa = 5.0\times 10^{-7} \left(\frac{\ln \Lambda _{\rm c}}{37}\right)^{-1}T^{5/2} f_{\rm c}, \end{equation} (37)in units of erg s−1 cm−1 K−1, and fc is the conductivity suppression factor because thermal conduction in the perpendicular directions to magnetic fields can be suppressed. The value of the suppression factor has been discussed by various theoretical arguments and estimated as fc ∼ 0.1 (e.g. Narayan & Medvedev 2001; Maron, Chandran & Blackman 2004). However, the precise value of the suppression factor is highly uncertain and depends on the nature of magnetic turbulence. It is known that there are two instabilities of magnetized plasma: the magneto-thermal instability (MTI; Balbus 2000) and the heat-flux-driven buoyancy instability (HBI; Quataert 2008). The MTI occurs when the temperature decreases with height and the HBI does in the opposite case. McCourt et al. (2011) have performed long-term and global MHD simulations in order to follow the evolution of both instabilities into the non-linear regime. They have found the following two results: (1) the MTI can drive strong turbulence leading to the isotropic configuration of magnetic fields (fc ∼ 0.3) and produce a large convective energy flux of $${\gtrsim}0.01 \rho c_{\rm s}^3$$ and (2) the HBI reorients the magnetic field and suppresses the conductive heat flux through the plasma. Since the value of the suppression factor is uncertain, we adopt fc = 0.1 as a fiducial one. Using our accretion solution inside the centrifugal radius, where the temperature follows T ≃ T∞f(γ)RB/r, the conductive luminosity (Lcond ≡ 4πr2Fcond) is estimated as   \begin{equation} L_{\rm cond} \simeq 6.4\times 10^{29}\ M_6 c_7^{-2} \left(\frac{f_{\rm c}}{0.1}\right) \left(\frac{r}{R_{\rm B}}\right)^{-5/2} \ {\rm erg\ s^{-1}}. \end{equation} (38)On the other hand, the convective luminosity is roughly given by   \begin{eqnarray} L_{\rm conv} &\simeq & \eta _{\rm conv} \left(\frac{\alpha }{0.01}\right)^\delta \dot{M}_{\rm B} c_\infty ^2, \nonumber \\ &\simeq & 1.1\times 10^{35}\ \rho _{-22}M_6^2 c_7^{-1} \left(\frac{\eta _{\rm conv}}{0.2}\right) \left(\frac{\alpha }{0.01}\right)^\delta \ {\rm erg\ s^{-1}}. \end{eqnarray} (39)Thus, the ratio of the two luminosities is   \begin{eqnarray} \frac{L_{\rm cond}}{L_{\rm conv}} &\simeq & 4.5\times 10^{-8} \left(\frac{f_{\rm c}}{0.1}\right) \left(\frac{\eta _{\rm conv}}{0.2}\right)^{-1} \left(\frac{\alpha }{0.01}\right)^{-\delta } \nonumber \\ &&\ \times \left(\frac{\dot{m}_{\rm B}}{10^{-3}}\right)^{-1} \left(\frac{T_\infty }{10^7\ {\rm K}}\right)^{-2} \left(\frac{r}{R_{\rm B}}\right)^{-5/2}. \end{eqnarray} (40)Therefore, we can estimate a characteristic radius where Lconv ≃ Lcond as   \begin{eqnarray} \frac{R_{\rm tr}}{R_{\rm B}} & \simeq & 1.5\times 10^{-3} \left(\frac{f_{\rm c}}{\eta _{\rm conv}}\right)^{2/5} \left(\frac{\alpha }{0.01}\right)^{-2\delta /5} \nonumber \\ && \ \times \left(\frac{\dot{m}_{\rm B}}{10^{-3}}\right)^{-2/5} \left(\frac{T_\infty }{10^7\ {\rm K}}\right)^{-4/5}. \end{eqnarray} (41)Inside this radius, thermal conduction would dominate the energy transport over convection, and the temperature profile would be flatter rather than the adiabatic scaling law (T ∝ r−1). The ratio of cs/vK would decline and the flow would resemble an optically thin viscous disc. Note that thermal conduction would hardly affect the dynamics of accretion flows at the vicinity of the BH (r < 100 RSch), where the conductive energy flux is limited to a few per cent of the saturated value of $${\sim } \rho c_{\rm s}^3$$ (Foucart et al. 2016, 2017). This fraction is consistent with that obtained in our solution at r ≳ Rtr, namely $$L_{\rm cond}/(4\pi r^2 \rho c_{\rm s}^3)\sim 1.5\times 10^{-2}(r/R_{\rm tr})^{-5/2}$$. To explore the nature of thermal conduction on the convective motions is left for future investigations. 4.3 Connection to the innermost region As discussed in Section 4.2, some physical processes, e.g. thermal conduction, can transport the energy outwards instead of convection. Once additional energy transport operates, if any, the temperature profile is no longer that for an adiabatic flow. We characterize the transition radius as Rtr = ξRB (ξ ≪ 1). The critical radius where Lconv ≃ Lcond corresponds to ξ = ξcond (≃10−3). Since the temperature follows T ∝ r−p (0 ≲ p < 1) within Rtr, the accretion disc becomes thinner and thus the density follows ρ ∝ r−3 + 3p/2, assuming that the accretion rate is constant through the disc as $$\dot{M}\simeq \xi (\alpha /0.01)^\delta \dot{M}_{\rm B}$$ (see equation 32). Inside the transition radius, the optical depth to electron scattering increases and is estimated as   \begin{eqnarray} \tau _{\rm es}(r)&=& \int _r^{R_{\rm tr}} \rho \kappa _{\rm es}\,{\rm d}r, \nonumber \\ &\simeq & \frac{40\xi ^{1/2}}{4-3p} \frac{\dot{m}_{\rm B}}{\beta } \left(\frac{c_\infty }{c}\right) \left(\frac{r}{R_{\rm tr}}\right)^{-2+3p/2}, \end{eqnarray} (42)where β ≪ 1 is assumed. We evaluate the optical depth at the innermost stable circular orbit (ISCO; r ∼ 3 RSch) for p = 0 as   \begin{eqnarray} \tau _{\rm es}^{\rm ISCO} \simeq 0.022 \left(\frac{\beta }{0.1}\right)^{-1} \left(\frac{\xi }{\xi _{\rm cond}}\right)^{5/2} \left(\frac{\dot{m}_{\rm B}}{10^{-3}}\right) \left(\frac{T_\infty }{10^7\ {\rm K}}\right)^{3/2}. \end{eqnarray} (43)Since the accretion flow at the ISCO is optically thin even in an isothermal case (p = 0), the gas is more likely to be optically thin for 0 < p < 1. In this paper, our numerical simulations do not treat thermal conductivity because the effect is subdominant in the computation domain of r ≥ rmin( > Rtr). We plan to numerically study the transition to the innermost solution in future work. 4.4 Radiation luminosities of BHs with very low accretion rates Finally, we estimate radiation luminosities produced from a low-density, radiatively inefficient accretion flow. Although our simulations do not treat the inner region (r ≲ Rtr), where energy transport by thermal conduction would dominate that by convection, we can infer the accretion rate on to the central BH from equations (32) and (41),   \begin{eqnarray} \frac{\dot{M}}{\dot{M}_{\rm Edd}}\simeq 1.5\times 10^{-6} \left(\frac{\alpha }{0.01}\right)^{3\delta /5} \left(\frac{T_\infty }{10^7\ {\rm K}}\right)^{-4/5} \left(\frac{\dot{m}_{\rm B}}{10^{-3}}\right)^{3/5}, \end{eqnarray} (44)where fc/ηconv = 1 is set. In order to estimate the radiation luminosity, we adopt a radiative-efficiency model for different values of $$\dot{M}/\dot{M}_{\rm Edd}$$, provided by axisymmetric numerical simulations of accretion flows around a supermassive BH with M• = 108 M⊙ on a small scale (r ≤ 100 RSch), including general relativistic MHD and frequency-dependent radiation transport (Ryan et al. 2017) (see also Ohsuga et al. 2009; Sa̧dowski et al. 2017). We note that their simulation results including radiation processes and electron thermodynamics, are no longer scale-free, i.e. the radiative efficiency we adopt probably depend on the choice of the BH mass. Nevertheless, it is worth demonstrating the radiation luminosity produced from a self-consistent, radiatively inefficient accretion flow, assuming that the radiative efficiency hardly depends on the BH mass. Fig. 17 shows the radiation luminosities estimated from our results (red square) as a function of the Bondi rate normalized by the Eddington luminosity. The radiation luminosity produced in the vicinity of the BH is extremely low because (1) the radiative efficiency decreases with the accretion rate in such low-accretion-rate regime ($$L\propto \dot{M}^{1.7}$$; see fig. 1 in Ryan et al. 2017) and (2) the net accretion rate at r ≲ RC is strongly suppressed by convective motions from the Bondi accretion rate $$\dot{M}_{\rm B}$$. We can approximate the radiation luminosity for $$\dot{M}/\dot{M}_{\rm Edd}\lesssim 10^{-3}$$ as   \begin{equation} \frac{L_{\rm bol}}{L_{\rm Edd}} \simeq 3\times 10^{-5} \left(\frac{\alpha }{0.01}\right)^\delta \frac{\dot{M}_{\rm B}}{\dot{M}_{\rm Edd}}, \end{equation} (45)or   \begin{equation} L_{\rm bol} \simeq 3\times 10^{-6} \dot{M}_{\rm B}c^2 \left(\frac{\alpha }{0.01}\right)^\delta . \end{equation} (46)We also present observational results of SgrA* and BHs in M31 and M87 in (blue symbols in Fig. 17). Those supermassive BHs in the nearby Universe are known as accreting BHs at low rates of $$\dot{M}_{\rm B}/\dot{M}_{\rm Edd}\sim 10^{-5}{\rm -}10^{-3}$$. Moreover, gas properties in the nuclei are studied since the angular sizes of their Bondi radii ( ≳ 1 arcsec) can be resolved well by Chandra observations (e.g. Baganoff et al. 2003; Garcia et al. 2010; Russell et al. 2015, and references therein).3 Our theoretical estimate agrees well with the observational results for SgrA* and M31. Although the observed luminosity of M87 is several times higher than our estimate, this would be because suppression of the accretion by convective motions in the CDAF regime becomes less efficient for $$\dot{M}_{\rm B}/\dot{M}_{\rm Edd}\gtrsim 10^{-3}$$, and/or because there is a level of contamination from the dominant emission due to the jet. Figure 17. View largeDownload slide Bolometric radiation luminosity produced by radiatively inefficient accretion flows (red square). The luminosities are calculated from our results given by equations (32) and (41), combining with a radiative efficiency provided by MHD simulations on a small scale (r ≤ 100 RSch), including general relativity and radiation transfer (Ryan et al. 2017). We also present observational results for SgrA* and BHs in M31 and M87 (blue asterisk). Figure 17. View largeDownload slide Bolometric radiation luminosity produced by radiatively inefficient accretion flows (red square). The luminosities are calculated from our results given by equations (32) and (41), combining with a radiative efficiency provided by MHD simulations on a small scale (r ≤ 100 RSch), including general relativity and radiation transfer (Ryan et al. 2017). We also present observational results for SgrA* and BHs in M31 and M87 (blue asterisk). Isolated BHs in our Galaxy can also accrete interstellar medium (ISM) at very low rates. Assuming typical density of the ISM and molecular clouds (n ∼ 1–100 cm−3) and the BH has a peculiar velocity of V ∼ 40–100 km s−1, the Bondi–Hoyle–Lyttleton accretion rate is as low as $$\dot{M}_{\rm B}/\dot{M}_{\rm Edd}\sim 10^{-8}{\rm -}10^{-5}$$ for M• ∼ 10 M⊙. Some of those isolated BHs in our Galaxy could be observed as high-energy sources (Fujita et al. 1998; Armitage & Natarajan 1999; Agol & Kamionkowski 2002; Ioka et al. 2017; Matsumoto, Teraki & Ioka 2018). However, since those previous studies assumed self-similar ADAF solutions, it is worth revisiting their arguments considering our self-consistent accretion solutions. Note that Perna et al. (2003) have discussed the same problem for neutron stars, and pointed out that the reduction of accretion rates from the Bondi rates can explain why isolated neutron stars accreting the ISM are rarely observed in X-rays. 5 SUMMARY We study radiatively inefficient rotating-accretion flows on to a BH with two-dimensional hydrodynamical simulations. We consider axisymmetric accretion flows and allow angular momentum transport adopting the α viscosity prescription. When the gas angular momentum is low enough to form a rotationally supported disc/torus within the Bondi radius (RB), we find a global steady accretion solution. The solution consists of three phases: (1) a rotational equilibrium distribution at r ∼ RB, where the density distribution follows ρ ∝ (1 + RB/r)3/2 and the gas is very subsonic; (2) a geometrically thick accretion torus at the centrifugal radius RC( < RB), where thermal energy generated by viscosity is transported via strong convection motions expected in CDAF; and (3) an estimated inner optically thin disc controlled by conduction and viscosity. Physical properties of the accretion flow in the intermediate region (2) agree with those in CDAFs (e.g. ρ ∝ r−1/2). In the CDAF solution, the gas accretion rate decreases towards the centre due to convection ($$\dot{M}\propto r$$), and the net accretion rate is strongly suppressed by several orders of magnitude from the Bondi accretion rate. We find that the net accretion rate depends on the viscous strength and the size of the innermost radius of the computational domain (rmin), and can be approximated as $$\dot{M}\sim (\alpha /0.01)^{0.6}(r_{\rm min}/R_{\rm B})\dot{M}_{\rm B}$$. This solution holds for low accretion rates of $$\dot{M}_{\rm B}/\dot{M}_{\rm Edd}\lesssim 10^{-3}$$ having minimal radiation cooling. In a hot plasma at the bottom of this solution (r < 10− 3 RB ≲ rmin), thermal conduction would dominate the convective energy flux. Since suppression of the accretion by convection ceases, the final BH feeding rate is found to be $$\dot{M}/\dot{M}_{\rm B} \sim 10^{-2}{\rm -}10^{-3}$$. This rate is as low as $$\dot{M}/\dot{M}_{\rm Edd} \sim 10^{-7}{\rm -}10^{-6}$$ inferred for SgrA* and the nuclear BHs in M31 and M87, and can explain the low accretion rates and low luminosities in these sources, without invoking any feedback mechanism. Finally, we briefly mention several effects which we do not treat in our simulations. First, stellar winds and X-ray emission from massive stars around the Bondi radius are neglected. These effects heat the gas to T ∼ keV and could drive outflows as in the Galactic centre around SgrA* (e.g. Najarro et al. 1997; Baganoff et al. 2003; Quataert 2004). Secondly, non-axisymmetric perturbations would lead to instability and allow additional angular momentum transport in the accretion flow (e.g. Papaloizou & Pringle 1984; Chandrasekhar 1961). In order to include these effects, we further need to conduct three-dimensional hydrodynamical simulations of accretion flows at large scales (e.g. Gaspari, Ruszkowski & Oh 2013). Thirdly, thermal conductivity in the inner region (r < Rtr ∼ 10−3 RB) plays an important role as an efficient process carrying the energy outwards instead of convection. In fact, the location of the transition radius would determine the net accretion rate on to the nuclear disc and the radiative luminosity produced in the vicinity of the BH as shown in Fig. 17. We need to investigate the detailed properties of radiatively inefficient accretion flows, in order to make a prediction for future observations with the Event Horizon Telescope.4 Acknowledgements We thank James Stone, Charles Gammie, Ramesh Narayan, Eliot Quataert, Lorenzo Sironi, Daniel Wang, Kengo Tomida, Kazumi Kashiyama, Kohei Ichikawa, Takashi Hosokawa, and Kazuyuki Sugimura for useful discussions. This work is partially supported by the Simons Foundation through the Simons Society of Fellows (KI), Simons Fellowship in Theoretical Physics (ZH), and NASA grant NNX15AB19G (ZH). 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. Numerical computations were carried out on Cray XC30 at the Center for Computational Astrophysics of the National Astronomical Observatory of Japan. Footnotes 1 For a spherically symmetric flow, the sonic point is estimated as Rsonic = [(5 − 3γ)/4]RB. Since for γ = 5/3, the sonic point is located at the origin, the inner boundary conditions would affect gas dynamics in the computational domain. Thus, we set a slightly smaller value of γ = 1.6, for which the sonic point can be resolved (Rsonic ≃ 0.05 RB). 2 The heating term due to convective viscosity and other higher order terms are neglected. Our purpose is not showing precise formulae but giving an order-of-magnitude estimate for the convective energy flux. 3 The bolometric luminosity of the M31 BH is estimated from the X-ray luminosity, assuming a bolometric correction factor of fbol = 10 for such a low-luminosity BH (Hopkins et al. 2007). 4 http://eventhorizontelescope.org REFERENCES Agol E., Kamionkowski M., 2002, MNRAS , 334, 553 CrossRef Search ADS   Armitage P. J., Natarajan P., 1999, ApJ , 523, L7 CrossRef Search ADS   Baganoff F. K. et al.  , 2003, ApJ , 591, 891 CrossRef Search ADS   Balbus S. A., 2000, ApJ , 534, 420 CrossRef Search ADS   Balbus S. A., Hawley J. F., 1991, ApJ , 376, 214 CrossRef Search ADS   Balbus S. A., Hawley J. F., 1998, Rev. Mod. Phys. , 70, 1 CrossRef Search ADS   Ball G. H., Narayan R., Quataert E., 2001, ApJ , 552, 221 CrossRef Search ADS   Bender R. et al.  , 2005, ApJ , 631, 280 CrossRef Search ADS   Blandford R. D., Begelman M. C., 1999, MNRAS , 303, L1 CrossRef Search ADS   Blandford R. D., Begelman M. C., 2004, MNRAS , 349, 68 CrossRef Search ADS   Bondi H., 1952, MNRAS , 112, 195 CrossRef Search ADS   Chandrasekhar S., 1961, Hydrodynamic and Hydromagnetic Stability , International Series of Monographs on Physics. Clarendon, Oxford Ciotti L., Ostriker J. P., 2001, ApJ , 551, 131 CrossRef Search ADS   Ciotti L., Ostriker J. P., Proga D., 2009, ApJ , 699, 89 CrossRef Search ADS   Fishbone L. G., Moncrief V., 1976, ApJ , 207, 962 CrossRef Search ADS   Foucart F., Chandra M., Gammie C. F., Quataert E., 2016, MNRAS , 456, 1332 CrossRef Search ADS   Foucart F., Chandra M., Gammie C. F., Quataert E., Tchekhovskoy A., 2017, MNRAS , 470, 2240 CrossRef Search ADS   Fujita Y., Inoue S., Nakamura T., Manmoto T., Nakamura K. E., 1998, ApJ , 495, L85 CrossRef Search ADS   Garcia M. R. et al.  , 2010, ApJ , 710, 755 CrossRef Search ADS   Gaspari M., Ruszkowski M., Oh S. P., 2013, MNRAS , 432, 3401 CrossRef Search ADS   Gebhardt K., Adams J., Richstone D., Lauer T. R., Faber S. M., Gültekin K., Murphy J., Tremaine S., 2011, ApJ , 729, 119 CrossRef Search ADS   Ghez A. M. et al.  , 2003, ApJ , 586, L127 CrossRef Search ADS   Hawley J. F., Balbus S. A., Stone J. M., 2001, ApJ , 554, L49 CrossRef Search ADS   Ho L. C., 2008, ARA&A , 46, 475 CrossRef Search ADS   Ho L. C., 2009, ApJ , 699, 626 CrossRef Search ADS   Hopkins P. F., Richards G. T., Hernquist L., 2007, ApJ , 654, 731 CrossRef Search ADS   Ichimaru S., 1977, ApJ , 214, 840 CrossRef Search ADS   Igumenshchev I. V., Abramowicz M. A., 1999, MNRAS , 303, 309 CrossRef Search ADS   Igumenshchev I. V., Abramowicz M. A., 2000, ApJS , 130, 463 CrossRef Search ADS   Igumenshchev I. V., Abramowicz M. A., Narayan R., 2000, ApJ , 537, L27 CrossRef Search ADS   Igumenshchev I. V., Narayan R., Abramowicz M. A., 2003, ApJ , 592, 1042 CrossRef Search ADS   Inayoshi K., Haiman Z., Ostriker J. P., 2016, MNRAS , 459, 3738 CrossRef Search ADS   Ioka K., Matsumoto T., Teraki Y., Kashiyama K., Murase K., 2017, MNRAS , 470, 3332 CrossRef Search ADS   Johnson B. M., Quataert E., 2007, ApJ , 660, 1273 CrossRef Search ADS   Kato S., Fukue J., Mineshige S., 2008, Black-Hole Accretion Disks – Towards a New Paradigm  –. Kyoto Univ. Press, Kyoto, Japan Kuiper R., Klahr H., Beuther H., Henning T., 2010, ApJ , 722, 1556 CrossRef Search ADS   Kuiper R., Klahr H., Beuther H., Henning T., 2011, ApJ , 732, 20 CrossRef Search ADS   Kuo C. Y. et al.  , 2014, ApJ , 783, L33 CrossRef Search ADS   Li J., Ostriker J., Sunyaev R., 2013, ApJ , 767, 105 CrossRef Search ADS   Machida M., Matsumoto R., Mineshige S., 2001, PASJ , 53, L1 CrossRef Search ADS   Maron J., Chandran B. D., Blackman E., 2004, Phys. Rev. Lett. , 92, 045001 CrossRef Search ADS PubMed  Matsumoto R., Tajima T., 1995, ApJ , 445, 767 CrossRef Search ADS   Matsumoto T., Teraki Y., Ioka K., 2018, MNRAS , 475, 1251 CrossRef Search ADS   McCourt M., Parrish I. J., Sharma P., Quataert E., 2011, MNRAS , 413, 1295 CrossRef Search ADS   McKinney J. C., Gammie C. F., 2004, ApJ , 611, 977 CrossRef Search ADS   Mignone A., Bodo G., Massaglia S., Matsakos T., Tesileanu O., Zanni C., Ferrari A., 2007, ApJS , 170, 228 CrossRef Search ADS   Milosavljević M., Bromm V., Couch S. M., Oh S. P., 2009, ApJ , 698, 766 CrossRef Search ADS   Najarro F., Krabbe A., Genzel R., Lutz D., Kudritzki R. P., Hillier D. J., 1997, A&A , 325, 700 Narayan R., Fabian A. C., 2011, MNRAS , 415, 3721 CrossRef Search ADS   Narayan R., Medvedev M. V., 2001, ApJ , 562, L129 CrossRef Search ADS   Narayan R., Yi I., 1994, ApJ , 428, L13 CrossRef Search ADS   Narayan R., Yi I., 1995, ApJ , 444, 231 CrossRef Search ADS   Narayan R., Igumenshchev I. V., Abramowicz M. A., 2000, ApJ , 539, 798 CrossRef Search ADS   Narayan R., Sa̧dowski A., Penna R. F., Kulkarni A. K., 2012, MNRAS , 426, 3241 CrossRef Search ADS   Ohsuga K., Mineshige S., Mori M., Kato Y., 2009, PASJ , 61, L7 CrossRef Search ADS   Ostriker J. P., Weaver R., Yahil A., McCray R., 1976, ApJ , 208, L61 CrossRef Search ADS   Papaloizou J. C. B., Pringle J. E., 1984, MNRAS , 208, 721 CrossRef Search ADS   Park M.-G., 1990, ApJ , 354, 64 CrossRef Search ADS   Park K., Ricotti M., 2011, ApJ , 739, 2 CrossRef Search ADS   Perna R., Narayan R., Rybicki G., Stella L., Treves A., 2003, ApJ , 594, 936 CrossRef Search ADS   Pringle J. E., 1981, ARA&A , 19, 137 CrossRef Search ADS   Proga D., 2007, ApJ , 661, 693 CrossRef Search ADS   Proga D., Begelman M. C., 2003, ApJ , 592, 767 CrossRef Search ADS   Quataert E., 2004, ApJ , 613, 322 CrossRef Search ADS   Quataert E., 2008, ApJ , 673, 758 CrossRef Search ADS   Quataert E., Gruzinov A., 2000, ApJ , 539, 809 CrossRef Search ADS   Roberts S. R., Jiang Y.-F., Wang Q. D., Ostriker J. P., 2017, MNRAS , 466, 1477 CrossRef Search ADS   Russell H. R., Fabian A. C., McNamara B. R., Broderick A. E., 2015, MNRAS , 451, 588 CrossRef Search ADS   Ryan B. R., Ressler S. M., Dolence J. C., Tchekhovskoy A., Gammie C., Quataert E., 2017, ApJ , 844, L24 CrossRef Search ADS   Sa̧dowski A., Wielgus M., Narayan R., Abarca D., McKinney J. C., Chael A., 2017, MNRAS , 466, 705 CrossRef Search ADS   Sano T., Inutsuka S.-I., Turner N. J., Stone J. M., 2004, ApJ , 605, 321 CrossRef Search ADS   Shakura N. I., Sunyaev R. A., 1973, A&A , 24, 337 Shapiro S. L., 1973, ApJ , 180, 531 CrossRef Search ADS   Sharma P., Quataert E., Stone J. M., 2008, MNRAS , 389, 1815 CrossRef Search ADS   Shcherbakov R. V., Baganoff F. K., 2010, ApJ , 716, 504 CrossRef Search ADS   Spitzer L., 1962, Physics of Fully Ionized Gases . Interscience, New York Stone J. M., Norman M. L., 1992, ApJS , 80, 753 CrossRef Search ADS   Stone J. M., Pringle J. E., 2001, MNRAS , 322, 461 CrossRef Search ADS   Stone J. M., Hawley J. F., Gammie C. F., Balbus S. A., 1996, ApJ , 463, 656 CrossRef Search ADS   Stone J. M., Pringle J. E., Begelman M. C., 1999, MNRAS , 310, 1002 CrossRef Search ADS   Walsh J. L., Barth A. J., Ho L. C., Sarzi M., 2013, ApJ , 770, 86 CrossRef Search ADS   Yuan F., Quataert E., Narayan R., 2003, ApJ , 598, 301 CrossRef Search ADS   Yuan F., Wu M., Bu D., 2012, ApJ , 761, 129 CrossRef Search ADS   APPENDIX: NUMERICAL TESTS In this appendix, we discuss the effects of numerical resolutions, outer boundary conditions, initial conditions, prescriptions of viscosity on our results. As a reference, we define our fiducial case with RC/RB = 0.1 and α = 0.01, where (1) the number of grid cells is Nr × Nθ = 512 × 256, (2) the initial condition is given by a rotating equilibrium density configuration by equation (17), and (3) the additional viscosity which would be driven by rotational instability is assumed (see equation 16). Fig. A1 shows the time evolution of the net accretion rates on to the BH (i.e. a sink cell at the centre) with different numerical resolutions: Nr × Nθ = 512 × 256 (fiducial, red), 256 × 256 (low resolution, green), and 1024 × 512 (high resolution, blue). For the three cases, the time-averaged accretion rate is $$\dot{M}/\dot{M}_{\rm B}=7.1\times 10^{-3}$$ (red), 7.9 × 10−3 (green), and 6.3 × 10−3 (blue), respectively. The estimated errors are at most ∼ 20 per cent. Figure A1. View largeDownload slide Time evolution of the net accretion on to a sink cell (in units of the Bondi rate) for RC/RB = 0.1 and α = 0.01. Simulation results with three different resolutions are shown: Nr × Nθ = 512 × 256 (fiducial, red), 256 × 256 (low resolution, green), and 1024 × 512 (high resolution, blue). Figure A1. View largeDownload slide Time evolution of the net accretion on to a sink cell (in units of the Bondi rate) for RC/RB = 0.1 and α = 0.01. Simulation results with three different resolutions are shown: Nr × Nθ = 512 × 256 (fiducial, red), 256 × 256 (low resolution, green), and 1024 × 512 (high resolution, blue). As seen in Figs 2 and 3, a significant fraction of the gas is outflowing towards the polar directions and the rate dominates the inflow rate outside the Bondi radius. In order to ensure stationarity of the accretion system, we study the time evolution of the outflow rate. Fig. A2 shows radial profiles of the outflow rates (i.e. vr > 0, see equation 20) for different elapsed times (6 ≤ t/tdyn ≤ 23). This component propagates outwards, accumulating the ambient gas which has a uniform distribution and no accretion initially at r > RB. However, the outflow rates inside the Bondi radius seems steady, and the rates even outside the Bondi radius have nearly converged (r ≲ 7 RB). Therefore, our long-term simulations allow us to ensure the stationarity of the accretion flow in the inner region. Figure A2. View largeDownload slide Time evolution of radial profiles of the outflow rates for RC/RB = 0.1 and α = 0.01. Different curves show the time evolution: t/tdyn = 6 (1, red), 10 (2, green), 14 (3, blue), 18 (4, magenta), and 23 (5, light blue). Although the outflow propagates, accumulating the ambient gas, the outflow rate inside the Bondi radius is steady. This result ensures the stationarity of the accretion flow in the inner region. Figure A2. View largeDownload slide Time evolution of radial profiles of the outflow rates for RC/RB = 0.1 and α = 0.01. Different curves show the time evolution: t/tdyn = 6 (1, red), 10 (2, green), 14 (3, blue), 18 (4, magenta), and 23 (5, light blue). Although the outflow propagates, accumulating the ambient gas, the outflow rate inside the Bondi radius is steady. This result ensures the stationarity of the accretion flow in the inner region. Fig. A3 shows the effect of initial conditions on the time evolution of the accretion rate. Red and green curve present the fiducial case and the result with a Bondi accretion solution added a constant angular momentum j∞ as the initial conditions. In the latter case, the accretion rate decreases rapidly by several orders of magnitude in the early stage due to the centrifugal force, begins to increase gradually at t ≳ 2 tdyn and approaches that in the fiducial case at t ≳ 6 tdyn. Therefore, this result clearly shows that the net accretion rate does not depend on the choice of the initial conditions. Figure A3. View largeDownload slide Time evolution of the net accretion on to a sink cell (in units of the Bondi rate) for RC/RB = 0.1 and α = 0.01. Red curve presents the fiducial case as shown in Fig. 1. Green one is the result where the simulation starts from the Bondi solution with a constant angular momentum j∞. Blue one is the result where the second term in equation (16), which provides additional viscosity led by rotational instability, is neglected. Figure A3. View largeDownload slide Time evolution of the net accretion on to a sink cell (in units of the Bondi rate) for RC/RB = 0.1 and α = 0.01. Red curve presents the fiducial case as shown in Fig. 1. Green one is the result where the simulation starts from the Bondi solution with a constant angular momentum j∞. Blue one is the result where the second term in equation (16), which provides additional viscosity led by rotational instability, is neglected. Finally, we study the effect of the additional viscosity given by the second term in equation (16), which would be led by rotational instability. Blue curve in Fig. A3 shows the accretion rate without the additional term, and almost identical to that in the fiducial case. Fig. A4 shows radial profiles of the specific angular momentum normalized by $$j_\infty (=\sqrt{\beta }R_{\rm B}c_\infty )$$ for RC/RB = 0.1 and α = 0.01. Solid (red) curve presents the result at t = 13 tdyn in our fiducial case. The profile has a Keplerian-like profile expected from a CDAF solution inside the centrifugal radius, $$j=g(\gamma )\sqrt{GM_\bullet r}$$ (dotted black), and a constant value outside the torus (r ≳ 2 RC). Two dashed curves show the cases without additional viscosity in the second term in equation (16) at two different elapsed time of t = 13 tdyn (green) and t = 16 tdyn (blue). Compared to the fiducial case, the profile at the same elapsed time (green curve) have an excess of the angular momentum from the boundary value of j∞ outside the torus. Since viscosity does not work outside the torus, the angular momentum is accumulated and the bump structure grows with time (see blue curve). As discussed in Section 2, the regions with negative gradients of the angular momentum (i.e. dj/dr < 0) are unstable and would produce turbulence, which allows additional angular momentum transport. In order to capture the physics of the instability, we need three-dimensional hydrodynamical simulations. However, we have confirmed that the treatment of the rotational instability does not affect our results. Figure A4. View largeDownload slide Radial profiles of the specific angular momentum normalized by $$j_\infty (=\sqrt{\beta }R_{\rm B}c_\infty )$$ for RC/RB = 0.1 and α = 0.01. Solid (red) curve presents the profile at t = 13 tdyn in our fiducial case, where the second term in equation (16) is included to provide additional viscosity led by rotational instability. Dashed curves show the results at t = 13 tdyn (green) and 16 tdyn (blue), respectively, in cases where the additional viscosity is neglected. Dotted (black) curve is that expected in a CDAF solution, $$j=g(\gamma )\sqrt{GM_\bullet r}$$. Figure A4. View largeDownload slide Radial profiles of the specific angular momentum normalized by $$j_\infty (=\sqrt{\beta }R_{\rm B}c_\infty )$$ for RC/RB = 0.1 and α = 0.01. Solid (red) curve presents the profile at t = 13 tdyn in our fiducial case, where the second term in equation (16) is included to provide additional viscosity led by rotational instability. Dashed curves show the results at t = 13 tdyn (green) and 16 tdyn (blue), respectively, in cases where the additional viscosity is neglected. Dotted (black) curve is that expected in a CDAF solution, $$j=g(\gamma )\sqrt{GM_\bullet r}$$. © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society

Journal

Monthly Notices of the Royal Astronomical SocietyOxford University Press

Published: May 1, 2018

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off