Dynamics near the inner dead-zone edges in a proprotoplanetary disk

Abstract

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ässer 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.

1 Introduction

The evolution of protoplanetary disks (PPDs) is controlled by angular momentum transfer. One of the important angular momentum transfer mechanisms is the magnetorotational instability (MRI; Velikhov 1959; Chandrasekhar 1961; Balbus & Hawley 1991). The turbulence driven by MRI transfers angular momentum radially. MRI has been extensively investigated in local shearing-box simulations (e.g., Hawley et al. 1995) and in global simulations (Armitage 1998; Hawley et al. 2011; Sorathia et al. 2012; Suzuki & Inutsuka 2014; Takasao et al. 2018; Zhu & Stone 2018; Jacquemin-Ide et al. 2021).

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 the 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 by 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 and 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 and 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 non-axisymmetric shadow on the outer disk.

In the global simulations containing the inner dead-zone edge performed by Dzyurkevich et al. (2010) and Lyra and 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 dead-zone edge is still unclear.

In this paper, we perform global three-dimensional non-ideal 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, the coefficients of which 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.

2 Numerical setup and methods

2.1 Basic equations

The resistive MHD equations with OR and AD are given by

and

where ρ is the gas density, |$\boldsymbol {v}$| is the gas velocity, |$\boldsymbol {B}$| is the magnetic field, ψ is the gravitational potential of the central star, |$\boldsymbol {J}$| is the current density, |$\boldsymbol {J}_\perp \equiv \boldsymbol {J} - \boldsymbol {J}\cdot \boldsymbol {B}/|\boldsymbol {B}|$| is the electric current density perpendicular to the local magnetic field, P is the gas pressure, |${\bf \textsf {T}}$| is the stress tensor, |$\rm{{\bf \textsf {I}}}$| is the identity matrix, and η_O and η_A are the diffusion coefficients for OR and AD, respectively. They are defined in sub-subsection 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 with a temperature distribution as given in equation (10). We note that a locally isothermal equation of state potentially drives the vertical shear instability (VSI; Nelson et al. 2013). In our simulations, the VSI does not grow because of insufficient simulation time, as discussed in sub-subsection 5.1.3.

2.2 Models and methods

2.2.1 Initial condition

The initial surface-density profile of the disk is given by

where Σ₀ = 500 g cm⁻² 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⁴ 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\pi R_*^2 \sigma _\mathrm{SB} T_*^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) = \sqrt{2}c_\mathrm{s,mid}/\Omega _\mathrm{K}(R)$|⁠, where |$c_\mathrm{s,mid}=\sqrt{k_\mathrm{B} T_\mathrm{mid}/\mu m_\mathrm{H}}$| is the sound speed at the mid-plane, μ = 2.33 is the mean molecular weight, m_H is the hydrogen mass, and |$\Omega _\mathrm{K}= \sqrt{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) = (μm_H|$/$|k_B)[ϵ_mid(r)²|$/$|2]r²Ω_K(r)² 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 |$\Sigma _r(r,\theta ) = \int _{r_\mathrm{min}}^r \rho (r^{\prime },\theta )dr^{\prime }$| is smaller than Σ_FUV, the gas temperature increases to T_atm(r) = (μm_H|$/$|k_B)[ϵ_atm(r)²|$/$|2]r²Ω_K(r)², where ϵ_atm(r) = 0.2(r|$/$|1 au)^1/4. The initial temperature profile smoothly connects between T_d and T_atm around Σ_r ∼ Σ_FUV as follows:

and

The values of Σ_FUV are highly uncertain, and Perez-Becker and Chian (2011) found that FUV layers have a thickness of ∼0.01–0.1 g cm⁻² in plane-parallel models. We adopt Σ_FUV = 0.3 g cm⁻². This is three times larger than the upper limit given in Perez-Becker and Chian (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_\mathrm{boost}^{-1}$| 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 β₀ = 10⁴. The initial condition of our simulation is displayed in figure 1.

The height of the base of the atmospheres at a given R is denoted by z = z_atm(R) above which the gas temperature is more than twice the mid-plane temperature (figure 1).

2.2.2 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, and Z_je is the charge and n_j 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 m_j 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 H₂ number density is 3 × 10⁻⁷. The thermal ionization degree (x_e)_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 |$Tf_\mathrm{boost}^{-1}$| in the Saha equation. When the gas temperature exceeds ∼1000 K, the thermal ionization of potassium provides a sufficiently 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 (x_e)_T is the ionization degree obtained from the Saha equation for potassium ionization and the second term indicates that the electrons are provided by photo-ionization 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, Nishi, and Umebayashi (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 |$Tf_\mathrm{boost}^{-1}$|⁠. For simplicity, the single size dust grain with a radius of |$a=0.1\, \mu$|m and an internal density of 2 g cm⁻³ is considered. In this paper, the dust-to-gas mass ratio is set to 10⁻⁴, 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, the ionization rate of which is given by

where Σ_CR = 96 g cm⁻² (Umebayashi & Nakano 1981) and |$\Sigma _\theta ^+(r,\theta )$| [|$\Sigma _\theta ^-(r,\theta )$|] 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 ζ₁ = 6.0 × 10⁻¹¹ s⁻¹, Σ_{X, a} = 3.6 × 10⁻³ g cm⁻², α = 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.

2.3 Basic properties of the diffusion coefficients

2.3.1 Definition of the dead zone boundaries

The non-ideal MHD effects weaken the growth rate of the MRI when the dimensionless Elsässer numbers

are smaller than unity, where |$v_\mathrm{A}=B/\sqrt{4\pi \rho }$| is the Alfvén speed (Sano & Miyama 1999; Blaes & Balbus 1994; Wardle 1999; Kunz & Balbus 2004; Bai & Stone 2011).

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 and Miyama (1999) and Bai and 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 sub-subsection 3.2.2.

Hereafter, the dead-zone boundaries for OR and AD are defined as Λ_O = 1 and Λ_A = 1, respectively.

2.3.2 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 R_MRI, 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 does not since Λ_O ∝ B² and Λ_A ∝ B⁰, where we use the fact that η_O ∝ B⁰ and η_A ∝ B² (sub-subsection 2.3.4). In figure 2b, we plot the radial profiles of Λ_O and Λ_A at |$t=500\, t_{\mathrm{K0}}$| when the MRI turbulence has been saturated. At R = R_MRI = 0.5 au, Λ_O becomes unity, and R_MRI 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 ionization 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 (sub-subsection 2.3.4).

2.3.3 Two-dimensional distributions of the diffusion coefficients

The two-dimensional distributions of η_O and η_A in the (R, z) plane at |$t=500\, t_{\mathrm{K0}}$| 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 ∼ R_MRI, while high-latitude regions are the first places 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 = z_atm(R) is well above the dead-zone boundaries.

2.3.4 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πen_e) ∝ B² (Salmeron & Wardle 2003), where β_i and β_e are the Hall parameters of ions and electrons, respectively, and n_e 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 sub-subsection 5.1.2. Development of sharp structures in the magnetic null due to AD is reported in Brandenburg and 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² is realized at the weak field limit β_e ≪ 1 (Xu & Bai 2016). For β_e ≳ 1, η_A is not necessarily proportional to B² 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 and Bai (2016).

Figure 4a shows how the radial profile of η_A at the mid-plane changes as the field strength is increased from the initial value B_ini(R). As B|$/$|B_ini increases, η_A does not increase in proportion to B², 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 strengths, the B dependence of η_A become weak, leading to inefficient formation of substructures due to AD.

2.4 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 & Stone 2008). The second-order Runge–Kutta–Legendre super-time stepping technique is used to calculate the magnetic diffusion processes (Meyer 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 r_in ≤ r ≤ r_out, 0 ≤ θ ≤ π, and 0 ≤ ϕ < 2π, where r_in = 0.1 au and r_out = 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 meridional 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. Each rectangle enclosed by the black lines in figure 5 consists of (N_r, N_θ) = (24, 16). Table 1 lists the models considered in this paper, where |$t_{\mathrm{K0}}=2\pi \sqrt{r_\mathrm{in}^3/GM_*}$| is the Kepler rotation period at the inner boundary of the computation box. The 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 (N_r, N_θ, N_z) = (1536, 1024, 1024). The disk scale height H(R) is resolved by 19 (R|$/$|R_MRI)^1/4 cells.

The LowRes run was conducted until |$2000\, t_{\mathrm{K0}}$| (table 1). As will be shown in section 1, the resolution of the LowRes run is not high enough to obtain the converged MRI turbulence. For comparison, we conduct the HighRes run, where the cells in the active zone are further refined in one more higher level (the lower panel of figure 5). The 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 the 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 the 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 the LowRes run, referring the results of the 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:

This density floor works only in the polar regions. In the other regions, the disk winds keep the densities greater than the floor value.

2.5 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(r_in), v_ϕ,ini], and v_r and v_θ are set to zero, where v_ϕ,ini is the initial toroidal velocity. We set the boundary conditions for B_r and B_ϕ so that B_r ∝ r⁻² and B_ϕ ∝ r⁻¹.

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_ϕ, B_r, and B_ϕ so that v_ϕ ∝ r^−1/2, B_r ∝ r⁻², and B_ϕ ∝ r⁻¹. The radial velocity is set by using the zero-gradient boundary condition only for v_r > 0, otherwise it is set to zero.

In order to mitigate the artificial effect of the inner boundary conditions, we introduce a buffer region in r_in ≤ r ≤ r_buf, where r_buf − r_in = 0.025 au is the width of the buffer region.

The treatment in the buffer region is similar to that used in Takasao et al. (2018). The poloidal velocities are damped according to the following equation:

where τ_d is the damping timescale, |$\tau _\mathrm{d}^{-1} = f_\mathrm{rad}(r) f_v(v_r) \tau _\mathrm{d0}^{-1}$|⁠,

and

f_d,min = 1, f_d,max = 100, and v_r,c = −0.05v_K,in. The other physical quantities are unchanged in the damping layer.

2.6 Normalization units

In this paper, all the physical quantities are normalized by the following three quantities: the length scale R₀ = 1 au, the time scale |$t_0 = \sqrt{R_0^3/GM_*}=0.1\:$|yr, and the surface density Σ₀ = 500 g cm⁻². The gas temperature is normalized by T₀ = (μm_H|$/$|k_B)(R₀|$/$|t₀)² = 6.3 × 10⁵ K. Hereafter, we use normalized physical quantities unless otherwise noted, except for R and z, which are shown in units of the astronomical unit.

2.7 Averaged quantities

We define the 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 ρ, ρv_z, |$B_\phi ^2/8\pi$|⁠, or B_R. For velocities and kinetic energies per unit mass, we take the density-weighted average

where a quantity A can be v_z or |$v_R^2/2$|⁠.

In order to eliminate stochastic features originating from the MRI turbulence, 〈A〉_ϕ is averaged over t₁ ≤ t ≤ t₂ 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 |z| ≤ z_b:

Using |$\langle \boldsymbol {v} \rangle$| and |$\langle \boldsymbol {B} \rangle$|⁠, we define the turbulent components of velocities and magnetic fields as variations from the ϕ-averaged values as follows:

where the turbulent components (⁠|$\delta \boldsymbol {v}$| and |$\delta \boldsymbol {B}$|⁠) satisfy |$\langle \rho \delta \boldsymbol {v} \rangle =0$| and |$\langle \delta \boldsymbol {B} \rangle =0$|⁠.

3 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.

3.1 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 = 2000t_K0 for the 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 R_MRI ≤ R ≤ R_η “the transition zone”, and it has different properties from the conventional dead zone. We call the outer region of R > R_η “the coherent zone”, and it has the same properties as the conventional dead zone.

The overall properties of the active, transition, and coherent zones are briefly summarized in sub-subsections 3.1.1, 3.1.2, and 3.1.3, respectively. In sub-subsection 3.1.4, we discuss the influence of the active zone on the outer regions.

3.1.1 Active zone (R ≤ R_MRI)

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 in appendix  2). 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 (figure 37d).

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 B_z are anti-correlated; the magnetic fluxes are concentrated in the density gaps. Possible formation mechanisms of these ring structures are investigated in sub-subsection 3.2.1 and are discussed in sub-subsection 5.1.1.

3.1.2 Transition zone (R_MRI < 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 R_MRI ≲ R ≲ R_η (sub-subsection 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 (sub-subsection 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 (sub-subsection 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 〈ρv_R〉 is averaged over time (sub-subsection 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 ∼ R_MRI 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 ≤ R_MRI) and the coherent zone (R ≥ R_η). Details will be investigated in sub-subsections 4.1 and 4.2.2. No B_z concentration occurs while the density peak forms around the outer-edge of the transition zone (figure 9d). This is because the magnetic flux drifts outward from the transition zone to the coherent zone as shown in sub-subsections 4.3.2 and 4.3.3.

3.1.3 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 −B_z∂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 z-symmetry adopted in the initial condition is broken (figure 7c). This is because η_O and η_A are so large inside the disk that a current sheet cannot exist inside the coherent zone (figure 2a). As a result, a current sheet is lifted either upward or 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 by 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 mid-plane at R > R_ch is amplified, resulting in the OR dead-zone 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. (2018, 2019) and also around the inner edge of the transition zone as discussed in sub-subsections 3.2.3. The magnetic tension force −B_z∂B_ϕ|$/$|∂z drives gas accretion at the mid-plane in R > R_ch (figure 8b).

3.1.4 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 sub-subsection 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 (sub-subsection 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 (sub-subsection 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=2000\, t_{\mathrm{K0}}$|⁠, 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 〈B_z〉 has a concentration around R ∼ R_ch. Concentrations in 〈B_z〉 propagate outward following quasi-periodic inversion of the sign of the toroidal magnetic field. Details will be investigated in sub-subsection 3.2.4.

3.2 Main findings

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

3.2.1 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 anti-correlated with the radial profile of the vertical magnetic flux averaged over the scale height 〈B_z〉_H; the rings (gaps) of Σ correspond to the gaps (rings) of 〈B_z〉_H.

The top and middle panels of figure 11 show the face-on views of the density and B_z averaged over |z| < H(R). At |$t=1000\, t_{\mathrm{K0}}$|⁠, 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=2000\, t_{\mathrm{K0}}$|⁠, 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=1000\, t_{\mathrm{K0}}$|⁠.

The resolution-dependence on the density distributions are shown in the bottom panels of figure 11. The density is lower for the HighRes run than for the LowRes run in the active zone because the HighRes run exhibits faster viscous evolution (figure 36 in appendix 1). Despite the difference in the density structures, ring structures are also found in the HighRes run. This suggests that if a long-term simulation was performed for the HighRes run, the ring structures would develop.

Several formation mechanisms of ring structures in the ideal MHD have 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 sub-subsection 5.1.1, we further discuss the origin of the ring structures in our simulations.

3.2.1.1 Viscous instability

Viscous instability in MRI-active disks has been investigated in Riols and Lesur (2019) and Suzuki (2023). Riols and Lesur (2019) assumed that the α parameter (Shakura & Sunyaev 1973) is proportional to |$\beta _z^{-p_\beta }$|⁠, where β_z = 8π〈〈P〉〉|$/$|〈〈B_z〉〉² is the plasma beta measured by the net vertical field 〈〈B_z〉〉 since α is anti-correlated with β_z in local shearing-box simulations (Hawley et al. 1995; Suzuki et al. 2010; Salvesen et al. 2016; Scepi et al. 2018), where 〈〈A〉〉 is a volume-average of A. The instability criterion derived by Riols and Lesur (2019) is p_β > 1. Suzuki (2023) extended the linear analysis in Riols and Lesur (2019) by considering the dependence on 〈〈ρ〉〉 and 〈〈B_z〉〉 separately, or |$\alpha \propto \langle \langle \rho \rangle \rangle ^{-p_\rho } \langle \langle B_z \rangle \rangle ^{-p_{B_z}}$|⁠. The instability criterion derived by Suzuki (2023) is p_ρ > 1. |$p_{B_z}$| does not contribute to the instability criterion, but it changes the growth rate. We note that the instability criteria derived in Riols and Lesur (2019) and Suzuki (2023) are consistent since |$\beta _z^{-p_\beta } \propto \langle \langle \rho \rangle \rangle ^{-p_\beta } \langle \langle B_z \rangle \rangle ^{2p_\beta }$|⁠.

We investigate whether our results satisfy the instability criteria given by Riols and Lesur (2019) and Suzuki (2023). The α parameter is defined by using the thermal pressure, Reynolds stress, and Maxwell stress averaged over |z| ≤ H as follows:

The plasma beta averaged over z ≤ H is defined as |$\langle \beta _{z} \rangle _H = 8\pi \langle P \rangle _H/\langle B_z \rangle _H^2$|⁠. The radial profiles of 〈α_Rϕ〉_H are fitted by the two fitting functions |$\alpha _\mathrm{ana,\beta }\propto (\langle \beta _z \rangle _H/\beta _0)^{-p_\beta }$| and |$\alpha _\mathrm{ana,\rho }\propto (\langle \rho \rangle _H/\rho _0)^{-p_\rho }$|⁠, ¹ where β₀ = 10⁴ and ρ₀ = ρ_mid,ini(R = 0.3 au). The best-fitting functions are shown in table 2. Figure 10c (figure 10d) compares 〈α_Rϕ〉_H and the best-fitting functions α_ana,β (α_ana,ρ) at both |$t=1000\, t_{\mathrm{K0}}$| and |$2000\, t_{\mathrm{K0}}$|⁠. Figures 10c and 10d show that both fitting functions α_ana,β and α_ana,ρ reproduce 〈α_Rϕ〉_H reasonably well.

Table 2 shows that the results of the LowRes run do not satisfy the Riols and Lesur (2019) instability criterion, while they do satisfy the Suzuki (2023) instability criterion since p_β ∼ 0.3–0.4 and p_ρ > 1.² This suggests that the viscous instability may contribute to the ring formations in our simulations if the α parameter is determined by the gas density.

Although a quantitative argument why p_ρ becomes greater than unity is still missing, an anti-correlation between 〈ρ〉_H and 〈B_z〉_H as seen in figures 11a and 11b may be the key to understanding the strong density dependence of the α parameter. Similar anti-correlations between the density and the vertical field were reported, for instance, in Bai and Stone (2014), Suriano et al. (2019), Jacquemin-Ide, Lesur, and Ferreira (2021), and Suzuki (2023). Using |$\alpha _\mathrm{ana,\beta } \propto \langle \beta _z \rangle _H^{-p_\beta } \propto \langle \rho \rangle _H^{-p_\beta } \langle B_z \rangle ^{2p_\beta }_H$| and assuming |$\langle B_z \rangle _H\propto \langle \rho \rangle _H^{-q_{B_z}}$|⁠, we obtain |$\alpha _\mathrm{ana,\beta }\propto \langle \rho \rangle _H^{-p_\beta (2q_{B_z}+1)}$|⁠. When 〈B_z〉_H and 〈ρ〉_H are anti-correlated (⁠|$q_{B_z}\\gt 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. The LowRes run shows that |$q_{B_z} \sim 0.43$| for |$t=1000\, t_{\mathrm{K0}}$| and ∼0.95 for |$t=2000\, t_{\mathrm{K0}}$|⁠. One obtains that |$-\partial \ln \alpha _\mathrm{ana,\beta }/\partial \ln \langle \rho \rangle _H = p_\beta (2q_{B_z}+1) \sim 0.76$| for |$t=1000\, t_{\mathrm{K0}}$| and ∼0.93 for |$t=2000\, t_{\mathrm{K0}}$|⁠, and both values are not consistent with p_ρ shown in table 2.

3.2.1.2 Wind-driven instability

Riols and 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 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 have a more 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 v_n stands for the velocity component perpendicular to the ±z_atm surfaces (figures 1 and 3), ρ_mid is the mid-plane density, and c_s,mid is the mid-plane sound speed.

Assuming that both C_w and α are anti-correlated with β_z and they follow the relations |$C_\mathrm{w}\propto \beta _z^{-p_\mathrm{w}}$| and |$\alpha \propto \beta _z^{-p_\beta }$|⁠, Riols and Lesur (2019) found that the wind-driven instability occurs when p_w > p_β. Since p_β is roughly equal to 0.3–0.4 in our simulations (figure 10c and table 2), p_w needs to be larger than 0.3–0.4 in order for the wind-driven instability to occur.

Figure 12a compares the radial profiles of C_w with the inverse of the time-averaged mid-plane plasma beta |$\langle \beta _{z,\mathrm{mid}} \rangle _t \equiv 8\pi \langle P \rangle _t(z=0)/\langle B_z \rangle _t^2(z=0)$|⁠. C_w is poorly correlated with 〈β_{z, mid}〉⁻¹, 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 |$\langle \beta _{z,\mathrm{mid}} \rangle _t^{-1}$| have local maxima, C_w is negative, indicating that the gas flows into the density gap regions rather than being ejected from the disk in the vertical direction. The vertical inflow into the density gap regions is evident from the arrows in figure 12b. Why are there no outflows from the gap regions where the magnetic field is strong? Riols and 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 z-axis (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| increases, suggesting that the gas is not continuously accelerated along the field lines.

3.2.2 Gas dynamics 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 sub-subsections 3.2.4 and 3.2.5. In this section, we investigate the quasi-steady structure of the transition zone by taking the time average. We will show that AD plays an important role in forming the hourglass poloidal magnetic fields and driving the mid-plane gas accretion around the inner edge of the transition zone (figure 6). Disk winds are launched from higher latitudes around the inner edge of the transition zone. In addition, we will also see that almost no gas accretion is driven in the transition zone. We will quantitatively discuss the gas flows driven by the torques in sub-subsections 4.1 and 4.2.

The left-hand 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 = R_MRI = 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. This is also clearly seen as an increasing trend of B_R with z for a constant B_z across the mid-plane in figure 14a, which presents the vertical distribution of each component of the magnetic field at R = R_MRI. The toroidal magnetic field dominates over the other components, and has a sharp gradient at the mid-plane.

3.2.2.1 Hourglass-shaped poloidal magnetic field

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 hourglass-shaped magnetic fields, we measure their contributions to |$-\langle \left(\nabla \times \boldsymbol {E}\right)_R \rangle$|⁠, which is equal to ∂〈B_R〉|$/$|∂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. Since B_z is almost constant in |z| ≲ 0.01 au, the mechanism that steepens the B_R gradient (figure 14a) creates an hourglass-shaped magnetic field. Figure 14b shows that the AD term tilts the poloidal magnetic field toward the radial direction of the disk (〈(∂B_R|$/$|∂t)_A〉 > 0 for z > 0, and 〈(∂B_R|$/$|∂t)_A〉 < 0 for z < 0), and the ideal MHD term does the same. The vertical profile of 〈B_R〉 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 in quite the opposite way. The AD term works as a diffusion of 〈B_R〉 amplified by the ideal MHD term.

3.2.2.2 Mid-plane gas accretion

The mid-plane accretion flow seen in 0.45 au ≲ R ≲ 0.55 au (figure 13c) is driven by the magnetic torque, −B_z∂B_ϕ|$/$|∂z, due to the pinched magnetic field by AD near the mid-plane. 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 〈(∂B_R|$/$|∂t)_A〉, far from the mid-plane the AD term works as diffusion of 〈B_ϕ〉 amplified by the ideal MHD term.

3.2.2.3 Disk wind launching

Just above the thin mid-plane gas accretion layer, the wind-like gas flows in the outward direction are driven. This is also caused by the magnetic torque of the coherent magnetic fields (Blandford & Payne 1982; Bai et al. 2016, see sub-subsection 4.1.2). 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-hand panels of figure 13 zoom out the region shown in the left-hand panels to cover the entire transition zone. Four streamlines are highlighted by 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_η as a result of 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).

3.2.2.4 Absence of gas accretion in the transition zone

Comparison between figures 13f and 7b shows that the radial mass flux inside the OR dead-zone disappears when the time average is taken. The radial mass flux existing inside the OR dead-zone in figure 7b originates from the sound waves generated from the MRI turbulence in the active zone (sub-subsection 3.2.5; Pucci et al. 2021).

3.2.3 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). We will present that both a steep spatial gradient of η_O and amplification of the toroidal electric field by AD are necessary for efficient outward transfer of the vertical flux. This flux transport is very efficient, and occurs in less time than seven rotations at R = R_MRI.

Figure 15 shows the time sequence of the magnetic field structure. At t = 50 t_K0, the poloidal magnetic field lines are almost vertical inside the AD dead zone at R > 0.6 au. After 1.8 rotations at R = R_MRI (t = 70 t_K0), 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 t_K0, the magnetic loop structure with its center 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 at R = 0.7 au are shown in figure 16. As shown in figure 15, 〈B_R〉 is amplified around the mid-plane during 50 ≤ t|$/$|t_K0 ≤ 70. At the same time, 〈B_ϕ〉 is also amplified. Interestingly, 〈B_z〉 decreases while 〈B_R〉 and 〈B_ϕ〉 are amplified, indicating a rapid transport of the vertical field.

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 |$\boldsymbol {E}_\mathrm{I}=-\boldsymbol {v}\times \boldsymbol {B}$| is the electric field in the ideal MHD, and |$\boldsymbol { E}_\mathrm{O}=\eta _{\mathrm{O}}\boldsymbol { J}$| and |$\boldsymbol { E}_\mathrm{A}=\eta _{\mathrm{O}}\boldsymbol { J}_\perp$| are the electric fields caused by OR and AD, respectively.

We investigate the radial transport of the vertical field. Equation (40) shows that the radial transport velocity of 〈B_z〉 is estimated by 〈E_ϕ〉|$/$|〈B_z〉. Since 〈B_z〉 is positive in most regions at the mid-plane, the sign of 〈E_ϕ〉 represents the direction of the vertical field transport. Figure 17 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 of 〈E_ϕ〉 in R ≲ 0.9 au is attributed to the 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\sim 50\, t_{\mathrm{K0}}$| to |$70\, t_{\mathrm{K0}}$|⁠, representing the rapid outward transport of 〈B_z〉 (figure 16c).

Figure 18 shows that an increase in E_ϕ results in outward magnetic flux transfer. We investigate how 〈E_ϕ〉_O is generated. Because |∂B_R|$/$|∂z| dominates over |∂B_z|$/$|∂R|,

the vertical distribution of B_R is critical for the vertical field transport.

The evolution of 〈B_R〉 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 t_K0, 〈(E_ϕ)_A〉 mainly contributes to 〈E_ϕ〉. Around the mid-plane, ∂〈(E_ϕ)_A〉|$/$|∂z is positive (negative) for z > 0 (<0), leading to amplification of 〈B_R〉 since the signs of ∂〈(E_ϕ)_A〉|$/$|∂z and 〈B_R〉 are the same (left-hand 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 18f). Since 〈η_A〉 is proportional to B² in this region, the gradient of 〈η_A〉 becomes steeper and steeper owing to the amplification of the magnetic field. As both 〈B_R〉 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 ≥ 70t_K0, 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 〈B_R〉.

The gradient of 〈B_R〉 with respect to z becomes steeper as the time passes, leading to the current sheet around the mid-plane becoming 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 and Zweibel (1994).

3.2.4 Quasi-periodic disturbances of magnetic fields beyond the inner dead-zone edge due to the MRI activity

In this section, we will show that the MRI activity drives quasi-periodic variations of the magnetic fields in the transition and coherent zones. The sign of B_ϕ reverses quasi-periodically. Anti-diffusion caused by AD generates a B_z condensation where the sign of B_ϕ is flipped. The magnetic field disturbances do not propagate beyond R = R_ch.

Figure 19 shows the radius–time diagrams of 〈B_z〉_H and 〈B_ϕ〉_H. As shown in sub-subsection 3.2.1, in the active zone (R ≤ R_MRI), the ring and gap structures develop in 〈B_z〉_H and 〈ρ〉_H.

In the transition zone (R_MRI < R < R_η), as shown in sub-subsection 3.2.3, the combination of OR and AD causes the magnetic flux to be transferred outward rapidly in the early evolution (⁠|$\sim\! 100\, t_{\mathrm{K0}}$|⁠). As a result, 〈B_z〉_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 dead-zone 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 a 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 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 〈B_ϕ〉_H variations occur quasi-periodically. The flip of the sign of 〈B_ϕ〉_H occurs almost synchronously in R_MRI ≲ R ≲ R_ch around |$t\sim 150\, t_{\mathrm{K0}}=13.4\, t_\mathrm{K}(R_{\mathrm{MRI}})$|⁠, |$\sim\! 520\, t_{\mathrm{K0}}= 46.5\, t_\mathrm{K}(R_{\mathrm{MRI}})$|⁠, |$\sim\! 1100\, t_{\mathrm{K0}}= 98.4\, t_\mathrm{K}(R_{\mathrm{MRI}})$|⁠, and |$\sim\! 1450\, t_{\mathrm{K0}}= 130\, t_\mathrm{K}(R_{\mathrm{MRI}})$|⁠.

In order to investigate the origin of the quasi-periodic time variation in 〈B_ϕ〉_H in R_MRI ≲ 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, the 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 Δ〈B_ϕ〉 ≡ 〈B_ϕ〉(z = −2.5H) + 〈B_ϕ〉(z = 2.5H) at R = 0.3 au. 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 that shows a correlation 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\sim\! 500\, t_{\mathrm{K0}}$| (figure 19b).

Figure 21a shows the time evolution of the radial profiles of 〈B_R〉, 〈B_ϕ〉, and 〈B_z〉 at the mid-plane. When the time passes from |$t=540\, t_{\mathrm{K0}}$| to |$600\, t_{\mathrm{K0}}$|⁠, the radius where the sign of 〈B_ϕ〉 is reversed moves outward while 〈B_R〉 and 〈B_z〉 do not change significantly. At |$t=600\, t_{\mathrm{K0}}$|⁠, a concentration of 〈B_z〉 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 |$\langle \eta \rangle /\Delta R \sim\! 0.1\:\mbox{au}\, t_{\mathrm{K0}}^{-1}$|⁠, indicating that the structure of 〈B_ϕ〉 moves 1 au per |$10\, t_{\mathrm{K0}}$|⁠. Although this speed is a few times larger than estimated from figure 21a, they are consistent based on the rough estimate.

The 〈B_z〉 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 |$-\boldsymbol {v}\times \boldsymbol {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 exhibits anti-diffusion behavior since OR provides pure diffusion of magnetic fields. As shown in figure 21d, the opposing signs 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 this triggers the concentration of 〈B_z〉.

Comparison between figures 19a and 19b shows that the concentrations of 〈B_z〉_H are located at the radii where the sign of 〈B_ϕ〉_H is flipped for R ≳ R_η. The outward propagation of the concentrations of 〈B_z〉_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 〈B_z〉_H.

3.2.5 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 figures 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, |$\langle {\cal M}_{R,\phi ,z} \rangle$| are as large as ∼0.1. Comparison between the Mach numbers and density in figure 22a shows that the velocity dispersion is anti-correlated with the density; it is larger in the gaps than in the rings.

Inside the transition zone (R_MRI < R < R_η), |$\langle \delta {\cal M}_R\rangle$|⁠, |$\langle \delta {\cal M}_\phi \rangle$|⁠, and |$\langle \delta {\cal M}_z\rangle$| behave differently. |$\langle \delta {\cal M}_R\rangle$| continues to decrease in the transition zone, and reaches ∼10⁻² at R = R_η. By contrast, after the rapid decrease at R ≈ R_MRI, |$\langle \delta {\cal M}_\phi \rangle$| and |$\langle \delta {\cal M}_z\rangle$| become almost constant in the outer part of the transition zone. Around R = R_η, all the components of the Mach numbers have a similar value of ∼10⁻².

In the coherent zone (R ≥ R_η), |$\langle \delta {\cal M}_R\rangle$| and |$\langle \delta {\cal M}_z\rangle$| are of similar magnitude and decrease with increasing radius. |$\langle \delta {\cal M}_\phi \rangle$| 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 |$\langle {\cal F}_R \rangle = \langle ({\cal F}_R)_+ \rangle - \langle ({\cal F}_R)_- \rangle$|⁠.

The radial distributions of |$\langle ({\cal F}_{R})_+ \rangle$|⁠, |$\langle ({\cal F}_R)_- \rangle$|⁠, and |$\langle {\cal F}_R \rangle$| are shown in figure 22b. In the active zone, |$\langle ({\cal F}_R)_+ \rangle$| and |$\langle ({\cal F}_R)_- \rangle$| are comparable; there are equal amounts of outward and inward propagating sound waves.

Outside the active zone, |$\langle ({\cal F}_R)_+ \rangle$| dominates over |$\langle ({\cal F}_R)_- \rangle$|⁠. This clearly shows that the sound waves transport more energy outward. Around R ∼ R_η, |$\langle ({\cal F}_R)_- \rangle$| is enhanced and becomes comparable to |$\langle ({\cal F}_R)_+ \rangle$| 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).

4 Detailed analyses on transfer of mass, angular momentum, and magnetic flux

4.1 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 transfer 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 z_atm.

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 t₁ ≤ t ≤ t₂, one obtains a prediction of the radial mass flux |$\dot{M}_\mathrm{pred}^{z_{\mathrm{b}}}$|⁠, which is given by integrating 2πR〈ρv_R〉 over −z_b ≤ z ≤ z_b as a function of radius as follows:

where

and

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 |$\dot{M}_\mathrm{pred}^{z_\mathrm{b}}$|⁠. The third term comes from the radial mass flux affected by the vertical mass flux.

4.1.1 Early evolution

First, we investigate the angular momentum transfer in the early evolution |$400\, t_{\mathrm{K0}}\le t \le 500\, t_{\mathrm{K0}}$|⁠. 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=450\, t_{\mathrm{K0}}$|⁠, indicating that both the transition and coherent zones have not reached quasi-steady 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.

4.1.1.1 Mass transfer in the active zone

The top and second rows of figure 23 show comparison between the measured and predicted mass transfer rates by integrating over |z| ≤ H and |z| ≤ z_atm, respectively. Equation (53) reproduces the actual mass transfer 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. For both torques, the Maxwell stress plays a more dominant role than the Reynolds stress does. 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 the 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 ≲ R_MRI). The outward radial mass flux around R ∼ 0.4 au is more significant in the HighRes than in the LowRes run.

When the upper layers of the disk are considered in the estimation of the mass transfer rate, the resolution-dependence becomes more significant. For the HighRes run, the mass transfer rate is not influenced by the upper layers significantly, or |$\dot{M}^H\sim \dot{M}^{z_{\mathrm{atm}}}$| (figures 23a and 23c). By contrast, the LowRes run shows that 〈W_ϕz〉 owing to large-scale magnetic fields in the upper layers makes |$\dot{M}_{z_{\mathrm{atm}}}$| positive in the inner active zone although |$\dot{M}^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 the LowRes run than the HighRes run.

4.1.1.2 Mass transfer around the inner edge of the transition zone

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 |$\dot{M}_{R\phi }$| and |$\dot{M}_{\phi z}$| at a given radius only from the first and second rows of figure 23. The color maps of the integrands of |$\dot{M}_{R\phi }$| and |$\dot{M}_{\phi z}$|⁠, or |${\cal J}_{R \phi }$| and |${\cal J}_{\phi z}$| [see equations (54) and (55)] are displayed in figure 24.

The bottom rows of figures 23 and 24, or the distributions of 〈ρv_R〉 and |${\cal J}_{R\phi } + {\cal J}_{\phi z}$|⁠, look quite similar, showing that the radial mass flux is determined by the local torque acting on the gas. Interestingly, figure 24 shows that small-scale substructures in |${\cal J}_{R\phi }$| and |${\cal J}_{\phi z}$| have complementary spatial distributions to each other, and their sum has a much smoother distribution.

As explained in sub-subsection 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 figure 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 = R_MRI because the gas accretion rate due to |${\cal J}_{\phi z}$| dominates over the outflow rate due to |${\cal J}_{R\phi }$|⁠. When the contribution of the upper layers (H ≤ |z| ≤ z_atm) is taken into account, the net mass flux |$\dot{M}^{z_{\mathrm{atm}}}$| is directed outward (figures 23c and 23d), leading to the mass supply from the active zone to the transition zone (subsection 4.2).

As one moves inward from R = R_MRI, the layer with |${\cal J}_{\phi z}\\lt 0$| is divided into two layers that sandwich 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.

4.1.1.3 Ring structures around the inner edge of the transition zone

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 the 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 (sub-subsection 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 the LowRes run, the inner density peak is slightly closer to the central star than for the HighRes run since the outward radial mass flux around R ∼ 0.4 au is smaller for the LowRes run (figures 23c and 23d).

4.1.2 Later evolution

Figure 25 is the same as figure 23 but showing the late evolution for the LowRes run. In the active zone, as shown in sub-subsection 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 |$\dot{M}^H\\gt 0$| and |$\dot{M}^{z_{\mathrm{atm}}}\\gt 0$| in the long-term evolution (figures 25a and 25b). The net mass transfer rates through the disk with |z| ≤ H and |z| ≤ z_atm are almost zero. One can identify the gas flow converging to the rings both in figures 25a and 25b. The fact that gas accretion does not occur around the disk mid-plane may be due to the ring structures. Jacquemin-Ide, Lesur, and Ferreira (2021) investigated an MRI disk (model SEp) similar to our active zone and found that the radial migration of the ring structures is not seen. However, the physical mechanism is unclear and needs to be investigated. We note that the insufficient resolution of the LowRes run may cause the almost-zero accretion rate. The gas accretion continues in the higher latitudes although it halts inside the disk. This is evidence that the mass transfer rate integrated over |z| ≤ 5z_atm (⁠|$\dot{M}^{5z_{\mathrm{atm}}}$|⁠) takes negative values except around the rings (see also figure 27).

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 ∼ R_MRI, |z| ≤ H), the mass is transferred inward by the ϕz torque due to the hour-glass-shaped magnetic fields. When the upper layers of the disk (|z| ≤ z_atm) is considered, the inward angular momentum transfer is almost compensated by the outward angular momentum transfer owing to the Rϕ torque. In subsection 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 (figures 26a and 26c).

For R_η ≤ R ≲ R_ch, surface gas accretion flows are seen in H < |z| < z_atm 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 transfer rate in the coherent zone is determined by 〈W_ϕz〉. Using the ratio f_ϕz = B_ϕ|$/$|B_z and equation (55), the mass transfer rate driven by magnetic braking is estimated to be

where |$\beta _z = 8\pi \rho _\mathrm{mid} c_\mathrm{s,mid}^2/B_{z}^2$|⁠. Equation (57) predicts that |$\dot{M} \sim\! -10^{-3}$| in the code units at R = 1 au and |$|\dot{M}|$| decreases in proportion to ∝R^−1/4 since Σ ∝ R⁻¹ [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⁴ (sub-subsection 3.2.3) and f_ϕz is estimated to be ∼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).

4.2 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 on to the disk in the vertical direction.

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₊ = R_MRI for the active zone, R₋ = R_MRI 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 = +z_atm(R) and z = −z_atm(R), respectively.

The time evolution of the total masses of the active, transition, and dead zones is given by

where the subscripts “ac,” “tr,” and “co” stand for the active, transition, and coherent zones, respectively, the mass increasing rates passing through the surfaces R = R₋ and R₊ are denoted by |$(\dot{M}_\mathrm{ac,tr,co})_{R-}$| and |$(\dot{M}_\mathrm{ac,tr,co})_{R+}$|⁠, respectively, and the mass increasing rates from the surfaces z = −z_atm and z = +z_atm are added and denoted as |$(\dot{M}_\mathrm{ac,tr,co})_z$|⁠. The signs of |$(\dot{M}_\mathrm{ac,tr,co})_{R-}$|⁠, |$(\dot{M}_\mathrm{ac,tr,co})_{R+}$|⁠, and |$(\dot{M}_\mathrm{ac,tr,co})_{z}$| are determined so that they are positive when the total mass inside the enclosed region increases owing to its contribution. Since the surface at R = R_MRI (R = R_η) is shared between the active and transition zones (the transition and coherent zones), |$(\dot{M}_\mathrm{ac})_{R+} = -(\dot{M}_\mathrm{tr})_{R-}$| [|$(\dot{M}_\mathrm{tr})_{R+} = -(\dot{M}_\mathrm{co})_{R-}$|] is satisfied. Figures 26a, 26c, and 26d show the time evolution of |$(\dot{M}_\mathrm{ac,tr,co})_{R+}$|⁠, |$(\dot{M}_\mathrm{ac,tr,co})_{R-}$|⁠, |$(\dot{M}_\mathrm{ac,tr,co})_z$|⁠, and |$(\dot{M}_\mathrm{ac,tr,co})_\mathrm{tot}$|⁠.

4.2.1 The active zone

Figure 26a shows that |$(\dot{M}_\mathrm{ac})_\mathrm{tot}$| is always negative, indicating that the total mass of the active zone decreases with time. The time evolution of |$(\dot{M}_\mathrm{ac})_{R-}$|⁠, |$(\dot{M}_\mathrm{ac})_{R+}$|⁠, and |$(\dot{M}_\mathrm{ac})_z$| can be divided into the early and late phases. The former shows quasi-periodic oscillations and the latter is characterized by a quasi-steady state (figure 37). The late phase corresponds to the development of the ring structures.

4.2.1.1 Early evolution

In the early phase (⁠|$t\lesssim 1000\, t_{\mathrm{K0}}$|⁠), the mass of the active zone is lost both vertically and radially with significant time variations. The fact that the mass of the active zone is lost through R = R_MRI [|$(\dot{M}_\mathrm{ac})_{R+}\\lt 0$|] is consistent with the positive mass transfer rate at R = R_MRI in figures 23c and 23d. In the radial directions, the mass of the active zone flows out more from R = R₊ = R_MRI than from R = R₋ = R_buf. The mass ejection rate by the disk wind |$(\dot{M}_\mathrm{ac})_z$| is comparable to |$(\dot{M}_\mathrm{ac})_{R+}$|⁠. This indicates that the vertical mass loss is as important as the radial mass transport.

Why does |$(\dot{M}_\mathrm{ac})_{R+}$| exhibit quasi-periodic variation even though the MRI disappears around R = R_MRI (figure 24)? The quasi-periodic variations generated in the atmospheres of the active zone disturb the magnetic field around R = R_MRI. Figure 26b shows the time variations of 〈B_ϕ〉 at z = 2H at two different radii. As a representative radius in the active zone, R = 0.3 au is chosen because the quasi-periodic variations caused by the MRI are seen in 0.2 au ≲ R ≲ 0.4 au. At R = 0.3 au, 〈B_ϕ〉 at z = 2H is positive for a longer time than it is negative because it is amplified from negative B_R that is generated by dragging of the magnetic fields toward the central star (Suzuki & Inutsuka 2014; Takasao et al. 2018). The amplified positive 〈B_ϕ〉 suddenly breaks and the sign of 〈B_ϕ〉 is flipped. Afterwards, the positive 〈B_ϕ〉 is amplified again. This cycle is repeated in a quasi-periodic manner in the active zone (figure 37). When the sign of 〈B_ϕ〉 reverses and reverts in the active-zone atmospheres, the disturbances propagate outwards and affect the magnetic fields around R = R_MRI. This can be seen in the fact that 〈B_ϕ〉(R = R_MRI, z = 2H) suddenly drops slightly after 〈B_ϕ〉(R = 0.3 au, z = 2H) becomes negative in figure 26b. The quasi-periodic variations in the magnetic fields at R = R_MRI produce the fluctuations of |$(\dot{M}_\mathrm{ac})_{R+}$|⁠.

Comparison between figures 26a and 26b indicates that the time variation of the vertical mass transfer rate |$(\dot{M}_\mathrm{ac})_z$| is also correlated with that of 〈B_ϕ〉 at R = 0.3 au. This is because the gas is ejected radially obliquely upward when the sign of 〈B_ϕ〉 returns to positive in the northern hemisphere.

The variations of |$(\dot{M}_\mathrm{ac})_{R-}$| are not quasi-periodic and correlate little with |$(\dot{M}_\mathrm{ac})_{R+}$| and |$(\dot{M}_\mathrm{ac})_{R+}$|⁠. The variations of |$(\dot{M}_\mathrm{ac})_{R-}$| are also likely to be produced by the MRI activity, but it may be affected by the buffer region and the inner boundary conditions (subsection 2.5).

4.2.1.2 Later evolution

The time evolution of |$(\dot{M}_\mathrm{ac})_{R+}$| transitions to a stationary state because the magnetic field around R = R_MRI is no longer subject to quasi-periodic disturbances from the atmospheres of the active zone (figure 26b). As mentioned in sub-subsection 4.1.2, |$(\dot{M}_\mathrm{ac})_{R+}$| approaches zero because the mid-plane gas accretion is almost compensated by the outflow in higher latitudes.

The vertical mass flow |$(\dot{M}_\mathrm{ac})_z$| mainly determines the time evolution of the mass of the active zone because |$(\dot{M}_\mathrm{ac})_{R-}$| and |$(\dot{M}_\mathrm{ac})_{R+}$| approach zero (figure 26a). Where does the mass ejected vertically from z = ±z_atm(R) go? Since the radial velocity in |z| > z_atm above the active zone is negative (the bottom panels of figures 23 and 25), the disk wind falls on to the center. In order to take into account mass accretion at high latitudes at R = R_buf, we measure the mass transfer rate |$(\dot{M}_\mathrm{ac})_\mathrm{buf}$| through the sphere r = R_buf, where we take into account only the regions where ρ > 10⁻², to remove the funnel regions. Figure 27 shows that |$(\dot{M}_\mathrm{ac})_\mathrm{buf}$| is comparable to |$(\dot{M}_\mathrm{ac})_z$| in all the times for LowRes and HighRes runs. This indicates that most mass ejected vertically falls on to the center (Takasao et al. 2018).

4.2.2 The transition zone

Figure 26c shows that the mass transfer rate |$(\dot{M}_\mathrm{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 coherent zone to the transition zone by the quasi-steady surface accretion flow.

Figure 26c shows that the vertical mass supply |$(\dot{M}_\mathrm{tr})_z$| stays almost zero. This does not mean that there is 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 = ±z_atm from R to R + dR that is given by

where the gas is ejected from the surfaces of z = ±z_atm when |$d\dot{M}_\mathrm{wind}/d\ln R \\gt 0$|⁠. Figure 28 shows that the mass is lost in R_MRI ≤ 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 of the gas outflowing from the inner part falls back within the transition zone owing to the failed disk wind (figure 13f).

4.2.3 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 26d). 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 |$(\dot{M}_\mathrm{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 |$(\dot{M}_\mathrm{co})_z\sim (\dot{M}_\mathrm{co})_{R-}$| (figure 26d).

Ferreira (1997) defines “the ejection index” as

which indicates a mass ejection efficiency because |$\dot{M}_\mathrm{wind}$| obtained from the radial integration of |$d\dot{M}_\mathrm{wind}/d\ln R$| increases more rapidly compared with the mass transfer rate |$\dot{M}^{z_{\mathrm{atm}}}$| as the radial extent considered in the integration for larger ξ.

Figure 28b shows that |$\dot{M}^{z_{\mathrm{atm}}} \sim d\dot{M}_\mathrm{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 λ ≡ (R_A|$/$|R₀)² are related, where R₀ is the radius at a radius where the disk wind is launched and R_A is the Alfvén radius of the radius, or ξ = [2(λ − 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 = z_atm, the plasma beta is β₀ exp [−(z_atm|$/$|H)²] ∼ 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.

4.3 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 = r_in 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 = r_in to R penetrate through the northern hemisphere at r.

From equation (63), one obtains

where R_out = r_out is the outermost cylindrical radius of the simulation box. Equation (64) shows how the total magnetic flux Ψ_rad(r_out) is distributed between Ψ_rad(r_in) and Ψ_mid(R_out).

The time evolution of Ψ_rad(r_in) and Ψ_mid(R_out) is shown in figure 29. Because Ψ_mid(R_out) is much larger than Ψ_rad(r_in) at the initial condition, figure 29 displays the deviations from the respective initial values, δΨ_rad(r_in) and δΨ_mid(R_out). In addition, their values are normalized by the initial mid-plane magnetic flux within R = R_MRI. Figure 29 shows that Ψ_rad(r_in) increases with time while Ψ_mid(R_out) decreases with time, keeping Ψ_rad(r_in) + Ψ_mid(R_out) almost constant.³ This clearly shows that the magnetic fluxes passing through the mid-plane accrete on to the inner edge. Comparison between the results of LowRes and HighRes shows that the magnetic flux accretes on to the inner boundary more rapidly for the HighRes than for the LowRes run.

The increasing rate of Ψ_rad(r_in) decreases with time except for a sudden increase around |$t=1000\, t_{\mathrm{K0}}$|⁠, which corresponds to the accretion of the innermost 〈B_z〉 concentration on to the inner boundary of the simulation box around |$t=1000\, t_{\mathrm{K0}}$| (figure 19). In the late evolution (⁠|$1500\, t_{\mathrm{K0}}\le t\le 2000\, t_{\mathrm{K0}}$|⁠), the magnetic flux at the inner boundary hardly changes while the gas accretion on to the inner boundary continues (figure 27). This decrease in the increasing rate of Ψ_rad(r_in) has been reported in Beckwith, Hawley, and Krolik (2009) and Takasao et al. (2018).

In sub-subsections 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 the drift speed of the magnetic flux v_B as follows:

(Guilet & Ogilvie 2012; Bai & Stone 2017).

4.3.1 The active zone

In this section, we will show that the magnetic flux drifts inward (outward) in inner (outer) regions of the active zone in the early phase evolution. In the later phase, the drift velocities oscillate between positive and negative with time in the inner active zone while the magnetic flux drifts inward in the outer active zone.

The trajectories of the magnetic flux and the color map of 〈B_z〉_H in figure 19a show that the flux concentrations become prominent for |$t\gtrsim 500\, t_{\mathrm{K0}}$|⁠, and the contours of Ψ_rad converges to the multiple 〈B_z〉 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 appendix 1). In each of the inner and outer regions, we average 〈v_B〉 in the radial direction as follows:

In order to investigate what determines the drift speed, we decompose 〈(E_ϕ)_I〉_H = 〈v_RB_z〉_H − 〈v_zB_R〉_H into the laminar component

and the turbulent component

where 〈(E_ϕ)_I〉_H = [〈(E_ϕ)_I〉_H]_lami + [〈(E_ϕ)_I〉_H]_turb is satisfied. Although the mass-weighted average is taken when we calculate the mean values of the velocity in subsection 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 〈v_{R, 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 |$(\widetilde{v_\mathrm{B}})_\mathrm{lami}$| and |$(\widetilde{v_\mathrm{B}})_\mathrm{turb}$|⁠, respectively.

In the early evolution (⁠|$t\le 500\, t_{\mathrm{K0}}$|⁠), 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 |$\widetilde{v_\mathrm{B}}$| 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⁻⁵v_K. Although the resolution-dependence of |$\widetilde{v_\mathrm{B}}$| is found in |$t\le 300\, t_{\mathrm{K0}}$|⁠, it disappears in 300 ≲ t|$/$|t_K0 ≲ 500.

The magnetic flux transport is not only determined by the dynamics inside the disk with |z| < H, but is also 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 on to 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(r_in) increases.

In the long-term evolution in the active zone, the turbulent and laminar contributions nearly cancel each other out and their summation, or |$\widetilde{v_\mathrm{B}}$| 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\ge 1000\, t_{\mathrm{K0}}$|⁠, indicating that there is almost no net magnetic transport. Such a quasi-steady 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, the phases of which are correlated with those of 〈B_ϕ〉 oscillations in the upper atmospheres. Unlike the LowRes run, the 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\sim\! 1500\, t_{\mathrm{K0}}$|⁠, beyond which the quasi-periodic MRI activities disappear (figure 37). Although |$\widetilde{v_\mathrm{B}}$| is positive on average in the early phase, the net drift velocities become around zero in 500 ≤ t|$/$|t_K0 ≤ 1000. In addition, unlike in the inner active zone, |$(\widetilde{v_\mathrm{B}})_\mathrm{lami}$| and |$(\widetilde{v_\mathrm{B}})_\mathrm{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\gtrsim 1500\, t_{\mathrm{K0}}$|⁠, |$(\widetilde{v_\mathrm{B}})_\mathrm{lami}$| and |$(\widetilde{v_\mathrm{B}})_\mathrm{turb}$| do not completely cancel each other out, and the magnetic flux drifts inward at |$\widetilde{v_\mathrm{B}}\sim\! 10^{-4}v_\mathrm{K}$|⁠. The magnetic flux drifts inward at a faster speed than the gas that accretes at a speed of ∼10⁻⁶v_K. Jacquemin-Ide, Lesur, and Ferreira (2021) conducted a simulation (model SEp) similar to ours and obtained a drift speed consistent with our results.

4.3.2 The transition zone

In sub-subsection 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 〈B_z〉_H ∼ 0 expands with time because the magnetic flux in the coherent zone drifts outward, as will be shown in sub-subsection 4.3.3. A typical drift speed of the outer edge that is identified at the radius where 〈B_z〉_H = 〈B_z,ini〉_H is estimated to about 10⁻³ v_K 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 〈v_B〉_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 〈B_z〉_H ∼ 0. However, the inward expansion would not proceed further because of rapid decreases in η_O and η_A (figure 2).

4.3.3 The coherent zone

Figure 33 shows the time evolution of the radial profiles of 〈v_B〉_H|$/$|v_K. For R ≲ 3.5 au, although 〈v_B〉_H|$/$|v_K is affected by the quasi-periodic disturbances from the active zone, the drift velocity is positive on average. Beyond R ≳ 3.5 au, 〈v_B〉_H is steady. When the time average over |$1000\, t_{\mathrm{K0}}\le t\le 2000\, t_{\mathrm{K0}}$| is taken, 〈v_B〉_{H, t} becomes almost constant, and is about 〈v_B〉_{H, t} ∼ 10⁻³v_K. This value is comparable to those in Bai and Stone (2017) and Bai (2017). Figure 34a shows that the drift velocity v_B = 〈E_ϕ〉_t|$/$|〈B_z〉_t does not depend on z significantly, indicating that the magnetic fields drift outward keeping their shapes. This feature is also consistent with Bai and 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 ∼ 〈v_RB_z〉_t is negative. Diffusion due to AD, which plays an important role in the surface accretion layers (Λ_A ∼ 1), tries to cause the magnetic flux to drift 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]|$/$|〈B_z〉 is comparable to the drift velocity in the region with Λ_O > 1.

5 Discussion

5.1 Structure formation

5.1.1 The active zone

As shown in figure 11, prominent multiple rings and flux concentrations are found in the active zone. They are not transient structures but are maintained at least until t|$/$|t_K0 = 2000. It should be noted that not all of the simulations conducted in previous studies have shown ring formations; 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 sub-subsection 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 understand the ring formation in the ideal MHD fully. 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.

5.1.2 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 sub-subsection 3.2.3, if a finite B_z exists in the initial condition, AD amplifies 〈B_R〉, 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 = R_MRI.

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. 2018, 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 to 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 (|B_R| ≫ |B_z|) at the mid-plane where the toroidal magnetic field is reversed. Then, magnetic reconnection is triggered, resulting in magnetic flux transport. Cui and 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_ϕ ≫ B_R, B_z. Béthune et al. (2017) pointed out that AD works as anti-diffusion. This behavior is also confirmed in our simulations in figures 21 and 34.

Although our simulations show the B_z 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². 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 concentrations. 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 ∼eight rotations at R = 4 au.

5.1.3 Vertical shear instability

Since the locally isothermal equation of state is used in our simulations (subsection 2.1), VSI may occur (e.g., Urpin & Brandenburg 1998; Nelson et al. 2013). In this section, we will show that our simulation time is too short to see non-linear growth of VSI in the transition and coherent zones.

Vertically global linear analyses show that the VSI modes can be divided into the two body modes (corrugation and breathing modes) that have large amplitudes near the mid-plane and the surface modes that emerge when the boundaries are set in the vertical direction (Nelson et al. 2013; Barker & Latter 2015; Umurhan et al. 2016). In the nonlinear stage, the fundamental corrugation mode dominates over the other modes (e.g., Nelson et al. 2013).

In all the regions, VSI does not appear to grow in our simulations. In the active zone, VSI is suppressed by magnetic fields (Latter & Papaloizou 2018; Cui & Lin 2021; Latter & Kunz 2022). Also in the transition and coherent zones, our simulations do not show the non-linear development of the corrugation modes.

To estimate the growth rate of VSI in the transition and coherent zones, we used the reduced model proposed by Nelson et al. (2013) with T ∝ R^−1/2 (also see Barker & Latter 2015; Umurhan et al. 2016). A significant difference from the previous studies (Nelson et al. 2013; Barker & Latter 2015; Umurhan et al. 2016) lies in the vertical extent where VSI occurs. In our simulations, for |z| ≳ H, Λ_O is greater than unity and the gas is marginally coupled with the magnetic field (figure 34c). It is reasonable to set the boundaries at |z| = H because the VSI is suppressed by the Lorentz force in |z| ≳ H. In addition, in the coherent zone, the surface gas accretion driven just above the OR dead-zone (|z| ∼ H) is likely to suppress the VSI (figure 34b). By contrast, the previous studies consider a larger vertical range (z > 5c_s|$/$|Ω). Unlike the previous studies (see figure 18 in Nelson et al. 2013), the limited vertical range in our setup (|z| ≤ H) makes the fundamental corrugation mode the most unstable of the body modes and reduces the growth rates of the higher-order body modes. This may be consistent with the fact that the velocity dispersion driven by the VSI is strongly suppressed when the vertical extent of the VSI unstable region is limited to |z| ≲ 2c_s|$/$|Ω (Fukuhara et al. 2023). The growth rate σ of the fundamental corrugation mode is insensitive to the radial wavelength when the radial wavelength is larger than O(ϵ²R). The growth rate is about ∼0.1 ϵΩ. The ratio of the growth timescale to the simulation time is about ∼σ⁻¹|$/$|t_end ∼ 0.5(R|$/$|1 au)^1.25 (table 1). This indicates that the simulation time is not long enough for the VSI to grow sufficiently nonlinear.

The surface mode, which develops near the vertical boundaries, is suppressed due to the limited vertical extent (|z| ≤ H) when the radial wavelength is less than ∼40(ϵ²R)⁻¹ in the reduced model. Although surface modes would emerge for higher k_R at higher growth rates, the radial resolution of our simulations is not sufficient to resolve them.

5.2 Implications for long-term evolution

Although we conducted the LowRes run until |$t=2000\, t_{\mathrm{K0}}= 40\:\mbox{yr}$|⁠, this is much shorter than a typical disk lifetime. In this section, we discuss the long-term evolution of the disk inferred from our results.

5.2.1 Gas evolution

First, we consider long-term evolution by assuming the magnetic flux is fixed. If the gas accretion rate shown in figure 27 is maintained (⁠|$\dot{M}\sim\! 5\times 10^{-4}$| in the simulation unit corresponding to |$\dot{M}=2.8\times 10^{-7}M_\odot \:\mbox{yr}^{-1}$|⁠), the depletion time-scale of the active zone is ∼4 × 10⁵ yr if there is no mass supply. This may be a formation mechanism of transition disks (e.g., van der Marel 2023).

What would the mass transport be between the transition and active zones? As discussed in sub-subsection 3.2.2 and subsection 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 |$\dot{M} = 10^{-3}$| in the simulation units (figure 25), a mass corresponding to the initial mass of the active zone is transferred to the transition zone in 2 × 10⁵ yr. The accumulated gas is not directly transferred to the active zone since gas accretion inside the transition zone does not occur because there is no net magnetic flux. Thus, the transferred mass forms a ring structure that has a surface density contrast of about 10³ if the radial width is ∼0.2 au at 2 × 10⁵ yr. In reality, the ring will become vortices via the Rossby-wave instability (Lovelace et al. 1999; Ono et al. 2016, 2018).

Gas accretion from the inner edge of the transition zone to the active zone may increase R_MRI when the surface density becomes low enough for cosmic rays and/or ionizing radiation to reach the mid-plane. If R_MRI gets close enough to R_η, the gas accumulated at the outer edge of the transition zone may be transferred to the active zone.

5.2.2 Magnetic flux evolution

The existence of the transition zone significantly affects 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 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 > R_MRI. If this is the case, the magnetic flux in R ≤ R_MRI decreases with time, leading to a decrease in the gas accretion rates on to the central star. As the process regulates the magnetic flux accretion on to the central star, it can be important to determine the magnetization of young stellar objects.

The outer edge of the region where B_z ∼ 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 B_z ∼ 0. If v_B 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 B_z ∼ 0, the inner edge of the coherent zone where the surface gas accretion occurs moves outward. If v_B 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 differences of about an order of magnitude 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.

5.3 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) models 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 3D 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, Flock, and Okuzumi (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⁻¹ (Blum & Wurm 2000; Wada et al. 2013). The dead zone in Ueda, Flock, and Okuzumi (2019) corresponds to the transition zone in our simulations. We note that Ueda, Flock, and Okuzumi (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 (sub-subsection 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 δv_tr does not satisfy |$\delta v_\mathrm{tr}/c_\mathrm{s} \sim \sqrt{\alpha _\mathrm{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 (δv_tr|$/$|c_s)². In terms of the velocity dispersion δv_tr, 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⁻¹.

Our results show that two rings form at both the edges of the transition zone, or R = R_MRI and R_η, although the setting of Ueda, Flock, and Okuzumi (2019) yields a single ring forming 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 number of the velocity dispersion |$\langle \delta {\cal M}_\mathrm{tot} \rangle$| 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⁻¹. 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⁴ would be preferred.

5.4 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 R_MRI is roughly determined by T ∼ 10³ K (Desch & Turner 2015), R_MRI is predicted to be

where κ_R is the Rossland mean opacity, and κ_R = 5 cm² g⁻¹ corresponds to the opacity of dust grains with |$0.1\, \mu$|m and a dust-to-gas mass ratio of 10⁻². In the 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.

Jankovic et al. (2021, 2022) investigated steady-state axisymmetric disk structures considering detailed physical and chemical processes. They found that the temperature structure is more complex than our adopted one [equation (10)]. This will affect the shape of the dead-zone boundaries. They derived a similar R_MRI as equation (71) even when they took into account the disk structures modified by convective instability.

Mori, Bai, and Okuzumi (2019) pointed out that R_MRI 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 the 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 subsection 5.2. Although criteria to form the transition zone are unclear, if the active zone exists stably, sub-subsection 3.2.3 shows that the formation of the transition zone requires η_O to increase rapidly across R ∼ R_MRI and η_A to be large enough to amplify B_R 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.

6 Summary

In this paper, we performed global non-ideal MHD simulations of the 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, R_MRI, R_η, and R_ch, are defined from the spatial distribution of η_O and η_A. The first radius R_MRI = 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 R_MRI < 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 < R_MRI) 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 (sub-subsections 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 (sub-subsection 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⁻⁴ v_K.

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 (subsections 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 (sub-subsection 4.3.3). Our simulations show that the latter dominates over the former, and the magnetic flux drift outwards at a speed of ∼10⁻³v_K.

3. The transition zone appears in R_MRI ≤ R ≤ R_η. Although it belongs to the conventional dead zone defined by Elsässer numbers, the physical properties of the transition zone 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 (sub-subsection 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 via 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 edges of the transition zone. The outflow launched around the inner edge of the transition zone fails to escape the disk but falls on to 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 a wider parameter search to investigate such structures systematically.

There are two main caveats in our study. One is that the Hall effect is neglected for simplicity. In our model setup, the diffusion coefficient of the Hall effect is as large as η_O and η_A near the dead-zone inner boundaries. We will include the Hall effect because it can change the disk structure significantly (Bai & Stone 2017; Béthune et al. 2017; Martel & Lesur 2022). The other is that we use the locally isothermal equation of state with a simple temperature profile [equation (10)]. The temperature profile is more complicated near the dead-zone inner boundaries when the radiative transfer is considered (e.g., Flock et al. 2017; Gressel et al. 2020). This can change the spatial distributions of the diffusion coefficients and the thermal evolution of the gas.

Acknowledgements

Numerical computations were carried out on supercomputer Fugaku provided by the RIKEN Center for Computational Science (Project ids: hp210164, hp220173), Oakforest-PACS at Joint Center for Advanced High Performance Computing (Project ids: hp190088, hp200046, hp210113), and Cray XC50 at the CfCA of the National Astronomical Observatory of Japan.

Funding

This work was supported in part by the Ministry of Education, Culture, Sports, Science and Technology (MEXT), Grants-in-Aid for Scientific Research, JP21H00056 (K.I.), JP16H05998 (K.T. and K.I.), JP21H04487, JP22KK0043 (K.T., K.I., and S.T.), JP22K14074 (S.T.), JP20H00182, JP23H00143, JP23K25923, JP23H01227 (S.O.), JP17H01105, JP21H00033, and JP22H01263 (T.K.S.). This research was also supported by MEXT as “Exploratory Challenge on Post-K computer” (Elucidation of the Birth of Exoplanets [Second Earth] and the Environmental Variations of Planets in the Solar System) and “Program for Promoting Researches on the Supercomputer Fugaku” (Toward a unified view of the universe: from large scale structures to planets, JPMXP1020200109).

Appendix 1. Resolution-dependence of the α parameter

We investigate whether the resolution of our simulations is sufficient to resolve the MRI. In our simulations, the scale height H(R) is resolved with N = 17(R|$/$|0.3 au)^1/4 for the LowRes run and N = 34(R|$/$|0.3 au)^1/4 for the HighRes run.

Noble, Krolik, and Hawley (2010) introduce the so-called quality factors, which are defined as the number of grid cells resolving the fastest MRI growing mode measured by amplified magnetic fields:

where |$v_{\mathrm{A},r,\theta ,\phi } = B_{r,\theta ,\phi }/\sqrt{4\pi \rho }$|⁠. Sorathia et al. (2012) show that in order to obtain the converged result Q_θ needs to be larger than 10–15 if Q_ϕ ∼ 10, and Q_θ can be smaller if Q_ϕ > 25 (see also Hawley et al. 2013).

Figure 35 shows the radial profiles of Q_r, Q_θ, and Q_ϕ averaged over 0 ≤ ϕ ≤ 2π and |z| < H both for LowRes and HighRes runs. The HighRes run barely satisfies the condition derived by Sorathia et al. (2012). By contrast, for the LowRes run, the fastest growing mode is not resolved in our simulation, although the MRI at longer wavelengths can be captured because Q_ϕ is around 10 for r > 0.2 au and Q_θ is around 3.

Figure 36 compares the time evolution of 〈α_Rϕ〉_H [equation (36)] between LowRes and HighRes runs. Both models show that the MRI turbulence is saturated after t ∼ 20t_K(R = 0.3 au). The saturation level of 〈α_Rϕ〉_H is about twice larger for HighRes than for LowRes runs, consistent with the fact that the LowRes run does not resolve the MRI turbulence sufficiently. The saturation level 〈α_Rϕ〉_H ∼ 2.8 × 10⁻² for the HighRes run is comparable to those found in local shearing-box simulations (e.g., Hawley et al. 1995; Sano et al. 2004). We confirmed that the HighRes run is expected to provide converged results by measuring the so-called quality factors (Noble et al. 2010).

Appendix 2. MRI turbulence in the active zone

Evolution of the MRI-driven turbulence is characterized by magnetic-field amplification owing to the channel flows followed by dissipation (Sano & Inutsuka 2001). The time evolution of the vertical structures at R = 0.3 au is shown in figure 37.

Figure 37a shows that the vertical structure of 〈B_ϕ〉 changes around |z| ∼ 2H (Suzuki & Inutsuka 2009). Around the mid-plane for |z| < 2H, the gas pressure dominates over the magnetic pressure, and the magnetic fields are amplified by the MRI dynamo activity. At |z| ∼ 2H, the magnetic pressure becomes comparable to the gas pressure (figure 37b), and the maximum growth scale of the MRI is comparable to H since it is expressed as 2πv_A|$/$|Ω₀. As a result, the magnetic fields are amplified by the channel flow that has a scale comparable to H, and are dissipated by magnetic reconnection. The t–z diagram of B_ϕ shows a clear so-called butterfly structure where B_ϕ flips its sign quasi-periodically and drifts toward upper atmospheres. The timescale of the quasi-periodic activity is about a few tens of rotations, which is consistent with local shearing-box simulations with stratified density structures.

Figure 37c shows |$\sqrt{\langle B_\phi \rangle ^2/\langle B_\phi ^2 \rangle _\phi }$|⁠, which is an indicator of the contribution of the coherent component of B_ϕ to the total toroidal field strength. The coherent field component dominates in the upper atmosphere for |z| > 2H while the field structure is controlled by the MRI turbulence in the disk (|z| < 2H).

The gas motion is affected by the magnetic field structures. As shown in figure 37d, R〈ρv_R〉 stochastically changes owing to the MRI turbulence for |z| < 2H, and it has a coherent quasi-periodic structure in the upper atmosphere. The time-variability of R〈ρv_R〉 is correlated with that of 〈B_ϕ〉. This suggests that the magnetic torque of the coherent magnetic field drifting upwards for |z| > 2H drives the coherent gas motion near the surface layers.

In the long-term evolution, figure 37 shows that the MRI dynamo at R = 0.3 au is somewhat weakened and the plasma β starts to increase when t ≳ 1300t_K0 corresponding to 250t_K(R = 0.3 au). This is related to the formation of ring structures shown in sub-subsection 3.2.1.

Footnotes

References