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.
keywords:Galaxies: nuclei – Galaxies: Seyfert – ISM: dust, extinction – Radiative transfer – Hydrodynamics – ISM: evolution.
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
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
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 ISM
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
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:
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 ().
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
(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 orbits
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
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 quantities
Concerning the other coordinates, we model a range of with
periodic boundary conditions in
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
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.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.
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
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.
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.
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.
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
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.
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 MRN
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
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.
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 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
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.
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).
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
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.
- pagerange: The effect of stellar feedback on the formation and evolution of gas and dust tori in AGN–References
- pubyear: 2008
- We assume a starburst duration of 40 Myrs in the following.
- 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.
- This is only the case, if the so-called double-degenerate scenario of two merging white dwarfs is realised in nature.
- Global orbit always refers to a Keplerian orbit, measured at the torus radius at pc and corresponds to a real time of approximately yr.
- A linear interpolation in time is used.
- Ghost cells are additional boundary cells, not belonging to the simulated physical domain, but necessary to implement boundary conditions.
- The dust model is labelled according to the initials of Mathis, Rumpl and Nordsieck.
- 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.
- Antonucci R., 1993, Ann. Rev. Astron. Astrophys., 31, 473
- Antonucci R. R. J., Miller J. S., 1985, Astrophys. J., 297, 621
- Balbus S. A., Hawley J. F., 1991, Astrophys. J., 376, 214
- Beckert T., Duschl W. J., 2004, Astron. Astrophys., 426, 445
- 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
- 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
- Bressan A., Fagotto F., Bertelli G., Chiosi C., 1993, Astron. Astrophys. Suppl. Ser., 100, 647
- Camenzind M., 1995, Reviews of Modern Astronomy, 8, 201
- Courant R., Friedrichs K. O., 1948, Supersonic flow and shock waves. Pure and Applied Mathematics, New York: Interscience, 1948
- Dalgarno A., McCray R. A., 1972, Ann. Rev. Astron. Astrophys., 10, 375
- 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
- de Avillez M., Breitschwerdt D., 2004, Astrophys. Space. Sci., 292, 207
- de Avillez M. A., Breitschwerdt D., 2005, Astron. Astrophys., 436, 585
- de Donder E., Vanbeveren D., 2003, New Astronomy, 8, 817
- Draine B. T., Lee H. M., 1984, Astrophys. J., 285, 89
- Elitzur M., 2006, New Astronomy Review, 50, 728
- 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
- Ferland G. J., 1993, Hazy, A Brief Introduction to Cloudy 84. University of Kentucky Internal Report, 565 pages
- 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
- 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
- Gerola H., Kafatos M., McCray R., 1974, Astrophys. J., 189, 55
- 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
- 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
- Hönig S. F., Beckert T., Ohnaka K., Weigelt G., 2006, Astron. Astrophys., 452, 459
- 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
- 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
- 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
- Joung M. K. R., Mac Low M.-M., 2006, Astrophys. J., 653, 1266
- Jungwiert B., Combes F., Palouš J., 2001, Astron. Astrophys., 376, 85
- Klahr H. H., 1998, PhD thesis, Friedrich-Schiller-Universität Jena
- Klahr H. H., Henning T., Kley W., 1999, Astrophys. J., 514, 325
- Kley W., 1989, Astron. Astrophys., 208, 98
- Königl A., Kartje J. F., 1994, Astrophys. J., 434, 446
- Krolik J. H., 2007, Astrophys. J., 661, 52
- Krolik J. H., Begelman M. C., 1988, Astrophys. J., 329, 702
- Kwok S., 2005, Journal of Korean Astronomical Society, 38, 271
- Laor A., Draine B. T., 1993, Astrophys. J., 402, 441
- 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
- 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
- Levermore C. D., Pomraning G. C., 1981, Astrophys. J., 248, 321
- Maciejewski W., 2006, astro-ph/0611259
- Marigo P., Bressan A., Chiosi C., 1996, Astron. Astrophys., 313, 545
- Mathis J. S., Rumpl W., Nordsieck K. H., 1977, Astrophys. J., 217, 425
- Mignone A., Bodo G., Massaglia S., Matsakos T., Tesileanu O., Zanni C., Ferrari A., 2007, Astrophys. J., Suppl. Ser., 170, 228
- Nayakshin S., Dehnen W., Cuadra J., Genzel R., 2006, Mon. Not. R. Astron. Soc., 366, 1410
- Nenkova M., Ivezić Ž., Elitzur M., 2002, Astrophys. J., Lett., 570, L9
- Nenkova M., Sirocky M. M., Ivezic Z., Elitzur M., 2008a, ArXiv e-prints, 806
- Nenkova M., Sirocky M. M., Nikutta R., Ivezic Z., Elitzur M., 2008b, ArXiv e-prints, 806
- 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
- Pier E. A., Krolik J. H., 1992, Astrophys. J., Lett., 399, L23
- Plewa T., 1995, Mon. Not. R. Astron. Soc., 275, 143
- Plummer H. C., 1911, Mon. Not. R. Astron. Soc., 71, 460
- Poncelet A., Perrin G., Sol H., 2006, Astron. Astrophys., 450, 483
- Prieto M. A., Maciejewski W., Reunanen J., 2005, The Astron. J., 130, 1472
- Raban D., Jaffe W., Röttgering H., Meisenheimer K., Tristram K. R. W., 2008, Astron. Astrophys.
- Risaliti G., Elvis M., Nicastro F., 2002, Astrophys. J., 571, 234
- 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.
- Roche P. F., Aitken D. K., Smith C. H., Ward M. J., 1991, Mon. Not. R. Astron. Soc., 248, 606
- Salpeter E. E., 1955, Astrophys. J., 121, 161
- Sanders D. B., Phinney E. S., Neugebauer G., Soifer B. T., Matthews K., 1989, Astrophys. J., 347, 29
- 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
- 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
- Schartmann M., , 2003, Modelle für Staubtori in Aktiven Galaktischen Kernen
- Schartmann M., 2007, PhD thesis, Max-Planck-Institute for Astronomy, Heidelberg, Germany
- Schartmann M., Meisenheimer K., Camenzind M., Wolf S., Henning T., 2005, Astron. Astrophys., 437, 861
- Schartmann M., Meisenheimer K., Camenzind M., Wolf S., Tristram K. R. W., Henning T., 2008, Astron. Astrophys., 482, 67
- Shapiro P. R., Field G. B., 1976, Astrophys. J., 205, 762
- 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
- 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
- Stolte A., Brandner W., Grebel E. K., Lenzen R., Lagrange A.-M., 2005, Astrophys. J., Lett., 628, L113
- Stolte A., Grebel E. K., Brandner W., Figer D. F., 2002, Astron. Astrophys., 394, 459
- Stone J. M., Mihalas D., Norman M. L., 1992, Astrophys. J., Suppl. Ser., 80, 819
- Stone J. M., Norman M. L., 1992a, Astrophys. J., Suppl. Ser., 80, 753
- Stone J. M., Norman M. L., 1992b, Astrophys. J., Suppl. Ser., 80, 791
- 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
- Sutherland R. S., Dopita M. A., 1993, Astrophys. J., Suppl. Ser., 88, 253
- Thompson T. A., Quataert E., Murray N., 2005, Astrophys. J., 630, 167
- 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
- Urry C. M., Padovani P., 1995, Publ. Astron. Soc. Pac., 107, 803
- van Leer B., 1977, Journal of Computational Physics, 23, 276
- Vázquez G. A., Leitherer C., 2005, Astrophys. J., 621, 695
- Wada K., Norman C. A., 2002, Astrophys. J., Lett., 566, L21
- 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
- Weidemann V., 2000, Astron. Astrophys., 363, 647
- Weingartner J. C., Draine B. T., 2001, Astrophys. J., 548, 296
- Wolf S., 2003, Computer Physics Communications, 150, 99
- Wolf S., Henning T., Stecklum B., 1999, Astron. Astrophys., 349, 839