{'fraction_non_alphanumeric': 0.031125038449707782, 'fraction_numerical': 0.020051522608428175, 'mean_word_length': 4.381440099317194, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Turbulence and Steady Flows in 3D Global Stratified MHD Simulations of Accretion Disks

We present full 2 Pi global 3-D stratified MHD simulations of accretion disks. We interpret our results in the context of proto-planetary disks. We investigate the turbulence driven by the magneto-rotational instability (MRI) using the PLUTO Godunov code in spherical coordinates with the accurate and robust HLLD Riemann solver. We follow the turbulence for more than 1500 orbits at the innermost radius of the domain to measure the overall strength of turbulent motions and the detailed accretion flow pattern. We find that regions within two scale heights of the midplane have a turbulent Mach number of about 0.1 and a magnetic pressure two to three orders of magnitude less than the gas pressure, while outside three scale heights the magnetic pressure equals or exceeds the gas pressure and the turbulence is transonic, leading to large density fluctuations. The strongest large-scale density disturbances are spiral density waves, and the strongest of these waves has m=5. No clear meridional circulation appears in the calculations because fluctuating radial pressure gradients lead to changes in the orbital frequency, comparable in importance to the stress gradients that drive the meridional flows in viscous models. The net mass flow rate is well-reproduced by a viscous model using the mean stress distribution taken from the MHD calculation. The strength of the mean turbulent magnetic field is inversely proportional to the radius, so the fields are approximately force-free on the largest scales. Consequently the accretion stress falls off as the inverse square of the radius.


Introduction
The magneto-rotational instability (MRI) is a powerful process to drive turbulence and angular momentum transport in protoplanetary disks, ultimately enabling the accretion of matter onto the central object Balbus & Hawley 1998).
There is vast literature studying this mechanism in local shearing box simulations with an ideal MHD description (Brandenburg et al. 1995;Hawley et al. 1995Hawley et al. , 1996Matsumoto & Tajima 1995;Stone et al. 1996;Sano et al. 2004).
The effect of non-ideal MHD on MRI, regarding the issue of resistive protoplanetary disks, was mostly studied in local box simulations (Blaes & Balbus 1994;Jin 1996;Sano et al. 2000;Fleming et al. 2000;Sano & Stone 2002a,b;Fleming & Stone 2003;Inutsuka & Sano 2005;Turner et al. 2007;Turner & Sano 2008;Turner et al. 2010). The various studies showed that at a certain level of resistivity the MRI will be suppressed.
Up to now there is no prescription of resistive profile in protoplanetary disks which applies for longer timescales. It is known that the dust grains control the ionization level in protoplanetary disk. The particle cross section and dust-to-gas ratio are most important parameters for defining the ionization level of the gas (Turner et al. 2006;Wardle 2007). Most studies of non-ideal MRI turbulence use a static dust distribution and neglect dust growth and evolution (Simon & Hawley 2009;Turner et al. 2006Turner et al. , 2007Dzyurkevich et al. 2010). But exactly in non-turbulent regions, the dust particles can grow, reducing quickly the cross section and so driving to better ionization levels enabling again MRI . In our study we focus on ideal MHD, which applies to sufficiently ionized disk regions depleted of small dust grains. This also applies to the innermost hot parts of protoplanetary accretion disk and even extended radial regions as expected for transitional disks (Chiang & Murray-Clay 2007). Here, an MRI turbulent ionization front starting at the inner rim of the disk propagates radially outward and the disk get evacuated from inside-out.
Recent results on the MRI obtained in local unstratified box simulations by Lesur & Longaretti (2007), Fromang et al. (2007), Simon & Hawley (2009) and Fromang (2010) show that occurrence and saturation level of MRI in zero net flux simulations is controlled by the magnetic Prandtl number 1 . Ideal MHD simulations including stratification found at least convergence for zero-net flux local simulations (Davis et al. 2010;Flaig et al. 2010). Here, the vertical resolution plays the key role for the convergence (Shi et al. 2010). Despite the simplicity of box simulations and their interesting results obtained over the last years, the study of radial extended structures are very restricted. For instance, the aspect ratio of the boxes is known to influence the saturation level of turbulence (Bodo et al. 2008). Besides local box simulations, global simulations of MRI have also been performed (Armitage 1998;Hawley & Krolik 2001;Steinacker & Henning 2001;Hawley 2001;Arlt & Rüdiger 2001;Fromang & Nelson 2006;Lyra et al. 2008;Fromang & Nelson 2009).
They confirmed the picture of a viscously spreading disk as a proxy for the action of MHD turbulence. Recently, the first global non-ideal MHD simulation , which included a radial dead-zone / active zone interface demonstrated the importance of the inner edge of the dead zone as a trap for planetesimals and even small planets.
So far only finite difference schemes as implemented in the ZEUS and Pencil codes have been used to perform global simulations. However, a Godunov code would have several advantages over a finite difference scheme. Without using any artificial viscosity the code solves the MHD Riemann problem and can better handle the supersonic MHD turbulence in corona regions of the disk. Several papers recognized the importance of Godunov-type shock capturing upwind schemes for future astrophysical simulations Mignone et al. 2007;Flock et al. 2010). In this project we use the Godunov code PLUTO (Mignone 2009) to 1 The ratio of viscosity vs. resistivity perform isothermal global MHD simulations of protoplanetary disks. In future work, we will switch on the total energy conservation property of this scheme and include radiation MHD to follow the temperature evolution in global disk.
An important issue in stratified global simulations is the relatively low resolution per scale height compared to what is possible in local box simulations. In addition, the extent of the azimuthal domain is often restricted to save computational time. The first 2π global disk simulations were performed by Armitage (1998), Hawley (2000) and Arlt & Rüdiger (2001) for a short period of time. But most of the global simulations were performed in restricted azimuthal domains like π/4 or π/2 (Fromang & Nelson 2006Dzyurkevich et al. 2010). By restricting the azimuthal domain one also restricts the largest possible mode in the domain. One of the goal of this paper is to investigate whether these modes affects the nonlinear state of the turbulence.
The standard viscous α−disk theory (Shakura & Sunyaev 1973) introduces an effective turbulent viscosity η = ρν = αP/Ω with the local thermal pressure P and the orbital frequency Ω, arising from undefined magnetic or hydrodynamic processes, transporting angular momentum outward and allowing mass accretion onto the star. Lynden-Bell & Pringle (1974) calculated the radial mass accretion rate and the radial accretion velocity for a 1D viscous disk model as a local function of surface density Σ and the value for α. Interestingly 2D viscous disk models (Kley & Lin 1992) showed the appearance of meridional outflows. Here, the mass flows radially outward near the midplane compensated by increased radial inflow at upper layers of the disk to allow for net-accretion, for α < 0.05. Much emphasis was given to this radial outflow and its role for the radial transport for grains and chemical species over large distances and relative short time scales (Keller & Gail 2004;Ciesla 2009).
In addition, we will investigate the onset of a vertical outflow as it was described in local box simulations (Suzuki & Inutsuka 2009;Suzuki et al. 2010) using a net flux vertical field.
Such outflows can be launched in the magnetized corona region of the MRI turbulent disk (Miller & Stone 2000;Machida et al. 2000). They could have an important effect on the dissipation timescales of accretion disks and may be related to jet production from accretion disks (Ferreira et al. 2006). An interesting property of MRI in stratified zero-net flux simulations is the emergence of a "butterfly" pattern, an oscillating mean azimuthal magnetic field with an period of 10 local orbits. It was found in many local MRI simulations, recently again by Davis et al. (2010), Gressel (2010), Flaig et al. (2010) and in global simulations by Sorathia et al. (2010) and Dzyurkevich et al. (2010). We indeed identify such a "butterfly" pattern in our global runs, which was suggested to be linked to magnetic dynamo action in accretion disks (Sorathia et al. 2010;Gressel 2010).
Our paper is organized in the following way. In Section 2 we will present our model setup and the numerical configuration. Section 3 contains the results Section 4 and 5 will provide a discussion, summary and outlook for our work.

Model setup
The setup follows closely the disk model which is presented by Fromang & Nelson (2006. We define the cylindrical radius with R = r sin (θ) with the spherical radius r and polar angle θ. The initial density, pressure and azimuthal velocity are set to be in hydrostatic equilibrium.
We choose an isothermal equation of state. The pressure is set to P = c 2 s ρ with c s = c 0 · 1/ √ R.
The azimuthal velocity follows For the initial velocities V R and V θ we use a white noise perturbation amplitude of V Init R,θ = 10 −4 c s . We start the simulation with a pure toroidal magnetic seed field with constant plasma beta β = 2P/B 2 = 25. The radial domain extends from 1 to 10 radial code units (CU) 2 with radial buffer zones from 1 to 2 CU and 9 to 10 CU. In the buffer zones we use a linearly increasing resistivity.
This damps the magnetic field fluctuations and suppresses boundary interactions, especially for the closed boundary runs. Our buffer zone follows mainly the ones used in global simulations by (Fromang & Nelson 2006Dzyurkevich et al. 2010). The θ domain is set to θ = π/2 ± 0.3, corresponding to ±4.3 scale heights. We calculated in total five disk models. Three models cover the complete 2π azimuthal domain (FC, FO and BO in Table 1) and two models are constrained to π/4 (PC, PO in Table 1). The simulation FO is also used for a test model FOR which is described later. The simulation BO has the best resolution. One subset of models has a closed boundary (FC and PC). Here we use a reflective radial boundary with a sign flip for the tangential magnetic fields and a periodic boundary condition for the θ direction. A second subset of models has an outflow boundary condition (FO,PO and BO). Here, we use a relaxation function in the radial buffer zones which reestablishes gently the initial value of density over a time period of one local orbit. In the buffer zones we set: ρ new = ρ − (ρ − ρ Init ) · ∆t/T Orbits . Our outflow boundary condition projects the radial gradients in density, pressure and azimuthal velocity into the radial boundary and the vertical gradients in density and pressure at the θ boundary. We ensure to have no inflow velocities. For an inward pointing velocity we mirror the values in the ghost cell to ensure no inward mass flux. The θ boundary condition for the magnetic field are also set up to be zero gradient, which approximates "force-free" -outflow conditions. We also ensure the force free 2 We refer to CU instead of a physical length unit because ideal MHD simulations without radiation transport are scale free. Thus our simulations could represent a disk from 1 to 10 AU as much as a disk from 0.1 to 1 AU. Only explicit dust physics and radiative transfer will introduce a realistic physical scale. character of the tangential components for the radial boundary by adjusting the 1/r profile in the magnetic field components in the ghost cells. The normal component of the magnetic field in the ghost cells is always set to have ∇ · B = 0. We set the CFL value to 0.33. Also higher CFL values were successfully tested and will be used for future calculations. We use a uniform grid with an aspect ratio of the individual cells at 5 CU of 1 : 0.67 : 1.74 (∆r : r∆θ : r∆φ). Using a uniform grid instead of a logarithmic grid, where ∆r/r is constant, has the disadvantage that it will reduce the accuracy in the sense that the inner part of the disk is poorly resolved, compared to the outer part of the disk: H(1AU)/∆r < H(10AU)/∆r.
However for the uniform grid, the relative broad radial inner buffer zone lies in the poorly resolved disk part and is excluded from analysis. The outer parts of the disk are, compared to a logarithmic grid with the same resolution, better resolved. Logarithmic grid requires a much smaller buffer zone, e.g. a logarithmic grid would place one third of the total number of grid cells in the first ninth of the domain, between 1 and 2 AU. Of course, using a uniform grid will always restrict the range of the radial domain and for more radially extended simulation the need of a logarithmic grid is mandatory.
For all runs we employ the second order scheme in PLUTO with the HLLD Riemann solver (Miyoshi & Kusano 2005), piece-wise linear reconstruction and 2 nd order Runge Kutta time integration. We treat the induction equation with the "Constrained Transport" (CT) method in combination with the upwind CT method described in . The detailed numerical configuration is presented in Flock et al. (2010). Our high resolution run BO was performed on a Blue-gene/P cluster with 4096 cores and was calculated for over 1.5 million time steps which corresponds to 1.8 million CPU hours.

Code Units vs. Physical Units
Isothermal ideal MHD simulations are scale-invariant. One has to define unit-variables to transform from code to cgs units. We can set three independent values to define our problem. This is gas density for which we choose for instance ρ u = 10 −10 g/cm 3 , the radial distance unity as length 1CU = 1AU and the Keplerian velocity v u = √ G · M ⊙ /l u with the gravitational constant G and the solar mass M ⊙ .
With those three quantities, we translate the values for our measured surface density and mass accretion rate into cgs units. Using this values we derive at 1 AU at the midplane a gas density of ρ = 10 −10 g/cm 3 with a Keplerian velocity of v K = 2.98 * 10 6 cm/s. With this the surface density becomes 524g/cm 2 at 1AU. Gas velocities and the Alfvénic speed are always presented in units of the sound speed for convenience.

Turbulent stresses
The α parameter relates the turbulent stresses to the local thermal pressure. For the calculation of the α values we measure the Reynolds and Maxwell stresses, which are the R − φ components of the respective stress tensors. The Reynolds stress is calculated as T R = ρv ′ φ v ′ R and the Maxwell stress as T M = B ′ φ B ′ R /4π with the turbulent velocity or magnetic fields, e.g.
The mean component of the velocity and magnetic field are always calculated only along the azimuthal direction because of the radial and vertical gradients in the disk. In our simulations, the amplitude of Maxwell stress is about three times the Reynolds stress. For the total α value we integrate the mass weighted stresses over the total domain The respective turbulence enhanced viscosity can now be represented as ν = αHc s with the height of the disk H and the sound speed c s .

Disk evolution
We first describe the typical evolution for an azimuthal MRI (  models. During the time period between 800 and 1200 inner orbits, we get a relatively constant α value of 5 · 10 −3 (model BO). We mark three different time stages of the turbulent disk evolution: In period I (0 to 800 inner orbits), the turbulence is not yet saturated. After a strong initial rise due to the net azimuthal field the turbulent decays to a level where self-sustained turbulence is possible, e.g., the loss of magnetic flux in the vertical direction is balanced by the generation of magnetic flux in the turbulent flow. The nature of this generation of magnetic fields can be an indication for dynamo action, e.g., an αΩ dynamo (Brandenburg & Subramanian 2005;Gressel 2010), but detailed studies of this effect will be subject to future work. In period II, we have a quasi steady state of at least 400 inner orbits. During this time period, we do all our analysis. In the period III, a comparison of the models becomes less useful. The models with lower resolution and open boundary (PO and FO) show a decreasing α-stress in time. Thus they are not useful for long time integrations past 1200 inner orbits. In the closed models, on the other hand, the magnetic flux cannot escape vertically (PC and FC) and therefore, turbulence does not decay.
On the contrary, turbulence even increases in these runs as the flux in the box cannot efficiently escape.
A closer view shows that the stress can oscillate locally on shorter time scales. In Fig. 1, bottom right, we plot the mass weighted stresses αρ at 4 AU over height and local orbits. The strength of stresses locally oscillates with a period of around 5 local orbits. The maxima in the stresses always appear first in the midplane and then propagate vertically. These oscillations in the stresses are connected to the B φ "butterfly" structures where the azimuthal mean field oscillates with a frequency of 10 local orbits (Fig. 13). Every change of sign in the mean B φ is now correlated with a minimum in the stresses, which both occur every 5 local orbits. At the same time this plot shows the importance of the stresses in a region up to 3 disk scale heights.
Although one can always define a total α-parameter in the disk, the spatial variations are enormous, especially in the vertical direction. The vertical α profile is plotted in Fig.1, bottom left. For model BO, the turbulent stress at the midplane increases from 2.0 · 10 −3 up to 8 · 10 −2 at 4 scale heights. The simulations with moderate resolution show significant lower values around the midplane due to the lack of resolution there. For the closed models, the stresses in the corona are artificially increased due to the periodic boundary.

Radial profile of turbulent stress
Beside the vertical profile of turbulent stress, which has been already studied in local box simulations, the radial profile of the turbulent stress can only be addressed in global simulations.
In Fig. 1, top right, we present the radial α profile, averaged between 300 and 1200 inner orbits.
In the inner buffer zones (1 -2 AU) the α values are practically zero because of the resistive damping. Starting from 2 AU α rises until it levels off at around 3 AU. From 3 to 8 AU we obtain a radial α profile which can be approximated by a √ r dependence. Beyond 9 AU α is again close to zero because of the damping applied there. We mark three regions in radius ( Fig. 1, top right, green lines). Region A, extending from 1 to 3 AU, is affected by the buffer zone. Region B, ranging from 3 to 8 AU shows the √ r slope. Region C, covering 8 to 10 AU, is again affected by the buffer zone. In the following analysis, we will therefore concentrate on region B.
In order to have a radial force-free accretion disk, fields have to drop radially as B ∝ r −1 (Fig. 12, top left). This was also observed for magnetic fields in galactic disks (Beck 2001) (Fig. 1). If the most important toroidal field follows ∝ r −1 , the radial Lorentz force vanishes: In the case of ∂ log ρ/∂ log r = −1.5 and ∂ log c s /∂ log r = −0.5 the α value, dominated by the Maxwell stresses will then scale as √ r, which is actually matching the value that we measure in our best resolved model (see Fig. 1, top right).

Mass loss
The models with open boundaries show a considerable mass loss over the course of our simulations. A vertical outflow removes a substantial amount of mass. The total mass loss over time is presented in Fig. 5. The mass loss is determined in space region B. The closed models FC and PC loose there mass due to radial mass movement. The open models loose there mass mainly due to the vertical outflow. We will discuss this outflow in section 3.5 . To check the possible impact of this mass loss onto the properties of the turbulence, we restarted run FO after 1000 inner orbits with the current velocity and magnetic field configuration but the initial density distribution. We call this model FOR (FO Restarted, Table 1). After restarting the simulation, the turbulence needs a couple of inner orbits to readjust the fully turbulent state. We compared the mean total α stress of the runs FO and FOR and found a comparable evolution of the α values. We measure α = 1.4 · 10 −3 for FOR and α = 1.3 · 10 −3 for FO in the time period from 1000 up to 1400 inner orbits. We conclude that the mass loss is not yet influencing the development and strength of turbulence.

Viscous disk models
The classical α viscous disk model should reproduce the radial mass flow as it occurs in global MHD simulations of MRI turbulent disks. Balbus & Papaloizou (1999) have argued that the mean flow dynamics in MRI turbulence follows the α prescription. To further test this supposition, we performed a series of 2D HD viscous comparison simulations with the PLUTO code for several of our 3D MHD runs. We use the same resolution and the same initial setup. Yet, the magnetic field evolution is now replaced by an explicit shear viscosity from the time averaged radial α profile ν(R) = α(R)Hc s obtained from the MHD simulations (see Fig. 1 and find a value very close to the time dependent viscous run ofṀ = 6 · 10 −9 M ⊙ yr , based on a surface density profile of Σ g = 524 · R −0.5 [g/cm 2 ] and our disk parameter H = 0.07 * R. In Fig. 2, top right, we show the time and azimuthal average of the accretion rate over radius and height. There is a dominant inward accretion at the midplane (Fig. 2, top right, red color). This result is in contrast to the viscous runs where we see the minimum of accretion and even a small outflow at the midplane (Kley & Lin 1992;Takeuchi & Lin 2002). After Takeuchi & Lin (2002) (eq. 8) there are several possibilities which could change the vertical profile of the radial velocity, and therefore the mean accretion flow. Radial and vertical gradients in the orbital frequency as well as a spatially varying α will affect the vertical profile of the meridional outflow. For the MHD simulations one has to include the vertical gradient as well as the time derivative of the orbital frequency. Fig. 2, bottom right, demonstrate the change of the orbital frequency with a period of around 50 local orbits at 5 AU.
Radial mass flow and surface density evolution have shown that we can fit our MHD global models with a viscous disk model as long as we use an α profile compatible with our MHD run.
Of course, the disk spreading that we observe in our MHD run is partly due to the existence of our radial buffer zones, in which not only the fields decay, but also the α stresses vanish. In a larger radial domain we can expect that also a larger region of the disk will get into a steady state of accretion. However, one could also argue that in a realistic protoplanetary accretion disk one will ultimately reach dead zones which behave similar as our buffer zones. In that sense the active part of our global disk is embedded between two dead zones.

Vertical outflow
In the previous section, we have seen that the MHD simulations point to the presence of an additional process besides turbulent "viscous" spreading that removes gas from the disk. Fig.   4, right, shows the initial vertical density profile at 5 AU (dotted line), the respective profile for the MHD model BO (solid line) and the viscous HD run (dashed line) after 1000 inner orbits (90 local orbits). For our model, the vertically integrated optical depth of τ = 1 is around 2.8 scale heights (Fig. 4, red line). For the calculation of the optical depth, we used Rosseland mean opacities (Draine & Lee 1984) with the temperature at 5 AU. The additional magnetic pressure as well as the vertical mass flow, present in the MHD run, generates higher gas density above 3 scale heights compared to the hydrostatic equilibrium. In Fig. 4, left, we plotted a snapshot of the azimuthally averaged velocity field, taken after 1000 inner orbits. The plot indicates a vertical outflow in the disk starting above 2 scale heights. We measured the angle between the cylindrical radial velocity V R and the vertical velocity V Z for the mean and turbulent components (Fig. 3,   left). The angle is measured with respect to the midplane axis (pointing to the star, see Fig. 3, left, V R = −1 and V Z = 0). The upper (red solid line) and the lower hemisphere (blue solid line) present similar profiles. From the midplane up to 1.8 scale heights, the turbulent velocity field is directed upwards but still pointing to the star. The low angle of 10 o for the mean velocities shows the gas motion pointing to the star and towards the midplane. At 1.8 scale heights the turbulent velocity is pointing away from the midplane and the star. Also the mean velocity angle changes quickly in the region between 1.6 and 2 scale heights to an outflow angle, e.g., steeper than 90 o . This region coincides with the region where the vertical outflow is launched (Suzuki et al. 2010).
Above 2 scale heights the angle between the turbulent and the mean velocity components stays above 90 o , leading a vertical outflow with a small radial outward component. The so-called dynamical evaporation time is the time to evacuate the gas completely from the disk assuming no supply of matter. In our model the value is slightly larger than 2000 local orbits (Fig. 3, right) which provides a confirmation for the vertical outflow obtained by local box simulations from Suzuki et al. (2010) with a vertical net flux field.
In Fig. 3, right, we plotted the vertical mass flow over height at 5 AU. The outflow starts at 2 scale heights and reaches mass fluxes of 10 −10 M ⊙ /yr at 5 AU (model BO, solid line). The influence of the pure outflow boundary is observed to cause the small outflow in the HD viscous run (dashed line). In the midplane region, the disk reestablishes the hydrostatic equilibrium due to the radial mass loss at the midplane (Fig. 6, bottom right, red solid line). This drives to a small mean vertical motions visible in the vertical velocity (Fig. 6, bottom right, green dotted line).
The gas does leave the grid with Mach numbers of only 0.5, which is significantly lower that the local escape velocity, which would be about Mach 20. Even the results indicate a stable vertical outflow, without including the sonic point and the Alfvénic point in the simulation it is not possible to make prediction about the flow, leaving or returning to the disk at larger radii. Thus the fate of the vertical outflow to be a disk wind or not will have to be determined in more detail in future simulations with a much broader vertically extent. For this study one would most probably need an additional vertical field in the corona, which could support additional propulsion effects like magneto-centrifugal acceleration.

Velocity analysis
Planet formation processes in circumstellar accretion disk are strongly dependent on the strength of the turbulence. Turbulence mixes gas and particles, diffuses or concentrates them and makes them collide (Ilgner et al. 2004;Johansen & Klahr 2005;Johansen et al. 2007;Brauer et al. 2008;Cuzzi et al. 2008;Carballido et al. 2010;Birnstiel et al. 2010). The property of MHD turbulence that is important for planet formation are the turbulent velocity and density fluctuations of the gas. The density fluctuations are around 10% and follow the results by Fromang & Nelson (2006). The spatial distribution of the fluctuating and mean part of the velocities is presented in A global picture of the total rms velocity is presented in Fig. 7. The 3D picture is taken after 750 inner orbits and shows again the different turbulent structures of the midplane and coronal region. There are also localized supersonic turbulent motions in the disk corona (Fig. 7, white color). Compared to the turbulent velocity, the mean velocities of the gas are two order of magnitude smaller. They show small but steady gas motions in the disk. The vertical dependence for the mean velocity (Fig. 6 right -bottom) shows the small inward motion (red solid line) as well as the change of r and θ-velocity components to an outflow configuration around 1.6 scale heights.

Kinetic spectra
Not all particles do couple alike to the turbulent gas flow. In fact particles have a sizedependent friction or stopping time (Weidenschilling 1977). This stopping time is also the time a particle needs to couple to the turbulent gas flow. Best coupled to turbulence are those particles, which have a coupling time shorter then turbulent correlation time. Particle collision velocities are maximized for particles whose stopping time coincides with the turbulent correlation time, e.g., the eddy turn over time. This means that particles of different sizes couple to different length scales of the turbulent spectrum. Therefore, a study of planet formation processes needs not only the mean turbulent velocity but also its spectral distribution.
In the global domain, only the k φ space of the spectra is without modifications accessible as only the φ direction is periodic in space. The investigation of the complete k-space for this model goes beyond the scope of this work. The classical Kolmogorov theory predicts the scaling of the energy spectra per wavenumber: Fig. 1) and height (±0.5 disk scale heights). For our spectra we use the azimuthal wavenumber m instead of k to be independent from radius: k = 2π/λ = m/R. In Fig.   6, bottom left, we plotted the energy power spectra E m = v(m) 2 · m with time and space averaged over ±0.5 scale heights around the midplane. In our models we do not observe Kolmogorov inertial like range, E m · m 2 ∼ m 1/3 . The 2π runs F0 and BO have most of the energy placed at m = 5. The high resolution model BO present a k −1.2 dependence, starting from m = 5 until m = 30. The π/4 run PO piles up the energy at its domain size (m = 8), reaching higher energy levels compared to the 2π models FO and BO. The velocity spectra for each component along the azimuth, plotted in Fig. 8, left, indicate that all velocity components have similar amplitude for the small scales and do not deviate by more than a factor of two at the largest scales. The radial velocities peak at m equals 5, but overall the entire spectrum above m=20 is essentially flat.
The peak at m=5 could be connected to the production of shear waves in the simulations. These shear or density waves are described in Heinemann & Papaloizou (2009). On top of the shorter time scale of MRI turbulence, these long "time scale" shear waves are visible in the contour plot of the radial velocity in the r − φ midplane (Fig. 8 right). The shear wave structures drive the radial velocity up to 0.3c s . In the velocity spectra we see the start of the dissipation regime at m = 30 − 40 for the high-resolution run BO. For the model BO, this corresponds to 26 or rather 19 grid cells per wavelength, which is still well resolved by the code .
Shearing waves are also visible in a r − θ snapshot of the velocity (Fig. 9, left). Here we plot the azimuthal velocity V φ − V K as contour color, over-plotted with the velocity vectors. Red contour lines show Keplerian azimuthal velocities. Super-Keplerian regions are important for dust particle migration. They reverse the radial migration of particles, leading to their efficient concentration and triggering parasitic instabilities in the dust layer, like the streaming instability leading potentially to gravoturbulent planetesimal formation Johansen et al. 2007). In our simulation these super-Keplerian regions are not completely axisymmetric, but have a large extension in the azimuthal direction of several scale heights. The variation of the orbital frequency over time and space, presented in Fig. 9, left, and

Magnetic field analysis
The azimuthal MRI generates a turbulent zero-net field configuration in the disk. Despite the loss of mass and magnetic flux through the vertical boundary there is no sign of decay for the highest resolution case BO (Fig. 1, top left, bottom right). We find well established turbulence. Fig. 9, right, presents a snapshot of the magnetic fields after 750 inner orbits. The r − θ components are shown as vectors with the azimuthal magnetic field as background color.

Magnetic energy spectrum
To understand the magnetic turbulence at the midplane, we investigated the spectral distribution of the magnetic energy. The magnetic energy power spectrum (Fig. 12 -bottom left) is plotted along the azimuthal direction with the same time and space average as for the kinetic energy power spectra. We plot the magnetic energy power spectra times the wave-number m · B 2 m /2P Init−5AU to show where most of the magnetic energy is located. For all runs, most of the magnetic energy is deposited in small scale magnetic turbulence. This was found in several recent MRI simulations, latest in local box simulations by Davis et al. (2010) and (Fromang 2010). The peak of the magnetic energy lies just above the dissipation regime. For the 2π model FO, the peak is located between m = 10 and m = 20, whereas for the high-resolution run BO this regime is shifted to m = 20 and m = 30. The spectra follows closely the m 1.0 slope until the dissipation regime is reached. The π/4 run does not resolve the scales where we observe this m dependence. In the restricted model PO most of the magnetic energy is again located at the scale of the domain size.

Convergence
The convergence of MRI is an important aspect in ongoing MRI research in local and global simulations. In local boxes, there was found convergence for the large scale turbulence between 32 and 64 grid cells per scale height (Davis et al. 2010). Due to the large domain in global simulations, it was up to now not feasible to reach such resolutions per scale height. Here the first resolution level is needed to reach a self sustaining turbulence, at least for simulations with a zero-net flux toroidal field (Fromang & Nelson 2006) Doing this, we expect only a weak decrease for the large scale modes, as presented in Davis et al. (2010) , Fig. 3.

Plasma beta
The overall strength of the magnetic fields is best analyzed by this plasma beta value β = 2P/B 2 . Fig. 10 presents a 3D picture of the logarithmic plasma beta for the r − θ components, taken at 750 inner orbits. The two-phase structure of the disk is again visible. The well established turbulence at the midplane has a broad distribution of high plasma beta values (Fig. 11 , top left).
In contrast, there are regions in the corona of the disk with plasma beta below unity (Fig. 10, black regions). The azimuthal and time averaged plasma beta at the midplane lies around 400 (Fig. 11, bottom right). In Fig. 11, top left, we plot the correlation of plasma beta over height in a scatter plot of all grid cells. We find the distribution of beta values to be very narrow in the disk corona (1-10) but on the other hand to be much broader (10 -10 4 ) around the midplane, but strongly peaked around β = 500. The value of plasma beta in the disk corona depends on several issues. A zero-net flux MRI turbulence with toroidal field produces lower magnetic fields in the corona. This was already shown in a very similar simulation by Fromang & Nelson (2006) (Fig.   8, solid line, model S2). In contrast, a vertical initial field produces a stronger turbulence level with plasma beta values below unity in the corona. The boundary condition also affects the values in the corona. A closed boundary condition, e.g. periodic in the vertical direction will accumulate large amount of magnetic flux in the corona and drive to a plasma beta value smaller then one (observed in model FC and PC). The small increase of plasma beta above 3 disk scale heights is connected to the vertical outflow and the increase of gas pressure and density in this area (Fig. 4,right). This effect has to be investigated in future work with a much broader vertical extent. Very high plasma beta values in the midplane (Fig. 11, top left) indicate reconnection. Two magnetic fields with different sign and comparable strength coming too close to each other, e.g., in the same grid cell, do reconnect. Such reconnections are visible in single grid cells with nearly no magnetic field. For our BO model, the reconnection zones reach plasma beta values up to 10 11 . The heating due to reconnection in those regions is not covered in our isothermal model, but shall be a subject for future studies.

Spatial distribution
As we already mentioned in section 3.2, the radial profile of the turbulent magnetic field has a direct effect onto the radial profile of the Maxwell stress in the α parameter. The dominant turbulent azimuthal magnetic field goes as 1/r, as shown in a azimuthal and time average in Fig.   12, top left. The saturated turbulent field is 4 times lower than the initial azimuthal field. All values are normalized to the initial gas pressure at 5 AU at the midplane and the radial profiles are again mass weighted. The vertical profile shows a constant distribution around ±2 scale heights from the midplane until it decreases with height (Fig. 12, top right). In contrast, the radial and θ component show a local minimum at the midplane with a peak of turbulent magnetic field slightly above 2 scale heights.
The turbulent magnetic fields are around 2 orders of magnitude larger than the mean fields. The vertical profiles of mean magnetic fields over height are presented in Fig. 12, bottom right.
The radial magnetic field is anti-symmetric to the midplane and correlated with the dominating azimuthal component. The distribution of mean magnetic fields are connected to the "butterfly" oscillations.

Butterfly structure
The butterfly pattern is a general property of MRI turbulence and was found in many local and global simulations, latest by Gressel (2010), Flaig et al. (2010), Sorathia et al. (2010) and Dzyurkevich et al. (2010) . The "butterfly" pattern becomes visible for the mean B φ evolution, plotted over disk height and time. In Fig. 13, bottom, we plotted the B φ component of the magnetic field averaged over a small radius (4 -5 AU) and over azimuth for model FO, left, and PO, right. We see a clear "butterfly" pattern in both models. This pattern is also visible in the total accretion stress with doubled period (Fig. 1, bottom right). In comparison, the π/4 run shows no systematic and more violent picture of the butterfly. The amplitudes are stronger and it has mixed symmetry (Fig. 13, bottom right). Also the total magnetic flux evolution shows these properties for model PO (Fig. 13, top right). The FO run presents a similar amplitude and period as the BO run. The effect of the narrow azimuthal domain on the mean fields will be investigated in a follow-up work. The reason of this butterfly structure and its role for the MRI is still under discussion. Recent studies show the connection to the MHD dynamo (Gressel 2010) and magnetic buoyancy (Shi et al. 2010).

Discussion
For a number of aspects, our Godunov method confirms results previously obtained with a finite difference method as presented by Fromang & Nelson (2006).
• A minimum amount of grid cells per scale height, which is about 25 grid cells, is needed to sustain turbulence, which was to be expected as both methods present similar numerical behavior . Otherwise the turbulent magnetic energy slowly decays in the nonlinear MRI evolution (Fromang & Nelson 2006). Our highest resolution model BO was able to sustain a constant level of turbulent stress for more than 400 inner orbits.
• The toroidal magnetic net flux is quickly lost via an open vertical boundary. Then, there is a oscillating zero-net flux field present in the disk.
• The disks show a two layer structure of turbulence.

α-stress evolution
We obtained a steady state α value of about 5 · 10 −3 , which is comparable to the results obtained in Fromang & Nelson (2006). The time averaged radial profile of α follows √ r. This profile can be explained by the choice of our radial pressure (density and temperature) profile, in combination with the resulting magnetic field profile which is force free, |B ′ φ | ∝ r −1 . For this magnetic field profile, the net radial magnetic force vanishes. Any quasi steady state of disk turbulence must display this profile, otherwise large scale radial readjustments in the density profile would occur. Both B ′ φ and B ′ r determine the Maxwell stress B ′ φ B ′ r to be ∼ 1/r 2 . For the chosen ∂ ln P/∂ ln R = −2.5 this results in α ∼ B ′ φ B ′ r /P ∼ √ r.
Of course, this profile is only valid for well-ionized accretion disk regions. The radial α profile in protoplanetary disks remains an open question. The ionization rate, and possibly the MRI activity, will be a function of radius and height (Semenov et al. 2004;Dzyurkevich et al. 2010).
Furthermore, the pressure scale height will vary with radius, and this also changes the MRI evolution. Both effects will lead to different saturation levels for MRI turbulence at different radii and thus also to different α values. We can already conclude that our global magnetic energy spectra as well as the effect of increasing resolution (Fig, 12, bottom left, FO to BO) look very similar to the results presented in local box simulations. The magnetic energy power spectra reveals that most of the magnetic energy is deposited in the small turbulent scales. For the model with a restricted azimuthal domain of π/4, the largest energy scale is always the domain size.

Turbulent velocity
Recently observed turbulent velocities in TW Hya and HD 163296 (Hughes et al. 2010) fit nicely to our computed Mach numbers of 0.1 and 0.4 in the midplane and corona of the disks.
The kinetic energy spectra as well as the velocity spectra along azimuth (k φ space) show a peak between m = 5 and 6 due to radial shear waves. Latest results in local box MRI simulations presented a k −3/2 slope for the kinetic energy power spectra (Fromang 2010). Similar as in Fromang (2010), our power-law fit applies only for a small range in k space. We do not find a Kolmogorov type slope of k −5/3 (Fig. 6). Further studies are needed, including the k r and k θ .
A Kolmogorov scaling was predicted for magnetic ISM turbulence by Goldreich & Sridhar (1995) and recently confirmed in numerical simulations by Beresnyak (2010). However, it only applies for the inertial range of incompressible isotropic turbulence. The driving of the turbulence via MRI, the anisotropy of the turbulent eddies, the geometry and rotation of the disk and the compressibility of the gas make it difficult to argue for a Kolmogorov scaling. We expected a spectrum to be more or less deviating from this simple law.

Two-phase disk structure
We observe that the accretion disks establish a two-phase structure: -The midplane region between ± two scale heights shows a pretty constant turbulent RMS velocity of about 10% of the local sound speed independent of radius or height. Part of the RMS velocity occurs due to global shear waves which have radial peak velocities of up to 30% of the local sound speed. The amplitude of the azimuthal fluctuations in the magnetic field is also independent of height (± two scale heights around the midplane) but develops a 1/r profile in radius. The midplane region shows a broad distribution of plasma beta values, β = 2P B 2 , with a mean value of about 500 and a full width at half maximum of two order of magnitude.
-In the coronal region, more than two scale heights above the midplane, the mean turbulent velocity reaches a Mach number of 0.5 with supersonic peaks up to 1.5. The mean magnetic fields decrease in this region with height. The disk corona shows a narrower distribution of the plasma beta values with most values between 1 and 10. Here the magnetic fields are buoyant, gas and fields are expelled from the disc. Relative high plasma beta (β > 1) in the corona have been reported in Fromang & Nelson (2006) for global models of AMRI with open boundary.
The magnetic flux escapes through the vertical boundary with a remaining zero-net flux in the computational domain. This leads to the weakly magnetized corona (below equipartition).

Vertical outflow
Our models show a MRI driven vertical outflow. Above 2 scale heights, the gas flow is directed vertically and radially outward, Fig. 3. The outflow velocity of the gas (measured at the vertical boundary) is still subsonic. The disk evaporation time was determined at 5 AU to 2000 local orbits. The launching region is located between 1.6 and 2 scale heights. This results matches values obtained in local box simulations (Suzuki & Inutsuka 2009;Suzuki et al. 2010) with a vertical net-flux field.
However, we are aware that a detailed study of the vertical outflows requires much broader vertical extended simulations to confirm that the gas is evacuated from the disk and not returning. These simulations should then include the sonic point or even the Alfvénic point to give further insight into disk-wind and disk-jet interacting regions.

Meridional flows
Our present work shows that the meridional outflow at the midplane is only present in HD simulations, e.g., in viscous simulations with an α value assumed to be constant in time and space.
For our MHD models, we find time variations of the orbital frequency of around 50 local orbits, which are not present in the viscous disk models and which prevent a steady radial outflow. A similar result, the absence of a meridional flow in global MHD simulations was recently found by Fromang et al. (2011).
However, we confirm the more general picture of a viscous disk and show that viscous disk models with a radial viscosity profile can reproduce successfully the radial mass flow rate in global MRI turbulent stratified disks. Clearly, the vertical mass flow cannot be described with such an HD model.

Mean field evolution
The azimuthal MRI is self-sustaining in our zero net flux simulations with open boundaries.
The fact that the total flux oscillates around zero could be due to the generation of a mean poloidal magnetic field by a turbulent toroidal field. We observe also an antisymmetric distribution of the mean magnetic fields in the upper and lower hemisphere which could be an indication for the action of an MHD dynamo in our global simulations.
The existence of an α-Ω MHD dynamo and its role for accretion disks was investigated by Brandenburg et al. (1995), Ziegler & Rüdiger (2000), Arlt & Rüdiger (2001), Brandenburg & von Rekowski (2007), Lesur & Ogilvie (2008) and Blackman (2010). The temporal oscillations of the mean azimuthal field, plotted over height and time (Fig. 13), generates a butterfly pattern. Latest results connect the butterfly pattern with a dynamo mechanisms (Gressel 2010). We present also a butterfly pattern with a period of 10 local orbits, independent of the azimuthal extent. Additionally, the butterfly structure is reflected in the temporal spatial fluctuations of the mean turbulent stresses with double period. A change of sign of the mean azimuthal field occurs every five local orbits, at the same time the α-stresses show a minimum ( Fig. 1, bottom right).
The magnetic energy as well as the mean field evolution have shown that the π/4 model does not capture the correct properties of the larger scale simulations. Hawley (2000) also studied full 2π and restricted π/2 models of accretion tori. However, a detailed study of the impact of different azimuthal domain extents is still needed and will be covered in future work.

Summary
We have performed full 2π 3D stratified global MHD simulations of accretion disks with the Godunov code PLUTO. Our chosen disk parameter represent well-ionized proto-planetary disk regions. We obtain a quasi steady state zero-net flux MRI turbulence after around 250 inner orbits.
• The second order Godunov scheme PLUTO including the HLLD Riemann solver presents a similar nonlinear MRI evolution as finite difference schemes. There is also a need of about 25 grid cells per pressure scale height to reach a self-sustaining MRI turbulence in global zero net flux azimuthal MRI simulations.
• We observe a total α parameter of about 5 · 10 −3 , which remains constant for at least 400 inner orbits and scales with √ r for our used pressure and density profile.
• The turbulent magnetic fields show a 1/r profile in radius, mainly visible in the dominating toroidal magnetic field. This configuration is force free in the sense that there exist no large scale net force on the gas. This profile determines the slope of the α parameter.
• The magnetic energy spectra is similar as in local box simulations. Most of the magnetic energy is placed in the smallest resolved turbulent scale.
• The kinetic energy spectra as well as the velocity spectra peak for an azimuthal wavenumber between m = 3 and 5 due to shear waves, driving the radial velocity up to a Mach number of 0.3. We do not find a Kolmogorov type scaling in the k φ space.
• The model with an azimuthal extent of only π/4 has most of the energy at the domain size and does not show the same mean field evolution.
• We observe a butterfly pattern with then local orbits independent of the azimuthal extent.
The butterfly period becomes also visible in the Maxwell stress with double period. The mean magnetic fields are antisymmetric for the two hemispheres.
• At the midplane (±2 disk scale heights), our turbulent RMS velocity presents a constant Mach number of 0.1 independent on radius. At the corona (> 2 disk scale heights), the turbulent velocity increases up to a Mach number of 0.5 at 4 scale heights.
• The turbulent magnetic fields at the midplane present a broad plasma beta distribution with a mean of about 500 ± one order of magnitude. In the corona the plasma beta is between unity and ten.
• The turbulent and the mean velocities are pointing vertically and radially outward in the disk corona (> 2 disk scale heights). We observe a steady vertical outflow for the open boundary models, dominating the radial accretion flow.
• We do not see a meridional flow pointing radially outward at the midplane in our MHD models. However, we reproduce our total radial mass flow in 2D viscous disk simulations with a radial dependent α-viscosity.

Outlook
This paper presents a huge data set of about 10 TBytes. This means we will continue to analyze the data for different goals. One study will deal with a closer investigation of dynamo properties in our global disk models. Another one will analyze the turbulent spectra in a better way to derive correlation times for the turbulence. We will also fill the parameter space with π and π/2 models to identify whether a subsection will be sufficient. Higher resolution is envisioned to reach resolution per scale height comparable to recent stratified local box simulations. Finally, our global MHD model will be the work horse for our future investigations of planet formation processes in circumstellar disks, like collisions of boulders, planetesimal formation and planet migration.
In future runs, we also plan to use non-ideal MHD to include more realistic magnetic Prandtl numbers and magnetic Reynolds numbers to understand the occurrence and saturation level of the turbulence. Improving the thermodynamics is also a must in future work, dealing with the proper ionization of the disk, like capturing p∆V terms or magnetic dissipation as heat input.
We thank Andrea Mignone for very useful discussions about the numerical setup. We thank Sebastien Fromang for the helpful comments on the global models and on the manuscript. We