A Tale of Two Grains: impact of grain size on ring formation via nonideal MHD processes

Studying the dynamics of protoplanetary disks is crucial for not only constructing comprehensive planet formation theory, but also understanding the fine features in recent high resolution observations. The 2014 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, are found to be common in protoplanetary disks in subsequent observations (e.g. Huang et al., 2018; Long et al., 2018).

Recent MHD simulations suggest that substructure formation can be achieved by non-ideal 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. 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 for 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 non-ideal 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 pre-calculated with dust distribution from global dust evolution including sintering effects (Okuzumi et al., 2016). The ionization structure and the non-ideal MHD diffusivities have sharp changes at the snow lines due to the change of 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 $\mathrm{CO_{2}}$ and $\mathrm{C_{2}H_{6}}$ snow lines. This dust-to-gas feedback requires much less dust to take effect, compared to hydrodynamic drag mechanism (e.g., Gonzalez et al., 2017).

However, the simple treatment of dust on magnetic diffusion is also one of the major caveats of current non-ideal MHD simulations. Grain-modified magnetic diffusivities in protoplanetary disks were first explored in the very small grain limit ($\sim$ nm size, or PAH scale). They are known for substantially reducing the level of ionization, killing magnetorotational instability (MRI)(Bai, 2011a; Wang et al., 2018). Charged grains, when very small, can be regarded as heavy ions in diffusion calculations, while its geometric size starts to take effect in larger ( sub-micron to micron-sized) grains. Xu & Bai (2016) systematically explored the dependence of magnetic diffusivities on magnetic field strength with various grain sizes, and found the dependence of magnetic diffusion onto field strength shifts over grain sizes.

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 $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., Hu et al., 2019; Béthune et al., 2017; Bai & Stone, 2017).

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

Similar to Hu et al. (2019), we use the radial profiles of dust surface density and particle size derived from a one-dimensional global dust evolution model having multiple snow lines by Okuzumi & Tazaki (2019). In this model, the grain size distribution follows a power-law with minimum and maximum grain sizes $a_{\mathrm{gr},\mathrm{min}}$ and $a_{\mathrm{gr},\mathrm{max}}$, respectively. The power law is such that the vertically integrated number density of grains per unit grain radius $a_{\mathrm{gr}}$ is proportional to $a_{\mathrm{gr}}^{-3.5}$. We fixed the minimum grain size $a_{\mathrm{gr},\mathrm{min}}$ as $0.1~{}{\rm\,\mu m}$, whereas the maximum size $a_{\mathrm{gr},\mathrm{max}}$ and total dust surface density $\Sigma_{d}$ are evolved by computing radial drift and collisional growth/fragmentation. The calculation takes into account the low stickiness of CO₂ ice (Okuzumi & Tazaki, 2019) and aggregate sintering (Okuzumi et al., 2016). In this 1-D 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 velocity dispersion of 1.7% of the sound speed. The gas surface density distribution for the dust evolution calculation is fixed in time and assumed to scale as $R^{-1}$, where $R$ is the distance from the central star, at $R<150~{}\textsc{au}$. We evolve the dust disk 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).

The profile of $\Sigma_{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 snowlines CO₂, C₂H₆, and CH₄ (located at $R=13,29,{\rm\ and\ }82~{}\textsc{au}$, respectively) experience enhanced collisional fragmentation and pile up there. Collisional fragmentation is also enhanced interior to the $\mathrm{H_{2}O}$ snow line (located at $4~{}\textsc{au}$) due to lack of water ice and results in another traffic jam of dust in this region.

The adopted dust’s total geometric cross-section for chemistry calculation is plotted in the lower panel of Figure 1. Multiplying average cross-section of grains of all sizes and total dust abundance (number density over Hydrogen nucleus) we get dust’s total geometric cross-section per Hydrogen nucleus. This number quantifies the contribution of dust particles to the whole chemical network. We make it unchanged in different models that are to be presented below. Note that we enhanced $\sigma_{\mathrm{gr}}/\mathrm{H}$  when $R<4~{}\textsc{au}$, to ensure a moderate ionization fraction at inner disk.

In a weakly ionized protoplanetary disk, the equation of motion for charged species (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,

where for the charged species indicated by $j$, $Z_{j}$ is its charge number (in unit charge $e$), $m_{j}$ is its mass, $\bm{v}_{j}$ is its drift velocity relative to the neutral background, $\gamma_{j}\equiv\langle\sigma v\rangle/(m+m_{j})$ ($m$ is the averaged particle mass of the neutrals) and $\langle\sigma v_{j}\rangle$ is its rate coefficient of momentum transfer with neutrals. The electric field in the frame comoving with the neutrals is denoted by $\bm{E}^{{}^{\prime}}$.

In the case of very tiny grains, the interaction between charged grains and neutral molecules are mediated by electric field from induced electric dipoles in the neutrals. This $r^{-4}$ electric potential makes the collision rate coefficient $\langle\sigma v\rangle_{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:

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

So the transition from electric potential dominated cross section to geometric cross section is at $\sim nm$ scale. Given $T$=100K, any single charged grain with $a_{\mathrm{gr}}>5.7\times 10^{-8}~{}{\rm\,cm}$ needs to consider geometric effect when calculating collision coefficient. This difference is shown in the Ohmic, Hall and Pederson conductivities:

Here the summation index $j$ runs through all charged species, with $n_{j}$ indicating the number density and charge of individual charged species. Since $\beta_{gr},\beta_{i}\ll 1$ stands true throughout most regions of our disk, the difference between ions and grains is negligible for Hall conductivity $\sigma_{\mathrm{H}}$. For Ohmic and Pederson 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 the three magnetic diffusivities are (Bai, 2011b; Wang et al., 2018).

In general, when we deal with protoplanetary disks, the diffusivities $\eta_{O}$ and $\eta_{A}$ increase with larger grains due to aforementioned analyses.

This letter studies PPDs that are subject to the impact of spatial distribution of dust grains. For brevity, we refer the readers to Wang et al. (2018) (WBG19 hereafter) for the details in (1) the methods of global full MHD simulations combined with consistent thermochemistry and ray-tracing radiation, (2) the initial conditions. 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_{\phi}=-B_{r}|_{t=0}$ (the initial value of $r$ component) to suppress magnetic instabilities there. (note that the magnetization here is $\sim 20$ times WBG19 in terms of absolute magnetic pressures and stresses. The plasma $\beta$ at midplane is $2\times 10^{4}$, 5 times smaller then WGB19) Other hydrodynamic and field components are identical to WBG19.

The dust grains are still treated as single-sized carbonaceous grains co-moving with the gas. The size is either $a_{\mathrm{gr}}=5~{}\mathrm{\AA}$ (Model S) or $a_{\mathrm{gr}}=10^{-2}~{}{\rm\,\mu m}$ (Model L). On the other hand, the variation of their properties is reflected by adjusting their total abundance so that the total geometric cross-section of dust grains per hydrogen nucleus, $\sigma_{\mathrm{gr}}/\mathrm{H}$, matches the value obtained from the dust evolution calculation (bottom panel of Figure 1), as a function of the cylindrical radius $R$ to the central star. The basic properties of our fiducial model (Model S) are summarized in Table 1. Maximum $\sigma_{\mathrm{gr}}/\mathrm{H}$ is $8\times 10^{-23}~{}{\rm\,cm}^{2}$ in both models, which equals to $7\times 10^{-6}$ in dust to gas mass ratio in Model S, and $1.4\times 10^{-4}$ in Model L, if we take the average density of a single grain as $2.25~{}{\rm\,g}~{}{\rm\,cm}^{3}$. 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 $\langle\sigma v\rangle_{\mathrm{gr}}$ that is two orders of magnitude larger than the grains in Model S. The behavior of very small grains are actually similar to ions, as the Hall parameter of ions is approximately $\beta_{i}\approx 3.3\times 10^{-3}(B_{\rm G}/n_{15})$. For larger charged grains, smaller $\beta_{gr}$ indicates weaker coupling to magnetic fields, which means that grains are more prone to frequent collision with neutrals.

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. In Model S, where grains have the same size as WBG19 ($5~{}\mathrm{\AA}$ PAH) , the variation of dust number density yields some effects on surface density beyond $10~{}\textsc{au}$, and there are no apparent features associated with $\sigma_{\mathrm{gr}}/\mathrm{H}$ bumps at $13~{}\textsc{au}$ and $28~{}\textsc{au}$. In Model L, the location of rings and gaps at these two snow lines are very close to Hu et al. (2019). From the non-ideal MHD point of view, the simulations prove that number density variation of PAH-scaled grains has sub-linear impacts on magnetic diffusivities. The two middle panels shows the Elsasser numbers of both Ohmic resistivity and Ambipolar diffusion of Model S, which have mild change at $13~{}\textsc{au}$ and $28~{}\textsc{au}$ snow lines, comparing to order of magnitude jump of $\sigma_{\mathrm{gr}}/\mathrm{H}$. The diffusivity profiles of Model L reflect the $\sigma_{\mathrm{gr}}/\mathrm{H}$ structures quite well. Diffusion strength increases up to 4 orders of magnitude at those two snow lines, where the total $\sigma_{\mathrm{gr}}/\mathrm{H}$ rises by about 10 times.

After 1500 orbits at the innermost radius, the overall disk structure of both simulations reaches a quasi-steady phase, as shown in Figure 3. 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 local sound speed. Above the wind base, the Alfvén surface is where the velocity in the poloidal plane equals to the poloidal Alfvén velocity $v_{A,p}=\sqrt{(v_{z}^{2}+v_{R}^{2})/4\pi\rho}$.

This work discusses the impact of grain sizes on the non-ideal MHD accretion processes and subsequently formation of sub-structures in PPDs. Dust grains are categorized into two populations by size, which have qualitatively different levels of impact onto the sub-structure formation in PPDs. PAH-sized grains behave as ions when charged, thus do not correspond to significant non-ideal MHD diffusivity features. Larger ($\gtrsim$sub-micron) charged grains interact with neutrals with 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 jump of the non-ideal MHD accretion rates would lose efficiency wherever PAHs are much more abundant (in terms of effective cross-section) than sub-micron-sized grains. Therefore, the ubiquity of PPD sub-structures suggests that PAHs could be relatively rare in absence of planet-induced structure formation. This inference accords with the lack of PAH signals detected in low mass embedded YSOs (Geers et al., 2009). In the dense regions where the temperature is relatively low, PAHs tend to freeze out and be retained by the surfaces of larger grains (e.g., Anderson et al., 2017).

Our non-ideal MHD simulations have connected ionization chemistry (especially dust-related processes) to sub-structure formation through accretion dynamics, but they still have a few caveats. 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 sub-structure 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 feedbacks in the sub-structure formation. Radial drift alone is not efficient even for the larger grains in Model L; those grains that concentrate rapidly at 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 hence transfer mass in dust rings to smaller grains. Thus, a complete multi-species 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; Tamfal et al., 2018; Kanagawa 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 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 with Grant No. NNX17AK40G. Simulations are carried out with the support from the Texas Advanced Computing Center (TACC) at The University of Texas at Austin through XSEDE grant TG-AST130002. L.W. acknowledges support from Center for Computational Astrophysics of Flatiron Institute. S.O. is supported by JSPS KAKENHI Grant Numbers JP16K17661, JP18H05438, and JP19K03926