Global simulations of magnetorotational turbulence I: convergence and the quasisteady state
Abstract
Magnetorotational turbulence provides a viable mechanism for angular momentum transport in accretion disks. We present global, three dimensional (3D), magnetohydrodynamic accretion disk simulations that investigate the dependence of the turbulent stresses on resolution. Convergence in the timeandvolumeaveraged stresstogaspressure ratio, , at a value of , is found for a model with radial, vertical, and azimuthal resolution of 1251, 27, and 12.5 cells per scaleheight (the simulation mesh is such that cells per scaleheight varies in the radial direction). The gas pressure dependence of the quasisteady state stress level is also examined using models with different scaleheighttoradius aspect ratio (), revealing a weak dependence of on pressure.
A control volume analysis is performed on the main body of the disk () to examine the production and removal of magnetic energy. Maxwell stresses in combination with the mean disk rotation are mainly responsible for magnetic energy production, whereas turbulent dissipation (facilitated by numerical resistivity) predominantly removes magnetic energy from the disk. Recasting the magnetic energy equation in terms of the power injected by Maxwell stresses on the boundaries of, and by Lorentz forces within, the control volume highlights the importance of the boundary conditions (of the control volume). The different convergence properties of shearingbox and global accretion disk simulations can be readily understood on the basis of choice of boundary conditions and the magnetic field configuration. Periodic boundary conditions restrict the establishment of largescale gradients in the magnetic field, limiting the power that can be delivered to the disk by Lorentz forces and by stresses at the surfaces. The factor of three lower resolution required for convergence in for our global disk models compared to stratified shearingboxes is explained by this finding.
keywords:
accretion, accretion disks  magnetohydrodynamics  instabilities  turbulence1 Introduction
For the astrophysically common process of mass accretion through a disk to be effective, outward angular momentum transport must occur (Shakura & Sunyaev 1973; Pringle 1981). In the past two decades it has become clear that selfsustaining magnetized turbulence driven by the magnetorotational instability (MRI) can play this role (Balbus & Hawley 1998).
Due to the highly nonlinear nature of magnetorotational turbulence, numerical simulations have become a common tool in its study. These simulations come in a number of flavours: unstratified shearingboxes (where the vertical component of gravity is neglected)(Hawley et al. 1995; Fromang & Papaloizou 2007; Fromang et al. 2007; Lesur & Longaretti 2007, 2011; Lesaffre et al. 2009; Latter et al. 2009; Heinemann & Papaloizou 2009; Simon et al. 2009; Guan et al. 2009; Bodo et al. 2011; Käpylä & Korpi 2011; Latter & Papaloizou 2012), stratified shearingboxes (Brandenburg et al. 1995; Stone et al. 1996; Miller & Stone 2000; Fleming et al. 2000; Brandenburg 2005; Johansen et al. 2009; Gressel 2010; Shi et al. 2010; Davis et al. 2010; Simon et al. 2011; Guan & Gammie 2011; Oishi & Mac Low 2011; Simon et al. 2012), unstratified global models (Hawley 2001; Armitage et al. 2001; Nelson & Gressel 2010; Sorathia et al. 2012), and stratified global models (Hawley 2000; Hawley & Krolik 2001; Arlt & Rüdiger 2001; Fromang & Nelson 2006, 2009; Beckwith et al. 2008; Lyra et al. 2008; Sorathia et al. 2010; O’Neill et al. 2011; Flock et al. 2011, 2012a; Noble et al. 2010; Beckwith et al. 2011; Hawley et al. 2011; McKinney et al. 2012; Parkin & Bicknell 2013). Shearingbox simulations focus on a local patch of an accretion disk whereas global simulations have the potential to study the entire radial (and vertical) extent of an accretion disk. Despite these numerous different approaches to modelling accretion disk turbulence, similarities exist in the magnetorotational turbulence that they exhibit. In general, there is an initial phase where the MRI develops and transient magnetic field amplification arises, following which the growth of stresses subsides and the disk settles into a quasisteady state (QSS).
There have been mixed results from simulations as to what sets the QSS stress level. The results of unstratified shearingbox simulations by Fromang & Papaloizou (2007) (see also Lesur & Longaretti 2007; Fromang et al. 2007; Simon et al. 2009; Guan et al. 2009; Fromang 2010; Käpylä & Korpi 2011) show that dissipation (i.e. resistivity and viscosity) dictates the QSS stress level. When this dissipation is purely numerical in origin, increasing the simulation resolution causes a reduction in the volume averaged stress in zeronet flux, unstratified shearingbox simulations. Fromang & Papaloizou (2007) argue that this occurs because magnetorotational turbulence always drives energy to the smallest resolved scale, thus removing energy from the larger (angular momentum transporting) eddies. Sorathia et al. (2012) have recently revisited this issue using unstratified global disks, revealing a contrasting result of converged stresses with increasing resolution. What then sets the QSS stress level? Vishniac (2009) has argued that stratification, if present, will affect the QSS stress, and it is indeed found that including stratification facilitates convergence (Davis et al. 2010; Shi et al. 2010; Oishi & Mac Low 2011). Furthermore, including a net flux field in unstratified shearingboxes enables convergence (e.g. Simon et al. 2009). Considering the aforementioned results, there is a clear indication that the choice of numerical setup and/or magnetic field configuration play crucial roles.
In this work we take the logical next step and investigate convergence in stratified global disk models, which is indeed found, but for lower resolutions than in equivalent shearingbox simulations. A complementary analysis of magnetic energy production leads us to conclude that boundary conditions have a profound influence on the QSS stress. The remainder of this paper is organised as follows: in § 2 we describe the simulation setup and diagnostics used in this investigation. In § 3 we examine the dependence of the saturated turbulent state on simulation resolution and disk scaleheight. The results from the application of a controlvolume analysis to the simulations are presented in § 4. We discuss our findings in the context of a unified interpretation for magnetorotational turbulence in different numerical setups in § 5 and close with conclusions in § 6.
2 The model
2.1 Simulation code
The timedependent equations of ideal magnetohydrodynamics are solved using the PLUTO code (Mignone et al. 2007) in a 3D spherical coordinate system. The relevant equations for mass, momentum, energy conservation, and magnetic field induction are:
(1)  
(2)  
(3)  
(4) 
Here is the total energy, is the internal energy, is the velocity, is the mass density, is the pressure, and is the magnetic energy. We use an ideal gas equation of state, , with an adiabatic index . The adopted scalings for density, velocity, temperature, and length are, respectively,
where is the speed of light, and the value of corresponds to the gravitational radius of a black hole.
The gravitational potential due to a central point mass situated at the origin, , is modelled using a pseudoNewtonian potential (Paczyńsky & Wiita 1980):
(5) 
Note that we take the gravitational radius (in scaled units), . The Schwarzschild radius, for a spherical
black hole and the innermost stable circular orbit (ISCO) lies at
. The term on the RHS of Eq (3) is an
adhoc cooling term used to keep the scaleheight of the disk
approximately constant throughout the simulations; without any
explicit cooling in conjunction with an adiabatic equation of state,
dissipation of magnetic and kinetic energy leads to an increase in gas
pressure and, consequently, disk scaleheight over time. The form of
used is identical to that of Parkin &
Bicknell (2013); further
details can be found in that paper
The PLUTO code was configured to use the fivewave HLLD Riemann solver of Miyoshi & Kusano (2005), piecewise parabolic reconstruction (PPM  Colella & Woodward 1984), limiting during reconstruction on characteristic variables (e.g. Rider et al. 2007), secondorder RungeKutta timestepping, and the upwind CONTACT Constrained Transport scheme of Gardiner & Stone (2008) (to maintain ) which includes transverse corrections to interface states. This configuration was found to be stable for linear MRI calculations by Flock et al. (2011).
The grid used for the simulations is uniform in the and directions and extends from and . A graded mesh is used in the direction which is uniform within and stretched between , where is the thermal disk scaleheight. For our fiducial model, gblsr, there are a total of 170 cells in the direction, of which 108 are uniformly distributed within , and the remaining 62 cells on the stretched sections between . Details of the grid resolutions used in the simulations are provided in Table 1. The adopted boundary conditions are identical to those used in Parkin & Bicknell (2013). Finally, floor density and pressure values are used which scale linearly with radius and have values at the outer edge of the grid of and , respectively.
2.2 Initial conditions
The simulations start with an analytic equilibrium disk which is isothermal in height (, where is the temperature) and possesses a purely toroidal magnetic field. The derivation of the disk equilibrium and a detailed description of the initial conditions can be found in Parkin & Bicknell (2013). In cylindrical coordinates (), the density distribution, in scaled units, is given by,
(6) 
where the pressure, , and the ratio of gastomagnetic pressure, is initially set to 20 in all models. For the radial profiles and we use simple functions inspired by the Shakura & Sunyaev (1973) disk model, except with an additional truncation of the density profile at a specified outer radius:
(7)  
(8) 
where sets the density scale, and are the radius of the inner and outer disk edge, respectively, is a tapering function (Parkin & Bicknell 2013), and and set the slope of the density and temperature profiles, respectively. In all of the global simulations , , , , and . In § 3, models with aspect ratios of and 0.1 are considered. These ratios are achieved by setting and , respectively. The rotational velocity of the disk is close to Keplerian, with a minor modification due to the gas and magnetic pressure gradients,
(9) 
where,
(10) 
The region outside of the disk is set to be an initially stationary, spherically symmetric, hydrostatic atmosphere. The transition between the disk and background atmosphere occurs where their total pressures balance. To initiate the development of turbulence in the disk, a low wavenumber, nonaxisymmetric Fourier mode is excited in the poloidal velocities with amplitude , where is the sound speed.
2.3 Diagnostics
A volumeaveraged value (denoted by angled brackets) for a variable is computed via,
(11) 
Similarly, azimuthal averages are denoted by square brackets,
(12) 
Time averages receive an overbar, such that a volume and time averaged quantity reads . Throughout this paper we concentrate on the region between and , where and (where is the sound speed). We define this region as the “disk body”.
The efficiency of angular momentum transport is typically quantified from the total stress,
(13) 
where the Reynolds stress tensor,
(14) 
and the Maxwell stress tensor,
(15) 
The largest contribution comes from the component of (Brandenburg et al. 1995; Hawley et al. 1995; Stone et al. 1996),
(16) 
where we have defined the perturbed flow velocity as
(17) 
Furthermore, we calculate the component of the Maxwell stress normalised by the magnetic pressure,
(18) 
We follow Noble et al. (2010) and Hawley et al. (2011) and utilize a “quality factor” to measure the ability of the simulations to resolve the wavelength of the fastest growing MRI mode, . Defining,
(19) 
where , and is the Alfvén speed, the “quality factor” is given by,
(20) 
where is the cell spacing in direction . The “resolvability”  the fraction of cells in the disk body that have (e.g. Sorathia et al. 2012)  is then defined as,
(21) 
where represents a cell.
Model  Resolution  

(, , )  ()  
gbllr  0.1  340,112,128  8.536  18  8 
gblsr  0.1  512,170,196  1251  27  12.5 
gblhr  0.1  768,256,256  1877  37  16 
gbllrla  0.1  420,140,70  10.545  20  4.5 
gblhrla  0.1  768,256,128  1877  37  8 
gblthin  0.05  512,170,320  625  27  10 
2.4 Fourier analysis
The simulation data is Fourier transformed in spherical coordinates to compute power spectra for different simulation variables. A detailed description of the method used is given in Appendix A. In brief, we define the Fourier transform of a function as,
(22) 
It then follows that the angleaveraged (in Fourier space) amplitude spectrum,
(23) 
where an asterisk () indicates a complex conjugate. The total power at a given wavenumber  the power spectrum  is then given by .To compute a power spectrum, we take 300 simulation checkfiles (equally spaced over the last , where is the orbital period at a radius ) and compute 30 powerspectra, each timeaveraged over . These 30 realisations are then averagedover to produce the final power spectrum.
Model  

gbllr  1231  0.57  0.27  0.67  131  395  14  12  0.043  0.39 
gblsr  1231  0.69  0.43  0.75  128  361  17  14  0.040  0.42 
gblhr  1231  0.78  0.56  0.80  123  332  18  15  0.039  0.42 
gbllrla  1831  0.31  0.07  0.25  655  2296  34  32  0.013  0.29 
gblhrla  1231  0.76  0.50  0.61  150  430  18  15  0.035  0.39 
gblthin  1231  0.52  0.44  0.75  120  416  15  13  0.044  0.41 
2.5 Summary of models
In Table 1 we list six simulations aimed at investigating the following points:

Convergence with resolution: Models gbllr, gblsr, and gblhr are low, standard, and high resolution variants, respectively, with identical cell aspect ratio and .

Importance of azimuthal resolution: Model gbllrla (gblhrla) is identical to gbllr (gblhr) with the exception of a lower azimuthal resolution (denoted by the affix “la”).

Scaleheight dependence: Models gblsr and gblthin have disk scaleheights of and 0.05, respectively. These models feature an identical number of cells per scaleheight in the vertical and azimuthal directions.




3 The quasisteady state
Following the initial transient phase of evolution, the disk settles into a QSS. In this section we examine the characteristics of this state for the different simulations. We list time and volume averaged parameter values pertaining to the steadystate turbulence in Table 2.



3.1 Resolution dependence
Convergence: gbllr, gblsr, and gblhr
The volume averaged stress normalised to gas pressure, , displays a dependence on resolution in the magnitude of the transient peak at (Fig. 1). Following this, steadily declines until , at which point the curves leveloff and a quasisteady state is reached. We find a timeaveraged value during the quasisteady state of for models gbllr, gblsr, and gblhr, indicating convergence with resolution. At the point of convergence, , where is the disk body plasma. This is in agreement with, but slightly higher than, the relation found for unstratified shearingbox simulations (Sano et al. 2004; Blackman et al. 2008; Guan et al. 2009). The volume averaged Maxwell stress normalised to the magnetic pressure, converges at a value of 0.42 (lower panel of Fig. 1), consistent with previous high resolution stratified shearing box (Hawley et al. 2011; Simon et al. 2012) and global disk models (Parkin & Bicknell 2013).
Convergence in is coincident with convergence in the resolvability (Eq 21)  the ability of the numerical grid to resolve the fastest growing MRI modes. Examining Fig. 2, one sees that the direction is the best resolved, followed by the radial direction, and then the direction. At the resolution of model gblsr, is clearly higher than , suggesting that convergence in global models is tied to the radial magnetic field. That convergence with resolution is more readily achievable for the radial magnetic field is illustrated by the relative magnetic field strengths: , , and . The converged value for is roughly 0.8, consistent with some fraction of the disk having weak magnetic fields (due to zeronet flux dynamo oscillations) which have corresponding values which are below the simulation resolution. Poloidal magnetic fields in our simulations appear relatively strong, with our models returning and compared to respective values of and for the highest resolution model in Hawley et al. (2011). We attribute this difference to the higher resolution used in our models.
Power spectra computed for density, kinetic energy, and total magnetic energy perturbations are shown in Fig. 3. A perturbation for a variable q is calculated by subtracting the azimuthal average, such that . Perturbations are used to reduce the influence of large scale radial and vertical gradients on the resulting power spectra. The position of the low wavenumber turnover between models gbllr, gblsr, and gblhr is consistent at , illustrating that the resolution of the largest physical structures is converged. The slope of the magnetic energy power spectrum is approximately . This apparent constant powerlaw slope suggests a selfsimilar transfer of energy from large to small scales which may indicate an inertial cascade, although it may also be due to the injection of energy by the MRI at all realisable scales (Fromang & Papaloizou 2007). The power spectra for magnetic energy and kinetic energy exhibit very similar shapes. On further inspection one sees that magnetic energy is slightly larger than kinetic energy on length scales of roughly a disk scaleheight (, where ), whereas they are approximately equal on the smallest length scales (). This differs from the stratified shearingbox simulations presented by Johansen et al. (2009), for which kinetic energy was found to dominate over magnetic energy at all but the very largest scales in the box. Comparing to the global models of Beckwith et al. (2011), we note that the low wavenumber turnover for magnetic and kinetic energy arises at a similar value (they find ). However, there is a considerable difference in the amplitudes of kinetic and magnetic energy fluctuations, where Beckwith et al. (2011) find the latter to be an order of magnitude lower than the former. The source of this difference is unclear. However, there are number of differences between our approach and that used by Beckwith et al. (2011) in the calculation of the Fourier transforms and the related power spectra. In calculating the value of a fluctuating quantities we have adopted a straightforward approach of subtracting an azimuthally averaged value of , whereas Beckwith et al. (2011) fit a twodimensional distribution in the radial and vertical directions, which they then subtract to determine . Another major difference is that we define a conventional spherical Fourier transform through Eq (22), whereas Beckwith et al. (2011) construct azimuthal averages of fluctuating quantities, define a normalized measure of spatial fluctuations in the coordinates and then define a Fourier transform in and treating and as pseudoCartesian coordinates (their equation (15)). A comparison of these two approaches and the implications for comparing computed accretion disk spectra with one another and also with textbook spectra for homogeneous turbulence, is beyond the scope of this paper.
The power spectra all display a pronounced turnover at high wavenumber, which depends on resolution, and which we interpret as the dissipation scale. Indeed the morphology of the steep, but slightly curved, step at the high wavenumber end of the magnetic energy power spectrum is indicative of a resolved separation between the Ohmic and viscous dissipation scales (see, e.g., Kraichnan & Nagarajan 1967).
Irrespective of resolution, most of the power in is on the largest length scales (i.e. at low wavenumber  Fig. 4). In fact, the relatively flat slope to the power spectrum indicates that a large amount of power is also contained in moderate length scales. The slope of the power spectrum changes at , becoming steeper, and indicating that smaller length scales contribute considerably less to the global stress. Therefore, although magnetic field correlation lengths demonstrate that MRIdriven turbulence is localised (Guan et al. 2009), we find evidence for angular momentum transport being dominated by larger length scales, of size .
Increasing the simulation resolution permits structure to occupy smaller spatial scales. This is illustrated by the simulation snapshots of shown in Fig. 5. As one progresses to higher resolution through models gbllr, gblsr, and gblhr the size of structures get progressively finer. Also, contrasts in the magnetic energy, which are particularly noticeable in the coronal region, become sharper at higher resolution. This equates to an increase in with resolution, which we examine in more detail in § 4.
In summary, convergence is achieved for a resolution of 1251 cells/ in radius, 27 cells/ in the direction, and 12.5 cells/ in the direction (model gblsr). This is considerably below the 64128 cells/ required for convergence in stratified shearing box simulations found by Davis et al. (2010), whereas the vertical resolution is comparable to the 25 cells/ necessary to produce sustained turbulence in the models of Fromang & Nelson (2006) and Flock et al. (2011). In § 5 we provide an explanation for this dramatic difference.
Influence of : gbllrla and gblhrla
When the azimuthal field is underresolved, turbulent activity dies out, as discussed by Fromang & Nelson (2006), Flock et al. (2011), and Parkin & Bicknell (2013). This effect can be seen in the stresses and resolvabilities computed for model gbllrla (a lower azimuthal resolution variant of gbllr  see Table 1), which we plot in Figs. 1 and 2. Repeating this experiment at higher resolution (models gblhr and gblhrla), one finds that even though the azimuthal field is barelyresolved (8 cells/ in the azimuthal direction for gblhrla), only a slightly lower value is obtained. Therefore, we find a similar dependence on azimuthal resolution to that discussed by Hawley et al. (2011), although this dependence appears to become less pronounced at higher resolution, and this is possibly due to compensation by the poloidal grid resolution. This indicates that low azimuthal resolution can, to some extent, be compensated for by higher poloidal resolution. However, based on these results it would seem advisable to adopt an aspect ratio close to unity. We also note that our and values are higher than in the models of Beckwith et al. (2011) and Hawley et al. (2011). We attribute this to the higher simulation resolution and lower cell aspect ratio used in our models  (see also the discussion in Fromang & Nelson 2006; Flock et al. 2011; Parkin & Bicknell 2013).
As discussed in the previous section, the power spectra in Fig. 3 display a turnover at high wavenumber corresponding to the dissipation scale. All simulations presented in this paper rely on numerical dissipation, hence one may anticipate that numerical resolution sets this scale and adopting a lower resolution in a certain direction may shift the dissipation scale to lower wavenumbers. Comparing the curves for models gbllr and gbllrla in Fig. 3, the slope in the magnetic energy power spectrum is steeper for gbllrla in the wavenumber range . This steeper slope is readily understood as a consequence of underresolving the fastest growing MRI modes in the direction  energy cannot be injected by the MRI if the modegrowth is not resolved. In relation to the discussion in the previous section, this implies that the powerlaw slope in the magnetic energy power spectrum of model gbllr, for example, is a consequence of magnetic energy injection by the MRI and not solely due to an inertial cascade of energy from large to small scales, consistent with results presented by Fromang & Papaloizou (2007) and Johansen et al. (2009) which illustrate the driving of magnetic energy to smaller scales by the MRI.
3.2 Other factors which might affect the saturated state
Gas pressure dependence: gblthin
Based on a suite of unstratified shearingbox simulations,
Sano et al. (2004) have reported a dependence of the QSS stress level on
the gas pressure. To examine whether this dependence exists in global
simulations, one can utilise simulations with different aspect ratios
because . Models gblsr and
gblthin
Initial magnetic field strength
The independence of the saturated state on the initial field strength has been demonstrated by Sano et al. (2004) (zeronet flux, unstratified shearingboxes), Guan & Gammie (2011) (stratified shearingboxes), and Hawley et al. (2011) (global stratified disks). In all cases, the dissipation, or expulsion, of the initial magnetic field configuration disconnects its influence from the QSS.
Initial perturbation to the disk
The growth rate of the nonaxisymmetric MRI depends on the wavenumber of the initial perturbation (Balbus & Hawley 1992). Faster magnetic field growth arises for higher wavenumbers provided they are subcritical. Therefore, the growth rate of the stresses in the disk will be higher for higher wavenumber perturbations. Parkin & Bicknell (2013) examined this in the context of global disks, finding that irrespective of the initially excited MRI mode, the onset of nonlinear (turbulent) motions in the disk erases the initial perturbation and leads to a statistically similar saturated state.
3.3 Astrophysical implications
Simulations gbllr, gblsr, and gblhr converge at which should be compared with values of commonly derived from relaxation times in postoutburst cataclysmic variables (CVs)(Smak 1999; King et al. 2007; Kotko & Lasota 2012). Although a discrepancy exists, we note that our converged values are consistent with those for isolated AGN disks (Starling et al. 2004). Attaining higher values for may, therefore, require the influence of a companion star to be included in simulations. Alternatively, largescale (net vertical flux) magnetic fields and/or large magnetic Prandtl numbers have been shown to yield larger stresses (Lesur & Longaretti 2007; Fromang et al. 2013; Bai & Stone 2013; Lesur et al. 2013).
4 Control volume analysis
In order to understand the global characteristics of our model accretion disks, we have preformed a control volume analysis of the magnetic energy budget. This involves evaluating terms in the magnetic energy equation integrated over a specific volume, and for this purpose we choose the disk body (defined in § 2.3). The boundaries of this control volume are open in the radial and vertical directions, and periodic in the azimuthal direction.
Previous shearing box simulations show a common characteristic of a consistent power on large scales when convergence is achieved (Simon et al. 2009; Davis et al. 2010). Thus the large scale power in the turbulent spectrum is important when assessing the convergence properties of accretion disk simulations. In turbulent gases (and in our simulations) most of the energy is contained in the largest scales so that our volumeintegrated approach to the energy budget is useful in understanding the characteristics of the low wave number part of the spectrum.
Another important aspect of the control volume analysis is that it provides a guide for constraining the various terms describing production, advection, volumetric changes and dissipation in analytic modeles of accretion disks (e.g. Balbus & Papaloizou 1999; Kuncic & Bicknell 2004). For example, this analysis informs us whether the vertical advection of turbulent energy is important compared to the production of turbulence in magnetised disks.



4.1 Magnetic energy evolution
To set the scene, we first examine the control volume averaged magnetic energy, (Fig. 6). The curves follow the general morphology of a rapid rise in during the initial transient phase, followed by a similarly rapid fall in magnetic energy which gradually flattens out as the quasisteady state () is reached. For gblsr, for example, the quasisteady state is reached after . Subsequently, there is a slow, but steady, decrease in magnetic energy.
We begin with the magnetic field induction equation, to which we add a term for numerical resistivity, , such that Eq (4) now reads,
(24) 
The motivation for introducing will become clear in the remainder of the paper. For now we merely note that the truncated order of accuracy of numerical finite volume codes (such as the PLUTO code used in this investigation) brings with it a truncation error which we interpret as a numerical resistivity and which we model with the additional Ohmic term in Eq (24). Taking the scalar product of with Eq (24) and rearranging terms gives,
(25) 
where , the fluid shear tensor,
(26) 
and the Maxwell stress tensor is given by Eq (15), and a subscript comma denotes partial differentiation. Next we expand , where is the perturbed velocity field in the rotating frame and,
(27) 
is the azimuthally averaged rotational velocity. This step allows the respective contributions to the terms in Eq (25) from the mean background disk rotation and the perturbed velocity field (in the rotating frame) to be inspected. Substituting Eq (27) into Eq (25) and integrating over a control volume with bounding surface , and using the relation,
(28) 
to separate shear and expansion terms (where a subscript semicolon indicates a covariant derivative) one arrives at,
(29) 
where,
(30)  
(31)  
(32)  
(33)  
(34)  
(35) 
and where the current density, , and the Maxwell stress tensor, , is given by (15). All terms featuring in Eqs (29)(35) are exact and can be explicitly calculated from the simulation data. The numerical resistivity, , is estimated by solving Eq (29) for and then solving Eq (35) for . To maintain consistency with the thirdorder spatial reconstruction used in the simulations, we compute terms appearing in Eqs (29)(35) to thirdorder accuracy using reconstruction via the primitive function (see Colella & Woodward 1984; Laney 1998, for further details).
Before proceeding to the results of the control volume analysis, a brief description of the terms and their respective meaning is worthwhile. The volume integrated rate of change of magnetic energy is given by . and are the production of magnetic energy by the shear in the mean disk rotation and the turbulent velocity field, respectively. corresponds to changes in magnetic energy due to expansion in the gas. is a surface term for the advection of magnetic energy in/out of the control volume by the turbulent velocity field. There are contributions to the surface integrals from the radial and direction  periodic boundaries in the azimuthal direction lead to a cancellation, and thus no contribution from those surfaces. (Note that the term vanishes because of the periodic boundary conditions in the azimuthal direction, therefore the mean disk rotation does not advect magnetic energy in/out of the control volume.) Finally, corresponds to numerical dissipation. We note that the value of is exact, as it is merely the remainder required to balance the magnetic energy equation (Eq 29). However, our determination of from is not exact in view of our assumed Ohmic form for the numerical resistive term in Eq (24). Nevertheless, we consider this estimate to be indicative of the actual numerical resistivity.
In Fig. 7 we plot the results of applying the control volume analysis to model gblsr. One immediately notices that magnetic energy production is dominated by (see also discussion in Kuncic & Bicknell 2004) and removal is predominantly via numerical dissipation, . There is nonnegligible magnetic energy removal by divergence in the velocity field (). Examining the directional contributions to this term, one finds roughly equal magnitudes for the , , and components. However, the poloidal contributions () are expansions, which remove magnetic energy, whereas the azimuthal contribution is compressive, thus being a source of magnetic energy. Flow divergence impacting on magnetic field evolution has also been observed in stratified shearingbox simulations by Johansen et al. (2009). The turbulent velocity field does not contribute greatly to magnetic energy production, as is demonstrated by the comparably small values for the curve. A negligible amount of magnetic energy appears to be advected out of the volume in the radial and vertical directions as shown by the curve for , consistent with stable magnetically buoyant (Parker) modes within (Shi et al. 2010). Therefore, although “butterfly” diagrams indicate quasiperiodic vertical magnetic field expulsion (e.g. Gressel 2010), it would seem that a much greater amount of energy is dissipated within the disk body. The rate of change of magnetic energy is relatively small compared to magnetic energy production by and dissipation by . Computing a timeaveraged value between orbits 1232, we find . Therefore, although exhibits quasisteady behaviour, is continually declining, but at a constantly decreasing rate.
Examining the dissipation term, , in more detail, one
finds that the first term in square brackets on the RHS of
Eq (35) is considerably smaller than the second. This
shows that dissipation is primarily powered by the current
density

The formulation of the magnetic energy equation used in Eq (29) allows one to distinguish the contributions from shearing and expansions in the disk. However, with a view to understanding the influence of the boundary conditions for the control volume on the magnetic energy, and on the power injected by Maxwell stresses and Lorentz forces, an alternative formulation may be used. To this end we recast Eq (29) as:
(36) 
where,
(37)  
(38)  
(39)  
(40)  
(41)  
(42) 
and where , , and are given by Eqs (30), (34), and (35), respectively. The rates of work done within the control volume by the Lorentz force, , in combination with the mean disk rotation, , and the turbulent velocity field (in the rotating frame), , are given by and , respectively. Similarly, the rates of work done on the surfaces of the control volume by combinations of the Maxwell stresses and and are, respectively, given by and .
In Fig. 9 we show the result of applying Eq (36) to model gblsr. Magnetic energy production, which was shown to be predominantly due to in Fig. 7, can now be attributed to the rates of work done by the mean disk rotation in combination with Maxwell stresses applied to the boundaries of the volume, , and Lorentz force acting within the volume, . In contrast, the turbulent velocity field acts to remove energy from the control volume, as shown by the terms and . Examining , which involves an integration over the surfaces of the control volume, one finds the magnitude of the radial surface terms to be much greater than from the vertical surfaces and, in particular, the inner radial surface dominates. Therefore, the rate of magnetic energy production is to a large extent due to the difference between the rates of work done on the radial surfaces of the control volume by Maxwell stresses, and by Lorentz forces within the volume.


Model  

gbllr  1231  770  
gblsr  1231  1328  
gblhr  1231  2139  
gbllrla  1831  1328  
gblhrla  1231  2081  
gblthin  1231  1695 
4.2 Numerical resistivity
Computing the numerical resistivity, , provides insight into the intrinsic dissipation arising from the simulation method, embodying the truncated order of accuracy present in commonly used numerical schemes. For example, in our present investigation we use thirdorder accurate spatial reconstruction and secondorder accurate timestepping. Model gblsr returns and , whereas, on the basis of the conclusions drawn by Fleming et al. (2000), Oishi & Mac Low (2011), and Flock et al. (2012b), sustained turbulence should not be observed for . This disagreement does not appear to be related to our approximation of a constant throughout the control volume, as tests computed for annuli with radial range and reveal variations in of only . However, it may be due to our assumption that numerical resistivity behaves like an Ohmic resistivity. Previous estimates of numerical resistivity (Fromang & Papaloizou 2007; Simon et al. 2009) adopt a Fourier analysis of the dissipation term whereby is derived from the high wavenumber end of the spectrum. These analyses reveal that numerical dissipation deviates from Ohmic at low wavenumbers. Our volume averaged values for provide an estimate which is biased towards the largescales, and thus may be higher than values at the small scales (i.e. the turbulent dissipation scale). We do note, however, that the thirdorder accurate spatial reconstruction used in our simulations may allow sustained turbulence at lower than the secondorder accuracy used by Flock et al. (2012b). Furthermore, shearingbox boundary conditions suppress terms in the magnetic energy equation that can supply/sustain large scale magnetic fields  the importance of global disk boundary conditions to magnetic field generation is discussed in more detail in § 5.3. Therefore, the largescale dynamo apparent in stratified global disks (Fig. 11  see also Arlt & Rüdiger 2001; Fromang & Nelson 2006; O’Neill et al. 2011) may operate effectively at low (Brandenburg 2009; Käpylä & Korpi 2011), meaning that global disks could exhibit sustained turbulence at lower than in a shearingbox (Fleming et al. 2000; Oishi & Mac Low 2011). For further discussion of numerical resistivity see Hirose et al. (2006) and Hawley et al. (2011).
In summary, there are two main reasons for the difference by a factor for the critical Reynolds number for the maintenance of turbulence in global disk models compared to the work of Flock et al. (2012b): (1) Differing mathematical approaches to the estimation of the resistivity  Flock et al. (2012b) include an Ohmic resistive term specifically in their simulations, whereas we estimate it using an Ohmic model for the numerical resistivity, (2) The Flock et al. simulations are spatially second order accurate whereas our simulations (and analysis) are spatially third order accurate.
Examining the results for models gbllr, gblsr, and gblhr in
Fig. 10 and Table 3, there is the
consistent trend that as the resolution is increased (and the cell
aspect ratio is kept fixed) the value of
decreases. For example, between models gbllr and gblsr the
resolution has been increased by a factor of 1.5 resulting in a
decrease in by a factor of 1.8. Based on the results
for these three simulations we find , where is the number of cells in the radial
direction, leading to an estimated resolution requirement of cells to achieve a magnetic Reynolds number, . This poses a significant computational
challenge
Somewhat surprisingly, model gbllrla displays a lower value for , and thus higher , than gbllr despite the former having a larger cell aspect ratio. Considering the lower level of turbulent activity in model gbllrla compared to gbllr (see Figs. 1 and 8), this shows that numerical resistivity scales with the turbulent motion of the magnetic field, i.e. a larger value of causes a larger net truncation error. In model gbllrla, turbulent activity wanes for , and simultaneously the value of dips (Fig. 10).
5 Boundary conditions and convergence
The question of convergence in simulation studies of magnetized accretion disk turbulence has been longstanding (e.g. Hawley et al. 1995; Stone et al. 1996; Sano et al. 2004; Fromang & Papaloizou 2007; Simon et al. 2009; Guan et al. 2009; Johansen et al. 2009; Davis et al. 2010; Hawley et al. 2011; Sorathia et al. 2012). It is clear that the development, or initial presence, of largescale magnetic field components is a vital ingredient in enabling convergence with increasing simulation resolution  see the discussion in the second paragraph of § 4. When present, largescale magnetic fields can replenish low wavenumber magnetic energy. This is a prerequisite for convergence since otherwise the reservoir of magnetic energy on the largest scales is drained by the turbulent cascading of magnetic energy to smaller scales. We show in this section that the simulation boundary conditions dictate whether large scale mean fields can grow and thereby promote convergence. In doing so, we consider three different classes of simulation: (1) Unstratified shearing boxes; (2) Stratified shearing boxes; (3) Global, stratified disks and (4) Global, unstratified disks. In the simplest case  the unstratified shearingbox with periodic boundary conditions  mean radial and vertical fields cannot readily evolve. When stratification is introduced the associated interface between the disk body and corona, relaxes the constraint on mean radial field growth such that an dynamo can operate effectively. In global models, mean fields grow relatively quickly, enabling largescale dynamo activity and magnetic energy replenishment. A key result of this analysis is that lower simulation resolution is required in stratified global models compared to shearingboxes because the largescale radial gradients enabled by open radial boundaries permit a larger magnitude contribution to the creation of magnetic energy from the Lorentz force terms. Thus, convergence is attained at lower resolutions in global models than in shearingboxes.
5.1 The magnetic energy balance in accretion disk turbulence
Irrespective of the specific setup (unstratified/stratified shearing box, global simulation) numerical simulations of magnetorotational turbulence exhibit common features. During the early phases of simulation evolution, the MRI develops and the subsequent magnetic field amplification causes a sharp rise in the magnetic energy, . Magnetic energy builtup during the initial transient growth phase supports optimal MRI growth and turbulent driving, which inturn dissipates magnetic energy via the resistivity and the current density, . Magnetic energy subsides, and a state is approached where magnetic field production and turbulent dissipation come into balance. This latter stage is the QSS.
Informed by the analysis in § 4, we write the steadystate magnetic energy evolution equation as (see Eq (36):
(43) 
where the separate Lorentz and magnetic stress terms for rotational and turbulent contributions are combined in the following terms:
(44)  
(45) 
where is the total (rotational plus turbulent) velocity, is the Maxwell stress tensor and the Lorentz force. From § 4.1 we have , so that Eq (43) for the magnetic energy balance now reads:
(46) 
Eq (46) states that the achievement of a QSS requires the rate of work done on the surfaces of the control volume by Maxwell stresses, and by Lorentz forces within the volume, to be balanced by dissipation. Magnetic field saturation and the quasisteady state are solutions to Eq (46).
5.2 Integrated form of the induction equation
In our analysis below, the interplay between the induction equation and the boundary conditions also plays an important role. We begin with the induction equation,
(47) 
where a numerical resistive term is included for consistency with the magnetic energy balance Eq (46). Integrating Eq (47) over a control volume, , with bounding surface , we have:
(48) 
where is the element of surface area and we have assumed to be approximately spatially constant.
The surface bounding the control volume as well as the boundary conditions on take several different forms depending upon the simulation  stratified/unstratified shearing box, stratified/unstratified global disk. However, in general we can use the coordinate convention introduced for shearing boxes (Hawley et al. 1995), adapting the mathematical analysis in each of the different cases. Therefore, the coordinates , and correspond to the radial, azimuthal, and vertical directions in the control volume, respectively. The surface integral then involves three separate integrals over the , and faces, which we denote by , and respectively. In our representation of the integrated induction equation we introduce the resistive flux,
(49) 
This represents a diffusion of magnetic field, resulting from resistive effects through the bounding surface . The integrated induction equation becomes for each coordinate:
(50)  
(51)  
(52)  
The surface integrals in Eqs (50)  (52) show the influence of the velocity and magnetic field values at the boundaries on the volume integrated field within the control volume.
5.3 Dependence of convergence on boundary conditions and magnetic field configuration
In the following sections we utilise our description of the magnetic energy balance combined with inferences from the induction equation to describe how the convergence properties of simulations with different numerical setups can be readily understood in terms of the respective boundary conditions and net magnetic field configuration.
Unstratified shearingbox
In this case the model is a periodic box with background shear applied via source terms in the momentum equation. The shearingbox method is used to represent a small patch of an accretion disk in a Cartesian coordinate system such that , , and correspond to the radial, azimuthal, and vertical directions, respectively; the corresponding lengths of each side of the box are , and (see Hawley et al. 1995, for further details). The boundaries of the control volume in this setup are the boundaries of the computational domain. For an unstratified shearingbox, the following shearingperiodic boundary conditions are applied in the radial (), azimuthal () and vertical () directions for all dynamical variables except the azimuthal velocity. Let be the shear parameter (=3/2 for a Keplerian disk), then the , and boundary conditions are:
(53)  
(54)  
(55) 
The exception is the azimuthal velocity, which satisfies the above  and  boundary conditions but whose boundary condition is:
(56) 
Applying these boundary conditions to Eq (46), and noting that , we have,
(57) 
where, as noted, refers to the inner radial boundary and is the corresponding element of surface area. Furthermore, noting that structures are typically elongated in the direction, then or , and retaining the largest terms (those linear in or ), one finds,
(58) 
The crucial feature of Eq (58) is that magnetic energy is produced by a combination of the component of the Maxwell stress at the radial boundary and Lorentz forces doing work within the volume; the Lorentz forces depend on the radial and vertical field components as well as the radial and vertical gradient in . However, the contributions from the second and third terms on the LHS of Eq (58) are negligible on large scales if there is zeronet radial and vertical magnetic field, and/or no radial or vertical gradient in