A Tale of Two Grains: Impact of Grain Size on Ring Formation via Nonideal Magnetohydrodynamic Processes

Studying the dynamics of protoplanetary disks (PPDs) is crucial for not only constructing comprehensive planet formation theory, but also understanding the fine features in recent high-resolution observations. The 2014 Atacama Large Millimeter/submillimeter Array ALMA campaign provided unprecedented details of the HL Tau disk (ALMA Partnership et al. 2015). Patterns that are broadly axisymmetric, such as multiple bright and dark rings in dust continuum emission, have been found to be common in protoplanetary disks (PPDs) in subsequent observations (e.g., Huang et al. 2018; Long et al. 2018).

Recent magnetohydrodynamic (MHD) simulations suggest that substructure formation can be achieved by nonideal MHD processes, including ohmic resistivity and ambipolar diffusion. Processes like redistribution of magnetic flux, direct feeding of avalanche accretion, and midplane magnetic reconnection (Suriano et al. 2017, 2018, 2019) can alter disk surface density. The ratio between gas surface density and vertical magnetic strength affects the local accretion rate. Regions with higher field-to-density ratios tend to drive stronger winds and accrete faster producing gaps; those with lower field-to-density ratios tend to accrete more slowly, allowing matter to accumulate and form dense rings (Suriano et al. 2017). Similar instability associated with disk wind has also been reported by Riols & Lesur (2019). With magnetic diffusion, the initial smooth magnetic flux is redistributed, causing radial density variation. The abrupt change in azimuthal velocity at disk surface twists vertical magnetic field to the extent that it efficiently brakes the accretion stream as it loses angular momentum. This is a rapid, runaway-like process similar to an avalanche (Kudoh et al. 1998; Suriano et al. 2017). The accretion stream feeds the local disk directly, hence, building a radially denser region, i.e., a ring. Magnetic reconnection is one of the major processes that form substructures observed in MHD disks with moderate ambipolar diffusion (Suriano et al. 2018; Hu et al. 2019). The accretion flow concentrated at the midplane drags field lines inward, pinching the magnetic field radially. Gradually, the highly pinched field lines, together with field lines dragged by infalling gas toward the midplane, eventually reconnect. The magnetic stress drives gas to the opposite direction at each side of the reconnected field loop, draining the material inside to form a gap.

Meanwhile, snow lines modify dust drift, growth, and fragmentation, causing surface density variation radially (Zhang et al. 2015; Okuzumi et al. 2016). The abundance of dust is an important factor in a disk's ionization, affecting the coupling between gas and magnetic fields (e.g., Wardle 2007; Bai 2011a). Near snow lines, we can expect the feedback from dust to nonideal MHD, then to disk accretion itself. In Hu et al. (2019), this idea was tested by incorporating snow-line-induced dust distribution and ionization change into MHD disk simulations. The ionization structure was precalculated with dust distribution from global dust evolution, including sintering effects (Okuzumi et al. 2016). The ionization structure and the nonideal MHD diffusivities undergo sharp changes at the snow lines due to the change in dust size there. This leads to a discontinuous accretion flow across the snow lines. With time, such a discontinuous accretion flow naturally produces gaps and adjacent rings at the CO₂ and C₂H₆ snow lines. This dust-to-gas feedback requires much less dust to take effect (∼10⁻² mass fraction) than a hydrodynamic drag mechanism like streaming instability (e.g., Gonzalez et al. 2017), where a dust-to-gas ratio at the order of unity is usually needed where a disk structure has formed. However, using only a spatially dependent and time-invariant magnetic diffusion profile (e.g., Hu et al. 2019) is also one of the major caveats of current nonideal MHD simulations. Grain-modified magnetic diffusivities in PPDs were first explored in the very small grain limit (approximately nanometer size, or on the polycyclic aromatic hydrocarbon (PAH) scale). They are known for substantially reducing the level of ionization, killing magnetorotational instability (Bai 2011a; Wang et al. 2019, hereafter WBG19. Charged grains, when very small, can be regarded as heavy ions in diffusion calculations, while its geometric size starts to take effect in larger (submicron- to micron-sized) grains. Thus, there is no simple functions to describe the relation between magnetic diffusivities and magnetic field strength when dust is included in disks, and the exact dependence also varies over different grain sizes. This is systematically explored in Xu & Bai (2016).

Another major caveat of most global MHD simulations is the uncertain disk temperature at the disk surface. Usually, a fixed or fast relaxation temperature profile is employed, and the initial density structure is derived based on hydrostatic equilibrium in the R–z plane. The transition between the hot atmosphere (also referred to as corona) and the cold, dense disk is either sharp (e.g., Suriano et al. 2017, 2018, 2019) or with some artificial smooth profile (e.g., Bai & Stone 2017; Béthune et al. 2017; Hu et al. 2019).

In this work, we aim to overcome these two shortcomings by utilizing a nonideal MHD simulation with self-consistent thermochemistry and radiative transfer (Wang & Goodman 2017; WBG19). The dust structure from global dust evolution is input into the simulations, while ohmic resistivity and an ambipolar diffusion profile can be obtained in real time based on the ionization structure from the thermochemistry calculations. This paper is organized as follows. Section 2 briefly describes the global dust evolution model and discusses the different roles played by small and large grains in magnetic diffusion. Section 3 summarizes the numerical methods and parameter choices for our simulations. Section 4 presents diagnostics of disk radial and vertical structures, focusing on relating chemistry to ionization and at last accretion. We summarize the paper in Section 5.

Similar to Hu et al. (2019), we use the Okuzumi & Tazaki (2019) radial profiles of dust surface density and particle size derived from a one-dimensional global dust evolution model having multiple snow lines. Putting a 1D dust profile into a global spherical disk simulation may look inconsistent, while in reality, the majority of dust grains settle down to the midplane, where nonideal MHD processes are most significant. The disk atmosphere is much more ionized due to its low density, so the presence of dust does not change its physical nature in any significant way. Thus, a 1D dust calculation is physically adequate to describe the radial dust structure close to the midplane.

The evolution of dust surface density Σ_d is described by

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{\rm{\Sigma }}}_{d}}{\partial R}+\displaystyle \frac{1}{R}\displaystyle \frac{\partial }{\partial R}\left({{Rv}}_{R}{{\rm{\Sigma }}}_{d}\right)=0.\end{eqnarray} \tag{ 1 }$$

Here, R is the distance to central star in cylindrical coordinates, v_R is the dust radial velocity caused by aerodynamic drag from gas. v_R at each location is calculated by selecting a representing size considering coagulation and fragmentation. The details of the single-sized approach is described in Sato et al. (2016). The calculation takes into account the low stickiness of CO₂ ice (Okuzumi & Tazaki 2019) and aggregate sintering (Okuzumi et al. 2016). In this 1D dust evolutionary model, the fragmentation threshold velocity in regions where aggregate sintering takes place is chosen to be 40% of the threshold for unsintered aggregates. We assume a disk of weak gas turbulence with a velocity dispersion of 1.7% of the sound of speed. The distribution of the gas surface density for the dust evolution calculation is fixed in time and assumed to scale as R⁻¹, where R is the distance from the central star, at R < 150 au. We evolve the dust disk 1D model for 1.3 Myr, until the computed total millimeter fluxes from the dust disk match the ALMA observations of the HL Tau disk (ALMA Partnership et al. 2015). To get a correct ionization, a range of grain sizes is reconstructed with the dust evolution calculation, which follows a power law with minimum and maximum grain sizes of ${a}_{\mathrm{gr},\min }$ and ${a}_{\mathrm{gr},\max }$, respectively, where ${a}_{\mathrm{gr},\max }$ is the representative size obtained from the dust evolution. The power law is such that the vertically integrated number density of grains per unit grain radius a_gr is proportional to ${a}_{\mathrm{gr}}^{-3.5}$. We fixed the minimum grain size ${a}_{\mathrm{gr},\min }$ as 0.1 μm.

This work aims to determine whether a given dust profile can change the gas structure by the means of a nonideal MHD process. This feedback via nonideal MHD can only be modeled in the global simulation, not in the 1D dust evolution calculation. To avoid unnecessary complexity and isolate nonideal MHD feedback from other physical processes, we chose the simplest gas structure, "the smooth disk," in the 1D dust evolutionary model. Generating similar dust structures to match observations is not the main purpose of this paper. The idea is that if we believe these dust structures follow snow lines, how will they change the gas structure when nonideal MHD is considered? This is partially answered in Hu et al. (2019), where a similar 1D dust evolution calculation is implemented in a 3D global nonideal MHD simulation. In this work, we study the nonideal MHD feedback in a more detailed and consistent way. It is the first time the diffusion coefficients are calculated with real-time temperature and magnetic fields, provided with different grain sizes. We realize that it is fully not self-consistent since the gas structure will subsequently change the dust structure, e.g., dust trapping at the pressure bump. Whether this feedback is positive or negative will be addressed in our next work, where we will solve dust dynamics simultaneously during the MHD simulation. However, we think understanding individual physical process (here, the dust's effect on the gas) is important for our ultimate understanding of dust and nonideal MHD processes, especially since such an aspect has not been studied before.

The profile of Σ_d is plotted in the upper panel of Figure 1 (also shown in the left, second bottom panel in Figure 3 of Okuzumi & Tazaki 2019). Because of sintering, dust particles behind the snow lines CO₂, C₂H₆, and CH₄ (located at R = 13, 29, and 82 au, respectively) experience enhanced collisional fragmentation and pile up there. Collisional fragmentation is also enhanced interior to the H₂O snow line (located at 4 au) due to lack of water ice and results in another traffic jam of dust in this region.

Zoom In Zoom Out Reset image size

The adopted dust's total geometric cross section for chemistry calculation is plotted in the lower panel of Figure 1. Multiplying the average cross section of grains of all sizes by the total dust abundance (number density over hydrogen nucleus), we get the dust's total geometric cross section per hydrogen nucleus (σ_gr/H hereafter). The roles that dusts play in our thermochemical networks can be categorized as follows: (1) absorb irradiative photons and yield free electrons via a photoelectric effect for h ν > W photons (W ≃ 4.4 eV is the surface work of carbonaceous grains), (2) affect gas thermodynamics by dust-gas thermal accommodation (this is the key process that maintains the thermal structure of PPDs), and (3) catalyze reactions (especially neutralization of charged chemical species). All three types of processes, especially (3), which determines the MHD-dominated dispersal of disks, are included in our thermochemical calculations. These processes depend on the per hydrogen geometric cross section indirectly (1) and directly ((2) and (3)). To make all models thermochemically consistent with the profile calculated in Okuzumi & Tazaki (2019), we keep this number invariant in different models within this paper. Note that we enhanced σ_gr/H by 10 times when R < 4 au to ensure a moderate ionization fraction at inner disk.

In a weakly ionized PPD, the equation of motion for all charged species (e⁻, ions, and Gr^±, whose inertia is usually negligible, see, e.g., Xu & Bai 2016) is set by the equilibrium between the Lorentz force and the collision with neutrals,

$$\begin{eqnarray}&&{Z}_{j}e\left({\boldsymbol{E}}^{\prime} +\displaystyle \frac{{{\boldsymbol{v}}}_{j}}{c}\times {\boldsymbol{B}}\right)={\gamma }_{j}\rho {m}_{j}{{\boldsymbol{v}}}_{j},\end{eqnarray} \tag{ 2 }$$

where for the charged species indicated by j, Z_j is its charge number (in unit charge e), m_j is its mass, v _j is its drift velocity relative to the neutral background, γ_j ≡ 〈σ v〉/(m + m_j ) (m is the averaged particle mass of the neutrals), and 〈σ v_j 〉 is its rate coefficient of momentum transfer with neutrals. The electric field in the frame co-moving with the neutrals is denoted by ${\boldsymbol{E}}^{\prime}$.

In the case of very tiny grains, the interaction between charged grains and neutral molecules are mediated by the electric field from induced electric dipoles in the neutrals. This r⁻⁴ electric potential makes the collision rate coefficient 〈σ v〉_j independent of temperature (e.g., Draine 2011). When grains get larger, their geometric cross section dominates their interactions with neutrals. This transition brought by the grain size is reflected by the Hall parameter, which is the ratio between the gyrofrequency of charged species under Lorentz force and their collision frequency with the neutrals:

$$\begin{eqnarray}&&{\beta }_{j}=\displaystyle \frac{| {Z}_{j}| {eB}}{{m}_{j}c}\displaystyle \frac{1}{{\gamma }_{j}\rho }.\end{eqnarray} \tag{ 3 }$$

In our simulation, we adopt the recipes in Bai (2011b, 2014) to calculate collision coefficients between charged grains and neutrals:

$$\begin{eqnarray}&&\begin{array}{l}\langle \sigma v{\rangle }_{\mathrm{gr}}=\max \left[1.3\times {10}^{-9}| {Z}_{\mathrm{gr}}| ,\right.\\ \quad \left.4\times {10}^{-3}{\left(\displaystyle \frac{{a}_{\mathrm{gr}}}{1\mu {\rm{m}}}\right)}^{2}{\left(\displaystyle \frac{T}{100{\rm{K}}}\right)}^{1/2}\right]\,{\mathrm{cm}}^{3}\,{{\rm{s}}}^{-1}.\end{array}\end{eqnarray} \tag{ 4 }$$

So the transition from an electric-potential-dominated cross section to a geometric cross section is at approximately the nanometer scale. Given T = 100 K, any single charged grain with a_gr > 5.7 × 10⁻⁸ cm needs to consider the geometric effect when calculating collision coefficient. This difference is shown in the ohmic, Hall, and Pedersen conductivities:

$$\begin{eqnarray}\begin{array}{rcl}{\sigma }_{{\rm{O}}} & = & \displaystyle \frac{{ec}}{B}\sum _{j}{n}_{j}| {Z}_{j}| {\beta }_{j},\\ {\sigma }_{{\rm{H}}} & = & \displaystyle \frac{{ec}}{B}\sum _{j}\displaystyle \frac{{n}_{j}{Z}_{j}}{1+{\beta }_{j}^{2}},\\ {\sigma }_{{\rm{P}}} & = & \displaystyle \frac{{ec}}{B}\sum _{j}\displaystyle \frac{{n}_{j}| {Z}_{j}| {\beta }_{j}}{1+{\beta }_{j}^{2}}.\end{array}\end{eqnarray} \tag{ 5 }$$

Here, the summation index j runs through all charged species, with n_j indicating the number density and charge of individual charged species. Since β_gr, β_i ≪ 1 stands true throughout most regions of our disk, the difference between ions and grains is negligible for Hall conductivity σ_H. For ohmic and Pedersen conductivities, the contribution of larger grains can be orders of magnitude smaller than very small grains and ions. Using these conductivities, the general expressions for ohmic diffusivity η_O, Hall diffusivity η_H, and ambipolar diffusivity η_A are (Bai 2011b; Wang et al. 2019)

$$\begin{eqnarray}\begin{array}{rcl}{\eta }_{{\rm{O}}} & = & \displaystyle \frac{{c}^{2}}{4\pi }\displaystyle \frac{1}{{\sigma }_{{\rm{O}}}},\quad {\eta }_{{\rm{H}}}=\displaystyle \frac{{c}^{2}}{4\pi }\displaystyle \frac{{\sigma }_{{\rm{H}}}}{{\sigma }_{{\rm{H}}}^{2}+{\sigma }_{{\rm{P}}}^{2}},\\ {\eta }_{{\rm{A}}} & = & \displaystyle \frac{{c}^{2}}{4\pi }\displaystyle \frac{{\sigma }_{{\rm{P}}}}{{\sigma }_{{\rm{H}}}^{2}+{\sigma }_{{\rm{P}}}^{2}}-{\eta }_{{\rm{O}}},\end{array}\end{eqnarray} \tag{ 6 }$$

In general, when we deal with PPDs, the diffusivities η_O and η_A increase with larger grains due to aforementioned analyses.

The nonideal induction equation is

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t} & = & {\rm{\nabla }}\times \left({\boldsymbol{v}}\times {\boldsymbol{B}}\right)\\ & & -\displaystyle \frac{4\pi }{c}{\rm{\nabla }}\times \left({\eta }_{{\rm{O}}}{\boldsymbol{J}}+{\eta }_{{\rm{H}}}{\boldsymbol{J}}\times {\boldsymbol{b}}+{\eta }_{{\rm{A}}}{{\boldsymbol{J}}}_{\perp }\right),\end{array}\end{eqnarray} \tag{ 7 }$$

where v is the gas velocity, B is the magnetic field, b = B /∣B∣ is the unit vector representing the field line direction, J is the current density vector, and J _⊥ is the current component perpendicular to the magnetic field.

This paper studies PPDs that are subject to the impact of spatial distribution of dust grains. We simulate the disk evolution with the combination of nonideal MHD effects using the higher-order Godunov MHD code Athena++, ray-tracing radiative transfer for high energy photons, and consistent thermochemistry. For each MHD evolution timestep, the nonequilibrium thermochemistry is coevolved in each zone throughout the simulation domain using a GPU-accelerated the semi-implicit method. In general, this numerical system is the same as that in WBG19, and we refer the reader for WBG19 and references therein for the details in the initial conditions and the overall setup of thermochemical reactions. For the boundary conditions, we inherited the setups in WBG19, except for the toroidal field above the disk region (viz. inside the wind region) at the inner radial boundary: we set B_ϕ = − B_r ∣_t=0 (the initial value of r component) to suppress magnetic instabilities there (note that the magnetization here is ∼20 times that of WBG19 in terms of absolute magnetic pressures and stresses). The plasma β at the midplane is 2 × 10⁴, 5 times smaller than that of WGB19. Other hydrodynamic and field components are identical to those of WBG19. We use the same β value as Hu et al. (2019), which is based on the value in Hasegawa et al. (2017). Such a value at the disk midplane will be able to reproduce both the accretion rate and dust settling observed in the HL Tau disk under the framework of magnetically driven disk accretion.

For simplicity, the dust grains are still treated as single-sized carbonaceous grains co-moving with the gas. The size is either a_gr = 5 Å (Model S) or a_gr = 10⁻² μm (Model L). Since the thermochemical properties of dusts are mainly determined by σ_gr/H (Section 2), in order to clearly exhibit the impact of dust sizes on the magnetic diffusivities and disk dispersal, we keep the σ_gr/H profile matching the bottom panel of Figure 1 for both models, as a function of the cylindrical radius R to the central star. Note that the number density of grains is changed accordingly. The basic properties of our fiducial model (Model S) are summarized in Table 1. Maximum σ_gr/H is 8 × 10⁻²³ cm² in both models, which equals 7 × 10⁻⁶ in dust-to-gas mass ratio in Model S, and 1.4 × 10⁻⁴ in Model L, if we take the average density of a single grain as 2.25 g cm⁻³. We note that the grains in the two models have the same total geometric cross section but couple to magnetic fields differently. The grains in Model L have a grain-neutral collision coefficient 〈σ v〉_gr that is two orders of magnitude larger than the grains in Model S. The behavior of very small grains is actually similar to ions, as the Hall parameter of ions is approximately β_i ≈ 3.3 × 10⁻³(B_G/n₁₅). For larger charged grains, smaller β_gr indicates weaker coupling to magnetic fields, which means that grains are more prone to frequent collision with neutrals.

Table 1. Properties of Model S (Section 3)

Item	Value
Radial domain	1 au ≤ r ≤ 100 au
Latitudinal domain	0.06 ≤ θ ≤ 1.5708
Resolution	${N}_{\mathrm{log}r}=384$, Nθ = 128
Stellar mass	1.0 M⊙
Mdisk(1 au ≤ r ≤ 100 au)	0.10 M⊙
Initial midplane density	8.03 × 1014(R/au)−2.2218 mp cm−3
Initial midplane plasma β	104
Initial midplane temperature	305(R/au)−0.57 K
Artificial heating profile a	305(R/au)−0.57 K
Luminosities [photon s−1]
7 eV ("soft" far-UV)	4.5 × 1042
12 eV (Lyman–Werner)	1.6 × 1040
3 keV (X-ray)	1.1 × 1038
Initial abundances [nX/nH]
H2	0.5
He	0.1
H2O	1.8 × 10−4
CO	1.4 × 10−4
S	2.8 × 10−5
SiO	1.7 × 10−6
Dust/PAH properties
agr	5 Å
σgr/H	Variable (Section 3)

Note.

^a The artificial heating profile indicates the temperature of dust in the midplane. Similar to WBG19, because the radiative transfer of diffuse infrared radiation field is not calculated in this paper, we adopt this profile of dust temperature as the floor of the dust temperature.

Download table as:  ASCIITypeset image

Figure 2 summarizes the main point of this paper: the impact of grain sizes on radial ionization structure and how rings and gaps form from it. The top panel shows that the two models have significant different free electron density when σ_gr/H ≲ 10⁻²³ cm. As we kept σ_gr/H unchanged between the two models, Model L has a lower grain number density. Because we only considered single electron capture/loss of a grain, the process Gr + e⁻ − → Gr⁻ is saturated in Model L when σ_gr/H ≲ 10⁻²³ cm, leaving more free electrons around. The diffusion structure is plotted in dimensionless Elsässer numbers of ohmic resistivity (${\rm{\Lambda }}={v}_{a}^{2}/\left({\eta }_{{\rm{O}}}{{\rm{\Omega }}}_{K}\right)$) and ambipolar diffusion ($\mathrm{Am}={v}_{a}^{2}/\left({\eta }_{{\rm{A}}}{{\rm{\Omega }}}_{K}\right)$) at the disk midplane. The ambipolar diffusivity η_A is roughly proportional to B², so its Elsässer number Am is relatively constant when the field strength changes during simulation, and is a better way to quantify the diffusion strength. In Model S, where grains have the same size as those in WBG19 (5 Å PAH), the variation of dust number density yields some effects on surface density beyond 10 au, and there are no apparent features associated with σ_gr/H bumps at 13 and 28 au. In Model L, the location of rings and gaps at these two snow lines are very close to those of Hu et al. (2019), where the same dust evolution model was used. The ionization structure (Elsässer number) is much simpler in that work. It is precalculated, only spatially dependent and fixed in the MHD simulation. The similarity of substructure produced in the model with large grains and Hu et al. (2019) indicated the robustness of the nonideal MHD-based structure forming mechanism. From the nonideal MHD point of view, the simulations prove that number density variation of PAH-scaled grains has sublinear impacts on magnetic diffusivities. The two middle panels show the Elsässer numbers of both ohmic resistivity and ambipolar diffusion of Model S, where only mild changes are seen at the 13 and 28 au snow lines. The diffusivity profiles of Model L reflect the σ_gr/H structures quite well. Diffusion strength increases up to four orders of magnitude at those two snow lines, where the total σ_gr/H rises by about 10 times.

Zoom In Zoom Out Reset image size

After 1500 orbits at the innermost radius, the overall disk structure of both simulations reaches a quasi-steady phase, as shown in Figure 3. Quasi-steady means the overall structure does not change much during the last couple hundred orbits. The disk is revolved 10 orbits at 30 au, enough to develop and stabilize gaseous structures within the area of interest. In MHD simulations, there are a lot of physical processes occurring at different timescales and spatial scales. When we are mostly interested in the formation of gas density structures, 1D gas surface density is a good way to check if it is a quasi-steady state. Figure 3 also gives a rudimentary view on how different structures are created between the two setups. The leftmost panels indicate the number density of the summation of all positive charge carriers except Gr⁺ and are similar in Model L and Model S. Their distinctions start to show when it comes to magnetic diffusion. The ambipolar Elsässer number Am has a generally smooth profile in Model S, while in Model L it has sharp transitions at the location of the snow lines, where the grain number density changes dramatically. These are magnetic field responses to these changes. Both vertical and azimuthal components form structures: B_ϕ is nearly cut off inside 4 au, and B_z reversed its direction at the snow-line radii. This reflected in how the two disks accrete. In Model S, accretion ( v _R < 0) is almost everywhere near the midplane, while in Model L, there are several decretion ( v _R > 0) zones that expand all the way down to the midplane. Around the boundaries between accretion and decretion, materials are concentrated or depleted, forming rings and gaps.

Zoom In Zoom Out Reset image size

In what follows, we recognize the disk surface as the location of wind launching, i.e., where both v_z and v_R are positive, and poloidal velocity $\sqrt{{v}_{z}^{2}+{v}_{R}^{2}}$ is more than 50% of the local sound speed. Above the wind base, the Alfvén surface is where the velocity in the poloidal plane equals the poloidal Alfvén velocity ${v}_{{\rm{A}},p}=\sqrt{({v}_{z}^{2}+{v}_{R}^{2})/4\pi \rho }$.

4.1. Thermochemistry

4.2. Dynamics and Kinematics of Accretions and Winds

In a (quasi-)steady accretion disk, the accretion is driven by (1) the radial gradient of T_{R ϕ} (R ϕ component of the magnetic stress tensor), and (2) the difference in the T_{z ϕ} stress (z ϕ component of the magnetic stress tensor) between the top (marked as top for the integration limit) and bottom disk surfaces (marked as bottom for the integration limit), namely,

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{\dot{M}}_{\mathrm{acc}}{v}_{{\rm{K}}}}{4\pi } & = & \displaystyle \frac{\partial }{\partial R}\left({R}^{2}{\displaystyle \int }_{\mathrm{bottom}}^{\mathrm{top}}{T}_{R\phi }\ {dz}\right)\\ & & +{R}^{2}{\left[{T}_{z\phi }\right]}_{\mathrm{bottom}}^{\mathrm{top}}.\end{array}\end{eqnarray} \tag{ 8 }$$

Here, ${\dot{M}}_{\mathrm{acc}}$ is accretion rate within the disk and v_K is the Keplerian orbital velocity. The first term on the right-hand side of Equation (8) (${\dot{M}}_{\mathrm{acc}}({T}_{R}\phi )$ in Figures 4 and 5) resembles a radial stress and can be characterized by the equivalent Shakura–Sunyaev α parameter,

$$\begin{eqnarray}&&\alpha ={\left[{\int }_{\mathrm{bottom}}^{\mathrm{top}}{p}_{\mathrm{gas}}{dz}\right]}^{-1}\times {\int }_{\mathrm{bottom}}^{\mathrm{top}}{T}_{R\phi }\ {dz}.\end{eqnarray} \tag{ 9 }$$

Quantities associated with these components of accretion rates are summarized in Figures 4 and 5. The accretion rate ${\dot{M}}_{\mathrm{acc}}$ is the integrated mass flux below the wind base, i.e., the dotted–dashed gray lines in Figure 3. The wind loss rate ${\dot{M}}_{\mathrm{wind}}$ is vertical mass flux (ρ v_z ) measured at the wind base (the equivalent of the disk surface, marked as wb in equations):

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{acc}}(\mathrm{Measured})=2\pi R{\int }_{\mathrm{bottom}}^{\mathrm{top}}{\rho }_{\mathrm{gas}}{v}_{R}{dz}\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{wind}}(R)=4\pi {\int }_{{{\rm{R}}}_{\min }}^{{\rm{R}}}{\rho }_{\mathrm{gas},\mathrm{wb}}{v}_{z,\mathrm{wb}}R^{\prime} {dR}^{\prime} .\end{eqnarray} \tag{ 11 }$$

Similar to Hu et al. (2019) and WGB19, the T_{z ϕ} stress (viz. the stress that drives wind) is the predominant factor of accretion in both models.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The second term on the right-hand side of Equation (8) represents the contribution of wind stresses (${\dot{M}}_{\mathrm{acc}}({T}_{z}\phi )$ in Figures 4 and 5, also written as ${\dot{M}}_{\mathrm{acc},\mathrm{wind}}$). There is a certain ratio between the accretion rate caused by wind stress and the local wind loss through the disk surface $d{\dot{M}}_{\mathrm{wind}}/d\mathrm{ln}R$. Along a field line, the cylindrical radius of wind base R_wb and the radius where it intersects with the Alfvén surface R_A is related to the ratio of the wind loss rate and the steady accretion rate caused by wind torque:

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{acc},\mathrm{wind}}=2\left[{\left(\displaystyle \frac{{R}_{{\rm{A}}}}{{R}_{\mathrm{wb}}}\right)}^{2}-1\right]\left(\displaystyle \frac{d{\dot{M}}_{\mathrm{wind}}}{d\mathrm{ln}R}\right).\end{eqnarray} \tag{ 12 }$$

The ratio R_A/R_wb is referred to as the "magnetic lever arm." For example, based on the field lines presented in the lower panel of Figure 3, the value of the lever arm is about 1.4 within 10 au; hence, the local accretion rate should be 1.92 times $d{\dot{M}}_{\mathrm{wind}}/d\mathrm{ln}R$ in Model S.

4.2.1. Model S

For Model S, the lower panel of Figure 3 shows a steady accretion below ∼2/3 the height of the wind base, and outflow dominates above that. Both the velocity streamlines and magnetic field lines in Figure 3 show apparent wind launching inside 10 au, while beyond 10 au, the outgoing flux is almost parallel to the disk surface. The strengths of both the poloidal and toroidal field components are significantly weaker than R < 15 au. Therefore, the wind torque is weaker at the outer region, and a lower accretion rate and slower evolution of the substructures are expected.

Inside the disk, the Elsässer numbers Am and Λ have no apparent feature observed at the 13 and 28 au snow lines due to the effects elaborated in Section 4.1. The accretion rate peaks (${\dot{M}}_{\mathrm{acc}}\simeq 5\times {10}^{-6}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$) near R ∼ 7 au. Beyond that, the accretion rate decreases with radius and means that less material is being transported into this region than the amount of material that is accreted away from the inner radius. The net loss of the material is clearly reflected by the resulting surface density profile (Figure 2). Though T_{z ϕ} -driven accretion also drops beyond 7 au, the decrease in the measured accretion rate between 7 and 20 au is mainly caused by T_{R ϕ} , as the dashed line between 9 au ≲ R ≲ 22 au shows the radial gradient of the T_{R ϕ} stress is countering the direction of accretion: it tries to move material outward rather than inward. A more direct way to understand how material is redistributed in the disk is radial differentiation of accretion $d{\dot{M}}_{\mathrm{acc}}/{dR}$ and local wind loss rate $d{\dot{M}}_{\mathrm{wind}}/{dR}$. In the bottom panel of Figure 4, the local wind loss rate is always positive, meaning that there is no material infall from the atmosphere. $d{\dot{M}}_{\mathrm{acc}}/{dR}\gt 0$ within 8 au is almost an order of magnitude higher than the wind loss rate per unit radius. The disk gains mass at R ≲ 8 au. Between 8 and 20 au, negative $d{\dot{M}}_{\mathrm{acc}}/{dR}$ means both accretion and wind is taking material away, leaving a wide gap in the surface density profile, as shown in Figure 2. This indicates that the toroidal field decreases faster than R⁻² in this range. Within 10 au, the accretion rate is about twice as that of $d{\dot{M}}_{\mathrm{wind}}/d\mathrm{ln}R$, which is consistent with the magnetic lever arm argument (1.96 as calculated in Equation (12)). The total wind loss rate, which is the maximum cumulative wind loss rate _wind(cumulative) is shown in Figure 4), and is about the same as the maximum accretion (at ∼8 au). Overall, the rate of accretion and wind loss are at the same level in Model S.

4.2.2. Model L

The outcomes of our simulation are quite different in Model L. At radii slightly above the 4, 13, and 28 au snow lines (where σ_gr/H rises significantly), both Am and Λ decrease to yield stronger diffusion (see also Section 4.1). Such features in magnetic diffusion modulate the way the disk accretes. Figure 3 presents an obvious characteristic of Model L: the absence of wind and wind-driven accretion at radii R ≲ 4 au. Immediately outside this snow-line radius, the concentration of magnetic fields causes excessive mass transfer rates in both wind and accretion, causing a gap at R ∼ 5 au. Qualitatively similar features are exemplified at snow lines at 13 and 28 au. Overall, the mass transport is much more chaotic than in Model S; wind and accretion structures below the Alfvén surface are not laminar. Similar to the findings in Hu et al. (2019), accretion and decretion (flowing outwards but the streamlines not reaching super-Alfvénic) regions are associated by poloidal field loops at the midplane, as shown in the upper-middle panel of Figure 3. From 12–15 au, the averaged net accretion flow in the disk is even negative (i.e., net decretion); a similar decretion was reported in Hu et al. (2019) at approximately the same location. Such differentiation is exemplified at the R ≃ 28 au snow line. In Figure 5, it can be seen that the driving force of accretion is different from that of Model S. Around the local density minima associated with these snow lines, the equivalent viscous α reaches its local maxima, as the magnetic and fluid fluxes are more turbulent there than in the neighbors. From the point of view of magnetic stress, ${\dot{M}}_{\mathrm{acc}}({T}_{R\phi })$ and ${\dot{M}}_{\mathrm{acc}}({T}_{z\phi })$ are at about the same scale in most parts of the disk. Beyond 6 au, radial stress ( T _{R ϕ} ) induced accretion starts to play an equally important role as wind ( T _{z ϕ} ).

4.2.3. How Not to Launch Wind

One of the most distinctive features of Model L is the absence of wind within 4 au (Figure 3). We can understand this by analyzing the vertical accretion structure. Since the mass accretion rate is ${\dot{m}}_{\mathrm{acc}}=-2\pi R{\int }_{\mathrm{bot}}^{\mathrm{up}}\rho {v}_{R}{dz}$, we can rewrite Equation (8) in derivative form so we can calculate the accretion flux (mass flux at the cylindrical R direction) at different vertical layers:

$$\begin{eqnarray}&&-\rho {v}_{R}=\displaystyle \frac{2}{{{Rv}}_{K}}\displaystyle \frac{\partial }{\partial R}\left({R}^{2}{{\boldsymbol{T}}}_{R\phi }\right)+\displaystyle \frac{2R}{{v}_{K}}\displaystyle \frac{\partial {{\boldsymbol{T}}}_{z\phi }}{\partial z}.\end{eqnarray} \tag{ 13 }$$

We show a vertical slice of mass flux and Maxwell stress in Model L at R = 3 au and R = 7 au in Figure 6, with R = 7 au being the reference. At R = 7 au, the disk has a "regular" structure as seen in the mass flux and density panel in Figure 3. This is reflected in the vertical accretion analysis. Accretion (v_R < 0) is confined within a relatively thin sheet (a little more than one scale height) in the midplane, and above that we see decretion (v_R > 0) and wind. Though wind is usually measured as mass flux in the z direction, it also has an R component, thus we can still perform wind analysis at the cylindrical R direction. Both T _{R ϕ} and T _{z ϕ} drive accretion at the midplane. Closer to the disk surface there is competition between T _{R ϕ} -driven decretion and T _{z ϕ} -driven accretion and T _{R ϕ} prevails. So there is a layer of decretion at the disk surface. Continuing up to the atmosphere, at the wind region, the R component of the wind is related to T _{z ϕ} (v_R > 0) and in most radii, it dominates over T _{R ϕ} (v_R < 0), which tends to drag wind back to the central star. Adding the effect of T _{z ϕ} and T _{R ϕ} together, the combined result is the predicted flux shown in the leftmost panels of Figure 6. Mass flux at both the wind region and the disk region is well predicted by a stress calculation, as the dashed lines (predicted) and solid lines (measured ρ v_R ) overlap a lot in Figure 6.

Zoom In Zoom Out Reset image size

At this point we can understand how not to launch wind at R = 3 au. Contribution from T _{R ϕ} and T _{z ϕ} cancels out in the atmosphere when z/R > 0.2. Both measured and predicted mass flux are almost null as shown in the upper leftmost panel in Figure 6. Magnetic field lines threading the disk surface need to be bent both vertically (opening angle outward, B_R > 0) and azimuthally (strong ϕ component) to lift material up. In the upper-middle panel of Figure 3, the red arrow at 4 au marks the field lines that are almost purely vertical, with near zero B_R and B_ϕ , which is an order of magnitude smaller. Looking at the field lines at the same location in Model S (lower-middle panel), we find that they all lean outwards with a significant positive R component. This is why disk wind launches normally in Model S but not in the innermost regions of Model L.

This work discusses the impact of grain sizes on nonideal MHD accretion processes and subsequently the formation of substructures in PPDs. Dust grains are categorized into two populations by size, which have qualitatively different levels of impact on the substructure formation in PPDs. PAH-sized grains behave as ions when charged, and thus do not correspond to significant features in nonideal MHD diffusivity. Larger (≳ submicron) charged grains collide with neutrals with a much bigger cross section, whose number density is related to diffusivities directly. Physical processes that affect the radial distribution of larger grains (e.g., by snow lines) lead to the formation of rings and gaps via the radial variation of accretion rates (dominating) and wind launching rates (secondary importance).

Structure formation due to the increase in the nonideal MHD accretion rates lost efficiency wherever PAHs were much more abundant (in terms of effective cross section) than submicron-sized grains. Therefore, the ubiquity of PPD substructures suggests that PAHs could be relatively rare in the absence of planet-induced structure formation. This inference accords with the lack of PAH signals detected in low mass embedded young stellar objects (Geers et al. 2009). The sublimation temperature of hydrocarbons is roughly around 400 K (e.g., Anderson et al. 2017), so when it solidifies in the disk (most of it actually), there is a high chance that a PAH grain may collide and coagulate with a larger dust grain. This way the effect of PAH in substructure formation can be ignored.

Our nonideal MHD simulations have connected ionization chemistry (especially dust-related processes) to substructure formation through accretion dynamics, but there are still a few caveats. When we expand the chemical calculation from PAH only to submicron-sized grains, the assumption that all grains are only singly charged starts to show its limits. To better account for multiple charged grains, the electric potential of an already charged grain needs to be considered. The variation in grain sizes feeds back to the rings' locations via snow lines, as the severity of traffic jams in accretion, in turn, depends on dust size and location. Therefore, the majority of dust grains redistributed by snow lines are not necessarily the group that dominates substructure formation via magnetic diffusion. In our current models, dust grains are still single-sized species co-moving with gas, which is still insufficient to account for the dust redistribution for the complete loop of feedback in the substructure formation. Radial drift alone is not efficient even for the larger grains in Model L; those grains that concentrate rapidly at the pressure maxima cannot dominate the balance of ionization as they are much less in number density. Such concentration of dust, in the meantime, can also result in more fragmentation and hence transfer mass in dust rings to smaller grains. Thus, a complete multispecies dust model needs to be implemented with multiple dust sizes. A viable method is to use the two-population (big and small) recipe of dust that involves coagulation and fragmentation processes for the transfer of mass between these two populations (e.g., Birnstiel et al. 2012; Kanagawa et al. 2018; Tamfal et al. 2018). In addition, particle-based dust grain models are critical to study the vertical distribution of the two-population dusts for their sedimentation, settling, and PAH freezing processes, which are likely to concentrate the grains near the equatorial plane and make the interior of the disk more susceptible to the sizes of grains.

X.H. and Z.Z. acknowledge support from the National Aeronautics and Space Administration through the Astrophysics Theory Program under grant No. NNX17AK40G. Simulations were carried out with the support of the Texas Advanced Computing Center (TACC) at The University of Texas at Austin through XSEDE grant TG-AST130002. X.H. acknowledges support from the Virginia Initiative on Cosmic Origins (VICO) via a Cosmic Origins Postdoctoral Fellowship. L.W. acknowledges support from Center for Computational Astrophysics of the Flatiron Institute. S.O. is supported by JSPS KAKENHI grant Nos. JP16K17661, JP18H05438, and JP19K03926.