RADIATION-DRIVEN FOUNTAIN AND ORIGIN OF TORUS AROUND ACTIVE GALACTIC NUCLEI

Based on our three-dimensional, multi-phase hydrodynamic model (Wada et al. 2009), we included the radiative feedback from the central source, i.e., radiation pressure on the dusty gas and the X-ray heating of cold, warm, and hot ionized gas. In contrast to previous analytical and numerical studies of radiation-dominated tori (e.g., Pier & Krolik 1993; Krolik 2007; Shi & Krolik 2008; Diamond-Stanic & Rieke 2012; Dorodnitsyn et al. 2011, 2012), we assume neither dynamical equilibrium nor geometrical symmetry. Here, we solve fully three-dimensional hydrodynamic equations, accounting for radiative feedback processes from the central source using a ray-tracing method. In a manner similar to the approach by Wada & Norman (2002), we account for self-gravity of the gas, radiative cooling for 20 K <T_gas < 10⁸ K, uniform UV radiation for photoelectric heating, and H₂ formation/destruction (Wada et al. 2009); however, we do not include the effect of supernova feedback in this study, in order to clarify the effect of the radiation feedback only. We plan to investigate the effects of supernova in the circumnuclear starburst in a future study.

We use a numerical code based on Eulerian hydrodynamics with a uniform grid (Wada & Norman 2001; Wada 2001; Wada et al. 2009). We solve the following equations numerically to simulate the three-dimensional evolution of a rotating gas disk in a fixed spherical gravitational potential under the influence of radiation from the central source.

$$\begin{eqnarray} {\partial \rho}/{\partial t} + \nabla \cdot (\rho \mbox{\boldmath ${v}$}) &=& 0, \end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray} {\partial (\rho \mbox{\boldmath ${v}$})}/{\partial t} + (\mbox{\boldmath ${v}$} \cdot \nabla)\mbox{\boldmath ${v}$}+{\nabla p} &=& - \rho \big(\nabla \Phi +\mbox{\boldmath ${f}$}_{{\rm rad}}^{r} \big), \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} && {\partial (\rho E)}/{\partial t} + \nabla \cdot [(\rho E+p)\mbox{\boldmath ${v}$}] -\rho \mbox{\boldmath ${v}$}\cdot \nabla \Phi = \rho \Gamma _{\rm UV}(G_0) \nonumber\\ &&\quad + \rho \Gamma _{\rm X} - \rho ^{2} \Lambda (T_{{\rm gas}}, f_{\rm H_2}, G_0), \end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray} \nabla ^2 \Phi _{\rm sg}&=& 4 \pi G \rho, \end{eqnarray} \tag{ 4 }$$

where $\Phi (\mbox{\boldmath ${x}$}) \equiv \Phi _{\rm ext}(r) + \Phi _{\rm BH}(r)+\Phi _{\rm sg}(\mbox{\boldmath ${x}$})$; $\rho,p,\, {\rm and} \; \mbox{\boldmath ${v}$}$ denote the density, pressure, and velocity of the gas, and the specific total energy $E \equiv |\mbox{\boldmath ${v}$}|^2/2+ p/(\gamma -1)\rho$, with γ = 5/3. We assume a time-independent external potential Φ_ext(r) ≡ −(27/4)^1/2[v²₁/(r² + a²₁)^1/2 + v₂²/(r² + a²₂)^1/2], where a₁ = 100 pc, a₂ = 2.5 kpc, v₁ = 147 km s⁻¹, v₂ = 147 km s⁻¹, and Φ_BH(r) ≡ −GM_BH/(r² + b²)^1/2, where M_BH = 1.3 × 10⁷ M_☉ (see Figure 1 in Wada et al. 2009 for the rotation curve). The potential caused by the BH is smoothed within r ∼ b = 4δ, where δ denotes the minimum grid size (=0.25 pc), in order to avoid too small time steps around the BH. In the central grid cells at r < 2δ, physical quantities remain constant.

We solve the hydrodynamic part of the basic equations using the advection upstream splitting method (Liou & Steffen 1993). We use 256³ grid points. The uniform Cartesian grid covers a 64³ pc ³ region around the galactic center (i.e., the spatial resolution is 0.25 pc). The Poisson equation, Equation (4), is solved to calculate the self-gravity of the gas using the fast Fourier transform and the convolution method with 512³ grid points along with a periodic Green's function. We solve the non-equilibrium chemistry of hydrogen molecules along with the hydrodynamic equations (Wada et al. 2009).

We consider the radial component of the radiation pressure:

$$\begin{eqnarray} &&\mbox{\boldmath ${f}$}_{{\rm rad}}^{r} = \int \frac{\chi _{\nu}\mbox{\boldmath ${F}$}^{r}_{\nu}}{c} d\nu, \end{eqnarray} \tag{ 5 }$$

where χ_ν denotes the total mass extinction coefficient due to dust absorption and Thomson scattering, i.e., χ_ν ≡ χ_{dust, ν} + χ_T. The radial component of the flux at the radius r, $\mbox{\boldmath ${F}$}^{r}_{\nu}$ is

$$\begin{eqnarray} &&\mbox{\boldmath ${F}$}^{r}_{\nu} \equiv \frac{L_{\nu}(\theta) e^{-\tau _{\nu}}}{4\pi r^{2}}\mbox{\boldmath ${e}$}_{r}, \end{eqnarray} \tag{ 6 }$$

where τ_ν = ∫χ_νρds.

The only explicit radiation source used here is an accretion disk whose size is five orders of magnitude smaller than the grid size in the present calculations. Therefore, we assumed that radiation is emitted from a point source. However, the radiation flux originating from the accretion disk is not necessarily spherically symmetric (e.g., Netzer 1987). In our study, we simply assume emission from a thin accretion disk without considering the limb-darkening effect, i.e., L_ν ≡ 2L_AGN|cos θ|, where θ denotes the angle from the rotational axis (z-axis; see also the discussion in Section 4).

In the simulations, we considered the radial component of the flux $\mbox{\boldmath ${f}$}_{{\rm rad}}^{r}$ only for the radiative heating and pressure. This approximation is not necessarily correct, especially with respect to optically thick dusty gas around the central few parsecs; however, beyond this region, the acceleration of the gas due to strong X-rays can be expected to be radial. Even for the central region spanning a few parsecs, Roth et al. (2012) showed that the gas acceleration is nearly radial using three-dimensional Monte Carlo radiative transfer calculations. It is noteworthy that since the hydrodynamics are calculated in a fully three-dimensional manner, the radiative acceleration may contribute to gas dynamics along any direction. Therefore, the resultant density and velocity structures could be very different from the symmetric ones (cf. Ohsuga & Umemura 2001). The optical depth τ_ν is calculated at every time step along a ray from the central source at each grid point, i.e., 256³ rays are used in the computational box.

In the energy equation (Equation (3)), we use a cooling function based on a radiative transfer model of photodissociation regions (Meijerink & Spaans 2005), denoted as $\Lambda (T_g, f_{\rm H_2}, G_0)$, assuming solar metallicity, which is a function of the molecular gas fraction f_H2 and the intensity of far-ultraviolet (FUV), G₀, (Wada et al. 2009). The radiative cooling rate below ≃ 10⁴ K is self-consistently modified depending on $f_{\rm H_2}$, which is determined for each grid cell. We include the effect of heating due to the photoelectric effect, Γ_UV. We assume a uniform FUV radiation field, with the Habing unit G₀ = 1000, originating in the star-forming regions in the central several tens of parsecs. In reality, both the density field and the distribution of massive stars should not be uniform, and G₀ should have a large dispersion with some radial dependence. However, the FUV radiation affects the chemistry of the ISM (mostly distribution of molecular clouds), although this does not significantly alter the dynamics of the circumnuclear ISM (Wada et al. 2009). The effects of the non-uniform radiation field with multiple sources on the gas dynamics is an important problem, but it is beyond the scope of this paper.

The temperature of the interstellar dust, T_dust, at a given position irradiated by the central source is calculated by assuming thermal equilibrium (e.g., Nenkova et al. 2008), and if T_dust > 1500 K, then we assume that the dust is sublimate. Here, we assume that the ISM in the central tens parsecs is optically thin with respect to the re-emission from the host dust. The SED of the AGN and the dust absorption cross section are taken from Laor & Draine (1993). We use a standard galactic dust model, which is the same as that used by Schartmann et al. (2010), with eight logarithmic frequency bins between 0.001 μm and 100 μm.

The effects of hard X-ray on the thermal and chemical structures of the ISM are significant in the region spanning subparsec distances to tens of parsecs (Pérez-Beaupuits et al. 2011). Furthermore, X-ray heating may be essential for galactic-scale feedback (Hambrick et al. 2011). Here, we consider heating by the X-ray, i.e., interactions between the high-energy secondary electrons and thermal electrons (Coulomb heating), and H₂ ionizing heating in the warm and cold gas (Maloney et al. 1996; Meijerink & Spaans 2005). The Coulomb heating rate Γ_{X, c} for the gas with the number density n is given by

$$\begin{eqnarray} &&\Gamma _{\rm X,c} = \eta n H_{\rm X}, \end{eqnarray} \tag{ 7 }$$

where η denotes the heating efficiency (Dalgarno et al. 1999; Meijerink & Spaans 2005), and the X-ray energy deposition rate H_X = 3.8 × 10⁻²⁵ξ erg s⁻¹ and the ionizing parameter $\xi = 4\pi F_{\rm X}/n = L_{\rm X}e^{-\int \tau _{\nu} d\nu}/n r^{2}$, with L_X = 0.1L_AGN. For optically thin hot gas with T ≳ 10⁴ K, the effects of Compton heating and X-ray photoionization heating are considered, and the net heating rate (Blondin 1994) is approximately

$$\begin{eqnarray} \Gamma _{\rm X,h} &=& 8.9\times 10^{-36} \xi (T_{\rm X} - 4T) \nonumber\\ &&+1.5 \times 10^{-21} \xi ^{1/4} T^{-1/2}(1-T/T_{\rm X}) \;\; {\rm erg}\, {\rm s}^{-1}\, {\rm cm}^{3},\quad \end{eqnarray} \tag{ 8 }$$

with the characteristic temperature of the X-ray radiation, T_X = 10⁸ K.

In order to prepare a quasi-steady initial condition without radiation feedback, we first run a model of an axisymmetric and rotationally supported thin disk with a uniform density profile (thickness of 2.5 pc) and a total gas mass of M_g = 6.68 × 10⁶ M_☉. Random density fluctuations, which are less than 1% of the unperturbed values, are added to the initial gas disk. After t ∼ 1.7 Myr, the thin disk settles in a flared shape with the total mass 6.64 × 10⁶ M_☉. Since we allow outflows from the computational boundaries, the total gas mass decreases during the evolution.

A free parameter in the present simulations is the luminosity of the AGN. Here, we explore models with two different Eddington ratios L_AGN/L_E = 0.1 and 0.01, where the Eddington luminosity L_E ≡ 4πGcm_pM_BH/σ_T. In each model, L_AGN stays constant during calculations, i.e., ∼3 Myr.¹

Figure 1 shows the onset of biconical outflows, cavities, and thick disk after the radiation from the central source is initiated. Radiation-driven outflows start forming in the inner region, and subsequently propagate outward. At t = 1.986 Myr, "back-flows" toward the disk plane appear outside the biconical outflows. These accretion flows interact with the disk, thereby leading to the formation of a geometrically thick disk that has a clear boundary with ρ ∼ 10² M_☉ pc⁻³ (Figure 1).

Zoom In Zoom Out Reset image size

The interactions of vertical flows with opposite directions are more clearly seen in Figure 2 (t = 2.146 Myr), which shows the distribution of v_z. There are several contact surfaces between upward and downward flows. As a result, turbulent motions that support the thick geometrical structure of the disk are generated. The complicated distribution of flows in the thick disk is also seen in the velocity field at the equatorial plane (left panel of Figure 2).

Zoom In Zoom Out Reset image size

The downward flows can be naturally expected due to the balance between radiation pressure and gravity. The ratio between the radiation pressure for the dusty gas and the gravity due to the central BH is

$$\begin{eqnarray} &&\frac{f_{{\rm rad}}}{f_{{\rm gravi}}} = \frac{L_{{\rm AGN}} \chi _{d} |\cos \theta |}{4\pi G c M_{{\rm BH}}}e^{-\tau _{d}}, \end{eqnarray} \tag{ 9 }$$

where τ_d denotes the optical depth for the dusty gas at a given point. Suppose that the gas density is a function of z, e.g., $n(z) = N_{{\rm col}}/\sqrt{2\pi}h e^{-z^{2}/2h^{2}}$. Figure 3 shows how the radiation-pressure-dominated region and gravity-dominated region are distributed around the AGN for the given gas density. Here, we assume the dust-to-gas mass ratio of 0.03 (Ohsuga & Umemura 2001), L_AGN = 10⁴⁴ erg s⁻¹, dust extinction coefficient χ_d = 10⁵ cm² g⁻¹, and the gas density at z = 0 as n₀ = 1000 cm⁻³. The four curves represent the condition f_rad/f_gravi = 1, with the vertical extent of the gas h = 4, 6, 8, and 10 pc, for which radiation pressure dominates in the inner side of the curve, therefore outflows are expected. On the other hand, the radiation pressure cannot prevent gas accretion outside the critical lines. In fact, accretion flows appear in the gravity-dominated domain as seen in Figures 1 and 2. Moreover, since the critical lines dividing the two domains are not purely radial (Figure 3), outflows with large angles from the z-axis may collide with the boundary of the cavity. This pushes the boundaries of the outflow outward, and the gases in the outflows change their directions and fall toward the center as indicated by the velocity vectors in Figure 2. This also causes the biconical cavities to become non-steady, and the vertical circulation of the flow, or the "fountain" of gas, circulates around the AGN on a scale of 10 pc. The kinetic energy of this circulating flow is partially used to generate the random motion in the disk and maintain its thickness. In other words, both radiation energy and gravity are the cause of origin of the internal turbulent motion in the disk.²

Zoom In Zoom Out Reset image size

Figure 4 shows that the density structures in a quasi-steady states differ depending on the luminosity of the central source. Here, "quasi-steady" means that the global morphology has reached an approximately steady state; that is, the phase diagram of the gas (Figure 6) does not significantly change after this point. The accretion rates to the central region are also approximately constant. This is achieved roughly after t = 3.0–3.5 Myr (Figure 9). Note that since the duration of the present simulation is much shorter than the accretion timescale (∼10⁹ Myr), the "quasi-steady" does not necessary mean the long-term stable.

Zoom In Zoom Out Reset image size

In the more luminous model (L_AGN/L_Edd = 0.1, upper panels), inhomogeneous outflows with ρ ∼ 0.1–1 M_☉ pc⁻³ and velocity of ∼100 km s⁻¹ are formed, and diffuse gas (ρ ∼ 0.1 M_☉ pc⁻³) extends for over 20 pc above the equatorial plane. On the other hand, diffuse and relatively stable outflows with ρ ∼ 10⁻² M_☉ pc⁻³ or less are formed in the less luminous model with L_AGN/L_Edd = 0.01. The gas temperature distributions of the two models (Figure 5) show that dense spiral arms and clumps in the equatorial plane exhibit a temperature of less than 1000 K, in contrast to the "warm" (∼10⁴–10⁵ K) torus. In the less luminous AGN model, an increased number of hot gas (≳ 10⁶ K) regions appear in the biconical outflow. This is because the radiative cooling is less effective due to reduced gas density in the region, and as a consequence, the X-ray is not attenuated.

Zoom In Zoom Out Reset image size

Figure 6 shows the effects of the radiative feedback on the phases of the gas. The dominant phase is the one around T_gas ∼ 8000 K. The dense and cold gas (ρ > 10² M_☉ pc⁻³ and T_gas < 1000 K) is not affected by the radiative feedback. The X-ray heating produces warm (T_gas < 1000 K) gas in the outer regions (i.e., high latitude) of the thick disk, as is expected from the X-ray dominated region model (Pérez-Beaupuits et al. 2011). Shock-heated gas with a temperature of several 10⁴ K also appears in the same density region.

Zoom In Zoom Out Reset image size

The structure of the torus does not show perfect symmetry for the rotational axis and the equatorial plane. Figure 7 shows four different vertical slices of density around the rotational axis of the model shown in the upper panel of Figure 4. The diffuse part of the torus (ρ ∼ 10²) is slightly tilted, thereby indicating that the outlying region of the disk is not completely settled. The shell-like outflows also show asymmetric distributions. It is interesting to note that these structures are caused by interactions between the radiation from the central symmetric source and the surrounding inhomogeneous gas, thereby demonstrating the importance of the three-dimensional treatment of radiation hydrodynamics in understanding the structure of AGNs.

Zoom In Zoom Out Reset image size

Obscuration due to the vertically extended gas in the central region over tens of parsecs and how it differs between models are shown in Figure 8, which shows the plot of the column density toward the central BH as a function of the viewing angle. If there is no radiation from the central source, then the opening angle for N ⩽ 10²³ cm⁻² is approximately 140°. The radiation feedback causes the thickness of the disk to increase, and the opening angles are approximately 60° and 100° for L_AGN/L_E = 0.1 and 0.01, respectively. It is noteworthy that the column density does not only depend on the viewing angle, but also on the azimuthal angles. The column density varies by a factor of 10 or more for a given viewing angle, thereby indicating the internal inhomogeneous structure of the disk. The structure of the disk is similar to that generated by SNe (Wada & Norman 2002; Wada et al. 2009). However, the present simulations imply that the geometrically and optically thick torus can be naturally formed due to gravitational energy in the central part with the aid of the radiation feedback in the case of the AGNs with Eddington ratios of 0.01-0.1.

Zoom In Zoom Out Reset image size

The radiation pressure caused by the central source usually has a "negative" effect on the mass accretion toward the center; however, our results indicate that mass accretion is not terminated by radiative feedback. Figure 9 shows the time evolution of the mass of the two models in four radial bins. The lower luminosity model (L_AGN/L_Edd = 0.01) shows a continuous net mass accretion in the central 8 pc region, and there is almost no net inflow beyond r ∼ 4 pc. The mass accretion in the central part is caused by gravitational torque due to the non-axisymmetric features near the equatorial plane of the central part of the disk (Figure 1).³

Zoom In Zoom Out Reset image size

The accretion rate to the central parsec increases with decreasing AGN luminosity: $\dot{M} = 1.7\times 10^{-3}\ M_{\odot}\, {\rm yr}^{-1}$ and $\dot{M} = 1.3\times 10^{-4}\ M_{\odot}\, {\rm yr}^{-1}$ for L_AGN/L_E = 0.01 and 0.1, respectively. For L_AGN = 0.1L_E, the net accretion rate is one order of magnitude smaller than the accretion rate required for constant AGN luminosity; the AGN luminosity is assumed to be constant during the calculations provided that the energy conversion efficiency is 0.1. On the other hand, the observed accretion rate in the low-luminosity AGN model is comparable to that required for the luminosity of the source. These results may suggest that luminous AGNs are intermittent, and their luminosity is not instantaneously determined by mass inflow toward r ∼ 1 pc.