Dynamics Near the Inner Dead-Zone Edges in a Proprotoplanetary Disk

We perform three-dimensional global non-ideal magnetohydrodynamic simulations of a protoplanetary disk containing the inner dead-zone edge. We take into account realistic diffusion coefficients of the Ohmic resistivity and ambipolar diffusion based on detailed chemical reactions with single-size dust grains. We found that the conventional dead zone identified by the Els\"asser numbers of the Ohmic resistivity and ambipolar diffusion is divided into two regions:"the transition zone"and"the coherent zone". The coherent zone has the same properties as the conventional dead zone, and extends outside of the transition zone in the radial direction. Between the active and coherent zones, we discover the transition zone, the inner edge of which is identical to that of the conventional dead zone. The transition zone extends out over the regions where thermal ionization determines diffusion coefficients. The transition zone has completely different physical properties than the conventional dead zone, the so-called undead zone, and the zombie zone. The combination of amplification of the radial magnetic field owing to the ambipolar diffusion and a steep radial gradient of the Ohmic diffusivity causes the efficient evacuation of the net vertical magnetic flux from the transition zone within several rotations. Surface gas accretion occurs in the coherent zone but not in the transition zone. The presence of the transition zone prohibits mass and magnetic flux transport from the coherent zone to the active zone. Mass accumulation occurs at both edges of the transition zone as a result of mass supply from the active and coherent zones.

Low ionization degrees in PPDs induce the non-ideal MHD effects and affect the growth of the MRI.There are three non-ideal MHD effects, Ohmic resistivity (OR), Hall effect (HE), and ambipolar diffusion (AD).The regions where non-ideal MHD effects suppress the MRI are called dead zones (Gammie 1996).Ionization fractions sufficient to drive the MRI are provided only in innermost and outermost radii and upper atmosphere of PPDs owing to thermal ionization, far-ultraviolet and X-ray radiation from the central star, and cosmic rays.
The angular momentum transfer in the inner regions (R < ∼ 10 au) of PPDs where OR plays an important role was investigated in Gammie (1996).He showed that so-called layered accretion is driven around the disk surfaces where MRI is active while the MRI is completely suppressed in the disk interior.Detailed linear analyses of the MRI modified by OR were done by Jin (1996) and Sano & Miyama (1999).MRI-driven turbulence above the dead zone was also studied in local simulations (Fleming & Stone 2003;Turner & Sano 2008;Hirose & Turner 2011) When AD operates in addition to OR, the MRI is suppressed even in the disk surface and the layered accretion disappears because AD works in higher latitudes than OR (Bai & Stone 2013).Instead of the layered accretion, the ϕz component of the magnetic stress owing to global coherent magnetic fields drives laminar surface gas accretion just above the dead zone in the inner region of PPDs where both OR and AD are important (Bai & Stone 2013;Bai 2013;Gressel et al. 2015).
The inner dead-zone edge, inside which the thermal ionization activates the MRI, has attracted attention as a dust accumulation site.In such a region, a local pressure bump is created because the outer part of the active zone moves outward owing to a negative turbulent viscosity gradient.A pressure bump at the inner dead-zone edge traps dust particles in the dead zone that are moving to the pressure bump (Kretke et al. 2009;Ueda et al. 2019).
The inner dead-zone edge has been investigated with global simulations taking OR into account.Dzyurkevich et al. (2010) found that a pressure bump develops at the inner dead-zone edge.Lyra & Mac Low (2012) pointed out that the Rossby-wave instability (Lovelace et al. 1999) induces the formation of vortices at the pressure bump.Flock et al. (2017) performed realistic non-ideal MHD simulations with radiation transfer, and found that a vortex forms at the inner dead-zone edge that casts a nonaxisymmetric shadow on the outer disk.
In the global simulations containing the inner dead-zone edge performed by Dzyurkevich et al. (2010) and Lyra & Mac Low (2012), the angular momentum of the disk is transferred by the MRI turbulence above the dead zone (Gammie 1996).They observed that the height-averaged mass accretion stress does not change significantly across the dead-zone edge.Although sudden decreases in the accretion stresses are found across the inner dead-zone edge in Flock et al. (2017), their range of the polar angle may not be wide enough for the MRI to occur above the dead zone.In addition, in their simulations, the spatial distributions of the OR coefficient are given by simplified analytic formulae with a sharp transition across the inner dead-zone edge.
Inclusion of AD changes the angular momentum transfer mechanism in the dead zone from the MRI turbulence.The MRI above the dead zone will be suppressed once AD is considered (Bai & Stone 2013;Gressel et al. 2015;Bai 2017).In a realistic situation where both OR and AD work, the angular momentum transfer mechanism changes from the MRI turbulence to the magnetic stress due to global coherent magnetic field across the inner dead-zone edge.When both OR and AD are considered, the spatial distribution of the accretion stress across the inner deadzone edge is still unclear.
In this paper, we perform global three-dimensional nonideal MHD simulations of a protoplanetary disk around an intermediate mass star to reveal the structure around the inner dead-zone edge by taking into account both OR and AD whose coefficients are given based on detailed chemistry.In the previous global simulations, a limited range of θ in the spherical polar coordinates of (r, θ, ϕ) is considered, which may affect surface gas accretion expected to be driven above the dead zone (Bai & Stone 2013) and launching of disk winds both from the active zone (Suzuki & Inutsuka 2009, 2014) and the dead zone (Bai & Stone 2013).Our simulation box covers a full solid angle (0 ≤ θ ≤ π and 0 ≤ ϕ < 2π) to capture the angular momentum transfer in the upper atmosphere and vortex formation.
This paper is organized as follows.We will describe the numerical setup and method in Section 2. The main results are presented in Section 3. In Section 4, we will present detailed analyses on transfer of the mass, angular momentum, and magnetic flux in the disk.Astrophysical implications are discussed in Section 5. Our findings are summarized in Section 6.

Basic Equations
The resistive MHD equations with OR and AD are given by ∂ρ ∂t + ∇ • (ρv) = 0 (1) and where ρ is the gas density, v is the gas velocity, B is the magnetic field, ϕ is the gravitational potential of the central star, J is the current density, J ⊥ ≡ J − J • B/|B| is the electric current density perpendicular to the local magnetic field, P is the gas pressure, T is the stress tensor, I is the identity matrix, ηO and ηA are the diffusion coefficients for OR and AD, respectively.They are defined in Section 2.2.2.For simplicity, the Hall effect is neglected in this work.
In this paper, for simplicity we use the locally isothermal equation of state where the gas temperature remains the same as the initial one (Suzuki & Inutsuka 2014;Zhu & Stone 2018).To achieve this, we solve the MHD equations including the energy equation with γ = 5/3, and the gas pressure is derived from the locally isothermal condition where the temperature distribution is shown in Equation (10).We note that a locally isothermal equation of state potentially drives the vertical shear instability (VSI) (Nelson et al. 2013).In the active zone, Latter & Papaloizou (2018) shows that the VSI is suppressed by magnetic fields in the ideal MHD limit (also see Cui & Lin 2021;Latter & Kunz 2022).In the dead-zone, OR can revive the VSI if the Elässer number is sufficiently low.Latter & Kunz (2022) show that the behaviours of AD are rather complicated because the VSI is coupled with AD shear-instability (Kunz 2008).In our simulations, at least in the regions around the inner dead-zone edge, which is focused on in this paper, the VSI does not appear to grow because the gas is significantly disturbed by the MRI turbulence in the active zone.The radial velocity dispersion is much larger than the vertical one, which is not expected in the VSI growth mode.By contrast, outer regions far from the inner dead-zone edge can suffer from the VSI.Our simulation time however is too short to see the growth the VSI in such outer regions.

Initial Condition
The initial surface density profile of the disk is given by Σini(R) = Σ0 R 1 au −1 (5) where Σ0 = 500 g cm −2 and R is the cylindrical radius.The disk mass integrated from R = 0 to R is given by The mid-plane gas temperature T mid is determined by assuming radiative equilibrium between the gas and radiation from the central star.We consider a central star with M * = 2.5 M⊙ of mass, R * = 2.5 R⊙ of radius, and T * = 10 4 K of effective temperature, which corresponds to a Herbig Ae star.The mid-plane gas temperature is given by where f boost is a coefficient which is set to 2 to increase the aspect ratio of the disk, L * = 4πR 2 * σSBT 4 * is the luminosity of the central star, and σSB is the Stefan-Boltzmann constant.The adopted stellar luminosity is typical in Herbig Ae stars with ages larger than 1 Myr although the luminosities of Herbig Ae stars with ages less than 1 Myr are highly uncertain (e.g., Montesinos et al. 2009).The scale height is given by H(R) = √ 2c s,mid /ΩK(R), where c s,mid = kBT mid /µmH is the sound speed at the midplane, µ = 2.33 is the mean molecular weight, mH is the hydrogen mass, and ΩK = GM * /R 3 is the Kepler angular velocity of the disk.The aspect ratio of the disk is given by The initial mid-plane density is given by We consider a similar temperature distribution as Bai (2017).The whole region is divided into a cold gas disk and warm atmospheres sandwiching the cold disk.Inside the disk, the temperature is given by T d (r) = (µmH/kB)(ϵ mid (r) 2 /2)r 2 ΩK(r) 2 which is identical to T mid (Equation ( 7)) at the mid-plane.In the warm atmosphere well above the disk, the gas is heated up mainly by the radiation from the central star.We define ΣFUV as the column density below which the FUV photons from the central star penetrate.If the radial column density Σr(r, θ) = r r min ρ(r ′ ,θ)dr ′ is smaller than ΣFUV, the gas temperature increases to Tatm(r) = (µmH/kB)(ϵatm(r) 2 /2)r 2 ΩK(r) 2 , where ϵatm(r) = 0.2(r/1 au) 1/4 .The initial temperature profile smoothly connects between T d and Tatm around Σr ∼ ΣFUV as follows: T (r, θ) = (µmH/kB)(ϵ(r) 2 /2)r 2 ΩK(r) 2  (10) and ϵ(r, θ) = ϵ mid (r) + ϵatm(r) 2 + ϵ mid (r) − ϵatm(r) 2 tanh 0.3 log ΣFUV Σr(r, θ) .
The values of ΣFUV are highly uncertain, and Perez-Becker & Chiang (2011) found that FUV layers have a thickness of ∼ 0.01 − 0.1 g cm −2 in plane-parallel models.We adopt ΣFUV = 0.3 g cm −2 .It is three times larger than the upper limit of that of Perez-Becker & Chiang (2011) because the radial column density is considered instead of the vertical column density (Bai 2017).
The aspect ratio for R < 1 au is too small to resolve the MRI structure.As explained below Equation ( 7), we artificially increase the aspect ratio by increasing the gas temperature by a factor of two, or f boost = 2.As mentioned later, for the calculation of ηO and ηA, T f −1 boost is used as the gas temperature.
The initial profiles of the density and ϕ component of the velocity are constructed using an iterative manner to satisfy equations ( 5) and (10) simultaneously in the hydrostatic equilibrium between the pressure-gradient force, gravitational force, and centrifugal force without the Lorentz force.The initial poloidal velocities are set to be zero.The initial magnetic field has only the poloidal components which have an hourglass-shape field given by Zanni et al. (2007).It is described by the ϕ component of the vector potential, where the coefficient C is determined so that the mid-plane plasma beta takes a constant value of β0 = 10 4 .The initial condition of our simulation is displayed in Figure 1.
A height of the base of the atmospheres at a given R is denoted by z = zatm(R) above which the gas temperature is more than twice the mid-plane temperature (Figure 1).

Diffusion Coefficients
The diffusion coefficients ηO and ηA are given by where σO, σH, and σP are the Ohmic, Hall, and Pedersen conductivities, respectively, and c is the speed of light.The three conductivities are given by and where the subscript "j" denotes the species of charged particles, Zje is the charge and nj is the density of the charged particle j.The Hall parameter βj is the ratio between the gyrofrequency of the charged particle j and its collision frequency with the neutral particles, and is defined as where mj is the mass of the charged particle j and γj is the rate coefficient for momentum transfer with the neutral particles through collisions.
The three conductivities are computed by the summation between the conductivities calculated with thermal and non-thermal ionization instead of calculating a chemical network including both thermal and non-thermal ionization, as follows: where i = (O, H, P) and the subscripts "T" and "NT" indicate the thermal and non-thermal ionization limits of the conductivities, respectively.Since the thermal ionization is important in the regions of sufficiently high temperatures, dust grains are assumed to be sublimated in the calculation of chemical reactions.We consider the thermal ionization of potassium, which has an ionization energy of 3.43 eV.The potassium abundance relative to the H2 number density is 3 × 10 −7 .The thermal ionization degree (xe)T is derived by the Saha equation where the electron densities are assumed to be equal to the potassium ion densities.As mentioned before, T is replaced with T f −1 boost in the Saha equation.When the gas temperature exceeds ∼ 1000 K, the thermal ionization of K provides a sufficient high ionization and the region becomes MRI active.
In order to take into account the photo-ionization of carbon and sulfur due to FUV photons in the disk atmosphere, we replaced the ionization degree calculated in the thermal ionization limit by where (xe) T is the ionization degree obtained from the Saha equation for potassium ionization and the second term indicates that the electrons are provided by photoionization of carbon and sulfur when Σr < ΣFUV (Bai 2017).
For the non-thermal ionization, we use the data table of the three conductivities based on a chemical reaction network of e − , H + , He + , C + , H + 3 , HCO + , Mg + in the gas phase and the charged dust grains using the methods described by Nakano et al. (2002) and Okuzumi (2009).The table is a function of ρ/B, T , and ρ/ζ, where ζ is the ionization rate that will be defined in Equation ( 20) and T is replaced with T f −1 boost .For simplicity, the single size dust grain with a radius of a = 0.1 µm and an internal density of 2 g cm −3 is considered.In this paper, the dust-to-gas mass ratio is set to 10 −4 , assuming that dust growth reduces the amount of dust grains smaller than the interstellar value ∼ 0.01.The following three kinds of non-thermal ionization sources are taken into account, The first source is the cosmic rays whose ionization rate is given by where ΣCR = 96 g cm −2 (Umebayashi & Nakano 1981), and is the column density integrated along the θ-direction from θ = 0 (θ = π) to θ at a given r.Secondly, we consider ionization due to the X-ray radiation from the central star using the fitting formula (Bai & Goodman 2009, based on the calculation done by Igea & Glassgold (1999)) as follows: where ζ1 = 6.0 × 10 −11 s −1 , ΣX,a = 3.6 × 10 −3 g cm −2 , α = 0.4, and Σr(θ,r) is the column density integrated along the r-direction from the innermost radius of the simulation box to r at a given θ.The third ionization source is provided by radioactive nuclei, (Finocchi & Gail 1997).To reduce the computational cost, ζ(r, θ) is calculated for the initial condition and fixed during the simulations.

Definition of the Dead Zone Boundaries
The non-ideal MHD effects affect the growth rate of the MRI when the dimensionless Elsässer numbers
In order for the MRI to occur at a given height z, the wavelength of the maximum growing mode measured locally should be smaller than the scale height (Sano & Miyama 1999).Using the expressions of the most unstable wavelengths λmax,O and λmax,A given by Sano & Miyama (1999) and Bai & Stone (2011), respectively, we confirmed that λmax,O > H and λmax,A > H are satisfied in most regions where ΛO < 1 and ΛA < 1, respectively, because of steep spatial gradients of ηO and ηA.This indicates that in the regions just outside the dead zone, the MRI turbulence is partly suppressed as seen in Section 3.2.2.
Hereafter, the dead-zone boundaries for OR and AD are defined as ΛO = 1 and ΛA = 1, respectively.

Spatial Distributions of the Diffusion
Coefficients at the Mid-plane We present the spatial distributions of ηO and ηA in the initial condition.Figure 2a shows the radial profiles of ηO and ηA at the mid-plane.Both resistivities show a similar behavior; they increase steeply as a function of the radius, reach maxima around R ∼ 1 au, and then decrease steeply.
Here we define three characteristic radii from the radial profiles of ηO, ηA, ΛO, and ΛA.The first characteristic radius is RMRI which is the inner edge of the OR dead zone at the mid-plane.As shown in Figure 2b, in the initial condition, the inner dead-zone edges for OR and AD are located in the region where the thermal ionization dominates over the non-thermal ionization.The inner dead-zone edge for OR moves outward until the MRI turbulence is saturated while that for AD do not since ΛO ∝ B 2 and ΛA ∝ B 0 , where we use the fact that ηO ∝ B 0 and ηA ∝ B 2 (Section 2.3.4).In Figure 2b, we plot the radial profiles of ΛO and ΛA at t = 500 tK0 when the MRI turbulence has been saturated.At R = RMRI = 0.5 au, ΛO becomes unity, and RMRI does not change after the MRI turbulence is saturated.Hereafter, the inner dead-zone edge for OR is simply called the inner dead-zone edge.
The second characteristic radius is R = Rη = 1.2 au around which ηO and ηA take the largest values and ΛO and ΛA take the lowest values (Figures 2a and 2b).This corresponds to the radius around which the main ioniza- tion process switches from the thermal to non-thermal ionization.
The third characteristic radius is R ch = 3 au beyond which the main negative charge carrier changes from dust grains to electrons (Section 2.3.4).

Two Dimensional Distributions of the Diffusion Coefficients
The two-dimensional distributions of ηO and ηA in the (R,z) plane at t = 500 tK0 are displayed in Figures 3a and  3b, respectively.ηO decreases from the mid-plane toward upper latitudes because the density decreases.The spatial distribution of ηA is different from that of ηO especially for R < 1 au, where ηA increases from the mid-plane toward upper latitude in the AD dead zone.This comes from the fact that AD is important for strong B and/or low ρ when dust grains do not play an important role in the resistivities.
The shapes of the dead-zone boundaries around the inner edges reflect the physical properties of OR and AD.As the radius increases from the inner boundary of the simulation box, OR makes the disk dead at the mid-plane first around R ∼ RMRI while high latitude regions are the first place which become dead for AD.
The vertical dead-zone boundaries for OR and AD lie between z = H(R) and z = 2H(R), and the AD dead-zone is more extended vertically than the OR dead-zone.The base of the warm atmospheres z = zatm(R) is well above the dead-zone boundaries.

The Dependence of the Diffusion Coefficients on the Field Strength
For OR, ηO is independent of B because B in the factor (ec/B) is cancelled with βj ∝ B in Equation ( 14).By contrast, ηA depends on B in a more complex manner than ηO (Xu & Bai 2016).If there are no dust grains, using the fact that βi ≪ βe, one finds that ηA = βi ×cB/(4πene) ∝ B 2 (Salmeron & Wardle 2003), where βi and βe are the Hall parameters of ions and electrons, respectively, and ne is the electron number density.
The dependence of ηA on the field strength is critical in the structure formation of magnetic fields as will be discussed in Section 5.1.2.Development of sharp structures in the magnetic null due to AD has reported in Brandenburg & Zweibel (1994).For such a structure to develop, ηA should increase with B so that diffusion near the magnetic null is much more inefficient than in strongly magnetized regions.
Existence of dust grains changes ηA significantly.The dependence ηA ∝ B 2 is realized at the weak field limit βe ≪ 1 (Xu & Bai 2016).For βe > ∼ 1, ηA is not necessarily proportional to B 2 if the contribution from dust grains to the number density of charged particles is non-negligible.The dependence of ΛA on the field strength for various grain sizes and dust-to-gas mass ratios are shown in Figure 4 of Xu & Bai (2016).
Figure 4a shows how the radial profile of ηA at the midplane changes as the field strength is increased from the initial value Bini(R).As B/Bini increases, ηA does not increase in proportional to B 2 , but the increase in ηA is saturated especially in the inner region of the dead zone Rη < ∼ R < ∼ R ch , where the dust grains are the main negative charge carrier.Figure 4b shows the field strengths where d ln ηA/d ln B = 1 and 0.5.For field strengths larger than these critical field strength, the B dependence of ηA become weak, leading to inefficient formation of substructures due to AD.

Methods
To solve the basic equations (1)-(3), we use Athena++ (Stone et al. 2020) which is a complete rewrite of the Athena MHD code (Stone et al. 2008).The HLLD (Harten-Lax-van Leer-Discontinuities) method is used as the MHD Riemann solver (Miyoshi & Kusano 2005).Magnetic fields are integrated by the constrained transport method (Evans & Hawley 1988;Gardiner et al. 2008).The second-order Runge-Kutta-Legendre super-time stepping technique is used to calculate the magnetic diffusion processes (Meyers et al. 2014), where the time step is limited so that it is not greater than 30 times the time step determined by OR and AD.
A spherical-polar coordinate system (r, θ, ϕ) is adopted in the simulation box rin ≤ r ≤ rout, 0 ≤ θ ≤ π, and 0 ≤ ϕ < 2π, where rin = 0.1 au and rout = 12.5 au.The radial cell width increases with radius by a constant factor so that the radial cell width is almost identical to the zenith one.
The static mesh refinement technique is used to resolve the disk region which requires high resolution.We adopt two mesh configurations as shown in Figure 5 angle enclosed by the black lines in Figure 5 consists of (Nr, N θ ) = (24, 16).Table 1 lists the models considered in this paper, where tK0 = 2π r 3 in /GM * is the Kepler rotation period at the inner boundary of the computation box.LowRes run corresponds to the lower resolution simulation where the meshes are refined by 5 levels.If the whole computation domain was divided by the finest cell, the resolution would be (Nr, N θ , Nz) = (1536, 1024, 1024).The disk scale height H(R) is resolved by 19 (R/RMRI) 1/4 cells.
LowRes run was conducted until 2000 tK0 (Table 1).As will be shown in Section 1, the resolution of LowRes run is not high enough to obtain the converged MRI turbulence.For comparison, we conduct HighRes run, where the cells in the active zone are further refined in one more higher level (the lower panel of Figure 5).HighRes run could only be calculated to a short time, 500 rotations at the inner boundary, which is long enough for the MRI turbulence to be a saturated state in the active zone.
We note that the resolution of LowRes run is high enough to drive turbulence in the active zone although the saturated ⟨α Rϕ ⟩H is underestimated by a factor of 2 compared with HighRes run that is expected to give the converged ⟨α Rϕ ⟩H (Appendix 1).In this paper, the long-term evolution is shown on the basis of LowRes run, referring the results of HighRes run to keep the resolution issue in mind.
To prevent the time step ∆t from being extremely small due to the Alfvén speed, we set the following spatially dependent density floor, ρ(r, θ, ϕ) = max ρ(r, θ, ϕ), 10 −6 ρ(R, z = 0) . (25) This density floor works only in the polar regions.In the other regions, the disk winds keep the densities greater than the floor value.

Boundary Conditions and Buffer Layer
In Athena++, the boundary conditions are applied by setting the primitive variables in the so-called ghost zones located outside of the computation domain.We impose the following inner boundary conditions.The density and pressure in the ghost cells are fixed to be the initial values, the toroidal velocity follows the rigid rotation v ϕ = min(RΩK(rin), v ϕ,ini ), and vr and v θ are set to zero, where v ϕ,ini is the initial toroidal velocity.We set the boundary conditions for Br and B ϕ so that Br ∝ r −2 and B ϕ ∝ r −1 .
The outer boundary conditions are as follows.For the density, pressure, v θ and B θ , zero-gradient boundary conditions are imposed.We set the boundary conditions for v ϕ , Br, and B ϕ so that v ϕ ∝ r −1/2 , Br ∝ r −2 , and B ϕ ∝ r −1 .
The radial velocity is set by using the zero-gradient boundary condition only for vr > 0 otherwise set to zero.
In order to mitigate the artificial effect of the inner boundary conditions, we introduce a buffer region in rin ≤ r ≤ r buf , where r buf − rin = 0.025 au is the width of the buffer region.
The treatment in the buffer region is similar to that of Takasao et al. (2018).The poloidal velocities are damped according to the following equation: where τ d is the damping timescale, and f d,min = 1, f d,max = 100, and vr,c = −0.05vK,in.The other physical quantities are unchanged in the damping layer.

Normalization Units
In this paper, all the physical quantities are normalized by the following three quantities, the length scale R0 = 1 au, the time scale t0 = R 3 0 /GM * = 0.1 yr, and the surface density Σ0 = 500 g cm −2 .The gas temperature is normalized by T0 = (µmH/kB)(R0/t0) 2 = 6.3 × 10 5 K. Hereafter, we use normalized physical quantities unless otherwise noted except for R and z that are shown in unit of the astronomical unit.

Averaged Quantities
We define averaged quantities used to analyse the simulation results in this section.In the subsequent sections, we present the data converted from the spherical polar coordinates (r, θ, ϕ) to the cylindrical coordinates (R, ϕ, z).
We define the data averaged over ϕ, which are denoted by using angle brackets, or ⟨A⟩ ϕ .Different weight functions are used for different physical quantities in taking averages.For quantities with the dimensions of mass density, momentum density, energy density, and field strength, we take the simple volume-weighted average where a quantity A can be ρ, ρvz, B 2 ϕ /8π, or BR.For velocities and kinetic energies per unit mass, we take the density-weighted average where a quantity A can be vz or v 2 R /2.In order to eliminate stochastic features originating from the MRI turbulence, ⟨A⟩ ϕ is averaged over t1 ≤ t ≤ t2 as follows: When the radial dependence of the disk structure is investigated, we take the following vertical average of the disk, where z b is a thickness over which ⟨A⟩ is averaged.
When the ϕ dependence of the physical quantities is investigated, we take the following vertical average over Using ⟨v⟩ and ⟨B⟩, we define turbulent components of velocities and magnetic fields as variations from the ϕaveraged values as follows; where the turbulent components (δv and δB) satisfy ⟨ρδv⟩ = 0 and ⟨δB⟩ = 0.

Main Results
Figure 6 summarises our findings about the gas dynamics and magnetic field properties in our disks.Details are provided in the subsequent sections.

Overall Features
In this section, we briefly describe overall features of our results.Figures 7 and 8 show the spatial distributions of the four ϕ-averaged variables in the inner region of R ≤ 2 au and outer region of R ≤ 6 au at t = 2000tK0 for LowRes run, respectively.We note that the outer region has not reached a quasi-steady state since 2000 rotations at the inner boundary corresponds to only 12 rotations at R = 3 au.Figure 9 shows the face-on color maps of the four vertically-averaged variables.From Figure 7, in the polar regions of the northern and southern hemispheres, we identify the so-called funnel magnetic fields that are connected with the inner boundary and the magnetic energy dominates over the thermal    energy.The radial size of this region increases with time owing to magnetic flux accretion (Beckwith et al. 2009;Takasao et al. 2019).The size of the funnel regions may be overestimated because in reality the funnel magnetic fields come from open fields around the central star, whose size is much smaller than the inner boundary in our simulations.We found that the conventional dead zone identified by ΛO ≤ 1 and ΛA ≤ 1 can be divided into two regions separated by R = Rη.We call the inner region of RMRI ≤ R ≤ Rη "the transition zone" which has different properties from the conventional dead zone.We call the outer region of R > Rη "the coherent zone" which has the same properties as the conventional dead zone.
The overall properties of the active, transition, and coherent zones are briefly summarized in Sections 3.1.1,3.1.2,and 3.1.3,respectively.In Section 3.1.4,we discuss the influence of the active zone on the outer regions.

Active Zone (R ≤ RMRI)
The physical properties of the MRI turbulence in the active zone are consistent with those found in the previous studies.We briefly summarise the physical properties of the MRI turbulence, and detailed descriptions are found in Appendix 2.
The vertical structure of the magnetic fields changes around |z| ∼ 2H.The MRI turbulence generates turbulent magnetic fields for |z| < 2H while the magnetic fields become coherent in the upper atmosphere (Suzuki & Inutsuka 2009).We observe a so-called butterfly structure in the t-z diagram of B ϕ at a fixed radius (Figure 37).The signs of the toroidal field change quasi-periodically and B ϕ drifts toward upper atmospheres.
The MRI turbulence drives gas accretion in the inner region of the active zone.The outer region expands outward by receiving the angular momentum from the inner region (Lynden-Bell & Pringle 1974).In the upper atmospheres (|z| > 2H), the magnetic torque of the coherent magnetic field drives coherent gas motion near the surface layers.
One of the striking features of our simulations is the occurrence of multiple ring structures in the density field in the long-term evolution, as shown in Figures 7a and 9a (also see Figure 11).The spatial distributions of the density and Bz are anti-correlated; the magnetic fluxes are concentrated in the density gaps.Possible formation mechanisms of these ring structures are investigated in Section 3.2.1 and are discussed in Section 5.1.1.

Transition Zone (RMRI < R < Rη)
We discover a distinct region, the transition zone, in this simulation.Although it was traditionally classified as a dead zone since ΛO < 1 and ΛA < 1 are satisfied in most regions, it has many characteristics not found in conventional dead zones.An interesting feature is found in the magnetic field structures in Figure 7d.The vertical magnetic fields almost completely disappear in RMRI < ∼ R < ∼ Rη (Section 3.2.3).Such a region has not been found in the previous studies (Lyra & Mac Low 2012;Dzyurkevich et al. 2010;Flock et al. 2017) because AD, which was neglected in their work, plays an essential role (Section 3.2.3).The gap in the vertical magnetic fields is almost concentric as shown in Figure 9d.
The disappearance of the vertical magnetic field suppresses surface gas accretion expected in the conventional dead zone (Section 3.2.2). Figure 7b shows disturbances in the radial mass flux although there are no turbulence driving mechanisms.These disturbances originate from the active zone (Pucci et al. 2021), and the net radial mass flux is extremely low when ⟨ρvR⟩ is averaged over time (Section 3.2.2).Thus, in the transition zone, gas accretion does not occur either around the mid-plane or on the surface of the disk.
Figure 9a clearly shows that a density peak appears at each edge of the transition zone, or R ∼ RMRI and R ∼ Rη.This structure is formed by a combination of no net radial gas motion in the transition zone and the gas supply from the active zone (R ≤ RMRI) and the coherent zone (R ≥ Rη).Details will be investigated in Sections 4.1 and 4.2.2.

Coherent Zone (R ≥ Rη)
The coherent zone with R ≥ Rη has properties consistent with those of the conventional dead zone found in the literature.The magnetic field lines are smooth and coherent both inside and outside the disk.
Just above the coherent zone, the toroidal magnetic field is amplified by the differential rotation because the magnetic field is relatively coupled with the gas.The magnetic tension force −Bz∂B ϕ /∂z extracts the angular momentum from the gas efficiently, driving surface gas accretion between the OR and AD dead-zone boundaries as shown in Figure 7b (Bai & Stone 2013;Gressel et al. 2015).The electric current sheet where the sign of ⟨B ϕ ⟩ is reversed is located not at the mid-plane but around the lower AD dead-zone boundary, indicating that the zsymmetry adopted in the initial condition is broken (Figure 7c).This is because inside the disk, ηO and ηA are so large that a current sheet cannot exist inside the coherent zone (Figure 2a).As a result, a current sheet is lifted either upward and downward to the height where ΛO is larger than unity (Bai & Stone 2013;Bai 2017).The off-mid-plane current sheet causes the surface gas accretion to be asymmetric with respect to the mid-plane because the magnetic torque exerted in the lower side is stronger than that in the upper side.
The behaviors of the magnetic field change around R ∼ R ch .For R < ∼ R ch , the current sheet is located at the lower disk surface as explained before.Beyond R ∼ R ch , ηO and ηA are small enough for the current sheet to remain at the mid-plane.A similar feature was reported in Lesur (2021) who demonstrated that the surface gas accretion (mid-plane gas accretion) occurs for lower (higher) ΛO. Figure 8c shows that the toroidal field around the midplane at R > R ch is amplified, resulting in the OR deadzone shrinking vertically toward the mid-plane.Around the mid-plane, the toroidal magnetic fields are amplified by AD.A similar amplification has been found in Suriano et al. (2018Suriano et al. ( , 2019) ) and also around the inner edge of the transition zone as discussed in Section 3.2.3.The magnetic tension force −Bz∂B ϕ /∂z drives gas accretion at the midplane in R > R ch (Figure 8b).

Influence of the Active Zone on the Transition and Coherent Zones
The MRI activity in the active zone affects both the velocity and magnetic fields in the outer regions.As mentioned in Section 3.1.2,spatial variations in the radial mass flux inside the disk seen in Figures 7b and 8b are caused by outward propagation of disturbances generated by MRI turbulence (Section 3.2.5).
The MRI activity induces quasi-periodic inversion of the sign of the toroidal magnetic field in the transition zone and the inner part of the coherent zone (Section 3.2.4).This leads to a quasi-periodic switching of the current sheet position between the top and bottom dead-zone boundaries.At t = 2000tK0, the toroidal field in the disk is negative (Figures 7c and 8c), but at another time it can be positive.The quasi-periodic disturbance of the field structures does not penetrate beyond R ∼ R ch .Comparison between Figures 8c and 8d shows that ⟨Bz⟩ has a concentration around R ∼ R ch .Concentrations in ⟨Bz⟩ propagate outward following quasi-periodic inversion of the sign of the toroidal magnetic field.Details will be investigated in Section 3.2.4.

Main Findings
In this section, we present detailed analyses about main findings in our simulations.

Ring Formation in the Active Zone
We observe the formation of multiple rings and gaps in the active zone as shown in the radial profiles of the column densities Σ (Figure 10a).Comparison between Figures 10a and 10b shows that the ring structures in Σ are anticorrelated with the radial profile of the vertical magnetic flux averaged over the scale height ⟨Bz⟩H ; the rings (gaps) of Σ correspond to the gaps (rings) of ⟨Bz⟩H .
The top and middle panels of Figure 11 show the face-on views of the density and Bz averaged over |z| < H(R).At t = 1000tK0, multiple-rings are found in the density map.
Thin magnetic flux concentrations are found around R ∼ 0.2 au although the magnetic flux concentrations are not as significant as the variations of the density map.At t = 2000tK0, the density contrasts between the rings and gaps become significant, and the magnetic flux concentrations at the density gaps are more prominent than those at t = 1000tK0.
The resolution dependence on the density distributions are shown in the bottom panels of Figure 11.Although the density is lower for HighRes run than for LowRes run, ring structures are also found in HighRes run.This suggests that if the long term simulation performed for HighRes run, the ring structures would develop.
Several formation mechanisms of ring structures in the ideal MHD has been proposed.In this section, we examine two mechanisms: the viscous instability (Lightman & Eardley 1974;Suzuki 2023) and the wind-driven ring formation (Riols & Lesur 2019).As shown below, the viscous instability is a possible mechanism of the ring formation found in our simulations while the wind-driven instability may not occur.In Section 5.1.1,we further discuss the origin of the ring structures in our simulations.
Although a quantitative argument why pρ becomes greater than unity is still missing, an anti-correlation between ⟨ρ⟩H and ⟨Bz⟩H seen in Figures 11a and 11b  .When ⟨Bz⟩H and ⟨ρ⟩H are anticorrelated (qB z > 0), the density dependence of α ana,β is apparently stronger than when α ana,β is assumed to be a function of ⟨βz⟩H .However this simple argument does not quantitatively explain our results.LowRes run shows that qB z ∼ 0.43 for t = 1000 tK0 and ∼ 0.95 for t = 2000 tK0.One obtains that −∂ lnα ana,β /∂ ln⟨ρ⟩H = p β (2qB z +1) ∼ 0.76 for t = 1000 tK0 and ∼ 0.93 for t = 2000 tK0, and both values are not consistent with pρ shown in Table 2.

Wind-driven Instability
Riols & Lesur (2019) proposed a wind-driven instability where disk winds destabilize disks when the amount of gas removed from gap regions due to winds is greater influence of the buffer layer (Section 2.5) and the transition zone.The least square methods are applied to the scattered data in the (ln⟨βz⟩ H , ln⟨α Rϕ ⟩ H , ) and (ln⟨ρ⟩ H , ln⟨α Rϕ ⟩ H ) planes. 2 The power-law index p β may be lower than those obtained in local shearing-box simulations that span from p β ∼ 0.5 to p β ∼ 1. Suzuki et al. (2010) and Okuzumi & Hirose (2011) found that p β ∼ 1, and the results of Sano et al. (2004)  than that supplied by viscous diffusion from ring to gap regions.Local simulations found that both the α parameter and mass loss rate depend negatively on the plasma beta; more strongly magnetized disks yield faster viscous diffusion and more efficient mass loss.Thus, in order for the wind-driven instability to occur, the mass loss rate should has a sensitive dependence on the plasma beta than the viscous diffusion rate.
We here define the normalized mass loss flux due to the disk wind averaged over time as where vn stands for the velocity component perpendicular to the ±zatm surfaces (Figures 1 and 3), ρ mid is the midplane density, and c s,mid is the mid-plane sound speed.
Assuming that both Cw and α are anti-correlated with βz and they follow the relations Cw ∝ β −pw z and α ∝ β Lesur (2019) found that the wind-driven instability occurs when pw > p β .Since p β is roughly equal to 0.3 − 0.4 in our simulations (Figure 10c and Table 2), pw needs to be larger than 0.3 − 0.4 in order for the winddriven instability to occur.
Figure 12a compares the radial profiles of Cw with the inverse of the time-averaged mid-plane plasma beta ⟨β z,mid ⟩t ≡ 8π⟨P ⟩t(z = 0)/⟨Bz⟩ 2 t (z = 0).Cw is poorly correlated with ⟨β z,mid ⟩ −1 , indicating that the wind-driven instability is not caused by the disk wind at least in our simulations.
Figure 12a shows that in the gap regions where the radial profiles of ⟨β z,mid ⟩ −1 t have local maxima, Cw is negative, indicating that the gas flows into the density gap regions rather than being ejected from the disk in the vertical direction.Why are there no outflows from the gap regions where the magnetic field is strong?Riols & Lesur (2019) pointed out the gas in the gap regions is ejected by "wind plumes" where the magnetic fields are strong enough to be coherent and are tilted with respect to the zaxis (see also Riols & Lesur 2018).Figure 12b shows that the poloidal magnetic field structures originating from the gap regions do not maintain a large tilt with respect to the z-axis as |z| increase, suggesting that the gas is not continuously accelerated along the field lines.

Disk Wind Launched Around the Transition-
zone Inner Edge and Non-existence of Gas Accretion in the Transition Zone The transition zone is disturbed by the influence from the turbulence in the active zone, and exhibits time variations as will be discussed in Sections 3.2.4 and 3.2.5.In this section, we investigate the quasi-steady structure of the transition zone by taking time average.The left panels of Figure 13 show the close-up views of the magnetic fields and velocities around the dead-zone inner edge.The MRI turbulence appears to be suppressed around R ∼ 0.45 au, which is slightly different from R = RMRI = 0.5 au defined by ΛO = 1.This is because MRI turbulence is partially suppressed even where ΛO is slightly greater than unity.
Around the mid-plane in 0.45 au < ∼ R < ∼ 0.55 au, hourglass-shaped poloidal magnetic fields develop.The vertical distribution of each component of the magnetic field at R = RMRI is shown in Figure 14a.The toroidal magnetic field dominates over the other components, and has a sharp gradient at the mid-plane.
AD plays a critical role in the formation of the hourglass-shaped poloidal magnetic fields in 0.45 au < ∼ R < ∼ 0.55 au.To investigate which of the ideal, OR, and AD terms is the most effective to produce the hourglassshaped magnetic fields, we measure their contributions to   −⟨(∇ × E) R ⟩, which is equal to ∂⟨BR⟩/∂t, at a given z.
First, we focus on the region of |z| < ∼ 0.01 au where the mid-plane gas accretion is seen in Figure 13c.Figure 14b shows that the AD term tilts the poloidal magnetic field toward the radial direction of the disk (⟨(∂BR/∂t)A⟩ > 0 for z > 0, and ⟨(∂BR/∂t)A⟩ < 0 for z < 0), and the ideal MHD term does the same.The vertical profile of ⟨BR⟩ is almost stationary by the balance between diffusion due to the OR term and amplification due to the AD and ideal MHD terms.Far from the mid-plane (|z| > ∼ 0.02 au), the AD term behaves quite the opposite.The AD term works as diffusion of ⟨BR⟩ amplified by the ideal MHD term.
AD also amplifies the toroidal fields near the mid-plane in 0.45 au < ∼ R < ∼ 0.55 au. Figure 14c is the same as Figure 14b but for ⟨B ϕ ⟩. Figure 14c shows that the signs of ⟨B ϕ ⟩ and ⟨(∂B ϕ /∂t)A⟩ are the same near the mid-plane (|z| < ∼ 0.01 au), indicating that the AD term amplifies ⟨B ϕ ⟩ and steepens its gradient.The vertical profile of ⟨B ϕ ⟩ is kept almost stationary by the OR term smoothing ⟨B ϕ ⟩.The ideal MHD term partially contributes to the steepening of ⟨B ϕ ⟩.In a similar way as ⟨(∂BR/∂t)A⟩, far from the midplane, the AD term works as diffusion of ⟨B ϕ ⟩ amplified by the ideal MHD term.
The magnetic torque due to −Bz∂B ϕ /∂z drives midplane gas accretion (Figures 13a and 13c).Just above the thin mid-plane gas accretion layer, the wind-like gas flows directing outward are driven.This is also caused by the magnetic torque of the coherent magnetic fields (Blandford & Payne 1982;Bai et al. 2016).The gas streaming lines at lower latitudes return to the mid-plane.As a result, meridional flows are formed in 0.45 au < ∼ R < ∼ 0.5 au and |z| < ∼ 0.05 au; the streamlines of the gas flows are circulated each in the north and south sides of the disk as shown in Figure 13c.
The right panels of Figure 13 zooms out the region shown in the left panels to cover the entire transition zone.Four streamlines are highlighted by colors with thick lines in Figure 13f.One can find that the lower two streamlines correspond to failed disk winds; the material does not reach the outer boundary of the simulation domain but falls back on to the disk surface in R < ∼ Rη by the central star gravity (Takasao et al. 2018).By contrast, the disk wind flowing from the higher latitudes (z = 0.073 au and 0.08 au, green and red) reaches the outer boundary of the simulations box.The directions of the winds are not parallel to the poloidal magnetic fields owing to AD (Figure 14b).

Disappearance of the Vertical Magnetic Flux
in the Transition Zone In this section, we investigate how the vertical flux disappears in the transition zone (Figures 7d and 9d).This flux transport is very efficient, and occurs in less time than 7 rotations at R = RMRI.
Figure 15 shows the time sequence of the magnetic field structure.At t = 50tK0, the poloidal magnetic field lines are almost vertical inside the AD dead zone at R > 0.6 au.After 1.8 rotations at R = RMRI (t = 70tK0), the poloidal magnetic fields are highly inclined and their directions are almost horizontal in 0.6 au < ∼ R < ∼ 0.8 au.Inside the transition zone, the radial magnetic fields are amplified.At t = 90 tK0, the magnetic loop structure whose center is at R = 0.62 au and z = 0 is formed.The time evolution shown in Figure 15 occurs on a dynamical timescale at R = 0.6 au.
In order to investigate the evolution of the magnetic fields, the vertical profiles of the ϕ-averaged magnetic fields What causes this magnetic field evolution?The answer can be obtained from the ϕ-averaged induction equations, which are given by and The electric field can be divided into three components, where EI = −v × B is the electric field in the ideal MHD, EO = ηOJ and EA = ηOJ ⊥ are the electric fields caused by OR and AD, respectively.The evolution of the vertical field is determined by Equation ( 40).In the early evolution, ⟨(E ϕ )O⟩ provides a dominant contribution to ⟨E ϕ ⟩ as will be shown in Figure 18.Because |∂BR/∂z| dominates over |∂Bz/∂R| the vertical distribution of BR is critical for the vertical field transport.
The evolution of ⟨BR⟩ is determined by the vertical structure of ⟨E ϕ ⟩ (Equation ( 38)).Around the mid-plane, ⟨(E ϕ ) I ⟩ is negligible since this region is inside the dead zones.For t ≤ 60 tK0, ⟨(E ϕ ) A ⟩ mainly contributes to ⟨E ϕ ⟩.Around the mid-plane, ∂⟨(E ϕ ) A ⟩/∂z is positive (negative) for z > 0 (< 0), leading to amplification of ⟨BR⟩ since the signs of ∂⟨(E ϕ ) A ⟩/∂z and ⟨BR⟩ are the same (the left panel of Figure 16).This downward-facing convex profile of ⟨(E ϕ ) A ⟩ is attributed to the profile of ⟨ηA⟩, which increases toward upper low density regions (Figure 16f).Since ⟨ηA⟩ is proportional to B 2 in this region, the gradient of ⟨ηA⟩ becomes steeper and steeper owing to the amplification of the magnetic field.As both ⟨BR⟩ and ⟨B ϕ ⟩ increase (Figure 16), the current density is enhanced around the mid-plane, and the electric field owing to OR becomes important.For t ≥ 70tK0, the downward-facing convex profiles disappear in the ⟨E ϕ ⟩ profile since the downward-facing convex part in ⟨(E ϕ ) A ⟩ is almost compensated by the upward-facing convex part in ⟨(E ϕ ) O ⟩.As a result, OR suppresses the amplification of ⟨BR⟩.
The gradient of ⟨BR⟩ with respect to z become steeper as the time passes, leading that the current sheet around the mid-plane becomes thinner around the central region where the magnetic fields are weak.Development of such sharp structures in the magnetic null by AD was previously reported in Brandenburg & Zweibel (1994).
Finally, we investigate the radial transport of the vertical field.Equation (40) shows that the transport velocity of ⟨Bz⟩ is estimated by ⟨E ϕ ⟩/⟨Bz⟩.Since ⟨Bz⟩ is positive in most regions at the mid-plane, the sign of ⟨E ϕ ⟩ indi- cates the direction of the vertical field transport.Figure 18 shows the time evolution of the radial profiles of the toroidal electric fields at the mid-plane.OR mainly contributes to ⟨E ϕ ⟩ while ⟨E ϕ ⟩I is negligible at the mid-plane where OR suppresses the MRI.AD plays a minor role in ⟨E ϕ ⟩ except for R > ∼ 1 au.A sharp increase in ⟨E ϕ ⟩ in R < ∼ 0.9 au is attributed to a strong radial dependence of ⟨ηO⟩ in R < Rη (Figure 2).The radial gradient of R⟨E ϕ ⟩ around R ∼ 0.8 au increases with time from t ∼ 50 tK0 to 70 tK0, representing the rapid outward transport of ⟨Bz⟩ (Figure 16).The rapid outward transport of the vertical field in the transition zone is due to steepening of the vertical gradient of BR around the mid-plane by AD as seen in Figures 16 and 17.

Quasi-periodic Disturbances of Magnetic
Fields Beyond the Inner Dead-zone Edge due to the MRI Activity The turbulence and dynamo activity in the active zone drives time variation in the transition and coherent zones.
Figure 19 shows the radius-time diagrams of ⟨Bz⟩H and ⟨B ϕ ⟩H .As shown in Section 3.2.1, in the active zone (R ≤ RMRI), the ring and gap structures develop in ⟨Bz⟩H and ⟨ρ⟩H .
In the transition zone (RMRI < R < Rη), as shown in Section 3.2.3, the combination of OR and AD causes the magnetic flux to be transferred outward rapidly in the early evolution ∼ 100tK0.As a result, ⟨Bz⟩H is almost zero in the transition zone.The gap structure is maintained at least during the simulation.
Before showing the evolution of ⟨B ϕ ⟩H , we recall the spatial structure of the magnetic fields found in Figures 7c and 13d.Just outside the inner edge of the transition zone, B ϕ is amplified above and below the mid-plane, and the current sheet where the sign of B ϕ is flipped is located around the mid-plane.As shown in Figure 14, the profile of B ϕ is determined by the balance between dissipation of B ϕ due to OR and B ϕ amplification due to AD and the induction term.For R > ∼ 0.7 au, strong magnetic diffusion moves the current sheet to either upper or lower AD deadzone boundaries (Figure 7c).
Figure 19b shows that ⟨B ϕ ⟩H just outside the inner edge of the transition zone changes their sign quasi-periodically.Since the position where ⟨B ϕ ⟩ = 0 is around the mid-plane, the quasi-periodic variations of ⟨B ϕ ⟩H are not caused by change of the current sheet position.The vertical profile of ⟨B ϕ ⟩ is not perfectly inversely symmetric with respect to the mid-plane, and the difference between |⟨B ϕ ⟩| for z > 0 and for z < 0 varies quasi-periodically.
Figure 19b shows that the variations of ⟨B ϕ ⟩H propagate outward rapidly in the radial direction while they do not propagate beyond R ∼ R ch .These In order to investigate the origin of the quasi-periodic time variation in ⟨B ϕ ⟩H in RMRI < ∼ R < ∼ R ch , we compare the time variations of B ϕ at three different radii, R = 0.3 au, 0.7 au, and 1.5 au.In the active zone, net toroidal field ⟨B ϕ ⟩ is almost zero for |z| ≤ 2H because of the MRI turbulence, and the quasi-periodic variations of ⟨B ϕ ⟩ are seen in higher latitudes as the butterfly structure in the time-z diagram (see Figure 37a).In order to quantify the time variations of ⟨B ϕ ⟩, we define Since the signs of ⟨B⟩ ϕ in z > 0 and z < 0 are opposite, ∆⟨B ϕ ⟩ is a measure of the antisymmetry of ⟨B ϕ ⟩ with respect to the mid-plane.
Figure 20 shows that the time variations of ∆⟨B ϕ ⟩ at R = 0.3 au are correlated with those of ⟨B ϕ ⟩H at R = 0.7 au.This means that the drift of the magnetic fields in the active zone disturbs the transition zone.The variations of ⟨B ϕ ⟩H in the transition zone propagate rapidly in the radial direction owing to efficient magnetic diffusion.This behavior can be seen in Figure 20  between ⟨B ϕ ⟩H at R = 1.5 au and ⟨B ϕ ⟩H at R = 0.7 au.
Next, we investigate the radial diffusion of the magnetic fields in the coherent zone (R > Rη).We focus on the second flip of ⟨B ϕ ⟩H starting around t ∼ 500 tK0 (Figure 19b).
Figure 21a shows the time evolution of the radial pro-files of ⟨BR⟩, ⟨B ϕ ⟩, and ⟨Bz⟩ at the mid-plane.When the time passes from t = 540 tK0 to 600 tK0, the radius where the sign of ⟨B ϕ ⟩ is reversed moves outward while ⟨BR⟩ and ⟨Bz⟩ do not change significantly.At t = 600 tK0, a concentration of ⟨Bz⟩ appears around R ∼ 2.1 au where the sign of ⟨B ϕ ⟩ is flipped.
The radial diffusion speed of ⟨B ϕ ⟩ is roughly estimated by ⟨ηO⟩ and ⟨ηA⟩, which are shown in Figure 21b.⟨ηA⟩ increases by increasing the field strength while ⟨ηO⟩ does not change in time significantly.Around R ∼ 1.2 au, ⟨ηO⟩ and ⟨ηA⟩ are both roughly ∼ 0.3.The spatial scale of the spatial variation of ⟨B ϕ ⟩ is around ∆R ∼ 0.5 au.The typical speed of the radial diffusion of the magnetic fields is estimated to be ⟨η⟩/∆R ∼ 0.1 au t −1 K0 , indicating that the structure of ⟨B ϕ ⟩ moves 1 au per 10 tK0.Although this speed is a few times larger than estimated from Figure 21a, they are consistent based on the rough estimate.
The ⟨Bz⟩ concentration is determined by R⟨E ϕ ⟩. Figure 21c shows the radial profiles of R⟨E ϕ ⟩ and the contributions of OR and AD to R⟨E ϕ ⟩ which are denoted by R⟨(E ϕ )O⟩ and R⟨(E ϕ )A⟩, respectively.The contribution from the ideal-MHD (induction) term −v × B is not plotted in Figure 21c because it is negligible.Figure 21c shows that R⟨(E ϕ )O⟩ and R⟨(E ϕ )A⟩ have opposite signs and nearly cancel each other.This indicates that AD behaves anti-diffusion since OR provides pure diffusion of magnetic fields.As shown in Figure 21d, the opposite sign of R⟨(E ϕ )O⟩ and R⟨(E ϕ )A⟩ comes from the fact that the sign of ⟨J ϕ ⟩ is opposite to ⟨(J ϕ ) ⊥ ⟩ since ⟨(E ϕ )O⟩ ∝ ⟨J ϕ ⟩. and ⟨(E ϕ )A⟩ ∝ ⟨J ϕ ⟩ ⊥ .This feature was previously pointed out by Béthune et al. (2017).Anti-diffusion owing to AD is not completely cancelled out by the normal diffusion owing to OR, and triggers the concentration of ⟨Bz⟩.
Comparison between Figures 19a and 19b shows that the concentrations of ⟨Bz⟩H are located at the radii where the sign of ⟨B ϕ ⟩H is flipped for R > ∼ Rη.The outward propagation of the concentrations of ⟨Bz⟩H slows down as the radius increases, and almost stall around R ∼ R ch .This is because ηA decreases rapidly beyond R = R ch , resulting in a rapid decrease of the propagation speed of the condensations of ⟨Bz⟩H .

Radial
Propagation of the Velocity Disturbances Originating from the MRI Turbulence The MRI turbulence driven in the active zone generates velocity disturbances that propagate outward and disturb the transition and coherent zones as seen in Figure 7b and  8b.
Figure 22a shows the Mach numbers of the velocity dispersion in the R, ϕ, and z directions at the mid-plane.
They are defined as Inside the active zone, ⟨M R,ϕ,z ⟩ are as large as ∼ 0.1.Comparison between the Mach numbers and density in Figure 22a shows that the velocity dispersion is anticorrelated with the density; it is larger in the gaps than in the rings.
In the coherent zone (R ≥ Rη), ⟨δMR⟩ and ⟨δMz⟩ are of similar magnitude and decrease with increasing radius.⟨δM ϕ ⟩ is smaller than the other two.
We investigate the radial propagation of the radial velocity fluctuations.The energy flux of the sound waves is estimated as where δρ ≡ ρ − ⟨ρ⟩.It can be divided into the energy flux propagating in the +R direction and that propagating in the −R direction They satisfy ⟨FR⟩ = ⟨(FR)+⟩ − ⟨(FR)−⟩.
Outside the active zone, ⟨(FR)+⟩ dominates over ⟨(FR)−⟩.This clearly shows that the sound waves transport more energy outward.Around R ∼ Rη, ⟨(FR)−⟩ is enhanced and becomes comparable to ⟨(FR)+⟩ although the outward flux is still larger than the inward flux.This indicates that a part of the outgoing waves are reflected by the density bump around R = Rη (Figure 22a).

Mass and Angular Momentum Transfer Throughout the Disk
In this section, we investigate the radial dependence of the mass accretion rate throughout the disk.The mass accretion rate at a given R is defined as where z b is a reference height, and it will take the value of either H or zatm.
We derive some equations to analyze the angular momentum transfer.We start from the continuity equation and the momentum conservation equation in the ϕ direction where they are averaged over 0 ≤ ϕ < 2π.
The Rϕ and ϕz components of the Reynolds stress tensor are decomposed into the laminar and turbulent parts, where i = (R, z) (Equation ( 35)).Combining the Maxwell stress tensor and the turbulent Reynolds stress tensor, we define the Rϕ and ϕz stresses as and respectively.Substituting Equations ( 50)-( 52) into Equations ( 48) and ( 49) and averaging Equations ( 48) and ( 49) over t1 ≤ t ≤ t2, one obtains a prediction of the radial mass flux Ṁ z b pred , which is given by integrating 2πR⟨ρvR⟩ over −z b ≤ z ≤ z b as a function of radius as follows: where where the terms with ∂/∂t are neglected.In the derivation of the above equations, ⟨v ϕ ⟩ = RΩK is assumed.We note that this approximation is not very accurate in the ring and gap regions where the rotation velocity deviates from the Keplerian profile due to the gas pressure gradient.The first and second terms correspond to the radial mass fluxes driven by the torques due to ⟨W Rϕ ⟩ and ⟨W ϕz ⟩, respectively.These two provide the dominant contribution to Ṁ z b pred .The third term comes from the radial mass flux affected by the vertical mass flux.

Early Evolution
First, we investigate the angular momentum transfer in the early evolution 400tK0 ≤ t ≤ 500tK0.The simulation time is long enough for the turbulence in the active zone to reach saturation (Figure 36), but not long enough to capture the secular evolution, such as the ring formation and flux concentration (Figure 11).At R = Rη, the gas rotates only 14 rotations at t = 450tK0, indicating that both the transition and coherent zones have not reached quasisteady states.In this section, we investigate the angular momentum transport mechanism especially in the active zone and around the inner edge of the transition zone.
The top and second rows of Figure 23 show comparison between the measured and predicted mass accretion rates by integrating over |z| ≤ H and |z| ≤ zatm, respectively.Equation ( 53) reproduces the actual mass accretion rates reasonably well for all the panels.
Next, we investigate the mass transfer in the active zone around the mid-plane (|z| ≤ H).Figures 23a and 23b show that the mass transfer is mainly driven by the Rϕ torque due to the MRI turbulence while the contribution from the ϕz torque is negligible.The active zone shows that the radial mass fluxes are negative at R < ∼ 0.25 au while they are positive in the outer region.This is a typical feature of viscous evolution of an isolated accretion disk where the outer region expands outward by receiving the angular momentum from the inner region where the gas accretes inwards.The resolution dependence of the radial mass flux is seen especially in the expanding outer region (0.25 au < ∼ R < ∼ RMRI).The outward radial mass flux around R ∼ 0.4 au is more significant in HighRes than in LowRes run.
When the upper layers of the disk is considered in the estimation of the mass accretion rate, the resolution dependence becomes more significant.For HighRes run, the mass accretion rate is not influenced by the upper layers significantly, or Ṁ H ∼ Ṁ z atm (Figures 23a and 23c).By contrast, LowRes run shows that ⟨W ϕz ⟩ owing to largescale magnetic fields in the upper layers makes Ṁz atm positive in the inner active zone although Ṁ H is negative there.The corresponding structure in Figure 23h is the outgoing gas upper layers in the active zone.This is because the magnetic fields are more coherent for LowRes run than HighRes run.
Next, we investigate the mass transfer around the inner edge of the transition zone.Since the angular momentum transfer mechanisms change significantly both radially and vertically, it is unclear which heights of the gas mainly contribute to ṀRϕ and Ṁϕz at a given radius only from the first and second rows of Figure 23.The color maps of the integrands of ṀRϕ and Ṁϕz , or J Rϕ and J ϕz (see Equations ( 54) and ( 55)) are displayed in Figure 24.
The bottom rows of Figures 23 and 24, or the distributions of ⟨ρvR⟩ and J Rϕ + J ϕz , look quite similar, showing that the radial mass flux is determined by the local torque acting on the gas.Interestingly, Figure 24 HighRes (400 ≤ t/t K0 ≤ 500) LowRes (400 ≤ t/t K0 ≤ 500) −0.004 −0.002 0.000 0.002 0.004 R 2 ρv R t  that small-scale substructures in J Rϕ and J ϕz have complementary spatial distributions to each other, and their sum has a much smoother distribution.
As explained in Section 3.2.2, the hourglass-shaped magnetic field with a steep gradient of B ϕ at the mid-plane is generated mainly by AD (Figure 14b).The middle row of Figures 24 clearly shows that the ϕz stress drives gas accretion at the mid-plane.By contrast, the outflows just above the mid-plane accreting layer identified in the bottom row of Figure 23 are driven by the Rϕ stress.
From Figures 23a and 23b, it is found that the gas within |z| ≤ H moves inward around R = RMRI because the gas accretion rate due to J ϕz dominates over the outflow rate due to J Rϕ .When the contribution of the upper layers (H ≤ |z| ≤ zatm) is taken into account, the net mass flux Ṁ z atm is directed outward (Figures 23c and 23d), leading to the mass supply from the active zone to the transition zone (Section 4.2).
As one moves inward from R = RMRI, the layer with J ϕz < 0 is divided into two layers that sandwitch the active zone in between (the middle row of Figure 24).Inside the active zone, the Rϕ stress moves the gas outward.Just above the two accreting layers, the outflows are driven mainly by the ϕz stress, although the Rϕ stress also contributes.
The surface density radial profiles are shown in Figures 23e and 23f.Roughly speaking, the gas is accumulated around the inner-edge of the transition zone by outward turbulent diffusion in the outer region of the active zone.
Looking more closely at the inner edge of the transition zone in Figure 23e, there are two density peaks around R ∼ 0.43 au and R ∼ 0.52 au for HighRes run.The inner peak is formed by the viscous expansion of the active zone while the outer peak is formed by the wind launched around the inner-edge of the transition zone (Section 3.2.2).The gap between the two peaks is located at the boundary between the expanding layer due to the MRI turbulence and the mid-plane accreting layer (Figure 23h).For LowRes run, the inner density peak is slightly closer to the central star than for HighRes run since the outward radial mass flux around R ∼ 0.4 au is smaller for LowRes run (Figures 23c  and 23d).25a and 25b).The net gas accretion flows through the disk with |z| ≤ H and |z| ≤ zatm are almost zero.One can identify the gas flow converging to the rings both in Figures 25a and 25b.As will be shown in Section 4.2.1, gas accretion in higher latitudes to the inner boundary continues in the late evolution.

Later Evolution
The behaviors around the inner edge of the transition zone shown in Figure 25 are similar to those shown in Figure 23.Inside the disk (R ∼ RMRI, |z| ≤ H), the mass is transferred inward by the ϕz torque due to the hourglass-shaped magnetic fields.When the upper layers of the disk (|z| ≤ zatm) is considered, the inward angular momentum transfer is almost compensated by the outward angular momentum transfer owing to the Rϕ torque.In Section 4.2, we will see that the absence of net angular momentum transport around the inner edge of the transition zone results in almost zero radial mass transfer from the active zone to the transition zone (Figure 26).
For Rη ≤ R < ∼ R ch , surface gas accretion flows are seen in H < |z| < zatm in Figure 25d.The surface gas accretion flows are anti-symmetric with respect to the mid-plane because the current sheet is not at the mid-plane but at the lower AD dead-zone boundary (see Figure 7c).Figures 25a and 25b show that this accretion flow is driven by the magnetic braking owing to the coherent magnetic field (⟨W ⟩ ϕz ).Since the outer region of the transition zone does not have gas accretion, the gas is accumulated around R = Rη (Figure 25c).
The mass accretion rate in the coherent zone is determined by ⟨W ϕz ⟩.Using the ratio f ϕz = B ϕ /Bz and Equation ( 55), the mass accretion rate driven by magnetic braking is estimated to be Ṁ = 8πR ΩK where βz = 8πρ mid c 2 s,mid /B 2 z .Equation ( 57) predicts that Ṁ ∼ 10 −3 in the code units at R = 1 au and Ṁ decreases in proportion to ∝ R −1/4 since Σ ∝ R −1 (Equation ( 5)), where we use the fact that the radial magnetic flux transport from the transition zone decreases βz in the coherent zone by a factor of two from the initial plasma beta 10 4 (Section 3.2.3)and f ϕz is estimated to ∼ 5 from the simulation result.
Beyond R ∼ R ch , the surface gas accretion flows disappear and mid-plane gas accretion is driven by the magnetic torque exerted at the mid-plane (Figures 25a and  25d) since ηO becomes small enough for the toroidal field to be amplified inside the disk (Figure 8c).

Time Evolution of the Masses of the Three Zones
In this section, we investigate how much mass is transferred between the active, transition, and coherent zones and is ejected from or falls onto the disk in the vertical directions.We measure the mass fluxes passing through surfaces enclosing each of the active, transition, and coherent zones.In each zone, the inner and outer radii are denoted by R− and R+, respectively ( R− = R buf and R+ = RMRI for the active zone, R− = RMRI and R+ = Rη for the transition zone, and R− = Rη and R+ = R ch for the coherent zone).The upper and lower boundaries in the vertical directions are z = +zatm(R) and z = −zatm(R), respectively.

The Active Zone
Figure 26a shows that ( Ṁac)tot is always negative, indicating that the total mass of the active zone decreases with time.The time evolution of ( Ṁac)R−, ( Ṁac)R+, and ( Ṁac)z can be divided into the early and late phases.The former shows quasi-periodic oscillations and the latter is characterized by quasi-steady state (Figure 37).The late  phase corresponds to the development of the ring structures.
In the early phase (t < ∼ 1300 tK0), the mass of the active zone is lost both vertically and radially with significant time variations.In the radial directions, the mass of the active zone flows out more from R = R+ = RMRI than from R = R− = R buf .The mass ejection rate by the disk wind ( Ṁac)z is comparable to ( Ṁac)R+.The mass transport rate through R = R+ is investigated in Section 4.2.2 in detail.This indicates that the vertical mass loss is as important as the radial mass transport.
In the later phase (t > ∼ 1300 tK0), ( Ṁac)R− and ( Ṁac)R+ approach zero.The vertical mass flow mainly removes the mass of the active zone.Where does the mass ejected vertically from z = ±zatm(R) go?Since the radial velocity in |z| > zatm above the active zone is negative (the bottom panels of Figure 23 and Figure 25), the disk wind falls onto the center.In order to take into account mass accretion at high latitudes at R = R buf , we measure the mass accretion rate ( Ṁac)buf through the sphere r = R buf , where we take into account only the regions where ρ > 10 −2 to remove the funnel regions.Figure 27 shows that ( Ṁac)buf is comparable to ( Ṁac)z in all the times for LowRes and HighRes runs.This indicates that most mass ejected vertically falls onto the center (Takasao et al. 2018).

The Transition Zone
The mass flux from the active zone to the transition zone ( Ṁtr)R− = ( Ṁac)R+ shows significant quasi-periodic time variations in the early phase of t < ∼ 1300tK0; time periods with large and small ( Ṁtr)R− are periodically repeated.The rate of the mass supply passing through R = R− = RMRI is positive most of the time.This is consistent with the mass transfer rate shown in Figures 23.The time variability of ( Ṁtr)R− is attributed to the MRI activity in the active zone.For t > ∼ 1300tK0, the quasi-periodic time variations in ( Ṁtr)R− disappear and ( Ṁtr)R− gradually decreases with time.This corresponds to the decrease in the MRI activity in the active zone shown in Figure 37.
Figure 26b shows that the mass transfer rate ( Ṁtr)R+ is positive and does not show a significant time variation.This is consistent with the fact that the gas is transferred from the dead zone to the transition zone by the quasisteady surface accretion flow.
Figure 26b shows that the vertical mass supply ( Ṁtr)z stays almost zero.This does not mean that there are no mass transfer at each radius.To investigate the radial profile of the mass ejection rate by the disk wind, in Figure 28 we plot the mass ejection rate from the surfaces of z = ±zatm from R to R + dR that is given by where the gas is ejected from the surfaces of z = ±zatm when d Ṁwind /d ln R > 0. Figure 28 shows that the mass is lost in RMRI ≤ R ≤ 0.6 au and is supplied in the outer region.This profile can be understood by the velocity field above the transition zone.Most amount of the gas outflowing from the inner part falls back within the transition zone owing to the failed disk wind (Figure 13f).

The Coherent Zone
Unlike the active and transition zones, the change rate of the coherent zone mass (Rη ≤ R ≤ R ch ) appears to approach zero in the long-term evolution (Figure 26c).The mass supply through R = R+ = R ch almost balances with the mass loss through R = R− = Rη and disk wind.We note that our calculation time may not be long enough for ( Ṁco)R+ to reach a quasi-steady state.About half of the mass supplied through R = R+ = R ch is ejected by the disk wind, and the remaining half is transferred to the transition zone since ( Ṁco)z ∼ ( Ṁco)R− (Figure 26c).Ferreira (1997) defines "the ejection index" as which indicates a mass ejection efficiency because Ṁwind obtained from the radial integration of d Ṁwind /d ln R increases more rapidly compared with the accretion rate Ṁ z atm as the radial extent considered in the integration for larger ξ.
Figure 28b shows that Ṁ z atm ∼ d Ṁwind /d ln R, indicating that significant mass loss occurs in our simulations.The ejection index is about unity, or ξ ∼ 1 (Equation ( 60)).By the angular momentum conservation and the induction equation, ξ and the magnetic lever arm λ ≡ (RA/R0) 2 are related, where R0 is the radius at a radius where the disk wind is launched and RA is the Alfvén radius of the radius, or ξ = [2 (λ − 1)] −1 (Pelletier & Pudritz 1992).From the fact that ξ ∼ 1, the magnetic lever arm is estimated to λ ∼ 1.5.The lever arm is much shorter than expected in a typical magneto-centrifugal wind (Ferreira 1997).Disk winds with short lever arms were reported; λ ∼ 1.4 (Béthune et al. 2017) and λ ∼ 1.15 (Bai 2017).Bai et al. (2016) pointed out the wind kinematics is mainly controlled by the plasma beta around the wind base.Supposing that the wind base is at z = zatm, the plasma beta is β0 × exp(−(zatm/H) 2 ) ∼ 20.For such a high β, the magneto-centrifugal winds (Blandford & Payne 1982) do not occur since the magnetic fields are too weak to corotate with the Kepler rotation at the wind base within the Alfvén radius.Instead, the thermal and magnetic pressure gradient forces accelerate the winds with the short lever arm λ ∼ 1.5.

Magnetic Flux Transport
We define the following two kinds of the magnetic fluxes.One is the magnetic flux threading the northern hemisphere at a given radius r that is given by The other is the magnetic flux threading the mid-plane from R = rin to R, The divergence-free condition in the northern hemisphere gives where R is identical to r at θ = π/2, because all the magnetic field lines passing through the northern hemisphere at the inner edge and the mid-plane from R = rin to R penetrate through the northern hemisphere at r. From Equation (63), one obtains where Rout = rout is the outermost cylindrical radius of the simulation box.Equation (64) shows how the total magnetic flux Ψ rad (rout) is distributed between Ψ rad (rin) and Ψ mid (Rout).
The time evolution of Ψ rad (rin) and Ψ mid (Rout) is shown in Figure 29.Because Ψ mid (Rout) is much larger than Ψ rad (rin) at the initial condition, Figure 29  the results of LowRes and HighRes shows that the magnetic flux accretes onto the inner boundary more rapidly for HighRes than for LowRes run.
The increasing rate of Ψ rad (rin) decreases with time except for a sudden increase around t = 1000tK0, which corresponds to the accretion of the innermost ⟨Bz⟩ concentration onto the inner boundary of the simulation box around t = 1000tK0 (Figure 19).In the late evolution (1500tK0 ≤ t ≤ 2000tK0), the magnetic flux at the inner boundary hardly changes while the gas accretion onto the inner boundary continues (Figure 27).This decrease in the increasing rate of Ψ rad (rin) has been reported in Beckwith et al. (2009) and Takasao et al. (2018).
In Sections 4.3.1,4.3.2, and 4.3.3,we will investigate the radius dependence of the magnetic flux transport in the active, transition, and coherent zones, respectively.The trajectories of the magnetic flux are shown by the contours of Ψ rad shown in Figure 19.It is clearly seen that the behaviors of the magnetic field transport are different between the three zones.We define a drift speed of the magnetic flux vB as follows: ⟨vB⟩ ≡ ⟨E ϕ ⟩/⟨Bz⟩ (65) (Guilet & Ogilvie 2012;Bai & Stone 2017).

The Active Zone
The trajectories of the magnetic flux and the color map of ⟨Bz⟩H in Figure 19a show that the flux concentrations become prominent for t > ∼ 500tK0, and the contours of Ψ rad converges to the multiple ⟨Bz⟩ concentrations.
First we investigate the magnetic flux transport before developing the flux concentrations.We divide the active zone into two regions, the inner region with R < 0.25 au and the outer region with 0.25 au ≤ R ≤ 0.4 au because the magnetic flux moves inward at R < 0.25 au.The outermost radius is set to R = 0.4 au because beyond R = 0.4 au the MRI turbulence is affected by OR (Figure 35).In each of the inner and outer regions, we average ⟨vB⟩ in the radial direction as follows: In order to investigate what determines the drift speed, we decompose ⟨(E ϕ )I⟩H = ⟨vRBz⟩H − ⟨vzBR⟩H into the laminar component and the turbulent component where Although the mass-weighted average is taken when we calculate the mean values of the velocity in Section 2.7, the volume-weighted average should be applied for the velocity when the flux transport is discussed.To distinguish the volume-weighted mean velocity, the subscript "v" is used, or ⟨vR,z⟩H,v.The radially-averaged drift velocities originating from the laminar and turbulent components of ⟨(E ϕ )I⟩ are defined by replacing ⟨(E ϕ )I⟩H with (⟨(E ϕ )I⟩H ) lami and (⟨(E ϕ )I⟩H ) turb in Equation ( 66), and are denoted by ( vB) lami and ( vB) turb , respectively.
In the early evolution (t ≤ 500 tK0), the magnetic flux moves inward in the inner active zone (R < 0.25 au) and moves outward in the outer active zone (0.25 au ≥ R ≤ 0.4 au) as shown in Figures 30a and 30b which display the time evolution of vB and the contributions from the laminar and turbulent components.In the inner active zone, the turbulent (E ϕ )I mainly drives inward drift of the magnetic flux at a speed of ∼ 5 × 10 −5 vK.Although the resolution dependence of vB is found in t ≤ 300 tK0, it disappears in 300 < ∼ t/tK0 < ∼ 500.The magnetic flux transport is not determined only by the dynamics inside the disk with |z| < H, but is affected by the upper atmospheres.The poloidal magnetic field structures are shown in Figure 31.Just below the boundaries (z ∼ ±0.05 au at r = 0.1 au) between the funnel regions and the atmospheres, the poloidal magnetic fields are dragged toward the inner boundary.After the dragged poloidal fields accrete onto the inner boundary, loop-like poloidal   fields penetrating the disk form (Beckwith et al. 2009).Their edges are at the inner boundaries in the southern and northern hemispheres.When the loop poloidal fields disappear through accretion and/or magnetic reconnection, the net flux at the inner boundary Ψ rad (rin) increases.
In the long-term evolution in the active zone, the turbulent and laminar contributions nearly cancel each other out and their summation, or vB remains low (Figure 30c).The laminar and turbulent toroidal electric fields tend to move the magnetic flux outward and inward, respectively.The drift velocities oscillate with time and takes both positive and negative values for t ≥ 1000 tK0, indicating that there is almost no net magnetic transport.Such a quasisteady state has been obtained in previous global simulations (Beckwith et al. 2009;Takasao et al. 2018).
In the early evolution of the the outer active zone, outward drift of the magnetic flux is driven mainly by ⟨(E ϕ )I⟩ turb .Both ⟨(E ϕ )I⟩ lami and ⟨(E ϕ )I⟩ turb exhibit quasi-periodic oscillations whose phases are correlated with those of ⟨B ϕ ⟩ oscillations in the upper atmospheres.Unlike LowRes run, HighRes run shows that the laminar component plays an important role in the outward flux transport because the outward velocity in the outer active zone is larger for HighRes than for LowRes (Figure 23a).
Figure 30d shows that the quasi-periodic oscillations continue until t ∼ 1500 tK0 beyond which the quasi-periodic MRI activities disappear (Figure 37).Although vB is positive on average in the early phase, the net drift velocities become around zero in 500 ≤ t/tK0 ≤ 1000.In addition, unlike in the inner active zone, ( vB) lami and ( vB) turb do not cancel each other out, and they oscillate in the same phase in the outer active zone.In the later phase with t > ∼ 1500 tK0, ( vB) lami and ( vB) turb do not completely cancel each other out, and the magnetic flux drifts inward at vB ∼ 10 −4 vK.The magnetic flux drifts inward at a faster speed than the gas that accretes at a speed of ∼ 10 −6 vK.Jacquemin-Ide et al. ( 2021) conducted a simulation (model SEp) similar to ours and obtained a drift speed consistent with our results.

The Transition Zone
In Section 3.2.3,we found that the rapid outward transport of the magnetic flux is caused by OR and AD.In this section, we investigate the long term evolution of the magnetic flux around the transition zone.Figure 32b shows that the outer edge of the region where ⟨Bz⟩H ∼ 0 expands with time because the magnetic flux in the coherent zone drift outward as will be shown in Section 4.3.3.A typical drift speed of the outer edge that is identified at the radius where ⟨Bz⟩H = ⟨Bz,ini⟩H is estimated to about 10 −3 vK from Figures 32a and 32b.This drift speed is almost identical to that in the coherent zone (Figure 33).
The flux concentration around R ∼ 0.5 au, where a density gap exists (Figure 25c), weakens with time because ⟨vB⟩H is positive around R ∼ 0.5 au (Figure 32a).Further time evolution would eliminate the flux concentration around R ∼ 0.5 au, resulting in inward expansion of the region where ⟨Bz⟩H ∼ 0. However, the inward expansion would not proceed furthermore because of rapid decreases in ηO and ηA (Figure 2).

The Coherent Zone
Figure 33 shows the time evolution of the radial profiles of ⟨vB⟩H /vK.For R < ∼ 3.5 au, although ⟨vB⟩H /vK is affected by the quasi-periodic disturbances from the active zone, the drift velocity is positive on average.Beyond R > ∼ 3.5 au, ⟨vB⟩H is steady.When the time average over 1000tK0 ≤ t ≤ 2000tK0 is taken, ⟨vB⟩H,t becomes almost constant, and is about ⟨vB⟩H,t ∼ 10 −3 vK.This value is comparable to those in Bai & Stone (2017) and Bai (2017).Figure 34a shows that the drift velocity vB = ⟨E ϕ ⟩t/⟨Bz⟩t does not depend on z significantly, indicating that the magnetic fields drift outward keeping their shapes.This feature is also consistent with Bai & Stone (2017) and Bai (2017).Since the magnetic field is almost completely decoupled with the gas inside the coherent zone where ΛO ≪ 1 and ΛA ≪ 1, the magnetic flux transport is determined above the coherent zone.We focus on the regions just above the OR dead zone where the surface gas accretion occurs (Figure 34b in H < ∼ |z| < ∼ 2H).The surface gas accretion tries to carry the magnetic flux inward since ⟨(E ϕ )I ⟩t ∼ ⟨vRBz⟩t is negative.Diffusion due to AD, which plays an important role in the surface accretion layers (ΛA ∼ 1), tries to drift the magnetic flux outward.Figure 34a shows that the diffusion due to AD slightly dominates over the accretion, resulting in the outward transport of the magnetic flux.
Above the surface accretion layers, the wind drags the vertical field outward (⟨(E ϕ )I ⟩ > 0).Unlike in the surface accretion layers, the diffusion due to AD is not sufficient to dominate over the wind drag because ΛA increases with z (Figure 34c).The drift velocity around the base of the wind (|z| ∼ 2.5H) is close to that in the surface accretion layers.
Inside the OR dead-zone boundary (⟨ΛO⟩t < 1), (⟨E ϕ ⟩)O and (⟨E ϕ ⟩)A almost cancel each other out, and the magnetic field structure inside the disk are adjusted so that {(⟨E ϕ ⟩)O + (⟨E ϕ ⟩)A} /⟨Bz⟩ is comparable to the drift velocity in the region with ΛO > 1.

The Active Zone
As shown in Figure 11, the prominent multiple rings and flux concentrations are found in the active zone.They are not transient structures but are maintained at least until t/tK0 = 2000.It should be noted that not all of the simulations conducted in previous studies have shown ring formation; some simulations (Bai & Stone 2014;Hawley 2001;Steinacker & Papaloizou 2002;Dzyurkevich et al. 2010;Jacquemin-Ide et al. 2021) form rings and some (Suzuki & Inutsuka 2014;Takasao et al. 2018;Zhu & Stone 2018) do not.
In Section 3.2.1,we examined two ring formation mechanisms proposed before, the viscous instability (Suzuki 2023) and the wind-driven instability (Riols & Lesur 2019).The viscous instability may be a possible ring structure formation mechanism in our simulations.However, our analyses are insufficient to fully understand the ring formation in the ideal MHD.The viscous instability found in Suzuki (2023) forms transient ring structures, unlike our results.The possibility of the wind-driven instability is rejected for our simulation, but we note that our result could be affected by the inner boundary condition.
We still do not know detailed mechanisms and what determines whether multiple-rings form or not.We will address the formation of ring structures in MRI active disks by performing global simulations with a wide range of parameters in the future.

The Transition and Coherent Zones
Previous studies investigating the inner dead-zone edge with global non-ideal MHD simulations consider only OR (Lyra & Mac Low 2012;Dzyurkevich et al. 2010;Flock et al. 2017).Our simulation shows that AD significantly changes the structure around the inner dead-zone edge.
Without AD, no rapid magnetic flux transport in the transition zone occurs.As shown in Section 3.2.3,if a finite Bz exists in the initial condition, AD amplifies ⟨BR⟩, leading to the enhancement of ⟨E ϕ ⟩O.Since ηO increases rapidly with R in the transition zone (Figure 2), the magnetic flux is efficiently transported outward within a few rotations at R = RMRI.
The previous studies focusing on outer disks show that AD triggers spontaneous magnetic flux concentration and form ring-like substructures in disks (Bai & Stone 2014;Béthune et al. 2017;Suriano et al. 2018Suriano et al. , 2019;;Cui & Bai 2021).Although the formation mechanism of the spontaneous flux concentration is still unclear, the previous papers suggest that a driving mechanism of the flux concentration is related with the fact that AD steepens the gradients of magnetic fields around magnetic nulls (Brandenburg & Zweibel 1994).Suriano et al. (2018) found that AD generates a pinched magnetic field (|BR| ≫ |Bz|) at the mid-plane where the toroidal magnetic field is reversed.Then, magnetic reconnection is triggered, resulting in magnetic flux transport.Cui & Bai (2021) identified that the direction of B ϕ is reversed in the radial direction where the magnetic flux is concentrated.Another property of AD that may form substructures in disks is that J ϕ,⊥ and J ϕ can have opposite sign when B ϕ ≫ BR,Bz.Béthune et al. (2017) pointed out AD works as anti-diffusion.This behavior is confirmed also in our simulations in Figures 21  and 34.
Although our simulations show the Bz concentrations around the regions where the sign of B ϕ is reversed, prominent gap structures do not form in the coherent zone.This may be caused by the B dependence of ηA.Significant flux concentrations require that ηA increases as ∝ B 2 .As shown in Figure 4, ηA is almost independent of B for Rη < ∼ R < ∼ R ch where the dominant negative charge carrier is dust grains.In addition to the Ohmic diffusion, this may be one of the reasons why flux concentration does not occur there.
For R > R ch , since the physical conditions are similar to those in Suriano et al. (2019), AD should generate the flux The reason why our results do not show the flux concentrations is that the simulation time is too short for them to develop; 2000 rotations at the inner boundary corresponds to ∼ 8 rotations at R = 4 au.

Implications for Long-term Evolution
Although we conducted LowRes run until t = 2000tK0 = 40 yr, it is much shorter than a typical disk lifetime.In this section, we discuss a lone-term evolution of the disk inferred from our results.

Gas Evolution
First, we consider a long-term evolution by assuming the magnetic flux is fixed.If the gas accretion rate shown in Figure 27 is maintained ( Ṁ ∼ 5 × 10 −4 in the simulation unit corresponding to Ṁ = 2.8 × 10 −7 M⊙ yr −1 ), the depletion time scale of the active zone is ∼ 4 × 10 5 yr if there is no mass supply.This may be a formation mechansim of transition disks (e.g., van der Marel 2023).
What would be the mass transport between the transition and active zones?As discussed in Sections 3.2.2 and 4.1, the direction of the mass transfer is determined by the two competing mass flows; the mid-plane accretion and the outflow in the high latitudes.If the mass of the active zone becomes sufficiently small, the mid-plane accretion from the transition to active zone will dominate over the outflows from the active zone to the transition zone.This leads to inward mass transfer from the transition zone to the active zone.
In the coherent zone, the surface gas accretion continues to accumulate the gas in the outer edge of the transition zone.If the gas accretion rate is kept to Ṁ = 10 −3 in the simulation units (Figure 25), the mass corresponding to the initial mass of the active zone is transferred to the transition zone in 2 × 10 5 yr.The accumulated gas is not directly transferred to the active zone since gas accretion inside the transition zone does not occur because of no net magnetic flux.Thus, the transferred mass forms a ring structure whose surface density contrast is about 10 3 if the radial width is ∼ 0.2 au at 2 × 10 5 yr.In reality, the ring will become vorticies by the Rossby-wave instability (Lovelace et al. 1999;Ono et al. 2016Ono et al. , 2018)).
Gas accretion from the inner edge of the transition zone to the active zone may increase RMRI when the surface density becomes low enough for cosmic rays and/or ionizing radiation to reach the mid-plane.If RMRI gets close enough to Rη, the gas accumulated at the outer-edge of the transition zone may be transferred to the active zone.

Magnetic Flux Evolution
Existence of the transition zone significantly affect the magnetic flux evolution in PPDs.From our results, the following important conclusion can be drawn; the magnetic fluxes cannot drift from the coherent zone to the active zone.
The existence of the transition zone affects on the evolution of the active zone.The magnetic flux initially possessed by the active zone falls to the central star or is transported to the transition zone.The magnetic flux transported to the transition zone does not return to the active zone.This indicates that no magnetic flux may be supplied to the active zone from R > RMRI.If this is the case, the magnetic flux in R ≤ RMRI decreases with time, leading to a decrease in the gas accretion rates onto the central star.As the process regulates the magnetic flux accretion onto the central star, it can be important to determine magnetization of young stellar objects.
The outer-edge of the region where Bz ∼ 0 moves following the drift velocity of the magnetic flux in the coherent zone as shown in Figure 32.We here define the transition zone as the region where Bz ∼ 0. If vB is positive as shown in our simulations and the previous simulations (Bai & Stone 2017;Bai 2017;Gressel et al. 2020), the outer edge of the transition zone moves outwards.Since the gas accretion will be suppressed in the region where Bz ∼ 0, the inner-edge of the coherent zone where the surface gas accretion occurs moves outward.If vB is negative, the outer edge of the transition zone will not move and the magnetic flux will be accumulated just outside the transition zone.
Long-term evolution of magnetic flux transport in the coherent zone remains elusive.The drift velocities have about an order of magnitude differences between the previous studies (Bai & Stone 2017;Bai 2017;Gressel et al. 2020).Lesur (2021) shows that the radial drift velocities strongly depend on the spatial distributions of ηO and ηA using simple self-similar solutions.Although their results show only outward drifts, we need to investigate how drift velocities depend on the disk parameters and the size distribution of dust grains.If the drift velocity can be negative in a parameter range, the magnetic flux will be accumulated at the outer edge of the transition zone.We will investigate how drift velocities depend on the disk parameters and the size of dust grains in the future.

Implication for Dynamics of Dust Grains
The inner dead-zone edge is thought to be a possible site where dust grains accumulate since a pressure bump is created owing to outward viscous expansion of the inner active zone (e.g., Kretke et al. 2009).So far, accumulation of dust grains at the inner dead-zone edge was investigated mainly by using one-dimensional (1D) model taking into account viscous evolution of disks and coupling of the gas and dust grains (Kretke et al. 2009;Pinilla et al. 2016;Ueda et al. 2019).The 1D models can compute the evolution for much longer periods than three-dimensional global disk simulations.
From our simulations considering realistic spatial distributions of ηO and ηA, we found that the transition zone is created with a finite width at the inner dead-zone edge.The gas accumulates from both inner and outer edges of the transition zone, producing the two ring structures at both the edges (Figure 22a).The structure found in this study is quite different from that predicted by the 1D models that show a single pressure bump because the 1D models assume a step-function-like distribution of the viscous parameter α and neglect the gas surface accretion in the dead zone.Ueda et al. (2019) derived a critical value of the α parameter α dead in the dead zone below which the dust pileup occurs, where v frag is the minimum collision speed that triggers fragmentation of dust grains.For silicate aggregates, v frags spans 1 − 10 m s −1 (Blum & Wurm 2000;Wada et al. 2013).The dead zone in Ueda et al. (2019) corresponds to the transition zone in our simulations.We note that Ueda et al. (2019) considers α dead as the turbulent viscosity.In our simulations, the velocity dispersion in the transition zone is generated by the sound waves originating from the active zone (Section 3.2.5).The α parameter estimated from the Reynolds and Maxwell stresses, αtr, cannot be substituted into α dead in Equation ( 69) because the velocity dispersion in the transition zone δvtr does not satisfy δvtr/cs ∼ √ αtr , which is valid in the MRI turbulence.Instead, α dead in Equation ( 69) should be considered as a measure of the random velocity dispersion, and is replaced with (δvtr/cs) 2 .In terms of the velocity dispersion δvtr, Equation ( 69) is rewritten as Equation ( 70) shows that in order for dust grains to accumulate without being radially diffused, the Mach number of the velocity dispersion in the transition zone is required to be lower than 0.02 for v frag = 1 m s −1 .Our results show that two rings form at both the edges of the transition zone, or R = RMRI and Rη although the setting of Ueda et al. (2019) yields a single ring forms at the inner edge of the dead zone.Thus, accumulation of dust grains potentially occurs in the two rings.Figure 22a shows that the Mach numbers of the velocity dispersion ⟨δMtot⟩ at the mid-plane is as large as 0.1 at the inner edge of the transition zone, and decreases with radius to reach ∼ 0.01 at the outer edge of the transition zone.Comparing our results with Equation (70) shows that accumulation of dust grains may be suppressed in the ring near the inner edge of the transition zone.In the ring at the outer edge of the transition zone, dust grains may accumulate if v frag is well above 1 m s −1 .The velocity dispersion in the transition zone is expected to be proportional to that in the active zone.To accumulate dust grains in the rings, plasma beta larger than 10 4 would be preferred.

Implications for Protoplanetary Disks Around Solar-type Stars
In this paper, PPDs around Herbig Ae stars were focused on.It is worth discussing implications for PPDs around solar-type stars from our results.
Since the inner dead-zone edge RMRI is roughly determined by T ∼ 10 3 K (Desch & Turner 2015), RMRI is predicted to be RMRI ∼ 0.6 au κR 5 cm 2 g −1 where κR is the Rossland mean opacity, and κR = 5 cm 2 g −1 corresponds to the opacity of dust grains with 0.1 µm and a dust-to-gas mass ratio of 10 −2 .In derivation of Equation ( 71), we adopt the standard α disk model (Shakura & Sunyaev 1973;Gammie 1996) where the vertical optical depth is set to κRΣ/2.Mori et al. (2019) pointed out that RMRI derived from Equation (71) assumes the existence of MRI turbulence.If no MRI turbulence exists, the disk is laminar.Joule heating is much less effective than viscous heating (Mori et al. 2019), and the gas temperature becomes comparable to that determined by the stellar irradiation (Chiang & Goldreich 1997).Thus, the MRI should be self-sustained in a sense that the MRI turbulent heating provides a sufficient amount of charged particles to drive the MRI.It is however unclear whether this self-sustaining cycle works stably.If the MRI is suppressed, the inner edge of the dead zone moves closer to the central star.
The existence of the transition zone should be critical in global disk evolution as discussed in Section 5.2.Although criteria to form the transition zone are unclear, if the active zone exists stably, Section 3.2.3shows that the formation of the transition zone requires ηO to increase rapidly across R ∼ RMRI and ηA to be large enough to amplify BR and B ϕ around the inner edge of the dead zone.Further investigations are required to reveal whether the transition zone exists and whether the active zone can exist stably by using non-ideal MHD numerical simulations with radiative transfer.

Summary
In this paper, we performed global non-ideal MHD simulations of inner regions of a protoplanetary disk around an intermediate star focusing on the behaviors of the inner dead zone boundaries by taking into account the Ohmic resistivity (OR) and ambipolar diffusion (AD).The radial extent spans from R = 0.1 au to R ∼ 10 au.
The three characteristic radii, RMRI, Rη, and R ch , are defined from the spatial distribution of ηO and ηA.The first radius RMRI = 0.5 au corresponds to the radius outside of which the MRI is suppressed around the mid-plane.The second radius Rη = 1.2 au shows the radius around which ηO and ηA are maximized.In the disk where R < Rη, the thermal ionization of metals determine the ionization degree.Beyond the last radius R ch = 3 au, the diffusion coefficients are small enough for the accretion flow to be driven at the mid-plane.
We found that the dead zone identified by ΛO < 1 and ΛA < 1 can be divided into two regions separated by Rη.The transition zone in RMRI < R < Rη is discovered in this paper and has different properties than those of the conventional dead zones while the coherent zone in R ≥ Rη shares the main characteristics with the conventional dead zone.
We summarize our work as follows: 1.The overall physical properties of the active zone (R < RMRI) are consistent with those found in the literature.
• Prominent ring structures develop in the long-term evolution (Figure 11), independent of resolution.
The vertical field strength is anti-correlated with the surface density.Although the driving mechanisms are uncertain, the viscous instability (Suzuki 2023) would be promising (Sections 3.2.1 and 5.1.1).The development of ring structures in MRI disks is an open question that requires further investigation.• The magnetic flux transport in the active zone depends both on time and radius (Section 4.3.1).After the ring structures and flux concentrations have developed, the magnetic flux drifts inward.A typical drift velocity in the active zone is about ∼ 10 −4 vK. 2. The physical properties of the coherent zone (R ≥ Rη) are also consistent with those found in the literature.
• The behaviors of the angular momentum transfer are characterized by R = R ch .For Rη < R < R ch , surface gas accretion is driven just above the OR dead zone.The mass loss rate of the disk wind is comparable to the mass accretion rate.For R ≥ R ch , ηO is low enough for the toroidal field to be amplified inside the disk, leading to mid-plane gas accretion (Sections 3.1 and 4.1).• The drift velocity of the magnetic flux is determined by the inward drag owing to the surface gas accretion and outward diffusion owing to AD (Section 4.3.3).Our simulations show that the latter dominates over the former, and the magnetic flux drift outwards at a speed of ∼ 10 −3 vK.
Although it belongs to the conventional dead zone defined by Elsässer nubmers, their physical properties are completely different from those in the conventional dead zone.Our findings regarding the transition zone are as follows: • The vertical magnetic flux is rapidly transported outward both by OR and AD, forming a magnetic gap where the vertical field is almost zero (Section 3.2.3).Without AD, the rapid flux transport would not occur.• Unlike the conventional dead zone (Bai & Stone 2013), the lack of the vertical field suppresses the surface gas accretion.Since neither MRI nor coherent magnetic fields extract the angular momentum, gas accretion is largely suppressed in the transition zone.
• The mass of the transition zone monotonically increases with time by the gas inflow both from the inner-and outer-edges of the transition zone during the calculated period.As a result, rings form in both the edges of the transition zone.The outflow launched around the inner edge of the transition zone fails to escape the disk but falls onto the outer part of the transition zone, resulting in zero net vertical mass flux.
Our simulations have demonstrated that the complicated structures formed by the non-ideal MHD effects affect the mass, angular momentum, and magnetic flux transport as well as dust accumulation.We are planning to perform wider parameter search to investigate such structures systematically.In addition, we will include the Hall effect in our future work.

Fig. 1 .
Fig. 1. (Left) Initial setup of the simulation.The color map indicates the initial density profile in the code units, and the black lines indicate the streamlines of the magnetic fields.The red (blue) line corresponds to the condition ΛO = 1 (ΛA = 1), which show the OR (AD) dead zone boundary.The magenta lines show z = ±zatm(R) above which the warm atmospheres exist.The right panel is the same as the left panel but it shows the zoom-in of the region enclosed by the black rectangle shown in the left panel.

Fig. 2 .
Fig. 2. (a) Radial profiles of the diffusion coefficients ηO and ηA at the midplane at t = 0 (solid) and t = 500 tK0 (dashed) when the MRI turbulence has been saturated.The diffusion coefficients are shown in the code units.The thin dashed and dotted lines show the diffusion coefficients at t = 0 derived by chemical reactions with thermal and non-thermal ionization, respectively.Panel (b) is the same as Panel (a) but for the Els ässer numbers ΛO and ΛA at the mid-plane at t = 0 and t = 500 tK0.For both panels, the three vertical dashed lines correspond to RMRI, Rη, and R ch from left to right.

Fig. 3 .
Fig. 3. Two-dimensional distributions of (a) ηO and (b) ηA in the code units at t = 500 tK0 when the MRI turbulence is saturated.The white dashed lines show the contours of ηO and ηA at 10 −4 and 10 −2 .The conditions of ΛO = 1 and ΛA = 1 are displayed by the red and blue lines in Panel (a) and (b), respectively.In each panel, the two black dashed lines correspond to z = H(R) and z = 2H(R), and the green dashed line show z = zatm(R) above which the warm atmospheres exist.

Fig. 4 .
Fig. 4. (a) Dependence of ηA on the field strength.The radial profiles of ηA at the mid-plane are displayed by changing the field strengths.The values of log 10 (B/Bini(R)) are labelled on the lines.The vertical axis is shown in the code units.The black line shows the radial profile of ηA in the initial condition.(b) The field strengths Bcri where d ln ηA/d ln B = 1.0 and 0.5.For field strengths larger than Bcri, the B-dependence of ηA becomes weak.

Fig. 5 .
Fig. 5.The mesh structures of LowRes and HighRes runs in the (r, θ) plane.Each cell corresponds to a mesh block with (Nr, N θ ) = (24, 16).The red and blue lines represent the boundaries of the dead zones for OR and AD (ΛO = 1 and ΛA = 1), respectively.The green line indicates the scale height z = ±H.

Fig. 6 .Fig. 7 .
Fig. 6.A schematic picture showing our main findings with respect to the time-averaged gas dynamics of the disk.The orange arrows represent the gas motion, and the blue lines show the poloidal magnetic field lines.

Fig. 8 .
Fig. 8.The same as Figure 7 but the radial range 0 ≤ R ≤ 6 au.The three vertical dashed lines show R = RMRI, Rη, and R ch from left to right.An animation of this figure is available.

Fig. 9 .
Fig. 9. Color maps of (a) the density, (b) the radial mass flux, (c) the toroidal magnetic field, (d) the vertical magnetic field averaged over |z| ≤ H at t = 2000 tK0.All the quantities are shown in the code units.In each panel, the three dashed circles correspond to RMRI, Rη, and R ch from the center outward.An animation of this figure is available.

Fig. 11 .
Fig. 11.Face-on views of the density (left) and Bz (right) averaged over |z| < H(R) at t/tK0 = 1000 (top) and 2000 (middle).(Bottom) Face-on views of the density averaged over |z| < H(R) at t = 500tK0 for LowRes (left) and HighRes (right) runs.All the quantities are shown in the code units.

Fig. 12 .
Fig. 12.(a) Normalized vertical mass loss rate in the active zone Cw defined in Equation (37) as a function of radius.The solid and dashed lines correspond to the radial profiles of Cw where Cw > 0 and Cw < 0, respectively.The thin dashed lines show the radial profiles of ⟨β z,mid ⟩ −1 t whose value is shown on the right vertical axis.The data are averaged over 900 ≤ t/tK0 ≤ 1000 and 1900 ≤ t/tK0 ≤ 2000.(b) Color maps of ⟨ρ⟩ at t = 2000 tK0.The black lines show the streamlines of the poloidal magnetic field averaged over ϕ.The two dotted magenta lines correspond to z = ±zatm.

Fig. 13 .
Fig. 13.Spatial Structures of R⟨B ϕ ⟩t (top),⟨B ϕ ⟩ 2 t /⟨B 2 ϕ ⟩t (middle), and R⟨ρv R ⟩t (bottom) in the transition zone.The left panels show the enlarged view around the inner dead-zone boundaries, and the right panels display the entire region of the transition zone.The results are averaged over 0 ≤ ϕ < 2π and 1500 ≤ t/tK0 ≤ 2000, and all the quantities are shown in the code units.In each panel, the red and blue lines represent the OR and AD dead-zone boundaries, respectively.The black lines show the streamlines of the poloidal magnetic fields (the top and middle panels) and poloidal velocity fields (the bottom panels).In the bottom panels, the streamlines of the poloidal velocity fields are not shown inside the OR dead-zone, and the region enclosed by the green dashed lines is used to evaluate the mass transfer in the transition zone (Section 4.2).In Panel (f), the four streamlines originating from (R, z) = (0.45 au, 0.06 au), (0.45 au, 0.07 au), (0.45 au, 0.073 au), and (0.45 au, 0.08 au) are shown by the thick lines.

Fig. 15 .
Fig. 15.Time sequence of the ϕ-averaged field structure near the dead zone boundary from t/tK0 = 50 to 90.In each panel, the colors show the R⟨B R ⟩ map, and the black lines correspond to the streamline of the poloidal magnetic fields.The values of R⟨B R ⟩ are shown in the code units.

Fig. 16 .
Fig. 16.Time evolution of the vertical profiles of ⟨B R ⟩, ⟨B ϕ ⟩, and ⟨Bz⟩ at R = 0.7 au from t/tK0 = 50 to 80 in an interval of 10.All the quantities are shown in the code units.

Fig. 19 .
Fig. 19.Radius-time diagrams of (a)⟨Bz⟩ H /⟨Bz,ini⟩ H and (b)R⟨B ϕ ⟩ H , where ⟨Bz,ini⟩ H is ⟨Bz⟩ H at t = 0.The vertical red dashed lines correspond to RMRI, Rη, and R ch from left to right.The black dashed lines show the contours of Ψ rad , which is defined in Equation (61).They correspond to the trajectory of the magnetic flux originated at each equally-spaced position at t = 0.All the quantities are shown in the code units.

Fig. 20 .
Fig. 20.Time variations of B ϕ at three different radii.The blue line shows∆⟨B ϕ ⟩ ≡ ⟨B ϕ ⟩(z = −2.5H)+ ⟨B ϕ ⟩(z = 2.5H) at R = 0.3 au, and ∆⟨B ϕ ⟩ is multiplied by 0.05 to facilitate comparison with the other two lines.The orange and green lines correspond to ⟨B ϕ ⟩ H at R = 0.7 au and 1.5 au, respectively.All the quantities are shown in the code units.

Fig. 22 .
Fig. 22. Radial propagation of the sound waves in the whole disk.(a) The Mach numbers in the R, ϕ, and z directions at the mid-plane as a function of R (the left vertical axis).The radial profile of the mid-plane density is shown by the gray line (the right vertical axis).(b) The radial profiles of the energy fluxes of the sound waves.The blue, orange, and green lines correspond to ⟨(F R )+⟩, ⟨(F R )−⟩, and ⟨F R ⟩ at the mid-plane, respectively.In all the panels, all the data are averaged over 1000 ≤ t/tK0 ≤ 2000, and the vertical dashed lines denote R = RMRI and Rη from left to right.All the quantities are shown in the code units.

Fig. 23 .
Fig. 23.Radial dependence of the mass accretion rate in the early stage for HighRes (left) and LowRes (right).All the quantities shown in this figure are averaged over 400tK0 ≤ t ≤ 500tK0.The first row compares the mass accretion rate (pink) integrated over |z| ≤ H and their predictions (black) given by Equation (53).The blue and orange lines correspond to Ṁ H Rϕ and Ṁ H ϕz , respectively.The second row is the same as the first row but the mass accretion rate is integrated over |z| ≤ zatm.The radial profiles of the surface densities are shown in the third row.In the fourth row, the color maps of R 2 ⟨ρv R ⟩t are displayed in the plane of (log 10 R, z/H), and the black solid and dashed lines correspond to z = ±H and z = ±zatm.In each panel, the two vertical dotted lines show R = RMRI and R = Rη.In the first and second row, to reduce significant spatial fluctuations of the mass accretion rates Ṁ H,z atm Rϕ , Ṁ H,z atm ϕz

Fig. 24 .
Fig. 24.Color maps of J Rϕ , J ϕz , and J Rϕ + J ϕz averaged over 400tK0 ≤ t ≤ 500tK0 for HighRes (left panels) and LowRes runs (right panels).The dotted lines correspond to z = ±zatm, respectively.All the quantities are shown in the code units.
Figure 25 is the same as Figure 23 but showing the late evolution for LowRes run.In the active zone, as shown in Section 3.2.1, the ring structures and flux concentrations develop.Unlike in the early evolution, the outer regions of the active zone (R ∼ 0.4 au) do not show outward mass flux Ṁ H > 0 and Ṁ z atm > 0 in the long-term evolution (Figures

Fig. 25 .
Fig. 25.The same as Figure 23 but the quantities are averaged over 1500 ≤ t/tK0 ≤ 2000 for LowRes run.In Panels (a) and (b), to reduce significant spatial fluctuations of the mass accretion rates Ṁ H,z atm Rϕ , Ṁ H,z atm ϕz

Fig. 26 .
Fig. 26.Mass transfer rate through the four surfaces of (a) the active, (b) transition, and (c) coherent zones for LowRes run.The blue, orange, and green lines correspond to Ṁac,tr,co R+ , Ṁac,tr,co R− , and Ṁac,tr,co z , respectively.The gray lines show Ṁac,tr,co tot .To reduce significant temporal fluctuations of the mass transfer rates, we take the 20tK0 centered moving average for them.All the quantities are shown in the code units.

Fig. 27 .
Fig. 27.Mass transfer in the active zone.Time evolution of ( Ṁac)R−, ( Ṁac)z, and ( Ṁac)buf for (a) LowRes and (b) HighRes runs.To reduce significant temporal fluctuations of the mass transfer rates, we take the 20tK0 centered moving average for them.All the quantities are shown in the code units.

Fig. 28 .
Fig. 28.Vertical mass outflow rate from R to R + dR as a function of radius.Panel (a) compares the results for LowRes and HighRes, and they are averaged over 400 ≤ t/tK0 ≤ 500.Panel (b) displays the time evolution of d Ṁwind /d ln R by showing the result averaged over three different time intervals, 500 ≤ t/tK0 ≤ 1000, 1000 ≤ t/tK0 ≤ 1500, and 1500 ≤ t/tK0 ≤ 2000.The gas accretion rate in the coherent zone is shown by the magenta line.All the quantities are shown in the code units.

Fig. 30 .
Fig. 30.Time evolution of the drift velocities averaged over R < 0.25 au (top) and 0.25 au ≤ R ≤ 0.4 au (bottom) for LowRes.The left and right panels show the early and whole evolution, respectively.For comparison, the results for HighRes are shown by the thin dashed line in the left panels.In each panel, vB lami , vB turb , and vB are shown, and are normalized by the Kepler speed vK averaged over the corresponding radius range.To reduce significant temporal fluctuations of the drift velocities, we take the 20tK0 centered moving average for them.

Fig. 31 .
Fig. 31.Poloidal magnetic fields structure near the inner boundary at t = 500tK0.The color map shows log 10 ⟨ρ⟩, and the streaming lines show the poloidal magnetic fields averaged over 0 ≤ ϕ < 2π.The two dashed lines correspond to z = ±zatm.The densities are shown in the code units.

Fig. 32 .
Fig. 32.Time Evolution of the radial profiles of ⟨vB⟩ H /vK and ⟨Bz⟩ H /⟨Bz⟩ H,ini from t = 1000tK0 to t = 2000tK0.In Panel (a), we do not plot ⟨vB⟩ H /vK in 0.6au ≤ R ≤ 1 au because scatters of the data are significant owing to almost zero ⟨Bz⟩ H .All the quantities are shown in the code units.

Fig. 33 .
Fig. 33.Time sequence of the radial profiles of the drift velocity of the magnetic flux ⟨vB⟩ H normalized by the Kepler speed from t = 1000tK0 to t = 2000tK0.The radial profiles of ⟨vB⟩ H averaged over 1000tK0 ≤ t ≤ 2000tK0 are shown by the red line.The vertical dashed lines correspond to R = Rη and R ch from left to right.

Fig. 34 .
Fig. 34.Vertical profiles of the quantities related to the magnetic field transport at R = 2 au.All the data are averaged over 1000tK0 ≤ t ≤ 2000tK0.(a) The vertical profiles of ⟨(E ϕ )I⟩t/⟨Bz⟩t, ⟨(E ϕ )O⟩t/⟨Bz⟩t, ⟨(E ϕ )A⟩t/⟨Bz⟩t, and ⟨(E ϕ )⟩t/⟨Bz⟩t.(b) The vertical profile of ⟨ρv R ⟩t.(c) The vertical profiles of ⟨ΛO⟩t and ⟨ΛA⟩t.In each panel, the vertical gray dashed lines correspond to the positions where ⟨ΛO⟩t = 1.All the quantities are shown in the code units.

Table 1 .
List of the models.The simulations are conducted up to t end . †
Salvesen et al. (2016)y,Salvesen et al. (2016)found that p β ∼ 0.5.It is unclear what causes p β in our simulation to be smaller than those in local shearing-box simulations, but we note that p β obtained from the spatial variations of ⟨α Rϕ ⟩ H need not necessarily consistent with p β obtained from a volume average of the α parameter in local shearing-box simulations.