3D Radiation Hydrodynamics of a Dynamical Torus

We use the public version of the N-body+hydrodynamics code GIZMO (Hopkins 2015) in pressure smoothed particle hydrodynamics (P-SPH) mode. GIZMO's P-SPH algorithm is a "modern" SPH method as opposed to the "traditional" T-SPH mode. P-SPH is a pressure–entropy scheme, incorporating artificial conductivity and switches for artificial viscosity. GIZMO is a descendant of GADGET (specifically GADGET-3; see Springel 2005), an N-body+SPH code in widespread use, and GIZMO shares many of the same features as GADGET, including file formats. GIZMO primarily differs from GADGET by including a new hydrodynamics scheme that differs significantly from SPH, but in our initial tests we found that this new algorithm sometimes developed extremely small time steps that essentially halted the simulation. It was not clear to us whether the cause of this issue was a bug in the particular public version of GIZMO we downloaded (a bug that may have been fixed in a more recent release), a genuine issue with the algorithm, or a problem in our particular model and its implementation, but recent results show that many of the results of cosmological simulations do not strongly depend on the hydrodynamics scheme (Hopkins et al. 2018). Hence, rather than attempting to disentangle this problem, we chose to run GIZMO in the simpler P-SPH mode, where we did not encounter such vanishingly small time steps.

Our simulations include self-gravity, which can produce gravitational instabilities, as we discuss in Section 3.1. Additionally, we include a static background potential, which we discuss in Section 2.6.

We summarize the key points of our method in the subsections below.

The previous simulations of the AGN torus region mentioned in Section 1 were performed with an Eulerian formulation on a fixed grid in two or three dimensions. We differ from the standard approach by using SPH, a Lagrangian formulation of hydrodynamics. Eulerian methods are limited by length resolution and tend to resolve high-density fluid poorly, while Lagrangian methods are limited by mass resolution and tend to resolve low-density fluid poorly. Eulerian methods suffer from an inherent numerical diffusivity (in addition to an explicit numerical viscosity), which causes fluid to spread out, especially at high velocities or when the velocity field is not aligned with the grid, producing errors in high-velocity outflows. On the other hand, SPH methods include an explicit numerical viscosity term to capture shocks, which dissipates energy and causes fluid to collapse, as well as suppressing fluid instabilities. There are various techniques to alleviate the problems with either method, such as adaptive mesh refinement and particle splitting, but these remain only partial solutions. Essentially, SPH and grid methods have complementary strengths and weaknesses, and by choosing to use a Lagrangian method, we are exploring the problem from an almost orthogonal perspective to existing studies. In this particular context, our model will not effectively resolve low-density winds in the ionization cone of the AGN. However, it will efficiently resolve the denser gas that contributes the most to the opacity and IR luminosity of the torus region.

Traditional SPH methods have a number of well-known problems with sharp density gradients, producing a spurious surface tension effect, and suppressing turbulence and the Kelvin–Helmholtz and Rayleigh–Taylor instabilities (Agertz et al. 2007). Several modifications to traditional SPH have been proposed to alleviate these effects (see Hopkins 2013, and references therein). The pressure–entropy formulation (Hopkins 2013; Saitoh & Makino 2013) explicitly calculates pressure from the smoothing kernel rather than from particle densities. This removes the spurious forces that result from sharp density contrasts, provided that the physical pressure gradient is smooth. An artificial viscosity is still required to capture shocks, where the pressure gradient is very steep. This method is inherently adiabatic in the sense that entropy is manifestly conserved in the equations of motion, while energy conservation is only subject to very small errors. Note that, with the appropriate choice of definitions, the pressure–entropy formulation is equivalent to a "pressure–energy" formulation.

The basic equations of motion to be solved are the compressible Euler equations. Mass is conserved trivially in SPH, and so we need only solve the momentum and energy conservation equations:

$$\begin{eqnarray}&&\displaystyle \frac{D{\boldsymbol{v}}}{{Dt}}=-\displaystyle \frac{{\rm{\nabla }}P}{\rho }+{\boldsymbol{g}}+{{\boldsymbol{a}}}_{{\boldsymbol{r}}},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{Du}}{{Dt}}=-\displaystyle \frac{P}{\rho }{\rm{\nabla }}\cdot {\boldsymbol{v}}+H,\end{eqnarray} \tag{ 2 }$$

where D/Dt is the convective derivative; P, ρ, u, and ${\boldsymbol{v}}$ are the pressure, density, internal energy, and velocity of the fluid, respectively; ${\boldsymbol{g}}$ is the acceleration from gravity; ${{\boldsymbol{a}}}_{{\boldsymbol{r}}}$ is the radiative acceleration; and H is the combination of all heating and cooling effects. The base equations without the source terms (${\boldsymbol{g}},{{\boldsymbol{a}}}_{{\boldsymbol{r}}}$, H) are discretized in the code following the pressure–entropy formulation in Hopkins (2013). The source terms are then added to the purely hydrodynamic rates, which are integrated with a leapfrog scheme. Heating and cooling rates can vary rapidly with temperature, and so to avoid short system time steps, we integrate heating and cooling over shorter substeps and apply the time-averaged rate to Du/Dt.

We assume an ideal equation of state, i.e.,

$$\begin{eqnarray}&&P=(\gamma -1)u\rho ,\end{eqnarray} \tag{ 3 }$$

where the adiabatic constant γ = 5/3. Temperature is only calculated for CLOUDY table lookup and for visualization, and it is calculated through

$$\begin{eqnarray}&&T=(\gamma -1)\displaystyle \frac{\mu {m}_{p}}{{k}_{{\rm{B}}}}u,\end{eqnarray} \tag{ 4 }$$

where k_B is the Boltzmann constant, m_p is the proton mass, and μ ≈ 0.122 is the mean molecular mass, i.e., ρ/n_H.

We have implemented a ray-tracing algorithm within GIZMO that calculates the optical depth between the central AGN and each gas particle. Ray-tracing, even from a single point, can be computationally expensive—a naive algorithm scales with the square of the number of particles, N², as there are N particles receiving radiation and N particles to potentially contribute extinction. Our algorithm uses an oct-tree structure to efficiently identify which particles intersect each ray and uses this to calculate the optical depth between the particle and the AGN. We use the algorithm of Revelles et al. (2000) for efficient oct-tree traversal.

We have performed a suite of CLOUDY models (Ferland et al. 2013, 2017) to tabulate the internal optical depth, heating, cooling, and radiation pressure of dusty gas as a function of the gas temperature, the gas density, the optical depth from the gas particle to the AGN, and the intensity of incident AGN radiation. The CLOUDY models also include a radiation field from bulge stars (see Section 2.4). For each particle, we identify the particles that intersect a ray from the AGN to the particle's center. For each intersecting particle, we calculate its contribution to the optical depth from its tabulated opacity and the intersecting column density, determined from the impact parameter through a cubic spline kernel of variable smoothing length.

That is, most of the details of radiative transfer are precalculated with CLOUDY, and only the optical depths from the AGN to each particle are calculated in the code. Essentially, the radiative transfer equation is reduced to simply

$$\begin{eqnarray}&&{\zeta }_{i}={\zeta }_{\mathrm{tab}}({T}_{i},{\rho }_{i},{F}_{i},{\tau }_{i}),\end{eqnarray} \tag{ 5 }$$

where ζ_i represents the radiative acceleration, heating/cooling rate, or opacity of particle i, and ${\zeta }_{\mathrm{tab}}({T}_{i},{\rho }_{i},{F}_{i},{\tau }_{i})$ is that value interpolated from the precalculated CLOUDY tables as a function of T_i, the particle temperature, ρ_i, the particle density, F_i, the unextinguished AGN flux, and τ_i, the optical depth between the AGN and the particle. F_i is defined below in Section 2.3.2. τ_i is calculated as

$$\begin{eqnarray}&&{\tau }_{i}=\displaystyle \sum _{j}{\kappa }_{j}{{\rm{\Sigma }}}_{j}({{\boldsymbol{r}}}_{i},{{\boldsymbol{r}}}_{\mathrm{AGN}}),\end{eqnarray} \tag{ 6 }$$

where the sum is over all particles that are close enough to perhaps intersect the ray between the AGN and particle i (as determined by the oct-tree algorithm), κ_j is the opacity (in area/mass units) of particle j as determined from the table in the previous step, and ${{\rm{\Sigma }}}_{j}({{\boldsymbol{r}}}_{i},{{\boldsymbol{r}}}_{\mathrm{AGN}})$ is the column density contributed by particle j to a ray from the AGN to the center of particle i, as calculated by the ray-tracing method.

2.3.1. Resolution

Most of the momentum deposition from the radiation pressure of the AGN is applied to the optically thin (τ ≲ 1) "skin" of the gas. If this is not resolved—that is, if the gas particles that receive AGN radiation are optically thick—then the wind will be suppressed. Instead of force being applied to an optically thin skin, which then peels away as a wind, this momentum would be distributed over the optical depth of the particle, diluting the acceleration of the gas. Furthermore, an optically thick particle can no longer be reasonably assumed to have a single value of temperature, opacity, or heating and cooling rates.

In the Lagrangian formulation, the central column density of an SPH particle is determined by its mass and smoothing length, where the smoothing length h is adaptive to ensure that each particle has enough neighbors. For a particle of fixed mass m_p, the column density is proportional to ${m}_{p}{h}^{-2}$. Hence, at any mass resolution, for a small enough h, the column density for a particle is high enough that the particle is optically thick. The adaptive smoothing length is effectively a density criterion, which roughly ensures that $h\propto {\rho }^{-1/3}/{m}_{p}$. If we define the threshold for an optically thick column density to be Σ_T, the cutoff density is proportional to ${{\rm{\Sigma }}}_{T}^{3/2}{m}_{p}^{-1/2}$, beyond which particles become optically thick.

Hence, while out-of-the-box SPH with adaptive smoothing and softening lengths can capture the dynamics of dense gas, such as that produced by gravitational collapse (although this becomes computationally expensive as densities become large), radiative processes will not be correctly treated. Additionally, CLOUDY does not always produce converged results at high densities.

We have taken several steps to alleviate these problems.

• 1.  
  We only calculate CLOUDY tables up to n = 10⁹ cm⁻³, and we extrapolate for higher n values, by assuming that the heating rate remains constant but that the cooling rate per unit mass is proportional to density (that is, the volumetric cooling rate is proportional to n²). This dense gas mostly occurs in dense cores, and the precise form of the extrapolation does not affect the nature of the wind, which is the focus of this study. To avoid overaccelerating distant gas, we also extrapolate the radiation pressure for fluxes below F = 10^1.25 erg s⁻¹ cm⁻² by assuming that the radiation pressure is proportional to the flux.
• 2.  
  We set the star formation threshold to n = 10¹⁰ cm⁻³, providing a sink for the very high density gas. Our star formation algorithm is explained more fully in Section 2.5.
• 3.  
  We generate a new set of tables that account for self-shielding. The optical depth of a particle depends on its opacity (as determined from the table), its mass, and its softening length (and the softening kernel). Particle mass and kernel shape are constant in the simulation, opacity is already determined from the table, and the smoothing length is very close to a monotonic function of density. Hence, this new table does not introduce any new independent variables and is only accurate for particles of a single fixed mass. We calculate the new table by dividing each particle into 2000 slices of equal column density (but not equal mass). These slices are thin enough that each is optically thin, and we can use the CLOUDY tables to determine the dependent properties (temperature, density, heating, cooling) of each slice, because they do not vary much over each slice. We then perform a mass-weighted average over the slices to determine the new table values for gas particles of that density.
• 4.  
  We initialize our disk to have a fairly low surface density, to reduce the number of high-density particles where self-shielding dominates.

To demonstrate the importance of self-shielding on individual particles, the radiative acceleration of particles of temperature T = 100 K from an AGN flux of 10⁶ erg⁻¹ s⁻¹ cm⁻² of various densities and particle masses are plotted in Figure 1, along with the typical smoothing lengths of particles of this mass resolution with a density of n_H = 10⁴ cm⁻³. At large particle masses and high densities, self-shielding means that the radiation is absorbed in a thin skin and then spread over the mass of the entire particle, suppressing the acceleration. If self-shielding was ignored and the entire particle received the acceleration that is received at its face, then the radiative acceleration would be massively overestimated instead. We choose a particle mass of 10⁻⁴ M_⊙, but note that even with this fine resolution, acceleration is significantly suppressed for densities greater than 10⁵ cm⁻³. As a basis of comparison with Eulerian simulations, we calculate the typical smoothing length of particles with a mass of 10⁻⁴ M_⊙ and a density of n_H = 10⁴ cm⁻³ to be about 10^−2.5 pc. To resolve the radiative acceleration of dense dusty gas, an Eulerian simulation must have a resolution at least this fine.

Zoom In Zoom Out Reset image size

As an illustration of this, we performed test simulations with resolutions from 10⁻⁵ to 10⁻¹ M_⊙. To reduce computational expense in these tests, we kept the total number of particles constant and the initial surface density constant, and so the outer edge of the initial disk is larger with increasing particle mass. Plots of the density distribution of the gas are in Figure 2. The biconical wind is only formed at the highest resolutions. This can be seen in the phase plots in Figure 3, where the hot wind phase is either not present or poorly resolved at m_p ≥ 0.001 M_⊙. High-density gas is also not resolved at higher resolutions, and much of the gas overcools. The outflow rates and mean wind inclinations are also dependent on resolution (Figure 4). We still see some resolution dependence at our finest test resolution, and so our simulations are not entirely converged with respect to resolution, but we find that this dependence weakens beyond a threshold resolution of ∼10⁻⁴ M_⊙. At this resolution, some of the gas still reaches densities high enough that their radiative acceleration would still be theoretically suppressed. However, the resolution dependence becomes weaker, as much of this very dense gas is shielded from the AGN by a large optical depth and receives very little radiative acceleration regardless of resolution. Overall, these tests support that a particle resolution coarser than our chosen resolution of 10⁻⁴ M_⊙ will not accurately capture the generation of a dusty wind by radiation pressure.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

2.3.2. Anisotropy of Radiation Field

Within the hydrodynamic simulation, the AGN flux is assumed to be anisotropic, following

$$\begin{eqnarray}&&F(r,\theta )=\displaystyle \frac{L}{4\pi {r}^{2}}f(\theta ),\end{eqnarray} \tag{ 7 }$$

where L is the luminosity of the AGN, r is the distance from the AGN, θ is the angle from the polar axis of the AGN, and the anisotropy function f(θ) is defined to be

$$\begin{eqnarray}&&f(\theta )=\displaystyle \frac{1+a\cos \theta +2a{\cos }^{2}\theta }{1+2a/3},\end{eqnarray} \tag{ 8 }$$

where we define $a=({\eta }_{a}-1)/3$, introducing the parameter η_a as the "anisotropy factor," equal to the ratio between the polar flux and the equatorial flux. This is a generalization of the classic angle dependence of Netzer (1987), to account for deviations from perfect limb darkening of a geometrically thin disk. It is effectively a superposition of an isotropic component with the classic limb-darkened disk. This chosen generalization produces only a small variation in polar flux with η_a, despite a large variation in equatorial flux.

The anisotropy factor η_a is not well constrained by observations and strongly depends on the model assumed for the accretion disk and corona. It is likely to vary by wavelength, as emission from the corona should be more isotropic than emission from the accretion disk. In this work, we make the simplification that η_a is independent of wavelength, as this allows us to assume a constant SED shape before extinction is included. We also assume that η_a is a free parameter and vary it between the simulations.

Accretion of dusty gas along the disk plane is possible at large ${\eta }_{a}$ and low Eddington factors γ_Edd, where the equatorial flux is low. However, in test simulations we found that the narrow region where dusty accretion is possible was not properly resolved, even at the very fine particle resolution of m_P = 10⁻⁴ M_⊙. The result was that radiation pressure would blow away most of the inner surface of the disk (see also Chan & Krolik 2016), but a few particles would be left behind in the disk plane—too few to be properly resolved in the SPH kernel. To avoid this problem, we ensure that η_a is low enough and γ_Edd high enough that dusty accretion is not possible. Gas, however, can avoid being blown away if the dust is destroyed/sublimated.

2.3.3. IR Radiation from the Dusty Gas

We use version c17.00 of the photoionization code CLOUDY to incorporate the detailed microphysics that determines the ionization, chemical, and thermal states through a slab of dust and gas exposed to an external radiation field. Our CLOUDY models take three parameters as input: the volumetric density of the cloud, the incident AGN radiation field, and the temperature of the gas. The models are in non-local thermal equilibrium, allowing the dust temperature to reach equilibrium with the radiation field and the gas temperature. Each run produces outputs as a function of hydrogen column density, up to 10²⁶ cm⁻².

It has been suggested that the properties of AGN winds do depend on the details of the driving SED (Dyda et al. 2017), but we focus on a single SED in this work. For the incident radiation field we assume an AGN SED with a big blue bump component at a temperature of T = 10^5.5 K and a logarithmic X-ray-to-UV flux ratio of −1.40. The slope of the big blue bump and the slope of the X-ray component are left at the CLOUDY default values of −0.5 and −1.0. This produces less inverse Compton scattering of X-ray photons to radio wavelengths than CLOUDY's built-in AGN continuum command and is the preferable method as per the CLOUDY handbook. We also include the cosmic microwave background, the cosmic-ray background, and a stellar background with an integrated intensity of 1000 Habing units, which accounts for the nuclear clusters found in the centers of galaxies. This stellar background intensity is taken from Wada (2012), who note that the exact value does not greatly affect the dynamics of the circumnuclear region (Wada et al. 2009).

We adopt the built-in ISM dust grain model in CLOUDY and use the command function sublimation to account for sublimation where applicable. The grain model includes graphites and silicates with sizes following the MRN distribution (Mathis et al. 1977), ranging from 0.005 to 0.25 μm. We do not include PAHs in our calculations, as there is an expectation that these very small grains are destroyed close to the AGN (but see Jensen et al. 2017). We define the dust temperature used in visualizations as the cross-section-weighted average over all dust sizes. For T > 10⁵ K, we assume that the dust has been destroyed by sputtering rather than sublimation and do not include dust at those high temperatures.

The table covers gas temperatures $\mathrm{log}T\,({\rm{K}})=1\mbox{--}8$, gas densities $\mathrm{log}n\,({\mathrm{cm}}^{-3})=0\mbox{--}9$, and AGN fluxes $\mathrm{log}I\,=1.2\mbox{--}9.2$, where I is in units of erg cm⁻² s⁻¹. However, for $\mathrm{log}n\,({\mathrm{cm}}^{-3})\gt 7$ and $\mathrm{log}T\,({\rm{K}})\gt 3$, the CLOUDY simulations did not converge. Such hot dense gas is rare, and so we extrapolate the table for the few gas particles in this region.

For a given set of the three input values, we compute a total of eight physical quantities: the heating and cooling rates, the temperature of the dust, the dust-to-gas mass ratio, the radiative acceleration due to the incident radiation, the absorption and scattering opacities, and the optical depth. The radiative acceleration, intensity-weighted opacities, and optical depth were extracted by calling repeatedly the save continuum emissivity command for each of the 5277 frequencies considered in CLOUDY¹ and using a post-processing routine to integrate over the resulting spectrum. The spectrum is sampled in the interval 1.296 × 10⁻⁸–9.080 × 10⁶ ryd with nearly logarithmically increasing widths.

Gravitational acceleration is calculated by the sum of an analytic background term (Section 2.6) and a self-gravity term. Brute-force gravitational acceleration could be calculated through

$$\begin{eqnarray}&&{{\boldsymbol{g}}}_{i}=\displaystyle \sum _{j}\displaystyle \frac{{{Gm}}_{j}{{\boldsymbol{r}}}_{{ij}}}{{r}_{{ij}}^{3}}+{{\boldsymbol{g}}}_{\mathrm{bg}},\end{eqnarray} \tag{ 9 }$$

where G is the universal gravitational constant, ${{\boldsymbol{r}}}_{{ij}}$ is the displacement between particles i and j, and ${{\boldsymbol{g}}}_{\mathrm{bg}}$ is the analytic background acceleration. A Barnes–Hut oct-tree (Barnes & Hut 1986) is used to reduce the size of the sum.

Cold, dense, self-gravitating gas tends to be susceptible to gravitational collapse. In the absence of a feedback mechanism, this collapse continues until stars are formed. The smoothing length of the gas is adaptive, and so we can follow this collapse to arbitrary densities without violating the Truelove criterion (Truelove et al. 1997). However, as we note in Section 2.3, heating, cooling, and radiation pressure must be extrapolated once the density is higher than we have tabulated. As we additionally note, there is also a threshold density beyond which an individual particle's optical depth is greater than unity and the important skin region is not resolved.

To prevent this in early tests, the surface density was set so that the Toomre parameter (Toomre 1964) was Q = 2 at all radii, and hence the gravitational instabilities that produce star-forming clumps are suppressed. In these lower-resolution tests, we found that radiation pressure nevertheless drove gas to high enough densities that it became Jeans unstable and underwent unimpeded collapse. To prevent this, we implemented a star formation routine that converts gas above a threshold density into stars at the freefall timescale. However, we found that the high resolution of our production runs is sufficient for shear and turbulence to halt the collapse, and star formation never occurs within the simulated disk. This implies that we would not expect significant star formation within the inner few parsec around an AGN, while it may occur at slightly larger distances.

We have also implemented a supernova (SN) feedback model, using a kinetic feedback method, which we use in some of our simulations. Our simulation time is shorter than the lifetimes of even the most massive stars, and so we cannot self-consistently calculate the rate and location of SNe based on resolved star formation. Instead, we place SNe at a fixed rate in an evenly weighted random location in the disk within a distance R_SN of the AGN center, where R_SN = 0.4 pc is a model parameter. The kinetic energy of an SN is ${\epsilon }_{\mathrm{SN}}{E}_{\mathrm{SN}}$, where ${E}_{\mathrm{SN}}={10}^{51}$ erg is the fiducial energy output of an SN and _SN is the fraction of energy that is converted into kinetic energy and is a model parameter, which we set to _SN = 0.1. SNe are assumed to occur instantaneously, and an SN event affects all n particles within r_SN by giving each particle a radially directed kinetic energy kick ${\epsilon }_{\mathrm{SN}}{E}_{\mathrm{SN}}/n$, where r_SN = 0.01 pc is a model parameter.

We model the initial conditions for the gas in the inner region of the torus as a disk of mass M = 69.0294 M_⊙ and thickness 0.02 pc, with an inner radius of 0.005 pc and an outer radius of 0.47 pc. The surface density is set to 100 M_⊙ pc⁻², giving an initial density of ρ = 5000 M_⊙ pc⁻³ (n_H ≈ 1.7 × 10⁵ cm⁻³). The gas disk has a uniform initial temperature of 1000 K, but heating and cooling are sufficiently rapid in our simulations that the evolution is not sensitive to the initial temperature. At our high resolution we cannot model the entire dusty region out to several tens of parsecs, but instead we concentrate on examining the inner region in detail. For comparison with a full-scale "torus," if the initial surface density of the disk was extended to a radius of 20 pc, it would have a total mass of 130,000 M_⊙. The gas is modeled by 690,294 particles, each with a mass of 10⁻⁴ M_⊙.

The imposed gravitational potential is a superposition of a softened Keplerian term to represent the black hole and a Hernquist bulge term to represent the potential of the galactic stars and dark matter. In practice, the Keplerian term dominates for most particle orbits. The form of the gravitational force law is

$$\begin{eqnarray}&&F(r)/G=\displaystyle \frac{{M}_{\mathrm{BH}}{r}^{2}}{{r}^{2}+{c}_{\mathrm{BH}}^{2}}+\displaystyle \frac{{M}_{{\rm{H}}}{\left(r/{c}_{{\rm{H}}}\right)}^{2}}{{\left(1+r/{c}_{{\rm{H}}}\right)}^{2}},\end{eqnarray} \tag{ 10 }$$

where M_BH and M_H are the masses of the black hole and Hernquist bulge, respectively, c_BH is the softening length for the black hole potential, and c_H is the scale length for the Hernquist bulge. We keep most of the parameters constant in this suite of simulations, setting M_H = 10⁹ M_⊙, ${c}_{\mathrm{BH}}={10}^{-4}$ pc, and c_H = 250 pc, but vary the black hole mass with γ_Edd so that the AGN luminosity is constant at 1.26 × 10⁴² erg s⁻¹. We set ${M}_{\mathrm{BH}}={10}^{6}\,{M}_{\odot }{({\gamma }_{\mathrm{Edd}}/0.02)}^{-1}$. The particles are given velocities such that the initial rotation curve is in equilibrium with the potential, although the thickness of the disk is primarily supported by pressure. Under strong radiation pressure, this means that the gas disk is not in equilibrium, but is pushed outward. Hence, when radiation pressure is strong, the simulations do not represent a long-term equilibrium, but rather the short-term response of a gas disk to the AGN radiation field.

We perform multiple simulations from γ_Edd = 0.02 to γ_Edd = 0.2. We set η_a = 10² in most of the simulations, except for the simulations with γ_Edd = 0.1, where we also perform simulations with η_a = 10⁰ (i.e., isotropic), η_a = 10¹, and η_a = 10³. Most simulations have a temperature floor of 10 K (which is never reached in the simulation), but in some we introduce higher-temperature floors of T = 30, 300, and 3000 K to see the effects of a thicker disk, as we might expect to be produced by the disk puffing itself up with IR radiation pressure. We also produce one simulation ("NoAgn") with no AGN radiation field, where gas is evolved adiabatically. We run each simulation for 9.78 kyr, and all systems either reach a steady state or are blowing out the gas by the end of the simulation.

We use a naming scheme that represents the radiation anisotropy and Eddington factor for each run, as well as the presence of SN feedback or an increased temperature floor. These simulation parameters are summarized in Table 1.

Table 1.  Summary of Run Parameters

			Tfloor	SN Rate	${F}_{p}/\langle F\rangle$
Name	γEdd	ηa	(K)	(Myr−1)
NoAGN	N/A	N/A	101	0	N/A
a2_e01	0.01	102	101	0	4.35
a2_e02	0.02	102	101	0	4.35
a2_e05	0.05	102	101	0	4.35
a2_e1	0.10	102	101	0	4.35
a2_e2	0.20	102	101	0	4.35
a2_e01_SN100	0.01	102	101	102	4.35
a2_e01_SN1000	0.01	102	101	103	4.35
a2_e01_T30	0.01	102	3 × 101	0	4.35
a2_e01_T300	0.01	102	3 × 102	0	4.35
a2_e01_T1000	0.01	102	103	0	4.35
a0_e1	0.10	100	101	0	1.00
a1_e1	0.10	101	101	0	3.33
a3_e1	0.10	103	101	0	4.49

Note. Here γ_Edd is the Eddington factor, η_a is the anisotropy factor (i.e., the ratio between the polar flux and the equatorial flux), T_floor is the temperature floor, and ${F}_{p}/\langle F\rangle$ is the ratio between the polar flux and the isotropic average flux.

Download table as:  ASCIITypeset image

All of the AGN runs follow a similar evolution, and so here we focus on the properties and evolution of Run a2_01 and only describe the other runs when variations appear. The evolution of the gas density and temperature for Run a2_e01 are plotted in Figure 5. The gas disk remains thin and cool, because the heating from the AGN is deposited only in a thin skin at the inner edge of the disk and is not reradiated inward. Gravitational instabilities are visible as ring and spiral patterns, but these do not greatly thicken the disk. In the absence of radiation pressure and of nonadiabatic heating and cooling, the gravitational potential causes the system to evolve into a flared disk (Figure 6) with no outflow.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A hot wind is formed. The thin skin region on the inner edge of the disk receives strong radiation pressure in addition to heating, causing the inner gas to be pushed outward as a warm wind. The thermal pressure of the gas pushes it vertically, thickening the outflow. Although the gas of this wind reaches temperatures of ∼10⁴ K, the dust component is significantly cooler and is not destroyed by sublimation or sputtering. The opacity of the gas therefore remains large, and the gas continues to be pushed outward. We expect that the inclusion of reradiation of IR radiation between gas particles would provide an additional vertical force through radiation pressure, and we will explore this in a future paper.

The inner radius of the disk slowly expands throughout the simulation. This is partially a result of the inner edge of the disk being ablated through heating and radiation pressure, but it is also a result of the net momentum deposition on the disk as a whole from radiation pressure. Most of the deposited momentum is carried off in the wind, but some is transferred to the bulk of the disk, pushing the disk outward. This effect is stronger when the intensity of radiation in the plane of the disk is stronger, especially in Run a0_e1. The rotation curve of the gas disk is initially in equilibrium with the gravitational potential and does not represent a self-consistent long-term state of the system, but rather the short-term response of a disk exposed to AGN radiation. We will explore more self-consistent models with gas accretion in a later paper.

To characterize the wind, we must precisely define which gas particles are considered to be "wind particles." In this work, we define the wind as all gas particles with velocities greater than the escape velocity from the gravitational potential at their position, within 5 pc of the simulation center. This cutoff radius is used to prevent gas expelled early in the simulation from dominating the results. We define the outflow rate ${\dot{M}}_{w}$ as the rate at which particles join the wind, plotted in the top row of Figure 7.

Zoom In Zoom Out Reset image size

Considering only the runs with η_a = 10², the runs with a greater Eddington factor show a higher outflow rate, quickly reaching a steady state. This suggests a fairly direct relationship between Eddington factor and outflow rate, and we plot the outflow rate at t = 2 kyr against Eddington factor in the bottom row of Figure 7. Here we find a power-law relationship between outflow rate and Eddington factor that is slightly sublinear, with a power index of ∼0.8.

More critically, we find that varying the anisotropy of the radiation field has a dramatic effect. It has sometimes been assumed that no radiation at all is emitted in the AGN plane, for at least part of the SED (such as in the UV, e.g., Wada et al. 2016), or conversely that the radiation field is anisotropic across the entire SED (e.g., Krolik 2007; Dorodnitsyn et al. 2016). We should expect limb darkening for radiation emitted by the AGN accretion disk, but we should not expect the radiation in the plane to be absolutely zero. The disk is not likely to be perfectly flat and perfectly thin, and we should expect a contribution of reflected radiation from the corona as well. We have varied the equatorial emission from our fiducial value of 1% of the polar emission (the a2_e* runs) to 0.1% in Run a3_e1 up to 100% in Run a0_e1. We find a strong dependence of outflow rate on anisotropy. This is particularly dramatic in a1_e1 and a0_e1. Here the more isotropic radiation field produces a very large radiation pressure on the inner surface of the disk, producing a strong outflow that blows away the gas disk entirely. In the bottom row of Figure 7 we plot the anisotropy factor against outflow rate and find a power-law trend with an index of −∼0.4. This demonstrates that decisions on the anisotropy of the AGN radiation field cannot be taken lightly, as sensitivity of the outflow rate to the anisotropy is of a similar order to its sensitivity to the Eddington factor.

We can see related effects in the evolution of the wind angle, which is defined as the mean latitude of all wind particles, plotted in the middle row of Figure 7. The wind angle is higher when the initial outflow rate is lower. This is caused by the interaction between thermal pressure, causing the wind to expand above and below the disk, and radiation pressure, pushing the gas outward. The thermal pressure of the wind particles and the disk does not vary greatly between the runs, despite the disk receiving different fluxes. This is the result of negative feedback from radiative cooling that causes the gas to be attracted toward discrete phases of temperature and density—it requires a large amount of energy to push the gas beyond T = 10⁴ K. However, there is no feedback effect to reduce radiation pressure until dust is destroyed, and in this regime radiative acceleration essentially increases in proportion to the flux on the inner disk edge. Thus, at lower anisotropies, the gas receives greater radiation pressure but still has a similar thermal pressure, and the wind is pushed outward more radially, producing a lower wind angle. Similarly, at higher Eddington factors, the gravitational potential is stronger while the radiation pressure is unchanged, and the outflows also become more equatorial.

The interaction between thermal and radiation pressure is also demonstrated in the simulations with varying temperature floor, in the right panel of Figure 7. As a reminder, we use the temperature floor as a means to approximate the effect of IR radiation pressure puffing up the dusty disk, which we aim at modeling directly in the upcoming paper. As the wind-launching gas is rapidly heated, the wind is not very sensitive to the previous temperature of the gas. At T_floor = 1000 K we do see an increase in the wind angle, due to thermal pressure producing a slightly thicker disk, but this effect is small.

SNe also only have a relatively small effect. Each SN is visible as a significant but short-lived peak in the outflow rate. In a2_e01_SN1000 the SN rate is high enough that the disturbance caused by SNe also causes the outflow rate to somewhat increase in between SNe, although the increase in mass outflow is smaller than that caused by lowering the anisotropy. The wind angle also becomes significantly more polar, although as we note in Section 3.4, this represents a small amount of gas thrown to high angles rather than a significant puffing up of the wind in general.

In these small-scale simulations, the coupled energy input of an SN (10⁵⁰ erg; see Section 2.5) is greater than the initial kinetic energy of the disk (∼10⁴⁹ erg). SNe therefore have an extremely small effect compared to their huge energy input. This is because the energy input is not efficiently coupled to the disk as a whole. Instead, as the energy deposition region is well resolved, it forms a bubble of hot gas and immediately escapes in a "chimney." Rather than boosting the wind and disk in general, a small amount of hot gas is ejected at high speeds.

We also note that, assuming 1 SN per ∼100 M_⊙ of star formation, the rate of 1000 SNe Myr⁻¹ is the equivalent of a star formation rate of about 0.1 M_⊙ yr⁻¹ concentrated within a ∼0.47 pc disk, or a star formation surface density of >4 × 10⁵ M_⊙ yr⁻¹ kpc⁻². This modest effect on the disk and wind is therefore the result of an extremely powerful nuclear starburst that we should not expect to be ubiquitous in obscured AGNs.

3.4.1. Column Density Profiles

AGN unification requires column densities with a large covering fraction, and so the column density as a function of angle above the disk plane is a critical constraint for AGN torus models. Although it is not possible to observe a single AGN from multiple angles, the general shape of the column density distribution can be inferred statistically from observing a large number of AGNs.

We calculate the column density for each angle by tracing rays from the center of the simulation outward and summing the contribution to column density across all particles that intersect the ray. We convert the mass surface density into a number surface density by assuming a mean molecular mass of ≈1.22. For each of 40 different inclinations, we trace 40 rays, equally spaced azimuthally. The mean, maximum, and minimum of each set of rays for all runs at t = 1.96 kyr and t = 9.78 kyr are plotted in Figure 8. We plot the earlier time because the gas disk is blown away at later times in some runs and the column densities become very small, and we plot the later time because the effects of SNe do not become apparent until later in the simulation.

Zoom In Zoom Out Reset image size

The simulations all have remarkably similar column density profiles, but there are some trends to be observed. Runs a0_e1 and a1_e1 have less material at very low inclinations as the disk is blown away. We see a related effect when comparing runs with different Eddington ratio. Here the column density slowly decreases with increasing Eddington ratio, except at ϕ > 20°, where the trend is perhaps reversed, although this represents only a very small amount of gas.

As above, the temperature floor has little effect. This indicates that the details of the hydrodynamics of the disk are not important for the wind—in the wind generation region, the gas properties are dominated by the radiation of the AGN and the intervening opacity, and small changes in disk thickness are unimportant.

SNe do have an effect at large inclinations, propelling small amounts of gas and dust to large heights about the disk. The effect here is not to greatly change the opening angle of the optically thick part of the wind, but to provide a small level of baseline extinction at more polar viewing angles.

3.4.2. Kinematics

Specific observational tracers are typically sensitive to gas or dust at some particular range of temperatures. To demonstrate what dynamical phases a particular observation might capture, phase plots of gas temperature, dust temperature, and radial velocity from the simulation center of Run a2_e01 at t = 1.96 kyr are plotted in Figure 9. The gas mostly lies in two distinct phases, with some gas in intermediate states—wind gas in the process of being heated. The bulk of the gas mass remains in the disk, with a radial velocity of $| {v}_{\mathrm{rad}}| \lesssim 10$ km s⁻¹, a gas temperature of hundreds of kelvin, and a dust temperature of tens of kelvin. Most of the rest of the gas mass is either transitioning into wind or already in the wind—a warm ionized phase with gas temperatures of just under 10⁴ K.

Zoom In Zoom Out Reset image size

The disk phase is clearly visible in the phase plots in Figure 9 as having low dust and gas temperatures. We can isolate this phase by applying cuts of either gas or dust temperature, as plotted in Figure 10. We choose T_d = 100 K and T_g = 1000 K as threshold temperatures and plot the velocity fields of gas above and below these thresholds in Figure 10. Both the gas and dust temperature cuts select qualitatively similar velocity fields, which is expected as there is a correlation between dust and gas temperature, at low gas temperature. Material with temperature below either cut shows strong rotation, while material above either cut is strongly outflowing. This suggests that observational tracers of hot/cold gas will respectively trace hot/cold dust and vice versa, and that tracers of hot material should observe outflowing gas and tracers of cold material should observe indicators of rotation. However, as we note below, this coupling is not exact.

Zoom In Zoom Out Reset image size

3.4.3. Mock Images

In Figure 11, we have plotted ray-traced mock images of Run a2_e01 from various angles. For each ray—corresponding to a pixel in each image—intersecting particles are sorted by distance along the line of sight. Summing from nearest to farthest, each particle contributes a flux equal to ${\sigma }_{{\rm{B}}}{f}_{d}{T}_{d}^{4}{e}^{-\tau }$, where σ_B is the Stefan–Boltzmann constant; f_d and T_d are the dust mass fraction and dust temperature, respectively, as calculated with the CLOUDY tables; and τ is the optical depth from all previous particles along the line of sight. Here we do not use the CLOUDY table opacities, as these are weighted by the flux of the AGN. Instead, for this visualization we use an opacity of 65.2 cm² g⁻¹, which is typical for ISM dust. The normalization of the plotted flux is arbitrary.

Zoom In Zoom Out Reset image size

The model images reproduce some features that are observed in IR interferometry (Hönig et al. 2012, 2013). We find a bright ring of hot dust in the center of the disk, as observed in near-IR. In the edge-on view, this is no longer visible, but is blocked by extinction from the disk, and most of the emission appears above and below the disk, representing the warm wind. This warm wind corresponds to the polar extended mid-IR emission, although in our simulations the wind does not extend as far in the polar direction as interferometry suggests owing to the lack of IR reradiation.

Recently, ALMA observations have produced high-resolution imagery of nearby AGNs on slightly larger than parsec scales (e.g., García-Burillo et al. 2016, 2017; Alonso-Herrero et al. 2018; Imanishi et al. 2018; Izumi et al. 2018), observing molecular lines from molecules such as HCN and CO. Using CLOUDY, we have produced tables of emissivities for two CO lines, two HCN lines, and three vibrational transitions of H₂. The CO and HCN lines lie within the range observable by ALMA, while the H₂ lines are observable by the VLT and VLTI instruments (H₂ 1−0 S(1) and H₂ 0−0 S(3)), or will be in the future with the James Webb Space Telescope ((H₂ 1−0 S(0), although only at low redshift). The emissivities as a function of dust temperature and gas density are plotted in Figure 12 for Run a2_e01 at t = 9.78 kyr, the end of the simulation. Generally, the lines are stronger at high gas densities and low dust temperatures, with a rapid transition from strong emissivity to zero emissivity as the molecules are destroyed in hotter, lower-density environments. For these lines, dust temperature is a proxy for the intensity of radiation. We note above that there is a correlation between dust temperature and gas temperature, but the molecular emission is more closely tied to dust temperature than to the gas temperature, and the emissivity cutoff is less clean in a gas temperature–density plot.

Zoom In Zoom Out Reset image size

All of the lines trace high-density gas with cool dust, but they differ in total strength and in the slope and location of the transition to nonemitting. The CO lines and the H₂ (1−0) S(1) and H₂ (0−0) S(3) vibrational lines retain emission in higher temperatures and lower densities than the HCN lines, with the H₂ (0−0) S(0) line lying somewhere in between. To show how this might appear in high-resolution observations, we have plotted the fluxes from three of these lines in Figure 13, along with a combined RGB image. The flux in each pixel is simply the sum of the column density multiplied by the tabulated emissivity per particle within one smoothing length of the line of sight through the pixel, and it does not include extinction effects. Here it can be see that the HCN lines trace material closest to the disk, including part of the outflowing wind. The CO line traces material farther from the disk, and the H₂ 1−0 S(1) vibrational line emphasizes the extended material even more. The expectation from these simulations is that the H₂ 1−0 S(1) lines show a higher velocity dispersion than the HCN and CO lines, which is consistent with observations of these lines (e.g., Hicks et al. 2009, and references therein).

Zoom In Zoom Out Reset image size