A JOINT MODEL OF THE X-RAY AND INFRARED EXTRAGALACTIC BACKGROUNDS. I. MODEL CONSTRUCTION AND FIRST RESULTS^*

The CXB was discovered almost half a century ago (Giacconi et al. 1962), and is now understood to be dominated by the integrated emission from gas being heated and accreted onto SMBHs at the centers of active galaxies. The pure blackbody spectrum of the cosmic microwave background places strong upper limits on the contributions to the CXB from diffuse hot intergalactic gas (10⁻⁴; Wright et al. 1994). Deep extragalactic surveys have now resolved a significant fraction (∼70%) of the background light at <10 keV into discrete sources (Hasinger et al. 1998; Mushotzky et al. 2000; Alexander et al. 2003; Worsley et al. 2005), with the resolved fraction decreasing with the increasing energy bands. The exact fraction also depends on the absolute normalization of the CXB spectrum whose dispersion reaches ∼30% at low energies (<10 keV) and ∼10% above 10 keV among different studies (e.g., Marshall et al. 1980; Gruber et al. 1999; Vecchi et al. 1999; Revnivtsev et al. 2003; Ajello et al. 2008; Churazov et al. 2007; Moretti et al. 2009; Türler et al. 2010).

The CXB spectrum peaks around 20–30 keV and drops toward both lower and higher energies, and its slope at <20 keV is harder than those of Seyfert 1 galaxies. In order to reproduce this slope and turnover around 20–30 keV, the general strategy is to invoke the SED of a type 1 AGN along with a varying H i column density (Setti & Woltjer 1989; Madau et al. 1994; Comastri et al. 1995), consistent with the unified scheme for AGNs (Antonucci 1993) where Seyfert 2 galaxies are obscured variants of Seyfert 1 galaxies and exhibit harder X-ray spectra. Today CXB models not only reproduce the CXB spectrum (1–100 keV), but also the X-ray number counts, redshift distributions, and spectral properties in both soft and hard bands (e.g., Treister & Urry 2005; Gilli et al. 2007; Treister et al. 2009). It is generally accepted that the majority of X-ray AGN activity occurs in the obscured phase and that the evolution of X-ray AGNs is luminosity dependent, with the density of Seyfert-like AGNs peaking around z = 1 and quasars peaking at z = 2–3. However, there are still important uncertainties in the derived parameters, such as the type 1/type 2 AGN ratio and its variation with luminosity and redshift, and the amount of Compton-thick objects. Such difficulties are at least partly caused by the fact that the current X-ray surveys at <10 keV are not sensitive to heavily extinguished AGNs, especially Compton-thick objects, and overall number counts remain limited.

The CIRB was first directly measured in the 1990s by COBE (for a review, see Hauser & Dwek 2001). The energy of CIRB is comparable to the cosmic UV/optical energy (Dole et al. 2006), while the ratio in the local universe is only one-third (Driver et al. 2008). This implies an increase in dusty sources with increasing redshift. The slope of the CIRB in the far-IR/submillimeter range is flatter than that of local starburst galaxies, which provides further evidence for the strong evolution of IR sources whose peak emission is redshifted toward longer wavelength, flattening the CIRB spectrum at the longer wavelength.

Unlike the CXB model where the H i column density is considered to be the main factor driving the variation of the SED, the CIRB model typically varies the SED mainly as a function of IR luminosity. They either adopt a fixed SED at a given luminosity (Chary & Elbaz 2001; Le Borgne et al. 2009), or further consider the scatter of the SED at a given luminosity (Chapman et al. 2003; Valiante et al. 2009), or assume different IR populations that exhibit different SEDs and occupy different luminosity ranges (Xu et al. 2001; Lagache et al. 2003, 2004; Rowan-Robinson 2009; Franceschini et al. 2010; Gruppioni et al. 2011; Béthermin et al. 2011). The SMBH radiation, with typically a much warmer IR SED compared to that of SF regions, is energetically important in mid-IR surveys: the contribution of AGNs to the mid-IR background reaches noticeable fractions, e.g., 10%–20% in the 1–20 μm range by Silva et al. (2004) or 2%–10% in the 3–24 μm by Treister et al. (2006); classifications (AGN versus SF dominated) of 24 μm sources indicate that the fraction of objects whose mid-IR emission is dominated by AGNs ranges from ∼40% at f_24 μm > 5 mJy and drops down with decreasing fluxes but remains around 10% down to the faintest end f_24 μm ∼ 0.08 mJy (Donley et al. 2008; Fu et al. 2010; Wu et al. 2010; Choi et al. 2011). However, the 24 μm number counts can be reproduced equally well by models with AGNs (Valiante et al. 2009; Franceschini et al. 2010; Gruppioni et al. 2011) and without AGNs (Béthermin et al. 2011; Le Borgne et al. 2009), reflecting the limitation of the CIRB synthesis models in decomposing SMBH and SF emission. Models without AGNs can reproduce the effect of AGNs in the mid-IR simply by modifying the SF mid-IR SED or by adjusting the slope of LFs and their redshift evolution.

We represent the galaxy total IR LF as a double power-law distribution with 10 free parameters:

$$\begin{eqnarray} \Phi (L_{\rm TIR}, z) &=& \frac{dN}{{d\,{\rm log}L}\times dV} = \frac{\Phi _{*}}{(L/L_{*})^{\gamma _{1}}+(L/L_{*})^{\gamma _{2}}}\nonumber \\ && \times \left({\rm Mpc}^{-3}\, {\rm log}\,L_{\odot }^{-1}\right). \end{eqnarray} \tag{ 1 }$$

We adopt a simple independent luminosity and density evolution. The evolution is described as a polynomial function of redshift:

$$\begin{equation} {\rm log}\,L_{*}(z) = {\rm log}\,L_{*, 0} + \eta \times k_{1, l} + \eta ^{2}\times k_{2, l}+\eta ^{3}\times k_{3, l} \hspace{12pt} \end{equation} \tag{ 2 }$$

$$\begin{equation} {\rm log}\,\Phi _{*}(z) = {\rm log}\,\Phi _{*, 0} + \eta \times k_{1, d} + \eta ^{2}\times k_{2, d}+\eta ^{3}\times k_{3, d}, \end{equation} \tag{ 3 }$$

where η = log(1 + z). For this component of the model, there are 10 free parameters, including four in the luminosity evolution (logL_{*, 0}, k_{1, l}, k_{2, l}, and k_{3, l}), four in the density evolution (logΦ_{*, 0}, k_{1, d}, k_{2, d}, and k_{3, d}), and two slopes (γ₁ and γ₂).

Individual objects in the model are assumed to be experiencing both SF and SMBH accretion activities, so that L_TIR = L^BH_TIR + L^SF_TIR. For the formula of the SMBH energy fraction, we assume a bounded Gaussian distribution for the logarithm of the SMBH energy fraction in the total IR band f^BH_TIR = (L^BH_TIR/L_TIR), as motivated by works of Valiante et al. (2009):

$$\begin{eqnarray} &&P^{\rm BH}\left({\rm log}f^{\rm BH}_{\rm TIR}\right) =\frac{dN}{d\,{\rm log}f^{\rm BH}_{\rm TIR}}\propto {\rm exp}{\left({-}\frac{ \left({\rm log}f^{\rm BH}_{\rm TIR} - {\rm log}f^{\rm BH}_{\rm *, TIR}\right)^{2} }{2(\sigma ^{\rm BH})^{2}}\right)} \nonumber \\ \end{eqnarray} \tag{ 4 }$$

$$\begin{equation} {\int _{\rm min}^{\rm max}}P^{\rm BH}\left({\rm log}f^{\rm BH}_{\rm TIR}\right) d{\rm log}f^{\rm BH}_{\rm TIR} = 1, \end{equation} \tag{ 5 }$$

where the boundary of logf^BH_TIR is set to have the minimum and maximum at −6.0 and 0.0, respectively (for details, see Section 4).

This Gaussian distribution is assumed to be redshift- and luminosity-dependent. The mean value of the distribution (logf^BH_*) is a function of both total IR luminosity and redshift:

$$\begin{equation} {} {\rm log}f^{\rm BH}_{*} = {\rm log}f^{\rm BH}_{*, z} + p_{f^{\rm BH}}\times {\rm log}\frac{L_{\rm TIR}}{L_{*, f^{\rm BH}}}, \end{equation} \tag{ 6 }$$

where

$$\begin{equation} {\rm log}f^{\rm BH}_{*, z} = {\rm log}f^{\rm BH}_{*, 0} + \eta \times k_{1, d}^{\rm BH} + \eta ^{2}\times k_{2, d}^{\rm BH}. \end{equation} \tag{ 7 }$$

The width of the distribution is only luminosity-dependent:

$$\begin{equation} \sigma ^{\rm BH} = \sigma _{0}^{\rm BH} + p_{\sigma }\times \log \frac{L_{\rm TIR}}{L_{*, f^{\rm BH}}}. \end{equation} \tag{ 8 }$$

At z = 0 and L_TIR = $L_{*, f^{\rm BH}}$, the SMBH energy fraction is fixed at f^BH_{*, 0}, which is based on the sample of local 1 Jy ULIRGs (Kim & Sanders 1998). The advantage of ULIRGs compared to lower luminosity objects is that the SMBH contributions are important so that the distribution of the SMBH energy fraction can be derived to a high accuracy. For a sample of 74 ULIRGs, we derived the SMBH IR luminosity from the observed 6.2 μm EW and 6.0 μm continuum luminosity as measured by Veilleux et al. (2009):

$$\begin{eqnarray} L_{6\,{\mu }{\rm m}}^{\rm BH}&=&{\rm max}[(1.0-{\rm EW}_{6.2\,{\mu }{\rm m}}/{\rm EW}_{\rm star\ formation})\nonumber \\ && \times L_{6\,{\mu }{\rm m}}\times \exp (\tau _{6\,{\mu }{\rm m}}), 0], \end{eqnarray} \tag{ 9 }$$

where EW_6.2 μm is the observed 6.2 μm equivalent width, L_6 μm is the observed 6.0 μm continuum luminosity, EW_{star formation} is the equivalent width of pure SF and taken to be 4.2 μm (Veilleux et al. 2009), and τ_6 μm is the effective optical depth at 6 μm and set to be 10% of the measured 9.7 μm silicate depth using the Small Magellanic Cloud (SMC) extinction curve (Pei 1992). Given that the star-forming EWs still have a small scatter (e.g., Wu et al. 2010; S. Stierwalt et al., in preparation), we set the maximum value of the derived SMBH radiation in Equation (9) to be zero. Note that the aromatic feature EW relies on the method to measure it, e.g., sometimes the wing emission of the feature is included through extrapolating the analytic profile. As long as the same method is adopted, the SMBH energy fraction from the above equation should not differ significantly. The L^BH_6 μm is then corrected to the total SMBH IR luminosity by a factor of 1.6 using our own SMBH SED. Combined with the measured total IR luminosity, the final distribution of the SMBH energy fraction is shown in Figure 1. Our assumption of a Gaussian distribution is more or less correct given the derived shape. The fit with a Gaussian gives logf^BH_{*, 0} = −1.21. The median luminosity of this sample gives log$L_{*, f^{\rm BH}}$ = 12.3. Note that the above values do not depend significantly on the adopted EW_{star formation} since EW_6.2 μm is small in the majority of the objects. However, we did not use the above distribution to fix σ^BH₀ but set it as a free parameter because the width measured in Figure 1 is broadened by the intrinsic scatter in EW_SF and luminosity correction L_TIR/L_6 μm. We note that the median fraction derived here is significantly lower than the energy fraction of roughly 35%–40% as derived by Veilleux et al. (2009). The main reason is that their numbers are given in the bolometric band while ours are integrated from 8 to 1000 μm. While the emission is similar between these two bands for ULIRGs, a factor of ∼5 needs to be applied for the SMBH radiation.

Zoom In Zoom Out Reset image size

There are five free parameters in total, including two describing redshift evolution of the mean fraction (k^BH_{1, d} and k^BH_{2, d}), one characterizing luminosity dependence of the mean fraction ($p_{f^{\rm BH}}$), one for the scatter (σ^BH₀) at the total IR luminosity of 10^12.3 L_☉, and another for its luminosity dependence (p_σ).

For the IR SED of the SF emission, we adopted the luminosity-dependent SED shapes and associated dispersions. This set of SEDs is well characterized in the local universe; the uncertainty lies at high-z where the difference is already observed but not fully quantified.

3.3.1. Local Infrared SED Templates

It is now established that the IR luminosity is not the main driver of the IR SED shape, i.e., at a given luminosity, the IR color shows a large scatter (Chapman et al. 2003; Hwang et al. 2010; Rujopakarn et al. 2011a; Elbaz et al. 2011). Fortunately, this scatter has been well quantified. We thus adopted the observed trend between IR luminosity and color (f_ν(IRAS60 μm)/f_ν(IRAS100 μm)) as well as the associated scatter as established by Chapman et al. (2003) for local IRAS galaxies. The probability distribution of IR color at a given star-forming IR luminosity is described by a Gaussian function:

$$\begin{equation} P^{\rm SF\,color} = \frac{dN}{dC}= \frac{1}{\sigma _{c}\sqrt{2\pi }} {\rm exp}\left(-\frac{(C-C_{0})^2}{2(\sigma _{c})^2}\right), \end{equation} \tag{ 10 }$$

where C = log (f_ν(60 μm)/f_ν(100 μm)), σ_C = 0.065, and

$$\begin{equation} {} C_{0} = C_{*} + {\gamma }\times {\rm log}\left(1+\frac{L_{\rm TIR}^{^{\prime }}}{L_{*, C}}\right) - {\delta }\times {\rm log}\left(1+\frac{L_{*, C}}{L_{\rm TIR}^{^{\prime }}}\right), \end{equation} \tag{ 11 }$$

C_* = −0.35, γ = 0.16, δ = 0.02, L_{*, C} = 5.0 × 10¹⁰ L_☉, and L^'_TIR is the star-forming total IR luminosity from 3 μm to 1100 μm. For a given color of (f_ν(IRAS60 μm)/f_ν(IRAS100 μm)), we used the infrared SED template library of Dale & Helou (2002, hereafter DH02), which sorts the IR SED by that color. The template was built empirically first from IRAS/Infrared Space Observatory (ISO) 3–100 μm observations (Dale et al. 2001) and then extended up to submillimeter wavelengths based on SCUBA and ISO data (Dale & Helou 2002). Spitzer data (3–160 μm) confirm its ability to cover the overall range of SED variations and describe the median trend between different colors below 200 μm (e.g., Dale et al. 2007; Wu et al. 2010). Recent Herschel far-IR observations further demonstrated its power (Rowan-Robinson et al. 2010; Boselli et al. 2010; Dale et al. 2012). Nevertheless, observations suggest that this set of templates likely lacks the dynamic range in the SED variations above a rest-frame of 200 μm. For example, DH02 templates do not produce objects with rest-frame f_ν(Herschel 500 μm)/f_ν(Spitzer 24 μm) > 6.5 as observed in Rowan-Robinson et al. (2010), although they bracket well the spread in f_ν(Spitzer 70 μm)/f_ν(Spitzer 24 μm). The templates also underrepresent the diversity of slopes in the mid-IR (Wu et al. 2010). It should be emphasized that all template libraries in the literature are based on broadband photometry, mainly using ISO and IRAS filters. However, the detailed SED shape at fine scales within a filter can actually affect model predictions for broadband surveys as the different parts of spectra are redshifted into a fixed wavelength filter.

3.3.2. IR SEDs At High-z

The sensitive infrared surveys carried out by Spitzer and Herschel are revealing the SED shapes of high redshift luminous sources to be significantly different from local galaxies at the same luminosity (e.g., Papovich et al. 2007; Rigby et al. 2008; Symeonidis et al. 2009; Muzzin et al. 2010; Hwang et al. 2010; Rujopakarn et al. 2011a, 2011b; Elbaz et al. 2011). Many high-z LIRGs and ULIRGs show IR SEDs similar to those of local, less luminous galaxies, i.e., high-z LIRGs and ULIRGs are colder compared to their local counterparts; IR luminosities do not correspond to unique properties of star-forming regions. To account for this fact, we evolve the L_{*, C} in Equation (11) through

$$\begin{equation} L_{*, C}(z) = L_{*, C}(z=0)(1+z)^{p_{c}}. \end{equation} \tag{ 12 }$$

It is, however, still difficult to constrain p_c because of the limited number of available broadband photometric bands, the dynamic range of the IR luminosities at a fixed redshift, and the observational selection biases. It is possible that the evolution in the SF SED only takes place at z ≳ 1.5 (Hwang et al. 2010; Elbaz et al. 2010), but it is still not certain whether or not the basic shape of Equation (11), i.e., a two power-law mean trend and a Gaussian distribution of the scatter, is valid at high-z. For example, Elbaz et al. (2011) described the high-z L_TIR/L_8 μm color as a Gaussian distribution plus a higher value tail that is independent of the IR luminosity. In this paper, we assume that the form of the local relationship holds at high-z and explore two possibilities of p_c, zero for no evolution and 2.0 for rapid evolution to bracket the range as discussed in the literature (e.g., Symeonidis et al. 2009; Muzzin et al. 2010; Hwang et al. 2010; Magnelli et al. 2009; Chapin et al. 2011). Note that we ran the model with various p_c values from 0.0 to 3.0 but failed to find that a non-zero p_c value can provide a better fit compared to the no SED evolution case. Because of this, the p_c = 2.0 is used to illustrate the effect of the SED evolution in the model's fitting.

3.3.3. X-Ray SED

Deep X-ray surveys, especially in the soft band, can probe emission from SF regions. We thus related the star-forming X-ray emission to the SF IR luminosity and constructed the X-ray spectrum to go with each IR SED template. In galaxies with active SF, the X-ray emission is dominated by a short-lived high-mass X-ray binary and is known to be related to the star formation rate (SFR; e.g., Grimm et al. 2003; Persic et al. 2004; Lehmer et al. 2010). In quiescent objects, the low-mass X-ray binary emission becomes important and X-ray luminosity is no longer a reliable SFR tracer. Nevertheless, IR emission also has a contribution from dust heated by long-lived stars, so it is reasonable to expect a relation between IR and X-ray even in quiescent objects. We thus compiled measurements from Persic et al. (2004) and Lehmer et al. (2010) but excluded AGN and AGN+SF composite objects. For a few common objects in both studies, data from the latter are used. The total IR luminosities unavailable in the former study are measured by our own based on the IRAS four-band fluxes through the Sanders & Mirabel (1996) equation. The final sample covers the IR luminosity from 10⁸ up to >10¹²L_☉, as shown in Figure 2. A two power-law fit to the data gives

$$\begin{eqnarray} {\rm log}L_{\rm 2\hbox{--}10\,keV}^{\rm SF} = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}-4.70 + 1.08\times {\rm log}L_{\rm TIR}^{\rm SF} & \!\! {\rm for }\; {\rm log}\frac{L_{\rm TIR}^{\rm SF}}{L_{\odot }} \,{\le}\, 10.31\\ \ms -0.59 + 0.68\times {\rm log}L_{\rm TIR}^{\rm SF} & \!\! {\rm for }\; {\rm log}\frac{L_{\rm TIR}^{\rm SF}}{L_{\odot }} \,{>}\, 10.31,\\ \end{array}\right.\nonumber \\ &&\hspace{-10pt} \end{eqnarray} \tag{ 13 }$$

where all luminosities are in units of solar luminosity. The 1σ scatters for trends at the low- and high-luminosity ends are 0.28 and 0.27 dex, respectively. The fit does not use only the one upper limit, but the result does not change if we include it. As discussed by Lehmer et al. (2010), the scatter between X-ray emission and SFR can be reduced if we include the stellar mass quantity, which, however, is not feasible for our model that does not have optical/near-IR components to infer stellar masses. Independent of the X-ray luminosity, the spectrum from 0.2 to 100 keV is fixed to be that of the star-forming galaxy M82 based on the BeppoSAX observation from Cappi et al. (1999) and Swift data from Baumgartner et al. (2012). We further assume no redshift evolution for the L^SF_X-ray − L_IR^SF relationship and X-ray spectrum. The lack of redshift evolution in the L^SF_X-ray − L_IR^SF relationship has been tentatively suggested by some works (Grimm et al. 2003; Symeonidis et al. 2011).

Zoom In Zoom Out Reset image size

The grand unified model of AGNs (Antonucci 1993) states that all AGNs are intrinsically the same and the emerging SED is only determined by the line of sight with respect to the dusty torus; the type 2 objects are the obscured version of the type 1 objects. Adopting this scenario, we first constructed the intrinsic AGN SED from type 1 data set and then modified it with extinction as characterized by the H i column density ($N_{{\rm H\,\scriptsize{I}}}$) to obtain the extinction-dependent SED. A reflected component in the X-ray is always added:

$$\begin{eqnarray} f_{\nu }^{\rm AGN\ observed} &=& f_{\nu }^{\rm AGN\ intrinsic}\times {\rm exp}(-\tau (\nu , N_{{\rm H\,\scriptsize{I}}})) \nonumber \\ && + {\rm reflected-component}, \end{eqnarray} \tag{ 14 }$$

where f^{AGN observed}_ν is the observed SED of AGNs, f^{AGN intrinsic}_ν is the intrinsic SED before extinction, and $\tau (\nu , N_{{\rm H\,\scriptsize{I}}})$ is the extinction curve.

The significant difference from the SF SED templates is that we only used the mean SMBH SED but did not consider its dispersion. This is partly to keep computations manageable but also because of the lack of studies that fully address the dispersion of the AGN SED from the X-ray all the way to the submillimeter as a function of the BH luminosity (Hao et al. 2012).

3.4.1. The IR/Submillimeter SED of Unobscured AGNs

The distinct IR continuum feature of the unobscured AGN is strong emission of dust heated to near the sublimation temperature (∼1000 K) by UV/optical radiation from the accretion disk. This feature serves as the basis of various mid-IR AGN selection techniques (Lacy et al. 2004; Stern et al. 2005; Armus et al. 2007; Donley et al. 2008; Wu et al. 2011). The shape of this feature has been well characterized with the emission peaking at 3–5 μm and the 1–30 μm slope α = 1 ± 0.3 (f_ν ∝ ν^−α), although the slope estimate depends on the definition of the photometric bands and their placement with respect to the silicate feature besides the intrinsic dispersion. Generally, the broadband photometric observations indicate slopes between 1.0 and 1.5 (e.g., Neugebauer et al. 1987; Haas et al. 2003). The Spitzer Infrared Spectrograph (IRS) observations of a large sample of local AGNs reveal smaller slopes for the continuum underlying the silicate feature. For example, the whole sample of 87 PG quasars shows a relatively universal slope with median and standard deviations of α_5–30 μm = 0.75 ± 0.35 (f_ν ∝ ν^−α) (Shi et al. 2007; also see Netzer et al. 2007; Maiolino et al. 2007). The Seyfert 1 objects in the complete 12 μm sample have α_{15–30 μm} = 0.85 ± 0.61 (Wu et al. 2009) where the slightly steeper slope and the larger dispersion are likely caused by the larger SF contamination in these low-luminosity AGNs. No obvious evidence has been found in the above studies for the dependence of the slope on the BH luminosity.

The most prominent spectral feature of the unobscured AGN SED is the silicate emission features at 9.7 and 18 μm as established unambiguously by Spitzer observations (e.g., Hao et al. 2005; Shi et al. 2006). It is important to include this feature in the SED template for our model as it can affect the broadband photometry. The strength of this feature depends weakly on the SMBH luminosity (Maiolino et al. 2007), with only a factor-of-two increase for an SMBH luminosity range of 10⁴. We thus adopt the same feature strength independent of the SMBH luminosity.

The most uncertain part of the SMBH IR SED lies above ∼30 μm where the SF emission becomes important or even dominates for quasars. One solution is to spatially separate the nuclear emission from SF regions in host galaxies; such observations are only feasible below 30 μm (e.g., Gandhi et al. 2009). Another challenging problem with this strategy is that we have no idea if the dusty torus itself is experiencing SF and how important that is. The current approach is to scale the SF with the observed aromatic feature and then subtract it off (e.g., Shi et al. 2007; Netzer et al. 2007). This leads to the conclusion that above ∼100 μm, all the observed emission typically comes from SF, and that the SED part then relies on the extrapolation from the short wavelength. In spite of significant caveats associated with the AGN SED above 30 μm, its effect on our model is limited. This is because the majority (90%) of the AGN IR emission (3–1000 μm) emerges below 30 μm; their far-IR contribution is negligible to the observed far-IR radiation.

In practice, we constructed the unobscured IR SED based on Spitzer observations of the whole PG sample that is selected based on the UV excess and are thus type 1 objects with little obscuration (Green et al. 1986). To construct the intrinsic SED, we first derived the median composite spectrum of Spitzer/IRS spectra (5–30 μm) and MIPS photometry (24, 70, and 160 μm) of all of the PG objects. With the composite spectrum, we measured the 11.3 μm aromatic feature and scaled it to derive the star-forming emission at the Spitzer wavelength based on the DH02 template. The SF emission contributes less than 15% at <30 μm and reaches ∼50% in the 70 μm filter while dominating the 160 μm filter emission. A continuous IR SED template is then constructed based on the SF-subtracted IRS spectrum and photometry in the 70 μm filter where all the aromatic features are also carefully removed. To fill the gap in the SED beyond 30 μm, as shown in Figure 3, the underlying continuum of the silicate feature is fitted with a double power law (purple dotted line) to the 5–8 μm IRS part and MIPS 70 μm photometry where the power-law index at long wavelengths is fixed to be 4 to mimic the modified blackbody with emissivity of 1.0 (νf_ν ∝ν⁴). A single power law connecting the IRS 30 μm emission and MIPS 70 μm photometry (dashed red line) is used to model the silicate emission between 30 and 38 μm. The resulting IR SED is quite consistent with the photometry-based SED derived from the Sloan Digital Sky Survey (SDSS) quasars by Richards et al. (2006), except for the much better characterized silicate features in our template.

Zoom In Zoom Out Reset image size

Extracting the intrinsic IR SED for a large sample of AGNs with a range of host galaxy contamination, extinctions, etc. would be quite difficult. Moreover, the spatial and spectral resolution data are not available to allow for a result comparable to that for star-forming galaxies with well characterized mean SED and associated scatter as a function of the luminosity (Dale & Helou 2002; Chapman et al. 2003). We thus adopt the mean PG quasar IR SED with no luminosity dependence or scatter. The recent work by Mullaney et al. (2011) claimed a luminosity dependence of the AGN intrinsic mid-IR/far-IR ratio in a small AGN sample. However, their low- and high-luminosity sub-samples have completely different extinctions, which could easily induce such a dependence. The studies of X-ray AGNs in the COSMOS field show no apparent luminosity and redshift dependence for their SED (Hao et al. 2012).

3.4.2. X-Ray SED As a Function of $N_{{\rm H\,\scriptsize{I}}}$

For the X-ray SED as a function of the H i column density, we basically adopted the work of Gilli et al. (2007). Since no scatter is included for the IR part in our model, we did not consider the scatter in the X-ray spectral slope either. While Gilli et al. (2007) argued that the scatter in the X-ray spectra is important to constrain the Compton-thick AGN population, this will not be the case for our model due to the fundamental difference in the way that the Compton-thick AGN abundance is constrained. Their CXB model derived the Compton-thick AGN by subtracting the Compton-thin AGN contribution to the CXB spectrum around 30 keV. The extrapolation of the Compton-thin AGN contribution from <10 keV to 30 keV certainly depends critically on the scatter of the power-law slope. Our model constrains the Compton-thick AGN through mid-IR data, which depends on the SMBH IR/X-ray relationship. The additional difference between their work and this study is that we included the reflected component for the quasars (e.g., Piconcelli et al. 2010). A summary is given as follows.

• 1.  
  The unobscured SMBH SED is composed of three components, a primary power law with a cutoff at high (>100 keV) energies, a reflected component and iron 6.4 keV emission line. The power-law component arises directly from the hot corona above the accretion disk. Its shape is characterized by a photon index of Γ = 1.9 (Nandra & Pounds 1994; Reeves & Turner 2000; Piconcelli et al. 2005; Mainieri et al. 2007) and a cutoff energy of 300 keV (Molina et al. 2006; Dadina 2008). For the latter, the uncertainty is still important. This component can be reflected off the far side of the accretion disk or dusty torus, hardening the emerging radiation. The reflected component depends on the incident spectrum slope and the observed inclination. It is modeled for Γ = 1.9 and the relative normalization is fixed to be 1.3, which is the average for type 1 lines of sight assuming a torus half-opening angle of 45°. The iron line is assumed to have a Gaussian profile with a width of 0.4 keV and an EW of 280 eV.
• 2.  
  For the obscured Compton-thin AGN, the intrinsic power law is the same as that for the unobscured one but the reflected component normalization is reduced to 0.88 given different viewing angles toward obscured nuclei. The extinction curve is the photoelectric absorption at solar abundances as characterized by Morrison & McCammon (1983). In total, we adopt $\log({N_{{\rm H\,\scriptsize{I}}}}/{\rm cm^{-2}})$ = [21.5, 22.5, 23.5] to sample the Compton-thin X-ray SEDs.
• 3.  
  For Compton-thick AGNs with $N_{{\rm H\,\scriptsize{I}}}$ < 10²⁵ cm⁻², the intrinsic power law is the same as that for the unobscured one while the normalization of the reflected component is fixed at 0.37. The extinction curve needs to consider the non-relativistic Compton scattering as modeled through the code of Yaqoob (1997).
• 4.  
  For Compton-thick AGNs with $N_{{\rm H\,\scriptsize{I}}}$ > 10²⁵ cm⁻², all transmitted photons are down-scattered by Compton recoil and only the reflected component is visible.

The relative normalization of the X-ray spectrum with respect to the IR part is derived through the modified version of the tight relationship between the 2–10 keV and 12.3 μm luminosities for 22 Seyfert galaxies in Gandhi et al. (2009). These galaxies are well resolved at 12.3 μm so that SF can be removed reliably. Galaxies in that work are a mix of Seyfert 1, Seyfert 2, and Compton-thick AGNs, while what we need here is the relationship only for type 1 AGNs, specifically PG quasars. At 12.3 μm, the silicate emission feature can be important. As shown in Figure 4, we thus modify the Gandhi et al. (2009) relationship by keeping the slope but modifying the normalization to produce the median X-ray/12.3 μm ratio of 34 PG quasars whose X-ray and IR data are derived from Cappi et al. (2006) and Shi et al. (2007), respectively. For these PG objects, the contribution from SF to the 12.3 μm emission is negligible as evidenced by the extreme low aromatic feature equivalent width (EW_{11.3 μm PAH} ≲ 0.05 μm). The derived relationship is

$$\begin{equation} {\rm log} \frac{L_{\rm 2\hbox{--}10\,keV}}{\rm 10^{43}\ {\rm erg\ s^{-1}}} = 0.90\, {\rm log}\frac{L_{12.3\,{\mu }{\rm m}}}{\rm 10^{43}\ {\rm erg\ s^{-1}}} - 0.45, \end{equation} \tag{ 15 }$$

where L_2–10 keV and L_12.3 μm are unobscured 2–10 keV and 12.3 μm luminosities, respectively. To evaluate the effect of this scatter on the model, we explored two additional variants of the model with the normalizations 0.2 dex higher and lower, respectively, where 0.2 dex is the 3σ dispersion of the mean ratio of the PG objects.

Zoom In Zoom Out Reset image size

3.4.3. H i Column Density Distribution

The H i column density distribution is key to understanding AGN physics and the CXB. Its proper characterization requires unbiased sample selections and secure measurements of H i column densities. In the local universe, optically selected samples offer the least bias compared to UV/X-ray (biased against heavily extinguished objects and low-luminosity AGNs) and IR (biased against unobscured objects and low-luminosity AGNs). Maiolino & Rieke (1995) have attempted to compile all Seyfert galaxies at m_B < 13.4. The Seyfert 2 to Seyfert 1 ratio is found to be as high as 4:1. This ratio is consistent with that found in the Palomar nuclear spectroscopic survey of local galaxies at m_B < 12.5 (Ho et al. 1997). Both studies dealt carefully with host galaxy light dilution and the result should be quite reliable. Both studies probed very low luminosity AGNs with bolometric luminosities around 10⁹ L_☉. Hao et al. (2005) instead found 1:1 for the Seyfert 2/Seyfert 1 ratio around L_bol ∼ 10⁹L_☉ and decreasing toward higher luminosities. The discrepancy can be partly due to the luminosity dependence but also that the SDSS spectra in Hao et al. (2005) with large aperture and relatively low quality suffer significant host dilution. By building a Seyfert 2 sample from Maiolino & Rieke (1995) and Ho et al. (1997), Risaliti et al. (1999) measured the H i column density distribution. Combining this with works of Cappi et al. (2006) for Seyfert 1 galaxies, we characterized the column density distribution with three bins:

$$\begin{eqnarray} &&P^{\rm N_{H\,\scriptsize{I}} } = \frac{dN}{d{\rm log}N_{{\rm H\,\scriptsize{I}}}} = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}0.5f_{\rm type\hbox{-}1}\; {\rm for }\; {\rm log}\frac{N_{{\rm H\,\scriptsize{I}}}}{\rm cm^{-2}} \in [0, 22] \\ \ms 0.5(1-f_{\rm type\hbox{-}1})(1-f_{\rm CT})\; {\rm for }\\ \ms \quad \times \log \frac{N_{{\rm H\,\scriptsize{I}}}}{\rm cm^{-2}} \in [22, 24] \\ \ms 0.5(1-f_{\rm type\hbox{-}1})f_{\rm CT}\; {\rm for }\\ \ms \quad \times \log \frac{N_{{\rm H\,\scriptsize{I}}}}{\rm cm^{-2}} \in [24, 26], \\ \end{array}\right. \end{eqnarray} \tag{ 16 }$$

where $P^{\rm N_{H\,\scriptsize{I}} }$ is the probability per logarithm of the H i column density and the unobscured AGN fraction f_type-1 is defined as the number of objects with $\log ({N_{{\rm H\,\scriptsize{I}}}}/{\rm cm^{-2}})$ < 22 to the total, while f_CT is defined to be the fraction of CT ($\log({N_{{\rm H\,\scriptsize{I}}}}/{\rm cm^{-2}})$ > 24) among the obscured AGNs (22 < $\log({N_{{\rm H\,\scriptsize{I}}}}/{\rm cm^{-2}})$ < 26).

Beyond the local universe, it is increasingly difficult to measure the H i column density distribution directly. Our model thus allows both luminosity and redshift evolution of f_type-1 and f_CT:

$$\begin{eqnarray} f_{\rm type\hbox{-}1} &=& f_{0, \rm type\hbox{-}1}\times \frac{1}{(1+z)^{\beta _{z, {\rm type\hbox{-}1}}}}\nonumber \\ && \times \exp \left(-\left({\rm log}L_{\rm TIR}^{\rm BH}-9.0\right)\times \beta _{l, {\rm type\hbox{-}1}}\right) \end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray} &&f_{\rm CT} = f_{0, \rm CT}\times \frac{1}{(1+z)^{\beta _{z, {\rm CT}}}}\times {\rm exp}\left(-\left({\rm log}L_{\rm TIR}^{\rm BH}-9.0\right)\times \beta _{l, {\rm CT}}\right), \nonumber \\ \end{eqnarray} \tag{ 18 }$$

where f_{0, type-1} is fixed to be 0.2 and f_{0, CT} is fixed to be 0.5 at L^BH_TIR = 10⁹ L_☉, in order to approximate the result of Maiolino & Rieke (1995), Ho et al. (1997), and Risaliti et al. (1999). We also impose a maximum and minimum fractions of 95% and 5%, respectively, for each column range (type 1, Compton-thin type 2, and Compton-thick). There are four free parameters in total to describe the H i column density distribution.

3.4.4. The UV/Optical Part of Unobscured SED

The composite UV/optical/near-IR SED is taken from Richards et al. (2006) based on the SDSS type 1 quasars. We combined this part with X-ray emission through the X-ray/optical luminosity relationship as established by Steffen et al. (2006)

$$\begin{equation} \alpha _{\rm ox}=-0.107\log \left(\frac{L_{\nu ,2500\ {\mathring{\rm{A}}}}}{\rm erg\ s^{-1}\ Hz^{-1}}\right)+1.739, \end{equation} \tag{ 19 }$$

where α_ox = −0.384log(L_{ν, 2500 Å}/L_{ν, 2 keV}).

3.4.5. The Extinction Curve at UV/Optical/IR Wavelength

Our model first assumes total IR LFs that are described by a double power law with four parameters (break luminosity and density along with the faint- and bright-end slopes). The IR LF is assumed to evolve with redshift in both break luminosity and density, each of which is a polynomial function of redshift with three parameters. Therefore, the total IR LF has a total of 10 parameters. The model then introduces the SMBH energy fraction in the total IR band to decompose the emission into SMBH and star-forming components. At a given IR luminosity, the distribution of SMBH energy fractions is assumed to be a Gaussian function with both mean and scatter depending on the luminosity, while the mean is further assumed to show redshift dependence. There are five free parameters in total, including two for the redshift dependence of the mean and one for its luminosity dependence, one for the scatter at the total IR luminosity of 10^12.3 L_☉, and one for its luminosity dependence. The mean fraction at 10^12.3 L_☉ is fixed based on local ULIRG objects. For each power source (star forming and SMBH), the model then converts the emission in the total IR band to other wavelengths using star-forming and AGN SEDs. For the former, luminosity-dependent empirical SEDs are used. For the latter, H i columns are added to extinguish the type 1 empirical SEDs. H i columns are divided into three ranges representing type 1, Compton-thin type 2, and Compton-thick, which are described by two parameters to give the relative frequencies of occurrence (the sum of which needs to be normalized). Each parameter has a local value fixed at the SMBH IR luminosity of 10⁹ L_☉, and is assumed to follow a simple power-law dependence on both luminosity and redshift. Therefore the H i columns have four parameters in total. The model has a total of 19 free parameters.

5.1.1. X-Ray Number Counts

5.1.2. Spitzer Number Counts

5.1.3. IRAS Number Counts

5.1.4. Herschel Number Counts

5.1.5. SCUBA/Aztec/Mambo Submillimeter Number Counts

5.2.1. X-Ray Redshift Distributions

5.2.2. IR/Submillimeter Redshift Distributions

We fitted four local LFs, including the local 15–55 keV LF (Ajello et al. 2012), the local 2–10 keV LF (Ueda et al. 2011), the IRAS 12 μm LF (Rush et al. 1993), and the IRAS 25 μm LF (Shupe et al. 1998).

For the reference model (the solid line), the differences between the data and the fit include 0.1 dex constant overfits of the 15–55 keV number counts, 0.2–0.3 dex underfits of the 2–10 keV counts at the bright end, the 0.2 dex underfits around 0.5–2 keV fluxes from 1.5 to 4.0 × 10⁻¹⁷ erg s⁻¹ cm⁻², 0.1–0.2 dex underfits around the peak 24 μm count in the flux density range of 0.2–0.5 mJy, 0.2 dex underfits and overfits at the 70 μm bright and faint ends, respectively. Unlike the 15–55 keV counts, the 17–60 keV ones are well produced by our model so that the discrepancy seen in 15–55 keV can be due to subtle differences in the X-ray SED around the range of 15–17 keV and 55–60 keV. It is yet interesting to note that, similar to ours, several other models (Gilli et al. 2007; Treister et al. 2009; Draper & Ballantyne 2010) overproduce the 15–55 keV counts by 0.1–0.2 dex, while their adopted SEDs are not exactly the same. The difference seen at the bright end of the 2–10 keV counts is likely caused by contamination of galaxy clusters in the data (Gilli et al. 1999). In 0.5–2 keV counts, the data clearly show a transition flux (3 × 10⁻¹⁷ erg s⁻¹ cm⁻²) below which normal galaxies other than AGNs start to dominate, causing the counts to go above the simple power-law extrapolation from the higher fluxes (Bauer et al. 2004; Lehmer et al. 2012), while the model also produces a similar feature but at fluxes ∼3 times lower. For the 24 μm peak counts, although the highest observed point is the result of the SWIRE data near the detection limit (Béthermin et al. 2010), the differences are seen over several adjacent data points. Since this peak is caused by aromatic features moving into the 24 μm filter around z = 1, a possible explanation is that z  ∼  1 star-forming galaxies have slightly increased aromatic-feature EWs. Differences seen in 70 μm are also noticeable, for which the subtle caveat in star-forming SEDs can also be the potential cause.

For the fast_evol_SED_SF variant (dashed line), additional differences between the fits and data, besides those for the reference model, are the significant underestimates of the 70, 100, and 160 μm counts above the flux density of the peak count (Figures 6(h)–(j)). The reason is that the redshift evolution in the SF SED results in a colder SF color f_60 μm/f_100 μm that decreases the rest-frame 30–100 μm emission. The behavior of the other two variants (low_IR2X_BH and high_IR2X_BH) in fitting the data are quite similar to that of the reference variant.

The redshift distributions are shown in Figure 7. The distributions at 2–10 keV are plotted in Figures 7(a)–(g). The data are from five fields including CDF-N, CDF-S, ECDF-S, COSMOS, and XMS that complement each other by probing from small deep field (0.12 deg²) to wide shallow area (3 deg²). The limiting fluxes range from >3 × 10⁻¹⁶ up to > 5 × 10⁻¹⁴ erg s⁻¹ cm⁻². Compared to the number counts, the observational uncertainties here are larger due to the redshift measurement errors or cosmic variances. For example, for the CDFN and CDFS (Figures 7(a) and (b), respectively) with the same field size and limiting X-ray flux, we see obvious differences in the data. In the soft X-ray (Figures 7(h)–(n)), the data extend from >5 × 10⁻¹⁷ to >7 × 10⁻¹⁴ erg s⁻¹ cm⁻² while the field area ranges from 0.12 to 80 deg². Again, cosmic variances and photo-z errors play roles in the differences among redshift distributions with similar flux cuts, emphasizing the limitations of the current available data. The redshift distributions in the IR bands are shown in Figures 7(o)–(u). Although the redshift measurements in general have high completeness, the incompleteness is always redshift-dependent, missing more objects at higher z. This seems to be consistent with the fact that our model produces more objects than the observations at the high-z end in several panels (e.g., Figures 7(m), (n), and (r)).

All variants of the model perform in similar ways. The fit to the 2–10 keV distributions is generally acceptable where the distribution with the brightest cut as shown in Figure 7(g) shows the largest difference between the fits and the data. In the 0.5–2 keV, the two distributions with the brightest flux cut (Figures 7(m) and (n)) have the worst fits. Models overproduce z > 2 sources, potentially caused by the fact that sources without redshift measurements lie at high-z. At 24 μm, for the two lowest flux cuts (Figures 7(o) and (p)), all variants of the model well reproduce two peaks around z = 1 and z = 2 when different aromatic features are redshifted into the 24 μm filter. The only difference is seen at z < 0.5, where the COSMOS field does not probe a significant cosmic volume. For the GOODS-N (Figure 7(q)), all the fits are systematically lower at z < 0.7. Here the cosmic variance or photo-z error should also be important; the GOODS-N distribution has a higher flux cut but shows higher densities at z < 0.5 compared to two COSMOS distributions. The match is excellent to the highest flux cut at z < 2 (5 mJy in 24 μm, Figure 7(r)) where the data quality is high (effectively 10 deg² field size and 98% spectroscopic completeness; Wu et al. 2010). For this distribution, sources with no redshift measurements are not included but likely lie at high-z according to their photo-z, which at least partly causes the difference between the model and the data at z >2. All the fits to redshifts of local IRAS 60 μm sources brighter than 5.24 Jy are excellent. All the fits to the SCUBA 850 μm and Aztec 1100 μm are acceptable but produce higher mean redshifts.

For the 15–55 keV LF (Figure 8(a)), all variants of the model produce 0.3–0.5 dex higher spatial densities than observations below the break luminosity. The fits to the 2–10 keV LF (Figure 8(b)) offer acceptable fits given the error bars. For the IRAS 12 μm as shown in Figure 8(c), all variants of the model produce slightly flatter bright-end slopes. At the faint end, all variants underestimate the observed LF. The local supercluster is known to boost the spatial density of faint sources, an effort that is not yet corrected in the 12 μm LF (Rush et al. 1993; Fang et al. 1998; Shupe et al. 1998). For the IRAS 25 μm LFs as shown in Figure 8(d), the fits are excellent.

All variants of the model, including the reference, fast_evol_SED_SF, low_IR2X_BH, and high_IR2X_BH as listed in Table 1, offer similar performance in fitting number counts, redshift distributions, and local LFs, with the fast_evol_SED_SF variant showing additional significant underestimates of the 70, 100, and 160 μm number counts above the flux density of the peak count. The reduced χ² ranges from 2.7 to 2.9. The contribution to the reduced χ² from systematics associated with the data is estimated to be at least 0.75 as shown in the following. For the number counts, we assumed negligible data systematic errors so that the difference between model predictions and data are due to the limitation of the model. For the local 12 μm LF (Figure 8(a)), the excess faint sources due to the local super cluster are not corrected, causing higher densities than our model's predictions. We quantitatively estimate its contribution to χ² by summing the square of the difference between observations and predictions divided by the observed errors at the 12 μm luminosity below 10⁹ L_☉. Cosmic variances or photo-z errors are evident for redshift distributions between different fields with similar flux cuts including at 2–10 keV, CDFN in Figure 7(a) versus CDFS in Figure 7(b), ECDFS in Figure 7(c) versus CDFN in Figure 7(d), COSMOS in Figure 7(e) versus CDFS in Figure 7(f); at 0.5–2 keV, CDFS in Figure 7(h) versus CDFN in Figure 7(i), and CDFN in Figure 7(k) versus COSMOS in Figure 7(l). In addition, they are evident in the 24 μm bands at z  <  0.3, COSMOS above 0.08 mJy in Figure 7(o) versus GOODS-N above 0.3 mJy in Figure 7(q). Although this pair has quite different flux cuts, the cosmic variance or photo-z errors must play dominant roles because GOODS-N with much higher flux cuts shows significantly higher source densities at z < 0.3 compared to the COSMOS field. Since the best fit to two Y values at fixed X is the mean of two Y values, independent of any model assumption, we estimated the minimum contribution to χ² for each pair as listed above by assuming the best fit to be the mean of two distributions in the pair. The sum of all the above contributions normalized by the degree of freedom gives 0.75, indicating that the true reduced χ² that our four variants of the model would achieve if data systematics was removed could be better than 2.0.

The comparison of the model variants and their performance in fitting the data does not favor strong SF SED evolution; on the other hand, direct comparisons with the observed high-z SED color do not provide definitive evidence about SED evolution. Figure 9 shows such comparisons where we sampled the model output to mimic the observational selection, which is crucial to eliminate selection biases (Magnelli et al. 2009; Chapin et al. 2011). However, it is impossible to ascertain that all selection biases have been properly represented. Figure 9(a) verifies that the model reproduces the local relationship (Equation (11)). The difference in the median color at a given IR luminosity is around 0.02, which is reasonable given the limited resolution (0.037) of the numerical calculation. Beyond z = 0, we tested the evolution of the SED shape by comparing to five samples, an SMG sample from Chapman et al. (2005), two star-forming samples from Magnelli et al. (2011), the Herschel GOODS sample from Elbaz et al. (2011), and a star-forming sample from Kirkpatrick et al. (2012). For the SMG sample, the radio emission in our model is derived through the radio–FIR relation (for details, see Dale & Helou 2002). As shown in Figure 9(b), although the scenario of no evolution in the SF SED offers a slightly better match to the 850 μm/1.4 GHz ratio of the SMG sample, these SMGs are radio selected for follow-up observations and thus subject to possible contamination from SMBH radio emission. In Figure 9(c), the tight empirical relationship between L_{observed 24 μm} and L_{observed 70 μm} of z = 1.3–1.8 star-forming galaxies is defined by stacked data and it favors strong evolution in the SF SED. For the stacked color of z =1.8–2.3 galaxies, the variants of no SED evolution perform better at the low-luminosity end while the variant with strong SED evolution does better at the high-luminosity end. The well-defined relationship between rest-frame 8 μm and total IR luminosities of the Herschel GOODS sample in Figure 9(e) is based on individual detections, but on those from different redshifts. Only at the high-luminosity end, corresponding to z ∼2, does the variant with strong SED evolution produce results closer to the observations. In Figure 9(f), the L_{obs 24 μm} versus L_{obs 250 μm} relationship is based on 24 μm sources (>0.1 mJy) with Spitzer mid-IR spectroscopic measurements at the redshift range of 0.3–2.5. We selected objects with inferred AGN fractions smaller than 20% to define the above relationship. The 250 μm detection rate is quite high (95%) so that the relationship should not bias against low 250 μm/24 μm sources. The model variants with no SED evolution better match the median of the distribution compared to the variant with SED evolution. At the bright 24 μm luminosity end that corresponds to z = 2, modeled galaxies with SED evolution are too cold to produce enough observed frame 250 μm (i.e., rest-frame 80 μm) radiation. We therefore conclude that although the high-z ULIRGs are colder than local counterparts as argued in the literature, their IR SEDs do not represent exactly those of local normal star-forming galaxies across the whole IR band (e.g., Kirkpatrick et al. 2012).

Zoom In Zoom Out Reset image size

As listed in Table 3 and shown in Figures 12 and 13, all variants of the model indicate that the median SMBH energy fraction in the total IR band increases with the IR luminosity, which is consistent with various studies in the local universe (e.g., Sanders & Mirabel 1996; Genzel et al. 1998; Armus et al. 2007; Desai et al. 2007; Veilleux et al. 2009; Petric et al. 2011) and at high-z (Kartaltepe et al. 2010), while the Gaussian width of the distribution increases with the decreasing luminosity. The overall trend of the median SMBH fraction with redshift shows a minimum around z ∼1–2. To test the robustness of this redshift trend, we first ran a model without redshift evolution (k^BH_{1, d}  ≡  0, k^BH_{2, d}  ≡  0) while all others are the same as those in the reference variant. The resulting best-fit reduced χ² increased by as much as 1.3, which offers strong evidence for the requirement of redshift evolution in the SMBH energy fraction. We then ran a variant of the model that only allows monotonic evolution (i.e., k^BH_{1, d} = free, k^BH_{2, d} ≡ 0). The increase in the best-fit reduced χ² is smaller but still significant by 0.4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 14, we compared the 6.2 μm EWs in our model to several observed distributions at different redshifts, which offers a stringent check for our SMBH energy fraction measurements. Here the modeled samples are defined in similar ways to those observed. All EWs are measured using the spline interpolations, with the pure SF median value around 0.6–0.7 μm. Note that our SF templates do not have scatters in the EW, and thus the modeled EWs do not have values in bins higher than EW = 0.6–0.7 μm and any value below EW = 0.6 μm is caused by SMBH contributions in our model. At logL_TIR = 12.3 and z = 0, the SMBH energy fraction in the model is fixed to that of 1 Jy ULIRGs, and is allowed to be luminosity and redshift-dependent. As shown in Figure 14(a) for ULIRGs, the model produces an almost flat trend, roughly consistent with the result of GOALS's ULIRGs (S. Stierwalt et al., in preparation) that additionally show excess objects in the lowest EW bin above a flat trend. For local LIRGs+ULIRGs, the distribution has a peak at the SF EW with a tail toward the low EW end, which is roughly consistent with the GOALS's result (S. Stierwalt et al., in preparation). Beyond the local universe, we compared the model's results to those of 5MUSES (Wu et al. 2010). At 0.05  <  z  <  0.3, the higher SF peak in the model compared to the observation is due to zero scatter in our SED templates. Overall, the SF peak plus a low EW tail as produced by the model is consistent with the observations. At 0.3  <  z  <  0.5, the model reproduces a bimodal EW distribution as observed, an SF peak plus a peak at the lowest bin with the latter higher than the former. For two higher redshift bins (0.5  <  z  <  1.0 and 1.0  <  z  <  2.0), the model reproduces a peak at the lowest bin with a tail toward the SF end as observed. In Figure 14, we also list the fraction of objects with EW < 0.4 for each observed or modeled distribution. This quantitative comparison further supports the consistency between the model's predictions and observations, with the former only slightly lower than the latter.

Zoom In Zoom Out Reset image size

The model characterizes the H i column density distribution through three bins, the type 1 (log$N_{{\rm H\,\scriptsize{I}}}/$cm² < 22), Compton-thin type 2 (22 < log$N_{{\rm H\,\scriptsize{I}}}/$cm² < 24), and Compton-thick AGNs (24 < log$N_{{\rm H\,\scriptsize{I}}}/$cm² < 26). Figure 15 shows the evolution of the type 1 and Compton-thick AGN fractions in all AGNs as a function of the SMBH luminosity and redshift. All four variants of the model require that the type 1 fraction increases with the BH luminosity but decreases with redshift, and the CT fraction has the opposite trend.

Zoom In Zoom Out Reset image size

The trend of the type 1 fraction with the SMBH luminosity is in general consistent with previous studies (e.g., Ueda et al. 2003; Hasinger 2004, 2008; La Franca et al. 2005; Barger et al. 2005; Treister et al. 2006; Ballantyne et al. 2006; Gilli et al. 2007). Although previous CXB models can explain the X-ray survey data either with the redshift dependence of the type 1 fraction (La Franca et al. 2005; Ballantyne et al. 2006) or without it (Ueda et al. 2003; Treister & Urry 2005; Gilli et al. 2007), the requirement of such a redshift evolution is highly preferred by our model, which is constructed to explain both the X-ray and IR data, consistent with recent studies of a large X-ray AGN sample (Hasinger 2008). Without redshift evolution (β_{z, type − 1} ≡ 0 & β_{z, CT} ≡ 0), the best-fit reduced χ² rises by 0.5, which is essentially caused by the increased brightness of high-z AGNs at both X-ray and 24 μm wavelengths. In the pure CXB model that only deals with the X-ray data, the increase in the apparent brightness of individual objects can be complemented by the decrease in the spatial number density in order to keep the number of objects above a given flux density unchanged. This is why the CXB model is not sensitive to redshift evolution of the type 1 fraction. However, the change of the AGN number density in our model is additionally constrained by the IR data. The redshift evolution of the break luminosity and density of the total IR LF is largely fixed throughout IR survey data. Although the change in the SMBH energy fraction can vary the AGN number densities, it also modifies the number density of the SF-dominated sources, especially the far-IR/submillimeter sources that lie at the bright end of the IR LF and thus are sensitive to the SMBH energy fraction. Without redshift evolution of the H i column density, the best-fit result shows that the redshift distributions at z > 1 of 24 μm sources brighter 5 mJy and 0.5–2 keV sources brighter than 10⁻¹⁴ erg s⁻¹ cm⁻² are overproduced by a factor of 10, and the redshift distribution of submillimeter sources systematically offsets by a factor of three along the Y-axis.

Our model requires a large number of the Compton-thick AGNs. At <10 keV, Compton-thick AGNs are extremely faint. The traditional CXB model first constrains the Compton-thin AGN through X-ray data at <10 keV and then add the Compton-thick AGN to match the CXB spectrum above 10 keV, especially the 20–30 keV peak (e.g., Ueda et al. 2003; Gilli et al. 2007). In this paper, we take advantage of the capability of the mid-IR data in probing these deeply buried AGNs thanks to significantly lower apparent gas-to-dust ratio. A test run of the model with zero Compton-thick AGN abundance shows that the reduced χ² would worsen by as much as 1.0. The main reason is that the relative fraction between Compton-thick and Compton-thin type 2 AGNs as a function of the flux density is different at 24 μm from that at 0.5–2 and 2–10 keV, so that it is able to break the degeneracy between Compton-thick and Compton-thin type 2 AGNs seen in pure CXB models. As shown in Figure 16, the Compton-thick AGNs are only important over the faintest one to two orders of magnitude flux ranges in X-ray bands, while Compton-thin type 2 AGNs are important over a much larger flux ranges. The resulting Compton-thick to Compton-thin type 2 AGN fraction increases with the decreasing flux density at both bands. Such similar behaviors at both soft and hard bands imply that if the model is forced to replace Compton-thick objects with Compton-thin type 2 AGNs in one band, the result in the other band would not significantly violate the observation. However, at 24 μm, the Compton-thick contribution is noticeable over a flux range as large as four orders of magnitude, which is quite similar to that of the Compton-thin type 2 AGN. The relative ratio between the two populations is not a strong function of the flux density, with a value roughly around 1:3. This means that if the model is forced to replace Compton-thick AGNs with Compton-thin AGNs at 24 μm, the result in X-ray bands would significantly overproduce objects at the high flux end but underproduce them at the faint flux end. This argument was verified by running a test model without a Compton-thick AGN population and the expected behavior was observed in the model output.

Zoom In Zoom Out Reset image size

As shown in Figure 17, the predicted total IR LFs are consistent with the observations over redshift and luminosity ranges where the data are available. For comparison, we converted the observed LFs at given filters to the total IR LFs using the reference variant. We first modeled the number density Φ_model(filter, z, νL_ν) at the same filter, redshift bin (z  ±  Δz), and luminosity bin (νL_ν  ±  ΔνL_ν) as the observed ones (Φ_obs(filter, z, νL_ν)). The νL_ν is then corrected to the total IR luminosity with luminosity corrections based on the probability-weighted average of all model SEDs within z  ±  Δz and νL_ν  ±  ΔνL_ν. This Φ_model/Φ_obs is taken as the difference of the observed and modeled LFs in the total IR band. For the local LFs, the IRAS data offer the largest dynamical range and are thus the best data set to test the validity of the model's LF. The 12 and 25 μm LFs are used to fit the model and the result is discussed in Section 6.3. Here we further compared our model to the total IR LF based on IRAS 60 μm selected sources (Sanders et al. 2003), and found a good match between the two (Figure 17). Consistency is excellent with other works of significantly smaller fields that cover the whole IR/submillimeter wavelength range. Above z > 0.1 up to z ∼ 3, current deep surveys are able to probe objects up to one order of magnitude below the break luminosity as shown in Figure 17. At short wavelengths (Wu et al. 2011; Rodighiero et al. 2010; Rujopakarn et al. 2010), no significant difference is found between the model and data up to z  ∼  2.5. Some deviations at the faint end are most likely caused by incompleteness corrections particularly associated with 1/V_max method (Page & Carrera 2000). At rest-frame 250 μm, the model is consistent with that of Eales et al. (2010) but the result of Dye et al. (2010) is systematically higher at 0.1  <  z  <  0. The total IR LF of Caputi et al. (2007) is systematically above our modeled LF; their work may overestimate the total IR luminosity due to K-corrections (also see Magnelli et al. 2011).

Zoom In Zoom Out Reset image size

Our model hypothesizes no evolution in either the faint- or bright-end slopes of the total IR LF. As shown in Figure 17, the fixed bright-end slope matches the available data up to z  ∼  3. For the faint-end slope, we did not see deviations from the data up to z  ∼  1 above that which the available data are quite limited to justify. With fixed faint- and bright-end slopes, the total IR LF as predicted by our model requires evolution in both luminosity and density. Although the detail relies on the evolution formula adopted, we found in general that the break luminosity (logL_*) increases with redshift but the break density (logΦ_*) decreases with redshift. As shown in Figure 18, for the evolution formula given in Equations (2) and (3), the break luminosity increases rapidly until z = 2 and then gradually reaches the peak at z = 4, while the break density shows a gradual increase until z = 0.4 and then monotonically decreases with redshift. To test the redshift out to which luminosity or density evolution is required, we have carried out several test runs of the model by modifying the evolution formula to force no evolution above redshifts of z = 0.5, 1.0, 1.5, 2.0, 2.5, and 3.0. If we adopt Δχ²/dof = 1 as an indicator of significantly worse fit, we found that the no-evolution scenario can be ruled out below z ∼ 2.5. Above z  ∼  3, the fit result could not differentiate evolution versus no evolution, which is not surprising since the current constraint on the z > 3 infrared and X-ray populations is still of low statistical significance as, shown in Figure 7. To evaluate the degeneracy between the luminosity and density evolution, we carried out test runs that halted density or luminosity evolution at selected redshifts. The result shows that the data are not able to break the degeneracy at z  ≳  1.5, which is consistent with Figure 17 showing that above z ∼ 1.5 the observed LF does not extend significantly below the break luminosity. However, no matter which evolution formula is adopted, the robust result is that the break luminosity always increases with redshift, while the break density overall decreases with redshift. Our result of the positive luminosity evolution is consistent with different literature studies while the negative trend in the density is not always favored. We note that our result of negative density evolution is inconsistent with studies in Le Floc'h et al. (2005), Pérez-González et al. (2005), and Rodighiero et al. (2010), but consistent with ones in Caputi et al. (2007) and Magnelli et al. (2011).

Zoom In Zoom Out Reset image size