MODELING MID-INFRARED DIAGNOSTICS OF OBSCURED QUASARS AND STARBURSTS

Gadget-2 is a TreeSPH (Hernquist & Katz 1989) code that conserves both energy and entropy (Springel & Hernquist 2002). The simulations include radiative heating and cooling as in Katz et al. (1996), which assumes ionization equilibrium with an ultraviolet background (Faucher-Giguère et al. 2009) and no metal line cooling. SF is modeled to reproduce the global Kennicutt–Schmidt law (Kennicutt 1998) via the volumetric relation $\rho _{\rm SFR} \propto \rho _{\rm gas}^{1.5}$ with a minimum density threshold n ∼ 0.1 cm⁻³. The SF law employed should be considered an empirically and physically motivated prescription to summarize physics we do not resolve. We do not track the formation of molecular gas or resolve individual molecular clouds, and this limitation is a key uncertainty in our RT step (see Section 2.3).

The structure of the ISM is modeled via a two-phase sub-resolution model in which cold, star-forming clouds are embedded in a diffuse, hot medium (Springel & Hernquist 2003) pressurized by supernova feedback that heats the diffuse ISM and evaporates the cold clouds (Cox et al. 2006b). All supernova energy is assumed to thermally heat the hot ambient medium, and no kinetic stellar wind kicks are used. Metal enrichment is calculated by assuming each gas particle behaves as a closed box with a yield y = 0.02. SMBH particles accrete via Eddington-limited Bondi-Hoyle accretion and deposit 5% of their luminosity to the nearby ISM as thermal energy (Springel et al. 2005; Di Matteo et al. 2005). The luminosity is computed from the accretion rate assuming 10% radiative efficiency, $L_{\rm bol}= 0.1\, \dot{M}c^2$.

In this paper, we focus on two merger simulations, summarized in Table 1, that represent starbursts with different characteristics. Our "highly obscured" simulation is a very massive, gas-rich merger meant to mimic very luminous starbursts at z ∼ 2–4, while the "less obscured" simulation is more representative of typical gas-rich mergers at z ∼ 0. In each case, we use two identical progenitor galaxies. The progenitors of the highly obscured merger are composed of an exponential disk with baryonic mass 4 × 10¹¹ M_☉ (60% of which is gas) and a dark matter halo of mass 9 × 10¹² M_☉ described by a Hernquist (1990) profile. The galaxy properties are scaled to z ∼ 3 following Robertson et al. (2006). The orbital parameters are identical to those of the "e" orbit of Cox et al. (2006a). In contrast, the progenitors of the weakly obscured example have lower mass (M_* ∼ 10¹¹ M_☉), are less gas-rich (40%) than the highly obscured example, and are scaled to z ∼ 0. Each progenitor galaxy is seeded with one SMBH particle.

In this work, we focus on signatures of AGN activity that normalize out the galaxy's mass or total luminosity, and therefore the salient differences between these two cases are the relative extents to which the merger triggers SF, AGN activity, and obscuration. A disadvantage with this approach is that we are restricted in the range of ULIRG scenarios we can probe, owing in part to limited computational resources. In this paper, we have chosen to explore two merger simulations that potentially reflect the extremes of IR-luminous galaxies, and we discuss further consequences of this limitation in Section 2.4.

We use Sunrise⁹ Version 3 in post-processing to calculate the SED observed from seven cameras distributed isotropically in solid angle every 10 Myr. For a full description of Sunrise see Jonsson (2006) and Jonsson et al. (2010), and for a summary of components essential for this work and a description of other specific choices, see Hayward et al. (2011). We describe specific choices for some of the RT parameters in Section 2.3.

2.2.1. Sources

2.2.2. ISM Structure, Metals, and Dust

2.2.3. Radiative Transfer

In order to gain additional physical insight, we consider different sets of assumptions for the RT stage, including two sub-resolution models for the dust distribution on scales below that resolved by the Gadget-2 simulation. For our highly obscured simulation, we will focus on the ISM treatment referred to in Hayward et al. (2011) as "multi-phase off." This choice assumes that the dust mass contained in both phases of the Springel & Hernquist (2003) ISM model is distributed uniformly across each resolution element, which is likely appropriate for gas-rich mergers in which the central regions are composed of dense, almost exclusively molecular gas. Throughout this work, we refer to this model as our "default ISM" treatment.

Alternatively, we can assume that gas in the cold, dense clouds of the Springel & Hernquist (2003) model has a negligible volume filling factor. This choice retains the dust in the "hot" phase, which is assumed to have a volume filling factor equal to unity. Specifically, the dust mass that corresponds to this diffuse gas phase of the Springel & Hernquist (2003) model is distributed uniformly across each resolution element, while the dust mass occupying the dense, cold clouds is ignored. In this case, we neglect attenuation and emission from both the cold clouds in which stars are formed and from other clouds photons encounter along the line-of-sight. This assumption thus gives a lower limit on the amount of attenuation. We refer to this assumption as our "alternate ISM" treatment, and apply it to both our highly and weakly obscured simulations.

The fraction of gas in the cold phase, and therefore the amount of dust ignored by the alternate ISM model, varies with time and position during our simulations. For the highly obscured merger, the amount of mass in the cold phase is ∼38% (∼27%) at the peak of the starburst (AGN). For the less obscured merger, this fraction is ∼28% at the peak starburst and ∼15% at peak AGN power. Gas which is at higher densities contains a larger fraction in the cold phase, and so during merger coalescence gas in the central few kpc achieves a cold phase fraction of ∼90%. Thus the alternate ISM assumption, discarding this fraction of dust in the cold clouds, more strongly affects the central regions than the galaxy as a whole. Our RT calculation implicitly assumes that this dust does not absorb any energy, and therefore it is not included in the IR emission calculation. This dust emission would primarily affect far-IR wavelengths, and so it should not affect the analysis presented here.

The two treatments of subresolution dust structure we employ should be considered two plausible and physically-motivated, yet uncertain, ISM models that encapsulate unresolved processes. With current simulations, it is only possible to parameterize our ignorance in this way. However, any model used to interpret ULIRG SEDs is limited by this uncertainty. These simulations are thus particularly useful for their ability to quantify the uncertainty caused by dust clumpiness.

For both mergers, we perform the RT calculations after multiplying the SMBH luminosity at all times by zero (AGN × 0), one (AGN × 1), and ten (AGN × 10), allowing us to manually adjust the AGN contribution to the SEDs. Note that the effect of AGN feedback, the thermal heating of the ISM surrounding the SMBH particles, is kept fixed at the fiducial level of Section 2.1. The AGN contribution to the mid-IR SED arises from separate self-consistent RT calculations with each of our assigned AGN luminosities, dust grain models, and ISM assumptions. This enables us to cleanly test how a given indicator depends on AGN luminosity for fixed geometry and galaxy ISM conditions. At each of these AGN strengths, we apply both ISM assumptions to the highly obscured merger and the "alternate ISM" treatment for the less obscured case. For these three sets of three simulations, we perform the RT using both the MW and SMC bar dust models.

Although the "high obscuration" and "low obscuration" simulations are meant to mimic gas rich starbursts in massive major mergers at z ∼ 3 and z ∼ 0, respectively (Table 1), these cases do not necessarily accurately reflect the real ULIRG population at a given redshift. The "high obscuration" example is, by choice, a particularly massive and gas-rich merger, and thus becomes a luminous hyper-LIRG that reaches a peak log L_IR ∼ 13 briefly, having log L_IR > 12 for ∼1 Gyr. This merger is consistent with extreme starbursts such as sub-millimeter galaxies (SMGs) and hot-dust ULIRGs at z ∼ 2–4 (Hickox et al. 2012; Hayward et al. 2012, 2013), and therefore reflects a typically observed source at these redshifts (note: this does not mean it is a typical galaxy). By contrast, the "weakly obscured" example is marginally a ULIRG, briefly reaching L_IR  ∼  10¹² L_☉ at merger coalescence under our default assumptions, consistent with the idea that local ULIRGs are almost exclusively ongoing mergers (Sanders & Mirabel 1996). Note that the initial gas fractions assigned in our progenitors, 60% and 40%, are high compared to comparable star-forming galaxies at these two epochs—this is to account for our neglecting cosmological accretion and gas recycling, and a fairer comparison is to the gas fractions at the time of merger, which are ∼30% and ∼20%.

Of course, real IR-luminous galaxy samples at a given redshift will draw from a wide range of scenarios. Therefore, we cannot hope to fully represent the ULIRG population in the present study. However, we will focus primarily on the goal of separating AGNs from SF activity in a small set of experiments defined by our two mergers with the RT variants described in Section 2.3, whose properties may more broadly span plausible situations. As an example, at z ∼ 3, there are some mergers with the same stellar mass as our highly obscured simulation but a lower gas fraction. In the Eddington-limited phase, such mergers may have higher AGN-to-SF fractions than our fiducial highly obscured calculation, so the AGN × 10 model may more accurately represent such cases. Similarly, our simulations at z ∼ 0 may not be a close match to the conditions in all local ULIRGs. Moreover, while the final mass of our accreting SMBHs depends primarily on the fraction of feedback energy coupled to the ISM (e.g., Springel et al. 2005), the growth history of a SMBH, and hence its luminosity at a given time, can depend sensitively on the adopted sub-resolution model.

Given these uncertainties, the goal of our experiments—in particular the AGN × 0, AGN × 1, and AGN × 10 models—is to cover our bases and controllably boost the luminosity so that we span a wide enough range in the ratio of AGN-to-SF activity, and not to comprehensively predict the behavior of observed ULIRGs. In the future, a more realistic treatment of the population may be possible as large suites of high-resolution simulations become more common. However, many uncertainties remain not only in the physical models of the ISM and AGN accretion on sub-resolution scales, but also in how to interface emission and dust attenuation from the central engine with the host galaxy in the RT stage (e.g., Jonsson et al. 2010). Therefore, we begin with an exploration of these few numerical experiments in a first attempt to explore the applicability of this modeling technique.

In Figure 1, we summarize the predicted SEDs of our mergers. We focus on four times during each merger spanning ∼500 Myr, and present the total mid-IR SEDs emerging in a particular direction. In the left SED columns, we show the results of our fiducial ISM assumption, but artificially vary the total luminosity of the two (and eventually one) SMBH particles. In the right SED columns, we explore our assumptions about the ISM and dust, a study that we analyze in more detail below in Section 4.2 (see also Section 2.3).

Zoom In Zoom Out Reset image size

Alongside the SEDs, we show high-resolution false-color composite images corresponding to the same viewing direction as the highly obscured SEDs. The left images present the system in the U, V, and J bands (Johnson & Morgan 1953), and the second in the 3.6, 4.5, and 8 μm channels of the Spitzer Infrared Array Camera (IRAC; Fazio et al. 2004). These images provide some insight into how these systems might be classified in terms of a merger-stage diagnostic—if this system is observed at high redshift, then the only time it is obviously a merger is at t ≲ t₁. This precedes the SMBH's peak luminosity and therefore the times at which it makes an obvious contribution to the mid-IR SED, a period spanning ∼1 Gyr.

We compute X-ray fluxes from SMBH particles and subsequent attenuation by the ISM following the approach described by Hopkins et al. (2005, 2006). This method uses the intrinsic quasar SED (Section 2.2 and Hopkins et al. 2007c), and calculates its extinction at X-ray wavelengths by applying the photoelectric absorption cross sections of Morrison & McCammon (1983), as well as Compton scattering cross sections, each scaled by metallicity.

A key difference between Hopkins et al. (2005) and the present work is that their assumptions correspond to the "alternate ISM" treatment we described in Section 2.3, where the cold clumps in the ISM, and hence often a significant fraction of the dust mass, are assumed to have a small volume filling factor. This may be true under many conditions, but may not be applicable to extremely gas-rich ULIRGs with abundant supplies of molecular gas. Therefore, we choose to use the same assumptions for X-ray attenuation that we use for the IR calculations described in Section 2.3. For our highly obscured merger simulation we calculate the column densities that attenuate the X-ray flux both ways: by discarding the cold phase mass ("alternate ISM"), and by keeping the cold phase mass ("default ISM"). For our less obscured merger simulation, we use only the "alternate ISM" model, discarding the cold phase mass.

For this paper, we focus on the AGN fraction, which we define as the ratio of SMBH luminosity to total luminosity across the wavelength range used by Sunrise: 0.09–10³ μm. We denote this quantity L_agn/L_bol. Our definition uses the intrinsic AGN power regardless of how it appears in the final observed SED (i.e., it is insensitive to the amount of attenuation), but it does not include AGN emission at X-ray wavelengths. This quantity, L_agn/L_bol, is not normally available for a given observed sample. Here, the SMBH accretion rate and L_agn/L_bol are available directly from the Gadget and Sunrise calculations. Thus we can directly evaluate the effectiveness of observed indicators, defined in Section 3, at estimating the AGN contribution utilizing self-consistent calculations of the galaxy's stars, dust, and SMBH.

We compute mid-infrared diagnostics based on the rest-frame spatially integrated Sunrise SEDs, which we store as specific luminosities, i.e., the energy through the camera per unit time and per unit frequency (or wavelength) interval: L_ν or L_λ. For simplicity and for analogy with observable quantities, we denote these as f_λ ≡ L_ν(λ), where λ is in microns. Owing to computational constraints, these SEDs have very low spectral resolution (λ/Δλ ∼ 10), so we limit consideration to approximations of several observationally-motivated estimators. We focus on two spectral dust features, r_9.7, r_7.7, defined by

$$r_{9.7}= \frac{f_{9.7}}{c_{9.7}}$$

$$r_{7.7}= \frac{f_{7.7}}{c_{7.7}},$$

where c_λ is the continuum level estimated by interpolating linearly between f_5.7 and f₁₅. The remaining indicators are ratios of f_λ at 1.6, 5.7, 15, and 30 μm. These points were chosen to match common mid-IR colors, in which a point is often chosen just below 6 μm to avoid the 6.2 μm PAH and 6 μm water-ice features. Unless otherwise specified, we present all SEDs and their derived quantities in the source's rest frame.

In Figure 2, we present an example rest-frame SED in the mid-IR and highlight the spectral features that we will consider here. In addition to the final SED, we plot the intrinsic SMBH SED that we assume, which exhibits a generally flat shape with an "IR bump." If this source were observed unobscured, or obscured only by an optically thin (at several  μm) dust column, then this would resemble a "power-law-like" AGN source with red near- and mid-IR colors. We see this situation in the t₂ and t₃ panels for the less obscured merger (right side) of Figure 1: the f_4.5/f_1.6 slopes depend strongly on the AGN strength, and a comparison of the "no dust" and "alternate ISM" cases indicates that dust attenuation is insignificant at wavelengths longer than a few microns. Furthermore, a near-IR slope (such as f_4.5/f_1.6) has been used to select AGN in wide-field IR surveys (see, e.g., Stern et al. 2005, 2012). Thus we will keep this feature in mind while analyzing the mid-IR properties of more highly obscured sources.

Zoom In Zoom Out Reset image size

By contrast, the starburst SED falls as lambda increases at 1 μm < λ < 10 μm, but also exhibits distinctive features from dust grains identified as PAHs superimposed at 3.3, 6.2, 7.7, 8.6, 11.3, 12.6, and 17 μm. For example, SEDs characteristic of starbursts are seen in the first row of SEDs in Figure 1 (i.e., the black curves at time t₁).

Observationally, weaker PAH and stronger continuum (i.e., lower PAH equivalent width—EW) effectively indicate a higher relative contribution from the AGN in dusty galaxies, especially ULIRGs (e.g., Laurent et al. 2000; Sturm et al. 2000; Tran et al. 2001). This diagnostic is supported by mid-IR fine structure line diagnostics implying an AGN-like radiation field in sources of lower PAH EW (Genzel et al. 1998; Armus et al. 2007).

More recently, such mid-IR diagnostic studies discovered sources that have a small PAH EW but a deep 9.7 μm silicate absorption feature (e.g., Houck et al. 2005), which seems to require a deeply embedded, centrally concentrated source (Levenson et al. 2006). These sources are believed to represent Compton-thick AGNs (Bauer et al. 2010; Georgantopoulos et al. 2011), although there is at least one known star forming dwarf without an AGN whose mid-IR spectrum fits the above criteria (Roussel et al. 2003).

In Section 5, we will compare our predictions with data for local and high-z starbursts with available low-resolution Spitzer IRS spectra. Specifically, we use the ULIRG sample of V09 and a sample of 192 24 μm-selected z ∼ 0.3–3 starbursts and obscured quasars. The redshifts and IRS spectra of the 24 μm-bright sample are presented in Yan et al. (2007) and Dasyra et al. (2009), while the full IR SEDs are compiled in S12.

Using the IRS spectra for both ULIRG samples, we estimate r_7.7, r_9.7, f₆, f₁₅, and f₃₀ in the same manner as for the simulated SEDs. However, the observed spectra have a more limited spectral coverage. The observed coverage of the IRS spectra varies between 5.2 and 38 μm, or 14–38 μm if only the Long-Low modules are used, as is the case for most of the higher-z sample. This means that for the local ULIRG sample, the required rest-frame 6–30 μm is covered by the available IRS spectra, while for the 24 μm-bright sample, the 15 μm and 30 μm continuum points are outside the IRS coverage above z ∼ 1.5. Below z ∼ 1.5, the 6 μm continuum point is often outside the IRS coverage (depending on whether or not the Short-Low module is available), and even the 7.7 μm point can be outside the IRS coverage below z ∼ 0.8. We compensate using the empirical SED model fits from S12, which necessarily introduces additional systematic uncertainty. These extrapolations are constrained below or above the IRS spectral coverage with IRAC 3.6, 4.5, 5.8, and 8 μm, and MIPS 70 μm photometric points, ensuring that the mid-IR colors obtained are reasonably accurate.

We choose to exclude the z ≲ 0.9 sources from this sample in order to avoid all cases where the 7.7 μm PAH feature is not covered by the IRS spectra. This step also removes all sources whose overall IR luminosities (L_3–1000) place them in the LIRG rather than ULIRG category. We also exclude the z > 2.5 sources, for which the 9.7 μm silicate absorption feature is too poorly defined, being only partially covered by the IRS spectra. This leaves a total of 118 sources. We explore determining the above quantities with and without interpolating the IRS spectra onto the much lower resolution simulated SED wavelength array and we find no significant difference in the overall results.

Lastly, in four cases, the silicate feature is saturated, making the measured r_9.7 value a lower limit. The saturation flag is triggered when the mean of the IRS spectra between 9.0 and 10.4 μm is less than or equal to the standard deviation in the same region.

In Section 6, we discuss the rest-frame near-IR as another means of diagnosing AGN power. For observational comparisons in this regime, we use Spitzer IRAC photometry for the high-redshift sample, and for the low-redshift ULIRGs, we interpolate from Two Micron All Sky Survey (2MASS) H- and K-band photometry (using total magnitudes). The unknown aperture corrections add some uncertainty here; however, ULIRGs are typically compact enough in the near-IR relative to the ∼2.5 arcsec 2MASS beam that aperture effects should be small (Surace & Sanders 1999; Surace et al. 2000). For several sources, we use 2.5–5 μm spectra (Imanishi et al. 2008; Sajina et al. 2009) from AKARI for more accurate near-IR fluxes.

In Figure 3, we explore how the mid-IR diagnostics vary with time for our two fiducial simulations: "highly obscured" and "less obscured." We show spectral diagnostics along with quantities of interest such as the overall IR luminosity, the SFR, and the AGN fraction, L_agn/L_bol. We assume our default ISM model and default AGN torus model, but we show the effect of varying these choices in Section 4.2 and Section 4.4, respectively. All curves shown are the median over the viewing angle. The scatter introduced by the uncertainty in the viewing angle is addressed in Section 4.3.

Zoom In Zoom Out Reset image size

To guide the discussion, we indicate the four times from the panels in Figure 1 with vertical lines for the highly obscured (left column) and less obscured (right column) mergers. The progenitors are still separated at time t₁, which precedes the period of final coalescence and peak L_IR by ∼200 Myr. The maximum SFR and AGN power occur at t ∼ t₂–t₃, after the galaxies have merged, but the central sources are still largely obscured by dust. At t₄, up to ∼300 Myr later, the SFR has plummeted but the AGN fraction is still elevated, and the remnants are on their way to being ellipticals.

Qualitatively, we find that our simulated mid-IR diagnostics behave largely as expected, but not always. Below, we address each of the diagnostics in turn.

The PAH strength, r_7.7, drops at the time of coalescence (t₂–t₃), reaching its lowest values at t₃ when the AGN fraction is highest. Pre-coalescence, the PAH strength is weakest for the model with strongest AGN (AGN×10). Post-coalescence, the PAH strength again increases but does not return to its pre-coalescence levels. No significant "coalescence" dip in the PAH strength is seen in the less obscured case. The level of absorption by silicate grains, r_9.7, is greatest when the AGN is most luminous, which is also when the relative strength of the PAH emission is weakest. These extrema are relatively sharp and confined to the period of final coalescence at ∼t₂–t₃, although the AGN luminosity fraction remains elevated for ≈1 Gyr. Overall, the shape of the PAH strength curve in both the highly obscured and less obscured cases does not bear a strong resemblance to the AGN fraction curve nor to the SFR curve.

In the highly obscured case, the f₁₅/f_5.7 slope is reddest when the AGN is most luminous, a somewhat counter-intuitive result since the 6 μm continuum is generally believed to be enhanced by AGN torus emission. This result likely owes to the fact that the time when the AGN fraction is greatest corresponds to the time of greatest obscuration, when the 9.7 μm silicate absorption feature is deepest. This obscuration would lead to significant absorption of the observed 6 μm continuum, reddening the f₁₅/f_5.7 color. However, this may result from the simulations not having sufficiently clumpy ISM to allow for lines of sight directly to the AGN torus. This slope is also reddened at t₂, the time of peak SFR activity, as expected due to enhanced H ii-region type emission. For the less obscured case, the f₁₅/f_5.7 color is reddened throughout the merger (t₁–t₄), likely the result of the enhanced SFR, but lacks a sharper reddening at the time of coalescence. This supports our view that the "coalescence" reddening seen in the "highly obscured" case is a direct result of galaxy dust attenuation.

Similarly, the f₃₀/f₁₅ slope is effectively constant for the duration of the merger for the less obscured case, but is significantly reddened at the time of coalescence in the highly obscured case. At a given time, this slope is bluer when the AGN source is more luminous (e.g., AGN×10), and in contrast to f₁₅/f_5.7, f₃₀/f₁₅ does not experience a sharp peak at t₃ for the most powerful AGN in the highly obscured merger.

Overall, our simulations indicate that the dependence of any of these mid-IR spectral diagnostics on the AGN fraction is time-dependent and non-linear and is affected by the properties of the host galaxy. The dependence on the overall level of obscuration, as parameterized by the silicate absorption depth, is to some degree stronger. For example, these indicators vary much less when there is less obscuration.

When the highly obscured and less obscured simulations have the same AGN fraction, the highly obscured case shows much sharper coalescence-stage SED reddening and deepening of the silicate feature. These effects on the SED are therefore the direct result of the host galaxy dust attenuation. This supports earlier observational evidence (Lacy et al. 2007; Sajina et al. 2007; Juneau et al. 2011) that in at least some AGN-dominated, deep silicate absorption sources, the absorption could arise in the host galaxy rather than the AGN torus.

In this section, we address the effect of varying both the dust structure in the ISM (clumpy versus smooth) and its composition (MW-like versus SMC-like). In Section 2.4, we described the two ISM treatments we apply to the highly obscured merger. Broadly speaking, the "default ISM" case assumes the dust and gas associated with both hot and cold ISM phases are smoothly distributed. The "alternate ISM" assumes the cold phase gas is in clumps sufficiently dense to have a negligible volume filling factor, so that the dust it contains does not contribute to the overall attenuation. In Figure 4, we explore these effects on the mid-IR SEDs.

Zoom In Zoom Out Reset image size

When switching to the alternate ISM assumption for the highly obscured simulation, the magnitude of r_9.7 is reduced by ∼0.3 dex, and the color f₃₀/f₁₅ by ∼0.7 dex, while r_7.7 and f₁₅/f_5.7 are mostly unchanged. In addition, the total IR luminosity (first row of Figure 4) is negligibly affected by this change. Therefore, the differences from this ISM structure assumption arise only because the source energy is redistributed in a different way, and not because of changes in the total amount of energy that is attenuated. This implies that there is less dust self-absorption in the alternate ISM model, which typically experiences less attenuation than the default case, so the source radiation that is reprocessed into the near-IR and mid-IR is not as often reprocessed again to longer wavelengths.

With SMC or MW dust, the IR SEDs at a given time during the starburst have similar slopes f₃₀/f₁₅ and f₁₅/f_5.7 (see also Figure 1). Furthermore, L_IR, a measure of the amount of attenuated light, is nearly identical in these cases. By contrast, the SEDs show very different spectral characteristics r_7.7 and r_9.7, reflecting the difference between MW and SMC grain composition. The SMC dust is assumed to contain relatively more silicate than carbonaceous grains, leading to virtually non-existent PAH emission r_7.7 and correspondingly stronger silicate absorption r_9.7.

In addition to the (evolving) source SEDs, the observed SED depends on the evolving relative distributions of sources and dust. This evolution changes the relative obscuration of AGN and stellar sources, affecting the behaviors of mid-IR indicators depending on whether the luminosity is dominated by SF or AGN, because the AGN is generally much more highly obscured than the stellar component. In the Appendix, we use toy models to confirm and further examine the behavior of the simulated diagnostics.

Figure 5 shows the time-dependence of variation in the mid-IR quantities with viewing angle, focusing on the highly obscured merger. The camera angle scatter is much lower (median S < 0.05 dex) for the less obscured merger. We show the same assumptions regarding AGN strength and ISM models from Figure 4, and we plot the scatter S, in dex, of the mid-IR quantities as a function of time. We define S(q) to be the normalized median absolute deviation (MAD; e.g., Mosteller & Tukey 1977; Beers et al. 1990) of the quantity q among the seven Sunrise viewing directions at a particular time during the simulation. First, we calculate the median M_q of q, and then the seven values |q − M_q|, the median of which we define to be the MAD of q. We define S(q) = MAD/0.67 so that if q is drawn from a Normal distribution with standard deviation σ, then S is an unbiased estimator of σ.

Zoom In Zoom Out Reset image size

At t ≲ t₂, S < 0.1 dex for all models, so the SED shape and spectral features are similar in all directions. Generally speaking, the scatter S reaches a maximum (∼0.2–0.4 dex) for conventional AGN diagnostics during t₃ to t₄. These times are also when the AGN luminosity fraction and IR luminosity are highest. Thus the increased viewing-angle scatter implies a large AGN contribution to the IR SED. This is because the AGN is effectively a point source, so the observed SED is sensitive to attenuation along a single line of sight, which can vary widely. Moreover, the sudden change in the extinction curve at 9.7 μm causes r_9.7 to have an exaggerated dependence on the line-of-sight attenuation. r_7.7 experiences much less viewing angle scatter than r_9.7, S(r_7.7) ∼ 0.0 versus S(r_9.7) ∼ 0.2 dex, because the extinction is the same for both the PAH features and 7.7 μm continuum. Turning off AGN emission leads to essentially no viewing angle variation in either quantity.

Figure 6 presents the SED from seven viewing angles at time t₃ (peak AGN) for our highly obscured merger, with AGN × 1 and fiducial ISM models. This gives us a clearer physical picture of what is driving changes to the mid-IR diagnostics during coalescence. The generally higher level of obscuration at these times leads to the enhanced silicate absorption, causing the dip in median r_9.7 at t₃ (Figure 3). However, at a fixed time during this period of very high attenuation, like the one shown in Figure 6, the silicate absorption feature is stronger along lines of sight in which the AGN is less obscured. In other words, the silicate feature only correlates positively with AGN attenuation until the level of attenuation is so extreme that the emergent stellar SED dominates the mid-IR emission associated with the AGN. This phenomenon is derived in the Appendix for a toy model of r_9.7.

Zoom In Zoom Out Reset image size

In Figure 7, we consider the effect of the intrinsic mid-IR AGN SED shape on the final observed SED. We show the manually boosted model (AGN×10) in order to maximize the possible effect. We demonstrate that the intrinsic AGN mid-IR spectrum has no impact on the observed mid-IR SED at t₂–t₃ of the highly obscured merger, during the starburst and AGN peak at final coalescence when the central densities are highest. We use our default AGN SED, the template by Hopkins et al. (2007c), as well as two templates by Nenkova et al. (2008). These templates reflect obscuration by clumpy tori inclined by i = 0° and i = 90°, spanning the range of derived mid-IR torus SED shapes.

Zoom In Zoom Out Reset image size

With our default ISM assumption, dust obscures this central source at wavelengths as long as λ ∼ 20–50 μm. Thus the flux at λ ∼ 2–5 μm, conventionally attributed to AGN tori, only indicates direct AGN emission when the central source is sufficiently unobscured (e.g., t₁ and t₄). In the less obscured simulation, the AGN torus SED has a large effect throughout the merger, as seen by the near-IR slopes of the SEDs in the left column of the "Less Obscured" panels in Figure 1.

In Section 4, we presented the time evolution of the observed mid-IR indicators and noted that their dependence on L_agn/L_bol is both time-dependent and non-linear. We showed that our mid-IR spectral diagnostics strongly depend on the dust model and can be highly scattered as a function of viewing angle. How could we hope to observationally discern the level of AGN in a given galaxy based on its mid-IR spectrum, given that typically no information on the merger stage is available, the dust model is even less constrained, and we will never have more than one viewing perspective? Here we attempt to do just that by using the simulation results as a diagnostic tool.

In Figure 8, we present mid-IR diagnostic diagrams that have been used in the literature for the purpose of disentangling the role of AGN and starburst activity in dusty galaxies. Specifically, using our definitions of the PAH strength, silicate feature depth, and mid-IR colors, we present slightly modified versions of the "Laurent" diagram (Laurent et al. 2000), the "Spoon" diagram (Spoon et al. 2007), the "Veilleux" diagram (Veilleux et al. 2009), and the "Lutz" diagram (Lutz et al. 2004). The last diagram presents the relation between the 6 μm luminosity and X-ray luminosity for unobscured or mildly obscured AGNs. The left-hand column of Figure 8 presents the highly obscured simulation, the middle column presents the less obscured simulation, and the right-hand panel presents real galaxy data that we discuss in more detail in the following section.

Zoom In Zoom Out Reset image size

For the simulation columns in Figure 8, we plot the IR quantities in the following way. Each panel is divided into a grid of 25 × 25 bins, each of which contains zero or more simulated points. These points include all seven viewing angles at each timestep in t₁–t₄, which are sampled every 10⁷ yr. We use our fiducial ISM settings and plot the three AGN strengths (×10, ×1, and ×0) in order to span a full range of L_agn/L_bol in our experiments. Bins containing zero of these points are white, and bins with one or more point are shaded according to the median value of L_agn/L_bol among the points it contains, as indicated by the color bar: darker/blue regions have L_agn/L_bol ∼ 1, while lighter/yellow regions have L_agn/L_bol ≲ 0.1. We plot histograms (on a logarithmic scale) along the edges of the first three rows to show the relative timescales of the model populations giving rise to these diagrams.

We caution that our plotting and selection here are not meant to directly characterize the observed populations of ULIRGs that we plot in the right column; these are included primarily as a reference and plausibility check. Instead, from this coding we seek to identify which regions, if any, isolate powerful AGNs during these mergers, regardless of the length of time the observed properties may be visible. A larger sample of models will be needed in order to fully constrain the expected observed distribution of L_agn/L_bol at each location in these diagrams, but we can begin to highlight some basic trends implied by the simulations.

Veilleux et al. (2009) analyzed the standard evolutionary scenario for ULIRGs (Sanders et al. 1988), which broadly speaking proceeds through three stages. In "Stage 1," the IR emission of a gas-rich merger is dominated by SF, leading to large PAH EWs and relatively small silicate absorption EWs. As the merger proceeds to "Stage 2," r_9.7 increases as the central source is buried in the gas-rich starburst, and the PAH emission decreases relative to the IR continuum that is boosted by the AGN, which remains sub-dominant to the starburst. Finally, at "Stage 3," the silicate absorption and PAH emission both subside as the starburst fades, the obscuring medium is cleared or consumed, and the ULIRG is completely or nearly AGN-dominated. Overall, Veilleux et al. (2009) concluded that AGNs contribute ∼35%–40% to the IR luminosity of ULIRGs, and that ULIRGs experience significant scatter from this general pattern.

In Figure 8, the highly obscured model points roughly fill the space spanned by observed sources. Our procedure does not mandate such consistency, and so we interpret this as one important check of the modeling technique. In particular, ignoring other diagnostics (such as SED slopes), this simulation reproduces the observed spread in r_9.7 and r_7.7. In Table 2, we summarize a few numerical results of our deeply obscured simulation. We divide the merger into three stages, t₁–t₂ (pre-coalescence), t₂–t₃ (coalescence), and t₃–t₄ (post-coalescence), and report the median values of L_agn/L_bol, r_9.7, and r_7.7 calculated in those periods.

Generally, the SED properties along the model timeline follow those of the merger stages by V09, albeit with an expectedly high amount of scatter. Local ULIRGs have a large amount of silicate absorption (r_9.7 ∼ 0.7), suggesting they may correspond to the modeled points at t₂ < t < t₃ where r_9.7 peaks. However, the AGN fractions in models depend sensitively on the dust model and AGN strength choices, and so it is unclear how well they may correspond to the ∼30% reported in V09. This value is consistent with our AGN × 1, MW dust simulation, but is much lower than the ∼70% from the "combined" model set, in which each of the AGN strengths are given equal weight. The more realistic choice is uncertain, but the AGN × 1 case has trouble reproducing observed objects with r_9.7 < −1. However, SMC dust with lower AGN strengths can also provide that feature, and observations suggest that SMC-like dust may be more appropriate during bright starburst or AGN phases, as small grains are preferentially depleted (Gordon et al. 1997; Hopkins et al. 2004). This may imply that SMC dust is more appropriate to use in this situation; in any case, it highlights the difficulty in determining the path of the AGN fraction in such objects if their dust grain populations are changing significantly.

In Figure 8, swaths of strong and weak AGN exist, but most of these projections do not isolate AGN in simple ways. For example, while it is true that sources with strong starbursts and weak AGNs typically have stronger PAH emission and redder f₃₀/f₁₅ as discussed in Veilleux et al. (2009), strong AGNs occupy all values of f₃₀/f₁₅, f₁₅/f_5.7, and r_9.7. Furthermore, the PAH strength depends sensitively on the dust model we assumed (see Figure 4), and the r_7.7 values of low L_agn/L_bol points depend on the obscuration level of the simulation.

In addition, we compare the relative fluxes emerging at hard X-ray and mid-IR wavelengths in the bottom row of Figure 8. We compute model X-ray fluxes as in Section 2.5. Bauer et al. (2010) found that z > 1 ULIRGs that are otherwise thought to be AGN-dominated (with low r_7.7) are undetected at X-ray wavelengths, suggesting that they are at least mildly Compton-thick. Moreover, they presented Compton-thick examples from the literature at z ∼ 0, which lie ∼2 dex below the relation for less obscured AGNs (cf. Hönig et al. 2010). Many model galaxies with powerful AGNs lie in the same region of this diagram, and are reflected in the dark blue squares mixed in with the lighter squares at L₆ ∼ 10¹¹–10¹² L_☉. These are the same sources in the low-r_7.7 (high r_9.7) tail of dark points in the upper panels, implying that this class of Compton-thick sources is a short-lived phase of our highly obscured model. The bulk of the simulation points have X-ray fluxes ∼0.5 dex below the gray Type 1 AGN band, reflecting our choice to not analyze completely unobscured models.

There are some properties for which the model values are not realized by observed sources. Generally speaking, the differences are most apparent in the mid-IR SED slopes f₃₀/f₁₅ and f₁₅/f_5.7. The most AGN-dominant model points are ≈0.5 dex redder in f₁₅/f_5.7 than any observed object. In bulk, the model f₃₀/f₁₅ values are larger by ≈0.3 dex than the observed distribution, and there are few model galaxies with log r_9.7 ∼ 0 and log r_7.7 ∼ 0, in contrast to those observed. However, this is sensitive to our dust grain model (Section 6.2), and may also reflect our not having modeled AGNs with no galaxy obscuration. Our calculations are necessarily uncertain, so in Section 6.3 below, we summarize some of the theoretical limitations and potential implications of these mismatches.

As we have seen in Section 4, when sufficiently obscured, a powerful AGN has an SED shape that is the same as a dusty pure starburst (see also Younger et al. 2009; Narayanan et al. 2010a). Figure 8 shows high L_agn/L_bol squares occupying a large fraction of the r_7.7, r_9.7, f₃₀/f₁₅, and f₁₅/f_5.7 values. We also saw that the model SEDs can change drastically on timescales ∼10⁸ yr. A further challenge is that the implied L_agn/L_bol, given a set of SED properties, depends sensitively on assumptions about the dust model and AGN strength. An ideal indicator is one that identifies powerful AGNs regardless of the host galaxy's content or state. An example is hard X-ray emission which can penetrate even very large optical depths, and in Section 6, we use our simulations to suggest other ideal indicators that have not been easily accessible with current data, but may be useful for future mid-IR studies with the JWST.

Ultimately, we seek an indicator that can be used to determine the AGN contribution to IR-luminous galaxies across the full range of obscuration by galaxy dust. Such a tool could simplify efforts to accurately track the growth of SMBHs and assess their impact on the ISM in diverse and evolving host galaxies.

Near-infrared colors effectively identify AGNs when their IR SED resembles a power-law, leading to enhanced emission beyond the stellar bump at ≳ 2 μm. This yields a high value of our example near-IR color f_4.5/f_1.6, for which the numerator is a proxy for torus emission and the denominator is a proxy for stellar mass. Such techniques have been used to successfully identify high-redshift AGN-dominated ULIRGs using Wide-field Infrared Survey Explorer photometry (e.g., Stern et al. 2012; Eisenhardt et al. 2012; Wu et al. 2012). This near-IR reddening operates in powerful AGNs for a wide variety of attenuation. However, as we saw in Figure 7, on short timescales, the near-IR SED shape can change from a power-law AGN to one that resembles a starburst-dominated ULIRG when the obscuration increases drastically. In these high-obscuration situations, r_9.7 is sensitive to the AGN power, and so combining these indicators into a generic diagnostic may be feasible.

In Figure 9, we plot two-dimensional diagrams with our near-IR color f_4.5/f_1.6 on the y-axis and r_9.7 on the x-axis. We have included our fiducial ISM assumptions for each mergers and all three AGN strengths: AGN × 10, AGN × 1, and AGN × 0. We note that these settings are chosen by hand in order to more broadly span the range of L_agn/L_bol produced by our calculations, and are justified by the uncertainty inherent in the simulation methods and limited range of cosmologically motivated models we consider (see, e.g., Section 2.4). Therefore, these model diagrams do not directly correspond with any realistic ULIRG populations; we plot two samples in the right panel primarily as a reference.

Zoom In Zoom Out Reset image size

The light/yellow shaded region corresponds to low AGN emission (L_agn/L_bol ≲ 0.1) and is confined to a small region (∼0.4 by 0.4 dex) of this parameter space centered at (log r_9.7, log f_4.5/f_1.6) ≈ (− 0.32, −0.5). In both the highly and less obscured simulations, the value log f_4.5/f_1.6 ∼ −0.5 corresponds to the SED shape of purely stellar emission, which can be seen in the gray dotted curves in Figure 1 at 1 μm < λ < 5 μm. The value log r_9.7 ∼ −0.32 appears to be the typical SED produced by ordinary (non-bursting) SF in our simulations. Numerically, this value may be set by our simplistic continuum estimate (Section 3.1), which is likely affected to some extent by the nearby PAH features. Both of these "zeropoints" do not depend strongly on the dust model, and are roughly the same in the highly and less obscured mergers.

Away from this point, the median AGN fraction varies in a semi-regular pattern. A redder f_4.5/f_1.6 implies a contribution at 4.5 μm from the relatively brighter AGN source, as expected from conventional hot dust/torus interpretations. The absorption diagnostic r_9.7 does not vary strongly in the less obscured case, as we have seen in Section 4.1, and so the locus of increasing L_agn/L_bol traces roughly a vertical line in the second column of Figure 9.

By contrast, in the highly obscured simulation, a significant area of the log r_9.7 ≲ −0.5 regions of the diagram is occupied by squares with L_agn/L_bol ≳ 0.1. In the lower left quadrant of each panel in the left column of Figure 9, the dust column is large enough to produce a significant silicate absorption feature and to attenuate the near-IR continuum associated with AGN-heated dust, leading to log f_4.5/f_1.6 ≲ 0. In the upper left quadrant, there are several squares with both a strong silicate feature and a red near-IR continuum. At a fixed L_agn/L_bol, these two quantities are correlated in the sense that increased silicate absorption implies more attenuation of the near-IR AGN-heated continuum relative to the starlight, giving a bluer f_4.5/f_1.6, leading to the slanted or semi-elliptical pattern in the areas of constant shading (L_agn/L_bol) in Figure 9. These qualities are the same when switching to the "alt. ISM" model for the highly obscured merger (not shown).

We seek to summarize these trends and estimate L_agn/L_bol directly from these variables. Visually, the areas of constant L_agn/L_bol form a semi-elliptical pattern that can be approximated by an analytic formula to potentially estimate L_agn/L_bol with low scatter. For MW dust, we consider the formula

$$D_{\rm MW} = \sqrt{ \left(\log r_{9.7}+ 0.32 \right)^2 + \left(\left[\log \frac{f_{4.5}}{f_{1.6}} + 0.5 \right] \frac{2}{3} \right)^2}.$$

The effectiveness of combination D_MW is evident in the top row of Figure 9: the map between color/shading (L_agn/L_bol) and these axes, f_4.5/f_1.6 and r_9.7, appears to be a well defined and slowly varying function of few parameters. In the Appendix, we confirm the rough behavior of this diagnostic from toy models.

We caution that we have analyzed only a handful of specially chosen models: AGN × 10, AGN × 1, and AGN × 0 in two mergers that likely occupy two very different extremes of the ULIRG population. Therefore, the shaded points in Figures 8 and 9 are likely to change with a larger and more representative sample of models. Nevertheless, it is encouraging that a clear trend in the f_4.5/f_1.6–r_9.7 plane exists for both our highly and less obscured mergers (in contrast to the diagrams in Figure 8), and for either the default or alternate ISM model in the highly obscured merger. This can be understood as a result of the relatively simple ways that the two axes of Figure 9 reflect AGN-associated emission and its obscuration. In Section 6.2, we show that this visual trend we identify in Figure 9 defines a tight relationship between L_agn/L_bol and D_MW at all times t₁–t₄ when the dust model is constrained.

However, we note that shape of the trend in the f_4.5/f_1.6–r_9.7 plane depends on our dust model: the bottom row of Figure 9 shows the same simulation using SMC dust instead of MW dust. Importantly, a nice mapping in the above sense still exists, it is just skewed toward higher values of r_9.7. Therefore, we can define a similar diagnostic,

$$D_{\rm SMC} = \frac{5}{8} \sqrt{ \left(\log r_{9.7}+ 0.32 \right)^2 + \left(\left[\log \frac{f_{4.5}}{f_{1.6}} + 0.5 \right] \frac{5}{4} \right)^2},$$

centered around the same point as D_MW, but with a different axis ratio and scaling.

In order to select between these estimators, we must first constrain the dust properties; in other words, we must jointly constrain the dust and L_agn/L_bol. One option is to use the PAH emission: assuming SMC dust, r_7.7 ∼ 0 for all model points instead of ∼0.3 for MW dust (see Figure 4). We demonstrate this procedure in Figure 10: while D_MW does not yield a low-scatter predictor of L_agn/L_bol when applied to a system with SMC dust, it does broadly select stronger AGN. Moreover, at D ≳ 0.2, sources with these different dust models separate based on their PAH strength, r_7.7, permitting either a direct choice between the formulae for D_MW and D_SMC or an intermediate case. This suggests the following procedure:

• 1.  
  Estimate D_MW.
• 2.  
  Estimate the PAH strength (here, r_7.7).
• 3.  
  Select dust model X based on the value obtained in step 2, by, e.g., checking whether at the D_MW value, r_7.7 lies closer to the MW or SMC points.
• 4.  
  Calculate D_X and use it to estimate L_agn/L_bol.

Zoom In Zoom Out Reset image size

With a known dust model, diagnostic D robustly predicts L_agn/L_bol regardless of the physical conditions of the ULIRG source. We compare the ability of diagnostics D_MW and D_SMC to predict L_agn/L_bol to that of hard X-ray fluxes in Figure 11.

Zoom In Zoom Out Reset image size

Using a least absolute deviation method, we fit the relation L_agn/L_bol = a + bQ for three quantities Q (D_MW with MW dust, and both D_MW and D_SMC with SMC dust) and each of our two ISM models for the highly obscured simulation, and compile the results in Table 3. To measure the scatter in the correlations' residuals, we used the same scale estimator S defined in Section 4.3, which is generally insensitive to outliers. In both cases, L_2–10 keV and each D predict L_agn/L_bol with roughly the same median scatter, S ∼ 10%. However, we find that L_2–10 keV experiences a failure rate (>3σ) two to three times higher than D in the highly obscured calculation. Both quantities predict L_agn/L_bol equally well under the alternate ISM treatment and for the less obscured candidate (not shown), with S < 10% and failure rates <1%.

Several simulated points at L_agn/L_bol ≈ 0.8 fail in both L_2–10 keV and D_MW. In these cases, corresponding to the most highly obscured lines of sight in our AGN × 10 calculation from Figure 6, the mid-IR is so heavily attenuated that no r_9.7 signature from the central source survives in the observed SED (see also Section 4.3).

Diagnostic D requires information about the galaxy's stellar mass (rest-frame λ  ∼  1–2 μm), limiting its immediate application to large new samples. However, the rest-frame wavelength range D requires will be covered by the JWST, an observatory whose Mid-Infrared Instrument (MIRI; e.g., Wright et al. 2004) permits measurements of r_7.7 and r_9.7 for sources at z ≲ 1. The tight correlation we find between D and L_agn/L_bol reflects the possible benefit in combining information about r_9.7 and r_7.7 with near-IR color diagnostics from JWST itself (e.g., Messias et al. 2012), or by matching to sources from other IR surveys (e.g., Stern et al. 2005, 2012). For the two example SEDs in Figure 12, whose known AGN fractions are both L_agn/L_bol = 0.33, this procedure estimates L_agn/L_bol = 0.33 ± 0.09 and L_agn/L_bol = 0.38 ± 0.09 (uncertainties are from the minimum theoretical scatter in the correlations) for the high and low obscuration examples, respectively.

Zoom In Zoom Out Reset image size

A concern with the present method is that we cannot be sure how our few idealized simulations fit into a cosmological context. We have mitigated this issue by presenting two rather different ULIRG candidates, by varying the relative power between the AGNs and SF, and by varying the assumptions acting on dust obscuration and emission. We are planning to build a library of model SEDs for a much wider variety of systems, enabling more robust studies across IR-luminous galaxy populations. Moreover, we have focused here on AGN signatures at mid-infrared wavelengths, and information at longer wavelengths (cf. Younger et al. 2009), such as the location of the far-IR peak, could further constrain the AGN fractions in such sources (e.g., S12). In an upcoming paper (E. Roebuck et al., in preparation), we include the far-IR peak predictions from these simulations in interpreting the properties of z ∼ 0.3–3 starbursts and quasars (the S12 sample).

The simplest interpretation of Figure 8 is that analogs of our "highly obscured" and "less obscured" examples are both realized in observed sources, in addition to both MW and SMC dust: all four models are required to cover the wide observed distributions of mid-IR color, PAH emission, and silicate absorption. A higher ratio of silicate to carbonaceous dust grains increases r_9.7 while shrinking r_7.7 (Figures 4), without changing f_4.5/f_1.6 (Figure 9), permitting estimates of L_agn/L_bol under a variety of physical conditions, dust models, and optical depths. Although different ISM or dust models can produce broadly similar mid-IR trends, a separate measure of the AGN power, such as X-ray flux, is essential to further constrain these models. Furthermore, the recent discovery of high-redshift, hyper-luminous, hot-dust ULIRGs (Eisenhardt et al. 2012; Wu et al. 2012) suggests that we may lack a key phase in these few simulations, in which feedback initiates a transition between scenarios represented by our "highly obscured" and "less obscured" SED models.

Several assumptions regarding the dust RT may lead to discrepancies between our model sources and real sources. In particular, from Figure 8 we have identified two failures to which these assumptions likely contribute: (1) the mid-IR colors of the most AGN-dominated highly obscured sources are redder than their real counterparts (although bluer at λ ≲ 5 μm; see Figure 9), and (2) we do not simulate enough sources with little signs of PAH emission or silicate absorption: r_9.7 ∼ 0 and r_7.7 ∼ 0.

The IR calculation in Sunrise distributes the energy absorbed by small PAH grains equally between emission in thermal equilibrium and in modes giving rise to the distinctive mid-IR PAH lines (Jonsson et al. 2010). Therefore, any change to this assumption or a more sophisticated calculation, such as thermally fluctuating small grains (e.g., Guhathakurta & Draine 1989), will affect the dust continuum shape and PAH strength. In addition, we do not currently model dust destruction, and our model for dust production is crude, both key areas where near-future progress may be tenable. By not allowing our radiation to destroy dust grains or alter the dust grain distribution (or changing more generally the gas distribution owing to realistic feedback), we may be preventing ourselves from obtaining truly "power-law" AGN SEDs with no mid-IR dust features that may arise after AGN feedback disrupts the ISM.

We have purposefully focused on how galaxy dust can affect the fixed SMBH source SED (with minor variations in Section 4.4), rather than how the physics on smaller scales play out in observations of ULIRGs. While high levels of galaxy dust can reproduce a number of the properties of obscured quasars, the evolution of the central regions during these times of maximum obscuration will depend on the behavior of gas structures very close to the SMBH. This behavior may set the minimum scatter to be observed in correlations between these indicators. Using multiscale simulations, Hopkins et al. (2012a) showed how gravitational instabilities propagate down to scales of ∼1–100 pc during large-scale gas inflows. They demonstrated that a number of properties of obscured AGNs attributed to a torus are explained by such dynamically-evolving accretion structures. Furthermore, these simulations suggest that the angular momentum vector of the small-scale "torus" can often be misaligned with the angular momentum of gas on galaxy scales, and the relative orientations can vary significantly in time (Hopkins et al. 2012b). Thus, the attenuation from the host galaxy may vary more significantly than our current simulations suggest. Therefore, we are motivated to develop methods of integrating knowledge of the impact of these different scales on the emergent SEDs as a means for further interpreting the evolution of IR-luminous galaxies.

Recent work on merger simulations at higher resolution (Hopkins et al. 2013) and with different hydrodynamics solvers (e.g., Arepo: Springel 2010; C. Hayward et al., in preparation) shows that the average global evolution in simulations like the ones we analyze here, and hence conclusions made based on them, are robust. However, such advances will enable more sophisticated treatments of feedback processes and hydrodynamics in gas with a complex phase structure and at smaller scales; therefore, future SED modeling could exploit these improvements to produce ever more detailed comparison libraries.