Infrared Emission from Cold Gas Dusty Disks in Massive Ellipticals

Dust grains are injected into the ISM by both AGB winds and supernovae. Grains are quickly coupled to gas via collisions, and so they are passively advected. Collisions with hot gas atoms can erode grains, while dust in cold gas can grow via accretion of metals. In parallel to the treatment of metals (Equation (4)), we implement these processes in a dust continuity equation (see also Hensley et al. 2014):

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial {\rho }_{\mathrm{dust}}}{\partial t}+{\rm{\nabla }}\cdot ({\rho }_{\mathrm{dust}}{\boldsymbol{v}})+{\rm{\nabla }}\cdot {{\boldsymbol{T}}}_{\mathrm{dust}}={\dot{\rho }}_{\star ,\mathrm{dust}}\\ \quad -{\dot{\rho }}_{\star ,\mathrm{dust}}^{+}+{\dot{\rho }}_{\mathrm{dust},\mathrm{gg}}-{\dot{\rho }}_{\mathrm{dust},\mathrm{gd}},\end{array}\end{eqnarray} \tag{ 6 }$$

where

$$\begin{eqnarray}&&{{\boldsymbol{T}}}_{\mathrm{dust}}=({\rho }_{\mathrm{dust}}/{\rho }_{\mathrm{gas}})\cdot {\boldsymbol{T}},\quad {\dot{\rho }}_{\star ,\mathrm{dust}}^{+}=({\rho }_{\mathrm{dust}}/{\rho }_{\mathrm{gas}})\cdot {\dot{\rho }}_{\star }^{+},\end{eqnarray} \tag{ 7 }$$

are the mass advection and sink terms due to the Toomre instability and star formation, respectively.

The mass density X_i of each element i is given by Equation (4). We partition this into a gas-phase density G_i and a density in dust grains D_i, i.e.,

$$\begin{eqnarray}&&{X}_{{\rm{i}}}={G}_{{\rm{i}}}+{D}_{{\rm{i}}}.\end{eqnarray} \tag{ 8 }$$

X_i and ${\rho }_{\mathrm{dust}}$ can be solved from the continuity Equations (4) and 6, respectively. To derive G_i and D_i, we assume for simplicity that all grains have a fixed composition. Table 1 lists the assumed mass fraction of each metal in dust $\langle {\tilde{D}}_{{\rm{i}}}/{\tilde{\rho }}_{\mathrm{dust}}\rangle$. Thus, D_i is given by

$$\begin{eqnarray}&&{D}_{{\rm{i}}}\equiv {\rho }_{\mathrm{dust}}\cdot \langle {\tilde{D}}_{{\rm{i}}}/{\tilde{\rho }}_{\mathrm{dust}}\rangle ,\end{eqnarray} \tag{ 9 }$$

with G_i following from Equation (8). This formulation allows us to track the time evolution and the spatial distribution of each chemical species in the gaseous ISM and in dust.

The source term ${\dot{\rho }}_{\star ,\mathrm{dust}}$ in Equation (6) describes the injection of dust grains into the ISM from dying stars, including both AGB winds and supernovae. Specifically,

$$\begin{eqnarray}&&{\dot{\rho }}_{\star ,\mathrm{dust}}={\dot{\rho }}_{\star }\cdot {R}_{\star ,\mathrm{dust}}+{\dot{\rho }}_{{\rm{I}}}\cdot {R}_{{\rm{I}},\mathrm{dust}}+{\dot{\rho }}_{\mathrm{II}}\cdot {R}_{\mathrm{II},\mathrm{dust}},\end{eqnarray} \tag{ 10 }$$

where ${R}_{\star ,\mathrm{dust}}$, ${R}_{{\rm{I}},\mathrm{dust}}$, and ${R}_{\mathrm{II},\mathrm{dust}}$ are the maximum ratios of dust mass to the total ejecta of AGB stars, SN Ia, and SN II, respectively, when the total available metals are transformed to the dust phase, e.g.,

$$\begin{eqnarray}&&{R}_{\star ,\mathrm{dust}}=1-{R}_{\star ,\mathrm{gas}}=\mathop{\min }\limits_{i=1...12}\left(\displaystyle \frac{{\dot{X}}_{\star ,{\rm{i}}}/{\dot{\rho }}_{\star }}{{D}_{{\rm{i}}}/{\rho }_{\mathrm{dust}}}\right).\end{eqnarray} \tag{ 11 }$$

That is, we assume that each of these sources injects as much dust of our prescribed composition—numerically, we "make" dust grains from the stellar metal yields until any of the metal ingredients runs out, then no more dust grains can be made according to the fixed dust chemical composition, i.e., each planetary nebula makes as much dust as its ejected metals permit. We calculate (${R}_{{\rm{I}},\mathrm{dust}}$, ${R}_{{\rm{I}},\mathrm{gas}}$) and (${R}_{\mathrm{II},\mathrm{dust}}$, ${R}_{\mathrm{II},\mathrm{gas}}$) in a similar way.

Finally, we track the mass exchange between the gaseous ISM and the dust by allowing grain growth ${\dot{\rho }}_{\mathrm{dust},\mathrm{gg}}$ and grain destruction ${\dot{\rho }}_{\mathrm{dust},\mathrm{gd}}$. These are described in more detail in the following sections.

In hot gas, impacting ions can erode dust grains on short timescales (Draine & Salpeter 1979). Thermal sputtering is expected to be the dominant grain destruction process in the hot ISM of an elliptical galaxy. Following Hensley et al. (2014), we adopt an analytic approximation for the sputtering rate that depends only on the gas density and temperature (Draine 2011). Letting a be the grain radius, T₆ the gas temperature in units of 10⁶ K, and ${n}_{{\rm{H}}}$ the proton number density in units of ${\mathrm{cm}}^{-3}$,

$$\begin{eqnarray}&&\ \dot{a}=-\displaystyle \frac{{10}^{-6}{n}_{{\rm{H}}}}{1+{T}_{6}^{-3}}\ \mu {\rm{m}}\ {\mathrm{yr}}^{-1}.\end{eqnarray} \tag{ 12 }$$

We assume that the grain size distribution is everywhere given by the Mathis–Rumpl–Nordsieck (Mathis et al. 1977) size distribution, i.e.,

$$\begin{eqnarray}&&\ \displaystyle \frac{{{dn}}_{d}}{{da}}=\displaystyle \frac{{{Aa}}_{\max }^{2.5}}{{a}^{3.5}}\end{eqnarray} \tag{ 13 }$$

where $({dn}/{da}){da}$ is the number density of grains with sizes between a and a+da. A is a normalization factor that can be determined from the total dust mass density ${\rho }_{\mathrm{dust}}$ via

$$\begin{eqnarray}&&{\rho }_{\mathrm{dust}}=\displaystyle \frac{8\pi }{3}A\left(1-\sqrt{{a}_{\min }/{a}_{\max }}\right){\rho }_{\mathrm{grain}}{a}_{\max }^{3},\end{eqnarray} \tag{ 14 }$$

where ${\rho }_{\mathrm{grain}}$ is the mass density of a dust grain, taken to be 3.5 g cm⁻³ as typical for a silicate grains (Draine & Lee 1984). Guided by Mathis et al. (1977) and Weingartner & Draine (2001), we adopt ${a}_{\min }=0.005\ \mu {\rm{m}}$ and ${a}_{\max }$ = 0.3 μm. The total destruction rate is obtained by integrating the mass destruction rate of grains of radius a over the size distribution, i.e.,

$$\begin{eqnarray}&&\ {\dot{\rho }}_{\mathrm{dust},\mathrm{gd}}=8\pi A{\rho }_{\mathrm{grain}}{\left(\displaystyle \frac{{a}_{\max }}{{a}_{\min }}\right)}^{0.5}{a}_{\max }^{2}\dot{a}.\end{eqnarray} \tag{ 15 }$$

This implementation of thermal sputtering implicitly assumes that the all of the gas in a cell can be described by a single temperature. However, cool, dense gas can be shock heated by supernovae, creating conditions favorable for nonthermal sputtering and grain–grain collisions at spatial scales not resolved by the simulation (Seab & Shull 1983). Grain destruction by supernova-driven shocks has been implemented in other studies by coupling the grain destruction timescale to the supernova rate (e.g., McKinnon et al. 2016; Zhukovska et al. 2016). Given the nature of the hot ISM of an elliptical galaxy and the necessarily phenomenological approach needed to incorporate this sub-grid physics, we do not consider this destruction mechanism in this work. However, as we discuss in Section 7.2, this may lead to an overestimate of the dust content in cool gas in the simulation.

When the gas density is sufficiently high, dust grains can grow from the accretion of metal atoms from the gas phase (Hensley et al. 2014, and references therein). Assuming a sticking efficiency f and average speed v_Z of these metal atoms, the grain growth rate due to collisions is $f{\rho }_{Z}{v}_{Z}4\pi {a}^{2}{n}_{d}$, where ${\rho }_{Z}$ is the mass density of the metals available for grain growth, i.e.,

$$\begin{eqnarray}&&{\rho }_{{\rm{Z}}}={\rho }_{\mathrm{gas}}\cdot \mathop{\min }\limits_{i=1...12}\left(\displaystyle \frac{{G}_{{\rm{i}}}/{\rho }_{\mathrm{gas}}}{{D}_{{\rm{i}}}/{\rho }_{\mathrm{dust}}}\right).\end{eqnarray} \tag{ 16 }$$

In this formulation, ${\rho }_{{\rm{Z}}}\to 0$ when any of the atomic constituents of dust (see Table 1) has been significantly depleted from the gas phase, thereby limiting the efficiency of grain growth.

Integrating over the grain size distribution, we obtain

$$\begin{eqnarray}&&\ {\dot{\rho }}_{\mathrm{dust},\mathrm{gg}}=f{\rho }_{Z}\displaystyle \frac{{c}_{s}}{\sqrt{{\mu }_{Z}}}8\pi {{Aa}}_{\max }^{2}{\left(\displaystyle \frac{{a}_{\max }}{{a}_{\min }}\right)}^{0.5}.\end{eqnarray} \tag{ 17 }$$

In the equation above, we have made the approximation ${v}_{Z}={c}_{s}/\sqrt{{\mu }_{Z}}$, where c_s is the sound speed of the gas, and ${\mu }_{Z}$ is the mean atomic mass of the metal atoms. For simplicity, we set constant f = 0.2 and ${\mu }_{Z}=16$.

It is known that dense ISM will be cooled down efficiently via inelastic collisions with dust grains (Ostriker & Silk 1973; Smith et al. 1996; Draine 2011). The cooling rate of gas per unit volume per unit time can be written as,

$$\begin{eqnarray}\begin{array}{rcl}{H}_{\mathrm{dust},\mathrm{collision}}^{\mathrm{gas}} & = & {\displaystyle \int }_{{a}_{\min }}^{{a}_{\max }}{da}\displaystyle \frac{{{dn}}_{d}}{{da}}\cdot 4\pi {a}^{2}\cdot \displaystyle \frac{1}{2}{\rho }_{\mathrm{gas}}{c}_{s}^{3}\\ & & \cdot \,\left[1+\displaystyle \frac{{T}_{4}}{1+{T}_{4}}{\left(\displaystyle \frac{{m}_{p}}{{m}_{e}}\right)}^{\tfrac{1}{4}}\right],\end{array}\end{eqnarray} \tag{ 18 }$$

where T₄ is the gas temperature in units of 10⁴ K. The factor in the brackets is a correction for the energy exchange considering the charge of dust grains and that electrons dominate the collisions when the gas temperature is high.

In the circumnuclear disk, both the gas density and dust density can be very high, which makes the process above very efficient—its timescale could be as short as ≲1 yr. Therefore, it is a stiff sink term in the energy equation of the gas dynamics, and it is numerically expensive to resolve the timescale explicitly. To overcome this difficulty, we assume thermal equilibrium when the timescale is shorter than the time step given by the Courant–Friedrichs–Lewy condition—we calculate the equilibrium temperature while assuming the thermal pressure is fixed. After that, we evaluate the energy loss according to the difference between the instantaneous gas temperature and the equilibrium temperature.

It is known that absorption by dust grains is the dominant sources of AGN obscuration, often exceeding the opacity due to electron scattering by orders of magnitude. The dust opacity in the UV, optical, and IR bands can be written roughly as,

$$\begin{eqnarray}\begin{array}{rcl}{\kappa }_{\mathrm{dust},\mathrm{opt}} & = & 300\times \left(\displaystyle \frac{{\rho }_{\mathrm{dust}}}{{\rho }_{\mathrm{gas}}}\right){\left(\displaystyle \frac{{\rho }_{\mathrm{dust}}}{{\rho }_{\mathrm{gas}}}\right)}_{\mathrm{MW}}^{-1}\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1},\\ {\kappa }_{\mathrm{dust},\mathrm{UV}} & = & 4{\kappa }_{\mathrm{dust},\mathrm{opt}},\quad {\kappa }_{\mathrm{dust},\mathrm{IR}}\ ={\kappa }_{\mathrm{dust},\mathrm{opt}}/150\end{array}\end{eqnarray} \tag{ 19 }$$

where we have scaled the opacities according to the dust-to-gas ratio normalized by the MW value. For simplicity, we only consider the absorption of UV photons from the central AGN by dust grains. The photon energy absorbed by dust grains per unit volume per unit time can be written as

$$\begin{eqnarray}&&{H}_{\mathrm{dust},\mathrm{abs}}^{\mathrm{AGN},\mathrm{UV}}={\rho }_{\mathrm{dust}}{\kappa }_{\mathrm{dust},\mathrm{UV}}\cdot \displaystyle \frac{{L}_{\mathrm{AGN},\mathrm{UV}}}{4\pi {r}^{2}}\exp (-{\tau }_{\mathrm{dust},\mathrm{UV}})\end{eqnarray} \tag{ 20 }$$

where ${L}_{\mathrm{AGN},\mathrm{UV}}$ is the AGN luminosity in the UV band. ${\tau }_{\mathrm{dust},\mathrm{UV}}$ is the optical depth of dust in the UV band along radial directions, i.e.,

$$\begin{eqnarray}&&{\tau }_{\mathrm{dust},\mathrm{UV}}={\int }_{{r}_{\mathrm{in}}}^{r}{\rho }_{\mathrm{dust}}{\kappa }_{\mathrm{dust},\mathrm{UV}}\cdot {{dr}}^{{\prime} }.\end{eqnarray} \tag{ 21 }$$

The calculation of the optical depth above is very computationally expensive because it could not be parallelized. For the sake of computational efficiency, we update the optical depth only every 100 time steps.

The dust grains gain momentum when they absorb the photons, which results in the radiation pressure on dust grains. Compared to electron scattering, the radiation pressure due to irradiation by AGNs could be very significant especially during bursts. Since we assume that the dust grains are coupled with gas dynamics, there will be a net acceleration to the gas as follows (Draine 2011):

$$\begin{eqnarray}&&{{\boldsymbol{f}}}_{\mathrm{rad},\mathrm{dust}}=\displaystyle \frac{{\rho }_{\mathrm{dust}}{\kappa }_{\mathrm{dust},\mathrm{UV}}}{{\rho }_{\mathrm{gas}}+{\rho }_{\mathrm{dust}}}\cdot \displaystyle \frac{{L}_{\mathrm{AGN},\mathrm{UV}}}{4\pi {r}^{2}c}\exp (-{\tau }_{\mathrm{dust},\mathrm{UV}})\cdot \hat{{\boldsymbol{r}}}.\end{eqnarray} \tag{ 22 }$$

At the inner boundary, we keep track of the dust accreted onto the black hole accretion disk. Since the temperature of the black hole accretion disk (even in the cold mode) is high enough so that no dust can survive, we return the metals in dust grains (which are accreted onto the center) back to the gas phase where they are ejected in the BAL winds. For simplicity, we assume the BAL wind is dust free, as it can be easily heated up by shocks, which will make it difficult for dust to survive.

At the outer boundary, we allow for inflows of low-metallicity material ($Z=0.15{Z}_{\odot };$ see Table 1) to mimic the CGM infall. For simplicity, we assume that the CGM is dust free.

Finally, we assume that the initial ISM is dust free as its temperature is high.

The simulation in this paper is designed for an isolated, massive elliptical galaxy, such as NGC 5129. The galaxy model is built using the procedure described in Gan et al. (2019b), and its main properties are summarized in Table 2 (see also Bentz & Manne-Nicholas 2018; S. Pellegrini et al. 2021, in preparation). The only difference is a radius-dependent description of the Satoh's k-parameter that is now parameterized as

$$\begin{eqnarray}&&k(r)={k}_{0}+({k}_{\mathrm{est}}-{k}_{0})\cdot \displaystyle \frac{r}{{r}_{0}+r},\end{eqnarray} \tag{ 23 }$$

where ${k}_{0}=0.42,{k}_{\mathrm{est}}=0.05$, ${r}_{0}=40$ kpc. These parameters produce a rather flat rotation curve, which approximates ${v}_{\mathrm{rot}}\sim 100$ km s⁻¹ when $r\lesssim {R}_{e}$, and then decreases outward in order to reproduce qualitatively the main features of the rotation profile of NGC 5129.

Table 2. Parameters of the Initial Galaxy Model

Galaxy Property	Value	Description
LB	$1.25\times {10}^{11}\ {L}_{B\odot }$	total stellar luminosity
${M}_{\star }$	$6.1\times {10}^{11}\ {M}_{\odot }$	total stellar mass
Re	11.36 kpc	projected half-light radius
Mg	$1.22\times {10}^{13}\ {M}_{\odot }$	total mass of the galaxy
${\sigma }_{\star }$	260 km s−1	1D velocity dispersion at 0.1Re
	0.37	ellipticity
${M}_{\mathrm{BH}}$	$6.1\times {10}^{8}\ {M}_{\odot }$	initial mass of central black hole

Download table as:  ASCIITypeset image

In our standard simulations, we start with a gas-free galaxy allowing infall and stellar evolution to provide all of the subsequent ISM. In our results section (Section 6.4), we will also describe the results of a perhaps more plausible initial state where the gas/star mass ratio at the beginning of the simulation is $\sim 30 \%$.

In our model, we have two physical channels for star formation, i.e., one via the Toomre instability, when the disk is gravitationally unstable (Toomre 1964; Goodman 2003; Jiang & Goodman 2011), and the other via radiative cooling when it is Jeans unstable. To ensure it is gravitationally unstable at the site of star formation, we allow star formation due to radiative cooling only when the local Jeans radius is smaller than the cell height. The Jeans radius reads,

$$\begin{eqnarray}&&{R}_{\mathrm{Jeans}}=3\times \left(\displaystyle \frac{{c}_{s}}{2\,\mathrm{km}\,{{\rm{s}}}^{-1}}\right)\sqrt{\displaystyle \frac{{10}^{3}\,{\mathrm{cm}}^{-3}}{{n}_{{\rm{H}}}}}\quad \mathrm{parsec}.\end{eqnarray} \tag{ 24 }$$

We have not made an allowance for shear or turbulence velocities within cells; these effects would reduce star formation in the innermost regions.

In the densest zone on the circumnuclear disk, the cooling radius of SN II feedback can be very small, i.e., the SN II heating will fade before it can affect its distant surroundings. Therefore, we allow SN II heating only when its cooling radius is larger than the cell height. The cooling radius of an SN II is approximately

$$\begin{eqnarray}&&{R}_{\mathrm{cool}}=5\times {E}_{51}^{1/5}{n}_{{\rm{H}}}^{1/5}\ {\left(\displaystyle \frac{{t}_{\mathrm{cool}}}{{10}^{3}\mathrm{yr}}\right)}^{2/5}\quad \mathrm{parsec},\end{eqnarray} \tag{ 25 }$$

where E₅₁ is the total energy output of an SN II explosion in units of 10⁵¹ erg. ${t}_{\mathrm{cool}}=\min ({t}_{\mathrm{low}},{t}_{\mathrm{high}})$ is the time of expansion of a supernova remnant before cooling. ${t}_{\mathrm{low}}$ and ${t}_{\mathrm{high}}$ are the estimates from the low- and high-temperature branches, respectively (see Kim & Ostriker 2015),

$$\begin{eqnarray}\begin{array}{rcl}{t}_{\mathrm{low}} & = & 4\times {10}^{4}\,{E}_{51}^{0.22}{n}_{{\rm{H}}}^{-0.55}\quad \mathrm{yr},\\ {t}_{\mathrm{high}} & = & 4\times {10}^{5}\,{E}_{51}^{0.38}{n}_{{\rm{H}}}^{-0.55}\quad \mathrm{yr}.\end{array}\end{eqnarray} \tag{ 26 }$$

In our previous model (e.g., Gan et al. 2019a), we assumed a top-heavy IMF for the star formation in the circumnuclear disk (Goodman & Tan 2004; Bartko et al. 2010; Lu et al. 2013; Palla et al. 2020), i.e., $\sim 60 \%$ of the mass in new stars is in massive stars, of which a large fraction is in binary systems. Thus, there is a significant probability for the massive stars to be kicked out of the circumnuclear disk, typically, with a runaway speed greater than 30 km s⁻¹. Finally, the runaway stars explode outside the circumnuclear disk and heat the ISM there. Therefore, the runaway stars may be an additional heating source for the ISM in the central galactic regions, which is a direct consequence of the star formation in the circumnuclear disk. Thus, following Eldridge et al. (2011; their Figure 2), we assume that 40% of the massive stars are runaway stars, and sample the escaping speed of the runaway stars according to the fitting formula below for the observed cumulative distribution function (CDF),

$$\begin{eqnarray}&&{v}_{\mathrm{runaway}}=30\mbox{--}100\mathrm{log}(1-\mathrm{CDF})\quad \quad \mathrm{km}\,{{\rm{s}}}^{-1},\end{eqnarray} \tag{ 27 }$$

i.e., no runaway stars have ${v}_{\mathrm{runaway}}\lt 30$ km s⁻¹, 50% of the stars have that ${v}_{\mathrm{runaway}}\lt 60$ km s⁻¹, 90% of the stars have ${v}_{\mathrm{runaway}}\lt 130$ km s⁻¹, etc. We set a fixed travel time for the runaway stars of 5 × 10⁶ yr, while the direction of the runaway stars is sampled randomly with equal probability over the solid angle of 4π.

Previously, we assumed that all SNe II inject a single characteristic explosion energy of 10⁵¹ erg s⁻¹. However, the explosion energy of an SN II should be dependent on the mass of its progenitor. In this paper, we adopt the mass-dependent SN II explosion energy in Morozova et al. (2018; their Figure 6), fitting it with the formula below,

$$\begin{eqnarray}&&\mathrm{log}({E}_{\mathrm{SNII}})=50.6+2.2\mathrm{log}(M/10{M}_{\odot })\quad \quad \mathrm{erg}\,{{\rm{s}}}^{-1}.\end{eqnarray} \tag{ 28 }$$

We now allow SN II kinetic feedback, i.e., we keep the vertical momentum component of the SN II ejecta, which corresponds to 1/16 of the total SN II energy, while the other 15/16 of the energy is still in form of thermal energy.

In Figure 3, we show the hydrodynamical properties of the circumnuclear disk, including its gas and dust components, at the end of the simulation. We can see that the gas component of the circumnuclear disk is dense and cold, which is ideal for star formation and dust grain growth. In the upper-mid panel, we see that in most of the volume, the gas is virialized, with a typical temperature around 1 keV. In the circumnuclear disk, the gas can be cooled far below 10⁴ K until it hits the temperature floor of 50 K, as we now allow for dust cooling of the hot gas. In the lower-left panel, we plot the dust density. In the lower-mid panel, we present the dust-to-gas mass ratio, which is up to ∼0.1. It is also much higher than the typical value in the solar neighborhood, which indicates that the dust depletion of metals can be very significant, though in the innermost part of the disk, AGN outbursts have efficiently destroyed all dust (Kawakatu et al. 2020).

Zoom In Zoom Out Reset image size

To demonstrate the effects of dust depletion of metals, in Figure 3, we plot the gas metallicity and total (gas and dust) metallicity in the upper-right and lower-right panels, respectively. We can see that the metallicities are significantly enhanced in the inner region, but there is a void in the spatial distribution of the gas metallicity in the outer disk, which is coincident with the peak of the dust-to-gas mass ratio (see the lower-mid panel in Figure 3). In the bulk of the hot X-ray emitting gas, the metallicity is below the solar value on average in the late times, but the BAL winds can have a highly super-solar metal abundance. The innermost and outermost regions are nearly dust free, the former because of AGN activity, the latter due to the low metallicity of cosmological inflow gas. We note that the circumnuclear disk is dusty, especially in its outer region where most of the metals are in dust grains. In Figure 4, we present the mass distribution of Fe in the forms of gas and dust grains. We can see that the mass density of the gas-phase Fe peaks in the inner region of the circumnuclear disk, and it is spread into the polar regions due to the AGN winds. From the right panel of Figure 4, we can see the significant dust depletion of Fe; almost all of the Fe is in dust grains in the outer disk.

Zoom In Zoom Out Reset image size

In Figure 5, we plot the surface density of the gas and dust components with blue and orange lines, respectively. The solid lines represent the results during an outburst at t = 3.87 Gyr, while the dashed lines show the results (at t = 13.64 Gyr) when the system is relatively quiescent. We can see from the figures above that the spatial distribution of dust overlaps with that of the gas component. But unlike the gas component, the dust density peaks in the outer disk rather than in the inner disk, as the gas temperature is relatively high in the inner disk due to the AGN feedback and supernova heating, which makes it difficult for the dust grains to survive.

Zoom In Zoom Out Reset image size

In Figure 6, we plot the total mass of dust grains and gaseous ISM in different ranges of temperature. It shows that the total dust mass is typically $5\times {10}^{7}{M}_{\odot }$, and the typical masses of the cold circumnuclear disk and the X-ray emitting hot ISM are ∼10⁹ and ${10}^{11}{M}_{\odot }$, respectively.

Zoom In Zoom Out Reset image size

Dust obscuration of photons is known to be very important in many astrophysical systems, especially at UV and optical wavelengths, and the dust opacity in AGN environments can dominate over electron scattering. In next subsection, we calculate the infrared emission from the dust, where the irradiation by the AGN and starlight is an important energy source for the dust emission. Therefore, it is useful to discuss below the dust obscuration of both the AGN radiation and the starlight in/through the circumnuclear disk.

Before we evaluate the dust absorption of starlight, we present in Figure 7 the cumulative star formation during the simulation. We can see that almost all of the star formation ($\sim 3.3\times {10}^{10}{M}_{\odot }$) is embedded in the dusty disk, where the gas is over-dense. The characteristic radius of the star-forming region is only ∼20 pc and would be difficult to detect except in the nearest elliptical systems. Note that: (i) in our star formation criterion (Section 5.2), we have not yet included the shear and turbulent velocity dispersion within cells. Allowing for this would reduce star formation in the innermost regions; and (ii) we do not include the dynamical and structural effects of the new stars, which may dominate the stellar mass and shift the dynamics significantly in the innermost region. One can expect that the new stars are an important source of UV photons that can be easily absorbed by dust grains, and thus can be effective in heating up the dust, and enabling the latter to emit in infrared.

Zoom In Zoom Out Reset image size

In Figure 8, we plot the dust optical depth for the UV photons emitted in the circumnuclear disk during the starburst at ${t}_{\mathrm{age}}=4.13$ Gyr. For simplicity, the optical depth is integrated along the $\theta \mbox{-}$ direction. Again, we can see that the dusty disk is optically thick to the starlight, i.e., most of the stellar radiation is expected to be absorbed within the disk, and only the "runaway stars" would be visible to distant observers. Of course, after the outbursts when the dust clears, an A-type spectrum would be seen from stars distributed as shown in Figure 7.

Zoom In Zoom Out Reset image size

In Figure 9, we plot the dust optical depth for the UV photons from the AGN (see Equation (21)) at redshift zero, where the UV optical depth is integrated along radial directions with fixed θ, since the AGN can be treated as a point source. Comparing to the lower-left panel of Figure 3, we find that most of the dusty disk is optically thick to the UV photons from the central AGN, i.e., all of the AGN emission is obscured behind the circumnuclear disk. The solid angle of the dusty disk is ∼2.36 sr, i.e., the covering factor is ∼19% of the 4π solid angle.

Zoom In Zoom Out Reset image size

In this section, we describe how we derive the infrared emission of dust from the simulation outputs. Following Hensley et al. (2014), we assume the dust emission in infrared is in equilibrium with its heating processes, i.e.,

$$\begin{eqnarray}\begin{array}{rcl}{P}_{\mathrm{dust}}^{\mathrm{lR}} & \equiv & \sigma {T}_{\mathrm{dust}}^{4}{\displaystyle \int }_{{a}_{\min }}^{{a}_{\max }}{da}\displaystyle \frac{{{dn}}_{d}}{{da}}\cdot 4\pi {a}^{2}\bar{Q}(a,{T}_{\mathrm{dust}})\\ & = & {H}_{\mathrm{dust},\mathrm{abs}}^{\mathrm{rad}}+{H}_{\mathrm{dust},\mathrm{collision}}^{\mathrm{gas}}\end{array}\end{eqnarray} \tag{ 29 }$$

where ${P}_{\mathrm{dust}}^{\mathrm{lR}}$ is the dust emission power per unit volume, ${T}_{\mathrm{dust}}$ is the dust temperature, and $\bar{Q}(a,{T}_{\mathrm{dust}})$ is the Planck-averaged emission efficiency (see Hensley et al. 2014, their Equation (33), where the Planck-averaged absorption cross sections for silicate grains are used; see also Schurer et al. 2009; Draine 2011). Therefore, provided the heating terms, the dust temperature can be calculated according to the equation above.

The heating sources are mainly contributed to by collisional heating ${H}_{\mathrm{dust},\mathrm{collision}}^{\mathrm{gas}}$ and radiative heating ${H}_{\mathrm{dust},\mathrm{abs}}^{\mathrm{rad}}$. The collisional heating is calculated according to Equation (18), and assumes equilibrium when the cooling timescale is short as in the hydrodynamical part of the simulation. Unfortunately, the assessment of the radiative heating of dust is much more complicated—there are two radiative heating sources for the dust grains, i.e., the AGN radiation and the stellar radiation,

$$\begin{eqnarray}&&{H}_{\mathrm{dust},\mathrm{abs}}^{\mathrm{rad}}={H}_{\mathrm{dust},\mathrm{abs}}^{\mathrm{AGN}}+{H}_{\mathrm{dust},\mathrm{abs}}^{\mathrm{stars}}.\end{eqnarray} \tag{ 30 }$$

The dust heating by AGN irradiation ${H}_{\mathrm{dust},\mathrm{abs}}^{\mathrm{AGN}}$ is calculated similarly to Equations (20) and (21), but with both the UV and optical photons. However, one can not assess ${H}_{\mathrm{dust},\mathrm{abs}}^{\mathrm{stars}}$ self-consistently without including sophisticated radiative transfer, and the latter is extremely computationally expensive. Fortunately, we find that the dust is usually optically thick to the local stellar radiation of the new stars, which makes it possible to allow for some simplification, i.e., we estimate ${H}_{\mathrm{dust},\mathrm{abs}}^{\mathrm{stars}}$ directly with the local stellar radiation from the new stars ${\dot{E}}_{\star ,\mathrm{rad}}^{+}$, assuming that all of the local stellar radiation is absorbed by dust.

$$\begin{eqnarray}&&{H}_{\mathrm{dust},\mathrm{abs}}^{\mathrm{stars}}\sim {\dot{E}}_{\star ,\mathrm{rad}}^{+}\sim {\dot{E}}_{\star ,\mathrm{UV}}^{+}+{\dot{E}}_{\star ,\mathrm{opt}}^{+},\end{eqnarray} \tag{ 31 }$$

where the stellar radiation of the new stars is mainly in the UV (${\dot{E}}_{\star ,\mathrm{UV}}^{+}$) and optical (${\dot{E}}_{\star ,\mathrm{opt}}^{+}$) bands.

Following Ciotti & Ostriker (2012), we compute the UV and optical luminosity per unit volume of the new stars according to the star formation history,

$$\begin{eqnarray}&&{\dot{E}}_{\star ,\mathrm{UV}}^{+}=\displaystyle \frac{{\epsilon }_{\mathrm{UV}}{c}^{2}}{{t}_{\mathrm{UV}}}{\int }_{0}^{t}{\dot{\rho }}_{\star }^{+}({t}^{{\prime} })\cdot {e}^{-\tfrac{t-{t}^{{\prime} }}{{t}_{\mathrm{UV}}}}\ {{dt}}^{{\prime} },\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}&&{\dot{E}}_{\star ,\mathrm{opt}}^{+}=\displaystyle \frac{{\epsilon }_{\mathrm{opt}}{c}^{2}}{{t}_{\mathrm{opt}}}{\int }_{0}^{t}{\dot{\rho }}_{\star }^{+}({t}^{{\prime} })\cdot {e}^{-\tfrac{t-{t}^{{\prime} }}{{t}_{\mathrm{opt}}}}\ {{dt}}^{{\prime} },\end{eqnarray} \tag{ 33 }$$

where ${\dot{\rho }}_{\star }^{+}$ is the star formation rate density (which is given by the MACER simulations). In the equations above, ${\epsilon }_{\mathrm{opt}}=1.24\,\times {10}^{-3}$, ${\epsilon }_{\mathrm{UV}}=8.65\times {10}^{-5}$, ${t}_{\mathrm{opt}}=1.54\times {10}^{8}$ yr, ${t}_{\mathrm{UV}}=2.56\,\times {10}^{6}$ yr are the efficiency and characteristic time of optical and UV emission, respectively.

Note that we ignore the dust opacity of the infrared emission itself in the calculations above, which could be important in the near or mid-near-infrared band (Hönig 2019), and we leave it for future work.

In Figure 10, we present the results from the data post-processing of the dust infrared emission. In the calculation, we assume the gas and dust are transparent to the infrared photons; therefore, the total infrared luminosity ${L}_{\mathrm{TIR}}$ of the galaxy is simply the sum of the dust infrared emission in each cell, making it possible to assess the contributions of dust emission by the AGN irradiation, starlight, and dust–gas collisional heating, separately (in the upper panels, (a) and (b), of Figure 10). In the lower panels, (c) and (d), we present the total infrared luminosity (black line).

Zoom In Zoom Out Reset image size

In Figure 11, we plot the spatial distribution of the dust emission during the AGN outburst at t = 2.44 Gyr. It shows that the dust emission during the outbursts is mostly due to the AGN irradiation of the innermost region and the surface of the disk where the optical depth to the central AGN is not obscured. As other evidence, in the bottom panel of Figure 12, we plot the scatter of the half-light radius of the dust infrared emission versus the total infrared luminosity. We can see that the half-light radius decreases moderately when the infrared luminosity increases, with typical values of approximately 30 pc at 2 × 10⁴⁴ erg s⁻¹. In the top panel of Figure 12, we plot the scatter of the infrared emission weighted dust temperature versus the total infrared luminosity, where a tight correlation is found. While most of the gas mass is in the X-ray emitting, extended hot component, most of the infrared radiation is emitted by the highly concentrated cold component.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 13, we present the histogram of infrared luminosity in dust temperature bins, in which the vertical axis is the total instantaneous infrared luminosity integrated over the galaxy within a given dust temperature bin. The colored lines are the results from three representative time snapshots, i.e., t = 2.44 Gyr during an AGN burst without starburst; t = 4.13 Gyr without AGN burst while with starburst; and t = 9.39 Gyr when the galaxy is quiescent, respectively (see also the samples in Figure 12 for more details). We can see that the dust temperature peaks around $\sim {10}^{2}$ K, and it is systematically higher when the infrared luminosity increases, especially during the AGN outbursts when the dust heating is dominated by the AGN irradiation. The dust temperature can reach ∼500 K in the innermost region near the central AGN, while the typical dust temperature in the circumnuclear disk is ∼30 K.

Zoom In Zoom Out Reset image size

In Figure 14, we present the duty cycle of the infrared luminosity from dust, i.e., the fraction of time when the infrared luminosity is above the given values. It shows that most of time the infrared luminosity is $\gt 3\times {10}^{43}$ erg s⁻¹, and it systematically decreases when the redshift decreases. The median value is ∼2 × 10⁴⁴ erg s⁻¹ (the blue line), or $\sim 1/2$ of that at late times (the red line).

Zoom In Zoom Out Reset image size

In Figure 15, we plot the histograms of the dust infrared emission. It shows that there is a clear redshift dependency of the dust infrared luminosity, i.e., the infrared luminosity decreases systematically as the redshift decreases. Again, this is due to the AGN irradiation (which contributes the most energy for dust emission), as there is more AGN activity at high redshift (see also Figures 1 and 2). In all of the histograms, we can also see a cutoff around the infrared luminosity ∼3 × 10⁴³ erg s⁻¹. We find in the simulation that the collisions between dust grains and the gaseous medium always contribute significant infrared emission, no matter whether there are AGN bursts and/or star formation, though the latter may boost the dust collisional heating via feedback processes. As a consequence, there is a "floor" luminosity due to the collisional heating (∼10⁴³ erg s⁻¹ in our simulated galaxy; see the green line in the top panel of Figure 10), which results in the cutoff in the histogram. The stellar irradiation by the embedded new stars also contributes a significant fraction of the energy for the dust infrared emission (see the orange line in the top panel of Figure 10).

Zoom In Zoom Out Reset image size

In this subsection, we study the dependence of our results on the initial conditions and spatial resolution.

In the fiducial model, we have adopted the previously standard assumption that, after formation, elliptical galaxies have a low gas fraction, and so for simplicity, we took a gas-free initial state with all subsequent gas content due to cosmological infall or mass loss from evolving stars. This fiducial model is the "high-res" model in Table 3, and it has an inner boundary of 2.5 pc. To determine the importance of our high spatial resolution, we include, as the "low-res" model, a simulation that is identical except for an inner boundary of 25 pc, which is still small compared to the resolution of most cosmological simulations.

Table 3. Simulation Outputs of the Model Galaxy ^a

Model	AGN Duty Cycle b			${M}_{\mathrm{BH}}$ c	$\mathrm{SF}$ d	$\left\langle \mathrm{SFR}\right\rangle$ e	$\left\langle {L}_{{\rm{X}}}\right\rangle$ f , l	$\left\langle {R}_{0.5}^{X}\right\rangle$ g	$\left\langle {L}_{\mathrm{TIR}}\right\rangle$ h	$\left\langle {R}_{0.5}^{{IR}}\right\rangle$ i	$\left\langle {M}_{\mathrm{gas}}^{\mathrm{cold}}\right\rangle$ j	${R}_{\mathrm{disk}}$ k
	${f}_{l\gt 0.1}$	${l}_{\mathrm{median}}^{\mathrm{energy}}$	${l}_{\mathrm{median}}^{\mathrm{time}}$	(${10}^{9}{M}_{\odot }$)	(${10}^{9}{M}_{\odot }$)	(${M}_{\odot }$ yr−1)	(1042erg s−1)	(kpc)	(1043erg s−1)	(kpc)	(${10}^{9}{M}_{\odot }$)	(kpc)
high-res	7.62%	0.011	0.00094	6.72	32.9	1.24	0.07	14.5	8.23	0.03	0.97	1.53
low-res	9.59%	0.025	0.00042	8.32	9.84	0.15	0.10	31.6	1.60	0.41	2.05	1.52
low-res+gas	37.0%	0.039	0.00026	25.6	26.4	0.33	0.33	19.5	5.39	0.40	3.06	1.67

Notes.

^a From left to right, the columns in the table above present the simulation results as follows: ^b AGN duty cycle, including the percentage of the cumulative radiation energy that the AGN emits when its Eddington ratio $l\equiv {L}_{\mathrm{AGN}}/{L}_{\mathrm{Edd}}$ is above 0.1 (i.e., ${f}_{l\gt 0.1}$), the median Eddington ratio above which the AGN emits half of its total energy (i.e., ${l}_{\mathrm{median}}^{\mathrm{energy}}$), and the median Eddington ratio above which the AGN spends half of the simulated cosmological time span (i.e., ${l}_{\mathrm{median}}^{\mathrm{time}}$). ^c the final black hole mass ${M}_{\mathrm{BH}}$ (the initial black hole is $0.61\times {10}^{9}{M}_{\odot }$ in all simulations). ^d total star formation. ^e averaged star formation rate ($\left\langle \mathrm{SFR}\right\rangle$). ^f averaged X-ray luminosity from the hot ISM ($\left\langle {L}_{{\rm{X}}}\right\rangle$). ^g averaged half-light radius of the X-ray emission ($\left\langle {R}_{0.5,{\rm{X}}}\right\rangle$). ^h averaged infrared luminosity from dust ($\left\langle {L}_{\mathrm{TIR}}\right\rangle$). ⁱ averaged half-light radius of the infrared emission ($\left\langle {R}_{0.5,\mathrm{IR}}\right\rangle$). ^j averaged mass of the cold gas with temperature lower than 100 K ($\left\langle {M}_{\mathrm{gas}}^{\mathrm{cold}}\right\rangle$). ^k size of the circumnuclear disk at the end of the simulations (${R}_{\mathrm{disk}}$). ^l All of the averages are made using the simulation data in the last 1 Gyr.

Download table as:  ASCIITypeset image

We see some dramatic changes in the lower-resolution model. The infrared luminosity is much higher, and the dusty disk is much smaller in the high-resolution model. Additionally, the new stars are much more numerous and concentrated when resolution is increased. The black hole final mass is less of a consequence, since the gas flowing in toward the center can be turned into stars before reaching the central black hole.

The implication is strong that a still higher spatial resolution (currently beyond technical reach) would further enhance these trends. It is possible however that heating from the central black hole would limit star formation within, as we see in Figure 3, $r\lt 10$ pc.

The other variation on our fiducial model is based on starting the simulation with a gas/star ratio of 30%, close to the value seen in the recent ALMA/SCUAB-2 survey of submillimeter galaxies (Dudzevičiūtė et al. 2020). That work (see their Figure 11(d)) indicates a gas/star ratio of roughly 20%–40% at redshift z = 2. To save computer time, we ran the simulation with 30% initial gas fraction, which was identical to the low-resolution model "low-res," labeling the high initial gas model "low-res+gas" in Table 3. We see several important changes. First, it is much more "bursty" with 37% of the AGN radiation energy emitted above $L/{L}_{\mathrm{Edd}}=0.1$. As expected, there is more star formation, more infrared emission, and a higher X-ray luminosity.

Figure 15 shows the distribution of infrared luminosities for the three simulations as a function of redshift range. Both higher resolution and more initial gas push up the infrared luminosities to 10⁴⁴ erg s⁻¹ at late times and 10⁴⁵ erg s⁻¹ at redshift above z = 2.

A model with both higher resolution and more initial gas is being studied and will be reported in a subsequent communication.

Extending to high redshift, Schreiber et al. (2018) fit models to stacked SEDs from deep CANDELS fields, including stacked Herschel and ALMA data. Although their model templates applied to the observations do not include emission from a dusty AGN torus, this should have little impact on their derived ${L}_{\mathrm{TIR}}$ and stellar masses. For their highest stellar mass bin ($11\lt \mathrm{log}{M}_{\star }\lt 11.5$), both data and models suggest stacked-average ${L}_{\mathrm{TIR}}$ ∼ 3, 7, 17, 29, and 51 $\times \ {10}^{44}$ erg s⁻¹ at z ∼ 0.5, 1, 1.5, 2.2, and 3, respectively. These ${L}_{\mathrm{TIR}}$ values are in excellent agreement with our predictions (see Figure 15).

Such high-luminosity massive galaxies are likely to evolve into the quiescent massive galaxies observed in the local universe. For example, using ALMA, Wang et al. (2019) recently discovered a population of very massive, optically faint, IR-luminous galaxies at z ∼ 4 with average ${L}_{\mathrm{TIR}}$ = $40\times {10}^{44}$ erg s⁻¹. This population has a clustering/bias factor consistent with the population expected to be the progenitors of massive galaxies in groups and clusters in the local universe, such as the one modeled here.

Results from Schreiber et al. (2018) suggest that at z > 1, massive galaxies dominate the luminous end of the total infrared (TIR) luminosity function (LF). Our model suggests high IR variability with a broad, flattened distribution over log ${L}_{\mathrm{TIR}}$. Measurements of the TIR (or far-IR) LF (e.g., Gruppioni et al. 2013; Koprowski et al. 2017) are relatively consistent with the range and shape of the distribution we infer. However these LFs are themselves discrepant on the bright end, possibly due to differing selection effects with respect to AGNs or galaxies with a warmer dust component (see, e.g., Gruppioni & Pozzi 2019).

At low redshifts, for $z\lt 0.5$, our model predicts a median infrared luminosity $\sim 1.2\times {10}^{44}$ erg s⁻¹, and 90% of the time, L ${}_{\mathrm{TIR}}\lt 4\times {10}^{44}$ erg s⁻¹ (see Figure 14). While luminosities fall just below the lower bound for the typical definition of luminous infrared galaxies, these luminosities are nevertheless comparable or even higher than those of less-massive star-forming galaxies. There is evidence that some relatively quiescent massive galaxies have TIR luminosities in this range.

Andreani et al. (2018) measured the bivariate luminosity and mass functions from a local sample from the Herschel Reference Survey. For galaxies with high stellar masses (e.g., log M_⋆ > 11.5), their bivariate function implies a broad range of TIR luminosities centered on (1–2) $\times {10}^{43}$ erg s⁻¹. Their analysis only includes galaxies classified as late-types, but the massive galaxies included in their analysis are likely to contain a significant spheroid component.

Surveys of IR emission from dust in early-types in the nearby universe typically contain only a handful of galaxies with high stellar masses ($\mathrm{log}{M}_{\star }\gt 11$). Amblard et al. (2014) extract physical properties from 221 early-type galaxies of which ∼50 have high stellar masses. The median ${L}_{\mathrm{dust}}$ of this subsample is 10⁴² erg s⁻¹, with five (10%) above 10⁴³ erg s⁻¹. Infrared AKARI observations of 260 early-type galaxies in the Atlas3D sample (Kokusho et al. 2017, 2019) do include several tens of high-mass galaxies, with a spread of ${L}_{\mathrm{TIR}}$ in the 10⁴²–10⁴³ erg s⁻¹ range, and occasionally higher. Interestingly, they note much wider scatter in warm dust luminosities (at a fixed stellar mass), and they conclude that this is likely due to variation in AGN activity across the sample. Galaxies with dusty nuclei and/or post-starburst signatures also show high IR luminosities. For example, in a sample of E+A galaxies, Smercina et al. (2018) measure ${L}_{\mathrm{TIR}}$ $\sim {10}^{43.5}$ and 10⁴⁴ erg s⁻¹ for two galaxies with high stellar mass.