Toward Horizon-scale Accretion Onto Supermassive Black Holes in Elliptical Galaxies

We present high-resolution, three-dimensional hydrodynamic simulations of the fueling of supermassive black holes in elliptical galaxies from a turbulent medium on galactic scales, taking M87* as a typical case. The simulations use a new GPU-accelerated version of the Athena++ AMR code, and span more than 6 orders of magnitude in radius, reaching scales similar to the black hole horizon. The key physical ingredients are radiative cooling and a phenomenological heating model. We find that the accretion flow takes the form of multiphase gas at radii less than about a kpc. The cold gas accretion includes two dynamically distinct stages: the typical disk stage in which the cold gas resides in a rotationally supported disk and relatively rare chaotic stages ($\lesssim 10\%$ of the time) in which the cold gas inflows via chaotic streams. Though cold gas accretion dominates the time-averaged accretion rate at intermediate radii, accretion at the smallest radii is dominated by hot virialized gas at most times. The accretion rate scales with radius as $\dot{M}\propto r^{1/2}$ when hot gas dominates and we obtain $\dot{M}\simeq10^\mathrm{-4}-10^\mathrm{-3}\,M_\odot\,\mathrm{yr^{-1}}$ near the event horizon, similar to what is inferred from EHT observations. The orientation of the cold gas disk can differ significantly on different spatial scales. We propose a subgrid model for accretion in lower-resolution simulations in which the hot gas accretion rate is suppressed relative to the Bondi rate by $\sim (r_\mathrm{g}/r_{\rm Bondi})^{1/2}$. Our results can also provide more realistic initial conditions for simulations of black hole accretion at the event horizon scale.


INTRODUCTION
Supermassive black holes (SMBHs) are ubiquitous in the nuclei of elliptical galaxies (Kormendy & Ho 2013). The fueling of these black holes involves a sizeable dynamic range from megaparsec to milliparsec scales and a variety of physical processes including radiative cooling, heating, turbulence, and magnetic fields.
It remains uncertain how SMBHs in galactic nuclei are fed in detail. For spherically symmetric adiabatic accretion of hot gas (Bondi 1952), the accretion begins from a characteristic (Bondi) radius, (r B = 2GM • /c 2 s,∞ , where G is the gravitational constant, M • is the black hole mass, and c s,∞ is the sound speed of ambient gas), where the gravitational energy due to the SMBH becomes comparable to the thermal energy of the gas. For many low-luminosity active galactic nuclei (AGN), it is believed that the quasi-spherical Bondi flow at large radii transitions at smaller radii to a hot, tenuous radiatively inefficient accretion flow (RIAF), with a radiative cooling time-scale longer than the accretion time-scale (Ichimaru 1977;Narayan & Yi 1995). However, if cooling plays an important role, the cold gas will contribute to the accretion and change the accretion flow substantially by forming a geometrically thin disk or chaotic cold accretion flow (Li & Bryan 2012;Sharma et al. 2012;Gaspari et al. 2013), in which the accretion rate is boosted by several orders of magnitude due to cooling of the hot gas.
Contemporary event horizon-scale modeling of accretion flows is limited in part by a dependence on ad hoc initial conditions.
Most previous general relativistic magnetohydrodynamic simulations (GRMHD) (Narayan et al. 2012;Porth et al. 2019;White et al. 2020) set the initial condition as a torus in hydrodynamic equilibrium with a specified angular momentum profile (Fishbone & Moncrief 1976;Kozlowski et al. 1978). However, the flow structure can strongly depend on the initial conditions, forming structures such as a magnetically arrested disk (MAD) (Narayan et al. 2003). More realistic initial and boundary conditions on galactic scales may help to construct a more consistent model of black hole accretion at smaller radii. Connecting the large and small scales in this way is also critical for developing more physical models of black hole growth and feedback in cosmological simulations that lack the physics or resolution to follow the gas to small radii.
Among massive galaxies, the giant elliptical galaxy M87 is of particular interest. The Event Horizon Telescope (Event Horizon Telescope Collaboration et al. 2019) presented the first image of event horizon scale plasma in the galaxy M87, deriving a BH mass of M • = (6.5 ± 0.7) × 10 9 M (consistent with prior stellar dynamics results; Gebhardt & Thomas 2009;Gebhardt et al. 2011). The corresponding gravitational radius is r g = GM • /c 2 = 3.1 × 10 −4 pc. The mass accretion rate near the event horizon is estimated to be (3 − 20) × 10 −4 M yr −1 (Event Horizon Telescope Collaboration et al. 2021). For M87, the temperature of the ambient hot plasma far from the black hole is ≈ 2 × 10 7 K (Russell et al. 2015), which gives a Bondi radius of r B ≈ 0.12 kpc ∼ 4 × 10 5 r g and a Bondi accretion rate ofṀ ≈ 0.05 M yr −1 , much larger than the accretion rate measured on horizon scales. In principle, cooling of hot gas on kpc scales and larger could generate an even larger inflow rate of ∼ 10 M yr −1 , although observations imply that heating (rather than gas inflow) largely offsets this radiative cooling (Fabian 2012a). In addition to the hot virialized plasma, HST observations show that there is a warm ionized T ∼ 10 4 K gas disk within ∼ 40 pc (e.g., Walsh et al. 2013). Despite these constraints on the gas properties ∼ 100 pc from the SMBH, the realistic structure and properties of the accretion flow in M87 and related systems are still uncertain. To self-consistently link the accretion rate at galactic scales to the event horizon, a dynamic range of ∼ 6 orders of magnitude must be covered.
This work is a first step towards building a multi-scale model of the accretion flow onto SMBHs in elliptical galaxies from galactic down to event horizon scales, involving different physics at different scales. Extending previous works (Li & Bryan 2012;Gaspari et al. 2013), our simulations span approximately six orders of magnitude in radius, reaching all the way down to radii pc. As much as possible, we set our initial conditions at large radii using observations (Gebhardt & Thomas 2009;Urban et al. 2011;Russell et al. 2015). The simulations connect the black hole accretion rate at small radii to the physical conditions in the galaxy at larger radii. More broadly, the problems that we address in this work include: the contribution of cold and hot gas to the accretion rate, the radial dependence of cold gas and hot gas, as well as the morphology, angular momentum, and dynamics of the accretion flow. The flow structure found here at small radii can also be used as more realistic initial conditions for simulations of the black hole accretion disk at event horizon scales, which will help to interpret Event Horizon Telescope observations of M87.
The rest of this article is organized as follows. In Section 2 we describe the physical model and numerical methods we adopt. Section 3 presents the results of our simulations. In Section 4, we discuss the implications of our results. We conclude in Section 5.

MODEL AND NUMERICAL METHODS
We perform hydrodynamic simulations using a performance portable version of the Athena++  code implemented using the Kokkos library (Trott et al. 2021). Athena++ provides a variety of reconstruction methods, Riemann solvers, and integrators for solving the equations of hydrodynamics. The static mesh refinement (SMR) in Athena++ enables us to achieve a high resolution and good performance over an extremely large dynamic range. We perform the simulations on Cartesian grids to avoid the severe time-step restriction and discontinuities at the pole inherent in 3D sphericalpolar coordinates. In addition, there is no a priori symmetry axis in our problem, making Cartesian grids particularly suitable.
In our simulations, we adopt the piecewise linear (PLM) reconstruction method, the Harten-Lax-van Leer-Contact (HLLC) Riemann solver, and the RK2 time integrator to solve the hydrodynamic equations. We also tested the piecewise parabolic (PPM) reconstruction method and the RK3 time integrator, and the results are very similar, although the higher-order methods tend to be less stable, leading to problematic cells with unphysically large temperatures or velocities. Thus our production runs adopt second-order methods. The source terms (gravity, cooling, and heating) are included by operator splitting.
In this study, we solve ∂ρ ∂t + ∇ · (ρv) = 0, ∂ρv ∂t + ∇ · (P I + ρvv) = ρg, where ρ is the gas density, v is the velocity, P is the gas pressure, with mean molecular weight µ = 0.618 by assuming the gas is fully ionized, g is the gravitational acceleration, E is the total energy density, with the gas adiabatic index γ = 5/3, Q − is the volume cooling rate, and Q + is the heating rate. The implementations of initial and boundary conditions, gravitational field, radiative cooling, and heating are presented below. . Free-fall time, gas inflow time (assuming a spherical Bondi accretion flow), and gas cooling time as functions of radius. On kpc scales, the gas cooling time is larger than the free-fall time by several times but smaller than the gas inflow time, enabling thermal instabilities to fully develop. The shaded region marks the Bondi radius rB. The violet lines illustrate the radial domain of the simulations and the arrows indicate that we restart the larger-scale simulation with a smaller inner radius.

Initial and Boundary Conditions
We investigate the accretion physics in an elliptical galaxy embedded in hot plasma. For the initial conditions of gas, we assume a flat cored entropy profile, following Martizzi et al. (2019), where the normalized radius is defined as x = r/r 0 . The gas pressure is where n = ρ/(µm H ) is the number density. Then the density profile is obtained by solving hydrostatic equilibrium, where Φ is the gravitational potential. In this work, we set r 0 = 2 kpc, K 0 = 7 keV cm 2 , and ξ = 1.1 to match the X-ray observations of the Virgo Cluster in the work of Urban et al. (2011). The normalization of the density profile is set as n(x = 1) = 0.1 cm −3 by taking the observed density fitting from the observations of M87 by Russell et al. (2015). By adopting these parameters, the radial profiles of density and temperature are similar to observations of the galaxy M87 as well as fairly general to be applicable to similar systems. Some important time-scales of the initial condition are shown in Figure 1. The free-fall time ranges from roughly 10 Myr at 10 kpc to 1 yr at 0.03 pc (100r g ). On kpc scale, the gas cooling time t cool (∼ 10 Myr) is larger than the free-fall time t ff but by no more than a factor of 10, i.e., t cool /t ff 10. Furthermore, assuming a Bondi accretion rate, the cooling time is smaller than the gas inflow time t inflow . This implies that cooling is dynamically important; indeed, as we shall see, thermal instability develops in this region.
The simulations adopt a cubic box of size 64 kpc in each direction and cover a radial domain from 32 kpc down to 30 mpc or smaller, using static mesh refinement to span 6 orders of magnitude in scale. The center of the potential is placed at the center of the simulation box. The root grid is a cube of 128 3 cells, which is divided into 4 × 4 × 4 mesh blocks, with each mesh block being a cube of 32 3 cells. At each level of mesh refinement, we double the resolution by refining the central 2 × 2 × 2 parent mesh blocks into 4×4×4 child mesh blocks. This enables us to retain a resolution of 1/64 < ∆x/x < 1/32 over the whole simulation domain. In the simulations, we typically employ 10 to 20 levels of static mesh refinement. In one case for the largest zoom-in simulation, we use 25 levels with inner radius r in = 1 mpc(≈ 3r g ) and the finest resolution of ∼ 0.05r g . In one high-resolution simulation for testing convergence, we double the resolution to 1/128 < ∆x/x < 1/64 everywhere.
For the outer boundary conditions, we fix the outer region of the box (r > r out = 32 kpc) to be the initial hydrostatic equilibrium solution. This allows us to enforce spherical symmetry at the outer boundaries. Since the outer radius r out is sufficiently large, the fixed outer boundaries have little influence on the flow at small scales. For the inner boundaries, we apply a vacuum sink region with r < r in by evacuating the gas within the region and resetting a negligible fixed density n sink = 5 × 10 −3 cm −3 , temperature T sink = 2 × 10 5 K, and zero velocity every time step. This sink region avoids artificial overpressure bounces and correctly reproduces the Bondi accretion solution with robust numerical stability.
To deal with the excessively severe time-step restriction at small scales, we first evolve the simulations with a relatively larger inner radius r in,old . Then we restart the simulation with a smaller inner radius r in,new . The new cells in the simulation domain are initialized with the values used for the boundary conditions. The evolution time for each run is 10 4 free-fall time at the inner boundary. In Figure 1, we illustrate the radial domain of the simulations with different r in . When restarting with a smaller inner radius, to avoid the sudden change of inner boundary conditions and the unphysical shocks therefrom, we gradually shrink the radius of the inner boundary by where the transition time-scale, t trans , is a few times longer than the local dynamical time at r in so that the change is nearly adiabatic. We tested the effects of the transition time on the accretion flow. The results show little difference for different t trans as long as t trans is long enough to avoid the strong unphysical shocks. The simulations are initialized with random perturbations to model turbulent motions, which is crucial in seeding thermal instability and shaping the accretion flow. We seed thermal instabilities using two types of perturbations, i.e., density perturbations (associated with lower angular momentum) and velocity perturbations (higher angular momentum). In the first type, we set an isobaric density perturbation, in which the velocity of the perturbations is relatively small. This initial perturbation means a relatively smaller dispersion of angular momentum of the gas. In the second type, we set a divergence-free Kolmogorov-like velocity perturbation. The amplitude of velocity is subsonic, with σ v ∼ 100 km s −1 , similar to the turbulent velocities in elliptical galaxies (de Plaa et al. 2012;Sanders & Fabian 2013). This perturbation induces a larger root-meansquare angular momentum, but still with a very small net angular momentum. The kinetic energy in the initial velocity perturbation is roughly 1% of the internal energy. We do not continuously drive the turbulence during the computation but only seed the velocity perturbation in the initial condition. In practice, we set the velocity perturbation only on large scales, with wavelength 3.2 kpc < λ < 8 kpc, and let it naturally cascade to smaller scales during the evolution. Note that the nonlinear cascade time ∼ λ/σ v ∼ 10 8 yrs which is comparable to the duration of our simulation.

Gravitational Field
We adopt a spherically symmetric gravitational field with the gravitational acceleration being where M (r), the mass enclosed in radius r, is the sum of the mass of the SMBH, stars, and dark matter, We do not include the self-gravity of the gas. For M87, the Event Horizon Telescope Collaboration et al. (2019) derived a SMBH mass of M • = (6.5 ± 0.7) × 10 9 M , corresponding to a gravitational radius of r g = 3.1 × 10 −4 pc ≈ 0.3 mpc, consistent with the dynamical inferred SMBH mass by Gebhardt & Thomas (2009). For the stars and dark matter, Gebhardt & Thomas (2009) presented radial profiles using observations of stellar light and kinematics. We model the mass of stars using a broken power law for density, similar to the NFW profile (Navarro et al. 1997), which gives a profile of mass as where we set r s = 2 kpc and M s = 3 × 10 11 M . For the dark matter, we also adopt the NFW profile with r DM = 60 kpc and M DM = 3 × 10 13 M . The potential at the inner boundary is softened by where we set s 2 = r 2 a exp −(r/r a ) 2 with r a = 0.5r in .

Cooling and Heating
The cooling function in Equation 3 takes the form of where Λ(T ) is the emissivity.
For T < 10 8.15 K, we use a linear interpolation of the tabulated cooling rate taken from Table 2 of Schure et al. (2009) for solar metallicity. We include this cooling source term by operator splitting forward Euler method.
Despite the presence of radiative cooling, X-ray observations of cluster cores do not show a sharp drop in temperature or the presence of X-ray lines from cooler gas (Vikhlinin et al. 2006), implying that the system is in global thermal equilibrium due to the presence of heating. The source of heating, unlike cooling, is more uncertain and difficult to model analytically. The major source of heating is probably AGN feedback, but with a possible contribution from galaxy mergers, stellar feedback, conductive heating, plasma viscous heating, etc. Similar to Sharma et al. (2012) and Gaspari et al. (2013), we assume an idealized phenomenological heating prescription, in which the global cooling is balanced by global heating with a heating rate approximately equal to the angle-averaged cooling rate at the same radius, where the normalization factor f Q = 0.98 in practice, slightly smaller than unity. This heating prescription provides a way to keep a global equilibrium but not a local equilibrium. Locally, due to thermal instability, some dense gas cools quickly and forms cold clumps and filaments, leading to a two-phase structure of the gas.
We tested the effects of different normalizations of the heating rate. Overall, we find that a higher heating rate tends to suppress the formation of cold gas and therefore decrease the accretion rate. Increasing the heating rate causes the accretion rate to gradually decrease from the pure cooling limit to the adiabatic limit. In fact, if we choose normalization such that the total heating rate equals the cooling rate (f Q = 1), there is little cold gas generated, especially during the later stage of the simulation. The heating prescription in Equation 16 occasionally produces artificial overheating when there is considerable cold gas, leading to the formation of unphysically high-temperature gas. So in practice, we do not include the cooling rate from the cold gas (T < 2×10 6 K) in evaluating the heating in Equation 16. We also artificially turn off the heating function for gas with T > 10 9 K. In the central region, heating plays a less significant role. It is unlikely that heating from AGN feedback operates the same on small scales as it does on larger scales in ellipticals and clusters. For example, jets appear to deposit most of their energy on scales of kpc or greater. For this reason, we elected to turn off heating within 1 pc of the BH, although we find that it does not make a big difference to the results. In practice, we begin heating after 5 Myr and smoothly increase it from zero to the value given in Equation 16 at 10 Myr in order to accelerate the formation of cold gas (similar to Gaspari et al. 2013).

Floors and Flux Corrections
Using Athena++ , we evolve the conservative variables (mass, momentum, and total energy) and treat the source terms by operator splitting. In practice, however, due to the combination of strong shocks, strong rarefaction, and strong cooling over the large dynamic range, the primitive variables of ρ and P can occasionally reach unphysical values, e.g., they can be negative. In addition, accurate integration of the cooling term generally requires a time step much smaller than the hydrodynamic time step, imposing significant restrictions on the computational speed. Lastly, the cold gas can readily cool to temperatures of 10 4 K, making the scale height of the disk too small to resolve, especially at small radii in the simulations. Thus instead of imposing a cooling time step limit, we apply floors when calculating the primitive variables from the conserved variables. For density, we apply both a fixed density floor of n floor = 10 −3 cm −3 and a radius-dependent density floor of 10 −3 (r/1 kpc) −1 cm −3 . For temperature, we apply a floor of T floor = 2 × 10 5 K. The temperature floors help to resolve the vertical structure of the cold disk so that the dynamical evolution within the disk can be more physically meaningful. In this way, we accelerate the computation significantly and ensure that the model remains stable. We tested the effects of floors by running simulations with a higher T floor . The density of the cold gas is smaller, but most results, such as accretion rate, are still similar.
Furthermore, we adopt the first-order flux correction algorithm (described by Lemaster & Stone 2009) to update problematic cells with unphysically large temperatures or velocities. In this algorithm, we first check the cells after the full update to identify the problematic cells in which the density is smaller than the density floors or the density ratio of adjacent cells is larger than a certain criterion (∼ 10 2 in practice). Then we roll back, replace the higher-order numerical fluxes at the interface of the problematic cells with spatially first-order fluxes computed using the Local Lax Friedrichs (LLF) solver and update the cells. This helps to improve the stability of the algorithm and reduces (though does not eliminate) the need for floors. In practice, we find the properties of accretion flow in the normal cells are not affected with or without this method. In general, the correction is required in less than 0.05% of the cells.
The parameters of the simulations and the device type for running them are summarized in Table 1. Note-The models of black hole accretion in this work. The bold simulations are the fiducial models.
(1) In the model labels, C stands for cooling, H stands for heating, D stands for density perturbation, V stands for velocity perturbation, and A stands for adiabatic. The digits following the character for the initial conditions indicate the number of refinement levels used in the model. (11) Restart time for the zoom-in simulation from the larger scale simulation. (12) The type of device on which we run the simulation. a The gravitational radius of the SMBH in M87 is rg = 0.31 mpc.
b A temperature floor of T floor = 10 4 K is used.
c These models are presented in the Appendix.
In this section, we present the results of our fiducial suite of simulations, which include optically thin radiative cooling, heating, and turbulence. In general, we find there are two-phase (hot and cold) accretion flows. The cold accretion flows show two distinct stages that emerge from the simulations: a disky accretion flow and a chaotic accretion flow. Most of the time the disky accretion flow is a better description of the cold gas dynamics. Figure 2 shows a volume rendering and zoomin of the central ∼ 1 kpc region for both of these two stages. The two snapshots are deliberately chosen to represent the typical structure of the chaotic (left) and disky (right) flows. We first describe in detail the properties of these accretion flows. Then we compare them with models with different levels of turbulence and idealized models with no cooling. Finally we present convergence tests.
Throughout the remainder of the paper, we define the gas with T < T cold ≡ 10T floor as cold gas, gas with T ≥ T hot ≡ 0.3T ini (r) as hot gas (where T ini (r) is the Figure 2. Three-dimensional rendering over a 2 × 2 kpc scale region centered on the black hole of chaotic cold accretion (left panel) and cold disk accretion (right panel) emerging in our simulations. This rendering is made using twelve "layers" evenly spaced logarithmically in number density between 10 −3 and 10 3 cm −3 . A zoom-in of a 200 × 200 pc scale region is shown in the bottom right in both images. Usually, a torus structure with size ranging from 100 pc to nearly 1 kpc is present, depending on the angular momentum. Even with chaotic accretion, there is generally a small torus-like structure within ∼ 10 pc. Occasionally, due to collisions of cold clouds and filaments, the torus disappears. initial equilibrium solution), and gas in between (T cold ≤ T < T hot ) as the intermediate gas.

Fiducial Suite: with Radiative Cooling, Heating, and Density Perturbation
In the fiducial suite of models, we include optically thin radiative cooling, radial shell averaged heating, and isobaric density perturbations (Models CHD10, CHD15..., and CHD20... in Table 1). With radiative cooling and heating producing approximate global thermal balance, the initial perturbations serve as the seed of local thermal instabilities, leading to the formation of a two-phase medium in about one cooling time.
We first evolve the large scale simulations for t evo = 300 Myr with a inner radius of r in = 30 pc. We pick several representative time points (30, 90, and 200 Myr, for CHD15-t030, CHD15-t090, and CHD15-t200, respectively) and restart the simulation with a smaller inner radius of r in = 1 pc, and run the simulations for t evo = 2 Myr with t trans = 1 Myr. As we shall describe below, the accretion flow at 90 Myr and 200 Myr is more disky, similar to a typical accretion disk, and accretion at 30 Myr is more chaotic without the presence of a torus or disk. Then we again pick representative times for the middle scale simulations (Models CHD15...), restart them using an even smaller in-ner radius r in = 30 mpc (≈ 100r g ) and run them for t evo = 20 kyr with t trans = 10 kyr (Models CHD20...).
In Figures 3 and 4, we present the two-dimensional z = 0 slices of density, temperature, and the Bernoulli parameter (E + P )/(ρΦ) − 1 with streamlines and their zoom-in from Models CHD20-t200 and CHD20-t030, respectively. Generally, a cold dense torus or accretion disk (similar to that shown in Figures 3) is formed during most of the time of our simulations. Depending on the level of turbulence, the size and orientation of the disk change frequently due to the continuous inflow of filaments and clouds of cool gas. The cold torus resembles a typical cold disk in equilibrium. The size of the torus is typically ∼ 200 pc, but can be larger or smaller, varying from tens to hundreds of parsec. At the same time, the orientation of the disk can differ greatly at different scales. The warps in the disk could be important for explaining the difference between the gas-based SMBH mass and the star/EHT-based mass in M87, which we will discuss in detail in Section 4.
Occasionally ( 10% of the time) the accretion flow is very chaotic and there is no disk structure at a larger scale (10 − 1000 pc). However, if we run a zoom-in simulation during those chaotic phases, we often find a tiny disk at smaller scales. This phase is shown in Figure 4.
1000 500 0 500 1000 . An example slice of density (top), temperature (middle), the Bernoulli parameter (E + P )/(ρΦ) − 1 (bottom), and streamlines of gas velocity through the z = 0 plane for the disk phase in the simulation with cooling, heating, and density perturbations (Model CHD20-t200, t = 20 kyr). From left to right are zoom-in slices covering 1 kpc, 30 pc, and 1 pc, respectively. There is a cold dense torus of size r ∼ 300 pc stretching down to the sink cell at very small scale (∼ 0.03 pc). Note the direction of the disk can differ at different scales. The hot gas is chaotic at the large scale but more laminar at the small scale.

Time Evolution
The top panel of Figure 5 shows the time evolution of mass accretion rate of Model CHD10. After about one cooling time (∼ 10 Myr), the accretion becomes chaotic and stochastic everywhere, with the mass accretion rate through the inner boundary r in = 30 pc varying on timescales of ∼ 1 Myr in a large range between the two simulations with spherically symmetric initial conditions: 0.1 M yr −1 for adiabatic (Bondi) accretion (Model A10) and 10 M yr −1 for accretion with cooling but no heating (Model C10), which is essentially set by mass cooling rate. The time-averaged accretion rate is approximately 1 M yr −1 during the simulation. The accretion rate contributed by the hot gas is around the Bondi rate, though it also varies modestly in time.
The mass-weighted angle-averaged value of the specific angular momentum, 1000 500 0 500 1000 suggesting that these flips in angular momentum might still be important even if we included MHD turbulence for angular momentum transport.
The mass and angular momentum are dominated by the cool disk in the inner region. The mass of cold gas within 100 pc varies between ∼ 10 6 M and 10 8 M (bottom panel of Figure 5). The contribution to the accretion rate from hot gas and cold gas varies in different cases, depending on the exact time and radius that we consider. In the long run, the time-averaged mass accretion rate may be set by the time-averaged rate at Cold, r < 100 pc Hot, r < 100 pc which cold gas is formed,Ṁ cool 1 M yr −1 since all the cold gas would finally be accreted in the absence of star formation or feedback. But most of the time the accretion rate is smaller than the time-averaged value since the formation of the disk or torus suppresses the accretion rate. The peaks of accretion rates close to the spherically symmetric cooling case ( 10 M yr −1 ) are often associated with the rapid change of angular momentum of the cold gas. Due to the cancellation of angular momentum, a considerable fraction of the cold gas can accrete to smaller radii. Figure 6 plots the same quantities as in Figure 5 but for Models CHD...-t030 and CHD...-t200. In the disk stage (t = 200 Myr, red lines), when we trace the accretion flow down to a smaller region by gradually shrinking the inner radius, the accretion rate decreases significantly. The contribution from cold gas drops considerably and remains small during most of the simulation, though there are some peaks in the accretion rate associated with the infall of small cold clouds. The hot gas accretion rate also decreases, but essentially follows a power lawṀ hot ∝ r 0.5 in (marked by the dotted line). We discuss the possible origin of this power-law in Section 4. The angular momentum and mass around the inner boundary is still dominated by the cold gas (middle and bottom panels). The angular momentum of the gas through the inner boundary is similar to the Keplerian value due to the motion of the cold disk. However, the direction of angular momentum in the inner region can differ greatly from that in the outer region. In the chaotic stage (t = 30 Myr, yellow lines), the accretion rate remains high at ∼ 1 pc, but drops significantly at ∼ 0.03 pc. This is associated with the formation of a small torus, as we will show below.  Cold, r < 1 pc Hot, r < 1 pc in . The gas is nearly Keplerian near the sink cell with angular momentum changing frequently. The chaotic accretion and the disky accretion show significant differences. The accretion rate in chaotic accretion tends to be higher than the disky cases. There is more cold gas on small scales for the chaotic cases.

Dynamics
As is shown in Figure 3, the gas appears in two distinct thermal phases: cold and hot phases with apparent boundaries. The cold gas is dense, remaining at a temperature similar to T floor , while the diffuse hot gas roughly keeps the virial temperature (the initial equilibrium temperature T ini ). The cold gas forms a torus within ∼ 300 pc. In Figure 4, the cold gas is more chaotic, but we still find a small torus on parsec scales (nearly face-on by coincidence). Generally, a small disk is formed easily in most cases, with sizes varying from tens of parsecs to hundreds of parsecs. When the disk is stable, the accretion rate tends to be low, similar to adiabatic accretion (Model A10). When some newlyformed cold clouds or filaments fall down to the central region, due to the cancellation of angular momentum, a considerable amount of gas flows into the sink cell. The accretion rate then increases significantly, approaching the spherically symmetric cooling case (Model C10). The direction of the disk in the inner region can be very different from the outer region. The hot flow is more chaotic at large scale but more laminar at small scale.
In the suite of models with initial density perturbation, the initial angular momentum dispersion is very small. The flow tends to be very similar to the cold chaotic accretion model Gaspari et al. (2013) on larger scales. However, a large torus still forms intermittently. On small scales, a disk is present most of the time.
When there is a disk, the scale height of the cold disk is related to the temperature (sound speed) of the gas by Thus it is extremely difficult to resolve the disk when the temperature is low and the radius is small. We set a temperature floor of 2 × 10 5 K, higher than the physical value of 10 4 K, so as to resolve the structure of cold 10 2 10 1 10 0 10 1 10 2 10 3 10 4 10 5 10 6 10 7 Time ratio , and CHD20-t200 (r < 1 pc, 20-30 kyr). The lighter lines show the variances of angle-averaged density and temperature within the averaging time range. The emission-weighted temperature is also shown in purple, where we use a range of 0.5 − 7 keV. For the hot gas, the variables follow power law in radius of ρ ∝ r −1 , T ∝ r −1 ,Ṁ ∝ r 1/2 . The cold gas is isothermal with a higher density and large inflow and outflow rate. The emission-weighted temperature profile is flat over a wide range of radii, implying that at small radii, the gas that dominates ∼ keV emission is not the hot viral gas but rather the intermediate gas with T ∼ keV.
gas well within the simulation domain. A more physical cold disk would be denser and retain a smaller scale height by a factor of ∼ 5. But the accretion rate and circularization of the cold gas should be similar, since they are basically determined by the angular momentum and inflow motion of the cold gas in the simulations, regardless of T floor . Going to small radii, the cold gas accretion is likely subdominant, so the impact of T floor might be relatively small in this case. Within ∼ 0.3 pc, the disk is aligned with the grid, which is perhaps due to the low resolution compared with the disk scale height. More generally, the formation of a torus is likely robust but the detailed thermodynamics, cooling, and heating properties of the disk might not be very reliable, especially on small scales. It remains to be investigated in the future about the exact effects of floors and resolutions on the accretion flow.

Radial Profiles
The time-and angle-averaged radial profiles of the gas and its two components (hot and cold) are shown in Figures 7 and 8. Throughout this work, we define the time and angle average of a variable A as where S is the area where we count the variable (the place of total, cold, and hot gas). The mass-weighted average is defined similarly by 10 4 pc 10 3 10 2 10 1 10 0 10 1 10 2 10 0 10 2 10 4 L [kpc km s 1 ] Figure 9. Mass distribution of gas density, temperature, and angular momentum at different radii of the disk stage (Model CHD20-t200). Dashed lines are initial density, initial temperature, and Keplerian angular momentum in the corresponding radius. There is a clearly bimodal distribution of density and temperature within ∼ 1 kpc, indicating that the gas is two-phase. The distribution of specific angular momentum has a long tail.
Note that we do not directly obtain a global steady solution over the whole simulation domain shown in Figures 7 and 8. This would require 10 4 times more computational time, which is unrealistic in the foreseeable future. Instead, the solution at r < 10 pc is essentially the steady solution with r in = 0.03 pc (CHD20-t200; which is run for 20 kyr), the solution at 10 < r < 1000 pc is the solution with r in = 1 pc (CHD15-t200; which is run for 2 Myr), and the solution with r > 1000 pc is the solution with r in = 30 pc (CHD10; which is run for 300 Myr).
As is shown in Figure 7, cold gas is generated within ∼ 2 kpc. The temperature profiles of the hot gas follow a T ∝ r −1 scaling, as expected from conservation of energy, while the cold gas remains at the floor value of ∼ 2 × 10 5 K. By contrast, the X-ray emission weighted temperature profile is rather flat (we use an emission range of 0.5 − 7 keV). This implies that the gas that dominates ∼ keV emission is not hot gas at the viral 10 1 10 0 (v 2 R + v 2 z ) 1/2 /v Kep (100 pc) Figure 10. Snapshot of azimuthally-averaged temperature (left) and density (right) for disky stage (Model CHD20-t200) in a coordinate system defined by the angular momentum of the cold gas. The velocity is normalized by the Keplerian velocity at 100 pc. The accretion flow is laminar with smooth streamlines. The accretion is mostly from the polar region and the hot gas both accretes to smaller radii and condenses onto the surface of the cold disk.  Figure 11. Angular profile of the azimuthally-averaged gas density (top), temperature, (middle), and rotational velocity for the disk stage (Model CHD20-t200) in a coordinate system defined by the angular momentum of the cold gas. There is a clear cold dense Keplerian disk at ∼ 100 pc. temperature but rather the gas with T ∼ keV, over a wide range of radii. Note that the somewhat cooler gas with T ∼ keV < T vir has higher density (at roughly constant pressure) which is why it dominates the keV luminosity. Observations of M87 (Russell et al. 2015) suggested there are two significant temperature components within 2 kpc and did not find evidence for a temperature increase within ∼ 0.25 kpc. The flat temperature profile is also reported in observations of the galaxy M84 (Bambic et al. 2022, in prep). This is consistent with our simulations.
The angle-integrated mass inflow and outflow accretion rates are defined aṡ andṀ We also define the net accretion rate byṀ (r) ≡ M in (r) −Ṁ out (r). For cold gas, both the gas inflow and outflow rates are large, but the net accretion rate through the inner boundary is relatively small. For the hot gas, the inflow rate follows a r 1/2 scaling. However there is no strong hot outflow, so the net accretion rate also follows a radial dependence of r 1/2 , rather than a constant. This could imply an accumulation of hot gas but as is shown in Figure 5, the hot gas mass does not change. So theṀ ∝ r 1/2 scaling implies that the hot gas is mixed with the cold gas. The density profile shows a scaling of ρ ∝ r −1 in Figure 7, which is different from the expected power law of r −3/2 for γ = 5/3 if the flow were adiabatic and spherically symmetric. This is consistent with scaling of hot gas accretion rate, sinceṀ ∼ ρvr 2 ∝ ρr 3/2 . The cooling time in Figure 7 shows that t cool /t ff ≡ E int / Q − /t ff > 10 for hot gas at 1 pc r 1 kpc. However, we note that the low-temperature tail of the hot gas can still cool efficiently due to radiative cooling and mixing with the cold gas. Actually, we find t cool /t ff ∼ 10 for the intermediate gas in this region (though not shown here).
The Bernoulli parameter in Figure 8 shows that, in contrast to the kpc scale, the hot gas is, on average, only slightly unbound inside 100 pc. As is shown by the snapshots in Figure 3, part of hot gas is marginally bound within 1 pc. The cold gas is bound and dominated by rotational energy. So overall, the hot gas is unbound in most of the volume, but some hot gas is slightly bound within ∼ 1 pc. The cold gas is completely bound with energy dominated by Keplerian rotation.
The radial profile of mass-weighted time-and angleaveraged absolute value of the specific angular momentum |L| and the three components, L x , L y , and L z , are shown in Figure 8. The hot gas is overall sub-Keplerian,  varies. Due to the lack of viscosity, inefficient transport of angular momentum leads to a small accretion rate from cold gas. The accretion rate is mainly from the hot gas. It is unclear whether magnetic fields will change the relative hot and cold gas accretion significantly; we plan to investigate this in future work. The radial dependence of net accretion rate shown in Figure 7 is consistent with the cold gas mass accumulation in a large region shown in the bottom panel of Figure 5. This implies that we may not be in a statistically steady state. More time averaging may be needed, though the current averaging time is already 10 4 times of the free-fall time at the inner radius. The accumulation of cold gas more likely reflects the fact that we do not model angular momentum transport in the cold disk, which likely precludes the existence of a global statistical steady state.

Phase Distribution of Hot Gas and Cold Gas
In analyzing the radial structure, we divided the gas into hot and cold phases. Here we discuss the mass distribution of the gas phases and the validity of the two-phase description. In Figure 9, we show the mass distribution of density, temperature, and angular momentum for different radii. There is clearly a bimodal distribution of density and temperature within ∼ 1 kpc, indicating the two-phase structure of the gas. The cold gas is generated at around 2 kpc and the two-phase medium extends to the central region. The mass of cold gas is estimated to be ∼ 10 6 − 10 8 M in the whole simulation domain. The hot gas has a virial temperature and density slightly lower than the initial solution. Despite the bimodal distribution, there is a considerable amount of intermediate gas between T cold and T hot . This broad temperature range could be important in interpreting and modeling ∼ keV X-ray emission (as indicated by the emission-weighted temperature in Figure 7).
A Keplerian cold disk is assembled around the black hole. The angular momentum distribution in Figure 9 reflects the sub-Keplerian hot gas and nearly Keplerian cold gas. For the hot gas, the tail of the angular momentum distribution follows dM/d log 10 L ∝ L 2 or equivalently, dM/dL ∝ L. in . The angular momentum approaches Keplerian, but it is still sub-Keplerian with a relatively large dispersion. There is little difference between simulations with density and velocity perturbations.

The Accretion Disk
In the typical disky stage, there is generally a welldefined disk or torus structure at smaller radii. To analyze the properties of the disk, it is convenient to define a new coordinate system (r, θ , φ ) in which the zdirection is aligned with the radial-and angle-averaged angular momentum of the cold gas.
We plot the φ -averaged density and temperature snapshots in Figures 10 and the angular profiles in 11. The results indicate a cold dense disk structure at ∼ 10 − 100 pc scale. The hot flow is laminar and falls down to the central sink region from the polar region or merges onto the cold disk by condensation. The cold disk closely follows Keplerian motion with a very small radial component. Most of the accretion is from hot gas in the polar region. The orientation of the disk changes within 10 pc so the inner disk is not described well in this coordinate system.
To investigate whether the cold disk would fragment into clumps and form stars due to spiral-mode gravitational instability, we calculate the Toomre parameter (Toomre 1964 where Σ is the surface density. Given the temperature floor of T floor = 2 × 10 5 K we use now, we find that the Toomre parameter of the disk reaches a minimum value of ∼ 20 at ∼ 100 pc. For T floor = 10 4 K, we would obtain Q min ∼ 4 if we suppose the surface density does not change, which is an appropriate assumption since the mass of the cold gas at ∼ 100 pc is basically determined by the infall and accumulation of the cold gas, regardless of T floor . In this case, the disk is marginally stable, but might be unstable if the surface density is higher or the temperature is even lower, which would lead to efficient fragmentation and star formation. In reality, in many cases the gas may be able to cool to below 10 4 K, i.e., to become molecular ∼ 100 K. This is observed in many ellipticals and groups and clusters (Davis et al. 2013(Davis et al. , 2017, though probably not in M87. This implies that the gas may fragment and form stars in some cases and that star formation and supernovae feedback could be important to include.

Chaotic Stage
The chaotic stage is obtained only in relatively rare cases ( 10% of the time) but is generally associated with a high accretion rate of cold gas. When a large amount of cold gas with negligible angular momentum is falling into the central region, the accretion flow is chaotic at smaller radii. The simulation at 30 Myr (Model CHD15-t030 and CHD20-t030) is a typical cold chaotic accretion flow, as is shown in Figures 4, 5, and 6. Here we discuss the chaotic case in more detail.
The radial profiles for the chaotic stage are shown in Figures 12 and 13  the accretion rate decreases by one order of magnitude. The disk direction within 0.3 pc is different from the torus at ∼ 1 pc. Similar to the disky stage, the disk is aligned with the grid at small scales, which is also possibly due to the low resolution compared with the disk scale height. The hot gas is more chaotic compared with the disky stage, with a strong hot outflow on small scale. Due to the significant cooling, the heating rate around 10 pc is also slightly higher, leading to a relatively higher temperature and Bernoulli parameter for the hot gas, which is possibly associated with the strong hot outflow. The density of hot gas is also higher, possibly due to intense mixing with the cold clouds and filaments. The emission-weighted temperature profile is still relatively flat over a large range. On smaller scales ( 1 pc), a cold disk forms and the morphology and radial profiles tend to be similar to the disky cases. Gaspari et al. (2013) proposed that chaotic cold accretion could boost the accretion rate by up to two orders of magnitude. Here we find that while the accretion is overall chaotic and turbulent on large scales, a cold dense disk is formed easily in most of our simulations, especially on small scales ( 100 pc). Even for the very chaotic stage, there is usually a small disk within 3 pc. The formation of the cold disk reduces the accretion rate toward the central region most of the time. We will discuss this in more detail in Section 4.

Accretion with Larger Angular Momentum
Here we discuss the models in which we add velocity perturbations to the initial conditions (Models CHV... , and momentum flux for the adiabatic simulations with density (blue) and velocity (green) perturbations. We average between 100-200 Myr for r > 30 pc, 1.5-2.5 Myr for 1 < r < 300 pc, and 10-20 kyr for r < 10 pc. The scalings are ρ ∝ r −1 ,Ṁin ∝ r 1/2 , anḋ pin ∝ r 0 .
in Table 1). These perturbations imply a larger range in the specific angular momentum |L|, although the net angular momentum L remains small. The resulting flows are generally similar to the fiducial suite, except that the size of the disk tends to be larger, typically 500 pc, The orientation and size of the torus are changed by the inflow of new cold gas with different angular momentum, but this occurs infrequently. The typical time-scale for this process is around tens of million to hundreds of million years.
The time evolution for Models CHV... is shown in Figure 14. The accretion around 100 Myr is suppressed because most of the cold gas is accumulated in the disk,  contributing little to the accretion rate. The inflow of cold clouds leads to an increase in the accretion rate at later times. Overall the accretion flows show similar behavior to the disk stage of the fiducial suite, except that the mass of cold gas is larger. When we trace the accretion flow down to small scales, the accretion rate follows a similar radial scaling as the fiducial suite. Accretion is also dominated by the hot gas on small scales. Although not shown here, we find that the flow morphology, radial profiles, and mass distribution show similar properties to the disky stage of the fiducial suite.

Adiabatic Accretion with Turbulence
In the fiducial suite of simulations, despite the presence of significant cooling, the cold disk does not contribute significantly to the accretion rate on small scales. Thus it is of interest to investigate the properties of adiabatic accretion with no cooling, so that cold gas is absent. Here we present adiabatic simulations with density and velocity perturbations (Models AD... and AV... in Table 1) and compare them with the fiducial suite that has cooling and heating and generates a multiphase medium.
We plot the time evolution of mass accretion rate and angular momentum in Figure 15. The accretion rate fluctuates, but overall is similar to the spherically symmetric case at ∼ 30 pc. The absolute value of specific angular momentum through the inner boundary is sub-Keplerian and small compared with the simulations with cooling. The net angular momentum L 2 x + L 2 y + L 2 z is even smaller, implying a relatively large dispersion ( |L| 2 − L 2 x − L 2 y − L 2 z ). As we restart the simulations with a smaller r in , the accretion rate follows the similar power law r 1/2 in as in the hot gas in the fiducial suite. Towards the inner region, the angular momentum also approaches Keplerian gradually, but it is still sub-Keplerian with a relatively large dispersion and never forms a Keplerian disk. As is shown by the snapshots plotted in Figure 16, though the fluctuations in the density and temperature are relatively small, the velocity field is highly turbulent. Most of the gas is unbound except at small scales. Despite many similarities with the fiducial suite, the streamlines in the adiabatic simulation are chaotic, remarkably different from the laminar streamlines in the disk stage of the fiducial suite.
We investigate the radial dependence of the accretion flow in the presence of moderate turbulence. The radial profiles are shown in Figure 17. We obtain a steady solution in r < 10 pc with r in = 0.03 pc, a steady solution in 1 < r < 300 pc with r in = 1 pc, and a steady solution in 30 < r < 10 4 pc with r in = 30 pc, as is shown clearly by the net accretion rate. Adiabatic accretion also shows similar scaling to the fiducial cases. The density clearly scales with radius by ρ ∝ r −1 . Combining the velocity scaling v r ∝ r −1/2 , this gives a power law ofṀ in ∝ r 1/2 and constant momentum fluxṗ in ∝ r 0 . In contrast to the disk stage of the fiducial suite, however, there is a strong outflow of hot gas which leads to this scaling. Thus though both fiducial suite and adiabatic simulations showṀ in ∝ r 1/2 , it seems to be for very different reasons. We return to this in Section 4.

Convergence Tests
To examine the convergence with resolution, we present a higher-resolution simulation with a root grid of 256 3 (Model CHD10-r256 in Table 1). In this case, the resolution is doubled everywhere compared with the fiducial simulation (Model CHD10). As is shown in Figure 18, the history of the accretion rate is very similar to the fiducial case. The morphology, dynamics, and radial profiles of the two simulations are also very similar. This further confirms the validity of our results.
To test whether our method of changing the inner boundary affects the outcome, we directly simulate the accretion flow with r in = 1 pc using 12 levels of mesh refinement with the finest resolution of ∆x = 0.12 pc for a long evolution time of ∼ 100 Myr (Model CHD12 in Table 1). This simulation required nearly one million core hours due to the extremely small time step. The results essentially confirm our results from the fiducial suite of simulations. We find that most of the time, the accretion rate at small scales stays low compared with the case with a larger r in . Occasionally the accretion rate is higher than the large-scale simulation due to the accumulation of cold gas. The hot gas accretion rate is also systematically lower than the fiducial case CHD10. In Figure 19, we plot the radial profile of accretion rate in the chaotic stage and compare it with the simulation CHD15-t030. Across the range of radii, the accretion rate only shows slight deviations from the model CHD15-t030. This indicates that the changing inner boundary strategy does not introduce significant systematic deviations at smaller radii.
Interestingly, we find a two-tori structure during the late time of this simulation, as shown by the z = 0 slice of density with streamlines in Figure 20. Before the old torus is depleted, a new torus pointing to a different direction forms on a larger scale due to the infalling of cold clouds with different angular momentum from the old torus. The size of the old torus is ∼ 20 pc and the new torus is ∼ 60 pc. This structure sustained for a long time (∼ 30 Myr). Similar to the warped disk found in the fiducial suite, this also implies that the direction of the accretion disk could differ greatly on different scales, though presumably for various reasons. The potential caveat is the grid alignment of inner torus, which is possibly due to the low resolution (r in /∆x min ≈ 8). Given that each torus is a massive structure, the inclusion of self-gravity may change the evolution of each torus and their interaction qualitatively.

Towards the Event Horizon
In Figure 21, we summarize the accretion rate through the inner boundary at different times in Model CHD as well as Models CHV, AD, and AV. In Table 2, we list the time-averaged mass accretion rate and the contribution from different phases of gas for Model CHD10 and CHD...-t200 (the typical disk stage). The radial dependence in Figure 21 is well described by the r 1/2 in scaling, especially when cold inflow is negligible (i.e., except for the yellow points for Model CHD at t=30Myr). The sink region boundary condition in our simulations means the absence of radial pressure support, leading to a negligible outflow rate through the inner boundary and increasing the net accretion rate. This effect is reasonable only if r in is the event horizon of the black hole. Unfortunately, covering such a large physical region is a formidable challenge due to the small time-steps required.
Fortunately, we can reasonably extrapolate the accretion rate down to smaller radii based on the selfsimilar radial dependence of the accretion flow in the simulations we have conducted. The effect of the inner boundary is to force v r (r in ) ≈ −c in ∝ r −1/2 in . Combining ρ ∝ r −1 , we expectṀ ∝ r 1/2 . As is shown in Figure 21, we find a mass accretion rate of 10 −1 −1 M yr −1 at 30 pc, 10 −2 −10 −1 M yr −1 at 1 pc,  Figure 21. Relationship between mass accretion rate and inner radius for different runs. The accretion rates are values measured with a time range of t = 3 Myr for rin = 30 pc, t = 1 Myr for rin = 1 pc, t = 10 kyr for rin = 30 mpc, and t = 100 yr for rin = 1 mpc (3rg). The simulations with the highest accretion rate (Fiducial, 30 Myr; yellow) are in the cold chaotic accretion phase. Overall, we find a universal power law ofṀ ∝ r 1/2 in over a large dynamic range from Bondi scale to horizon scale, spanning three to five orders of magnitude. The normalization of the accretion rate vs. radius is similar for simulations not in the cold, chaotic accretion phase. and 10 −3 − 10 −2 M yr −1 at 0.03 pc. The exception is the chaotic case (fiducial simulation at 30 Myr, Models CHD...-t030), which is relatively rare, and which produces a larger accretion rate at small radii. By extrapolating our simulation results down to the event horizon, we obtain an accretion rate ∼ 10 −4 −10 −3 M yr −1 . We find a fairly good agreement in the simulations. The accretion rate is a thousand times smaller than the Bondi accretion rate prediction. It is remarkably consistent withṀ ∼ (3 − 20) × 10 −4 M yr −1 from EHT observations (Event Horizon Telescope Collaboration et al.

2021).
As a test of this scaling of the accretion rate with inner boundary radius, we furthermore ran one simulation with r in = 1 mpc (3r g !) and 25 levels of mesh refinement (Model CHD25-t200 in Table 1) to explore the accre- 1.0 pc 2-3 Myr 0.034 M yr −1 36% 24% 40% 35% 7.9% 57% CHD20-t200 0.03 pc 20-30 kyr 0.010 M yr −1 0.1% 28% 72% 0.0% 5.2% 95% Note-The mass accretion rate and the role of different phases of gas for Models CHD10 and CHD..-t200. Here fṀ is the mass fraction contributed by the corresponding gas phase and ft is the time fraction that the corresponding phase dominates the accretion.
tion around the event horizon. This result is shown in Figure 21, and yields a accretion rate ∼ 10 −3 M yr −1 , consistent with the extrapolation. We note that in this case, the resolution is insufficient to capture the dynamics of the cold disk and we neglect crucial physics (general relativity, magnetic field, radiation, collisionless effects, etc) that we plan to include in the future. We also do not aim to give a precise prediction of the accretion rate around the event horizon using this simulation. But nonetheless, it demonstrates that the scalinġ M ∼Ṁ Bondi (r g /r B ) 1/2 might be a good prescription for the accretion rate on horizon scales. The approach used in this paper, when applied to future MHD simulations which we plan to run, will enable us to perform MHD simulations over the entire dynamic range of the accretion flow. This will have the potential to benefit more realistic initial or boundary conditions for GRMHD simulations on horizon scales, instead of an idealized torus with specified angular momentum.

Physics of Hot Gas Accretion
It remains unclear how to self-consistently link the accretion rate and dynamics at the Bondi scale to the event horizon for hot gas accretion. For adiabatic accretion flows, there exist three possible solutions assuming selfsimilarity (Gruzinov 2013): 1. the spherically symmetric Bondi flow, with ρ ∝ r −3/2 and constant mass flux; 2. the Convection Dominated Accretion Flow (CDAF) with ρ ∝ r −1/2 , characterized by strong entropy inversion and constant energy flux withṀ ∝ r; 3. a turbulent accretion flow identified by constant momentum flux with ρ ∝ r −1 andṀ ∝ r 1/2 .
As is discussed in Section 4.1, the third scaling oḟ M ∝ r 1/2 is an excellent description of the accretion flow in our simulations in the presence of moderate turbulence. It is worth noting that there is some evidence that this scaling is universal if the accretion is nearly adiabatic and modestly turbulent. Similar scaling is reported in various contexts and with varying physics, including some MHD and GRMHD simulations of radiatively inefficient black hole accretion (Pang et al. 2011;White et al. 2020), fueling of Sagittarius A* via the stellar winds of Wolf-Rayet stars within the central parsec of the galactic center (Ressler et al. 2018(Ressler et al. , 2020, Bondi-Hoyle-Lyttleton accretion in supergiant X-ray binaries (Xu & Stone 2019), wandering massive black holes in the outskirts of galaxies (Guo et al. 2020), and adiabatic super-Eddington accretion (Hu et al. 2022).
It is somewhat surprising that our simulations with multiphase gas, heating, and cooling produce a similar accretion rate ∝ r 1/2 in to our adiabatic simulations of Bondi accretion with turbulence. The reason this is surprising is that the hot gas dynamics is very different in these two cases. There are no strong outflows of hot gas when cooling produces a significant cold disk. Instead, we find mixing of hot and cool gas that depletes the amount of hot gas at smaller radii. The streamlines of the hot gas with cooling ( Figure 3) and in the adiabatic simulations ( Figure 16) highlight the large differences in the hot gas physics: in the former case the hot gas settles relatively smoothly onto the cold disk while in the latter case the hot gas is much more turbulent. The fundamental explanation in most of the literature for a reduction in the accretion rate relative to the Bondi rate (Ṁ ∝ r p for some 1 > p > 0) is that outflows suppress the accretion rate at small radii (Blandford & Begelman 1999), with some other conserved quantity setting the density profile and thusṀ (r). This is essentially the argument forṀ ∝ r 1/2 orṀ ∝ r in the previous literature (e.g., Narayan et al. 2000;Gruzinov 2013). However, in our simulations with multiphase gas,Ṁ is suppressed but not for this reason. Given this different physics, we do not see any particularly compelling reason to expecṫ M (r) to be similar in our simulations with and without cooling, but yet they seem to be; this is an interesting puzzle.

Physics of Cold Gas Accretion
In our simulations, cold gas accretion includes two distinct stages: the typical disk stage and the relatively rare chaotic stage. We find that while the accretion is overall chaotic and turbulent, there is always a cold dense disk or torus in the inner hundreds of parsec (or, in the most extreme cases, in the inner few pc). However, when it is present the cold disk does not generally contribute significantly to the accretion rate at small radii, which is dominated by the hot gas. This is partially by construction, since we do not include magnetic fields or viscosity in our simulations. As a result, the accretion rate of cold gas typically decreases with decreasing inner radius because the cold gas accumulates in a disk.
Consistent with previous results (Li & Bryan 2012;Gaspari et al. 2013), we find that there can occasionally be a significant accretion rate of cold gas at intermediate radii. The driver of this accretion is the large fluctuations in the angular momentum of the cold gas at large radii and the resulting cancellation of angular momentum via frequent collisions between the inflowing cold gas and the disk. A series of analytical calculations (King & Pringle 2006, 2007 postulated that BHs may grow through a series of randomly oriented accretion events analogous to what we find in our simulations. It is possible that the random fluctuations in inflowing angular momentum that promote occasional cold gas accretion may still be important even in the presence of turbulence and magnetic fields in the cold gas disk. This is because the viscous time in a thin disk is relatively long compared to the time over which the inflowing angular momentum changes in our simulations. On the other hand, large scale magnetic stresses can remove angular momentum from the inflowing cold gas much more efficiently than local turbulence; we will study the role of these large-scale magnetic stress in future work. One important feature of the cold disk in our simulations is that the orientation of the disk can differ significantly on different scales. Furthermore, we find the co-existence of two tori with different sizes and orientations in one simulation covering a large radial dynamic range (Model CHD12 in Table 1). The warps that we find could be important for explaining the difference between the gas-based SMBH mass in M87 of 3 × 10 9 M and the star/EHT-based mass of 6×10 9 M . The reason is that the gas-based mass inference is sensitive to inclination (e.g., Walsh et al. 2013). If there is a warp and the inclination varies on the scale of the observed disk, that could bias the mass inferences from gas kinematics.

Accretion of Hot vs. Cold Gas
AGN feedback is believed to play a key role in galaxy evolution and is crucial for cosmological simulations to produce a realistic population of massive galaxies, groups, and galaxy clusters (Fabian 2012b;Heckman & Best 2014;Somerville & Davé 2015). Our work has the potential to provide a more physical and accurate subgrid model of AGN feedback in groups and clusters by linking the accretion rate at a larger scale to the accretion rate on horizon scales. In cosmological simulations the Bondi accretion rate is typically utilized to model how the black hole accretes gas from its surroundings (Somerville & Davé 2015;Pillepich et al. 2018). By contrast, accretion of cold gas is the key in current highresolution simulations of how AGN feedback balances cooling in clusters (Li & Bryan 2012;Gaspari et al. 2013, etc). The latter models are essentially a limit cycle in which cold gas accretion generates feedback that suppresses cooling for a period of time (a few cooling times of hot gas), then a new burst of cold gas accretion sets in. We find, by contrast, that the accretion at small radii is dominated by the hot gas most of the time (see Table 2). In future work it will be critical to understand how the interplay between hot gas and cold gas accretion changes with the inclusion of magnetic fields, which can further drive accretion of the cold gas. This will clarify whether cold or hot gas dominates the time-averaged accretion and feedback in massive galaxies. Thermal conduction and viscosity may also be important because they can alter the dynamics of the hot gas and the mixing between cool and hot phases that is important in setting the hot gas accretion rate in our simulations (Nayakshin 2004).

Comparison to Previous Work
The results presented here on the interplay between heating, cooling, and accretion in massive galaxies and galaxy clusters are closely related to previous simulations by Li & Bryan (2012) and Gaspari et al. (2013) (among others). Relative to these earlier works, our main improvements include: • We focus on connecting physics on multiple scales, adopting a very small inner boundary radius, generally spanning six orders of magnitude down to 0.03 pc, with one test case down to ∼ 3r g .
• We achieve a high resolution even around the sink region, which is crucial in resolving the cold gas.
• We run the simulation for a very long time to investigate the statistically steady state. Li & Bryan (2012) studied the accretion of gas with modest rotation from Mpc scales down to pc scales. They found that the cooling gas rapidly forms a thin disk, as we do, but did not trace the evolution of the disk for a long time. Gaspari et al. (2013) carried out simulations linking the physics at tens of kpc in clusters to sub-pc scales. Overall, they found that chaotic cold accretion was more prominent/important than we find in our simulations. Part of this could be due to the relatively short duration of their simulations, which was ∼ 40 Myr, only a few t cool . We find that chaotic accretion is more prominent at early times comparable to these in our simulations, while the disk phase dominates later as the angular momentum of the cold gas continues to accumulate. Another difference is that we only include turbulence (which generates angular momentum fluctuations) in the initial conditions. Gaspari et al. (2013) adopted continuous turbulence driving, which might promote chaotic accretion by changing the angular momentum of the inflowing gas. Simulations with more realistic turbulence driven by AGN, substructure, and/or stellar feedback would be valuable to better understand the cold gas dynamics on small scales in massive galaxies.

SUMMARY
We have presented a series of high-resolution, threedimensional hydrodynamic simulations of the fueling of SMBHs in elliptical galaxies, taking M87 as a typical case. The key physical ingredients are radiative cooling and a phenomenological heating model. The simulations span more than 6 orders of magnitude in radius, reaching all the way down to sub-parsec scales with one test simulation reaching to near the event horizon.
Our main conclusions are: • The accretion flow takes the form of multiphase gas at radii less than about a kpc. Most of the time accretion at the smallest radii is dominated by hot virialized gas. Nonetheless, simple adiabatic simulations of hot gas accretion do not reproduce the dynamics found in our full simulations: the hot gas inflow rate in the latter is strongly influenced by condensation onto the cold disk formed by radiative cooling.
• Cold gas accretion includes two dynamically distinct stages: the typical disk stage in which the cold gas resides in a rotationally supported disk and relatively rare chaotic stages in which the disk is broken up and the cold gas inflows via chaotic streams. Disks are formed easily in most of our simulations, especially on small scales. Chaotic cold accretion is obtained only in a fine-tuned parameter space.
• The mass accretion rate scales with radius aṡ M net ∝ r 1/2 in when hot gas dominates. By extrapolating the accretion rate down to the event horizon, we obtainṀ 10 −4 − 10 −3 M yr −1 , similar to (3−20)×10 −4 M yr −1 from EHT observations. The same scaling ofṀ ∝ r 1/2 can explain both observations and simulations of the Galactic Center (Ressler et al. 2020). This prescriptionṀ ∝ r 1/2 is thus an attractive option for connecting hot gas properties at large radii (in observations or cosmological simulations) to the accretion rate near the event horizon. We suggest that this is a more physical and accurate subgrid model of SMBH fueling by hot gas than the standard Bondi accretion rate assumption typically used in large-scale cosmological simulations.
• The relative contribution of hot and cold gas to the total accretion rate varies on different spatial scales and on different time-scales. Most of the time, accretion at small radii is primarily due to hot gas, and only in rare periods does cold gas dominate. It is possible, however, that cold gas accretion could still dominate the time-averaged accretion rate, but simulations with more physics in the cold gas (e.g., magnetic fields) are needed to quantitatively determine this.
• The orientation of the gas disk can differ significantly on different spatial scales, primarily due to the inflow of gas with a wide range of angular momentum. Such warps in the disk could be important for explaining the difference between the gas-based SMBH mass in M87 of 3 × 10 9 M and the star/EHT-based mass of 6 × 10 9 M (because the gas-based SMBH mass is sensitive to the inclination of the gas disk on scales that are not well resolved observationally, and thus to the presence of warps).
Because of insufficient resolution and the absence of crucial physics (e.g., magnetic fields), our results on the internal dynamics of the cold gas are likely less reliable than our results on the dynamics of the hot gas; inclusion of anisotropic conduction and viscosity would further bolster the realism of our treatment of the hot gas. MHD simulations with these more realistic physics ingredients will be carried out in the near future. Using these results, we can also provide more realistic initial and boundary conditions for GRMHD simulations of accretion on event horizon scales.

A.1. Adiabatic Spherical Accretion
First, we present our tests of adiabatic spherical accretion (Bondi accretion). For a point source, the Bondi accretion rate isṀ where c s,∞ is the sound speed at infinity, ρ ∞ the density at infinity, and the factor q γ (γ) = 1 4 2 5 − 3γ with q γ = 1/4 for γ = 5/3. At small radii, we have a self-similar solution, In this study, if the gas is adiabatic, the accretion is similar to a Bondi accretion, but the potential is from the mass of a point source (the SMBH) as well as the stars and the dark matter. In addition, there is no fixed sound speed and density at infinity. The accretion rate is a function of the density and temperature of the gas, which could be different at different radii. However, due to the assumption of a flat entropy core (ρ ∝ c 2/(γ−1) s = c 3 s ), the Bondi rate is constant within the core radius. In this test (Model A13), we employ 13 levels of mesh refinement and set r in = 0.3 pc. We evolve the simulation for ∼ 30 Myr, which is about the dynamic time at 10 kpc, 100 Bondi time, and 10 6 times the free-fall time at the inner boundary of 0.3 pc. The radial profiles of density, temperature, and mass accretion rate are shown by dashed lines in Figure 22. The accretion rate is similar to the Bondi accretion rate ∼ 0.05 M yr −1 within ∼ 300 pc. The radial profiles are very stable, following scaling of ρ ∝ r −3/2 , T ∝ r −1 inside Bondi radius. The profiles are nearly steady over the simulation, without significant change. The density and temperature profiles are very close to the initial condition. The inflow velocities are subsonic over the simulation domain. The flat entropy core is maintained during the simulation. The snapshots of density and temperature are shown in Figure 23. There are only tiny (nearly invisible) asymmetries due to the Cartesian mesh. This test builds a steady accretion flow similar to the Bondi solution in a large domain from ∼ 20 kpc down to the sub-parsec scale spanning 5 orders of magnitude.

A.2. Cooling Spherical Accretion
Radiative cooling drastically changes the accretion flow when adiabaticity is violated. Though the free-fall time of the gas is smaller than the gas cooling time by a factor of 5 or higher, the ratio of cooling time and gas inflow time can be smaller. As is shown in the adiabatic case, the gas cooling time is much shorter than the inflow time at ∼ 1 kpc. In this test (Model C13), we include radiative cooling and set a temperature floor of T floor = 10 4 K. The floor is different from the fiducial suite because there is no disk to resolve. Thus we can safely set a more physical temperature. Similar to the adiabatic case, we also employ 13 levels of mesh refinement and set r in = 1 pc.
The radial profiles are plotted by solid lines in Figure 22. With cooling, the accretion rate increases exponentially and finally saturates at ∼ 20 Myr. It is enhanced by a factor of ∼ 100, approaching ∼ 10 M yr −1 . 10 1 10 0 10 1 10 2 10 3 10 4 n [cm 3 ] 10 5 10 6 10 7 10 8 10 9 T [K] The gas within ∼ 100 pc cools down to the temperature floor of 10 4 K in ∼ 20 Myr. There is a cooling flow outside ∼ 300 pc with t cool ∼ t inflow . The gas density within ∼ 1 kpc increases by up to two orders of magnitude, following a power law of r −3/2 in the whole simulation domain. The cold gas is quasi-free-fall within ∼ 100 pc. The sonic point is at ∼ 300 pc. The compressional heating does not influence the accretion flow. For the hot accretion flow, the time ratio t cool /t ff increases towards smaller radii. But if cooling develops in about one cooling time, it changes the boundary condition substantially and thus we have t cool /t ff 1. The Bondi solution is not appropriate in this circumstance.
The snapshots in Figure 23 show the condensation of the galactic core. The collapse is symmetric and monolithic. We do not find thermal instabilities present, though there are tiny asymmetries. The giant core continues to expand slowly in size. The flow becomes supersonic at the boundary of the cold nucleus and the hot ambient gas, consistent with the region where the gas is in a fast cooling regime (Λ(T ) ∝ T −1/2 , T 10 6 K). Even if the initial sonic point is within the Bondi radius, the mass accretion rate would be so large pushing the sonic point r s > r B . The solutions are in good numerical stability with very low oscillations both in space and in time. We note that there are still some small asymmetries in smaller radii.

A.3. Turbulent Cold Accretion
It is worthwhile to investigate the accretion flow in the existence of turbulence and cooling but no heating, which sets the upper limit of chaotic cold accretion to some extent. The tests are Models CD10 and CV10 in Table 1. The setup is identical to Models CHD10 and CHV10 except that we do not add any artificial heating.
The snapshots and time evolution are shown in Figures 24 and 25. The mass accretion rate is fluctuating around the spherically symmetric cooling accretion rate. We still find that there is usually a small torus within ∼ 100 pc with the angular momentum of the gas around the inner boundary being Keplerian. The accretion flow is more chaotic compared with the simulation with heating. The disk changes its orientation more frequently. There is little difference between the two types of initial perturbations. Due to the absence of heating, there is more cold gas formed. The accretion rate is therefore higher.