Global Three-Dimensional Simulations of Outer Protoplanetary Disks with Ambipolar Diffusion

The structure and evolution of protoplanetary disks (PPDs) are largely governed by disk angular momentum transport, mediated by magnetic fields. In the most observable outer disk, PPD gas dynamics is primarily controlled by ambipolar diffusion as the dominant non-ideal magnetohydrodynamic (MHD) effect. In this work, we study the gas dynamics in outer PPDs by conducting a set of global 3D non-ideal MHD simulations with ambipolar diffusion and net poloidal magnetic flux, using the Athena++ MHD code, with resolution comparable to local simulations. Our simulations demonstrate the co-existence of magnetized disk wind and turbulence driven by the magneto-rotational instability (MRI). While MHD wind dominates disk angular momentum transport, the MRI turbulence also contributes significantly. We observe that magnetic flux spontaneously concentrate into axisymmetric flux sheets, leading to radial variations in turbulence levels, stresses, and accretion rates. Annular substructures arise as a natural consequence of magnetic flux concentration. The flux concentration phenomena show diverse properties with different levels of disk magnetization and ambipolar diffusion. The disk generally loses magnetic flux over time, though flux sheets could prevent the leak of magnetic flux in some cases. Our results demonstrate the ubiquity of disk annular substructures in weakly MRI turbulent outer PPDs, and imply a stochastic nature of disk evolution.


INTRODUCTION
Protoplanetary disks (PPDs) are the birth place of planets, where dust and ice particles coagulate into planetesimals, seeding further growth into terrestrial planets and cores of gas giants. PPDs accrete onto the central protostar at typical rate of ∼ 10 −8±1 M⊙ yr −1 (e.g., Hartmann et al. 1998;Herczeg & Hillenbrand 2008), with typical lifetime of several million years (Haisch et al. 2001;Mamajek 2009). The accretion and dispersal of PPDs are governed by disk gas dynamics, which is also the key to understanding many physical processes in planet formation. In particular, the mechanisms of angular momentum transport shape the structure and govern the global evolution of the disk, and also determine the level and nature of disk turbulence, thereby affecting virtually all stages of planet formation (Armitage 2010).
In the recent years, the Atacama Large Millimiter/submillimeter Array (ALMA) has been revolutionary that uncovers a diverse range of disk substructures (ALMA Partnership et al. 2015;Andrews 2020), particularly annular substructures (e.g Isella et al. 2016;Andrews et al. 2016;van Boekel et al. 2017; Avenhaus et al. ⋆ E-mail: cc795@cam.ac.uk † E-mail: xbai@tsinghua.edu.cn 2018; Long et al. 2018;Huang et al. 2018a), spirals (e.g Benisty et al. 2015;Pérez et al. 2016;Huang et al. 2018b), and crescent-shaped azimuthal asymmetries (e.g van der Marel et al. 2013van der Marel et al. , 2020. These features have triggered extensive research on their origins. They could be interpreted as signposts of unseen planets embedded in the disk, reflecting planet formation processes caught in action, while they may simply reflect local concentration of dust that could potentially become precursors to planet formation.
Below, we give a brief overview on the above two aspects of recent development before introducing the objectives of this study.

Angular Momentum Transport
Magnetic fields are known to play a leading role in driving angular momentum transport in PPDs (e.g., see recent reviews by Turner et al. 2014;Weiss et al. 2021). The driving mechanisms can be broadly divided into two categories -radial transport of angular momentum mediated by the turbulence through the magnetorotational instability (MRI; Chandrasekhar 1961;Balbus & Hawley 1991) and vertical transport of angular momentum through magnetized disk winds (Blandford & Payne 1982).
Conventionally, the MRI turbulence is considered as the primary mechanism to transport angular momentum, which requires the perfect coupling between disk gas and magnetic fields. However, the extremely weakly ionized gas in PPDs undermines this coupling (Gammie 1996;Wardle 2007;Bai 2011a), which is measured by three non-ideal magnetohydrodynamic (MHD) effects -Ohmic resistivity, Hall effect, and ambipolar diffusion (AD). Among them, Ohmic resistivity and AD are dissipative and effectively stabilize the MRI such that in inner disks the MRI is suppressed by Ohmic resistivity at the midplane (e.g., Fleming & Stone 2003;Ilgner & Nelson 2008;Oishi & Mac Low 2009) and by AD at the disk surface Bai 2013;Gressel et al. 2015). In outer disks, AD is the dominant nonideal MHD effect. The MRI is found to be damped (but not suppressed) by AD at the midplane, but is expected to become more vigorous at the disk surface due to substantially boosted ionization levels by stellar far-ultraviolet (FUV) radiation (Perez-Becker & Chiang 2011;Simon et al. 2013a,b;Bai 2015).
With the MRI turbulence being largely suppressed or damped, several purely hydrodynamic instabilities have been identified and systematically explored (e.g. Klahr & Bodenheimer 2003;Nelson et al. 2013;Marcus et al. 2013;Klahr & Hubbard 2014;Lin & Youdin 2015;Cui & Lin 2021). The onset and outcome of these instabilities sensitively depend on disk thermodynamics, despite that the resulting turbulence level remains relatively low, insufficient to account for the bulk accretion rate in PPDs (see Lyra & Umurhan 2019 for a recent review).
In the meantime, it is revealed that magnetized disk winds play a major role in driving disk accretion. This has been robustly demonstrated in the inner disk regions in both local Bai 2013Bai , 2014Lesur et al. 2014) and global simulations (Gressel et al. 2015;Bai 2017;Wang et al. 2019;Gressel et al. 2020), where all three non-ideal MHD effects are relevant and dominate in different vertical heights. As a prerequisite, disk winds requires the disk be threaded by net vertical magnetic flux, which is likely inherited from the star formation process (see, e.g., Zhao et al. 2020 for a review). The efficiency of angular momentum transport is then directly linked to the amount of magnetic flux threading the disk.
In outer disks, however, the situation remains ambiguous. Although high resolution local simulations can appropriately resolve the MRI (Simon et al. 2013a,b;Bai 2015), the role played by disk winds is uncertain for two reasons. First, local framework ignores disk curvature, hence cannot distinguish between inward and outward radial directions, which is particularly problematic when the disk is turbulent . Second, the disk surface is artificially truncated by the vertical boundaries that is not able to accommodate the wind propagation (Fromang et al. 2013, but see Riols et al. 2016). On the other hand, existing global simulations are either in 2D thus cannot capture the MRI (e.g., Bai & Stone 2017;Béthune et al. 2017;Suriano et al. 2018;Riols et al. 2020), or in 3D with resolution too coarse to properly resolve the MRI (e.g., Béthune et al. 2017;Suriano et al. 2019;Hu et al. 2019). Although these studies show efficient winddriven angular momentum transport, it is unclear that to what extent their results will be undermined due to the lack of third dimension and/or resolution. Therefore, a fundamental question left unanswered is whether magnetized disk wind can still be launched in outer PPDs, and the relative importance angular momentum by the MRI and the wind.

Annular Substructures
Annular structures are one of the most common type of disk substructures (see review by Andrews 2020). They are most commonly considered as a consequence of planet-disk interactions (e.g. Goldreich & Tremaine 1979;Lin & Papaloizou 1986, 1993, with significant recent development to interpret disk observations (e.g. Bryden et al. 1999;Nelson et al. 2000;Zhu et al. 2011;Kley & Nelson 2012;). This interpretation is recently corroborated by additional signatures in gas kinematics, suggesting the presence of Jupitermass planet responsible of carving deep gaps in a few systems (Teague et al. 2018;Pinte et al. 2019).
While the planet scenario is especially compelling, it does not necessarily apply to all annular substructures observed. In particular, rings and gaps are also discovered during the embedded Class 0/I phase (Sheehan & Eisner 2018;Nakatani et al. 2020;Segura-Cox et al. 2020;Sheehan et al. 2020) at early stages of disk formation, requiring either rapid formation of massive planets, or alternative, planet-free interpretations.
Indeed, there is extensive recent literature on alternative origins of disk rings. To name but a few, one set of proposed mechanisms is related to the condensation fronts of various volatile species, which mark a transition in dust properties (e.g., Zhang et al. 2015;Okuzumi et al. 2016;Hu et al. 2019;Owen 2020). Another set of scenarios is associated with gas-dust coupling, including secular gravitational instability (Shariff & Cuzzi 2011;Youdin 2011;Takahashi & Inutsuka 2014;Tominaga et al. 2018Tominaga et al. , 2019, dust-driven viscous ring-instability (Wünsch et al. 2005;Dullemond & Penzlin 2018), and self-induced dust traps (Drazkowska et al. 2016;Gonzalez et al. 2017).
Forming annular substructures generally requires a redistribution of disk mass, and hence disk angular momentum. Since magnetism dominates the angular momentum transport in PPDs, any non-smoothness in angular momentum transport associated with magnetic fields must directly lead to the formation of annular substructures. This leads to another broad category of mechanisms to generate rings and gaps. These include MRI-induced zonal flows (e.g., Johansen et al. 2009;Dittrich et al. 2013;Simon & Armitage 2014), gas density accumulation at the dead-zone boundaries (Varnière & Tagger 2006;Flock et al. 2015), and spontaneous magnetic flux concentrations (Bai & Stone 2014;Béthune et al. 2017;Suriano et al. 2018Suriano et al. , 2019Riols & Lesur 2019;Riols et al. 2020) We note that zonal flows and magnetic flux concentration are intimately related (Bai & Stone 2014), which are the main focus of this work. It was found in local and global disk simulations that through the development of the MRI and launching of magnetized disk winds, the global distribution of poloidal magnetic flux threading the disk tends to become non-uniform and spontaneously concentrate into quasiaxisymmetric flux sheet. Regions with stronger magnetic flux transports angular momentum more efficiently and become depleted (i.e., gaps), whereas mass tends to accumulate in regions in between flux sheets (i.e., rings). However, the aforementioned works were conducted either as local simulations, which again suffer from numerical artifacts, or global simulations that significantly under-resolve the MRI. It is yet unclear whether magnetic flux concentration is a robust phenomenon, how magnetic flux concentrates at first place, hence the characteristics of the resulting features.

This Work
In this paper, we conduct global 3D non-ideal MHD simulations of the outer PPDs. Our simulations include net poloidal magnetic flux, with domain size opening up to the polar region to fully accommodate wind launching, while in the meantime feature a numerical resolution comparable to those achieved in local shearing box simulations. As an initial study, we avoid complications from evaluating the diffusion coefficients self-consistently from ionization chemistry and thermodynamics. Instead, our simulations are designed to be scalefree to systematically scan the parameter space by prescribing the strength of AD and level of magnetization. In doing so, we aim to clarify the dominant mechanism governing angular momentum transport in outer PPDs, and in the meantime, we will verify how disks naturally form annular substructures from magnetic flux concentration.
The paper is organized as follows. In §2, we describe numerical methods and simulation setup. In §3, we list diagnostic quantities to facilitate the analysis of simulation results. The results of fiducial models including gas dynamics and magnetic flux concentration phenomenon are detailed in §4. We conduct parameter study on the strength of ambipolar diffusion and magnetic field strengths in §5. Implications of our results and mechanism of magnetic flux concentration are further discussed in §6. Finally, we conclude in §7.

Dynamical Equations
We use the grid-based high-order Godunov MHD code Athena++ (Stone et al. 2020) to carry out the simulations in this work. Athena++ is a rewrite from the predecessor Athena code (Gardiner & Stone 2005, 2008Stone et al. 2008) in C++ language with a comprehensive set of physics and numerical capabilities, and significantly improved performance thanks to its modular design and the use of dynamical execution model. We solve the equations of non-ideal MHD (as implemented in Bai & Stone 2017) using the van Leer time integrator, the HLLD Riemann solver, and the piecewise linear reconstruction. Super time-stepping is employed (based on Alexiades et al. 1996, described in Simon et al. 2013b) to speed up calculations for diffusive non-ideal MHD physics, as was extensively used in our previous works. The dynamical equations in the conservative form read Here, v, ρ, and P are gas velocity, density, and pressure, and B is the magnetic field vector. The moduli of velocity and magnetic vectors are v = |v| and B = |B|. The stress tensor M is defined as where P * = P * I with I being the identity tensor, P * = P + PB is the total pressure, and PB = B 2 /8π denotes the magnetic pressure. The total energy density is E = ǫ+ρv 2 /2+ B 2 /8π, where ǫ is the internal energy density and is related to the gas pressure by an ideal gas equation of state P = (γ − 1)ǫ. We adopt adiabatic index γ = 7/5 for the molecular gas in the outer disks. The cooling term Λc in the last equality will be elaborated in §2.2. The gravitational potential of the central star is implemented as a source term given by Φ = −GM/r, with stellar mass M . The components of non-ideal MHD terms manifest in the induction Equation (3). The electric field in the rest fluid frame is where Ohmic and ambipolar diffusivities are denoted by ηO and ηA. The current density is given by J = c∇ × B/4π, where we express the component of J that is perpendicular to the magnetic field by J ⊥ = −(J ×b) ×b, and the unit vector of magnetic field is denoted byb = B/B. The Poynting flux associates to the non-ideal MHD is written as The above governing equations are written in Gaussian units, whereas in code units the 4π factors are absorbed to the definition of magnetic pressure such that the magnetic permeability µ0 = 1. Our simulations are conducted in spherical-polar coordinates (r, θ, φ). In code units, we choose GM = R0 = 1, with R0 a reference radius at the inner radial boundary. The domain spans over r ∈ [1, 100] in radius, in polar angles over θ ∈ [0, π], and in azimuthal angles over φ ∈ [0, π/4]. The radial and theta domain encompasses adequate dynamical range to fully accommodate the launching of magnetized disk wind. We also consider cylindrical coordinates, with R = r sin θ and z = r cos θ, to better describe simulation setup and diagnostics.

Disk Model
The disk model is scale-free and is largely identical to that of Bai & Stone (2017) and Cui & Bai (2020). The density profile comprises of separate radial and theta components where ρ0 is the midplane density at r = R0, and we set qD = 2. The function f (θ) is computed by solving hydrostatic equilibrium based on a prescribed equilibrium temperature profile Teq, which is specified as where T0 = v 2 K0 = GM/R0 is the midplane temperature at r = R0. The aspect ratio is ǫ(θ) = H/r, where H = cs/Ω is the pressure scale height, and c 2 s = P/ρ is the isothermal sound speed. The relation between f (θ) and ǫ(θ) can be found in Bai & Stone (2017); Cui & Bai (2020) by solving hydrostatic equilibrium. We choose the temperature profile ǫ(θ) to be which transitions from a disk value ǫ d = 0.1 to a surface value ǫw = 0.5 smoothly. We thus quote H d = 0.1R as the disk scale height. The transition peaks at θtrans = 3.5H d /r and is mostly achieved within a width of 0.05 radian, with δθ = |θ − π/2| − θtrans the angle that deviates from θtrans. Such choice mimics additional heating from UV and X-rays in the disk atmosphere (e.g., Glassgold et al. 2004;Walsh et al. 2010). The resulting equilibrium rotation profile is given by while radial and meridional velocities are zero at equilibrium. On the other hand, to trigger the MRI, we apply random noise to these velocity components with an amplitude of 5% of the local sound speed. We initialize magnetic fields only for poloidal component by prescribing an azimuthal vector potential (Zanni et al. 2007): The poloidal field is computed by B = ∇ × A φ , allowing preservation of divergence-free condition. In the above, Bz0 is the midplane field strength at R0. The parameter m reflects the degree to which initial poloidal field lines bend. We choose m = 0.5, corresponding to modestly bent initial fields. At the midplane, we obtain a pure vertical field The strength of the poloidal fields is characterized by plasma β, the ratio of gas to magnetic pressure, which is constant at the disk midplane in our prescription. We choose a midplane value of β0 = 10 4 in most of our simulation models. We relax disk temperature to initial temperature through the Λc term on the right hand side of Equation (4). After each time step over a time interval ∆t, we adjust gas temperature T to T + ∆t with where the relaxation time τ is prescribed as a fraction of the local Keplerian orbital period, τ (R) ∝ P orb ≡ 2π/ΩK(R). To avoid the excitation of potential hydrodynamic instabilities in the disk, primarily vertical shear instability (VSI), we choose τ = P orb in the disk region. We further make τ smoothly transition to 0 in the wind zone (i.e., locally isothermal) by strictly setting to equilibrium values at each time step. This allows the system to better relax to a quasi-steady state based on our numerical experiments.

Elsässer Numbers
The dimensionless Elsässer numbers quantify the strength of non-ideal MHD effects. For AD, it is given by where vA = B 2 /4πρ is the Alfvén speed. Given an ionization fraction, ambipolar diffusivity typically scales as ηA ∝ B 2 /ρ 2 , and hence Am does not depend on the field strength.
In our simulations, Am is set to a constant throughout the bulk disk (Table 1), and transitions to Am = 100 from the disk zone towards the wind zone to mimic the elevated ionization level above the disk surface mainly by stellar Xrays and far-UV (FUV, Perez-Becker & Chiang 2011). The transition of Am follows the same form as that of temperature, since they share the same physical origins. Correspondingly, we consider the transition height, characterized by θtrans = 0.35, to be the "wind base", as identified in earlier studies Gressel et al. 2015).
Besides AD, we also incorporate Ohmic resistivity near the inner boundary for pure numerical considerations. The resulting magnetic dissipation allows poloidal magnetic fields to stay straight, which helps stabilize the inner boundary. The value of midplane resistivity is fixed to ηO = 0.02csH or 0.03csH at R = R0 and is gradually relaxed in radial direction to approach approximately zero at R = 2. In θ direction, ηO transitions from the midplane value to zero around 1.2 radian above/below midplane (near the pole) over a thickness of ∼ H d /r, and we strictly require both AD and resistivity to vanish at two active cells from the pole.
The Hall effect can be marginally important at the midplane region up to ∼ 30 − 60 AU (Bai 2015). Nevertheless, it was found in local simulations to have only minor impact on angular momentum transport. Hence, as a first study, we do not include the Hall effect to avoid potential complications.

The Mesh and SMR
In spherical-polar coordinates, the most severe limitation to simulation timestep is the small azimuthal grid size in the polar region. Taking the advantage of hierarchical mesh structure in Athena++, we use static mesh refinement (SMR), which allows us to use coarse grid in the polar region, and in the meantime to zoom in the disk zone to achieve adequate resolution.
All 3D simulations carried out in this paper are equipped with the same grid structure. At the root level, we set 144 × 64 × 12 active cells in (r, θ, φ). The radial domain spans from r = 1 to r = 100 with logarithmic grid spacing. The grid spacing in θ-direction increases by a constant factor of 1.05 Figure 1. Left: the 3D mesh that covers the entire simulation domain. Right: a zoom-in mesh in R, z-plane of inner disks. The meshes in both panels are plotted by every 4 cells in radial direction and every 2 cells in θ direction for illustrative purposes. Note that due to visualization artifacts, straight line segments in θ and φ should be considered as arcs.
per root grid cell from the midplane towards the two poles, yielding a resolution at the midplane about five times finer than that at the pole. The grid spacing in the φ direction is uniform.
On top of the root level, we set three levels of refinement. The mesh structure is shown in Figure 1. The left panel illustrates the mesh of the entire simulation domain, and the right panel zooms in to the inner parts of the disk. In the radial direction, the finest level spans over r ∼ 3 − 15, where our data analysis are conducted within this range. In θ-direction, the finest level n = 3 covers ±H, the n = 2 level occupies ±2.5H, n = 1 level fills ±6H. Overall, this grid setting enables us to achieve a resolution of 32 cells per H in r, 30 cells per H in θ at the disk midplane, and 12 cells per H in φ in the n = 3 level.

Boundary Conditions
We choose a "fixed-state" type of inner boundary condition, which we find is able to sufficiently stabilize the inner boundary so that it does not strongly interfere with the bulk computational domain (see further discussions in Section 6.1). Specifically, we impose temperature to be equal to their equilibrium values from Equation (9). Density is extrapolated from the last active zone, assuming ρ ∝ r −qD . The angular velocity is set to the minimum between equilibrium v φ and the angular velocity of solid body rotation ΩK(R0)R. The other two velocity components are vr = v θ = 0. The addition of resistivity mentioned in Section 2.3 further helps stabilize the inner boundary by creating less turbulent environments and straightening field lines. In some simulations which yield strong turbulence, we further enforce the poloidal magnetic field lines penetrating through the inner boundary to be tied to its initial position (no sliding) by resetting the φ− electric field E φ = 0 at the inner boundary.
Other boundary conditions are standard. At outer radius, the hydrodynamic variables are extrapolated from the last active zone, assuming ρ ∝ r −q D , T ∝ r −q T and v φ ∝ r −1/2 . Radial and meridional velocities vr and v θ are copied directly from the last active zone, except setting vr = 0 when vr < 0 (standard outflow boundary condition). Magnetic field variables of the inner and outer ghost cells are set based by Our θ-domain directly reaches the poles, and we employ polar wedge boundary conditions (Zhu & Stone 2018). Lastly, boundary condition in φ-direction is periodic. Table 1 lists our simulations models. Model names with an 'E' indicates those we enforce E φ = 0 at the inner boundary. In this work, two physical parameters are under investigation, the ambipolar Elsässer number Am and disk magnetization β0. The parameter spaces cover β0 ∈ {10 3 , 10 4 } and Am ∈ {0.1, 0.3, 1, 3, 10}. In addition, all simulations are run up to 3000P0, where P0 is the orbital period at the inner boundary R0 = 1. We choose fiducial model to be Am1 with Am = 1 and β0 = 10 4 . This is due to the fact that the Ambipolar Elsässer number is found to be on the order of unity towards outer PPD regions (Bai 2011a,b), and controls the onset of MRI with a threshold ∼ 1 (e.g., Bai & Stone 2011). Moreover, the strength of disk magnetization that yields an accretion rate similar to the typically observed values is on the order of β0 = 10 4 (Simon et al. 2013a;Bai 2015).

Mass Accretion
In this section, we present some relevant diagnostics that will assist data analyses in result sections. In cylindrical coordinates, the φ−component of momentum equation (2) can be cast to expresses the conservation of angular momentum: where j = v φ R ≈ vKR is the specific angular momentum, ρdz is the disk surface density, angle brackets · denote azimuthal and temporal averages. We truncate the disk at z ± = ±0.32R, which is slightly smaller than the expected position of tan(0.35)R. This is chosen because we impose the transition over a finite thickness (instead of a sharp jump) of ∼ ǫ d /2 = 0.05 in θ and we do observe that wind can be launched at a lower height.
The accretion rate is defined aṡ and the Rφ and zφ components of the stress tensor are given by where δvR and δv φ are the radial and azimuthal velocity fluctuations. The hydrodynamic and magnetic parts in the above are known as the Reynolds and Maxwell stresses, respectively, which we denote as T Rey and T Max . In Equation (16), the first term corresponds to the change of angular momentum due to the time evolution of surface density, the second term is the advection of angular momentum from the accretion flow, the third and last terms are the driving forces of disk angular momentum transport, giving the fluxes of angular momentum in the radial and vertical dimensions, respectively.
We note that from the continuity equation, we have whereṀ loss (R) is the cumulative disk mass loss rate within radius R, and it can be more conveniently expressed as where note that we use log R to ensure the dimension in the above equation to match that of the accretion rate. By substituting (19) to Equation (16), we obtain 2πRΣ ∂j ∂t +Ṁacc ∂j ∂R Here, only the Maxwell stress is retained for vertical transport of angular momentum. In quasi-steady state, we anticipate the specific angular momentum j ∼ vK R to be stationary in space and time, with ∂j/∂R = vK /2. Therefore, in our standard analysis, we divide the above equation by vK /2, so that the second term becomesṀacc if j = vK R. This approach, which we refer to as "original" decomposition, directly illustrates the expected contribution to accretion from viscous stress (3rd term) and wind (last term). There is also the first term, which reflects the time-variation in specific angular momentum. We call it time-dependent correction, and it can arise when disk forms ring-like substructures which modifies the rotation profile.
In the meantime, as the formation of disk substructure makes rotation profile deviate from Keplerian, implies that ∂j/∂R does not necessarily be equal to vK /2. To obtain the actual contribution from these physical process toṀacc, one needs to, on top of the "original" decomposition, multiply this equation by a dimensionless correction factor Corr(R) ≡ (vK /2)(∂j/∂R) −1 . We call this approach "actual" decomposition.
The radial transport of angular momentum is mainly mediated by turbulence in the disk, and is characterized by the classic dimensionless α parameter (Shakura & Sunyaev 1973), Over the disk column, the vertically averaged α is obtained by The vertical transport of angular momentum is mediated by the magnetized disk winds. Assuming magnetized winds drive disk accretion in steady state with j = vK R, the local mass flux can be written by (Wardle 2007 where fast accretion occurs at steep vertical gradient of B φ . This relation is derived in laminar state, while we will examine our results under turbulent conditions by taking averages using ρvR , Bz and B φ .

Units Conversion
Since the code units are scale-free, here we present the conversion of code units to physical units for accretion rateṡ M and magnetic field strength B. Given a radius in code units, we aim to find the corresponding mass accretion rate and magnetic field strength in realistic PPDs. Consider at a physical radius of 30 AU, the surface density is 10 g cm −2 as for minimum-mass solar nebula (Weidenschilling 1977). The mass accretion rate in physical units is expressed bẏ and the magnetic field strength is given by In the above, subscripts "phys" and "code" represent quantities measured in physical units and code units, respectively. These results will be utilized later for discussions.

Characterizing Rings and Gaps
Later in §4.2 and §5, we will see the disk form annular substructures with variations in surface density profile. For the ease of analysis, we define the amplitude of surface density variation (depth δΣ) of a ring and a gap as which is the difference between the local maximum and minimum of surface density divided by the minimum. The width wΣ of a ring or gap is measured as the half distance of two adjacent local maxima in surface density Σ, Note that our definition of the depth and width are purely based on gas surface density, which can be directly obtained from simulation data. In observations, since the surface density is not directly observable, the gap depth is usually defined as the difference between the intensity at the gap location and the peak value of its outer ring, and the gap width is as the full width at half depth (e.g. Long et al. 2018). The definitions of gap depth and width in surface density also vary in the literature (see Section 3.2 in Dong & Fung 2017). In this paper, we employ simple definitions of the depth and width in order to give a semi-quantitative description of our simulation results instead of directly comparing with observed characteristics of rings and gaps.

THE FIDUCIAL MODEL
In this section, we focus on the fiducial model Am1, following by detailed analyses on the overall disk evolution ( §4.1), magnetic flux concentration ( §4.2), magnetic flux transport ( §4.2.2), turbulence and stress levels ( §4.3), angular momentum transport ( §4.4), and azimuthal features ( §4.5). We present time in units of the innermost orbital time P0 ≡ P orb (R0) = 2π throughout the result sections. Based upon the above discussion, we choose to analyze data over a time domain t ∈ [2600, 3000]P0 and focus on a spatial domain R ∈ [5, 10] if not otherwise noted. We quote that 2600P0 is approximately 233 and 82 local orbits at R = 5 and R = 10, respectively, and 3000P0 corresponds to 268 and 95 local orbits. The MRI is well developed to a radius of R ∼ 15 over the time interval selected.

Overall Disk Evolution
The time evolution of the fiducial Model Am1 is illustrated in Figure 2. We show snapshots of radial mass fluxes (top), poloidal (middle) and toroidal (bottom) magnetic field strength at different innermost orbital times t = 0, 1000, 2000, 3000P0. The overlaid black contours represent equally spaced lines of constant magnetic flux Φ, which represents the cumulative magnetic flux from the north pole, calculated by or equivalently, by the divergence free condition ∇ · B = 0, where Br and B θ are azimuthally averaged, and Φ(r0, θ) is the magnetic flux at the inner radial boundary. At t = 0P0, we see the smooth profile of poloidal field imposed throughout the disk on top of initial velocity perturbations. The large-scale poloidal field threading the disk results in launching of magnetized winds from surface layers. The Keplerian rotation winds up the poloidal field to generate the toroidal magnetic field, and it dominates the magnetic field strength over the entire course of the simulation. At t = 1000P0, for a large extent of the disk (R 6), toroidal fields are symmetric but possess opposite signs about the midplane, generating a thin current layer there and driving the gas accretion. The accretion flow can be clearly seen in the top panel for radial mass flux, where efficient inward transport of mass is discernible at the midplane, consistent with previous global 2D simulations (e.g. Suriano et al. 2018;Bai & Stone 2017;Gressel et al. 2020). The MRI turbulence sets in progressively from inner regions towards larger radii. At the inner part of the disk (R 6), the MRI turbulence leads to the disruption of current sheet at the midplane.
As the system evolves, at t = 2000P0 and t = 3000P0, the disk becomes fully MRI turbulent up to a radius of R ∼ 10 and 15, as seen from the disruption of the midplane current sheet. Field strength is still dominated by the toroidal component in the bulk disk. As the mean toroidal fields flip at arbitrary vertical locations, the resulting accretion flow becomes randomly distributed in the disk. We also see that the distribution of poloidal magnetic flux within the disk becomes non-smooth, as opposed to the initial smooth profile. Despite the disk is fully turbulent, the field configuration in the wind region remains smooth.
In the left two columns of Figure 3, we show diagnostic quantities at the end of simulation, averaged over the last 400P0, including plasma β at z = 0 and z = 3H, gas and magnetic pressure, and three components of magnetic field strength. We also note that the MRI is mainly driven by the vertical field, as found in previous local simulations (e.g., Bai & Stone 2011;Simon et al. 2013a,b). These radial and vertical profiles will be discussed in the subsequent subsections to facilitate the analysis.

Magnetic Flux Concentration and Evolution
In the second panel of Figure 2, we see that within the disk, the vertical magnetic flux accumulates, assembling into magnetic flux sheets at separated radii. This phenomenon, known as magnetic flux concentration, has been reported in a variety of MHD simulations for PPDs (e.g. Bai Figure 2. Snapshots of azimuthally averaged radial mass flux and magnetic fields for Model Am1 at t = 0, 1000, 2000, 3000P 0 from left to right. Top panels: snapshots of radial mass flux ρvr, overlaid with green arrows of velocity vectors. Middle panels: snapshots of scaled meridional magnetic field RB θ . Bottom panels: snapshots of scaled toroidal magnetic field RB φ . Overlaid black curves are equally spaced contours of poloidal magnetic flux. Dashed curves mark an opening angle of 0.35 above and below the midplane, which correspond to the transition from non-ideal MHD dominated disk zone to ideal MHD dominated wind zone. Riols et al. 2020). We note that, however, these existing studies are either local simulations, or 2D and 3D global simulations without properly resolving the MRI. Thus, our simulations can serve to provide robust evidence of the magnetic flux concentration phenomena, which will be analyzed in further details in this subsection and §5.2.

Riols & Lesur 2019;
In Figure 4, we show the evolution of Bz/Bz0 at the disk midplane in a space-time diagram, where this quantity is azimuthally and vertically (±1H) averaged. 2 The overlaid white contours delineate linearly equally spaced lines of constant magnetic flux Φ(r, π/2) at the midplane. From Figure  5 10 2 and Figure 4, we see that the development of magnetic flux sheets is associated with regions where the MRI turbulence is developed. As such regions expand towards larger radii over time, so do regions of magnetic flux concentration. At 3000P0, discernible flux sheets emerge at R ≃ 5.5, 8.5, 11, 13, with high level of flux concentration. In the bulk disk, the peak vertical field strength reaches a factor of ∼ 4 compared to its initial value, whereas in between the flux sheets it can be reduced to close to zero. This is also manifested in the radial profile of plasma β in the top left panel of Figure 3, where local minimum of β are obtained at the locations of flux concentration. We also see from the space-time diagram that the flux sheets are dynamically evolving in our fiducial Model Am1, where thicker sheets can split into multiple thinner sheets, and thin sheets can also merge.
In our fiducial model Am1, the typical thickness of the flux sheets is ∼ 0.5H, and the flux sheets are typically separated by ∼ 2 − 5H. Local shearing-box simulations that resolve the MRI turbulence show comparable thickness of the flux sheets (Bai & Stone 2014;Bai 2015;Riols & Lesur 2018). However, the separation between flux sheets is more ambiguous in these shearing-box simulations due to their local nature, where flux sheets are either too closely spaced with a separation of ∼ 0.5H in unstratified simulations (Bai & Stone 2014), or form only one single flux sheet over a simulation domain of several H (Bai 2015). Such ambiguity is resolved in our global simulations.
In Figure 5, we further show the radial surface density profile Σ(R) at t = 0, 600, 1200, 1800, 2400, 3000P0, as well as radial profiles of Bz/Bz0 at the midplane for all models. Focusing on the fiducial model Am1, we see that the overall surface density profile deviates from the initial profile, and develops multiple density maxima and minima. In particular, the surface density minima are co-located with where magnetic flux is concentrated, in agreement with previous local and global simulations (e.g. Bai  . Surface densities Σ (first row and third row) as a function of radius at t = 0, 600, 1200, 1800, 2400, 3000P 0 from light to dark colors of all simulation models. Vertical magnetic fields normalized by initial values B Z /B Z0 (second row and forth row) as a function of radius at t = 1200P 0 (red) and t = 3000P 0 (black). Béthune et al. 2017;Suriano et al. 2019;Riols et al. 2020). The width of the rings or gaps is wΣ ∼ 1.5−2.5H. The surface density variations, estimated based on Equation (27), span over δΣ ∼ 15 − 50% for Model Am1. We defer further discussion on the physical mechanisms that produce density contrast and magnetic flux concentration phenomena in §4.4.1 and §6.2.
We note that the inner disk is gradually depleted over time (see §4.4), which in reality may affect the ionization and thermal structure of the disk. Nevertheless, we neglect such further complications in this study, and our simulations are terminated before depletion becomes drastic.

Toroidal magnetic field evolution
We see from the middle panels of Figure 2 that the poloidal fields in the disk region, especially the flux sheets, are not purely vertical, but become locally tilted in a random manner. As the sign of toroidal field anti-correlates with the sign of the radial field component, the mean toroidal fields show arbitrary sign flips across the disk. The orientation of these tilts are dynamically evolving together with the flux sheets. This is accompanied by sign changes in toroidal field. In the second column of Figure 3, we see that the toroidal magnetic fields dominate over the poloidal components across the disk, with sign flips of B φ also identified as a function of radius at the disk midplane.
In Figure 6, we show the time evolution of the vertical profiles of azimuthaly-averaged toroidal magnetic fields B φ at r = 5.5. The space-time diagram is often employed in the context of MRI dynamo actions, where the mean toroidal fields flip every 10 orbital periods in ideal MHD simulations of the MRI with zero net vertical magnetic flux (e.g., Davis et al. 2010;Shi et al. 2010;Flock et al. 2011), but the dynamo actions are weakened or even suppressed by net vertical magnetic flux with or without the presence of non-ideal MHD effects Simon et al. 2013a,b). The pattern in the space-time diagram of B φ seen in Model Am1 is clearly not a standard MRI dynamo, but rather just reflects the orientation of the tilt in poloidal fields. The sign flip at about 1850P0 is mainly associated with a change from BR > 0 to BR < 0 in the bottom half of the disk. Comparing to local simulations with similar setup (Simon et al. 2013a), our results share similar irregular sign changes in the bulk disk, but lacks dynamo-like sign flips in the FUV layer (θ > θtrans). This is likely due to global conditions imposed by wind-launching from the disk surface, fixing the tilt of poloidal field and hence sign of B φ .

Magnetic Flux Transport
Besides magnetic flux concentration, we further examine global trends of flux evolution. This is known as magnetic flux transport. As disk evolution is largely controlled by the amount of magnetic flux threading the disk, flux transport governs the long-term evolution of the disk at a more fundamental level. The white contour lines in Figure 4 show the time evolution of constant magnetic flux surface. It can be seen that after a rapid initial relaxation, the magnetic flux is slowly transported radially outward at a steady rate. Later on, as flux concentration develops, the location of constant flux surface tracks the location of the flux sheets where most of the flux assembles.
To quantify the radial transport of poloidal magnetic flux in Model Am1, we show the evolution of Φ/Φ0, where Φ0 = Φ|t=0, for all simulation models at R = 10 in Figure 7. It illustrates the fraction of magnetic flux enclosed within R = 10 compared to initial flux distribution. As we have discussed, the MRI is fully developed at this radius after t = 2000P0. Before this time, the enclosed magnetic flux steadily declines, corresponding to outward flux transport. Assuming the flux transport outward with such a steady rate, it requires ∼ 800 local orbits to lose half of the initial magneitc flux enclosed within R = 10 for the fiducial model. We have verified by performing additional axisymmetric 2D simulations with identical simulation setup and found the same rate of flux transport. This rate is much lower than that reported in Bai & Stone (2017) (only about ∼ 10%), which we speculate is related to the different simulation setup. In particular, wind in Bai & Stone (2017) is launched from a lower height with θtrans = 0.3, and it is found that the rate of flux loss is sensitive to θtrans under the current setup, where small transition angles correspond to fast magnetic flux losses (Yang & Bai, to be submitted).
After t = 2000P0, magnetic flux transport is affected by flux concentration. For our fiducial model Am1, we see that the rate of flux transport stalls after t ∼ 2400P0. This is associated with the fact that despite exhibiting dynamical evolution at earlier time, the locations of major flux sheets appear to be relatively steady. As a large fraction of magnetic flux is trapped in flux sheets, they serve to protect the system from further losing magnetic flux. On the other hand, we find a diversity of behavior from other models, to be discussed in Sections 5.2 and 5.3.

Turbulence and Stress Levels
To parameterize turbulence levels induced by the MRI, we compute the root mean square of velocity fluctuations, where angle bracket v denotes azimuthal averaging. The third column, top panel of Figure 3 presents the mean δv/cs as a function of radius, where the spatial average of δv/cs is performed over ∆R ∈ [0.95R, 1.05R] and ∆z ∈ [−3H, 3H], and the time average is over the last 400P0. The turbulence is dominated by radial and azimuthal components, with strengths of δv ∼ 0.02cs for the radial extent 5 < R < 10. The turbulence levels gradually decline over R > 13 as the MRI has not yet fully saturated in this region. The middle and right panels show the vertical profile of δv at R = 5.5 (with flux concentration) and R = 6.5 (with less flux concentration), averaged over ∆R ∈ [0.95R, 1.05R] and the last 400P0. The turbulence strength is on the order of 0.01cs towards the disk midplane and rise towards high latitude approaching 0.1cs at the surface layers at both radii. This may suggest a turbulent FUV layer as seen in local shearing-box simulations (Simon et al. 2013b). The difference between regions with and without flux concentrations is moderate, which can be up to a factor of ∼ 2 in δv/cs. The forth column of Figure 3 show the Reynolds and Maxwell stresses defined right below Equation (22). The quantities are temporally and spatially averaged in the same way as for turbulence levels δv/cs. We see that Maxwell stress dominates the Reynolds stress over the entire radial extent. The dimensionless α parameter is a few times 10 −4 for Reynolds stress and ∼ 1 − 2 × 10 −3 for Maxwell stress in the radial domain 5 < R < 10. Similar to turbulence levels, the difference between regions with and without flux concentrations of stresses is also moderate, up to a factor of ∼ 2 − 3 in α parameter.

Angular Momentum Transport
We discuss angular momentum transport and mass loss in the fiducial model in Figure 8. Following Section 3, we show in the top panel the "original" decomposition of terms in Equation (21), which reflects the torques exerted by various transport mechanisms. The data are averaged over the last 300P0, except for the time-dependent term, which is computed from the finite difference between data from the first and second halves over this period, and gives an approximate estimate on the sign and magnitude of the term over a timescale of ∼ 150/R 3/2 local orbits. At each radii, we have further smoothed the data over ±0.5H.
The accretion rate is on the order of 10 −4 in code units, and is generally consistent with that expected for wind-driven  (21), as well as the mass loss rate. Bottom panels: actual accretion rateṀacc by multiplying a correction factor Corr(R) to the top panel. See §3 for definitions. Dashed curves denote negative accretion rates. (25), this accretion rate corresponds to 10 −8 M⊙ yr −1 in physical units for typical outer disks. We see that the dominant contribution is from vertical transport by the wind, but radial transport, as well as the time-dependent term, also play a significant role, and can reach a level comparable to wind-driven accretion at certain locations. Moreover, Maxwell stress is found to be a major contribution to the viscously-driven accretion, and the Reynolds stress has typically negligible contribution at all locations, consistent with expectations discussed in Section 4.3.

accretion. Recalling Equation
There is very significant mass loss in our simulations, with dṀ loss /d log R comparable or even exceeding the accretion rate. Such significant mass loss is generally expected for magneto-thermal disk winds (Bai et al. 2016), where thermal pressure plays an important role and has been observed in a number of simulations of the kind (e.g., Bai & Stone  2017; Cui & Bai 2020;Rodenkirch et al. 2020). On the other hand, the exact mass loss rate is sensitive to the thermal chemistry in the wind launching region, where our treatment is considerably simplified, and more realistic simulations incorporating these physics tend to yield more mild mass loss rates Gressel et al. 2020). Therefore, when discussing wind mass loss, we focus on the trends instead of absolute values. The significant mass loss can alter the surface density profile, potentially leading to the depletion of the inner disk (e.g., Suzuki et al. 2016), and this is indeed seen in our simulations. Nevertheless, we will not pursue further discussion on this density depletion due to unrealistic mass loss in our simulations.

Substructure formation
Magnetic flux concentration leads to variations in winddriven accretion and mass losses. We see that the wind stress BzB φ is directly correlated with magnetic flux distribution when comparing with the second panel of Figure 5, where regions with flux concentration has on average stronger stress. However, the radial variation in the BzB φ stress is not as dramatic as magnetic flux concentration seen in the midplane. This is readily visible from the middle row of Figure 2, where magnetic flux becomes less concentrated in the surface layer. The distinction of magnetic flux concentration in the ADdominated regime from the ideal MHD regime was observed in local shearing-box simulations (Bai & Stone 2014), where only the former leads to thin flux sheets (see also our discussion in Section 5.2). Therefore, we expect the thickening of flux sheets in the wind zone is likely associated with the transition into the ideal MHD regime.
Interestingly, we notice that peaks in the measured mass loss rate profile do not coincide, but tend to anti-correlate with those in the stress profile. This is somewhat counter-intuitive as one generally expects stronger net vertical field to yield both higher wind-driven accretion and mass loss rates. This expectation is based on the assumption that gas is well coupled to magnetic field and hence its streamlines follows poloidal magnetic fields. As we implement the transition from non-ideal MHD to ideal MHD regime over a relatively broad (∼ H) layer, we do observe departures between flow streamlines and poloidal field lines over the transition region. More systematic studies are needed to examine how the mass loss profile depends on the smoothness of this transition.
The formation of ring-like substructures further complicates the process of angular momentum transport. They are more strongly exhibited in radial transport of angular momentum, as well as the time-dependent correction term. Because radial transport of angular momentum depends on the radial gradient of the Rφ stress, radial variations in surface density as well as the α parameter seen in Figure 3 give rise to frequent sign changes in its contribution to gas accretion. In addition, the mass accumulation and/or dilution in rings/gaps also changes the level of radial pressure gradient and hence rotation profiles, which is at the expense of angular momentum and leads to the first term in Equation (21), i.e., the time-dependent correction. Contribution from these terms, coupled with ring-like substructure, give rise to additional variations in the total accretion rates, as can be readily tracked in Figure 8. These terms could be considered as non-linear back-reaction to disk angular momentum transport in response to driving from magnetic-flux concentration.
The bottom panel of Figure 8 shows the actual accretion rate profiles, as well as contributions from individual terms. This is obtained by correcting for the ∂j/∂R factor. We see that after this correction, there are more frequent radial variations at relatively small amplitudes in wind-driven accretion rates, but they barely show any correlation with the profile of poloidal magnetic flux. In fact, the amplitudes of rings and gaps are dynamically evolving, reflecting the complex interplay among different transport mechanisms together with the aforementioned "back-reaction" due to formation of disk substructure. Altogether, we see that despite substantial magnetic flux concentration in the bulk disk, the disk rarely form deep gaps.

Vertical profiles of the accretion flow
In the left panel of Figure 9, we examine the vertical structure of the accretion flow by showing the azimuthally and temporally averaged radial mass flux profiles at R = 5.5 for the fiducial Model Am1. The quantities are averaged over last 200P0. We see there are both accretion and decretion flows of comparable level of mass flux present at different heights, thus the net accretion flow that we discussed earlier results from cancellations from such flows. By comparing the results with Equation (24), we see that the mass fluxes are generally consistent with predictions suggesting that the overall flow structure is primarily driven by the torques exerted by magnetic forces. In particular, it is set by the vertical gradient of the toroidal field, and following the discussion in Section 4.2.1, the locations of the accretion/decretion flows are random and dynamically-evolving.
We also see that the mass flux profiles show departures from predicted by Equation (24). Such departures can be substantial, sometimes a factor of a few in the peaks. This reflects the role played by the MRI turbulence. Following our previous discussion, contributions to such departures include radial transport and time-dependent correction, which is in part associated with formation of disk substructures.
The right panel of Figure 9 illustrates the radial velocity normalized by local sound speed as a function of height, at radii of R = 5.5 and R = 6.5, showing alternating radial velocity on the order of ∼ 10 −2 − 10 −1 within ±2H. We anticipate that such alternating radial gas flows driven by disk winds can effectively enhance the dust radial diffusion due to the weak MRI turbulence alone, as shown in Hu & Bai (2021). It would be interesting to explore this subject further by incorporating a dust component, as the level of radial diffusion in dust, together with the profile of the gas rings, sets the width of dust rings in mm-continuum observations .

Azimuthal Features
The discussions above are based upon azimuthally averaged simulation data. It remains to address whether these features are largely axisymmetric or involve additional nonaxisymmetric structures.
In Figure 10, we show snapshots in (R − φ) plane of surface density, scaled meridional and azimuthal magnetic fields, and vorticity at the midplane. In the second panel, the meridional field B θ is mainly axisymmetric, indicating that the vertical magnetic flux is indeed confined in flux sheets rather than an artifact of azimuthal averaging. Likewise, the toroidal field B φ shown in the third panel is also highly axisymmetric, which exhibits sign changes in R, similar to sign changes in z and time as previously discussed. We note that toroidal fields reverse polarity sharply, leading to steep gradients at certain radii. This phenomenon could be attributed to the sharp structures resulted from the anisotropic diffusive properties of AD (Brandenburg & Zweibel 1994).
The formation of ring-like substructures may induce vortex formation by Rossby wave instability (RWI) (Lovelace et al. 1999;Li et al. 2000). In the last panel, we exmamine the z−component of midplane vorticity, defined as where axisymmetry and Keplerian rotation give ωz/ΩK = 0.5. We note that no discernible vortex formation is observed, indicating the absence of the RWI. Ono et al. (2016) calculated the fastest growth rates as a function of amplitudes and widths of the surface density variation (see their Figure 6).
Taking the Gaussian bump model for surface density variations in their paper, and given a width of wΣ ∼ 1.5 − 2.5H for our fiducial Model Am1, the smallest amplitude of surface density variation to trigger the RWI is on the order of ∼ 100 − 200%, which is achieved with the azimuthal wavenumber m = 1. Seen in our Figure 5, fiducial model Am1 results in variations of δΣ ∼ 15 − 50%. Despite our simulations have a reduced azimuthal domain, the critical surface density variations found in Ono et al. (2016) is even higher for larger azimuthal wavenumber m. This is consistent with the lack of vortices in the last panel.

PARAMETER STUDY
In this section, we present the results for Model B3, Am0.1, Am0.3, Am3E, Am10E, where §5.1 focuses on flow properties including turbulence levels, stresses, and angular momentum transport, §5.2 discusses magnetic flux concentration and transport, and §5.3 details each individual simulation model.

Flow properties
Similar to Figure 3 for the fiducial model Am1, we show in Figure 11 and Figure 12 the vertical and radial profiles of turbulence levels and stresses for the rest of the simulation models. The averaging is taken in the same way as in Figure 3 (see §4.3). In Figure 11, we see that the turbulence levels increase with larger Am and stronger magnetization (smaller β0). The values of δv/cs gradually raise from ∼ 10 −2 for Model Am0.1 to ∼ 10 −1 for Model 10E as well as for Model B3. Correspondingly, the Reynolds stresses α Rφ increase from ∼ 10 −4 to ∼ 10 −3 , and the Maxwell stresses from ∼ 10 −3 to ∼ 10 −1 . Unlike the fiducial model Am1, we do not observe significant boost of turbulence levels at the wind region for the rest of simulation models, but rather, turbulence level is largely flat over height.
Observations on turbulent motions through non-thermal line broadening infer a turbulent velocity δv < 0.04−0.06cs at different heights (Flaherty et al. 2017). The vertical profiles of δv/cs shown in Figure 3 and Figure  Am0.1 Am0.3 Am3E Am10E Figure 11. Turbulence levels δvr , δv θ , δv φ , normalized by local sound speed cs, as a function of height z at r=5.5 (top row) and of cylindrical radius R (bottom row) for Models B3, Am0.1, Am0.3, Am3E, Am10E. models with Am 1 for β0 = 10 4 are favored, which place less stringent constraints than in Simon et al. (2018). Figure 13 shows radial profiles of accretion rates with contributions from Equation (21), similar to the top panel of Figure 8. Angular momentum transport in these models is dominated by magnetized winds, consistent with the fiducial model. Accretion rates are largely set by the net vertical magnetic flux, and is similar among models with β0 = 10 4 . Model B3 possesses an accretion rate with a factor of ∼ 10 higher as expected, since accretion rate scales quadratically with field strength. Contributions to accretion rates from individual terms also show radial variations akin to magnetic flux concentration. The mass loss rates remain substantial, with ratio of d logṀ loss /d log R to accretion rates reaching a few. This ratio is smaller in Model B3, as is expected in disk wind theory for stronger magnetization (Bai et al. 2016). Radial transport by Maxwell stress and time-dependent correction show frequent sign changes that are associated with disk substructures. In sum, these models share qualitatively similar features to those in the fiducial model ( §4.4).

Magnetic flux evolution
All of these simulation models show features of spontaneous magnetic flux concentration (Figure 4 and Figure 14), though their properties, such as the amount of flux sheets emerged, rates of outward flux transport, widths of rings, and surface density variations, are diverse. All models show the formation of rings and gaps in surface density profiles due to the flux concentration ( Figure 5). The widths of rings or gaps span over wΣ ∼ 1 − 5H, and the surface density variations cover δΣ ∼ 30% − 200% for different simulation models.
The space-time diagrams of toroidal magnetic field B φ at r = 5.5 are shown in Figure 6. Similar to the fiducial model Am1, the sign in B φ mainly reflects the orientation of poloidal field through the disk. We note that the toroidal field shows symmetry about the midplane, only flipping sign once across the disk for Model B3, Am3E, and Am10E, while Model Am0.1 and Am0.3 show more stochastic sign changes. This trend could be related to the fact that for stronger net poloidal magnetic fields or weaker non-ideal MHD effects, B φ is more strongly amplified, which is more stiff against reversal.
The overall transport of magnetic flux, as examined in Figure 7, show a two-stage evolution (except for Model Am10E which quickly turns into strong turbulence), with an initial steady decline before the development of the MRI and magnetic flux concentration. Interestingly, after the development of the MRI turbulence, for Model Am1 and Model B3, the flux sheets become stationary and cease to drift outwards. For Model Am0.1 and Model Am3E, magnetic flux keeps drifting outwards, but at an averaged rate comparable to the first stage. Model Am0.3, on the other hand, exhibits faster outward transport with flux concentration. The lack of trend in our results likely reflects the complexity in magnetic flux transport in MRI-turbulent disks.

Notes on individual simulations
In this subsection, we discuss special features in each individual simulation, owing to the diversity of our simulation results.
Model B3 -with β0 = 10 3 , this model is most similar to several previous global disk simulations that report magnetic flux concentration (Béthune et al. 2017;Suriano et al. 2018;Riols et al. 2020). Comparing to previous global simulations (e.g. Suriano et al. 2019), the flux sheets in this simulation are more incoherent and short-lived since the MRI turbulence is resolved. We speculate that it is related to the corrugation of midplane current sheet seen in Figure 6 that oscillates over time. The flux sheets are also thicker than those in the fiducial model Am1, which may be due to stronger turbulent diffusion that likely act against flux concentration.
Three prominent thick flux sheets are present at the end of simulation as seen in Figure 14, resulting in the formation of three gaps in Figure 5. The gaps located at R 6 have widths wΣ ∼ 3H in surface density, and a surface density variable of δΣ ∼ 20% − 70%. We note that a thick flux sheet is developed at R ∼ 4 − 6, driving a wide deep gap. As this region is close to the inner boundary, the gas surface density is significantly depleted over simulation timescale, attributed to the fast accretion and significant mass loss. Hence, the disk is more heavily magnetized. This depletion with trapped magnetic flux merits more in-depth study in the future.
Model Am0.1 -With Am = 0.1, this model is borderline since the MRI transitions to the ambipolar shear instability in the limit of Am ≪ 1 (Kunz 2008;Pandey & Wardle 2012). In Figure 14, we see that two major concentrated flux sheets are developed, locating at 5 < R < 6 and 10 < R < 15. These  two sheets are almost stationary, promoting efficient mass depletion in these regions, leading to relatively deep gaps. Seen in Figure 5, the surface density shows two bumps with widths of wΣ ∼ 4 − 5H and varies by δΣ ∼ 80% − 200%. We have also experimented with different inner boundary conditions (e.g., whether to enforce E φ = 0) while keeping Am = 0.1, finding that the small Am simulations tend to generate fewer but widely-separated flux sheets regardless of the numerics. We also note that in Figure 13, regions with strong magnetic flux concentration show reduced or negative mass loss rate. This is related to the fact that the fiducial location where we set as the wind base is still occupied by the accretion flow (with negative vertical velocity) in this model.
Model Am0.3 -The flux concentration phenomenon present in this model shares similarities with Model Am1 and Model Am3E (Figure 14), where multiple concentrated flux sheets are emerged. Two discernible bumps seen in surface density ( Figure 5) are associated with two sheets with stronger flux concentration located at R = 6 and R = 11. The corresponding widths of rings and gaps are wΣ ∼ 3−5H, and the surface density varies by δΣ ∼ 80% − 130%. In Figure 4, the flux sheets show fastest migration towards large radii among all simulation models, accompanied by rapid flux transport as seen in Figure 7.
Model Am3E -Seen in Figure 4 and Figure 14, the flux sheets in this model are more ordered and particularly steady comparing to other models. These thin flux sheets correspond to a narrower width of wΣ ∼ 1 − 2H in surface density, and a surface density variable of δΣ ∼ 30% − 60%. The magnetic flux is slowly leaking outward comparing to the fiducial model ( Figure 7).
Model Am10E -Model Am10E is close to the ideal MHD regime and is characterized by strongest turbulence among all the other models. As a result, while magnetic flux concentration can be identified, it appears flocculent in the space-time diagram and does not lead to appreciable substructures (Figure 4). The outward transport of flux is fast. The associated surface density does not display prominent annular substructures ( Figure 5). In this model, the features of MRI dynamo are still not observed, in contrast with local shearing-box simulations for which Am = 10 simulation results in a period of dynamo flipping of B φ on the order of 50 local orbits (Simon et al. 2013b).

Numerics
Grid and resolution -We emphasize that conducting 3D simulations with sufficiently high resolution to resolve the MRI turbulence is essential to properly characterize the gas dynamics of outer PPDs. We have experimented with 2D simulations with the same resolution, and the results indicate that the system is steady with no evidence of MRI for the long term evolution. Our results can also be compared with recent global 3D simulations (e.g. Suriano et al. 2019;Hu et al. 2019) that employed coarser resolutions (∼ 9 cells per H), where the flows were considered to be almost laminar without addressing the development of MRI turbulence. On the other hand, the condition to properly resolve the MRI turbulence in non-ideal MHD is yet to be established, especially given the wide range of Am values explored here.
If we adopt the quality factor employed in ideal MHD simulations (Hawley et al. 2011) (see their Equation 10), and replace the most unstable wavelength λMRI with the fitting formula of Bai & Stone (2011) (see their Equation 22), we find the condition in the azimuthal dimension is well satisfied with Qy ≫ 10, and that Qz ≈ 10 for most runs, indicating turbulence is marginally resolved in the vertical dimension.
One caveat of this work arises from the limited φ−domain employed. More extended φ−domain would accommodate more low-m MRI modes to emerge, while restricted φ−domain may artificially amplify turbulence, based on ideal MHD simulations of the MRI turbulence (Flock et al. 2012). At present stage, it is computationally infeasible to scan a large parameter space of Am as in this paper while employing an extended azimuthal domain. We leave extended φ−domain global simulations for a future study.
Inner boundary condition -We have taken special care in setting the inner radial boundary conditions. As we focus on the gas dynamics in the outer PPDs, regions inside the inner boundary corresponds to the inner disk, which is expected to be threaded by poloidal magnetic flux and launches its own disk winds. As this region is artificially truncated by our inner boundary, one expects there is poloidal magnetic flux anchored to this boundary, and would ideally design the boundary conditions that reflect such conditions: it should accommodate mass injection mimicking inner disk outflows at high latitudes, as well as the conventional "outflows" which represents accretion flow in the disk region. Our choice of a "fixed-state" type boundary condition represents a compromise that allows artificial launching of outflows from the high-latitude region of the rotating inner boundary, despite of having some mass accumulation in the disk region of the inner boundary which becomes less realistic after long-term evolution. This is to be compared with the conventional outflowtype boundary condition, which tends to pull gas into the boundary from high altitude, and further attracts and destabilizes poloidal field lines near the inner boundary, although they would have been part of a more stable wind launched from the putative disk regions inside the inner boundary. In ideal MHD simulations, it leads to coronal accretion (e.g., Zhu & Stone 2018), and as field lines penetrating to the inner boundary can extend to large distances, they may substantially affect the overall simulation results. On the other hand, we note that coronal accretion still develops when modeling the inner boundary to mimic a weakly magnetized protostar (Takasao et al. 2018). Such differences highlights the importance of inner boundary conditions in global disk simulations, and merits more systematic studies in the near future.
Thermodynamics and chemistry -We have opted for simple thermodynamics, prescribing the temperature, ambipolar diffusivity, and thermal relaxation timescales. The excessive mass loss via MHD winds observed in all models could be possibly due to these crude prescriptions, and more detailed treatment including ionization thermochemistry and radiative transfer (Bai 2017;Wang et al. 2019;Gressel et al. 2020) shall be considered to properly track the wind launching, acceleration, and propagation.

Mechanisms of magnetic flux concentration
A handful of mechanisms have been proposed to explain the spontaneous magnetic flux concentration seen in global simu-lations. Previously, Bai & Stone (2014) speculated that in the ideal MHD regime, the recurrent action of the MRI channel flows makes oppositely directed fields induced by the channel flows reconnect, producing magnetic loops and separating mass from magnetic flux. Suriano et al. (2018) found that at the disk midplane, toroidal magnetic field reverses polarity, promoting pinching of poloidal field lines and fast accretion. The reconnection of the pinched field acts to separate mass from magnetic flux. Lastly, Riols & Lesur (2019) demonstrated a linear secular instability where initial radial density perturbations lead to transport of gas and magnetic fluxes towards the gap via viscous stress. The excess of poloidal magnetic flux in the gap extracts more mass so that initial density perturbations are reinforced. The presence of the MRI turbulence and diversity of flux concentration phenomena in our simulations complicate the quest for concentration mechanisms. In the following, we give a qualitative phenomenological discussion on the possible causes focusing on the last two scenarios.
Our Model B3 shares similarity with Suriano et al. (2018). Specifically, by closely examining the simulation data, we find small magnetic loops presumably generated by the magnetic reconnection are emerged at the current layers in between the segregated poloidal field lines. For the rest of our simulation models which all possess β0 = 10 4 , however, we do not observe such magnetic loops at the midplane, especially because that most of these models lack persistent midplane current sheet.
The secular wind instability requires, firstly, turbulence to be built in the disk in order to yield viscous stress. Secondly, the magnetic field shall be well-coupled to the gas such that magnetic flux is brought towards gaps together with gas. Lastly, stronger magnetization should lead to an excess of ejected mass through disk winds. Our simulations do possess excessive mass loss, but it is far from straightforward to test whether magnetic flux is brought to the gap by the viscous flow towards filling in the gap. One possible way to testify this mechanism is by increasing the transition heights of the disk θtrans, for which mass loss rate by wind can be reduced, and one can inspect if the spontaneous ring formation is weakened.

Implications for annular substructures
Our simulation results reveal the fact that the spontaneous ring formation through magnetic flux concentration is robust over a wide range of parameter space, including different disk magnetization β0 and ambipolar diffusion strengths. Our results also suggest that magnetically-induced substructures in the MRI turbulence show a diversity of properties. We have seen that the process is stochastic, and the spacing of the substructures to some extent sensitively depends on simulation parameters. Besides the simulations shown in this paper, we have also experimented with other treatments for the buffer zone near the inner radial boundary, including the strengths of resistivity buffer and enforcing E φ = 0. Given the same physical parameters, these simulations share qualitatively similar disk evolutionary patterns and properties of flux sheets, such as stresses and turbulence levels as well as separations and field strengths of the flux sheets. Nevertheless, they differ case by case in terms of exact time and loca-tions in developing these concentrated flux sheets, and their migrations over time.
Another important finding from our results is that the substructures are typically shallow, with surface density contrast being order unity or less, and dynamically evolving. For instance, in our Model B3, Am1 and Am3E, the locations of rings and gaps in the radial surface density profiles are timedependent ( Figure 5). Since observations on ring-like substructures are through dust, it is unclear whether dusts could respond timely. In other words, such shallow substructures may evolve faster, a few tens of local orbital times, than the dust trapping timescales, precluding the formation of deep rings and gaps.
Observationally, by assuming steady state with dust concentration balanced by turbulent diffusion, one can constrain the level of turbulence α/St by separately estimating the width of dust and gas ring width, where St is the dust Stokes number, or dimensionless stopping time characterizing the level of coupling with gas. From the DSHARP disk sample (Andrews et al. 2018), it was found that α/St ∼ 0.1, implying non-negligible level of disk turbulence depending on the not-well-constrained dust Stokes number Rosotti et al. 2020), with upper limit of up to α ∼ 10 −2 . Level of turbulence in our simulations easily satisfies these constraints, and in the meantime, our simulation results call for caution in this approach because the steady-state assumption may not hold.
While rings and gaps are commonly interpreted as a result of interactions between disks and the embedded planets (e.g. Pinte et al. 2020), by resolving key disk microphysics, our simulations reveal that the mechanism drives disk angular momentum transport and evolution is inherently non-smooth, and could result in additional, likely shallower disk substructures. Our magnetic mechanism might be a plausible explanation for substructures found in young Class 0/I disks, for which the evolutionary stage of these systems might be too early for Jupiter-mass planets to form (Sheehan & Eisner 2018;Nakatani et al. 2020;Segura-Cox et al. 2020;Sheehan et al. 2020). For more matured, Class II disks, Jennings et al. (2021) employed a superresolution technique and discovered finer disk substructures in the DSHARP disk samples, including shallow annular substructures in the outer disks. Our results could be a possible interpretation for such shallow substructures.

Implications for B-field observations
Observing magnetic fields in PPDs is challengin. Recent attempts have been made through the Zeeman effect, obtaining only upper limits (Vlemmings et al. 2019;Harrison et al. 2021). These upper limits constrain the instantaneous lineof-site magnetic field strength averaged over the line-emitting region, typically in the outer disk surface layers. Models used to constrain Zeeman observations commonly make simplified assumptions about magnetic field geometry (Brauer et al. 2017;Mazzei et al. 2020), which in our simulations is dominated by the toroidal component. However, the toroidal field randomly reverses polarity in space and time within the bulk disk, which could potentially lead to cancellations of the Zeeman signal and hence complicates the interpretation of observations. On the other hand, we note that for the widely-considered molecule of the Zeeman effect, CN mainly traces disk heights of z/r ∼ 0.2 − 0.4 (Cazzoletti et al. 2018;Teague & Loomis 2020). These locations might be sufficiently high so that the sign of toroidal field does not flip, as in our ficucial model Am1, seen in the second column, top panel of Figure 3. Nevertheless, future works with more realistic treatment of thermodynamics in the wind launching region shall be conducted in order to assess whether there can be sign flips of B φ in the CN-emitting regions.
Another approach to probe disk magnetic field is through paleomagnetism, by recovering remnant field recorded in meteorites when they formed in the solar nebula (Weiss et al. 2021), usually believed to be in the disk midplane region. One major uncertainty in paleomagnetic measurements is that the radial locations where meteorites parent bodies formed are largely unknown, leaving a broad range of plausible radii in the solar nebula. Comparison to theories in the literature have relied on models that assume the disk magnetic field strength is smooth, and link field strength to expected disk accretion rates (e.g. Wardle 2007;Bai & Goodman 2009;Weiss et al. 2021). Fortunately, although poloidal magnetic fields become spontaneously concentrated into flux sheets, the total field strength dominated by toroidal field is relatively smooth, as seen in the second column of Figure 3, which validates this approach.
By scaling numerical units in the simulations to physical units of magnetic fields using Equation (26), we find a magnetic field strength of ∼ 0.01 G for typical outer disk conditions around ∼ 30 AU from our simulation region between R ∼ 5 − 10. We can verify the basic principles used in paleomagnetism, which connects field strength to accretion rates. Within this radii, the measured accretion rates translate to ∼ 2 − 3 × 10 −8 M⊙ yr −1 based on Equation (25). In the wind driven accretion scenario, we employ Equation (3) of Weiss et al. (2021), applying m ≈ 4 (ratio of toroidal field between midplane and wind base), f ′ ≈ 2 (ratio of B φ and Bz at the wind base), and R = 30 AU, we find a magnetic field strength of ∼ 0.008 G, in good agreement with the actual measurements. Both of these results are also consistent with constraints placed by paleomagnetic measurements for the Jupiter family comet 67P, expected to form in the region between 15 to 45 AU, of a local field strength < 0.03 G (Biersteker et al. 2019).

CONCLUSIONS
In this paper, we study the gas dynamics and magnetic field evolution in outer regions of PPDs. To this end, we conduct a set of global 3D non-ideal MHD simulations with ambipolar diffusion, achieving numerical resolutions comparable to those in local shearing box simulations. Our simulations are designed to be scale-free, and we perform parameter study on the strengths of AD and levels of disk magnetization. We summarize our main findings as follows.
(i) Our simulations demonstrate the co-existence of magnetized disk winds and MRI turbulence. While MHD winds dominate the disk angular momentum transport, MRI turbulence also makes considerable contributions.
(ii) For the fiducial model Am1 with Am = 1 and β0 = 10 4 , the turbulence level is on the order of δv ∼ 0.03cs, and the vertically-integrated α parameter is ∼ 10 −3 , dominated by the Maxwell stress. Turbulence and stress levels increase with weaker AD or stronger magnetization, achieving up to δv ∼ 0.1cs, α ∼ 10 −2 for Models Am10E and B3, and down to δv ∼ 0.01cs, but still α ∼ 10 −3 for Model Am0.1.
(iii) Spontaneous concentration of poloidal magnetic flux into largely axisymmetric flux sheets are observed in all of our simulation models, leading to radial variations in turbulence levels, stresses, and mass accretion rates. With angular momentum transport being inherently non-smooth due to magnetic flux concentration, annular substructures arise as a natural consequence.
(iv) In the fiducial model, the flux sheets have typical thickness of ∼ 0.5H, and are dynamically evolving. The widths of the resulting rings or gaps are 1.5 − 2.5H in surface density, with surface density variations span over 15 − 50%, while they show more diverse properties as magnetization and AD Elsässer number vary.
(v) The disk generally loses magnetic flux over time, though flux sheets may serve to prevent the loss of magnetic flux, as in Model B3 and Model Am1.
Our results show that the contrast in surface density variation due to the magnetic flux concentration is on the order of unity or even less, and they evolve as flux sheets stochastically migrate radially in disks. This implies that the magneticallyinduced annular substructures are typically shallow, dynamically evolving, and show diverse properties. It may account for some of the observed disk substructures, such as young disks in Class 0/I phase, for which the evolutionary stage may be too early to form Jovian planets, as well as finer substructures found in super-resolution images of DSHARP disk samples (Jennings et al. 2021).
In surveying a wide range of parameter space, we have chosen relative small azimuthal domain size to reduce computational cost. Larger domain size may be needed to fully capture the global properties of the MRI turbulence (e.g., Flock et al. 2012;Suzuki & Inutsuka 2014). It will also lower the threshold in surface density variation to trigger the RWI (Ono et al. 2016), which may potentially drive vortex formation. In addition, more self-consistent treatment of ionization and thermodynamics/radiative-transfer is needed towards more realistic study and characterization of the gas dynamics in the outer PPDs. These extensions will be explored in future works.