Global MHD Simulations

# Global Evolution of an Accretion Disk with Net Vertical Field: Coronal Accretion, Flux Transport, and Disk Winds

Zhaohuan Zhu11affiliation: Department of Physics and Astronomy, University of Nevada, Las Vegas, 4505 South Maryland Parkway, Las Vegas, NV 89154, USA , AND James M. Stone22affiliation: Department of Astrophysical Sciences, 4 Ivy Lane, Peyton Hall, Princeton University, Princeton, NJ 08544, USA
###### Abstract

accretion, accretion disks - dynamo - diffusion - magnetohydrodynamics (MHD) - instabilities - turbulence - protoplanetary disks -

## 1. Introduction

Accretion and outflow have been observed in a wide range of astrophysical systems, ranging from protoplanetary disks (Hartmann, 1998) to supermassive black holes (Begelman et al., 1984). Two mechanisms are likely to drive accretion in most systems: MHD turbulence or a magnetized disk wind. Turbulence can lead to net stress within the disk that transports angular momentum radially, while a disk wind can lead to a stress at the disk surface which carries angular momentum away vertically. It is thought the main mechanism driving turbulence is the magnetorotational instability (MRI, Balbus & Hawley 1991, 1998), while in a Newtonian potential the main mechanism which produces a wind is the magnetocentrifugal effect associated with vertical fields (Blandford & Payne, 1982). Since many astrophysical disks are poorly ionized (e.g. protoplanetary disks), the effects of non-ideal MHD on MRI (e.g. the review by Turner et al. 2014b) and disk winds (Konigl, 1989; Wardle & Koenigl, 1993; Königl et al., 2010; Salmeron et al., 2011) is also important in these systems.

Outflow/wind/jet launching mechanisms require large scale poloidal magnetic fields. In the magnetocentrifugal wind model, the magnetic field is anchored in the rotating disk. If the poloidal component of the field makes an angle of less than 60 with the disk surface, the centrifugal force can break the potential barrier and accelerates matter outwards, leading to outflow. The stress at the wind base launches the outflow and, at the same time torques the disk, leading to accretion below the wind region. Such a picture has been confirmed by numerical simulations with prescribed poloidal magnetic fields (e.g. the review by Pudritz et al. 2007).

Net vertical magnetic fields also play an important role in MRI turbulence. It is known that the turbulent stress increases as the net field increases (Hawley et al., 1995), thus net field promotes both outflow and disk accretion.

To address the relative importance of turbulence and a disk wind in driving accretion, we need to rely on numerical simulations. However, simultaneously resolving small scale turbulence and capturing the large scale disk wind is challenging. In order to resolve the MRI, most previous simulations with net vertical flux only study a small patch of the disk using the shearing box approximation (Suzuki & Inutsuka, 2009; Bai & Stone, 2013; Fromang et al., 2013). However, the resulting wind in such simulations can flow in either radial direction (Bai & Stone, 2013; Lesur et al., 2014), and the outflow rate drops dramatically when a taller box has been used (Fromang et al., 2013). Global simulations of turbulence within an accretion disk with net vertical fields have been carried out recently (Suzuki & Inutsuka, 2014; Gressel et al., 2015) and strong outflows and high disk accretion rates have been observed. However, these simulations have limited radial range, only cover a wedge around the disk midplane , and have not been evolved long enough for a true steady-state to emerge. On the other hand, global (usually two-dimensional) simulations of disk winds that do not capture the MRI in the disk interior have been reported (Stone & Norman, 1994; Ouyed & Pudritz, 1997a; Krasnopolsky et al., 1999, 2003; Kato et al., 2002; Porth & Fendt, 2010; Ramsey & Clarke, 2011). In such work, either only the wind region has been studied or explicit resistivity has been assumed in the disk region as a ”sub-grid” of the MRI (e.g. Casse & Keppens 2002, 2004; Fendt & Čemeljić 2002; Zanni et al. 2007; Tzeferacos et al. 2009, 2013). In this work, we carry out global three-dimensional simulations which capture both the MRI and disk wind self-consistently by using both mesh-refinement and special polar boundary conditions.

An important question that can only be addressed in global simulations is the net rate of transport of vertical magnetic field. Net vertical fields are both advected by the large-scale flow, and diffuse due to turbulence. In the phenomenological model by van Ballegooijen (1989) disk accretion and field diffusion is modeled using turbulent viscosity and resistivity, however the result depends sensitively on the relative amplitude of both. As clearly summarized in Guan & Gammie (2009), the evolution of the poloidal fields is governed by

 ∂tAϕ= −vRBz−ηR∂zBR+η∂RBz ∂tAϕ= vR1R∂R(RAϕ)+η∂2zAϕ+η∂R[1R∂R(RAϕ)]. (1)

if only radial motion is considered. In Equation 1, is the azimuthal component of the vector potential. On the right hand side, the first term is advection by the accretion flow. Since in the viscous model, this term is roughly . The second term is the vertical diffusion of radial field and is roughly . If we assume the large scale field enters the disk at an angle of , as in the disk wind model, we have . The third term is radial diffusion, roughly . Thus, the first and second terms dominate, and they balance each other when or the turbulent magnetic Prandtl number (Pr) . However, local shearing box simulations suggest that Pr1 (Guan & Gammie, 2009; Lesur & Longaretti, 2009; Fromang & Stone, 2009), implying that large-scale fields will diffuse outwards faster than the inward advection (Lubow et al., 1994). To maintain large scale magnetic fields, either turbulent diffusion is significantly reduced (Spruit & Uzdensky, 2005; Rothstein & Lovelace, 2008) or inward accretion is increased. An important insight is provided by Beckwith et al. (2009) (see also simulations by Stone & Norman 1994) which observed that magnetic flux is mainly transported in the corona of the disk and magnetic fields are pinched within the corona. Thus, the above 1-D estimates may not apply in more complex multi-dimensional field geometries. Since the radial distribution of global magnetic fields also determines the collimation of disk winds (Anderson et al., 2005; Pudritz et al., 2006), the processes which determine and maintain large scale magnetic fields in disks are essential for sustaining outflow or even accretion.

In §2, the theoretical framework on describing turbulence and disk wind is presented. Our numerical method is introduced in §3. The results are presented in §4. After a short discussion in §5, we conclude in §6.

## 2. Theoretical Framework

Which is more important for disk accretion: turbulence or outflow? If we average the angular momentum equation in the azimuthal direction and integrate it in the vertical direction, we can derive

 ∂∫⟨ρvϕ⟩dz∂t =−1R2∂∂R(R2∫(⟨ρvRδvϕ⟩−⟨BRBϕ⟩)dz) −12πR2∂Rvk˙Macc∂R−(⟨ρvzvϕ⟩−⟨BzBϕ⟩)∣∣∣zmaxzmin (2)

or

 ∂∫⟨ρδvϕ⟩dz∂t =−1R2∂∂R(R2∫(⟨ρvRδvϕ⟩−⟨BRBϕ⟩)dz) −˙Macc2πR2∂Rvk∂R−(⟨ρvzδvϕ⟩−⟨BzBϕ⟩)∣∣∣zmaxzmin (3)

for the perturbed quantities (). The symbol denotes that the quantity has been averaged over the direction 111 has been assumed to be constant along . Without this assumption, there will be an additional term related to ., and . We refer to and as the radial Reynolds and Maxwell stress respectively, to distinguish them from the vertical stresses . If we normalize the stresses with pressure, we can identify different contributions to the total stress as

 αRϕ,Rey=⟨ρvRδvϕ⟩/⟨p⟩andαRϕ,Max=−⟨BRBϕ⟩/⟨p⟩ αzϕ,Rey=⟨ρvzδvϕ⟩/⟨p⟩andαzϕ,Max=−⟨BzBϕ⟩/⟨p⟩.

Stresses and the parameters in spherical-polar coordinate can be defined in similar ways. Since the vertically integrated stress determines the disk accretion rate as in Equation 3, we can define the vertically integrated parameter as

 αint=∫TRϕdzΣc2s. (4)

where is the sum of both radial Reynolds and Maxwell stress. If we choose so that =0 and assume that the magnetic stress dominates at the disk surface, Equation 3 can be written as

 ˙Macc=−2π∂Rvk/∂R(∂∂R(R2αRϕ,intΣc2s)−R2⟨BzBϕ⟩∣∣∣zmaxzmin) (5)

Equation 5 suggests that disk accretion is due to both within the disk and exerted at the disk surface. Normally, the stress is from disk turbulence, and the stress at the surface is from a magnetocentrifugal disk wind. If both the and stresses have similar values, the second term on the right will be larger than the first term by a factor of which can be quite large for a thin disk . On the other hand, with vigorous turbulence the internal stress may be larger than the surface stress, thus it is unclear if disk accretion is driven by turbulence or wind when both processes are present.

## 3. Method

We solve the magnetohydrodynamic (MHD) equations in the ideal MHD limit using Athena++ (Stone et al. 2017, in preparation). Athena++ is a newly developed grid based code using a higher-order Godunov scheme for MHD and the constrained transport (CT) to conserve the divergence-free property for magnetic fields. Compared with its predecessor Athena (Gardiner & Stone, 2005, 2008; Stone et al., 2008), Athena++ is highly optimized for speed and uses a flexible grid structure that enables mesh refinement, allowing global numerical simulations spanning a large radial range. Furthermore, the geometric source terms in curvilinear coordinates (e.g. in cylindrical and spherical-polar coordinates) are carefully implemented so that angular momentum is conserved to machine precision, a crucial feature to enable the angular momentum budget analysis as presented in §4.1.

Since a disk wind flows radially, we adopt a spherical-polar coordinate system (, , ) for our simulations, which should minimize the effects of the domain boundary on the wind properties (Ustyugova et al., 1999). We implement a special polar boundary in the direction allowing our simulation domain to extend all the way to the pole, which is different from previous similar simulations where a hole was carved out close to the pole. This boundary condition prevents the loss of magnetic fields at the previously carved-out pole, which is crucial to study transport of magnetic fields in disks. The details on the implementation are given in the Appendix.

Although we adopt spherical-polar coordinates for the simulations, we transform some quantities to cylindrical coordinates to study physical processes in disks. In this paper, we use (, , ) to denote positions in cylindrical coordinates and (, , ) to denote positions in spherical polar coordinates. In both coordinate systems represents the azimuthal direction (the direction of disk rotation). Our simulation domain extends from to with from 0 to . extends from 0 to except for the thin disk case for which is from 0 to .

The initial density profile at the disk midplane is

 ρ0(R,z=0)=ρ0(R0,z=0)(RR0)p, (6)

while the temperature is assumed to be constant on cylinders

 T(R,z)=T(R0)(RR0)q. (7)

Hydrostatic equilibrium in the plane requires that (e.g. Nelson et al. 2013)

 ρ0(R,z)=ρ0(R,z=0)exp[GMc2s(1√R2+z2−1R)], (8)

and

 vϕ(R,z)=vK[(p+q)(csvϕ,K)2+1+q−qR√R2+z2]1/2, (9)

where is the isothermal sound speed, , and . A local isothermal equation of state is assumed during the simulation 222To achieve this, we have used the adiabatic equation of state with =1.4, but instantaneous cooling is applied at each timestep. . Using Equations 8 and 9, the density and velocity will become infinite at the pole. To avoid this, we use in the above equations.

We adopt a density floor which varies with position

 ρfl=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩ρfl,0(RR0)p(1z2)ifR≥rminρfl,0(rminR0)p(1z2)ifR3rminρfl,0(RR0)p(1z2)(5−2r−rminrmin)(4rmin−Rrmin+1)ifR

We choose this density floor so that it becomes small at the disk atmosphere and becomes large close to the inner radial boundary.

Initial magnetic field is assumed to be vertical. To maintain , we use the vector potential to initialize the magnetic fields. We set

 Aϕ=⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩12×R×B0(rminR0)mifR≤rminB0Rm0Rm+1m+2+B0rm+2minRm0R(12−1m+2)ifR>rmin (11)

where . Thus, the vertical magnetic fields are

 (12)

With this setup, the plasma at the disk midplane beyond is a constant.

We choose and in our simulation so that the disk surface density . Then is 1. The time unit in the simulation is . The grid is uniformly spaced in ln(), , with 13664128 grid cells in the domain of [ln(0.1), ln(100)][, ][0, 2] at the root level. Thus, equals 0.052 at the root level. We use an open boundary condition in the radial direction, the polar boundary condition in the direction, and periodic boundary conditions in the direction.

We have carried out three simulations, one with and initial plasma =1000 at the disk midplane (the fiducial case), one with and (the weaker field case), and one with and (the thin disk case). For the simulations with , three levels of refinement have been adopted towards the disk midplane, and at the midplane is 0.0065. Since =0.1 at , the disk scale height is resolved by 15 grids at the finest level. For the thin disk case, 4 levels of mesh refinement have been used so that the disk scale height is still resolved by 15 grids. To save computational cost, in the thin disk case only extends from 0 to . The and case is our fiducial model, since the MRI is better resolved at the midplane and the simulation runs longer. We will mainly present results for our fiducial case, but other cases will be discussed briefly in §4.4.

In ideal MHD, the wavelength for the fastest growing linear MRI mode satisfies

 λ=2π√16/15|vA|/Ω=9.17β−1/2H (13)

(Hawley et al. 1995). With =1000, is 0.29 H so that the most unstable wavelength is resolved by 4.5 grid cells at . The density floor parameter is for the fiducial case and for the other cases.

To demonstrate that our fiducial case has reached the necessary resolution to capture MRI, we plot the azimuthally averaged quality factor ( and ) for and in Figure 2. The quality factor is defined as the number of grid cells that resolve the fastest MRI growing mode (Noble et al., 2010),

 Qθ=λθ/(rΔθ) (14) Qϕ=λϕ/(rsinθΔϕ) (15)

where and . The quantities and are the Alfven velocity calculated using the and components of the magnetic field, and is the Keplerian angular velocity at the disk midplane. To get converged results for MRI turbulence, Sorathia et al. (2012) have shown that if , needs to be , and if , while can be smaller ( 6 in their Figure 8). As shown in Figure 2, the quality factor decreases towards the disk midplane despite the fact that it has been boosted by a factor of 8 towards the midplane due to mesh refinement. is larger than 10 in the whole disk (50 at ), and is around 10 at R=1. Thus, our simulation should have adequate resolution to capture MRI turbulence except at the midplane of the inner region ().

## 4. Results

After running for t=45.6 , which is equivalent to 1442 orbits at the inner boundary, our fiducial disk reaches a steady state within R3, i.e. the inner factor of in radius, as evident in Figure 3. A steady state throughout the whole disk cannot be established with our simulation setup since we do not supply new material at the outer boundary. As shown in the panel of Figure 3, the MRI saturates at later times at larger radii. By the end of our simulation, the inner disk reaches a steady state and the disk mass accretion rate is almost a constant within R3, shown in the panel of Figure 3.

Once a quasi-steady state is reached at t=42 , the surface density profile follows , and the midplane profile follows . The midplane stress then follows . The integrated follows , so that the vertically integrated stress follows and the disk accretion rate is a constant with radii. With these values and Equation 5, we can derive =-0.0038, which is consistent with the direct measurement in the upper right panel of Figure 3. We can see that the vertically integrated is much larger than the midplane by a factor of 10. This is because the stress is almost uniform vertically, or even increases towards the disk surface, while the density drops sharply towards the disk surface. If the stress is perfectly uniform vertically to some disk heights, the vertically integrated stress should follow and the midplane and the vertically integrated should have the same slope. But the stress panel in Figure 3 suggests that the vertically integrated stress is much steeper with , which is because the stress is larger at the inner disk’s atmosphere where the fields are pinched as shown in §5.1. Thus, the midplane and the vertically integrated have different slopes. Since 0.5 at R=1, 42 is already longer than the viscous timescale () at R=1, explaining why the disk has a constant accretion rate.

While the disk is losing mass through accretion, it is also losing magnetic field. The net vertical magnetic field decreases with time, shown in the panel of Figure 3. The plasma calculated with the net vertical magnetic field increases from to at the inner disk (the lower center panel). On the other hand, the ratio between the net vertical field and the disk surface density does not vary much over time (the panel). Although it is tempting to interpret this as flux freezing during the accretion process, we will show later that the global field transport is much more complicated.

### 4.1. Accretion Structure

Although the whole disk accretes inwards, the accretion flow has a complicated structure within the disk as shown in Figure 4. The green curves in the figure are the velocity field lines and the white curves are where the plasma =1. The most striking feature is that the disk accretes through the disk atmosphere where the disk is magnetically dominated () as shown in the rightmost panel. Such surface accretion has been seen as early as Stone & Norman (1994) and recently been studied in GRMHD simulations by Beckwith et al. (2009) where it was termed ”coronal accretion”.

The fast coronal accretion at the disk surface carries the magnetic fields inwards, thus pinching magnetic fields at the disk surface. The magnetic field lines are shown in Figure 5. The field lines change from the vertical direction at the disk midplane to the horizontal direction at the disk surface, and then they reverse their directions higher above. The net fields along the height at are shown in Figure 6. In this work, we define the corona region (the yellow shaded region) as the region which is magnetically dominated () and has a negative radial velocity (). Away from the midplane, first becomes negative and then positive due to the radial drag from the surface accretion. Most importantly, the resulting horizontal field lines are stretched azimuthally due to the Keplerian shear, similar to the growth of azimuthal fields in the linear phase of MRI. At , the negative close to the midplane is stretched to the positive and the positive at the upper corona region is stretched to the negative (also shown in the color map of Figure 5). Such net is crucial for maintaining coronal accretion and launching disk wind. Net is always positive but it gets amplified at the upper boundary of the corona region since the field lines are pinched there. Such net field geometry exerts strong stresses or torques to the disk. Since and have opposite signs as shown above, the magnetic stress, -, is positive from the midplane to the wind region. The stress - has similar shapes as - since is always positive. At , the stress changes from positive values in the lower corona region to negative values at the wind base. These magnetic stresses are shown in Figure 7. In the corona region, both and stresses are mostly due to the mean field as shown by the blue curves.

Figures 6 and 7 show disk quantities along the disk height () at . The corona region is again marked as the yellow shaded area. We can see that the corona region has a supersonic inflow velocity, reaching 2 and transports a significant amount of mass inwards. The density reaches a plateau in this magnetically supported corona region. Figure 6 also suggests that the disk is sub-Keplerian in the corona region. Considering the corona region is magnetically dominated, the sub-Keplerian motion could be due to that magnetic fields in the corona connect to the midplane at larger radii. In other words, the midplane magnetically breaks the corona, as described in detail below. The panel in figure 7 also reveals that, on average, the midplane transports mass outwards from t=35 to 42 . The wind region which is beyond the corona region carries little mass outwards since the density is low there (more discussion on mass loss due to the wind is in §4.2).

To understand how magnetic stresses can lead to the coronal accretion, we analyze the angular momentum budget due to these stresses in different regions. Figure 4 suggests that the boundary between the corona and wind region is around z1.5 R, and the boundary between the midplane and corona region is around z0.1 R. Thus, the midplane, corona, and wind regions are better represented by wedges in the spherical-polar coordinate instead of slabs in cylindrical coordinates. Thus, we rewrite and analyze the angular momentum equation under the spherical-polar coordinate coordinate. After inserting the mass conservation equation and dividing the equation by , we get

 ∂⟨ρδvϕ⟩∂t =−1r3∂(r3⟨Trϕ⟩)∂r−⟨ρvr⟩r∂rvk∂r −1rsinθ2∂(sin2θ⟨Tθϕ⟩)∂θ−⟨ρvθ⟩rsinθ∂(sinθvk)∂θ (16)

where

 Trϕ≡ρvrδvϕ−BrBϕ Tθϕ≡ρvθδvϕ−BθBϕ

We then write Equation 16 as

 ∂⟨ρδvϕ⟩∂t=mrϕ+˙mr+mθϕ+˙mθ. (17)

The left term is the change of the azimuthal momentum. The first term on the right () is the radial gradient of the stress. At the disk midplane, this term is associated with disk turbulence. However, within the wind region, this term accelerates the wind since the magnetic fields are aligned with the radial direction. The second term on the right () is the angular momentum carried by the radial accretion flow. It is the radial accretion term. When it is positive, the disk accretes inwards. The third term () is the gradient of the stress. It can also be considered as the torque between different layers in the disk. The forth term () is the angular momentum loss due to the flow in the direction, such as disk wind.

These different terms are integrated over wedges as shown in Figure 8. The wedges are chosen in a way that they represent different disk regions at . At the disk midplane (upper left panel), the term is negative at R=1, suggesting that the flow is outwards. With our setup, we can see that the term balances the addition of two larger terms: and . At some other radii, the balance of two large terms can also lead to positive (or inward accretion). It is also possible that with some other disk parameters, the midplane could have different flow directions. As shown in Figure 3, follows at the disk midplane and it is almost a constant with height (Figure 9). Using the definition of , is which is negative. The negative term means that the stress tries to make the disk accrete inwards. However, the stress term () which is mainly due to net magnetic fields is positive, trying to make the disk flow outwards. The positive value of can be understood from Figure 9, where becomes large moving away from the midplane due to the azimuthal field stretch. Then, the term () is positive, driving disk outwards. We could also regard the term as the magnetic breaking. It torques the midplane to move outwards but the surface inwards. Eventually, at the disk midplane, the stress term wins over the stress term at some radii and the disk flows outwards at R=1.

In the corona region, both terms due to and stresses are negative, driving the disk to accrete inwards. We note that the term that is vertically integrated with a weight of sin is basically the difference between the stresses at the upper and bottom surfaces of the corona region. Since is negative and positive at the upper and bottom surfaces of the corona region, the integrated term is negative. Thus, the stress torques the disk surface to flow inwards and the disk midplane to flow outwards, again like magnetic breaking. At R=1, the coronal accretion rate is 10 times larger than the midplane outflow rate.

For the disk as a whole, including both the midplane and the corona, the large stress within the disk cannot lead to the overall disk accretion. Only the stress exerted at the upper and lower surface of the disk can torque the disk, leading to accretion. The lower right panel in Figure 8 shows that the total term is negative due to the stress exerted by the wind to the disk surfaces. Thus, the wind torque contributes to the disk accretion. However, the stress term is smaller than the stress term by a factor of 20 by examining the lower right panel in more detail. Thus, 95% of the disk accretion is due to the stress. Such stress is from MRI turbulence at the midplane and the global magnetic fields in the corona (the T panel in Figure 7).

To understand how the coronal accretion is determined by different components of stresses, we fit simple curves for both and stresses and use Equation 16 to calculate the accretion rate at different heights. We fit with the constant of 0.00032 and with

 Tθϕ=⎧⎪⎨⎪⎩0.00025×(θ−0.77)ifθ<1.4−0.00096×(θ−π/2)if1.4<θ<1.740.00025×(θ−2.37)if1.74<θ

shown as the red curve in the middle panel of Figure 9. Then with the assumptions that the disk is in a steady state, term is negligible, and , we can derive the accretion rate as

 ρvr=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩−0.00073−0.0005−0.001×ctgθ×(θ−0.77)ifθ<1.4−0.00073+0.0019+0.0038×ctgθ×(θ−π/2)if1.4<θ<1.74−0.00073−0.0005−0.001×ctgθ×(θ−2.37)if1.74<θ

shown as the red curve in the right panel of Figure 9. We can clearly see that the stress torques the disk midplane to flow outwards and the disk surface to flow inwards. The integrated mass flux (with the weight of sin) at is thus -2. Since the disk region extends from =0.69 to =2.45, we can derive the integrated flux to be -0.00135.

If we only use the stress, we can derive as -0.00073 and the vertically integrated disk mass flux as -0.00073. With and plugged in, the mass flux due to stress is -0.00128. From this simple model, 95 of disk accretion is due to the stress and from stress exerted at the disk surface, which is consistent with the value by examining the components in the lower right panel of Figure 8.

Since both turbulence and net magnetic fields can contribute to stresses, we would like to know their relative importance. We separate the magnetic stress into ()() where is the net field that has been averaged over the direction. Thus the azimuthally averaged stress can be divided into the stress due to net fields and the stress due to turbulence . The stresses are shown in the upper panels of Figure 10 which reveals that the stress is dominated by the MRI turbulence at the disk midplane. Previous work has established that in MRI-resolved models is 1/4 of , and the ratio between the Maxwell stress and the magnetic pressure (Hawley et al., 1995). In our simulation, the MRI turbulent midplane also satisfies these relationships (Figure 11). Figure 10 also shows that the stress is dominated by the net fields in the corona region, and by the turbulence again at the transition between the corona and wind regions. Since the stress determines the radial inflow, we conclude that both turbulence and net fields contribute to the radial accretion. On the other hand, the net fields dominate the stress in the wind and corona region and contributes significantly to the stress even at the disk midplane. Since the stress determines the internal flow structure within the disk, we conclude that the net fields determine the vertically-sheared motion.

Finally, we would like to know the total mass flux in different regions. Thus, we cut three wedges in the simulation domain and measure the radial mass fluxes, shown in Figure 12. We can see that the coronal accretion dominates. At , the midplane outflow rate is of the coronal accretion rate. However, at some other radii, the disk flows inwards at the midplane. In any case, the net flow rate at the midplane is always much smaller than the inflow rate in the corona region. The wind region has a very small mass flow rate. At larger distances, our cut for the wind region includes part of the disk region so that the mass flow rate does not reflect the true wind loss rate. A more detailed analysis on the outflow properties is presented in the next subsection.

### 4.2. Wind Region

Although the torque by the wind at the disk surface only accounts for 5% of the disk accretion rate, the disk wind is launched at all the radii in the disk and it can be directly probed by observations (Bjerkeli et al., 2016). Thus, we will study its properties in more detail.

To properly study the disk outflow/wind, the simulation domain needs to be larger than the Alfven surface and the fast magnetosonic surface (Blandford & Payne, 1982). This is satisfied in our simulation as shown in Figure 4 and 5 where the dashed curves are the Alfven surface while the dotted curves are the fast magnetosonic surface. Alfven surface is where the velocity in the poloidal plane equals the poloidal Alfven velocity . Fast magnetosonic surface is where equals the fast magnetosonic speed where . In the disk’s atmosphere is much larger than , so that .

For a steady axisymmetric outflow, there are four conserved quantities along the magnetic and velocity field lines (Weber & Davis, 1967; Blandford & Payne, 1982). In a steady flow, the induction equation becomes =0. If and are separated into the poloidal and the toroidal components ( and ), it has been shown that and are in the same direction, and using mass conservation equation we have the first constant

 k=ρvpBp, (18)

which is the mass load parameter, affecting the dynamical properties of the wind (Ouyed & Pudritz, 1999). Using the same equations, the second constant can also be derived

 ω=Ω−kBϕρR. (19)

With the angular momentum equation, we have the third constant

 l=R(vϕ−Bϕk) (20)

which is the specific angular momentum of the wind. In a barotropic fluid, we can use Bernoulli’s equation to derive the fourth constant

 e=12v2+Φ+h+B2ρ−B⋅vk (21)

where is . We can also write in some other ways using the first three constants above

 e =12v2+Φ+h+BϕBϕρ−Bϕvϕk (22) =12v2+Φ+h−RBϕωk, (23)

Although our disk is locally isothermal which is not barotropic, we will see that is much smaller than other terms (which means the wind is “cold” ) and we can still treat as a constant.

Using the azimuthally averaged velocity structure, we derive velocity streamlines and the four conserved quantities along the streamlines. We find that the conserved quantities are not constants. For example, the values of and have a factor of 2 peak at at t=42 . When we check these quantities over time, we notice that the peak in originates from the disk surface and propagates outwards. This is not surprising since the outflow is launched from the turbulent disk. Previous simulations have also observed the episodic disk wind (Ouyed & Pudritz, 1997b).

Thus, we average all primitive quantities over time, trying to smooth the unsteady wind and study the statistical properties of the wind. The averaged disk structure and the conserved quantities are shown in Figure 13. The conserved quantities are derived using the azimuthally and time averaged primitive quantities. The constant is calculated using the component of and . The conserved quantities are almost constants. This is encouraging, suggesting that we may still be able to use the traditional steady wind solution to study the statistical properties of the wind generated from the turbulent disk. Accordingly, all the relationships based on these four conserved quantities are also satisfied in our simulations. For example, is satisfied. We can verify this by noting that , (derived from in the lower right panel of Figure 13 and =0.35 at r=10 in Figure 14, or we can directly read the value of by following the streamlines in Figure 4 ), and in Figure 13. Thus, the wind lever arm is 10 at R=1.

Another useful relationship that should also be satisfied in our simulation is (Anderson et al., 2003), which reduces to

 J=12v2+Φ+h−ωRvϕ. (24)

At the wind base, , and at the distance far away from the base . Thus

 −3/2v2K∼12v2p,∞−ωR∞vϕ,∞, (25)

Since is much larger than and both and are , this can be reduced to (Anderson et al., 2003)

 v2p,∞2R∞vϕ,∞=√GM∗R30 (26)

which relates the wind properties at large distances to the launching point. This relationship is particularly useful in observations to estimate the position of the launching point by using quantities far away from the disk (Bjerkeli et al., 2016). We can test this relationship from our simulations. All primitive quantities () along the field are shown in Figure 14. If we take as far enough away from the wind base, we have , , and (based on Figure 13). Thus, derived from Equation 26 is , similar to the real launching point .

Some other wind properties can also be observed in Figure 14. Clearly, the density and magnetic fields drop off quickly along the streamline. The poloidal field is stronger than the toroidal field at the disk surface. But after the Alfven surface at , the toroidal field is stronger than the poloidal field. At the wind launching point, 1/ calculated using the net field reaches 1000 and 240 for the component. Bai et al. (2016) pointed out that at the wind base is one important parameter, besides the global field geometry and wind temperature, to determine the wind properties. In our particular case, the wind is highly magnetized when it is launched. Beyond the launching point, both and decreases. But decreases faster than . At the Alfven point and are roughly equal and the wind is dominated by beyond the Alfven point. On the other hand, the poloidal velocity increases along the streamline while the azimuthal velocity is almost a constant. Eventually the poloidal velocity reaches the terminal velocity 4 (, Pudritz et al. 2007, and 3 for wind launching from =1 in our simulation). Knowing and , we can estimate the dimensionless parameter

 μ≡ρvpvϕB2p (27)

to see if the wind is “light” or “heavy”, where all the quantities are estimated at the wind base. Using our measured values, is suggesting that the wind is very light. Anderson et al. (2005) has shown that the wind is still strongly collimated with such small .

Since the wind that is launched from different disk positions has different properties, we apply the same approximation to streamlines that are launched from other disk positions as shown in Figure 16. We can see that the wind is launched high in the disk around , beyond the corona region. The launching points are also overplotted as crosses in Figure 15. We can see that the material starts to flow outwards beyond . But the region between and sometimes flow inwards as the corona and sometimes flows outwards as the wind. Only the region beyond always flows outwards as the wind so that the conserved quantities are constant beyond . The two conserved quantities and are almost the same for different streamlines (their values are shown in Figure 13), while and . at the wind base roughly follows . The angle between the velocity vector and the vertical direction is 0.36-0.45 (). Although this angle is smaller than 30 required for launching disk wind from the disk midplane, the wind in our simulation is actually launched from high above the disk atmosphere (). If we use the marginally stable equipotential surface

 ϕ(R,z)=−GMR0[12(RR0)2+R0(R2+z2)1/2], (28)

we can derive that the critical angle to launch the wind at is only 17. Thus, field lines in our simulations are tilted enough to launch disk winds.

We can also use quantities along different streamlines in Figure 16 to estimate the total mass loss rate from the wind. If the wind is in a steady state, the mass loss rate between two streamlines will be a constant at different . Thus 2 between two streamlines will be a constant along . Since the separation between two streamlines does not change dramatically along r, we expect that should be roughly a constant. As shown in the upper middle panel of Figure 16, only changes by a factor of 2 from the launching point (black crosses) to the domain boundary (red crosses). We also know that wind that is launched at larger disk radii has a larger opening angle, as shown in the lower right panel. Thus, we can integrate over the opening angle at the domain boundary to derive the total mass loss rate from two sides of the the disk region [, ],

 ˙Mloss=4πr2∫θroθrisinθρvrdθ. (29)

where and are the positions of the streamlines that are originated from and . If we take , at , and at , we can estimate the from two sides of the disk region at [0.5,5]. Considering the total disk accretion rate is , the wind loss rate from this disk region is only 0.4% of the disk accretion rate.

We can also use the mass flux at the disk surface to estimate the two sided mass loss rate using

 ˙Mloss=∫4πRρvzdR. (30)

The panel in Figure 16 suggests that . Thus, , which is consistent with the estimation above.

We would also like to know how much angular momentum is carried away by the wind and if the wind torque is consistent with our direct measurement in §4.1. At the wind base, the stress term in Equation 2 can be written as

 ⟨ρvzvϕ⟩−⟨BzBϕ⟩=ρvzlR=14π∂˙Mloss∂RR2AR2ω (31)

using the conserved constants 333Note that the Reynolds stress in the wind uses instead of . But the stress at the disk surface won’t change much if we calculate the Reynolds stress using since that the magnetic stress dominates. . Figure 16 suggests that is around at the wind base around . Figure 13 suggests that . Thus the total stress at the wind base is around , which is roughly consistent with the direct measurement in Figure 7. This confirms that the wind torque is very small compared with the radial stress and the wind torque contributes little () to the disk accretion.

### 4.3. How to Maintain Global Magnetic Fields?

An important question to address is what is the rate of inward transport of vertical magnetic flux, and how is it related to the mass accretion rate? To answer this question, we plot near the inner boundary and at in Figure 17. The total magnetic flux through the sphere at is . In previous 3-D MHD simulations which have not covered the polar region, magnetic fields can be lost at the boundary close to the pole. However, our simulations cover the full 4 sphere, magnetic fields cannot be lost at the poles. If no field is being accreted, the total flux should remain a constant. In Figure 17, we can see that, after the initial relaxation, at the inner boundary increases with time, which suggests that magnetic fields are accreted to the central star while mass is being accreted. Simulations with longer timescale are needed to see if the accumulation of flux will saturate at some point. We also observe that magnetic field is strongest at the equator, and not the pole. However, we caution that a density floor is applied in the polar regions which may affect migration of the field from the equator.

Another important question regarding net magnetic fields is how to maintain such a large scale field in an accreting disk. The balance between the field advection and diffusion will determine the global field strength (Okuzumi et al., 2014; Takeuchi & Okuzumi, 2014). As discussed in the introduction, if Pr1, large-scale fields will diffuse outwards faster than the inward advection and the disk quickly loses magnetic flux (Lubow et al., 1994).

To understand how magnetic fields are maintained, we can write down the induction equation using the magnetic vector potential,

 ∂A∂t=v×B−η×∇×B (32)

where is the resistivity due to turbulence.

The poloidal field is determined by , which is

 ∂Aϕ∂t=−vRBz+vzBR−η∂BRR∂z+η∂Bz∂R. (33)

When the field is steady, the advection of fields (the first two terms on the right side) is balanced by the field diffusion due to turbulence. These four terms on the right side of the equation are plotted in Figure 18. In the wind region, the poloidal components of the velocity and magnetic vectors are parallel to each other so that . These two advection terms are balanced by each other. To calculate the turbulent diffusion terms, we assume the turbulent resistivity is 0.6 (where ) below the wind region and zero in the wind region, shown as the blue curve in the lower left panel. In the corona region, the radial inflow carries fields inwards (). With our choice of turbulent resistivity, this inward motion is balanced by the turbulence diffusion (the black curves in the top panels are almost zero.). Diffusion due to is more important close to the wind region, while diffusion due to is more important close to the disk midplane.

On the other hand, we can derive directly using equation 33 by assuming that the magnetic field structure has reached a quasi-steady state. In other words, fields advection and diffusion occur at a much shorter timescale than the secular evolution of the fields. Thus, by setting =0, we have

 η=−vRBz+vzBR∂BRR∂z−∂Bz∂R. (34)

Using this equation, the derived is shown in the lower right panel of Figure 18. This profile is quite similar to the profile, consistent with local simulations that . The quasi-steady state with seems to be in contradiction to previous studies that the disk loses magnetic fields due to the turbulent diffusion if . However, the inflow velocity in our simulation (e.g. 0.2 at R=1) is faster than (e.g. 0.03 from Figure 18) based on the viscous theory. This is due to the internal stress which creates the vertical shear. Another reason for the quasi-steady state is that, since the disk extends to , the vertical diffusion of the field is on the order of instead of . Thus, the faster inflow and the slower diffusion allow the disk to maintain a global field even if .

To confirm that the disk has evolved to a quasi-steady state when the secular evolution of the global field is much slower than the diffusion and advection, we checked the field strength in the simulation. From t=200 to t=420, at changes by at most 0.002. Then, is , which is two orders of magnitude smaller than other terms.

Besides studying the field structure within the disk, we would like to know the secular evolution of the global fields. As Ogilvie & Livio (2001) and Okuzumi et al. (2014) point out, the poloidal field evolution is determined by the balance between the conductivity-weighted radial velocity and conductivity-weighted resistivity. If we define the effective viscosity as the product of the conductivity-weighted radial velocity and R, we can calculate the effective Prandtl number as the ratio between the effective viscosity and the conductivity-weighted resistivity 444This Pr is different from the definition in Takeuchi & Okuzumi (2014) by a factor of ., and it can be reduced to

 Preff=12∫R−Rvr(z)η(z)dz. (35)

We assume that the disk extends from z=-R to R. Using time and azimuthally averaged and assuming , we calculate . If we use the averaged in the lower right panel of Figure 18, we calculate , although this value has a large uncertainty since the derived is very uncertain. Nevertheless, the effective Prandtl number is on the order of unity.

### 4.4. Different Net Fields and H/r

To explore how our results depend on the imposed magnetic fields strength, we have carried out simulations with initial but keeping 555 We have also carried out a case with . But that case behaves very differently from our fiducial case. The accretion rate is so high that the disk quickly loses mass and becomes magnetically dominated. We leave the discussion for the strongly magnetized disk to a later paper.. The velocity and magnetic fields structure are shown in Figure 19. Velocity streamlines are shown in the left panel and the magnetic streamlines are shown in the right panel. The color map in the right panel shows . We can clearly see that coronal accretion is still present in this case. The magnetically dominated corona still extends to . Similar to our fiducial case, field lines are pinched at the coronal and weak winds are launched.

The radial disk structure is presented in Figure 20. The disk evolves much slower than our fiducial case since the accretion is less efficient with instead of -0.005 in our fiducial case. Similar to our fiducial case, the vertically integrated () is significantly larger than the midplane , suggesting that most accretion occurs at the disk surface. On the other hand, is one order of magnitude smaller than our fiducial case, implying that is proportional to the net field strength.

We also did the angular momentum budget analysis as the fiducial case. Similar to our fiducial case, most of the accretion occurs at the surface. Less than 2% of disk accretion is due to the wind torque. Figure 21 shows that the inflow is still supersonic in the corona region, and some material is transported outwards at the disk midplane. The stress at the disk midplane is smaller than the value in the fiducial case by a factor of 3, consistent with local shearing box simulations (Hawley et al., 1995) that is proportional to the initial . On the other hand, the stress at the wind base is smaller than the value in our fiducial case by a factor of 10. Thus, the disk wind seems to play a less important role in the disk threaded by a weaker field.

To see if the coronal accretion picture will hold for thin disks, we have also tried one case with and . As shown in Figure 22, the coronal accretion still dominates the disk accretion. The corona still extends to . On the other hand, the disk accretion rate is -0.002, which is similar to the weak field case but smaller than the fiducial case (Figure 23). This lower accretion rate is mainly due to a smaller in the disk since is similar to the value in the fiducial case. Figure 21 also suggests that the stress at the wind base is 6 times smaller than the fiducial case. Considering the total accretion rate is only 2-3 times smaller than the fiducial case, the wind seems to play a less important role in thinner disks.

## 5. Discussion

### 5.1. Meridian Circulation

How mass is transported in an accretion disk is important not only for understanding star and planet formation, but also for explaining components of primitive meteorites, or chondrites, in our solar system (Cassen, 1996). The refractory inclusions in chondrites are formed at 1400-1800 K (Grossman, 2010). Such high temperature environment exists at the inner disk within 1 AU. In order to explain their presence in chondrites and even in comets (e.g. Simon et al. 2008), outward mass transfer is needed. In viscous disks, mass can flow outwards at the disk midplane, so-called “meridian circulation” (e.g. Takeuchi & Lin 2002). Such outward mass transfer could explain the refractory inclusions in meteorites (Ciesla, 2007; Hughes & Armitage, 2010). However, this“meridian circulation” is not supported in MHD simulations with net toroidal magnetic fields (Fromang et al., 2011).

In this paper, a “meridian circulation” pattern is found in our MHD simulations with net vertical fields. However, the driving mechanism in our simulations is entirely different from the traditional “meridian circulation” in viscous disks.

Takeuchi & Lin (2002) have shown that, in viscous disks with the stress and , the radial velocity at the disk midplane is positive whenever 3p+2q+60. With the normal disk parameters of q=-1/2 and p=-2.25, 3p+2q+6 is -1.75 and the disk flows outwards at the midplane due to viscous stresses. At larger , becomes negative and the disk accretes inwards. Fromang et al. (2011) have shown that such meridian circulation is mostly due to stresses. Even if