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

THE EFFECT OF AGN FEEDBACK ON THE INTERSTELLAR MEDIUM OF EARLY-TYPE GALAXIES: 2D HYDRODYNAMICAL SIMULATIONS OF THE LOW-ROTATION CASE

LUCA CIOTTI¹, SILVIA PELLEGRINI¹, ANDREA NEGRI², AND JEREMIAH P. OSTRIKER³,⁴

¹Department of Physics and Astronomy, University of Bologna, via Ranzani 1, I-40127, Bologna, Italy

²CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98bis bvd Arago, F-75014 Paris, France

³Department of Astronomy, Columbia University, 550 W. 120th Street, New York, NY 10027, USA

⁴Princeton University Observatory, Princeton, NJ 08544, USA

Received 2016 July 25; revised 2016 November 15; accepted 2016 November 27; published 2017 January 16

ABSTRACT

We present two-dimensional hydrodynamical simulations for the evolution of early-type galaxies containing central massive black holes (MBHs), starting at an age of ≈2 Gyr. The code contains accurate and physically consistent radiative and mechanical active galactic nucleus (AGN) wind feedback, with parsec-scale central resolution. Mass input comes from stellar evolution; energy input includes Type Ia (SNIa) and II supernovae and stellar heating; star formation (SF) is included. Realistic, axisymmetric dynamical galaxy models are built solving the Jeans' equations. The lowest mass models ($M_{*} = 8 \ 10^{10} M_{\odot}$) develop global outflows sustained by SNIa heating, ending with a lower amount of hot gas and new stars. In more massive models, nuclear outbursts last to the present epoch, with large and frequent fluctuations in nuclear emission and from the gas ($L_{X}$). Each burst lasts ~$10^{7.5}$ years, during which cold, inflowing, and hot, outflowing gas phases coexist. The $L_{X}-T_{X}$ relation for the gas matches that of local galaxies. AGN activity causes positive feedback for SF. Roughly half of the total mass loss is recycled into new stars ($\Delta M_{*}$), just ≈3% of it is accreted on the MBH, the remainder being ejected from the galaxy. The ratio between the mass of gas expelled to that in new stars, the load factor, is ≈0.6. Rounder galaxy shapes lead to larger final MBH masses, $\Delta M_{*}$, and $L_{X}$. Almost all of the time is spent at very low nuclear luminosities, yet one quarter of the total energy is emitted at an Eddington ratio >0.1. The duty-cycle of AGN activity is approximately 4%.

Key words: galaxies: elliptical and lenticular, cD – galaxies: evolution – galaxies: ISM – quasars: supermassive black holes – X-rays: galaxies – X-rays: ISM

1. INTRODUCTION

The relationship of quasi-stellar object (QSO) activity at high z with the surrounding interstellar medium (ISM) of the host galaxy, including the possibility of triggering or quenching star formation (SF), is one of the currently most debated and still unsettled topics in the field of galaxy evolution. The massive black holes (MBHs) at the centers of massive elliptical galaxies were in place already when the universe was ≈1 Gyr old (e.g., Madau & Rees 2001; Alvarez et al. 2009; Wu et al. 2015), and fast, large scale, and massive outflows, driven by QSOs, are supposed to transform young, star-forming galaxies into “red and dead” spheroids (Sturm et al. 2011; Cano-Díaz et al. 2012; Faucher-Giguère & Quataert 2012; Zubovas & King 2012; Feruglio et al. 2015; King & Pounds 2015). But QSO and SF activities are not likely to suddenly stop at high z. Concerning the origin of QSO activity, it appears that the most luminous active galactic nucleus (AGN) phases, preferentially found at z > 2, may be connected to direct accretion of cold gas and to mergers (e.g., Di Matteo et al. 2005; Kazantzidis et al. 2005; Dubois et al. 2012), while less luminous and lower z AGNs seem to be driven by other processes, unrelated to the merging phenomenon in general (Cisternas et al. 2011; Schawinski et al. 2012; Treister et al. 2012; Kocevski et al. 2012; see Heckman & Best 2014 for a review).

A likely possibility is that the stellar mass losses normally produced during stellar evolution cyclically feed a central gas inflow (Norman & Scoville 1988), and then trigger QSO activity for isolated early-type galaxies (ETGs), that are often considered “dead.” Indeed, these losses represent a major source of mass for the ISM (e.g., Ciotti et al. 1991). Also, the

high metallicity of the observed outflows from low- and high-z galaxies provides evidence that the fuel source for these flows is highly processed gas, not “cold flows” accreted from outside (Cooksey et al. 2010; Fox 2011; Lehner et al. 2013). Numerical simulations describing the evolution of stellar mass losses in ETGs provide clear support for the conjecture that they trigger QSO activity at epochs closer than z ~ 2 (Ciotti et al. 2010). Galaxy–black hole coevolution under the effects of these internal processes has been generally termed “secular evolution” (e.g., Heckman & Best 2014).

In a number of previous works, Ciotti, Ostriker, and coworkers investigated the relationship between the secular evolution of the stellar population in ETGs, QSO activity, and the resulting feedback action on the host galaxy, since z ~ 2, after the MBH and the stellar population have terminated their major growing phase.<fnref n="5" /> With high-resolution hydrodynamical simulations in spherical symmetry, they studied in detail the mechanical and radiative feedback effects induced by accretion of stellar mass losses. The physics of feedback was modeled on the observed dominant processes: radiative output and broad absorption line (BAL) wind output. In fact, the emitted photons impart energy and momentum to the ISM via electron scattering, photoionization, scattering due to atomic resonance lines, and absorption by dust grains; in addition, accretion drives BAL winds that convey mass, momentum, and energy to the ISM surrounding the nucleus (as expected: Silk & Rees 1998; King 2003; Bieri et al. 2016 and as intensively

<footnote n="5">⁵ For a standard cosmology $z = 2$ corresponds to an age of the universe of 3.3 Gyr, an epoch at which massive ETGs are thought to have completed the bulk of their SF processes and MBH formation.</footnote>
observed: Reichard et al. 2003; Greene et al. 2011; Arav et al. 2013; Liu et al. 2013; Carniani et al. 2015; McElroy et al. 2015). The simulations consider the corresponding cooling and heating functions, including photoionization plus Compton scattering, and solved the radiative transport equations, also in the presence of dust; they allowed for mass and energy inputs from stellar winds, and Type Ia and Type II supernovae (SNIa and SNII); and they considered mechanical feedback due to nuclear wind and SF induced by accretion.

These simulations covered length scales from ~5 pc to ~200 kpc, and timescales from ~10² years (or less) to 10¹⁰ years; thus, all the relevant length and timescales were resolved (from the Bondi accretion radius to tens of optical effective radii), and the accretion rates, as well as the effect of AGN feedback on the gas over the whole galaxy, were self-consistently determined. This, together with the complex but exhaustively implemented input physics, was a specific and very important feature of the simulations, crucial to establish what is the exact accretion rate, and then the feedback effects on the final MBH mass, on the ISM, and on SF. In other numerical studies (e.g., Di Matteo et al. 2005, 2008; Booth & Schaye 2009; Choi et al. 2015; Sijacki et al. 2015) the mass accretion rate remains “unresolved,” and is set from recipes or algorithms, due to the lack of required spatial extent and resolution, or due to difficulties in resolving the gas mass distribution intrinsic to the numerical modeling (e.g., in SPH codes). For example, cosmological simulations suffer inevitable limits due to numerical resolution, and the mass accretion rate is given by recipes generally based on the Bondi rate (Bondi 1952), or ad hoc prescriptions based on it (with all the associated uncertainties; e.g., Hobbs et al. 2012; Curtis & Sijacki 2016; Korol et al. 2016). Also, the effect of the MBH on its surroundings has often been modeled by injecting thermal energy into the ISM again following recipes for its amount and distribution (typically with the goal of reproducing observed relations or general properties).

The simulations in spherical symmetry of Ciotti, Ostriker, and coworkers showed that in medium-high mass ETGs the resulting evolution is highly unsteady (e.g., Ciotti et al. 2010). At early times (starting from $z \sim 2$) major accretion episodes caused by cooling flows trigger AGN flaring, with duty cycles small enough to account for the small fraction of massive galaxies observed to be in the QSO phase, when the accretion luminosity approaches the Eddington luminosity. At low redshift, the majority of models are characterized by smooth, very sub-Eddington mass accretion rates. At the end of the evolution, the mass of the MBH is limited to the range of masses observed today, even though the mass lost by the stars is roughly two orders of magnitude larger than the MBH masses observed in local ETGs. Note that MBH heating alone has been shown to be insufficient by itself to avoid long-lasting and massive inflows toward the galactic center at early times, but when coupled with the SNIa's heating, it becomes very efficient in sustaining galaxy degassing and preventing large mass accumulation in the central regions. During the evolution roughly half of the mass lost by stars gets ejected in SNIa driven winds, and half falls to the center and ends accreted or in starbursts. This series of simulations in spherical symmetry was followed by an investigation in two dimensions, again applied to spherical galaxies (Novak et al. 2011, 2012; Gan et al. 2014). A two-dimensional (2D) implementation of the same feedback physics described above showed that MBH accretes

some of the infalling gas and expels a conical wind; and that the cool shells, forming at 0.1–1 kpc from the center, are Rayleigh–Taylor unstable to fragmentation, leading to a somewhat higher accretion rate, and less effective feedback.

In the current work we study the radiative and mechanical (due to AGN winds) feedback effects on the ISM by improving the treatment of our previous 2D works in two main respects. First, the galaxy models underlying the ISM evolution are more realistic and accurate than ever previously: they are axisymmetric, allow for various degrees of flattening, and include, in addition to the MBH, a generalized de Vaucouleurs stellar profile coupled with a Navarro–Frenk–White (NFW; Navarro et al. 1997) dark matter (DM) halo; the halo contributes less dark mass than the stellar mass within one effective radius, as observed recently for ETGs (e.g., Cappellari et al. 2015); the stellar kinematics is determined solving the Jeans equations for the total mass (MBH + stars + dark halo) and a chosen orbital distribution; all galaxy parameters are determined in order to keep the models on the main observed scaling laws. Second, we consider a secularly evolving stellar population input, i.e., a secularly decreasing stellar mass loss rate and SNIa explosion rate, distributed over the galaxy. With these two improvements implemented in our 2D simulations with well resolved length and timescales, we have addressed the following questions. What is the effect of feedback on the MBH mass? What is the increase of the MBH mass due to accretion at epochs more recent than $z \sim 2$? Is it plausible that the $M_{BH}-\sigma$ relation was already in place at $z \sim 2$? What is the effect of radiative feedback and BAL winds on the ISM? Are mass losses and AGN activity connected with residual SF episodes at late times (i.e., after most of the SF has completed)? Is feedback responsible for more or less SF? Do we correctly predict the properties of the circumgalactic medium?

The paper is organized as follows: Section 2 describes the galaxy models; Section 3 presents the numerical code and the hydrodynamical equations it solves, the inputs to them, and the implementation of the feedback and SF physics; Section 4 presents the results of the simulations; Section 5 summarizes the main conclusions.

## 2. GALAXY MODELS

We briefly summarize here the main characteristics of the galaxy models, with respect to their stellar population properties and evolution (Section 2.1), and their internal dynamics (Section 2.2).

### 2.1. Stellar Population

In ETGs the gas is lost by evolved stars mainly during the red giant, asymptotic giant branch, and planetary nebula phases. These losses originate ejecta that initially have the velocity of the parent star, then interact with the mass lost from other stars or with the hot ISM, and mix with it (Mathews 1990; Parrott & Bregman 2008; Bregman & Parrott 2009). Thus, stellar winds are heated to X-ray temperatures by thermalization of the kinetic energy of collisions between stellar ejecta, as will be presented in Section 2.2 below. Far infrared observations allow for measurements of the stellar mass loss rate for the whole galaxy ($\dot{M}_{*}$), giving an average rate in reasonable agreement with theoretical predictions (Athey et al. 2002). According to single burst stellar population synthesis models (Maraston 2005), the trend of $\dot{M}_{*}$ with time, for solar metal
abundance, after an age of ≳2 Gyr, can be approximated as

$$ \dot{M}_{*}(t) = 10^{-12} A \times M_{*} t_{12}^{-1.3} (M_{\odot} \text{yr}^{-1}), \quad (1) $$

where $M_*$ is the galactic stellar mass in solar masses at an age of 12 Gyr, $t_{12}$ is the age in units of 12 Gyr, and $A = 2.0$ or 3.3 for a Salpeter or Kroupa IMF (the latter is adopted here; see also Pellegrini 2012). The relation above agrees well with previous theoretical estimates (e.g., Mathews 1989; Ciotti et al. 1991).

Also SNIa explosions provide mass and heat to the ISM, and the total mass loss rate of a stellar population is $\dot{M}(t) = \dot{M}_{*}(t) + \dot{M}_{\text{SN}}(t)$, where the mass input due to SNIas is $\dot{M}_{\text{SN}}(t) = 1.4 M_{\odot} R_{\text{SN}}(t)$. Here $R_{\text{SN}}(t)$ (in yr$^{-1}$) is the evolution of the explosion rate with time, and each SNIa ejects 1.4 $M_{\odot}$. In models of SNIa explosions past a burst of SF (Greggio 2010), $R_{\text{SN}}(t)$ experiences a raising epoch during the first ≃1 Gyr, and then decreases slowly with a timescale of the order of 10 Gyr, down to the present day observed rate. A parameterization of the rate after the peak, in number of events per year, is

$$ R_{\text{SN}}(t) = 0.16(H_0/70)^2 \times 10^{-12} L_B t_{12}^{-s} (\text{yr}^{-1}), \quad (2) $$

where $H_0$ is the Hubble constant in units of km s$^{-1}$ Mpc$^{-1}$, $L_B$ is the present epoch B-band galaxy luminosity in $L_{B,\odot}$, and $s$ characterizes the secular evolution; when $t_{12} = 1$, Equation (2) gives the rate for local ETGs in recent measurements (e.g., Mannucci et al. 2005; Maoz et al. 2011). For the rate in Equation (2) and $H_0 = 70$ km s$^{-1}$ Mpc$^{-1}$, one obtains $\dot{M}_{\text{SN}}(12 \text{ Gyr}) = 2.2 \times 10^{-13} L_B M_{\odot} \text{yr}^{-1}$, that is almost ~100 times smaller than the “quiescent” stellar mass loss rate $\dot{M}_{*}(12 \text{ Gyr}) \approx 2 \times 10^{-11} L_B M_{\odot} \text{yr}^{-1}$ given above. Recent estimates of the slope $s$ agree with a value around $s \simeq 1$ (Sharon et al. 2010; Maoz et al. 2011).

The heating rate provided by SNIa explosions $L_{\text{SN}}(t)$ is the product of the kinetic energy injected by one event ($E_{\text{SN}} \approx 10^{51}$ erg) times the rate $R_{\text{SN}}(t)$, and times an efficiency factor. In the code we adopt an efficiency of 0.85, an assumption that is not unreasonable for hot diluted gas (see also Section 3.1). Then at most $L_{\text{SN}}(t) = E_{\text{SN}} R_{\text{SN}}(t)$, that is

$$ L_{\text{SN}}(t) = 5.1(H_0/70)^2 \times 10^{30} L_B t_{12}^{-s} \quad (\text{erg s}^{-1}). \quad (3) $$

Of course, another major source of ISM heating is provided by the central MBH, as will be discussed in Section 3 below.

## 2.2. Dynamical Structure

We consider here a subset of the large suite of models built for the study of gas flows in galaxies of various shapes and internal kinematics of Negri et al. (2014b). These are axisymmetric galaxy models, where, in addition to the central MBH (of initial mass $M_{\text{BH},0} = 10^{-3} M_*$), there is a stellar component characterized by different intrinsic flattenings, and a spherical DM halo. The stellar density is described by the ellipsoidal deprojection (Mellier & Mathez 1987) of the de Vaucouleurs (1948) law

$$ \rho_{*}(R, z) = \rho_{0}\zeta^{-0.855}\exp(-\zeta^{1/4}), \quad (4) $$

with

$$ \rho_0 = \frac{M_* b^{12}}{16\pi q R_{e0}^3 \Gamma(8.58)}, \quad \zeta = \frac{b^4}{R_{e0}} \sqrt{R^2 + \frac{z^2}{q^2}}, \quad (5) $$

where $(R, \phi, z)$ are the cylindrical coordinates, $b \simeq 7.67$, $R_{e0}$ is the projected half-mass radius (effective radius) when the galaxy is seen face-on,<fnref n="6" /> and the parameter $q \lesssim 1$ controls flattening, so that the minor axis is aligned with the z-axis. For the simulations we consider $q$ values of (1, 0.6, 0.3), corresponding to E0, E4, and E7 galaxies when seen edge-on, respectively. For the DM halo we adopt an untruncated NFW (Navarro et al. 1997) profile

$$ \rho_h(r) = \frac{\rho_{\text{crit}} \delta_c r_h}{r(1 + r/r_h)^2}, \quad (6) $$

where $\rho_{\text{crit}} = 3H^2/8\pi G$ is the critical density for closure, and

$$ \delta_c = \frac{200}{3} \frac{c^3}{\ln(1+c) - c/(1+c)}, \quad c \equiv \frac{r_{200}}{r_h} \quad (7) $$

and $r_{200}$ is the radius of a sphere of mean interior density of $200\rho_{\text{crit}}$. We refer to the DM mass enclosed within $r_{200}$ as the halo mass $M_h$.

All the relevant dynamical properties of the models were computed with a code built for this purpose (Posacki et al. 2013). Starting from an axisymmetric density distribution $\rho_*(R, z)$ produced by a two-integral phase-space distribution function, the code solves the Jeans equations in cylindrical coordinates, and computes the velocity fields of the stars, the total potential $\Phi_{\text{tot}}(R, z)$ due to all components (stars, dark halo, MBH), and the vertical and radial forces. The radial and vertical velocity dispersions are equal ($\sigma_R = \sigma_\phi \equiv \sigma$), and the only non-zero streaming motion is in the azimuthal direction ($\bar{v}_\phi$). To set the latter, we adopt the Satoh (1980) $k$-decomposition $\bar{v}_\phi^2 = k^2(\bar{v}_\phi^2 - \sigma^2)$, from which the azimuthal velocity dispersion is recovered as $\sigma_\phi^2 \equiv \bar{v}_\phi^2 - \bar{v}_\phi^2 = \sigma^2 + (1 - k^2)(\bar{v}_\phi^2 - \sigma^2)$, where $0 \le k \le 1$. For $k = 1$ the galaxy is an isotropic rotator, while for $k = 0$ no net rotation is present, and all the flattening is due to $\sigma_\phi$. In general, $k$ can be a function of $(R, z)$, and more complicated (realistic) velocity fields can be realized (Ciotti & Pellegrini 1996; Negri et al. 2014a, 2014b). In any case, $k$ is bounded from above by the function $k_{\text{max}}^2(R, z) = \bar{v}_\phi^2 / (\bar{v}_\phi^2 - \sigma^2)$. In this work we restrict ourselves to the $k \simeq 0$ case, and then in practice only stellar random motions are thermalized (see also Section 4 for more discussion on the role and then the adopted values of the $k$ parameter).

The only free parameter of the stellar distribution is $\sigma_{e8}$, the aperture luminosity-weighted velocity dispersion within $R_e/8$, from which the galaxy luminosity is recovered from the Faber–Jackson relation, and the size $R_e$ from the size–luminosity relation (Desroches et al. 2007). From the stellar mass-to-light ratio, fixed at that of a 12 Gyr old stellar population with a Kroupa initial mass function (IMF), the stellar mass $M_*$ is derived (see Posacki et al. 2013 for more details). By assigning a given $\sigma_{e8}$, we build a spherical galaxy, that we call the “progenitor.” The free parameters of the DM halo are determined by the need to reproduce the assumed $\sigma_{e8}$, by fixing $M_h/M_* \simeq 20$ (Behroozi et al. 2013), and imposing that the DM fraction $f_{\text{DM}}$ within a sphere of radius $R_e$ keeps well below unity (Cappellari et al. 2015). These constraints produce $r_h \simeq 2R_e$, $22 \lesssim c \lesssim 41$, and $f_{\text{DM}} \simeq 0.6$ for the spherical progenitors. The flattened descendants of each progenitor have the same circularized $R_e$ as the spherical progenitor, when seen

<footnote n="6">⁶ For an edge-on view, the circularized effective radius is $R_e = R_{e0} \sqrt{q}$.</footnote>
<table><caption>**Table 1**
Structural Parameters of the Galaxy Models</caption><thead><tr><th>Name<br>(1)</th><th>L<sub>B</sub><br>(10<sup>11</sup>L<sub>B,☉</sub>)<br>(2)</th><th>R<sub>e</sub><br>(kpc)<br>(3)</th><th>M<sub>*</sub><br>(10<sup>11</sup>M<sub>☉</sub>)<br>(4)</th><th>M<sub>h</sub><br>(10<sup>11</sup>M<sub>☉</sub>)<br>(5)</th><th>M<sub>BH,0</sub><br>(10<sup>8</sup>M<sub>☉</sub>)<br>(6)</th><th>σ<sub>e8</sub><br>(km s<sup>-1</sup>)<br>(7)</th><th>f<sub>DM</sub><br>(8)</th><th>c<br>(9)</th></tr></thead><tbody><tr><td>E4<sup>180</sup></td><td>0.18</td><td>3.26</td><td>0.81</td><td>16.20</td><td>0.81</td><td>160</td><td>0.62</td><td>41</td></tr><tr><td>E7<sup>180</sup></td><td>0.18</td><td>3.26</td><td>0.81</td><td>16.20</td><td>0.81</td><td>137</td><td>0.73</td><td>41</td></tr><tr><td>E4<sup>210</sup></td><td>0.32</td><td>4.57</td><td>1.54</td><td>30.80</td><td>1.54</td><td>187</td><td>0.62</td><td>35</td></tr><tr><td>E7<sup>210</sup></td><td>0.32</td><td>4.57</td><td>1.54</td><td>30.80</td><td>1.54</td><td>155</td><td>0.66</td><td>35</td></tr><tr><td>E4<sup>250</sup></td><td>0.65</td><td>7.04</td><td>3.35</td><td>67.00</td><td>3.35</td><td>223</td><td>0.63</td><td>28</td></tr><tr><td>E7<sup>250</sup></td><td>0.65</td><td>7.04</td><td>3.35</td><td>67.00</td><td>3.35</td><td>184</td><td>0.67</td><td>28</td></tr><tr><td>E4<sup>300</sup></td><td>1.38</td><td>11.8</td><td>7.80</td><td>160.00</td><td>7.80</td><td>267</td><td>0.66</td><td>22</td></tr><tr><td>E7<sup>300</sup></td><td>1.38</td><td>11.8</td><td>7.80</td><td>160.00</td><td>7.80</td><td>221</td><td>0.68</td><td>22</td></tr></tbody></table>

**Note.** (1) Model name: the letters and numbers identify the shape (E4 or E7), the superscript gives the $\sigma_{e8}$ of the progenitor. (2) Luminosity in the B-band. (3) Effective radius (edge-on view). (4) Total stellar mass. (5) Total DM mass. (6) Black hole initial mass. (7) Luminosity-weighted stellar velocity dispersion within a circular aperture of radius $R_e/8$, for an edge-on view. (8) Ratio of the DM mass to the total mass enclosed within a sphere of radius $R_e$. (9) Concentration parameter of the NFW halo. See Section 2 for more details.

edge-on, thus they expand with decreasing $q$ ($\rho_* \propto \sqrt{q}$); their $\sigma_{e8}$, as a consequence, also decreases with respect to that of the progenitor (while $M_*$ and $L_B$ remain the same). In the flattening procedure the DM halo is maintained fixed to that of the progenitor (see Posacki et al. 2013 for a more comprehensive model description$^7$).

We consider here four spherical progenitors, with $\sigma_{e8}$ values of 180, 210, 250, and 300 km s$^{-1}$, and corresponding stellar masses $M_*$ of 0.81, 1.54, 3.35, and $7.80 \times 10^{11} M_⊙$. The simulations are run for eight flattened descendants: for each progenitor, they have intrinsic flattenings with $q = 0.6$ (E4 shape) and $q = 0.3$ (E7 shape). Table 1 lists all the relevant parameters characterizing these eight galaxy models.

### 3. FEEDBACK AND SF IN THE HYDRODYNAMICAL EQUATIONS

We describe here the hydrodynamical equations that are solved to evolve the gas flow (Section 3.1), and the input physics for them, with particular regard to the implementation of the SF process (Section 3.2) and feedback effects (Sections 3.3 and 3.4).

#### 3.1. The Hydrodynamical Equations

The hydrodynamical equations that are numerically integrated in spherical coordinates, by a significantly updated version of the Novak et al. (2011) code, are the following:

$$ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) = \dot{\rho}_{\text{Ia}} + \dot{\rho}_{*} + \dot{\rho}_{\text{II}} - \dot{\rho}_{\text{SF}} + \dot{\rho}_{\text{w}}, \quad (8) $$

$$ \begin{aligned} \rho \frac{\partial \mathbf{u}}{\partial t} + \rho (\mathbf{u} \cdot \nabla) \mathbf{u} = & - \nabla p - \rho \nabla \Phi_{\text{tot}} - \nabla p_{\text{rad}} \\ & + (\dot{\rho}_{\text{Ia}} + \dot{\rho}_{*} + \dot{\rho}_{\text{II}}) \\ & \times (\mathbf{v}_{*} - \mathbf{u}) + \dot{\rho}_{\text{w}} (\mathbf{v}_{\text{w}} - \mathbf{u}), \end{aligned} \quad (9) $$

$$ \begin{aligned} \frac{\partial E}{\partial t} + \nabla \cdot (E \mathbf{u}) = & -p \nabla \cdot \mathbf{u} + H - C + \dot{E} \\ & - \dot{E}_{\text{SF}} + \frac{\dot{\rho}_{\text{w}}}{2} \|\mathbf{v}_{\text{w}} - \mathbf{u}\|^2 \end{aligned} \quad (10) $$

where $\rho$, $\mathbf{u}$, $E$, $p$ are respectively the mass density, velocity, internal energy density, and pressure of the gas; $\mathbf{v}_*$ is the streaming velocity of the stellar component; $\mathbf{v}_\text{w}$ is the AGN wind velocity. $H$ is the volumetric heating rate due to radiative feedback, and $C$ is the bolometric cooling rate per unit volume (see Section 3.4). The ISM is a fully ionized monoatomic gas, with $p = (\gamma - 1)E$, $\gamma = 5/3$, and solar composition ($\mu = 0.62$). The mass density rates $\dot{\rho}_{\text{Ia}}$ and $\dot{\rho}_*$ describe the mass injection from the old stellar population, i.e., from SNIa explosions and from normal stars (Section 2.1); $\dot{\rho}_{\text{SF}}$ and $\dot{\rho}_{\text{II}}$ describe the mass sink due to SF, and the mass input due to the SNII produced by SF (see Section 3.2); $\dot{E}_{\text{SF}}$ is the energy sink due to SF (Section 3.2). $\dot{E}$ is the total energy injection rate due to the old and new stellar populations, produced by the thermalization of the kinetic energy of SNIa and SNII explosions ($\dot{E}_{\text{Ia}}$ and $\dot{E}_{\text{II}}$ respectively), and by the thermalization of the relative motions between stars and the ISM (existing at the moment of injection). $\dot{E}$ is then given by

$$ \dot{E} = \dot{E}_{\text{Ia}} + \dot{E}_{\text{II}} + \frac{\dot{\rho}_{\text{Ia}} + \dot{\rho}_{*} + \dot{\rho}_{\text{II}}}{2} [\|\mathbf{v}_{*} - \mathbf{u}\|^2 + \text{Tr}(\boldsymbol{\sigma}^2)], \quad (11) $$

where $\dot{E}_{\text{Ia}} = \dot{\rho}_{\text{Ia}} \vartheta_{\text{SNIa}} E_{\text{SN}} / (1.4 M_\odot)$, with $E_{\text{SN}} = 10^{51}$ erg and $1.4 M_\odot$ being respectively the kinetic energy and ejected mass of one SNIa event, and $\vartheta_{\text{SNIa}}$ is the thermalization efficiency, for which we adopt the value of 0.85, as a plausible one for a low density and hot medium (see Mathews 1989; Tang & Wang 2005); $\dot{E}_{\text{II}}$ is the SNII energy injection rate, and is calculated considering the overlap in time of subsequent SF episodes (as for $\dot{\rho}_{\text{II}}$ below; see Section 3.2, and Negri et al. 2015). $\dot{E}_{\text{II}}$ also depends on the thermalization efficiency, that we take equal to $\vartheta_{\text{SNIa}}$. Finally, $\text{Tr}(\boldsymbol{\sigma}^2)$ is the trace of the stellar velocity dispersion tensor.

<footnote>$^7$ The models considered here belong to the class of “edge-on built” descendants of Posacki et al. (2013), while the latter authors also built the class of “face-on built” descendants, that have the same $R_e$ as the spherical progenitor when seen face-on.</footnote>
The other symbols in Equations (8)–(10) describe the mass, energy, and momentum source terms due to radiative and mechanical feedback, and are described in Sections 3.3 and 3.4 below. In particular, the terms $\dot{\rho}_w$ and $\mathbf{v}_w$ are given in Equations (25)–(28). $\nabla p_{\text{rad}} = (\nabla p_{\text{rad}})_{\text{photo}} + (\nabla p_{\text{rad}})_{\text{es}}$ is the total radiative pressure gradient for the radiation coming from the accreting MBH. Most of the input physics concerning feedback and SF is the same as in Ciotti et al. (2010), with a few modifications described in detail below.

## 3.2. Star Formation

SF and consequent SNII production are treated as in Novak et al. (2011) and Negri et al. (2015). SF is implemented by subtracting gas from the grid, and this mass sink corresponds to an SF rate per unit volume given by

$$ \dot{p}_{\text{SF}} = \frac{\eta_{\text{SF}} \rho}{\tau_{\text{SF}}}, \qquad (12) $$

where $\rho$ is the local gas density, $\eta_{\text{SF}}$ is the SF efficiency, for which we adopt the value of 0.1, and

$$ \tau_{\text{SF}} = \max(\tau_{\text{cool}}, \tau_{\text{dyn}}), \qquad (13) $$

where

$$ \tau_{\text{cool}} = \frac{E}{C}, \quad \tau_{\text{dyn}} = \min(\tau_{\text{jeans}}, \tau_{\text{rot}}), \qquad (14) $$

with

$$ \tau_{\text{jeans}} = \sqrt{\frac{3\pi}{32G\rho}}, \quad \tau_{\text{rot}} = \frac{2\pi r}{v_c(r)}. \qquad (15) $$

Here $G$ is the Newtonian gravitational constant, $r$ is the distance from the galaxy center, $\tau_{\text{rot}}$ is an estimate of the radial epicyclic period, and $v_c(R)$ is the galaxy circular velocity in the equatorial plane. The energy and momentum sinks associated with SF are then

$$ \dot{E}_{\text{SF}} = \frac{\eta_{\text{SF}} E}{\tau_{\text{SF}}}, \quad \dot{\mathbf{m}}_{\text{SF}} = \frac{\eta_{\text{SF}} \mathbf{m}}{\tau_{\text{SF}}} = \dot{p}_{\text{SF}} \mathbf{u}, \qquad (16) $$

where $E$ and $\mathbf{m}$ are the internal energy and momentum density of the ISM.

SF removes mass, momentum, and energy from the grid, but also injects new mass and energy from SNII explosions. For each SF episode, assuming that the new stars form with a Salpeter IMF, the mass returned in SNII events is 20% of the new star mass in that episode; the SNII mass source term $\dot{\rho}_{\text{II}}$ at each time $t$ comes from considering that a given SF episode generates SNIIs that inject mass (at a rate exponentially declining on a timescale $\tau_{\text{II}} = 2 \times 10^7$ years), and that during the evolution of that episode other episodes may take place, forming other SNIIs that in turn eject mass into the ISM. The same considerations are taken into account to compute the SNII energy injection rate $\dot{E}_{\text{II}}$ (see Negri et al. 2015 for more details on how $\dot{\rho}_{\text{II}}$ and $\dot{E}_{\text{II}}$ are computed).

## 3.3. Mechanical Feedback

Following the general, self-consistent treatment in Ostriker et al. (2010), the basic quantities involved in the

implementation of mechanical feedback are written as

$$ \dot{M}_{\text{BH}} = \frac{\dot{M}_{\text{in}}}{1 + \eta}, \qquad (17) $$

$$ \dot{M}_{\text{out}} = \eta \dot{M}_{\text{BH}}, \qquad (18) $$

$$ L_w = \epsilon_w \dot{M}_{\text{BH}} c^2, \qquad (19) $$

$$ \dot{p}_w = \dot{M}_{\text{out}} v_w, \qquad (20) $$

where $\dot{M}_{\text{BH}}$ is the mass accretion rate on the MBH, $\dot{M}_{\text{out}}$ is the mass outflow rate in the conical wind, $\dot{M}_{\text{in}}$ is the mass inflow rate at the first active radial grid, $\epsilon_w$ is the efficiency of generating mechanical energy with an AGN wind, $v_w = \|\mathbf{v}_w\|$ is the modulus of the AGN wind velocity, and $\eta \equiv 2\epsilon_w c^2 / v_w^2$. With the adoption of this simplified scheme, the sub-grid physics near the MBH (described in Ciotti & Ostriker 2012) is not activated here. Then, in the present simulations, all gas that flows in eventually either flows onto the MBH or back into the simulation grid as a conical wind, without passing through a circumnuclear star-forming disk.

Equations (17)–(20) above are adopted by most authors to treat AGN feedback as a process comprising both infall and outflow; however, typically $\eta = 0$ is adopted, implicitly assuming $v_w \to \infty$, so that $\dot{M}_{\text{out}}$ and $\dot{p}_w$ are neglected, and $L_w$ and $\dot{M}_{\text{BH}}$ may be overestimated. For example, if we adopt $\epsilon_w = 5 \times 10^{-3}$ as many authors have done (e.g., Springel et al. 2005; Johansson et al. 2009), and $v_w = 10^4 \text{ km s}^{-1}$ (e.g., Moe et al. 2009), then $\eta = 9$ and all the normally neglected effects may in fact be dominant: the bulk of the inflowing mass ($\dot{M}_{\text{in}}$) may be ejected in a broad-line disk wind, and the effects of the mass ($\dot{M}_{\text{out}}$) and momentum ($\dot{p}_w$) input deposited in the ambient gas may dominate over the energy input ($L_w$), which may be largely radiated away. Thus, it is important to consider consistently the effects of including mass, energy, and momentum conservation when $\eta > 0$; our treatment is consistent with observations of BAL quasars and the radiative mechanisms that drive them (Ciotti et al. 2010).

In general, $\eta$ is fixed for given $\epsilon_w$ and $v_w$. The wind efficiency $\epsilon_w$ is not known very well, neither from observations nor from detailed simulations. The best estimates might be in the range $10^{-3} > \epsilon_w > 3 \times 10^{-4}$ (Proga et al. 2000; Proga & Kallman 2004; Krongold et al. 2007; Kurosawa et al. 2009). Here we consider models where $\epsilon_w$ and $v_w$ are allowed to depend on the mass accretion rate, and are described by the following laws (see Novak et al. 2011 for a discussion on these assumptions):

$$ \epsilon_w = \frac{\epsilon_{w0} A_w \dot{m}}{1 + A_w \dot{m}}, \quad v_w = \frac{v_{w0} A_w \dot{m}}{1 + A_w \dot{m}}, \qquad (21) $$

where $A_w = 1000$, $\dot{m}$ is the dimensionless mass accretion rate, i.e., normalized to the Eddington mass accretion rate (see below), and the two functions $\epsilon_w$ and $v_w$ saturate to $\epsilon_{w0}$ and $v_{w0}$, respectively, for large accretion rates. We choose the constant values $\epsilon_{w0} = 10^{-4}$ and $v_{w0} = 10^4 \text{ km s}^{-1}$ (as observed for the outflow velocity in UV absorption lines of BAL AGNs; e.g., Reichard et al. 2003; Gibson et al. 2009 in ionized emission of high-redshift quasars, e.g., Liu et al. 2013; Zakamska et al. 2016; in nearby Seyfert galaxies, e.g., Fischer et al. 2013; and in molecular outflows, Feruglio et al. 2015; Tombesi et al. 2015). Independently of the mechanical feedback model, the radiative luminosity of the AGN is given
by

$$L_{BH} = \epsilon_{EM} \dot{M}_{BH} c^2, \qquad (22)$$

where the electromagnetic efficiency $\epsilon_{EM}$ is given by the advection dominated accretion flow inspired formula (Narayan & Yi 1995), to reproduce even the very low nuclear luminosities typically observed for local MBHs (e.g., Pellegrini 2005):

$$\epsilon_{EM} = \frac{\epsilon_0 A_{EM} \dot{m}}{1 + A_{EM} \dot{m}}, \qquad (23)$$

and $A_{EM} = 100$ and $\epsilon_0 = 0.125$. The dimensionless mass accretion rate is

$$\dot{m} \equiv \frac{\dot{M}_{BH}}{\dot{M}_{Edd}} = \frac{\epsilon_0 \dot{M}_{BH} c^2}{L_{Edd}}, \qquad (24)$$

where $L_{Edd}$ is the Eddington luminosity. With the settings above, $\eta$ has a minimum value of 0.18 at the highest accretion rates, and increases for decreasing $\dot{m}$ (for example, $\eta = 0.36$ for $\dot{m} = 10^{-3}$).

Our previous 1D simulations used a prescription based on pressure balance between the outgoing wind and the ambient gas, to compute how the mass, energy, and momentum are radially distributed. The present 2D simulations instead inject the desired mass, energy, and momentum into the innermost radial cells, and self-consistently compute the radial transport of these quantities. In the 2D simulations, the total mass, energy, and momentum injected into the ISM by the AGN wind are calculated as described by Equations (18)–(20) for a given MBH accretion rate $\dot{M}_{BH}$; thus, a specification of the angular dependence of the properties of the AGN conical wind ($\dot{\rho}_w$, $\dot{m}_w$ and $\dot{E}_w$) is required. For this dependence we adopt

$$\dot{\rho}_w = \dot{M}_{\text{out}} \frac{\delta(r - r_0)}{r^2 \sin \theta} f(\theta), \qquad (25)$$

where $\delta$ is the Dirac delta-function, $r_0$ is the first gridpoint, $\theta$ is the angle from the z-axis, measured clockwise, and

$$f(\theta) = \frac{(n + 1) \sin \theta |\cos \theta|^n}{4\pi}. \qquad (26)$$

The anisotropic momentum source in Equation (9) is given by

$$\dot{m}_w = \dot{\rho}_w v_w e_r = \dot{\rho}_w \frac{\delta(r - r_0)}{r^2 \sin \theta} f(\theta) e_r \qquad (27)$$

where the second equality follows from using Equations (21) and (26). Finally, the energy injection due to the AGN wind, from Equations (20) and (26), is

$$\dot{E}_w = \frac{1}{2} \dot{\rho}_w v_w^2 = L_w \frac{\delta(r - r_0)}{r^2 \sin \theta} f(\theta). \qquad (28)$$

Note that $v_w$ is independent of $\theta$, and the anisotropy in the momentum and energy injection is due to the anisotropic mass injection (Equation (26)). The wind injects mass, momentum, and energy within two symmetric cones above and below the equatorial plane. Following Novak et al. (2011), we adopt $n = 2$, that corresponds to $\approx 40^{\circ}$ for the half-opening angle of each of the two cones enclosing half of the mass, momentum, and energy injected by the wind in each of the two half-spaces. The two “effective” cones

correspond to $\approx 1/5$ of the total solid angle.<fnref n="8" /> Note however that, from Equation (26), the conical wind acts over the whole spherical solid angle (reducing to zero on the equatorial plane only); so that the figures given are just indicative, and the region affected by the wind is larger. These figures agree with observations of the fraction of obscured and unobscured AGNs under the assumption that the two populations are made up of a single population of objects that differ only in viewing angle (e.g., Tombesi et al. 2010; Liu et al. 2015; Bae & Woo 2016).

### 3.4. Radiative Feedback: Heating and Cooling

Radiative heating and cooling are computed by using the formulae in Sazonov et al. (2005), which describe the net heating/cooling rate per unit volume of a plasma in photoionization equilibrium with a radiation field characterized by the average quasar spectral energy distribution (Sazonov et al. 2005, 2008), whose associated spectral temperature is $T_c \approx 2$ keV. In particular, Compton heating and cooling, bremsstrahlung losses, line and continuum heating, and cooling are taken into account. The net gas energy change rate per unit volume for $T \gtrsim 10^4$ K is given by

$$H - C \equiv n^2 (S_1 + S_2 + S_3), \qquad (29)$$

where $n$ is the hydrogen number density, and positive and negative terms are grouped together in the heating ($H$) and cooling ($C$) functions (all quantities are expressed in cgs system). The bremsstrahlung losses are given by

$$S_1 = -3.8 \times 10^{-27} \sqrt{T}. \qquad (30)$$

The Compton heating and cooling is given by

$$S_2 = 4.1 \times 10^{-35} (T_c - T) \xi, \qquad (31)$$

and $\xi$ is the ionization parameter. The sum of photoionization heating, line and recombination continuum cooling is

$$S_3 = 10^{-23} \frac{a + b (\xi/\xi_0)^c Z}{1 + (\xi/\xi_0)^c Z_\odot}, \qquad (32)$$

where the almost perfect linear dependence on metallicity is explicit, and

$$a = - \frac{18}{e^{25(\log T - 4.35)^2}} - \frac{80}{e^{5.5(\log T - 5.2)^2}} - \frac{17}{e^{3.6(\log T - 6.5)^2}}, \qquad (33)$$

$$b = 1.7 \times 10^4 T^{-0.7}, \quad c = 1.1 - \frac{1.1}{e^{T/1.8 \times 10^5}} + \frac{4 \times 10^{15}}{T^4}, \qquad (34)$$

$$\xi_0 = \left( \frac{1.5}{T^{0.5}} + \frac{1.5 \times 10^{12}}{T^{2.5}} \right)^{-1} + \frac{4 \times 10^{10}}{T^2} \left[ 1 + \frac{80}{e^{(T-10^4)/1.5 \times 10^3}} \right]. \qquad (35)$$

<footnote n="8">⁸ The half-opening angle $\theta_q$ of each cone containing a fraction $q$ of the emitted mass has $\cos \theta_q = (1 - q)^{1/(n+1)}$, and the fraction of the hemispheric solid angle $2\pi$ covered by the top of each cone is $1 - \cos \theta_q$.</footnote>
Gas temperatures are bounded from below by the adopted atomic cooling curve, that has an exponential cutoff below $10^4$ K.

In this work we do not consider the effects of dust, i.e., the radiation momentum transfer associated with dust absorption (cf. Debuhr et al. 2011; Ciotti & Ostriker 2012; Hensley et al. 2014; Bieri et al. 2016). Moreover, since the effects of dust absorption and reprocessing are not considered here, the radiation pressure due to the reprocessing of light emitted by stars formed during the evolution is also omitted. For the radiation pressure we consider electron scattering and the force exerted by absorption of AGN photons by atomic lines (see Equation (55) in Ciotti & Ostriker 2007). The integration scheme used here for the radiation transport is not the full version in Novak et al. (2012, their Section 3), but the simplified optically thin case, in order to speed up the simulations. However, we stress that the integration is done radius by radius, and along each radius the ISM density and temperature vary as prescribed by the hydrodynamics, so that heating and cooling are not spherically symmetric, but depend on the 2D ISM properties at each time step. In the case of a spherically symmetric ISM, the equation to be integrated would be

$$ \frac{dL_{\text{BH,photo}}^{\text{eff}}(r)}{dr} = -4\pi r^2 H, \qquad (36) $$

where $L_{\text{BH,photo}}^{\text{eff}}(r)$ is the effective accretion luminosity at $r$, and the equation is solved with central boundary condition $L_{\text{BH,photo}}^{\text{eff}}(r=0) = L_{\text{BH}}$ given by Equation (22). The force per unit mass due to photoionization + Compton opacity would be expressed as

$$ (\nabla p_{\text{rad}})_{\text{photo}} = - \frac{\rho \kappa_{\text{photo}} L_{\text{BH,photo}}^{\text{eff}}(r)}{c 4\pi r^2} e_r, \qquad (37) $$

where

$$ \begin{aligned} \kappa_{\text{photo}} &= -\frac{1}{\rho(r)L_{\text{BH,photo}}^{\text{eff}}(r)} \frac{dL_{\text{BH,photo}}^{\text{eff}}(r)}{dr} \\ &= \frac{4\pi r^2 H(r)}{\rho(r)L_{\text{BH,photo}}^{\text{eff}}(r)}. \end{aligned} \qquad (38) $$

The radiation pressure due to electron scattering would be

$$ (\nabla p_{\text{rad}})_{\text{es}} = - \frac{\rho \kappa_{\text{es}} L_{\text{BH}}}{c 4\pi r^2} e_r, \qquad (39) $$

where $\kappa_{\text{es}} = 0.35 \text{ cm}^2 \text{ g}^{-1}$. Finally, Equations (31) and (32) depend on the ionization parameter, that would be given by

$$ \xi = \frac{L_{\text{BH,photo}}^{\text{eff}}(r)}{n(r)r^2}. \qquad (40) $$

In summary, in the present simulations we maintained the treatment of radiation feedback at a minimum level of complexity, restricting ourselves to the optically thin case. We plan to run more realistic simulations with the full equations of radiative transport activated, including the radiation reprocessing (from X-ray, to UV, to optical, to IR), and non negligible galactic rotation, with the consequent formation of a cold and optically thick equatorial disk (e.g., see Negri et al. 2014a, 2014b).

## 4. THE SIMULATIONS

We employed our modified version of the parallel ZEUS code (Hayes et al. 2006), in a 2D axisymmetric configuration, with a radially logarithmic grid in spherical coordinates $(r, \theta)$ of $128 \times 32$ meshpoints, spanning from 2.5 pc to 250 kpc. Reflecting boundary conditions are set along the z-axis, while at the outer edge of the computational domain the fluid is free to flow out.

Negri et al. (2014a, 2014b) showed that in a galaxy with substantial ordered rotation (without AGN feedback) the gaseous halo is almost co-rotating with the stars, and angular momentum conservation leads to the formation of a cold, rotationally supported, star-forming equatorial disk of kpc-scale (Negri et al. 2015), thus preventing any substantial accretion on the central MBH. These massive disks are expected to be gravitationally unstable, fragment, and consequently transport material to the galaxy center (Bertin & Lodato 2001; Hopkins & Quataert 2011). At the present stage we do not account for these processes, and we restrict ourselves to the low-rotation case. However, in axisymmetric systems some rotation is numerically needed to prevent gas from unphysically sticking on the z-axis. Different recipes have been figured out to solve this problem (Novak et al. 2011; Li et al. 2013); here, we take advantage of the Jeans solver that allows us to tune the Satoh k-parameter (Section 2.2). In detail, we determine $k(r, \theta)$ in order to have negligible but non-zero ordered rotation over the main galaxy body, with angular momentum of the stars never exceeding $J_0$, the specific angular momentum of the circular orbit at the first radial gridpoint, in the gravitational field of the MBH. In this way the gas is allowed to enter the first gridpoint, and so be accreted on the MBH. In practice, we consider a low rotation regime, where the centrifugal barrier prevents gas from sticking onto the z-axis, but allows for accretion down to the innermost radial gridpoint; such a rotation field is built by defining the Satoh parameter as follows:

$$ k(r, \theta) = \frac{\eta_{\text{rot}} J_0}{\sqrt{J_0^2 + R^2 v_{\text{IS}}^2}}; \qquad v_{\text{IS}}^2 = \overline{v_{\varphi}^2} - \sigma^2, \qquad (41) $$

where $R = r \sin\theta$ is the distance from the z-axis, $v_{\text{IS}}(r, \theta)$ is the ordered stellar velocity if the galaxy model were an isotropic rotator, and we fix $\eta_{\text{rot}} = 0.9$.

We computed various X-ray properties of the gas flows, with the X-ray luminosity in the 0.3–8 keV band and the X-ray-emission-weighted temperature respectively defined as

$$ L_X = \int \varepsilon_X dV, \qquad T_X = \frac{1}{L_X} \int T \varepsilon_X dV, \qquad (42) $$

where $\varepsilon_X$ is the thermal emissivity in the 0.3–8 keV band of a hot, collisionally ionized plasma (see Negri et al. 2014a for more details), and the integrals are computed over the volume of interest. We also calculated the X-ray surface brightness maps and maps of projected temperatures (as detailed in Pellegrini et al. 2012). We also calculated the following general properties of the models: the duty cycle ($D$), defined as the percentage of time spent at a $l \equiv L_{\text{BH}}/L_{\text{Edd}} > 0.05$; the time at which half of the MBH radiation energy has been emitted ($t_L$); the time at which half of the new stellar mass has been created ($t_M$); and the half-mass radius of the new stars formed until the end of the simulation ($r_M$).
<table><caption>**Table 2**
Simulations Results</caption><thead><tr><th>Name<br>(1)</th><th>ΔM<sub>BH</sub><br>(10<sup>8</sup> M<sub>☉</sub>)<br>(2)</th><th>ΔM<sub>BH</sub><br>(M<sub>BH,0</sub>)<br>(3)</th><th>l<sub>0.5</sub><br>(4)</th><th>D<br>(%)<br>(5)</th><th>t<sub>L</sub><br>(Gyr)<br>(6)</th><th>ΔM<sub>*</sub><br>(10<sup>8</sup> M<sub>☉</sub>)<br>(7)</th><th>SFR<br>(M<sub>☉</sub> yr<sup>-1</sup>)<br>(8)</th><th>t<sub>M</sub><br>(Gyr)<br>(9)</th><th>r<sub>M</sub><br>(kpc)<br>(10)</th><th>L<sub>X</sub><br>(10<sup>40</sup> erg s<sup>-1</sup>)<br>(11)</th><th>T<sub>X</sub><br>(keV)<br>(12)</th></tr></thead><tbody><tr><td colspan="13"><strong>Rad+mech</strong><br>(FF)</td></tr><tr><td>E4<sup>180</sup></td><td>0.15</td><td>0.18</td><td>0.12</td><td>5.11</td><td>4.56</td><td>14.51</td><td>0.04</td><td>5.38</td><td>9.09</td><td>0.02</td><td>0.62</td></tr><tr><td>E7<sup>180</sup></td><td>0.06</td><td>0.07</td><td>0.12</td><td>4.73</td><td>3.27</td><td>7.21</td><td>0.01</td><td>4.66</td><td>10.23</td><td>0.007</td><td>0.61</td></tr><tr><td>E4<sup>210</sup></td><td>1.54</td><td>1.01</td><td>0.03</td><td>1.60</td><td>5.04</td><td>77.60</td><td>0.28</td><td>5.25</td><td>1.68</td><td>1.31</td><td>0.47</td></tr><tr><td>E7<sup>210</sup></td><td>1.06</td><td>0.69</td><td>0.03</td><td>1.17</td><td>4.83</td><td>64.82</td><td>0.13</td><td>5.22</td><td>1.68</td><td>0.57</td><td>0.47</td></tr><tr><td>E4<sup>250</sup></td><td>6.31</td><td>1.90</td><td>0.04</td><td>2.94</td><td>5.08</td><td>170.16</td><td>0.31</td><td>4.97</td><td>2.44</td><td>3.72</td><td>0.77</td></tr><tr><td>E7<sup>250</sup></td><td>4.91</td><td>1.48</td><td>0.04</td><td>2.18</td><td>5.24</td><td>155.71</td><td>0.44</td><td>5.24</td><td>2.29</td><td>2.88</td><td>0.63</td></tr><tr><td>E4<sup>300</sup></td><td>23.19</td><td>3.00</td><td>0.06</td><td>3.43</td><td>5.92</td><td>342.59</td><td>1.02</td><td>5.04</td><td>4.22</td><td>23.57</td><td>1.02</td></tr><tr><td>E7<sup>300</sup></td><td>20.33</td><td>2.63</td><td>0.05</td><td>3.03</td><td>5.91</td><td>309.85</td><td>2.95</td><td>5.16</td><td>4.76</td><td>47.62</td><td>0.92</td></tr><tr><td colspan="13"><strong>Mech</strong><br>(MF)</td></tr><tr><td>E4<sup>180</sup></td><td>0.16</td><td>0.20</td><td>...</td><td>...</td><td>...</td><td>14.05</td><td>0.02</td><td>5.19</td><td>11.51</td><td>0.006</td><td>0.64</td></tr><tr><td>E7<sup>180</sup></td><td>0.06</td><td>0.07</td><td>...</td><td>...</td><td>...</td><td>5.63</td><td>0.01</td><td>4.17</td><td>9.09</td><td>0.004</td><td>0.60</td></tr><tr><td>E4<sup>210</sup></td><td>1.80</td><td>1.18</td><td>...</td><td>...</td><td>...</td><td>68.44</td><td>0.33</td><td>5.83</td><td>2.93</td><td>1.57</td><td>0.48</td></tr><tr><td>E7<sup>210</sup></td><td>1.37</td><td>0.90</td><td>...</td><td>...</td><td>...</td><td>61.55</td><td>0.26</td><td>5.42</td><td>2.59</td><td>1.08</td><td>0.45</td></tr><tr><td>E4<sup>250</sup></td><td>7.29</td><td>2.20</td><td>...</td><td>...</td><td>...</td><td>152.91</td><td>2.30</td><td>5.68</td><td>3.97</td><td>12.39</td><td>1.06</td></tr><tr><td>E7<sup>250</sup></td><td>5.74</td><td>1.73</td><td>...</td><td>...</td><td>...</td><td>136.96</td><td>0.88</td><td>5.91</td><td>3.74</td><td>4.89</td><td>0.67</td></tr><tr><td>E4<sup>300</sup></td><td>22.98</td><td>2.98</td><td>...</td><td>...</td><td>...</td><td>278.46</td><td>0.76</td><td>5.46</td><td>6.81</td><td>19.48</td><td>1.09</td></tr><tr><td>E7<sup>300</sup></td><td>20.96</td><td>2.71</td><td>...</td><td>...</td><td>...</td><td>267.40</td><td>0.50</td><td>6.14</td><td>6.81</td><td>11.86</td><td>1.03</td></tr><tr><td colspan="13"><strong>None</strong><br>(NOF)</td></tr><tr><td>E4<sup>180</sup></td><td>7.82</td><td>9.76</td><td>...</td><td>...</td><td>...</td><td>23.59</td><td>0.25</td><td>5.35</td><td>1.47</td><td>0.30</td><td>0.31</td></tr><tr><td>E7<sup>180</sup></td><td>4.83</td><td>6.03</td><td>...</td><td>...</td><td>...</td><td>15.36</td><td>0.01</td><td>4.43</td><td>1.66</td><td>0.005</td><td>0.60</td></tr><tr><td>E4<sup>210</sup></td><td>28.17</td><td>18.47</td><td>...</td><td>...</td><td>...</td><td>53.94</td><td>0.26</td><td>5.69</td><td>0.43</td><td>1.23</td><td>0.49</td></tr><tr><td>E7<sup>210</sup></td><td>18.30</td><td>12.00</td><td>...</td><td>...</td><td>...</td><td>46.45</td><td>0.28</td><td>5.41</td><td>0.57</td><td>0.79</td><td>0.46</td></tr><tr><td>E4<sup>250</sup></td><td>84.20</td><td>25.38</td><td>...</td><td>...</td><td>...</td><td>115.42</td><td>0.53</td><td>5.71</td><td>0.43</td><td>4.70</td><td>0.67</td></tr><tr><td>E7<sup>250</sup></td><td>69.65</td><td>21.00</td><td>...</td><td>...</td><td>...</td><td>103.21</td><td>0.46</td><td>5.60</td><td>0.46</td><td>3.59</td><td>0.64</td></tr><tr><td>E4<sup>300</sup></td><td>221.23</td><td>28.66</td><td>...</td><td>...</td><td>...</td><td>245.76</td><td>1.24</td><td>5.88</td><td>0.50</td><td>20.28</td><td>0.94</td></tr><tr><td>E7<sup>300</sup></td><td>195.85</td><td>25.37</td><td>...</td><td>...</td><td>...</td><td>224.55</td><td>1.09</td><td>5.92</td><td>0.53</td><td>16.88</td><td>0.92</td></tr></tbody></table>

**Note.** (1) Name of the galaxy model, following the nomenclature of Table 1. (2) Total MBH accreted mass. (3) Percent variation of the MBH mass, with respect to the initial MBH mass $M_{BH,0}$. (4) Value of the Eddington ratio $l = L_{BH}/L_{Edd}$ with respect to which the MBH energy is emitted equally above and below $l$. (5) Duty cycle ($D$) defined as the ratio between the time spent by the MBH at $l > 0.05$ and the total simulation time. (6) Time at which half of the total MBH radiation energy, emitted over 2–13 Gyr, has been emitted, measured since the birth of the galaxy ($t = 0$). (7–8) Total mass of stars produced and star formation rate (SFR) at 13 Gyr. (9) Time at which half of the new stellar mass, produced over the simulation time-lapse of 2–13 Gyr, has been created, measured since the birth of the original stellar population. (10) Radius containing half of the stars produced by $t = 13$ Gyr. (11) ISM luminosity within $5 R_e$ in the 0.3–8 keV band. (12) Emission-weighted temperature within $5 R_e$, and in the 0.3–8 keV band, calculated as in Equation (42).

We assume that each galaxy at the beginning of the simulation is 2 Gyr old, and is depleted of gas due to the intense high SF occurring during the initial stages of its evolution. The simulations follow the galaxy evolution for the subsequent 11 Gyr, thus ending at a galaxy age of 13 Gyr.

## 5. RESULTS

We present here the main results for the whole set of eight galaxy models in Table 1. For each galaxy model we ran three simulations: one without any feedback from the MBH (we refer to these models as NOF models), one with mechanical feedback only (MF models), and one where the feedback is radiative plus mechanical, due to an AGN conical wind (full-feedback FF models). Due to the simplified treatment of radiation feedback, we ran models with only radiative feedback just for a few exploratory cases (briefly discussed in Section 6).

The most interesting quantities at the end of the 24 simulations are listed in Table 2. Rather than presenting the specific evolution of the gas flow in each of the various models, in the following we focus on a representative model, and then we discuss the overall results across the whole set, as a function of galaxy mass, shape, and type of feedback. We focus on the results concerning the hot gas (Section 5.1), the MBH growth (Section 5.2), the newly formed stars (Section 5.3), and the duty cycle (Section 5.4).

### 5.1. Gas Evolution

During quiescent phases the flow is characterized by a central inflowing region, and an external outflowing one, as found in our previous studies (e.g., Ciotti et al. 2010 for spherical models; Negri et al. 2014b for 2D simulations without AGN feedback). The main effect of varying the galaxy mass is
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Time evolution of the main gas flow properties for the E4$^{210}$ galaxy model (see Tables 1 and 2). From left to right, plotted quantities are the X-ray luminosity and emission weighted temperature, in the 0.3–8 keV band, within $5R_e$, and the mass of the hot gas and the SF rate within the whole numerical grid. Blue, red, and black lines refer to no feedback (NOF), mechanical feedback only (MF), and full feedback (FF) models. In the bottom panels, a zoom in time is made between 10 and 10.4 Gyr.</figcaption></figure>

that the size of the central inflow is larger for more massive galaxies, with consequences for the gas content, the average gas density and cooling time, the mass accretion rate, and related quantities (as discussed below). At fixed galaxy mass, a change in the galaxy shape produces a lower $L_X$ for flatter morphologies, because the outflowing region becomes larger, as obtained for 2D models without feedback (Negri et al. 2014b).

We start by considering a representative E4 model of average mass (E4$^{210}$ in Table 1). This choice is motivated by the fact that the morphological E4 type of ETGs is observed to be more common than the E7 one, and also that our galaxy models (that are essentially non-rotating) are more realistic for E4 galaxies than for E7 ones. Figure 1 shows the time evolution of some quantities describing the hot gas phase in the E4$^{210}$ galaxy: $L_X$ and $T_X$ calculated within $5R_e$, the hot gas mass within the whole grid, and the SFR. The NOF model shows a smooth evolution in $L_X$, $T_X$, hot gas mass $M_h$ (the mass with $T > 10^6$ K), and SFR; in particular, in each panel of Figure 1, the NOF curves represent a sort of average lower envelope for the largely fluctuating behavior typical of feedback models. The present-epoch values of $L_X$, $T_X$, $M_h$, and SFR are not significantly different between the NOF, MF, and FF cases; this is true for the whole set of models, at all galaxy masses and shapes (see also Table 2). Therefore, we can conclude that the global properties of the hot gas are not significantly affected by AGN feedback. However, we note that after each outburst $L_X$ of feedback models, when calculated within $(1-2)R_e$, there is a drop that reaches one order of magnitude below $L_X$ of the NOF models, i.e., that is much larger than those in Figure 1. The

drops are the natural consequence of the clearing of the gas in the central regions.

We also varied the opening angle of the nuclear wind, changing from the adopted $n = 2$ to $n = 0$ (spherical emission) up to $n = 4$. Note that the conical wind acts over the whole spherical solid angle (reducing to zero on the equatorial plane only), independently of the value of $n$. In fact, the final results are not significantly different.

Figure 2 shows the $L_X$ versus $T_X$ relation for all models, at the present epoch, in comparison with the observed $L_X$ and $T_X$ values recently measured for the gas only, using Chandra data (Kim & Fabbiano 2015; Goulding et al. 2016). The figure shows how the NOF, MF, and FF models occupy similar regions in the plot, and also reveals a remarkable agreement between models and real ETGs. In particular, the large spread in $L_X$ at low luminosities is well reproduced, due to mostly outflowing ETGs at lower galaxy masses. Also well reproduced is the observed average trend for $L_X \gtrsim 10^{40}$ erg s$^{-1}$ (note that most of the highest $L_X$ observed are due to central ETGs in groups, Goulding et al. 2016).

The time evolution of $L_X$, $T_X$, $M_h$, and SFR is similarly rich in large fluctuations in models more massive than the E4$^{210}$ one, while it generally becomes smooth and slowly declining in less massive models, where a global outflow is established within the present epoch. Massive ETGs, then, may appear to spend much of their lifetime at very high $L_X$, $T_X$, and SFR values (Figure 1). This aspect can be checked quantitatively by computing the gas “duty-cycle,” as the fraction of time, during a chosen range of time, spent by a certain gas quantity at a level above a chosen average value for the NOF case. These fractions turn out to be very low. For example, for the gas $L_X$,
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Relation between hot gas luminosity $L_X$ and temperature $T_X$, calculated at the present epoch, for the models in Table 2. Also shown with error bars are hot gas luminosities and temperatures for observed ETGs, recently calculated using *Chandra* data for the whole emission region, and after subtraction of resolved point-sources and unresolved stellar emission (from Kim & Fabbiano 2015, in green, and Goulding et al. 2016, in blue). Most ETGs with $L_X > 10^{41}$ erg s$^{-1}$ are central galaxies in groups.</figcaption></figure>

one can compute the fraction of time spent above $L_X = 3 \times 10^{41}$ erg s$^{-1}$, that represents a sort of upper value for the observed gas emission of normal ETGs in the local universe<fnref n="9" /> (Kim & Fabbiano 2015), over the past 3 or 5 Gyr. For both choices of lookback time, the $E4^{250}$ and $E7^{250}$ FF models spend ~1% of the time with $L_X > 3 \times 10^{41}$ erg s$^{-1}$, and lower mass models spend <1% of their time above it (the very massive $E4^{300}$ and $E7^{300}$ models have $L_X$ precisely of the order of $2-3 \times 10^{41}$ erg s$^{-1}$ in the past few Gyr; see Table 2). Of course, the fraction of time during which disturbances in the gas remain can be larger than this, as already discussed for spherical models by Pellegrini et al. (2012).

The large fluctuations in the gas properties are produced by nuclear outbursts. For illustrative purposes, we have selected one representative outburst, well isolated in time, taking place at 6.85 Gyr for the $E4^{180}$ FF model; in Figures 3–7 we show the evolution of the main hydrodynamical quantities during this outburst. The leftmost panels in these figures show the quiescent time closest to the outburst, when the gas properties still have a smooth distribution. In the middle panels, corresponding to subsequent times, the gas reaches its peak in emission, and cold and dense fingers are approaching the galactic center (Figures 3 and 4), where outflow and inflow regions coexist (Figure 5); these cold filaments are mixed with hot and low density regions already created by the developing outburst. The outburst is then fading; there is still some outflowing material from the nucleus, and some hot, low density outflowing material at a radius of ≲1 kpc. The rightmost panels show again a quiescent state, at a time of 6.95 Gyr, i.e., $10^8$ years after the start of the burst; the galaxy

still has a lower density and higher temperature gas than right before the outburst. Close to the nucleus one can now notice the heating effect of the fading AGN wind, as a slightly hotter, lower density bi-conical region (Figures 3 and 4); the wind itself is outflowing in a small hourglass region above and below the nucleus (Figure 5), and imparts a tangential velocity to the gas surrounding it (Figure 6).

Figure 7 shows the SFR during the outburst: SF is very active in the cold filaments close to the center; at the end of the burst, one can see a lower SFR region close to the nucleus, due to the AGN wind. This feature nicely agrees with an observed spatial anti-correlation between Hα emission, a tracer of SF, and line emission from powerful outflows in quasars (Cano-Díaz et al. 2012; Carniani et al. 2015; Cresci et al. 2015). This observation is considered important evidence for negative feedback. However we find that overall the feedback effect on SF is positive (Section 5.3). Figure 8 shows finally the X-ray surface brightness and projected temperature maps weighted with the emission over 0.3–8 keV (Section 4), corresponding to the two central times during the outburst of Figures 3–7 ($t = 6.85$ and $t = 6.86$ Gyr). A very bright central region of ~100 pc radius is apparent in the brightness map; sharp and very hot arcs can be seen in the temperature map. Both features are very transient, but should be detectable with *Chandra*'s high angular resolution in nearby galaxies; similar features were found in 1D simulations discussed in Pellegrini et al. (2012).

We note that the bursting activity is almost continuous and lasts for the whole lifetime of massive models, at odds with the activity shown by 1D hydrodynamical simulations that was characterized by peaks well isolated in time and decreasing in frequency at later times (Ciotti & Ostriker 2012). This major difference is due to obvious geometrical reasons: in spherical symmetry the cold shell, produced by the snowplow phase of AGN-launched shock waves, cannot fragment; as a result, the shock waves are very efficient in clearing the inner regions of the galaxy from the gas. Instead, 2D hydrodynamics allows for Rayleigh–Taylor instabilities of the cold shell that breaks “permitting cold fingers of material to accrete onto the BH” (Ciotti & Ostriker 2001), while hot gas escapes from the center at the same time. This behavior was also found in spherical models simulated with the 2D code (Novak et al. 2012); the possibility of multiphase and cold gas accretion onto an MBH has been intensively investigated recently with a variety of numerical simulations (e.g., Pizzolato & Soker 2010; Barai et al. 2012; Nayakshin & Zubovas 2012; Gaspari et al. 2013).

Finally, during a burst episode we find no major differences between the FF and MF cases in the behavior of the hydrodynamical quantities on the galactic scale; the time evolution of central or nuclear properties, instead, as the gas emission during its peaks, or the nuclear emission, is more structured in the FF than in the MF case. This results in a longer duration, on average, of the outbursts in FF models, as shown in the lower panels of Figure 1, that gives an example of this difference in time evolution, with a zoom in between 10 and 10.4 Gyr. This property, already shown by spherical models (Ciotti et al. 2009, 2010), has interesting consequences for the cumulative SF of MF and FF models (Section 5.3).

## 5.2. The Black Hole Mass Growth

Figure 9 shows the adopted initial $M_{BH,0}-M_*$ relation (Section 2.2), the $M_{BH}-M_*$ relation derived from dynamical

<footnote n="9">⁹ Note that this is a conservative assumption; for example, central group or cluster ETGs can show much larger $L_X$ values.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Map of the gas number density (in cm⁻³) for the E4¹⁸⁰ FF model, at four representative times in correspondance with an outburst at 6.85 Gyr. From left to right: *t* = 6.84 Gyr (smooth appearance, immediately before the outburst), *t* = 6.85 Gyr (close to the peak in emission; cold and dense fingers are approaching the galactic center, where outflow and inflow regions coexist; these are mixed with hot and low density regions already created by the outburst; see the next Figure 4 for the mentioned features in the temperature, and Figure 5 for those in the velocity), *t* = 6.86 Gyr (the outburst is fading), and *t* = 6.95 Gyr (the main outburst effects have vanished, the galaxy is left with less dense gas than in the leftmost panel, before the outburst; a less dense bi-conical region close to the nucleus is evident, produced by the AGN wind). The solid lines represent the optical isophotes.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Map of the gas temperature (in K) for the same E4¹⁸⁰ FF model during the outburst at 6.85 Gyr of the previous Figure 3, at the same representative times. From left to right: *t* = 6.84 Gyr (smooth appearance, immediately before the outburst), *t* = 6.85 Gyr (close to the peak in emission; cold fingers are approaching the galactic center, while very hot regions have already been created by the outburst), *t* = 6.86 Gyr (the outburst is fading; there is still some hot material outflowing at a radius of ~1 kpc), and *t* = 6.95 Gyr (the main outburst effects have vanished; the heating effect of the conical wind of the fading AGN is visible as an inner slightly hotter bi-conical region).</figcaption></figure>

studies of well observed local ETGs (McConnell & Ma 2013), and the $M_{BH}$ values at the end of the simulations, when the MBH mass has grown; different symbols distinguish the galaxy shapes and types of feedback. The first thing to notice is that

models with feedback produce a final $M_{BH}-M_*$ relation that agrees well with the observed one to within the uncertainties. This agreement is not fulfilled by NOF models that end with overwhelmingly large $M_{BH}$ masses (not shown in the figure,
<figure><img src="image_6.png" /><figcaption>**Figure 5.** The same as in the previous Figures 3 and 4, for the map of the radial component ($u_r$) of the gas velocity (in km s$^{-1}$). Redder regions are outflows, bluer regions indicate inflows. From left to right: $t = 6.84$ Gyr (immediately before the outburst), $t = 6.85$ Gyr (close to the peak in emission; outflow (hot) and inflow (cold) regions coexist around the galactic center), $t = 6.86$ Gyr (the outburst is fading; there is still some outflowing gas from the center), and $t = 6.95$ Gyr (the main outburst effects have vanished, the conical wind of the fading AGN becomes visible as an hourglass feature at the nucleus).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Map of the meridional tangential component ($u_\theta$) of the gas velocity (in km s$^{-1}$) for the outburst of the previous Figures 3–5. Redder regions are moving clockwise, bluer ones are moving counter-clockwise. From left to right: $t = 6.84$ Gyr (smooth appearance, immediately before the outburst), $t = 6.85$ Gyr (close to the peak in emission), $t = 6.86$ Gyr (the outburst is fading), and $t = 6.95$ Gyr (the main outburst effects have vanished; note close to the nucleus the tangential motions imparted to the flow by the AGN conical wind).</figcaption></figure>

but see Table 2), thus demonstrating once again how feedback is needed to keep the MBH masses reasonable, even after ETGs have become “red and dead.”

MBH masses measured locally and those of feedback models are consistent at all galactic $M_*$ and shapes; in fact, the differences between the final $M_{BH}$ values at fixed $M_*$, due to shape and type of feedback, are minor. This shows that the

effect on the MBH mass growth of mechanical feedback is very important (even though this does not imply that radiative feedback alone has minor effects). Looking deeper into Figure 9, however, one can notice some trends at fixed $M_*$. E7 models show lower $M_{BH}$ than E4 ones, a result due to the lower gas binding energy of flatter galaxies (Ciotti & Pellegrini 1996; Posacki et al. 2013), with a consequent larger outflow
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Map of the SFR density (in $M_{\odot}$ kpc$^{-3}$ yr$^{-1}$) for the E4$^{180}$ FF model during the outburst of the previous Figures 3–6. From left to right: $t = 6.84$ Gyr (immediately before the outburst), $t = 6.85$ Gyr (close to the peak in emission), $t = 6.86$ Gyr (the outburst is fading), and $t = 6.95$ Gyr (the main outburst effects have vanished).</figcaption></figure>

region with respect to the central inflowing one, and larger effectiveness of AGN feedback, leading to a lower accreted mass. For the same galaxy shape, FF models end with lower $M_{BH}$ than MF ones (Figure 9), a result of their better capability of preventing gas from accreting. As a future investigation, it would be interesting to study whether the final MBH masses keep close to those measured locally even when starting from different (lower) values for $M_{BH,0}$, and more generally what is the effect of the initial (albeit at $z \sim 2$) value of $M_{BH,0}$ on the flow evolution and on the final $M_{BH}$. However, simulations with smaller $M_{BH,0}$ are more time-consuming, if, as done in all our studies, one wants to resolve the Bondi radius, and the even smaller radius at which the thermal energy associated with the AGN spectral temperature equals the particle's gravitational energy in the MBH potential (Ciotti & Ostriker 2001).

The increase in $M_{BH}$ for feedback models is different for the different galaxy masses: it ranges from a figure of the order of ~10% for the lowest mass models (E$^{180}$) to roughly ~100% (the MBH mass ~doubles) for the E$^{210}$ models, to roughly triple for the E$^{250}$ models, to almost quadruple for the very massive E$^{300}$ galaxies (see col. 3 in Table 2, and Figure 9, right panel). Therefore, MBHs in more massive ETGs not only start larger and are still larger at the end of the simulations, but are also able to grow more in percentage with respect to their initial mass. This is because more massive MBHs have larger Eddington luminosities, and, most importantly, they reside in ETGs where the gas is more bound per unit mass; then, more massive ETGs have larger accreting gas mass per unit stellar mass, or per unit MBH mass. At fixed $M_*$, the right panel of Figure 9 shows the same trends shown in the left panel: the MBH mass increases more, even with respect to the initial $M_{BH,0}$, in E4 models than in E7 ones, and in MF models than in FF ones. Such a mass increase may appear large, but in fact is remarkably small: in the absence of AGN feedback and SNIa assisted galactic winds, the mass increase of $M_{BH,0}$ would be larger by up to two orders of magnitude (Table 2). It is a

valuable property of the present models that the average increase in $M_{BH,0}$ since $z \sim 2$ is just of a factor of a few (2–4), as found for the average growth history of MBHs with given starting masses (e.g., Marconi et al. 2004).

### 5.3. The Newly Formed Stars and the Circumgalactic Medium

Figure 10 (left panel) shows the mass in newly formed stars at the end of the simulations for models with feedback. $\Delta M_*$ ranges from $\simeq 6 \times 10^8 M_{\odot}$ to $\simeq 3 \times 10^{10} M_{\odot}$, from the low to the high galaxy masses. The lowest mass models show significantly lower $\Delta M_*$ because during their evolution the gas is mostly outflowing, they experience little accreting mass and little possibility of SF. At fixed $M_*$, rounder models produce slightly larger $\Delta M_*$ than flatter ones, due to their larger capability to retain the gas, that can eventually be converted into stars. MF and FF models produce similar $\Delta M_*$, at any fixed $M_*$ and galaxy shape, with larger $\Delta M_*$ in the FF case. The better efficiency of FF models to produce SF is a consequence of the richer structure and longer duration of each of their nuclear outburst episodes, which prevent prompt accretion, provide the gas more time for its cooling, and lead to a better capability to keep cold and dense gas far from the nucleus (see Section 5.1, and Figure 1, where an FF outburst shows more fluctuations in $L_X$, and a larger duration of the feedback action on the gas, than an MF outburst). Note that the final $\Delta M_*$, at fixed $M_*$ and galaxy shape, increases from NOF to MF to FF models (Table 2, except for the two NOF E$^{180}$ models), which proves definitively that feedback has the overall effect of making SF more efficient, i.e., it has a positive action with respect to SF. Also, the behavior of the $\Delta M_*$ growth is opposite to that of the MBH growth with respect to feedback: $\Delta M_*$ increases from NOF to FF models, while the opposite is true for $\Delta M_{BH}$ (cf. Figures 9 and 10). Both $\Delta M_*$ and $\Delta M_{BH}$, instead, are larger for rounder models.
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Surface brightness maps in the 0.3–8 keV band (upper panels), and projected temperature maps (lower panels), for the same model in the outburst of the previous Figures 3–7, at the two central times in those panels: *t* = 6.85 Gyr (the outburst has just started), and *t* = 6.86 Gyr (the outburst is fading). Solid contours indicate the optical isophotes.</figcaption></figure>

It is interesting to compare the relative increase of $\Delta M_*$ and $\Delta M_{BH}$ (values in Table 2). In NOF models, $\Delta M_{BH}$ and $\Delta M_*$ are not largely different, and MBH growth is much favored over SF with respect to feedback models. For feedback models, $\Delta M_{BH}$ ranges from $\simeq 6 \times 10^6 M_\odot$ to $\simeq 2 \times 10^9 M_\odot$, thus $\Delta M_*$ is always much larger than $\Delta M_{BH}$ at all galaxy masses. There is a trend though: $\Delta M_{BH}/M_{BH,0}$ (Figure 9) increases clearly with $M_*$, while $\Delta M_*/M_*$ keeps quite flat or even decreases (Figure 10, right panel). Thus, $\Delta M_{BH}/\Delta M_*$ increases with $M_*$, and there is relatively more MBH growth than SF in more massive galaxies. The reason for this trend is that the mass flowing to the central region, where SF takes place (see below), is converted in new stars more efficiently than in MBH mass at low galaxy masses, while the opposite works in the most massive galaxies. In more massive ETGs, in fact, the feedback can displace the gas out to distances from the MBH that are relatively lower than in less massive ones, consequently the gas has less time to fall back to the center, while giving origin to SF

(as revealed by an inspection of the evolution of the hydrodynamical properties of the flow). One final remark here concerns the Magorrian relation $M_{BH}-M_*$: as the stellar mass does not increase significantly after $z \sim 2$ (Figure 10, right panel), and $M_{BH}$ can increase up to a factor of 4 (Figure 9, right panel), this implies that after $z \sim 2$ the Magorrian relation is expected to shift upwards in massive ETGs, just as a consequence of the passive evolution of the stellar population. The extent of the shift is small, though, within the scatter and uncertainties of the local and low redshift observed relations (e.g., McConnell & Ma 2013; Schulze & Wisotzki 2014).

It is also interesting to compare $\Delta M_*$ with the total stellar mass losses from the beginning to the end of the simulation, i.e., from an age of 2 Gyr for the stellar population to that of 13 Gyr. From an integration of the rate in Equation (1), these losses amount to $\sim 10\%$ of $M_*$. From Figure 10 we see that $\Delta M_*$ is $\simeq 5\%$ of $M_*$ (or lower for the least massive models). Thus, roughly half of the integrated stellar mass losses within a galaxy since $z \sim 2$ goes into new stars. The remaining part of the mass losses goes for a minor fraction into $\Delta M_{BH}$ (that reaches at most $3 \times 10^{-3} M_*$, and then at most 3% of the integrated stellar mass losses, see Figure 9), and for a major fraction is ejected from the galaxy, due to SN heating and the further help of the AGN feedback action. We recall that AGN feedback produces an increase in the ejected mass from the galaxy, but has never been found capable of producing a global/major outflow by itself during passive galaxy evolution after $z \sim 2$. To better quantify this increase, we can evaluate the amount of gas residing outside a reference radius of $5R_\odot$ (that we consider lost by the galaxy), at the end of the simulations, for models with and without feedback. This amount is larger for feedback models than for NOF models by a percentage ranging from $\sim 20\%$ to $\sim 40\%$, from the largest to the smallest galaxy masses. We can also compute the “load factor,” defined as the ratio between the gas ejected from the galaxy during the evolution, and the mass in new stars $\Delta M_*$; such a factor is useful to establish the possible role of AGN feedback in adding material to the circumgalactic medium, after $z \sim 2$. For feedback models, the load factor is $\simeq 0.6$, except for the lowest mass (E$^{180}$) galaxies, that eject a larger amount of gas, and for which then the load factor is $\simeq 3$. For NOF models, the factor is slightly larger ($\simeq 0.7$), but not so different, because they eject less mass, but also form fewer stars; a markedly lower factor ($\simeq 1$) than in feedback models is instead shown by the E$^{180}$ NOF models, where the energy input by the feedback can be significant in clearing the gas from the galaxy. Note that these load factors are similar to those ($\sim 0.7$) recently determined thanks to background quasar lines of sight passing near star-forming galaxies (e.g., Schroetter et al. 2016).

Figure 11 shows the epoch at which half of $\Delta M_*$ is formed ($t_M$ from Table 2). This epoch (measured since the birth of the original stellar population, i.e., 2 Gyr before the start of the simulation) keeps between 5 and 6 Gyr for most models, and tends to be larger for MF than for FF models. Thus, half of the new stars created during the life period of 2–13 Gyr form quite early, during the first 3–4 Gyr. Figure 1 shows the time evolution of the SFR for the E4$^{210}$ model, and Table 2 gives the instantaneous value of the SFR at the end of the simulations. These present-epoch rates typically keep below $\simeq 0.5 M_\odot \text{ yr}^{-1}$, and compare well with the current rates observed for molecular gas-rich ETGs of the local ATLAS$^{3D}$ sample, which range from 0.01 to $3 M_\odot \text{ yr}^{-1}$, with a median value of 0.15 $M_\odot \text{ yr}^{-1}$.
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Left panel: final values of the MBH masses vs. $M_*$ (given in Table 1). Plotted are the 16 models with feedback in Table 2: E4 in black, E7 in red; open symbols are MF models, full symbols are FF ones. The positions of the models are compared with the initial $M_{BH,0}-M_*$ relation (dotted line), and the observed $M_{BH}-M_*$ relation (solid line), derived for dynamically measured stellar masses of local ETGs (McConnell & Ma 2013). The vertical bar shows the estimated intrinsic scatter in log log $M_{BH}$ for the plotted scaling relation. Right panel: the percental increase in the MBH mass, for the same feedback models on the left, with the same meaning of symbols.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Left panel: final mass in the newly formed stars $\Delta M_*$ (from Table 2) vs. $M_*$ (from Table 1), for models with feedback (E4 in black, E7 in red). Plotted are the 16 models with feedback in Table 2; MF and FF models are shown with open and full symbols, respectively. Right panel: the percentage increase in the stellar mass, for the same models on the left, with the same meaning of the symbols.</figcaption></figure>

(Davis et al. 2014). Larger SFRs than ≈0.5 $M_\odot$ yr⁻¹ are shown by models caught in an outburst (Figure 1), or in the most massive models, that host more gas, and have a larger AGN activity; note that the ATLAS³D survey is not representative of the most massive models in Table 2.

Another interesting aspect for comparison with observations is *where* SF mainly takes place. Figure 12 shows meridional sections of the ratio between the density in newly formed stars, at the end of the simulations, and that in the original stellar population, for the E4 shape (results are similar for the E7 case). SF is very low in the lower mass model (E4$^{180}$), that experiences an almost global outflow over its whole lifetime. SF is instead evident in the larger mass E4$^{250}$ model: it forms a

nuclear stellar disk in the NOF case,<fnref n="10" /> while it has a roughly spherical distribution, peaking within ≈1 kpc radius in the cases with feedback. The extent of the SF is measured more quantitatively by $r_M$, the radius including half of the new stars at the end of the simulation (in Table 2). Typically $r_M \simeq R_c/2$, for feedback models, except for the least massive ones, where SF is very low but extended. In contrast, $r_M$ is much smaller for NOF models, indicating how most of SF takes place at the center, i.e., the gas flows to the central regions before having time to start SF. In feedback models, instead, as already noted,

<footnote n="10">¹⁰ This disk is the result of the stellar streaming imposed in the central regions for numerical reasons, and discussed in Section 4, see Equation (38).</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Epoch (measured since the birth of the original stellar population of the galaxy, i.e., 2 Gyr before the start of the simulation) at which half of the final mass in the newly formed stars $\Delta M_{*}$ is formed ($t_{M}$ in Table 2), vs. $M_{*}$ (in Table 1); E4 models are plotted in black, E7 in red. Plotted are the 16 feedback models in Table 2; MF and FF models are shown with open and full symbols, respectively.</figcaption></figure>

the gas is compressed from time to time by AGN activity, and kept at larger distances by the feedback action, which produces longer times available for SF, that is then favored. The feedback models also show an X-shaped structure in SF, inclined by $\sim 45^{\circ}$, more evident in the MF case. This is related with the aperture of the conical wind: at the contact surface between the wind region and the external gas, the density is increased, with a consequent increase in SF. Note that in the simulations the newly formed stars remain in the position where they are born, while in reality they move away from that place. From this point of view, it is not clear whether these features could be evident in real ETGs; however, curiously, X-shaped features have been observed in the morphology of bulges (although on a much larger scale; e.g., Ness & Lang 2016), and in scattered light produced by illuminated dusty cones in quasars hosting winds (e.g., Obied et al. 2016).

When seen in the surface brightness profile, these new stars produce a central cusp, as observed in ETGs at high angular resolution (see Kormendy et al. 2009 for a review on these nuclear cusps or cuspy cores). Even coreless galaxies do not have featureless power-law profiles, but, rather, they show central extra light above the inward extrapolation of the outer Sérsic profile. Kormendy et al. (2009) and Hopkins et al. (2009) suggested that the extra light is produced by starbursts fed by gas dumped inward during dissipative mergers. The origin of the observed cusps could also reside in the evolutionary phenomena investigated here (see also Ciotti & Ostriker 2007; Ciotti 2009).

### 5.4. The Evolution of $\dot{M}_{BH}$, $L_{BH}$, and the Duty Cycle

Figure 13 shows the time evolutions of $\dot{M}_{BH}$ (Equation (18)) and $L_{BH}$ (Equation (23)) that are relevant to quantify the accretion history of the MBH for the representative model E4$^{210}$ in Figure 1. In feedback models, the average accretion rate slowly decreases with time, in pace with the declining rate of mass input from stars. The accretion episodes may extend throughout the galaxy lifetime, with $\dot{M}_{BH} \gtrsim 10^{-4} M_{\odot} \text{ yr}^{-1}$, as

for more massive models, or for a shorter time, and with smaller $\dot{M}_{BH}$, as for lower mass models (E$^{180}$).

A major difference with respect to 1D hydrodynamical simulations, where major bursts are well separated in time (Ciotti & Ostriker 2012), is that the bursting activity is here almost continuous (as described in Section 5.1). In spite of this difference between the bursting activity of 1D and 2D hydrodynamics, the resulting duty-cycle is not so different (both in spherical and flat galaxies). In order to quantify this important aspect, we considered the time spent by each model at any chosen Eddington ratio $l$, and with $l$ above/below a chosen threshold. In particular, for the E4 and E7 FF models, Figure 14 shows the fraction of total simulation time, and the fraction of total radiative energy, respectively, spent and emitted at different values of $l$. For all models the distribution of times has a long, almost flat tail extending to the lowest $l < 10^{-10}$, and peaks at an $l$ that increases with the galaxy mass: the peak is located at $l \simeq 10^{-3}$, and reaches $l \simeq 10^{-2}$ for the most massive models (while becoming less and less pronounced). All distributions drop sharply above $l \simeq 10^{-2}$, so that all models spend very little time above $l \simeq 10^{-1}$. An average, realistic galaxy, as could be described by the E4$^{210}$ FF model, spends 75% of its time since $z \sim 2.2$ below $l \approx 10^{-3}$. Note however how the fraction of the energy emitted at high Eddington ratios is significant. Even though the nuclei are usually very faint, all models typically emit 25% of the energy above $l = 0.1$, as shown by the right column in Figure 14.

Figure 15 shows the percentage of the total simulation time (11 Gyr) and of the total emitted energy, respectively, spent and emitted above and below each value of $l$, for the E4 FF models (the same figure for the E7 FF models, not shown, presents the same trends, both in shapes and in normalizations). Estimates based on 1D models, and with a luminosity threshold of $L_{BH} = L_{\text{Edd}}/30$, gave a figure of a few percent for the fraction of time spent above this threshold (Ciotti et al. 2010). Remarkably, we can confirm this estimate, with a percentage of time spent above $l = 1/30$ ranging from 3 to 5%, going from the E4$^{210}$ to the E4$^{300}$ models (Figure 15, left panels, solid lines).

For reference, Table 2 gives the fraction of time spent above $l = 0.05$ ($D$); the $l$-value with respect to which the MBH energy is emitted equally above and below $l$ ($l_{0.5}$); and the time $t_{L}$ at which half of the total MBH radiation energy, emitted over 2–13 Gyr, has been emitted, measured since the birth of the galaxy (i.e., 2 Gyr before the start of the simulation).

The duty-cycles resulting from the present set of simulations, in the form shown in Figures 14 and 15, represent a useful diagnostic that could be used in observational works, and phenomenological models, that try to derive the growth rate of each MBH mass, the evolution of their Eddington ratio, the total luminosity emitted as a function of $l$, and the duty cycle of AGN activity (e.g., Marconi et al. 2004; Trakhtenbrot & Netzer 2012; Caplar et al. 2015; Schulze et al. 2015). We stress that our predictions concern massive non-rotating ETGs, with $M_{BH} \gtrsim 10^8 M_{\odot}$, since $z \sim 2.2$, i.e., after galaxy formation has ended. How much the rotation of a galaxy could eventually affect these findings will be addressed in a forthcoming paper.

## 6. SUMMARY AND CONCLUSIONS

In this work we have followed the evolution of hot gas flows in ETGs with central MBHs, using 2D hydrodynamical simulations where the most accurate and physically consistent description of AGN feedback (both radiative and mechanical,
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Meridional sections of the ratio between the density in newly formed stars at the end of the simulations and that in the original stellar population. Solid lines show constant densities for the original stellar population. From left to right the panels refer to: the FF E4$^{180}$ model, and to the NOF, MF, and FF E4$^{250}$ models. SF is very low in the lower mass model that experiences an almost global outflow over its whole lifetime; SF is instead significant in the larger mass galaxy, and it forms a nuclear stellar disk in the NOF case, while it has a roughly spherical distribution in the cases with feedback (see Section 5.3 for more details).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 13.** Time evolution of the mass accretion rate on the MBH ($\dot{M}_{\text{BH}}$, left), and the radiated accretion luminosity ($L_{\text{BH}}$, right), for the E4$^{210}$ model (for which Figure 1 shows the gas evolution). Blue, red, and black lines refer to NOF, MF, and FF models, respectively; in the right panel the nuclear luminosity is shown only in the FF case. The lower panels show the same zoom in time as in Figure 1.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure 14.** Left panels: percentage of the total simulation time (11 Gyr) spent at the values of the Eddington ratio on the x-axis, for the E4 FF models. Right panels: for the same models on the left, percentage of the total energy emitted at the Eddington ratio on the x-axis. In each panel, for the E4 models only, the vertical dashed lines mark the Eddington ratios below which the model spends 25%, 50%, and 75% of the total time (left), or below which 25%, 50%, and 75% of the total energy is emitted (right).</figcaption></figure>

due to AGN winds) is implemented, the spatial resolution at the center is parsec-scale, and the underlying galaxy models are the most realistic dynamical models adopted so far in this kind of study. The mass and energy input from the stellar population are secularly evolving, according to the prescriptions of stellar

evolution theory for the stellar mass losses, and to the predictions of progenitors evolution and to observations, for the declining SNIa rate. SF is implemented via a simple scheme based on physical arguments shown to reproduce well the Kennicutt-Schmidt relation (Negri et al. 2015). The galaxy
<figure><img src="image_16.png" /><figcaption>**Figure 15.** Left panels: percentage of the total simulation time (11 Gyr) spent above (solid line) and below (dashed line) the values of the Eddington ratio on the x-axis, for the E4 (in black) and E7 (in red) FF models. Right panels: for the same models on the left, percentage of the total energy emitted above (solid line) and below (dashed line) the Eddington ratio on the x-axis. In each panel, the vertical dashed lines mark the Eddington ratios below which the model spends 25%, 50%, and 75% of the total time (left), or below which 25%, 50%, and 75% of the total energy is emitted (right).</figcaption></figure>

models are axisymmetric, and the Jeans equations provide detailed internal stellar kinematics on which stellar kinematical heating is based. The stellar density profile follows a deprojected (ellipsoidal) Sérsic law, and the main observables ($L, R_e, \sigma_{e8}$) are chosen in order for the galaxy models to lie on the main scaling laws. The spherical DM halo has a radial profile predicted by cosmological simulations (NFW), and is normalized to account for a dark mass within $R_e$ lower than the

stellar mass, as observed. A few limitations remain in the present work, though. First of all, significant galactic rotation is not included; rotation is known to lead to an almost corotating gaseous halo (Negri et al. 2014b), which in turn produces a massive rotating equatorial cold disk, of ~kpc scale, that is expected to be unstable, with consequent disposal of fresh gas on the central MBH. A new timescale is then introduced, determined by the global stability of the disk, and the amount
of gas accreted on the MBH depends also on the SF taking place in the disk. Second, the radiative transport is here calculated by using the simplified (and computationally much faster) scheme described in Novak et al. (2012). Also not included are the effects related with the presence of dust (radiative transport, reprocessing of the radiation from the new stellar population, sputtering, as described in Hensley et al. 2014). Third, the circumnuclear sub-grid accretion disk, implemented in some of our previous works (e.g., Ciotti & Ostriker 2012) is not considered. Finally, the presence of a jet that seems relevant at recent epochs is still to be included. For the limitations of the adopted radiative feedback, we ran just a few test models with only radiative feedback. These showed similar properties to the MF and FF models (e.g., in the fluctuations in time of $L_X$, $T_X$, the SFR, etc.); the mass accreted on the MBH is instead much larger, yet still much lower than in the NOF models.

The main results of this work are as follows.

(1) $L_X$ and $T_X$ of the gas at the present epoch, for models with feedback, reproduce well those observed for ETGs similar to our galaxy models, at all galaxy masses; in particular, the observed large range of $L_X$ at the low $T_X$ is accounted for, thanks to the prevalence of outflows in lower mass galaxies. At fixed galaxy mass, a change in the galaxy shape produces a lower $L_X$ for flatter morphologies, because the outflowing region becomes larger, as obtained for 2D models without feedback (Negri et al. 2014b). The evolution of the gas when feedback is present is characterized by large and frequent fluctuations in $L_X$ and $T_X$, except for the lowest mass models explored here, where these fluctuations stop after the first few Gyr, due to the onset of a global outflow sustained by SNIa heating.<fnref n="11" /> During the past 3–5 Gyr, the fraction of time spent above a value representative of the largest gas emission observed for normal ETGs in the local universe ($L_X = 3 \times 10^{41} \text{ erg s}^{-1}$), is ≲1%. After each outburst, $L_X$ within (1–2)$R_e$ is lower by ≳1 order of magnitude than the same quantity for models without AGN feedback.

(2) Each major burst triggered by accretion on the MBH is made by several smaller bursts, and lasts ~few × 10⁷ yr; we stress that this timescale is not imposed a priori, but it results from the simulations. The 2D hydrodynamics and the ellipsoidal shape of the galaxy models produce a complicated gas evolution in the inner ~2–3 kpc region, where a cold and inflowing gas phase coexists with a hot and outflowing one. Bursts from mechanical or full feedback are qualitatively similar, but temporally more isolated in the purely mechanical case. The X-ray surface brightness map shows a very bright central region (on a scale of the order of ~100 pc), including sharp and very hot arcs, as a very transient feature. The development and fading out of an outburst are also similar for different galaxy shapes. We expect though that the addition of galactic rotation will change this.

(3) The mass in newly formed stars $\Delta M_*$, for feedback models, ranges from $6 \times 10^8$ to $3 \times 10^{10} M_\odot$, from the low to the high galaxy masses. Since the models are basically non-rotating, SF takes place in an almost spherical region; the radius including half of the new stars is $r_M \simeq R_e/2$, except for the least massive models, where SF is very low but extended.

Half of the mass in the new stars forms within 5–6 Gyr from the birth of the original stellar population.

(4) At fixed $M_*$, $\Delta M_*$ is larger for the full feedback case than the mechanical one, that in turn produces a larger $\Delta M_*$ than without feedback. This proves that feedback has the overall effect of making SF more efficient, at least in the phase following the major galaxy formation period studied here, when SF and MBH accretion are fueled by passive stellar evolution. These results confirm previous findings of spherical models about positive feedback (Ciotti & Ostriker 1997; e.g., see also Ishibashi & Fabian 2012; Nayakshin & Zubovas 2012); thus geometrical effects allowed for by 2D simulations and a flat galaxy shape do not invalidate the global picture.

(5) $\Delta M_*/M_*$ is of the order of 0.04–0.05, excluding the lowest mass models that show lower values (≈0.01), since most of their injected stellar mass loss is expelled in an outflow. On average, excluding the lowest mass models, roughly half of the total mass losses since an age of ≃2 Gyr ends recycled into new stars. The other half goes for a minor fraction into $\Delta M_{BH}$ (for feedback models), and for a major fraction is ejected from the galaxy, mostly due to SN heating. AGN feedback produces an increase (of the order of 20%–40%) in the ejected mass from the galaxy, during the passive galaxy evolution after $z \sim 2$. The “load factor,” defined as the ratio between the gas displaced outside $5R_e$ during the evolution, and the mass in new stars $\Delta M_*$, is ≃0.6 for feedback models, and slightly larger for NOF models (due to their weaker ability to form new stars). The lowest mass (E¹⁸⁰) galaxies eject a larger amount of gas, and for them the load factor is larger (≈3 for models with AGN feedback). We note that AGN outflows are often considered responsible for clearing gas out of massive galaxies at the end of the epoch of galaxy formation, which requires an efficiency of energy transfer to the ISM much larger than shown by our simulations (that refer to the passive evolution phase). It could be that AGN feedback at the epoch of galaxy formation, with a denser and more optically thick gas, and a shallower galactic potential well, was stronger.

(6) Final $M_{BH}$ values of feedback models are consistent with those measured in the local universe, at all $M_*$ considered here, while they are overwhelmingly large without feedback; the latter is then fundamental to keep the MBH masses at the values observed, after $z \sim 2$. In feedback models just ≃3% of the total injected mass from stars is accreted on the MBH; the average increase in $M_{BH}$ since $z \sim 2$ is just of a factor of a few (2–4). The MBHs in more massive galaxies grow more than in lower mass ones, during the later phase of (small) MBH growth studied here. At fixed $M_*$, flatter models show lower $M_{BH}$, due to the lower gas binding energy, and FF models end with lower $M_{BH}$ than MF ones, thanks to their better capability of preventing gas accretion. $\Delta M_{BH}/\Delta M_*$ ranges from $8 \times 10^{-3}$ to $7 \times 10^{-2}$, from the smallest to largest galaxy masses. The implementation of rotation could change this trend, by making MBH growth less favored, and SF in the cold disk more favored.

(7) Half of the radiative output from the AGN is emitted within $t_L \approx 5$ Gyr since the birth of the galaxy, roughly the same value of the timescale $t_M$ within which half of the new stars form; this is another indication that the two activities are linked. We notice a slight but clear dependence on galaxy mass in the time-shift between these two timescales: $t_L$ precedes $t_M$ in lower mass galaxies, while the reverse is true in larger mass ones.

<footnote n="11">¹¹ In these low mass ETGs the MBH accretes almost steadily, from a hot atmosphere, at a very sub-Eddington rate, as shown previously (e.g., Ciotti & Ostriker 2012).</footnote>
(8) The outburst activity is almost continuous, and extends to the present epoch, for more massive galaxies, or stops a few Gyr ago, for the least massive ones considered here. Almost all of the time is spent at very low nuclear luminosities (75% of the simulation time is spent at Eddington ratios $l \lesssim 10^{-3}$), the energy is emitted at high Eddington ratios, with one quarter of the total energy emitted at $L_{BH}/L_{Edd} > 0.1$. The resulting duty-cycles of AGN activity, estimated as the fraction of time spent above $L_{Edd}/30$, range from 3 to 5%.

We note that recently Cheung et al. (2016) identified a population of local ETGs (that they called “red geysers”) showing bisymmetric outflow-like structures in their ionized gas emission line maps. The authors argue that these are low-luminosity AGN-driven winds, and that they are occurring in roughly 5%–10% of the red sequence at moderate masses. This finding would prove how the AGN wind mechanical feedback continues down to late times, as predicted by our modeling.

We thank G. Novak for providing an initial version of the code used here, and S. Posacki for providing the dynamical properties of the underlying galaxy models. We also thank the referee for useful comments. L.C., S.P., and A.N. were supported by the MIUR grant PRIN 2010–2011, project “The Chemical and Dynamical Evolution of the Milky Way and Local Group Galaxies,” prot. 2010LY5N2T. AN acknowledge funding from the European Research Council under the European Community’s Seventh Framework Programme (FP7/2007–2013 Grant Agreement no. 614199, project BLACK).

REFERENCES

Alvarez, M. A., Wise, J. H., & Abel, T. 2009, ApJL, 701, L133

Arav, N., Borguet, B., Chamberlain, C., et al. 2013, MNRAS, 436, 3286

Athey, A., Bregman, J., Bregman, J., Temi, P., & Sauvage, M. 2002, ApJ, 571, 272

Bae, H.-J., & Woo, J.-H. 2016, ApJ, 828, 97

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

Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57

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

Bieri, R., Dubois, Y., Rosdahl, J., et al. 2016, MNRAS, in press (arXiv:1606.06281)

Bondi, H. 1952, MNRAS, 112, 195

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

Bregman, J. N., & Parrott, J. R. 2009, ApJ, 699, 923

Cano-Díaz, M., Maiolino, R., Marconi, A., et al. 2012, A&A, 537, L8

Caplar, N., Lilly, S. J., & Trakhtenbrot, B. 2015, ApJ, 811, 148

Cappellari, M., Romanowsky, A. J., Brodie, J. P., et al. 2015, ApJL, 804, L21

Carniani, S., Marconi, A., Maiolino, R., et al. 2015, A&A, 580, 102

Cheung, E., Bundy, K., Cappellari, M., et al. 2016, Natur, 533, 504

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

Ciotti, L. 2009, Natur, 460, 333

Ciotti, L., D'ERcole, A., Pellegrini, S., & Renzini, A. 1991, ApJ, 376, 380

Ciotti, L., & Ostriker, J. P. 1997, ApJL, 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 Hot Interstellar Matter in Elliptical Galaxies, ed. D.-W. Kim & S. Pellegrini (Berlin: Springer), 8

Ciotti, L., Ostriker, J. P., & Proga, D. 2009, ApJ, 699, 89

Ciotti, L., Ostriker, J. P., & Proga, D. 2010, ApJ, 717, 708

Ciotti, L., & Pellegrini, S. 1996, MNRAS, 279, 240

Cisternas, M., Jahnke, K., Inskip, K., et al. 2011, ApJ, 726, 57

Cooksey, K. L., Thom, C., Prochaska, J. X., & Chen, H-W. 2010, ApJ, 708, 868

Cresci, G., Marconi, A., Zibetti, S., et al. 2015, A&A, 582, A63

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

Davis, T. A., Young, L. M., Crocker, A. F., et al. 2014, MNRAS, 444, 3427 de Vaucouleurs, G. 1948, AnAp, 11, 247

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

Desroches, L.-B., Quataert, E., Ma, C.-P., & West, A. A. 2007, MNRAS, 377, 402

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

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

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

Faucher-Giguère, C., & Quataert, E. 2012, MNRAS, 425, 605

Feruglio, C., Fiore, F., Carniani, S., et al. 2015, A&A, 583, A99

Fischer, T. C., Crenshaw, D. M., Kraemer, S. B., & Schmitt, H. R. 2013, ApJS, 209, 1

Fox, A. J. 2011, ApJ, 730, 58

Gan, Z., Yuan, F., Ostriker, J. P., et al. 2014, ApJ, 789, 150

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

Gibson, R. R., Brandt, W. N., Gallagher, S. C., & Schneider, D. P. 2009, ApJ, 696, 924

Goulding, A. D., Greene, J. E., Ma, C.-P., et al. 2016, ApJ, 826, 167

Greene, J. E., Zakamska, N. L., Ho, L. C., & Barth, A. J. 2011, ApJ, 732, 9

Greggio, L. 2010, MNRAS, 406, 22

Hayes, J. C., Norman, M. L., Fiedler, R. A., et al. 2006, ApJS, 165, 188

Heckman, T. M., & Best, P. N. 2014, ARA&A, 52, 589

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

Hobbs, A., Power, C., Nayakshin, S., & King, A. R. 2012, MNRAS, 421, 3443

Hopkins, P. F., Cox, T. J., Dutta, S. N., et al. 2009, ApJS, 181, 135

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

Ishibashi, W., & Fabian, A. C. 2012, MNRAS, 427, 2998

Johansson, P. H., Naab, T., & Burkert, A. 2009, ApJ, 690, 802

Kazantzidis, S., Mayer, L., Colpi, M., et al. 2005, ApJL, 623, L67

Kim, D.-W., & Fabbiano, G. 2015, ApJ, 812, 127

King, A. 2003, ApJL, 596, L27

King, A., & Pounds, K. 2015, ARA&A, 53, 115

Kocevski, D. D., Faber, S. M., Mozena, M., et al. 2012, ApJ, 744, 148

Kormendy, J., Fisher, D. B., Cornell, M. E., & Bender, R. 2009, ApJS, 182, 216

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

Krongold, Y., Nicastro, F., Elvis, M., et al. 2007, ApJ, 659, 1022

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

Lehner, N., Howk, J. C., Tripp, T. M., et al. 2013, ApJ, 770, 138

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

Liu, G., Arav, N., & Rupke, D. N. 2015, ApJS, 221, 9

Liu, G., Zakamska, N. L., Greene, J. E., Nesvadba, N. P. H., & Liu, X. 2013, MNRAS, 436, 2576

Madau, P., & Rees, M. J. 2001, ApJL, 551, L27

Mannucci, F., Della Valle, M., Panagia, N., et al. 2005, A&A, 433, 807

Maoz, D., Mannucci, F., Li, W., et al. 2011, MNRAS, 412, 1508

Maraston, C. 2005, MNRAS, 362, 799

Marconi, A., Risaliti, G., Gilli, R., et al. 2004, MNRAS, 351, 169

Mathews, W. G. 1989, AJ, 97, 42

Mathews, W. G. 1990, ApJ, 354, 468

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

McElroy, R., Croom, S. M., Pracy, M., et al. 2015, MNRAS, 446, 2186

Mellier, Y., & Mathez, G. 1987, A&A, 175, 1

Moe, M., Arav, N., Bautista, M. A., & Korista, K. T. 2009, ApJ, 706, 525

Narayan, R., & Yi, I. 1995, ApJ, 452, 710

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493

Nayakshin, S., & Zubovas, K. 2012, MNRAS, 427, 372

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

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

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

Ness, M., & Lang, D. 2016, AJ, 152, 14

Norman, C., & Scoville, N. 1988, ApJ, 332, 124

Novak, G. S., Ostriker, J. P., & Ciotti, L. 2012, MNRAS, 427, 2743

Novak, G. S., Ostriker, J. P., & Ciottik, L. 2011, ApJ, 737, 26

Obied, G., Zakamska, N., Wylezalek, D., & Liu, G. 2016, MNRAS, 456, 2861

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

Parrott, J. R., & Bregman, J. N. 2008, ApJ, 681, 1215

Pellegrini, S. 2005, ApJ, 624, 155

Pellegrini, S. 2012, in Hot Interstellar Matter in Elliptical Galaxies, ed. D.-W. Kim & S. Pellegrini (Berlin: Springer), 21

Pellegrini, S., Ciotti, L., & Ostriker, J. P. 2012, ApJ, 744, 21

Pizzolato, F., & Soker, N. 2010, MNRAS, 408, 961

Posacki, S., Pellegrini, S., & Ciotti, L. 2013, MNRAS, 433, 2259

Proga, D., & Kallman, T. R. 2004, ApJ, 616, 688

Proga, D., Stone, J. M., & Kallman, T. R. 2000, ApJ, 543, 686
Reichard, T. A., Richards, G. T., Hall, P. B., et al. 2003, AJ, 126, 2594

Satoh, C. 1980, PASJ, 32, 41

Sazonov, S., Krivonos, R., Revnivtsev, M., Churazov, E., & Sunyaev, R. 2008, A&A, 482, 517

Sazonov, S. Y., Ostriker, J. P., Ciotti, L., & Sunyaev, R. A. 2005, MNRAS, 358, 168

Schawinski, K., Simmons, B. D., Urry, C. M., et al. 2012, MNRAS, 425, L61

Schroetter, I., Bouché, N., Wendt, M., et al. 2016, ApJ, 833, 39

Schulze, A., Bongiorno, A., Gavignaud, I., et al. 2015, MNRAS, 447, 2085

Schulze, A., & Wisotzki, L. 2014, MNRAS, 438, 3422

Sharon, K., Gal-Yam, A., Maoz, D., et al. 2010, ApJ, 718, 876

Sijacki, D., Vogelsberger, M., Genel, S., et al. 2015, MNRAS, 452, 575

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

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

Sturm, E., González-Alfonso, E., Veilleux, S., et al. 2011, ApJL, 733, L16

Tang, S., & Wang, Q. D. 2005, ApJ, 628, 205

Tombesi, F., Meléndez, M., Veilleux, S., et al. 2015, Natur, 519, 436

Tombesi, F., Sambruna, R. M., Reeves, J. N., et al. 2010, ApJ, 719, 700

Trakhtenbrot, B., & Netzer, H. 2012, MNRAS, 427, 3081

Treister, E., Schawinski, K., Urry, C. M., & Simmons, B. D. 2012, ApJL, 758, L39

Wu, X., Wang, F., Fan, X., et al. 2015, Natur, 518, 512

Zakamska, N. L., Hamann, F., Pâris, I., et al. 2016, MNRAS, 459, 3144

Zubovas, K., & King, A. 2012, ApJL, 745, L34