Three-dimensional radiative transfer modeling of AGN dusty tori as a clumpy two-phase medium
We investigate the emission of active galactic nuclei (AGN) dusty tori in the infrared domain. Following theoretical predictions coming from hydrodynamical simulations, we model the dusty torus as a 3D two-phase medium with high-density clumps and low-density medium filling the space between the clumps. Spectral energy distributions (SED) and images of the torus at different wavelengths are obtained using 3D Monte Carlo radiative transfer code SKIRT. Our approach of generating clumpy structure allows us to model tori with single clumps, complex structures of merged clumps or interconnected sponge-like structure. A corresponding set of clumps-only models and models with smooth dust distribution is calculated for comparison. We found that dust distribution, optical depth, clump size and their actual arrangement in the innermost region, all have an impact on the shape of near- and mid-infrared SED. The m silicate feature can be suppressed for some parameters, but models with smooth dust distribution are also able to produce a wide range of the silicate feature strength. Finally, we find that having the dust distributed in a two-phase medium, might offer a natural solution to the lack of emission in the near-infrared, compared to observed data, which affects clumpy models currently available in the literature.
keywords:galaxies: active – galaxies: nuclei – galaxies: Seyfert – radiative transfer.
According to the unification model, the different appearance of type 1 and type 2 active galactic nuclei (AGN) is only a matter of orientation (Antonucci, 1993; Urry & Padovani, 1995). This hypothesis postulates a supermassive black hole and its accretion disc as the central engine. The accretion disc is the source of the strong X-ray emission and UV/optical continuum. Superimposed on the continuum are the broad emission lines, coming from gaseous clouds moving at high velocities, the so-called broad-line region (BLR). This central region is surrounded by a geometrically and optically thick toroidal structure of dust and gas with an equatorial visual optical depth much larger than unity. Viewed edge-on, this dusty torus blocks the radiation coming from the accretion disc and BLR. In this case an UV/optical bump and broad emission lines are absent and the object appears as a type 2 AGN. In the case when the line of sight does not cross the dusty torus, both accretion disc and BLR are exposed, giving rise to a strong UV/optical continuum and broad emission lines, and the object is classified as a type 1 active galaxy. Observed polarized nuclear emission in type 2 sources, scattered by electrons and tenuous dust (Antonucci & Miller, 1985; Pier et al., 1994; Packham et al., 1997), supports the unification model. It proves that an active galactic nucleus is present even though no direct emission from accretion disc is observed. The toroidal geometry also explains several other observables such as, the presence of ionizing cones (Pogge, 1988, 1989; Tadhunter & Tsvetanov, 1989) and high hydrogen column densities in the X-ray, usually associated with type 2 sources (e.g. Shi et al., 2006).
The dusty torus is expected to absorb the incoming accretion disc radiation and re-emit it in the infrared domain. Thus, further support for the existence of a dusty torus comes from the observed mid- to far-infrared bump and silicate feature at m in the spectral energy distribution (SED) of AGNs. In type 1 sources, hot dust in the inner region can be seen directly and the feature is expected to be detected in emission. Recent mid-infrared observations obtained with the Spitzer satellite confirm the silicate emission feature in AGNs (Siebenmorgen et al., 2005; Hao et al., 2005). In type 2 objects, the dust feature is observed in absorption (e.g., Jaffe et al., 2004) due to obscuration by the cold dust. Additional evidence of dusty tori comes from the interferometric observations. Using speckle interferometry, the nucleus of NGC 1068 was successfully resolved in the -band (Wittkowski et al., 1998) and in the -band (Weigelt et al., 2004). This resolved core is interpreted as dust close to the sublimation radius. Tristram et al. (2007) reported VLTI interferometric observations with strong evidence of a parsec scale dust structure in the Circinus galaxy. Kishimoto et al. (2011) reported indication of a partial resolution of the dust sublimation region in several type 1 AGNs using the Keck interferometer.
In order to prevent the dust grains from being destroyed by the hot surrounding gas, Krolik & Begelman (1988) suggested that the dust in the torus is organized in a large number of optically thick clumps. However, due to the difficulties in handling clumpy media and lack of computational power, early work was conducted by adopting a smooth dust distribution. Several authors explored different radial and vertical density profiles (Pier & Krolik, 1992, 1993; Granato & Danese, 1994; Efstathiou & Rowan-Robinson, 1995; van Bemmel & Dullemond, 2003; Schartmann et al., 2005; Fritz et al., 2006). The first effort of developing the formalism for the treatment of clumpy media was undertaken by Nenkova et al. (2002, 2008a, 2008b). They utilized a 1D radiative transfer code to compute the SED of a single irradiated clump. In a second step a statistical generalization is made to assemble the SED of the torus. They claim that only clumpy tori are able to reproduce the observed properties of the silicate feature. However, Dullemond & van Bemmel (2005) modeled the torus as a whole, using a 2D radiative transfer calculations. They adopted a geometry with axial symmetry and modeled clumps in the form of rings around the polar axis. They made a direct comparison of models with clumps and corresponding smooth dust distributions and concluded that there is no evidence for a systematic suppression of the silicate emission feature in the clumpy models. Hönig et al. (2006), with an upgrade of their model in Hönig & Kishimoto (2010), adopted a similar method as Nenkova et al. (2002), but they employed a 2D radiative transfer code and took into account various illumination patterns of clumps. They reported that a suppression of the silicate emission feature strongly depends on the random distribution and density of the clumps in the innermost region. More recently Schartmann et al. (2008) presented a 3D radiative transfer models of clumpy tori. Their findings are in agreement with those by Hönig et al. (2006). On the other hand, using the models developed by Fritz et al. (2006), Feltre et al. (2011) found that a smooth distribution of dust is also capable of reproducing the observed variety of the silicate feature strength. A further study of the silicate feature and its properties, such as the appearance in emission in some type 2 objects and the apparent shifting toward long wavelengths in some objects, is presented in Nikutta et al. (2009).
A problem which the obscuring torus hypothesis faced from the beginning was formation of the dynamically stable structure and maintenance of the required scale-height. Krolik & Begelman (1988) presented a scenario according to which the scale-height is maintained through regular elastic collisions between the clumps (see also Beckert & Duschl, 2004). In the case of a continuous dust distribution, Pier & Krolik (1992), followed by Krolik (2007), suggested that radiation pressure within the torus may be enough to support the structure. Wada & Norman (2002) (with model update in Wada et al., 2009) performed a 3D hydrodynamical simulations, taking into account self-gravity of the gas, radiative cooling and heating due to supernovae. They found that such a turbulent medium would produce a multiphase filamentary (sponge-like) structure, rather then isolated clumps. A scenario where the effects of stellar feedback from a nuclear cluster play a major role is discussed in Schartmann et al. (2009).
In this paper we present our modeling of 3D AGN dusty tori. We model the torus as the whole, using the 3D Monte Carlo radiative transfer code SKIRT. We take a step further toward a more realistic model by treating the dusty torus as a two-phase medium, with high density clumps and low density medium filling the space between them. We calculated SEDs and images of the torus for a grid of parameters. Our approach allows us to, for each two-phase model, generate a clumps-only model (with dust distributed to the clumps exclusively, without any dust between them) and a smooth model with the same global physical parameters. Our aims are (a) to investigate the influence of different parameters on model SEDs and their observable properties, (b) to put to a test reports that the observed SEDs in the mid-infrared domain unambiguously point to a clumpy structure of dusty tori; if that is indeed the case, a comparison of clumpy and smooth models should show a systematic difference of their observable properties, such as the strength of the silicate feature.
The paper is organized as follows. In Section 2 we give the description of the radiative transfer code and present the physical details of our model. In Section 3 we discuss how different parameters affect modeled SEDs, analyze their observable properties and compare two-phase, clumps-only and smooth models. Finally, in Section 4 we outline our conclusions.
2.1 The radiative transfer code
We have used the radiative transfer code SKIRT (Baes et al., 2003, 2011) for the modelling of AGN dusty tori. SKIRT is a 3D continuum radiative transfer code based on the Monte Carlo algorithm (Cashwell & Everett, 1959; Witt, 1977), initially developed to study the effect of dust absorption and scattering on the observed kinematics of dusty galaxies (Baes & Dejonghe, 2001, 2002; Baes et al., 2003). It has been extended with a module to self-consistently calculate the dust emission spectrum under the assumption of local thermal equilibrium – LTE (Baes et al., 2005a). This LTE version of SKIRT has been used to model the dust extinction and emission of galaxies and circumstellar discs (Baes et al., 2010; De Looze et al., 2010; Vidal & Baes, 2007).
Similar to most modern Monte Carlo codes (e.g. Gordon et al., 2001; Wolf, 2003; Niccolini, Woitke, & Lopez, 2003; Bianchi, 2008), the SKIRT code contains many deterministic elements which makes the Monte Carlo technique orders of magnitude more efficient than the naive Monte Carlo recipe. These include the peeling-off technique (Yusef-Zadeh, Morris, & White, 1984), continuous absorption (Lucy, 1999; Niccolini, Woitke, & Lopez, 2003), forced scattering (Cashwell & Everett, 1959; Witt, 1977) and smart detectors (Baes, 2008). For the current simulations, we use the technique of frequency distribution adjustment presented by Bjorkman & Wood (2001) and critically discussed by Baes et al. (2005b). This technique ensures that at each moment during the simulation, the wavelength distribution from the photon packages emitted from the cell are in agreement with the cell’s current temperature. The main advantage of this technique is that no iteration is necessary.
2.2 Model description
Dust distribution and properties
We approximate the spatial dust distribution around the primary continuum source with a conical torus (i.e. a flared disc). Its characteristics are defined by (a) the half opening angle – defining also the maximum height extent to which the dust is distributed –, (b) the inner and outer radius, and respectively, and (c) the parameters describing dust density distribution, and (see below). A schematic representation of the adopted geometry is given in Fig. 1. For the inner radius of the dusty torus we adopted the value of pc, at wich the dust grains are heated to the temperature of K, according to the prescription given by Barvainis (1987):
where is the bolometric ultraviolet/optical luminosity emitted by the central source, expressed in units of erg s and is the sublimation temperature of the dust grains given in units of 1500 K.
We describe the spatial distribution of the dust density with a law that allows a density gradient along the radial direction and with polar angle, as the one adopted by Granato & Danese (1994):
where and are coordinates in the adopted polar coordinate system (see Fig. 1).
The dust mixture consists of separate populations of graphite and silicate dust grains with a classical MRN size distribution (Mathis, Rumpl & Nordsieck, 1977):
where the size of grains, , varies from to m for both graphite and silicate. The normalization factors for size distribution are and for graphite and silicate, respectively (Weingartner & Draine, 2001). Optical properties are taken from Laor & Draine (1993) and Li & Draine (2001).
The dust is distributed on a 3D cartesian grid composed of a large number of cubic cells. The dust density is constant within each cell. The standard resolution for our simulations is 200 cells along each axis ( cells in total). However, to properly sample the dust properties, an increase in the torus size requires an increase of the resolution of the computational grid as well. Thus, to simulate a torus twice the size of the ‘standard model’, one needs to employ a grid with cells along each axis, that is, cells in total.
In the case of a smooth density distribution, the axial symmetry in our model reduces the 3D radiative transfer computations to a 2D problem, with a significant gain both on the computational time and the memory usage. However, the prescription we use for generating clumpy models demands a 3D cartesian grid. Therefore, such a grid was used throughout all our simulations, in order to avoid any possible differences due to the adoption of different grids. To ensure that the adopted resolution is sufficient to properly sample the dust, for each simulation we compare the actual and expected values of (a) the face-on and edge-on central surface density and (b) the total dust mass.
The emission for all models was calculated on an equally spaced logarithmic wavelength grid ranging from to m. A finer wavelength sampling was adopted between and m, in order to better resolve the shape of and m silicate features. Each simulation is calculated using photon packages.
Spectral energy distribution of the primary source
The primary continuum source of dust heating is the intense UV-optical continuum coming from the accretion disc. A very good approximation of its emission is a central, point-like energy source, emitting isotropically. Its SED is very well approximated by a composition of power laws with different spectral indices in different spectral ranges. The adopted values are:
and the resulting SED is plotted on Fig. 2. These values have been quite commonly adopted in the literature, and come from both observational and theoretical arguments (see e.g., Schartmann et al., 2005). We have anyway verified that changes in the shape of the primary source SED affect only very marginally the infrared re-emission. For the bolometric luminosity of the primary source we adopted the value of (see e.g., Davis & Laor, 2011).
As mentioned above, an isotropic emission of the central source is commonly adopted in the literature (e.g. Schartmann et al., 2005; Hönig et al., 2006; Nenkova et al., 2008a). However, the disk emission is inevitably anisotropic (see, for example, Kawaguchi & Mori, 2011, and references therein). Therefore, we have performed additional calculations assuming anisotropic radiation of the central source, in order to investigate the resulting influence on the model SEDs. Radiation flux () from a unit surface area of an optically thick disk toward a unit solid angle at the polar angle of decreases with an increasing according to the formula given by Netzer (1987):
where the first term represents the change in the projected surface area, and the the latter represents the limb darkening effect. For simplicity, in our calculations we take into account only the first term. As the accretion disk radiation varies with , the dust sublimation radius also depends on it:
where is the dust sublimation radius estimated assuming isotropic emission. As a result, the inner edge of the torus is (a) closer to the central source than suggested by the Eq. 1 and (b) the structure of the edge is concave (Kawaguchi & Mori, 2010, 2011). Also, Kawaguchi & Mori (2010) found that the dust sublimation radius can decrease down to , all the way to the outermost radius of accretion disk. However, due to the numerical constraints, that is, the minumum size of the dust cell in the computational grid we are currently limited to, we cannot allow the dust to be placed all the way to the primary source. Instead, we are forced to put a limit on the minimum allowed dust sublimation radius at pc (). We discuss the influence of the anisotropic central source radiation on the dusty tori model SEDs in the Sec. 3.4; throughout the rest of the paper, the isotropic emission is assumed.
Two-phase medium: the approach
Models of emission in which the dust is uniformly, smoothly distributed within the toroidal region are obtained starting from Eq. 2. The density gradient is determined by the two parameters, and . The total amount of dust is fixed based on the equatorial optical depth at m (). While for the smooth models distributing the dust is straightforward, for the clumpy model this process is non-trivial. As hydrodynamical simulations of Wada & Norman (2002) demonstrated, dust in tori is expected to take the form of a multiphase structure, rather than isolated clumps. Therefore, we adopted the approach which allows us to generate such a medium. We start from the corresponding smooth models, i.e. the one with the same global parameters, and then apply the algorithm described by Witt & Gordon (1996) to generate a two-phase clumpy medium. According to this algorithm, each individual cell in the grid is assigned randomly to either a high- or low-density state by a Monte Carlo process. The medium created in such a way is statistically homogeneous, but clumpy. The filling factor determines the statistical frequency of the cells in the high-density state and can take values between and . For example, a filling factor of represents a case of rare, single high-density clumps in an extended low-density medium. The process for the clump generation is random with respect to the spatial coordinates of the individual clumps themselves. Thus, as the filling factor is allowed to increase, the likelihood that adjoining cells are occupied by clumps increases as well. This leads to the appearance of complex structures formed by several merged clumps. For filling factor values , clumps start to form an interconnected sponge-like structure, with low-density medium filling the voids. We form larger clumps by forcing high-density state into several adjoining cells.
To tune the density of the clumps and the inter-clump medium, we use the ‘contrast parameter’, defined as the ratio of the dust density in the high- and low-density phase. This parameter can be assigned any positive value. For example, setting the contrast to one would result in a smooth model. Setting extremely high value of contrast () effectively puts all the dust into the clumps, without low-density medium between them. An example of dust density distributions at the meridional plane for smooth, two-phase and clumps-only models is given in Fig. 3.
2.3 Parameter grid
In this section we present the adopted values of parameters used to calculate a grid of models for our analysis.
For the inner and outer radius of the torus we adopted the values and pc, respectively. The half opening angle of the torus – – is fixed to for all of our model realizations. All models are calculated at 7 different line-of-sight inclinations, namely , , , , , and , where represents a face-on (type 1) AGN and an edge-on (type 2) AGN. Inclinations between and (dust-free lines of sight) were omitted since their SED shows no appreciable difference with respect to the full face-on view. The equatorial optical depth takes the values , , , . We note here that this is the optical depth of the initial smooth model, before the dust is redistributed to make the clumpy one (see Sec. 2.2.3). Thus, the exact value along the given line of sight will vary due to the random distribution of the clumps. The parameters defining the spatial distribution of the dust density (Eq. 2) are , and , , , .
Defining the relative clump size, , as the ratio of the outer radius of the torus over the clump size:
we explored the clumpy models for two different clump sizes, pc and pc, that is, and , respectively. The number of clumps in the former case is , and each clump occupies one grid cell. In the latter case there are clumps, each one being composed of grid cells. The adopted filling factor values are in the case , and in the case models. Both values allow us to have single, as well as clusters of high-density clumps immersed into a low-density medium. The contrast between high- and low-density phases is fixed at .
We generated three sets of models with the same global physical parameters: (a) models with the dust distributed smoothly, (b) models with the dust as a two-phase medium and (c) models with a contrast parameter set to an extremely high value (), effectively putting all the dust into the high-density clumps. We will refer to these models as ‘smooth’, ‘two-phase’ and ‘clumps-only’, respectively. For clumpy models (both two-phase and clumps-only) we generated another set of models with the same parameters, but with a different random distribution of clumps.
For each model we calculated SEDs and images of torus for all the points in the wavelength grid. Calculated flux is scaled for the torus distance of 10 Mpc from the observer. The parameter grid is summarized in Table 1.
|, , ,|
|, , ,|
|Size of clumps||pc, pc|
|Inclination||, , , , , ,|
3 Results and discussion
In this section we discuss the influence of different parameters on the general shape of the SEDs of the two-phase models and analyze their observable properties. The following analysis refers to the two-phase models with . We will address how the properties of these models compare to properties of models with , clumps-only, and smooth models in Sections 3.9 and 3.10.
3.1 SED dependence on the viewing angle
As it was demonstrated in early works (e.g. Granato & Danese, 1994), the SED of a dusty torus depends on the viewing angle. In Fig. 4 we show the SED dependence on the inclination of the torus. There is a clear distinction between cases of dust-free lines of sight (, ) and those that pass through the torus (, , ). For the adopted value of half opening angle (), this transition occurs at inclination . The most notable difference is found shortward of m. In the case of dust-free lines of sight, we directly see the radiation coming from the accretion disc, while in the case of dust-intercepting paths most of the radiation is absorbed and scattered at different wavelengths. This is especially pronounced in those systems where the density remains constant with polar angle. In the case of a non-constant density, the transition from face-on to edge-on view is smoother. An exception is the case of very low optical depths, where it is possible to directly ‘see’ the central source even when viewed edge-on. Another important difference between dust-free and dust-intercepting lines of sight is the m silicate feature, which is expected to appear in emission in the former case and in the absorbtion in the latter. The properties of this feature is discussed in detail in Section 3.5. Images of the torus at different wavelengths, for the type 1 and type 2 lines of sight are shown in Fig. 5. From the figure it is clear that, at shorter wavelengths, it is the radiation from the inner (and hotter) region that dominates. At longer wavelengths, the emission arises from the dust placed further away. Thus, the size of torus is wavelength dependent. In Fig. 6 we present the total SED and its thermal and scattered components, along with the primary source emission, for face-on and edge-on view. As it can be seen from the figure, the thermal component is dominant in the mid- and far-infrared part of SED and its shape is similar for both type 1 and type 2 orientations. The shape and amount of scattered component is quite different; in the edge-on view it determines the total SED at shorter wavelengths, while in the face-on view it is negligible compared to the primary source emission.
3.2 SED dependence on the filling factor and contrast
As described in Section 2.2.3, the two parameters that determine the characteristics of the two-phase medium are the filling factor and the contrast. The filling factor determines the percentage of grid cells in a high-density state. Models with low values for the filling factor (e.g. ) represent systems with rare, single high-density clumps in extended low-density medium. As the filling factor increases, the number of clumps will increase as well, forming clusters of clumps, or even single, interconnected sponge-like structure. This is illustrated in Fig. 7, where we show dust density distributions at the meridional plane for different filling factors. Fig. 8 shows SEDs of models for different filling factors, compared with SED of the corresponding smooth model. From this Fig. we see that, as the filling factor increases, the overall mid-infrared emission increases as well. For filling factor value of , in the face-on view, the silicate feature is attenuated. As the filling factor increases, clumps start to form sponge-like structures, more and more resembling a smooth dust distribution, and the strength of the silicate feature approaches the strength in the corresponding smooth model. In the edge-on, a filling factor value of produces silicate features in weaker absorption than in the corresponding smooth model. As the filling factor increases, the strength of the silicate feature approaches the strength of the feature in the smooth models.
The ‘contrast’ parameter sets the density ratio between the high- and low-density phases. Fig. 9 shows the model SED dependence on this parameter. In the face-on view, for increasing contrast, both the hot dust emission ( m) and the strength of the silicate feature decrease. From the same figure we also see that, for higher contrast values, the peak of the silicate feature is slightly shifted toward longer wavelengths. In the edge-on view, the silicate feature in absorption gets weaker with increasing contrast.
3.3 SED dependence on the random distribution of clumps
The shape and overall near- and mid-infrared emission strongly depend on the distribution of dust in the innermost region. Changing the random arrangement of clumps, along with choosing a particular line of sight, can affect the resulting SED significantly, as illustrated in Fig. 10. As described in Section 2.2.3, the process for clump generation is random with respect to the spatial coordinates of the individual clumps. As a consequence, adjoining cells can be occupied by individual clumps, forming complex structures of several connected clumps. In models with a higher concentration of clumps in the innermost region, due to the shadowing effect, the absorption is increased and silicate feature is suppressed.
This characteristic imports a degree of degeneracy in the features of the SEDs, which will be less directly dependent on the physical input parameters. Even though the spatial position of the clumps is not related to the physical properties of dusty tori, their re-arrangement has a clear impact on the infrared emission. It is, in some way, mimicking a change in the optical depth, which might appear either to increase or decrease, depending on the clumps re-arrangement, especially in the innermost regions.
Some random arrangements of the clumps have interesting repercussions. Because of clumpiness, the difference between the SED of type 1 and 2 objects is not truly an issue of orientation; it is rather a matter of probability of directly viewing the main energy source of the AGN (Nenkova et al., 2008b). As a result, type 1 sources can be detected even from what are typically considered as type 2 orientations. Such a scenario provides an explanation for the few Seyfert galaxies with type 1-like optical spectra whose m SED resembles that of a type 2 AGN (Alonso-Herrero et al., 2003). Conversely, if a clump happens to obscure the central engine from an observer, that object would be classified as type 2 irrespective of the viewing angle. In such cases, the clump may move out of the line-of-sight, creating a clear path to the nucleus and a transition to a type 1 spectrum. Such transitions between type 1 and type 2 line spectra have been observed in a few sources (see Aretxaga et al., 1999, and references therein).
3.4 Influence of anisotropic primary source radiation on model SED
As described in Sec. 2.2.2, an isotropic source emission is commonly adopted in the radiative transfer modeling of dusty tori; however, the accretion disk emission is actually anisotropic. In this section, we discuss the influence of anisotropic source radiation on the model SEDs. Dependence of the accretion disk radiation on the direction is taken according to Eq. 5 and the corresponding change of the dust sublimation radius according to Eq. 6. In Fig. 11 we present the resulting model SEDs if the anisotropic radiation of the primary source is assumed (dotted line) and compare them to the corresponding SEDs obtained in the case of the isotropic source (solid line) for different inclinations. SEDs were calculated for the inclinations between and with the step of ; for the clarity of the Figure, only SEDs for three inclinations, , and are shown.
We found that, when anisotropy of the central source is assumed, the IR SED can indeed change, resulting in a lower emission, though roughly keeping the same shape. This is a logical consequence coming from the fact that, for a given bolometric luminosity of the accretion disc, an anisotropic source whose characteristics are those as described above, is emitting more power in the dust-free region: the overall result is a less luminous torus. The excess of emission shortward of m is seen in the dust-free lines of sight, because, at these wavelengths the primary source contribution is still significant. Again, as expected from the properties of radiative transfer (Ivezić & Elitzur, 1997), we found that the shape and the features of the SED (e.g. the m feature) are not affected. Therefore we conclude that our analysis in the rest of the paper is not affected by the isotropic approximation for the central source radiation.
3.5 The silicate feature strength
As it was mentioned above, an important characteristic in the infrared part of an AGN SED is the so-called silicate feature. This silicate feature is caused by Si-O stretching modes, giving rise to either emission or absorption band, peaking at m. All of the early models were dealing with the following issue: while they were properly predicting it in absorption in type 2 objects –in agreement with what was indeed observed–, observations from that period were not supporting the models’ prediction of a silicate feature in emission in type 1 AGN. In fact, one of the main issues driving the development of clumpy models, aimed at addressing this discrepancy between models and observations. Later observations performed by Spitzer with its infrared spectrometer IRS, showed that for a number of type 1 object this feature is indeed observed in emission, partially solving this issue. Recently, Hony et al. (2011) reported the detection of a very strong m feature in emission. On the other hand, Fritz et al. (2006) showed that smooth models are also able to properly reproduce the observed emission in this range. Furthermore, the comparative study performed by Feltre et al. (2011) showed that clumpy and smooth dust distributions are equally able to reproduce both observed broad-band SEDs and mid-infrared Spitzer spectra.
The strength of the m feature can be characterized by the dimensionless parameter , the natural logarithm of the peak-over-continuum ratio (Pier & Krolik, 1992; Granato & Danese, 1994). The continuum is defined by a power law connecting the fluxes at and m. assumes positive values for a feature in emission and negative ones if it is in absorption.
In a face-on view, takes values in the range . The silicate feature is present in a strong emission in the models with lower optical depths (). Models with an optical depth of are showing a wider range of intensities, most of them of moderate strength, with a few cases of strong or weak emission. The strength of the feature in models with high optical depth () is also showing a wide range of intensities, but with overall lower values, and is significantly attenuated in some cases.
For the majority of the edge-on models with optical depths of and the silicate feature is in absorption, with . Models with low optical depths ( and ) do not provide enough dust to absorb the silicate feature and, in this case, it is present in emission even in the edge-on view. Fig. 12 shows SED dependence on the optical depth. To further illustrate dependence of the strength of the silicate feature on the different parameters, in Fig. 13 we plot its intensity, , as a function of the optical depth (), of the dust distribution parameters ( and ) and of the clumps size . For these calculations, the following values of the parameters were chosen: , , , and then each of these parameters was varied while all the others were kept constant.
3.6 SED width
Following Pier & Krolik (1992) and Granato & Danese (1994), the width of the SED, , is defined as the logarithmic wavelength interval in which the power emitted in the infrared is more than one third of the peak value. For a black body this parameter has a value of , while in the observed spectra its value is always larger than . The vast majority of model SEDs, both in the face-on and edge-on views, have a width spanning the range. SEDs with widths are produced by models with optical depths of and because (a) these are the models with the larger amounts of dust and (b) the high values of the optical depth provide a better shielding of the primary source, allowing colder dust temperatures causing, in turn, the broadening of the SED. A small fraction of model SEDs have . These widths are almost exclusive to models with optical depth of and a density law parameter , which produce the silicate feature in very strong emission. Since the maximum of the infrared emission often coincides with the peak of the silicate feature, such models produce lower values.
Another parameter that affects the SED width is the size of the torus. Increasing the radius, while keeping the optical depth constant, means that the amount of dust in the outer (and colder) regions increases. As these regions emit in the far-infrared, an increase in the radius makes the SED wider. For the same reason, will increase with the total amount of dust, that is, with the optical depth (see Fig. 12).
The edge-on orientations produce wider SEDs than the face-on ones, with almost 50% of them having . This is because in the edge-on view the silicate feature is usually in absorption. As a result, the peak of the infrared emission decreases, leading to a wider SEDs. Furthermore, in the edge-on view the received radiation is mainly coming from the outer regions that contribute to the far-infrared emission.
3.7 Isotropy of the infrared emission
Following Dullemond & van Bemmel (2005) we define the isotropy parameter, , as the ratio of the total integrated infrared flux in an edge-on view over the total integrated infrared flux in a face-on view. Larger value of implies there is more isotropy.
Anisotropy in the infrared emission is expected in all systems with torus-like geometry. This is because in the face-on view we have a direct view of the primary source and the inner, hotter region of the torus, while in the edge-on view they are obscured. The values of strongly depend on the optical depth. Models with a low optical depth are almost perfectly isotropic: models with produce and in models with takes values around . Models with a higher optical depth have anisotropic emission, with most values being around and for optical depths of and , respectively. The lowest value in our models is .
3.8 The peak of the infrared emission
Another important feature characterizing the infrared SED of AGNs is the wavelength at which it peaks. We measure this quantity in our model SEDs expressed in . The majority of the models in our grid peak around m, more or less corresponding to centre of the silicate band. A small fraction of models has its maximum in the to m range: all of these models have a high optical depth (either 5 or 10). In the face-on view, almost all models peak at m. In the edge-on view, models exhibiting the silicate feature in emission (i.e. models with low optical depths of ) also peak at m, due to the prominence of the m feature in emission in low-optical depth systems and lower dust content. Edge-on models with higher optical depths peak beyond m.
3.9 Comparison of two-phase and smooth models
In this section we investigate the differences between the models with homogeneous dust distribution (smooth models) and models with dust as a two-phase medium. In order for this comparison to be as consistent as possible, for each two-phase model we have generated its corresponding smooth configuration, using the same global physical parameters. Furthermore, we have generated two different sets of two-phase models using a relative clump size (see Eq. 7) value of and , respectively: in the latter case, the clumps are eight times bigger than in the former.
We found that two-phase models with (small clumps) tend to have a less pronounced emission in the m range, when compared to the smooth ones. If we compare the intensity of the m silicate feature, we find virtually no difference in type 1 view, while slightly lower absorption are measured in two-phase SEDs for type 2 lines of sight. As expected, the dust distributed in a large number of small clumps, embedded in a smooth, homogeneous medium, will closely resemble the characteristics of a smooth SED.
Two-phase models with bigger clumps () are showing more difference compared to both smooth and models. In the face-on view, they tend to have even less pronounced emission and also a different, flatter, slope in the m range. Depending on the parameters of the dust distribution and on the optical depth, the silicate feature is in general less pronounced. This behavior can be attributed to the shadowing effect caused by the clumps in the innermost region, where the dust is hotter and the feature is produced. In the edge-on view, the silicate band absorption is less deep compare to both smooth and -models because we are able to penetrate – at least partly – between the clumps deeper into the torus. Fig. 14 presents a comparison of SEDs of typical models with the smooth and the two two-phase dust distributions for the two clumps sizes.
For face-on view, although both smooth and two-phase models are able to produce almost the same range of values of the silicate feature strength, two-phase models tend in general to produce attenuated emission compared to those produced by the corresponding smooth models. The majority of both smooth and two-phase model SEDs have their infrared emission maximum around m. However, while no smooth model peaks beyond m, there are several two-phase models that peak around m. This is because the two-phase models tend to produce an attenuated silicate emission feature, and when it is very weak or absent, the peak of the emission is shifted toward the longer wavelengths. In the edge-on view, the two-phase models produce a weaker silicate absorption feature, with the lowest strength around . The smooth models produce a deeper silicate feature, with the strength value reaching a minimum of .
Two more characteristics which are of interest when comparing smooth and two-phase models, are the isotropy of the infrared emission and the SED width (see sec. 3.7 and 3.6 for definitions). Both two-phase and smooth models produce a similar range of values of the isotropy parameter . However, compared individually, two-phase models are more isotropic than the smooth ones. Regarding the SED width, , we found that clumpiness does not have a profound effect on this parameter.
In Fig. 15 we present plots of SEDs covering our standard parameter grid, for three characteristic inclinations (, , ). This figure illustrates how SEDs of smooth, two-phase and clumps-only models compare to each other and evolve with the different parameters, i.e. inclination, optical depth and the two parameters determining the dust distribution. In models with low optical depth, the silicate feature appears in a strong emission and the difference between smooth and clumpy models is marginal. With increasing optical depth the difference is increasing as well. Also, the difference between smooth and clumpy models is greater in the cases of constant dust density with polar angle ( in Eq. 2), and non-constant dust density in the radial direction ().
3.10 Comparison of two-phase and clumps-only models
As it can be seen from Fig. 15, the major difference between SEDs of two-phase and clumps-only models arises in the near-infrared range and mainly for face-on view. At these wavelengths, most of the two-phase models with type 1 inclination have a flatter SED when compared to the corresponding clumps-only models. This difference is caused by the presence of the smooth component in which the clumps are embedded. Dust in this component, exposed to the radiation field of the central source, can reach high temperatures and will give rise to higher luminosity in the m range.
Regarding the m silicate feature, we do not find any significant difference between the two dust configurations: depending on the parameters, in clumps-only models it could be slightly attenuated compared to the one in the two-phase models, but the difference is in the most cases marginal. A similar behaviour can be observed in SEDs of edge-on views, in which the smooth low density component is responsible for additional absorption, so the silicate feature is slightly deeper in the two-phase models. The dissimilarities between the SEDs in these two dust configurations increase as the optical depth increases: from models with the lowest value (), where the SEDs are identical, to models with the highest value () which display the most evident differences. The difference is the most pronounced in the cases of constant dust density with polar angle (), and non-constant dust density in the radial direction ().
It is very interesting to note that such a behaviour of the near- and mid-infrared SED of the two-phase dust distribution, would overcome an issue that seems to be common to the most clumpy models currently available in the literature. Exploiting the model of Nenkova et al. (2008b), Mor et al. (2009) fit a sample of mid-infrared spectra of 26 luminous quasar, finding the need of an extra hot-dust component, which they add to the clumpy torus SED, in order to properly reproduce the shorter wavelengths part of the Spitzer spectrum. The addition of this hot dust, whose emission is represented by a black-body with a temperature of about the sublimation limit of graphite, is required by the lack of emission from the adopted clumpy model at these wavelengths. More recently Deo et al. (2011), adopted the same clumpy model to reproduce a combination of observed broad-band photometry and the mid-infrared spectrum of 26 high redshift type-1 quasars. Similarly to Mor et al. (2009), the adopted clumpy models are not able to simultaneously reproduce the intensity of the silicate feature and the near-infrared continuum emission: models that would properly fit the continuum were overestimating the silicate feature emission. An analogous problem was also spotted by Vignali et al. (2011), when using the same clumpy models to fit the observed photometry and IRS spectrum of a type-2 quasar. Adopting the clumpy models developed by Hönig et al. (2006), to fit both photometry and mid-infrared spectroscopy data, Polletta et al. (2008) reach similar conclusions.
As we have shown, a torus model with the dust distributed in a two-phase medium, has a more pronounced (‘hotter’) emission in the m range while displaying, at the same time, a silicate feature whose intensity is almost identical to that of the corresponding clumps-only model.
3.11 Other results in the literature
Making a detailed comparison of our modeling approach with models previously developed in the literature and their results, is quite a tricky task, and is beyond the scope of our work. Furthermore, what we describe in this paper is a model which would ideally put itself in between smooth and clumpy models approach, and it is hence not directly comparable to any of formerly published work. In this section we give a very brief description, which is by no means meant to be exhaustive, of some of the aforementioned works, limiting ourselves to models that consider a clumpy dust distribution.
The exploitation of radiative transfer codes to model AGN IR emission, taking into account the clumpy nature of dust surrounding the central source, includes at present only few works: Nenkova et al. (2002, 2008a, 2008b), Dullemond & van Bemmel (2005), Hönig et al. (2006); Hönig & Kishimoto (2010), Schartmann et al. (2008) and Kawaguchi & Mori (2011). Each of these works exploits different techniques and approximations.
In their series of works, Nenkova et al. used the radiative transfer code DUSTY (Ivezić & Elitzur, 1997) to solve the radiative transfer equation for the single clouds, that where modeled as 1-D dust slab. The final torus SED was obtained by adding the emission from different slabs at different viewing (phase) angles, after statistically weighting them. They find that 5 to 15 clumps in a equatorial line of sight, each with an optical depth in the range , are successful in reproducing the observed characteristics of AGNs. Models in which the clouds are more concentrated at shorter distances from the central source, i.e. with a radial distribution following a power low, are favoured.
Another approach was followed by Dullemond & van Bemmel (2005), who exploit a 2-D Monte-Carlo code, in which the clumps where modeled as concentrical rings. Once the temperature of the dust is known throughout all the cells, the torus SED is calculated by means of ray tracing techniques. Starting from models with dust continuously distributed, they calculate the respective clumpy models, finding that it is not possible to use observed infrared data to distinguish between the effects due to the two different distributions.
The model developed by Hönig et al. (2006) and its further development (Hönig & Kishimoto, 2010), also adopts a Monte-Carlo technique to solve the radiative transfer problem, calculating the SEDs for various phase angles, for each cloud, setting up in this way a database of clumps emission. They consider that both the clouds optical depth and their size (radius) are related to their distance from the central source. The clouds are then randomly displaced, according to a spatial distribution function, and the torus SED is calculated by summing the emission of directly and non-directly illuminated clumps. This approach allows them to also study the dependence of the dust SED on the arrangement of the clouds: relevant differences are found especially for intermediate angle lines of sight.
Monte-Carlo, coupled to ray tracing techniques, are used by Schartmann et al. (2008), who are not using any prescription for the dust distribution which is instead computed from the equilibrium between the gravitational potential and pressure forces. They explore the effect of the clouds filling factor, of changing the dust mass, of the clump size and their positions. Again, their analysis of the SED for different arrangements of the clumps, shows non-negligible differences which tend to be the highest for edge-on views. The case of a non–isotropically emitting central source, whose emission is varying according to a law, was also studied by Schartmann et al. (2005), but their results are not directly comparable to ours since in their case the dust was continuously distributed.
In this paper we have investigated the infrared emission of AGN dusty tori. Following theoretical predictions coming from hydrodynamical simulations, we modeled the dusty torus as a 3D two-phase medium with high-density clumps and low-density medium filling the space between the clumps. We employed a 3D radiative transfer code based on the Monte Carlo technique to calculate SEDs and images of torus at different wavelengths. We calculated a grid of models for different parameters and analyzed the properties of the resulting SEDs. For each two-phase model we have calculated two corresponding models with the same global physical parameters: a clumps-only model and a model with a smooth dust distribution. For both two-phase and clumps-only models, another set is generated keeping all the parameters constant but varying the random distribution of the clumps. From the analysis of the SED properties and comparison of the corresponding models, we conclude the following:
The SED at near- and mid-infrared wavelengths is determined by the conditions of dust the innermost region of the torus: different random distributions of the clumps may result in the very different SEDs in otherwise identical models.
The shape of the silicate feature is not only a function of inclination. Optical depth, dust distribution parameters, clump size and actual arrangement of the clumps, all have an impact on the appearance of the silicate feature. Low optical depth tori produce silicate feature in a strong emission. Models with high-density clumps occupying the innermost region will have the emission feature attenuated due to the shadowing effects.
The clump size has a major impact on the SED properties. SEDs of the clumpy models with small clumps ( or clump size of pc) are very similar to the ones obtained by a homogeneous distribution of the dust. The silicate feature in absorption in these models is shallower and they tend to have less near-infrared emission than the corresponding smooth models. However, the silicate feature in emission is not suppressed. Clumpy models with bigger clumps ( or clump size of pc) are showing more differences compared to both small clump and smooth models. The silicate feature in absorption in these models is even less deep and they have less near-infrared emission than the small clump and smooth models. The silicate feature in emission is in general less pronounced. We stress that suppression strongly depends on the dust distribution parameters. The effect is the most notable in the case of a non-constant density in the radial direction and constant density in the polar direction (, ); as is allowed to increase the effect is weaker or even absent.
Although the silicate emission feature could be suppressed in the clumpy models for certain parameters, the smooth models are able to reproduce almost the same range of the silicate feature strength. Our analysis shows that, overall, when considering characteristics of the silicate feature, models with the three dust configurations (smooth, two-phase, clumps-only) are not distinguishable.
Low density dust, smoothly distributed between the clumps in the two-phase model, significantly contributes to the near-infrared emission in type 1 view. This is the main difference with respect to the clumps-only models that typically show a deficiency in this range. This peculiar characteristic of the two-phase models might represent a possible solution to a similar issue found when fitting observed SED with currently available clumpy models from the literature.
This work will be extended and the parameter grid will be progressively improved; models with different chemical composition of the dust, different torus and clumps sizes and their spatial distributions will be further explored. SEDs, in the form of ascii files are available on the following address: https://sites.google.com/site/skirtorus/. Images, in the form of fits files are available upon request.
We thank the anonymous referee for useful comments and suggestions. This work was supported by the European Commission (Erasmus Mundus Action 2 partnership between the European Union and the Western Balkans, http://www.basileus.ugent.be) and by the Ministry of Education and Science of Serbia through the projects ‘Astrophysical Spectroscopy of Extragalactic Objects’ (146001) and ‘Gravitation and the Large Scale Structure of the Universe’ (146003). The European Commission and EACEA are not responsible for any use made of the information in this publication.
- pagerange: Three-dimensional radiative transfer modeling of AGN dusty tori as a clumpy two-phase medium–Three-dimensional radiative transfer modeling of AGN dusty tori as a clumpy two-phase medium
- pubyear: 2011
- Alonso-Herrero A., Quillen A. C., Rieke G. H., Ivanov V. D., Efstathiou A., 2003, AJ, 126, 81
- Antonucci R., 1993, ARA&A, 31, 473
- Antonucci R. R. J., Miller J. S., 1985, ApJ, 297, 621
- Aretxaga I., Joguet B., Kunth D., Melnick J., Terlevich R. J., 1999, ApJL, 519, L123
- Baes M., 2008, MNRAS, 391, 617
- Baes M., Dejonghe H., 2001, ApJ, 563, L19
- Baes M., Dejonghe H., 2002, MNRAS, 335, 441
- Baes M., et al., 2003, MNRAS, 343, 1081
- Baes M., Dejonghe H., Davies J. I., 2005a, AIPC, 761, 27
- Baes M., Stamatellos D., Davies J. I., Whitworth A. P., Sabatini S., Roberts S., Linder S. M., Evans R., 2005b, NewA, 10, 523
- Baes M., et al., 2010, A&A, 518, L39
- Baes M., Verstappen J., De Looze I., Fritz J., Saftly, W., Vidal Pérez E., Stalevski M., Valcke S., 2011, ApJS, accepted, arXiv:1108.5056v1 [astro-ph.CO]]
- Barvainis R., 1987, ApJ, 320, 537
- Beckert T., Duschl W. J., 2004, A&A, 426, 445
- Bianchi S., 2008, A&A, 490, 461
- Bjorkman J. E., Wood K., 2001, ApJ, 554, 615
- Cashwell, E. D., & Everett, C. J., 1959, A Practical Manual on the Monte Carlo Method for Random Walk Problems, Pergamon, New York
- Davis S. W., Laor A., 2011, ApJ, 728, 98
- De Looze I., et al., 2010, A&A, 518, L54
- Deo, R. P., Richards, G. T., Nikutta, R., Elitzur, M., Gallagher, S. C., Ivezić, Ž., & Hines, D. 2011, ApJ, 729, 108
- Dullemond C. P., van Bemmel I. M., 2005, A&A, 436, 47
- Efstathiou A., Rowan-Robinson M., 1995, MNRAS, 273, 649
- Feltre A., Hatziminaoglou E., Fritz J., Franceschini A., in preparation
- Fritz J., Franceschini A., Hatziminaoglou E., 2006, MNRAS, 366, 767
- Gordon K. D., Misselt K. A., Witt A. N., Clayton G. C., 2001, ApJ, 551, 269
- Granato G. L., Danese L., 1994, MNRAS, 268, 235
- Hao L., Spoon H. W. W., Sloan G. C., Marshall J. A., Armus L., Tielens A. G. G. M., Sargent B., van Bemmel I. M., Charmandaris V., Weedman D. W., Houck J. R., 2005, ApJL, 625, L75
- Hony, S., et al. 2011, A&A, 531, A137
- Hönig S. F., Beckert T., Ohnaka K., Weigelt G., 2006, A&A, 452, 459
- Hönig S. F., Kishimoto M., 2010, A&A, 523, A27+
- Ivezić, Z., & Elitzur, M. 1997, MNRAS, 287, 799
- Jaffe W., Meisenheimer K., Röttgering H. J. A., Leinert C., Richichi A., Chesneau O., Fraix-Burnet D., Glazenborg-Kluttig A., Granato G., Graser U., Heijligers B., Köhler R., Malbet F., et al. M., 2004, Nature, 429, 47
- Kawaguchi, T., & Mori, M. 2010, ApJL, 724, L183
- Kawaguchi, T., & Mori, M. 2011, arXiv:1107.0678
- Kishimoto M., Hönig S. F., Antonucci R., Barvainis R., Kotani T., Tristram K. R. W., Weigelt G., Levin K., 2011, A&A, 527, A121+
- Krolik J. H., 2007, ApJ, 661, 52
- Krolik J. H., Begelman M. C., 1988, ApJ, 329, 702
- Laor A., Draine B. T., 1993, ApJ, 402, 441
- Li A., Draine B. T., 2001, ApJ, 554, 778
- Lucy L. B., 1999, A&A, 344, 282
- Mathis J. S., Rumpl W., Nordsieck K. H., 1977, ApJ, 217, 425
- Mor, R., Netzer, H.,& Elitzur, M. 2009, ApJ, 705, 298
- Nenkova M., Ivezić Ž., Elitzur M., 2002, ApJL, 570, L9
- Nenkova M., Sirocky M. M., Ivezić Ž., Elitzur M., 2008, ApJ, 685, 147
- Nenkova M., Sirocky M. M., Nikutta R., Ivezić Ž., Elitzur M., 2008, ApJ, 685, 160
- Netzer, H. 1987, MNRAS, 225, 55
- Niccolini G., Woitke P., Lopez B., 2003, A&A, 399, 703
- Nikutta R., Elitzur M., Lacy M., 2009, ApJ, 707, 1550
- Packham C., Young S., Hough J. H., Axon D. J., Bailey J. A., 1997, MNRAS, 288, 375
- Pier E. A., Antonucci R., Hurt T., Kriss G., Krolik J., 1994, ApJ, 428, 124
- Pier E. A., Krolik J. H., 1992, ApJ, 401, 99
- Pier E. A., Krolik J. H., 1993, ApJ, 418, 673
- Pogge R. W., 1988, ApJ, 328, 519
- Pogge R. W., 1989, ApJ, 345, 730
- Polletta, M., Weedman, D., Hönig, S., Lonsdale, C. J., Smith, H. E., & Houck, J. 2008, ApJ, 675, 960
- Schartmann M., Meisenheimer K., Camenzind M., Wolf S., Henning T., 2005, A&A, 437, 861
- Schartmann M., Meisenheimer K., Camenzind M., Wolf S., Tristram K. R. W., Henning T., 2008, A&A, 482, 67
- Schartmann M., Meisenheimer K., Klahr H., Camenzind M., Wolf S., Henning T., 2009, MNRAS, 393, 759
- Shi Y., Rieke G. H., Hines D. C., Gorjian V., Werner M. W., Cleary K., Low F. J., Smith P. S., Bouwman J., 2006, ApJ, 653, 127
- Siebenmorgen R., Haas M., Krügel E., Schulz B., 2005, A&A, 436, L5
- Tadhunter C., Tsvetanov Z., 1989, Nature, 341, 422
- Tristram K. R. W., Meisenheimer K., Jaffe W., Schartmann M., Rix H., Leinert C., Morel S., Wittkowski M., Röttgering H., Perrin G., Lopez B., Raban D., Cotton W. D., Graser U., Paresce F., Henning T., 2007, A&A, 474, 837
- Urry C. M., Padovani P., 1995, PASP, 107, 803
- van Bemmel I. M., Dullemond C. P., 2003, A&A, 404, 1
- Vidal E., Baes M., 2007, BaltA, 16, 101
- Vignali, C., et al. 2011, MNRAS, 1088
- Wada K., Norman C. A., 2002, ApJL, 566, L21
- Wada K., Papadopoulos P. P., Spaans M., 2009, ApJ, 702, 63
- Weigelt G., Wittkowski M., Balega Y. Y., Beckert T., Duschl W. J., Hofmann K., Men’shchikov A. B., Schertl D., 2004, A&A, 425, 77
- Weingartner J. C., Draine B. T., 2001, ApJ, 548, 296
- Witt A. N., 1977, ApJS, 35, 1
- Witt A. N., Gordon K. D., 1996, ApJ, 463, 681
- Wittkowski M., Balega Y., Beckert T., Duschl W. J., Hofmann K., Weigelt G., 1998, A&A, 329, L45
- Wolf S., 2003, CoPhC, 150, 99
- Yusef-Zadeh F., Morris M., White R. L., 1984, ApJ, 278, 186