Resolution Dependence of Magnetorotational Turbulence in the Isothermal Stratified Shearing Box

TL;DR

This study investigates how the resolution affects the properties of magnetorotational turbulence in stratified shearing boxes, revealing that shear stress and magnetic correlation length depend on resolution and do not converge.

Contribution

It provides a detailed resolution study of MRI-driven turbulence in stratified shearing boxes, highlighting the resolution dependence of key turbulence metrics.

Findings

Shear stress scales as N^{-1/3} with resolution.

Magnetic correlation length decreases as N^{-1/2}.

Results are consistent with previous studies but show lack of convergence.

Abstract

Magnetohydrodynamic (MHD) turbulence driven by the magnetorotational instability can provide diffusive transport of angular momentum in astrophysical disks, and a widely studied computational model for this process is the ideal, stratified, isothermal shearing box. Here we report results of a convergence study of such boxes up to a resolution of $N = 256$ zones per scale height, performed on blue waters at NCSA with ramses-gpu. We find that the time and vertically integrated dimensionless shear stress $\overline{\alpha} \sim N^{-1/3}$, i.e. the shear stress is resolution dependent. We also find that the magnetic field correlation length decreases with resolution, $\lambda \sim N^{-1/2}$. This variation is strongest at the disk midplane. We show that our measurements of $\alpha$ are consistent with earlier studies. We discuss possible reasons for the lack of convergence.

Figures

Tables

Table 1. Table 1: Convergence Properties of MHD Disk Turbulence Models

Geometry	Stratified	Net	Dissipation	Isothermal	Convergent	Maximum	References
		Flux				Resolution
local	no	zero	ILES	yes	no	$256/H$	(1) (2) (3) (4)
local	no	zero	explicit	yes	yes	$512/H$	(5)
local	no	mean	ILES	yes	yes	$256/H$	(6) (2) (3)
local	no	mean	explicit	yes	yes	$800/H$	(7)
local	no	mean	ILES	no	yes	$64/H$	(8)
local	yes	zero	ILES	yes	this work	$256/H$	(9) (10) (11)
local	yes	zero	ILES	no	unclear	$64/H$	(12) (13)
local	yes	zero	explicit	yes	unclear	$128/H$	(14) (15)
local	yes	mean	ILES	yes	unclear	$48/H$	(16) (17)
local	yes	zero	ILES	no	unclear	$64/H$	(18) (19)
global	no	zero	ILES	yes	yes	$480\times 1920\times 128$	(20)
global	yes	zero	ILES	no	unclear	$768\times 256\times 256$	(21) (22) (23)
global	yes	mean	ILES	no	unclear	$288\times 128\times 128$	(24) (25) (26)

Table 2. Table 2: Model parameters

Label	N = Zones/$H$	$t_0\left({\Omega }^{-1}\right)$	$\Delta t\left({\Omega }^{-1}\right)$
r32	32	1800	300
r64	64	2100	300
r128	128	2400	300
r256	256	2648	288

Table 3. Table 3: Model results

Label	$\overline{\alpha }$	$\overline{⟨-B^x{B}^y⟩}$	$\overline{⟨\rho v^x\delta v^y⟩}$	$\Delta M/M_0$	${\sigma }_{\alpha }/\overline{\alpha }$	${\lambda }_{\mathrm{minor}}$	${\lambda }_{\mathrm{major}}$	${\theta }_{\mathrm{tilt}}$
r32	0.039	0.0061	0.0017	0.69%	0.24	0.12	0.61	16.0°
r64	0.034	0.0053	0.0015	0.65%	0.37	0.085	0.40	17.8°
r128	0.025	0.0039	0.0011	0.55%	0.26	0.060	0.27	18.6°
r256	0.019	0.0029	0.0008	0.40%	0.23	0.043	0.20	19.0°

Table 4. Table 4: RMS and standard deviation of net magnetic fluxes present for each run. Vertical magnetic flux is zero, conserved to machine precision.

Label	${\overline{⟨B^x⟩}}_{\mathrm{RMS}}$	${\sigma }_{⟨B^x⟩}$	${\overline{⟨B^y⟩}}_{\mathrm{RMS}}$	${\sigma }_{⟨B^y⟩}$	${\overline{\alpha }}_{\mathrm{NF}}$
r32	$7.0\times {10}^{-4}$	$4.2\times {10}^{-4}$	$3.6\times {10}^{-2}$	$1.5\times {10}^{-2}$	$1.1\times {10}^{-2}$
r64	$4.0\times {10}^{-4}$	$3.8\times {10}^{-4}$	$1.8\times {10}^{-2}$	$1.6\times {10}^{-2}$	$5.6\times {10}^{-3}$
r128	$4.3\times {10}^{-4}$	$3.2\times {10}^{-4}$	$2.0\times {10}^{-2}$	$1.4\times {10}^{-2}$	$6.2\times {10}^{-3}$
r256	$3.7\times {10}^{-4}$	$3.5\times {10}^{-4}$	$1.1\times {10}^{-2}$	$1.0\times {10}^{-2}$	$3.4\times {10}^{-3}$

Equations

$$Re\equiv \frac{c_{s}H}{\nu }R{e}_M\equiv \frac{c_{s}H}{\eta }P{r}_M\equiv \frac{\nu }{\eta }.$$

$$\left(x,y,z\right)=\left(R-R_0,R_0\left(\varphi -\Omega t-{\varphi }_0\right),z\right).$$

$$\frac{\partial \rho }{\partial t}$$

$$\frac{\partial \rho \mathbf{v}}{\partial t}$$

$$-2\rho \Omega {\hat{\mathbf{e}}}_z\times \mathbf{v}+\rho \nabla \left(-\frac{3}{2}{\Omega }^2{x}^2+\frac{1}{2}{\Omega }^2{z}^2\right),$$

$$\frac{\partial \mathbf{B}}{\partial t}$$

$$\nabla \cdot \mathbf{B}=0.$$

$$\rho$$

$$\mathbf{v}$$

$$c_s=1\Omega =1{\rho }_0=1$$

$$N=\frac{H}{\Delta x}$$

$$⟨f⟩\equiv \frac{\int dxdydzf}{\int dxdydz},$$

$$[f]\equiv \frac{\int dxdyf}{\int dxdy},$$

$$\overline{f}\equiv \frac{\int dtf}{\int dt}.$$

$$\alpha \equiv \frac{⟨\rho v^x\delta v^y-B^x{B}^y⟩}{⟨P⟩}.$$

$$T^{ij}(\delta x,\delta y,z,t)\equiv \frac{1}{[B^2]}[\delta B^i(x,y,z,t)\delta B^j(x+\delta x,y+\delta y,z,t)].$$

$$\alpha (t)={\alpha }_0+\sum _j{a}_j\cos ({\omega }_{j}t+{\varphi }_j);{\omega }_j=\frac{2\pi j}{T}$$

$$P(a_j)d{a}_j=\exp \left(-\frac{a_j^2}{2⟨a_j^2⟩}\right)\frac{d{a}_j}{(2\pi ⟨a_j^2⟩{)}^{1/2}}$$

$$⟨{\sigma }^2⟩=⟨(\alpha -{\alpha }_0{)}^2⟩=\frac{1}{2}{\int }_0^{\infty }\frac{d\omega }{2\pi }{P}_{\omega }.$$

$$\xi (\tau )\equiv \frac{1}{T}{\int }_0^{T}dt\alpha (t)\alpha (t+\tau )-{\alpha }_0^2=\frac{1}{2}{\int }_0^{\infty }\frac{d\omega }{2\pi }{P}_{\omega }\cos \omega \tau .$$

$$m=\frac{1}{\Delta t}{\int }_0^{\Delta t}dt(\alpha -{\alpha }_0)$$

$$m=\frac{1}{\Delta t}{\int }_0^{\Delta t}dt\sum _j{a}_j\cos ({\omega }_{j}t+{\varphi }_j)=\sum _j{a}_j\frac{1}{{\omega }_j\Delta t}\left\{\sin ({\omega }_j\Delta t+{\varphi }_j)-\sin {\varphi }_j\right\}$$

$$⟨m^2⟩=\sum _j⟨a_j^2⟩\frac{1}{{\omega }_j^2\Delta t^2}\left(1-\cos ({\omega }_j\Delta t)\right)\to \int \frac{d\omega }{2\pi }\frac{P_{\omega }}{{\omega }^2\Delta t^2}\left(1-\cos (\omega \Delta t)\right).$$

$$P_{\omega }\propto (1+(\omega /{\omega }_0{)}^2{)}^{-p/2}.$$

$$P_{\omega }=\frac{\Gamma (p/2)8\sqrt{\pi }}{\Gamma ((p-1)/2)}\frac{{\sigma }^2}{{\omega }_0}\frac{1}{(1+{\omega }^2/{\omega }_0^2{)}^{p/2}}.$$

$$\xi (\tau )=\frac{2^{3/2-p/2}}{\Gamma ((p-1)/2)}{\sigma }^2({\omega }_0|\tau |{)}^{(p-1)/2}{K}_{(p-1)/2}({\omega }_0|\tau |)$$

$$⟨(\alpha -\overline{\alpha }{)}^2⟩={\int }_0^{\infty }\frac{d\omega }{2\pi }{P}_{\omega }\left(\frac{1}{2}+\frac{\cos (\omega \Delta t)-1}{{\omega }^2\Delta t^2}\right).$$

$$⟨m^2⟩=\frac{{\sigma }^2}{{\omega }_0\Delta t}2\left(1-\frac{1-e^{-{\omega }_0\Delta t}}{{\omega }_0\Delta t}\right).$$

$$⟨(\alpha -\overline{\alpha }{)}^2⟩={\sigma }^2\left(1-\frac{2}{{\omega }_0\Delta t}+\frac{2(1-e^{-{\omega }_0\Delta t})}{{\omega }_0^2\Delta t^2}\right).$$

Full text

\WarningFilter

revtex4-1Repair the float

1 Introduction

Astrophysical disks form in galaxies and around black holes, neutron stars, white dwarfs, main sequence stars, and planets because the angular momentum of the parent plasma is approximately conserved while kinetic energy in noncircular or non-coplanar motion is easily dissipated and radiated away. Disk evolution is therefore governed by angular momentum transport, which can take the form of external torques (e.g. from magnetized winds) or internal transport.

Diffusive internal transport of angular momentum has been fruitfully described with the phenomenological anomalous viscosity, or $\alpha$, model (Shakura & Sunyaev 1973; Lynden-Bell & Pringle 1974), which attributes transport to localized turbulence. No general driver of turbulence in non-self-gravitating Keplerian disks was known until the discovery by Balbus & Hawley (1991) of the magnetorotational instability (MRI), a local, linear instability of weakly magnetized disks. Subsequent nonlinear numerical studies convincingly demonstrated that the MRI leads to turbulence and outward angular momentum transport (see the review of Balbus & Hawley, 1998). Later work has uncovered purely hydrodynamical instabilities including the zombie vortex (Marcus et al. 2015, but see Lesur & Latter 2016), the vertical shear instability (Urpin 2003, Nelson et al. 2013), the baroclinic instability (Klahr & Bodenheimer, 2003; Petersen et al., 2007a, b; Lesur & Papaloizou, 2010), and convective overstability (Klahr & Hubbard 2014). Nonetheless, MRI-driven turbulence remains the leading candidate for driving disk evolution in many astrophysical settings.

Our paper probes numerical convergence of MHD turbulence in a particular disk model. By convergence, we mean resolution and dissipation-scale independence in average quantities like the angular momentum flux. We begin by reviewing the various classes of numerical models used to study MHD turbulence in disks, and describing the claims of convergence or nonconvergence made for each class.

Numerical simulations of disk turbulence can be divided into local and global models. In a local model (or shearing box; Goldreich & Lynden-Bell 1965, Hawley et al. 1995), the equations of motion are expanded to lowest order in the ratio of the scale height $H$ to the local radius $r$ in a co-orbiting Keplerian frame. Differential rotation manifests as a linear shear flow. The shearing box boundary conditions then make it possible to model the disk in a shear-periodic, rectangular box. The local model is highly symmetric and cannot, for example, distinguish between the inward and outward directions (it is symmetric under a rotation by $\pi$ around the $z$ axis). In a global model, by contrast, one simulates some radial range within a disk without requiring $H\ll r$. Global models do not have the inward-outward symmetry of the local model.

The vertical ($z$) component of gravitational acceleration in the local model is $-\Omega^{2}z$, where $\Omega\equiv$ orbital frequency. Unstratified local models turn off the vertical component of gravity, begin with a uniform vertical density profile, and typically use periodic vertical boundary conditions. Stratified local models turn on the vertical component of gravity, begin with a $z$-dependent vertical density profile, and use a variety of vertical boundary conditions.

For most boundary conditions, local simulations conserve one or more components of the mean magnetic field. For example, unstratified local models with periodic vertical boundary conditions conserve the mean vertical and toroidal field if the mean radial field vanishes.111A nonvanishing mean radial field is conserved, but it causes the toroidal field to vary linearly in time. See Hawley et al. (1995). Numerical investigations show that the mean field can have a profound effect on the saturated turbulent state, so we need to distinguish between zero mean field models, where all the currents that sustain the field are contained within the simulation volume and can therefore decay, and mean field models, where one or more components of the field is fixed by the boundary conditions and cannot decay.

Turbulence leads to dissipation. In explicit dissipation models (or direct numerical simulations) dissipation is incorporated directly in the model, for example by a scalar viscosity $\nu$ and resistivity $\eta$ that are dimensionlessly parameterized by the Reynolds numbers and their ratio, the magnetic Prandtl number:

[TABLE]

In implicit large eddy simulation (ILES) models there is no explicit dissipation, and dissipation is provided by the numerical scheme through truncation error at the grid scale. Notice that for ILES models run with a conservative scheme, lost kinetic and magnetic energy is entirely captured as plasma thermal energy. In this sense reconnection can be “included” in an ILES model, although the reconnection rate and dynamics may be incorrect.

The consequences of using ILES to study high Reynolds number hydrodynamic turbulence are fairly well understood (e.g. Sagaut 2006): if there is sufficient dynamic range (large enough zone number) then the character of dissipation at small scales has little influence on turbulent structures at large scales. It is large-scale structures that often determine the flow properties of greatest astrophysical interest, such as turbulent momentum flux. The consequences of using ILES to study high Reynolds numbers magnetohydrodynamic (MHD) turbulence are less well understood (Miesch et al. 2015). It is fair to say that many disk simulators (including us) have frequently assumed that with enough dynamic range MHD ILES would converge to the astrophysically relevant high Reynolds numbers result (but see Lesur & Longaretti 2007, Fromang et al. 2007, Longaretti & Lesur 2010, Meheut et al. 2015, Walker et al. 2016).

Finally, disk simulations can be subdivided according to their treatment of heating and cooling of the plasma. Direct simulation of the interaction of radiation with matter has, until recently, been expensive in comparison to available computational resources. Most disk simulations have therefore used simplified treatments of plasma thermodynamics, with phenomenological cooling and heating, or assumed an isothermal equation of state with pressure $P=\rho c_{s}^{2}$, and $c_{s}$ constant in time and space. Isothermal models are relevant to disks heated by external illumination, such as disks around compact objects at many gravitational radii, where the thermal timescale can be short compared to the dynamical timescale.

Local models also depend on the box dimensions which are purely numerical parameters. Changes in box sizes are known to produce qualitative changes in shearing box models (e.g. Simon et al. 2012, Shi et al. 2016). Even the largest domains find correlations on the scale of the box, at least in the corona (Guan & Gammie 2011). Two related questions emerge. Does the shearing box model converge as the box sizes goes to infinity? Does shearing box evolution match global behavior as the box size goes to infinity? These questions are challenging to answer numerically.

Much is now understood about convergence of the gross, time-averaged properties of MRI-driven turbulence (e.g. $\alpha$) in every corner of the five dimensional disk model parameter space: local/global, stratified/unstratified, mean/zero net field, ILES/explicit dissipation, isothermal/nonisothermal. A summary of previous calculations emphasizing convergence is given in Table 1.

Zero net field, local, unstratified, isothermal, ILES models are particularly interesting: Fromang & Papaloizou (2007) showed that these models are nonconvergent (see also Pessah et al. 2007), and this has been independently confirmed (Simon et al. 2009, Guan et al. 2009). With $N$ the number of resolution elements along one axis, with zone aspect ratios fixed, nonconvergence appears as $\alpha\propto N^{-1}$ (but see Bodo et al. 2011) and magnetic correlation length $\lambda\propto N^{-1}$ (i.e. correlation length is proportional to zone size). But this is not the full story: Shi et al. (2016) have recently found convergence if the vertical extent of the model is large compared to the radial extent. In this case MHD turbulence excites waves that travel vertically, and this may be connected to the butterfly oscillations seen in stratified models. However, the connection between these tall boxes and traditional unstratified (and stratified) shearing boxes is still uncertain, and we consider it premature to change the relevant conclusion for convergence in Table 1.

Unstratified models converge, however, if either explicit dissipation (Fromang 2010, but see Bodo et al. 2011) or a mean magnetic field (Simon et al. 2009, Guan et al. 2009) are added. When a mean field is added $\alpha$ increases proportional to the mean field strength (Hawley et al. 1995, Salvesen et al. 2016).

What about stratified models? One might think that stratification would lead to magnetically driven convection which could organize the field on the scale of the convective eddies, leading to convergence. But the numerical evidence for convergence of zero net field, local, stratified, isothermal ILES models is contradictory. The work of Davis et al. (2010), using the athena code, is consistent with convergence, while the work of Bodo et al. (2014), using the pluto code, shows a sharp drop in Maxwell stress at the highest accessible resolution of $200/\sqrt{2}\simeq 141$ zones per scale height. The question of convergence for stratified, isothermal ILES models is particularly pressing because they are sometimes used to interpret observations in both local (e.g. Simon et al. 2015) and global (e.g. Mościbrodzka et al. 2009, Flock et al. 2015) forms.

This paper therefore returns to study the convergence of zero net field, local, stratified, isothermal ILES models at high resolutions made newly accessible by the combination of NCSA’s blue waters machine and the ramses-gpu code. In Section 2 we present the physical model and numerical method. Section 3 contains the results of our calculations. Section 4 discusses the implications of our results and future directions. Section 5 concludes.

2 Model

2.1 Governing Equations

The local model expands the equations of motion to lowest nontrivial order around a Keplerian orbit at $R=R_{0},\phi=\phi_{0}+\Omega t,z=0$ and defines the local Cartesian coordinates

[TABLE]

In the local model for a Keplerian disk the equations of ideal MHD, with an isothermal equation of state ($P=\rho c_{s}^{2}$; $P\equiv$ pressure, $\rho\equiv$ density, $c_{s}\equiv$ sound speed, which is assumed constant), are

[TABLE]

where $\boldsymbol{v}\equiv$ velocity in the local frame and $\boldsymbol{B}\equiv$ magnetic field, subject to the constraint

[TABLE]

Equation (4) includes Coriolis and tidal forces. Notice that there is no explicit dissipation (resistivity or viscosity) and that $R_{0}$ does not appear in the governing equations.

For $\boldsymbol{B}=0$ these equations admit the equilibrium

[TABLE]

Here $H=c_{s}/\Omega$. Notice that others (e.g. Davis et al. 2010, Bodo et al. 2014) define the scale height as $\sqrt{2}c_{s}/\Omega$. This implies that their resolution should be multiplied by $1/\sqrt{2}$ for comparison with ours. The initial conditions for our model are the unmagnetized equilibrium (7)-(8), seeded with a uniform toroidal field at $|z|<2H$; ${\bf B}=0$ elsewhere. The initial plasma $\beta\equiv 2P/B^{2}=50$ at the midplane.

Hereafter we set

[TABLE]

which together imply that $H=1$ and the surface density $\Sigma=\sqrt{2\pi}$. The mass, length, and time units are thus $\rho_{0}H^{3}$, $H$, and $\Omega^{-1}$, respectively. Occasionally we reinsert these units for clarity.

For the $x$ and $y$ boundaries we use shearing box boundary conditions (see Hawley et al., 1995). With these boundary conditions the model is translation-invariant in the $x-y$ plane, and also invariant under rotations by $\pi$ around the $z$ axis. In addition, the vertical magnetic flux $\Phi_{z}\equiv\int dxdyB_{z}$ (integral taken over the entire $x-y$ domain at any $z$) is conserved. Our initial conditions have $\Phi_{z}=0$, and our model extends over $-L_{x}/2<x<L_{x}/2$, with $L_{x}=3$, and over $-L_{y}/2<y<L_{y}/2$, with $L_{y}=4$.

At least three different $z$ boundary conditions have been used for stratified shearing boxes. Beginning with Stone et al. (1996), many have used outflow boundary conditions. Davis et al. (2010) used periodic boundary conditions, which have the advantage that all three components of the mean magnetic field are conserved in exchange for altering the domain topology. Several authors (Brandenburg et al. 1995, Bodo et al. 2014) adopt impenetrable, stress-free boundaries that set $\partial v_{x}/\partial z=\partial v_{y}/\partial z=v_{z}=0$ and $B_{x}=B_{y}=\partial B_{z}/\partial z=0$ (or the equivalent conditions on the magnetic vector potential). The effect of boundary condition choice has not been systematically studied at modern resolution, although Stone et al. (1996) found that an artificial resistive layer at $2<|z|<3$ did not affect midplane dynamics significantly, and Oishi & Mac Low (2011) demonstrate similar behavior in three runs that differ only by choice of vertical boundary conditions. For finite thermal diffusivity, Gressel (2013) find a significant change in energy transport between outflow and impenetrable vertical boundaries.

We chose outflow boundary conditions and a large vertical extent to minimize the influence of the vertical boundaries on the body of the disk. Formally, outflow boundary conditions are $\partial{\bf B}/\partial z=0$ and $\partial{\bf v}/\partial z=0$, and $-(1/\rho)\partial(c_{s}^{2}\rho)/\partial z-\Omega^{2}z=0$, consistent with hydrostatic equilibrium. Our model extends over $-L_{z}/2<z<L_{z}/2$ with $L_{z}=12$.

Outflow boundary conditions imply that fluid mass in the computational domain is not conserved. The characteristic outflow timescale $\tau_{\mathrm{out}}\equiv\Sigma/\dot{\Sigma}$. Assuming sonic outflow at the boundaries, $\dot{\Sigma}\approx 2\rho(|z|=L_{z}/2)c_{s}$. Using the density profile fit from Guan & Gammie (2011) that takes account of magnetic support of the disk atmosphere, $\tau_{\mathrm{out}}\approx 6\times 10^{4}\Omega^{-1}$. This is long compared to our integration times. Outflow boundary conditions also imply that the radial and toroidal magnetic flux are not conserved.

The domain size $(L_{x},L_{y},L_{z})=(3,4,12)H$ may influence the turbulent state. Guan & Gammie (2011) and Simon et al. (2012) provide numerical evidence that angular momentum transport and variability may depend on structures large compared to $H$, but such large domains are currently inaccessible at our target resolution. Simon et al. (2012) demonstrate a transition to anomalous behavior as $L_{x}$ goes from $2$ to $0.5$. The minimum $L_{x}$ that avoids these pathologies is known to be less than $L_{x}=3$; the results of Davis et al. (2010) suggest that this minimum is less than $L_{x}=\sqrt{2}$.

Finally, the integration time $\Delta t$ should be long enough that average values for $\alpha$ and other quantities of interest can be measured with reasonable signal to noise. Our typical integration time is $\approx 300\Omega^{-1}$ (see the Appendix for a discussion of measurement errors).

2.2 Numerical Methods

We integrate the model with ramses-gpu (Fromang et al., 2006; Kestener et al., 2010, 2014), a modern astrophysical MHD code with support for GPU acceleration222Freely available: http://www.maisondelasimulation.fr/projects/RAMSES-GPU/html/. ramses-gpu is a second-order finite volume MUSCL scheme. Fluxes are evaluated with the HLLD approximate Riemann solver (Miyoshi & Kusano 2005). The constraint $\nabla\cdot{\bf B}=0$ is preserved via constrained transport with face-centered magnetic fields (Evans & Hawley 1988).

Numerical resolution is characterized by

[TABLE]

i.e. the number of zones per scale height in the radial direction. We take $\Delta x:\Delta y:\Delta z=1:2:1$, so this is also the number of zones per scale height in the vertical direction, and twice the number of zones per scale height in the azimuthal direction. Hawley et al. (2011) showed that for MRI growth the azimuthal direction is typically better resolved than the vertical direction by a factor of a few in shearing boxes, as did Parkin & Bicknell (2013). Guan et al. (2009) showed that the autocorrelation function of the magnetic field in unstratified, isothermal shearing box models is anisotropic and approximately in the ratio $1:4:1$ in the radial, azimuthal, and vertical directions, suggesting that near the midplane the $y$ direction is slightly better resolved than $x$ and $z$ in our model.

The mean azimuthal velocity $\overline{v}_{y}=-(3/2)\Omega x$. Truncation error depends on the velocity of the fluid with respect to the grid, and therefore if $\overline{v}_{y}$ is the dominant component of the velocity field the truncation error will vary systematically with $x$. This problem can be solved by using orbital advection (also known as “FARGO”; Masset 2000) for the MHD equations (Johnson et al. 2008, and references therein). We do not use orbital advection, but the shear velocity at the edge of our boxes is only $1.5c_{s}$, and we have checked that the Maxwell and Reynolds stresses do not vary significantly with $x$.

We start preliminary models from smooth initial conditions. These were seeded with white noise, with $\delta v_{i}\sim 0.01c_{s}$ and $\delta\rho\sim 0.01\rho_{0}$ to excite a spectrum of unstable modes. We used late-time snapshots from these models to initialize our production models. Each run at resolution $2N$ was initialized with a snapshot from the final (or near-final) state of a model with resolution $N$ using a divergence-free prolongation operator (Fromang et al., 2006). While this avoids running high resolution models through an initial transient phase (and allows our model to forget the initial net azimuthal magnetic flux), it does introduce a potential bias by coupling the initial state of one simulation to the final (or near-final) state of a lower resolution model.

Stratified shearing box models have high Alfvén speeds in the upper atmosphere ($v_{A}\sim\rho^{-1/2}$), which via the Courant condition can demand an impractically small timestep. This is a standard problem in numerical MHD, and can be solved by applying a density floor, or re-introducing a displacement current that limits the Alfvén speed to a maximum speed (Boris, 1970). In shearing box models, Miller & Stone (2000) used a version of the Boris fix with speed of light $v_{A,max}=(1,4,8)c_{s}$. Guan & Gammie (2011), by contrast, impose a density floor of $10^{-5}\rho_{0}$. We impose a density floor such that $v_{A}<v_{A,max}=10c_{s}$. Our $v_{A,max}$ is higher than the expected $v_{A}$ at $z=6$ (as deduced from the fit to averaged stratified shearing box properties of Guan & Gammie 2011) but small enough to limit the integration to a practical timestep.

In characterizing the saturated state, we use the following averages: an average over volume

[TABLE]

an average over $x$ and $y$

[TABLE]

and an average over time

[TABLE]

The height-integrated Shakura-Sunyaev $\alpha$ parameter is

[TABLE]

This definition does not depend explicitly on box size. It is the average used in height-integrated disk evolution models (e.g. King et al. 2007) for comparison with observation.

3 Results

We consider four models, marching forward in linear resolution by factors of two from $N=32$ to $N=256$. Each model is started using late-time data from the preceding lower-resolution model. All share a common coordinate time $t$. The runs and their linear resolution, initial time $t_{0}$, and duration $\Delta t$ are given in Table 2. We define $t^{\prime}$ for each run as $t-t_{0}$. Poloidal slices from all four resolutions are shown in Figure 1.

For $L_{z}=12H$ about $0.5\%$ of the disk mass is lost per $300\Omega^{-1}$ after accounting for mass added via the density floor (see Table 3; $M_{0}$ is the mass of the disk at the start of that run).

We now turn to the effects of resolution on one- and two-point statistics of the saturated state. Section 3.1 considers volume- and area-averaged quantities over the domain. Section 3.2 presents the correlation function of the magnetic field.

3.1 Space and Time Averages

Does $\overline{\alpha}$ depend on resolution? Figure 2 shows $\alpha$ as a function of time and resolution. Average values are given in Table 3. Interestingly, the stress monotonically decreases with resolution and there is no evidence for convergence. The resolution dependence is well fit by $\overline{\alpha}\propto N^{-1/3}$.

How large are the error bars on our estimate of $\overline{\alpha}$, and is the observed variation with $N$ significant? We assume that $\alpha(t)$ is a stationary process with mean $\overline{\alpha}$ and variance $\sigma_{\alpha}^{2}$. We provide evidence in the Appendix that the fluctuations in $\alpha(t)$ decorrelate over large time intervals for a long-integration-time $N=32$ model, and that the correlation time $\tau_{c}\Omega\approx 63$. A measurement of $\overline{\alpha}$ averaged over some interval $T$ therefore consists of approximately $\sim T/\tau_{c}$ independent measurements, and one expects an rms error in evaluating $\overline{\alpha}$ of $\approx\sigma_{\alpha}(T/\tau_{c})^{-1/2}$ (see Fig 4 of Longaretti & Lesur, 2010, which implies $\tau_{c}\Omega\sim 10$ in an unstratified local model).

In the Appendix we work out the relation between $\sigma,\tau_{c}$ and the rms error in evaluating $\overline{\alpha}$ for a class of model power spectra, assuming $\alpha(t)$ is a Gaussian process.333The PDF of $\alpha$ is not consistent with a Gaussian. The PDF of $\log\alpha$ is consistent with a Gaussian. The analysis in the appendix does not change if carried out for $\log\alpha$ instead of $\alpha$ For a fit to the $N=32$ run power spectrum, these imply that the expected rms error is $\approx 0.6\sigma_{\alpha}/\alpha_{0}\approx 0.17$, assuming that $\sigma_{\alpha}/\alpha_{0}$ is independent of $N$, consistent with Table 3. This can be compared to $\overline{\alpha}(N)/\overline{\alpha}(2N)-1\approx 0.25$. Therefore the observed trend over a factor of $8$ in $N$ and $\approx 2$ in $\overline{\alpha}$ is significant. A naive estimate of the probability that $d\log\overline{\alpha}/d\log N\geq 0$ gives $\approx 3\%$.

The run of magnetic field energy density $[E_{B}](z,t)=\left[{\bf B}\cdot{\bf B}/2\right]$ for all runs is shown in Figure 3. Evidently the “butterfly” or dynamo oscillations, which are independent, quasi-periodic enhancements in magnetic energy density on either side of the disk, followed by buoyant rise of magnetic field through $z\sim 2H$, are present at all resolutions.

Does the time-averaged vertical structure of the disk change with resolution? Figure 4 shows $x$,$y$-averaged quantities for all runs averaged over time. Also shown are fits to $\rho$ and $E_{B}$ from Guan & Gammie (2011), who study boxes of lower resolution but greater radial and azimuthal extent than we do here.444The fit is $\rho=0.93\rho_{0}\exp(-z^{2}/(2H^{2}))$ for $|z|<2.55H$ and $\rho=0.036\rho_{0}\exp(-(|z|/H-2.55)/0.44)$ otherwise, and $E_{B}=0.012\rho_{0}c_{s}^{2}$ for $|z|<2.55H$ and $0.012\rho_{0}c_{s}^{2}\exp(-(|z|/H-2.55)/0.64)$ otherwise. The density profile is consistent with an exponential profile (rather than Gaussian) at large $|z|$, with scale height $0.44H$. The magnetic energy density is also consistent with an exponential profile at large $|z|$, but with scale height $0.64H$. $E_{B}$ has a feature close to the vertical boundaries, perhaps caused by field lines breaking as they intersect the boundary (Miller & Stone, 2000).

The top right panel of Figure 4 shows the $t$,$x$,$y$-average of total stress. Little variation is seen at large $|z|$, and monotonic decrease of stress with resolution is seen near the midplane. Notice, however, that as resolution increases the structure of averaged stress develops a local minimum around $z=0$ and a local maximum around $|z|\sim 2H$.

3.2 Magnetic Field Correlations

Earlier work (Guan et al., 2009) has shown that the magnetic field correlation length (defined below) scales as $N^{-1}$ in zero-net-field, unstratified local ILES models where $\overline{\alpha}\sim N^{-1}$. How does the characteristic size of structures in MHD disk turbulence change with $N$ for our stratified models?

The dimensionless magnetic field autocorrelation tensor is

[TABLE]

The dimensionless scalar magnetic autocorrelation function $\xi_{B}\equiv Tr(T^{ij})$. Evidently $\xi_{B}(\delta x=0,\delta y=0)=1$. We consider only $\xi_{B}$; $\xi_{v}$ and $\xi_{\rho}$ contain comparatively larger contributions from the compressive disturbances evident in Figure 1 (see also Beckwith et al., 2011).

First we average $\xi_{B}(\delta x,\delta y,z,t)$ over $|z|<2H$ and $t$, as did Davis et al. (2010). The result is shown in Figure 5. The correlation function is an ellipse swept back by the shear into a trailing spiral structure. The shape and orientation of the ellipse do not change significantly with resolution, but the scale of the correlation ellipse drops monotonically as resolution is increased.

Next we average $\xi_{B}(\delta x,\delta y,z,t)$ over time and over bins in $z$ of width $\Delta z=0.5H$, then fold around the midplane. We then evaluate the second moments of $\xi_{B}(z)$ in the contiguous region around $\delta x=\delta y=0$ where $\xi_{B}>0$. The eigenvectors of this moment tensor define a major and minor axis with major axis tilted at a small angle $\theta$ to the $y$ axis. The correlation lengths $\lambda_{\mathrm{minor}}$ and $\lambda_{\mathrm{major}}$ are defined as the distance along each eigenvector at which $\xi_{B}=1/e$. The shape of the correlation departs from an exponential at both small and large scales, although the correlations at large scale are weak and hard to measure accurately (although they must be present, as Guan & Gammie (2011) have shown that butterfly oscillations are coherent over large boxes). The correlation length is the outer scale of disk turbulence.

Figure 6 shows $\lambda_{\mathrm{minor}}(|z|)$, $\lambda_{\mathrm{major}}(|z|)$, and $\theta(|z|)$. All depend on height. The tilt rises toward $\approx 19\degree$ for $|z|<2.5H$. It declines out to $4.5H$ and then rises again toward the boundary (this rise may signal the influence of boundary conditions or the density floor). The major axis correlation length converges toward $\approx 0.5H$ for $|z|>3H$, but is monotonically decreasing with $N$ at $z=0$. The minor axis correlation length is also monotonically decreasing with $N$ at $z=0$, and rises steadily with a bump at $\approx 3H$ toward the boundaries.

Figure 7 shows explicitly the resolution sensitivity $d\log\lambda/d\log N$ for the minor and major axis correlation lengths, along with the resolution sensitivity $d\log\overline{[w_{r\phi}]}/d\log N$ of the shear stress, as a function of $|z|$. Both correlation lengths are sensitive to resolution at the midplane, and far less sensitive (perhaps converged) at higher altitude. At the midplane, both correlation lengths scale as $N^{-1/2}$. $\overline{[w_{r\phi}]}$ exhibits a similar trend, especially for $|z|\lesssim 3H$.

Does this mean the outer scale of turbulence is unresolved, even at our highest resolution? Figure 6 also shows $\lambda_{\mathrm{minor}}(|z|)$, $\lambda_{\mathrm{major}}(|z|)$ in units of $\Delta x$ in the right panels. Above $|z|=3H$ even the minor axis is very well resolved, with in excess of $30$ zones per correlation length. At the midplane $\lambda_{\mathrm{minor}}(N=32)/\Delta x\approx 3$ and $\lambda_{\mathrm{minor}}(N=256)/\Delta x\approx 10$. This differs from the nonconvergence observed in unstratified, zero-net-field ILES models, where $\lambda/\Delta x$ are independent of $N$; here, the outer scale is better resolved as resolution increases.555The ratio of correlation length to resolution $\lambda/\Delta x$ is related to, but not exactly the same as, the quality factor $Q\equiv\lambda_{MRI}/\Delta x$, where $\lambda_{MRI}\propto v_{A}/\Omega$ is a characteristic wavelength for the MRI (Sano et al., 2004; Noble et al., 2010; Hawley et al., 2011). The ratio of the two ratios is $\propto M\equiv v_{A}/\Omega\lambda$, which is the Alfvén Mach number of MRI-driven turbulence at the correlation length. Walker et al. (2016) demonstrated that in their unstratified models $M$ is approximately constant in MRI driven turbulence. In our simulations $M$ varies by a factor $\sim 2$ inside the disk.

3.3 Evolution of Net Magnetic Flux

Our choice of boundary conditions permit evolution of $\langle B^{x}\rangle$ and $\langle B^{y}\rangle$. How important is the mean field in driving the evolution?

The RMS and standard deviation of $\langle B^{y}\rangle$ and $\langle B^{x}\rangle$ are given in Table 4. Evidently $\langle B^{y}\rangle\gg\langle B^{x}\rangle$. We can estimate the effect of the mean field on $\overline{\alpha}$ using the saturation predictor of Hawley et al. (1995) for an unstratified shearing box with a net toroidal field666We emphasize that this predictor is for unstratified models; how well it recovers the behavior of stratified models is uncertain. We also use the mean field through the box as input; locally, the net field may vary.: $B^{2}\sim 4\pi\sqrt{16/15}(0.012)\rho_{0}L_{y}v_{A}\Omega$, where $v_{A}$ is the Alfvén speed associated with the rms net toroidal field. Then using $\overline{\alpha}_{NF}\approx\int dz(1/(2\beta))/\int dzP=0.25B^{2}(5H)/\sqrt{2\pi}$ (where $5H$ comes from assuming $B^{2}=\mathrm{const}$ for $|z|<2.5H$, $B^{2}=0$ else), we find a predicted $\overline{\alpha}$ associated with the mean azimuthal field that is, for all models, at least a factor of $3$ smaller than the measured $\overline{\alpha}$ (and nearly an order of magnitude for r256). This suggests that the boundary conditions are not controlling the saturation.

The mean field sensed locally by the turbulence may still control $\alpha$ locally. To illustrate this point, Figure 8 shows a sample estimate of a local mean field: the azimuthal field averaged over sheets at constant $z$. This fluctuates in sign, so to avoid cancellation we take the time average of the absolute value of this mean field. The resulting mean field is an order of magnitude larger than $\langle B_{y}\rangle$, which the unstratified box saturation predictor suggests would produce an $\overline{\alpha}$ comparable to what is measured. In sum: a localized mean field may play an important role in controlling the outcome, but the mean field over the entire computational domain does not.

4 Discussion

Our simulations have yielded several unexpected dependences on resolution: (1) $\overline{\alpha}\sim N^{-1/3}$ (2) $\lambda_{\mathrm{minor}}\propto\lambda_{\mathrm{major}}\sim N^{-1/2}$ in the midplane of the disk, and nearly $\sim N^{0}$ in the corona (3) The total stress, scaling similarly to $\lambda$ with $N$, develops a local maximum at $|z|\sim 2H$ as resolution increases.

Surprisingly, we do not see convergence of the time-averaged, vertically integrated shear stress $\overline{\alpha}$ for resolution up to $N=256$ zones per scale height in stratified, isothermal, local ILES models. This is broadly consistent with Bodo et al. (2014) and in tension with the results of Davis et al. (2010). Both Davis et al. (2010) and Bodo et al. (2014) find a plateau in $\overline{\alpha}$ between $\approx 45$ and $\approx 90$ zones per $H$. We do not find evidence for this behavior, but the plateau could be hidden in our measurement errors due to finite run time and finite computational volume.

Are our results consistent with earlier work? To compare, we need to convert to common units and a common measurement of stress, for which we will use $\alpha$ as defined in eq.(14).

Davis et al. (2010) report volume-and-time averaged stresses in units of the midplane pressure. This is equivalent to volume-averaged stress in our units. Notice that Davis et al. (2010) define $H=\sqrt{2}c_{s}/\Omega$. Then for $N\approx(23,$ $45,$ $91)$ their volume averaged stress (see their Table 1) is $(0.0149,$ $0.0093,$ $0.0092)$. Converting to vertically integrated stress (multiply by $4\sqrt{2}$) and dividing by the vertically integrated gas pressure ($\sqrt{2\pi}$ in our units), we find (using our definition of $\alpha$) $\overline{\alpha}=(0.034,$ $0.021,$ $0.021)$.

Bodo et al. (2014) also define $H=\sqrt{2}c_{s}/\Omega$, and set $c_{s}=1/\sqrt{2}$, $\rho_{0}=1,$ and $\Omega=1$, so their unit of stress is a factor of $2$ larger than ours. They consider models with $N\approx(23,$ $45,$ $91,$ $141)$. Since they do not report time-averaged stresses, we will estimate these from their Figure 2. We estimate that the volume integrated maxwell stress in their units is $\simeq(0.022,$ $0.017,$ $0.017,$ $0.01)$. We convert this vertically integrated stress to our units (multiply by $2\sqrt{2}$; the factor of $2$ is for the stress unit and the factor of $\sqrt{2}$ is for the length unit), multiply by $1.25$ to incorporate an assumed $25\%$ Reynolds stress contribution, then divide by the vertically integrated pressure ($\sqrt{2\pi}$) in our units to find $\overline{\alpha}\simeq(0.031,$ $0.024,$ $0.024,$ $0.014)$.

To facilitate comparison, at a resolution of $N=(32,$ $64,$ $128,$ $256)$ we find $\overline{\alpha}=(0.039,$ $0.034,$ $0.025,$ $0.019)$. These results are shown in Figure 9. The overall offset of the Davis et al. and Bodo et al. series from ours is significant, but may be explained in part by the larger vertical extent of our models. The algorithms used also differ, possibly yielding different effective resolutions, and of course the vertical boundary conditions also differ. Nevertheless, it is reassuring that all simulations lead to values of $\overline{\alpha}$ that are within $1\sigma$ of our results. Indeed, least squares power-law fits to the Davis et al. and Bodo et al. series yield slopes ($-0.35$ and $-0.37$, respectively) consistent with ours ($-0.36$) and the relationship $\overline{\alpha}\sim N^{-1/3}$.

The correlation function in the $x-y$ plane is approximately ellipsoidal and characterized by the major axis length, minor axis length, and the “tilt angle” between the major axis and the $y$ axis. The tilt angle $\theta_{\mathrm{tilt}}(N=256)\sim 19\degree$ is consistent with Davis et al. (2010) who find $\theta_{\mathrm{tilt}}(N\approx 91)\sim 18\degree$. The increase in $\theta_{\mathrm{tilt}}$ with resolution was also reported by Guan et al. (2009). Although Davis et al. (2010) do not quote a value for $\lambda_{\mathrm{minor}}$, visual inspection of their $\xi_{B}(N\approx 64)$ slice yields a value comparable to what we find at similar resolution.

The sensitivity of stress to $N$ depends on height (see Figure 7). The midplane shear stress decreases with $N$ at a rate that is inconsistent with convergence, but the stress at $|z|\gtrsim 2H$ is much less sensitive to $N$ and convergence is not excluded by our limited time- and volume-sampled data. One consequence of this is that a local minimum develops in the total stress at $z=0$ and a local maximum develops at $|z|\simeq 2H$. A qualitatively similar local maximum in the stress is observed in stratified shearing boxes with self-consistent thermodynamics, at least when they are radiation pressure-dominated (Hirose et al. 2009, Jiang et al. 2016). This effect appears to be due to a convective process which also significantly enhances $\alpha$ in these models (e.g. Hirose et al., 2014).

Are our simulations run long enough? From a long-duration, low-resolution simulation we measured a correlation time of $\approx 60\Omega^{-1}$ (this is slightly shorter than the $90\Omega^{-1}$ correlation time seen in the $N=90$ run of Davis et al. (2010) 777We thank S. Davis for kindly providing us with the data.), and our assessment of the error bars on $\bar{\alpha}$ relies on this measurement. Stratified shearing box models frequently give an impression of order-unity enhancements in $\alpha$ (“bursts”) separated by long intervals, and rare bursts could change the correlation time. Our data are not sufficient to assess whether this impression is statistically well grounded or not. If it is, then the bursts might correspond to long-timescale power in the power spectrum of a Gaussian process that is undetectable in a short simulation, or non-Gaussianity associated with the flares. There is, however, no evidence for non-Gaussianity in our data; the probability distribution for $\log{\alpha}$, for example, is consistent with Gaussian. There is also no evidence to changes in the variance of $\log(\alpha)$ with $N$; the relative variance, shown in Table 3, shows no systematic trend.

Why no convergence? The cause may lie either with our numerical realization of the stratified isothermal zero-net-flux ILES shearing box model (A), or with assumptions made by the model itself (B). We have assembled an incomplete list of possible explanations:

(A1) The nonconvergence is physical and $\alpha_{MRI}\rightarrow 0$ in isothermal astrophysical disks with vanishing mean field. Although we cannot rule this out, it seems inconsistent with the result of Fromang (2010) for an unstratified model with explicit scalar viscosity and resistivity that converges to nonzero $\alpha$, albeit only for $P_{m}=4$.

(A2) The apparent nonconvergence is a consequence of a combination of statistical errors associated with a finite sampling time and an initial transient that results from using resolution $N/2$ data to initialize resolution $N$ models. Our analysis (see Appendix) suggests, however, that even though $\overline{\alpha}$ has a long correlation time this is improbable.

(A3) The nonconvergence is an artifact of the limited size of the model. Fluctuations in $\alpha(t)$ will depend on the volume of the simulation. Naively, they would scale as $1$ over the square root of the number of correlation volumes. But there is coupling between correlation volumes via large-scale magnetic fields and this is connected to the butterfly oscillations. Furthermore, it is already known that in unstratified, local simulations the imposition of a mean field causes an ILES model to converge. Ultimately, it must be that turbulence is locally unable to distinguish between uniform fields and magnetic fields that have structure on a sufficiently large scale. Perhaps our models are simply too small to see this sufficiently large scale, and so they are analogous to the zero mean field unstratified models that do not converge.

The interaction between small and large scale fields has been explored by Sorathia et al. (2012) who measured the net magnetic flux in local regions of global unstratified ILES models. They found distributions of $\langle B^{x}\rangle$, $\langle B^{y}\rangle$, and $\langle B^{z}\rangle$ inconsistent with zero, with the linear MRI growth associated with these mean fields typically being well resolved in their simulations.

(A4) The model will converge at higher $N$, and we are simply not in the high resolution limit yet. The magnetic field correlation length is $\sim 10\Delta x$ in our highest resolution models, so there is only a dynamic range of $\simeq 2$ between the outer scale $\lambda$ and the dissipation scale.

(A5) The model will not converge, with $\alpha\sim\lambda\sim N^{-1}$, in the complete absence of a mean magnetic field. Although we cannot account for the $\overline{\alpha}$ we measure from the net flux through our computational domain, our estimate is based on a fit to results from unstratified models. Mean fields in stratified boxes may behave differently. They may, for example, be playing a stronger role than we estimate in the magnetically-dominated corona, contributing to our near-convergence of $\lambda$ above $|z|\sim 3H$. Nonetheless, Davis et al. (2010) do maintain zero net flux in a stratified model, and their results are inconsistent with $\alpha\sim N^{-1}$.

(B1) The nonconvergence is an artifact of our use of an ILES model. In models in which the numerical resolution and Reynolds numbers are increased together, there is numerical evidence that both unstratified and stratified models converge (Fromang 2010, Simon et al. 2011). It would be interesting to know whether this extends to larger $N$ and the large $Re_{M}$, large $Re$ limit relevant to astrophysical disks. There is also numerical evidence that computational models of the solar dynamo depend strongly on the dissipation model (see, e.g. Charbonneau, 2014, for a review).

(B2) The nonconvergence is an artifact of the absence of consistent vertical energy transport by radiation and convection. It is now known that convective disks in models with consistent treatment of energy transport exhibit enhanced $\alpha$ (Hirose et al. 2014), which may enhance the amplitude of dwarf nova outbursts. The convective process may aid convergence (Bodo et al. 2015). It is not yet clear how well converged the energetically consistent models are; current model have $N\sim O(64)$.

(B3) The nonconvergence is an artifact of the symmetry of the local model. The local model is invariant under translations in the plane of the disk, and invariant under rotations by $\pi$ around the $z$ axis. The incorporation of higher order terms in $H/r$ would break these symmetries and might qualitatively change the outcome. There is limited numerical evidence for convergence in unstratified global models (Sorathia et al., 2012), although with a tendency for $\overline{\alpha}$ (and hence $\dot{M}$ and $\beta^{-1}$) to increase with resolution (see also Shiokawa et al., 2012; Hawley et al., 2011, 2013).

We are unsure which (if any) of these explanations is correct, but all except (A1) are amenable to future numerical investigation.

What are the implications of nonconvergence? It is difficult to say without testing the hypotheses above with new numerical simulations. For example, if A4 is correct (insufficient resolution) then current lower resolution models may yield $\bar{\alpha}$ to within a factor of two. On the other hand, if A5 is correct ($\alpha$ is zero without a mean field) then the result would have profound implications for our understanding of disk structure and evolution, which would presumably be controlled by the generation and transport of large-scale magnetic field. No matter what the explanation for the nonconvergence seen here, future disk simulations need to be tested carefully for convergence.

5 Conclusion

The isothermal stratified zero-net-flux shearing box is a minimal model with zero physical parameters for the turbulent saturation of the magnetorotational instability and is thus central to accretion disk theory. We have attempted to sort out apparently conflicting reports of convergence in the literature using the ramses-gpu code on blue waters to probe convergence at an unprecedented resolution of $N=256$ zones per scale height.

Our results imply that existing local and perhaps global zero-mean-field ILES models of disks are, at best, underresolved. We have found that $\overline{\alpha}\sim N^{-1/3}$. This is not convergent, but it differs from the sharp nonconvergence identified by Fromang & Papaloizou (2007) in unstratified ILES models, with $\overline{\alpha}\sim N^{-1}$.

We have also compared our results to earlier work by Davis et al. (2010) and Bodo et al. (2014). These earlier calculations are consistent with our to within the error bars, and all show a similar trend with resolution. Like Bodo et al. (2014) and unlike Davis et al. (2010), our models do not conserve net toroidal magnetic flux. Although first estimates suggest the net flux present in our model is not controlling our results, this remains an uncertainty in performing comparisons. Box size effects may also confound comparisons.

We have reviewed possible physical and numerical causes of this nonconvergence. All of these are amenable to further numerical investigation when sufficient computational resources are available. One implication is clear, however: simulations of MHD turbulence in disks need to be tested carefully for convergence, and the attendant uncertainties need to be allowed for when weighing the results.

It is a pleasure to thank G. Bodo, M. Chandra, S. Davis, J. Dolence, V. Paschalidis, J. Simon, and J. Stone for useful discussions. The numerical calculations presented here were performed on the Blue Waters supercomputer at NCSA. BRR was supported by an Illinois Distinguished Fellowship and by NSF grant AST-1333612. Work at Los Alamos National Laboratory was done under the auspices of the National Nuclear Security Administration of the US Department of Energy. CFG’s work was also supported in part by a Romano Professorial Scholar appointment, a Simons Fellowship in Theoretical Physics, and a Visiting Fellowship at All Souls College, Oxford.

Appendix A Measurement Error Estimates with a Gaussian Process Model

Shearing box simulations estimate the true, long-term average $\alpha_{0}$ from $\overline{\alpha}$ measured over a finite time $\Delta t$. How long is long enough?

Note that in this section, $\langle\rangle$ denotes an expectation value for consistency with previous literature on Gaussian random fields, rather than the volume average of Equation 11. Suppose $\alpha$ has a correlation time $\tau_{c}$ and variance $\sigma^{2}$. Then our intuition is that $\langle(\overline{\alpha}-\alpha_{0})^{2}\rangle$should be proportional to $\sigma^{2}\tau_{c}/\Delta t$, i.e. the rms error averaged over many realizations of $\alpha$ should scale as one over the square root of the number of correlation times. But with what coefficient?

We can estimate $\langle(\overline{\alpha}-\alpha_{0})^{2}\rangle$ for a Gaussian process with known power spectrum over some long but finite time $T$. That is,

[TABLE]

The sum is taken only over $\omega_{j}>0$, $\phi_{j}$ is uniformly distributed in $[0,2\pi)$ (random phase) and $a_{j}$ is Gaussian distributed:

[TABLE]

The power spectrum $P_{\omega}\equiv T\langle a_{j}^{2}\rangle$. In the limit that $T$ is large the modes are closely spaced and $\sum_{j}\rightarrow\int d\omega T/(2\pi)$. The expected variance in $\alpha$ over the interval $T$ is

[TABLE]

where the factor of $2$ comes from phase averaging. $P_{\omega}$ is independent of $T$ if $\sigma^{2}$ is fixed.

The autocorrelation function is

[TABLE]

The error in estimating $\alpha_{0}$ from a finite interval $\Delta t$ is

[TABLE]

Expanding and integrating,

[TABLE]

Then

[TABLE]

To go further we need to know the power spectrum.

We consider model power spectra that decorrelate on long timescales, so that $P_{\omega}\propto\omega^{0}$ for $\omega$ small, and scale as a power law at high frequency. A suitable model is

[TABLE]

Evidently if the process is stationary then $p>1$. The power spectrum can be normalized by $\sigma^{2}$:

[TABLE]

Then

[TABLE]

where $K_{n}$ is a modified Bessel function of the second kind. It is easy to show that $\xi(0)\rightarrow\sigma^{2}$.

We estimate $\sigma^{2}$ from data taken over an interval $\Delta t$. This estimate is biased because it does not include contributions to the variance from low frequency components. The expected value of $\sigma^{2}$ sampled over time $\Delta t$ is

[TABLE]

If $p$ and $\omega_{0}$ are known then this expression can be used to produce an debiased estimate of $\sigma^{2}$.

An auxiliary $N=32$ run with $\Delta t=2000\Omega^{-1}$ has a power spectrum consistent with $p\sim 2$. For this special case, $\xi(\tau)=\sigma^{2}e^{-\omega_{0}|\tau|}$, $\tau_{c}=\omega_{0}^{-1}$, and

[TABLE]

Consistent with expectations, this scales as $\sigma^{2}\tau_{c}/\Delta t$. Also,

[TABLE]

The auxiliary $N=32$ run has $\tau_{c}\Omega\simeq 60$, so $\omega_{0}=60^{-1}\Omega$, which we will assume is independent of $N$. Runs in Table 2 with $\Delta t=300\Omega^{-1}$ therefore have $\omega_{0}\Delta t\simeq 5$. Runs in Table 2 also have $\sigma_{\alpha}/\overline{\alpha}\simeq 0.25$. Then (A13) implies the debiased $\sigma_{\alpha}/\alpha\simeq 0.30$. Combined with (A12), we find $\langle m^{2}\rangle^{1/2}=0.17\alpha$. This implies that the one-sigma error is small compared to the total change in $\overline{\alpha}$ over a factor of $8$ in $N$.