2D hydrodynamical simulations of the BAL wind feedback

# The Effect of the AGN Feedback on the Interstellar Medium of Early-Type Galaxies: 2D Hydrodynamical Simulations of the Low-Rotation Case.

Luca Ciotti11affiliationmark: , Silvia Pellegrini11affiliationmark: , Andrea Negri22affiliationmark: , Jeremiah P. Ostriker33affiliationmark: 44affiliationmark: 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
###### Abstract

We present 2D hydrodynamical simulations for the evolution of hot gas flows in early-type galaxies with central massive black holes (MBHs), starting from an age of Gyr; the code has an accurate and physically consistent description of radiative and mechanical (due to AGN winds) feedback, and a parsec-scale resolution at the center. The mass input for the flow comes from stellar mass losses, and the energy input includes Type Ia supernova and stellar heating; the flow can form stars. Realistic, axisymmetric dynamical models for the galaxies are built by solving the Jeans’ equations. We find that the lowest mass models explored () develop a global outflow sustained by SNIaâs heating after a few Gyr, and then end with a significantly lower amount of hot gas and new stars (). In more massive models, instead, nuclear outbursts last up to the present epoch, with large and frequent fluctuations in the emission of the nucleus and of the gas (). Each burst lasts few yr, during which, in the inner kpc region, cold, inflowing, and hot, outflowing gas phases coexist. The relation for the gas reproduces well that of local galaxies. AGN activity determines a positive feedback on star formation. Roughly half of the total mass losses ends recycled into new stars, just % of it is accreted on the MBH, and most of the remaining part is ejected from the galaxy; the ratio between the mass of gas displaced outer of 5 during the evolution and that in new stars , the load factor, is . A rounder galaxy shape corresponds to a larger final MBH mass, , and . Half of the radiative output from the AGN is emitted within Gyr from the start of the simulation, roughly the same timescale within which half of the new stars form; 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 . The duty-cycles of AGN activity ranges from 3 to 5%.

galaxies: elliptical and lenticular, cD – galaxies: evolution – quasars: supermassive black holes – galaxies: ISM – X-rays: galaxies – X-rays: ISM
slugcomment: Submitted, July 25, 2016

## 1. Introduction

The relationship of the QSO activity at high with the surrounding ISM of the host galaxy, including the possibility of triggering or quenching star formation, is one of the currently most debated and still unsettled topics in the field of galaxy evolution. The massive black holes (hereafter MBHs) at the centers of massive elliptical galaxies were in place already when the universe was 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, Feruglio et al. 2015, King & Pounds 2015). But the QSO and starformation activities are not likely to suddenly stop at high . Concerning the origin of the QSO activity, it appears that the most luminous AGN phases, preferentially found at , 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 AGNs seem to be driven by other processes, unrelated with the merging phenomenon, in general (Cisternas et al. 2011, Treister et al. 2012, Schawinski 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 the QSO activity also for isolated early-type galaxies (hereafter ETGs), that are often considered “dead”. Indeed, these losses represent a major source in mass for the ISM (e.g., Ciotti et al. 1991); also, the high metallicity of the observed outflows from low and high- 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 provided clear support for the conjecture that they trigger the QSO activity at epochs closer than (Ciotti et al. 2010). The 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, the QSO activity, and the resulting feedback action on the host galaxy, since , after the MBH and the stellar population have terminated their major growing phase111For a standard cosmology 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 star formation process, and MBH formation.. 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 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 broad absorption line (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 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 considered the corresponding cooling and heating functions, including photoionization plus Compton scattering, and solved the radiative transport equations, also in presence of dust; they allowed for mass and energy inputs from stellar winds, and Type Ia and Type II supernovae (hereafter respectively SNIa and SNII); and they considered the mechanical feedback due to the nuclear wind and star formation (hereafter SF) induced by accretion.

These simulations covered length scales from pc to kpc, and timescales from yr (or less) to yr; thus, all the relevant length and time-scales 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 all the required spatial extent and resolution, or to difficulties in resolving the gas mass distribution intrinsic to the numerical modeling (e.g., in the 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. 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 the medium-high mass ETGs the resulting evolution is highly unsteady (e.g., Ciotti, Ostriker & Proga 2010). At early times (starting from ) 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 the MBH heating alone has been shown to be not sufficient by itself to avoid long-lasting and massive inflows towards the galactic center at early times, but when coupled with the SNIa’s heating, it becomes very efficient in sustaining the 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 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.11 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: 1) 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 NFW (Navarro et al. 1997) dark matter 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. 2) 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 time scales, 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 ? Is it plausible that the relation was already in place at ? 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 a follows: Sect. 2 describes the galaxy models; Sect. 3 presents the numerical code and the hydrodynamical equations it solves, the inputs to them, the implementation of the feedback and SF physics; Sect. 4 presents the results of the simulations; Sect. 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 (Sect. 2.1), and their internal dynamics (Sect. 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, Parriott & Bregman 2008, Bregman & Parriott (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 Sect. 2.2 below. Far infrared observations allow for measurements of the stellar mass loss rate for the whole galaxy (), 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 with time, for solar metal abundance, after an age of Gyr, can be approximated as:

 ˙M⋆(t)=10−12A×M⋆t−1.312(M⊙yr−1), (1)

where is the galactic stellar mass in solar masses at an age of 12 Gyr, is the age in units of 12 Gyrs, and 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’s explosions provide mass and heat to the ISM, and the total mass loss rate of a stellar population is , where the mass input due to SNIa’s is . Here (in yr) is the evolution of the explosion rate with time, and each SNIa ejects . In models of SNIa’s explosions past a burst of SF (Greggio 2010), experiences a raising epoch during the first 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

 RSN(t)=0.16(H0/70)2×10−12LBt−s12(yr−1), (2)

where is the Hubble constant in units of km s Mpc, is the present epoch B-band galaxy luminosity in , and characterizes the secular evolution; when , Eq. 2 gives the rate for local ETGs in recent measurements (e.g. Mannucci et al. 2005, Maoz et al. 2011). For the rate in Eq. 2 and km s Mpc, one obtains yr, that is almost times smaller than the ”quiescent” stellar mass loss rate yr given above. Recent estimates of the slope agree with a value around (Maoz et al. 2011, Sharon et al. 2010).

The heating rate provided by SNIa’s explosions is the product of the kinetic energy injected by one event ( erg) times the rate , and times an efficiency factor. In the code we adopt an efficiency of 0.85, an assumption that is not unreasonable for the hot diluted gas (see also Sect. 3.1). Then at most , that is:

 LSN(t)=5.1(H0/70)2×1030LBt−s12(ergs−1). (3)

Of course, another major source of ISM heating is provided by the central MBH, as will be discussed in Sect. 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 ), there is a stellar component characterized by different intrinsic flattenings, and a spherical dark matter (DM) halo. The stellar density is described by the ellipsoidal deprojection (Mellier & Mathez 1987) of the de Vaucouleurs (1948) law:

 ρ⋆(R,z)=ρ0ζ−0.855exp(−ζ1/4), (4)

with

 ρ0=M⋆b1216πqR3e0Γ(8.58),ζ=b4Re0√R2+z2q2, (5)

where are the cylindrical coordinates, , is the projected half mass radius (effective radius) when the galaxy is seen face-on222For an edge-on view, the circularized effective radius is ., and the parameter controls the flattening, so that the minor axis is aligned with the axis. For the simulations we consider values of (1, 0.6, 0.3), corresponding to E0, E4 and E7 galaxies when seen edge-on. For the DM halo we adopt an untruncated NFW (Navarro et al. 1997) profile:

 ρh(r)=ρcrit δcrhr(1+r/rh)2, (6)

where is the critical density for closure, and

 δc=2003c3ln(1+c)−c/(1+c),c≡r200rh (7)

and is the radius of a sphere of mean interior density of 200. We refer to the DM mass enclosed within as to the halo mass .

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 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 due to all components (stars, dark halo, MBH), and the vertical and radial forces. The radial and vertical velocity dispersions are equal (), and the only non-zero streaming motion is in the azimuthal direction (). To set the latter, we adopted the Satoh (1980) -decomposition , from which the azimuthal velocity dispersion is recovered as , where . For the galaxy is an isotropic rotator, while for no net rotation is present, and all the flattening is due to . In general, can be a function of , and more complicated (realistic) velocity fields can be realized (Ciotti & Pellegrini 1996; Negri et al. 2014a,b). In any case, is bounded from above by the function . In this work we restrict to the case, and then in practice only stellar random motions are thermalized (see also Sect. 4 for more discussion on the role and then the adopted values of the parameter).

The only free parameter of the stellar distribution is , the aperture luminosity-weighted velocity dispersion within , from which the galaxy luminosity is recovered from the Faber–Jackson relation, and the size 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, the stellar mass is derived (see Posacki et al. 2013 for more details). By assigning a given , 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 , by fixing (Behroozi et al. 2013), and imposing that the DM fraction within a sphere of radius keeps well below unity (Cappellari et al. 2015). These constraints produce , , and for the spherical progenitors. The flattened descendants of each progenitor have the same circularized as the spherical progenitor, when seen edge-on, thus they expand with decreasing (); their , as a consequence, also decreases with respect to that of the progenitor (while and 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 description333The 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 as the spherical progenitor when seen face-on.).

We consider here four spherical progenitors, with values of 180, 210, 250 and 300 , and corresponding stellar masses of 0.81, 1.54, 3.35 and 7.80. The simulations are run for eight flattened descendants: for each progenitor, they have intrinsic flattenings with (E4 shape) and (E7 shape). Table LABEL:tab:params lists all the relevant parameters characterizing these eight galaxy models.

## 3. Feedback and star formation in the hydrodynamical equations

We describe here the hydrodynamical equations that are solved to evolve the gas flow (Sect. 3.1), and the input physics for them, with particular regard to the implementation of the SF process (Sect. 3.2) and feedback effects (Sects. 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:

 ∂ρ∂t+∇⋅(ρu)=˙ρIa+˙ρ⋆+˙ρII−˙ρSF+˙ρw, (8) ρ∂u∂t+ρ(u⋅∇)u=−∇p−ρ∇Φtot−∇prad++(˙ρIa+˙ρ⋆+˙ρII)(v⋆−u)+˙ρw(vw−u), (9) ∂E∂t+∇⋅(Eu)=−p∇⋅u+H−C+˙E−˙ESF+˙ρw2∥vw−u∥2 (10)

where , , , are respectively the mass density, velocity, internal energy density, and pressure of the gas; is the streaming velocity of the stellar component; is the AGN wind velocity. is the volumetric heating rate due to radiative feedback, and is the bolometric cooling rate per unit volume (see Sect. 3.4). The ISM is a fully ionized monoatomic gas, with , , and solar composition (). The mass density rates and describe the mass injection from the old stellar population, i.e., from SNIa’s explosions and from normal stars (Sect. 2.1); and describe the mass sink due to SF, and the mass input due to the SNII produced by SF (see Sect. 3.2); is the energy sink due to SF (Sect. 3.2). 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 ( and respectively), and by the thermalization of the relative motions between stars and the ISM (existing at the moment of injection). is then given by:

 ˙E=˙EIa+˙EII+˙ρIa+˙ρ⋆+˙ρII2[∥v⋆−u∥2+Tr(σ2)], (11)

where , with erg and being respectively the kinetic energy and ejected mass of one SNIa event, and 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); is the SNII energy injection rate, and is calculated considering the overlap in time of subsequent SF episodes (as for below; see Sect. 3.2, and Negri et al. 2015). also depends on the thermalization efficiency, that we take equal to . Finally, is the trace of the stellar velocity dispersion tensor.

The other symbols in eqs. (8)–(10) describe the mass, energy and momentum source terms due to radiative and mechanical feedback, and are described in Sects. 3.3 and 3.4 below. In particular, the terms and are given in eqs. (25)–(28). 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

Star formation, and the 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 a SF rate per unit volume given by:

 ˙ρSF=ηSFρτSF, (12)

where is the local gas density, is the SF efficiency, for which we adopt the value of 0.1, and

 τSF=max(τcool,τdyn), (13)

where

 τcool=EC,τdyn=min(τjeans,τrot), (14)

with

 τjeans=√3π32Gρ,τrot=2πrvc(r). (15)

Here is the Newtonian gravitational constant, is the distance from the galaxy center, is an estimate of the radial epicyclic period, and is the galaxy circular velocity in the equatorial plane. The energy and momentum sinks associated with star formation are then:

 ˙ESF=ηSFEτSF,˙mSF=ηSFmτSF=˙ρSFu, (16)

where and 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 at each time comes from considering that a given SF episode generates SNII’s that inject mass (at a rate exponentially declining on a timescale yr), and that during the evolution of that episode other episodes may take place, forming other SNII’s that in turn eject mass into the ISM. The same considerations are taken into account to compute the SNII energy injection rate (see Negri et al. 2015 for more details on how and 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:

 ˙MBH=˙Min1+η, (17)
 ˙Mout=η˙MBH, (18)
 Lw=ϵw˙MBHc2, (19)
 ˙pw=˙Moutvw, (20)

where is the mass accretion rate on the MBH, is the mass outflow rate in the conical wind, is the mass inflow rate at the first active radial grid, is the efficiency of generating mechanical energy with an AGN wind, is the modulus of the AGN wind velocity, and . 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 starforming disk.

Equations (17)â(20) above are adopted by most authors to treat AGN feedback as a process comprising both infall and outflow; however, typically is adopted, implicitly assuming , so that and are neglected, and and may be overestimated. For example, if we adopt as many authors [e.g., Springel et al. (2005), Johansson et al. (2009), McCarthy et al. (2010)], and (e.g., Moe et al. 2009), then , and all the normally neglected effects may in fact be dominant: the bulk of the inflowing mass () may be ejected in a broad-line disk wind, and the effects of the mass () and momentum () input deposited in the ambient gas may dominate over the energy input (), which may be largely radiated away. Thus, it is important to consider consistently the effects of including mass, energy, and momentum conservation when ; our treatment is consistentwith observations of BAL quasars and the radiative mechanisms that drive them (Ciotti et al. 2010).

In general, is fixed for given and . The wind efficiency, , is not known very well, neither from observations nor from detailed simulations. The best estimates might be in the range (Proga et al. 2000; Proga & Kallman 2004; Krongold et al. 2007; Kurosawa et al. 2009). Here we consider models where and 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):

 ϵw=ϵw0Aw˙m1+Aw˙m,vw=vw0Aw˙m1+Aw˙m, (21)

where , is the dimensionless mass accretion rate, i.e., normalized to the Eddington mass accretion rate (see below), and the two functions and saturate to and , respectively, for large accretion rates. We choose the constant values and  km/s (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, Tombesi et al. 2015, Feruglio et al. 2015). Independently of the mechanical feedback model, the radiative luminosity of the AGN is given by:

 LBH=ϵEM˙MBHc2, (22)

where the electromagnetic efficiency 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):

 ϵEM=ϵ0AEM˙m1+AEM˙m, (23)

and and . The dimensionless mass accretion rate is

 ˙m≡˙MBH˙MEdd=ϵ0˙MBHc2LEdd, (24)

where is the Eddington luminosity. With the settings above, has a minimum value of 0.18 at the highest accretion rates, and increases for decreasing (for example, for ).

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 eqs. (18)-(20), for a given MBH accretion rate ; then, the specification of the angular dependence of the properties of the AGN conical wind (, and ) is required. For this dependence we adopt:

 ˙ρw=˙Moutδ(r−r0)r2sinθf(θ), (25)

where is the Dirac delta-function, is the first gridpoint, is the angle from the -axis, measured clockwise, and:

 f(θ)=(n+1)sinθ|cosθ|n4π. (26)

The anisotropic momentum source in eq. (9) is given by:

 ˙mw=˙ρwvwer=˙pwδ(r−r0)r2sinθf(θ)er (27)

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

 ˙Ew=12˙ρwv2w=Lwδ(r−r0)r2sinθf(θ). (28)

Note that is independent of , and the anisotropy in the momentum and energy injection is due to the anisotropic mass injection [eq. (26)]. Following Novak et al. (2011), we adopt , so that the half-opening angle enclosing half of the energy is . In terms of solid angle, this means that the wind is visible from of the available viewing angles. This fraction is in agreement 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. 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 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  K is given by:

 H−C≡n2(S1+S2+S3), (29)

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

 S1=−3.8×10−27√T. (30)

The Compton heating and cooling is given by

 S2=4.1×10−35(Tc−T)ξ, (31)

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

 S3=10−23a+b(ξ/ξ0)c1+(ξ/ξ0)cZZ\sun, (32)

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

 a=−18e25(logT−4.35)2−80e5.5(logT−5.2)2−17e3.6(logT−6.5)2, (33)
 b=1.7×104T−0.7,c=1.1−1.1eT/1.8×105+4×1015T4, (34)
 ξ0=(1.5T0.5+1.5×1012T2.5)−1+4×1010T2[1+80e(T−104)/1.5×103]. (35)

Gas temperatures are bounded from below by the adopted atomic cooling curve, that has an exponential cutoff below K.

In this work we do not consider the effects of dust, i.e., the radiation momentum associated with dust absorption (cf. Debuhr et al. 2011; Bieri et al. 2016); however, the effect of dust has been shown not to alter significantly the gas evolution (Hensley et al. 2014). 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 eq. 55 in Ciotti & Ostriker 2007). The integration scheme used here for the radiation transport is not the full scheme of Novak et al. (2012, their Sect. 3), but the simplified version in their Appendix B, in order to speed up the simulations. For example, in the case of spherical symmetry, the equation to be integrated would be:

 dLeffBH,photo(r)dr=−4πr2H, (36)

where is the effective accretion luminosity at , and the equation is solved with central boundary condition given by equation (22). The force per unit mass due to photoionization+Compton opacity can be expressed as

where

 κphoto=−1ρ(r)LeffBH,photo(r)dLeffBH,photo(r)dr==4πr2H(r)ρ(r)LeffBH,photo(r). (38)

The radiation pressure due to electron scattering is

where .

Finally, equations (31) and (32) depend on the ionization parameter, that is given by:

 ξ≡LeffBH,photo(r)n(r)r2. (40)

## 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 of meshpoints, spanning from 2.5 pc to 250 kpc. Reflecting boundary conditions were set along the -axis, while at the outer edge of the computational domain the fluid is free to flow out.

Negri et al. (2014a,b) 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 to the low-rotation case. However, in axysimmetric systems some rotation is numerically needed to prevent gas from unphysically sticking on the -axis. Different recipes have been figured out to solve this problem (Novak et al. 2011; Li, Ostriker & Sunyaev 2013); here, we take advantage of the Jeans solver that allows us to tune the Satoh -parameter (Sect. 2.2). In detail, we determine in order to have negligible but non-zero ordered rotation over the main galaxy body, with angular momentum of the stars never exceeding , 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 -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,θ)≡ηrotJ0√J20+R2v2IS;v2IS≡¯¯¯¯¯v2φ−σ2, (41)

where is the distance from the -axis, is the ordered stellar velocity if the galaxy model were an isotropic rotator, and we fix .

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

 LX=∫εXdV,TX=1LX∫TεXdV, (42)

where 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 (), defined as the percentage of time spent at a ; the time at which half of the MBH radiation energy has been emitted (); the time at which half of the new stellar mass has been created (); and the half-mass radius of the new stars formed until the end of the simulation ().

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 star formation 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 Tab. 1. For each galaxy model we ran three simulations: one without any feedback from the MBH (we refer to these models as to 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). The most interesting quantities at the end of the twenty-four simulations are listed in Tab. 2. Rather than presenting the specific evolution of the gas flow in each one 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 the hot gas (Sect. 5.1), the MBH growth (Sect. 5.2), the newly formed stars (Sect. 5.3), and the duty cycle (Sect. 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 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 the related quantities (as discussed below). At fixed galaxy mass, a change in the galaxy shape produces a lower for flatter morphologies, because the outflowing region becomes larger, as obtained for 2D models without feedback (Negri et al. 2014b).

We start considering a representative E4 model of average mass (E4 in Tab. 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 galaxy: and calculated within 5, the hot gas mass within the whole grid, and the SFR. The NOF model shows a smooth evolution in , , hot gas mass (the mass with K), and SFR; in particular, in each panel of Fig. 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 , , 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 Tab. 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 of feedback models, when calculated within , shows drops that reach one order of magnitude below of NOF models, i.e. that are much larger than those in Fig. 1. The drops are the natural consequence of the clearing of the gas in the central regions.

Figure 2 shows the vs. relation for all models, at the present epoch, in comparison with the observed and values recently measured for the gas only, using 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 a remarkable agreement between models and real ETGs. In particular, the large spread in at low luminosities is well reproduced, as due to mostly outflowing ETGs at lower galaxy masses. Also well reproduced is the observed average trend for erg s (note that most of the highest observed are due to central ETGs in groups, Goulding et al. 2016).

The time evolution of , , and SFR is similarly rich in large fluctuations in models more massive than the E4 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 , and SFR values (Fig. 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 , one can compute the fraction of time spent above erg s, that represents a sort of upper value for the observed gas emission of normal ETGs in the local universe444Note that this is conservative assumption; for example, central group or cluster ETGs can show much larger values. (Kim & Fabbiano 2015), over the past 3 or 5 Gyr. For both choices of lookback time, the E4 and E7 FF models spend % of the time with erg s, and lower mass models spend % of their time above it (the very massive E4 and E7 models have precisely of the order of erg s in the past few Gyr; see Tab. 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 FF model; in Figs. 37 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 (Figs. 3 and 4), where outflow and inflow regions coexist (Fig. 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 kpc. The rightmost panels show again a quiescent state, at a time of 6.95 Gyr, i.e., yr 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 (Figs. 3 and 4); the wind itself is outflowing in a small hourglass region above and below the nucleus (Fig. 5), and imparts a tangential velocity to the gas surrounding it (Fig. 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, Cresci et al. 2015, Carniani et al. 2015); this observation is considered an important evidence for negative feedback, however we find that overall the feedback effect on SF is positive (Sect. 5.3). Figure 8 shows finally the X-ray surface brightness and projected temperature maps weighted with the emission over 0.3–8 keV (Sect. 4), corresponding to the two central times during the outburst of Figs. 3–7 (t=6.85 and t=6.86 Gyr). A very bright central region of hundred of 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 ’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 variance 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 is escaping 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 on to a MBH has been intensively investigated recently with a variety of numerical simulations (e.g. Pizzolato & Soker 2010; Barai, Proga & Nagamine 2012; Nayakshin & Zubovas 2012; Gaspari, Ruszkowski & Oh 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 Fig. 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, Ostriker & Proga 2009, 2010), has interesting consequences for the cumulative SF of MF and FF models (Sect. 5.3).

### 5.2. The black hole mass growth

Figure 9 shows the adopted initial relation (Sect. 2.2), the relation derived from dynamical studies of well observed local ETGs (McConnell & Ma 2013), and the values at the end of the simulations (when the MBH mass has grown, while remains constant); different symbols distinguish the galaxy shapes and types of feedback. The first thing to notice is that models with feedback produce a final relation that agrees well with the observed one, within the uncertainties. This agreement is not fulfilled by NOF models that end with overwhelmingly large masses (not shown in the figure, but see Tab. 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 and shapes; in fact, the differences between the final values at fixed , 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 the radiative feedback alone has minor effects). Looking deeper into Fig. 9, however, one can notice some trends at fixed . E7 models show lower than E4 ones, a result due to the lower gas binding energy of flatter galaxies (Ciotti & Pellegrini 1996, Posacki et al. 2013), with the consequent larger outflow region with respect to 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 than MF ones (Fig. 9), a result of their larger 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 , and more in general what is the effect of the initial (albeit at ) value of on the flow evolution and on the final . However, simulations with smaller 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 , for feedback models, is different for the different galaxy masses: it ranges from a figure of the order of % for the lowest mass models (E), to roughly % (the MBH mass doubles) for the E models, to roughly triplicate for the E models, to almost quadruplicate for the very massive E galaxies (see col. 3 in Tab. 2, and Fig. 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 , the right panel of Fig. 9 shows the same trends shown in the left panel: the MBH mass increases more, even with respect to the initial , 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 absence of AGN feedback and SNIa’s assisted galactic winds, the mass increase of would be larger by up to two orders of magnitude (Tab. 2). It is a valuable property of the present models that the average increase in since is just of a factor of few (), as found for the average growth history of MBHs with given starting mass (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. ranges from to , from the low to the high galaxy masses. The lowest mass models show significantly lower because during their evolution the gas is mostly outflowing, they experience little accreting mass and little possibility of SF. At fixed , rounder models produce slightly larger 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 , at any fixed and galaxy shape, with larger in the FF case. The larger efficiency of FF models to produce SF is a consequence of the richer structure and longer duration of each of their nuclear outburst episodes, that prevent prompt accretion, provide the gas more time for its cooling, and lead to a larger capability to keep cold and dense gas far from the nucleus (see Sect. 5.1, and Fig. 1, where an FF outburst shows more fluctuations in , and a larger duration of the feedback action on the gas, than a MF outburst). Note that the final , at fixed and galaxy shape, increases from NOF to MF to FF models (Tab. 2, except for the two NOF E models), which proves definitely 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 growth is opposite to that of the MBH growth, with respect to feedback: increases from NOF to FF models, while the opposite is true for (cfr. Figs. 9 and 10). Both and , instead, are larger for rounder models.

It is interesting to compare the relative increase of and (values in Tab. 2). In NOF models, and are not largely different, and the MBH growth is much favored over the SF, with respect to feedback models. For feedback models, ranges from to , thus is always much larger than , at all galaxy masses. There is a trend, though: (Fig. 9) increases clearly with , while keeps quite flat or even decreases (Fig. 10, right panel). Thus, increases with , 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 : as the stellar mass does not increase significantly after (Fig. 10, right panel), and can increase up to a factor of 4 (Fig. 9, right panel), this implies that after 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 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 eq. (1), these losses amount to % of . From Fig. 10 we see that is % of (or lower for the least massive models), thus, roughly half of the integrated stellar mass losses within a galaxy since goes into new stars. The remaining part of the mass losses goes for a minor fraction into (that reaches at most , and then at most 3% of the integrated stellar mass losses, see Fig. 9), and for a major fraction is ejected from the galaxy, due to SNe heating and the further help of the AGN feedback action. We recall that the AGN feedback produces just an increase in the ejected mass from the galaxy, but has never been found capable of producing a global/major outflow by itself during the passive galaxy evolution after . To better quantify this increase, we can evaluate the amount of gas residing outside a reference radius of 5 (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 % to %, 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 ; such a factor is useful to establish the possible role of AGN feedback in adding material to the circumgalactic medium, after . For feedback models, the load factor is , except for the lowest mass (E) galaxies, that eject a larger amount of gas, and for which then the load factor is . For NOF models, the factor is slightly larger (), but not so different, because they eject less mass, but also form less stars; a markedly lower factor () than in feedback models is instead shown by the E 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 () 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 is formed ( from Tab. 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 Gyr form quite early, during the first Gyr. Figure 1 shows the time evolution of the SFR for the E4 model, and Tab. 2 gives the instantaneous value of the SFR at the end of the simulations. These present-epoch rates typically keep below yr, and compare well with the current rates observed for molecular gas-rich ETGs of the local ATLAS sample, which range from to yr, with a median value of yr (Davis et al. 2014). Larger SFR than yr are shown by models caught in an outburst (Fig. 1), or in the most massive models, that host more gas, and have a larger AGN activity; note that the ATLAS survey is not representative of the most massive models in Tab. 2.

Another interesting aspect for a 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), that experiences an almost global outflow over its whole lifetime. SF is instead evident in the larger mass E4 model: it forms a nuclear stellar disk in the NOF case555This disk is the result of the stellar streaming imposed in the central regions for numerical reasons, and discussed in Sect. 4, see eq. (38)., while it has a roughly spherical distribution, peaking within kpc radius, in the cases with feedback. The extent of the SF is measured more quantitatively by , the radius including half of the new stars at the end of the simulation (in Tab. 2). Typically , for feedback models, except for the least massive ones, where SF is very low but extended. In contrast, 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, 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 the SF, inclined by , 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. 2015).

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 Sersic 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 ˙MBH, LBH and the duty cycle

Figure 13 shows the time evolution of (eq. 18), and (eq. 23), that are relevant to quantify the accretion history of the MBH, for the representative model E4 in Fig. 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 Myr, as for more massive models, or for a shorter time, and with smaller , as for lower mass models (E).

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 , and with above/below a chosen threshold. In particular, for the E4 and E7 FF models, Fig. 14 shows the fraction of total simulation time, and the fraction of total radiative energy, respectively spent and emitted at different values of . For all models the distribution of times has a long, almost flat tail extending to the lowest , and peaks at an that increases with the galaxy mass: the peak is located at , and reaches for the most massive models (while becoming less and less pronounced). All distributions drop sharply above , so that all models spend very little time above . An average, realistic galaxy, as could be described by the E4 FF model, spends 75% of its time since below . 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 , as shown by the right column in Fig. 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 , 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 , 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 ranging from 3 to 5%, going from the E4 to the E4 models (Fig. 15, left panels, solid lines).

For reference, Table 2 gives the fraction of time spent above (); the -value with respect to which the MBH energy is emitted equally above and below (); and the 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 (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 Figs. 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 , and the duty cycle of AGN activity (e.g., Marconi et al. 2004, Trakhtenbrot & Netzer 2012, Schulze et a. 2015, Caplar et al. 2015). We stress that our predictions concern massive non-rotating ETGs, with M, since , i.e. after the 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, 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 studies. 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’s rate. Star formation is implemented via a simple scheme based on physical arguments shown to reproduce well the Kennicutt-Schmidt relation (Negri et al. 2015). The galaxy models are axisymmetric, and the Jeans equations provide the detailed internal stellar kinematics on which the stellar kinematical heating is based. The stellar density profile follows a deprojected (ellipsoidal) Sersic law, and the main observables () are chosen in order for the galaxy models to lie on the main scaling laws. The spherical dark matter halo has a radial profile predicted by cosmological simulations (NFW), and is normalized to account for a dark mass within 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 the consequent disposal of fresh gas on the central MBH. A new time-scale 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. Secondly, the radiative transport is here calculated by using the simplified (and computationally much faster) scheme described in Novak et al. (2012), that has been shown to produce very accurate results when compared to the exact numerical solution. Also not included are all 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, Novak et al. 2012). 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.

The main results of this work are as follows.

1) and 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 at the low 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 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 and , 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’s heating666In 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). 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 ( erg s), is %. After each outburst, within (1–2) is lower by 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 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 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 one hundred 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 , for feedback models, ranges from to , 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 , 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 , is larger for the full feedback case than the mechanical one, that in turn produces a larger 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 2007; e.g. see also Nayakshin & Zubovas 2012, Ishibashi & Fabian 2012); thus geometrical effects allowed for by 2D simulations and a flat galaxy shape do not invalidate the global picture.

5) is of the order of 0.040.05, excluding the lowest mass models that show lower values (), 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 Gyr ends recycled into new stars. The other half goes for a minor fraction into (for feedback models), and for a major fraction is ejected from the galaxy, mostly due to SNe heating. AGN feedback produces just an increase (of the order of %) in the ejected mass from the galaxy, during the passive galaxy evolution after . The “load factor”, defined as the ratio between the gas displaced outer of 5 during the evolution, and the mass in new stars , is for feedback models, and slightly larger for NOF models (due to their lower ability to form new stars). The lowest mass (E) galaxies eject a larger amount of gas, and for them the load factor is larger ( for models with AGN feedback).

6) Final values of feedback models are consistent with those measured in the local universe, at all considered here, while they are overwhelmingly large without feedback; the latter is then fundamental to keep the MBH masses at the values observed, after . In feedback models just % of the total injected mass from stars is accreted on the MBH; the average increase in since is just of a factor of few (). 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 , flatter models show lower , due to the lower gas binding energy, and FF models end with lower than MF ones, thanks to their larger capability of preventing gas accretion. ranges from to , from the lowest to the largest galaxy masses. The implementation of rotation could change this trend, by making the MBH growth less favored, and SF in the cold disk more favored.

7) Half of the radiative output from the AGN is emitted within Gyr since the birth of the galaxy, roughly the same value of the timescale 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: precedes in lower mass galaxies, while the reverse is true in larger mass ones.

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 ), but the energy is emitted at high Eddington ratios, with one quarter of the total energy emitted at . The resulting duty-cycles of AGN activity, estimated as the fraction of time spent above , ranges 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. LC, SP and AN 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.

## References

• Alvarez (2009) Alvarez M.A., Wise J.H., Abel T. 2009, ApJ, 701, L133
• Arav (2013) Arav, N., Borguet, B., Chamberlain, C., et al. 2013, MNRAS, 436, 3286
• Athey (2002) Athey, A., Bregman, J., Bregman, J., Temi, P., Sauvage, M. 2002, ApJ, 571, 272
• Bae (2016) Bae, H.-J., & Woo, J.-H. 2016, ApJ, in press (arXiv:1606.05348)
• Barai (2012) Barai, P., Proga, D., Nagamine, K. 2012, MNRAS, 424, 728
• Behroozi (2013) Behroozi, P. S., Wechsler, R. H., Conroy C. 2013, ApJ, 770, 57
• Bertin (2001) Bertin G., Lodato G. 2001, A&A, 370, 342
• Bieri (2016) Bieri, R., Dubois, Y., Rosdahl, J., et al. 2016,MNRAS, in press (arXiv:1606.06281)
• Bondi (1952) Bondi H. 1952, MNRAS, 112, 195
• Bregman (2009) Bregman, J.N., Parriott, J.R. 2009, ApJ, 699, 923
• Booth (2009) Booth, C. M., Schaye, J. 2009, MNRAS, 398, 53
• Cano (2012) Cano-Díaz M., Maiolino R., Marconi A., et al. 2012, A&A 537, L8
• Caplar (2015) Caplar N., Lilly S.J., Trakhtenbrot B. 2015, ApJ811, 148
• Cappellari (2015) Cappellari, M., Romanowsky, A. J., Brodie, Jean P., et al. 2015, ApJ, 804, L21
• Cappellaro (1999) Cappellaro, E., Evans, R., Turatto, M. 1999, A&A, 351, 459
• Carniani (2015) Carniani S., Marconi A., Maiolino R., et al. 2015 A&A, 580, 102
• Cheung (2016) Cheung E., Bundy K., Cappellari M., et al. 2016, Nature, 533, 504
• Choi (2015) Choi, E., Ostriker, J.P., Naab, T., Oser, L., Moster, B. P. 2015, MNRAS, 449, 4105
• Ciotti (2009) Ciotti L. 2009, Nature, 460, 333
• Ciotti (1991) Ciotti, L., DâErcole, A., Pellegrini, S., Renzini, A. 1991, ApJ, 376, 380
• Ciotti (1996) Ciotti, L., Pellegrini, S. 1996, MNRAS, 279, 240
• Ciotti (1997) Ciotti, L., Ostriker, J.P. 1997, ApJ, 487, L105
• Ciotti (2001) Ciotti, L., Ostriker, J.P. 2001, ApJ, 551, 131
• Ciotti (2007) Ciotti, L., Ostriker, J.P. 2007, ApJ, 665, 1038
• Ciotti (2009) Ciotti, L., Ostriker, J.P., Proga, D. 2009, ApJ, 699, 89
• Ciotti (2010) Ciotti, L., Ostriker, J.P., Proga, D. 2010, ApJ, 717, 708
• Ciotti (2012) Ciotti, L., Ostriker, J.P. 2012, in Hot Interstellar Matter in Elliptical Galaxies, Kim D.-W., Pellegrini S., eds, Astrophysics and Space Science Library, Vol. 378. Springer-Verlag, Berlin, p. 8
• Cisternas (2011) Cisternas, M., Jahnke, K., Inskip, K., et al. 2011, ApJ, 726, 57
• Cooksey (2010) Cooksey, K.L., Thom, C., Prochaska, J.X., Chen, H-W. 2010, ApJ, 708, 868
• Cresci (2015) Cresci, G. et al. 2015, A&A, 582, A63
• Curtis (2016) Curtis, M., Sijacki, D. 2016, MNRAS, submitted (arXiv:1606.02729)
• Davis (2014) Davis, T.A., Young, L.M., Crocker, A.F., et al. 2014, MNRAS, 444, 3427
• Debuhr (2011) Debuhr, J., Quataert, E., Ma, C.-P. 2011, MNRAS, 412, 1341
• Desroches (2007) Desroches, L.-B., Quataert, E., Ma, C.-P., West, A. A. 2007, MNRAS, 377, 402
• deVaucouleurs (1948) de Vaucouleurs, G. 1948, Annales dâAstrophysique, 11, 247
• Dimatteo (2005) Di Matteo, T., Springel, V., Hernquist, L. 2005, Nature, 433, 604
• Dimatteo (2008) Di Matteo, T., Colberg, J., Springel, V., Hernquist, L., Sijacki, D. 2008, ApJ, 676, 33
• Dubois (2012) Dubois,Y., Devriendt, J., Slyz, A., Teyssier, R. 2012, MNRAS, 420, 2662
• Faucher (2012) Faucher-GiguÃ¨re, C., Quataert, E. 2012, MNRAS425, 605
• Feruglio (2015) Feruglio, C., Fiore, F., Carniani, S., et al. 2015, A&A, 583, A99
• Fischer (2013) Fischer, T.C., Crenshaw, D. M., Kraemer, S.B., Schmitt, H.R. 2013, ApJS, 209, 1
• Fox (2011) Fox, A.J. 2011, ApJ, 730, 58
• Gaspari (2013) Gaspari, M., Ruszkowski, M., Oh, S.P. 2013, MNRAS, 432, 3401
• Gibson (2009) Gibson, R.R., Brandt, W.N., Gallagher, S. C., Schneider, D.P. 2009, ApJ, 696, 924
• Goncalves (2008) Goncalves, T.S., Steidel, C.C., Pettini, M. 2008, ApJ, 676, 816
• Greene (2011) Greene, J.E., Zakamska N.L., Ho L.C., Barth, A.J. 2011, ApJ732, 9
• Greggio (2010) Greggio, L. 2010, MNRAS, 406, 22
• Goulding (2016) Goulding, A.D., Greene, J.E., Ma, C.-P., et al. 2016, in press on ApJ(arXiv:1604.01764)
• Hayes (2006) Hayes, J.C., Norman, M.L., Fiedler, R.A., et al. 2006, ApjS, 165, 188
• Heckman (2014) Heckman, T.M., Best, P.N. 2014, ARA&A, 52, 589
• Hensley (2014) Hensley, B.S., Ostriker, J.P., Ciotti, L. 2014, ApJ, 789, 78
• Hopkins (2009) Hopkins, P.F., Cox, T.J., Dutta, S.N., et al. 2009, ApJS, 181, 135
• Hopkins (2011) Hopkins, P.F., Quataert E. 2011, MNRAS, 415, 1027
• Ishibashi (2012) Ishibashi, W., Fabian, A.C. 2012, MNRAS, 427, 2998
• Johansson (2009) Johansson, P.H., Naab, T., Burkert, A. 2009, ApJ, 690, 802
• Kazantzidis (2005) Kazantzidis, S., Mayer, L., Colpi, M., et al. 2005, ApJ, 623, L67
• Kim (2015) Kim, D.-W., Fabbiano, G., 2015, ApJ, 812, 127
• King (2003) King, A. 2003, ApJ, 596, L27
• King (2015) King, A., Pounds, K. 2015, ARA&A, 53, 115
• Kocevski (2012) Kocevski, D.D., Faber, S.M., Mozena, M., et al. 2012, ApJ, 744, 148
• Kormendy (2009) Kormendy, J., Fisher, D.B., Cornell, M.E., Bender, R. 2009, ApJS, 182, 216
• Korol (2016) Korol, V., Ciotti, L., Pellegrini, S. 2016, MNRAS, 460, 1188
• Krongold (2007) Krongold, Y., Nicastro, F., Elvis, M., et al. 2007, ApJ, 659, 1022
• Kurosawa (2009) Kurosawa, R., Proga, D., Nagamine, K. 2009, ApJ, 707, 823
• Lehner (2013) Lehner, N., Howk, J.C., Tripp, T.M., et al. 2013, ApJ, 770, 138
• Li (2013) Li, J., Ostriker, J.P., Sunyaev, R. 2013, ApJ, 767, 105
• Liu (2013) Liu, G., Zakamska, N. L., Greene, J. E., Nesvadba, N. P. H., Liu, X. 2013, MNRAS, 436, 2576
• Liu (2015) Liu, G., Arav, N., Rupke, D.N. 2015, ApJS, 221, 9
• Magorrian (1998) Magorrian, J. et al., 1998, AJ, 115, 2285
• Mannucci (2005) Mannucci, F., Della Valle, M., Panagia, N., Cappellaro, E., Cresci, G., Maiolino, R., Petrosian A., Turatto M. 2005, A&A, 433, 807
• Maoz (2011) Maoz, D., Mannucci, F., Li W., Filippenko, A.V., Della Valle, M., Panagia, N. 2011, MNRAS, 412, 1508
• Maraston (2005) Maraston, C. 2005, MNRAS, 362, 799
• Marconi (2004) Marconi, A., Risaliti, G., Gilli, R., Hunt, L.K., Maiolino, R., Salvati, M. 2004, MNRAS, 351, 169
• Mathews (1989) Mathews, W.G. 1989, AJ, 97, 42
• Mathews (1990) Mathews, W.G. 1990, ApJ, 354, 468
• Mccarth (2010) McCarth, I.G., Schaye, J., Ponman, T.J., et al. 2010, MNRAS, 406, 822
• Mcconnell (2013) McConnell, N.J., Ma, C.-P. 2013, ApJ, 764, 184
• Mcelory (2015) McElroy, R., Croom, S.M., Pracy, M., et al. 2015, MNRAS, 446, 2186
• Mellier (1987) Mellier, Y., Mathez, G. 1987, A&A, 175, 1
• Moe (2009) Moe, M., Arav, N., Bautista, M.A., Korista, K.T. 2009, ApJ, 706, 525
• Narayan (1995) Narayan, R., Yi, I. 1995, ApJ, 452, 710
• Navarro (1997) Navarro, J.F., Frenk, C.S., White, S.D.M. 1997, ApJ, 490, 493
• Nayakshin (2012) Nayakshin, S., Zubovas, K. 2012, MNRAS, 427, 372
• Negri (2014a) Negri, A., Ciotti, L., Pellegrini, S. 2014a, MNRAS, 439, 823
• Negri (2014b) Negri, A., Posacki, S., Pellegrini, S., Ciotti, L. 2014b, MNRAS, 445, 1351
• Negri (2015) Negri, A., Pellegrini, S., Ciotti, L. 2015, MNRAS, 451, 1212
• Ness (2016) Ness, M., Lang, D. 2016, AJ, 152, 1
• Norman (1988) Norman, C., Scoville, N. 1988, ApJ, 332, 124
• Novak (2011) Novak, G.S., Ostriker, J.P., Ciottik L. 2011, ApJ, 737, 26
• Novak (2012) Novak, G.S., Ostriker, J.P., Ciotti, L. 2012, MNRAS, 427, 2743
• Obied (2016) Obied, G., Zakamska, N., Wylezalek, D., Liu, G. 2016, MNRAS, 456, 2861
• Ostriker (2010) Ostriker, J.P., Choi, E., Ciotti, L., Novak, G.S., Proga, D. 2010, ApJ, 722, 642
• Parriott (2008) Parriott, J.R., Bregman, J.N. 2008, ApJ, 681, 1215
• Pellegrini (2005) Pellegrini, S. 2005, ApJ, 624, 155
• Pellegrini (2012a) Pellegrini, S. 2012, in Hot Interstellar Matter in Elliptical Galaxies, Kim D.-W., Pellegrini S., eds, Astrophysics and Space Science Library, Vol. 378. Springer-Verlag, Berlin, p. 21
• Pellegrini (2012b) Pellegrini, S., Ciotti, L., Ostriker, J.P. 2012, ApJ, 744, 21
• Pizzolato (2010) Pizzolato, F., Soker, N. 2010, MNRAS, 408, 961
• Posacki (2013) Posacki, S., Pellegrini, S., Ciotti, L. 2013, MNRAS, 433, 2259
• Proga (2000) Proga, D., Stone, J.M., Kallman, T.R. 2000, ApJ, 543, 686
• Proga (2004) Proga, D., Kallman, T.R. 2004, ApJ, 616, 688
• Reichard (2003) Reichard, T.A., Richards, G.T., Hall, P.B., et al. 2003, AJ, 126, 2594
• Satoh (1980) Satoh, C. 1980, PASJ, 32, 41
• Sazonov (2005) Sazonov, S.Y., Ostriker, J.P., Ciotti, L., & Sunyaev, R.A., 2005, MNRAS, 358, 168
• Sazonov (2008) Sazonov, S., Krivonos, R., Revnivtsev, M., Churazov, E., Sunyaev, R., 2008, A&A 482, 517
• Schroetter (2016) Schroetter, I., et al. 2016, submitted to ApJ(arXiv:1605.03412)
• Schulze (2014) Schulze, A., Wisotzki, L. 2014, MNRAS438, 3422
• Schulze (2015) Schulze, A., Bongiorno, A., Gavignaud, I., et al. 2015, MNRAS, 447, 2085
• Sharon (2010) Sharon, K., Gal-Yam, A., Maoz, D., et al. 2010, ApJ, 718, 876
• Sijacki (2015) Sijacki, D., Vogelsberger, M., Genel, S., et al. 2015, MNRAS, 452, 575
• Silk (1998) Silk, J., Rees, M. J. 1998, A&A, 331, L1
• Springel (2005) Springel, V., Di Matteo, T., Hernquist, L. 2005, MNRAS, 361, 776
• Sturm (2011) Sturm, E., González-Alfonso, E., Veilleux, S., et al. 2011, ApJ, 733, L16
• Tang (2005) Tang, S., Wang, Q.D., 2005, ApJ, 628, 205
• Tombesi (2015) Tombesi, F., Meléndez, M., Veilleux, S., Reeves, J.N., González-Alfonso, E., Reynolds, C.S., 2015, Nature, 519, 436
• Trakhtenbrot (2012) Trakhtenbrot, B., Netzer, H. 2012, MNRAS, 427, 3081
• Treister (2012) Treister, E., Schawinski, K., Urry, C.M., Simmons, B.D. 2012, ApJ, 758, L39
• Wu (2015) Wu, X., Wang, F., Fan, X., et al. 2015, Nature, 518, 512
• Zakamska (2016) Zakamska, N.L., Hamann, F., Pâris, I., et al. 2016, MNRAS, 459, 3144