DUSTY TORI OF LUMINOUS TYPE 1 QUASARS AT z ∼ 2

In the unified model of active galactic nuclei (AGNs; Antonucci 1993; Urry & Padovani 1995), the dust torus is a region immediately outside the accretion disk where dusty clouds are no longer sublimated by the radiation from the central engine. The dust torus reprocesses the incident ultraviolet (UV)/optical radiation from the accretion disk and this energy emerges in the near- and mid-infrared (IR) bands. Richards et al. (2006) presented panchromatic spectral energy distributions (SEDs) for 259 type 1 quasars selected from the Sloan Digital Sky Survey⁶ (SDSS; York et al. 2000). These quasar SEDs constructed from broadband photometry are remarkably similar over a large range in both luminosity and redshift. However, Richards et al. (2006) noted small differences in the 1.3–8 μm range between optically luminous and optically dim quasars. Gallagher et al. (2007) investigated this further and found that the 1–8 μm spectral index (α_ν) is strongly anti-correlated with IR luminosity in type 1 quasars. More luminous quasars have bluer 1–8 μm slopes. Further, they noted a tight linear correlation between the UV continuum luminosity and the IR luminosity for these quasars. This suggested that the observed near-IR emission at 3 μm in the SEDs of many type 1 objects is driven by the dust reprocessing of the intrinsic optical/UV continuum from the accretion disk, as had been noted previously (Rees et al. 1969; Neugebauer et al. 1979; Edelson & Malkan 1986; Barvainis 1987; Sanders et al. 1989). As a recent example, the near-IR emission is clearly visible in the spectrum of Mrk 1239 (Rodríguez-Ardila & Mazzalay 2006).

Theoretical work on the response of accretion disks to radiation and hydromagnetic pressure suggests that outflow of matter is associated with all accretion disks in the form of a wind coming off the surface of the disk (Konigl & Kartje 1994; Murray & Chiang 1995; Proga et al. 2000). The dusty torus itself may be the outermost part of this accretion disk wind close to the equator of the system (Konigl & Kartje 1994; Elitzur & Shlosman 2006). Disk winds have a natural dependence on luminosity through radiation pressure, and this begs the question: "Is the structure of the dusty torus related to the physics of the accretion disk?" The need for proper radiation transfer treatment of clumps in dusty tori was recognized in pioneering early studies (Krolik & Begelman 1988; Pier & Krolik 1992; Rowan-Robinson 1995), and was fully developed by Nenkova et al. (2002). More recently, Nenkova et al. (2008a) presented their model in detail (denoted by Clumpy hereafter).

Significant effort has been invested in understanding the torus dust distributions with various groups favoring both clumpy and smooth dust density distributions (Nenkova et al. 2002, 2008a; Dullemond & van Bemmel 2005; Schartmann et al. 2005, 2008; Fritz et al. 2006; Hönig et al. 2006). The primary difference between clumpy and smooth models is that of the dust temperature distributions (see Figure 3 of Schartmann et al. 2008). While in smooth density models, the temperature steadily declines with radius from the inner wall, clumpy models can show a range of temperatures at different distances from the central source. This effect occurs primarily due to the shadowing effect from the finite size of clouds. The inner faces of clouds are directly exposed to radiation from the central source, and are hence hotter, while their outer faces are much cooler. And because of clumpiness, clouds farther out in radius can still have their inner faces exposed directly to radiation from the central source.

The effective optical depth in a clumpy torus is a function of the number density of clouds in the central regions of the torus. This important model construction has resulted in better fits to both low-resolution Spitzer spectra (Mor et al. 2009; Nikutta et al. 2009), and high-resolution interferometric observations of dusty tori in NGC 1068 (Jaffe et al. 2004) and Circinus (Tristram et al. 2007).

Clumpy models appear to be the most promising set of models with a wide range of applications to both AGNs (e.g., Mason et al. 2006) and merger-driven ultraluminous infrared galaxies (Levenson et al. 2007). Other notable models that employ clumps arranged in a disk-like geometry include Schartmann et al. (2008) and Hönig et al. (2006). For example, Polletta et al. (2008) employed clumpy torus models from Hönig et al. (2006) to fit their optically obscured but IR-bright sources at high redshift.

Clumpy models show changes in their near-IR continua based on the average number of clouds (N₀) encountered along a radial equatorial ray (see Figure 6 in Nenkova et al. 2008b). Using Spitzer mid-IR spectroscopy of high-redshift quasars it is then possible to constrain the parameters of their dusty tori. While Spitzer archives are rich in observations of low-redshift Seyfert galaxies, they are deficient in high-redshift observations of radio-quiet quasars at the peak of the quasar activity in the universe. In this paper, we present such observations as obtained with the Infrared Spectrograph (IRS) on board Spitzer. Our goals include (1) presenting high-quality mid-IR quasar spectra covering rest-frame 2–12 μm for comparison to the low-redshift templates already available (e.g., Hao et al. 2005; Weedman et al. 2005; Buchanan et al. 2006; Glikman et al. 2006; Shi et al. 2006; Schweitzer et al. 2006); and (2) testing the validity of Clumpy torus models by fitting the observed spectra with model SEDs. Using good-quality IRS spectra we hope to model the 10 μm region properly and constrain Clumpy torus parameters for these luminous quasars.

The properties of the sample and reduction process of the IRS spectra are presented in Section 2. The IRS spectra and SEDs of the sample are discussed in Section 3. Clumpy torus models are summarized in Section 4, and Section 5 presents results of model fits to UV-to-mid-IR SEDs. Results are summarized in Section 7. In all calculations, we assume a standard cosmology with H₀ = 71 km s^-1 Mpc⁻¹, Ω_M = 0.27, and Ω_Vac = 0.730.

2.1. The Sample

2.2. Data Reduction

We obtained the basic calibrated data (BCD) products processed with the standard Spitzer IRS pipeline (version S18.7.0) from the Spitzer Science Center (SSC) archive. We cleaned the BCD images using the IRSCLEAN software package to fix rogue pixels using SSC supplied masks, and a weak thresholding of the pixel histogram. We co-added the multiple data collection event image files into one image for each module, spectral order, and the "nod" position (e.g., SL, first order, first nod position) using the fair-coadd option in SMART. We differenced the co-added images from the opposite "nod" positions to remove the sky background. The spectra were extracted using the optimal extraction option within the SMART package. All the image combining and spectrum extraction operations were carried out using SMART (Higdon et al. 2004). We also checked our extractions using the SPICE program. We obtained an average S/N of ∼6–10 for four archival spectra from program 3046, ∼10–13 for spectra from program 50087, and ∼25–35 for the spectra from proposal 50328. These S/N estimates were commensurate with the pre-determined configuration of each observation. Figure 1 displays the observed spectra plotted along with SEDs.

Zoom In Zoom Out Reset image size

In general, the spectra show two features peaking at ∼3 and ∼10 μm (see Figure 1, features are marked by vertical dashed lines) in νL_ν units. The IR spectral index (α_ν) from 3–8 μm ranges from −0.49 to −1.82, with a median of −0.86. The 10 μm emission feature is the well-known 10 μm silicate emission feature due to the Si–O stretching mode of the silicate molecule. This emission feature was well known in stellar spectra for a long time (e.g., Little-Marenin & Little 1988), but has only recently been detected in quasar spectra (Siebenmorgen et al. 2005; Hao et al. 2005; Sturm et al. 2005) due to the sensitive spectroscopy and broad wavelength coverage possible with Spitzer.

The weakness of the 10 μm emission feature in IR spectra of local type 1 AGNs had motivated suggestions of the presence of different chemical compositions and/or size distributions of dust grains (Laor & Draine 1993; Maiolino et al. 2001). Instead, Clumpy models of Nenkova et al. (2002) make use of the clumpy nature of the dusty medium to improve model fits to the 10 μm region. However, as we will see below, different sublimation temperatures and radii for graphite and silicate grains remain an important issue to be resolved in torus models.

The emission peaking between 2 and 4 μm can be attributed to the blackbody emission from dust close to its sublimation temperature (Rees et al. 1969; Davidson & Netzer 1979; Barvainis 1987), which is typically expected to be T ≳ 1500 K for graphite dust. This hot dust emission has long been expected based on broadband IR data (Sanders et al. 1989). Measurement of the strength of this feature relative to longer wavelength mid-IR emission is important because it can give constraints on the inclination of the torus assuming a disk-like configuration (Pier & Krolik 1993; Murayama et al. 2000). Recent advances in near-IR ground-based spectroscopy have led to observations of the near-IR bump in Mrk 1239 (Rodríguez-Ardila & Mazzalay 2006) and NGC 4151 (Riffel et al. 2009).

Figure 1 shows the SEDs constructed using the photometric data points from Tables 2, 3, and 4. Also overplotted for each object is the IRS spectrum along with the mean quasar SED template from Richards et al. (2006) scaled to the SDSS i-band luminosity for each object. While the mean SED captures the overall trend quite well, individual spectra reveal significant differences from the mean SED. Objects with similar UV luminosities can have different relative IR powers (∼0.3 dex). Obscured sources (e.g., SDSS J142730, program 50087) are significantly more IR luminous than sources with similar observed UV luminosities (e.g., SDSS J172522, program 50087) that are probably not as strongly obscured based on their optical SDSS spectra. This trend is reflected in the mean SEDs constructed by Richards et al. (2006).

We use the Clumpy torus models from Nenkova et al. (2008a) to fit the complete SEDs. The models are constructed by assuming an intrinsic AGN SED that heats the dust clouds (see Figure 4 of Nenkova et al. 2008a). We do not consider the effects of a different intrinsic AGN SED. This effect was partially studied by Nenkova et al. (2008a; see their Figure 12), and it is expected that the SED longward of 1 μm should not change significantly. Clumpy models contain a standard Galactic mix of silicates (53%) and graphite (47%) dust grains. We have not explored changes in composition and size distribution of dust grains (see, e.g., Laor & Draine 1993), and contributions from species other than silicates (e.g., Markwick-Kemper et al. 2007). These areas should be addressed by future work on torus models.

The Clumpy torus model is realized as a collection of individual molecular clumps/clouds arranged in a toroidal structure around the central accretion disk. In reality, this region is likely to be a continuous extension of the outer accretion disk (Elitzur & Shlosman 2006). The primary parameters of the Clumpy torus model are described below.

• 1.  
  N₀: It is the average number of clouds along a radial equatorial ray in a given model. It represents the normalization of a Gaussian distribution of clouds around the equatorial plane. The total number of clouds intersecting a given equatorial ray is different for different lines of sight. The intrinsic AGN continuum can escape along many different lines of sight, and the observed mid-IR 10 μm silicate features can be seen in emission even for lines of sight close to the equator. The total effective optical depth to the continuum source is thus a function of the number of clouds along the line of sight and optical depth of each cloud.
• 2.  
  τ_V: Each of the clouds/clumps is assumed to have the same optical depth τ_V in the V band. Assuming standard Milky Way dust extinction with R_V = 3.1, A_V/τ_9.7μm = 18.0 ± 1.0 (Roche & Aitken 1984; Whittet 2003), and A_V = 1.086 τ_V; only when τ_9.7μm ⩾ 1 or τ_V ≳ 16.5, we will notice the effects of self-absorption on the 10 μm feature.
• 3.  
  Y: The radial extent of the torus, Y, is the ratio of the outer (R_o) to the inner radius (R_d) of the torus. The inner radius depends on the onset of dust sublimation due to the incident UV radiation from the accretion disk (Barvainis 1987). See also Equation (1) in Nenkova et al. (2008b). The radial extent Y of the torus decides the IR turnover at long wavelengths (λ ≳ 30 μm).
• 4.  
  q: The clouds are distributed along the radius with a power-law distribution (r^−q) parameterized with the exponent "q." For q>1, the clumps are concentrated closer to R_d. When the clumps are packed closer to R_d, the resultant IR SED is dominated by the emission from dust close to its sublimation temperature, and there is little long-wavelength mid-IR emission. The corresponding width of the SED (Pier & Krolik 1993) in this case is also small.
• 5.  
  σ: The torus angular width, σ, is the width of the Gaussian distribution of clumps around the equatorial plane. Thick tori (large σ) generate redder 3–8 μm continua (in λF_λ units).
• 6.  
  i: The models produce the IR SED longward of ∼1 μm for each inclination i from 0° (face-on) to 90° (edge-on) in steps of 10°.

The torus models are constructed using the radiative transfer code, Clumpy (Nenkova et al. 2002). The tabulated SEDs for different parameters are accessible from the Clumpy project Web site.⁷

Clumpy dust density distributions differ from smooth distributions in one important aspect: in smooth dust distributions the temperature is uniquely determined by the distance from the source of radiation. While this is also roughly true for clumpy distributions, the presence of lines of sight with different dust columns allows both hot and cold temperature regions to co-exist at similar radial distances. This leads to a greater dependence of the output SED on N₀, τ_V, and q. The primary motivating factor for considering clumpy models for the torus is interferometric observations of local AGNs (Jaffe et al. 2004) which constrain the tori to be physically small (R_o ≲ a few pc).

To fit our data with the Clumpy torus models, we adopt the procedure developed by Nikutta et al. (2009). We analyze the distributions of best-fitting Clumpy torus parameters for each quasar in our sample. We consider the following grid of parameters:

• 1.  
  q = 0.0–3.0, in steps of 0.5
• 2.  
  N₀ = 1–15, in steps of 1
• 3.  
  τ_V = 5, 10, 20, 30, 40, 60, 80, 100, 150
• 4.  
  Y = 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100
• 5.  
  σ = 15–70, in steps of 5
• 6.  
  i = 0–90, in steps of 10.

In all, there are ∼1.25 million possible combinations of model parameters. Very large values of Y and τ_V, present in the original model grid in Nikutta et al. (2009), are excluded here as the objects under study are type 1 quasars with, in most cases, silicate 10 μm feature in emission.

Each model is scaled and fitted such that the overall fitting error E is minimized. We adopt Equation (1) from Nikutta et al. (2009) shown below:

$$\begin{equation} E = \frac{1}{N}\sqrt{\sum ^{N}_{i=1}\left(\frac{ F_{\rm AGN} \cdot f_{i}^{m} - f_{i}^{\rm {obs}}}{\sigma _{i}} \right).} \end{equation} \tag{ 1 }$$

Here, $f_{i}^{\rm {obs}} = \lambda _{i} F_{i}$ are the observed SED data points that are interpolated at model grid points denoted by λ_i (see below for why we take this approach), f^m_i are the corresponding model SED points, and σ_i are the 1σ errors on λ_iF_i. The scaling of the model, F_AGN, provides a measure of the IR luminosity of the Clumpy torus, which can be converted into an estimate of the bolometric luminosity of the system.

For each parameter, we construct a discrete distribution of values by selecting a sample of well-fit models. For each model, the fitting error E is computed from Equation (1). The model with the minimum value of fitting error, E_Min, is considered to be the best-fit model. Further, a relative error, E_r = 100 × |E − E_Min|/E_Min, for each model is constructed. Models that differ by 10% from the minimum value E_Min are considered to represent the distribution of parameter values that best represents the data for a given quasar. For each parameter, we consider the mode of the distribution of parameter values as the most probable value of the parameter for a given quasar. Note that the best-fitting value may not be the most probable one. The 90% confidence intervals for a parameter are also computed.

The model SEDs are scaled and fitted to the data SEDs constructed from the photometric data and the IRS spectrum. We attempt the fitting procedure for all 1.25 million model SEDs, and record their respective relative error E_r. Parameter distributions are then constructed where the acceptance criteria to form the samples are E_r∼ 10%, 20%, and 30%. We find that the distributions gradually become flatter or uniform as the relative error criterion is relaxed. Thus, a narrower distribution suggests a better constrained parameter value.

To measure how well a parameter is constrained we use the discrete Kullback–Leibler (KL) divergence. The KL divergence measures the similarity between two histograms (or discrete distributions) of identical sampling k. The KL divergence is written as

$$\begin{equation} D_{\rm KL} = \sum _{k} Po_{k} * \log _{2}(Po_{k}/Pr_{k}) / \log _{2}(N) \nonumber \end{equation} \tag{ 2 }$$

with N being the number of sampled bins, Pr the prior distribution of parameter values (uniform in this analysis), and Po the posterior distribution of parameter values (histogram of "accepted" parameter values). The normalization log₂(N) ensures that D_KL = 1.0 when all accepted models happen to have a parameter value within a single bin. A D_KL value close to 1 indicates a better constrained parameter.

Figure 2 shows the parameter distributions corresponding to three sources from our sample for brevity. For each parameter there are three figures from left to right corresponding to the three sources: SDSS J100401.27+423123.0, SDSS J160950.72+532909.5, and SDSS J142730.19+324106.4. SDSS J100401 shows strong near-IR emission, and has a high UV luminosity (Figure 1). SDSS J160950 is weak in the UV and also has weak 3 and 10 μm features; SDSS J142730 has a deep 10 μm absorption feature. This source shows a power-law optical/UV continuum in its SDSS spectrum, but its emission lines are absorbed (and may be a BALQSO), which is consistent with its mid-IR nature. For each source, we show the distribution of Clumpy parameter values that forms by accepting models that fall within 10% of the best-fit model.

Zoom In Zoom Out Reset image size

It should be noted that we are fitting the entire SED from UV to mid-IR, and the IRS range is more densely sampled than the photometry. To avoid problems due to uneven spectral sampling, the data SED was resampled to the wavelength grid of the models. This is not a significant issue since we are interested in fits to the broadband features of the SED such as the optical AGN power law, the near-IR bump, and the 10 μm feature. The Clumpy model SED is better sampled near 9.7 μm than elsewhere; thus, this improves the fit to the 10 μm region without biasing the fit to be weighted more by the 2–8 μm continuum. Another important point to be noted is that the selection of the model (AGN+TORUS) SED is also constrained by the optical/UV portion of the data SED. While we do not investigate changes in the intrinsic AGN SED, by including fits to the optical SDSS photometry, we are preferentially selecting model SEDs that satisfy consistent flux density scaling in both UV and mid-IR regime at the same time.

5.1. Model Fits using the Clumpy SED

Initially, we used only the Clumpy model SED to fit the data SED. The best-fit values of the parameters for each model are given in Table 5. Example model fits are shown in the left-hand panels of Figure 3. The best-fitting models of the entire sample have N₀ ∼ 1, τ_V∼ 20–100, q ≳ 2, σ ∼ 15, and i∼ 60–80 (see Table 5 and dark bars in Figure 2). The radial extent Y of the torus is unconstrained with parameter distributions nearly flat over the sampling grid.

Zoom In Zoom Out Reset image size

Table 5. Best-fit Clumpy Torus Parameters without an Additional Blackbody Component

SDSS ID	q	N0	τV	Y	σ	i	χ2ν	EMin
095047.47+480047.3	3.0	1	60.0	5	15	90	2.3023	0.2086
100401.27+423123.0	3.0	1	80.0	5	30	60	7.5916	0.3655
103931.14+581709.4	2.0	2	40.0	5	15	80	1.0076	0.1420
104114.48+575023.9	3.0	1	30.0	5	30	80	1.7045	0.1829
104155.16+571603.0	1.0	1	60.0	5	25	80	4.4097	0.2914
104355.49+562757.1	1.0	3	20.0	5	15	90	7.0194	0.3577
105001.04+591111.9	3.0	2	150.0	5	35	40	7.2815	0.3710
105153.77+565005.7	2.5	1	40.0	5	35	70	4.6752	0.3029
105447.28+581909.5	2.5	2	80.0	10	30	60	2.4313	0.2124
105951.05+090905.7	3.0	1	100.0	5	25	90	0.7722	0.1220
132120.48+574259.4	3.0	1	40.0	5	15	90	8.0009	0.3889
142730.19+324106.4	0.5	15	10.0	5	35	80	2.3793	0.2161
142954.70+330134.7	3.0	1	60.0	5	35	70	0.8621	0.1301
143102.94+323927.8	2.0	4	30.0	60	15	80	0.8203	0.1245
143605.07+334242.6	3.0	1	80.0	10	30	60	1.8428	0.1866
151307.75+605956.9	3.0	1	100.0	5	20	70	8.4899	0.3899
160004.33+550429.9	3.0	4	20.0	60	20	80	3.3760	0.2574
160950.72+532909.5	3.0	1	100.0	5	30	40	4.8608	0.3089
161007.11+535814.0	0.0	1	40.0	5	60	40	2.5430	0.2234
161238.27+532255.0	2.5	2	10.0	40	15	0	38.6937	0.8475
163021.65+411147.1	1.5	2	10.0	10	15	0	90.6623	1.9043
163425.11+404152.4	3.0	1	40.0	5	30	70	1.0808	0.1429
163952.85+410344.8	2.0	1	150.0	5	20	20	1.0997	0.1455
164016.08+412101.2	3.0	2	60.0	5	15	70	5.2334	0.3205
172522.06+595251.0	3.0	3	20.0	5	15	80	1.4518	0.1688

Notes. Descriptions of Clumpy torus parameters: q: index of the radial distribution (r^−q) of clouds; N₀: average number of clouds along radial equatorial rays; τ_V: optical depth through a single cloud at optical wavelengths; Y: the ratio of outer to inner (sublimation) radius of the torus.; σ: the angular width of the torus in degrees; i: inclination of line-of-sight of the observer in degrees; χ²_ν and E_Min provide measures of how well the best-fit model fits the observed data. Typically, χ²_ν close to 1 and smaller values of E_Min indicate a better fit.

Download table as:  ASCIITypeset image

We find that models with N₀ ∼ 1, q ∼ 3, τ_V ≲ 10, Y ∼ 5, and σ ∼ 15 show peaked 10 μm silicate emission features for all values of the inclination of the line of sight. For τ_V ≲ 15, all wavelengths longer than ∼1.5 μm have τ < 1, and the dust emission is optically thin (Nenkova et al. 2008b). In this case, the SED simply follows the shape of the dust absorption coefficients, which decreases rapidly at longer wavelengths in the mid-IR. The observed spectra should then have blue 3–8 μm continua, which is indeed the case for luminous objects like SDSS J100401 and SDSS J151307 (both from program 50328), as can be seen in Figure 1 (bottom row of panels).

Further, q is well constrained in the case of single-component models to a high value of 2–3 in the case of most objects. This suggests a steep radial distribution of clumps, with most clumps concentrated close to R_d. Nenkova et al. (2008a) show that whenever q ≳ 2, Y is fundamentally unconstrained. As most clumps are closer to R_d in this case, the absolute size of the torus does not matter; the output SEDs from tori of all sizes look the same. On the other hand, for sources with q ∼ 0, the clump distribution is flatter/spatially extended, and Y can be constrained much better for such sources as cooler temperatures contribute at longer wavelengths.

Increasing N₀, τ_V, and/or Y causes the SED to become redder in the 2–8 μm wavelength range, and the overall flux density peak shifts to longer mid-IR wavelengths (due to the Wien displacement law). The increasing N₀ and τ_V essentially increases the obscuration due to the torus, and leads to increased contribution from the cooler parts of the clouds. This effect can be seen by comparing best-fit values of N₀ for SDSS J142730 (Table 5) with the rest of the sample. Larger N₀ at smaller τ_V and small Y apparently still produce deep absorption features. A larger Y has a similar effect if q ≲ 1, as clouds are more spread out radially, and hence cooler. Thus, detecting a blue SED in the 2–8 μm range suggests small N₀, τ_V, and Y, along with a radially steep distribution (q ≳ 1) of clouds. This conclusion, however, comes with a caveat: while it is clear that the near-IR emission is generated by the dust close to its sublimation point, the strong silicate emission features predicted by the Clumpy models with these parameter configurations are not observed.

The near-IR emission is fitted well by Clumpy models with N₀ ∼ 1, τ_V ∼ 5, and q>1; the 10 μm feature profiles are not well fit by the same models. The model 10 μm profiles are more peaked than observed profiles, which are broad and shallow. We note that this uncertainty about the origin of the near-IR emission in torus models was also encountered previously in the study by Pier & Krolik (1993), where they also had to employ an additive blackbody component to represent the near-IR contribution separate from their mid-IR torus component. Even in smooth density models, where dust temperatures are functions of radial distance from the source, use of a common sublimation temperature for graphite and silicate dust leads to this effect. Using different sublimation radii for different grain populations is computationally expensive, which could explain some of these discrepancies.

Fitting UV/optical continuum and mid-IR together highlights the need for an additional blackbody component (see the left panels of Figure 3). Polletta et al. (2008, see their Figure 1) also came to similar conclusions in their effort to fit high-z extremely obscured sources with clumpy torus models from Hönig et al. (2006). This appears to be a common problem to all clumpy models constructed so far.

5.2. Additive Near-IR Blackbody Emission

The 10 μm silicate emission feature gets broader and weaker with increasing τ_V, N₀, σ, and Y. None of our objects show peaked 10 μm silicate emission profiles in the spectra, indicating that the hot dust generating the near-IR emission is depleted in silicate dust, and that the 10 μm region receives contributions from multiple "colder than sublimation temperature" sources which likely make the feature broader and weaker. Right-hand panels in Figure 3 show the fits of silicate features in the presence of an extra blackbody component.

In most objects, the feature either peaks around 9.7 μm (SDSS J100401) or is mostly flat (SDSS J160950). In some cases, there is a well-defined plateau from 9.7 to 11.4 μm (SDSS J151307, last rightmost panel in Figure 3). We find that with the Clumpy+Blackbody fits, the models mostly reproduce the observed shapes within the errors of the observations, with the exception of emission around 11.3 μm. This suggests the presence of dust species other than silicates in these quasar spectra (see also Hao et al. 2005; Sturm et al. 2005; Markwick-Kemper et al. 2007). We find that with the exception of excess flux around 11 μm, the silicate features in 14 out of 25 sources are fitted well.

Another issue in fits to the 10 μm features is the observed shift of the feature peak in quasar spectra (see, e.g., Figure 3 of Hao et al. 2005). Whether this shift is a real effect is still uncertain; however, we note that radiation transfer in clumpy media as demonstrated by the fits in this paper may explain the varied shapes and apparent shift of the feature peak.

We present Spitzer/IRS observations of a sample of optically luminous type 1 quasars at z ∼ 2. Their rest-frame 2–12 μm IR spectra show two prominent features peaking at ∼3 and 10 μm. The 10 μm feature is the 10 μm silicate emission feature, commonly observed in Spitzer observations of other type 1 AGNs (Hao et al. 2005; Siebenmorgen et al. 2005; Sturm et al. 2005). The 3 μm bump is the expected signature of the hottest thermal dust emission from the inner region of the dust torus. There is a strong correlation between the optical/UV and IR luminosities (Gallagher et al. 2007), and the detection of this near-IR bump in a sample of optically luminous high-redshift quasars shows that the optical/UV continuum heats the dust in the inner torus, which then radiates in the thermal near- to mid-IR.

We fit the spectra and the UV-to-mid-IR SED with Clumpy torus models (Nenkova et al. 2008a). This is the first time such fits have been attempted to spectroscopically confirmed high-z quasars with near-IR data. We considered two different approaches. In the first case, we use the Clumpy model SED. These Clumpy torus models provide good fits to the 2–8 μm part of the spectrum, if we only fit data longward of 1 μm. Models with average number clouds along a radial equatorial ray (N₀) ∼1, optical depth through each cloud (τ_V) ≲10, and a radial distribution of clouds (r^−q) described by a power-law exponent (q) ∼3 fit IRS spectra (not complete SEDs) with a strong hot-dust bump very well. The q ∼ 3 values suggests that the hot dust component is more centrally concentrated as expected. However, the 10 μm silicate emission features of these models show strongly peaked profiles, and the 10 μm features in the observed spectra are more broad and flat. This problem can be partially removed by fitting the entire SED from UV to mid-IR; using this long lever-arm, the Clumpy model SED is consistently weaker than the observed SED in the 1–7 μm range (see the left panels of Figure 3), highlighting the lack of additional near-IR contribution in the models, if both UV and IR data are fitted together.

To accurately model the 10 μm silicate emission features and remove the above inconsistency, we considered the Clumpy+Blackbody model where we fit the spectra and the SED with a linear combination of a hot dust blackbody and a Clumpy model. In these fits, the clumpy models provide good fits to the 10 μm region, while the blackbody contributes more strongly to the region between 2 and 8 μm. Use of the additional blackbody leads to a stronger contribution of the Clumpy model to the far-IR emission. Whether this is a real effect may be tested using far-IR facilities like Herschel.

We compared the IR properties of this sample to the low-redshift PG quasar sample (z ∼ 0.1) from the Spitzer archives, and find that the primary difference in the 2–8 μm range between low- and high-redshift samples is the absolute luminosity. There are, however, significant object-to-object differences in the 10 μm silicate emission features, which point to real differences in the dust structure of their tori. In a few cases, such as SDSS J142945, the 9.7 μm peak of the silicate feature appears shifted to longer wavelengths. Just as other observations have noted the presence of different dust species (Hao et al. 2005; Sturm et al. 2005; Markwick-Kemper et al. 2007), we note a feature around 11.3 μm in some sources that may be due to crystalline silicates (Markwick-Kemper et al. 2007).

The 10 μm feature shapes in 14 out of 25 objects are well reproduced by Clumpy models; the agreement is weak in other cases mostly due to the lack of a clear emission feature. Presence of additional dust species also seems to contribute to this issue. More work is necessary to connect the near-IR emission with the rest of the torus structure. The lack of near-IR contribution in the torus models with clumpy media (in general) appears to be rooted in not considering the balance of amounts of silicate and graphite grains as a function of distance from the source.

However, we find that the near-to-mid IR SED analysis is a powerful tool for distinguishing between different distributions of q, N₀, and τ_V in Clumpy models. Observing a blue 3–8 μm continuum indicates that the source is compact (q>1) with N₀ ∼ 1. A redder continuum may require a more extended (q < 1) distribution of clumps with N₀ ∼ 10 and τ_V ∼ 30. Further, improvements in fits using the complete UV-to-mid-IR SED suggests the importance of using UV/optical data if available. Further far-IR data where the contribution from cold dust associated with star formation in the host galaxy of the quasar may be dominant (Netzer et al. 2007) are also important. The radial extent of the torus (Y) is constrained by the location of the far-IR turnover in the IR SED; however, contribution from cold dust in the host galaxy is also dominant in the same region, disentangling these contributions will be interesting (see, for example, Hatziminaoglou et al. 2010).

In a Clumpy torus, the probability of viewing the AGN as a type 1 object depends more strongly on N₀ and τ_V than on the inclination to the line of sight i. Using multi-component models decreases this sensitivity of the model SED to parameters like N₀. This is observed in the number of accepted models in Table 7; even for objects with S/N ∼25 (SDSS J100401, SDSS J151307), the number of accepted models is ≲1000. The argument in favor of Clumpy+Blackbody models is that they represent the complete data range better, and adding a blackbody component improves the fits to the 10 μm region (right panels in Figure 3), even in the case of objects like SDSS J142730 that should be dominated by the Clumpy model alone.

Addition of the blackbody component to represent the near-IR emission does not by itself represent a failure of Clumpy models, but suggests that more detailed treatment of the origin of the near-IR emission is required. The composite grain approximation assumed in radiative transfer calculations (DUSTY; Ivezić et al. 1999) may lead to stronger 10 μm features than would be generated in the actual dust sublimation transition region. This effect is also seen in models of Schartmann et al. (2005) that use the standard Mathias–Rumpl–Nordsieck dust grain mixture, and obtain strong 10 μm emission features in their SEDs. As the models fits in this paper show, Clumpy models can reproduce the 10 μm shapes adequately. Differences in number density of dust grains of different sizes and compositions with distance from the continuum source likely contribute to the nature of near-IR emission. This dust sublimation region may also be spread out over an extended region rather than in a thin AGN-facing layer of the cloud as assumed in Clumpy models. Future clumpy torus models should consider both these effects to properly model the near- to mid-IR SEDs of active galaxies.

We thank the anonymous referee for their comments that significantly improved this paper. This work is based on observations and archival data from the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA. Support for this work was provided by NASA through awards (RSA's 1353801 and 1365236) issued by JPL/Caltech. G.T.R. was supported in part by an Alfred P. Sloan Research Fellowship. S.C.G. thanks the National Science and Engineering Research Council of Canada for support. The IRS was a collaborative venture between Cornell University and Ball Aerospace Corporation funded by NASA through the Jet Propulsion Laboratory and Ames Research Center. SMART was developed at Cornell University and is available through the Spitzer Science Center at Caltech. This research has also made use of NASA's Astrophysics Data System Bibliographic Services.

Facility: Spitzer - Spitzer Space Telescope satellite