Dynamics and formation of obscuring tori in AGNs

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.


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(Tristram et al., , 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: M torus = 10 5 M (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;Chan and Krolik, 2017;Dorodnitsyn et al., 2016), 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. Section 3 and Section 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.

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, θ wind = 30 • , are shown in Fig.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.
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, M cl = M torus /(N − 1): where r i = (x, y, z) is the vector of a cloud from the central mass normalized to the major radius of the torus (R); a i is the vector of the acceleration of the i-th cloud acquired from all of the clouds of the torus M torus and from the central mass M BH . 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 R cl = R. We choose the physical parameters corresponding to the case of NGC1068: M BH = 10 7 M and M torus = 10 5 M = 0.01M BH . Since the method of parallel calculations is used, the number of particle in N-body problem must be N = 2 n . In presented simulations we adopt n = 14 (or n = 13), thus the number of clouds comes N cl = N − 1 = 16, 383 (or N cl = 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, |E − E 0 | = 5 · 10 −6 , where E 0 is the initial value of the total energy. We use the unity system: M BH = G = R = 1. 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 ρ = ± x 2 + y 2 (right). Initial conditions: 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 Fig. 1. (One orbital period corresponds T = 30, 000 years for M BH = 10 7 M 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 (N 1 , N 2 ), the same velocity dispersion of the clouds in the torus will be reached for the second system in the time interval t 2 = t 1 N 2 /N 1 , where t 1 is the time interval for the first system (Bannikova et al., 2012).

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 (Fig.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. 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 Fig.1.
Knowing from N-body simulations the coordinates and velocity components, we can calculate the orbital elements of each cloud. Fig.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 (Fig. 2 left, right) do not essentially differ from the initial ones, while there is a little difference in the distribution for eccentricity (Fig. 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.

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(Gallimore et al., , 2001(Gallimore et al., , 2004Greenhill 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, M torus = 10 5 M (García- Burillo et al., 2016), which corresponds to 0.01M BH . It is seen from Fig.3 (left) that the motion at the inner boundary of the torus (black and blue points at Fig.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 (Fig.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 Fig.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.

CONDITION OF OBSCURATION
Here we will determine the probability of obscuration of the central source by clouds of spherical shape with the radius R cl = · 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 S cl i = π( /r i ) 2 . 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 Fig. 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 . Fig.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 = 10 5 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).
The width of the Gaussian function, σ, depends on the initial half-opening angle of the torus. The right panel of Fig.4 shows the dependence of the obscuration coefficient on the angle between the lineof-sight and the equatorial plane. The fitting of function, k obs (θ), 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.

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 of the clouds 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 (Fig. 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.