Dust Clumping in Outer Protoplanetary Disks: The Interplay among Four Instabilities

Pinghui Huang; Xue-Ning Bai

doi:10.3847/2041-8213/adcebb

1. Introduction

Planetesimal formation marks the initial stage of planet formation, yet it remains one of the least understood phases (E. Chiang & A. N. Youdin 2010). This is not only due to the difficulty in directly observing the process, but also because of significant uncertainties in the theoretical models. It is widely accepted that planetesimal formation proceeds in a two-stage process: first, dust particles concentrate into dense clumps, and then these clumps undergo gravitational collapse due to their own self-gravity. The first stage is the most critical, where it is generally understood as a result of the streaming instability (SI; J. Goodman & B. Pindor 2000; A. N. Youdin & J. Goodman 2005) or dust trapping in pressure maxima (S. J. Weidenschilling 1980; P. Barge & J. Sommeria 1995). However, both scenarios sensitively depend on the gas dynamics in the protoplanetary disk (PPD).

PPDs are believed to be (weakly) turbulent. This turbulence not only influences the overall structure of the disks but also introduces diffusive effects that strongly affect dust dynamics. In the vertical direction, the turbulence properties determines the vertical profile of dust particles (J. N. Cuzzi et al. 1993; A. N. Youdin & Y. Lithwick 2007), while in the horizontal direction, many gas dynamical processes can lead to the formation of pressure bumps or anticyclonic vortices that can efficiently trap dust (J. Bae et al. 2023). Moreover, the SI is mostly studied under idealized conditions without external turbulence. There has been controversy on the role of external turbulence. On one hand, turbulence is expected to suppress the SI based on linear theory (K. Chen & M.-K. Lin 2020; O. M. Umurhan et al. 2020) and simulation of driven (forced) turbulence (D. A. Gole et al. 2020; J. Lim et al. 2024). However, the situation can be very different under more realistic gas dynamics.

Here we focus on the outer region of PPDs (at a distance of ≳10 au). On one hand, PPDs are weakly ionized, leading to weak coupling between gas and magnetic field, which is mainly governed by ambipolar diffusion (AD) in outer disk conditions. The magnetorotational instability (MRI; S. A. Balbus & J. F. Hawley 1998), a powerful mechanism for driving turbulence in most accretion disks, is expected to be damped by AD, leading to weak-to-modest turbulence (X.-N. Bai & J. M. Stone 2011; J. B. Simon et al. 2013; C. Cui & X.-N. Bai 2021). On the other hand, the vertical shear instability (VSI; R. P. Nelson et al. 2013; A. J. Barker & H. N. Latter 2015) is expected to operate thanks to the relatively short cooling timescale in outer PPDs, generating a moderate level of turbulence. We note that in reality, the MRI and VSI can coexist (C. Cui & X.-N. Bai 2022), with turbulence exhibiting characteristic properties of both. In this work, we focus on the VSI turbulence and its interplay with dust dynamics.

Recently in P. Huang & X.-N. Bai (2025; hereafter HB25), we investigated the dust dynamics in both 2D and 3D VSI turbulence in a suite of global simulations, including dust back-reaction. Our 2D simulations demonstrated that SI can coexist with VSI in axisymmetric setups, where weak pressure bumps generated by 2D VSI turbulence may help form additional dust clumps, consistent with the findings of U. Schäfer et al. (2020), U. Schäfer & A. Johansen (2022), and U. Schäfer et al. (2025). In 3D, the VSI turbulence can produce zonal flows that trigger the Rossby wave instability (RWI; R. V. E. Lovelace et al. 1999; H. Li et al. 2000, 2001), leading to the formation of anticyclonic vortices. With dust, such vortices can substantially enhance dust trapping (S. Richard et al. 2016; N. Manger & H. Klahr 2018; N. Manger et al. 2020). Our study also showed that dust feedback significantly affects 3D VSI turbulence, where dust buoyancy/mass loading suppresses VSI corrugation modes, allowing dust to settle deeper, creating more favorable conditions for the SI to trigger planetesimal formation. However, the resolution in our earlier 3D simulations was insufficient to resolve the SI and capture dust clumping. We expect that with higher resolution, the SI and/or dusty RWI (H. Liu & X.-N. Bai 2023) should operate in the dust layer under 3D VSI turbulence.

In this study, we extend the 3D global simulation of HB25 by incorporating several additional levels of mesh refinement. This allows us to, for the first time, simultaneously resolve the SI while accommodating the VSI turbulence in a 3D global setup. We expect that this setup will enable us to study the interplay of the VSI and RWI with the SI toward a better understanding of how dust clumping and planetesimal formation occur under more realistic gas dynamics in the outer PPDs. Moreover, a thin dust layer can be further subject to the Kelvin–Helmholtz instability (KHI), as a result of strong vertical shear in azimuthal velocity at the interface of the dust layer (M. Sekiya 1998; A. Johansen et al. 2006). This is expected to generate additional stirring to sustain the dust layer thickness. As a result, our simulation allows us to simultaneously observe the VSI, SI, RWI, and KHI.

The outline of this Letter is as follows. We describe the numerical setup in Section 2. In Section 3, we present the results of the 3D global simulation with mesh refinement. In Section 4, we provide the conclusions and further discussions. In HB25, we have already investigated the interplay between VSI and SI, as well as VSI and RWI. For readers who may not be familiar with these topics, we provide a brief introduction about VSI, SI, and RWI in the Appendix.

2. Numerical Setups

We use the multifluid dust module (P. Huang & X.-N. Bai 2022) in the Athena++ code (J. M. Stone et al. 2020) to conduct a 3D global simulation in spherical-polar coordinates (r − θ − ϕ) with mesh refinement for both gas and dust. This simulation explores the interaction between VSI and SI, with RWI and KHI emerging as byproducts. We name this simulation “3D-four instabilities” (3D-FIs) to reflect the interplay among these instabilities.

The computational domain spans r ∈ [1, 3], θ ∈ [π/2 − 0.4, π/2 + 0.4], and ϕ ∈ [0, π]. The setup follows the 2D and 3D simulations described in HB25, using a locally isothermal equation of state, $P={c}_{\,\rm{s}\,}^{2}(R){\rho }_{\,\rm{g}\,}$ , where P and c_s are gas pressure and sound speed, respectively. The radial profiles of sound speed, as well as initial gas and dust densities (ρ_g,init, ρ_d,init), are given by

$\begin{eqnarray}&&{c}_{\,\rm{s}\,}^{2}(R)={c}_{\,\rm{s,0}\,}^{2}{\left(\tfrac{R}{{R}_{0}}\right)}^{-\tfrac{1}{2}},\end{eqnarray} \tag{ 1 }$

$\begin{eqnarray}&&{\rho }_{\,\rm{g,init}\,}(R,z)={\rho }_{0}{\left(\tfrac{R}{{R}_{0}}\right)}^{-\tfrac{3}{2}}\exp \left[\displaystyle \frac{GM}{{c}_{\,\rm{s}\,}^{2}}\left(\displaystyle \frac{1}{r}-\displaystyle \frac{1}{R}\right)\right],\end{eqnarray} \tag{ 2 }$

$\begin{eqnarray}&&{\rho }_{\,\rm{d,init}}=0.01{\rho }_{\rm{g,init}\,},\end{eqnarray} \tag{ 3 }$

where GM ≡ 1 for stellar gravity, and we take ρ₀ ≡ 1, c_s,0 = 0.08, and R₀ ≡ 1 in this simulation. A single dust species with a constant Stokes number St ≡ T_sΩ_K = 0.1 is considered, where T_s is the dust stopping time. Governing equations and initial/boundary conditions are detailed in Sections 2.1 and 2.2 of HB25. We use four stages of refinement to achieve three refinement levels, optimizing memory use and load balancing across computational nodes. The activation timing of mesh refinement does not significantly affect the final results.

We adopt the same “3D-feedback” (3D-FB) setup following Table 1 of HB25 as the root-grid (level 0) configuration for the “3D-FIs” model in this study. The root-grid resolution for both runs is 512 × 256 × 768 cells in the r–θ–ϕ directions. This corresponds to an effective resolution of approximately 40 × 28 × 22 ≃ 29³ cells per ${H}_{\,\rm{g}\,}^{3}$ at r = 1.5. There is no initial difference in the numerical setup between the “3D-FB” and “3D-FIs” models. Dust is initially released, and dust feedback is included in both cases. However, in “3D-FB” of HB25, mesh refinement was not enabled, and thus the SI in “3D-FB” could not be resolved. From 300 to 400 orbits of “3D-FIs,” we activate the first level of mesh refinement (doubling the resolution) to resolve the thin dust layer near the midplane. Subsequent refinement transitions proceed sequentially: from level 1 to level 2, and then from level 2 to level 3. During the transition from level 2 to level 3, mesh refinement first focuses on the region r = 1.2−1.8, followed by refinement in the region r = 1.8−2.4. After full three-level mesh refinement in “3D-FIs,” the effective resolution reaches 318 × 227 × 173 ≃ 232³ cells per ${H}_{\,\rm{g}\,}^{3}$ at r = 1.5 and θ = π/2 (midplane), which is sufficient to capture the unstable modes of SI.

The finest grids of “3D-FIs” cover r ∈ [1.2, 2.4], θ ∈ [π/2 − 0.025, π/2 + 0.025], and ϕ ∈ [0, π]. Each snapshot with full three-level refinement involves 1.2 billion grid points. Due to the computational cost of three refinement stages for “3D-FIs” (∼1.5 × 10⁴ CPU*hours per orbit), we limit it to 100 orbits between 470 and 570 orbits.

For diagnostics, we use cylindrical coordinates (R−ϕ−z) to analyze and present selected results of “3D-FIs,” following the methods in Section 2.3 and Appendix A of HB25. Consistent with HB25, final outcomes remain unaffected by whether dust is introduced initially or after VSI turbulence saturation.

3. Numerical Results and the Coexistence of VSI, SI, RWI, and KHI

In this section, we present the main simulation results. In particular, we closely compare the results of our new “3D-FIs” model with the “3D-FB” model in HB25, following similar analysis procedures and diagnostics, highlighting the new physics brought by the higher resolution thanks to mesh refinement.

3.1. Overview of Simulation Results

We start from Figure 1, which illustrates the overall features of the “3D-FIs” simulation from the snapshot at 550 orbits. The VSI generates strongly anisotropic turbulence throughout the entire simulation domain, characterized by strong vertical motion dominated by the breathing modes (HB25). This leads to stronger dust settling than in 2D, and enhanced dust–gas density ratio at the midplane region. In the meantime, the dusty RWI stems from the zonal flows produced by the VSI, which produces numerous small vortices, many of which are dust clumps that contain high dust–gas density ratios (ρ_d/ρ_g ≳ 100) at the midplane (see Section 3.2).

Figure 1. Refer to the following caption and surrounding text. — **Figure 1.** A 3D schematic diagram of the snapshot at 550 orbits for the “3D-FIs” model. The two vertical slices represent the azimuthally averaged vertical Mach number ${\langle {v}_{\,\rm{g},z}/{c}_{\rm{s}\,}\rangle }_{\phi }$ , the lower and upper semicircles show the dust–gas density ratio at the midplane ${({\rho }_{\,\rm{d}}/{\rho }_{\rm{g}\,})}_{\,\rm{midplane}\,}$ , and the vertical vorticity at the midplane ${({\omega }_{z})}_{\,\rm{midplane}\,}$ , respectively. For visualization purposes, the radial domain is limited to r ∈ [1.2, 2.4].
Download figure:
Standard image High-resolution image

**Figure 1.** A 3D schematic diagram of the snapshot at 550 orbits for the “3D-FIs” model. The two vertical slices represent the azimuthally averaged vertical Mach number ${\langle {v}_{\,\rm{g},z}/{c}_{\rm{s}\,}\rangle }_{\phi }$ , the lower and upper semicircles show the dust–gas density ratio at the midplane ${({\rho }_{\,\rm{d}}/{\rho }_{\rm{g}\,})}_{\,\rm{midplane}\,}$ , and the vertical vorticity at the midplane ${({\omega }_{z})}_{\,\rm{midplane}\,}$ , respectively. For visualization purposes, the radial domain is limited to r ∈ [1.2, 2.4].
Download figure:
Standard image High-resolution image

The overall properties of the VSI turbulence are largely unaffected by the refinement of the dust layer. This can be seen from the left panel of Figure 2, which shows the spacetime (R − t) plot of the vertical Mach number at the midplane, ${\langle M{a}_{z,\,\rm{midplane}\,}\rangle }_{\phi }$ . It closely resembles that of the “3D-FB” model in HB25, which shows “trains” of inward-propagating inertial waves, while the transition between two neighboring wave zones propagates outward at the group velocity of the inertial waves. No abrupt change is observed in the wave patterns as an additional refinement level is introduced.

Figure 2. Refer to the following caption and surrounding text. — **Figure 2.** The spacetime (R − t) plots of the azimuthally averaged vertical Mach number at the midplane (left: ${\langle M{a}_{z,\,\rm{midplane}\,}\rangle }_{\phi }\equiv {\langle {v}_{\,\rm{g},z,\rm{midplane}}/{c}_{\rm{s}\,}\rangle }_{\phi }$ ); and the maxima of the dust–gas density ratio (right: $\max {({\rho }_{\,\rm{d}}/{\rho }_{\rm{g}\,})}_{\phi ,z}$ ) for the “3D-FIs” model are shown. The gray dashed–dotted lines divide the stages where the maximum level in the dust layer is achieved. Level 0 spans 0–300 orbits, level 1 spans 300–400 orbits, level 2 spans 400–470 orbits, and level 3 spans 470–570 orbits.
Download figure:
Standard image High-resolution image

**Figure 2.** The spacetime (R − t) plots of the azimuthally averaged vertical Mach number at the midplane (left: ${\langle M{a}_{z,\,\rm{midplane}\,}\rangle }_{\phi }\equiv {\langle {v}_{\,\rm{g},z,\rm{midplane}}/{c}_{\rm{s}\,}\rangle }_{\phi }$ ); and the maxima of the dust–gas density ratio (right: $\max {({\rho }_{\,\rm{d}}/{\rho }_{\rm{g}\,})}_{\phi ,z}$ ) for the “3D-FIs” model are shown. The gray dashed–dotted lines divide the stages where the maximum level in the dust layer is achieved. Level 0 spans 0–300 orbits, level 1 spans 300–400 orbits, level 2 spans 400–470 orbits, and level 3 spans 470–570 orbits.
Download figure:
Standard image High-resolution image

On the other hand, dust clumping is highly sensitive to resolution, as shown in the right panel of Figure 2, which illustrates the maxima of the dust–gas density ratios, $\max {({\rho }_{\,\rm{d}}/{\rho }_{\rm{g}\,})}_{\phi ,z}$ . The clumping phenomenon occurs right after turning on the first level of mesh refinement, and every further refinement leads to a leap in $\max {({\rho }_{\,\rm{d}}/{\rho }_{\rm{g}\,})}_{\phi ,z}$ without sign of convergence (see further discussion in Section 3.2). All dust clumps slightly migrate inward. According to the linear eigenvalue analysis, the most unstable modes for SI with St = 0.1 exhibit growth rates of $0.1\,\unicode{x02013}\,0.4\ {\Omega }_{\,\rm{K}\,}^{-1}$ at the finest resolution used in this study (see Figure 1 in A. Youdin & A. Johansen 2007). When full three-level mesh refinement is enabled, the effective resolution reaches ≳200³ cells per ${H}_{\,\rm{g}\,}^{3}$ , allowing these fastest-growing modes to be well resolved. As a result, strong dust concentration emerges rapidly, within just a few orbits, consistent with the predicted SI growth rates.

More quantitatively, Figure 3 represents various parameters in the “3D-FIs” simulation, with those from the “3D-FB” model in HB25. These parameters include: the temporal evolution of the Reynolds stress, ${\langle {\alpha }_{R\phi }\rangle }_{R,z\in \pm {H}_{0}}$ (as defined in Equation 14 of HB25), the maxima of the dust–gas density ratios, $\max {({\rho }_{\,\rm{d}}/{\rho }_{\rm{g}\,})}_{R,\phi ,z}$ , and the radially averaged dust–gas scale height ratio, ${\langle {H}_{\,\rm{d}}/{H}_{\rm{g}\,}\rangle }_{R}$ (Equation 22 of HB25). Their initial evolution is identical to that of the “3D-FB” model during the first 300 orbits. When mesh refinement is activated gradually after 300 orbits, stronger Reynolds stress emerges near the midplane, indicating that additional instabilities (e.g., SI, dusty RWI, and KHI) develop within the global VSI turbulence. The progressively higher dust–gas density ratios ( $\max ({\rho }_{\,\rm{d}}/{\rho }_{\rm{g}\,})=100\sim 1000$ ) with increasing resolution suggest that the SI begins to be resolved and generates dust clumps. Although not yet reaching convergence, the dust–gas density ratios within these clumps can well exceed the critical value ${({\rho }_{\,\rm{d}}/{\rho }_{\rm{g}\,})}_{\,\rm{crit}\,}\sim 100$ (Equation 24 in HB25) in typical outer disk conditions after reaching our second finest resolution, potentially allowing dust to gravitationally collapse into planetesimals with the aid of self-gravity. Given that the properties of the VSI turbulence largely remain similar in both “3D-FB” and “3D-FIs,” this reduction in H_d/H_g in “3D-FIs” is associated with dust clumping itself, which mostly occurs in the midplane region.

Figure 3. Refer to the following caption and surrounding text. — **Figure 3.** Top left: temporal evolution of the radially and vertically averaged Reynolds stress, ${\langle {\alpha }_{R\phi }\rangle }_{R,z\in \pm {H}_{0}}$ , where H₀ = 0.08. Top right: temporal evolution of the maximum dust–gas density ratio, $\max {({\rho }_{\,\rm{d}}/{\rho }_{\rm{g}\,})}_{R,\phi ,z}$ . Bottom left: temporal evolution of the radially averaged dust–gas scale height ratio, ${\langle {H}_{\,\rm{d}}/{H}_{\rm{g}\,}\rangle }_{R}$ . Bottom right: cumulative distribution function (CDF) of dust–gas density ratios, computed by counting the number of grid cells, $P[({\rho }_{\,\rm{d}}/{\rho }_{\rm{g}\,})\gt {({\rho }_{\,\rm{d}}/{\rho }_{\rm{g}\,})}_{\,\rm{threshold}\,}]$ . The red solid lines represent the properties of the “3D-FB” model in HB25 without mesh refinement (measured at the same period as we measure the level 0 CDF), while the black solid lines correspond to the “3D-FIs” model. The gray dashed lines indicate the (rough) critical dust–gas density ratios required for planetesimal formation. The blue, green, cyan, and black lines in the CDF plot indicate different refinement stages of the “3D-FIs” model. Radially averaged values are calculated within the range R ∈ [1.2, 2.4]. Similar to Figure 2, the gray dashed–dotted lines indicate the stages with the maximum refinement level.
Download figure:
Standard image High-resolution image

**Figure 3.** Top left: temporal evolution of the radially and vertically averaged Reynolds stress, ${\langle {\alpha }_{R\phi }\rangle }_{R,z\in \pm {H}_{0}}$ , where H₀ = 0.08. Top right: temporal evolution of the maximum dust–gas density ratio, $\max {({\rho }_{\,\rm{d}}/{\rho }_{\rm{g}\,})}_{R,\phi ,z}$ . Bottom left: temporal evolution of the radially averaged dust–gas scale height ratio, ${\langle {H}_{\,\rm{d}}/{H}_{\rm{g}\,}\rangle }_{R}$ . Bottom right: cumulative distribution function (CDF) of dust–gas density ratios, computed by counting the number of grid cells, $P[({\rho }_{\,\rm{d}}/{\rho }_{\rm{g}\,})\gt {({\rho }_{\,\rm{d}}/{\rho }_{\rm{g}\,})}_{\,\rm{threshold}\,}]$ . The red solid lines represent the properties of the “3D-FB” model in HB25 without mesh refinement (measured at the same period as we measure the level 0 CDF), while the black solid lines correspond to the “3D-FIs” model. The gray dashed lines indicate the (rough) critical dust–gas density ratios required for planetesimal formation. The blue, green, cyan, and black lines in the CDF plot indicate different refinement stages of the “3D-FIs” model. Radially averaged values are calculated within the range R ∈ [1.2, 2.4]. Similar to Figure 2, the gray dashed–dotted lines indicate the stages with the maximum refinement level.
Download figure:
Standard image High-resolution image

3.2. Dust Clumps and Vortices

In this subsection, we analyze in further detail the dust dynamics and the properties of dust clumping. In addition to the interplay between VSI and (dusty) RWI as studied in HB25, we highlight how the picture is enriched by the SI.

We start by looking at the horizontal and vertical distributions of the dust–gas density ratio, as shown in Figure 4 at the ϕ = π/2 and θ = π/2 (z = 0, midplane) slices over different times. From the horizontal distributions, we see that as the grid gets refined, we clearly identify the formation of dust clumps. Such clumps become much more densely populated toward higher resolution, and the individual clumps get more and more compact. From the vertical distribution, we see that most of the dust settles into very thin layers. There are regions where dust shows larger scale heights with corrugation patterns. They are associated with the transition between two neighboring inertial wave zones that propagate at the group velocity, as identified and discussed in Figure 10 of HB25. Interestingly, this region shows very little dust clumping even with much increased resolution.

Figure 4. Refer to the following caption and surrounding text. — **Figure 4.** The dust–gas density ratios ρ_d/ρ_g on the R–z (ϕ = π/2) section (top four panels) and R–ϕ (θ = π/2, z = 0, midplane) section (bottom four panels) for the “3D-FIs” model at various times are shown. From left to right, from top to bottom in the top four panels, the snapshots are taken at 300, 400, 450, and 550 orbits. The bottom four panels correspond to the same times as the top four panels. Different times and the maximum mesh levels are marked in the top-left corners of the top four panels. The black dashed rectangles indicate the zoomed-in regions shown in Figure 5.
Download figure:
Standard image High-resolution image

The dust clumps are associated with dusty RWI vortices. In Figure 5, we show zoomed-in views of the vertical vorticity and dust–gas density ratios at the midplane. Before activating mesh refinement (prior to 300 orbits), the dust–gas density ratio profiles resemble those of the “3D-FB” model. As refinement progresses, narrower dust rings form, which eventually break into smaller dust clumps, in accordance with the development of the dusty RWI (H. Liu & X.-N. Bai 2023). Indeed, these dust rings and clumps are associated with regions of negative vertical vorticity, indicating that they correspond to dust concentrations within tiny dusty vortices (mild pressure bumps). We note that there is no one-to-one correspondence between gas vortices and dust clumps; there are, in general, many more gas vortices. Moreover, only a small fraction (1% ∼ 2%) of dust clumps have density ratios exceeding the critical value, as highlighted with white contours.

Figure 5. Refer to the following caption and surrounding text. — **Figure 5.** The zoomed-in plots of gas vertical vorticity at the midplane ${({\omega }_{z})}_{\,\rm{midplane}\,}$ (top four panels), and dust–gas density ratios at the midplane ${({\rho }_{d}/{\rho }_{g})}_{\,\rm{midplane}\,}$ (bottom four panels) at various snapshots. The order of snapshots is the same as Figure 4. The white spots and contours indicate the regions where the local dust–gas density ratios exceed the critical ratio (ρ_d/ρ_g)_crit = 98 (Equation (24) in HB25).
Download figure:
Standard image High-resolution image

**Figure 5.** The zoomed-in plots of gas vertical vorticity at the midplane ${({\omega }_{z})}_{\,\rm{midplane}\,}$ (top four panels), and dust–gas density ratios at the midplane ${({\rho }_{d}/{\rho }_{g})}_{\,\rm{midplane}\,}$ (bottom four panels) at various snapshots. The order of snapshots is the same as Figure 4. The white spots and contours indicate the regions where the local dust–gas density ratios exceed the critical ratio (ρ_d/ρ_g)_crit = 98 (Equation (24) in HB25).
Download figure:
Standard image High-resolution image

The vortices show vertical structures. Figure 6 displays slices at constant θ of the vertical vorticity ω_z at 550 orbits. The midplane dust layer is the most densely populated with tiny vortices. The vortex structures at 1 ∼ 2H_g above the midplane are similar, showing small vortices similar to those reported in HB25, while the vorticity increases over height. At higher altitudes (θ = π/2 + 4H₀/R₀), density waves induced by VSI surface modes and larger vortices are observed.

Figure 6. Refer to the following caption and surrounding text. — **Figure 6.** The vertical vorticity ω_z at various θ-angles is shown at the snapshot of 550 orbits for the “3D-FIs” model, similar to Figure 12 in HB25. We select θ = π/2 (z = 0, midplane), π/2 + H₀/R₀, π/2 + 2H₀/R₀, and π/2 + 4H₀/R₀, where H₀/R₀ = 0.08 in this simulation. Note that the mesh level at θ = π/2 is 3, while the mesh levels at other θ-angles are 0 (root-grid level).
Download figure:
Standard image High-resolution image

Also shown in Figure 3 are the cumulative distribution functions (CDFs) of the dust–gas density ratio exceeding specific threshold values, $P[({\rho }_{\rm{d}}/{\rho }_{\rm{g}})\gt {({\rho }_{\rm{d}}/{\rho }_{\rm{g}})}_{\rm{threshold}}]$ , obtained by counting the number of cells. At level 0, the CDF agrees with the “3D-FB” model presented in HB25.⁵ With full three-level refinement, the tail of the CDF extends to ρ_d/ρ_g ≳ 1000. At this stage, the ratio of the mass of dust clumps with $({\rho }_{\,\rm{d}}/{\rho }_{\rm{g}\,})\gt {({\rho }_{\,\rm{d}}/{\rho }_{\rm{g}\,})}_{\,\rm{crit}\,}$ to the total dust mass is approximately M_clumps,total/M_dust,total ∼ 0.02–0.03. Each clump contains roughly 10⁻⁶ ∼ 10⁻⁵ of the total dust mass. Despite reaching a very high resolution sufficient to resolve the SI, we do not find a sign of convergence of the CDF with increasing resolution. Note that even for pure SI, convergence in dust clumping behavior is known to require very high resolution, especially in 3D. Our study represents a step forward as the first study to investigate dust clumping under a self-consistent turbulent environment in a 3D global setup. Our results here mainly suggest that dust clumping and planetesimal formation can be more efficient than quoted above. Our results also call for independent verification from particle-based code, as they may show different convergence properties in certain problems (L. Krapp et al. 2019; P. Huang & X.-N. Bai 2022).

Based on the analysis above, we now discuss the role of SI on dust clumping. We first note that zonal flows from the VSI are not strong enough to overcome the background pressure gradient in the PPD (as discussed in HB25), thus, the condition is favorable to promote the SI and dust clumping (e.g., X.-N. Bai & J. M. Stone 2010a). From HB25, with low resolution, we also found that dust feedback already substantially enhances dust concentration. On the other hand, the association of dust clumps with vortices is suggestive of the dusty RWI (H. Liu & X.-N. Bai 2023), although existing studies assume a pressure maxima and are only 2D. We thus favor the interpretation of a combination of the SI and dusty RWI operating under such conditions. Finally, given that vortices are all very small and get smaller with increasing resolution, we do not expect to see the operation of the SI within vortices with migrating dust.

3.3. KHI

With dust settling into a thin layer in the disk midplane, the gas azimuthal velocity is accelerated toward Keplerian by the strong dust feedback, while the gas above and below the dust layer is largely unaffected and still moves at non-Keplerian speed by the global (usually negative) pressure gradient. This gas azimuthal velocity transition leads to strong vertical shear in v_g,ϕ and can be subject to the KHI, which develops at the azimuthal–vertical (ϕ−z) plane (M. Sekiya 1998; A. N. Youdin & F. H. Shu 2002; E. Chiang 2008; X.-N. Bai & J. M. Stone 2010a; A. T. Lee et al. 2010a, 2010b; K. Gerbig et al. 2020; D. Sengupta & O. M. Umurhan 2023). One useful criterion for the onset of KHI is the Richardson number Ri, which can be defined as (E. Chiang 2008; X.-N. Bai & J. M. Stone 2010a)

$\begin{eqnarray}&&Ri\equiv -\displaystyle \frac{GMz}{{\rho }_{\,\rm{eff}}{r}^{3}}\displaystyle \frac{\partial {\rho }_{\rm{eff}\,}/\partial z}{{(\partial {v}_{\,\rm{g}\,{,}\phi }/\partial z)}^{2}}.\end{eqnarray} \tag{ 4 }$

For tightly coupled dust (with a Stokes number St ≲ 0.1), the effective density is given by ρ_eff = ρ_d + ρ_g. In the absence of rotation, it is known that the KHI develops for Ri < 0.25 (J. N. Cuzzi et al. 1993; M. Sekiya 1998; A. N. Youdin & F. H. Shu 2002; A. Johansen et al. 2006; A. T. Lee et al. 2010a, 2010b). When considering the stabilizing effects of the Coriolis force, E. Chiang (2008) suggested that empirically, the critical value of Ri is about 0.1, which we adopt here.

To assess the potential activation of the KHI, we perform both a visual inspection of the dust layer at specific radial locations (r = 1.3, 1.5) based on data at 550 orbits and a more quantitative analysis of the Richardson number (Ri) profile at the corresponding positions. The results are presented in Figure 7.

Figure 7. Refer to the following caption and surrounding text. — **Figure 7.** The top two panels display the radial profiles of the dust–gas density ratio for the “3D-FIs” model at 550 orbits, taken at r = 1.3 and 1.5. The middle two panels show the corresponding radial–azimuthal profiles of the Richardson number Ri (Equation (4)), where we apply a moving-average with window width 0.5H_g along the azimuthal direction (indicated by the black solid segment in the bottom-right corners). White contours in the top four panels outline regions where the dust–gas density ratio exceeds unity. The bottom panel presents vertical profiles of the azimuthally averaged values of Ri and the minimum Ri from the moving average at the two radial locations.
Download figure:
Standard image High-resolution image

The dust layer does exhibit distinct wave-like patterns in the azimuthal–vertical (ϕ−z) plane at different radii, reminiscent of the activation of the KHI. When we compute the Richardson number profile by averaging data over the entire azimuthal domain, we observe that the Richardson number remains large (Ri ≫ 0.1) beyond the dust layer, as expected. While it rapidly decreases toward the dust layer, the values remain well above 0.1 (see the bottom panel of Figure 7). However, noticing that dust across the azimuth does not necessarily communicate, it is more reasonable to measure the local Richardson number through a moving window of finite width along the ϕ direction. When taking the window size to be 0.5 H_g, which is more than 10 times the mean dust layer thickness, we see that the Richardson number can go below the critical value of 0.1 (the green and yellow regions). Interestingly, most such regions coincide with the edges of the white contours, where the dust–gas density ratio is unity. Such contours can be considered as a marker for the dust–gas interface, suggesting that the KHI can indeed be triggered near such interfaces.

Interestingly, the dust density at the dust–gas interface exhibiting Ri < 0.1 is systematically lower (ρ_d/ρ_g ≲ 1), and no dust clumps are observed in these regions. To further elucidate the role of KHI, we notice from Figure 2 that when averaged over azimuth, the dust “rings” persist over hundreds of orbits. On the other hand, dust in these narrow rings is spread into many small clumps within individual dusty vortices. Further examination of Figure 5 and neighboring snapshots reveals that the azimuthal dust distribution varies substantially from snapshot to snapshot (separated by one innermost orbit). This suggests that many of the clumps likely form and dissolve in transient, potentially due to the stirring (at least in part by) the KHI, although more detailed investigation is beyond the scope of this work.

4. Conclusions and Discussion

In this work, we conduct the first 3D global simulation of dust dynamics in PPDs subject to VSI with sufficient resolution to resolve the SI. With an initial dust–gas mass ratio of 0.01, we observe strong dust clumping and the complex interplay among the VSI, RWI, SI, and KHI.

1.
Among the four instabilities, the VSI sets the bulk level of turbulence and governs the overall gas dynamics, despite that the corrugation mode is suppressed, leaving the breathing mode as the dominant dust stirrer.
2.
The VSI generates weak zonal flows that do not necessarily form pressure maxima, allowing the SI to develop. Meanwhile, zonal flows facilitate the development of dusty RWI, leading to the formation of small dust clumps inside RWI vortices. Dust clumping is likely the outcome of both instabilities.
3.
The midplane layer is characterized by densely populated small gas vortices, though only a small fraction (1% ∼ 2%) of such vortices host dust clumps exceeding the Roche density.
4.
With increasing resolution, more dust participates in dust clumping, forming higher-density clumps that do not show evidence of numerical convergence. About 2% ∼ 3% of dust resides in regions of clumps subject to gravitational collapse at our highest resolution.
5.
The dust layer is further affected by the KHI, which develops in localized azimuthal regions. It likely produces wave-like structures at the dust–gas interface along the azimuth, and contributes to the temporary disruption of dust clumps.

Our work demonstrated that when accounting for dust feedback, dust clumping, likely triggered by the SI, can persist on top of 3D VSI, RWI, as well as KHI turbulence. It suggests that under a realistic environment, turbulence is not necessarily a diffusive and obstructive process for dust concentration, and dust back-reaction plays a crucial role in assisting dust clumping. On the other hand, recall that dust clumping under the MRI turbulence likely occurs only in pressure maxima in addition to dust back-reaction (Z. Xu & X.-N. Bai 2022), a condition different from what we observe here. Altogether, these works reveal the rich physics in dust–gas interaction and highlight the importance of self-consistent treatment of the gas turbulent environment, which is essential for a comprehensive understanding of dust dynamics and planetesimal formation.

Due to limited computational resources, we conducted only one single 3D global simulation in this Letter with a set of parameters that we consider typical for general PPDs. Further investigations should explore more parameters such as dust size and abundance, as well as variations in the disk model.

We made several simplifying assumptions in terms of the physics incorporated, leaving room for future improvements. First, we simplify thermodynamics by fixing disk temperature with a locally isothermal equation of state. More self-consistent treatment of radiation transport can alter the dominant mode of the VSI (J. D. Melon Fuksman et al. 2024; S. Zhang et al. 2024), which, in the meantime, may also trigger the convective overstability (H. Klahr 2024). Intermediate cooling rates could also affect SI growth rates (M. Lehmann & M.-K. Lin 2023) and the decay of RWI vortices (J. Fung & T. Ono 2021). Second, we did not incorporate the magnetic field, which likely dominates disk angular momentum transport by launching MHD winds, together with the MRI turbulence (C. Cui & X.-N. Bai 2021). The MRI can coexist with the VSI depending on disk ionization (C. Cui & X.-N. Bai 2022). The disk wind configuration is subject to magnetic flux concentration to drive strong zonal flows (X.-N. Bai & J. M. Stone 2014; S. S. Suriano et al. 2018) and generate the RWI on its own (C.-Y. Hsu et al. 2024). Third, self-gravity was neglected, but it becomes critical for understanding the initial mass function of planetesimals (A. Johansen et al. 2015; J. B. Simon et al. 2016). Finally, this study considered a single dust species, but dust size distribution can significantly affect disks' cooling and hence the VSI turbulence (Y. Fukuhara & S. Okuzumi 2024), as well as dust trapping in RWI vortices (Y.-P. Li et al. 2020), and potentially the properties of the SI as well (L. Krapp et al. 2019; C.-C. Yang & Z. Zhu 2021; Z. Zhu & C.-C. Yang 2021; J. Matthijsse et al. 2025). It should also be noted that dust coagulation can proceed efficiently in dust clumps, promoting further dust clumping and planetesimal formation (R. T. Tominaga & H. Tanaka 2023; K. W. Ho et al. 2024). All of these effects are potentially fruitful avenues for future research.

Acknowledgments

We thank Orkan Umurhan, Andrew Youdin, Hui Li, Eugene Chiang, Min-Kai Lin, Chao-Chin Yang, Shengtai Li, Can Cui, Rixin Li, and Shangjia Zhang for helpful discussions. This work is supported by the National Science Foundation of China under grant Nos. 12233004, 12325304, and 12342501, and the China Manned Space Project with No. CMS-CSST-2021-B09. P.H. acknowledges the Canadian Grant NSERC ALLRP 577027-22.

Software: Athena++ (J. M. Stone et al. 2020; P. Huang & X.-N. Bai 2022).

Appendix: Theoretical Background

A.1. VSI

In PPDs, the density profile is primarily set by angular momentum transport and stellar gravity, while the temperature profile is predominantly regulated by stellar irradiation (E. I. Chiang & P. Goldreich 1997). Notably, the contours of constant density and pressure are generally not parallel, making PPDs essentially nonbarotropic (∇ρ_g × ∇P ≠ 0). This leads to the presence of weak vertical shear (A. J. Barker & H. N. Latter 2015):

$\begin{eqnarray}&&2R{{\rm{\Omega }}}_{\rm{g}}\displaystyle \frac{{\rm{\partial }}{{\rm{\Omega }}}_{\rm{g}}}{{\rm{\partial }}z}=-\displaystyle \frac{\left({\rm{\nabla }}{\rho }_{\rm{g}}\times {\rm{\nabla }}P\right)}{{\rho }_{\rm{g}\,}^{2}}\cdot \hat{\phi }.\end{eqnarray} \tag{ A1 }$

When subjected to rapid thermal cooling (the thermal relaxation timescale is much smaller than Keplerian dynamic timescale, ${t}_{\,\rm{cool}}\ll {\Omega }_{\rm{K}\,}^{-1}$ ), the strong stabilizing vertical buoyancy within disks is diminished (M.-K. Lin & A. N. Youdin 2015). This triggers the VSI (V. Urpin & A. Brandenburg 1998; V. Urpin 2003; R. Arlt & V. Urpin 2004; R. P. Nelson et al. 2013), which is an axisymmetric (R–z plane) and centrifugal instability analogous to the Goldreich–Schubert–Fricke instability in differentially rotating stars (P. Goldreich & G. Schubert 1967; K. Fricke 1968).

The VSI taps the free energy from vertical shear in disks and generates modest and anisotropic turbulence characterized by the “body modes.” They result from the destabilization of the inertial waves and are usually dominated by strong vertical oscillation, which levitates dust particles (M. H. R. Stoll & W. Kley 2016; M. Flock et al. 2020; C. P. Dullemond et al. 2022). Furthermore, VSI can generate weak zonal flows and vortices from secondary instabilities (RWI), facilitating dust trapping/concentration (S. Richard et al. 2016; M. Flock et al. 2017; N. Manger & H. Klahr 2018; M. Flock et al. 2020; N. Manger et al. 2020).

A.2. SI

In PPDs, dust particles undergo radial drift toward higher gas pressure (S. J. Weidenschilling 1977). When considering the dust back-reaction to the gas, this motion becomes linearly unstable to axisymmetric perturbations, leading to the SI (J. Goodman & B. Pindor 2000; A. N. Youdin & J. Goodman 2005; A. Youdin & A. Johansen 2007). It can also be interpreted as a resonant drag instability (P. F. Hopkins & J. Squire 2018; J. Squire & P. F. Hopkins 2018, 2020), resulting from the resonance between the dust streaming motion (drift velocity w) and gas epicyclic motion (epicyclic frequency κ), i.e., ${\boldsymbol{w}}\cdot \hat{k}=\kappa$ , where $\hat{k}$ is the projected wavevector. The free energy of SI comes from the streaming motion due to radial pressure gradient, characterized by

$\begin{eqnarray}&&\eta \equiv \displaystyle \frac{1}{2}{\left(\displaystyle \frac{{c}_{\,\rm{s}}}{{v}_{\rm{K}\,}}\right)}^{2}\displaystyle \frac{d\,{\mathrm{ln}}\,P}{d\,{\mathrm{ln}}\,R},\end{eqnarray} \tag{ A2 }$

and it is often nondimensionized by Π ≡ ηv_K/c_s, where c_s and v_K are the gas sound speed and Keplerian speed.

In the nonlinear regime, gas drag, dust feedback, and the Coriolis force can create a positive loop that amplifies dust density fluctuations, leading to strong dust overdensities (A. Johansen & A. Youdin 2007; J. Squire & P. F. Hopkins 2018; N. Magnan et al. 2024). Considering vertical stratification, the SI can yield significant dust clumping depending on the dust abundance, size distribution and radial pressure gradient (X.-N. Bai & J. M. Stone 2010a, 2010b; D. Carrera et al. 2015; R. Li & A. N. Youdin 2021; J. Lim et al. 2025; D. Ostertag & M. Flock 2025), and when further incorporating self-gravity, such dust clumps can further collapse to directly form planetesimals (A. Johansen et al. 2007; J. B. Simon et al. 2016; R. Li et al. 2018).

The interplay between VSI and SI has recently been investigated in U. Schäfer et al. (2020), U. Schäfer & A. Johansen (2022), U. Schäfer et al. (2025), as well as our work (HB25), in a 2D axisymmetric setting. It was found that SI exhibits the capacity to generate more robust dust clumping when compared to scenarios where SI acts alone.

A.3. RWI

The RWI is a local nonaxisymmetric instability (R. V. E. Lovelace et al. 1999; H. Li et al. 2000; T. Ono et al. 2016), which occurs in regions where there is a local maximum of the key function ${\mathscr{L}}$ known as “entropy modified inversed vortensity”:

$\begin{eqnarray}&&{\mathscr{L}}\equiv \displaystyle \frac{{\Sigma }_{\,\rm{g}\,}{s}^{2/\gamma }}{2\left({\rm{\nabla }}\times {{\boldsymbol{v}}}_{\,\rm{g}\,}\right)\cdot \hat{z}},\end{eqnarray} \tag{ A3 }$

where Σ_g, γ, v_g, and s are gas surface density, adiabatic index, gas velocities, and specific gas entropy, respectively. The specific entropy is defined as $s\equiv {p}_{\,\rm{g}}/{\Sigma }_{\rm{g}\,}^{\gamma }$ , and p_g is the vertical integrated pressure. A local maximum of ${\mathscr{L}}$ could occur in locations characterized by anomalous radial shear at orbital frequency ∂Ω_g/∂R compared to background Keplerian shear, typically associated with strong radial pressure variations. Examples of such locations include gap edges created by planets (M. de Val-Borro et al. 2007; M.-K. Lin & J. C. B. Papaloizou 2011; Z. Zhu et al. 2014), zonal flows (A. Johansen et al. 2009; S. Richard et al. 2016), mass infall (J. Bae et al. 2015; A. Kuznetsova et al. 2022), boundary layer (Z. Fu et al. 2023) and dusty rings edges (W. Fu et al. 2014; A. Pierens et al. 2019; H.-F. Hsieh & M.-K. Lin 2020; P. Huang et al. 2020; C.-C. Yang & Z. Zhu 2020; H. Liu & X.-N. Bai 2023; K. Chan & S.-J. Paardekooper 2024). Since the anomalous radial shear is important, RWI can be considered as a special Kelvin–Helmholtz type instability acting on the horizontal (R−ϕ) plane for accretion disks.

The nonlinear development of the RWI in disks usually leads to the production of multiple anticyclonic vortices, which may eventually merge into one big vortex (H. Li et al. 2001; H. Meheut et al. 2012a, 2012b, 2012c; T. Ono et al. 2018). The vortices can launch density waves (G. Bodo et al. 2005; H. Meheut et al. 2010; P. Huang et al. 2019; X. Ma et al. 2025), which lead to their migration (H. Li et al. 2001; S.-J. Paardekooper et al. 2010; T. Ono et al. 2018). A most attractive property of RWI vortices is that their anticyclonic nature is characterized by a pressure maxima at the center, thus they are expected to be able to concentrate dust particles and facilitate planetesimal formation (P. Barge & J. Sommeria 1995). Recent observations made by the Atacama Large Millimeter/submillimeter Array indicate that azimuthal asymmetries in the dust continuum are not uncommon and may be attributed to the trapping of dust within RWI vortices (N. van der Marel et al. 2013, 2016; R. Dong et al. 2018; H. Yang et al. 2023).