Dusty Structure Around Type-I Active Galactic Nuclei:
Clumpy Torus Narrow Line Region and Near-Nucleus Hot Dust
We fitted Spitzer/IRS spectra of 26 luminous QSOs in attempt to define the main emission components. Our model has three major components: a clumpy torus, dusty narrow line region (NLR) clouds and a blackbody-like dust. The models utilize the clumpy torus of Nenkova et al. (2008) and are the first to allow its consistent check in type-I AGNs. Single torus models and combined torus-NLR models fail to fit the spectra of most sources but three component models adequately fit the spectra of all sources. We present torus inclination, cloud distribution, covering factor and torus mass for all sources and compare them with bolometric luminosity, black hole mass and accretion rate. The torus covering factor and mass are found to be correlated with the bolometric luminosity of the sources. We find that a substantial amount of the radiation originates from a hot dust component, which likely situated in the innermost part of the torus. The luminosity radiated by this component and its covering factor are comparable to those of the torus. We quantify the emission by the NLR clouds and estimate their distance from the center. The distances are times larger than the dust sublimation radius and the NLR covering factor is about 0.07. The total covering factor by all components is in good agreement with the known AGN type-I:type-II ratio.
Subject headings:infrared: galaxies – galaxies: active – galaxies: nuclei – quasars: general
An obscuring dusty structure surrounding an accreting massive black hole (BH) is believed to be a common feature of most active galactic nuclei (AGNs). Since the obscuration of the central region is anisotropic, sources with small inclination angles to the line of sight (face-on) would be classified as type-I AGNs and those with large inclination angles (edge-on) would be classified as type-II AGNs. In this picture, the bulk of the radiation from the central engine is absorbed by the obscuring structure and re-emitted mainly in mid-infrared (MIR) wavelengths. The central obscuration is not necessarily a single component structure and the exact nature of its different components remains an open question and provides the motivation to the present work.
A main component of the obscuring structure is believed to be a dusty torus. The MIR spectral energy distribution (SED) of such torus depends on its dimensions and geometry, the density distribution and the dust grain properties. Initial attempts to model such tori assumed smooth density distributions (e.g. Pier & Krolik 1992; 1993; Granato & Danese 1994; Efstathiou & Rowan-Robinson 1995; van Bemmel & Dullemond 2003; Schartmann et al. 2005). This can explain part of the SED but falls short of providing realistic MIR spectra. For example, Silicate emission has been detected in Type-2 AGNs (e.g. Sturm et al. 2006; Teplitz et al. 2006) despite of the the fact that absorption features are expected in edge-on tori. Several studies suggested that the dusty medium should be clumpy (e.g. Krolik & Begelman 1988; Rowan-Robinson 1995; Nenkova et al. 2002; Tristram et al. 2007; Thompson et al. 2009). The recent works of Nenkova et al. (2002) and Nenkova et al. (2008a; 2008b; hereafter N08) offer the required formalism to calculate the SED of such tori and a framework to constrain some of their properties.
Other AGN components can contribute to the observed MIR spectrum of AGNs. Some of this emission may originate farther from the central radiation source, at distances exceeding the dimensions of the torus. Broad-band 10 m imaging of several nearby AGNs suggests extended MIR continuum source (Cameron et al. 1993; Bock et al. 2000; Tomono et al. 2001; Radomski et al. 2003; Packham et al. 2005). Dusty clouds in the narrow line region (NLR) may be the source of such radiation (Schweitzer et al. 2008; hereafter S08). Thus, a significant contribution at 10–30 m due to components not related to the torus, must be considered.
|PG 0157+001 (Mrk 1014)||0.1630||779||44.68||2,3,4,5,16,19,20,21|
|PG 0050+124 (IZw1)||0.0611||273||44.30||6,7,9,10,12,16,22,23,24|
References. – (1) Rudy et al. (1982); (2) 2MASS magnitudes (Jarrett et al. 2000); (3) Neugebauer et al. (1987); (4) Guyon et al. (2006); (5) Glikman et al. (2006); (6) Sanders et al. (1989); (7) Elvis et al. (1994); (8) Neugebauer et al. (1979); (9) Balzano & Weedman (1981); (10) Rieke (1978); (11) McAlary et al. (1983); (12) Stein & Weedman (1976); (13) Surace et al. (2001); (14) Hayland & Allen (1982); (15) Gavazzi & Boselli (1996); (16) Veilleux et al. (2006); (17) Matsuoka et al. (2005); (18) Sitko et al. (1982); (19) Scoville et al. (2000); (20) Surace et al. (1999); (21) Imanishi et al. (2006); (22) Spinoglio et al. (1995); (23) Jarret et al. (2003); (24) Allen (1976)
Another component not necessarily related to the torus is hot dust emission in the immediate vicinity of the central engine. Minezaki et al. (2004) reported delayed and correlated variations between the V and the K band emission in the nucleus of the Seyfert 1 galaxy NGC4151. The measured lag between the V and K bands, days, lead to the conclusion that the near infrared emission is dominated by thermal radiation from hot dust 0.04 pc from the center. This result was recently supported by the work of Riffel et al. (2009) who fitted the near infrared spectrum of NGC 4151 with blackbody (BB) spectrum representing emission from hot dust in the inner region of the torus. Other studies used similar models to fit the SED of more luminous AGNs (e.g. Edelson & Malkan 1986; Barvainis 1987; Kishimoto et al. 2007). The study of Suganuma et al. (2006) indicates that luminous and variable K-band emission is common in several nearby Seyfert-1 galaxies. More generally, observational support for powerful 1–3 m emission is known for decades (e.g., Hyland & Allen 1982; McAlary et al. 1983; Neugebauer et al. 1987). It is not clear that such emission is consistent with the torus dust emission, which is expected to peak at longer wavelengths. Furthermore, part of this radiation requires dust temperature that exceeds the sublimation temperature of silicate type grains suspected to be responsible for the bulk of torus MIR radiation.
The high quality spectra made available by the IRS spectrometer on-board the Spitzer Space Observatory (Houck et al. 2004) allows the detailed analysis of the MIR SED of a large number of AGNs. Here we fit the observed MIR spectra of 26 PG QSOs, already investigated by S08, using more realistic three component models made of a clumpy torus, dusty NLR clouds and very hot dust clouds. In §2 we describe the observational data. In §3 we detail our model and the fitting procedure and in §4 we present the results of this procedure. Our main findings are summarized and discussed in §5.
2. Sample Selection Spitzer Observations and Data Reduction
Our sample consists of all AGNs in the QUEST Spitzer spectroscopy project (PID 3187, PI Veilleux). It is described in detail in Schweitzer et al. (2006; hereafter S06), Netzer et al. (2007) and S08. Most of the objects are Palomar-Green (PG) QSOs (Schmidt & Green 1983) and are taken from Guyon (2002) and Guyon et al. (2006). The luminosity range is where stands for at rest wavelength 5100Å. Radio loudness and infrared excess are typical of these properties in low redshift PG QSOs. The QUEST sample, while representing the Guyon et al. (2006) sample, includes fewer sources and thus is not complete. This was explained in S06. Thus, several of the correlations discussed below may differ from those in the entire sample in a way which is difficult to asses until a larger data set, with the same spectral coverage, is obtained. Here we aims at modeling the entire 2-35 m SED hence we omit PG 1307+085 (which appeared in the S06 sample) due to lack of observations in the 5.2-8.7 m range (Spitzer SL1 mode). Table 1 lists names, redshifts, and for all 26 QSOs in our sample.
The Spitzer observations and the data reduction of the QUEST sample are detailed in S06 and are summarized here for completion. IRS spectra for all objects were taken in both low resolution (5-14 m) and high resolution (10-37 m) modes. The standard slit widths of 3”.6 to 11”.1 include flux from the QSOs hosts and the vicinity of the AGNs. Data reduction starts from the basic calibrated data (BCD) provided by the Spitzer pipeline. Specially developed IDL-based tools are then used for removing outlying values for individual pixels and for sky subtraction. The SMART tool (Higdon et al. 2004) is used for extraction of the final spectra. More details are given in S06.
We have supplemented the Spitzer spectra with near infrared (NIR) data obtained from the literature. The data are obtained from various sources, in particular Neugebauer et al. (1987) and the 2MASS extended and point source catalogs (Jarrett et al. 2000). For some of the sources there are several photometric measurements. For these, we average all measurements in each band and used the standard deviation as the flux error. For single observations, the flux error is estimated as 20%. These uncertainties must be underestimated since many of the NIR data were taken about 20 years before the Spitzer observations and all sources are variables. For a list of JHKL data references see table 1.
We computed the bolometric luminosity of all sources using their and a simple bolometric correction factor, BC, to convert it to . The best values of BC are discussed in Marconi et al. (2004), Netzer et al. (2007), and Netzer (2009). The approximation used here is taken from Netzer (2009) and is given by BC, where / erg s. For most objects the value of is obtained from the observations of Boroson & Green (1992). For PG 1244+026, PG 1001+054, and PG 0157+001 we use spectra from the SDSS (York et al. 2000) data release seven (DR7, Abazajian et al.2008) that were measured in the way described in Netzer and Trakhtenbrot (2007).
There are two uncertainties associated with the use of . The first is source variability, which is an important effect since the optical and MIR spectroscopy were separated by many years. We estimate this uncertainty to be a factor of . The second uncertainty involves the approximation used for BC. We estimate this uncertainty to be 30%. Both uncertainties affect the derived model parameters such as torus covering factor and NLR distance. Because of this, we do not attach great importance to specific values in specific sources obtained from our best acceptable models. On the other hand, the sample is large enough to enable a significant analysis of its mean properties since the larger of these effects, due to source variability, is changing in a random way. The uncertainty due to the estimate of BC is smaller but more problematic since the expression we use can under-estimate or over-estimate for all sources. This can introduce a systematic difference in torus and other component properties, e.g. a smaller covering factors for all sources. We return to this point in §5.1.
3. Spectral Modeling
3.1. Spectral Components
The aim of the present work is to fit the 2-35 m spectra of the sources in our sample. For this we combine models of three different physical components; a dusty torus, a dusty NLR and a very hot dust component. An additional starburst component representing emission from the host galaxy, is also present but is not part of the fit procedure. This component is subtracted from the observed spectra using the MIR spectrum of M82 in a way similar to the one used in S06. The reality of the subtraction can be judged from the remaining emission or absorption in the wavelengths corresponding to the strong PAH features (e.g. 6.2 or 7.7 m). More details are given in S06 and in §3.2. This section describes the spectral properties of the three components and the following sections give detailed account of the fitting procedure.
The first component represents a dusty torus surrounding the central energy source. AGN unification schemes suggest that the central engine is surrounded by dusty, optically thick, toroidal structure (e.g., Krolik & Begelman 1988; Antonucci 1993). Earlier torus models assumed a smooth density distribution (Pier & Krolik 1992; Rowan-Robinson 1995; van-Bemmel & Dullemond 2003). They all suffer from various incompletions and their agreement with MIR spectral observations is generally poor. The more recent clumpy torus models of N08 represent a significant improvement and are the basis for the present study.
The fundamental difference between clumpy and smooth density distributions is that the former allows the radiation to propagate freely between different regions of the optically thick medium. The clumpy dust distribution results in the coexistence of clouds with a range of dust temperatures at the same distance from the central radiation source. They also allow clouds at large distances to be exposed to the direct AGN continuum. In contrast, smooth density distribution models associate a certain dust temperature with a certain distance from the central source, thus limiting the range of acceptable SEDs. These differences allow the clumpy torus models to have a range of spectral properties not accessible in smooth density distribution torus models. N08 describe the formalism and the detailed radiative transfer used to calculate the MIR spectrum of clumpy dusty tori and our use of the torus models follows the same procedure.
The first parameter of the clumpy torus model is the inner radius of the cloud distribution that is set to the dust sublimation radius . This corresponds to a dust sublimation temperature that depends on the grain properties and mixture. Here we adopt two different sublimation radii appropriate for two types of grains. The first is
where . This gives the innermost radius where graphite dust can survive. It corresponds to a sublimation temperature of about 1800K and is the one used in S08. N08 adopted a somewhat larger distance, which is more appropriate for silicate type grains with a sublimation temperature of 1500K. This is given by
In the following we specify which of the two is used for each purpose. The clumpy torus model requires six additional parameters:
The visual (5500Å) dust optical depth of a single cloud, (all clouds are assumed to have the same ).
The mean number of clouds along a radial equatorial line, .
The ratio between the outer dimension of the torus and , Y.
The inclination angle of the torus with respect to the line of sight, .
The torus width parameter , which is analogues to its opening angle.
A parameter q that specifies the radial power-law distribution of the clouds, i.e. , where is the number of clouds.
N08 also define an additional parameter that is set by the above free parameters. This is the probability that light from the central source will escape the obscuring structure at a given angle without interacting with the clouds. Assuming that individual clouds are optically thick
where . By integrating over all angles and subtracting from 1, we obtain the probability of absorption by the torus.
This is also the fraction of obscured objects in a random sample (e.g. the fraction of type-II AGNs out of the entire AGN population). is equivalent to the “real” (geometrical) covering factor of the torus, and does not depend on the inclination angle since it reflects the global torus geometry. In other words, is the ratio between the total torus luminosity and . Due to the anisotropy of the torus radiation, differs from the apparent covering factor of the torus deduced from the ratio between its observed (angle dependent) luminosity and the bolometric luminosity of the central source. This apparent covering factor can be represented by
where is the angle dependent monochromatic luminosity of the torus. For a more detailed description of all torus parameters, as well as other model properties and assumptions, see N08.
The emission and absorption properties of the torus depend on the dust grains composition and other properties. All models consider here assume spherical dust grains with MRN size distribution (Mathis, Rumple & Nordsieck 1977). The dust has a standard galactic mix of 53% silicates and 47% graphite. The optical properties of graphite are taken from Draine (2003) and for silicates from Ossenkopf, Hennig & Mathis (1992) (OHM). The OHM dust mixture produces better agreement with observations of the 10 and 18 m silicate features (Sirocky et al. 2008) and is the only one used in the present work.
The second component of the model represents a collection of dusty NLR clouds. The motivation for this component is explained in S08 where it was shown that such a component can contribute, significantly, to the MIR flux of luminous AGNs. The properties assumed here for these clouds are similar to the ones used in S08. We assume constant column density clouds with . We further assume constant hydrogen density of 10 cm, solar composition and galactic dust-to-gas ratio. The important physical parameters for this component are the cloud-central source distance (which determines the dust temperature), the incident SED and the dust column density. The assumed gas density and column density only serve to define the emission from this component using convenient parameters such as the ionization parameter. The NLR component is assumed to be concentrated in a thin spherical shell with a small covering fraction. These assumptions (that are reviewed later) mean that a single thin shell with the chosen properties is contributing to any 3-component model.
The calculations of the emitted IR continuum of the clouds, as well as their emission line spectrum, are obtained with the photoionization code ION (Netzer 2006 and references therein). The assumed SED of the central continuum is the one described in S08 and the resulting NLR spectra are similar to those shown in S08 Fig. 1. The clouds are optically thin to MIR radiation and hence no radiative transfer is required for the calculations of this emission. Obviously, this is not the case for the UV part of the spectrum where the opacity is much larger and where the transfer approximation we use is the one built into ION.
The third component is a single BB representing hot dust emission. The need for such a component has been noticed in earlier works, e.g. in S08. This is very clear from the data, where a rise towards short wavelengths can be seen from the IRS spectra even without the supplemented NIR photometry. The simplified BB assumption is not entirely consistent with hot dust clouds illuminated from one side since for large dust optical depth, there must be a temperature gradient across such clouds. For small dust optical depth, the dust temperature is more uniform but the cloud may be partly transparent to the incident radiation. For the purpose of the present study, we assume that all the incident optical-UV radiation is absorbed by this component and the dust temperature is constant. We discuss this component further in §4.3 and §5.3 below.
3.2. Fitting Method
As shown in S08 and explained in detail in Netzer et al. (2007), the starburst contribution to the MIR spectrum can be significant, especially at long wavelengths. To remove this contribution, we follow the procedure of S08 who fitted and subtracted a “nominal” M82 spectrum from all spectra prior to the model fitting. We use the ISO-SWS mid-IR spectrum of M82 from Sturm et al. (2000) and the scaling factors of S08 that were obtained from the intensity of the strong PAH emission features that are clearly seen in some of our objects. Although this component is subtracted prior to the model fitting, its normalization introduces another degree of freedom to the procedure. The remaining spectra are assumed to represent the intrinsic AGN continuum.
The fitting procedure minimizes the modified ,
where is the number of degrees of freedom, is the starburst subtracted luminosity deduced from the observed monochromatic flux density at rest wavelength , and is the combined emission due to all three components,
where and are fitting parameters. The bolometric luminosity of the source uniquely determine the total emergent flux of the torus component through the radiative transfer calculation. Since an independent measure of the bolometric luminosity is available for these type-I sources, the normalization of the torus component is not a free parameter of the fitting procedure (i.e. the first term in Eq. 7 is set by the torus parameters but the total energy emitted by the torus is fixed.
To reduce the uncertainty, we first bin the observed spectra into 100 bins of equal energy widths and take in each the mean observed flux. The torus, NLR and BB models are then interpolated to the same energy grid. The binning of the data practically removes all emission lines and smooth the spectrum, allowing to focus on the broad band continuum shape. Some remaining peaks may result in larger local values but since none of the models contains lines, the effect would be the same for all fits.
According to S08, most of the errors in the IRS spectra are systematic and are approximately 10% of the flux. We add this systematic error to the standard deviation of the data points in each energy bin to get the value of in Eq. 6. This is taken to represent the error on the observed binned flux. The use of this affects the derived such that its absolute value is different from the one normally used in statistical analysis. However, the relative value of can still be used to discriminate between different models and will be treated as such in the remaining of the paper.
The fitting algorithm computes values for all possible combinations of torus, NLR and BB models. There are six free parameters in the torus model that describe its geometrical properties (see section 3.1). The radial extent of the torus, Y, ranges from 5 to 60 . The total number of clouds along radial equatorial lines ranges from 1 to 10. The optical depth of individual cloud, , ranges from 5 to 100. The torus width parameter, , ranges from to . The index of the power law radial distribution of clouds, , only accepts one of the three integer values, 0, 1 or 2.
Finally, the torus inclination angle, , ranges from to . In the NLR model there are two free parameters, the cloud distance and the normalization factor . The distance is changed in steps of 0.075 dex between 1 pc to 850 pc for a source with . For the hot BB, there are also two free parameters, temperature and a normalization factor . The range in temperatures is 800–1800 K.
4.1. Torus only models
We first examine the possibility that single torus models can explain the observed spectrum. This is important in order to evaluate the general properties of clumpy tori models and has never been attempted on large Spitzer -IRS samples. In general, allowing for a torus-only model results in large and poor fit to almost all spectra. On average, values of such fits are about an order of magnitude larger than those obtained by using all three components. None of the torus models could fit well the full observed spectrum of any of the objects. The best torus-only models are those fitted to the spectra of PG 1626+554 ,PG 1004+130 and PG 0838+770. These objects are fitted reasonably well around the silicate emission features and at long wavelengths but fail at m where they show clear flux deficiencies (see left panel of Fig. 1). Single torus fits to all other objects are significantly worse. This is illustrated in the right panel of Fig. 1 that shows the best torus-only fit to the spectrum of PG 1440+356. Not only that the model poorly fits the silicate features region, it all shows insufficient flux both at shorter and longer wavelengths.
4.2. Combined torus-NLR models
Next we examine the quality of two component models. In general, the fits improved compared with torus-only models especially over the wavelength range 6 m. An example is shown in Fig. 2 where we reproduce the best two component fit to the spectrum of PG 1411+442. In this case, the fit is within 10% of the observed flux for all 5 m. Clear deviations are still found at the short wavelength part of the spectrum. We were able to obtain satisfactory fits over the 6–35 m range for 17 out of the 26 objects using two component models. The values for these fits are, on average, a factor 6 smaller than the single torus fit values. These results are comparable to those shown in S08 where a combination of several BBs and NLR components was fitted to the 5–35 m spectrum of all sources. Only in PG 1626+554, which is not part of the above group of 17 objects, the torus-only fit is superior to the two-component fit (i.e. the NLR component is not required). For the remaining 8 objects, the values of the two-component fits are only about a factor 2 smaller than the values obtained for torus-only fits. In these 8 fits, none of the main spectral features are well fitted.
In all cases fitted with by two component models, a significant fraction of the flux shorter of 6 m is still unaccounted for. The mean missing fractional (i.e. the difference in flux between model and observations over the 2–35 m range) is about 40%.
|Object||Torus-only fit||2-component fit||3-component fit|
4.3. Three component models
The missing flux in the short wavelength part of the spectrum in most of our sources is the main motivation for our third, hot dust component. This was not discussed in S08 who fitted only the m part of the spectrum.
We have investigated two possibilities: The first is emission due to hot dust with the same grain composition used for the other components. This can represent additional dusty clouds in the innermost part of the torus that are not part of the clumpy structure of N08. Including such a component, obtained from the single cloud models of N08, resulted in poorer agreement between model and observations. The reason is that such a component produces strong silicate emission features, especially the broad 10 m feature, that completely change the SED. The apparent strength of the 10 m and the 18 m features in pure emission models represents the average temperature of the dust. The introduction of very hot silicate dust gives the impression of a very strong 10 m feature compared with a weak 18 m feature, which clearly contradicts the observations. Since the hot dust luminosity is very significant, the variations in the intensities of the other components cannot compensate for this increase. Thus, using hot silicate dust we could not successfully fit the observations.
The second possibility is a hot BB. The inclusion of the this component improves, significantly, the quality of most fits. It allows for a better fit at both short and long wavelengths since the torus and NLR components are no longer needed to account for the short wavelength part of the spectrum. For the same reason, the models come closer to the observations at 10 m and 18 m. The values obtained with the three-component fits are, on average, an order of magnitude smaller than the torus-only fits and about a factor of 3 smaller compared with the two-component fits. values for all objects and all types of fits are listed in table 2 (note again that these values do not represent the formal values). Given all three components, we obtain fits that are within 20% of the observations over the entire wavelength range. In the wavelength range of the silicate emission features ( m), 19 of our best fit models deviates by no more than 5% from the data. In the remaining 7 sources, the deviation over this wavelength range is at most 10%. A more realistic model for the short wavelength component is a collection of hot graphite-only clouds (see §5.3). Such clouds have no silicate features in their spectra. Model fitting including such a component is deferred to a future publication.
Fig. 3 shows the best fit three-component models for two representative cases, B2 2201+31A and PG 1617+175. The top panel of each diagram shows the best fit model (in red) and the observed spectrum (in black). We also show individual components: torus (dashed line), NLR (thick black line) and hot BB (dash-dotted line). Asterisks represent the K and L photometry. In the bottom and middle panels of each diagram we show the quality of the fit in each wavelength bin.
Having obtained satisfactory fits for all sources, we have calculated the median contribution, over the entire sample, of each of the components. This is done at every wavelength and is shown in Fig. 4 in such a way that the sum of all components at each wavelength bin is 1. The diagram shows that the BB component dominates the spectrum below 4 m while the torus dominates between 5-25 m. The NLR component has a non-negligible contribution (%) above 16 m, including the 18 m emission feature.
4.4. Fit quality
Each of the fits is assigned a quality flag set by four criteria, the calculated over the entire 2–35 m range and partial values for three shorter segments of the spectrum. Regarding the first, since the errors are not normally distributed, the best has a different meaning from the commonly used and setting a standard confidence level is not applicable. Given this, we decided to consider only fits with values not exceeding the minimum by more than a factor of 1.25. This conservative narrow range of acceptable was chosen to ensure that it contains, beside the best fit, only acceptable fits. Indeed all values within this range seem very good and in most cases, there are also a few acceptable fits outside of this range. However, the distribution of the model parameters for the acceptable fits is not affected much by these outliers.
The three additional partial criteria are aimed to identify those fits where a certain part of the spectrum is not fitted well even though the global is within the above mention range. For this, we examine the goodness of the fit in three different parts of the spectrum by calculating values for the ranges m, m and m. We used a similar approach (in this case a factor of 1.5 relative to the best ) and reject fits that are outside of this range.
Having fixed the limits, we selected the fits that pass all four criteria (one global and three partial limits). There are thousands of possible torus models and a much larger number of possible combinations hence, the number of accepted fits is large, typically 10-50 per sources. All these are inspected visually to verify the success of the automatic selection process. The groups of acceptable models, for all sources, is the basis for the following analysis.
Two objects in our sample, PG 0157+001 and PG 0050+124, have total covering factors larger than unity. This can perhaps be explained by source variability. To check this for the entire sample we compare the Boroson & Green (1992) data with more recent observations from the SDSS. New spectra are available for 8 objects in our sample, 5 of which exhibit an increase in luminosity by a factor , two show a decrease of about the same factor and one shows no significant change. Of the two objects with a covering factor larger than unity, only PG 0157+001 has recent SDSS spectra. A comparison with the older observations show no significant change in luminosity between the two epochs. Under-estimation of the star formation contribution to the mid-infrared range of the spectrum may also result in an over estimation of the covering factor. Obviously, we could not fit the spectra of these two objects with adequate quality, using the measured value of their bolometric luminosity (see the discussion regarding the implications of related uncertainties in §5.1). To proceed, we assume that flux variation is the cause for the apparent discrepancy and artificially lower the bolometric luminosity of these objects so that their covering factors become 1. This allows us to fit their spectra in a way similar to the other sources and get adequate fits and reasonable model parameters. However, we do not include them in the overall statistics and all parameter distributions and mean values presented here do not include these objects.
Having found the best three-component models for all individual spectra, we now discuss the properties of these components as well as their distribution in the sample.
5.1. Torus properties
We examine the distribution of the different torus parameters for each of the objects in our sample. Some of the parameters have a narrow distribution around a mean value and in some cases only a single acceptable value. Other parameters exhibit a broader, more uniform distribution. In Fig. 5 we present a typical example of the acceptable torus parameters for one source, PG 2349-014. As seen in this case, , , and have very narrow distributions for all acceptable torus models. ,, and have broader distributions around their mean values. The values obtained for individual objects from their best fit models are listed in table 3 The parameter distributions of all sources were combined by giving each value within the acceptable range in a certain source its relative weight in the distribution. The results of this combination are general parameter distributions that are shown in Fig 6. Inspection of the various parts of this figure lead to the following conclusions:
The mean torus width parameter is . None of the fits requires the largest allowed of . This value determines, more than any of the other parameters, the geometrical covering factor () of the torus and is consistent with the expected type-I:type-II ratio (see below).
The mean radial extent of the torus, Y, has a broad distribution with a mean value of 31. The range in Y implies torus outer sizes, which range from 1 to 35 pc using and 3 to 90 pc using .
The average number of clouds along an equatorial ray, , has a broad distribution with a mean of 5 clouds. The distribution is practically uniform for .
The parameter that specifies the radial power-law distribution of the clouds, q, has a mean value of 1. This parameter is related to the anisotropy of the torus radiation. As q decreases the torus radiation becomes less isotropic. The SED shape is affected less by q and hence we do not attach great significance to its specific value.
The mean optical depth of a cloud is =58. This parameter, again, has a broad distribution covering the range . Since the requirement for large MIR optical depth is , it is not surprising that the torus models are not very sensitive to the exact value of beyond this value.
The torus inclination angle shows a broad distribution for all with a mean of . This again is consistent with the assumption that the direct line-of-sight to the center of type-I AGNs is almost completely free of obscuring material.
We also checked for various correlations between those parameters and found no such correlation. This suggest that torus models with different parameters can result in similar SEDs.
As explained, SED related uncertainties, in particular the bolometric luminosity of the sources, make it difficult to determine the exact torus properties in individual objects. There is a relatively large range in some of the parameters (e.g. and ) and small changes in can result in a significant change in those parameters. For example, lower bolometric luminosity, due to changes in between the optical and the Spitzer observations, requires lower and higher inclination angle. In principle this may introduce non-negligible uncertainties on (and other parameters) in individual sources. However, the average sample properties are less likely to depend on such uncertainties and the discussion below focus on this aspect of the work. Below we examine three of the parameters, , , and using Eq. 3. We set two of the parameters to their mean values and investigate the dependence of on variations of the remaining free parameter. This is equivalent to the investigation of a generic torus that represents the entire sample of 26 AGNs.
|Real total CF||0.0033(-)||0.024(-)||0.0056(-)||0.045(-)||…|
|Apparent total CF||0.0011(-)||0.0063(-)||0.0013(-)||0.014(-)||…|
Note. – Plus and minus signs indicate positive and negative correlations, respectively. Results for are not listed.
First we consider the torus inclination angle as the free parameter. Using Eq.3 together with the mean values found here, the probability that light from the central region will escape such a structure, and the source will be classified as type-I AGN, is about 75%. Setting and to their mean values (5 and , respectively) and changing the inclination angle, we find that for and for . We therefore expect that for type-I objects, the inclination angle should not drop below . As can be seen in Fig.6, the distribution of the inclination angle is indeed consistent with random inclination angles in the range with no preference to a specific angle.
The second important parameter is the mean number of clouds, . Setting and to their mean values and changing , we find that remains large even for very high values of . Thus, the exact value of is not very important provided the inclination angle is not too large. This is the reason why has the broad distribution seen in Fig.6.
Next we tested the acceptable range of the torus width parameter, , by fixing all other parameters to their mean values. We find that is sensitive to such changes and falls rapidly for . This is the reason that the distribution falls rapidly for . Torus with high value of would have low and thus a low probability of being classified as type-I AGN even if the inclination angle is zero (a face-on torus).
The combined effect of the above parameters is, perhaps, a more meaningful test. In particular, the combination of and determine the geometrical covering factor of the torus, , which is combined, later, with the equivalent property of the other components. This issue is discussed in §5.4.
The above procedure is not applicable to the visual optical depth parameter. This parameter does not affect or but influence the SED shape, in particular the apparent strength of the silicate emission features. The mean visual optical depth for the entire sample is =58 and the lowest value is . Thus, realistic clumpy torus models require large IR optical depths to explain the observed spectrum. The exact value of is not very important once the clouds are thick enough.
The torus size found here is consistent with upper limits set by several high resolution observations of nearby AGNs (e.g, Soifer et al. 2003; Radomski et al. 2003; Jaffe et al. 2004). This is in contrast to previous theoretical works involving smooth density distribution torus models. In these models, the torus is required to be much larger (up to few hundreds of pc) in order to posses the large amounts of cool dust necessary for producing the observed MIR emission. This requirement arises because in smooth density distribution models, the dust temperature is uniquely determined by the distance from the center. In contrast, the clumpy torus model enables emission from cool dust clouds that are much closer to the center thus able of producing the required MIR emission with much smaller dimensions.
Perhaps the more interesting results are the correlations of the torus parameters with the physical properties of the AGNs, in particular BH mass (), bolometric luminosity and normalized accretion rate, . To study these correlations we have estimated these quantities for all sources using the procedure described in Netzer & Trakhtenbrot (2007). In this procedure (the “virial” mass determination, see Eq. 1 in Netzer & Trakhtenbrot 2007), and FWHM(H) are combined to obtain , and is obtained by using the assumed bolometric correction factor BC. Table 4 lists p-values for Pearson and Spearman-rank correlation tests for the different torus parameters against those properties. These values represent the probability that the observed correlation occurred by chance. We regard as significant all correlations where and indicate with plus and minus signs positive and negative correlations, respectively. Out of the above geometrical parameters, only and show a significant correlation with . This is consistent with a ’receding torus’ assumption (Lawrence 1991) and discussed in §5.4.
5.2. NLR properties
The important parameters of the NLR component are its distance from the center, which determines the dust temperature, and its covering factor, which determines the fraction of MIR flux emitted by the clouds ( in Eq.7).
The photoionization models used to calculate the NLR emission are all of the same constant density () and column density ( cm). As explained earlier, the clouds are optically thin in the MIR range and their distance from the central source uniquely determine the NLR dust temperature. The assumed gas density plays no important role and serves only as a convenient way to determine the photoionization model parameters (in the actual calculations, the model parameters are density and ionization parameter). The dust column density is determined by the assumed gas composition and depletion. The present calculations use solar composition with abundances smaller than assumed in several recent NLR models such as the constant pressure models of Groves et al. (2006). This again is of little importance since the dust column density is the only important parameter, unless an unusual grain mixture is used.
Table 5 lists NLR distances for the 26 best models. Radii are given in pc as well as in units of the sublimation distances, and . For the graphite sublimation radius, the NLR distances span the range with a mean value of 707. This translates to with a mean value of 280. The distribution of in units of the graphite sublimation radius is shown in figure 7.
The mean distance found here is about a factor of 4 larger than the values found by S08. The reason for the larger distances is mostly due to the introduction of real torus models compared with the combination of BBs used in S08. The clumpy torus models contribute a certain amount to the observed silicate emission and hence the NLR contribution, which was the only one with such features considered in S08, is reduced. This can also be seen when comparing the NLR covering factor found here and the one found in S08 (see Table 6). The mean covering factor for the NLR component in S08 is a factor 2 larger than the one found here which is 7%. We consider the present values more appropriate to the objects in question (see §5.4)
The slope is consistent with 0.5, which is the one predicted for simple single-cloud photoionization models and a constant gas.
The dusty NLR dimensions can be compared with earlier works on NLR size such as Bennert et al. (2002), Schmitt et al. (2003), Netzer et al. (2004) and Baskin and Laor (2005). These papers use two different methods to estimate the NLR size. The Bennert et al. and Schmitt et al. works are based on detailed HST imaging of AGN samples that are of the same redshift range as the one considered here. Our dust temperature-based dimensions are smaller by a factor of compared with these imaging-based estimates. This is not surprising since our estimates are biased toward those NLR clouds that have the largest covering factor for the hottest dust clouds. The discrepancy suggests that real NLRs, i.e. those with a range of dust temperatures contribute more to the longer wavelength part of the spectrum. Furthermore, Netzer et al. (2004) suggested that the Bennert et al. (2002) estimates are considerably larger than “real” NLR sizes for two reasons. The first is that the dependence must break down at some intermediate luminosity since extrapolating this relation to high luminosity objects would result in NLR dimensions that exceed the host galaxy dimension. The second reason is that Bennert et al. (2002) measured image sizes that encompasses more than 98% of the detectable [O iii] emission. A 90% or 95% intensity limits can lead to much smaller dimensions.
Baskin & Laor (2005) used a different method to estimate the distances to 40 NLRs in the Boroson & Green (1992) sample. They fitted the observed intensities of the [O iii] 5007, [O iii] 4363 and H narrow emission lines using two different models. The first, assumed a single uniform, [O iii] 5007 emitting region. This resulted in 111In this expression and the following ones, we converted the continuum luminosity given by Baskin & Laor (2005) into bolometric luminosity. and a very large mean gas density of . These dimensions are smaller than those found by Bennert et al. (2002). The second model assumed two emitting regions: a low, constant gas density () region where most of the [O iii] 5007 originates and a high, constant gas density () region where most of the [O iii] 4363 originates. For the low density region they found and for the high density region . The [O iii] 5007 dimensions are comparable to those of Bennert et al. (2002). It is clear that single-zone and even two-zone NLR models are highly simplified. Moreover, dimension deduced from a narrow line intensity is uncertain since the gas density is changing across the NLR and in general . Thus, it is not surprising that the NLR size obtained from such estimates spans a very large range.
The M82 starburst template subtracted from the spectra could, in principle, have an effect on the calculated NLR properties, e.g. a different template could have larger contribution towards shorter wavelengths (30 m). Consequently, the NLR component would have less weight and different spectral shape that would result in larger NLR distances. We regard this possibility as an additional uncertainty on the determination of the NLR distance.
In conclusion, our hot dust measure of the NLR size is based on a method different than the direct imaging method and the photoionization modelling method. While further discussion of those discrepancies is beyond the scope of the present paper, the main conclusion is that simple, dusty NLR clouds (or single NLR shells) at the distances found above, in combination with the two other model components, can adequately fit the 2–35 m spectrum of the QUEST AGNs.
5.3. Hot BB properties
The wavelength range is dominated by the hot BB component of our model (see Fig. 4). On average, the luminosity of this component is about 40% of the total luminosity and is comparable to the luminosity radiated by the clumpy torus. The mean BB temperature is 1400 K and the distribution of temperatures is shown in Fig. 9. While such a large contribution is evident in all of the sources, the physical origin of this component is not yet clear. Clumpy tori of the type considered here cannot produce more than one luminosity bump. Since most of the torus radiation peaks at 5–20 m, the short wavelength excess cannot result from this structure. The conclusion is that the assumed BB must be a separate component.
We have considered several possible explanations for the hot BB component. The first is a contribution from old stellar population near the galactic center. To examine this, we use to estimate the host galaxy luminosity using the scaling relations of Lauer et al. (2007; see their Eq. 3). To convert visual magnitudes to K–band magnitudes we use (V-K)=3.2 typical of giant elliptical galaxies (e.g. Grasdalen 1980). This K–band luminosity was compared to the luminosity of the BB found here. In all cases, the BB component is about a factor 10 more luminous than the derived K–band luminosity of the host galaxy. Thus, such a stellar population would not have sufficient luminosity to account for the observed 2–7 m flux.
The space distribution of such an old stellar population is another limitation. Broad band monitoring of several nearby, intermediate luminosity AGN clearly show K-band variations on time scales of weeks to months (e.g., Suganuma et al. 2006). This has been interpreted as changes in the temperature of the innermost dust, just beyond the sublimation radius. Given the great similarity between the 2–35 m spectrum of the present sample and the spectra of those variable AGNs, we conclude that the origin of the K and the L-band emission, in the present sample and the above AGN sample, is the same. Such dimensions are clearly too small for the assumed stellar population. Needless to say, any K or L-band variations are inconsistent with stellar population properties.
Another strong source of radiation is starburst activity in the host galaxy. The IR emission from starburst, however, is expected to peak at much larger wavelengths. S06 found a strong correlation between the strength of the PAH feature at 7.7 m and the continuum luminosity at 60 m. The same trend is found in starburst dominated ultraluminous infrared galaxies (ULIRGs; see Fig.4 in S08). Thus, luminosity produced by starburst in the QUEST sample dominates the far-IR part of the spectrum and cannot account for the hot dust component.
Another possibility that has already been mentioned is that the new component represents emission by very hot dust, which is not a part of the clumpy torus structure. As explained in §4.3, we could not fit the observations assuming single hot dust clouds of the same grain composition as the torus. The reason is that such clouds produce extremely strong 10 m emission. This is not the case for pure graphite dust that has no prominent features at MIR wavelengths. Graphite dust grains can survive at the BB temperatures found here and provide a suitable explanation for the BB component. The location of this dust can be inside the silicate sublimation radius but still outside the dust-free broad line region. A full discussion of this possibility, including realistic hot graphite grains spectra, is beyond the scope of this paper.
Finally we note that the above suggested location, close or even inside the silicate dust sublimation radius, is also problematic from the point of view of the derived covering factors. Our calculations assume no obscuration of one component by another yet such a location must obscure some of the radiation impinging on the torus. The location of these clouds imply that the radiation incident upon the torus should change and include part of the re-emitted radiation from the BB component. Given the unknown properties of the hot BB component, we cannot take this into account within the framework of the present work.
5.4. Covering factors
A major assumption of this work is that the entire MIR spectrum, after starburst subtraction, is reprocessed AGN radiation. This can be used to deduce the covering factor of the central source by the three components (see also S08 and references therein; Maiolino et al. 2007). Regarding the BB and the NLR components, this factor is obtained directly from the comparison of their deduced luminosities with . The covering factor of the torus, , is different because of its non-isotropic radiation. For a given model, is calculated from the torus parameters. To demonstrate these differences, we compare in Table 6, , , and the covering factors of the NLR and the BB components. The total covering factor can be defined in two ways: 1. Apparent total covering factor, defined as . This is similar to the number used in all earlier investigations of this type. Since the Spitzer-IRS spectra extend only to 35 m, we use the mean “intrinsic” AGN SED of Netzer et al. (2007), obtained from the starburst subtracted QUEST spectra, to calculate the ratio between and . This ratio is 1.072 and thus we multiply of each object by this factor to get its apparent total covering factor. 2. Real total covering factor, calculated by summing together and the NLR and BB covering factors. The various covering factor distributions are shown in Fig. 10.
The mean value of is 0.27. This value is lower than the , for which the mean value is 0.34. The difference is due to the anisotropy of the torus radiation and the fact that, for the present sample of type-I sources, the torus inclination angle is small. The mean covering factor of the NLR component is 0.07. This is in agreement with other estimates that are based on narrow emission line imaging and line equivalent width measurements, which suggest NLR covering factor % (e.g., Baskin & Laor 2005; Netzer & Laor 1993). The mean covering factor of the BB component is 0.23. The mean real total covering factor in our sample is 0.57, which is slightly lower than the mean apparent total covering factor of 0.59.
|Object||Real total CF||Apparent total CF||Real torus CF ()||NLR CF||BB CF|
The dependence of the different covering factors on bolometric luminosity is seen in Fig. 11. Looking at the entire sample, we find no significant correlation of bolometric luminosity with any of the covering factors. However, two of the sources, PG 1448+273 and PG 1244+026, seem to be affecting the entire correlation in a way that their removal from the sample makes the luminosity CF correlations significant. These sources are marked with different symbols on the diagram. The reason why these objects should be omitted from the sample is not clear. We note, however, that these are the sources with the lowest thus may represent different properties. In particular, they are likely to be the highest amplitude variables and the highest accretion rate BHs thus their derived may be wrong. Visual inspection of the bottom left panel shows that, without the two sources, the geometrical covering factor (+BB) decreases with and its highest value, at the low luminosity end, is roughly 0.7, in agreement with typical assumed distributions between type-I and type-II AGN in the local universe.
Several earlier studies found clear anti-correlation between covering factor and (e.g., Wang et al. 2005; Richards et al. 2006; Maiolino et al. 2007; Treister et al. 2008). The decrease in covering factor translates, within the unification scheme, to a decreasing fraction of obscured AGNs as a function of luminosity. All earlier studies were based on the total MIR emission and took no account of the various contributers to the flux and the fact that the geometrical covering factor cannot be obtained, directly, from L(IR)/. Our results are consistent with these studies but take into account, more accurately, the exact torus geometry and the needed translation between observed IR flux and geometrical structure. They show that of the derived covering factor is due to NLR clouds and properly account for the inclination angle distribution within the clumpy torus scenario. Taking the values found here we conclude that for there is roughly 1:1 ratio between type-I and type-II sources (note that the NLR contribution is subtracted from the total covering factor). Both Maiolino et al. (2007) and Treister et al. (2008) find a somewhat larger value for this bolometric luminosity. However, estimates of this ratio based on X-ray surveys are similar to the one found here and consistent across the entire luminosity range in our sample (e.g., Hasinger et al. 2004; Treister & Urry 2006).
The physical mechanism responsible for the decrease of covering factor with is still undetermined. One possibility is a receding torus mentioned above, where higher luminosity implies larger dust sublimation distance, and hence an obscuring structure that is located farther away from the center. In this scenario, the obscuring structure must have a constant height. In the clumpy torus model this scenario corresponds to a decrease of with . Indeed, we find that is anti-correlated with the bolometric luminosity.
We also find a positive correlation between the inclination angle of the torus and the bolometric luminosity. This is also consistent with a receding torus. An obscuring structure located at a larger distance, in this scenario, implies a larger solid angle through which an unobscured view to the center is possible. Thus, the luminosity of the source determines the possible range of inclination angles in which the model would be consistent with a type-I AGN SED. Since there is no preference to a specific angle within that range, the correlation is due to the expansion of the possible range in with increasing luminosity.
5.5. Torus mass
Given the torus parameters, we can estimate the torus mass from
where is the column density of a single cloud in units of cm obtained from and the assumed dust-to-gas ratio. is a function of Y and has the values of 1, and for 1, 2, and 0, respectively.
The torus mass shows a clear correlation with both the bolometric luminosity and the mass of the central black hole (see Fig. 11). These correlations are expected since the bolometric luminosity appears in the torus mass calculation and, in general, . No correlation is found between the torus mass and . We also checked for correlation between and and found none.
The mass of the torus ranges from to for the silicate dust sublimation radius. The ratio between the torus and the BH masses is ranges from to 0.08 and does not show any correlation with the source luminosity.
Assuming a mass-independent BH accretion efficiency of , we can calculate the time it would take for the torus mass to be completely accreted by the BH. This is an extremely short duration of about . This time scale is similar for all sources and is independent of torus mass, source luminosity and any of the other properties. The typical e-folding time for a BH mass of 10 and the above is about . This mean that, if all accretion to the center goes through the torus, there are more than 100 short accretion episodes per one e-folding growth time.
The above result indicates at least two different possibilities: 1. A constant replenishment of the torus material by larger scale mass inflow through the plane of the torus. Such a scenario, requires a constantly changing torus structure. 2. Matter accreting onto the BH does not originate in the torus. Such a scenario has been suggested by Elitzur & Shlosman (2006), who assumed that mass from the host galaxy arrives directly into the accretion disk. Consequently, wind from the accretion disk forms the toroidal obscuration. As long as the large scale mass infall is sufficiently large and continuous, the wind sustains a steady outflowing clumpy structure, which is the clumpy torus discussed in this work. In this case, the small mass of the torus is not directly related to the growth time of the BH.
- Abazajian et al. (2008) Abazajian, K., & Sloan Digital Sky Survey, f. t. 2008, arXiv:0812.0649
- Allen (1976) Allen, D. A. 1976, ApJ, 207, 367
- Antonucci (1993) Antonucci, R. 1993, ARA&A, 31, 473
- Balzano & Weedman (1981) Balzano, V. A., & Weedman, D. W. 1981, ApJ, 243, 756
- Barvainis (1987) Barvainis, R. 1987, ApJ, 320, 537
- Baskin & Laor (2005) Baskin, A., & Laor, A. 2005, MNRAS, 358, 1043
- Bennert et al. (2002) Bennert, N., Falcke, H., Schulz, H., Wilson, A. S., & Wills, B. J. 2002, ApJ, 574, L105
- Bock et al. (2000) Bock, J. J., et al. 2000, AJ, 120, 2904
- Boroson & Green (1992) Boroson, T. A., & Green, R. F. 1992, American Institute of Physics Conference Series, 254, 584
- Cameron et al. (1993) Cameron, M., Storey, J. W. V., Rotaciuc, V., Genzel, R., Verstraete, L., Drapatz, S., Siebenmorgen, R., & Lee, T. J. 1993, ApJ, 419, 136
- Draine (2003) Draine, B. T. 2003, ApJ, 598, 1017
- Edelson & Malkan (1986) Edelson, R. A., & Malkan, M. A. 1986, ApJ, 308, 59
- Efstathiou & Rowan-Robinson (1995) Efstathiou, A., & Rowan-Robinson, M. 1995, MNRAS, 273, 649
- Elitzur & Shlosman (2006) Elitzur, M., & Shlosman, I. 2006, ApJ, 648, L101
- Elvis et al. (1994) Elvis, M., et al. 1994, ApJS, 95, 1
- Gavazzi & Boselli (1996) Gavazzi, G., & Boselli, A. 1996, Astrophysical Letters Communications, 35, 1
- Glikman et al. (2006) Glikman, E., Helfand, D. J., & White, R. L. 2006, ApJ, 640, 579
- Granato & Danese (1994) Granato, G. L., & Danese, L. 1994, MNRAS, 268, 235
- Grasdalen (1980) Grasdalen, G. L. 1980, Objects of High Redshift, 92, 269
- Groves et al. (2006) Groves, B., Dopita, M., & Sutherland, R. 2006, A&A, 458, 405
- Guyon (2002) Guyon, O., 2003, PhD Thesis, Univ. Paris VI
- Guyon et al. (2006) Guyon, O., Sanders, D. B., & Stockton, A. 2006, ApJS, 166, 89
- Hasinger (2004) Hasinger, G. 2004, Nuclear Physics B Proceedings Supplements, 132, 86
- Higdon et al. (2004) Higdon, S. J. U., et al. 2004, PASP, 116, 975
- Houck et al. (2004) Houck, J. R., et al. 2004, Proc. SPIE, 5487, 62
- Hyland & Allen (1982) Hyland, A. R., & Allen, D. A. 1982, MNRAS, 199, 943
- Imanishi et al. (2006) Imanishi, M., Dudley, C. C., & Maloney, P. R. 2006, ApJ, 637, 114
- Jaffe et al. (2004) Jaffe, W., et al. 2004, Nature, 429, 47
- Jarrett et al. (2000) Jarrett, T. H., Chester, T., Cutri, R., Schneider, S., Skrutskie, M., & Huchra, J. P. 2000, AJ, 119, 2498
- Jarrett et al. (2003) Jarrett, T. H., Chester, T., Cutri, R., Schneider, S. E., & Huchra, J. P. 2003, AJ, 125, 525
- Kishimoto et al. (2007) Kishimoto, M., Hönig, S. F., Beckert, T., & Weigelt, G. 2007, A&A, 476, 713
- Krolik & Begelman (1988) Krolik, J. H., & Begelman, M. C. 1988, ApJ, 329, 702
- Lauer et al. (2007) Lauer, T. R., et al. 2007, ApJ, 662, 808
- Lawrence (1991) Lawrence, A. 1991, MNRAS, 252, 586
- Maiolino et al. (2007) Maiolino, R., Shemmer, O., Imanishi, M., Netzer, H., Oliva, E., Lutz, D., & Sturm, E. 2007, A&A, 468, 979
- Marconi et al. (2004) Marconi, A., Risaliti, G., Gilli, R., Hunt, L. K., Maiolino, R., & Salvati, M. 2004, MNRAS, 351, 169
- Mathis et al. (1977) Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425
- Matsuoka et al. (2005) Matsuoka, Y., Oyabu, S., Tsuzuki, Y., Kawara, K., & Yoshii, Y. 2005, PASJ, 57, 563
- McAlary et al. (1983) McAlary, C. W., McLaren, R. A., McGonegal, R. J., & Maza, J. 1983, ApJS, 52, 341
- Minezaki et al. (2004) Minezaki, T., Yoshii, Y., Kobayashi, Y., Enya, K., Suganuma, M., Tomita, H., Aoki, T., & Peterson, B. A. 2004, ApJ, 600, L35
- Nenkova et al. (2002) Nenkova, M., Ivezić, Ž., & Elitzur, M. 2002, ApJ, 570, L9
- Nenkova et al. (2008) Nenkova, M., Sirocky, M. M., Ivezić, Ž., & Elitzur, M. 2008, ApJ, 685, 147
- Nenkova et al. (2008) Nenkova, M., Sirocky, M. M., Nikutta, R., Ivezić, Ž., & Elitzur, M. 2008, ApJ, 685, 160
- Netzer & Laor (1993) Netzer, H., & Laor, A. 1993, ApJ, 404, L51
- Netzer & Trakhtenbrot (2007) Netzer, H., & Trakhtenbrot, B. 2007, ApJ, 654, 754
- Netzer et al. (2004) Netzer, H., Shemmer, O., Maiolino, R., Oliva, E., Croom, S., Corbett, E., & di Fabrizio, L. 2004, ApJ, 614, 558
- Netzer et al. (2007) Netzer, H., et al. 2007, ApJ, 666, 806
- Netzer et al. (2007) Netzer, H., Lira, P., Trakhtenbrot, B., Shemmer, O., & Cury, I. 2007, ApJ, 671, 1256
- Netzer (2006) Netzer, H. 2006, ApJ, 652, L117
- Netzer (2009) Netzer, H. 2009, ApJ, 695, 793
- Neugebauer et al. (1979) Neugebauer, G., Oke, J. B., Becklin, E. E., & Matthews, K. 1979, ApJ, 230, 79
- Neugebauer et al. (1987) Neugebauer, G., Green, R. F., Matthews, K., Schmidt, M., Soifer, B. T., & Bennett, J. 1987, ApJS, 63, 615
- Ossenkopf et al. (1992) Ossenkopf, V., Henning, T., & Mathis, J. S. 1992, A&A, 261, 567
- Packham et al. (2005) Packham, C., Radomski, J. T., Roche, P. F., Aitken, D. K., Perlman, E., Alonso-Herrero, A., Colina, L., & Telesco, C. M. 2005, ApJ, 618, L17
- Pier & Krolik (1992) Pier, E. A., & Krolik, J. H. 1992, ApJ, 401, 99
- Pier & Krolik (1993) Pier, E. A., & Krolik, J. H. 1993, ApJ, 418, 673
- Radomski et al. (2003) Radomski, J. T., Piña, R. K., Packham, C., Telesco, C. M., De Buizer, J. M., Fisher, R. S., & Robinson, A. 2003, ApJ, 587, 117
- Richards et al. (2006) Richards, G. T., et al. 2006, ApJS, 166, 470
- Rieke (1978) Rieke, G. H. 1978, ApJ, 226, 550
- Riffel et al. (2009) Riffel, R. A., Storchi-Bergmann, T., & Mcgregor, P. J. 2009, arXiv:0903.2798
- Rowan-Robinson (1995) Rowan-Robinson, M. 1995, MNRAS, 272, 737
- Rudy et al. (1982) Rudy, R. J., Levan, P. D., & Rodriguez-Espinosa, J. M. 1982, AJ, 87, 598
- Sanders et al. (1989) Sanders, D. B., Phinney, E. S., Neugebauer, G., Soifer, B. T., & Matthews, K. 1989, ApJ, 347, 29
- Schartmann et al. (2005) Schartmann, M., Meisenheimer, K., Camenzind, M., Wolf, S., & Henning, T. 2005, A&A, 437, 861
- Schmidt & Green (1983) Schmidt, M., & Green, R. F. 1983, ApJ, 269, 352
- Schmitt et al. (2003) Schmitt, H. R., Donley, J. L., Antonucci, R. R. J., Hutchings, J. B., Kinney, A. L., & Pringle, J. E. 2003, ApJ, 597, 768
- Schweitzer et al. (2006) Schweitzer, M., et al. 2006, ApJ, 649, 79
- Schweitzer et al. (2008) Schweitzer, M., et al. 2008, ApJ, 679, 101
- Scoville et al. (2000) Scoville, N. Z., et al. 2000, AJ, 119, 991
- Sirocky et al. (2008) Sirocky, M. M., Levenson, N. A., Elitzur, M., Spoon, H. W. W., & Armus, L. 2008, ApJ, 678, 729
- Sitko et al. (1982) Sitko, M. L., Stein, W. A., Zhang, Y.-X., & Wisniewski, W. Z. 1982, ApJ, 259, 486
- Soifer et al. (2003) Soifer, B. T., Bock, J. J., Marsh, K., Neugebauer, G., Matthews, K., Egami, E., & Armus, L. 2003, AJ, 126, 143
- Spinoglio et al. (1995) Spinoglio, L., Malkan, M. A., Rush, B., Carrasco, L., & Recillas-Cruz, E. 1995, ApJ, 453, 616
- Stein & Weedman (1976) Stein, W. A., & Weedman, D. W. 1976, ApJ, 205, 44
- Sturm et al. (2000) Sturm, E., Lutz, D., Tran, D., Feuchtgruber, H., Genzel, R., Kunze, D., Moorwood, A. F. M., & Thornley, M. D. 2000, A&A, 358, 481
- Sturm et al. (2006) Sturm, E., Hasinger, G., Lehmann, I., Mainieri, V., Genzel, R., Lehnert, M. D., Lutz, D., & Tacconi, L. J. 2006, ApJ, 642, 81
- Suganuma et al. (2006) Suganuma, M., et al. 2006, ApJ, 639, 46
- Surace & Sanders (1999) Surace, J. A., & Sanders, D. B. 1999, ApJ, 512, 162
- Surace et al. (2001) Surace, J. A., Sanders, D. B., & Evans, A. S. 2001, AJ, 122, 2791
- Teplitz et al. (2006) Teplitz, H. I., et al. 2006, ApJ, 638, L1
- Thompson et al. (2009) Thompson, G. D., Levenson, N. A., Uddin, S. A., & Sirocky, M. M. 2009, ApJ, 697, 182
- Tomono et al. (2001) Tomono, D., Doi, Y., Usuda, T., & Nishimura, T. 2001, ApJ, 557, 637
- Treister & Urry (2006) Treister, E., & Urry, C. M. 2006, ApJ, 652, L79
- Treister et al. (2008) Treister, E., Krolik, J. H., & Dullemond, C. 2008, ApJ, 679, 140
- Tristram et al. (2007) Tristram, K. R. W., et al. 2007, A&A, 474, 837
- van Bemmel & Dullemond (2003) van Bemmel, I. M., & Dullemond, C. P. 2003, A&A, 404, 1
- Veilleux et al. (2006) Veilleux, S., et al. 2006, ApJ, 643, 707
- Wang et al. (2005) Wang, J.-M., Zhang, E.-P., & Luo, B. 2005, ApJ, 627, L5
- York et al. (2000) York, D. G., et al. 2000,AJ, 120, 1579