- 1Institute of Radio Astronomy, National Academy of Science of Ukraine, Kharkiv, Ukraine
- 2Institute of Astronomy, V. N. Karazin Kharkiv National University, Kharkiv, Ukraine
We considered the evolution of a self-gravitating clumpy torus in the gravitational field of the central mass of an active galactic nucleus (AGN) in the framework of the N-body problem. The initial conditions take into account winds with different opening angles. Results of our N-body simulations show that the clouds moving on orbits with a spread in inclinations and eccentricities form a toroidal region. The velocity of the clouds at the inner boundary of the torus is lower than in a disk model that can explain the observed rotation curves. We discuss the scenario of torus formation related with the beginning of the AGN stage.
1. Introduction
A dusty torus is an important structural element of an active galactic nucleus (AGN). In the framework of the unified scheme, the observational properties of AGNs (of type 1 and 2) are explained by the different orientation of the torus relative to an observer. This scheme was applied to Sy-galaxies (Antonucci, 1993) and generalized to other classes of AGNs (Urry and Padovani, 1995). Direct observations of obscuring tori exist only for a the nearby Sy-galaxies NGC 1068 (Jaffe et al., 2004; Raban et al., 2009), Circinus (Tristram et al., 2007, 2014).
The observed spectral energy distributions (SEDs) in the IR correspond to the model of a clumpy thick torus with a Gaussian distribution of clouds in its cross-section (Toroidal Obscuration Region) (Nenkova et al., 2008a,b). The radiative transfer model of a clumpy torus with 3D cloud distributions was considered in (Hönig and Kishimoto, 2010; Stalevski et al., 2012). In the other radiative transfer models the IR emission is explained by two components: disk and wind (Hönig and Kishimoto, 2017 and references therein). ALMA observations allowed to estimate the mass of the torus in NGC 1068: (García-Burillo et al., 2016). It means that self-gravity of the torus can influence the motion of the clouds in it. In fact, the dynamics of the matter in the torus shows non-circular motions which can be related to its self-gravitating properties.
One of the main problem concerns the explanation of the geometrical thickness of the torus. In such a torus the vertical velocity component of the clouds must be of the same order of magnitude as the orbital velocity. Several mechanisms were offered for the solution of this problem. The geometrical thickness of the torus can be explained by IR radiation pressure (Krolik, 2007; Dorodnitsyn et al., 2016; Chan and Krolik, 2017), by turbulent motions (Schartmann et al., 2010), or by starburst processes (Wada, 2012; Wada et al., 2016). Other models propose clumpy or dusty winds as an obscuring region (for example, Elitzur and Shlosman, 2006). Indeed, many AGNs show the presence of outflows (winds) which can appear due to influence of the radiation pressure or of a magnetic field (Proga and Kallman, 2004; Netzer, 2013, and references therein).
Our main idea is that the geometrical thickness of torus can be achieved by the motion of clouds in inclined orbits (Bannikova et al., 2012). This assumption is quite natural because there is the external accretion of matter which replenishes the central region of AGN. On the other hand, since we consider a discrete medium in the torus, the orbital plane of each cloud must pass through the central mass. In this case, the toroidal structure may form due to the motion of the clouds in inclined orbits.
In our previous papers (Bannikova et al., 2012), we used special initial conditions (Keplerian torus) for the investigation of the torus evolution (see also Bannikova, 2015). In a Keplerian torus, the distribution of particles by eccentricities and inclinations obeys some law and all major semi-axes equal a major radius of the torus. Here we continue to investigate the properties and stability of a self-gravitating clumpy torus in the framework of N-body problem for more general initial conditions, taking into account the wind cones with different opening angles. In section 2 we discuss the initial conditions and present the result of N-body simulation. Sections 3 and 4 are devoted to the main results on the cloud dynamics in the torus and to the conditions of obscuration (section 5). A discussion of possible scenario of torus formation is provided in section 6.
2. N-body Simulation of A Torus: Initial Conditions
We consider, as initial conditions, the random distribution of clouds over all orbital elements: eccentricity, inclination, major semi-axis, and three angles. Such an initial distribution is more general than in Bannikova et al. (2012). To form a toroidal structure, an anisotropy in the distribution of the clouds associated with winds in AGNs is needed. In these wind cones, the clouds acquire additional momentum due to radiation pressure and may overcome the gravitational forces of the central mass (and torus), leaving the system. This fact is accounted here by a simple assumption: the clouds from two opposite polar opening angles are excluded. In this case, the half-opening angle of the wind is a parameter influencing the resulting equilibrium cross-section of the torus and the obscuration condition. Three projections of the initial distribution of the clouds for the half-opening angle of the winds, , are shown in Figure 1, (top panels). We will use in the following a value of the half-opening angle of the torus, θ0 = π/2−θwind, and show the results of simulations for different values of θ0 in section 5.
Figure 1. Distribution of the clouds in the torus for the initial state (top panels) and after 1, 000 average orbital periods (bottom panels): projections on the equatorial plane (left), on the meridional plane (center), and with all the clouds gathered in one meridional plane (right). Initial conditions: Mtorus = 0.01MBH, N = 16, 384, ϵ = 0.01, .
N-body problem is reduced to the numerical integration of the equations of motion taking account the gravitational field of the central mass and (N − 1) clouds of constant mass, Mcl = Mtorus/(N − 1):
where ri = (x,y,z) is the vector of a cloud from the central mass normalized to the major radius of the torus (R); ai is the vector of the acceleration of the i-th cloud acquired from all of the clouds of the torus Mtorus and from the central mass MBH. A softening length ϵ in the N-body problem allows us to avoid unlimited increasing of the gravitational forces by collisions of particles and can be interpreted as the radius of the cloud Rcl = ϵR. We choose the physical parameters corresponding to the case of NGC1068: and . Since the method of parallel calculations is used, the number of particle in N-body problem must be N = 2n. In presented simulations we adopt n = 14 (or n = 13), thus the number of clouds comes Ncl = N − 1 = 16, 383 (or Ncl = 8, 191).
N-body simulations were carried out using a technology of parallel calculations with GPU (CUDA). We used the Euler method with a step 0.001 to solve equation (1). In this case the total energy of the system E is a constant with a good accuracy, , where E0 is the initial value of the total energy. We use the unity system: MBH = G = R = 1.
The result of the simulation shows that the torus cross-section is changed and achieves its equilibrium state after a few hundred orbital periods. The distribution of clouds after 1, 000 orbital periods are shown in the bottom panels of Figure 1. (One orbital period corresponds T = 30, 000 years for and R = 1 pc.) It is seen that the clouds are spread along z-axis as compared to the initial state, and that the cloud density increases towards the center of the torus cross-section. Indeed, the torus potential in the N-body problem might be divided into two summands (Bannikova et al., 2012): a regular part which is related with a smooth potential of the torus, and an irregular one which is due to the gravitational interactions of the clouds. The regular potential leads to an increase of the cloud density towards the center of the torus cross-section, but the irregular forces between clouds tend to stretch it. As the result, the self-gravitating torus is geometrically thick, which is needed for the obscuration condition (see section 5).
We can suggest that the behaviour of the system for a larger number of clouds N will not essentially differ from the obtained results. For example, if two toroidal systems have the same mass but different numbers of clouds (N1, N2), the same velocity dispersion of the clouds in the torus will be reached for the second system in the time interval Δt2 = Δt1N2/N1, where Δt1 is the time interval for the first system (Bannikova et al., 2012).
3. Distribution of the Clouds in A Torus
At equilibrium, the clouds are distributed in such a way that they form the toroidal structure: the number of clouds is exponentially decreasing along the z-axis (Figure 1, bottom panels), and the distribution of the clouds in the torus cross-section is Gaussian. This distribution is similar to that obtained by (Nenkova et al., 2008a,b) from the analysis of the spectral energy distribution (SED) in the IR and was also used in 3D radiative transfer model of a clumpy torus (Hönig and Kishimoto, 2010). In our simulation, such a distribution is produced by the gravitational interaction between all the clouds and the central black hole.
Knowing from N-body simulations the coordinates and velocity components, we can calculate the orbital elements of each cloud. Figure 2 shows a comparison of the cloud distribution in the initial state (red curves) and after 1,000 average orbital periods for two values of the softening length ϵ = 0.01 (green curves) and ϵ = 0.005 (blue curves). It can be seen that the resulting distributions for the inclination and the major semi-axis (Figure 2, left, right) do not essentially differ from the initial ones, while there is a little difference in the distribution for eccentricity (Figure 2, center). This result shows that the initial state is near to the equilibrium one, and therefore, these distributions are statistically similar. Note that this torus evolution differs from the result presented in Bannikova et al. (2012), where the initial condition (Keplerian torus) was far from equilibrium and the distribution of clouds was noticeably modified by the evolution of the system. It is seen that the softening length does not influence the resulting distribution of the clouds because the equilibrium state, achieved due to self-gravity, depends on the mass (or the volume density) of the torus.
Figure 2. Histograms of the orbital elements for the initial condition (red) and after 1, 000 average orbital periods for inclination (left), eccentricity (center), and major semi-axis (right), for ϵ = 0.01;0.005 (green, blue). All other parameters correspond to Figure 1.
4. Dynamics of the Clouds in Dusty Torus
Rotation curves allow us to understand the dynamics of clouds in a dusty torus from observations of the megamaser emission. Indeed, the conditions for the formation of water maser emission (λ = 1.35 cm) appear at the inner region of the torus. Rotation curves were obtained only for a few nearby Sy-galaxies, including NGC1068 (Gallimore et al., 1996, 2001, 2004; Greenhill et al., 1996), and demonstrated that the matter in the torus is in sub-keplerian motion. Some models were proposed to explain such a motion and the rotation curve in NGC1068, considering the torus in the approximation of a self-gravitating disk with the mass comparable to the central mass (Huré, 2002; Lodato and Bertin, 2003). This does not agree with the mass of the torus obtained from ALMA observations, (García-Burillo et al., 2016), which corresponds to 0.01MBH.
It is seen from Figure 3 (left) that the motion at the inner boundary of the torus (black and blue points at Figure 3, left) is sub-keplerian. Indeed, the inner boundary of the torus is formed by the clouds that move in orbits with different eccentricities and inclinations (Figure 3, right). These clouds can pass through the apocenter or occupy some arbitrary (temporary) position on these orbits. A spread of the cloud velocities at the inner boundary of the torus can explain the observational data usually assigned to a turbulent motion. The clouds have high values of the velocity in the equatorial plane near the supermassive black hole (orange points at Figure 3, left). These clouds move on orbits with large values of the eccentricity passing exactly through the pericenter; they could ultimately feed the accretion disk.
Figure 3. Distribution of clouds in the torus. All particles are gathered in the meridional plane (ρ, z), where . Left panel: colors mark the deviations of the cloud velocity in the torus (V) from a keplerian velocity for the case of a disk (Vk): (V/Vk − 1). Right panel: colors mark the values of the eccentricity of the cloud orbits.
5. Condition of Obscuration
Here we will determine the probability of obscuration of the central source by clouds of spherical shape with the radius Rcl = ϵ · R, where the major radius of the torus is R = 1. It is convenient to use spherical coordinates, (r, θ, ϕ), and consider a sphere of unit radius r = R. We divide the sphere by the coordinate θ into n uniform parts (bins) whose length on the unit sphere of visibility is h = π/n. The ring length on the sphere for each angle θi is determined by the expression: l(θ) = 2 π ·cosθ. Thus, an area of a ring on the sphere can be written as:
Using the results of N-body simulations, we construct a histogram of the cloud distribution with respect to the angle θ between the equatorial plane and the line-of-sight with a cell size equal to h. The area of the projection of the cloud on the sphere is . Let us suppose that the obscuration area equals the sum of all the cloud areas in the bin:
Then we determine an obscuration coefficient as the ratio of these areas:
It is seen from Figure 4 that the obscuration coefficient depends on the inclination angle of the torus θ which is well fitted by the Gaussian function:
where σ(θ0) is the width of the torus half-opening angle, and the amplitude A characterizes the number of the clouds in the equatorial line-of-sight for the certain value of ϵ. Figure 4 shows that the amplitude A has a linear dependence on the cloud number N and a square dependence of the size ϵ (see also Bannikova et al., 2012). The amplitude values for N = 16, 384 is A = 1.26, while for N = 8, 192 is A = 0.68. So, we can predict that for N = 105 the number of clouds in the equatorial line-of-sight is about 7.7 (for ϵ = 0.01), which is consistent with the estimations of radiative transfer models (Nenkova et al., 2008a,b; Hönig and Kishimoto, 2010).
Figure 4. Dependence of the obscuration coefficient on the angle between the line-of-sight and the equatorial plane. Left panel: values and the corresponding best-fitting curves are obtained from simulations for , N = 16, 384 and ϵ = 0.01, 0.005; N = 8, 192 and ϵ = 0.01. Right panel: each curve corresponds to the different initial half-opening angle of the torus: .
The width of the Gaussian function, σ, depends on the initial half-opening angle of the torus. The right panel of Figure 4 shows the dependence of the obscuration coefficient on the angle between the line-of-sight and the equatorial plane. The fitting of function, kobs(θ), gives the values of σ(45°) = 24.5°, σ(60°) = 32.4°, σ(75°) = 41.9°, and the amplitude A takes the values 1.35, 1.26, and 1.16.
6. Discussion
N-body simulations show that a torus like the one observed in NGC1068 can stay thick, if its clouds initially have a random distribution of the orbital elements and anisotropy in two polar directions. The clouds in such a toroidal structure move in inclined orbits with a spread in eccentricity, and the equilibrium state of the torus corresponds to Gaussian density distribution, which satisfies the obscuration conditions and the observed SED in the IR.
The considered initial distribution of the clouds may be a consequence of the evolutionary processes in AGNs. It may suggest that, at the first stage, the supermassive black hole and the accretion disk are embedded in a quasi spherical distribution of dust (Liu and Zhang, 2011). An example could be the system IRAS16399-0937, a galaxy whose core is immersed in quasi spherically distributed optically thick clouds (Sales et al., 2015). The beginning of the active stage may lead to an increase of the wind energy and to the anisotropy in the cloud distribution. Within the wind cones, the clouds acquire additional impulse against the gravitational forces due to radiation pressure. The dusty clouds located outside of the wind cones are unaffected by the wind and continue to move in inclined and eccentric orbits. These clouds can form the thick toroidal distribution (Figure 1, bottom) which plays the role of an obscuring structure in AGNs.
These simulations do not take into account the effects of dissipation, which will influence the distribution of the clouds and the stability of the torus. The dissipation can be related to the collisions of the clouds and their heating. More frequent collisions will occur near the center of torus cross-section, where density is higher. The heating of the clouds in this region can boost the radiation pressure which may be an additional factor compensating the dissipation effects. The clouds near the supermassive black hole will move to the center due to dissipation and eventually feed the accretion disk. The clouds moving in the inner (and outer) boundary of the torus will rarely collide because their number decreases exponentially. Thus, it can be assumed that the toroidal structure can be conserved for certain values of the dissipation coefficient. Note, that the main mechanisms which were proposed to explain the geometrical thickness of the torus may work in such a dynamical clumpy model.
Author Contributions
All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Acknowledgments
We thank Massimo Capaccioli for the useful discussions, and the referees for helpful comments that improved the manuscript. One of us (EB) likes to thank the Organizing Committee of the conference “Quasars at all cosmic epochs” for financial support to present this work.
References
Antonucci, R. (1993). Unified models for active galactic nuclei and quasars. Annu. Rev. Astron. Astrophys. 31, 473–521. doi: 10.1146/annurev.aa.31.090193.002353
Bannikova, E. Y. (2015). Clouds distribution in obscuring tori of active galactic nuclei. Radio Phys. Radio Astron. 20, 191–204. doi: 10.15407/rpra20.03.191
Bannikova, E. Y., Vakulik, V. G., and Sergeev, A. V. (2012). N-body simulation of a clumpy torus: application to active galactic nuclei. Month. Notices R. Astron. Soc. 424, 820–829. doi: 10.1111/j.1365-2966.2012.21186.x
Chan, C.-H., and Krolik, J. H. (2017). Geometrically thick obscuration by radiation-driven outflow from magnetized tori of active galactic nuclei. Astrophys. J. 843:58. doi: 10.3847/1538-4357/aa76e4
Dorodnitsyn, A., Kallman, T., and Proga, D. (2016). Parsec-scale accretion and winds irradiated by a quasar. Astrophys. J. 819:115. doi: 10.3847/0004-637X/819/2/115
Elitzur, M., and Shlosman, I. (2006). The AGN-obscuring torus: The end of the “Doughnut” paradigm? Astrophys. J. Lett. 648, L101–L104. doi: 10.1086/508158
Gallimore, J. F., Baum, S. A., and O'Dea, C. P. (2004). The parsec-scale radio structure of NGC 1068 and the nature of the nuclear radio source. Astrophys. J. 613, 794–810. doi: 10.1086/423167
Gallimore, J. F., Baum, S. A., O'Dea, C. P., Brinks, E., and Pedlar, A. (1996). H 2O and OH masers as probes of the obscuring torus in NGC 1068. Astrophys. J. 462:740. doi: 10.1086/177187
Gallimore, J. F., Henkel, C., Baum, S. A., Glass, I. S., Claussen, M. J., Prieto, M. A., et al. (2001). The nature of the nuclear H2O masers of NGC 1068: reverberation and evidence for a rotating disk geometry. Astrophys. J. 556, 694–715. doi: 10.1086/321616
García-Burillo, S., Combes, F., Ramos Almeida, C., Usero, A., Krips, M., Alonso-Herrero, A., et al. (2016). ALMA resolves the torus of NGC 1068: Continuum and molecular line emission. Astrophys. J. Lett. 823:L12. doi: 10.3847/2041-8205/823/1/L12
Greenhill, L. J., Gwinn, C. R., Antonucci, R., and Barvainis, R. (1996). VLBI imaging of water maser emission from the nuclear torus of NGC 1068. Astrophys. J. Lett. 472:L21. doi: 10.1086/310346
Hönig, S. F., and Kishimoto, M. (2010). The dusty heart of nearby active galaxies. II. From clumpy torus models to physical properties of dust around AGN. Astron. Astrophys. 523:A27. doi: 10.1051/0004-6361/200912676
Hönig, S. F., and Kishimoto, M. (2017). Dusty winds in active galactic nuclei: reconciling observations with models. Astrophys. J. Lett. 838:L20. doi: 10.3847/2041-8213/aa6838
Huré, J.-M. (2002). Origin of non-keplerian motions of masers in NGC 1068. Astron. Astrophys. 395, L21–L24. doi: 10.1051/0004-6361:20021445
Jaffe, W., Meisenheimer, K., Röttgering, H. J. A., Leinert, C., Richichi, A., Chesneau, O., et al. (2004). The central dusty torus in the active nucleus of NGC 1068. Nature 429, 47–49. doi: 10.1038/nature02531
Krolik, J. H. (2007). AGN obscuring tori supported by infrared radiation pressure. Astrophys. J. 661, 52–59. doi: 10.1086/515432
Liu, Y., and Zhang, S. N. (2011). Dusty torus formation by anisotropic radiative pressure feedback of active galactic nuclei. Astrophys. J. Lett. 728:L44. doi: 10.1088/2041-8205/728/2/L44
Lodato, G., and Bertin, G. (2003). Non-keplerian rotation in the nucleus of NGC1068: evidence for a massive accretion disk? Astron. Astrophys. 398, 517–524. doi: 10.1051/0004-6361:20021672
Nenkova, M., Sirocky, M. M., Ivezić, Ž., and Elitzur, M. (2008a). AGN dusty tori. I. Handling of clumpy media. Astrophys. J. 685, 147–159. doi: 10.1086/590482
Nenkova, M., Sirocky, M. M., Nikutta, R., Ivezić, Ž., and Elitzur, M. (2008b). AGN Dusty Tori. II. Observational implications of clumpiness. Astrophys. J. 685, 160–180. doi: 10.1086/590483
Netzer, H. (2013). The Physics and Evolution of Active Galactic Nuclei. Cambridge, UK: Cambridge University Press.
Proga, D., and Kallman, T. R. (2004). Dynamics of line-driven disk winds in active galactic nuclei. II. effects of disk radiation. Astrophys. J. 616, 688–695. doi: 10.1086/425117
Raban, D., Jaffe, W., Röttgering, H., Meisenheimer, K., and Tristram, K. R. W. (2009). Resolving the obscuring torus in NGC 1068 with the power of infrared interferometry: revealing the inner funnel of dust. Month. Notices R. Astron. Soc. 394, 1325–1337. doi: 10.1111/j.1365-2966.2009.14439.x
Sales, D. A., Robinson, A., Axon, D. J., Gallimore, J., Kharb, P., Curran, R. L., et al. (2015). An embedded active nucleus in the OH megamaser galaxy IRAS16399-0937. Astrophys. J. 799:25. doi: 10.1088/0004-637X/799/1/25
Schartmann, M., Burkert, A., Krause, M., Camenzind, M., Meisenheimer, K., and Davies, R. I. (2010). Gas dynamics of the central few parsec region of NGC 1068 fuelled by the evolving nuclear star cluster. Month. Notices R. Astron. Soc. 403, 1801–1811. doi: 10.1111/j.1365-2966.2010.16250.x
Stalevski, M., Fritz, J., Baes, M., Nakos, T., and Popović, L. Č. (2012). 3D radiative transfer modelling of the dusty tori around active galactic nuclei as a clumpy two-phase medium. Month. Notices R. Astron. Soc. 420, 2756–2772. doi: 10.1111/j.1365-2966.2011.19775.x
Tristram, K. R. W., Burtscher, L., Jaffe, W., Meisenheimer, K., Hönig, S. F., Kishimoto, M., et al. (2014). The dusty torus in the Circinus galaxy: a dense disk and the torus funnel. Astron. Astrophys. 563:A82. doi: 10.1051/0004-6361/201322698
Tristram, K. R. W., Meisenheimer, K., Jaffe, W., Schartmann, M., Rix, H.-W., Leinert, C., et al. (2007). Resolving the complex structure of the dust torus in the active nucleus of the Circinus galaxy. Astron. Astrophys. 474, 837–850. doi: 10.1051/0004-6361:20078369
Urry, C. M., and Padovani, P. (1995). Unified schemes for radio-loud active galactic nuclei. Publ. Astron. Soc. Pac. 107:803. doi: 10.1086/133630
Wada, K. (2012). Radiation-driven fountain and origin of torus around active galactic nuclei. Astrophys. J. 758:66. doi: 10.1088/0004-637X/758/1/66
Keywords: active galactic nuclei (AGN), quasars, Seyfert galaxies, N-body problems, gravitation
Citation: Bannikova EY and Sergeyev AV (2017) Dynamics and Formation of Obscuring Tori in AGNs. Front. Astron. Space Sci. 4:60. doi: 10.3389/fspas.2017.00060
Received: 30 August 2017; Accepted: 28 November 2017;
Published: 12 December 2017.
Edited by:
Mauro D'Onofrio, Università degli Studi di Padova, ItalyReviewed by:
Omaira Gonzàlez Martín, Instituto de Radioastronomía y Astrofísica, MexicoAlessandro Ballone, Osservatorio Astronomico di Padova (INAF), Italy
Copyright © 2017 Bannikova and Sergeyev. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Elena Yu. Bannikova, bannikova@astron.kharkov.ua