Wind-driven Accretion in PPDs

Wind-driven Accretion in Protoplanetary Disks — I: Suppression of the Magnetorotational Instability and Launching of the Magnetocentrifugal Wind

Xue-Ning Bai11affiliation: Department of Astrophysical Sciences, Peyton Hall, Princeton University, Princeton, NJ 08544 22affiliation: Institute for Theory and Computation, Harvard-Smithsonian Center for Astrophysics, 60 Garden St., MS-51, Cambridge, MA 02138 33affiliation: Hubble Fellow & James M. Stone11affiliation: Department of Astrophysical Sciences, Peyton Hall, Princeton University, Princeton, NJ 08544

We perform local, vertically stratified shearing-box MHD simulations of protoplanetary disks (PPDs) at a fiducial radius of 1 AU that take into account the effects of both Ohmic resistivity and ambipolar diffusion (AD). The magnetic diffusion coefficients are evaluated self-consistently from a look-up table based on equilibrium chemistry. We first show that the inclusion of AD dramatically changes the conventional picture of layered accretion. Without net vertical magnetic field, the system evolves into a toroidal field dominated configuration with extremely weak turbulence in the far-UV ionization layer that is far too inefficient to drive rapid accretion. In the presence of a weak net vertical field (plasma at midplane), we find that the magnetorotational instability (MRI) is completely suppressed, resulting in a fully laminar flow throughout the vertical extent of the disk. A strong magnetocentrifugal wind is launched that efficiently carries away disk angular momentum and easily accounts for the observed accretion rate in PPDs. Moreover, under a physical disk wind geometry, all the accretion flow proceeds through a strong current layer with thickness of that is offset from disk midplane with radial velocity of up to times the sound speed. Both Ohmic resistivity and AD are essential for the suppression of the MRI and wind launching. The efficiency of wind transport increases with increasing net vertical magnetic flux and the penetration depth of the FUV ionization. Our laminar wind solution has important implications on planet formation and global evolution of PPDs.

Subject headings:
accretion, accretion disks — instabilities — magnetohydrodynamics — methods: numerical — planetary systems: protoplanetary disks — turbulence

1. Introduction

Protoplanetary disks (PPDs) are gaseous disks surrounding protostars. The gas in PPDs are found to be rapidly accreting to the protostar with accretion rate of yr, with typical disk lifetime of about 1-10 Myrs (e.g., Hartmann et al., 1998; Sicilia-Aguilar et al., 2006). Despite large number of observational programs aiming at revealing the structure, composition and evolution of PPDs (see Williams & Cieza, 2011 and references therein), two crucial theoretical questions on the gas dynamics of PPDs remain poorly understood: What is the level of turbulence in PPDs? How efficient is angular momentum transport in PPDs? The answer to these questions are essential to understanding the structure and evolution of the PPDs, as well as a series of processes in planet formation. In particular, grain growth (e.g., Birnstiel et al., 2010; Zsom et al., 2010), transport of solids (e.g., Garaud, 2007; Hughes & Armitage, 2010) are both sensitive to the radial structure of PPDs and level of turbulence. Current models for planetesimal formation such as the streaming instability (Johansen et al., 2009; Bai & Stone, 2010a, b) and gravitational instability (Youdin, 2011) generally favor weak turbulence and small radial pressure gradient. Moreover, turbulent mixing and disk winds have a significant influence on the disk chemistry (Semenov et al., 2006; Heinzeller et al., 2011). When planets have formed, planet migration via planet-disk interaction also depends on the disk radial profile and diffusion processes (e.g., Paardekooper et al., 2010; Baruteau et al., 2011).

1.1. The Current Understanding of Accretion in PPDs

We are mainly interested in the T-Tauri (class II) phase of PPDs when the envelope infall has ended and the entire disk is visible. At this stage the disk is in general not massive enough for gravitational instability to take place (Zhu et al., 2010). In this paper, we focus on magnetic mechanisms.

Most regions of PPDs111The region that is close to the inner edge of PPDs is sufficiently hot (K) due to direct illumination by the star that thermal ionization of Alkali species Na and K will provide sufficient ionization and the gas behave as ideal MHD (Umebayashi & Nakano, 1988), which is not the concern of this paper. are very weakly ionized (Hayashi, 1981; Igea & Glassgold, 1999), hence the gas dynamics is strongly affected by non-ideal magnetohydrodynamics (MHD) effects due to the finite gas conductivity, which include Ohmic resistivity, Hall effect and ambipolar diffusion (AD). All three effects are relevant and important in PPDs (Wardle, 2007; Bai, 2011a). Generally speaking, Ohmic resistivity dominates in dense regions with weak magnetic field (e.g., midplane in the inner region of PPDs), AD dominates in tenuous regions with strong magnetic field (e.g., disk surface and outer region of PPDs), while the Hall regime lies in between.

It is widely believed that PPDs are turbulent as a result of the magnetorotational instability (MRI, Balbus & Hawley, 1998). The MRI turbulence transports angular momentum radially within the disk that allows the majority of the materials to be accreted onto the protostar while a small fraction of mass disperses away. Non-ideal MHD effects in PPDs strongly modify the behavior of the MRI. Currently, most studies of the MRI in PPDs take into account only the effect of Ohmic resistivity, and it is found that in the inner region of PPDs (about 0.5-5 AU), the disk midplane is too weakly ionized for the MRI to operate (i.e., the dead zone, Gammie, 1996), while the disk surface is still prone to the MRI and should be turbulent (i.e., the active layer). A large number of numerical simulations have been conducted either in the local shearing-box framework or using full global approach to study and characterize the gas dynamics of the active layer and the dead zone, as well as exploring their physical consequences (e.g., Fleming & Stone, 2003; Turner et al., 2007; Turner & Sano, 2008; Ilgner & Nelson, 2008; Oishi & Mac Low, 2009; Dzyurkevich et al., 2010; Hirose & Turner, 2011; Flaig et al., 2012; Flock et al., 2012). These studies show that the boundary between the active layer and the dead zone is characterized by the Ohmic Elsasser number (see equation (6) ). Moreover, the dead zone is not completely “dead” in the sense that sound waves injected from the active layer bounce back and forth and give rise to some small Reynolds stress. Nevertheless, angular momentum transport is largely dominated by the MRI turbulence in the active layer, and the strength of the turbulence (considering Ohmic resistivity only) appears to be able to drive rapid accretion consistent with observations.

Comparatively, the effects of non-ideal MHD terms other than Ohmic resistivity on the MRI in PPDs are less well understood. They are extensively studied in the linear regime without vertical stratification (see Wardle, 1999; Balbus & Terquem, 2001; Wardle & Salmeron, 2012 for the Hall effect, Blaes & Balbus, 1994; Kunz & Balbus, 2004; Desch, 2004 for the effect of AD, and Pandey & Wardle, 2012 for a general study), and with vertical stratification (Salmeron & Wardle, 2005, 2008). In the non-linear regime, so far all numerical simulations (all using the local shearing-box approach) focus on individual non-ideal MHD effects such as MRI with Ohmic and Hall terms (Sano & Stone, 2002a, b) and MRI with AD (Hawley & Stone, 1998; Bai & Stone, 2011). Most of the simulations are vertically unstratified. These simulations provide useful criteria for the MRI to be self-sustained in the non-linear regime and the results were applied in the framework developed by Bai (2011a) to estimate the efficiency of MRI-driven angular momentum transport in PPDs. It was shown that even in the most optimal scenario, the MRI-driven accretion rate falls below typical observed rate by about an order of magnitude at the inner disk around 1 AU. The main reason is that the strength of the MRI turbulence is expected to be dramatically reduced in the conventional “active layer” of the disk once AD is taken into account. Similar conclusions were also drawn from Perez-Becker & Chiang (2011a, b), and from Mohanty et al. (2013); Dzyurkevich et al. (2013) for different stellar masses. Relatively large accretion rate can be achieved in the outer region of the disk under optimal magnetic field geometry, and with the assistance of tiny grains (Bai, 2011b).

An alternative scenario for describing the gas dynamics in PPDs is the picture of magnetocentrifugal wind (Blandford & Payne, 1982; Pudritz & Norman, 1983): outflowing gas from accretion disks can be centrifugally accelerated along magnetic field lines when the inclination angle of the poloidal field is above (relative to the disk normal)222The X-wind model (Shu et al., 1994) is also magnetocentrifugal in nature, with the wind launched near the inner edge of the disk. We are interested in the wind launched from radially extended region in PPDs.. The magnetocentrifugal wind scenario has also been extensively explored with global simulations. Early global simulations treat the disk as a boundary condition (i.e., razor-thin) with axisymmetry, and prescribe the rate of outflow from the disk (Ouyed & Pudritz, 1997; Krasnopolsky et al., 1999, 2003). These simulations demonstrated the robustness of the magnetocentrifugal acceleration and collimation, and further found that the flow structure is sensitive to the prescribed rate of mass loading from the disk, and may lead to episodic formation of jets (Ouyed et al., 1997; Krasnopolsky et al., 1999; Anderson et al., 2005). More recent simulations (most of which are two-dimensional) that do resolve the disk generally rely on artificially prescribed and excessively large diffusion (that is unjustified) to prevent rapid magnetic flux accumulation near the central object as mass accretes (Kato et al., 2002; Casse & Keppens, 2002, 2004), and the resulting wind properties largely depend on the prescribed resistivity profile in the disk (Zanni et al., 2007).

In reality, the wind launching process is governed by the microphysics within the disk, with mass loading rate determined by requiring that the flow smoothly passes the slow magnetosonic point (Wardle & Koenigl, 1993; Li, 1995; Ogilvie & Livio, 2001; Ogilvie, 2012). It was found that launching a laminar disk wind generally requires a relatively strong vertical background magnetic field with around equipartition strength at disk midplane (Wardle & Koenigl, 1993; Ferreira & Pelletier, 1995), though weaker field is possible in the presence of strong ambipolar diffusion (Li, 1996). The vertical magnetic field can not be too strong, which would make it difficult to be bent by the disk material, and would result in substantial sub-Keplerian rotation that hinders wind-launching (Shu et al., 2008; Ogilvie, 2012). More detailed study of the wind-launching criteria and representative solutions were presented in Königl et al. (2010) and Salmeron et al. (2011) where all the non-ideal MHD effects were taken into account. It appears that the MRI and the laminar wind scenarios tend to mutually exclude each other (Salmeron et al., 2007): MRI operates when the vertical magnetic field strength is well below equipartition strength at the disk midplane, while launching an magnetocentrifugal wind requires the vertical field strength to be around equipartition. The strong magnetic field in the magnetocentrifugal wind scenario would drive very efficient accretion, making the disk much more tenuous (with very small surface density) than that in the MRI scenario for a standard accretion disk (Combet & Ferreira, 2008).

To briefly summarize, MRI generally requires the background vertical field to be weak, and has difficulty in accounting for the rapid accretion rate in PPDs, especially in the inner region of AU. The alternative picture of magnetocentrifugal wind generally requires the presence of strong background vertical field of equipartition strength that either drives accretion that is too rapid or results in a tenuous disk, but the microphysics of the wind launching process still requires more realistic treatment.

1.2. This Work

We conduct vertically stratified shearing-box simulations of a local patch of a PPD that include a realistic profile of both Ohmic resistivity and AD coefficients. The diffusion coefficients are interpolated from a pre-computed look-up table in real simulation time based on the gas density, temperature (fixed) and ionization rate (calculated from the density profile). We have not included the Hall effect (which is numerically more challenging and demanding), and this is the first step towards understanding the non-ideal MHD of PPDs beyond Ohmic resistivity. This is the first time that Ohmic resistivity and AD are simultaneously included in vertically stratified simulations to study the gas dynamics of PPDs that incorporates the disk microphysics in the most realistic manner.

All our simulations except one include a weak vertical magnetic field, which is likely to be more realistic for a local patch of a PPD. The other reason for including a vertical field is that such a field geometry is more favorable for the MRI to operate in the presence of AD (Bai & Stone, 2011). We note that most vertically stratified shearing-box simulations to date adopt a zero net vertical magnetic flux field geometry (e.g., Stone et al., 1996; Shi et al., 2010; Davis et al., 2010) and we will show that such field geometry would make the strength of the MRI turbulence diminishingly small due to AD. Including net vertical magnetic flux in shearing-box simulations places strong demands on the robustness of numerical algorithms, especially in the magnetic dominated disk corona (Miller & Stone, 2000). Recently such simulations are successfully performed by several groups, mostly in the ideal MHD regime, and it was shown that the inclusion of net vertical magnetic field always leads to an outflow from an MRI-turbulent disk (Suzuki & Inutsuka, 2009; Fromang et al., 2012; Bai & Stone, 2013). In the context of PPDs, it was found that including a net vertical magnetic flux does not change the basic picture of layered accretion (again, considering Ohmic resistivity only), but stronger net vertical flux leads to stronger MRI turbulence in the active layer and reduces the vertical extent of the dead zone, as well as stronger outflow (Okuzumi & Hirose, 2011). Although these results suggest the simultaneous existence of the MRI and magnetocentrifugal wind, Bai & Stone (2013) pointed that the outflow from MRI active region is unlikely to be directly connected to a global magnetocentrifugal wind due to the MRI dynamo and symmetry considerations. Our simulations in this paper will further address the potential connection of the disk outflow to an magnetocentrifugal wind in the context of PPDs, and we arrive at positive conclusions, in contrast with the ideal MHD case studied by Bai & Stone (2013).

More specifically, we consider a fiducial model that corresponds to a minimum-mass solar nebular at 1 AU, assuming solar abundance of chemical composition and (in mass) of m sized grains for the chemistry calculation. We find that although the initial condition is unstable to the MRI (with net vertical magnetic field much weaker than equipartition), the disk rapidly adjusts to a new laminar configuration that is stable to the MRI. The new laminar state is characterized by an outflow launched by the magnetocentrifugal mechanism, and the outflow can achieve a physical wind geometry (poloidal streamlines at the top and bottom sides of the disk bend towards the same radial direction) by having a thin current layer where the horizontal field flips. The magnetocentrifugal wind launched in this scenario can efficiently transport angular momentum to account for the observed PPD accretion rate while without being too efficient to deplete the disk (as in the conventional wind scenario). For clarity, we focus on the physical properties of the new laminar solution and together with a parameter study all at 1 AU in this paper, extension of the results to other disk radii will be presented in a companion paper.

This paper is organized as follows. In Section 2, we describe the numerical method and the overall setup of our simulations. Three contrasting simulations are presented in Section 3, highlighting the new laminar wind solution. In Section 4, we study the physical properties of the new laminar wind solution in detail and discuss the wind launching mechanism. In Section 5, we conduct a thorough parameter study to further explore the properties of the laminar wind. We discuss the robustness and implications of our wind solution in Section 6, together with the conclusion.

2. Simulation Setup

2.1. Formulation

We use the Athena MHD code (Stone et al., 2008) to study the gas dynamics in PPDs and perform three-dimensional numerical simulations. We consider a local patch of a PPD and adopt the conventional shearing-box approach (Goldreich & Lynden-Bell, 1965). MHD equations are written in a Cartesian coordinate system in the corotating frame at a fiducial radius with Keplerian frequency :


where is the total stress tensor


is the identity tensor, , , are the gas density, velocity and pressure respectively, is the magnetic field, are unit vectors pointing to the radial, azimuthal and vertical directions respectively, where is along the direction. Note that the equations are written in units such that magnetic permeability is 1. Vertical gravity is included to account for density stratification. Periodic boundary conditions are used in the azimuthal direction, while the radial boundary conditions are shearing periodic as usual. Shearing-box source terms (Coriolis force and tidal gravities) have been readily implemented in Athena (Stone & Gardiner, 2010), which uses an orbital advection scheme that splits the system into an advective part for the background shear flow , and a fluctuating part with velocity fluctuation :


The advection scheme not only accelerates the calculation by permitting larger time steps, but also improves the accuracy by making the truncation error independent of radial location.

We use an isothermal equation of state , where is the isothermal sound speed. In reality, multiple radiative processes such irradiation and scattering are involved in PPDs, making the disk surface layer substantially hotter (e.g., Hirose & Turner, 2011). The hotter disk surface affects the properties disk wind, however, the main problems addressed in this paper are largely magnetic: the suppression/survival of the MRI and wind launching are largely controlled by non-ideal MHD effects, where thermodynamics only plays a minor role.

Being weakly ionized, PPDs are not perfectly conducting, which is reflected in the non-ideal MHD terms in the induction equation (e.g., Wardle, 2007; Bai, 2011a)


where is the current density, denotes unit vector along , subscript “” denotes the vector component that is perpendicular to , and are the Ohmic, Hall and the ambipolar diffusivities. Note that represents the velocity for the bulk of the gas (i.e., neutrals), while tracer amount of charged species provides conductivity and gives rise to non-ideal MHD effects.

The magnetic diffusivities depend on the number density of the charged species, and are characterized by the dimensionless Elsasser numbers, defined as


for Ohmic, Hall and AD respectively, where is the Alfvén velocity. In the absence of grains, we have , and being independent of (see next subsection for details). Generally speaking, self-sustained MRI turbulence requires these Elsasser numbers to be greater than 1.

Ohmic resistivity and ambipolar diffusion (AD) have been implemented in Athena (Davis et al., 2010; Bai & Stone, 2011). Furthermore, we have implemented super time-stepping that substantially accelerates the calculation, which we have described in detail in Simon et al. (2013). Although our simulations do not include the Hall term, we do evaluate and assess its importance in our analysis.

We use natural unit in our simulations, where , . The initial density profile is taken to be Gaussian


with being the midplane gas density in code unit, and is the thermal scale height. We perform simulations with both zero and non-zero net-vertical magnetic flux. For simulations with zero net vertical flux, the initial field is purely vertical given by , where is the radial size of the simulation box, and parameterized by the midplane plasma . For simulations with net-vertical flux, we add a uniform vertical field (parameterized by midplane plasma ) on top of the sinusoidally varying component.

Physically, we consider a patch of the PPD using the minimum-mass solar nebular (MMSN, Weidenschilling, 1977; Hayashi, 1981) disk model, with the surface density given by g cm, and temperature given by K, where is the radius to the central star (whose mass is fixed at ) measured in astronomical unit (AU). The mean molecular weight of the neutrals is taken to be from which the sound speed km s (at 1 AU) hence other quantities can be easily evaluated. In particular, in our code unit corresponds to field strength of Gauss. These physical scales are needed to normalize the magnetic diffusivities (next subsection) to code unit.

We use outflow boundary condition in the vertical direction which copies the density, velocity and magnetic fields in the boundary cells to the ghost zones, with the density attenuated following the Gaussian profile to account for vertical gravity. In the case of mass inflow, the vertical velocity is set to zero at the ghost zones. A density floor of (in code unit) is applied to avoid numerical difficulties at magnetic dominated (low plasma ) regions. We have checked that horizontally averaged densities in the saturated states of all our simulations are always well above the density floor333Except for Run OA-nb-F3 where a density floor of is applied, which will be discussed in Section 3.2.. Moreover, the use of outflow boundary condition no longer conserves mass in the simulation box. We compensate for the mass loss so that a steady state can be achieved, following the same procedure as Bai & Stone (2013): the density of each cell is modified by the same proportion at the end of each time step so that the total mass in the simulation box remains the same. For most of our runs, the mass change over the duration of our simulations (if mass conservation were not enforced) is only a tiny fraction of the total mass.

2.2. Calculation of Magnetic Diffusivities

The magnetic diffusivities are evaluated based on the chemistry calculation. Instead of evolving the chemical network in real time, as done by a number of previous works (Turner et al., 2007; Turner & Sano, 2008; Ilgner & Nelson, 2008), we assume equilibrium chemistry (similar to Hirose & Turner, 2011), because the recombination time has been shown to be much shorter than the dynamical time scale (Bai, 2011a). We adopt a complex chemical reaction network (Ilgner & Nelson, 2006; Bai & Goodman, 2009; Bai, 2011a) that is based on the UMIST database (Woodall et al., 2007). In our fiducial model considered in this paper, we fix the elemental composition to be solar, with well-mixed m grains whose abundance is in mass (about 0.01 solar, corresponding to substantial grain growth and settling). We follow the same procedure and methodology described Section 3.4 in Bai (2011a) to evolve the network for years, with further details provided in Sections 3.2-3.5 of Bai & Goodman (2009). Given the chemical composition, variable parameters of the network include gas density , gas temperature , and the ionization rate , which are scanned to give a complete coverage of the parameter space relevant to PPDs. The outcome of the scan is a look-up table of magnetic diffusivities that is read into the code so that , and at each grid cell can be evaluated by interpolation in real simulation time. Since we adopt an isothermal equation of state, is fixed, our look-up table is essentially two-dimensional ( and ).

The ionization rate in the disk depends on the column density to the disk surface. We calculate the horizontally averaged vertical density profile in real simulation time, from which a column density profile can be reconstructed, an approach similar to previous works (e.g. Turner et al., 2007). The sources of ionization include radioactive decay, cosmic ray and stellar X-ray ionizations, with prescriptions given in Section 3.2 of Bai (2011a), with fixed X-ray luminosity of erg s and X-ray temperature of keV. In addition, we consider the effect of far-ultraviolet (FUV) ionization. According to Perez-Becker & Chiang (2011b), FUV photons almost completely ionize tracer species such as C and S and give ionization fraction of the order with penetration depth g cm depending on the effectiveness of dust attenuation and self-shielding 444There are large uncertainties associated with the FUV ionization. For example, it has been noted that the FUV photons may be shielded if a dusty wind is launched from the inner disk near the dust sublimation radius (Bans & Königl, 2012). Moreover, photochemistry in the FUV layer plays an important role in determining the molecular composition and and ion abundances (Walsh et al., 2012), and particularly, the FUV penetration depth may be overestimated in Perez-Becker & Chiang (2011b) due to the conversion of C and S into molecular ions which facilitates recombination (Ádámkovics et al., 2011, Al Glassgold, private communication). Our treatment should be considered as a first approximation that captures basic physics rather than deep into the details.. For simplicity, we assume an ionization fraction of in the form of carbon in the FUV layer, whose column density is chosen by default as g cm. Within the FUV layer, the magnetic diffusivities expressed in the form of Elsasser numbers under the MMSN disk model are found to be


where is ratio of magnetic pressure to midplane gas pressure. A smooth transition (across about 4 grid cells) on the magnetic diffusivities from the FUV ionization layer at to the X-ray/cosmic-ray dominated ionization layer (based on our chemistry calculations) at is applied. As Ohmic resistivity plays essentially no role in the low density region of the FUV layer, we simply do not reset Ohmic resistivity in this layer.

We calculate the magnetic diffusivities from the number density of all charged species at the end of the chemical evolution following Section 2 and 3.5 of Bai (2011a). Note that the Ohmic resistivity is independent of magnetic field strength , while Hall and ambipolar diffusivities do depend on . In the grain-free case, we have , and . In this case, we can simply fit the proportional coefficients, and respectively, and put them into the look-up table. In the presence of small grains, a situation studied in detail in Bai (2011b), () is proportional to () when is sufficiently weak or sufficiently strong respectively, while is roughly proportional to () at some intermediate field strength. In this case, we include in the look-up table the two proportional coefficients () at weak and strong field regimes from the fitting respectively, together with a transition field strength so that


By comparing with Figure 1 of Bai (2011b), we see that the transition field strength corresponds to the situation that the ion gyro-frequency equals its collision frequency with the neutrals (or the ion Hall parameter equals one), which can be directly calculated given the gas density.

Being diffusive processes, large Ohmic resistivity near the disk midplane and strong AD in the tenuous disk corona would significantly limit the code efficiency even super time-stepping is used to accelerate the calculation. In the Ohmic regime, the diffusive time step scales as where is the minimum grid spacing. In practice, we add a cap to the Ohmic resistivity so that in code unit . Since the magnetic field strength near the disk midplane never reaches equipartition in all our simulations, the Elsasser number at the disk midplane is always smaller than , well below the threshold value . Also, this cap value of makes the diffusion time scale much smaller than the dynamical time scale, hence captures the basic effect of strong diffusion at disk midplane even if resistivity is much higher in reality. Note that our resistivity cap is much larger than most previous works (thanks to the use of super time-stepping), where the cap was of the order in natural unit (e.g., Fleming & Stone, 2003; Okuzumi & Hirose, 2011). In the AD regime, the time step scales as ( is the local ratio of gas to magnetic pressure). In the same spirit, we apply a floor to so that in every grid cell. This floor value of is again sufficiently small so that it does not make the otherwise stable field configuration unstable to the MRI (see Figure 16 of Bai & Stone, 2011), and it retains the effect of strong diffusion.

3. Fiducial Simulations and Results

In this section, we present three benchmark simulations with very similar initial setup but evolve into dramatically different states. All three simulations are three-dimensional (3D) shearing-box with vertical stratifications, located at 1 AU in a MMSN disk, adopting a chemistry model with well-mixed m grains with abundance of ( solar), and a FUV column density g cm. All simulations are run for about 150 orbits ().

In the first simulation (Run O-b5), only Ohmic resistivity is included, which aims at modeling the conventional picture of layered accretion. We have also included a net vertical magnetic flux of . A sinusoidally varying vertical field of is also added as initial condition to avoid the strong channel flows (but has no effect on the saturated state of the system). The simulation box size is in the radial (), azimuthal () and vertical () dimensions respectively, with a computational grid of cells. Our simulations have relatively high resolution in and (24 cells per , or 34 cells if one defines the scale height to be , as in a number of works) to properly resolve the MRI turbulence (Davis et al., 2010; Sorathia et al., 2012; Bai & Stone, 2011), and have relatively large simulation box to capture the mesoscale structures of the MRI turbulence (Simon et al., 2012).

In the second simulation (Run OA-nobz), both Ohmic resistivity and AD are included, and we adopt a zero net vertical magnetic flux configuration with . Again, the choice of has no effect on the saturated state of the system. After experimenting with several initial setups, we find that due to the extremely weak level of turbulence, the gas near the disk surface has very limited magnetic support hence drops very rapidly, and we must reduce the box height to (and cells) so that the gas density at vertical boundaries is not too small which would severely reduce the numerical time step. Furthermore, we started the simulation with density floor , but then found that the density at vertical boundaries is always at the level of the density floor. We thereby keep reducing the density floor until when we no longer find artificial features near the vertical boundaries.

Finally, we conduct a contrasting simulation from the above two cases (Run OA-b5). The initial setup of the simulation is exactly the same as in the first simulation (Run O-b5) except that AD is included.

Figure 1 illustrates the initial profile of the Ohmic, Hall and ambipolar Elsasser numbers for our runs O-b5 and OA-b5 (similar but not exactly the same for run OA-nobz due to the different magnetic configuration). Also shown is the initial profile of plasma . From the disk midplane to surface, the dominant non-ideal MHD effects are Ohmic resistivity, Hall effect and AD respectively for the initial field configuration, and MRI unstable region is located at around where all Elsasser numbers are greater than 1 and plasma is well above 1. The initial FUV ionization front is located at about . The density floor of is applied to regions beyond , hence the and curves flattens out (this artifact will disappear as the system evolves).

The most common accretion diagnostics is the component of the Reynolds and Maxwell stresses, which measures the local rate of radial transport of angular momentum


where the over bar indicates horizontal averaging. The total rate of radial angular momentum transport is characterized by the parameter (Shakura & Sunyaev, 1973)


In Figure 2 we show the time evolution of the horizontally averaged Maxwell stress for the three fiducial simulations. Albeit for similar initial setup, the three runs show very distinctive characteristics as the systems evolve into saturated/steady states. The results from the three simulations will be analyzed in detail in the following three subsections.

For clarity, we further provide a list of all our fiducial simulation runs and their parameters in Table 1. In particular, we introduce the letter “S” for runs with very small horizontal domain, where the simulations are essential one-dimensional, and letter “E” for runs with enforced even- symmetry (see Section 4). We use the label “bn” to denote plasma for the vertical background field. All other simulations are run for about 200 orbits ().

Figure 1.— Initial Elsasser number profile for Ohmic resistivity (red solid), Hall diffusivity (green dashed) and ambipolar diffusivity (blue dash-dotted) for our fiducial runs O-b5 and OA-b5 . Note the cap on and around the disk midplane and on beyond about . Also shown is the initial profile for plasma (black solid). Note that we show the Hall Elsasser number () while it is not included in the simulations.
Run Diffusion Box Size Section
O-b5 Ohm 3.1
OA-nobz Ohm, AD 3.2
OA-b5 Ohm,AD 3.3, 4
S-OA-b5 Ohm,AD 4.4
E-OA-b5 Ohm,AD 4.4.1
S-OA-b5-12H Ohm,AD 4.5
S-OA-b5-20H Ohm,AD 4.5
S-OA-b5-24H Ohm,AD 4.5

MMSN disk model at 1AU, X-ray luminosity of ergs s and temperature keV, well mixed m grains with abundance of , penetration depth of g cm for the FUV ionization are assumed for all runs.

Table 1Summary of All Fiducial Simulations.
Figure 2.— The space-time plot for the (horizontally averaged) vertical profiles of the Maxwell stress in the three fiducial runs O-b5 (top), OA-nobz (middle) and OA-b5 (bottom). The colors are in logarithmic scales (in code unit). White contours in the upper panel correspond to vertical Elsasser number . The discontinuous transitions in the middle panel correspond to the sudden reduction of the density floor which attains the final value of at .

3.1. The Ohmic-Resistivity-Only Run

Run O-b5 quickly develops into turbulence. From the upper panel of Figure 2, the separation between the highly turbulent active layer and the more or less quiescent midplane region is clearly seen. We further show the time and horizontally averaged vertical profiles of density, magnetic pressure, Maxwell and Reynolds stresses in Figure 3. The time averages are performed from onward.

The MRI turbulence generates strong magnetic field that buoyantly rises, and the disk surface quickly forms a strongly magnetically dominated corona beyond (Miller & Stone, 2000), as shown in the bottom panel of Figure 3. The gas density in the corona follows an exponential profile due to strong magnetic support, instead decreasing with height as an Gaussian in the hydrostatic case. The velocity in the coronal regions is highly supersonic, with turbulent kinetic energy exceeding the gas pressure beyond . A strong outflow is launched from the active layer of the disk as has been studied by Suzuki et al. (2010), who also found that the mass outflow rate scales roughly linearly with the vertical net magnetic flux. In our simulation, the time and horizontally averaged mass outflow rate is found to be about from each side of the box, comparable with the measurements in Suzuki et al. (2010) with the same . We note that the mass outflow rate is not a well-characterized quantity in shearing-box (Bai & Stone, 2013), these measurements should only be taken as a reference.

Figure 3.— Vertical profiles of the Maxwell stress and Reynolds stress (upper panel), as well as gas/magnetic pressure and kinetic energy (lower panel) in the fiducial run with only Ohmic resistivity (O-b5).

The disk midplane is too resistive to become MRI active, with the boundary of the active layer well described by the Elsasser number criterion: , where is the vertical component of the Alfvén velocity. The Maxwell stress remains very small at the midplane for the first 20 orbits. However, the midplane magnetic field is then gradually amplified, while the flow in the midplane remains more or less laminar. We note that the Maxwell stress in the midplane is even larger than most regions in the active layer, which contrasts with some previous simulations with a similar setup, where the Maxwell stress becomes very small near the midplane (Suzuki et al., 2010; Hirose & Turner, 2011; Okuzumi & Hirose, 2011). The main reason for the difference, while somewhat counterintuitive, lies in the usage of a much larger resistivity cap in our simulations555See last paragraph of Section 2. In fact, we start the simulation with a smaller resistivity cap of in code unit. At , the cap is raised to its final value , and the Maxwell stress at disk midplane rises shortly afterwards.. The large resistivity at the disk midplane in our simulations strongly suppresses electric current, leaving the horizontal magnetic field to be almost constant across the midplane (see the bottom panel of Figure 3). Therefore, the midplane field strength is largely set by the field strength in the active layer. Note that although this is phenomenologically similar to the “undead” zone proposed by Turner & Sano (2008), it is conceptually very different. We note that it is essential for the resistivity cap to be large enough so that the gas and magnetic field become decoupled at the cap value, and we have further tested that the simulation results are independent of the resistivity cap once its value is greater than or so.

By contrast, the Reynolds stress as well as kinetic energy clearly have a large dip in the undead region between . There is still random motion near the midplane due to the sound waves launched from the base of the active layer (Fleming & Stone, 2003; Oishi & Mac Low, 2009), though the velocity amplitude is at least one order of magnitude smaller than that in the active layer.

Integrating the profiles of the Maxwell and Reynolds stresses using (12), we obtain the Shakura-Sunyaev parameter , and . Assuming steady state accretion, one can express the accretion rate for a MMSN disk model as


Therefore, in the absence of AD, and at the fiducial location AU, we obtain an accretion rate of yr, large enough to account for the observed accretion rates in most T-Tauri stars.

Given the usage of more realistic resistivity profiles and chemistry models, as well as the much larger resistivity cap enabled by the super time-stepping technique, our run with pure Ohmic-resistivity deserves more discussion on its own right. Nevertheless, we only consider our run O-b5 as a test and reference case since our main point of interest is the effect of AD on the gas dynamics of PPDs, which makes a dramatic difference from the conventional picture of PPD accretion.

3.2. Zero Net Vertical Flux Run with both Ohmic Resistivity and AD

Figure 4.— Same as Figure 3, but for zero net vertical flux run with both Ohmic resistivity and AD (OA-nobz).Note the different scales.

Setting the initial density floor of , we find that the MRI turbulence sets in at the beginning, then the simulation gets stuck with relatively strong magnetic field (dominantly toroidal) accumulating near the vertical boundaries with relatively little activities. The flux does not escape and provides very little magnetic support due to the flat profile, and the density near the vertical boundaries stays at the value of . This situation implies that the gas near vertical boundaries tends to fall back towards disk midplane, but this is prevented by the outflow boundary condition. Therefore, we gradually lower the density floor over the course of the simulation and find that the phenomena of artificial magnetic flux accumulation and density cutoff at disappears when is brought down to , which is achieved at . The whole process is reflected in Figure 2. From this time onward, the system evolves smoothly and no artifacts are seen from the vertical boundaries. Also note that setting such small density floor of at the beginning would lead to dramatically small simulation timestep which is not quite realistic, hence gradual reduction of is necessary, although one has to be patient to get rid of the artifacts.

We see from Figure 2 that in the saturated state, the disk only has extremely weak level of MRI turbulence in a very thin layer between about H. Note that the FUV ionization front is located at , beyond which we find and the gas behaves almost as in the ideal MHD regime (since density profile is still largely Gaussian, one can still use Figure 1 as a reference for the profile of ), hence the MRI operates in this region. Below , the value of is of order unity or smaller, and there is no evidence of turbulent activity. This is consistent with unstratified simulation results of Bai & Stone (2011), where it was found that the MRI is suppressed for in the absence of net vertical magnetic flux.

There is some residual Maxwell stress around the disk midplane as a result of initial conditions that slowly diminishes over the course of the simulation. We extract the profiles of various physical quantities and show them in Figure 4, where the time average is performed between and . The MRI-active region exhibits as bumps in the stress plot, as well as the bumps in the kinetic energy plot. The magnetic field strength roughly stays constant through the entire disk, which is likely set by the turbulent activities in the active zone: the value of plasma (ratio of gas pressure to magnetic pressure) reaches order unity within the turbulent layer. Moreover, we find that the magnetic field is predominantly toroidal in the entire disk.

From the simulation, we find the Shakura-Sunyaev parameter using (12) to be , and . For such small level of stress, we find that the resulting accretion rate is only about yr, which is three orders of magnitude too small compared with observations.

We can compare this result with optimistic predictions of the MRI-driven accretion rate using the semi-analytical framework of Bai (2011a). Using this framework, we first extract the density profiles and the profiles of and . Assuming constant magnetic field strength across the MRI-active layer, the MRI-driven accretion rate can be expressed as (Equation (28) of Bai, 2011a)


where the integral is performed in regions where the MRI is permitted based on the criteria (20) of Bai (2011a). We scan the field strength to maximize , which gives yr as an upper limit of the MRI-driven accretion rate. We see that this value is a factor of more than larger than our simulation result. The main reason that our simulation yields an accretion rate that is significantly smaller than theoretical expectations is that the field geometry is not optimal: in both the ideal MHD (e.g., Hawley et al., 1995; Bai & Stone, 2013) and non-ideal MHD (e.g., Fleming et al., 2000; Bai & Stone, 2011), one finds stronger turbulence in the presence of a net vertical magnetic flux, and the optimistic estimate can possibly be achieved only when a optimal field geometry is realized.

This simulation has also demonstrated the physical significance of AD which drastically reduces the accretion rate compared with the Ohmic-only case. Although such a comparison may be unfair since our Ohmic-only simulation contains net vertical magnetic flux, Ohmic-only simulations with similar setup and zero net vertical flux generally yield total stress level that is only slightly smaller than ours (e.g., Turner et al., 2007; Ilgner & Nelson, 2008). Clearly, with AD taken into account, zero net vertical magnetic field geometry is far from being capable of driving rapid accretion in the inner region of PPDs.

3.3. Net Vertical Flux Run with Both Ohmic Resistivity and AD

Figure 5.— Vertical profiles of various quantities in the fiducial run with both Ohmic resistivity and AD and net vertical magnetic flux of (OA-b5). Upper left: gas pressure and magnetic pressure. Lower left: Elsasser numbers for Ohmic resistivity (), Hall term () and AD (), together with the plasma . Note that we evaluate while the Hall effect is not included in the simulations. Upper right: three components of gas velocity, where the bold green curve is for the vertical velocity. The Alfvén points are indicated as black dots, while open circles mark the base of the outflow. Lower right: three components of the magnetic field, where the bold green curve is for vertical field.

The inclusion of net vertical magnetic flux in our fiducial run with both Ohmic resistivity and AD (OA-b5) makes the field geometry more favorable for the MRI as we originally expect. The bottom panel of Figure 2 illustrates the time evolution of Maxwell stress. The initial field configuration is unstable to the MRI, which gives rise to channel-flow like behaviors and initiates some turbulent activities in the surface layer. However, we find that quite surprisingly, the system then quickly relaxes into a non-turbulent state in about 10 orbits with the MRI suppressed completely. The laminar configuration is then maintained for the remaining of the simulation time666To justify the validity of this result, particularly that it is not due to an unrealistic initial condition, we restart from the end of run O-b5 with AD turned on and find that also in about 10 orbits of time, the system settles to the same laminar state..

As the system settles to a completely laminar state, we extract the exact vertical profiles of various physical quantities, as shown in Figure 5. The magnetic field is strongest within about of the disk midplane, where the field is essentially constant due to the large resistivity and AD. There is no distinction between active layer and dead zone as the entire disk is laminar. The disk becomes magnetically dominated beyond . The FUV ionization front is located at about as seen in the sharp increase of and profiles, which also corresponds to the point where the gas density starts to deviate from Gaussian and follow an exponential profile.

Being a 3D time-dependent simulation, the fact that the system reaches a steady laminar state suggests the stability (particularly, against the MRI) of the configuration. This stability can be qualitatively understood using the criterion based on MRI simulations that include individual non-ideal MHD effects. In regions near the disk midplane where Ohmic resistivity dominates, the Ohmic Elsasser number is below one within about , too small for the MRI to operate (Turner et al., 2007). Beyond this region where AD is the dominant non-ideal MHD effect, which requires weak magnetic field for the MRI to operate, the magnetic field is too strong, as judged from Figure 16 of Bai & Stone (2011). Even the FUV ionization increases substantially beyond , the disk has already become magnetically dominated () in these regions. The suppression of the MRI can be understood as a result of magnetic field amplification in the surface layer during the initial growth from the MRI-unstable configuration. The MRI is quenched once the field becomes too strong for it to operate because of AD.

The fact that MRI is suppressed in the disk surface layer implies that the hypothesis of Bai (2011a) does not always hold, where it was assumed that the magnetic field can be amplified by the MRI to maximize the efficiency of the MRI. In other words, we see that the magnetic field is amplified to much greater extent that the MRI is suppressed. Nevertheless, this result is not inconsistent with the framework of Bai (2011a) since it only estimates the upper limit of MRI-driven accretion rate, and complete suppression of the MRI simply means the rate is zero.

The most prominent feature of the laminar state is a strong outflow that leaves the vertical boundaries. All three components of the velocities exceed the sound speed relative to the background Keplerian flow, with reaches about at the vertical boundaries, and an outflow mass loss rate of from each side of the simulation box. The Alfvén critical point, given by , is contained within our simulation box and is indicated in the upper right panel of Figure 5. Moreover, according to the conventional definition (Wardle & Koenigl, 1993), we define the base of the outflow at the location where the azimuthal velocity starts to become super-Keplerian, and the it will become more obvious later about this definition (see Figure 8) .

Our fiducial run OA-b5 suggests that the structure of the disk has a pure one-dimensional (1D) profile. To verify this result, we perform a new simulation (run S-OA-b5) with everything kept the same except that the horizontal domain is reduced to resolved by cells. We will refer to this type of runs as “quasi-1D” simulations. Such small box is obviously too small to study the MRI, but we find that except for different initial evolution of the original MRI-unstable configuration, the system relaxes to exactly the same laminar state with a strong outflow as in the fiducial run OA-b5. Moreover, we have further checked that although the initial evolution involves horizontal variations, such variations vanish as the system relaxes to the final steady state777However, if we start from a pure one-dimensional profile, the system will remain one-dimensional but does not relax to a steady state.. Therefore, we conclude that the inclusion of AD makes the disk structure purely 1D.

In sum, we have seen that the three fiducial simulations differ with each other in only one piece of physics, while show dramatically different characteristics. These results best demonstrate the importance of including AD in the study of gas dynamics in PPDs, as well as the significant role played by the magnetic field geometry. They strongly suggest that the MRI operates extremely ineffectively in the inner region of PPDs. In the mean time, the new laminar solution from our run OA-b5 strongly points to the magnetocentrifugal wind as the most promising driving mechanism for PPD accretion. The richness of the new findings deserves detailed study, and we devote the next section to address the nature of the new laminar wind solution. Moreover, the 1D nature of the new solution allows us to perform the simulations using very small horizontal domains which tremendously reduces the computational cost.

4. Nature of the Wind Solution

In this section, we do detailed analysis of our new laminar wind solution from run OA-b5. Table 2 summarizes the main physical properties of our fiducial laminar wind run together with a number of companion runs as they will be elaborated in the text.


4.69 6.10

4.63 6.22

4.56 6.23

4.45 5.33

4.70 7.12

4.83 7.67

: normalized Maxwell stress in the disk zone; : wind stress at wind base; : total mass loss rate (from both sides of the disk); : location of the wind base; : location of the Alfvén point; : rate of work done at the shearing-box boundaries; : energy loss rate due to Poynting flux; : kinetic part of the energy loss rate. All quantities are in natural unit ().

Table 2Wind Properties from All Fiducial Simulations.

4.1. Field Line Geometry and Wind Launching

We first consider the geometry of poloidal magnetic field lines, from which much insight can be gained on the wind launching mechanism. Using the magnetic field vectors from our run OA-b5, we integrate a poloidal field line from the disk midplane all the way to the vertical boundary of our simulation box, and show the results in Figure 6. In the mean time, we overplot the direction of velocity vectors as red arrows.

Figure 6.— The poloidal field line geometry in our fiducial run OA-b5 (blue solid line). Overplotted are the unit vectors of the poloidal gas velocity (red arrows). The location of the wind launching point, the plasma point, the FUV ionization front and the Alfvén point are indicated (black dash-dotted). Also marked is the location at the base of the wind (green dashed).

The field line is straight within about from the midplane due to the extremely large resistivity, where the gas and magnetic field are essentially decoupled as we discussed in the previous section. The field lines start to bend once the gas become partially coupled to the magnetic field, characterized by Ohmic as well as AD Elsasser numbers and exceeding unity, which occurs at , as can be seen from the bottom left panel of Figure 5. We label this point as the launching point in Figure 6. Beyond this point, the azimuthal magnetic field decreases rapidly with height, creating large current density in the radial direction. Together with the vertical field, we obtain the following force balance equation in the azimuthal direction


which states that the Lorentz force is balanced by the Coriolis force. This explains the increase of radial velocity in the upper right panel of Figure 5 beyond about , as well as the direction of the arrows in Figure 6 at the same locations. In this region, AD is the dominant non-ideal MHD effect. The radial motion of the gas makes the velocity vector deviate from the direction of the magnetic field, which drags and bends the magnetic field lines towards the same direction via the ion-neutral drag. Correspondingly, increases in strength.

The field lines tend to become straight again as gas density decreases and the flow becomes magnetically dominated with total plasma . Beyond this point, the gas only has limited effect on the field line geometry. Eventually, at the base of the wind (located at ), defined as the location where the azimuthal velocity exceeds the Keplerian velocity, the magnetocentrifugal acceleration starts to operate. Throughout this paper, we define the region between as the disk zone, while regions beyond as the wind zone.

We have also indicated the location of the FUV ionization front, beyond which the gas behaves more or less as ideal MHD due to the large ionization fraction. We see that the poloidal velocity of the gas is aligned with the poloidal magnetic field beyond the FUV front, as expected for an ideal MHD wind. Below the FUV front, gas velocity deviates from the direction of the magnetic fields as a result of AD and Ohmic resistivity. In this fiducial run with , the location of the FUV front overlaps with the base of the wind, which we note is just a coincidence. The importance of the FUV ionization will further be discussed in Section 5.2.

4.2. Conservation Laws

A laminar magnetized disk wind is characterized by the conservation of mass, angular momentum and energy along poloidal streamlines and magnetic field lines (which are aligned with each other in ideal MHD). They are very useful for diagnosing the mechanism for wind launching and acceleration (e.g., Pelletier & Pudritz, 1992; Spruit, 1996). In the shearing-box framework, the counterparts of these conservation laws are derived by Lesur et al. (2012). It was found that poloidal gas streamlines do not necessarily follow exactly the poloidal field lines, which introduces new terms in the conservation laws. Nevertheless, in the absence of strong magnetic field, the deviations are small (as we checked in our simulations) therefore, we just consider conventional terms in the treatment.

Figure 7.— Angular momentum (left) and energy (right) conservation along a streamline in the disk surface from to for the laminar wind solution in our fiducial run OA-b5. For the energy plot, the various terms in Equation (18) are represented by: (black dashed), (blue solid), (magenta solid), (green solid) and (red solid). For the angular momentum plot, the various terms in Equation (17) are represented by: (black dashed), (blue solid), (red solid).

Starting from mass conservation, we have


at each side of the disk. In practice, because we constantly add mass to the simulation box to maintain steady state, we find that does not become constant until beyond for our run OA-b5.

Specific angular momentum is conserved along streamlines, as long as is constant:


where is the fluid part of the specific angular momentum. The partition between and describes the angular momentum exchange between gas and magnetic field.

Energy conservation is expressed by the Bernoulli invariant in the ideal MHD regime, and we expect it to strictly hold only beyond the FUV ionization front (at ). It reads


where represents the specific energy along a streamline, with the four terms denoting kinetic energy, enthalpy, potential energy and the work done by the magnetic torque, respectively, and


Note that the full velocity (rather than ) enters , and the potential energy for the shearing-box is


To test these conservation laws, we integrate the velocity profiles shown in Figure 5 to obtain the streamline. We start from where approaches the constant value and trace the streamlines to the vertical boundary , setting at (note that by construction the conservation laws do not depend on the choice of the zero point of ). This range covers the most interesting regions of the disk where wind launching and acceleration take place.

In Figure 7 we show the vertical profiles of individual terms in the specific angular momentum and energy for streamlines obtained from run OA-b5. On the left panel, we see that the black dashed line flattens beyond , indicating angular momentum conservation. This is consistent with our expectations ( is conserved once approaches constant). It is clear that the increase of gas angular momentum () is compensated by the magnetic torque. On the right panel, we see that the black dashed line flattens beyond , indicating energy conservation. This is again consistent with our expectations ( is conserved once the gas obeys ideal MHD, i.e., beyond the FUV ionization front). The rise of the blue curve in the energy plot indicates flow acceleration, which is mostly compensated by the reduction of potential energy (red). This is the basic picture of magnetocentrifugal acceleration, which is not surprising since the bending angle of the field lines relative to disk normal in the wind zone is about , well above the threshold of . In this frame of reference, the magnetic torque plays a minor role, although its effect and the centrifugal potential are interchangable by shifting the zero point of (Bai & Stone, 2013). Thermal pressure gradient only plays a diminishing role in the acceleration process.

4.3. Angular Momentum Transport

Figure 8.— The vertical profiles of the components of the Reynolds (red dashed) and Maxwell (blue solid) stresses in our fiducial run OA-b5. The disk is divided into the wind zone and disk zone, given by .

Angular momentum transport can be achieved by radial redistribution, due to the component of the stress tensor as discussed in the beginning of Section 3, as well as by vertical extraction, due to the component of the stress tensor


where subscripts corresponds to the location for the top and bottom sides of the disk, which we choose to be at the base of the wind (the reason will become clear shortly). The torque from the outflow can be obtained by simply multiplying by the radius (which is unspecified in shearing-box).

The total rate of angular momentum transport, assuming steady-state accretion in the disk zone, can be obtained by generalizing equation (13) to include the wind torque, which gives


where we have assumed from symmetry considerations, which will further be discussed in the next subsection. We see that if and are comparable with each other, than vertical transport would be more efficient than radial transport by a factor of about . Since we are considering thin accretion disks (where shearing-box approximation applies), it is in general much easier for disk wind to transport angular momentum compared with turbulent transport.

We start by examining the radial transport of angular momentum. Although radial transport is usually a result of turbulence which gives correlated fluctuations of velocity and magnetic field in the radial and azimuthal direction, ordered velocity/field structure also produce radial transport. In Figure 8, we show the vertical profiles of the components of the Reynolds and Maxwell stresses. Clearly, even there is no turbulence, a non-zero Maxwell stress given by ordered magnetic fields dominates the radial transport in the disk zone. This corresponds to the undead zone scenario discussed in Turner & Sano (2008). We also see that contribution from the Reynolds stress in the disk zone is completely negligible. Integrating the Maxwell stress across the disk zone, we find .

Meanwhile, Figure 8 best demonstrates the advantage for dividing the profile of our the laminar wind solution into the disk zone and the wind zone: the Reynolds stress is by definition zero at the base of the wind, and increases rapidly in the wind zone as a result of radially-outward and super-Keplerian motion consequent of the magnetocentrifugal acceleration; the Maxwell stress has a prominent peak at the base of the wind, and decreases in the wind zone, as magnetic energy is converted to the kinetic energy of the outflowing gas.

For the vertical transport, we measure the component of the stress tensor at the base of the wind. Again, our choice of the wind base shows its advantages: by definition, Reynolds stress is zero, and one only needs to consider the Maxwell stress, which we find is . Hereafter, we will refer to the component of the Maxwell stress at the base of the wind as “wind stress”. We note that since where is constant, from Figure 5 we see that varies smoothly within the disk and should peak between . Although wind stress eventually drives accretion, a small fraction of the gas flow undergoes “decretion” to bend magnetic field lines (between the launching point and the wind base, see Figure 6). This outwardly directed gas flow is properly accounted for by measuring the wind stress at rather than at .

Applying the MMSN disk model into Equation (22), we can estimate the accretion rate driven by the combined effect of radial and vertical transports


where is the accretion rate measured in yr. We see that radial transport by the Maxwell stress (from ordered magnetic field) in the disk zone can drive accretion for about yr. This number is already significantly larger than the case for run OA-nobz with zero net vertical magnetic flux, yet still too small to account for the typical accretion rate for PPDs. The accretion driven by wind transport, on the other hand, reaches yr, which is sufficient to account for the observed accretion rate for most PPDs.

We see that our laminar wind solution simultaneously solves two problems facing the conventional MRI scenario and the conventional wind scenario described in Section 1. First, it efficiently drives disk accretion that matches the rates inferred from observations, which can hardly be achieved in the MRI scenario. Second, it launches the magnetocentrifugal wind with only a very weak net vertical magnetic field, in contrast with the conventional understanding that requires equipartition field at disk midplane. Therefore, the rate of wind transport in our scenario is only moderate and does not result in rapid depletion of the disk.

4.4. Symmetry and Strong Current Layer

Figure 9.— Cartoon illustration of the geometry of wind magnetic field lines in shearing-box simulations. The shearing-box approximation ignores all radial gradients (except the shear), hence does not contain information about the location of the central object. The natural wind geometry launched from shearing-box has the “odd-” symmetry where the field lines on the top and bottom of the box bend toward opposite directions (left). In a physical situation, where the location of the central object is fixed, the field lines on the top and bottom sides of the box should bend to the same direction away from the central object.

The previous subsection demonstrates the success of the magnetocentrifugal wind scenario as the the driving force of disk accretion. However, one problem was ignored in the previous discussion. The laminar wind solution we obtained obeys the “odd-” symmetry:


As a result, the wind field lines on the top and bottom sides of our simulation box bend to opposite radial and toroidal directions, as illustrated on the left of the Cartoon picture in Figure 9. However, physically, one expects the field lines on the top and bottom sides of the box to bend toward the same direction, pointing away from the central star, as illustrated on the right of Figure 9. A particularly disturbing fact is that in the case with odd- symmetry, the “wind” does not transport angular momentum, since at the top and bottom sides of the disk have the same sign and value, and cancel out. The wind angular momentum transport mechanism works only for a physical wind geometry.

We note that all radial gradients except the azimuthal shear are neglected in the shearing-box approximation, hence one can not tell whether the location of the central star is in the “inner” or “outer” (radial) side of the box. In other words, the two radial sides of the box are symmetric. In our simulations, we find that the direction that the wind field lines bend as we start from the initial configuration is totally random, and the top and bottom sides of the disk evolve independently888We have performed more than one copies of each 3D and the quasi-1D fiducial simulations that verifies the randomness.. This suggests that the chances to find a solution with the unphysical odd- symmetry and with a physical geometry are equal. In fact, our 1D run S-OA-b5 results in a solution with the physical geometry.

In Figure 10, we show the profiles of the magnetic field and velocity field from our run S-OA-b5-F3 that have the physical wind geometry. We see that this physical solution is NOT a solution under the “even-” symmetry:


which is almost exclusively used for constructing semi-analytical local wind solutions (Wardle & Koenigl, 1993; Ogilvie & Livio, 2001; Königl et al., 2010; Salmeron et al., 2011). Very interestingly, the physical wind solution follows exactly the odd- symmetry solution we obtained before, as can also be seen from Table 2 that all physical quantities measured in this run agree with those in our 3D fiducial run OA-b5 within . In particular, the wind stress is almost exactly the same. The only difference is that the horizontal field lines (and horizontal velocities) are flipped at one side of the box to achieve the physical wind geometry. This flipping is mediated by a sharp transition of the horizontal fields, which exhibits as a strong current layer. The strong current layer is not located at the midplane, but is located at about (or at , and the selection is random). Therefore, the physical solution does not have any symmetry about the disk midplane.

Figure 10.— The vertical profiles of the magnetic fields (left) and velocities (right) in the laminar solution that has the physical wind geometry (run S-OA-b5, see Section 4.4). The inset on the right panel shows the zoomed-in view of the velocity profiles at the strong current layer. This is essentially how accretion proceeds in our wind-driven scenario.

The location of the strong current layer roughly corresponds to where both the Ohmic and AD Elsasser numbers and become greater than , as seen from the bottom left panel of Figure 5. The reason for such large offset from the disk midplane is that magnetic diffusion near the disk midplane is so strong that the gas and magnetic field are essentially decoupled. Correspondingly, the magnetic field lines must be straight and current is excluded, and the strong current layer can exist only when magnetic diffusion becomes weaker.

The right panel of Figure 10 demonstrates the velocity profile of the physical wind solution. The strong current layer exhibits as a seemingly small feature at about in this plot, which is in linear scale. As we zoom into this region shown in the inset, we find that the gas has a large inflow velocity of about at the location of the strong current layer. The large inflow is directly the consequence of the wind stress: the Maxwell stress exerted at the base of the wind is released at the strong current layer, driving a large inflow which is essentially how accretion proceeds. More specifically, it is the balance between the Lorentz force and Coriolis force that leads to the large inflow and accretion in the strong current layer (see equation (15)).

To further demonstrate the effectiveness of accretion in the current layer, we note that for a MMSN disk, in order to have accretion to approach yr, a bulk radial drift velocity of about is required


In the physical wind solution, we find that only a tiny fraction (about ) of disk mass is contained in the strong current layer, but it drifts at very large velocities (). The combined effect is an efficient accretion that easily accounts for the typical accretion rates observed in PPDs. More specifically, by integrating the radial mass flux () across the disk zone, we find the mean radial velocity of the gas to be . Plugging to Equation (26), this would lead to an accretion rate of yr, which agrees very well with the estimate in Section 4.3.

Finally, we note that the strong current layer is not a thin current sheet, since at its location there is still substantial magnetic diffusion (mainly AD). It has a thickness of about , which is well resolved in about cells in our quasi-1D simulations. We have found that the strong current layer is stable in our quasi-1D fiducial run S-OA-b5. The stability of the strong current layer in 3D requires justification in the more geometrically appropriate global simulations999We have repeated our 3D fiducial run OA-b5 with a different random seed which evolves to a physical wind solution with a strong current layer. We find that the strong current layer is maintained for about 100 orbits but eventually escapes the simulation box and the odd- symmetry is recovered. This is likely to be due to the intrinsic limitations of the shearing-box, where curvature terms are neglected and hence favors the odd- symmetry solution. Moreover, as the fast magnetosonic point is beyond our simulation box, the stability of the strong current layer can also be affected from outside of the simulation domain..

4.4.1 Wind Solution from Even- Symmetry

Figure 11.— Same as Figure 10, but for run E-OA-b5 where an even- symmetry (25) is enforced (see Section 4.4.1). The inset on the right panel shows the zoomed-in view of the radial velocity profile in the disk zone.

We have discussed that the odd- symmetry is the most natural symmetry in shearing-box but it is unphysical for a disk wind. A physical wind solution can be achieved by flipping the horizontal field and stream lines at one side of the disk, which yields a solution that contains a strong current layer offset from the disk midplane. Because the two solutions have exactly the same properties except the flip, it is natural to classify them as the same solution. In this subsection, we further consider a solution that obeys even- symmetry.

We design a new 3D simulation E-OA-b5 which is exactly the same as our fiducial run OA-b5 except that we use only half of the shearing-box, with the lower boundary located at the disk midplane. At the lower boundary we use the standard reflection/conduction boundary condition so that even- symmetry is enforced. We find that after a brief period of the MRI, the system also settles into a pure laminar state, with velocity and magnetic field profiles shown in Figure 11 (where we have recovered the other side of the disk based on even- symmetry). Under even- symmetry, the horizontal magnetic field at the midplane is enforced to be zero. Due to the large resistivity (which can not support a significant current), the field is kept very close to zero within where both Ohmic and AD Elsasser numbers are much less than . Beyond , the horizontal field increases, and then beyond , the field configuration recovers to almost exactly the same solution we obtained under odd- symmetry.

We note that the field configuration similar to our even- symmetry solution has previously been reported by Li (1996) and Wardle (1997), and it was discussed that in the presence of low conductivity in the midplane, wind solution can be constructed with midplane magnetic field pressure much smaller than the gas thermal pressure. Our solution confirms their findings and demonstrates it robustness and uniqueness since it is constructed in a self-consistent and evolutionary manner.

The velocity profiles from our even- symmetry solution are exactly the same as the odd- symmetry solution in the wind zone (except the flip). In the disk zone, we see that the radial velocity has two dips located at . The dips are much shallower (only about ) than the dip due to the strong current layer in our run S-OA-b5, but they are much wider. Integrating the over height, we obtain the radial mass flux to be , which is almost exactly the same as in run S-OA-b5. This is not surprising since the properties of the wind in the two solutions are exactly the same.

In sum, we find that the suppression of the MRI and wind launching is inevitable in the inner region of PPDs around 1AU, and the property of the solution in the wind zone is independent of the large-scale field geometry (or symmetry).

4.5. Outflow

Strong outflows are launched in our fiducial simulations, and we have measured that the rate of mass outflow from each side of the box is about from the quasi-1D run and from the 3D run. In Table 2, we also list total mass loss rate,


as the sum of the mass loss rate from the two vertical boundaries for our fiducial runs. Although mass loss is not significant for the duration of our simulations, the value is relatively high that without replenishing to the disk, the mass loss timescale is only , which amounts to only about years. Moreover, the measured wind mass loss rate is in fact comparable to the mass accretion rate driven by the wind itself, which is inconsistent with the observationally inferred ratio that the mass loss rate is only about of the accretion rate (Cabrit et al., 1990; Hartigan et al., 1995).

We note that, however, the magnetocentrifugal wind is intrinsically a global phenomenon. A global wind solution is not fully determined until all critical points (or surfaces), namely, the slow magnetosonic, Alfvén and fast magnetosonic points are passed, where the fast magnetosonic point is far beyond the extent of our simulation box. This leads to one extra degree of freedom in the wind solution. Moreover, the global structure of the wind, and the location of the critical points, are set by the interplay between the microphysics in the disk, and the large-scale field structure (i.e., the magnetic flux distribution over the entire disk), where the latter is not reflected in local shearing-box simulations. Therefore, while useful for studying wind-launching, the shearing-box approximation has its limitations and the measured rate of outflow from our simulations should only be taken as reference and it may significantly overestimate the mass outflow rate in reality.

To further clarify the discussion above, we perform three additional quasi-1D simulations, changing the box size to , and respectively, and the three runs are labeled as S-OA-b5-12H, S-OA-b5-20H and S-OA-b5-24H. We measure a series of quantities discussed before and also list the results in Table 2. We find that as we increase the height of the simulation box, the location of the Alfvén point moves to higher altitudes. In the mean time, the mass loss rate decreases. In fact, these two quantities are connected, and in shearing-box, the definition of the Alfvén point indicates (Bai & Stone, 2013)


where denotes the value at the Alfvén point. Therefore, higher location of the Alfvén point indicates lower density, hence smaller mass loss rate.

The trend of deceasing mass outflow rate with vertical box size reflects one limitation to shearing-box: the gravitational potential increases quadratically with vertical height, which increasingly hinders the gas from escaping. The outflow in our simulations is possible because gravitational potential is cut off at the vertical boundaries. This limitation of shearing-box has been noticed and discussed in a number of recent works on ideal MHD simulations of the MRI (Fromang et al., 2012; Lesur et al., 2012). Because the depth of the shearing-box potential is about the same as the real depth of the potential at , Bai & Stone (2013) suggest that the mass loss rate can be roughly estimated by performing shearing-box simulations with a box height of , or to study the height dependence of and extrapolate the result to . In MMSN, we have at 1 AU (note that is not specified in the shearing-sheet approximation). Such a thin disk would require a shearing-box height of , which is numerically unpractical. Based on our three quasi-1D runs with to , the trend can be roughly described by . Extrapolating this trend, we obtain a rough estimate of mass loss rate to be in code unit. This level of brings the mass loss time scale to about years, which is about an order of magnitude longer than our original naive estimate. Although still somewhat short, the mass reservoir in the outer disk (which contains much more mass than that in the inner disk) may be sufficient to feed the inner disk to compensate for the mass loss.

It is also interesting to note that although the mass outflow rate is largely affected by the vertical box size, the angular momentum transport depends on the vertical box size very weakly. For box sizes of , , and , the values of at the base of wind are found to be , , and respectively, hence increasing box size leads to more or less the same (but slightly stronger) wind-driven accretion.

In sum, the trend that the mass outflow rate decreases with box size and mass accretion rate increases with box size implies that in a real system, the ratio of mass accretion rate to mass loss rate would be a factor of several larger than that obtained in our simulations, and would be more consistent with observations.

4.6. Energetics

Although our simulations assume an isothermal equation of state where total energy is not conserved, it is important to verify that our simulation results are energetically feasible and to examine the energy balance for real situations.

In shearing-box simulations, the energy is injected due to the work done by the shearing-box boundaries, the rate of which is given by (Bai & Stone, 2013)


where the integral over is performed within (the disk zone). Energy is lost through the open vertical boundaries, given by


where is the unit vector pointing away from the disk in the vertical direction, and we sum over contributions from the top and bottom of the wind base. The first terms in the bracket represent the kinetic energy loss and the work done by the mass outflow, and the second term represents energy loss from the Poynting flux. We have ignored the energy loss term associated with internal energy loss due to the mass outflow, since we keep feeding mass to the system which balances this term exactly.

In Table 2 we also list the values of , and , normalized in natural units. We see that the sum of and comprises of about of the total work done by the shear , hence energy conservation is not violated. The extra energy is likely to be radiated away in real systems.

4.7. Radial Drift of Magnetic Flux

Figure 12.— The toroidal EMF profile in our fiducial run OA-b5. Blue dashed: inductive EMF (); red dash-dotted: AD EMF (); black solid: the total EMF (). The Ohmic EMF is very close to and is not plotted. The EMFs are normalized to and are in natural unit of the sound speed .

Globally, poloidal magnetic flux drifts radially at velocity in the presence of a net toroidal electromotive force (EMF, )


Negative would lead to accumulation of magnetic flux in the inner disk while positive would make the disk lose magnetic flux to large scales. Moreover, should be constant with height so that magnetic flux drifts uniformly across the disk. In global models, the radial drift velocity is determined by the radial distribution of magnetic flux (Teitler, 2011). In local models, it is often chosen as a free parameter in local shearing-box models of disk winds (Wardle & Koenigl, 1993) due to the extra degree of freedom discussed in Section 4.5. Our simulations in principle have one extra degree of freedom because the fast magnetosonic point is not contained in the simulation box. We expect that the value of in our simulations has direct correspondence with this extra degree of freedom.

The toroidal component of the EMF in our simulations is the sum of the inductive EMF (), Ohmic EMF () and AD EMF (), which reads


In Figure 12 we show the radial profiles of the three toroidal EMFs and their sum, all normalized to . We see that in most regions of the disk, the inductive EMF is essentially exactly balanced by the AD EMF. We did not plot the Ohmic EMF because it is always negligibly small due to the very small current near the midplane and the very small resistivity in the upper layer. The sum of all EMFs, is consistent with zero at all locations of the disk. The little sharp peaks of the total EMF at around are purely numerical features due to the sharp variations of the AD coefficient and current density at the FUV ionization front.

Therefore, the wind solution we obtained from shearing-box simulations is a zero radial drift solution , which corresponds to a stationary magnetic configuration in the global picture. A direct consequence of is that poloidal field lines and poloidal streamlines follow each other in the absence of non-ideal MHD terms, as clearly shown in Figure 6. The fact that in our simulations is likely to be a result of our vertical boundary condition, where there is no vertical gradient of magnetic fields. We note that in the ideal MHD shearing-box simulations of wind launching process by Lesur et al. (2012), significant deviation between poloidal field lines and streamlines are found (see their Figure 3). This is likely to be a result of a different vertical boundary condition: poloidal field is forced to be vertical at ghost zones, resulting in a current sheet at the vertical boundary.

5. Parameter Study

To assess the importance of various physical effects on the properties of the laminar wind solution, we conduct a parameter study by surveying four parameters and compare the results from our fiducial model. All simulations in this parameter study are quasi-1D. First, we vary the vertical net magnetic flux, and consider , and . These runs are labeled as S-OA-b, where b denotes . Second, we vary the penetration depth of the FUV ionization, and consider g cm and g cm, and the run labels are attached with “F1” and “F10” respectively. Third, we consider the effect of grain abundance, and consider the case with grain-free chemistry (whose run label is attached with “nogr”) and the case with grain abundance (run label attached with “gr01”, indicating depletion factor of 0.1), 10 times the value in our fiducial run. Finally, we vary the surface density of our disk model, and consider and . In the mean time, we maintain the same field strength in physical unit, thus in these two runs are set to and respectively. These two runs are attached with “M3” and “M03”, respectively. A list of all these runs are given in Table 3.

S-OA-b3 1 0.03
S-OA-b4 1 0.03
S-OA-b6 1 0.03
S-OA-F1 1 0.01
S-OA-F10 1 0.1
S-OA-nogr 1 0.03 0
S-OA-gr01 1 0.03
S-OA-M3 3 0.1
S-OA-M03 0.3 0.1
Table 3Summary of All Simulations in the Parameter Study.

3.95 8.00

3.92 7.47

4.63 6.22

4.52 5.52

4.17 7.11

4.20 5.52

4.55 6.30

4.92 5.86

4.89 6.23

4.38 6.09

Note: we have duplicated run S-OA-b5 as a standard reference.

Table 4Wind Properties from All Simulations in the Parameter Study.

All the additional quasi-1D simulations are run to