Modeling the Strongest Silicate Emission Features of Local Type 1 AGNs

We use the latest version of the AGN catalog of Véron-Cetty & Véron (2010) and the sample of low-luminosity type 1 AGN in Mason et al. (2012) to select AGN that have redshifts z < 0.1 and available 2D low-resolution IRS/Spitzer spectra in the CASSIS database (v6., Lebouteiller et al. 2011). We select those objects with spectra extracted as point-sources that cover the 5–35μm spectral range. The first criterion allows us to probe the MIR emission from the central region of the AGN, and the second one allows us to study both the 10 and 18 μm silicate features. As an exception, we include NGC 3998, a strong silicate emitter (e.g., Mason et al. 2012), which has a low-resolution spectrum in the ∼7.5–14.5 μm range and a high-resolution spectrum in the ∼14–35 μm range. In order to ensure that the emission in the IRS/Spitzer spectra is mostly dominated by dust heated by the AGN, we decompose them into their stellar, ISM, and AGN emission (dust heated by AGN) components at ∼5–15 μm, using the spectral decomposition tool DeblendIRS (Hernán-Caballero et al. 2015). For more details on the spectral decomposition, we refer the reader to Appendix A. We select only those objects (67) for which the spectral contribution due to AGNs is >80% (see Figure 1). Hereafter, we refer to the spectra of these objects as the AGN-dominated IRS/Spitzer spectra. Note that we refer here to the original spectra, not to the component obtained from the decomposition.

Zoom In Zoom Out Reset image size

We obtained the reduced low-resolution (R ∼ 60–127) IRS/Spitzer spectra from the CASSIS database (v6., Lebouteiller et al. 2011). The spectra include the SL1 (λ ∼ 7.4–14.5 μm) and SL2 (λ ∼ 5.2–7.7 μm) modules with a slit-width of 36, and the LL1 (λ ∼ 19.9–39.9μm) and LL2 (λ ∼ 13.9–21.3 μm) modules with a slit-width of 105 (Houck et al. 2004; Werner et al. 2004). We use the final stitched spectra between 5 and 35 μm. The different module spectra were stitched by scaling the LL and SL1 flux modules to the shortest module SL2 flux.

To measure the silicate emission features, we start by interpolating a local continuum between both sides of the emission-line complexes. We inspect each spectrum visually and choose three bands in which to measure the continuum. Band 1 is located at the blue extreme of the 10 μm silicate feature, while band 3 is located at the red extreme of the 18 μm silicate feature, and band 2 is located between the two silicates features. Bands 1, 2 , and 3 are located within the 7.8–8.5, 13.5–14.0, and 20.0–21.0 μm ranges.

In order to better define uncertainties in the continuum definition, which might be affected by the presence of Polycyclic Aromatic Hydrocarbon (PAH) molecular emission lines around 7 μm, and in some cases by the artificial "teardrop"⁸ feature around 14 μm, we trace fiducial mean values of each continuum band by bootstrapping on the measured fluxes (hatched pink regions in Figure 2). We generate 100 continuum values, considering the uncertainties, between bands 1 and 2, and between bands 2 and 3. We randomly associate shorter and longer wavelength mean continuum values to generate linear continua below the silicate features. We note that fitting a spline continuum gives similar results. The dark blue solid lines in Figure 2 are bootstrapped local continua derived from the continuum bands and they define the regions where we measure the strength of the silicate features.

Zoom In Zoom Out Reset image size

Considering that the peak of the silicate features vary from object to object, we choose to measure the strength of the features at the wavelengths where they peak. However, for simplicity, we still call them the 10 and 18 μm silicate features. The silicate strength is, hence, defined as the silicate peak relative to the continuum, at the wavelength where the silicate feature peaks (see, e.g., Hao et al. 2007), according to the equation:

$$\begin{eqnarray}&&{\mathrm{Si}}_{p}=\mathrm{ln}\displaystyle \frac{{f}_{p}(\mathrm{spectrum})}{{f}_{p}(\mathrm{continuum})}.\end{eqnarray} \tag{ 1 }$$

In Table 1, we list the mean and the 68% intervals for the objects with the largest silicate feature strengths. In Appendix B, we list the silicate emission strengths for the full sample. On average, we find that type 1 AGN have a 10 μm silicate strength ${\mathrm{Si}}_{10.3\mu {\rm{m}}}={0.13}_{-0.37}^{+0.15}$ that peaks at ${10.3}_{-0.9}^{+0.7}$ μm, and a 18 μm silicate strength ${\mathrm{Si}}_{17.3\mu {\rm{m}}}={0.14}_{-0.06}^{+0.06}$ that peaks at 17.3 μm.

Table 1.  Basic Properties of Objects in the Si–s Sample and Their Silicate Feature Strengths Measured from the IRS/Spitzer Spectra

Name	Activitya	z	$\mathrm{log}{L}_{X(2-10\mathrm{keV})}$	Ref.	λp	Si${}_{10\mu {\rm{m}}}$	λp	Si${}_{18\mu {\rm{m}}}$
			erg s−1		(μm)		(μm)
NGC 7213	LINER	0.0058	42.2	6	10.7 ± 0.1	0.52 ± 0.05	17.3 ± 0.1	0.23 ± 0.04
PG 2304+042	Sy1.2	0.0420	43.4*	9	10.3 ± 0.1	0.30 ± 0.06	17.3 ± 0.1	0.33 ± 0.08
PKS 0518-45	LINER	0.0420	44.9	8	10.4 ± 0.2	0.31 ± 0.05	17.6 ± 0.1	0.18 ± 0.04
PG 0844+349	Sy1/QSO	0.0640	43.7	1	10.3 ± 0.1	0.33 ± 0.06	16.7 ± 0.1	0.18 ± 0.04
PG 1351+640	Sy1.5/QSO	0.0882	43.1	6	9.8 ± 0.1	0.52 ± 0.03	18.4 ± 0.1	0.14 ± 0.03
PG 2214+139	Sy1.0/QSO	0.0658	43.8	1	10.2 ± 0.1	0.29 ± 0.04	17.8 ± 0.2	0.14 ± 0.04
PG 0804+761	Sy1/QSO	0.1000	44.5	1	9.9 ± 0.1	0.36 ± 0.03	17.3 ± 0.1	0.11 ± 0.04
OQ208	Sy1.5	0.0766	40.8	3	10.2 ± 0.1	0.34 ± 0.05	16.5 ± 0.1	0.12 ± 0.04
NGC 4258	LINERb	0.0015	40.9	5	11.0 ± 0.1	0.29 ± 0.05	17.3 ± 0.1	0.17 ± 0.04
NGC 3998	LINER	0.0035	41.2	4	10.7 ± 0.1	0.36 ± 0.01	16.3 ± 0.1*	0.23 ± 0.02

Note. Columns 1 and 2 list the name and type of activity, column 3 gives the redshift, and columns 4 and 5 list the logarithm of the intrinsic hard X-ray luminosity and the corresponding reference, respectively. Columns 6 and 7 list the wavelength where the $10\,\mu {\rm{m}}$ silicate emission feature peaks, and the $10\,\mu {\rm{m}}$ silicate strength, and columns 8 and 9 are are the wavelength and strength of the 18μm silicate emission feature, respectively. *This value needs to be used carefully, since we are using the high angular resolution spectrum to fix the band 3 and measure the 18μm silicate feature strength. References:^aNED and Véron-Cetty & Véron (2010), ^bMason et al. (2012). *Estimated from (15 to 150)+keV X-ray luminosity assuming a spectral power law with an index α = 1.8. References for hard X-ray luminosity: 1: Zhou & Zhang (2010); 2: Sambruna et al. (2011); 3: Ueda et al. (2005); 4: Younes et al. (2011); 5: Cappi et al. (2006); 6: Bianchi et al. (2009); 7:Brightman & Nandra (2011); 8: Winter et al. (2012); 9: Tueller et al. (2010); 10: Cusumano et al. (2010).

Download table as:  ASCIITypeset image

In most cases, our measurements, which were performed using the AGN-dominated Spitzer spectra, produce similar results to previous works (Hao et al. 2005; Sirocky et al. 2008; Thompson et al. 2009; Mendoza-Castrejón et al. 2015) (see Appendix B). Additionally, we compare our measurements of the silicate strengths with the values obtained using deblendIRS and find that they are similar within the uncertainties (see Appendix A).

We build our final sample by selecting only those objects that show the strongest silicate emission features (see Table 1). From the original sample of 67 local type 1 AGNs, we select those objects with the largest 10 μm silicate strength (${\sigma }_{{\mathrm{Si}}_{10\mu {\rm{m}}}}\gt 0.28$, see Figure 3). The final sample is composed of 10 objects: six Seyfert (Sy) galaxies, and four low-ionization nuclear emission-line region (LINER) galaxies. Four of the Seyfert galaxies are also classified as PG QSOs. Note that we are not selecting the sample according to the type of AGN. We list the type only to give information about their basic properties.

Zoom In Zoom Out Reset image size

Hereafter, we will refer to this sample as the type 1 AGN strong Silicate selected sample (Si–s). The sample spans a range of hard (2–10 keV) X-ray unobscured luminosities between ∼10⁴¹ and ∼10⁴⁵ erg s⁻¹.

This model is represented by a flared disk delimited by the inner and outer torus radius (e.g., Efstathiou & Rowan-Robinson 1995). The inner radius (R_min) is defined by the sublimation temperature of dust grains (1500 K) under the strong radiation produced by the central engine. The size of the torus is determined by the outer radius (Y = R_out/R_in, where Y is a free parameter), and the angular width of the torus Θ (see the cartoon (a) in Figure 4). This model assumes that the dust is composed of graphite grains (53%) with sizes from 0.005 to 0.25 μm, and silicate grains (47%) with sizes from 0.025 to 0.25 μm. The grain sizes are distributed according to Mathis et al. (1977). For the two species, they used the scattering and absorption coefficients given by Laor & Draine (1993). The dust is illuminated by an isotropic central point-like emitting source, which is represented by a broken power law of the form λL_λ ∝ λ^α, with α = 1.2 if 0.001 < λ < 0.03 μm, α = 0 if 0.03 < λ < 0.125 μm, and α = 0.5 if 0.125 < λ < 20.0 μm. Other parameters of the model are the viewing angle i, the polar (γ) and radial (β) indices of the gas density distribution $\rho (r,{\rm{\Theta }})\propto {r}^{\beta }{e}^{-\gamma \times \cos ({\rm{\Theta }})}$ within the torus, and the optical depth τ_{9.7 μm}. For a more complete description of this model, see Fritz et al. (2006).

Zoom In Zoom Out Reset image size

This type of model assumes a central point-like emitting source surrounded by a toroidal distribution of clouds. The emission of the central source is characterized by a broken power law of the form λf_λ ∝ λ^1.2 for λ ≤ 0.01 μm, ∝λ^−0.5 for $0.1\mu {\rm{m}}\leqslant \lambda \leqslant 1\mu {\rm{m}}$, ∝λ⁻³ for λ ≥ 1 μm, and a constant power law between 0.01 and 0.1 μm. Due to its clumpy nature, the central source can directly heat the dust in the inner region of the torus and the dust located at several sublimation radii from the central source. This model assumes standard Galactic spherical dust grains (standard ISM) composed of graphites (47%, Draine 2003) and silicates (53%, Ossenkopf et al. 1992), with a power-law distribution of grain sizes (∝a^−3.5), where a_min = 0.05–0.01μm and a_max = 0.25 μm, respectively. This model has six free parameters: the viewing angle i, the number of clouds along the equatorial ray N₀, the angular width σ, the radial extend Y = R_outer/R_inner, the index of the radial distribution of clouds q, and the optical depth per cloud τ_V (see cartoon (b) in Figure 4). For a more complete description of this model, see Nenkova et al. (2008a, 2008b).

In this case, the toroidal distribution of dusty clumps surrounds a central point-like emitting source, which is described by a broken power law of the form λf_λ ∝ λ for λ < 0.03 μm, ∝λ⁻³ for λ > 3 μm, constant for λ between 0.03 and 0.3 μm, and ∝λ^−4/3 for λ between 0.3 and 3 μm. They modeled the torus following a 3D Monte Carlo radiative transfer simulation. These kinds of simulations fail to properly sample optically thick surface regions with enough grid cells so that each cell is optically thin, which results in underestimating the emission temperature and resulting in a smaller number of model clouds with respect to Nenkova et al. (2008a, 2008b). The dust is composed of three components: a standard ISM component, a standard ISM component with larger grains (a_min = 0.1μm and a_max = 1.0 μm), and intermediate to larger grains (a_min = 0.05μm and a_max = 0.25 μm) with 70% graphite and 30% silicates. The free parameters that describe this model are the viewing angle i, the number of clouds along the equatorial ray N₀, the half-opening angle Θ, the index of the radial distribution a, and the optical depth τ_V (see cartoon (b) in Figure 4). For a more complete description of this model, see Hönig & Kishimoto (2010).

This type of model is based on the parameterization of the clumpy torus model of Hönig & Kishimoto (2010), but instead of adding a blackbody component to take into account the NIR emission, they include a set of different sublimation temperatures for silicate and graphite dust that results in more emission from graphite located in the inner edge of the torus. In this way, when the dust is heated to temperatures >1200 K, smaller silicate grains are destroyed leaving only graphite grains that can be heated up to 1900 K. In the innermost radius only, grains with a size between 0.075 and 1 μm survive.

The clouds are distributed according to a radial power law ∝r^a, where a is the power-law index and r the distance from the black hole in units of the sublimation radius r_sub. They also add a polar outflow, modeled as a hollow cone that can be characterized by the radial distribution of dust clouds in the wind a_w, the half-opening angle of the wind (Θ_w), and the angular width (σ_Θ). Other parameters are the number of clouds along the equatorial ray N₀, and the scale height in the vertical Gaussian distribution of clouds h in the disk, (see cartoon (c) in Figure 4). For a more complete description of this model, see Hönig & Kishimoto (2017).

This kind of models assumes a distribution of high-density dusty clumps embedded in a low-density smooth dusty component. This assumption produces both weaker silicate features and a pronounced NIR emission. They assume that the accretion disk in the nucleus radiates as a broken power law of the form λL_λ ∝ λ^α, where α = 1.2 for a spectral range of 0.001 ≤ λ ≤ 0.01 μm, α = 0 for 0.01 < λ ≤ 0.1 μm, α = –0.5 for 0.1 < λ ≤ 5 μm, and α = –3 for 5 < λ ≤ 50 μm. The dust is distributed following a law that allows a density gradient along the radial (r) and polar (θ) directions, inside a flare disk defined by the inner (R_in), outer radii (R_out), and half-opening angle. The inner radius is defined by the sublimation temperature of 1500 K for an average dust grain size of 0.05 μm. They assumed a standard ISM dust composition with optical properties from Laor & Draine (1993) and Li & Draine (2001). For a more complete description of this model, see Stalevski et al. (2012, 2016).

In this section, we explore how well the dusty torus models reproduce the central wavelength and strength of both 10 and 18 μm silicate features, and the NIR (α_NIR) and MIR (α_MIR) spectral indexes. In order to make a proper comparison, the synthetic and observed central wavelength and strength of both the 10 and 18 μm silicate features are measured following the same methodology described in Section 2.3, and fixing the bands to the sides of the silicate features between 7 and 7.5 μm, 14–15 μm, and 25–26 μm. The synthetic and observed spectral indexes α_NIR and α_MIR are measured between 5.5 and 7.5 μm, and between 7.5 and 14.0 μm, respectively, according to the following definition ${\alpha }_{\mathrm{2,1}}=-\mathrm{log}({f}_{\nu }({\lambda }_{2})/{f}_{\nu }({\lambda }_{1}))/\mathrm{log}({\lambda }_{2}/{\lambda }_{1})$, with λ₂ > λ₁ (see, e.g., Buchanan et al. 2006).

In Figures 5 and 6, we plot the wavelength of the peak of the both 10 and 18 μm silicate features and the strengths of both silicate features as predicted by the models and as observed in the AGN-dominated IRS/Spitzer spectra of the Si–s sample. Additionally, we color-coded the objects according to their bolometric luminosity, which we estimate using the hard X-ray luminosity and the relation derived by Marconi et al. (2004) and Alexander & Hickox (2012).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Since the models are probabilistic in nature, we compare the envelope of measurements performed on the models to the measurements of the Si–s AGN sample. While the Fritz 06, Stalev16, and Nenkova08 models' envelope covers the wavelength space (Figure 5) where our Si–s AGN measurements lie, the Hoenig10 and Hoenig17 models miss one and two objects, respectively, with the lowest value central wavelength of the 18 μm silicate feature. Curiously, these objects are of low bolometric luminosity (OQ 208 and NGC 4258).

For the silicate strengths (Figure 6), we note that the range of synthetic values sampled by the models mostly match the observations. However, we also note that the Hoenig10 and Hoenig17 models show a narrower range of values with respect to the other models. Additionally, the Hoenig17 model never predicts both silicates in absorption. Similar results were reported for a larger sample of AGNs (González-Martín et al. 2019b). Finally, all of the models tend to produce extremely prominent silicate emission features that have not been observed.

In Figure 7, we compare the synthetic and observed silicate strength of the 10 μm silicate feature with the α_NIR and α_MIR spectral indexes. We note that the range of synthetic values sampled by the Fritz06, Hoenig10, and Hoenig17 models mostly matches the observations, independently of the bolometric luminosity. However, the Nenkova08 and Stalev16 models miss several of the NIR and MIR spectral indexes observed. Note that NGC 3998 is excluded from the plots in Figures 5, 6, and 7 because of the short spectral range (∼7.5–14.5 μm) covered by the low-resolution IRS/Spitzer spectrum of this object.

Zoom In Zoom Out Reset image size

We apply the torus models using the computational spectral fitting tool XSPEC, which is part of the HEASOFT⁹ software. These models were recently converted into the XSPEC format in order to fit data in a similar way (see Section 2 in González-Martín et al. 2019a). We also use the set of synthetic stellar and empirical H ii components (ISM) previously converted into XSPEC format (González-Martín et al. 2019a). The former corresponds to a stellar population of 10¹⁰ yr and solar metallicity from the stellar spectral libraries of Bruzual & Charlot (2003), while the empirical H ii components are average starburst templates from Smith et al. (2007).

We model the AGN-dominated IRS/Spitzer spectrum of each object in the Si–s sample, using each one of the four torus family models and the disk+outflow models. Additionally, we add a stellar synthetic and/or H ii empirical component to investigate in which cases adding one or both components really improves the fitting. We also add foreground extinction to the torus models by using the ZDUST component (Pei 1992).

In detail, we start converting the IRS/Spitzer spectra into the XSPEC format, loading the data into XSPEC, removing the parts of the spectra dominated by emission lines, and modeling the spectra assuming the following component combinations:

• 1.  
  AGN dust emission: torus models or disk+outflow model
• 2.  
  AGN + stellar
• 3.  
  AGN + H ii
• 4.  
  AGN + stellar + H ii

For each step, we save the reduced ${\chi }_{\mathrm{red}}^{2}$, the parameters of the model with their uncertainties, and the emission contribution of each component to the total emission of the IRS/Spitzer spectrum between 5.5 and 30 μm. In those cases where ${\chi }_{\mathrm{red}}^{2}\gt 2$, we reported the case as "non-modeled." We perform this procedure for each torus model, resulting in 20 spectral fits per object. Figure 8 shows, for each torus model, the number of spectral fits obtained, using one, two, or three components. Those cases in which none of the component combinations (1–4) are able to model the spectrum with a ${\chi }_{\mathrm{red}}^{2}\lt 2$ are called "non-modeled" spectral fits (see Appendix C). In the next section, we investigate which model best fits the peak and strength of both silicate emission features in the Si–s sample.

Zoom In Zoom Out Reset image size

In most cases, we find that the same object can be modeled assuming one (AGN), two (AGN+stellar or AGN+H ii), or three components (AGN+Stellar+H ii) with a ${\chi }_{\mathrm{red}}^{2}\sim 1$. For each object and torus model, we use the statistical f-test, which allows us to evaluate in which cases the addition of a new component improves the fit from a statistical point of view. We add a new component when the f-test probability is <10⁻⁴. The computational tool xspec includes the ftest command line, which allows us to calculate the f-statistic and its probability, when new and old values of χ² and the degrees of freedom (dof) are provided. As an example, in Figure 9, we show the best fit of the IRS/Spitzer spectrum of PG 0844+349 for each torus model. For this particular object, the Fritz06, Nenkova08, and Stalev16 models need an additional stellar component to fit the spectrum, while the Hoenig10 and Hoenig17 models need also the H ii component. Additionally, for this object, we observe that clumpy and disk+outflow models produce the flattest residuals within the uncertainties and smaller ${\chi }_{\mathrm{red}}^{2}$ (see bottom panel in Figure 9 and Table 3).

Zoom In Zoom Out Reset image size

Table 3.  Best Fits from Smooth Models

			Component Contributions						Component Contributions
Name	${\chi }_{\mathrm{red}}^{2}$	Best Combination	AGN	Stellar	H ii	Name	${\chi }_{\mathrm{red}}^{2}$	Best Combination	AGN	Stellar	H ii
			(%)	(%)	(%)				(%)	(%)	(%)
Smooth torus models of Fritz et al. (2006)						Disk+outflow models of Hönig & Kishimoto (2017)
PG 2304+042	1.92	Torus+Stellar	93	7	0	PG 2304+042	0.57	Disk+outflow+Stellar	93	7	0
PG 0844+349	1.20	Torus+Stellar	93	7	0	PKS 0518-45	0.31	Disk+outflow	100	0	0
PG 1351+640	1.52	Torus+Stellar	97	3	0	PG 0844+349	0.13	Disk+outflow+Stellar+H ii	87	5	8
PG 2214+139	1.52	Torus+Stellar	96	4	0	PG 2214+139	0.20	Disk+outflow+Stellar+H ii	95	1	4
PG 0804+761	0.63	Torus+Stellar	94	6	0	PG 0804+761	0.49	Disk+outflow+H ii	95	0	5
NGC 3998	1.91	Torus+Stellar	100	0	0	OQ 208	1.78	Disk+outflow+Stellar	98	2	0
						NGC 4258	1.88	Disk+outflow+H ii	62	0	38
	.	.	.	.	.	NGC 3998	0.39	Disk+outflow	100	0	0
Clumpy torus models of Nenkova et al. (2008a, 2008b)						Two-phase media torus models of Stalevski et al. (2016)
NGC 7213	1.11	Torus+Stellar	96	4	0	PG 2304+042	1.99	Torus+Stellar	93	7	0
PG 2304+042	1.25	Torus+Stellar	92	8	0	PKS 0518-45	1.65	Torus	100	0	0
PKS 0518-45	1.29	Torus+Stellar	97	3	0	PG 0844+349	0.81	Torus+Stellar	92	8	0
PG 0844+349	0.32	Torus+Stellar	91	9	0	PG 2214+139	0.91	Torus+Stellar+H ii	82	10	8
PG 1351+640	1.92	Torus+Stellar	99	1	0	PG 0804+761	1.11	Torus+Stellar	90	10	0
PG 2214+139	1.92	Torus+Stellar	87	13	0
PG 0804+761	1.60	Torus+Stellar	89	11	0	.	.	.	.	.	.
OQ 208	0.32	Torus+Stellar+H ii	91	3	6	.	.	.	.	.	.
NGC 4258	0.78	Torus+Stellar+H ii	66	4	30	.	.	.	.	.	.
NGC 3998	0.32	Torus+Stellar	81	19	0	.	.	.	.	.	.
Clumpy torus models of Hönig & Kishimoto (2010)						.	.	.	.	.	.
PG 2304+042	1.18	Torus	100	0	0	.	.	.	.	.	.
PKS 0518-45	0.46	Torus+Stellar+H ii	89	6	5	.	.	.	.	.	.
PG 0844+349	0.20	Torus+Stellar+H ii	80	10	10	.	.	.	.	.	.
PG 2214+139	0.71	Torus+Stellar+H ii	76	9	14	.	.	.	.	.	.
PG 0804+761	1.79	Torus+Stellar+H ii	76	9	15	.	.	.	.	.	.
NGC 4258	1.64	Torus+Stellar+H ii	67	3	29	.	.	.	.	.	.
NGC 3998	0.40	Torus	100	0	0	.	.	.	.	.	.

Note. Column 1 lists the name of the object, column 2 the ${\chi }_{\mathrm{red}}^{2}$, column 3 the combination of components that best fit the IRS/Spitzer spectra, and columns 4–6 the percentage contribution of each component.

Download table as:  ASCIITypeset image

In Figure 10, we plot the histogram of the fractional contribution of each spectral component. The AGN component dominates the emission in all cases, which is expected due to our selection of AGN-dominated IRS/Spitzer spectrum sources.

Zoom In Zoom Out Reset image size

Fritz06 models are able to reproduce the IRS/Spitzer spectra in seven of ten objects. The AGN component contributes more than 93%. The Nenkova08 models reproduce the IRS/Spitzer spectra of ten objects, in most cases with an AGN contribution >81%, except for NGC 4258, for which the AGN component contributes 66%. For the Hoenig10 models, we find that four of the objects can be modeled with an AGN contribution >80%, and three other (PG 2214+139, PG 0804+761, and NGC 4258) with a contribution between 67 and 76%. The Hoenig17 models are able to fit seven of the objects with an AGN contribution >87%, and NGC 4258 with an AGN component that contributes >62%. Finally, in the case of two-phase medium models of Stalevski et al. (2016) only five of the ten objects can be reproduced with a contribution of the AGN > 82%.

Irrespective of the model used, NGC 4258 always needs a large (∼20%–38%) contribution from the H ii component. It is possible that in this case, the spectral decomposition tool DeblendIRS is not perfectly separating the different spectral contributions at MIR. We also note that in all cases, the stellar component is necessary to take into account the emission in the bluer extreme of the spectrum, while the H ii component is necessary to take into account the emission in the redder extreme of the spectrum. A similar result was also found by González-Martín et al. (2019a) for a large sample of AGNs.

We obtain upper (or lower) limits of the resulting free parameters and covering factors. The covering factor is defined as 1–P_esc, where P_esc is the probability that a photon emitted in the central engine is able to escape without being absorbed by the torus (in the case of the smooth torus and two-phase medium torus models) or by a dusty cloud in clumpy torus models (Nenkova08, Hoenig10). In the case of the Hoenig17 models, the covering factor derived is the sum of the geometrical covering factors of the disk and the covering factor of the outflow (see González-Martín et al. 2019b).

To obtain well-constrained parameters, a detailed modeling that includes NIR and far-infrared data is necessary (see, e.g., Ramos Almeida et al. 2014). However, simultaneously modeling the NIR and MIR components of the IRS/Spitzer spectrum of type 1 AGNs has been a challenge (e.g., Mor et al. 2009; Hernán-Caballero et al. 2015; Martínez-Paredes et al. 2017). In this work, we assume the NIR as a stellar component. Our purpose is to show if any of the proposed models are able to explain both the peak and shape of the strongest silicate emission features. For this purpose, we only need to check that the range of values of the covering factors obtained from modeling the IRS/Spitzer spectra with each model is within the range of values expected for type 1 AGNs (see, e.g., Alonso-Herrero et al. 2011; Ramos Almeida et al. 2011; Ichikawa et al. 2015; Mateos et al. 2016; González-Martín et al. 2017; Martínez-Paredes et al. 2017).

We note that on average, the smooth models produce a dusty torus with small angular widths, with low and high viewing angles; although, the angular width is poorly constrained. The dusty clumpy torus models of Nenkova08 produce both large and low viewing angles and a range of angular widths from low (15°) to high (70°) values, and a number of clouds along the equatorial ray that are in general ≲7 clouds, resulting in escape probabilities ≳40%. The clumpy models of Hoenig10 produce values of the viewing angle that range from nearly 30° to 80°, angular widths around 55°, and number of clouds in the range 2.5–10.0. The viewing angles in disk+outflow models of Hoenig17 that we find are between 0° and 50°, and the angular widths are between 30° and 45°. Although in three cases, we obtain lower limits of the angular width, indicating that this could be larger. The two-phase models of Stalev16 produce on average large viewing angles (∼80), only in one case we obtain a viewing angle around 10°, and the angular width is around 80°, resulting in very obscured AGNs (see Tables 6–10, and Figure 17 in Appendix C).

In general, Fritz06 and Hoenig17 models produce lower covering factors than the Nenkova08, Hoenig10, and Stalevski16 models. The Fritz06 and Hoenig17 models produce covering factors around 0.2, although they range from 0 to 1. The Nenkova08 model produces covering factors around 0.6 with a range from 0.3 to 1.0. These values are consistent with the range of values obtained by Martínez-Paredes et al. (2017) for a sample of PG QSOs using Nenkova08 models. Feltre et al. (2012) compare smooth and clumpy torus models of Nenkova08 and found that both torus model families produce similar MIR continuum shapes for different model parameters. In the Hoenig10 and Stalevski16 models, the covering factors are large (around 0.8), probably due to the fact that, in these models, the angular width of the torus tends to be larger, which leads to a more obscured AGN.

Using the spectral residuals for all modeled and non-modeled cases, we calculate the average spectral residuals for each torus model. The vertical black solid lines in Figure 11 represents the mean wavelength where the silicate features peak in the IRS/Spitzer spectra, and the black dotted lines their 1σ intervals.

Zoom In Zoom Out Reset image size

In order to discuss qualitatively the similarities and differences between the residuals of the models, we divide the spectral range into three parts in Figure 11. These parts are the region bluewards of 10 μm, between 10 and 18 μm, and redwards 18 μm. In part (a), we observe that on average, around 5 μm, the Fritz06, Nenkova08, Hoenig10, and Stalev16 models are the worst at reproducing the bluer extreme of the spectra within the uncertainties. However, at longer wavelengths, the Fritz06 and Stalevski16 models are the least accurate at reproducing the shape of the spectra. The Hoenig17 models best reproduce the spectra at all wavelengths within the uncertainties. In part (b), Nenkova08 show the flattest residual, while in the redder extreme, the Fritz06 and Stalevski16 models show the largest residuals, while the Hoenig10 and Hoenig17 models shows flatter residuals. In part (c), all models show similar residuals, although Nenkova08 shows the flattest residual within the uncertainties.

We also observe that the Fritz06 and Stalevski16 models underestimate the strength of the 10 and 18 μm silicate features, while the Nenkova08, Hoenig10, and Hoenig17 models best reproduce the peak and the shape of both features. In general, we note that the Hoenig17 and Nenkova08 models show flatter residuals, resulting in the models that best reproduce the shape and peak of the strong silicate features observed in these objects.

Some objects deserve particular attention. For instance, in Figure 19, we can see that for NGC 7213, we require the stellar component to fit the bluer extreme of the spectrum but that both silicate peaks are still underestimated by the Nenkova08 models.

PG 1351+640 is modeled only by the smooth and clumpy models of Nenkova08. In both cases it is necessary to add the stellar component, which contributes 3% and 1% in the case of smooth (Fritz06) and clumpy models (Nenkova08 and Hoenig10), respectively. OQ208 is modeled only by the Nenkova08 and Hoenig17 models. In this case, the residuals from the Nenkova08 model are flatter than the residuals from the Hoenig17 models, although in the case of the Nenkova08 model, it is necessary to include the H ii component with a contribution of 6%, in addition to the stellar component (3%) in order to fit the redder extreme of the spectrum. In contrast, the Hoenig17 model is able to reproduce the entire spectral range requiring only a small contribution from the stellar component (2%). The remaining objects are modeled successfully by the Nenkova08, Hoenig10, and Hoenig17 models.

We divide the objects in the Si–s sample into three groups according to their bolometric luminosities and BH accretion rates. For each group and torus model, we combine the residual obtained from fitting the AGN-dominated Spitzer spectrum with the components C1 (AGN), C2 (AGN+Stellar), C3 (AGN+H ii), and C4 (AGN+Stellar+H ii). In Figure 12, we plot the average residual for each group. We note that for the first group, all models are unable to produce a flat residual around the 10 μm silicate feature. But, for the second group, which covers a larger range of bolometric luminosities and BH accretion rates, all of the residuals become flatter. At the largest bolometric luminosities of the third group, all models show the flattest residuals. These results show that all models fail in reproducing the central wavelength of the 10 μm silicate feature in the objects with lower bolometric luminosities, as we see in Figure 5 in Section 4.1. Additionally, we note that the Hoenig17 models always produce the flattest residual around 5 μm for low, intermediate, and high luminosities.

Zoom In Zoom Out Reset image size

Based on our analysis, the clumpy models of Nenkova08 (${\chi }_{\mathrm{red}}^{2}\sim 1.03$), Hoenig10 (${\chi }_{\mathrm{red}}^{2}\sim 0.89$), and Hoenig17 (${\chi }_{\mathrm{red}}^{2}\sim 0.70$) produce a better fitting of the IRS/Spitzer spectra than the smooth (${\chi }_{\mathrm{red}}^{2}\sim 1.36$) and two-phase medium (${\chi }_{\mathrm{red}}^{2}\sim 1.20$) models. This result is in agreement with the previous evidence that the surrounding medium around the AGN should be clumpy (see Ramos Almeida & Ricci 2017, and references therein). Indeed, Mendoza-Castrejón et al. (2015) found that isolated type 1 AGNs have a clumpy dust distribution, while interacting type 1 AGNs can have both clumpy or smooth dusty distributions. Additionally, we find that all models produce flatter residuals in high-luminosity AGNs than in low-luminosity AGNs, probably due to all models better sampling the central wavelength and both 10 and 18 μm silicate features of high-luminosity AGNs. However, high-luminosity Hoenig17 models produce the flattest residuals along the entire spectral range.

Mor et al. (2009) used IRS/Spitzer spectra of a sample of 26 nearby QSOs to constrain the clumpy models of Nenkova08. They argue that in order to model the entire spectral range between 5 and 30 μm, they need to add two more components, one that takes into account the emission produced by the dust in the narrow-line region (Schweitzer et al. 2008), and another one that takes into account the emission produced by dust close to the AGN, not directly related with the dusty torus (e.g., Minezaki et al. 2004; Kishimoto et al. 2007; Riffel et al. 2009). However, Alonso-Herrero et al. (2011) showed, for a sample of local Seyferts, that when high angular resolution data at NIR and MIR is used, it is not necessary to add any additional components in order to constrain the clumpy models of Nenkova08; although, for some type 1 Seyferts galaxies, they found that the NIR emission is underpredicted by the fitted SED. In a previous work, we used the starburst-subtracted IRS/Spitzer spectra from ∼7.5 to 15μm plus NIR high angular resolution data from the Near Infrared Camera and Multi-Object Spectrometer on Hubble Space Telescope for a sample of 20 nearby QSOs in order to constrain the clumpy torus models of Nenkova08. Martínez-Paredes et al. (2017) found that including the spectral range between 5 and 8μm resulted in a poor fitting of the 10 μm silicate emission feature. In this work, we find that the AGN-dominated Spitzer spectra of the Si–s type 1 AGN can be fitted by Nenkova08 and/or Hoenig10 models by adding a stellar component that takes into account the bluer spectral range, and in some cases, the H ii component, in order to improve the fitting in the redder spectral range. A similar result was also found for a large sample of AGNs (González-Martín et al. 2019a, 2019b).

Hernán-Caballero et al. (2017) found, for a sample of 85 QSOs, that a superposition of two blackbodies between 1.7 and 8.4 μm, with temperatures around 1000 K for the hot blackbody component, and 400 K for the warm blackbody component, can fit the Spitzer spectra between 0.1 and 10 μm. They argue that an additional hotter component of dust is necessary to reproduce the excess emission at 1–2 μm. On the other hand, Lyu et al. (2017) argues that the strong silicate emission features observed in dust-deficient PG QSOs can be explained assuming a reduced height scale of the warm dust, i.e., the dust between the sublimation zone and the outer region of cold dust, allowing for less interception of the radiation from the accretion disk with the inner dust, resulting in a decrease of the MIR continuum, but keeping the heating of the outer dust responsible for the silicate emission features invariant (see Figure 20 in Lyu et al. 2017).

In general, the Nenkova08 models produce a better fit than the Hoenig10 ones. They assume the same geometry, but the former uses a standard dust composition, while the latter includes standard ISM dust plus standard ISM with large grains and also a composition of dust mostly dominated by graphite. The Hoenig10 models assume this dust composition in order to take into account the observational suggestion that the dust composition in AGNs deviates from standard ISM dust (Suganuma et al. 2006; Kishimoto et al. 2007). One of the clear improvements in the Hoenig17 models is that they allow for the existence of different dust compositions in different parts of the surrounding dusty structure (disk+outflow).

The chemical composition of the dust has been largely studied for the ISM (see Henning 2010, and references therein), and to a lesser extent in AGNs (see Lyu et al. 2014, and references therein). However, determining the exact chemical composition of the dust has been challenging. Srinivasan et al. (2017) studied the dust composition of a large sample of PG QSOs with a redshift z < 0.5, which showed the 10 μm silicate feature in emission. They found that the dust is mostly composed of amorphous oxides and silicates, plus a small fraction in crystalline form. However, this small fraction is nearly four times larger than the last upper limit (≲2.2) reported by Kemper et al. (2005) for the ISM, and more similar to the upper limit (≲5) previously reported by Li & Draine (2001) for the ISM.

In general, we find that neither the peak nor the shape of the silicate features of the AGN-dominated IRS/Spitzer spectra of the Si–s type 1 AGN are perfectly reproduced by the models. However, we note that for each object, either the Hoenig17 or the Nenkova08 model produces the smallest ${\chi }_{\mathrm{red}}^{2}$ and flattest residuals. Hönig & Kishimoto (2017) point out that their dust composition explains both the observed small NIR reverberation mapping and interferometric sizes with dust sublimation physics. They argue that this combination of disk+outflow clumpy models is able to reproduce the 3–5μm bump in type 1 AGNs, and preserve the MIR bump produced by the wind. Indeed, according to this model, a standard ISM composition of the dust in the wind would be responsible for the emission of the silicate features. García-González et al. (2017) found that the Hoenig17 (disk+outflow) models predict MIR slopes (between 8.1 and 12.5 μm) and silicate strengths at 10 μm that are in agreement with the values observed in type 1 AGNs. Particularly, they noted that when clouds are more concentrated toward the inner region of the torus (see, e.g., Hönig & Kishimoto 2010; Ramos Almeida et al. 2011; Ichikawa et al. 2015; Martínez-Paredes et al. 2017), the MIR spectral indices are flatter, and the silicate features are stronger than those observed in Seyferts and QSOs.