RESISTIVITY-DRIVEN STATE CHANGES IN VERTICALLY STRATIFIED ACCRETION DISKS

In all of our simulations, ν and η are parameterized in terms of Reynolds numbers. At a scale height, H, away from the radial center of the shearing box, the fluid velocity is |v_y| ∼ qHΩ ∼ c_s. Thus, the sound speed is a representative velocity for the fluid, and we define the Reynolds numbers as

$$\begin{equation} Re \equiv \frac{c_{\rm s}H }{\nu }. \end{equation} \tag{ 5 }$$

Similarly we define the magnetic Reynolds number,

$$\begin{equation} Rm \equiv \frac{c_{\rm s}H}{\eta }, \end{equation} \tag{ 6 }$$

and their ratio, the magnetic Prandtl number,

$$\begin{equation} P_{\rm m}\equiv \frac{\nu }{\eta } = \frac{Rm}{Re}. \end{equation} \tag{ 7 }$$

Because resistivity affects the MRI directly, another useful dimensionless quantity is the Elsasser number

$$\begin{equation} \Lambda \equiv \frac{v_{\rm A}^2}{\eta \Omega }. \end{equation} \tag{ 8 }$$

Λ can be computed on a zone-by-zone basis, but it is also helpful to calculate an average Λ for a simulation. Thus, we can write the characteristic ${\rm Alfv\acute{e}n}$ speed in the direction i via the averaged magnetic field in that direction:

$$\begin{equation} v_{{\rm A}i} = \sqrt{\frac{\big\langle B_i^2\big\rangle }{\langle \rho \rangle }}, \end{equation} \tag{ 9 }$$

where the angled brackets denote a volume average and the subscript i = (x, y, z) depending on the direction of interest. In the above definition of Λ, all three components are used in calculating the ${\rm Alfv\acute{e}n}$ speed. As the Elsasser number approaches unity, resistivity becomes more dominant, stabilizing the MRI. Simon & Hawley (2009) found that turbulence decayed in net toroidal field, unstratified shearing boxes when the Elsasser number was ≲100.

Although we focus on physical dissipation, there is numerical dissipation as well. One way to characterize the effects of finite resolution is to compare the characteristic MRI wavelength, λ_MRI = 2πv_A/Ω, to the grid zone size. Noble et al. (2010) defined this "quality factor" Q as

$$\begin{equation} Q_i \equiv \frac{\lambda _{{\rm MRI},i}}{\Delta x_i} = \frac{2 \pi v_{{\rm A}i}}{\Omega \Delta x_i}. \end{equation} \tag{ 10 }$$

Again, v_Ai is defined via Equation (9). For Q_i ≲ 6, the growth of the MRI can be reduced (Sano et al. 2004) and the MRI is underresolved, although this number has some uncertainty and should be taken only as an estimate.

One can also calculate a volume-averaged ${\rm Alfv\acute{e}n}$ speed rather than an ${\rm Alfv\acute{e}n}$ speed determined by the averaged field, i.e.,

$$\begin{equation} v_{{\rm A}i} = \left\langle \frac{|B_i|}{\sqrt{\rho }}\right\rangle. \end{equation} \tag{ 11 }$$

For convenience, most of our analysis uses Equation (9), as the volume-averaged history data and the one-dimensional, horizontally averaged quantities are routinely computed at high time resolution for use in creating, e.g., space–time plots. As a check, we have calculated the volume-averaged ${\rm Alfv\acute{e}n}$ speed via Equation (11) for several hundreds of orbits in a few of our vertically stratified simulations (the volume average is done within 2H of the mid-plane). We found that v_A from Equation (11) is roughly 1–1.5 times larger than that from Equation (9).

We have run shearing box simulations with different field geometries, dissipation values, and resolutions, both with and without vertical gravity. Here, we describe the parameters and initial conditions used in these two types of simulations. The results from the unstratified simulations are presented in Section 3 and the results from the stratified simulations are in Section 4.

2.2.1. Sine Z Unstratified Simulations

The first set of simulations is unstratified, zero-net magnetic flux shearing boxes. This is the same type of problem as investigated in Fromang et al. (2007) and Simon & Hawley (2009). The initial magnetic field is

$$\begin{equation} {\mbox{\boldmath {$B$}}} = \sqrt{2 P_o/\beta } {\rm sin}(2 \pi x/L_x) \hat{\mbox{\boldmath {$z$}}} \end{equation} \tag{ 12 }$$

with β = 400. The isothermal sound speed is c_s = 0.001, corresponding to an initial gas pressure P_o = 10⁻⁶ with initial density ρ_o = 1. The orbital velocity of the local domain is Ω = 0.001. For these simulations, we define the scale height to be

$$\begin{equation} H \equiv \frac{c_{\rm s}}{\Omega } = 1, \end{equation} \tag{ 13 }$$

which is a slightly different definition than that in the stratified box (by a factor of $\sqrt{2}$). The size of the box is L_x = 1H, L_y = 4H, and L_z = 1H. We run several resolutions in order to study convergence, and at each resolution, we study four different cases of Rm and P_m values. See Table 1 for the parameters of these runs. The labeling scheme of the runs refers to resolution, field geometry, and dissipation values, e.g., 16SZRm12800Pm16 corresponds to 16 zones per H, "SZ" for the "Sine Z" geometry, and Rm = 12800, P_m = 16. Since, as we will see, Rm is the critical parameter in many of our simulations, we choose to include it as part of the run label, differing from the convention of related works (e.g., Fromang et al. 2007; Simon & Hawley 2009).

Table 1. Unstratified Simulations

Label	Re	Rm	Pm	Resolution	α	Description
				(zones per H)
16SZRm12800Pm16	800	12800	16	16	0.011	Zero-net flux
32SZRm12800Pm16	800	12800	16	32	0.033	Zero-net flux
64SZRm12800Pm16	800	12800	16	64	0.042	Zero-net flux
128SZRm12800Pm16a	800	12800	16	128	0.046	Zero-net flux
16SZRm12500Pm4	3125	12500	4	16	0.0043	Zero-net flux
32SZRm12500Pm4	3125	12500	4	32	0.013	Zero-net flux
64SZRm12500Pm4	3125	12500	4	64	0.015	Zero-net flux
128SZRm12500Pm4a	3125	12500	4	128	0.013	Zero-net flux
16SZRm6250Pm1	6250	6250	1	16	decay	Zero-net flux
32SZRm6250Pm1	6250	6250	1	32	decay	Zero-net flux
64SZRm6250Pm1	6250	6250	1	64	decay	Zero-net flux
16SZRm25600Pm2	12800	25600	2	16	0.0077	Zero-net flux
32SZRm25600Pm2	12800	25600	2	32	0.010	Zero-net flux
64SZRm25600Pm2	12800	25600	2	64	0.0078	Zero-net flux
32FTNum	...	...	...	32	0.021b	Flux tube, num. dissipation, βy = 100
32FTNumβ1000	...	...	...	32	0.020c	Flux tube, num. dissipation, βy = 1000
32FTNumβ10000	...	...	...	32	0.018c	Flux tube, num. dissipation, βy = 10000
32FTRm800Pm0.5	1600	800	0.5	32	decay	Restarted from 32FTNum
32FTRm3200Pm0.5	6400	3200	0.5	32	0.0094	Restarted from 32FTNum
32FTRm3200Pm2	1600	3200	2	32	0.018	Restarted from 32FTNum
32FTRm3200Pm4	800	3200	4	32	0.028	Restarted from 32FTNum
32FTRm6400Pm4	1600	6400	4	32	0.029	Restarted from 32FTNum
32FTRm6400Pm8	800	6400	8	32	0.041	Restarted from 32FTNum
32FTRm1600Pm1β1000	1600	1600	1	32	decay	Restarted from 32FTNumβ1000
32FTRm3200Pm2β1000	1600	3200	2	32	0.015d	Restarted from 32FTNumβ1000
32FTRm3200Pm2β10000	1600	3200	2	32	decay	Restarted from 32FTNumβ10000

Notes. ^aThese runs were taken from Simon & Hawley (2009). ^bTime averaged from orbit 20 to 150. ^cTime averaged from orbit 20 to 110. ^dTime averaged from orbit 120 to 400.

Download table as:  ASCIITypeset image

The MRI is seeded with random perturbations to the density and the velocity components introduced at the grid scale. The amplitude of the density perturbations is δρ = 0.01 and the amplitude of the velocity perturbations is (1/5)δρc_s for each component (a different perturbation is applied for each component). We do not employ orbital advection in these simulations (see the Appendix). All simulations are run to 400 orbits, except for the runs in which the turbulence decays and also 32SZRm12800Pm16 and 32SZRe12500Pm4, which were run to 289 orbits and 246 orbits, respectively. We also include the results from previous, higher resolution simulations (Simon & Hawley 2009) in our analysis.

2.2.2. Flux Tube Unstratified Simulations

The second set of unstratified simulations contain a net toroidal field and are initialized with the twisted azimuthal flux tube of Hirose et al. (2006), with minor modifications to the dimensions and β values. The initial conditions are designed to match those used in the vertically stratified simulations of Section 2.2.3, so we define H to be that of a stratified, isothermal disk,

$$\begin{equation} H = \frac{\sqrt{2} c_{\rm s}}{\Omega }. \end{equation} \tag{ 14 }$$

The isothermal sound speed c_s = 7.07 × 10⁻⁴, corresponding to an initial value of the gas pressure of P_o = 5 × 10⁻⁷. With Ω = 0.001, the value for the scale height is H = 1. As in the SZ simulations, random perturbations are added to the density and velocity components.

The initial toroidal field, B_y, is given by

$$\begin{equation} B_y = \sqrt{\frac{2 P_o}{\beta _y} - \big(B_x^2+B_z^2\big)} \end{equation} \tag{ 15 }$$

and we have run simulations with a toroidal field β of β_y= 100, 1000, and 10000. The poloidal field components, B_x and B_z, are calculated from the y component of the vector potential:

$$\begin{equation} A_y = \left\lbrace \begin{array}{@{}ll}- \sqrt{\frac{2 P_o}{\beta _p}}\frac{H}{2 \pi } \left[1 + {\rm cos}\left(\frac{2 \pi r}{H}\right)\right] & \quad \mbox{if \ $r < \frac{H}{2}$} \\ 0 & \quad \mbox{if \ $r \ge \frac{H}{2}$}, \end{array} \right. \end{equation} \tag{ 16 }$$

where $r = \sqrt{x^2 + z^2}$ and β_p = 1600 is the poloidal field β value.

The labeling convention is the same as in Section 2.2.1, but with "FT" for "Flux Tube" (see Table 1). We run simulations both with and without physical dissipation. Those simulations with no physical dissipation are labeled with "Num" for "Numerical dissipation," and for the β_y = 1000 (10000) cases, we append β1000 (β10000) to the end of the run label. The domain size is L_x = 2H, L_y = 4H, and L_z = 1H, and the resolution is 32 zones per H. Orbital advection is employed in these calculations (see the Appendix). 32FTNum was run to 150 orbits, and both 32FTNumβ1000 and 32FTNumβ10000 were run to 110 orbits.

The runs with physical dissipation are initiated from the turbulent state (at t = 100 orbits) of the corresponding β_y value run with only numerical dissipation. These simulations were all run out to 220 orbits.

2.2.3. Vertically Stratified Simulations

In vertically stratified shearing boxes, the initial density corresponds to an isothermal hydrostatic equilibrium

$$\begin{equation} \rho (x,y,z) = \rho _o {\rm exp}\left(-\frac{z^2}{H^2}\right), \end{equation} \tag{ 17 }$$

where ρ_o = 1 is the mid-plane density and H is the scale height in the disk, as defined in Section 2.2.2. A density floor of 10⁻⁴ is applied to the physical domain as too small a density leads to a large ${\rm Alfv\acute{e}n}$ speed and a very small timestep. Furthermore, numerical errors in energy make it difficult to evolve regions of very small plasma β. All other parameters and initial conditions are identical to the corresponding FT runs of Section 2.2.2, except that the toroidal field is set to zero where the poloidal field strength is zero. All the stratified runs use β_y = 100 for the initial toroidal field strength.

The domain size for all vertically stratified simulations is L_x = 2H, L_y = 4H, and L_z = 8H, and we have run most simulations at 32 grid zones per H, with a few at 64 zones per H. Orbital advection is employed in all of these runs (see the Appendix). The runs are listed in Table 2. The labeling scheme is the same as for the unstratified simulations, but we omit the "FT" label as all the stratified runs use the flux tube initial conditions. We study the effect of physical dissipation by restarting the equivalent resolution, numerical-dissipation-only run at t = 100 orbits.

Table 2. Vertically Stratified Simulations

Label	Re	Rm	Pm	Resolution	Integration time	Description
				(zones per H)	tstop–tstart (orbits)
32 Num	...	...	...	32	1058	Num. dissipation
32Rm800Pm0.5	1600	800	0.5	32	263	...
32Rm3200Pm0.5	6400	3200	0.5	32	863	...
32Rm3200Pm2	1600	3200	2	32	1082	...
32Rm3200Pm4	800	3200	4	32	487	...
32Rm6250Pm1	6250	6250	1	32	337	...
32Rm6400Pm8	800	6400	8	32	325	...
32Rm6400Pm4	1600	6400	4	32	488	...
32Rm3200Pm2_By+	1600	3200	2	32	584	By added at 150 orbits
32ShearBx	...	...	...	32	45	Net Bx within mid-plane
64Num	...	...	...	64	159	Num. dissipation
64Rm800Pm0.5	1600	800	0.5	64	108	...
64Rm3200Pm0.5	6400	3200	0.5	64	83	...
64Rm3200Pm2	1600	3200	2	64	80	...
64Rm3200Pm4	800	3200	4	64	84	...
64Rm6400Pm4	1600	6400	4	64	80	...

Download table as:  ASCIITypeset image

We begin with a resolution study of unstratified simulations with physical dissipation. Figure 1 shows α as a function of resolution for three different P_m values in the SZ simulations. The N_x = 128 data listed in the table and plotted in the figure were taken from Simon & Hawley (2009). Note that those runs were similar to, but not the same as those done for this paper. Here, the grid zone size is equal in all directions, Δx = Δy = Δz, whereas the Simon & Hawley (2009) runs had Δx = Δz = 0.4Δy. Furthermore, the calculations done in Simon & Hawley (2009) were performed with the Roe method for the Riemann solver, in contrast to the HLLD solver used here.

Zoom In Zoom Out Reset image size

For the SZ simulations that have sustained turbulence, α appears to be converging with resolution. More specifically, by 32 grid zones per H, α is within a factor of ∼1.4 of the corresponding value at 128 zones per H. In contrast, α in the P_m = 16 and P_m = 4 runs increases by about a factor of three going from 16 to 32 zones per H. The P_m = 2 case shows a much smaller change and the turbulence dies out for all resolutions with Re = 6250 and P_m = 1 in agreement with the higher resolution simulations of Fromang et al. (2007).

Next, we explore the influence of physical dissipation terms on the FT initial field configuration, using 32 zones per H. The time history of the volume-averaged stress for these runs is displayed in Figure 2. There is a clear dependence on the dissipation parameters and on P_m in particular. For large enough resistivity (i.e., low Rm), the turbulence decays; the critical Rm value is ∼1000, in agreement with the higher resolution simulations of Simon & Hawley (2009). In Figure 3, we plot the time-averaged α values, averaged from orbit 120 to the end of the simulation, versus P_m for these runs (asterisks). Note that 32FTRm3200Pm4 and 32FTRm6400Pm4 have different Rm values, but the same P_m and nearly the same saturation level. We also plot α from the higher resolution simulations of Simon & Hawley 2009 (see their Table 3). The dashed lines are linear fits to the data in log–log space. Assuming δ in α ∝ P^δ_m we find δ = 0.54 for 32 grid zones per H, and δ = 0.33 for 128 grid zones per H; the P_m dependence is steeper at the lower resolution. Furthermore, all of the α values for the higher resolution simulations are larger than the corresponding lower resolution simulations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Both the zero-net flux and net toroidal flux results suggest that moderate resolutions, i.e., 32 grid zones per H, may be sufficient to capture the general effects of changing ν and η, at least for the range of Re, Rm, and P_m values considered here. This is not to say that everything is converged at this resolution. Indeed, Figure 3 shows a noticeable resolution effect. Full convergence likely requires a sufficiently high resolution to ensure that the effective numerical dissipation scale is below the viscous and resistive dissipation scales (see, e.g., Fromang et al. 2007; Simon et al. 2009; Simon & Hawley 2009), but the general dependence of α on dissipation parameters appears to be captured even at 32 zones per H. This conclusion is also supported by the recent results of Flaig et al. (2010), which show that ∼30 zones per H in the vertical direction may be sufficient to characterize the turbulence in their stratified shearing box simulations. This is an important point, as available computational resources limit us to this resolution in carrying out the comprehensive study of physical dissipation effects presented in this work.

Because the fastest growing MRI wavelength is proportional to v_A, the background field strength can play a significant role in the outcome of a shearing box simulation. For example, Longaretti & Lesur (2010) found that the dependence of angular momentum transport on Rm and P_m becomes steeper for weaker background (vertical) fields in unstratified shearing boxes.

Here, we consider the effects of different toroidal field strengths on the FT initial condition in the unstratified shearing box. We consider β_y = 1000 and β_y = 10000 runs both with and without physical dissipation (see Table 1 and Section 2). The time history of the volume-averaged stresses for these simulations is shown in Figure 4. For β_y = 1000, the turbulence survives at Rm = 3200 but dies at Rm = 1600, whereas for β_y = 10000, the turbulence dies at Rm = 3200. It would seem that for each increase in β_y by a factor of 10, the critical Rm value increases by roughly a factor of 2; it becomes easier to kill off the MRI with a lower resistivity as the background toroidal field is weakened.

Zoom In Zoom Out Reset image size

The time- and volume-averaged Elsasser numbers, calculated via Equation (8), are 2.3 and 137 for the β = 1000, Rm = 1600 and 3200 models, respectively. The time average is from orbit 140 to the end of the calculation, and since the turbulence decays, the average Λ for Rm = 1600 agrees with that given by the initial toroidal field value. For the β = 10000, Rm = 3200 model, the Elsasser number is 0.45, again, calculated from either the initial toroidal field strength or from Equation (8) after the turbulence has completely decayed. In both cases the initial poloidal field has β_p = 1600, corresponding to an initial poloidal Λ on order unity for these Rm values. These results are consistent with the Λ values calculated in Simon & Hawley (2009); Λ ≲ 100 leads to decay.

The time- and volume-averaged Q_y values are Q_y = 34 and Q_y = 33 for 32FTNumβ1000 and 32FTNumβ10000, respectively. The time average is done from orbit 20 to 110. The toroidal field MRI appears to be quite well resolved in both simulations. Averaging over the same time interval, we find that Q_z = 7.3 for 32FTNumβ1000 and Q_z = 6.9 for 32FTNumβ10000; the vertical field MRI is marginally resolved. Note that these Q values are calculated via the saturated state of these runs. Indeed, the initial Q values are substantially lower and either marginal or underresolved. It is not clear which Q values are a better representation of how well resolved the MRI actually is; the initial values correspond to the background field which ultimately drives the MRI, but other modes may become important in driving the MRI in the fully nonlinear state. In any case, despite the low initial Q values, sustained turbulence develops when explicit dissipation terms are not included, suggesting that at least some MRI modes are resolved. It is only when resistivity is turned on that decaying turbulence is observed. The initial low Q values may introduce uncertainty for the specific values of the critical resistivity, but the relationship between critical resistivity and initial field strength is likely to hold nevertheless.

We now turn to a series of simulations that investigate the effects of physical dissipation in vertically stratified, isothermal disks. The various combinations of runs are summarized in Table 2. Some time- and volume-averaged values of various quantities are given in Table 3. The baseline calculations without physical dissipation are done at two resolutions: 32Num which uses 32 grid zones per H and 64Num which uses 64. 32Num is run to a total integration time of 1058 orbits. It has sustained turbulence at a level of α = 0.028, calculated via Equation (18), time averaging from orbit 20 until the end of the calculation and volume averaging over all x and y and for |z| ⩽ 2H. In model 64Num, α = 0.022, where the time average is from orbit 20 until the end of the simulation at 159 orbits.

One particularly useful diagnostic is the space–time diagram of horizontally averaged quantities. Figure 5 shows space–time plots of horizontally averaged B_y and pressure-normalized total stress for a 100 orbit period in the 64Num simulation. From this, we see that B_y changes sign periodically and rises above the equatorial plane. The behavior of B_y is roughly mirror symmetric around the equator and has a period of ∼10 orbits at both resolutions. This has been seen in previous vertically stratified shearing boxes (e.g., Brandenburg et al. 1995; Stone et al. 1996; Hirose et al. 2006; Guan & Gammie 2011; Shi et al. 2010; Gressel 2010; Davis et al. 2010), with a variety of initial fields, numerical resolutions and codes, and it has even been observed in global simulations (e.g., Fromang & Nelson 2006; Dzyurkevich et al. 2010; O'Neill et al. 2010). It is apparently a generic property of MRI-driven turbulence in the presence of vertical gravity.

Zoom In Zoom Out Reset image size

The white contours in the top panel of Figure 5 indicate where the gas β value switches from greater than to less than unity. For |z| ≳ 2–2.5H, β < 1, except for some regions very near the vertical boundaries where β>1 as a result of the absence of magnetic field. The region around |z| = 2H is also where the fluid becomes completely buoyantly unstable. Using the criterion of Newcomb (1961) as outlined in Guan & Gammie (2011), the gas is buoyantly stable if

$$\begin{equation} \left|\frac{d\rho }{dz}\right| > \left|\frac{\rho ^2 g}{\gamma \! P}\right|, \end{equation} \tag{ 19 }$$

where γ = 1 because the gas is isothermal. According to this criterion applied to the simulation data, the fluid is buoyantly unstable for |z| ≳ 2H. Consistent with this, the space–time plot shows that the field rises faster for |z| ≳ 2H. For |z| ≲ 2H, there are regions of instability as well as stability, and it is within this region that the MRI is active as indicated by the presence of significant stress. Indeed, the total stress drops off rapidly near |z| ∼ 2H. It appears that the marginal buoyancy stability coupled with the MRI-induced turbulence leads to a slower rise of magnetic field until |z| ∼ 1.5–2H, where the gas then becomes completely buoyantly unstable. These results are consistent with the recent ZEUS calculations of Guan & Gammie (2011) with large radial extent as well as with Shi et al. (2010) using a version of ZEUS that includes radiation physics and total energy conservation. Compare, e.g., the top panel of Figure 5 to Figure 6 in Shi et al. (2010).

Figure 6 shows the volume-averaged toroidal field within |z| ⩽ 0.5H as a function of time. The rms field fluctuations averaged within this region are roughly a factor of ∼4 times larger; the turbulent fluctuations dominate over the mean, but not significantly so. The temporal oscillations in the mean field are apparent. The right figure shows the temporal power spectrum, revealing the dominant ∼10 orbit period. The oscillation amplitude appears to be modulated on longer timescales, ranging from tens to hundreds of orbits. Furthermore, the averaged radial field appears to exhibit the same 10 orbit cyclic behavior as the toroidal field, but with a slight temporal lag, as shown in Figure 7. The behavior of the mean radial and toroidal fields resembles the α–Ω dynamo model derived in Guan & Gammie (2011). Specifically, they write simplified evolution equations for volume-averaged field components,

$$\begin{equation} \frac{d \langle B_y\rangle }{dt} = -q \Omega \langle B_x\rangle - \frac{|v_{\rm A}|}{2H}\langle B_y\rangle + \frac{\alpha _1}{2 H} \langle B_x\rangle, \end{equation} \tag{ 20 }$$

$$\begin{equation} \frac{d \langle B_x\rangle }{dt} = - \frac{|v_{\rm A}|}{2H}\langle B_x\rangle - \frac{\alpha _2}{2 H} \langle B_y\rangle. \end{equation} \tag{ 21 }$$

The first term on the right hand side of Equation (20) is the background shear acting on the radial field. The second term is the loss of toroidal field due to buoyant rise, estimated to have a characteristic buoyant velocity equal to the toroidal field ${\rm Alfv\acute{e}n}$ speed. The third term is the α-dynamo term coupling B_x to the evolution of B_y. Equation (21) is similar except there is no shear term and the toroidal and radial field components have been flipped with respect to Equation (20). Also, in general, α₁ ≠ α₂. Note that while we follow the notation and general approach of Guan & Gammie (2011) here, Gressel (2010) has also constructed an α–Ω model, which incorporates a helicity-based dynamo mechanism and reproduces the observed behavior of the horizontally averaged field components.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Guan & Gammie (2011) numerically integrated this set of equations and found a solution that looks strikingly similar to the red and black curves in Figure 7. The value of α_1,2 sets the oscillation frequency in the mean field; α₁ = α₂ = −0.01ΩH reproduces the 10 orbit variability observed in stratified shearing boxes. Here, we numerically integrate Equation (20) using the simulation data for 〈B_x〉 (the red curve) and the initial condition for 〈B_y〉 taken from 〈B_y〉 at t = 0. We have set α₁ = 0 so that the 〈B_y〉 evolution is controlled entirely by shear and buoyancy. The result is the blue curve in the figure. The agreement between the actual evolution of 〈B_y〉 and the α–Ω dynamo model suggests that the evolution of the mean toroidal field within the mid-plane region is almost completely controlled by the shearing of radial field and the buoyant removal of the generated toroidal field.

The remaining question then is, what creates the mean radial field? In the dynamo model, the α₂ term is essential, but what mechanism is responsible for α₂ ≠ 0? The most likely candidate is MRI turbulence; turbulent fluctuations create electromotive forces (EMFs) that generate poloidal field (e.g., Brandenburg et al. 1995; Davis et al. 2010; Gressel 2010), but the details of this are still not well understood.

In Figure 8, we plot the time averaged and horizontally averaged vertical distributions for various quantities in 32Num. The time average runs from orbit 100 to the end of the calculation. The figure shows that the stress drops off rapidly near |z| ∼ 1.5H. The shape of the distribution is generally the same for both Maxwell and Reynolds stresses, with the Maxwell stress always greater than the Reynolds stress by a factor that varies from 2.3 to 5.6 depending on z; this factor is ∼4 when averaged over all z, in agreement with unstratified simulations (e.g., Hawley et al. 1995). The gas pressure maintains an approximately Gaussian distribution consistent with hydrostatic equilibrium, while the magnetic pressure is subthermal and relatively flat for all |z| ≲ 2H. The magnetic field becomes superthermal for |z|>2H, although its magnitude continues to decrease with height. These results are consistent with previous studies of isothermal disks (e.g., Stone et al. 1996; Miller & Stone 2000; Guan & Gammie 2011). Interestingly, the vertical structure of the turbulence is also consistent with local simulations containing radiation physics (Hirose et al. 2006, 2009; Krolik et al. 2007; Flaig et al. 2010) and global simulations (Hawley et al. 2001; Fromang & Nelson 2006), though we do not observe the double peak profile in the stress as seen in the simulations of Hirose et al. (2009) and Flaig et al. (2010).

Zoom In Zoom Out Reset image size

Figure 9 examines the three-dimensional structure of the magnetic field in the fully turbulent gas using streamline integration for 32Num at t = 100 orbits. The field strength is denoted by color rather than by field line density. Within |z| ≲ 2H, the field is primarily toroidal but has a smaller scale, tangled structure in the x and z directions. Very near the vertical boundaries, however, the field appears to develop larger poloidal components. Utilizing snapshots taken throughout the evolution of 32Num, we find that this is typical of the saturated state, except for t = 550 orbits, in which the field near the boundaries is primarily vertical. Figure 16 of Hirose et al. (2006) shows the field in a shearing box with a vertical domain twice as large as 32Num. In their simulation, the field at |z| = 4H is mainly toroidal while the field structure at |z| = 4H in 32Num resembles the field structure near |z| = 8H in their simulation. This suggests that the vertical outflow boundary conditions are influencing the field structure very near the boundaries. Away from the vertical boundaries, however, the magnetic field in 32Num appears to have a very similar structure to that in Hirose et al. (2006).

Zoom In Zoom Out Reset image size

Having established the baseline simulations without physical dissipation, we now turn to the main focus of our present work: the effect of ν and η on vertically stratified MRI turbulence.

The Maxwell and Reynolds stress time evolution for the simulations performed with 64 zones per H is shown in the left panel of Figure 10. Runs 64Rm3200Pm0.5, 64Rm800Pm0.5, and 64Rm3200Pm4 have decreasing turbulence levels, while turbulence is sustained at a statistically constant value in the remaining simulations. Furthermore, 64Rm800Pm0.5 undergoes alternating periods of low and high stress, though the overall trend is downward with time.

Zoom In Zoom Out Reset image size

The right panel of Figure 10 is the stress evolution for the equivalent simulations with 32 zones per H, shown for a much longer time period than the 64 zone runs. There is considerable variability on long timescales. 32Rm3200Pm0.5 in particular exhibits alternating periods of low and high stress levels, occurring on ∼100–200 orbit timescales. The more viscous and resistive run, 32Rm800Pm0.5, shows similar variability but on a much shorter timescale of ≳10 orbits. Turbulence in 32Rm3200Pm2 appears to have leveled off at a smaller value without any indication of regrowth to the higher state. Both 32Rm6400Pm4 and 32Rm3200Pm4 remain at relatively high stress levels, which are very similar between the two runs, likely a result of having the same P_m.

Turning to the high and low states in 32Rm3200Pm0.5, the left panel of Figure 11 shows the vertical profile for the horizontally averaged total stress during a high turbulence state (solid line, time averaged from 500 to 570 orbits) and a low turbulence state (dashed line, time averaged from 700 to 770 orbits). The high state has considerable stress within the mid-plane region (|z| ≲ 2H), whereas in the low state, the stress is peaked near |z| ∼ 2H. Note that there is still a nonzero stress in the mid-plane region even during the low state. The Reynolds stress is nonzero and relatively flat throughout the mid-plane, whereas the Maxwell stress drops off dramatically at the mid-plane. What Maxwell stress is present in the low states seems to derive from the presence of residual radial and toroidal field. (This effect was also noted by Turner et al. 2007 in the dead zone regions of their shearing box calculations.) One characteristic of MRI turbulence is the ratio of the Maxwell stress to the magnetic pressure, which is typically on order 0.4 (Hawley et al. 1995; Blackman et al. 2008; Table 3). The vertical profile of this ratio is shown in the right panel of Figure 11. This ratio is between 0.3 and 0.4 for |z| ≲ 1.5H in the high state and near |z| ∼ 2H in the low state; in the mid-plane region during the low state, the ratio is much smaller and approaches zero.

Zoom In Zoom Out Reset image size

These differences between the high and low states are seen in the other simulations that exhibit this variability (except for 32Rm800Pm0.5, in which the variability occurs on such a short timescale that temporal averaging becomes difficult). In summary, during the high state, the MRI is fully active, producing turbulent stresses for |z| ≲ 2H. In the low state, the turbulence in the mid-plane region has died down, leaving weaker Reynolds stresses. MRI-driven turbulence is present, however, in narrow regions near |z| ∼ 2H. This is similar to the behavior seen in Fleming & Stone (2003) where active layers above and below the mid-plane drive Reynolds stresses in the equatorial dead zone.

Is the variability between high and low turbulence an artifact of using a relatively low resolution in these calculations? Several pieces of evidence suggest that this is not the case. First of all, both resolutions for Re = 1600, P_m = 0.5 show the same variable stress behavior on ≳10 orbit timescales. Second, even at 32 zones per H, dissipation coefficients play a significant role in determining the stress level, as shown above. The fact that the low-resolution, unstratified simulations show constant saturation levels for sufficiently small η, whereas vertically stratified simulations exhibit this variability for the same parameters, suggests that the variability is a direct result of adding in vertical gravity. Third, this variability was seen in the simulations of Davis et al. (2010), which were run at a higher resolution of 64 zones per H. We note, however, that while Davis et al. (2010) used the same numerical algorithm as in this work, they used a different initial magnetic field configuration and vertical boundary conditions.

Finally, we examine the Q values (Equation (10)) using an ${\rm Alfv\acute{e}n}$ speed defined by Equation (9) and a volume average for all x and y and for |z| ⩽ 0.5H. Figure 12 shows Q_y and Q_z as a function of time for 32Rm6400Pm4 and 64Rm6400pm4. The results suggest that the toroidal field MRI is quite well resolved, but that the vertical field MRI may be only marginally resolved for the lower resolution simulation. The higher resolution Q values are roughly a factor of two larger than the lower resolution Q values, which is simply a result of the change in the grid zone size; the turbulent saturation level is roughly the same between the two resolutions. This result, coupled with the somewhat low Q_z value, suggests that the vertical field MRI may not be playing a particularly significant role in setting the saturation level in these simulations.

Zoom In Zoom Out Reset image size

If this high–low variability is indeed a physical effect, what is its origin? First, consider the space–time diagram of the horizontally averaged B_x and B_y for several simulations. Figure 13 shows the first 200 orbits of 32Rm800Pm0.5. The turbulence level decreases dramatically from the beginning (see also Figure 10). This is not too surprising considering that the same Rm, P_m values cause rapid decay of the turbulence in the unstratified case; the resistivity is large enough to damp the MRI. The space–time plots show that after this decay, there is a residual magnetic field left within the mid-plane region of the disk. In particular, within |z| ≲ 0.5H, there is a net horizontally averaged B_x < 0 and B_y < 0 around t = 110 orbits. The average B_x within this region remains constant for awhile, and B_y increases due to the shear of B_x, eventually flipping to B_y>0. By t ∼ 130 orbits, the turbulence has re-emerged, and the average B_y rises to larger |z|. The resistivity then kills off the MRI again, leading to another period of B_x shearing into B_y before the next period of increased turbulence.

Zoom In Zoom Out Reset image size

Model 32Rm3200Pm0.5 also experiences alternating high and low turbulence states, but it is not immediately clear why the turbulence within the mid-plane region should decay since the unstratified run with the same dissipation terms has sustained turbulence. As noted earlier and in Simon & Hawley (2009), the critical Rm value below which the turbulence decays in unstratified shearing boxes is ∼1000, but here Rm = 3200. The major difference between the stratified and unstratified simulations is that with vertical gravity, the net magnetic field within a localized region of the domain can change due to buoyancy, whereas with unstratified boxes the net toroidal field remains constant. Loss of flux can raise the critical Rm value (numerical resolution is also likely to have an effect); indeed, the unstratified simulations initialized with β = 1000 and β = 10000 fields have demonstrated that the critical Rm depends upon the background field strength. The loss of net flux in a stratified box could similarly change the critical Rm value.

In this simulation, the average toroidal field within the mid-plane region, 〈B_y〉, oscillates around zero with a period of 10 orbits. Thus, every 10 orbits or so, 〈B_y〉 is conceivably weak enough for resistivity to kill the turbulence. However, the turbulence remains for many of these 10 orbit periods. Furthermore, averaging B_y within some vertical distance from the mid-plane erases information about the field structure there, e.g., 〈B_y〉 might be small but there could still be strong toroidal fields of opposite polarity close to z = 0. The point is one cannot necessarily expect the turbulence to decay away strictly whenever 〈B_y〉 drops below a certain (small) value.

The 〈B_y〉 oscillation amplitude appears to be modulated by a longer timescale, more on the order of ∼100 orbits (see, e.g., Figure 6). This behavior is present in all simulations with and without physical dissipation. This is also roughly the same timescale over which the system switches between low and high states in the Rm = 3200 simulations. Comparing the time evolution of 〈B_y〉 for |z| ⩽ 0.5H with the evolution of the total stress shows that the minima in the oscillation amplitude are generally correlated with the decay of mid-plane turbulence. One exception is near 200 orbits in 32Rm3200Pm0.5 in which 〈B_y〉 becomes rather small, but the turbulence remains active, though relatively weak compared to the fully active state. This correlation implies that if the mean toroidal field near the mid-plane remains sufficiently small for some time, resistivity can overwhelm the MRI and cause decay.

Evidently, there exists a critical Rm below which the turbulence experiences long-timescale variability; this critical value is Rm < 6000. We carry out two more stratified simulations with Rm ∼ 6000 but different P_m values in order to further test this hypothesis. The first simulation is 32Rm6400Pm8; thus, the Rm value is 6400, and P_m is relatively large. The mid-plane turbulence is sustained over a long time, nearly 330 orbits, without any sign of decay. The dissipation coefficients of the second simulation, 32Rm6250Pm1, are chosen to match the relatively high Rm, low P_m simulation that decays in the zero-net flux shearing box (see Fromang et al. 2007; Simon & Hawley 2009). This simulation also remains in the high state for nearly 330 orbits.

If the dissipation is large enough to cause MRI-driven turbulence to decay, what leads to its reactivation after ≳100 orbits of time? Figure 14 shows the space–time plot of B_x and B_y for a 300 orbit period in 32Rm3200Pm0.5 during which the mid-plane turbulence dies out and eventually returns. For clarity, we also plot the volume-averaged horizontal field, 〈B_x〉 and 〈B_y〉, where the average is done for all x and y and for |z| ⩽ 0.5H. During the low state, there is a small net radial field that remains in the mid-plane region and generates B_y through shear. The sign of 〈B_x〉 reverses and with it the sign of ∂〈B_y〉/∂t. The last such switch in the low state occurs around 750 orbits after which 〈B_y〉 continues to grow up until 790 orbits. At this point, the 10 orbit period oscillations resume.

Zoom In Zoom Out Reset image size

Figures 13 and 14 suggest that it is the growth of B_y due to shear that periodically reactivates the MRI in the mid-plane. Indeed, during the low states, the poloidal fields are sufficiently weak that the most unstable wavelengths of the radial and vertical field MRI are very underresolved; Q_x ≲ 1 and Q_z ≲ 1 for both 32Rm800Pm0.5 and 32Rm3200Pm0.5. (Q was calculated as a function of time using Equations (9) and (10).) Considering the toroidal field, however, we find that Q_y is well above the marginal resolution limit when the mid-plane turbulence starts to decay. In the low states of 32Rm800Pm0.5, Q_y ∼ 10–20, occasionally dropping to Q_y = 6. The typical Q_y values for the low states of 32Rm3200Pm0.5 are similar but somewhat smaller. Of course, when studying the behavior of the MRI in a situation where the Q values are small, one is in a regime where numerical resolution is likely to matter. The specific behavior of the disk in the low state might be different at higher resolution, but shearing of radial into toroidal field is itself not so dependent on resolution. Thus, Q may well play a role in determining when the MRI is reactivated, but not in how it is reactivated.

As a further test, we carry out two additional experiments. First, we take the state of the gas in 32Rm3200Pm2 at t = 150 orbits when the stress levels are decreasing and the average of B_y within |z| ⩽ 0.5H, 〈B_y〉 = 5.9 × 10⁻⁶, is relatively small compared to the oscillation amplitude of 〈B_y〉 in the high state, which is ∼5 × 10⁻⁵. We restart this simulation and add a net B_y = 8.9 × 10⁻⁵ into the region |z| ⩽ 0.5H, which corresponds to a toroidal β ≈ 126 (using β defined with the initial mid-plane gas pressure P_o). This run is 32Rm3200Pm2_By+ (see Table 2). Figure 15 shows the subsequent evolution of the stress along with the stress evolution of 32Rm3200Pm2. Not only does the mid-plane turbulence return, but the system undergoes episodic transitions between low and high states on ∼100 orbit timescales, as in 32Rm3200Pm0.5.

Zoom In Zoom Out Reset image size

In a second experiment, we initialize a stratified shearing box with all the same parameters as in 32Num but with an initially very weak radial magnetic field. Specifically, for |z| ⩽ 0.5H, $B_x = -\sqrt{2 P_o/\beta _x}$, where β_x = 10⁶. This field strength is very underresolved; the Q_x value is 0.2 and the radial field MRI will not be significant. Figure 16 shows the space–time diagrams of horizontally averaged B_x and B_y. The weak radial field is sheared into toroidal field that grows linearly in time until it reaches a sufficient strength to activate the MRI. The onset of MRI turbulence is rapid, and once the MRI sets in, the subsequent behavior is very similar to the other vertically stratified MRI simulations; there are rising magnetic field structures, dominated by the toroidal component, and the period of oscillations in the mean field is ∼10 orbits. Note that this experiment is not an exact imitation of the low state in our simulations; there is no MRI activity near |z| ∼ 2H initially. However, this calculation demonstrates that shear amplification of a weak radial field can eventually lead to turbulence and move the system to the high state.

Zoom In Zoom Out Reset image size

All of the simulations in which turbulence in the mid-plane sets in after a period of quiescence show the presence of a net radial field in the mid-plane during the low state. 32Rm3200Pm2, however, is the only simulation that does not show the re-emergence of the mid-plane MRI, despite over 1000 orbits of integration. An examination of the mid-plane region (up to |z| ≲ H) in the low state of this run shows that the residual radial field is weaker than in the low states of the other simulations. If this radial field would remain constant in time, the toroidal field would continually strengthen to the point of reactivating the MRI. However, 〈B_x〉 continues to change sign even in the absence of turbulence, though with a period of many hundreds of orbits. That is, 〈B_x〉 oscillates about zero but with a very small amplitude. 〈B_y〉 similarly oscillates around zero and never reaches a sufficient amplitude to reactivate the MRI, and the simulation remains in the low state.

In fact, 〈B_x〉 oscillates about zero in all of our simulations, even in the low states (see, e.g., Figure 14). The various space–time diagrams suggest that the horizontal field migrates toward the mid-plane region from near |z| ∼ 2H where MRI activity is ongoing. The oscillations in the mid-plane field components are a direct result of the oscillations previously observed in vertically stratified MRI turbulence. It is not entirely clear, however, what causes the field to migrate to the mid-plane. Does it diffuse there or is it carried downward by some sort of large-scale flow? The shearing box calculations of Turner et al. (2007) and Turner & Sano (2008), which include Ohmic resistivity via a treatment of ionization chemistry, suggest that turbulent diffusion from the active layers is responsible for the field migration. While our simulations include a simpler prescription for resistivity, it is conceivable that a similar process is at work here. These issues will be addressed in future work.

To summarize the results for transitions between high and low states, we find that for sufficiently small Rm, MRI turbulence within the mid-plane region can decay away, in part because of the loss of net flux due to buoyancy. Residual net radial field remains and subsequently creates toroidal field from shear. Once this toroidal field reaches a sufficiently large amplitude, the MRI is reactivated, and the turbulence is sustained for some duration until it decays again, repeating the pattern. In other words, the α–Ω dynamo, seen in stratified simulations without resistivity, continues to operate in resistive simulations where the mid-plane MRI is suppressed. The dynamo is key to the reactivation of MRI-driven turbulence.

This behavior appears to be independent of the P_m values we have used, except for near Rm = 3200; 32Rm3200Pm4 remains in the high state, whereas 32Rm3200Pm2, 32Rm3200Pm2_By+, and 32Rm3200Pm0.5 do not. However, in the higher resolution runs, 64Rm3200Pm4 and 64Rm3200Pm0.5 show decay but 64Rm3200Pm2 appears to be sustained (though again, these simulations were not integrated very far). While P_m may play a role here, the line between sustained and highly variable stress levels is unlikely to be hard and fast, and many factors probably contribute to the nature of the turbulence near this Rm value.

Comparing our results with those of Davis et al. (2010), we note that the largest Rm used in their simulations was Rm = 3200, consistent with the largest critical Rm observed in our simulations. Second, an examination of the space–time data from their simulation with Rm = 1600, P_m = 2 (kindly provided by the authors) shows the same behavior as we have observed here; a net radial field remains within the mid-plane region after decay, shearing into toroidal field, from which the MRI is reactivated within this region.

Interestingly, Turner et al. (2007) briefly note that the toroidal field in the dead zone of some of their highly resistive simulations can occasionally become strong enough to reactivate turbulence there, after which the MRI is rapidly switched off again due to Ohmic dissipation. These simulations are consistent with our highest resistivity (Rm = 800) simulations, in which turbulence occurs in short outbursts rather than through a transition to a long-lived high state, as is the case at lower resistivities.

Lastly, we examine the topology of the magnetic field in the shearing box during the low state. Figure 17 shows the equivalent information as Figure 9, but for orbit 550 of 32Rm3200Pm2_By+. For |z| ≳ 2H, the field remains mainly toroidal but with noticeable excursions from being purely toroidal, resembling the field structure in this region during the high state. Within 2H, the field is almost completely toroidal, and any small poloidal field present within this region is not visible in this image. We also examined the azimuthally averaged poloidal field structure in several snapshots of this run. We found that the structure of the field was different, depending on which snapshot we examined. At some times, the poloidal field within 2H is almost completely radial, with very little vertical field. At other times, the vertical and radial fields are comparable in size such that the field takes on a more loop-like structure.

Zoom In Zoom Out Reset image size

In this section, we investigate how P_m affects angular momentum transport in simulations that do not exhibit the high–low variability. Figure 18 shows the volume-averaged stress levels for these simulations normalized by the gas pressure. The volume average is done for all x and y and for |z| ⩽ 2H. The figure shows a general increase in the turbulence levels with P_m, but there is also significant time variability in the stress. As a result, the curves overlap at certain times, much more so than for unstratified turbulence (see Figure 2).

Zoom In Zoom Out Reset image size

Figure 19 displays the time-averaged α values as a function of P_m for the unstratified (left panel) and stratified simulations (right panel). The time average for the unstratified simulations is the same as in Figure 3, from orbit 120 to the end of the calculation, and for the stratified simulations, this average is done from orbit 150 until the end of the simulation. The error bars denote one standard deviation about the temporal average of the numerator in α.³ While there is a clear dependence of α on P_m in the stratified simulations, it is not as steep as in the unstratified simulations. Taking a linear fit in log–log space and assuming α ∝ P^δ_m, δ = 0.54 for the unstratified runs (from Section 3) and δ = 0.25 for the stratified calculations. In addition, the time variability is significantly larger in the stratified runs than in the unstratified case.

Zoom In Zoom Out Reset image size

Figure 20 is the vertical profile of the time averaged and horizontally averaged total stress for the sustained turbulence runs. The time average was done from orbit 150 to the end of each simulation. Increasing P_m increases the stress at nearly all z, and in all cases the stress drops off dramatically above |z| ∼ 1.5H, consistent with the baseline simulations. Run 32Rm6400Pm8 appears to have a sharper peak in the stress profile near z = 0, compared to a flatter stress profile for |z| ≲ 1.5H in the other simulations. Using different temporal averaging windows for the averaged stress profiles produces the same general result; the stress increases with P_m and 32Rm6400Pm8 has a sharper peak near the mid-plane. In some cases, however, the stress does not necessarily increase monotonically with P_m at |z| ≳ 1.5H. Finally, we examined the vertical profile of the same quantities as in Figure 8. We found that these profiles in the sustained turbulence simulations are all very similar to each other and to 32Num; the general vertical structure of the disk does not appear to be sensitive to P_m.

Zoom In Zoom Out Reset image size

In summary, where the MRI is operating within a stratified disk, increasing P_m leads to an increase in stress, consistent with the behavior seen in unstratified simulations.

The algorithms for solving the shearing box Equations (1)–(3) and the specific implementation of those algorithms within the Athena code are described in Stone & Gardiner (2010); here we provide a brief summary of the numerical methods detailed in that paper.

The lefthand sides of the equations are solved via the standard Athena CTU flux-conservative algorithm (Stone et al. 2008). The gravitational and Coriolis source terms in the momentum equation are evolved via a combination of an unsplit method consistent with CTU and Crank–Nicholson differencing constructed to conserve energy exactly within epicyclic motion (Gardiner & Stone 2005; Stone & Gardiner 2010).

For our simulations in which the radial size of the shearing box is a scale height, H, the velocity is initialized with a background shear flow given by

$$\begin{equation} v_y = -q \Omega x. \end{equation} \tag{ A1 }$$

For sufficiently large x domains, this velocity can become supersonic (orbital velocities are, in general, supersonic in accretion disks). For large x domains, the Courant limit on the timestep can become significant. Furthermore, the background shear flow can lead to a systematic change in truncation error with radial position in the box, which in turn introduces purely numerical features in the radial density and stress profiles (Johnson et al. 2008). Hence, for large radial domain simulations we utilize an orbital advection scheme, which subtracts off this background shear flow and evolves it separately from the fluctuations in the fluid quantities (Masset 2000; Johnson et al. 2008; Stone & Gardiner 2010).

In the shearing box, the y boundary condition is periodic. The x boundary condition is shearing periodic (Hawley et al. 1995; Stone & Gardiner 2010): quantities are reconstructed in the ghost zones from appropriate zones in the physical domain that have been shifted along y to account for the relative shear from one side of the box to the other. In Athena, this reconstruction step is performed on the fluid fluxes, and then the ghost zone fluid variables are updated via these reconstructed fluxes. The order of this reconstruction matches the spatial reconstruction in the physical grid, e.g., third-order reconstruction of the ghost zone fluxes is done when the PPM spatial reconstruction is employed. Furthermore, the y momentum is adjusted to account for the shear across the x boundaries as fluid moves out one boundary and enters at the other. The y component of the EMF is reconstructed at the radial boundaries to ensure precise conservation of the net vertical magnetic flux (Stone & Gardiner 2010).

In unstratified shearing box simulations the z boundary is periodic. This same boundary condition has been employed in stratified boxes as well (e.g., Davis et al. 2010). For our simulations, we use an outflow boundary condition instead. The density ρ is extrapolated into the ghost zones based upon an isothermal, hydrostatic equilibrium, using the last physical zone as a reference. Therefore, for the upper vertical boundary, the ρ value in grid cell k is

$$\begin{equation} \rho (k) = \rho (ke) {\rm exp}\left(-\frac{z(k)^2-z(ke)^2}{H^2}\right), \end{equation} \tag{ A2 }$$

where ke refers to the last physical zone at the upper boundary. An equivalent expression holds for the lower vertical boundary. This extrapolation provides hydrostatic support against the opposing vertical gravitational forces, which are also applied in the ghost zones. All velocity components, B_x, and B_y are copied into the ghost zones from the last physical zone assuming a zero slope extrapolation. If the sign of v_z in the last physical zone is such that an inward flow into the grid is present, v_z is set to zero in the ghost zones. Finally, the ghost zone values of B_z are calculated from B_x and B_y to ensure that ${\bf \nabla } \cdot {\mbox{\boldmath {$B$}}} = 0$ within the ghost zones.

Viscosity and resistivity are added via operator splitting; the fluid variables updated from the main integration (Appendix A.1) are used to calculate the viscous and resistive terms. The viscosity term is calculated via the divergence of the viscous stress tensor, Equation (4), and the resistive term is included as an additional EMF within the induction Equation (3). This formulation allows us to discretize the viscous and resistive terms in a flux-conservative and constrained-transport manner, consistent with the Athena algorithm. Note that the resistive contribution to the y EMF must also be reconstructed at the shearing-periodic boundaries in order to preserve B_z precisely.

The addition of viscosity and resistivity places an additional constraint on the timestep,

$$\begin{equation} \Delta t = {\rm MIN}\left(\Delta t_{\rm CTU}, C_o \frac{\Delta ^2}{4 \nu }, C_o \frac{\Delta ^2}{4 \eta }\right), \end{equation} \tag{ A3 }$$

where Δt_CTU is the timestep limit from the main integration algorithm (see Stone et al. 2008), C_o = 0.4 is the CFL number, and Δ is the minimum grid spacing, Δ = MIN(Δx, Δy, Δz). Note that in most of our simulations ν and η are sufficiently small so that they do not limit the timestep.