<figure><img src="image_1.png" /></figure>

HAL
open science

# Black hole feeding and feedback: the physics inside the 'sub-grid'

Andrea Negri, Marta Volonteri

► To cite this version:

Andrea Negri, Marta Volonteri. Black hole feeding and feedback: the physics inside the 'sub-grid'. Monthly Notices of the Royal Astronomical Society, 2017, 467 (3), pp.3475-3492. 10.1093/mn-ras/stx362. hal-01554715

HAL Id: hal-01554715

https://hal.science/hal-01554715v1

Submitted on 23 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Black hole feeding and feedback: the physics inside the ‘sub-grid’

A. Negri* and M. Volonteri

Institut d'Astrophysique de Paris, UPMC Univ Paris 6 et CNRS, UMR 7095, Sorbonne Universités, 98 bis bd Arago, F-75014 Paris, France

Accepted 2017 February 8. Received 2017 January 17; in original form 2016 October 14

## ABSTRACT

Black holes (BHs) are believed to be a key ingredient of galaxy formation. However, the galaxy-BH interplay is challenging to study due to the large dynamical range and complex physics involved. As a consequence, hydrodynamical cosmological simulations normally adopt sub-grid models to track the unresolved physical processes, in particular BH accretion; usually the spatial scale where the BH dominates the hydrodynamical processes (the Bondi radius) is unresolved, and an approximate Bondi–Hoyle accretion rate is used to estimate the growth of the BH. By comparing hydrodynamical simulations at different resolutions (300, 30, 3 pc) using a Bondi–Hoyle approximation to sub-parsec runs with non-parametrized accretion, our aim is to probe how well an approximated Bondi accretion is able to capture the BH accretion physics and the subsequent feedback on the galaxy. We analyse an isolated galaxy simulation that includes cooling, star formation, Type Ia and Type II supernovae, BH accretion and active galactic nuclei feedback (radiation pressure, Compton heating/cooling) where mass, momentum and energy are deposited in the interstellar medium through conical winds. We find that on average the approximated Bondi formalism can lead to both over- and underestimations of the BH growth, depending on resolution and on how the variables entering into the Bondi–Hoyle formalism are calculated.

**Key words:** accretion, accretion discs – black hole physics – hydrodynamics – methods: numerical.

## 1 INTRODUCTION

Black holes (BHs) are believed to play a crucial role in galaxy formation and evolution. The correlations between BH mass and host galaxy properties (e.g. Magorrian et al. 1998; Ferrarese & Merritt 2000; Gebhardt et al. 2000; Gültekin et al. 2009; Kormendy & Ho 2013; McConnell & Ma 2013; Kormendy 2016) suggest that their respective evolutions are closely tied. However, two decades after the pioneering work of Magorrian et al. (1998) the impact of BHs formation and growth on the host galaxies evolution is not fully understood yet.

The BH-galaxy coevolution is believed to rely on active galactic nuclei (AGN) feedback, the prodigious energetic emission from the BH accretion disc that, from millipc scales, is able to influence the surrounding medium up to kpc scales (for a review see Fabian 2012). The mutual interaction between BHs and their host galaxies, and in particular the BH accretion processes, has been studied both analytically and numerically. Various models have been proposed to describe BH accretion physics and the relative feedback (Silk & Rees 1998; Fabian 1999; King 2005). Although they have been successful in explaining, for instance, the relationship between BH mass and the galaxy velocity dispersion, $M_{BH}-\sigma$, analytical models

are forced to assume simplified geometries and neglect important physical processes such as mass inflow and star formation.

In the last 10 years hydrodynamical cosmological simulations including BH growth and feedback have become an invaluable tool to shed light on galaxy formation and evolution (e.g. Di Matteo, Springel & Hernquist 2005; Springel, Di Matteo & Hernquist 2005; Robertson et al. 2006; Sijacki et al. 2007; Di Matteo et al. 2008; Booth & Schaye 2009; Schaye et al. 2010; Debuhr, Quataert & Ma 2011; Degraf, Di Matteo & Springel 2011; Teyssier et al. 2011; DeGraf et al. 2012; Dubois et al. 2014; Hirschmann et al. 2014; Vogelsberger et al. 2014b; Khandai et al. 2015; Feng et al. 2016; Volonteri et al. 2016). These simulations are often at the edge of modern supercomputers computational capabilities, and they have to overcome numerous computational challenges. Perhaps the most important is the involved spatial range: a proper simulation of both the large-scale structure of the Universe and the BH accretion scales would require a dynamical range spanning from the Schwarzschild radius of a supermassive BH – $r_s \approx 10^{-6}$ pc for a mass of $10^7 M_{\odot}$ – to a scale of hundreds of Mpc. Bridging all the involved scales is currently computationally infeasible, and it will remain so for the near future. As a consequence, cosmological simulations including AGN feedback have to resort to so-called sub-grid models in the case of BH accretion and feedback.

Perhaps the most commonly used sub-grid accretion models are based on the Bondi–Hoyle–Lyttleton accretion solution (Hoyle &

<footnote>* E-mail: negri@iap.fr</footnote>
Lyttleton 1939; Bondi & Hoyle 1944; Bondi 1952), first proposed by Springel et al. (2005) and then widely adopted in the literature. The main reason of its popularity relies in its simplicity: the Bondi problem describes the spherical stationary inflow of a perfect, non-viscous, non-self-gravitating gas on a BH, subject to polytropic transformations and in absence of any type of feedback. Given these hypotheses, the BH accretion rate can be written as

$$ \dot{M}_{\text{B}} = \frac{4\pi\lambda_{\text{c}}G^2 M_{\text{BH}}^2 \rho_{\infty}}{c_{s_{\infty}}^{3}}, \qquad (1) $$

where

$$ \lambda_{\text{c}} = \frac{1}{4} \left( \frac{2}{5-3\gamma} \right)^{\frac{5-3\gamma}{2(\gamma-1)}}, \qquad (2) $$

G is the gravitational constant, $M_{\text{BH}}$ is the BH mass, $\rho_{\infty}$ and $c_{s_{\infty}}$ are the density and sound speed at great distance (formally infinite) from the BH and $\gamma$ is the polytropic transformation index (Bondi 1952). Thus, with a simple formula BH mass accretion can be quantified.

However, Bondi accretion suffers from various problems if applied to realistic environments. First, the accretion of a turbulent, non-uniform gas flow is not considered (see Barai, Proga & Nagamine 2011, 2012; Park & Ricotti 2012; Gaspari, Ruszkowski & Oh 2013, and reference therein). Secondly, Bondi accretion spherical symmetry implies that gas angular momentum is negligible, which is unlikely to be true in galaxy formation simulations. Gas with angular momentum will circularize before reaching the innermost stable orbit around the BH, creating an accretion disc. Extensive numerical works addressed this issue, by either analytically modelling the angular momentum impact on accretion (Hopkins & Quataert 2011; Rosas-Guevara et al. 2015; Tremmel et al. 2016) or by directly resolving the rotating accretion flow on small scales (Debuhr et al. 2010; Power, Nayakshin & King 2011; Park & Ricotti 2012; Li, Ostriker & Sunyaev 2013; Gaspari, Brighenti & Temi 2015; Curtis & Sijacki 2016; Hopkins et al. 2016; Anglés-Alcázar et al. 2017). Moreover, BHs are not isolated objects as assumed in the Bondi problem, but they usually reside at the bottom of the galactic potential well. The additional pull due to the galactic gravitational field could provide additional material to the BH with respect to the standard Bondi solution, where this additional contribution is neglected (for a fully analytical extended Bondi solution considering the halo gravitational pull see Korol, Ciotti & Pellegrini 2016).

Another level of complexity is then added by further physical mechanisms not considered in the Bondi solution, above all gas cooling, AGN feedback and star formation. Even evaluating the rather simple right-hand-side of equation (1) in a numerical simulation is not trivial: where should $\rho_{\infty}$ and $c_{s_{\infty}}$ be calculated and in which way? In addition, a comparison between Eulerian (grid) and Lagrangian (smooth particle hydrodynamics, hereafter SPH) codes is far from being trivial due to the structural difference in the fluid discretization. Finally, from the numerical point of view, Bondi accretion suffers from resolution dependent issues (see Section 1.1).

Despite the intrinsic difficulty of the Bondi approach in estimating realistic accretion rates in a non-ideal environment, it has been adopted by a large part of the numerical simulation community, providing results that broadly agree with observational constraints on BH growth and AGN evolution. Thus, in the next section we summarize the main sub-grid processes based on Bondi accretion and their pitfalls.

## 1.1 Black hole accretion in hydrodynamical cosmological simulations

All hydrodynamical cosmological simulations employ a sub-grid recipe for BH accretion. The first implementation of BH accretion based on the Bondi solution in a hydrodynamical galaxy-scale simulation has been proposed by Springel et al. (2005). The BH accretion rate was computed as $\dot{M}_{\text{BH}} = \alpha \dot{M}_{\text{B}}$, where $\alpha$ is a dimensionless constant named boost parameter with a typical value of 100,<fnref n="1" /> and $\rho_{\infty}$, $c_{s_{\infty}}$ and $v_{\text{BH}}$ were obtained using a kernel weighted SPH estimator. The accretion rate is then capped at the Eddington limit, thus super-Eddington accretion is not allowed by construction.<fnref n="2" /> This is the accretion rate employed by the Magneticum simulation (Hirschmann et al. 2014).

Booth & Schaye (2009) motivate the adoption of the boosted accretion by noting that at the resolutions of cosmological simulations, the Bondi radius (the scale where the BH dominates the hydrodynamics) and the properties of the interstellar medium (ISM) cold phase are not resolved. Thus, they argue that the temperature of the material accreted on the BH would be overestimated, and its density underestimated, leading to a very low accretion rate in equation (1) if large values of $\alpha$ are not adopted. However, they also discuss the limits of these justifications: when the BH is in a state of a hot accretion, which typically shows small deviations from spherical symmetry, the Bondi assumptions are fulfilled, and a large boost factor is no longer needed. Moreover, the combination of the Eddington limiter and the constant $\alpha$ factor can present unexpected consequences. They show that by assuming $\alpha = 100$, BHs with a mass > $10^6$ M$_{\odot}$ reach a sub-Eddington accretion only when the gas density has become extremely low, almost two orders of magnitude less than the star formation threshold. In other words, BH accretion a constant boosting will be sensitive to the ISM density (and thus more easily regulated by AGN feedback) only when the BH surroundings have been almost completely emptied of gas. To minimize these effects, they proposed a density-dependant boost factor, where $\alpha = 1$ if $n_{\text{gas}}/n_{\text{sf}} < 1$ and $\alpha = (n_{\text{gas}}/n_{\text{sf}})^\beta$ otherwise, where $n_{\text{gas}}$ is the gas number density, $n_{\text{sf}} = 0.1 \text{ cm}^{-3}$ is the gas density threshold for star formation and $\beta = 2$; the accretion is then limited at the Eddington luminosity, assuming a radiative efficiency of 10 per cent.

Following Springel et al. (2005) and Booth & Schaye (2009), the boosted Bondi accretion has been adopted (and modified) by a large number of authors. Pelupessy, Di Matteo & Ciardi (2007) do not follow a constant $\alpha$ prescription; instead $\dot{M}_{\text{B}}$ is the sum of both the hot and cold phase Bondi accretion rates, calculated with boosting factors equal to their respective mass fractions. This scheme is adopted by the Massiveblack (DeGraf et al. 2012), Massiveblack II (Khandai et al. 2015) and BLUETIDES simulations (Feng et al. 2016; Di Matteo et al. 2016). A similar approach is employed by Steinborn et al. (2015), where they calculate the physical quantities of equation (1) using 295 neighbours, and calculating separately the accretion rate for the cold and the hot gas, with $\alpha = 100$ and 10,

<footnote n="1">¹ The value of $\alpha$ is not explicitly mentioned in the earlier works of hydrodynamical cosmological simulations, but see table 2 of Booth & Schaye (2009) for a list of the adopted boost parameters.</footnote>

<footnote n="2">² We do not further discuss here the validity of capping the mass accretion rate, see Volonteri, Silk & Dubus (2015) and Inayoshi, Haiman & Ostriker (2016) for an overview. Note that Lusso & Ciotti (2011) describe the transition from the pure hydrodynamical to the optically thin radiation-dominated regime by self-consistently solving the Bondi problem including electron scattering.</footnote>
respectively (this is the default BH sink algorithm of the GADGET-3-x code). Finally, Choi et al. (2012) apply the average directly on $\dot{M}_{BH}$ with $\alpha = 100$ calculated over 64 neighbours, in order to reduce the dependence on the numerical kernel. The actual accretion of the particles around the BH is performed in a stochastic fashion. Only the particles inside the Bondi radius are accreted, with a probability proportional to the amount of overlap between the particle and the Bondi radius (if the Bondi radius is unresolved, it is set equal to the particles softening length). An additional probability taking into account the free fall time for each particles is considered. The same prescription is used in Choi et al. (2015) without $\alpha$ boosting.

Some of these prescriptions have been adopted as the default accretion algorithm in numerical codes. The GADGET-3-OWLS code implements the Booth & Schaye (2009) accretion. This modified version of GADGET-3 is employed in the OWLS (Schaye et al. 2010), COSMO-OWLS (Le Brun et al. 2014) and BAryons and HALoes of MAssive Systems (BAHAMAS) simulations (McCarthy et al. 2017). The prescription has been modified in the Evolution and Assembly of GaLaxies and their Environments (EAGLE) simulation (Schaye et al. 2015) to take into account the effects of gas angular momentum on BH accretion (adopted from Rosas-Guevara et al. 2015); at variance with the majority of cosmological simulations, the boost factor is no longer employed.

The implementation in CHANGA modifies the Bondi formalism used in GASOLINE (Bellovary et al. 2010) to include the suppressive effect of gas rotational support (Tremmel et al. 2016). In both cases, the accretion rate is computed as the sum of the Bondi accretion rate of each of the 32 individual gas particles nearest the BH, rather than simply averaging the gas quantities to insert them in equation (1). A boost factor may or may not be included.

BHs in the AREPO code (Springel 2010) are collision-less sink particles whose accretion rate is calculated following equation (1), where a volume weighted mean is performed (Vogelsberger et al. 2013; Sijacki et al. 2015). However, Vogelsberger et al. (2013) note that when no star-forming gas is present in the immediate vicinity of the BH, the boosted Bondi prescription can lead to an overestimate of the mass accretion rate. Thus, in this regime they lower the accretion rate based on a pressure criterion. This configuration has been employed in the Illustris simulation (Vogelsberger et al. 2014b), while Curtis & Sijacki (2015) proposed a new super-Lagrangian refinement to reach the Bondi radius near a BH.

Moving to adaptive mesh refinement, in RAMSES (Teyssier 2002) BHs accrete following equation (1) where the right-hand-side is composed by mean values calculated by disseminating clouds particles around the BH then weighted using a kernel (Dubois et al. 2012). The accretion rate is then boosted following the Booth & Schaye (2009) prescription. In the Horizon-AGN simulation (Dubois et al. 2014), the mean quantities are computed in a mass weighted fashion (Dubois, private communication). New sink schemes are proposed in Bleuler & Teyssier (2014). All the above mentioned sub-grid models are Eddington limited. In the ENZO code (Kim et al. 2011; Bryan et al. 2014), Bondi accretion is assumed, and the density in equation (1) is extrapolated at the Bondi radius, no boost factor is taken into account.

Despite the details of the BH accretion scheme not always being stated explicitly in the literature, two clear approaches emerge: volume weighting (the SPH kernel weighted estimators are similar to volume weighting) and mass weighting for the calculation of the Bondi rate. Finally, there have been recent proposals that forgo the Bondi approximation in favour of angular-momentum driven sub-grid implementations (e.g. see Debuhr et al. 2010; Anglés-Alcázar et al. 2017).

All the sub-grid algorithms presented above rely on the common assumption that at the resolutions of cosmological simulations the ISM in the surroundings of a BH is smoothed, so its density and temperature are respectively under and overestimated, leading to an accretion rate lower than in reality. However, the boost factor $\alpha$ is poorly constrained, being problem and resolution dependent. In addition, comparative studies such as Wurster & Thacker (2013) and Elahi et al. (2016) employ BH sub-grid recipes in their simulations, while a comparison with a non-parametrized accretion would be ideal. High-resolution studies, where the Bondi radius is resolved and the accretion rate is computed in a non-parametrized way, only simulate the innermost regions of a galaxy (up to the first few hundred of pc, see Kurosawa & Proga 2009a,b; Kurosawa, Proga & Nagamine 2009; Barai et al. 2011; Hopkins & Quataert 2011; Barai et al. 2012; Hopkins et al. 2016) or multiple simulations on different scales of a single object (Hopkins & Quataert 2010). Multi-scale studies in presence of feedback has been performed by Ciotti, Ostriker and collaborators (Ciotti & Ostriker 1997, 2001, 2007; Ostriker et al. 2010; Novak, Ostriker & Ciotti 2011; Ciotti et al. 2017, and references therein), but no comparison with Bondi accretion has been performed. One should not be surprised in finding that at a particular time the $\dot{M}_B$ under- or overestimate the real accretion rate, e.g. when large amounts of gas very far from the BH produces a large $\dot{M}_B$, while the BH can actually be in a state of low accretion. However, due to the non-linear nature of hydrodynamics and in presence of additional physical processes such as radiative cooling, AGN feedback and star formation it is extremely difficult to quantify a priori both the BH growth with time and the global effect of adopting the Bondi solution as a proxy for BH accretion.

Our aim is to probe two types of Bondi algorithms (using mass and volume weighted quantities) and their global effects on the BH mass by comparing 2.5D hydrodynamical simulations of an isolated galaxy at different resolutions to sub-parsec runs with non-parametrized accretion (calculated as the mass flux on the BH, see Section 2.3). In our reference run we simulate radii from 0.1 pc to 250 kpc, thus resolving the Bondi radius for the adopted BH mass at all the possible gas temperatures (Section 2.2), in presence of cooling, star formation, Type Ia and Type II supernovae and AGN feedback (radiation pressure, Compton heating/cooling, broad absorption line winds).

The paper is organized as follows: in Section 2, we describe the code input physics and our initial conditions. In Section 3, we present the results of our simulations, while in Section 4, we discuss our results and we compare them to other recent works. Finally, in Section 5, we draw our conclusions.

## 2 THE SIMULATIONS

We use a customized, parallel version of the hydrodynamical code ZEUS (Stone & Norman 1992; Hayes et al. 2006). For the most part the detailed input physics is described in Novak et al. (2011); Novak, Ostriker & Ciotti 2012 hereafter N11 and N12, here we briefly summarize the main implementations and differences with respect to N11 and N12.

### 2.1 The input physics

The mass, momentum and energy injections/sinks from stellar evolution, Type Ia and II supernovae, star formation and radiative cooling (bremsstrahlung and lines cooling, with an exponential cut-off below $10^4$ K) are the same as N11. In particular, the star
formation recipe can reproduce the Kennicutt–Schmidt relation (Kennicutt 1998) with an adopted efficiency equal to 0.1.

Radiative feedback from the BH is taken into account by considering photo-ionization heating, Compton heating/cooling in the energy equation and radiation pressure due to electron scattering and photo-ionization in the momentum equation. Given a BH accretion rate $\dot{M}_{BH}$, the total BH luminosity is

$$L_{BH} = \epsilon_{EM} \dot{M}_{BHC} c^2, \quad (3)$$

where c is the speed of light and $\epsilon_{EM}$ is the electromagnetic efficiency. We here adopt a dependence of $\epsilon_{EM}$ on the accretion rate that captures the transition from optically thick and geometrically thin, radiatively efficient accretion discs (Shakura & Sunyaev 1973), to optically thin, geometrically thick, radiatively inefficient advection dominated accretion flows (Narayan & Yi 1994):

$$\epsilon_{EM} = \frac{\epsilon_0 A \dot{m}}{1 + A \dot{m}}, \quad \dot{m} \equiv \frac{\dot{M}_{BH}}{\dot{M}_{Edd}} = \frac{\epsilon_0 \dot{M}_{BHC} c^2}{L_{Edd}}, \quad (4)$$

here $\dot{m}$, $\dot{M}_{Edd}$ and $L_{Edd}$ are the dimensionless BH and Eddington mass accretion rates, respectively, with $A = 100$ and $\epsilon_0 = 0.125$, while $L_{Edd} = 4\pi G M_{BH} m_p c / \sigma_T$ is the Eddington luminosity, where $m_p$ and $\sigma_T$ are the proton mass and the Thomson cross-section.

Mechanical feedback is included as a broad absorption line (BAL) wind originating from the BH accretion disc. In N11, the BAL wind is modelled as a kinetic wind following Ostriker et al. (2010), which guarantees that the mass, energy and momentum carried by the wind are self-consistent. In this approach, the mass outflow (towards the BH) at the inner radial grid is transferred into a sub-grid accretion disc model, then a fraction of the inflow is injected back into the computational grid as a BAL wind, while the rest is accreted on the BH. This model can be formalized as:

$$\dot{M}_{BH} = \frac{\dot{M}_{in}}{1 + \eta}, \quad (5)$$

$$\dot{M}_{out} = \eta \dot{M}_{BH}, \quad (6)$$

$$L_w = \epsilon_w \dot{M}_{BHC} c^2, \quad (7)$$

$$\dot{p}_w = \eta \dot{M}_{BH} v_w, \quad (8)$$

$$\eta \equiv 2\epsilon_w c^2 / v_w^2, \quad (9)$$

$$\epsilon_w = \frac{\epsilon_{w0} A_w \dot{m}}{1 + A_w \dot{m}}, \quad (10)$$

where $\dot{M}_{in}$ is the mass inflow rate at the innermost radial boundary, $M_{out}$, $\epsilon_w$, $L_w$, $\dot{p}_w$ and $v_w$ are the BAL wind injection mass rate, efficiency, mechanical luminosity, injection momentum rate, injection velocity, respectively; in this work we consider $v_w = 10^4$ km s$^{-1}$, $\epsilon_{w0} = 10^{-4}$ and $A_w = 1000$.

We model the BAL wind as a direct injection, at the first radial grid, of mass and momentum within two symmetric cones above and below the equatorial plane, with a fixed half-opening angle of 45° containing half of the energy emitted, i.e. this means that the wind is visible from ~1/5 of the entire solid angle. Equations (6)–(8) provide the total injection rates, which are distributed on the innermost radial grid following the prescription of Ciotti et al. (2017) (see their equations 25–28). By defining the wind thermalized when its velocity drops below five times the velocity dispersion of the halo (150 km s$^{-1}$), in the region outside the double cone with semi-aperture of 45 deg, the thermalization radius is typically ~10 pc,

with a maximum at 40 pc, while in the polar region is ~100 pc, with maximum of 300 pc.

We stress that AGN feedback is a combination of radiative (via photo-ionization heating and Compton heating/cooling) and kinetic forms, instead of the thermal feedback (injection of thermal energy with an efficiency parameter setting the radiation–gas coupling) usually employed in cosmological simulations (see also Barai et al. 2014). At variance with N12 in this work we do not consider the treatment of dust formation and destruction, and the associated radiation pressure effects (see also Hensley, Ostriker & Ciotti 2014), neither radiation pressure due to stellar light in different bands (see e.g. Ciotti & Ostriker 2012).

## 2.2 Numerical grid

The code employs a 2D axisymmetric, computational grid in spherical coordinates $(r, \theta)$. The mesh is not dynamically refined, however the non-uniform spacing allows to concentrate more gridpoints near the centre of the simulations box. 2D axisymmetric simulations offer the best tradeoff between realism and computational requirements. Neglecting the dependence from the azimuthal angle drastically decreases the computational effort; at the same time, Rayleigh–Taylor instabilities can fragment a cold shell, thus forming a clumpy, turbulent multiphase medium, and the BH can simultaneously accrete and eject mass in a non-spherical shape, which is not allowed in 1D runs. However, in absence of non-axisymmetric torques an infalling rotating gas blob will circularize without feeding the BH. For this very reason, we neglect a large level of rotation, in order to accrete on the central object.

In this work, we adopt different radial spacing for different sets of simulations, whereas the $\theta$ setup is maintained fixed. For our reference models (see Sections 2.3 and 3.1), the mesh consists of 512 × 30 gridpoints, spanning from 0.1 pc to 250 kpc in the radial direction (with the spacing growing following a geometric progression), with a central resolution of 0.1 pc. The grid in the $\theta$ direction is regularly sampled and to avoid the coordinate singularity on the axis of symmetry we exclude a narrow double-cone of semi-aperture $\pi \times 10^{-4}$ rad. In this configuration, we obtain an almost constant 0.1 pc resolution in the first 10 pc from the centre, and then resolution grows linearly with radius. As mentioned above, some simulation sets adopt a different, low-resolution grid (3, 30, 300 pc), while the innermost and outermost radial grid are maintained at 0.1 pc and 250 kpc, respectively.

In hydrodynamical cosmological simulations BH accretion is commonly treated as sub-grid physics. Due to the large dynamical range required to resolve properly galaxy evolution in a cosmological context, one has to account for extremely different spatial scales. A crucial scale for the galaxy-BH co-evolution is the Bondi radius, formally defined as $r_B \equiv GM_{BH}/c_s^2$, where the BH dominates the hydrodynamical evolution of the ISM. The temperature dependence is apparent: given $M_{BH}$, the higher the gas temperature, the smaller the Bondi radius. Thus, a numerical simulation that aims at following the ISM during accretion events should properly resolve the Bondi radius for all the temperatures of the inflowing gas. In our case, with the high-resolution grid we resolve spatial scales corresponding to temperatures less than $10^8$ K. Note that while the highest temperature allowed in the code derives from the complete thermalization of the BAL wind (~2.5 × $10^9$ K), this hot gas will be outflowing by construction. As a consequence, the temperature upper limit for the inflowing gas is the Compton temperature $T_c = 2 \times 10^7$ K, so we actually resolve the Bondi radius for all the (inflowing gas) temperatures.
For all the simulations in this work, reflecting boundary conditions were set along the $\theta$ boundary (with the gas rotational velocity flipped in sign); at the outer edge of the spherical simulated box the fluid is free to flow out of the computational grid. The inner boundary conditions are a special case: they are set on reflection when we adopt Bondi accretion (see Section 2.3), while in the case of non-parametrized accretion we adopt outflow conditions (off of the computational grid towards the centre of the simulation) reported in N11. As already noted in Prasad, Sharma & Babul (2015), the cold gas can stick to the $\theta$ boundary when reflecting boundary conditions are applied, thus inducing an unphysical polar accretion flow. To avoid this phenomenon, they exclude the first 8 grid zones near each pole in the calculation of the BH accretion rate. However, this can still lead to an enhanced star formation confined in the polar direction due to the cold gas presence therefore we tuned the initial ISM angular momentum distributions near the poles to minimize cold ISM accumulation along the $\theta$ boundary (see below).

## 2.3 BH accretion

We employ different types of accretion models. For our high resolution, reference simulations, we employ a non-parametrized accretion scheme, 'flux' hereafter: the mass accretion rate $\dot{M}_{in}$ is simply the flux out of the grid and towards the BH calculated on the innermost radial grid, which is at a distance of 0.1 pc from the BH. This allows the code to treat the ISM hydrodynamics self-consistently, since the gas is accreted on the BH if and only if is actually crossing our innermost cell, well below the Bondi radius.

To quantitatively study the effect of the most common sub-grid accretion prescriptions adopted in hydrodynamical cosmological simulations, we implemented four different types of Bondi accretion schemes. For each scheme we define an accretion radius $r_{acc}$, that includes the volume where we calculate $\rho_{\infty}$ and $c_{s\infty}$ in equation (1) as a volume or mass weighted average. Two accretion schemes employ $\dot{M}_{in} = \dot{M}_{B}$, and other two cap accretion at the accretion rate corresponding to the Eddington luminosity assuming a radiative efficiency of 10 per cent, i.e. $\dot{M}_{in} = \min(\dot{M}_{B}, \dot{M}_{Edd})$. Thus, the four schemes differ in the way $\rho_{\infty}$ and $c_{s\infty}$ are calculated and whether an Eddington limit is applied. In order to ensure mass conservation, when the BH accretes a certain mass of gas $M$ in a timestep, we drain the same amount of gas from the volume inside $r_{acc}$. For each cell the rate of mass depletion $\dot{M}$ is given by: $\dot{M} = \dot{M}_{in}M_{cell}/M_{tot}$, where $M_{cell}$ and $M_{tot}$ are the gas mass inside a cell and the gas mass enclosed by $r_{acc}$, respectively. For sake of brevity from now on we will assume that a Bondi simulation employs mass weighting, when not explicitly stated.

## 2.4 The initial conditions

The total gravitational potential is an isothermal sphere with a velocity dispersion $\sigma = 150$ km s$^{-1}$, accounting for the stellar and dark matter distributions, plus a point mass potential for the central BH. The gas is not self-gravitating.

The stellar distribution is described by a Jaffe (1983) profile with a half-light radius of 5.3 kpc, a total mass $M_{*} = 3 \times 10^{10}$ M$_{\odot}$, and the mass-to-light ratio is assumed to be 5.8 in solar unit in the B band, independently of the position (the structural and dynamical properties of these galaxy models are discussed in detail in Ciotti, Morganti & de Zeeuw 2009). The initial BH mass is $M_{BH} = 10^{-3} M_{*}$, thus placing the galaxy near the Magorrian relation (Magorrian et al. 1998; Kormendy & Ho 2013).

The initial ISM density is a spherical distribution described by

$$ \rho(r) = \frac{\mu m_{p} n_{0}}{(1 + \tilde{r}^{2})^{\beta}}, \quad \tilde{r} \equiv r/r_{c}, \qquad (11) $$

where $m_p$, $\mu$, $r_c$, $n_0$, $\beta$ are the proton mass, mean molecular weight, core radius, central number density and slope of the ISM distribution, respectively. We adopt $\mu = 0.62$ (solar metallicity), for simplicity $r_c = 5.3$ kpc (consistent with the stellar distribution) and $\beta = 1.3$ to obtain a density similar to the critical baryon density in the last radial grid when $n_0 = 0.01$ cm$^{-3}$ is adopted (see Section 3). Then, the temperature profile is calculated by solving the equation of hydrostatic equilibrium for the ISM with a zero pressure boundary condition at the outer edge of the grid.

In principle, in 2D axisymmetric simulations, gas and stars can have an arbitrary rotational velocity $u_{\phi}$ and $v_{\phi}$ as a function $(r, \theta)$. Negri, Ciotti & Pellegrini (2014a) and Negri et al. (2014b) show that a rotating stellar population is able to drag the ISM, thus producing a gaseous halo almost corotating with the stars. However, angular momentum conservation in rotating galaxies would lead to the formation of a cold, rotationally supported disc of kpc-scale (Negri et al. 2014a,b; Negri, Pellegrini & Ciotti 2015), thus preventing any substantial accretion on the central BH. In real objects, these massive discs would be prone to gravitational instabilities that can transport angular momentum outwards and so drive gas on the central object (Bertin & Lodato 2001; Hopkins & Quataert 2011). Yet, this kind of instability cannot be captured in 2D axisymmetric simulations without a prescription for angular momentum transport. In order to allow accretion on the BH, in this work we limit ourself in the low rotation regime, by constraining the centrifugal barrier radius to reside inside the innermost radial grid. The local specific angular momentum of the circular orbit in the meridional plane is given by

$$ j(r, \theta) = \sin^2 \theta \sqrt{r^2 \sigma^2 + G M_{BH,0} r}. \qquad (12) $$

Thus, we impose the following rotational field $u_{\phi}$ ($v_{\phi}$) for the gaseous (stellar) component:

$$ \frac{1}{u_{\phi}} = \frac{d}{\sigma R} + \frac{1}{f \sigma} + \frac{R}{j(r_{\min}, \pi/2)}, \qquad (13) $$

where $R = r \sin \theta$ is the distance from the axis, $M_{BH,0}$ is the BH mass at the beginning of the simulation, $d = 10$ pc, $f = 1$, $j(r_{\min}, \pi/2)$ is the specific angular momentum of the circular orbit at radius $r_{\min}$ in the equatorial plane and $r_{\min}$ is the radius of the first grid. The resulting profile gives solid body rotation at $R < d$, and a constant angular momentum distribution at large radii, while the $f$ term guarantees that the velocity field never exceeds $\sigma$. As mentioned in Section 2.2, if a gas blob hits the $\theta$ boundary, it remains 'attached' to it, thus leading to un unphysical BH accretion and star formation density. In order to avoid this numerical issue, we modified the stellar rotation field $v_{\phi}$ near the $\theta$ boundary, by taking it as the maximum between equation (13) and $5\psi(r, \theta)j(r, \theta)/R$, where

$$ \psi(r, \theta) = \frac{1}{8} \left[ 2 + \tanh \frac{\theta_1 - \theta}{\sigma_{\theta}} + \tanh \frac{\theta_1 + \theta - \pi}{\sigma_{\theta}} \right] \qquad (14) $$

$$ \times \left[ 1 + \tanh \frac{r_2 - r}{\sigma_r} \right] \left[ 1 + \tanh \frac{r - r_1}{\sigma_r} \right], \qquad (15) $$

where $\theta_1 = 0.15$ rad, $r_1 = 70$ pc and $r_2 = 300$ kpc. This formulation provides a smooth transition from the biconical region described by equation (15) and the rest of simulation box, where the rotation law is provided by equation (13).
<table><caption>**Table 1.** Mass weighted simulations.</caption><thead><tr><th>Accretion<br>(1)</th><th>Resolution<br>(pc)<br>(2)</th><th>r<sub>acc</sub><br>(pc)<br>(3)</th><th>n<sub>0</sub><br>(cm<sup>-3</sup>)<br>(4)</th><th>ΔM<sub>BH</sub><br>(10<sup>6</sup> M<sub>☉</sub>)<br>(5)</th><th>ΔM<sub>*</sub><br>(10<sup>6</sup> M<sub>☉</sub>)<br>(6)</th><th>⟨SFR⟩<br>(M<sub>☉</sub> yr<sup>-1</sup>)<br>(7)</th><th>M<sub>gas</sub><br>(10<sup>10</sup> M<sub>☉</sub>)<br>(8)</th></tr></thead><tbody><tr><td>B300_S_H<sup>a</sup></td><td>0.1</td><td>300</td><td>1</td><td>2772.13</td><td>418.17</td><td>10.06</td><td>36.10</td></tr><tr><td>B30_S_H</td><td>0.1</td><td>30</td><td>1</td><td>15.52</td><td>2092.68</td><td>41.85</td><td>35.99</td></tr><tr><td>B3_S_H</td><td>0.1</td><td>3</td><td>1</td><td>2.27</td><td>2841.47</td><td>56.83</td><td>36.84</td></tr><tr><td>B300_S_L</td><td>0.1</td><td>300</td><td>0.01</td><td>12.86</td><td>1.05</td><td>0.02</td><td>5.39</td></tr><tr><td>B30_S_L</td><td>0.1</td><td>30</td><td>0.01</td><td>8.86</td><td>4.68</td><td>0.10</td><td>5.39</td></tr><tr><td>B3_S_L</td><td>0.1</td><td>3</td><td>0.01</td><td>1.30</td><td>16.71</td><td>0.33</td><td>0.37</td></tr><tr><td>B300_C_H</td><td>300</td><td>300</td><td>1</td><td>171.89</td><td>186.10</td><td>3.72</td><td>30.96</td></tr><tr><td>B30_C_H</td><td>30</td><td>30</td><td>1</td><td>100.22</td><td>1085.32</td><td>21.71</td><td>31.18</td></tr><tr><td>B3_C_H</td><td>3</td><td>3</td><td>1</td><td>1.24</td><td>3059.29</td><td>61.19</td><td>31.98</td></tr><tr><td>B300_C_L</td><td>300</td><td>300</td><td>0.01</td><td>17.87</td><td>0.44</td><td>0.01</td><td>0.38</td></tr><tr><td>B30_C_L</td><td>30</td><td>30</td><td>0.01</td><td>2.62</td><td>0.30</td><td>0.01</td><td>0.38</td></tr><tr><td>B3_C_L</td><td>3</td><td>3</td><td>0.01</td><td>0.56</td><td>12.22</td><td>0.24</td><td>0.37</td></tr><tr><td>BE300_S_H</td><td>0.1</td><td>300</td><td>1</td><td>32.01</td><td>1157.94</td><td>23.16</td><td>35.98</td></tr><tr><td>BE30_S_H</td><td>0.1</td><td>30</td><td>1</td><td>7.47</td><td>2240.76</td><td>44.82</td><td>36.76</td></tr><tr><td>BE3_S_H</td><td>0.1</td><td>3</td><td>1</td><td>1.44</td><td>2121.02</td><td>42.42</td><td>36.63</td></tr><tr><td>BE300_S_L</td><td>0.1</td><td>300</td><td>0.01</td><td>31.63</td><td>1.17</td><td>0.02</td><td>5.39</td></tr><tr><td>BE30_S_L</td><td>0.1</td><td>30</td><td>0.01</td><td>5.39</td><td>4.96</td><td>0.10</td><td>0.37</td></tr><tr><td>BE3_S_L</td><td>0.1</td><td>3</td><td>0.01</td><td>1.15</td><td>15.41</td><td>0.31</td><td>5.38</td></tr><tr><td>BE300_C_H</td><td>300</td><td>300</td><td>1</td><td>30.67</td><td>1787.89</td><td>35.76</td><td>30.42</td></tr><tr><td>BE30_C_H</td><td>30</td><td>30</td><td>1</td><td>32.25</td><td>2180.58</td><td>43.61</td><td>30.70</td></tr><tr><td>BE3_C_H</td><td>3</td><td>3</td><td>1</td><td>0.76</td><td>1103.63</td><td>22.07</td><td>30.00</td></tr><tr><td>BE300_C_L</td><td>300</td><td>300</td><td>0.01</td><td>17.05</td><td>1.09</td><td>0.02</td><td>0.38</td></tr><tr><td>BE30_C_L</td><td>30</td><td>30</td><td>0.01</td><td>1.87</td><td>0.35</td><td>0.01</td><td>0.38</td></tr><tr><td>BE3_C_L</td><td>3</td><td>3</td><td>0.01</td><td>0.56</td><td>24.04</td><td>0.48</td><td>0.37</td></tr></tbody></table>

Note. (1) Simulation name, containing the accretion scheme (F, B, BE and BV for flux, Bondi, Bondi capped at Eddington rate and Bondi volume weighted, respectively), accretion radius (3, 30 and 300 pc), grid type (S, C for standard high resolution and coarse grid, respectively) and central initial gas number density (H or L for $n_0 = 1, 0.01 \text{ cm}^{-3}$, respectively). (2) Grid resolution. (3) Accretion radius for B and BE simulations. (4) Initial gas central number density. (5) Gas mass accreted on the BH. (6) Stellar mass formed during the run. (7) Mean SFR calculated over the entire run. (8) Final gas mass.

<table><caption>**Table 2.** Flux simulations.</caption><thead><tr><th>Accretion<br>(1)</th><th>Resolution<br>(pc)<br>(2)</th><th>r<sub>acc</sub><br>(pc)<br>(3)</th><th>n<sub>0</sub><br>(cm<sup>-3</sup>)<br>(4)</th><th>ΔM<sub>BH</sub><br>(10<sup>6</sup> M<sub>☉</sub>)<br>(5)</th><th>ΔM<sub>*</sub><br>(10<sup>6</sup> M<sub>☉</sub>)<br>(6)</th><th>⟨SFR⟩<br>(M<sub>☉</sub> yr<sup>-1</sup>)<br>(7)</th><th>M<sub>gas</sub><br>(10<sup>10</sup> M<sub>☉</sub>)<br>(8)</th></tr></thead><tbody><tr><td>F_S_H</td><td>0.1</td><td>—</td><td>1</td><td>0.53</td><td>2523.00</td><td>50.46</td><td>36.66</td></tr><tr><td>F_S_L</td><td>0.1</td><td>—</td><td>0.01</td><td>0.49</td><td>14.60</td><td>0.29</td><td>0.36</td></tr></tbody></table>

Note. Same as Table 1.

## 3 RESULTS

Here, we present the main results of the simulations, summarized in Tables 1–4. As reference runs, in Section 3.1 we consider two different runs employing the flux accretion scheme and 0.1 pc grid resolution, but having different initial gas density ($n_0 = 1; 0.01 \text{ cm}^{-3}$). For each flux simulation we produced two different families with a sub-grid Bondi scheme: (1) a pure Bondi accretion, as described in Section 2.3 and (2) a Bondi accretion capped at the Eddington limit. In each of these families we keep the 0.1 pc hydrodynamical grid, but consider three different values of the accretion radius ($r_{acc} = 3, 30, 300$ pc). We retain the 0.1 pc resolution hydrodynamical grid of the flux simulations in order to isolate the effects of the Bondi accretion scheme (Section 3.2) from grid effects. In addition, the mechanical feedback is deposited in the innermost radial grids as in the flux runs, so the ratio of mass/momentum/energy injection and the injection radius is independent of $r_{acc}$. This setup

is similar to the BH accretion implementation of Curtis & Sijacki (2015); we discuss few tailored simulations employing their accretion scheme (see Section 4). Bondi high-resolution simulations with density and sound speed calculated as volume weighted quantities are presented in Section 3.4. Finally, we employ two additional families of runs, with identical initial conditions and BH accretion schemes, but with a central spatial resolution equal to the accretion radius $r_{acc}$ explored in the above families (Section 3.3).

Each simulation is run for 100 Myr; the evolution of the main quantities is rescaled in time in the figures, so that the first burst of BH accretion (at $10^{-2}$ the BH Eddington ratio) is at 2 Myr, ending with a presented evolution of 50 Myr. This time-scale is longer than the Salpeter time ($\simeq 45$ Myr) and the crossing time inside the effective radius ($\simeq 35$ Myr, upper limit in absence of AGN feedback). The simulations do not show a significant change in the trends after this time. We avoid running the isolated galaxies for a much longer time-scale because isolated galaxies lack the cosmic

<footnote>a Run terminated at $t = 36$ Myr.</footnote>
<table><caption>**Table 3.** Simulations without AGN feedback.</caption><thead><tr><th>Accretion<br>–<br>(1)</th><th>Resolution<br>(pc)<br>(2)</th><th>r<sub>acc</sub><br>(pc)<br>(3)</th><th>n<sub>0</sub><br>(cm<sup>-3</sup>)<br>(4)</th><th>ΔM<sub>BH</sub><br>(10<sup>6</sup> M<sub>☉</sub>)<br>(5)</th><th>ΔM<sub>*</sub><br>(10<sup>6</sup> M<sub>☉</sub>)<br>(6)</th><th>⟨SFR⟩<br>(M<sub>☉</sub> yr<sup>-1</sup>)<br>(7)</th><th>M<sub>gas</sub><br>(10<sup>10</sup> M<sub>☉</sub>)<br>(8)</th></tr></thead><tbody><tr><td colspan="8"><strong>Mass weighted</strong></td></tr><tr><td>BE300_C_H</td><td>300</td><td>300</td><td>1</td><td>60.20</td><td>2088.61</td><td>41.77</td><td>29.35</td></tr><tr><td>BE30_C_H</td><td>30</td><td>30</td><td>1</td><td>56.96</td><td>5702.33</td><td>114.04</td><td>29.26</td></tr><tr><td>BE3_C_H</td><td>3</td><td>3</td><td>1</td><td>56.65</td><td>7330.36</td><td>146.54</td><td>29.79</td></tr><tr><td>BE300_C_L</td><td>300</td><td>300</td><td>0.01</td><td>13.87</td><td>4.11</td><td>0.08</td><td>0.36</td></tr><tr><td>BE30_C_L</td><td>30</td><td>30</td><td>0.01</td><td>43.48</td><td>43.96</td><td>0.87</td><td>0.35</td></tr><tr><td>BE3_C_L</td><td>3</td><td>3</td><td>0.01</td><td>44.34</td><td>44.25</td><td>0.88</td><td>0.35</td></tr><tr><td colspan="8"><strong>Volume weighted</strong></td></tr><tr><td>BEV300_C_L</td><td>300</td><td>300</td><td>0.01</td><td>7.22</td><td>73.47</td><td>1.47</td><td>0.30</td></tr><tr><td>BEV30_C_L</td><td>30</td><td>30</td><td>0.01</td><td>49.43</td><td>179.47</td><td>3.79</td><td>0.35</td></tr><tr><td>BEV3_C_L</td><td>3</td><td>3</td><td>0.01</td><td>234.47</td><td>2.40</td><td>0.05</td><td>0.35</td></tr><tr><td>BEV300_C_H</td><td>300</td><td>300</td><td>1</td><td>60.33</td><td>2268.57</td><td>45.37</td><td>29.33</td></tr><tr><td>BEV30_C_H</td><td>30</td><td>30</td><td>1</td><td>57.85</td><td>21600.51</td><td>432.01</td><td>29.18</td></tr><tr><td>BEV3_C_H</td><td>3</td><td>3</td><td>1</td><td>61.06</td><td>5731.05</td><td>114.62</td><td>29.04</td></tr></tbody></table>

*Note.* Same as Table 1.

<table><caption>**Table 4.** Volume weighted simulations.</caption><thead><tr><th>Accretion<br>–<br>(1)</th><th>Resolution<br>(pc)<br>(2)</th><th>r<sub>acc</sub><br>(pc)<br>(3)</th><th>n<sub>0</sub><br>(cm<sup>-3</sup>)<br>(4)</th><th>ΔM<sub>BH</sub><br>(10<sup>6</sup> M<sub>☉</sub>)<br>(5)</th><th>ΔM<sub>*</sub><br>(10<sup>6</sup> M<sub>☉</sub>)<br>(6)</th><th>⟨SFR⟩<br>(M<sub>☉</sub> yr<sup>-1</sup>)<br>(7)</th><th>M<sub>gas</sub><br>(10<sup>10</sup> M<sub>☉</sub>)<br>(8)</th></tr></thead><tbody><tr><td>BV300_S_H</td><td>0.1</td><td>300</td><td>1</td><td>0.17</td><td>2584.31</td><td>52.49</td><td>36.20</td></tr><tr><td>BV30_S_H</td><td>0.1</td><td>30</td><td>1</td><td>0.25</td><td>2715.80</td><td>55.10</td><td>36.54</td></tr><tr><td>BV3_S_H</td><td>0.1</td><td>3</td><td>1</td><td>0.58</td><td>1405.10</td><td>28.52</td><td>35.19</td></tr><tr><td>BV300_S_L</td><td>0.1</td><td>300</td><td>0.01</td><td>0.02</td><td>29.81</td><td>0.90</td><td>5.38</td></tr><tr><td>BV30_S_L</td><td>0.1</td><td>30</td><td>0.01</td><td>0.26</td><td>28.75</td><td>0.63</td><td>5.39</td></tr><tr><td>BV3_S_L</td><td>0.1</td><td>3</td><td>0.01</td><td>0.47</td><td>27.95</td><td>0.56</td><td>5.38</td></tr><tr><td>BV300_C_H</td><td>300</td><td>300</td><td>1</td><td>1.37</td><td>7489.82</td><td>149.80</td><td>28.73</td></tr><tr><td>BV30_C_H</td><td>30</td><td>30</td><td>1</td><td>36.10</td><td>4341.17</td><td>86.82</td><td>30.34</td></tr><tr><td>BV3_C_H</td><td>3</td><td>3</td><td>1</td><td>0.64</td><td>2493.44</td><td>49.87</td><td>31.73</td></tr><tr><td>BV300_C_L</td><td>300</td><td>300</td><td>0.01</td><td>0.31</td><td>39.46</td><td>0.79</td><td>0.35</td></tr><tr><td>BV30_C_L</td><td>30</td><td>30</td><td>0.01</td><td>0.51</td><td>24.67</td><td>0.49</td><td>0.36</td></tr><tr><td>BV3_C_L</td><td>3</td><td>3</td><td>0.01</td><td>0.88</td><td>9.53</td><td>0.19</td><td>0.37</td></tr><tr><td>BEV300_C_H</td><td>300</td><td>300</td><td>1</td><td>1.37</td><td>7587.78</td><td>151.75</td><td>28.70</td></tr><tr><td>BEV30_C_H</td><td>30</td><td>30</td><td>1</td><td>9.17</td><td>5893.80</td><td>117.88</td><td>29.70</td></tr><tr><td>BEV3_C_H</td><td>3</td><td>3</td><td>1</td><td>0.63</td><td>2502.28</td><td>50.84</td><td>31.91</td></tr><tr><td>BEV300_C_L</td><td>300</td><td>300</td><td>0.01</td><td>0.31</td><td>39.46</td><td>0.79</td><td>0.35</td></tr><tr><td>BEV30_C_L</td><td>30</td><td>30</td><td>0.01</td><td>0.57</td><td>29.46</td><td>0.59</td><td>0.36</td></tr><tr><td>BEV3_C_L</td><td>3</td><td>3</td><td>0.01</td><td>0.37</td><td>9.01</td><td>0.18</td><td>0.37</td></tr></tbody></table>

*Note.* Same as Table 1.

context that modulates the gas content and galaxy structure, via gas accretion from filaments and mergers.

## 3.1 Flux accretion simulations

We present here our two reference simulations, characterized by a high-resolution grid and non-parametrized accretion. Despite the different initial gas fractions, the hydrodynamical evolution of both models follows a common picture.

Fig. 1 shows a typical accretion event: cold, star forming gas is accreting along the equatorial plane, while the BAL wind begins to inject mass, momentum and energy in the ISM along the polar directions. Accretion continues until the bipolar outflow creates a hot, elongated expanding bubble of size ≈100 pc in the direction of least resistance. BH accretion is thus suppressed (Fig. 2). However,

AGN feedback is not powerful enough to evacuate all the gas present in the galaxy, and after ~10⁴–10⁵ yr the external cold gas is able to penetrate the hot bubble and to accrete on the BH, thus producing an irregular accretion-feedback-hot bubble cycle that lasts for the entire evolution and self-regulates the BH accretion. Fig. 3 (green line) shows even more clearly the effects of the erratic accretion and AGN feedback cycle: BH accretion rate and star formation rate (SFR) time evolution are highly fluctuating (green line, note that this is a moving average over 1 Myr) during the entire evolution. Remarkably, the BH self-consistently limits its own accretion rate to sub-Eddington values via AGN feedback (note that in the flux simulations we do not employ any capping on the BH accretion rate), producing a BH mass increment of ~4 × 10⁵ M☉ (see Table 2).

The top-left panel of Fig. 3 shows the temporal evolution of BH mass, BH accretion rate, stellar mass and SFR for the flux
<figure><img src="image_4.png" /><figcaption>**Figure 1.** From left to right: meridional section of density (in cm⁻³, left-hand panel), temperature (in K, middle panel), SFR density (in M☉ yr⁻¹ kpc⁻³) and radial velocity (in km s⁻¹, right-hand panel) for the flux accretion simulation with n₀ = 0.01 cm⁻³ at t = 38.7 Myr. The figure presents a representative accretion event: the cold, infalling gas is feeding the BH, while the BAL wind begins to inject mass, energy and momentum in the polar direction.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Same quantities as Fig. 1, but at the end of the accretion event, at *t* = 48.7 Myr. The BAL wind created a central hot bubble, temporarily shutting down any accretion on the BH.</figcaption></figure>

simulation having $n_0 = 10^{-2} \text{ cm}^{-3}$. After a few cooling times, a cold, infalling shell is formed in both simulations at $r \gtrsim 1$ kpc. The shell quickly fragments due to the Rayleigh–Taylor instability, producing a multiphase ISM characterized by a hot, diffuse component plus clumpy cold gas. This leads to a clumpy, cold accretion that efficiently feeds the BH (for a comparison between 1D and 2D simulations see Novak et al. 2011).

The counterpart at high gas density (top-right panel in Fig. 3) shows a much higher SFR, with a final stellar mass about two orders of magnitude higher, but, strikingly, the BH accretion rate and global growth are not boosted by the higher gas fraction, staying stable at a fraction about few per cent of the Eddington rate. AGN feedback seems to be more effective at regulating BH growth than star formation in the galaxy.

## 3.2 Bondi simulations with a high-resolution grid

Moving from flux accretion to the Bondi accretion runs, we show here the set of simulations employing the high-resolution grid. BH accretion is computed from equation (1) having an accretion radius $r_{\text{acc}} = 3, 30$ and $300$ pc (see Section 2.3), and the accretion radius is resolved respectively by 23, 98 and 208 grid elements in the radial direction.

The top-left panel in Fig. 3 presents the time evolution of the Bondi models having $n_0 = 10^{-2} \text{ cm}^{-3}$. The most striking feature is the apparent trend of $\Delta M_{\text{BH}}$ with $r_{\text{acc}}$, in contrast with the usual assumptions in the literature that an unboosted Bondi formalism underestimates the accretion rate. The larger the accretion radius,
Bondi mass weighted simulations

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Flux and Bondi mass weighted models evolution with time. Each panel is divided in four quadrants showing BH mass (top left), BH accretion rate (top right), stellar mass formed (bottom left) and SFR (bottom right) as a function of time. BH accretion rate and SFR are shown as a moving average over 1 Myr. The top (bottom) row shows high (low) resolution simulations, while low (high) gas density simulations are presented in the left (right) column. The 300 pc run in the top-left panel stopped due to the timestep becoming too short. The flux accretion run never exceeds the Eddington limit. To exclude the formation time of the gas cold shell, all the plots have been shifted in time so that the cold accretion mode (marked when $\dot{M}_{BH}/\dot{M}_{Edd}$ reaches the value of 2 for the first time) takes place at 2 Myr.</figcaption></figure>

the more the BH mass accretion is overestimated with respect to flux accretion (see also Table 1).

A common feature of all the runs employing the high-resolution grid is the concurrent formation of the cold shell at the same time of the flux simulation. The BH in the 300 pc run (red line in the top-left panel of Fig. 3) begins to accrete when the cold shell falls inside the accretion radius, reaching immediately a state of super-Eddington accretion. This produces a strong BAL wind that reduces the amount of cold gas inside $r_{acc}$ to $6 \times 10^6 M_{\odot}$. Following the

same accretion-feedback-hot bubble cycle of the flux simulation, a hot expanding bubble is formed, and star formation is quenched. This cycle holds until the BAL wind is finally able to evacuate all the cold gas from the innermost 10 kpc at $t = 32$ Myr. From then onwards, no cold gas is formed and the BH is in a hot, sub-Eddington accretion mode.

The 30 pc run shows a similar evolution of $\Delta M_{BH}$ with respect to the 300 pc simulation; $\dot{M}_{BH}$ fluctuates from super-Eddington accretion to strongly sub-Eddington, while $\dot{M}_{BH}$ in the 3 pc case is
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Time evolution of the BH luminosity over gas mass inside the accretion radius (in cgs units) for high-resolution, low-density simulations (top-left panel of Fig. 3).</figcaption></figure>

oscillating at higher frequency, and accretion remains overall lower than in the two previous cases (the Bondi 3 pc model accretes at super-Eddington rate at $t \simeq 2$ Myr, but this phase is not visible in figure due to the moving average).

AGN feedback is much more effective for smaller $r_{\text{acc}}$ since the volume to be cleared from cold material is smaller. The amount of cold gas inside $r_{\text{acc}}$ for the 30 and 3 pc runs is smaller than the 300 pc simulation, having less than $10^4$ and $10 M_{\odot}$ of cold gas inside $r_{\text{acc}}$ for the 30 and 3 pc model, respectively (note that in the last case the BH is mainly in a state of hot accretion). In the 300 pc case, accretion proceeds until $\dot{M}_{\text{BH}}$ becomes so large that feedback sweeps away all the cold material in one very strong event, producing a fast wind and fierce heating. Instead, in the 30 and 3 pc runs feedback is more effective in removing the cold gas at $r < r_{\text{acc}}$ at early times; hence $\dot{M}_{\text{BH}}$ is prevented from growing to high values.

This is the physical reason of the unexpected positive correlation between accreted mass and accretion radius: with large accretion radii the AGN feedback has to sweep away or heat a larger amount of gas to prevent accretion, thus, the feedback efficiency depends on $r_{\text{acc}}$. In Fig. 4, we quantify the efficiency with the ratio between $L_{\text{BH}}$, which regulates AGN feedback, and the gas mass inside $r_{\text{acc}}$, which has to be displaced or heated to stop the accretion. A clear trend with $r_{\text{acc}}$ is present: even if the 3 pc simulation has the least luminous AGN (and thus the weaker feedback), it shows the highest ratio of $L_{\text{BH}}/M_{\text{gas}}$. AGN feedback is more effective in clearing the accretion sphere for simulations with small $r_{\text{acc}}$.

The trend of BH accreted mass and AGN feedback strength with $r_{\text{acc}}$ at fixed resolution affects also the amount of gas mass converted into stars. On the one hand in simulations with small $r_{\text{acc}}$ feedback limits BH accretion therefore more cold material outside $r_{\text{acc}}$ is left to fuel star formation. On the other hand, for the largest $r_{\text{acc}}$ (the 300 pc run), even the cold gas at large radii is heated at $t \simeq 3$ Myr by the very strong AGN feedback, prompting a global reduction of cold gas mass and star formation.

To complete the set of high-resolution simulations, we present in Fig. 3 (top-right panel) the Bondi runs having $n_0 = 1 \text{ cm}^{-3}$. They resemble their low-density counterparts, in particular the BH mass increases with increasing $r_{\text{acc}}$. The physical mechanism is the same, at large $r_{\text{acc}}$ the AGN is less efficient in preventing cold gas accretion, producing BH masses up to 100 times the final mass of the flux BH simulation. The stellar mass shows the opposite trend: the larger $r_{\text{acc}}$, the smaller $\Delta M_*$, due to larger AGN luminosities at larger $r_{\text{acc}}$. The effects of AGN feedback on SFR are minimized in the high-density case, because star formation occurs mainly in the massive outer cold shell, at kpc-scales, where gas is sufficiently far that the heating of the AGN feedback is unimportant, but the density is high enough to form stars. The introduction of the Eddington limit does not affect the SFR, while the BH final mass is reduced up to a factor of $\simeq 90$ for the 300 pc case (see Table 1).

### 3.3 Bondi low-resolution grid simulations

In this section, we present the simulation families employing progressively low-resolution hydrodynamical grids, with central resolution of 3, 30 and 300 pc. For each grid, $r_{\text{acc}}$ is resolved by one cell (in the radial direction) therefore the mechanical feedback is injected in an area that grows with $r_{\text{acc}}$.

#### 3.3.1 Simulations without AGN feedback

Given the highly non-linear behaviour of the hydrodynamics, and the complex interplay between accretion and feedback, we first performed a subset of simulations without AGN feedback (see Table 3) by setting $\epsilon_{\text{EM}}$ and $\epsilon_{w0} = 0$ in equations (4) and (10).

As in the high-resolution runs, a cold shell is formed concurrently in all the models due to the cooling of the hot halo. However, the cold shells fragment in BE3_C_L and BE30_C_L, producing cold clumps that have an easier way to the BH than the main shell body. Instead, in the BE300_C_L case the shell proceeds towards the centre in its integrity due to the lower resolution, delaying accretion with respect to the former model.

The lack of feedback has a major impact on the global evolution of all the systems. As soon as the cold shell falls inside the accretion radius, a cooling flow is established until all the cold material present in the galaxy has been accreted. The abundance of cold gas and the absence of a powerful heating mechanism prevent the formation of the accretion-feedback-hot bubble cycle, producing a smooth evolution of BH accretion rate and a steeper stellar density profile in the central 100 pc.

Fig. 5 presents the time evolution of $\rho_{\infty}$ and $c_{s\infty}$ at $r < r_{\text{acc}}$, the trend usually assumed when $\alpha$ boosting is applied: the higher the resolution, the denser and colder the gas. This can be explained by the fact that at $r \gtrsim r_{\text{acc}}$ the gas is better resolved in high-resolution runs. As a consequence, smaller and smaller clumps of gas can be resolved during the infall, thus producing denser and colder clouds that prompt higher values of $\rho_{\infty}$ when they cross $r_{\text{acc}}$ and ending up with higher Bondi accretion rates.

#### 3.3.2 Simulations with AGN feedback

The introduction of AGN feedback entirely modifies the scenario of the previous section. In Fig. 6, we can contrast the (mass weighted) density and sound speed that enter the Bondi formalism with and without AGN feedback. The mean gas temperature inside the accretion radius is higher for smaller resolutions, while the gas densities are all very similar, after the cold gas shell has reached $r_{\text{acc}}$. In
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Mass weighted density and sound speed inside $r_{\text{acc}}$ for all the models without AGN feedback (see Table 3). Top and bottom panels refer to low- and high gas density simulations, respectively.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Mass weighted density and sound speed inside $r_{\text{acc}}$ for the B3\_C\_L, B30\_C\_L and B300\_C\_L model (top row), and B3\_C\_H, B30\_C\_H and B300\_C\_H model (bottom row). The corresponding simulations employing Eddington limited accretion present the same evolution.</figcaption></figure>

the 3 pc resolution run the density is much lower than in the case without AGN feedback, while in the 300 pc runs the densities are comparable. The physical reason can be again identified in the dependence of the AGN relative strength on the injection radius: the smaller the volume into which AGN feedback is injected, the more efficient it becomes, i.e. the energy input per unit gas mass in the accretion region is higher the smaller is $r_{\text{acc}}$ (cf. Fig. 4). We can conclude that the common assumption that low-resolution grids lead to an underestimate of $\dot{M}_{\text{BH}}$ is not verified in presence of resolved and realistic AGN feedback.

The BH and stellar evolution of the low-density, low-resolution runs are presented in the bottom-left panel of Fig. 3, showing the same trends of the high-resolution runs, due to the analogous physical evolution of the gaseous phase: the cold shell is formed, it accretes on the BH<fnref n="3" /> and the accretion-feedback-hot bubble cycle typical of the high-resolution runs is established. The 3 pc resolution run shows high-frequency oscillations in $\dot{M}_{\text{BH}}$, where the 30 and 300 pc runs present a much smoother accretion history. In all cases AGN feedback is able to halt the accretion on short time-scales.

The trend of $\Delta M_{\text{BH}}$ increasing with the size of $r_{\text{acc}}$ is in contrast with respect to the no AGN simulations, while it matches the scenario of high-resolution runs (Fig. 3), where $\Delta M_{\text{BH}}$ increases with the accretion radius. The interplay between BH feeding and AGN feedback at different resolutions leads to an accreted BH mass overestimated up to ~30 times with respect to the flux simulation, while stellar mass can be underestimated by a factor of ~15.

The same trend of $\Delta M_{\text{BH}}$ increasing at decreasing resolution appears in the high-density, low-resolution runs shown in the bottom-right panel of Fig. 3 and in the bottom panel of Fig. 6. Here, the mass accretion becomes largely super-Eddington in the 30 and 300 pc runs due to the large gas mass present in the accretion radius. On the other hand, the cold shell, while falling towards the BH, is heated by the AGN radiation. In the 3 pc runs (both at low and high density), the accretion is mainly sub-Eddington, but at larger $r_{\text{acc}}$, $L_{\text{BH}}$ is high enough to heat the cold gas still outside the accretion radius, thus keeping low the SFR.

With degrading grid resolution, at early times, before the BH start accreting and AGN feedback to operate, star formation is higher in the better resolved grid, as densities and temperatures are respectively higher and lower. As soon as AGN feedback is unleashed, star formation gets highly suppressed.

<footnote n="3">³ Note that the B300\_C\_L shows a delayed accretion like its counterpart without AGN feedback due to resolution effects, see Section 3.3.1.</footnote>
Bondi volume weighted simulations

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Same as Fig. 3, for the Bondi $n_0 = 10^{-2}$ and $1 \text{ cm}^{-3}$ simulations employing volume weighted quantities in equation (1) (left and right figure, respectively).</figcaption></figure>

3.4 Volume weighted Bondi runs

An alternative to mass weighting to compute $\rho_\infty$ and $c_{s\infty}$ is to use
volume weighted quantities (e.g. Vogelsberger et al. 2013, 2014a,
and SPH codes employing kernel weighted quantities). We present
here a family of high-resolution simulations employing a volume
weighted average of $\rho_\infty$ and $c_{s\infty}$ in equation (1). All the remaining
parameters (e.g. the accretion radius $r_{\text{acc}}$), initial conditions and
grids are unchanged from Section 3.2.

Fig. 7 shows the BH and stellar mass evolution, compared with the
corresponding flux simulations, for both the high- and low-density
runs. In both cases, at variance with the mass weighted simulations

of Fig. 3 (top panels), the accreted mass is smaller for large $r_{\text{acc}}$, since the computed Bondi accretion rate is sensitive to volume effects through $\rho_\infty$ and $c_{s\infty}$, while at small radii the simulations converge to the flux runs. This can be explained by the hot gas naturally having a filling factor larger than the cold gas inside the accretion radius therefore the hot phase is the main actor in volume weighted quantities. As a consequence, in the BV300\_S\_L (where the hot gas filling factor is larger and thus its influence is maximized) the BH is in a state of hot accretion for most of the evolution, until a significant volume of cold gas crosses $r_{\text{acc}}$, prompting accretion (still sub-Eddington, see $t \simeq 2$ Myr in Fig. 7, left-hand panels); the subsequent feedback then removes the remaining cold gas and
restores the hot accretion mode. As expected, in the high-density counterpart the deviation from the flux accretion at large $r_{\text{acc}}$ is minimized by the large amount of gas that can cools down, thus enhancing the accretion rate by increasing the cold gas filling factor.

SFR and $\Delta M_*$ in low-density cases show a weaker dependence on $r_{\text{acc}}$ compared to the mass weighted simulations (cf. top rows of Figs 3 and 7), since the main situ of star formation is localized in the cold shell, which is formed at $r \gtrsim 1$ kpc. Moreover, all the simulations spend almost all their evolution at low accretion rates (with mean values of $\dot{M}_{\text{BH}}/\dot{M}_{\text{Edd}}$ ranging from $4 \times 10^{-4}$ to $10^{-2}$), thus the corresponding AGN feedback is unable to significantly affect star formation. The similarities between runs with different $r_{\text{acc}}$ are even less noticeable in the high gas density case (Fig. 7, right-hand panels), due to the more massive cold shell, where the majority of star formation takes place.

The scenario is modified when the grid is coarser. On the one hand, in low-resolution runs without AGN feedback the evolution of gas density and temperature inside the accretion radius shows the same trend of mass weighted simulations (cf. left-hand and right-hand panels of Fig. 5), with lower density values attributable to volume weighted mean. On the other hand, when AGN feedback is present, the models with $r_{\text{acc}} = 3$ present a global evolution of $\Delta M_{\text{BH}}$ and $\Delta M_*$ similar to the flux accretion runs. Instead, at larger $r_{\text{acc}} = 300$ feedback is less efficient in stopping the accretion, producing a larger accretion rate than in the high-resolution runs. This phenomenon compensates the reduced accretion rate typical of volume weighted Bondi accretion, ending in a BH mass of the same order of magnitude of the flux run for the $r_{\text{acc}} = 300$ pc, whereas in the mass weighted case it was overestimated by a factor of 100. The BV30\_C\_H presents a peculiar evolution, characterized by an oscillating BH accretion rate at $t > 20$ Myr. This is due to a large blob of cold gas entering in the accretion radius, reducing $c_s$ and thus increasing the accretion rate. In addition, the general time evolution of BH and stellar mass is smoother with respect to the low-resolution mass weighted models, being dominated by hot gas, and lacking the strong, super-Eddington peaks triggering strong feedback phases.

## 3.5 Eddington rate limiter?

In the runs discussed above, the BH can grow theoretically without limitations. However, in cosmological simulations BH accretion is usually capped at the Eddington luminosity, by assuming $\dot{M}_{\text{BH}} = \min(\dot{M}_{\text{B}}, \dot{M}_{\text{Edd}})$. Thus, we simulate the same set of galaxies presented in the previous sections by employing Bondi accretion (mass and volume weighted), capped at the Eddington rate, i.e. assuming a maximal accretion rate corresponding to the Eddington luminosity for a radiative efficiency of 10 per cent. For the majority of the models this additional constraint does not significantly modify the global picture (see Tables 1 and 3 and 4). The reduced accretion rate produces a less energetic feedback that is not able to evacuate the accretion volume, maintaining the BH accretion rate at the Eddington limit for the entire evolution. Smaller BH masses are thus produced with respect to the Bondi runs, with the solely exception of the BE300\_S\_L model, due to the fact that in its Bondi counterpart spends half of its evolution in a sub-Eddington phase (cf. with the red line in top-left panel of Fig. 3). Given the low accretion rates involved, the volume weighted simulations capped at the Eddington rate reported in Table 4 do not show any significant difference with respect to the Bondi accretion runs (except for the BV30\_C\_H run).

## 4 DISCUSSION AND COMPARISON WITH PREVIOUS WORKS

Our analysis of the Bondi accretion in high- and low-resolution simulations has revealed how BH feeding and AGN feedback influence each other, and that the exact implementation affects the results. Usually, cosmological simulations justify the adoption of a boost parameter $\alpha$ to artificially increase Bondi accretion with two main arguments: (1) large-scale simulations do not have the necessary resolution nor the physics (e.g. molecular cooling) to resolve the cold gas, thus the gas temperature near the AGN may be overestimated and (2) the Bondi radius is normally unresolved therefore the ISM density may be underestimated.

We show that the above assumptions, while they remain true in simulations without AGN feedback, are not always verified in its presence. Our reference model employing flux accretion with a 0.1 pc resolution grid display moderately low BH accretion. This is due to feedback, which is able to evacuate the BH surroundings, halt accretion and keep the BH in a self-regulation regime. A cycle of cooling, accretion, feedback and hot bubble is created, that ultimately leads to a cold, irregular accretion. This accretion-feedback loop is not new in the literature: using 1D simulations of Ciotti & Ostriker (1997, 2001, 2007) found that the flickering of the AGN is due to cold shells of gas accreting on the BH. Higher frequency oscillations in $L_{\text{BH}}$ are produced in 2 and 3D simulations, where the gas cold shell naturally fragments, producing a rain of cold gas clumps that can dramatically increase $\dot{M}_{\text{BH}}$ (this accretion mode has been dubbed Cold Chaotic Accretion in Gaspari et al. (2013) and Gaspari, Temi & Brighenti 2017, see also Ciotti et al. 2017).

We explored the parameter space by varying the accretion algorithm, the accretion radius and the numerical resolution of the simulations. Comparison between different models can be performed by either keeping constant the resolution, and comparing simulations with different $r_{\text{acc}}$, or keeping fixed $r_{\text{acc}}$ and considering different resolutions, in order to isolate pure resolution effects. Finally, one could consider comparing simulations with different resolutions and $r_{\text{acc}}$, but with a specified accretion algorithm. The results are summarized in Fig. 8, and we will discuss them in the following.

### 4.1 Comparison at constant resolution with various $r_{\text{acc}}$

From Fig. 9, by keeping fixed the resolution (filled and empty circles), a clear trend emerges: in the case of mass weighted accretion (left panel), $\Delta M_{\text{BH}}$ increases and $\Delta M_*$ decreases at increasing $r_{\text{acc}}$. When $r_{\text{acc}}$ is much larger than the Bondi radius (red filled circles) the BH mass is overestimated up a factor of 100. In volume weighted simulations, instead, the trend is opposite: $\Delta M_{\text{BH}}$ and $\Delta M_*$ are respectively under and overestimated with respect to the flux case, at increasing $r_{\text{acc}}$.

This difference is due to the complex interplay between accretion and feedback. We show that the feedback efficiency in clearing the accretion region is a function of the accretion radius. On the one hand, in the mass weighted case, which is dominated by accretion of cold gas, at large $r_{\text{acc}}$ AGN feedback is not effective in evacuating the accretion region, since the mass accumulated is too large to be pushed outside $r_{\text{acc}}$ or heated up. This leads to a continuous accretion that proceeds until $L_{\text{BH}}$ becomes so large that feedback generates a hot bubble that shuts off accretion. However, at this point the BH has already increased its mass by a considerable amount with respect the flux accretion run. In this regime, the BH accretion rate should be de-boosted to recover the evolution of the flux accretion simulation.
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Time evolution for the 300 pc (left) and 3 pc Bondi (right) runs having $n_0 = 10^{-2} \text{ cm}^{-3}$ (top) or $n_0 = 1 \text{ cm}^{-3}$ (bottom), compared with flux accretion (simulations capped at the Eddington rate are not included).</figcaption></figure>

On the other hand, the volume average is essentially dominated by the hot and cold gas filling factors and not by their mass. Thus, at large $r_{\text{acc}}$, where the filling factor of the hot gas dominates over that of the cold gas, the BH is in a state of hot accretion, leading to a lower BH final mass with respect to the flux case (see empty circles in Fig. 10). In this regime, that is what is normally assumed in literature, boosting the accretion rate when large $r_{\text{acc}}$ are used goes in the right direction.

This picture is strengthened when an Eddington capping is applied to Bondi accretion: in the mass weighted case AGN feedback is further hindered by the capping, leading to continuous accretion for the entire evolution, whereas in the volume weighted case there is no sensitive difference.

## 4.2 Comparison at constant $r_{\text{acc}}$ with various resolutions

In considering the pure effect of resolution, we can compare high- and low-resolution simulations sharing the same $r_{\text{acc}}$ (represented by squares and circles of the same colour in Fig. 9). In the mass weighted case no clear trend is visible,<fnref n="4" /> except at low $r_{\text{acc}}$ (black symbols), where $\Delta M_{\text{BH}}$ is lower at low resolution, although the two simulations are very similar. Surprisingly, in volume weighted runs $\Delta M_{\text{BH}}$ increases at decreasing resolution with constant $r_{\text{acc}}$, showing an evolution similar to the flux case (with the notable exception of the BV30_C_H).

In principle one would expect that the lower the resolution, the worse the ability of the simulation to capture the properties of the

gas feeding the BH, located, in reality, on sub-pc scales. This result is caused by a combination of two factors cancelling each other out. On the one hand, volume weighting tends to *underestimate* BH accretion, because of the hot gas filling factor dominance, and this behaviour is enhanced at large $r_{\text{acc}}$ (see also Section 4.1). On the other hand, at low resolution and large $r_{\text{acc}}$, AGN feedback is less effective in clearing out cold gas. The combination of these effects produces reasonable BH and stellar masses, although stellar mass may be overestimated in gas-rich galaxies (see Fig. 8). This surprising compensation may explain why, statistically, BH and galaxy properties are well reproduced in large-scale, relatively low-resolution cosmological simulations (see the discussion in Volonteri et al. 2016).

Finally, as expected, there is convergence between Bondi and flux accretion when the former is calculated using information in the immediate surroundings of the BH, at a distance of the scale of the Bondi radius. As shown in Figs 8 and 10, the spread in BH and galaxy mass is within a factor of 2 for all runs with an accretion radius of 3 pc, regardless of the numerical implementation and of the galaxy gas content.

## 4.3 Comparison with previous works

Our results can be compared with other recent papers on BH accretion and feedback using high-resolution simulations. Spherical Bondi accretion on supermassive BH has been studied in Barai et al. (2011), by performing high-resolution simulations with the 3D code GADGET-3. The spatial scales involved range from 0.1 to 20–200 pc, with a softening of 0.005–0.02 pc (gas self-gravity is not considered), thus the Bondi radius is resolved. Compton heating and cooling are included, but momentum feedback is not. The BH accretes all particles crossing the minimum radius of the

<footnote n="4">⁴ Note that the BV30_C_L simulation shows a much smaller $\Delta M_*$ with respect to what the trend of the other models suggest due to a very early, and short-lived phase of strong AGN feedback, suppressing star formation. See Fig. 3, bottom-left panel, $t \sim 2$ Myr.</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Accreted BH mass versus stellar mass for the simulations of Tables 1, 2 and 4, except Eddington limited runs. Circles and squares represent runs with standard and coarse grid, respectively, while filled and open symbols refer to mass and volume weighted simulations. Crossed symbols are Eddington limited runs. The accretion radius is mapped by colours: red, blue and black represent $r_{\text{acc}} = 300, 30$ and $3$ pc. Green stars represent the flux simulations (which employ the standard $0.1$ pc resolution grid). All the simulations with $\Delta M_{\star} < 10^{8} \text{ M}_{\odot}$ (vertical dotted line) have $n_{0} = 10^{-2} \text{ cm}^{-3}$, $n_{0} = 1 \text{ cm}^{-3}$ otherwise. The models converge towards the flux runs when the Bondi radius is progressively better resolved. Mass weighted simulations always overestimate $\Delta M_{\text{BH}}$, high-resolution volume weighted simulations underestimate it, while for low-resolution volume runs a clear trend is not present.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Accreted BH mass versus stellar mass for all the simulations of Tables 1, 2 and 4. Symbols are the same as in Fig. 9.</figcaption></figure>

simulation and the BH luminosity is artificially kept constant during the evolution. While Barai et al. (2011) clearly show that an SPH code can follow the steady evolution of an adiabatic Bondi solution without major issues, in presence of radiative heating and cooling the picture is modified. At intermediate $L_{\text{BH}}$ non spherical cold clumps can be formed via thermal instability, while at high $L_{\text{BH}}$ an outflow is created, reducing the accretion rate of a factor of $10^3$ (see their fig. 13).

The formation of a multiphase gas is further investigated in Barai et al. (2012) using the same setup of the previous work. They

confirm that formation of cold gas clumps in a spherical setup is a consequence of thermal instabilities acting on small perturbations of the initial conditions. These cold filaments accrete on the BH and trigger large fluctuations on short time-scales in the accretion rate (like our reference model). The mass inflow rate ratio between the cold and the hot phase is calculated to constrain the boost parameter $\alpha$ commonly employed in large-scale simulations. They find a value $\le 1-5$ at the Bondi radius, much lower respect the usual factor $\alpha = 100$ used in cosmological simulations, meaning that they do not find an important unresolved fraction of gas hidden in the cold phase. Finally, they argue that $\alpha$ can be resolution dependent. We confirm their results and extend them in the case of more realistic AGN feedback, galaxy model and in presence of star formation. The additional feedback mechanisms in our work are able to produce a multiphase gas in almost every circumstance, while they obtain a multiphase medium and a non-smooth BH accretion rate only in case of high BH luminosity. In addition, we find that the Bondi rate is not a good proxy for accretion when the Bondi radius is not resolved.

A comparable work regarding explored spatial scales and resolutions but with a more realistic feedback is presented in Hopkins et al. (2016). They perform 3D simulations of the first 100 pc around an isolated BH with the GIZMO code, assuming a force softening of $\epsilon = 0.02$ pc, including gas self-gravity, rotation, radiative cooling down to 10 K and star formation. The Compton heating/cooling and AGN feedback are the same of the present work.<fnref n="5" /> Moreover, the BH accretes immediately any gas particle

<footnote n="5">⁵ The adopted AGN feedback efficiencies are higher in Hopkins et al. (2016). For our models $\eta = \dot{M}_{\text{out}}/\dot{M}_{\text{BH}} = 0.4$ and $v_{\text{w}} = 10^4 \text{ km s}^{-1}$, while Hopkins et al. (2016) adopt much higher values of $\eta$ (up to 60, named $\beta$ in their table 1), and in general lower outflow velocities (500–5000 km s$^{-1}$).</footnote>
having an apocentric radius <2.8ε. They compare this accretion scheme (that can represent the meshless counterpart of the flux accretion in grid codes) with two different types of Bondi approaches, namely $\dot{M} = 4\pi G^2 M_{BH}^2 \rho_{\text{gas}} c_s^{-3}$ and $\dot{M} = 4\pi G^2 M_{\text{enc}}^2 (<R)\rho_{\text{gas}}(c_s^2 + (V_{\text{gas-BH}}^2))^{-3/2}$, where the latter takes into account the mass outside the BH as a source of gravity for the gas outside the centre of the galaxy (however, it is not made explicit how $\rho_{\text{gas}}$ and $c_s$ are calculated). Both schemes fail to reproduce the correct inflow rate in simulations without AGN feedback by overestimating the BH accretion rate of a factor of $10^8$ and $10^3$, respectively (see fig. 6 in Hopkins et al. 2016, note that the Bondi accretion is computed, but not used to calculate feedback). The accretion of cold and molecular gas is regulated by gravitational torques in gaseous disc, while Bondi accretion assumes a spherical gas distribution with no rotation, thus the low $c_s$ of the rotating gas produces a large accretion rate in the Bondi formalism.

A direct comparison between the different models in Hopkins et al. (2016) with our results is not straightforward, due to the low rotation regime in our 2.5D simulations (see Section 2.4). However, from the point of view of BH accretion, these two studies are complementary: accretion in Hopkins et al. (2016) is completely regulated by the gravitational torques of the gas+stars disc, thus confirming and quantifying the inadequacy of the Bondi accretion in presence of rotation. On the other hand, we show that the Bondi accretion does not predict the right accretion rate even in absence of angular momentum. This can depend on resolution, feedback, volume and algorithm employed to calculate density and sound speed in the Bondi formula, producing both an over- and underestimate of the accretion rate. When BH winds are activated, Hopkins et al. (2016) find that AGN feedback can evacuate the polar regions with high velocity outflows ($5-30 \times 10^3 \text{ km s}^{-1}$), suppressing star formation and producing a BH accretion rate with an extreme time variability, in the same way mechanical feedback in our reference model is able to clear the innermost 100 pc from cold gas, thus producing an extremely fluctuating BH accretion rate that is maintained below the Eddington limit in a self-consistent way.

Curtis & Sijacki (2015) attempt to bridge the scale of the Bondi radius to the resolutions typical of galactic scale simulations with a new super-Lagrangian refinement in the Voronoi-mesh AREPO code. The cells inside the smoothing length of the BH are refined almost down to the Bondi radius, and the Bondi rate is then computed from the mass weighted mean of density and sound speed over the cells. They apply this novel refinement scheme to 3D simulations of an isolated Milky Way like object, by employing different types of feedback (thermal and/or momentum, both isotropic and non-isotropic). The simulations with super-Lagrangian refinement show somewhat lower rates of gas accretion than those without (note that this resolution test is performed with isotropic thermal feedback only), due to the higher temperature generated by feedback in the higher resolution run. While this result agrees with our work, the difference in the final BH mass between different refinements found in Curtis & Sijacki (2015) is small. We designed two tailored simulations with a setup similar to Curtis & Sijacki (2015). We adopted their accretion model with $r_{\text{acc}} = 700$ pc, a 1 pc resolution grid (see their fig. 7) and $n_0 = 0.01 \text{ cm}^{-3}$, $M_{\text{BH},0} = 10^6 \text{ M}_\odot$. The two runs employ a mass weighted Bondi accretion and a flux accretion, respectively. The results are in line with our main set of high-resolution simulations presented in Section 3.2. The mass weighted Bondi accretion run overestimates the BH accreted mass with respect to the flux accretion at the same resolution and at 0.1 pc resolution by a factor ≃5 and 25, respectively. The BH in the Bondi run accretes at Eddington rate until feedback is able to

evacuate the accretion radius from cold gas therefore its evolution is quite smooth, in contrast to the simulations using flux accretion. We find that the main trends resemble the Bondi mass weighted high-resolution runs presented in Section 3.2, where the BH accretion is always overestimated when compared to the flux accretion. This is due to the poor feedback efficiency in stopping the accretion, while in the flux accretion runs the BH feeding is much more irregular, and feedback is able to halt accretion more effectively. This can explain in part the BH mass rapid increase in fig. 12 of Curtis & Sijacki (2015) when momentum feedback only is employed. The scenario where AGN feedback is not effective in halting the accretion at large $r_{\text{acc}}$, due to the large amount of cold gas inside the accretion radius, is thus confirmed. The super-Lagrangian refinement presented in Curtis & Sijacki (2015) is promising, given its ability to dynamically resolve the region around BHs, and indeed in such a scheme the increased resolution can be used directly to apply a flux accretion to overcome the issues of Bondi accretion.

## 5 CONCLUSIONS

We employ 2.5D hydrodynamical simulations of an isolated galaxy to test when and whether Bondi accretion is a good proxy for BH feeding. In cosmological simulations the cold gas phase and the Bondi radius are not resolved, thus the ISM temperature is expected to be overestimated and its density underestimated. To correct for these effects, a boost parameter $\alpha$ is often employed to increase the BH accretion rate. The large amount of sub-grid recipes based on Bondi accretion present in the literature motivated us to study how well BH accretion and the subsequent AGN feedback on galactic scales can be captured.

We include AGN radiative and mechanical feedback together with a 0.1 pc resolution grid. Our simulation setup consists of an isolated galaxy with a numerical grid (in the reference case) spanning from 0.1 pc to 250 kpc, thus resolving the Bondi radius for all the temperatures. Our reference runs implement flux accretion on the innermost grid. We compare the reference case with simulations employing Bondi accretion at different resolutions, and with different tweaks of the implementation (mass versus volume weighting). The main results are summarized below:

(i) The reference runs show a self-regulated accretion that never exceeds the Eddington limit, without any capping explicitly included. The hydrodynamics is characterized by irregular and highly variable cycles of cooling, accretion, AGN feedback.

(ii) When mass weighting is employed to estimate gas density and sound speed in the Bondi formalism in high-resolution simulations (where the grid still resolves the Bondi radius), the BH accreted mass is overestimated up to a factor of 100. This is due to the complex interplay between feeding and feedback. The feedback efficiency in clearing the BH surroundings is a function of the accretion radius. At large accretion radii the AGN feedback is less effective in evacuating the accretion region, since the gas mass accumulated is too large to be swept away or heated up, leading to a smoother accretion and in the end to a larger BH mass.

(iii) In simulations with progressively lower resolution and larger $r_{\text{acc}}$, in absence of AGN feedback, the ISM temperature is overestimated and its density underestimated. However, the introduction of feedback in simulations with different resolutions leads to an overestimate of the BH accretion due to poor feedback efficiency at lower resolutions. Capping the accretion rate at the Eddington limit only enhances this trend since it does not allow the BH to reach the luminosity needed to effectively halt accretion.
(iv) Volume weighted Bondi accretion underestimates the BH accreted mass when the high-resolution grid is employed. BH accretion is dominated by hot gas, as a result of the volume weighted averaging. In this scenario, the use of a boost parameter is justified, although the exact value of the boost factor is resolution and problem dependent.

(v) Unexpectedly, reducing the resolution in volume weighted Bondi simulations, BH and stellar mass are similar to the 'flux' case. This occurs because both the accretion rate and AGN feedback are underestimated, and compensate each other.

(vi) In general, we show that convergence between Bondi and flux accretion is reached only when the accretion radius resolves the Bondi radius.

## ACKNOWLEDGEMENTS

We thank Luca Ciotti and Jeremiah Ostriker for useful discussion and comments, Greg Novak for providing an initial version of the code used here and the anonymous referee for helpful comments and contributions. AN thanks Tilman Hartwig, Muhammad Latif for interesting discussions, Rebekka Bieri and Mélanie Habouzit for help with the Horizon Cluster setup. AN and MV acknowledge funding from the European Research Council under the European Community's Seventh Framework Programme (FP7/2007-2013 Grant Agreement no. 614199, project 'BLACK'). This work has made use of the Horizon Cluster funded by Institut d'Astrophysique de Paris. We thank Stephane Rouberol for running smoothly this cluster for us.

## REFERENCES

Anglés-Alcázar D., Davé R., Faucher-Giguère C.-A., Özel F., Hopkins P. F., 2017, MNRAS, 464, 2840

Barai P., Proga D., Nagamine K., 2011, MNRAS, 418, 591

Barai P., Proga D., Nagamine K., 2012, MNRAS, 424, 728

Barai P., Viel M., Murante G., Gaspari M., Borgani S., 2014, MNRAS, 437, 1456

Bellovary J. M., Governato F., Quinn T. R., Wadsley J., Shen S., Volonteri M., 2010, ApJ, 721, L148

Bertin G., Lodato G., 2001, A&A, 370, 342

Bleuler A., Teyssier R., 2014, MNRAS, 445, 4015

Bondi H., 1952, MNRAS, 112, 195

Bondi H., Hoyle F., 1944, MNRAS, 104, 273

Booth C. M., Schaye J., 2009, MNRAS, 398, 53

Bryan G. L. et al., 2014, ApJS, 211, 19

Choi E., Ostriker J. P., Naab T., Johansson P. H., 2012, ApJ, 754, 125

Choi E., Ostriker J. P., Naab T., Oser L., Moster B. P., 2015, MNRAS, 449, 4105

Ciotti L., Ostriker J. P., 1997, ApJ, 487, L105

Ciotti L., Ostriker J. P., 2001, ApJ, 551, 131

Ciotti L., Ostriker J. P., 2007, ApJ, 665, 1038

Ciotti L., Ostriker J. P., 2012, in Kim D.-W., Pellegrini S., eds, Astrophysics and Space Science Library, Vol. 378, Hot Interstellar Matter in Elliptical Galaxies. Springer-Verlag, Berlin, p. 83

Ciotti L., Morganti L., de Zeeuw P. T., 2009, MNRAS, 393, 491

Ciotti L., Pellegrini S., Negri A., Ostriker J. P., 2017, ApJ, 835, 15

Curtis M., Sijacki D., 2015, MNRAS, 454, 3445

Curtis M., Sijacki D., 2016, MNRAS, 463, 63

Debuhr J., Quataert E., Ma C.-P., Hopkins P., 2010, MNRAS, 406, L55

Debuhr J., Quataert E., Ma C.-P., 2011, MNRAS, 412, 1341

Degraf C., Di Matteo T., Springel V., 2011, MNRAS, 413, 1383

DeGraf C., Di Matteo T., Khandai N., Croft R., Lopez J., Springel V., 2012, MNRAS, 424, 1892

Di Matteo T., Springel V., Hernquist L., 2005, Nature, 433, 604

Di Matteo T., Colberg J., Springel V., Hernquist L., Sijacki D., 2008, ApJ, 676, 33

Di Matteo T., Croft R. A. C., Feng Y., Waters D., Wilkins S., 2016, MNRAS, preprint (arXiv:1606.08871)

Dubois Y., Devriendt J., Slyz A., Teyssier R., 2012, MNRAS, 420, 2662

Dubois Y. et al., 2014, MNRAS, 444, 1453

Elahi P. J. et al., 2016, MNRAS, 458, 1096

Fabian A. C., 1999, MNRAS, 308, L39

Fabian A. C., 2012, ARA&A, 50, 455

Feng Y., Di-Matteo T., Croft R. A., Bird S., Battaglia N., Wilkins S., 2016, MNRAS, 455, 2778

Ferrarese L., Merritt D., 2000, ApJ, 539, L9

Gaspari M., Ruszkowski M., Oh S. P., 2013, MNRAS, 432, 3401

Gaspari M., Brighenti F., Temi P., 2015, A&A, 579, A62

Gaspari M., Temi P., Brighenti F., 2017, MNRAS, 466, 677

Gebhardt K. et al., 2000, ApJ, 539, L13

Gültekin K. et al., 2009, ApJ, 698, 198

Hayes J. C., Norman M. L., Fiedler R. A., Bordner J. O., Li P. S., Clark S. E., ud-Doula A., Mac Low M.-M., 2006, ApJS, 165, 188

Hensley B. S., Ostriker J. P., Ciotti L., 2014, ApJ, 789, 78

Hirschmann M., Dolag K., Saro A., Bachmann L., Borgani S., Burkert A., 2014, MNRAS, 442, 2304

Hopkins P. F., Quataert E., 2010, MNRAS, 407, 1529

Hopkins P. F., Quataert E., 2011, MNRAS, 415, 1027

Hopkins P. F., Torrey P., Faucher-Giguère C.-A., Quataert E., Murray N., 2016, MNRAS, 458, 816

Hoyle F., Lyttleton R. A., 1939, Proc. Camb. Phil. Soc., 35, 405

Inayoshi K., Haiman Z., Ostriker J. P., 2016, MNRAS, 459, 3738

Jaffe W., 1983, MNRAS, 202, 995

Kennicutt, R. C., Jr, 1998, ApJ, 498, 541

Khandai N., Di Matteo T., Croft R., Wilkins S., Feng Y., Tucker E., DeGraf C., Liu M.-S., 2015, MNRAS, 450, 1349

Kim J.-h., Wise J. H., Alvarez M. A., Abel T., 2011, ApJ, 738, 54

King A., 2005, ApJ, 635, L121

Kormendy J., 2016, in Laurikainen E., Peletier R. F., Gadotti D. A., eds, Astrophysics and Space Science Library, Vol. 418, Galactic Bulges. Springer International Publishing, Switzerland, p. 431

Kormendy J., Ho L. C., 2013, ARA&A, 51, 511

Korol V., Ciotti L., Pellegrini S., 2016, MNRAS, 460, 1188

Kurosawa R., Proga D., 2009a, MNRAS, 397, 1791

Kurosawa R., Proga D., 2009b, ApJ, 693, 1929

Kurosawa R., Proga D., Nagamine K., 2009, ApJ, 707, 823

Le Brun A. M. C., McCarthy I. G., Schaye J., Ponman T. J., 2014, MNRAS, 441, 1270

Li J., Ostriker J., Sunyaev R., 2013, ApJ, 767, 105

Lusso E., Ciotti L., 2011, A&A, 525, A115

Magorrian J. et al., 1998, AJ, 115, 2285

McCarthy I. G., Schaye J., Bird S., Le Brun A. M. C., 2017, MNRAS, 465, 2936

McConnell N. J., Ma C.-P., 2013, ApJ, 764, 184

Narayan R., Yi L., 1994, ApJ, 428, L13

Negri A., Ciotti L., Pellegrini S., 2014a, MNRAS, 439, 823

Negri A., Posacki S., Pellegrini S., Ciotti L., 2014b, MNRAS, 445, 1351

Negri A., Pellegrini S., Ciotti L., 2015, MNRAS, 451, 1212

Novak G. S., Ostriker J. P., Ciotti L., 2011, ApJ, 737, 26 (N11)

Novak G. S., Ostriker J. P., Ciotti L., 2012, MNRAS, 427, 2734 (N12)

Ostriker J. P., Choi E., Ciotti L., Novak G. S., Proga D., 2010, ApJ, 722, 642

Park K., Ricotti M., 2012, ApJ, 747, 9

Pelupessy F. I., Di Matteo T., Ciardi B., 2007, ApJ, 665, 107

Power C., Nayakshin S., King A., 2011, MNRAS, 412, 269

Prasad D., Sharma P., Babul A., 2015, ApJ, 811, 108

Robertson B., Bullock J. S., Cox T. J., Di Matteo T., Hernquist L., Springel V., Yoshida N., 2006, ApJ, 645, 986

Rosas-Guevara Y. M. et al., 2015, MNRAS, 454, 1038

Schaye J. et al., 2010, MNRAS, 402, 1536

Schaye J. et al., 2015, MNRAS, 446, 521

Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337
Sijacki D., Springel V., Di Matteo T., Hernquist L., 2007, MNRAS, 380, 877

Sijacki D., Vogelsberger M., Genel S., Springel V., Torrey P., Snyder G. F., Nelson D., Hernquist L., 2015, MNRAS, 452, 575

Silk J., Rees M. J., 1998, A&A, 331, L1

Springel V., 2010, MNRAS, 401, 791

Springel V., Di Matteo T., Hernquist L., 2005, MNRAS, 361, 776

Steinborn L. K., Dolag K., Hirschmann M., Prieto M. A., Remus R.-S., 2015, MNRAS, 448, 1504

Stone J. M., Norman M. L., 1992, ApJS, 80, 753

Teyssier R., 2002, A&A, 385, 337

Teyssier R., Moore B., Martizzi D., Dubois Y., Mayer L., 2011, MNRAS, 414, 195

Tremmel M., Karcher M., Governato F., Volonteri M., Quinn T., Pontzen A., Anderson L., 2016, MNRAS, preprint (arXiv:1607.02151)

Vogelsberger M., Genel S., Sijacki D., Torrey P., Springel V., Hernquist L., 2013, MNRAS, 436, 3031

Vogelsberger M. et al., 2014a, MNRAS, 444, 1518

Vogelsberger M. et al., 2014b, Nature, 509, 177

Volonteri M., Silk J., Dubus G., 2015, ApJ, 804, 148

Volonteri M., Dubois Y., Pichon C., Devriendt J., 2016, MNRAS, 460, 2979

Wurster J., Thacker R. J., 2013, MNRAS, 431, 2513

This paper has been typeset from a *TEX/LATEX* file prepared by the author.