2.1. How standard disc theory is used

The ultimate goal of accretion theory is to explain and predict observations of accretes, most fundamentally, the observed flux, radiation spectra and temporal evolution. Accretion engines include not only emission from the disc but contributions from large-scale outflows, extended coronae or reprocessed emission from the central star, as in the case of YSO (Calvet et al. Reference Calvet, D’Alessio, Watson, Franco-Hernández, Furlan, Green, Sutter, Forrest, Hartmann, Uchida, Keller, Sargent, Najita, Herter, Barry and Hall2005). The ionization complexity, the mixture of material in different states of matter and the potential role of self-gravity make YSO accretion discs possibly the most challenging. A separate set of complexities for black hole accreting discs includes the effects of radiation pressure, general relativity for the most rapid rotators and the associated fact that outflows can be powered by magnetically mediated extraction of the rotational energy of the back hole. The latter is a separate energy reservoir from that of the accretion, although the latter would supply the magnetic field.

Figure 1.

Cyg X-1 flux spectrum as an example of how engine models are patched together to match observations. Typically the optically thick continuum contribution to the disc spectrum is computed from an SS73 type model while the optically thin coronal spectrum is separately computed from a hot optically thin inverse-Compton dominated plasma. The two are then added in the relative combination needed to match the observations. The reflected spectrum is self-consistently computed given the coronal spectrum and assumptions about the disc geometry. What is missing is a self-consistent theory that predicts, from first principles, the relative contribution of disc, corona and jet without having to patch together separate spectra empirically.

But even for an accretor without an additional energy source from a spinning central gravitater, the accretion itself can be a source of an optically thick disc, optically thin corona and outflows. Presently, modellers commonly construct spectral models by pasting together empirically weighted contributions from these ingredients at the fractions needed to match observations. Figure 1 shows this for Cygnus X-1, a well-known X-ray microquasar (Gierliński & Zdziarski Reference Gierliński, Zdziarski, Poutanen and Svensson1999). The three contributions needed to match the spectrum are the continuum optically thick disc emission, direct optically thin coronal emission and disc-reprocessed coronal emission. The reprocessed coronal emission is determined self-consistently, but the relative fractions of direct coronal and disc emissions are tuned to match the observations. A more developed theory of accretion discs would avoid this empirical tuning. That is, given the engine size (and perhaps spin) parameters along with the accretion rate, a developed theory of accretion would allow us to directly predict the relative contributions of the emitting components, and thus the overall spectra, from first principles. We would ultimately like such a theory to also predict the time evolution of these relative fractions as the accretion rate evolves.

2.2. Standard disc theory does not distinguish local versus non-local transport

The standard SS73 ‘viscosity’ formalism is elegantly suited for capturing the local radial transport of angular momentum within the disc but the inclusion of non-local and vertical transport requires further development. Among the important legacies of SS73 is thus the opportunity for further work to build a more inclusive, but still practical framework.

To explain this further we emphasize that the SS73 formalism makes sense only as a mean field theory. This follows because all dynamical quantities such as velocity and density are taken to be axisymmetric, and thus a function only of radius. But if angular momentum is transported by a turbulent viscosity, axisymmetry cannot apply on arbitrarily small scales. Thus the equations for standard disc theory must emerge from mean field equations. To derive the theory, dynamical quantities such as velocity, magnetic field, density and pressure must be divided into mean and ‘fluctuating’ quantities, with suitably chosen averages of the latter vanishing. While this procedure is not formally shown in SS73, a first principles approach using different averaging procedures (Balbus & Hawley Reference Balbus and Hawley1998; Kuncic & Bicknell Reference Kuncic and Bicknell2004) reveals that the turbulent viscosity emerges from a specific closure approximation to the correlated stresses associated with magnetic and velocity fluctuations.

To see more explicitly the origin and limitations of turbulent viscosity, note that the rate of change of vertical angular momentum for a magnetohydrodynamic disc can be obtained by multiplying the azimuthal ${\it\phi}$ component of momentum density by the radius in cylindrical coordinates. This gives (Balbus & Hawley Reference Balbus and Hawley1998)

(2.1) $$\begin{eqnarray}\displaystyle & & \displaystyle \partial _{t}({\it\rho}rv_{{\it\phi}})+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\left[r{\it\rho}v_{{\it\phi}}\boldsymbol{v}-r\frac{\boldsymbol{B}B_{{\it\phi}}}{4{\rm\pi}}+r\left(P_{th}+\frac{B_{r}^{2}+B_{z}^{2}}{8{\rm\pi}}\right)\hat{\boldsymbol{e}}_{{\it\phi}}\right]\nonumber\\ \displaystyle & & \displaystyle \quad -\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\left[\frac{r{\it\nu}_{V}\hat{\boldsymbol{e}}_{{\it\phi}}}{3}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{v}+r^{2}{\it\nu}_{V}\boldsymbol{{\rm\nabla}}\frac{v_{{\it\phi}}}{r}\right]=0,\end{eqnarray}$$

where ${\it\rho}$ is the density, $P_{th}$ is the thermal pressure, $r$ is the radius, $\boldsymbol{B}=(B_{r},B_{{\it\phi}},B_{z})$ is the magnetic field, $\boldsymbol{v}=(v_{r},v_{{\it\phi}},v_{z})$ is the velocity, ${\it\nu}_{V}$ is the dynamic viscosity and $\hat{\boldsymbol{e}}_{{\it\phi}}$ is a unit vector in the ${\it\phi}$ direction. We assume radiation pressure is subdominant.

The sum of the square bracketed quantities is the total angular momentum flux. We are interested in the transport of angular momentum radially and vertically but not azimuthally. Since the divergence of $\hat{\boldsymbol{e}}_{{\it\phi}}$ zero, any terms multiplying $\hat{\boldsymbol{e}}_{{\it\phi}}$ (including the contributions to the second and third terms from the ${\it\phi}$ vector components) make no contribution to the radial or vertical fluxes. The terms proportional to ${\it\nu}_{V}$ are microphysical dissipation terms and we assume here that they are negligible for high Reynolds number and magnetic Reynolds number astrophysical contexts. Removing the ${\it\phi}$ terms and the microscopic flux terms leaves the vertical and radial fluxes associated with macroscopic flows, namely

(2.2) $$\begin{eqnarray}F_{z}=r{\it\rho}[v_{{\it\phi}}v_{z}-v_{Az}v_{A{\it\phi}}],\end{eqnarray}$$

and

(2.3) $$\begin{eqnarray}F_{r}=r{\it\rho}[v_{{\it\phi}}v_{r}-v_{Ar}v_{A{\it\phi}}],\end{eqnarray}$$

where we have defined the Alfvén speed as $\boldsymbol{v}_{A}=(v_{Ar},v_{A{\it\phi}},v_{Az})=((B_{r}/\sqrt{4{\rm\pi}{\it\rho}}),(B_{{\rm\pi}}/\sqrt{4{\rm\pi}{\it\rho}}),(B_{z}/\sqrt{4{\rm\pi}{\it\rho}}))$ . So far, all quantities can depend on $r,{\it\phi}$ and $z$ and it is only after a formal averaging that we obtain a reduction of independent variables. Different procedures involving some combination of ensemble, temporal or spatial averages can be used. For compressible flows, a convenient approach is to use density-weighted time averages and density-weighted azimuthal averages (Kuncic & Bicknell Reference Kuncic and Bicknell2004). Averages of density are taken to be time and azimuthally averaged so that splitting ${\it\rho}$ into mean and fluctuating components gives

(2.4) $$\begin{eqnarray}{\it\rho}=\overline{{\it\rho}}+{\it\rho}^{\prime },\end{eqnarray}$$

where

(2.5) $$\begin{eqnarray}\langle {\it\rho}^{\prime }\rangle =\overline{{\it\rho}^{\prime }}=0\quad \text{and}\quad \langle {\it\rho}\rangle =\overline{{\it\rho}}=\frac{\displaystyle \int _{0}^{t_{M}}\!\int \!{\it\rho}(r,{\it\phi},z,t+t^{\prime })\,\text{d}{\it\phi}\,\text{d}t^{\prime }}{2{\rm\pi}t_{M}},\end{eqnarray}$$

where $t_{M}$ is a long time scale compared to turbulent turnover time scales. Averages of velocities, Alfvén speeds and their higher order moments, are taken to be density-weighted azimuthal and time averages and indicated by a ‘tilde’. For example, for the velocity we have

(2.6) $$\begin{eqnarray}\tilde{\boldsymbol{v}}\equiv \langle \boldsymbol{v}\rangle =\frac{\displaystyle \int _{0}^{\text{T}}\!\int \!{\it\rho}(r,{\it\phi},z,t+t^{\prime })\boldsymbol{v}(r,{\it\phi},z,t+t^{\prime })\,\text{d}{\it\phi}\,\text{d}t^{\prime }}{2{\rm\pi}T\,\overline{{\it\rho}}}\end{eqnarray}$$

with the assumption that the density and velocity fluctuations are uncorrelated, i.e. $\langle {\it\rho}^{\prime }\boldsymbol{v}^{\prime }\rangle =0$ . For the combination of density and two velocities for example, this averaging procedure then gives

(2.7) $$\begin{eqnarray}\langle {\it\rho}v_{r}v_{z}\rangle =\overline{{\it\rho}}\tilde{v}_{r}\tilde{v}_{z}+\overline{{\it\rho}}\langle {v^{\prime }}_{r}{v^{\prime }}_{z}\rangle ,\end{eqnarray}$$

where we have made the additional approximation that triple correlations vanish to obtain the last term. Then $\langle {\it\rho}{v^{\prime }}_{r}{v^{\prime }}_{z}\rangle =\overline{{\it\rho}}\langle {v^{\prime }}_{r}{v^{\prime }}_{z}\rangle$ after using (2.4).

Applying the above averaging procedure to (2.2) and (2.3) gives

(2.8) $$\begin{eqnarray}\langle F_{z}\rangle =r\overline{{\it\rho}}[\tilde{v}_{{\it\phi}}\tilde{v}_{z}-\tilde{v}_{A{\it\phi}}\tilde{v}_{Az}+\langle {v_{{\it\phi}}}^{\prime }v_{z}^{\prime }\rangle -\langle v_{A{\it\phi}}^{\prime }v_{Az}^{\prime }\rangle ],\end{eqnarray}$$

and

(2.9) $$\begin{eqnarray}\langle F_{r}\rangle =r\overline{{\it\rho}}[\tilde{v}_{{\it\phi}}\tilde{v}_{r}-\tilde{v}_{A{\it\phi}}\tilde{v}_{Ar}+\langle {v_{{\it\phi}}}^{\prime }v_{r}^{\prime }\rangle -\langle v_{A{\it\phi}}^{\prime }v_{Ar}^{\prime }\rangle ].\end{eqnarray}$$

These mean fluxes are functions of $r$ , $z$ and coarse grained in $t$ .

To further simplify without completely eliminating the $z$ dependence, we integrate separately over hemispheres. This additional coarse graining is indicated by the subscript ‘ $\overline{{\it\rho}}$ ’. Thus $\langle \tilde{v}_{z}\rangle _{\overline{{\it\rho}}}\equiv \int _{0}^{H/2}\overline{{\it\rho}}\tilde{v}_{z}\,\text{d}z/({\it\Sigma}/2)$ , where ${\it\Sigma}/2=\int _{0}^{H/2}\overline{{\it\rho}}\,\text{d}z$ and $H$ is the disc thickness defined by a density scale height. We also assume that the three quantities $\tilde{v}_{{\it\phi}}$ , $\tilde{v}_{r}$ and $\tilde{v}_{A{\it\phi}}$ are independent of $z$ . The results of the $z$ integration of (2.8) and (2.9) are then

(2.10) $$\begin{eqnarray}\langle F_{z}\rangle _{\pm }=\frac{{\it\Sigma}r}{2}[\tilde{v}_{{\it\phi}}\tilde{v}_{z\pm ,\overline{{\it\rho}}}-\tilde{v}_{A{\it\phi},}\tilde{v}_{Az\pm ,\overline{{\it\rho}}}+\langle v_{{\it\phi}}^{\prime }v_{z}^{\prime }\rangle _{\pm ,\overline{{\it\rho}}}-\langle v_{A{\it\phi}}^{\prime }v_{Az\pm ,\overline{{\it\rho}}}^{\prime }\rangle _{\pm ,\overline{{\it\rho}}}],\end{eqnarray}$$

and

(2.11) $$\begin{eqnarray}\langle F_{r}\rangle _{\pm }=\frac{{\it\Sigma}r}{2}[\tilde{v}_{{\it\phi}}\tilde{v}_{r}-\tilde{v}_{A{\it\phi}}\tilde{v}_{Ar\pm ,\overline{{\it\rho}}}+\langle v_{r}^{\prime }v_{{\it\phi}}^{\prime }\rangle _{\pm ,\overline{{\it\rho}}}-\langle v_{Ar}^{\prime }v_{A{\it\phi}}^{\prime }\rangle _{\pm ,\overline{{\it\rho}}}],\end{eqnarray}$$

where the $\pm$ subscripts indicate the values in the ‘ $+$ ’ (upper) and ‘ $-$ ’ (lower) hemispheres, but there is otherwise no $z$ dependence left. In each hemisphere, the quantities remain functions of radius and coarse grained time. The stresses presented in this form offer a unifying focal point for a variety of generalizations of SS73, and keeping the hemispheres averaged separately is relevant for tracking any pseudo scalars (such as helicities) or fluxes into coronae that might play a dynamical role.

To see SS73 emerge from these equations, we take $\tilde{v}_{z}=\tilde{v}_{Az}=\tilde{v}_{Ar}=0$ , and ignore vertical stresses from fluctuations. Then, only the radial flux from (2.10) matters. If we further assume the $\pm$ contributions are equal (symmetric across the mid plane) then

(2.12) $$\begin{eqnarray}\langle F_{r}\rangle =\langle F_{r}\rangle _{+}+\langle F_{r}\rangle _{-}={\it\Sigma}r[\tilde{v}_{{\it\phi}}\tilde{v}_{r}+T_{r{\it\phi}}],\end{eqnarray}$$

where $T_{r{\it\phi}}$ is now the stress combination from fluctuations given by

(2.13) $$\begin{eqnarray}T_{r{\it\phi}}=\langle v_{r}^{\prime }v_{{\it\phi}}^{\prime }\rangle -\langle v_{Ar}^{\prime }v_{A{\it\phi}}^{\prime }\rangle \simeq -{\it\nu}_{T}r\partial _{r}{\it\Omega}\sim q{\it\alpha}_{ss}c_{s}^{2},\end{eqnarray}$$

where the penultimate similarity (Lynden-Bell & Pringle Reference Lynden-Bell and Pringle1974) represents the closure that converts a sum of second order fluctuation correlations into a term linear in the angular velocity shear. The latter similarity follows from setting ${\it\nu}_{T}\equiv {\it\alpha}_{ss}c_{s}H$ , where ${\it\alpha}_{ss}$ is the dimensionless SS73 ‘viscosity’ coefficient, along with ${\it\Omega}\simeq c_{s}/H$ in hydrostatic equilibrium and ${\it\Omega}\propto r^{-q}$ , where $q=3/2$ for Keplerian flow. Although SS73 warned that ${\it\alpha}_{ss}$ could be a function of radius, it has become common practice to employ the simplest model in which ${\it\alpha}_{ss}$ is a constant.

We can now recover a standard form of the accretion rate balance equation in the steady state. Applying Gauss’s theorem to (2.1) in cylindrical coordinates with only a radial derivative, and integrating the flux divergence when vertical stresses are ignored implies that $r\langle F_{r}\rangle$ is constant. If we further assume that stresses from fluctuations vanish at the inner boundary $r=r_{0}$ , then using $r\langle F_{r}\rangle |_{r=r_{0}}=r\langle F_{r}\rangle$ in (2.12) along with ${\dot{M}}=-2{\rm\pi}{\it\Sigma}r\tilde{v}_{r}$ for the steady state accretion rate, we obtain

(2.14) $$\begin{eqnarray}-\frac{{\dot{M}}r_{0}^{2}{\it\Omega}_{0}}{2{\rm\pi}}=-\frac{{\dot{M}}r^{2}{\it\Omega}}{2{\rm\pi}}+r^{2}{\it\Sigma}T_{r{\it\phi}}.\end{eqnarray}$$

Thus

(2.15) $$\begin{eqnarray}T_{r{\it\phi}}=\frac{{\dot{M}}{\it\Omega}}{2{\rm\pi}{\it\Sigma}}\left(1-\frac{r_{0}^{2-q}}{r^{2-q}}\right),\end{eqnarray}$$

for $\tilde{v}_{{\it\phi}}={\it\Omega}r$ and ${\it\Omega}={\it\Omega}_{0}(r/r_{0})^{-q}$ . A positive stress implies positive accretion for $q<2$ . For $q>2$ , this equation alone would suggest excretion but for this regime epicyclic frequencies become imaginary and the orbits and disc are unstable, violating the validity of the steady state assumed above and invalidating the excretion interpretation. This is also why the shearing box model produces runaway mean velocities for $q>2$ and is unsuitable for that regime (Nauman & Blackman Reference Nauman and Blackman2014).

The above derivation highlights that a physically consistent derivation of the standard ${\it\alpha}_{ss}$ viscous disc accretion model from the MHD equations requires the following: (i) formal averaging that leaves the mean field quantities dependent only on radius; (ii) dropping vertical transport terms associated with large-scale fields; (iii) dropping vertical transport terms associated small-scale MHD fluctuations; (iv) dropping radial transport terms associated with large-scale magnetic fields; and (v) including the radial transport from fluctuations by replacing the sum of their stress contributions by a viscous term linear in the mean angular velocity gradient. There is a significant opportunity to improve this standard theory by relaxing even one of the above assumptions. The importance of large-scale fields, as evidenced by jets or coronae, often invalidates at least assumption (ii) above. In addition, vertical stresses associated with fluctuations can be smaller than the radial stresses by a factor $H/r$ but still contribute equally to the torque when both radial and vertical torques are present, violating assumption (iii). Radial transport by large-scale loops is excluded by assumption (iv), which may also be a significant restriction (see also Rüdiger Reference Rüdiger1987). Note that even if stresses from large-scale structures might be written as proportional to $c_{s}^{2}$ by analogy to that of viscous transport, the dimensionless coefficient can be different for that contribution, corresponding to large-scale structures which fuel coronae and outflows compared to local viscous transport which supplies heating within the disc.

The standard SS73 closure modelling the transport as purely viscous is simple and elegant, but there have been early hints that it may be misleading when applied to MRI driven turbulence (Balbus & Hawley Reference Balbus and Hawley1998). There have been specific attempts to improve the closure (Ogilvie Reference Ogilvie2003; Pessah, Chan & Psaltis Reference Pessah, Chan and Psaltis2006), and more recent direct evidence that the MRI stress does not agree with the SS73 formulation (Pessah, Chan & Psaltis Reference Pessah, Chan and Psaltis2008; Nauman & Blackman Reference Nauman and Blackman2015). As discussed further below, our interpretation of numerical evidence further suggests that the transport from the MRI is actually dominated by large scales.

How does previous work on accretion disc theory connect to relaxing assumptions (i)–(v)? Blandford & Payne (Reference Blandford and Payne1982) and Königl (Reference Königl1989) exemplify models for which the primary angular momentum transport occurs via large-scale magnetic fields. Such models implicitly include (i), (ii) and (iv), but not (iii) and (v). Blandford & Begelman (Reference Blandford and Begelman1999) and Quataert & Narayan (Reference Quataert and Narayan1999) are examples of practical models which include mass loss to reduce disc luminosity, and thus implicitly the role of vertical transport, but the associated stresses are not formally derived from a closure nor are the outflows emergent from equations that include large-scale magnetic field dynamics. The models are rooted in the viscous formalism of assumption (v), with the additional assumption that in collisionless plasmas, the dissipation may heat primarily ions not electrons, with the coupling of the two being limited primarily by Coulomb collisions. Understanding the physics of this dissipation (Blackman Reference Blackman1999) and energy transfer is itself an important area of active research at the small-scale collisionless frontier of accretion theory (Quataert, Heinemann & Spitkovsky Reference Quataert, Heinemann and Spitkovsky2014), but is outside of the scope of the present paper.

Rüdiger, Elstner & Schultz (Reference Rüdiger, Elstner and Schultz1993), Campbell & Caunt (Reference Campbell and Caunt1999) and Campbell (Reference Campbell2000) developed models for accretion that transfer angular momentum via large-scale ${\it\alpha}{-}{\it\Omega}$ dynamo-produced fields, and showed that this contribution to transport can be comparable to that from turbulent transport of the SS73 type. Campbell (Reference Campbell2003) solved for the vertical structure and radial structure in such discs. Taken together, these papers maintain assumption (v) above but relax the others. They should be revisited with an eye toward reducing the complexity and incorporating new developments from large-scale dynamo theory (see §§ 4.2 and 4.3). A separate important step towards a first principles improvement to mean field disc theory is Kuncic & Bicknell (Reference Kuncic and Bicknell2007) which adopts (i), (ii), (iv) and (v) but drops (iii) by including some vertical transport. Campbell & Caunt (Reference Campbell and Caunt1999) and Campbell (Reference Campbell2000, Reference Campbell2003) did not compute the observed spectral influence of the vertical transport and Kuncic & Bicknell (Reference Kuncic and Bicknell2004, Reference Kuncic and Bicknell2007) did so partially but only by removing it from the contribution to the disc spectra. Overall, there remains much opportunity for further work.

2.3. Averaging and precision of disc theory

That standard ${\it\alpha}_{ss}$ accretion disc theory is a mean field theory means that the precision of the predictive power is limited in spatial, temporal and spectral resolution. Without quantifying this precision, and recognizing that it is finite, a disagreement between theory and observed data can be misinterpreted to imply the need for new physics rather than as a consequence of statistical fluctuations within the predicted precision error.

To obtain the standard link between accretion and stress in the previous section, we considered a time average and a $z$ -integration without averaging in $r$ . The longer the duration over which the time average is taken, the lower the predicted error in precision. In comparing theory with data, it is important to know the time scales over which the data are binned to be sure that appropriate theoretical averages are being used to compare to appropriately averaged data. The smaller the scales of turbulent fluctuations compared to the smoothing scales, the more agreement (and thus freedom in choice) there is between temporal, ensemble and quasi-local spatial averages. However, if fluctuations are sufficiently large, then there may be a mismatch between predictions of the theory and observations taken over very short integration times, even if the theory would correctly explain the data averaged over much longer times. The case of YSO is particularly noteworthy because data against which theory is often compared are integrated for time scales much smaller than an orbit time. Some work on this mismatch between the statistically expected deviations from the predicted mean field theory has been done (Blackman Reference Blackman1998; Blackman, Nauman & Edgar Reference Blackman, Nauman and Edgar2010) but more is needed to address the non-local contribution of emission to a given wavelength.

A related subtlety is that although variability and deviations from standard predictions of axisymmetric theory may arise from systematic phenomena such as warps (Hartnoll & Blackman Reference Hartnoll and Blackman2000; Tremaine & Davis Reference Tremaine and Davis2014), when the theory itself uses a viscous model for transport, self-consistency requires an assessment that the predicted structures are large compared to scales of the statistical fluctuations over which the mean field theory is averaged.

2.4. Local radial transport is insufficient: vertical and non-local transport must be included

Although most accretion disc modelling presumes a viscosity transport model and a constant ${\it\alpha}_{ss}$ , observations and simulations challenge this minimalist paradigm. Evidence for its insufficiency emerges from: (i) observations indicating non-local radial and vertical transport via jets, outflows and coronae and the important role of large-scale fields in plausible jet models; (ii) inadequacy of total stress magnitudes derived from local shearing box MRI simulations if taken at face value; (iii) dominance of large-scale stresses in local and global simulations; (iv) contradiction between transport computed from simulations and the dependence on shear predicted from the constant ${\it\alpha}_{ss}$ viscosity model; and (v) evidence that ${\it\alpha}_{ss}$ depends on radius in some global simulations (Penna, Narayan & Sa̧dowski Reference Penna, Narayan and Sa̧dowski2013). We discuss points (i)–(iv) further below.

4.1. Understanding the MRI in shearing boxes as a transfer of linear momentum

A widely used physical explanation for the linear MRI appeals to two masses tethered by a spring separated by an initial radial displacement in a background rotation profile of outwardly increasing angular momentum (Balbus & Hawley Reference Balbus and Hawley1998). For a weak enough spring, the inner mass loses angular momentum to the outer mass and will continue to move toward loci of lower background angular momentum, namely, inward. Complementarily, the outer mass moves toward regions of higher background angular momenta, namely outward. By considering the equations for each mass separately and taking their difference, the expression for the separation between the masses can be obtained with exponential growth being the signature of the instability. A similar approach can be used to argue that a flow with an outwardly decreasing background angular momentum gradient should be stable to the MRI. Such stability is consistent with recent generalizations Shakura & Postnov (Reference Shakura and Postnov2014).

Because the Cartesian shearing box has an infinite radius of curvature, the angular momentum is ill-defined therein. But the instability still occurs because it manifests in the linear momentum equations. Angular momentum transport is a secondary consequence that arises for a finite radius of curvature. Although these points follow simply from the known equations, the standard physical explanation of the MRI given in the previous paragraph with angular momentum commonly presented without clarification for a shearing box. For clarification, we discuss the tethered spring analogue in more detail below.

To apply the viscous closure (2.13) to that of a shearing box, we shift to the local Hill frame Cartesian coordinates for a disc co-rotating with local angular speed ${\it\Omega}_{0}={\it\Omega}(r=r_{0})$ , assumed to represent the range $r_{0}-x\leqslant r\leqslant r_{0}+x$ , with $x\ll r_{0}$ . Ignoring curvature terms, the mean $y$ -velocity in the box equals

(4.1) $$\begin{eqnarray}\tilde{v}_{{\it\phi}}-{\it\Omega}_{0}r_{0}=r({\it\Omega}-{\it\Omega}_{0})\simeq xr\partial _{r}{\it\Omega}|_{r=r_{0}}=-xq{\it\Omega}\end{eqnarray}$$

to lowest order in $x$ so that

(4.2) $$\begin{eqnarray}\langle T_{r{\it\phi}}\rangle ={\it\alpha}c_{s}Hq{\it\Omega}_{0}\end{eqnarray}$$

in place of (2.13). From the right side of (4.1) we see that the centreline of the box $x=0$ represents the co-rotation radius whose mean azimuthal ( $y$ -direction) velocity vanishes.

In the absence of pressure gradients, the Hill frame equations capturing the MRI instability (Balbus Reference Balbus2003) in the presence of a vertical field of strength $B_{z}$ , subjected to displacement ${\bf\xi}=(x(t)\text{e}^{\text{i}kz},y(t)\text{e}^{\text{i}kz})$ about the vertical line $(x=0,y=0)$ can be written

(4.3) $$\begin{eqnarray}\ddot{x}-2{\it\Omega}{\dot{y}}=-(K_{A}-T)x,\end{eqnarray}$$

and

(4.4) $$\begin{eqnarray}\ddot{y} +2{\it\Omega}{\dot{x}}=-K_{A}y,\end{eqnarray}$$

where $T=-r\,\text{d}{\it\Omega}^{2}/\text{d}r=+2q{\it\Omega}^{2}$ ( $q=3/2$ for Keplerian) is the coefficient of the surviving tidal force in the Hill frame; $K_{A}=(k\tilde{v}_{zA})^{2}$ , where $\tilde{v}_{zA}$ is the Alfvén speed associated with the vertical field; the second terms on the left sides come from the Coriolis force. The appearance of $K_{A}$ arises as a consequence of magnetic tension force because the linear fluctuation in the field ${\it\delta}\boldsymbol{B}$ satisfies $\partial _{t}{\it\delta}\boldsymbol{B}\sim \text{i}kB{\it\delta}\boldsymbol{v}$ from the induction equation, so that ${\it\delta}\boldsymbol{B}\sim \text{i}kB{\it\delta}\boldsymbol{v}t\sim \text{i}kB{\bf\xi}$ on time integrating. The curvature force appearing on the right side of the momentum equation is then $(\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}/4{\rm\pi}{\it\rho}){\it\delta}\boldsymbol{B}\sim -K_{A}{\bf\xi}$ with the magnetic pressure having been absorbed into the total pressure, assumed to have a vanishing gradient.

Since the $z$ dependence dropped out of (4.3) and (4.4), the relevant motion takes place in the $x,y$ plane. The equations can then be interpreted as describing the evolution of the position of a small parcel of fluid at position ${\bf\xi}$ tethered by a spring to a fixed point at the origin $x,y=0$ . Given that $K_{A}$ and $T$ are constants, the right hand side of the above equations then represents two directional components of a spring force per unit mass with distinct spring constants $K_{A}-T$ and $K_{A}$ , respectively. The spring constant in (4.4) has a positive sign, such that positive displacement strengthens the force opposing the displacement. However, if $T>K_{A}$ then (4.3) would have a negative spring constant, so the associated force would exacerbate the separation of the fluid parcel from the origin as the separation increases in the $x$ direction. When interpreted as a runaway displacement of the fluid parcel from its initial position, instability requires $T>K_{A}$ and stability requires $T\leqslant K_{A}$ as found by solving the two equations. (Note that for $q<0$ , the azimuthal velocity increases outward. Then $T<0$ , and the effective spring constant in (4.3) would be positive, implying stable oscillations.)

Equations (4.3) and (4.4) can also be adjusted to describe the displacement of two equal density fluid parcels mutually attached by a spring. To do so, we must construct two sets of equations similar to (4.3) and (4.4), one for each of the two masses, and then subtract them to obtain the equations of evolution of their separation. However, note that while the tidal force has the same form as in the above equations, the force from the $K_{A}$ terms depend on the separation between the two masses. The equations for each of the two parcels thus take the form

(4.5) $$\begin{eqnarray}\ddot{x}_{i}-2{\it\Omega}{\dot{y}}_{i}=Tx_{i}-K_{A}(x_{i}-x_{j}),\end{eqnarray}$$

and

(4.6) $$\begin{eqnarray}\ddot{y} _{i}+2{\it\Omega}{\dot{x}}_{i}=-K_{A}(y_{i}-y_{j}),\end{eqnarray}$$

where $i=1,j=2$ represent the equations of motion for parcel 1 and $i=2,j=1$ represents the equations for parcel 2. Subtracting the respective $x$ and $y$ equations for parcel 2 from those of parcel 1, and letting $({\it\zeta}_{x},{\it\zeta}_{y})=(x_{1}-x_{2},y_{1}-y_{2})$ gives

(4.7) $$\begin{eqnarray}\ddot{{\it\zeta}}_{x}-2{\it\Omega}\dot{{\it\zeta}}_{y}=-(2K_{A}-T){\it\zeta}_{x},\end{eqnarray}$$

and

(4.8) $$\begin{eqnarray}\ddot{{\it\zeta}}_{y}+2{\it\Omega}\dot{{\it\zeta}}_{x}=-2K_{A}{\it\zeta}_{y}.\end{eqnarray}$$

Instability as interpreted as separation of the two connected masses requires $T>2K_{A}$ .

Neither of the above explanations of the instability criteria require explicit mention of the angular momentum. The instability is evident in the linear momentum equation and the stress that it produces in the shearing box limit represents a transfer of linear momentum. For a finite radius of curvature, angular momentum transport then becomes a physical consequence.

4.2. Stress spectra from shearing boxes: the dominance of large scales

4.2.1. Stress spectra from fluctuations

The simplifications and limitations of the shearing box are significant compared to real discs, but as with the triply periodic box commonly used in forced MHD turbulence studies, the shearing box should be thought of as its own well-defined computational physics experiment. If we can understand how fields and flows evolve in such suitably chosen minimalist systems, we will hopefully gain insight into key elements of real systems.

How well are the basic assumptions of the SS73 paradigm confirmed or refuted by shearing box studies? This is germane given that stress magnitudes from shearing box simulations are widely quoted as inputs for spectral models. Shearing boxes were designed with the intuition that angular momentum transport is a primarily local and radial transport phenomena. In this context we can ask: (i) What scales dominate the stresses? (ii) Do the radial stresses depend linearly on shear as predicted by the SS73 framework with constant ${\it\alpha}_{ss}$ ? As we will see below, the answers are: (i) large (ii) no.

We first determine the scale dependence of stresses from from stratified shearing box simulations by computing the $k$ -dependent cumulative stress from fluctuations (Nauman & Blackman Reference Nauman and Blackman2014). This is computed by integrating from the minimum wavenumber up to $k$ and determining how large $k$ must be in the integral to capture most of the stress in the box. We have used the ATHENA code and solved the vertically stratified MHD equations without explicit diffusivities, for an isothermal equation of state. The boundary conditions are shear periodic in the radial $x$ direction and periodic in the azimuthal $y$ and vertical $z$ directions (with gravity thus taken to vanish at vertical boundaries). For the results shown in figure 2, we employed a box domain size of $4H\times 4H\times 8H$ where $H$ is the density scale height, with a resolution of 96 zones per $H$ and an initial plasma ${\it\beta}$ (thermal to magnetic pressure ratio) equal to 100 everywhere and an initially toroidal mean field which decreases in strength vertically like the density.

We obtain the cumulative stresses from fluctuations as a function of wavenumber by first computing 2-D Fourier transforms $f(\boldsymbol{k},z)=f(k_{x},k_{y},z)$ of a physical quantity $F(\boldsymbol{x})$ (which was either the magnetic field, off-diagonal magnetic stress, momentum density or off-diagonal Reynolds stress) in variables $x$ and $y$ . For each quantity we then computed three different integrals of the quadratic norm: (i) circle average, $|f(k,z)|^{2}=\int |f(\boldsymbol{k}^{\prime })|^{2}{\it\delta}(|\boldsymbol{k}^{\prime }|-k)\,\text{d}\boldsymbol{k}^{\prime }$ , where $|\boldsymbol{k}|=\sqrt{k_{x}^{2}+k_{y}^{2}}$ ; and (ii) $y$ -integrated, $|f(k_{x},z)|^{2}=\int \!|f(\boldsymbol{k})|^{2}\,\text{d}k_{y}$ , (iii) $x$ -integated $|f(k_{y},z)|^{2}=\int \!|f(\boldsymbol{k})|^{2}\,\text{d}k_{x}$ . To eliminate $z$ , we then averaged the resulting 1-D spectra over $z$ ranges specified by either ‘disc’ or ‘corona’ as determined by the height above which the saturated ${\it\beta}$ is above or below unity, respectively. The cumulative stress and energy spectra then represent the integral of a given spectrum up to wavenumber $k$ divided by the integral over the full range of $k$ , that is:

(4.9) $$\begin{eqnarray}\left.\begin{array}{@{}l@{}}Q_{f}(k)=\frac{\displaystyle \int _{1}^{k}|f(k^{\prime })|^{2}2{\rm\pi}k^{\prime }\,\text{d}k^{\prime }}{\displaystyle \int _{1}^{k_{max}}|f(k^{\prime })|^{2}2{\rm\pi}k^{\prime }\,\text{d}k^{\prime }},\\ Q_{f}(k_{x})=\frac{\displaystyle \int _{0}^{k_{y}}|f(\boldsymbol{k})|^{2}\,\text{d}k_{y}\left/\int _{}^{}\right.\,\text{d}k_{y}}{\displaystyle \int _{0}^{k_{y,max}}|f(\boldsymbol{k})|^{2}\,\text{d}k_{y}\left/\int _{}^{}\right.\text{d}k_{y}},\\ Q_{f}(k_{y})=\frac{\displaystyle \int _{0}^{k_{x}}|f(\boldsymbol{k})|^{2}\,\text{d}k_{x}\left/\int _{}^{}\,\text{d}k_{x}\right.}{\displaystyle \int _{0}^{k_{x,max}}|f(\boldsymbol{k})|^{2}\,\text{d}k_{x}\left/\int _{}^{}\right.\,\text{d}k_{x}}.\end{array}\right\}\end{eqnarray}$$

The stresses computed as above do not include contributions for which both $k_{y}=0$ and $k_{x}=0$ , though each can separately vanish. We discuss the $k_{y}=k_{x}=0$ contributions in § 4.2.2.

Figure 2 shows that, for both coronae and disc regions, more than half of the radial and vertical stresses are accounted for by contributions from the first few wavenumbers. We found that this conclusion is independent of the box size for the domain sizes studied. The cumulative stresses are dominated by correlations from the large scales of the box. Even though stresses are nearly constant, as the box size and resolution change (a necessary condition for convergence), the upper integration bound on the wavenumber needed to capture 50 % of the cumulative stress does not increase as the box domain size increases. This suggests that the local simulations are insufficient to capture significant large-scale contributions to the transport. That in turn highlights a limitation of the shearing box; the box size cannot be increased indefinitely without having to incorporate curvature terms or vertical shear (McNally & Pessah Reference McNally and Pessah2014) to ensure relevance when scaled to a real system. That shearing boxes suggest non-local transport in the sense described is further supported by measurements from global simulations of energy spectra (Flock et al. Reference Flock, Dzyurkevich, Klahr, Turner and Henning2012; Suzuki & Inutsuka Reference Suzuki and Inutsuka2014) and stress spectra (Beckwith, Armitage & Simon Reference Beckwith, Armitage and Simon2011; Sorathia et al. Reference Sorathia, Reynolds, Stone and Beckwith2012; Parkin & Bicknell Reference Parkin and Bicknell2013).

While energy spectra and stresses for both disc and coronae (respectively distinguished by plasma ${\it\beta}$ greater than or less than unity) in cases with net flux and without explicit dissipation are dominated by large scales in both local and global simulations and the magnetic energy spectra generally do not show a turnover at large scales (Nauman & Blackman Reference Nauman and Blackman2014), the zero net flux, unstratified simulations of (Fromang Reference Fromang2010) with explicit dissipation do seem to show a turnover in the magnetic and kinetic energy spectra. The wavenumber of the magnetic peak seems to increase with increasing magnetic Reynolds number. A concern about the magnitude of the stress in previous unstratified zero flux simulations however is that they did not seem to converge Pessah, Chan & Psaltis (Reference Pessah, Chan and Psaltis2007), although Fromang (Reference Fromang2010) argues for convergence. Controversy still lingers over the question of convergence of the stress magnitude, even for stratified simulations (Davis et al. Reference Davis, Stone and Pessah2010; Bodo et al. Reference Bodo, Cattaneo, Mignone and Rossi2014). Stress spectra, in addition the stress magnitude, should be assessed in convergence studies.

Figure 2.

Modified from Nauman & Blackman (Reference Nauman and Blackman2014) (to include $(0,k_{y})$ and $(k_{x},0)$ modes): the fractional spectral power for several quantities as indicated in a Keplerian shearing box simulation (see text for parameters) versus a dimensionless measure of the wavenumber. Black curves are for corona $|z|\geqslant 2H$ and blue curves are for disc $|z|<2H$ and the legend shown applies to both colour sets of curves. The wavenumbers are of the form $k=2{\rm\pi}n/L$ where $L=4H$ is the domain size in the $x$ or $y$ direction and $n$ ranges from 1 to 96. The $x$ -axes are $kH/2{\rm\pi}=n/4$ . The total stresses are dominated by the Maxwell stresses and so lessons learned from their plots also apply for the total. Both coronae and disc show that the fractional power for all quantities is well over 50 % already for the first few wave modes. Fractional power is calculated using (4.9) using $|f(\boldsymbol{k})|^{2}=\{|B(\boldsymbol{k})|^{2},|\sqrt{{\it\rho}}v(\boldsymbol{k})|^{2},|B_{x}B_{y}(\boldsymbol{k})|^{2},|B_{y}B_{z}(\boldsymbol{k})|^{2}\}$ for magnetic energy, kinetic energy, $xy$ Maxwell stress and $yz$ Maxwell stress, respectively. For energy $Q_{f}$ gives the fractional power, but for the stresses it gives fractional power of stress squared, so we plotted $\sqrt{Q_{f}}$ in those cases. The fractional power distribution in $k_{y}$ is significantly different from that of $k_{x}$ and $k$ for all quantities, highlighting anisotropy.

4.2.2. Comparing stresses from mean fields to those from fluctuations

The analysis of the previous section does not include contributions to the stress from $k_{x}=k_{y}=0$ modes at fixed $z$ . We refer to these as stresses from the mean field. Here we compare the contribution from these stresses to those from the fluctuations. As the measure of stresses and energies from mean magnetic fields we compute $T_{M,ij}(z)=\langle B_{i}\rangle \langle B_{j}\rangle$ where the brackets indicate $x,y$ averages, with the product being left as a function of $z$ . We then average the product over $z$ and orbit times $101<t<220$ to obtain the total contribution $T_{M,ij}$ . For the measure of stress from fluctuations, we first compute the contribution from the total field $T_{ij}(z)\equiv \langle B_{i}B_{j}\rangle$ , averaged over $z$ and $101<t<220$ as above to obtain $T_{ij}$ , and then subtract the total contribution from the mean fields to obtain ${T_{F}}_{ij}=T_{ij}-T_{M,ij}$ .

Using the above procedure for a $q=3/2$ case from a run with 32 zones/ $H$ and domain size $H\times 2H\times 4H$ (Nauman & Blackman Reference Nauman and Blackman2015), we find that the ratio of the stresses from the mean field to that of the total to be $T_{M,xy}/T_{xy}\simeq 0.24$ and the ratio of magnetic energy from the mean field to that of the total to be $0.47$ . These indicate a significant contribution from the mean field, bolstering the overall message that large-scale fields and transport are important. The simulation from which we derived this result has no extended corona and so the stress ratio likely underestimates the importance of mean fields: figure 3 shows that already at the base of the corona the mean stresses are dominant.

Figure 3.

Colour gradient plot of the ratio of the Maxwell stress from the $x$ and $y$ averaged mean fields to to the $x$ and $y$ averaged total Maxwell stress $T_{M,xy}(z)/T_{xy}(z)$ as a function of $z$ and $t$ . The data come from a vertically stratified, vertically periodic simulation of domain size $H\times 2H\times 4H$ with 32 zones $\text{h}^{-1}$ , shear parameter $q=3/2$ and an initial toroidal field with initial plasma ${\it\beta}=100$ . The plot shows the dominance of the mean components to the stress near the top of the box. Values in excess of unity can arise because the stress is a signed quantity and the denominator $T_{xy}(z)$ includes some non-vanishing oppositely signed contributions from cross terms between mean and fluctuations that do not vanish. This simulation has no extended corona where plasma ${\it\beta}<1$ but the particular importance of the large-scale field where the density and saturated ${\it\beta}$ decrease near the top of the box is evident.

4.3. More deviations from transport as viscosity in MRI sims

Simulations in shearing boxes reveal other important deviations from the traditional constant ${\it\alpha}_{SS}$ prediction such as the dependence on the shear profile index $q$ . Equation (4.2) shows that the SS73 formalism combined with the assumption of constant ${\it\alpha}$ predicts a stress linear in $q$ , but the dependence on shear has been found to be strongly inconsistent with such a linear dependence from both unstratified (Pessah et al. Reference Pessah, Chan and Psaltis2008) and stratified (Abramowicz, Brandenburg & Lasota Reference Abramowicz, Brandenburg and Lasota1996; Nauman & Blackman Reference Nauman and Blackman2015) simulations.

In addition, stress and magnetic energy spectra are anisotropic, being elongated in the $y$ direction and more so on large scales than small scales (see figure 2). This suggests that a nonlinear cascade ensues with turnover time scales that progressively decrease with scale compared to the shear time scale, the latter being constant on all scales.

The subtleties in anisotropy and scale dependence further highlight that stress magnitude should not be the sole criterion for convergence in MRI simulations. Convergence studies should also address the stress spectra and anisotropy. There remains opportunity for a systematic assessment of their convergence as a function of initial conditions, presence or absence of explicit dissipation, magnetic Prandtl number, resolution and local versus global simulations.

4.4. Saturation of MRI and tilt angle

The saturation mechanism of MRI stress in shearing boxes is still uncertain and may even differ from that of global simulations or real discs. We have emphasized that saturated stress spectra as a function of resolution and domain size should be evaluated in addition to stress magnitude as part of generalized convergence studies. However, an interesting nugget of universality (for simulations that are not radiation dominated) is a nearly constant value of $0.1\lesssim {\it\alpha}_{ss}{\it\beta}\lesssim 0.4$ . Blackman, Penna & Varnière (Reference Blackman, Penna and Varnière2008). Here the dimensionless measure of the stress associated with fluctuations in the SS73 formalism is ${\it\alpha}_{ss}=T_{xy}/P_{th}$ and ${\it\beta}=P_{th}/P_{mag}$ , the ratio of thermal to magnetic pressure. (for radiation dominated accreters this would be generalized to include radiation pressure and thus the constancy of ${\it\alpha}_{ss}{\it\beta}$ might be replaced by ${\it\alpha}_{ss}{\it\beta}_{tot}$ where ${\it\beta}_{tot}\sim (P_{rad}+P_{part})/P_{mag}$ ). As we show below, the near constancy of ${\it\alpha}_{ss}{\it\beta}$ is equivalent to a constancy of the tilt angle between the azimuthal and radial fluctuating field components, and can be compared to the associated tilt angle of the mean field if planar averages are considered for the latter. Computed as box averages, stresses from mean fields would not survive in the absence of a net flux, however a given $z$ , mean stresses can be important, as seen in figure 3. Below, we estimate and compare the tilt angles for fluctuations and planar averaged mean fields.

4.4.1. Tilt angle from fluctuating fields

Because the magnetic stress dominates the Reynolds stress we have for the box averaged contribution from fluctuations

$$\begin{eqnarray}{\it\alpha}_{ss}\simeq \frac{\langle b^{2}\rangle }{P}\frac{\langle b_{x}b_{y}\rangle }{\langle b^{2}\rangle }={\it\beta}^{-1}\frac{\tan {\it\theta}}{(1+\tan ^{2}{\it\theta}+b_{z}^{2}/b_{x}^{2})}\sim {\it\beta}^{-1}\tan {\it\theta},\end{eqnarray}$$

where $\tan {\it\theta}=b_{y}/b_{x}$ , and we assume $|b_{z}/b_{x}|\ll 1$ and $\tan {\it\theta}\ll 1$ for the last similarity (justified later).

The constancy of tilt angle across simulations can be explained as a consequence of the fact that the toroidal field is amplified above the radial field by shear over a time scale approximately equal to a correlation time scale (Blackman et al. Reference Blackman, Penna and Varnière2008). Roughly, $|b_{y}|\sim |b_{x}|(1+q{\it\Omega}{\it\tau}_{c})$ implying that $\tan {\it\theta}\sim (1+q{\it\Omega}{\it\tau}_{c})^{-1}$ . Nauman & Blackman (Reference Nauman and Blackman2015) showed that the tilt angle is also independent of the shear parameter $q$ using ATHENA because an increase or decrease in shear is compensated by a corresponding decrease or increase in the relevant correlation time. Whether the source of the reduced correlation time is physical or a consequence of the numerical method (such as the remap time scale across the radial boundary) needs further investigation.

4.4.2. Tilt angle from mean fields and comparison to that from fluctuations

Compared to that of the fluctuating fields, we expect the tilt angle from mean fields to be larger if sustained by a large-scale dynamo. We explain this prediction and compare it to the data here.

In shearing box simulations, mean fields can be computed as horizontally averaged fields left as a function of $z$ . The large-scale dynamo would, via the electromotive force, sustain the radial ( $x$ ) field whilst the azimuthal ( $y$ ) field would be stretched by shear. Although fluctuating fields can be stretched coherently in a correlation time, the mean field could be sheared coherently over the longer vertical diffusion time. For mean fields we therefore have $\langle B_{y}\rangle /\langle B_{x}\rangle \sim q{\it\Omega}t_{df}$ where $t_{df}=({\it\alpha}_{ss}c_{s}Hk_{z}^{2})^{-1}$ is the vertical diffusion time of the mean field mode of wavenumber $k_{z}$ . For $c_{s}\sim {\it\Omega}H$ in hydrostatic equilibrium, we then have

(4.10) $$\begin{eqnarray}\tan {\it\theta}_{M}=\langle B_{x}\rangle /\langle B_{y}\rangle \sim ({\it\alpha}_{ss}/q)(k_{z}H)^{2}\sim {\it\alpha}_{ss}/q,\end{eqnarray}$$

for $k_{z}H\sim 1$ . Empirically, ${\it\alpha}_{ss}\propto q/(2-q)$ (Abramowicz et al. Reference Abramowicz, Brandenburg and Lasota1996; Pessah et al. Reference Pessah, Chan and Psaltis2008; Nauman & Blackman Reference Nauman and Blackman2015) so $\tan {\it\theta}_{M}\simeq 1.5{\it\alpha}_{ss,0}/(2-q)$ , where ${\it\alpha}_{ss,0}\equiv {\it\alpha}_{ss}(q=3/2)$ . Combining the latter expression with that for the fluctuating field above we predict that

(4.11) $$\begin{eqnarray}\frac{\tan {\it\theta}_{M}}{\tan {\it\theta}}=\frac{1.5{\it\alpha}_{ss,0}(1+q{\it\Omega}{\it\tau}_{c})}{(2-q)}\ll 1,\end{eqnarray}$$

for ${\it\alpha}_{ss}\sim 0.01$ . The predicted tilt angle ratio from (4.11) is roughly consistent with that found in simulations: for the run shown in figure 3, the tilt angle for the mean field is $\tan {\it\theta}_{M}=(\langle B_{x}\rangle \langle B_{y}\rangle )/\langle B\rangle ^{2}\sim 0.03$ while the tilt angle from everything but the mean gives $(\langle B_{x}B_{y}\rangle -\langle B_{x}\rangle \langle B_{y}\rangle )/\langle B^{2}\rangle \sim 0.26$ . The ratio of these two is 0.14. The predicted ratio from (4.11) for ${\it\alpha}_{ss,0}=0.01$ and $q=3/2$ is 0.06, so the results are not too far off. The stresses from planar averaged mean fields can be significant at given $z$ even if the box averaged stresses appear to show tilt angles more consistent with those for fluctuations.

For large-scale dynamos operating in galaxies similar questions about the tilt angle of the fluctuating and the mean fields can be considered. However, turbulent diffusion is likely primarily due to supernovae and not to the MRI in this case. Extracting the shear dependence must then include consideration of how the turbulent diffusion does or does not correlate with the shear. This circumstance is distinct from the MRI, but is useful to consider when generally probing the predictions of large-scale dynamos. Recently, Van Eck et al. (Reference Van Eck, Brown, Shukurov and Fletcher2015) found that the shear dependence of tilt angle of the mean field in galaxies did not have a shear dependence that was easily captured by simple models.

5.1. Growth of magnetic field from shear-driven turbulence without Coriolis force

In the absence of a Coriolis force, Keplerian shearing boxes are unstable to the development of flow turbulence, but the earliest simulations of this type did not show an accompanying growth of magnetic energy Hawley, Gammie & Balbus (Reference Hawley, Gammie and Balbus1996). This highlighted the fact that the mere existence of turbulence does not guarantee field growth, but the reason for this absence of growth was not fully determined. Is there something fundamental about the turbulence generated by this kind of shear that prevents field growth at all magnetic Reynolds numbers or is it just that e.g. the effective magnetic Reynolds number of the simulation had not reached the critical value for field growth? The difference between these two explanations is important because conventional wisdom from minimalist statistical methods of small-scale magnetic field growth from turbulent flows at high enough Reynolds number would predict small-scale dynamo action. We have revisited this issue with higher resolution simulations and our preliminary results show that at higher resolution, Keplerian shear without the Coriolis force does in fact grow turbulent magnetic energy, even though lower resolution simulations do not. While determination of the critical magnetic Reynolds number is still a work in progress for us. The result for simulations without explicit resistivity or viscosity is shown in figure 4.

In addition, assessing whether large-scale fields (i.e. fields for which the flux of suitably defined averages remains coherent over spatial or temporal scales that are large compared to those of fluctuations) are formed from self-generated linear shear driven turbulence is important in helping to identify the minimum conditions needed for large-scale dynamos. Previous studies have found that imposed forcing of non-helical turbulence combined with linear shear leads to large-scale field growth Yousef et al. (Reference Yousef, Heinemann, Schekochihin, Kleeorin, Rogachevskii, Iskakov, Cowley and McWilliams2008), possibly best interpreted as a ‘fluctuating’ dynamo ${\it\alpha}$ effect Heinemann, McWilliams & Schekochihin (Reference Heinemann, McWilliams and Schekochihin2011).

Figure 4.

Unstratified simulations using ATHENA showing that linear Keplerian shearing boxes without the Coriolis force are both unstable to generation of turbulence and that at high enough resolution, the magnetic field is amplified. The initial conditions of the simulations were zero net magnetic flux with an initial magnetic field $B_{0}=\sin kz\sqrt{2P_{0}/{\it\beta}}$ and an initial plasma ${\it\beta}=P_{th}/P_{mag}=1600$ . The domain size used in $(x,y,z)$ is $H\times 2{\rm\pi}H\times H$ . Time units are $2{\rm\pi}/{\it\Omega}$ , with ${\it\Omega}=10^{-3}$ . The left panel shows the magnetic energy (dashed) and Maxwell stress (solid) and the right panel shows the kinetic energy (dashed) and the Reynolds stress (solid). In each panel the black lines are for a resolution of 16 zones/ $H$ and blue lines are for $32~\text{zones}~\text{h}^{-1}$ . For both resolutions, the right panel shows that the kinetic energy and Reynolds stress grow. The left panel shows however, that only for the higher resolution case do the magnetic energy and Maxwell stress grow. Previous simulations (Hawley et al. Reference Hawley, Gammie and Balbus1996) did not use high enough resolution to obtain field growth.

5.2. What are the non-local implications of shear periodic boundaries and absence of vertical differential rotation?

Because the shear ( $y$ -velocity) in a shearing box has different signs at the inner and outer shear-periodic radial ( $x$ ) boundary, while all other velocities and magnetic fields are periodic, quantities that depend on an odd power of the shear velocity such as the electric field, Poynting flux and even the magnetic helicity flux (Hubbard & Brandenburg Reference Hubbard and Brandenburg2011; Hubbard et al. Reference Hubbard, McNally, Oishi, Lyra and Mac Low2014) can be non-periodic. For the case with an initially radial mean field for example, an inflow of Poynting flux can non-locally add a comparable amount of energy to that supplied locally by the imposed shear (Hubbard et al. Reference Hubbard, McNally, Oishi, Lyra and Mac Low2014). Even if there is no mean field across the box, fluctuations in velocity or magnetic field that depend linearly on the background shear velocity are not guaranteed to be periodic.

What consequences do such fluxes have? The larger the boxes, the longer the time scales over which local dynamics in the box can be studied without the influence of the boundary term (Regev & Umurhan Reference Regev and Umurhan2008), but the larger the boxes, the more unphysical it is to ignore curvature terms and vertical differential rotation (McNally & Pessah Reference McNally and Pessah2014) when interpreting the results for a real system.

5.3. Are shearing boxes with rotation and very large Reynolds numbers nonlinearly stable?

The question of whether rotating quasi-Keplerian shear flows can develop or sustain turbulence has lingered in both theoretical and laboratory contexts of Taylor–Couette flows (the latter being complicated by Ekman circulation) (Longaretti Reference Longaretti2002; Lesur & Longaretti Reference Lesur and Longaretti2005; Balbus Reference Balbus2011; Ji Reference Ji, Bonanno, de Gouveia Dal Pino and Kosovichev2011; Paoletti & Lathrop Reference Paoletti and Lathrop2011; Schartman et al. Reference Schartman, Ji, Burin and Goodman2012). We can however distinguish two sub-questions: (i) Does the Coriolis force always stabilize an already turbulent flow? (ii) Does the Coriolis force prevent the onset of turbulence in the first place? If the answer to the first question were yes, then the answer to the second question would also likely be yes, but if the answer to the first question is no, then the answer to the second question would still be undetermined. While hydrodynamic modal analytic stability analyses and hydrodynamic shearing box simulations to date that include the Coriolis force have been found to be stable, there are physical reasons why instability may not yet be ruled out.

First, modal analysis is insufficient to prove the absence of nonlinear instability. The most relevant example is Couette flow, which is known to be linearly stable in a modal analysis, but unstable in the laboratory. Moreover, psuedospectral, non-modal methods also point to instability Trefethen et al. (Reference Trefethen, Trefethen, Reddy and Driscoll1993). Second, because a linear shear flow without the Coriolis force is unstable, the question arises as to whether there is a range of scales, even in a rotating flow, where the influence of the Coriolis force on the flow is so weak that the flow appears to exhibit a purely linear shear with respect to nonlinear energy transfer. The time scale over which the Coriolis force acts on all scales is ${\it\Omega}^{-1}$ where ${\it\Omega}$ is the angular speed of the rotating frame, but a nonlinear cascade could proceed on a time scale $l/v(l)$ where $l$ is the eddy scale and $v(l)$ is the velocity on that scale. This is particularly relevant to question (i) above as if nonlinear fluctuations or turbulence are already present, we might expect a dynamic range where the turnover times are much shorter than the rotation time. To capture this regime may require very high Reynolds numbers, but there could be a subset of the dynamic range above the dissipation range where the effect of the Coriolis force is small (Longaretti Reference Longaretti2002).

A suggestive feature of the simulations in this context is that for the MRI unstable case, figure 2 shows that the turbulence incurs a nonlinear cascade which is more anisotropic on large scales than small scales. As alluded to earlier, this suggests a nonlinear cascade proceeding such that decreasing scales have relatively shortened nonlinear cascade time scales with respect to the shear time scale, which is independent of scale. In any case, the possibility that a subset of the dynamic range is less influenced by the Coriolis force cannot be identified from simple modal linear theory as evidenced by (4.3); the condition for stability $K_{A}-T<0$ is independent of scale in such a treatment.

5.4. Vertical stresses

Vertical stresses in shearing boxes are small compared to radial stresses and even negative (Miller & Stone Reference Miller and Stone2000). For our run used to make figure 3 for example, we find the box averaged ratio of vertical Maxwell stress from fluctuations to radial Maxwell stress from fluctuations to be $(\langle v_{y}v_{z}\rangle -\langle b_{y}b_{z}\rangle )/(\langle v_{x}v_{y}\rangle \langle b_{x}b_{y}\rangle )=-0.05$ . For real discs, it is unlikely that vertical stresses are negative. The prevalence of outflows, jets and coronae suggest otherwise. The absence of vertical stratification in shearing boxes is a limitation (McNally & Pessah Reference McNally and Pessah2014) and global simulations will be an essential part of the enterprise to obtain realistic estimates of vertical stresses to inform mean field models.

6.1. Including vertical and non-local transport in mean field disc models

As we have seen, the closure that produces the standard SS73 mean field model does not include non-local (large-scale) radial stresses or any vertical stresses. Kuncic & Bicknell (Reference Kuncic and Bicknell2004, Reference Kuncic and Bicknell2007) include mass loss and vertical angular momentum transport to an outflow using a mean field theory. However, they do not separately include vertical stresses from large and small scales as they do not include large-scale magnetic fields. In computing the effect of the vertical transport on the observed radiation spectrum, they compute the reduction in disc emission rather than the specific spectral contributions of corona, jet or wind.

Ultimately, we desire a dynamical theory that governs the relative contributions to transport from small-scale and large-scale magnetic fields and predicts the resultant spectral contributions of disc, corona and jet spectra. If the magnetic field is amplified in situ within the disc, then only structures large enough to overcome turbulent shredding Blackman & Pessah (Reference Blackman and Pessah2009) give rise to coronae. The fraction can be determined if the magnetic energy and stress spectra in the disc can be determined. This provides a guiding principle towards quantifying the fraction of magnetic energy and stress that dissipate in coronae. Some relaxation of this coronal field to larger scales can provide global fields suitable for jets. Studying the processes by which large scale fields in accretion discs grow and evolve must be a component of this effort.

6.2. Role of large-scale dynamos in MRI simulations and transport

In § 4 we discussed that transport emerging from the MRI is non-viscous and that large-scale fluctuations and mean fields make a significant contribution to the stress. Large-scale mean fields emerge in MRI unstable simulations with large enough vertical domains, whether stratified, unstratified, local or global. In shearing boxes or global disc simulations, vertically-dependent mean fields are evident after averaging over radius ( $x$ ) and azimuth ( $y$ ), and commonly exhibit cycle periods of ${\gtrsim}10$ orbits (Brandenburg et al. Reference Brandenburg, Nordlund, Stein and Torkelsson1995; Davis et al. Reference Davis, Stone and Pessah2010; Lesur & Ogilvie Reference Lesur and Ogilvie2010; Guan & Gammie Reference Guan and Gammie2011; Simon, Hawley & Beckwith Reference Simon, Hawley and Beckwith2011; Sorathia et al. Reference Sorathia, Reynolds, Stone and Beckwith2012; Suzuki & Inutsuka Reference Suzuki and Inutsuka2014). In cylindrical global simulations with conducting walls, radially varying large-scale fields are identified after vertically and azimuthally averaging (Ebrahimi & Bhattacharjee Reference Ebrahimi and Bhattacharjee2014). Understanding the mechanism of the underlying large-scale dynamo, the extent to which the stress spectra are converged and the ratio of stress from mean fields to that from the total all need further work.

Whenever a large-scale dynamo is present, a key unifying property is the presence of an electromotive force (per unit charge) $\overline{\boldsymbol{{\mathcal{E}}}}=\langle \boldsymbol{v}\times \boldsymbol{b}\rangle$ . The time evolution of the large-scale field depends on the curl of the EMF. In mean field formalism, solving for the time evolution of the large-scale field requires identifying a turbulent closure for the EMF in terms of the mean fields in the problem (namely, the mean magnetic field and mean velocities) analogous to the need to identify an appropriate closure for the stresses discussed in § 2. For weakly compressible flows, expansion of the EMF in terms of the mean field and its derivatives reveals that large-scale field growth is typically facilitated by a mean magnetic field aligned EMF such that $\overline{\boldsymbol{{\mathcal{E}}}}\boldsymbol{\cdot }\overline{\boldsymbol{B}}\neq 0$ . The importance of $\overline{\boldsymbol{{\mathcal{E}}}}\boldsymbol{\cdot }\overline{\boldsymbol{B}}$ also highlights the utility of tracking the temporal and spatial evolution of magnetic helicity (a topological measure of magnetic field line linkage) because $\langle \overline{\boldsymbol{{\mathcal{E}}}}\boldsymbol{\cdot }\overline{\boldsymbol{B}}\rangle$ can be written as the sum of a time derivative of the magnetic helicity density associated with the large-scale magnetic field, a spatial divergence of large-scale magnetic helicity flux and a resistive term associated with the large-scale current helicity. This same $\overline{\boldsymbol{{\mathcal{E}}}}\boldsymbol{\cdot }\overline{\boldsymbol{B}}$ can also be written as the sum of a time derivative of mean small-scale helicity density plus a spatial divergence of mean small-scale magnetic helicity flux and a resistive term associated with the mean small-scale current helicity (see Brandenburg & Subramanian Reference Brandenburg and Subramanian2005; Blackman Reference Blackman2015 for reviews).

The connection between magnetic helicity and large-scale astrophysical dynamos was first evident in the spectral model of helical MHD turbulence of Pouquet, Frisch & Leorat (Reference Pouquet, Frisch and Leorat1976). They demonstrated an inverse transfer growth of large-scale magnetic helicity for which the driver is the difference between kinetic and current helicites. Kleeorin & Ruzmaikin (Reference Kleeorin and Ruzmaikin1981) presented an equation that couples the small-scale magnetic helicity to the mean electromotive force, but the time evolution was not studied. The spectral work of Pouquet et al. (Reference Pouquet, Frisch and Leorat1976) was re-derived as a time-dependent mean field theory without shear Blackman & Field (Reference Blackman and Field2002) and with shear Blackman & Brandenburg (Reference Blackman and Brandenburg2002)

The most important properties of magnetic helicity underlying its role in dynamo theory are: (i) for a given mircrophysical resistivity, the total magnetic helicity is better conserved than magnetic energy for typical MHD turbulent spectra (Blackman Reference Blackman2004); (ii) for a given magnetic helicity, the energy of the field configuration is minimized when the magnetic field relaxes to the largest scale available subject to the boundary conditions. This can be approximately understood as follows: for a given wavenumber, magnetic energy in a helical field is ${\sim}kH_{M}$ where $H_{M}$ is the magnetic helicity, and so for fixed $H_{M}$ , magnetic energy would decrease as $k$ decreases; (iii) magnetic helicity can change rapidly on large scales and on small scales with opposite signs even when their sum total evolves only slowly on resistive time scales; and (iv) magnetic helicity can flow through boundaries and can be exchanged spatially between sectors within a closed system.

During mean field growth, the field-aligned EMF is substantial but for a turbulent system in a saturated steady state, the growth terms and turbulent decay terms may nearly cancel to a residual value that exactly compensates the microphysical resistive terms. When divergence terms are allowed (either by the global boundary conditions, or because the chosen averaging procedure involves retention of a subset of local spatial dependences), the EMF can be sustained by helicity fluxes. The importance of helicity fluxes to eject small-scale magnetic helicity to maintain fast cycle periods and fast growth rates and avoid catastrophic quenching in stellar and galactic systems has been emphasized (Blackman & Field Reference Blackman and Field2000a ,Reference Blackman and Field b ; Vishniac & Cho Reference Vishniac and Cho2001; Blackman & Brandenburg Reference Blackman and Brandenburg2003; Brandenburg & Subramanian Reference Brandenburg and Subramanian2005; Shukurov et al. Reference Shukurov, Sokoloff, Subramanian and Brandenburg2006; Vishniac Reference Vishniac2009; Ebrahimi & Bhattacharjee Reference Ebrahimi and Bhattacharjee2014). Recent observations of the sun (Zhang et al. Reference Zhang, Moss, Kleeorin, Kuzanyan, Rogachevskii, Sokoloff, Gao and Xu2012; Pipin & Pevtsov Reference Pipin and Pevtsov2014) suggest that both signs of helicity emerge in both hemispheres, a circumstance that helps resolve earlier concerns about which sign dominates Seehafer (Reference Seehafer1990) and how this should be interpreted. When compared with simulations of large-scale dynamos that use the same gauge as used to extract the vector potential from the observations (Pipin & Pevtsov Reference Pipin and Pevtsov2014), the scale-dependent signs of helicity measured from observations seem to agree with expected predictions (Blackman & Brandenburg Reference Blackman and Brandenburg2003). That helicity flux can sustain or drive the growth of large-scale fields in magnetically dominated fusion plasma relaxation or ‘laboratory plasma dynamos’ has long been studied (Strauss Reference Strauss1985; Bhattacharjee & Hameiri Reference Bhattacharjee and Hameiri1986) and there are basic conceptual commonalities between laboratory and astrophysical dynamos despite the different plasma conditions.

The method of averaging can influence the specific form of terms contributing to $\overline{\boldsymbol{{\mathcal{E}}}}$ and its interpretation. The same system averaged in different ways might lead to different explanations of the differently computed mean magnetic fields (Blackman Reference Blackman2015). Mean fields in shearing boxes are typically computed as $x,y$ planar averages leaving quantities as a function of $z$ (see § 4), whereas mean fields in MRI simulations in global cylindrical systems with conducting boundaries have invoked azimuthal and vertical ( ${\it\phi},z$ ) averages (Ebrahimi & Bhattacharjee Reference Ebrahimi and Bhattacharjee2014).

Some have modelled the large-scale cycle periods in shearing box simulations using traditional ${\it\alpha}{-}{\it\Omega}$ type dynamos, with parameterized coefficients (Brandenburg & Donner Reference Brandenburg and Donner1997; Brandenburg Reference Brandenburg, Abramowicz, Björnsson and Pringle1998; Simon et al. Reference Simon, Hawley and Beckwith2011). Taking a step beyond this approach, Gressel (Reference Gressel2010) identified that when the large-scale field evolution is modelled by an ${\it\alpha}{-}{\it\Omega}$ dynamo, the current helicity associated with fluctuations – as previously calculated for stratified sheared rotators Rüdiger & Kichatinov (Reference Rüdiger and Kichatinov1993) – seems to correlate well with the time averaged behaviour of ${\it\alpha}$ empirically determined from the simulations. In addition, Gressel (Reference Gressel2010) found that the time-dependent dynamical quenching predicted using an ${\it\alpha}$ effect represented by the difference between kinetic and current helicities is promising to explain the time-dependent large-scale dynamo. Others have begun to investigate in more detail the role and form of the EMF in MRI unstable flows (Vishniac Reference Vishniac2009; Käpylä & Korpi Reference Käpylä and Korpi2011; Ebrahimi & Bhattacharjee Reference Ebrahimi and Bhattacharjee2014; Ebrahimi & Blackman Reference Ebrahimi and Blackman2015; Squire & Bhattacharjee Reference Squire and Bhattacharjee2015a ,Reference Squire and Bhattacharjee b ). There is much to be done on this front, and to connect with the basic question of what minimal ingredients are needed for large-scale dynamos, as alluded to above in § 5.1. Different circumstances, even if not minimal, likely lead to different important contributions to the EMF.

The presence and importance of large-scale dynamos in accretion simulations is not entirely surprising because accretion disc theory and large-scale dynamo theory are actually artificially separated components of what really should be captured by a unified mean field theory. Accretion disc theory has traditionally ignored the evolution of the large-scale fields whilst mean field dynamo theory has ignored the dynamics of transport. Ultimately, the two must operate together. There have been some steps along this path of unification (Rüdiger et al. Reference Rüdiger, Elstner and Schultz1993; Campbell & Caunt Reference Campbell and Caunt1999; Campbell Reference Campbell2000; Moss, Shukurov & Sokoloff Reference Moss, Shukurov and Sokoloff2000; Rekowski, Rüdiger & Elstner Reference Rekowski, Rüdiger and Elstner2000; Moss & Shukurov Reference Moss and Shukurov2004). These papers, along with Kuncic & Bicknell (Reference Kuncic and Bicknell2004), provide important starting points but have not incorporated modern lessons learned from the MRI or magnetic helicity dynamics. Efforts to include all of these ingredients are beginning e.g. Vishniac (Reference Vishniac2009), Ebrahimi & Bhattacharjee (Reference Ebrahimi and Bhattacharjee2014) but we are far from the ‘practical model’ stage in predicting observed spectra and temporal evolution. Again, there is much opportunity for further work.

6.3. Global simulations and a unified mean field accretion/dynamo theory

Evidence for significant non-local transport by large-scale fields is also evident in global simulations De Villiers & Hawley (Reference De Villiers and Hawley2003), Fromang & Nelson (Reference Fromang and Nelson2006), Beckwith et al. (Reference Beckwith, Armitage and Simon2011), Sorathia et al. (Reference Sorathia, Reynolds, Stone and Beckwith2012), Penna et al. (Reference Penna, Narayan and Sa̧dowski2013), Suzuki & Inutsuka (Reference Suzuki and Inutsuka2014). How the transport in such magnetized high magnetic Reynolds number discs is proportioned between disc, jets, outflows and coronae and what determines this proportion is not yet understood.

A complicating feature is the interaction between disc and stellar or blackhole magnetosphere (Ghosh & Lamb Reference Ghosh and Lamb1978; Küker, Henning & Rüdiger Reference Küker, Henning and Rüdiger2003; Matt & Pudritz Reference Matt and Pudritz2005; Perna et al. Reference Perna, Bozzo and Stella2006; Penna, Sa̧owski & McKinney Reference Penna, Sa̧owski and McKinney2012; Romanova et al. Reference Romanova, Ustyugova, Koldoba and Lovelace2012). If the spin of the central object such as a neutron star or black hole acts as a significant source of input energy, then the jet or coronae need not be entirely sourced by the accretion. The interaction undoubtedly influences jet formation, disc truncation and the coronal energy supply, depending on the compactness, spin and magnetic field of the central object.

There is evidence from numerical simulations Penna et al. (Reference Penna, Narayan and Sa̧dowski2013) that rapidly spinning black hole engines launch jets consistent with the mechanism of Blandford & Payne (Reference Blandford and Payne1982) being dominant. And, as long as the simulations run long enough, emergence of the central jet appears to be independent of whether the initial field configuration consists of self-contained loops in the disc or more open initial large-scale poloidal fields. The processes by which the field in the disc seeds the MRI, evolves and amplifies in the central regions (via some combination of advection and dynamo action), and subsequently relaxes to large scales to form jets all suggest that tracking magnetic helicity evolution will provide useful insight. Conceptually, a generic set of processes seems to be occurring; field amplification in a flow-dominated disc produces structures with large enough scale to buoyantly rise into the corona and avoid turbulent shredding. Once in the corona, some of these structures relax to form very large fields in the low plasma ${\it\beta}$ regime corona. Generalizations of the simulations of Penna et al. (Reference Penna, Narayan and Sa̧dowski2013) to include a wider range of initial configurations with even smaller scale initial fields would be of interest.

If global simulations ultimately reveal that a diversity of initial states lead to a much less diverse set of end states, that is a helpful guide for mean field theories that couple accretion, field growth, coronae and jets. Some use of mean field theory as an input to global simulations is emerging Stepanovs, Fendt & Sheikhnezami (Reference Stepanovs, Fendt and Sheikhnezami2014), Sa̧dowski et al. (Reference Sa̧dowski, Narayan, Tchekhovskoy, Abarca, Zhu and McKinney2015) but also important is the other direction – informing mean field models with outputs from global simulations. For most of the long history of simulations of jets (Pudritz et al. Reference Pudritz, Hardcastle and Gabuzda2012), the magnetic fields have been imposed and the disc treated as a boundary condition but the new generation of simulations is evolving beyond this limitation. A similar evolution between connecting interior to exterior is emerging in solar dynamo simulations as well (Nelson et al. Reference Nelson, Brown, Brun, Miesch and Toomre2011; Nelson & Miesch Reference Nelson and Miesch2014).

6.4. Quantifying the precision of mean field models

In § 2 we explained that standard axisymmetric accretion theory makes sense only as a mean field theory and thus has a finite precision. In fact, both mean field dynamo and accretion theory are subject to the need to quantify this finite precision. There has been only limited work toward this goal Blackman (Reference Blackman1998), Blackman et al. (Reference Blackman, Nauman and Edgar2010). The precision of the theory itself depends on the size and temporal scales of the fluctuations over which averaging must take place to construct the mean field theory. The larger these fluctuations, the less precise the theory. There is a wide range of ratios of observational integration times versus dynamical times in accreting systems, but characterization of this imprecision is particularly important when the integration times are very short such as for YSO infrared observations. Apparent disagreements between observations and theory must be evaluated with theoretical imprecision kept in mind. By analogy, offsets between the emergent spin axis and magnetic axis in stars or planets with dynamos need not require a separate systematic physical explanation if such deviations are consistent simply with the predicted statistical imprecision in the mean field dynamo theory itself.

6.5. Developing an analogue to the HR diagram for accretion discs?

If we are to ultimately understand disc evolution to a level that matches stellar evolution, we will need a characterization of observations analogous to that of the Hertzsprung–Russell diagram for stars. This should provide a systematic way to document the evolution of discs through various states of luminosity, coronal and outflow activity. For accretion discs, the effective surface temperature varies outward across the surface of the disc, although the emission from the inner-most regions likely dominates the luminosity. Among the features to capture in a diagram might be how far the thermal component of the disc extends inward, a measure of the fraction of non-thermal versus thermal emissions (spectral hardness), a measure of the jet contribution and a characterization of variability. For a given type of accretion engine, the evolution through these states may depend primarily on the accretion rate. For microquasars, the long-observed evolution between states are being characterized by such diagrams Zhang (Reference Zhang2013), Kylafis & Belloni (Reference Kylafis, Belloni, Contopoulos, Gabuzda and Kylafis2015). Generalizing such diagrams to other types of accretion engines is important.

Black hole and stellar accreters have a range of masses whereas white dwarf and neutron star masses do not vary much. For rapidly spinning neutron stars, black holes and classical novae, there are additional energy sources beyond that of the accretion itself. All such properties should be ultimately evident in the ensemble of evolutionary paths on appropriate state diagrams that provide the HR diagram analogue.