Rings and gaps produced by variable magnetic disk winds and avalanche accretion streams: I. Axisymmetric resistive MHD simulations

TL;DR

This study uses axisymmetric resistive MHD simulations to explore how variable magnetic disk winds and avalanche accretion streams create rings and gaps in protoplanetary disks, influencing planet formation processes.

Contribution

It demonstrates how different accretion modes driven by magnetic fields and resistivity produce observable disk structures like rings and gaps, advancing understanding of disk evolution.

Findings

Fast avalanche streams lead to dense rings through continuous feeding.

Regions with low mass-to-flux ratios form gaps, high ratios form dense rings.

Both wind-driven and stream-driven accretion modes produce observable disk features.

Abstract

Rings and gaps are being observed in an increasing number of disks around young stellar objects. We illustrate the formation of such radial structures through idealized, 2D (axisymmetric) resistive MHD simulations of coupled disk-wind systems threaded by a relatively weak poloidal magnetic field (plasma-$\beta \sim 10^3$). We find two distinct modes of accretion depending on the resistivity and field strength. A small resistivity or high field strength promotes the development of rapidly infalling `avalanche accretion streams' in a vertically extended disk envelope that dominates the dynamics of the system, especially the mass accretion. The streams are suppressed in simulations with larger resistivities or lower field strengths, where most of the accretion instead occurs through a laminar disk. In these simulations, the disk accretion is driven mainly by a slow wind that is typically accelerated by the pressure gradient from a predominantly toroidal magnetic field. Both wind-dominated and stream-dominated modes of accretion create prominent features in the surface density distribution of the disk, including rings and gaps, with a strong spatial variation of the magnetic flux relative to the mass. Regions with low mass-to-flux ratios accrete quickly, leading to the development of gaps, whereas regions with higher mass-to-flux ratios tend to accrete more slowly, allowing matter to accumulate and form dense rings. In some cases, avalanche accretion streams are observed to produce dense rings directly through continuous feeding. We discuss the implications of ring and gap formation driven by winds and streams on grain growth and planet formation.

Tables

Table 1. Table 1: Model parameters for all simulation runs. Note: the diffusivity parameters D 0 subscript 𝐷 0 D_{0} and D m , 0 subscript 𝐷 𝑚 0 D_{m,0} are measured at the inner edge of the disk, and the simulation ref-x2grid has a higher resolution in the disk region (see Section  2.3 ).

	$ϵ$	$\beta$	$\eta \left[{\mathrm{cm}}^2{s}^{-1}\right]$	$D_0$	$D_{m,0}$
reference	0.05	${10}^3$	$2.5\times {10}^{15}$	0.16	3.6
D_4	0.05	${10}^3$	$6.25\times {10}^{14}$	0.04	0.89
D4	0.05	${10}^3$	${10}^{16}$	0.64	14
beta_3	0.05	$3.33\times {10}^2$	$2.5\times {10}^{15}$	0.16	2.1
beta3	0.05	$3.0\times {10}^3$	$2.5\times {10}^{15}$	0.16	6.2
beta100	0.05	${10}^5$	$2.5\times {10}^{15}$	0.16	36
t4	0.1	${10}^3$	${10}^{16}$	0.16	3.6
ref-x2grid	0.05	${10}^3$	$2.5\times {10}^{15}$	0.16	3.6

Equations

$$\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho \mathbf{v}\right)=0,$$

$$\rho \frac{\partial \mathbf{v}}{\partial t}+\rho \left(\mathbf{v}\cdot \nabla \right)\mathbf{v}=-\nabla P+\mathbf{J}\times \mathbf{B}/c-\rho \nabla {\Phi }_g,$$

$$\frac{\partial \mathbf{B}}{\partial t}=\nabla \times \left(\mathbf{v}\times \mathbf{B}-\eta \nabla \times \mathbf{B}\right),$$

$$\frac{\partial e}{\partial t}+\nabla \cdot \left(e\mathbf{v}\right)=-P\nabla \cdot \mathbf{v}.$$

$${\rho }_d(r,\theta )={\rho }_0{\left(\frac{r}{r_0}\right)}^{-\alpha }\exp \left(-\frac{{\cos }^2\theta }{2{ϵ}^2}\right),$$

$$P_d(r,\theta )={\rho }_d{c}_s^2,$$

$$v_{\varphi }=v_K\sqrt{1-(1+\alpha ){ϵ}^2}.$$

$${\rho }_c(r)={\rho }_0{ϵ}^2(1+\alpha )\exp \left[-\frac{{\cos }^2{\theta }_0}{2{ϵ}^2}\right]{\left(\frac{r}{r_0}\right)}^{-\alpha }={\rho }_{c,0}{\left(\frac{r}{r_0}\right)}^{-\alpha },$$

$$P_c(r)={\rho }_c{v}_K^2/(1+\alpha ).$$

$$\Psi (r,\theta )=\frac{4}{3}{r}_0^2{B}_{p,0}{\left(\frac{r\sin \theta }{r_0}\right)}^{3/4}\frac{m^{5/4}}{{\left(m^2+{\cot }^2\theta \right)}^{5/8}},$$

$$B_r=\frac{1}{r^2\sin \theta }\frac{\partial \Psi }{\partial \theta },$$

$$B_{\theta }=-\frac{1}{r\sin \theta }\frac{\partial \Psi }{\partial r}.$$

$$\Delta \dot{M}=\rho \Delta A\left(\frac{{\mathbf{v}}_p\cdot {\mathbf{B}}_p}{B_p}\right)=\rho \Delta \Psi \left(\frac{{\mathbf{v}}_p\cdot {\mathbf{B}}_p}{B_p^2}\right).$$

$$\Delta M=\int \rho \Delta Ads=\Delta \Psi \int (\rho /B_p)ds,$$

$$t_{\text{dep}}=\frac{\Delta M}{\Delta \dot{M}}=\int (\rho /B_p)ds{\left(\frac{\rho {\mathbf{v}}_p\cdot {\mathbf{B}}_p}{B_p^2}\right)}^{-1}.$$

$$t_{\text{acc}}=r/v_r.$$

$${\alpha }_{\text{eff}}=-v_r{v}_K/c_s^2,$$

Full text

††pubyear: 2017††pagerange: Rings and gaps produced by variable magnetic disk winds and avalanche accretion streams: I. Axisymmetric resistive MHD simulations–References

1 Introduction

Disks around young stellar objects (YSOs) are the birthplaces of planets. Planet formation and evolution thus depend critically on the properties of these disks. Given the large number of super-Earths and mini-Neptunes discovered by Kepler at distances between $\sim 0.1$ and $\sim 1$ au from their host stars (Winn & Fabrycky, 2015), there is strong motivation for studying the structure and evolution of the inner ($\lesssim 1$ au) parts of protostellar disks. This will be the focus of our investigation.

The inner circumstellar disk plays an important role in launching outflows. Jets and winds are ubiquitously observed in YSOs and have a long history of observational and theoretical study (see Frank et al. 2014 for a review and references therein). It is thought that such outflows are driven by rapidly rotating magnetic field lines, although where the outflow-driving field lines are anchored remains unclear. One school of thought posits that the field lines are trapped at the inner edge of the disk, giving rise to so-called ‘X-winds’ (Shu et al., 2000). Another suggests that they are distributed over a more extended region of the inner disk, driving a ‘disk wind’ (Konigl & Pudritz, 2000). Part of the outflow may also be attributed to an accretion-enhanced stellar wind (e.g., Matt & Pudritz 2008). Some evidence favoring disk winds comes from the observed gradients in the line-of-sight velocity measured across optical jets (e.g., Bacciotti et al. 2000) and molecular outflows (e.g., Launhardt et al. 2009). If such gradients arise from rotation in a magnetocentrifugal wind (Blandford & Payne, 1982), the implied outflow launching region would extend to au scales or larger (e.g., Anderson et al. 2003; Ray et al. 2007; Bjerkeli et al. 2016; however, see Lee et al., submitted, for SiO observations of HH 212 that provide evidence for X-winds). The outflow would carry angular momentum away from the disk in a launching region that is directly relevant to the formation of the terrestrial planets in the Solar system and a considerable fraction of the exoplanets discovered by Kepler.

The idea of magnetic wind-driven disk evolution has been discussed in the literature for a long time (e.g., Pudritz & Norman 1983; Wardle & Koenigl 1993; Combet & Ferreira 2008; Königl et al. 2010; Bai 2016; Suzuki et al. 2016, and references therein). One way to quantify the effects that winds have on disks is to construct global, coupled wind-disk solutions to the MHD equations. Many of the early investigations along this line adopted simplifying self-similarity assumptions (e.g., Wardle & Koenigl 1993; Li 1995; Ferreira 1997; Salmeron et al. 2011). These semi-analytic treatments have been important for illuminating the basic mechanics of the coupled system, including how a small fraction of the accreting material is peeled off of the disk surface and ejected along the rapidly rotating field lines as a wind, how the wind removes angular momentum from the disk and drives it to accrete, and how the accretion and ejection processes depend on the non-ideal MHD effects (Ohmic dissipation, ambipolar diffusion, and the Hall effect) that are expected to occur in the lightly ionized YSO disks. Self-similar solutions require that all physical quantities, such as the density, magnetic field strength, and the coefficients for non-ideal MHD effects, must scale with the radius as power-laws with specific indices. This requirement limits the applicability of the self-similar solutions to realistic systems. We seek to relax the self-similarity ansatz through a series of MHD simulations that incorporate non-ideal effects.

Simulations of coupled wind-disk systems have been carried out by several groups (e.g., Kato et al. 2002; Casse & Keppens 2002, 2004; Zanni et al. 2007; Murphy et al. 2010; Lii et al. 2012; Sheikhnezami et al. 2012; Čemeljić et al. 2013; Stepanovs & Fendt 2014, 2016). These simulations confirmed and extended the earlier semi-analytic results, finding self-consistent wind-disk solutions for a wide range of disk masses and magnetic field distributions. However, the focus of such simulations has typically been on the launching of outflows from the disk and on their propagation to large distances. We are thus motivated to start a long-term program to investigate the wind-disk system that will eventually include all three non-ideal MHD effects using the ZeusTW MHD code (Krasnopolsky et al., 2011; Li et al., 2011), focusing in particular on the structure of the disk in such a system. As a first step, we start with the simplest case of a resistive wind-disk system under the assumptions of axisymmetry (2D) and reflection symmetry across the disk midplane.

Resistive disk winds under such an idealized geometry have been studied by several groups, especially Stepanovs & Fendt (2016) (herein, SF16), who were able to determine wind properties (such as its mass loss rate, energy, and angular momentum) for disks with a wide range of magnetic field strengths. The emphasis of such work has been on quasi-steady state wind-disk solutions, which enable the computation of several well-known conserved quantities along each magnetic flux surface (e.g., Spruit 2010) that can be checked against the semi-analytic solutions. However, a steady state can only be obtained when outward diffusion of the magnetic field in the disk is balanced by inward advection of the field, which occurs only under rather restrictive conditions for the prescribed magnetic diffusivity. Under more general conditions, the magnetic winds remain highly variable, driving non-steady disk accretion even for 2D (axisymmetric) systems; they are likely to be exacerbated in 3D simulations. Such non-steady state phenomena are the focus of this investigation. We find that the disk structure can be completely reshaped from its initial state by even a weak initial magnetic field, often with rapidly accreting streams developing near the disk surface and rings and gaps developing near the disk midplane.

The accretion streams are a form of the magnetorotational instability (MRI) channel flows (Balbus & Hawley, 1991; Hawley & Balbus, 1992; Goodman & Xu, 1994). They are the result of the runaway magnetic braking of an infalling stream, which is linked to more slowly rotating material at larger radii. We term them ‘avalanche accretion streams’, motivated by the work of Kudoh et al. (1998) who find similar structures in ideal MHD simulations of thick AGN disks (see also Stone & Norman 1994; Matsumoto et al. 1996; Kudoh et al. 2002; Beckwith et al. 2009). We find that the accretion streams become more numerous and more important dynamically as the disk material becomes better coupled to the magnetic field. Their rapid formation and disruption forms a thick envelope of chaotic infall and outflow motions above the disk, which is intimately related to the so-called ‘coronal accretion’ found recently by Zhu & Stone (2017) (herein, ZS17) in global (3D) ideal MHD simulations of thin disks with net vertical magnetic flux (see also Beckwith et al. 2009). These investigations highlight the possible importance of vigorous accretion in a vertically extended structure outside the dense disk on disk evolution.

The formation of rings and gaps on sub-au scales strongly influences how dust grains are trapped and thus affects how the grains grow into planetesimals and ultimately planets. Although our simulations are limited to the inner ($\sim 0.01-1$ au scale) disk, the same mechanism of ring and gap formation through variable disk winds and avalanche accretion streams should operate at larger radii as well, where rings and gaps have now been observed in a number of disks (ALMA Partnership et al., 2015; Andrews et al., 2016; Zhang et al., 2016; Isella et al., 2016). The 2D simulations presented in this study serve as an illustration of these generic mechanisms, although they are likely to generate rings and gaps more easily than 3D simulations, because of the assumed axisymmetry. Nevertheless, they provide a starting point for future explorations that will include more detailed disk microphysics and less restrictive geometry.

The rest of the paper is organized as follows. In Section 2, we describe the simulation setup. In Section 3, we present the results of a reference simulation. In Section 4, we analyze how the simulation outcome depends on three key dimensionless parameters that control the disk magnetic field strength, resistivity and disk thickness. In Section 5, we discuss this work in the context of previous studies and examine its implications. The results are summarized in Section 6.

2 Problem setup

2.1 MHD equations

We use the ZeusTW code (Krasnopolsky et al., 2010) to solve the time-dependent, resistive magnetohydrodynamic (MHD) equations in axisymmetric spherical coordinates ($r,\theta,\phi$). The equations solved are as follows:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The current density is $\bm{J}=(c/4\pi)\nabla\times\bm{B}$, the internal energy is $e=P/(\gamma-1)$, and $\eta$ is the resistivity. The remaining parameters have their usual definitions. When referring to cylindrical coordinates, we will use the notation ($R,\phi,z$) such that $R=r\sin{\theta}$ and $z=r\cos{\theta}$.

2.2 Initial conditions

The simulation domain is separated into two regions: a thin, cold, rotating disk orbiting a 1 M⊙ central source at the grid origin and a stationary, hot corona above the disk. We choose $\gamma=1.01$ so that the disk and corona regions are locally isothermal.

2.2.1 Disk

The geometrically thin disk is characterized by the dimensionless parameter $\epsilon=h/r=c_{s}/v_{K}\ll 1$, where $h$ is the disk scale height, $c_{s}$ is the sound speed, and $v_{K}$ is the Keplerian speed. The disk is limited to the equatorial region where the polar angle $\theta\in{[\pi/2-\theta_{0},\pi/2]}$, with disk (half) opening angle set to $\theta_{0}=\arctan(2\epsilon)$. The disk density takes the form of a radial power law multiplied by a Gaussian function of $z/r=\cos\theta$,

[TABLE]

as dictated by hydrostatic balance. The subscript ‘0’ refers to values on the disk midplane at the inner radial boundary. For all simulations shown this paper, we use $\alpha=3/2$. The choice of power-law exponent is consistent with sub-millimeter observations of Class II sources that find surface density power-law exponents of $0.4-1.1$ (Andrews et al., 2010). The disk pressure is set as

[TABLE]

with $c_{s}=\epsilon v_{K}$. The radial pressure gradient causes the equilibrium rotational velocity $v_{\phi}$ to be slightly sub-Keplerian as

[TABLE]

2.2.2 Corona

We require that the hydrostatic corona is initially in pressure balance with the disk surface. This constraint sets the density drop from the disk surface to the corona as $(1+\alpha)\epsilon^{2}$. Therefore, the coronal density and pressure are

[TABLE]

[TABLE]

The coronal sound speed is then $c_{s,c}=v_{K}/\sqrt{1+\alpha}$.

2.2.3 Magnetic field

To ensure that the magnetic field is divergence-free initially, we set the magnetic field components using the magnetic flux function $\Psi$ as in Zanni et al. (2007),

[TABLE]

where the $B_{p,0}$ is the poloidal field strength on the midplane at $r_{0}$, and the parameter $m$ controls the bending of the field. Since varying $m$ from 0.1 to 1 has little effect on the long-term disk or wind magnetic field structure (Stepanovs & Fendt, 2014), we use $m=0.5$ for all simulations presented in this work. The initial magnetic field components are then calculated as

[TABLE]

[TABLE]

2.2.4 Resistivity

For illustrative purposes in this initial study, we adopt a spatially and temporally constant resistivity $\eta(\bm{r},t)=\eta$; more refined treatments, including one based on detailed ionization calculation, are postponed to future investigations. This implementation differs from the resistivity profile used in other disk wind simulations where $\eta\propto hv_{A}$ ($v_{A}$ is the Alfvén speed) inside the disk and $\eta=0$ outside (e.g., Stepanovs & Fendt 2014). To quantify the relative importance of the prescribed resistivity, we define a dimensionless magnetic diffusivity parameter $D\equiv\eta/(hc_{s})$ (where $h$ is the disk scale height and $c_{s}$ is the disk sound speed) as in Li (1995), similar to the well-known $\alpha$-parameter of Shakura & Sunyaev (1973). This diffusivity parameter is initially proportional to $r^{-1/2}$ with $D_{0}=0.16$ at the inner edge of the disk in the reference simulation, corresponding to an ionization fraction of $x_{e}=2.3\times 10^{-12}$ at 1 au (Blaes & Balbus, 1994). This is not far from the minimum ionization fraction due to the radioactive decay of 26Al alone ($x_{e}\approx 1.9\times 10^{-13}$ according to Armitage 2011). The resistivity adopted in our reference simulation is large enough to suppress the fastest growing MRI mode near the disk midplane (Jin, 1996). However, the MRI may still operate, especially near the disk surface where the plasma-$\beta$ is lower and the critical resistivity needed for MRI suppression (which is inversely proportional to $\beta$) is higher.

2.3 Grid

The equations are solved for $r\in{[0.02,10]}$ au and $\theta\in{[0,\pi/2]}$, where $\theta=0$ corresponds to the polar axis and $\theta=\pi/2$ to the disk midplane. We choose $r_{0}=0.02~{}\rm{au}\sim 4~{}\rm{R}_{\odot}$ as the inner radius because it is approximately the inner radius of a protostellar disk truncated by the magnetosphere of a classical T Tauri star (Hartmann et al., 2016). It yields an inner orbital period $t_{0}\approx 1$ d. For this initial study, we focus on only one hemisphere and assume that the system has reflection symmetry across the disk midplane. This assumption will be relaxed in future investigations. In the radial direction, a ‘ratioed’ grid is used, such that $dr_{i+1}/dr_{i}$ is constant and $r_{i+1}=r_{i}+dr_{i}$. The grid spacing at the inner edge is set such that $dr_{0}=2.3r_{0}d\theta$. The grid is uniform in the $\theta$-coordinate for most simulations presented in the work, with a resolution of $n_{r}\times n_{\theta}=400\times 360$. This results in 24 grid cells from the disk midplane to surface (at two disk scale heights) in the reference run. The fastest growing (ideal MHD) MRI wavelength would be marginally resolved at the disk midplane with $\lambda=3.3rd\theta$ initially (Balbus & Hawley, 1991).

In the simulation called x2-grid, the $\theta$-grid is separated into two regions: a uniform grid in the disk region ($\theta\in{[\pi/3,\pi/2]}$), and a ratioed grid in the corona ($\theta\in{[0,\pi/3]}$), with $d\theta$ matched at the boundary and increasing towards the polar axis. The uniform region contains 240 cells, thereby decreasing $d\theta$ by a factor of two compared to the reference run. The non-uniform grid section has 120 cells, so the total number of cells in the $\theta$-direction remains 360. We decrease the value of $dr_{0}$ in order to have the ratio of $\frac{dr_{0}}{r_{0}d\theta}$ be the same in all runs. The value of $d\theta$ increases to $1.3^{\circ}$ at the polar axis for this ratioed grid. In this case, the fastest growing (ideal MHD) MRI wavelength would be resolved by about 7 cells.

2.4 Boundary conditions

The inner radial boundary is separated into two regions. For $\theta\in{[0,\pi/2-\arctan(4\epsilon)]}$, mass is injected into the simulation domain with $\rho=\rho_{c,0}$ and $v=10~{}\text{km}~{}\text{s}^{-1}$ (Krasnopolsky et al., 1999). The injection boundary is used to lower the Alfv́en speed in polar region close to the inner boundary and thereby increase the simulation time step, as the region would quickly be evacuated by gravitational infall if the standard outflow boundary was used. The total mass that enters the simulation domain through the inner boundary is small compared to both the initial mass in the simulation and the mass that is eventually carried through the outer boundary via the disk wind. It thus has little influence on the dynamics of the simulation. The remaining section of the inner radial boundary, $\theta>\pi/2-\arctan(4\epsilon)$, and the outer radial boundary both use the standard outflow condition, as usual in Zeus codes.

The reflection boundary condition is used on the polar axis ($\theta=0$) and the equatorial plane ($\theta=\pi/2$). The $\phi$-component of the magnetic field is set to zero on the polar axis. The assumed reflection symmetry across the midplane ($\theta=\pi/2$) demands that $B_{\phi}=0$ at this boundary. We also set $B_{\phi}$ to vanish on the inner radial boundary because it is not rotating.

3 Reference model

We start by discussing the result of a ‘reference’ simulation. It is used to highlight the salient features of the coupled wind-disk system, and serves as a benchmark against which other simulations with different parameters will be compared in Section 4. Of the parameters that will be changed in Section 4, the reference model has $\epsilon=0.05$, an initial magnetic field strength characterized by plasma-$\beta$ of $10^{3}$ on the disk midplane (corresponding to 0.27 G at 1 au), and a dimensionless magnetic diffusivity parameter of $D_{0}=0.16$. This simulation runs for 1820 inner orbital periods ($\sim 5$ yr). Fig. 1 shows the initial conditions for the system and snapshots at three representative times.

3.1 Global evolution

The initially weak poloidal magnetic field is wound up by differential rotation between the (nearly) Keplerian disk and the static corona, which inflates a bubble of strong toroidal magnetic field that expands outward against the coronal material. This winding operates fastest near the inner edge of the disk, where the orbital period is the shortest. As a result, the outflow propagates with the highest speed near the polar axis. By $t=150t_{0}$ (Fig. 1(b)), the outflow has reached the outer edge of the computation domain along the axis, but remains confined by the initial coronal material away from the polar region. By $t=1200t_{0}$ (Fig. 1(c)), most of the initial coronal material has been completely swept out of the computational domain. Beyond this time, the effects of the initial corona on the coupled wind-disk system should be relatively small, and for the inner part of the system (where much of our analysis is focused; see Section 3.3.3 below), the initial corona ceases to be important at an even earlier time. The wind-disk system shows large spatial and temporal variability throughout the simulation. This variability is reflected in Fig. 1(d) at $t=1800t_{0}$; its appearance is broadly similar to that in panel (c), but there are important differences such as the appearance of low-density wind ‘channels’ that are more prominent in the outflow of the former than the latter. The variability of the outflow is intimately tied to the structures that develop in the disk.

3.2 Outflow

Even though the initial magnetic field at the disk midplane is rather weak ($\beta=10^{3}$), it is still able to drive a powerful outflow. Unlike in the disk, the pressure in the outflow is dominated by the magnetic field rather than the thermal gas, as illustrated in Fig. 2(a). This dominant magnetic pressure is provided mostly by the toroidal component $B_{\phi}$, as shown in panel (b) where the ratio of the toroidal and poloidal components $|B_{\phi}/B_{p}|$ is plotted. This ratio is spatially inhomogeneous, with a filamentary appearance. The dominance of the toroidal field component suggests that the outflow is driven mostly by the magnetic pressure gradient, as is expected for a relatively weak initial magnetic field (e.g., SF16) or a heavy mass loading (e.g., Anderson et al. 2005). The exceptions are a narrow region near the polar axis and two filaments at roughly $\theta=10^{\circ}$ and $80^{\circ}$. The magnetic field lines that run through the polar axis region are connected to the inner radial boundary, which is assumed to be non-rotating; the lack of a toroidal field component here is therefore expected. The other two filaments are physically meaningful and correspond to the two low-density ‘channels’ that are prominent in the density map of Fig. 1(d). In these channels, the thermal pressure is completely dominated by the magnetic pressure, with a plasma $\beta$ less than $10^{-3}$ (see Fig. 2(a)). In other words, their field lines are much more lightly mass-loaded than in the rest of the outflow. As a result, the magnetic field is able to accelerate the mass to a much higher speed than in the denser surrounding regions (see Fig. 2(c)).

The distribution of the poloidal velocity is plotted in panel (c) of Fig. 2. The outflow in the outer low-density channel reaches a speed of $\sim 50$ km s*-1*, whereas that in the inner channel is much faster, reaching up to 200 km s*-1*. The outflow close to the axis can also reach a relatively high speed of $\sim 50-100$ km s*-1*. However, over most of the simulation volume, especially to the right of the inner channel, the poloidal outflow speed is rather low, typically a factor of two or more smaller than the the rotation speed and well below the local free-fall speed (see Fig. 2(d)). Except for the channels and polar axis region, the magnetic and velocity fields of the slowly-expanding, low plasma-$\beta$ outflow are dominated by their toroidal components.

Even though the outflow is slow in general, there is disk material that gets accelerated to a high speed. In Fig. 3, we plot the mass flux of the outflowing material with a radial velocity greater than 100 km s*-1* as a function of time through a hemisphere of radius of $r=1$ au. It is clear that the fastest component of the outflow is highly variable in time. The mass outflow rate routinely changes by an order of magnitude on timescales of $t\lesssim 100t_{0}$, and can sometimes drop by four orders of magnitude on similar timescales. This variability may be essential for generating internal shocks that are required to keep the outflow heated to a relatively high temperature and visible through, e.g, optical forbidden lines (e.g., Shang et al. 2002). The average mass loss rate of this fast component is approximately $5\times 10^{-8}$ M⊙ yr*-1*, which is at the high end of the mass loss rate inferred in the classical T Tauri jets (e.g., Frank et al. 2014). This fast, massive outflow component is remarkable in light of the fact that the disk is only weakly magnetized, at least initially with a midplane plasma-$\beta$ of $10^{3}$.

The fast component is only a relatively minor component of the outflow produced in the reference simulation. This is illustrated in Fig. 3, where the mass loss rate is plotted as a function of time for the slow (1 km s*-1*$<v_{r}<$ 10 km s*-1*) and intermediate-speed (10 km s*-1*$<v_{r}<$ 100 km s*-1*) components of the outflow. Both components have mass loss rates of order $10^{-6}$ M⊙ yr*-1*, which is at least an order of magnitude higher than the mass flux of the fast component. Such a massive outflow is expected to strongly affect – indeed, control – the disk dynamics, as we demonstrate in Section 3.3.1.

3.3 Disk-wind connection

Fig. 4 shows the density distribution and magnetic field lines for the final time frame of the reference simulation, as in panel (d) of Fig. 1 and 2, but now zoomed in on the inner part of disk where the disk structure is most strongly affected by the outflow. As expected, the modification of disk structure starts near the inner edge, where the rotation period is the shortest. By $t=1800~{}t_{0}$ (Fig. 4), three distinct regions have been established: (1) an inner region ($r<0.1$ au) with a highly variable mass distribution that contains prominent dense rings and gaps, (2) an intermediate radius region (0.1 au $<r<0.5$ au) with a much smoother mass density distribution, and (3) an outer region ($r>0.5$ au) where the disk ‘puffs up’ due to the presence of a dense surface layer – the ‘avalanche accretion stream.’ We label these regions Zone I, II, and III, respectively, and discuss each in detail below. We save the discussion of the most variable region, Zone I, until Section 3.3.3.

3.3.1 Zone II: Steady wind-driven accretion

We start our discussion with the most laminar, intermediate radius region, Zone II. The location of this region drifts radially outward as the simulation progresses. By $t=1800t_{0}$ (Fig. 4), it is located roughly between 0.1 and 0.5 au. The magnetic field line intersecting the midplane at a radius of 0.3 au is marked by a dashed line in Fig. 4 (left panel). Fig. 5(a) plots the density distribution along that field line, with the location of the sonic point marked. For comparison, we also plot the density profile expected for an isothermal (thin) disk with the temperature set to that at the midplane. Clearly, the presence of the magnetic field and the launching of a magnetized wind have significantly changed the vertical density profile of not only the disk itself, but also the region surrounding the disk. Specifically, the vertical density gradient is steeper near the disk surface compared to an isothermal (non-magnetic) profile. This difference is explained by a decreasing temperature along the path and by an increasing toroidal component of the magnetic field, $B_{\phi}$, which reaches a maximum near 0.02 au before decreasing again (Fig. 5(b)). Indeed, the toroidal field becomes so strong near the disk surface that it dominates the thermal pressure (see Fig. 5(f) and discussion below), and generates a downward magnetic pressure force that compresses the disk significantly. Above the disk surface, the magnetic pressure force has a positive radial component conducive to launching a wind.

The transition from disk to wind can be seen most clearly in the distribution of the radial component of the velocity $v_{r}$ along the field line (Fig. 5(c)). The transition occurs approximately $\sim 0.03$ au above the disk midplane. Material below this point moves radially inward with a speed of $\sim$0.2 km s*-1*, while material above this point is gradually accelerated outward. Beyond the sonic point the flow speed becomes more variable, because the wind has yet to reach a steady state. Even though the wind speed becomes supersonic at large distances, it remains well below the rotational speed except toward the edge of the simulation domain (see Fig. 5(d)). In fact, $v_{\phi}$ is close to the speed needed to be rotationally supported against the gravitational pull from the central star in the cylindrical radial ($R$) direction (also shown in Fig. 5(d) for comparison). The implication is that, just like the disk, the slowly expanding wind is almost entirely supported in the $R$-direction by rotation. The difference lies in the vertical direction; whereas the disk is supported against gravity (and magnetic compression) mostly by the thermal pressure gradient, the wind is supported mostly by the magnetic pressure gradient. This difference is highlighted in Fig. 5(e), which shows that the ratio of the vertical component of the Lorentz force to the gravitational force changes sign at the disk surface/wind base and reaches unity just before the sonic point. The available evidence points to a magnetically dominated wind (Fig. 5(f)) that is slowly lifted out of the (deep) gravitational well near the disk surface by the magnetic pressure gradient, consistent with the ‘magnetic tower’ picture of Lynden-Bell (1996).

Despite being much slower than the Keplerian speed at its footpoint, the wind launched from 0.3 au strongly affects the disk structure and dynamics. First, it rapidly depletes the disk material. To estimate a local disk mass depletion time, we define a magnetic flux tube containing an infinitesimally small magnetic flux $\Delta\Psi$ around the field line. The cross-sectional area of the flux tube is $\Delta A=\Delta\Psi/B_{p}$, where $B_{p}$ is the poloidal magnetic field. The rate of wind mass flux along the magnetic flux tube is

[TABLE]

This is to be compared with the amount of mass contained inside the flux tube

[TABLE]

where the integration is along the field line, and is dominated by the mass inside the disk. The local disk mass depletion time due to wind mass loss is then

[TABLE]

In a steady state, the quantity $\Delta\dot{M}$ should be constant along a field line. However, the wind does not reach a strict steady-state, and the non-steadiness introduces some variation to that combination, especially at large distances. We obtain a depletion time of $t_{\text{dep}}\approx 9.6$ yr, which is 58 times the local orbital period at the footpoint. This time scale is to be compared with the accretion time scale for the disk material, which is roughly

[TABLE]

For an average accretion speed of $\sim 0.2$ km s*-1* at 0.3 au, the accretion timescale is $t_{\text{acc}}\approx 8.6$ yr, or 52 times the local orbital period. The fact that the accretion time is comparable to the depletion time means that a significant fraction of the disk will be ejected in the wind.

One way to quantify the disk accretion rate is through an effective $\alpha$ parameter, defined as

[TABLE]

which, for $v_{r}\approx 0.2$ km s*-1*, $v_{K}=54$ km s*-1*, and $c_{s}=2.3$ km s*-1*at 0.3 au, yields $\alpha_{\text{eff}}\approx 2$. This is much larger than the effective $\alpha$ typically obtained from turbulent MRI simulations ($\alpha\sim 1/\beta$, Hawley et al. 2015), which highlights the dynamical importance of the slow disk wind launched even by a relatively weak initial poloidal magnetic field, in that it can still drive rapid disk accretion.

3.3.2 Zone III: Avalanche accretion stream and slow midplane decretion

For the outer disk region ($r>0.5$ au), we plot the physical properties of the gas and magnetic field along a representative field line that starts from a radius of $\sim 1$ au on the disk midplane (Fig. 6). Panel (a) of Fig. 6 shows the density distribution as a function of the distance along the field line. We can see that near the midplane, the density profile is close to that expected for a (non-magnetic) isothermal disk (dashed line), indicating that the magnetic field is not important in the vertical force balance. Moving closer to the disk surface (where the magnetic pressure dominates the thermal pressure), magnetic compression causes the density to drop significantly below the isothermal value, as in Zone II. The difference is that the density in Zone III increases sharply again before dropping off precipitously. This density spike marks the ‘avalanche accretion stream’ that is prominent in Fig. 4. It is the defining characteristic of Zone III, not only in the density structure but also in the magnetic field structure.

From Fig. 4, it is clear that the dense stream occurs where the magnetic field line pinches severely, with a sharp reversal of the radial field component across it. This reversal shows up clearly in Fig. 6(b), where all components of the magnetic field are plotted. The poloidal magnetic field line near the midplane first bows outward, with a small but positive radial component, as in Zone II (see Fig. 6(b)). However, it is forced to bend sharply inward (see Fig. 4) by the rapidly infalling stream, producing a large (negative) value for $B_{r}$. The toroidal field component increases from the midplane towards the stream, dominating the other two field components before reversing direction across the stream (Fig. 6(b)).

As the material in the stream falls radially inward, it is orbiting faster relative to material both above and below it, as shown in Fig. 6(d), where we plot $v_{\phi}$ along the representative field line. Simple geometric considerations show that the faster rotation at the tip of the sharply pinched field line will twist the inward-directing field into a toroidal field of positive sign and the outward-directing field into a toroidal field of the opposite sign, as shown in Fig. 6(b). The twisted magnetic field geometry efficiently brakes the stream as it loses angular momentum to both the wind and the disk material that is magnetically linked to the stream at larger radii. Both $B_{r}$ and $B_{\phi}$ are amplified as the fast-rotating stream accretes, which leads to a stronger magnetic braking torque. The end result is a run-away collapse of the stream similar to the MRI, and hence the name ‘accretion avalanche’ (Kudoh et al., 1998). At the time shown in Fig. 6, the infall speed of the stream is supersonic and reaches about 2 km s*-1* (only the sonic point in the wind above the the stream is shown in the figure). The bulk of the disk material below the stream moves much more slowly and in the opposite direction, as shown in Fig. 4 (left panel). The contrast between the fast accreting stream and slowly expanding disk is quantified in the radial velocity profile plotted in Fig. 6(c). This ‘decretion’ is caused by the disk material receiving angular momentum from the infalling stream and is opposite of the disk motion in Zone II (compare to Fig. 5(c)).

Although the stream infalls supersonically, its infall speed is still much smaller than the rotation speed $v_{\phi}$ (Fig. 6(d)). Indeed, the rotation speed is close to the value needed for the centrifugal force to balance the gravity in the cylindrically radial direction inside both the disk and the stream, indicating that both are close to being rotationally supported. A supersonic wind is launched above the stream, reaching speeds of $\sim$50 km s*-1*, which is higher than that reached along the field line originating from 0.3 au (see Fig. 5(d)). This speed difference is mainly due to lower mass loading, which allows the wind to experience greater acceleration despite rotating more slowly near the footpoint. Nevertheless, both the velocity and magnetic field are still dominated by the toroidal component over a large fraction of the wind, which is broadly similar to the wind driven from Zone II. The thermal pressure of the wind here is also completely dominated by the magnetic pressure, again similar to the wind from Zone II. Therefore, the formation of the stream appears to modify, but not completely disrupt, the wind launching and acceleration.

3.3.3 Zone I: Rings and gaps in inner disk

The innermost disk region ($r<0.1$ au) is most strongly modified from its initial state during the simulation because it has the shortest rotation period and is located where the (poloidal) magnetic flux tends to accumulate through disk accretion. It is also the most variable region of the disk-wind system. The spatial variation of mass distribution is evident in Fig. 4, where the right panel shows a ‘face-on’ view of the axisymmetric (one hemisphere) surface density normalized by the initial surface density distribution, \Sigma_{i}=810~{}\text{g~{}cm{}^{-2}}~{}(r/r_{0})^{-1/2}, to highlight features at large radii. There is a striking contrast between the dense inner rings and their surrounding gaps; in fact, most of the mass of the inner disk is concentrated in these rings. Through the course of the simulation up to five rings are formed in Zone I. The first is formed just outside of the inner boundary and is quickly accreted after 200 inner orbital periods. The remaining rings all form in an inside-out manner; once one ring and gap are formed (with the gap located at a larger radius than the ring), another ring develops at a slightly larger radius than the first gap. Once formed, the multiple rings and gaps drift slowly inward maintaining their relative order.

Although the magnetic field topology near the rings and gaps is irregular, winds are still launched from Zone I. The winds that are launched from the gaps accelerate material to very high velocities ($v_{r}>100$ km s*-1*), both because very little mass gets loaded onto the field lines anchored in the low-density gaps and because the vertical magnetic field strength peaks in the gaps (allowing for rigid rotation of the field lines out to a larger radius). This is quantified in Fig.7(a), where the column density as a function of radius at a representative time ($t=1250~{}t_{0}$) is plotted. One might naively expect the vertical component of the magnetic field to be concentrated in the high density rings, because the mass-to-flux ratio, $\Sigma/|B_{\theta,\rm{mid}}|$, would be conserved in the strict ideal MHD limit. However, there is a tendency for the vertical field to be weaker in the dense rings than in the gaps. This anti-correlation is shown in Fig. 7(b), where $-B_{\theta,\rm{mid}}$ is plotted versus radius ($-B_{\theta,\rm{mid}}=|B_{\theta,\rm{mid}}|$ since the magnetic field points upward at the midplane). It is clear that the dense rings are less magnetized relative to both the initial disk value and the neighbouring gaps. The contrast between the rings and gaps is particularly striking for the distribution of the mass-to-flux ratio shown in Fig. 7(c).

Since the innermost rings are always observed to form first, it is natural to relate their formation to the inside-out development of the disk wind (i.e., at any given time, the disk wind is further developed at smaller disk radii because of their shorter dynamical times). The wind drains angular momentum from the disk, forcing the disk material to accrete and drag magnetic flux along with it. In the ideal MHD limit, the accumulation of mass at a given location would lead to a corresponding pile-up of magnetic flux at the same location. But in the resistive MHD of this simulation, the field lines diffuse away from the region of magnetic flux concentration. This leads to two key steps in the formation of a ring and gap. First, the magnetic flux begins to drop in the region where surface density starts to spike. The decreasing dynamical importance of the magnetic field in this region increases the accretion timescale of the dense ring, effectively creating a trap for mass accretion from larger radii. Second, the magnetic flux begins to accumulate just outside the growing overdensity. Now the low mass-to-flux region can efficiently drain disk material of angular momentum, quickly moving any remaining mass to smaller radii where it is added to the ring. The fast accretion of the material in the gaps onto the nearly stationary rings is shown in Fig. 7(d), where the radial component of the velocity at the disk mid-plane is plotted. This interplay of magnetic flux redistribution, angular momentum removal through magnetic wind, and disk accretion naturally leads to a configuration of a high mass-to-flux region (ring) just inside of a low mass-to-flux region (gap; Fig. 7(c)). Once an over-dense region is formed, its survival allows for the development of similar patterns at larger radii.

3.4 Magnetic diffusivity and the disk-wind structure

To summarize, there are three distinct zones that develop in the reference simulation. The dimensionless diffusivity parameter $D=\eta/(hc_{s})$ decreases with radius as $D\propto r^{-1/2}$. Zone I is thus the most magnetically diffusive region initially, as measured by $D$. It is also the region where the most prominent rings and gaps form, because magnetic diffusion is needed for the redistribution of magnetic flux relative to the matter. In Zone II the magnetic diffusion is such that the outward diffusion balances the inward advection of the magnetic field, thereby establishing a laminar, quasi-steady disk-wind system. The least magnetically diffusive region, Zone III, develops a dense surface accretion stream that dominates the mass accretion in that region and drives the disk material below it to expand. These results suggest that the level of magnetic diffusion is a key factor in controlling the structure and dynamics of the coupled disk-wind system.

4 Parameter study

In order to determine how robust the basic features of the reference simulation are, especially the formation of rings, gaps and avalanche accretion streams, we performed several additional simulations in which we varied three dimensionless model parameters: $D_{0}$, $\beta$, and $\epsilon$. These characterize the disk magnetic diffusivity, field strength, and disk thickness (or temperature), respectively (see Table 1).

We find that all simulations show characteristics of at least one of the three Zones described in the reference simulation and all simulations develop variable winds (see Fig. 9). Almost all simulations show enough spatial variation in the surface density that rings and gaps are readily apparent, and many of the simulations show dense avalanche accretion streams developing near the disk surface. In some cases the inner disks appear to be entirely dominated by rapid formation and break-up of short-lived accretion streams, forming a vertically extended ‘envelope’ above the disk characterized by chaotic infall and outflow motions. In the following subsections we discuss the effects of each of the three dimensionless parameters in turn, focusing in particular on the lower diffusivity (D$\_$4), stronger field (model beta$\_$3), and higher temperature (t4, or thicker disk) cases (illustrated in Fig. 8 and Fig. 9).

4.1 Resistivity

We start with model D_4, where all parameters are the same as in the reference case except for the (constant) resistivity $\eta$, which is reduced by a factor of four (thereby reducing $D_{0}$ by the same amount). The reduction in $\eta$ makes the magnetic field lines better coupled to the material in the disk. We expect this to facilitate the formation of avalanche accretion streams, as in Zone III of the reference simulation, and this is indeed the case. Movies of simulation D_4 (see the supplementary material in the online journal for animated versions of Fig. 8) show that avalanche accretion streams form continuously, starting as early as $\sim 10$ inner orbital periods. Many of the streams are pushed up high into the disk corona by the growing toroidal magnetic fields beneath them (as in the case of the stream in Zone III of the reference simulation, see Fig. 6(b)). They are eventually disrupted as they fall radially inward and towards the midplane. The constant formation and disruption of the (often elevated) avalanche accretion streams leads to a thick, clumpy envelope above the disk that is highly inhomogeneous in density, velocity, and magnetic field. This stream-produced envelope can be seen clearly in Fig. 8(a), where its maximum vertical extent is comparable to the cylindrical radius and is positioned between the denser equatorial disk and the more tenuous polar wind region. While the (poloidal) field lines in the polar wind region are well organized, those in the envelope are disordered with severe bunching of field lines in some places and looping (due to reconnection) in others. The disordered field is a reflection of the chaotic motions inside the envelope characterized by simultaneous infall and outflow. The motion of the envelope is best seen in Fig. 8(b), where the spatial distribution of mass flux per unit polar angle, $d\dot{M}/d\theta=2\pi r^{2}\rho v_{r}\sin{\theta}$, is plotted. The mass accretion in the envelope is dominated by the fragments of disrupted avalanche accretion streams (blue in middle row of Fig. 8). They are mixed together with packets of outflow (red) that sometimes extend to the disk midplane. The picture is reminiscent of the MRI-driven channel flows that are seen, e.g., in the thick disk simulations of Kudoh et al. (2002), and they share the same physical origin – the exchange of angular momentum between magnetically connected material located at different radii. In our thin-disk case, the difference is that the alternating pattern of inflow and outflow occurs mostly in the envelope above the disk rather than inside the disk. This dynamic envelope is a new feature that formally does not exist in the reference simulation, although it is intimately related to the avalanche accreting stream of Zone III.

In addition to the envelope, a wind is also launched in the low-resistivity case (model D_4; see also Suzuki & Inutsuka 2009; Fromang et al. 2013). The wind can be seen in Fig. 8(a), especially in (but not limited to) the polar region. The wind is highly variable, as illustrated in Fig. 9 (blue line), where the mass loss rate for the fastest wind component ($v_{r}>100$ km/s) is plotted as a function of time. The mass loss rate is comparable to that in the reference case, especially in the second half of the simulation. The highly variable wind and the chaotic envelope create radial structures in the disk. Fig. 8(c) shows a face-on view of the surface density normalized by its initial distribution. Several rings and gaps are clearly visible, and they are quantified in Fig. 10(b), where the surface density is plotted as a function of radius. The development of rings and gaps starts at small radii, and gradually spreads to large radii, because the dynamical time increases with radius. Some of the rings appear to change quickly, while others are more stationary. These features are present throughout the simulation, suggesting that they are a robust characteristic of the coupled disk-envelope-wind system. We find that the surface density anti-correlates with the poloidal field strength on the disk, with the dense rings typically less magnetized than the gaps, as in Zone I of the reference simulation. This indicates that magnetic flux redistribution through either the relatively small resistivity or turbulent reconnection, is playing a role in creating the rings and gaps. Since the disk is tightly connected to the stream-dominated envelope, the streams are also expected to play an important role in gap and ring formation. This role is difficult to quantify precisely, however, because of the chaotic flow pattern in the disk and envelope. A cleaner case of stream-induced ring formation will be presented in the higher-temperature case below (see Section 4.3). We have also carried out an ideal MHD simulation with the same parameters except for $\eta=0$, and found results broadly similar to this low-resistivity case.

In contrast to D_4, the more resistive model D4 is much more laminar. A wind is launched steadily throughout the run over most of the disk radii that have had time for at least one orbit to take place, as in the smooth Zone II of the reference run. There are no avalanche accretion streams in the simulation. A new feature of model D4 is that, because of its high resistivity, the wind from the innermost part of the disk becomes weaker with time as the magnetic flux in the region diffuses to larger radii. By $t=1800t_{0}$, the region of weakened wind launching extends to a radius $\sim 0.2$ au (Fig. 10(e)). In this region, the rate of angular momentum removal by the wind is significantly reduced, allowing mass to accumulate there. This results in a large plateau in the surface density profile followed immediately by a gap at 0.2 au, with a drop in surface density of one order of magnitude. Model D4 therefore provides another example of gap formation through a wind of varying strength at different radii, in addition to Zone I of the reference simulation. At 2200 inner orbital periods, the mass in the disk within $r<0.2$ au has doubled, and the magnetic flux has dropped well below its initial value. The magnetic field remaining in the dense plateau is too weak to break the surface density into rings and gaps, unlike Zone I of the reference case. The drastically different behaviors of models D_4, D4, and the reference model underscores the key role of the magnetic diffusivity in determining the structure and dynamics of the coupled disk-wind system.

4.2 Magnetic field strength

We will focus first on model beta_3, where the (vertical) magnetic field is $\sqrt{3}$ times stronger than that of the reference case at the disk midplane (i.e., $\beta$ is decreased by a factor of three). One expects a stronger magnetic field to drive a more powerful disk wind, and this is indeed the case. Fig. 9 shows that the wind mass loss rate is typically well above that of the reference case. It is also clear from Fig. 9 that the disk wind in this case is just as variable as, if not more variable than, that in the reference case.

The flow pattern in and around the disk in model beta_3 is more similar to that in the least resistive model D_4 than to the reference case. Specifically, it is dominated by the constant formation and disruption of dense avalanche accretion streams (like that seen in Zone III of the reference case); the quasi-laminar region (Zone II) has all but disappeared. These changes are expected because the magnetic torque exerted by a stronger field removes angular momentum from the disk more efficiently, allowing the accreting streams to develop earlier and at smaller radii. The first prominent stream develops within 150 inner orbital periods (earlier than the reference run but later than model D_4) and terminates near $r=0.2$ au. By 500 inner orbital periods the inner disk is almost completely restructured by the interaction of multiple avalanche streams, as in model D_4, again leading to the formation of a thick envelope where matter moves rather chaotically both inward and outward (see Fig. 8(d)). Compared with model D_4 at similar times, the number of avalanche streams is smaller, the stream-dominated envelope is thinner (typically reaching only $z\sim 0.5R$ compared to $\sim R$), and the wind above the envelope is more massive (compare Fig. 8(e) and (b)). In addition, the expansion of disk material near the midplane is more clearly visible in model beta_3; for example, about half of the disk inside 0.6 au is decreting at the time shown in Fig. 8(e). Nevertheless, the flow patterns in the two cases are fundamentally similar: both are dominated by the avalanche accretion streams, as they are not suppressed by a large enough magnetic diffusivity. We will return to this important point toward the end of the subsection.

The strong, highly variable wind and the stream-dominated envelope create structures in the disk surface density distribution, which can be seen in Fig. 8(f) and Fig. 10(c). There are numerous rings and gaps at the time shown ($t=1650t_{0}$). Their development starts near the inner boundary and spreads outward with time, as the disk material at an increasingly larger radius becomes affected by the magnetic wind launching and the formation of streams. We find an anti-correlation between the surface density and vertical magnetic field strength in the midplane (as in model D_4 and Zone I of the reference case), although the interpretation of the anti-correlation is again complicated by the chaotic motions in the disk and envelope.

The weaker field model, beta3, is dominated by a quasi-laminar region beyond 0.1 au (similar to Zone II of the reference run), and lacks a region with avalanche streams. In fact, not a single stream is formed throughout the duration of the simulation. Within a radius of 0.1 au, a wind appears only intermittently and is less efficient in carrying away angular momentum compared to the wind launched beyond this radius. This allows mass accumulation at small radii, forming the high surface density plateau, as in model D4 (Fig. 10(f)). This plateau is adjacent to a prominent gap at $r=0.1$ au, which separates the variable inner and steady outer wind launching region, similar to the boundary between Zone I and Zone II in the reference case. This difference is further illustrated by the weakest magnetic field model, beta100. In this case, an outflow is still launched, but is unable to escape the simulation domain; instead, it inflates a bubble up to $z\sim 1$ au that is dominated by the toroidal magnetic field with $\beta\approx 5\times 10^{-2}$. The disk remains laminar and mostly unchanged from its initial state.

Taken together, the four models with different field strengths but the same resistivity (beta100, beta3, reference, and beta_3) show a clear trend: as the field becomes stronger, the disk-wind system becomes less laminar, with the avalanche accretion streams becoming increasingly more important. This is the same trend that we see in Section 4.1, where the importance of the avalanche accretion streams increases as the resistivity decreases for a given field strength. The similarity indicates a deep connection between the effects of a stronger magnetic field and a lower resistivity.

The connection can be understood with the help of a second diffusivity parameter $D_{m}\equiv\eta/(hv_{A})$, where $v_{A}$ is the Alfvén speed; it is the inverse of the standard Lundquist number in plasma physics, $S\equiv hv_{A}/\eta$, and is related to the first diffusivity parameter as $D_{m}=\sqrt{\beta/2}\ D$. It measures the importance of the resistive magnetic diffusion on magnetically induced motions. Obviously, a lower resistivity would lead to a better coupling between the magnetic field and the matter. Perhaps less obvious is that a stronger magnetic field would also lead to a better field-matter coupling. This is because a stronger field induces a faster motion over a given length-scale, which leaves less time for the field lines to diffuse resistively relative to the matter. In this sense, the diffusivity parameter $D_{m}$ provides a better indicator for how the coupled disk-wind system behaves: the smaller $D_{m}$ is, the more important the avalanche accretion streams would become. Conversely, the streams are suppressed by a large enough $D_{m}$ (analogous to the suppression of the MRI). This explanation applies not only to models with different parameters, but also across different regions of the reference model. For example, it naturally explains why an avalanche accretion stream develops only at the largest radius (Zone III) of the reference simulation, where the parameter $D_{m}$ is the lowest because it initially decreases with radius as $r^{-1/2}$.

4.3 Disk thickness/temperature

In the simulation t4, the disk aspect ratio $\epsilon$ is increased by a factor of two, thereby increasing the initial disk temperature by a factor of four (see Table 1). To keep $\beta$ and $D_{0}$ (and $D_{m,0}$) the same as the reference run, we increase the initial poloidal field strength $B_{p,0}$ by a factor of two and the resistivity $\eta$ by a factor of four. In this simulation, there is an extended period of time where a stable disk wind is launched within a radius of $r<0.5$ au. Immediately outside of this is an avalanche accretion stream that terminates at $r\approx 0.5$ au. The mass deposited by the stream grows fast enough to form a ring (see Fig. 8(g) and (i)), while the disk region beneath the stream moves radially outward (Fig. 8(h)). This provides the cleanest evidence that a ring can be formed directly from a stream. This ring formation process was also present in the reference model, where it increased the surface density at $r=0.5$ au by a factor of five relative to its neighbouring gap (at a slightly larger radius). The effect is even more noticeable in model t4, with the surface density contrast between the ring and gap increased to a factor of 40. We again see an anti-correlation between the surface density and magnetic flux for the rings and gaps created via an avalanche accretion stream. Naturally, material being deposited at the end of the accretion stream will increase the disk mass locally, but the magnetic flux will not increase here because the mass accretion through the stream is nearly parallel to the field that confines the flow. The magnetic flux is reduced in the region under the stream that has been moving radially outward.

4.4 Resolution

In order to study the effects of the grid resolution on the simulations, we increase the disk resolution by a factor of two in the $\theta$-coordinate direction in model ref-x2grid. The grid now contains 240 uniformly-spaced cells from $\theta=\pi/3$ to $\pi/2$ and 120 non-uniform cells from $\theta=\pi/3$ to the polar axis ($\theta=0$) with a constant ratio for the widths of adjacent cells (see Section 2.3). In order to match the cell aspect ratio with all previous simulations, we also decrease the value of $dr_{0}$ by two ($dr_{0}$ is the width of the first cell on the $r$-grid). This simulation has the same qualitative evolution as the reference run, and we can separate the disk evolution into three zones with radial delineations matching those described for the reference run. The primary difference seen between the reference run and the high resolution run is in the number of rings/gaps formed in Zone I. There are up to three rings at any one time in the high resolution run, but by $t=1650t_{0}$ they have accreted to form one large ring outside the inner boundary at $r=0.03$ au. The remaining zones follow similar behavior to the reference run, indicating that the qualitative behavior of the disk-wind system – especially the disk-wind variability and the ring, gap, and stream formation – are independent of the simulation resolution.

5 Discussion

5.1 Comparison with other work

5.1.1 Steady disk wind

Although one of the focuses of this work is on the variability of winds and their impact on forming radial structures in disks, previous studies have focused largely on steady winds such as those seen in Zone II of the reference simulation. Many authors aim to explore the parameter space of disk properties in an effort to connect them to properties of the wind (e.g., Tzeferacos et al. 2009; Murphy et al. 2010; Sheikhnezami et al. 2012; Stepanovs & Fendt 2014; SF16). In these works, the resistivity often takes the form $\eta\propto hv_{A}$ inside the disk and zero outside. The scaling of this expression is chosen such that the inward accretion and outward diffusion of magnetic flux are able to reach a quasi-steady state, with a fraction of the disk material getting launched into the wind, thereby driving disk evolution through angular momentum removal. Such steady-state solutions, while important for illuminating the physics of wind-launching and for making connection with analytic work, are not guaranteed in nature. There is no apriori reason why the Lundquist number $S\equiv\eta/(hv_{A})$ should be spatially constant or of order unity, especially in the inner disks where the electron number density varies rapidly with radius (caused by collisional ionization of the alkali metals and/or thermionic emission from dust grains; Umebayashi & Nakano 1981; Desch & Turner 2015). We have shown that, in the simplest illustrative case of a constant resistivity, a steady-state is not reached in general.

SF16 have shown that winds can be launched for disks with plasma-$\beta$ from $10^{0.7}$ to $10^{3.5}$. They have also confirmed that there is a critical magnetization below which the launching mechanism transitions from the classical magnetocentrifugal mechanism of Blandford & Payne (1982) to the magnetic tower mechanism driven by the toroidal magnetic pressure gradient (Lynden-Bell, 1996). This transition occurs near $\beta\approx 30$. Our simulations typically have $\beta\approx 10^{3}$ (specifically in Zone II of the reference simulation, $\beta=10^{2.36}$ at $r=0.3$ au after $t=1800t_{0}$), so our finding that the wind is in the magnetic tower regime (see Section 3.3.1) is consistent with previous work. This is also consistent with the 1D analytic model of Bai et al. (2016), who found that the disk wind tends to be driven by a magnetic pressure gradient unless the Alfvén speed near the disk surface is much larger than the local Keplerian speed. Given the small magnetic lever arm in this regime, the mass loading rate is high and the wind speed remains low in the steady part of the wind.

5.1.2 Rings and gaps

In the era of ALMA, rings and gaps are being observed in an increasing number of disks, including, e.g., HL Tau (ALMA Partnership et al., 2015), TW Hya (Andrews et al., 2016), and HD 163269 (Isella et al., 2016). The most commonly invoked explanations for such features are that they are cleared by planets (Dong et al., 2015) or are the result of the condensation of abundant volatile compounds (i.e., snowlines; Zhang et al. 2016). In almost all of the disks simulated in this work, gaps are opened with surface density contrast ratios of at least one order of magnitude. The fact that gaps can be created purely through MHD processes is interesting in its simplicity. The mechanism requires only a magnetized disk with an accretion rate that varies as a function of radius. In our simulations, one way to achieve this variation is through an MHD disk wind that carries away differing amounts of angular momentum from different radial locations in the disk. The magnetic torque acting on the surface of the disk then results in a mass accretion rate of $\dot{M}_{\text{acc}}(R)\approx RB_{\theta}B_{\phi}/\Omega_{K}$ (Bai et al. 2016, see their equation 19). Therefore, a magnetic disk wind can open a gap (or create a ring) so long as $B_{\theta}B_{\phi}$ reaches a local maximum (or minimum) somewhere in the disk. In our particular problem setup, the rings and gaps are formed within a few au of the central mass because only such inner regions have short enough orbital periods to evolve through multiple orbits in the duration of the simulation. However, this mechanism should operate at larger disk radii as long as the magnetic disk wind is able to redistribute angular momentum on these scales as well (Bjerkeli et al., 2016).

Several other MHD mechanisms have been proposed for the creation of radial surface density maxima in disks. For example, the simulations of Flock et al. (2015) and Ruge et al. (2016) have shown that a surface density bump is formed inside the MRI dead zone of a disk, while a gap is opened up outside the dead zone due to MRI-driven mass accretion. The authors find an anti-correlation between the magnetic field and surface density (similar to that shown in Fig. 7) as the magnetic field accumulates in the MRI-active gap outside the dead zone. A related phenomenon are the so-called ‘zonal flows’ (Johansen et al., 2009). Zonal flows occur when a radial pressure gradient is balanced by the Coriolis force, leading to alternating radial bands of sub and super-Keplerian flows (see Armitage 2015 for a recent review). These flows develop in MRI disk simulations from large-scale variations in the Maxwell stress ($B_{r}B_{\phi}$) and result in an anti-correlation between the magnetic pressure and the mass density (Johansen et al., 2009; Kunz & Lesur, 2013; Bai & Stone, 2014). In global disk simulations, zonal flows and zonal (magnetic) fields can be generated by the Hall effect because the additional Hall term in the induction equation can act to radially confine the vertical magnetic flux (Béthune et al., 2016b). The rings and gaps that form in our simulations share some characteristics with the zonal flows (e.g., both show an anti-correlation between surface density and magnetic pressure), but their formation mechanisms are quite distinct.

In addition to variable magnetic winds, our simulations show that rings can also form directly out of rapidly accreting streams. This ring formation mechanism is observed most clearly in the simulation with a larger initial disk thickness (or temperature; model t4, Section 4.3), and to a lesser extent in the reference simulation (Zone III; Section 3.3.2). The streams most likely play a role in forming the rings and gaps observed in several other simulations as well, especially in the lower resistivity (D_4) and stronger field (beta_3) cases, where the disk evolution is dominated by a stream-induced envelope. The concept of avalanche accretion streams itself is not new; Kudoh et al. (1998) found such streams in their 2D simulations of thick AGN accretion disks (where the streams were termed ‘accretion avalanches’; see also Stone & Norman 1994; Matsumoto et al. 1996; Beckwith et al. 2009). What we have shown here is that such avalanche streams tend to develop in regions of low dimensionless magnetic diffusivity $D_{m}$ (or high Lundquist number $S=1/D_{m}$) and that they can deposit material at fixed locations quickly enough to form rings in some cases or generate enough spatial variation in mass accretion to produce multiple rings and gaps in others.

5.2 Magnetic diffusion and wind-dominated vs. stream-dominated accretion

We found two distinct modes of accretion in coupled disk-wind systems initially threaded by a relatively weak, ordered poloidal magnetic field. They are controlled by the dimensionless magnetic diffusivity parameter $D_{m}$ (or the Lundquist number $S=1/D_{m}$). Low values of $D_{m}$ promote the formation of avalanche accretion streams, leading to a ‘stream-dominated’ mode of accretion. Higher values of $D_{m}$ suppress stream formation, leading to a more steady wind that drives accretion in a more laminar disk, i.e., a ‘wind-dominated’ mode of accretion. If $D_{m}$ becomes too large, the strong wind dies out quickly as the wind-launching magnetic field weakens due to diffusion.

The zero diffusivity limit has recently been explored in the 3D ideal MHD simulations of ZS17 that came to our attention near the conclusion of this investigation. Their study has several aspects in common with the work presented here; for example, both concentrate on relatively thin disks (disk thickness $5-10$ per cent of the radius), both adopt a large-scale poloidal magnetic field initially corresponding to a typical plasma-$\beta$ of $10^{3}$ on the disk midplane, and both employ a spherical-polar coordinate system. The main differences are that (1) their simulations are fully 3D, whereas ours are 2D axisymmetric, and (2) their simulations are formally ideal MHD, whereas we included explicit resistivity. Despite these differences, there is qualitative agreement regarding some key features of the simulations. One of the most intriguing features of our 2D simulations is the development of fast accreting avalanche streams near the disk surface. ZS17 also observes the simultaneous fast accretion at high altitudes and slow expansion on the midplane, where the fast accretion occurs in a thick layer above the disk extending up to $\sim 1-1.5$ times the local radius (which they term ‘disk corona’). We believe that their disk corona is physically equivalent to the thick envelope found in our low diffusivity ($D_{m,0}$) cases. In our case, the thick envelope is produced by the repeated formation and disruption of multiple avalanche accretion streams, as can be seen in Fig. 8 (also see movies in the online journal). Although it is not discussed explicitly in their work, we suspect that a similar mechanism is responsible for forming the corona in their simulations as well, although it may be harder to identify rapidly infalling streams in 3D simulations than in 2D ones. The avalanche streams tend to drag the field lines into a radial configuration, which are then wound up by differential rotation. This naturally produces the $r-\phi$ magnetic stress that is found by ZS17 to dominate the coronal accretion. In any case, our simulations add weight to the emerging picture that a chaotic, fast accreting ‘corona’ or envelope plays a crucial role in shaping the structure and evolution of thin disks threaded by relatively weak ($\beta\sim 10^{3}$), open magnetic field lines, especially when the magnetic diffusivity is low.

In more magnetically diffusive disks with larger values of $D_{m}$, the situation can be quite different. On one hand, the explicit diffusivity can enable outward field diffusion by itself, allowing for the possibility of a steady-state balance between the outward diffusion and the inward advection of field lines by mass accretion. On the other hand, it can suppress the MRI in the disk and thus remove the turbulent diffusivity that could also enable such a state. This field advection-diffusion balance is illustrated in Zone II of the reference run (Fig. 5), where the disk connects smoothly to the wind. One might be tempted to identify the subsonic region at the base of the wind as the ‘corona’, which has a similar height ($z\sim R$) in the simulations of ZS17. However, this subsonic wind base is expanding slowly, and the fast accretion (with an effective $\alpha$ parameter of $\sim 2$) is limited to the high density disk proper (see Fig. 5(c)). Most of the accretion in the more diffusive cases is through the disk and is driven by the magnetic torque from the wind; this is in contrast to the ideal MHD simulations of ZS17, where the wind plays only a minor role in driving the accretion. Wind-driven accretion appears to hold for other 2D non-ideal MHD simulations of coupled disk-wind systems that include either resistivity (e.g., SF16) or ambipolar diffusion and the Hall effect (Bai & Stone, 2017). Whether it stays true in 3D or not remains to be determined. In any case, our results, together with those of ZS17 and others, suggest that explicit magnetic diffusion from non-ideal MHD effects plays a key role in determining the extent to which the accretion is dominated by a wind or by a thick ‘corona’/stream-induced envelope.

We speculate that the two modes of accretion may be unified in the following sense. For a thin disk magnetized with a relatively weak poloidal magnetic field, a toroidal magnetic field is quickly generated. This would naturally ‘puff up’ the disk in the vertical direction, creating a thick envelope that is supported by the magnetic pressure gradient vertically (see also Hirose & Turner 2011) and by rotation in the cylindrical $R$-direction. This rotationally supported envelope can become unstable to the MRI, just as the disk, especially when the dimensionless magnetic diffusivity $D_{m}$ is low. Indeed, its larger vertical extent makes the envelope more prone to the MRI than the thin disk itself. When the diffusivity $D_{m}$ is high, the MRI is suppressed in the envelope, which allows the magnetic pressure gradient-driven expansion that produced the envelope to continue smoothly to larger distances, forming a slow wind.

5.3 Implications of ring and gap formation on dust dynamics and planet formation

It is well-understood that an outward radial pressure gradient in a disk causes the gas to rotate at a sub-Keplerian velocity. Consequently, grains that orbit at the Keplerian speed would experience a headwind, lose angular momentum, and migrate inward (Weidenschilling, 1977). If, however, the radial pressure gradient is reversed due to a ‘bump’ in the disk surface density, this local maximum is able to halt the inward radial drift of particles and trap them there. In 3D, these radial bumps in the disk surface density can lead to the growth of the Rossby wave instability (RWI), whereby a high pressure nonaxisymmetric vortex grows exponentially, creating an azimuthal dust trap (Lovelace et al., 1999). Azimuthally asymmetric features are observed in disks and can be explained by azimuthal variations in the gas-to-dust ratio and grain size segregation (Casassus et al., 2013). The formation of radial and azimuthal pressure traps are critical for stopping solid particles from spiraling into the star on rapid timescales, thereby allowing the onset of grain growth and the eventual formation of planetesimals.

The emergence of multiple radial density maxima in the simulations presented here shows the potential importance that MHD disk winds and avalanche streams can have on growing the seeds of planet formation. This is especially true in the inner ($\lesssim 1$ au) regions of the disks, where jets and winds of young stellar objects are most likely launched and where the most common type of planets (super-Earths/mini-Neptunes) are located. Indeed, Sun-like stars have a 50 per cent chance of having a compact system of small (non-giant) planets within 1 au (Winn & Fabrycky, 2015). While the importance of migration for this population of planets is still unknown, it seems plausible that many of these planets may have formed at their current locations. If this is the case, an efficient mechanism is needed for trapping dust grains on sub-au scales, especially during the relatively early phases of star formation when most of the mass is processed through the inner disk. Rings and gaps produced by variable magnetic winds and avalanche streams that drive rapid accretion potentially provide an opportunity for grains to accumulate and grow early in the disk’s life. This mechanism will likely not be as efficient at later times when YSO jets and winds are observed to be less powerful and the accretion rate is lower. More work is needed to quantify the effects that the disk structures created by magnetically driven winds and avalanche streams have on the dynamics and growth of dust grains.

5.4 Future directions

The initial set of idealized simulations presented in this work is aimed at illustrating the basic processes through which rings and gaps are produced by variable winds and avalanche accretion streams. In future refinements we will extend the simulations to a less restrictive geometry. Preliminary 2D (axisymmetric) simulations with two hemispheres show qualitatively similar results to those presented here. Full 3D simulations, such as those of ZS17 but with non-ideal MHD effects (see Suzuki & Inutsuka 2014; Gressel et al. 2015; Béthune et al. 2016a), will be needed to determine whether the formation of rings and gaps is artificially enhanced by the axisymmetric geometry and whether the rings formed in 2D simulations are unstable to the Rossby wave instability. Also, 3D simulations are needed to determine whether the accretion streams will survive in 3D since the streams are physically related to the MRI channel flows, which are known to behave differently in 2D and 3D (Balbus & Hawley, 1998; Latter et al., 2009, 2015). Just as important, there is a need to include a detailed calculation of the ionization structure of the disk and the non-ideal MHD effects, especially the Hall effect, which is expected to be important near the disk midplane on the au scale (Balbus & Terquem, 2001; Kunz & Balbus, 2004). Post-processing the simulations through radiative transfer calculations is needed to make predictions for observable quantities (e.g., variability in near infrared emission and profiles of optical forbidden lines) that can be compared directly with observations (e.g., Cody et al. 2014; Simon et al. 2016). The long term goal is to include a self-consistent treatment of dust grains in the full 3D non-ideal MHD coupled disk-wind simulations. This would allow us to quantify the effects that wind and stream-induced disk structures have on grain trapping and growth, as well as any potential back-reaction that the grains might have on the gas dynamics of disk accretion, stream formation, and wind launching. The rings and gaps found in our starting simulations provide added impetus to explore the interplay between dust dynamics and winds and streams in greater depth.

6 Conclusions

We have presented the results of 2D (axisymmetric) resistive MHD simulations of coupled disk-wind systems with a range of disk parameters (resistivity, magnetic field strength, and temperature), focusing on geometrically thin disks. We find that the structure and dynamics of the disk-wind system strongly depend on the dimensionless magnetic diffusivity parameter $D_{m}\equiv\eta/(hv_{A})$ and that interesting disk features, including rings and gaps, are naturally produced. Specifically, we find that:

There are two distinct modes of accretion depending on the dimensionless parameter $D_{m}$. Disks with low values of $D_{m}$, from either a small resistivity or high field strength, tend to develop fast ‘avalanche accretion streams’. The rapid formation and disruption of such streams often leads to a clumpy, thick envelope above the disk that dominates the dynamics of the system (e.g., models beta_3 and D_4; see Fig. 8), although a highly variable wind is still launched above the envelope. This envelope appears equivalent to the thick disk ‘corona’ found independently by ZS17 in their 3D ideal MHD simulations. In both cases, the disk below the corona/envelope is often expanding radially outward. The streams (and the thick clumpy envelope they produce) are suppressed in simulations with larger values of $D_{m}$ (from either a large resistivity or low field strength, e.g., models beta3 and D4). In these more diffusive (larger $D_{m}$) simulations, most of the accretion occurs through a laminar thin disk rather than the thick clumpy envelope, and the disk accretion is driven mainly by a magnetic wind. 2. 2.

Both wind-dominated and stream-dominated accretion create prominent features in the surface density distribution, especially rings and gaps. The wind-driven ring and gap formation is illustrated most clearly in the innermost region (Zone I) of the reference simulation, where there is substantial redistribution of magnetic flux relative to the mass in the disk that is enabled by the resistivity (Fig. 7). Regions with lower mass-to-flux ratios tend to drive stronger winds and accrete faster, producing gaps; those with higher mass-to-flux ratios tend to accrete more slowly, allowing matter to accumulate and form dense rings. Another clear illustration of wind-driven gap formation can be found in models D4 and beta3, where a strong wind is driven from the outer part of the disk but not from the more magnetically diffusive inner part, creating a deep gap between them (Fig. 10(e) and (f)). The stream-driven ring formation is illustrated most clearly in the thicker disk model (t4), where a stream feeds a prominent ring at a roughly constant radius (Fig. 8, right column). Rings formed this way also have high mass-to-flux ratios. Multiple rings and gaps are formed in other, more complicated cases, especially those with stream-induced envelopes (model D_4 and beta_3). It is likely that both magnetic winds and avalanche accretion streams play a role in their formation, although the relative importance of the two mechanisms is hard to quantify due to the complexity of the flow pattern inside and above the disk. 3. 3.

Powerful winds are launched despite the fact that the magnetic field in the disk is rather weak initially (corresponding to a typical plasma-$\beta\sim 10^{3}$). In the reference simulation where the wind is analyzed in detail, we find that the bulk of the wind is heavily mass-loaded and accelerated by the magnetic pressure gradient to relatively low speeds (a few \times 10~{}\text{km~{}s{}^{-1}}). There are, however, lightly mass-loaded regions that are accelerated magnetocentrifugally to speeds exceeding 100~{}\text{km~{}s{}^{-1}}, comparable to the jet speeds observed in young stellar objects. The magnetic wind can remove angular momentum from the disk efficiently, leading to disk accretion with an effective $\alpha$ parameter up to order unity. Our simulations add weight to the notion of wind-driven disk evolution, especially in the presence of a suitable level of magnetic diffusivity.

Rings and gaps produced in circumstellar disks by magnetic winds and avalanche accretion streams have important implications on the dynamics and growth of dust grains and ultimately planet formation. The local pressure maxima associated with the rings would act to stop the radial drift of solid particles, possibly trapping them long enough to enable enhanced grain growth that facilitates planetesimal formation. This may be especially important in the inner (i.e., few tenths of an au) disk regions where the largest population of planets reside, as seen by Kepler. The simulations presented in this work lay the foundation for future explorations of these and other aspects of the coupled disk-wind systems.

Acknowledgements

We thank Dom Pesce, Zhaohuan Zhu, and William Béthune for useful comments. This work is supported in part by NSF AST-1313083 and NASA NNX14AB38G. The simulations in this work were carried out on the Rivanna computer cluster at the University of Virginia.