EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 710 (June 2026) A&A, 710 (2026) A21 Full HTML
Open Access
Issue
A&A
Volume 710, June 2026
Article Number A21
Number of page(s) 16
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/202659809
Published online 28 May 2026

© The Authors 2026

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model.

Open Access funding provided by Max Planck Society.

1 Introduction

Dust in protoplanetary disks plays a central role in planet formation and serves as one of the most direct observational tracers of disk dynamics. Observations with the Atacama Large Mil-limeter/submillimeter Array (ALMA) have revealed that dust emission commonly exhibits substructures, including annular gaps and rings (ALMA Partnership 2015; Andrews et al. 2016, 2018; Isella et al. 2016; Long et al. 2018, 2019; Cieza et al. 2019; Pérez et al. 2019; Bosschaart et al. 2026), central cavities (van der Marel et al. 2018, 2021, 2022; Pérez et al. 2020; Francis & van der Marel 2020; Benisty et al. 2021; Facchini et al. 2026), winding spiral arms (Pérez et al. 2016; Dong et al. 2018; Huang et al. 2018; Kurtovic et al. 2018; Rosotti et al. 2020; Yoshida et al. 2025), and non-axisymmetric crescents and clumps (Fukagawa et al. 2013; van der Marel et al. 2013; Isella et al. 2018; Casassus et al. 2019; Long et al. 2022; Ribas et al. 2024; Wölfer et al. 2025). The ubiquity of these features suggests that disks are intrinsically structured and dynamically active on planet-forming scales, potentially reflecting instabilities from coupled dust and gas. Because ALMA observations predominantly trace millimeter-sized grains, interpreting these substructures requires a detailed understanding of how marginally coupled dust is transported and concentrated within disks.

The dynamics of millimeter-sized grains is particularly important near the disk midplane, where dust back-reaction onto the gas is enhanced by vertical settling and where planetesimal formation is expected to initiate. However, grain growth beyond millimeter sizes faces several well-known barriers: rapid inward drift removes large grains on short timescales, high-velocity collisions lead to fragmentation, and inefficient sticking causes particles to bounce off rather than to grow (Weidenschilling 1977; Blum & Wurm 2008; Brauer et al. 2008; Birnstiel et al. 2010; Güttler et al. 2010; Zsom et al. 2010, 2011; Drazkowska et al. 2023; Birnstiel 2024). Overcoming these barriers requires mechanisms that locally retain and concentrate dust, for example through trapping at pressure maxima or via instability-driven accumulation processes (Pinilla et al. 2012; Lyra & Lin 2013; Dullemond et al. 2018; Lesur et al. 2023). The observed substructures can therefore serve as key signposts of regions where dust accumulates efficiently enough to trigger planetesimal formation.

A leading mechanism for concentrating marginally coupled dust into dense structures is the streaming instability (SI; Youdin & Goodman 2005; Johansen & Youdin 2007). The SI belongs to the broader class of resonant drag instabilities, which arise when dust streams through gas and the dust-gas drift resonates with a natural wave mode of the gas (Squire & Hopkins 2018a,b). Simulations demonstrate that SI can produce strong dust clumping and associated turbulence for millimeter-sized grains, and that its efficiency depends sensitively on the local dust-to-gas ratio and stopping time (Bai & Stone 2010; Carrera et al. 2015; Yang et al. 2017; Lim et al. 2024, 2025; Ostertag & Flock 2025). These properties make SI a natural bridge between disk dynamics and the onset of planetesimal formation.

Despite its effectiveness under favorable conditions, the growth and saturation of SI are sensitive to the gas dynamics. Turbulent diffusion and viscous damping can weaken or suppress SI by dispersing dust overdensities and disrupting the resonance that drives its growth (Chen & Lin 2020; Umurhan et al. 2020; Lim et al. 2024). This sensitivity highlights the need to consider SI in conjunction with other hydrodynamic processes that contribute to disk turbulence and angular momentum transport.

In regions with weak magnetic activity, hydrodynamic instabilities can provide an important source of turbulence. One such instability is the vertical shear instability (VSI), a baroclinic instability associated with the radial temperature gradient in disks with a short thermal relaxation time (Nelson et al. 2013; Barker & Latter 2015; Lin & Youdin 2015). Simulations show that VSI generates turbulence and large-scale vertical gas motions that originate in the surface layers and extend toward the midplane, thereby modifying dust settling, concentration, and the conditions under which SI develops (Stoll & Kley 2014, 2016; Lin & Youdin 2017; Lin 2019).

Recent studies have begun to explore the coexistence of SI and VSI in global disk models, revealing that dust back-reaction and hydrodynamic turbulence can mutually influence instability growth, saturation, and dust-layer structure. Using two-dimensional axisymmetric disk models with dust treated as Lagrangian particles, Schäfer et al. (2020) showed that the outcome of SI-VSI interaction depends sensitively on their relative growth and saturation timescales: when VSI develops first, VSI-driven turbulence dominates the dust layer and inhibits settling, whereas when SI grows concurrently, dust back-reaction enhances midplane turbulence and produces denser dust concentrations. Building on this work, Schäfer & Johansen (2022) showed that VSI-induced pressure perturbations can act as seeds for SI, enabling strong dust clumping at lower dust-to-gas ratios and smaller grain sizes than would be required for SI operating in isolation. More recently, Schäfer et al. (2025) quantified dust diffusion driven by SI and VSI both separately and jointly, showing that VSI produces highly anisotropic, scale-dependent diffusion associated with large-scale motions, while SI contributes an intrinsic diffusion level that broadens the midplane dust layer. In parallel, Huang & Bai (2025a,b) investigated SI-VSI coupling using both two- and three-dimensional multi-fluid simulations, highlighting the emergence of zonal flows, vortices, and secondary instabilities that further regulate dust concentration. In particular, Huang & Bai (2025b) report that the outcome of SI-VSI interaction is independent of whether SI or VSI grows first, contrary to Schäfer et al. (2020). Collectively, these studies demonstrate that SI and VSI can mutually influence each other’s saturation, turbulence properties, and dust morphology, motivating further exploration of their coupled behavior under more general disk conditions.

Existing studies of SI-VSI coexistence have largely focused on inviscid conditions. However, in realistic disk environments, SI and VSI are unlikely to be the only sources of angular momentum transport: disk winds, nonideal magnetohydrodynamic effects, and other processes can also contribute to the background viscosity (Bai & Stone 2013; Bai 2013; Bai et al. 2016; Béthune et al. 2017; Lesur et al. 2023). Observational constraints suggest that α is approximately between 10−4 and 10−3, encompassing the processes discussed above (Rosotti 2023). Such additional turbulence can modify the development of instabilities; in particular, VSI is sensitive to turbulent diffusion, which can weaken its growth or shift unstable modes to smaller scales (Barker & Latter 2015; Lin & Youdin 2015). Further investigation is therefore required to better understand how these instabilities operate under more realistic disk conditions.

Motivated by these considerations, we investigated the coupled evolution of SI and VSI in viscous protoplanetary disks using two-dimensional axisymmetric models. We focused on the dynamics and morphology of marginally coupled millimetersized dust near the midplane. Our study explores how viscosity modifies the coupled evolution of SI and VSI, how dust substructures emerge under these conditions, and what morphological signatures can arise from the combined action of drag-driven clumping and vertical-shear turbulence. By extending SI-VSI studies into viscous regimes, our work provides a step toward linking instability-driven dust dynamics with the diversity of substructures observed in protoplanetary disks.

The structure of this work is as follows. We first describe the numerical methods and simulation setup in Sect. 2, and present our results on dust dynamics and substructure formation in viscous environments in Sect. 3. We discuss the implications of our findings for disk observations in Sect. 4, and conclude with a summary of our main results in Sect. 5.

2 Method

We considered a two-dimensional, axisymmetric protoplanetary disk model composed of gas and dust orbiting a central star of mass M*. We use {r, θ} to denote the spherical radius and polar angle, and {R, Z} for the cylindrical radius and height. Both coordinate systems are centered on the star. The superscript “ini” denotes initial values. The subscript “0” denotes evaluations in the midplane at the reference radius (R = R0, Z = 0).

2.1 Model of gas

We assumed an isothermal equation of state P=ρgcs2,Mathematical equation: P = \rho_{\rm g}c_{\rm s}^2,(1)

where {P, ρg, cs} are the pressure, volumetric density, and sound speed of the gas, respectively. We assumed a vertically isothermal, time-independent temperature profile Tg(R)∝(R/R0)−q.Mathematical equation: T_{\rm g}(R) \propto (R/R_0)^{-q}.(2)

Here, q = 3/7 is adopted from Chiang & Goldreich (1997), leading to a moderately flared disk with an aspect ratio of the pressure-supported gas disk given by hg(R)≡Hg/R=hg0(R/R0)2/7.Mathematical equation: h_{\rm g}(R) \equiv H_{\rm g}/R = h_{\rm g0}(R/R_0)^{2/7}.(3)

In the above, we assume hg0 = 0.1, Hg(R)=cs/ΩKMathematical equation: H_{\rm g}(R) = c_{\rm s}/\Omega_{\rm K}(4)

is the pressure scale height of the gas, ΩK(R)=vK/R=GM⋆/R3Mathematical equation: \Omega_{\rm K}(R) = v_{\rm K}/R = \sqrt{GM_\star / R^3}(5)

is the Keplerian angular velocity, and G is the gravitational constant. The disk’s self-gravity and magnetic fields are neglected.

2.2 Model of dust

We considered a single species of dust modeled as a pressureless fluid. The dynamical coupling between gas and dust is parameterized by the Stokes number St=τsΩK,Mathematical equation: {\rm St} = \tau_{\rm s}\Omega_{\rm K},(6)

where τs is the particle stopping time characterizing the frictional drag force between gas and dust. For dust grains that remain marginally coupled to the gas, we assumed the gasdust interaction is in the Epstein regime. Assuming a fixed size (s) and internal density (ρo), the stopping time then becomes (Weidenschilling 1977) τs=ρ∘s/ρgcs,Mathematical equation: \tau_{\rm s} = \rho_\circ s / \rho_{\rm g}c_{\rm s},(7)

and in practice, it is prescribed by τs=ρg0iniρgcs0csSt0iniΩK0.Mathematical equation: \tau_{\rm s} = \frac{\rho_{\rm g0}^{\rm ini}}{\rho_{\rm g}}\frac{c_{\rm s0}}{c_{\rm s}}\frac{{\rm St}_{\rm 0}^{\rm ini}}{\Omega_{\rm K0}}.(8)

We adopted a fiducial reference Stokes number St0ini=10−3Mathematical equation: ${\rm St}_0^{\rm ini} = 10^{-3}$, which corresponds to approximately 100-μm-sized grains at R = 100 au in protoplanetary disks such as HL Tau1,2. We also explore the effect of varying the Stokes number in Sect. 3.5.

2.3 Basic equations

The hydrodynamic equations for gas and dust are given by ∂ρg∂t+∇⋅(ρgVg)=0Mathematical equation: &\frac{\partial \rho_{\rm g}}{\partial t} + \nabla \cdot (\rho_{\rm g} \textit{\textbf{V}}_{\rm g}) = 0 \\[0.5em](9) ∂Vg∂t+Vg⋅∇Vg=−∇Φ+ϵτs(Vd−Vg)−1ρg∇P+1ρg∇⋅TMathematical equation: &\begin{aligned} \frac{\partial \textit{\textbf{V}}_{\rm g}}{\partial t} + \textit{\textbf{V}}_{\rm g} \cdot \nabla \textit{\textbf{V}}_{\rm g} = - &\nabla \Phi + \frac{\epsilon}{\tau_{\rm s}}(\textit{\textbf{V}}_{\rm d} - \textit{\textbf{V}}_{\rm g}) \\[0.5em] - &\frac{1}{\rho_{\rm g}} \nabla P + \frac{1}{\rho_{\rm g}} \nabla \cdot \mathcal{T} \end{aligned}(10) ∂ρd∂t+∇⋅(ρdVd)=0Mathematical equation: &\frac{\partial \rho_{\rm d}}{\partial t} + \nabla \cdot (\rho_{\rm d} \textit{\textbf{V}}_{\rm d}) = 0 \\[0.5em](11) ∂Vd∂t+Vd⋅∇Vd=−∇Φ−1τs(Vd−Vg).Mathematical equation: &\frac{\partial \textit{\textbf{V}}_{\rm d}}{\partial t} + \textit{\textbf{V}}_{\rm d} \cdot \nabla \textit{\textbf{V}}_{\rm d} = - \nabla \Phi - \frac{1}{\tau_{\rm s}}(\textit{\textbf{V}}_{\rm d} - \textit{\textbf{V}}_{\rm g}).(12)

Here, {ρg, Vg} and {ρd, Vd} are the volumetric density and velocity vector of gas and dust, respectively. Φ(r)=−GM⋆/rMathematical equation: \Phi(r) = -GM_\star/r(13)

is the gravitational potential from the star, and ϵ≡ρd/ρgMathematical equation: \epsilon \equiv \rho_{\rm d}/\rho_{\rm g}(14)

is the dust-to-gas density ratio. The viscous stress tensor is given by T=ρgν[∇Vg+(∇Vg)†−23I∇⋅Vg],Mathematical equation: \mathcal{T} = \rho_{\rm g}\nu\left[\nabla\textit{\textbf{V}}_{\rm g} + (\nabla\textit{\textbf{V}}_{\rm g})^{\dagger} - \frac{2}{3}\textit{\textbf{I}}\,\nabla \cdot\textit{\textbf{V}}_{\rm g}\,\right],(15)

where I is the identity tensor, and the kinematic viscosity ν is parameterized as ν=αcsHg.Mathematical equation: \nu = \alpha c_{\rm s} H_{\rm g}.(16)

In our models, α is a dimensionless parameter for viscosity (Shakura & Sunyaev 1973) that is spatially constant unless otherwise specified (e.g., in Sect. 3.1). The value of α varies between models, and we primarily focused on disk models with α between 10−5 and 10−4. But we also explore the effect of varying α in Sect. 3.3. Dust back-reaction onto the gas is included in our models, while dust diffusion in the turbulent gas is intentionally neglected in order to isolate the dust stirring effect driven by dust back-reaction3.

2.4 Gas and dust initialization

The volumetric density of the gas (ρg) is initialized using the vertically integrated surface density profile Σgini(R)=Σg0ini[(R/R0)+(R/Rc)5]−1,Mathematical equation: \Sigma_{\rm g}^{\rm ini}(R) = \Sigma_{\rm g0}^{\rm ini}\left[(R/R_0)+(R/R_{\rm c})^5\right]^{-1},(17)

where Σg0iniMathematical equation: $\Sigma_{\rm g0}^{\rm ini}$ can be arbitrary for our non-self-gravitating disk models, and the corresponding volumetric density is ρgini=Σgini2πHg×exp⁡(−Z22Hg2).Mathematical equation: \rho_{\rm g}^{\rm ini} = \frac{\Sigma_{\rm g}^{\rm ini}}{\sqrt{2\pi}H_{\rm g}} \times \exp{\left(-\frac{Z^2}{2H_{\rm g}^2}\right)}.(18)

The last term of Eq. (17) provides a steeper decrease in the gas density profile beyond a cut-off radius Rc, which is inspired by recent rotation curve studies of protoplanetary disks (e.g., Martire et al. 2024). With Rc = 2R0, the deviation between Eq. (17) and an ordinary power-law profile Σg ∝ R−1 becomes noticeable (about a factor of 2) starting from 2R0, and reaches a factor of roughly 20 at 5R0.

The volumetric density of the dust (ρd) is initialized via the dust-to-gas density ratio ϵini≡ρdini/ρgini=ϵ0ini×δR×δZ.Mathematical equation: \epsilon^{\rm ini} \equiv \rho_{\rm d}^{\rm ini}/\rho_{\rm g}^{\rm ini} = \epsilon_0^{\rm ini} \times \delta_R \times \delta_Z.(19)

Here, δR is a radial-direction taper given by δR={(11−ζR−RminR0−Rmin)2forR<[R0−ζ(R0−Rmin)],(11−ζR−RmaxR0−Rmax)2forR>[R0−ζ(R0−Rmax)],1.0for anywhere else,Mathematical equation: \delta_R = \begin{cases} \left(\dfrac{1}{1-\zeta}\dfrac{R - R_{\rm min}}{R_0 - R_{\rm min}}\right)^2 \\[1.5em] \qquad\qquad\qquad {\rm for}\; R < [R_0-\zeta(R_0 - R_{\rm min})], \\[0.5em] \left(\dfrac{1}{1-\zeta}\dfrac{R - R_{\rm max}}{R_0 - R_{\rm max}}\right)^2 \\[1.5em] \qquad\qquad\qquad {\rm for}\; R > [R_0-\zeta(R_0 - R_{\rm max})], \\[0.5em] 1.0 \qquad\qquad\;\:\:\text{for anywhere else}, \end{cases}(20)

where Rmin and Rmax are the inner and outer radial boundaries of the computational domain (specified in Sect. 2.5), and ζ = 3/4, such that ρdni gradually becomes zero when approaching the radial boundaries. Meanwhile, δZ is a vertical-direction taper given by δZ=exp⁡(−Z2Hg2−Hd22Hg2Hd2),Mathematical equation: \delta_Z = \exp{\left(-Z^2\frac{H_{\rm g}^2 - H_{\rm d}^2}{2H_{\rm g}^2 H_{\rm d}^2}\right)},(21)

such that ρdni retains a vertically Gaussian profile with a dust scale height Hdini=0.1HginiMathematical equation: $H_{\rm d}^{\rm ini} = 0.1H_{\rm g}^{\rm ini}$. We adopted this pre-settled dust layer to reduce computational cost by bypassing the initial dustsettling phase, which is not expected to affect the subsequent dust evolution in our models. With a fiducial value ϵ0ini=1Mathematical equation: $\epsilon_0^{\rm ini} = 1$, this setup corresponds to a disk model with a total dust-to-gas mass ratio (i.e., metallicity) of Z ≈ 0.1. We also explore the effect of varying Z in Sect. 3.4.

The azimuthal velocities are initialized to Vg,ϕini=RΩK(R)1−2η,Mathematical equation: & V_{\rm g,\phi}^{\rm ini} = R\Omega_{\rm K}(R)\sqrt{1-2\eta},(22) Vd,ϕini=rΩK(r),Mathematical equation: & V_{\rm d,\phi}^{\rm ini} = r\Omega_{\rm K}(r),(23)

where η=12[(p+q)hg2+q(1−Rr)]Mathematical equation: \eta = \frac{1}{2}\left[(p+q)h_{\rm g}^2 + q\left(1 - \frac{R}{r}\right)\right](24)

is a term from the radial gradient of gas pressure, with p=−dlogρg|Z=0dlogR=−dlogΣgdlogR+dlogHgdlogR=(R/R0)+5(R/Rc)5(R/R0)+5(R/Rc)5+97=5(R/R0)4+255(R/R0)4+25+97.Mathematical equation: \begin{split} p & = -\frac{{\rm d\,log}\,\rho_{\rm g}|_{Z=0}}{{\rm d\,log}\,R} = -\frac{{\rm d\,log}\,\Sigma_{\rm g}}{{\rm d\,log}\,R} + \frac{{\rm d\,log\,{\it H}_{g}}}{{\rm d\,log}\,R} \\[0.5em] & = \frac{(R/R_0)+5(R/R_{\rm c})^5}{(R/R_0)+\phantom{5}(R/R_{\rm c})^5} + \frac{9}{7} \\[0.5em] & = \frac{5(R/R_0)^4 + 2^5}{\phantom{5}(R/R_0)^4 + 2^5} + \frac{9}{7}. \end{split}(25)

Both radial and vertical velocities of gas and dust are initialized to zero. We note that this initialization does not correspond to a strict steady state of the coupled gas-dust system (Kanagawa et al. 2017). In the absence of explicit dust diffusion, the vertical equilibrium of the dust layer cannot be maintained, as dust continues to settle and therefore immediately deviates from the initial state once the simulation begins. Consequently, even an exact steady-state initialization would not be preserved. We therefore adopted this setup as a convenient starting configuration and do not expect the velocity initialization to significantly affect the subsequent dust evolution.

2.5 Numerical tool

Our disk models are evolved by the multi-fluid hydrodynamic code FARGO3D (Benítez-Llambay & Masset 2016; Benítez-Llambay et al. 2019). We adopted an axisymmetric domain in spherical coordinates, spanning 0.1R0 to 5.0R0 in r, and between π/2 ± atan(2hg0) in θ (i.e., roughly ±2Hg in the vertical direction). Our fiducial grid resolution is Nr × Nθ = 8,000 × 800 with logarithmic spacing in r and uniform spacing in θ. This corresponds to approximately 200 cells per Hg in both directions. We also explore the effect of varying the grid resolution in Appendix A.

We adopted a wave-killing (i.e., Stockholm) boundary condition for gas in the radial direction. The inner (respectively, outer) damping zone ends (begins) at the radius where the Keplerian velocity is 1.5 times slower (faster) than that of the inner (outer) radial boundary. At the boundaries, the gas density is assumed to be in hydrostatic equilibrium, and the dust density is assumed to be symmetric. The meridional velocities of both gas and dust are assumed to be zero at the boundaries, except that the dust is allowed to leave the computational domain through the inner radial boundary.

3 Result

In this section we first describe the dust dynamics in two of our fiducial setups and then explore a broader parameter space to understand how dust dynamics and the resulting substructures depend on dust and gas properties.

3.1 Dust distribution in the model with α = 10−4

Figure 1 shows the distribution of dust-to-gas density ratio at 400 reference orbital periods4 in the disk model with α = 10−4 and St0 = 10−3. For the first time, we demonstrate that strong dust back-reaction from marginally coupled grains in the midplane can lead to the formation of spatially episodic, shuttlecockshaped dust substructures that feature a denser, more vertically settled “head”, and a less dense but more vertically extended “tail”.

The head of the shuttlecock exhibits a peak dust-to-gas density ratio of typically 1 ≲ max(ρd/ρg) ≲ 3, and therefore does not qualify as strong clumping. This is consistent with expectations, as our adopted metallicity Z = 0.1 lies below the threshold required for strong clumping due to SI (Lim et al. 2024): log⁡Zcrit(α,τs)=0.15×log⁡α×log⁡α+1.18×log⁡α−0.24×log⁡α×log⁡τs−1.48×log⁡τs.Mathematical equation: \begin{split} \label{eq:z_crit} \log \mathcal{Z}_{\rm crit}(\alpha, \tau_{\rm s}) = \: &0.15\times\log\alpha\times\log\alpha \,+ 1.18\times\log\alpha \\ - \: &0.24\times\log\alpha\times\log\tau_{\rm s} - 1.48\times\log\tau_{\rm s}. \end{split}(26)

For dust grains with St = 3 × 10−3 at R = 2R0, the critical value is Zcrit ≳ 0.17 in our setup5.

To quantify the vertical extent of the dust layer, we used the mass-weighted root mean square vertical height Hdrms=∑(mdZ2)∑mdMathematical equation: H_{\rm d}^{\rm rms} = \sqrt{\frac{\sum(m_{\rm d}Z^2)}{\sum m_{\rm d}}}(27)

for the dust scale height relative to the midplane, and the weighted standard deviation Hdstd=∑(mdZ2)∑md−(∑(mdZ)∑md)2Mathematical equation: H_{\rm d}^{\rm std} = \sqrt{\frac{\sum(m_{\rm d}Z^2)}{\sum m_{\rm d}} - \left(\frac{\sum(m_{\rm d}Z)}{\sum m_{\rm d}}\right)^2}(28)

for the dust scale height relative to the center of the dust layer. The measurements are shown in Fig. 2a. We find that Eqs. (27) and (28) yield similar results, meaning that the dust layer is close to being centered at the disk midplane. We further performed fits with a Gaussian centered at the midplane, but the result deviates from HdrmsMathematical equation: $H_{\rm d}^{\rm rms}$ and HdstdMathematical equation: $H_{\rm d}^{\rm std}$, meaning that the dust layer does not match a Gaussian profile. Therefore, we adopted HdstdMathematical equation: $H_{\rm d}^{\rm std}$ as our fiducial measurement of dust scale height hereafter. Without explicit dust diffusion, the dust scale height reaches Hd/Hg ≳ 0.15 at the head of the shuttlecocks and Hd/Hg ≳ 0.25 at the tail. Notably, these ratios are comparable to, or even above, the expected value resulting from the balance between vertical settling and turbulent diffusion (e.g., Dubrulle et al. 1995; Youdin & Lithwick 2007; see Fig. 2): HdHg=(1+Stα1+2St1+St)−1/2.Mathematical equation: \frac{H_{\rm d}}{H_{\rm g}} = \left(1 + \frac{\rm St}{\alpha}\frac{\rm 1+2St}{\rm 1+St}\right)^{-1/2}.(29)

However, we emphasize that this value should be regarded only as a reference estimate, since explicit dust diffusion is not included in our models and a true settling-diffusion equilibrium is therefore not expected to be established.

Why is the substructure shuttlecock-shaped?

The streamlines of gas velocity are overplotted in the bottom panel of Fig. 1. The streamlines exhibit radially alternating, vertically extended gas flows that appear to compress the head toward the midplane and stretch the tail in the opposite direction. However, these alternating gas flows may either cause the shuttlecock or be its consequence. On the one hand, the gas flows might be the breathing modes of VSI that are symmetric about the midplane, and thus the cause of the shuttlecock, provided that the midplane-crossing corrugation modes are suppressed by dust back-reaction (Lin & Youdin 2017; Lin 2019; Schäfer et al. 2020; Huang & Bai 2025a). On the other hand, the gas flows might be the result of SI in the midplane, as has also been observed in Fig. 9 of Schäfer et al. (2020), and thus the consequence. To distinguish between these two scenarios, we first performed one control simulation with a higher viscosity of α = 10−3, which was accompanied by two additional runs in which α was allowed to vary between 10−4 and 10−3 along the vertical direction. The results are shown in Fig. 3.

In the models where α varies spatially, the vertical profiles of α are given by α1(θ)=1(9ψθ−8)×10−4Mathematical equation: \alpha_1(\theta) &= \phantom{1}(9\psi_\theta-8) \times 10^{-4} \\(30) α2(θ)=(19−9ψθ)×10−4,Mathematical equation: \alpha_2(\theta) &= (19-9\psi_\theta) \times 10^{-4},(31)

where ψθ=ψ+(θ)+ψ−(θ)∈(1,2)Mathematical equation: $\psi_\theta = \psi^+(\theta) + \psi^-(\theta) \in (1,2)$ is a window function connected by two mirrored logistic sigmoids ψ+(θ)=1−[1+exp⁡(−θ−π/2+10δθδθ)]−1Mathematical equation: \psi^+(\theta) = \phantom{1\,-} &\left[1 + \exp{\left(-\frac{\theta-\pi/2+10\delta_\theta}{\delta_\theta}\right)}\right]^{-1} \\(32) ψ−(θ)=1−[1+exp⁡(−θ−π/2−10δθδθ)]−1,Mathematical equation: \psi^-(\theta) = 1\,- &\left[1 + \exp{\left(-\frac{\theta-\pi/2-10\delta_\theta}{\delta_\theta}\right)}\right]^{-1},(33)

with δθ = arctan(0.1 × hg0). This enables α to transition smoothly within a vertical range around Z ∼ Hg (see the last panel in Fig. 3), thereby not introducing sharp gradients of α close to the midplane that would impact dust dynamics there.

With a globally α = 10−3 prescription (Fig. 3a), the shuttlecock-shaped substructure disappears as a result of the suppression of both SI and VSI. The shuttlecock structure does not reappear when α decreases to 10−4 only in the disk atmosphere (panel b), indicating that VSI, if present in the globally α = 10−4 model (panel d), is unlikely to be responsible for the formation of the shuttlecock. In fact, as we show in Fig. 5 and Appendix B, VSI remains inactive when α ≥ 10−4 in the disk atmosphere. In contrast, the shuttlecock re-emerges when α decreases to 10−4 in the midplane while remaining large in the atmosphere (panel c). This behavior suggests that the shuttlecock is the result of SI, i.e., strong dust back-reaction from marginally coupled dust grains in the midplane.

Thumbnail: Fig. 1 Refer to the following caption and surrounding text. Fig. 1

Dust-to-gas density ratio at t = 400 T0 from the model with α = 10−4 and St0 = 10−3. Dashed curves in the upper panel denote contours for 25% of the gas scale height, and the arrow below denotes the Stokes number of dust grains in the midplane. Gas streamlines are overplotted in the lower panel, with the line color denoting the vertical gas velocity.
Thumbnail: Fig. 2 Refer to the following caption and surrounding text. Fig. 2

Radial profiles of the dust-to-gas scale height ratio at t = 400 T0 from models with α = 10−4 (panel a) and α = 10−5 (panel b), both with St0 = 10−3. The dust scale heights are measured using the root mean square (HdrmsMathematical equation: $H_{\rm d}^{\rm rms}$, Eq. (27)), the standard deviation (HdstdMathematical equation: $H_{\rm d}^{\rm std}$, Eq. (28)), or by fitting a Gaussian profile. The dashed curve indicates a reference-only expected ratio assuming explicit dust diffusion, as given by Eq. (29). The radial extent is the same as in Figs. 1 and 4.
Thumbnail: Fig. 3 Refer to the following caption and surrounding text. Fig. 3

Dust-to-gas density ratio at t = 400 T0 from models with different vertical prescriptions of the α parameter. Panel a : model in which α = 10−3. Panel b : model in which α decreases from 10−3 in the midplane to 10−4 in the atmosphere. Panel c: model in which α increases from 10−4 in the midplane to 10−3 in the atmosphere. Panel d: model in which α = 10−4. Last panel: Corresponding vertical profiles of α, given by Eqs. (30)-(33). Dashed curves in panels b and c denote contours for one gas scale height. The shaded regions indicate where the transition of α occurs and are bounded by aspect ratios evaluated at the inner and outer radial boundaries.

3.2 Dust distribution in the model with α = 10−5

Similar to Fig. 1, Fig. 4 presents the dust-to-gas density ratio, but for the model with α = 10−5. In contrast to the higher-viscosity case, the shuttlecock-shaped substructures are absent in this model. Instead, the lower viscosity allows clearer signatures of VSI to emerge, particularly its corrugation modes, as illustrated by the streamlines in panel b. The vertical gas motions reach typical Mach numbers of around 0.1, comparable to those reported for VSI-driven flows in inviscid models (Schäfer et al. 2020; Huang & Bai 2025b). These vertical motions corrugate the dust layer about the midplane, leading to discrepancies between HdrmsMathematical equation: $H_{\rm d}^{\rm rms}$ and HdstdMathematical equation: $H_{\rm d}^{\rm std}$ in Fig. 2b.

The presence of VSI is further supported by the distribution of gas specific angular momentum shown in Fig. 5. In the α = 10−5 model, pronounced vertical bands are evident in the disk’s atmosphere, consistent with VSI-driven angular momentum mixing (Melon Fuksman et al. 2024). In contrast, such features are absent in the α = 10−4 model, supporting our interpretation that the VSI is suppressed at higher viscosity.

Despite the lower viscosity, which provides more favorable conditions for SI, the maximum dust-to-gas density ratio typically remains below 3, although our adopted metallicity is above the threshold Zcrit(α,τs) = 0.048. Nevertheless, compared with the α = 10−4 model, a larger fraction of the disk midplane remains dense in the α = 10−5 case. This suggests that even a midplane with ρd/ρg ≳ 1 in a low-viscosity environment is insufficient to suppress the corrugation modes of the VSI.

In terms of the dust-layer thickness measured by HdstdMathematical equation: $H_{\rm d}^{\rm std}$, the α = 10−5 model yields Hd/Hg ≳ 0.25 (see Fig. 2b), comparable to the α = 10−4 case. However, this value is significantly larger than the prediction from Eq. (29) for α = 10−5.

3.3 Dust distribution in models with different α

Given the distinct dynamical behaviors identified in Sects. 3.1 and 3.2, we now compare the dust distributions across these models, together with an additional run adopting an even lower viscosity of α = 10−6. The comparison is presented in Fig. 6.

The results indicate that shuttlecock-shaped substructures only arise in models where SI dominates the dust dynamics in an otherwise relatively undisturbed disk midplane. These structures disappear in higher-viscosity models where SI is suppressed, as well as in lower-viscosity models where active VSI interacts with SI in the midplane.

For models with α ≤ 10−5, the dust layer exhibits a qualitatively different morphology from those in higher-α models: the midplane becomes vertically corrugated due to VSI, while smaller-scale dust clumps produced by SI are still present. The primary distinction between panels (c) and (d) is the increased prominence of vertical dust streaks in the lowest-viscosity model, reflecting the stronger VSI activity with lower α.

3.4 Dust distribution in models with different metallicities

In our fiducial models, we adopted a relatively high metallicity of Z = 0.1 to ensure strong dust back-reaction. However, such an elevated metallicity may not be typical of protoplanetary disks. It is therefore important to assess how dust dynamics and the resulting substructures depend on the dust load. To this end, we performed two additional simulations at a lower metallicity, Z = 0.01. The results are shown in Fig. 7.

A comparison between Figs. 7a and 7c shows that the shuttlecock-shaped substructures persist even when the dust load is reduced to Z = 0.01. However, the maximum dust-to-gas density ratio barely reaches unity, reflecting the lower metallicity. In addition, the vertical thickness of the dust layer decreases due to weaker SI, resulting from reduced dust back-reaction.

Comparing panels b and d, we find that while SI dominates the midplane dust dynamics in the Z = 0.1 case, the weaker back-reaction at Z = 0.01 allows VSI to dominate instead. As a result, the settled dust layer at the midplane is strongly disrupted and effectively destroyed.

Thumbnail: Fig. 4 Refer to the following caption and surrounding text. Fig. 4

Similar to Fig. 1, but from the model with α = 10−5. The colorbar for gas streamlines has been adjusted to fit stronger vertical gas flows in this lower-viscosity disk model.
Thumbnail: Fig. 5 Refer to the following caption and surrounding text. Fig. 5

Specific angular momentum of gas at t = 400 T0 from models with α = 10−4 (panel a) and α = 10−5 (panel b), both with St0 = 10−3. The radial extent is the same as in Figs. 1 and 4.

3.5 Dust distribution in models with different Stokes numbers

So far, we have explored dust dynamics by varying α and Z, the two key parameters entering the SI-clumping threshold Zcrit(α, τs). However, we have not yet examined how dust dynamics and the resulting substructures depend on the grains’ Stokes number (i.e., stopping time). To this end, we performed two additional simulations with St0 = 10−2. The results are shown in Fig. 8, while measurements of the dust scale height and the maximum dust-to-gas density ratio are presented in Fig. 9.

A comparison between Figs. 8a and 8c shows that the shuttlecock-shaped substructures persist when the Stokes number increases to St0 = 10−2. However, their vertical extent is significantly reduced due to the more efficient settling of larger grains. Nevertheless, the tail of the shuttlecock still exhibits a dust scale height larger than the value predicted by Eq. (29) (see Fig. 9a). With the adopted metallicity Z = 0.1 exceeding the threshold Zcrit(α,τs) = 0.054 for St = 1.4 × 10−2 at R = 1.3R0, and given the enhanced settling of larger grains, the maximum dust-to-gas density ratio increases to typically 3 ≲ max(ρd/ρg) ≲ 10. However, this still falls short of the strong clumping regime.

Comparing Figs. 8b and 8d, the clumping efficiency becomes substantially higher, and the dust concentrations are more pronounced. Owing to the intrinsic turbulence generated by active SI in the midplane, the dust scale height in the α = 10−5 model is globally larger than the value predicted by Eq. (29) (see Fig. 9b). In contrast to the smaller-grain cases, the midplane dust layer does not exhibit a prominent VSI-driven corrugated morphology, as larger grains are less tightly coupled to the gas. Furthermore, since the adopted metallicity Z = 0.1 is well above the corresponding threshold Zcrit(α,τs) = 0.026 for St = 1.4 × 10−2 at R = 1.3R0, SI operates in a strongly nonlinear regime. The maximum dust-to-gas density ratio reaches typical values of 20 ≲ max(ρd/ρg) ≲ 50, potentially exceeding the Hill density of the disk (Klahr & Schreiber 2020): ρHill=94πM⋆R3.Mathematical equation: \rho_{\rm Hill} = \frac{9}{4\pi}\frac{M_\star}{R^3}.(34)

For dust clumps at 100 au orbiting a solar-mass star, pHill = 4.3 × 10−13g/cm3. Given a typical midplane gas density of ρg = 1.6 × 10−14g/cm3 (corresponding to Σg = 6.0g/cm2 and Hg = 10 au), this implies that a dust-to-gas density ratio ρd/ρg ≳ 27 is required for gravitationally bound dust clumps. Our measured peak values therefore approach, and may exceed, this threshold. Meanwhile, we note this result does not contradict the conclusion of Ostertag & Flock (2025), who found that SI struggles to produce clumps exceeding the Hill density. The difference mainly reflects the orbital radii considered: the clumps in that work were studied near 10 au, whereas our clumps are located at 100 au. Since ρHill ∝ R−3 increases inward more steeply than the disk gas density, much stronger dust concentration is required at smaller radii. Using our disk model, a dust-to-gas density ratio of ρd/ρg ≳ 139 would be required at 10 au to reach the Hill density. The two results are therefore consistent and mainly reflect the strong radial scaling of the Hill density.

Thumbnail: Fig. 6 Refer to the following caption and surrounding text. Fig. 6

Dust-to-gas density ratio at t = 400 T0 from models with different α values. Dashed curves denote contours for 25% of the gas scale height.

4 Discussion

4.1 Limitations of the numerical setup

The VSI arises in baroclinic disks where the angular velocity of gas varies with height owing to the radial temperature gradient (Nelson et al. 2013). Studies have shown that the fastestgrowing VSI modes are typically localized in the disk’s surface layers, where the vertical shear is strongest, before extending toward the midplane (Stoll & Kley 2014; Barker & Latter 2015). Because of this intrinsically vertical character, previous numerical studies commonly adopt vertically extended domains, often spanning several gas scale heights (e.g., ±4-5Hg), to capture both the surface and body modes. The treatment of vertical boundary conditions is also known to influence VSI development, as wave reflection and transmission at the disk surfaces can modify growth rates and the nonlinear saturation amplitude (Wu et al. 2024).

In our simulations, we adopted a vertically limited computational domain together with reflective vertical boundary conditions chosen for computational efficiency and numerical stability. We acknowledge that this restricted vertical extent may limit the VSI’s full development, particularly for the fastestgrowing surface modes. As a result, the saturated VSI amplitude, and hence the level of VSI-driven vertical stirring, may be underestimated in our models. This limitation should be borne in mind when interpreting the competition between VSI and SI and its influence on dust morphology. Future work employing more vertically extended domains and carefully designed boundary conditions will be necessary to fully capture the global VSI dynamics and its coupling to dust evolution.

We would also want to discuss the potential role of numerical diffusion in shaping the shuttlecock-shaped substructures identified in Sect. 3.1. Through private communications (e.g., Ostertag et al., in prep), we learned that SI in inviscid disk simulations evolved with relatively diffusive Riemann solvers can produce qualitatively similar dust morphologies, raising the question of whether the shuttlecock features might be numerical artifacts. We argue, however, that numerical diffusion is unlikely to dominate the behavior in our models. First, the vertical structure of the shuttlecocks, including both the head and the tail, exhibits a dust scale height larger than the reference value predicted by Eq. (29) (see Fig. 2) and is resolved by ≳100 numerical cells. In the radial direction, the separation between neighboring heads also exceeds the dust scale height. Such spatial scales are therefore not expected to be strongly affected by numerical diffusion.

Moreover, as shown in Appendix D, our cross-code comparison using the PLUTO code with both the Roe Riemann solver (Roe 1981; Mignone et al. 2007), one of the least diffusive solvers commonly used (Toro 2009), and the more diffusive Harten-Lax-van-Leer (HLL) solver yields results consistent with those obtained with FARGO3D. This agreement, independent of the choice of Riemann solver, further supports the robustness of the shuttlecock structures.

We nevertheless acknowledge another limitation of our numerical setup: although gas viscosity is included, an explicit dust diffusion term is not implemented. This is because the public versions of all FARGO3D, PLUTO, and ATHENA++ codes currently lack a dust module that can simultaneously incorporate dust diffusion and consistently account for dust back-reaction onto the gas. Future studies would benefit from numerical frameworks capable of treating both processes self-consistently, such as the recently developed but not yet open-sourced dust fluid module for ATHENA++ (Huang & Bai 2022), which will enable more rigorous investigations of diffusion effects on coupled gas-dust dynamics.

Thumbnail: Fig. 7 Refer to the following caption and surrounding text. Fig. 7

Dust-to-gas density ratio at t = 400 T0 from models with different dust loads. The left and right columns correspond to models with α = 10−4and 10−5, respectively. Upper panels: Results from the fiducial metallicity Z = 0.1. Lower panels: Models with Z = 0.01. Dashed curves denote contours for 25% of the gas scale height.
Thumbnail: Fig. 8 Refer to the following caption and surrounding text. Fig. 8

Dust-to-gas density ratio at t = 400 T0 from models with different reference Stokes numbers. The left and right columns correspond to models with α = 10−4 and 10−5, respectively. Upper panels: Results from the fiducial reference Stokes number St0 = 10−3. Lower panels: Models with St0 = 10−2. Dashed curves denote contours for 25% of the gas scale height, and arrows below denote the Stokes number of dust grains in the midplane. The spatial domains differ between the upper and lower panels due to the different radial drift timescales.
Thumbnail: Fig. 9 Refer to the following caption and surrounding text. Fig. 9

Radial profiles of the dust-to-gas scale height ratio (left y-axis, blue curves) and the maximum dust-to-gas density ratio (right y-axis, orange curves) at t = 400 T0 from models with α = 10−4 (panel a) and α = 10−5 (panel b), both with St0 = 10−2. The dust scale heights are measured using the standard deviation (HdstdMathematical equation: $H_{\rm d}^{\rm std}$, Eq. (28)), and the dashed curve indicates a reference-only expected ratio assuming explicit dust diffusion, as given by Eq. (29). The radial extent is the same as in Figs. 1 and 4.

4.2 Outlooks

Future work should adopt a more vertically extended domain to better capture VSI development and more rigorously assess its interaction with SI in viscous environments. Extending the models to three dimensions will also be essential to explore non-axisymmetric structures, such as zonal flows, vortices, and possible secondary instabilities arising from SI-VSI interactions. In addition, implementing numerical methods that consistently account for both dust diffusion and dust back-reaction will

enable a more controlled investigation of the effects of diffusion on dust substructures.

We will also assess the feasibility of observing structures generated by VSI and SI and quantify their expected observational signatures. The potential observations will rely on dust grain alignment and the resulting net linear polarization of the thermal re-emission radiation, explicitly accounting for competing processes, in particular self-scattering (Lietzow-Sinjen et al. 2025). Finally, we will derive quantitative observational requirements and compare them with the specifications of current and next-generation submillimeter and millimeter interferometric facilities.

5 Conclusion

In this work, we investigated the radial and vertical distribution of dust in protoplanetary disks by performing twodimensional, axisymmetric gas-dust simulations that include dust back-reaction and parameterized viscosity. Motivated by ALMA observations revealing ubiquitous dust substructures, we focused on understanding how the interplay between SI, VSI, vertical settling, and radial drift shapes the morphology and concentration of marginally coupled dust in disks. The results of our study can be summarized as follows:

  • In models with moderate viscosity α = 10−4, we identify a characteristic shuttlecock-shaped dust substructure composed of a dense, vertically settled “head” and a vertically extended, diffuse “tail” (Fig. 1). This morphology originates from the SI (Fig. 3), triggered by the back-reaction of marginally coupled dust concentrated near the midplane;

  • For dust grains with 5 × 10−3 ≲ St ≲ 10−2, the dust-togas scale height ratio reaches ≳0.25 in the shuttlecock tails (Fig. 2a), while the dust-to-gas density ratio attains ≲3 in the heads (Fig. 1). For larger dust grains with 1.5 × 10−2 ≲ St ≲ 2.5 × 10−2, the scale height ratio reaches ≳0.1 and the density ratio increases up to ≲10 (Fig. 9a). In both regimes, the maximum dust scale height exceeds the reference value predicted by the balance between turbulent diffusion and vertical settling (Eq. (29)), indicating strong vertical stirring induced by the SI. However, the maximum dust densities remain below the critical Hill density required for planetesimal formation (Eq. (34));

  • In lower-viscosity models (α = 10−5), the shuttlecock morphology disappears, as the VSI becomes sufficiently vigorous to disrupt dust structures near the midplane (Fig. 4). We find that even a midplane dust-to-gas density ratio of order unity is insufficient to suppress VSI when a small but finite background viscosity is present;

  • The enhanced instabilities in lower-viscosity models lead to dust scale heights much larger than those predicted by Eq. (29). For grains with St ≳ 10−2, the dust-to-gas density ratio can reach up to 20-50 (Fig. 9b), potentially exceeding the Hill density and thus approaching conditions favorable for gravitationally bound clump formation.

Acknowledgements

We thank the anonymous reviewer for the helpful comments. We also thank Pinghui Huang for insightful discussions on disk instabilities and dust dynamics. This work was supported by the Deutsche Forschungsgemeinschaft (DFG) under Grant No. 544937803. J.B. further gratefully acknowledges funding from COST Action CA22133 PLANETS, supported by COST (European Cooperation in Science and Technology). The simulations were carried out on computing clusters hosted by the Max Planck Computing and Data Facility (MPCDF). The data underlying this article will be made available upon reasonable request to the corresponding author.

References

  1. ALMA Partnership (Brogan, C. L., et al.) 2015, ApJ, 808, L3 [Google Scholar]
  2. Andrews, S. M., Wilner, D. J., Zhu, Z., et al. 2016, ApJ, 820, L40 [Google Scholar]
  3. Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJ, 869, L41 [NASA ADS] [CrossRef] [Google Scholar]
  4. Bai, X.-N. 2013, ApJ, 772, 96 [Google Scholar]
  5. Bai, X.-N., & Stone, J. M. 2010, ApJ, 722, 1437 [NASA ADS] [CrossRef] [Google Scholar]
  6. Bai, X.-N., & Stone, J. M. 2013, ApJ, 769, 76 [Google Scholar]
  7. Bai, X.-N., Ye, J., Goodman, J., & Yuan, F. 2016, ApJ, 818, 152 [Google Scholar]
  8. Barker, A. J., & Latter, H. N. 2015, MNRAS, 450, 21 [NASA ADS] [CrossRef] [Google Scholar]
  9. Benisty, M., Bae, J., Facchini, S., et al. 2021, ApJ, 916, L2 [NASA ADS] [CrossRef] [Google Scholar]
  10. Benítez-Llambay, P., & Masset, F. S. 2016, ApJS, 223, 11 [Google Scholar]
  11. Benítez-Llambay, P., Krapp, L., & Pessah, M. E. 2019, ApJS, 241, 25 [Google Scholar]
  12. Béthune, W., Lesur, G., & Ferreira, J. 2017, A&A, 600, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. Birnstiel, T. 2024, ARA&A, 62, 157 [Google Scholar]
  14. Birnstiel, T., Dullemond, C. P., & Brauer, F. 2010, A&A, 513, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Blum, J., & Wurm, G. 2008, ARA&A, 46, 21 [CrossRef] [Google Scholar]
  16. Booth, A. S., & Ilee, J. D. 2020, MNRAS, 493, L108 [NASA ADS] [CrossRef] [Google Scholar]
  17. Bosschaart, Q., Guerra-Alvarado, O. M., Van Der Marel, N., & Mulders, G. D. 2026, A&A, 708, A143 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Brauer, F., Dullemond, C. P., & Henning, T. 2008, A&A, 480, 859 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  19. Carrera, D., Johansen, A., & Davies, M. B. 2015, A&A, 579, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  20. Casassus, S., Marino, S., Lyra, W., et al. 2019, MNRAS, 483, 3278 [Google Scholar]
  21. Chen, K., & Lin, M.-K. 2020, ApJ, 891, 132 [Google Scholar]
  22. Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368 [Google Scholar]
  23. Cieza, L. A., Ruíz-Rodríguez, D., Hales, A., et al. 2019, MNRAS, 482, 698 [Google Scholar]
  24. Dong, R., Liu, S.-y., Eisner, J., et al. 2018, ApJ, 860, 124 [Google Scholar]
  25. Drazkowska, J., Bitsch, B., Lambrechts, M., et al. 2023, in Protostars Planets VII, 534, 717 [Google Scholar]
  26. Dubrulle, B., Morfill, G., & Sterzik, M. 1995, Icarus, 114, 237 [Google Scholar]
  27. Dullemond, C. P., Birnstiel, T., Huang, J., et al. 2018, ApJ, 869, L46 [NASA ADS] [CrossRef] [Google Scholar]
  28. Facchini, S., Curone, P., Benisty, M., et al. 2026, ApJ, 998, L16 [Google Scholar]
  29. Francis, L., & van der Marel, N. 2020, ApJ, 892, 111 [Google Scholar]
  30. Fukagawa, M., Tsukagoshi, T., Momose, M., et al. 2013, PASJ, 65, L14 [Google Scholar]
  31. Güttler, C., Blum, J., Zsom, A., Ormel, C. W., & Dullemond, C. P. 2010, A&A, 513, A56 [Google Scholar]
  32. Harten, A., Lax, P. D., & van Leer, B. 1983, SIAM Rev., 25, 35 [Google Scholar]
  33. Huang, P., & Bai, X.-N. 2022, ApJS, 262, 11 [NASA ADS] [CrossRef] [Google Scholar]
  34. Huang, P., & Bai, X.-N. 2025a, ApJ, 986, L13 [Google Scholar]
  35. Huang, P., & Bai, X.-N. 2025b, ApJ, 986, 76 [Google Scholar]
  36. Huang, J., Andrews, S. M., Pérez, L. M., et al. 2018, ApJ, 869, L43 [Google Scholar]
  37. Isella, A., Guidi, G., Testi, L., et al. 2016, Phys. Rev. Lett., 117, 251101 [NASA ADS] [CrossRef] [Google Scholar]
  38. Isella, A., Huang, J., Andrews, S. M., et al. 2018, ApJ, 869, L49 [NASA ADS] [CrossRef] [Google Scholar]
  39. Johansen, A., & Youdin, A. 2007, ApJ, 662, 627 [Google Scholar]
  40. Kanagawa, K. D., Ueda, T., Muto, T., & Okuzumi, S. 2017, ApJ, 844, 142 [NASA ADS] [CrossRef] [Google Scholar]
  41. Klahr, H., & Schreiber, A. 2020, ApJ, 901, 54 [NASA ADS] [CrossRef] [Google Scholar]
  42. Kurtovic, N. T., Pérez, L. M., Benisty, M., et al. 2018, ApJ, 869, L44 [Google Scholar]
  43. Lesur, G., Flock, M., Ercolano, B., et al. 2023, in Protostars Planets VII, 534, 465 [Google Scholar]
  44. Lietzow-Sinjen, M., Reissl, S., Flock, M., & Wolf, S. 2025, A&A, 703, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  45. Lim, J., Simon, J. B., Li, R., et al. 2024, ApJ, 969, 130 [NASA ADS] [CrossRef] [Google Scholar]
  46. Lim, J., Simon, J. B., Li, R., et al. 2025, ApJ, 981, 160 [Google Scholar]
  47. Lin, M.-K. 2019, MNRAS, 485, 5221 [CrossRef] [Google Scholar]
  48. Lin, M.-K., & Youdin, A. N. 2015, ApJ, 811, 17 [NASA ADS] [CrossRef] [Google Scholar]
  49. Lin, M.-K., & Youdin, A. N. 2017, ApJ, 849, 129 [NASA ADS] [CrossRef] [Google Scholar]
  50. Long, F., Pinilla, P., Herczeg, G. J., et al. 2018, ApJ, 869, 17 [Google Scholar]
  51. Long, F., Herczeg, G. J., Harsono, D., et al. 2019, ApJ, 882, 49 [Google Scholar]
  52. Long, F., Andrews, S. M., Zhang, S., et al. 2022, ApJ, 937, L1 [NASA ADS] [CrossRef] [Google Scholar]
  53. Lyra, W., & Lin, M.-K. 2013, ApJ, 775, 17 [Google Scholar]
  54. Martire, P., Longarini, C., Lodato, G., et al. 2024, A&A, 686, A9 [Google Scholar]
  55. Melon Fuksman, J. D., Flock, M., & Klahr, H. 2024, A&A, 682, A140 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  56. Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228 [Google Scholar]
  57. Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, MNRAS, 435, 2610 [Google Scholar]
  58. Ostertag, D., & Flock, M. 2025, A&A, 695, L13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  59. Pérez, L. M., Carpenter, J. M., Andrews, S. M., et al. 2016, Science, 353, 1519 [Google Scholar]
  60. Pérez, S., Casassus, S., Baruteau, C., et al. 2019, AJ, 158, 15 [Google Scholar]
  61. Pérez, S., Casassus, S., Hales, A., et al. 2020, ApJ, 889, L24 [Google Scholar]
  62. Pinilla, P., Birnstiel, T., Ricci, L., et al. 2012, A&A, 538, A114 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  63. Ribas, Á., Clarke, C. J., & Zagaria, F. 2024, MNRAS, 532, 1752 [NASA ADS] [CrossRef] [Google Scholar]
  64. Roe, P. L. 1981, J. Comput. Phys., 43, 357 [Google Scholar]
  65. Rosotti, G. P. 2023, New Astron. Rev., 96, 101674 [NASA ADS] [CrossRef] [Google Scholar]
  66. Rosotti, G. P., Teague, R., Dullemond, C., Booth, R. A., & Clarke, C. J. 2020, MNRAS, 495, 173 [Google Scholar]
  67. Schäfer, U., & Johansen, A. 2022, A&A, 666, A98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  68. Schäfer, U., Johansen, A., & Banerjee, R. 2020, A&A, 635, A190 [Google Scholar]
  69. Schäfer, U., Johansen, A., & Flock, M. 2025, A&A, 694, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  70. Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
  71. Squire, J., & Hopkins, P. F. 2018a, MNRAS, 477, 5011 [Google Scholar]
  72. Squire, J., & Hopkins, P. F. 2018b, ApJ, 856, L15 [NASA ADS] [CrossRef] [Google Scholar]
  73. Stoll, M. H. R., & Kley, W. 2014, A&A, 572, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  74. Stoll, M. H. R., & Kley, W. 2016, A&A, 594, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  75. Toro, E. F. 2009, Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction, 3rd edn. (Heidelberg: Springer Berlin) [Google Scholar]
  76. Umurhan, O. M., Estrada, P. R., & Cuzzi, J. N. 2020, ApJ, 895, 4 [Google Scholar]
  77. van der Marel, N., van Dishoeck, E. F., Bruderer, S., et al. 2013, Science, 340, 1199 [Google Scholar]
  78. van der Marel, N., Williams, J. P., Ansdell, M., et al. 2018, ApJ, 854, 177 [Google Scholar]
  79. van der Marel, N., Birnstiel, T., Garufi, A., et al. 2021, AJ, 161, 33 [Google Scholar]
  80. van der Marel, N., Williams, J. P., Picogna, G., et al. 2022 [arXiv:2204.08225] [Google Scholar]
  81. Weber, P., Benítez-Llambay, P., Gressel, O., Krapp, L., & Pessah, M. E. 2018, ApJ, 854, 153 [NASA ADS] [CrossRef] [Google Scholar]
  82. Weidenschilling, S. J. 1977, MNRAS, 180, 57 [Google Scholar]
  83. Wölfer, L., Barraza-Alfaro, M., Teague, R., et al. 2025, ApJ, 984, L22 [Google Scholar]
  84. Wu, Y., Yu, C., & Cui, C. 2024, MNRAS, 534, 948 [Google Scholar]
  85. Yang, C.-C., Johansen, A., & Carrera, D. 2017, A&A, 606, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  86. Yoshida, T. C., Nomura, H., Doi, K., et al. 2025, Nat. Astron, 1 [Google Scholar]
  87. Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459 [Google Scholar]
  88. Youdin, A. N., & Lithwick, Y. 2007, Icarus, 192, 588 [NASA ADS] [CrossRef] [Google Scholar]
  89. Ziampras, A., Sudarshan, P., Dullemond, C. P., et al. 2025, MNRAS, 536, 3322 [NASA ADS] [CrossRef] [Google Scholar]
  90. Zsom, A., Ormel, C. W., Güttler, C., Blum, J., & Dullemond, C. P. 2010, A&A, 513, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  91. Zsom, A., Ormel, C. W., Dullemond, C. P., & Henning, T. 2011, A&A, 534, A73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
1

Since the global ρg profile does not significantly deviate from its initial condition and cs is time-independent, we hereafter use St0 to refer to St0iniMathematical equation: ${\rm St}_0^{\rm ini}$ for brevity.

2

Assuming a disk mass of 0.2 M⊙ and an outer disk radius of 150 au (Booth & Ilee 2020).

3

A second motivation is that the dust diffusion prescription in the FARGO3D code cannot account for dust back-reaction and therefore does not strictly conserve the angular momentum of the gas (Weber et al. 2018), complicating a consistent treatment of dust feedback in turbulent flows.

4

The reference orbital period is given by T0 = 2π/ΩK0.

5

Lim et al. (2024) did not include explicit dust diffusion, but they accounted for dust self-gravity. Consequently, the threshold for strong clumping in our setup may be even higher.

Appendix A Numerical resolution test

To assess the robustness of our results against numerical resolution, we conducted a series of resolution tests in which the grid resolution was doubled separately in the radial and vertical directions for both the α = 10−4 and α = 10−5 models. The results are shown in Fig. A.1. Owing to the computational cost of higher-resolution simulations, these models were evolved for a shorter duration, t = 100 T0, which is sufficient to capture the initial development of the dust substructures. For the α = 10−4 models, we find good agreement in the overall dust morphology across all resolutions. In the α = 10−5 models, the higher-resolution runs, particularly the one with doubled vertical resolution, exhibit more pronounced midplane corrugations due to better-resolved VSI activity. Nevertheless, the fiducial simulation reproduces the key features of the dust distribution. These tests indicate that our results are numerically converged at the fiducial resolution and are not artifacts of grid discretization in either the radial or vertical direction.

Thumbnail: Fig. A.1 Refer to the following caption and surrounding text. Fig. A.1

Dust-to-gas density ratio at t = 100 T0 from models with different numerical resolutions. The left and right columns correspond to models with α = 10−4 and 10−5, respectively. Middle row: Models with the fiducial resolution (i.e., ~200 cells per Hg in both directions). Top and bottom rows: Models with two times higher resolution than fiducial in the r direction (top) and θ direction (bottom).

Appendix B VSI activity in models with different α

Thumbnail: Fig. B.1 Refer to the following caption and surrounding text. Fig. B.1

Specific angular momentum of gas at t = 400 T0 from models with different vertical prescriptions of the α parameter. Panels a, d, and g: Models in which α = 10−3, 10−4, and 10−5, respectively. Panel b: Model in which α decreases from 10−3 in the midplane to 10−4 in the atmosphere. Panel c: Model in which α increases from 10−4 in the midplane to 10−3 in the atmosphere. Panel e: Model in which α decreases from 10−4 in the midplane to 10−5 in the atmosphere. Panel f : Model in which α increases from 10−5 in the midplane to 10−4 in the atmosphere. Last panel: Corresponding vertical profiles of α.

To examine the activity of VSI in our models and its role in shaping the dust layer, we performed two additional simulations in which α varies vertically between 10−4 and 10−5, similar to the setup explored in Sect. 3.1. The resulting gas specific angular momentum is shown in Fig. B.1, and the corresponding dust-to-gas density ratio is presented in Fig. B.2. We observe more pronounced signatures of VSI activity only in panels e and g of Fig. B.1, which correspond to models with α = 10−5 in the disk atmosphere. This behavior is consistent with our earlier conclusion that the shuttlecock-shaped substructure is primarily driven by SI rather than VSI. Meanwhile, the comparison between panels d and e of Fig. B.2 shows that the midplane dust layer becomes more corrugated when α is reduced in the disk atmosphere, consistent with stronger VSI activity. The comparison between panels f and g further shows that when α decreases in the disk midplane, the dust layer can form more clearly defined shuttlecock-shaped substructures (panel f ), although the morphology can also be significantly modulated by VSI-induced corrugations (panel g).

Thumbnail: Fig. B.2 Refer to the following caption and surrounding text. Fig. B.2

Similar to Fig. B.1 but for the dust-to-gas density ratio. Dashed curves in panels b, c, e, and f denote contours for one gas scale height. The shaded regions indicate where the transition of α occurs and are bounded by aspect ratios evaluated at the inner and outer radial boundaries.

Appendix C Dust substructures in VSI-suppressed models

To further confirm that the shuttlecock-shaped dust substructures are driven primarily by the SI rather than the VSI, we performed an additional “toy-model” simulation at half the fiducial resolution in both dimensions, in which the temperature profile was modified to suppress VSI activity. Specifically, we imposed a radially constant temperature profile by setting the disk flaring index to 0.5. This removes the vertical shear in the gas rotation profile and thereby suppresses the VSI. The resulting dust-to-gas density ratio is shown in Fig. C.1. We find that the shuttlecock-shaped substructure still forms in this VSI-suppressed model, demonstrating that the SI alone can generate such dust morphologies in the absence of the VSI. This result further supports our conclusion that the dust substructures observed in our simulations are driven primarily by the SI rather than the VSI.

Thumbnail: Fig. C.1 Refer to the following caption and surrounding text. Fig. C.1

Dust-to-gas density ratio at t = 200 T0 from models with our fiducial temperature profile T ∝ r3/7 (panel a) and with a radially isothermal prescription (panel b). VSI is therefore completely suppressed in the latter case, and the resulting dust substructure is solely shaped by SI.

Appendix D Cross-code validation test

To further verify the robustness of the results from the α = 10−4 model, we performed cross-code comparisons using the PLUTO code (Mignone et al. 2007) and its dust-fluid module (Ziampras et al. 2025) with the same physical setup described in Sect. 2. The results of this comparison are shown in Fig. D.1.

To explore the potential impact of numerical diffusion, we carried out two PLUTO simulations with different numerical configurations. In the setup designed to minimize numerical diffusion, we adopted fourth-order parabolic reconstruction, third-order Runge-Kutta time integration for the gas, first-order implicit time integration for the dust, and the Roe Riemann solver (Roe 1981; Toro 2009). In the more diffusive configuration, we instead employed the HLL Riemann solver (Harten et al. 1983), replaced the parabolic reconstruction with linear, second-order reconstruction, and used second-order Runge-Kutta time integration for the gas.

Despite these differences, we find good agreement among the three models. In particular, the formation of the shuttlecockshaped substructures and the magnitude of dust concentration are consistently reproduced. Minor quantitative differences appear in the detailed morphology and small-scale structures, likely due to differences in numerical methods. However, these variations do not affect the qualitative behavior or the main conclusions. This comparison demonstrates that our results are not specific to a particular numerical implementation and are robust across different simulation frameworks.

Thumbnail: Fig. D.1 Refer to the following caption and surrounding text. Fig. D.1

Dust-to-gas density ratio at t = 200 T0 from models evolved with FARGO3D (panel a), the less diffusive Roe solver in PLUTO (panel b), and the more diffusive HLL solver in PLUTO (panel c).

All Figures

Thumbnail: Fig. 1 Refer to the following caption and surrounding text. Fig. 1

Dust-to-gas density ratio at t = 400 T0 from the model with α = 10−4 and St0 = 10−3. Dashed curves in the upper panel denote contours for 25% of the gas scale height, and the arrow below denotes the Stokes number of dust grains in the midplane. Gas streamlines are overplotted in the lower panel, with the line color denoting the vertical gas velocity.
In the text
Thumbnail: Fig. 2 Refer to the following caption and surrounding text. Fig. 2

Radial profiles of the dust-to-gas scale height ratio at t = 400 T0 from models with α = 10−4 (panel a) and α = 10−5 (panel b), both with St0 = 10−3. The dust scale heights are measured using the root mean square (HdrmsMathematical equation: $H_{\rm d}^{\rm rms}$, Eq. (27)), the standard deviation (HdstdMathematical equation: $H_{\rm d}^{\rm std}$, Eq. (28)), or by fitting a Gaussian profile. The dashed curve indicates a reference-only expected ratio assuming explicit dust diffusion, as given by Eq. (29). The radial extent is the same as in Figs. 1 and 4.
In the text
Thumbnail: Fig. 3 Refer to the following caption and surrounding text. Fig. 3

Dust-to-gas density ratio at t = 400 T0 from models with different vertical prescriptions of the α parameter. Panel a : model in which α = 10−3. Panel b : model in which α decreases from 10−3 in the midplane to 10−4 in the atmosphere. Panel c: model in which α increases from 10−4 in the midplane to 10−3 in the atmosphere. Panel d: model in which α = 10−4. Last panel: Corresponding vertical profiles of α, given by Eqs. (30)-(33). Dashed curves in panels b and c denote contours for one gas scale height. The shaded regions indicate where the transition of α occurs and are bounded by aspect ratios evaluated at the inner and outer radial boundaries.
In the text
Thumbnail: Fig. 4 Refer to the following caption and surrounding text. Fig. 4

Similar to Fig. 1, but from the model with α = 10−5. The colorbar for gas streamlines has been adjusted to fit stronger vertical gas flows in this lower-viscosity disk model.
In the text
Thumbnail: Fig. 5 Refer to the following caption and surrounding text. Fig. 5

Specific angular momentum of gas at t = 400 T0 from models with α = 10−4 (panel a) and α = 10−5 (panel b), both with St0 = 10−3. The radial extent is the same as in Figs. 1 and 4.
In the text
Thumbnail: Fig. 6 Refer to the following caption and surrounding text. Fig. 6

Dust-to-gas density ratio at t = 400 T0 from models with different α values. Dashed curves denote contours for 25% of the gas scale height.
In the text
Thumbnail: Fig. 7 Refer to the following caption and surrounding text. Fig. 7

Dust-to-gas density ratio at t = 400 T0 from models with different dust loads. The left and right columns correspond to models with α = 10−4and 10−5, respectively. Upper panels: Results from the fiducial metallicity Z = 0.1. Lower panels: Models with Z = 0.01. Dashed curves denote contours for 25% of the gas scale height.
In the text
Thumbnail: Fig. 8 Refer to the following caption and surrounding text. Fig. 8

Dust-to-gas density ratio at t = 400 T0 from models with different reference Stokes numbers. The left and right columns correspond to models with α = 10−4 and 10−5, respectively. Upper panels: Results from the fiducial reference Stokes number St0 = 10−3. Lower panels: Models with St0 = 10−2. Dashed curves denote contours for 25% of the gas scale height, and arrows below denote the Stokes number of dust grains in the midplane. The spatial domains differ between the upper and lower panels due to the different radial drift timescales.
In the text
Thumbnail: Fig. 9 Refer to the following caption and surrounding text. Fig. 9

Radial profiles of the dust-to-gas scale height ratio (left y-axis, blue curves) and the maximum dust-to-gas density ratio (right y-axis, orange curves) at t = 400 T0 from models with α = 10−4 (panel a) and α = 10−5 (panel b), both with St0 = 10−2. The dust scale heights are measured using the standard deviation (HdstdMathematical equation: $H_{\rm d}^{\rm std}$, Eq. (28)), and the dashed curve indicates a reference-only expected ratio assuming explicit dust diffusion, as given by Eq. (29). The radial extent is the same as in Figs. 1 and 4.
In the text
Thumbnail: Fig. A.1 Refer to the following caption and surrounding text. Fig. A.1

Dust-to-gas density ratio at t = 100 T0 from models with different numerical resolutions. The left and right columns correspond to models with α = 10−4 and 10−5, respectively. Middle row: Models with the fiducial resolution (i.e., ~200 cells per Hg in both directions). Top and bottom rows: Models with two times higher resolution than fiducial in the r direction (top) and θ direction (bottom).
In the text
Thumbnail: Fig. B.1 Refer to the following caption and surrounding text. Fig. B.1

Specific angular momentum of gas at t = 400 T0 from models with different vertical prescriptions of the α parameter. Panels a, d, and g: Models in which α = 10−3, 10−4, and 10−5, respectively. Panel b: Model in which α decreases from 10−3 in the midplane to 10−4 in the atmosphere. Panel c: Model in which α increases from 10−4 in the midplane to 10−3 in the atmosphere. Panel e: Model in which α decreases from 10−4 in the midplane to 10−5 in the atmosphere. Panel f : Model in which α increases from 10−5 in the midplane to 10−4 in the atmosphere. Last panel: Corresponding vertical profiles of α.
In the text
Thumbnail: Fig. B.2 Refer to the following caption and surrounding text. Fig. B.2

Similar to Fig. B.1 but for the dust-to-gas density ratio. Dashed curves in panels b, c, e, and f denote contours for one gas scale height. The shaded regions indicate where the transition of α occurs and are bounded by aspect ratios evaluated at the inner and outer radial boundaries.
In the text
Thumbnail: Fig. C.1 Refer to the following caption and surrounding text. Fig. C.1

Dust-to-gas density ratio at t = 200 T0 from models with our fiducial temperature profile T ∝ r3/7 (panel a) and with a radially isothermal prescription (panel b). VSI is therefore completely suppressed in the latter case, and the resulting dust substructure is solely shaped by SI.
In the text
Thumbnail: Fig. D.1 Refer to the following caption and surrounding text. Fig. D.1

Dust-to-gas density ratio at t = 200 T0 from models evolved with FARGO3D (panel a), the less diffusive Roe solver in PLUTO (panel b), and the more diffusive HLL solver in PLUTO (panel c).
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.