X-RAY PROPERTIES EXPECTED FROM ACTIVE GALACTIC NUCLEUS FEEDBACK IN ELLIPTICAL GALAXIES

Supermassive black holes (SMBHs) at the centers of bulges and elliptical galaxies play an important role in the processes of galaxy formation and evolution (e.g., Cattaneo et al. 2009), as testified by remarkable correlations between host galaxy properties and the SMBH masses (e.g., Magorrian et al. 1998; Ferrarese & Merritt 2000; Gebhardt et al. 2000; Graham et al. 2001; Kormendy et al. 2009) and as supported by many theoretical studies (e.g., Silk & Rees 1998; Haiman et al. 2004; Merloni et al. 2004; Sazonov et al. 2005; Di Matteo et al. 2005; Hopkins et al. 2006; Somerville et al. 2008; Lusso & Ciotti 2010). An important aspect of the coevolution process is the radiative and mechanical feedback by the accreting SMBH onto the galactic interstellar medium (ISM) that is continuously replenished by normal stellar mass losses, at a rate of the order of ≈1 M_☉ yr⁻¹ in a medium-mass galaxy. In the absence of some form of heating (and stripping from the intracluster medium in the case of cluster galaxies), this ISM would develop a flow directed toward the galactic center, accreting ≳ 1 M_☉ yr⁻¹ in a process similar to a "cooling flow." Instead, in low-mass galaxies, type Ia supernova (SNIa) heating is able to sustain a low-luminosity, global galactic wind (e.g., Ciotti et al. 1991; David et al. 1991; Pellegrini & Ciotti 1998), and the central SMBH is in a state of permanent, highly sub-Eddington hot accretion (Ciotti & Ostriker 2012).

Therefore, in medium- to high-mass galaxies, heating is required by the following empirical arguments: (1) the large amount of gas lost by the passively evolving stellar population during the galaxies' lifetime is not observed (e.g., Peterson & Fabian 2006), and just ≲1% of the mass made available by stars is contained in the masses of present epoch SMBHs; (2) bright active galactic nuclei (AGNs), as would be expected given the predicted mass accretion rate, are not commonly seen in the spheroids of the local universe (e.g., Fabian & Canizares 1988; Pellegrini 2005). Since over a large part of the galaxies' extent, and for a large fraction of their lifetime, the ISM cooling time t_cool is much lower than the galaxy age, one early quasar phase cannot be the solution to the cooling flow problem. The solution requires either steady heating or heating with bursts on a timescale Δt ≈ t_cool. The AGN feedback is not the only heating mechanism acting on the ISM of elliptical galaxies (in addition to the thermalization of stellar mass losses and the supernova heating); thermal conduction and cosmic rays may also be present. These latter two forms of heating have been mainly proposed to halt cooling flows on cluster scales (e.g., Parrish et al. 2009; Guo & Oh 2008), but appear not as crucially relevant on galactic scales. For example, the hot coronae of the two D galaxies in the Coma cluster require highly attenuated thermal conduction (Vikhlinin et al. 2001); and a large survey of early-type galaxies in hot clusters shows that, to avoid evaporation of the retained gas after stripping, thermal conduction must be suppressed by at least two orders of magnitude compared to the Spitzer value (Sun et al. 2007). On the other hand, the energy in cosmic rays (not from the AGN) cannot exceed that from supernovae, as these are believed to be their main source. Therefore, the AGN feedback heating appears a most natural solution. Unfortunately, on a purely theoretical ground, how much radiative and mechanical energy and momentum output from the SMBH can effectively interact with the surrounding ISM, and what the resultant SMBH masses are, is difficult to establish.

Recently, the interaction of the SMBH output with the inflowing gas has been studied with high-resolution one-dimensional (1D) hydrodynamical simulations in a series of papers (Ciotti & Ostriker 2007; Ciotti et al. 2009b, 2010, hereafter Papers I and III; Shin et al. 2010a, 2010b; Ostriker et al. 2010; Jiang et al. 2010), that are currently being extended to two-dimensional (2D) treatments (Novak et al. 2011). These simulations implement a physically based detailed treatment of the radiative energy and momentum input from the SMBH into the ISM, consistent both with observed average AGN spectra and theoretical calculations of radiation transport; they also include star formation and a modeling of the mechanical energy and momentum feedback from AGN winds. We refer the reader to the papers listed above for a detailed description of the feedback physics and effects on the ISM (see also Ciotti & Ostriker 2012); here we recall that the combined effects of radiative and mechanical feedback produce recurrent AGN burst phases accompanied by star formation in cold shells of ISM on the kiloparsec scale, temporally spaced apart by longer phases of relative quiescence. A cycle repeats with the galaxy seen alternately as an AGN/starburst for a small fraction of the time and as a "normal" elliptical for much longer intervals. Accretion-fueled feedback thus proves effective in suppressing long lasting cooling flows and in maintaining SMBH masses within the range observed today, since the gas is mostly lost in outflows or consumed in starbursts. Note that a major role in producing global degassing, and in regulating the flow evolution, is also played by the SNIa's heating. Remarkably, our feedback models show a two-stage regime of feedback-modulated star formation. First, the final mass of newly formed stars in the kiloparsec-scale cold shells originated by the AGN feedback during the outbursts is significantly larger than in pure cooling flow models. From this point of view, the AGN activity induces star formation. However, the AGN activity also terminates star formation because, at the end of each major accretion episode, feedback is strong enough to empty the galaxy of gas on the kiloparsec scale. Finally, during the quiescent, low-luminosity phases of hot SMBH accretion, the smooth feedback helps keep the ISM hot, thus delaying a new cooling catastrophe and the associated star formation. It is tempting to identify these low-luminosity phases as the predicted counterparts of recent observational results that star formation is suppressed when the AGN is in the low-luminosity state (e.g., Schawinski et al. 2009; Kaviraj et al. 2011).

The previous papers, with the exception of Pellegrini et al. (2009), were mainly dedicated to the study of the accretion physics and the feedback effects. Here instead we focus on the appearance that the models would have if observed in the X-ray band, both in quiescence and during an outburst of activity. We concentrate on two models (named B2₀₂ and B3₀₂) extracted from the suite of cases presented in Paper III (Table 1 therein), characterized by a mechanical feedback efficiency dependent on the Eddington-scaled accretion luminosity. B2₀₂ and B3₀₂ were considered particularly successful, since their input parameters agree with previous theoretical studies or observations (as, e.g., for the AGN wind opening angle, and the peak value of mechanical efficiency, that are in accord with those estimated from 2D and 3D numerical simulations; Proga et al. 2000; Proga & Kallman 2004; Benson & Babul 2009), and, at the same time, their final properties (as the mass fraction of a younger stellar population, SMBH mass, etc.) are in reasonable accord with observations. The only difference in the input physics of B2₀₂ and B3₀₂ is the maximum value of the mechanical efficiency. Since in the simulations the treatment of feedback is physically based, not tuned to reproduce observations, any agreement or discrepancy of the resulting model properties with X-ray observations is relevant to improving our understanding of the SMBH-ISM coevolution, putting further constraints on the input ingredients, and possibly telling us what additional physics may be important in the problem. However, we stress that the models describe an isolated galaxy, where ram-pressure stripping (in the case of cluster galaxies) and intracluster medium pressure confinement (in the case of group or cluster central galaxies) are not taken into account (see Shin et al. 2010b).

The main observational signatures investigated here include the nuclear and gaseous emissions, and the ISM temperature and brightness profiles in the quiescent phases and before, during, and after a burst. Particular attention is paid to the appearance and detectability of central hot bubbles, with diameters of about a hundred parsecs, that are produced by the models during outbursts, and to various kinds of disturbances in the hot ISM. The analysis is performed in the soft and hard X-ray bands, also after unsharp-masking. These predictions are relevant for their observational consequences since the high angular resolution of the Chandra satellite has allowed us to obtain the best definition ever for the hot gas properties of the galaxies of the local universe by separating the contributions of stellar sources and hot gas, and the emission coming from different spatial regions within galaxies (e.g., Fabbiano 2012). In particular, in several elliptical galaxies various kinds of hot gas disturbances have been detected, likely resulting from nuclear activity (e.g., Finoguenov & Jones 2001; Jones et al. 2002; Forman et al. 2005; Machacek et al. 2006; O'Sullivan et al. 2007; Million et al. 2010). At the same time, nuclear emission values have been detected down to very low luminosities, comparable to those of X-ray binaries (e.g., Loewenstein et al. 2001; Gallo et al. 2010; Pellegrini 2010). Note that cosmological hydrodynamical simulations have already been attempted to study the effects of feedback from an AGN on the evolution of the galaxies, though mainly residing in clusters (e.g., Springel et al. 2005; Johansson et al. 2009; Sijacki et al. 2009; Teyssier et al. 2011; Hambrick et al. 2011; Kim et al. 2011; Kaviraj et al. 2011); the present attempt is the first of a detailed comparison of the X-ray properties resulting from simulations covering the past ∼10 Gyr of evolution, and sampling the galactic scale, with those of real galaxies.

The paper is organized as follows. Section 2 summarizes the main evolutionary phases of the representative models (B3₀₂ and B3₀₂), Section 3 describes how the observational properties of the models are derived, Section 4 presents a comparison of the nuclear luminosities with observations, Section 5 discusses the evolution of the X-ray luminosity and emission-weighted temperature of the ISM, and Section 6 presents the projected temperature and surface brightness profiles at representative times during quiescence and during a nuclear outburst. Finally, in Section 7 we summarize and discuss the main results.

The basic ideas behind the present class of models for feedback-modulated accretion flows have been introduced in Ciotti & Ostriker (1997, 2001), and Ostriker & Ciotti (2005), and developed in detail in the papers listed in the Introduction; a comprehensive recent discussion is given in Ciotti & Ostriker (2012). Here we focus on models B2₀₂ and B3₀₂ from Paper III (where the description of the numerical code and the input physics is given). The initial parameter values and the main final properties of the models are given in Table 1.

Table 1. The Representative Models

Model	Mw	〈w〉	〈EM〉	log ΔMBH	log ΔM*	log ΔMw	log Mgas	log LeffBH, opt/LEdd
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)
B202	10−3	2.0 × 10−5	0.105	8.74	9.74	10.27	9.68	−5.13
B302	3 × 10−4	1.2 × 10−5	0.133	9.05	10.22	10.31	9.34	−5.43

Notes. Relevant model properties at an age of 12 Gyr; masses are in units of solar mass M_☉, and luminosities in erg s⁻¹. The value of ^M_w is reached when L_BH ⩾ 2L_Edd, and the maximum radiative efficiency is set to 0.2. Columns 3 and 4 give the accretion weighted values of the mechanical and radiative efficiencies. ΔM_BH is the total accreted SMBH mass, ΔM_* is the total stellar mass formed during the evolution, ΔM_w is the total amount of ISM lost at 10R_e, and M_gas is the amount of gas inside 10R_e. L^eff_{BH, opt} is the fiducial SMBH luminosity in the optical as would be seen at infinity after absorption, with L^eff_{BH, opt} = 0.1L_BH at the first grid point (see Paper III for details).

Download table as:  ASCIITypeset image

The two models refer to an isolated elliptical galaxy placed on the Fundamental Plane, with a projected central stellar velocity dispersion σ = 260 km s⁻¹, a total B-band luminosity L_B = 5 × 10¹⁰L_B☉, and an effective radius R_e = 6.9 kpc. The stellar density profile is described by a Jaffe (1983) law, and the dark halo profile is such that the total (stellar+dark) mass density profile scales as ρ∝r⁻² at large radii; all relevant dynamical properties used in the code are discussed in Ciotti et al. (2009a). The dark-to-visible mass ratio is one within R_e, and the resulting stellar mass-to-light ratio is M_*/L_B = 5.8. Finally, the standard mass loss rate predicted by stellar evolution theory, and SNIa's rate declining with time t as t^−1.1, are assumed (e.g., Pellegrini 2012). The initial SMBH mass is set to 10⁻³ of the initial stellar mass M_*, i.e., it is 2.9 × 10⁸ M_☉. The simulations begin at a galaxy age of ∼2 Gyr (that is, a redshift z ∼ 2, the exact value depending on the epoch of elliptical galaxy formation, usually put at z ≳ 2).

The efficiency for producing radiation (Soltan 1982; Yu & Tremaine 2002) of material accreting on the SMBH at a rate of $\dot{M_{\rm BH}}$ is

$$\begin{equation} \epsilon =0.2\times {A \dot{m} \over 1+A\dot{m} }, \end{equation} \tag{ 1 }$$

where $\dot{m}=\dot{M_{\rm BH}}/\dot{M}_{\rm Edd}$ is the Eddington-scaled accretion rate. Thus, for A = 100, as adopted for the models, one has ∼ 0.2 at large mass accretion rates $\dot{m} \gg 0.01$, and declining as for radiatively inefficient accretion flows (RIAFs; Narayan & Yi 1994), as $\epsilon \sim 20 \dot{m}$, for $\dot{m}\lesssim 0.01$ (see also Section 4 for additional comments on this choice). Evidence for a transition at $\dot{m}\sim 0.01$ has been found recently by Trump et al. (2011). The mechanical feedback implemented is that of the Broad Line Region winds (leading to outflow velocities of ≃ 10⁴ km s⁻¹, similar to that observed by, e.g., Chartas et al. 2003 and Crenshaw et al. 2003), and reproduces the main features of the numerical modeling by Proga (2003). In particular, the mechanical efficiency scales with the Eddington ratio l = L_BH/L_Edd (where $L_{\rm BH}=\epsilon \dot{M_{\rm BH}}c^2$ is the instantaneous bolometric accretion luminosity, and $L_{\rm Edd}=\epsilon _0 \dot{M}_{\rm Edd} c^2$, with ₀ = 0.2; Paper III), reaching a maximum value ^M_w of 3 × 10⁻⁴ (for B3₀₂), and 10⁻³ (for B2₀₂), when l ≳ 2; also, the aperture solid angle of the conical nuclear wind increases with increasing l. Note that the values of the mechanical efficiency in Columns 2 and 3 of Table 1 are to be contrasted with the generally higher fixed efficiency of 5 × 10⁻³ commonly adopted in the literature (e.g., Hopkins et al. 2005; Di Matteo et al. 2005; Johansson et al. 2009). The mechanical output of a nuclear jet is also computed but not added to the hydrodynamical equations and it will be inserted in a future work.

The evolution of the gas flows is obtained by integrating the time-dependent (1D) Eulerian equations of hydrodynamics, with a logarithmically spaced and staggered radial grid, extending from 2.5 pc from the central SMBH to 250 kpc. It is most important that the resolution is high enough that the inner boundary is within the Bondi radius (Bondi 1952); if this is not ensured, the accretion rate will be calculated incorrectly. Typical values of the Bondi radius (for the central gas temperature and density of elliptical galaxies) range between 10 and 100 pc (e.g., Pellegrini et al. 2007b; Pellegrini 2010). Thus "Bondi accretion" is not assumed in the simulations; the correct, time-dependent accretion rate is computed from the hydrodynamical equations.

The code self-consistently derives the source and sink terms of mass, momentum, and energy associated with the evolving stellar population (stellar mass losses, SNIa events); the temporal steepening of the stellar velocity dispersion within the sphere of influence of the SMBH as a consequence of its growth; the star formation during nuclear starbursts; and finally accretion and SMBH feedback. Needless to say, the code conserves mass, energy, and momentum (Ostriker et al. 2010). Gas heating and cooling are calculated for a photoionized plasma in equilibrium with an average quasar spectral energy distribution (as detailed by Sazonov et al. 2005), and the resulting radiation pressure and absorption/emission are computed and distributed over the ISM from numerical integration of the radiative transport equation. The effects of radiation pressure on dust due to the starburst luminosity in the optical, UV, and infrared are also considered. A circumnuclear accretion disk is modeled at the level of "sub-grid" physics, and a set of differential equations describing the mass flow on the disk, its star formation rate, mass ejection, and finally SMBH accretion are solved at each time-step.

The resulting evolution of L_BH for B2₀₂ and B3₀₂ is shown in Figure 1. These models have fairly standard values of the input parameters, and their general properties in the X-ray band are typical of the 1D models investigated in Paper III. At the beginning, the galaxy is replenished by gas produced by the mass return from the evolving stars. Soon AGN outbursts develop, due to accretion of this gas, accompanied by star formation, and followed by degassing and a precipitous drop of the nuclear accretion rate. The outbursts are separated by long time intervals during which the galaxy is replenished again by gas from the stellar mass losses.

Zoom In Zoom Out Reset image size

The behavior of the gas during an outburst is almost independent of the specific burst episode considered. The outburst precursor is the off-center growth of a thin shell of dense gas (at a radius of ∼0.5–1 kpc) that progressively cools below X-ray-emitting temperatures and falls toward the center; compression of the enclosed gas follows and a central burst is triggered even before the cold shell reaches the center.⁴ The physical origin of the first cold shell (e.g., Ciotti & Ostriker 2007; Papers I and III) is a typical field instability. At variance with galaxy models with a flat stellar density profile in their central region, in cuspy models as those used here, the thermally unstable region, where the first cold shell originates at the beginning of each major burst, is located off-center, due to a combination of the galaxy density and stellar heating profiles. In less than a million years, a radiative shock originates from the center and quickly (in ≈10⁶ yr) produces an outward-moving shell that collides with the original shell falling in. Reflected shock waves carry fresh material for accretion, and produce new sub-bursts and new cold shells, hosting vigorous star formation (see the Introduction). This leads to the rich temporal structure of each outburst event, especially visible in model B3₀₂ (bottom panel, Figure 1).

Eventually, the cold material left after star formation (Ciotti & Ostriker 2007) is accreted in a final burst, a major shock leaves behind a very hot and dense center, and causes a substantial galaxy degassing. In general, while radiative effects mainly work on the kiloparsec scale, mechanical feedback from the AGN winds is more concentrated and affects the ISM on the ∼100 pc scale (see Figure 11 in Paper III). During the next few tens of million years, the gas cools, resumes its subsonic velocity, the density starts increasing again due to stellar mass losses, and the cycle repeats. At late epochs, the flow finally sets in a state of steady, hot, and low-luminosity accretion.

The observational model properties considered in Paper III were the nuclear bolometric, optical, and UV luminosities and the X-ray luminosity of the ISM within 10R_e. The latter was evaluated fiducially just by cutting the bolometric gas emission below a threshold temperature of 5 × 10⁶ K. In the following, we calculate the X-ray luminosity of the nucleus (L_{BH, X}, Section 4), and, with a more detailed and realistic procedure, the total luminosity and emission-weighted temperature of the hot gas (L_X and 〈T_X〉, Section 5); we also calculated the temperature and the surface brightness profiles during quiescent times and during an outburst (Section 6). We briefly describe below how the observational gas properties are calculated from the numerical outputs for the gas density and temperature.

The X-ray emission of the different model components is calculated over the energy range of 0.3–8 keV (the Chandra sensitivity band), and also in two separate bands, 0.3–2 keV ("soft") and 2–8 keV ("hard"), typically used in studies of the nuclear and gaseous properties. In practice, at any given time, the volume gas luminosity is calculated from the gas density and temperature distribution on the numerical grid as

$$\begin{equation} L_{\rm X}=4\pi \int _0^{\infty }{\cal {E}}(r) r^2 dr, \end{equation} \tag{ 2 }$$

where the emissivity is given by ${\cal {E}}(r) = n_{\rm e}(r) n_{\rm H}(r) \Lambda [T(r), Z]$, n_e and n_H are the number densities of electrons and hydrogen, and Λ(T, Z) is the cooling function. The cooling function is calculated over the two energy intervals by means of the radiative emission code APEC, valid for hot plasmas at the collisional ionization equilibrium (Smith et al. 2001), as available in the XSPEC package for the analysis of the X-ray data. For simplicity we choose constant abundance at the solar value, and the solar abundance ratios of Grevesse & Sauval (1998), which is consistent with observed gas metallicities (e.g., Loewenstein & Davis 2010). In order to speed up the analysis, we derived with APEC the values of Λ, for each energy band, for a large set of temperatures in the range 0.1–16 keV; then we obtained two very accurate non-linear fits of the tabulated values (with maximum deviations <1%, see Ciotti & Pellegrini 2008). These fits were used to compute the integral in Equation (2), and in every other integration where the emissivity is needed. For example the emission-weighted temperature for the whole galactic volume was computed as

$$\begin{equation} \langle T_{\rm X}\rangle = {4\pi \over L_{\rm X}}\int _0^{\infty }{\cal {E}}(r) T(r) r^2 dr. \end{equation} \tag{ 3 }$$

The surface brightness profile Σ(R), the emission-weighted projected temperature profile T_p(R), and the associated emission-weighted aperture temperature profile T_a(R), were obtained by numerical integration of the simulation outputs as

$$\begin{equation} \Sigma (R)=2\int _R^{\infty }{{\cal {E}}(r) r\over \sqrt{r^2-R^2}}dr, \end{equation} \tag{ 4 }$$

$$\begin{equation} T_{\rm p}(R) = {2\over \Sigma (R)} \int _R^{\infty }{T(r){\cal {E}}(r) r\over \sqrt{r^2-R^2}}dr, \end{equation} \tag{ 5 }$$

$$\begin{equation} T_{\rm a}(R) ={\int _0^R T_{\rm p}(R^{\prime })\Sigma (R^{\prime })R^{\prime } dR^{\prime } \over \int _0^R \Sigma (R^{\prime }) R^{\prime } dR^{\prime } }. \end{equation} \tag{ 6 }$$

The accuracy of the integrations above is verified by checking that the surface integral of Σ(R) over the whole grid recovers the same luminosity calculated via Equation (2), and that T_a(∞) = 〈T_X〉 within a few percent (Ciotti & Pellegrini 2008). The surface integral of Σ(R) is also used to compute the gas emission within the optical effective radius R_e, and in Equation (6) to compute the average temperature within the optical effective radius.

In order to highlight local departures from the mean ISM brightness profile, and to evidence major brightness features that could be revealed by observations, "fluctuation" profiles have also been created. These have been constructed with a technique similar to the so-called unsharp masking, frequently adopted in observational analysis (e.g., Fabian et al. 2003). In practice, the brightness profiles Σ(R) have been convolved with a 2D Gaussian of dispersion σ:

$$\begin{equation} {\rm PSF}={e^{-{R^2\over 2\sigma ^2}}\over 2\pi \sigma ^2}, \end{equation} \tag{ 7 }$$

so that the resulting surface brightness profile can be written as

$$\begin{equation} \Sigma _{\rm obs}(R)=\int _0^{\infty } I_0\left({RR^{\prime }\over \sigma ^2}\right) \,\, {e^{-{R^{\prime 2}+R^2\over 2\sigma ^2}} \over \sigma ^2} \,\, \Sigma (R^{\prime })R^{\prime }dR^{\prime }, \end{equation} \tag{ 8 }$$

where I₀ is the zeroth-order modified Bessel function of the first kind. In the analysis of the simulations, the integral above is solved numerically after a careful choice of σ. As expected, a too large σ produces an almost featureless profile, while a too small σ reproduces the unprocessed profile. After some attempts, it turned out that, in order to highlight local features, the optimal choice is that of a σ equal, at each gridpoint, to the sum of the lengths of the immediately preceding and subsequent grid intervals. The "unsharp-masked" profile is then defined in a natural way as

$$\begin{equation} \Sigma _{\rm UM}(R)\equiv {\Sigma (R)\over \Sigma _{\rm obs}(R)}-1. \end{equation} \tag{ 9 }$$

Figure 1 shows the time evolution of the nuclear bolometric accretion luminosity L_BH for B2₀₂ and B3₀₂ (whose input parameters differ only for the maximum value of the mechanical efficiency ^M_w, that is, respectively, 10⁻³ and 3 × 10⁻⁴; Table 1). Strong intermittencies at an earlier epoch, with L_BH reaching the Eddington value, become rarer and rarer with time, as the mass return rate from the stellar population declines, until a smooth, hot, and very sub-Eddington accretion phase establishes. The different mechanical efficiency is responsible for the sharp bursts in model B2₀₂, and the more time-extended and structured bursts in model B3₀₂ (Paper III). Toward the present epoch, at a galaxy age of 12 Gyr, the mass accretion rate on the SMBH for both models is $\dot{M_{\rm BH}}\approx 0.02\, M_{\odot }$ yr⁻¹, which translates into an Eddington-scaled accretion rate $\dot{m} \simeq 1.7\times 10^{-3}$ and $\dot{m} \simeq 1.2\times 10^{-3}$, respectively, for B2₀₂ and B3₀₂. The value of $\dot{m}$ of B3₀₂ is a bit smaller than for B2₀₂ because of its larger final M_BH (Table 2), a consequence of the weaker mechanical feedback. At the present time, and during the interburst periods, accretion has then entered the RIAF regime, and the radiative efficiency is ≃ 0.02; the nuclear bolometric luminosity is L_BH ≃ 2.4 × 10⁴³ erg s⁻¹ for both models, and the corresponding Eddington ratios are l ≃ 2 × 10⁻⁴ (B2₀₂) and l ≃ 10⁻⁴ (B3₀₂); see Table 2.

These results agree with the observation that in the local universe massive SMBHs are mostly radiatively quiescent, and the fraction of them at luminosities approaching their Eddington limit is negligible (e.g., Ho 2008). For example, in the sample of nuclei of the Palomar Spectroscopic Survey of northern galaxies, a nearly complete sample, magnitude-limited at B_T ⩽ 12.5 mag, ∼50% of ellipticals show detectable emission line nuclei,⁵ but mostly at low level (L_Hα < 10⁴⁰ erg s⁻¹; Ho 2008). For this sample, the nuclear bolometric luminosity (L_{bol, nuc}) was derived from the observed nuclear X-ray emission L_{X, nuc}, using the correction L_{bol, nuc}/L_{X, nuc} = 15.8 for the 2–10 keV band (Ho 2009). It was found that elliptical galaxies span a large range of L_{bol, nuc}, from 10³⁸ to 10⁴³ erg s⁻¹, with a median value of L_{bol, nuc} ≃ 1.7 × 10⁴⁰ erg s⁻¹ and a mean value of 4.6 × 10⁴¹ erg s⁻¹; the Eddington-scaled luminosity l has a median value of l = 1.2 × 10⁻⁶ (mean l = 1.2 × 10⁻⁵), with l < 10⁻³ for elliptical galaxies, and extending down to l = 10⁻⁸. The representative models tend therefore to lie on the upper end of the observed distributions of L_{bol, nuc} and l, a result that probably remains true regardless of the uncertainty in the bolometric correction from the 2–10 keV band. We will return to this point (already noticed in Papers I and III) in the conclusions.

Another way of comparing the model L_BH with observed values, is to estimate L_{BH, X} of the models, and compare it directly with observed L_{X, nuc} values. Accurate L_{X, nuc} measurements in a large number have been derived recently for elliptical galaxies of the local universe, based on Chandra data (e.g., Gallo et al. 2010; Pellegrini 2010); in the 0.3–10 keV band, L_{X, nuc} ranges from ≳ 10³⁸ erg s⁻¹ to ∼10⁴² erg s⁻¹, with most of L_{X, nuc}/L_Edd values as low as 10⁻⁶–10⁻⁸. In particular, for the final SMBH masses of the models (Table 2), it is observed that 10³⁸ ≲ L_{X, nuc} (erg s⁻¹) ≲ 10⁴², and 10⁻⁹ ≲ L_{X, nuc}/L_Edd ≲ 10⁻⁴, both from the sample of the Hubble Virgo Cluster Survey (Gallo et al. 2010), and that of 112 early-type galaxies within ∼60 Mpc (Pellegrini 2010). The 0.3–10 keV nuclear luminosity L_{BH, X} of the models at the present epoch can be derived adopting a correction factor appropriate for the spectral energy distribution of a radiatively inefficient accretion flow, i.e., L_{BH, X} ≲ 0.2L_BH for low-luminosity AGNs (Mahadevan 1997). This gives L_{BH, X} ≲ 5 × 10⁴² erg s⁻¹, and L_{BH, X}/L_Edd ≲ 4.8 × 10⁻⁵ (Table 2). Also in the X-ray band, then, the model nuclear luminosity tends to be larger than that typically observed; L_{BH, X}/L_Edd is within the observed range (in its upper part). All this may suggest that in real galaxies an additional mechanism is further reducing the mass available for accretion, as could be provided by a nuclear jet and/or a wind from an RIAF (Narayan & Yi 1994; Blandford & Begelman 1999). In the latter case, only a fraction of the gas within the accretion radius actually reaches the SMBH; the binding energy released by the accreting gas is transported radially outward to drive away the remainder in the form of a wind. Of course, at low accretion luminosities, the AGN wind mechanical feedback is not important, so that a change in its efficiency does not fix the nuclear overluminosity of the quiescent states. On the other hand, adopting a high fixed value of the mechanical feedback efficiency to simulate the role of jets (e.g., Section 2; Di Matteo et al. 2005) produces a full degassing of the models and then a much too low ISM X-ray luminosity (see also Ostriker et al. 2010).

Alternatively, it may be that the quadratic dependence of l on $\dot{m}$ sets in at a higher value of $\dot{m}$ than adopted here (Section 2), for which it starts at $\dot{m}\sim 10^{-2}$. Within the current uncertainties, we might have chosen the constant A = 10 rather than 100; in this way, the quadratic dependence would have set in at l ≲ 0.1 rather than l ≲ 0.01. This would have reduced the late time nuclear luminosity by a factor of ≃ 9, as confirmed by a supplementary run⁶ of models B2₀₂ and B3₀₂ with A = 10.

Finally, another interesting—albeit more difficult—comparison with observational results can be done using the duty cycle. The latter can be calculated as the fraction of the total time spent by the AGN in the "on" state, defined by a luminosity greater than L_Edd/30 in a given band, over some temporal baseline (Paper III).⁷ In so doing, the duty cycle turns out to be a small number (Table 2); for example, for a temporal baseline ending at the present epoch (i.e., at 13 Gyr in Figure 1), the duty cycle of model B2₀₂ is zero when starting from 9 Gyr (at a redshift z ≈ 0.45, for a flat universe with H₀ = 71 km s⁻¹ Mpc⁻¹, Ω_M = 0.27, Ω_Λ = 0.73), as no burst occurs after ≃ 7.5 Gyr; when starting from 6 Gyr (z ≈ 1), it is ≃ 6.3 × 10⁻³, 3.2 × 10⁻³, and 3.0 × 10⁻³, respectively, in the bolometric, optical (absorbed), and UV (absorbed) bands. For model B3₀₂, in the same bands, the duty cycle is, respectively, ≃ 4.9 × 10⁻², 1.7 × 10⁻², and 7.1 × 10⁻³ (when starting at 9 Gyr), and ≃ 4.8 × 10⁻², 1.8 × 10⁻², and 8.6 × 10⁻³ (when starting at 6 Gyr). The duty cycle decreases from the bolometric, to the optical, to the UV bands, due to the different values of the opacity in these bands.

These duty-cycle values are broadly consistent with the fraction of active galaxies measured in observational works. For example, Greene & Ho (2007) estimated the (mass dependent) number of active galaxies, using broad-line luminosities from Sloan Digital Sky Survey (SDSS) DR4, for galaxies with z < 0.352 (age of the universe ≳ 10 Gyr); statistically speaking, the fraction of active systems can be interpreted as a duty cycle for SMBHs in a given mass range. Greene & Ho reported duty-cycle values of the order of 4 × 10⁻³ for SMBHs of masses of 10⁷ M_☉, declining at increasing mass. Similar duty-cycle values of ∼2 × 10⁻³, decreasing at increasing SMBH mass, were reported by Heckman et al. (2004). The duty cycles of the models tend to be larger than those observed; however, the comparison is limited by the small number of models considered, and the fact that the only way to compute duty cycles different from zero is to extend the temporal baseline back in the past. A more consistent comparison with observations requires an increased data set of models, and computing the duty cycle for the last 2–3 Gyr. Clearly this procedure will reduce the duty cycle as the models are characterized by a declining nuclear activity. Again, the computed duty cycle also would have been reduced significantly if we had raised the threshold for RIAF-like behavior of the radiative efficiency to l = 0.1.

We describe here the time evolution of the global thermal luminosity and temperature of the ISM.

5.1. Luminosity Evolution

The top panels of Figure 2 show the time evolution of the total gas emission L_X for models B2₀₂ (left) and B3₀₂ (right). Red lines refer to the soft (0.3–2 keV) band, and blue lines to the hard (2–8 keV) band; for reference, the scaled-down L_BH (black dotted line) is also shown.

Zoom In Zoom Out Reset image size

The luminosity evolution of the gas for the two models is qualitatively similar, with peaks in L_X coinciding with those in L_BH. The sudden increase of the gas emission during outbursts is due to the increase in temperature and density in the central galactic regions (≃ 10²–10³ pc) caused by radiative gas heating (Compton and photoionization) and by compression due to direct and reflected shock waves produced by mechanical and radiative feedback which are associated with the AGN and the starburst. For short times, most of the luminosity in the peaks of L_X originates from a very small region at the galactic center (≈100 pc), thus it is observationally indistinguishable from the luminosity of the nucleus (see also Section 6.2). The hard emission oscillates in phase with the soft one, and presents the same overall behavior, but remains at a level almost two orders of magnitude lower. Hard emission during the outburst, as shown in Figure 2, would be detectable with Chandra, if centrally concentrated (see also Section 6.2 below). However, hard emission during quiescent times would be difficult to distinguish from the contribution of unresolved binaries, even with Chandra, if extended (e.g., Pellegrini et al. 2007a; Trinchieri et al. 2008).

A comparison of the peaks in L_X and L_BH reveals differences and analogies. While L_BH shows sharp and sudden spikes at the outbursts (increasing by two or more orders of magnitude in ≈10⁶–10⁷ yr), and is almost constant between them, L_X increases slowly but steadily between outbursts when the galaxy is replenished by the stellar mass losses. The peaks in L_X become broader with increasing time, but not weaker; for example, the increase of L_X during the last major outburst of B2₀₂ is the largest one, with the largest amount of gas heated and then removed from the galaxy in its entire life. Instead, when the burst ends, L_X has the same abrupt decrease as L_BH, due to the density drop following the final (and usually strongest) sub-burst in each major accretion event. Another similarity is that both L_BH and L_X show sharper and "cleaner" bursts in B2₀₂ than in B3₀₂: more radiatively dominated feedback bursts (as in B3₀₂) are richer in temporal substructure, because it takes longer for the cold shells to be destroyed, thus more star formation and SMBH accretion occur (Paper III).

The quiescent values of L_X during the past few Gyr remain at the same level of ∼10⁴⁰ erg s⁻¹ for model B2₀₂, and decrease by a small factor of ≲ 2 to reach a present epoch L_X = 2 × 10³⁹ erg s⁻¹ for B3₀₂. Previous compilations of observed L_X values (Fabbiano et al. 1992; O'Sullivan et al. 2001; Diehl & Statler 2007; Mulchaey & Jeltema 2010) for early-type galaxies of the local universe, residing in all kinds of environments (from the field to groups to clusters such as Virgo and Fornax), show a range of L_X from 10⁴⁰ erg s⁻¹ up to 10⁴³ erg s⁻¹, at a B-band optical luminosity of L_B = 5 × 10¹⁰L_{B, ☉} as the model galaxy. L_X values larger than a few× 10⁴¹ erg s⁻¹ belong to galaxies at the centers of groups, clusters, or subclusters, for which an important contribution from the intragroup or intracluster medium, or a confining action, is likely (e.g., Renzini et al. 1993; Brighenti & Mathews 1998; Brown & Bregman 2000; Helsdon et al. 2001). However, the final L_X of B2₀₂ and B3₀₂ is on the lower end of the range of observed values; this indicates that degassing is important in the models, and for many real cases it must be impeded. In the simulations this could be obtained, for example, by adding the external pressure from an outer medium (e.g., Vedder et al. 1988), and it will be considered in future works.

5.1.1. L_X − L_B and L_X − L_BH

Real galaxies show a wide range of L_X values, and the observed L_X variation has remained a subject of debate for years (e.g., Fabbiano 1989; Ciotti et al. 1991; White & Sarazin 1991; Pellegrini 2012). Thus, we check here whether the L_X variation in the models during their evolution can (partly) account for the large observed variation. For this check we considered the range of hot gas emission for the largest sample of early-type galaxies of the local universe (O'Sullivan et al. 2001), after having excluded AGN-dominated cases, and central dominant cluster or group galaxies. Only galaxies with optical luminosities within a range close to that of the model galaxy have been considered (i.e., from log L_B = 10.5 to log L_B = 10.8, for which there are 43 galaxies). The discrete stellar sources' contribution estimated by O'Sullivan et al. (2001) has been removed. The distribution of the observed X-ray emission values so obtained (N, normalized to the total number of galaxies) is then compared with that built for the soft X-ray emission of the models during the epoch from 2 to 12 Gyr (Figure 3). Each model emission value enters the histogram with the fraction (again N) of the chosen epoch during which it is present (here then N is normalized to the total chosen epoch); the hypothesis underlying this comparison is that, statistically, an observed galaxy can be caught in any one of the phases shown in the past 2–12 Gyr by the representative models.

Zoom In Zoom Out Reset image size

Figure 3 shows that L_X of the models keeps within the observed range, but it does not significantly exceed ∼10⁴¹ erg s⁻¹, while a fraction of galaxies populates this region. Model B2₀₂ has a ∼constant L_X in the past few Gyr, so that it populates mostly a couple of bins; model B3₀₂ (which experiences more outbursts) reaches a wider coverage of the observed L_X range, but its L_X distribution extends more to lower L_X values than to larger ones, due to a larger overall degassing. Larger L_X for the models can be obtained when considering that real galaxies of similar L_B can have different values of the central stellar velocity dispersion and effective radius. These differences determine a variety of flow evolution, and then of L_X (Ciotti & Ostriker 2012). Following this idea, Figure 3 (bottom left) shows the histogram of one possible variation to model B3₀₂, that with σ = 280 km s⁻¹; this indeed reaches larger L_X values. Finally, the bottom right panel shows the average of the histograms of the three models, and shows how this average reproduces the observed histogram reasonably well. Overall, this analysis shows that the large observed variation in L_X at similar optical L_B has another contributing factor, nuclear activity, in addition to those already put forward.

Another useful diagram for an observational comparison is given by the relationship between L_BH and L_X; this is shown in Figure 4, considering all the available temporal outputs, for model B2₀₂. The analogous figure for model B3₀₂ is very similar, with the only difference being the broader temporal extension of the bursts. One can recognize the interburst times, when the galaxy is replenished by stellar mass losses, as periods with L_BH almost constant, and L_X increasing (which produces almost vertical lines). During outbursts, L_X and L_BH abruptly increase and then decrease, which produces loops running clockwise on the right of each vertical line; these loops are occupied for a very short time. As time proceeds, the vertical lines (the interburst periods) shift toward lower L_BH. The final quiescent period is described by a line (the most crowded with points) where L_BH and L_X both decrease, though L_BH decreases faster than L_X. During the bursts L_BH usually peaks at a later epoch than L_X, since the maximum L_X is reached at the cold shell formation, when L_BH is still increasing due to the increasing density of the hot gas inside the shell. When L_BH reaches a maximum, the gas density and temperature have usually been modified and produce a lower L_X.

Zoom In Zoom Out Reset image size

The more the outbursts are confined at earlier epochs, the more the final weak correlation between L_X and L_BH extends with time; the more the outbursts extend toward the present epoch, the more a trend is expected for L_X to increase with L_BH with a large scatter around it (as produced by the vertical lines, where galaxies reside most of the time, and by the loops). For a sample of early-type galaxies of the local universe, the relationship between L_X and the nuclear emission L_{X, nuc} indeed shows a weak trend, dominated by a large scatter (Pellegrini 2010); in the present framework, such an observation could be explained with many galaxies still being in the phases made of the vertical rise followed by the loop. The comparison cannot be pushed farther, though, since the L_BH values—when converted to an X-ray band—are typically larger than the observed L_{X, nuc} (Section 4).

5.2. Temperature Evolution and L_X − T_a

Another important global property of the ISM typically observed is the emission-weighted aperture temperature T_a (Equation (6)), here calculated within the optical R_e, and within 10R_e. The time evolution of these temperature diagnostics is shown in the bottom panels of Figure 2, in parallel with the luminosity evolution; red and blue lines again refer to the 0.3–2 keV and 2–8 keV bands, respectively. Note that T_a(10R_e) in the two bands is indistinguishable from the corresponding global emission-weighted temperature 〈T_X〉 (which is calculated but not shown in Figure 2), as expected given that the density profile is steeply decreasing outward (Paper III); thus T_a(10R_e) is a good proxy for 〈T_X〉. Also, temperatures computed for the whole 0.3–8 keV band (not shown) are always very close to those weighted with the 0.3–2 keV emission, except for those short burst times during which there is a very hot gas component. As for L_X, also for T_a the temperature peaks of B3₀₂ are significantly more structured in time than those of B2₀₂.

Three main characteristic features of T_a common to this class of models can be pointed out. The first is the complex behavior of T_a during outbursts, with variations going in opposite directions for the two bands. This is due to the coexistence of hot (the central bubble) and cold (the radiative shells) ISM phases during the bursts. The sharp and high peaks in the hard band (blue) correspond to the onset of very hot regions at the center, while the decrements in the soft band (red) are due to a dense cold shell created immediately before the major burst, and to cold gas accumulated by the passage of radiative shock waves produced by the outburst (see also Section 6.1).

A second, observationally relevant feature is that in each band T_a(R_e) is higher than T_a(10R_e), which is especially evident in the interburst periods (Figure 2). For example, in the quiescent phases, both models show in the soft band a similar T_a(10R_e) ≃ 0.4–0.5 keV, while T_a(R_e) ∼ 0.7 keV. This is explained by the radial temperature distribution decreasing outward (Section 6.1).

Finally, for model B2₀₂ Figure 5 shows the relationship between L_X and T_a(10R_e), both calculated for the 0.3–8 keV band; such a diagram is often produced in observational works, for galaxies of all σ (e.g., Pellegrini 2011). During the long interburst epochs, L_X increases with little variation of T_a. In the short burst episodes, L_X and T_a reach a maximum, though at slightly different times: at each burst, the density in the central regions decreases when the shock expands, with a decline in L_X; after the shock has crossed the inner regions, and increased the post-shock gas temperature, there is the slightly delayed maximum in T_a. Then L_X and T_a follow a clockwise loop on the right, reach back to the original T_a value, move left of it, and follow a counterclockwise smaller loop. Then the cycle repeats with L_X increasing and T_a almost constant. During the last few Gyr, L_X remains at ∼10⁴⁰ erg s⁻¹ and T_a ∼ 0.5 keV.

Zoom In Zoom Out Reset image size

We now compare these results with observations. In the largest sample of global, soft-X-ray emission-weighted temperatures derived from ROSAT observations 〈T_X〉 is in the range ∼0.4–0.8 keV for galaxies with σ ≃ 260 km s⁻¹ as for the model galaxy (O'Sullivan et al. 2003); T_a(10R_e) of the models during quiescence, in the soft band, is within this range, though on its lower end. However, various factors tend to bias high these observed temperatures as incomplete subtraction of the hard stellar emission due to binaries or hard AGN emission; in addition, many of the sample galaxies reside in high density environments with possible contamination from the hotter group/cluster medium and a temperature profile that is commonly rising outward (e.g., Diehl & Statler 2008; Nagino & Matsushita 2009), a behavior opposite to that of the models which refer to isolated galaxies (Section 6.1).

Chandra observations, with large sensitivity and much higher angular resolution (∼1''), allowed for an accurate subtraction of all the AGN and stellar source contributions from the total emission, giving measurements of the gas temperature of unprecedented accuracy (e.g., Boroson et al. 2011). For example, Boroson et al. derived global, 0.3–8 keV emission-weighted temperatures for a few galaxies with σ ∼ 260 km s⁻¹ ranging from 0.3 to 0.6 keV; T_a(10R_e) of the models agrees well with this result, falling in the middle of the observed range. Coming to temperature estimates for more central regions, Athey (2007) derived 0.35–8 keV emission-weighted temperatures within R_e, for 53 galaxies with Chandra observations. For a selection of 20 galaxies with L_B similar to that of the model galaxy, from log L_B(L_{B, ☉}) = 10.5 to 10.8, the average temperature is 0.61 ± 0.03 keV (calculated weighting each measurement with its uncertainty). Figure 6 shows a full comparison between the distribution N of these temperatures and of T_a(R_e), calculated for the 0.3–8 keV band, during the evolution of B2₀₂ and B3₀₂ (N for the observed values and the models is calculated as in Section 5.1.1). The model T_a(R_e) tends to be concentrated in the upper half of the observed distribution; this result, if confirmed with a larger set of simulated galaxies, indicates that heating in the central galactic region of the models may be too efficient.

Zoom In Zoom Out Reset image size

In conclusion, the global temperatures of the models fall in the middle of the range of values recently observed, while the model T_a(R_e) reproduces easily the larger observed values and less easily the lower ones.

We consider here the radial profiles of the ISM temperature and surface brightness, as they would be revealed for the models by observations. For simplicity we limit the discussion to model B2₀₂ whose sharp bursts allow for an easier presentation; the results are substantially similar for B3₀₂. The profiles are constructed using Equations (4)–(9) both during the quiescent phases, which occupy most of the ISM evolution, and during an outburst. The recurrent burst phases are temporally limited, but represent a central aspect of the models, thus they are devoted special attention. Statistically, the feedback features should be present, and possibly revealed by current X-ray observations, in ≃ 5%–10% of the isolated galaxies with L_B similar to that of the models (Section 6.2 below). In the following we present snapshots of the projected profiles during quiescence, at an age of 3, 6.5, 9, and 12 Gyr, and centered on the last outburst at 7.5 Gyr.

6.1. Temperature Profiles

The emission-weighted projected temperature profiles T_p(R) during quiescent interburst times are smooth, with the temperature monotonically decreasing for increasing radius (Figure 7). From an age of 6 Gyr onward, the temperature keeps at ≃ 1 keV at a radius of ≃ 100 pc, and at ∼0.4 keV at a radius of ≃ 30 kpc. Figure 7 (right panel) shows the corresponding aperture temperature profiles T_a(R) calculated in bins (i.e., in Equation (6) the numerator and denominator are evaluated over radial bins), chosen to reproduce the typical binning used for Chandra observations of nearby galaxies (Humphrey et al. 2006; Diehl & Statler 2007).

Zoom In Zoom Out Reset image size

Temperature profiles with negative gradients had already been found for gas inflows in steep galactic potentials without feedback, due to compressional heating (e.g., Pellegrini & Ciotti 1998). As typical of models without a central SMBH, in that case the temperature also had a central drop. In the present computations, instead, the combined effects of the gravitational field of the SMBH, and of the high injection temperature of the stellar mass losses (a consequence of the stellar velocity dispersion that is enhanced by the presence of the SMBH, within its sphere of influence), keep the gas temperature increasing approaching the center, even outside the burst episodes (see also Pellegrini 2012). A temperature profile that keeps smoothly decreasing toward the center for a long time, as expected in classical cooling flows, is never shown by the model runs.

With increasing time, the value of the central temperature does not increase significantly after 6 Gyr, since it is mainly determined by the gravitational effects of the SMBH, whose mass remains approximately constant. The external T_p(R) values instead steadily increase with time since they are determined by the SNIa's heating; the latter has a secular trend that produces a time-increasing specific heating of the mass return rate (i.e., $L_{\rm SNIa}/\dot{M}_*\propto t^{0.2}$, where L_SNIa is the energy injected per unit time by SNIa's supernova explosions, and $\dot{M}_*$ is the stellar mass loss injected per unit time; e.g., see Pellegrini 2012).

We now focus on the evolution during the last outburst (Figure 8). As typical, starting from the unperturbed profile, a shell of denser gas is created, in this case at a radius of ≃ 800 pc; it is particularly evident as a dip in the otherwise monotonically decreasing profile (see the −2 Myr red line in the top left panel). The shell falls toward the center, progressively cooling, so that the soft X-ray emission-weighted temperature decreases (Figure 2, bottom panels). The first AGN burst is produced before the shell reaches the center, and the resulting outward-moving shock heats the gas in the inner regions of the galaxy on a timescale of ≈1 Myr. The subsequent snapshot in Figure 8, at +6 Myr after the first burst, shows the high central temperature typical of the outburst phase: T_p(R) ∼few keV within a radius of ∼50–100 pc.

Zoom In Zoom Out Reset image size

When the shock enters the radiative snow-plow phase, the associated secondary dense shell interacts with the first one still falling and a series of direct and reflected shock waves are produced accompanied by accretion events and star formation. The profiles between −2 Myr and +18 Myr (not shown) are very irregular and consist of a series of dips propagating outward with the temperature in the central region quickly and alternately increasing and decreasing. As already noted for Figure 2, both hotter and colder regions are continuously created, the hottest ones within ∼1 kpc (mainly due to shocks), and the coldest ones due to the cooling gas in the shells.

Each sub-burst is stronger than the previous one; finally, a last shock from the center concludes the phenomenon, halting star formation and leaving behind a very hot nucleus (the +66 Myr red line), while emptying the rest of the galaxy (see also the Σ(R) profile at +66 Myr in Figure 10). Then, during the following ≃ 100 Myr, the gas resumes the temperature typical of quiescent times (the +202 Myr green line). The inner very hot phase lasts for ≲ 0.1 Gyr.

During the bursts, cosmic rays are shock-accelerated, and the inner regions look similar to a gigantic supernova remnant. Only thermal X-rays are considered here while those due to synchrotron emission from accelerated particles due to shocks were not computed (see Jiang et al. 2010).

The profiles in the right panels of Figure 8 show the binned aperture temperature profiles for the same times as for the left panels. The binning smears the small-scale features, but the major distinctive characteristics, as the temperature dip when the first shell forms, and the large temperature drop outside the center, at +18 Myr and +66 Myr, are still present.

6.1.1. Comparison with Observed Profiles

Negative radial gradients, as shown by the model temperature during quiescence (Figure 7), are common among ellipticals, as revealed most recently by Chandra observations (e.g., Kim & Fabbiano 2003; Humphrey et al. 2006; Sansom et al. 2006; Fukazawa et al. 2006; Diehl & Statler 2008; Nagino & Matsushita 2009). In a large collection of temperature profiles (Diehl & Statler 2008; Athey 2007), those cases with negative gradients show a temperature that roughly halves from the innermost bin (that in general extends out to a radius of 0.5-a few kpc from the center) to the outer galactic region; this is roughly the same behavior shown by the models. Observed temperatures range between 0.5 and 1 keV at the innermost radial bin, where the model temperature is 0.8–1 keV (Figure 7, right panels). Interestingly, the galaxies with a negative gradient reside in the field or in less dense environments, as the outer Virgo regions (Matsushita 2001; Fukazawa et al. 2006; Diehl & Statler 2008), and all have⁸L_X < few × 10⁴⁰ erg s⁻¹, both characteristics that match those of the models.

More complex temperature profiles are also common, and could correspond to phases of the AGN activity. For example, the so-called hybrid profiles show a central negative gradient, until the temperature reaches a minimum, and an outer positive gradient; for example in NGC 1316, NGC 4552, and NGC 7618 the temperature has a minimum of ∼0.4–0.5 keV at ∼few kpc, and rises both going toward the center (of ∼0.2–0.3 keV) and going outward (Diehl & Statler 2008). Also in the models, after each major burst, when the last shock is moving outward and fading, leaving a hot, rapidly cooling core, there is a drop in the temperature profile at a radius of ∼1 kpc or more (Figure 8, bottom panels, black and red profiles). However, in the models the temperature reaches >1 keV at the center, larger than the observed values (the comparison also depends on the binning used for the profiles, though). Interestingly, preliminary results from 2D simulations (Novak et al. 2011), show lower central temperatures than in Figure 8 during outbursts due to the fragmentation of the cold shell which causes a lower gas compression while falling to the center.

Other observed profiles stay roughly isothermal or are roughly flat out to ∼R_e and then increase outward (e.g., Humphrey et al. 2006; Diehl & Statler 2008; Nagino & Matsushita 2009). Such a positive outer gradient is shown by galaxies in high-density environments, suggesting the influence of circumgalactic hot gas; also, these galaxies have typically a larger hot gas content than the models. Environmental confining effects, currently not included in the models, are expected to increase the gas luminosity and produce outer positive gradients (see, e.g., Sarazin & White 1987; Vedder et al. 1988).

It has been proposed that galaxies could behave according to a scenario where weak radio AGNs distribute their heat locally and host negative inner temperature gradients, whereas more luminous radio AGNs heat the gas more globally through a jet or rising bubbles and produce a flat profile or a positive gradient (Diehl & Statler 2008). Our models during quiescence could correspond to the weak AGN phase; also, after an outburst, the temperature profile can be rather flat (excluding the innermost bin) out to a few kiloparsecs, and show a positive gradient outward of this radius (see the red line of Figure 8, bottom right panel). However, without a confining environment, an exploratory investigation conducted by us shows that the kinetic heating of a bright radio AGN could cause a major degassing rather than just a reversal of the temperature gradient from negative to flat or positive. A gas-rich environment providing confinement for the galactic coronae is likely a necessary condition to observe a bright (extended) radio source, and the confinement in turn also helps produce a positive temperature gradient in the outer galactic regions. Both aspects (the jet and a dense environment) will be implemented in models in the future.

In conclusion, the addition of a jet to the simulations could have positive effects if it heats the gas outside ∼R_e, and reduces the accretion rate and then L_BH during the stationary hot accretion phase (Section 4); it should not increase the temperature at the center, though, since this is already on the upper end of those observed even during quiescence.

6.2. Brightness Profiles

Figure 9 shows the evolution of the X-ray surface brightness profile of the gas Σ(R) for B2₀₂ at the same quiescent times of the temperature profiles in Figure 7 for the soft (0.3–2 keV) and hard (2–8 keV) bands. In the hard band, Σ(R) is always comparable to or lower than a profile representing the unresolved stellar emission due to low-mass X-ray binaries, calculated for a deep (≳ 200 ks) Chandra pointing of a galaxy of the same optical luminosity as B2₀₂ and distance ≲ 20 Mpc. Thus, hard emission during quiescent times would be difficult to distinguish from the contribution of unresolved binaries, even with Chandra. In both bands Σ(R) becomes flatter in shape with increasing time, since the emission level decreases mostly in the central galactic regions (within R_e). This important effect is produced by the nuclear outbursts which remove gas from the center.

Zoom In Zoom Out Reset image size

Analogous to the discussion of the temperature profiles, Figure 10 shows Σ(R) just before, during, and after the last major nuclear burst at 7.498 Gyr, with times counted from the first accretion event. At −2 Myr the formation of the off-center cold shell produces the characteristic feature of a sharp decrease of Σ in the hard band and a bright rim in the soft band. The subsequent curves show the shock moving outward after the major burst and the presence of very hot gas at the center revealed by the central peak in the hard band. The disturbances in the Σ profiles due to the shells launched by the repeated sub-bursts are much more visible in the soft than in the hard band, as particularly apparent at +18 Myr. At this time, note the presence of a hot center surrounded by a denser and colder shell, producing a sharp peak in Σ(R) at R = 600–700 pc and a sharp dip in T_a(R) in Figure 8. The final two times (+66 Myr and +202 Myr) show the result of the degassing caused by the passage of the shock waves: the gas density is low and subsonic perturbations remain at a radius of ≳ 10 kpc. A different representation of the profiles during the burst phases is given by Figure 11, where Σ_UM(R) resulting from the unsharp masking procedure (Equation (9)) is shown. Fluctuations in brightness that can be clearly distinguished are those due to the cold shell (time −2 Myr), and, after the burst, the negative off-center region in Σ_UM(R) delimited by a sharp positive peak (times +6, +18, +66 Myr); the latter resembles what is commonly called a hot gas "cavity."

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Note that, when implemented in 2D simulations, the same input physics adopted for the present class of models give conical hot gas outflows from the nucleus during outbursts: hot gas is ejected along the symmetry axis so that elongated "cavities" (i.e., regions with a gas density decrement) are created (Novak et al. 2011). In addition, in their study of a model very similar to B3₀₂, Jiang et al. (2010) found after a burst taking place at 6.5 Gyr, a cavity filled with radio-emitting particles of ∼4.4 kpc in size, detectable during the first ∼10 Myr after the burst. Therefore, we expect that cavities in the hot gas such as those seen as off-center minima in the profiles of Figures 10 and 11 at times from 6 to 66 Myr after the burst should be fairly bright in the radio and would have two essentially hollow conical lobes in X-rays. They should even show non-thermal X-ray emission of the type seen in the Crab nebula (Jiang et al. 2010). Finally, taking the results of the present analysis and the estimates on the radio detectability of Jiang et al. (2010) at face value, hot gas cavities seem more long-lasting than their radio detectability.

Coming to the observability of the predicted features, one major property of the models is the decrease of Σ(R) in the central galactic regions produced by the nuclear outbursts with respect to models without feedback. Interestingly, bright ellipticals imaged with Chandra (e.g., Loewenstein et al. 2001) show a brightness profile that is quite flat within the central ∼1 kpc, a feature impossible to reproduce with pure inflow models, while it resembles the profile of the "pre-burst" phase (black lines in Figure 10, left panels) or at the end of a burst (+202 Myr). To better illustrate this point, Figure 12 shows the different shape of Σ(R) for B2₀₂ and for a model with the same L_B, σ and R_e, but without the AGN feedback, and the sole SNIa's heating included. The two Σ(R) profiles considered refer to the present quiescent epoch, yet the difference in steepness is clear. Figure 12 also shows the observed Σ(R) for an elliptical at the periphery of the Virgo cluster, with L_X close to that of B2₀₂; the agreement between model and observation is very good.

Zoom In Zoom Out Reset image size

Another major prediction is given by the disturbances in the profile during an outburst; these stay above the level of the unresolved stellar emission for a deep observation of a nearby galaxy with Chandra, as shown by Figure 10. The central spike in Σ(R) during the high-temperature and high-density phase at the center, is confined within ∼100 pc, and then likely to remain a central unresolved feature even in galaxies observed at the high angular resolution of Chandra. Disturbances such as shells and ripples farther out in the galaxy last ≲ 0.2 Gyr and are more likely to be observed. Given the typical duration of these features and the presence of three to four major outbursts during the whole evolution, statistically they should be present and possibly revealed by current X-ray observations in ≃ 5%–10% of the galaxies with L_B similar to that of the model galaxy.

In fact, many nearby galaxies show a disturbed appearance, as revealed by studies based on Chandra data (e.g., Finoguenov & Jones 2001; Forman et al. 2005; Soria et al. 2006; Diehl & Statler 2007; Nulsen et al. 2009; Baldi et al. 2009; Dunn et al. 2010). Many of the best-studied gas-rich galaxies show decrements in the X-ray surface brightness map, identified as cavities formed when AGN jets inflate radio lobes and displace surrounding gas; in many cases the cavities are filled with radio plasma and surrounded by armlike features sometimes classified as shocks. The hot gas disturbances have then been generally attributed to jet activity. As discussed above, we expect that if we took two cones from our solutions, they would give "lobes." These would be accompanied by shocks at the edges and would be filled with radio-emitting particles.

There are also a few galaxies without currently evident extended radio emission, but with signs of an outflow and hot central gas, such as NGC 4552 (Machacek et al. 2006). This galaxy shows a weak core radio source unresolved by the Very Long Baseline Array, and in the (unsharp masked) X-ray image two conspicuous ringlike features at 1.3 kpc from the galaxy center, surrounding two cavities; these features have been found to be consistent with shocked gas driven outward by recent nuclear activity, as in a bipolar nuclear outflow. The gas temperature in the central ∼100 pc of the galaxy is 1 ± 0.2 keV, hotter than elsewhere in the galaxy, suggesting that we may be directly observing the reheating of the galaxy ISM by the outburst (Machacek et al. 2006). These characteristics resemble those predicted by the models for the temperature and the surface brightness during the afterburst phase (Figures 8, 10, and 11).

The hot gas properties of massive ellipticals with regard to luminosity and temperature and their spatial distributions, allow us to derive insight into the hot gas evolutionary status and its link with the host galaxy. Since the gas cooling times are short compared to the galaxy age, it is now commonly accepted that repeated cooling catastrophes have occurred in the past, accompanied by central starbursts and AGN outbursts. The interest in a better understanding of this phenomenon is obvious, as it is directly linked to galaxy formation and evolution and to the growth of the central SMBH. Unfortunately, a complete theoretical picture is still missing so that modeling coupled with a close comparison with observations is crucial in this field. These feedback events must leave signatures on the X-ray properties of the galaxies; indeed, the observed temperature and brightness profiles often cannot be fit with smooth profiles, such as those predicted for galactic winds or cooling flows (e.g., Sarazin 2012; Statler 2012). In this work we have calculated the observational properties of two galaxy models in the X-ray band, representative of a class of high-resolution 1D simulations of physically based models for feedback-modulated accretion in isolated galaxies (taken from Paper III). The feedback physics include the combined effects of radiation and AGN winds. The observational properties derived provide good matches to those observed in general for the local universe, and account for a few otherwise puzzling observed properties; on the other hand, they also evidence the need for improvements or additions to the input physics. The main results of the present investigation are the following.

• 1.  
  After an evolution of ≃ 10 Gyr, the models are typically in a permanent quiescent phase. The bolometric nuclear emission is very sub-Eddington (l ≃ 10⁻⁴), within the range observed for l, though the most frequently observed values are somewhat lower (l ≈ 10⁻⁵). Unfortunately, uncertainties in the bolometric correction which apply to observed nuclear luminosities, appropriate for a spectral energy distribution at low emission levels, do not allow us to make a stringent comparison between modeled and observed values. However, the nuclear X-ray emission L_{BH, X} of the models, estimated as L_{BH, X} ≲ 0.2L_BH (as appropriate for low-luminosity nuclei) and L_{BH, X}/L_Edd, also tend to lie on the upper end of that observed. Thus in real galaxies an additional mechanism may be at work to further reduce the mass available for accretion; this could be provided by the mechanical feedback of a (nuclear) jet, and/or by a wind from an RIAF. Alternatively, the switch from disk to RIAF behavior (that is $l\propto \dot{m}^2$) should occur at a larger l than assumed here (e.g., at l ≃ 0.1 rather than l ≃ 0.01). Ram-pressure stripping effects, however, cannot produce a reduction in the accretion rate because in the quiescent hot accretion regime the accreting mass comes mainly from the stellar mass losses within the central ∼100 pc of the galaxy, which are not affected by stripping (Shin et al. 2010b).
• 2.  
  The X-ray luminosity of the ISM oscillates in phase with the nuclear luminosity, though with much broader peaks; at the present epoch, L_X lies at the lower end of the large observed range for galaxies of L_B similar to that of the models. The degassing/heating in the models may then be too efficient, or a larger/more concentrated gravitating mass, or a confining external medium, are needed. However, when the gas luminosities during the whole evolution are considered, the observed L_X range is better reproduced. This is even truer if an additional model of slightly larger σ is also included. Part of the observed large variation in L_X for galaxies of a given L_B could then be contributed by nuclear activity.
• 3.  
  The average ISM temperature is within the observed large range for the model σ; when estimated within R_e, the model temperatures better reproduce the upper half of those observed. Modifications to the models such as the addition of a jet or an external medium, as suggested in the previous points, should then not increase the average temperature within R_e to be viable.
• 4.  
  During quiescence, the profiles of the gas temperature and brightness resemble those observed for many local galaxies. Especially remarkable is the lack of the steep brightness profile shape typical of inflowing models due to the frequent removal of gas from the galactic central regions and to the heating provided by mechanical feedback (which is always present, even during quiescent phases, though at a very low level). The models show negative temperature gradients, as is common for isolated galaxies; the addition of a jet or a confining agent should change the temperature profiles into a flat or outwardly increasing profile, as also frequently observed for galaxies in rich environments.
• 5.  
  During outbursts, disturbances are predicted in the temperature and brightness profiles; the ISM resumes the smooth appearance of steady and low-luminosity hot accretion onto the SMBH on a timescale of ≲ 200 Myr. The most conspicuous variations with respect to smooth profiles are within current detection capabilities and could correspond to (part of) the widespread disturbances observed in galaxies of the local universe. In particular, shocked hot gas should be seen at the galactic center (within ∼100 pc), possibly not resolved; this would be a certain sign of prior AGN activity. These hot bubbles could be revealed by emission of cosmic rays in a structure similar to a gigantic supernova remnant. Preliminary results from 2D runs show bipolar nuclear outflows, which should be seen as conical cavities extending from the galactic center, and may be called jet-like features.
• 6.  
  The duty cycle of nuclear activity is of the order of a few × (10⁻³ to 10⁻²), depending on the assumed mechanical feedback efficiency; in general, a burst cycle lasts for ≈10⁷ yr. These duty-cycle values are broadly consistent with the fraction of active galaxies measured in observational works, though reported values for the local universe are somewhat lower for the SMBH mass of the models. In order to make a more consistent comparison with observations, the data set of models should be increased and the duty cycle computed only over the past 2–3 Gyr; this will reduce the duty cycle, as the models are characterized by a declining nuclear activity.
• 7.  
  The duty cycle of perturbances in the ISM is of the order of 5%–10% from their average number and duration for galaxies of L_B similar to that of the model. This duty cycle likely increases with galaxy mass because an outburst has a greater impact in less massive (and less gas-rich) systems, which are then "on" for a shorter time (Ciotti & Ostriker 2012). The ISM duty cycle and its trend with galaxy mass compare reasonably well with preliminary estimates obtained from a large sample of hot gas coronae in elliptical galaxies observed with Chandra (Nulsen et al. 2009): the fraction of galaxies with X-ray cavities in the hot gas is ≲ 10% when L_X < 10⁴¹ erg s⁻¹ (as for the models) and reaches ∼25% in the most luminous ones.
• 8.  
  Two diagnostic planes have been constructed. In the first one, the nuclear luminosity L_BH and the ISM luminosity L_X are followed during the whole model evolution. The models populate a wedge region, which should then be occupied when observing a large set of galaxies. The other plane shows the evolution of L_X versus the average gas temperature T_a; here the most populated region is that of a large L_X variation (factor of ∼10) for T_a keeping between 0.4 and 0.6 keV.

Clearly, a larger set of models is to be explored in order to better establish the final gas properties (such as gas content, nuclear and ISM duty cycle, etc.) to be compared with those of a statistically large sample. For example, a general expectation is that changes of the galaxy properties have an impact on the number of nuclear outbursts: depending on many model parameters (supernova rate, central σ, dark matter amount and distribution, and even external pressure due to an intragroup or intracluster medium), bursts could take place even toward the present epoch or be confined to the early epoch (Ciotti & Ostriker 2012).

L.C. and S.P. are supported by the grant MIUR PRIN2008.