What Is The Numerically Converged Amplitude of MHD Turbulence in Stratified Shearing Boxes?
aa r X i v : . [ a s t r o - ph . H E ] D ec What Is The Numerically Converged Amplitude of MHDTurbulence in Stratified Shearing Boxes?
Jiming Shi and Julian H. KrolikDepartment of Physics and Astronomy, Johns Hopkins University, Baltimore, MD 21218andShigenobu HiroseInstitute for Research on Earth Evolution, JAMSTEC, 3173-25 Showamachi,Kanazawa-ku, Yokohama, Kanagawa 263-0001, JapanReceived ; accepted 2 –
ABSTRACT
We study the properties of the turbulence driven by the magnetorotationalinstability (MRI) in a stratified shearing box with outflow boundary conditionsand an equation of state determined by self-consistent dissipation and radiationlosses. A series of simulations with increasing resolution are performed withina fixed computational box. We achieve numerical convergence with respect toradial and azimuthal resolution. As vertical resolution is improved, the ratio ofstress to pressure increases slowly, but the absolute levels of both the stress andthe pressure increase noticeably. These results are in contrast with those of pre-vious work on unstratified shearing boxes, in which improved resolution causeda diminution in the magnetic field strength. We argue that the persistence ofstrong magnetic field at higher resolution found in the stratified case is due tobuoyancy. In addition, we find that the time-averaged vertical correlation lengthof the magnetic field near the disk midplane is ≃ Subject headings: accretion, accretion disks — MHD — radiative transfer — methods:numerical
1. INTRODUCTION
The physical mechanism for transfer of angular momentum through disk gas is believedto be the magnetorotational instability (MRI) (Balbus & Hawley 1991). The MRI is alocal linear instability, driven by exchange of angular momentum along magnetic field linesthreading material at different distances from the central object. Its growth rate is of orderthe orbital angular frequency, and the fastest-growing wavelength is approximately thepropagation length of Alfv´en waves in one orbital period (Balbus & Hawley 1998). Butwhat mechanisms limit the exponential growth of the MRI and determine the saturationamplitude of the MRI-driven turbulence is still a big problem to be resolved.In the context of hydrodynamic turbulence, small scale dynamics do not significantlyaffect large scale behavior. In magnetohydrodynamics, however, and in particular in thesomewhat special circumstances of the anisotropic turbulence expected in accretion disksdue to their characteristic orbital shear, this may not be the case. Given the inadequacyof analytic methods, numerical simulations are required for quantitative studies. To reachthe smallest possible length scales, the best approach is not to study an entire disk, butonly an annulus of limited azimuthal extent. A “shearing box” approximation, in whichsuch an annulus is stretched into rectangular geometry, periodic boundary conditions areapplied along the azimuthal direction, and sheared-periodic boundary conditions appliedalong the radial direction, works well (Hawley et al. 1995). Recent unstratified shearingbox simulations (without a vertical component of gravity and imposing periodic boundaryconditions in the vertical direction) have shown that the saturation level of the magneticstress appears to converge to zero as the grid cell width becomes infinitesimal when thereis no net vertical flux (Fromang & Papaloizou 2007). In these simulations, the correlationlength of the turbulent magnetic field also became smaller and smaller as the grid scaleshrank. Only if some small but non-zero dissipation is included (Lesur & Longaretti 2007; 4 –Fromang et al. 2007), is convergence achieved, with the field having some small finite value.This result has led to considerable discussion and puzzlement. One possible explanationis that it is due to an unphysical approximation. Because the vertical gravity is proportionalto altitude from the midplane, it has been thought that unstratified shearing boxes area good model for the midplane region of the disk, where gravity should be very weak.However, eliminating gravity also eliminates a physical length scale, the vertical pressurescale height. The absence of any physical length scale other than the grid scale mayexplain the decay of the correlation length with diminishing grid scale, and the decay of themagnetic field strength as well. In this paper, we wish to test whether the dependence ofmagnetic stress on grid scale changes when vertical gravity, along with the physical lengthscale it introduces, are included. We will also add a further degree of realism by explicitlycomputing the temperature by balancing dissipation with radiative cooling.We are not the first to explore MRI behavior in a stratified shearing box; the shearingbox with vertical stratification has been studied and developed since the 1990s. Earlystratified disk study has already shown very different vertical structures compared withthe unstratified model (Brandenburg et al. 1995; Stone et al. 1996). Even though thesesimulations extended only a few scaleheights from the midplane, where the structure didresemble that of an unstratified homogeneous disk, at high altitudes, the gas density wasmuch lower and the magnetic field dominated the pressure. In those simulations, a simpleequation of state, either isothermal or adiabatic, was adopted. No cooling, or only simpleoptically-thin thermal relaxation, was considered. Due to the numerical difficulties ofdefining an outflow boundary above an approximately hydrostatic fluid, less realistic verticalboundary conditions were often used: for instance, stress-free and reflecting boundaryconditions were used in Brandenburg et al. (1995), and outflow boundaries with extraresistive layers below them in Stone et al. (1996). By including the displacement current in 5 –the equation of motion, Miller & Stone (2000) introduced an Alfv´en speed limiter in theirsimulations in order to lengthen the time-step, providing a way to more computationallyefficient long-term simulations. Again, the equation of state was a simple isothermal one,and no radiative transfer effects were studied. Using flux-limited diffusion (FLD) to solvethe radiation transfer problem, Turner (2004) performed the first illustrative calculation ofa vertically stratified disk segment that included dissipation and radiation effects. However,in this calculation, energy was not completely conserved: only magnetic energy losses werecaptured into heat. Full energy conservation was first achieved by Hirose et al. (2006), whosimulated a gas pressure dominated disk annulus. To smooth the field when it crosses theoutflow boundaries, they added a small amount of artificial resistivity into the ghost cells.They also applied the FLD approximation to describe the radiative transfer within thedisk. Recently, vertically stratified disk segments with both comparable radiation and gaspressures (Krolik et al. 2007; Blaes et al. 2007) and radiation pressure much larger than gaspressure (Hirose et al. 2009) have also been studied. In these papers, a similar techniquewith thin diffusive layers extended into the problem volume near the top and bottomboundaries was tested and implemented. In all these stratified simulations, significantcontrasts between stratified and unstratified disks were found. For instance, magneticbuoyancy leads to a highly magnetized ‘corona’ (e.g., Miller & Stone 2000; Turner 2004;Blaes et al. 2007) that is completely absent in unstratified disks.As already mentioned, it is the objective of this paper to study numerical convergenceof a density-stratified shearing box in which energy is conserved and radiation transfer istaken into account. To do so, we performed a set of numerical simulations with increasingresolution, but we did not include any explicit diffusivity other than the small amount nearthe boundaries and the von Neumann-Richtmyer bulk viscosity in compressive regions tothermalize kinetic energy in shocks. 6 –The paper is organized as follows. In §
2, we describe the calculational method andinitial setup of the simulations. The main results are presented in §
3. In §
2. CALCULATIONAL METHOD2.1. BASIC MODEL
We adopt the code described in Hirose et al. (2006) as our basic tool for the calculation.The complete updated (with Compton scattering included) equation set can be found in § . v y = − (3 / x in the background. We carefully includethe Coriolis force, gravitational tidal forces and the vertical component of the gravity inthe momentum equation. As the magnetic Reynolds number is usually large in accretiondisks, the ideal MHD limit is a simple but reasonable choice to describe the magnetic field.The radiation field is described by the FLD approximation. For simplicity, the opacityis thermally averaged. The gas and radiation exchange both momentum and energy viaThomson scattering and free-free absorption. Energy exchange via Compton scatteringis also included although it contributes little under the conditions we examine here. Tocomplete the equation set, we assume an adiabatic index γ = 5 / − times the initial midplane density; our energy floor and velocity capare the same as in Hirose et al. (2006). As discussed in previous work (e.g. Hirose et al.2006, 2009), the artificial energy injection due to those three limiting values is negligiblecompared to the total energy content. 7 – As a test case, we chose to repeat the physical problem first investigated byHirose et al. (2006), a stratified shearing box in which, averaged over 50 orbits ofwell-developed turbulence, the gas pressure was approximately five times greater than theradiation pressure. Our primary reason for using this as our test-case is that it is the onlyexample in the published literature of a gas-dominated stratified shearing box studied withself-consistent thermodynamics. When the thermal state of such a shearing box is foundself-consistently, it is fully characterized by only two parameters: the rotation rate at thecenter of the box Ω = 5 .
90 s − and the surface density Σ = 9 . × g cm − , producingan electron scattering optical thickness 3 . × . In terms of a Shakura-Sunyaev model,those parameters correspond to a disk annulus at a radial distance r = 300 r g ( r g = GM/c is the gravitational radius) from a central black hole with mass M = 6 . M ⊙ , an accretionrate that would yield a total luminosity of 0 . α = 0 .
03. The predicted effectivetemperature is T = 4 . × K. In the gas pressure-dominated limit, the scale height(half-thickness) of the disk is H = 3 . × cm. We choose the characteristic scale height H as our length unit.Our inital condition is similar to the one used in Hirose et al. (2006) except that weassumed a dissipation profile dF/dz ∝ d Σ / /dz instead of ∝ d ln Σ /dz . This form forthe dissipation profile is more similar to the time-averaged dissipation profile found inprevious simulations, so it helps the simulation to pass the transient phase quickly, butdoes not affect the later stages. The initial disk is in approximate hydrostatic and radiativeequilibrium below the photosphere. Outside it, the flux is constant, and the gas density isset to the density floor value. The initial configuration of the magnetic field is a twistedazimuthal flux tube of cicular cross section, centered at x = z = 0 and having a radius of 8 –0 . H . The tube has uniform interior field strength B = 2 . × G, corresponding to4% of the initial box-averaged gas pressure plus radiation pressure. The maximum poloidalcomponent of the field strengh is 5 . × G, while the net vertical flux is zero. Althoughthis initial magnetic field has both net azimuthal flux and net helicity, neither quantity ispreserved, and therefore neither provides any long-term constraint on the evolution of themagnetic field. The fastest-growing vertical MRI wavelength in the midplane in the initialstate is 1 . × cm, which is resolved with 8 . c s ≡ p (4 E/ γp ) /ρ , where E is the radiation energy, p is the thermal pressure and ρ is the gas density.All simulations in this paper start with the same computational box, which extends 2 H in the radial direction, 8 H along the orbit, and 8 H on either side of the midplane. We firstperformed a standard run, denoted as STD32, with moderate resolution: 32 × ×
256 cells( x × y × z ) with constant cell size ∆ x = ∆ z = H/
16 = ∆ y/
2. Next we carried out threeruns that had doubled resolution in only one of the three directions; they are labelled: X64for 64 × × × ×
256 and Z512 for 32 × × . × min (∆ x , ∆ y , ∆ z ) / ∆ t , to zero at 1 H into thecomputational domain, i.e. 16 cells into the problem volume for STD32, X64 and Y128,and 32 cells for Z512 and DBLE. Here ∆ t is the time step in simulations.The box height is designed to be large enough to contain most of the gas for long-termevolutions. However, the surface mass density of the disk segment can vary due to matteradded by the density floor or, more importantly here, gas outflow. In this paper, wepermitted variations of the surface mass density only up to 2 . .
25 to prevent outflow when the matter swells .When doing so, we kept the resized box at the same resolution as before, i.e., ∆ x ,∆ y and∆ z were unchanged. When the new cells were initialized, they were given density andenergy equal to the floor values. The velocities were set to zero except for the backgroundshear. The transverse components of the magnetic field were copied from the values inthe vertically-aligned cells on the old boundaries , while the perpendicular component wascalculated by enforcing the divergence-free constraint. The layers with artificial resistivitywere resized by the same factor of 1.25. We configured the resizing this way to keepthe surface mass density of the box nearly constant while not introducing large pressuregradients or energy injection in the extended zones (the added energy was <
1% of thetotal). We also tested whether this size change alone causes any noticeable effects on theresults. A test run was started at t = 90 orbits of STD32 with its height enlarged by a factorof 1 .
25 and evolved 50 more orbits. It was then compared with the data from t = 90 to 140orbits in STD32. They are statistically similar: the time- and volume-averaged stresses and Once the loss of surface density becomes greater than 2 . ≤
3% from the standard. Note that there is no significantgas outflow in STD32 during this test period, so the box resizing must be responsible forany changes. In all runs, we chose a duration long enough that the staturation level appearsto be quasi-steady. Depending on the simulation, this required between ∼
120 and 300orbits, which amount to ∼
20 –40 cooling times.
3. RESULTS
We are primarily interested in time-averaged quantities representative of conditions ina steady-state disk. To determine the appropriate time-averaging period, we examined howrapidly the effects of our initial conditions and transient response die away. As Figures 1,2, and 3 show, in every simulation the various contributions to the stress reach levelscharacteristic of the statistical steady-state by ∼
10 orbits from the beginning. Topologicalproperties of the initial magnetic field are erased rapidly: the net azimuthal flux changessign on ∼ ∼ We show the time evolution of the stresses and energies of STD32 in Fig.1. STD32has the same resolution and nearly identical initial setup as the run in Hirose et al. (2006),which makes it easy to compare them statistically (note, however, that STD32 ran for fivetimes as long). STD32 appears to have two stages after the transient decay, which separatearound 130 orbits. During the first stage, the maximum Maxwell stress is about 3 timesas great as the minimum, which is nearly the same as the gas pressure-dominated run in 11 –Hirose et al. (2006). The ratio of box-integrated stress to box-integrated total pressure(including both the gas and radiation pressures) is 0 . .
016 in their paper). In the second stage, thepeak value of the stress doubles, and the range of the fluctuations is roughly a factor of 4.A quick look at the radiation energy and gas energy plots in Fig.1 reveals that although thedisk starts as a gas pressure-dominated system, it gradually evolves toward a situation witha larger ratio of radiation to gas pressure: the time- and volume-averaged ratio hh p rad i / h p ii is ∼ . ∼ . ≃ . − .
5. Here the first (inner) hi represents a volume average and the second (outer)denotes a time average. Comparison with the previous simulation with comparable gasand radiation pressure (Krolik et al. 2007) is helpful. In that run, the variation of thestress is ∼
6, which is slightly bigger than that of the second stage of STD32; the nominaltime-averaged α -parameter is ≃ .
03, while it is a bit smaller, ∼ .
02, for STD32; the ratioof radiation pressure to gas pressure varies over the range ≃ . −
2, which is beyond that ofSTD32. Despite large fluctuations, STD32 clearly achieves a quasi-stable stage for the last150 orbits. All other four runs in this paper also show the feature of increasing h p rad i / h p i ,and they are terminated when a quasi-stable stage like the one in STD32 is reached. The saturation states of the other four runs are illustrated in Fig.2. Clearly thesaturation level of X64 and Y128 (left two panels in the graph) are similar while the Z512and DBLE (right two panels) have relatively higher mean values. The offset between theleft two and the right two begins right after the transient decay. It grows even larger afterthe first 60 orbits in both Z512 and DBLE: at that point, the radiation energy becomescomparable with the gas energy. 12 –To study convergence quantitatively, we need not only temporal and spatial averagesof the stress, but proper ways to normalize it. Here we choose the averaging time to befrom the end of the transient phase (10 orbits) to the end of the run; the volume-integralof the stress is used as the spatial average. Unlike the nonstratified case, the pressuresin our simulations show both consistent spatial gradients, depending on height from themidplane, and significant trends over time. We therefore employ three different methods ofnormalization, and examine the quality of numerical convergence in each case.We first normalize the different kinds of stress and energy by the initial volume-averagedtotal pressure P (sum of gas pressure and radiation pressure): P = 9 . × ergs cm − in our simulations. This normalization definition is extensively used in unstratified shearingbox simulations, but with the total pressure replaced by thermal pressure alone. The time-and box-averaged values are given in Table 1. Scanning across each line, one can see whichquantities are sensitive to resolution; in general, convergence has clearly been reached withregard to x and y cell size, but not with respect to ∆ z . The normalized Maxwell stress isconstant at ≃ .
03 when the resolution in the radial or azimuthal directions increases, butits value is almost doubled when the vertical resolution is raised by a factor of 2. Similarly,the magnetic energy and turbulent kinetic energy rise by about a factor of two when thevertical cell count is doubled, but are independent of the horizontal cell dimensions. Bycontrast, in unstratified simulations, when the resolution improves, the saturation leveleither decreases toward zero with a zero net-field configuration (Fromang & Papaloizou2007; Simon et al. 2009) or increases weakly for mean azimuthal field models (Guan et al.2009). We have net azimuthal field in our simulations, but it is not fixed, and even the signof the net azimuthal flux changes.A second useful normalization standard is the horizontal average of the time-dependenttotal pressure in the midplane. The time-averaged values for the stresses and energies 13 –normalized in this fashion are listed in Table 1, too. They depend on resolution in a wayvery similar to the ones using the absolute normalization except that their values are almostone order of magnitude smaller. Again, in this normalization there is no dependence onhorizontal resolution, but increasing resolution in z leads to larger values. For example, theMaxwell stress for DBLE is almost twice that of STD32.Considering that P is just an arbitrary initial guess for the total pressure, and P (0)does not reflect the properties of the whole box, it is more physical to normalize the energiesand stresses to the simultaneous volume-averaged total pressure, i.e., h P i ≡ h p rad + p i .The time evolution curves of the stress ratios using this normalization are plotted in Fig.3.Compared to the absolute stresses shown in Fig.1 and 2, the stress ratios after normalizationshow considerably smaller peak-to-peak variations. Normalized in this way, there is also asignificantly smaller increase when the resolution along the vertical direction increases. Thetime averaged values for this normalization are the third group in Table 1. The normalizedMaxwell stress is ∼ .
02 for STD32, X64 and Y128, and increases only to ≃ .
03 for Z512and DBLE. This sort of normalization is the best of the three to use for estimating theShakura-Sunyaev α parameter because it makes use of the actual volume-integrated totalpressure, and we see that with the best resolution employed here, we are approachingconvergence in defining its value. We emphasize, however, that one can speak of a singlevalue for this number only in terms of a particular location in the disk and after both avertical integration and a time average that encompasses many thermal times.To summarize this section, we find that numerical convergence with respect toresolution in the x and y , but not z , directions has been achieved for the absolute values ofstress and energy in a stratified shearing box. Increasing z resolution at the level we havereached leads to rising absolute values of stress and pressure. On the other hand, we comeclose to reaching convergence with respect to all three sorts of resolution for the ratio of 14 –stresses and energies to the time-dependent total pressure. In the next subsection, we showthat certain detailed features of the magnetic field show similar convergence properties andcast light on why stratified shearing boxes differ from unstratified. In order to investigate the effect of stratification on the magnetic field, we considerthe vertical correlation length of the field. As the azimuthal component dominates inour simulations, we restrict our attention to B y . We present two 3-D snapshots of thedistribution of azimuthal field strength in STD32 in Fig. 4. In the plot, color contours ofthe field strength are mapped onto the surfaces of the box. For contrast, we choose a pair oftimes, one ( t = 73 orbits) when the magnetic energy is low and one ( t = 250 orbits) whenit is high. We find the field distribution below ± H of height is distinct from that above:within that distance of the midplane, the field is turbulent, whereas at higher altitudes it ismuch more regular. Far from the midplane, long filamentary regions of relatively smooth B y with the same sign extend more than 1–2 H beyond the disk core. The field becomesdisordered again above ∼ ± H , where the sign of B y sometimes flips. The same featuresare consistently observed in all our other runs. Stratification clearly has a strong effect onthe qualitative organization of the magnetic field.This fact leads us to ask if stratification also influences the field structure of themidplane region. Let us first calculate the physical scale height by taking h c s / Ω i . The timeaveraged scale height is then denoted as H hereafter. This scale height is a better measureof the physical scale length of the disk than our guessed unit of length, H . The physicalscale height varies within the range ≃ . H for STD32, with time average 1 . H (similarvariations and average values are also found in X64 and Y128, see Table 1); its range ofvariation is slightly wider, ≃ . . H in Z512 and DBLE, and the average is also a bit 15 –larger, ∼ . H .We now calculate the vertical correlation length of B y in the disk core, which we defineas | z | ≤ H , i.e. within roughly half a physical scale height of the midplane. This choiceallows us to compare our results with those of Fromang & Papaloizou (2007), whose boxextended one physical scale height vertically. Note that our DBLE run has exactly the sameresolution as their STD64 (both are H /
64 vertically). The two-point correlation function of B y at t = t is defined as: C z ( B y ; t , l z ) ≡ L x L y Z Z dxdy R B y ( t , x, y, z ) B y ( t , x, y, z − l z ) dz R B y ( t , x, y, z ) dz , (1)where C z ( B y ; t , l z ) denotes the vertical correlation of B y at separation l z as a functionof time, and L x and L y are the box sizes along the x and y directions. FollowingFromang & Papaloizou (2007), we define an integrated correlation length as the integralover different separations: λ z, int ( B y ; t ) = Z C z ( B y ; t , l z ) dl z . (2)Note that in Fromang & Papaloizou (2007) the correlation length is calculated only in the x − z plane. We find that our value of the integrated correlation length would change byonly ∼
3% (see Table 1) if we had used their definition instead of ours. We also calculatethe correlation length defined as the full width at half maximum (FWHM) of C z , i.e., λ z, F W HM ( B y ; t ) = FWHM of C z ( B y ; t , l z ) . (3)We show the time histories for several definitions of λ z in STD32 and DBLE in the toppanel of Fig. 5. As this figure shows, both the integrated and the FWHM correlation lengthvary with time, but λ z, int has considerably larger fluctuations than λ z, F W HM . A close lookat the correlation function explains the discrepancy. In the middle of the same graph, threeinstantaneous correlation functions, at t = 73, 97 and 25 orbits in DBLE are presented. At 16 – t = 73 orbits, C z ( B y ) clearly shows positive correlation even at large separations, makingthe integrated correlation length λ z, int as large as 0 . H , about twice the λ z, F W HM found atthe same time. On the other hand, at t = 25 orbits, the correlation function turns negativewhen | l z | > . H , and thus a dip of λ z, int ( ∼ . H ) is observed at t = 25 orbits, when λ z, F W HM is ∼
60% large. The two measures of λ z become comparable to each other whenthe wings of C z ( B y ) asymptote to zero at larger separations. For example, at t = 97 orbits,when this occurs, λ z, F W HM is similar to, but slightly less than λ z, int . We can test whether λ z, int or λ z, F W HM is the more realistic measure of the correlation length by looking atsnapshots of B y , such as those shown under each correlation plot. There is no apparentchange in the lengthscales of magnetic field features between t = 25 and 97 orbits, despitethe factor of two change in the intergrated correlation length between those times. TheFWHM correlation lengths at these two times are nearly the same. Thus, we claim that theFWHM definition is a better measure of the correlation length than the integrated version:it corresponds more closely to one’s visual impression and varies less in time. The largefluctuations in λ z, int appear to be due to sensitivity to the tails of C z .As shown by the time averaged values listed in Table 1, λ z, F W HM decreases slightly inmagnitude when the resolution is doubled. We find h λ z, F W HM i t ∼ . H for STD32, butdrops to ∼ . H once the z − cell number is doubled. This change corresponds to a 30%decrease if we scaled the length with the physical scale height H : h λ z, F W HM i t ∼ . H for STD32 and ∼ . H for DBLE. As we mentioned before, the DBLE run possesses thesame resolution as STD64 in Fromang & Papaloizou (2007), in which the time averagedcorrelation length using their definition (an integrated one) is 0 . H ; with the samedefinition, the correlation length of DBLE is ∼ . H , approximately triple theirs.In sum, we find comparatively large scale ( & H ) structures in B y more than ≃ H from the midplane, and even in the region quite close to the midplane, features several 17 –times larger than found in unstratified simulations. We also find that the FWHM definitionof the vertical correlation length may be more useful than the integrated one. In terms of itsdegree of numerical convergence, λ z is similar to the ratio of stress to simultaneous pressure:at the resolution scales achieved in these simulations, it appears to be close to convergence,but not quite there, particularly with respect to resolution in the vertical direction. As discussed in the last subsection, large filament-like structures of magnetic fieldemerge above the core of the disk and extend vertically for another 1 − H . These featureshave also been observed in many previous simulations (e.g., Turner 2004; Hirose et al.2006). In Figure 6, we show a space-time diagram of the horizontally-averaged azimuthalcomponent of the field from a 50-orbit segment of STD32 . There are roughly ten episodesof field upwelling during this sample, corresponding to a period of ∼ B y alternates in successive events. The upward pattern speed is ∼ . H perorbit at the base of the plume ( ∼ ± H ), accelerating to ∼ − H per orbit near the topof the disk ( & ± H ). These events are almost symmetrical about the midplane except forsome small phase offsets and intensity variations.Computing the simple hydrodynamic Brunt-V¨ais¨al¨a frequency for our data wouldsuggest that the disk is always stable against buoyancy-driven instability below ∼ H , As the default dumping rate for STD32 is only one time per orbit, to obtain this data wecarried out a special high dump rate run by restarting STD32 from t = 100 orbits, runningfor 50 orbits, and dumping physical quantities every 0.1 orbits. These data therefore do notfollow STD32 exactly, but are physically and statistically equivalent. In that sense, we stillcall it STD32 here. 18 –and remains mostly stable at higher altitudes, including the upwelling regions. However,including magnetic forces leads to a very different result.There are two magnetic buoyancy modes to consider: the interchange mode andthe undulatory Parker instability. The former mode is not present in our simulations:evaluating its dispersion relation with horizontally averaged data shows that the disk isstable at nearly all times and in nearly all locations. To demonstrate that the undulatorymode does seem to be responsible for these repeated episodes of magnetic upwelling, webegin by writing down expressions for the square of the generalized magnetic Brunt-V¨ais¨al¨afrequency for this mode in two limits: fast and slow radiation diffusion (Tao & Blaes 2009,in preparation): N mag,fast = gv A c i + v A d ln | B | dz , (4) N mag,slow = N hyd + gv A c t d ln | B | dz . (5)Here N hyd is the hydrodynamic Brunt-V¨ais¨al¨a frequency for a thermally coupled gas andradiation mixture (Blaes et al. 2007), g = z Ω is the vertical component of gravity, v A is theAlfv´en speed, and | B | is the magnitude of the field strength. The isothermal sound speed inthe gas is c i = ( p/ρ ) / , and the total adiabatic sound speed is c t = [Γ ( p + E/ /ρ ] / , where p and E again are gas pressure and radiation energy density, and Γ is Chandrasekhar’sgeneralized adiabatic constant (Chandrasekhar 1967). The “slow diffusion” limit describesthe case in which the growth rate of the instability exceeds the photon diffusion rateso that photons are dynamically well-coupled to the fluid; in the “fast diffusion” limit,photons diffuse rapidly compared to the instability growth rate, so that there is no radiationpressure response to the mode and the perturbations in the gas are isothermal. Note thathydrostatic equilibrium with no magnetic tension forces is assumed in both expressions for N mag . Evaluating the growth rate using the wavelength of the fastest-growing mode in therapid diffusion limit (equation A15 of Blaes et al. (2007)), we find that, for regions above 19 – ∼ ± . H , the fast diffusion limit is the appropriate one, while at lower altitudes the slowdiffusion limit generally applies. In approximate terms, N mag can then be treated as theproper combination of the frequency squared under those two limits. The fastest-growingwavelength is always well confined in the box and well resolved numerically in most regionsexcept near the midplane.We plot the zero-frequency contour of the combined N mag (black curves) in thespace-time diagram of Figure 6. Instability takes place only outside the contour lines.Inside ± H , i.e. the core disk region, this mode is generally stable, although sometimesonly marginally so. Small unstable patches exist in the core region, and most of them areelongated in a way suggesting buoyancy, but they exist only briefly. The magnitude ofthe growth rate in this region is small (see panel b of Fig. 7), and the implied rise speedsare relatively slow. These small episodes of unstable buoyancy may help explain why wefind larger scale features in our stratified simulations than have been found in unstratifiedsimulations, but the fact that the wavelength of the fastest-growing modes exceeds the sizeof the box suggests that this explanation probably cannot completely answer the questionof the origin of our larger correlation length. On the other hand, N mag is negative almosteverywhere at altitudes above ∼ ± (2–3) H . At these higher altitudes, the undulatory Parkermode is almost always unstable with a linear growth rate ∼ (1–2)Ω.In order to show the operation of this mechanism in greater detail, we study a magnifiedview of a small domain of the spacetime diagram. Data from the region within the whiterectangle in Figure 6 are displayed on a larger scale in Figure 7. As discussed, the risingpattern is roughly symmetric about the midplane, so by choosing a sample from the upperhalf of the disk we should not lose any interesting physics. The time range of this sampleis 10 orbits and contains almost two complete buoyancy episodes. In this figure, we plotnot only B y and N mag , but also three characteristic speeds: the gas vertical velocity v z , 20 –the magnetosonic speed v ms , and the Poynting flux speed v Pf . The last one is defined as v Pf ≡ πS z /B , where S z represents the z component of the Poynting flux. Thus, positive v Pf means field energy is transported upward and corresponds to field ascending. ThePoynting flux speed can also be written as v Pf = v z − ( v · B ) B z / | B | ; this form shows how,even in the MHD limit, the gas and magnetic velocities can differ to the extent that thefluid can slide along field lines.Panels a, c, and e of Figure 7 present an enlightening contrast: although the magneticfield spacetime diagram consistently shows features in B y rising upward from near themidplane all the way to the top of the box, the v Pf and v z illustrations demonstrate thatconsistent upward field and fluid motion begins only at ≃ H above the midplane. In fact,the bottom edge of the true upwelling region corresponds closely to the lower boundary ofthe Parker instability region. Thus we see that the apparent motion of regions of strong fieldin the disk core and lower corona is a wave or pattern motion almost completely divorcedfrom genuine mass or field transport. In fact, the pattern speeds observed in the simulationare close to the Alfv´en speed, suggesting that the waves in question are either Alfv´enor slow magnetosonic waves. The latter mode appears to be particularly important forshort-wavelength, high-frequency variability, for which a distinct anti-correlation betweenmagnetic and gas pressure fluctuations can be seen.As already remarked, genuine rising motion of magnetic field and gas is closelyassociated with the onset of magnetic buoyancy instability. Further evidence thatundulatory Parker mode instability is the origin of upwelling is provided by the fact that theregions of greatest growth rate (the dark blue strips in the N mag panel) coincide with theleading edges of regions where both v Pf and v z are large. In other words, the accelerationto high upward velocity takes place where the instability grows most strongly.Although true upward motion certainly takes place at altitudes several scale heights 21 –from the midplane, it is in an absolute sense relatively slow. Typical speeds for both v Pf and v z are ∼ (0 . . H Ω. Both are subsonic, only ∼ O (0 . v ms (note that the units inthe v ms panel are 10 cm s − rather than the 10 cm s − of the v Pf and v z panels). Theslowness of these speeds demonstrates that even though much of the corona is buoyantlyunstable, the net accelerations are never more than a fraction of gravity.The relationship between v Pf and v z changes sharply between the disk core and thecorona, as can be seen in panel f of Figure 7, which shows their ratio. In the disk core, theyare entirely uncorrelated, consistent with our conclusion that there are no genuine flows ofeither in that part of the disk. On the other hand, in the coronal region, | z | > H , althoughnot strictly proportional, more often than not they vary together, with v Pf ≃ (1–2) v z .Where the undulatory Parker mode grows, field lines develop upward bends and can carrygas upward with them. However, the curvature of the field lines permits the gas to slipdiagonally downward relative to the field, an effect explaining why the mean upward speedfor gas is smaller than that for field.To sum up, we find that the undulatory Parker instability is the major driver of themagnetic upwelling that has been so often noted in stratified shearing box simulations.Most of the time, it is marginally stable in the core region ( | z | . H ), and unstable above it.Consistent upward mass and field motions begin only at the lower boundary of the Parkerinstability, and both of them are one order of magnitude slower than the magnetosonicspeed. The upward velocity of gas motion is generally half of the field flux velocity as thegas may slide towards the field line valleys. 22 – The turbulence produced by the MRI is frequently described in the literature as“horizontal” because the initial work on shearing boxes found that the kinetic energy ofturbulent motions in the x − and y − directions, corresponding to the radial and azimuthaldirections if the box were placed in an actual disk, was substantially greater than that inthe z − , or vertical, direction. For example, Hawley et al. (1995) found that in their fiducialsimulation of a three-dimensional, but unstratified, shearing box (in this case, with initialmagnetic field uniform and vertical), h ρv x i ≃ h ρ ( δv y ) i ≃ h ρv z i , where δv y is the y-velocityafter subtracting out the shearing motion. An initially azimuthal field also led to primarilyhorizontal motions, although with somewhat larger amplitude in the radial than in theazimuthal direction. Similarly, in a stratified shearing box that extended ± √ H fromthe midplane, Stone et al. (1996) found that the turbulence was predominantly azimuthal: h ρ ( δv y ) i ≃ h ρv x i ≃ h ρv z i , with little dependence on whether the equation of state wasisothermal or adiabatic or on whether the initial field was vertical or azimuthal.However, the situation changes when stratified shearing boxes with greater verticalextent are simulated, and the periodic boundary conditions used by Stone et al. (1996)are replaced by outflow boundary conditions. Miller & Stone (2000) reported that in theirsimulations, in which the box extended to ± √ H from the midplane, the amplitude ofvertical motion near the midplane ( | z | ≤ √ H ) was almost as great as the amplitudeof the azimuthal turbulent motions: h ρ ( δv y ) i ≃ . h ρv z i , nearly independent of whetherthe initial field was vertical or azimuthal. Their result for the ratio h ρv z i / h ρv x i was moreambiguous: its value was ≃ . ≃ . ± H of the midplane, we have 23 – h ρ ( δv y ) i ≃ . h ρv x i ≃ . h ρv z i . In the highest-resolution simulation, the azimuthalmotions grow slightly relative to both the radial and vertical motions, but overall, all threevelocity components remain comparable in magnitude: h ρ ( δv y ) i ≃ . h ρv x i ≃ . h ρv z i .That such a qualitative contrast in the geometry of turbulent motions is associated with thechange from unstratified to stratified boxes strongly suggests that its physical origin can beidentified with the principal physical contrast between the two cases: vertical gravity, whichin this context drives magnetic buoyancy.
4. CONCLUSIONS
We have shown that in a stratified shearing box, unlike an unstratified one, themagnetic saturation level does not diminish with increasing grid resolution. Instead, atthe resolutions used in our simulations, the Maxwell stress is independent of the x and y resolution. When the vertical resolution increases, both the stress and the pressure grow.Despite the nearly linear increase of both magnetic stress and total pressure with resolution( h B r B φ i , h P i ∝ (∆ z ) − ), the ratio of stress to pressure is almost converged at this resolution(as also found in the isothermal simulations of Davis et al. (2009)), increasing only slowlywith finer vertical resolution.The field structure of a stratified box is quite different from that of an unstratified box.At high altitude ( | z | > H ), the field is dominated by its azimuthal component and formssmoothly distributed filamentary structures with vertical thicknesses & H . In the core ofthe disk, the region unstratified boxes are thought to mimic, the field is turbulent, but withsignificantly larger scale features than found in unstratified simulations: for simulationswith identical resolution, the vertical correlation length of B y is ∼ ≃
30% when 24 –the resolution is doubled.There are several ways to understand the contrast between magnetic field behavior inthe stratified and unstratified cases. One is through dimensional analysis. Without verticalgravity, the ratio c s / Ω, although well-defined, loses its physical significance. Consequently,the only remaining significant lengthscale in the problem is the gridscale. As a result, thetwo lengthscales determined by the turbulence, the field correlation length and v A / Ω, bothtrack the gridscale and approach zero as the resolution grows finer. On the other hand,with vertical gravity, c s / Ω plays an explicit role in the turbulent dynamics, and both thefield correlation length and v A / Ω can become associated with it. A related argument hasbeen made by Vishniac (2009).A second approach ties the result more directly to gravitational dynamics. As wehave shown, buoyancy plays an important role in the character of MRI-driven MHDturbulence in stratified shearing boxes. Outside the midplane region ( | z | ≥ ± (2–3) H ), N mag is consistently negative, a signal that buoyant regions are continually being created.Consequently, at any given time there is nearly always a rising magnetic filament in thatregion. Even inside the midplane region, vertical motions, presumably excited by thebuoyantly-driven pressure fluctuations on its edges, are greater than they are in unstratifiedsituations. In notable contrast to studies of both unstratified shearing boxes and stratifiedshearing boxes with periodic vertical boundary conditions and containing only a few verticalscale heights, we find that the turbulent motions are fully three-dimensional, with thekinetic energy in vertical motion comparable to that in either of the horizontal directions.Because fully three-dimensional motion is essential to dynamo action (e.g., as reviewedby Cowling (1981) and first emphasized in this context by Brandenburg et al. (1995)), thebuoyant excitation of vertical motion is essential to maintaining the magnetic field: as fluidelements rise and fall, they create vertical field from horizontal. The pressure scale height 25 –sets the characteristic scale for vertical motions, linking the magnitude of the field (withproportionality constant ∼ √ πρ Ω) to this characteristic lengthscale (see Vishniac (2009)for a proposed scaling relation). At saturation, field amplification by dynamo action isbalanced by two field loss mechanisms: dissipation and expulsion of field from the box. Aspreviously shown by Hirose et al. (2006), the former strongly dominates the latter whenallowance is made for photon energy losses.We have also used these simulations to identify the dynamical mechanism responsiblefor the upward magnetic motions commonly seen in previous stratified simulations: it is theundulatory Parker instability. This mode is marginally stable within the core region, butis almost always unstable outside that region. True upward motions of both magnetic fieldand fluid begin at the instability boundary, but are one order of magnitude slower thanthe magnetosonic speed. Because the gas can slide along field lines and fall into field linevalleys, its mean upward velocity is generally only about half the Poynting velocity.Although we have allayed fears that the true converged state of MRI-driven MHDturbulence is zero stress, the question of numerical convergence of the stratified case remainsopen. We understand neither why the stress and pressure are strongly dependent on vertical(but not horizontal) cell size nor at what resolution this dependence may weaken. Theorigin of the large fluctuations as a function of time observed in box-integrated quantitiesin these simulations is also unclear, although there are hints that their magnitude may haveto do with the radial width of the box (e.g., Fromang & Stone 2009).Several physical questions also remain unanswered: The pressure scale height certainlysets a physical lengthscale, but can we understand more specifically how it determines thevertical correlation length in the disk core’s MHD turbulence? It seems plausible that thefield buoyancy leads to larger scale features in the core region, as the field senses the verticaloutflow boundaries via the magnetic upwelling, but exactly how is that communicated 26 –to the midplane? Lastly, one of the most striking features of the disk corona is thequasi-regular alternation in sign of B y : what causes this alternation and what determinesits characteristic lengthscale and timescale?We thank Omer Blaes and Shane Davis for very useful discussions. The authorsacknowledge the Texas Advanced Computing Center (TACC) at The University of Texas atAustin for providing HPC resources that have contributed to the research results reportedwithin this paper. Part of the work was also performed on the Johns Hopkins HomewoodHigh-Performance Computing Center cluster. This work was partially supported by NASAATP Grant NNG06GI68G and by NSF Grant AST-0507455. 27 – REFERENCES
Balbus, S.A. & Hawley, J.F. 1991, 376, 214Balbus, S.A., & Hawley, J.F. 1998, Rev. Mod. Phys., 70, 1Blaes, O., Hirose, S. & Krolik, J.H. 2007, ApJ,664, 1057Brandenburg,A. et al. 1995, ApJ, 446, 741Chandrasekhar, S. 1967, An Introduction to the Study of Stellar Structure (New York:Dover)Fromang, S. & Papaloizou, J. 2007, A&A, 476, 1113Fromang, S. et al. 2007, A&A, 476,1123Fromang, S. & Stone, J.M.,astro-ph 09064422Guan, X. et al. 2009, ApJ, 694, 1010Hawley, J.F., Gammie, C.F., & Balbus, S.A. 1995, ApJ, 440, 742Hirose, S., Krolik, J.H., & Stone, J.M. 2006, ApJ, 640, 901Hirose, S., Krolik, J.H., & Blaes, O. 2009, ApJ, 691, 16Krolik, J.H., Hirose, S., & Blaes, O. 2007, ApJ, 664, 1045Lesur, G. & Longaretti, P.-Y. 2007, MNRAS, 378, 1471Miller, K.A. & Stone, J.M. 2000, ApJ, 534, 398Shakura, N.I. & Sunyaev, R.A. 1973, A&A, 24, 337Simon, J.B., Hawley, J.F. & Beckwith, K. 2009, ApJ, 690, 974 28 –Stone, J.M., Hawley, J.F., Gammie, C.F., & Balbus, S.A. 1996, ApJ, 463, 656Tao, T. & Blaes, O. 2009, in preparationTurner, N.J. 2004, ApJ, 605, 45Vishniac, E.T. 2009, ApJ, 696, 1021Cowling, T.G. 1981, ARAA, 19, 115Davis, S.W., Stone, J.M. & Pessah, M. 2009, arXiv 0909.1570This manuscript was prepared with the AAS L A TEX macros v5.2. 29 – S t r ess × Ω / ( e r g c m − o r b i t − ) E n e r g y ( e r g c m − o r b i t − ) Fig. 1.—
Volume-integrated stresses (upper panel) and energies (lower panel) in STD32 as afunction of time.
Top panel : Maxwell stress (red), Reynolds stress (green) and the sum of thesetwo (black).
Bottom panel : internal energy (black), radiation energy (blue), 10 times the magneticenergy (red), 10 times the kinetic energy(green) and total energy (grey dash-dotted). The quantitiesare all vertically integrated and the stress is multipled by 3Ω /
30 – S t r ess × Ω / ( e r g c m − o r b i t − ) S t r ess × Ω / ( e r g c m − o r b i t − ) Fig. 2.—
Volume-integrated stresses in Y128(upper left panel), X64(lower left panel), Z512(upperright panel) and DBLE(lower right panel) as a function of time. The color curves are coded as inthe top panel of Fig.1.
31 – S t r ess / < P > S t r ess / < P > Fig. 3.—
Volume-integrated stresses in ratio to simultaneous volume-averaged total pressure forSTD32, Y128, X64 (left, from top to bottom), Z512 and DBLE (right, from top to bottom). Colorcode is the same as in Fig.2. −4 • −2 • • • B y (G)t= 73 orbits −1 0 1−4−2024 −1 0 1−505−4−2024 −505 −4 • −2 • • • B y (G)t= 250 orbits −1 0 1−4−2024 −1 0 1−505−4−2024 −505 Fig. 4.—
Snapshots of the B y distribution on the box surfaces at t = 73 orbits and t = 250 orbitsin STD32.
32 – λ z ( B y ) ( H ) −1.0 −0.5 0.0 0.5 1.0l z (H)−0.20.00.20.40.60.81.0 C z ( B y ;t = ) −1.0 −0.5 0.0 0.5 1.0l z (H)−0.20.00.20.40.60.81.0 C z ( B y ;t = ) −1.0 −0.5 0.0 0.5 1.0l z (H)−0.20.00.20.40.60.81.0 C z ( B y ;t = ) −1.0 −0.5 0.0 0.5 1.0x (H)−1.0−0.50.00.51.0 z ( H ) −1.0 −0.5 0.0 0.5 1.0x (H)−1.0−0.50.00.51.0 z ( H ) −1.0 −0.5 0.0 0.5 1.0x (H)−1.0−0.50.00.51.0 z ( H ) −2 • −1 • • • B y (G)t=73 orbits t=97 orbits t=25 orbits Fig. 5.—
Top : Radially averaged vertical correlation length of B y for STD32 (black solid for λ z, F W HM and blue dashed for λ z, int ) and DBLE (red solid for λ z, F W HM and green dash-dottedfor λ z, int ) as a function of time; Middle : The correlation function at t = 73, 97 and 25 orbits ofDBLE; Bottom : B y contours from DBLE in y = 0 plane at t = 73, 97 and 25 orbits.
33 –Fig. 6.—
Spacetime diagram of the azimuthal component of the magnetic field. The contour lineson top delineate where and when the magnetic Brunt-V¨ais¨al¨a frequency squared (black) changesfrom positive (inside the contours) to negative (outside the contours) and the plasma β (white)increases from less than unity (outside the contours) to greater (inside the contours). The whiterectangle shows the section magnified in Fig.7.
34 –Fig. 7.—
Magnified view of a portion of the spacetime diagram for STD32: (a) azimuthal magneticfield strength; (b) magnetic Brunt-V¨ais¨al¨a frequency squared; (c) vertical velocity of the Poyntingflux; (d) magnetosonic speed; (e) vertical gas velocity; (f) ratio of the Poynting flux velocity to thevertical gas velocity. See legend on the top of each panel for units. The black and white contourlines are the same as in Fig.6.
35 –Table 1. Properties of the simulations in the paper.
STD32 X64 Y128 Z512 DBLEresolution 32 × ×
256 64 × ×
256 32 × ×
256 32 × ×
512 64 × × max / Σ (%) 2.5 2.0 2.0 1.1 1.4 H = hh c S / Ω ii ( H ) 1 .
65 1 .
67 1 .
67 2 .
00 1 . hh t cool ii (orbits) 8 . . . . . h λ z, int ( B y ) i t ( H ) 0 .
31 0 .
35 0 .
34 0 .
35 0 . h λ z, int ( B y ; y = 0) i t ( H ) 0 .
31 0 .
36 0 .
33 0 .
35 0 . h λ z, F W HM ( B y ) i t ( H ) 0 .
26 0 .
25 0 .
25 0 .
23 0 . hh− B x B y / π ii /P .
029 0 .
030 0 .
031 0 .
064 0 . hh ρv x δv y ii /P . × − . × − . × − .
015 0 . hh B / π ii /P .
11 0 .
11 0 .
11 0 .
24 0 . hh B x / π ii /P . × − . × − .
010 0 .
025 0 . hh B y / π ii /P .
095 0 .
090 0 .
096 0 .
21 0 . hh B z / π ii /P . × − . × − . × − .
014 0 . hh ρδv / ii /P .
029 0 .
031 0 .
031 0 .
061 0 . hh ρv x / ii /P .
012 0 .
012 0 .
012 0 .
024 0 . hh ρδv y / ii /P .
010 0 .
011 0 .
011 0 .
023 0 . hh ρv z / ii /P . × − . × − . × − .
014 0 . hh− B x B y / π i /P (0) i . × − . × − . × − . × − . × − hh ρv x δv y i /P (0) i . × − . × − . × − . × − . × − hh B / π i /P (0) i .
014 0 .
014 0 .
015 0 .
040 0 . hh B x / π i /P (0) i . × − . × − . × − . × − . × − hh B y / π i /P (0) i .
012 0 .
012 0 .
013 0 .
036 0 . hh B z / π i /P (0) i . × − . × − . × − . × − . × − hh ρδv / i /P (0) i . × − . × − . × − . × − . × −
36 –Table 1—Continued
STD32 X64 Y128 Z512 DBLE hh ρv x / i /P (0) i . × − . × − . × − . × − . × − hh ρδv y / i /P (0) i . × − . × − . × − . × − . × − hh ρv z / i /P (0) i . × − . × − . × − . × − . × − hh− B x B y / π i / h P ii .
020 0 .
021 0 .
021 0 .
029 0 . hh ρv x δv y i / h P ii . × − . × − . × − . × − . × − hh B / π i / h P ii .
075 0 .
073 0 .
076 0 .
13 0 . hh B x / π i / h P ii . × − . × − . × − .
011 0 . hh B y / π i / h P ii .
064 0 .
062 0 .
065 0 .
11 0 . hh B z / π i / h P ii . × − . × − . × − . × − . × − hh ρδv / i / h P ii .
020 0 .
021 0 .
021 0 .
028 0 . hh ρv x / i / h P ii . × − . × − . × − .
011 0 . hh ρδv y / i / h P ii . × − . × − . × − .
010 0 . hh ρv z / i / h P ii . × − . × − . × − . × − . × − hh− B x B y / π i / h B / π ii .
27 0 .
28 0 .
28 0 .
21 0 . hh− B x B y / π i / h ρv x δv y ii . . . . . hh− B x B y / π ii / hh B / π ii .
26 0 .
27 0 .
28 0 .
27 0 . hh− B x B y / π ii / hh ρv x δv y ii . . . . ..