{'fraction_non_alphanumeric': 0.03898863274195262, 'fraction_numerical': 0.0190906193562095, 'mean_word_length': 4.238813446126447, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Geometrically thick obscuration by radiation-driven outflow from magnetized tori of active galactic nuclei

Near-Eddington radiation from active galactic nuclei (AGNs) has significant dynamical influence on the surrounding dusty gas, plausibly furnishing AGNs with geometrically thick obscuration. We investigate this paradigm with radiative magnetohydrodynamics simulations. The simulations solve the magnetohydrodynamics equations simultaneously with the infrared (IR) and ultraviolet (UV) radiative transfer (RT) equations; no approximate closure is used for RT. We find that our torus, when given a suitable sub-Keplerian angular momentum profile, spontaneously evolves toward a state in which its opening angle, density distribution, and flow pattern change only slowly. This"steady"state lasts for as long as there is gas resupply toward the inner edge. The torus is best described as a mid-plane inflow and a high-latitude outflow. The outflow is launched from the torus inner edge by UV radiation and expands in solid angle as it ascends; IR radiation continues to drive the wide-angle outflow outside the central hole. The dusty outflow obscures the central source in soft X-rays, the IR, and the UV over three quarters of solid angle, and each decade in column density covers roughly equal solid angle around the central source; these obscuration properties are similar to what observations imply.

Whatever dynamical processes operate, they must explain the structure of these parsec-scale dusty regions as detected by mid-infrared (MIR) interferometry in nearby AGNs. Some AGNs reveal two components, one elongated in the equatorial direction and the other in the polar direction (e.g., Jaffe et al. 2004; Raban et al. 2009;Tristram et al. 2014). The equatorial component is understood as the torus inner edge (e.g., Jaffe et al. 2004), while the polar component could come from optically thin dust in the polar regions under direct illumination by the central source (Hönig et al. 2012). Polar emission is also seen in MIR imaging on 100 pc scales (e.g., Braatz et al. 1993). Circinus is a particularly striking example: MIR emissions on ∼ 1 pc and ∼ 100 pc are aligned with each other (Asmus et al. 2016), and also with the edge of the ionization cone (Packham et al. 2005;Asmus et al. 2016). This observation suggests that dust is concentrated along the edge of the cone (Hönig et al. 2012;Asmus et al. 2016), which could be the case if radiation pressure on dust drives a dusty outflow from the inner surface of the torus (Hönig et al. 2012).
Several authors have tried to explain dusty outflows. Roth et al. (2012) performed Monte Carlo radiative transfer (RT) on a static distribution of dusty gas and found that radiative acceleration can overcome gravity. Keating et al. (2012) constructed analytic solutions of a magnetocentrifugal wind where radiation pressure on dust provides additional driving. Detailed radiative hydrodynamics (RHD) simulations have also been performed. The simulations by Wada et al. (2016) considered momentum deposition by ultraviolet (UV) radiation on dust, but ignored momentum transfer between infrared (IR) radiation and dust. Dorodnitsyn et al. (2016) took the opposite approach by including only momentum coupling between IR radiation and dust. CK16 was the first attempt to capture the effects of both IR and UV radiative support in one single simulation. Despite their significant differences, all these studies agree that, for a torus orbiting a central object of mass M ∼ 10 7 M radiating at an Eddington ratio L/L E ∼ 0.1, the typical mass loss rate due to radiation-driven dusty outflow is 0.1 M yr −1 Ṁ 1 M yr −1 . This is not a feeble outflow: A torus of Thomson optical depth τ T located just outside the dust sublimation surface can be depleted by mass loss in This means a torus cannot remain in a steady state for many orbits unless there is a constant resupply of gas from galactic scales through the torus down to the inner edge (see also Krolik & Begelman 1988;Roth et al. 2012;Schartmann et al. 2014;Wada 2015). Magnetohydrodynamic (MHD) stresses generated by the magnetorotational instability (MRI) are widely accepted as the mechanism for angular momentum transport and mass accretion in geometrically thin disks (e.g., ; such stresses conceivably regulate accretion in geometrically thick disks as well. The MRI grows in the ideal MHD condition, which holds even when the gas is weakly ionized (Blaes & Balbus 1994;Gammie 1996). Such low levels of ionization can be sustained in the torus interior by X-rays (Neufeld & Maloney 1995) since the energy produced by AGNs in X-rays is ∼ 0.05 times that in the UV (e.g., Zamorani et al. 1981). Indeed, magnetic field of strength ∼ 10 mG has been detected on 30 pc scales in the nucleus of NGC 1068 (Lopez-Rodriguez et al. 2015). It is therefore critical that we understand what implications magnetic field has on torus dynamics.
This article reports the extension of our RHD simulations (CK16) to radiative magnetohydrodynamics (RMHD). Our code marries the finite-volume MHD code Athena (Stone et al. 2008) with time-dependent IR (Jiang et al. 2014) and time-independent UV (CK16) RT modules. It is uniquely capable of concurrently solving the MHD and RT equations without resorting to arbitrary closures for RT; this property is crucial for correctly treating gasradiation interaction in systems with optical depths comparable to unity (Davis et al. 2012).
In addition, we quantitatively explore the effect of varying the rotational profile of the torus. An important lesson is that a torus with Keplerian rotation cannot withstand irradiation for a long time ( §3.2), but a torus with sub-Keplerian rotation may maintain a steady inflow-outflow morphology and survive under mass loss for multiple orbits ( §4.2). The kinematics and obscuration properties of the radiation-driven outflow in this "steady" state generally agree with observations. We recount our simulation setup in §2, report our results in § §3 and 4, and discuss them in §5.

METHODS
Our simulations are based on CK16; the only differences are in our initial condition ( §2.3.1) and in how we conduct the simulations ( §2.3). We adopt cylindrical coordinates (R, φ, z) and define the spherical radius r ≡ (R 2 + z 2 ) 1/2 for convenience. Quantities are normalized by the fiducial quantities listed in Table 1. The surface on which the UV optical depth τ UV to the central source is unity is called the "inner surface," whereas the part of this surface near the mid-plane is called the "inner edge."

Magnetohydrodynamics
The equations of ideal MHD are ∂ρ ∂t + ∇ · (ρv) = 0, Here ρ, v, and p are gas density, velocity, and pressure. The magnetic field is B, whose unit is chosen such that magnetic permeability is unity. Gas temperature, total pressure, and total energy density are T = p/(ρR ideal ), p * = p + 1 2 B 2 , and E = 1 2 ρv 2 + p/(γ − 1) + 1 2 B 2 respectively, where R ideal and γ are the specific ideal gas constant and the ratio of specific heats. The gravitational potential of the central mass is Φ(r) = −GM/r. The energy and momentum source terms due to radiation are S e IR,UV and S m IR,UV , to be defined in §2.2. Finally, the isotropic rank-two tensor is denoted by I.

Radiative transfer
The torus is illuminated in the UV by the innermost regions of an accretion disk at the origin. The angular distribution of the UV radiative flux is modulated by a number of physical effects: Geometrical projection and limb darkening due to a scattering atmosphere favor emission in the polar direction, while relativistic boosting, beaming, and lensing enhance the equatorial radiative flux. For simplicity, we adopt a spherically symmetric central source in our simulations.
The propagation of IR radiation is handled by the timedependent RT module of Athena, which solves the timedependent RT equation on a large number of grid rays. To first order in v/c, where c is the speed of light, the mixed-frame time-dependent RT equation for IR radiation interacting with gray material reads (Jiang et al. 2014) 1 c The specific intensity integrated over the IR in the observer frame is I IR (n); its lowest three angular moments are J IR , H IR , and K IR , from which the IR radiation energy density and flux follow as E IR = (4π/c)J IR and F IR = 4πH IR . The frequency-integrated black-body mean intensity is B(T ) = ca SB T 4 /(4π), where a SB is the radiation constant. The comoving absorption and scattering cross sections per mass in the IR are κ IR and σ IR respectively, andĉ is the reduced speed of light (CK16; see also Gnedin & Abel 2001;Skinner & Ostriker 2013). Taking the zeroth and first angular moments of Equation (6) yields 1 c The remaining piece to specify in Equations (6)- (8) is H 0 IR , the first angular moment of the IR specific intensity in the fluid frame. It is related to the angular moments in the observer frame by a Lorentz transformation (Mihalas & Weibel-Mihalas 1984): The time-dependent RT module of Athena introduces artifacts along grid rays in optically thin regions; therefore, we have written a time-independent long-characteristics RT module specifically for UV radiation from the central source. The module computes the UV radiation energy density in any cell as where L UV is the luminosity of the central source, while τ UV and τ * UV are the UV optical depths from the source to the cell and across the cell respectively. The energy and momentum source terms due to UV radiation are then with κ UV being the comoving absorption cross section per mass in the UV. The chief sources of opacity in our system are dust absorption and electron scattering, which we model as In these fitting formulae, T ds ≈ 1500 K is the dust sublimation temperature (e.g., Rees et al. 1969;Rieke & Lebofsky 1981;Barvainis 1987), T hi ≈ 4013 K is the temperature at which hydrogen atoms in local thermodynamic equilibrium at a number density of 10 4 cm −3 are collisionally half-ionized, and κ T ≈ 0.397 cm 2 g −1 is the Thomson scattering cross section per mass. The dust opacities are normalized to Thomson as κ IR /κ T = 20 andκ UV /κ T = 80; the parameters governing the transition between opacity regimes are ∆ ds = 0.05 and ∆ hi ≈ 0.196.

Simulation strategy
MHD turbulence stirred by the MRI saturates after several tens of orbits, but RMHD simulations are computationally prohibitively expensive on that timescale. One solution is to carry out the simulation in stages; in order of execution, they are the MHD stage, the UV-RMHD stage, and the IR-RMHD stage. Each stage is documented at length below.

MHD stage
In the MHD stage, we follow a geometrically thick, gassupported torus in pure MHD to saturation. The rationale is that if IR radiation does behave like a pressure in optically thick regions, we should be able to replace gas pressure by IR radiation pressure at a later stage with minimal changes to the geometry of the torus. We begin by constructing the initial condition for the MHD stage. Assuming v(R) = (GM/R p ) 1/2 (R/R p ) 1−qê φ and p = Kρ Γ , the density of an axisymmetric hydrostatic torus is uniquely given by (Papaloizou & Pringle 1984) with maximum ρ = ρ p at (R, z) = (R p , 0). The free parameters of the initial condition are the radial coordinate R p of the density maximum, the shear parameter q, the polytropic constant K, the polytropic index Γ, and the constant on the left-hand side. Note that we require 1.5 < q < 2 for the torus to have finite height and be stable; indeed, a geometrically thick torus supported by isotropic pressure must have a sub-Keplerian rotational profile ( §2.3.3). A geometrically thick torus with bounded radial and vertical extent can be obtained with R p = r 0 , ρ p = ρ 0 , q = 1.9, K = 0.2p 0 ρ −Γ 0 , and Γ = 4. The surface of the torus is ρ(R, z) = 0; solving this equation for z = 0 yields R ≈ 0.545 r 0 and R ≈ 4.72 r 0 .
Overlaid on the torus is a poloidal loop of magnetic field derived from the vector potential A ∝ max(ρR ξ − Cρ 0 r ξ 0 , 0)ê φ ; we choose an exponent ξ = 0.4 so that the plasma betas at the inner and outer surfaces of the torus are similar, and a cutoff C = 0.4 so that all field lines are properly confined within the torus. The proportionality constant is selected to make the ratio of volume-integrated gas pressure to volume-integrated magnetic energy equal to 1000. Gas pressure is also perturbed at the 0.01 level to seed the MRI. We set up the ambient material around the torus as in CK16, but here its length and density scales are R amb = r 0 andρ amb = 10 −5 ρ 0 respectively, and its sound speed is v φ (R amb ).
The configuration is evolved until MHD turbulence has fully saturated, at t = 200 t 0 . MHD stresses have increased the specific angular momentum j ≡ Rv φ at (R, z) = (2 r 0 , 0) from ≈ 0.76 to ≈ 0.81 times Keplerian, and steepened its mid-plane radial profile from j(R) ∝ R 0.1 to Meanwhile, the torus outer edge spreads radially outward and the torus loses mass through all boundaries, but eventually the torus stabilizes. To make up for the mass loss, we multiply ρ 1/2 , p 1/2 , and B at the end of the MHD stage by the same factor everywhere, determined as follows. The Thomson optical depth averaged over all sightlines contained within a mid-plane wedge of aspect ratio m 1 is where ϑ(·) is the step function and φ min,max are the coordinates of the azimuthal boundaries of the simulation domain. The factor is chosen such that τ T of a wedge with m = 0.1 matches some desired value, which we shall give in §2.3.3.

UV-RMHD stage
We now switch over to RMHD. When the central source is turned on and radiation starts pushing on the gas at the inner surface, the torus is no longer in equilibrium. Since we would like a steady state, the simulation must be run long enough that any transients excited at the inner surface have time to propagate radially outward, away from the inner surface where the dynamics is the most interesting. The UV-RMHD stage accomplishes this: We run the MHD solver with the UV RT module but not the IR RT module. The most expensive step of the UV RT module is ray casting, which is performed just once before the simulation starts, thus the amortized cost is low. In addition to being a preprocessing step for the IR-RMHD stage, the UV-RMHD stage is also valuable for studying dynamics driven exclusively by UV radiation. The UV RT module needs to be modified specifically for this stage. First, gas temperature in this stage is unsuitable for computing κ UV with Equation (14) because gas temperature in a gas-supported torus is virial, whereas gas temperature in a realistic IR-supported torus reflects the balance between radiative absorption and re-emission. Since UV radiative acceleration is strong only in the central hole, we estimate what the gas temperature there may be in the IR-RMHD stage, use that expressly for κ UV in the UV-RMHD stage, while keeping the actual gas temperature unchanged. As we shall justify in §4.7, such an opacity temperature would be this means κ UV is a function not of gas temperature, but of position. Second, gas should convert most of the energy it receives from UV radiation to the IR, but since it cannot do so without the IR RT module, Equation (11) overestimates the energy actually imparted to the gas. The correct value of S e UV should simply be the rate of work done by UV radiation, that is, S e UV ≡ v · S m UV . Note that although Wada (2012) also considered the case where UV radiation from the central source deposits momentum, not energy, their energy equation does not include a similar term. Third, we assume gas velocity vanishes in Equation (12), thereby ignoring the minute effect of Lorentz transformation. Fourth, the time step is arbitrarily multiplied by 0.25 to account for the fact that S m UV GMρ/r 2 at τ UV 1. We reset t to zero when the UV-RMHD stage begins; consequently, all times reported below are reckoned from the beginning of the stage.
2.3.3. Reduction of angular momentum in UV-RMHD stage MHD stresses establish j(R) in realistic tori. Simulating this process all the way to steady state is impossible in the IR-RMHD stage because of computational cost; it is also impossible in either the UV-RMHD stage or the IR-RMHD stage because MHD stresses redistribute angular momentum over tens of orbits, but our choice of κ UV implies that the radiation-driven outflow drains all the mass from an isolated torus in a couple orbits, and our simulations do not provide continuous mass resupply. Although the torus cannot reach a formal steady state in our simulations, it may nonetheless exhibit an approximate "steady" state wherein its inner edge stays close to the dust sublimation surface and its morphology is qualitatively the same over time. This "steady" state ends if UV radiation does enough positive work over time to gravitationally unbind the torus, but the amount of work required depends on j(R). We may obtain a relatively longlasting "steady" state by reducing j(R) of the MHD stage output before forwarding it to the UV-RMHD stage, but we must first determine what j(R) produces the longest "steady" state. Dorodnitsyn et al. (2011) likewise concluded that geometrically thick tori must have sub-Keplerian rotation, but on the basis of maintaining dynamical equilibrium in the spherically radial direction. Their argument can be recast more generally: Isotropic IR radiation pressure provides vertical and radial support simultaneously, hence a geometrically thick, IR-supported torus can only be in radial balance if rotation is sub-Keplerian. The same logic applies whenever geometrical thickness is ascribed to some isotropic pressure, be it gas pressure, radiation pressure, magnetic pressure, or velocity dispersion. Similar to ours, the latest simulations by Dorodnitsyn et al. (2016) employed an initial condition whose angular momentum profile is shallower than Keplerian.
Irradiation strengthens this argument. A geometrically thick torus, by definition, intercepts a sizable fraction of the radiation from the central source. Some of the radiative momentum absorbed by the torus may be carried away in an outflow, but whatever left behind constitutes a radially outward force, which could be strong enough to counteract gravity if L/L E τ T . Consider an unirradiated torus in which rotation precisely balances gravity; in other words, rotation is sub-Keplerian only to the extent isotropic pressures, if present, compel it to be. This torus cannot remain in equilibrium when the central source is turned on; to do so, its rotation must be even more sub-Keplerian. Moreover, radiation does positive work on outward-moving gas; if this energy is not advected away in its entirety by the outflow, the torus will increase in mechanical energy and be unbound eventually. The torus can stay in place only if it is replenished with mass to compensate for the outflow, and if this mass has sub-Keplerian rotation to offset the gain in mechanical energy.
For ease of parameterization, we multiply j(R) from the MHD stage by α(R/r ds ) −β , where 0 < α ≤ 1, β ≥ 0, and the dust sublimation radius is defined by (CK16) A parameter study decides the optimal values of α and β. While the computational cost of the UV-RMHD stage is merely a few percent of the IR-RMHD stage, extensive sampling of the parameter space is still unfeasibly expensive; we have therefore tested 11 pairs of parameters, as shown in the left panel of Figure 1. All runs have the same gravitational potential energy E grav at t = 0 but different kinetic energy E kin ; consequently, a useful parameter is b ≡ 1 + E kin /E grav , the binding energy normalized to the negative of the gravitational potential energy. Conscious effort is expended to ensure each run has a value of b similar to that of at least one other run. We choose τ T = 1 for the parameter study because this optical depth lies within the observed range (e.g., Risaliti et al. 1999), but its exact value is immaterial as long as the torus is optically thick to UV radiation. We shall demonstrate in §3.2 that α = 0.8 and β = 0.25 grant the torus the longest "steady" state. With parameters thus fixed, it remains to choose a snapshot of the UV-RMHD stage for passing on to the IR-RMHD stage. We impose two criteria on such snapshot. First, in the UV-RMHD stage, gas falls radially inward due to decreased rotational support, then rebounds and relaxes upon encountering the centrifugal barrier. The infalling and relaxing regions are separated by an outward-propagating shock; we consider a snapshot eligible only if the shock has moved sufficiently far away from the inner surface. Second, to ensure that the torus would not be blown away immediately in the IR-RMHD stage, we stipulate that the torus survive in the UV-RMHD stage beyond the selected snapshot for two more orbits at the inner edge. A more massive torus is able to withstand UV irradiation longer, allowing ourselves greater freedom in picking a snapshot that obeys both requirements; therefore, we conduct yet another run in the UV-RMHD stage with the optimal α and β, but with τ T = 2. The snapshot we opt for in the additional run is t = 50 t 0 .
This exercise provides us with a "steady"-state torus in the UV-RMHD stage; of course, there is no guarantee that it would remain so in the following IR-RMHD stage.

IR-RMHD stage
As we advance to the IR-RMHD stage, we reinstate IR radiative support by replacing gas pressure with a combination of gas and IR radiation pressure under thermal equilibrium; in other words, if gas temperatures before and after the replacement are T 1 and T 2 respectively, and the isotropic IR specific intensity in the fluid frame after the replacement is I 0 2 , then we demand ρR ideal T 1 = ρR ideal T 2 + 1 3 a SB T 4 2 and a SB T 4 2 = (4π/c)I 0 2 . The degree to which gas pressure is replaced is quantified by T 2 /T 1 = [1+ 1 3 a SB T 4 2 /(ρR ideal T 2 )] −1 . We clearly have R ideal T 2 ∼ c 2 s , where c s is the sound speed after the replacement. Since a SB T 4 2 ∼ ρv 2 φ , and since c 2 s v 2 φ if the torus is to be IR-supported and not gas-supported, we have T 2 /T 1 ∼ c 2 s /v 2 φ 1. Several comments are in order. First, the replacement preserves pressure, not energy or momentum, because we are interested in how the torus is supported. Second, the replacement does not promise exact force balance in the inertial frame; in fact, considering that the gas pressure tensor in the fluid frame is isotropic while the IR radiation pressure tensor in the same frame is ellipsoidal, there is no trivial transformation from one kind of pressure to another that would secure force balance everywhere. Third, the assumption that IR radiation pressure acts like a gas pressure is valid only in optically thick regions, so we are not justified to perform the replacement within the central hole; nevertheless, since the gas there is optically thin to the central source, the steady-state temperature profile quickly establishes itself no matter what the initial temperature is.
At the beginning of the stage, the radial Thomson and IR optical depths along the mid-plane are ≈ 1.59 and ≈ 23.1; the mass in the quarter-circle simulation domain ( §2.4) is ≈ 8.83 × 1 2 πρ 0 r 3 0 , compared to ≈ 2.10 × 1 2 πρ 0 r 3 0 for the initial condition of the RHD torus (CK16). We again reset t to zero. Now that the torus is IR-supported, we enable both IR and UV RT modules and study whether it can self-consistently stay so. The simulation is conducted in this stage to t = 14 t 0 . The final radial Thomson and IR optical depths along the mid-plane are ≈ 2.34 and ≈ 37.9; radiation-driven mass loss results in a final mass of ≈ 6.62 × 1 2 πρ 0 r 3 0 .

Simulation parameters and domain
The simulation domain spans [0.3 r 0 , 9.9 r 0 ] × [− 1 4 π, 1 4 π] × [−5 r 0 , 5 r 0 ] in (R, φ, z) in all stages. A large radial extent is needed to capture the extended flow after the radial expansion of the torus in the MHD stage ( §2.3.1); the geometrical thickness of the torus demands a similarly large vertical extent. We pick the number of grid cells to be 480 × 60 × 500 in (R, φ, z), large enough to resolve both MRI-driven turbulence (e.g., Hawley et al. 2013) and the UV absorption layer at the inner surface. The number of grid rays per cell is 168.
To make contact with our previous simulations (CK16), we let L UV /L E = 0.1 in both UV-RMHD and IR-RMHD stages. Numerical artifacts can appear if rotation does not exactly cancel gravity; such artifacts are smoothed out when c s /v φ O(0.1), where c s and v φ are the gas sound speed and orbital speed respectively at the inner edge (CK16). We set up our simulations such that the gas-radiation equilibrium temperature at the inner edge is always ≈ T ds , hence c s is independent of M. Since v 2 φ ∝ M/r ds ∝ M/L 1/2 UV , where we used r ds from Equation (20), we have this means a constraint on c s /v φ is also a constraint on M.
Here we choose c s /v φ = 0.1 as in CK16, corresponding to M ≈ 0.758 M . Its actual value in the IR-RMHD stage ranges from c s /v φ ≈ 0.15 at (R, z) = (r ds , 0) to c s /v φ ≈ 0.5 at (R, z) = (3.5 r 0 , 0) (see also §4.3); therefore, gas pressure remains a minor contributor to support where gas is most dense ( §4.2), and our results are not qualitatively affected by the small value of M. Quantitative changes may occur at larger M ( §5.3), but Equation (21) shows that the ratio of simulated to realistic value of M is much smaller than the analogous ratio for (c s /v φ ) −1 .

Scaling properties
Let us examine the scaling properties of Equations (2)-(5) in the three stages. In the MHD stage, the radiative source terms are zero; if we adopt a system of normalization in which v 2 0 = GM/r 0 , it is clear that the dimensionless equations are independent of M. The UV-RMHD stage introduces the UV radiative source terms in their modified forms ( §2.3.2). Since the (Table 1), and the normalization of S e UV is v 0 times that, the dimensionless equations remain independent of M. This means we are not committed to a particular value of M in either stage, and we may simply scale our results as needed to match any M.
The situation is very different in the IR-RMHD stage. The normalization of S m IR,UV is ρ 0 v 2 0 /r 0 as before, but now S m IR,UV has additional terms beyond zeroth order in v 0 /c; worse still, the normalization of S e IR,UV is now c/v 0 times that of S m IR,UV . The introduction of a fixed velocity scale c therefore breaks scalability in all equations except Equation (2) and the leading order of Equation (3). Because v 4 0 = GMκ T a SB T 4 ds , a choice of v 0 /c is equivalent to a choice of M, which we made in §2.4.

UV-driven dynamics
All runs except run A involve an initial suppression of angular momentum ( §2.3.3) and evolve in qualitatively similar fashion. We illustrate this general behavior with run G, displayed in the top row of Figure 2.
The torus quickly settles into a "steady" state ( §2.3.3). Because rotation alone provides insufficient support against gravity, gas falls radially inward at speeds ∼ v φ and converges toward the inner edge. After passing through a shock, the gas joins a lump at 0.6 r 0 R r 0 and |z| 0.3 r 0 , highlighted in Figure 3. UV radiation opens up the central hole as expected, but only weakly; all runs considered, the UV half-opening angle, defined as the angle between the axis and the τ UV = 1 surface, finds equilibrium somewhere between ∼ 0.24 rad and ∼ 0.38 rad.
Wind launching by UV radiation is bursty. When gas is shot out from the inner edge, it is accelerated along the τ UV = 1 surface by UV radiation, and at the same time pushed horizontally outward beyond the τ UV = 1 surface by centrifugal and radiative accelerations. Such gas excites a weak shock that propagates outward into the slower-moving gas at τ UV 1, a shock we call a palisade. Figure 3 has palisades at |z|/R 0.5 and τ UV 1; Figure 4 shows a schematic palisade. Palisades are found exclusively above a certain height because the lump stops outward motion at low latitudes.
The angle ψ between the palisades and the mid-plane is fixed by three parameters: the angle the τ UV = 1 surface makes with the mid-plane, χ ∼ 1.3 rad; the characteristic speed of UV-launched gas, v ∞ ∼ (GM/r ds ) 1/2 (L UV /L E ) 1/2 (κ UV /κ T ) 1/2 (CK16); and the shock propagation speed in the gas-supported torus, which is the sound speed c s ∼ (GM/r ds ) 1/2 . The angles in Figure 4 for our parameters ( § §2.2 and 2.4), we have c s /v ∞ ≈ 0.354, so the solution to the equation is ψ ≈ 1.6 rad. Palisades therefore look vertical in our simulations, but they need not be so in simulations with different parameters. Palisades are important because they represent a sizable portion of the outflow: They have v r > 0 even at τ UV 1, and most of the mass outflow through the vertical boundaries in the "steady" state is in fact at τ UV 1. Furthermore, their immediate adjacency to infalling gas implies inflow and outflow could interact, or even regulate each other. Figure 2 also shows (GMρ/r) −1 ( 1 2 ρv 2 + p), which is related to the Bernoulli constant; roughly speaking, gas with this parameter below or above unity is gravitationally bound to or unbound from the system respectively. The palisades are bound, and UV radiation cannot unbind them since they are at τ UV 1, so the outflow in the palisades cannot continue indefinitely outward.

Parameter study of reduction of angular momentum
One requirement for "steady" state is that the inner edge stays near the dust sublimation surface ( §2.3.3). We therefore define lifetime as the time when the intersection of the azimuthally averaged τ UV = 1 surface and the mid-plane passes outside  (20)) is the dotted black contour around the origin, and the red contour traces the surface on which τ UV = 1. All quantities are normalized to fiducial units (Table 1). Top row: Run G in the UV-RMHD stage at t = 40 t 0 . Bottom row: IR-RMHD stage at t = 14 t 0 . First and second columns: Gas density is presented on linear and logarithmic scales respectively as blue intensities (see color bars along the top edge). Third to fifth columns: Colors represent v R , v z , and (GMρ/r) −1 ( 1 2 ρv 2 + p) respectively (see color bars along the top edge).  R = r ds for the last time; we use the last time because the position of the inner edge fluctuates at the beginning and its initial crossings of the dust sublimation surface are of little import. It is unsurprising that lifetime increases with b, as evidenced by the right panel of Figure 1. Variation in α or β mostly creates scatter around this trend.
All runs have the same L UV , while the UV covering fraction C UV , defined as the solid angle around the central source with UV optical depth greater than unity, changes by at most a few percent over time and from run to run; therefore, the rate of UV momentum deposition also varies by a similar amount. We might expect lifetime to increase linearly with b, but that thought is not corroborated by the right panel of Figure 1. A plausible explanation is that the rate of work done by UV radiation is v · (ρκ UV F UV /c), so halting infall in fact raises binding energy at a rate proportional to −v r . A torus with lower j(R) has decreased radial support and faster inflow, hence UV radiation is less effective at unbinding it. Figure 4. Schematic diagram of a palisade in poloidal section. Gas is launched from the inner edge O along the inner surface OA at characteristic speed v ∞ ; as it travels, it excites a shock AB that propagates outward into the torus at speed c s . Palisades appear vertical because ψ ≈ 1 2 π rad for our parameters.
Lifetime diminishes for b 0.9 because such large binding energy is due to a j(R) so far below Keplerian that UV radiation cannot prevent a large portion of the gas from falling all the way through the inner-radial boundary of the simulation domain. Consequently, runs H and K cannot represent obscuring tori.
Of the remaining runs, run G stands out with the longest lifetime; we therefore select its parameters for the additional run that eventually goes on to the IR-RMHD stage ( §2.3.3). Incidentally, run G at t = 0 has approximately flat mid-plane j(R). Figure 5 shows that j(R) stays flat in the mean throughout the UV-RMHD stage, yet its increasing jaggedness suggests that the torus is moving away from a steady state, namely, that of the MHD stage. Because specific angular momentum is nearly homogeneous at t = 0, this must be due to angular momentum redistribution by either MHD stresses or non-axisymmetric pressure fluctuations. It is probable, but by no means certain, that the UV-RMHD stage has its own steady state, and the torus must pass through a disturbed state to reach it, but our simulations are not long enough for this to happen.
Our parameter study qualitatively corroborates our claim in §2.3.3 that a steady-state irradiated torus rotates more slowly than one unirradiated. The exact amount by which rotation is slower can depend on many factors, likely including L UV /L E , which we vary and which sets both the mass loss and the energy gain rates.

RESULTS OF IR-RMHD STAGE
Our principal results come from the IR-RMHD stage. The first part of this section deals with the initial transient phase ( §4.1), but the rest is dedicated to the "steady" state: the internal structure of the torus ( §4.2); the relative importance of different forces ( §4.3); the flow resulting from these forces ( §4.4); the character of the radiation-driven, high-latitude, wide-angle outflow ( §4.5); the distribution of IR radiation ( §4.6); the distribution of temperature ( §4.7); and the distribution of magnetic field ( §4.8). . Plot of specific axial angular momentum, weighted by density and averaged azimuthally and vertically over |z| ≤ 0.5 r 0 , against radial coordinate in run G in the UV-RMHD stage. The lines plot snapshots 5 t 0 apart, going from t = 0 at the bottom to t = 50 t 0 at the top. Upward shifts of 0.1 r 0 v 0 per line are added for legibility. All quantities are normalized to fiducial units (Table 1).

Transient behavior
The central hole opens up quite dramatically at the beginning of the IR-RMHD stage, in contrast to the UV-RMHD stage ( §3.1). The difference can be explained by considering radial force balance in the plane, say, z ∼ 3 r 0 , near where it intersects the τ UV = 1 surface. In the UV-RMHD stage, the locus of intersection fluctuates over time but is generally within 0.5 r 0 R r 0 . An approximate balance obtains in that stage at R 2 r 0 between radially outward centrifugal and UV radiative accelerations, inward gravitational acceleration, and largely inward gas pressure acceleration. Gravitational and UV radiative accelerations play minor roles here because their radial components, which are ∝ R/(R 2 + z 2 ) 3/2 , peak at |z|/R = √ 2 for any given |z|, but |z|/R √ 2 for the situation under discussion. Just before the IR-RMHD stage, most of the gas pressure is swapped for IR radiation pressure. An inward acceleration is accordingly removed, but it is not replaced because IR radiative flux diffuses outward from τ UV 1; centrifugal acceleration is no longer opposed, so the central hole widens. As the UV half-opening angle increases, the radial component of UV radiative acceleration at τ UV ∼ 1 becomes stronger, hence UV radiation begins to participate in expanding the central hole as well.
Analogous to previous simulations (CK16), we observe a strong chevron-shaped transient propagating radially outward from the inner surface as the central hole opens up. The transient can be regarded as dividing the torus into two parts: The gas outside has not fully responded to the change induced by turning on the IR RT module, whereas the gas inside has undergone at least partial relaxation and is evolving toward a "steady" state.

"Steady"-state behavior
The torus achieves a "steady" state at t 6 t 0 that lasts until the end of the simulation; at the average mass loss rate during the "steady" state, it would survive for ≈ 21 orbits. The much longer lifetime compared to CK16 is primarily because the current torus has 4 times the mass.
The bottom row of Figure 2 portrays the torus at a time when the transient has left the plotted volume entirely. Four structures with distinctive and persistent morphologies emerge; Figure 3 points out where they are. The structures are most recognizable in the first and second panels of Figure 2, but the other panels provide kinematic information that helps demarcate them. One is the UV-launched wind at τ UV 1 (CK16). Because the second panel resembles a top-down view of a bird flying toward the left, we name the other three by analogy with avian anatomy: The head refers to the very dense region enclosed by 0.6 r 0 R 1.4 r 0 and |z| 0.3 r 0 , the body is the somewhat less dense region at 1.3 r 0 R 4.2 r 0 and −0.6 r 0 z 0.4 r 0 , and the wings are the density ridges parallel to the τ UV = 1 surface at τ UV 1 and |z| 0.5 r 0 . We shall see in §4.4 that these structures are not hydrostatic, but merely parts of a global flow that retain their shapes throughout the simulation. The body does not lie entirely along the mid-plane because MHD turbulence breaks the symmetry about it. The head is denser than the lump in the UV-RMHD stage ( §3.1) by a factor of ∼ 2; although the head and the body are the densest parts of the torus, they take up 30 % of the total mass owing to their small volumes.
Although the torus is already in "steady" state, the wind is launched in bursts because density is not smoothly distributed at the inner edge of the head. This irregularity gives rise to the complex density structure at τ UV 1. Because the position of the τ UV = 1 surface is easily influenced by the presence of trace amounts of dusty gas in the central hole, the "steady"-state UV half-opening angle fluctuates between ≈ 0.61 rad and ≈ 0.69 rad at t 6 t 0 (see also §5.2). The average angle, ≈ 0.65 rad, is tantalizingly close to arctan(1/ √ 2) ≈ 0.62 rad. We speculate this is because the cylindrically radial component of UV radiative acceleration attains its maximum in any horizontal slice at |z|/R = √ 2.
The "steady"-state inner surface is corrugated in the azimuthal direction, similar to our simulations in CK16. Although the corrugation in CK16 increased rapidly in radial span, the one here remains bounded within 0.62 r 0 R 0.75 r 0 . We speculate that the corrugation grows only after the inner edge has started moving away from the dust sublimation surface. This speculation is supported by the rapid growth of the corrugation in every run of the UV-RMHD stage, but only during this outward recession phase, whereas the inner edge in the IR-RMHD stage never separates from the dust sublimation surface. We further hypothesize that if a realistic torus could maintain a true steady state as a result of mass resupply ( §2.3.3), its inner edge may stay close to the dust sublimation surface and hence not have strong corrugation.

Forces
The central unanswered question about obscuring tori is the nature of the forces supporting them vertically against gravity. Our simulation allows all relevant forces to be measured. We quantitatively define the advective, centrifugal, gravitational, gas, magnetic, IR, and UV forces as Note that f adv contains f cent . If we restrict ourselves to dynamics in the poloidal plane, the difference f adv − f cent can be understood in two equivalent ways: either as the rate of change of local momentum due to poloidal advection, or as the force f ram arising from ram pressure. The individual forces combine to form the Eulerian and "Lagrangian" non-gravitational forces: The Eulerian non-gravitational force is easy to grasp: If f E + f grav vanishes, then the flow is time-steady, but gas may still accelerate along streamlines. To interpret the "Lagrangian" nongravitational force, we consider the force equation for a gas packet: Since the R-and z-components of the left-hand side are ρ[dv R /dt − v φ (dφ/dt)] and ρ(dv z /dt) respectively, the same components of f L + f grav are ρ(dv R /dt) and ρ(dv z /dt). The "Lagrangian" non-gravitational force therefore tells us how the gas packet moves in the poloidal plane; the special case of f L + f grav having zero poloidal projection means that the gas packet does not move poloidally. Figure 6 compares non-gravitational against gravitational forces in radial and vertical directions. Only forces that explain support are included: The radial components of f gas and f mag are not shown because they, ignoring signs, are 0.1 and 0.05 times gravity in the body; similarly, the vertical component of f mag is omitted due to it being on average weaker than most other forces. Although the torus is asymmetric about the mid-plane throughout the simulation, it must on average be symmetric in the long run; therefore, we fold each quantity in the figure about the mid-plane by averaging the quantity with its vertical reflection. We then smooth out fluctuations by averaging over the interval 6 t 0 ≤ t ≤ 14 t 0 in which the torus is in "steady" state. The figure also shows density contours: The head and the body are the vertically extended and flat structures about the mid-plane at 0.5 r 0 R 1.8 r 0 and 1.8 r 0 R 3.9 r 0 respectively, whereas the wings refer to the region directly above them at τ UV 1.
Gravity is closely matched both radially and vertically by f E in a sizable region that encompasses the body and the densest parts of the head. This implies the velocity field in the region changes little over time, so structures there have time-steady morphologies, which is what we claimed in §4.2. All forces contribute to create the approximate equality between f E and gravity, with the relative significance of each force varying by position. The head and the body are inflated in the vertical direction by both gas and IR radiation pressure; the same can be said of the lower parts of the wings directly above the body, in the triangular region bounded by 2 r 0 R 3 r 0 , |z| 0.5 r 0 , and |z|/R 0.7. The head appears taller than the body because IR and UV radiation push the inner edge up into a wind. UV and IR radiation are the primary drivers of outflow in the wind and in the upper parts of the wings at τ UV 1 respectively. These are exactly the same observations we made about our previous simulations (CK16), but the tori in those simulations were never in a "steady" state long enough for concrete remarks to be made.
In the radial direction, f L only partially supports the body against gravity, so gas accelerates radially as it falls inward through the body. In the vertical direction, we interpret this as advection bringing gas with downward momentum to the body, or equivalently, as vertically collapsing gas exerting downward ram pressure on the body. The origin of this collapsing gas will be discussed in §4.4. Ram pressure squeezes the body vertically, making the vertical extent of the body smaller than its sound speed would otherwise suggest.
The vertical components of f adv , f gas , f mag , f IR , and f UV in the head, normalized by the negative of the vertical component of f grav , are ≈ −0.01, ≈ 0.42, ≈ 0.02, ≈ 0.44, and ≈ 0.18 respectively; analogous quantities in the body are ≈ −0.17, ≈ 0.67, ≈ 0.15, ≈ 0.35, and ≈ 0.00. Clearly vertical force balance prevails in the body, so gas motion is nearly horizontal; moreover, gas pressure is chiefly responsible for counteracting downward forces. We shall extrapolate these results to realistic AGNs in §5.3.

Streamlines
Figure 7 depicts streamlines in regions of "steady" state, where the flow timescale is shorter than the simulation duration. The global inflow-outflow is now manifest. The head, body, wind, and wings are not hydrostatic; rather, gas from, say, 3 r 0 R 4 r 0 and r 0 |z| 2 r 0 migrates first to the body, then to the head, and finally to the wind or the wings. The four structures are recognizable during the entire "steady" state because they retain qualitatively similar shapes ( §4.2) even as gas passes through them. Let us follow the marked streamline in Figure 7. Gas and IR radiation together furnish less radial and vertical support against gravity in the IR-RMHD stage than gas with artificially elevated pressure ( §2.3.1) did in the UV-RMHD stage; indeed, Figure 6 shows f L to be weaker than gravity in both directions at the starting point of the streamline. Weaker support, combined with the fact that gas was already infalling in the UV-RMHD stage ( §3.1), means that gas streams inward radially and vertically, as demonstrated by the bottom row of Figure 2 and by Figure 7. This gas piles onto the upper surface of the body, generating ram pressure ( §4.3). The reaction by the body on the gas results in a shock, visible in Figure 2, which removes the vertical component of velocity from the gas; this is why f L in Figure 6 points strongly upward throughout the body.  (20)) is the dotted black contour around the origin, the red contour traces the surface on which τ UV = 1, and gray contours plot ρ/ρ 0 from 10 −5 to 1 in logarithmic steps of 10 0.5 . All quantities are normalized to fiducial units (Table 1).  (20)) is the dotted black contour around the origin, red contours trace the surfaces on which τ UV = 1 and τ UV = 4, and light gray contours plot ρ/ρ 0 from 10 −5 to 1 in logarithmic steps of 10 0.5 . Gas in the body moves horizontally inward because f E nearly balances gravity in the vertical direction ( §4.3). This journey is arrested when gas reaches the irregular surface separating the head from the body. We witness in Figure 2 that a shock along this interface strips the gas of much of its remaining radial component of velocity. Flow in the head is much slower than in the body, but Figure 7 reveals that gas on the whole is still moving horizontally inward to the inner edge. As soon as it gets there, it is ejected by UV radiation.

Wide-angle outflow
The streamlines in Figure 7 suggest that while part of the gas leaving the inner edge becomes the wind at τ UV 1, part of it recedes to τ UV 1 and blends into the wings. Outflow occurs through both structures; what makes the wings special is that they stretch across a broader latitude range than the wind, and they also host a more massive outflow: The "steady"-state mass, momentum, and kinetic energy loss rates at τ UV > 1 are respectively ∼ 6.1, ∼ 1.8, and ∼ 0.85 times those at τ UV < 1. IR radiation diffusing outward from the central hole through the torus is vital for driving this wide-angle outflow. One piece of evidence supporting this claim is in Figure 6: The upper parts of the wings at τ UV 1 and |z|/R 0.7 experience outward IR radiative acceleration much stronger than gravity. Another is in the top panel of Figure 7, which plots the spherically radial component of gas momentum along the vertical boundaries of the simulation domain. The curve peaks at τ UV ≈ 4.87, close to τ UV = κ UV /κ IR = 4, the surface of unit IR optical depth from the origin.
It should be remembered that this division of the outflow into wind and wings is only for our cognitive convenience. The curve in the top panel of Figure 7 has no perceptible discontinuity at τ UV ∼ 1; IR and UV radiation work in tandem to power a continuous outflow across a large solid angle, from the wind at τ UV 1 to the wings at τ UV 1.
The upper parts of the wings in Figure 2 appear to be gravitationally unbound. Given our limited simulation domain, we cannot say definitively if the outflow in the wings would reach infinity. If the infalling part of the torus has a flaring shape, as it ostensibly does in Figure 2, and if gas velocity in the wings does not make a large enough angle with the mid-plane, then the outflow may eventually run into the inflow. The outflow may lose energy in shocks, become bound, and merge with the inflow, thus creating a circulation of gas in the torus.
Wings are the direct analogues of palisades in the UV-RMHD stage ( §3.1). Much that applies to palisades carries over to wings: Both are similarly located in the torus, both have gas flowing outward only above a certain latitude, and both contain small-scale density inhomogeneities caused by bursty wind launching and consequent shocks. The two structures are nevertheless very different dynamically: Palisades are not propelled by IR radiation, move much more slowly, and are gravitationally bound.
We saw in §3.1 how the angle between the palisades and the mid-plane can be estimated from just three parameters; here we carry out the same analysis mutatis mutandis. While the gassupported torus in the UV-RMHD stage had c 2 s ∼ GM/r ds ( §3.1), the torus in the IR-RMHD stage has c 2 s ∼ R ideal T ds GM/r ds . It is obvious from Figure 4 that reducing c s /v ∞ brings ψ closer to χ; in the limit of c s /v ∞ 0.1, as is applicable to our simulation, Equation (22)

IR radiation
We define the normalized IR radiative flux and radiation energy density as [L UV /(4πr 2 )] −1ê r · F IR and [L UV /(4πr 2 )] −1 (cE IR ) respectively; the former quantity is unity if all UV radiation were converted to the IR, and if IR radiative flux were spherically symmetric. Figure 8 plots the quantities measured on two surfaces: Blue curves are for a sphere of radius r = 3 r 0 cutting through the body; orange curves are for the outer-radial and vertical boundaries of the simulation domain. Our simulation domain is large enough for the normalized IR radiative flux to asymptote at large distances. The wiggles along the curves are the consequence of the IR RT module directing radiation into preferred directions in optically thin regions ( §2.2) and are not physical. The wiggles do not die out with distance, but fluctuations about the mean are at worst 20 %.
First consider the blue curves. The solid curve is 0.1 times the dashed curve at low latitudes; this means IR radiation is  Figure 9. Zoom-in of the azimuthally averaged and vertically folded poloidal plane in the IR-RMHD stage at t = 14 t 0 . Solid and dotted blue contours show gas and IR radiation temperatures; dash-dotted gray contours show the modeled temperature given by Equation (35). Temperatures on these contours go down from T ds in steps of 0.1 T ds as one moves away from the origin. Only parts of the IR and modeled temperature contours are drawn to reduce clutter. The upper and lower dashed red contours are the edge-on and face-on IR photospheres ( §5.1). Gas density is presented on a logarithmic scale as gray intensities (see color bar along the right edge), the dust sublimation surface r = r ds (Equation (20)) is the dotted black contour around the origin, and the solid red contour traces the surface on which τ UV = 1. All quantities are normalized to fiducial units (Table 1). fairly isotropic inside the body. In addition, the value of the solid curve at low latitudes is 0.05 times that at high latitudes; this is indicative of the degree to which the optically thick torus concentrates IR radiative flux into the polar direction (CK16).
Next we compare the solid curves. The curves almost coincide at latitudes 0.4; this suggestsê r · F IR ∝ r −2 in the central hole, the wind, and the wings, from as small a radius as r = 3r 0 outward. Both curves approach ∼ C IR /(1 − C IR ) near the axis, in agreement with our previous result (CK16); here C IR is the IR covering fraction, defined as the fraction of sightlines toward the central source with IR optical depth above unity. The blue curve falls below the orange at latitudes 0.4, but this is simply because the head and the body have limited radial extent, so outgoing IR radiation can diffuse around and reach regions behind them.

Temperature
Gas and IR radiation temperature contours in Figure 9 are very close to each other at τ UV 1. They are not spherical because the dashed curves in Figure 8 increase with latitude; in our particular simulation, the contours are strikingly vertical from the mid-plane almost up to the τ UV = 1 surface. The dashed orange curve in Figure 8 has greater variation over latitude than the blue; accordingly, contours further from the origin are less spherical.
Temperature in the central hole has a spherically symmetric distribution, with deviations only at the angles into which the IR RT module concentrates IR radiation ( §2.2). It can be modeled by considering the balance of energy emitted and absorbed by a dust grain of radius a in thermal equilibrium: Following CK16, the IR radiative flux is F IR ≈ C UV /(1 − C IR ) × L UV /(4πr 2 ). The UV absorption term contains the UV radiative flux corrected for extinction; such correction is unnecessary for the IR because τ UV 1 automatically implies τ IR 1, where τ IR is the IR optical depth from the central source. Rearranging, we get Using the operational definition of C IR and the simplification C UV ≈ C IR from CK16, we plot the modeled temperature in Figure 9; the model is excellent at τ UV 1. Equations (19) and (35) agree if τ UV 1 and F IR = 0, hence the opacity temperature in the UV-RMHD stage is consistent with the actual temperature in the IR-RMHD stage.

Magnetic field
Figure 10 graphs gas-only plasma beta β g ≡ p/( 1 2 B 2 ) and total plasma beta β t ≡ (p + 1 3 E IR )/( 1 2 B 2 ) at t = 0 and t = 14 t 0 . The range of β g varies little across space and time. The same can almost be said of β t , except that its value at t = 14 t 0 is much higher in the central hole and the upper parts of the wings than the rest of the simulation domain. This is because E IR is a few times higher in the central hole than in other parts of the torus ( §4.6), while p is several orders of magnitude lower; it then follows from the definitions of the plasma betas that β t β g . The spatial distributions of β g and β t are virtually identical apart from the overall normalization, and apart from the central hole. The temporal constancy of β g and β t means that complex gas motion does not perceptibly modify the MHD saturation state, at least not within our finite simulation time. Additionally, Figure 10 suggests that only the wings, where outflow drags out field lines, have large-scale order in the magnetic field, but even there the field does not point uniformly inward or outward, and neighboring regions can have fields in opposite directions.
We define the density-weighted plasma betas and alpha parameter as where ∆v φ is the perturbation of v φ about its azimuthal average. We see β g rise steadily from ∼ 4.2 to ∼ 4.6 over 0 t 10 t 0 and more rapidly to ∼ 7.5 at t = 14 t 0 , while β t drops from ∼ 27 to ∼ 19 over 0 t 2 t 0 and holds still at ∼ 16 thereafter. In addition, α climbs from ∼ 0.019 to ∼ 0.044 over 0 ≤ t ≤ 14 t 0 . The MHD simulations cited by Hawley et al. (2011) have β ∼ 15 and α ∼ 0.02, while the RMHD simulations by Hirose et al. (2009) have α ∼ 0.02; since our β t and α are close to these values, MHD turbulence is likely at saturation.
We see large spatial variation in plasma betas in Figure 10 because field fluctuations are of the same order as the mean. The associated velocity fluctuations, having Mach numbers ∼ 0.1, are much smaller than the global flow speed, so gas motion remains well ordered throughout the simulation, and the smooth time-averaged streamlines in Figure 7 closely resemble the flow at any one time.

DISCUSSION
In this section, we consider the observational implications of the IR-RMHD stage ( § §5.1 and 5.2), extrapolate results from that stage to parameters applicable to realistic tori ( §5.3), and compare our simulation with other torus models with outflows ( §5.4).

Observed temperature profiles in IR-RMHD stage
Interferometric observations of tori are analyzed by fitting ellipsoidal blackbodies of various sizes and temperatures to visibilities; the result can be interpreted as a crude temperature profile (e.g., Tristram et al. 2007). To facilitate comparison of the "steady"-state torus in the IR-RMHD stage with observations, it is useful to locate its IR photospheres as seen by observers both face-on and edge-on, and to determine its observed temperature profiles. Here we adopt the gray-opacity approximation of our simulation, that is, κ IR /κ T = 20 wherever T T ds ( §2.2). We ignore the possibility that rarefied gas at τ UV 1 could be much hotter than in our simulation because we do not treat its physics accurately, and because its Thomson optical depth along the vertical sightlines described below is 5 × 10 −3 .
The lower dashed red contour in Figure 9 shows the face-on IR photosphere of the torus, that is, where the IR optical depth as measured vertically from the upper-vertical boundary equals unity. The photosphere is approximately horizontal, separating the wings above from the head and the body below. Vertical sightlines at R r ds always intersect the photosphere; the photospheric temperature has a gradient ∝ R −0.73 starting from ≈ T ds at R ≈ r ds . The observed IR emission includes the contribution from all gas above the photosphere; however, gas at τ UV 1 does not significantly modify the temperature of the IR radiation from the photosphere because temperature contours inside this gas are close to vertical ( §4.7). According to Figure 9, if we were to enlarge the simulation domain vertically, the part of the photosphere at R 3 r 0 would barely change because sightlines from the upper-vertical boundary would simply encounter additional amounts of rarefied gas at τ UV 1, but the part at R 3 r 0 would shift upward noticeably because sightlines would cut through more of the denser gas at τ UV 1. This vertical displacement of the photosphere should not qualitatively change the observed temperature gradient since temperature contours are vertical ( §4.7).
The upper dashed red contour shows the edge-on IR photosphere, that is, where the IR optical depth as measured horizontally from the outer-radial boundary equals unity. Horizontal sightlines at |z| 2 r 0 intersect the vertical portion of the photosphere at R ∼ 5 r 0 ; the fact that this portion has temperature ∼ 0.2 T ds over its entirety may be germane to the observed constant temperature of ∼ 300 K at 2 pc scales in Circinus  (20)) is the dotted black contour around the origin, and the red contour traces the surface on which τ UV = 1. All quantities are normalized to fiducial units (Table 1).
intersect the photosphere where its tangent |dR/dz| is large, so the observed temperature on these sightlines rises rapidly with |z|; however, this part of the photosphere presents only a small projected area to edge-on views. Finally, horizontal sightlines at |z| 3.4 r 0 do not intersect the photosphere at all; the IR emission along these sightlines is likely weaker because they are optically thin, but hotter because they pass through gas with higher average temperature than sightlines at |z| 3.4 r 0 . If we were to expand the simulation domain radially, the photosphere would necessarily move outward, but the magnitude of the shift depends on the unknown gas distribution outside the simulation domain. Nevertheless, since density generally falls off with height, the presence of additional opacity at greater radii should only boost |dR/dz| of the photosphere, so we would still expect a sharp increase in observed temperature at some |z|.

X-ray, IR, and UV obscuration in IR-RMHD stage
The left panel of Figure 11 plots Thomson, IR, and UV optical depths against latitude for one snapshot in the IR-RMHD stage. To fit the curves into the same scale, IR and UV optical depths have been divided byκ IR /κ T = 20 andκ UV /κ T = 80 respectively ( §2.2); the actual optical depths can be discerned by means of the dotted lines marking the normalized optical depths at which sightlines become optically thick. The IR and UV curves are indistinguishable because κ IR and κ UV have identical temperature dependence; the two curves are very close to the Thomson curve because κ IR and κ UV depend weakly on temperature in most of the torus ( §2.2). The thinness of the shaded regions establishes that optical depth in our torus is largely independent of azimuth.
The right panel plots dΩ/d log 10 τ, the solid-angle coverage around the origin as a function of logarithmic normalized optical depths, during the entire "steady" state of the IR-RMHD stage. The similarity of the histograms is again due to the weak temperature dependence of κ T , κ IR , and κ UV . The peak at larger optical depth corresponds to sightlines through both the head and the body, while the peak at smaller optical depth is due to sightlines through the head only.
We define obscuration in X-rays by τ T = 0.01: Neutral gas of such column density is optically thick to X-rays at 2 keV. Likewise, we define obscuration in the IR and UV by τ IR = 1 and τ UV = 1 respectively. For our assumed mid-plane Thomson optical depth of ∼ 2 ( §2.3.4), the soft X-ray, IR, and UV covering fractions in the "steady" state of the IR-RMHD stage fluctuate within the ranges 0.78 C soft 0.83, 0.71 C IR 0.73, and 0.77 C UV 0.82 respectively; these ranges are close to one another, and also to the observed fraction of type-2 AGNs (e.g., Lawrence & Elvis 2010). Such broad coverage of the central source is achieved by combining the head, body, wings, and wind. In contrast, hard X-rays up to a few hundred keV are blocked along Compton-thick sightlines with N H 10 24 cm −2 or τ T ≥ 1; only sightlines traversing both head and body have such high optical depths, so they are concentrated near the midplane, taking up a solid angle of 0.15 C hard 0.28. This covering fraction coincides with the fraction of Compton-thick AGNs, observed to be ∼ 0.2 to ∼ 0.3 (e.g., Ricci et al. 2015;Koss et al. 2016).
Density profiles along sightlines to the central source serve as diagnostics of the statistics of density fluctuations, sometimes called "clumping" in this context. We graph these profiles in Fig-ure 12, leaving out regions not in a "steady" state ( §4.1). Apart from their overall scales, all profiles are similar in that density along a sightline decreases outward, but with fluctuations. These fluctuations are a blend of turbulent structures and, for latitudes 0.4 rad, density ridges in the wings ( §4.2); their density contrast goes from 50 % at low latitudes up to a factor of a few at high latitudes. These fluctuations, with their irregular shapes and small amplitudes, do not at all resemble the symmetrical, isolated clumps posited in many phenomenological RT models (e.g., Nenkova et al. 2002;Hönig et al. 2006;Nenkova et al. 2008;Schartmann et al. 2008;Heymann & Siebenmorgen 2012;Roth et al. 2012;Stalevski et al. 2012). In addition, on sightlines with latitudes 0.4 rad, the spacing between maxima stretches with latitude because the wings are parallel to the inner surface ( §4.5).

Extrapolating to realistic AGN tori
Our simulations cannot employ parameters that apply to realistic tori due to numerical reasons; instead, we must adopt a smaller central mass M, or equivalently, a higher c s /v φ (CK16), as well as a lower UV opacity κ UV ( §2.2). Moreover, our simulations consider just a few values for the Thomson optical depth τ T , and completely ignore photoionization and Compton heating. It is imperative that we understand how our results may change with c s /v φ , κ UV , and τ T , and how they may be modified by photoionization and Compton heating.

Mass loss rate and outflow speed
Mass, momentum, and kinetic energy loss rates are primary observables of our simulations. To extrapolate loss rates from our simulations to realistic tori, we introduced in CK16 a simple analytic model of an outflow of unit UV optical depth powered by UV radiation pressure. The fiducial mass loss rate and outflow speed of the model are given by Equations (25) and (26) in that article, which we duplicate below: We carry this model over to our current simulations. The model strictly pertains only to the UV-driven outflow in the UV-RMHD stage, but it can be adapted to describe the IR-driven part of the outflow in the IR-RMHD stage because that part also penetrates to unit IR optical depth ( §4.5). If we multiply L UV by ∼ C IR /(1 − C IR ) in Equations (39) and (40) to account for repeated scattering of IR radiation inside the central hole (CK16) and replace κ UV by κ IR , we finḋ The M-and κ UV -dependences of Equations (39)-(42) form the basis of our extrapolation. In the remainder of this discussion, we take R in = r ds with r ds from Equation (20), and set C IR ≈ 0.72 as found in the "steady" state of the IR-RMHD stage ( §5.2).  Figure 12. Plot of gas density in the IR-RMHD stage at t = 14 t 0 along sightlines to the central source at fixed azimuth φ = 0 and varying latitudes as indicated in the legend. The dashed and solid portions of each curve are at τ UV < 1 and τ UV > 1 respectively. All quantities are normalized to fiducial units (Table 1).
For the "steady" state of the UV-RMHD stage, the mass loss rate through the vertical boundaries is measured directly, while the outflow speed is obtained by comparing momentum and kinetic energy loss rates with the mass loss rate; the two values are then scaled by Equations (39) and (40). We obtain to their counterparts in that article, to wit, ∼ 0.1 M yr −1 and ∼ 5000 km s −1 . The temporal fluctuations of the two quantities in our current simulations are by ≈ 11 % and ≈ 26 % respectively. For the "steady" state of the IR-RMHD stage, we remove the contribution due to UV radiation by subtracting from each loss rate its average value in the UV-RMHD stage. The mass loss rate and outflow speed derived from these IR-only loss rates are then scaled by Equations (41) and (42). We obtain with fluctuations by ≈ 21 % and ≈ 16 % respectively. The IR-driven outflow has a greater mass loss rate but a smaller outflow speed than the UV-driven outflow; this is because the former has higher density and occupies a larger solid angle ( §4.5). Both IR-and UV-driven outflows fall within the observed ranges of mass loss rates and outflow speeds of UV and X-ray warm absorbers, which are respectively ∼ 10 −4 to ∼ 10 M yr −1 (e.g., Anderson & Kraft 1969;Crenshaw et al. 1999;Kaastra et al. 2000;Kaspi et al. 2000) and ∼ 100 to ∼ 2000 km s −1 (e.g., Blustin et al. 2005;Zhang et al. 2011;Crenshaw & Kraemer 2012).
In view of the torus not being in a formal steady state ( §4.2), it is remarkable that the mass loss rate and outflow speed of the combined IR-and UV-driven outflow in the "steady" state of the IR-RMHD stage fluctuate only by ≈ 11 % and ≈ 9 %. This constancy suggests that the properties of the outflow may be determined solely by the conditions at the inner edge, which does remain at approximately the same location throughout the "steady" state.

Morphology
The degree to which gas pressure lends support against gravity is ∼ (c s /v φ ) 2 , which according to Equation (21) is ∝ M −1/2 . In contrast, the typical IR radiation energy density in the torus is always ∼ a SB T 4 ds ; IR radiation pushes on the gas as it escapes the torus, resulting in a typical acceleration of ∼ a SB T 4 ds /(ρr ds ), which is independent of M insofar as ρr ds is constant. In other words, as M increases, gas pressure support weakens while IR radiation pressure support remains roughly the same.
Vertical support against gravity in the "steady" state of the IR-RMHD stage is due to different forces in different places: UV radiation powers the wind; IR radiation accelerates the wings ( §4.5); gas and IR radiation pressures contribute equally to inflate the head ( §4.3); the same pressures prop up the body, but gas pressure contributes roughly twice as much as IR radiation pressure ( §4.3). Based on the argument in the previous paragraph, we expect when M assumes realistic values, the wind and the wings would remain equally well supported against gravity, hence neither their dimensions nor their density fluctuations ( §5.2) would change substantially in character. The head and the body would lose respectively ≈ 1 2 and ≈ 2 3 of their vertical support and become geometrically thinner. Because gas and radiation influence the distribution of each other, there may not be a strict proportionality between c s /v φ and the aspect ratio of either head or body, but the greater importance of IR radiation in the head suggests that its height should scale more slowly with c s /v φ than that of the body.
The latitude range of the wind at τ UV 1 is likely governed by κ UV ; expressly, a higher κ UV implies a more slender wind along the inner surface (CK16).
The properties of all four structures may vary with the typical torus density, which is reflected in τ T . The density of the wind and the wings should be independent of density elsewhere; this is because the rate of momentum delivery to the torus by UV radiation, the characteristic size of the outflow, and the outflow speeds as given by Equations (40) and (42) do not depend on density. The defining attribute of the head is its vertical thickness relative to the body ( §4.2), maintained in part by IR and UV radiation ( §4.3); if the torus is denser, the penetration depth of radiation is reduced, so the head may appear less radially extended. The situation with the body is less certain due to positive feedback: An increase in density makes the body less permeable to IR radiation, weakening not only IR radiative support, but also gas support since gas temperature is tied to IR temperature ( §4.7); the torus becomes thinner and even denser as a result (see also Roth et al. 2012). The feedback loop may be limited by other physical effects that support the body vertically, such as Compton heating ( §5.3.6). Contrarily, an initial decrease in density triggers the feedback loop to run the opposite way, which causes the body to become less and less dense.

Temperature
Energy balance and temperature in realistic tori are controlled by external illumination; therefore, we expect IR radiation to diffuse outward from the inner edge, and temperature contours near the mid-plane to be not far from vertical, no matter what c s /v φ is. In contrast, if internal dissipation at the mid-plane were to dominate, IR radiation would diffuse from the mid-plane, which means temperature contours would make a sharp angle at the mid-plane. Indeed, temperature contours in our simulation are quite vertical far above the head and the body ( §4.7). We should not interpret this as a sign that temperature contours are strictly vertical at any c s /v φ , but as a suggestion that their shape enjoys relative independence from that of the head and the body, so they likely remain vertically extended even as c s /v φ is reduced.
Assuming such, we could ask how the observed face-on and edge-on temperature profiles ( §5.1) change as a function of c s /v φ . The face-on IR photosphere follows the outline of the head and the body. With smaller c s /v φ , both structures would be flatter; by virtue of the verticality of temperature contours near the mid-plane, a radially outward temperature gradient should always be observed. As for the temperature profile of the edge-on IR photosphere, the jump at large altitudes is a result of density in the wings diminishing with height; as long as wings are geometrically thick, the jump should remain.

Obscuration
The obscuration properties of the torus are a direct consequence of its density distribution. The soft X-ray, IR, and UV covering fractions should not vary strongly with c s /v φ or τ T ; this is because the primary obscurers are the wings and the wind ( §5.2), which should remain geometrically thick at lower c s /v φ and equally dense at any τ T ( §5.3.2). However, the Compton-thick fraction would decline with smaller c s /v φ because the obscurers here are the head and the body ( §5.2), both of which would be thinner at reduced c s /v φ ( §5.3.2).
Another important aspect of AGN obscuration is the distribution of observed column densities. Stated in terms of simulation variables, the observed dΩ/d log 10 τ is generally flat for τ T 0.01, with a slight rise toward higher τ T (e.g., Gilli et al. 2007, and references therein). Our torus already has flat dΩ/d log 10 τ for τ 0.4 ( §5.2), and we get an even better agreement with observations if we consider how dΩ/d log 10 τ would change under smaller c s /v φ or different τ T .
When c s /v φ is lowered, gas would gather toward the midplane, and both head and body would take up a smaller solid angle around the origin; as a consequence, the slender peaks in the right panel of Figure 11, which reflect the densest parts of the torus, would move to the right, and the area under them would diminish. There may still be the same gentle roll-off to a plateau on the left side due to the wings and the wind, but the plateau would be higher because the total area under the histogram is conserved. The histogram would therefore become flatter overall, which means dΩ/d log 10 τ would be practically constant over a large span of log 10 τ.
Alternatively, realistic tori may have a wide range of τ T extending up to a few. Since our torus already has τ T ∼ 2 ( §2.3.4), we are mostly interested in cases of smaller τ T , which have leftshifted peaks. When dΩ/d log 10 τ histograms for tori of different τ T are stacked, offset peaks would add up to form a flatter distribution; in other words, the near constancy of dΩ/d log 10 τ over log 10 τ could be a natural consequence of realistic tori having a broad distribution of τ T .

Photoionization
Our simulation is best at treating IR RT at τ IR 1, or where the column density from the central source is 8 × 10 22 cm −2 ; it is less reliable inside the central hole since we omit UV and X-ray photoionization. Here we estimate how reinstituting these effects may modify our torus.
Photoionization raises gas temperature in the central hole (Krolik & Begelman 1986;Krolik & Kriss 2001), and the higher gas pressure could compress the torus Dorodnitsyn et al. 2016). Although our simulation does not explicitly constrain pressure from photoionized gas in the central hole, we can estimate its effect by computing the ionization parameter where L ion is the ionizing luminosity, and p t is the sum of gas and IR radiation pressure (Krolik et al. 1981). AGN photoionization calculations indicate that Ξ locks in at ∼ 10 if, as here, a cool gas reservoir is present (Krolik & Kriss 2001). Taking L ion ∼ L UV , we find that Ξ in the head, body, wings, and wind are ∼ 1, ∼ 2, ∼ 2.5, and ∼ 1.2 respectively. This means our torus is too pressurized to be confined by photoionized gas; rather, its geometrical thickness is limited primarily by gravity.
Hotter gas can also destroy dust grains through sputtering (Draine & Salpeter 1979;Tielens et al. 1994), thus reducing the effectiveness of wind driving through radiation pressure on dust. The sputtering timescale at gas temperature of ∼ 10 5 K is ∼ 100 [n H /(10 6 cm −3 )] −1 [a/(0.1 µm)] yr, where n H and a are the number density of hydrogen atoms and the radius of dust grains respectively (Draine & Salpeter 1979;Tielens et al. 1994), and it decreases sharply at higher temperatures until ∼ 10 7 to ∼ 10 7.5 K. This is not very different from the time it takes for the UV-driven wind to escape from the inner edge at R = R in to infinity: The dust sublimation radius r ds is given by Equation (20), and the wind speed is defined in CK16. The dust content of the outflow therefore depends sensitively on how quickly cold gas evaporated from the inner surface rises in temperature. Photoionization could impart momentum to the wind through lines, similar to the accretion disk wind of Proga et al. (2000) and the radiation-accelerated magnetocentrifugal wind of Everett (2005). Heating by photoionization can also produce a thermally driven wind (Krolik & Begelman 1986;Balsara & Krolik 1993;Krolik & Kriss 2001;Blustin et al. 2005). This suggests some kind of complementarity between outflow-driving mechanisms: In places where Ξ 1, gas temperature is 10 4.5 K (Krolik et al. 1981), so gas may stay dusty long enough for radiation pressure to evacuate it from the AGN; in places where Ξ 1, gas temperature is near Compton (Krolik et al. 1981), so gas may already have enough thermal pressure to expel itself.

Compton heating
We ignore volumetric Compton heating in our simulation, but X-rays at luminosities a few times weaker than in the UV can produce sufficient IR radiation deep inside the torus to induce important changes to the distributions of density and specific angular momentum, particularly near the mid-plane (Shi & Krolik 2008). As a matter of fact, since the mid-plane is optically thick in the IR but only marginally so in X-rays, X-rays could be the deciding factor in how concentrated gas is near the mid-plane. A positive feedback loop could exist, wherein Compton heating injects the energy requisite to seed a vertical IR radiative flux, which lowers the density near the mid-plane and allows IR radiation from the inner edge to enter the torus and support it (see also Roth et al. 2012).
Compton heating most certainly changes gas and IR temperatures at τ UV 1. If the torus is Compton thin in all directions and has approximately the same density everywhere, we would expect Compton heating to be uniform, hence temperature contours would be vertical near the mid-plane as if X-rays were absent ( §5.3.3). But if the torus transitions from Compton thick to Compton thin with increasing latitude, then Compton heating would occur predominantly near the mid-plane, at a rate diminishing with distance from the central source and with effect resembling internal dissipation ( §5.3.3); furthermore, if sufficient energy is deposited at the mid-plane that the secondary IR radiative flux has noticeable effect on vertical support, temperature contours may no longer be vertical near the mid-plane.

Relation of our model to literature
This article does not report the first attempt at understanding torus dynamics. On the topic of outflows alone, several authors, we included, have written extensively on how outflows could account for torus phenomenology. We situate the present work relative to past ones by summarizing the contrasts.

Comparison with our previous work
The approaches taken in CK16 and here for the UV-RMHD stage ( §2.3.3) are quite complementary. In the earlier simulations, we fixed a sub-Keplerian j(R) and searched for L UV such that the torus is reasonably long-lived; we concluded that L UV should be small to allow under-supported gas to flow toward the inner edge, but not so small that the gas contracts into the dust sublimation surface. In our current simulations, we choose L UV and determine j(R) that best promotes torus longevity ( §3.2); we discover that the torus can survive longer if j(R) is sub-Keplerian and gas feeds the inner edge, but there is a limit to how far below Keplerian we can go before the torus shrinks to smaller than the dust sublimation surface. Königl & Kartje (1994) proposed that obscuration in AGNs could be provided by a magnetocentrifugal wind, which, intuitively speaking, is a centrifugally driven outflow guided by open magnetic field lines. Their model presupposes a razor-thin accretion disk as the mass source for the wind; our simulation allows a given amount of gas to evolve without prescribing how mass resupply occurs. Their model also requires the specification of several parameters for the wind, among which are the launch radius from the accretion disk, the spherically radial density profile, and the conserved mass flux and angular momentum along the streamline, whereas our simulation fixes only the initial angular momentum profile ( §2.3.3). Most importantly, their model posits the existence of a large-scale, dynamically dominant magnetic field; in contrast, we assume only that the MRI amplifies magnetic field contained entirely within the gas. We further find that radiation alone can lift gas from the inner edge into a high-latitude outflow ( §4.5). Despite our torus being magnetized, this outflow is not a magnetocentrifugal wind because meandering loops of magnetic field in the outflow are too weak to exert much force; instead, they are passively dragged out by gas motion ( §4.8). Königl & Kartje (1994) acknowledged the importance of IR and UV radiation pressure on dust, but Everett (2005) was the first to study how an AGN-like spectrum photoionizes a magnetocentrifugal wind and injects momentum through atomic absorption. He put shielding gas of arbitrary column density into the model to prevent overionization of the wind, but omitted dust, which could also influence photoionization by removing UV radiation. Keating et al. (2012) expanded the work of Everett (2005) by adding dust opacity and momentum transfer from dust absorption. Neither model considers pressure on dust from reprocessed IR radiation, which transmits momentum and energy deposited by UV radiation through the torus, and which our simulation identifies as critical for maintaining an outflow at τ UV 1 ( §4.5).

Comparison with magnetocentrifugal wind models
Elitzur & Shlosman (2006; see also Kartje et al. 1999) examined another variation on the magnetocentrifugal wind, one in which dusty gas is clumped; the authors assumed clumps are individual entities and touched on how they may be magnetically confined. Our simulation sheds light on both issues. The first and second panels of the bottom row of Figure 2 depict the inhomogeneous density distribution in the wings and the wind, which obscure at high latitudes ( §5.2). We find wedges of various densities in the wind, density ridges in the wings, and even a hook-shaped feature near the top of the second panel. The density profiles in Figure 12 illustrate the same point in a different way. Our simulation therefore underlines the point that the common picture of spherical, well-separated clumps must not be taken too seriously. Furthermore, density perturbations in our torus are not static or stationary structures confined externally by gas or magnetic pressure, or internally by self-gravity; they are imprinted in the wind by bursty wind launching ( §4.2), and in the wings by bursts of newly launched, faster gas shocking with slower gas further out ( §3.1). Density perturbations are ephemeral; only by virtue of their frequent recurrence at the same place with the same morphologies do they become consistent features of the torus. Wada (2012;2015; see also Schartmann et al. 2014;Wada et al. 2016) explored through simulations the idea of a gas fountain powered by the central source through UV radiative acceleration, X-ray photoionization heating, and Compton heating. The neglect of heating by UV radiation from the central source in these simulations precludes the treatment of thermal IR radiation. When gas in these simulations moves out of the central hole, it receives reduced radiative acceleration and falls back to the mid-plane; this could be because these simulations presume UV radiation is concentrated in the polar direction, and because they ignore IR radiation, which transports momentum and energy to τ UV 1 and thus extends the outflow into that region ( §4.5). Wada et al. (2016) combined in one simulation turbulence generated by supernovae (Wada & Norman 2002) with driving by the central source, but stars cannot make gas geometrically thick on parsec scales (Krolik & Begelman 1988  : Their simulations have IR radiation driving a wide-angle outflow, while ours have IR and UV working in concert to achieve the same effect ( §4.5). However, there are also important differences between the two sets of simulations. Their simulations postulate a razor-thin accretion disk as a mass source for the outflow, the characterization of which introduces additional free parameters not self-consistently determined by the simulation; we avoid this by putting all the mass in the simulation domain right at the start and letting it develop structures on its own. The outflow in their simulations originates from across the entire accretion disk because their mid-plane boundary condition assures so; in our simulation, where no such boundary condition is assumed, gas is launched into the outflow exclusively at the torus inner edge ( §4.2). Finally, the IR-driven outflow in their simulations is a failed wind that apparently falls back to the mid-plane. Our simulation with both IR and UV radiation tells a different story: The UV-driven outflow at τ UV 1 is gravitationally unbound, and at least part of the IR-driven outflow at τ UV 1 is unbound ( §4.5).

Comparison with photoionization-driven models
The later simulations by  and Dorodnitsyn et al. (2016) consider how the central source deposits momentum through UV radiation, and momentum and energy through X-rays; they diverge markedly from the simulations by  and from ours. In both their simulations and ours, IR radiation is created when the central source heats dust; the difference is that in our simulation, heating through UV radiation is concentrated at the inner edge (CK16), whereas in their simulations, heating is more widespread because it is due to X-rays, not UV radiation. In their simulations, UV radiation transfers momentum solely through lines; this implies that the authors were looking at a different situation from ours, namely, one with a dust-free central hole, but whether dust is present depends on how fast photoionization raises gas temperature to Compton ( §5.3.5). Lastly, photoionization heating in their simulations produces a hot atmosphere that envelopes and vertically squeezes their tori; therefore, IR radiation can only push an outflow along, not far above, the mid-plane. Our torus, in comparison, has high enough gas and IR radiation pressure that photoionized gas is not the main determinant of its geometrical thickness ( §5.3.5).
The recent simulations by Namekata & Umemura (2016) study how dusty gas interacts with UV radiation and X-rays from the central source, and with reprocessed IR radiation; dust absorption, Compton recoil, and a variety of chemical processes transfer energy and momentum between gas and radiation. Their simulations always produce geometrically thin structures and thus do not explain vertical support in tori. There are several reasons that might explain this thinness. Their initial condition is geometrically thin and therefore has Keplerian rotation, but Keplerian rotation is inconsistent with a long-lived, geometrically thick structure ( §2.3.3). Their central source radiates zero flux in the mid-plane, so little radiation enters the gas; in addition, their geometrically thin structure is poorly resolved in the vertical direction at small radii. Lastly, a photoionized and Compton-heated atmosphere surrounds their geometrically thin structure and confines the cooler gas to the mid-plane.

Summary of comparisons
In sum, our simulation offers a new perspective on torus dynamics. With the bare minimum of physics, to wit, momentum and energy coupling between gas, IR radiation, and UV radiation, our simulation demonstrates that radiation on its own can propel an outflow far above the mid-plane; a mass-loading mechanism and a strong magnetic field steering gas to high latitudes, as in a magnetocentrifugal wind, are unnecessary ( §4.5). Moreover, our simulation distinguishes itself from gas fountain models by showing that IR radiation is pivotal in delivering momentum and energy to τ UV 1, thereby driving a wide-angle outflow ( §4.5). Although photoionization can augment the outflow in the central hole ( §5.3.5) and Compton heating can modify the shape of structures deep inside the torus ( §5.3.6), neither is likely to change the fundamental character of the radiation-driven outflow.
Our simulations in CK16 and here also highlight two aspects of tori that have not always received the attention they deserve. First, efforts to understand the observed IR spectrum of AGNs often assume torus gas is in clumps, and these clumps are almost invariably taken to be spherical, discrete, and pressure-confined (e.g., Elitzur & Shlosman 2006). In contrast, our current simulation suggests that fleeting, irregular density perturbations can arise simply from cool gas accelerating to high speeds due to radiation pressure, then shocking with slower gas ( §4.5). Detailed RT calculations will be needed to assess whether such density inhomogeneities produce strong enough far-infrared emission and a shallow enough silicate feature to match observations. Second, our previous and current simulations both portray the torus as a flow-through system: Gas is conveyed inward along the mid-plane; the majority of this gas is expelled in the outflow, and only a small fraction is captured by the central mass to fuel the generation of UV radiation. This view has already been taken by Krolik & Begelman (1988), but our work makes clear that sub-Keplerian rotation is necessary for maintaining a steady-state torus in the presence of strong radiation pressure ( §2.3.3), and we have provided an estimate of the mass resupply rate requisite for steady state ( §5.3.1).

CONCLUSIONS
We have performed three-dimensional, time-dependent RMHD simulations of AGN tori featuring quality RT and simultaneous evolution of gas and radiation. For the first time, our torus achieves a "steady" state lasting for more than an orbit at the inner edge, and potentially for much longer. This "steady" state is defined as the torus having constant overall morphology ( §4.2). It is obtained by reducing the angular momentum profile before the simulation starts ( §2.3.3), which raises the total binding energy and thus allows the torus to survive longer under UV radiation doing positive work ( §3.2).
The existence of a "steady" state is significant: While tori in previous simulations could not endure UV irradiation for more than two orbits at the inner edge (CK16), our current simulation demonstrates that a torus with the right parameters can indeed remain in a steady state for multiple orbits ( §3.2). Granted that our torus cannot formally reach equilibrium owing to our choice of κ UV (CK16) and to computational cost, we can already learn much from its approximate "steady" state that would conceivably carry over to the true steady state. Moreover, the ability to study the torus in a quasi-stationary state boosts our confidence in separating the "steady" state ( §4.2) from transitory behavior ( §4.1).
We perceive four "steady"-state structures in the torus, namely, head, body, wings, and wind ( §4.2). Vertical support against gravity is dominated by gas pressure in the head and the body, IR radiation pressure in the wings, and UV radiation pressure in the wind ( §4.3). By inspecting the "steady"-state flow and the forces driving it, we realize that these structures are not hydrostatic. Instead, due to insufficient support against gravity, most gas falls toward the inner edge through the body and the head; as gas reaches the inner edge, it flies outward on UV radiative acceleration ( §4.4). The outflow is initially directed at high latitudes, but it spreads out in solid angle once it climbs above the head and the body ( §4.5). The part remaining in the central hole is the wind; it is propelled by IR and UV radiation. The part expanding beyond the central hole is the wings; it is powered by IR radiation ( §4.5). The four structures are simply regions that hold on to their shapes as gas flows through them.
The study of forces clarifies the subtle role IR radiation plays in torus dynamics: It opens up the central hole ( §4.1), partially supports the body and the lower parts of the wings in the vertical direction ( §4.3), and drives an outflow in the wings where UV radiative acceleration fails ( § §4.5 and 5.3.1).
It is reassuring that most statements pertaining to RHD tori are valid here as well: Gas and IR radiation have equal temperature inside the optically thick torus ( §4.7), the torus focuses IR radiation toward the axis ( §4.6), and the outflow has mass loss rate and speed consistent with observations ( §5.3.1). The strong resemblance between RHD and RMHD tori suggests that the effect of magnetic field over timescales as short as a few orbits is small ( §4.8). It bears reiterating that the influence of magnetic field on realistic tori is felt only over many orbits, as MHD stresses redistribute angular momentum and thereby set the steady-state angular momentum profile.
Observational predictions can be more easily made for a torus in "steady" state. When seen face-on, the temperature profile of our torus should follow the radially outward temperature gradient of the body, which is T ∝ R −0.73 in our simulation; when seen edge-on, a jump in temperature should be seen at high altitudes ( §5.1). The tenuous wings and wind obscure the central source in soft X-rays, the IR, and the UV, while the dense head and body also stop hard X-rays. The soft X-ray, IR, and UV covering fractions are all approximately three quarters, which is close to the observed fraction of type-2 AGNs ( §5.2). Furthermore, if we assume AGNs have a finite range of mid-plane column densities, then our torus also naturally explains why the distribution of observed AGN gas columns over logarithmic column density is flat ( §5.3.4).
The torus around a given central mass M is governed by three important parameters: the Eddington ratio L UV /L E of the central source, the Thomson optical depth τ T of the torus, and the angular momentum profile j(R) of the same (CK16). The first two parameters are fixed in our simulations by observational constraints ( § §2.3.3 and 2.4); the only freedom we have is with j(R) ( §2.3.3). Yet, with practically no fine-tuning, our torus naturally arrives at a "steady" state typified by a high-latitude, wide-angle outflow whose obscuration properties agree fairly well with observations. Such outflow therefore deserves serious consideration as a model for geometrically thick obscuration in AGNs.