# Turbulent Dust-trapping Rings as Efficient Sites for Planetesimal Formation

Ziyan Xu, Xue-Ning Bai

<figure><img src="image_1.png" /></figure>

## To cite this version:

Ziyan Xu, Xue-Ning Bai. Turbulent Dust-trapping Rings as Efficient Sites for Planetesimal Formation. The Astrophysical journal letters, 2022, 937, 10.3847/2041-8213/ac8dff . insu-03854058

HAL Id: insu-03854058

https://insu.hal.science/insu-03854058

Submitted on 15 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
# Turbulent Dust-trapping Rings as Efficient Sites for Planetesimal Formation

Ziyan Xu<fnref n="1" />,<fnref n="2" />,<fnref n="3" /> ID and Xue-Ning Bai<fnref n="4" />,<fnref n="5" /> ID

Received 2022 July 8; revised 2022 August 23; accepted 2022 August 29; published 2022 September 15

## Abstract

Recent observations of protoplanetary disks (PPDs) at submillimeter wavelengths have revealed the ubiquity of annular substructures that are indicative of pebble-sized dust particles trapped in turbulent ringlike gas pressure bumps. This major paradigm shift also challenges the leading theory of planetesimal formation from such pebbles by means of the streaming instability, which operates in a pressure gradient and can be suppressed by turbulence. Here, we conduct 3D local shearing box nonideal magnetohydrodynamic simulations of dust trapping in enforced gas pressure bumps, including dust backreaction. Under a moderate level of turbulence generated by the magnetorotational instability with ambipolar diffusion, which is suitable for outer disk conditions, we achieve quasi-steady states of dust trapping balanced by turbulent diffusion. We find strong dust clumping in all simulations near the gas pressure maxima, reaching a maximum density well above the threshold for triggering gravitational collapse to form planetesimals. A strong pressure bump concentrates dust particles toward the bump's center. With a weak pressure bump, dust can also concentrate in secondary filaments off the bump's center, due to dust backreaction, but strong clumping still occurs mainly in the primary ring around the bump's center. Our results reveal dust-trapping rings to be robust locations for planetesimal formation in outer PPDs, while they may possess diverse observational properties.

*Unified Astronomy Thesaurus concepts: Protoplanetary disks (1300); Planetesimals (1259); Planet formation (1241); Magnetohydrodynamics (1964)*

## 1. Introduction

Planetesimals are ~1–1000 km sized bodies in protoplanetary disks (PPDs) and the building blocks of planets. Following the growth of dust grains to millimeter–centimeter sizes in PPDs, planetesimals are believed to form from such grown dust particles, representing the intermediate stage of planet formation. This stage is also considered to be the most difficult, as surface forces can no longer stick particles together upon collisions, nor can self-gravity, which dominates among solids of larger sizes (Chiang & Youdin 2010). Generally, planetesimal formation is understood as a two-phase process, as particles first concentrate into high-density clumps, followed by their gravitational collapse. What drives the formation of particle clumps thus becomes the key to understanding planetesimal formation. Regardless of the specific mechanism, particles of millimeter–centimeter sizes strongly interact with gas via aerodynamic drag, and the gas dynamics in PPDs, especially the disk turbulence, must play an important role.

The leading theory of planetesimal formation is the streaming instability (SI; Youdin & Goodman 2005) resulting from the two-way drag force between gas and dust in the presence of a background pressure gradient in PPDs. Hybrid gas–dust simulations have found that the nonlinear evolution of the SI leads to dust clumping (e.g., Johansen et al. 2009; Bai & Stone 2010a; Carrera et al. 2015; Yang et al. 2017), provided

that the dust abundance Z (the ratio of dust to gas surface density) exceeds some threshold (typically ≳0.02; but see Li & Youdin 2021). While highly successful, the vast majority of SI simulations have assumed no external turbulence. However, recent studies have found that even a modest level of external disk turbulence suppresses the SI (Chen & Lin 2020; Gole et al. 2020; Umurhan et al. 2020), challenging the SI paradigm.

The outer disk regions (≫20 au) are the most accessible to spatially resolved disk observations. At submillimeter wavelengths (where millimeter-sized dust radiates most efficiently), the prevalence of ringlike substructures has recently been firmly established (e.g., Andrews et al. 2018a; Long et al. 2018), which are widely believed to be indicative of dust trapping in pressure bumps (Dullemond et al. 2018).<fnref n="6" /> Dust trapping in local pressure maxima was first demonstrated by Whipple (1972), and the formation mechanisms of pressure bumps can be planet–disk interactions (see the review in Paardekooper et al. 2022) or have nonplanetary origins, such as magnetohydrodynamic (MHD) zonal flows (e.g., Suriano et al. 2018; Riols et al. 2020; Cui & Bai 2021). While there is direct and indirect evidence pointing to weak disk turbulence (Pinte et al. 2016; Flaherty et al. 2017, 2018; Teague et al. 2018), modeling dust ring widths suggests a modest level of turbulence (Dullemond et al. 2018; Rosotti et al. 2020). Together, these observational constraints lead to an emerging picture that most millimeter-sized dust particles in outer PPDs are most likely concentrated in weakly to modestly turbulent

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ Univ Lyon, Univ Lyon 1, ENS de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, F-69230, Saint-Genis-Laval, France
ziyan.xu@ens-lyon.fr</footnote>

<footnote n="2">² Kavli Institute for Astronomy and Astrophysics, Peking University, Yiheyuan 5, Haidian Qu, 100871 Beijing, People's Republic of China</footnote>

<footnote n="3">³ Department of Astronomy, Peking University, Yiheyuan 5, Haidian Qu, 100871 Beijing, People's Republic of China</footnote>

<footnote n="4">⁴ Institute for Advanced Study, Tsinghua University, 100084 Beijing, People's Republic of China; xbai@tsinghua.edu.cn</footnote>

<footnote n="5">⁵ Department of Astronomy, Tsinghua University, 100084 Beijing, People's Republic of China</footnote>

<footnote n="6">⁶ Unless otherwise noted, “pressure bumps” in this paper refer to dust-trapping pressure bumps, i.e., they are strong enough to overcome the background pressure gradient and to generate local pressure maxima.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** A snapshot of a 3D visualization of the dust spatial distribution in our simulation, with a strong gas bump (run Z2H) at the time of $T = 1380\Omega_K^{-1}$ after inserting dust particles. The dust density is mapped to the brightness, where the brighter regions indicate higher dust densities. Also shown in the figure are the vertically averaged gas density, projected to the $x-y (R-\phi)$ plane, and the vertically and azimuthally averaged radial profiles of the dust and gas densities, with the amplitudes of the profiles being rescaled for better illustration.</figcaption></figure>

and largely axisymmetric gas pressure bumps. The strongly enhanced local dust abundance in such pressure bumps yields favorable conditions for planetesimal formation. However, given that the SI is unexpected to operate without a pressure gradient, and can be suppressed by external turbulence, it is unclear whether it will further lead to dust clumping and hence planetesimal formation.

To address the feasibility of planetesimal formation in the more realistic environment of turbulent dust rings, we conduct hybrid gas-particle simulations in the standard shearing box framework that mimics a local patch of the disk. It is anticipated that the magnetorotational instability (MRI) will be a primary mechanism for driving turbulence in the outer PPDs. Due to a very low level of disk ionization, the MRI is weakened by ambipolar diffusion (AD), as the dominant nonideal MHD effect in the outer disk. The turbulence in our simulations is self-consistently generated by the MRI with AD. We further implement a novel forcing prescription to drive a pressure bump in the center of our simulation box, with essentially no impact on the disk turbulence. Our simulations generalize from our earlier work (Xu & Bai 2022; hereafter, Paper I), by imposing and enforcing the pressure bumps, allowing us to achieve a steady-state configuration with dust trapping toward pressure bumps, balanced by dust diffusion in realistic disk turbulence, in a self-consistent and well-controlled manner. They differ from the recent works of Carrera et al. (2021, 2022) and Carrera & Simon (2022), who conducted simulations using a similar procedure to enforce pressure bumps, but without external turbulence, meaning that they were not able to achieve such a steady-state balance of dust trapping and diffusion as has been envisioned in recent disk observations (Dullemond et al. 2018; Rosotti et al. 2020).

This letter is organized as follows. We first introduce the implementation of the pressure bumps and the simulation setup

in Section 2, before presenting our main results in Section 3. In Section 4, we further analyze the results, focusing on dust-clumping conditions and ring dynamics, and this is followed by a comparison with recent studies and a discussion of observational implications. We summarize and discuss future prospects in Section 5.

## 2. Method and Simulations

We use the Athena MHD code (Stone et al. 2008) to perform 3D hybrid gas-particle local shearing box simulations, considering nonideal MHD with AD. The formulation and basic simulation setup is similar to that of Paper I, which we only describe briefly here. The $x, y, z$ coordinates in a local shearing box correspond to the radial, azimuthal, and vertical directions in a global disk, as illustrated in Figure 1. We use an isothermal equation of state, with sound speed $c_s = 1$ and local Keplerian frequency $\Omega_K = 1$ in code units, which define our length unit of disk scale height $H = c_s/\Omega_K = 1$. The gas density, velocity, and pressure are denoted by $\rho, v$, and $P = \rho c_s^2$, respectively. The dust particles are characterized by their dimensionless stopping time $\tau_s \equiv t_s \Omega_0$ and the height-integrated dust-to-gas mass ratio $Z$ for the background gas surface density. Our simulations focus on the midplane region, and hence we neglect vertical gravity in the gas (unstratified), but vertical gravity is included for the particles, to allow dust settling. We impose a net vertical magnetic flux, characterized by $\beta_0 = 12,800$, the ratio of the background gas pressure to the magnetic pressure of the net vertical field. The AD is characterized by the Ambipolar Elsässer number $Am$, where we fix $Am = 2$, which is appropriate for typical outer disk conditions.
## 2.1. Implementation of the Gas Pressure Bump

On top of the simulation setup in Paper I, we further implement a gas pressure bump by using a prescribed forcing procedure. Let $\rho_0 = 1$ be the constant background gas density. The forcing will redistribute the density, and we aim to achieve a Gaussian-like density (and hence pressure) bump in our simulation box of the following form:

$$ \frac{\rho(x)}{\rho_0} = 1 + A_{\text{bump}} \exp \left[ - \frac{x^2}{2w_{\text{bump}}^2} - 0.2 \left( \frac{x^2}{2w_{\text{bump}}^2} \right)^4 \right], \quad (1) $$

where $A_{\text{bump}}$ is the dimensionless amplitude of the pressure bump and $w_{\text{bump}}$ is the effective width of the pressure bump. The first term in the exponential gives the overall Gaussian profile, while the second term damps the profile toward the domain edges, to satisfy the shearing-periodic radial boundary condition.

The radial pressure gradient associated with this bump must be balanced by an azimuthal flow:

$$ \frac{v'_{y, \text{bump}}(x)}{c_s} = \frac{H}{2} \frac{d}{dx} \ln \rho(x), \quad (2) $$

where the prime ' in $v_y$ denotes deviation from the Keplerian velocity $v_K = -(3/2)\Omega_K x$. Note that in our simulations, as we focus on regions near the pressure maxima, we do not impose any background pressure gradient.

To generate and sustain the bump, we apply an azimuthal forcing in our simulation, so that the azimuthal velocity relaxes to the desired profiles within a timescale of $t_{\text{relax}} = 60\Omega_K^{-1}$:

$$ \Delta v'_{y}(x) = [v_{y, \text{bump}}(x) - \bar{v}'_{y}(x)][1 - \exp(-\Delta t/t_{\text{relax}})], \quad (3) $$

where $\bar{v}'_y(x)$ is azimuthally and vertically averaged. Our forcing procedure is analogous to that of Carrera et al. (2021), but with several differences. First, our forcing acts on azimuthal velocity (momentum) only, thus it is mass conserving. It directly mimics the torque exerted on the disk that is responsible for making the bump (whether it be of planetary origin or nonplanetary origin, as in Cui & Bai 2021), and the desired density profile is a direct consequence from this forcing. Second, our forcing is homogeneously applied in the y and z directions, as opposed to relaxation at individual grid cells, which introduces local friction. Third, our relaxation time is much longer, as opposed to $\sim 1\Omega_K^{-1}$. These facts minimize the impact of our forcing procedure on disk turbulence, and especially avoid artificial damping at small scales.<fnref n="7" />

## 2.2. Simulation Procedures and Parameters

A list of our simulation runs is given in Table 1. We explore the strengths of the gas pressure bump, the presence/absence of dust feedback, different levels of solid abundance, and simulation resolutions. The simulation setup and parameters are described as follows.

<table><caption>**Table 1**
List of Simulation Runs</caption><thead><tr><th>Run</th><th>Z</th><th>N<sub>x</sub> × N<sub>y</sub> × N<sub>z</sub></th><th>A<sub>bump</sub></th><th>H<sub>d</sub>/H</th><th>w<sub>d</sub>/H</th></tr></thead><tbody><tr><td>Z0</td><td>1e-30</td><td>512 × 512 × 64</td><td>0.5</td><td>0.067</td><td>0.55</td></tr><tr><td>Z2</td><td>0.02</td><td>512 × 512 × 64</td><td>0.5</td><td>0.014</td><td>0.16</td></tr><tr><td>Z0H</td><td>1e-30</td><td>1024 × 1024 × 128</td><td>0.5</td><td>...</td><td>...</td></tr><tr><td>Z1H</td><td>0.01</td><td>1024 × 1024 × 128</td><td>0.5</td><td>0.025</td><td>0.25</td></tr><tr><td>Z2H</td><td>0.02</td><td>1024 × 1024 × 128</td><td>0.5</td><td>0.017</td><td>0.16</td></tr><tr><td>Z0w</td><td>1e-30</td><td>512 × 512 × 64</td><td>0.25</td><td>0.075</td><td>0.83</td></tr><tr><td>Z1w</td><td>0.01</td><td>512 × 512 × 64</td><td>0.25</td><td>0.031</td><td>0.87</td></tr><tr><td>Z2w</td><td>0.02</td><td>512 × 512 × 64</td><td>0.25</td><td>0.022</td><td>0.74</td></tr></tbody></table>

**Note.** All of the simulations have a box size of $8H \times 16H \times H$, with $\beta_0 = 12,800$, $Am = 2$. The gas bump width is $w_{\text{bump}} = 2H$, with no background pressure gradient. The particle numbers are two per cell, on average, with the particle stopping time being $\tau_s = 0.1$. All of the parameters in this table are for the initial setup, except $H_d$ and $w_d$, which are measured from simulation.

We first run simulations without dust particles for a time of $t_0 = 120\Omega_K^{-1}$, for the MRI to grow and fully saturate into turbulence, following our previous work (Paper I). The gas pressure bump is initiated at the beginning of the simulation, with a bump width of $w_{\text{bump}} = 2H$. We explore two cases, strong and weak bumps, with bump amplitudes of $A_{\text{bump}} = 0.5$ and 0.25, separately, where the weak-bump runs are marked by the presence of “w” in the run name.<fnref n="8" /> Our choice of bump width generally agrees with observational constraints (Dullemond et al. 2018; Rosotti et al. 2020), and both the bump width and the amplitudes in our simulation are reasonably consistent with the results of Cui & Bai (2021).

Dust particles are then inserted into the simulation box at the time of $t_0 = 120\Omega_K^{-1}$, without backreaction to the gas, with random horizontal positions and a vertical distribution following a Gaussian profile $\propto \exp(-z^2/2H_{d0}^2)$, with $H_{d0} = 0.02H$. We define $T = T_{\text{sim}} - t_0$, where $T_{\text{sim}}$ is the actual simulation time, and we use $T$ instead of $T_{\text{sim}}$ throughout this paper. We use a single particle species with a fixed stopping time $\tau_s = 0.1$, corresponding to particles of millimeter–centimeter sizes in the midplane regions of the outer disk in standard disk models (e.g., Chiang & Youdin 2010). As the dust particles passively follow the gas, the settling balanced by turbulent diffusion in the vertical direction causes them to maintain a certain scale height (see Paper I), while the balance between the trapping toward pressure maxima and turbulent diffusion in the radial direction allows them to concentrate into a ring with a certain ring width $w_d$ (see Figure 2), approaching a steady state after $T = 600\Omega_K^{-1}$. Next, we gradually introduce dust feedback to the gas, by gradually increasing the dust abundance to the desired value, over a time interval of $30\Omega_K^{-1}$. We consider two global dust abundances, $Z=0.01$ and $Z=0.02$, and these runs are labeled by “Z,” followed by a number indicating the $Z$ values. In all simulations, there are on average two particles per cell, ensuring sufficient particle statistics for studying dust dynamics.

We choose a simulation box size of $8H \times 16H \times H$ in the x, y, and z directions, respectively. We choose a fiducial

<footnote n="7">⁷ Our choice of relaxation time is generally consistent with the typical timescale for magnetic flux evolution (and thus bump evolution) indicated in Cui & Bai (2021). We further tested $t_{\text{relax}} = 30\Omega_K^{-1}$ and $120\Omega_K^{-1}$, respectively (not shown in the paper), and the results are not sensitive to the change of $t_{\text{relax}}$ in our simulation in this paper.</footnote>

<footnote n="8">⁸ We have confirmed that these bump amplitudes are below the threshold for triggering the Rossby wave instability (RWI) for all azimuthal wavenumbers under global settings (Ono et al. 2016), and we anticipate that the presence of the MRI turbulence will also make our pressure bump less prone to the RWI.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Time evolution of the dust maximum density (upper panel), the dust ring width (middle panel), and the dust scale height (lower panel) for all our simulation runs. The dust maximum densities for the no-feedback runs (Z0, Z0H, and Z0w) are scaled to $Z = 0.02$.</figcaption></figure>

resolution of 64/H (in the x and z directions; the resolution in the y direction is half this). We further conducted simulations with a higher resolution of 128/H for comparison, labeled with “H” at the end of the run names.

### 3. Results

An overview of our simulations is shown in Figure 1, featuring a 3D visualization of the dust spatial distribution for our run Z2H at time $T = 1380\Omega_K^{-1}$. Our forcing procedure generates a wide and steady gas bump, with additional density fluctuations (i.e., density waves) as a result of the MRI turbulence (Heinemann & Papaloizou 2009). Clearly, the dust concentrates at the center of the gas bump, its width being much narrower than the gas bump width. Besides the settling and trapping, the dust distribution is highly nonsmooth and clumpy under background turbulence.

In this section, we present our main simulation results. We measure the maximum dust density $\rho_{d,max}$, normalized by the background gas density $\rho_0$, as a proxy for dust clumping. The dust densities in the no-feedback runs (Z0 and Z0H) are rescaled to $Z = 0.02$ for better comparisons. For reference, we consider $\rho_{d,max} \gtrsim 200\rho_0$ as a threshold condition for triggering gravitational collapse, and hence planetesimal formation, which is applicable for typical outer disk conditions (see Equation (16) and Section 5.1 in Paper I). We further measure the dust ring width $w_d$ and dust scale height $H_d$, defined as the rms value of the radial and vertical coordinates for all dust particles. The

time history of $\rho_{d,max}/\rho_0$, $w_d$, and $H_d$ from all simulations is shown in Figure 2. The time-averaged values of $w_d$ and $H_d$ are presented in Table 1, averaged within $T = 800–1500\Omega_K^{-1}$. These are accompanied by Figure 3, where we show the projected dust densities in the x-y and x-z planes for the Z0, Z2H, Z02, and Z2w runs in the saturated states. We find that the results for the strong-bump and weak-bump cases are distinct, which we discuss separately below.

#### 3.1. The Strong-bump Case: Dust Clumping in a Single Thin Ring

The top panel of Figure 2 shows that the dust density presents a modest concentration in the absence of dust feedback, with $\rho_{d,max}/\rho_0$ reaching ~100, approaching yet barely reaching the threshold for triggering planetesimal formation. This is due to a combination of dust trapping and intrinsic fluctuations. Reading from the second and third panels of Figure 2, we find a dust ring width of $w_d \sim 0.55H$ and a scale height of $H_d \sim 0.067H$ (comparable to ~0.08H, as found in Paper I, without an enforced pressure bump). This yields a midplane density enhancement factor of ~50 at the bump's center. Indeed, as shown in Figure 3, the mean midplane density at the bump's center for run Z0 (renormalized assuming $Z=0.02$) is about $\bar{\rho}_d \sim 1$. We do not observe obvious dust clumping from the averaged dust densities along the vertical and azimuthal directions, and thus attribute an additional enhancement factor of ~50–100 in $\rho_{d,max}$ from
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Snapshots of the projected dust densities in the x-z (upper four panels) and x-y (lower four panels) planes for the runs Z0, Z2H, Z0w, and Z2w at time $T = 1380\Omega_K^{-1}$, when the secondary filament in Z2w is clearly seen. The dust density in the x-z plane is azimuthally averaged over the whole domain, while in the x-y plane, the dust density is vertically averaged around the midplane ($z = \pm 0.05 H$), to highlight the dust structures.</figcaption></figure>

localized density fluctuations. The results are largely indepen-
dent of resolution, as ρ<sub>d,max</sub>, w<sub>d</sub>, and H<sub>d</sub> are consistent with each
other between run Z0 and Z0H, within fluctuation levels.

Turning on dust feedback leads to a rapid boost in ρ<sub>d,max</sub>/ρ<sub>0</sub>,
well exceeding 10<sup>3</sup> (and hence the Roche density) for all runs
with Z = 0.01 (Z1H) and Z = 0.02 (Z2 and Z2H), which is
therefore expected to trigger planetesimal formation. Visible
bright spots are also seen in the bottom panel of Figure 3 for
run Z2H, suggesting that such clumping is physical instead of
random Poisson fluctuations. There is a trend of higher
resolution leading to a higher level of dust clumping, where,
for runs Z1H and Z2H, ρ<sub>d,max</sub>/ρ<sub>0</sub> ≳ 5000, an enhancement
factor of ~50 compared to the no-feedback case. We have

further examined the cumulative distribution function (CDF;
see Appendix A), which reveals that the dust density is indeed
substantially boosted, with a significant fraction (a few tens of
percent) residing in regions that are over the Roche density.

Somewhat surprisingly, with dust feedback, we observe a
dramatic reduction of the ring width to w<sub>d</sub> ≈ 0.2H, ~3 times
narrower than that without feedback. Moreover, the thickness
of the dust layer is also reduced substantially, by a factor of
3-4. Examining Figure 3 reveals that this reduction is not
homogeneous: there is a strong/spiky dust concentration and
more substantial settling at the bump's center (see also the top
left panels of Figures 3 and 4 in Paper I), which reduces the
measured w<sub>d</sub> and H<sub>d</sub> toward lower values. We note that the
reduction of $H_d$ was also reported in Paper I, where it was mainly attributed to the reduction of the turbulence correlation time. However, in the presence of the pressure bump, $H_d$ is further reduced by a factor of nearly ~2. The dramatic reduction of the ring width, as well as the further promotion of dust settling, are likely the result of strong dust mass loading in the bump center (see the further discussions in Section 4.2 and Appendix C). The results are convergent with resolution, as the values of $w_d$ and $H_d$ are consistent with each other between runs Z2 and Z2H.

## 3.2. The Weak-bump Case: Diverse Ring Properties

The dust ring properties are the most prominent feature of the weak-bump case compared to the strong-bump case. The rest of the results, especially the dust clumping and settling, are similar to the strong-bump case.

The middle panel of Figure 2 shows that without dust feedback, $w_d$ in the weak-bump case is generally higher compared to the strong-bump case, by a factor of up to ~2. This is easily accounted for by the reduced gas pressure gradient on the two sides of the gas bump, resulting in weaker dust trapping. This is also seen in Figure 3, which shows a smooth and radially extended distribution of dust from run Z0w, together with azimuthally extended stripes akin to the MRI density waves.

When dust feedback is turned on, interestingly, the dust ring width $w_d$, measured from the rms radial coordinates, remains similar for all cases. However, when examining the dust spatial distribution in Figure 3, we find that instead of the dust concentrating into a single ring at the bump center, the dust dynamics in the weak bump are more complicated, potentially forming multiple fine-scale filaments within the gas bump. In the case of run Z2w, we observe two dust rings, where the primary ring is close to the bump center, and a secondary ring is more than *H* away from the bump center, in addition to a few minor filaments. The formation of this more complex ring/ filamentary structure is likely caused by the interplay between dust feedback and gas dynamics within the bump. In addition, these structures generally remain during slow evolution over tens to hundreds of $\Omega_K^{-1}$. The details of the formation and evolution of (multiple) fine-scaled ring structures will be further discussed in Section 4.2 and Appendix B.

# 4. Discussion

## 4.1. Clumping Conditions

For both the strong and weak bumps, sufficient dust clumping is seen to trigger planetesimal formation. Here, we further discuss the clumping conditions, to gain insight into the clumping mechanism. In doing so, we show in the top four panels of Figure 4 the radial profiles of the vertically and azimuthally averaged dust density, the gas density, and the Keplerian-subtracted azimuthal gas velocity $v'_y$, at the very midplane of our simulation (±1 cell). The latter two are useful for quantifying the pressure gradient, and we further show the time evolution of the radial profiles of $v'_y$ in four representative runs in the bottom panels.

**Dust feedback.** Dust feedback is required for dust clumping in our simulations, as no clumping is seen in our runs without feedback (Z0 and Z0w). As seen in the upper panels of Figure 4, the midplane dust densities in the no-feedback cases generally exhibit a Gaussian-like radial profiles, conforming to

the standard expectation of concentration balancing turbulent diffusion (Dullemond et al. 2018). With the inclusion of dust feedback, the radial density profiles substantially deviate from being Gaussian-like, showing density cusps at the locations of the dust rings, where the dust density is enhanced by a factor of nearly 100 relative to the no-feedback counterparts. The cusps are associated with the locations of the dust clumping, and could be consequences of the dust clumping itself.

**Solid abundance.** Dust clumping occurs in all simulations for a mean solid abundance Z (averaged over the simulation domain) as low as 0.01. Note that due to dust trapping, the effective solid abundance (i.e., the vertically integrated dust-to-gas mass ratio) in the bump region is much higher. We may estimate the enhancement factor to be ~4–8 in the weak-bump and strong-bump cases, based on the value of $w_d$, corresponding to $Z \gtrsim 0.04$ around the bump center in all our simulations. Therefore, the overall solid abundance Z in the bulk disk appears to be less relevant, as dust trapping can easily bring a sufficient amount of dust to the bump region.

**Local pressure maxima.** In the strong-bump cases (runs Z1H and Z2H), dust clumping primarily occurs at the central cusp, which is marked by the dashed line in Figure 4. We further see from the radial profiles of $\rho_g$ and $v'_y$ that this location exactly coincides with the location of the gas pressure maxima, corresponding to $v'_y = 0$. In the weak-bump case, the primary ring close to the bump center is also located at the pressure maxima, but in this case it has widened into a plateau with a width of ~0.5H, as seen from the gas density profile, as well as an extended region with $v'_y \approx 0$. The cusp where the strong dust clumping occurs is located within the plateau, again consistent with the requirement of the local pressure maxima, or at least a zero pressure gradient region, being a necessary condition for strong dust clumping. These findings are consistent with the results reported in Paper I, where the pressure bumps were not enforced, but generated by the MRI zonal flows (although they were somewhat unrealistic in local unstratified simulations).

**Weak clumping in the secondary ring.** For the weak-bump cases of Z1w and Z2w, in addition to the primary ring near the bump center, an offcentered “secondary” ring is identified at $x \sim -1.2H$. We show in Appendix A that the dust density CDF in the primary ring is similar to that in the strong-bump case, with strong clumping. However, the dust CDF of the secondary ring is significantly different. The mean dust density is more than a factor of ~10 lower, and only a very small fraction (\lesssim 3%) of particles reside in the regions that are over the Roche density (which exhibit as the individual spots that are seen in the bottom right panel of Figure 3), indicating that the clumping is much less efficient. Examining the density (and hence pressure) profiles in the top right panel of Figure 4, the location of the secondary ring is clearly not at the pressure maxima.

Overall, for both the strong- and weak-bump cases, dust clumping is tightly related to the presence of pressure maxima, as long as dust backreaction is included. While secondary rings/filaments can form for the weak-bump case, they generally reside beyond the pressure maxima, and do not show prominent/efficient clumping.

## 4.2. Interplay between Dust Feedback and the Gas Bump

**Bump/ring profile.** In addition to dust clumping, the interaction between the dust and gas also influences the gas
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Upper four panels: radial profiles of the vertically and azimuthally averaged dust density ($\rho_d$; the solid curves in the upper two panels), the gas density ($\rho_g$; the dashed curves in the upper two panels), and the Keplerian-subtracted azimuthal gas velocity ($v'_y$; the solid curves in the lower two panels), at the midplane (±1 cell from $z=0$), averaged between $T = 800–1500\Omega_K^{-1}$. The vertical dashed lines mark the locations of the dust rings (the peaks of dust density), with colors corresponding to those used for the radial profiles of each run. Lower four panels: time evolutions of the radial profiles of $v'_y$, between $T = 800–1500\Omega_K^{-1}$.</figcaption></figure>

bump properties, leading to diversity in the dust structures, as we discuss here.

Dust feedback changes the radial profile of the gas bump. This can be seen in the yellow region in the time evolution map of $v'_y$ at the midplane in Figure 4, by comparing the upper and lower panels. In the strong-bump case, as well as around the primary ring in the weak-bump case, the dust feedback leads to a plateauing of the Keplerian rotation at the bump's center, in

comparison to that in the runs Z0 and Z0w. This is consistent with the previous finding in Huang et al. (2020).

The secondary dust ring, despite not occurring at the pressure maxima, can be sustained for at least dozens of orbits, without drifting toward the bump's center. Examining the profile of $v'_y$, we find that its value shows a dip at the secondary dust filament, implying near-Keplerian rotation. This is the result of feedback from the dusty filaments modifying the
rotational profile of the background gas. It likely occurs in our weak-bump runs because the gas pressure gradient is not strong enough to force the heavily loaded dust filaments to drift to the bump's center, resulting in dust trapping.<fnref n="9" /> Such a location shares similarities with the pressure maxima, in that there is almost no radial or azimuthal drift between the gas and dust, and is likely responsible for the weak clumping discussed in Section 4.1. These secondary rings also evolve over a timescale of tens to hundreds of orbits, as we discuss further in Appendix B.

**Turbulent diffusion.** The enhanced dust mass loading in the turbulent dust rings also affects the properties of the background MRI turbulence. With feedback included, we find that the gas rms velocities in both the radial and vertical directions are largely reduced in the dust ring (see the details in Appendix C). Such reductions can sufficiently explain the results reported in Section 3.1, where the dust rings in the strong-bump case become much narrower when dust feedback is included, and the dust becomes even more settled than that reported in Paper I.

### 4.3. Comparisons with Other Recent Studies

Recently, Carrera et al. (2021, 2022) and Carrera & Simon (2022) have conducted a series of hydrodynamic simulations of the SI, imposing a pressure bump by means of a similar (albeit more dissipative) forcing procedure as ours, but without external turbulence. However, the pressure bumps in most of their simulations are not sufficiently strong to overcome the background pressure gradient. As a result, particles are accumulated (but not trapped) in regions with lower pressure gradients, making such regions the most favorable areas for clumping and planetesimal formation (e.g., Bai & Stone 2010b). In other words, their simulations have identified favorable conditions for dust clumping by the conventional SI in the presence of a radial-dependent pressure gradient. Our simulations, on the other hand, focus on planetesimal formation in dust-trapping rings, i.e., the pressure maxima. The problem is only well-defined in the presence of external turbulence, as expected from observations (Dullemond et al. 2018; Rosotti et al. 2020), otherwise dust would become indefinitely concentrated toward the ring center over the trapping timescale (this issue may also apply to the analytic work of Auffinger & Laibe 2018), forming planetesimals along the way. Overall, the physical picture of our simulations is distinctly different from the works of Carrera et al. (2021, 2022), and is likely to better reflect the outer disk conditions. Together with Paper I, we show that for the MRI turbulence, the conditions for planetesimal formation require the presence of pressure maxima, which is again different from the conventional SI scenario.

The outer disk may also be subject to turbulence driven by the vertical shear instability (VSI; Nelson et al. 2013). Using 2D global simulations, Schäfer et al. (2020) found that the SI can coexist with the VSI, and that the VSI even facilitates the development of the SI and dust clumping, likely due to the overdensities induced by the VSI. We note that the

overdensities from the VSI itself are also insufficient to overcome the background pressure gradient. This is different from the case with the MRI, where we found that dust clumping requires the presence of pressure maxima, suggesting that dust clumping conditions likely depend on the nature of the turbulence. Lehmann & Lin (2022) find that the presence of an initial pressure bump (without forcing) further facilitates dust concentration, though their resolution is insufficient to capture dust clumping. Moreover, the VSI leads to vortex formation in 3D simulations (Manger et al. 2020), and hence dust clumping in dust-laden vortices (e.g., Raettig et al. 2021), representing another alternative for planetesimal formation in turbulent disks. It is worth noting that the VSI and MRI can coexist (Cui & Bai, 2022, submitted), and the turbulence is dominated by the MRI/VSI when the ambipolar Elsässer number $Am$ is large/small compared to unity. It is likely that the clumping conditions still depend on the nature of the turbulence (whether due to the VSI or the MRI), although more investigations are needed to understand the clumping mechanisms, and how the dust dynamics are affected when both instabilities are present.

### 4.4. Observational Implications

**Dust clumping and fine-scale substructures.** There are multiple observational diagnostics for probing the turbulent dust rings. First, as such dust rings are clumpy, there can be both optically thick and optically thin regions within the ring. Atacama Large Millimeter/submillimeter Array (ALMA) observations at submillimeter wavelengths can hardly resolve such regions, but the clumpiness could be reflected in the spectral indices, as suggested by Scardoni et al. (2021). Further analysis of our results with radiative transfer and comparisons with ALMA observations may provide evidence of dust clumping in the observed dust rings.

Second, we have seen the formation of secondary rings/filaments in our weak-bump simulations. These fine-scale substructures are separated on scales $\lesssim H$, but are all located within a single gas pressure bump. We note that we generally do not expect ringlike substructures in the gas component on scales $\lesssim H$ (Dullemond et al. 2018), whereas observational evidence of such fine-scale ringlike substructures is mounting (Jennings et al. 2022), with there being hints of multiple rings within previously identified rings. Our results provide a natural explanation for the formation of multiple sub-$H$-scaled substructures within a single gas bump. For instance, they may provide an alternative explanation for the compact double rings observed in HD 169142 at 57 and 64 au, separately, which have previously been interpreted as the consequence of a migrating protoplanet (Pérez et al. 2019).

Third, we note that the dust rings from our simulations are not entirely axisymmetric, with fine-scale substructures along the azimuth. Although our simulations are local and do not cover the entire azimuth, we notice that fine-scale low-contrast asymmetries have been observed in some systems (e.g., DM Tau; Hashimoto et al. 2021) and may be common among gapped PPDs (van der Marel et al. 2021).

**Disk turbulence.** The dust ring width has recently been used to constrain the level of turbulence (the $\alpha$ parameter) in disks (Dullemond et al. 2018; Rosotti et al. 2020). However, as discussed in Section 4.2, disk turbulence in the presence of a dust ring can be both anisotropic and nonuniform, due to dust feedback.

<footnote n="9">⁹ One might argue that the formation of the secondary dust rings seen in runs Z1w and Z2w could be a consequence of the special timing when we turn on the dust feedback, due to preexisting local density enhancements from the MRI turbulence (e.g., zonal flows). We thus further performed a simulation where the dust feedback was turned on at a much later time $T_{sim} = 1200\Omega_K^{-1}$. We still found similar secondary dust rings/filaments in this case.</footnote>
More specifically, the turbulent diffusion in the radial direction in our simulations is stronger overall than that in the vertical direction, leading to a large aspect ratio of $w_d/H_d \sim 10$ for the single dust ring in run Z2H, which is five times the aspect ratio of the gas bump inserted into the simulation. This effect is even more extreme for the weak-bump case, with $w_d/H_d \sim 30$ for Z2w, due to the formation of multiple dust filaments. This flat morphology of the dust ring is consistent with the recently reported highly settled disk with relatively extended rings in Oph 163131 (Villenave et al. 2022), likely suggesting different levels of turbulence in the radial and vertical directions.

In addition, dust feedback weakens the turbulent diffusion at the dust ring center (see Section 4.2). It makes the dust ring narrower in our strong-bump case, but the presence of additional offcentered filaments increases the overall $w_d$ (if they are unresolved), which largely cancels out the effect of the reduced ring width. Thus, there is a degeneracy between the turbulence level, dust feedback, and bump amplitude in determining the dust ring width. A single $\alpha$ parameter is not necessarily a good indicator of the disk turbulence level, and special care is needed when interpreting the level of disk turbulence from observations of dust rings.

*Prospects for future observations.* Overall, our simulations suggest that the dynamics of dust-trapping rings can be constrained by the presence or absence of fine-scaled substructures. Simultaneously constraining the gas bump and characterizing the fine-scale dust substructures is crucial for determining such expectations. High-resolution observations searching for kinematic features (e.g., ExoALMA)<fnref n="10" /> will provide essential information about the properties of the gas pressure bump, while identifying sub-*H*-scaled substructures will require observations by ALMA with the longest baseline for nearby disks, and also potentially by future upgrades of ALMA (e.g., Burrill et al. 2022), ngVLA (e.g., Andrews et al. 2018b), and/or SKA (Ilee et al. 2020).

## 5. Conclusions and Perspectives

In this paper, we conduct the first controlled experiments of dust dynamics and clumping in MRI-turbulent pressure bumps. Our simulations include a realistic level of MRI turbulence, which is applicable to AD-dominated outer disks, and a carefully designed forcing scheme, to achieve a self-consistent balance between dust trapping and turbulent diffusion. Our main findings include:

1. At a solid abundance compatible to solar, we find robust and efficient dust clumping in radial pressure maxima, with a sizable fraction of the dust mass residing in clumps that are over the Roche density. Secondary dust filaments may form outside the pressure maxima, but they show only weak or no clumping.

2. Dust feedback substantially affects the ring properties. Strong mass loading in dust rings reduces turbulent diffusion, making the ring narrower, with stronger dust settling. Weaker pressure bumps can result in the formation of metastable secondary dust filaments near the bump.

Our work strongly supports the understanding that dust-trapping rings are robust sites for planetesimal formation. Our

simulations do not include self-gravity, and it has yet to be examined how efficiently the dust in dense clumps is converted into planetesimals. This could result from a competition between self-gravity and turbulent diffusion (Gerbig et al. 2020), though the internal turbulent diffusivity is uncertain (but see Klahr & Schreiber 2020, 2021) and requires further study.

Our results also have important observational implications. In particular, depending on the strength of the gas bump, as well as the dust abundance, dust feedback leads to a diversity of dust ring morphologies that may evolve over time. Conventional estimates of the level of turbulence, based on ring width and dust thickness, do not necessarily reflect the true levels of disk turbulence. Moreover, our results suggest the presence of fine-scale substructures, which may be resolvable at the highest resolution by ALMA and the future ngVLA.

As a first study, we only consider a single dust species with a fixed stopping time $\tau_s=0.1$. The natural next step is to incorporate a dust size distribution. This will allow us to further examine how particles of different sizes participate in dust clumping, and to yield more realistic observational signatures. Moreover, the pressure bumps in our shearing box simulations rely on an artificial forcing prescription. Future work should consider more realistic bump formation scenarios, such as magnetic flux concentration and planet-induced gas bumps, in the global context.

We thank Daniel Carrera, Kees Dullemond, Guillaume Laibe, and Chao-Chin Yang for useful discussions, and Greg Herczeg for overseeing the completion of this work. This work is supported by the National Key R&D Program of China, grant No. 2019YFA0405100, and the China Manned Space Project, grant No. CMS-CSST-2021-B09. ZX acknowledges the support of the ERC CoG project PODCAST No. 864965 and NSFC project No. 11773002. The numerical simulations were conducted on the Orion cluster at the Department of Astronomy, Tsinghua University, and on TianHe-1 (A) at the National Supercomputer Center in Tianjin, China.

## Appendix A
### Dust Clump Properties

In order to further analyze the dust properties in dust clumps, we present in Figure A1 the CDF as a function of dust density for runs Z2H and Z2w. The local dust density $\rho_d$ is evenly divided into 100 bins in logarithmic space between 0.01$\rho_0$ and 5000 $\rho_0$, and the CDF is obtained by measuring the probability of a dust particle residing in a region with a local dust density higher than the value of each bin.

The dust density is substantially boosted for both the strong-bump and weak-bump cases when feedback is included, as compared to the black curves showing the corresponding no-feedback runs (Z0 and Z0w). The CDFs are clearly time-dependent, but after about a hundred local orbits, a significant fraction of the dust particles ($\gtrsim 30\%$) are found to reside in regions that are $\gtrsim 200$ times the background gas density, reaching the Roche density. Although self-gravity is not included in our simulations, such a high fraction implies that planetesimal formation could be highly efficient.

For the weak-bump case, we further divide the domain in half at $x = -0.5$, to separate the primary and secondary rings. We see that strong clumping occurs mostly in the primary dust ring at the bump's center ($x > -0.5 H$; see the lower right panel), with the CDF being very similar to that of the

<footnote n="10">¹⁰ https://www.exoalma.com</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure A1.** CDFs of the dust density for runs Z2H and Z2w, as labeled in the titles of each panel. The upper panels show CDFs for the whole simulation domain, whereas the bottom two panels measure the CDFs of run Z2w for dust particles residing in the $x < -0.5$ H and $x > -0.5$ H regions separately, in order to investigate the clump properties in the primary and secondary rings. The CDFs are measured by averaging within different time intervals in units of $\Omega_K^{-1}$, as marked by the different colors (see the legend). The black curves in each panel shows corresponding runs with no feedback, i.e., Z0 and Z0w, averaged from $T = 720\Omega_K^{-1}$ to $1520\Omega_K^{-1}$.</figcaption></figure>

strong-bump case. In the secondary ring ($x < -0.5$ H; lower left panel), the overall dust density is only boosted by a small factor ≲5, with very little clumping. Only a couple of dust clumps are present at a later time, which appear as outliers in the CDF, representing a tiny fraction (≈3%) of particles. The natures of these clumps and their potential to form planetesimals are deserving of future studies.

## Appendix B
### Dust Ring Evolution over the Long Term

In order to examine the lifetime and long-term evolution of the dust rings in our simulation, especially the secondary rings/filaments in the weak-bump case, we continue to run Z1w and Z2w all the way to $T = 2880\Omega_K^{-1}$ ($T_{sim} = 3000\Omega_K^{-1}$), and show the evolution of the dust ring widths and dust ring profiles in Figure B1.

For run Z1w, shown in the upper panel of Figure B1, the overall ring width decreases after $T \sim 1500\Omega_K^{-1}$. This can also be clearly seen in the lower panel of Figure B1, showing the evolution of the dust radial distribution, where the secondary ring spreads and later merges into the primary ring. On the other hand, the secondary ring in run Z2w is more long-lived, persisting over hundreds of orbits. This is likely due to the higher inertia in the dust ring, which withstands potential disruption by turbulent diffusion (see also Appendix C) and/or radial drift. It is also in line with our results in Paper I, where dust feedback enhances the MRI zonal flows with more extended lifetimes. Moreover, we see that additional secondary rings/filaments may form at offcenter positions, and even the primary ring itself has made a small shift in its position over the duration of the long-term simulation (while still residing at the bump's center). These results indicate that the presence of secondary rings/filaments is likely generic, and that they dynamically evolve over a timescale of tens to a few hundreds orbits.
<figure><img src="image_9.png" /><figcaption>**Figure B1.** Upper panel: time evolution of the dust ring width $w_d$ as a function of time $T = 720\Omega_K^{-1}$ to $2880\Omega_K^{-1}$, featuring later stages of the runs Z1w and Z2w. Lower two panels: time evolution of the azimuthally and vertically averaged radial profiles of the dust density at the midplane ($\pm 1$ cell from $z=0$), for runs Z1w and Z2w, for the same time interval as the upper panel.</figcaption></figure>

## Appendix C
### Turbulent Diffusion in Dust Rings

Figures C1 and C2 show the azimuthal- and time-averaged gas rms velocities in the radial and vertical directions, i.e., $\sqrt{\langle v_x^2 \rangle}$ and $\sqrt{\langle v_z^2 \rangle}$, for the strong- and weak-bump cases, respectively. They serve as a proxy for estimating the strength of the turbulent diffusion. For run Z0, without dust feedback, the gas rms velocities in both directions are approximately

uniform, with $\sqrt{\langle v_x^2 \rangle} \sim 0.07c_s$ and $\sqrt{\langle v_z^2 \rangle} \sim 0.016c_s$ (the results from run Z0w are largely the same). With $w_d/w_{bump} = 0.28$ and $H_d/H = 0.067$, our results are consistent overall with the expectations of $w_d/w_{bump} \sim \sqrt{\langle v_x^2 \rangle/(\Omega_K \tau_s)}$ and $H_d/H \sim \sqrt{\langle v_z^2 \rangle/(\Omega_K \tau_s)}$, modulo correction factors from the eddy time (Youdin & Lithwick 2007; Zhu et al. 2015, Paper I). They are also generally consistent with the turbulence parameter $\alpha \sim 10^{-3}$ measured in our simulation.$^{11}$

<footnote n="11">11 The MRI turbulence level in our simulations can be characterized by the $\alpha$ parameter (Shakura & Sunyaev 1973), and is obtained by summing over the Maxwell and Reynold stresses normalized by thermal pressure. Our simulation without dust feedback (Z0) gives time-averaged $\alpha = 1.08 \times 10^{-3}$, consistent with previous studies (Zhu et al. 2015; Xu et al. 2017, Paper I).</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure C1.** The azimuthally averaged spatial distribution of the dust density (the upper panels of both plots), the gas rms velocity in the radial direction (the middle panels), and in the vertical direction (the lower panels), time-averaged within $T = 800–1500\Omega_K^{-1}$, for runs Z0 (left) and Z2H (right). Note that this figure is not to the real aspect ratio.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure C2.** The same as Figure C1, but for the runs Z0w and Z2w.</figcaption></figure>

With feedback included in run Z2H, we find that the radial gas rms velocity is largely reduced in the dust ring, which is clearly visible in Figure C1. We obtain $\sqrt{\langle v_x^2 \rangle} \sim 0.07c_s$ by averaging within the whole domain, which is consistent with the value measured in Z0. In the dust ring, however, we find $\sqrt{\langle v_{x,ring}^2 \rangle} \sim 0.02c_s$, which is about three times lower than that outside the ring. This reduction in $\sqrt{\langle v_x^2 \rangle}$ generally matches the reduction in ring width reported earlier. In the vertical direction, we note that in Paper I we found that dust feedback leads to a strong reduction in the turbulent eddy time, which promotes settling. In this work, we find that, in addition, with more enhanced dust mass loading, the vertical rms velocity $\sqrt{\langle v_z^2 \rangle}$ is also reduced. Although the reduction is less dramatic than that in the radial direction, it is sufficient to account for the additional

settling found in Figure 2. The reduction in both $\sqrt{\langle v_x^2 \rangle}$ and $\sqrt{\langle v_z^2 \rangle}$ is also seen in run Z2w (see the right panels of Figure C2).

ORCID iDs

Ziyan Xu https://orcid.org/0000-0002-2986-8466

Xue-Ning Bai https://orcid.org/0000-0001-6906-9549

References

Andrews, S. M., Terrell, M., Tripathi, A., et al. 2018a, ApJ, 865, 157

Andrews, S. M., Wilner, D. J., Macías, E., Carrasco-González, C., & Isella, A. 2018b, in ASP Conf. Ser., 517, Science with a Next Generation Very Large Array, ed. E. Murphy (San Francisco, CA: ASP), 137

Auffinger, J., & Laibe, G. 2018, MNRAS, 473, 796

Bai, X.-N., & Stone, J. M. 2010a, ApJ, 722, 1437

Bai, X.-N., & Stone, J. M. 2010b, ApJL, 722, L220
Burrill, B. P., Ricci, L., Harter, S. K., Zhang, S., & Zhu, Z. 2022, ApJ, 928, 40

Carrera, D., Johansen, A., & Davies, M. B. 2015, A&A, 579, A43

Carrera, D., & Simon, J. B. 2022, ApJL, 933, L10

Carrera, D., Simon, J. B., Li, R., Kretke, K. A., & Klahr, H. 2021, AJ, 161, 96

Carrera, D., Thomas, A. J., Simon, J. B., et al. 2022, ApJ, 927, 52

Chen, K., & Lin, M.-K. 2020, ApJ, 891, 132

Chiang, E., & Youdin, A. N. 2010, AREPS, 38, 493

Cui, C., & Bai, X.-N. 2021, MNRAS, 507, 1106

Dullemond, C. P., Birnstiel, T., Huang, J., et al. 2018, ApJL, 869, L46

Flaherty, K. M., Hughes, A. M., Rose, S. C., et al. 2017, ApJ, 843, 150

Flaherty, K. M., Hughes, A. M., Teague, R., et al. 2018, ApJ, 856, 117

Gerbig, K., Murray-Clay, R. A., Klahr, H., & Baehr, H. 2020, ApJ, 895, 91

Gole, D. A., Simon, J. B., Li, R., Youdin, A. N., & Armitage, P. J. 2020, ApJ, 904, 132

Hashimoto, J., Muto, T., Dong, R., et al. 2021, ApJ, 911, 5

Heinemann, T., & Papaloizou, J. C. B. 2009, MNRAS, 397, 64

Huang, P., Li, H., Isella, A., et al. 2020, ApJ, 893, 89

Ilee, J. D., Hall, C., Walsh, C., et al. 2020, MNRAS, 498, 5116

Jennings, J., Booth, R. A., Tazzari, M., Clarke, C. J., & Rosotti, G. P. 2022, MNRAS, 509, 2780

Johansen, A., Youdin, A., & Mac Low, M.-M. 2009, ApJL, 704, L75

Klahr, H., & Schreiber, A. 2020, ApJ, 901, 54

Klahr, H., & Schreiber, A. 2021, ApJ, 911, 9

Lehmann, M., & Lin, M. K. 2022, A&A, 658, A156

Li, R., & Youdin, A. 2021, ApJ, 919, 107

Long, F., Pinilla, P., Herczeg, G. J., et al. 2018, ApJ, 869, 17

Manger, N., Klahr, H., Kley, W., & Flock, M. 2020, MNRAS, 499, 1841

Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, MNRAS, 435, 2610

Ono, T., Muto, T., Takeuchi, T., & Nomura, H. 2016, ApJ, 823, 84

Paardekooper, S.-J., Dong, R., Duffell, P., et al. 2022, arXiv:2203.09595

Pérez, S., Casassus, S., Baruteau, C., et al. 2019, AJ, 158, 15

Pinte, C., Dent, W. R. F., Ménard, F., et al. 2016, ApJ, 816, 25

Raettig, N., Lyra, W., & Klahr, H. 2021, ApJ, 913, 92

Riols, A., Lesur, G., & Menard, F. 2020, A&A, 639, A95

Rosotti, G. P., Teague, R., Dullemond, C., Booth, R. A., & Clarke, C. J. 2020, MNRAS, 495, 173

Scardoni, C. E., Booth, R. A., & Clarke, C. J. 2021, MNRAS, 504, 1495

Schäfer, U., Johansen, A., & Banerjee, R. 2020, A&A, 635, A190

Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 500, 33

Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, ApJS, 178, 137

Suriano, S. S., Li, Z.-Y., Krasnopolsky, R., & Shang, H. 2018, MNRAS, 477, 1239

Teague, R., Henning, T., Guilloteau, S., et al. 2018, ApJ, 864, 133

Umurhan, O. M., Estrada, P. R., & Cuzzi, J. N. 2020, ApJ, 895, 4

van der Marel, N., Birnstiel, T., Garufi, A., et al. 2021, AJ, 161, 33

Villenave, M., Stapelfeldt, K. R., Duchêne, G., et al. 2022, ApJ, 930, 11

Whipple, F. L. 1972, in From Plasma to Planet, Proc. of the Twenty-First Nobel Symp., ed. A. Elvius (New York: Wiley Interscience), 211

Xu, Z., & Bai, X.-N. 2022, ApJ, 924, 3

Xu, Z., Bai, X.-N., & Murray-Clay, R. A. 2017, ApJ, 847, 52

Yang, C. C., Johansen, A., & Carrera, D. 2017, A&A, 606, A80

Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459

Youdin, A. N., & Lithwick, Y. 2007, Icar, 192, 588

Zhu, Z., Stone, J. M., & Bai, X.-N. 2015, ApJ, 801, 81