Hints on the Gradual Resizing of the Torus in AGNs through Decomposition of Spitzer/IRS Spectra

We have only included spectra observed with both the short-low (SL) and long-low (LL) modules to guarantee the full Spitzer/IRS coverage ($\sim 5\mbox{--}30\,\mu {\rm{m}}$). CASSIS and SINGS provide flux- and wavelength-calibrated spectra. However, the observations using data from both the SL and LL spectral modules suffer from mismatches due to telescope pointing inaccuracies or due to different spatial resolutions of the IRS orders. This is not corrected in the final products given by CASSIS and SINGS. We therefore scaled each spectra to the immediate prior (in wavelength range) to overcome such effects. Thus, our flux level is scaled to the level of the shortest wavelengths, which is the order with the highest spatial resolution (3.6 arcsec). This guarantees that the flux level is scaled to the best spatial resolution that Spitzer can provide. Note, however, that this does not solve the problem related to the fact that each IRS module might be seeing a different region due to a different spatial resolution and/or slight changes on the position of the spectrum. However, this is the best that we can do. Finally, the spectra are shifted to the rest frame according to the redshifts of the objects (see Col. 3 in Tables 1–3).

We have included in our analysis mid-infrared spatially resolved images taken with CanariCam/GTC using the filter "Si6" centered at 11.5 μm. These observations are part of proprietary data of a sample of faint and Compton-thick LINERs observed with CanariCam/GTC (proposal ID GTC10-14A, P.I. González-Martín). Some of these images were already used by González-Martín et al. (2015) to compare their nuclear fluxes with Spitzer/IRS fluxes. Note that here we make a more sophisticated treatment of the images to better isolate the nuclear emission from its extended emission (see Section 4). The full sample contains 19 LINERs and it will be the subject of subsequent publications focused on the nuclear (J. Masegosa et al. 2017, in preparation; see Masegosa et al. 2013 for preliminary results on the nuclear flux) and extended emission (as an example, see the analysis of the extended emission of NGC 835 presented in González-Martín et al. 2016). Here we have used the nuclear flux taken from the CanariCam data in 12 objects in common with the Spitzer/IRS sample of LINERs. The summary of the observations used in this paper is reported in Table 4.

Table 4.  CanariCam Observations

	Target Observations				Standard Star Observations				Results
Name	Date	ObsID	Expt.	Config.	Name	ObsID	Expt.	${\mathrm{FWHM}}_{\mathrm{PSF}}$	${\mathrm{FWHM}}_{{\rm{N}}}$	${\mathrm{Flux}}_{{\rm{T}}}$	${\mathrm{Flux}}_{{\rm{N}}}$
	(Y-M-D)		(s)		(HD)		(s)	(arcsec)	(arcsec)	(mJy)	(mJy)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)
NGC 315	2014 Sep 22	1048	927	10/10/90/-180	4502	1045	66	0.34-0.44	0.24–0.29	28 ± 5	9 ± 2
	2014 Sep 22	1050	927	10/10/90/-180	4502	1045	66	0.34–0.44
NGC 835	2014 Sep 24	1281	993	10/10/90/-180	11353	1287	66	0.31–0.36	0.25-0.39	32 ± 5	27 ± 4
	2014 Sep 24	1283	993	10/10/90/-180	11353	1287	66	0.31–0.36
NGC 2685*	2013 Jan 04	7636	661	10/10/90/-180	73108	7640	83	0.30–0.32	⋯	⋯	⋯
	2013 Jan 04	7638	661	10/10/90/-180	73108	7640	83	0.30–0.32
NGC 2655	2013 Jan 04	7642	661	10/10/90/-180	73108	7640	83	0.30–0.32	0.26-0.27	29 ± 5	8 ± 1
	2013 Jan 04	7644	661	10/10/90/-180	73108	7640	83	0.30–0.32
UGC 05101	2013 Jan 01	7270	617	10/10/90/-180	86378	7268	77	0.34–0.43	0.28-0.37	96 ± 15	11 ± 2
NGC 4321	2015 Apr-03	2700	927	10/10/90/-180	108381	2698	66	0.24–0.36	0.26-0.26	15 ± 3	3.0 ± 0.6
	2015 Apr-04	2702	927	10/10/90/-180	108381	2698	66	0.24–0.36
NGC 4486	2015 Apr-06	3117	927	16/16/0/-180	108985	3113	66	0.36–0.42	0.28-0.33	14 ± 2	6 ± 1
	2015 Apr-06	3119	927	16/16/0/-180	108985	3113	66	0.36–0.42
MRK 266NE	2013 Jan 01	7276	617	10/10/90/-180	120933	7272	77	0.32–0.35	0.37-0.46	34 ± 5	16 ± 3
MRK 266SW	2013 Jan 01	7276	617	10/10/90/-180	120933	7272	77	0.32–0.35	0.31-0.37	94 ± 15	11 ± 2
UGC 08696	2013 Jan 01	7274	308	10/10/90/-180	120933	7272	77	0.33–0.35	0.48-0.57	61 ± 10	25 ± 4
NGC 6251	2012 Sep 25	5990	617	12/12/0/-180	144204	5988	77	0.40–0.52	0.25-0.38	18 ± 3	7 ± 1
	2013 Aug 24	3063	662	16/16/0/0	144204	3061	66	0.25–0.26
IRAS 17208-0014	2013 Jul 19	1897	596	10/10/90/-180	153210	1901	66	0.31–0.36	0.28-0.42	87 ± 14	4 ± 1
	2013 Jul 19	1899	596	10/10/90/-180	153210	1901	66	0.31–0.36

Note. The object marked with an asterisk (NGC 2685) was not detected with CanariCam. The column called "Config." corresponds to the nod throw (keyword NODTHROW, in units of arcseconds), chop throw (keyword CHPTHROW, in units of arcseconds), instrument position angle (keyword INSTRPA, in degrees), and chop position angle (keyword CHPPA, in degrees; written as NODTHROW/CHPTHROW/INSTRPA/CHPPA).

Download table as:  ASCIITypeset image

CanariCam uses a Raytheon $320\times 240$ pixel Si:As detector that covers a field of view (FOV) of 26 × 19 arcsec on the sky with a pixel scale of 0.0798 arcsec. Standard mid-infrared chopping–nodding techniques were used to remove the time-variable sky background, the thermal emission from the telescope, and the detector 1/f noise. The employed chopping and nodding throws, and chop and nod position angles are reported in Table 4 (Col. 5).

Images of point-spread function (PSF) standard stars were obtained in the same filter immediately after the science target to accurately sample the image quality and allow for flux calibration of the target observations. Table 4 includes the name (Col. 6), integrating time (Col. 8), and the FWHM of the standard stars associated with each target (Col. 9, representing the FWHM of the PSF at the time of the observations). To compute it, we have fitted a 2D Gaussian to the standard star observations. The two numbers given in Col. 9 in Table 4 show the minor and major widths of this Gaussian fit. Any point-like source detected in our images should show an FWHM contained between these two values. However, as will be discussed below, the sky at mid-infrared is highly variable and often the conditions of the sky at the time the target was observed might have slightly changed from those when the standard star was observed. Thus, although this number can give a rough estimate of the image quality, it cannot be taken as a strict limit to the FWHM of a point-like source.

Each observing block was processed using the pipeline RedCan (González-Martín et al. 2013), which is able to produce flux-calibrated images and wavelength- and flux-calibrated spectra for CanariCam/GTC and T-ReCS/Gemini low-resolution data. The combination of the different observing blocks for the same source (when available; see Table 4) were made after flux calibration using Python routines.

We have compiled all the high spatial resolution mid-infrared images associated with our Spitzer/IRS sample. For that purpose we have used the atlas of mid-infrared observations reported in Asmus et al. (2014). It contains 895 observations of 253 AGNs taken with 8 m class telescopes up to 2014.

This complements our study in three ways: (1) mid-infrared observations centered at wavelengths other than the 11.5 μm of the current CanariCam/GTC-observed LINERs, (2) mid-infrared observations of other LINERs not observed in our campaign with CanariCam/GTC, and (3) mid-infrared observations of objects included in our comparison samples. We have retrieved 285 mid-infrared observations for 18 LINERs, 10 type-1 Seyferts, 23 type-2 Seyferts, and two Starbursts contained in the atlas by Asmus et al. (2014). Among the LINERs, only five of them (UGC 05101, NGC 4486, MRK 266SW, MRK 266NE, and NGC 6251) are in common with the CanariCam/GTC sample, although not exactly with the same filter because the CanariCam/GTC sample was selected, avoiding objects already observed with 8 m class telescopes.

Some of the high spatial resolution images in our sample show extended structures together with the nuclear point-like source (Asmus et al. 2015). In order to better isolate the nuclear component, we have developed a code able to decompose both emissions. This procedure is based on the idea that both emissions can be roughly fitted with 2D Gaussians. We have used a 2D Gaussian fit included in the package satrapy within Python, which allows the normalization of the Gaussian, the width along the major and minor axes, and the angle in which the major axis is located to be varied. First, we trim the image in a box of $40\times 40$ pixels, centered at the position of the source, which is wide enough to contain all of the extended emission for all the objects in our sample but sufficiently small to guarantee that the procedure avoids any artificial or real structure away from our target. Then, we followed several steps until the results converge to the final solution. (1) We fit the image using a single 2D Gaussian. At this stage, the width of the Gaussian is fixed to the width of the standard star associated with the target for CanariCam data (FWHM reported in Col. 9 of Table 4) or with the major axis of the FWHM of the Gaussian fit reported by Asmus et al. (2014). This is considered as an initial guess for the nuclear component. (2) The fitted Gaussian is subtracted from the original images, producing a first guess of the extended structure. This extended structure is fitted with another Gaussian, which is centered at the position of the first Gaussian but now allowing the widths of the Gaussian to vary. (3) This Gaussian fit to the extended emission is now subtracted from the original image, producing a new guess for the nuclear component. At this stage, the process starts over in (1), using this new guess for the nuclear component as the input image. This process continues until the residuals are within three standard deviations over the background of the image.

As an example of the result of this process, Figure 2 shows the case for Mrk 266SW. Panel (a) shows the original image where both point-like and extended emission can be clearly spotted. Panels (b) and (c) correspond to the best Gaussian fit to the extended and nuclear emissions, respectively. Panels (d) and (e) show the resulting "extended" and "nuclear" images, respectively, computed as the original image minus the Gaussian best fit for the nuclear component, and the original minus the Gaussian best fit for the extended emission. Panel (f) (inset in panel (a)) shows the residuals of the final fit. This method nicely isolates the nuclear from extended emission.

Zoom In Zoom Out Reset image size

The resulting minor and major axes of the FWHM are recorded in Col. 10 in Table 4 for the CanariCam images and in Col. 3 in Table 5 for other archival data. In the case of CanariCam data, we can compare the final FHWM of the nuclear component with that of the PSF (as traced by the standard star). Although in most cases the FWHM of the nuclear component is consistent with that of the PSF, it is clear that in some cases the point-like source detected in the target image shows a smaller or larger FWHM than that of the standard star (e.g., UGC 08696 or NGC 315). When the FWHM of the nuclear component is larger than the PSF of the standard star, it is plausible that the source is partially resolved. However, it could also be due to changes on the conditions of the observation throughout the night. Due to the faintness of our targets, the exposure times are rather long. Thus, a delay of one hour or more between the observation of the target and that of the standard star is usual.

Table 5.  Sample of Archival High-resolution Mid-IR Imaging Results

Name	Wavelength	${\mathrm{FWHM}}_{{\rm{N}}}$	${\mathrm{Flux}}_{{\rm{N}}}$
	($\mu {\rm{m}}$)	(arcsec)	(mJy)
NGC 1052	10.6	0.5-0.4	133 ± 24
	11.5	0.4–0.5	144 ± 26
	12.5	0.6–0.5	172 ± 31
	9.8	0.4–0.4	112 ± 20
	7.8	0.6–0.4	65 ± 11
	8.6	0.4–0.5	62 ± 11
	11.9	0.4–0.4	125 ± 21
	18.7	0.5–0.5	307 ± 51
	18.3	0.6–0.5	315 ± 53
	8.7	0.8–0.6	63 ± 10
NGC 1097	12.3	0.4–0.4	24.9 ± 4.2
	8.6	0.3–0.6	25.1 ± 4.2
	11.2	0.6–0.7	25.7 ± 4.5
	11.2	0.4–0.5	26.6 ± 4.4
	11.9	0.5–0.4	22.1 ± 3.9
	11.9	0.3–0.4	21.8 ± 3.7
	18.7	0.5–0.5	47.3 ± 7.9
	18.3	0.5–0.6	45.1 ± 7.7
	10.8	0.4–0.5	24.5 ± 4.1
	11.7	0.5–0.5	23.3 ± 4.0

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

The nuclear fluxes (reported in Col. 11 of Table 4 and in Col. 4 of Table 5) are computed by performing aperture photometry in the nuclear images (i.e., after subtracting the extended emission), using an aperture radius 2.5 times the major width of the 2D Gaussian fit for the nuclear component. This ensures that over 97% of the flux is contained within this aperture and, at the same time, avoids any residuals of the extended emission contributing to the nuclear flux. We computed the error as the quadratic sum of the flux calibration uncertainty plus the error due to the S/N of each observation, but the errors are fully dominated by the former one. They are assumed to be 15% of the flux. Note that calibration errors need to be included especially when combining data sets from different facilities. Alonso-Herrero et al. (2016) computed the flux calibration in a sample of CanariCam/GTC observations of AGNs, finding that the mean calibration error is 11%, close to our estimate (see also Díaz-Santos et al. 2010; Ramos Almeida et al. 2011). We also compared the nuclear fluxes reported here with those reported by a Gaussian fit in Asmus et al. (2014). All of them are consistent within the errors. However, our errors are larger than those reported by Asmus et al. (2014). Note that our final purpose is to constrain the torus component when decomposing the Spitzer/IRS spectra. A less restrictive limit could translate into a larger contribution for this component. Since we are studying the plausible disappearance of the torus, this less restrictive constraint yields to a more conservative result on the disappearance of the torus (see Section 4.2).

We used the model-independent spectral decomposition called DeblendIRS presented by HC15 to decompose the Spitzer/IRS spectra. This code uses a set of IRS spectra as templates for purely stellar, interstellar (dominated by polycyclic aromatic hydrocarbons, PAHs), and AGN-dominated components. Here in after, we call these three components stellar, ISM, and torus, respectively. Note that we refer to the AGN component as the torus component since the torus is the dominant source of AGN continuum in the mid-infrared. The algorithm computes the marginalized probability distribution of physical parameters, from comparison of the observed data with all the models in the library, using Bayesian inference.

We initially have used as templates the same library presented by HC15. However, in their work, they focused on the decomposition in the range between 5 and $15\,\mu {\rm{m}}$. For that reason, some of the templates did not cover the full Spitzer/IRS wavelength range. Our purpose is to decompose the full range covered by IRS spectra as much as possible in order to study the entire Spitzer spectra (and to compare with clumpy models). For that reason, we have removed the templates with redshift above 0.2. All together, we have removed 76 sources from the torus template list and 2 from the ISM list. This guarantees that all the templates can be used to decompose our local sample with spectra in the range 5.5–29 μm. We have chosen this range in order to maximize the number of templates that we are able to use, also maximizing the range covered. Finally, we have also removed from the AGN template library all the sources included in our current analysis (12 AGNs, mainly type-1 and type-2 Seyferts). The final template list contains 101 torus-, 59 ISM-, and 19 stellar-dominated spectra.

We have used the following constraints on the torus component to improve the uncertainty on the torus component as our aim is to study if the torus is present in our sample. First, we have used high spatial resolution data (described in Sections 3.2 and 3.3) to put an upper limit to the torus component flux at the specific observed wavelengths. Note that we are able to include these constraints because we know that the torus component is point like at both Spitzer and high-resolution imaging data. We used the flux measured in the nuclear component plus three times the error. Second, we have used the X-ray luminosity to set an upper limit to the torus component at $12\,\mu {\rm{m}}$, using the well-established relation between the $12\,\mu {\rm{m}}$ and 2–10 keV X-ray luminosities for AGNs (Asmus et al. 2015). We have used the $12\,\mu {\rm{m}}$ to X-ray relation found by Asmus et al. (2015) because it is the most recent relation with the largest number of sources. Since X-ray luminosities can show short- and long-term variations, which can be as high as a factor of 10 (e.g., González-Martín & Vaughan 2012), we have assumed the error on the 2–10 keV X-ray luminosity to be a factor of 10. Then, we have used this value to estimate another upper limit to the $12\,\mu {\rm{m}}$ flux of the torus component. It is worth noticing that the X-ray to mid-infrared relation has been tested down to X-ray luminosities of ${L}_{{\rm{X}}}(2\mbox{--}10\,{\rm{keV}})\sim {10}^{41}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ (Asmus et al. 2015). However, our current sample reaches X-ray luminosities down to ${L}_{{\rm{X}}}(2\mbox{--}10\,{\rm{keV}})\sim {10}^{38}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$. Therefore, this relation might not apply to the very faint end of the luminosity function of AGNs. Indeed, one of the goals of this work is to test the plausible disappearance of the torus, although little is known of the very low-luminosity sources. In this case, we expect the $12\,\mu {\rm{m}}$ flux to be an upper limit to the actual contribution of the torus. Thus, using these limits as upper limits does not bias our results. Furthermore, this upper limit is very useful because it is available for all the objects in our sample.

We have used version 1.2 of the deblendIRS¹¹ code presented by HC15. It improves over previous versions in that it allows the flux priors to constrain the torus component to be included. For each object, we have decomposed its Spitzer/IRS spectrum, using the upper limit on the $12\,\mu {\rm{m}}$ flux (obtained from the 2–10 keV X-ray flux) and the mid-infrared fluxes obtained from high spatial resolution images. Figure 3 shows the best fit for the spectrum of the type-2 Seyfert galaxy NGC 5728 as an example.

Zoom In Zoom Out Reset image size

Tables 1–3 include ${C}_{\mathrm{torus}}$ (Col. 6), ${C}_{\mathrm{stellar}}$ (Col. 7), and ${C}_{\mathrm{ISM}}$ (Col. 8) for the LINER, Seyfert, and Starburst samples, respectively. These values refer to the percentage of each component contributing to the 5–15 $\mu {\rm{m}}$ wavelength range. The Appendix shows the decomposition and posterior distributions for the full sample. The median values and 16th–84th percentiles for each class are listed in Table 6. Starbursts nicely group into ${C}_{\mathrm{torus}}\lt 1.3$%, demonstrating the good performance of this method to decompose Spitzer spectra (see Table 6). Only NGC 3367 shows ${C}_{\mathrm{torus}}\gt 5 \%$ (8.4%). This object might actually host an AGN (see González-Martín et al. 2015). Note that we have also confirmed that the change of the spectral range used for the decomposition to 5–15 μm has no impact in our results.

Table 6.  Median and 16th and 84th Percentiles (in Parentheses) of the Torus, Stellar, and ISM Distributions

	${C}_{\mathrm{torus}}$	${C}_{\mathrm{stellar}}$	${C}_{\mathrm{ISM}}$
	(%)	(%)	(%)
Starbursts	0.1 (0.0–1.3)	4.3 (0.7–19.8)	93.3 (79.2–97.3)
LINER1	6.4 (0.6–36.7)	24.3 (16.2–60.4)	40.1 (22.2–74.3)
LINER2	3.5 (0.1–24.0)	36.0 (1.8–76.7)	51.6 (17.9–85.3)
Seyfert 1	77.2 (70.5–82.2)	14.3 (11.4–16.6)	9.4 (6.1–15.2)
Seyfert 2	58.2 (23.0–80.2)	10.9 (4.1–16.2)	30.7 (12.0–57.9)

Download table as:  ASCIITypeset image

Figure 4 shows a diagram representing the relative contribution of each component for the objects in our sample. Objects with 100% of a single component are located at one of the corners of this diagram (marked with their names in the diagram). Objects with 0% of one of the components would be located in the side of the triangle opposite to its corner. Interestingly, the diagram is not equally populated. There are no objects with ${C}_{\mathrm{ISM}}\sim 0 \%$, indicating that all the objects in our sample show a non-negligible ISM contribution. It is worth noting that each class can be clearly differentiated using their contributors. Starbursts are mainly dominated by ISM components. Type-1 Seyferts are fully dominated by the torus component with less than 30% of the stellar component, and with ${C}_{\mathrm{ISM}}\lt 20 \%$. Type-2 Seyferts have a wide range of torus and ISM components, with less than 30% of the stellar component. Finally, LINERs show less than 50% of the torus component with a wide range of both stellar and ISM components. The small number of type-1 LINERs prevents any firm conclusion on the difference between type-1 and type-2 LINERs to be reached.

Zoom In Zoom Out Reset image size

The DeblendIRS code computes the root-mean-square error (RMSE) of the final fit for each object (included in Tables 1–3, Col. 9). We use these numbers to investigate when the final fit is good enough to represent the data. We have selected as bad fits those with $\mathrm{RMSE}\gt 0.3$ (HC15). Ten objects show RMSE above that limit; two type-2 Seyferts (NGC 3393 and NGC 4945), two Starbursts (NGC 925 and NGC 3184), a type-1 LINER (NGC 4450), and five type-2 LINERs (IRAS 17208-0014, III Zw 035, IRAS 14348-1447, NGC 4125, and IRAS 12112+0305). Note that large RMSE values also imply a larger error on the estimates. We marked these sources in the plots in the following analysis. According to this criterion (RMSE $\lt \,0.3$), DeblendIRS successfully fits the spectra in 90% of the sample. We have not found any trend on the quality of the fit according to the AGN classification. However, most of the objects with large RMSE are actually ULIRGs with deep silicate features (as already discussed by HC15).

Although there is a trend showing that less luminous objects have smaller ${C}_{\mathrm{torus}}$, the relation between the bolometric luminosity and ${C}_{\mathrm{torus}}$ is not linear (coefficient of correlation $r=0.6$). Figure 5 (left panel) shows the expected luminosity below which the torus would disappear (Elitzur & Shlosman 2006) as a red dashed vertical line. Indeed, objects below that limit tend to show small ${C}_{\mathrm{torus}}$. Above that limit, there is a wide range of percentages of torus contributions ${C}_{\mathrm{torus}}$. In order to quantify this (and also to define groups within the plot for subsequent analysis), we have used the clustering Affinity Propagation (AP) method (Frey & Dueck 2007) to look for groups in this diagram. Clustering analysis is aimed at discovering the underlying clusters in the data points according to their similarities. The AP method, in particular, is based on the concept of message passing between data points. The advantage of this method compared to other clustering methods (e.g., k-means) is that AP does not require the number of clusters as an input to the algorithm. The AP method divides the sample into groups and chooses one object as the representative of its group.

We have applied the AP method to the pair of data $[{L}_{\mathrm{bol}}$, ${C}_{\mathrm{torus}}]$ for the full sample using the Python routine "AffinityPropagation" within the package scikit-learn.¹² We used "euclidean" as the affinity method, which uses the negative squared Euclidean distance between points. We set the maximum number of iterations ${\max }\_{iter}$ to 1000 although the actual number of iterations needed to converge was only 43. The routine iterates until the number of estimated clusters does not change for a selected number of iterations ${convergence}\_{iter}$. We set this number of iterations to ${convergence}\_{iter}=15$. Note that the ${convergence}\_{iter}$ parameter is a different restriction than ${\max }\_{iter}$ parameter. However, we have noticed that changes in this ${convergence}\_{iter}$ parameter do not affect the results as long as ${convergence}\_{iter}\gt 1$. This means that the groups in the sample were found after the second iteration.

Using the AP method, the points in Figure 5 (left panel) can be classified into five groups. Col. 10 in Tables 1–3 show which group each object belongs to. Table 7 gives the pair of positions for the representative member of each group and the median (and the 25th–75th percentiles as the width of the distribution) of the objects belonging to each group. The locus of the representative members and their median values are shown in the small inset within Figure 5 (left panel) as squares and black crosses, respectively. The main difference between groups is the percentage of the torus contribution, ${C}_{\mathrm{torus}}$. Among them, only one group is below the line of ${L}_{\mathrm{bol}}={10}^{42}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ and the median ${C}_{\mathrm{torus}}$ for them is consistent with zero (${C}_{\mathrm{torus}}\lt 1.4 \%$; see Table 7). Furthermore, only three objects are consistent with this group but showing ${L}_{\mathrm{bol}}\gt {10}^{42}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ (namely NGC 4945, UGC 05101, and IRAS 14348-1447, shown in the small panel in Figure 5, left panel). Among them, all but UGC 05101 have RMSE $\gt \,0.3$, indicating poor decompositions. Thus, the behavior seen in Figure 5 is consistent with a wide range of ${C}_{\mathrm{torus}}$ above ${L}_{\mathrm{bol}}={10}^{42}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ and negligible ${C}_{\mathrm{torus}}$ below that limit.

Table 7.  Overall Results of the Affinity Propagation (AP) Clustering Method

Group	Representative		Median		Morph.	Slit Width
	Member		of the Group			log(pc)
	$\mathrm{log}({L}_{\mathrm{bol}})$	${C}_{\mathrm{torus}}$	$\mathrm{log}({L}_{\mathrm{bol}})$	${C}_{\mathrm{torus}}$
1	40.8	1.3	40.5 (40.1, 41.3)	0.6 (0, 1.4)	2 (0, 5)	2.6 (2.4, 2.8)
2	42.1	18.5	42.5 (41.9, 43.0)	18.8 (14.1, 23.0)	1 (−1, 3)	2.8 (2.6, 3.1)
3	42.6	39.8	43.5 (43.0, 43.7)	41.2 (37.4, 43.9)	2 (−3, 2)	3.1 (2.9, 3.4)
4	43.7	58.8	43.5 (43.2, 43.7)	60.6 (57.4, 64.8)	1 (0,1)	2.8 (2.7, 3.0)
5	44.1	78.0	44.0 (43.3, 45.0)	78.4 (74.9, 82.3)	1(−2, 4)	3.0 (2.7, 3.2)

Note. The confidence levels computed for the median values are the 25th–75th percentiles; they are not an error on the estimates.

Download table as:  ASCIITypeset image

It is also worth mentioning the results of this classification method compared to the optical classes included in this analysis (see Col. 8 in Tables 1–3). All the Starbursts are classified as Group 1 (i.e., ${C}_{\mathrm{torus}}\lt 1.4 \%$; see Table 7). This is fully expected under the assumption that all these sources are non-AGNs. However, Group 1 also contains four type-2 Seyferts, four type-1 LINERs, and 25 type-2 LINERs. Eight out of the ten type-1 Seyferts are classified within Group 5 (i.e., ${C}_{\mathrm{torus}}\simeq 80 \%$; see Table 7). Only one LINER (NGC 1052) and 11 type-2 Seyferts are associated with Group 5. Except NGC 1052, all LINERs belong to groups below Group 3. In the case of Seyferts, there is a complete mix of groups.

The measurement of ${C}_{\mathrm{torus}}$ could depend on the host galaxy properties. We have collected the morphological types and B magnitudes¹³ for the galaxies in our sample and studied whether the distribution changes with the AP groups. We have computed the 25th–75th percentiles of each distribution for morphological types and B magnitudes (these numbers are recorded in Table 7). There is an increase on the B total magnitude from Group 1 to Groups 4–5, although all of them agree with each other within the percentiles. Furthermore, more distant objects could also include more dust within the slit of Spitzer/IRS. We also computed the 25th–75th percentiles of each distribution of "Slit widths" (see Table 7). There is no correlation between the AP groups and the physical portion of the galaxy included in the IRS slit width.

We have investigated if the percentage of ISM contribution (${C}_{\mathrm{ISM}}$) might affect our results. ${C}_{\mathrm{ISM}}\gt 50 \%$ for Starbursts. This is consistent with the idea that the PAH features dominating the ISM are produced in the photodissociating region associated with a star-forming region. Starbursts, essentially H II galaxies, are dominated by star-forming regions all over the SED, and therefore, dominated by the PAH features at the mid-infrared. Among the AGN classes, both LINERs and Seyferts are spread in a wide range of ${C}_{\mathrm{ISM}}$. This is consistent with the idea that ${C}_{\mathrm{ISM}}$ depends on the inner star formation (∼1 kpc). Indeed, we have studied whether this contribution depends on AGN luminosity, finding no correlation at all (r = 0.18).

We have examined how Figure 5 (left panel) changes if we consider only objects with ${C}_{\mathrm{ISM}}\lt 50 \%$ (see Figure 5, right panel). Now the tendency toward a large ${C}_{\mathrm{torus}}$ for larger AGN bolometric luminosities is much clearer ($r=0.8$). This correlation is still better ($r=0.9$) if we select objects with ${C}_{\mathrm{torus}}\,\lt$ 50% (see the solid orange line in Figure 5, right panel). Furthermore, it is worth noting that objects with ${C}_{\mathrm{torus}}\,\gt$ 50% do not follow the same relation (only two objects are consistent within the error bars). In fact, they tend to have a larger ${C}_{\mathrm{torus}}$ than that expected from the extrapolation of the linear fit of objects with ${C}_{\mathrm{torus}}\lt 50 \%$ (dotted orange line in Figure 5, right panel).

We constructed an average spectrum for each of the five groups found using the AP method (see Figure 6, top panel). The dispersion shown as a shadow area is the 1σ uncertainty of the mean. All the spectra are scaled to their flux at $12\,\mu {\rm{m}}$ before averaging. As in Figure 5 (right panel), we have excluded objects with ${C}_{\mathrm{ISM}}\gt 50 \%$ to avoid those spectra where the ISM is dominating the observed emission.¹⁴ The slope of the spectra gradually becomes bluer as we move to lower ${C}_{\mathrm{torus}}$ (i.e., as we move from Group 4 and 5 to Group 1). We interpret this as stellar dominance for the observed first AP classes (${C}_{\mathrm{torus}}\lt 40 \%$). Our decomposition method allows us to study the average spectrum by isolating the torus component for each AP group. For that purpose, we have subtracted the stellar and ISM components from the Spitzer/IRS spectra. The average spectra are shown in Figure 6 (bottom panel). We do not show in this plot the average spectrum for Group 1 (${C}_{\mathrm{torus}}\lt 1.4 \%$) because, after removing stellar and ISM spectra, the residuals are negligible in this group (see Table 7). Now the average spectra are much more similar among the four groups. The main difference is the presence of strong emission lines in Group 2. Indeed, while collisional lines such as [S iv], [Ne ii], [Ne iii], and [O iv] are present in all the categories, it seems that there is an enhancement of the ${{\rm{H}}}_{2}$ molecular lines in Groups 2 and 3 (${C}_{\mathrm{torus}}\sim 20 \%$ and ${C}_{\mathrm{torus}}\sim 40 \%$, respectively).

Zoom In Zoom Out Reset image size

It is worth mentioning that Group 2 (2, 3, and 4) show the ${{\rm{H}}}_{2}$ molecular line at $28\,\mu {\rm{m}}$ (PAH feature at 11.3 μm) in absorption. This indicates a slight oversubtraction of those emission lines. However, note that all the other ${{\rm{H}}}_{2}$ molecular lines are shown in emission and that the other PAH features (e.g., $8.7\,\mu {\rm{m}}$) are not seen in either emission or absorption. Thus, we believe this oversubtraction depends on the particular conditions of the lines. Perhaps this oversubtraction of the ${{\rm{H}}}_{2}$ molecular lines could indicate a different composition of the clouds. In the case of PAH features, it is even harder to evaluate since changes in the composition of PAHs could yield to a different profile of this line.

We have fitted the subtracted average spectra (Figure 6, bottom panel) with the Clumpy torus models (Nenkova et al. 2008a, 2008b) using Bayesian inference to estimate probability distributions for the parameters of the torus model. Note that emission and absorption features have been removed from the spectral coverage since these models only account for the shape of the continuum. The BayesClumpy code uses the Metropolis–Hastings Markov Chain Monte Carlo algorithm for sampling the posterior distribution function (see Asensio Ramos & Ramos Almeida 2009 for more details on the code). The parameters of the torus involved in the fitting are the half opening angle of the torus $\sigma$, the outer radius of the torus as a function of the inner radius of the torus $Y={R}_{\mathrm{out}}/{R}_{\mathrm{inner}}$, the number of clouds in the equatorial plane ${N}_{{\rm{o}}}$, the slope of the power-law distribution of clouds with respect to the angle from the equatorial plane $q$, the optical depth of the clouds ${\tau }_{\nu }$, and the inclination between the line of sight and the equatorial plane of the torus $i$. Using the best-fit parameters, BayesClumpy marginalizes over the AGN geometrical covering factor, ${f}_{C}$:

$$\begin{eqnarray}&&{f}_{c}=1-{\int }_{0}^{{90}^{^\circ }}{e}^{-{N}_{o}{e}^{-{({90}^{^\circ }-i)}^{2}/{\sigma }^{2}}}\cos (i){di}.\end{eqnarray} \tag{ 2 }$$

The inner radius of the torus is linked to the AGN bolometric luminosity because it depends on the radius at which the dust sublimates (see Barvainis 1987; Nenkova et al. 2008a, 2008b). Thus, the outer radius of the torus can be computed using the bolometric luminosity of the AGN and the parameter ${\rm{Y}}$. The height of the torus H can also be computed as $H=({R}_{\mathrm{out}}\,+{R}_{\mathrm{in}})\sin (\sigma )/2$.

We did not fit Group 1 because after subtracting the stellar and ISM components, only residuals are left. Figure 7 shows the best fit for Groups 2, 3, 4, and 5. We failed to find a good fit for Group 2. Although the average spectrum of Group 2 resembles Groups 3, 4, and 5, below $10\,\mu {\rm{m}}$ it decreases with a slope too steep and it increases too quickly above $25\,\mu {\rm{m}}$ for any of the Clumpy models. We believe that the extra contribution above $25\,\mu {\rm{m}}$ could be a residual of the ISM component. The deficit of emission compared to the model below $10\,\mu {\rm{m}}$ could be due to extra extinction or due to a distribution and/or composition of the dusty region not predicted by the Clumpy models described by Nenkova et al. (2008a, 2008b).

Zoom In Zoom Out Reset image size

We were able to successfully fit Groups 3, 4, and 5, although the spectral fitting for Group 3 is visually worse than that for Groups 4 and 5. The resulting parameters are listed in Table 8. All the parameters are in the range found for AGNs fitted with Clumpy models (Nenkova et al. 2008b; Ramos Almeida et al. 2009, and references therein). Torus sizes recovered by the models are fully in agreement with mid-infrared interferometric observations (Burtscher et al. 2013).

Table 8.  Parameters of the Clumpy Models Obtained using the Code BayesClumpy

Param.	Group 3	Group 4	Group 5
(${C}_{\mathrm{torus}}$)	(40%)	(60%)	(80%)
			All	Type 1	Type 2
σ	${52}_{-11}^{+9}$	${48}_{-13}^{+12}$	${38}_{-14}^{+19}$	${36}_{-11}^{+15}$	${35}_{-7}^{+17}$
Y	${13.2}_{-2.0}^{+2.8}$	${14.6}_{-2.0}^{+4.2}$	${21.0}_{-3.1}^{+7.7}$	${27}_{-8}^{+26}$	${22.0}_{-2.8}^{+5.3}$
Rin (pc)	0.025	0.089	0.14	0.14	0.14
Rout (pc)	${0.33}_{-0.05}^{+0.07}$	${0.82}_{-0.11}^{+0.23}$	${2.36}_{-0.34}^{+0.84}$	${3.0}_{-0.9}^{+2.8}$	${2.42}_{-0.58}^{+0.31}$
H (pc)	${0.14}_{-0.05}^{+0.05}$	${0.34}_{-0.14}^{+0.14}$	${0.76}_{-0.45}^{+0.45}$	$\sim 0.9$	$\sim 0.7$
${N}_{{\rm{o}}}$	${6.3}_{-1.5}^{+2.9}$	${6.5}_{-1.7}^{+3.3}$	${7.0}_{-2.2}^{+3.6}$	${6.7}_{-2.8}^{+3.8}$	${8.5}_{-2.5}^{+3.5}$
${\rm{q}}$	${0.39}_{-0.25}^{+0.45}$	${0.61}_{-0.40}^{+0.72}$	${0.72}_{-0.46}^{+0.68}$	${1.26}_{-0.56}^{+0.46}$	${0.56}_{-0.36}^{+0.58}$
${\tau }_{\nu }$	${71}_{-18}^{+25}$	${110}_{-25}^{+21}$	${64}_{-24}^{+27}$	${42}_{-16}^{+32}$	${43}_{-17}^{+21}$
i	${58}_{-32}^{+18}$	${60}_{-24}^{+16}$	${65}_{-26}^{+13}$	${35}_{-21}^{+22}$	${73}_{-14}^{+9}$
${{\rm{f}}}_{C}$	${0.74}_{-0.24}^{+0.13}$	${0.67}_{-0.24}^{+0.16}$	${0.48}_{-0.29}^{+0.26}$	${0.42}_{-0.18}^{+0.22}$	${0.45}_{-0.14}^{+0.27}$

Note. R_in is estimated from the bolometric luminosity, R_out is estimated from Y R_in, ${\rm{H}}(\mathrm{pc})$ is estimated using R_out, R_in, and $\sigma$, and ${{\rm{f}}}_{C}$ is estimated using the other parameters with Equation (2). Group 1 was not fitted. We did not find a good fit for Group 2 (see text). We have also fitted the average spectrum of type-1 and type-2 Seyferts included in Group 5 (see the text).

Download table as:  ASCIITypeset image

Rather than focus on the actual numbers obtained for each group, it is more relevant to look for differences in the measured parameters among the three groups. The number of clouds in the equatorial plane (${N}_{{\rm{o}}}$) and the viewing angle (i) is similar for the three groups. However, when we split Group 5 into type-1 and type-2 AGNs¹⁵ (see the last two columns in Table 8), we naturally recover the dependence on the viewing angle (as expected under the unified model). There is a marginal trend for an increase of the slope of the cloud distribution q with the distance from the equatorial plane, the height of the torus, and a decrease on the half opening angle $\sigma$ from Group 3 to Group 5. The latter results in a marginal detection on higher geometrical covering factors for lower luminosity AGNs.

The parameter that is changing the most among the three groups is the outer size of the torus (R_out), increasing from Group 3 to Group 5 (i.e., for ${C}_{\mathrm{torus}}$ from 40% to 80%). Alonso-Herrero et al. (2011) did not find a different size of the torus, fitting 13 AGNs to the very same Clumpy models used here. However, there are two main reasons why these results are not in contradiction. First, the range of luminosities covered here is much larger than those reported in their sample. Indeed, the sizes of the torus and bolometric luminosities reported in their publications are consistent with all of them being mainly in Group 5. Second, we have found this tendency after averaging objects classified into the same group using the AP method while they used individual Clumpy fits to derive their conclusions. Since the torus half opening angle ($\sigma$) is similar and the size of the torus (R_out) is smaller when moving from Group 5 to Group 3, the covering factor (the half opening angle of the torus) of the AGN is larger (smaller) for Group 3 compared to Group 5. Note that far-infrared flux observations would be more sensitive to the torus extent independently (Ramos Almeida et al. 2011). Thus, the addition of these measurements could confirm the trend found here for the AP groups (i.e., with ${C}_{\mathrm{torus}}$).

Motivated by the ${{\rm{H}}}_{2}$ molecular emission lines detected in Groups 2 and 3 found using the AP method, we have done a search for all the ${{\rm{H}}}_{2}$ molecular lines found in the Spitzer/IRS spectra of the sample. Due to the coverage of our spectra, we have been able to look for the transitions S(0) $28.22\,\mu {\rm{m}}$, S(1) $17.03\,\mu {\rm{m}}$, S(2) $12.27\,\mu {\rm{m}}$, S(3) $9.67\,\mu {\rm{m}}$, S(4) $8.02\,\mu {\rm{m}}$, S(5) $6.91\,\mu {\rm{m}}$, and S(6) $6.11\,\mu {\rm{m}}$. We have built a code using Python routines (within scipy) to automatically detect these lines. We have forced a fit to a Lorentzian profile to each of the proposed lines. We allowed the amplitude, center, and width of the Lorentzian profile to vary. We have run this fitting 200 times using Monte Carlo simulations to estimate the error of these fittings from the error of each spectrum. We consider a line to be detected if the amplitude of the line is above three times its error and the center of the line is consistent with the expected wavelength of the line within the width of the Lorentzian. Note that we have performed such a detection in the Spitzer/IRS spectra and not after the stellar and ISM components were subtracted.

We have not detected S(0) in any of the objects of our sample. Furthermore, none of the ${{\rm{H}}}_{2}$ molecular lines are detected in any of the Starbursts and type-1 Seyferts included in our comparison sample. S(5) and S(6) lines are detected among both LINERs and type-2 Seyferts. S(4) is detected only in five objects, three in Group 1 and two in Group 5. S(1), S(2), and S(3) are only detected in Groups 1, 2, and 3 (i.e., ${C}_{\mathrm{torus}}\lt 40 \%$).