© The Authors 2026

1 Introduction

Streaming instability (Youdin & Goodman 2005; see Squire & Hopkins 2020 for a review) is the model that is currently favored with respect to the formation of planetesimals in a protoplanetary disk. This instability arises from the relative drift between dust particles (mm or cm in size) and gas with mutual aerodynamic coupling, concentrating dust into filaments during a linear growth phase, where dust clumps can then form via a nonlinear growth phase (Abod et al. 2019; Li & Youdin 2021). When the dust density in a clump reaches the Roche density, the clump collapses under its own self-gravity, forming 10–100 km sized bodies, known as planetesimals (Klahr & Schreiber 2020).

However, there is evidence to suggest that the radial drift of dust relative to gas is limited to a fraction of the disk’s lifetime, likely within the first million years (Myr). Indeed, extrasolar protoplanetary disks remain rich in dust throughout their lifetime of ~5 Myr and their dust radius does not seem to shrink over time (Najita & Bergin 2018), suggesting that dust radial drift has been blocked (Birnstiel 2024). Moreover, the Solar System presents two generations of planetesimals. The first generation, including the parent bodies of iron meteorites and other achondrites, did form in the first Myr (Spitzer et al. 2021); however, the parent bodies of chondritic meteorites formed much later, namely, between 2 and 4 Myr after the birth of the protoplanetary disk (e.g., Neumann et al. 2024). Despite this, there are pairs of achondrites-chondrites (e.g., aubrites and enstatite chondrites, and the achondrites NWA 6704 and NWA 011 together with CR chondrites) that have identical nucleosynthetic isotopic anomalies. Given the disk was radially heterogeneous (e.g., Kleine et al. 2020), this implies that dust remained blocked and confined during the time lapse separating the two planetesimal formation events. Thus, the very framework of the streaming instability (dust radial drift) most likely does not apply to planetesimals forming in a late disk.

An efficient mechanism for trapping dust in the disk is a pressure bump. A local maximum in an otherwise monotonically decreasing radial pressure profile of the gas acts as a dust trap by halting the inward drift and allowing dust grains to accumulate. Protoplanetary disk observations by the Atacama Large Millimeter/submillimeter Array (ALMA) have revealed that concentric dust rings are a ubiquitous feature among these observed disks (see, e.g., Andrews et al. 2018) and are nearly all confined within pressure maxima (Stadler et al. 2025). The question is then whether these dust rings can also act as the nursery for dust clumps and be the eventual birth place of planetesimals. To answer this question, we must consider the dust-gas dynamics within the dust rings and the hydro-instabilities that may be active and capable of concentrating dust into denser regions.

Although the dust is effectively trapped, the formation of these rings does not immediately trigger planetesimal formation. Turbulence in the gas prevents dust from settling sufficiently toward the midplane, a necessary condition for gravitational collapse. Even modest turbulence levels (with an α-viscosity around l0⁻⁴) can inhibit efficient sedimentation. The same is true for the radial concentration at the pressure maximum (Dullemond et al. 2018). Thus, a large amount of dust might need to accumulate over time into the ring, before planetesimal formation becomes possible.

This highlights the importance of modeling dust rings with realistic levels of turbulence in order to fully understand how hydrodynamic instabilities behave within and near pressure bumps. While pressure bumps can effectively trap dust and suppress classical instabilities like the Kelvin-Helmholtz instability at their centers, the overall picture remains complex. The limited effectiveness of the streaming instability—particularly at the pressure maximum where relative dust-gas drift vanishes—suggests that additional mechanisms are needed to form dust clumps in these regions (see, e.g., Auffinger & Laibe 2018; Carrera & Simon 2022). This motivates the exploration of alternative instabilities that might operate under realistic disk conditions and potentially drive dust clumping within the dust rings. One such candidate is the dusty Rossby wave instability (DRWI), first coined by Liu & Bai (2023), and investigated in Cui et al. (2025).

These authors identified two distinct modes of the DRWI. Type I is similar to the classic Rossby wave instability (RWI), but involving dust as well. This mode leads to the formation of a single large gas vortex that traps and concentrates dust—a process that has been well studied. However, the viability of type I DRWI faces several challenges in realistic disk conditions.

First, the classic RWI requires a relatively sharp pressure bump, whereas milder pressure bumps are expected to be more common in protoplanetary disks. Second, as previously mentioned, observations show that dust rings are a ubiquitous feature of these disks, while crescent-shaped asymmetries, which are thought to be signatures of gas vortices, are comparatively rare. The widespread presence of dust rings also suggests that they are long-lived structures. In contrast, as type I DRWI leads to vortex formation, it disrupts and destroys the dust ring on timescales that are much shorter than the lifetime of dust rings inferred from observations.

This motivates a more in-depth examination of the newly identified type II DRWI, which presents a promising alternative pathway to planetesimal formation. One key advantage is that it can operate in relatively mild pressure bumps (including many that remain stable to the standard RWI), making it a more broadly applicable and plausible scenario. Additionally, this instability has the potential to concentrate dust into dense clumps, which may eventually form planetesimals, while preserving the overall dust ring structure.

In this work, after presenting our computational methods (Section 2) and disk model (Section 3), we reproduce the results of Liu & Bai (2023) on the type II DRWI and extend them to longer timescales (Section 4). We then extend the analysis to 3D, beginning with a viscous disk with α = 10⁻⁴ (Section 5), followed by an inviscid disk (Section 6). In the inviscid case, we consider a setup with a finite dust reservoir (no inflow) and a scenario with continuous dust inflow toward the pressure bump, with the latter investigated using a high-resolution 2D r – z simulation. We show that it is an inefficient mechanism for increasing the local dust density, particularly in the inviscid regime. Moreover, as discussed above, sustained dust inflow over the disk lifetime is unlikely. For these reasons, in Section 7, we study an alternative mechanism for planetesimal formation in a pressure bump (i.e., during gas removal). In particular, we investigate how the maximal dust density should scale with the overall dust-to-gas mass ratio in the ring, which should eventually lead to a gravitational instability of the dust clumps as gas is removed. The conclusions follow in Section 8.

2 Numerical methods

We considered a global protoplanetary disk model treated as a non-self-gravitating mixture of gas and dust whose motion is described by the Navier–Stokes equations. We used spherical coordinates (r, φ, θ), where r is the radial distance from the host star of mass M_* (which is at the origin of the coordinate system), φ the polar angle measured from the z-axis (the colatitude), and θ is the azimuthal coordinate measured from the x-axis. The midplane of the disk is located at the equator $\varphi =\frac{\pi }{2}$. The gravitational potential of the central star acting on the disk is Φ = −GM_*/r. We did not consider indirect terms that arise from the primary acceleration due to disk’s gravity (Crida et al. 2025).

We used the fargOCA (Lega et al. 2014) code, which is a grid-based code that solves the hydrodynamic equations with a time-explicit method, using operator splitting and upwind techniques. The code has been recently reorganized to run both on CPUs and GPUs with the addition of dust dynamics (Miniussi et al., in prep.).

We model the gas as a fluid with volume density, ρ_g, and velocity, u = (u_r, u_θ, uφ), with u_θ = r sin(ϕ)(ω + Ω_f), where ω is the azimuthal angular velocity in a frame rotating around the z-axis at angular velocity, Ω_f. The dust is modeled as a nonviscous pressureless fluid described by a volume density, ρ_d, and a velocity, v. Dust experiences a drag force from the collisions with gas that can be expressed (per unit volume) as $$F{d}_{}=-{\rho }_d\frac{{\Omega }_k}{S_t}\left(v-u\right),$$

where Ω_k/S_t is the friction timescale, which characterizes the time needed for dust grains with the Stokes number, S_t, to adjust their velocity to a change of gas velocity. The gas feels a reaction force or feedback of sign opposite to drag force exerted on the dust, which is proportional to the amount of dust, $$F_g=-{\rho }_d\frac{{\Omega }_k}{S_t}\left(v-u\right).$$

The continuity equation for the gas is expressed as $$\frac{\partial {\rho }_g}{\partial t}+\nabla \cdot \left({\rho }_{g}u\right)=0.$$

The Navier-Stokes equations for the radial momentum, Y_r = ρ_gu_r, the polar momentum, Υ_φ = ρ_gru_φ, and the angular momentum, Υ_θ = ρ_gr sin (φ)u_θ, are expressed as $$\{\begin{array}{lll}\frac{\partial Y_r}{\partial t}+\nabla \cdot \left(Y_{r}u\right)\hfill & =\hfill & {\rho }_g\left[\frac{u_{\varphi }^2}{r}+\frac{u_{\theta }^2}{r}-\frac{\partial \Phi }{\partial r}+\frac{1}{{\rho }_g}\left(f_r-\frac{\partial P}{\partial r}\right)+\frac{{\rho }_d}{{\rho }_g}\frac{{\Omega }_k}{S_t}\left({\upsilon }_r-u_r\right)\right],\hfill \\ \frac{\partial Y_{\varphi }}{\partial t}+\nabla \cdot \left(Y_{\varphi }u\right)\hfill & =\hfill & {\rho }_{g}r\left[\frac{u_{\theta }^2\text{\hspace{0.17em}cot}\left(\varphi \right)}{r}-\frac{1}{r{\rho }_g}\left(f_{\varphi }-\frac{1}{r}\frac{\partial P}{\partial \varphi }\right)+\frac{{\rho }_d}{{\rho }_g}\frac{{\Omega }_k}{S_t}\left({\upsilon }_{\varphi }-u_{\varphi }\right)\right],\hfill \\ \frac{\partial Y_{\theta }}{\partial t}+\nabla \cdot \left(Y_{\theta }u\right)\hfill & =\hfill & {\rho }_{g}r\text{\hspace{0.17em}sin}\left(\varphi \right)\left[\frac{1}{{\rho }_g}\left(f_{\theta }-\frac{1}{r\text{\hspace{0.17em}sin\hspace{0.17em}}\varphi }\frac{\partial P}{\partial \theta }\right)+F_{\theta }+\frac{{\rho }_d}{{\rho }_g}\frac{{\Omega }_k}{S_t}\left({\upsilon }_{\theta }-u_{\theta }\right)\right].\hfill \end{array}$$(1)

The function f = (f_r, fφ, f_θ) is the divergence of the stress tensor (see, for example, Tassoul 1978). The gravitational potential, Φ, depends on r and therefore contributes to a radial acceleration in the above equations. The term F_θ is the forcing term defined in Section 3. The fluid equations are closed by the definition of an equation of state (EoS). In this paper, we consider a locally isothermal disk with the pressure, $P=c_s^2{\rho }_g$, where the sound speed, c_s, is given by c_s = Ω_kH_g with the pressure scale height, H_g = h₀r, and h₀ is the disk aspect ratio.

The continuity equation for the dust is expressed as $\frac{\partial {\rho }_d}{\partial t}+\nabla \cdot \left({\rho }_{d}v\right)=\nabla \cdot j$, with $$j=D_d{\rho }_g\nabla \left(\frac{\rho d}{\rho g}\right),$$(2)

where D_d is the dust diffusion coefficient (Liu & Bai 2023): D_d = ρ_g/(ρ_g + ρ_d)D_g with $D_g=\alpha H_g^2{\text{Ω}}_k$ and α is the Shakura–Sunyaev turbulent viscosity parameter related to the disk’s viscosity, ν, by ν = αc_sH. This formula for the dust diffusion coefficient can be simply derived by assuming that for a turbulent fluid, the turbulent kinetic energy is $1/2{\rho }_g{u}_{turb}^2$. However, if we add ρ_d, the energy becomes $1/2\left({\rho }_g+{\rho }_d\right){u^{\prime }}_{turb}^2$, where ${u^{\prime }}_{turb}$ is the new turbulent velocity. For the energy to be the same (the only source of energy comes from the turbulence of the gas, the dust is just an inertia), we have $${u^{\prime }}_{turb}=\sqrt{{\rho }_g/\left({\rho }_g+{\rho }_d\right)}{u}_{turb}.$$

In other words, the turbulent velocity is quenched by a factor $\sqrt{{\rho }_g/\left({\rho }_g+{\rho }_d\right)}$. We notice that the diffusion coefficient is related to the turbulent velocity, u_turb, by the formula $D_g=u_{turb}^2/{\Omega }_k$ so that the dust diffusion coefficient is quenched by a factor ρg/(ρg + ρd).

The Navier-Stokes equations for the dust’s radial momentum, J_r = ρ_dv_r, the polar momentum, J_φ = ρ_drvφ, and the angular momentum, J_θ = ρ_dr sin(φ)V_θ, is expressed as $$\{\begin{array}{lll}\frac{\partial J_r}{\partial t}+\nabla \cdot \left(J_{r}v\right)\hfill & =\hfill & {\rho }_d\left[\frac{{\upsilon }_{\varphi }^2}{r}+\frac{{\upsilon }_{\theta }^2}{r}-\frac{\partial \Phi }{\partial r}-\frac{{\Omega }_K}{\text{St}}\left({\upsilon }_r-u_r\right)\right],\hfill \\ \frac{\partial J_{\varphi }}{\partial t}+\nabla \cdot \left(J_{\varphi }v\right)\hfill & =\hfill & {\rho }_{d}r\left[\frac{{\upsilon }_{\theta }^2\text{\hspace{0.17em}cot}\left(\varphi \right)}{r}-\frac{\partial \Phi }{\partial r}-\frac{{\Omega }_K}{\text{St}}\left({\upsilon }_r-u_r\right)\right],\hfill \\ \frac{\partial J_{\theta }}{\partial t}+\nabla \cdot \left(J_{\theta }v\right)\hfill & =\hfill & -{\rho }_{d}r\sin \left(\varphi \right)\left[\frac{{\Omega }_K}{\text{St}}\left({\upsilon }_{\theta }-u_{\theta }\right)\right].\hfill \end{array}$$(3)

3 Disk model

We consider disks with background density profile $${\Sigma }_b\left(r\right)={\text{Σ}}_0{\left(r/r_0\right)}^{-p},$$(4)

where Σ₀ is the gas surface density at r/r₀ = 1 and p is the power law slope, while r₀ is the unit of distance in astronomical units (au). In all of our simulations, we applied the following parameters: Σ₀ = 10⁻⁴, p = 1, and aspect ratio, h₀ = 0.05. The code units are defined such that G = M_* ≡ M_⊙ = 1. We adopt r₀ = 1 au throughout this work, although the results may be rescaled to any physical value of r₀. The unit of time is (r₀/au)^3/2/(2π) yr. We denote by P_orb the orbital period at r = r₀. For simplicity, in the following, we omit the r₀ term. In the 3D model, we define the background volume density from hydrostatic equilibrium in cylindrical coordinates (R, z) as $${\rho }_b\left(R,z\right)=\frac{{\text{Σ}}_0}{h_0\sqrt{2\pi }}{R}^{-p-1}\text{exp}-{\left(\frac{z}{2H_g}\right)}^2,$$(5)

where z = r cos(φ) and R = r sin(ϕ), with R ~ r in the thin disk approximation used in this paper.

3.1 Density bump: Initialization and forcing term

Following a similar approach to Liu & Bai (2023), we modeled the gas density bump as a Gaussian profile and maintained it by applying a torque via a forcing term, F_θ, in the azimuthal direction that balances viscous diffusion. We introduced the density bump in an otherwise monotonically decreasing radial gas density profile. Overall, we are primarily interested in the 3D model, however, we start our discussion with a simulation involving the 2D (r, θ) case. Therefore, we want to start by defining the density bump and the forcing term in the two cases. More precisely, in the case of a 2D vertically integrated disk, we set an initial surface density profile of the form, $${\Sigma }_g\left(r\right)={\text{Σ}}_b\text{\hspace{0.17em}exp}\left[A\text{\hspace{0.17em}exp}\left(-\frac{{\left(r-r_{bump}\right)}^2}{2\Delta w^2}\right)\right],$$(6)

where A is the bump amplitude and Δw is the bump width. The position of the bump maximum is r_bump. In 3D, the gas density profile takes the following form, with the same bump definition as in 2D (Eq. (6)), expressed as $${\rho }_g\left(r,\varphi \right)={\rho }_b\text{\hspace{0.17em}exp}\left[A\text{\hspace{0.17em}exp}\left(-\frac{{\left(r-r_{bump}\right)}^2}{2\Delta w^2}\right)\right].$$(7)

To maintain the density bump against viscous diffusion, we have introduced a forcing term in the azimuthal component of Eq. (1). Since the pressure has no azimuthal dependence, to compensate for the viscous diffusion the forcing term is $$F_{\theta }=-\frac{1}{d}{f}_{\theta },$$(8)

with the function d corresponding to the surface density, Σ_g, or the volume density, ρ_g, according to the dimension of the problem and f_θ is the azimuthal component of the divergence of the stress tensor. In the 3D case f_θ is given by $$f_{\theta }=\frac{1}{r\text{\hspace{0.17em}sin\hspace{0.17em}}\varphi }\left[\frac{\text{sin\hspace{0.17em}}\varphi }{r}\frac{\partial }{\partial r}\left(r^2{\tau }_{{\theta }_r}\right)+\frac{\partial }{\partial \varphi }\left({\tau }_{\varphi \theta }\text{\hspace{0.17em}sin\hspace{0.17em}}\varphi \right)+\frac{\partial {\tau }_{\theta \theta }}{\partial \theta }\right]+\frac{{\tau }_{r\theta }}{r}+\frac{{\tau }_{\theta \varphi }\text{\hspace{0.17em}cot\hspace{0.17em}}\varphi }{r},$$(9)

with the components of the viscous stress tensor from Tassoul (1978) simplified with the condition u_r = 0 and u_ϕ = 0, $$\{\begin{array}{lllll}{\tau }_{\theta r}\hfill & =\hfill & {\tau }_{r\theta }\hfill & =\hfill & {\rho }_{g}vr\frac{\partial }{\partial r}\left(\frac{u_{\theta }}{r}\right),\hfill \\ {\tau }_{\varphi \theta }\hfill & =\hfill & {\tau }_{\theta \varphi }\hfill & =\hfill & {\rho }_{g}v\frac{\text{sin\hspace{0.17em}}\varphi }{r}\frac{\partial }{\partial \varphi }\left(\frac{u_{\theta }}{\text{sin\hspace{0.17em}}\varphi }\right),\hfill \\ {\tau }_{\theta \theta }\hfill & \hfill & \hfill & =\hfill & 2{\rho }_{g}v\left(\frac{1}{r\text{\hspace{0.17em}sin\hspace{0.17em}}\varphi }\frac{\partial u_{\theta }}{\partial \theta }\right).\hfill \end{array}$$(10)

Considering that u_θ has no azimuthal dependence, we have τ_θθ = 0 and the forcing term is finally expressed as $$F_{\theta }=-\frac{1}{{\rho }_g}\left[\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2{\tau }_{\theta r}\right)+\frac{{\tau }_{r_{\theta }}}{r}+\frac{1}{r\sin \varphi }\frac{\partial }{\partial \varphi }\left({\tau }_{\varphi \theta }\text{\hspace{0.17em}sin\hspace{0.17em}}\varphi \right)+\frac{{\tau }_{\theta \varphi }\text{\hspace{0.17em}cot\hspace{0.17em}}\varphi }{r}\right].$$(11)

By neglecting the terms depending on the colatitude in Eq. (11), we have the forcing term in the 2D case, $$F_{\theta }=-\frac{1}{{\text{Σ}}_g}\left(\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2{\tau }_{\theta r}\right)+\frac{{\tau }_{r\theta }}{r}\right),$$(12)

where ${\tau }_{\theta r}={\tau }_{r\theta }={\Sigma }_{g}vr\frac{\partial }{\partial r}\left(u_{\theta }/r\right)$.

We built an equilibrium disk by considering the initial radial velocity, u_r = 0, and the azimuthal velocity from the centrifugal balance, $$u_{\theta }^2/r={\text{Ω}}_k^{2}r+\frac{1}{d}\frac{\partial P}{\partial r},$$(13)

and $P=c_s^{2}d$. In the 3D case, we set the initial polar velocity to zero.

3.2 Initial conditions for dust

We recall that the gas rotates at a slightly sub-Keplerian speed and the relative difference with respect to the Keplerian speed, V_k, is indicated by η = (V_k − u_θ)/V_k, which can be computed as $$\eta =-\frac{1}{2}{\left(\frac{H}{r}\right)}^2\frac{\partial \text{\hspace{0.17em}log\hspace{0.17em}}P}{\partial \text{\hspace{0.17em}log\hspace{0.17em}}r},$$(14)

with $\eta \propto h_0^2$ in usual power law disks. The radial dust velocity component, υ_r, is then initialized (Takeuchi & Lin 2002) as $${\text{v}}_r=\frac{u_r}{S_t^2+1}-\frac{2S_t}{S_t^2+1}\eta {\text{V}}_K.$$(15)

We remark that the gas velocity u_r is expected to be small. Neglecting u_r if the pressure gradient is negative then η > 0 and the dust moves toward the star, while in the case of a positive pressure gradient (η < 0) the situation is reversed. In the presence of a pressure maximum (η = 0), the dust can be trapped. The initial azimuthal velocity is $${\text{v}}_{\theta }=u_{\theta }-\frac{1}{2}{S}_t{\text{v}}_r,$$(16)

while the initial polar velocity is zero as for the gas. The dust density is initialized as a fraction, ϵ₀, typically of a few percent of the gas density.

3.3 Boundary conditions

We provide the values of dust and gas components in the ghost cells (for the radial and vertical directions) according to values of the closest neighbor or first active cell. We considered a global disk with periodic azimuthal boundary conditions. For the gas components in the radial direction, we used the classical prescription (de Val-Borro et al. 2006) of the evanescent boundary condition.

For the dust component, we implemented an open inner radial boundary, allowing dust to drift out of the computational domain. Specifically, dust with negative radial velocity in the first active cells is advected into the ghost cells. To prevent a dust pile-up, we copied the density from the first active cells into the ghost cells. The azimuthal velocity component of the first active cells was extrapolated using a Keplerian profile within the ghost cells. At the outer radial boundary, the azimuthal velocity of the last active cells is similarly extrapolated using a Keplerian profile in the ghost cells. For the dust density and radial velocity, we considered two distinct boundary conditions, detailed below.

• No inflow : the outer part of the simulation domain gets depleted according to radial dust drift towards the pressure bump. This setup is intended to mimic the case where a barrier in the disk beyond the simulation domain halts dust flow toward the considered pressure bump so that only the local dust can concentrate there. Here, the radial velocity in the ghost cells is set to zero.

• Inflow: we impose a constant dust inflow by setting the ghost cell’s radial velocity to the theoretical value (Eq. (15)) and resetting its density to the initial profile. This simulates a continuous replenishment of the dust reservoir in the simulation domain via advection and diffusion. In 3D simulations, we additionally account for dust settling by renormalizing the vertical density distribution in the ghost cells at each time step, using the vertical profile of the first active cell.

At the vertical boundaries, we enforced a no-flow condition by setting the vertical velocity component to zero for both gas and dust, preventing any inflow or outflow from the domain. All other quantities in the ghost cells are directly copied from the nearest active cells.

4 2D simulation results

The first goal of this paper is to reproduce the 2D local shearing-box results of Liu & Bai (2023), hereafter LB23, using 2D global disk simulations. We begin by describing the initial disk setup and then present our findings for type II DRWI, highlighting similarities and differences with LB23.

We begin with an equilibrium disk with a homogeneous global dust-to-gas density ratio, ϵ₀ = 0.02, and viscosity, α = 10⁻⁴, and initialized a gas density bump centered at r_bump ≡ r₀ = 1, with an amplitude of A = 0.8 and width of Δw/h₀ = 1.5 (the same density bump parameters and viscosity used in LB23’s fiducial type II run), which is expected to be stable to the classic RWI (Chang et al. 2023). The radial simulation domain spans the interval [0.6, 1.4] and we consider (N_r × Ν_φ) = (1024 × 3072) grid cells.

The established pressure bump reverses the sign of η (Eq. (14)) and traps the inward drifting dust and forms the dust ring, locally increasing the dust-to-gas ratio within the ring. The initial condition of our disk differs from that of LB23, as they begin from a preestablished dust ring in equilibrium which has sufficient dust-to-gas ratio to trigger type II DRWI, once perturbed. We instead impose dust inflow from the outer radial boundary (according to Section 3.3) to gradually increase the dust-to-gas ratio, ensuring that at the pressure bump, we can eventually reach the critical dust-to-gas ratio required for the dust ring to become unstable. Following LB23, we characterize the dust enrichment at the pressure bump via $$f_g{}_{\text{min}}=\text{min}\left(\frac{{\overline{\text{Σ}}}_g}{{\overline{\text{Σ}}}_g+{\overline{\text{Σ}}}_d}\right),$$(17)

where ${\overline{\Sigma }}_g$ and ${\overline{\Sigma }}_d$ are azimuthally averaged quantities and the minimum is taken over the simulation domain and occurs near the pressure maximum. We find that the type II DRWI is triggered once f_gmin decreases below 0.5 (see Fig. 2). LB23 found a threshold at f_gmin = 0.536. Given the differences in the simulation set-up, the results are in satisfactory agreement.

Figure 1 shows the evolution of the dust ring after the onset of the type II DRWI. After the ring becomes unstable, dust clumps eventually form, which can be seen in the panels at 3000 and 3350 P_orb. The panel at 3350 P_orb is repeated with the dust clumps circled in black. We identified the dust clumps by applying an algorithm to find the local dust density maxima; then, the maxima exceeding three times the azimuthal average are considered to be dust clumps. Although LB23 note the transient nature of the dust clumps in their simulations, we extended their analysis by examining the behavior and evolution of these clumps when dust is continuously supplied to the ring. Figure 2 presents the evolution of the number of dust clumps in the ring (top panel), the maximum dust clump strength (middle panel), and the evolution of f_gmin. We defined the clump strength as $\text{max}\left[{\Sigma }_d/{\overline{\Sigma }}_d\left(r_{\text{clump}}\right)\right]$, with ${\overline{\Sigma }}_d\left(r_{\text{clump}}\right)$ being the azimuthal average of the dust density at the clump’s location. We find that both the number of dust clumps and the maximum density of the clumps increase as more dust is collected in the ring (i.e., as f_gmin ^decreases).

The instability reaches saturation near 4500 P_orb, after which the number and strength of the dust clumps continue to fluctuate, but show no significant increase. The number of dust clumps averages at 20–25 (Fig. 2 top panel) and the dust clump strength around 10 (Fig. 2 middle panel), corresponding to Σ_d/Σ_b ≈ 30, where Σ_b is defined in Eq. (4). The maximum dust clump strength we find is close to 18, which corresponds to Σ_d/Σ_b ≈ 34.

Although our 2D simulations produce dust clumps, assessing whether they reach the Hill density, the threshold for gravitational collapse, is difficult for two reasons. First, a fully 3D treatment is required to measure the volume density of dust on the midplane (although some estimates can be used). Second, because our code does not include self-gravity, all resulting densities scale with the value of Σ_b we assume. Nevertheless, it is important at this stage to extend the analysis to 3D simulations, as different behaviors and instabilities may appear.

Fig. 1 Snapshots in the 2D type II DRWI run of the dust density in the azimuthal-radial plane. The time is annotated on the top left of each panel. The panels are colored in a power-law scale and the dust surface density is plotted with respect to the initial background gas density profile (Eq. (4)). The panel at 3350 Porb is repeated to the right of the dotted black line, with the identified dust clumps circled in black.

Fig. 2 Evolution of the dust clumps in the 2D type II DRWI run. Top panel: number of identified dust clumps at each time. Middle panel: maximum density of the dust clumps at each time, normalized by the average dust density at the clump’s radial location. Bottom panel: evolution of the minimum gas-to-dust ratio.

5 3D simulation results: Viscous disk

Here, we expand on the work by LB23 using 3D global disk simulations. The vertical domain spans the φ interval [1.54,1.60] (φ = π/2 being the midplane), which is ±0.6h₀ about the mid-plane (much larger than the scale height of the dust layer), and the radial domain extends over [0.8, 1.2]. The full 2π range in θ is considered. Additional parameters for all 3D simulations performed in this work are provided in Table 1. To begin, the setup is identical to the 2D case, but with an initial homogeneous global dust distribution of ϵ₀ = 0.05. Since 3D simulations are computationally expensive, we have used this unusually large ϵ₀ value to reach the critical fg_min for DRWI type II in fewer orbital periods than in 2D. The dust is then allowed to settle self-consistently and drift radially under the hydrodynamics of the code. Again, we impose dust inflow from the outer radial boundary, as described in Section 3.3. The gas density bump in 3D is established according to Eq. (7) and enforced according to Eq. (11).

To thoroughly study DRWI and the evolution of the dust ring in 3D, we performed a series of tests, first examining the case of a viscous disk for a comparison with the previous 2D results and then considering a disk with vanishing viscosity.

5.1 Stratified viscous disk

Our first test considered a 3D vertically stratified viscous disk. Figure 3 shows the dust density after 1790 P_orb. No instability or dust clumping develops, even though f_gmin decreases well below the critical value required to trigger type II in 2D, and we find the same result in tests with a larger Stokes number St = 0.1. This suppression is likely a consequence of the full 3D structure: dust settling concentrates dust into a thin midplane layer, leaving the upper regions dust-poor and unable to develop or sustain the instability. Viscous coupling between vertical layers then allows the unperturbed gas above and below the midplane to damp and smooth out any perturbations that would otherwise develop in the dust-rich midplane, preventing the growth of the instability throughout the column. We tested this interpretation in the following section using a 3D vertically unstratified test.

5.2 Unstratified viscous disk

To simulate an unstratified disk, we removed the vertical component of gravity on both dust and gas and, thus, the gas is initialized with a homogeneous vertical density profile and the dust-to-gas ratio is the same at every z, so that each layer is identical. Dust settling does not occur because of the suppression of the vertical component of the gravitational force.

Figure 4 shows that the dust ring in the unstratified case exhibits stronger perturbations than in the stratified run, even developing dust clumps and closely resembling the behavior seen in the 2D configuration (Figure 1). Because dust is present throughout the vertical column, all layers become unstable, allowing perturbations to grow and dust clumps to form. This confirms our interpretation that, in the stratified disk, as dust settling concentrates solids in the midplane, it leaves the upper layers dust-poor and stable; viscous coupling between these layers then damps perturbations arising at the midplane, ultimately suppressing the instability.

Fig. 3 Same as Fig. 1, but depicting the volume density of dust on the midplane in a 3D simulation with the same pressure bump parameters and gas viscosity as in the 2D model of Section 4.

6 3D simulation results: Inviscid disk

Given the type II DRWI is suppressed in 3D in a disk with α = 10⁻⁴, we wanted to test what happens in the limit of an inviscid disk. For this purpose, we set α = 0. Of course, as all numerical codes, the simulation is still affected by some numerical viscosity, but various tests suggest that it is weaker than α = 10⁻⁵ (Lega et al. 2021). Importantly, in the inviscid case, no diffusion is prescribed on the dust evolution. The only diffusion that is eventually present is that induced self-consistently by the gas dynamics. Furthermore, the function f_θ (Eq. (9)) is now identically null.

To investigate the role of dust supply, we considered two inviscid simulations that differ in their dust reservoir. In the first case, the disk contains a finite amount of dust with no inflow from the outer boundary, representing a situation in which dust has already drifted inward and accumulated in the pressure bump. Accordingly, we adopt a relatively high initial dust-to-gas ratio of ϵ₀ = 0.05. In the second case, we allow for continuous dust inflow from the outer disk, starting from a standard interstellar value of ϵ₀ = 0.01.

Fig. 4 Same as Fig. 3, but for the 3D unstratified simulation.

6.1 No dust inflow

We began with a setup containing a finite amount of dust (no inflow). In the absence of turbulent diffusion driven by α-viscosity, dust is expected to settle into a thin, high-density midplane layer and to form correspondingly a narrower, sharper ring than in the viscous case. In this scenario, dust settling and radial contraction alone may be sufficient to reach the critical dust-to-gas ratio at the midplane required for instability within the ring.

Figure 5 shows the evolution of the dust ring. Without viscous diffusion, multiple dust subrings emerge and interact dynamically. These subrings gradually merge, decreasing in number while concentrating the dust into a few dominant structures. Ultimately, two long-lived rings remain, spanning a finite radial width of w_d = 0.02. Moreover, neither of these rings become razor-thin.

Figure 6 illustrates the relationship between the dust-density peaks (brown) marking the locations of the rings and the radial profile of η (pink). The solid pink curve gives the pressure gradient, while the dashed pink curve provides the true value of η, defined as η = (υ_kep − υ_θ)/υ_θ. The latter is simply the former quenched by the factor ρ_d/ρ_g, as expected and, thus, they have the same zero points. The zeros of η (indicated by the horizontal dashed grey line) identify sign reversals and those with positive slope coincide with the positions of the dust rings. In this high dust-to-gas ratio regime, the strengthened dust back-reaction, amplified by radial dust pile-ups and the radial dependence of the drift velocity, modifies the gas density profile and, thus, the pressure structure. Although the gas is only weakly compressible, this is sufficient, particularly when η is small, to generate new zeros in the η profile, where additional dust rings form. This happens when ρ_d/ρ_g approaches 10. We discuss this in more detail in Section 6.3.

We also checked that the inversions in the sign of η occur also in a simulation with an adiabatic equation of state (instead of isothermal), with a cooling timescale of 100 orbital periods, consistent with the disk at 1 au. This is because the little compression needed to reverse the sign of η produces negligible heating in the gas. Finally, we performed additional simulations with intermediate values of α = 10⁻⁵ and 10⁻⁶. We find that dust subrings associated with η sign inversions develop at α = 10⁻⁶, whereas they are suppressed at α = 10⁻⁵, with the evolution resembling that of the viscous case. This suggests that the transition occurs for α between 10⁻⁶ and 10⁻⁵.

6.2 No dust inflow: 2D r – z validation case

Before proceeding with the next set of simulations, we performed a validation test to verify we would recover the same results (i.e., formation of multiple axis-symmetric rings) in a 2D r – z simulation. 2D r – z simulations are typical of most studies of the streaming instability (e.g., Li & Youdin 2021) and offer the great advantage of enhancing resolution within the same computation time with respect to a full 3D simulation.

Figure 7 shows the evolution of the dust surface density in a 2D r – z disk (left) and in the fully 3D disk (right). Both simulations use the same set of parameters and resolution in the radial and vertical directions for a proper comparison, even though the 2D simulations can be run at much higher resolution. The results are qualitatively very similar, demonstrating that the 2D r – z configuration captures the essential behavior observed in 3D. Therefore, we adopted the 2D setup for the next simulation, described in Section 6.3, so that we could push the resolution as high as possible.

6.3 Dust inflow

Next, we examined the inviscid evolution under continuous dust inflow, starting from a more realistic initial dust-to-gas ratio of ϵ₀ = 0.01. We adopted a 2D r – z axisymmetric disk, implemented using a single grid cell in θ, and we used the same parameters as in the previous inviscid run (Section 6.1). The resolution is increased by a factor of 10 in both the radial and vertical directions. Figure 8 shows the resulting evolution with time versus the radial distance to the star of, respectively: the vertically integrated density ratio, Σ_d/Σ_g (second panel), η computed from the pressure gradient at the midplane (third panel), and the radial velocity of the dust at the midplane υ_r (bottom panel). In the top panel of Fig. 8, we show the vertical distribution of the dust at the end of the simulation, at t = 1500P_orb. In the second panel, we observe filaments of high density which correspond to the white regions along z = 0 in the top panel.

6.3.1 Multiple ring formation and streaming instability

Although the simulation starts from a lower dust-to-gas ratio, dust progressively drifts toward and accumulates at the pressure bump and the resulting back-reaction becomes strong enough to reshape the gas profile. As in the no inflow case, this generates multiple new locations where η = 0 (Fig. 8 third panel). The continuous supply of dust steadily modifies the gas density structure, partially flattening and shifting the original pressure bump. The partial flattening can be deduced from the second panel, which shows that the radial distribution of dust first contracts (in the first ~ 300 orbital periods) and then expands again. The shifting of the main pressure bump is shown in the third panel where the main boundary between the blue and red regions (η < 0 and >0 respectively) moves from ~ 0.98 to ~ 0.965 in the course of the simulation. This occurs because dust drift modifies the radial gas velocity through back-reaction, thereby reshaping the gas density and pressure profile. Remember that the azimuthal stress term f_θ (Eq. (9)) vanishes because ν = 0.

As the pressure structure evolves, dust is able to drift inward of the original bump, transforming the ring into a nonuniform multi-ring configuration with localized regions of enhanced dust-to-gas ratio. New filaments continuously form and sometimes merge, each associated with new η = 0 points. The dust radial drift velocity is essentially zero in the filaments (see the predominant white color in the bottom panel of Fig. 8), meaning that dust is trapped in the filament and co-moving with it. However, because of turbulent fluctuations of the velocity (the red and blue dots in the vicinity of the main white filaments), some dust can diffuse out of the filament and then starts to drift again until it is captured in a new filament.

The first two panels of Fig. 8 show a pattern reminiscent of the streaming instability, which is expected as dust approaches the pressure maximum (initially located at r = 0.983), slows down and accumulates on its outer side, progressively enhancing the local dust-to-gas ratio until the threshold for the streaming instability can be reached. The difference from the classic streaming instability is that as soon as pronounced filaments appear, η changes sign, as indicated by the appearance of blue stripes in the red domain of the third panel of Fig. 8. The theory of the streaming instability is usually developed for an incompressible fluid, but in reality the gas is weakly compressed by the dust filaments, inducing oscillations of η around its well defined positive value. In the vicinity of the pressure bump, η is small, and the weak compression of the gas is enough to force it to change sign. Once this happens, the dust is trapped in each filament (see bottom panel, where the dust has zero radial velocity within the filaments).

Finally, in the top panel of Fig. 8 we observe that the dust layer becomes increasingly vertically extended at radii r > 1. In an inviscid disk, we would expect the dust to settle into an extremely thin midplane layer, yet this is not what was produced by the simulation. Although vertical stirring could, in principle, arise from the Kelvin-Helmholtz instability (KHI), the Richardson-number criterion is not satisfied due to the small value of η. To identify the mechanism responsible for this vertical “puffing up,” we examined the dust distribution and velocities at t = 1000 P_orb, as detailed in the following subsection. We also include a comparison between the effective α derived from the vertical dust structure and that inferred from the radial dust distribution.

The resulting radial broadening of the dust distribution prevents ρ_d /ρ_g from increasing indefinitely, implying that additional processes may still be required to raise the dust-to-gas ratio to densities (Hill density) high enough for gravitational collapse. One such mechanism is gas evaporation, which occurs later in the evolution of the disk. We explore this phase of the disk’s evolution and its implications on planetesimal formation in Section 7.

Fig. 5 Same as Fig. 3, but for the 3D inviscid (α = 0) simulation without dust inflow.

Fig. 6 Radial profiles for the 3D inviscid simulation without dust inflow of the pressure (green), dust density (brown), and η (pink). The solid pink curve is the pressure gradient, while the dashed pink curve is the true value of η.

Fig. 7 Comparison between the dust ring evolution in the 2D r – z inviscid disk versus the full 3D disk with the same radial and vertical resolution.

6.3.2 Vertical versus radial dust profiles

The vertical dust distribution and gas dynamics at t = 1000 P_orb are shown in Figure 9. As described above, we observe that the vertical dust profile varies with radius, becoming increasingly puffed up outside the main pressure bump. Within the dust rings, the dust settles into a thinner midplane layer than in the surrounding regions, yet still maintains a finite vertical width. We see alternating patterns of positive and negative vertical gas motions (Fig. 9 third column), which are reminiscent of the vertical shear instability (VSI). However, we tested that VSI is not responsible for the vertical stirring of the dust. Indeed, in a simulation without dust, the magnitude of the velocities, representing the strength of the VSI, was three orders of magnitude lower than in the simulation with dust.

In addition, an adiabatic simulation with slow cooling produces equivalent results. We think the strong vertical motion of the gas in the simulation with dust arises because the dust back-reaction, which becomes more significant when there is a high dust-to-gas ratio, tends to compress the gas in the radial direction. Indeed, we observed locations of radial compression in the radial gas velocity profile at the midplane. These compression points can be seen in the second column of Figure 9, where the velocity profile transitions from red (outward flowing gas) to blue (inward flowing gas). Because the gas is only mildly compressible by the dust (indeed, the theory of the streaming instability is traditionally developed in the framework of an incompressible fluid), the radial compression of the gas induces an upward vertical flow away from the midplane, as in the streaming instability model (Squire & Hopkins 2020), entraining and lifting the dust with it. The resulting pattern (Fig. 9, third column) consists of alternating red–blue vertical bands, while the spikes in dust density are aligned with the blue (upward flow) regions. The final panel, which overlays dust density on the vertical gas velocity, clearly demonstrates this correlation: dust spikes are co-located with zones of upward-moving gas, coinciding with the locations of strongest radial gas compression (e.g., Fig. 9 bottom right panel at r = 0.978, 0.99).

From the vertical dust distribution, we can compute an effective α that induces a sufficient vertical diffusion and compare this to the turbulent α computed from the Reynold’s stress tensor, T^Rey. This is usually done by fitting a Gaussian to the vertical dust density profile, from which the vertical height of the dust layer H_d can be extracted and then be used to compute α from $H_d=H_g\sqrt{\alpha /\text{St}}\sqrt{{\rho }_g/\left({\rho }_g+{\rho }_d\right)}$. Outside of the dust ring near r = 1.04 we find H_d ≈ 10⁻³ and ρ_g/(ρ_g + ρ_d) ≈ 0.98, from which we find α ≈ 10⁻⁵ .Within the dust ring between r = 0.975–0.985, we find an average H_d on the order of ~ 10⁻⁴ and an average α on the order of ~10⁻⁶.

We can then compute α from the rϕ component of Reynold’s stress tensor, $T_{r\phi }^{Rey}$, using $\alpha =T_{r\phi }^{Rey}/P$, where P is the gas pressure. The value of α is then vertically averaged over the dust-containing region. Near r = 1.04, we find α ~ 10⁻⁵ and within the ring, we find an average α on the order of ~10⁻⁶. Both results are in agreement with the effective α values computed from the vertical dust height. As such, we can conclude that it is the gas turbulence induced by the dust-back reaction (as it is orders of magnitude smaller if the dust is removed) that limits the settling of the dust.

Finally, we also computed α from the vertical component of the dispersion velocity using the formula ${\upsilon }_z^{rms}=\sqrt{\alpha }{c}_s$ where ${\upsilon }_z^{rms}$ is the root mean square of the vertical velocity and c_s is the sound speed. Outside the ring we find α ≈ 10⁻⁴, whereas inside the ring we find α ≈ 10⁻⁵. Interestingly, both values are an order of magnitude larger than the corresponding values computed above.

We turn our attention to the radial width of the ring. The top panel in Fig. 8 shows that the ring is composed of multiple filaments, each at location of η = 0. An observer with limited resolution would not be able to detect the individual filaments and would instead interpret this as a single ring with a broad width of 0.025 in r/r₀. Interpreting this as a result of turbulent diffusion and using the formula $w_{ring}=\Delta w\sqrt{{\alpha }_{ring}\text{/St}}$, where Δw is the width of the gas profile, we would then conclude that α_ring ~ 3 × 10⁻³. Performing the same exercise for the inviscid case without dust inflow (Section 6.1), and considering the final panel in Figure 5, we see that the broad width of the dust ring would appear as 0.02 in r/r₀, resulting in an α_ring ~ 2 × 10⁻³.

These values are close to those found in Dullemond et al. (2018) for several rings observed in protoplanetary disks. This large value of α, compared to that obtained from the vertical distribution of dust in protoplanetary disks (Villenave et al. 2022), have prompted the discussion that turbulence might be anisotropic (i.e., stronger in the radial direction than in the vertical direction). Our results suggest that the radial width of rings may not be due to turbulence but to the existence of unresolved filaments separated by a few 10⁻² in relative units. Thus, turbulence in disks does not need to be anisotropic.

Fig. 8 Evolution of the high resolution 2D r – z inviscid disk with ϵ0 = 0.01 and dust inflow. In the top panel, we plot the r – z view of the dust to gas volume density at the final integration time t = 1500 Porb. The lower panels show the evolution with time versus r/r0 of, respectively: the vertically integrated density ratio, Σd/Σg (second panel), the η values computed from midplane pressure gradient (third panel) and the midplane dust’s radial velocity υr (bottom panel). In the third panel, a transition from blue to red indicates an η sign inversion and thus a dust trapping location.

Fig. 9 Vertical dust structure as a function of radius for the high-resolution 2D r – z inviscid disk simulation with ϵ0 = 0.01 and dust inflow at t = 1000 Porb. The top panel shows a broader radial region, while the bottom panel shows a zoomed-in view of the dust ring region. In the top panel, the outer edge of the dust ring is indicated with a dashed red line. Each set of panels (from left to right) shows: the dust density, radial gas velocity, vertical gas velocity (both normalized by the azimuthal Keplerian velocity), and dust density overlaid with the vertical gas velocity.

7 Late planetesimal formation: Gas evaporation

In addressing planetesimal formation, we recall that our code does not include self-gravity with respect to either the gas or the dust. Thus, even if we were to reach the Hill density, we would not observe the formation of a self-gravitating clump of dust. Moreover, it has been shown that the real dust dynamics diverge from that of a simulation neglecting self-gravity when the dust density reaches a fraction of a few of ρ_Hill (Johansen et al. 2007). A simulation without self-gravity (as done here) only gives information about the evolution of the dust-to-gas ratio and not about the absolute values of either of these quantities. This is why in Figures 3, 4, and 5, we plot ρ_d/ρ_b, where ρ_b is the density that the gas would have in the absence of the pressure bump (Eqs. (4) and (5)). With a typical gas density at 1 AU of 10⁻⁹ g/cm³, the maximal dust density in Figure 5 would reach 4 × 10⁻⁷ g/cm³, which is the Hill density. However, it is possible that at a later time in the evolutionary sequence of the disk, the gas was already partially depleted (i.e., ρ_b ≪ 10⁻⁹ g/cm³) or the total amount of dust available (relative to the gas) was smaller than that considered in the simulation. In both these cases, the Hill density would still be out of reach. Thus, in this work, we consider the late evolutionary phase in which gas is gradually removed through photoevaporation. Under the assumption that at late times all of the available dust has been trapped in multiple rings, gas removal might become the dominant mechanism for enhancing the local dust-to-gas ratio and dust volume density within each ring, offering a possible explanation for the delayed formation of chondrites. We expect gas removal to increase the volume density of the dust by promoting both further sedimentation of dust toward the midplane and radial concentration, for the reasons explained below. In this study, we also investigated the possibility of triggering instabilities within the ring, which could further promote dust clumping. This was explored in both a viscous and inviscid disk, where we removed gas via a uniform exponential decay with a characteristic timescale of 200 orbital periods.

7.1 Viscous disk

We begin by considering a viscous disk, for which we first present a derivation of the expected scaling of the dust density within the ring, ρ_d,ring, with respect to the dust-to-gas mass ratio, M_d/M_g. In first considering the vertical profile, the dust scale height is given by $$H_d=H_g\sqrt{\frac{\alpha }{\text{St}}\frac{{\rho }_g}{{\rho }_g+{\rho }_d}.}$$(18)

We can also approximate ρ_g/(ρ_g + ρ_d) ≈ ρ_g/ρ_d, which is valid for ρ_d ≫ ρ_g, which is the correct regime once enough gas is removed. Then, we can re-write Eq. (18) as $$H_d=H_g\sqrt{\frac{\alpha }{\text{St}}\frac{{\text{Σ}}_g}{{\text{Σ}}_d}\frac{H_d}{H_g},}$$(19)

and solving for H_d, we find $$H_d=H_g\frac{\alpha }{\text{St}}\frac{{\text{Σ}}_g}{{\text{Σ}}_d}.$$(20)

Next, we compute the dust ring density at the midplane, $${\rho }_{d,\text{ring,0}}=\frac{{\text{Σ}}_d}{\sqrt{2\pi }{H}_d}=\frac{{\text{Σ}}_d}{\sqrt{2\pi \left(\alpha /\text{St}\right)}{H}_g}\left[\frac{1}{\sqrt{\alpha /\text{St}}}\frac{{\text{Σ}}_d}{{\text{Σ}}_g}\right].$$(21)

The first fraction term is the midplane density in the regime of low dust-to-gas ratio and the second term inside the brackets is the enhancement factor, which is proportional to Σ_d/Σ_g. The formulae in Eqs. (20) and (21) show that H_d and ρ_d depend on Σ_d/Σ_g, which increases as gas is removed. The reason is that the effective α is α[ρ_g/(ρ_d + ρ_g)]. If we had not considered the quenching of alpha by the dust-to-gas ratio, H_d and ρ_d would not have appeared to change during gas removal, which is unphysical (i.e., in the absence of gas, there is no reason for the dust to have a scale height). We stress that in our simulations, we kept St constant during gas removal. In reality, the dust velocity dispersal also scales as $\sqrt{\alpha \left[{\rho }_g/\left({\rho }_d+{\rho }_g\right)\right]\text{St}}$ c_s, so that in a fragmentation-limited regime, St should increase as ρ_g/(ρ_d + ρ_g) decreases. This effectively enhances dust sedimentation.

Considering the radial profile, the equation for Σ_d at the center of the dust ring is given by $2\pi r{\Sigma }_d=M_d/\left(\sqrt{2\pi }{w}_d\right)$, where M_d is the mass of the dust in the ring and w_d is its Gaussian width. The same is true for the gas, given by $2\pi r{\Sigma }_g=M_g/\left(\sqrt{2\pi }{w}_g\right)$. At the center of the ring, we then have $$\frac{{\text{Σ}}_d}{{\text{Σ}}_g}=\frac{M_d}{M_g}\frac{w_g}{w_d}.$$(22)

Plugging this into Eq. (21), we obtain the dependence of the dust density on the mass ratio. The radial contraction of the dust ring width is also expected to depend on M_d/M_g, following the same approach presented for the vertical contraction; however, there are two important differences between radial contraction and midplane settling. The first is that the locations at which η = 0 are not in steady state and shift radially over time, so the dust continuously readjusts to follow these moving trapping points; however, the midplane obviously does not move. The second is that, while the settling timescale is proportional to 1/(ΩSt), the radial contraction timescale is proportional to δr/(Stηυ_k), where δr is the distance from the pressure bump. The factor δr/η ≫ 1, so the contraction timescale is much longer than the sedimentation one. The gas may be removed faster than the dust has the time to radially adjust. In summary, although some degree of radial contraction is still expected, this time-dependent structure prevents a simple analytical description of the ring width. As a result, while the midplane dust density in the ring is expected to scale as a power of (M_d/M_g), the exact exponent cannot be determined without the expected scaling for the radial width evolution.

To investigate this behavior, we initiated the gas-removal phase from an intermediate stage of our previous 3D stratified viscous disk (Section 5.1) and by this time, a well-developed dust ring has already formed. We also repeat the same gas evaporation simulation in a 2D r – z disk at an identical resolution, which provides a controlled comparison to isolate the role of the azimuthal dimension and to highlight additional effects that may arise in the full 3D evolution, such as clumping.

Figure 10 shows the evolution of the maximum midplane dust density within the ring, plotted as a function of the total dust-to-gas mass ratio in the simulation domain for both the 2D (orange) and 3D (blue) cases. The evolution can be divided into three distinct stages, each characterized by a different scaling. The 2D and 3D evolution closely follow each other during the first two stages and diverge only in the final stage. In the first stage, with a slope of m = 1.02 (dashed line), gas removal enhances dust sedimentation towards the midplane (Eq. (20)), while the radial width of the ring also contracts. Radial contraction contributes a fraction of approximately 0.4 of the density increase, while vertical settling accounts for the remaining 0.6. The latter is less than expected from Eq. (21), but this is because we are already close to vertical resolution. Indeed, in the second stage, the thickness of the dust layer becomes limited by the vertical resolution, preventing further settling, so density growth is driven only by radial contraction, resulting in a weaker scaling of m = 0. 25 (dotted line). During this phase, subrings also begin to form within the main dust ring.

In the final stage, the 2D and 3D cases clearly diverge. The dust ring develops additional substructures, forming two distinct subrings, and radial contraction effectively ceases. This is clear in the 2D case, where we see that the dust density saturates. In contrast, the 3D simulation shows renewed growth, with a steeper scaling of m = 0.35 (dash-dotted line) This enhancement is driven by the development of azimuthal asymmetries within the rings, indicating the formation of dust clumps, a process that is absent in the 2D case. Once sufficient gas has been removed, the dust forms clumps that continue to grow in density, such that we expect the Hill density to eventually be reached and gravitational collapse to occur.

Fig. 10 Maximum dust density in the ring as a function of the dust-togas mass ratio during gas removal for both a 2D (orange) and 3D (blue) viscous disk (α = 10−4). The upper x-axis indicates orbital periods at the reference radius of r = 1. Linear best-fit lines are shown for three stages of the gas evaporation process, with their slopes indicated next to each line. The slopes of the first two stages are consistent between the 2D and 3D case and then diverge in the final stage, where there is clumping in 3D.

Fig. 11 Same as Fig. 5, but showing the evolution of the dust density as gas is removed from the disk. The dust density is plotted with respect to the Hill density and the total dust-to-gas mass ratio at each given time is indicated at the top of the panel. The orbital period in the top-left corner is given relative to the start of the original inviscid run with no inflow in Section 6.1.

7.2 Inviscid disk

Next, we sought to investigate gas evaporation in an inviscid disk. In this case, we initiated the gas-removal phase from the previous inviscid run without dust inflow (Section 6.1) at t = 1240P_orb, at which point the disk already exhibits well-developed dust ring substructures, but no clumps. Our analysis of the subsequent evolution of the dust density focuses on the dominant subring, namely, the most massive of the two principal subrings that eventually form in Figure 5.

Figure 11 shows the evolution ofthe dustrings as gas is gradually removed. Within the ring, the thickness of the dust layer is already limited by the vertical resolution, so initially the growth of dust density is due only to radial contraction. Considering the dominant inner subring, we find that the ring width w_d contracts as w_d ∝ (M_d/M_g)^−0.19, until we reach log(M_d/M_g) ≈ 0.4−0.5, after which further contraction is limited by the radial resolution.

Eventually, we also observe the formation of dust clumps. To quantify the growth rate of the clumping, we compute the ratio, at the midplane, of the maximum dust density within the ring to the azimuthally averaged dust density at the center of the ring. We find that this ratio scales as (M_d/M_g)^0.5, indicating that there are azimuthal asymmetries within the ring that grow faster than the average dust density and thereby signaling the formation of dust clumps. We note that clumps have not formed at the corresponding times in the simulation shown in Fig. 5, indicating that gas removal is primarily responsible for dust clumping. These clumps would eventually reach Hill density and as in the viscous case, they are expected to undergo gravitational collapse as a result.

In summary, we observed dust clumping in both viscous and inviscid cases, with the maximum dust density growing as (M_d/M_g)^0.35 and (M_d/M_g)^0.5, respectively. In the viscous case, the late-stage evolution becomes similar to that of the inviscid case as gas is progressively removed. This convergence arises because gas depletion drives the disk toward the inviscid limit: even for a fixed α-viscosity, the effective dust diffusivity is reduced by the factor ρ_g/(ρ_g + ρ_d). In both cases, the dust clumps are therefore expected to undergo gravitational collapse.

8 Conclusions

In this work, we reproduced and extended the 2D shearing-box results of Liu & Bai (2023) in a global 2D disk simulation. In a relatively mild pressure bump with the same parameters as those characterizing LB23, we recovered the type II DRWI and associated formation of dust clumps. By continuously supplying dust to the ring, we were able to extend their analysis by characterizing the evolution of these clumps, finding an average of 20-25 clumps and reaching a maximum dust-to-gas surface density ratio of Σ_d/Σ_b, ≈ 34.

In contrast, we find that the DRWI is suppressed in fully 3D disks with the same viscosity (i.e., α = 10⁻⁴ in our study, just as in LB23): no dust clump develops and the dust maintains a Gaussian distribution in r and z around the ring center. We attribute this suppression to dust sedimentation, which concentrates solids at the midplane, while leaving the upper and lower layers dust-poor and dynamically stable. Viscous coupling between vertical layers then acts to damp perturbations that arise at the midplane. This interpretation is supported by the unstratified 3D simulation, in which dust is present at all heights and the disk exhibits stronger perturbations reminiscent of the 2D behavior.

In inviscid 3D disks, we find that the dust ring evolves differently. Multiple subrings form because the back-reaction of the dust on the gas changes the gas density enough to form new pressure maxima around the original (single) maximum. When dust is continuously supplied, instead of enhancing the dust density in the already existing rings, it modifies the overall gas density profile, partially flattening and shifting the global pressure bump. This allows for the formation of more numerous subrings, over a more extended radial region. As a result, the dust density may remain everywhere lower than the threshold for gravitational instability. We stress that our code does not account for self-gravity, so the results are independent of the gas and dust densities, but dependent only on their ratio. Thus, the question of whether the Hill density is ultimately reached in such cases depends on the initial assumed density.

This result reveals that the radial width of the dust rings imaged by ALMA in protoplanetary disks can arise from the presence of multiple, closely spaced subrings that remain unresolved at observational resolutions. When interpreted as a single broad ring, the contrast with the vertical thinness of the dust layer suggests a much larger turbulent diffusion in the radial direction than in the vertical direction (Dullemond et al. 2018). Therefore, our results suggest that this inferred anisotropy of turbulence might not be real.

Finally, we investigate the late evolutionary stages of disks by modeling gas loss, while noting that a detailed photoevaporation model is not included. In inviscid disks undergoing gas depletion, we find growing azimuthal asymmetries within the dust rings, consistent with the formation of dust clumps. In viscous disks, the evolution eventually resembles the inviscid case, including clump formation, once sufficient gas has been removed; as gas is depleted, the effective α in the dust diffusion coefficient is reduced and, thus, the disk approaches the inviscid limit.

Overall, our results suggest that planetesimal formation within pressure bumps is favored in very low-viscosity disks and can be triggered during late stages of disk evolution, when sufficient gas has been removed by photoevaporation. The amount of gas removal required depends on the gas viscosity. This picture is consistent with the inferred late formation times of the parent bodies of chondritic meteorites in the Solar System.

Acknowledgements

AM and MT are grateful for support from the ERC advanced grant HolyEarth No. 101019380. This work was granted access to the HPC resources of IDRIS and CINES under the allocation A0180407233 made by GENCI. EL and MT wish to thank Alain Miniussi for maintenance and re-factorization of the code fargOCA.

References

All Tables

All Figures

Fig. 1 Snapshots in the 2D type II DRWI run of the dust density in the azimuthal-radial plane. The time is annotated on the top left of each panel. The panels are colored in a power-law scale and the dust surface density is plotted with respect to the initial background gas density profile (Eq. (4)). The panel at 3350 Porb is repeated to the right of the dotted black line, with the identified dust clumps circled in black.

In the text

	Fig. 2 Evolution of the dust clumps in the 2D type II DRWI run. Top panel: number of identified dust clumps at each time. Middle panel: maximum density of the dust clumps at each time, normalized by the average dust density at the clump’s radial location. Bottom panel: evolution of the minimum gas-to-dust ratio.
In the text

	Fig. 3 Same as Fig. 1, but depicting the volume density of dust on the midplane in a 3D simulation with the same pressure bump parameters and gas viscosity as in the 2D model of Section 4.
In the text

	Fig. 4 Same as Fig. 3, but for the 3D unstratified simulation.
In the text

	Fig. 5 Same as Fig. 3, but for the 3D inviscid (α = 0) simulation without dust inflow.
In the text

	Fig. 6 Radial profiles for the 3D inviscid simulation without dust inflow of the pressure (green), dust density (brown), and η (pink). The solid pink curve is the pressure gradient, while the dashed pink curve is the true value of η.
In the text

	Fig. 7 Comparison between the dust ring evolution in the 2D r – z inviscid disk versus the full 3D disk with the same radial and vertical resolution.
In the text

Fig. 8 Evolution of the high resolution 2D r – z inviscid disk with ϵ0 = 0.01 and dust inflow. In the top panel, we plot the r – z view of the dust to gas volume density at the final integration time t = 1500 Porb. The lower panels show the evolution with time versus r/r0 of, respectively: the vertically integrated density ratio, Σd/Σg (second panel), the η values computed from midplane pressure gradient (third panel) and the midplane dust’s radial velocity υr (bottom panel). In the third panel, a transition from blue to red indicates an η sign inversion and thus a dust trapping location.

In the text

Fig. 9 Vertical dust structure as a function of radius for the high-resolution 2D r – z inviscid disk simulation with ϵ0 = 0.01 and dust inflow at t = 1000 Porb. The top panel shows a broader radial region, while the bottom panel shows a zoomed-in view of the dust ring region. In the top panel, the outer edge of the dust ring is indicated with a dashed red line. Each set of panels (from left to right) shows: the dust density, radial gas velocity, vertical gas velocity (both normalized by the azimuthal Keplerian velocity), and dust density overlaid with the vertical gas velocity.

In the text

Fig. 10 Maximum dust density in the ring as a function of the dust-togas mass ratio during gas removal for both a 2D (orange) and 3D (blue) viscous disk (α = 10−4). The upper x-axis indicates orbital periods at the reference radius of r = 1. Linear best-fit lines are shown for three stages of the gas evaporation process, with their slopes indicated next to each line. The slopes of the first two stages are consistent between the 2D and 3D case and then diverge in the final stage, where there is clumping in 3D.

In the text

	Fig. 11 Same as Fig. 5, but showing the evolution of the dust density as gas is removed from the disk. The dust density is plotted with respect to the Hill density and the total dust-to-gas mass ratio at each given time is indicated at the top of the panel. The orbital period in the top-left corner is given relative to the start of the original inviscid run with no inflow in Section 6.1.
In the text