Stellar feedback in AGN tori

The effect of stellar feedback on the formation and evolution of gas and dust tori in AGN


Recently, the existence of geometrically thick dust structures in Active Galactic Nuclei (AGN) has been directly proven with the help of interferometric methods in the mid-infrared. The observations are consistent with a two-component model made up of a geometrically thin and warm central disk, surrounded by a colder, fluffy torus component. Within the framework of an exploratory study, we investigate one possible physical mechanism, which could produce such a structure, namely the effect of stellar feedback from a young nuclear star cluster on the interstellar medium in centres of AGN. The model is realised by numerical simulations with the help of the hydrodynamics code TRAMP. We follow the evolution of the interstellar medium by taking discrete mass loss and energy ejection due to stellar processes, as well as optically thin radiative cooling into account. In a post-processing step, we calculate observable quantities like spectral energy distributions and surface brightness distributions with the help of the radiative transfer code MC3D. The interplay between injection of mass, supernova explosions and radiative cooling leads to a two-component structure made up of a cold geometrically thin, but optically thick and very turbulent disk residing in the vicinity of the angular momentum barrier, surrounded by a filamentary structure. The latter consists of cold long radial filaments flowing towards the disk and a hot tenuous medium in between, which shows both inwards and outwards directed motions. With the help of this modelling, we are able to reproduce the range of observed neutral hydrogen column densities of a sample of Seyfert galaxies as well as the relation between them and the strength of the silicate 10 m spectral feature. Despite being quite crude, our mean Seyfert galaxy model is even able to describe the SEDs of two intermediate type Seyfert galaxies observed with the Spitzer Space Telescope.

Galaxies: nuclei – Galaxies: Seyfert – ISM: dust, extinction – Radiative transfer – Hydrodynamics – ISM: evolution.

1 Introduction

Within the so-called Unified Scheme of Active Galactic Nuclei (Antonucci, 1993; Urry & Padovani, 1995), a torus-like geometrically thick gas and dust component is coined to explain two observed classes of Active Galactic Nuclei (AGN) by one intrinsically unique AGN model. This torus is illuminated by the accretion disk, surrounding the central black hole. As the dust opacity peaks in the UV/optical wavelength range, where also the accretion disk emits maximally, most of its light is absorbed and reemitted in the form of a pronounced peak in the mid-infrared (Sanders et al., 1989), which is frequently observed in AGN. This gives rise to the separation into type 2 sources, where the torus is viewed edge-on and type 1 sources, where a face-on view allows the visibility of the accretion disk, which shows up in form of the so-called Big Blue Bump in the spectral energy distribution (SED) of AGN. This geometrical unification idea was first proposed by Antonucci & Miller (1985), after broad emission lines – arising from fast moving gas close to the central black hole (the so-called Broad Line Region) – have been detected in polarised light in the Seyfert 2 galaxy NGC 1068, which was explained by scattering by dust and electrons above the torus opening. In direct light, type 2 sources only exhibit narrow emission lines, originating from slower gas motions further out, where they are less affected by the gravitational potential. This is the so-called Narrow Line Region, located beyond the opening of the torus funnel. The first direct evidence for a dusty torus came from interferometric observations with the help of the MID-infrared Interferometer (MIDI), which revealed geometrically thick dust distributions on parsec scale in the two Seyfert galaxies NGC 1068 (Jaffe et al., 2004; Poncelet et al., 2006; Raban et al., 2008) and the Circinus galaxy (Tristram et al., 2007).

Up to now, no conclusive, physical model for the distribution of gas and dust within these tori exists. This is due to the fact that very complicated and yet not fully understood astrophysical processes are thought to happen in the centres of AGN. Additionally, their large distances do not allow direct imaging observations. The major problem for the persistence of geometrically thick gas and dust distributions is to obtain stability of the vertical scale height against gravity. So far, several models have been put forward: Krolik & Begelman (1988) proposed that the torus is made up of clumps which possess supersonic random velocities, maintained by transferring orbital shear energy with the help of sufficiently elastic collisions between the clumps (see also Beckert & Duschl, 2004). To reach this elasticity, high magnetic field strengths are necessary. As these tori are made up of a multiphase mixture of gas and dust with cold, warm and hot components, a clumpy structure also helps to prevent the dust from being destroyed by hot surrounding gas. Further evidence for clumpiness or a filamentary structure of the cold component – in this case mainly for the distribution of neutral gas – comes from X-ray measurements of the absorbing column density distribution (Risaliti et al., 2002) and their combination with measurements of the strength of a spectral feature of a dust ingredient (see also discussion in Sect. 6.3). Very recently, Tristram et al. (2007) found hints for clumpiness in the dust torus of the Circinus galaxy by means of interferometric observations in the mid-infrared. Wada & Norman (2002) performed hydrodynamical simulations of the interstellar medium in centres of active galaxies under the assumption of starburst conditions. A very high supernova type II rate in a narrow sheet around the midplane (following in situ star formation in the densest region) is able to puff up an initially rotationally supported thin disk. In an alternative scenario, no torus is needed, but the necessary obscuration is given by dusty clouds, which are embedded into a hydromagnetic disk wind (Königl & Kartje, 1994 and discussion in Elitzur, 2006). The most recent approach comes from Krolik (2007), following an idea of Pier & Krolik (1992). His idealised analytical calculations show that the scale-height of AGN tori can be stabilised against gravity with the help of infrared radiation pressure.

A second approach towards modelling gas and dust structures of AGN is via radiative transfer calculations of specific dust distributions, which are motivated by simplified astrophysical scenarios. The most up-to-date simulations (Nenkova et al., 2002; Hönig et al., 2006; Schartmann et al., 2008; Nenkova et al., 2008a, b) feature clumpy dust distributions. From these calculations, we get some rough idea of the geometrical distribution of dust within AGN from comparison with high resolution spectral (e.g. NACO or Spitzer) as well as interferometric observations with MIDI (the MID-infrared Interferometer).

The goal of this work is to explore the effects of stellar feedback from a nuclear stellar cluster on the formation and evolution of AGN tori in terms of an exploratory study. In this scenario, the torus is mainly built up by gas loss due to stellar evolution. Energetic supernova explosions provide an effective structuring mechanism. To this purpose, we deploy three-dimensional hydrodynamic simulations. Eventually, this model should account for the obscuration as well as the feeding of the central source. In a second step, we convert the gas density distribution into a dust distribution, which we then use for continuum radiative transfer calculations, which can later be compared to observational results. In this work, we focus on low-luminosity AGN (Seyfert galaxies), as they are more abundant in the local universe and, therefore, their nuclear regions are accessible via current generation interferometric instruments.

After having discussed the main physical assumptions of our work in Sect. 2, we will briefly present the numerical method we use (Sect. 3). In Sect. 4, we describe the results of this exploratory study in terms of the evolution of our standard torus model. Furthermore, we characterise its final state and explain the results from a study of varying several model parameters. After a critical discussion (Sect. 5), we compare our torus simulations to recent observations of Seyfert galaxies (Sect. 6) in terms of the spectral energy distributions of dust re-emission, as well as gas extinction column densities and summarise our work (Sect. 7).

2 Physical assumptions of our model

The main assumption of our model is that a young star cluster exists in those nuclei of AGN, which possess a torus. Evidence for this was given recently by observations of a number of Seyfert galaxies observed with adaptive optics techniques (Davies et al., 2007) and for the special case of NGC 1068 with the help of HST/NICMOS images (Gallimore & Matthews, 2003 and references therein), one of the closest – and therefore best studied – Seyfert 2 galaxies. We model this nuclear stellar cluster in form of a Plummer profile (Plummer, 1911) of the gravitational potential and assume that it was built up during a short-duration starburst3, in concordance with the findings of Davies et al. (2007), which forms a so-called Coeval Stellar Population (CSP). The first phase of evolution of the stellar population is expected to be very violent, such that supernova type II explosions might evacuate most of the nuclear region from gas and dust. During this phase, we do not expect that a stable gas and dust distribution (torus) exists. According to our population synthesis modelling with the starburst 99-code (Leitherer et al., 1999; Vázquez & Leitherer, 2005), this violent phase lasts for approximately 40 million years. Within about the same time period, double systems of lower mass stars have been able to form. According to models which include binary systems (de Donder & Vanbeveren, 2003), the rate of supernova type Ia explosions – arising from so-called double degenerate systems – begins to rise steeply and starts to dominate the input of energy. At the same time, mass injection should be dominated by the emission of planetary nebulae from stars in the intermediate mass range (from 1.5 to 8, Kwok, 2005).

2.1 Mass input

Stars of all masses lose material during their lifetime in form of winds of different strengths, mainly happening within short periods (e.g. at the tip of the Red Giant Branch). For the case of high-mass stars, also supernova type II explosions have to be taken into account. The time dependent mass loss rate of a CSP can in principle be calculated with the following integration over the initial stellar mass for a given (time dependent) mass loss rate of individual stars and initial mass function (IMF) :


As both functions are not very well known, a number of simplifying assumptions have to be applied. We use the model of Jungwiert et al. (2001), which summarises the various phases of mass loss in one delta peak for the intermediate mass stars (IMS) and the high mass stars (HMS) and two delta peaks for low mass stars (LMS). For the IMS, the delta peak represents the AGB (Asymptotic Giant Branch) wind peak and for the HMS the SNe explosions plus winds. For the LMS, it represents the wind loss at the tips of the RGB (Red Giant Branch) and AGB. The IMF is modelled according to Scalo (1998) as a broken power law. Using the stellar initial-final mass relation and the stellar lifetime-mass relation from the Padua stellar evolutionary models (Bressan et al., 1993; Marigo et al., 1996), an approximate, normalised mass loss rate function can be derived:


where the subscript n denotes that it is normalised to the initial mass of the CSP. This leads to a mass loss rate of (normalised to a total stellar mass of ) at the beginning of our simulations (after 40 Myr). Recent observations tend to find rather top heavy IMFs in the centres of galaxies. This means that the number of massive stars is increased compared to low-mass stars. It has been shown for the centre of the Milky Way (Stolte et al., 2002, 2005; Nayakshin et al., 2006; Paumard et al., 2006), the Andromeda galaxy (Bender et al., 2005) as well as more distant galaxies (Scalo, 1990; Elmegreen, 2005). Therefore, we use a five times larger value of in our standard model, which is kept constant over the whole evolution of the simulation. For simplification reasons, we assume that all of the mass is injected in form of planetary nebulae.

Numerically, mass injection means a source term in the continuity equation as well as the momentum equation, when taking the ejection velocity of the gas into account. According to the mass injection rate, the number of planetary nebulae (PNe) to inject4 into the whole model space is calculated. Residual mass is injected into the model space within the subsequent time step. The position of the ejecting star is chosen randomly following the distribution of stellar mass. The ejection velocity of the PN is chosen such that it matches the velocity dispersion (assumed to be independent of position) plus the rotational velocity of the stars in the central star cluster. As we cannot follow the evolution of single stars explicitly, all of the ejected planetary nebulae have the same gas mass , which corresponds to the mean mass by weighting the initial-final mass relation (given by Weidemann, 2000) with the Initial Mass Function (IMF), which in this mass regime is well approximated by a Salpeter-like function (Salpeter, 1955). This leads to an average mass of per planetary nebula, which we distribute homogeneously over cells surrounding the randomly chosen position.

A detailed modelling of the temperature structure of a planetary nebula is beyond the reach of our current resolution, as we are not able to resolve the small-scale physical processes. As it will finally thermally merge with the surrounding gas on short time-scales, we choose the temperature to be the same as the gas in the cells where the PNe are injected.

2.2 Energy input

Our simulations start after roughly 40 Myr, when the violent phase of supernova II explosions (e.g. a mass star lives only about 5 Myr) has ended. At roughly the same time, the rate of SN Ia explosions begins to dominate the energy input into the ISM5 – a result found in population synthesis modelling including binary systems (de Donder & Vanbeveren, 2003). According to today’s theoretical understanding, SN Ia are thought to be thermonuclear disruptions of accreting white dwarfs, which reach the Chandrasekhar mass (), being the maximum mass which can be supported by electron degeneracy pressure. The determination of SN Ia rates for the considered nuclear stellar system in our simulations is complicated. Observations have to rely on low number statistics. Additionally, only the stellar content of the whole galaxy can be taken into account. A recent publication by Sullivan et al. (2006) describes the SN Ia rate for young and old stellar populations and finds that it can be well represented by the sum of (old population of stars) and (young population of stars). The rate for the old stellar population is normalised to 1  of stars and the latter rate to a star formation rate of 1 .

For an assumed duration of the starburst of , this results in a normalised SN Ia rate of roughly for our modelling. Due to the large uncertainties in theoretical derivation of supernova rates and as it is doubtful whether the observational derivations are appropriate for our simulations, we consider the SN Ia rate a more or less free parameter within our models. In an alternative starburst scenario for nuclei of Seyfert galaxies, Wada & Norman (2002) use a much larger rate of roughly 100 times our SN-rate, namely 1 SN detonation per year, restricted to a narrow strap around the midplane of their model space. This is motivated by violent star formation within a dense gas disk in the galactic midplane, accounting for SN II explosions. After a SN-rate parameter study, a value of was chosen for our standard model. It is assumed to be constant in our simulations, as our current period of time evolution (1.2 Myrs) is short compared to the expected evolutionary time scales of AGN tori, of the order of 100 Myrs. From the supernova rate for the whole model space, the number of supernovae exploding in each timestep can be determined. When a supernova is launched, the position is determined randomly, but according to the underlying stellar distribution. The energy is injected into a single cell in the form of thermal energy. Additionally, a mass of is added to this cell, also ensuring numerical stability. We use an efficiency of , reducing the input of thermal energy to erg. This factor accounts for energy losses in the early phases of evolution of the supernovae, which we cannot model at the current resolution. The problem of radiating away most of the introduced thermal energy due to vigorous cooling in very dense regions before a supernova shell can form is mitigated by introducing a cooling delay of several timesteps. Direct modelling of an expanding shell as initial condition for a supernova explosion is not feasible at the resolution of our current simulations. These shock fronts efficiently dissipate their kinetic energy when interacting with dense blobs or filaments or dense shock fronts of neighbouring supernova explosions, transferring their kinetic energy back to thermal energy. Similar procedures of energy input have been applied to many different astrophysical scenarios. For example Wada & Norman (2002) also introduce the complete amount of supernova energy into a single cell in form of thermal energy within a starbursting gaseous disk with even higher ambient densities compared to our simulations. Purely thermal energy is also introduced by Joung & Mac Low (2006) for the case of the interstellar medium in our own galaxy, where the mean gas density is much lower compared to our computations. Given the aims of this exploratory study, we therefore think that our energy injection mechanism is an appropriate treatment of this feedback process.

2.3 Optically thin cooling

Figure 1: Cooling curve (Plewa, 1995). The hatched regions denote three temperature regimes: cold neutral medium (CNM), warm neutral medium (WNM) and hot ionised medium (HIM). The thick red line segments mark the thermally stable parts of the cooling curve. Dotted lines refer to cooling curves used by other groups, see text.

A mixture of neutral and ionised gas mainly cools via the transformation of kinetic energy to radiation by collisional processes. This also means that the efficiency of cooling is a sensitive function of the constituents of the gas. The radiation can only escape, if the medium is optically thin for these wavelengths, otherwise cooling is inefficient.

Fig. 1 shows the effective radiative cooling curve applied in our simulations (Plewa, 1995), which was calculated with the help of the Cloudy code (Ferland, 1993). For the case of our cooling curve, a purely collisionally ionised optically thin plasma with solar abundances was used.

Following Joung & Mac Low (2006), such a gas is thermally stable at temperatures, where the slope in logarithmic scale is steeper or equal to one. For the case of our cooling function, we find five distinct stable regions: K, K, K, K and K. They are marked with the thick red line segments in Fig. 1. Overplotted as yellow dotted lines are the often used cooling functions from Dalgarno & McCray (1972) – where only the part below K is shown – with ionisation fractions of 0.1 for the upper curve and 0.01 for the lower one. The violet dotted line results from equilibrium ionisation cooling, calculated by Sutherland & Dopita (1993) for the case of solar metallicity. As can be seen from this, the cooling curves are in relatively good agreement in the high temperature regime (above the ionisation temperature of hydrogen at around K). At temperatures lower than K, the values of the assumed cooling curves are most uncertain, as in this temperature regime, the approximation of thermal equilibrium is usually not valid. Therefore, large deviations exist between the different calculational methods. It is also a matter of debate, which ionisation fraction to assume for the low temperature part. We decided to use 0.1 (black curve and upper yellow curve in Fig. 1) in concordance with de Avillez & Breitschwerdt (2004), whereas Joung & Mac Low (2006) assume an ionisation fraction of 0.01 (lower yellow curve). For temperatures lower than 10 K, we set the cooling rate to zero.

Similar as in Joung & Mac Low (2006), we make a division into three temperature regimes: the Cold Neutral Medium (CNM) up to 175 K, where the first stable part of the cooling curve ends, the Warm Neutral Medium (WNM) in the intermediate temperature regime up to the ionisation temperature of hydrogen and the Hot Ionised Medium (HIM) for higher temperatures (compare to Fig. 1).

2.4 Parameters of our standard model

The parameters we use for our standard model are summarised in Table 1. In the past subsections, we already motivated the supernova as well as the mass injection rate. For the remaining parameters, short explanations are given in this paragraph.

The underlying gravitational potential is made up of two components:

  • The potential of the nuclear star cluster, which is modelled according to a Plummer distribution:

  • The Newtonian point-like potential of the nuclear black hole:

Parameter Value Parameter Value
25 pc
5 pc SNR
50 pc
165 km/s

Mass of the black hole (), normalisation constant of the stellar potential (), initial gas mass (), cluster core radius (), torus radius (), outer radius (), stellar velocity dispersion (), exponent of the angular momentum distribution (), initial gas temperature (), normalised mass injection rate (), mass of a single injection (), supernova rate (SNR) and adiabatic exponent ().

Table 1: Parameters of our standard model.

A core radius =25 pc, a black hole mass and a normalisation factor for the stellar distribution of is used. The latter leads to a total integrated stellar mass within the core radius of . Our initial condition comprises of a TTM-model similar to the one discussed in Schartmann et al. (2005), but with a Plummer potential (as given above) and the parameters listed in Table 1. The angular momentum of the stellar content, the gas in the initial condition as well as the injected gas is given by the following distribution:


with a constant and a torus radius of , where the torus radius is defined as the minimum of the effective potential. The turbulent pressure of the TTM-model is substituted by thermal pressure. For the case of our standard model, the velocity dispersion of the assumed dust clouds in the TTM-model of km/s then corresponds to an initial temperature of = K. The mass loss rate was chosen to be , which means a mass injection of solar masses of gas per year, normalised to one solar mass of stars (see discussion in Sect. 2.1). A SN-rate of supernovae per year and normalised to in stellar mass was used (see Sect. 2.2 for further discussion). We apply the equation of ideal gas as a closing condition with an adiabatic index of an ideal monoatomic gas (). The gas has solar abundances in our simulations, with a mean molecular weight, which changes from 0.6 (fully ionised gas with solar abundances) to 1.2, following a smooth transition at the ionisation temperature of hydrogen. Most of the simulations discussed in this paper ran for an evolutionary period of 10 orbits at the torus radius (global orbits6), corresponding to yr. Note that this is a tiny fraction of the age of the stellar population. Therefore, evolutionary effects can be ignored.

3 Numerical method

For our simulations, we use the 3D radiation hydrodynamics code TRAMP (Kley, 1989; Klahr, 1998; Klahr et al., 1999), which uses similar numerical schemes as the Zeus code (Stone & Norman, 1992a, b; Stone et al., 1992). It solves the equations of ideal hydrodynamics with the help of an Operator Splitting technique on a spherical grid that is logarithmically spaced in the radial direction. Using this method, one has always an optimal resolution and at the same time more or less boxy e. g. equi-spaced grid cells. For the advection step, the monotonic transport scheme of van Leer (1977) is applied, which is a spatially second-order-accurate upwind method. The radiation transport in TRAMP is treated in the flux-limited diffusion approximation using the flux limiters of Levermore & Pomraning (1981). Additional source terms for mass input (Sect. 2.1) and energy input (Sect. 2.2) are treated in a separate step. TRAMP discretises the hydrodynamic equations in the finite volume description and, therefore, locally and globally conserves advected quantities. The force step is calculated using a finite difference discretisation. Shocks are treated by introducing a non-linear artificial viscosity, which smears out discontinuities over several cells and, thereby, results in the correct generation of entropy in the shock, the correct shock velocities and suppresses post-shock oscillations. It is treated as an additional contribution to the pressure. The non-linearity ensures that the artificial viscosity is large in shocks, but negligible elsewhere.

3.1 Domain decomposition and boundary conditions

In order to gain resolution, we split the whole computational domain into three parts: the inner one ranges from to pc, the middle one from to pc and the outer part from to pc (logarithmic grid). Towards the centre, the size of the single cells decreases and, additionally, the obtained velocities increase due to the presence of the central gravitational potential and the conservation of angular momentum of infalling material. This leads to ever smaller timesteps, in order to prevent information to be transported along more than one cell during a single timestep (Courant-Friedrichs-Levy condition, Courant & Friedrichs, 1948; Ritchmyer & Morton, 1967). Hence, we are able to calculate a much longer time period for the outer domain with the same amount of CPU usage and the same grid size, compared to the inner domain. On the other hand, further in, processes act on much shorter time-scales and only a shorter time evolution is needed in order to obtain a steady state. Therefore, we calculate 10 global orbits for the outer domain, but only the last orbit for the middle domain and one tenth of an orbit for the innermost part. This leads to a comparable amount of CPU time of approximately one weak per domain.

In the first step, the outer domain is calculated. The physical state of material at the position of the outer boundary of the middle model is stored in a file every 0.1% of a global orbit. Outflow boundary conditions are used radially inwards as well as outwards. This means we allow for the outflow of matter from the grid and minimize the reflection of waves at the grid edge, while not allowing for inflow, thus preventing the boundary cells from generating a numerical instability. During the calculation of the middle domain in the second step, the inner radial boundary is set to outflow, whereas the outer boundary is set by the afore stored quantities7, which can result in in- or outflowing material, thanks to the wide overlap of the domains. At the same time, physical quantities advected across the position of the outer boundary of the inner domain are stored. In the last step, the same is done for the innermost model space. For the final quantitative analysis and the subsequent radiative transfer calculations, the three domains are combined again. This procedure is a reasonable approximation, as in most of the simulations we observe infall of gas towards the inner region. Outflowing material can partly be taken into account, as the domains have a quite large overlapping region.

Concerning the other coordinates, we model a range of with periodic boundary conditions in direction. In -direction, where we model from 4 degrees to 176 degrees, as well as for the inner- and outermost boundary in radial direction, outflow boundary conditions are applied in a sense that mass is allowed to flow out of the integration domain, but no material can get back into it. This means that the values of the ghost cells 8 are set such, that the normal derivatives of the boundary values equal to zero for all quantities, presenting a smooth continuation of the flow through the boundary. Furthermore, we only allow outward directed motion, otherwise the normal velocity components are set to zero. For each of the domains, we are able to calculate 75 cells in radial, 94 cells in and 49 cells in direction, leading to cell sizes between 0.006 pc in the innermost part of the whole domain and 1.5 pc at the outer edge.

3.2 Initial condition

A hydrodynamically stable initial condition for the case of our effective potential – made up of the nuclear stellar cluster, the black hole and the additional centrifugal potential due to orbital rotation – is given by the isothermal Turbulent Torus Model (TTM, Camenzind, 1995, Schartmann et al., 2005). A conical volume around the rotation axis with a half opening angle of needs to be cut out, because the funnel region of this field configuration is dynamically unstable and would cause unphysical behaviour of the state vector.

4 Results

Figure 2: Temporal evolution of the density in our standard model, given in a meridional slice of the torus. Shown are four stages of the torus at 0.25, 1.0, 5.0 and roughly 10.0 global orbits. The scaling is logarithmic.
Figure 3: State of the torus after roughly ten global orbits. Shown are from left to right the density distribution, the temperature distribution and the pressure distribution of a meridional slice through the 3D data cube. All images are displayed in logarithmic scaling.

4.1 The evolution of our standard model

Within the first few percent of a global orbit, single injected blobs of gas are visible, mainly within the central 25 pc, the core radius of the stellar distribution (Fig. 2a). Spherical cavities are formed due to the explosion of supernovae. While these explosions predominantly form shock fronts, the planetary nebulae interact in two different ways, as we could see in separate simulations of single events (see Schartmann, 2007): In the subsonic case (e.g. hot surrounding medium), colliding streams tend to stick together and form long filaments in perpendicular direction from the initial direction. The supersonic interaction of two planetary nebulae (e.g. in cold medium) will cause the formation of a common shock front, which surrounds the point of interaction. With ongoing evolution, the single blobs collide to create a filamentary, net-like structure, enhanced and further shaped by supernova explosions (Fig. 2b and c). This process is also amplified by the onset of cooling instability in the dense regions of the filaments. Cooling also leads to loss of thermal pressure and, therefore, material tends to fall towards the minimum of the potential well, located at the torus radius (5 pc) within the equatorial plane of the torus. Around this radius, the filamentary structure opens out into a geometrically thin, but optically thick and very turbulent disk (Fig. 2d).

This process is comparable to the so-called galactic fountain process (Shapiro & Field, 1976) visible in galactic disks. The more tenuous the gas, the more effective are the supernova explosions in pushing away material. This mechanism finally leads to the large-scale, filamentary structure visible in Fig. 2d.

As the implemented optically thin cooling scales proportional to the gas density squared, only the dense filaments cool on short time-scales, whereas the interclump-medium is effectively heated by the energy injection of the supernova explosions. Accordingly, the gas velocity can appreciably differ in various temperature regimes, as will be discussed in Sect. 4.3.

Figure 4: Phase diagrams for the total model space of our standard model after an evolution of ten global orbits. The hatched regions indicate thermally stable regimes according to the cooling curve (see thick red line segments in Fig. 1).

In order to ensure that the memory of the initial condition is completely lost, we start calculating the middle part of the computational domain after an evolution of nine global orbits of the outer domain. Already after a fraction of the simulated one global orbit, a geometrically relatively thin, but optically very dense disk forms in the equatorial plane in the vicinity of the minimum of the potential (see Fig. 6 and discussion in Sect. 4.4).

The final state of our standard model after ten global orbits is depicted in terms of density, temperature and pressure within a meridional plane in Fig. 3. A multiphase medium has formed (Sect. 4.2). The temperature and the density distribution are complementary in the sense that high density regions possess low temperature due to rapid cooling and hot regions possess low density as material has been blown away by supernova explosions. It is important to note that the pressure distribution (right panel in Fig. 3) is far away from pressure equilibrium, as it spans more than five orders of magnitude.

4.2 Multiphase medium

Fig. 4 shows phase diagrams for the total model space of our standard model (all three spatial subregimes) after ten global orbits. The hatched regions indicate thermally stable parts of the cooling curve (see Sect. 2). Capital letters (A-F) denote the six bumps in the probability density function for gas temperature (upper left panel). C and D as well as the small excess E can be directly related to the stable regions of the cooling curve, because time spent in each temperature range is proportional to the inverse of the cooling function (Gerola et al., 1974). At these temperatures, gas cools comparatively slow and therefore assembles in these regimes. The explanation for the maximum B is identical. The minimum in-between A and B evolves, because gas at these temperatures is mostly located in dense clumps. As cooling depends on the square of the density, it can cool quite fast, even if the cooling curve indicates thermally stable behaviour. Peak F can be explained by constant energy input from supernova explosions, which allows gas to reside even in the thermally unstable regime (see also discussion in de Avillez & Breitschwerdt, 2005). Looking at the mass distribution of the various temperature phases (lower left panel), a similar curve is obtained, but tilted, as most of the mass is concentrated in the cold temperature regime (dense filaments and disk), while the hot gas fills most of the volume (supernova blown cavities).

Plotted with respect to the gas density (second column of Fig. 4), three local maxima are visible (G,H,I). Feature I is caused by the very dense disk component, as it is dominant in mass, but unimportant in volume. This material had enough time to cool to temperatures below 175 K and, therefore, corresponds to maximum A and B in the left panel. G can be accounted to regions, where supernova explosions happened recently, contributing to a large volume fraction, but only little to the total mass. In the left column, this corresponds to feature F. The bump (H) in between arises from material in the large scale filaments, which is in the process of settling down and thereby cooling.

The distribution of mass and volume with respect to pressure is shown in the right column of Fig. 4. As already evident from the two-dimensional cuts in Fig. 3, no pressure equilibrium is reached. This fact is in qualitative agreement with recent observations of Jenkins & Tripp (2007) in our local neighbourhood, who find variations in the pressure of the diffuse, cold neutral medium in a range of  KK  by using C i absorption lines, measured in the ultraviolet spectra of roughly 100 hot stars, taken from the HST archive. Instead of a single and sharp peak (expected for pressure equilibrium), we obtain three local maxima (right panel in Fig. 4). The system tries to reach an equilibrium value (located around peak K). Peak L can be assigned to regions with recent supernova action, leading to the highest pressure values of the simulation and filling a large volume. This is the same component causing the peaks F and G. We interpret deviations from a global equilibrium pressure value towards small pressure (peak J) as a sign of ongoing thermal instability within the fast cooling, densest regions of the flow (see also discussion in Joung & Mac Low, 2006). This interpretation is backed by Fig. 3, where cold and dense regions correspond to low pressure states and hot and tenuous regions to high pressure values.

In terms of volume and mass filling factors, the following picture evolves: Most important in mass is the CNM with 60.4% of the total mass, mostly residing in the disk component and the inner parts of the long filament, which had enough time to cool to temperatures of the CNM and occupies the smallest fraction in volume (only 5.4%). Second most important (22.6% in mass) is the WNM, made up of material in the process of cooling down. Only 16.9% of the mass is in the HIM, which occupies most of the space (83.3% of the volume). Only 11.3% of the whole model space is occupied by the WNM.

4.3 Dynamical state of the torus

Figure 5: Radial velocities of the outer model space between eight and ten global orbits, separated into temperature phases. Negative velocities refer to infalling motions. The blue line represents the free fall velocity (), scaled with a factor of , the dashed line symbolises the location of the core radius of the stellar distribution.
Figure 6: Formation of a turbulent, geometrically thin disk close to the midplane. Shown is the logarithm of the density in units of g/ as a cut through the midplane (upper row) and through the central meridional plane (lower row). Time is given in global orbits. Labels are given in pc.

The build-up of the disk and torus component can also be quantified by means of radial velocity distributions, as given in Fig. 5. The velocity is averaged over all angles and and separated into the three temperature components. For an integration time between eight and ten global orbits, the data are plotted every 0.05 global orbits (single dots). The solid lines are temporal averages for the whole range of displayed data. For comparison, we show the free fall velocity of a test particle within the total gravitational potential of the black hole and the stellar distribution, scaled down by a factor of three. It is a remarkable feature of our simulations that the slope of the velocity distributions of all three components are similar to free fall.

Clearly, the hot tenuous medium is able to provide the largest thermal pressure support, which slows down the infall towards the mass centre. In contrast to this, cold material sinks much faster to the bottom of the potential well, within the long, dense filaments.

At the border line between the phases with different flow structures, Kelvin-Helmholtz and Rayleigh-Taylor instabilities are expected to occur, but are not visible in our simulations due to the lack of sufficient resolution. They are supposed to increase turbulence and promote further mixing. Due to their high density, the filaments are very stable and single supernova explosions are not able to disrupt them. In them, the CNM reaches radial velocities of almost 30% of the free fall velocity, the WNM approximately 15% and the HIM 10%. This is a further indication that energy in form of supernovae or due to the planetary nebulae injection is mainly dispensed into the hot medium, which is also visible in the larger scatter of the radial velocities of the hot component.

4.4 Formation of a turbulent disk

As already discussed in Sect. 4.1, in the vicinity of the torus radius, a geometrically rather thin, but optically very dense disk forms. Fig. 6 shows the temporal evolution of density for four snap-shots of the middle domain. Its dynamical evolution is mainly determined by mass inflow through the outer radial boundary. This material reaches the middle part both in the form of broad streams as well as in blobs and filaments of material and gets structured by supernova explosions and planetary nebulae injection. As more and more material falls towards the minimum of the potential, a turbulent, fluffy disk emerges between the inner boundary and roughly six or seven parsec. We note the remarkably filamentary structure of the disk, showing strings of material, elongated in azimuthal direction (for reasons identical to those already discussed for the outer domain). The difference is that the average flow is here directed in azimuthal direction and not in radial direction.

Altogether, the evolution can be summarised as a collapse from a geometrically thick distribution (the torus component) to a geometrically thin disk. The initial torus is mainly (thermal) pressure supported and, therefore, collapses due to radiative cooling processes into dense filaments, as cooling and heating of the gas is locally unbalanced. This finally leads to the formation of the thin disk in the vicinity of the minimum of the gravitational potential, close to the angular momentum barrier. Thermal pressure does not play a role anymore, as the disk has cooled to temperatures of less than 100 K and reaches in some regions even our minimal temperature of 10 K. The scale height is now determined by random motions of the filamentary disk and rotation. This hydrodynamical turbulence leads to angular momentum transfer outwards and accretion towards the centre.

4.5 Parameter studies

In order to investigate the dependencies of our model, we conducted several parameter studies, which are briefly presented here (see Schartmann, 2007 for a more detailed discussion).

  • Mass and energy input rate: Changing the mass injection rate and the supernova rate, we found that a complementary behaviour is obtained. Both, with an increasing supernova rate or a decreasing mass-loss rate, the inward bound flow is more and more turned into an outflow, as demonstrated for example in Fig. 7. This has also been shown in 1D hydrodynamical simulations of the interstellar medium in elliptical galaxies by Gaibler et al. (2005). The transition radius between inflowing and outflowing hot gas increases with increasing mass loss rate or decreasing supernova rate. Inside the transition radius, the curves flatten for models with larger mass input rate. This is caused by the additional turbulent pressure component, provided by the PNe injection of gas blobs with random velocities.

    Figure 7: Comparison of radial velocities within the mass loss rate study for the HIM between global orbit eight and ten. Lines represent time-averaged (global orbits eight to ten) values.
  • Cooling rate: Cooling is a sensitive function of metallicity. The larger the metallicity, the larger the cooling rate. We varied the cooling rate to account for the unknown metallicity of the gas in AGN tori. The most striking difference appears in the broadening of the filaments, when increasing the cooling rate, as gas of even lower densities is able to effectively cool to CNM temperatures.

  • Mass of the central stellar cluster: Varying the (generally unknown) mass content in stars is of interest, as it first determines the gravitational potential in the outer part beyond the influence of the black hole. Second, it changes the amount of supernova input and mass input in equal measure. A higher stellar mass produces a deeper potential well and therefore accelerates the gas to higher velocities and the transition radius between inflow and outflow of the hot temperature phase increases. For the case of a higher total stellar mass, the evolving structures are more compact and a more massive disk is able to form at the angular momentum barrier.

  • Various initial conditions: In order to check the dependence of our results on the initial condition, we started one simulation with an initially empty model space (possessing a homogeneous density distribution with the lower density threshold of our simulations –  g/). As expected, this has no visible consequence on the evolution of the gaseous structures and a steady state in terms of total gas mass is already reached after less than 3 global orbits.

    Also starting with cold ( K) initial conditions, where the initial torus is supported by turbulent pressure on small scales yields identical results. This small-scale turbulence decays within a fraction of an orbit and is dissipated into heat. Additionally, the supernova input quickly heats the tenuous medium outside-in towards the minimum temperature at the torus radius. Therefore, after a short time, we again get a very similar behaviour as in our standard model.

5 Discussion

It needs to be emphasized that within our model, the large-scale torus component is not a long-living steady state, but it can only exist as long as mass and energy is injected from the stellar cluster. After switching off energy and mass input, most of the large scale structure will disappear within approximately 1.5 global orbits. The turbulent disk in the vicinity of the angular momentum barrier is the only remaining structure. In reality, this shutdown will happen naturally due to the evolution of the nuclear star cluster, but on much longer time-scales. It means that – within our modelling – the existence of tori is closely linked to the existence and evolution of a stellar cluster in the centre of a galaxy. This fits well to the observations by Davies et al. (2007), which reveal a link between nuclear starbursts and AGN activity with a certain time delay in-between. Furthermore, it needs to be stressed that the torus scale height in our modelling is mainly determined by the spatial distribution of the stars.

The total simulated time of our runs is of the order of 10 global orbits, which corresponds to  Myr. This is sufficiently long to avoid memory effects from the initial condition and at that point, a global dynamical equilibrium concerning the total mass and total energy is reached for the outer domain. In the middle and inner domain, the calculated time evolution of one and 0.1 global orbits is not long enough to reach a global dynamical equilibrium state. However, the technical limitation of this study did not allow a longer evolution.

5.1 Limitations of our simulations and future work

This paper presents an exploratory study of a self-consistent hydrodynamical model of gaseous and dusty AGN tori. Due to the simplifications required by the efficiency of our code, some important problems and astrophysical effects remain open:

  • Longer time evolution: With our serial code we are limited in temporal evolution to about 1.2 Myrs, while typical evolutionary times of such systems are of the order of 10 to 100 Myrs. Therefore, we currently port our routines to the parallel PLUTO code (Mignone et al., 2007).

  • Radiation contribution from central source and stars: The effect of the radiation of the central accretion disk is twofold: First of all, it contributes to the heating of the dust and gas. Second, its radiation pressure plays an important role in the dynamics of AGN dust tori, as was recently shown by Krolik (2007). With the help of our detailed radiative transfer calculations, we compared the accelerations due to radiation pressure of the dust reemission with the gravitational accelerations in vertical direction for each cell within the computational domain of the snapshot after roughly 10 global orbits, using a flux mean dust opacity in combination with a grey approximation. In this estimation, infrared radiation pressure significantly contributes only close to the midplane in the inner few parsecs from the heating source, which is in accordance with the detailed calculations of Krolik (2007). Within the large scale filaments, the ratio between the accelerations due to the infrared radiation pressure and the gravitational acceleration in vertical direction amounts to a few percent only, except for the direct illuminated parts of the filaments, where it can also exceed the gravitational accelerations and will therefore drive an outflow along the torus axis, widening the (so far too narrow) opening angle in our simulations. However, detailed radiative transfer calculations during the whole course of the hydrodynamical simulations are necessary to finally clarify this question. Opposed to Thompson et al. (2005), who study radiation pressure-supported starburst disks, we are currently interested in the phase following an efficient starburst. By then, the stellar luminosity has already dropped by a large factor. Giving the low stellar luminosity densities derived from the observations of a Seyfert galaxy sample discussed in Davies et al. (2007), the stellar luminosity seems to be only a small fraction of the luminosity of the central accretion disk at this evolutionary stage. Therefore, we currently neglect its contribution to the dynamics, as well as to the heating of the dust within our postprocessing radiative transfer calculations.

  • Magnetic fields: Weak magnetic fields in combination with shear due to differential rotation cause magneto-rotational instability (Balbus & Hawley, 1991), which provides a means to transport angular momentum and thus enables accretion towards the core (that is through the inner boundary of our computational domain). Probably coexisting strong magnetic fields in a spatially distinct region cause hydromagnetic outflows, which might be preferentially directed along the torus axis and thereby widening the opening angle of the torus.

  • Self gravity and star formation: At the current state of our modelling, self-gravity of the gas is not taken into account. However, examining the Jeans-criterion for each cell of our model space, we find that part of the inner turbulent disk, as well as a small fraction of the filamentary structure is Jeans unstable and will form stars. This is consistent with the findings of Nayakshin et al. (2006) and Paumard et al. (2006). In the snapshot after 10 orbits, the unstable mass amounts to 0.3% in volume or 15% in mass respectively. However, it should be emphasized again that our purely hydro- dynamical code is not capable to treat the accretion flow towards the core correctly, leading to a pile-up of cold material in the disk which does not describe the steady state density. Subsequently, supernova type II explosions might inflate or even partly disrupt the evolving dense and geometrically thin turbulent disk. However, one should note again that detailed simulations including self gravity will be subject of our future work.

  • Further mass input: There will also be material transported into the simulated domain of the active galactic nucleus from further out of the galaxy. Prominent examples are dusty mini-spirals, of which some have been discovered by NACO (Prieto et al., 2005) or HST-ACS observations (Hubble Space Telescope – Advanced Camera for Surveys). Simulations (Maciejewski, 2006) show that with the help of such nuclear spirals, enough gas inflow can be produced in order to power luminous nearby AGN.

    Figure 8: Images of our standard model at m. Shown are the inclination angles , , and for a common azimuthal angle of with a logarithmic colour scale. Labels are given in pc.
  • Detailed data comparison: After including part of the discussed additional features, one should do a detailed comparison with observational data. Starting from the simulated surface brightness distributions, we are able to calculate visibilities in order to compare to recent mid-infrared interferometric observations. Furthermore, SINFONI observations (e. g.  Davies et al., 2007) of Seyfert nuclei and starbursting galaxies will give us important constraints for our modelling, e. g. in terms of the nuclear stellar cluster, as well as fuelling of the central active region.

6 Comparison with observations

The only possibility to directly resolve the dust structure in nuclei of Seyfert galaxies so far is with the help of the MID-infrared Interferometric instrument MIDI at the VLTI. Thereby, tori were found in the closest Seyfert nuclei NGC 1068 (Jaffe et al., 2004; Raban et al., 2008) and Circinus (Tristram et al., 2007). Remarkably, both objects are well modelled by a warm, more or less spherical component and a less extended, hotter elongated component, which might be interpreted as a disk seen edge-on. The latter component coincides with an edge-on disk observed in maser emission (Greenhill et al., 2003 for the case of the Circinus galaxy). Some of these interpretations qualitatively match well with the results of our simulations (see discussion in Tristram et al., 2007). It has to be taken into account that our modelling is meant to be suitable for an average Seyfert galaxy.

6.1 Observing our standard model

The basic procedure of our effort of simulating the observational properties of dust tori is subdivided into several steps: First of all, hydrodynamic simulations are carried out with the help of the TRAMP code as described above. In a postprocessing step, the resulting gas density distribution is converted into the corresponding dust density distribution with the help of a constant gas-to-dust-ratio of 250 applied to those patches of the gas distribution, where temperatures are below the sublimation temperatures () of the various grain species. Our dust model takes graphite (K) as well as silicate grains (K) into account. This mixture leads to a prominent spectral feature at 9.7 m and a less pronounced feature at 18.5 m. A slightly higher gas-to-dust ratio compared to the local galactic value of 160 (Sodroski et al., 1994) was applied, as we expect it to increase towards the centres of active galaxies, as dust is likely to be destroyed in the harsh environment close to the energy source. The resulting dust density distributions are then fed into the 3D radiative transfer code MC3D (Wolf et al., 1999; Wolf, 2003). Here, in a two-step approach, first the temperature distribution is calculated, that is then used to simulate observable quantities like spectral energy distributions (SEDs) or images at various wavelengths. Grain-size dependent sublimation (as described in Schartmann et al., 2008) is only considered during the radiative transfer calculations and, therefore, only for dust in the close vicinity of the AGN. Here, we split the dust opacity model into 6 different grains, made up of the smallest (m) and largest (m) grain of the MRN9-size distribution (Mathis et al., 1977) and the three grain species: astronomical silicate and the two orientations of graphite grains (Draine & Lee, 1984; Laor & Draine, 1993; Weingartner & Draine, 2001). As radiative transfer calculations are very time and memory consuming, the hydrodynamic grid – an assembly of all three domains – has to be mapped on to a slightly smaller grid. A resolution of 114 grid cells in radial direction (instead of 181 for the composed hydrodynamic model), 61 grid cells in direction (instead of 94 in the hydrodynamic models) and 128 in -direction (instead of 49 per quadrant in our hydrodynamic models) is currently feasible. We use 96 distinct wavelengths which are logarithmically distributed between 0.001 and 2000 m, with an enhanced resolution in the infrared regime and especially around the silicate features. The gas temperature from the hydrodynamic simulation is only used to decide whether dust is present at a certain grid point () or not (). As we do not follow the dust temperature during the hydrodynamical simulation explicitly, we set the dust temperature to zero prior to the radiative transfer calculation and assume that the dust distribution is solely heated by the central AGN. We assume a bolometric luminosity of the accretion disk of , with a radiation characteristic, preferentially emitting perpendicular to the plane of the disk.

Fig. 8 shows an inclination angle study for images at m of the final stage (see Fig. 3) of our standard model. Shown is the dust re-emission and the central radiation of the almost face-on case ()10, , and the edge-on case () for an azimuthal angle of . Even at an inclination angle of , the central source still possesses the highest surface brightness at these -angles, but it vanishes behind the dense disk component at larger inclination angles (). The dense dust disk is the second brightest characteristic for our models. As expected, it does not appear smooth, but a filamentary structure is visible, as well as a slight nuclear spiral structure (only visible in the high-resolution electronic version), caused by the differential rotation of the disk. Due to the large optical depth at intermediate inclination angles, the torus shows an asymmetric structure with respect to the midplane, as only the upper, directly illuminated funnel wall is visible. The same holds true for the inclination. Here, as well as in the edge-on case, shadows of dense – in most cases radially inwards moving – filaments are visible. The anisotropic radiation source (preferentially emitting along the axis) also contributes to the vertically elongated bright regions and enhances the dark lanes in the midplane.

In Fig. 9, SEDs of the standard model are presented in terms of pure dust re-emission SEDs. Due to the turbulent structure of the dense disk component close to the midplane and the filaments and blobs of high density gas within or close to the funnel region (compare to Fig. 6), the model for the reduction of the silicate feature in emission as described in Schartmann et al. (2008) works fine and only a moderate emission band with a strength of 0.31 is visible at m. For the following analysis, we define the silicate feature strength in accordance with Shi et al. (2006) as


where is the flux at the wavelength of the silicate feature and is the spline fitted continuum, anchored at wavelengths between 5.0 and 7.5 m and from 25.0 to 40.0 m. With increasing inclination angle, it changes to moderate absorption: -0.10 for , -0.41 for and -0.54 for . This is remarkable, as the same average optical depth within the midplane (!!) would cause much deeper absorption – which has never been observed – for the case of a torus with continuous dust distribution. These small feature depths are caused by the fact that the silicate emission feature is produced within an elongated volume along the funnel region. Thus, the appearance of the silicate feature is due to a mixture of emission and absorption and, accordingly, existing channels of low absorption towards the observer (though leading to silicate emission features) can also contribute to the final morphology of the silicate feature.

Figure 9: SEDs of our standard model. Shown are the inclination angles , , and (top to bottom) for a common azimuthal angle of .

For a detailed discussion of the accuracy and limitations of these kind of radiative transfer calculations, we refer the reader to Schartmann et al. (2008). The most severe problems occur in the midplane, where the largest optical depths are reached, which lead to an overestimation of the temperature due to steep gradients. Here, a higher spatial resolution would be desirable. However, we are confident that our generic results are robust, since a resolution study for the case of SEDs yielded almost identical results, as also described in Schartmann et al. (2008).

6.2 Comparison with Spitzer spectra

Figure 10: Comparison of our standard model (solid line – ) with IRS Spitzer observations of NGC 4151 (Weedman et al., 2005) and Mrk 841 (H. W. W. Spoon, private communication).

Figure 10 shows a comparison of our standard hydrodynamic model with observations of two Seyfert galaxies with the Infra-Red Spectrograph (IRS) onboard the Spitzer Space Telescope, in high-sensitivity, low resolution mode. The solid curve refers to an inclination angle of towards our standard model. Shown in colour are the mid-infrared SEDs of the two intermediate type (Sy 1.5 according to NED11, the NASA/IPAC Extragalactic Database) Seyfert galaxies NGC 4151 (blue solid line, Weedman et al., 2005) and Mrk 841 (green, H. W. W. Spoon, private communication). Both are scaled, in order to obtain fluxes similar to our standard model, as we are only interested in a comparison of the shape of the curves. The -inclination angle of our simulated SED shows a self-absorbed silicate feature, yielding a good approximation of the green curve in the respective wavelength range. As we are applying pure dust radiation transfer, we are not able to produce emission lines in our modelled SEDs and we can only compare the underlying continuum flux. Deviations between observations and model are visible in a shift of the maximum of emission of our model towards slightly longer wavelengths and a feature around m. The latter arises most likely from crystalline forsterite in emission (H. W. W. Spoon, private communication), which is not included in our dust model.

As our SEDs do not result from a fitting procedure, this is a very noteworthy result. Especially because the physical parameters we use are meant to represent a typical Seyfert galaxy. The sharp decrease of the flux towards longer wavelengths results from missing cold dust in the outer part.

6.3 Silicate feature strength H i column density relation

Compared to the previously discussed infrared wavelength range, the opacity drops steeply in the X-ray range. Therefore, a certain amount of X-ray photons can escape and provides thus the possibility to determine the hydrogen column density on the line of sight observationally.

For the case of our hydrodynamic models, H i column densities are directly integrated along the line of sight. Only cells are considered, where the temperature is lower than the ionisation temperature of hydrogen. However, in our continuous models to which we will compare our new hydrodynamical models in the following, we cannot take this into account. Therefore, we assume a constant gas-to-dust ratio of 250 (see discussion in Sect. 6.1), in order to calculate column densities from the continuous dust distributions. The strength of the silicate feature is derived as given in formula 6 and is taken at the wavelength, where the maximum distance between the SED and the fitted continuum is reached, as it can be appreciably shifted for the case of self-absorbed emission.

Figure 11: Comparison of resulting silicate feature strengths and column densities of our hydromodels (standard model plus models from a supernova and mass loss rate study) given as stars and the continuous TTM-models (standard model and models from a dust mass study, see Schartmann et al., 2005) given as triangles. The filled stars denote our hydrodynamical standard model. Colour denotes inclination angle (blue – , green – , magenta – , red – , black – intermediate angles of the continuous standard model). The solid lines are linear fits to the observational data shown in Shi et al., 2006 (black – all objects, magenta – Seyfert galaxy sample). The two outlying objects NGC 7213 and PG 1351+640 are given additionally. Ellipses are drawn to guide the eye.

With the help of this procedure, our simulations can be directly compared to the observational data published in Shi et al. (2006). The Shi et al. (2006) sample comprises of 85 AGN of various types, for which Spitzer InfraRed Spectrograph (IRS, Houck et al., 2004) observations or ground-based measurements from the Roche et al. (1991) sample (four objects) exist, as well as X-ray data. H i column densities are obtained from X-ray spectra taken from the Chandra data archive directly for 8 objects. For further 77 objects, column densities were collected from the literature. A correlation was found between the silicate feature strength (see equation 6) and the X-ray absorption column density in a sense that low H i columns correspond to silicate emission and high columns correspond to silicate absorption:


This is shown as the black line in Fig. 11. When Shi et al. (2006) determine the correlation between the silicate feature strength and the X-ray column density data separately for the various types of AGN, a large scatter between them is found, as is also visible in the data range for each type of AGN. Indicated in magenta is the correlation for Seyfert galaxies alone, given by:


In Fig. 11, these correlations are compared to our simulations. Colours of the symbols indicate inclination angle: blue – , green – , magenta – and red – . For the case of inclination angles without any dust along the line-of-sight, we assigned a value of , taking foreground gas into account. Our hydrodynamic models are given as stars and the triangles denote results from our continuous TTM-model simulations (Schartmann et al., 2005). Overlayed ellipses roughly mark the distribution of symbols for the two classes of simulations. Excluding the values at , which are only due to foreground extinction, one can see that the hydrodynamical models agree quite well with the observed correlations, whereas the continuous models show a much too steep dependence. Although the emission features in the Seyfert 1 case are rather too pronounced, they are in good agreement with the two outlying galaxies NGC 7213 and PG 1351+640. Furthermore, the silicate feature changes from emission to absorption within a very small range of column densities and it is also evident that column densities much higher than cannot be reached, as they would cause unphysically high absorption depths of the silicate feature (as already discussed in Schartmann et al., 2008). A further concentration of the dust density towards the minimum of the potential within the midplane might cause a flattening of this distribution, but then, too narrow infrared bumps will be produced, see Schartmann (2003). Obviously, the observations presented in Shi et al. (2006) are in favour of clumpy or filamentary models instead of continuous matter distributions, concerning the models we investigated so far. At least, models are needed which feature high column densities in combination with moderate emission as well as absorption features. In the hydrodynamical models presented in this paper, this is possible due to the clumpy and filamentary nature of the dust distribution, which enables the coexistence of almost free lines of sight towards the inner region, where the silicate feature in emission is produced, together with highly obscured lines of sight. The highest column densities are reached for lines of sight through the dense and turbulent disk in our models. Furthermore, a remarkable consequence of this is the very tight linear correlation of the continuous models and the broad distribution for the case of our (clumpy) hydrodynamical models. Therefore, the broad scatter in the observational data might really be taken as a case for clumpy models, as claimed in Shi et al. (2006). The observed shallow slope of the discussed relation has also been investigated by Levenson et al. (2007) on basis of radiative transfer models of simplified geometries, representing dust clouds and continuous dust distributions. Whereas the latter can yield arbitrarily deep silicate absorption features in their simulations, single dust clouds can not, as long as their dimensions are much smaller than their distance to the illuminating source. With the help of these simulations they conclude that the tori of Seyfert galaxies must be clumpy, as their mean SED shows only moderate emission as well as absorption features and ULIRGs (Ultraluminous Infrared Galaxies) must possess nuclear sources, which are deeply embedded into a geometrically and optically thick and smooth distribution of dust, as their mean SED features deep silicate absorption (Hao et al., 2007). Similar conclusions have been drawn with the help of a clumpy torus model already discussed in section 1 by Nenkova et al. (2002) and Nenkova et al. (2008a, b).

7 Summary

We present a physical model for the evolution of gaseous tori in Seyfert galaxies, in which the stellar evolution of a young nuclear star cluster is predominantly responsible for setting up a torus-like dust (and gas) distribution. The star cluster provides both, (i) mass injection, which is mainly due to planetary nebulae and (ii) injection of energy from supernova explosions. Within the course of the simulations, a highly dynamical multiphase medium evolves, which is comparable to high resolution models of the supernova-driven turbulent ISM in galaxies (Breitschwerdt & de Avillez, 2007; Joung & Mac Low, 2006). There, blast waves from supernovae explosions sweep up the ISM by transferring kinetic energy to the gas, which leads to the subsequent formation of thin shells. These shells interact with one another and with already existing filaments. Radiative cooling then produces cold clumps. In our scenario, the ISM is constantly refilled by stellar mass loss processes distributed according to the stars in the nuclear star cluster. The main structuring agent of the gas and dust in the torus is the interplay between mass input, supernova blast waves locally sweeping away the ISM, the streams of cold gas towards the centre as well as cooling instability.
Within this scenario, cold dense clouds and filaments, where dust can exist, form naturally. The implemented optically thin cooling reduces thermal pressure and, therefore, material streams towards the minimum of the potential. This forms long radial filaments, with lots of substructure due to interaction processes. In agreement with MHD simulations of the ISM by Breitschwerdt & de Avillez (2007), we find absence of pressure equilibrium of the temperature phases. This is the case, as mixing due to supersonic motions and turbulence as well as cooling processes act on shorter time-scales than relaxation processes. Although no pressure equilibrium can be reached, the flow is in a global dynamical equilibrium state. The characteristics of this equilibrium delicately depend on the driving parameters – that are energy input (supernova rate) and mass input rate, as well as the gravitational potential, the cooling curve and rotation. We find completely different velocities and streaming patterns for the various temperature phases. Whereas cold material generally tends to move towards the centre, the radial velocity of the hot ionised medium shows both inflow (in the inner region) and outflow (further out). Depending on the ratio between supernova rate and mass loss rate, the inflow/outflow transition shifts towards smaller radii for higher SN-to-mass injection ratios and towards larger radii, when mass input dominates over supernova heating. For very high supernova rates, we observe solely outflow solutions for the hot gas component. This solution is likely to dominate the initial supernova type II phase shortly after the birth of the stellar population.

The resulting density and temperature structure of the gas is subsequently fed into the radiative transfer code MC3D in order to obtain observable quantities (SEDs and surface brightness distributions in the mid-infrared), which are then compared to recent observations. A comparison with spectral energy data from the Spitzer satellite shows that, despite of its severe limitations in resolution and astrophysics considered within this exploratory study, our hydrodynamical standard model describes one class of Spitzer SEDs fairly well, without the necessity of any finetuning of parameters. For example the two intermediate type Seyfert galaxies NGC 4151 as well as Mrk 841 belong to this class (see Fig. 10). The success of the comparison with the silicate feature strength to Hi column density relation for the hydrodynamical models and the failure of the continuous TTM-models leads to the important conclusion that the presence of a disk-like density enhancement within the equatorial plane in combination with a fluffy structure surrounding it is needed in order to get agreement with observational data. Otherwise, large X-ray column densities cannot simultaneously be obtained together with moderate silicate absorption feature strengths at viewing angles within the midplane. The latter are guaranteed by a mixture of almost dust-free and heavily obscured lines of sight towards the hot inner part of the torus close to the funnel. Therefore, we see a much weaker silicate absorption feature. Due to these remarkable similarities of our exploratory study and observations, we feel encouraged to implement more physical effects into our simulations and to study the central gas and dust structure over a longer evolutionary time of the nuclear star cluster.


We would like to thank the anonymous referee for many useful suggestions to improve the paper, H.W.W. Spoon for providing us the reduced Spitzer SED of Mrk 841 prior to publication, as well as K.R.W. Tristram and M. Krause for helpful discussion. S. W. was supported by the German Research Foundation (DFG) through the Emmy Noether grant WO 857/2.


  1. pagerange: The effect of stellar feedback on the formation and evolution of gas and dust tori in AGNReferences
  2. pubyear: 2008
  3. We assume a starburst duration of 40 Myrs in the following.
  4. As the main mass input at the time of our simulations is given by the ejection of PNe, we further refer to this mass input mechanism as the planetary nebulae injection, although all mass loss effects mentioned above are included.
  5. This is only the case, if the so-called double-degenerate scenario of two merging white dwarfs is realised in nature.
  6. Global orbit always refers to a Keplerian orbit, measured at the torus radius at pc and corresponds to a real time of approximately  yr.
  7. A linear interpolation in time is used.
  8. Ghost cells are additional boundary cells, not belonging to the simulated physical domain, but necessary to implement boundary conditions.
  9. The dust model is labelled according to the initials of Mathis, Rumpl and Nordsieck.
  10. An inclination angle of is used in order to investigate only lines of sight within the computational domain of the hydrodynamical simulations as well as within the torus of our previously investigated 2D radiative transfer simulations we compare to in Sect. 6.3.


  1. Antonucci R., 1993, Ann. Rev. Astron. Astrophys., 31, 473
  2. Antonucci R. R. J., Miller J. S., 1985, Astrophys. J., 297, 621
  3. Balbus S. A., Hawley J. F., 1991, Astrophys. J., 376, 214
  4. Beckert T., Duschl W. J., 2004, Astron. Astrophys., 426, 445
  5. Bender R., Kormendy J., Bower G., Green R., Thomas J., Danks A. C., Gull T., Hutchings J. B., Joseph C. L., Kaiser M. E., Lauer T. R., Nelson C. H., Richstone D., Weistrop D., Woodgate B., 2005, Astrophys. J., 631, 280
  6. Breitschwerdt D., de Avillez M. A., 2007, in Elmegreen B. G., Palous J., eds, IAU Symposium Vol. 237 of IAU Symposium, Dynamical evolution of a supernova driven turbulent interstellar medium. pp 57–64
  7. Bressan A., Fagotto F., Bertelli G., Chiosi C., 1993, Astron. Astrophys. Suppl. Ser., 100, 647
  8. Camenzind M., 1995, Reviews of Modern Astronomy, 8, 201
  9. Courant R., Friedrichs K. O., 1948, Supersonic flow and shock waves. Pure and Applied Mathematics, New York: Interscience, 1948
  10. Dalgarno A., McCray R. A., 1972, Ann. Rev. Astron. Astrophys., 10, 375
  11. Davies R. I., Mueller Sánchez F., Genzel R., Tacconi L. J., Hicks E. K. S., Friedrich S., Sternberg A., 2007, Astrophys. J., 671, 1388
  12. de Avillez M., Breitschwerdt D., 2004, Astrophys. Space. Sci., 292, 207
  13. de Avillez M. A., Breitschwerdt D., 2005, Astron. Astrophys., 436, 585
  14. de Donder E., Vanbeveren D., 2003, New Astronomy, 8, 817
  15. Draine B. T., Lee H. M., 1984, Astrophys. J., 285, 89
  16. Elitzur M., 2006, New Astronomy Review, 50, 728
  17. Elmegreen B. G., 2005, in de Grijs R., González Delgado R. M., eds, ASSL Vol. 329: Starbursts: From 30 Doradus to Lyman Break Galaxies The Initial Mass Function in Starbursts. p. 57
  18. Ferland G. J., 1993, Hazy, A Brief Introduction to Cloudy 84. University of Kentucky Internal Report, 565 pages
  19. Gaibler V., Camenzind M., Krause M., 2005, in Merloni A., Nayakshin S., Sunyaev R. A., eds, Growing Black Holes: Accretion in a Cosmological Context Evolution of the ISM in elliptical galaxies and black hole growth. pp 66–67
  20. Gallimore J. F., Matthews L., 2003, in ASP Conf. Ser. 290: Active Galactic Nuclei: From Central Engine to Host Galaxy NICMOS Observations of the Nuclear Star Cluster of NGC 1068. p. 501
  21. Gerola H., Kafatos M., McCray R., 1974, Astrophys. J., 189, 55
  22. Greenhill L. J., Booth R. S., Ellingsen S. P., Herrnstein J. R., Jauncey D. L., McCulloch P. M., Moran J. M., Norris R. P., Reynolds J. E., Tzioumis A. K., 2003, Astrophys. J., 590, 162
  23. Hao L., Weedman D. W., Spoon H. W. W., Marshall J. A., Levenson N. A., Elitzur M., Houck J. R., 2007, Astrophys. J., Lett., 655, L77
  24. Hönig S. F., Beckert T., Ohnaka K., Weigelt G., 2006, Astron. Astrophys., 452, 459
  25. Houck J. R., Roellig T. L., van Cleve J., Forrest W. J., Herter T., Lawrence C. R., Matthews K., Reitsema H. J., Soifer B. T., Watson D. M., Weedman D., Huisjen M., Troeltzsch J., Barry D. J., Bernard-Salas J., Blacken C. E., Brandl B. R., Charmandaris V., Devost D., Gull G. E., Hall P., Henderson C. P., Higdon S. J. U., Pirger B. E., Schoenwald J., Sloan G. C., Uchida K. I., Appleton P. N., Armus L., Burgdorf M. J., Fajardo-Acosta S. B., Grillmair C. J., Ingalls J. G., Morris P. W., Teplitz H. I., 2004, Astrophys. J., Suppl. Ser., 154, 18
  26. Jaffe W., Meisenheimer K., Röttgering H. J. A., Leinert C., Richichi A., Chesneau O., Fraix-Burnet D., Glazenborg-Kluttig A., Granato G.-L., Graser U., Heijligers B., Köhler R., Malbet F., Miley G. K., Paresce F., Pel J.-W., Perrin G., Przygodda F., Schoeller M., Sol H., Waters L. B. F. M., Weigelt G., Woillez J., de Zeeuw P. T., 2004, Nature, 429, 47
  27. Jenkins E. B., Tripp T. M., 2007, in Elmegreen B. G., Palous J., eds, IAU Symposium Vol. 237 of IAU Symposium, New results on the distribution of thermal pressures in the diffuse ISM. pp 53–56
  28. Joung M. K. R., Mac Low M.-M., 2006, Astrophys. J., 653, 1266
  29. Jungwiert B., Combes F., Palouš J., 2001, Astron. Astrophys., 376, 85
  30. Klahr H. H., 1998, PhD thesis, Friedrich-Schiller-Universität Jena
  31. Klahr H. H., Henning T., Kley W., 1999, Astrophys. J., 514, 325
  32. Kley W., 1989, Astron. Astrophys., 208, 98
  33. Königl A., Kartje J. F., 1994, Astrophys. J., 434, 446
  34. Krolik J. H., 2007, Astrophys. J., 661, 52
  35. Krolik J. H., Begelman M. C., 1988, Astrophys. J., 329, 702
  36. Kwok S., 2005, Journal of Korean Astronomical Society, 38, 271
  37. Laor A., Draine B. T., 1993, Astrophys. J., 402, 441
  38. Leitherer C., Schaerer D., Goldader J. D., Delgado R. M. G., Robert C., Kune D. F., de Mello D. F., Devost D., Heckman T. M., 1999, Astrophys. J., Suppl. Ser., 123, 3
  39. Levenson N. A., Sirocky M. M., Hao L., Spoon H. W. W., Marshall J. A., Elitzur M., Houck J. R., 2007, Astrophys. J., Lett., 654, L45
  40. Levermore C. D., Pomraning G. C., 1981, Astrophys. J., 248, 321
  41. Maciejewski W., 2006, astro-ph/0611259
  42. Marigo P., Bressan A., Chiosi C., 1996, Astron. Astrophys., 313, 545
  43. Mathis J. S., Rumpl W., Nordsieck K. H., 1977, Astrophys. J., 217, 425
  44. Mignone A., Bodo G., Massaglia S., Matsakos T., Tesileanu O., Zanni C., Ferrari A., 2007, Astrophys. J., Suppl. Ser., 170, 228
  45. Nayakshin S., Dehnen W., Cuadra J., Genzel R., 2006, Mon. Not. R. Astron. Soc., 366, 1410
  46. Nenkova M., Ivezić Ž., Elitzur M., 2002, Astrophys. J., Lett., 570, L9
  47. Nenkova M., Sirocky M. M., Ivezic Z., Elitzur M., 2008a, ArXiv e-prints, 806
  48. Nenkova M., Sirocky M. M., Nikutta R., Ivezic Z., Elitzur M., 2008b, ArXiv e-prints, 806
  49. Paumard T., Genzel R., Martins F., Nayakshin S., Beloborodov A. M., Levin Y., Trippe S., Eisenhauer F., Ott T., Gillessen S., Abuter R., Cuadra J., Alexander T., Sternberg A., 2006, Astrophys. J., 643, 1011
  50. Pier E. A., Krolik J. H., 1992, Astrophys. J., Lett., 399, L23
  51. Plewa T., 1995, Mon. Not. R. Astron. Soc., 275, 143
  52. Plummer H. C., 1911, Mon. Not. R. Astron. Soc., 71, 460
  53. Poncelet A., Perrin G., Sol H., 2006, Astron. Astrophys., 450, 483
  54. Prieto M. A., Maciejewski W., Reunanen J., 2005, The Astron. J., 130, 1472
  55. Raban D., Jaffe W., Röttgering H., Meisenheimer K., Tristram K. R. W., 2008, Astron. Astrophys.
  56. Risaliti G., Elvis M., Nicastro F., 2002, Astrophys. J., 571, 234
  57. Ritchmyer R. D., Morton K. W., 1967, Difference methods for initial-value problems. Interscience Tracts in Pure and Applied Mathematics, New York: Interscience, 1967, 2nd ed.
  58. Roche P. F., Aitken D. K., Smith C. H., Ward M. J., 1991, Mon. Not. R. Astron. Soc., 248, 606
  59. Salpeter E. E., 1955, Astrophys. J., 121, 161
  60. Sanders D. B., Phinney E. S., Neugebauer G., Soifer B. T., Matthews K., 1989, Astrophys. J., 347, 29
  61. Scalo J., 1990, in Fabbiano G., Gallagher J. S., Renzini A., eds, ASSL Vol. 160: Windows on Galaxies Top-Heavy IMFs in Starburst Galaxies. p. 125
  62. Scalo J., 1998, in Gilmore G., Howell D., eds, ASP Conf. Ser. 142: The Stellar Initial Mass Function (38th Herstmonceux Conference) The IMF Revisited: A Case for Variations. p. 201
  63. Schartmann M., , 2003, Modelle für Staubtori in Aktiven Galaktischen Kernen
  64. Schartmann M., 2007, PhD thesis, Max-Planck-Institute for Astronomy, Heidelberg, Germany
  65. Schartmann M., Meisenheimer K., Camenzind M., Wolf S., Henning T., 2005, Astron. Astrophys., 437, 861
  66. Schartmann M., Meisenheimer K., Camenzind M., Wolf S., Tristram K. R. W., Henning T., 2008, Astron. Astrophys., 482, 67
  67. Shapiro P. R., Field G. B., 1976, Astrophys. J., 205, 762
  68. Shi Y., Rieke G. H., Hines D. C., Gorjian V., Werner M. W., Cleary K., Low F. J., Smith P. S., Bouwman J., 2006, Astrophys. J., 653, 127
  69. Sodroski T. J., Bennett C., Boggess N., Dwek E., Franz B. A., Hauser M. G., Kelsall T., Moseley S. H., Odegard N., Silverberg R. F., Weiland J. L., 1994, Astrophys. J., 428, 638
  70. Stolte A., Brandner W., Grebel E. K., Lenzen R., Lagrange A.-M., 2005, Astrophys. J., Lett., 628, L113
  71. Stolte A., Grebel E. K., Brandner W., Figer D. F., 2002, Astron. Astrophys., 394, 459
  72. Stone J. M., Mihalas D., Norman M. L., 1992, Astrophys. J., Suppl. Ser., 80, 819
  73. Stone J. M., Norman M. L., 1992a, Astrophys. J., Suppl. Ser., 80, 753
  74. Stone J. M., Norman M. L., 1992b, Astrophys. J., Suppl. Ser., 80, 791
  75. Sullivan M., Le Borgne D., Pritchet C. J., Hodsman A., Neill J. D., Howell D. A., Carlberg R. G., Astier P., Aubourg E., Balam D., Basa S., Conley A., Fabbro S., Fouchez D., Guy J., Hook I., Pain R., Palanque-Delabrouille N., Perrett K., Regnault N., Rich J., Taillet R., Baumont S., Bronder J., Ellis R. S., Filiol M., Lusset V., Perlmutter S., Ripoche P., Tao C., 2006, Astrophys. J., 648, 868
  76. Sutherland R. S., Dopita M. A., 1993, Astrophys. J., Suppl. Ser., 88, 253
  77. Thompson T. A., Quataert E., Murray N., 2005, Astrophys. J., 630, 167
  78. Tristram K. R. W., Meisenheimer K., Jaffe W., Schartmann M., Rix H.-W., 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, Astron. Astrophys., 474, 837
  79. Urry C. M., Padovani P., 1995, Publ. Astron. Soc. Pac., 107, 803
  80. van Leer B., 1977, Journal of Computational Physics, 23, 276
  81. Vázquez G. A., Leitherer C., 2005, Astrophys. J., 621, 695
  82. Wada K., Norman C. A., 2002, Astrophys. J., Lett., 566, L21
  83. Weedman D. W., Hao L., Higdon S. J. U., Devost D., Wu Y., Charmandaris V., Brandl B., Bass E., Houck J. R., 2005, Astrophys. J., 633, 706
  84. Weidemann V., 2000, Astron. Astrophys., 363, 647
  85. Weingartner J. C., Draine B. T., 2001, Astrophys. J., 548, 296
  86. Wolf S., 2003, Computer Physics Communications, 150, 99
  87. Wolf S., Henning T., Stecklum B., 1999, Astron. Astrophys., 349, 839
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description