Exploring the Mid-infrared SEDs of Six AGN Dusty Torus Models. I. Synthetic Spectra

Barvainis (1987) was the first to propose that the excess of emission observed at ∼3 μm in active galactic nuclei (AGN), compared to the disk emission at optical/UV frequencies, was probably due to multiple-temperature dust components. Since then, many attempts have been made to solve the radiative transfer equations to produce the spectral energy distribution (SED) due to dust.⁸ Pioneering works on the subject were Krolik & Begelman (1988), Pier & Krolik (1992), Granato & Danese (1994), Stenholm (1994), Efstathiou & Rowan-Robinson (1995), and Nenkova et al. (2002).

Initially, for the sake of simplicity, most authors used smooth dust distributions with different radial and vertical density profiles (e.g., Pier & Krolik 1992; Granato & Danese 1994; Efstathiou & Rowan-Robinson 1995; van Bemmel & Dullemond 2003; Schartmann et al. 2005; Fritz et al. 2006). However, before the development of the first geometrical torus models, it was known that dust was probably organized into clouds because it would be very difficult for the dust to survive in the region otherwise. The dust should be moving with random velocities of hundreds of kilometers per second to reproduce the geometrical thickness required to obscure a substantial solid angle around the central source. If the dust were distributed homogeneously, the temperature would be of the order of ∼10⁶ K, preventing the dust from surviving (e.g., Krolik & Begelman 1988; Tacconi et al. 1994). Several radiative transfer models were developed to account for 2D or 3D clumpy dust distributions (Nenkova et al. 2008a, 2008b; Hönig & Kishimoto 2010; Hönig et al. 2010). A mix of smooth plus clumpy distributions have also been proposed (Stalevski et al. 2012; Siebenmorgen et al. 2015). More recently, motivated by the discovery of a polar dust contribution to the mid-infrared SEDs, a more complex scenario has been proposed to explain the infrared nuclear emission of Seyfert galaxies (e.g., Hönig et al. 2013; Asmus et al. 2016). Hönig & Kishimoto (2017) produced a model that includes a compact, geometrically thin disk in the equatorial region of the AGN and an extended, elongated polar structure, which is cospatial with the outflow region of the AGN on larger scales (see also the recently published semi-empirical SED libraries for quasars by Lyu & Rieke 2018). They claim that this model reproduces the 3–5 μm bump observed in many type 1 AGN that other models could not reproduce. Using infrared interferometry of 23 Seyfert galaxies (for which 7 sources have three or more well-spaced directions in the (u, v) plane, the point-like source contributes less than 70%, and the uncertainties are below 10%), this mid-infrared polar emission has been detected so far in six sources: five by López-Gonzaga et al. (2016) and an additional one reported by Leftley et al. (2018). We give a short overview of the six models used in this paper in Section 2 (see also Ramos Almeida & Ricci 2017, for further details of different AGN dust models).

Apart from the obvious differences in the torus geometry and smoothness/clumpiness, these models also invoke different grain sizes and chemical compositions. Fritz et al. (2006) considered typical silicate and graphite grain sizes of a_G = 0.005–0.25 μm and a_G = 0.025–0.25 μm, respectively. Nenkova et al. (2008a) assumed that dust grains are spherical, with size distribution from Mathis et al. (1977) and a standard interstellar medium (ISM) composition of 53% silicates and 47% graphites. Hönig & Kishimoto (2010) included the standard ISM composition, although they also used ISM large grains (sizes between 0.1 and 1 μm), dominated by intermediate to larger graphite grains.

Minor effort has been put into comparing the various models or disentangling which of them better reproduces the data. Most comparisons are confronting parent models developed by the same group (Schartmann et al. 2008; García-González et al. 2017). Although some papers have tried to qualitatively compare different models (e.g., van Bemmel & Dullemond 2003), Feltre et al. (2012) is the only work comparing the smooth model by Fritz et al. (2006) and the clumpy model by Nenkova et al. (2008b) using matched parameters. They found that the behavior of the silicate features and spectral slopes below ∼7 μm are different.

The question we would like to answer in this paper is whether the data allow us to distinguish which of the proposed models better describes AGN infrared SEDs. For this purpose we created a set of synthetic spectra using six of these well-known models (see Section 2) using different instruments. A comparison between models and real AGN mid-infrared spectroscopic data is performed by González-Martín et al. (2019, hereafter Paper II). The paper is organized as follows. Section 2 gives a brief summary of the dusty models used throughout this paper. Section 3 describes how to convert multiparameter models to the spectral fitting tool environment XSPEC (Arnaud 1996), how to create synthetic spectra using the most relevant current and future mid-infrared instrumentation, and how to fit these synthetic spectra. The main results are included in Section 4. These results are discussed in Section 5, and the paper ends with a report of the main findings in Section 6.

This paper is devoted to comparing as many dust SED libraries as possible. Table 1 shows the compilation of radiative codes with publicly available SEDs. All of them show slightly different dust distributions (Col. (2)) and compositions (Col. (3)). Furthermore, they map a different range and/or set of parameters (see Col. (6)). There are other available radiative calculations to model the dusty structure of AGN, but there is no available SED library to test them. Furthermore, note that there are semi-empirical SED libraries that are also beyond the scope of this paper (e.g., Lyu & Rieke 2018). Below we describe the SED libraries used. We quote the abbreviation used for each model throughout this paper within brackets.

• 1.  
  Smooth torus model by Fritz et al. (2006) [Fritz06] (see also Feltre et al. 2012). Radiative transfer code used to produce the SED of a simple toroidal geometry consisting of a flared disk that can be represented as two concentric spheres, delimiting respectively the inner and the outer torus radius, having the polar cones removed. They considered typical silicate and graphite grain sizes of a_G = 0.005–0.25 μm and a_G = 0.025–0.25 μm, respectively. This model uses different sublimation radii for graphite and silicate grains. The free parameters of this model are the viewing angle toward the toroidal structure, i, the opening angle, σ, the exponent of the logarithmic azimuthal density distribution, Γ, the exponent of the logarithmic radial profile of the density distribution, β, the outer radius of the torus, R_max, compared to the inner one (expressed as Y = R_max/R_min), and the edge-on optical depth at 9.7 μm, τ_{9.7 μm}. The size of the torus is defined by the outer radius, R_max, and the opening angle, σ. The inner radius is defined by the sublimation temperature of dust grains under the influence of the strong nuclear radiation field (T_sub = 1500 K). Its incident radiation spectrum is defined by a combination of power laws such as those described by Granato & Danese (1994). The SED library contains 24,000 SEDs in the 0.001–1000 μm wavelength range.
• 2.  
  Clumpy torus model by Nenkova et al. (2008b) [Nenkova08]⁹ (see also Nenkova et al. 2008a). They developed a formalism that accounts for the concentration of dust in clumps or clouds, referred to as clumpy, to describe the nature of the AGN torus. Nenkova et al. (2008a) assumed that dust grains are spherical, with size distribution from Mathis et al. (1977) and a standard Galactic mix of 53% silicates and 47% graphites (i.e., standard ISM). The model assumes a toroidal distribution that depends on the viewing angle toward the torus, i, the inner number of clouds, N₀, the half-opening angle of the torus, σ, the outer radius of the torus, R_out (scaled to the inner radius, i.e., Y = R_out/R_in), the slope of the radial density distribution, q, and the optical depth (i.e., at 0.55 μm) of the individual clouds, τ_v. The inner radius is fixed to the dust sublimation radius (T_dust = 1500 K). They produced a library including 1,247,400 SEDs in the 0.001–1000 μm wavelength range.
• 3.  
  Clumpy torus model by Hönig & Kishimoto (2010) [Hoenig10] (see also Hönig et al. 2006, 2010). Radiative transfer model of three-dimensional clumpy dust tori using optically thick dust clouds and a low torus volume filling factor. They included an improved handling of the diffuse radiation field in the torus, which is approximated by a statistical approach and the possibility of different dust compositions and grain sizes. They synthesized three compositions: the standard ISM (47% graphites and 53% silicates), ISM large grains (grains between 0.1 and 1 μm in size), and Gr dominated (dominated by intermediate to larger graphite grains; 70% graphites and 30% silicates). The parameters of this library of SEDs are the viewing angle, i, the number of clouds along an equatorial line of sight, N₀, the half-opening angle of the distribution of clouds, θ, the radial dust cloud distribution power-law index, a, and the opacity of the clouds, τ_cl. The outer torus radius, R_max, is fixed to the inner radius as R_max = 150 R_min, where R_min is set to the dust sublimation radius. This library includes 1680 SEDs in the 0.01–36,000 μm wavelength range.
• 4.  
  Two-phase (clumpy + smooth) torus model by Siebenmorgen et al. (2015) [Sieben15]. They assumed that the dust near the AGN is distributed in a torus-like geometry, which can be described as a clumpy medium or a homogeneous disk, or a combination of the two. They include an isothermal disk that is embedded in a clumpy medium, including parsec-sized dust clouds passively heated. Moreover, this model is based on the idea that a significant part of the mid-infrared emission appears to come from the ionization cones (Braatz et al. 1993; Hönig et al. 2013). Therefore, in this model dust exists also in the polar region. The dust particles considered are fluffy and have higher submillimeter emissivities than grains in the diffuse ISM. The SED library depends on the viewing angle, i, the inner radius, R_in, the volume filling factor of the clouds, η, the optical depth (i.e., at 0.55 μm) of the individual clouds, τ_cl, and the optical depth of the disk midplane, τ_disk. The outer radius up to where dust exists and their distribution function are kept constant as R_out = 170R_in. They produced a library including 3600 SEDs in the 0.0005–500 μm wavelength range.
• 5.  
  Two-phase (clumpy + smooth) torus model by Stalevski et al. (2016) [Stalev16] (see also Stalevski et al. 2012). They model the dust in a torus geometry with a two-phase medium, consisting of a large number of high-density clumps embedded in a smooth dusty component of low density. This model is claimed to reproduce the near-infrared excess around 5 μm, producing at the same time attenuated silicate features (Stalevski et al. 2012). One of the improvements of the current set of models described by Stalevski et al. (2016) compared to the previous modeling (i.e., Stalevski et al. 2012) is the inclusion of anisotropic emission. This implies that the inner radius is not constant but shows the same dependency with the polar angle shown by the anisotropy. Dust is also distributed along the radial and polar directions according to a power-law distribution (based on the same prescription given by Fritz et al. 2006). The dust chemical composition is set to a mixture of silicate and graphite grains. The parameters of the model are the viewing angle toward the observer, i, the ratio between the outer and the inner radius of the torus, Y = R_out/R_in, the half-opening angle of the torus, σ, the indices that set the dust density gradient with the radial p and polar q distribution of dust, and the 9.7 μm average edge-on optical depth, τ_{9.7 μm}. The fraction of total dust mass in clumps compared to the total dust mass is set to M_cl = 0.97, and the inner radius of the torus is set to R_in = 0.5 pc. This SED library contains 19,200 SEDs in the 0.001–1000 μm wavelength range.
• 6.  
  Clumpy disk and outflow model by Hönig & Kishimoto (2017) [Hoenig17]. This model is built on the suggestion that the dusty gas around the AGN consists of an inflowing disk and an outflowing wind. In practice, this model consists of clumpy disk-like models (following that described by Hönig & Kishimoto 2010) plus a polar outflow. Common parameters for disk and wind are the viewing angle, i, and the number of clouds in the equatorial plane, N₀. The disk is also governed by the exponent of the radial distribution of clouds, a, and the optical depth of individual clouds, τ_cl, which is fixed to τ_cl = 50. The outflow is modeled as a hollow cone and characterized by three parameters: the index of the dust cloud distribution power law along the wind, a_w, the half-opening angle of the wind, θ, and the angular width of the hollow wind cone, σ. Finally, a wind-to-disk ratio, f_wd, defines the ratio between the number of clouds along the cone and N₀. The library includes 132,300 SEDs in the 0.01–36,000 μm wavelength range.

A caveat to take into account when comparing all these models is that they are not produced to fully cover similar dusty structures. Indeed, this is impossible when producing SEDs with completely different structural components. The exception is the parameter space from [Fritz06] and [Nenkova08] because the former libraries were produced to match similar physical parameters of the torus (see Feltre et al. 2012). However, we rely on the assumption that the parameter space for all the models is set by the authors to better match observational evidence of the AGN dusty structure. In that sense, the comparison shown in this paper is still valid. We defer the study of the adequacy of the parameter space by comparison to an AGN sample with Spitzer/IRS spectroscopic data to Paper II.

We devoted this paper to testing whether the models can be self-constrained and whether they can be differentiated in their spectral shape among each other. For that purpose we converted the SED libraries to multiparametric models within the spectral fitting tool XSPEC (Section 3.1) and produced a set of synthetic spectra (Section 3.2) to test these models.

3.1. Model Conversion into XSPEC Format

3.2. Synthetic Spectra in XSPEC Format

We create synthetic spectra to test the models. We simulated mid-infrared GTC/CanariCam, mid-infrared low-resolution Spitzer/IRS, near-infrared James Webb Space Telescope (JWST)/NIRSpec, and mid-infrared JWST/MIRI spectra to compare them against models. We do not attempt to test models against photometric data because the large variety of filters and their combinations add complexity to the analysis (see, e.g., Ramos Almeida et al. 2014). XSPEC is able to simulate spectra with the fakeit task using the instrument response¹² and a model. The flx2xsp task, within HEASOFT, reads a text file containing one or more spectra and their errors and writes out a standard XSPEC pulse height amplitude (PHA¹³ ) and response files.

We convert into XSPEC format a real spectrum in order to get the response file to produce the synthetic spectra. Note that this response file does not depend on the actual spectrum of the object used but on the instrumental setup. The information on the spectrum is recorded in a different file, while the response file contains only the sensitivity/performance of the instrument used. We used the N-band (∼7–13 μm) and Q-band (∼17–23 μm) T-ReCS/Gemini spectra for NGC 1068 downloaded from the Gemini archive¹⁴ and reduced as a point-like source with the RedCan pipeline (González-Martín et al. 2013). Note that this is representative of the expected wavelength coverage and sensitivity of GTC/CanariCam because CanariCam and T-ReCS are twin instruments. T-ReCS is no longer available for the community, whereas CanariCam is currently available. GTC/CanariCam mid-infrared simulated spectra are also equivalent to other ground-based facilities (e.g., VLT/VISIR or Michelle/Gemini) owing to their similarities in the wavelength range and sensitivities.

As for the Spitzer/IRS response, we obtained the Spitzer/IRS low-resolution spectra for NGC 1052 using the Combined Atlas of Sources with Spitzer IRS Spectra (CASSIS¹⁵ ), which provides fully calibrated 1D spectra for both low-resolution and high-resolution Spitzer instruments.

The JWST/MIRI and JWST/NIRSpec spectra were simulated using the JWST exposure time calculator (ETC¹⁶ ) tool. This tool provides simulated JWST data for any uploaded spectrum. We created MIRI integral field unit (IFU) data cubes for the 12 wavelength ranges (three channels and four wavelength ranges per channel) and NIRSpec IFU data cubes using G140H, G235H, and G395H filters. We uploaded one of the SEDs produced by [Nenkova08] to produce the simulated data cube for NGC 1068 using a point-like source as the proposed model distribution. In order to get flux-calibrated spectra, we also produced the same data cubes for a flat spectrum (i.e., with a known constant flux across the wavelength range) to produce a wavelength-dependent count-to-flux conversion. We created the simulated spectra by multiplying the simulated [Nenkova08] SED by this count-to-flux conversion.

A schematic view of the procedure to produce and fit synthetic spectra can be seen in Figure 1. Our synthetic spectra include several instrumental setups as a combination of four instruments (GTC/CanariCam, Spitzer/IRS, JWST/MIRI, and JWST/(MIRI + NIRSpec) and four sensitivities (F(12 μm) = 30 mJy, 100 mJy, 300 mJy, and 10 Jy). We use the six AGN dust models alone and a combination of AGN dust models and circumnuclear contributors (stellar and/or ISM contributors). The stellar contribution is set to 10%, 50%, and 100% of the AGN at 5 μm and 100% of the AGN at 10 μm. The ISM contribution is set to 10%, 50%, and 100% of the AGN at 30 μm. We also explore a combination of 10%, 50%, and 100% of stellar and ISM contribution at 5 and 30 μm, respectively. Note that only the F(12 μm) = 300 mJy was used when the circumnuclear contributors were added to the synthetic spectra. We produce 96,000 SEDs using random realization of the parameters for the synthetic spectra. This includes 1000 SEDs per model and instrumental setup. For the only AGN synthetic spectra we also produce 147,744 SEDs using three fixed values per model, including 243, 729, and 6561 SEDs for the five-, six-, and eight-parameter models. These three values are chosen as 2/5, 3/5, and 5/5 of the parameter space. Moreover, we generate 29,600 SEDs using random realizations of the parameters when AGN dust models are combined with circumnuclear contributors. This includes 100 SEDs per model and 100 additional SEDs for [Fritz06] and [Hoenig17] to explore whether the limited number of SEDs affects our results. Overall we produced 273,344 synthetic SEDs and ∼1.5 million spectral fits. This took over 3 months of computational time in three 64-core, one 32-core, and two 8-core dedicated servers.¹⁷

Zoom In Zoom Out Reset image size

We use the synthetic spectra to study the spectral shapes produced by these models (Section 4.1), to test the accuracy of the parameter determination per model (Section 4.2), whether these models can be distinguished among them (Section 4.3), and how the recovery of the parameters is affected by unresolved circumnuclear contributors to the spectra (Section 4.4).

4.1. Spectral Shapes

We produced a set of 1000 spectra per model using random parameters. The synthetic SEDs with fixed parameters will be used to ensure a good coverage of the parameter space. For each of them, we used the fakeit task to produce synthetic spectra for GTC/CanariCam, Spitzer/IRS, JWST/MIRI, and JWST/(MIRI + NIRSpec) instruments that can be fitted within XSPEC. We scaled the normalization parameter to four 12 μm fluxes (∼30 mJy, 100 mJy, 300 mJy, and 10 Jy) to simulate different signal-to-noise ratio (S/N) (equivalent to S/N ∼ 3–10, 20–30, 40–60, and 100–150, respectively). Therefore, we produced 12,000 SEDs per model using random realizations (1000 SEDs per model, four instruments, and three S/N levels) and 9234 SEDs using fixed steps for the six models studied (times the 12 instrumental setups; see Figure 1). These SEDs are used to study the accuracy of the parameter determination in Section 4.2.

We compute the following spectral parameters to compare the spectral shapes provided by these models, following previous works on the infrared spectra of AGN (e.g., Hernán-Caballero et al. 2015; García-González et al. 2017; Hönig & Kishimoto 2017), although optimized to the spectral range provided by Spitzer/IRS. Note that this analysis has been focused on the synthetic spectra provided by Spitzer/IRS for consistency with the analysis performed in Paper II.

• 1.  
  Spectral slopes: We computed three spectral slopes of the form $\alpha =-\mathrm{log}({F}_{\nu }({\lambda }_{2})/{F}_{\nu }({\lambda }_{1}))$/$\mathrm{log}({\lambda }_{2}/{\lambda }_{1})$ (where λ₂ > λ₁). Note that, under this nomenclature, negative (positive) values mean that the flux increases (decreases) with wavelength. We defined α_NIR, α_MIR, and α_FIR as the slopes evaluated at [λ₁, λ₂] equal to [5.5, 7.5], [7.5, 14], and [25, 30] μm, respectively.
• 2.  
  Silicate features strength: We also computed the silicate feature strength using the formula ${{\rm{Si}}}_{\lambda }=-\mathrm{ln}({F}_{\nu }(\lambda )/{F}_{\nu }(\mathrm{continuum}))$, for the two silicate features located at ∼9.7 and ∼18 μm. Silicate features in emission (absorption) show negative (positive) Si_λ. We anchor the continuum at the ∼9.7 μm (∼18 μm) silicate feature using a linear fit to the continuum at 7–7.5 μm and 14–15 μm (14–15 μm and 25–26 μm).

Figures 2–4 show α_MIR versus α_NIR, α_FIR versus α_MIR, and Si_{18 μm} versus Si_{9.7 μm}, respectively. Table 2 shows the range of values for each parameter and model. Synthetic spectra produced with fixed steps (small circles) always show spectral parameters well within those obtained using the synthetic spectra produced with random parameters (large circles). This ensures that 1000 random realizations are enough to cover the spectral shapes covered by the SEDs.

Zoom In Zoom Out Reset image size

Table 2.  Range of the Spectral Parameters per Model

	αNIR	αMIR	αFIR	Si9.7 μm	Si18 μm
F06	[−3.5, −0.5]	[−3.0, −1.0]	[−4.0, 1.0]	[−1.5, 4.0]	[−0.5, 0.6]
N08	[−7.0, −1.0]	[−7.5, −1.0]	[−5.0, 0.0]	[−1.5, 1.5]	[−0.6, 0.4]
H10	[−7.0, 0.0]	[−4.0, 0.0]	[−1.0, 3.0]	[−1.0, 1.0]	[−0.6, 0.2]
S15	[−8.0, −1.0]	[−8.0, −2.0]	[−6.0, 0.0]	[−1.0, 5.0]	[−0.4, 0.4]
S16	[−3.5, −1.0]	[−3.0, −1.0]	[−2.5, 0.5]	[−1.5, 2.5]	[−0.4, 0.3]
H17	[−8.0, 0.0]	[−4.0, 0.0]	[−1.0, 2.0]	[−0.8, 0.8]	[−0.6, 0.0]

Note. F06: [Fritz06]; N08: [Nenkova08]; H10: [Hoenig10]; S15: [Sieben15]; S16: [Stalev16]; and H17: [Hoenig17].

Download table as:  ASCIITypeset image

Figures 2 and 3 show that all the models overlap in a range of α_NIR and α_MIR with ∼[−3, −1.0] and in a range of α_FIR with ∼[−1, 0], except for [Sieben15], which lacks α_MIR in the range [−2, −1] (see also Table 2). [Nenkova08] ([Sieben15]) shows values of these slopes as low as α_NIR ∼ −7 (α_NIR ∼ −8), α_MIR ∼ −7.5 (α_MIR ∼ −8), and α_FIR ∼ −5 (α_FIR ∼ −6). [Hoenig10] and [Hoenig17] also extend α_NIR and α_MIR to low values with α_NIR ∼ −8 and α_MIR ∼ −4 but also to larger values up to α_NIR ∼ 0 and α_MIR ∼ 0. Some of these characteristics might explain the need for nonrealistic stellar and ISM components by the data when using some of the models (see Paper II).

Zoom In Zoom Out Reset image size

The strength of the silicate features also differs depending on the model used (see Figure 4). All of them overlap in the range of Si_{9.7 μm} = [−1, 1] and Si_{18 μm} = [−0.4, 0], i.e., producing relatively weak emission features. However, [Hoenig10] and [Hoenig17] almost do not include SEDs with 18 μm absorption features, while [Fritz06] contains SEDs with 18 μm silicate features with strengths as high as Si_{18 μm} = 0.6. Similarly, [Fritz06] and [Sieben15] also contain SEDs with absorption 9.7μm silicate features with strengths as high as Si_{18 μm} = 4–5. These plots and Table 2 could be used to study the adequacy of the models to AGN mid-infrared spectra (see also Paper II).

Zoom In Zoom Out Reset image size

4.2. Parameter Determination

We automatically fitted the synthetic SEDs used in the previous section with the same model used to produce them. This allows us to study how accurate is the parameter determination for these instruments. Note that in order to fit [Sieben15] using the JWST/(MIRI + NIRSpec) combination of instruments, we needed to rule out the spectral range below 4 μm because their SEDs show a chaotic behavior in the near-infrared domain, which produces a failure in the chi-squared procedure.

Figure 5 shows the results obtained for [Nenkova08] with f(12 μm) = 300 mJy. Each panel shows the estimated versus simulated values for one of the parameters of the model. The results using random parameters for the four instrumental setups are included in each panel using different colors and symbols in the bottom panels. The top panels show the histogram of the resulting values obtained for the synthetic spectra produced using fixed parameter steps (the simulated parameters are marked with vertical dashed lines). The best estimates of the parameters are obtained for Y, q, and τ_v. This is clearly seen in the small dispersion around the simulated value using either random or fixed parameter values. Ramos Almeida et al. (2014) found that Y is sensitive to wavelengths above ∼20 μm while q and τ_v can be restricted through the silicate features. However, i and σ are more sensitive to the near- to mid-infrared slope. This might explain why Y, q, and τ_v are better recovered than i and σ.

Zoom In Zoom Out Reset image size

The best way to analyze the parameter determination is to produce the posterior distribution function (pdf) per parameter. However, producing pdf's requires ∼5 minutes per parameter and spectral fit. Since 250,000 fits requires longer than 2 yr computer time, we performed this calculation only for the Spitzer/IRS data in Paper II. The results obtained in Paper II are consistent with the analysis performed here. In the synthetic spectra we estimate the accuracy on the parameter determination using as the error the difference between the computed and simulated values compared to the parameter space range. We consider that a parameter is well determined if this error is within 15% of the parameter space on average. Figure 6 shows the average percentage error compared to the parameter space range per model and instrumental setup. The 15% of the parameter space is highlighted as the first dashed circle. The results for different S/N simulated spectra are linked with cyan solid, blue long-dashed, and red short-dashed lines. Except for [Hoenig17], all the parameters are recovered well within 15% of error for a source with f(12 μm) = 100 mJy. However, sources with lower fluxes (green-filled area in Figure 6) show errors larger than 15% of the parameter space. Note that the results are irrespective of the use of synthetic spectra using random initial parameters (solid lines) or fixed parameter values (dashed lines). We can still recover four (i, h, a, and a_w) out of the eight parameters for [Hoenig17] using any of the instrumental setups. The large errors obtained for [Hoenig17] are associated with the large number of parameters involved in this model.

Zoom In Zoom Out Reset image size

It is also worth mentioning that when the angular width of the torus is a free parameter on the models (i.e., [Fritz06], [Nenkova08], [Hoenig10], [Stalev16], and [Hoenig17]), it is the less constrained parameter, except for [Stalev16]. Several authors have pointed out that the viewing angle of the torus is very difficult to constrain without near-infrared data, at least for [Nenkova08] (Ramos Almeida et al. 2011, 2014; Alonso-Herrero et al. 2011). We found that the viewing angle is well constrained for all the models analyzed, even using GTC/CanariCam data, although this can be improved with JWST/(MIRI+NIRSpec) spectra. This could be explained by the fact that we use here full-coverage mid-infrared spectra while they used photometry and N-band spectroscopy. Moreover, the inclusion of host galaxy dilution might be key to understanding why near-infrared data are needed (see Section 4.4). [Sieben15] is particularly well determined with less than 10% error for all the instrumental setups, except for the faintest sources with f(12 μm) ∼ 30 mJy.

In a broad comparison among instruments, an increment on the errors below 5% is obtained for GTC/CanariCam data compared to other instruments. The difficulty in the determination of parameters using GTC/Canaricam data is related to the limited wavelength range. Furthermore, Spitzer/IRS spectra are better suited to getting a good determination of the models than JWST/MIRI, and only a marginal improvement is obtained when using JWST/(MIRI+NIRSpec). This is due to the larger wavelength range of the Spitzer/IRS spectra (5–38 μm) compared to JWST/MIRI (5–30 μm). However, we emphasize here that Spitzer/IRS spectra include a larger portion of the galaxy compared to JWST/MIRI or ground-based GTC/CanariCam, imposing a further limitation on the determination of the torus parameters (see Section 4.4).

4.3. Model Discrimination

We investigate in this section whether the spectrum from a particular model can be distinguished from the others. We used for this purpose the 6000 synthetic Spitzer/IRS spectra (1000 spectra per model) obtained with a 12 μm flux of 300 mJy (same as those used in Section 4.2; see Figure 1).

Each synthetic spectrum is fitted with the six models, including the one used to produce the synthetic spectrum. Figure 7 shows the residuals (left) and distribution of χ²/dof (right) obtained when fitting the 1000 synthetic spectra produced with [Nenkova08]. The distributions of reduced χ² show that the only model able to obtain a good fit for these synthetic spectra is its own model. The main differences are (1) the 9.7 and 18 μm silicate features; (2) deficit (for [Fritz06]), excess (for [Hoenig10] and [Hoenig17]), or both (for [Sieben15]) on the flux below 7 μm; and (3) steep slopes above ∼25 μm for [Hoenig10] and [Hoenig17]. Figures 11 and 12 show the results for all the simulated spectra (one panel per model). In general, all the models can be distinguished based on the χ² statistics and residuals except for the comparisons [Hoenig10] versus [Hoenig17] and [Fritz06] versus [Stalev16]. This is somewhat expected because the CAT3DWIND model reported by Hönig & Kishimoto (2017) is based on the CAT3D model reported by Hönig & Kishimoto (2010) and the two media torus model produced by Stalevski et al. (2016) uses similar geometry and dust composition to that used by Fritz et al. (2006). Nevertheless, it is worth noticing that [Stalev16] synthetic spectra are more easily recovered with the [Fritz06] model than the other way around. This is due to the largest parameter space covered by [Fritz06] compared to [Stalev16]. Note that [Stalev16] has enough spectral coverage to reproduce real data (see Paper II). In general, the main differences are the same, i.e., a different slope below ∼7 μm, a different shape for the 10 and 18 μm silicate features, and a steeper slope above ∼25 μm for [Hoenig10] and [Hoenig17] compared to [Fritz06], [Nenkova08], [Sieben15], and [Stalev16].

Zoom In Zoom Out Reset image size

We also repeat our analysis for a source with f(12 μm) ∼ 10 Jy and 100 mJy. In the former, the models are so different that we could not converge to a good solution. In the case of a source with f(12 μm) ∼ 100 mJy, our results are similar to those obtained for a source with f(12 μm) ∼ 300 mJy. Thus, except for the combinations [Hoenig10] versus [Hoenig17] and [Fritz06] versus [Stalev16], models can be distinguished based solely on Spitzer/IRS spectra. Note, however, that we also repeated the analysis with the f(12 μm) ∼ 30 mJy source. At this low S/N limit models are indistinguishable irrespective of the instrumental setup or model used.

4.4. Host Galaxy Dilution

In previous sections, we showed that mid-infrared spectra are useful to discern which is the best model and to constrain the parameters of the models for relatively low flux AGN (f(12 μm) > 100 mJy). However, a question remains open: the effect of the dilution by external contributors (i.e., host galaxy) on the estimated parameters. In order to analyze that, we have included 100 synthetic spectra with stellar components, ISM components, and a combination of the two per AGN model used. Note that we also produced 100 additional synthetic spectra for [Fritz06] and [Hoenig17] to ensure that the relatively low number of SEDs does not affect our result. The ISM component is taken from Smith et al. (2007), which are averaged Starburst templates in the ∼5–160 μm wavelength range for different 6.2, 7.7, 11.3, and 17 μm PAH feature strengths (see their Figure 13). Note that the results for ISM dilution using JWST/(MIRI+NIRSpec) are not accurate because these templates lack spectral information at near-infrared wavelengths. The stellar component corresponds to a stellar population of 10¹⁰ yr and solar metallicity from the stellar libraries provided by Bruzual & Charlot (2003). We set the parameters of the ISM component to the first ISM template for simplicity. The stellar component contributes mostly to the shortest wavelengths (<10 μm), and the ISM component contributes to the entire mid-infrared wavelength range, with the highest contribution at longer wavelengths (>20 μm).

These libraries have been converted into XSPEC following the same procedure explained in Section 3.1. We produced four combinations for each synthetic spectrum using a stellar contribution scaled to 10%, 50%, and 100% of the torus flux at 5 μm and 100% of the torus flux at 10 μm. Note, however, that we include 100% contribution of the stellar component at mid-infrared wavelength for completeness, although it is well documented in the literature that this percentage is never reached at mid-infrared (stellar contribution always below 50%–60% for Sy2s and almost negligible for Sy1s; see Rodriguez Espinosa et al. 1987; Dultzin-Hacyan et al. 1988; Dultzin-Hacyan & Benitez 1994; Dultzin-Hacyan & Ruano 1996). Similarly, we produced three combinations for each synthetic spectrum using an ISM contribution scaled to 10%, 50%, and 100% of the torus flux at 30 μm (see Figure 1). Finally, we also produced 100 synthetic spectra using a combination of both ISM and stellar components as follows: (1) 10% of stellar contribution at 5 μm + 10% of ISM contribution at 30 μm; (2) 50% of stellar contribution at 5 μm + 50% of ISM contribution at 30 μm; and (3) 100% of stellar contribution at 5 μm + 100% of ISM contribution at 30 μm. The dusty torus spectra are simulated with f(12 μm) ∼ 300 mJy (i.e., the intermediate S/N ∼ 40–60 studied in Sections 4.2 and 4.3). We then repeated the analysis performed in Section 4.2 to study how good is the determination of the parameters.

Note that the simulated fractional contributions of host galaxy (both stellar and ISM) are well recovered with less than 5% error. Figures 8–10 show the average percentage error (compared to the parameter range) per AGN dust model and instrumental setup, after including dilution by the stellar, ISM, and stellar + ISM component in the synthetic spectra, respectively. Our results are robust since tests increasing the number of synthetic spectra do not change them. Better results are obtained when using Spitzer/IRS and JWST/(MIRI + NIRSpec) compared to GTC and JWST/MIRI. However, this is due to the lack of near-infrared coverage of the ISM templates. However, it is worth mentioning that the superb spatial resolution obtained with GTC/CanariCam allows us to isolate the nuclear component from host contributors much better than Spitzer/IRS (see below). Indeed, a large portion of the Swift/BAT sample with Spitzer/IRS spectra is highly contaminated by circumnuclear contributors (see Paper II).

Zoom In Zoom Out Reset image size

We were able to constrain all the parameters (roughly within 15% error) with up to 50% of stellar contribution at 5 μm or 50% of the ISM contribution at 30 μm for all the models (except the [Hoenig17] model and GTC/CanariCam instrumental setup; see Figures 8 and 9). [Nenkova08] and [Hoenig10] can further constrain the parameters with less than 15% error for 100% of the stellar component at 5 μm. [Sieben15] is able to constrain the parameters with less than 15% uncertainty even for 100% of the stellar component at 10 μm. [Hoenig17], due to a large number of parameters, shows the largest percentage errors even for 10% of the stellar component at 5 μm or 10% of the ISM component at 30 μm. The parameters are still within 15% uncertainty only for [Nenkova08], [Hoenig10], and [Sieben15] (excluding in this case GTC/CanariCam) when the stellar and ISM contributions are combined up to 10%–50% (see Figure 10). The parameters for [Fritz06] and [Stalev16] could also be estimated including up to 50% of combined stellar and ISM contributions except for σ and γ for [Fritz06] and p for [Stalev16] (with all the instrumental setups except GTC/CanariCam).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The impact of galaxy dilution on the parameter estimate is stronger than the instrumental setup. In general, 100% of the stellar component at 10 μm doubles the error on the resulting parameters, while 100% of the stellar component at 5 μm introduces a 50% additional error on the resulting parameters. Similarly, 100% of the ISM component at 30 μm doubles the error obtained in the resulting parameters. The combination of 100% of the stellar + ISM components includes twice the error found with 100% of the stellar component and similar error to that found when including 100% of the ISM component. Thus, infrared high spatial resolution spectra are key to studying the AGN dust by restricting the host galaxy contribution to the lowest level. JWST/MIRI or ground-based GTC/CanariCam spectra are very useful to estimate the torus parameters because they better isolate the AGN from the host galaxy. This is clearly seen in Figure 9 if we compare, for instance, the error on the parameter estimate of GTC/CanariCam or JWST spectra and 50% of ISM at 30 μm (data points linked with solid lines) with Spitzer/IRS spectra and 100% of ISM at 30 μm (data points linked with long-dashed lines). Thus, high spatial resolution GTC/CanariCam spectra still play an important role in the parameter estimate until JWST is able to give both spectral coverage and spatial resolution.

Radiative transfer models have proven successful in reproducing the infrared SED of AGN (e.g., Fritz et al. 2006; Netzer et al. 2007; Ramos Almeida et al. 2009; Alonso-Herrero et al. 2011; González-Martín et al. 2015, 2017; Ichikawa et al. 2015; Martínez-Paredes et al. 2015, 2017; Netzer 2015; Fuller et al. 2016; Audibert et al. 2017; Hönig & Kishimoto 2017). These models differ on the dust chemical composition, distribution, morphology, or dynamics. We attempt to discuss for the first time how a large variety of these models could be distinguished based on SED fitting. We have found that four out of the six models could be distinguished based on the slope below ∼7 μm, the slope above ∼25 μm, and the silicate features. Indeed, one of the main controversial aspects of the AGN model unification regarding the dust emission is associated with the silicate features. Most type 1 Seyferts exhibit a rather weak emission feature (Hönig et al. 2010), and AGN in general lack deep 9.7 μm silicate absorption features (Hao et al. 2005, 2007). These weak features are naturally produced when large graphite grains are dominating the dust composition (Hönig & Kishimoto 2010). The peak wavelength is above 10.2 μm in about 65% of the silicate emission features, whereas the shift of the 9.7 μm absorption feature is small (see also Nikutta et al. 2009; Hönig et al. 2010; Hatziminaoglou et al. 2015). Nenkova et al. (2002) proposed that these findings might be explained by a clumpy distribution of dusty clouds because after the inclusion of effects as directly illuminated clumps and clouds illuminated by others, these features are never deep (see also Levenson et al. 2006; Spoon et al. 2007; Sirocky et al. 2008). Nevertheless, different sublimation temperatures of the silicate and graphite grains can also produce these weak features using smooth models (Fritz et al. 2006; Schartmann et al. 2008; Hönig et al. 2010). Indeed, we also found that the smooth toroidal model [Fritz06] and the smooth + clumpy toroidal model [Stalev16] produce a deep absorption feature at 9.7 μm compared to the other models (see Figure 4).

The first question that can be tackled with our analysis is whether we can distinguish between smooth and clumpy torus models. Dullemond & van Bemmel (2005) compared their own radiative transfer models of smooth and clumpy tori, finding that, despite the distinct nature of the models, it was not possible to distinguish between them based on the SEDs (similar results are found by Stalevski et al. 2012). Schartmann et al. (2008) also compared smooth and clumpy torus models developed by their own group both based on the radiative transfer code MC3D. They found that the main difference is that a clumpy medium allows the central source to heat the clouds at large radii. The only work comparing models not based on similar prescriptions for the dust composition is presented by Feltre et al. (2012). They compared the smooth model [Fritz06] and the clumpy model [Nenkova08B], finding that even with matched parameters they do not produce similar SEDs. Due to the different chemical composition in [Fritz06] compared to [Nenkova08], the behavior of the silicate features is quite distinctive (Feltre et al. 2012). They can also be distinguished throughout the near-infrared slopes because the clumpy torus model [Nenkova08] is sensitive to a much wider range of near-infrared slopes than the smooth torus model [Fritz06] (see Figure 2). Interestingly, this narrow range of slopes is also found for the two-phase model [Stalev16]. Feltre et al. (2012) suggested that the difference in the near-infrared slopes might be due to the different slope of the primary emission from the accretion disk. However, all the models except [Fritz06] and [Stalev16] show a wider range for the near-infrared slopes. Note, however, that these steep near-infrared slopes found for [Nenkova08], [Hoenig10], [Sieben15], and [Hoenig17] are not required for the AGN Spitzer/IRS spectra analyzed in Paper II, which might indicate an oversampling of some of the parameter spaces.

We can also attempt to answer whether models with similar cloud distributions produce similar SEDs. Using identical cloud distributions, our synthetic spectral analysis shows that the two clumpy models [Nenkova08] and [Hoenig10] can be distinguished by looking to their silicate feature residuals, which are mainly attributed to different cloud compositions (by Ossenkopf et al. 1992; Draine & Lee 1984, for [Nenkova08] and [Hoenig10], respectively). This can also be seen in Figure 4. [Nenkova08] can produce 18 μm silicate features in absorption (i.e., positive values), while [Hoenig10] can only produce these features in emission. After exploring the resulting spectral shapes for different ranges on the parameters, we found that these differences cannot be attributed to unmatched parameter space for these SED libraries. Indeed, Feltre et al. (2012) showed that, even though clumpy and smooth dust models produce different SEDs, most of the differences arise from the model assumptions and not from the dust distribution. This might also explain why [Hoenig10] and [Hoenig17], despite their major differences in overall components and morphology of the distribution of dust, are undistinguishable when it comes to SED fitting, as shown here. This also happens for the parent models [Fritz06] and [Stalev16]. Therefore, it seems that the dust composition needs to be better explored because it might have a major impact on the shape of the silicate features.

The clumpy torus model [Hoenig10] and its disk+wind version [Hoenig17] naturally produce flat near- and mid-infrared slopes (see Figure 2), a lack of steep far-infrared slopes (see Figure 3), and a narrower range of silicate feature strengths (see Figure 4) compared to other models. These main differences are fully consistent with previous results. Hönig et al. (2010) pointed out that the silicate features and the near-infrared slopes are some of the main differences among models (see also Feltre et al. 2012; Stalevski et al. 2012; Ramos Almeida et al. 2014; Hönig & Kishimoto 2017). The smooth torus model [Fritz06] and its two-phase version [Stalev16] can also be distinguished from the others (see above). However, parent models do not always produce the same spectral features. García-González et al. (2017) found that [Hoenig10] produces stronger near-infrared emission and bluer mid-infrared spectral slopes than its previous version.

Finally, another important result of the present analysis is that the parameters for each model can be determined (within less than 15% of the parameter space). This is true for any model except [Hoenig17] for a source with a 12 μm flux above ∼100 mJy and for any instrumental configuration except for GTC/CanariCam (for which sources with 12 μm flux larger than ∼300 mJy are required). However, the isolation of the AGN mid-infrared emission is key; more than 50% of stellar contribution at 5 μm, 50% of ISM contribution at 30 μm, or a combination of 10% of stellar and 10% of ISM can double the error on the parameters. This is the main reason to require high-resolution infrared observations, such as those provided by the future JWST.

We have investigated a set of six SED libraries with the aim of reproducing the dust infrared emission of AGN. We produced synthetic spectra for GTC/CanariCam, Spitzer/IRS, JWST/NIRSpec, and JWST/MIRI instrumental setups and four sensitivities (equivalent to 30 mJy < F_{12 μm} < 10 Jy or 3 < S/N < 150). We fitted them with the set of models and using a general parameterization that includes three slopes and the strength of the silicate features. The main results are as follows:

• 1.  
  Each model can be distinguished from the others based on the chi-squared statistics. The exceptions are the comparison between [Hoenig10] versus [Hoenig17] and [Fritz06] versus [Stalev16]. This is probably due to the fact that [Hoenig17] ([Stalev16]) is based on [Hoenig10] ([Fritz06]), and therefore they used similar prescriptions for the dust composition.
• 2.  
  In general, residuals show that the main differences among the models are the slopes of the spectra below ∼7 μm and above ∼25 μm and the shape of the 10 and 18 μm silicate features. We also show that these models can be distinguished based on the continuum slopes and the silicate feature strengths.
• 3.  
  The parameters can be well determined within 15% error (compared to their parameter space) for all the models using either Spitzer/IRS, JWST/MIRI, or JWST/(MIRI+NIRSpec). This is true except for [Hoenig17], which has a large number of free parameters. In this case the inclusion of JWST/(MIRI+NIRSpec) data and the selection of the brightest AGN at mid-infrared are needed to improve the resulting parameter determination. We do not see any significant improvement by including near- to mid-infrared data for other models. However, this might be relevant for host galaxy decontamination. Slightly worse results are obtained when using GTC N and Q bands.
• 4.  
  Dilution plays an important role in the parameter determination. Either more than 50% of stellar contribution compared to the AGN contribution at 5 μm, more than 50% of ISM compared to the AGN contribution at 30 μm, or a combination of 50% of each component prevents the parameter determination.

Infrared dust models (in the form of either a torus, disk, or wind) could be distinguished and their parameter constrained using the spectral fit to the 5–30 μm wavelength range. However, high spatial resolution data are required to isolate this AGN emission from host galaxy dilution. For that, current ground-based mid-infrared (e.g., GTC/CanariCam) and future JWST data are required. In Paper II, we fit these models to a sample of 110 AGN selected from the Swift/BAT survey with available Spitzer/IRS spectra.

We thank the anonymous referee for his/her comments and suggestions, which have improved significantly the results of this research. This research is mainly funded by the UNAM PAPIIT project IA103118 (PI OG-M). M.M.-P. acknowledges support by KASI postdoctoral fellowships. I.M. and J.M. acknowledge financial support from the research project AYA2016-76682-C3-1-P (AEI/FEDER, UE). J.M.R.-E. acknowledges support from the Spanish Ministry of Science under grants AYA2015-70498-C2-1 and AYA2017-84061-P. I.G.-B. acknowledges financial support from the Spanish Ministry of Science and Innovation (MICINN) through projects PN AYA2015-64346-C2-1-P and AYA2016-76682-C3-2-P. I.M. and J.M. acknowledge financial support from the State Agency for Research of the Spanish MCIU through the "Center of Excellence Severo Ochoa" award for the Instituto de Astrofísica de Andalucía (SEV-2017-0709). D.E.-A. acknowledges support from a CONACYT scholarship. D.-D. acknowledges PAPIIT UNAM support from grant IN113719. This research has made use of dedicated servers maintained by Jaime Perea (HyperCat at IAA-CSIC); Alfonso Ginori González, Gilberto Zavala, and Miguel Espejel (Galaxias, Posgrado04, and Arambolas at IRyA-UNAM); and Daniel Díaz-González (IRyAGN1 and IRyAGN2). All of them are gratefully acknowledged.

Figures 11 and 12 show the residuals of a ${\chi }^{2}$ distribution for each set of simulated spectra (as in Figure 7 of Nenkova et al. 2008b).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size