© The Authors 2026

1. Introduction

In recent years, high-resolution observations of protoplanetary disks by the Atacama Large Millimeter/Submillimeter Array (ALMA) have revealed striking dust substructures, in the form of rings, gaps, spirals, and asymmetric crescent-shaped features (Bae et al. 2023). Observations in the submillimeter (sub-mm) dust continuum of dust crescents suggest that large dust particles are being trapped in the disk midplane (van der Marel et al. 2013; Pérez et al. 2018). Determining which mechanism is primarily responsible for the observed dust asymmetries remains a significant challenge.

The most widely accepted physical mechanism for generating gas vortices is the Rossby wave instability (RWI; Lovelace et al. 1999; Li et al. 2000, 2001). This instability is triggered by a local minimum in the radial profile of the potential vorticity (or vortensity), defined as $\eta =\frac{{\kappa }^2{S}^{-2/\gamma }}{2{\mathrm{\Omega }}_k{\mathrm{\Sigma }}_{\mathrm{gas}}}$, where κ is the epicyclic frequency, which describes the oscillation frequency of a fluid element undergoing small radial displacements about a circular orbit and is given by $\kappa =R^{-3}\frac{\partial (R^4{\mathrm{\Omega }}_k^2)}{\partial R}$; Ω_k is the Keplerian angular velocity; Σ_gas is the gas surface density; S = P/Σ_gas^γ denotes the entropy, with P as the gas pressure; and γ is the adiabatic index. This quantity depends on the disk’s density, temperature, and rotational structure and is governed by the balance of radial forces including the centrifugal, Coriolis, and pressure gradient forces. The RWI leads to the formation of anticyclonic vortices at planetary gap edges, provided the planet is sufficiently massive–comparable to a Jupiter- or Neptune-mass planet (de Val-Borro et al. 2007; Lyra et al. 2009a). Vortices may also arise near disk boundaries where sharp pressure or viscosity transitions generate long-lived pressure maxima, often observed as dusty rings in ALMA observations (Lin 2014). Additionally, RWI-induced vortices can form in disks perturbed by a binary companion on an eccentric orbit (i.e., Price et al. 2018).

The RWI can also emerge in turbulent disk environments driven by the vertical shear instability (VSI; Nelson et al. 2013). The VSI is a hydrodynamical (HD) instability in accretion disks that arises from vertical gradients in rotational velocity, leading to velocity perturbations and turbulence in both the vertical and radial directions (Barker & Latter 2015; Stoll & Kley 2016; Lyra & Umurhan 2019). Richard et al. (2016) first demonstrated that small vortices can form in disks susceptible to the VSI. Later, Latter & Papaloizou (2018) confirmed that the Kelvin-Helmholtz (KH) instability, acting as a parasitic instability in VSI modes, is responsible for the small-scale vortices. Richard et al. (2016) also suggested that secondary instabilities, such as the RWI or the elliptical instability (E.I.; Lesur & Papaloizou 2009), may ultimately determine the fate of these vortices, an idea later summarized by Lyra & Umurhan (2019). Building on subsequent vortex studies in VSI scenarios, Manger & Klahr (2018) showed that the VSI can generate large-scale, long-lived vortices – azimuthally elongated and vertically aligned – as a secondary instability triggered by the RWI (see also Manger et al. 2020; Pfeil & Klahr 2021). While such vortices have been commonly reported in VSI simulations under locally isothermal configuration, recent high-resolution simulations by Shariff & Umurhan (2024), with up to 70 cells per scale height, were the first to show that large-scale vortices are likely absent due to the enhanced dynamical range of kinetic energy cascading to smaller scales. This resolution effect was later confirmed by Lesur et al. (2025), who also reported the absence of long-lived midplane vortices and attributed this to the action of the E.I. These results highlight the ongoing uncertainties regarding how resolution and gas-cooling timescales influence vortex formation in VSI-driven turbulence.

It remains unclear which mechanism initiates the observed dust asymmetries – whether the RWI acts as a secondary instability in VSI-dominated disks, or a massive planet first carves a gap that subsequently triggers a vortex responsible for the dust crescent. The HD simulations by Huang et al. (2019) have shown that spiral arms generated by RWI-driven vortices exhibit surface density contrasts comparable to those produced by subthermal mass planets but are too weak to be detected in scattered light, suggesting that the prominent near-infrared (NIR) spirals observed in systems such as SAO 206462 and MWC 758 are more likely planet-induced. To distinguish between these mechanisms and their observational manifestations, it is essential to analyze the associated gas dynamics.

Vortices have been widely investigated as efficient dust traps in protoplanetary disks, where local pressure maxima can concentrate solids and potentially trigger planetesimal formation (Johansen et al. 2004; Lyra et al. 2009a). In particular, RWI vortices have been studied in the context of planet–disk interaction simulations that include Lagrangian dust particles, and their ability to trap dust until gravitational collapse leads to the formation of bodies with masses up to that of a super-Earth has been confirmed (Lyra et al. 2009a; Fu et al. 2014; Lin 2014; Bae et al. 2016; Ma et al. 2025). Vortices formed at the edges of planet-induced gaps can become significantly elongated in azimuth when the planet grows over sufficiently long timescales, a feature that may help distinguish them observationally from other asymmetry-generating mechanisms (Hammer et al. 2017, 2019). The lifetime of such vortices can extend to several thousand orbits when dust coagulation and fragmentation are taken into account (Fu et al. 2014). While previous studies have explored VSI vortices in gas-only simulations (Richard et al. 2016; Manger & Klahr 2018; Manger et al. 2021; Pfeil & Klahr 2021), only a few works (Flock et al. 2020; Blanco et al. 2021; Huang & Bai 2025) have conducted 3D simulations demonstrating that VSI-driven turbulence can concentrate dust into RWI vortices and small-scale clumps that are related to the KH instability.

From an observational standpoint, the recent exoALMA Large Program (Teague et al. 2025) has presented new high-resolution gas observations including four of the disks that exhibit prominent dust crescents in the submillimeter continuum. These observations achieve a high angular resolution of 100 mas (approximately 14 au at the typical distances of the sources) and a spectral resolution of 26 m s⁻¹, enabling the detection of subtle structures and motions that offer critical insights into dust asymmetries and their role as physical processes governing planet formation. In Wölfer et al. (2025), the authors analyzed the line-of-sight velocities of ¹²CO (J = 3 − 2) and ¹³CO (J = 3 − 2) and identified deviations from purely Keplerian rotation around the four disks with dust crescents. Vortex diagnostics in simulations reveal distinct azimuthally asymmetric red-blue velocity patterns which have been used to probe deviations from Keplerian motion (e.g., Robert et al. 2020). In their HD simulation work, vortices formed at the edge of gas-depleted regions produced apparent spiral arm features in the gas line-of-sight velocity maps. However, the interpretation of such features in observations remains challenging due to the complexity of gas dynamics, where multiple processes (e.g., other disk instabilities or binary interactions) may obscure or mimic vortex signatures.

In this paper, we investigate the influence of VSI-driven, self-sustained vortices on dust distribution and gas kinematics and their potential observational signatures in protoplanetary disks. In Section 2, we describe the 3D simulation setup for both gas and dust. Section 4 presents the results and we focus on the locations within the disk where particles become concentrated and trapped inside the vortices (see Appendix A for the RWI vortex criterion). We used a single time snapshot to test whether large-scale dust concentrations at similar locations motivated by those observed in the MWC 758 disk (Dong et al. 2018) can arise from vortices triggered by the RWI in a turbulent disk driven by the VSI. However, we are not specifically concerned with matching the models solely to the properties of MWC 758. In general, we aim to link observed submillimeter asymmetries to vortex-induced dust trapping from our models. We demonstrate how post-processing the simulation data allowed us to diagnose dust asymmetries around vortices, analyze gas kinematic perturbations from Keplerian motion, and predict how these features would appear in ALMA observations. In Section 6, we evaluate the feasibility of explaining observed dust asymmetries in disks based on our model predictions. Finally, our conclusions are summarized in Section 7.

2. Methods

We performed 3D simulations of the gas in the protoplanetary disk based on the principles of mass, momentum, and energy conservation using the PLUTO¹ code (Mignone 2007). The simulations are carried in spherical coordinates (r, θ, ϕ) using a uniform spaced grid in θ and ϕ, and logarithmic spaced grid in r. The fluxes in each cell are computed using the parabolic reconstruction method using fourth order with the Harten, Lax, Van Leer (HLLC) solver. For the time integration we use the third-order Runge-Kutta algorithm with the Courant-Friedrichs-Lewy (CFL) number set to 0.3. The grid extends from 0.5 R₀ to 2 R₀ in the radial direction, where R₀ = 1 in code units, and from π/2 − 0.35 to π/2 + 0.35 (i.e., ±3.5H) in the meridional direction. The analysis are presented in cylindrical coordinates (R, ϕ, Z) and are scaled to 40 au in physical units. This numerical setup was chosen to achieve a resolution of approximately 37 cells per scale height in all directions (Table 1). Our disk configuration follows the equilibrium solutions of the density and rotational velocity profile in cylindrical coordinates based on Nelson et al. (2013)

$$\begin{array}{cc}\hfill \rho (R,Z)& ={\rho }_0{\left(\frac{R}{R_0}\right)}^{p}exp\left(\frac{-Z^2}{2H^2}\right),\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill \mathrm{\Omega }(R,Z)& ={\mathrm{\Omega }}_{\mathrm{k}}{[(p+q){\left(\frac{H}{R}\right)}^2+(1+q)-\frac{\mathit{qR}}{\sqrt{R^2+Z^2}}]}^{1/2}.\hfill \end{array}$$(2)

The variables p and q corresponds to the exponents in the radial profiles of temperature and midplane density, respectively. The initial midplane density is given by ${\rho }_0=\frac{{\mathrm{\Sigma }}_0}{\sqrt{2\pi }H/R{R}_0}$, where Σ₀ = 30 g cm⁻² is the initial surface density at the reference radius R₀, and H/R is the aspect ratio (see Table 1). The scale height profile of the gas is H ∝ R^{(q + 3)/2}.

The gas is modeled as an ideal gas, with internal energy e related to temperature and density via the caloric equation of state. The pressure is given by P = ρe/(γ − 1), where γ = 1.4 is the adiabatic index appropriate for diatomic molecular hydrogen. The gas pressure is defined as $P=c_{\mathrm{s}}^2\rho$, where $c_{\mathrm{s}}\propto {\left(\frac{R}{R_0}\right)}^q$ is the radially varying isothermal sound speed. The temperature profile is related to the gas pressure based on c_s as T ∝ c_s² and the sound speed is related to the gas pressure scale height as c_s = HΩ_k.

We used a simple thermal relaxation or cooling relaxation scheme in which the gas pressure evolves toward an equilibrium state over a characteristic cooling timescale, $t_{\mathrm{cool}}={\beta }_{\mathrm{cool}}{\mathrm{\Omega }}_k^{-1}$. The thermal relaxation is represented by the time evolution of the pressure dependent of the local density and the local sound speed

$$\begin{array}{c}\hfill \frac{\mathit{dP}}{\mathit{dt}}=-\frac{P-\rho c_{s,\mathrm{ini}}^2}{{\beta }_{\mathrm{cool}}{\mathrm{\Omega }}_k^{-1}},\end{array}$$(3)

which is relaxed to its initial value described by c_s, ini. The initial disk conditions adopt values typical of T Tauri systems, with q = −0.5, which has been shown to favor VSI growth (Manger et al. 2021), and p = −2.25, a choice that ensures a constant mass accretion rate is maintained (see Equation (13) in Manger et al. 2021), consistent both with theoretical expectations (Lynden-Bell & Pringle 1974) and with observational constraints (e.g., Wolf et al. 2003; Sauter et al. 2009). We adopt a cooling time of ${10}^{-4}{\mathrm{\Omega }}_k^{-1}$, which is sufficiently short to permit the development of the VSI in the outer disk (Table 1).

We adopted the same boundary conditions as in Manger & Klahr (2018) to avoid wave reflections and minimized boundary effects at the inner and outer edges of the domain. Specifically, within buffer zones of width ΔR = 1.0H and Δθ = 0.05, the density and velocity components were damped toward their initial values when the velocity normal to the boundary points inward into the domain. We used periodic boundaries in the azimuthal direction. Within these buffer zones, the damping timescale was set relative to the local orbital period as

$$\begin{array}{c}\hfill \frac{d{v}_x}{\mathit{dt}}=-\frac{v_x-v_{x,0}}{{\tau }_{\mathrm{damp}}}·f^2.\end{array}$$(4)

Here, v_x represents any of the primitive variables, such as density or velocity components. The damping timescale is defined as ${\tau }_{\mathrm{damp}}=0.1\frac{2\pi }{{\mathrm{\Omega }}_{\mathrm{k}}}$, and the damping factor is given by $f=\frac{R-R_{\mathrm{boun}}}{\mathrm{\Delta }R}$ in the radial direction, and by $f=max(\frac{R-R_{\mathrm{boun}}}{\mathrm{\Delta }R},\frac{\theta -{\theta }_{\mathrm{boun}}}{\mathrm{\Delta }\theta })$ in the meridional direction. The parameters R_boun and θ_boun denote the radial and meridional locations of the boundaries of the damping layer within the buffer zone.

3. Dust

In this paper, we are interested in understanding where dust-loading concentration occurs in VSI-active protoplanetary disks and determining qualitatively whether it has observational implications. From the 2D (R, Z) model of the density and rotational velocity described in Section 2, we introduced three populations of Lagrangian particles with sizes of 100 μm, 500 μm, and 1 mm in the simulation. We reset the dust component midway through the simulation runtime. This was necessary because, by approximately 200 orbits, the initial set of dust particles had already drifted inward before fully developed vortices – whose full development occurs beyond 250 orbits – could effectively trap incoming particles. This approach allowed us to investigate and analyze the role of vortices in dust trapping more thoroughly. The particle dynamics is described using a Lagrangian dust module (Mignone et al. 2019) in the PLUTO code using a spherical coordinates system. These dust grains are assumed to be homogeneous, compact, and spherical with a dust density of, ρ_dust = 1.7 g cm⁻³. We do not include the effects of the dust-to-gas feedback (see, e.g., studies of dust-to-gas feedback in VSI simulations by Schäfer et al. 2020; Schäfer & Johansen 2022). We assume that dust particles are sufficiently small compared to the gas mean free path, allowing treatment within the Epstein regime (Epstein 1924). In this regime, the Stokes number is given by St = t_sΩ_k, where t_s is the stopping time–the timescale over which a dust grain adjusts to gas motion, depending on its size. It dictates the time that it takes for a dust particle to be slowed down by the gas flow and in terms of the disk properties is defined as

$$\begin{array}{c}\hfill t_s=\frac{\mathrm{St}}{{\mathrm{\Omega }}_{\mathrm{k}}}=\frac{a{\rho }_{\mathrm{dust}}}{\rho (R,Z)H{\mathrm{\Omega }}_{\mathrm{k}}}=\frac{a{\rho }_{\mathrm{dust}}}{\rho (R,Z)c_s},\end{array}$$(5)

where a is the dust particle size, which is kept constant throughout the dynamical evolution of the disk (see Appendix C for the Stokes number values at different scale heights in the disk for each grain size). Aside from the Lagrangian particle sizes, which are considered large in our grain sample, we include small grains between 0.1 μm and 30 μm discretized in 10 bins important to consider when post-processing the simulation data in the RADCM3D code (see Section 5). These small grains are not modeled as Lagrangian particles in the VSI simulations, as they remain well coupled to the gas. For the dust mass distribution we follow the approach as in Ruge et al. (2016) where the dust-size distribution is defined based on the MRN (Mathis et al. 1977) distribution with a power-law of −3.5. The MRN distribution fits well for the larger particles that dominate the intensity at the chosen wavelengths (see Section 5). The small grains are very well-coupled with the gas and spatially distributed uniformly. We inserted a million of particles, and for the mentioned distribution, we obtained dust-to-gas mass ratio values for the simulation run: M_small, d/M_gas = 0.0003, M_{100 μm, d}/M_gas = 0.0012, M_{500 μm, d}/M_gas = 0.0026, and M_1 mm, d/M_gas = 0.0036.

4. Results

To assess the time-integrated turbulence in the disk, Figure 1 (top) shows the time evolution of the volume-weighted Reynolds stress-to-pressure ratio, α_r, ϕ, that is radially and azimuthally averaged

$$\begin{array}{c}\hfill {\alpha }_{r,\varphi }=\frac{\int T_{r,\varphi }dV}{\int PdV}=\frac{\int \rho v_r^{\prime }{v}_{\varphi }^{\prime }dV}{\int \rho c_s^{2}dV}·\end{array}$$(6)

Fig. 1. Reynolds stress-to-pressure ratio, αr, ϕ. The top panel shows the time-averaged αr, ϕ, where the yellow line represents the convergence value, ⟨αr, ϕ⟩ = 3 × 10−4. For reference, the vertical dotted line corresponds to the outer local orbit at 0.25 R0 used by Lesur et al. (2025) in their analysis, as shown in their Figure 1. The bottom-left panel displays the time series of αr, ϕ in the r-direction, averaged over azimuth in the 3D VSI-active disk. The bottom-right panel presents the radial Reynolds stress-to-pressure ratio over the time series and averaged in the azimuthal directions.

Here v_r′ and v_ϕ′ are the velocity perturbations calculated by subtracting the azimuthally averaged value at each r and θ. The saturation phase occurs around 100 orbits, during which α_r, ϕ converges to a value of 3 × 10⁻⁴, and the turbulence remains stable until approximately 500 orbits. We use a single snapshot at 350 orbits that shows the two large-scale dust crescents at similar locations motivated by those observed in the MWC 758 disk (Dong et al. 2018). Figure 1 (bottom) shows the time series of α_r, ϕ in the midplane, where an increasing trend is observed over time.

4.1. Vorticity and dust concentration

To identify dynamical gas structures acting as dust traps, we generated midplane vorticity maps in the R − ϕ plane. The vertical component of the vorticity, normalized by the Keplerian frequency, is given by

$$\begin{array}{c}\hfill \frac{{\omega }_z}{{\mathrm{\Omega }}_k}={(\mathrm{\nabla }\times \mathit{v})}_z=\frac{1}{{\mathrm{\Omega }}_k}(\frac{1}{r}\frac{\partial r{v}_{\theta }}{\partial r}-\frac{1}{r}\frac{\partial v_r}{\partial \varphi }),\end{array}$$(7)

where v is the gas velocity field (Lovelace et al. 1999). Figure 2 shows the resulting midplane vorticity map normalized by Ω_k, overlaid with contours of the dust-to-gas ratio (log₁₀ε = −1.0 and −2.0) for 1 mm particles at the midplane simply to illustrate where these larger particles tend to concentrate. Two large-scale anticyclonic vortices are evident at distinct azimuthal locations and they exhibit strong particle accumulation toward their centers.

Fig. 2. Vorticity of the gas in the midplane in units of the keplerian frequency at 350 orbits. Overplotted are two dust-to-gas mass ratio values of 1 mm particles in the midplane.

To demonstrate the evolution of vortices that are self-sustained within the VSI-unstable disk, we generated a time-resolved vorticity residual map (Figure 3). This was constructed by applying a box filter to the midplane vertical vorticity, ω_z, using a filter width of 200 grid cells in the azimuthal direction and 30 cells radially (approximately two pressure scale heights). The filtering suppresses small-scale fluctuations, enhancing coherent vortex structures. We then computed the residual as the deviation of the smoothed local vorticity from its azimuthal average within each annulus, allowing us to isolate vortex signatures from large-scale azimuthal flows.

Fig. 3. Time evolution of the minimum vorticity residual as a function of radius, illustrating the radial position and migration of vortices in the VSI-unstable disk. The horizontal dashed lines outline the final locations of the two large-scale vortices.

The inner vortex forms near a local pressure maximum (see Appendix D for more information about the pressure perturbations) at the inner damping zone (see Section 2), where VSI turbulence is intentionally suppressed to produce a low-viscosity transition region. Our vorticity residual analysis (Figure 3) shows that the vortex persists away from the damping region, indicating that it is most likely a self-sustained structure within the VSI flow. The outer vortex becomes fully developed within 60 au after approximately 250 orbits and gradually migrates inward at a rate of ∼0.04 au per orbit (≈9 × 10⁻⁴ R₀ per orbit) because of the steep radial pressure gradient (see more in Manger & Klahr 2018). We notice a vorticity field perturbation appearing around ∼100 orbits that seems to be associated with the origin of the outer vortex’s trailing structures, but no obvious vortex signatures appear in the snapshots between 50 and 100 orbits. Because its formation occurs after 250 orbits, when most 1 mm dust particles have already drifted inward, it captures fewer particles than the inner vortex. Its aspect ratio, χ ≳ 7, is roughly half that of the inner vortex. Additional small-scale vortex trails are also visible in Figure 3. These smaller structures can transiently concentrate particles–particularly those with sizes of 500 μm and 1 mm–but we do not investigate them further in this study, as our analysis focuses on the largest, long-lived vortices that dominate the disk dynamics.

We additionally verified that vortices are self-sustained in VSI disks using a model with the same density and temperature gradients but a different cooling parameter, ${\beta }_{\mathrm{cool}}{\mathrm{\Omega }}_k^{-1}=5\times {10}^{-3}$ (see Appendix B). The vorticity residual evolution in this model also exhibits coherent vortex signatures forming away from the damping zone, confirming that vortex generation is not restricted to the fiducial boundary setup.

4.2. Azimuthal and vertical velocity fields

To examine how the vortices modify the local gas dynamics, we analyzed the azimuthal and vertical velocity perturbations. Figure 4 shows the azimuthal and vertical gas velocity perturbations together with the vertically averaged dust distribution for the same snapshot as in Figure 2. The inner and outer vortices are traced by concentrations of 1 mm and 100 μm particles, respectively (indicated by arrows). Both vortices produce strong azimuthal velocity deviations from Keplerian rotation, exhibiting characteristic red–blue dipolar patterns that indicate anticyclonic motion. The outer vortex, located near 50 au, shows the clearest retrograde flow signature, with sub-Keplerian motion on one side and super-Keplerian motion on the other. These perturbations persist from the midplane to the disk surface (> 2H), though they become more pronounced at higher altitudes (Figure 4, bottom panel). Vertical velocities remain weak near the vortices but increase toward the surface, where the flow becomes more asymmetric and turbulent. The presence of velocity perturbations caused by vortices significantly alters the surrounding gas kinematics.

Fig. 4. Pertubational azimuthal velocity (top and bottom-left panels) and vertical velocity (top and bottom middle panels) of the disk at the midplane (z = 0, top panels) and close to the disk surface (z > 2H, bottom panels). To make it easier to identify the two large-scale vortices in the disk, we show the dust distribution vertically averaged in the top right panel. The black arrows point at the location of the inner and outer vortices.

4.3. Dynamical trajectory of particles inside a vortex

We further analyzed the trajectory of randomly selected particles, representative of a 100 μm-sized and 1 mm sized grain, to assess its velocity profile after being captured by the outer vortex (see Appendix E for their spatial 3D trajectories). To reduce computational expense, the analysis focuses on a limited orbital interval between 400 and 460 orbits, once the particles are already confined within the vortex. In Figure 5, we present the particle’s velocity components over time for a 100 μm-sized grain (left) and 1 mm sized grain (right). In both particles’ profile, we find that the vertical velocity no longer exhibits the characteristic upward and downward symmetry typically expected in VSI-dominated regions, as shown in Figure 6 and Figure 7 in Flores-Rivera et al. (2025), whereas the radial and perturbed azimuthal velocities largely retain their sinusoidal behavior. Additionally, we observe that the vertical velocity becomes subdominant, with a smaller amplitude compared to the radial component.

Fig. 5. Velocity profile of a single 100 μm-sized particle (left) and a 1 mm-sized particle (right) within the vortex.

We also conducted a test analysis of the trajectories of a 100 μm- and a 1 mm-sized particles trapped within the inner vortex. In both particles’ cases, the vertical velocity amplitudes are nearly a factor of two lower than those observed for particles in the outer vortex. Nevertheless, we observe a consistent trend in which the sinusoidal symmetry of the vertical velocity is broken and remains subdominant relative to the radial component.

4.4. Spatial distribution of the dust

To examine the spatial distribution of solids, Figure 6 shows the dust density for three grain sizes in the ϕ–R and Z–R planes at the same snapshot as in Figure 2. The 100 μm particles remain broadly distributed and well coupled to the gas, whereas the 1 mm grains occupy smaller disk radii and are concentrated within the inner vortex. The intermediate 500 μm particles are trapped in both the inner and outer vortices, where the later captures primarily the 100 μm and 500 μm grains after most 1 mm particles have drifted inward. Notably, small-scale vortices outside the pressure maxima also contribute to localized dust concentrations, especially for 1 mm particles (Figure 6, right panel).

Fig. 6. Dust density distribution. Top panels: Dust density distribution vertically averaged at 350 orbits in the ϕ − R plane. The large grains are 100 μm (blue), 500 μm (green), and 1 mm (orange). The small grains spanning 0.1–30 μm are plotted as a black shade in the disk and are assumed to be well coupled with the gas. Bottom panels: Same dust density distribution as the top panels but viewed in the Z − R plane and averaged in the azimuthal direction.

The vortices are more difficult to visualize in the vertical distribution of particles (Figure 6, bottom row). The vertical dust distribution of the 1 mm particles shows that these particles are more settled toward the midplane (Z/R < 0.1), while 500 μm grains exhibit slightly greater vertical extent. At the outer vortex location, the 500 μm particles form slightly more compact radial and vertical concentration, whereas the 100 μm grains appear more dispersed. In contrast, the inner vortex is less prominent in the vertical dust profile, with the 1 mm distribution appearing relatively flat. Figure 7 confirms these trends, showing azimuthally averaged midplane dust-to-gas mass ratios: the 1 mm and 500 μm particles exhibit the highest enhancements, consistent with the vortex trapping locations.

Fig. 7. Dust-to-gas mass ratio at the midplane at 350 orbits. The dotted lines represent the initial dust-to-gas mass ratio value for each grain size.

To quantify the vertical distribution of solids, we measured the vertical position, Z, of all particles radially averaged between 40 and 50 au, after 300 orbits. Figure 8 shows the normalized vertical dust distributions for each grain size (solid lines), with Gaussian fits overplotted (dashed lines). At 50 au, we find $H_d^{1\mathrm{mm}}/R=0.0325\pm 0.0003$, $H_d^{500\mathrm{\mu }\mathrm{m}}/R=0.0404\pm 0.0002$, and $H_d^{100\mathrm{\mu }\mathrm{m}}/R=0.0562\pm 0.0001$, yielding dust-to-gas scale height ratios of H_d/H = 0.32, 0.40, and 0.56, respectively.

Fig. 8. Mean vertical dust distribution, Z, averaged over the radial range 40–50 au and over time after 300 orbits. Colored lines show the mean normalized distributions for 1 mm (orange), 500 μm (green), and 100 μm (blue) grain sizes. The standard deviation of each distribution corresponds to the dust scale height (Hd) for the respective grain size.

The time evolution of the dust scale height, shown in Figure 9, reveals that 500 μm and 1 mm particles gradually settle and drift inward, while the 100 μm grains remain more vertically extended. After ∼350 orbits, both larger grain sizes exhibit strong fluctuations in H_d, with additional variations around ∼420 orbits. We note that the strong fluctuations seen after ∼450 orbits are not caused by an artificial inverse cascade or a breakdown of the turbulence. Instead, they arise when the large outer vortex migrates into the 40–50 au measurement region, temporarily dominating the dust distribution and producing non-axisymmetric vertical oscillations. These oscillations arise from vertically stirred, compact dust concentrations within vortices, as seen in the R–Z and R–ϕ projections (bottom panels of Figure 9).

Fig. 9. Dust scale height over time for all particle species. The dust scale height is radially averaged from 40 au to 50 au. Each marker represents the averaged dust scale height calculated by taking the rms of the vertical position of the particles for all sizes per orbit. Strong vertical fluctuations are observed for the 1 mm and 500 μm grain sizes after approximately 424 orbits. This corresponds to a regime where the dust distribution is characterized by localized and compact concentrations within vortices (see the orange region in the bottom three panels at three different snapshots).

During these fluctuation phases, the vertical dust distribution deviates from a Gaussian profile, forming either a sharply peaked layer or a bimodal structure. In such cases, Gaussian fits fail to capture the true vertical extent, making the rms of the particle vertical positions a more robust diagnostic of H_d. After ∼460 orbits, the 1 mm particles have largely drifted inward, leaving only small remnants trapped in the vortices, visible as narrow, non-Gaussian dust layers (orange regions in Figure 9).

5. Synthetic images at different wavelength bands

To generate the intensity maps of the dust continuum emission of our disk model, we used RADMC-3D² (Dullemond et al. 2012). The stellar parameters used are characteristic of a T Tauri star to explore the general physical conditions under which such features may arise. We used T_★ = 4000 K, M_★ = 0.5 M_⊙, and R_★ = 2 R_⊙.

We generated the dustkapscatmat_dsharp_*.inp files that contains Mie-based scattering matrices for compact, spherical grains made of 36% water, 17% astronomical silicates, 3% troilite, and 44% refractory organics by volume (Birnstiel et al. 2018). The opacity values are made by the DSHARP_OPAC package using the Mie code of Bohren and Huffman (BHMIE; Bohren & Huffman 1983) embedded in a python script. The opacities are computed for all the dust populations considered here, ranging from 0.1 μm to 1 mm, under the standard power-law distribution. We adopted anisotropic scattering by setting scattering_mode = 2 in the radmc3d.inp file, which enabled the Henyey–Greenstein approximation using the scattering opacities and their corresponding anisotropy parameters provided in the dust opacity file.

The synthetic images of our modeled protoplanetary disk, are shown in Figure 10 (top panels) at four different ALMA wavelength bands: 340 μm (Band 9), 650 μm (Band 8), 870 μm (Band 7), and 1.3 mm (Band 6). As a test, we also generated synthetic images at longer wavelengths, including 3.0 mm (the lowest-frequency end of the Next Generation Very Large Array, ngVLA), 7.0 mm, and 10 mm (the highest-frequency end of the ngVLA). At these wavelengths, the two bright crescent-like features remain prominent, but the emission appears more compact, consistent with dust trapping in the inner and outer vortices identified in previous figures. The rest of the disk becomes fainter, as the ngVLA bands are most sensitive to grains between 5 and 15 mm in size – particles not included as Lagrangian components in our simulations – making it noteworthy that the modeled 1 mm grains still produce sufficient flux to weakly reveal the vortices. The dust crescents are most prominent in Bands 8, 7, and 6. At shorter wavelengths, such as Band 9, the brightness contrast of the crescents is considerably reduced compared to the longer-wavelength bands. Moreover, the outer crescent is no longer visible, and spiral arms are also absent, consistent with the lack of spirals previously seen in the dust density distribution.

Fig. 10. Synthetic images of the simulated disk at 0° inclination. The top row shows the brightness of the simulated disk after radiative transfer for four different ALMA wavelength bands (columns). The bottom row shows the convolution of the disk images using a circular beam with a FWHM of 12 au (σbeam = 5 au). Overall, the presence of two asymmetric dust features are clearly seen at the inner ring and in the outer parts of the dust continuum emission. The inner one was found to be at the same location as the dust concentration.

The bottom panels of Figure 10 show the synthetic intensity maps at the same wavelengths bands but convolved with a synthetic 2D circular Gaussian beam of σ_beam = 5 au (equal to a circular beam with a full width at hald maximum (FWHM) of 12 au, similar to the synthesized beam used in Dong et al. (2018) for MWC 758). For visualization purposes, the color scale in each panel was adjusted to enhance the contrast between the different ALMA bands. The previously sharp and intricate substructures observed in the unconvolved images such as multiple narrow rings, crescents, and azimuthal asymmetries, are now smoothed out due to the finite beam size. Despite the convolution, prominent features remain visible in Bands 8, 7, and 6, where the two asymmetric crescents are still present. In Band 9, the inner dust crescent appears over-brightened but takes on a more ring-like morphology due to the lower brightness contrast at this wavelength, while the outer crescent essentially disappears. For the ngVLA wavelength bands, these features become increasingly diffuse and faint at longer wavelengths (≥3.0 mm) and are nearly undetectable at 10 mm, consistent with reduced dust emission and beam dilution of small-scale structures. Compared to the original unconvolved images, this highlights how observational resolution and sensitivity significantly affect the ability to detect and characterize such fine-scale substructures in vortices. However, our synthetic images do not include noise which is another important component to consider in order to accurately determine if the outer vortex, for example, is resolved or not. We also understand that while the absolute contrast in dust emission may vary with the seeding time, the qualitative morphology, particularly the asymmetric structures associated with the vortices, will remain robust.

5.1. Brightness contrast and aspect ratio

There are two main physical properties we can analyze from our disk imaging that can be directly compared with those estimated by ALMA observations from disks with dust crescents. The dust crescents in disks observed in submillimeter by ALMA are identified by their large brightness contrast. We define the brightness contrast ratio as $\frac{I_{\nu ,\max }}{I_{\nu ,\min }}$ of the dust crescent as a function of radial location. Using the same time snapshot as in previous analyses, we calculated the brightness contrast of the brightest inner dust crescent, adopting the 85% peak intensity level for all ALMA bands shown in Figure 10, except for Band 9. Because Band 9 exhibits a much lower brightness contrast and appears as a nearly smooth, ring-like structure with little evidence of substructure, we instead used a 95% threshold; using 85% would no longer capture a crescent-shaped feature but rather a continuous ring. We find contrast ratios of 4.23 (Band 6), 3.93 (Band 7), 3.57 (Band 8), and 2.49 (Band 9), indicating that the contrast decreases toward shorter wavelengths. From Band 9 to Band 6, I_ν, max drops by ∼21%. As shown in Figure 10, the contrast between the crescent-shaped features and the surrounding disk increases with observing wavelength, while the overall intensity declines slightly at longer wavelengths, reflecting both the preferential tracing of larger, more concentrated dust grains and the reduced thermal emission.

The second physical property is the aspect ratio, χ defined earlier as the ratio between the azimuthal extent to the radial extent of the dust crescent, as a function of its radial location. For this, we fit ELLIPSEMODEL() from the SKIMAGE.MEASURE module in Python which uses the least-square geometric fitting method to estimate the parameters of an ellipse from the X and Y pixels points in the image data. The ellipse model uses the spatial distribution of the dust crescent pixels corresponding to 85% of I_ν, max for all bands (and 95% for Band 9), and computes the best-fit ellipse by minimizing the algebraic distance between the image points and the conic representation of an ellipse. We use the lengths of the semi-major and semi-minor axes as the azimuthal and radial extent for χ. We estimate the aspect ratios for the four ALMA wavelength bands to be: 5.1 (Band 9), 4.0 (Band 8), 3.8 (Band 7), and 3.6 (Band 6). By dividing the estimated semi-major axis by σ_beam, we find that the inner vortex is well resolved in Band 6 (ratio ∼3.6) and in Band 9 (ratio ∼2.6). Our estimates of the brightness contrast and aspect ratio are further compared and discussed in the context of the observed estimations in Section 6.4.

5.2. Effects due to disk inclination

So far, we have analyzed the observational features of our model assuming a disk inclination of 0°. In Figure 11, we explore how the appearance of the dust crescents changes with increasing inclination. We find that the crescent becomes progressively more elongated and can resemble a ring-like structure at higher inclinations. This projection effect has been reported in previous studies (Doi & Kataoka 2021; Villenave et al. 2025), where optically thin, vertically extended, and radially narrow dust rings exhibit azimuthal over brightness asymmetries w. In particular, emission is enhanced along the projected semi-major axis because the line-of-sight intersects a longer path through the vertically extended dust layer, while the semi-minor axis samples a shorter vertical column and therefore appears fainter. Similar inclination-dependent azimuthal smearing of dust substructures has also been discussed in recent observational analyses (e.g. Ribas et al. 2024). As seen in the highest-inclination case in Figure 11, this geometric effect can partially wash out crescent contrasts. This result is encouraging, as many disks observed in ALMA Bands 6–7 have moderate to high inclinations, implying that vortex-induced crescents may remain undetected in some systems due purely to projection effects.

Fig. 11. Dust crescent at different disk inclinations. The black contour represents 85% of the peak intensity for the ALMA band 6 frequency. The elongation of the dust crescent increases with disk inclination due to an imaging projection effect. This effect is also evident in the other ALMA wavelength bands.

We further analyze the radial intensity profiles of the images for all disk inclination cases to identify any potential intensity bumps at the location of the dust crescent. For the 0° inclination case, which shows the clearest intensity features, the maximum peak corresponds to the location of the pressure bump as in Figure D.1 (bottom). While the intensity generally decreases with radius, a second peak at the location of the outer vortex is not evident, despite the presence of minor fluctuations. Similarly, the radial intensity profiles for disks at higher inclinations exhibit a smoother distribution.

5.3. Velocity residuals around the dust asymmetries

To post process the simulation data in the gas and explore observable features to compare with ALMA observations, we used RADMC-3D and performed a line transfer of a line velocity of interest assuming local thermodynamic equilibrium (LTE) conditions and following the methodology presented in Barraza-Alfaro et al. (2021). We used the ¹²CO (J = 2 − 1) line transition centered at 230.538 GHz. For the molecular number density, we assumed a constant fraction of ¹²CO relative to H₂ of 1 × 10⁻⁴. The radiative properties of the CO molecular data was done with the Leiden LAMDA database³ (Schöier et al. 2005). We also considered the 3D structure of all velocity components and assumed that the dust temperature is equal to the gas temperature. The synthetic data cubes have a total bandwidth of 6 km s⁻¹ distributed over 40 channels, each with a spectral resolution of 0.15 km s⁻¹, as produced by RADMC-3D.

Figure 12 displays the moment 1 maps of the velocity residuals, v₀, obtained by subtracting the background Keplerian rotation profile from the simulated gas velocity field for the same four disk inclinations shown in Figure 11. These maps were computed using the BETTERMOMENTS⁴ Python package (Teague & Foreman-Mackey 2018), which applies a quadratic fitting method centered on the brightest pixel of the point-spread function to more precisely estimate the line centroids of spectral profiles. Once the pure Keplerian velocity is subtracted from the line-of-sight velocity maps, significant deviations from Keplerian rotation emerge across all inclinations. For the lowest inclination (first panel), concentric ring-like features of sub-Keplerian residuals dominate the disk. At the location of the outer vortex (Y = 50 au), a characteristic quasi-dipolar red-blue pattern is evident, indicating faster-moving gas (blue) following the disk’s rotation, and slower gas (red) on the opposite side moving against it. In contrast, Keplerian deviations near the inner vortex are more difficult to discern; however, as the inclination increases, this red-blue pattern becomes increasingly apparent in both vortex regions. Outside the vortices, the concentric ring-like Keplerian deviations remain somewhat symmetric, whereas in the vortex-hosting regions, this symmetry is clearly disrupted. We do not observe clear spiral features in the velocity residuals near the vortex regions. Future studies should analyze the velocity residuals of CO isotopologues to more thoroughly investigate the gas dynamics surrounding vortices.

Fig. 12. Velocity residuals of the simulated disks at different disk inclinations. The same beam size used in Figures 10 and 11 was applied to the convolved images (bottom panels).

6. Discussion

6.1. Model limitations

In terms of thermodynamics, the simulation employs a simplified cooling prescription that is short enough to allow the growth of the VSI. We also conducted additional simulation runs using different cooling times to investigate the behavior of vortices in dust trapping under both optimal and marginal VSI conditions; these results will be presented in a forthcoming paper. The most recent state-of-the-art studies addressing the vertical thermal structure in VSI-unstable disks incorporate radiative transfer into radiation hydrodynamics (RHD) models. For instance, Melon Fuksman et al. (2024a) evolve the frequency-integrated radiation energy flux, while Zhang et al. (2024) implement angle-dependent (anisotropic) radiative transfer to more accurately model radiation transport in both optically thin and thick regimes. Both studies find that the VSI results in a quiescent midplane, with turbulence primarily confined to the disk’s upper layers. These approaches have been consistently applied in VSI-unstable disk models to explore gas–dust dynamics under varying gas cooling prescriptions, opacities, and molecular emissivities, concluding that VSI growth is strongest away from the midplane, at intermediate heights where the vertical shear is largest (typically Z/R ≈ 0.1 − 0.2; Flock & Mignone 2021; Melon Fuksman et al. 2024b). More advanced radiative transfer schemes, such as those used by Muley et al. (2023), further resolve the energy exchange between gas, dust, and radiation, underscoring their potential for future applications in modeling hydrodynamic instabilities.

We limited our analysis to processes occurring prior to dust aggregation and the initial formation of planetary embryos. While this study demonstrates that vortices can generate dust asymmetries in submillimeter continuum emission within turbulent protoplanetary disks, further VSI simulations are necessary to investigate the longevity of these vortices under the influence of dust backreaction. The formation of planetesimals is favored in scenarios where the SI develops following the onset of the VSI in an isothermal configuration, with strong clumping occurring for dust concentrations as low as 0.5% (St = 0.1) or 1% (a = 3 mm) (Schäfer et al. 2020; Schäfer & Johansen 2022). However, by changing the equation of state, gas buoyancy can affect the VSI in the disk midplane (Lin & Youdin 2017; Lin 2019).

6.2. Resolution effect

It is worth noting that a recent high-resolution study by Shariff & Umurhan (2024) and Lesur et al. (2025), reaching up to 70 and 250 cells per scale height, respectively, reported no vortices in the midplane of their VSI simulations using an effectively isothermal setup. Instead, these simulations only show small-scale, isotropic fluctuations below the disk scale height. They suggest that this absence of vortices could be linked to the elliptical instability (E.I.; Lesur & Papaloizou 2009), not present in lower resolution runs of previous works with 18-50 cells per H; e.g., Manger & Klahr 2018; Richard et al. 2016; Huang & Bai 2025), including the present work (37 cells per H). We note that the simulations in Lesur et al. (2025) report evolution in inner orbits measured at R = 0.25 R₀, which corresponds to < 187 local orbits for the longest resolution which is significantly shorter than the timescale over which vortices emerge in our models. Regarding the resolution used in Shariff & Umurhan (2024), it is difficult to directly compare their resolution with ours because they employ a Padé filter, a high-order numerical smoothing operator that selectively damps unresolved small-scale fluctuations near the grid scale while preserving the resolved large-scale flow; its operation also requires an accompanying shock-viscosity treatment. We encourage future studies to investigate the occurrence of long-lived vortices in high-resolution models across a wide range of cooling timescales, since vortices may persist at longer cooling timescales (Richard et al. 2016, Sengupta et al., in prep.)

We investigate the kinetic energy spectrum of our simulations by extracting a smaller box region of the domain centered at 1.25 in code units, taking r ∈ [1.1, 1.4] in code units, and θ ∈ [1.57, 1.81], and ϕ ∈ [3.14, 3.68] in radians. To perform a fast Fourier transform (FFT) on the simulation data, the grid points must be equally spaced so that the transformed fields satisfy the periodicity assumptions of the FFT. Because the radial grid in our case is non-uniform, we first resample it onto a uniform grid using linear interpolation. We then apply a Hanning window (i.e., Lesur et al. 2025) to the selected box to reduce edge effects in case the boundaries are not periodic. After windowing, we multiply the simulation variables by the window and compute the Fourier transform of the momentum components in each direction, followed by the Fourier transform of the velocities. The resulting energy spectrum is defined as the sum of each velocity Fourier component multiplied by the complex conjugate of its corresponding mass-flux Fourier component. The 3D spectrum is obtained by binning the Fourier-space energy into spherical shells using evenly spaced wavenumber bins, summing the total energy, E(k), in each shell, dividing by the shell thickness to obtain an energy density, and finally normalizing the spectrum so that its integral over all wavenumbers equals unity. This methodological procedure follows the approach described in Appendix B of Melon Fuksman et al. (2024b).

From Fig. 13 we observe the energy spectrum together with two reference slopes: a Kolmogorov-like power-law slope (green) and a much steeper power-law slope (red). The Kolmogorov-like regime is indicative of a turbulent cascade; however, in anisotropic, rotating, and vertically stratified flows such as VSI turbulence, this slope does not uniquely imply a forward cascade from large scales. Energy may instead be injected near the fastest-growing or marginal VSI wavelengths and redistributed both upscale and downscale. By estimating the scale-dependent Rossby number, we find that Ro < 0.1 near the spectral kink and decreases toward both larger and smaller scales, indicating that rotation increasingly constrains the flow away from the injection scale. This behavior is consistent with quasi-2D rotating turbulence, in which energy is preferentially transferred upscale while enstrophy cascades downscale, naturally explaining the coexistence of a Kolmogorov-like slope and a steeper high-k regime (see also Manger & Klahr 2018). Here k = 2π/λ denotes the Fourier wavenumber associated with spatial scale λ. The two vertical gray lines mark the characteristic VSI wavelengths in the Fourier space, expressed in units of kH. Structures with wavelength equal to the disk scale height (λ = H) appear at kH = 2π, while the fastest-growing and marginal modes located at kH = 40 (dotted-line) and kH = 20 (dashed-line), respectively. We theoretically calculate the fastest-growing VSI wavelength using the prescription from Umurhan et al. (2016), further supported by the numerical analysis of Shariff & Umurhan (2024), where $\frac{\lambda }{H}=\pi |q|{\left(\frac{H}{R}\right)}^2$ = 0.16. Using the local scale height at the spectrum center, we obtain λ = 0.020 with the marginal mode at 2λ. Converting to Fourier space, this yields the characteristic wavenumbers marked in vertical gray lines in Fig. 13. These characteristic scales lie close to the transition toward the dissipation range, where the spectrum steepens sharply. This proximity suggests that part of the VSI-injected energy may be removed numerically rather than fully cascading through a resolved inertial range, although the presence of both Kolmogorov-like and steeper slopes indicates a turbulence regime strongly influenced by shear, rotation, and vertical stratification.

Fig. 13. Energy spectrum for a small box extracted from the simulation domain (blue curve) and for the full domain (orange curve) at 400 orbits. Overplotted are two power laws, k−5/3 and k−5 illustrating the Kolmogorov regime and a steeper dissipation regime, respectively. The dotted and dashed vertical gray lines mark the wavenumbers corresponding to the marginal and fastest-growing VSI modes, based on the prescription from Umurhan et al. (2016). The dash-dotted vertical line denotes to the mode whose wavelength equals the local disk scale height λ = H, corresponding to kH = 2π.

By multiplying $\frac{\lambda }{H}$ by the number of cells per scale height shows that the fastest-growing modes are resolved within roughly 6 grid points, meaning that the power injected at that spatial scale can accumulate and contribute to the generation of the vortices. The fact that the energy injection occurs close to the dissipation range implies that VSI energy is dissipated locally and turbulence remains underdeveloped at small scales (large k). However, in our case the vortices remain self-sustaining due to continuous VSI driving throughout the simulation time. Thus, early damping of VSI motions effectively prolongs vortex survival relative to what would occur in a fully resolved turbulent cascade.

6.3. Dust vertical diffusion

The VSI drives turbulent vertical gas motions that play a key role in determining the extent of dust vertical diffusion in the disk. We estimate the degree of vertical stirring, α_z, from the directly fitted value of H_d/R_50 au obtained from the dust grain distribution (Figure 8) and the Stokes number for all grain sizes at the midplane (see Appendix C). The vertical stirring is defined as

$$\begin{array}{c}\hfill {\alpha }_z=\mathrm{St}{[{\left(\frac{H}{H_d}\right)}^2-1]}^{-1},\end{array}$$(8)

(Dubrulle et al. 1995) yielding α_z = 1.7 × 10⁻³ with St = 1.4 × 10⁻² for the 1 mm grains, α_z = 1.3 × 10⁻³ with St = 7.1 × 10⁻³ for the 500 μm grains, and α_z = 6.4 × 10⁻⁴ with St = 1.8 × 10⁻³ for the 100 μm grains. Because dust backreaction is absent, the gas turbulence is identical for all grain sizes; however, the inferred α_z values reflect how each Stokes number responds to that same flow. In this sense, α_z should be interpreted as a grain-dependent measure of vertical stirring efficiency rather than a unique turbulent viscosity of the gas. We note, however, that the Stokes number entering Eq. (8) implicitly depends on the characteristic turbulence correlation or eddy turnover timescale, τ_eddy. In VSI turbulence this timescale is not well constrained. The commonly adopted definition St = t_stopΩ_k assumes that the eddy turnover time is of order Ω_k⁻¹. More generally, St = t_stop/τ_eddy, such that variations in τ_eddy would modify the inferred stirring efficiency. Therefore, the α_z values derived here should be regarded as effective diagnostics under this standard assumption rather than absolute measures of the intrinsic gas diffusivity.

The dust scale heights we measured (see Section 4.4 and Figure 9) are comparable to the vertical wavelengths of the VSI modes. This is consistent with the nonlinear VSI dynamics described by Latter & Papaloizou (2018), who showed that vertical roll-up of VSI modes produces coherent overturning motions capable of vertically stirring dust on these same scales.

Our estimate of the dust scale height, H_d^1mm/R_50 au = 0.032, is very similar to the value reported in RHD VSI simulations by Flock et al. (2020), H_d = 0.034, for the same grain size. Consistently, we obtain α_z/St = 0.12 for the 1 mm particles, close to the value of 0.10 reported by Flock et al. (2020), indicating that turbulent stirring dominates over vertical settling in both cases. Since the VSI strength is highly sensitive to initial disk conditions and the adopted thermodynamic treatment – particularly the cooling times with radius (e.g., Richard et al. 2016; Manger & Klahr 2018; Pfeil & Klahr 2021) – further exploration of dust distributions under different initial setups is needed.

Villenave et al. (2025) proposed that the VSI is not active in disks with well-settled dust scale heights, corresponding to α_z/St ≤0.1, a threshold motivated by the finding of Stoll et al. (2017) that α_z = 650α_r under isothermal conditions. In 2D isothermal VSI simulations, Flores-Rivera et al. (2025) – using the same disk profile (p = −1.5, q = −1) as Stoll et al. (2017) – found a dust scale height of H_d = 0.033 for 1 mm particles, with α_z/St ≈0.12, consistent with a regime of efficient VSI-driven mixing. This agrees with the trend shown in Figure 5 of Villenave et al. (2025), where α_z/St values below 0.1 correspond to weak vertical stirring and diminished VSI activity. Notably, their fits were obtained for 88 μm grains, comparable to the 100 μm case in our study, for which we find α_z/St ≈0.36, again in line with their results. By contrast, Flock et al. (2020), who included radiative transfer in their VSI simulations, reported α_z/St ratios closer to the marginal threshold proposed by Villenave et al. (2025). These comparisons highlight that VSI strength is highly sensitive to both disk structure and thermodynamic treatment, underscoring the need for further simulations with dust that systematically explore a wider range of initial conditions and realistic physics to refine this threshold.

We further investigated the eddy timescale governing particle vertical diffusion. This analysis also provides an empirical estimate of the turbulence correlation time, which is otherwise uncertain in VSI-driven flows and implicitly assumed in the Stokes number definition used above. From the vertical settling-diffusion equilibrium definition (Dubrulle et al. 1995; Youdin & Lithwick 2007)

$$\begin{array}{c}\hfill {\left(\frac{H_d}{H}\right)}^2=\frac{D_z}{\mathrm{St}{c}_{s}H}\end{array}$$(9)

where D_z is the dust vertical diffusion related to the eddy turnover time defined as Fromang & Papaloizou (2006)

$$\begin{array}{c}\hfill D_z=⟨v_z^2⟩{\tau }_{\mathrm{eddy}}\to {\tau }_{\mathrm{eddy}}{\mathrm{\Omega }}_k=\frac{\mathrm{St}{c}_s^2}{⟨v_z^2⟩}{\left(\frac{H_d}{H}\right)}^2.\end{array}$$(10)

Outside the vortices, we evaluated the averaged vertical velocity of particles over the radial range 40–50 au, time-averaged after 300 orbits. For 1 mm grains, we use the averaged Stokes number between (Z/H_R) = 0 and (Z/H_R) = 1 which is 0.0185 and the vertical velocity of the particle is $⟨v_z^2⟩/c_s^2=7.2\times {10}^{-4}$, the estimated eddy timescale for a particle inside a vortex is τ_eddyΩ_k = 2.63. For comparison, Stoll & Kley (2016) reported a value of 0.2, while Flock et al. (2020) obtained 20 based on α_z = 5.4 × 10⁻³ and $⟨v_z^2⟩/c_s^2=2.8\times {10}^{-4}$. Our results are consistent with the theoretical framework of Sengupta et al. (2024), who argued that anisotropic turbulence naturally produces Schmidt numbers different from unity. By inferring the eddy turnover time directly from the dust-diffusion properties in our VSI simulations, we obtain τ_eddyΩ_k∼ 1–3, placing our system the Ω_L > Ω_k regime identified in their work.

Inside the vortex, the vertical diffusion (Equation (8)) is estimated to be α_z, vortex = 5.8 × 10⁻⁴ using a dust scale height of H_d/R = 0.02 or (H_d/H)² = 0.04 (see Figure 9). The corresponding eddy timescale is 1.17, calculated from a value of $⟨v_z^2⟩/c_s^2=6.3\times {10}^{-4}$ (Figure 5). The α_z, vortex value is a factor of about three times lower than previously determined, more global value of α_z = 1.7 × 10⁻³ for the 1 mm grain size case. These very low values suggest that vortices preferentially enhance diffusion in the radial and azimuthal directions rather than vertically.

6.4. Dust crescents and disk properties

Dust crescents are readily identifiable in submillimeter observations due to their characteristic azimuthal brightness variations within the dusty ring. To date, a total of 19 dust crescents have been identified across 13 disks (Bae et al. 2023), suggesting that they may be relatively rare features in protoplanetary disks. However, as shown in Figure 11, when exploring crescent morphology under varying disk inclinations, we find that the apparent asymmetry of the dust crescent changes significantly with increasing inclination due to geometric projection effects discussed in Section 5.2. At high inclinations, the crescent-like morphology can become so distorted that it no longer resembles a crescent. Additionally, brightness contrast varies with inclination; higher inclinations reduce the line-of-sight interception of photons through the semi-minor axis, making that region appear fainter. We also predict that dust crescents traced by large particles in the disk midplane are unlikely to be significantly affected by scattering phase functions, as scattering is not expected to play a major role at such depths.

Given the limited number of dust crescents identified in the literature, little is known about how their size and brightness depend on global disk properties. The aspect ratio of dust crescents is thought to be a key parameter for constraining the physical mechanisms responsible for their formation. Although our estimates of the aspect ratio at 0° disk inclination fall within the range reported for observed dust crescents at similar radial locations (see Figure 5 in Bae et al. 2023), both the radial widths (ranging from 7 to 80 au) and azimuthal extents (40–190 au) exhibit significant variation. Consequently, there appears to be no clear correlation between aspect ratio and radial location, especially when projection effects due to disk inclination introduce further degeneracies. Bae et al. (2023) also report that the occurrence rate of crescents is approximately 30% for disks with M_disk/M_⊙ < 0.05. Our simulated aspect ratios are consistent with this trend, as our model disk has a mass of 0.038 M_⊙. Regarding crescent multiplicity, three disks to date – MWC 758 (Dong et al. 2018), HD 143006 (Benisty et al. 2018), and HD 139614 (Muro-Arena et al. 2020) – have been found to host more than one crescent, consistent with our predictions that multiple vortices can coexist in a single disk.

Other works by Flock et al. (2020) and Blanco et al. (2021) have explored and characterized the observational predictions and substructures in the distribution of millimeter-sized and small particles in 3D RHD models of the VSI. Both studies demonstrate that a vortex generated by the VSI can trap dust and create a localized azimuthal bump in the dust continuum emission, particularly at the inner edge of a planet-induced gap. The resulting dust concentration is more compact and produces a higher brightness contrast for larger grains, leading to stronger emission features at longer wavelengths – consistent with our findings. Additionally, they found that VSI-induced dust crescents exhibit radial widths of ∼2 au, which may serve as observational signatures of the VSI in protoplanetary disks.

Whether arising from the RWI or the KH instability, vortices are expected to be ubiquitous in both VSI-active and non-VSI disk environments. To guide observers, we propose that high-angular-resolution observations (≥14 au, as in the exoALMA program) of velocity residuals can help distinguish between vortices triggered by planets and those arising independently. If spiral features are detected in velocity residuals or in NIR scattered-light images, the RWI vortex could likely be excited by a giant planet carving a gap.

6.5. Dust-to-gas feedback effect

There have been relatively few studies exploring the numerical effects of dust-induced feedback on gas dynamics within vortices. Johansen et al. (2004) provided the first 3D demonstration that anticyclonic vortices can efficiently trap dust in protoplanetary disks, highlighting their potential role as long-lived sites of planetesimal formation and paving the way for connecting vortex theory with future ALMA observations. Meheut et al. (2012b) notably found that the vertical velocity structure within a vortex differs at higher disk altitudes. In their simulations, small particles (< 2 cm) experienced vertical lifting due to a positive feedback effect on the vertical velocity field. Fu et al. (2014) investigated dust back-reaction in a planet-induced vortex and found that when the local dust-to-gas mass ratio approached unity, a dynamical instability could be triggered within the vortex. This instability altered the vorticity and led to a reduced vortex lifetime. In contrast, Surville et al. (2016), using very high-resolution simulations (N_r × N_ϕ = 2048 × 4096), concluded that vortices could survive up to 1000 orbits even when accounting for dust back-reaction, eventually transitioning into ring-like structures that may become susceptible to the SI and enable rapid planetesimal formation (see also Raettig et al. 2015). Similarly, Miranda et al. (2017) performed HD simulations including dust feedback and produced synthetic observations, finding that while dust feedback may weaken vortices and slow dust accumulation, the dust crescents remain detectable for thousands of orbits, unlike in the planet-induced vortex scenario. They attribute this to the continual mass accumulation at the gap edge, which replenishes the vortex and prevents its dissipation. More recently, Huang & Bai (2025) conducted 3D high-resolution VSI simulations (with 230 cells per scale height) and found evidence of dust clumping, with approximately 2% of the dust mass trapped in vortices potentially undergoing gravitational collapse. They also identified small-scale dust clumps forming in azimuthal vortices within shearing zonal flows, which may be associated with Kelvin–Helmholtz (KH) instabilities. These findings emphasize that VSI-driven turbulence is not purely diffusive, but instead enables interactions with multiple hydrodynamic instabilities. In our work, we similarly confirm the interplay of three key instabilities: the VSI, which dominates the global gas dynamics; the RWI, which supports long-lived, large-scale vortices; and the KH instability, which emerges at the shear interfaces of zonal flows.

6.6. Implications to planet formation

Vortices are widely considered prime sites for planetesimal formation due to their efficiency in trapping pebbles, thereby enhancing the local dust-to-gas mass ratio and enabling the growth of planetary embryos (Johansen et al. 2004; Lyra et al. 2008). Lyra et al. (2009a) conducted global 2D simulations incorporating self-gravity and found that solid material accumulated within vortices and subsequently grew through collisions and pebble accretion (Lambrechts & Johansen 2012), ultimately can form Earth-mass planets. Later, Lyra et al. (2009b) demonstrated that vortices could be sustained in 3D through the RWI at gap edges where a viscosity transition occurs due to a Jupiter-mass planet opening a gap. When a Jupiter-mass planet grows over a sufficiently long timescale (e.g., > 700 orbits), the vortex triggered at the gap edge can become highly elongated (extending over approximately 180° in azimuth), which appears to be a promising characteristic for distinguishing dust asymmetries produced by a vortex formed due to planet-induced gap opening (Hammer et al. 2017, 2019). The persistence of RWI-induced vortices in three dimensions was subsequently confirmed by multiple studies (Meheut et al. 2010, 2012b,a,c; Lin 2012a,b, 2013; Lyra et al. 2015). In our case with no self-gravity, the vortices are able to trap a fraction of the Moon mass around their center.

Vortices offer a compelling solution to the drift barrier in planet formation (Birnstiel et al. 2010). While micrometer-sized dust particles can coagulate into millimeter-sized pebbles, further growth is inhibited by rapid radial drift toward the central star. This drift, dependent on particle size (Weidenschilling 1977), arises because gas in the disk, supported by pressure, rotates at sub-Keplerian speeds, causing dust particles to experience drag and spiral inward. However, since vortices modify the local pressure distribution in the disk by creating pressure maxima, they act as dust traps that draw grains toward their centers of rotation, preventing them from drifting inward and being lost. The dust density enhancement at the vortex (Barge & Sommeria 1995; Klahr & Bodenheimer 2006; Meheut et al. 2012b) provides favorable conditions for gravitational collapse or the streaming instability (Raettig et al. 2015, 2021), to ultimately promote planetesimal formation.

6.7. Vortices as a storage reservoir for CV chondrites

In cosmochemistry, the initial solid inventory of the protoplanetary disk is represented by chondritic meteorites–fragments of primitive, undifferentiated asteroids primarily composed of chondrules and fine-grained dust that formed in other regions of the disk before being accreted onto the parent body (Davis 2014). The accretion ages of these chondrites are determined by modeling the thermal histories of their parent bodies, driven by heating from the decay of ²⁶Al (Grimm & McSween 1993). These thermal profiles are compared with the formation ages and peak temperatures of secondary minerals such as fayalite (Fe₂SiO₄) and complex carbonates ([Ca,Mg,Mn]CO₃). These minerals serve as valuable indicators of aqueous alteration because they form through fluid-assisted reactions on the parent bodies, thereby providing a lower limit for the accretion age of the chondrites. Using the Mn-Cr dating method (Doyle et al. 2015), the formation age of fayalite in CV chondrites has been determined to be 4.2 ± 0.8 million years after the formation of CV Calcium-Aluminium-Inclusions (CAIs). These CAIs mark the reference time (t₀) for the first solids in the solar system, as established by the U-corrected Pb-Pb dating method (Connelly et al. 2012). Based on this lower limit, the modeled accretion age of the CV chondrite parent body is approximately 2.5 million years after CAI formation. This time gap between solar system formation and the accretion of CV parent bodies aligns with the timescales over which particles can be retained in vortices before dissipation. We therefore propose that vortices acted as storage reservoirs, accounting for both the preservation of early solar system materials (CV CAIs) and the delayed accretion of CV chondrules into their parent bodies.

7. Conclusions

We studied the 3D dynamics of dust in turbulent protoplanetary disks, examining how particles are trapped and retained within vortices self-sustained by the VSI. We performed 3D simulations using Lagrangian particles of three different sizes (1 mm, 500 μm, and 100 μm) to assess how vortices capture particles of varying sizes. Our main findings are summarized below:

1. We showed that vortices can trap dust particles in disks dominated by VSI-driven turbulence. Our simulations show the formation of multiple vortices, including small-scale structures likely induced by shearing azimuthal flows from the KH instability. Two large-scale vortices became fully developed and stable after approximately 250 orbits. The inner vortex efficiently captured most of the 1 mm particles, that underwent rapid radial drift and accumulated at the inner pressure maxima. Similarly, the 500 μm particles were trapped by both large-scale vortices. In contrast, the 100 μm particles were not captured by the inner vortex, as its efficient trapping of millimeter-sized particles occurs more rapidly, whereas the smaller grains are more dispersed in the surrounding area. This may explain their absence in NIR scattered light observations. However, these smaller particles are more effectively trapped by the outer vortex, which initially formed at the outer edge of the disk and slowly migrated inward; as a result, it lacks an associated gap in its vicinity.

2. We find that the VSI drives significant vertical diffusion, with our estimated α_z values for 1 mm, 500 μm, and 100 μm dust grains consistent with previous work. The ratio α_z/St ≈ 0.12 for 1 mm grains suggests a regime of efficient turbulent mixing. Inside the vortex, we found the α_z to be about three times lower than the globally determined α_z consistent with a low eddy timescale suggesting that vortices preferentially enhance diffusion in the radial and azimuthal directions rather than vertically.

3. We did not observe spiral structures associated with the two large-scale vortices in the gas velocity residuals or in the distribution of any particle species. Such features are more commonly linked to planet formation in a gap, where spiral wakes are generated in the disk. However, when analyzing the perturbed gas pressure, a trailing spiral-like feature became apparent at both the disk midplane and the disk surface.

4. By analyzing particle trajectories within vortices, we found that turbulent vertical motions driven by the VSI are reduced by up to a factor of two compared to regions where VSI turbulence is fully developed. This reduction implies that particle collision velocities remain low enough to prevent fragmentation, enabling growth through sticking and coagulation. Consequently, particles may reach sizes conducive to gravitational collapse and the formation of planetary embryos, consistent with Carrera et al. (2025), who demonstrated that pebbles drifting into vortices experience reduced turbulence via their drag force on the gas, thereby facilitating further growth. Future simulations should incorporate dust growth processes and dust backreaction to further assess the longevity of vortices and their role in planetesimal formation.

5. At high disk inclinations, dust crescents are affected by projection effects, that distort their asymmetry and reduce the brightness contrast along the dusty ring. This is important to consider, as many disks observed with ALMA have inclinations greater than 30°, making it challenging to identify dust crescents formed by vortices. Therefore, more than the 13 disks with observed dust crescents (Bae et al. 2023) may in fact host dust vortices. Whether the associated RWI vortices are formed by the VSI or by planetary companions, the presence of spirals observed in NIR scattered light in some systems (e.g., MWC 758 and HD 135344B) suggests that planetary perturbations may also play a role.

In this work, we have demonstrated that vortices may serve as initial sites of dust trapping. These initial sites of dust trapping by vortices could eventually lead to the formation of the first generation of planetesimals and planetary embryos in protoplanetary disks. Recent high-resolution studies (Lesur et al. 2025; Shariff & Umurhan 2024) have further highlighted the sensitivity of VSI-driven structures to numerical resolution, suggesting that small-scale turbulence, rather than coherent vortices, may dominate under effectively isothermal conditions. This emphasizes the need to explore the nonlinear regime of the VSI across a broader range of resolutions and cooling timescales. Future work should therefore include significantly higher-resolution simulations to determine whether the apparent absence of vortices at high resolution reflects a genuine physical effect or a limitation of spatial resolution.

Acknowledgments

We thanks the referee for the helpful discussions that improved the quality of the manuscript. L. Flores-Rivera and M. Lambrechts acknowledge this work is funded by the European Research Council (ERC Starting Grant 101041466-EXODOSS). A.J. is funded by the Carlsberg Foundation (Semper Ardens: Advance grant FIRSTATMO). We also thank Marcelo Barraza, for the helpful discussion about the velocity residuals. To David Melon Fuksman for helpful discussions on generating and interpreting the kinetic-energy spectrum. To Elishevah van Kooten and Jean Bollard for their input regarding to the application to cosmochemistry.

References

Appendix A: Vortex criterion

The Rossby Wave Instability (RWI) was first derived by Lovelace et al. (1999) and later expanded upon by Li et al. (2000, 2001) as a global instability that arises at local extrema in pressure or vortensity. It operates in regions where the function ℒ exhibits an extremum. This critical function, which serves as the criterion for the onset of the RWI, is defined as

$$\begin{array}{c}\hfill \mathcal{L}=\frac{{\mathrm{\Sigma }}_{\mathrm{gas}}}{{\omega }_z}{\left(\frac{P}{{\mathrm{\Sigma }}_{\mathrm{gas}}^{\gamma }}\right)}^{2/\gamma }\end{array}$$(A.1)

where Σ_gas is the surface gas density, ω_z is the vertical component of the vorticity, P is the vertically integrated gas pressure, and γ = 1.4 is the adiabatic index. We compute ℒ using azimuthally and vertically averaged profiles of density, pressure, and velocity and normalized it with their initial profiles. Figure A.1 (left) shows the RWI criterion at at 350 and 463 orbits. The ℒ/ℒ₀ curve displays multiple global extrema, indicating that large regions of the disk are susceptible to the RWI. Two prominent extrema are identified, and they align radially with the locations of the two large-scale vortices (indicated by vertical dashed gray lines).

Fig. A.1. Rossby instability criteria at the extrema in the disk.

The Rossby number, Ro, can also be used to identify regions exhibiting anticyclonic rotation. It quantifies the ratio of Coriolis to inertial forces in the Navier–Stokes equations and is defined as $\mathrm{Ro}={(\mathrm{\nabla }\times \mathit{v}-{\mathit{v}}_k)}_z/{\mathrm{\Omega }}_k$, where the velocity field is measured relative to Keplerian motion and averaged vertically. Figure A.1 (right) presents the spatial distribution of the Rossby number, highlighting azimuthally elongated vortical structures with negative Rossby numbers. These localized minima signal the presence of anticyclonic vortices, where the dynamics are dominated by a balance between pressure and Coriolis forces.

Appendix B: Vorticity residual for a different model

We also present results from a different simulation, with the same physical disk parameters as presented in the main text, but a different cooling time, $t_{\mathrm{cool}}=5\times {10}^{-3}{\mathrm{\Omega }}_k^{-1}$ (Flores-Rivera et al., in prep). Applying the same methodology used to construct the vorticity residual over time as in Figure 3, Figure B.1 reveals multiple vortex trails at this higher cooling rate. We observe multiple vortices originating between 40 au and 60 au, far from the boundary conditions; this suggests that such features form spontaneously in VSI-dominated disks simulated at the same resolution as in the present case. A caveat remains that at higher resolutions, with less numerical diffusion, the elliptical instability may inhibit vortex formation (Lesur & Papaloizou 2009; Lesur et al. 2025).

Fig. B.1. Time evolution of the minimum vorticity for a model with identical physical disk parameters but a modified cooling time, $t_{\mathrm{cool}}=5\times {10}^{-3}{\mathrm{\Omega }}_k^{-1}$.

Appendix C: Stokes number in the disk

The Stokes number at the midplane (Z/H_R = 0) is St = 1.4×10⁻² for a = 0.1 cm at 50 au. At one scale height (Z/H_R = 1), the Stokes number increases to St = 2.3×10⁻² for the same particle size. This implies that a one-scale-height variation leads to a change in Stokes number by less than a factor of two. Similarly, the Stokes number for particles of size 0.05 cm and 0.01 cm at 50 au is 7.1×10⁻³ and 1.4×10⁻³ at the midplane, respectively. This is further supported when considering a turbulent fragmentation velocity of 1 m s⁻¹ and α = 3 × 10⁻⁴ (adopted from our VSI turbulence convergence value, see Fig. 1 top), which yields a fragmentation-limited Stokes number of 1.8×10⁻³, corresponding to a particle size of ≃0.02 cm at the midplane at 50 au (see Eq. 33 in Birnstiel 2024).

These fragmentation-limited Stokes number values apply to the outer vortex region that is centered at ∼48.5 au, where the estimated α-viscosity ranges from ∼7 × 10⁻⁵ from the midplane up to 2 H. For the inner vortex, the α-viscosity is similar (3–7 × 10⁻⁴ over the same height range), but differences in scale height and sound speed affect the Stokes number. We find that the fragmentation-limited Stokes number in the vortices is 8 × 10⁻³ for α_{vortex, outer} = 7 × 10⁻⁵ for the outer vortex and, for the inner vortex it is one order of magnitude lower using α_{vortex, inner} = 7 × 10⁻⁴.

Fig. C.1. Stokes number in the disk. The different colored contour lines in each panel represents different Stokes number at different heights in the disk.

Appendix D: Pressure fluctuations around vortices

We examined the perturbed gas pressure in the disk, defined as $\delta P_{\mathrm{gas}}=\frac{P_{\mathrm{gas}}}{P_{\mathrm{gas},0}}$, where P_gas, 0 is the initial pressure. As shown in the top panels of Figure D.1, the strongest gas pressure perturbations occur at the locations of the inner and outer vortices. The perturbed gas pressure map reveals an arc-like feature extending toward the lower right quadrant of the disk, seemingly originating from the outer vortex in the upper left, suggestive of a trailing spiral or filamentary structure driven by vortex-induced disturbances. The radial profile, shown in the bottom panel of Figure D.1, displays steep pressure gradients at the radial positions of both vortices, while the root-mean-square (rms) velocity, v_rms (plotted in red), exhibits local minima at the same locations. Similar steep pressure gradients were reported in Stadler et al. (2025) for the sources that exhibit dust crescents in the submillimeter continuum emission.

Fig. D.1. Perturbed gas pressure structure in the midplane in logarithmic scale (top left) and in the surface (top right) of the global domain. Bottom panel shows the radial maximum gas pressure (blue curve) and the rms of the vertical velocity (red curve), vrms. At the location of the inner vortex and outer vortex we observe the pressure bumps following by a reduced vrms demonstrating that vortices are weak turbulence features in the disk.

Appendix E: Dynamical trajectory of the particle

We investigated the dynamical behavior of a particle once it is trapped by a vortex. A randomly selected particle, identified as a 1 mm-sized grain and a 100 μm-sized grain, was used to trace its 3D trajectory, as shown in Figure E.1. The scatter markers, connected by light gray line segments, indicate the particle’s position at each orbit. The apparent jumps in the trajectory arise from the azimuthal angle wrapping from π to −π as the vortex completes full revolutions around the disk. In the θ–r plane, both particles follow looping paths within the vortex, while the vortex itself gradually migrates inward.

Fig. E.1. Dynamical trajectory of a 100 μm-sized (left) and a 1 mm-sized (right) particle inside the outer vortex.

All Tables

All Figures

Fig. 1. Reynolds stress-to-pressure ratio, αr, ϕ. The top panel shows the time-averaged αr, ϕ, where the yellow line represents the convergence value, ⟨αr, ϕ⟩ = 3 × 10−4. For reference, the vertical dotted line corresponds to the outer local orbit at 0.25 R0 used by Lesur et al. (2025) in their analysis, as shown in their Figure 1. The bottom-left panel displays the time series of αr, ϕ in the r-direction, averaged over azimuth in the 3D VSI-active disk. The bottom-right panel presents the radial Reynolds stress-to-pressure ratio over the time series and averaged in the azimuthal directions.

In the text

	Fig. 2. Vorticity of the gas in the midplane in units of the keplerian frequency at 350 orbits. Overplotted are two dust-to-gas mass ratio values of 1 mm particles in the midplane.
In the text

	Fig. 3. Time evolution of the minimum vorticity residual as a function of radius, illustrating the radial position and migration of vortices in the VSI-unstable disk. The horizontal dashed lines outline the final locations of the two large-scale vortices.
In the text

Fig. 4. Pertubational azimuthal velocity (top and bottom-left panels) and vertical velocity (top and bottom middle panels) of the disk at the midplane (z = 0, top panels) and close to the disk surface (z > 2H, bottom panels). To make it easier to identify the two large-scale vortices in the disk, we show the dust distribution vertically averaged in the top right panel. The black arrows point at the location of the inner and outer vortices.

In the text

	Fig. 5. Velocity profile of a single 100 μm-sized particle (left) and a 1 mm-sized particle (right) within the vortex.
In the text

Fig. 6. Dust density distribution. Top panels: Dust density distribution vertically averaged at 350 orbits in the ϕ − R plane. The large grains are 100 μm (blue), 500 μm (green), and 1 mm (orange). The small grains spanning 0.1–30 μm are plotted as a black shade in the disk and are assumed to be well coupled with the gas. Bottom panels: Same dust density distribution as the top panels but viewed in the Z − R plane and averaged in the azimuthal direction.

In the text

	Fig. 7. Dust-to-gas mass ratio at the midplane at 350 orbits. The dotted lines represent the initial dust-to-gas mass ratio value for each grain size.
In the text

	Fig. 8. Mean vertical dust distribution, Z, averaged over the radial range 40–50 au and over time after 300 orbits. Colored lines show the mean normalized distributions for 1 mm (orange), 500 μm (green), and 100 μm (blue) grain sizes. The standard deviation of each distribution corresponds to the dust scale height (Hd) for the respective grain size.
In the text

Fig. 9. Dust scale height over time for all particle species. The dust scale height is radially averaged from 40 au to 50 au. Each marker represents the averaged dust scale height calculated by taking the rms of the vertical position of the particles for all sizes per orbit. Strong vertical fluctuations are observed for the 1 mm and 500 μm grain sizes after approximately 424 orbits. This corresponds to a regime where the dust distribution is characterized by localized and compact concentrations within vortices (see the orange region in the bottom three panels at three different snapshots).

In the text

Fig. 10. Synthetic images of the simulated disk at 0° inclination. The top row shows the brightness of the simulated disk after radiative transfer for four different ALMA wavelength bands (columns). The bottom row shows the convolution of the disk images using a circular beam with a FWHM of 12 au (σbeam = 5 au). Overall, the presence of two asymmetric dust features are clearly seen at the inner ring and in the outer parts of the dust continuum emission. The inner one was found to be at the same location as the dust concentration.

In the text

	Fig. 11. Dust crescent at different disk inclinations. The black contour represents 85% of the peak intensity for the ALMA band 6 frequency. The elongation of the dust crescent increases with disk inclination due to an imaging projection effect. This effect is also evident in the other ALMA wavelength bands.
In the text

	Fig. 12. Velocity residuals of the simulated disks at different disk inclinations. The same beam size used in Figures 10 and 11 was applied to the convolved images (bottom panels).
In the text

Fig. 13. Energy spectrum for a small box extracted from the simulation domain (blue curve) and for the full domain (orange curve) at 400 orbits. Overplotted are two power laws, k−5/3 and k−5 illustrating the Kolmogorov regime and a steeper dissipation regime, respectively. The dotted and dashed vertical gray lines mark the wavenumbers corresponding to the marginal and fastest-growing VSI modes, based on the prescription from Umurhan et al. (2016). The dash-dotted vertical line denotes to the mode whose wavelength equals the local disk scale height λ = H, corresponding to kH = 2π.

In the text

	Fig. A.1. Rossby instability criteria at the extrema in the disk.
In the text

	Fig. B.1. Time evolution of the minimum vorticity for a model with identical physical disk parameters but a modified cooling time, $t_{\mathrm{cool}}=5\times {10}^{-3}{\mathrm{\Omega }}_k^{-1}$.
In the text

	Fig. C.1. Stokes number in the disk. The different colored contour lines in each panel represents different Stokes number at different heights in the disk.
In the text

Fig. D.1. Perturbed gas pressure structure in the midplane in logarithmic scale (top left) and in the surface (top right) of the global domain. Bottom panel shows the radial maximum gas pressure (blue curve) and the rms of the vertical velocity (red curve), vrms. At the location of the inner vortex and outer vortex we observe the pressure bumps following by a reduced vrms demonstrating that vortices are weak turbulence features in the disk.

In the text

	Fig. E.1. Dynamical trajectory of a 100 μm-sized (left) and a 1 mm-sized (right) particle inside the outer vortex.
In the text