Global Simulations of the Inner Regions of Protoplanetary Disks with Comprehensive Disk Microphysics

The gas dynamics of weakly ionized protoplanetary disks (PPDs) is largely governed by the coupling between gas and magnetic fields, described by three non-ideal magnetohydrodynamical (MHD) effects (Ohmic, Hall, ambipolar). Previous local simulations incorporating these processes have revealed that the inner regions of PPDs are largely laminar accompanied by wind-driven accretion. We conduct 2D axisymmetric, fully global MHD simulations of these regions ($\sim1-20$ AU), taking into account all non-ideal MHD effects, with tabulated diffusion coefficients and approximate treatment of external ionization and heating. With net vertical field aligned with disk rotation, the Hall-shear instability strongly amplifies horizontal magnetic field, making the overall dynamics dependent on initial field configuration. Following disk formation, the disk likely relaxes into an inner zone characterized by asymmetric field configuration across the midplane that smoothly transitions to a more symmetric outer zone. Angular momentum transport is driven by both MHD winds and laminar Maxwell stress, with both accretion and decretion flows present at different heights, and modestly asymmetric winds from the two disk sides. With anti-aligned field polarity, weakly magnetized disks settle into an asymmetric field configuration with supersonic accretion flow concentrated at one side of disk surface, and highly asymmetric winds between the two disk sides. In all cases, the wind is magneto-thermal in nature characterized by mass loss rate exceeding the accretion rate. More strongly magnetized disks give more symmetric field configuration and flow structures. Deeper far-UV penetration leads to stronger and less stable outflows. Implications for observations and planet formation are also discussed.


INTRODUCTION
Planet formation takes place in protoplanetary disks (PPDs) surrounding young stars. Composed of gas and dust, PPDs offer rich observational diagnostics that help constrain the physical scenarios of planet formation. With typical lifetime of a few Myrs (Haisch et al. 2001), PPDs are known to be rapidly accreting, with typical accretion rate of ∼ 10 −8 M ⊙ yr −1 (e.g., see Hartmann et al. 2016 for an updated review). The accretion phenomenon is closely related to jets and outflows that are ubiquitous among young stellar objects (e.g., see Frank et al. 2014 for a recent review). The solid materials are primarily probed by the dust thermal emission at sub-millimeter (mm) as well as in scattered light at infrared wavelengths (e.g., see Williams & Cieza 2011).
With the advent of the Atacama Large Millimeter/sub-millimeter Array (ALMA) and extreme adaptive optics systems such as SPHERE/VLT and the Gemini Planet Imager, PPDs have revealed rich substructures (e.g., ALMA Partnership 2015;Nomura et al. 2016;Pérez et al. 2016;Isella et al. 2016;Ginski et al. 2016;de Boer et al. 2016), even down to AU scales . Furthermore, astrochemistry has emerged to provide more refined information about the physical environments of PPDs with implications for planetary composition (e.g., Öberg et al. 2011;Henning & Semenov 2013;Qi et al. 2013;Cleeves et al. xbai@cfa.harvard.edu

2016).
Theoretically, the gas dynamics of PPDs, especially the local and global disk structure as well as internal flow structure (such as the level of turbulence), plays a crucial role in almost all aspects of planet formation (Armitage 2011). This is because small dust grains are coupled with the gas aerodynamically, whereas larger bodies are coupled with the disk gravitationally. In particular, dust grains always migrate towards higher pressure, leading to radial drift and particle trapping at pressure maxima, and these theoretical predictions have found observational support (e.g., Pinilla et al. 2012;Birnstiel & Andrews 2014;Zhang et al. 2016). For planets that form within the disk lifetime, planet-disk interaction leads to planet migration, and its direction and rate sensitively depend on the radial gradients of various disk quantities, as well as the level of turbulence (e.g., see Baruteau et al. 2016 for a recent review).

Current Understandings of PPD Gas Dynamics
The central question on PPD gas dynamics lies in the mechanism of angular momentum transport, which shapes the disk structure and drives disk accretion and global evolution. Angular momentum can be transported radially (viscous accretion), mainly mediated by turbulence or large-scale magnetic stress, or vertically, mediated by a magnetized disk wind (see Turner et al. 2014 for a recent review). In either scenario, magnetic field is believed to play an essential role, as we briefly discuss below, and this is further supported from paleomagnetic studies of the Semarkona chondrite (Fu et al. 2014).
In the absence of magnetic fields, a number of mechanisms have been studied, such as the verticalshear instability (Nelson et al. 2013;Stoll & Kley 2014;Lin & Youdin 2015), the convective overstability (Lesur & Papaloizou 2010;Lyra & Klahr 2011;Klahr & Hubbard 2014;Lyra 2014), and the zombie vortex instability (Marcus et al. 2013(Marcus et al. , 2014. Nevertheless, besides the fact that these instabilities all require certain thermodynamic conditions to operate, the resulting level of turbulence is typically weak, with Shakura-Sunyaev α reaching at most ∼ 10 −3 , which is too small to account for the accretion rates in the bulk disk population (e.g. Andrews et al. 2009Andrews et al. , 2010. Other mechanisms, such as the gravitational instability (Gammie 2001;Rafikov 2009), and spiral-density waves driven by envelop infall , can provide significant angular momentum transport, although only at the very early stages of PPD evolution.
With magnetic fields, the key microphysical processes involve determining how well magnetic fields are coupled with the gas. The ionization of PPDs largely relies on external sources such as cosmic-rays and X-rays, leading to extremely low level of disk ionization with vertically stratified ionization structure (Sano et al. 2000;Ilgner & Nelson 2006;Bai & Goodman 2009). As a result, magnetic fields are no longer frozen in to the bulk gas as in ideal magnetohydrodynamics (MHD), introducing three non-ideal MHD effects: Ohmic resistivity, the Hall effect and ambipolar diffusion (AD). The three effects control the gas dynamics in different ways, and at fixed field strength, the dominant effect transitions from resistivity to the Hall effect, and to AD as density decreases (Wardle 2007;Bai 2011).
Conventionally, the magnetorotational instability (MRI, Balbus & Hawley 1991) has been considered as the dominant mechanism to drive disk accretion. However, the MRI is strongly affected by non-ideal MHD effects. Linear modes are damped by resistivity and AD (Blaes & Balbus 1994;Jin 1996;Kunz & Balbus 2004), whereas the Hall effect modifies the dispersion relation depending on the polarity of vertical field threading the disk (Wardle 1999;Balbus & Terquem 2001). Taking only Ohmic resistivity into account, the picture of layered accretion has been established (Gammie 1996;Fleming & Stone 2003;Turner & Sano 2008;Oishi & Mac Low 2009), where the MRI is suppressed by resistivity in the midplane region of the inner disk (i.e., the densest region where resistivity dominates), while it operates in the much better ionized, low-density disk surface, driving viscous accretion through the surface layer.
The conventional picture of layered accretion no longer holds when AD is further taken into account. As AD dominates towards low density regions, the MRI is found to be almost entirely suppressed in the inner region of PPDs ( 15 AU, Bai & Stone 2013;Bai 2013;Gressel et al. 2015). The MRI is substantially damped in the low-density outer disk (Simon et al. 2013b,a;Bai 2015), whereas it can, however, operate in the disk surface thanks to far-UV (FUV) ionization (Perez-Becker & Chiang 2011). In other words, layered accretion is likely more applicable to the outer instead of the inner disk. As the MRI is suppressed or damped, efficient angular momentum transport requires the disk to be threaded with net vertical magnetic field, and a magnetized disk wind is likely the primary mechanism to drive disk accretion.
The inclusion of the Hall effect further makes the gas dynamic depend on the polarity of the net vertical field. In particular, when vertical field is aligned with the disk rotation axis, horizontal components of the field are amplified due to the Hall-shear instability (Kunz 2008), which leads to enhanced radial transport of angular momentum by large-scale magnetic stresses Bai 2014Bai , 2015Simon et al. 2015). In the opposite case of anti-aligned vertical field, horizontal field is reduced towards zero. In both cases, disk winds likely remain the dominant mechanism to drive disk accretion (Bai 2014). In the outer disk where the MRI is damped, the Hall effect can enhance/reduce turbulence depending on polarity (Sano & Stone 2002a,b;Bai 2015;Simon et al. 2015), though not substantially.
MHD disk winds have been studied extensively in the literature, ranging from global self-similar solutions (e.g., Blandford & Payne 1982;Li 1995;Ferreira 1997) to local solutions that match to such solutions (e.g., Wardle & Koenigl 1993;Königl et al. 2010;Salmeron et al. 2011), to global simulations (e.g., Krasnopolsky et al. 1999;Pudritz et al. 2006;Zanni et al. 2007;Tzeferacos et al. 2009). The wind properties are now well known to depend mainly on the strength and distribution of the magnetic flux threading the disk, and on the mass loading, with the latter mainly being controlled by disk physics and thermodynamics. We emphasize that previous studies are not directly applicable to PPDs because the main disk microphysics (all three non-ideal MHD effects with realistic ionization structure) and thermodynamics were not properly taken into account. Most studies considered vertical field strength near equipartition (to avoid the development of the MRI), and the resulting wind-driven accretion rate would be orders of magnitude larger than typical PPD accretion rates. On the other hand, realistic local simulations have demonstrated that a weak vertical field can naturally sustain wind launching, and drives disk accretion at desired accretion rates (Bai & Stone 2013;Bai 2013).

Outstanding Issues
These recent works have revealed rich disk physics resulting from the non-ideal MHD effects, highlighting the importance of incorporating realistic disk microphysics for studying PPD gas dynamics. However, the aforementioned works are mostly local simulations. One exception is Gressel et al. (2015), where the simulations were radially global, yet vertically local, and the Hall effect was not included. Another exception is the very recent work of Béthune et al. (2017), which we will discuss in more detail in Section 8.2. Three outstanding issues remain to be worked out and clarified.
First, wind kinematics. Local simulations fail to cover the full depth of the gravitational potential well of the central star, and the wind mass loss rate has been found to depend on the vertical height of the simulation box (Fromang et al. 2013;Bai & Stone 2013). This effect can be understood from a different point of view: the wind is artificially truncated by the imposed boundary conditions. To overcome this limitation, the computational domain must be sufficiently extended so that wind flow passes major critical points (i.e., wind velocity exceeds sonic/Alfvén speed) and loses causal connection with the disk surface.
In anticipation of this work, we developed a semianalytical theory for MHD disk winds from PPDs in Bai et al. (2016) that allows the flows to pass all critical points to arrive at unique solutions. We have also adopted an approximate treatment of the thermodynamics of the disk wind to mimic FUV/X-ray heating, which has conventionally been considered to drive photoevaporation (as a pure thermal wind, see Alexander et al. 2014 for a review). The wind solutions, which we call magneto-thermal disk winds, indicate that PPDs lose mass from disk winds at a rate comparable to winddriven accretion rate, and wind kinematics is most sensitive to poloidal field strength as well as how deep FUV/X-ray can penetrate into (and hence heat and ionize) the disk. These predictions remain to be verified and calibrated through realistic global simulations.
Second, the symmetry issue. A physical wind geometry requires that poloidal field lines bend away from the central star, which further requires that toroidal field must change sign across the disk. The toroidal field gradient is also directly associated with the torque exerted by the disk wind, which drives the accretion flow. In Bai & Stone (2013) and Bai (2013), we found in local shearing-box simulations that in the inner disk, the flip of toroidal field occurs at a location that is offset from the midplane, resulting in symmetry breaking, which is later confirmed in global simulations (Gressel et al. 2015). Including the Hall effect with aligned vertical field, the Hall-shear instability amplifies the toroidal field near the midplane so strongly that in local simulations, toroidal field of one sign overwhelms and no flip can be sustained in the simulation box (Bai 2014). While we speculated that this is an artifact of vertical boundary condition, global simulations are essential to resolve this issue. This is also crucial to determine the global field configuration and flow structures in PPDs.
Third, the origin and evolution of magnetic flux. As discussed earlier, the presence of net poloidal magnetic flux is essential to drive disk accretion. We further showed more quantitatively that global disk evolution is primarily governed by the amount of flux threading the disks, and its radial distribution ). Therefore, transport of magnetic flux in PPDs is an even more fundamental question. In Bai & Stone (2017), we conducted preliminary studies of magnetic flux transport in PPDs, and emphasized the unique roles played by the Hall effect and AD that has been overlooked in the literature (Lubow et al. 1994;Okuzumi et al. 2014;Guilet & Ogilvie 2014). We found that when the disk is laminar, magnetic flux is systematically transported outward in a polarity-dependent manner, with rate of transport being faster in the anti-aligned case. We also noted that the exact rate of transport can be sensitive to the details of the disk ionization structure, which we treated very roughly, and more realistic simulations are needed to quantitatively characterize the rate of magnetic flux transport in PPDs.

This Work
We aim to conduct global simulations of PPDs that incorporate the most realistic disk microphysics. A unique aspect of our simulations is that we properly resolve the most important disk microphysics, and in the mean time the computational domain extends all the way to near the polar region, which is essential to accommodate the launching and propagation of MHD disk winds, as well as to accommodate magnetic flux evolution. Our simulations include all three non-ideal MHD effects with realistic treatment of disk ionization chemistry, together with approximate treatment of disk thermodynamics. These simulations are made possible thanks to the newly developed Athena++ MHD code (Stone et al. in preparation). We have implemented all non-ideal MHD effects (Bai & Stone 2017), and carefully designed simulation setup that circumvent difficulties especially associated with boundary conditions. With these simulations, we are able to address the aforementioned three major issues simultaneously, offering the most realistic PPD simulations to date. This paper is organized as follows. In Section 2, we provide detailed descriptions of the numerical method and simulation setup. Main diagnostics are discussed in Section 3. We focus on three fiducial simulations and analyze the results of each simulation in detail in Sections 4-6. A brief parameter study is conducted in Section 7. The results are discussed in broader contexts in Section 8, and we summarize and conclude in Section 9.

METHOD
We use the newly developed grid-based higher-order Godunov MHD code Athena++ (Stone et al., in preparation), which is the successor of the widely used Athena MHD code (Gardiner & Stone 2005, 2008Stone et al. 2008), to carry out global simulations of PPDs in this work. Athena++ works for curvilinear coordinate systems (here we use spherical-polar coordinates), where geometric source terms are properly implemented that guarantees angular momentum conservation. It also employs flexible grid spacings, allowing simulations to be performed over large dynamical ranges. All three nonideal MHD terms have been implemented (Bai & Stone 2017), which are the key microphysical processes for our simulations.

Dynamical Equations
Using Athena++, we solve the MHD equations in conservative form, including non-ideal MHD effects where ρ, v, and P are gas density, velocity and pressure, P * = P + B 2 /8π is total pressure, P * ≡ P * I with I being the identity tensor, B is magnetic field, with b ≡ B/B being the unit vector along the field direction. Total energy density is given by E = P/(γ − 1) + ρv 2 /2 + B 2 /8π, where γ is the adiabatic index, and J is the current density, with J ⊥ = −(J × b) × b being the component of J that is perpendicular to the magnetic field. We specify the static gravitational potential of the protostar as Φ = −GM * /r, with M * = M ⊙ . Magnetic diffusivities are represented by η O , η H and η A for Ohmic resistivity, the Hall effect and ambiploar diffusion (AD), which depend on ionization chemistry to be described in Section 2.4. Thermodynamics is mainly controlled by a thermal relaxation term Λ c in Equation (4). Note that in this equation, we have not included the Poynting flux (and hence heating) from non-ideal MHD terms. While not fully self-consistent, ignoring this contribution is mainly for convenience because we only treat disk thermodynamics approximately by artificially relaxing disk temperature to a target temperature in relatively short timescales through the Λ c term. This approach renders the discrepancy largely irrelevant. More details will be given in Section 2.3.
The above equations are written in c.g.s. units, which will be used consistently in this paper. In the simulations, factors of 4π are absorbed into the definition of B so that magnetic permeability is µ B = 1. The equations are solved in spherical-polar coordinates (r, θ) in 2D with axisymmetry. For convenience, we also use cylindrical coordinates (R, z) in this work to facilitate analysis.

Basic Disk Model
Motivated from observations of PPDs (e.g., Andrews et al. 2009Andrews et al. , 2010, we consider a powerlaw disk surface density profile with exponential cutoff where R AU is radius normalized to 1 AU, R c is the outer radius of the disk beyond which the disk surface density cuts off. We take Σ 0 = 500 g cm −2 with power-law index q S = 1. We are mainly interested in the inner regions of PPDs with R ∼ 1 − 20 AU, where the disk is expected to be largely laminar (Bai 2013), and simply set R c = 30 AU. We note that while Σ should be defined by integrating gas density over the vertical column, in practice, we calculate Σ by integrating over θ at fixed r = R for convenience: it makes little difference as long as the disk is thin. The temperature of the bulk disk is taken to follow the minimum-mass solar nebular scaling (MMSN, Weidenschilling 1977;Hayashi 1981) where we take the power-law index q T = 1/2. Disk temperature sets the isothermal sound speed, given by c s (R) 2 = P/ρ = k B T (R)/µm p , where µ = 2.34 is the mean molecular weight in the bulk disk (molecular gas), k B is the Boltzmann constant, m p is proton mass. The disk scale height is given by where Ω K (R) = GM/R 3 is the Keplerian angular velocity. With q T = 1/2, the disk is flared, with disk aspect ratio ǫ d given by For an MMSN disk, we have T 0 = 280K and ǫ(R = 1AU) = 0.034. In our simulations, we slightly enlarge the disk thickness with ǫ d (R = 1AU) = 0.045 in order to adequately resolve the disk with available computational resources (see Section 2.5 for details). Assuming disk temperature is vertically isothermal, the gas density is given by solving the vertical hydrostatic equilibrium In the radial direction, radial pressure gradient modifies the rotation profile to yield 1

Thermodynamics
The bulk of the PPD is heated by the thermal radiation from the protostar, with temperature approximately given by to Equation (6). The bulk disk is connected to a tenuous disk atmosphere, which is subject to significant heating from higher energy radiation, especially far-UV (FUV) and X-rays, and reaches much higher temperatures (e.g., Glassgold et al. 2004;Walsh et al. 2012). The thermodynamics of the disk atmosphere is much more complex, and detailed modeling would involve UV and X-ray radiative transfer calculations coupled with photochemistry (e.g. Woitke et al. 2016;Haworth et al. 2016). While there are also significant uncertainties, particularly related to the unknown abundance of very small grains, the general result is that there is a rapid temperature transition as high-energy photons are exhausted.
Based on the above, we adopt an approach following the simplified treatment in semi-analytical models of Bai et al. (2016) and Bai (2016) (see Figure 1 in these two papers). The disk is divided into a cold disk zone, and a warm atmosphere, where the division is set by the depth that high-energy radiation (especially FUV) is able to penetrate into the disk. Let Σ FUV be this penetration depth (in units of column density), which is treated as a simulation parameter and is assumed to be a constant. In the simulations, we trace radial rays from the central star in spherical grid to obtain the column densities traversed by the rays, denoted by Σ r (r, θ). For convenience, we express gas temperature in terms of the local disk aspect ratio ǫ ≡ H/R, with 1 We have ignored the exponential surface density cutoff at large R in deriving this equation.
We assume that within the penetration column of FUV radiation, the disk is heated to a temperature that corresponds to a constant local ǫ = ǫ a ≡ 0.2. 2 Beyond this column density, the disk aspect ratio returns to ǫ d given by (7). We then prescribe a smooth but rapid transition to join the two limits, given by (11) To some extent, this treatment is a realization of the semi-analytical global disk evolution model of .
Equations (10) and (11) define the target disk temperature T c (r, θ). In the simulations, we relax the temperature to T c with a simple cooling prescription where the relaxation time τ is set by . (13) Basically, we relax the gas in orbital timescale in the disk zone. This prescription avoids the development of vertical shear instability, which requires fast cooling (Nelson et al. 2013). The relaxation time is gradually reduced to about 1/10 orbital time in the atmosphere, which we find is necessary to sufficiently heat the gas to approach T 0 . We take the adiabatic index γ = 5/3 in our simulations. This is more appropriate for the atomic gas in the atmosphere/wind zone. In reality, the gas in the system transitions to become largely molecular in the bulk disk, leading to an abrupt change in γ at the disk surface. As a caveat, this may result in additional temperature variations in the transition region, and may affect wind kinematics near the wind base. While this is not captured in our treatment, thermal-chemical calculations in hydrostatic disks generally found monotonic increase of temperature from disk to atmosphere (e.g., (Walsh et al. 2012)). Moreover, we do not find appreciable difference in the disk dynamics by using a constant γ = 7/5.
Overall, thermodynamics is treated very approximately in several aspects mentioned previously. In this work, we mainly focus on the role of non-ideal MHD effects on the overall gas dynamics (see the next subsection), which likely play a dominant role governing the overall disk angular momentum transport and flow structure. The simplified treatment of thermodynamics allows us explore these aspects in a more controlled manner. On the other hand, we expect that our treatment of thermodynamics at least captures the most essential ingredients, with Σ FUV being the main controlling parameter.

Ionization, Chemistry and Non-ideal MHD Effects
The strength of non-ideal MHD effects is determined by the ionization degree, or more precisely by the abundance of all charge carriers, a result of disk chemistry 2 In Bai et al. (2016), we find that the wind properties is sensitive to Σ FUV , whereas the detailed temperature structure in the wind zone is of less importance. initiated by the ionization processes. We focus on regions not too close to the protostar so that the disk is not sufficiently hot to trigger thermal ionization ( 800K, Desch & Turner 2015). Below we discuss the nonthermal ionization processes and the calculation of the magnetic diffusivities.

Ionization Rates
The main ionization sources in the bulk disk include cosmic-rays, X-rays. We follow standard prescriptions, where the ionization rates are given as a function of column densities. In addition to tracing radial rays to obtain Σ r (r, θ), we further trace θ−rays from the upper and lower poles towards the disk at constant r, and define two column densities Σ top θ (r, θ) and Σ bot θ (r, θ). While these rays are not straight, the two column densities reach physically meaningful values only close to the bulk disk, where the rays are largely vertical in geometrically thin disks considered here.
For cosmic-ray ionization, the ionization rate is given by (Umebayashi & Nakano 1981) (14) where Σ CR = 96 g cm −2 . We note that the rate of cosmic-ray ionization rate bears large uncertainties (e.g., McCall et al. 2003;Cleeves et al. 2013a), though it would not change the fact that the midplane region of the inner disk is extremely weakly ionized, suppressing the MRI (Bai 2011).
For X-ray ionization, we adopt the fitting formula of Bai & Goodman (2009), based on calculations of Igea & Glassgold (1999). We assume an X-ray luminosity of L X = 10 30 erg s −1 , and use the fitting coefficients at X-ray temperature T X = 3 keV, which gives The first term accounts for the direct absorption of the X-rays along radial rays, with ζ 1 = 6.0 × 10 −11 s −1 , Σ X,a = 3.6 × 10 −3 g cm −2 , and α = 0.4. We note that for the direct absorption component, the Σ X,a value in the original fitting formula corresponds to the vertical instead of radial column density. We here multiply Σ X,a by a geometric factor of 5 to account for the conversion. 3 The second term describes the scattered X-rays that penetrate deeper, with ζ 2 = 1.0 × 10 −14 s −1 , Σ X,s = 1.7g cm −2 , β = 0.65. In addition, we include ionization from the decay of short-lived radionuclides by adding a constant ionization rate of ξ SLR = 6.0 × 10 −19 s −1 .
Note that the rate is expected to be in the range of ∼ (1 − 10) × 10 −19 s −1 and is higher at smaller radii and early times (Umebayashi & Nakano 2009;Cleeves et al. 2013b).

Disk Chemistry and Diffusivity Table
In the bulk disk, ionization-recombination equilibrium is typically achieved within dynamical time (Bai 2011), especially in the presence of grains. In this case, the ionization degree, and hence magnetic diffusivities, are functions of the ionization rate discussed in the previous subsection, as well as gas density and temperature. We make magnetic diffusivity tables from chemistry calculations described below.
We use a complex chemical reaction network developed in Bai & Goodman (2009) and Bai (2011), based on the work of Ilgner & Nelson (2006), with 175 gas-phase species. Since Bai (2014), we extract gas-phase reactions from 2012 UMIST database described in McElroy et al. (2013), adopting the updated rate coefficients from the new database. In total, there are 2147 gas-phase reactions, including four ionization reactions and one ad hoc reaction to account for H 2 formation. 4 In addition, a population of single-sized grains are also included in the network, with maximum grain charge set to ±3. We choose grain size a = 0.1µm and mass fraction of f = 10 −4 , where the total surface area is comparable to those obtained by more realistic grain coagulation/fragmentation calculations (Birnstiel et al. 2011). For a given set of parameters, we start from singleelement species and evolve the network species for 3×10 6 years, which is sufficient to reach chemical equilibrium over a wide range of parameter space. Even the abundance of some species still show secular evolution trends, the ionization fraction converges well before the end of the calculation.
Magnetic diffusivities are calculated following standard formulas (Wardle 2007;Bai 2011). Ohmic resistivity is always independent of magnetic field strength. For Hall and ambipolar diffusivities, we have η H ∝ B and η A ∝ B 2 when field is either weak or very strong. Complex dependence on B may be present at intermediate field strength in the presence of charged grains (Xu & Bai 2016). Nevertheless, as we studied in detail in Xu & Bai (2016), unless small grains are very abundant (say 0.1µm grains with f = 10 −2 ), the weak field limit is always satisfied in practice. Therefore, it suffices to assume η O , η H /B and η A /B 2 , as we adopt in the table.
Further complications in determining the magnetic diffusivities may arise from the non-linear Ohm's law Mori & Okuzumi 2016) that we have not accounted for. Nevertheless, the theory has only been worked out for Ohmic resistivity, which operates in the presence of very strong current, while it is less likely to be relevant in disks that are largely laminar, as we have in this work. 4 We have made two changes in the chemistry calculations compared with our previous work (e.g., Bai 2011Bai , 2014. First, in case when temperature falls out of the range of validity of the fitting formula provided in in the data base, we now still use the fitting formula as if it remained valid. See Section 3.1 of Xu & Bai (2016) for further explanation. Second, for electron-grain collisions, we now adopt a constant sticking coefficient of se = 0.3, instead of directly calculating se described in the Appendix of Bai (2011). This change follows from Ivlev et al. (2016), see also Weingartner & Draine (2001).

FUV Ionization in the Disk Atmosphere
High-energy radiation, not only heat the disk atmosphere, but also significantly boost its ionization level. In particular, FUV can fully ionize atomic carbon and sulfer, raising the ionization fraction to x e ≡ (n e /n) = 10 −5 to 10 −4 (Perez-Becker & Chiang 2011), with a sharp transition in x e at the FUV ionization front. Here, n e and n are the number densities of the electrons and neutrals, respectively. Note that our chemical network does not include photo-reactions to account for FUV ionization, but we mimic this effect by setting the ionization fraction due to FUV to be The resulting non-ideal MHD diffusion coefficients are evaluated according to (applicable when x e ≫ charged grain abundance) where m n = µm H is the mean molecular mass, σv e ≈ 8.3 × 10 −9 (T /100K) 1/2 cm 3 s −1 , and σv i ≈ 2.0 × 10 −9 cm 3 s −1 are coefficients of momentum exchange in electron-neutral and ion-neutral collisions (Draine 2011). We calculate the diffusivities both from the diffusivity table, as well as from (18) with x e given by (17), and set the diffusivities to be the ones with smaller values, which guarantees a smooth transition from the disk zone to the atmosphere.
In practice, we further boost x e,FUV by a factor g given by The sole purpose of this factor is to ensure that the gas behave in the ideal MHD regime throughout the wind zone (otherwise AD would become progressively more important, while this is not the case from more realistic calculations as shown in Figure 9 of Walsh et al. 2012). The value of Σ FUV is uncertain, and is particularly sensitive to the abundance of very small grains. Perez-Becker & Chiang (2011) found Σ FUV ∼ 0.01 − 0.1 g cm −2 in their 1D calculations. The value they quote has been converted to a vertical column density with a geometric factor of 0.3. Since we attenuate FUV along radial rays, we expect a range of Σ FUV ∼ 0.03 − 0.3 g cm −2 to be more appropriate, and we take Σ FUV = 0.03 g cm −2 as fiducial.

Simulation Setup
We set the radial grid to span from r in = 0.6 AU to r out = 60 AU with logarithmic grid spacing. The θ−grid extends from the midplane all the way to near the poles (leaving only a 2 • cone at each pole), which we find is essential to properly accommodate the MHD disk wind. The θ−grid is concentrated around the disk which guarantees adequate resolution, where ∆θ increases by a constant factor per cell from midplane to pole, with a con-trasting factor of 3.5 between minimum and maximum ∆θ. The full grid size is 1152 × 512 in (r, θ), so that at R = 3 AU where disk aspect ratio ǫ d ≈ 0.06, we achieve a grid resolution of 15 cells per H d in r, and 20 cells per H d in θ.
We initialize the disk from the hydrostatic solution (8) and (9), where disk temperature is set to be vertically isothermal with ǫ(r, θ) = ǫ d (R). A density floor of ρ(r) = 10 −8 ρ mid (r) is set to prevent excessive density drop in the disk atmosphere in the initial condition. This density floor is sufficiently small so that upon achieving quasisteady state, the density in the disk atmosphere is well above the floor value thanks to disk wind launching.
In the simulations, we trace radial rays and θ−rays to compute Σ r (r, θ) and Σ top,bot θ (r, θ). This operation involves global communications, and is executed only at a time interval of 0.5Ω −1 K (r in ) (corresponding to at least a few thousand timesteps) so that it does not affect the overall code performance.
We first run the simulations for 300Ω −1 K (r in ) without magnetic fields, which allows the disk atmosphere to be heated to desired temperature according to (11), and the gas density to adjust to a new equilibrium state where analytical solutions are not available. Right afterwards, we apply an external poloidal field using a vector potential generalized from Zanni et al. (2007) where m is a parameter that controls how much poloidal fields are bent, with m → ∞ giving a pure vertical field. We choose m = 1, and we have verified that the results are insensitive to the choice of m. Poloidal fields are obtained from B = ∇ × (A φφ ) in a way that guarantees ∇ · B = 0. At the midplane, we have B = B z0ẑ (r/R 0 ) −(α+qT )/2 , maintaining constant ratio of gas to magnetic pressure, defined by plasma β 0 . Fiducially, we choose β 0 = 10 5 , appropriate for the inner region of PPDs (Bai & Stone 2013;Bai 2013), but we also consider stronger fields in Section 7.1.
In simulations including the Hall effect with aligned field polarity, we find that the outcome of the simulation depends on initial conditions. This issue is discussed in more detail in Appendix A. For simulations shown in the main text, we have modified the simulation setup to obtain more consistent result (see Section 5 for details).
An important aspect of the simulations is to properly control the inner boundary condition. We note that the flow near the polar region should originate from a part of the disk located inside the inner boundary, whose dynamics is beyond the reach of the simulation. A standard outflow-type boundary condition would violate causality, and can become unstable in the presence of magnetic fields which further interfere with the wind flow in the main computational domain. After experimenting with a number of options, we adopt a boundary condition that is close to a fixed state, which we find can better constrain the flow structure near the inner boundary and minimize its effect to the rest of the simulation domain. We fix the density profile as (8) with a density floor, v r = v θ = 0, and v φ is set to be the minimum of the initial v φ (9) and AU . The simulation domain extends from 0.6-60 AU in radius.
Ω K (r in )R. Gas temperature is set based on (11). We also set a buffer zone between r = r in and r = 1.5r in , where we linearly reduce all magnetic diffusivities with radius to zero, and damp gas poloidal velocities on local orbital timescale. This approach prevents accretion into the inner boundary and would lead to some mass accumulation within the buffer zone. We then deplete the gas in the buffer zone over a timescale of 10 4 Ω −1 K (r in ). Also note that the duration of our simulations is relatively short, and we do not observe significant modifications of the inner disk structure.
The rest of the boundary conditions are straightforward. The outer radial boundary follows from standard outflow boundary prescriptions, where hydrodynamic variables are copied from the last grid zone assuming ρ ∝ r −2 , v φ ∝ r −1/2 , with v r and v θ unchanged except that we set v r = 0 in case of inflow. Magnetic variables in the inner/outer ghost zones are copied from the nearest grid zone assuming B r ∝ r −2 and B φ ∝ r −1 , with B θ unchanged. Reflection boundary conditions are applied in the θ−boundaries.

List of Runs and Parameters
We list all our simulation runs in Table 1. Most parameters are fixed as described in previous subsections, and we only vary two parameters: disk magnetization (parameterized by plasma β 0 , fiducially 10 5 ), and FUV penetration depth Σ FUV (fiducially 0.03g cm −2 ). Time is measured in unites of Ω −1 0 ≡ Ω K (r in ) −1 in our simulations. For r in = 0.6 AU, we have Ω −1 0 = 0.074 yr. We will focus on our fiducial runs, labeled as "Fid±", where the +/− signs correspond to simulations with poloidal field aligned/anti-aligned with disk rotation. For comparison, we also conduct a run "Fid0", where we turn off the Hall effect. These simulations are run for more than 2000 orbits at innermost radius for detailed analysis. We then vary one parameter at a time, and for each variation, three runs labeled by "±0" are performed as in the fiducial case. They are run for shorter amount of time but are sufficient to illustrate the dominant features.

DIAGNOSTICS
In this section, we discuss major diagnostics to be employed to analyze our simulation results.

Elsasser Numbers
The strength of the non-ideal MHD effects are conveniently measured by dimensionless Elsasser numbers., defined as where v A ≡ B/ √ 4πρ is the Alfvén speed. Non-ideal MHD terms are considered strong if the Elsasser numbers are of order unity of less. For Ohmic resistivity, Λ < 1 is generally sufficient to suppress the MRI (Turner et al. 2007;Ilgner & Nelson 2008). For AD, Am < 1 can suppress or damp the MRI depending on vertical field strength (Bai & Stone 2011). Since only Am is independent of field strength. For the Hall effect, a field strength independent measure is the Hall length l H , defined as which is the generalization of ion inertial length in weakly ionized plasmas (Kunz & Lesur 2013). Strong Hall effect is characterized by l H H.

Angular Momentum Transport and Disk Flow
Structure Angular momentum transport is mainly mediated by magnetic stresses. Ignoring hydrodynamic processes, the equation of angular momentum transport in cylindrical coordinates can be written aṡ whereṀ acc ≡ −2πR are Maxwell stresses, and ±z b mark the vertical coordinates that separate the disk and atmosphere. Overlines represent time and azimuthal averages, and we have assumed Keplerian rotation.
Physically, the first term on the right hand side corresponds to radial transport of angular momentum. In our simulations, the disks are largely laminar, and this term is dominated by large-scale fields that wind up into spirals, corresponding to magnetic braking. By convention, we define the Shakura & Sunyaev (1973) which is dimensionless measure of the stress. Note that the accretion rate is related to the radial gradient of T Rφ . The second term on the right hand side corresponds to vertical transport of angular momentum by magnetized disk winds. We can normalize T zφ by the midplane gas pressure. We note that given similar field strength, vertical transport is more efficient than radial transport by a factor of ∼ R/H (Wardle 2007;Bai & Goodman 2009). The B z B φ stress drives accretion by exerting a torque on the disk, and the torque density is proportional to its vertical gradient. Because B z ≈ constant in a thin disk, the wind-driven local accretion velocity is given by We see that the accretion mass flux ρv R is directly proportional to toroidal field gradient. This is the most important relation for understanding the flow structure in our simulations.

Wind Kinematics
In steady state and axisymmetry, a magnetized disk wind is characterized by a series of conservation laws along field lines (e.g., Spruit 1996), as long as the gas is well coupled with magnetic field (i.e., ideal MHD). These include the conservation of mass angular velocity of magnetic flux surface and specific angular momentum subscript p denotes the poloidal component. Here, k, ω and l are conserved along poloidal field lines, and Ω ≡ v φ /R. In practice, we normalize these quantities to k 0 = 4πρ mid v K /B z0 , Ω K , and Ω K R 2 0 , where v K , Ω K and B z0 are defined at the wind launching radius R 0 . If the equation of state is barotropic (ours is not), energy conservation can also be expressed explicitly. We will test the above three relations in our simulations.
The most important wind diagnostics is the mass loss rate. We defineṀ wind (R) as the cumulative wind mass loss rate within radius R. Locally, we quote the mass loss rate per logarithmic radius as An important concept is the Alfvén radius R A , where for the wind flow originating from radius R 0 , R A is the radius of the point along the field line where the poloidal flow velocity v p equals to the poloidal Alfvén velocity v Ap = B p / √ 4πρ. The local wind mass loss rate is closely related to wind-driven accretion rate by (Ferreira & Pelletier 1995;Bai et al. 2016) Therefore, the location of the Alfvén point provides an alternative and more convenient measure of the wind mass loss rate. The ratio λ ≡ (R A /R 0 ) 2 is defined as the magnetic lever arm.

Magnetic Flux Transport
The basic physics of magnetic flux transport due to non-ideal MHD effects has been studied in detail in Bai & Stone (2017), and we do not pursue analysis in . Left three panels: snapshots of magnetic field configuration represented by equallyspaced contours of poloidal magnetic flux and color (scaled toroidal field RB φ ) at t =1200, 3000 and 15000 Ω −1 0 . Rightmost panel: radial mass flux ρvr (rescaled by r q D R 1/2 at the last snapshot (t=15000Ω −1 0 = 1109 yrs), overlaid with poloidal flux contours (black) and velocity vectors (green arrows). Thick white dashed lines in the four panels mark the FUV ionization front, and the magneta contours in the left three panels mark the Alfvén surface. Note that the simulation domain extends to r = 60 AU.
as much detail as was done there. In most occasions, we simply measure the total magnetic flux enclosed within radius r at the midplane and follow its evolution.

SIMULATION
We start from the Hall-free simulation before introducing further complications owing to the Hall effect. In Figure 1, we show snapshots of magnetic field configurations from run Fid0. Overall, the system quickly settles into a laminar configuration in approximately steady state in 10-15 local orbital time, launching a disk wind. To a large extent, the system is symmetric about the midplane (except for regions between R ∼ 2 − 5 AU). We will further discuss the field configuration in Section 4.2.
The quasi-steady state configuration allows us to choose two characteristic radii, 2 AU and 10 AU, and analyze the overall gas dynamics in further detail. In Figures 2 and 3, we show the vertical profiles of main diagnostic quantities at the last simulation snapshot (t = 15000Ω −1 0 ≈1109 yrs), including density, temperature, Elsasser numbers, pressure, Maxwell stress, magnetic fields and velocity fields. The results are discussed in Sections 4.1, 4.2 and 4.3 from different perspectives. Section 4.3 further addresses the overall disk angular momentum transport and examines the disk flow structure. We analyze wind kinematics in Section 4.4.

Magnetic Diffusivities
From Figure 2, we see that at both R = 2 and 10 AU, the FUV front is located around z = ±5H d . Our thermodynamic scheme nicely maintains constant disk aspect ratio ǫ = ǫ d (R) within the disk zone, and allows it to smoothly rise to ǫ ∼ 0.2 beyond the FUV front. The density decreases with height much more slowly near and beyond the the FUV front, owing to magnetic pressure support (see Figure 3).
The right panels of Figure 2 show the Elsasser number profiles. Note that while we have switched off the Hall effect, we can still show the Hall Elsasser numbers and the Hall length l H . Also, l H and the AD Elsasser number Am profiles, being independent of field strength, remain largely unchanged over the course of the simulation. This also holds in simulations with the Hall term turned on (we thus do not repeat similar plots in later discussions).
We first focus on the ambipolar Elsasser number Am. The FUV front leads to the most sharp increase of Am to near ∼ 100, making ideal MHD a good approximation in the wind zone. Below the FUV front, the Am profiles display a number of wiggles, corresponding to contributions from individual ionization sources (direct X-ray absorption, X-ray scattering, cosmic-rays).
At 2 AU, resistivity and the Hall effect are the two dominant non-ideal MHD effects at the midplane, with Elsasser number orders of magnitudes below 1. This corresponds to the conventional Ohmic dead zone (Gammie 1996), with extremely weak level of ionization with magnetic field largely decoupled with the gas. The Hall effect dominates between z ∼ 1 − 3H d , and AD takes over to dominate at disk upper layers. At 10 AU, Ohmic resistivity becomes largely irrelevant in the entire vertical column. The Hall effect dominates the midplane region, and AD dominates beyond z ∼ ±2H d . Note that Hall length l H well exceeds H d at R = 2 AU in the midplane region, whereas it drops below H d at 10 AU, reflecting that the Hall effect weakens towards larger disk radii.
The Elsasser number profiles obtained here are largely consistent with those obtained in local simulations for the inner disk using similar parameters (Bai & Stone 2013; Bai 2013). We also note a flat Am profile ∼ 1 at 10 AU within the FUV front. This fact also approximately holds towards larger radii.

Magnetic Field Configuration
The vertical structure of magnetic fields in the disk can be seen from Figure 3. Overall, the field configuration is consistent with previous (vertically-)local simulations (e.g., Figure 11 of Bai & Stone 2013, Figure 6 of Gressel et al. 2015).
Initially, the poloidal fields bend radially outwards, 5 generating oppositely-directed toroidal fields above/below the midplane, which build up magnetic pressure. As seen from the Figure, magnetic pressure dominates over gas pressure beyond about z = ±4H d . The vertical gradient of toroidal field exerts a torque to the disk, driving radial flows according to (25). The radial flow is an accretion flow within a certain height about the disk midplane before |B φ | reaches its maxima. Beyond the maxima, |B φ | slowly decreases with height (driving a weak decretion flow), and its associated magnetic pressure gradient directly drives the disk wind (see Section 4.4).
Toroidal field is always the dominant field component.
In steady state, the generation of toroidal field (from radial field) is mainly balanced by Ohmic/ambipolar dissipation in the bulk disk, whereas in the wind zone, it is simply balanced by advection in the wind flow. The poloidal field configuration, especially the level field lines bend, is mainly set by the radial flow structure through the vertical extent of the disk (which transitions from accretion to decretion towards surface, and drives field line bending), Ohmic/ambipolar dissipation (which works to straighten the field in the bulk disk), as well as outflow advection (towards the surface).

The Symmetry Issue
Launching MHD wind with physical symmetry requires toroidal field to change sign across the disk, which is directly connected to driving the accretion flow (Bai & Stone 2013). This is achieved in different ways at small and large radii.
At small radii (e.g., R 2 AU), the extremely weak level of ionization means that there are very few charge carriers to sustain current, which minimizes the vertical gradient of toroidal field. In Bai & Stone (2013), we found that toroidal field of one sign reaches maxima has a flat profile across the midplane, it then flips sharply at a few scale heights at the other side the midplane (as the gas becomes better coupled to the field). This is observed in early stages of evolution (i.e., second panel in Figure  1). Later, on the other hand, midplane toroidal field decreases, and toroidal field peaks at a similar heights both above and below the midplane (i.e., third panel in Figure  1, and central top panel of Figure 3). These are similar to the "belt" structures observed in Gressel et al. (2015).
At large radii (e.g., R 5 AU), the midplane region becomes better coupled with magnetic field (Elsasser number ∼ 1), allowing the flip of toroidal field to take place right at the midplane. This leads to reflection symmetry about the midplane, as seen in Figure 1, and the central bottom panel of Figure 3. This result is consistent to earlier local studies of Bai (2013).
In between R ∼ 2 − 5 AU, the midplane region of the disk is dominated by patches of toroidal field in opposite signs. These patches interact with each other, leading to some secular evolution of the system. We do not find signatures of unstable MRI channel modes develop, as in some runs in Gressel et al. (2015). Instead, we simply interpret these phenomena as inherent to the transition from the extremely poorly coupled regime at small radii to the marginally coupled regime at larger radii. 6

Angular Momentum Transport and Flow Structure
Magnetized disk wind is the dominant mechanism of disk angular momentum transport. The vertical distribution of wind-driven accretion flow directly results from toroidal field gradient, shown in Equation (25). This is verified in the top panels of Figure 4 for R = 2 and 10 AU, respectively. At R = 2 AU, there are mainly two accretion layers above and below the midplane at around z = ±1.5H d , where maximum toroidal field gradient develops, with slightly different mass fluxes (owing to slight   The radial flow velocity can be found in the insets on the right panels of Figure 3. The thickness of the accreting layer is typically a good fraction of a scale height, and is well resolved in the simulation. We note that beyond the accreting layer towards the surface, the gas is directed radially outward (i.e., decretion) because the toroidal field gradient reverses. Nevertheless, combined with the rapid density drop, the surface layer decretion carries a negligible fraction of mass flux, as seen in Figure  4.
The transition to ideal MHD regime at the FUV front shows distinct features in magnetic field and flow structures, as can be found in Figure 3, and they are mostly consistent with previous local simulations (Bai & Stone 2013;Bai 2013). In particular, we have defined the wind base as the location where the gas flow transitions from being sub-Keplerian (below) to super-Keplerian (above) (Wardle & Koenigl 1993), and the location is found to largely coincide with the FUV front (Bai & Stone 2013;Gressel et al. 2015). However, in our global simulations, we find that such transition never takes place. Instead, v φ almost never exceeds Keplerian, and simply decreases towards larger height. This is related to the nature of magneto-thermal disk winds Bai et al. (2016), to be fur-ther discussed in Section 4.4. On the other hand, the FUV front does correspond to a local maxima in v φ , which may still be considered as a reasonable way to define the wind base z b , as we adopt here. Similar situation holds for simulations with the Hall effect 7 .
The bottom panel of Figure 4 shows the radial profiles of accretion and mass loss rates, together with accretion rates computed from (23), separating the contribution from wind-driven transport from radial transport of angular momentum. Overall, the mass accretion rate is approximately constant over radius, and is around 2 × 10 −8 M yr −1 . Note that over the course of our simulation, there has been very little evolution of gas surface density. The flat accretion rate profile is largely a result of proper choice of the initial magnetic flux distribution (i.e., constant plasma β 0 ).
We again see that disk wind accounts for almost the entire accretion process in the disk, whereas radial transport of angular momentum due to T Rφ from the Maxwell stress is about one order of magnitude smaller. The Shakura-Sunyaev α measured from our simulation is small, ranging from α ∼ 10 −4 at R ∼ 1 AU, to α ∼ 4 × 10 −4 beyond R = 5 AU. Such small value can already be inferred from the left panels of Figure 3.
At intermediate radii between ∼ 2 − 5 AU, the sys-   (25). The horizontal dashed line marks zero mass flux to guide the eye. Bottom panel: radial profile of mass accretion rate (black), mass loss rate per logarithmic radius (blue solid), cumulative mass loss rate (blue dashed), predicted wind-driven accretion rate (red solid), and the predicted viscouslydriven accretion rate (red dashed). All calculated from the last snapshot of the simulation. tem exhibits slightly larger accretion rates, which is related to enhanced magnetic activities associated with the secular evolution of toroidal field patches. We also see from the rightmost panel of Figure 1 that the disk gas shows complex radial flow structures. In these regions, (23) is no longer applicable to predict accretion rates due to time variability, but overall, the accretion flow still largely correlates with strong toroidal field gradients at the boundaries of the toroidal field patches, where most of the torque is exerted. Figure 4 further demonstrates that wind mass loss rate, calculated from (29), is excessive. The mass loss rate per logarithmic radii already exceeds the accretion rate, whereas the cumulative mass loss rate reaches as large as 4 times the accretion rate at R = 20 AU. This is again a consequence of magneto-thermal disk winds, which we focus on in this subsection.

Wind Kinematics
To analyze the wind kinematics, we again choose R 0 = 2 and 10 AU, and trace poloidal field lines from the midplane all the way to the boundary of our simulation domain. We measure various diagnostic quantities along the field lines and show the results in Figure 5.
We see that the conservation laws outlined in Section 3.3 are generally well satisfied beyond the FUV front where the gas behaves approximately in the ideal MHD regime. The mass flux along field lines is substantial. With k ∼ 10 −5 k 0 and H/R ∼ 0.05 − 0.08 (at R = 2 − 10 AU), the corresponding local disk depletion timescale is typically only a few thousand orbits. We also find that ω ≈ Ω K (R 0 ) as the angular velocity of magnetic flux surface. The value of l 2Ω K (R 0 )R 2 0 , indicating that the wind carries less than twice the specific angular momentum at the origin, which is consistent with the fact that local wind mass loss rate is comparable to winddriven accretion rate.
Poloidal velocity constantly increases along the field line within the simulation domain. We see in Figure 1 that the Alfvén surface is located relatively close to the disk, and is only slightly beyond the FUV front. The fast magnetosonic point, defined when poloidal velocity equals to the poloidal fast magnetosonic velocity v 2 , is not contained in the simulation domain, as seen from Figure 5. Note that the fast point is typically at very large distances in wind theory . In practice, we have tested that containing fast magnetosonic point is not crucial to wind kinematics 8 By field line tracing at R 0 = 2 AU and 10 AU, we find R A /R 0 ≈ 1.36 − 1.40. According to Equation (30), this implies that dṀ wind /d ln R ≈ 0.52 − 0.59. On the other hand, the actual mass loss rate appears to be a factor of ∼ 3 higher, as seen from Figure 4 This apparent discrepancy will be discussed in Section 4.4.1.
The heavily loaded wind in our simulations is mainly driven by the toroidal magnetic pressure gradient. To show this, we decompose the poloidal forces following Section 3.1.1 of Bai et al. (2016) where the three terms correspond to thermal pressure gradient, the net centrifugal force, and the Lorentz force from toroidal magnetic pressure gradient. Note that we define the net centrifugal force as the excess of centrifugal force over gravitational acceleration. We see in Figure 6 that this force is always negative, meaning that corotation is far from being enforced to drive centrifugal fling, as in the conventional Blandford & Payne (1982) picture. Instead, acceleration is dominated by magnetic pressure gradient. This is because poloidal fields in PPD winds are too weak to enforce corotation. They are thus wound up by differential rotation, developing strong toroidal 8 While working on the semi-analytical wind model of Bai et al. (2016), we have also solved time-dependent MHD wind equations along prescribed poloidal field lines (unpublished). We find that the Alfvén point is quickly settled even the flow is far from reaching the fast magnetosonic point, and the Alfvén radius established early on is almost identical with the final steady-state solution.  fields. We also see from Figure 5 that v φ falls off approximately as R −1 , and that |B φ |/B p well exceeds 1 in the wind zone. These results are all consistent with the conclusions in Bai et al. (2016).

Physics Behind Excessive Wind Mass Loss
From Figure 4, we estimate dṀ wind /Ṁ acc ∼ 1.5, which translates to R A /R 0 ∼ 1.15 based on Equation (30), indicating extremely small lever arm. On the other hand, as mentioned earlier, when taking R 0 as the radius of the field origin at the midplane, we find R A /R 0 ∼ 1.4. This apparent inconsistency is resolved by noting that when computing the lever arm, R 0 should be defined as the radius of the wind base. More appropriately, it should be taken to be the radius of the FUV front R FUV (see Section 4.3). In fact, from Figure 5, we find exactly R A /R FUV ≈ 1.15 at both R 0 = 2 AU and R 0 = 10 AU.
Another apparent inconsistency arises when comparing the mass loss rate with semi-analytical theory of Bai et al. (2016). In our fiducial run, the poloidal Alfvén velocity and the sound speed at the wind base are found to be around 0.1 − 0.15v K (R FUV ), very close to the fiducial parameter values adopted in Bai et al. (2016) Figure 1, but for run Fid+ with all three non-ideal MHD terms included and vertical field aligned with disk rotation. Note that this run is restarted from the Hall-free run Fid0 at t = 3000Ω −1 0 ∼ 222 yr, after which the Hall term is turned on in an inside-out manner over a period of 5 local orbits. For reference, at R = 10 and 20 AU, the Hall term is fully included after 380 yrs and 669 yrs. See Section 5 for more details.

(taken
to be 0.1v K ). However, although it was pointed out there that PPD wind is heavily loaded, the predicted mass loss rate is about an order of magnitude smaller than measured in our simulation.
Two factors contribute to the excessive mass loss rate in our simulation. First, poloidal field strength B p appears to drop at a rate comparable to or faster than R −2 , as seen from the bottom right panel of Figure 5. In Bai et al. (2016), it was found that the wind lever arm is sensitive to how rapidly B p decreases with R. The fiducial model adopted there assumes B p ∝ R −1 near the disk, and transitions to B p ∝ R −2 at larger distances, giving R A /R 0 ∼ 2.3. On the other hand, assuming B p ∝ R −2 dramatically reduces the lever arm with R A /R 0 ∼ 1.5. Our simulation result suggests even faster decrease of B p , which is likely related to the fact that consecutive poloidal field lines are collimated at different levels, where field lines originating from smaller radii are more collimated.
The second, and probably more important factor lies in the angular velocity of magnetic flux surface ω. We have shown that ω ≈ Ω K (R 0 ). However, it is more appropriate to normalize ω to Ω K (R FUV ) (i.e., at the wind base). For the two field lines shown in Figure 5, we find ω ≈ 1.25 − 1.3Ω K (R FUV ). This is substantially larger than the range of values considered in Bai et al. (2016), who considered the ratio in the range of 0.95 − 1.05. Despite the offset, a clear trend was identified that higher ω leads to heavier mass loading (see their Figure 10). 9 We may further ask why ω deviates substantially from Ω K at the wind base. The reason is that the field is still largely anchored to the disk at radius R 0 , thus rotating at ∼ Ω K (R 0 ). In the disk upper layers (before reaching the FUV front) where X-rays are the dominant ionization 9 While we can repeat the calculations done in Bai et al. (2016) using higher ω, we end up violating the assumptions made there: the slow magnetosonic point is found to be well within the FUV front (where ideal MHD no longer applies). source, the coupling between gas and field is marginal, giving Am ∼ 1. Such marginal coupling allows poloidal fields to bend, reaching the FUV front at a larger radius R FUV . This makes the field lines rotate faster than the local Keplerian speed, enhancing wind mass loading.

FIDUCIAL SIMULATION IN THE ALIGNED CASE
Following the discussion in Appendix A, for simulations including the Hall effect with aligned field polarity, we restart from simulations from the Hall-free run at a certain time t 0 , and then turn on the Hall term in an inside-out manner. This is motivated from evolutionary (i.e., disk formation) considerations detailed in Appendix A. More specifically, we choose t 0 = 3000Ω −1 0 (∼ 222 yr) 10 , and turn on the Hall term according to [(t − t 0 )/t H ] 4 for t < t 0 + t H , where t H is set to 5 local orbits. After t = t 0 + t H , the Hall term is completely included. The simulation is run for sufficient amount time, measuring more than 10 local orbits at R = 20 AU.

Overall Evolution and Magnetic Field
Configuration In Figure 7, we show the time evolution of magnetic field configuration from run Fid+. Once the Hall term is turned on, horizontal components of the magnetic field are quickly amplified due to the Hall-shear instability (HSI, Kunz 2008;Lesur et al. 2014;Bai 2014). An MHD disk wind is always launched in the presence of net poloidal field, maintaining a physical geometry that poloidal field lines bend radially outward on both sides of the disk. This physical field geometry requires toroidal field to change sign across the disk (Bai & Stone 2013).
In the presence of HSI, we see that the disk can be clearly divided into two zones. The inner zone shows The above phenomena are closely related to the development of the HSI. Globally, the HSI is associated with radial transport of poloidal magnetic flux along the direction of the Hall drift (or electron-ion drift in a dustfree gas), which stretches poloidal field lines towards a radially-elongated configuration. The radial field is then sheared by differential rotation in the disk to further amplify the toroidal field (Bai & Stone 2017). The direction of the Hall-drift is along the electric current, which is mainly due to the vertical gradient of toroidal field, whose growth feeds back to the HSI. In the inner zone, we see that poloidal field lines are highly inclined in the midplane region, which is the source of strong toroidal field there. Similarly, poloidal field lines in the outer zone also have a significant radial component, which strongly contrasts with the Hall-free run Fid0. The flip of toroidal field is associated with a kink in poloidal field.
Generally speaking, the transition from the inner to the outer zone is related to the transition from Halldominated to AD-dominated regime in the bulk disk. This transition can already be traced from the Elsasser profiles shown in Figure 2 for run Fid0. For run Fid+ we find that around 12-20 AU, the Hall term and AD have comparable strength in the midplane region, whereas the strength of the Hall term drops rapidly towards the surface. In practice, there are more subtle issues, which are further discussed in Appendix B (e.g., comparison with the simulations of Bai & Stone 2017 which show fully symmetric solutions). We also notice that in later stages of our simulation, there is a segregation of magnetic flux in the inner few AU of the simulation box. This is partly related to the asymmetry in the inner zone: poloidal field lines are highly inclined, and hence for the same field line, it reaches the disk surface (where the wind is launched) at different radii in the upper and lower sides of the disk. However, this asymmetry is broken by the presence of the inner boundary. Once a field line is attached to the inner boundary, it loses disk support and becomes isolated. In reality, however, it should be connected to some part of the disk inside the inner boundary. In our simulations, we find that the inner boundary gradually attracts magnetic flux from the lower side of the disk where the pinched poloidal field line first reaches the inner boundary, and then the entire field line is accreted to the inner boundary, building up magnetic flux there. While this flux segregation phenomenon might be real to a certain extent, it is clearly affected by the inner boundary, and hence is not very trustable.
In the rest of the discussion, we analyze the result at the end of our simulation, focusing on regions characteristic of the asymmetric inner zone (∼ 8 − 12 AU), and the more symmetric outer zone (∼ 16 − 20 AU). We avoid regions affected by the flux segregation phenomenon (within a few AU).

Angular Momentum Transport and Flow Structure
In this subsection, we discuss the mechanisms that drive disk angular momentum transport in run Fid+, and the associated flow structure. In Figure 8, we show the vertical profiles of major diagnostic quantities at two representative radii R = 10 AU and 18 AU for the inner and outer zones. We first discuss the results at these two locations separately before analyzing the global diagnostics.

Vertical Structure at the Inner Zone
The inner zone is characterized by a strong toroidal magnetic field B φ whose strength peaks at the midplane. Due to the HSI, we see from Figure 8 that B φ is amplified to ∼ 60 times the net vertical field. Note that even with such significant amplification, midplane magnetic pressure only reaches a few percent of gas pressure (plasma β ∼ 30). This is consistent with previous local simulations (Bai 2014(Bai , 2015, and for given vertical field, the amplification factor also depends on the assumed grain abundance (e.g., see Lesur et al. 2014;Simon et al. 2015;Xu & Bai 2016 where different grain abundances are adopted).
The vertical profile of B φ is approximately symmetric about the midplane before changing sign at one side of the disk surface. This configuration largely determines the vertical profiles of the bulk flow velocity according to Equation (25). In the left panel of Figure 9, we find excellent agreement between the radial mass flux measured in the simulation and expectation from Equation (25).
To elaborate, as the radial mass flux is largely determined by the vertical gradient of B φ , it becomes clear that the bulk gas flows inward at one side of the disk, and outward at the other side. Additional radial mass flux occurs at the position where B φ flips(at one side of the disk surface). The absolute mass fluxes in the three regions are comparable (given the B φ profile). On the other hand, because the midplane region is much denser than the surface, radial flow velocity there is typically small (as seen in the top right panel of Figure 8), on the order of 1% of the sound speed. The flow velocity in the surface, on the other hand, is very significant, and is on the order of the sound speed.
We note that the net wind-driven accretion rate is determined by the difference in the wind stress −B z B φ at the top and bottom disk surfaces. With B z approximately constant within the disk, net accretion rate is largely set by B φ at the surface. While |B φ | is a factor of several stronger in the midplane, it does not yield additional mass flux, but leads to the radial inflow-outflow pattern whose mass fluxes largely cancel each other. The net accretion flux mainly results from the surface layer where B φ flips.

Vertical Structure at the Outer Zone
The outer zone is characterized by an approximately symmetric field configuration where B φ flips at the midplane. Significant amplification of B φ due to the HSI still occurs near the midplane, with amplification factor of up to ∼ 40 as seen in bottom middle panel of Figure 8 at the radius of 18 AU, corresponding to a plasma β of ∼ 100 at midplane. The amplification factor will decrease towards larger radii as the Hall effect gets weaker.
The vertical profile of B φ again leads to a very unusual radial flow structure, mainly determined from Equation (25), as illustrated in the right panel of Figure 9. More specifically, the accretion flow is concentrated at the midplane strong current layer where B φ flips. On the other hand, because |B φ | maximizes right outside of the strong current layer, the drop in |B φ | leads to radial outflows in these regions, both above and below the midplane. Because these regions are not far from the (dense) midplane, both outflow velocities, as well as the midplane accretion flow velocity, are relatively small. At the radius of R = 18 AU, they are on the order of ∼ 1% of the sound speed, as can be traced from the bottom right panel of Figure 8, where c s /v K ≈ 0.09.
The radial outflows above and below the midplane partially cancel the midplane accretion mass flux. The net wind-driven accretion rate is again determined by the difference in the wind stress −B z B φ at the top and bottom disk surfaces, which amounts to roughly 20 − 30% of the midplane accretion rate at R = 18 AU.

Radial Profiles of Accretion Rates
To further analyze the mechanism of angular momentum transport, we show in Figure 10 the radial profiles of accretion and outflow rates. We again treat the FUV front, marked as black dashed lines in Figure 7, as the wind base that separates the bulk disk and the wind zone. We further show accretion rates derived from Equation (23), separating contributions from the disk wind and the laminar Maxwell stress.
We first analyze wind-driven accretion. As noted earlier, towards the end of the simulation, there is a deficit Note that mass outflow rates from the top and bottom sides of the disk are generally different. Middle: radial profiles of mass accretion rates derived from Equation (23), separating contributions from the wind (black) and the laminar Maxwell stress (viscously-driven, solid for accretion, dashed for "decretion".). Right: radial profiles of normalized wind stress |T zφ |/P mid and Shakura-Sunyaev α parameter [from Equation (24) of magnetic flux in the first few AU of the domain, leading to very low wind-driven accretion rates. We discard this region since it is likely related to the limitations of inner boundary conditions. Wind-driven accretion rates reaches about 2 × 10 −8 M yr −1 beyond about 10 AU, with normalized wind stress several times of 10 −4 . This is comparable and slightly larger than the Hall-free case, which is also consistent with local shearing-box simulation result (Bai 2014(Bai , 2015.
The total accretion rate as seen in the left panel of Figure 10, however, differ significantly from winddriven accretion rate. This is largely owing to contributions from viscously-driven accretion from the laminar Maxwell stress, as we discuss below.
From the rightmost panel of Figure 10, we see that the α value peaks at around 7−10 AU, which is right outside the region deficit of magnetic flux. Towards larger radii, the α value decreases, which is related to the fact that the Hall effect weakens. The typical α value reaches a few times 10 −3 , which is consistent with previous local simu-lations in this region (Bai 2015;Simon et al. 2015). Such α values (greater than T zφ /P mid by a factor of several) already suggests that viscously-driven accretion rate can be significant compared with wind-driven accretion rate, as discussed in previous local simulations Bai 2014).
We note that viscously-driven accretion rateṀ acc,V depends on the radial gradient of α and other disk properties. For our adopted radial surface density profile, we haveṀ acc,V ∝ d(αΣT R 2 )/dR ∝ d(αR 1/2 )/dR. Therefore, if α decreases with radius more rapidly than R −1/2 , the laminar stress would drive an "decretion" flow instead of accretion. This appears to be the case beyond ∼ 10 AU in our simulation, which results in significant reduction of net accretion rate with increasing radius. With this accretion rate profile, steady-state accretion is not possible, although the duration of our simulations is too short to show significant surface density evolution.
We emphasize that whether the laminar stress leads to accretion or decretion is also affected by the over-X.-N. Bai all density and temperature profiles, and our simulation Fid+ can be considered as one realization at the given disk model and magnetic flux distribution. While there can be many other possibilities, one general trend likely holds. In regions where the Hall effect is important and magnetic flux distribution is approximately uniform, α likely decreases with radius, and hence reduces or even reverses viscously-driven accretion rate. Given that the value of α is significant, the consequence of this effect on global disk evolution can be very profound, which requires further investigations in the future.

Wind Kinematics
Overall, beyond a few AU (where magnetic flux has not evolved significantly), the local mass outflow rate from our run Fid+ is comparable to that from run Fid0, reaching ∼ 10 −8 M yr −1 , as seen from the left panel of Figure 10. For this particular run, mass outflow rate there well exceeds accretion rate, due to the reduction of the latter from viscously-driven decretion.
The asymmetry in the inner zone also leads to an asymmetry in the outflow rate: the bottom side loses mass slower than the top side. In the outer zone near 20 AU where symmetry is retained, mass loss rate from the bottom side catches up and approaches that from the top side. The wind kinematics in the outer zone is similar to that discussed in the Hall-free case. Below we focus on the wind kinematics of the asymmetric inner zone.
In Figure 11, we choose R 0 = 9 AU, and trace poloidal field lines from the midplane both towards the upper and lower sides of the disk, and measure various diagnostics along the field lines similar to those done in Figure 5. The asymmetry is already evident by looking at the location of the Alfvén points, where the Alfvén radii are clearly larger at the lower side of the disk. This can also be tracked directly from Figure 7 over a broader range of radii. The smaller mass loss rate at the bottom side of the disk is consistent with larger Alfvén radii, due to Equation (30). We can also see from the middle panel of Figure 11 that toroidal velocity in the wind drops faster in the upper side than in the lower side, again consistent with the fact that wind from the upper side is more heavily loaded, as discussed in Bai et al. (2016).
To understand why wind in the lower side loses mass slower, we note that because poloidal field lines at 9 AU are significantly inclined, the wind bases at the top and bottom sides of the disk along this field line are located at different cylindrical radii. We see from the third panel that the angular velocity of the field line is approximately the same at the upper and lower sides and is close to the Keplerian frequency at radius R 0 . 11 Because wind in the lower side of the disk is launched from smaller radii, thus ω is smaller when normalized to the local Keplerian frequency at the wind base. Similarly as discussed before, this will lead to a larger Alfvén radius based on Bai et al. (2016).

FIDUCIAL SIMULATION IN THE ANTI-ALIGNED CASE
In the anti-aligned case, we find that the outcome of the simulation is insensitive to initial conditions, and hence the simulations are performed in the normal way as described in Section 2.5. The simulation is run for 24000Ω −1 0 ≈ 1775 yrs.

Overall Evolution and Magnetic Field Configuration
In Figure 12, we show the time evolution of magnetic field configuration from run Fid−. The initial stage of the evolution is very similar to those found in Bai & Stone (2017). Namely, instead of field amplification due to the HSI, horizontal components of the field are reduced towards zero in the midplane region. Overall, reflection symmetry across the midplane is preserved in this initial stage, and magnetic flux is transported outward due to the Hall drift. The rate of flux transport is rapid. By the time of 6000Ω −1 0 ≈ 443 yrs, the inner ∼ 3 AU is largely depleted of magnetic flux (two other field lines in the Figure are attached to the inner boundary), corresponding to a timescale of ∼ 100 local orbits. This rate is comparable to (by order-of-magnitude) while a factor of a few (∼ 3) slower than the rate of flux transport reported in Bai & Stone (2017) for the same given parameter (β 0 = 10 5 ). 12 We also note that at t = 6000Ω −1 0 ≈ 443 yrs, in between R ∼ 5−12 AU, the orientation of the poloidal field in the disk upper layers points radially inward and outward in an oscillatory manner (with time and radius). This leads to patches of toroidal fields of alternating polarities. This is closely related to the phenomenon observed in earlier local simulations, e.g., Figure 7 of Bai (2014) and Figure 9 of Bai (2015), where more detailed explanations were given. In brief, with anti-aligned vertical field, the disk becomes more susceptible to the MRI in localized region in the disk where the Hall Elsasser number is close to unity. In local simulations, such behavior may persist with time (e.g., in the aforementioned figures), making it ambiguous to interpret its global consequences. On the other hand, towards larger radii, toroidal field of one sign eventually overwhelms, terminating the oscillatory behavior, as seen in Figure  11 of Bai (2015). This is exactly what we observe in the global simulation Fid−, illustrated in the third panel of Figure 12. Once this sign of B φ is established, the inner region of the disk is quickly affected and settle to the new asymmetric configuration. This pattern also propagates outward, and by the end of our simulation, regions up to ∼ 20 AU are being affected. Interestingly, under the new asymmetric field configuration, outward transport of magnetic flux appears to be stalled, or at least significantly slowed down.
Towards the end of the simulation, the overall field configuration is well established within R ∼ 12 AU. While toroidal field is the dominant field component in the bulk disk, its strength is only modest, in fact the mean field strength smaller than the Hall-free case. Toroidal field then flips sharply at one (lower) side of the disk slightly below the FUV front, creating a strong current layer. This is where most of the accretion flow is concentrated (see the last panel), and the accretion flow is supersonic (see next subsection). In the upper side of the disk, we find that the poloidal field still show oscillatory behavior (mainly beyond ∼ 5 AU) for same reason discussed earlier. This leads to some slow motions in the midplane. Moreover, the deficit of magnetic flux within ∼ 3 AU means that the supersonic surface accretion flow suddenly stops at ∼ 3 AU. In reality, it plunges into the inner regions, causes strong disturbances, and also leads to rapid density variation in the outflows. This further affects the location of FUV fronts at larger radii, leading to variabilities near the entire strong current layer. We caution that because segregation of magnetic flux between the inner boundary and a few AU is the main cause of this behavior, it may be subject to the limitations in setting the inner boundary conditions discussed earlier.

Angular Momentum Transport and Flow Structure
In this section, we first focus on a representative radius R = 9 AU to analyze disk vertical structure, and address the mechanism of angular momentum transport in run Fid−.

Vertical Structure
In Figure 13, we show the vertical profiles of major diagnostic quantities at R = 9 AU. Without amplifying the horizontal field by the HSI, midplane field is much weaker than in run Fid+, with plasma β ∼ 10 3 . The disk maintains a small Maxwell stress −B R B φ in the bulk disk, corresponding to α ∼ 10 −4 .
The flow structure is again mainly determined by the vertical gradient of B φ base on Equation (25). From the last panel of Figure 13, we see that radial velocity in the midplane region is largely negligible. Essentially all the accretion flow is concentrated in the strong cur- rent layer. We note that the location of the strong current layer is very close to the FUV front where non-ideal MHD effects (dominated by AD) are greatly reduced, making the layer thin (but is well resolved by more than 10 cells). The sharpness of the strong current layer leads to very rapid radial inflows, where the velocity reaches ∼ 15% Keplerian speed. This is nearly twice the sound speed in that region. This layer is likely unstable to the Kelvin-Helmholtz instability, as found in surface current layers in some of Gressel et al. (2015)'s Hall-free simulations. In our case, despite the flow being supersonic, it is not straightforward to further discuss the stability of this layer due to the disturbances in this layer discussed in the end of the previous subsection.

Radial Profiles of Accretion Rates
Similar as in the aligned case, we show in Figure 14 the radial profiles of accretion and outflow rates. Overall, the accretion rate is around 10 −8 M yr −1 . Accretion is primarily wind-driven, as can be seen in the middle panel. Wiggles in the radial profiles of wind-driven accretion rate/wind stress are related to the oscillatory behavior discussed earlier. Similarly, the instantaneous α profile also exhibits wiggles. The corresponding viscouslydriven accretion rate shows accretion-decretion oscilla-tions, which time-averages to much smaller values than the absolute values shown in the Figure. The typical α values reach 2 − 3 × 10 −4 over a wide radial range, which is much smaller than in run Fid+, as expected.

Wind Kinematics
From Figure 14, we see there is a striking contrast between mass outflow rates in the upper and lower sides of the disk. Mass loss from the upper side is similar to those found in the grain-free case, whereas from the lower side, the mass loss rate is about two orders of magnitude smaller! This is related to the location of Alfvén surfaces, where we can identify from the third panel of Figure 12 that at the lower side of the disk, it is located at much larger distance than its counterpart at the upper side.
To understand this difference, we choose a characteristic radius R 0 = 8 AU, and trace poloidal field lines from the midplane towards both the upper and lower sides of the disk. In order to minimize the effect of disturbances to the field lines in the bottom side of the disk, we have averaged the data between t = 18000−24000Ω −1 0 (1331 − 1775 yrs) before tracing the field lines. In Figure 15, we show various diagnostics along the field lines similarly as in Figure 11. Note that even after time averaging, standard conservation laws are satisfied only ap- proximately, especially there are more deviations from the lower side of the wind (third panel).
We note that due to the large accretion velocity in the thin strong current layer, field lines strongly pinch radially inward there. Because this layer lies right below the FUV front, the wind base in the lower side of the disk has a radius that is smaller than R 0 . On the other hand, the wind base in the upper side is located at a radius larger than R 0 . The relatively contrast between the two wind base radii is large, amounting to a factor of ∼ 1.4 difference. Note that the angular velocity of the field lines ω at the upper and lower sides of the disk are approximately the same, and are comparable to Keplerian at radius R 0 (third panel). This means that there is a large difference in ω when being normalized to the wind base radii at the upper and lower sides of the field line. For similar reasons as discussed before, Alfvén radius must be larger at the lower side of the disk based on Bai et al. (2016), and hence much smaller mass loss rate.

PARAMETER STUDY
In this section, we briefly explore the role of net vertical field strength and the depth of FUV penetration, only focusing on different features exhibited in these simulations.

Disk Magnetization
With stronger magnetization (β 0 = 10 4 ), disk evolution proceeds much faster, and hence it suffices for shorter run time. Same as before, runs B40 and B4− are set up as described in Section 2.5, while run B4+ is restarted from run B40 as described at the beginning of Section 5. In Figure 16, we show the magnetic field configuration at the end of each simulation. The radial profiles of accretion and mass loss rates from these runs are shown in Figure 17.

The Hall-free Case
We start by discussing the Hall-free run B40. Overall, we see from the middle panel of Figure 16 that the field configuration is very similar to run Fid0. With stronger net vertical field, the disk becomes more stable and is fully symmetric about the midplane within the simulation domain. The toroidal field flips at the midplane beyond ∼ 4 AU, where the accretion flow is concentrated. Within that radius, toroidal field reduces to close to zero and the accretion flow splits into two branches above and below the midplane.
With stronger magnetization, wind-driven accretion rate increases by about an order of magnitude. Wind mass loss rate also increases, but by a smaller factor (∼ 3). Therefore, the wind becomes more lightly loaded, which is consistent with the fact that Alfvén surface is located higher than that in run Fid0. This trend is also consistent with expectations from semi-analytical theory .
The most interesting feature from this simulation is that because of higher wind mass loss rate, the wind becomes denser and FUV radiation is almost completely shielded from reaching the disk surface. This situation is distinct from earlier studies which attribute efficient wind launching largely to FUV ionization (which brings the gas in the disk surface layer to the ideal MHD regime) and associate the wind base with the FUV front.
To further understand the nature of this wind, we show in Figure 18 the vertical profiles of the Elsasser numbers and velocity components at a characteristic radius of R = 6 AU. In the wind zone, the dominant non-ideal MHD effect is AD, with its Elsasser number Am increases smoothly from ∼ 1 at about z = ±4H d to of order ∼ 10 or higher in the bulk wind zone. Without FUV, the wind column is instead mainly ionized by the X-rays (the direct absorption component). The smooth profile in Am also leads to a smooth flow structure, where we see that the v φ profile varies monotonically from midplane to surface, instead of having additional peaks as in run Fid0 (see Figure 3). We also find that in the wind zone, the conservation laws are not satisfied exactly, but still approximately.
We caution that the simulations in this work are   not designed to handle the situation with FUV being shielded, and the physical condition in the wind zone can be far from being realistic. In particular, without our prescribed heating beyond the FUV front, the wind zone is even colder than the disk zone due to rapid expansion. This helps reduce the recombination rate to achieve higher ionization level. In reality, heating from X-rays (e.g., Glassgold et al. 2004;Ercolano et al. 2009) and ambipolar diffusion (Garcia et al. 2001) can be substantial. Moreover, equilibrium chemistry may no longer apply, and even the prescription of X-ray ionization rate can become inaccurate in this region. More careful calculations are needed to better understand the wind properties in this regime (Panoglou et al. 2012;Wang & Goodman 2017). Overall, this simulation mainly demonstrates that assisted by X-ray ionization, MHD disk winds can still be launched when FUV radiation is shielded.

The Aligned Case
In the aligned case, we see from Figure 16 that again, the system relaxes to a state with an asymmetric inner zone within R ∼ 12 AU, and a more symmetric outer zone. 13 The overall mass wind-driven accretion rate and mass loss rate are comparable to the Hall-free run. With enhanced mass loss, FUV is again shielded and wind launching occurs in the AD dominated surface layer. Partly because of this, the vertical extent of the toroidal field patch (of single sign) in the asymmetric inner zone is smaller (instead of extending to FUV front). Correspondingly, the wind mass loss rates from the top and bottom sides of the disk are similar.
One general trend in the more strongly magnetized disks is that magnetic field amplification factor through the HSI is smaller. At 10 AU, we recall that midplane toroidal field in run Fid+ is amplified to ∼ 60 times the vertical field. In our run B4+ and at the same radius, we measure that midplane toroidal field is only amplified to ∼ 25 times the vertical field. Correspondingly, the laminar Maxwell stress increases more slowly than the increase in the wind torque. As a result, wind-driven accretion becomes more dominant, as can be seen in Figure  17. Comparing with the fiducial run Fid+ (Figure 10), radial transport of angular momentum clearly becomes less prominent in run B4+.

The Anti-aligned Case
With stronger net vertical field, the initial field evolution is similar to run Fid−, with horizontal field reduced to close to zero at the midplane maintaining reflection symmetry, and magnetic flux is rapidly transported outward. After about 1000 innermost orbits (∼ 460 years), the surface layer around a few AU becomes less stable, leading to symmetry breaking with toroidal field of a single sign dominates, with accretion proceeding at one side of the surface layer, again similar to run Fid−. However, this state is very short lived. Due to significant mass loss, the FUV front is again far from the disk surface, Fig. 19.-The first three panels are the same as Figure 16, but for simulations with higher FUV penetration depth (Σ FUV = 0.3g cm −2 ), run FUV+ (left), FUV0 (middle) and FUV− (right). Note that we have slightly time averaged the data spanning 60Ω −1 0 in making these plots. The rightmost panel shows the density structure at the last snapshot (plotted as ρR q D , no time average) of run FUV0. Green arrows indicate velocity vectors. and wind launching proceeds in the AD-dominated region. The system then settles to a state within ∼ 10 AU shown in the rightmost panel of Figure 16, where the strong current layer (and hence the accretion flow) wiggles around the midplane region. With these wiggles, we find that rapid loss of magnetic flux in this region is stalled. The region characterized by the wiggled strong current layer slowly moves outward over time, and regions beyond ∼ 10 AU are not yet affected by the end of the simulation (characterized by full symmetry across the midplane with rapid outward flux transport). While more detailed analysis is beyond the scope of this work, we note that in this state, the overall disk dynamics is much more symmetric and mass loss rate from the top and bottom sides of the disk are very similar.

FUV Penetration Depth
In Figure 19, we show the magnetic field configuration at the end of each simulation with deeper FUV penetration (Σ FUV = 0.3g cm −2 ). The radial profiles of accretion and mass loss rates from these runs are shown in Figure 20.
Overall, the evolution of magnetic field configuration share many similarities with the fiducial simulations. Without the Hall term, symmetry across the midplane is roughly preserved, with toroidal field being the dominant field component. With the Hall term in the aligned case, the HSI strongly amplifies the horizontal field, leading to an asymmetric inner zone and a symmetric outer zone with a smooth transition in between. In the antialigned case, horizontal field is first reduced towards zero in the midplane, followed by symmetry breaking. In the end, one sign of B φ dominates the bulk disk, with accretion flow concentrated at one side of the surface where this B φ flips.
With the Hall effect, the outflow from the top and bottom sides of the disk show very significant asymmetry for both aligned and anti-aligned cases. In particular, in the anti-aligned case, owing to the depletion of magnetic flux between the inner boundary to about R ∼ 4 AU (which may be unrealistic) and deep FUV penetration, the bottom side of the disk surface where B φ flips shows complex evolution (likely due to the MRI), and starts to tangle the poloidal magnetic flux by the end of our simulation (and the measured "outflow" rate becomes negative in some regions). Further investigation of this issue is desirable but beyond the scope of this work. The discussion of outflows in the next subsection, on the other hand, remains applicable to the top side of the disk for run FUV−.

Turbulent Outflow and Shielding
One important influence of larger Σ FUV is that it makes the outflow turbulent and drives much more significant mass loss. Three main features are worth discussing.
First of all, we find that in all three runs, the disk surface layer becomes unstable within radius r ∼ 2 AU, and launches episodic outflows. This is because for given vertical field strength, the well ionized disk surface layer near the FUV front becomes less magnetized, making it more susceptible to the MRI. As an example, the rightmost panel of Figure 19 shows the density and velocity structure from the last snapshot of the FUV0 run, which are clearly indicative of vigorous and turbulent outflow activities originating from the inner regions.
Second, beyond ∼ 3 AU, we find that the disk is relatively stable. This is mainly because that the unsteady outflow launched from smaller radii is so dense that it substantially shields the FUV radiation. In fact, the FUV front in run FUV0 at the distance of a few AU is located higher than that in run Fid0, despite that the former has much larger Σ FUV ! As a result, the disk surface layer is stable against the MRI.
Third, the bulk wind-driven accretion rate and wind mass loss rate in all three runs all exceed those in the fiducial case. Moreover, the ratio of mass loss to accre- tion rate is higher. This result first appears natural, and agrees with semi-analytical studies as the natural outcome of deeper FUV penetration ). This argument likely applies only in the innermost ∼ 2 AU (although the wind is episodic in this region), but does not apply to the outer region because the FUV front is located even higher than in the fiducial runs. The main reason for the enhanced accretion and mass loss rates is that the ram pressure from the heavily loaded episodic inner wind pushes the poloidal field lines and bends them further. This is similar to reducing the θ angle in the Bai et al. (2016) wind model. Although this parameter was only very briefly explored there, the trend is that more inclined field tends to develop stronger toroidal field, leading to both enhanced accretion and outflow rates, with the latter being more pronounced. This is consistent with what we observe in Figures 19 and 20.
Moreover, we notice that in many cases, the Alfvén radius is located within the FUV front. Namely, the development of the wind proceeds in the presence of strong non-ideal MHD effect. We have already discussed this phenomenon in the previous subsection. It again adds more complications to the dynamics of the system. On the other hand, we have also seen that semi-analytical theory is still useful that helps interpret the basic trend.
Finally, we caution that because part of the inner disk is MRI unstable, our 2D simulations are unable to fully characterize the flow properties, and hence the properties of the episodic winds launched from these regions. Full 3D investigations are necessary to resolve the gas dynamics more self-consistently.

Magnetic Flux Evolution
Given the very pronounced differences in simulations with different initial poloidal field strengths, it is clear that the overall disk evolution is largely controlled by the strength and radial distribution of poloidal magnetic flux. Following from Section 1.2, we now discuss the global evolution of magnetic flux based on our fiducial simulations. For each run, we compute the magnetic flux function Φ B,mid from Equation (31), and follow its evolution at three representative radii (different set of radii for different runs, since the dynamics and field configuration are different in each case). In Figure 21, we show the time evolution of dΦ B (t) ≡ Φ B,mid (r, t) − Φ B,mid (r, t 0 ), the amount of magnetic flux that has been transported through the chosen locations since time t 0 (time when magnetic flux is introduced, or when the Hall term is turned on in the case of run Fid+). The values are normalized by 2πB z0 R 2 , where B z0 is the initial vertical field strength, so that one can easily estimate the timescale of flux transport.
In the Hall-free run Fid0, we see that at R ∼ 2 AU, magnetic flux is consistently transported. From the figure, we estimate the timescale of flux transport to be ∼ 10 4 yr, translating to a speed of v B ∼ 3 × 10 −5 v K . This is more than an order of magnitude slower than the rate found in Bai & Stone (2017). As speculated there, the rate of flux transport is sensitive to the vertical diffusivity profile. On the other hand, the timescale of 10 4 year is still too short compared with disk lifetime. At R = 4 AU, flux is also transported outwards despite some irregularities, which is related to the secular evolution of toroidal field patches in Figure 1 that precludes accurate measurements of flux transport rate. At R = 10 AU, on the other hand, we do not find obvious signs of outward flux transport (despite small variations) within ∼ 1500 yrs, amounting to ∼ 50 local orbits of evolution. If the same rate of flux transport measured at ∼ 2 AU applies here, we would expect dΦ B /(2πB z0 R 2 ) ∼ 0.01 by ∼ 1200 yrs, yet this is not achieved. Therefore, the rate of flux transport at larger radii is even slower.
In the aligned case, run Fid+, we see that at all three radii, magnetic flux follow a pattern of being transported inward first (dΦ B > 0), followed by outward transport. This is consistent with the findings in Bai & Stone (2017), as a result of the HSI, followed by outward diffusion (at R = 18 AU, outward diffusion just starts to develop by the end of the simulation). There is some anomalous trend at R = 5 AU towards later stages, where we have discussed in Section 5.1 that might be an artifact due to inner boundary conditions. None of the three representative radii have achieved a quasi-steady state in magnetic flux evolution, which again makes it difficult to assess the overall rate of flux transport.
In the anti-aligned case, run Fid−, we see that at the beginning, magnetic flux is systematically transported outward at all radii, at a rate that is several times faster than the Hall-free case. This is overall consistent with the findings in Bai & Stone (2017). However, after the symmetry across the midplane is broken, a sudden change in the rate of transport is induced. At R = 5 AU, we find that flux is even transported slowly inward, whereas in between at R = 9 AU, flux evolution comes to a stall. We note that in this field configuration, the vertical profile of toroidal field is mostly flat, and hence substantially reduces the radial Hall-drift (proportional to ∂B φ /∂z). This is likely the main reason for the reduction and even reversal in flux transport. More detailed analysis about the direction and steady-state rate of flux transport is beyond the scope of this work, but is worth pursuing in the near future.
Overall, we have seen that while the initial stages of magnetic flux evolution in our simulations are similar to those found in Bai & Stone (2017), more complex behaviors are found as the system develops more complex field configurations and flow structures, resulting from more realistic prescriptions of disk microphysics.

Comparison with Other Works
As this work was in preparation, Béthune et al. (2017) (hereafter BLF17) conducted similar types of simulations of PPDs, and reported a variety of phenomena related to disk angular momentum transport, wind launching, etc. In methodology, there are two major differences between our simulations and theirs. First, the simulation domain in BLF17 extends to about ±60 • above/below the midplane, truncating at least part of the wind launched from the disk, as well as some of the magnetic flux in the simulation box. Our simulations do not suffer from this limitation. Second, our simulation domain covers a factor of 100 in radius as opposed to 10, which allows us to comfortably follow the launching and propagation of disk winds. There are several other differences at implementation level. To list a few, our simulations use more realistic tabulated magnetic diffusivities based on a complex chemical reaction network containing dust grains, as opposed to analytical diffusivity prescriptions mimicking the grain-free case in BLF17. We consider a flaring disk geometry, and use ray-tracing to estimate the radial and vertical disk column densities, whereas BLF17 considered flat disks (constant H d /R), and did not account for the radial column. The treatment of the transition from the disk zone to the atmosphere is very different. The transition in BLF17 is located at prescribed and constant latitudes, and is generally closer to midplane (by more than one scale height) than ours, leading to more unstable surface layer and more significant mass loss. Moreover, most of the BLF17 simulations are significantly more strongly magnetized, leading to accretion rates that are at least an order of magnitude higher than our fiducial runs.
Compared with a large variety of behaviors found in BLF17, we find a much more unified set of behaviors that are unique to the aligned and anti-aligned field polarities. In particular, BLF17 found "non-accreting" cases where B φ vanishes in the disk corona region, with meridional circulations but no wind-driven accretion. We do not observe this behavior: wind is always found leading to net accretion, on top of which there are meridional flows in the aligned cases as explained in Section 5.2. We speculate that the non-accreting solutions in BLF17 may be related to their over-constraining domain size. With much larger dynamical range and longer simulations, we have also identified and clarified regimes where the disk/wind structure become symmetric or asymmetric with respect to the midplane. Moreover, the diversity of behaviors in BLF17 is also likely related to the sensitivity to initial conditions in the aligned case, as we discuss in Appendix A. By mimicking conditions of disk formation, we have (at least partially) avoided this problem.
BLF17 also considered "cold" wind and "warm" wind, where the wind region is heated to very different temperatures. Most of the behaviors are consistent with the semi-analytical magneto-thermal wind framework of Bai et al. (2016). In particular, we note that keeping all other parameters fixed, while the "warm" wind simulations lead to much stronger acceleration due to thermal driving, the wind mass loading (or mass loss rate) remains similar to the "cold" wind. In other words, wind mass loss rate is largely determined by the conditions at the wind launching region, but not subsequent thermal accelerations. BLF17 described the wind from several of their warm simulations with different properties as "magneto-thermal". In Appendix C, we aim to systematically clarify the nomenclature on PPD winds and the corresponding phenomenology.
Finally, we comment that we do not observe the magnetic flux concentrations and zonal flows in our 2D simulations. The phenomenon identified and explained in BLF17 is found in highly strongly magnetized (and hence accretion rate is several orders of magnitude higher) disks with net vertical flux β 0 = 10 2 , a regime not explored in our simulations. We also comment that the mechanism is also likely related to the sharpening of flux concentration observed in Bai (2015).

Implications for Planet Formation and Disk
Evolution The most important implication of this work on planet formation is from the complex flow structures. In particular, in simulations with aligned field polarity, the presence of both accretion and decretion flows at different heights in the bulk disk is completely unexpected based on conventional models (especially viscous evolution models) of accretion disks. It poses very interesting questions on how it would affect the transport of solids, and subsequent stages of planet formation. The systematic outward motion in a substantial fraction of the bulk disk may be an important source for large-scale mixing. Evidence for such radial mixing in the solar system has been mounting, especially based on the findings of crystalline silicates in comets (e.g., Bockelée-Morvan et al. 2002 and references therein), and more directly from samples collected from comet 81P/Wild 2 (Brownlee et al. 2006;Nakamura et al. 2008). Similar evidence has been found in nearby PPDs (van Boekel et al. 2004;Watson et al. 2009). Viscous diffusion with large-scale radial flows have been commonly involved to explain such largescale mixing (e.g., Keller & Gail 2004;Ciesla 2009;Hughes & Armitage 2010), as well as variations in various isotopic ratios in the solar system, such as the D/H ratio in water (Jacquet & Robert 2013;Yang et al. 2013;Albertsson et al. 2014). Our simulation results offer a first-principle demonstration of the large-scale flow structure that differ substantially from the conventional picture. With typical radial flow velocity of the order ∼ 1% sound speed, and given that the sub-Keplerian velocity in the inner regions of PPDs is around ∼ 5% of sound speed, this means that decretion flow in the bulk disk can overcome radial drift for particle Stokes numbers St 0.1. This flow thus has the potential to transport mm-sized particles to ∼ 30 − 40 AU scale (if our results can be generalized to outer radii) in a standard MMSN disk. More detailed calculations are necessary to further demonstrate its feasibility.
Global evolution of PPDs is determined by the transport of angular momentum, both radially and vertically, as well as mass loss. Our simulations with aligned field polarity also show dramatic radial variations in accretion rate due to significant contribution from the laminar Maxwell stress as a result of the HSI (whose rate and even direction depends on the radial gradient of the stress). This fact, together with significant mass loss rate, implies that global disk evolution is highly com-plex, and it is unclear whether a steady state can ever be achieved. Uncertainties in magnetic flux transport and evolution discussed earlier add further complications.
The discussions above mainly focused on the aligned case. PPD gas dynamics in the anti-aligned case shows completely different behaviors. It is thus very likely that planet formation takes very different pathways in these two cases, though the details remain to be filled in upon better understandings of long-term disk evolution (but see a toy model by Simon 2016).
Finally, we speculate that the solar nebular was once threaded by poloidal field with aligned polarity, for two reasons. First, large-scale outward radial flows are only present in simulations in the aligned case. While outward transport is also possible in the anti-aligned case through diffusion, the weak level of turbulence is unlikely to lead to efficient large-scale mixing. Second, recent paleomagnetic measurement of the Semarkona meteorite revealed a strong magnetic field of ∼ 0.5 Gauss (Fu et al. 2014), 14 presumably corresponding to the asteroid belt region in the midplane. We note that that in the aligned case, radial transport of angular momentum is comparable to wind contributions. For typical accretion rate of 10 −8 M yr −1 , ∼ 0.5G is exactly the expected field strength at ∼ 2 − 3 AU scale from radial transport of angular momentum (see Wardle 2007;Bai & Goodman 2009). Comparing Figures 7 and 12, we see that the midplane field in the anti-aligned case is typically a factor of 4 − 5 times lower, inconsistent with the paleomagnetic measurement.

Connection to Disk Observables
One important prediction from this work is that disk winds from the inner region of PPDs are likely asymmetric between the two sides. While wind signatures seem to be ubiquitous among T Tauri disks (Hartigan et al. 1995;Natta et al. 2014;Simon et al. 2016) based on optical-infrared forbidden line (blue-shifted) observations, these observations typically see winds only from one side, and the observations themselves already bare large uncertainties in constraining wind kinematics. Jets from young-stellar-objects often show asymmetric signatures between the jet and counter-jet (e.g., Hirth et al. 1994;Woitas et al. 2002;Hartigan & Hillenbrand 2009;Liu et al. 2012), although it is less clear whether such asymmetry extends to the lower-velocity wind components. Recently, ALMA has revealed large-scale molecular outflows in several sources (e.g., Klaassen et al. 2013;Salyk et al. 2014;Bjerkeli et al. 2016), which all show complex spatial and velocity structures, some of which are only one-sided. Environmental effects (e.g., envelope, foreground, tidal interaction with binary) may be important contributing factors, but it is unclear whether some are caused by intrinsic asymmetry during the wind launching processes. Very encouragingly, Klaassen et al. (2016) derived the kinematics of disk winds from the HL Tau disk (despite the systematics), and found dramatically differences (by ∼an order of magnitude) between the redshifted and blueshifted sides.
The complex flow structures in the bulk disk found in our aligned simulations reach systematic radial velocities of up to a few percent of the sound speed. We note that careful modeling of ALMA data has already enabled level of turbulence to be constrained at a comparable precision (Flaherty et al. submitted). We thus expect that such systematic flow structures that depart from Keplerian rotation to be potentially detectable. Moreover, some specific accreting layers in the disk surface, in both aligned and anti-aligned cases, have accretion velocities near or exceed the sound speed. Note that these layers are very thin and contain only a very small fraction of disk mass. Detecting such flow structure would provide smoking-gun evidence of our simulation predictions, but it is also very challenging, since it would require specific tracers whose optical depth τ ∼ 1 surface is right in the vicinity of the thin accreting layer.

SUMMARY
In this work, we conducted the most comprehensive/realistic global simulations of the inner regions of PPDs to date, that have incorporated all non-ideal MHD effects coupled to steady-state chemistry with dust grains, as well as proper ray tracing schemes to calculate disk ionization and control thermodynamics. All simulations include net poloidal magnetic flux, which is an essential ingredient to launch MHD disk winds and drive disk accretion. We have largely focused on a set of fiducial simulations, which give accretion rates on the order of 10 −8 M yr −1 , but also briefly explored the role of poloidal field strength and FUV penetration depth. Our main findings from the fiducial simulations are as follows.
• The bulk disk is largely laminar, launching an MHD disk wind that drives disk accretion. The wind is magneto-thermal in nature, launched by magnetic pressure gradient, with very strong mass loss rate that is comparable or larger than winddriven accretion rate.
• In the aligned case, the Hall shear instability (HSI) strongly amplifies horizontal field, making the outcome dependent on initial field configuration. Mimicking realistic initial conditions, we find that the disk is divided into an asymmetric inner zone (within ∼ 10 AU) and a more symmetric outer zone that are smoothly connected. Accretion and decretion flows at the ∼ 1% of sound speed coexist in the bulk disk, determined by the vertical gradient of toroidal field.
• In the aligned case, both MHD wind and the laminar Maxwell stress contribute at comparable level to disk accretion. However, since the latter contribution depends on its radial gradient, making local accretion rates sensitive to radial disk structure, and the accretion process is likely non-steady.
• In the anti-aligned case, the disk may achieve a symmetric state that rapidly loses magnetic flux, or (more likely) an asymmetric state that retain magnetic flux. In the latter case, accretion is predominantly wind-driven, with the bulk accretion flow located at one side of the disk surface at transsonic to supersonic speed. The resulting disk wind is highly asymmetric, with most mass loss at the opposite side of the accretion flow.
In addition, increasing poloidal field strength enhances mass accretion rate more than enhancing mass loss rate. By contrast, increasing FUV penetration enhances mass loss rate more than accretion rate. Moreover, the mass loss rates found in these additional simulations are sufficiently high to substantially shield the incoming FUV radiation. We find wind launching can still operate in the X-ray ionization dominated disk surface layer when FUV is largely shielded, with gas marginally coupled to the magnetic field through ambipolar diffusion. We raised several outstanding issues in PPD gas dynamics in Section 1.2, and our global simulations have largely clarified the issues on wind kinematics and symmetry. Our simulations have also found complex behaviors in magnetic flux evolution, but in general the disks appear to be able to retain magnetic flux, or lose flux much more slowly than found in idealized simulations (Bai & Stone 2017).
Our simulation results also have major implications on planet formation, especially that the complex flow structure may transform our understandings on how solids are transported in disks. It also calls for a reassessment of other stages of planet formation and migration. Based on the results, we further speculate that the solar nebula was originally threaded by poloidal fields aligned with disk rotation.

Limitations and Future Directions
Given the richness of the results from our fiducial simulations, we have only very briefly explored the parameter space. We expect the simulation results to be representative, and the physics we have explained to be widely applicable. Further parameter exploration may lead to additional variations and complications. In particular, the gas dynamics of the bulk disk can be affected by grain abundance, as well as ionization rates, especially the X-ray properties of the protostar.
Several of our simulations (especially those with antialigned polarity and deeper FUV penetration) show signs of turbulence. Moreover, the stability of the strong current layer from the HSI in the aligned case also requires further investigation. Future 3D simulations are necessary to properly characterize their properties.
In addition, we have only treated thermodynamics in the bulk disk and the wind zone very roughly. Fully self-consistent calculations would require coupling radiative transfer and chemistry (including photo-chemistry) with dynamics (as in some photo-evaporation simulations, e.g., Owen et al. 2010;Wang & Goodman 2017). These treatments are essential to better determine wind kinematics and compare with observations, and call for future investigations.
Finally, we have only explored the inner region of the disk (∼ 0.6 − 20 AU). Future explorations should also focus on other radial ranges, including the innermost re-X.-N. Bai gion (e.g., Flock et al. 2017) where most exoplanets are found, and the outer disk regions which are more accessible with spatially-resolved observations. Moreover, simulations would also benefit from using further larger domain size to potentially capture the fast magnetosonic points, and to cover broader dynamical ranges. I thank the referee for a prompt and detailed report, and acknowledge support from Institute for Theory and Computation, Harvard-Smithsonian Center for Astrophysics. Computations for this work are performed on the Hydra cluster managed by the Smithsonian Institution, and on Stampede at the Texas Advanced Computing Center through XSEDE grant TG-AST140001.

A. DEPENDENCE ON INITIAL CONDITION IN THE ALIGNED CASE
In this appendix, we address the dependence on initial conditions in simulations with all three non-ideal MHD effects included and initial vertical field aligned with disk rotation. In Figure 22, we show the time evolution of magnetic field configuration from a modified version of Fid+, where all three non-ideal MHD terms are turned on from the beginning (instead of turning on the Hall term from inside out). In other words, the setup is the same as run Fid− except that field polarity is flipped.
We immediately notice that during the evolution, there are patches in the disk that possess strong toroidal field with opposite and alternating signs. The presence of these patches is a result of the Hall-shear instability (HSI), as discussed in Kunz (2008); Lesur et al. (2014) and Bai (2014) via local analysis and simulations. The global manifestation of the HSI was studied in Bai & Stone (2017), and is briefly reviewed in Section 5.1.
The nature of the HSI dictates that the formation of such discrete patches depends on initial condition, and can be stochastic. For instance, a random perturbation in vertical field can create regions of radial field with opposite signs, with each region growing their own HSI, eventually creating a pair of relatively strongly magnetized patches. Indeed, we have also run simulations with different field geometries (controlled by the m parameter in Equation (20)), and found qualitatively similar outcomes except that these patches are distributed differently.
The nature of the HSI also dictates that these strongly magnetized patches (but still not sufficiently strong to become magnetically dominated), once formed, do not annihilate with neighboring patches with opposite sign of toroidal field. Instead, each patch attempts to expand its "territory" through the growth of HSI. At the boundaries, the balance between growth and resistivie/ambipolar dissipation allows these patches to survive and evolve slowly. Because the flow structure in the disk strongly depends on the vertical gradient of toroidal field described by Equation (25), we see in the rightmost panel of Figure 22 that strong radial flows are induced at patch boundaries.
The sensitive dependence on the initial conditions may suggest that in the aligned case, the overall magnetic field structure, and hence the flow structure in the disk is somewhat unpredictable. While this is entirely possible and worth further investigations, it also suggests that we need to consider setting up the initial conditions in a more realistic way that mimics the initial stage of disk formation.
Formation of PPDs follows from protostellar core collapse, and it is well known that non-ideal MHD effects play an important role throughout the processes of core collapse and disk formation (see Li et al. 2014 for a review, and more recent works by Tomida et al. 2015;Tsukamoto et al. 2015;Wurster et al. 2016). We note that these processes are accompanied by rapid increase of gas density, and the development of more rapid rotation that winds up poloidal field into toroidal field. According to the relative ordering among the three non-ideal MHD effects, AD is the dominant effect at the beginning, the Hall term gradually picks up as density increases. The key notion here is that before the Hall effect becomes important, the system has already evolved under AD with well-defined sense of (differential) rotation to produce toroidal field with ordered vertical structures. As density builds up further, the Hall effect enters and modifies the field structure through the HSI.
The discussion above motivates us to adopt the procedure described in Section 5 to set up the initial conditions. The key is to allow the system to develop some initial toroidal field structures before introducing the Hall effect, so that the HSI can grow on top of these pre-existing structures. This procedure, when applied to run Fid+, is very successful in minimizing the number of discrete magnetized patches to result in sufficiently simple field structure, which also appears more physically reasonable.
However, the same procedure appears less successful for some other runs in our parameter study (Section 7). This again reflects the sensitive dependence of the outcome on initial conditions, but it is also true that our new procedure of setting up the initial condition is still far from representing realistic conditions of disk formation. In view of these results, we are still argue for the generality of the field structure obtained in run Fid+, although future works are needed to reach firm conclusions.

B. FURTHER DISCUSSIONS ON THE SYMMETRY ISSUE IN THE ALIGNED CASE
Our simulation results resolve a number of puzzles found in previous local simulations. First, previous local shearingbox simulations of Lesur et al. (2014) and Bai (2014) always found that in the inner disk, toroidal field of a single sign overwhelms the entire box after the development of the HSI even the vertical box extends into the FUV layer. The failure for B φ (and hence B R ) to flip across the disk means that while the system launches an outflow, there is no net transport of angular momentum. Our results show that toroidal field does flip, and a fully global setup is essential which avoids the influence of artificial boundary conditions in local simulations.
Second, local simulations of Bai (2014Bai ( , 2015 also found that towards outer radii ( 5 AU), flipping the toroidal field within the simulation box is possible, but it was unclear how this field configuration connects to the inner radii. Our simulations demonstrate that they join smoothly.
In brief, the inner zone with strong asymmetry is related to strong Hall effect. Development of such asymmetry is unavoidable in these simulations because within a few AU, resistivity is so strong that a flipped field configuration is not sustainable. This is the reason why simulations of Bai & Stone (2017) also dominated by the Hall effect yet does not produce the asymmetric field configuration.

C. NOMENCLATURE OF PPD DISK WINDS
In this appendix, we aim to clarify the nomenclature of disk winds from PPDs, in light of recent development of magnetized wind theory and global disk simulations. Broadly speaking, disk winds can be thermally-driven (by thermal pressure), magnetically-driven (by the Lorentz force), or driven by radiation pressure (e.g., Proga 2007). For PPDs in T-Tauri phase, radiation pressure is largely negligible due to their low luminosities (e.g., Cabrit 2007). We focus on thermal and magnetic effects.
Thermally-driven PPD wind is generally referred to as photoevaporation. It results from external heating by highenergy photons (UV to X-ray), either from the central protostar, or nearby massive stars. Extensive literature has focused on this type of disk winds (see Alexander et al. 2014 for a review), and the outcome is very sensitive to details of the heating and cooling processes.
Magnetically-driven wind has two flavors.
• Magnetocentrifugal wind: strong poloidal field lines anchored to the disk enforces the outflow to corotate with the wind foot-point (wind base), leading to centrifugal acceleration as viewed from the corotating frame when field geometry is favorable. It is directly related to the Blandford & Payne (1982) scenario.
• Wind driven by magnetic pressure gradient: poloidal field lines are too weak to enforce corotation, and get winded up to build up toroidal field. The wind is launched by vertical gradient of magnetic pressure from the toroidal field. It is sometimes referred to as magnetic tower flows (Lynden-Bell 1996, 2003. Viewed in the observer's frame, the Lorentz force is the driving force in both cases, again from the pressure gradient and tension force from the toroidal field (Spruit 1996). Nevertheless, it is physically intuitive to distinguish the two regimes, as is widely adopted in other contexts such as star formation (e.g., Seifried et al. 2012;Tomida et al. 2013). The centrifugally-driven wind typically corresponds to large Alfvén radius (long lever arm), and the wind is lightly loaded. The opposite regime corresponds to small Alfvven radius (short lever arm), with a heavily loaded wind. Fixing other conditions, one generally smoothly transitions from magnetic pressure gradient driven wind to centrifugally driven wind as poloidal field strength increases ). Conventionally, magnetically-driven wind models/simulations usually assume a cold gas flow and hence thermal pressure plays a negligible role throughout the process (with a few exceptions, such as Casse & Ferreira 2000). In principle, both magnetic and thermal effects can contribute to the launching and acceleration the wind flow.
We call a disk wind "magneto-thermal" when the wind properties are affected both the field strength and thermodynamics. This term can be considered to be broadly defined, encompassing the aforementioned wind driving mechanisms as long as both magnetic and thermal effects play a role. One can imagine that by applying poloidal fields to a pure thermal disk wind and increasing the field strength, the wind will transition to being driven by magnetic pressure gradient, and eventually become centrifugally-driven. In general, PPD winds are magneto-thermal, because poloidal magnetic fields are essential to drive disk accretion, and strong external heating is inevitable that can drive a thermal wind on its own. We have shown in Bai et al. (2016) as well as in this paper that magnetic pressure gradient is the main wind launching mechanism.
Many examples of magneto-thermal winds are explored in idealized models of Bai et al. (2016), where the wind is assumed to be launched in the ideal MHD regime from the wind base at the disk surface that is well separated from the poorly ionized disk main body, with a barotropic equation of state. The reality can always be more complex, as some of our simulations and the ones in BLF17 illustrate. Here we list two scenarios where wind properties are strongly modified by thermal effects, which can be considered as key characteristics of magneto-thermal disk winds: • Wind launching takes place where magnetic pressure is not much larger than thermal pressure, as studied in Bai et al. (2016). The resulting wind is typically very heavily loaded, with very small lever arm. Sometimes it even violates the requirement that the lever arm λ > 3/2 for a cold MHD wind (Casse & Ferreira 2000), as we have observed in several cases in our fiducial simulations, as well as in BLF17.
• The wind is launched magnetically, but is subsequently accelerated by thermal pressure from external heating. This case includes the hot wind simulations in BLF17. We studied a class of models of this type in Bai et al. (2016), finding that despite the subsequent wind acceleration from external heating, wind lever arm and hence wind mass loss rate is largely unaffected. This is also confirmed in the BLF17 simulations.
Note that the above two scenarios do not necessarily exclude each other. In addition, we have also found that when FUV is shielded by the wind itself (as we find in the case of deeper FUV penetration and stronger magnetization), the wind can be launched, and achieve super-Alfvénic velocities within the non-ideal MHD layer. Similar results are found in BLF17. These can be considered as extensions of the first scenario, though gaining more quantitative understandings is much less straightforward. Finally, the terms "MHD disk wind", or "magnetized disk wind", generally refer to disk winds that are launched magnetically without specifically referring to thermal effects. We consider these terms to be more broad and inclusive, and can be applied to magnetized PPD winds in general.