© The Authors 2026

1. Introduction

Supermassive black holes (SMBHs) are widely believed to exist at the centers of galaxies. In nearby galaxies, most of active galactic nuclei (AGNs) are low-luminosity AGNs (LLAGNs) with luminosities ranging from 10⁻⁹ to 10⁻¹ L_Edd, where L_Edd is the Eddington luminosity (Ho 2008, 2009). Hot accretion flows around SMBHs are generally employed to explain the observed properties of LLAGNs. Hot accretion flows are characterized by lower density and higher temperature, which causes them to be radiatively inefficient. The energy dissipated in hot accretion flows is advected into SMBHs along with the accreting matter. These flows are commonly referred to as advection-dominated accretion flows (ADAFs, Ichimaru 1977; Rees et al. 1982; Narayan & Yi 1994, 1995; Abramowicz et al. 1995).

In recent years, a series of numerical studies have investigated the hot accretion flows in detail (Stone et al. 1999; Igumenshchev & Abramowicz 1999, 2000; Hawley et al. 2001; Stone & Pringle 2001; Hawley & Balbus 2002; Yuan et al. 2012). One of the significant findings is the identification of strong winds launched from the radiatively-inefficient hot accretion flows (Yuan et al. 2012, 2015). Stone et al. (1999) implemented the two-dimensional hydrodynamic (HD) simulations of non-radiative hot accretion flows and found that the mass inflow rate diminishes as the radius decreases. Yuan et al. (2012) confirmed that the wind-induced mass loss results in a reduction of the mass inflow rate as the accretion gas moves inward. Based on observations from the Chandra X-ray Observatory, Wang et al. (2013) confirmed the existence of winds originating from hot accretion flows surrounding the SMBH in our galactic center (Sgr A^*).

For further investigating the properties of winds, Yuan et al. (2015) performed the trajectory analysis of the virtual Lagrangian particles based on a 3D general relativistic magnetohydrodynamic (MHD) simulation of non-radiative hot accretion flows around a Schwarzschild BH. They found that, at small scale (r < 10³r_s, with r_s being Schwarzschild radius), the radial profile of the mass flux of winds can be approximately expressed as ṁ_wind ≈ ṁ_BH (r/20r_s), where ṁ_BH is the mass accretion rate at the BH horizon. This result suggests that most of the winds are produced at large radii within the region of r < 10³r_s. When radiative cooling is considered, the strength of winds depends on the accretion rate (Bu & Gan 2018). Based on the 2D HD simulations, Bu & Gan (2018) found that radiative cooling becomes important with the increase of accretion rate. When the accretion rate increases to a value, winds become very weak. Their study did not consider the effect of magnetic field. However, if magnetic field is taken into account, the wind strength may be changed due to the presence of the magnetic-driving winds (Blandford & Payne 1982; Cao 2011; Li & Begelman 2014). The aforementioned studies primarily investigate the generation of winds on small scales (r < 10³r_s), i.e. on the region close to the BH. Therefore, the important questions are whether or not the winds can be produced at large scale (r >  10³r_s) and how the accretion rate influences the wind strength.

The accretion processes around a BH are generally believed to extend out to the Bondi radius (Bondi 1952). For the spherically symmetric flow without rotation and radiation, within the Bondi radius, the gas is gravitationally bound. For a realistic accretion flow, due to the presence of gas angular momentum and radiative cooling, the accretion flow can not be adequately described by the Bondi solution. A dynamical range extending from the Bondi radius to the vicinity of a BH is excessively large. It is difficult to implement numerical simulation on this dynamical range. A feasible approach is to independently implement numerical simulations at small scales (the vicinity of a BH) and large scales (near the Bondi radius) to investigate the dynamics of accretion flows. Several previous works have implemented HD and MHD simulations to study the hot accretion flows at large scales (Li et al. 2013; Bu et al. 2016a,b; Inayoshi et al. 2018, 2019; Mosallanezhad et al. 2022). Bu et al. (2016a,b) and Inayoshi et al. (2018) studied the radiatively-inefficient accretion flows around the sub-Bondi radius. Li et al. (2013) implemented HD simulations to study the effect of radiative cooling on the accretion flows around the sub-Bondi radius. They found that strong radiative cooling promotes accretion. In the HD simulations, the α-viscosity prescription is usually employed to mimic the anomalous shear stresses for angular momentum transport and the heating effects resulting from the work done by stresses. However, HD simulation results depends on the assumption of α. In reality, the transport of angular momentum in accretion flows should be driven by Maxwell stresses arising from MHD turbulence, which is induced by the magnetorotational instability (MRI; Balbus & Hawley 1998). Therefore, it is necessary to perform MHD simulations on the accretion flows around the Bondi radius.

In this paper, we perform MHD simulations to investigate the effects of radiative cooling on the dynamics and wind production in hot accretion flows around the sub-Bondi radius. In Section 2, we describe the method and model setup. In Section 3, we analyze the main results. In Sections 4 and 5, we present a discussion and a summary of the results, respectively.

2. Method

2.1. Basic equations and numerical method

To describe the BH hot accretion flows around the sub-Bondi radius, we employed the MHD equations as follows:

$$\begin{array}{cc}& \frac{\partial \rho }{\partial t}+\mathrm{\nabla }·(\rho \mathit{v})=0,\hfill \end{array}$$(1)

$$\begin{array}{cc}& \frac{\partial (\rho \mathit{v})}{\partial t}+\mathrm{\nabla }·[\rho \mathit{v}\mathit{v}-\mathit{B}\mathit{B}+(p_{\mathrm{g}}+\frac{{\mathit{B}}^2}{2})]=-\rho \mathrm{\nabla }\mathrm{\Psi },\hfill \end{array}$$(2)

$$\begin{array}{cc}& \frac{\partial E_{\mathrm{t}}}{\partial t}+\mathrm{\nabla }·[(E_{\mathrm{t}}+p_{\mathrm{g}})\mathit{v}-\mathit{B}(\mathit{v}·\mathit{B})]=-\rho \mathit{v}·\mathrm{\nabla }\mathrm{\Psi }-Q_{\mathrm{brem}}^{-},\hfill \end{array}$$(3)

$$\begin{array}{cc}& \frac{\partial \mathit{B}}{\partial t}-\mathrm{\nabla }\times (\mathit{v}\times \mathit{B})=0.\hfill \end{array}$$(4)

Here, ρ, p_g, v, and B¹ are the density, pressure, velocity, and magnetic field, respectively. E_t is the total energy density, E_t = e + ρv²/2 + B²/2, where e is the internal energy. We adopted an equation of state of ideal gases, p_g ≡ (γ − 1)e ≡ ρκ_BT_g/μm_p, and set γ = 5/3, where μ, κ_B, m_p, and T_g are the mean molecular weight, the Boltzmann constant, the proton mass, and the gas temperature, respectively. In this paper, we set μ = 0.6. The gravitational potential (Ψ) of the center BH is described using the pseudo-Newtonian potential introduced by Paczyńsky & Wiita (1980), i.e. Ψ(r) = − GM/(r − r_s), where G is the gravitational constant, M is the mass of the central BH, and r_s ≡ 2GM/c² is the Schwarzschild radius. $Q_{\mathrm{brem}}^{-}$ is the bremsstrahlung cooling term given by Li et al. (2013):

$$\begin{array}{c}\hfill Q_{\mathrm{brem}}^{-}={\alpha }_f{r}_e^2{m}_e{c}^3{n}^2(32/3){(2/\pi )}^{1/2}{\left(\frac{{\kappa }_{\mathrm{B}}{T}_{\mathrm{g}}}{m_e{c}^2}\right)}^{1/2},\end{array}$$(5)

where α_f is the fine structure constant, r_e is the classical electron radius, c is the speed of light, and n is the number density.

We set the BH mass to be M = 10⁸ M_⊙. Assuming the axisymmetric condition, i.e ∂/∂ϕ = 0, we employed the PLUTO code to numerically solve Equations (1)–(4) in a spherical coordinates (r, θ, ϕ). The PLUTO code is a Godunov-type code solving the HD and MHD equations (Mignone et al. 2007, 2012). The bremsstrahlung cooling term in equation (3) was considered as a source term, which was solved implicitly to update the gas temperature and internal energy in the PLUTO code.

The computation domain is from r_in = 2 × 10⁻²r_B to r_out = 5r_B in radial direction and from θ₁ = 0 to θ₂ = π in angular direction. Here, r_B is the Bondi radius, which is determined based on the gas temperature at infinity T_∞. We set T_∞ ∼ 1.7 × 10⁶ K, and then r_B equals 10⁶r_s. In the radial direction, we employ 256 zones with the radial size ratio of (△r)_i + 1/(△r)_i = 1.0177. In the angular direction, 104 zones with the angular size ratio of (△θ)_j + 1/(△θ)_j = 0.98106 are distributed over 0 < θ < π/2 and 104 zones with the angular size ratio of (△θ)_j + 1/(△θ)_j = 1.01930 are distributed over π/2 < θ < π. The smallest angular size is △θ = 0.00480956 at θ = π/2. Axisymmetric boundary conditions are employed at the poles and outflow boundary conditions are employed at the inner and outer radial boundaries.

2.2. Initial condition

Initially, we assumed that a torus in rotating equilibrium with constant specific angular momentum is embedded in a non-rotating, low-density medium. We followed Kato et al. (2004) to determine properties of the torus. The specific angular momentum of the torus was set to be l = (GMr₀³)^1/2/(r₀ − r_s), where r₀ (=r_B) is the radius of the torus center. For the initial torus, a polytropic relation between pressure (p_t) and density (ρ_t) is $p_{\mathrm{t}}=K{\rho }_{\mathrm{t}}^{1+1/n}$, where K is a constant and n was set to be equal to 3. The initial density and pressure distributions of the torus are given by

$$\begin{array}{c}\hfill {\rho }_t={\rho }_0{(1-\frac{\gamma }{c_{\mathrm{s},0}^2}\frac{\stackrel{\sim }{\mathrm{\Psi }}-{\stackrel{\sim }{\mathrm{\Psi }}}_0}{n+1})}^n,\end{array}$$(6)

and

$$\begin{array}{c}\hfill p_t={\rho }_0\frac{c_{\mathrm{s},0}^2}{\gamma }{\left(\frac{{\rho }_t}{{\rho }_0}\right)}^{1+\frac{1}{n}},\end{array}$$(7)

where c_s, 0² is the sound speed at the torus center, $\stackrel{\sim }{\mathrm{\Psi }}$ (≡Ψ + (l/rsinθ)²/2) is the effective potential, and ${\stackrel{\sim }{\mathrm{\Psi }}}_0$ is the effective potential at torus center, respectively. We can give the gas temperature (T_g, 0) at the torus center to determine c_s, 0² (≡γk_BT₀/μm_μ) and then determine the initial density and pressure of the torus. Here, we set T_g, 0 = 2 × 10⁵ K. For the initial ambient medium, we assume the initial density ${\rho }_a={10}^{-6}{\rho }_0{(\frac{r_{\mathrm{in}}}{r})}^4$ and the initial pressure $p_a=\frac{GM{\rho }_a}{r}$. The density floor is also set to be equal to ρ_a.

2.3. Magnetic field setup

The magnetic field topology at the Bondi radius is observationally unknown. According to Han (2017), although the magnetic field may exhibit disorder on small scales, it is expected to be ordered on galactic scales. The introduced weak magnetic field in numerical simulations can grow under MRI perturbations. Beckwith et al. (2008) investigated the impact of the initial magnetic field topology on accretion flows. They employed three different topological configurations, such as dipole, quadrupole, and toroidal fields. Their results indicate that once the accretion flow reaches a quasi-steady state, the initial topology of the magnetic field does not significantly affect its primary characteristics. Yuan & Narayan (2014) also pointed out that the initial topology of the magnetic field exerts minimal influence on the characteristics of the accretion flows once they have developed into a quasi-steady state. However, it may have a certain impact on the properties of the jet. Generally, the dipole field generates a strong and steady jet, whereas the quadrupole field produces a weaker and intermittent jet. In contrast, the toroidal field yields almost no jet at all. Here, we initially chose a vector potential to introduce a weak seed magnetic field.

Similar to the magnetic field configuration proposed by Penna et al. (2013), we constructed the initial magnetic fields that are purely poloidal, i.e. only the toroidal component (A_ϕ) of magnetic potential is nonzero. The toroidal component is then

$$\begin{array}{c}\hfill {\mathit{A}}_{\mathit{\varphi }}\propto \{\begin{array}{ccc}\hfill q\sin \left(\frac{N_{\mathrm{L}}\pi (g(r)-g(r_{\mathrm{beg}})}{g(r_{\mathrm{end}})-g(r_{\mathrm{beg}})}\right)& & q>0.0\hfill \\ \hfill 0.0& & \mathrm{otherelse}\hfill \end{array},\end{array}$$(8)

where

$$\begin{array}{c}\hfill \{\begin{array}{cc}\hfill g(r)& =r^{2/3}+\frac{15}{8}{r}^{-2/5}\hfill \\ \hfill q& ={(cos\theta )}^{N_{\mathrm{q}}}{\sin }^3\theta max(\rho -{\rho }_{\mathrm{cut}},0)\hfill \end{array}.\end{array}$$(9)

Here, r_beg and r_end are the inner and outer boundaries of the magnetized region, respectively. ρ_cut is the cut-off density and we set ρ_cut = 10⁻⁴ρ₀. N_L is a parameter to control the number of the poloidal magnetic loops along the r direction. N_q is a parameter that determines whether a magnetic loop is divided into two mirror-symmetric loops relative to the equator or not. When N_q = 0, only one loop is distributed along the θ direction. When N_q = 1, a loop is divided into two loops that are mirror symmetric relative to the equator.

For each magnetic loop, we specify the scale factor in Equation (8) to determine the magnetic field strength. The scale factor is determined by the ratio (β₀ = p_g/p_m) of the gas pressure to the magnetic pressure (p_m) at the location of maximum magnetic pressure over each loop. In our runs, we set N_L = 4, N_q = 1 and β₀ = 1000. The configuration of magnetic fields are shown in Figure 1.

Fig. 1. Distribution of the initial magnetic field (lines) and plasma β0 (colors). Solid and dashed lines represent clockwise and counterclockwise magnetic field lines, respectively.

3. Result

Table 1 summarized our simulation runs. In our simulations, we can adjust the accretion rate by varying the density (ρ₀) at the center of the initial torus. As shown in Table 1, the density increases from 10⁻²⁴ g cm⁻³ in run R24 to 10⁻²² g cm⁻³ in run R22. For the purpose of comparison, we do not incorporate radiative cooling in the simulation of run N24. In the simulations of runs R22–R24, we initially let the runs without radiative cooling evolve from the initial state to a quasi-steady state. Subsequently, we incorporate radiative cooling into these runs and then permit them to evolve continuously to achieve a quasi-steady state once again. In general, we incorporate radiative cooling into these runs after t = 2t_B, where t_B is the Keplerian orbit period at r = r_B. In the following analysis, we perform time-averaging over the time interval of 4−5 t_B. During this time interval, we gather 400 data files with equal time intervals.

In magnetized accretion flows, the MRI is the primary mechanism for angular momentum transport. To ensure that the MRI is adequately captured in numerical simulations, it is essential to resolve its fastest growing MRI mode. The characteristic wavelength of the fastest growing MRI mode is given by λ_MRI = 2πV_A/Ω, where V_A is the Alfvén speed and Ω is the angular velocity. Following Hawley et al. (2011, 2013), we define the MRI quality factors in the radial and polar directions as Q_r = 2πV_A, r/Ω Δr and Q_θ = 2πV_A, θ/Ω rΔθ, where V_A, r and V_A, θ are the radial and polar components of the Alfvén velocity, respectively, and Δr and rΔθ represent the local grid size in the corresponding directions. These quality factors quantify the number of grid cells available to resolve the fastest growing MRI mode. Previous studies have suggested that maintaining Q_r ≳ 10 and Q_θ ≳ 10 is sufficient to sustain MRI-driven turbulence throughout the nonlinear regime (Hawley et al. 2011, 2013). In our simulations, within the region that achieves a quasi-steady state, the MRI quality factors exceed the threshold value of 10 throughout the majority of the domain. This demonstrates that our simulations are sufficiently well-resolved to effectively capture MRI turbulence, thereby ensuring the reliability of the derived accretion dynamics and wind properties.

In all the runs, accretion flows achieve a quasi-steady state within the region of r < ∼0.3r_B, where the net accretion rate remains approximately constant (see the dotted line in Figure 11). Therefore, in the following analysis, we focus on the accretion flows within the quasi-steady state region. For the purpose of analysis, we employed two distinct methods for averaging physical quantities. The time-averaged value of a given quantity (a) in each grid (r, θ) is defined as follows:

$$\begin{array}{c}\hfill {⟨a⟩}_{\mathrm{t}}=\frac{{\int }_{t_1}^{t_2}adt}{{\int }_{t_1}^{t_2}dt}·\end{array}$$(10)

The time- and angular averaged value of a given quantity (a) at each radius is defined as follows:

$$\begin{array}{c}\hfill ⟨⟨a⟩⟩=\frac{{\int }_{t_1}^{t_2}{\int }_{{\theta }_1}^{{\theta }_2}asin\theta d\theta dt}{{\int }_{t_1}^{t_2}{\int }_{{\theta }_1}^{{\theta }_2}sin\theta d\theta dt}·\end{array}$$(11)

Here, the time interval spans from t₁ = 4t_B to t₂ = 5t_B. For evaluating the angular averaged value around the equatorial plane, the integration limits for the polar angle were defined as θ₁ = 86° and θ₂ = 94°.

3.1. Structure of accretion flows

Figure 2 shows the distribution of density and velocity for runs R22–R24 and N24. The black solid line in Figure 2 marks the density scale height (H) of accretion flows. The density scale height is defined as the vertical distance at which the density drops by a factor of e relative to the density at the equatorial plane. Here, we define the accretion disk as the accretion flows within the density scale height. From Figure 2, we can see that outflows are generated as matter escapes from the accretion disk. In order to compare the density scale height of various runs, we plot Figure 3, which shows the radial dependence of the ratio (H/r) of density scale height to radius. In run R24 with radiative cooling, the density scale height is comparable to that observed in run N24 without radiative cooling. Both runs exhibit similar values for H/r, ranging from approximately 0.3 to 0.6. This indicates that radiative cooling is very weak in R24. In run R23, H/r is slightly lower, approximately 0.2 − 0.5, due to radiative cooling. In run R22, the density scale height is the lowest among all the runs, and beyond 0.1 r_B the H/r value is lower than 0.2, which implies that the accretion flows undergo a significant contraction due to strong radiative cooling. Therefore, as the accretion rate increases from run R24 to run R22, the radiative cooling gradually strengthens, resulting in the accretion flows gradually contracting.

Fig. 2. Distribution of density and velocity field. In each panel, the left side presents a snapshot at 4.0tB, while the right side displays the time-averaged distribution. Colors represent the density, measured in units of the density at the center of the initial torus. The streamlines reflect the velocity field on the poloidal plane. The black solid line marks the scale height of accretion flows. The scale height is defined as the height at which the density drops by a factor of e relative to the density at the equatorial plane. Panel (A): Run R22. Panel (B): Run R23. Panel (B): Run R24. Panel (D): Run N24.

Fig. 3. Radial dependence of the ratio (H/r) of density scale height to radius. The black, blue, red, and green lines correspond to runs R22, R23, R24, and N24, respectively.

Figure 4 shows the radial profiles of the time-averaged and angle-averaged density near the equatorial plane. Within r = ∼0.07r_B, we can employ a power-law function ρ ∝ r^−p to describe the radial profiles of the density near the equatorial plane. In Table 1, we listed the values of the power-law index p in a function of ρ ∝ r^−p for all the runs. As radiative cooling strengthens, the value of p decreases from 0.37 in run R24 to 0.05 in run R22. For run N24 without radiative cooling, p = 0.50. This indicates that radiative cooling leads to a flatter density distribution in the inner region of the accretion flows. Outside 0.07r_B, the density distribution exhibits different behaviors depending on the specific run. In runs N24 and R24, the density slightly decreases as the flows move inward from ∼0.3r_B to ∼0.2r_B. In the range where r ∼ 0.07 − 0.2r_B, a slight increase inward in density is observed. In runs R23 and R22, the density rapidly decreases inward from 0.3r_B to ∼0.07r_B. Notably, the decrease inward in density is more pronounced in run R22, which is characterized by stronger cooling effects. This implies that radiative cooling reduces the efficiency of material transfer from a larger radius (∼0.3r_B) to a smaller radius (0.1r_B) becomes low. We will discuss the reason for this phenomenon in Section 3.5.

Fig. 4. Radial profiles of the time-averaged and angle-averaged density near the equatorial plane. The angle averaging is implemented over the angle range of 84° to 96°. The black, blue, green, and red lines correspond to runs R22, R23, R24, and N24, respectively.

Figure 5 shows the time-averaged distribution of the plasma-β value and the poloidal mean magnetic field for runs R22 to R24 and N24. The time-averaged plasma-β² is defined as the ratio of the time-averaged gas pressure to the magnetic pressure of the time-averaged magnetic field, i.e., β = 2⟨p_g⟩_t/⟨B⟩_t², where ${⟨\mathit{B}⟩}_{\mathrm{t}}\equiv {⟨B_{\mathrm{r}}⟩}_{\mathrm{t}}\widehat{\mathit{r}}+{⟨B_{\theta }⟩}_{\mathrm{t}}\widehat{\mathit{\theta }}+{⟨B_{\varphi }⟩}_{\mathrm{t}}\widehat{\mathit{\varphi }}$ is the time-averaged magnetic field. In run N24 without radiative cooling, the gas pressure around the equatorial plane is stronger than the magnetic pressure, resulting in a relatively weak magnetic field in this region. In runs R24 and R23, which are characterized by weak radiative cooling, similar results are also observed. In contrast, run R22 with strong radiative cooling exhibits a relatively stronger magnetic field around the equatorial plane. In the region near the equatorial plane beyond ∼2r_B, the magnetic pressure even exceeds the gas pressure. To further evaluate the relative strengths of the toroidal and poloidal magnetic fields, Figure 6 presents the distribution of the strength ratio between these two types of magnetic fields (⟨B_ϕ⟩_t²/(⟨B_r⟩_t² + ⟨B_θ⟩_t²)). In all runs, except for the polar region where the poloidal magnetic field is predominant, the toroidal magnetic field dominates in most regions. Around the equatorial plane, the toroidal magnetic energy accounts for 60% to 90% of the total magnetic energy (refer to the dashed line in Figure 20).

Fig. 5. Distribution of the plasma β and poloidal magnetic fields. Colors represent the ratio (β) of the time-averaged gas pressure to the magnetic pressure of the time-averaged magnetic field. The solid and dashed lines reflect the magnetic fields on the poloidal plane. The purple line indicates the interface at which β = 1. Panel (A): Runs R22 (left) and R23 (right). Panel (B): Runs R24 (left) and N24 (right).

Fig. 6. Distribution of the strength ratio between the toroidal and poloidal components of the mean magnetic field. The purple line indicates the interface where the toroidal magnetic field strength equals the poloidal magnetic field strength. Panel (A): Runs R22 (left) and R23 (right). Panel (B): Runs R24 (left) and N24 (right).

In the following, we further analyze the reasons for the strengthening of the magnetic fields near the equatorial plane in run R22. The MRI-induced turbulence plays an important role in magnetic field amplification. In run R22, the initial magnetic fields are consistent with those observed in other runs in terms of relative strength and topology. Additionally, the growth time of MRI is also the same in all the runs. It is reasonable to believe that the MRI-induced turbulence leads to a similar effect of magnetic field amplification in all the runs. In addition to the magnetic field amplification resulting from MRI-induced turbulence, the magnetic field strength in accretion flows can be further enhanced by two primary processes: local contraction and field advection by converging flows. For run R22, on the one hand, when radiative cooling weakens the gas pressure support in the accretion flows, the accretion flows undergo significant vertical contraction. On the other hand, magnetic field lines can be advected inward along with the accreting material, transporting magnetic flux from the outer regions to the inner regions and thereby strengthening the field. To quantitatively distinguish between these two mechanisms, we refered to Tchekhovskoy et al. (2011) and Narayan et al. (2012) to introduce the normalized magnetic flux:

$$\begin{array}{c}\hfill \dot{\mathrm{\Phi }}=\frac{\mathrm{\Phi }}{\sqrt{{\dot{M}}_{\mathrm{acc}}}},\end{array}$$(12)

where $\mathrm{\Phi }=\sqrt{4\pi }{⟨{\displaystyle {\int }_{86°}^{94°}}|B_{\varphi }|rdrd\theta ⟩}_t$ is the time-averaged toroidal magnetic flux around the equatorial plane and ṁ_acc is the net accretion rate (see Table 1). As previously mentioned, the toroidal magnetic field is the primary component near the equatorial plane; therefore, we only consider the magnetic flux of the toroidal magnetic field here. In Equation (12), the denominator can reflect the contribution of converging flows to transporting magnetic flux. If only the magnetic flux transport by converging flows is taken into account, both the numerator and denominator in Equation (12) increase simultaneously. A comparatively small value of $\dot{\mathrm{\Phi }}$ indicates that converging flows play a significant role in the transport of magnetic flux. In contrast, a comparatively large value of $\dot{\mathrm{\Phi }}$ suggests that the impact of converging flows on magnetic flux transport is relatively minor. Figure 7 displays the radial profiles of $\dot{\mathrm{\Phi }}$ for all runs. For runs N24, R24 and R23, the radial profiles of $\dot{\mathrm{\Phi }}$ are similar. This suggests that the magnetic flux transport by converging flows plays a similar role in strengthening magnetic fields for the three runs. In the region of r = ∼0.2 − 0.3r_B, the value of $\dot{\mathrm{\Phi }}$ in run R22 is significantly greater than that observed in the other three runs. Comparatively, in run R22, the impact of converging flows on magnetic flux transport is weaker than in the other three runs in the region of r = ∼0.2 − 0.3r_B. This region also corresponds to the area with the lowest disk scale height, as shown in Figure 3, and to the region of high density depicted in Figure 4. Therefore, for run R22, it is reasonably inferred that the increase in magnetic field strength is predominantly attributed to the contraction of accretion flows induced by radiative cooling, rather than the magnetic flux transport caused by converging flows.

Fig. 7. Radial profiles of normalized magnetic flux around equatorial plane. The black, blue, red, and green lines correspond to runs R22, R23, R24, and N24, respectively.

3.2. Heating and cooling

Figure 8 shows the radial profiles of the time- and angle- averaged temperature in the vicinity of the equatorial plane. As shown in Figure 8, the gas temperature in runs R23, R24 and N24 exhibits similar profiles. The gas temperature in run R23 is slightly lower than that in runs R24 and N24 in the region of r >  0.15r_B. For run R22, due to strong radiative cooling, the gas temperature is significantly lower than in other runs in the region. The distribution of gas temperature is attributed to cooling and heating. To analyze the energetics of accretion flows, we follow Yuan & Narayan (2014) to define the advection factor (f) of accretion flows as follows:

$$\begin{array}{c}\hfill f=\frac{Q_{\mathrm{adv}}}{Q^{+}}=1-\frac{Q^{-}}{Q^{+}}·\end{array}$$(13)

Fig. 8. Radial profiles of the time-averaged and angle-averaged temperature near the equatorial plane. The angle averaging is implemented over the angle range of 84° to 96°. The black, blue, green, and red lines correspond to runs R22, R23, R24, and N24, respectively.

Here, Q⁺ is the heating rate, Q⁻ is the cooling rate given by Equation (5), and Q_adv is the advection rate of the gas internal energy.

In magnetic fluids, heating primarily results from the dissipation of magnetic energy. The magnetic energy can be converted into thermal energy through reconnection in current sheets. In non-conservative numerical schemes, an explicit artificial resistivity term is often incorporated into numerical codes to implement the dissipation of magnetic energy and to effectively capture heating at current sheets (e.g. Stone & Pringle 2001). In conservative numerical schemes, the energy within computational cells can be exchanged between kinetic, magnetic, and thermal forms. In the code of conservative numerical schemes, the inclusion of an artificial resistivity term is unnecessary for capturing heating at current sheets. Numerical simulations in the shearing box have demonstrated that the heating rate equals the work done by the Maxwell and Reynolds stresses at the box boundaries (Hawley et al. 1995; Sano et al. 2004; Gardiner & Stone 2005; Jiang et al. 2013). Jiang et al. (2014) employed this result to investigate the heating rate in super-Eddington accretion flows. Due to the conservative nature of the numerical schemes employed in this study, here we followed Jiang et al. (2014) to calculate the heating rate. The local heating rate (Q⁺) is approximately given by

$$\begin{array}{c}\hfill Q^{+}=1.5{\mathrm{\Omega }}_{\mathrm{k}}W,\end{array}$$(14)

where Ω_k is the Keplerian angular velocity, and W is the Rϕ component of the combined Reynolds and Maxwell stresses in cylindrical coordinates. In spherical coordinates, the stresses are written as

$$\begin{array}{c}\hfill W=-(B_{r}sin\theta +B_{\theta }cos\theta )B_{\varphi }+\rho (v_{r}sin\theta +v_{\theta }cos\theta )\delta v_{\varphi },\end{array}$$(15)

where $\delta v_{\varphi }=v_{\varphi }-\overline{v_{\varphi }}$ is the difference between v_ϕ and the time-averaged v_ϕ.

The advection factor (f) defined in Equation (13) is often used to measure the relative importance of advection (Yuan & Narayan 2014). Figure 9 shows the radial profiles of the time-averaged and angle-averaged f near the equatorial plane. In runs R23 and R24, f ≈ 1, which implies that the advection cooling plays a significant role in the cooling process. For run R22, within 0.1r_B, f ≈ 1 and then the advection cooling is a dominant cooling term. Outside 0.1r_B, f is less than 1 and even drops to negative values at some radii. This implies that radiative cooling plays a significant role in the energetics outside 0.1r_B for run R22, and at some radii the radiative cooling rate exceeds the heating rate.

Fig. 9. Radial profiles of the time-averaged and angle-averaged advection factor f near the equatorial plane. The angle averaging is implemented over the angle range of 84° to 96°. The black, blue, green, and red lines correspond to runs R22, R23, R24, and N24, respectively.

For the accretion flow near a BH, when Q_adv ≈ Q₊ ≫ Q₋ (i.e. f ≈ 1), the accretion flow is radiatively inefficient and then becomes advection-dominated (Narayan & Yi 1994, 1995). As the accretion rate increases, f < 1 and then the accretion flow becomes a luminous hot accretion flow (Yuan 2001). Bu & Gan (2018) simulated the accretion flows with different accretion rates, and found that when radiative cooling is strong, the cooling rate near the equatorial plane exceeds the heating rate at some radii.

3.3. Convective stability

For a non-rotating flow, an inward increase of entropy can be used as a criterion for convective instability. For a rotating flow, we can employ the Høiland criterion, which is given by (Begelman & Meier 1982; Yuan et al. 2012; Bu et al. 2016b)

$$\begin{array}{c}\hfill (\mathrm{\nabla }s·\mathit{dr})(\mathit{g}·\mathit{dr})-\frac{2\gamma v_{\varphi }}{R^2}[\mathrm{\nabla }(v_{\varphi }R)·\mathit{dr}]dR<0,\end{array}$$(16)

where R = rsinθ is the cylindrical radius, $\mathit{dr}=dr\widehat{r}+rd\theta \widehat{\theta }$ is the displacement vector, s = ln(p_g)−γln(ρ), and $\mathit{g}=-\mathrm{\nabla }\mathrm{\Psi }+(v_{\varphi }^2/R)\widehat{R}$ is the effective gravity. When the accretion flow satisfies Equation (16), it becomes convectively unstable. Figure 10 shows the result obtained from Equation (16) based on simulation data at 5t_B. In Figure 10, the red solid lines mark the density scale height of accretion flows and also indicate the disk surface of accretion flows, as shown in Figure 2. The black regions are convectively unstable, while the white regions are convectively stable. We calculated the proportion of the area occupied by the stable regions between the two red lines for each snapshot, and subsequently averaged those proportions during the time interval of 4 to 5 t_B. Table 1 gave the averaged proportion. For runs N24 and R22–R24, the stable regions account ∼55−62%. This implies that the accretion flows around the sub-Bondi radius are marginally stable in convectively stability. For the accretion flows around the BH, Yuan et al. (2012) employed the Høiland criterion to find that they are convectively stable. In comparison to the accretion flows around the BH, the accretion flows around the sub-Bondi radius promote the occurrence of the convective instability.

Fig. 10. Convective stability analysis. The result is obtained according to Equation (16) based on simulation data at 4.5tB. The red solid lines mark the density scale height of accretion flows, as shown in Figure 2. Panel (A): Runs R22 (left) and R23 (right). Panel (B): Runs R24 (left) and N24 (right). The black regions are unstable.

3.4. Outflows

Following Stone et al. (1999), we define the mass inflow rate (ṁ_in) and outflow rate (ṁ_out) as:

$$\begin{array}{cc}& {\dot{M}}_{\mathrm{in}}(r)=-2\pi r^2{\displaystyle {\int }_0^{\pi }}\rho \min (v_r,0)sin\theta \mathrm{d}\theta ,\hfill \end{array}$$(17)

$$\begin{array}{cc}& {\dot{M}}_{\mathrm{out}}(r)=2\pi r^2{\displaystyle {\int }_0^{\pi }}\rho \max (v_r,0)sin\theta \mathrm{d}\theta .\hfill \end{array}$$(18)

Then, the net mass accretion rate (ṁ_net) is given by

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{net}}(r)={\dot{M}}_{\mathrm{in}}(r)-{\dot{M}}_{\mathrm{out}}(r).\end{array}$$(19)

We first calculate the mass inflow and outflow rates and the net accretion rate from each snapshot and then time-average them.

Figure 11 shows the radial profiles of the mass inflow and outflow rates and the net accretion rate for runs R22–R24 and N24. We can see that the net accretion rate remains approximately constant within ∼0.3 r_B, implying that all the runs achieve a quasi-steady state within this region. As also shown in Figure 11, the mass inflow and outflow rates decrease inward, similar to the behavior observed in the simulations of accretion flows around the BH (Stone et al. 1999; Yuan et al. 2012). A power-law function of radius can be employed to approximately describe the radial profiles of ṁ_in and ṁ_out. We employ the two functions ${\dot{M}}_{\mathrm{in}}(r)\propto r^{s_{\mathrm{in}}}$ and ${\dot{M}}_{\mathrm{out}}(r)\propto r^{s_{\mathrm{out}}}$ to fit the radial profiles of the inflow and outflow rate, respectively, within the radial range of 0.03−0.3 r_B. The values obtained from the fitting are given in Table 1. According to Table 1, for runs R23, R24, and N24, the value of s_in ranges from ∼0.76 to 1.29, while s_out ranges from ∼1.34 to 1.68. For run R22, the values of s_in and s_out are relatively large, being 1.64 and 1.94, respectively. For hot accretion flows near the BH (i.e. at small radii), Yuan et al. (2012) give s_in ∼ 0.4 − 0.75. In comparison to the hot accretion flows near the BH, for the accretion flows around the sub-Bondi radius, the mass inflow rate decreases inward more rapidly.

Fig. 11. Radial profiles of the time-averaged (from 4tB to 5tB) mass inflow rates (solid lines), outflow rates (dashed lines) and the net accretion rates(dotted lines). The black, blue, red, and green lines correspond to runs R22, R23, R24, and N24, respectively.

In this paper, the matter moving outward is defined as outflows, while the outward-moving matter with the Bernoulli parameter (Be) greater then zero is defined as winds. In magnetic flows, the Bernoulli parameter is given by Zhu & Stone (2018):

$$\begin{array}{c}\hfill \mathrm{Be}=\frac{{\mathit{v}}^2}{2}+\mathrm{\Psi }+h+\frac{B_{\varphi }{B}_{\varphi }}{\rho }-\frac{B_{\varphi }{v}_{\varphi }}{k}·\end{array}$$(20)

Here, $h=\frac{\gamma p}{(\gamma -1)\rho }$, and $k=\rho \frac{\sqrt{v_r^2+v_{\theta }^2}}{\sqrt{B_r^2+B_{\theta }^2}}$. For winds, Be >  0 implies that they have enough energy to overcome the gravitational potential well and escape to infinity. The inward decrease of the mass inflow rate was attributed to turbulent motion or mass loss caused by winds (Begelman 2012; Narayan et al. 2000). In the accretion flows with magnetic fields, turbulent motion may be induced by MRI or/and convective instability (Stone et al. 1999; Stone & Pringle 2001). Based on numerical simulations, Narayan et al. (2012) and Yuan et al. (2012) found that the hot accretion flows around the BH are convectively stable in the presence of magnetic fields. Yuan et al. (2015) further found that winds are present in the hot accretion flows and play a key role in the inward decrease of the mass inflow rate. Bu & Gan (2018) found that radiative cooling significantly weakens the winds in the hot accretion flows.

For the hot accretion flows around the sub-Bondi radius (i.e. all the runs in this work), we plot Figure 12 to compare the radial profile of the mass flux of both outflows and winds. Here, the mass flux of winds was calculated by only considering the contributions from the outflows with Be >  0. The total mass outflow rate encompasses the contributions from all outflows, including the outflows with Be >  0 as well as those with Be ≤ 0. In Figure 12, it is evident that the mass flux of winds decreases as the flows move inward. This trend is consistent with the inward decrease in the mass flux of outflows. We also plot Figure 13 to indicate the ratio of the mass flux of winds to the total mass outflow rate. As shown in Figure 13, for run N24, the winds account for 50%−70% of the total mass outflow rate within the quasi-steady state region. For run R22 with strong radiative cooling, the contribution of winds to the mass outflow rate is comparatively lower than in other runs. Specifically, in the region of r >  0.07r_B, the contribution of winds falls below 50%. This implies that the mass outflow rate of winds is significantly reduced due to strong radiative cooling.

Fig. 12. Radial profile of the mass flux of both outflows and winds. Solid lines represent the mass flux of outflows, while dashed lines denote the mass flux of winds. The black, blue, red, and green lines correspond to runs R22, R23, R24, and N24, respectively.

Fig. 13. Ratio of the mass flux of the winds (defined as outflows with Be >  0) to the total mass outflow rate. The black, blue, red, and green lines correspond to runs R22, R23, R24, and N24, respectively.

In Figure 14, it was observed that for the region of r < ∼0.08r_B, the mass flux of winds within the disk is larger than that of winds outside the disk. Conversely, in the region of r >  ∼0.08r_B, the mass flux of winds within the disk is lower than (or comparable to) that of winds outside the disk. For run N24 without radiative cooling, the mass flux of winds within and outside the disk is approximately comparable. As radiative cooling becomes stronger, the winds within the disk is suppressed. In run R22, for the region of r >  ∼0.1r_B, the mass flux of winds within the disk is significantly lower than the winds outside the disk.

Fig. 14. Ratio of the mass flux of winds within disk to the mass flux of winds outside disk. The black, blue, red, and green lines correspond to runs R22, R23, R24, and N24, respectively.

To determine the contribution of winds within and outside the disk to the mass outflow rate, we calculated the two ratios. These ratios were presented in Figure 15. In Figure 15, the upper panel shows the ratio of the mass flux of winds within the disk to the total mass outflow rate within the disk. The lower panel shows the ratio of the mass flux of winds outside the disk to the total mass outflow rate outside the disk. As shown in the upper panel, for run N24 without radiative cooling, the contribution of the winds within disk account for more than 50 percent of the mass outflow rate within the disk at most radii. As radiative cooling gradually becomes stronger, the contribution of the winds within the disk to the mass outflow rate gradually diminishes. For run R22 with strong radiative cooling, the contribution of the winds within the disk to the mass outflow rate is lower than 40 percent in the region where r >  0.06r_B, and drops below 15 percent for r >  0.1r_B. The lower panel also shows that, for run N24, the contribution of the winds outside the disk account for more than 60 percent of the mass outflow rate outside the disk at most radii. For run R22, outside the disk, the contribution of winds to the mass outflow rate is significantly lower than that in other runs. In run R22, for the region where r >  ∼0.1r_B, the contribution of winds to the mass outflow rate outside the disk decreases to less than 50%. Therefore, for runs with weak radiative cooling, the mass outflow rate outside the disk is dominated by winds. Additionally, radiative cooling also reduces the mass outflow rate of winds outside the disk, particularly at large radii.

Fig. 15. Ratio of winds within and outside the disk to the mass outflow rate. The black, blue, red, and green lines correspond to runs R22, R23, R24, and N24, respectively. Panel (A): the ratio of the mass flux of winds (defined as outflows with Be >  0) within the disk to the total mass outflow rate within the disk. Panel (B): the ratio of the mass flux of winds outside the disk to the total mass outflow rate outside the disk.

In addition to the winds that escape to infinity, the outflows that do not reach infinity should be attributed to the outward-moving matter in turbulence. For run N24 without radiative cooling and run R24 with weak radiative cooling, the mass flux of winds exceeds that of the outward-moving matter within turbulent eddies. Therefore, in the two runs, the inward decrease of mass inflow rate is mainly attributed to winds. As radiative cooling gradually becomes stronger, the winds weaken, leading to a reduced influence on the inward decrease of the mass inflow rate, while the role of turbulence becomes increasingly important. For run R22 with strong radiative cooling, in the region where r < 0.05r_B, the contributions of winds and turbulent motions to the inward decrease of mass inflow rate are roughly equal. However, in the region where r >  0.05r_B, the inward decrease of mass inflow rate is mainly attributed to turbulent motions. Panel A in Figure 2 shows a snapshot for run R22. As shown in this snapshot, turbulence is significant.

In hot accretion flows, both gas pressure and magnetic forces are commonly considered to play significant roles in the mechanisms that drive winds. Blandford & Payne (1982) suggested the magneto-centrifugally driven mechanism to drive the thin disk winds. Their model necessitates the presence of a large-scale poloidal magnetic field. However, in hot accretion flows, the magnetic field is relatively disordered. Yuan et al. (2012) also proposed an alternative mechanism for wind production, known as the ‘micro’ Blandford–Payne mechanism. This mechanism indicates that in small-scale turbulent magnetic fields, the combined effects of centrifugal force, magnetic pressure gradient, and gas pressure gradient can drive winds. In the ‘micro’ Blandford-Payne mechanism, the centrifugal force is primarily responsible for accelerating the winds, while the gas pressure gradient force and the magnetic pressure gradient force assist in altering the direction of the outflows or play a comparable role compared with the centrifugal force. To investigate the mechanisms driving winds, we examined the forces within accretion flows in two runs: run N24 without radiative cooling and run R22 with strong radiative cooling. In Figure 16, we presented the time-averaged forces at r = 0.15r_B on the disk surface (as indicated by the black line in Figure 2) for runs R22 and N24. In the two runs, the centrifugal force plays an important role in driving the winds. In run R22, the net force is directed outward and is weaker than that in run N24. In run N24, the net force is also directed outward and slightly upward. Consequently, the winds in run R22 are significantly weaker than those in the other runs. Therefore, our analysis indicates that the mechanism responsible for driving winds within the accretion flow around the sub-Bondi radius is the ‘micro’ Blandford–Payne mechanism.

Fig. 16. Force analysis at r = 0.15rB on the disk surface (as indicated by the black line in Figure 2). The black, blue, red, green, and yellow colors correspond to the net force, gravitational force, centrifugal force, pressure gradient force, and magnetic pressure gradient force, respectively. The arrows point in the direction of each force, and their lengths represent the magnitude of the force. All arrows are scaled identically, with the unit being the magnitude of the gravitational force. Panel (A): Run R22. Panel (B): Run N24.

3.5. Inflows

Based on HD simulations, Li et al. (2013) investigated the effects of radiative cooling on hot accretion flows around the Bondi radius. Their results indicated that as the accretion rate increases, strong radiative cooling promotes accretion. Here, we used MHD simulations to further investigate the effects of radiative cooling on inflows. If the radiative cooling term is not taken into account, the dynamical equations can become density free. Consequently, for the cases without radiative cooling, the results of numerical simulations can be applied to any units of density, and then the accretion rate is directly proportional to the adopted density in simulations. In our simulations, the density of accretion flows is scaled by the density at the center of the initial torus, i.e. ρ₀. In Figure 17, we plot the dependence of the net accretion rate on ρ₀. The red points correspond to the cases without radiative cooling, whereas the black points correspond to the cases with radiative cooling. For run R24, radiative cooling slightly reduces the net accretion rate due to weak radiative cooling. For run R22 with strong radiative cooling, the net accretion rate is reduced by nearly an order of magnitude compared to the runs without radiative cooling. This suggests that strong radiative cooling suppresses accretion. This result is not consistent with that given by Li et al. (2013). In HD simulations, an α-viscosity model is adopted to drive angular momentum transport and accretion. In the α-viscosity prescription, the effectiveness of angular momentum transport depends on the specific angular momentum flux driven by the anomalous viscosity shear stress. This shear stress is directly proportional to the kinematic viscosity (ν). In the study by Li et al. (2013), they assumed a constant, non-zero value for kinematic viscosity (ν). This assumption led to uniform effectiveness in angular momentum transport in all their simulations. In cases of significant radiative cooling, the net accretion rate increases as a result of a reduction in the mass outflow rate caused by the strong cooling effect (refer to Figure 6 in Li et al. 2013). However, in MHD simulations, angular momentum transport is driven by Maxwell and Reynolds stresses. To understand the reduction in the mass accretion rate observed in run R22, we examine the effective viscosity parameter (α) and the specific angular momentum flux in MHD simulations.

Fig. 17. Net accretion rates corresponding to different unit densities, with black and red points representing the cases with radiative cooling switched on and off, respectively.

In accretion flows, Maxwell stress (S_m) and Reynolds stress (S_r) are responsible for the angular momentum transfer. We then followed Jiang et al. (2019) to calculate the radial profiles of the r − ϕ stress as S_m = ⟨⟨−B_rB_ϕ⟩⟩ and S_r = ⟨⟨ρv_rv_ϕ⟩⟩ − ⟨⟨ρv_r⟩⟩⟨⟨v_ϕ⟩⟩. The effective α can be defined as α_m ≡ S_m/⟨⟨P_tot⟩⟩, α_r ≡ S_r/⟨⟨P_tot⟩⟩, and then α ≡ α_m + α_r. Here, P_tot is the total pressure, which encompasses both the gas pressure and the magnetic pressure. The specific angular momentum flux caused by Maxwell and Reynolds stresses is defined as

$$\begin{array}{c}\hfill t_{\mathrm{r}\varphi }=\frac{\alpha r⟨⟨P_{\mathrm{tot}}⟩⟩}{⟨⟨\rho ⟩⟩}·\end{array}$$(21)

In Figures 18 and 21, we plotted the radial profiles of the effective α and the specific angular momentum flux, respectively.

Fig. 18. Radial profiles of the effective α. The black, blue, red, and green lines correspond to runs R22, R23, R24 and N24, respectively.

According to Figure 18, in the runs (R24 and R23) with weak radiative cooling, the effective α varies from 0.16 to 0.42, with the average value of 0.26. In run R22 with strong radiative cooling, the effective α is higher compared to other runs within the region where r ≲ 0.06r_B. Conversely, for regions where r >  0.08r_B, the effective α is lower than those observed in the other runs. In Figure 19, the effective α_m is is significantly greater than the effective α_r. As a result, the effective α is primarily determined by the effective α_m. This suggests that Maxwell stresses contribute significantly to the effective α. The r − ϕ component of Maxwell stresses depends on the radial and toroidal components of magnetic fields. As shown in Figure 20, the magnetic fields are dominated by the toroidal component. In the region where r >  ∼0.1r_B, the radial magnetic field is relatively weak in run R22, leading to a reduction in both the effective α_m and α values.

Fig. 19. Radial profiles of the effective αm (solid lines) and αr (dashed lines). The black, blue, red, and green lines correspond to runs R22, R23, R24 and N24, respectively.

Fig. 20. Radial profiles of the percentage of the total magnetic field energy accounted for by the radial (solid lines) and toroidal (dashed lines) components of the magnetic field, respectively. The black, blue, red, and green lines correspond to runs R22, R23, R24 and N24, respectively.

Figure 21 shows that for run R22 with strong radiative cooling, the specific angular momentum flux driven by Maxwell and Reynolds stresses is significantly reduced in the region of r >  ∼0.08r_B. This is because the specific angular momentum flux is influenced not only by the effective α parameter but also by the temperature of accretion flows. In run R22, radiative cooling reduces the temperature of accretion flows (see Figure 8). A significant reduction in the specific angular momentum flux at r = ∼0.4r_B results in the significant decrease in the mass inflow rate (refer to Figure 11). Consequently, as shown in Figure 4, the radial density profile at the equator in run R22 differs significantly from those observed in other runs with weak radiative cooling. Due to the significant decrease in the mass inflow rate, the net accretion rate in run R22 is significantly reduced compared to the cases with weak radiative cooling (refer to Figure 17). Therefore, accretion in run R22 is inhibited due to significant radiative cooling.

Fig. 21. Radial profile of the stress tensor exerted on per unit mass for different runs. The black, blue, red, and green lines correspond to runs R22, R23, R24 and N24, respectively. rB and vk, B denote the Bondi radius and the Keplerian velocity at the Bondi radius, respectively.

A low-temperature, high-density accretion disk may become gravitationally unstable, which may result in fragmentation and star formation. Gravitational instability can be assessed using Toomre’s criterion (Toomre 1964):

$$\begin{array}{c}\hfill Q=\frac{C_s\mathrm{\Omega }}{\pi G\mathrm{\Sigma }},\end{array}$$(22)

where Σ is the surface density of the disk. When Q ≥ 1, the disk is gravitationally stable. When Q < 1, the disk is gravitationally unstable. Gravitational stress and Reynolds stress, associated with the gravitational instability-driven turbulence, can provide angular momentum transport in accretion flows. When Q ≪ 1, the disk is expected to be fragment. We computed the Q values and found Q ≫ 1 for all the runs, which indicates that the disks are gravitationally stable and angular momentum transport is still attributed to Maxwell stresses associated with MHD turbulence driven by the MRI.

4. Discussions

For the hot accretion flows near a BH, synchrotron radiation, Compton scattering, and two-temperature effects play significant roles. Nevertheless, these effects are not incorporated into our simulations. The reasons are as follows.

(1) When the electron temperature is higher than ∼10¹⁰ K, synchrotron radiation becomes important (refer to Figure 1 in Fragile & Meier 2009). In our simulations, the gas temperature is below 10⁸ K, which is insufficient to produce a significant population of relativistic electrons. Therefore, we did not take into account synchrotron radiation.

(2) For the low-luminosity AGNs, Compton scattering plays an important role. The hot accretion flows near a BH can produce a great number of high-energy photons, which can heat the gas at parsec and sub-parsec scales. The significance of this effect depends on the luminosity (L) of the hot accretion flows. Yuan & Li (2011) pointed out that when the luminosity exceeds 2% of Eddington luminosity, the gas at parsec and sub-parsec scales is heated to sufficiently high temperatures. As a result, the gas becomes unbound by the gravity of the central BH, thereby halting the accretion process. After the accretion process ceases, there is a significant decrease in luminosity, rendering Compton heating negligible. Subsequently, the accretion process resumes, leading to intermittent activity. In our simulations, the accretion rate remained consistently low, with a maximum value of 2.5 × 10⁻⁶ Ṁ_Edd. On the other hand, our computational inner boundary was positioned at a considerable distance from the BH. Subsequently, it was unknown how much matter falls into a BH, resulting in difficulty in accurately calculating the luminosity of the hot accretion flows. According to Wang et al. (2013), in the matter initially captured by the central BH in Sgr A^★, less than 1% reaches the innermost region around Sgr A^★. Therefore, we believed that for our runs, their luminosity is very low and thus Compton scattering effects are negligible.

(3) It is widely accepted that hot accretion flows surrounding a BH exhibit a two-temperature structure, attributed to the weak Coulomb coupling between ions and electrons in the very hot plasma. However, according to Narayan & Yi (1995; refer to Figure 4 in their paper), the temperature gap between ions and electrons diminishes with increasing radial distance and becomes negligible beyond 10³ r_s. Therefore, in our simulations, two-temperature effects are negligible.

In our simulations, we employed a torus with high angular momentum as the initial condition. Although the specific angular momentum of gases near the Bondi radius (∼11 parsec) is observationally unknown, it is reasonable to infer that the angular momentum of gases in many galaxies is likely not small. Recent kinematic surveys have indicated that approximately 80% of early-type galaxies (ETGs) are classified as regular rotators (Emsellem et al. 2011). These galaxies exhibit oblate axisymmetric shapes, which suggest the presence of underlying disk-like components. Yoon et al. (2018) implemented simulations of elliptical galaxies. In their simulations, the radial range extends from 2.5 parsec to 250 kiloparsec. In simulations of galaxy evolution, gases are introduced through various mechanisms, including AGN mass feedback, stellar mass loss, and supernova explosions. Their results indicate that the gases within a radius of r < 100 pc have high angular momentum, close to the Keplerian specific angular momentum (refer to Figure 8 in Yoon et al. 2018). Gan et al. (2019) further took into account the infall of the circumgalactic medium, and obtained similar results (refer to Figure 6 in Gan et al. 2019).

For all the runs, we considered only the gravitational potential of a non-spinning BH. Since our computational domain is located far from the BH, the effects of BH spin become negligible. However, near the BH, BH spin plays an important role and significantly affect the production and properties of winds. Previous studies have discussed the influence of BH spin on winds in accretion flows (Curd & Narayan 2023; Aktar et al. 2024). These winds, launched from the inner region, may propagate outward to the Bondi radius, where they may interact with the interstellar medium on larger scales. Nevertheless, this effect has not yet been investigated in the current study, but it remains an interesting topic for future investigation.

5. Summary

The accretion processes surrounding a SMBH in AGNs are thought to span a region extending from the vicinity of the black hole to the Bondi radius. It is difficult to implement numerical simulations on such a large region. Most previous numerical simulations have primarily concentrated on the region close to the BH. These studies have indicated that strong winds can be produced within hot accretion flows near the BH and that radiative cooling plays an important role in modulating the strength of these winds. In this paper, we performed 2D MHD simulations to investigate the impact of radiative cooling on hot accretion flows around the sub-Bondi radius. We exclusively considered radiative cooling through bremsstrahlung radiation. As the mass accretion rate increases, radiative cooling gradually becomes stronger, leading to the contraction of accretion flows and a reduction in the thickness of the accretion disk. In the case of strong radiative cooling, the significant contraction of accretion flows leads to a substantial increase in magnetic pressure within the accretion disk. The radial profiles of the mass inflow and outflow rates can be described by ${\dot{M}}_{\mathrm{in}}(r)\propto r^{s_{\mathrm{in}}}$ and ${\dot{M}}_{\mathrm{out}}(r)\propto r^{s_{\mathrm{out}}}$, respectively. For the runs with weak cooling, the exponent s_in is approximately within the range of 0.7 to 1.3, and the exponent s_out is approximately within the range of 1.2 to 1.7. For the run with strong radiative cooling, the inflow and outflow profiles are steeper, with s_in ≈ 1.6 and s_out ≈ 1.9. Similarly, the radial profile of density near the equatorial plane can be expressed as ρ(r)∝r^−p. For the runs with weak radiative cooling, the exponent p is approximately within the range of 0.2 to 0.5, while for the runs with strong radiative cooling, the density profile is flatter, with p ≈ 0.05.

We employed the Høiland criterion to analyze the convective stability of the accretion disk. We find that within the accretion disk (defined as the accretion flows within the density scale height), the region of convection stability accounts for ∼55−62%. This suggests that the hot accretion flows around the sub-Bondi radius are marginally convectively stable.

We defined the outflows with Be >  0 as winds that escape to infinity. In the runs with weak radiative cooling, the winds contribute to ∼50−70% of the mass outflow rate, and then they play a significant role in the inward decrease of the mass inflow rate. In the runs with strong radiative cooling, the mass outflow rate of winds is significantly reduced, particularly at large radii, and then the inward decrease of the mass inflow rate is primarily attributed to turbulence driven by MRI and convection. Radiative cooling suppresses the production of winds.

Furthermore, in the region of r ≳ 0.08r_B, the effective viscosity parameter (α) in the runs with strong radiative cooling is significantly lower than in those with weak radiative cooling. As a result, accretion processes are suppressed under conditions of strong radiative cooling. This result differs from the results obtained from HD simulations.

Acknowledgments

We thank the anonymous referee for the valuable suggestions, which were very helpful in improving the manuscript. This work is supported by the Natural Science Foundation of China (grant 12347101). D. Bu is supported by the Natural Science Foundation of China (grant 12192220, 12192223).

References

Appendix A: Assessing the outer boundary sensitivity and the impact of simulation dimensionality

We implemented run N24b to assess the sensitivity to the outer radius of the computational region, and run N24c to evaluate the impact of simulation dimensionality. In the following analysis, we perform time-averaging over the time interval of 2−4 t_B. During this time interval, we gather 200 data files with equal time intervals. Run N24b is an alternative version of run N24a. Compared to run N24a, run N24b has a larger outer boundary. The outer boundary of run N24b is set at 10r_B. In run N24b, the grid resolution is consistent with that of Run N24a within r < 5r_B. In the radial range extending from r = 5r_B to 10r_B, we employ 68 uniformly distributed zones. Run N24c is a 3D version of run N24a. In the r and θ directions, runs N24c and N24a have the same number of zones. In the ϕ direction, we employ 64 uniformly-spaced zones for run N24c. Compared to run N24a, runs N24b and N24c exhibit qualitatively consistent properties, as evidenced by the data presented in Table A.1. Figure A.1 shows the radial profiles of the net accretion rate, inflow rate, and outflow rate for comparing runs N24a (black lines) and N24b (blue lines). Figure A.2 also presents the radial density profiles at the equatorial plane. Figures A.1 and A.2 demonstrate that the radial profiles of these quantities are comparable in runs N24a, N24b and N24c. This suggests that our simulation results are not sensitive to the position of the outer boundary, provided the outer boundary is situated at a substantial distance from the initial torus, and also indicates that the axisymmetric 2D simulations are sufficient to capture the key features of the accretion flow dynamics in our study.

Fig. A.1. Radial profiles of the time-averaged mass inflow rates (solid lines), outflow rates(dashed lines) and net accretion rates(dotted lines). The black, blue and red lines correspond to runs N24a, N24b and N24c, respectively.

Fig. A.2. Radial profiles of the time-averaged and angle-averaged density near the equatorial plane. The angle averaging is implemented over the angle range of 84° to 96°. The black, blue and red lines correspond to runs N24a, N24b and N24c, respectively.

All Tables

All Figures

	Fig. 1. Distribution of the initial magnetic field (lines) and plasma β0 (colors). Solid and dashed lines represent clockwise and counterclockwise magnetic field lines, respectively.
In the text

Fig. 2. Distribution of density and velocity field. In each panel, the left side presents a snapshot at 4.0tB, while the right side displays the time-averaged distribution. Colors represent the density, measured in units of the density at the center of the initial torus. The streamlines reflect the velocity field on the poloidal plane. The black solid line marks the scale height of accretion flows. The scale height is defined as the height at which the density drops by a factor of e relative to the density at the equatorial plane. Panel (A): Run R22. Panel (B): Run R23. Panel (B): Run R24. Panel (D): Run N24.

In the text

	Fig. 3. Radial dependence of the ratio (H/r) of density scale height to radius. The black, blue, red, and green lines correspond to runs R22, R23, R24, and N24, respectively.
In the text

	Fig. 4. Radial profiles of the time-averaged and angle-averaged density near the equatorial plane. The angle averaging is implemented over the angle range of 84° to 96°. The black, blue, green, and red lines correspond to runs R22, R23, R24, and N24, respectively.
In the text

Fig. 5. Distribution of the plasma β and poloidal magnetic fields. Colors represent the ratio (β) of the time-averaged gas pressure to the magnetic pressure of the time-averaged magnetic field. The solid and dashed lines reflect the magnetic fields on the poloidal plane. The purple line indicates the interface at which β = 1. Panel (A): Runs R22 (left) and R23 (right). Panel (B): Runs R24 (left) and N24 (right).

In the text

	Fig. 6. Distribution of the strength ratio between the toroidal and poloidal components of the mean magnetic field. The purple line indicates the interface where the toroidal magnetic field strength equals the poloidal magnetic field strength. Panel (A): Runs R22 (left) and R23 (right). Panel (B): Runs R24 (left) and N24 (right).
In the text

	Fig. 7. Radial profiles of normalized magnetic flux around equatorial plane. The black, blue, red, and green lines correspond to runs R22, R23, R24, and N24, respectively.
In the text

	Fig. 8. Radial profiles of the time-averaged and angle-averaged temperature near the equatorial plane. The angle averaging is implemented over the angle range of 84° to 96°. The black, blue, green, and red lines correspond to runs R22, R23, R24, and N24, respectively.
In the text

	Fig. 9. Radial profiles of the time-averaged and angle-averaged advection factor f near the equatorial plane. The angle averaging is implemented over the angle range of 84° to 96°. The black, blue, green, and red lines correspond to runs R22, R23, R24, and N24, respectively.
In the text

	Fig. 10. Convective stability analysis. The result is obtained according to Equation (16) based on simulation data at 4.5tB. The red solid lines mark the density scale height of accretion flows, as shown in Figure 2. Panel (A): Runs R22 (left) and R23 (right). Panel (B): Runs R24 (left) and N24 (right). The black regions are unstable.
In the text

	Fig. 11. Radial profiles of the time-averaged (from 4tB to 5tB) mass inflow rates (solid lines), outflow rates (dashed lines) and the net accretion rates(dotted lines). The black, blue, red, and green lines correspond to runs R22, R23, R24, and N24, respectively.
In the text

	Fig. 12. Radial profile of the mass flux of both outflows and winds. Solid lines represent the mass flux of outflows, while dashed lines denote the mass flux of winds. The black, blue, red, and green lines correspond to runs R22, R23, R24, and N24, respectively.
In the text

	Fig. 13. Ratio of the mass flux of the winds (defined as outflows with Be >  0) to the total mass outflow rate. The black, blue, red, and green lines correspond to runs R22, R23, R24, and N24, respectively.
In the text

	Fig. 14. Ratio of the mass flux of winds within disk to the mass flux of winds outside disk. The black, blue, red, and green lines correspond to runs R22, R23, R24, and N24, respectively.
In the text

Fig. 15. Ratio of winds within and outside the disk to the mass outflow rate. The black, blue, red, and green lines correspond to runs R22, R23, R24, and N24, respectively. Panel (A): the ratio of the mass flux of winds (defined as outflows with Be >  0) within the disk to the total mass outflow rate within the disk. Panel (B): the ratio of the mass flux of winds outside the disk to the total mass outflow rate outside the disk.

In the text

Fig. 16. Force analysis at r = 0.15rB on the disk surface (as indicated by the black line in Figure 2). The black, blue, red, green, and yellow colors correspond to the net force, gravitational force, centrifugal force, pressure gradient force, and magnetic pressure gradient force, respectively. The arrows point in the direction of each force, and their lengths represent the magnitude of the force. All arrows are scaled identically, with the unit being the magnitude of the gravitational force. Panel (A): Run R22. Panel (B): Run N24.

In the text

	Fig. 17. Net accretion rates corresponding to different unit densities, with black and red points representing the cases with radiative cooling switched on and off, respectively.
In the text

	Fig. 18. Radial profiles of the effective α. The black, blue, red, and green lines correspond to runs R22, R23, R24 and N24, respectively.
In the text

	Fig. 19. Radial profiles of the effective αm (solid lines) and αr (dashed lines). The black, blue, red, and green lines correspond to runs R22, R23, R24 and N24, respectively.
In the text

	Fig. 20. Radial profiles of the percentage of the total magnetic field energy accounted for by the radial (solid lines) and toroidal (dashed lines) components of the magnetic field, respectively. The black, blue, red, and green lines correspond to runs R22, R23, R24 and N24, respectively.
In the text

	Fig. 21. Radial profile of the stress tensor exerted on per unit mass for different runs. The black, blue, red, and green lines correspond to runs R22, R23, R24 and N24, respectively. rB and vk, B denote the Bondi radius and the Keplerian velocity at the Bondi radius, respectively.
In the text

	Fig. A.1. Radial profiles of the time-averaged mass inflow rates (solid lines), outflow rates(dashed lines) and net accretion rates(dotted lines). The black, blue and red lines correspond to runs N24a, N24b and N24c, respectively.
In the text

	Fig. A.2. Radial profiles of the time-averaged and angle-averaged density near the equatorial plane. The angle averaging is implemented over the angle range of 84° to 96°. The black, blue and red lines correspond to runs N24a, N24b and N24c, respectively.
In the text