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

TL;DR

This study presents comprehensive global ideal MHD simulations of thin accretion disks with net vertical magnetic fields, revealing the dominance of coronal accretion, flux transport dynamics, and weak disk winds, advancing understanding of magnetic influences in accretion processes.

Contribution

First global ideal MHD simulations of thin accretion disks with net vertical magnetic fields, demonstrating coronal accretion, flux transport balance, and weak disk winds, with implications for astrophysical accretion systems.

Findings

Inward accretion occurs mainly in magnetically dominated upper disk regions.

A steady-state magnetic field geometry is established through flux transport and turbulent diffusion.

Weak, episodic disk winds are launched but contribute minimally to angular momentum transport.

Abstract

We report new global ideal MHD simulations for thin accretion disks (with thermal scale height H/R=0.1 and 0.05) threaded by net vertical magnetic fields. Our computations span three orders of magnitude in radius, extend all the way to the pole, and are evolved for more than one viscous time over the inner decade in radius. Static mesh refinement is used to properly resolve MRI. We find that:(1) inward accretion occurs mostly in the upper magnetically dominated regions of the disk, similar to the predictions from some previous analytical work and the "coronal accretion" in previous GRMHD simulations. Rapid inflow in the upper layers combined with slow outflow at the midplane creates strong $R\phi$ and $z\phi$ stresses in the mean field; the vertically integrated $\alpha\sim 0.5-1$ when the initial field has $\beta_{0}=10^3$ at the midplane. (2) A quasi-static global field geometry is established in which flux transport by inflows at the surface is balanced by turbulent diffusion. The field is strongly pinched inwards at the surface. A steady-state advection-diffusion model, with turbulent magnetic Prandtl number of order unity, reproduces this geometry well. (3) Weak unsteady disk winds are launched at $z/R\sim1$ with the Alfven radius $R_{A}/R_{0}\sim3$. Although the wind is episodic, the time averaged properties are well described by steady wind theory. Wind is not efficient at transporting angular momentum. Even with $\beta_{0}=10^3$, only 5% of the angular momentum transport is driven by torque from the wind, and the wind mass flux from the inner decade of radius is only $\sim$ 0.4% of the mass accretion rate. With weaker fields or thinner disks, the wind contributes even less. (4) Most of the disk accretion is driven by the $R\phi$ stress from the MRI and global magnetic fields. Our simulations have many applications to astrophysical accretion disk systems.

Figures

Equations

$${\partial }_t{A}_{\varphi }=$$

$$\frac{\partial \int ⟨\rho v_{\varphi }⟩dz}{\partial t}$$

$$\displaystyle-\frac{1}{2\pi R^{2}}\frac{\partial Rv_{k}\dot{M}_{acc}}{\partial R}-\left(\langle\rho v_{z}v_{\phi}\rangle-\langle B_{z}B_{\phi}\rangle\right)\bigg{|}_{z_{min}}^{z_{max}}$$

$$\frac{\partial \int ⟨\rho \delta v_{\varphi }⟩dz}{\partial t}$$

$$\displaystyle-\frac{\dot{M}_{acc}}{2\pi R^{2}}\frac{\partial Rv_{k}}{\partial R}-\left(\langle\rho v_{z}\delta v_{\phi}\rangle-\langle B_{z}B_{\phi}\rangle\right)\bigg{|}_{z_{min}}^{z_{max}}$$

$${\alpha }_{R\varphi ,Rey}=⟨\rho v_R\delta v_{\varphi }⟩/⟨p⟩\mathrm{and}{\alpha }_{R\varphi ,Max}=-⟨B_R{B}_{\varphi }⟩/⟨p⟩$$

$${\alpha }_{\varphi z,Rey}=⟨\rho v_z\delta v_{\varphi }⟩/⟨p⟩\mathrm{and}{\alpha }_{\varphi z,Max}=-⟨B_z{B}_{\varphi }⟩/⟨p⟩.$$

$${\alpha }_{int}=\frac{\int T_{R\varphi }dz}{\Sigma c_s^2}.$$

$$\dot{M}_{acc}=-\frac{2\pi}{\partial Rv_{k}/\partial R}\left(\frac{\partial}{\partial R}\left(R^{2}\alpha_{R\phi,int}\Sigma c_{s}^{2}\right)-R^{2}\langle B_{z}B_{\phi}\rangle\bigg{|}_{z_{min}}^{z_{max}}\right)$$

$${\rho }_0(R,z=0)={\rho }_0(R_0,z=0){\left(\frac{R}{R_0}\right)}^p,$$

$$T(R,z)=T(R_0){\left(\frac{R}{R_0}\right)}^q.$$

$${\rho }_0(R,z)={\rho }_0(R,z=0)\exp \left[\frac{GM}{c_s^2}\left(\frac{1}{\sqrt{R^2+z^2}}-\frac{1}{R}\right)\right],$$

$$v_{\varphi }(R,z)=v_K{\left[(p+q){\left(\frac{c_s}{v_{\varphi ,K}}\right)}^2+1+q-\frac{qR}{\sqrt{R^2+z^2}}\right]}^{1/2},$$

$$\rho_{fl}=\left\{\begin{array}[]{l}\rho_{fl,0}\left(\frac{R}{R_{0}}\right)^{p}\left(\frac{1}{z^{2}}\right)\quad{\rm if}\quad R\geq r_{min}\\ \rho_{fl,0}\left(\frac{r_{min}}{R_{0}}\right)^{p}\left(\frac{1}{z^{2}}\right)\quad{\rm if}\quad R<r_{min}\;{\rm and}\;r>3r_{min}\\ \rho_{fl,0}\left(\frac{R}{R_{0}}\right)^{p}\left(\frac{1}{z^{2}}\right)\left(5-2\frac{r-r_{min}}{r_{min}}\right)\left(4\frac{r_{min}-R}{r_{min}}+1\right)\\ \quad\quad\quad\quad\quad\quad\quad{\rm if}\quad R<r_{min}\;{\rm and}\;r<3r_{min}\\ \end{array}\right.$$

$$A_{\phi}=\left\{\begin{array}[]{l}\frac{1}{2}\times R\times B_{0}\left(\frac{r_{min}}{R_{0}}\right)^{m}\quad{\rm if}\quad R\leq r_{min}\\ \\ \frac{B_{0}}{R_{0}^{m}}\frac{R^{m+1}}{m+2}+\frac{B_{0}r_{min}^{m+2}}{R_{0}^{m}R}(\frac{1}{2}-\frac{1}{m+2})\quad{\rm if}\quad R>r_{min}\end{array}\right.$$

$$B_{z}=\left\{\begin{array}[]{l}B_{0}\left(\frac{r_{min}}{R_{0}}\right)^{m}\quad{\rm if}\quad R\leq r_{min}\\ \\ B_{0}\left(\frac{R}{R_{0}}\right)^{m}\quad{\rm if}\quad R>r_{min}\,.\end{array}\right.$$

$$\lambda =2\pi \sqrt{16/15}|v_A|/\Omega =9.17{\beta }^{-1/2}H$$

$$Q_{\theta }={\lambda }_{\theta }/(r\Delta \theta )$$

$$Q_{\varphi }={\lambda }_{\varphi }/(r\sin \theta \Delta \varphi )$$

$$\frac{\partial ⟨\rho \delta v_{\varphi }⟩}{\partial t}$$

$$-\frac{1}{r\sin {\theta }^2}\frac{\partial ({\sin }^2\theta ⟨T_{\theta \varphi }⟩)}{\partial \theta }-\frac{⟨\rho v_{\theta }⟩}{r\sin \theta }\frac{\partial (\sin \theta v_k)}{\partial \theta }$$

$$T_{r\varphi }\equiv \rho v_r\delta v_{\varphi }-B_r{B}_{\varphi }$$

$$T_{\theta \varphi }\equiv \rho v_{\theta }\delta v_{\varphi }-B_{\theta }{B}_{\varphi }$$

$$\frac{\partial ⟨\rho \delta v_{\varphi }⟩}{\partial t}=m_{r\varphi }+\dot{m_r}+m_{\theta \varphi }+\dot{m_{\theta }}.$$

$$\displaystyle T_{\theta\phi}=\left\{\begin{array}[]{l}0.00025\times(\theta-0.77)\quad{\rm if}\quad\theta<1.4\\ -0.00096\times(\theta-\pi/2)\quad{\rm if}\quad 1.4<\theta<1.74\\ 0.00025\times(\theta-2.37)\quad{\rm if}\quad 1.74<\theta\\ \end{array}\right.$$

$$\displaystyle\rho v_{r}=\left\{\begin{array}[]{l}-0.00073-0.0005-\\ \quad 0.001\times{\rm ctg}\theta\times(\theta-0.77)\quad{\rm if}\quad\theta<1.4\\ -0.00073+0.0019+\\ \quad 0.0038\times{\rm ctg}\theta\times(\theta-\pi/2)\quad{\rm if}\quad 1.4<\theta<1.74\\ -0.00073-0.0005-\\ \quad 0.001\times{\rm ctg}\theta\times(\theta-2.37)\quad{\rm if}\quad 1.74<\theta\\ \end{array}\right.$$

$$k=\frac{\rho {\mathbf{v}}_{\mathbf{p}}}{{\mathbf{B}}_{\mathbf{p}}},$$

$$\omega =\Omega -\frac{k{B}_{\varphi }}{\rho R}.$$

$$l=R(v_{\varphi }-\frac{B_{\varphi }}{k})$$

$$e=\frac{1}{2}{v}^2+\Phi +h+\frac{B^2}{\rho }-\frac{\mathbf{B}\cdot \mathbf{v}}{k}$$

$$e$$

$$=\frac{1}{2}{v}^2+\Phi +h-\frac{R{B}_{\varphi }\omega }{k},$$

$$J=\frac{1}{2}{v}^2+\Phi +h-\omega R{v}_{\varphi }.$$

Full text

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 $R-\phi$ stress within the disk that transports angular momentum radially, while a disk wind can lead to a $z-\phi$ 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; Bai & Stone, 2013; Bai, 2017; Simon et al., 2013a, b) 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 60o 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 $z\in[-0.5R,0.5R]$, 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

[TABLE]

if only radial motion is considered. In Equation 1, $A_{\phi}$ is the azimuthal component of the vector potential. On the right hand side, the first term is advection by the accretion flow. Since $v_{R}\sim\nu/R$ in the viscous model, this term is roughly $\nu B_{z}/R$. The second term is the vertical diffusion of radial field and is roughly $\eta B_{R}/z$. If we assume the large scale field enters the disk at an angle of $\sim 45^{o}$, as in the disk wind model, we have $B_{R}/z\sim B_{z}/H$. The third term is radial diffusion, roughly $\eta B_{z}/R$. Thus, the first and second terms dominate, and they balance each other when $\nu/\eta\sim R/H$ or the turbulent magnetic Prandtl number (Pr) $\sim R/H$. However, local shearing box simulations suggest that Pr$\sim$1 (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 2-D treatment of the magnetic flux transport is desired. Lovelace et al. (2009); Guilet & Ogilvie (2012, 2013) have carried out such study of the magnetic flux transport in 2-D and found that the fast flux transport at the disk surface may indeed sustain a steady magnetic field configuration, potentially solving the too-fast field diffusion problem. Since the radial distribution of global magnetic fields also determines the collimation of disk winds (Ogilvie & Livio, 2001; 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

[TABLE]

or

[TABLE]

for the perturbed quantities ($\delta v_{\phi}=v_{\phi}-v_{k}$). The symbol $\langle\rangle$ denotes that the quantity has been averaged over the $\phi$ direction 111$v_{k}$ has been assumed to be constant along $z$. Without this assumption, there will be an additional term related to $\dot{M}_{loss}\partial v_{k}/\partial z$., and $\dot{M}_{acc}=2\pi R\int\rho v_{R}dz$. We refer to $\langle\rho v_{R}\delta v_{\phi}\rangle$ and $-\langle B_{R}B_{\phi}\rangle$ as the radial Reynolds and Maxwell stress respectively, to distinguish them from the vertical stresses $T_{\phi z}$. If we normalize the stresses with pressure, we can identify different contributions to the total stress $\alpha$ as

[TABLE]

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

[TABLE]

where $T_{R\phi}$ is the sum of both radial Reynolds and Maxwell stress. If we choose $v_{k}$ so that $\langle\rho\delta v_{\phi}\rangle$=0 and assume that the magnetic stress dominates at the disk surface, Equation 3 can be written as

[TABLE]

Equation 5 suggests that disk accretion is due to both $T_{R\phi}$ within the disk and $T_{\phi z}$ exerted at the disk surface. Normally, the $R-\phi$ stress is from disk turbulence, and the $z-\phi$ stress at the surface is from a magnetocentrifugal disk wind. If both the $R-\phi$ and $z-\phi$ stresses have similar values, the second term on the right will be larger than the first term by a factor of $R/z$ 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 ($r$, $\theta$, $\phi$) for our simulations, which should minimize the effects of the domain boundary on the wind properties (Ustyugova et al., 1999). 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 ($R$, $\phi$, $z$) to denote positions in cylindrical coordinates and ($r$, $\theta$, $\phi$) to denote positions in spherical polar coordinates. In both coordinate systems $\phi$ represents the azimuthal direction (the direction of disk rotation).

Our simulation domain extends from $r_{min}=0.1$ to $r_{max}=100$ with $\theta$ from 0 to $\pi$. $\phi$ extends from 0 to $2\pi$ except for the thin disk case for which $\theta$ is from 0 to $\pi/4$. We implement a special polar boundary in the $\theta$ 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.

The initial density profile at the disk midplane is

[TABLE]

while the temperature is assumed to be constant on cylinders

[TABLE]

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

[TABLE]

and

[TABLE]

where $c_{s}=\sqrt{p/\rho}$ is the isothermal sound speed, $v_{K}=\sqrt{GM_{*}/R}$, and $H=c_{s}/\Omega_{K}$. A local isothermal equation of state is assumed during the simulation 222To achieve this, we have used the adiabatic equation of state with $\gamma$=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 $R=\max(R,r_{min})$ in the above equations.

Considering the density at the disk’s atmosphere is orders of magnitudes smaller than the density at the disk midplane, we adopt a density floor which varies with position

[TABLE]

This density floor becomes small at the disk atmosphere and becomes large close to the midplane and the inner radial boundary. When $\rho_{fl}$ gets too small ( smaller than $10^{-9}$ for $\beta_{0}=1000$ cases or $3\times 10^{-9}$ for the $\beta_{0}=10^{4}$ case), we choose $10^{-9}$ or $3\times 10^{-9}$ as the density floor.

Initial magnetic field is assumed to be vertical. To maintain $\nabla\cdot\bf{B}=0$, we use the vector potential $\bf{A}$ to initialize the magnetic fields. We set

[TABLE]

where $m=(p+q)/2$. Thus, the vertical magnetic fields are

[TABLE]

With this setup, the plasma $\beta$ at the disk midplane beyond $r_{min}$ is a constant.

We choose $p=-2.25$ and $q=-1/2$ in our simulation so that the disk surface density $\Sigma\propto R^{-1}$. Then $\rho_{0}(R_{0},z=0)$ is 1. The time unit in the simulation is $T_{0}=2\pi/\Omega(R=1)$. The grid is uniformly spaced in ln($r$), $\theta$, $\phi$ with 136$\times$64$\times$128 grid cells in the domain of [ln(0.1), ln(100)]$\times$[[math], $\pi$ ]$\times$[0, 2$\pi$] at the root level. Thus, $\Delta r/r$ equals 0.052 at the root level. We use an open boundary condition in the radial direction (but mass inflow into the domain is not allowed), the polar boundary condition in the $\theta$ direction (Appendix), and periodic boundary conditions in the $\phi$ direction.

We have carried out three simulations, one with $(H/R)_{R=R_{0}}=0.1$ and initial plasma $\beta$=1000 at the disk midplane (the fiducial case), one with $(H/R)_{R=R_{0}}=0.1$ and $\beta=10^{4}$ (the weaker field case), and one with $(H/R)_{R=R_{0}}=0.05$ and $\beta=10^{3}$ (the thin disk case). The density floor parameter $\rho_{fl,0}$ is $10^{-6}$ for the fiducial case and $3\times 10^{-7}$ for the other cases. For the simulations with $(H/R)_{R=R_{0}}=0.1$, three levels of refinement have been adopted towards the disk midplane, and $\Delta r/r$ at the midplane is 0.0065. Since $H/R$=0.1 at $R=1$, 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, $\phi$ in the thin disk case only extends from 0 to $\pi/4$. The above choice of $H/R$ covers the normal thickness range of protoplanetary disks. The choice of the field strength is based on both numerical and analytical considerations. MRI is sufficiently resolved numerically with a stronger field ($\beta_{0}=1000$), while analytical work (Guilet & Ogilvie, 2014; Okuzumi et al., 2014) suggests a weaker net field threading protoplanetary disks ($\beta_{0}\sim 10^{4}-10^{7}$). We will mainly present results for our fiducial case ($(H/R)_{R=R_{0}}=0.1$ and $\beta_{0}=1000$ ) since the MRI is better resolved at the midplane and the simulation runs longer, but other cases will be discussed in §4.4.

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

[TABLE]

(Hawley et al. 1995). With $\beta_{0}$=1000, $\lambda$ is $\sim$ 0.29 H so that the most unstable wavelength is resolved by 4.5 grid cells at $R=1$.

To demonstrate that our fiducial case has reached the necessary resolution to capture MRI, we plot the azimuthally averaged quality factor ($\langle Q_{\theta}\rangle$ and $\langle Q_{\phi}\rangle$) for $B_{\theta}$ and $B_{\phi}$ 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),

[TABLE]

where $\lambda_{\theta}=2\pi\sqrt{16/15}|v_{A,\theta}|/\Omega_{K}$ and $\lambda_{\phi}=2\pi\sqrt{16/15}|v_{A,\phi}|/\Omega_{K}$. The quantities $v_{A,\theta}$ and $v_{A,\phi}$ are the Alfven velocity calculated using the $\theta$ and $\phi$ components of the magnetic field, and $\Omega_{K}$ is the Keplerian angular velocity at the disk midplane. To get converged results for MRI turbulence, Sorathia et al. (2012) have shown that if $Q_{\phi}\approx 10$, $Q_{z}$ needs to be $\gtrsim 10-15$, and if $Q_{\phi}\gtrsim 25$, while $Q_{z}$ can be smaller ($\sim$ 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. $Q_{\phi}$ is larger than 10 in the whole disk ($>$50 at $R=1$), and $Q_{\theta}$ is around 10 at R=1. Thus, our simulation should have adequate resolution to capture MRI turbulence except at the midplane of the innermost region ($R<0.3$).

4. Results

After running for t=45.6 $T_{0}$, which is equivalent to 1442 orbits at the inner boundary, our fiducial disk reaches a steady state within R$\sim$3, i.e. the inner factor of $\sim 30$ in radius, as evident in the $\dot{M}$ panel of Figure 3 which shows that the disk accretion rate is almost a constant within R$\sim$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. The outer disk is still evolving. As shown in the $\alpha_{total,mid}$ panel of Figure 3, the MRI in the outer disk saturates at later times.

Once a quasi-steady state is reached at t=42 $T_{0}$ , the surface density profile follows $\sim 0.12R^{-0.6}$, and the midplane $\alpha_{r\phi}$ profile follows $\sim 0.1R^{0.5}$. With these profiles and the given temperature profile, the midplane stress then follows $\propto R^{-1.85}$. The vertically integrated $\alpha_{r\phi}$ follows $\sim 0.5R^{-0.4}$, leading to the profile of $R^{-1.5}$ for the vertically integrated stress. With these values and Equation 5, we can calculate $\dot{M}$=-0.0038, which is a constant with radii and consistent with the direct $\dot{M}$ measurement in the upper right panel of Figure 3. We can see that the vertically integrated $\alpha$ is much larger than the midplane $\alpha$ by a factor of $\sim$ 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 $R\phi$ stress should follow $R^{-0.85}$ and the midplane $\alpha$ and the vertically integrated $\alpha$ should have the same slope. But the stress panel in Figure 3 suggests that the vertically integrated $R\phi$ stress is much steeper with $R^{-1.5}$, which is because the stress is larger at the inner disk’s atmosphere where the fields are pinched as shown in §5.1. Since $\alpha_{int}\sim$ 0.5 at R=1, 42 $T_{0}$ is already longer than the viscous timescale ($R^{2}\Omega/(\alpha_{int}c_{s}^{2})$) at R=1, explaining why the inner 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 $B_{\theta,mid}$ panel of Figure 3. The plasma $\beta$ calculated with the net vertical magnetic field increases from $10^{3}$ to $10^{5}$ 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 $B/\Sigma$ 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 1 and 4. Figure 1 suggests that the disk surface is highly filamentary, especially in the corona region, similar to Suriano et al. (2017). Thus, to study the statistical properties of the disk, we average the quantities over the azimuthal direction, as shown in Figure 4. The green curves in the figure are the velocity field lines and the white curves are where the plasma $\langle\beta\rangle$=1. The most striking feature is that the disk accretes through the disk atmosphere where the disk is magnetically dominated ($\beta<1$) 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”. Analytical works by (Guilet & Ogilvie, 2012, 2013) have also seen such surface accretion when the turbulent viscosity is considered in a constant H/R disk.

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 $R=1$ are shown in Figure 6. The yellow shaded region is the corona region which we define as the region that is magnetically dominated ($\beta<1$) and has a negative radial velocity ($v_{r}<0$). The corona region separates the disk into three regions, with the wind region above and the disk midplane region below. Figure 6 shows that, away from the midplane, $B_{R}$ 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 $z>0$, the negative $B_{R}$ close to the midplane is stretched to the positive $B_{\phi}$ and the positive $B_{R}$ at the upper corona region is stretched to the negative $B_{\phi}$ (also shown in the color map of Figure 5). Net $B_{z}$ 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, crucial for maintaining coronal accretion and launching disk wind. For example, since $B_{\phi}$ and $B_{R}$ have opposite signs as shown above, the $R\phi$ magnetic stress, -$B_{R}B_{\phi}$, is positive from the midplane to the wind region, and the radial gradient of this stress leads to overall disk accretion. The $z-\phi$ stress -$B_{z}B_{\phi}$ has similar shapes as -$B_{\phi}$ since $B_{z}$ is always positive. At $z>0$, the stress changes from positive values in the lower corona region to negative values at the wind base, leading to vertically sheared motion which will be discussed below. These magnetic stresses are shown in Figure 7. We note that the net field dominates in the coronal region. Both $R\phi$ and $\phi z$ stresses are mostly from the mean field as shown by the blue curves.

Figures 6 and 7 show disk quantities along the disk height ($z$) at $R=1$. The corona region is again marked as the yellow shaded area with the wind region above and the disk midplane below. We can see that the corona region has a supersonic inflow velocity, reaching $\sim$ 2 $c_{s}$ 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 $<\rho V_{r}>$ panel in figure 7 also reveals that, on average, the midplane region transports mass outwards from t=35 to 42 $T_{0}$. 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 z$\sim$1.5 R, and the boundary between the midplane and corona region is around z$\sim$0.1 R. Thus, the midplane, corona, and wind regions are better represented by $\theta$ wedges in the spherical-polar coordinate instead of $z$ slabs in cylindrical coordinates. Thus, we rewrite and analyze the angular momentum equation under the spherical-polar coordinate. After inserting the mass conservation equation and dividing the equation by $r$, we get

[TABLE]

where

[TABLE]

We then write Equation 16 as

[TABLE]

The left term is the change of the azimuthal momentum. The first term on the right ($m_{r\phi}$) is the radial gradient of the $r-\phi$ 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 ($\dot{m_{r}}$) 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 ($m_{\theta\phi}$) is the $\theta$ gradient of the $\theta-\phi$ stress. It can also be considered as the torque between different layers in the disk. The forth term ($\dot{m_{\theta}}$) is the angular momentum loss due to the flow in the $\theta$ direction, such as disk wind.

These different terms that are integrated over $\theta$ wedges are shown in Figure 8. The wedges are chosen in a way that they represent different disk regions at $R=1$. Due to the symmetry of the problem, we will only focus on regions at $z>0$ or $\theta<\pi/2$. We can see that $m_{r\phi}$ ($\sim-1/r^{3}\partial(r^{3}T_{r\phi})/\partial r$) is negative in all these regions, thus trying to drive the disk inwards. The radial profile of $T_{r\phi}$ is shown in Figure 3, and it follows $r^{-1.85}$ at the disk midplane and it is almost a constant with height (Figure 9). With the definition of $m_{r\phi}$ in Equation 17, $m_{r\phi}$ is thus $-1.15T_{r\phi}/r$ which is negative.

At the disk midplane (upper left panel), the $\dot{m_{r}}$ term is negative at R=1, suggesting that the flow is outwards. With our setup, we can see that the $\dot{m_{r}}$ term balances the addition of two larger terms: $m_{r\phi}$ and $m_{\theta\phi}$. Although the $m_{r\phi}$ term is negative as shown above, the $\theta-\phi$ stress term ($m_{\theta\phi}\sim-\partial T_{\theta\phi}/\partial\theta$), which is mainly from net magnetic fields, is positive, trying to make the midplane flow outwards. The positive value of $m_{\theta\phi}$ can be understood from Figure 9, where $T_{\theta\phi}$ becomes large moving away from the midplane due to the azimuthal field stretch. We could also regard the $m_{\theta\phi}$ term as the magnetic breaking to the surface by the midplane. It tries to torque the surface to flow inwards and adds the momentum to the midplane driving midplane outwards. Eventually, at the disk midplane, the $\theta-\phi$ stress term wins over the $r-\phi$ stress term and the disk flows outwards at $R\sim 1$. We note that the balance of two large terms not necessarily leads to the inward accretion. At some other radii shown in Figure 8 or possibly in disks with some other parameters, the midplane can flow inwards.

In the corona region, both terms due to $r-\phi$ and $\theta-\phi$ stresses are negative, driving the disk to accrete inwards. We note that the $m_{\theta\phi}$ term that is vertically integrated with a weight of sin$\theta^{2}$ is basically the difference between the $\theta-\phi$ stresses at the upper and bottom surfaces of the corona region. Since $T_{\theta\phi}$ is negative and positive at the upper and bottom surfaces of the corona region, the integrated $m_{\theta\phi}$ term is negative. Thus, the $\theta-\phi$ 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 $\sim$10 times larger than the midplane outflow rate.

Although it is interesting to observe the different flow patterns at different vertical heights, we are more interested in the disk’s overall accretion rate as a whole. For the whole disk, including both the midplane and the corona, the large $\theta\phi$ stress within the disk cannot lead to the overall disk accretion. Only the $\theta\phi$ 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 $m_{\theta\phi}$ term is negative. Thus, the wind torque indeed contributes to the disk accretion. However, the wind torque (the $\theta-\phi$ stress term) is smaller than the $r-\phi$ 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 $r-\phi$ stress. Such stress is from MRI turbulence at the midplane and the global magnetic fields in the corona (the TR,ϕ panel in Figure 7).

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

[TABLE]

shown as the red curve in the middle panel of Figure 9. Then with the assumptions that the disk is in a steady state, $\dot{m_{\theta}}$ term is negligible, and $T_{r\phi}\propto r^{-1.85}$, we can derive the accretion rate as

[TABLE]

shown as the red curve in the right panel of Figure 9. We can clearly see that the $T_{\theta\phi}$ stress torques the disk midplane to flow outwards and the disk surface to flow inwards. The integrated mass flux (with the weight of sin$\theta^{2}$) at $r=1$ is thus $\sim-0.00073(\theta_{max}-\theta_{min})$-2$(T_{\theta\phi}(\theta_{max})-T_{\theta\phi}(\theta_{min}))$. Since the disk region extends from $\theta_{min}$=0.69 to $\theta_{max}$=2.45, we can derive the integrated flux to be -0.00135.

If we only use the $T_{r\phi}$ stress, we can derive $\rho v_{r}$ as -0.00073 and the vertically integrated disk mass flux as $\sim$ -0.00073$(\theta_{max}-\theta_{min})$. With $\theta_{max}$ and $\theta_{min}$ plugged in, the mass flux due to $T_{r\phi}$ stress is -0.00128. From this simple model, 95$\%$ of disk accretion is due to the $T_{r\phi}$ stress and $5\%$ from $T_{\theta\phi}$ 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 $B_{i}B_{j}$ into ($\langle B_{i}\rangle_{\phi}+\delta B_{i}$)($\langle B_{j}\rangle_{\phi}+\delta B_{j}$) where $\langle B\rangle_{\phi}$ is the net field that has been averaged over the $\phi$ direction. Thus the azimuthally averaged stress $\langle B_{i}B_{j}\rangle_{\phi}$ can be divided into the stress due to net fields $\langle B_{i}\rangle_{\phi}\langle B_{j}\rangle_{\phi}$ and the stress due to turbulence $\langle\delta B_{i}\delta B_{j}\rangle_{\phi}$. The stresses are shown in the upper panels of Figure 10. For the $r-\phi$ stress, we can see that the turbulent stress dominates at the midplane and the transition between the corona and wind regions, while the net field stress dominates in the corona and wind regions. Previous work has established that, in MRI turbulent models, $\alpha_{Rey}$ is $\sim$1/4 of $\alpha_{Max}$, and the ratio between the Maxwell stress and the magnetic pressure $\alpha_{mag}=T_{Max}/P_{b}\sim 0.45$ (Hawley et al., 1995). In our simulation, the MRI turbulent midplane also satisfies these relationships (Figure 11). Since the $r-\phi$ stress determines the radial inflow, the importance of both turbulence and net fields at different layers leads us to conclude that both processes contribute to the radial accretion.

On the other hand, the $\theta-\phi$ stress is dominated by net fields in both the wind and corona region and it is as important as the turbulent stress at the disk midplane. Since the $\theta-\phi$ 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 $R=1$, the midplane outflow rate is $10\%$ of the coronal accretion rate. However, we also note that, 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 $V_{p}\equiv\sqrt{\langle v_{r}\rangle^{2}+\langle v_{\theta}\rangle^{2}}$ equals the poloidal Alfven velocity $V_{Ap}\equiv\sqrt{(\langle B_{r}\rangle^{2}+\langle B_{\theta}\rangle^{2})/\langle\rho\rangle}$. Fast magnetosonic surface is where $V_{p}$ equals the fast magnetosonic speed $V_{F}\equiv\sqrt{1/2({{V_{A}}^{2}+c_{s}}^{2})+1/2|V_{A}^{2}-c_{s}^{2}|}$ where $V_{A}\equiv\sqrt{(\langle B_{r}\rangle^{2}+\langle B_{\theta}\rangle^{2}+\langle B_{\phi}\rangle^{2})/\langle\rho\rangle}$. In the disk’s atmosphere $V_{A}$ is much larger than $c_{s}$, so that $V_{F}\sim V_{A}$.

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 $\nabla\times(\bf{v}\times\bf{B})$=0. If $\bf{v}$ and $\bf{B}$ are separated into the poloidal and the toroidal components ($\mathbf{v}=\mathbf{v_{p}}+\Omega R\hat{\phi}$ and $\bf{B}=\bf{B_{p}}+\bf{B_{\phi}}$), it has been shown that $\bf{v_{p}}$ and $\bf{B_{p}}$ are in the same direction, and using mass conservation equation we have the first constant

[TABLE]

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

[TABLE]

With the angular momentum equation, we have the third constant

[TABLE]

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

[TABLE]

where $h$ is $\int dp/\rho$. We can also write $e$ in some other ways using the first three constants above

[TABLE]

Although our disk is locally isothermal which is not barotropic, we will see that $h$ is much smaller than other terms (which means the wind is “cold” ) and we can still treat $e$ 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 $k$ and $\omega$ have a factor of 2 peak at $r\sim 20$ at t=42 $T_{0}$. When we check these quantities over time, we notice that the peak in $k$ 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 $k$ constant is calculated using the $r$ component of $\bf{v_{p}}$ and $\bf{B_{p}}$. 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, $l=R_{A}^{2}\omega$ is satisfied. We can verify this by noting that $l\sim 14$, $R_{A}\sim 3$ (derived from $r_{A}=10$ in the lower right panel of Figure 13 and $\theta$=0.35 at r=10 in Figure 14, or we can directly read the value of $R_{A}$ by following the streamlines in Figure 4 ), and $\omega\sim 1.5$ in Figure 13. Thus, the wind lever arm $\lambda\equiv(R_{A}/R_{0})^{2}$ is $\sim$ 10 for field lines launched at $R_{0}=1$.

Another useful relationship that should also be satisfied in our simulation is $J=e-\omega l$ (Anderson et al., 2003), which reduces to

[TABLE]

At the wind base, $J\sim-3/2v_{K}^{2}$, and at the distance far away from the base $v_{p,\infty}>>v_{\phi,\infty}$. Thus

[TABLE]

Since $v_{p}$ is much larger than $v_{K}$ and both $\Omega$ and $\omega$ are $\sim\Omega_{K}$, this can be reduced to (Anderson et al., 2003)

[TABLE]

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 ($\rho,v_{p},v_{\phi},B_{\phi}$) along the field are shown in Figure 14. If we take $r=30$ as far enough away from the wind base, we have $v_{p}=4$, $v_{\phi}=1.2$, and $R=6$ (based on Figure 13). Thus, $R_{0}$ derived from Equation 33 is $\sim 1.1$, similar to the real launching point $R_{0}=1$.

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 $r\sim 10$, the toroidal field is stronger than the poloidal field. At the wind launching point, 1/$\beta$ calculated using the net field reaches 1000 and 240 for the $B_{p}$ component. Bai et al. (2016) pointed out that $\beta$ 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 $B_{p}$ and $B_{\phi}$ decreases. But $B_{p}$ decreases faster than $B_{\phi}$. At the Alfven point $B_{p}$ and $B_{\phi}$ are roughly equal and the wind is dominated by $B_{\phi}$ 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 ($\sim 2^{1/2}\Omega_{0}R_{A}$, Pudritz et al. 2007 with $R_{A}\sim$3 for wind launching from $R_{0}$=1 in our simulation). Knowing $v$ and $B$, we can estimate the dimensionless $\mu$ parameter

[TABLE]

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

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 $z\sim 1.5R$, 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 $z\sim R$. But the region between $z\sim R$ and $z\sim 1.5R$ sometimes flow inwards as the corona and sometimes flows outwards as the wind. Only the region beyond $z\sim 1.5R$ always flows outwards as the wind so that the conserved quantities are constant beyond $z\sim 1.5R$. We find that the two conserved quantities $l$ and $k$ for other streamlines are almost the same as the streamlines originated from $R=1$, while $\omega$ scales roughly $\sim 1.5\Omega_{K}$ and $e$ is roughly $\sim 10\omega$ as shown in Figure 16. $B_{p}$ at the wind base roughly follows $R^{-1.5}$. The angle between the velocity vector and the vertical direction is 0.36-0.45 ($20^{o}-25^{o}$). Although this angle is smaller than 30o required for launching disk wind from the disk midplane, the wind in our simulation is actually launched from high above the disk atmosphere ($z\sim 1.5R$). If we use the marginally stable equipotential surface

[TABLE]

we can derive that the critical angle to launch the wind at $z=1.5R_{0}$ is only 17o. Thus, field lines in our simulations are tilted enough to launch disk winds from the disk surface.

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 $r$. Thus 2$\pi r^{2}\int{\rm sin}\theta\rho v_{r}{\rm d}\theta$ between two streamlines will be a constant along $r$. Since the $\theta$ separation between two streamlines does not change dramatically along r, we expect that $r^{2}{\rm sin}\theta\rho v_{r}$ should be roughly a constant. As shown in the upper middle panel of Figure 16, $r^{2}{\rm sin}\theta\rho v_{r}$ 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 $r^{2}{\rm sin}\theta\rho v_{r}$ over the opening angle at the domain boundary to derive the total mass loss rate from two sides of the the disk region [$R_{in}$, $R_{out}$],

[TABLE]

where $\theta_{ri}$ and $\theta_{ro}$ are the $\theta$ positions of the streamlines that are originated from $R_{in}$ and $R_{out}$. If we take $r^{2}{\rm sin}\theta\rho v_{r}\sim 10^{-5}$, $\theta_{ri}=0.05$ at $R_{in}=0.5$, and $\theta_{ro}=0.2$ at $R_{out}=5$, we can estimate the $\dot{M}_{loss}\sim 2\times 10^{-5}$ from two sides of the disk region at [0.5,5]. Considering the total disk accretion rate is $5\times 10^{-3}$, 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

[TABLE]

The $\rho v_{z}$ panel in Figure 16 suggests that $\rho v_{z}\sim 2\times 10^{-6}r^{-1.5}$. Thus, $\dot{M_{out}}\sim 2\times 10^{-5}$, which is consistent with the estimation above.

We would also like to know how much angular momentum is carried away by the wind if we derive this wind torque using the measured wind properties, and if this calculated torque is consistent with our direct measurement in §4.1. At the wind base, the $\phi z$ stress term in Equation 2 can be written as

[TABLE]

using the conserved constants 333Note that the Reynolds stress in the wind uses $v_{\phi}$ instead of $\delta v_{\phi}$. But the stress at the disk surface won’t change much if we calculate the Reynolds stress using $\delta v_{\phi}$ since that the magnetic stress dominates. . Figure 16 suggests that $\rho v_{z}$ is around $2\times 10^{-6}$ at the wind base around $R=1$. Figure 13 suggests that $l\sim 10$. Thus the total stress at the wind base is around $2\times 10^{-5}$, 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 ($\sim 5\%$) to the disk accretion.

4.3. How to Maintain Global Magnetic Fields?

Since both MRI and disk wind depend on the strength of net magnetic fields, an important question to understand disk accetion is what is the rate of inward transport of net vertical magnetic flux. To answer this question, we plot $B_{r}$ both near and far away from the inner boundary at $r=0.12$ and $r=10$ in Figure 17. The total magnetic flux through the sphere at $r$ is $\int 2\pi r^{2}{\rm sin}\theta\langle B_{r}\rangle{\rm d}\theta$. In previous 3-D MHD simulations which have not covered the polar region, magnetic fields can be lost at the $\theta$ boundary close to the pole. However, our simulations cover the full 4$\pi$ 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, $B_{r}$ 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 Pr$\sim$1, 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 in our simulations, we can write down the induction equation using the magnetic vector potential,

[TABLE]

where $\eta_{turb}$ is the resistivity due to turbulence.

The poloidal field is determined by $A_{\phi}$, which is

[TABLE]

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 ${\bf v_{p}}\times{\bf B_{p}}=0$. The two advection terms ($v_{R}B_{z}$ and $v_{z}B_{R}$) are balanced by each other. To calculate the turbulent diffusion terms, we assume the turbulent resistivity is $\eta_{turb}=$0.6 $\nu$ (where $\nu=T_{total,R\phi}/(1.5\rho\Omega)$) below the wind region and zero in the wind region, shown as the blue curve in the lower left panel. In the upper corona region ($|z/R|>1$), $v_{p}$ and $B_{p}$ are still parallel to each other so that the two advection terms are balanced by each other. But at the lower corona region ($|z/R|<1$), the radial inflow carries fields inwards ($-v_{R}B_{z}<0$). With our choice of turbulent resistivity, this inward motion is balanced by the turbulence diffusion, especially the $dB_{R}/Rdz$ term there (the black curves in the top panels are almost zero.). Thus, the turbulent diffusion seems to balance the advection very well with the assumption that $\eta_{turb}/\nu$ is on the order of unity.

On the other hand, we can derive $\eta_{turb}$ directly using equation 40 if we assume that the magnetic field structure has reached a quasi-steady state, which is evident in our simulations. In other words, fields advection and diffusion occur at a much shorter timescale than the secular evolution of the fields. Thus, by setting $\partial_{t}A_{\phi}$=0, we have

[TABLE]

Using this equation, the derived $\eta_{turb}$ is shown in the lower right panel of Figure 18. This profile is quite similar to the $\nu$ profile, consistent with local simulations that $\eta_{turb}\sim\nu$. The quasi-steady state with $\eta_{turb}\sim\nu$ seems to be in contradiction to previous studies that the disk loses magnetic fields due to the turbulent diffusion if $\eta_{turb}\sim\nu$. However, two factors affect this. First, the inflow velocity in our simulation (e.g. 0.2 at R=1) is faster than $v_{R}\sim\nu/R$ (e.g. 0.03 from Figure 18) based on the viscous theory. This is due to the internal $\phi z$ stress which creates the vertical shear. Second, since the disk extends to $z\sim R$, the vertical diffusion of the field is on the order of $\eta B_{R}/R$ instead of $\eta B_{R}/H$. Thus, the faster inflow and the slower diffusion allow the disk to maintain a global field even if $\eta\sim\nu$. We want to caution that the curves in Figure 18 have large uncertainties (e.g. a factor of 2) since we need to average over 70 snapshots and also over the azimuthal direction to derive the mean fields out of the turbulent fields, and then we need to calculate the derivatives for these non-smooth curves to get the vector potential.

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, $B_{\theta}$ at $R=1$ changes by at most 0.002. Then, $\partial A_{\phi}/\partial t\sim B_{\theta}R/\Delta t$ is $\sim 10^{-5}$, 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 Preff is different from the definition in Takeuchi & Okuzumi (2014) by a factor of $C_{u}$., and it can be reduced to

[TABLE]

We assume that the disk extends from z=-R to R. Using time and azimuthally averaged $v_{r}$ and assuming $\eta(z)=T_{R\phi}/(1.5\rho\Omega)$, we calculate $Pr_{eff}\sim 2.7$. If we use the averaged $\eta(z)$ in the lower right panel of Figure 18, we calculate $Pr_{eff}\sim 10$, although this value has a large uncertainty since the derived $\eta(z)$ 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 $\beta_{0}=10^{4}$ but keeping $(H/R)_{R=R_{0}}=0.1$ 555 We have also carried out a case with $\beta_{0}=100$. 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 $B_{\phi}$. We can clearly see that coronal accretion is still present in this case. The magnetically dominated corona still extends to $z\sim 1.5R$. 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 $\dot{M}\sim-0.002$ instead of -0.005 in our fiducial case. Similar to our fiducial case, the vertically integrated $\alpha$ ($\alpha_{int}$) is significantly larger than the midplane $\alpha$, suggesting that most accretion occurs at the disk surface. On the other hand, $\alpha_{int}$ is one order of magnitude smaller than our fiducial case, implying that $\alpha_{int}$ 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 $T_{R\phi}$ stress at the disk midplane is smaller than the value in the fiducial case by a factor of $\sim$3, consistent with local shearing box simulations (Hawley et al., 1995) that $\alpha$ is proportional to the initial $v_{A}$. On the other hand, the $T_{\phi z}$ stress at the wind base is smaller than the value in our fiducial case by a factor of $\sim$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 $(H/R)_{R=R_{0}}=0.05$ and $\beta=1000$. As shown in Figure 22, the coronal accretion still dominates the disk accretion. The corona still extends to $z\sim R$. On the other hand, the disk accretion rate is $\sim$-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 stress associated with a weaker field. Even though the midplane $\beta_{0}$ is the same as the fiducial case, the 4 times smaller gas pressure means that the initial magnetic field is weaker by a factor 2, which leads to a weaker stress.

smaller $c_{s}$ in the disk considering $\alpha_{int}$ is similar to the value in the fiducial case. Figure 21 also suggests that the $\phi z$ stress at the wind base is $\sim$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 too.

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 $\sim$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” ( Urpin 1984; Takeuchi & Lin 2002; Jacquet 2013; Philippov & Rafikov 2017). Such outward mass transfer could explain the refractory inclusions in meteorites (Ciesla, 2007; Hughes & Armitage, 2010). However, this “meridian circulation” pattern with the midplane going out and surface going in is not supported in MHD simulations with net toroidal magnetic fields (Flock et al., 2011; Fromang et al., 2011). Fromang et al. (2011) found that disk material is going out at all disk heights.

In this paper, a “meridian circulation” pattern is found in our MHD simulations with net vertical fields ( a similar finding is reported in Suzuki & Inutsuka 2014 despite that their surface inflow is very close to their $\theta$ boundary). 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 $T_{R\phi}=\rho\nu R\partial\Omega/\partial R$ and $T_{\phi z}=\rho\nu R\partial\Omega/\partial z$, the radial velocity at the disk midplane is positive whenever 3p+2q+6$<$0. 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 $z$, $v_{R}$ becomes negative and the disk accretes inwards. Fromang et al. (2011) have shown that such meridian circulation in the viscous disk is due to the $R-\phi$ stress, while the $z-\phi$ stress actually tries to drive the midplane inwards. When $T_{\phi z}$ stress is completely ignored, the radial velocity is

[TABLE]

where $\nu=\alpha_{visc}c_{s}h$ 666$\alpha_{R\phi,visc}$ defined in this way is smaller than $\alpha_{R\phi}$ defined in our paper ($\alpha_{R\phi}=T_{R\phi}/\rho c_{s}^{2}$) by a factor of 1.5 . With q=-1/2 and p=-2.25, the disk still flows outwards at the midplane and flows inwards at the surface.

However, for the coronal accretion presented here, the vertically sheared motion is mostly due to the anomalous $z-\phi$ stress instead of the $T_{R,\phi}$ stress as in viscous disks. In most MHD simulations, the turbulent $T_{R,\phi}$ stress is almost uniform along the vertical direction instead of being proportional to the density (Fromang et al., 2011). In Figure 24, we can see that the stress is even higher at the atmosphere at our inner disk, and the stress is mostly due to the mean fields at the atmosphere (Figure 7). For an order-of-magnitude estimate, we assume that the stress is uniform vertically, and the $\alpha$ parameter at the disk midplane vary radially as $\alpha=\alpha_{0}(R/R_{0})^{\gamma}$. Then the disk accretion rate is

[TABLE]

based on Equation 25 in Fromang et al. (2011). In our fiducial case with $\gamma$=0.5, $p$=-1.85, $q$=-0.5, the whole disk should flow inwards. Our angular momentum analysis in Figure 8 confirms that the $r-\phi$ stress term is negative, trying to drive the disk to accrete inwards. However, Figure 8 also reveals that the $z-\phi$ (or $\theta-\phi$) stress transports the angular momentum from the disk surface to the midplane, or in other words the midplane magnetically breaks the surface, leading to the surface coronal accretion.

We have summarized how such $z-\phi$ stress term is produced and how it feeds back to the disk accretion flow in Figure 25. When the disk accretes at the surface, the magnetic fields at the surface are dragged inwards with the flow. This produces the negative $B_{R}$ component close to the midplane at $z>0$ (the lower left panel in Figure 6). Due to the Keplerian shear, the inward magnetic fields is twisted in the azimuthal direction, producing positive $B_{\phi}$ close to the midplane (the lower right panel in Figure 6). With the positive net $B_{z}$, the magnetic stress $-B_{\phi}B_{z}$ is negative at $z>0$ and $\partial T_{\phi z}/\partial z$ (or $\partial T_{\theta\phi}/\partial\theta$) is negative, driving the midplane to flow outwards (Equation 16). At the top of the corona region, $B_{R}$ is positive and it is twisted to generate negative $B_{\phi}$ due to the Keplerian shear. Thus $B_{\phi}$ decreases from positive to negative within the corona region, thus producing positive $\partial T_{\phi z}/\partial z$, which leads to the inward surface accretion. Such vertically sheared motion then feeds back to the field geometry. This picture is quite similar to the linear growth phase of MRI, but operating at the global scale.

To make the connection between the $\theta-\phi$ stress (Figure 10) and the classical viscous disk model, where

[TABLE]

using the relationship $\partial\Omega/\partial z\sim-\Omega_{K}z/R^{2}$ at the midplane from Equation 9, we can calculate $\alpha_{vis,\theta\phi}$ as $T_{\theta\phi}/(P(\theta-\pi/2))$ which is similar to $T_{\phi z}R/(Pz)$. Figure 10 shows that $T_{\theta\phi}$ is negative at z$>$0 and positive at z$<$0. Thus, $\alpha_{vis,\theta\phi}$ is negative, which is qualitatively different from the viscous model. In the viscous model, the positive $\alpha_{vis,\theta\phi}$ tries to drive the midplane flow inwards while $\alpha_{vis,r\phi}$ tries to drive the midplane flow outwards. In our MHD models, the negative $\alpha_{vis,\theta\phi}$, due to the mechanism illustrated in Figure 25, drives the midplane flow outwards. The magnitude of the $\alpha_{vis,\theta\phi}$ in our model is shown as the dashed curve in the lower right panel of Figure 10, which is comparable to $\alpha_{r\phi}$.

Overall, the vertically sheared flow is mostly determined by the internal $z-\phi$ (or $\theta-\phi$) stress from net magnetic fields. However, the $z-\phi$ stress contributes little to the disk overall accretion since this stress at the corona upper surface due to the wind torque is weak. As described in §4, the inward accretion of the disk is mostly due to the $r-\phi$ stress (illustrated in Figure 25). In the corona region, the $r-\phi$ stress is provided by the net fields, while at the disk midplane it is from MRI turbulence.

5.2. Accretion Mechanism and Wind vs Turbulence

If we are only interested in the net disk accretion rate, we can vertically integrate the angular momentum equation and it is the $z-\phi$ stress at the disk surface that matters (Equation 5). If we plug Equation 5 into the mass conservation equation

[TABLE]

and assume $v_{k}\propto R^{1/2}$, we have

[TABLE]

Inserting Equation 38, we get

[TABLE]

where $(R_{A}/R)^{2}\equiv\lambda$ is the wind lever arm. The terms on the right are accretion due to the radial stress gradient, the accretion due to the wind torque, and the mass loss rate due to the wind. Since $\omega\sim\Omega_{K}$, the accretion rate led by the wind torque is

[TABLE]

In our fiducial run, $\alpha_{R\phi,int}$ is 0.5 at $R=1$ and can be larger at the inner disk (Figure 3). The accretion rate due to the radial stress gradient is on the order of $-0.005$, the accretion due to the wind torque is $-2.5\times 10^{-4}$, and the mass loss from a large disk region [0.5,5] is $2\times 10^{-5}$. Thus, we can see that the importance of the terms decreases on their orders at the right side. Most disk accretion is led by the radial stress. Such stress is from the MRI turbulence at the disk midplane and the net global field at the disk corona. The wind carries the angular momentum away but it only leads to 5% of the accretion. The direct mass loss by the wind is another order of magnitude smaller.

We can also relate the accretion due to the wind and the wind mass loss rate by using Equation 49 and $\lambda\sim 10$ for the wind that launches at $R=1$. Considering $\partial\dot{M}_{loss}/\partial R=4\pi R\rho v_{z}$ and $\rho v_{z}\sim 2\times 10^{-6}$ (Figure 16), we can derive the mass accretion rate due to the wind is $\sim 3\times 10^{-4}$, consistent with our direct measurement.

Overall, the accretion mechanism in our simulations is very different from most previous works. The accretion is mostly driven by neither the turbulence nor the wind. Instead, it is mainly driven by the large scale net magnetic fields in the atmosphere, with some contribution from MRI turbulence at the midplane. The net field is sustained by the interplay between the accretion and diffusion. Large scale net magnetic fields have been found in some previous local MHD simulations. Turner & Sano (2008) has suggested that such fields can lead to accretion in the “undead zone”. We can also make an analogy between the field structure in our simulations with channel flows in magnetized disks (Goodman & Xu, 1994; Latter et al., 2010) where the stretch between different channel layers creates large scale magnetic fields which transport angular momentum.

5.3. Time and Spatial Variability

Due to MRI turbulence, both disk accretion and wind launching are stochastic. To asses the time variability of various quantities in our fiducial run, the space time diagram and the time evolution of various quantities are shown in Figures 26 and 27. We can see that the flow at the disk midplane is not always outwards. There seems to be a cyclic behavior on a timescale of 1.5 T0. Both $B_{r}$ and $B_{\phi}$ have similar cyclic patterns. However, different from zero field or net toroidal field simulations (Stone et al., 1996), the net $B_{\phi}$ field does not exhibit negative and positive alternation during one cycle. Instead $B_{\phi}$ is dominated by the stress from the coronal accretion and always positive at $z>0$ and negative at $z<0$. Previous shearing box simulations suggests that the cyclic dynamo behavior disappears when $\beta_{0}<1000$ (Bai & Stone, 2013) or $\beta_{0}<10$ (Salvesen et al., 2016). The $\beta$ calculated by the net $B_{z}$ in our simulation is slightly larger than 1000 during the quasi-steady state. However, our disk is also threaded by a strong net $B_{\phi}$ due to the azimuthal stretch of $B_{R}$ and it is unclear how the net $B_{\phi}$ will affect the cyclic dynamo.

5.4. Compared with Previous work

Our simulation indicates that the corona region extends to z$\sim$1.5 R and beyond which the magnetocentrifugal wind is launched. This suggests that many previous simulations do not have enough vertical range to capture the disk wind. Being unable to capture the wind region could explain why the mass loss rate in Fromang et al. (2013) is reduced by one order of magnitude when the shearing box extends from 5 to 10 H above the midplane (even $10$H is within our corona region). The tallest box simulation in Fromang et al. (2013) suggests that the disk depletion time is 700 orbits. Our disk patch at $R\sim$1 has a depletion timescale of $\Sigma/(\rho v_{z})\sim 0.7/2\times 10^{-6}\sim 5\times 10^{4}$ orbits, much longer than the depletion timescale in shearing box simulations. This implies that shearing box simulation is inadequate for studying disk wind.

The wind exerts a torque to the upper surface of the corona region. Since the corona region is magnetically dominated, it can accrete inwards as a whole. The corona region shares similarities with the sub-Keplerian disk region proposed by Wardle & Koenigl (1993) where the disk is confined by magnetic stresses.

The vertically sheared disk motion in our simulation has been observed in 3-D global MHD simulations by Suzuki & Inutsuka (2014). However, their simulations only extend to $z\sim 0.5R$, well within our corona region. Thus, it is unclear if the surface inflow they observed is from part of the corona, due to the wind as they suggested, or due to the imposed boundary. The coronal accretion in our simulation is very similar to the results in GRMHD simulation by Beckwith et al. (2009) and earlier simulations by Stone & Norman (1994), despite the very different numerical setups. Stone & Norman (1994) make an analogy with the channel mode of MRI. However, the processes on how the flow motion determines the global fields and how the global fields sustain the flow motion have not been studied in detail. In this paper, we have presented a careful analysis of the angular momentum budget and have studied the interplay between flow motion and global fields. We also point out that although the $z-\phi$ stress drives coronal accretion, it is the $R-\phi$ stress that determines the total accretion rate. Furthermore, although Beckwith et al. (2009) suggests that the parameterization of local turbulent viscosity and turbulent resistivity may not be sufficient to describe the field transport, our simulations agree quite well with R-z 2-D disk models (espeically Guilet & Ogilvie 2012, 2013) that still use local turbulent viscosity and resistivity but consider R-Z 2-D mean field and flow dynamics.

Although our simulations are qualitatively similar to the surface inflow picture suggested by these analytical models, the simulations do not support some assumptions made in some of these models. Rothstein & Lovelace (2008) and Lovelace et al. (2009) assume conducting non-turbulent surface layers, while the simulations imply a turbulent transition region between the coronal and the wind region. On the other hand, our simulations agree with Guilet & Ogilvie (2012) and Guilet & Ogilvie (2013) surprisingly well. Our $\alpha$ profile (Figure 10) is similar to that assumed by Guilet & Ogilvie (2013) (GO, the red curve in their Figure 1) and the resulting accretion structures are similar (the comparison is in Figure 28). For example, our maximum surface inflow velocity is 2 $c_{s}$ while it is $c_{s}$ in GO with H/R=0.1. The $v_{\phi}$ in our coronal region decreases to 5 $c_{s}$, while it decreases to 3 $c_{s}$ in GO. The maximum $B_{R}/B_{z,0}$ is 1 in our model and 2 in GO. The maximum $B_{\phi}/B_{z}$ is 10 in both models. Despite the surprising similarity, there are some notable differences. First, the maximum inflow occurs at 10 H in our models while it is 4 H in GO. We think such difference is due to that the coronal region is magnetically supported and the density there is different from the gaussian profile assumed by GO. Second, we have the wind region beyond the coronal region, which is not the case in GO due to the specific boundary condition assumed. Future analytical calculations with these considerations are desired.

5.5. Observables

Our simulations have observational implications for cataclysmic variables (CVs) and protoplanetary disks.

By fitting the light curves of CVs, previous studies have constrained that the $\alpha$ parameter is $\sim$0.1-0.3 during the outburst state (e.g. Mineshige & Osaki 1983; Smak 1984, 1999; Meyer & Meyer-Hofmeister 1984; Cannizzo 1993; Kotko & Lasota 2012). Since the disk should be fully ionized during the outburst state, MRI is thought to be the main angular momentum transport mechanism in such disks (Gammie & Menou, 1998). However, MRI in zero net field simulations can only generate $\alpha\sim 0.01$ which is much weaker than the $\alpha$ constrained from observations. There are two main solutions to this problem:1) MRI is enhanced by thermal convection (Hirose et al., 2014); 2) the disk is threaded by net vertical magnetic fields. The second solution has met several challenges (King et al., 2007): a relatively strong net vertical field (e.g. $\beta\sim 10$) is required to generate $\alpha\sim$0.1 based on local shearing box simulations, a strong disk wind will be generated if the disk is threaded by large fields, and the global field can be lost in the disk. On the other hand, our simulations suggest that we can generate vertically integrated $\alpha\sim 0.5$ even with a weak net field ($\beta\sim 10^{3}-10^{4}$). Part of this large stress is due to the geometry of the global magnetic fields. A very weak wind is launched in the disk and the net magnetic fields are maintained in our simulations. Thus, a weak net vertical field which could be carried by the inflow remains one solution to the large $\alpha$ in outbursting CV disks.

One puzzle for disks in Herbig stars is that their near-IR fluxes are several times larger than those predicted from radiation hydrodynamic models (Mulders & Dominik, 2012; Flock et al., 2016). Since most near-IR flux comes from the inner dust rim where dust sublimates, this discrepancy could be due to a high disk atmosphere that is magnetically supported (Turner et al., 2014a). Our simulations support this picture where the magnetic dominated corona extends to z$\sim$R. Figure 29 shows the position where the disk’s vertically integrated density is 0.001 and 0.01. For comparison, the two-sided surface density at $R=1$ is 0.1. As long as the $\tau=1$ surface is above the 0.01 curve, the corona will play an important role on the spectrum energy distribution. If we scale our fiducial simulation by choosing the length unit as 0.1 AU, the central star mass as 1 solar mass, and the disk accretion rate (-0.005 code unit) as $-10^{-8}\rm M_{\odot}\,yr^{-1}$, we can calculate the unit of the surface density as 8.9 g cm*-2*. Thus, the disk surface density at $R=1$ is 0.18.9$\sim$ 1 g cm-2*. If we assume that the Rosseland mean opacity is 10 $cm^{2}g^{-1}$, the disk’s optical depth is 10. Thus, the blue curves have optical depths of 0.1 and 1. If the disk has an accretion rate of $10^{-7}\rm M_{\odot}\,yr^{-1}$, the highest blue curve will correspond to the $\tau=1$ surface, which is twice higher than the $\tau=1$ surface in the hydrodynamical model. We note that our disk’s aspect ratio at R=1 is larger than the aspect ratio of a real disk at 0.1 AU which is normally $\sim$ 0.03. Although MHD simulations with realistic thermal dynamics is needed for addressing this problem properly in future, our simulations suggest that magnetically supported corona may play an important role in explaining the strong near-IR flux in Herbig Ae/Be stars.

Another puzzle our simulations may shed light on is the fast inflow in transitional disks. Transitional disks have moderate accretion rates but low surface density (Espaillat et al., 2014; van der Marel et al., 2016). This implies a large $\alpha$ in the disk or even supersonic inflow. Fast inflow could also explain the twist of channel maps observed in some transitional disks (Rosenfeld et al., 2014; Pineda et al., 2014; Casassus et al., 2015; van der Plas et al., 2016). Rosenfeld et al. (2014) derive that the inflow in HD 142527 approaches the infall velocity or the disk is warped. In our simulation, the inflow velocity in the corona region approaches 0.2 vK. On the other hand, ambipolar diffusion, which should operate in low density regions, may also lead to fast inflow (Wang & Goodman, 2016). MHD simulations with ambipolar diffusion included are desired in future to solve this puzzle.

Strong outflows have been observed in FU Orionis system. Since these disks are fully ionized within $\sim$AU (Zhu et al., 2007), we should be able to directly compare our simulations with observations. Recent ALMA high angular resolution observations (Zurlo et al., 2017) reveal a wide hourglass shape outflow with an outflow velocity of several km/s. Our outflow velocity is at a comparable rate. As shown in Figure 30, the terminal velocity of the outflow ranges from 0.5 to 3. If we assume that $R=1$ in our simulation corresponds to the ionized disk size ($\sim$1 AU) and the central star mass is 0.3 $M_{\odot}$ (Zhu et al., 2007), the terminal velocity in our simulations ranges from 8 km/s to 48 km/s. On the other hand, our simulations cannot explain some observables. For example, (Calvet et al., 1993) have estimated an outflow rate of 10${}^{-5}\rm M_{\odot}\,yr^{-1}$ for FU Ori which has an accretion rate of $2\times 10^{-4}\rm M_{\odot}\,yr^{-1}$ . In our simulations, the outflow rate from R=0.5 to 5 is only 0.4% of the accretion rate, while observations suggest that the outflow rate is 5% of the accretion rate. One solution is that the disk is threaded by a stronger net vertical magnetic fields. Strong magnetic fields have been observed in FU Orionis system (Donati et al., 2005). Another solution is that the wind strength will be much stronger in a thicker disk whose gravitational potential is smaller at the wind launching points so that more mass can escape the disk.

6. Conclusions

We have carried out global ideal MHD simulations to study accretion disks threaded by net vertical magnetic fields. Static mesh refinement has been adopted at the disk midplane to capture the growth of magnetorotational instability (MRI), and special boundary conditions have been used to prevent the loss of magnetic fields at the polar region.

For our fiducial case which has an initial field of $\beta=1000$ at the midplane, after running for 1442 orbits at the inner edge, the accretion flow reaches a steady state from R=0.1 to 3. The vertically integrated $\alpha$ follows $R^{-0.4}$, reaching almost 1 at the inner disk. Due to this $\alpha$ profile, the disk surface density follows $R^{-0.6}$ which is shallower than $R^{-1}$ in a viscous disk having a constant $\alpha$.

The disk exhibits a complicated accretion pattern with supersonic inflow at the corona region and little inflow or even outflow at the midplane. The corona region is magnetically dominated and the inflow velocity can reach $\sim 2c_{s}$. Although the inward surface flow is filamentary, it drags magnetic fields inwards, pinching fields at the disk surface. The Keplerian shear stretches the radial fields into azimuthal fields and creates large $z-\phi$ stress between the midplane and the corona. Such $z-\phi$ stress torques the corona inwards and the midplane outwards, thus sustaining the vertically sheared flow motion. If we treat the disk as a whole and neglect the internal sheared motion, only 5% of the disk accretion is due to the wind torque. 95% of the disk accretion is driven by the radial $r-\phi$ stress. Such stress is from MRI turbulence at the disk midplane and large scale net magnetic fields at the disk atmosphere. However, disk wind can play a more important role when non-ideal MHD effects have weakened MRI at the disk midplane (Bai, 2017; Béthune et al., 2017).

Even with such strong net vertical fields, a very weak disk wind is launched beyond the corona region at $z\sim 1.5R$. Although the wind is episodic, the time averaged wind properties can be fully described by the four conserved quantities from the steady wind theory. The launching angle (20o-25o) is smaller than 30o required to launch wind from the disk midplane, but it is enough to launch wind from $z\sim 1.5R$. The wind is highly magnetized and very weakly loaded. The mass loss rate from R=0.5 to 5 is only 0.4% of the disk accretion rate.

When a weaker net field has been applied ($\beta_{0}=10^{4}$) or a thinner disk has been considered, the wind torque accounts for less of the total accretion. Supersonic inflow at the disk surface also occurs and the picture of coronal accretion remains the same. The corona still extends to $z\sim 1.5R$. $\alpha_{int}$ is significantly smaller when a weaker field is applied. However, for a thinner disk which is also threaded with $\beta_{0}=10^{3}$ as in our fiducial case, $\alpha_{int}$ is similar to the value in the fiducial case.

Magnetic fields are accreted to the central star along with mass. The global magnetic field geometry in the disk is sustained by the inward accretion and turbulence. $Pr\sim$1 still seems to hold throughout the disk. Such steady field configuration is possible since 1) the surface inflow is faster than the viscous flow and 2) the magnetic structure in the disk, which extends to $z\sim R$, is much thicker than the thermal scale height so that magnetic fields diffuse much slower. Our simulations show excellent agreement with some previous analytical studies which consider magnetic flux transport in 2-D turbulent disks (e.g. Guilet & Ogilvie 2012 and Guilet & Ogilvie 2013).

Our simulations may shed light on some astrophysical problems. The large $\alpha$ in our simulations may be applicable to cataclysmic variables. The puffed corona may help to explain the high near-infrared flux in Herbig stars. Fast surface inflow may be the cause for the fast inflow in transitional disks. The self-consistent accretion and outflow geometry may be applied to FU Orionis systems. Finally, the vertically sheared mass transport may play an important role on transporting chondrite components in protoplanetary disks.

All simulations are carried out using computer supported by the Princeton Institute of Computational Science and Engineering, and the Texas Advanced Computing Center (TACC) at The University of Texas at Austin through XSEDE grant TG- AST130002. Z. Z. acknowledges support from the National Aeronautics and Space Administration through the Astrophysics Theory Program with Grant No. NNX17AK40G and support from the Sloan Research Fellowship. We thank Kengo Tomida and Christopher J. White for their contributions to the Athena++ code. We thank Jerome Guilet and Gordon Ogilvie for sharing their figures with us. We also thank Satoshi Okuzumi, Roman Rafikov, and Sebastien Fromang for very helpful discussions and comments. Finally, we thank the referee for very helpful suggestions. To prevent the loss of magnetic fields at the polar region, we implement a special polar boundary in the $\theta$ direction. The implementation of this boundary depends on the domain size in the $\phi$ direction. We will first describe this boundary condition when the domain extends over the full $2\pi$ in the $\phi$ direction. In this case, at the pole, the ghost zones in the $\theta$ direction at (ri, $\theta_{j=-1}$, $\phi_{k}$) overlap with the active zones at (ri, $\theta_{j=0}$, $2\pi-\phi_{k}$) where $-1$ and 0 refer to the ghost and active zone. Thus, we first assign all cell-centered and face-centered quantities in the ghost zones using the quantities from the corresponding active zones at the opposite $\phi$ position. Due to the positive direction adopted in the spherical-polar system, cell-centered $v_{\theta}$, $v_{\phi}$ and face-centered $B_{\theta}$, $B_{\phi}$ need to flip their signs when they are copied to the ghost zones. Since the first active zone in the $\theta$ direction basically has a zero area at the pole, we do not use the CT scheme to update $B_{\theta}$ at the pole. Instead $B_{\theta}$ there is the average of $B_{\theta}$ from the second active zone and the first ghost zone. Finally, $E_{r}$ at the pole is taken as the average of $E_{r}$ from all the grids touching the pole and the averaged value is shared among all these grids. If the domain only extends over a wedge in the $\phi$ direction, all above steps are the same except the first step where quantities at the ghost zones are directly copied from the active zones at the same $\phi$ slice instead of copying from the $2\pi-\phi$ slice. For mesh-refinement , the prolongation and restriction methods in Toth & Roe (2002) can preserve both $\nabla\cdot\bf{B}$ and $\nabla\times\bf{B}$ in Cartesian coordinates. We adopt the same methods for mesh-refinement with spherical-polar coordinates. Theoretically, we could follow the steps in Toth & Roe (2002) to design algorithms preserving both $\nabla\cdot\bf{B}$ and $\nabla\times\bf{B}$ under the spherical-polar coordinate system. However, due to the additional geometrical factor in $\nabla\cdot\bf{B}$ and $\nabla\times\bf{B}$, the final restriction function is very complicated. Thus, we only require the restriction to satisfy $\nabla\cdot\bf{B}$=0, which is crucial for CT. It can be shown that Equation (8)-(12) in Toth & Roe (2002) can guarantee $\nabla\cdot\bf{B}$=0 even in cylindrical and spherical-polar coordinates as long as $\nabla\cdot\bf{B}$=0 initially. Note that it cannot conserve $\nabla\cdot\bf{B}$ if $\nabla\cdot\bf{B}$ is not zero initially, which is different when this method is used in Cartesian coordinates.