Saturation of the magnetorotational instability in the unstratified shearing box with zero net flux: convergence in taller boxes

Abstract

Previous studies of the non-linear regime of the magnetorotational instability in one particular type of shearing box model – unstratified with no net magnetic flux – find that without explicit dissipation (viscosity and resistivity) the saturation amplitude decreases with increasing numerical resolution. We show that this result is strongly dependent on the vertical aspect ratio of the computational domain L_z/L_x. When L_z/L_x ≲ 1, we recover previous results. However, when the vertical domain is extended L_z/L_x ≳ 2.5, we find the saturation level of the stress is greatly increased (giving a ratio of stress to pressure α ≳ 0.1), and moreover the results are independent of numerical resolution. Consistent with previous results, we find that saturation of the magnetorotational (MRI) in this regime is controlled by a cyclic dynamo which generates patches of strong toroidal field that switches sign on scales of L_x in the vertical direction. We speculate that when L_z/L_x ≲ 1, the dynamo is inhibited by the small size of the vertical domain, leading to the puzzling dependence of saturation amplitude on resolution. We show that previous toy models developed to explain the MRI dynamo are consistent with our results, and that the cyclic pattern of toroidal fields observed in stratified shearing box simulations (leading to the so-called butterfly diagram) may also be related. In tall boxes the saturation amplitude is insensitive to whether or not explicit dissipation is included in the calculations, at least for large magnetic Reynolds and Prandtl number. Finally, we show MRI turbulence in tall domains has a smaller critical Pm_c, and an extended lifetime compared to L_z/L_x ≲ 1 boxes.

accretion, accretion discs, dynamo, instabilities, MHD, turbulence, methods: numerical

INTRODUCTION

The local shearing box model has proved to be very useful for studying the non-linear regime of the magnetorotational instability (MRI; Hawley, Gammie & Balbus 1995). It is useful to classify such models into four types, depending on whether or not the domain contains net magnetic flux, and whether or not the vertical component of gravity is included (producing a vertically stratified density profile). Generally speaking, the results of numerical simulations that explore three of these four types of shearing box models can be summarized as follows.

• Unstratified shearing box simulations with net flux show sustained MHD turbulence and numerically converged values of the stress and angular momentum transport. The saturated stress-to-pressure ratio α varies widely (∼10⁻³–0.1) depending on the net magnetic field strength (Hawley et al. 1995; Sano et al. 2004; Guan et al. 2009; Simon, Hawley & Beckwith 2009).

• Vertically stratified shearing box simulations with no net flux also show strong MRI driven turbulence with typical α ∼ 10⁻³–10⁻², as well as dynamo activity that leads to a quasi-periodic pattern of alternating toroidal field (Brandenburg et al. 1995; Stone et al. 1996; Gressel 2010; Simon et al. 2013a). The stress is independent of numerical resolution (Davis, Stone & Pessah 2010; Shi, Krolik & Hirose 2010), although recently Bodo et al. (2014) have claimed at very high resolution this may no longer be the case. Further investigation of this issue is required.

• Vertically stratified shearing box simulations with net flux show a wide range of behaviour depending on the field geometry and strength (Stone et al. 1996; Miller & Stone 2000). For example, sustained turbulence is observed with weak toroidal fields, while powerful outflows that depend on the field strength are produced in the case of net vertical fields, and the disc can show complex interplay between MRI and buoyancy (Parker) instabilities for strong toroidal fields (Suzuki & Inutsuka 2009; Guan & Gammie 2011; Simon, Hawley & Beckwith 2011; Bai & Stone 2013a; Fromang et al. 2013; Lesur, Ferreira & Ogilvie 2013; Simon et al. 2013b).

The final case of an unstratified shearing box with no magnetic flux, introduced by Hawley, Gammie & Balbus (1996) to study the MHD dynamo driven by the MRI, shows intriguing behaviour. In this case simulations find that the saturated level of stress decreases as the numerical resolution is increased (Fromang & Papaloizou 2007; Guan et al. 2009; Simon et al. 2009; Bodo et al. 2011). When physical dissipation (viscosity and resistivity) is included in the model, convergence of the stress with resolution is recovered (Fromang et al. 2007). Numerous works have also shown that the ratio of these two dissipation length-scales seems to control the saturation level of stress driven by MRI turbulence (Fromang et al. 2007; Simon et al. 2009; Fromang 2010), even when net flux is included (Lesur & Longaretti 2007; Simon & Hawley 2009; Meheut et al. 2015).

Since the important work of Fromang & Papaloizou (2007) there has been significant effort to investigate MHD dynamo action in the special case of unstratified no net flux shearing box. Almost all of this work is based on spectral methods, so that explicit dissipation (both viscosity and resistivity) must be included. For example, Lesur & Ogilvie (2008b) found that long-lived cyclic dynamo activity occurs in an incompressible MHD model in which the aspect ratio of the computational domain is (L_x : L_y : L_z) = (1 : 4 : 2), and these authors developed a toy dynamical model to explain the cycles. More recently, the work of Rincon, Ogilvie & Proctor (2007), Rincon et al. (2008) has explored the subcritical dynamo mechanism in more detail, and in particular the role of the magnetic Prandtl number (ratio of viscosity to resistivity) in determining the outcome. Recent work (Riols et al. 2013, 2015) reveals the non-linear dynamics of the MRI dynamo is more complex than expected.

Despite this progress, the question remains why, in compressible MHD with no explicit dissipation, does the saturation amplitude of the MRI depend so strongly on resolution? One hint may come from the extended vertical domain used in the study of MRI dynamos in incompressible MHD, e.g. Lesur & Ogilvie (2008b). Most previous studies in compressible MHD adopt a standard computational domain that spans one thermal scaleheight H in the radial direction, and uses a vertical aspect ratio L_z/L_x = 1 and a toroidal aspect ratio L_y/L_x ≳ 4. In contrast, studies of shear-driven dynamos have found that a large vertical aspect ratio is required for vigorous dynamo action (Yousef et al. 2008), consistent with the results of Lesur & Ogilvie (2008b).

In this work, we present new studies of the saturation of the MRI in compressible MHD in the no-net flux unstratified shearing box, both with and without explicit dissipation, to investigate the role of the size and aspect ratio of the computational domain on the results. Similar to the results reported in incompressible MHD, we find qualitatively different behaviour when the vertical aspect ratio of the domain is large, L_z/L_x ≳ 2.5. In such case vigorous dynamo action produces a strong, cyclic toroidal magnetic field that greatly increases the level of turbulence and stress. Moreover the saturation amplitude of the turbulence is independent of numerical resolution, and independent of whether explicit dissipation is included in the model (at least for large magnetic Reynolds numbers and magnetic Prandtl numbers greater than one).

In reality, the unstratified shearing box with no net flux has very little relevance for real astrophysical discs, because it is impossible for each local patch in a global disc to maintain zero-net-flux for all time. Even if overall the disc has no net flux, large-scale field loops on the scale of the vertical size of the disc will impart smaller local patches with a time-varying net flux. Moreover, in many astrophysical plasmas the dominant non-ideal MHD effect is not simply Ohmic dissipation, but rather ambipolar diffusion and/or the Hall effect (Wardle 1999; Salmeron & Wardle 2003). Significant progress has been made recently on investigating the role of ambipolar diffusion (Hawley & Stone 1998; Bai & Stone 2011, 2013b; Simon et al. 2013a,b) and the Hall effect (Sano & Stone 2002; Kunz & Lesur 2013; Lesur, Kunz & Fromang 2014; Bai 2014, 2015) on the saturation of the MRI in all four kinds of shearing box models described earlier, and the results differ significantly from models which include only Ohmic resistivity. Nevertheless, it is of interest to understand dynamo action in the unstratified no net flux shearing box (with and without resistivity) if only because it represents such a simple well-posed model. In addition, we show in this paper that the cyclic dynamo action observed in stratified shearing box models (resulting in the so-called butterfly diagram of toroidal field) may be related to that observed in more realistic unstratified domains.

The structure of the paper is as follows: we first describe the equations we solve and the numerical methods adopted in Section 2; the main results are then presented in Section 3; a short discussion of the dynamo action observed in our results follows in Section 4; and finally in Section 5, we summarize and conclude.

METHODS

Equations solved and code description

We solve the compressible MHD equations adopting the ‘shearing box’ approximation. In a Cartesian reference frame corotating with the disc at fixed orbital frequency |$\Omega \hat{\boldsymbol {z}}$|⁠, the equations solved are as follows:

\begin{eqnarray} \frac{\mathrm{\partial} \rho }{\mathrm{\partial} t} + \nabla \cdot (\rho \boldsymbol {v}) = 0 \,, \end{eqnarray}

(1)

\begin{eqnarray} \frac{\mathrm{\partial} \rho \boldsymbol {v}}{\mathrm{\partial} t} + \nabla \cdot \left(\rho \boldsymbol {v}\boldsymbol {v} + \boldsymbol {T} \right) = - 2\rho \Omega \hat{\boldsymbol {z}} \times \boldsymbol {v} +2q\rho \Omega ^2 x \hat{\boldsymbol {x}} \,, \end{eqnarray}

(2)

\begin{eqnarray} \frac{\mathrm{\partial} \boldsymbol {B}}{\mathrm{\partial} t} - \nabla \times \left(\boldsymbol {v}\times \boldsymbol {B} \, - \,\eta \nabla \times \boldsymbol {B}\right) = 0 \,, \end{eqnarray}

(3)

where |$\hat{\boldsymbol {x}}$| refers to the radial direction, ρ is the mass density, |$\boldsymbol {v}$| is the velocity, q = 3/2 is the Keplerian shear parameter, |$\boldsymbol {B}$| is the magnetic strength, and η is the Ohmic resistivity. The total stress tensor |$\boldsymbol {T}$| is defined as

\begin{equation} \boldsymbol {T} = \left(P+\frac{\boldsymbol {B}\cdot \boldsymbol {B}}{8\pi } \right)\boldsymbol {I} - \frac{\boldsymbol {B}\boldsymbol {B}}{4\pi } \, - \,\Pi \,, \end{equation}

(4)

where |$\boldsymbol {I}$| is the identity tensor, |$P=\rho c_{\rm {s}}^2$| is the gas pressure, and c_s is the isothermal sound speed. The viscous stress Π can be expanded as

\begin{equation} \Pi _{ij} = \rho \nu \left(\frac{\mathrm{\partial} v_i}{\mathrm{\partial} x_j}+\frac{\mathrm{\partial} v_j}{\mathrm{\partial} x_j} -\frac{2}{3}\delta _{ij}\nabla \cdot \boldsymbol {v}\right) \,, \end{equation}

(5)

where ν is the kinematic viscosity. Note that in equation (2), the vertical tidal acceleration |$\Omega ^2 z\hat{\boldsymbol {z}}$| is omitted, that is we are studying the unstratified shearing box. We include explicit viscosity and resistivity terms in some of our calculations, as discussed in Section 3.5. However, in order to investigate the convergence of stress with numerical resolution without dissipation, most of simulations are ideal MHD.

We use the athena (Stone et al. 2008; Stone & Gardiner 2010) MHD code for our numerical simulations. We adopt the CTU integrator, third-order piecewise parabolic reconstruction with characteristic tracing in the primitive variables, and the Roe Riemann solver as the basic algorithms. The orbital advection scheme is used to increase the efficiency of calculation (Masset 2000; Johnson, Guan & Gammie 2008). As the shearing box boundary conditions can cause mismatch of the integral of fluxes over the two radial faces due to remap (Gressel & Ziegler 2007), using orbital advection can reduce this mismatch and improve the conservation (Stone & Gardiner 2010). We ensure the conservation of vertical magnetic flux to machine precision by treating the shearing box boundary conditions carefully with a remapping method described in Stone & Gardiner (2010, Section 4).

Initial conditions and run setup

We initialize the disc density with a uniform distribution ρ = ρ₀ = 1 within the box. We set the unit of time Ω⁻¹ = 1 and the unit of length H = c_s/Ω = 1, and therefore the sound speed c_s and initial gas pressure |$P_0=\rho _0 c_{\rm {s}}^2$| are unity as well. The initial magnetic field is |$\boldsymbol {B}=B_0\sin (2\pi x/L_{x})\hat{\boldsymbol {z}}$|⁠, where B₀ is specified through the plasma |$\beta _0\equiv P_0/(B_0^2/8\pi )=100$| for all runs. With this configuration, the net magnetic flux is zero in all three directions. We also ran several models in which the initial magnetic field geometry and strength were varied to check that our results do not depend on the precise form of the initial conditions, provided the net flux is zero. We adopt periodic boundary conditions in the azimuthal (y) and vertical (z) directions, and shearing periodic in the radial (x) direction (Hawley et al. 1995).

We carry out multiple sets of simulations with various box sizes L_x ∈ {0.25, 0.5, 1.0} H, vertical to radial aspect ratios L_z/L_x = 0.5–12 and numerical resolutions (Δx)⁻¹ = (Δy)⁻¹ = (Δz)⁻¹ ∈ {32, 64, 128, 256} H⁻¹. Parameters for all our runs are all listed in Table 1. We do not vary the toroidal aspect ratio L_y/L_x in this work, instead we fix L_y/L_x = 4 for all runs. Provided L_y/L_x ≳ 1, it has been found that different sized domains in the toroidal direction primarily affect the spectrum and amplitude of spiral density waves excited by the turbulence and amplified by shear Heinemann & Papaloizou (2009), however such waves do not dominate the saturation level of stress.

Table 1.

Simulation parameters and results.

Namea	αtotb	αMb	αRb	ΔTc	ΔTavgd	Tcyclee	Δtf	|$\left\langle \!\left\langle \frac{B^2}{8\pi P_0}\right\rangle \!\right\rangle _t^{g}$|	|$\left\langle \!\left\langle \frac{\rho v^2}{2 P_0} \right\rangle \!\right\rangle _t^{g}$|	|$\left\langle \!\left\langle \frac{\overline{B}^2}{8\pi P_0} \right\rangle \!\right\rangle _t^{g}$|	|$\left\langle \!\left\langle \frac{\overline{B}^2_y}{8\pi P_0}\right\rangle \!\right\rangle _t^{g}$|
	(×0.01)	(×0.01)	(×0.01)	(orbits)	(orbits)	(orbits)	(orbits)
x1y4z0.5r32	0.63	0.43	0.2	200	150	|$6.4^{+3.6}_{ }$|	0.3	0.0098	0.0067	0.0006	0.0005
x1y4z0.5r64	0.43	0.31	0.12	200	150	|$7.6^{+0.2}_{ }$|	0.5	0.0066	0.0046	0.0005	0.0002
x1y4z0.5r128	0.19	0.14	0.05	200	150	|$7.8^{ }_{-1.8}$|	0.9+0.1	0.0029	0.0021	0.0002	0.0001
x0.5y2z0.5r32	0.44	0.35	0.09	300	200	|$6.2^{+3.8}_{ }$|	0.3	0.0084	0.0035	0.0014	0.0013
x0.5y2z0.5r64	0.32	0.25	0.07	300	200	|$11.0^{ }_{-6.0}$|	0.5−0.1	0.0055	0.0025	0.0004	0.0004
x0.5y2z0.5r128	0.19	0.15	0.04	300	200	|$9.6^{+1.4}_{ }$|	0.9+0.1	0.0032	0.0015	0.0001	0.0001
x1y4z1r32	1.12	0.85	0.27	300	200	|$5.2^{+4.8}_{ }$|	0.5	0.0190	0.0085	0.0017	0.0016
x1y4z1r64	0.54	0.41	0.13	300	200	|$7.6^{+5.8}_{-0.2}$|	|$1.0^{+0.1}_{-0.1}$|	0.0088	0.0041	0.0004	0.0004
x1y4z1r128	0.30	0.26	0.07	200	100	|$15.4^{ }_{-2.4}$|	|$2.1^{+0.1}_{-0.1}$|	0.0050	0.0023	0.0001	0.0001
x0.5y2z1r32	0.89	0.73	0.16	300	200	|$6.6^{+0.4}_{-0.2}$|	0.6+0.1	0.0176	0.0073	0.0042	0.0041
x0.5y2z1r64	0.57	0.47	0.10	300	200	|$8.4^{+0.2}_{-1.0}$|	|$0.9^{+0.1}_{-0.1}$|	0.0104	0.0047	0.0013	0.0013
x0.5y2z1r128	0.29	0.24	0.05	250	150	|$12.6^{+4.2}_{ }$|	|$1.6^{+0.3}_{-0.2}$|	0.0053	0.0024	0.0004	0.0004
x1y4z2r32	1.81	1.44	0.37	300	200	|$9.0^{+2.8}_{}$|	1.1+0.1	0.0333	0.0144	0.0049	0.0047
x1y4z2r64	0.86	0.68	0.17	300	200	|$12.6^{+4.0}_{-2.2}$|	|$1.8^{+0.2}_{-0.1}$|	0.0152	0.0068	0.0010	0.0009
x1y4z2r128	0.39	0.31	0.08	300	150	|$13.0^{+6.0}_{-0.4}$|	|$1.6^{+0.1}_{-0.1}$|	0.0068	0.0031	0.0002	0.0002
x1y4z2.5r32	2.27	1.81	0.46	300	200	|$10.6^{+2.8}_{-0.2}$|	1.4+0.1	0.0422	0.0181	0.0067	0.0065
x1y4z2.5r64	1.30	1.05	0.25	250	150	|$12.6^{+0.8}_{-0.4}$|	1.9+0.1	0.0239	0.0105	0.0022	0.0021
x1y4z2.5r128	1.31	1.08	0.22	200	100	|$20.0^{ }_{-3.4}$|	|$2.0^{+0.1}_{-0.5}$|	0.0202	0.0087	0.0020	0.0020
x1y4z3r32	3.11	2.51	0.60	300	200	|$13.4^{+2.0}_{-1.6}$|	1.6+0.2	0.0618	0.0251	0.0146	0.0141
x1y4z3r64	6.09	5.08	1.01	270	170	|$27.4^{ }_{-2.4}$|	|$2.0^{+0.3}_{-0.3}$|	0.1585	0.0508	0.0628	0.0614
x1y4z3r128	9.13	7.79	1.34	150	100	|$33.2^{ }_{ }$|	|$2.5^{+0.3}_{-0.4}$|	0.2496	0.0739	0.1048	0.1028
x1y4z3.5r32	3.50	2.83	0.68	300	200	|$25.0^{+6.1}_{-4.4}$|	|$1.9^{+0.1}_{-0.1}$|	0.0753	0.0281	0.0234	0.0228
x1y4z3.5r64	6.35	5.29	1.05	300	200	|$27.8^{+5.6}_{-2.8}$|	|$2.7^{+0.3}_{-0.3}$|	0.1630	0.0529	0.0626	0.0612
x0.25y1z1r64	0.42	0.35	0.07	200	100	|$19.6^{+7.0}_{}$|	|$2.9^{+0.2}_{-0.3}$|	0.0118	0.0035	0.0061	0.0060
x0.25y1z1r128	0.49	0.41	0.08	300	200	|$19.0^{+3.2}_{-6.4}$|	|$2.6^{+0.1}_{-0.2}$|	0.0125	0.0041	0.0052	0.0051
x0.25y1z1r256	0.51	0.43	0.08	300	200	|$27.7^{ }_{-5.0}$|	|$3.4^{+0.6}_{-0.7}$|	0.0161	0.0051	0.0065	0.0064
x0.5y2z2r32	1.44	1.17	0.27	300	200	|$17.4^{+4.6}_{-0.8}$|	|$2.4^{+0.2}_{-0.2}$|	0.0370	0.0117	0.0177	0.0173
x0.5y2z2r64	1.83	1.51	0.32	300	200	|$25.0^{+8.4}_{ }$|	|$3.2^{+0.2}_{-0.1}$|	0.0489	0.0151	0.0224	0.0219
x0.5y2z2r128	2.13	1.80	0.33	160	100	|$25.0^{+5.0}_{ }$|	|$3.1^{+0.6}_{-0.3}$|	0.0633	0.0183	0.0289	0.0284
x1y4z4r32	5.16	4.17	0.99	300	200	|$22.8^{+10.6}_{ }$|	|$3.0^{+0.2}_{-0.5}$|	0.1301	0.0417	0.0565	0.0553
x1y4z4r64	7.85	6.49	1.36	300	200	|$30.8^{+2.6}_{-8.6}$|	|$3.6^{+0.3}_{-0.3}$|	0.2562	0.0649	0.1374	0.1353
x1y4z4r128	8.41	7.14	1.26	200	100	|$30.0^{ }_{ }$|	|$3.3^{+0.5}_{-0.7}$|	0.2730	0.0714	0.1331	0.1310
x1y4z4r128pm8	9.25	7.90	1.35	120	100	|$19.4^{ }_{ }$|	|$4.4^{0.7}_{-1.1}$|	0.2716	0.0590	0.0767	0.0752
x1y4z4r32pm4	4.90	3.99	0.91	300	200	|$22.8^{+5.8}_{-0.6}$|	|$2.6^{+0.6}_{-0.7}$|	0.1188	0.0399	0.0387	0.0378
x1y4z4r64pm4	7.11	5.90	1.21	300	200	|$30.0^{ }_{-10.0 }$|	|$3.7^{+0.2}_{-0.6}$|	0.2053	0.0597	0.0989	0.0972
x1y4z4r128pm4	8.53	7.24	1.29	200	100	|$22.6^{+4.8}_{ }$|	|$3.8^{+0.3}_{-0.3}$|	0.2307	0.0604	0.0942	0.0924
x1y4z4r128pm2	5.39	4.51	0.88	300	200	|$32.8^{ }_{ }$|	|$3.6^{+0.5}_{-1.1}$|	0.1421	0.0439	0.0464	0.0455
x1y4z4r128pm1	–	–	–	300	–	–	–	–	–	–	–
x2y8z8r32	20.06	16.08	3.98	300	200	|$34.4^{ }_{-2.4}$|	|$3.9^{+0.6}_{-0.6}$|	0.6387	0.1701	0.3401	0.3351
x2y8z8r64	21.52	17.53	3.99	200	150	|$47.0^{+19.6}_{-11.0}$|	|$3.9^{+0.5}_{-0.7}$|	0.6823	0.1791	0.3451	0.3401
x1y4z6r32	5.93	4.79	1.14	300	250	|$47.6^{+2.4}_{-7.4}$|	|$4.8^{+0.3}_{-0.3}$|	0.1558	0.0477	0.0725	0.0710
x1y4z6r64	7.60	6.30	1.31	280	180	|$37.0^{+3.0}_{-2.2}$|	|$6.9^{+0.5}_{-0.4}$|	0.2390	0.0632	0.1228	0.1207
x0.25y1z2r64	0.42	0.34	0.08	300	200	|$33.4^{+6.6}_{ }$|	|$7.0^{+0.4}_{-0.5}$|	0.0146	0.0034	0.0091	0.0090
x0.25y1z2r128	0.48	0.40	0.08	250	150	|$50.0^{+3.4}_{-1.4}$|	|$9.2^{+1.4}_{-1.9}$|	0.0211	0.0040	0.0151	0.0150
x0.25y1z2r256	0.52	0.44	0.08	200	100	|$59.6^{+19.2}_{ }$|	|$9.6^{+2.2}_{-1.1}$|	0.0531	0.0044	0.0143	0.0141
x0.5y2z4r32	1.53	1.24	0.28	300	200	|$40.0^{ }_{-11.4}$|	|$6.0^{+0.2}_{-0.1}$|	0.0451	0.0124	0.0250	0.0245
x0.5y2z4r64	2.04	1.68	0.36	300	200	|$44.3^{ }_{-11.0}$|	|$8.8^{+1.4}_{-0.4}$|	0.0648	0.0168	0.0359	0.0353
x0.5y2z4r128	2.13	1.79	0.35	200	150	|$61.4^{+5.0}_{-11.4}$|	|$12.0^{0.4}_{-1.8}$|	0.0918	0.0179	0.0590	0.0583
x1y4z8r32	6.10	4.88	1.22	300	200	|$55.0^{+2.2}_{-7.2}$|	|$9.7^{+0.8}_{-2.3}$|	0.2261	0.0488	0.1428	0.1409
x1y4z8r32by	5.97	4.78	1.19	300	200	|$50.0^{+6.2}_{ }$|	|$8.3^{+0.7}_{-1.1}$|	0.2048	0.0478	0.1225	0.1207
x1y4z8r64	7.21	5.92	1.29	200	100	|$50.0^{+3.4}_{ }$|	|$7.6^{+1.1}_{-2.5}$|	0.3833	0.0595	0.1850	0.1827
x1y4z8r128	7.92	6.66	1.26	200	100	|$62.0^{+20.0}_{ }$|	|$9.2^{+3.0}_{-1.3}$|	0.3790	0.0680	0.2491	0.2470
x1y4z10r32	5.68	4.54	1.13	300	200	|$57.2^{+9.6}_{-7.2}$|	|$11.3^{+0.5}_{-0.4}$|	0.1949	0.0454	0.1178	0.1160
x1y4z10r64	7.16	5.91	1.25	300	200	|$80.0^{ }_{-7.0}$|	|$12.5^{+1.6}_{-2.0}$|	0.2444	0.0591	0.1357	0.1338
x1y4z12r32	6.29	5.04	1.25	300	200	|$96.8^{+4.0}_{ }$|	|$15.7^{+2.8}_{-1.0}$|	0.2000	0.0504	0.1142	0.1122
x1y4z12r64	7.00	5.73	1.27	300	200	|$100.0^{ }_{-5.0}$|	|$16.9^{+1.1}_{-3.1}$|	0.3349	0.0574	0.2306	0.2285

Notes.^aName convention: xn ym zk denotes box dimension (L_x,L_y,L_z)=(n,m,k)H; ‘r32’ means typical resolution 32/H, ‘r64’ and ‘r128’ are 64/H and 128/H respectively; ‘pm8’, ‘pm4’,‘pm2’, and ‘pm1’ denote resistive runs with magnetic Prandtl number Pm = 8,4,2, and 1 for fixed Reynolds number Re = 3125; ‘by’ stands for an initial azimuthal magnetic field |$\boldsymbol {B}(t=0)=B_0 \sin (2\pi x/L_{x})\hat{\boldsymbol {y}}$| instead of the vertical configuration |$\boldsymbol {B}(t=0)=B_0\sin (2\pi x/L_{x})\hat{\boldsymbol {z}}$| used in most runs.

^b α_M and α_R are the Maxwell and Reynolds stresses normalized by |$\rho _0 c_{\rm s}^2$|⁠; α_tot is the sum of these two.

^cDuration of the simulation.

^dLast ΔT_avg orbits chosen for time average.

^eThe period of the oscillating magnetic field cycle. Measured based on the power spectra of the largest vertical mode (k = 2π/L_z) of |$\overline{B_y}(t)$|⁠. We separate the time sequence to two parts, and the errors are derived from the range of T_cycle over two different parts.

^fThe time lag between the azimuthal and radial mean field defined in Section 4.1 and used in equation (13) for calculating the phase lag.

^gThe total magnetic energy 〈〈B²/8πP₀〉〉_t, total kinetic energy 〈〈ρv²/2P₀〉〉_t, total mean field energy |$\langle \langle \overline{B}^2/8\pi P_0\rangle \rangle _t$|⁠, and azimuthal mean field energy |$\langle \langle \overline{B}_y^2/8\pi P_0\rangle \rangle _t$|⁠.

Diagnostics

In order to facilitate analysis and obtain statistical properties of our simulations, we define a few ways to average physical variables. We first define a volume (box) average:

\begin{equation} \langle X \rangle \equiv \frac{\int \,X \,\mathrm{d}x \mathrm{d}y \mathrm{d}z}{\int \mathrm{d}x \mathrm{d}y \mathrm{d}z} \,. \end{equation}

(6)

We also define a time average

\begin{equation} \langle X \rangle _t \equiv \frac{\int \,X \,\mathrm{d}t}{\int \mathrm{d}t} \,. \end{equation}

(7)

The time average is generally applied over ∼200 orbits to eliminate chaotic fluctuations and achieve meaningful statistics (Winters, Balbus & Hawley 2003). Using these averages, we write the time-averaged Maxwell and Reynolds stresses as

\begin{equation} \alpha _{\rm M}\equiv \langle \langle -\frac{B_{x}B_{y}}{4\pi P_0}\rangle \rangle _{t} \,, \alpha _{\rm R}\equiv \langle \langle \frac{\rho v_{x} \delta v_{y}}{P_0}\rangle \rangle _{t} \end{equation}

(8)

respectively. In Table 1 we also list the total internal stress, α_tot ≡ α_R + α_M. We use these stresses as a measure of the saturated level of MRI turbulence throughout this paper.

As some of our simulations develop large-scale coherent field structures within the box, we therefore introduce a horizontal average:

\begin{equation} \overline{X} \equiv \frac{\int \,X\, \mathrm{d}x \mathrm{d}y}{\int \mathrm{d}x \mathrm{d}y} \,. \end{equation}

(9)

This average is useful for decomposing the total magnetic field |$\boldsymbol {B}$| into mean and fluctuating parts. i.e.

\begin{equation} \boldsymbol {B} = \overline{\boldsymbol {B}} + \boldsymbol {b} \,, \end{equation}

(10)

where |$\boldsymbol {b}$| denotes the turbulent (small-scale) fluctuating field.

RESULTS

Saturation in tall versus short boxes

We begin by studying saturation in the standard computational domain, (L_x, L_y, L_z) = (1, 4, 1) H. We begin with a numerical resolution of 32/H. We refer to this simulation as x1y4z1r32, where our naming convention gives both the box dimensions (in H) and resolution. We then double this resolution twice (simulations x1y4z1r64 and x1y4z1r128) to reproduce previous results. The top panels of Fig. 1 show 3D snapshots of the magnetic energy and density in these runs. As the resolution improves from left to right, the magnetic energy is reduced, and turbulence becomes weaker. At the highest resolution (x1y4z1r128 in the right-hand panel) the density field is dominated by k_z = 0 shearing waves (Heinemann & Papaloizou 2009), and there is little indication of turbulence.

Figure 1.

Snapshots (at t = 150 orbits) from standard box (top row, L_x = L_y/4 = L_z = 1 H) and tall box (bottom row, L_x = L_y/4 = L_z/4 = 1 H) simulations. Shown in each case are magnetic energy (far half) and density (near half) distributions. From left to right, the resolution is 32/H, 64/H, and 128/H. Turbulence becomes weaker in the standard box case as resolution increases; while it maintains constant amplitude in the tall box runs. Note the colour scale of the magnetic energy is logarithmic, and the colour bars in the bottom and top rows are different.

We then repeat these runs in a taller box, (L_x, L_y, L_z) = (1, 4, 4) H, varying the resolution from 32/H for run x1y4z4r32, to 64/H for run x1y4z4r64, and to 128/H for run x1y4z4r128. Snapshots of magnetic energy and density for these runs are shown in the bottom row of Fig. 1. In contrast to the standard box case, we see no signs of weakening turbulence as resolution improves from left to right. The typical magnetic energy rises slightly from x1y4z4r32 to x1y4z4r64, and stays roughly at the same level as the resolution doubles again in x1y4z4r128. The typical density fluctuation follows the similar trend, and are dominated by small-scale turbulent fluctuations. At the same time, both the density and magnetic field show a large amplitude vertical mode with wavelength about half the vertical size of the box.

To quantify the resolution effects, we measure the volume-averaged total stress α_tot (the sum of the Maxwell and Reynolds stresses normalized with thermal pressure defined in Section 2.3) and the results are shown in Fig. 2. In the standard box runs (the top-left panel), after a short transient growth (∼30 orbits), the total stresses reaches a quasi-steady state that last hundreds of orbits through the end of the simulations. The time-averaged stresses are therefore well defined and are measured over the last 100–200 orbits (see Table 1 for the exact number used for averages). Comparing the saturation levels at different resolutions, we confirm the stress decreases as the resolution improves, roughly by a factor of two every time the resolution is doubled. In particular, α_tot drops from 10⁻³ for x1y4z1r32 to 5 × 10⁻⁴ when the resolution is doubled, and drops by another factor of ∼2 when 128/H is used. This behaviour is identical to that previously reported (e.g. Fromang & Papaloizou 2007; Simon et al. 2009; Guan et al. 2009; Bodo et al. 2011).

Figure 2.

The volume-averaged stresses for a standard box ((L_x, L_y, L_z) = (1, 4, 1) H, left-hand panel) and a tall box ((L_x, L_y, L_z) = (1, 4, 4) H, right) with various resolutions: 32/H(black solid),64/H(red dotted) and 128/H(green dashed). In contrast to the small box, the vertically elongated box achieves good convergence.

Remarkably, the decrease in stress with resolution that appears in the standard box is not reproduced in the taller box simulations. In the bottom-left panel of Fig. 2, we show that the volume-averaged stress does in fact converge with numerical resolution for runs x1y4z4r32 through x1y4z4r128. In these cases, α_tot ≃ 0.08 is reached once the resolution exceeds 32/H.

Fig. 2 also demonstrates that the amplitude of the converged stress in the tall box is significantly greater than those in the standard domain. For example, the time-averaged α_tot ≃ 0.079 in x1y4z4r64, almost 15 times larger than α_tot ≃ 0.005 in x1y4z1r64, while α_tot ≃ 0.08 in x1y4z4r128, some ∼30 times greater than that of x1y4z1r128. The magnetic energy in simulations performed in the tall box is also much greater than in the standard box runs; typically 〈B²/8πP₀〉 ≃ 0.27 which is a factor of ≳50 time bigger (this is also evident in Fig. 1 where different colour scales must be used for the two cases). The large values of the stress achieved in taller boxes is also of interest in that they are close to the values inferred from observations suggesting that in fully ionized accretion disc (e.g. dwarf novae in outburst) α_tot ∼ 0.1–0.4, alleviating some of the concerns raised by King, Pringle & Livio (2007) regarding the discrepancy between values suggested by observations and measured in simulations. (But see Hirose et al. (2014), in which they found convection-enhanced α above 0.1 when the disc is near the hydrogen ionization transition using a stratified shearing box.)

Further important insights into the differences between the saturated state in standard and tall boxes are revealed by the Fourier power spectrum. The right column of Fig. 3 shows the magnetic energy density power spectra calculated over either spherical shells of constant k (⁠|$k B_k^2$|⁠) or along a constant k_z plane (⁠|$k_z B_{k_z}^2$|⁠) for both the standard (left) and tall box (right). The definition of both can be find in Davis et al. (2010, section 2.1). The spectra in the standard box peak at an intermediate wavenumber that increases with the resolution, and with an overall amplitude that decreases with resolution (but with a shape that is unchanged). In stark contrast, however, in the tall box the power spectra peak (or plateau) at the lowest wavenumbers (k ∼ 1–5) at all resolutions. Such spectra are much more reminiscent of turbulence driven at a large outer scale, with an energy cascade to high k independent of numerical resolution. These spectra suggest there is a resolved characteristic outer length-scale in the tall box simulations. We note the spectrum and its dependence on resolutions of those tall box runs are also very different than those using standard short boxes with explicit dissipation (e.g. fig. 3 in Fromang 2010) due to the lack of a large-scale dynamo in the short boxes. However, a large-scale dynamo, as we will show in the next subsection, is present in our tall box simulations.

Figure 3.

The power spectra of magnetic density energy for a standard box ((L_x, L_y, L_z) = (1, 4, 1) H, (left two panels) and a tall box ((L_x, L_y, L_z) = (1, 4, 4) H, right two panels) with various resolutions: 32/H(black solid),64/H(red dotted) and 128/H(green dashed), where |$B_k^2$| is spherical shell integrated value and |$B_{k_z}^2$| is integrated along constant k_z plane. In contrast to the small box case, where the self-similar power spectra decays with resolution, it achieves good convergence at low wavenumber in the tall box.

A large-scale dynamo

The images of magnetic energy and density shown in Fig. 1, and the power spectra in Fig. 3, both indicate that large-amplitude and large-scale vertical structure appears in simulations in tall boxes. In this section, we show this structure is associated with a large-scale dynamo action that is triggered in the vertically extended box. Fig. 4 plots a space–time diagram of the horizontally averaged azimuthal magnetic field in both standard and tall boxes. Patches of strong |$\overline{B}_y$| that extend over vertical regions of size ∼H are evident in run x1y4z4r32, with the sign of the field alternating quasi-periodically over ∼10–20 orbits. The pattern become even more clear in the panels corresponding to runs x1y4z8r32, x1y4z10r32 and x1y4z12r32 (see especially the middle right panel in Fig. 4). Clearly, a large-scale dynamo must be operating that generates such strong ordered toroidal field. The 3D structure of the magnetic field at three different times during one reversal cycle in run x1y4z4r32 is explored further in Fig. 5. Large-scale (k_x = k_y = 0 and k_z = 1) patterns in the azimuthal magnetic field form along the azimuthal direction. Spatially, the field lines flip signs (from red to blue) with height to conserve the zero-net-flux initial condition. Their orientations also reverse from one state of coherent structures (top left at t = 200 orbits) to another (top right at t = 220 orbits). This result is similar to the dynamo cycle reported in Lesur & Ogilvie (2008b) (see their fig. 3). However, we emphasize that no explicit dissipation (viscosity or resistivity) is required to capture the dynamo, nor is it required to see a converged level of stress.

Figure 4.

Space–time diagrams of the horizontally averaged azimuthal field |$\overline{B}_y$| for various runs as labelled. Runs with L_z < 8 H are duplicated 8 (for x1y4z1r32), 4 (for x1y4z2r32), and 2 (for x1y4z4r32) times, so that the images have the same aspect ratio of run x1y4z8r32. The individual black box shows the actual box we simulate. Large and structured azimuthal B-field appears when L_z > 2 H.

Figure 5.

Top row: large-scale magnetic field in x1y4z4r32 run reverses between t = 200 orbits (top left) and 220 orbits (top right); t = 210 at the middle show more tangled field during this transition. Bottom row: snapshots of large-scale field structures in other tall box runs (x1y4z8r32, x1y4z10r32, and x1y4z12r32 from left to right). In all panels, tubes are field lines in the box with starting points randomly distributed; colour bar describes the value of |$B_y/\sqrt{4\pi \rho _0}c_{\rm {s}}$| to show the orientation and strength of magnetic field in azimuth.

The large-scale dynamo maintains a time-averaged mean field |$\langle \langle \overline{B}^2/8\pi P_0\rangle \rangle _t\simeq 0.06$| for x1y4z4r32 and ≃ 0.13 for x1y4z4r64 and x1y4z4r128 runs, which amount to ∼40–50 per cent of the total magnetic energy in the box (see Tables 1 and 2), close to a state of equalized mean |$\overline{\boldsymbol {B}}$| and turbulent (⁠|$\boldsymbol {b}$|⁠) field strength. In contrast, runs in the standard box do not exhibit strong cyclic dynamo behaviour (see top-left panel of Fig. 4), and the mean field is rather weak, |$\langle \langle \overline{B}^2/8\pi P_0\rangle \rangle _t \lesssim 0.1 \langle \langle b^2/8\pi P_0\rangle \rangle _t$| for x1y4z1r32, and it drops further down to ≃2 per cent of the turbulent magnetic energy in run x1y4z1r128.

Table 2.

Time- and volume-averaged magnetic-field-related quantities.

	x1y4z1r128	x1y4z4r128	x1y4z8r128
〈〈B2/8πP0〉〉t	0.0050	0.2730	0.3790
|$\langle \langle B_x^2/8\pi P_0\rangle \rangle _t$|	0.0007	0.0352	0.0357
|$\langle \langle B_y^2/8\pi P_0\rangle \rangle _t$|	0.0041	0.2208	0.3273
|$\langle \langle B_z^2/8\pi P_0\rangle \rangle _t$|	0.0003	0.0170	0.0160
|$\langle \langle \overline{B}^2 /8\pi P_0\rangle \rangle _t$|	0.0001	0.1331	0.2491
|$\langle \langle \overline{B}_x^2/8\pi P_0\rangle \rangle _t$|	0.0	0.0020	0.0022
|$\langle \langle \overline{B}_y^2/8\pi P_0\rangle \rangle _t$|	0.0001	0.1311	0.2470
|$\langle \langle \overline{B}_z^2/8\pi P_0\rangle \rangle _t$|	0.0	0.0	0.0
〈〈b2/8πP0〉〉t	0.0049	0.1399	0.1299
|$\langle \langle {b}_x^2/8\pi P_0\rangle \rangle _t$|	0.0007	0.0332	0.0335
|$\langle \langle {b}_y^2/8\pi P_0\rangle \rangle _t$|	0.0040	0.0897	0.0803
|$\langle \langle {b}_z^2/8\pi P_0\rangle \rangle _t$|	0.0003	0.0170	0.0160
〈〈 − BxBy/4πP0〉〉t	0.0023	0.0714	0.0680
|$\langle \langle - \overline{B}_x\overline{B}_y/4\pi P_0\rangle \rangle _t$|	0.0	0.0115	0.0110
〈〈 − bxby/4πP0〉〉t	0.0023	0.0599	0.0570
〈〈ρv2/2P0〉〉t	0.0028	0.0707	0.0706
|$\langle \langle \rho v_x^2/2P_0 \rangle \rangle _t$|	0.0014	0.0226	0.0226
|$\langle \langle \rho \delta v_y^2/2P_0 \rangle \rangle _t$|	0.0009	0.0365	0.0369
|$\langle \langle \rho v_z^2/2P_0 \rangle \rangle _t$|	0.0005	0.0116	0.0111
〈〈ρvxδvy/P0〉〉t	0.0007	0.0126	0.0131
〈〈(δρ/ρ0)2〉1/2〉t	0.0312	0.1990	0.2243

Why does dynamo action produce a larger value for the stress which does not vary with numerical resolution? The large-scale vertical patches of azimuthal magnetic field produced by the dynamo act locally as a region with net toroidal flux. Thus, each L_z ∼ H patch acts as an unstratified shearing box with net toroidal field. As is already known, shearing boxes with net flux produce saturated stress which is independent of resolution (Guan et al. 2009; Simon & Hawley 2009). As a result of (locally) strong magnetic field, the saturated Maxwell stress in x1y4z4r128 is α_M ≃ 0.0714, ∼31 × greater than that of x1y4z1r128. Of this total, ∼16 per cent is due to the correlated mean field |$-\overline{B}_x\overline{B}_y/4\pi$|⁠, while the majority is still from the correlation between the perturbed field components − b_xb_y/4π (since |$\overline{B}_x$| is still small). In contrast, in the standard-sized box almost all of the stress is associated with the perturbed field.

We find the stress-to-energy ratio, α_mag ≡ 〈〈 − 2B_xB_y〉〉_t/〈〈B²〉〉_t ≃ 0.27 in the tall box, smaller than in the standard box case in which α_mag ≃ 0.46. A diminished α_mag is also observed in simulations with strong net azimuthal flux (run Y8 in Hawley et al. 1995), where the imposed azimuthal mean field resembles a subsection of our box which contains |$\overline{B}_y$| of the same sign.

Varying the aspect ratio L_z/L_x: a parameter survey

We now investigate the effect of varying the aspect ratio using L_z/L_x ∈ {0.5, 1, 2, 2.5, 3, 3.5, 4, 6, 8, 10, 12} with fixed size L_y = 4H and L_x = H in the horizontal dimensions. For each of L_z/L_x, we also vary the resolution to study numerical convergence. We first plot the α_tot in Fig. 6. This figure clearly demonstrates what may be our most important result. The data fall into two groups which show distinctly different behaviour; the two groups are separated by the vertical dashed line at L_z/L_x ≃ 2.5. For aspect ratios ≲ 2.5 (to the left of the dotted line), the stress decays linearly with the vertical size L_z. Moreover, at any given aspect ratio, the stress decreases with increasing numerical resolution as shown by the decreasing amplitude of the black (32/H), red (64/H) and green (128/H) points. Clearly the standard box, with L_z/L_x = 1 falls in this group. In contrast, for aspect ratios L_z/L_x ≳ 2.5 (to the right of the dotted line), the saturated stress associated with MRI turbulence appears independent of the vertical size L_z, approaching α_tot ∼ 0.1. Moreover, numerical convergence in the value of the stress is achieved for aspect ratios L_z/L_x > 4. For example, for L_z/L_x = 8, the stress values are nearly identical for resolutions of 32/H through 128/H. The behaviour of the two groups in this figure simply reflect the emergence of dynamo action in runs with a large vertical extent, which controls the stress. For example, the L_z/L_x = 4, 8, 10, and 12 cases in Figs 4 and 5 show strong cyclic azimuthal magnetic field driven by the underlying dynamo mechanism.

Figure 6.

The sum of Maxwell and Reynolds stress for various aspect ratios (L_z/L_x) and resolutions for L_x = L_y/4 = 1 H boxes. The vertical dotted line separates the ‘diverging’ and ‘converging’ regions: to the left of this line, the stress decays with increasing resolution; the stress converges to ∼0.08 for different L_z/L_x on the right.

For all runs, the Maxwell stress dominates the total turbulent angular momentum transport; the ratio of the Maxwell to the Reynolds stress varies with different L_z/L_x in a similar way (but at reduced amplitude) as the α_tot. In Fig. 7, we find α_M/α_R rises gradually from ∼3 to 4 as L_z/L_x is increased from 0.5 to 2. The ratio then levels off, ranging between 4 and 6. This is greater than the typical values previously reported in the literature with Keplerian shear (∼3–4) (Hawley et al. 1995; Abramowicz, Brandenburg & Lasota 1996; Stone et al. 1996; Hawley, Balbus & Winters 1999; Sano et al. 2004; Pessah, Chan & Psaltis 2006). As the strong coherent field structures tend to eliminate strong velocity fluctuations and therefore reduce the Reynolds stress in those regions (similar to those observed in the corona regions of stratified boxes (Brandenburg et al. 1995; Miller & Stone 2000)), the box averaged α_M/α_R therefore rises above the previously reported values.

Figure 7.

The ratio of the Maxwell stress and the Reynolds stress as a function of L_z/L_x and numerical resolution. The vertical dotted line at L_z/L_x = 2.5 same as Fig. 6. The Maxwell stress always dominates the Reynolds in all runs; the ratio increases by a factor of ∼2 in the small box regime, while stays roughly constant in the tall box cases. The tall box runs have greater α_M/α_R than previous measurements marked as the shaded region.

The ratio between the Maxwell stress and magnetic energy, α_mag (see definition in Section 3.2), is usually adopted as a measure of sufficient resolution to capture MRI turbulence (e.g. α_mag = 0.3–0.4 in Blackman, Penna & Varnière 2008; Hawley, Guan & Krolik 2011). We find consistent results for our standard boxes (L_z/L_x ≲ 2.5), in which a relatively constant α_mag ∼ 0.47 is obtained, in spite of that no convergence is achieved for the stress. Surprisingly, this ratio falls with (L_z/L_x)^−1/2 in our tall box runs where converged stresses are found. As the stress stays roughly constant in those tall boxes, it is the increase of magnetic energy (or 1/β) as (L_z/L_x)^1/2 that drives this scaling.

The dynamo effect in tall boxes alters the underlying magnetic field structure. As discussed in Section 3.1, most of the stress comes from the correlation of the small-scale field, and |$\langle -B_{x}B_{y}\rangle \sim \langle -b_{x}b_{y}\rangle \gg \langle -\overline{B}_{\rm x}\overline{B}_{\rm y}\rangle$| holds true for both groups. However, the magnetic energy in smaller boxes is mostly from the azimuthal component of the turbulent field, i.e. |$\langle B^2\rangle \sim \langle b_{y}^2 \rangle \gg \langle \overline{B}_{\rm y}^2\rangle$|⁠; while |$\langle B^2\rangle \sim \langle \overline{B}_{\rm y}^2\rangle \gg \langle b_{y}^2\rangle$| in tall boxes. As a result, α_mag roughly measures |b_x/b_y| in the standard box cases, but it traces a very different quantity, |$\left|b_{x}b_{y}\right|/\overline{B}_{\rm y}^2$| in the tall boxes. The effects of increasing the aspect ratio in an already elongated box would only introduce a stronger mean magnetic field. For instance, we find the total magnetic energy rises from 〈B²/8πP₀〉 ≃ 0.27 to 0.37 when comparing x1y4z4r128 and x1y4z8r128 in Table 2, and the energy increase mostly goes into |$\langle \overline{B}_{\rm y}^2/8\pi P_0\rangle$| so that the latter increases from ≃ 0.13 to 0.25 while the other magnetic field components stay constant.

In a recent review, Blackman & Nauman (2015) pointed out the stress-to-energy ratios for the pure mean and turbulent field could be different. Following their proposal, we define

\begin{equation} \alpha _{\rm mag}^{(m)}\equiv -2\frac{\langle \overline{B}_{x}\overline{B}_{y}\rangle }{\langle \overline{B}^2\rangle } \,; \alpha _{\rm mag}^{(t)}\equiv \frac{-2 \langle b_{x}b_{y}\rangle }{\langle B^2\rangle -\langle \overline{B}^2\rangle } \end{equation}

(11)

for the mean field and turbulent field, respectively. In the bottom panel of Fig. 8, we find the ratio of the turbulent field (larger symbols) still reaches a constant |$\alpha _{\rm mag}^{(t)}\sim 0.47$|⁠, while the mean field counterpart |$\alpha _{\rm mag}^{(m)}\ll \alpha _{\rm mag}^{(t)}$|⁠, and decreases linearly with L_z/L_x. These are consistent with the results found in Blackman & Nauman (2015).

Figure 8.

Top: α_mag ≡ −2〈〈B_xB_y〉〉_t/〈〈B²〉〉_t for various L_z/L_x and resolutions. Again, the vertical dotted line delineates the short and tall box regimes. The horizontal dashed line shows that α_mag stays constant in the small box regime; the dashed lines at L_z/L_x > 2.5 follows the (L_z/L_x)^−1/2 power law. Bottom: stress-to-energy ratios computed with pure mean field (small symbols) and turbulent field (large symbols) as defined in 11. The turbulent ratio keeps constant |$\alpha _{\rm mag}^{(t)}\sim 0.47$| (the horizontal dashed line); while its mean field counterpart |$\alpha _{\rm mag}^{(m)}\propto (L_{z}/L_{x})^{-1}$| (the dashed line at L_z/L_x > 2.5), much smaller than |$\alpha _{\rm mag}^{(t)}$|⁠.

Varying the box size with fixed aspect ratio

For given box size L_y = 4H, L_x = H in the horizontal domain, we find the stresses converge in the tall box runs (L_z/L_x ≳ 2.5) owing to the emergence of large-scale azimuthal magnetic field produced by dynamo action that acts as a local non-zero net flux. Since the mean field sets an extra length-scale that is directly related to the box size, it is also important to see how does the stress depend on the box size itself. Therefore we have performed an additional series of runs in which we vary L_x, but keep the aspect ratio L_z/L_x fixed.

In general, we find the amplitude of the saturated stress still con-verges to a constant value independent of numerical resolution as long as L_z/L_x > 2.5 (see Table 1). However, the value to which the stress converges depends on the box size, as found previously for shearing box simulations with net azimuthal field (Hawley et al. 1995; Guan et al. 2009). As shown in Fig. 9, both α_tot and 〈B²/8πP₀〉, measured from runs with aspect ratio of 1 : 4 : 4 scale as |$\propto L_{x}^2$|⁠. This differs from the saturation predictor used in Hawley et al. (1995), in which a linear relation ∝ L_y is reported. However, their relation applied to relatively weak (β ∼ O(10²)) external azimuthal field is very different than the very strong (β ∼ O(10)), dynamo generated oscillating azimuthal field observed in our runs.

Figure 9.

The box size dependence of the stress (left) and magnetic energy (right) for L_z/L_x = 4 runs. We note that for both panels, the x-axes are (L_x/H)² instead of L_x. The dashed lines shows the ∝ (L_x/H)² scaling. A strong linear correlation between the stress/energy and (L_x/H)² is obtained which indicates a way to renormalize α and 〈B²〉 using the squared length-scale |$L_{x}^2$|⁠.

In fact, the unstratified box does not have any intrinsic length-scales other than the box size L_x, cell width Δx, and (for compressible flows) the sonic scale c_s/Ω. If we assume that compressible effects can be ignored (we return to this point below), then based on a dimensional analysis the stress can be normalized by ρ₀(qΩL_x)² (Fromang & Papaloizou 2007; Guan et al. 2009), giving stress |$\propto L_{x}^2$| as reported above. As long as the simulation is resolved L_x/Δx ≳ 32, we find negligible dependence of the stress on numerical resolution, e.g. Fig. 9. We find similar scaling with size for other quantities such as magnetic and kinetic energy, and Reynolds stress. In addition, the scaling |$\propto L_{x}^2$| applies to runs using other values of the aspect ratio, e.g. L_z/L_x = 8 as tabulated in Table 1.

In Fig. 10, we plot the Fourier power spectra of five different runs computed with different numerical resolutions and physical box size (both L_x = H and L_x = 0.5H), but normalizing both the spectra and wavenumber according to the above scaling. Interestingly, all runs contain the same amplitude and slope at small k (the ‘inertial range’), but extend to successively larger k as the resolution is increased. Moreover, run x0.5y2z2r64 recovers the results of run x1y4z4r32 identically, and run x0.5y2z2r128 resembles x1y4z4r64 as well. Since these two pairs of runs have the same range of scales L_x/Δx, we expect them to show similar behaviour provided there are no intrinsic length-scales in the model. The result clearly indicates that indeed the properties of MRI turbulence in the unstratified shearing box model depend only on dimensionless wavenumber kL. As a result, similar properties found in previous sections for tall boxes of L_x = H are also present in tall boxes with different L_x.

Figure 10.

Rescaled magnetic energy power spectra for boxes with same aspect ratio L_z/L_x = 4. The close matches between x1y4z4 and x0.5y2z2r64, x1y4z4r64 and x0.5y2z2r128 indicate that the underling turbulence is intrinsically the same but only rescaled to a different box size.

Returning to Fig. 9, the stress measured in our simulations deviates slightly from the (L_x/H)² power law scaling at α ≳ 0.1, which occurs in the largest boxes. We speculate this is due to compressibility effects (Sano et al. 2004). Such large values of α are associated with very strong (nearly sonic) turbulence. For example, the rms density fluctuation for runs with large L_x become as large as 〈(δρ/ρ₀)²〉^1/2 ∼ 0.25 in x2y8z8r64. Compressibility strongly damps the turbulence, and prevents the turbulent magnetic field from growing even stronger. Obviously compressibility effects must play a role at some point, as turbulence alone cannot generate α ≫ 1.

Explicit dissipation

Previous work has shown that if explicit dissipation is included in the unstratified shearing box model with no net flux using the standard box size L_z/L_x = 1, then a converged value for the magnetic stress can be achieved (Fromang et al. 2007). The saturation level of the stress in this case is highly sensitive to the magnetic Prandtl number Pm. It increases almost linearly with Pm in the range 4 ≤ Pm ≤ 16; while turbulence is completely suppressed once the magnetic Prandtl number is below some critical value Pm_c ∼ 2–4 (Fromang et al. 2007; Simon & Hawley 2009). Further study suggests this behaviour shares similar origin to supertransient behaviour in chaotic systems, and that the lifetime of the turbulent active phase of the MRI increases exponentially with magnetic Reynolds number for fixed kinetic Reynolds number (Rempel, Lesur & Proctor 2010; Riols et al. 2013). In addition, studies with net flux and explicit dissipation have also shown a somewhat weaker dependence on the magnetic Prandtl number (Lesur & Longaretti 2007; Simon & Hawley 2009), with saturation levels appearing to reach asymptotic values in the Pm ≪ 1 limit (Meheut et al. 2015). These results have been used to argue it may be necessary to include explicit dissipation in all simulations of the non-linear regime of the MRI, even at very large magnetic Reynolds number and when the magnetic Prandtl number Pm ∼ 1.

However, we have shown above that for large aspect ratios (L_z/L_x ≳ 1), converged values of the stress are achieved in ideal MHD. To test whether the saturated stress is affected by the inclusion of explicit dissipation, we have repeated simulations in a tall box (L_x, L_y, L_z) = (1, 4, 4) H but with four different combinations of kinematic viscosity and Ohmic resistivity that give |$(\rm {Re},\rm {Pm}) = (1600,7.8125)$| for run x1y4z4r128pm8, (3125, 4) for x1y4z4r128pm4, (3125, 2) for x1y4z4r128pm2, and (3125, 1) for x1y4z4r128pm1, with resolution of 128/H. Parameters are chosen to match previous simulations in Fromang et al. (2007) and Simon & Hawley (2009). A resolution of 128/H ensures that the small-scale dissipation is well resolved in the parameter space we explored (Simon & Hawley 2009). The results of all four runs are listed in Table 1. We note all runs are simulated with the same initial conditions as described in Section 2.2 except x1y4z4r128pm8, which is restarted from t = 50 orbits of run x1y4z4r128pm2 to reduce the computational cost.

We find for Pm ≳ 4 (run x1y4z4r128pm4 and x1y4z4r128pm8), sustained stress and magnetic energy are achieved in tall boxes. Moreover, the saturated stress and magnetic energy for these runs are not much different from those without explicit dissipation. The time- and volume-averaged total stress is ≃0.085 for run x1y4z4r128pm4, and 0.092 for x1y4z4r128pm8, comparing to 0.084 for run x1y4z4r128 without explicit dissipation (see Fig. 11 and Table 1). We find, in Fig. 12, both the volume-averaged stresses and magnetic energy spectra are quite similar to the tall box runs without explicit dissipation as shown in right columns of Figs 2 and 3. It would seem that for taller boxes, in which the non-linear regime and saturation is controlled by dynamo action at least in the Pm ≳ 4 regime, including explicit dissipation has little effect on the results. We note that the saturated stress, α_tot ∼0.085 owning to the large-scale mean field, t is much greater than the values reported in run SZRe3125Pm4 of Simon & Hawley (2009) (∼0.013), and in run 128Re3125Pm4 of Fromang et al. (2007) (∼0.009). We also note that the shape of the power spectrum is different than the standard short box simulations with explicit dissipation (Fromang 2010) as the large-scale dynamo observed in our tall box runs is completely absent in this latter case.

Figure 11.

Left: time evolution of volume-averaged total stress for different magnetic Prandtl number Pm. Right: time-averaged total stress-to-pressure ratio as a function of Pm. All curves are calculated using resolution 128/H. We set α_tot = 0 for Pm = 1 case (run x1y4z4r128pm1), as turbulent transport vanishes after 250 orbits. The saturated stress level is sensitive to the magnetic Prandtl number when Pm ≲ 4; however, this dependence becomes less significant when Pm ≳ 4, unlike the linear dependence found in standard boxes (Fromang et al. 2007; Simon & Hawley 2009).

Figure 12.

Similar to Figs 2 and 3 but with explicit dissipation: Re = 3125 and Rm = 12500. The results are very similar to the runs without explicit dissipation in Figs 2 and 3.

When Pm = 2 (run x1y4z4r128pm2), the stress in Fig. 11 appears to be more bursty compared to higher Prandtl number runs and those without explicit dissipation. The large fluctuations might be a result of competition between supertransient decay and dynamo growth (Thaler & Spruit 2015). Time averaging over the last 200 orbits, we find α_tot ∼ 0.05, almost twice smaller than high Pm runs. The same parameters have also been explored previously using a standard box, see run 128Re3125Pm2 in fig. 8 of Fromang et al. (2007) and SZRe3125Pm2 in table 1 of Simon & Hawley (2009), in contrast to our tall box run x1y4z4r128pm2, no sustained MHD turbulence is observed in either of these cases.

The only run in which we find turbulence eventually decays is Pm = 1 (run x1y4z4r128pm1 in Fig. 11). After about 250 orbits, both the Maxwell stress and magnetic energy abruptly drop by several orders of magnitude and the flow becomes laminar. This finding confirms the existence of a critical Prandtl number Pm_c for a given kinetic Reynolds number even when the box is tall (L_z/L_x > 1 and L_y/L_x = 4). Together with the results of Pm = 2 run, it indicates a smaller critical value, 1 ≲ Pm_c ≲ 2, comparing to Pm_c = 2–4 found in standard boxes (Fromang et al. 2007; Simon & Hawley 2009). As we mainly varies the Re_M, we caution the readers that the results might possibly reflect the dependence of saturation amplitude of the MRI on another dimensionless number, the Lundquist number |$\rm {Lu}\equiv v_{\rm Az}^2/\eta \Omega$| (Sano & Stone 2002). We find the time-averaged Lu ≃ 332, 105 and 11 for run x1y4z4r128pm4, x1y4z4r128pm2 and x1y4z4r128pm1,¹ getting closer to the critical values reported in previous studies with vertical net flux (Lesur & Longaretti 2007; Pessah, Chan & Psaltis 2007; Turner, Sano & Dziourkevitch 2007; Masada & Sano 2008).

The result of Rempel et al. (2010, see their fig. 5), based on the statistics of MRI turbulence in a standard (L_x: L_y: L_z) = (1: π: 1) box, would predict a characteristic decay time ∼10 shear time units for Pm = 1 and 2 with the same Re = 3125 as used here. However, we find sustained turbulence that lasts more than several thousand shear time units in our Pm = 2 run; and the active time for our Pm = 1 run is ∼200 orbits, or 2000 shear time units, much longer than the prediction. Moreover, Fromang et al. (2007) find, when changing Re from 12 500 to 25 000, the turbulence lifetime of Pm = 1 increases from ≲40 orbits to ∼100 orbits. This would suggest a decaying time even shorter than 40 orbits for Re = 3125 case, therefore much smaller than our result. Clearly, the aspect ratio of the computational domain may play a big role in determining the dynamical lifetime of MRI turbulence (Riols et al. 2015) in the unstratified, no net flux shearing box.

To sum up, we find the Prandtl number dependence of the turbulent stress in our tall box runs is very different from that found earlier using a standard box (see right-hand panel of Fig. 11). When Pm ≥ 4, the saturated stress level is insensitive to the inclusion of explicit dissipation. We find converged turbulent stress for runs with and without explicit dissipation. The saturated stress shows stronger dependence when Pm < 4 as it drops linearly from Pm = 4 to 2. By choosing a taller box which promotes dynamo action, we find the dynamical lifetime of MHD turbulence is extended, and the critical magnetic Prandtl number is reduced compared to the standard smaller box simulations.

MODELING THE DYNAMO

Our previous results indicate that shearing box simulations of the MRI with a large aspect ratio, dynamo action produces strong ordered toroidal fields on vertical scales L_x. In this section we explore models that might explain this dynamo action.

Dynamo cycle period

An important property of the dynamo is that at any given vertical location z, the toroidal magnetic field is cyclic, and the period of these cycles is an important clue to the mechanism of dynamo action. As illustrated in Fig. 4, cyclic patterns of the mean field |$\overline{B}_y$| becomes strong and regular for those boxes with L_z ≳ 8 H. It is relatively easy to identify the cycle period from the space–time diagrams. However, for smaller boxes, e.g. L_z = 4 H as shown in Fig. 4, the cycle is less regular and it is more difficult to extract a single value of the period T_cycle. Thus, we measure the cycle period based on the power spectrum density (PSD) of the largest vertical mode (k_z = 2π/L_z) of |$\overline{B}_y (t)$| in Fourier space, i.e.

\begin{equation} {\rm PSD\, of }\, \widetilde{B}_y(t),\, {\rm where }\, \widetilde{B}_y (t) \equiv \int \!\overline{B}_y(z;t) e^{-i k_{\rm z} z} \mathrm{d}z \,. \end{equation}

(12)

We can also estimate errors from the range of T_cycle as measured from two independent time sequences. We have applied this measurement to all simulations regardless of the box size, but note the physical meaning of T_cycle in simulations with a small aspect ratio (L_z/L_x < 2) is less clear due to the absence of dynamo action. We list all values of T_cycle measured in this way in Table 1 for reference.

Once dynamo is at action, there is in general a phase lag between the radial mean field |$\overline{B}_x$| and the azimuthal |$\overline{B}_y$|⁠. In order to obtain the phase shift, we first apply filter on our mean field and get the largest vertical mode |$\overline{\boldsymbol {B}}^{(1)}(z,t)$|⁠, where ‘(1)’ denotes k_z = 2π/L_z only. We then calculate the cross-correlation between two time sequences, |$\overline{B}_x^{(1)}$| and |$\overline{B}_y^{(1)}$| at given z. The exact time delay Δt between radial and azimuthal mean field are the lag which gives the most negative value of correlation. In Table 1, we list the vertically averaged Δt, and the errors from the variation of Δt at different height. The azimuthal mean field therefore lags behind the radial field by

\begin{equation} \phi = \pi - 2\pi \frac{\Delta t}{T_{\rm {cycle}}} \end{equation}

(13)

in phase.

In the top panel of Fig. 13, we plot the measured T_cycle for all L_x = 1 H runs. In general, we find the cycle period is much longer than the dynamical (orbital) time, ranging from several to tens of orbital periods. The measured T_cycle becomes longer as L_z/L_x is increases with a slope slightly shallower than unity. In contrast to the aspect ratio dependence, T_cycle is not very sensitive to the box size. As shown in the middle panel of Fig. 13, T_cycle for L_z/L_x = 4 scatters between ∼20 and 30 orbits, while it varies around ∼50 orbits for L_z/L_x = 8 with ≲50 per cent fluctuations. The phase lag between radial and azimuthal mean field does not depend on box size either (see Table 1). It only weakly depends on the aspect ratio, roughly ∝ (L_z/L_x)^−0.2 as shown in the bottom panel of Fig. 13.

Figure 13.

Top: dynamo cycle periods for different L_z/L_x with fixed L_x/H = 1, dashed lines show ∝ L_z/L_x power law. Middle: dynamo periods versus box sizes (small symbols for L_z/L_x = 4, and larger symbols for L_z/L_x = 8). The dotted horizontal lines at T_cycle = 27 and 48 orbits show the constancy of cycle period with respect to various box sizes. Bottom: phase lag between radial and azimuthal mean field as calculated in equation (13), dashed lines follow ∝ (L_z/L_x)^−0.2 scaling. The dotted vertical lines at L_z/L_x = 2.5 separate the standard and tall boxes. The error bars of T_cycle are the variations measured with two independent time sequences; the error bars of phase lag show the max/min values given the variations of the measured T_cycle and time delay Δt.

A toy model

Following the model proposed in Lesur & Ogilvie (2008b), we also try to fit the dynamo cycles observed in our simulations with the following non-linear model:

\begin{equation} \mathrm{\partial} _t\overline{B}^{(1)}_x= \gamma \overline{B}^{(1)}_y(t-t_r)\frac{\left| \overline{B}^{(1)}_y(t-t_r)\right|-B_r}{B_r}\,, \end{equation}

(14)

\begin{equation} \mathrm{\partial} _t\overline{B}^{(1)}_y= -q\Omega \overline{B}^{(1)}_x- \beta \overline{B}^{(1)}_y(t-t_r)\,, \end{equation}

(15)

where |$\overline{B}^{(1)}_x$| and |$\overline{B}^{(1)}_y$| are horizontally averaged radial and azimuthal field that are filtered to conserve only the largest vertical mode (k_z = 2π/L_z, represented with the superscript ‘(1)’). Here t_r characterizes the time delay between the EMF and large-scale magnetic field. When the mean field |$|\overline{B}^{(1)}_y|$| is smaller than B_r, the γ term effectively amplifies the |$\overline{B}^{(1)}_x$| field, which in turn lead to |$\overline{B}^{(1)}_y$| growth via shear. In the opposite case, when |$|\overline{B}^{(1)}_y|$| exceeds B_r, the γ term starts to damp the radial field, and the azimuthal field is thus reduced via the β term. In the fitting, we choose t_r = 2/|qΩ| = 4Ω⁻¹/3 as suggested by linear analysis (Lesur & Ogilvie 2008a,b), and |$B_r/(\sqrt{4\pi \rho _0}c_{\rm {s}})=0.3$|⁠, a factor ∼2.5 greater than used in Lesur & Ogilvie (2008b) to better fit the mean field amplitude observed in our simulations. We then vary γ and β to match the observed cycle period and long-term amplitude. Empirically, γ determines the cycle frequency, and β prevent the solution from diverging. We focus on fitting the results for run x1y4z8r32 as it exhibits a well organized cyclic pattern. Without loss of generality, we fit the filtered |$\overline{B}^{(1)}_x$| and |$\overline{B}^{(1)}_y$| at z = −3 H. With γ = 2 × 10⁻⁴ and β = 4 × 10⁻⁴, the model is shown as the red curves in Fig. 14. For times after 100 orbits, the model provides a reasonable fit to the spatial filtered simulation data. The model predicts the correct dynamo period and relative strength between |$\overline{B}^{(1)}_x$| and |$\overline{B}^{(1)}_y$|⁠. However the phase shift between |$\overline{B}^{(1)}_x$| and |$\overline{B}^{(1)}_y$| is π/2, slightly smaller than the simulation data, which is ϕ ≃ (0.65 ± 0.08)π. This phase lag is also different than what is found in Herault et al. (2011) ( ∼ π), but is close to 3π/4, the phase shift of a marginally excited dynamo wave in a α-Ω dynamo (Brandenburg & Subramanian 2005). As the values of B_r and t_r are uncertain, there are likely to be degeneracies in the model parameters. Recent work has explored this model as a magnetic analogue to the ‘shear current’ effect proposed by Rogachevskii & Kleeorin (2004) in which the non-isotropic contribution of the turbulent diffusion can drive large-scale magnetic field growth in non-helical turbulent shearing flows. Interestingly, Squire & Bhattacharjee (2015b) recently find negative off-diagonal turbulent diffusivity (a negative η_yx term which drives exponential growth of B_x, see their equations (2) and (3), or equation (16) in this paper for its definition) by direct measurement in simulations assuming a simple closure model,

\begin{equation} {\mathcal {E}}_i = \alpha _{ij}\overline{B}_j - \eta _{ij}\epsilon _{jkl}\mathrm{\partial} _k\overline{B}_l\,, i,j,l \in \lbrace x,y\rbrace \,, k=z \end{equation}

(16)

where the turbulent EMF, |$\mathcal {E}\equiv \overline{v\times b}$|⁠, and mean field |$\overline{\boldsymbol {B}}$| can be measured from the simulation data, and ϵ_jkl is the Levi–Civita permutation symbol. After setting α_xx = α_yy, α_yx = 0, η_xx = η_yy and η_xy = 0, they find an average η_yx ∼ −10⁻⁴H²Ω. Following their procedure, i.e. computing |$\langle \mathcal {E}_i M\rangle$| from equation (16) for each of |$M = (\overline{B}_x, \overline{B}_y, \mathrm{\partial} _z\overline{B}_x, \mathrm{\partial} _z\overline{B}_y)$| and solving the resulting matrix equations at each time point, we also find negative η_yx ≃ −1.6 × 10⁻³H²Ω by time averaging over the dynamo's growth phase t = 10–50 orbits in x1y4z8r32.

Figure 14.

A toy model fit (red, based on equations 14 and 15) to our unstratified x1y4z8 run (black). The top panel shows the largest vertical mode of mean |$\overline{B}_y$| at z = −3 H; the bottom panel shows the corresponding radial field. Dotted lines in both panels mark zero amplitude level.

Alternatively, we can also solve equation (16) at any given z directly assuming time-independent transport coefficients (but a function of height) with multiple linear regression method. The results are shown in Fig. 15. All αs fluctuate about zero and are close to zero when averaging over z. The diagonal η_xx ∼ 5 × 10⁻³H²Ω are dissipations of magnetic field. The non-diagonal term η_yx ∼ −10⁻³H²Ω for all z is consistent with that of Squire & Bhattacharjee (2015b), indicates a dynamo growth. We find this negative η_yx behaviour does not change with non-zero αs and ηs, or unequal α_xx and α_yy, as long as the diagonal ηs are kept to be the same (η_xx = η_yy). We also obtain similar values of η_yx in the well established dynamo cyclic phase using data at t = 100–200 orbits.

Figure 15.

Turbulent transport coefficients measured via simulation at t = 10–50 orbits, assuming closure model like equation (16). We further set α_xx = α_yy, α_yx = 0, η_xx = η_yy and η_xy = 0 in this fitting. Square symbols are the time-averaged coefficients at evenly spaced disc height. The black dashed lines show the vertical averages and the grey dashed lines are set to zero to help read this plot. Clearly, negative η_yx is obtained for all z indicating an antidiffusive growth of magnetic field.

However, comparing a time dependent η_yx throughout a dynamo cycle² shows no evidence for sign change when the mean field |$|\overline{B}_y|$| surpasses some critical magnetic field strength as designed in the toy model, although we do see the filtered EMF |$\mathcal {E}_y^{(1)}$| changes sign over a cycle with the projection method used in Lesur & Ogilvie (2008b, see their eqautions 12 and 13) and after some smoothing. Further study is still required to explain the inconsistency we find here, and identify the dynamo mechanisms ultimately, but is beyond the scope of this paper.

Connection to the stratified shearing box

It is interesting to compare the dynamo cycles in our unstratified shearing box simulations with those observed in stratified discs. Simulations using the stratified shearing box always find strong dynamo cycles in the toroidal field (e.g. Brandenburg et al. 1995; Stone et al. 1996; Ziegler & Rüdiger 2000; Gressel 2010; Oishi & Mac Low 2011; Käpylä et al. 2013). As the unstratified shearing box discs approximates the midplane of a stratified discs, there may be some relation between the dynamo behaviour seen in both cases.

We start by analysing an existing stratified simulation data (run STD32) first published in Shi et al. (2010). Due to buoyancy, the mean magnetic field is expelled towards the disc surface producing a regular pattern known as the butterfly diagram (see fig. 6 of Shi et al. 2010) reproduced here in Fig. 16. By extracting data within ±H of the midplane, we find a clear dynamo cycle with a period ∼10 orbits as shown in the bottom panel of Fig. 16. At any given time, the mean azimuthal field is roughly uniform in z, and has one sign throughout the subvolume near the midplane (to preserve the constraint of no net flux, field of the opposite sign is located above and below the midplane). We can also model the dynamo in the stratified shearing box with the toy model described in Section 4.2. Again, we adopt t_r = 4/(3Ω) and |$B_r/(\sqrt{4\pi \rho _0}c_{\rm {s}}) = 0.3$|⁠, where ρ₀ and c_s are the midplane density and sound speed in the stratified disc. We fit the vertically averaged |$\overline{B}_x$| and |$\overline{B}_y$| with parameter values of γ = 10⁻³ and β = 0.014 as shown in Fig. 17. Similar to the unstratified case, the toy model captures the cycle period and relative strength of the field components correctly, but slightly underestimates the relative phase between |$\overline{B}_x$| and |$\overline{B}_y$| (⁠|$\phi /\pi \simeq 0.6^{+0.1}_{-0.2}$|⁠).

Figure 16.

Top: space–time diagram of |$\overline{B}_y$| in a stratified box STD32 (Shi et al. 2010). Bottom: a blow-up view of the left to show the main body of the disc within ± H.

Figure 17.

A toy model fit (based on equations 11 and 12) to the stratified disc STD32 (see Fig. 16; Shi et al. 2010). The fit is in red, the vertically averaged mean B_y (top) and B_x (bottom) within ± H are in black. Dotted lines in both panels mark zero amplitude level.

The overall similarity between the dynamo in the unstratified and stratified discs suggests the same mechanism could act in both cases. We speculate that the core of the stratified disc resembles our tall box runs. Extrapolating the linear scaling found in the top panel of Fig. 13 to L_z/L_x ≲ 2 gives T_cycle ∼ 10 orbits which matches the cycle period of the stratified shearing box very well. Assuming similar correlation between EMF and mean field as described in Section 4.2 (Squire & Bhattacharjee 2015b,a), we find again a negative off-diagonal resistivity η_yx ∼ −10⁻³H²Ω (Fig. 18), which favours an antidiffusive dynamo model similar to the unstratified case. However, we find, unlike the unstratified case, an α-effect is also present in the stratified box. The diagonal α_yy is antisymmetric with respect to the midplane. In the upper half of the disc, it is negative below 4 H, and systematically positive above 4 H. The negative α-effect (within the upper half of the disc) is generally attributed to the dynamo growth in previous stratified simulations and the negative sign is caused by the influence of the field buoyancy and/or strong shear in accretion discs (e.g. Brandenburg et al. 1995; Rüdiger & Pipin 2000). The positive α-effect has also been recently proposed to explain the dynamo (Gressel 2010; Gressel & Pessah 2015). Together with the η-effect, it seems to suggest a combination of direct and indirect dynamo mechanisms (Blackman & Tan 2004; Gressel 2010) in stratified shearing box. Further investigations are required to identify whether one or both mechanisms could exist and further cause the dynamo cycles.

Figure 18.

Similar to Fig. 15, turbulent transport coefficients measured via stratified disc as in Fig. 16 are plotted assuming only η_xx = η_yy. In each panel, different colours represent different coefficients as labelled on its y-axis. We see similar negative η_yx as in the unstratified disc within ±4 H, also a negative (positive) α_yy below (above) 4 H in the upper half of the disc, antisymmetric with respect to the midplane.

CONCLUSIONS

In this paper, we have studied MRI turbulence in the unstratified shearing box model with no net flux, and in particular how the size and aspect ratio of the computational domain affects the saturation amplitude and stress. We have performed simulations both with and without explicit dissipation. Previous simulations of this particular shearing box model have shown that the saturation amplitude of the MRI depends on numerical resolution. However, we have shown that the geometry of the computational domain also strongly affects the result. Our main results are summarized as follows.

• The amplitude of the stress and turbulence driven by the MRI converge to values independent of numerical resolution if the aspect ratio of the computational domain L_z/L_x > 2.5, even without explicit dissipation. These values are proportional to the box size |$L_{x}^2$|⁠. Only when the aspect ratio L_z/L_x < 2 do these levels depend on resolution.

• When L_x = H and L_z/L_x > 2.5 the converged value of the saturated stress is greatly increased compared to previous values reported for L_z/L_x = 1, and can be as large as α ≳ 0.1. These values are similar to those required by observations of dwarf novae discs in outburst, perhaps negating some of the concerns expressed in King et al. (2007).

• The Fourier power spectra of runs with large aspect ratio are very different from those with L_z/L_x ∼ 1. In particular, most of the power in the magnetic field is at small k, independent of resolution.

• For the limited range of parameter values explored in this paper, we find the saturation level of the MRI is independent of whether explicit dissipation (viscosity and resistivity) is included or not, provided that the magnetic Reynolds number is large. Turbulence can be suppressed at small Re_M, corresponding to small magnetic Prandtl number Pm ≲ 4, however simulations utilizing a taller box show a relatively smaller critical Prandtl number and much longer lifetime of turbulence than those with L_z/L_x ∼ 1.

• Cycles of strong large-scale toroidal magnetic field with alternating sign on scales of a few L_x are generated in domains with L_z/L_x > 2.5 via a dynamo. This large-scale magnetic field sets and sustains the turbulence independent of numerical resolution.

• Some aspects of the cyclic dynamo can be modelled with anisotropic turbulent resistivity (Lesur & Ogilvie 2008b), although the physical mechanisms driving the dynamo are unclear. Similar dynamo cycles observed in a stratified shearing box simulation suggest the stratified and unstratified discs might share common dynamo mechanisms.

As discussed in the Introduction, the unstratified zero-net-flux shearing box model is unlikely to have much relevance to real astrophysical discs. Instead, it is likely that different local (L ∼ H) patches of the disc are threaded by a broad range of net field strengths, even if the total field threading the disc is zero. Thus, global models of astrophysical discs are best constructed from an ensemble of net flux shearing box simulations. Nevertheless, studies of the unstratified zero-net-flux shearing box are of interest in order to study MHD dynamo action in a simple, well-posed model. Moreover, we have shown that some aspects of the dynamo observed in stratified shearing boxes may be related to properties of the dynamo discussed here. Thus, further study of the MRI dynamo in domains with large aspect ratio, in both compressible and incompressible MHD, and both with and without explicit dissipation are warranted.

We thank Steve Balbus, Eric Blackman, Charles Gammie, Geoffroy Lesur and Jonathan Squire for thoughtful comments on an early draft of this paper. We also thank Amitava Bhattacharjee, Fatima Ebrahimi, Sebastien Fromang and Geoffroy Lesur for encouraging discussions. We thank Shane Davis for sharing the code to compute the power spectrum. This work was supported in part by the National Science Foundation under grant PHY-1144374, ‘A Max-Planck/Princeton Research Center for Plasma Physics’ and grant PHY-0821899, ‘Center for Magnetic Self-Organization’. Resources supporting this work were provided by the Princeton Institute of Computational Science (PICSciE) and Engineering and Stampede at Texas Advanced Computing Center (TACC), The University of Texas at Austin through XSEDE grant TG-AST130002.

1

Take time average over turbulent active phase t < 250 orbits.

2

It is calculated again by fitting moment equations of the closure 16 at each time point (without time average) using |$\mathcal {E}$| and |$\boldsymbol {B}$| values directly from the simulation. We can also fit equation (16) with filtered data, e.g. by keeping only the first few wavenumbers. The resulting η_yx can switch sign, however is so noisy that we can hardly retrieve any meaningful results from it.