---
abstract: 'Early in the study of viscous accretion disks it was realized that energy transfers from distant sources must be important, not least because the flow at the disk midplane in the bulk of the disk is likely outwards, out of the gravitational potential well. If the source of the viscosity is powered by accretion, such as in the case of the magneto-rotational instability, such distant energy sources must lie in the innermost regions of the disk, where accretion occurs even at the midplane. We argue here that modulations in this energy supply can alter the accretion rate on dynamical, rather than far longer viscous, time scales. This means that both the steady state value of and fluctuations in the inner disk’s accretion rate, depending on the details of the inner boundary condition and occurring on the inner disk’s rapid evolution time, can affect the outer disk. This is particularly interesting because observations have shown that disk accretion is not steady (e.g. EX Lupi type objects). We also note that the power supplied to shearing boxes is set by the boxes themselves rather than the physical energy fluxes in a global disk. That is, their saturated magnetic field is not subject to the full set of energy constraints present in an actual disk. Our analysis suggests that large scale radial transport of energy has a critical impact on the evolution and variability of accretion disks.'
author:
- |
    Alexander Hubbard, Colin P. McNally, Jeffrey S. Oishi,\
    Wladimir Lyra and Mordecai-Mark Mac Low
title: Radial Stresses and Energy Transport in Accretion Disks
---

Introduction
============

Accretion disks surround many objects in the universe, from young stars, to compact binary companions, to supermassive black holes. Their observed accretion flows require outwards transport of angular momentum. The magnetic coupling of regions at different radii in accretion disks ranks high among the usual suspects for angular momentum transport mechanisms. This transport can occur either radially through the disk, such as with the magnetorotational instability [MRI, @1998RvMP...70....1B], or vertically through a disk wind [@2000prpl.conf..759K]. Understanding vertical angular momentum transport through disk winds depends on understanding the disk’s vertical boundaries. Similarly, understanding radial angular momentum transport requires understanding the radial boundary conditions. In this paper, we argue that this is true not only near to but also far from those boundaries. While that may seem a trivial statement, in practice it means that even distant boundary conditions must be taken into account by numerical simulations intended to study the magnitude and direction of any accretion flows.

The common assumption of local control implicitly places any transition between accretion and decretion far from any region explicitly considered. Simulating regions far from any physically distinct radius further allows the use of the shearing box approximation, which goes so far as to erase the distinction between inwards and outwards directions. @1999ApJ...521..650B [hereafter BP99] found that the mean flow dynamics of MRI-driven MHD turbulence were a local phenomenon, justifying the use of alpha-disk models and shearing boxes.

It has long been known that steady state, vertically integrated, viscous, regions of accretion disks require an outwards energy flux in the standard theory [@1974MNRAS.168..603L][^1]. That theory defines $r_b$ as the position of the boundary layer where the differential rotation goes to zero [@2002apa..book.....F]. Protostars rotate below the Keplerian frequency, so $r_b$ must exist for disks that extend to the surface of their protostar, although it may not exist in strongly magnetized disks [@2002apa..book.....F]. It was shown early in the study of viscous disks that dissipation at radii $r > 9r_b/4$ thermalizes significantly more energy than is released by the accretion flow there [@1974MNRAS.168..603L; @1998RvMP...70....1B; @2002apa..book.....F], with the balance made up by an energy flux built up within $ 9r_b/4$, and depleted over the entire outer disk.

Viscous disk theory predicts that disks decrete, rather than accrete, in regions where the stresses that drive accretion and angular momentum transport fall off faster with radius than $r^{-2}$ [@1974MNRAS.168..603L BP99]. In standard, self-similar disk models, this decretion condition is easily satisfied in disk midplanes where the pressure decreases rapidly with radius. It takes vertical integration through viscous disks to reliably find net accretion. This has led to the concept of meridional circulation, in which a decreting midplane is more than balanced by accreting surface layers . Such behavior has been demonstrated in numerical models with a fixed anomalous viscosity defined by a constant Shakura-Sunyaev $\alpha$ parameter .

These layered accretion structures arise in uniformly viscous disks because the disk thickness increases with radius in such disks. While the gas density and pressure in the midplane decrease with radius, they instead increase with cylindrical radius far enough above the midplane that the disk surface is encountered at finite radius. This occurs at increasing radius with height because the tilted disk surface reaches higher altitudes at larger radii. This transition in behavior is lost in slab geometries with radially constant scale height such as the shearing sheet approximation. In classical viscous disk theory the stress remains proportional to the gas pressure, so there exists a dividing height between decretion and accretion that occurs close to the surface, where the radial pressure gradient goes from falling off faster to falling off slower than $r^{-2}$ [e.g. @1994ApJ...423..736R].

Decreting regions of disks must be net importers of energy, since decretion boosts material out of a potential well. These regions still dissipate energy locally through non-ideal effects such as viscosity or resistivity, just as accretion regions do [@1974MNRAS.168..603L], but that just adds to the energy deficit that must be balanced by importing energy from inner regions. However, the source of that energy must be determined by the nature of the stresses driving accretion and decretion. If those stresses are powered by the accretion process itself and occur in the horizontal $r$-$\phi$ plane, as expected for magnetorational instability [@1991ApJ...376..214B; @1998RvMP...70....1B], We argue that this means that the midplane outflow of meridional circulation is powered not by the accreting upper layers, but rather by the inner edge of the midplane itself. Similar radial energy flows occur at every altitude, with the accreting inner edge of the disk at each altitude forming the accreting upper layer described by the meridional circulation picture.

BP99 argued that the mean flow dynamics of MRI-driven, MHD turbulence can be treated as a local phenomenon because the dynamical equations describing the flow are themselves local: the time derivatives at a specific spatial point of quantities such as the mean velocity or the mean magnetic field depend only on their values and spatial derivatives at that point. However, even distant regions of a disk will communicate on long enough time scales. This means that locality is more accurately treated as a time scale condition: on time scales shorter than some critical value, the disk can be treated as local, but it cannot be treated as local for longer time scales.

We show that the strength of the radial energy flows required to maintain the disk structure determine that critical time scale. The magnetic fields that provide the Maxwell stresses in MRI active disks have only a very small energy density, comparable to the accretion or decretion power integrated only over a few orbits. Order unity variations in the energy fluxes due to variability in the inner disk can therefore, everywhere within a few local orbits, cause order unity variations in the stresses, and hence the local accretion or decretion rates. This motivates us to extend the time dependent analysis in BP99 by including the time derivative of the stresses.

In Section \[sec\_angmom\] we revisit the basic equations for viscous disks. In Section \[Sec\_Fluxes\] we discuss the energy fluxes present and estimate the dynamical time-scale associated with the stresses. In Section \[sec:shear-box-appr\] we delve more deeply into the energy flux considerations in the context of the shearing sheet approximation, and discuss the difference between the case where stresses are powered by accretion energy (which includes most cases of turbulent viscosity) vs. the cases where the stresses are not powered by accretion energy (including many cases of microphysical viscosity). In Section \[sec:simul-accr-disks\] we discuss links to numerical simulations, both existing results and suggestions for future design and diagnostics. In Section \[sec:BP99\] we place our results in context as an extension of the analysis of BP99, and we conclude in Section \[sec\_conclusions\].

Angular Momentum Transport {#sec_angmom}
==========================

Azimuthal Lorentz Forces {#sec:azim-lorentz-forc}
------------------------

The action of the magnetic field on the gas to transfer angular momentum occurs through azimuthal Lorentz forces ${\mbox{\boldmath $F$} {}}_{L,\phi} ={\mbox{\boldmath $J$} {}}\times {\mbox{\boldmath $B$} {}}|_{\phi}$. In cylindrical coordinates, the magnetic field is $${\mbox{\boldmath $B$} {}}= B_r {{\bm{\hat{e}}_r}}+ B_\phi {{\bm{\hat{e}}_\phi}}+ B_z {{\bm{\hat{e}}_z}},$$ and the current density $$\begin{aligned}
\mu_0{\mbox{\boldmath $J$} {}}=& {\mbox{\boldmath $\nabla$} {}}\times {\mbox{\boldmath $B$} {}},  \\
 =&\left[ \frac{1}{r} \partial_\phi B_z -\partial_z B_\phi\right] {{\bm{\hat{e}}_r}}\nonumber\\
 & +\left[\frac{\partial B_r}{\partial z} - \frac{\partial B_z}{\partial r} \right] {{\bm{\hat{e}}_\phi}}\nonumber\\ 
 & + \left[\frac{1}{r}\partial_r(r B_\phi) -\frac{1}{r}\partial_\phi B_r \right] {{\bm{\hat{e}}_z}}.\end{aligned}$$ The azimuthal component of the Lorentz force can be usefully expanded and rearranged: $$\begin{aligned}
\mu_0 F_{L,\phi} =& \mu_0 \left( J_z B_r - J_r B_z\right) \label{eq:lorentz-terms}\\
= & \frac{1}{r}\left[\partial_r (r B_\phi)\right] B_r  - \frac{1}{r} (\partial_\phi B_r) B_r \nonumber  \\
& \qquad  - \frac{1}{r }(\partial_\phi B_z )B_z + (\partial_z B_\phi) B_z\\ 
=& -\frac{1}{2 r} \partial_\phi B_r^2 -\frac{1}{2 r} \partial_\phi B_z^2 +\frac{1}{r^2} \partial_r (r^2B_\phi B_r)  \nonumber\\
& \quad   -\frac{B_\phi}{r} \partial_r (r B_r)  + \partial_z(B_\phi B_z) - B_\phi \partial_z  B_z \\
=& -\frac{1}{2r} \partial_\phi \left( B_r^2 + B_z^2 \right) + \frac{1}{r^2}\partial_r (r^2 B_\phi B_r) \nonumber\\
&  + \partial_z (B_\phi B_z)+B_\phi\left[ -\frac{1}{r} \partial_r (r B_r) - \partial_z B_z \right].  \label{Reduced_Max_0}\end{aligned}$$ We can now invoke ${\mbox{\boldmath $\nabla$} {}}\cdot {{\bm B}}=0$ to clarify the component terms of the azimuthal Lorentz force: $$\begin{aligned}
\mu_0 F_{L,\phi}&= \underbrace{\frac{1}{r^2} \partial_r (r^2 B_\phi B_r)}_\text{Radial Stress} + \underbrace{\partial_z (B_\phi B_z)}_\text{Vertical Stress} \nonumber \\
 & \qquad \qquad \quad + \underbrace{\frac{1}{2r} \partial_\phi \left( B_\phi^2 - B_r^2 -B_z^2\right)}_\text{Anisotropic magnetic pressure}.
\label{Reduced_Max}\end{aligned}$$ The pressure term azimuthally averages to zero due to periodicity .

Accretion Stresses {#sec:accretion-stresses}
------------------

The components of Equation (\[Reduced\_Max\]) are derivatives of the radial and vertical components of the Maxwell stress tensor[^2] $$M_{ij} \equiv \mu_0^{-1}B_i B_j,$$ The evolution of angular momentum in disks includes contributions from both the magnetic Maxwell stress and the hydrodynamical Reynolds stress $$R_{ij} \equiv \rho v_i v_j,$$ often combined in a total stress tensor $$T_{ij}  \equiv R_{ij} - M_{ij}.\label{eq_Tdef}$$ We emphasize here that in the above equations we use the full, not fluctuating, velocities and magnetic fields. We can write the time evolution of the angular momentum density in terms of the stresses as $$\partial_t \langle \rho r v_\phi \rangle_\phi  + \frac
1r \partial_r  \langle r^2 T_{\phi r}\rangle_\phi + r \partial_z  \langle T_{\phi
  z}\rangle_\phi=0 \label{stress_0}.$$ where we have denoted averaging performed over the dimension $i$ with the notation $\langle \dots \rangle_{i}$. The azimuthal averaging here again eliminates the azimuthal pressure term. (Equation (\[stress\_0\]) is often vertically averaged as well.) We name stress terms such as $M_{\phi r}$ “horizontal” while terms such as $M_{\phi z}$ are labelled “vertical”.

### Alternate Definitions of the Stresses {#sec_rey_decomp}

It is common to decompose the velocities and magnetic field into mean (generally azimuthal averages) and fluctuating terms, marked here with overbars and primes respectively, while taking the azimuthal density fluctuations to vanish. In that case, Equation (\[eq\_Tdef\]) becomes $$\langle T_{ij}\rangle_\phi  = \rho {\overline}{v}_i{\overline}{v}_j - \mu_0^{-1} {\overline}{B}_i {\overline}{B}_j+\langle\rho v_i' v_j' \rangle_\phi -
\mu_0^{-1} \langle B_i'B_j'\rangle_\phi.\label{eq_Tdef2}$$

When moving from Equation (\[eq\_Tdef2\]) to Equation (\[stress\_0\]), two terms are often treated differently: $\rho {\overline}{v}_i {\overline}{v}_j$ is often separated from the stress, and written as the radial advection of angular momentum; while the term $\mu_0^{-1} {\overline}{B}_i{\overline}{B}_j$ is often neglected in determining the  $\alpha$ parameter.

### Independence of Torques from Density Gradients

Past studies, such as BP99, have chosen to perform a density weighted vertical and azimuthal average of Equation (\[stress\_0\]). The horizontal stress term is given by $$\int dz \left\langle \frac 1r \partial_r \left(r^2 T_{\phi r}\right)\right\rangle_{\triangle r,\phi}= \frac 1r \partial_r \left( r^2 \varSigma W_{r\phi}\right)
 \label{stress_1},$$ where $\triangle r$ represents a small interval in $r$, the disk surface density $\varSigma$ is given by $$\varSigma = \int dz \langle \rho
  \rangle_{\triangle r,\phi},$$ and $$W_{r \phi} \equiv \varSigma^{-1} \int dz \langle  T_{\phi r}
\rangle_{\triangle r, \phi} \label{W_def}.$$ For a truly viscous, azimuthally symmetric disk, the $T_{\phi r}$ term in Equation (\[stress\_0\]) is: $$\int dz\left(\frac1r \partial_r \langle r^2 T_{\phi r}\rangle_{\triangle r\phi}\right)=-\frac 1r \partial_r \left( r^3 \nu \varSigma \partial_r \varOmega \right),$$ where $\nu$ is the kinematic viscosity [@1974MNRAS.168..603L]. This averaging serves the purpose of showing that the horizontal Maxwell and Reynolds stresses act like a viscosity if we equate $$W_{r \phi} \sim -\nu r \partial_r \varOmega,$$ and induces no mathematical error. However, it can lead to confusion: a naive reading of Equation (\[stress\_1\]) suggests that a spatially varying surface density $\varSigma$ can create torques from a Maxwell stress with radially constant $r^2 M_{\phi r}$, in contradiction to Equation (\[stress\_0\]). Equation (\[W\_def\]) makes clear, however, that $W_{\phi r}$ is defined by dividing by $\varSigma$, so $W_{\phi r}$ and $\varSigma$ cannot be varied independently.

Condition for Accretion {#subsec:direction}
-----------------------

Equation (\[Reduced\_Max\]) tells us that the effect of the azimuthal force deriving from $M_{\phi r}$ depends on the sign of the product of $\partial_r (r^2 B_\phi B_r)$ with the angular velocity $\varOmega$. If $$\varOmega \partial_r (r^2 B_\phi B_r) > 0,$$ then the Lorentz force exerts a torque on the disk aligned with the orbital rotation, resulting in decretion, while if it is negative, the Lorentz force exerts a torque acting against the rotation, resulting in an accretion flow.

Driving an accretion flow then requires $$\label{Criterion}
\partial_r \left[r^2 (-M_{\phi r}) \right]>0$$ (from Eq. \[Reduced\_Max\]). This can be rephrased as the statement that accretion requires that the magnitude of $M_{\phi r}$ have a radial dependence shallower than $r^{-2}$.

Stresses Proportional to Pressure {#sec_T_propto_p}
---------------------------------

If the only physical quantities of the background disk are its density $\rho$ and sound speed $c_s$, then the the stresses must scale with the thermal pressure, the only combination of the parameters with the correct dimensions: $\rho_0 c_s^2$ (the standard case of an $\alpha$ disk, ). This means that $$M_{\phi r} \propto \rho_0 c_s^2 \propto H \varSigma  r^{-3}, \label{ZNF_scaling}$$ where we have used $\varOmega^2 \propto r^{-3}$. For such a region of a Keplerian disk to accrete, the radial dependence of $H \varSigma$ must exceed $r^{+1}$. Observations suggest that the surface density dependence for real accretion disks cannot be much shallower than $\varSigma \propto r^{-1}$ [@2003MNRAS.341..501S; @2007prpl.conf..555D; @2010ApJ...723.1241A], so to produce accretion the minimum disk flaring must be $H\propto r^2$. The standard disk model with constant opening angle and $\varSigma \propto r^{-1}$ has $H\propto r$, so it must have horizontal Maxwell stresses that drive decretion instead of accretion.

Of course, Equation (\[ZNF\_scaling\]) is local, and the scaling applies only to the midplane. Vertically integrating it leads to $$\int_z dz M_{\phi r} \propto H^2 \varSigma r^{-3},$$ which will generally fall off slower than $r^{-2}$, recovering the net accretion flow of viscous disk theory. Such a scenario is known [@2002ApJ...581.1344T] to have a meridional circulation with a decreting midplane as above counterbalanced by accreting surface layers.

### Stresses Not Proportional to Pressure

Especially in the case of magnetized accretion flows, more physics is expected, which provides additional parameters that can control the stresses. For example, the stresses in the MRI case likely depend on any net background vertical magnetic field, which has no reason to be proportional to the midplane pressure, and the stresses can certainly depend on non-ideal effects. Further, in the MRI case any estimation $M_{\phi r} \propto \rho c_s^2$ applies only in regions where the magnetic fields are locally generated, and simulations have shown that the upper layers of accretion are actually coronas powered by magnetic field generated closer to the midplane . This means that $M_{\phi r}$ will not be proportional to the thermal pressure precisely in the upper layers where standard meridional circulation models predict accretion flows, which may help explain why ’s MRI simulations saw both decreting midplanes and surface layers. Nonetheless, in the absence of dominant imposed physics, standard disk radial scalings imply decreting midplanes.

Energy flows {#Sec_Fluxes}
============

It is puzzling that decretion flows are predicted, and seen in numerical simulations. Where is the origin of the energy required to drive such a flow out of a gravitational potential? If a volume embedded in the disk is locally accreting, then the local loss of gravitational potential energy can provide the energy dissipated during turbulent angular momentum transport. If it is locally decreting, though, then some distant energy source must provide the energy for both the decretion flow and the associated turbulent dissipation.

Energy Transport by Poynting Flux {#sec_E_transport}
---------------------------------

Radial fluxes of energy (in the magnetic case, Poynting fluxes) are known to be dynamically significant in disks accreting through horizontal viscous stresses [@1974MNRAS.168..603L; @1998RvMP...70....1B; @2002apa..book.....F]. We can calculate the time evolution of the magnetic energy density: $$\partial_t \frac{B^2}{2 \mu_0} = {\mbox{\boldmath $B$} {}}\cdot \left[{\mbox{\boldmath $\nabla$} {}}\times \left( -{\mbox{\boldmath $E$} {}}\right)\right], \label{E0}$$ where the electric field ${\mbox{\boldmath $E$} {}}= -{\mbox{\boldmath $v$} {}}\times {\mbox{\boldmath $B$} {}}+ \eta {\mbox{\boldmath $J$} {}}$, with $\eta$ being the resistivity.

Noting that $$\begin{aligned}
{\mbox{\boldmath $B$} {}}\cdot \left[{\mbox{\boldmath $\nabla$} {}}\times \left( -{\mbox{\boldmath $E$} {}}\right)\right] &= -{\mbox{\boldmath $E$} {}}\cdot \left({\mbox{\boldmath $\nabla$} {}}\times {\mbox{\boldmath $B$} {}}\right) -{\mbox{\boldmath $\nabla$} {}}\cdot \left( {\mbox{\boldmath $E$} {}}\times {\mbox{\boldmath $B$} {}}\right)  \\
&=\left({\mbox{\boldmath $v$} {}}\times {\mbox{\boldmath $B$} {}}\right) \cdot {\mbox{\boldmath $J$} {}}-\eta {\mbox{\boldmath $J$} {}}^2 -{\mbox{\boldmath $\nabla$} {}}\cdot {\mbox{\boldmath $S$} {}},\end{aligned}$$ where ${\mbox{\boldmath $S$} {}}= \mu_0^{-1} {\mbox{\boldmath $E$} {}}\times {\mbox{\boldmath $B$} {}}$ is the Poynting flux, we can rewrite Equation (\[E0\]) as $$\partial_t \frac{B^2}{2 \mu_0} = -{\mbox{\boldmath $\nabla$} {}}\cdot {\mbox{\boldmath $S$} {}}- \eta {\mbox{\boldmath $J$} {}}^2 - {\mbox{\boldmath $F$} {}}_L \cdot {\mbox{\boldmath $v$} {}}, \label{E1}$$ where ${\mbox{\boldmath $F$} {}}_L \equiv {\mbox{\boldmath $J$} {}}\times {\mbox{\boldmath $B$} {}}$ is the Lorentz force. Equation (\[E1\]) states that, up to resistive terms, any energy taken from the kinetic flow by the Lorentz force ($-{\mbox{\boldmath $F$} {}}_L \cdot {\mbox{\boldmath $v$} {}}$) goes into the magnetic field (either locally, or, through the Poynting flux, elsewhere) and vice-versa. In a decretion flow azimuthal forces are aligned with the orbital motion, torquing up the disk. Therefore in Equation (\[E1\]) decretion flows imply ${\mbox{\boldmath $F$} {}}_L \cdot {\mbox{\boldmath $v$} {}}>0$. To balance Equation (\[E1\]) the power required to torque up the disk must come from either the local magnetic energy through the term $\partial_t B^2$, or from a deposit of magnetic energy generated elsewhere and transported by the Poynting flux (${\mbox{\boldmath $\nabla$} {}}\cdot {\mbox{\boldmath $S$} {}}$).

Calculating the Poynting flux using only the orbital velocity ${\mbox{\boldmath $v$} {}}= r \varOmega {{\bm{\hat{e}}_\phi}}$ we find $$\begin{aligned}
{\mbox{\boldmath $S$} {}}=& \left(-{\mbox{\boldmath $v$} {}}\times {\mbox{\boldmath $B$} {}}\right)\times {\mbox{\boldmath $B$} {}}/\mu_0   \label{SS_Orb} \\
&= -r \varOmega M_{\phi r} {{\bm{\hat{e}}_r}}- r \varOmega M_{\phi z} {{\bm{\hat{e}}_z}}+ r \varOmega \left(B_r^2 + B_z^2\right) {{\bm{\hat{e}}_\phi}}. \nonumber\end{aligned}$$ Taking the divergence of Equation (\[SS\_Orb\]) and averaging in azimuth we find $$\langle {\mbox{\boldmath $\nabla$} {}}\cdot {\mbox{\boldmath $S$} {}}\rangle_{\phi} = -\frac 1r \partial_r \left(r^2 \varOmega \langle M_{\phi r}\rangle_\phi\right) - \partial_z \left(r \varOmega
M_{\phi z}\rangle_\phi\right). \label{divS_Orb}$$ Equation (\[SS\_Orb\]) shows that horizontal stresses (such as $M_{\phi r}$) are associated with radial energy transport, while vertical stresses are associated with vertical energy transport. In disks with meridional circulation accreting through horizontal stresses, energy released in surface layers at a radius $r$ therefore does not travel vertically to power the decreting midplane at that same radius.

The horizontal component of Equation (\[divS\_Orb\]) can be approximated by assuming a power-law behavior $$M_{\phi r} \propto r^{-s}$$ to find $$\langle {\mbox{\boldmath $\nabla$} {}}\cdot {\mbox{\boldmath $S$} {}}\rangle_{\phi} =(1/2-s )  \varOmega  \langle M_{\phi r} \rangle_{\phi}.
\label{Poy1}$$ For comparison to this value, the local power released by accretion or required by decretion is approximately $$\frac 32 \varOmega  \langle M_{\phi r} \rangle_{\phi,t}.
 \label{eq:poy2}$$ Because $s$ is of order unity, [[Eq. (\[Poy1\])]{}]{} and [[Eq. (\[eq:poy2\])]{}]{} are of the same order, except in the special case of $s=1/2$. Thus, the Poynting flux extracts an order unity fraction of the accretion power from accreting regions, and injects similar amounts in decretion regions (again, except in the special case of $s=1/2$).

While we focus on the case of decretion to demonstrate the absolute importance of long distance energy transport, the order unity role played by the Poynting (or by similar arguments the kinetic energy) flux means that long distance energy transport also plays a crucial role in determining the energy density of the fluctuating fields that drive Maxwell (or Reynolds) stresses in disks accreting due to horizontal stresses powered by accretion energy. Accordingly, the saturated strength of those fields everywhere but the innermost edge of an accretion disk can be modulated by the rate at which energy is supplied by the inner disk, but only on the time scale required for that modulation to change the stresses. That time scale for the stresses to change sets a locality criterion: on longer time scales the accretion rate in the outer disk depends on the inner disk.

Inner Critical Radii and Disk Layering
--------------------------------------

![ Structure of a disk dominated by horizontal Maxwell stresses. In each layer, decretion occurs beyond the height-dependent cylindrical radius $r_c(z)$, where the radial dependence of the Maxwell stress drops below the critical value $\partial_r[r^2(- M_{\phi r}(z))] = 0$. Blue arrows in the bottom half show the direction of the accretion and decretion flow. For horizontal stresses other than the Maxwell stress, the structure is the same, with the flux and stress replaced as appropriate. Yellow arrows in the top half indicate the outward going Poynting flux, which grows though the accretion region, and deposits energy in the outer, decretion region, although the transition between the growing and shrinking fluxes occurs at $9r_b(z)/4$, which need not be the same as $r_c(z)$ and defines an alternate critical surface.[]{data-label="fig_fig1"}](fig1.pdf){width="\columnwidth"}

As we have reviewed above, in steady state viscous disk models, the vertically integrated radial mass flux is directed inwards while the midplane flow is directed outwards. This occurs because the much stronger radial gradients found at high altitudes drive fast inward flows that overwhelm the midplane’s slower outward flows. Horizontal stresses drive horizontal energy fluxes (Sect. \[sec\_E\_transport\]). Because the mass and energy fluxes are horizontal, this leads to a picture of accretion disks as made up of layered, nearly horizontal slabs rather than adjacent annuli (see Figure \[fig\_fig1\]). In this picture, each individual layer at a height $z$ above the midplane has its own critical radius $r_c(z)$ that defines the boundary between accretion and decretion. In a vertically isothermal disk in hydrostatic equilibrium, the density is $$\rho(r,z) = \rho_0(r) e^{-z^2/2H^2(r)},$$ where $\rho_0$ is the midplane density and $H$ the scale height. Taking the radial derivative we find $$\partial_r \ln \rho(r,z) = \partial_r \ln \rho_0 +\frac{z^2}{H^2(r)} \partial_r \ln H(r).
\label{flaring_eq}$$ The density, then, can only be taken to be a power-law in $r$ if the term capturing the radial dependence on $H$ is negligible compared to the radial dependence on the midplane density, which requires $z^2/H^2 \ll 1$. As noted in Section \[sec\_T\_propto\_p\], disk models with stresses strictly proportional to pressure predict decretion where the disk can behave in a scale-free, power law manner, i.e. where the radial dependence on H in Equation (\[flaring\_eq\]) can be neglected. This implies that for each height $z$, the critical radius separating accretion and decretion will be near the radius $r_{pl}(z)$ outside of which the disk, at the height $z$, can be approximated as a powerlaw in radius. While the precise location of $r_{pl}$ will depend on the details of the disk, Equation (\[flaring\_eq\]) gives the constraint $H(r_{pl}(z)) > z$.

Further, standard viscous disk theory [@2002apa..book.....F] for vertically integrated disks invokes another critical radius $r_b$ at which $\partial_r \varOmega=0$. Outside of $$r > \frac 94 r_b$$ viscous accretion flows locally dissipate more energy than provided by the accretion power. However, while vertical hydrostatic equilibrium can mathematically be extended to infinity, in practice accretion disks will generally be truncated vertically by stellar or disk winds, or the presence of infalling envelope material. It follows that each altitude $z$ has its own $r_b(z)$ where $\partial_r \varOmega(r,z)=0$. Accordingly, a similar layered picture results from a analysis of viscous disks performed layer-by-layer, as opposed to the more common vertically-integrated approach. We can construct the surface $$r(z)=\frac 94 r_b(z), \label{94r_star}$$ noting that it divides an outer region centered on the midplane that is a net consumer of energy from an inner, surface layer that is a net exporter of energy. (Because disks generally increase in height with radius, their upper surfaces are also their inner edges.) The energetics of global disks depend fundamentally on the material in the surface layer defined in Equation (\[94r\_star\]), because in viscous disk theory the accretion power below that surface at larger radii is inadequate to power the dissipation [@2002apa..book.....F].

Available Accretion Energy {#sec:power-decr-flows}
--------------------------

Can accretion at the inner edge of each layer of the disk provide enough energy to drive the outward motions expected further out? The local surface density of the energy released per unit time by accretion in a vertical slab of a Keplerian disk is $$Q_{acc}= \int_z  \frac 12 \rho v_r r \varOmega^2 dz, \label{acc_power_0}$$ where the factor of one half accounts for half of the gravitational potential energy being converted into orbital kinetic energy. Azimuthally integrating Equation (\[acc\_power\_0\]) for a given slab, we find the total power released at a radius $r$ to be $$\begin{aligned}
\nonumber \int_\phi Q_{acc} 2\pi r d \phi &= &
\frac 12 \left(2 \pi r \varSigma v_r r \varOmega^2\right)\\
& =& \frac 12 \dot{M} r \varOmega^2=\frac 12 \dot{M} r_f \varOmega_f^2 \left(\frac{r}{r_f}\right)^{-2}.
\label{acc_power_1}\end{aligned}$$ where $\dot{M}$ is the accretion rate (so that $\dot{M} < 0$ implies decretion), and $r_f$ is a fiducial reference position. In Equation (\[acc\_power\_1\]) $\Sigma$ is the surface density and $v_r$ is the density-weighted mean radial velocity of the slab.

Let us divide our slab into two annuli at $r_f$, which we choose to be where $v_r=0$, such that there is an inner accretion flow $\dot{M}_1>0$ with an inner edge at $r_0$, and an outer outwards flow $\dot{M}_2<0$ out to $r_1$, with $r_1 \gg r_f \gg r_0$. Then if the flows are to be powered by the accretion, net gravitational energy must be released by the system. The accretion power provided by the inner annulus is $$P_{in} = \int_{r_0}^{r_f} \frac {\dot{M}_1 r_f^3 \varOmega_f^2 r^{-2} }{2} dr= \frac 12 \dot{M}_1 \varOmega_f^2 \frac{r_f^2}{r_0} (r_f-r_0),$$ while the accretion power provided by the outer annulus is (negative because $\dot{M}_2<0$) $$P_{out} = \int_{r_f}^{r_1} \frac 12 \dot{M}_2 r_f^3 \varOmega_f^2 r^{-2} dr= \frac 12 \dot{M}_2 \varOmega_f^2 \frac{r_f^2}{r_1} (r_1-r_f).$$ The condition that net gravitational potential energy be released by this flow, allowing for accretion-powered accretion stresses, is then $$P_{in} + P_{out} >0,$$ or $$\left| \frac{\dot{M}_1}{\dot{M}_2} \right|> \frac{r_0}{r_1} \left(\frac{r_1-r_f}{r_f-r_0}\right) \simeq \frac {r_0}{r_f},$$ where the inequality holds for the broad annuli we assume, with $r_f/r_0 \gg 1$ and $r_1/r_f \gg 1$. Inner annuli accreting at a rate $\dot{M}_1$ therefore release adequate energy to drive broad regions of decretion with $|\dot{M}_2| \sim \dot{M}_1$, although this is only a necessary requirement, and the details of the energy transport also matter.

Time Scales of Energy Transfer {#subsec_time_scales}
------------------------------

The divergence of the Poynting (and analogous hydrodynamical energy) fluxes is comparable to the accretion power and can drive broad outflowing regions. Therefore, we need to consider the impact that changes in the Poynting flux, due, for example, to fluctuations or outbursts in the inner disk, will have on the stresses. We can quantify this by following the analysis of BP99. They derived the local steady state rate of energy dissipation into heat (their Eq. 37) starting from the full energy evolution equation: $$\frac{\partial \phantom{t}}{\partial t}\left(\frac 12 \rho v^2 + \rho \Phi+\frac{B^2}{2 \mu_0}\right) + {\mbox{\boldmath $\nabla$} {}}\cdot {{\mbox{ ${\cal F}$} {}}_E}= -{\mbox{ ${\cal Q}$} {}}_e
\label{EQ_QQQ}$$ (their Eq. 31). To find the steady state rate, they assume that the time derivatives are negligible, thus identifying the divergence of the energy flux with the local energy dissipation rate ${\mbox{ ${\cal Q}$} {}}_e$ (positive ${\mbox{ ${\cal Q}$} {}}_e$ means energy is deposited as heat). Note that the quantity we define as ${\mbox{ ${\cal Q}$} {}}_e$ is the volume density of the energy dissipation; so the variable $Q_e$ defined by BP99, the total energy dissipated per unit time per unit surface area, is $$Q_e = \int_z dz {\mbox{ ${\cal Q}$} {}}_e.$$ That is twice the total power that must be radiated from each disk surface element. Note that $Q_e$ is the surface power actually dissipated into heat and radiated away, while $Q_{acc}$ was the accretion power surface density.

In the case of accretion flows driven by Maxwell or turbulent Reynolds stresses, the accretion or decretion power is mediated by the energy density of the magnetic field or of the turbulent flows. If there is a mismatch between the local accretion power and the divergence of the energy fluxes, energy must be injected into or extracted from the local magnetic or turbulent energy densities (see Equation \[E1\] for the magnetic case). For a slab, those azimuthally averaged energy surface densities are $$\langle E_{T} \rangle_{\phi} = \int_z dz \left\langle \frac 12 \rho w^2 +\frac{B^2}{2\mu_0}\right\rangle_\phi,
\label{E_T_def}$$ where we recall that the velocity has been decomposed into a mean orbital $\bar{v}$ and a fluctuating $w$.

The terms that survive averaging are related to the stresses: we have $R_{ij} \equiv \rho w_i w_j$ and $M_{ij} \equiv \mu_0^{-1} B_i B_j$. Accordingly we can write $$E_T =  f \Sigma |W_{r \phi}|
\label{c_def}$$ where $f>1$ measures the difference between the squared amplitudes of the magnetic and fluctuating velocity fields $B^2$ and $w^2$ and their horizontal correlators $|B_\phi B_r|$ and $|w_\phi w_r|$. The parameter $f$ relates the stresses to an energy density, and so is a relative of the Shakura-Sunyaev $\alpha$. Indeed, [@2011ApJ...738...84H] defined $$\alpha_{mag} = -2 \frac{B_r B_\phi}{B^2},$$ finding $\alpha_{mag} \sim 0.3$. We can therefore estimate $f^{-1} \sim \alpha_{mag}/2 \sim 0.1$ for MRI active disks.

We can find the accretion power surface density in terms of the stresses from Equation (28) of BP99, $$P_{acc} = \Sigma R \varOmega_K^2 \langle u_r \rangle_\rho \sim 2 g \varOmega \Sigma W_{r\phi}$$ where we have assumed that $\Sigma R^2 W_{r \phi} \propto r^{-g}$. We see then that changes in the divergence of the energy fluxes cause changes in the local stresses and hence accretion rates on a time scale $$\left| \frac{E_T}{P_{acc}}\right| =\left| \frac{f}{2g\varOmega} \right|.$$ It immediately follows that local regions in accretion disks are only buffered against changes in the incoming energy fluxes for orbital to tens of orbital time scales. This means that the locality time scale for disk accretion is only a few dynamical times, far shorter than the global viscous evolution time. This is most clearly the case in decreting regions, where the energy that powers the decretion must come from somewhere, but even in accreting regions, putting energy into or removing energy from the magnetic field or the turbulence will alter the stresses. Note that this time evolution of $E_T$ was set to zero in BP99 even when they considered the non-steady state to derive their Equation 46, as we discuss in detail in Sect. \[sec:BP99\].

The Shearing Box Approximation {#sec:shear-box-appr}
==============================

The shearing box approximation neglects curvature terms and radial gradients of the hydrostatic disk background distribution of density and temperature. The curvature terms are consistently neglected by taking the horizontal size of the shearing box to be small compared to the radial position $r_0$. This flattening of the differential operators allows replacement of the cylindrical coordinate system with a local, Cartesian one. To arrive at a system where shear periodic boundary conditions can be applied in the radial direction, it is also necessary to neglect a disk’s background radial density and temperature gradients. For the hydrostatic equilibrium pressure $P_h$ that characterizes the global hydrostatic disk structure, this approximation about some radius $r_0$ is $$\begin{aligned}
P_h(r,z) = P_{h}(r_0,z) +  \mathcal{O}\left(\frac{r-r_0}{r_0}\right) \  \label{eq_phapprox}\end{aligned}$$ . This approximation has significant consequences a few scale heights $H$ above the midplane, no matter how cold and thin the disk is, because the approximation cannot capture a radially varying $H$. In particular, meridional flows cannot be driven once this approximation has been made.

The shearing box as commonly employed does, however, provide an internally consistent model for a shearing disk-like flow, if not a strictly valid asymptotic approximation to a section of an accretion disk in hydrostatic equilibrium. Given a shearing box, the underlying conservation properties of magnetohydrodynamics are preserved, up to specific effects of work done at the boundaries [@1995ApJ...440..742H; @2014arXiv1406.4864M]. As such, the shearing box has and can be used to understand the local dynamics in a shear flow closely analogous to that of a disk. However, the mass flux through a vertical surface drawn in a disk cannot be determined from this shearing box model alone, due to the neglect of the radial derivatives of the background density and temperature structure and the resulting elimination of any meridional flows.

The Poynting Flux in the Shearing Box {#subsec_Energy_fluxes}
-------------------------------------

Consider an incompressible ideal unstratified shearing box. where the velocity $\bm{v}$ has been decomposed into ${{\bm{\bar{v}}}}+ \bm{w}$ with the background shear flow $\bm{\bar{v}} = -q\varOmega x {{\bm{\hat{e}}_y}}$, $\bm{w}$ the fluctuation, and $q=3/2$ in the case of Keplerian rotation with $\varOmega \propto r^{-3/2}$. The governing equations are: $$\begin{aligned}
\partial_t \bm{w} &+ {{\bm{\bar{v}}}}\cdot \nabla \bm{w} + \bm{w}\cdot \nabla \bm{w} = \nonumber\\
& -2\varOmega {{\bm{\hat{e}}_z}}\times \bm{w} +q\varOmega w_x{{\bm{\hat{e}}_y}}+\frac{1}{\mu_0 \rho} (\nabla\times{{\bm B}}) \times{{\bm B}},  \\
\partial_t {{\bm B}}&= \nabla\times\left[ (\bm{w}+q\varOmega x{{\bm{\hat{e}}_y}}) \times {{\bm B}}\right],\end{aligned}$$ with $\nabla\cdot \bm{w}=0$.

An exact solution of these governing equations is $$\begin{aligned}
\bm{w}(t) &= 0, \\
{{\bm B}}(t) &=  \sin(k_z z) \left( {{\bm{\hat{e}}_x}}+ q\varOmega t {{\bm{\hat{e}}_y}}\right), \label{eq_path_Bt} \\
\mu_0 {{\bm J}}(t) &= \nabla \times {{\bm B}}(t) = k_z \cos (k_z z) \left(-q \varOmega t {{\bm{\hat{e}}_x}}+ {{\bm{\hat{e}}_y}}\right),\end{aligned}$$ which has zero net magnetic field at all times. The Lorentz force $${\mbox{\boldmath $F$} {}}_L = {{\bm J}}\times {{\bm B}}= \frac{k_z}{\mu_0}  \sin(k_z z) \cos(k_z z) \left(1 +q^2 \varOmega^2\right) {{\bm{\hat{e}}_z}}$$ merely captures the magnetic pressure. The fluid is incompressible, so no fluid is accelerated, the kinetic energy of the flow does not evolve and no work is done by the flow.

The only time varying energy in the system is magnetic energy, which evolves according to: $$\begin{aligned}
\frac{\partial}{\partial t} \left(\frac{B^2 }{2 \mu_0} \right)&= -\nabla \cdot \left[\frac{1}{\mu_0}\left(-{\mbox{\boldmath $v$} {}}\times {\mbox{\boldmath $B$} {}}+\eta{\mbox{\boldmath $J$} {}}\right)\times{\mbox{\boldmath $B$} {}}\right] \nonumber\\
&\qquad -\eta {\mbox{\boldmath $J$} {}}^2 - ({\mbox{\boldmath $J$} {}}\times {\mbox{\boldmath $B$} {}}) \cdot {\mbox{\boldmath $v$} {}}.
\label{eq_Bener_evolution}\end{aligned}$$ In the ideal limit ($\eta=0$), this has the solution $$\begin{aligned}
\frac{B^2(t)}{2 \mu_0} & = \frac{\sin^2(k_z z) \left(1+ q^2 \varOmega^2 t^2\right) }{2 \mu_0} \ ,\end{aligned}$$ which can also be obtained directly from Equation \[eq\_path\_Bt\].

The rate of energy density growth is then $$\begin{aligned}
\frac{\partial}{\partial t} \frac{B^2(t)}{2\mu_0} & = \frac{1}{\mu_0} \sin^2(k_z z) q^2 \varOmega^2 t \ . \label{eqpathsb_dtmagenergy}\end{aligned}$$ No work is done inside the volume, so the energy must be provided by some flux from outside the volume, namely the Poynting flux ${{\boldsymbol{S}}}$.

For this exact solution, $$\begin{aligned}
{\mbox{\boldmath $S$} {}}&= \frac{1}{\mu_0}  q\varOmega x  \sin^2 (k_z z)  \left[  -  q \varOmega t {{\bm{\hat{e}}_x}}+  {{\bm{\hat{e}}_y}}\right]\, .
\label{eqpathsb_poynting}\end{aligned}$$ The Poynting flux deposits energy into the field at the rate $$\begin{aligned}
-\nabla \cdot \bm{S} &= \frac{ q^2 \varOmega^2  t}{\mu_0}  \sin^2(k_z z) \label{eqpathsb_divs}\ ,\end{aligned}$$ which is the same rate that the magnetic energy grows (Equation \[eqpathsb\_dtmagenergy\]).

We reemphasize that the shear does no work in this scenario, as the associated Lorentz forces are directed in the ${{\bm{\hat{e}}_z}}$ direction, perpendicular to the fluid velocity in the ${{\bm{\hat{e}}_y}}$ direction, so the fluid does not encounter any resistance to its movement. Therefore, even though the energy growth rate contains a factor of the shear rate $q\varOmega$ (Equation \[eqpathsb\_dtmagenergy\]), the role of the shear is limited to setting the field geometry, which yields the Poynting flux (Equations \[eqpathsb\_poynting\] and \[eqpathsb\_divs\]).

Energy Fluxes in Shearing Boxes
-------------------------------

Starting with the ideal MHD energy flux and generalizing BP99 (Eq. 32), one arrives at $$\begin{aligned}
{{\mbox{ ${\cal F}$} {}}_E}= \bm{v} \left( \frac{1}{2} \rho v^2 + \rho \Phi + P \right) +\frac{{{\bm B}}}{\mu_0} \times (\bm{v} \times {{\bm B}}) - \bm{\sigma'}(\bm{v}) \cdot \bm{v}\end{aligned}$$ which includes the kinetic energy, potential energy, thermal energy, magnetic energy, and viscous fluxes. We use the same notation as in the previous section, but in addition, $P$ is the thermal pressure and $\Phi$ is the gravitational potential. The viscous stress tensor is $$\begin{aligned}
\sigma'_{ik}(v) = \nu\rho \left(\frac{\partial v_i}{\partial x_k} +\frac{\partial v_k}{\partial x_i} -\frac{2}{3}\delta_{ik}\frac{\partial v_l}{\partial x_l}\right) + \zeta \delta_{ik} \frac{\partial v_l}{\partial x_l},\end{aligned}$$ where $\nu$ is the viscosity, and $\zeta$ the bulk viscosity [@LandauAndLifshitz].

In a shearing box, the $y$ and $z$ components of ${{\mbox{ ${\cal F}$} {}}_E}$ are periodic, but the $x$ component is not: $$\begin{aligned}
{{\mbox{ ${\cal F}$} {}}_E}\cdot {{\bm{\hat{e}}_x}}=
   & -q\varOmega x \rho  w_x w_y +w_x \left( \frac{1}{2} \rho (\bar{v}^2 +\bm{w}^2  )  + \rho \Phi + P \right) \nonumber \\
   &  +q\varOmega x \frac{1}{\mu_0} B_x B_y + \frac{1}{\mu_0} w_x B^2 - \frac{1}{\mu_0} B_x w_k B_k \nonumber\\
   &  +q\varOmega x \nu\rho \left[ \frac{\partial (-q\varOmega x)}{\partial x} \right]   +q\varOmega x \sigma'_{xy}(\bm{w})  \nonumber\\
   & -  w_y  \nu\rho \left[ \frac{\partial (-q\varOmega x)}{\partial x} \right]   -  w_i \sigma'_{x i} (\bm{w})
   \label{F_e_dot_e_x}\end{aligned}$$ Hence, the divergence of the energy flux contains the terms $$\begin{aligned}
\nabla \cdot {{\mbox{ ${\cal F}$} {}}_E}=
   & -q\varOmega  \rho  w_x w_y +q\varOmega  \frac{1}{\mu_0} B_x B_y -\nu\rho q^2 \varOmega^2  + ...
   \label{stress_terms}\end{aligned}$$ This expression is essentially a differential version of equation (8) from @1995ApJ...440..742H. These are the energy sources that result from the non-periodic nature of the energy flux in the shearing box. Note that in vacuum, kinetic energy fluxes are zero and Poynting fluxes reduce to radiation fluxes. This means that the energy fluxes can have physically motivated boundary conditions when generalized from shearing boxes to global simulations with disks embedded in vacuum.

Viscosity-driven Energy Flux in the Shearing Box
------------------------------------------------

We can usefully compare shearing boxes behaving viscously due to Maxwell, Reynolds (or even microphysical viscous) stresses to shearing boxes with imposed constant viscosity $\nu$. A viscous shearing box flow with $\bm{w}=0$ also has a non-zero divergence of the energy flux [@1974MNRAS.168..603L], but this is accompanied by work done against the viscous friction force. To characterize this, it is useful to state the evolution equation for kinetic energy in an incompressible viscous hydrodynamic flow $$\begin{aligned}
\frac{\partial}{\partial t}\left( \frac{1}{2} \rho v^2\right) &= -\nabla \cdot \left[\rho \bm{v} \left(\frac{1}{2} v^2 +\frac{p}{\rho} \right) -\bm{v} \cdot \bm{\sigma'}\right] - \sigma'_{ik} \frac{\partial v_i}{\partial x_k} \ .
\label{eq_Kener_evolution}\end{aligned}$$ In the steady $\bm{w}=0$ shearing box flow without Reynolds stresses, this reduces to a balance between the kinetic energy source from the divergence of the kinetic energy flux, and the dissipation of kinetic energy to heat from the work done against viscous friction: $$\begin{aligned}
-\nabla \cdot \left[ -\bm{v} \cdot \bm{\sigma'}\right] &=\sigma'_{ik} \frac{\partial v_i}{\partial x_k} \ ,\\
-\nabla \cdot \left[ -\bm{{{\bm{\bar{v}}}}} \cdot \bm{\sigma'}(\bm{{{\bm{\bar{v}}}}}) \right]  &=  \sigma'_{yx}(\bm{{{\bm{\bar{v}}}}}) \frac{\partial {{\bm{\bar{v}}}}}{\partial x} \ , \\
 -\nu \rho q^2 \varOmega^2 &=  -\nu \rho q^2 \varOmega^2 \ .\end{aligned}$$

The dynamics of the energy flux that results from constant viscosity $\nu$ in a shearing box hence differ strongly from that which results from the shearing of magnetic fields. This is because the energy consumed by the shearing of magnetic fields can remain as magnetic energy, increasing the Maxwell stresses and driving further accretion. By contrast, in the viscous case the work done against viscous friction immediately thermalizes any energy left over from the radial energy flux, which includes the accretion or decretion power. The assumption of constant viscosity means that this energy plays no further role in the accretion process. This can be seen by contrasting the magnetic energy evolution (Eq. \[eq\_Bener\_evolution\]) to the kinetic energy evolution (Eq. \[eq\_Kener\_evolution\]), where in the kinetic energy case, the source of kinetic energy due to viscosity in the first term is exactly balanced by the final term, which removes this energy to heat at the same time.

Simulations of Accretion Disks {#sec:simul-accr-disks}
==============================

Applicability of Shearing Boxes {#sec_interpretation_of_sb}
-------------------------------

Shearing boxes may lead to unrealistic results when used to model local regions in accretion disks controlled by horizontal Maxwell stresses $M_{\phi r}$. This occurs because the Poynting flux into the simulation domain is unphysically set by the boundaries of the simulation domain, with energy being supplied at the boundaries of the periodic volume. In a given domain, the energy in the magnetic field, and hence the magnetic stress itself, depends on the energy supply. Equivalent constraints also apply to horizontal Reynolds stresses and their accompanying hydrodynamic energy fluxes. Our conclusions about $M_{\phi r}$ are, however, primarily based on overall energy considerations, and so do not address the extent to which local turbulent properties can still be addressed by local models.

Other recent work has indeed called that into question, though. One example is that the MRI appears to expand to the artificially constrained box scale [@2012MNRAS.422.2685S; @2012ApJ...748...79Y; @2014MNRAS.441.1855N], again arguing that the local treatment must be carefully interpreted. Another example is a comparison by of the growth of isolated MHD perturbations in the center of a periodic, 2D, shearing box with the growth of the same perturbations in an otherwise identical domain with free-slip, conducting walls on the radial boundary. When the boxes had width $L_x = \pi$, such that the radial boundaries were near the initial perturbations, the shearing box generated spurious energy in the box. However, when they moved the radial edges away by increasing the width of the box to $L_x = 2 \pi$, they found much better agreement between the shearing box and the wall-bounded flow. Their interpretation of this phenomenon was that the shearing box walls were adding energy to the perturbations. We interpret this as being due to the implied Poynting flux (Sect. \[subsec\_Energy\_fluxes\]).

Nonetheless, there is evidence that the small scale behavior of the MRI can be modeled in restricted domains. @2012ApJ...749..189S examined the behavior of MRI in restricted subdomains of a global model, and found that the correlations of the turbulent fluctuations remained similar. Nonmodal analysis of the most energetic structures over finite growth times in a laboratory-like global MRI configuration also suggests that the shearing waves are a dominant feature in that global problem, and are captured well in the local shearing box approximation [@2014arXiv1407.4742S].

 Existing Global Simulations {#sec:glob-simul-exist}
----------------------------

In general, global simulations could show either accretion or decretion locally at each radius. Further, the boundary conditions can force the flow. and find both accretion and decretion in global models, with accretion at small radii and decretion at larger radii. Additionally, searched for and did not detect a meridional circulation in an MRI active disk. They found in this case that the variation of stress is not proportional to the pressure in the vertical direction as would occur in a solution with vertically constant alpha, explaining their lack of detection of the expected circulation.

@2014ApJ...784..121S performed global simulations with initial net vertical fields. These simulations are especially interesting because they combine net vertical field with a radial temperature gradient that forces vertical shear. Their analysis of radial flows [@2014ApJ...784..121S Section 5.4] is consistent with our analysis. Indeed they estimate that accretion occurs in the interior section of their disk, and decretion occurs in the exterior section as we predict.

In these models, where some mass flows outwards and to higher total energy, the flow must be supplied with energy liberated from some other location where accretion occurs, or from the boundaries. Indeed, in each case, such locations exist, consistent with our findings.

Comparing Local and Global Simulations
--------------------------------------

[@2010ApJ...712.1241S] compared the flux-stress relation from a global simulation with the results for shearing boxes from [@2007ApJ...668L..51P]. They found a qualitative agreement, but also a marked discrepancy, which they suggested was due to inadequate resolution in the global simulations. This resolution difficulty was also noted by [@2011ApJ...738...84H] who attempted to derive resolution criteria from shearing boxes for use in global simulations.

From a different angle, [@2012ApJ...749..189S] compared small sub-subdomains of a small subdomain of a global simulation with each other and with the original subdomain (but did not invoke the shearing sheet approximation). They found that there was general agreement between the sub-volumes, confirming that local simulations with physical boundary conditions can be used to study small scale phenomena below the scale of their subdomain. Overall, our argument is consistent with published numerical simulations.

Diagnostics {#subsec-cons}
-----------

We have noted that shearing boxes are not subject to energy constraints present in accretion disks; and that global simulations of accretion disks may need to contend with locality time scales of a few tens of orbits. This raises the question of how meaningful our arguments are for global simulations: while the locality time scale we find is just a few orbits, this only allows, but does not force, the disk to behave on those time scales. We propose two diagnostics to measure the effect in practice.

Firstly, the outer-disk growth rate of magnetic instabilities such as the MRI will be partially slaved to the inner, faster growing magnetic field. Radial fluxes originating in the rapidly evolving inner disk should increase the outer disk’s MRI growth rate beyond that following from a radially local analysis. There are suggestions of this behavior in the literature . Secondly, the MRI has strongly time varying stresses and we expect the magnetic stresses can be correlated across large radial extents with a modest lag (on order of the orbital timescale). There are tantalizing hints of such behavior: figure 5 of shows the azimuthal magnetic field at $6$ AU being correlated with that at $4$ AU with a lag of only about $100$ years, or about $10$ local orbits. However, the correlations of the stresses across radius and time are far less clear (their figure 9). Further numerical simulations with more focused diagnostics could better quantify these stress correlations.

Energy Flux Boundary Conditions {#sec_choice_BC}
-------------------------------

If the Poynting flux at the inner edge of the simulation domain is not zero, then it is either being imposed by the boundary conditions or set by the simulation volume itself. Unless the boundary condition for the Poynting flux is somehow physically determined (e.g. by a stellar magnetosphere), the energy made available to the system will not be physically controlled, and may be significantly larger or smaller than the value in a simulation that included the entire disk. If, on the other hand, the Poynting flux is set to zero at a physical location within the accretion disk, then the energy available to the simulation volume will be significantly below the correct value because the accretion energy liberated inside of this location is not being transported any further into the domain. This is the case for the common choice of absorbing boundary conditions, which resistively destroy the magnetic field in buffer zones at the inner and outer edges of the simulation domain, forcing the Poynting flux to zero.

In vacuum, the Poynting flux reduces to the radiative flux, justifying radial boundaries positioned outside the disk with zero Poynting flux boundary conditions. At the outer edge this is reasonable: the evolution time scale is long, and the outer edge of the simulation domain can be put far from the disk, so that the outer edge of the disk will not cross the boundary during the lifetime of the simulation. At the inner edge however, the edge of the computational domain cannot be put arbitrarily far from the disk edge, and the disk evolution time is fast, so the inner edge of the disk will rapidly overrun the boundary. Once that occurs, the appropriate outwards energy flux entering the simulation domain from accreting gas cannot be tracked.

We need therefore to capture the structure of the inner edge of the disk directly. This is possible in simulations that include the surface of the central object, such as black hole accretion simulations with horizon-penetrating coordinates and computational domains that extend inside the horizon [ex: @2004ApJ...611..977M]. For circumstellar systems where the accretion disc extends smoothly to the stellar surface, this surface provides the appropriate location for the inner radial boundary [@2002MNRAS.330..895A; @2002ApJ...571..413S]. For circumstellar disks around magnetized stars, though, we suggest instead placing the inner boundary at the interface between the disk and the stellar magnetosphere: energy released by matter accreting in the magnetosphere will travel along the stellar magnetic field lines rather than entering the disk. In practice, this means using a conducting porous boundary that allows both mass flow and magnetic flux detachment. Note that this latter requires resolving the boundary layer that will form as the disk magnetic field piles up against the stellar field, as is done for example by in two-dimensional models, and @2012MNRAS.421...63R in three-dimensional models. These models indeed show strong time-variability, which our results suggest may be required to maintain accretion.

An alternative radial magnetic boundary condition may be to bracket the simulation volume with dead zones dominated by non-ideal effects. If the non-ideal effects are adequate to decouple an inner active region of an accretion disk from the outer active region, then the Poynting flux will not be dynamically important in the outer region. However, in stratified disks with active surface layers, this is unlikely to be the case, and recent results have shown that the Hall effect may allow large scale magnetic fields to develop even in dead zones . Nonetheless, conveniently, simulations using dead-zones as magnetic boundaries can check their self-consistency post-hoc by evaluating the radial Poynting flux in those self-same dead-zones, determining whether or not they are dynamically significant.

Extending BP99 {#sec:BP99}
==============

In their classic work, BP99 wrote the evolution of the non-thermal energy density in a viscous disk in their Equation (31), reproduced below: $$\begin{aligned}
\frac{\partial}{\partial t} &\left( \frac 12 \rho v^2 +  \rho \Phi +\frac{B^2}{8\pi}\right)  \nonumber \\
&+ {\mbox{\boldmath $\nabla$} {}}\cdot \left[{\mbox{\boldmath $v$} {}}\left(\frac{1}{2}\rho v^2 + \rho \Phi +P\right) +\frac{{\mbox{\boldmath $B$} {}}}{4\pi}\times({\mbox{\boldmath $v$} {}}\times{\mbox{\boldmath $B$} {}}) \right] \nonumber \\
&=P {\mbox{\boldmath $\nabla$} {}}\cdot {\mbox{\boldmath $v$} {}}- \eta_V \left(\delta_i v_j\right)\left(\delta_i v_j \right) -\frac{\eta_B}{4\pi} |{\mbox{\boldmath $\nabla$} {}}\times {\mbox{\boldmath $B$} {}}|^2\, , \end{aligned}$$ where ${\mbox{\boldmath $v$} {}}= R \varOmega {{\bm{\hat{e}}_\phi}}+{\mbox{\boldmath $u$} {}}$ is the complete velocity field, and $\varOmega$ is purely determined by the gravitational potential $\Phi$ (so ${\mbox{\boldmath $u$} {}}$ includes terms deriving from radial pressure support and does not have a zero average). The square brackets contain the energy flux and Gaussian CGS units are used. In this section we follow the notation of BP99 and use $R$ for the cylindrical radius.

When considering non-steady state cases (and after vertically integrating and averaging) they arrived at their Equation (46) for the energy dissipated into heat per unit surface area per unit time, also reproduced below: $$-Q_e =  \Sigma W_{R\phi} \frac{d\varOmega}{d\ln R}\, ; \label{Eq_46_base}$$ see also our Section \[subsec\_time\_scales\] and Equation (\[EQ\_QQQ\]). In the analysis adopted by BP99, Equation (\[Eq\_46\_base\]) is accurate to leading order. It is instructive to consider the full expression, with terms of all orders retained. In deriving this expression we use the same procedure as in that work, with two exceptions. Firstly, for a given quantity $X$ we average using $$\begin{aligned}
\langle X\rangle_\rho &\equiv
   \frac{1}{2\pi R \Sigma}  \int_{-\infty}^{+\infty} \int_0^{2\pi}  \rho X R d\phi dz \, , \label{ouravg}\\
 \Sigma &\equiv   \frac{1}{2\pi R}  \int_{-\infty}^{+\infty} \int_0^{2\pi} \rho R d\phi dz\, .\end{aligned}$$ Secondly we define a cylindrical part -spherical part decomposition of $\Phi$ $$\begin{aligned}
\Phi = \Phi_C(R) + \Phi_S(R,Z)\, , \quad \Phi_S(R,0) = 0\, .\end{aligned}$$

After significant algebra using the simplifications laid out by BP99, but no approximations, we arrive at a version of Equation (\[Eq\_46\_base\]) with all terms retained: $$\begin{aligned}
-Q_e =&  \Sigma   W_{R\phi}   \frac{\partial \varOmega }{\partial \ln R} + {\mbox{ ${\cal A}$} {}}+ {\mbox{ ${\cal B}$} {}}+ {\mbox{ ${\cal C}$} {}}+ {\mbox{ ${\cal D}$} {}}+ {\mbox{ ${\cal E}$} {}}\label{Eq_46_full}\end{aligned}$$ $$\begin{aligned}
 {\mbox{ ${\cal A}$} {}}=& \frac{\partial }{\partial t}  \Sigma \left\langle \frac{1}{2}u^2 + \frac{B^2}{8\pi\rho} \right\rangle_\rho  \\
                       =&   \Sigma \frac{\partial }{\partial t}   \left\langle \frac{1}{2}u^2 + \frac{B^2}{8\pi\rho} \right\rangle_\rho +
                           \left\langle \frac{1}{2}u^2 + \frac{B^2}{8\pi\rho} \right\rangle_\rho \frac{\partial }{\partial t}  \Sigma \nonumber \\
 {\mbox{ ${\cal B}$} {}}=& \frac{\partial }{\partial t}  \Sigma R\varOmega \langle  u_\phi  \rangle_\rho \\
 {\mbox{ ${\cal C}$} {}}=&- \frac{1}{R}\frac{\partial }{\partial R} R
       \Sigma \left\langle \frac{1}{2} u^2 u_R 
      + \frac{ B_z^2 }{4 \pi \rho} u_R  + \frac{  B_R B_z }{4 \pi \rho} u_z  \right. \nonumber \\
  &  \qquad \qquad \qquad \qquad \left. +\frac{B_\phi^2}{4 \pi \rho} u_R  
                            + \frac{B_R B_\phi}{4 \pi\rho} u_\phi   \right\rangle_\rho \\
  {\mbox{ ${\cal D}$} {}}=& - \frac{1}{R}\frac{\partial }{\partial R} R \Sigma \left\langle \frac{ P}{\rho} u_R \right\rangle_\rho  \\
  {\mbox{ ${\cal E}$} {}}=&   \frac{\partial }{\partial t} \Sigma \langle \Phi_S \rangle_\rho
    - \frac{1}{R}\frac{\partial }{\partial R} R \Sigma \langle \Phi_S u_R\rangle_\rho\, . \end{aligned}$$ This expression is exact for the averaging as specified by Equation (\[ouravg\]), however it can also be taken as an approximation when expressed in terms of the mean flow variables by extending the averaging over an appropriate small radial interval $\Delta R$; these details are discussed in Appendix \[appendix\_averaging\].

It is difficult to determine the order at which small perturbed quantities enter the terms on the right-hand side of Equation (\[Eq\_46\_full\]). The perturbed velocity $u$ is no longer asymptotically small once turbulence has saturated; we do not *a priori* know the value of the ratio of perturbed velocity components $\langle u_{\phi} u_r\rangle_\rho/\langle u^2\rangle_\rho$ but it could easily be less than a tenth; and, crucially, we do not know the time and length scales on which values vary.

However, we merely need to argue that effects beyond those captured in Equation (\[Eq\_46\_base\]) can arise. To do that, we continue by adding only the first term ${\mbox{ ${\cal A}$} {}}$ to Equation (\[Eq\_46\_base\]). We neglect the second term of ${\mbox{ ${\cal A}$} {}}$, proportional to $\partial_t
\Sigma$, and ${\mbox{ ${\cal B}$} {}}$, because we expect the surface density, pressure and hence also $\langle u_\phi \rangle_\rho$ to vary only slowly in time, on the long, viscous time scale. We can neglect ${\mbox{ ${\cal C}$} {}}$ and ${\mbox{ ${\cal D}$} {}}$ in favor of ${\mbox{ ${\cal A}$} {}}$ if background and fluctuating quantities vary more slowly in space than in time, i.e. if the fluctuating terms vary over length and time scales $\ell$ and $\tau$ such that $$\frac{\ell}{\tau} \gg u.$$ Finally, neglecting ${\mbox{ ${\cal E}$} {}}$ implies that the disk is thin enough that the spherical component of the potential $\Phi_S$ can be neglected.

Under those conditions, we arrive at an approximate form of Equation (\[Eq\_46\_full\]): $$-Q_e \approx \Sigma   W_{R\phi}   \frac{\partial \varOmega }{\partial \ln R}
   +\Sigma \frac{\partial }{\partial t}   \left\langle \frac{1}{2}u^2
   + \frac{B^2}{8\pi\rho} \right\rangle_\rho.
\label{Eq_46_corr}$$ We further recall from Equation (\[E\_T\_def\]) that the energy in those reservoirs is related to the stress $W_{R\phi}$. Hence, the energy density in the rightmost term of Equation (\[Eq\_46\_corr\]) can be estimated as $f \Sigma W_{R\phi}$ by invoking Equation (\[c\_def\]). Neglecting the rightmost term of Equation (\[Eq\_46\_corr\]) therefore amounts to neglecting the time derivative of the stresses. The term is however vital if considering physical processes such as the linear growth phase of MRI, where the stresses grow on a timescale $\sim \varOmega^{-1}$.

It is further important because in Section \[subsec\_time\_scales\] we argued that the energy contained in $E_T$ can only sustain the power $Q_e$ for a few orbits. That means that a significant change in the incoming energy fluxes must change $E_T$ on a time scale of only a few times $\varOmega^{-1}$. As a result, the two terms on the right-hand side of Equation (\[Eq\_46\_corr\]) are of the same order, and the rightmost term must be kept if time variation in the stresses is contemplated. Future numerical experiments correlating stresses at different radii in global models will be needed to move beyond these order-of magnitude estimates and determine the actual speed of propagation of changes in the energy fluxes, and hence stresses. To diagnose the energetics, the various terms of Equation (\[Eq\_46\_full\]) can be individually tracked in numerical simulations.

Conclusions {#sec_conclusions}
===========

It has been recognized since the work of [@1974MNRAS.168..603L] that large swathes of disks that accrete through horizontal stresses, including by mechanisms such as the MRI, do not release adequate gravitational potential energy to power their dynamics. This can be seen most clearly in the context of meridional circulation and viscously spreading disks: any parcel of gas flowing away from the central object cannot be releasing gravitational potential energy. Indeed it is only in a narrow inner annulus that the accretion power is greater than the dissipation of energy into heat. In regions where the accretion power is inadequate, the energy evolution equation must be balanced by long distance energy fluxes, such as Poynting flux traveling radially through an MRI active disk, or stellar irradiation powering a thermal instability. Further, the mismatch between the local accretion power and energy dissipation is within an order of magnitude of the total local energy density in the fields generating the stresses, in the case of turbulent or magnetic stresses.

We have described two important consequences for accretion disks whose accretion stresses are powered by the accretion energy itself, including MRI active disks. Firstly, because the outer disk dissipates more energy than it provides, and imports the remainder from the inner disk, the outer disk’s behavior depends on the inner disk’s accretion rate. It follows that the energy flux, produced on the fast, dynamical time scale of the inner disk, must be capable of modulating the stresses, and hence accretion rate even in the outer disk with its slower dynamical time scale. In the case of global disks, this means that dynamics in the outer disk (which in viscous disk models needs to import energy even though it is accreting, when vertically integrated) should depend on the energy provided by the inner disk.

In more practical terms, for numerical simulations this means that energy fluxes through the radial boundaries can be important to first order in setting stresses and accretion rates. This behavior cannot be represented in shearing box simulations, where the box sets its own boundary conditions and hence energy fluxes. Shearing boxes are therefore decoupled from a constraint on the available energy that drives the magnetic and turbulent velocity fields exerting stresses on the disk. In the case where the energy fluxes play major roles, temporal changes in the fluxes should cause any stress-exerting turbulence present to temporarily enter forced-turbulence or decaying-turbulence modes. Even in steady states, the local power deficit means that the forced-turbulence picture may apply.

Secondly, the energy fluxes are large compared with the local energy densities associated with the accretion stresses in the sense that the local energy density can only power horizontal accretion stresses for a few orbits. If the time scale were much longer, for example the global disk evolution time scale, then the importance of the energy fluxes would be mitigated: the time scale for the outer disk to adjust itself to the energy fluxes from the inner disk would be too long to matter. However, because the time scale is short, the outer disk must rapidly adjust to changes in the inner disk, which in turn has short dynamical time scales. This short timescale may help explain observations of time-variable accretion flows in protoplanetary disks such as FUors and EXors [@2014arXiv1401.3368A].

In sum, we have demonstrated that understanding and accurately treating long distance energy fluxes in disks is vital in modeling and understanding accretion flows. Changes in the energy fluxes in the outer disk due to changes in the accretion rate in the inner disk impact the amplitude of the outer disk’s stresses on dynamical, rather than viscous, time scales. Additionally, we have suggested how to diagnose these phenomena in numerical simulations. In particular, we suggest both correlating stresses at different radii with time lags to determine the speed with which and the extent to which the stresses in the outer disk depend on the stresses (and hence energy release and energy fluxes) in the inner disk; and tracking the full set of contributions to the local energy dissipation (the terms of Equation \[Eq\_46\_full\]) and comparing their measured magnitudes to empirically determine which are dominant.

We thank P. Armitage, E. Blackman, O. Gressel, T. Heinemann, G. Lesur, S. Balbus and many other colleagues, including participants in the Facebook group Circumstellar Disks and Planet Formation, for comments that improved the manuscript. The research leading to these results has received funding from the People Programme (Marie Curie Actions) of the European Union’s Seventh Framework Programme (FP7/2007-2013) under REA grant agreement 327995 (CPM), and the U.S. National Science Foundation under CDI grant AST08-35734 and AAG grant AST10-09802 (WL, AH, JSO, and M-MML). AH received additional support from NASA OSS grant NNX14AJ56G. WL carried out his part of the work at the Jet Propulsion Laboratory, under contract with the California Institute of Technology, with support from the National Aeronautics and Space Administration (NASA) through the Sagan Fellowship Program executed by the NASA Exoplanet Science Institute.

[60]{} natexlab\#1[\#1]{}

, S. M., [Wilner]{}, D. J., [Hughes]{}, A. M., [Qi]{}, C., & [Dullemond]{}, C. P. 2010, [, 723, 1241](http://dx.doi.org/10.1088/0004-637X/723/2/1241)

, P. J. 2002, [, 330, 895](http://dx.doi.org/10.1046/j.1365-8711.2002.05152.x)

, M., [[Á]{}brah[á]{}m]{}, P., [Dunham]{}, M. M., [et al.]{} 2014, ArXiv e-prints, [[arXiv:1401.3368 \[astro-ph.SR\]]{}](http://arxiv.org/abs/1401.3368)

, X.-N. 2014, [, 791, 137](http://dx.doi.org/10.1088/0004-637X/791/2/137)

, S. A. 2003, [, 41, 555](http://dx.doi.org/10.1146/annurev.astro.41.081401.155207)

, S. A., & [Hawley]{}, J. F. 1991, [, 376, 214](http://dx.doi.org/10.1086/170270)

, S. A., & [Hawley]{}, J. F. 1998, [Reviews of Modern Physics, 70, 1](http://dx.doi.org/10.1103/RevModPhys.70.1)

, S. A., & [Papaloizou]{}, J. C. B. 1999, [, 521, 650](http://dx.doi.org/10.1086/307594)

, E. G., & [Pessah]{}, M. E. 2009, [, 704, L113](http://dx.doi.org/10.1088/0004-637X/704/2/L113)

, G., [Mignone]{}, A., [Cattaneo]{}, F., [Rossi]{}, P., & [Ferrari]{}, A. 2008, [, 487, 1](http://dx.doi.org/10.1051/0004-6361:200809730)

, F. 2007, [Science, 318, 613](http://dx.doi.org/10.1126/science.1147273)

, F. 2009, [Icarus, 200, 655](http://dx.doi.org/10.1016/j.icarus.2008.12.009)

, S. W., [Stone]{}, J. M., & [Pessah]{}, M. E. 2010, [, 713, 52](http://dx.doi.org/10.1088/0004-637X/713/1/52)

, C. P., [Hollenbach]{}, D., [Kamp]{}, I., & [D’Alessio]{}, P. 2007, Protostars and Planets V, 555

, N., [Flock]{}, M., [Turner]{}, N. J., [Klahr]{}, H., & [Henning]{}, T. 2010, [, 515, A70](http://dx.doi.org/10.1051/0004-6361/200912834)

, T., & [Stone]{}, J. M. 2003, [, 585, 908](http://dx.doi.org/10.1086/345848)

, M., [Dzyurkevich]{}, N., [Klahr]{}, H., & [Mignone]{}, A. 2010, [, 516, A26](http://dx.doi.org/10.1051/0004-6361/200912443)

, M., [Dzyurkevich]{}, N., [Klahr]{}, H., [Turner]{}, N. J., & [Henning]{}, T. 2011, [, 735, 122](http://dx.doi.org/10.1088/0004-637X/735/2/122)

, J., [King]{}, A., & [Raine]{}, D. J. 2002, [Accretion Power in Astrophysics: Third Edition]{} (Cambridge)

, S., [Lyra]{}, W., & [Masset]{}, F. 2011, [, 534, A107](http://dx.doi.org/10.1051/0004-6361/201016068)

, S., & [Nelson]{}, R. P. 2006, [, 457, 343](http://dx.doi.org/10.1051/0004-6361:20065643)

, P., & [Lynden-Bell]{}, D. 1965, , 130, 125

, O. 2010, [, 405, 41](http://dx.doi.org/10.1111/j.1365-2966.2010.16440.x)

, J. F., [Gammie]{}, C. F., & [Balbus]{}, S. A. 1995, [, 440, 742](http://dx.doi.org/10.1086/175311)

, J. F., [Guan]{}, X., & [Krolik]{}, J. H. 2011, [, 738, 84](http://dx.doi.org/10.1088/0004-637X/738/1/84)

, A. L. H., & [Armitage]{}, P. J. 2010, [, 719, 1633](http://dx.doi.org/10.1088/0004-637X/719/2/1633)

, C., & [Gail]{}, H.-P. 2004, [, 415, 1177](http://dx.doi.org/10.1051/0004-6361:20034629)

, H., & [Hubbard]{}, A. 2014, [, 788, 21](http://dx.doi.org/10.1088/0004-637X/788/1/21)

, W., & [Lin]{}, D. N. C. 1992, [, 397, 600](http://dx.doi.org/10.1086/171818)

, A., & [Pudritz]{}, R. E. 2000, Protostars and Planets IV, 759

, L. D., & [Lifshitz]{}, E. M. 1959, [Fluid mechanics]{} (Pergamon Press)

, G., [Kunz]{}, M. W., & [Fromang]{}, S. 2014, [, 566, A56](http://dx.doi.org/10.1051/0004-6361/201423660)

, D., & [Pringle]{}, J. E. 1974, , 168, 603

, J. C., & [Gammie]{}, C. F. 2004, [, 611, 977](http://dx.doi.org/10.1086/422244)

, C. P., [Hubbard]{}, A., [Yang]{}, C.-C., & [Mac Low]{}, M.-M. 2014, [, 791, 62](http://dx.doi.org/10.1088/0004-637X/791/1/62)

, C. P., & [Pessah]{}, M. E. 2014, ArXiv e-prints, [[arXiv:1406.4864 \[astro-ph.EP\]]{}](http://arxiv.org/abs/1406.4864)

, K. A., & [Stone]{}, J. M. 2000, [, 534, 398](http://dx.doi.org/10.1086/308736)

, F., & [Blackman]{}, E. G. 2014, [, 441, 1855](http://dx.doi.org/10.1093/mnras/stu706)

, E. R. 2014, [, 438, 2513](http://dx.doi.org/10.1093/mnras/stt2379)

, M. E., [Chan]{}, C.-k., & [Psaltis]{}, D. 2007, [, 668, L51](http://dx.doi.org/10.1086/522585)

, O., & [Gitelman]{}, L. 2002, [, 396, 623](http://dx.doi.org/10.1051/0004-6361:20021492)

, O., & [Umurhan]{}, O. M. 2008, [, 481, 21](http://dx.doi.org/10.1051/0004-6361:20078413)

, M. M., [Ustyugova]{}, G. V., [Koldoba]{}, A. V., & [Lovelace]{}, R. V. E. 2012, [, 421, 63](http://dx.doi.org/10.1111/j.1365-2966.2011.20055.x)

, M., [Bodenheimer]{}, P., & [Bell]{}, K. R. 1994, [, 423, 736](http://dx.doi.org/10.1086/173853)

, N. I., & [Sunyaev]{}, R. A. 1973, , 24, 337

, J., [Krolik]{}, J. H., & [Hirose]{}, S. 2010, [, 708, 1716](http://dx.doi.org/10.1088/0004-637X/708/2/1716)

, A. 1988, Acta Astronomica, 38, 21

, J. B., [Beckwith]{}, K., & [Armitage]{}, P. J. 2012, [, 422, 2685](http://dx.doi.org/10.1111/j.1365-2966.2012.20835.x)

, E., & [Goodman]{}, J. 2003, [, 341, 501](http://dx.doi.org/10.1046/j.1365-8711.2003.06431.x)

, K. A., [Reynolds]{}, C. S., & [Armitage]{}, P. J. 2010, [, 712, 1241](http://dx.doi.org/10.1088/0004-637X/712/2/1241)

, K. A., [Reynolds]{}, C. S., [Stone]{}, J. M., & [Beckwith]{}, K. 2012, [, 749, 189](http://dx.doi.org/10.1088/0004-637X/749/2/189)

, J., & [Bhattacharjee]{}, A. 2014, ArXiv e-prints, [[arXiv:1407.4742 \[astro-ph.HE\]]{}](http://arxiv.org/abs/1407.4742)

, A., & [Papaloizou]{}, J. C. B. 2002, [, 571, 413](http://dx.doi.org/10.1086/339892)

, T. K., & [Inutsuka]{}, S.-i. 2014, [, 784, 121](http://dx.doi.org/10.1088/0004-637X/784/2/121)

, T., & [Lin]{}, D. N. C. 2002, [, 581, 1344](http://dx.doi.org/10.1086/344437)

, W. M., & [Gail]{}, H.-P. 2007, [, 463, 369](http://dx.doi.org/10.1051/0004-6361:20065794)

, O. M., & [Regev]{}, O. 2004, [, 427, 855](http://dx.doi.org/10.1051/0004-6361:20040573)

, V. A. 1984, , 28, 50

, C.-C., [Mac Low]{}, M.-M., & [Menou]{}, K. 2012, [, 748, 79](http://dx.doi.org/10.1088/0004-637X/748/2/79)

, C., & [Ferreira]{}, J. 2013, [, 550, A99](http://dx.doi.org/10.1051/0004-6361/201220168)

Alternative Averaging {#appendix_averaging}
=====================

To read Equation (\[Eq\_46\_full\]) as an approximation expressed in terms of the mean flow variables one exchanges the averaging of Equation (\[ouravg\]) for BP99’s choice of averaging which includes averaging over an appropriate, small, radial interval $\Delta R$ given by: $$\begin{aligned}
\langle X\rangle_\rho &\equiv
   \frac{1}{2\pi \Sigma \Delta R}\int_{-\infty}^{+\infty} \int_{R-\Delta R/2}^{R+\Delta R/2}\int_0^{2\pi} \rho X d\phi dR' dz\ , \label{BP99avg}\\
 \Sigma &\equiv   \frac{1}{2\pi \Delta R}\int_{-\infty}^{+\infty} \int_{R-\Delta R/2}^{R+\Delta R/2}\int_0^{2\pi} \rho d\phi dR' dz\, . \label{BP99sigma}\end{aligned}$$ The averaging schemes are very similar, and the resulting expression Equation (\[Eq\_46\_full\]) we arrive at can be read using this averaging as approximation. When replacing the averaging in Equation (\[Eq\_46\_full\]) with the form Equation (\[BP99avg\]) the approximations implicitly invoked are of the form $$\begin{aligned}
\int_{R-\Delta R/2}^{R+\Delta R/2} R' X dR' \approx R \int_{R-\Delta R/2}^{R+\Delta R/2} X dR' \, . \label{BP99DeltaRapprox}\end{aligned}$$ Thus, the approximation is valid for $\Delta R/R \ll 1$.

[^1]: See the footnote on page 611 of that article: “For a Newtonian point mass \[...\] the power liberated at radii between $R$ and $R+dR$ is three times larger than the energy generated there.”

[^2]: A common alternative definition for the stress is $\mu_0M_{ij} \equiv B_iB_j -
      \sfrac{1}{2}\,\delta_{ij} B^2$, but this only matters for the diagonal (pressure) components.
