GRAIN PHYSICS AND INFRARED DUST EMISSION IN ACTIVE GALACTIC NUCLEUS ENVIRONMENTS

It is well-established that the brightest active galactic nuclei (AGNs) and the most massive black holes reside in giant elliptical galaxies (Dunlop 2004). It is also accepted that the environment provided by the galaxy, particularly the abundance of hot metal-rich gas, aging stars, and heavily depleted dust, profoundly affects the supermassive black hole (hereafter SMBH) evolution, and in turn the host system is modified by radiative and mechanical feedback (Begelman et al. 1984; Norman & Scoville 1988; Ciotti & Ostriker 1997, 2012; Ostriker et al. 2010; Pellegrini et al. 2012; Vogelsberger et al. 2013, and references therein). Indeed, the interplay between the evolution of the black hole and the evolution of its host galaxy has been frequently invoked in order to explain the observed correlation between the properties of the black hole and the galaxy (e.g., Springel et al. 2005; Ostriker & Ciotti 2005; Cattaneo et al. 2009; Debuhr et al. 2011). Putting the questions we seek to answer in concrete terms: if the black hole in a galaxy such as M87 were to erupt as a quasar, what would we see? What mechanisms would govern the interaction between the accreting gas of the galaxy and the erupting black hole? What sources will contribute to the observed radiation?

A key ingredient in answering each of the above questions is the presence of dust in the galaxy. The potential importance of dust can be readily seen for several reasons. First, quasars emit at or near their Eddington limit, which classically is set by the Thomson scattering opacity of 0.4 cm² g⁻¹. By comparison, the opacity of dust grains to UV photons exceeds 1000 times this value, so the radiative forces on dust can easily overwhelm other dynamical drivers during outbursts. Second, it is widely accepted that a significant fraction of quasars are obscured. However, the giant elliptical galaxies which host AGN typically have very small optical depths, implying little dust. Therefore, not only does the dust govern the observed radiation from the galaxy, but the dust abundance itself may be a dynamic quantity that evolves in parallel with quasar "on" and "off" phases. In this work, we seek to understand the power of dust to influence the nature, evolution, and observational signatures of giant elliptical galaxies by introducing the processes that create and destroy dust within galaxies into our simulations.

One of the primary links between an AGN and its host galaxy is the immense radiative output of the AGN. The radiative feedback mechanisms between the central black hole and the accreting gas greatly influence the gas dynamics and thus have been the focus of much study (Binney & Tabor 1995; Ciotti & Ostriker 1997; Thompson et al. 2005; Ciotti & Ostriker 2007, hereafter CO07). Because dust grains efficiently absorb UV radiation and re-radiate in the IR, the presence of grains can dramatically affect the gas dynamics (Siebenmorgen & Heymann 2012, and references therein). In particular, infalling shells of cold gas (a common feature at the onset of interstellar medium (ISM) cooling episodes and detected in giant ellipticals Werner et al. 2014) can be supported by radiation pressure, slowing accretion, and altering the subsequent evolution. In turn, this evolution is related to the so-called "positive feedback" mode in which feedback from the AGN enhances star formation (CO07, Nayakshin & Zubovas 2012, and references therein). Observations have hinted at the presence of radiation-supported cold gas around AGN (Vasudevan et al. 2013) though in the majority (≃ 95%) of cases the gas is below its effective Eddington limit (Raimundo et al. 2010). Nevertheless, the latter authors note the inevitability of radiation-supported gas and its possible role in AGN feedback. Radiation pressure on dust grains may be instrumental in driving galactic winds (Coker et al. 2013), though there is disagreement as to whether this could be a substantial effect (Socrates & Sironi 2013).

In a recent series of papers (Ciotti et al. 2009b; Shin et al. 2010; Ciotti et al. 2010), one-dimensional (1D) hydrodynamical simulations were used to study the mechanisms of AGN feedback. The dust physics was implemented in a very simple, phenomenological way. As sputtering effectively destroys dust in regions of hot gas (Draine & Salpeter 1979), these previous papers approximated the sputtering by reducing the dust to gas ratio in hot gas by roughly two orders of magnitude relative to Galactic values in accord with observations of dust in elliptical galaxies. In practice, a factor inversely proportional to the gas temperature multiplies the fiducial absorption coefficient in the UV, optical, and IR.

In this paper, we significantly improve this aspect of the input physics. We seek to clarify the role of dust grains in AGN feedback by developing and numerically implementing the basic equations for dust production and destruction. While a similar but simpler version of the dust physics described here was implemented in the two-dimensional (2D) simulations of Novak et al. (2012), we present a more generalized treatment and focus principally on the effects of different models of dust abundance on the hydrodynamical evolution of the simulation, both in its galaxy-scale properties and accretion physics. Additionally, we discuss the ability or inability of models to reproduce the observational properties of giant elliptical galaxies in the infrared during both quiescent and accretion phases.

For the present hydrodynamical simulations we use an updated version of the 1D code used in our previous studies, with the relevant modifications detailed in Section 3. In particular, a major upgrade has been made in the solution of the radiative transfer equation. We note that in a recent related paper (Novak et al. 2012) the new "two-streams" approach is also tested in the context of 2D hydrodynamical simulations.

We restrict our focus to the secular evolution of the galaxy. While processes such as major and minor mergers or inflows of cold gas from the intergalactic medium (IGM) certainly occur and can account for a variety of observations, we are most interested in which aspects of the evolution of these galaxies can be accounted for by purely secular processes, i.e., driven by well-understood stellar evolutionary processes. We find that the AGN duty cycle, its time dependence, the black hole mass, and the IR luminosity are all naturally explained with our purely secular model.

The paper is organized as follows: in Section 2 we describe the dust grain physics employed in our simulations as well as a discussion of two-stream radiative transfer. In Section 3 we describe the suite of simulations conducted under different models for the grain abundance. In Section 4, we summarize the results of the simulations and compare against observations to select the most viable models. We also discuss the effects of dust on the properties of AGN outbursts. We then discuss the most viable models in the context of observations in Section 5, focusing on both observational signatures predicted by our model as well as how current observations constrain the survival, growth, and radiative properties of dust in elliptical galaxies with AGN. We summarize the implications of these comparisons in Section 6.

In this section we summarize the dust physics implemented in the new version of the code and used for the hydrodynamic simulations.

To aid comparisons between models, we define the factor D to describe how much the dust is depleted relative to what is observed in the Galaxy, i.e.,

$$\begin{equation} \ {D} \equiv \frac{Z_{\rm MW}}{Z} \left(\frac{\rho _d}{\rho }\right) \left(\frac{\rho _d}{\rho }\right)_{\rm MW}^{-1}, \end{equation} \tag{ 1 }$$

where ρ and ρ_d are the gas and dust mass densities, respectively, and Z and Z_MW are the metallicities of the simulated galaxy and the Milky Way, respectively. This is motivated by the fact that the dust to gas ratio of a galaxy should scale linearly with the metallicity. A depletion factor of one therefore indicates that the fraction of metals incorporated in dust grains is the same as in the Milky Way. We adopt Z_MW = Z_☉, a Galactic dust to gas ratio of 0.01, and Z/Z_MW = 4/3 for all times in the simulations.

2.1. Grain Opacity

Following CO07, we adopt the dust opacity values

$$\begin{equation} \ \kappa _{\rm Op} = 300 \frac{{\rm cm}^2}{\rm g} \left(\frac{Z}{Z_{\rm MW}}\right) \times {\rm D}, \quad \kappa _{\rm UV} = 4 \kappa _{\rm Op}, \quad \kappa _{\rm IR} = \frac{\kappa _{\rm Op}}{150}, \end{equation} \tag{ 2 }$$

i.e., dependent on the dust to gas ratio in the galaxy relative to the Milky Way. The bands correspond to the bands used by Sazonov et al. (2004) to compute the broadband AGN output and are listed in Table 1. Note that in CO07,

$$\begin{equation} \ {D} = \frac{1}{1+T_4}, \end{equation} \tag{ 3 }$$

where T₄ is the gas temperature in units of 10⁴ K and the metallicity taken to be Z_☉. For coronal X-ray emitting gas in a typical elliptical galaxy, this corresponds to a depletion factor of ≃ 10⁻². In the following, we present a method for computing D in a much more physically consistent way than Equation (3).

Table 1. Bandpasses

Band	Energy	AGN Output Fraction
X-ray	E > 2 keV	0.1
UV	13 <E < 2 keV	0.35
Op	1 <E < 13 eV	0.25
IR	E < 1 eV	0.3

Note. The assumed energy limits for the bandpasses used in this work.

Download table as:  ASCIITypeset image

2.2. Dust Abundance

The primary sources of dust in giant elliptical galaxies include the massive winds ejected by dying asymptotic giant branch (AGB) stars during thermal pulses (for example as planetary nebulae) and by supernova explosions, as well as grain growth in the metal-rich ISM. The ultimate fate of these winds is not fully understood, but it is commonly accepted that some form of mixing and thermalization with the preexisting ISM takes place on short timescales (e.g., Bregman & Parriott 2009). Therefore, it is natural to expect that the dust is transported through the galaxy by the large scale gas flows of early type galaxies (Kim & Pellegrini 2012), where it is sputtered by the hot ISM and, in cold gas, allowed to grow via collisions with metal atoms. Moreover, when dust-bearing gas forms new stars, the dust is removed from the ISM and sequestered into the stars.

In the previous simulations discussed in the Introduction, these effects were addressed qualitatively by assuming that the cold gas was grain-rich and the hot gas was grain-poor. In practice, this was implemented by approximating the dust to gas ratio at each radius by Equation (3). However, this prescription changes the dust abundance instantaneously with temperature and does not allow for the transport of grains from one radius to another.

2.2.1. One Component Model

The new treatment addresses both issues by implementing a grain continuity equation

$$\begin{equation} \ \frac{\partial \rho _{d}}{\partial t} + \nabla \cdot (\rho _d {\bf v}) = S_+ - S_-, \end{equation} \tag{ 4 }$$

where ρ_d is the mass density of dust grains. Note that v is the gas velocity given by the hydrodynamic code since we assume that the grains are coupled to the motion of the gas. Although drift relative to the gas will occur due to gravitational forces and anisotropic radiation fields, such drift will be mitigated by Coulomb drag and drag from the local magnetic field since the grains will be charged. The net drifts will in general be small relative to the Eulerian velocities. Finally S₊ and S₋ are the dust source and sink terms, respectively. As discussed below, each of the two functions S₊ and S₋ is in turn given by the sum of two terms. We term this the "One Component" model as all dust in this model is assumed to be mixed with the ISM. We describe a "Two Component" model in Section 2.2.2 that considers mixed and unmixed dust separately.

In the Milky Way, the dust to gas ratio in outflows of oxygen-rich AGB stars has been observed to be ≃ 0.0063 (Knapp 1985; Kemper et al. 2003). Since an early type galaxy can be metal-enriched relative to the Milky Way, we adopt a fiducial dust to gas mass ratio of 0.01 in these outflows, and the source term is simply

$$\begin{equation} \ S_{+, {\rm inj}} = 0.0133 \dot{\rho }_{*}, \end{equation} \tag{ 5 }$$

where $\dot{\rho }_{*}$ is the gas released by stars as described in CO07 Equation (13), and the numerical factor is dimensionless.

In the particularly hot ISM of elliptical galaxies, grains are rapidly sputtered. To compute the dust destruction rate due to sputtering, we use the relation

$$\begin{equation} \ \dot{a} = -\frac{10^{-6} n_{\rm H}}{1 + T_6^{-3}}\ \mu {\rm m}\ {\rm yr}^{-1} \end{equation} \tag{ 6 }$$

where a is the grain radius, T₆ is the gas temperature in units of 10⁶ K, and n_H is the proton number density in units of cm⁻³. This expression is a good approximation for graphite and silicate grains in gas with 10⁵ < T < 10⁹ K (Draine 2011); note that $\dot{a}$ is independent of a. Empirically, the size distribution of dust grains above ≃ 50 Å can be approximated by the standard Mathis–Rumpl–Nordsieck (MRN) distribution (Mathis et al. 1977) despite destruction and creation processes. Therefore, we assume that the MRN distribution is valid at all times, i.e.,

$$\begin{equation} \ \frac{{\rm d}n_d}{{\rm d}a} = \frac{H a_{\rm max}^{2.5}}{a^{3.5}} \end{equation} \tag{ 7 }$$

is the number density of grains with size between a and a +da, where H is a normalization constant, in principle dependent on time and position in the galaxy. In the following, we assume a_min = 0.005 μm and a_max = 0.3 μm. The total density in grains at a given radius and time is then

$$\begin{equation} \rho _{d} \simeq \frac{8 \pi }{3} H \rho _{{\rm grain}} a_{{\rm max}}^3, \end{equation} \tag{ 8 }$$

where ρ_grain is the internal density of a dust grain and we neglected the factor of $1 - \sqrt{a_{{\rm min}}/a_{{\rm max}}}$. We take ρ_grain = 3.5 g cm⁻³ which is a standard value for silicate grains. The total destruction rate is obtained by computing the mass destruction rate of grains of radius a and then integrating over the distribution. It follows that

$$\begin{equation} \ \dot{\rho }_{d,gd} = 8 \pi H \rho _{\rm grain} \left(\frac{a_{\rm max}}{a_{\rm min}}\right)^{0.5}a_{\rm max}^2\dot{a}, \end{equation} \tag{ 9 }$$

and from Equations (8) and (9), we can define a grain destruction frequency

$$\begin{equation} \ \nu _{gd} \equiv \frac{|\dot{\rho }_{d,gd}|}{\rho _{d}} = 3 \left(\frac{a_{{\rm max}}}{a_{\rm min}}\right)^{0.5} \frac{|\dot{a}|}{a_{\rm max}}. \end{equation} \tag{ 10 }$$

Therefore, the sputtering term in Equation (4) is

$$\begin{equation} \ S_{-,gd} = \nu _{gd} \rho _d. \end{equation} \tag{ 11 }$$

Note that in principle, in AGN environments, where high energy photons can ionize grains, the sputtering time can be altered by the effects of grain charging. However, Weingartner et al. (2006) find that the effect at r = 100 pc for an L_BH = 10⁴⁶ erg s⁻¹ quasar is negligible above ≃ 10⁶ K, even for large ionization parameters. Thus, this effect may be safely neglected.

In cold gas, metal atoms are able to collide with dust grains and stick. If these metal atoms have a probability f of sticking and are moving with average speed v_Z, then the source term due to these collisions is fρ_Zv_Z4πa²n_d, where ρ_Z = f_Zρ is the mass density of the metals and f_Z is the mass fraction of gas available for making grains. If we assume that the Milky Way has included all such materials in grains already, then f_Z for the Milky Way would just be its observed dust to gas ratio of 0.01. Since the simulated galaxy has a metallicity 4/3 greater than the Milky Way, we take f_Z to be 0.0133. However, grain growth cannot continue after the metals in the gas have been used up, so we replace f_Z with f_Z − ρ_d/ρ to disallow growth beyond the available metal atoms. Integrating this over the grain size distribution, we obtain

$$\begin{equation} \ \dot{\rho }_{d,gg} = f \rho \frac{c_s}{\sqrt{\mu _Z}} 8 \pi \left(f_Z - \frac{\rho _{d}}{\rho } \right) H a_{\rm max}^2 \left(\frac{a_{{\rm max}}}{a_{\rm min}}\right)^{0.5}, \end{equation} \tag{ 12 }$$

where we have made the approximation $v_Z = c_s/\sqrt{\mu _Z}$, with c_s being the sound speed in the gas and μ_Z the mean atomic mass of the metal atoms. Thus the grain growth frequency is

$$\begin{equation} \ \nu _{gg} \equiv \frac{\dot{\rho }_{d,gg}}{\rho _d} = \frac{3 f \rho c_s}{\rho _{\rm grain} a_{\rm max}\sqrt{\mu _Z}} \left(f_Z - \frac{\rho _d}{\rho } \right) \left(\frac{a_{\rm max}}{a_{\rm min}}\right)^{0.5}. \end{equation} \tag{ 13 }$$

Since grains are made primarily from carbon, oxygen, magnesium, silicon, and iron, we adopt μ_Z = 16 corresponding to oxygen and in agreement with the mean atomic mass in current grain models (Draine 2011, Table 23.1). Following Clayton & Wickramasinghe (1976), who modeled grain growth in the 10⁴ K gas of an expanding nova shell, we take f = 0.2, but note that there is considerable uncertainty in the surface chemistry of grains, particularly in the high temperature environment of an elliptical galaxy.

The ratio of grain growth by collisions to sputtering is a temperature dependent function given by

$$\begin{equation} \ \frac{\nu _{gg}}{\nu _{gd}} = 1.18\times 10^{-6} \left(1 + T_6^{-3}\right) \sqrt{T} \end{equation} \tag{ 14 }$$

for our assumed galaxy parameters and approximating f_Z ≫ ρ_d/ρ. Figure 1 plots the growth and destruction times as a function of temperature for gas with n_H = 0.1 cm⁻³. Grain growth is negligible in all but the coldest gas.

Zoom In Zoom Out Reset image size

Finally, grains will be removed from the ISM when the gas containing those grains form stars. The dust destruction due to star formation is

$$\begin{equation} \ S_{-,{\rm SF}} = \frac{\rho _d}{\rho _g} \dot{\rho }_*^+, \end{equation} \tag{ 15 }$$

where $\dot{\rho }_*^+$ is the star formation rate (SFR) at radius r as described in Section 3.

Combining the equations in this section, we now have the expression for the dust source and sink terms:

$$\begin{equation} \ S_+ - S_- = S_{+, inj} + S_{+, gg} - (S_{-, gd} + S_{-, SF}). \end{equation} \tag{ 16 }$$

By evolving Equation (4) for each position and time, we can evaluate the D function in Equation (1) with the computed value of ρ_d.

2.2.2. Two Component Model

The one component approach assumes that dust mixes with gas instantaneously after its creation in stellar outflows. However, stellar ejecta contains dust not yet mixed with the ambient ISM. Although dust-laden, these compact planetary nebulae will not contribute significantly to the total infrared opacity of the galaxy. Therefore, the mixing process introduces a lag time between dust creation and its consequent effects on the radiative feedback in the galaxy.

To estimate the mixing time, we consider a planetary nebula of mass ΔM expanding with velocity v₁ relative to the parent star into a surrounding medium of density ρ_ext and sound speed c_s. We further define v_rel as an estimate of the relative velocity of the star with respect to the surrounding ISM, so that a fiducial value is the local (1D) velocity dispersion of the stars. We denote the internal density as ρ₁ and the radius of the nebula as r₁. The nebula expands until it comes into pressure equilibrium with the surrounding medium, i.e.,

$$\begin{equation} \ v_1^2 \rho _1 = P_{{\rm ext}}. \end{equation} \tag{ 17 }$$

Hence the time to reach pressure equilibrium t_eq is

$$\begin{equation} \ t_{eq} = \left(\frac{3 \Delta M}{4 \pi v_1 P_{{\rm ext}}}\right)^{1/3}. \end{equation} \tag{ 18 }$$

Defining the Mach number ${\cal M} \equiv ({v_{{\rm rel}}}/{c_s})$, the ram pressure experienced by the expanding wind is

$$\begin{equation} \ P_{{\rm ext}} = \rho _{{\rm ext}} \left(v_{{\rm rel}}^2 + c_s^2\right) = v_{{\rm rel}}^2 \rho _{{\rm ext}} (1 + {\cal M}^{-2}). \end{equation} \tag{ 19 }$$

Let t_fr be the time it takes the planetary nebula to encounter a gas mass equal to its own mass, thereby fragmenting and mixing it. Then,

$$\begin{equation} \ \frac{\Delta M}{\pi r_1^2} = t_{fr} \rho _{{\rm ext}}\ {\rm max}(v_1, v_{{\rm rel}}), \end{equation} \tag{ 20 }$$

where typically v_rel ≫ v₁. Using Equation (19), t_fr can be expressed as

$$\begin{equation} \ t_{fr} = \left(\frac{16}{9\pi }\right)^{1/3} (1 + {\cal M}^{-2})^{2/3} \left(\frac{v_{{\rm rel}}}{v_1}\right)^{1/3} \left(\frac{\Delta M}{v_1^3 \rho _{{\rm ext}}}\right)^{1/3}. \end{equation} \tag{ 21 }$$

In addition to fragmentation, we must also consider evaporation due to thermal conduction. Following Draine (2011) Equation (34.17),

$$\begin{equation} \ t_{ev} = 1.6\times 10^3\ {\rm yr} \left(\frac{\Delta M}{1\ M_\odot }\right) \left(\frac{r_{ev}}{1\ {\rm pc}}\right)^{-1} \left(\frac{T}{10^7\ {\rm K}}\right)^{-2.5}, \end{equation} \tag{ 22 }$$

where T is the gas temperature and r_ev is the radius of the nebula when it evaporates. Note that we have assumed that the Coulomb logarithm ln Λ = 30. In our implementation of these equations, ΔM = 0.1 M_☉, and v₁ = 10 km s⁻¹.

There are three distinct regimes to consider. First, consider the case in which t_ev < t_eq, i.e., the nebula evaporates before it expands to pressure equilibrium. Then the radius r_ev = v₁t_ev. Solving for t_ev, which in this case is the mixing time, we obtain

$$\begin{equation} t_{\rm mix,1} = 8.4\times 10^3\ {\rm yr}\ \left(\frac{T}{10^7\ K}\right)^{-5/4}. \end{equation} \tag{ 23 }$$

However, if t_eq < t_ev, then r_ev = r₁. Once at pressure equilibrium, the nebula can mix via fragmentation or evaporation, depending on which is faster. In the former case,

$$\begin{eqnarray} \ t_{\rm mix,2} &=& 1.5\times 10^3\ {\rm yr}\ \left(\frac{T}{10^7\ K}\right)^{-5/2} \left(\frac{v_{\rm rel}}{300\ {\rm km {\rm \,s}^{-1}}}\right)^{2/3} \nonumber \\ && \times (1 + {\cal M}^{-2})^{1/3}. \end{eqnarray} \tag{ 24 }$$

However, if t_{mix, 2} > t_fr, we are in the third regime where the mixing time t_{mix, 3} = t_fr:

$$\begin{eqnarray} t_{\rm mix, 3} &=& 1.9\times 10^6\ {\rm yr} \left(\frac{v_{\rm rel}}{300\ {\rm km\,{\rm s}^{-1}}}\right)^{1/3} \left(\frac{n_{\rm H}}{0.01\ {\rm cm}^{-3}}\right)^{-1/3} \nonumber \\ && \times (1 + {\cal M}^{-2})^{2/3}. \end{eqnarray} \tag{ 25 }$$

In summary, the mixing time is defined as

$$\begin{equation} t_{\rm mix} = \left\lbrace \begin{array}{@{}lr@{}}\displaystyle t_{\rm mix, 1}, &t_{ev} < t_{eq};\\ \displaystyle {\rm min}(t_{\rm mix, 2}, t_{\rm mix, 3}) &t_{ev} > t_{eq}. \end{array} \right. \end{equation} \tag{ 26 }$$

Equipped with a mixing time, we may now modify the continuity equations for gas and dust by distinguishing between the planetary nebula (PN) and diffuse ISM phases. We assume that no dust growth or destruction occurs in the PN phase.

$$\begin{equation} \frac{\partial \rho _{g,{\rm ISM}}}{\partial t} + \nabla \cdot (\rho _g {v}) = \dot{\rho }_{\rm II} - \dot{\rho }_*^+ + \dot{\rho }_{\rm w} + \frac{\rho _{g,PN}}{t_{\rm mix}} \end{equation} \tag{ 27 }$$

$$\begin{equation} \frac{\partial \rho _{g,PN}}{\partial t} = \dot{\rho }_* - \frac{\rho _{g,PN}}{t_{\rm mix}} \end{equation} \tag{ 28 }$$

$$\begin{eqnarray} \frac{\partial \rho _{d,{\rm ISM}}}{\partial t} + \nabla \cdot (\rho _d {v}) &=& 0.0133 \dot{\rho }_{\rm II} + \nu _{gg} \rho _{d,{\rm ISM}} - \nu _{gd} \rho _{d,{\rm ISM}} \nonumber \\ &&-\frac{\rho _{d,{\rm ISM}}}{\rho _{g,{\rm ISM}}} \dot{\rho }_*^+ + \frac{\rho _{d,PN}}{t_{\rm mix}} \end{eqnarray} \tag{ 29 }$$

$$\begin{equation} \frac{\partial \rho _{d,PN}}{\partial t} = 0.0133 \dot{\rho }_* - \frac{\rho _{d,PN}}{t_{\rm mix}}, \end{equation} \tag{ 30 }$$

where $\dot{\rho }_{\rm II}$ is the gas source term associated with the young stellar population via Type II supernovae and $\dot{\rho }_{\rm w}$ is the source term associated with winds from the circumnuclear disk. In the limit of small but constant t_mix, we recover the gas continuity equation of Ciotti & Ostriker (2012), Equation (4.73).

2.3. Grain Temperature

The physics presented in the previous section directly influences the hydrodynamical evolution of the models. Here we present additional physics needed to compute observational properties of the models, the other focus of this work.

By numerical integration of the radiative transfer equations in the one-stream approximation (see Section 2.4), we compute the total radiation density in each shell. We approximate the total radiation absorbed in each radial shell by dust as the difference in luminosity at the base and end of the shell, i.e.,

$$\begin{equation} \ \Delta L = L_{\rm eff, inner} - L_{\rm eff, outer}, \end{equation} \tag{ 31 }$$

where the effective luminosities include contributions from both optical and UV bands. We assume this luminosity is radiated by a population of grains in that shell, all with steady-state temperature T_d, since the steady-state temperature of a dust grain is nearly size-independent (see, e.g., Draine 2011, Equations (24.19) and (24.20)). Imposing that the total dust emission in a shell of volume V is equal to the computed ΔL, the relation between T_d and ΔL in a given shell is:

$$\begin{equation} \ \Delta L = V \sigma T_d^4 \int _{a_{{\rm min}}}^{a_{{\rm max}}} \! {\rm d}a \frac{{\rm d}n_d}{{\rm d}a} 4 \pi a^2 Q(a, T_d), \end{equation} \tag{ 32 }$$

where Q(a, T_d) is the Planck-averaged emission efficiency. Following the power-law prescription for the silicate Planck-averaged emission efficiency of Draine (2011) Equation (24.15) at low T_d and approximating the high T_d behavior of Q(a, T_d)/a as a constant, we obtain:

$$\begin{equation} Q(a,T_d) = \left\lbrace \begin{array}{@{}lr@{}}\displaystyle \left(\frac{a}{0.1 \,\mu {\rm m}}\right) 1.3 \times 10^{-6}\ T_d^2, & \, T_d < 164; \\ \displaystyle \left(\frac{a}{0.1 \,\mu {\rm m}}\right) 3.5 \times 10^{-2}, & \, T_d > 164; \end{array} \right. \end{equation} \tag{ 33 }$$

where T_d is in Kelvin. We note that the steady-state temperature for graphitic grains does not differ substantially from that of silicate grains (Draine 2011, Equation (24.20)), allowing us to focus on silicates for specificity and simplicity. By inserting Equation (33) into Equation (32), some algebra shows that the equilibrium dust temperature is given by

$$\begin{equation} T_d = \left\lbrace \begin{array}{@{}lr@{}}\displaystyle 5.97 \left(\frac{\rho _{{\rm grain}} \Delta L}{\rho _d V}\right)^\frac{1}{6},& T_d < 164; \\ \displaystyle 1.14 \left(\frac{\rho _{{\rm grain}} \Delta L}{\rho _d V}\right)^\frac{1}{4} ,& T_d > 164; \end{array} \right. \end{equation} \tag{ 34 }$$

where all quantities are in cgs. Note that the integral in Equation (32) is monotonic in T_d, therefore ensuring only one branch of Equation (34) is selected for given input values of ΔL, V, and ρ_d. Because this calculation neglects the stochastic heating of small grains, the derived grain temperature should be considered as a characteristic temperature for the far infrared (FIR) dust emission.

Finally, we define the luminosity-weighted dust temperature 〈T_d〉 to be

$$\begin{equation} \langle T_d\rangle = \frac{1}{L_{\rm IR}} \int T_d(r) \left(\frac{\Delta L (r)}{\Delta r}\right) {\rm d}r, \end{equation} \tag{ 35 }$$

which provides an estimate of the temperature of the dust producing the observed IR emission.

2.4. One-Stream Radiative Transfer in Spherical Symmetry

The integration scheme for the radiative transfer equation has been improved with respect to Ciotti et al. (2010). The full description of the new scheme is given in Novak et al. (2012). In particular, we adopt the Simplified Radiation Transport in their Appendix B, which we describe briefly below. In Novak et al. (2012) the full equations of radiative transfer were solved by using a relaxation method, and it was shown that the following approximation works remarkably well for the present problem.

The radiation transport equations for the black hole radiation are particularly simple because all UV and optical photons emitted from the black hole will necessarily be outgoing photons assuming that the scattering opacity is negligible. This yields the relation

$$\begin{equation} \ \frac{dL_{\rm eff, BH}}{dr} = - \rho \kappa _i L_{\rm eff, BH}, \end{equation} \tag{ 36 }$$

where κ_i is the dust opacity in band i and the effective black hole luminosity L_{eff, BH} is the outgoing black hole luminosity that would be seen by an observer at radius r, i.e., after absorption.

We make the approximation that all absorption in the UV is due to dust. The photoionization opacity of the gas competes with the dust when the neutral fraction is above ∼10⁻³ for Galactic dust-to-gas ratios. However, in 10⁷ K gas, the neutral fraction is of order 10⁻⁸ due to collisional ionization alone (see, e.g., Draine 2011, Equations (14.39) and (14.43)). During AGN "on" phases, a photoionizing luminosity of 10⁴⁶ erg s⁻¹ from the AGN is able to maintain a steady-state neutral fraction of ∼10⁻⁷ in a dense (n_H = 10³ cm⁻³) cloud with T = 10⁴ K at r = 100 pc. In both cases, the gas opacity is negligible compared to the dust.

Following Sazonov et al. (2004), we assume the AGN radiates 10% of its energy in the X-ray, 35% in the UV, 25% in the optical, and 30% in the IR. The IR value includes contribution from a subgrid dusty accretion torus which we leave in place even in our "No Dust" model.

For other radiation sources, however, the equation is complicated by the fact that photons may be emitted inward toward the center of the galaxy and, providing the optical depth is low enough, re-emerge on the other side as an outgoing photon. To account for this, we approximate the probability that an emitted photon will be outgoing, either initially or by re-emerging to the same radius on the other side, by the function Ψ which is defined as

$$\begin{equation} \ \Psi \equiv 1 - \frac{0.5}{1+{\rm exp}(-\tau)}\frac{r_1^2}{{\rm max}\left(r_1^2,r^2\right)}, \end{equation} \tag{ 37 }$$

where τ is the optical depth from r to infinity and r₁ is the radius at which τ = 1. Using this parameterization, the radiative transport equation for the outgoing stellar radiation L_{eff, *} is given by

$$\begin{equation} \ \frac{dL_{\rm eff, *}}{dr} = 4\pi r^2 \Psi \dot{E}_{i} - \rho \kappa _i L_{\rm eff, *}, \end{equation} \tag{ 38 }$$

where $\dot{E}_{i}$ is energy radiated by stars per unit volume per unit time at radius r in band i. In practice, $\dot{E}_i$ is dependent upon time, radius, and the local gas density and is partitioned into the optical and UV bands through use of characteristic emission efficiencies and timescales for each band (Ciotti & Ostriker 2012, Equations (4.24) and (4.25)).

The treatment of the radiation pressure has remained unchanged from CO07 other than the method of computing the dust opacity, and includes radiation pressure on gas from electron scattering and X-ray photoionization as well as the radiation pressure on dust.

In addition to the changes detailed in Novak et al. (2012), we also include an updated prescription for the optical depth of a radial shell. If a shell is optically thick, only a portion of the shell will experience a force from the radiation pressure. To achieve the proper limiting behavior, we modify the optical depth in a band i of a given shell $\tau _i^{\prime }$ in the following way to obtain a τ_i for use in calculations of effective luminosities and radiation pressure:

$$\begin{equation} \ \tau _i \equiv 1 - {\rm e}^{-\tau _i^{\prime }}. \end{equation} \tag{ 39 }$$

We present a suite of 1D hydrodynamical simulations of the coevolution of a giant elliptical galaxy and its central supermassive black hole. The simulation begins after the initial starburst that produced the majority of the galaxy's stellar mass, leaving the galaxy with no remaining gas. Cooling flow instabilities in the secondary gas from stellar evolution primarily drive accretion onto the central SMBH, which leads to the production of nuclear and galactic winds. Mechanical and radiative feedback from the AGN, Type Ia and Type II supernovae, stellar radiation, and thermalization from stellar mass losses are all explicitly considered.

For simplicity, we restrict our simulations to the class of Type A models described in Ciotti et al. (2010). In these models, the opening angle of the broad line region (BLR) wind and the mechanical efficiency _w are independent of the accretion luminosity. _w = 10⁻⁴ is a factor of two below the mechanical efficiency assumed in many of the treatments of AGN feedback, such as Di Matteo et al. (2005), but similar to the value found most appropriate when winds are included (see, e.g., Choi et al. 2013). We note that including the momentum of the outgoing wind makes a given energy input far more effective (Choi et al. 2012). We choose to use the A class for the purpose of this study as the accretion physics is cleaner than the more intricate B class of models, whose efficiency increases with increasing Eddington ratio (Ciotti et al. 2009b), and the role of the dust is consequently easier to disentangle.

For ease of comparison, all of the dynamical properties relevant for the simulations is the same as in Ciotti et al. (2009b), i.e., a Jaffe stellar distribution plus a dark matter halo so that the total density profile is proportional to 1/r². The total stellar mass of 3 × 10¹¹ M_☉ and effective radius R_e = 6.9 kpc result in a central velocity dispersion 260 km s⁻¹. Dynamical properties of the model are given in Ciotti et al. (2009a). The initial mass of the central SMBH is fixed to M_BH = 10⁻³M_* as in previous papers, therefore approximately following the Magorrian relation. In practice, all of the evolutionary phases of galaxy formation leading to the establishment of the Magorrian relation are not considered. Accretion onto the BH is computed from the full hydrodynamic equations rather than assuming Bondi accretion or other approximate treatments. It is mediated by a circumnuclear accretion disk whose balance equations are integrated as subgrid physics (Ciotti & Ostriker 2012). Each simulation employs 240 cells with the innermost gridpoint at 2.5 pc and the outermost at 208 kpc. As in the previous papers, for simplicity we assume standard outflow boundary conditions at the grid outer boundary and use a dynamic time resolution based on the physical timescales in the galaxy. However, we increase the time resolution by an additional factor of 10 relative to previous work.

The treatment of the physics for the stellar component of the galaxy, including stellar evolution, Type Ia and Type II Supernovae, and star formation, as well as the hydrodynamical equations are fully described in Ciotti & Ostriker (2012), and we outline it briefly here. The SFR at a specific radius r is given by the equation

$$\begin{equation} \dot{\rho }_*^+ = \frac{\eta _{\rm form} \rho }{\tau _{\rm form}}, \end{equation} \tag{ 40 }$$

where η_form is an efficiency coefficient dependent on the local gas temperature and having typical values between 0.03 and 0.4 (Cen & Ostriker 2006) and τ_form is the maximum of the gas cooling time and the dynamical time. The gas cools via Compton cooling, bremsstrahlung, and both line and continuum cooling as estimated by the formulae given in Sazonov et al. (2005). These formulae are unmodified by the inclusion of dust.

In summary, the models are in all respect identical to previous models with the exception of a better treatment of dust, an improved numerical integration of the radiative transfer (see also Novak et al. 2012), and increased spatial and temporal resolution. However, for the same input physics and previous dust treatment, the results are nearly identical to previous ones. Our A₂ model refers to the precise implementation of the same model in Ciotti et al. (2010) as we use Equation (3) to model the dust depletion.

We introduce five variants of the A₂ model—and thus six models in all, with each variant utilizing a different prescription for the dust abundance and distribution. These models are summarized in Table 2, where they are listed in the approximate order of increasing dust to gas ratio at the end of the simulation.

Table 2. Summary of Models

Model	Depletion	ΔMw	ΔM*	Mgas	LX	$\frac{M_{\rm dust}}{M_{\rm gas}}$	$\frac{L_{\rm IR}}{L_{\rm Op*}}$	τOp	ΔMBH	eBol	eUV	eIR	Nburst	fduty
		(M☉)	(M☉)	(M☉)	(erg s−1)				(M☉)		eOp	eX
${\rm A}_2^{{\rm ND}}$	0	10.43	9.40	9.04	38.15	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	8.97	0.11	0.040 0.029	0.034 0.011	81	−2.21
A2	1/1+T4	10.42	9.26	9.14	38.26	−4.19	−4.75	−5.00	8.65	0.11	0.016 0.013	0.073 0.011	51	−2.45
${\rm A}_2^{{\rm -4}}$	0.75 × 10−2	10.29	9.38	9.91	40.49	−4.00	−3.68	−3.31	8.64	0.11	0.023 0.020	0.058 0.011	46	−2.25
${\rm A}_2^{\rm CE}$	Computed	10.32	9.36	9.83	40.47	−3.63	−3.15	−2.01	8.67	0.11	0.013 0.010	0.079 0.011	57	−2.42
${\rm A}_2^{\rm CE2}$	Computed	10.21	9.98	10.10	40.57	−2.54	−2.65	−1.93	9.01	0.12	0.007 0.005	0.092 0.012	86	−2.20
${\rm A}_2^{-3}$	0.75 × 10−1	10.24	10.06	9.38	39.03	−3.00	−3.30	−2.92	8.91	0.12	0.006 0.005	0.091 0.012	76	−2.27
${\rm A}_2^{{\rm MW}}$	1	10.34	8.19	9.80	39.48	−1.88	−2.07	−1.68	7.79	0.10	0.017 0.020	0.056 0.010	14	−3.14

Notes. The models are arranged roughly by dust content, from lowest to highest. All quantities are the values attained at the end of the simulation, which in all cases represents a quiescent giant elliptical galaxy. Depletion is the ratio of dust to metals in the model relative to the dust to metal ratio in the Galaxy (Equation (1)); ΔM_BH is the total mass accreted by the black hole; ΔM_w is the total mass ejected as a galactic wind; ΔM_* is the total mass of new stars; e_i ≡ ΔE_i/ΔM_BHc² is the total energy emitted by the black hole in band i (as seen from infinity) divided by the energy equivalent of the black hole mass growth; L_IR/L_Op* is the ratio of the IR luminosity from dust and the effective optical luminosity from stars; τ_Op is the optical depth in the optical band; L_X is the X-ray luminosity of the ISM in the galaxy; N_burst is the number of burst events; and f_duty is the fraction of time spent with L_Bol > L_Edd/30. All quantities except e_i and N_burst given as log₁₀.

Download table as:  ASCIITypeset image

In the first model, $A_2^{\rm ND}$, we consider a galaxy completely devoid of dust, i.e., ρ_d/ρ = 0 at all radii (equivalently, the depletion factor D in Equation (1) is fixed to zero).

$A_2^{\rm MW}$, with a dust to gas ratio equal to that of the Milky Way scaled to the metallicity of our galaxy (Z = 4/3 Z_MW), i.e., ρ_d/ρ = (4/3) × 10⁻², is the other extreme model. In this maximum dust model, D = 1.

We have two additional models in which the dust to gas ratio is a fixed number independent of time and position. First is $A_2^{-4}$, in which ρ_d/ρ = 10⁻⁴ at all radii (D = 0.75 × 10⁻²). This is motivated by recent Herschel observations (Smith et al. 2012) of the dust masses of 62 early type galaxies and scaled to our assumed stellar mass of 3 × 10¹¹ M_☉. In interest of spanning the viable range of dust to gas ratios, we also introduce $A_2^{-3}$ in which ρ_d/ρ = 10⁻³ (D = 0.75 × 10⁻¹).

Our most sophisticated models embody the suite of physics for grain production and destruction outlined in Section 2.2 to compute the dust mass density at each radius and the resulting dust opacity. The $A_2^{\rm CE}$ model employs the "One Component" formalism of Section 2.2.1 while $A_2^{\rm CE2}$ the "Two Component" formalism of Section 2.2.2.

We begin by comparing the overall behavior of all models in Table 2. Our purpose is two-fold: first to understand the effects of different treatments of the dust to gas ratio. In particular, Section 4.1 is dedicated to the effects during AGN bursts. Second, we select the subset of models that best corresponds to observations, which we discuss in detail in Section 5.

The first column of Table 2 shows the mass ejected as a galactic wind, illustrating that the bulk of the mass produced by stellar evolution is ejected as a galactic wind (see also Figure 2). Such galactic winds in our model are supported by thermalization of stellar motion and in particular by heating provided by Type Ia supernovae. Due to the time dependence of the supernovae and star formation, the specific heating rate increases with time.

Zoom In Zoom Out Reset image size

Earlier work by Renzini et al. (1993) and Ciotti & Ostriker (2001) has shown that Type Ia supernovae are capable of driving winds from the outer parts of elliptical galaxies but have little effect on the inner ∼1 kpc region. Within this radius, feedback from the central AGN prevents continual infall and can drive material to radii where supernova winds dominate.

The correlations we observe with dust abundance are tied to this assisting role of the black hole. AGN feedback is more effective when there is more dust due to increased radiation pressure, and thus models with high dust content, such as $A_2^{\rm MW}$, are able to drive out more mass in winds early in the simulations during periods of intense bursting. However, the dust abundance can have the opposite effect at late times since the black hole is not able to accrete dusty gas as effectively and bursting may stop. Consequently, the $A_2^{\rm MW}$ model has relatively little wind at late times whereas models with little dust continue ejecting mass throughout the duration of the simulation. Additionally, larger black holes have higher Eddington luminosities and are thus more effective in driving winds, and the black holes grow more in models with little dust. Thus, dust-rich models tend to eject mass in winds at early times more so than dust-poor models, whereas dust-poor models have significantly more winds at late times.

The next column reports the mass of new stars formed over the simulation. Note that this value is always intermediate between ΔM_BH and ΔM_w. This fact has important cosmological implications as it clearly shows how Type Ia supernovae are responsible for the metal pollution of the IGM since the bulk of the gas is ejected, not locked into new stars. As already described in CO07, AGN feedback has competing effects on star formation, acting as both positive feedback during bursts and as negative feedback at the end of each burst. This leads to the surprising result that pure cooling flow models may form fewer stars than models with AGN feedback. It is also known that, at least in 1D models, the bulk of star formation happens in a region of about a few hundred parsecs in size where cold shells are formed by recurrent cooling instabilities and shocks induced by AGN feedback. Therefore, we expect a correlation between the number of bursts and the number of new stars formed (see Table 2).

Table 3. IR Duty Cycle

Model	z: 3–2.5		z: 2.5–2		z: 2–1.5		z: 1.5–1		z: 1–0.5		z: 0.5–0
${\rm A}_2^{\rm ND}$	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
A2	0.23	0.08	1.14	0.71	0.76	0.49	0.16	0.09	0.03	0.01	0.29	0.22
${\rm A}_2^{-4}$	0.31	0.10	0.22	0.07	0.75	0.46	1.33	0.47	0.04	0.01	0.07	0.05
${\rm A}_2^{\rm CE}$	0.44	0.02	0.39	0.00	0.34	0.00	0.89	0.58	0.75	0.52	0.04	0.02
${\rm A}_2^{\rm CE2}$	0.44	0.00	0.39	0.01	0.39	0.00	2.03	1.20	11.60	1.27	0.06	0.03
${\rm A}_2^{-3}$	0.46	0.13	0.36	0.12	0.22	0.06	0.18	0.06	0.15	0.05	5.19	0.91
${\rm A}_2^{\rm MW}$	0.50	0.36	0.34	0.22	0.11	0.07	0.04	0.02	0.00	0.00	0.00	0.00

Notes. Each cell contains the percentage of time during a given redshift range that the galaxy has an IR luminosity comparable to LIRGs (10¹¹ L_☉, left) and ULIRGs (10¹² L_☉, right). In most models, phases of intense IR output have petered out below redshift ≃ 1 and in nearly all by ≃ 0.5.

Download table as:  ASCIITypeset image

In the next column, we report the total amount of gas in the simulation. Overall, the gas masses are consistent with observations of galaxies with comparable velocity dispersions (Canizares et al. 1987; Kim & Pellegrini 2012). The correlation between gas mass and dust abundance is tied to both the ability of dust to prevent accretion and the more nuanced effects of dust on galactic winds discussed above. The A₂ and $A_2^{\rm ND}$ models, which have settled into an outflow state by the end of the simulation, have the least gas.

In the next column is the final X-ray luminosity of the hot gaseous corona of our models obtained by integrating the gas emissivity in the 0.38 keV band within the volume of 10 effective optical radii. For all models, luminosities are in the observed range (Boroson et al. 2011). The low luminosities of models $A_2^{\rm ND}$ and A₂, in conjunction with their low total gas mass, demonstrate that these models are in a global wind phase at the present time. The luminosities of the other models are consistent with inflow/partial wind states. At the end of the simulation, all models are in a state of hot, low-luminosity accretion.

In the next column, we give the dust to gas ratio over the galaxy at the end of the simulation. The models are ordered as expected from the dust physics. By construction, $A_2^{\rm ND}$ is inconsistent with observations since it has no dust, and $A_2^{\rm MW}$ has too much dust relative to observed giant ellipticals. Additionally, $A_2^{-3}$ and $A_2^{\rm CE2}$ are on the high end of what would be expected (see Section 5.1). However, the enhanced dust content of the $A_2^{\rm CE2}$ model is not directly related to the dust treatment, but rather is due to a large star formation episode following the last AGN bust (see bottom panel of Figure 3). More extensive exploration of the $A_2^{\rm CE2}$ model is needed to determine if these star formation episodes are a generic feature of the model.

Zoom In Zoom Out Reset image size

The next column reports the ratio of the IR dust emission from reprocessed optical and UV radiation from stars and the black hole to the total effective luminosity of stars in the optical band. This ratio varies widely between models, and is thus an important observational diagnostic. The ratio has the expected behavior—as the dust abundance increases, the IR luminosity increases and the optical luminosity decreases due to absorption. Thus, the ratio should increase with increasing dust, which is the observed behavior. This trend is also evident in the next column, which gives the optical depth in the optical band to the center of the galaxy. For elliptical galaxies, Smith et al. (2012) find a ratio of FIR to B-band luminosity −2.5 < log L_FIR/L_B < −1.5 with a number of upper limits at the lower end of the range. We stress that our L_IR/L_Op* does not correspond exactly, but we can still make some useful comparisons. If we include upper limits, all models are in agreement. However, when restricting the comparison to detections, $A_2^{\rm ND}$ and A₂ are clearly ruled out and $A_2^{-4}$ and $A_2^{\rm CE}$ are only marginally consistent.

In summary, all models produce acceptable results from a hydrodynamic point of view. We can however exclude models based on their dust content—$A_2^{\rm ND}$ and $A_2^{\rm MW}$ clearly have too little and too much dust, respectively, and there is tension between observations of dust in elliptical galaxies the high dust content of models $A_2^{-3}$ and $A_2^{\rm CE2}$.

Now we discuss the energetic aspects of black hole accretion. The next column illustrates a factor of ≃ 10 spread in black hole growth by the end of the simulations. The mass growth is strongly correlated with the dust to gas ratio of the galaxy, with low dust models having more black hole growth. Dust grains, which have large UV absorption cross-sections, absorb UV photons and thus momentum from the luminous black hole. This radiative momentum in turn props up the gas, retarding its rate of accretion. Thus, the presence of dust tends to screen the black hole from accreting gas. Similarly, radiative feedback is able to more effectively terminate bursting events in models with more dust.

A simple check on the validity of a given model is whether the final black hole mass is consistent with the M_BH–σ relation or whether the final black hole mass is too large. Stellar evolution over a cosmological time releases an amount of gas into the galaxy equal to ≃ 30% of the initial stellar mass. If more than ≃ 1% of this gas were to be accreted, the M_BH–σ relation would be violated. The black hole growth in our models never exceeds ≃ 10⁹ M_☉ (see Figure 2, top panel), which preserves the Magorrian relation we assumed at the outset of the simulation

The next three columns of Table 2 give the integrated effective luminosity in the indicated band in units of ΔM_BHc². We recall that the adopted the electromagnetic efficiency of our simulations is ADAF-like, declining at low accretion rates and saturating to a prescribed value at high accretion rates. We use a saturation value of 0.125 (see CO07 Equation (33)). Since the values of e_Bol are roughly constant and near the saturation value, the bulk of accretion must occur at high accretion rates independent of the dust treatment. However, the distribution into different bands is sensitive to dust due to opacity effects.

e_UV indicates the amount of dust during periods of high quasar luminosity, so it is naturally maximal in the no dust model $A_2^{\rm ND}$. The optical output tends to follow the UV in its overall behavior. The hard X-ray output is very similar in all models since the dust plays no part in its transmission. However, this component is slightly lower in the $A_2^{\rm MW}$ model since the maximal dust model emits a larger fraction of its energy at low Eddington ratios where the overall radiative efficiency is lower in the A type of models.

Two of the energy output columns allow us to discriminate cleanly among the models, eliminating those having observational properties inconsistent with known data. One important ratio is that of the total quasar electromagnetic output to the observed AGN optical, as inferred by e_Bol/e_Op. This ratio, the "bolometric correction," has been classically estimated to be in the range of 5 to 10 (Soltan 1982; Yu & Tremaine 2002). Richards et al. (2006) created composite spectral energy distributions (SEDs) of 249 quasars using photometry from Spitzer and the Sloan Digital Sky Survey. They measured a mean ratio of bolometric luminosity to total optical luminosity (integrated from 0.1 to 1 μm) of 2.9 ± 1.5, with values ranging between 1.8 and 19, and a bolometric correction to the 5100 Å flux of 10.3 ± 2.1. Most of our models fall comfortably within the 5-10 range, though the $A_2^{\rm CE2}$ and $A_2^{-3}$ models have a bolometric correction exceeding 20.

Additionally, Richards et al. (2006) report an integrated IR flux between 1 and 100 μm for their quasar sample, with no corrections made for the ISM of the host galaxy. The ratio of the mean bolometric luminosity to the integrated IR luminosity is 2.58 ± 0.75, and IR to optical ratio of 1.3 ± 1.5. These ratios spanned a range of 1.1 to 5.8 and 0.36 to 18, respectively. With the exception of the $A_2^{\rm ND}$ model, all models have total IR (e_IR) exceeding total optical (e_Op) by a factor greater than two and as much as 8. $A_2^{\rm CE2}$, $A_2^{-4}$, and $A_2^{\rm MW}$ have values closer to the mean. For the radiation output from the AGN itself, we have implicitly assumed a total IR to optical ratio of 1.2. Deviations from this value are due entirely to processing by the galaxy.

In Figures 4 and 5, we present the evolution of the effective optical luminosities of the black hole and stars, the total IR luminosity (L_IR from Table 2 plus a contribution from the central black hole), the Eddington fraction, and the black hole mass. The two figures consider separately models with constant dust to gas ratios and those where this ratio varies with time and radius. The top panel of Figure 5 shows the evolution of the A₂ model taken from Ciotti et al. (2009b) with the improvements detailed in Section 3. The time evolution of each simulation has some variation from model to model, but the AGN activity of all models declines with cosmic time (see also Table 3). This decline demonstrates how the main driver of secular evolution is the relative importance of mass injection (declining as ≈t^−1.4) and supernova heating (declining as ≈t⁻¹), so that the specific heating of the galaxy declines and galaxies develop a global wind. The sharpness of the bursts is due to the use of the A family of models, which have sharper bursts and shorter duty cycles than the B family. All differences above these general trends are due to the treatment of dust.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The duty cycles shown in the last column of Table 2 are somewhat shorter than the 0.01 typical of current observations. As previously discussed, 1D simulations have inherently less steady accretion due to the inability of the gas to fragment. Additionally, the A class of models has routinely produced short duty cycles due to its fixed efficiency for driving winds resulting in short duration bursts. In contrast, the more intricate B class of models typically produced higher duty cycles (Ciotti et al. 2009b).

Figure 6 shows that in all models, most of the energy is emitted at or above the Eddington limit. However, the amount of time spent at a given fraction of Eddington varies substantially among the dust models considered here, with very low duty cycles being typical. A very small fraction of the time in all models is spent above L_Edd. To gauge how much time each model spends in a quiescent phase, we also plot the amount of time spent below a given fraction of Eddington. Ho (2009) finds that roughly 50% of AGN have L_BH/L_Edd < 10⁻⁵. As the dust content of the models goes down, the time spent at high Eddington fraction increases. This supports the idea that gas is more easily able to stream to the center of the galaxy in low dust models, resulting in sharp luminous bursts. In contrast, high dust models require more gradual buildup of cold dense shells of infalling gas before being able to overcome the radiative pressure exerted by the central black hole.

Zoom In Zoom Out Reset image size

Aird et al. (2012) find that the probability density function of finding a galaxy with a specific Eddington ratio is well-described by a power law between Eddington fractions of 10⁻⁴ and 1. We plot the corresponding cumulative distribution function for comparison in Figure 6 and find again that our 1D models spend too little time at Eddington ratios of 10⁻⁴–10⁻¹. However, the data is in rough agreement with the simulations for very high Eddington ratios.

Taken together, these aspects of the models clearly identify those that are unphysical. The $A_2^{\rm ND}$ model produces too little IR emission while the $A_2^{\rm MW}$ produces too much. The standard A₂ model at late times has very little IR output relative to optical, suggesting that it has too little dust given its stellar mass. $A_2^{-4}$, which has similar levels of depletion relative to the Galactic dust to metals ratio (D = 0.75 × 10⁻²), has an order of magnitude greater IR output as it also has roughly an order of magnitude more gas and therefore dust, bringing its value of L_IR/L_Op* closer to the observational value. However, the $A_2^{-4}$ model does not have adequate dust to produce significant obscuration, which is inconsistent with high obscured fractions. Like the $A_2^{-4}$ model, both models employing the continuity equation treatment of the dust abundance have a reasonable amount of IR emission, in addition to having sensible values for both the black hole mass and the duty cycle. In conclusion (and perhaps not surprisingly), the $A_2^{\rm CE}$ and $A_2^{\rm CE2}$ models pass the preliminary screenings better than the other models, and we will focus on these models in Section 5.

4.1. Burst Behavior

Due to the relevance of the black hole accretion physics, we now discuss the burst behavior. Each burst begins with the formation of a cooling gas shell at ≲ 1 kpc from the center of the galaxy. In 1D simulations, this cold shell starts to fall toward the center and compresses the gas interior to it. As the gas density is increased, the black hole luminosity also increases. As soon as the black hole reaches ≃ 0.01L_Edd, pre-heating instabilities appear and the accretion becomes unstable with shock waves propagating toward the falling cold shell. Fresh material is carried to the black hole by reflected shock waves. The gas in the cold shell is compressed and star formation is induced. However, the piling up of cooling material from outside the shell pushes the cold material to the center. The accretion of this material produces a large final accretion event that quenches star formation.

This general evolution is naturally affected by gas opacity, which determines how well radiation pressure works against the falling shell. Therefore, it is not surprising that the details of each burst change with the different dust treatments.

Figure 7 demonstrates the effects of changing the dust content of the gas on the burst dynamics. One burst episode was selected from each model near 3 Gyr, then scaled such that the maximum L_BH/L_Edd occurs at Δt = 0. We note that the short few Myr duration of the bursts in these models are a feature of the A family of models, and that the more complicated B models have burst durations of ≃ 10 Myr. While the common epoch for the burst ensures some level of consistency in the galaxy evolution among models, the $A_2^{\rm ND}$ model has undergone significantly more black hole growth by this time than the other models, which must be taken into account when interpreting the results. For each model, we plot a number of relevant quantities that change through the burst—the X-ray luminosity of the hot ISM in the 0.38 keV band, the optical depth to the center of the galaxy in the optical band, the dust luminosity L_IR, 〈T_d〉 as described in Equation (35), and the SFR.

Zoom In Zoom Out Reset image size

Overall, the burst evolution follows the qualitative picture of the hydrodynamics given at the beginning of this section irrespective of dust treatment. The black hole luminosity rises rapidly followed by a decline due to the expansion of gas in the central region. Coincident with the peak in black hole luminosity are peaks in both L_IR and the SFR. Following this positive feedback on star formation, the SFR drops due to the AGN feedback.

There are small but important differences in the evolution among the models. The trends are best illustrated by the models with constant dust to gas ratios, as the changes are often monotonic with this ratio. For instance, the dustier models have a faster decline in black hole luminosity after the peak. The black hole is more effective in pushing gas away in models with more dust, which slows accretion. Similarly, the drop in SFR is monotonic in dust to gas ratio since feedback is faster and more effective in dustier models.

We now ask how this picture is modified when the dust abundance is treated in a more realistic way. The simplest physical treatment is the A₂ model where the dust depletion is a simple function of temperature. The bursting behavior of this model is illustrated in the right-hand side of Figure 7. Prior to the burst, this model looks very much like the $A_2^{\rm MW}$ model since it has high values of τ_Op and L_IR. After the burst, however, the AGN heats the gas in the galaxy, which, due to the temperature-dependent dust to gas ratio, instantaneously destroys the dust. Indeed, the A₂ model closely resembles the $A_2^{\rm ND}$ model following the peak, notably in the slightly enhanced duration of the burst and its relative lack of suppressed star formation.

In the $A_2^{\rm CE}$ and $A_2^{\rm CE2}$ models, the dust must form and be destroyed on more realistic timescales. In the cold shells, the decreased temperatures allow for grain growth via collisions. From Figure 1, the grain growth time in a 10⁴ K shell of cold gas with density 10³ cm⁻³ is ≃0.01 Myr, short enough to ensure the shell is dusty. However, if the shells do not reach these densities, the growth time can become long compared to the infall time, rendering the dust unable to affect the dynamics. This is in contrast to A₂ in which cold gas would by assumption be immediately restored to MW-like grain abundances. Indeed, the right panel of Figure 7 has little evidence for enhanced dust for either the $A_2^{\rm CE}$ or $A_2^{\rm CE2}$ models, nor does the total dust mass plotted in Figure 3 show any evidence for enhancement during prior to ≃5 Gyr despite many bursts. However, the existing dust is able to affect the dynamics in ways comparable to the $A_2^{-3}$ and $A_2^{-4}$ models, notably more effective AGN feedback leading to shorter bursts and suppressed star formation. The mixing time in the $A_2^{\rm CE2}$ model does not appear to have noticeable effects on the hydrodynamics on the timescale of this burst.

By inspection of Figures 4 and 5, it is obvious that not all bursts are as sharp as that expanded in Figure 7. In general, a series of bursts culminates in a stronger final burst with considerably more time structure. These episodes are easily identified in Figures 4 and 5 as the thickest bands. We recall that in the B family of models, not discussed in this paper, the majority of bursts are of this kind. These bursts have longer duration during which a significant amount of material is accreted and are usually followed by long periods of quiescence for the galaxy. Models with little or no dust have correspondingly less radiative feedback from the AGN, and are thus characterized by short, clean bursts. In contrast, the dustier models have more long duration bursts as more material is allowed to build up and then accrete.

For illustration, in Figure 8 we expand the burst around 5 Gyr in model $A_2^{\rm CE}$ and give the Eddington fraction, effective optical luminosity from the black hole, τ_Op, L_IR, 〈T_d〉, and the SFR through the burst. Note that the time axis in the left panel spans 150 Myr, while the right panel has the same 4 Myr span as in Figure 7.

Zoom In Zoom Out Reset image size

In this model, the dust optical depth has sufficient time to rise above unity before being stopped by the dust-destroying AGN luminosity due to the buildup of a large, dense shell. Once the large structure is able to collapse, a cavity forms as the newly-fueled central AGN is able to drive out gas, dropping the accretion rate to effectively zero. In this time, gas will again accumulate until it becomes cool and dense enough to accrete. The galaxy oscillates between these modes for tens of Myr before returning to the equilibrium configuration, as illustrated in the left panel of Figure 8. Due to the prolonged existence and extreme density of the cold shell, grain growth becomes important, with spikes of grain growth evident in Figure 3 and the dust to gas ratio saturated in the cold shell evident in Figure 9. These dramatic bursts illustrate the close interplay between grain growth in cold shells and the radiation pressure that supports them.

Zoom In Zoom Out Reset image size

While the two models employing the dust continuity equation have similar overall behavior, the "two stream" approach results in more bursts of this nature. This can be attributed to the lag time between dust creation and mixing, which makes the AGN feedback less effective in the early stages of the burst before mixing can occur.

Each simulation discussed above can be assessed by its ability to reproduce the observed properties of elliptical galaxies containing supermassive black holes. Additionally, we can assess the importance of dust in the determination of each of the observational characteristics we present by analyzing the variation in these quantities among the simulations.

Of course, 1D simulations cannot adequately describe the observational signature of a galaxy viewed, e.g., along the jet axis. However, for most orientations, a 1D approach is sufficient to model the gas and dust intercepted by the line of sight. Comparisons to actual observations must be made with care bearing these limitations in mind.

5.1. Dust in Quiescent Early Type Galaxies

It is well-established that early type galaxies harbor very little dust due to rapid sputtering of grains in hot gas, with Clemens et al. (2010) putting an upper limit on grain lifetimes of 46 ± 25 Myr. However, Spitzer and Herschel have enabled study of the dust that is present and are providing important clues on the origin of that dust. Here we summarize some recent results on the dust in elliptical galaxies and compare with our simulations.

The Herschel Virgo Cluster Survey detected dust emission in 46 of 910 early type galaxies (ETGs) in their sample (di Serego Alighieri et al. 2013), with total dust masses ranging between 7  ×  10⁴ and 1.1 × 10⁷ M_☉. They further note that these masses are greater than expected for a passively evolving galaxy, and cite a potential external origin for the dust.

The Herschel Reference Survey performed a similar study on 62 ETGs, detecting dust in 31 (Smith et al. 2012). They too find that the dust masses exceed predictions for passively evolving galaxies after accounting for sputtering in hot gas, and posit that the excess dust may be the result of mergers. Both studies find a lack of correlation between the dust mass and stellar mass, casting doubt on the hypothesis that the dust originates solely from stellar outflows.

Using far-infrared Spitzer data, Temi et al. (2007a) analyzed the SEDs of 46 elliptical galaxies, finding large (∼100) variations in 70 μm and 160 μm luminosity even for ellipticals with the same B-band luminosity. Six galaxies showed extended 70 μm emission that was in excess of what would be predicted by dust production and sputtering rates. Further, none of the galaxies showed evidence of recent mergers and indeed some had quite old stellar populations. Observing dust emission as extended as 5–10 kpc, the authors suggest that the dust has been buoyantly transported out from a dusty nuclear region on a timescale less than the sputtering time.

Martini et al. (2013) use Spitzer observations of 38 ETGs to conclude that ETGs without dust lanes tend to have less than 10⁵ M_☉ of dust. Additionally, like di Serego Alighieri et al. (2013) and Smith et al. (2012), there is a large scatter in the inferred dust mass at a fixed stellar mass. Like Temi et al. (2007a), they conclude that mergers cannot alone count for the excess dust as the expected merger rate is too slow relative to the dust destruction time. They propose instead that grain growth can occur in externally accreted cold gas, with the enhanced lifetimes of the dust in the cold gas sufficient to explain the excess.

In our continuity models, which do not include any non-secular processes such as mergers or accretion of cold gas from the IGM, dust growth is able to occur in the cold gas produced by cooling flow instabilities. The presence of such gas is attested by multi-wavelength observations of giant ellipticals, revealing a cold ISM component (Werner et al. 2014). In the $A_2^{\rm CE}$ model in particular, the dust to gas ratio at z = 0 is a typical 10⁻⁴ while the dust mass is ≃ 10⁶ M_☉, values in accord with (Smith et al. 2012). Thus, our most detailed model is able to reconcile observations with theoretical estimates of dust production and destruction rates.

Additionally, the $A_2^{\rm CE}$ and $A_2^{\rm CE2}$ models predict that the distribution of dust in a quiescent galaxy (see Figure 9, left panel) that is concentrated within the inner 100 pc and then sharply declining. Although the right panel of Figure 9 is a snapshot of the radial profile in the midst of a complex burst, the AGN luminosity at that precise time (see Figure 8) is very low, and thus this object too would be interpreted as quiescent. The dust distribution of the galaxy at this time is markedly different, with high dust to gas ratios seen out to ≃ 10 kpc scale, similar to what is observed by Temi et al. (2007a, 2007b).

Finally, our models also anticipate a large variation in L_IR while L_Op* remains relatively fixed. In Figure 10 we give the histogram of the infrared dust luminosity in equally-spaced time intervals over the simulation. In both the $A_2^{\rm CE}$ and $A_2^{\rm CE2}$ models, the typical dust luminosity varies between ≃ 10⁴⁰ and 10⁴² erg s⁻¹.

Zoom In Zoom Out Reset image size

The dust distribution predictions of our $A_2^{\rm CE}$ and $A_2^{\rm CE2}$ models lends itself to a simple observational test. Because the dust abundance declines sharply with radius, we find the ratio of the half radii of the IR emission from dust and X-ray emission from the hot ISM to be ≃ 0.2 in these models. In contrast, this ratio has a value of ≃ 1 in models with constant dust to gas ratios.

5.2. Dust in Galaxies with AGN

5.3. Obscured Fraction

5.4. Star Formation Rate

As already described in CO07, all of our models predict a period of AGN-induced star formation, the so-called "positive feedback" (see also Ishibashi & Fabian 2012; Zubovas et al. 2013), with star formation occurring within the inner few hundred parsecs in the galaxy. Following this period, star formation is quenched to below pre-burst levels. Though the interplay is complex, it is evident that the black hole accretion rate (BHAR) and the SFR are closely entwined.

Chen et al. (2013) sought evidence of a BHAR–SFR relationship by studying the average BHARs of AGN as determined by their X-ray luminosity and looking for correlations with the SFR as inferred from the IR luminosity. Due to the intense variability of AGN on timescales short compared to star formation time, averaging is emphasized as painting a clearer picture of the relationship. They find that

$$\begin{equation} \log (L_X{\rm [erg s^{-1}]}) = 30.37 + 1.05\log (L_{\rm IR}/L_\odot), \end{equation} \tag{ 41 }$$

for their best-fit model, which analyzed galaxies in the redshift range 0.25 < z < 0.8 and with SFRs 0.85 < log SFR/M_☉ <2.56. Converting to BHAR and SFR, they obtain

$$\begin{equation} \log {\rm BHAR} = -3.72 + 1.05\log {\rm SFR}. \end{equation} \tag{ 42 }$$

To compare this result with our most physical simulated galaxies, in Figure 11 we make the same cuts in redshift and SFR and consider the BHAR and SFR in the simulation at equally-spaced times. Since the observational data do not have objects with L_X > 10⁴⁴ erg s⁻¹, and since these objects are likely to be obscured in our simulations at variance with the 10²⁰ cm⁻² column density assumed by Chen et al. (2013), we removed all points with L_X > 10⁴⁴ erg s⁻¹.

Zoom In Zoom Out Reset image size

Indeed, there is a strong linear correlation in all models between the BHAR and SFR. However, the points cluster more closely to the line BHAR = SFR/500, which Chen et al. (2013) derived from the M_BH–M_bulge relations of Marconi et al. (2004) than to the observations of Chen et al. (2013). Nevertheless, the slopes appear consistent. The L_X-L_IR plot varies significantly from the BHAR–SFR plot for our models. This could be partially due to rapid variations in the X-ray luminosity at relatively constant L_IR, to which the $A_2^{\rm CE2}$ model would be particularly susceptible given its delayed dust mixing. Averaging the data in 50 Myr time bins (since we cannot average over an ensemble of galaxies) brings the simulations into reasonable agreement with the observational data as shown in Figure 11, although the uncertainties are large.

It must also be noted that our models consider only AGN-induced star formation, and thus by neglecting star formation induced by other processes, e.g., mergers, we are likely under-predicting the total SFR. Secondly, we are modeling a single galaxy with a single velocity dispersion, not an ensemble of galaxies, so a quantitatively exact comparison is beyond the scope of this work. These caveats notwithstanding, AGN-induced star formation appears at least roughly consistent with the observed SFR–BHAR correlation.

By implementing a more physically-based dust treatment into 1D hydrodynamical simulations of the evolution of massive elliptical galaxies, we are able to link the computed IR emission from the galaxy during various stages of secular evolution with observations of IR emission. These models are capable of attaining LIRG and ULIRG-like phases of high IR emission without needing to invoke non-secular processes.

Despite the differing assumptions on dust abundance, the simulated galaxies illustrated a remarkably robust mass budget—in each simulation, the vast majority of the gas in the galaxy was expelled in outflows, about 10% was turned into stars, a few percent was accreted onto the central black hole, and a few percent remained as gas. The black hole growth is consistent both with current determinations of the Magorrian relation and the empirical fact that quasar "on" periods decline in frequency with decreasing redshift.

Our most physical dust models are able to reconcile the low observed dust abundance of quiescent galaxies (dust to gas ratios of ≃ 10⁻⁴) with presence of heavily obscured quasars through grain growth and reduced sputtering rates in cold gas. Additionally, optically-thick gas was able to oscillate between accretion and outflow phases for tens of Myr, resulting in sustained periods of large IR luminosity consistent with LIRGs and ULIRGs. However, at variance with Debuhr et al. (2011), τ_IR never exceeds unity and the momentum imparted to the dust gas never exceeds L_BH/c.

We identify two distinct types of AGN bursts common to all models—short-duration optically thin bursts that eventually culminate to a single large, complex burst that is largely optically thick. A clear prediction of this work is the presence of infrared emission from ≃ 100 K dust grains in the inner ≃ 1 kpc of massive galaxies during AGN bursts.

The presence of dust grains in accreting gas was also found to impact the star formation processes in the galaxy—AGN feedback and consequent quenching of star formation was enhanced in models with more dust. Similarly, dusty models also accrete less gas and have shorter duration bursts. Irrespective of our dust treatment, we find periods of "positive feedback" on star formation in which AGN activity precipitates a brief period of active star formation.

An inherent limitation of 1D simulations is the inability to account for fragmentation of gas. The influence of dust in this case, particularly in its role of preventing gas from accreting, has yet to be determined using a detailed physical prescription for the dust abundance. The formalism laid out in this work can be easily generalized to higher dimensional simulations, and given the importance of dust not only in the dynamics but also observational signatures, doing so may shed additional light on evolution of these galaxies.

We thank Bruce Draine, Jenny Greene, Jill Knapp, and Greg Novak for helpful discussions and the anonymous referee for detailed feedback that significantly improved the quality of this paper. B.H. acknowledges support from the National Science Foundation Graduate Research Fellowship under Grant No. DGE-0646086. L.C. is supported by the grant MIUR 2008, and PRIN MUIR 2010-2011, project "The Chemical and Dynamical Evolution of the Milky Way and Local Group Galaxies" 2010LY5N2T.