DISSIPATION AND VERTICAL ENERGY TRANSPORT IN RADIATION-DOMINATED ACCRETION DISKS

Much of the power radiated by the accretion flow in luminous states of X-ray binaries and active galactic nuclei necessarily originates from the vicinity of the central compact object where the gravity well is deepest. If the bolometric luminosity is anywhere close to Eddington, radiation pressure must dominate the thermal pressure of the accreting plasma in these regions, and understanding the physics of radiation-dominated accretion is therefore central to any explanation of how these sources work in their brightest states. Classical accretion disk theory (Shakura & Sunyaev 1973; Novikov & Thorne 1973) models the flow as geometrically thin and optically thick, and assumes that angular momentum is transported within the flow by internal stresses due to turbulence. The average internal stress is postulated to scale with thermal, and therefore mostly radiation, pressure. Accretion power is assumed to be locally dissipated as heat, which is then transported vertically outward by photon diffusion.

This model has been questioned over the years on a number of grounds. First and foremost, if one takes its assumptions literally, then the resulting equilibrium structure is unstable to both thermal and inflow ("viscous") instabilities (Lightman & Eardley 1974; Shibazaki & Hōshi 1975; Syunyaev & Shakura 1975; Shakura & Sunyaev 1976). However, these assumptions may not be valid. Bisnovatyi-Kogan & Blinnikov (1977) point out that if the turbulent dissipation rate per unit mass is constant with height in the radiation pressure-dominated regime, then the mass density would be independent of height inside the photospheres, a situation that would clearly be convectively unstable.⁵ The argument is simple and worth repeating. Hydrostatic equilibrium in a radiation-dominated disk implies that the vertical profile of radiation flux F(z) with height z obeys κF/c = Ω²z, where κ is the opacity (which is very nearly constant as Thomson scattering dominates the flux mean), c is the speed of light, Ω is the local angular velocity, and we have assumed Newtonian gravity for simplicity. If radiative diffusion dominates the heat transport, then equilibrium requires that the flux divergence dF/dz be equal to the dissipation rate per unit volume Q, giving Q = cΩ²/κ ≡ Q^⋆, which is constant.⁶ If the dissipation rate per unit mass is also assumed to be constant, then the density must be constant, implying convective instability. It should be noted, however, that vertical energy transport by convection under the conditions of constant dissipation per unit mass need not alter the thermal instability of the disk, as diffusive transport of energy could still be dominant (Shakura et al. 1978).

Energy may also be transported vertically in other ways. Rather than completely dissipating locally, magnetic energy in the turbulence could be transported outward by buoyancy in the form of Poynting flux, perhaps to be dissipated outside the disk photosphere in a corona (Galeev et al. 1979). It has even been argued that magnetic buoyancy would limit the magnetic energy density in a radiation-dominated plasma to be at most the gas pressure, thereby producing an accretion stress that scales with the gas pressure alone (Sakimoto & Coroniti 1981, 1989; Stella & Rosner 1984). Such a stress would eliminate the thermal and inflow instabilities that plague the standard model.

Unfortunately, this work was built on an incomplete foundation because it lacked an understanding of the physical nature of angular momentum transport in accretion disks. One could therefore only guess the vertical profile of dissipation per unit mass and pretend that the disk can be adequately modeled by a time-averaged, steady-state vertical structure. Moreover, arguments concerning the efficacy of magnetic buoyancy had to make assumptions about the rate of magnetic field generation, as well as the geometry of magnetic field lines.

Since the discovery of the relevance of magnetorotational (MRI) turbulence to accretion disks (Balbus & Hawley 1991, 1998; Hawley & Balbus 1991), it has become possible to explore these ideas in detail with numerical simulations. The most relevant simulations thus far have used a shearing box that incorporates the vertical tidal gravitational field of the central object (Brandenburg et al. 1995; Stone et al. 1996). The accretion stress in such simulations arises self-consistently from correlated magnetic and velocity fluctuations within the turbulence itself. Using such stratified shearing boxes with an isothermal equation of state, Miller & Stone (2000) found that the majority of the magnetic energy generated in MRI turbulence was (numerically) dissipated locally. Nevertheless, a quarter of the magnetic energy generated was vertically transported outward by buoyancy, forming a strongly magnetized corona outside a weakly magnetized structure near the disk midplane. Miller & Stone (2000) used an isothermal equation of state, however, and did not include the possibility of diffusive radiation transport outward along the vertical temperature gradients.

Inclusion of such transport would also introduce new dissipation physics, as compressible waves should be damped by radiative diffusion (Agol & Krolik 1998). This process is entirely analogous to the Silk damping of acoustic perturbations in the early universe (Silk 1967, 1968) and to the radiative damping of stellar pulsation modes (Cox 1980). It is particularly interesting because such dissipation can easily be resolved numerically in radiation MHD simulations. This contrasts sharply with our complete inability to resolve the microscopic scales on which viscous and resistive dissipation damps fluid and magnetic fluctuations.

Shearing-box simulations of MRI turbulence incorporating thermodynamics and radiation transport have been possible for some time now. The first such simulations neglected vertical gravity and explored the properties of the turbulence in the presence of radiation transport, treated numerically using flux-limited diffusion (Turner et al. 2002, 2003). If photon diffusion is rapid enough, and the magnetic pressure exceeds the pressure in the gas alone, then MRI turbulence can become extremely compressible with strong density fluctuations. These fluctuations are highly dissipative: net PdV work is done on the plasma over time, indicating an irreversible conversion of mechanical energy into internal energy (Turner et al. 2002, 2003). This result confirmed the suggestion of Agol & Krolik (1998), although the specific character of the fluctuations being damped was not entirely clear.

The first radiation MHD shearing-box simulation of MRI turbulence with vertical gravity was published by Turner (2004). The simulation was radiation pressure dominated and exhibited no evidence of thermal instability. Moreover, the time-averaged vertical entropy profile was stable to hydrodynamic convection. On the other hand, the simulation was not fully energy-conserving, nor was it able to incorporate the photospheres within the simulation domain. Substantial mass loss occurred during the simulation and this was suggested as a possible cause of the absence of any exponential thermal runaway. An attempt was also made to measure the radiative damping using the time-averaged PdV work, but the result was necessarily uncertain because this work can also be used to excite vertical mechanical motions when vertical gravity is present.

The technical issues with the Turner (2004) simulation were solved by Hirose et al. (2006). Radiation transfer was still treated by flux-limited diffusion, but a new diffusion solver permitted several improvements. The quasi-periodic radial boundary conditions appropriate to shearing boxes could be imposed properly. Low densities could now be handled, allowing the photospheres to be incorporated within the simulation domain. Most importantly, a total energy scheme was implemented so that grid-scale losses of magnetic and kinetic energy were fully captured as heat in the gas, and the code accurately conserves energy. While dissipation is therefore still numerical, it may nonetheless mimic local energy flow to microscopic dissipation scales at high wavenumbers in the turbulent cascade. Radiation and gas exchange momentum through Thomson scattering and free–free opacity, and exchange energy through free–free absorption and emission. Energy exchange by Compton scattering was later incorporated into the code in Hirose et al. (2009b).

Using this improved code, we have succeeded in establishing long-lived thermal equilibria in which gas pressure dominates (Hirose et al. 2006), gas and radiation pressure are comparable (Krolik et al. 2007), and radiation pressure dominates (Hirose et al. 2009b). The radiation-dominated state is thermally stable, even with no mass loss, and this is due to the fact that thermal pressure lags turbulent stress on timescales of order the thermal time (Hirose et al. 2009b). We have also established that the stress and total pressure are correlated (see also Ohsuga et al. 2009), but only on timescales longer than the thermal time, reviving the question of inflow instability (Lightman & Eardley 1974; Hirose et al. 2009a). Neither the gas density nor the dissipation rate per unit mass is constant with height in the radiation-dominated (or any other) regime (Hirose et al. 2009b), and the time- and horizontally averaged structures are stable to hydrodynamic convection (Krolik et al. 2007). However, the outer layers are generally magnetically supported (Hirose et al. 2006) and exhibit Parker instability dynamics (Blaes et al. 2007).

There remain a number of important questions about the thermodynamics of radiation-dominated accretion disks. How does one properly calculate the contribution of radiative damping to the overall dissipation when PdV work can also be used to excite vertical mechanical motions? What exactly are the compressive motions that are being dissipated by radiative diffusion? What happens if the local dissipation rate per unit volume exceeds the radiation pressure-dominated hydrostatic value of Q^⋆ = cΩ²/κ? What controls the relative shares of radiative diffusion, advection, and Poynting flux in vertical energy transport, and how does this depend on radiation to gas pressure ratio and height? Can there be significant coherent energy flux in the form of, e.g., vertical acoustic waves, and if so, how are these waves excited and how much energy do they transport and dissipate? Even if the time- and horizontally averaged vertical structure is convectively stable, can local and transient buoyancy lead to significant energy transport?

The goal of this paper is to answer these questions on the basis of detailed analysis of simulation data. In Section 2, we provide a brief overview of the radiation-dominated simulations we have analyzed. In Section 3, we analyze global energetics, first deriving and discussing the total energy conservation equation in Section 3.1, and then turning to the first law of thermodynamics in Section 3.2. There we show that work done by pressure is increasingly important at high radiation-to-gas pressure ratios, and this work is associated both with radiative damping and excitation of vertical mechanical motions. We identify two important classes of radiation pressure fluctuation in Section 4 and show how radiative damping acts upon them in Section 5. We then present the detailed results on the nature of vertical advective energy transport in Section 6. In Section 7, we discuss how our findings give rise to a more dynamic view of the thermal physics of a radiation-dominated disk. Finally, we summarize our conclusions in Section 8. We provide some mathematical background on trapped vertical modes that modulate the mechanical work and advective energy transport in an Appendix.

The radiation MHD equations solved in our simulations are discussed by Hirose et al. (2006, 2009b), but we list them here again as they are the basis for most of the equations we derive elsewhere in the paper:

$$\begin{equation} {\partial \rho \over \partial t}+{\nabla }\cdot (\rho {\bf v})=0 \end{equation} \tag{ 1 }$$

$$\begin{eqnarray} {\partial \over \partial t}(\rho {\bf v})+\nabla \cdot (\rho {\bf v}{\bf v})&=& {-}\nabla p-\nabla \cdot {\sf q}+{1\over 4\pi }(\nabla \times {\bf B})\times {\bf B}\nonumber\\ && +{\bar{\kappa }^{\rm R}\rho \over c}{\bf F}+{\bf f}_{\rm SB} \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} {\partial e\over \partial t}+\nabla \cdot (e{\bf v})&=&{-}p\nabla \cdot {\bf v} -{\sf q}:\nabla {\bf v} -(aT^4-E)c\bar{\kappa }^{\rm P}_{\rm ff}\rho\nonumber\\ && -cE\kappa _{\rm es}\rho {4k_{\rm B} (T-T_{\rm rad})\over m_ec^2}+\tilde{Q} \end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray} {\partial E\over \partial t}+\nabla \cdot (E{\bf v})&=&{-}{\sf P}:\nabla {\bf v} +(aT^4-E)c\bar{\kappa }^{\rm P}_{\rm ff}\rho\nonumber\\ && +cE\kappa _{\rm es}\rho {4k_{\rm B} (T-T_{\rm rad})\over m_ec^2}-\nabla \cdot {\bf F} \end{eqnarray} \tag{ 4 }$$

$$\begin{equation} {\partial {\bf B}\over \partial t}=\nabla \times ({\bf v}\times {\bf B}) \end{equation} \tag{ 5 }$$

$$\begin{equation} {\bf F}=-{c\lambda \over \bar{\kappa }^{\rm R}\rho }\nabla E. \end{equation} \tag{ 6 }$$

Here, ρ is the density, v is the fluid velocity, B is the magnetic field, p is the pressure in the gas, ${\sf q}$ is a diagonal tensor associated with the artificial bulk viscosity adopted by the code (Stone & Norman 1992), e = 3p/2 is the internal energy density in the gas, T is the gas temperature, E is the radiation energy density, ${\sf P}$ is the radiation pressure tensor, T_rad ≡ (E/a)^1/4 is the effective temperature of the radiation, a is the radiation density constant, F is the radiation flux, k_B is Boltzmann's constant, m_e is the electron mass, κ_es is the electron scattering opacity, $\bar{\kappa }_{\rm ff}^{\rm P}=\bar{\kappa }_{\rm ff}^{\rm P}(\rho,e)$ is the Planck mean free–free opacity, $\bar{\kappa }^{\rm R}=\bar{\kappa }^{\rm R}(\rho,e)$ is the Rosseland mean opacity including electron scattering and free–free contributions,⁷ λ is a flux limiter equal to 1/3 in the optically thick limit, and $\tilde{Q}$ is the dissipation rate per unit volume required to maintain total energy conservation due to grid-scale losses of magnetic and kinetic energy. More details on how many of these quantities are actually computed in the code can be found in Hirose et al. (2006, 2009b).

The vector f_SB represents the gravitational and inertial forces in the local shearing-box frame, which is rotating at fixed angular velocity Ω with respect to an inertial reference frame:

$$\begin{equation} {\bf f}_{\rm SB}=-2\rho \Omega \hat{\bf z}\times {\bf v}-\rho \nabla \phi, \end{equation} \tag{ 7 }$$

where ϕ ≡ −3Ω²x²/2 + Ω²z²/2 is the effective gravitational potential, and x, y, z are Cartesian spatial coordinates along the radial, azimuthal, and vertical directions, respectively.

Table 1 summarizes most of the numerical parameters of the six simulations we analyze in this paper, including the surface mass density Σ and scale height H of the initial condition. Each simulation had the same angular velocity Ω = 190 rad s⁻¹, corresponding to a radius of 30GM/c² around a 6.62 M_☉ Schwarzschild black hole. The grids all had 48 zones in the x (radial) direction, 96 zones in the y (azimuthal) direction, and 896 zones in the z (vertical) direction. The simulations were all initialized with a weak magnetic field in a twisted azimuthal flux tube geometry in the midplane regions of the domain. MRI turbulence is well established by 10 orbits into the simulation, and we ran each simulation a further 250–600 orbits. This is much longer than the thermal time of the simulations, the average of which ranges from 13 orbits for 1112a to 24 orbits for 0519b, as shown in Table 2. All the simulations reach approximate thermal equilibria with continued long timescale fluctuations. Radiation pressure dominates the thermal pressure in all these simulations, and the time- and box-averaged ratios of radiation to gas pressure are also listed in Table 2. More information about these simulations can be found in Hirose et al. (2009b, 2009a).

Table 2. Summary of Time-averaged Physical Properties of Simulations

Simulation	HP/Ha	tthermal/Orbitsb	$\langle {P_{\rm rad,av}\over P_{\rm gas,av}}\rangle$c	$\langle {\rho _{\rm max}\over \rho _{\rm min}}\rangle$d
0211b	0.55	16 ± 3.2	62 ± 16	2.3 ± 0.56
0519b	0.81	24 ± 5.2	70 ± 21	2.6 ± 1.3
1112a	0.71	13 ± 2.6	6.6 ± 2.4	1.9 ± 0.54
1126b	0.88	16 ± 2.9	10 ± 2.7	2.1 ± 0.81
0520a	1.3	22 ± 5.3	14 ± 4.7	2.2 ± 0.79
0320a	1.5	21 ± 3.8	6.6 ± 1.6	1.7 ± 0.33

Notes. ^aThe thermal pressure scale height, defined as H_P ≡ ∫^∞ _{− ∞}P_therm(z)dz/[2P_therm(0)], where P_therm(z) is the time-averaged vertical profile of gas plus radiation pressure, in units of the fiducial scale height H used in the simulation grid. ^bThe time average of the thermal time, defined as the ratio of instantaneous thermal energy content in the simulation domain to the instantaneous horizontally averaged radiative flux emerging from the top and bottom faces. Errors indicate one standard deviation in the time average. ^cThe time average of the ratio of box-averaged radiation pressure E/3 to box-averaged gas pressure p. Errors indicate one standard deviation in the fluctuations of the ratio about the time average. ^dThe time average of the ratio of maximum to minimum density at the midplane z = 0. Errors indicate one standard deviation in the fluctuations of the ratio about the time average.

Download table as:  ASCIITypeset image

The shearing-box equations of motion as implemented in our code conserve total energy to high accuracy (Hirose et al. 2006). Energy originates in the work done by turbulent stresses on the shearing radial walls of the box and ultimately escapes from the top and bottom of the box, largely in the form of photons. Long-term equilibrium requires that nearly all the work done by the walls be dissipated and that all the heat generated be vented. In this section, we analyze this energy flow and dissipation in detail as a function of height.

3.1. Total Energy Conservation

Equations (1)–(5) and (7) together imply a total energy conservation equation

$$\begin{equation} {\partial \over \partial t} ({\newmathcal {E}}_{\rm mech} + {\newmathcal {E}}_{\rm therm}) + \nabla \cdot ({\bf {F}}_{\rm mech}+ {\bf {F}}_{\rm therm})={\bf v} \cdot \left(\nabla \cdot {\sf {P}}+ {\bar{\kappa }^{\rm R}\rho \over c}{\bf {F}}\right). \end{equation} \tag{ 8 }$$

The grid-scale numerical heating $\tilde{Q}$ has been lost, as our total energy numerical scheme is designed to capture grid-scale losses of magnetic and kinetic energy and convert them into heat, which is promptly used to create photons. There are also numerical energy sources and sinks that should be present on the right-hand side of this equation due to density and internal energy floors and a velocity cap (see the Appendix of Hirose et al. 2006), but these are negligible and we ignore them here.

The energy density in Equation (8) has two parts. The first is mechanical, being the sum of kinetic, effective gravitational potential, and magnetic energy densities:

$$\begin{equation} {\newmathcal E}_{\rm mech}\equiv {1\over 2}\rho v^2+\rho \phi +{B^2\over 8\pi }. \end{equation} \tag{ 9 }$$

The second is thermal, being the sum of the gas and radiation internal energy densities:

$$\begin{equation} {\newmathcal E}_{\rm therm}\equiv e+E. \end{equation} \tag{ 10 }$$

The energy flux vector also has two similar pieces. The mechanical energy flux is the sum of kinetic energy flux, effective gravitational potential energy flux, flux of work done by artificial viscosity and by gas and radiation pressures, and Poynting flux:

$$\begin{equation} {\bf F}_{\rm mech}\equiv {1\over 2}\rho v^2{\bf v}+\rho \phi {\bf v}+ {\sf q}\,{\bm \cdot}\, {\bf v}+p{\bf v}+{\sf P}\cdot {\bf v}+ {c\over 4\pi }{\bf E}\times {\bf B}, \end{equation} \tag{ 11 }$$

where E = −v × B/c is the electric field in ideal MHD. The thermal energy flux is the sum of gas and radiation internal energy advection and heat transport by radiative diffusion,

$$\begin{equation} {\bf F}_{\rm therm}\equiv (e+E){\bf v}+{\bf F}. \end{equation} \tag{ 12 }$$

The right-hand side of Equation (8) is an artificial set of energy source and sink terms that result from our use of flux-limited diffusion to handle radiation transport. These terms exactly cancel in the optically thick limit. They would also cancel to lowest order in v/c if we were using the full radiation momentum equation (Equation (9) of Stone et al. 1992) instead of the diffusion Equation (6). These terms contribute negligibly to the overall energy balance in our simulations, and we ignore them from now on.

Assuming an equilibrium has been established over long timescales, the energy conservation Equation (8) can be rewritten to show that the local divergence of the time and horizontal average of the vertical energy flux is given by the local rate at which work is being done on the fluid by the shearing walls at that height z, i.e.,

$$\begin{equation} {d\over dz}(\langle F_{{\rm mech},z}\rangle + \langle F_{{\rm therm},z}\rangle)={3\over 2}\Omega \langle \tau _{xy}(z)\rangle, \end{equation} \tag{ 13 }$$

where angle brackets denote horizontal and time averages at a particular height z. For example,

$$\begin{equation} \langle F_{{\rm mech},z}(z)\rangle ={1\over \Delta tL_xL_y}\int _0^{\Delta t} dt \int _{-L_x/2}^{L_x/2}dx\int _{-L_y/2}^{L_y/2}dy F_{{\rm mech},z}, \end{equation} \tag{ 14 }$$

where Δt is the simulation duration minus the first 10 orbits when the MRI was still in its growth phase. All time averages presented in this paper use this time interval.

The quantity 〈τ_xy(z)〉 is the time-averaged xy component of the Reynolds and Maxwell stresses at height z, plus a negligibly small contribution from radiation viscosity,

$$\begin{equation} \langle \tau _{xy}(z)\rangle \equiv {1\over \Delta tL_y}\int _0^{\Delta t}dt \int _{-L_y/2}^{L_y/2}dy\left(\rho v_x\delta v_y-{B_xB_y\over 4\pi }+P_{xy}\right). \end{equation} \tag{ 15 }$$

Here, δv_y ≡ v_y + 3Ωx/2 is the perturbation of azimuthal velocity from the average shear flow and the y-integral is done either on the inner or outer radial wall of the box (Hawley et al. 1995). When vertically integrated, the relative contributions of Maxwell and Reynolds stresses are remarkably constant from simulation to simulation: they are 85% and 15%, respectively, to within 1% for all six simulations. The radiation viscosity contributes at most 10⁻⁴ of the total vertically integrated stress.

Figure 1 shows the time-averaged vertical profiles of the dominant contributions to the vertical thermal and mechanical energy fluxes in simulations 1112a (one of the two lowest radiation to gas pressure ratio simulations we consider in this paper) and 0519b (the highest radiation to gas pressure ratio). As in the case of gas pressure-dominated simulations (Hirose et al. 2006) and simulations with comparable gas and radiation pressure (Krolik et al. 2007), radiative diffusion is the dominant process of vertical energy transport in simulation 1112a. However, we now see that at the highest levels of radiation pressure support simulated thus far (0519b), radiation advection is just as important in the midplane regions.

Zoom In Zoom Out Reset image size

Figure 2 depicts the time-averaged vertical profiles of the various contributions to the divergences of the vertical thermal and mechanical energy fluxes and compares their sum to the vertical profile of stress times rate of strain. Equation (13) is accurately satisfied by all the simulations. The radiation diffusion flux divergence approximately matches the cΩ²/κ_es value required by hydrostatic equilibrium, presumably because departures from this equilibrium would result in very fast readjustments on the dynamical timescale. All the remaining flux divergence components are all positive in the midplane regions and negative further out. The energy injected by the stresses on the walls in the midplane regions is significantly larger than the critical value which can be carried away by radiative diffusion and still maintain vertical hydrostatic equilibrium. This excess energy is therefore transported outward by the other forms of energy flux, among which radiation advection is the most important. This process is generally completely ignored in standard accretion disk models.

Zoom In Zoom Out Reset image size

3.2. The First Law of Thermodynamics: PdV Work and Radiative Damping

By combining Equations (3) and (4), energy conservation may also be described in terms of the first law of thermodynamics,

$$\begin{equation} {\partial {\newmathcal E}_{\rm therm}\over \partial t}+ \nabla \cdot {\bf F}_{\rm therm}=Q-p\nabla \cdot {\bf v}- {\sf P}:\nabla {\bf v}. \end{equation} \tag{ 16 }$$

Here, $Q\equiv \tilde{Q}-{\sf q}:\nabla {\bf v}$ is the local dissipation rate per unit volume, being a sum of grid-scale conversions of magnetic and kinetic energy into heat as well as energy dissipated by artificial viscosity. (The artificial viscosity dissipation is approximately 3/4 of the grid-scale loss rate of kinetic energy in all of our simulations and has a similar time- and horizontally averaged vertical profile.) Again assuming an equilibrium has been established over long timescales, Equation (16) implies that

$$\begin{equation} {d\over dz}\langle F_{{\rm therm},z}\rangle = \langle Q-p\nabla \cdot {\bf v}-{\sf P}:\nabla {\bf v}\rangle. \end{equation} \tag{ 17 }$$

In other words, the local divergence of the average vertical thermal energy flux is equal to the local dissipation per unit volume plus the rate at which gas and radiation pressure work is being done on the plasma per unit volume.

Equation (17) appears to be familiar, but its form hides some subtleties peculiar to radiation-dominated disks. These new effects emerge when comparing the two sides of this equation. Figure 3 depicts the vertical profiles of the advection and radiative diffusion contributions to the divergence of the vertical thermal flux, as well as the numerical dissipation profile, for all six of our simulations. Because our numerical dissipation Q does not include damping by radiative diffusion, one might have expected it to be smaller than the divergence of the thermal flux. In fact the opposite is true, with the largest discrepancy occurring in the most radiation pressure-dominated simulation 0519b. Moreover, the absence of an explicit place for radiative damping suggests that something is missing because that process should contribute to the total dissipation rate.

Zoom In Zoom Out Reset image size

Figure 4 shows the same vertical profiles of thermal flux divergence, but a dissipation rate corrected by the PdV work. With this adjustment, the time-averaged first law of thermodynamics (17) is accurately satisfied in each of the simulations. From this fact, we reach two important conclusions. First, there is a significant conceptual flaw in classical time-steady accretion disk models. These models assume that the divergence of the diffusive radiation flux completely defines the left-hand side of the first law of thermodynamics, while the local dissipation rate is the only contribution to its right-hand side. But we have just seen that in radiation-dominated disks both of these simplifications are wrong: radiation advection must be included with radiative diffusion, and the PdV work terms are important. Second, radiative damping is actually included in the PdV work terms. In our scheme, the dissipation associated with it is conveyed by the diffusion equation, but the energy it dissipates flows into the gas via the pressure work terms in the first law of thermodynamics.

Zoom In Zoom Out Reset image size

It is perhaps helpful to make more explicit the actual dissipation associated with radiative diffusion. The first law of thermodynamics (16) can be combined with the radiative diffusion Equation (6) to derive an equation for the evolution of the entropy per unit mass s of the gas plus radiation mixture. Restricting consideration to optically thick regions for simplicity, this equation is

$$\begin{equation} \rho \left({\partial s\over \partial t}+{\bf v}\cdot \nabla s\right)+ \nabla \cdot \left({{\bf F}\over T}\right)={4acT\over 3\bar{\kappa }^{\rm R}\rho } \left(\nabla T\right)^2+{Q\over T}. \end{equation} \tag{ 18 }$$

The last term on the left-hand side is the divergence of the entropy flux due to radiative diffusion. The first and second terms on the right-hand side are sources of entropy (dissipation) due to radiative diffusion and grid-scale losses plus artificial viscosity, respectively. Assuming an equilibrium has been established on long timescales, this equation becomes

$$\begin{equation} \langle \rho {\bf v}\cdot \nabla s\rangle +{d\over dz}\left<{F_z\over T}\right>= \left<{4acT\over 3\bar{\kappa }^{\rm R}\rho }\left(\nabla T\right)^2\right>+ \left<{Q\over T}\right>. \end{equation} \tag{ 19 }$$

Advective and diffusive transport of entropy is therefore balanced by dissipation. The dissipation associated with radiative damping of fluctuations in the turbulence is in the first term on the right-hand side, but this term also includes entropy generation due merely to the photon diffusion down the background average vertical temperature gradient. An analogous situation holds in stars with radiative envelopes. Under static conditions, the outward luminosity at every radius in the envelope is constant. Because the temperature at the base of the envelope is higher than the temperature of the photosphere, the entropy leaving the photosphere exceeds the entropy entering the base of the envelope. The source of this entropy increase is the dissipative nature of photon diffusion itself and is not associated with dissipative release of energy that must then be transported away to establish thermal equilibrium.

Because of this non-energy releasing dissipation, we cannot use Equation (19) to calculate the true dissipative heating due to radiative damping of fluctuations. Instead, we must try and use the pressure work terms in the first law of thermodynamics (17). As we noted above, it is through these terms that energy flows into the gas from turbulent fluctuations and ultimately dissipates by radiative diffusion.

However, somewhat surprisingly, the time-averaged pressure work corrections $-\langle p\nabla \cdot {\bf v}\rangle -\langle {\sf P}:\nabla {\bf v}\rangle$ are typically negative, in contrast to the case of shearing boxes without vertical gravity studied by Turner et al. (2002, 2003). This sign change indicates that the plasma does net work through expansion even though these corrections must include radiative damping, which would be positive (it is dissipative). The fact that the overall pressure work corrections are negative means that the gas loses more thermal energy by driving expansion away from the midplane than it gains by radiative damping.

Fortunately, these two contributions to the pressure work corrections (pumping expansion and radiative damping) can be separated analytically. Considering only the optically thick regions, write the thermal pressure P_therm = p + E/3 and fluid velocity as P_{therm, av} + δP_therm and v_av + δv, respectively. Here, the "av" subscript denotes horizontal average and δP_therm and δv therefore have zero horizontal average by definition. The divergence of δv also has zero horizontal average, as one can show by integrating by parts and using the shearing-box boundary conditions. Hence the horizontal average of the thermal pressure work done on the plasma under optically thick conditions is

$$\begin{eqnarray} &&{1\over L_xL_y}\int _{-L_x/2}^{L_x/2} dx\int _{-L_y/2}^{L_y/2}dy \left(-p\nabla \cdot {\bf v}-{\sf P}:\nabla {\bf v}\right)= -P_{\rm therm,av}\nonumber\\ &&\quad \times{\partial \over \partial z}v_{{\rm av},z}(z,t) +{1\over L_xL_y}\!\!\int _{-L_x/2}^{L_x/2} dx\!\int _{-L_y/2}^{L_y/2}dy (-\delta P_{\rm therm}\nabla \cdot \delta {\bf v}).\nonumber\\ \end{eqnarray} \tag{ 20 }$$

After time averaging, the first term on the right-hand side represents pumping of vertical mechanical motions. This term will be negative if the plasma undergoes vertical expansion in a horizontally averaged sense, i.e., the plasma will do net work. We will henceforth call this term the mechanical pumping term.

The second term on the right-hand side of Equation (20) is the radiative damping. For a sinusoidal adiabatic fluctuation, this term would vanish identically in the time average, because δP_therm and ∇ · δv are 90° out of phase. Radiative diffusion removes this cancellation by making δP_therm lead ∇ · δv by a little more than 90°. For example, if the plasma is locally compressed nearly adiabatically, the temperature rises and radiative diffusion to the cooler surrounding regions increases. This diffusion is fastest as the point of maximum compression is approached, while at the same time the rate of compression is slowing down. The temperature therefore starts to drop just before the point of maximum compression, i.e., the maximum in temperature is reached before the point of maximum compression. The same holds on the expansion cycle, where the minimum in temperature is reached before the point of maximum expansion. Hence there is a greater than 90° lead in δP_therm relative to ∇ · δv.⁸ This phase offset results in a net positive time average of −δP_therm∇ · δv. The plasma has net positive work done on it, and that excess work is dissipated by radiative diffusion. In principle, the mechanical pumping term of Equation (20) may still include a small positive contribution from the radiative damping of purely vertical acoustic waves, but we show in Section 6.1 below that this contribution is negligible.

Figures 5 and 6 depict the time-averaged vertical profiles of stress times rate of strain and the different contributions to dissipation and mechanical work for simulations 1112a and 0519b, respectively. We did not save high time resolution data over the entire simulation duration to enable us to directly compute the last term on the right-hand side of Equation (20), so the radiative damping profiles in these figures were computed from subtracting the vertical profile of the first term from the vertical profile of the left-hand side. We also neglected gas pressure work here as it is very small in these simulations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Two features are worth noting about these profiles. First, even after time averaging, there remain spatial fluctuations in the total mechanical work profile as well as the profile of mechanical pumping. Those fluctuations are completely absent in the radiative damping profile, which is as smooth as the other (magnetic and kinetic) dissipation profiles. The radiative damping profile is also very similar in shape to these other dissipation profiles. Second, when compared to the numerical dissipation, the radiative damping and mechanical pumping, as well as the net total pressure work, are clearly relatively more important in 0519b (the simulation with the highest radiation to gas pressure ratio) than in 1112a. As much as 22% of the work done by the shearing walls ends up being dissipated by radiative diffusion in simulation 0519b.⁹

Figure 7 illustrates this trend of increasing relative importance of the pressure work corrections with growing radiation to gas pressure ratio. A larger rate of radiative damping presumably requires larger density fluctuations; to test this supposition, we have measured the time average of the ratio of maximum to minimum density at the midplane in each of the simulations, and the results are listed in Table 2. Density fluctuations in the midplane regions do indeed become slightly larger with increasing radiation to gas pressure. Figure 8 shows the fractional pressure work as a function of the time-averaged density contrast at the midplane. The large horizontal error bars reflect the large variations in the density contrast, but the radiative damping points of this figure are clearly consistent with the trend observed in non-stratified shearing boxes by Turner et al. (2003, see the top panel of their Figure 7).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

There is considerable structure in the temporal behavior of the pressure work terms. A 10 orbit segment of the time dependence of the vertical integral of the stress times rate of strain, the dissipation terms, and the pressure work terms for simulation 0519b is shown in Figure 9. The mechanical pumping work (lower gray dashed curve) shows clear oscillatory behavior on timescales of order the orbital period. This is completely absent in the radiative damping (upper gray dashed curve), which instead clearly exhibits much higher frequency variability. Both also exhibit much longer timescale variation.

Zoom In Zoom Out Reset image size

Some aspects of this behavior can be immediately understood by Fourier transforming the time dependence and plotting the temporal power spectrum of the vertically integrated mechanical pumping and radiative damping terms. The result is shown in Figure 10. The large oscillations seen in the mechanical pumping in Figure 9 are reflected in a series of discrete sharp peaks. These peaks represent standing vertical acoustic waves that are trapped in the box. Despite their prominence, we show below in Section 6.1 that they actually contribute negligibly to the energetics of the disk. The radiative damping power spectrum exhibits much higher frequency power, including a set of high frequency, extremely sharp spikes. This is consistent with the high-frequency variability seen in the time domain of Figure 9. To better understand this variability, we must first examine the types of fluctuation that can produce radiative damping, to which we now turn.

Zoom In Zoom Out Reset image size

It is radiation pressure (i.e., temperature) fluctuations that give rise to photon diffusion and therefore radiative damping. In this section, we discuss two distinct types of such fluctuations that clearly play a role in our simulations. We will also see in Section 6 that understanding both radiation advection and energy transport by Poynting flux will likewise be aided by a prior understanding of radiation pressure fluctuations.

In classical MHD theory (in which radiation pressure is negligible), linear compressible waves are classified into "fast" and "slow" modes. One way to understand qualitatively this distinction is to note that the magnetic and gas pressure perturbations are exactly in phase in the former mode and exactly out of phase in the latter; it is the partial cancellation of the pressure perturbation in the slow mode that causes its low propagation speed. Our simulations do not exactly correspond to this categorization in two ways: the fluctuations are often nonlinear; and radiation pressure is both significant and only imperfectly coupled to the fluctuations (that is what radiative damping is all about, of course). Nonetheless, this conceptual division remains useful because large total (gas plus radiation plus magnetic) pressure fluctuations at a given length scale tend to have much higher frequency than those with small total pressure fluctuations. The large total pressure fluctuations therefore propagate on an approximately stationary background set by the slowly evolving small total pressure fluctuations. We will call the former "acoustic waves" and the latter "isobaric waves" in order to stress this distinction.

4.1. Acoustic Waves

Agol & Krolik (1998) suggested that acoustic waves (i.e., fast magnetosonic waves) would be the dominant contributor to radiative damping in radiation-dominated accretion disks. We have already seen that standing vertical acoustic waves are excited in the box. In addition, strong nonaxisymmetric acoustic waves are pervasive in all our simulations. These waves are almost certainly stochastically excited by the MRI turbulence itself (Heinemann & Papaloizou 2009a, 2009b). In these waves, gas density, radiation pressure, and magnetic pressure fluctuations all oscillate in phase. Figure 11 depicts a snapshot of various quantities at the z = 0 midplane of simulation 0519b at 200.3 orbits. A wave pattern is clearly evident in most of the fluid quantities shown, though it is cleanest in the total pressure (upper left). Figure 12 depicts a vertical slice of the midplane regions at the same time. Particularly in the total pressure (upper left), the wave is once again evident in the vertical yellow and blue stripes in the midplane region |z| ≲ 0.2–0.3H and primarily propagates in the horizontal direction.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To better understand these nonaxisymmetric waves, it is helpful to project the spatial variation of the z = 0 midplane fluid quantities onto the natural set of basis vectors exp [i(k_x(t)x + k_yy)] of the shearing box, where

$$\begin{equation} k_x(t)={2\pi n_x\over L_x}+{3\over 2}\Omega k_y t, \end{equation} \tag{ 21 }$$

$$\begin{equation} k_y={2\pi n_y\over L_y}, \end{equation} \tag{ 22 }$$

and n_x and n_y are the integer quantum numbers (Hawley et al. 1995; Heinemann & Papaloizou 2009b). The largest azimuthal wavelength (|n_y| = 1) waves always have the most power. Figure 13 shows the real part of the total pressure Fourier amplitude as a function of time for n_y = 1 and various values of n_x. In agreement with Heinemann & Papaloizou (2009b), waves with different values of the n_x quantum number reach high amplitude at distinct times (in declining n_x order) when they swing from leading to trailing, resulting in a series of wave pulses. The separation in time between these pulses is determined entirely by the time interval between successive epochs at which the shearing radial boundaries become exactly periodic in the radial direction. This shear time of the box is δT_s = 2L_y/(3ΩL_x), i.e., 3π/4 ≃ 2.36 inverse orbital periods. We purposely chose the 200.3 orbits epoch for Figures 11 and 12 as it corresponds to the time of peak amplitude for one of these wave pulses (n_x = 9, n_y = 1).

Zoom In Zoom Out Reset image size

4.2. Isobaric Waves

The nonaxisymmetric acoustic wave pattern is evident in the spatial distribution of density, shown in the top right panel of Figure 11, but this pattern is markedly perturbed by shorter length scale fluctuations. Among the most prominent are rarefied regions, e.g., near (x = 0.05, y = −0.5), that are correlated with regions of low radiation pressure, low gas pressure, and high magnetic pressure (bottom left to right panels of Figure 11, respectively).

These are a second kind of radiation pressure fluctuation, one in which the magnetic pressure oscillates with very nearly the same amplitude as the sum of gas and radiation pressure, but with opposite phase. As a result, the total pressure hardly changes. This near-cancellation of pressure fluctuations is characteristic of slow magnetosonic modes, in which the cancellation is exact to ∼O(β⁻¹) when the plasma β ≫ 1. When such modes are placed in a rotating shear flow whose rotation rate declines outward, they become the MRI instability. We might naturally expect to see many such features in our shearing-box simulations, and some can be expected to grow to nonlinear amplitude.

Figure 14 shows correlation plots between magnetic pressure fluctuations, density fluctuations, and thermal (i.e., gas plus radiation) pressure fluctuations in the z = 0 midplane of simulation 0519b near the 200 orbit epoch. Most cells (i.e., most of the area as indicated by the red and yellow regions in all three panels) house weakly negative magnetic fluctuations, and the thermal pressure and density fluctuations in these cells approximately obey the expected adiabatic relation for acoustic waves: δP_therm/P_therm = Γ₁δρ/ρ, where Γ₁ ≃ 4/3 is the first generalized adiabatic index for a gas and radiation mixture (Chandrasekhar 1967). However, in the upper right plot of magnetic pressure perturbation versus thermal pressure perturbation, these acoustic waves oscillate horizontally back and forth about a mean thermal pressure perturbation that is set by a background of isobaric fluctuations. This is indicated by the fact that the contours of cell-counts stretch diagonally across this plot, showing that at large amplitude the magnetic pressure and the total thermal (mostly radiation) pressure are anti-correlated. In addition, the magnitudes of the two fluctuations are similar, so that |δP_mag + δP_therm| is in general considerably smaller than |δP_mag| + |δP_therm|. The isobaric modes also exhibit a clear anti-correlation between magnetic pressure and density, presumably because the density is positively correlated with the gas pressure component of the thermal pressure.

Zoom In Zoom Out Reset image size

While most of the volume is occupied by weak magnetic field, Figure 15 shows that the regions of enhanced magnetic pressure associated with the isobaric modes are numerous enough to dominate the total magnetic energy budget in the vicinity of the midplane. As much as half the magnetic energy can be located in regions where the magnetic pressure is more than twice the local horizontal average. The maximum enhancement ratio in a single horizontal slice can be as high as ∼10 at certain times and places.

Zoom In Zoom Out Reset image size

Any phenomenon involving radiation pressure fluctuations (acoustic or isobaric) can be dissipated (i.e., lose its mechanical energy to heat, with associated entropy production) via photon diffusion, which enters our formalism in the diffusion equation. In fact, we find that both types of radiation pressure fluctuation contribute significantly to the radiative damping, although the isobaric perturbations dominate.

One way of seeing this is to examine the integrated radiative damping as a function of fractional magnetic and thermal pressure fluctuations. This is shown in the upper left plot of Figure 16 for the 200–202 orbit interval at the z = 0 midplane in simulation 0519b. (The corresponding incidence of fractional magnetic and thermal pressure fluctuations themselves was shown in the upper right-hand plot of Figure 14.) Two distinct regions in this plane contribute to the radiative damping, one associated with large positive magnetic fluctuations stretching up the isobaric locus, and one associated with slightly negative magnetic fluctuations that can occur with both positive and negative thermal pressure fluctuations that are mostly due to acoustic waves.

Zoom In Zoom Out Reset image size

Integrating the data of the two-dimensional distribution function over all thermal pressure perturbation values results in the distribution shown in the upper right plot of Figure 16. Because weak magnetic fields (negative magnetic pressure perturbations) dominate the spatial volume in the midplane, radiative damping from negative magnetic pressure perturbations here is primarily due to acoustic waves. In contrast, positive magnetic pressure perturbations are mostly associated with the strongly magnetized, nonlinear isobaric fluctuations. The contributions to radiative damping from these two types of fluctuation are comparable: 35% for acoustic waves, and 65% for isobaric fluctuations for the 200–202 orbit epoch in simulation 0519b, for a total midplane radiative damping rate of 9 × 10¹⁴ erg cm⁻³ s⁻¹.

Because the acoustic waves are much faster oscillating fluctuations, they are presumably responsible for the rapid variability seen in the time dependence of the box-integrated radiative damping (upper dashed gray curve) of Figure 9. In fact, the origin of the high-frequency spikes in the power spectrum of the box-integrated radiative damping shown in Figure 10 is now clear. All of these spikes are higher order harmonics of a fundamental frequency of ≃ 2.36 inverse orbital periods. This frequency is exactly the inverse of the shear time of the box and reflects the time interval between successive pulses of the dominant |n_y| = 1 nonaxisymmetric waves (see Figure 13).

The measured contributions to the radiative damping can also be estimated analytically from the observed fluctuation amplitudes, at least for the nonaxisymmetric fast waves, which are not too nonlinear. All the fluctuations that we observe in the midplane regions have length scales such that the photon diffusion time is much longer than the sound crossing time. When radiation pressure dominates both magnetic and gas pressure, the exponential damping rate of the amplitude of a linear, plane wave fast mode (acoustic wave) in this slow diffusion limit is given by

$$\begin{equation} \Gamma _{\rm acoustic}\simeq {k^2c\over 6\bar{\kappa }^{\rm R}\rho }, \end{equation} \tag{ 23 }$$

where k is the wavenumber of the mode. Physically, acoustic waves in a radiation-dominated plasma are simply damped at the rate that photons diffuse across a wavelength. The acoustic radiative damping rate can therefore be estimated by multiplying the average wave energy density by twice this damping rate (twice because energy is proportional to the amplitude squared), i.e.,

$$\begin{eqnarray} Q_{\rm rad,acoustic}&=&\left<{(\delta P_{\rm tot})^2\over 2c_{\rm t}^2\rho }+ {1\over 2}\rho (\delta v)^2\right>(2\Gamma)= {(\delta P^2_{\rm tot,max})^2\over \rho c_{\rm t}^2}\Gamma\nonumber\\ &=&\left({k^2c \over 6\bar{\kappa }^{\rm R}}\right)\left({3c_{\rm t}\over 4} \right)^2 \left({\delta P_{\rm tot,max}\over P_{\rm tot}}\right)^2, \end{eqnarray} \tag{ 24 }$$

where c_t ≃ 6 × 10⁸ cm s⁻¹ is the total sound speed in the radiation-dominated plasma (see Figure 23). From Figure 11, we estimate a typical pressure amplitude of δP_{tot, max}/P_tot ∼ 0.06 and a wavenumber of k ∼ 2π/(0.3H). For these numbers, we find Q_{rad, acoustic} ∼ 3 × 10¹⁴ erg cm⁻³ s⁻¹, very close to our measured acoustic contribution near this epoch of 35% of 9 × 10¹⁴ erg cm⁻³ s⁻¹.

The short two-orbit epoch we have analyzed in detail in this section is not atypical for this simulation: the midplane radiative damping rate of 9 × 10¹⁴ erg cm⁻³ s⁻¹ at this time is close to the midplane radiative damping rate time-averaged over the entire ∼400 orbits of the simulation (upper gray dashed curve of Figure 6). We reiterate here how significant these numbers are (as much as 22% of the work done by the walls when integrated over the entire box) for the overall energy budget at these high levels of radiation pressure support.

Although the contributions of acoustic waves and isobaric fluctuations to radiative damping are comparable near the midplane, plots similar to Figure 16 at fixed heights more than ≃ 0.5H away from the midplane clearly show that isobaric fluctuations become dominant at these higher altitudes. At these locations, the two-dimensional distribution of radiative damping as a function of magnetic and thermal pressure closely follows the isobaric locus for both positive and negative magnetic perturbations, in contrast to the midplane distribution (upper left plot of Figure 16). This fact is consistent with the patterns of fluctuations that we observe as a function of height in Figure 12. In that figure, we see that the vertical yellow and blue stripes near the midplane associated with nonaxisymmetric acoustic waves in total pressure (upper left plot) extend out only to |z| ∼ 0.2–0.3H. Beyond that, isobaric fluctuations of both signs of thermal and magnetic pressure dominate the volume (lower left and right plots).

We now turn to examine two possible origins of the vertical radiation advection energy flux: transport by trapped vertical waves and buoyancy.

6.1. Vertical Epicyclic and Acoustic Waves

Figure 17 shows the temporal power spectrum of the vertical radiation advection energy flux at every height z in simulation 0519b. A hierarchy of vertical modes is apparent at frequencies at and above the orbital frequency. All but the lowest frequency of these modes are identical to those causing the oscillations and discrete frequency peaks in the mechanical pumping shown in Figures 9 and 10.

Zoom In Zoom Out Reset image size

The lowest order mode at the orbital frequency is just a vertical epicyclic oscillation, excited by vertical asymmetries in momentum losses through the vertical boundaries of the box. The vertical velocity v_z in this mode is the same at all heights, and the mode simply uniformly displaces the entire plasma up and down in the box. It produces zero compression or rarefaction and therefore does no pressure work, which explains its complete absence in the temporal power spectra of pressure work terms shown in Figure 10. Because it bodily displaces trapped photons up and down, it modulates the radiation advection energy flux, but it produces zero radiation advection through the plasma.

The next lowest frequency mode is a vertical breathing mode, and higher and higher frequency modes clearly represent standing acoustic waves with increasing numbers of vertical velocity nodes. The vertical velocities of these modes also modulate the radiation advection, but the net time-averaged flux would be zero if these modes were purely adiabatic. However, radiative damping of these modes does produce a net radiation advection energy flux through the plasma over time. With a small amount of analytic work (see the Appendix for details), we can directly calculate this flux and compare it to the measured value.

The strongest mode is the breathing mode, and we measured its velocity amplitude as follows. First, we computed unwindowed, unbinned power spectra of the vertical velocity time series at each height in the box. We then summed the power at each height over the measured mode line profile in these power spectra (corresponding to the lowest peak in the upper green and black curves of Figure 10) and fit the resulting power at each height in the midplane regions with the square of the expected mode eigenfunction (linear for the breathing mode). Finally, we used Parseval's theorem to determine the actual velocity amplitude of the mode. The results were ≃ (1–7) × 10⁶ cm s⁻¹ (roughly 1% of the thermal sound speed) across the simulations at one pressure scale height from the midplane.

The damping rate of the breathing mode is (see Equation (A16))

$$\begin{equation} \Gamma _{\rm rad}\sim {c\Omega ^2\over \kappa _{\rm es} E_0}, \end{equation} \tag{ 25 }$$

where E₀ is the radiation energy density at the midplane. Not surprisingly, this is approximately the reciprocal of the nominal thermal time of the disk.¹⁰ Because the breathing mode is the longest vertical wavelength acoustic mode, it should be damped at a rate given roughly by the rate at which photon diffusion can cool the disk.

The radiative damping of the breathing mode causes phase changes between the velocity and pressure perturbations, which result in a secular flux of mechanical energy carried by the mode. This energy flux is given by (see Equations (A19) and (A20) in the Appendix)

$$\begin{equation} \left<\delta v\left[\Delta P_{\rm therm}+ \Delta \left({B^2\over 8\pi }\right)\right] \right>={3z\Gamma _{\rm rad}A_{\rm b}^2\over 14\Omega ^2H_P^2} \left({4E\over 9}+{B^2\over 4\pi }\right){,} \end{equation} \tag{ 26 }$$

where A_b is the velocity amplitude of the breathing mode at one pressure scale height H_P from the midplane.

Evaluating Equations (25) and (26), we find a peak radiation advection flux of only 7 × 10¹⁶ erg cm⁻² s⁻¹ at the breathing mode amplitude observed in the 0519b simulation; the other simulations have comparable fluxes, although the ones with lower radiation to gas pressure ratios are smaller by factors of several. However, Figure 1 shows that the actual time-averaged peak radiation advection flux in 0519b is more than five orders of magnitude larger: (1–2) × 10²² erg cm⁻² s⁻¹. The breathing mode amplitude would have to be at least ∼400 times larger, i.e., several times larger than the sound speed, in order to produce that large a flux. Similar discrepancies exist in all six simulations, and we therefore conclude that the acoustic waves do not contribute significantly to the time-averaged outward radiation advection. Instead, they simply modulate it at high frequencies.

Similar calculations can be used to estimate the dissipation rate due to radiative damping of the breathing mode. For example, in simulation 0519b, this is ≃ 4 × 10¹⁰ erg cm⁻³ s⁻¹ at the midplane, orders of magnitude smaller than the actual rates of radiative damping and pressure work that we find in the midplane regions of this simulation (Figure 6). The same is true for all the other simulations, and we conclude that the standing vertical acoustic waves contribute negligibly to radiative damping and the time-averaged mechanical pumping work, although they do significantly modulate the latter.

6.2. Buoyancy as the Driver of Advection

We have just seen that, while the radiation advection is modulated by vertical acoustic modes, these waves do not contribute significantly to the net radiation advection flux. If we smooth over this rapid modulation, a clear pattern emerges as shown in Figure 18. The radiation advection flux is modulated on long (∼5 orbit) timescales, and this modulation is in fact highly correlated with a similar modulation in the smoothed spacetime pattern of vertical Poynting flux shown in Figure 19. These spacetime patterns of Poynting flux have been seen in all previous vertically stratified MRI simulations, both with and without explicit thermodynamics (Brandenburg et al. 1995; Stone et al. 1996; Miller & Stone 2000; Hirose et al. 2006; Davis et al. 2010; Shi et al. 2010), but their origin has remained a mystery.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We argue here that both the radiation advection and the vertical Poynting flux are due to buoyancy in the region inside |z| ∼ H. That there is any buoyancy at all in the midplane regions is at first surprising, as the time- and horizontally averaged pressure and density profiles are linearly stable to buoyant instabilities in the midplane regions. Indeed, when the square of the hydrodynamic Brunt–Väisälä frequency,

$$\begin{equation} N^2=g\left({1\over \Gamma _1}{d\ln P_{\rm therm}\over dz}- {d\ln \rho \over dz}\right){,} \end{equation} \tag{ 27 }$$

is computed using horizontally and time-averaged pressure and density profiles, it is positive everywhere off the midplane in all our simulations, indicating a hydrodynamically stable average vertical entropy profile.

Magnetic fields can still cause buoyancy instabilities in the form of the interchange and undulatory Parker modes. The interchange mode is stable provided the ratio of magnetic field strength to mass density does not decrease outward too fast (e.g., Acheson 1979). In our simulations, the time-averaged profiles have d/dzln (B/ρ) > 0 everywhere off the midplane, so the interchange mode is linearly stable.

This leaves the undulatory Parker modes, which are typically linearly unstable in the surface layers (Blaes et al. 2007). In the midplane regions, where radiative diffusion is slow, the linear stability criterion (Newcomb 1961) can be expressed in terms of the requirement that the square of the magnetic Brunt–Väisälä frequency be positive. This can be defined as

$$\begin{equation} N_{\rm mag}^2\equiv g\left(-{g\over c_{\rm t}^2}-{d\ln \rho \over dz}\right) =N^2+{gv_{\rm A}^2\over c_{\rm t}^2}{d\ln B\over dz}, \end{equation} \tag{ 28 }$$

where c_t ≡ (Γ₁P_therm/ρ)^1/2 is the adiabatic sound speed in the gas plus radiation mixture, v_A is the Alfvén speed, and the last equality follows if the field is purely horizontal with no vertical tension forces. The solid lines in Figure 20 show that the midplane regions are stable to the undulatory Parker modes by this criterion. If radiative diffusion is fast enough to suppress temperature fluctuations, the undulatory Parker stability criterion (Gilman 1970) can instead be written in terms of a modified magnetic Brunt–Väisälä frequency:

$$\begin{equation} N_{\rm mag,r}^2\equiv {gv_{\rm A}^2\over \big(c_{\rm i}^2+v_{\rm A}^2 \big)} {d\ln B\over dz}, \end{equation} \tag{ 29 }$$

where c_i ≡ (p/ρ)^1/2 is the isothermal sound speed in the gas. The dashed lines in Figure 20 show the profiles of N²_{mag, r}. Hence even if radiative diffusion were fast in the midplane regions, they would still be stable simply because the time- and horizontally averaged magnetic pressure typically increases outward in these regions.

Zoom In Zoom Out Reset image size

The midplane regions, viewed in terms of horizontally averaged quantities, are therefore linearly stable to buoyant perturbations. Nevertheless, finite amplitude perturbations with three-dimensional structure within the turbulence can still be locally buoyant. As we have already seen, magnetosonic slow modes create regions of high magnetic pressure and low gas density. Here they have nonlinear amplitude, but retain these qualitative characteristics. Because they have the same pressure as other locations at that altitude, but have lower density, they are, of course, buoyant. Such low density, high magnetic pressure fluctuations are clearly present throughout the region within ±0.7H of the midplane shown in Figure 12.

Unlike the breathing modes, these fluctuations do carry enough energy flux to explain the advection. Figure 21 shows the two-dimensional distribution of radiation advection flux at z = +0.5H in simulation 0519b, averaged over 200–208 orbits, as a function of magnetic pressure and density perturbations. Integrating over all densities and magnetic fields to compute the total net radiation advection gives 9 × 10²¹ erg cm⁻² s⁻¹, which is of course consistent with the horizontally and time-averaged radiation advection at this height around this epoch (Figure 18). It is also consistent with the average radiation advection flux that we see at this height over the entire duration of the simulation (Figure 1).

Zoom In Zoom Out Reset image size

Most importantly, it is clear that magnetic pressure and density continue to be anti-correlated above the midplane, just as they are in the midplane. Moreover, outward radiation advection tends to be associated with low density and high magnetic pressure, and inward radiation advection tends to be associated with high density and low magnetic pressure. This can also be seen if we integrate over all densities (the distribution in the upper right plot) or all magnetic pressures (the distribution in the lower plot). Thus, the advection is clearly associated with these local regions of buoyancy.

A corollary of this finding is that the regions responsible for the bulk of upward radiation advection actually have radiation energy density somewhat smaller than other regions at the same altitude. In other words, this process, although it bears some resemblance to classical convection, differs strikingly in that it carries heat upward in fluid elements that are initially cooler than their surroundings. Instead of being due to a higher entropy per unit mass (which in fact they have in spite of being cooler than their surroundings), their buoyancy is driven by a higher magnetic pressure per unit mass.

The fact that cooler than average buoyant fluid elements can still give rise to a net outward energy flux can also be seen by the following simple argument.¹¹ In any given horizontal plane, divide the radiation energy density, mass density, and vertical velocity into their horizontal average plus perturbations (just as we did with other variables to separate radiative damping from mechanical pumping in the pressure work terms in Section 3.2 above): E = E_av + δE, ρ = ρ_av + δρ, and v_z = v_{z, av} + δv_z. We observe zero net vertical mass flux in the simulations, so v_{z, av} = −(δρδv_z)_av/ρ_av. The horizontal average of the radiation advection energy flux is therefore given by

$$\begin{equation} {(Ev_z)_{\rm av}\over E_{\rm av}}=\left[\delta v_z \left({\delta E\over E_{\rm av}}-{\delta \rho \over \rho _{\rm av}}\right) \right]_{\rm av}. \end{equation} \tag{ 30 }$$

The fractional density fluctuations are generally much larger in magnitude than the fractional temperature (or radiation pressure or energy density) fluctuations, as we can see in Figures 11 and 12. It is the sign of these density fluctuations that give rise to an overall outward mean radiation advection energy flux. Put another way, the underdense regions give rise to a net outward horizontal average velocity, and because the temperature fluctuations are so small, this net outward velocity times the average radiation energy density is approximately the net outward radiation advection flux.

While the rising magnetized fluid parcels start out cooler than their surroundings, they nevertheless eventually become hotter than their surroundings, and it is for this reason that they are able to transfer heat from the hot midplane regions to the cooler regions away from the midplane. We can see this qualitatively through the following linearized treatment. Consider a bundle of (mostly horizontal) field lines that rises an infinitesimal height Δz from some initial height z. We take the gas pressure to be negligible and assume that radiative diffusion is slow so that the plasma inside the bundle undergoes an approximately adiabatic change (as shown by Figure 23, in this region the diffusion speed on a length scale 0.1H is about comparable to the typical upward advective speed in this region). We further assume that the mean gradient in the total pressure (radiation plus magnetic) is in hydrostatic balance with gravity. Pressure equilibrium with the surroundings then implies that the radiation pressure P^⋆_rad of the plasma inside the bundle will change by

$$\begin{equation} \Delta P_{\rm rad}^\star ={-\rho g\Delta z\over 1+3P_{\rm mag}^\star / 2P_{\rm rad}^\star }, \end{equation} \tag{ 31 }$$

where P^⋆_mag is the initial magnetic pressure inside the bundle. On the other hand, the radiation pressure in the local background will change by

$$\begin{equation} \Delta P_{\rm rad}={-\rho g\Delta z\over 1+{\it dP}_{\rm mag}/{\it dP}_{\rm rad}}, \end{equation} \tag{ 32 }$$

where dP_mag/dP_rad is a derivative following background pressures as they vary with height. We therefore find that the rising fluid parcel will cool less rapidly with height than the background provided

$$\begin{equation} \left({P_{\rm mag}\over P_{\rm mag}^\star }\right)\left({P_{\rm rad}^\star \over P_{\rm rad}}\right){d\ln P_{\rm mag}\over d\ln P_{\rm rad}}<{3\over 2}. \end{equation} \tag{ 33 }$$

This is nearly always true in our simulations, as the bundle starts out more magnetized and cooler than its surroundings, and the background magnetic pressure generally falls less rapidly with height than the radiation pressure (see, e.g., Figure 16 of Hirose et al. 2009b). In fact, near the midplane, the magnetic pressure generally rises with height, while the radiation pressure falls. For example, in the advective region of 0519b (|z| ≲ 1.5H), $d\ln P_{\rm mag}\over d\ln P_{\rm rad}$ rises outward from ≃ −1.5 very near the midplane to a plateau at ≃ 0.5. Thus, rising magnetized fluid parcels eventually become hotter than their surroundings, at which point heat will be transferred outward to the surroundings by radiative diffusion.

Not surprisingly, regions of strong Poynting flux are even more strongly associated with high magnetic field, low mass density isobaric fluctuations, as shown in Figure 22. This figure is constructed in a fashion exactly analogous to Figure 21 illustrating radiation advection. Supporting our argument that the two are fundamentally the same, the patterns seen in this figure are very similar to those seen in the radiation advection figure. The principal contrast is that, unlike the radiation advection case, there is very little downward Poynting flux; it is almost exclusively upward. This follows, of course, from the fact that the upward-moving isobaric fluctuations have greater than average magnetic energy density, whereas they have smaller than average radiation energy density.

Zoom In Zoom Out Reset image size

6.3. Characteristic Speed of Advective Motions

Because the outward advection fluxes are associated with regions of enhanced magnetic field, we can define a horizontally averaged vertical advection speed by using the z-component of the Poynting vector S = c(E × B)/(4π),

$$\begin{eqnarray} v_{\rm adv} &\equiv& 4\pi \int _{-L_x/2}^{L_x/2} dx \int _{-L_y/2}^{L_y/2} dy S_z(x,y,z)\Big/{}{}\int _{-L_x/2}^{L_x/2} dx\nonumber\\ &&\times \int _{-L_y/2}^{L_y/2} dy B^2(x,y,z). \end{eqnarray} \tag{ 34 }$$

We compare the time-averaged vertical profile of this speed in simulation 0519b to the time-averaged vertical profiles of the thermal sound speed c_t and the Alfvén speed in Figure 23.¹² The advection speed is several percent of the sound speed in the near-midplane region where advection is most important, but rises outward, reaching Mach numbers a little less than unity in the outer layers of the disk. This outward increase in v_adv is consistent with the characteristic buoyant rise time of 5–10 orbits for magnetic features: because H ∼ c_t/Ω, such a characteristic time implies a mean upward Mach number ∼0.1–0.2.

Zoom In Zoom Out Reset image size

Figure 23 also shows the characteristic one-dimensional speed of radiation diffusion v_diff over a length scale of ℓ = 0.1H, the characteristic size of the buoyant fluctuations (Figure 12). We define this speed as this length scale divided by the radiation diffusion time, i.e.,

$$\begin{equation} v_{\rm diff}\equiv {c\over 3\ell \bar{\kappa }^{\rm R}\rho }. \end{equation} \tag{ 35 }$$

The diffusion speed increases rapidly away from the midplane as the density drops. It is comparable to the advection speed in the midplane regions and then greatly exceeds the advection speed at heights greater than ≃ H. This coincides with the heights above which radiative diffusion finally becomes dominant over radiation advection in this simulation (bottom panel of Figure 1).

The average advection velocity profile also allows us to estimate the work done by the plasma in driving radiation advection. Figure 23 shows that at z = H, v_adv ≃ 10⁷ cm s⁻¹ in simulation 0519b. The average radiation pressure at that height is ≃ 8 × 10¹⁴ dyne cm⁻². Hence the mechanical pumping work associated with radiation advection at this height is ≃ −(E/3)dv_adv/dz ∼ −(E/3)v_adv/H ∼ −2 × 10¹⁵ erg cm⁻³ s⁻¹. This agrees very well with the actual mechanical pumping in Figure 6 (lower gray dashed curve). We conclude that the vertical expansion motions associated with radiation advection are largely responsible for the mechanical pumping work done by the plasma.

We have now answered a number of the questions posed at the beginning of this paper on the thermodynamics of radiation pressure-dominated disks. Central to our discussion is the existence of high magnetic energy density, low gas density regions that manifest within the turbulence as nonlinear slow magnetosonic modes. We named these regions "isobaric fluctuations" because the total pressure perturbation within them (thermal plus magnetic) is quite small. Higher than average magnetic pressure is associated with lower than average thermal pressure. The isobaric fluctuations are nonlinear, in that they are associated with order-unity magnetic intermittency. The magnetic energy density in these fluctuations is a few times the volume mean, and they constitute a significant fraction of the magnetic energy in the turbulence. On the other hand, these regions are rarely of sufficient size and coherence that they can be thought of as "flux tubes."

We have elucidated the nature and role of radiative damping. To supply the energy that is ultimately dissipated by photon diffusion, thermal pressure does work on the fluid in association with compressive fluctuations. This can be distinguished from work done by the fluid in driving vertical mechanical motions by subtracting off appropriate horizontal averages in pressure and velocity (Equation (20)). The compressive fluctuations responsible for radiative damping are a combination of isobaric fluctuations and acoustic waves. In the midplane regions of our most radiation-dominated simulation, we find that the dissipation associated with radiative damping is roughly 65% due to isobaric modes and 35% to acoustic waves. At higher altitude, isobaric fluctuations dominate the dissipation associated with radiative damping.

It is noteworthy that the time-averaged spatial distribution of radiative damping is similar in shape to the time-averaged spatial distribution of magnetic and kinetic energy dissipation (Figures 5 and 6). Similarly, the box-averaged radiative damping is well correlated with the magnetic and kinetic energy dissipation rates as a function of time: averaging over the short timescale fluctuations associated with acoustic waves in the radiative damping (upper dashed gray) curve of Figure 9 makes this history's relation to the magnetic (blue) and kinetic energy (green) dissipation curves clear. In addition, as we have already noted, a significant fraction of the total magnetic energy is found in regions where the magnetic energy density is of order twice the mean, places where strong isobaric fluctuations have led to field concentration. We therefore suspect that isobaric fluctuations play a significant role in the grid-scale magnetic and kinetic energy dissipation. In the magnetic case, this is perhaps not surprising, as the large concentrations of magnetic field must be associated with sharp gradients of magnetic field, as well as high electric current densities.

Surprisingly, we have also seen that when the ratio of radiation pressure to gas pressure is large enough, the total dissipation rate due to radiative damping can become comparable to the grid-scale loss of magnetic energy. If this trend continues to still larger radiation pressure, it will mean that we will have the unique advantage of being able to resolve the principal scale of dissipation in astrophysical MHD turbulence. That will finesse a great many of the uncertainties and difficulties that stand in the way of a clear interpretation of MHD turbulence in many other contexts.

Indeed, there is growing evidence that the saturated state of MRI turbulence may depend on the microscopic viscosity and resistivity of the plasma (Fromang et al. 2007; Lesur & Longaretti 2007; Simon & Hawley 2009; Davis et al. 2010). In particular, much attention has been paid to the sustainability and properties of the turbulence for different magnetic Prandtl numbers (the dimensionless ratio of viscosity to resistivity), because this sets the relative sizes of the velocity and magnetic dissipation scales. The inner, radiation-dominated regions of optically thick, geometrically thin accretion disks around black holes generally have magnetic Prandtl number greater than unity if the viscosity and resistivity are dominated by Coulomb scattering (Balbus & Henri 2008). In other words, the viscous scale, where velocity fluctuations are damped, exceeds the resistive scale where magnetic fluctuations are damped. Insofar as radiative damping also suppresses velocity fluctuations, albeit only those which are compressive, it reinforces this ordering: velocity fluctuations are damped on larger spatial scales than magnetic fluctuations in radiation-dominated accretion disks.

Perhaps the most important question that we posed at the beginning of this paper was how a radiation-dominated disk handles a dissipation rate that exceeds the critical value of Q^⋆ = cΩ²/κ associated with hydrostatic and radiative equilibrium. As shown in Figures 5 and 6, such an excess always occurs in the midplane regions. We find that the radiative diffusion flux is indeed constrained by hydrostatic equilibrium to have a divergence equal to Q^⋆, and excess energy flux is carried outward by radiative advection, fluid motion carrying trapped photons. A smaller (by ∼β⁻¹) energy flux is conveyed electromagnetically by the same fluid motions. At higher altitudes, the vertical energy flux becomes almost entirely diffusive radiation flux, simply because it becomes progressively harder for small fluid elements to trap photons as the density falls.

The existence of substantial vertical energy flux carried by something other than radiative diffusion or Poynting flux is a significant addition to the repertory of accretion disk physics. This is the first time that advected radiation flux has been shown to be important in disks.

We emphasize that these advected energy fluxes are not classical convection, as the horizontal and time-averaged entropy profiles are stably stratified. This is also true of the horizontal and time-averaged magnetic profiles, at least in the midplane regions. These advected fluxes are also not in any way related to the old argument that radiation-dominated disks have constant density and are therefore convectively unstable. As we discussed in the introduction, this argument rested on the assumption that the dissipation per unit mass is constant. This is simply not true, nor is the density constant with height.

Local buoyancy is the mechanism responsible for driving both radiation advection and upward Poynting flux, and this buoyancy arises from the low densities of the nonlinear, highly magnetized isobaric fluctuations. Because the lower than average thermal pressure within these fluctuations is mostly radiation pressure, these buoyant fluid parcels are cooler than their surroundings, at least initially. Nevertheless, they advect radiation internal energy outward because their magnetic pressure makes them buoyant. The work done by the plasma in driving this local vertical expansion explains the mechanical pumping work that we separated from the radiative damping in the total PdV work.

Although these buoyant regions play an important role in the energy budget in very radiation dominated disks, they are also important for a different reason in disks with weaker radiation content. The same isobaric fluctuations can also explain the characteristic upwelling of magnetic flux seen in all stratified disk simulations. Radiation forces per se are unnecessary: the only prerequisite is that there be thermal pressure fluctuations accompanied by proportionate density fluctuations that are balanced by opposite sign fluctuations in the magnetic pressure. Data from our earlier gas-dominated simulations show essentially the same dynamics as the data from radiation-dominated simulations that have played the primary role in this paper's analysis.

The quasi-periodic buildup of magnetic flux in the midplane regions is also associated with azimuthal field reversals, clearly indicating some sort of dynamo activity (Brandenburg et al. 1995; Johansen et al. 2009; Davis et al. 2010; Gressel 2010; O'Neill et al. 2010; Shi et al. 2010). The fact that the resulting upwelling motions are energetically significant is the first example that we know of in astrophysics of a magnetic dynamo that is important for the global energetics of the medium. For example, the solar dynamo is in part driven by the outward energy transport associated with convection in the Sun. In the case of a radiation-dominated accretion disk, the dynamo creates buoyant high magnetic field regions which then play a significant role in energy transport by advecting radiation (and magnetic field) outward. We have seen that this radiation advection is just as important as radiative diffusion in the midplane regions for the most radiation-dominated simulation we have conducted so far.

A final outstanding question is how the turbulence "knows" that it needs to produce a heating rate that is of order the critical value Q^⋆ required by hydrostatic and radiative equilibrium. We have seen that, in fact, when the overall radiation to gas pressure ratio is high, the dissipation rate generally exceeds Q^⋆ in the midplane regions, and the excess heat is carried outward by radiation advection, not radiative diffusion. But is this radiation advection therefore somehow regulated to match the requirements of overall thermal equilibrium? We suspect that the answer is no, at least not directly.

The nonlinear isobaric fluctuations are an inherent feature of intermittency in MRI turbulence and are related (somehow) to the dynamo process producing quasi-periodic reversals in the azimuthal field. They are inherently buoyant, and therefore the resulting rate of radiation advection is simply proportional to how much radiation energy density is present in the fluid to be carried upward with the magnetic field. As we just discussed, these fluctuations do appear to be associated with grid-scale dissipation that helps determine the local radiation energy density, but that is a slow process. Dissipation, which is still mostly magnetic in nature, has to occur over timescales comparable to the thermal time to build up enough heat to significantly change the radiation energy density, simply because the volume-averaged magnetic energy density is so small.¹³

Instead, any small excess radiation energy density that appears because the average dissipation exceeds Q^⋆ results in an excess vertical diffusive radiation flux. The resulting excess radiation pressure force quickly (on the fast sound crossing timescale) produces a vertical expansion of the medium. As a consequence, work is done by the fluid and this lowers the radiation energy density, restoring hydrostatic equilibrium. This mechanism may in part be responsible for the excitation of the standing vertical acoustic waves. It may very well be that hydrostatic equilibrium cannot be so maintained if the turbulent dissipation becomes significantly stronger. Outflows may appear at levels of radiation dominance we have not yet explored in our simulations.

Ultimately, the mechanism that sets the midplane dissipation to be in excess of Q^⋆, and the level of ambient radiation energy density that is then advected outward by buoyant magnetic field, is the level of turbulence in the medium. However, what sets the saturation level of MRI turbulence is still an open question.

Finally, we remark that photon bubbles have been proposed as an additional piece of physics that may alter the vertical transport of energy at high luminosities (Begelman 2006). The characteristic length scale of these instabilities is ∼c²_i/g, where c_i is the isothermal gas sound speed and g = Ω²|z| is the acceleration due to gravity (Gammie 1998; Begelman 2001; Blaes & Socrates 2003; Turner et al. 2005). In simulation 0519b, this is less than the size of our grid zones at distances more than ≃ 0.5H away from the midplane. Clearly, these instabilities cannot be resolved in our simulations, and their possible impact on the thermodynamics remains a topic for future study.

The thermodynamics of a radiation-dominated disk differs significantly from that envisaged in standard static models of accretion disks. In addition to microscopic dissipation associated with the MRI turbulent cascade, radiative damping of compressive fluctuations plays an increasingly significant role as the radiation to gas pressure ratio increases. The compressive fluctuations associated with this damping consist of fast magnetosonic waves as well as nonlinear isobaric fluctuations arising from the turbulence itself. This damping is energetically significant: on the order of tens of percent of the overall dissipation at the highest levels of radiation pressure support that we have simulated. It is also numerically resolvable, in contrast to microscopic resistive and viscous dissipation.

Buoyancy of the highly magnetized, low density isobaric fluctuations is responsible in part for the butterfly diagram that is always seen in simulations of MRI turbulence with vertical gravity. This buoyancy is intrinsically three dimensional in nature; even though the horizontally averaged structure is very stably stratified in the midplane regions, localized concentrations of magnetic field generated by the turbulent dynamo are underdense and still produce outward advection of magnetic field. When the plasma is radiation dominated, these buoyant motions become significant for the overall energetics of the plasma. Outward advection of photons becomes comparable to radiative diffusion, and the associated vertical expansion work must be included in balancing vertical heat transport and dissipation.

This work was supported by NSF grants AST-0707624 and AST-0507455, and by a Grant-in-Aid for Scientific Research (No. 20340040) from the Ministry of Education, Culture, Sports, Science and Technology of Japan. Numerical computations were carried out partially on the Cray XT4 at the Center for Computational Astrophysics at the National Astronomical Observatory of Japan, and on the SX8 at the Yukawa Institute for Theoretical Physics at Kyoto University. Our work has benefited from conversations with numerous people, including Steve Balbus, Shane Davis, John Hawley, Tobias Heinemann, Ken Henisey, Gordon Ogilvie, Chris Reynolds, Prateek Sharma, Jim Stone, Ted Tao, Neal Turner, and Ethan Vishniac. We also thank the referee, Jeremy Goodman, for detailed constructive criticism that greatly improved the paper. Part of this work was completed while O.B. was visiting the Kavli Institute for Astronomy and Astrophysics in Beijing.

In this appendix, we outline the derivations of the damping times and vertical energy fluxes of the standing vertical acoustic waves discussed in Section 6.1.

We consider a vertically stratified equilibrium with a background shear flow v_y = −3Ωx/2. Apart from v_y, all other fluid quantities are constant in the horizontal direction. We assume that the radiation and gas exchange heat sufficiently rapidly that T_rad = T everywhere. We also take the medium to be sufficiently optically thick that the radiation pressure tensor ${\sf P}$ only has diagonal elements E/3, and radiation transport is diffusive with $\bar{\kappa }^{\rm R}\simeq \kappa _{\rm es}={\rm constant}$. These thermodynamic assumptions are excellent approximations in the regions of the simulations where the acoustic waves are evident. Finally, we approximate the equilibrium magnetic field as being purely azimuthal and having only vertical gradients.

The equilibrium vertical structure is given by hydrostatic equilibrium,

$$\begin{equation} {d\over dz}\left(P_{\rm therm}+{B^2\over 8\pi }\right)=-\rho \Omega ^2z, \end{equation} \tag{ A1 }$$

radiative equilibrium,

$$\begin{equation} {dF\over dz}=Q \end{equation} \tag{ A2 }$$

and diffusive radiation transport,

$$\begin{equation} F=-{c\over 3\kappa _{\rm es}\rho }{{\it dE}\over dz}. \end{equation} \tag{ A3 }$$

Here, P_therm ≡ p + E/3 is the total thermal pressure. Note that we are including turbulence only insofar as it contributes energy dissipation Q.

We now consider vertically propagating, longitudinal waves on this equilibrium. The Eulerian perturbed fluid velocity δv(z, t) = ∂ξ/∂t, where ξ is the vertical Lagrangian displacement of each fluid element. The linearized continuity equation is

$$\begin{equation} \Delta \rho =-\rho {\partial \xi \over \partial z}, \end{equation} \tag{ A4 }$$

where, here and elsewhere, Δ ≡ δ + ξ(d/dz) refers to Lagrangian perturbations. The linearized vertical momentum equation is

$$\begin{equation} \rho {\partial \delta v\over \partial t}=-{\partial \delta P_{\rm therm} \over \partial z}-\delta \rho \Omega ^2 z-{1\over 4\pi }{\partial \over \partial z} (B\delta B). \end{equation} \tag{ A5 }$$

The artificial viscosity term has been lost as it is nonlinear in the perturbation velocity (Stone & Norman 1992). The linearized flux-freezing equation is

$$\begin{equation} \delta B=-{\partial \over \partial z}(B\xi). \end{equation} \tag{ A6 }$$

Finally, the linearized total internal energy equation may be written as

$$\begin{equation} {\partial \over \partial t}\left({\Delta P_{\rm therm}\over P_{\rm therm}}\right) -\Gamma _1{\partial \over \partial t}\left({\Delta \rho \over \rho }\right)= {\Gamma _3-1\over P_{\rm therm}}\left(\delta \tilde{Q}- {\partial \delta F\over \partial z} \right), \end{equation} \tag{ A7 }$$

where Γ₁ and Γ₃ are the usual generalized adiabatic exponents for a mixture of gas and radiation at the same temperature (Chandrasekhar 1967), both of which are close to 4/3 for our radiation-dominated simulations. This equation may be integrated in time to give

$$\begin{equation} {\Delta P_{\rm therm}\over P_{\rm therm}} -{\Gamma _1\Delta \rho \over \rho }= {\Gamma _3-1\over P_{\rm therm}}\int ^t\left(\delta \tilde{Q}- {\partial \delta F\over \partial z}\right)dt^\prime. \end{equation} \tag{ A8 }$$

Assuming the mode periods are short enough that there is little energy exchange between the modes and the turbulence, little turbulent dissipative heating, and little radiative damping, we can for the moment set $\delta \tilde{Q}=\delta F=0$ in the equations. Assuming a time dependence ∝exp (− iωt), Equations (A4)–(A6) and (A8) may then be combined to give a single equation in δv:

$$\begin{equation} {\partial \over \partial z}\left[\left(\Gamma _1P_{\rm therm}+{B^2\over 4\pi }\right) {\partial \delta v\over \partial z}\right]+\rho (\omega ^2-\Omega ^2)\delta v=0. \end{equation} \tag{ A9 }$$

The vertical epicyclic mode is the simplest general solution to this equation, with ω = Ω and δv being spatially constant. The breathing mode is the next simplest. Provided Γ₁ is spatially constant (which is true if gas or radiation dominates the thermal pressure), and the magnetic energy density is much less than the thermal pressure (an assumption that typically fails only in the low density surface layers), then the hydrostatic equilibrium Equation (A1) guarantees that ω = (1 + Γ₁)^1/2Ω and δv∝z satisfy Equation (A9).

Like the breathing mode, higher order modes are acoustic in nature, but their frequencies and eigenfunctions depend more sensitively on the details of the equilibrium vertical structure (Okazaki et al. 1987; Lubow & Pringle 1993; Blaes et al. 2006).

We now consider radiative damping of these modes. Multiplying Equation (A5) by δv, and then combining with Equations (A4), (A6), and (A8), we finally obtain after some algebra the following wave energy conservation equation:

$$\begin{eqnarray} &&{\partial \over \partial t}\left[{1\over 2}\rho \delta v^2+{1\over 2}\Gamma _1 P_{\rm therm} \left({\Delta \rho \over \rho }\right)^2+{1\over 2}\rho \xi ^2\Omega ^2+{(\Delta B)^2 \over 8\pi }\right] \nonumber\\ &&\qquad\quad\quad+{\partial \over \partial z}\left[\delta v\Delta P_{\rm therm}+\delta v \Delta \left({B^2\over 8\pi }\right)\right]\nonumber\\ &&\qquad\quad =-{(\Gamma _3-1)\over \rho }{\partial \Delta \rho \over \partial t}\int ^t\left(\delta \tilde{Q}- {\partial \delta F\over \partial z}\right)dt^\prime. \end{eqnarray} \tag{ A10 }$$

The right-hand side of Equation (A10) represents sources and sinks of wave energy due to perturbed turbulent dissipation and radiative losses. True microscopic viscosity and resistivity would also contribute additional terms on the right-hand side arising from additional terms in the perturbed momentum and flux-freezing equations, but our numerical dissipation scheme has no such terms. In reality they would be present, but likely even more important than these would be terms representing turbulent fluctuations, including turbulent viscosity and resistivity, which we have ignored here. However, we believe the damping to be dominated by radiative diffusion.

Assuming that the damping rate is small, we can estimate it by computing a work integral for the modes (e.g., Cox 1980). Approximating all the perturbed quantities on the right-hand side as being perfectly periodic at the mode period Π = 2π/ω_R, where ω_R ≡ Re(ω), we can time average Equation (A10) over this mode period and integrate the right-hand side by parts. Then integrating over all height, assuming negligible wave energy fluxes leaving the top and bottom boundaries, we finally obtain for the exponential damping rate Γ ≡ Im(ω) of the mode amplitude

$$\begin{equation} \Gamma =-{1\over 2\Pi {\newmathcal E}_{\rm mode}}\int dz(\Gamma _3-1) \int _0^\Pi dt{\Delta \rho \over \rho } \left(\delta \tilde{Q}-{\partial \delta F\over \partial z}\right), \end{equation} \tag{ A11 }$$

where

$$\begin{eqnarray} {\newmathcal E}_{\rm mode}&\equiv& \int dz \bigg[{1\over 2}\rho \delta v^2+{1\over 2}\Gamma _1 P_{\rm therm} \bigg({\Delta \rho \over \rho }\bigg)^2\nonumber\\ &&+{1\over 2}\rho \xi ^2\Omega ^2+{(\Delta B)^2 \over 8\pi }\bigg] \end{eqnarray} \tag{ A12 }$$

is the mode energy.

Because δv and therefore ξ are spatially constant for the epicyclic mode, Equation (A4) implies that Δρ = 0, i.e., fluid elements preserve their density under this mode. The work integral therefore vanishes, and there is no damping of this mode. This is as one would expect: the only way to change this mode is to exert an external force on the system or eject momentum through the vertical boundaries.

The breathing mode and higher order acoustic modes are damped, however. The linearized radiative diffusion equation is

$$\begin{equation} \delta F=-F{\delta \rho \over \rho }-{c\over 3\kappa _{\rm es}\rho } {\partial \delta E\over \partial z}. \end{equation} \tag{ A13 }$$

If we consider a radiation-dominated equilibrium, and neglect gas and magnetic pressure contributions to the mode dynamics, then Equations (A5) and (A13) give a simple relationship between the perturbed radiative flux and the perturbed velocity,

$$\begin{equation} \delta F={c\over \kappa _{\rm es}}{\partial \delta v\over \partial t}. \end{equation} \tag{ A14 }$$

Setting δv = A_b(z/H_P)cos ω_Rt for the breathing mode, where A_b is the velocity amplitude at one pressure scale height H_P and ω_R = (7/3)^1/2Ω, we find a radiative damping rate of

$$\begin{equation} \Gamma _{\rm rad}={c\Omega ^2\over 4\kappa _{\rm es}H_PE_0}\int dz, \end{equation} \tag{ A15 }$$

where E₀ is the midplane radiation energy density, and the integral goes over the entire equilibrium, or at least that portion which is optically thick. Clearly, this integral can be set equal to the pressure scale height H_P times some numerical constant. For example, for an n = 3 isentropic polytrope, we obtain

$$\begin{equation} \Gamma _{\rm rad}={315\over 256}\left({c\Omega ^2\over \kappa _{\rm es} E_0} \right)\sim {c\Omega ^2\over \kappa _{\rm es}E_0}. \end{equation} \tag{ A16 }$$

The existence of this damping produces nonzero time averages in the wave energy fluxes because it changes the phase difference between the velocity perturbation δv and the pressure perturbations ΔP_therm and Δ(B²/8π). In the absence of an exact solution for the non-adiabatic eigenfunctions, we compute these pressure perturbations using the adiabatic equations. Setting δv = A_b(z/H_P)e^−Γtcos ω_Rt, then the instantaneous energy fluxes are

$$\begin{eqnarray} \delta v\Delta P_{\rm therm}&\simeq& \delta v{\Delta E\over 3}= -{4EzA_{\rm b}^2\over 9\omega _{\rm R}H_P^2} e^{-2\Gamma t}\nonumber\\ &&\times\left(\cos \omega _{\rm R}t\sin \omega _{\rm R}t- {\Gamma \over \omega _{\rm R}}\cos ^2\omega _{\rm R}t\right)+ {\newmathcal O}\left({\Gamma ^2\over \omega _{\rm R}^2}\right)\nonumber\\ \end{eqnarray} \tag{ A17 }$$

and

$$\begin{eqnarray} \delta v\Delta \left({B^2\over 8\pi }\right)&=&{-}{B^2zA_{\rm b}^2\over 4\pi \omega _{\rm R}H_P^2} e^{-2\Gamma t}\nonumber\\ &&\times\left(\cos \omega _{\rm R}t\sin \omega _{\rm R}t- {\Gamma \over \omega _{\rm R}}\cos ^2\omega _{\rm R}t\right)\nonumber\\ &&+ {\newmathcal O}\left({\Gamma ^2\over \omega _{\rm R}^2}\right). \end{eqnarray} \tag{ A18 }$$

The first term in each of these expressions periodically changes sign, and a time average of this term over one mode period results in positive or negative values depending on the phase of the time-interval chosen. This term therefore produces no secular energy flux. The second term does, and the resulting time-averaged fluxes are, to lowest order in Γ/ω_R,

$$\begin{equation} \langle \delta v\Delta P_{\rm therm}\rangle ={2Ez\Gamma A_{\rm b}^2\over 21\Omega ^2H_P^2} \end{equation} \tag{ A19 }$$

and

$$\begin{equation} \left<\delta v\Delta \left({B^2\over 8\pi }\right)\right>={3B^2z\Gamma A_{\rm b}^2\over 56\pi \Omega ^2H_P^2}. \end{equation} \tag{ A20 }$$

Note that these two equations are completely independent of whether the decay rate Γ was due to radiation damping or some other process.