Global MHD disks I

Global simulations of magnetorotational turbulence I: convergence and the quasi-steady state


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 time-and-volume-averaged stress-to-gas-pressure ratio, , at a value of , is found for a model with radial, vertical, and azimuthal resolution of 12-51, 27, and 12.5 cells per scale-height (the simulation mesh is such that cells per scale-height varies in the radial direction). The gas pressure dependence of the quasi-steady state stress level is also examined using models with different scaleheight-to-radius 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. Re-casting 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 shearing-box 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 large-scale 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 shearing-boxes is explained by this finding.

accretion, accretion disks - magnetohydrodynamics - instabilities - turbulence

1 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 self-sustaining magnetized turbulence driven by the magnetorotational instability (MRI) can play this role (Balbus & Hawley 1998).

Due to the highly non-linear nature of magnetorotational turbulence, numerical simulations have become a common tool in its study. These simulations come in a number of flavours: unstratified shearing-boxes (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 shearing-boxes (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). Shearing-box 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 quasi-steady state (QSS).

There have been mixed results from simulations as to what sets the QSS stress level. The results of unstratified shearing-box 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 zero-net flux, unstratified shearing-box 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 shearing-boxes 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 shearing-box 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 scale-height. The results from the application of a control-volume 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 time-dependent 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:


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 pseudo-Newtonian potential (Paczyńsky & Wiita 1980):


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 ad-hoc cooling term used to keep the scale-height 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 scale-height over time. The form of used is identical to that of Parkin & Bicknell (2013); further details can be found in that paper3.

The PLUTO code was configured to use the five-wave HLLD Riemann solver of Miyoshi & Kusano (2005), piece-wise parabolic reconstruction (PPM - Colella & Woodward 1984), limiting during reconstruction on characteristic variables (e.g. Rider et al. 2007), second-order Runge-Kutta time-stepping, 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 scale-height. For our fiducial model, gbl-sr, 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,


where the pressure, , and the ratio of gas-to-magnetic 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:


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,




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, non-axisymmetric Fourier mode is excited in the poloidal velocities with amplitude , where is the sound speed.

2.3 Diagnostics

A volume-averaged value (denoted by angled brackets) for a variable is computed via,


Similarly, azimuthal averages are denoted by square brackets,


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,


where the Reynolds stress tensor,


and the Maxwell stress tensor,


The largest contribution comes from the component of (Brandenburg et al. 1995; Hawley et al. 1995; Stone et al. 1996),


where we have defined the perturbed flow velocity as4 , with . Normalising by the gas pressure defines the -parameter,


Furthermore, we calculate the component of the Maxwell stress normalised by the magnetic pressure,


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,


where , and is the Alfvén speed, the “quality factor” is given by,


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,


where represents a cell.

Model Resolution
(, , ) ()
gbl-lr 0.1 340,112,128 8.5-36 18 8
gbl-sr 0.1 512,170,196 12-51 27 12.5
gbl-hr 0.1 768,256,256 18-77 37 16
gbl-lr-la 0.1 420,140,70 10.5-45 20 4.5
gbl-hr-la 0.1 768,256,128 18-77 37 8
gbl-thin 0.05 512,170,320 6-25 27 10
Table 1: List of global simulations.

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,


It then follows that the angle-averaged (in Fourier space) amplitude spectrum,


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 power-spectra, each time-averaged over . These 30 realisations are then averaged-over to produce the final power spectrum.

gbl-lr 12-31 0.57 0.27 0.67 131 395 14 12 0.043 0.39
gbl-sr 12-31 0.69 0.43 0.75 128 361 17 14 0.040 0.42
gbl-hr 12-31 0.78 0.56 0.80 123 332 18 15 0.039 0.42
gbl-lr-la 18-31 0.31 0.07 0.25 655 2296 34 32 0.013 0.29
gbl-hr-la 12-31 0.76 0.50 0.61 150 430 18 15 0.035 0.39
gbl-thin 12-31 0.52 0.44 0.75 120 416 15 13 0.044 0.41
Table 2: List of time averaged quantities from the global simulations. (second column) is the time interval over which time averaging was performed.

2.5 Summary of models

In Table 1 we list six simulations aimed at investigating the following points:

  • Convergence with resolution: Models gbl-lr, gbl-sr, and gbl-hr are low, standard, and high resolution variants, respectively, with identical cell aspect ratio and .

  • Importance of azimuthal resolution: Model gbl-lr-la (gbl-hr-la) is identical to gbl-lr (gbl-hr) with the exception of a lower azimuthal resolution (denoted by the affix “-la”).

  • Scale-height dependence: Models gbl-sr and gbl-thin have disk scale-heights of and 0.05, respectively. These models feature an identical number of cells per scale-height in the vertical and azimuthal directions.

Figure 1: The time evolution of (upper) and (lower) in the global models, where time is in units of the orbital period at a radius of , . (For comparison, , therefore roughly 370 inner disk orbits are covered.) Details pertaining to the models are listed in Table 1 and corresponding time averaged results are given in Table 2.
Figure 2: Resolvability (see Eq 21) of the MRI in the global simulations in the (upper), (middle), and -directions (lower).
Figure 3: Angle-averaged energy spectrum calculated from time-averaged simulation data showing density (top), kinetic energy (middle), and magnetic energy (lower). The dotted lines are representative power-law slopes (see § 3). The horizontal axis is in units of .
Figure 4: Angle-averaged energy spectrum for the stress-to-gas-pressure ratio, , calculated from time-averaged simulation data. The horizontal axis is in units of .

3 The quasi-steady 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 steady-state turbulence in Table 2.

Figure 5: 3D volume rendering showing the ratio of magnetic pressure to gas pressure () for models gbl-lr (top), gbl-sr (middle), and gbl-hr (lower). A wedge has been excised from the upper hemisphere of the disk to expose the disk mid-plane.

3.1 Resolution dependence

Convergence: gbl-lr, gbl-sr, and gbl-hr

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 level-off and a quasi-steady state is reached. We find a time-averaged value during the quasi-steady state of for models gbl-lr, gbl-sr, and gbl-hr, 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 shearing-box 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 gbl-sr, 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 zero-net 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 gbl-lr, gbl-sr, and gbl-hr 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 power-law slope suggests a self-similar 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 scale-height (, where ), whereas they are approximately equal on the smallest length scales (). This differs from the stratified shearing-box 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 two-dimensional 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 pseudo-Cartesian 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 turn-over 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 MRI-driven 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 gbl-lr, gbl-sr, and gbl-hr 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 12-51 cells/ in radius, 27 cells/ in the -direction, and 12.5 cells/ in the -direction (model gbl-sr). This is considerably below the 64-128 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 : gbl-lr-la and gbl-hr-la

When the azimuthal field is under-resolved, 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 gbl-lr-la (a lower azimuthal resolution variant of gbl-lr - see Table 1), which we plot in Figs. 1 and 2. Repeating this experiment at higher resolution (models gbl-hr and gbl-hr-la), one finds that even though the azimuthal field is barely-resolved (8 cells/ in the azimuthal direction for gbl-hr-la), 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 turn-over 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 gbl-lr and gbl-lr-la in Fig. 3, the slope in the magnetic energy power spectrum is steeper for gbl-lr-la in the wavenumber range . This steeper slope is readily understood as a consequence of under-resolving the fastest growing MRI modes in the -direction - energy cannot be injected by the MRI if the mode-growth is not resolved. In relation to the discussion in the previous section, this implies that the power-law slope in the magnetic energy power spectrum of model gbl-lr, 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: gbl-thin

Based on a suite of unstratified shearing-box 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 gbl-sr and gbl-thin5, which feature and 0.05, respectively, and an identical number of cells per scale-height in the vertical and azimuthal directions (Table 1), return comparable values of , with a slightly higher value for the latter (Table 2). Based on the two disk aspect ratios we have explored, there is a weak dependence of on gas pressure. This lack of dependence stems from the similarity in values for the disk body plasma-, (Table 2). Essentially, even though gas pressure is higher in gbl-sr compared to gbl-thin, the relative strength of the Maxwell stresses is very similar.

Initial magnetic field strength

The independence of the saturated state on the initial field strength has been demonstrated by Sano et al. (2004) (zero-net flux, unstratified shearing-boxes), Guan & Gammie (2011) (stratified shearing-boxes), 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 non-axisymmetric MRI depends on the wavenumber of the initial perturbation (Balbus & Hawley 1992). Faster magnetic field growth arises for higher wavenumbers provided they are sub-critical. 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 non-linear (turbulent) motions in the disk erases the initial perturbation and leads to a statistically similar saturated state.

3.3 Astrophysical implications

Simulations gbl-lr, gbl-sr, and gbl-hr converge at which should be compared with values of commonly derived from relaxation times in post-outburst 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, large-scale (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 volume-integrated 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.

Figure 6: Volume averaged magnetic energy, , as a function of time.
Figure 7: Comparison of terms pertaining to the control volume analysis for model gbl-sr (see § 4.1). The upper and lower panels show results for time intervals 0-32 and , respectively. Rates of change of energy are plotted in units of - to convert to code units divide the values by . Note the difference in scale between the plots.

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 quasi-steady state () is reached. For gbl-sr, for example, the quasi-steady 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,


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 re-arranging terms gives,


where , the fluid shear tensor,


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,


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,


to separate shear and expansion terms (where a subscript semi-colon indicates a covariant derivative) one arrives at,




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 third-order spatial reconstruction used in the simulations, we compute terms appearing in Eqs (29)-(35) to third-order 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 gbl-sr. 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 non-negligible 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 shearing-box 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 quasi-periodic 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 time-averaged value between orbits 12-32, we find . Therefore, although exhibits quasi-steady 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 density6, . In Fig. 8 we show . There is a striking similarity between the morphology of the curves in this plot with those for (Fig. 6), suggesting an intimate link between the evolution of the magnetic energy, dissipation driven by a turbulent magnetic field, and magnetic energy production (also demonstrated by and in Fig. 7). Comparing time-averaged values for from models gbl-lr, gbl-sr, and gbl-hr, we find very little difference. Therefore, as convergence with resolution is achieved for , the level of dissipation also converges. We elaborate on the above points in § 5 in the context of a unified description for the observed evolution in accretion disk simulations.

Figure 8: Integral of the current density squared, , over the control volume. This term indicates the level of turbulent activity in the magnetic field and is the dominant contributor to .

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 re-cast Eq (29) as:




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 gbl-sr. 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.

Figure 9: Comparing different terms from the one-zone disk body model for model gbl-sr. The top and bottom panels show results over the time intervals 0-32 and , respectively. Rates of change of energy are plotted in units of - to convert to code units divide the values by . Note the difference in scale between the plots.
Figure 10: Numerical resistivity, , computed from the control volume analysis (see § 4.1 and 4.2).
gbl-lr 12-31 770
gbl-sr 12-31 1328
gbl-hr 12-31 2139
gbl-lr-la 18-31 1328
gbl-hr-la 12-31 2081
gbl-thin 12-31 1695
Table 3: Time averaged numerical resistivity, , and magnetic Reynolds number, . is estimated using Eq (29) and . (second column) is the time interval over which time averaging was performed.

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 third-order accurate spatial reconstruction and second-order accurate time-stepping. Model gbl-sr 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 large-scales, and thus may be higher than values at the small scales (i.e. the turbulent dissipation scale). We do note, however, that the third-order accurate spatial reconstruction used in our simulations may allow sustained turbulence at lower than the second-order accuracy used by Flock et al. (2012b). Furthermore, shearing-box 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 large-scale 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 shearing-box (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 gbl-lr, gbl-sr, and gbl-hr 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 gbl-lr and gbl-sr 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 challenge7. Note that our derived scaling for numerical resistivity is identical to that found by Simon et al. (2009) for unstratified net flux shearing-boxes.

Somewhat surprisingly, model gbl-lr-la displays a lower value for , and thus higher , than gbl-lr despite the former having a larger cell aspect ratio. Considering the lower level of turbulent activity in model gbl-lr-la compared to gbl-lr (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 gbl-lr-la, 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 long-standing (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 large-scale 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, large-scale magnetic fields can replenish low wavenumber magnetic energy. This is a pre-requisite 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 shearing-box 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 large-scale 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 shearing-boxes because the large-scale 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 shearing-boxes.

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 built-up during the initial transient growth phase supports optimal MRI growth and turbulent driving, which in-turn 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 steady-state magnetic energy evolution equation as (see Eq (36):


where the separate Lorentz and magnetic stress terms for rotational and turbulent contributions are combined in the following terms:


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:


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 quasi-steady 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,


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:


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,


This represents a diffusion of magnetic field, resulting from resistive effects through the bounding surface . The integrated induction equation becomes for each coordinate:


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 shearing-box

In this case the model is a periodic box with background shear applied via source terms in the momentum equation. The shearing-box 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 shearing-box, the following shearing-periodic 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:


The exception is the azimuthal velocity, which satisfies the above - and - boundary conditions but whose -boundary condition is:


Applying these boundary conditions to Eq (46), and noting that , we have,


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,


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 zero-net radial and vertical magnetic field, and/or no radial or vertical gradient in