On the obscuration of the SMBH

On the fate of the matter reinserted within young nuclear stellar clusters

Filiberto Hueyotl-Zahuantitla, Jan Palouš, Richard Wünsch Astronomical Institute, Academy of Sciences of the Czech Republic, Boční II 1401, 141 31 Prague, Czech Republic. Guillermo Tenorio-Tagle22affiliationmark: and Sergiy Silich22affiliationmark: Instituto Nacional de Astrofísica Optica y Electrónica, AP 51, 72000 Puebla, México filibert@asu.cas.cz

This paper presents a hydrodynamical model describing the evolution of the gas reinserted by stars within a rotating young nuclear star cluster (NSC). We explicitly consider the impact of the stellar component to the flow by means of a uniform insertion of mass and energy within the stellar cluster. The model includes the gravity force of the stellar component and a central supermassive black hole (SMBH), and accounts for the heating from the central source of radiation and the radiative cooling of the thermalized gas. By using a set of parameters typical for NSCs and SMBHs in Seyfert galaxies our simulations show that a filamentary/clumpy structure is formed in the inner part of the cluster. This “torus” is Compton thick and covers a large fraction of the sky (as seen from the SMBH). In the outer parts of the cluster a powerful wind is produced, that inhibits the infall of matter from larger scales and thus the NSC-SMBH interplay occurs in isolation.

AGN galaxies: nuclear starburst — black holes — hydrodynamics: accretion — methods: numerical

1 Introduction

The origin of the obscuring matter in active galactic nuclei (AGN) is one of the main challenges in modern astrophysics. The Unified scheme for AGNs requires of a dusty torus to explain the two main types: non-obscured (type 1) and obscured (type 2) AGN (Antonucci 1993; Urry & Padovani 1995). The model claims that the highly energetic central activity of AGNs is powered by mass accretion onto a supermassive black hole (SMBH) with mass , and that the distinction between different types of AGNs is simply determined by the orientation of the torus. The origin, structure, and dynamics of the torus remains as one of the key unsolved problems in AGN physics, which presumably is also related to the SMBH feeding and feedback. Some models assume uniform gas distribution for the obscuring torus, whose thickness is supported by IR radiation pressure (e.g. Pier & Krolik 1992; Krolik 2007, and references therein). Other models assume that the circumnuclear obscuring medium is clumpy, see for example Krolik & Begelman (1988) and Nenkova et al. (2002). However, the mechanism supporting the vertical thickness (in homogeneous models) and motions of the clumps (in clumpy models) is still under debate. Elitzur & Shlosman (2006) suggested that the torus is simply an outflow of dusty and optically thick clumps coming from the accretion disc . Konigl & Kartje (1994) proposed a model where a magneto-centrifugal wind is responsible for the obscuration. Dorodnitsyn et al. (2011) proposed a model in which the obscuration is produced by a dusty wind driven by infrared radiation pressure from a dense torus. The torus as a wind model does not suffer from the vertical structure problem, but the origin of the wind is still unclear. Nayakshin et al. (2012) suggest that the obscuring matter comes from fragmentation of solid bodies (asteroids, comets and terrestrial-like planets) in the vicinity of the SMBH. Wada (2012), based on three dimensional simulations including radiative feedback from the AGN, describe the formation of a turbulent torus from the interaction of back-flows in a bipolar fountain, starting from a preexisting rotationally supported thin disk. Li et al. (2012) provided numerical simulations for the accretion flows with angular momentum which is sufficient to inhibit the accretion if the viscous processes are negligible, and to form a torus from the centrifugal gas. Most of the above models are dedicated mainly to explain the mechanism for vertical support of the torus and they are based on the assumption that a cold gaseous structure already exists. They also neglect the effects that the stellar feedback may provide to the inflow onto the SMBH.

Wada & Norman (2002) and Schartmann et al. (2009) suggested that clumpy tori form with gas reinserted by stars within massive nuclear clusters during the late (post-supernovae) evolution. This model was motivated by the fact that stellar activity in the vicinity of the central SMBH has been found in a variety of Hubble type galaxies (see, for example, Filippenko & Ho 2003; González Delgado et al. 2008; Seth et al. 2008; Kormendy & Bender 2009, and references therein). Such nuclear starbursts are compact ( pc) and have masses in the range (Davies et al. 2007; Watabe et al. 2008; Seth et al. 2010). Some of them show complicated star formation histories (Walcher et al. 2006) and present evidence for global rotation, up to 45 km s, and stellar velocity dispersion of up to 120 km s. The Schartmann et al. (2009) model seems to be in agreement with Davies et al. (2007) and Wild et al. (2010) who claim that the accretion rates and thus the AGN luminosities rise rapidly at the late stage of the nuclear starburst evolution. However, this model does not explain the coincidence of luminous AGNs with young (ages less than 40 Myr) nuclear starbursts, as it is the case of NGC 1097 (see Figure 11 in Davies et al. 2007). Here we show that massive and compact young NSCs with a central SMBH can form filamentary/clumpy gaseous tori, and that the sizes of such tori are in the range inferred by near-infrared observations of Seyfert galaxies (Jaffe et al. 2004; Prieto et al. 2005) and with the sizes derived from their spectral energy distributions (SEDs, Alonso-Herrero et al. 2011, and references therein). These tori are consistent also with the compact sizes of the Compton thick medium estimated from X-ray observations of Seyfert 2 galaxies, type 2 AGNs (Risaliti et al. 1999).

The paper is organized as follows. Section 2 deals with all the ingredients of the hydrodynamical model. Details regarding the numerical scheme are given in section 3. Section 4 describes the hydrodynamic solution which leads to the formation of the torus. We start the discussion with a model without the central radiation field. Then we include the central source of radiation and calculate column densities and possible obscuration fraction of the torus. We also compare analytic estimates of the centrifugal barrier with the numerical results. Section 5 discusses the impact of the NSC wind on the host galaxy and section 6 presents our conclusions.

2 The physical model

We consider the accretion flow onto a SMBH embedded into a young stellar cluster which rotates like a solid body. The key point for our model is that the matter reinserted within young, massive, and compact star clusters could evolve in a catastrophic cooling regime which is very different from the Chevalier & Clegg (1985) adiabatic solution (see Silich et al. 2004; Tenorio-Tagle et al. 2007; Wünsch et al. 2008, hereafter S04, TT07 and W08, respectively). In the Chevalier and Clegg solution the thermalization of the matter ejected by massive stars inside the cluster leads to a high central pressure with an outward pressure gradient that steadily accelerates the gas from zero velocity at the center (i.e., the stagnation point is at the center) to the sound speed at the cluster edge. The reinserted matter exits the cluster as a free wind approaching its terminal velocity, which is twice the sound speed at the cluster border. If radiative cooling is considered, the gas in the central zones of massive clusters cools down and eventually becomes thermally unstable. This is because the average density of the gas increases linearly with the cluster mass while the cooling rate inside the cluster volume grows as a square function of the stellar cluster mass. Consequently, the central pressure drops and the gas cannot be accelerated outwards. The stagnation point moves out of the cluster center and the solution becomes bimodal: within the stagnation volume the thermal instability leads to mass accumulation, while in the outer parts of the cluster a stationary cluster wind is established. The size of the stagnation zone becomes larger as one considers more massive clusters where strong radiative energy losses favor the frequent generation of cold parcels of gas (see TT07 for the semi-analytic procedure for calculating the stagnation radius). This leads to mass accumulation and eventually to star formation and thus to a positive star formation feedback (see Tenorio-Tagle et al. 2005).

In the presence of a central SMBH the stagnation point is always out from the center (Silich et al. 2008; Hueyotl-Zahuantitla et al. 2010, here after S08 and HZ10, respectively) and thus the flow in the vicinity of the central SMBH is always bimodal. Matter inserted within the stagnation volume forms an accretion flow whereas the mass inserted between the stagnation radius and the cluster edge drives the cluster wind. In such cases, the stagnation radius, , is defined by the balance between the outward pressure gradient (strongly affected by radiative cooling) and the gravity force that makes the accretion flow very different from the classic Bondi solution. All matter reinserted within the stagnation zone remains bound to the cluster and thus defines the upper limit to the mass accretion rate onto the SMBH. In thermally unstable bimodal cases, as radiative cooling becomes more important, the stagnation point moves to a larger radius increasing substantially the mass accretion rate. Here by using typical values for the masses and sizes of NSCs, and masses of the SMBH in Seyfert galaxies, we explore the formation of the torus from the mass inserted within a rotating young NSCs with a central SMBH.

Our physical model consists of a spherically-symmetric young NSC of radius and mass with a homogeneous distribution of stars and with a central SMBH of mass . It accounts for the gravity pull from the stellar component and from the black hole. We consider a constant injection of mechanical energy () and mass () within the cluster volume via SNe II and stellar winds. Matter is inserted within the star cluster with a finite angular momentum given as the solid-body rotation around the polar axis: , where is the angular frequency, and are the radial and polar coordinates, respectively, and is the unit vector in the direction . We assume that the rotation velocity of the star cluster is small compared to the dispersion velocity of individual stars and thus we disregard the star cluster flattening due to its rotation. Radiative cooling is one of the main ingredients of the model and it is considered in all the computational domain.

One of the main features of AGNs is their strong emission of ionizing radiation. A proper treatment of the emission from the central source requires to know the intensity and the spectral energy distribution of the radiation field, and how it propagates through the ambient medium. Here we consider only the X-ray luminosity , since such photons can penetrate deeper into high density regions shielded from UV photons. We consider a constant Eddington ratio during the calculations ( erg s). This value is in the range 0.01–0.1 used by Wada (2012). In this way, the model includes Compton and X-ray heating from the central source and the radial component of the acceleration due to radiation pressure. These quantities were explicitly calculated using the ray-tracing method described in Appendix A.

The model implicitly accounts for shock heating by assuming full thermalization of stellar winds and SNe kinetic energy due to random interactions of the ejecta from massive stars that leads to gas temperatures of a few K. The reinserted gas at such temperature and relatively high density is not in thermal equilibrium.

In all calculations we assume a maximum rotation velocity along the equator 50 km s at the star cluster edge. The energy and mass deposition rates ( and ) relate to the wind adiabatic terminal speed () by . The second expression arises when one scales the average results from Starburst99 (SB99, Leitherer et al. 1999). The wind adiabatic terminal speed is an input parameter in the model, we assume in all cases that it is constant and equal to 1000 km s. Note that this value is 2.5 times lower than the average adiabatic wind terminal speed (see Wünsch et al. 2011) for an instantaneous starbursts during the first 40 Myr. Therefore, km s implicitly means additional mass loading to the flow, which may be due to evaporation and destruction of preexisting high density molecular clouds and filaments, and/or evaporation of circumstellar discs forming low mass stars, see for example Stevens & Hartwell (2003); Melioli & de Gouveia Dal Pino (2004). We assume that mass loading is five times that inserted by massive stars.

3 The numerical approach

The numerical models presented here are based on the finite-difference Eulerian hydrodynamic code ZEUS-3D version 3.5 (Stone & Norman 1992; Clarke 2010). All calculations were performed in spherical coordinates in 2D, with symmetry along the - direction. The set of hydrodynamic equations is


and it is closed by the equation of state , where is the internal energy density and the adiabatic index. The total internal energy is . The mass and energy deposition rates per unit volume are and , respectively. The parameter represents mass loading. A small value of means that the mass in winds of individual stars is loaded by additional mass from the parental cluster. The magnitude of the local acceleration due to gravity is , where is the mass enclosed within a sphere of radius . The terms and in equations (2) and (3), respectively, represent the radial acceleration due to radiation pressure and the heating rate per unit volume due to the X-rays from the central source, the method used to calculate these terms is described in Appendix A. The cooling rate per unit volume is , where is the gas number density, is the proton mass, and is the cooling function, which depends on temperature and metallicity . In all cases solar metallicity is assumed. To compute we use an approximate treatment for the ionization degree: we take the values and as the mean mass per particle for ionized and neutral gas respectively.

The model accounts for extremely fast cooling in all the computational domain and is considered in the calculation of the time step. The cooling rate is computed with an updated routine of W08 that uses, for temperatures above the ionization temperature of hydrogen, the Raymond & Cox cooling function tabulated by Plewa (1995), and for T K the Koyama & Inutsuka (2002) cooling function: erg cm s, where erg s. The minimum allowed temperature in the simulations is assumed to be T=100 K. Note that according to Joung & Mac Low (2006) the gas is thermally stable at temperatures where the slope of the cooling curve in the space is , what led Schartmann et al. (2009) to identify five stable zones in the Plewa (1995) cooling function.

The flow was modeled following the prescription given in W08, which explicitly considers a continuous replenishment of mass and internal energy in all cells within the starburst volume at rates and , respectively. The inserted mass is subject to the gravity force of the SMBH as well as from the NSC, like in HZ10, however here the mass is injected with an angular momentum that corresponds to the solid body rotation of the cluster. In summary, the procedure applied to each cell within the cluster volume at every time step is:

  1. The radial velocity of the flow is updated according to . where and are the local acceleration due to gravity and radiation pressure.

  2. The density and total energy in a given cell are saved to and .

  3. The mass is inserted so that , where is the injected mass per unit volume. The mass is inserted with rotation velocity , along the direction, assuming solid-body rotation for the star cluster: , where is the angular frequency, is the distance from the center, and the angle from the polar axis.

  4. The velocity is corrected so that the momentum is conserved: , where the components of the velocity vector v are and . is the unit vector in the direction .

  5. The internal energy is corrected to conserve the total energy .

  6. The new energy is inserted in a form of internal energy .

  7. The AGN-heating is included by increasing the internal energy in each grid cell by .

In steps 3 and 6, is a random number from the interval (-1,1) generated each time step it is used, and is the relative amplitude of the noise. The inclusion of the noise is necessary to break the spherical symmetry imposed by the initial conditions, an analysis of the effects of such noise is presented in W08 in the case of star cluster winds. Here we used their recommended value in all our simulations.

3.1 Initial and boundary conditions

After an initial relaxation period, the solution reaches a steady state with a quasi-stationary wind blowing from the outer parts of the cluster and mass accumulation at an approximately constant rate in the inner region. To make the transition as short as possible, the models start from the stationary adiabatic wind solution of a star cluster with mass , adiabatic wind terminal speed km s, and with radius in each case as given in Table 1. The boundary conditions are set open at both -boundaries and reflecting at both -boundaries. We used the scaled grid option in and a uniform one in . The computational domain and the number of zones in each direction were selected such that , which conserves the shape of the zones and provides a higher resolution closer to the center.

No. Rad.
() () (pc) (km s) (pc)
1 1 10 50 N 9.2 2 1.95 0.76 0.92
2 1 10 50 Y 9.2 2 1.95 0.93 0.86
3 1 10 50 Y
4 1 40 50 Y 31.3 10.5 9.7 0.71 0
5 10 40 50 Y

Note. – Summary of the models. The first column marks the label of each model, where model 2 is our reference model; and are the masses of the NSC and the SMBH respectively; is the radius of the star cluster; is the maximum rotation velocity at the star cluster surface; the labels Y and N, indicate whether a model includes the radiation from the central source or not. is the averaged value of the stagnation radius measured in the simulations at a time when the solution is steady; is the radius of the centrifugal barrier, where the subscripts and denote the results from the numerical simulations and the analytic prediction (Eq. 7) respectively, numerically it corresponds to the time averaged radius of the point of maximum density; indicates the time averaged of the fraction of the sky covered by column densities larger than cm, and gives the fraction of the sky covered only by cold gas ( 1500 K) with column density larger than cm.

Table 1: Selected models and results from the simulations.

4 Results

We have calculated a set of models (see Table 1) with NSC and SMBH parameters in the range given by the observations (see for example Figure 19 in Seth et al. (2010)). Models 1 and 2 have a star cluster radius = 10 pc and mass , respectively, both with a SMBH. In model 1, the radiation from the central source is not considered. Hereafter we refer to model 2 as the reference model. In model 3 a star cluster is 10 times less massive than in the first two cases. Models 4 and 5 consider more extended clusters with = 40 pc, and SMBHs of and , respectively. The computational domain extends radially from 0.5 pc to 15 pc in the case of models 1 – 3, and from 1 pc to 60 pc for models 4 and 5. The axial extent goes from to radians in all calculations. The complete set of relevant parameters for the simulated models are given in Table 1.

4.1 The hydrodynamic solution

4.1.1 Case without a central source of radiation (model 1)

We start with a description of model 1, in which the radiation from the central source is not considered. According to TT07 and W08 such NSC evolves in a catastrophic cooling regime. Here it is shown that the inclusion of gravity of the NSC+SMBH and the angular momentum of the inserted matter lead to the formation of a filamentary/clumpy torus.

Figure 1: Filamentary torus resulting from model 1, snapshots at 1 Myr. The torus is composed by a collection of cold filaments and a dense core located at the centrifugal barrier, at 2 pc. In the left column, from top to bottom, the plotted quantities are logarithms of density, temperature and pressure divided by the Boltzmann constant . In the right column, from top to bottom, the panels corresponds to radial velocity, tangential velocity, and logarithm of angular momentum.

Initially, the average gas density in the central part grows rapidly due to the mass deposition by the stellar cluster. Radiative cooling is enhanced because of its squared dependence on density and thus temperature in the densest zones drops. As soon as temperature decreases to approximately K, free-bound and bound-bound transitions become the major cooling processes. Consequently, when the temperature in some region approaches K the cooling rate increases steeply and the thermal instability starts to operate. This lowers the temperature to 10 K, and decreases the pressure by three orders of magnitude from that ( dyne cm) in the hot ( K) gas. The cold regions are compressed by the surrounding hot gas into dense filaments/clumps. After an initial relaxation period, the stagnation radius , which is defined by radiative cooling, remains almost constant. Note that in our model this radius determines the amount of gas that flows inwards and that later accumulates close to the center and that is different from the Bondi radius which is defined by the mass of the central black hole and by the sound speed in the surrounding ISM.

Figure 1 presents frames of the distribution of the hydrodynamical variables in the whole computational domain. The density distribution presents mainly two-phases: hot (few 10 K) gas with low density, and cold (T=10 K) gas in the densest zones in a form of filaments and clumps. The thermal pressure gradient within the stagnation volume ( pc) is not high enough to push the cold gas out from the cluster volume and instead the cold parcels of gas begin to stream toward the center because of the force of gravity. Due to angular momentum conservation, the rotation velocity of the inflowing gas increases as it approaches the center, to about 200 km s at pc. Such fast rotation prevents the gas to flow further inwards and favors the accumulation of mass around the centrifugal barrier. We identify the collection of cold filaments and the dense core at the centrifugal barrier as the torus responsible for the obscuration of the central SMBH. While all this happens in the central part of the cluster, a stationary cluster wind reaching a terminal speed of about 800 km s is well sustained above in all simulated cases. Such winds may prevent the inflow of matter from larger scales onto the NSC. This suggests that NSC-SMBH may evolve in isolation when the feedback from massive stars is highly active.

4.1.2 Reference model (model 2)

The input parameters in this model are the same as in model 1 except that in this case we consider the effect of the X-ray radiation from the central SMBH. Figure 2 presents a sequence of frames of density, temperature, pressure, radial and tangential components of velocity for the reference model. For comparison, the panels in the right most column display the distribution of the corresponding variables at 1 Myr for the case without central source of radiation (model 1). As one can note, the stagnation radius remains the same in both models. This implies that the same amount of mass is accumulated in the central part of the cluster. However, the inner zone is affected by the central radiation field, which prevents the gas from the fast cooling at the very center (1 pc). Therefore in the reference model the gas there remains hot (at K, see the second row in Figure 2). However, the flow at larger radii (1 pc 3 pc) is still dominated by cooling and a substantial amount of cold gas is accumulated there forming a filamentary/clumpy torus with a dense core supported by rotation at the centrifugal barrier. A fraction of this gas is ionized by the SMBH radiation, so, it remains warm (1500 K K). Note that the main difference between models 1 and 2 is due to the thermal pressure of the ionized gas in the central zone but not due to the radiation pressure.

Figure 2: Time sequence of the distribution of the hydrodynamic variables in the reference model. All panels show the plane, with the axis running horizontally. Each column displays frames at the time indicated on the top, for each variable indicated at the end of each row. As a comparison, the last column shows the corresponding variables for model 1 at Myr.

Figure 3 shows the multi-phase medium that results from the simulation in the reference model. We identify five regions in the vs plane which represent different components of the flow. Region 1 corresponds to the wind, whose temperature drops due to the expansion and radiative cooling. One can also find in this region temperatures that coincide with two stable regions (around and K) of the Plewa (1995) cooling function noticed by Schartmann et al. (2009). Region 2 presents hot (K) rarefied gas resulting from shock-shock collisions and radiative heating. Region 3 results from the radiative heating of the surfaces of dense clumps in the torus. This region does not exist in the case without central source of radiation. Region 4 corresponds to the collection of filaments which tend to settle in one of the stable branches of the cooling function. Eventually these cool down to the minimum allowed temperature in the simulation due to the increase of density as gas moves towards the center. Some high density ( cm ) zones are heated up to K by the AGN radiation. These zones are also not present in model 1 which does not include the central source of radiation. Region 5 corresponds to the core of the torus at the centrifugal barrier where most of the mass is concentrated. The majority of the gas is in the cold phase (T1500 K; see Figure 4).

Figure 3: Phase diagram for the standard model. Temperature against number density at Myr. The diagram was gridded into cells with the masses of the points depicted in color. We identify five components of the flow: 1- wind, 2- hot thermalized gas, 3- heated gas in the torus, 4- filaments and clumps, and 5- cold dense core of the torus. The gas tends to settle at stable regions in the cooling curve, in particular at around K, however the squared dependence of the cooling rate on density leads to lower temperatures.

4.1.3 Other models

In model 3, the stellar cluster is 10 times less massive than that in the reference model. Therefore, the density of the inserted gas is an order of magnitude lower. In this case, the central source of radiation heats up all the gas within the cluster volume up to few K and prevents the formation of thermally unstable zones (clumps). Therefore the reinserted matter does not form torus in this case.

In model 4, the stellar cluster has the same mass as in the reference model, however the cluster is more extended: pc. Therefore, the density of the inserted mass is also lower compared to the reference model. Nevertheless, thermal instabilities occur in the densest regions where the accumulated mass forms a torus. However, in this case the torus is composed only of warm gas. On the other hand, some clumps formed close to eventually join the cluster wind and leave the cluster, reducing the mass accumulation rate.

In model 5 we consider an extended cluster as in model 4 but, in this case the SMBH is 10 times more massive and therefore more energetic than in all previous cases. The strong radiation field keeps the matter within the whole cluster hot what does not allow to form a torus, however a powerful wind as in model 3 is generated. Such cases resemble adiabatic calculations given the impact of radiation.

4.2 Column density and obscuration fraction

In the AGN models it is supposed that a torus, uniform or clumpy, blocks the light coming from the accretion disc. The amount of obscuring gas is usually quantified by the column density, i.e. the number of particles per unit area along the line of sight . The optical/UV radiation is strongly attenuated above cm. If , the opacity is high enough to block even X-ray photons, in such a case the AGN is said to be Compton-thick (Comastri 2004, and references therein). There is observational evidence that suggests that a large fraction of AGNs in the local universe are obscured by Compton thick gas (Maiolino et al. 1998; Risaliti et al. 1999; Matt 2000) and that most of them are associated with Seyfert 2 galaxies (Risaliti et al. 1999).

Figure 4 displays for the reference model the column density at different evolutionary times as a function of the viewing angle as seen from the central SMBH. The line of sight along the equator corresponds to . Different colors correspond to different gaseous components. Note that the column density of the warm gas (blue line) is high enough (cm) to block a large fraction of the X-ray radiation, and thus may turn the AGN Compton thick. The cold component (green line) presents gaps and high variability in its covering angle, which implies that UV photons can escape from the torus through the holes in the neutral gas and an observer can see eventually directly to the center. Such events offer a natural explanation to the “mutation” of optically classified Seyfert 2 to Seyfert 1, and vice versa. Aretxaga & Terlevich (1994) and Aretxaga et al. (1999) give examples of such kind of objects.

Figure 4: Column densities as seen from the central SMBH for the reference model. The red line displays the total column density along each line of sight. The magenta line shows the column density for hot gas, T K. The blue line represents the column density for warm phase, 1500 K T K. The column densities of cold matter (T 1500 K) are shown by green lines. The cold gas does not cover all the sky.

Figure 5 presents the fraction of the sky covered by two different column densities as a function of time for models 1, 2 and 4: (red lines) represents the fraction of the sky covered by total column densities cm, (black lines) shows the fraction of the sky covered by cold gas with column density cm. In the case of model 1, the time average values are and . In the reference model, and . Note that is larger in the reference model compared to model 1. The same tendency was found by Wada (2012) for Eddington ratios in the range 0.01-0.1, where more luminous AGNs are obscured over larger solid angles. The fraction in the reference model is reduced due to ionization by the central source. In the case of model 4 the torus is composed only of warm gas, with average after 1.5 Myr.

Figure 5: Fraction of the sky covered by optically thick gas to optical/UV and X-ray radiation as a function of time. Black lines represent the fraction of the sky, , covered by cold gas in the torus with column densities cm. Red lines shows the fraction of the sky, , covered by total column densities cm.

4.3 Comparison with analytic predictions

The thermally unstable gas inserted within the stagnation radius is attracted by the gravity towards the cluster center. Due to its angular momentum it accumulates around the centrifugal barrier, where the rotation balance the gravitational attraction. Here we give the analytic formula for such radius and show that it is in a good agreement with the numerical results.

The radius where the mass accumulates, , is determined by the angular momentum of the matter inserted within the stagnation radius and by the central gravitational potential of the SMBH + NSC. Here we neglect the effect of the radiation pressure. The value of is defined by radiative cooling, which depends on the mass and compactness of the cluster (see TT07 & W08). Thus, in order to estimate the position of the centrifugal barrier we need to know for a given rotation velocity of the star cluster. In the following we derive an analytic relation for by assuming a star cluster in solid-body rotation.

Let us consider the total mass inserted within the stagnation zone at some instant. Then the specific angular momentum of a rotating parcel of gas is , where is the projection of on the equatorial plane, i.e., the distance to the rotation axis. An integration over the mass within the stagnation volume , gives the average specific angular momentum inserted within at some instant:


Then by considering  constant and using spherical coordinates one obtains:


One expects the accumulation of mass at the centrifugal barrier, i.e., at the radius , where corresponds to the Keplerian velocity of the gas orbiting the mass . This leads to the algebraic equation:


The physical solution of equation (6) is:






Thus, is the key parameter to estimate the centrifugal barrier, because it determines the amount of mass and angular momentum inserted within the stagnation zone for a given set of parameters: , , , and . In our calculations, is self-consistently determined. It splits the cluster into two zones: the inner one where the reinserted matter is accumulated and the outer one where the star cluster wind is formed. It is worth to note that is almost independent on the SMBH mass and may have a non-zero value even if . This makes our approach different from all the modified Bondi models used so far to estimate the size of accretion discs or torus, see for example Ulrich (1976); Proga & Begelman (2003); Krumholz et al. (2005); Inogamov & Sunyaev (2010), where the solution depends on the size of the so-called Bondi radius which is a function of the SMBH mass and the sound speed in the surrounding ISM.

Figure 6 shows the analytic predictions for the centrifugal barrier as a function of the maximum rotation velocity of the NSC, i.e., the rotation velocity at and . Different lines correspond to NSCs of different radii, all of them with a central SMBH of . The squares represent the average values in the case of the reference model ( = 10 pc) and model 4 ( = 40 pc). In all cases the mass of the NSC corresponds to that in the reference model (). As one can expect increases monotonically with the assumed rotation velocity of the cluster. Note that in the range of parameters here used, if one considers more extended clusters at a fixed , the absolute value of the stagnation radius is larger, but the ratio is smaller. Therefore, the specific angular momentum inserted is higher and the mass accumulates at a larger distance from the center. If the NSC is very extended, the impact of cooling is less important and even the absolute value of gets smaller, and eventually the NSC could be in a quasi-adiabatic regime111See, Silich et al. (2004) and Tenorio-Tagle et al. (2007) for a discussion on the threshold energy which separates star clusters evolving in the catastrophic cooling regime from those evolving in a quasi-adiabatic regime.. In such cases tends to a very small value (S04, TT07, W08) and is mainly defined by the gravitational potential (S08 and HZ10).

Figure 6: Analytic prediction for the centrifugal barrier. The radius where mass accumulates as a function of the maximum rotation velocity of star clusters of mass and central SMBH of mass , for different star cluster radius. The squares represent the reference model (model 2) and model 4, see Table 1.
() () (pc) ( yr)
1 1 10 9.2 31.2 7.3 22.1 1.8
2 1 10 9.2 31.2 7.2 23.9 0
4 1 40 31.3 31.2 15.9 14.6 0

Note. – Mass flow rates for models with typical parameters of NSC and SMBH in Seyfert type galaxies. A substantial amount of the total inserted mass () remains locked within the stagnation radius (), and accumulates at a rate in the torus. Note that a good fraction of the total injected mass blows out in the NSC wind at rate . The inflow rate through the inner boundary of the computational domain is denoted by .

Table 2: Mass accumulation rate.

4.4 Mass accumulation rate

In all simulations the hydrodynamic solution reaches a steady state. In the case of our reference model it happens after 0.1 Myr and at about four times longer time in the case of 40 pc clusters. From then onwards the matter inserted in the region flows through and leaves the cluster as a stationary wind which stops the income of matter from a large scale in the galaxy. On the other hand, the mass that remains locked within the stagnation volume streams toward the center and accumulates around the centrifugal barrier practically at a constant rate. Table 2 presents the corresponding rates.

Figure 7 presents for models in Table 2, the absolute and relative (normalized to the star cluster mass input rate) values of the rates of mass deposition within the volume of the cluster (), mass accumulation in the torus (), mass carried away by the wind (), and mass inflow towards the center through the inner zone of the computational domain (). In Model 1: of the total inserted mass accumulates in the torus, leaves the cluster as a wind. About 6% of the inserted mass escape through the inner zone of the computational domain. In the reference model (model 2): 77% of the total inserted mass goes into the torus and 23 % goes into the cluster wind. Model 4: approximately one half of the inserted mass leaves the cluster as a wind, the rest accumulates in the torus. In this case some clumps escape from the cluster, producing peaks in (dashed line) with the corresponding response in (solid line). Note that in models 2 and 4, radiation pressure prevents the inflow of mass through the inner boundary. The actual value of the accretion rate onto the SMBH is beyond the scope of this paper as the inflow of gas to the central black hole is inhibited by angular momentum when the viscous processes are negligible, Li et al. (2012). Note that in all cases the rate of mass accumulation is given by .

The mass of the torus at a given time can be estimated from the average . For example, at 1 Myr it reaches in the reference model, and in the case of model 4 for the same period. Due to the additional mass loading considered in our models, the mass of the torus grows substantially. Such torus is gravitationally unstable. An estimate of the Toomre parameter, , for the obscuring structure in our reference model in the region centered at the centrifugal barrier ( pc) with  pc, sound speed 1 km s, Keplerian rotation at frequency s and surface density g cm results in . Therefore, the mass accumulation should lead to a continuous star formation, which may amplify the effect of the stellar feedback in the nuclear region of the host galaxy.

Figure 7: Mass deposition, accumulation and outflow rates. Left and right axis show scales of absolute and relative quantities, respectively. Dotted lines represent , thick solid lines display , dashed lines shows , and thin solid lines represent . The average values are given in Table 2.

5 Impact of the NSC wind on the host galaxy

In all cases presented in Table 1, the star cluster wind is sufficiently powerful as to significantly re-structure the host galaxy ISM leading perhaps to a thick ring along the plane of the galaxy, and to a super galactic wind along the host galaxy symmetry axis (as in Tenorio-Tagle & Munoz-Tunon 1997, 1998).

A simple estimate of the wind power can be obtained from its ram pressure at the starburst edge. This is in all cases many orders of magnitude larger than the typical ISM pressure in our Galaxy ( dyne cm). For example, for models 1 and 2, dyne cm; and about half this value in model 3. In models 4 and 5, 1.4 and 2.5 dyne cm, respectively. Such values are comparable with the ram pressure of the freely falling gas, where is the free fall velocity, in the torus formation model of Hopkins et al. (2012) with the inflow rate yr ( dyne cm for models with 10 pc radius and dyne cm for models with 40 pc). Thus the NSC winds have to build up super bubbles and probably super galactic winds preventing, in most cases, the falling of the ISM onto the NSC.

6 Conclusions

We present 2D radiation hydrodynamic simulations of the gas reinserted by stars of a rotating young NSC with a central SMBH. Our model considers explicitly the impact which the stars from a NSC provide on the accretion flow. The model includes constant mass and energy deposition from stars and assumes that the mass is inserted with non zero angular momentum. It takes into consideration gravity from the central SMBH and from the NSC, and accounts for radiative cooling and the heating from a central isotropic source of X-ray radiation.

Here we have shown that the combined effect of gravity from the SMBH+NSC and the angular momentum of the inserted mass results in the formation of a dense structure (torus) inside the NSC, well within the stagnation radius , defined by radiative cooling. The torus is only a few parsecs across, filamentary/clumpy, with a core at the centrifugal barrier. It is composed of gas in two phases: a cold phase (T K), where dust can survive as close as a couple of parsecs from the SMBH, and a warm phase (1500 KTK), maintained at this temperature by heating from the central source of radiation. The torus is Compton thick and covers a large fraction of the sky, more than 80% in our reference model. This obscuring structure is embedded into a low density hot gas.

Note that models developed by Wada & Norman (2002) and Schartmann et al. (2009) and that are here discussed, lead to the flow of the cold gas towards the central zone of the host galaxy. The inflow of molecular gas in the inner pc of the Seyfert galaxy NGC 4051 was detected by Riffel et al. (2008) who suggested that the inflow occurs due to the spiral arms that reach the nucleus of the galaxy. Grier et al. (2012) presented the evidence for the inflowing gas in the broad line regions of four AGNs: Mrk 335, Mrk 1501, 3C 120 and PG2130+099. Sim et al. (2012) showed that the parsec-scale inflows do not result in significant absorption features in the X-ray spectra since the ionization degree of the infalling gas is high. Therefore, the lack of such observations does not rule out our model.

The accreting mass accumulates in the central region almost at a constant rate, resulting after some time in a very massive torus. As soon as it becomes gravitationally unstable a second generation of stars may form there leading to the formation of a stellar torus, and thus matter would be continuously reprocessed into stars at a rate dictated by the mass accumulation.

A necessary, but not sufficient, condition for the formation of the torus is that the matter reinserted within the NSC evolves in a thermally unstable regime. However, the formation of the torus may be prevented by the strong central source of radiation as it is the case in our models 3 and 5.

In all cases a powerful cluster wind is established outside the stagnation radius. Such winds can inhibit the income of gas from larger scale in the galaxy. This suggests that during the starburst phase, when massive stars dominate the NSC feedback to the host galaxy ISM, the NSC-SMBH interplay occurs in isolation.

The authors express their thanks to the anonymous referee, whose constructive comments and proposal have helped improve this paper. This study has been supported by the Czech Science Foundation grant 209/12/1795 and by the project RVO: 67985815; the Academy of Sciences of the Czech Republic and CONACYT-México research collaboration under the project 17048: Violent star formation; the CONACYT - México, research grants 131913 and 167169. H-Z. F. wishes to express his thanks to CONACYT-México for additional support through grants 162184 and 186720.

Appendix A The heating rate and acceleration due to the central radiation field

We consider only the high frequency band of an AGN spectrum, which allows to take the advantage of a parametric form of the Compton () and X-ray () heating functions given by Blondin (1994). By means of a ray tracing we calculate the optical depth at each radius , where is the inner boundary of the computational domain, is the local gas density, and is the radial length of a grid cell. We assume that the attenuation in the ionized gas is dominated by Thompson scattering, i.e. the opacity is cm g, and two orders of magnitude higher (see Proga et al. 2000) for the neutral gas.

The optical depth is used to compute the X-ray flux , which is required to calculate the heating rate and the acceleration due to the radiation pressure. The amount of energy emitted in X-rays depends on the SED and the total luminosity of the source. For the average AGN SED given in Korista et al. (1997), one gets that about 8% of the total luminosity is emitted in X-rays. Here we take this fraction and assume that the total luminosity corresponds to the Eddington limit for a given SMBH, therefore, we use .

Once is known, the local acceleration due to the radiation pressure in equation (2) is computed as , where is the speed of light. Then the radial velocity of the flow is corrected by .

The heating rate per unit volume is , it is included in equation (3) by increasing the internal energy in each grid cell by . Explicitly, and , both in units of erg s cm. Such functions depend on the temperature of the gas, the characteristic temperature of 10 keV X-ray radiation, and on the ionization parameter (Tarter et al. 1969). The last parameter characterizes states of ionization equilibrium and depends on the local flux and the number density of particles within a grid cell.


  • Alonso-Herrero et al. (2011) Alonso-Herrero, A., Ramos Almeida, C., Mason, R., et al. 2011, ApJ, 736, 82
  • Antonucci (1993) Antonucci, R. 1993, ARA&A, 31, 473
  • Aretxaga et al. (1999) Aretxaga, I., Joguet, B., Kunth, D., Melnick, J., & Terlevich, R. J. 1999, ApJ, 519, L123
  • Aretxaga & Terlevich (1994) Aretxaga, I. & Terlevich, R. 1994, in IAU Symposium, Vol. 159, Multi-Wavelength Continuum Emission of AGN, ed. T. Courvoisier & A. Blecha, 438
  • Blondin (1994) Blondin, J. M. 1994, ApJ, 435, 756
  • Chevalier & Clegg (1985) Chevalier, R. A. & Clegg, A. W. 1985, Nature, 317, 44
  • Clarke (2010) Clarke, D. A. 2010, ApJS, 187, 119
  • Comastri (2004) Comastri, A. 2004, in Astrophysics and Space Science Library, Vol. 308, Supermassive Black Holes in the Distant Universe, ed. A. J. Barger, 245–+
  • Davies et al. (2007) Davies, R. I., Müller Sánchez, F., Genzel, R., et al. 2007, ApJ, 671, 1388
  • Dorodnitsyn et al. (2011) Dorodnitsyn, A., Bisnovatyi-Kogan, G. S., & Kallman, T. 2011, ApJ, 741, 29
  • Elitzur & Shlosman (2006) Elitzur, M. & Shlosman, I. 2006, ApJ, 648, L101
  • Filippenko & Ho (2003) Filippenko, A. V. & Ho, L. C. 2003, ApJ, 588, L13
  • González Delgado et al. (2008) González Delgado, R. M., Pérez, E., Cid Fernandes, R., & Schmitt, H. 2008, AJ, 135, 747
  • Grier et al. (2012) Grier, C. J., Peterson, B. M., Horne, K., et al. 2012, ArXiv e-prints
  • Hopkins et al. (2012) Hopkins, P. F., Hayward, C. C., Narayanan, D., & Hernquist, L. 2012, MNRAS, 420, 320
  • Hueyotl-Zahuantitla et al. (2010) Hueyotl-Zahuantitla, F., Tenorio-Tagle, G., Wünsch, R., Silich, S., & Palouš, J. 2010, ApJ, 716, 324
  • Inogamov & Sunyaev (2010) Inogamov, N. A. & Sunyaev, R. A. 2010, Astronomy Letters, 36, 835
  • Jaffe et al. (2004) Jaffe, W., Meisenheimer, K., Röttgering, H. J. A., et al. 2004, Nature, 429, 47
  • Joung & Mac Low (2006) Joung, M. K. R. & Mac Low, M.-M. 2006, ApJ, 653, 1266
  • Konigl & Kartje (1994) Konigl, A. & Kartje, J. F. 1994, ApJ, 434, 446
  • Korista et al. (1997) Korista, K., Baldwin, J., Ferland, G., & Verner, D. 1997, ApJS, 108, 401
  • Kormendy & Bender (2009) Kormendy, J. & Bender, R. 2009, ApJ, 691, L142
  • Koyama & Inutsuka (2002) Koyama, H. & Inutsuka, S.-i. 2002, ApJ, 564, L97
  • Krolik (2007) Krolik, J. H. 2007, ApJ, 661, 52
  • Krolik & Begelman (1988) Krolik, J. H. & Begelman, M. C. 1988, ApJ, 329, 702
  • Krumholz et al. (2005) Krumholz, M. R., McKee, C. F., & Klein, R. I. 2005, ApJ, 618, 757
  • Leitherer et al. (1999) Leitherer, C., Schaerer, D., Goldader, J. D., et al. 1999, ApJS, 123, 3
  • Li et al. (2012) Li, J., Ostriker, J., & Sunyaev, R. 2012, ArXiv e-prints
  • Maiolino et al. (1998) Maiolino, R., Salvati, M., Bassani, L., et al. 1998, A&A, 338, 781
  • Matt (2000) Matt, G. 2000, A&A, 355, L31
  • Melioli & de Gouveia Dal Pino (2004) Melioli, C. & de Gouveia Dal Pino, E. M. 2004, A&A, 424, 817
  • Nayakshin et al. (2012) Nayakshin, S., Sazonov, S., & Sunyaev, R. 2012, MNRAS, 419, 1238
  • Nenkova et al. (2002) Nenkova, M., Ivezić, Ž., & Elitzur, M. 2002, ApJ, 570, L9
  • Pier & Krolik (1992) Pier, E. A. & Krolik, J. H. 1992, ApJ, 401, 99
  • Plewa (1995) Plewa, T. 1995, MNRAS, 275, 143
  • Prieto et al. (2005) Prieto, M. A., Maciejewski, W., & Reunanen, J. 2005, AJ, 130, 1472
  • Proga & Begelman (2003) Proga, D. & Begelman, M. C. 2003, ApJ, 582, 69
  • Proga et al. (2000) Proga, D., Stone, J. M., & Kallman, T. R. 2000, ApJ, 543, 686
  • Riffel et al. (2008) Riffel, R. A., Storchi-Bergmann, T., Winge, C., et al. 2008, MNRAS, 385, 1129
  • Risaliti et al. (1999) Risaliti, G., Maiolino, R., & Salvati, M. 1999, ApJ, 522, 157
  • Schartmann et al. (2009) Schartmann, M., Meisenheimer, K., Klahr, H., et al. 2009, MNRAS, 393, 759
  • Seth et al. (2008) Seth, A., Agüeros, M., Lee, D., & Basu-Zych, A. 2008, ApJ, 678, 116
  • Seth et al. (2010) Seth, A. C., Cappellari, M., Neumayer, N., et al. 2010, ApJ, 714, 713
  • Silich et al. (2008) Silich, S., Tenorio-Tagle, G., & Hueyotl-Zahuantitla, F. 2008, ApJ, 686, 172
  • Silich et al. (2004) Silich, S., Tenorio-Tagle, G., & Rodríguez-González, A. 2004, ApJ, 610, 226
  • Sim et al. (2012) Sim, S. A., Proga, D., Kurosawa, R., et al. 2012, MNRAS, 426, 2859
  • Stevens & Hartwell (2003) Stevens, I. R. & Hartwell, J. M. 2003, MNRAS, 339, 280
  • Stone & Norman (1992) Stone, J. M. & Norman, M. L. 1992, ApJS, 80, 753
  • Tarter et al. (1969) Tarter, C. B., Tucker, W. H., & Salpeter, E. E. 1969, ApJ, 156, 943
  • Tenorio-Tagle & Munoz-Tunon (1997) Tenorio-Tagle, G. & Munoz-Tunon, C. 1997, ApJ, 478, 134
  • Tenorio-Tagle & Munoz-Tunon (1998) Tenorio-Tagle, G. & Munoz-Tunon, C. 1998, MNRAS, 293, 299
  • Tenorio-Tagle et al. (2005) Tenorio-Tagle, G., Silich, S., Rodríguez-González, A., & Muñoz-Tuñón, C. 2005, ApJ, 628, L13
  • Tenorio-Tagle et al. (2007) Tenorio-Tagle, G., Wünsch, R., Silich, S., & Palouš, J. 2007, ApJ, 658, 1196
  • Ulrich (1976) Ulrich, R. K. 1976, ApJ, 210, 377
  • Urry & Padovani (1995) Urry, C. M. & Padovani, P. 1995, PASP, 107, 803
  • Wada (2012) Wada, K. 2012, ApJ, 758, 66
  • Wada & Norman (2002) Wada, K. & Norman, C. A. 2002, ApJ, 566, L21
  • Walcher et al. (2006) Walcher, C. J., Böker, T., Charlot, S., et al. 2006, ApJ, 649, 692
  • Watabe et al. (2008) Watabe, Y., Kawakatu, N., & Imanishi, M. 2008, ApJ, 677, 895
  • Wild et al. (2010) Wild, V., Heckman, T., & Charlot, S. 2010, MNRAS, 405, 933
  • Wünsch et al. (2011) Wünsch, R., Silich, S., Palouš, J., Tenorio-Tagle, G., & Muñoz-Tuñón, C. 2011, ApJ, 740, 75
  • Wünsch et al. (2008) Wünsch, R., Tenorio-Tagle, G., Palouš, J., & Silich, S. 2008, ApJ, 683, 683
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 minimum 40 characters and the title a minimum of 5 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