© The Authors 2024

1. Introduction

Radial drift stands out as a pivotal process in the evolution of dust within planet-forming disks, playing a crucial role in planet formation. Dust drift arises from the sub-Keplerian rotation of the gas, exerting drag forces that slow down the inherently Keplerian motion of dust. The resulting depletion of dust angular momentum triggers its radial displacement toward the central star on a timescale significantly shorter than the disk’s lifetime. This process is seemingly in contradiction with the observed presence of dust in much older disks. Moreover, planet formation is expected to occur over timescales on the order of a million years, and needs a large amount of dust left in the disk to form the observed massive objects. However, this timescale is significantly longer than the radial drift timescale, leading to an inconsistency in planet formation theories, known as the radial drift barrier (Weidenschilling 1977; Nakagawa et al. 1986). Consequently, planetesimal formation is believed to be localized in specific regions where radial drift is arrested, and dust concentrates.

Recent observational evidence suggests that disk turbulence exhibits lower amplitude than previously believed (Flaherty et al. 2015, 2018; Villenave et al. 2022). This is consistent with the saturated nonlinear state obtained from various hydrodynamical instabilities, such as the vertical shear instability (Arlt & Urpin 2004; Nelson et al. 2013) or the convective overstability (Klahr & Hubbard 2014; Lyra 2014), rather than situations displaying strong turbulence (e.g., through the magnetorotational instability, Balbus & Hawley 1991). A weak turbulence modifies our understanding of turbulent particle concentration. In particular, rotation and shear become nonnegligible compared to the strength of turbulence, resulting in a reshape of the turbulent eddies. In this context, turbulence cannot be treated as homogeneous and isotropic. Our previous work (Gerosa et al. 2023) demonstrates that rotation, through the Coriolis force, can counterbalance the ejection of particles from turbulent eddies. This phenomenon fosters the clustering of dust and potential planetesimal formation within anticyclonic eddies, in contrast with earlier studies that focused on lower rotation rates (Cuzzi et al. 2001; Pan et al. 2011). Furthermore, it shows that Keplerian turbulence cannot be simplistically approximated as a purely diffusive process acting on particles.

Turbulence is frequently invoked for its ability to transport gas angular momentum, often modeled as a viscosity (Shakura & Sunyaev 1973; Lesur 2021). This paradigm leads to gas radial motion, consequently amplifying dust radial drift through drag in an accreting disk (Takeuchi & Lin 2002). However, such a mean-field approach neglects the local modifications that turbulent eddies can impose on particle trajectories. The nuanced understanding of non-isotropic turbulence and its interplay with dust dynamics, as elucidated in Gerosa et al. (2023), prompts an investigation into whether Keplerian turbulence additionally alters dust drift. Turbulent gas could potentially slow down or even halt the radial drift of dust, thus overcoming the need for an azimuthally extended pressure bump, which is the conventional solution to this barrier (Whipple et al. 1972; Pinilla et al. 2012).

The aim of this study is to investigate the impact of gas turbulence on the radial and azimuthal drift of dust. We employ 2D numerical simulations of forced Keplerian turbulence within a shearing box and track the dynamics of drifting Lagrangian particles in this flow, characterized by the absence of mean radial velocity. Our findings reveal how turbulence plays a crucial role in diminishing the radial drift of particles, while concurrently enhancing their azimuthal drift.

2. Model and numerical methods

2.1. Gas

We assumed that the gas is incompressible; thus, its dynamics are governed by a divergence-free velocity field v that satisfies the Navier–Stokes equation, with the inclusion of the external gravitational potential from the star. Additionally, we considered the gas flow to be two-dimensional. We considered a small region located at distance r₀ from the star, rotating at a rate Ω along with the gas, and we employed the periodic shearing box approximation. If (r, θ) are the polar coordinates centered onto the star, we defined the local coordinates x = r − r₀ and y = r₀θ. The total gas velocity solves as

$$\begin{array}{c}\hfill {\partial }_t\mathit{v}+\mathit{v}·\mathbf{\nabla }\mathit{v}=-\frac{1}{{\rho }_{\mathrm{g}}}\mathbf{\nabla }p+\nu {\mathbf{\nabla }}^2\mathit{v}-2\mathbf{\Omega }\times \mathit{v}+3{\mathrm{\Omega }}^{2}x{\mathit{e}}_x+\mathit{f}-\mathit{D},\end{array}$$(1)

where ρ_g and ν represent the gas mass density and kinematic viscosity, respectively. The computational domain has size L_x × L_y, with an aspect ratio of L_y/L_x = 4 to limit spurious geometrical effects at large Ω. The fluctuating velocity of the gas $\mathit{u}=\mathit{v}+\frac{3}{2}\mathrm{\Omega }x{\mathit{e}}_{\mathit{y}}$ is assumed periodic in the frame distorted by the mean shear. We conducted numerical simulations with 256 × 1024 collocation points using the open-source pseudo-spectral solver SNOOPY (see Lesur & Longaretti 2005).

The turbulent flow was maintained in a statistically steady state by introducing both a forcing term, f, and a large-scale dissipation $\mathit{D}=\gamma (\mathit{v}+\frac{3}{2}\mathrm{\Omega }x{\mathit{e}}_{\mathit{y}})$. The prescribed force is random, Gaussian, homogeneous, and isotropic in the sheared frame, with a zero mean, white noise in time, and with spatial correlations concentrated at large scales. This implementation allowed us to investigate a generalized form of 2D turbulence. Although the resulting flow exhibited some similarities with gas subjected to the subcritical baroclinic instability (Lesur & Papaloizou 2010), investigating the transition to a turbulent state through instabilities is beyond the scope of this paper. The forcing has a fixed amplitude, so that the mean injection rates of kinetic energy ε_I and enstrophy η_I are prescribed. This specifies a forcing length-scale ℓ_f = (ε_I/η_I)^1/2 and timescale τ_f = (η_I)^−1/3, which remain constant across all runs. We can also define the small eddy turnover time τ_ω = ⟨ω²⟩^−1/2, where ω = ∂_xu_y − ∂_yu_x is the vorticity of the gas turbulent fluctuations. Using this timescale, we nondimensionalize the rotation rate Ω, defining the Rossby number Ro = 1/(τ_ωΩ). It serves as a measure of turbulence strength relative to rotation. For our study, we selected values of Ro in the range of [0.1, 10] (see Appendix A for details on the simulation parameters). We note that at small Ro, rotation becomes predominant over turbulence. The effects of such weak turbulence on particle dynamics have been understudied until now.

2.2. Dust

The position and velocity (R_p, U_p) of a particle in the rotating and sheared frame solve as

$$\begin{array}{cc}& \frac{\mathrm{d}{\mathit{R}}_{\mathrm{p}}}{\mathrm{d}t}={\mathit{U}}_{\mathrm{p}}+\frac{3}{2}\mathrm{\Omega }{R}_{\mathrm{p},x}{\mathit{e}}_y,\hfill \end{array}$$(2)

$$\begin{array}{cc}& \frac{\mathrm{d}{\mathit{U}}_{\mathrm{p}}}{\mathrm{d}t}=-\frac{1}{{\tau }_{\mathrm{p}}}[{\mathit{U}}_{\mathrm{p}}-\mathit{u}({\mathit{R}}_{\mathrm{p}},t)]-2\mathbf{\Omega }\times {\mathit{U}}_{\mathrm{p}}+\frac{3}{2}\mathrm{\Omega }{U}_{\mathrm{p},x}{\mathit{e}}_y-ϵ{\mathit{e}}_x.\hfill \end{array}$$(3)

The equation governing the motion of dust particles is derived in Appendix B. The first term on the right-hand side of Eq. (3) comes from the drag with the gas. It involves the particle stopping time τ_p defined in Epstein regime as τ_p = (8/π)^1/2ρ_pa/(ρ_gc_s), with ρ_p denoting the particle material mass density, a its radius, and c_s the speed of sound. We nondimensionalize τ_p by defining the Stokes number St = τ_p Ω. The subsequent term appearing in Eq. (3) accounts for the Coriolis force. The last term represents the adjustment for the difference in azimuthal velocity between gas and dust in the shearing box framework. It involves the parameter $ϵ=({\mathrm{\Omega }}_{\text{K}}^2-{\mathrm{\Omega }}^2)r_0$, where Ω_K denotes the particle Keplerian rotation rate. When scaled relative to the two other fixed parameters quantities of the simulation, namely τ_f and ℓ_f, which prescribe the level of gas turbulence, we can define the dimensionless drift parameter $\stackrel{\sim }{ϵ}=ϵ{\tau }_{\mathrm{f}}^2/{\ell }_{\mathrm{f}}$. We adopt the fiducial value $\stackrel{\sim }{ϵ}=0.1$ to examine the effects of weak but nonnegligible turbulence on dust drift. Nonetheless, we extensively explore the influence of this parameter on our results in the discussion, considering values in the range [10⁻², 10]. We neglect both particle-particle interactions and their feedback onto the gas.

3. Results

An overview of the results is depicted in Fig. 1. The strength of rotation with respect to turbulence is increased along the vertical axis, while the horizontal axis represents the particle Stokes number, varying from small to large grains. The orange background zone highlights the region where dust forms point clusters within anticyclones, whereas in the gray background area dust fills the whole space. In the other regions of the parameter space, dust particles tend to concentrate on filaments. For a more detailed discussion of distinct dust concentration regimes, we refer to Gerosa et al. (2023) or to the highlights given in Appendix C. In the lower region, indicative of strong turbulence and/or slow rotation, dust particles tend to be expelled from eddies. Conversely, in the upper region, dust predominantly resides within anticyclones. The transition between these two behaviors can be quantitatively characterized (see Appendix D), and is identified in Fig. 1 by the dashed line. To determine whether turbulence diminishes the radial drift of dust, we compare its mean speed in our turbulent runs to that in the laminar case. In the absence of turbulence, the particle radial drift velocity is given by $U_{\mathrm{p},x}^{\mathrm{NT}}=-{\tau }_{\mathrm{p}}ϵ/(1+{\tau }_{\mathrm{p}}^2{\mathrm{\Omega }}^2)$ (as detailed in Appendix E). We therefore define the relative velocity $\mathrm{\Delta }{U}_{\mathrm{p},x}/U_{\mathrm{p},x}^{\mathrm{NT}}=(⟨U_{\mathrm{p},x}⟩-U_{\mathrm{p},x}^{\mathrm{NT}})/U_{\mathrm{p},x}^{\mathrm{NT}}$, where ⟨ ⋅ ⟩ denotes the average over all particles and times. This quantity measures the variation in dust radial drift velocity attributable to turbulence. A negative (positive) value indicates a decrease (increase) in the radial drift of dust compared to the laminar case. It is observed in Fig. 1 that the radial drift velocity is significantly reduced at large Ro⁻¹ and small Stokes numbers. This indicates that, for these parameters, turbulence substantially slows down dust radial drift. It is noteworthy that, for Ro⁻¹ = 3.29, the plot is the result of an average over 12 runs. Because of particles clustering to a point in phase space, many runs are needed here to achieve sufficiently high statistics.

Fig. 1. Variation in dust radial drift velocity due to turbulence for $\stackrel{\sim }{ϵ}=0.1$ as a function of the inverse Rossby number and the Stokes number.

3.1. Radial drift

The mean radial drift velocity of particles is presented in Fig. 2. The dashed line corresponds to the nonturbulent case, $U_{\mathrm{p},x}^{\mathrm{NT}}$. For small Stokes numbers, the radial drift velocity exhibits a substantial reduction for almost all the Rossby numbers, diminishing with faster rotation rates (larger Ro⁻¹). This can be clearly identified on this figure as the absolute value of the mean radial velocities are smaller than in the laminar case. For Ro⁻¹ > 0.7, dust lies longer in anticyclones, where it gets locked to the gas radial velocity¹. The turbulent decrease in dust drift reaches up to an order of magnitude in some cases (i.e., circles for which $\mathrm{\Delta }{U}_{\mathrm{p},x}/U_{\mathrm{p},x}^{\mathrm{NT}}=-90\%$ in Fig. 1). The yellow circles in Fig. 1, on the other hand, correspond to a complete cessation of radial drift due to turbulence, with particles clustered to a point and trapped inside an anticyclone. For Ro⁻¹ = 0.49 particles mostly reside between eddies, yet their radial drift is still slowed down by turbulence. For Ro⁻¹ = 0.10, instead, small particles drift faster toward the star, as usually assumed for strong turbulence. We delve into the interpretation of these last two phenomena in Sect. 4. Finally, for all Rossby numbers, ⟨U_p, x⟩ converges at large Stokes numbers to $U_{\mathrm{p},x}^{\mathrm{NT}}$, indicating large solids weakly coupled to the turbulent flow.

Fig. 2. Mean radial drift velocity as a function of the Stokes number, for various values of the Rossby number and for $\stackrel{\sim }{ϵ}=0.1$. The dotted line represents the radial drift velocity particles would have in a laminar gas flow.

3.2. Azimuthal drift

Figure 3 shows the mean azimuthal drift velocity ⟨U_p, y⟩ of dust particles. The laminar drift velocity in the azimuthal direction is given by $U_{\mathrm{p},\mathit{y}}^{\mathrm{NT}}={\tau }_{\mathrm{p}}^2\mathrm{\Omega }ϵ/2(1+{\tau }_{\mathrm{p}}^2{\mathrm{\Omega }}^2)$ (see Appendix E), plotted as a dashed line. Turbulence enhances the azimuthal drift of small dust particles, while for large St, it remains unaltered. The values of ⟨U_p, y⟩ clearly appear larger at slower rotation rates (small Ro⁻¹). It is important to note that the magnitude of the strongest azimuthal drift velocity represents only a small fraction of the Keplerian velocity (⟨U_p, y⟩≈10⁻⁴r₀Ω for Ro⁻¹ = 0.10 and St = 0.1). Therefore, it would not currently be detectable from observational data. Conversely, the azimuthal drift is reduced at very high rotation rates and intermediate Stokes numbers. This is again consistent with the mechanism of dust trapping inside anticyclones, resulting in particle cluster velocities governed by the mean flow.

Fig. 3. Mean azimuthal drift velocity as a function of the Stokes number, for various values of the Rossby number and for $\stackrel{\sim }{ϵ}=0.1$. The dotted line represents the azimuthal drift velocity particles would have in a laminar gas flow.

4. Discussion

4.1. Preferential sweeping

Preferential sweeping is a well-known process that induces variations in particle speed, due to turbulence in the gas phase. This phenomenon has been extensively studied and demonstrated in the context of droplets or aerosols settling in the atmosphere (Maxey 1987; see also Bec et al. 2024 for a recent review). In such scenarios the interaction between gravity and turbulence leads to an acceleration of particle settling as particles tend to preferentially sample down-welling regions of the turbulent flow. We anticipate a similar mechanism for dust particles moving within turbulent protoplanetary disks, where gravity drives radial drift of the dust toward the central star. Figure 2 shows, particularly for Ro⁻¹ = 0.10, that the radial drift is increased by turbulence for intermediate Stokes numbers. A simple analysis using cellular eddies, as depicted in Fig. 4a, further elucidates this point. Cyclonic eddies (rotating in the direction of the disk) are shown in red, while anticyclonic eddies are shown in blue. It is evident that dust particles exclusively sample left-welling regions, attracted by the star’s gravity, resulting in increased radial drift velocities. When inspecting snapshots from our numerical simulations (Figs. 4b and d), the dynamics of dust appears to be more intricate. Here, particles predominantly sample left-welling regions, while also undergoing azimuthal drift. During this process, they are once again subject to preferential sweeping. Consequently, they experience an additional acceleration in the azimuthal direction, as seen for small enough Stokes numbers in Fig. 3.

Fig. 4. Vorticity ω, normalized by τω, and dust particle position Rp (black dots) for St ≈ 0.1. The snapshots in a) and c) were obtained for a simple cellular flow, while from our direct numerical simulations at $\stackrel{\sim }{ϵ}=0.1$ and Ro−1 = 0.10 for b) and Ro−1 = 0.49 for d).

4.2. Reduction of radial dust drift

Figure 2 shows that, for Ro⁻¹ > 0.10, the radial velocity of dust particles is decreased by turbulence. Two different regimes should be considered here. For Ro⁻¹ > 0.7, particles mostly remain trapped inside anticyclones. Therefore, the reduction in their drift can be readily explained: particles evolve with the radial velocity of eddies, which is slower than the dust velocity. On the other hand, at Ro⁻¹ = 0.49, where particles are instead ejected from eddies, counterintuitively in regard to the effect of preferential sweeping on dust speed, our findings are in contrast with the typical enhancement of dust speed caused by standard preferential sweeping. The reason for this radial drift reduction lies in the Coriolis force, which induces an interchange of azimuthal and radial velocities:

$$\begin{array}{c}\hfill f_{\mathrm{C}}=-2\mathbf{\Omega }\times {\mathit{U}}_{\mathrm{p}}=2\mathrm{\Omega }{U}_{\mathrm{p},y}{\mathit{e}}_x-2\mathrm{\Omega }{U}_{\mathrm{p},x}{\mathit{e}}_y.\end{array}$$(4)

The decrease in radial drift velocities by turbulence can therefore be easily explained. While preferential sweeping would normally enhance the radial velocity of dust, it also increases U_p, y that, when entering the Coriolis term, can counteract the acceleration of radial drift. This effect can even prevail, particularly at high rotation rates, resulting in an effective decrease in dust radial drift velocities. Moreover, as the rotation rate increases, the shear renders the turbulence increasingly non-isotropic (Gerosa et al. 2023). In particular, eddies become strongly stretched azimuthally, as can be seen when comparing Fig. 4b (Ro⁻¹ = 0.10) and Fig. 4d (Ro⁻¹ = 0.49). This, in turn, affects the dynamics of dust. To better understand this mechanism, in Figs. 4a and c, we illustrate dust dynamics in a model flow consisting of cellular eddies with two different aspect ratios to represent vortex stretching due to shear. In both cases, particles follow the preferential sweeping scenario, solely sampling left-welling regions. However, with elongated eddies, dust paths are predominantly azimuthally oriented, thus resulting in a significant reduction in particle mean radial motion.

4.3. Estimation of the turbulent parameter α

Turbulent intensity is classically defined as the ratio of the amplitude of the turbulent fluctuations u of the gas to its mean velocity. This quantity is the most relevant for estimating the effects of turbulence on particle dynamics. From its squared value, we can compute the turbulent α parameter as

$$\begin{array}{c}\hfill \alpha =\frac{⟨|\mathit{u}{|}^2⟩}{{(\mathrm{\Omega }H)}^2}.\end{array}$$(5)

The modifications of dust velocities, due to preferential sweeping, and of particle concentrations can be estimated through this value. Considering H = 10 l_f, the α parameter can be as large as 10⁻² for Ro⁻¹ = 0.10. At this rotation rate, turbulence is ejecting particles from eddies and accelerating their drift. Its effect on particle dynamics is therefore similar to what is already known in the literature for strong MRI-like turbulence (Yang et al. 2018). On the other hand, at Ro⁻¹ = 6.19 we obtain α = 10⁻⁵. At these low turbulence levels, not extensively studied in the literature, turbulent anticyclones concentrate dust and slow down drift. In this case, the Coriolis force due to rotation is able to counteract the diffusive power of eddies. We therefore observe a substantial difference in dust dynamics in weak turbulence compared to the purely laminar case (Ro⁻¹ = ∞), where particles would be randomly distributed and would drift radially at a considerably faster velocity.

When instead considering the angular momentum transport due to turbulence (Shakura & Sunyaev 1973), we can compute a viscous α using the Reynolds stress ⟨u_x u_y⟩ (in place of ⟨|u|²⟩). For this parameter we find values that are two orders of magnitude lower than those given above. This indicates a small amount of turbulent momentum transport in our simulations, in line with recent observations that tend to infer low turbulent viscosity in disks.

4.4. Dependence on the drift parameter $\stackrel{\sim }{ϵ}$

The shearing box approximation holds true only when the box is located at a distance r₀ from the star much greater than the largest scales of the simulation. This requires r₀ ≫ l_f. Consequently, we obtain the condition $ϵ\gg ({\mathrm{\Omega }}_K^2-{\mathrm{\Omega }}^2)l_{\text{f}}$. The acceptable values of $({\mathrm{\Omega }}_K^2-{\mathrm{\Omega }}^2)$ are determined from

$$\begin{array}{c}\hfill v_y^2/v_{\mathrm{K}}^2=1-n{c}_{\mathrm{s}}^2/v_{\mathrm{K}}^2,\end{array}$$(6)

where n is the index of the power law describing pressure as a function of radius. In the literature the typical disk properties yield $0.996<v_y^2/v_{\text{K}}^2<0.999$. For our analysis we make the arbitrary choice $v_y^2/v_{\text{K}}^2=0.998$. Therefore, for slow rotation rates, $\stackrel{\sim }{ϵ}$ can be as small as 10⁻⁴, while faster rotations necessitate considering a larger drift parameter (e.g., ${\stackrel{\sim }{ϵ}}_{\min }\approx 0.2$ for Ro⁻¹ = 10). In Fig. 5 the black hatched region shows the inaccessible values of $\stackrel{\sim }{ϵ}$ for the shearing box approximation as a function of the inverse Rossby number. Choosing a larger (smaller) value for $v_y^2/v_{\text{K}}^2$ shifts the inaccessible region to smaller (larger) values of the drift parameter $\stackrel{\sim }{ϵ}$.

Fig. 5. Phase diagram of the drift parameter vs. the inverse Rossby number. The highlighted parameters present the results of the dust concentration and radial drift for St = 0.1.

In our exploration of the parameter space, we investigated various values of $\stackrel{\sim }{ϵ}$ within the accessible range. Figure 5 presents the results of dust concentration and radial velocity variation as a function of $\stackrel{\sim }{ϵ}$ and Ro⁻¹ for St = 0.1. It demonstrates that, for $\stackrel{\sim }{ϵ}\le 1$, dust is ejected from eddies, and its velocity is enhanced at small Ro⁻¹. Conversely, dust tends to cluster more readily and experiences a decrease in radial drift for Ro⁻¹ > 0.7. The results at $\stackrel{\sim }{ϵ}=1$ further corroborate the findings presented in this paper, as they exhibit the same behavior as those at $\stackrel{\sim }{ϵ}=0.1$, while being much farther away from the inaccessible region (see Appendix F). However, at $\stackrel{\sim }{ϵ}=10$, for which turbulence becomes a secondary mechanism compared to drift, a higher rotation rate is required for the slowdown of dust and particle concentration to happen. It is important to note that the decrease in dust radial drift is quantitatively smaller for the same Ro and St at larger $\stackrel{\sim }{ϵ}$. Therefore, a small value of the drift parameter is essential for the reduction of radial drift due to turbulence to significantly impact the long-term evolution of dust in disks.

5. Conclusions

In this letter we addressed the issue of dust particle drift in a turbulent gas, usually considered as a problem in theories of planetesimal formation. We demonstrated that turbulent eddies can reduce the radial drift velocity of small dust particles through two distinct mechanisms: preferential sweeping of particles between elongated eddies at slow rotation rates and concentration within anticyclones for faster rotations. In the latter scenario, dust particles may even come to a complete halt in their radial drift. Concurrently, turbulence tends to enhance the azimuthal drift of dust in most cases, although azimuthal velocities of particles are reduced when they aggregate inside anticyclonic eddies.

The reduction of dust radial drift has significant implications, potentially reconciling the masses of dust disks inferred from observations with the predicted mass loss rates from theory, particularly when turbulence prolongs the timescale of dust drift by an order of magnitude. Furthermore, our findings suggest that a turbulent region in the disk could act as a traffic jam, potentially explaining the formation of certain observable substructures, such as those detected by telescopes like ALMA. Slowing down radial drift may also facilitate planetesimal formation by assisting dust particles to overcome the radial drift barrier and by reducing collision speeds, both crucial factors for a favorable outcome of their interactions.

In our simulations we made several assumptions. The shearing box approach is relevant for our local numerical simulations, focusing specifically on dust particles dynamics. The turbulence intensity and considered scales are also small enough for density perturbations to be weak in the box (Méheut et al. 2015), therefore supporting the incompressible approximation. The choice of two-dimensional simulations is justified for small Rossby numbers, since numerous studies have indicated a two-dimensionalization of the flow at high rotation rates (Yeung & Zhou 1998; Biferale et al. 2016). Additionally, for large Ro, the 2D approximation remains valid due to the highly stratified nature of protoplanetary disks (Cambon 2001). Finally, in this study we did not consider the back-reaction of dust on gas or self-interactions between particles. While these effects can be significant when dust particles are concentrated, they are often considered secondary aspects when diluted. However, dust feedback can have strong consequences on both particle clustering (Johansen & Youdin 2007) and dust radial drift (Dipierro et al. 2018), topics that are left for future work.

Acknowledgments

This work was supported by the “Programme National de Physique Stellaire” (PNPS) and “Programme national de planétologie” (PNP) of CNRS/INSU co-funded by CEA and CNES, and Observatoire de la Côte d’Azur. This work was supported by the French government, through the UCA^JEDI Investments in the Future project managed by the National Research Agency (ANR) under reference number ANR-15-IDEX-01. The authors are grateful to the OPAL infrastructure and the Université Côte d’Azur’s Center for High-Performance Computing for providing resources and support.

References

Appendix A: Run parameters

Table A.1 presents the parameters used in the direct numerical simulations with the spectral code Snoopy.

We conducted seven distinct runs, each time varying the rotation rate Ω, which resulted in specific values for the eddy turnover time τ_ω (as listed in Table A.2), and, consequently, for the Rossby number.

Finally, each of these runs is seeded with 48 sets of 10⁴ Lagrangian particles. These families are associated with 4 different values of ϵ (0.1, 1, 10, 100) and 12 values of τ_p (0.008, 0.014, 0.025, 0.044, 0.080, 0.144, 0.260, 0.470, 0.849, 1.533, 2.769, and 5.000).

Appendix B: Shearing box and drift

We consider dust particles in Keplerian rotation, thus rotating with a rotation rate Ω_K faster than the gas, which rotates with the shearing box at a sub-Keplerian speed Ω. Considering the box centered at a distance r₀ from the star, the equation for passive particles is

$$\begin{array}{cc}\hfill \frac{d{\mathit{V}}_{\mathrm{p}}}{\mathit{dt}}=& -\frac{1}{{\tau }_{\mathrm{p}}}[{\mathit{V}}_{\mathrm{p}}-\mathit{v}({\mathit{R}}_{\mathrm{p}},t)]-2\mathbf{\Omega }\times {\mathit{V}}_{\mathrm{p}}+(2{\mathrm{\Omega }}_K^2+{\mathrm{\Omega }}^2)R_{\mathrm{p},x}{\mathit{e}}_x\hfill \\ \hfill & -({\mathrm{\Omega }}_K^2-{\mathrm{\Omega }}^2)r_0{\mathit{e}}_x.\hfill \end{array}$$(B.1)

For numerical reasons, it is easier to compute the particle velocity in the frame where the background negative shear has been subtracted. Therefore, we perform the change of variables U_p, y = V_p, y + (3/2)Ωx_p giving

$$\begin{array}{cc}\hfill \frac{d{U}_{\mathrm{p},x}}{\mathit{dt}}=& -\frac{1}{{\tau }_{\mathrm{p}}}[U_{\mathrm{p},x}-u_x({\mathit{R}}_{\mathrm{p}},t)]+2\mathrm{\Omega }{U}_{\mathrm{p},y}+(2{\mathrm{\Omega }}_K^2+{\mathrm{\Omega }}^2)R_{\mathrm{p},x}\hfill \\ \hfill & -({\mathrm{\Omega }}_K^2-{\mathrm{\Omega }}^2)r_0-3{\mathrm{\Omega }}^2{R}_{\mathrm{p},x},\hfill \\ \hfill \frac{d{U}_{\mathrm{p},y}}{\mathit{dt}}=& -\frac{1}{{\tau }_{\mathrm{p}}}[U_{\mathrm{p},y}-u_y({\mathit{R}}_{\mathrm{p}},t)]-\frac{1}{2}\mathrm{\Omega }{U}_{\mathrm{p},x}.\hfill \end{array}$$(B.2)

Finally, if R_p, x ≪ r₀, we obtain $2({\mathrm{\Omega }}_K^2-{\mathrm{\Omega }}^2)x_{\text{p}}\ll ({\mathrm{\Omega }}_K^2-{\mathrm{\Omega }}^2)r_0$, so that

$$\begin{array}{c}\hfill \frac{d{\mathit{U}}_{\mathrm{p}}}{\mathit{dt}}=-\frac{1}{{\tau }_{\mathrm{p}}}[{\mathit{U}}_{\mathrm{p}}-\mathit{u}({\mathit{R}}_{\mathrm{p}},t)]-2\mathbf{\Omega }\times {\mathit{U}}_{\mathrm{p}}+\frac{3}{2}\mathrm{\Omega }{U}_{\mathrm{p},x}{\mathit{e}}_y-ϵ{\mathit{e}}_x,\end{array}$$(B.3)

with $ϵ=({\mathrm{\Omega }}_K^2-{\mathrm{\Omega }}^2)r_0$.

Appendix C: Dust clustering

The clustering properties of particles in a turbulent Keplerian flow can vary significantly depending on the flow and dust parameters. In a previous paper where no drift was considered (Gerosa et al. 2023), we identified three possible turbulent concentration paths:

1. Ejection from eddies and concentration on filamentary structures in between, for small rotation rates and dust sizes.

2. Distribution in spirals inside the anticyclones, but avoiding their core, for intermediate rotation rates.

3. Formation of point clusters in the cores of anticyclones, for large rotation rates and small dust sizes.

When the drift is considered, the picture remains qualitatively the same. As an example, we show in Fig. C.1 the concentration of dust particles over time, for $\stackrel{\sim }{ϵ}=0.1$, St = 1, and Ro⁻¹ = 3.29. Due to the Coriolis force, dust initially gathers in anticyclonic eddies (Fig. C.1a). Over time, each of these clouds coalesces into a point cluster in the core of an anticyclone (Fig. C.1b). Clusters eventually merge to form a single cluster (Fig. C.1c). In the extreme case of $\stackrel{\sim }{ϵ}=10$, however, the process of dust clustering is halted by drift, which carries particles away from anticyclones on a timescale too fast for dust to concentrate efficiently within them.

Fig. C.1. Vorticity ω, normalized by τω, and dust particle position Rp (black dots) for $\stackrel{\sim }{ϵ}=0.1$,St = 1, Ro−1 = 3.29 and a) t = 200τω, b) t = 2000τω, c) t = 20000τω. The clusters in b) and c) are highlighted with pink circles.

Appendix D: Okubo-Weiss parameter

The Okubo-Weiss parameter estimates the properties of a gas flow, helping to identify its local structures. It is defined as

$$\begin{array}{c}\hfill OW=-{({\partial }_x{u}_x)}^2-{\partial }_x{u}_y{\partial }_y{u}_x.\end{array}$$(D.1)

When positive, it describes rotating elliptic structures of the flow (e.g., eddies). For straining hyperbolic regions, the Okubo-Weiss parameter instead assumes negative values. This parameter has been measured along the particles path and its averaged value, as a function of Ro⁻¹ for various Stokes number, is shown in Fig. D.1. We note that ⟨OW(R_p)⟩ assumes negative values at small Ro⁻¹, while it become positive at high rotation rates. This means that, independently of the Stokes number, particles in a slow rotating gas flow predominately lie between eddies. Instead, for Ro⁻¹ ≳ 0.7, dust resides longer inside anticyclonic eddies.

Fig. D.1. Okubo-Weiss parameter computed at the particle position as a function of the inverse Rossby number, for various Stokes numbers and $\stackrel{\sim }{ϵ}=0.1$

Appendix E: Mean drift velocity

An expression for the mean particle velocity can be obtained by averaging Eq. (B.3) assuming statistical stationarity of turbulent fluctuations along particle paths. This leads to the following expressions for the mean radial and azimuthal particle velocities:

$$\begin{array}{cc}\hfill ⟨U_{\mathrm{p},x}⟩& =\frac{⟨u_x({\mathit{R}}_{\mathrm{p}},t)⟩+2{\tau }_{\mathrm{p}}\mathrm{\Omega }⟨u_y({\mathit{R}}_{\mathrm{p}},t)⟩-{\tau }_{\mathrm{p}}ϵ}{1+{\tau }_{\mathrm{p}}^2{\mathrm{\Omega }}^2},\hfill \\ \hfill ⟨U_{\mathrm{p},y}⟩& =-\frac{{\tau }_{\mathrm{p}}\mathrm{\Omega }⟨u_x({\mathit{R}}_{\mathrm{p}},t)⟩+⟨u_y({\mathit{R}}_{\mathrm{p}},t)⟩-{\tau }_{\mathrm{p}}^2\mathrm{\Omega }ϵ}{2(1+{\tau }_{\mathrm{p}}^2{\mathrm{\Omega }}^2)}.\hfill \end{array}$$(E.1)

In the absence of gas turbulence (u = 0), these equations yield the mean particle velocity $U_{\mathrm{p}}^{\mathrm{NT}}$ in a laminar flow.

It should be noted that without drift (ϵ = 0), even if the gas flow is turbulent, the average gas velocity at particle position ⟨u(R_p, t)⟩ would vanish, leading to ⟨U_p⟩ = 0. On the other hand, when drift is present, as in our case, the mean gas velocity at the particle position does not vanish, highlighting the phenomenon of preferential sweeping. This mechanism indeed arises from a preferential sampling of specific zones and structures of the flow along particle paths, modifying the mean radial and azimuthal drifts of the particles.

Appendix F: Radial drift at $\stackrel{\sim }{ϵ}=1$

Figure F.1 shows the relative increase in dust radial drift due to turbulence $\mathrm{\Delta }{U}_{\mathrm{p},x}/U_{\mathrm{p},x}^{\mathrm{NT}}=(⟨U_{\mathrm{p},x}⟩-U_{\mathrm{p},x}^{\mathrm{NT}})/U_{\mathrm{p},x}^{\mathrm{NT}}$ for $\stackrel{\sim }{ϵ}=1$. We note that the results closely resemble those presented in Fig. 1 for $\stackrel{\sim }{ϵ}=0.1$. A subtle reduction of radial drift can be noted at Ro⁻¹ = 6.19 and small St compared to Ro⁻¹ = 3.29, indicating a reversal in trend. In this case, turbulence is too weak to efficiently counteract dust drift.

Fig. F.1. Variation in dust radial drift velocity due to turbulence for $\stackrel{\sim }{ϵ}=1$ as a function of the inverse Rossby number and the Stokes number.

All Tables

All Figures

	Fig. 1. Variation in dust radial drift velocity due to turbulence for $\stackrel{\sim }{ϵ}=0.1$ as a function of the inverse Rossby number and the Stokes number.
In the text

	Fig. 2. Mean radial drift velocity as a function of the Stokes number, for various values of the Rossby number and for $\stackrel{\sim }{ϵ}=0.1$. The dotted line represents the radial drift velocity particles would have in a laminar gas flow.
In the text

	Fig. 3. Mean azimuthal drift velocity as a function of the Stokes number, for various values of the Rossby number and for $\stackrel{\sim }{ϵ}=0.1$. The dotted line represents the azimuthal drift velocity particles would have in a laminar gas flow.
In the text

	Fig. 4. Vorticity ω, normalized by τω, and dust particle position Rp (black dots) for St ≈ 0.1. The snapshots in a) and c) were obtained for a simple cellular flow, while from our direct numerical simulations at $\stackrel{\sim }{ϵ}=0.1$ and Ro−1 = 0.10 for b) and Ro−1 = 0.49 for d).
In the text

	Fig. 5. Phase diagram of the drift parameter vs. the inverse Rossby number. The highlighted parameters present the results of the dust concentration and radial drift for St = 0.1.
In the text

	Fig. C.1. Vorticity ω, normalized by τω, and dust particle position Rp (black dots) for $\stackrel{\sim }{ϵ}=0.1$,St = 1, Ro−1 = 3.29 and a) t = 200τω, b) t = 2000τω, c) t = 20000τω. The clusters in b) and c) are highlighted with pink circles.
In the text

	Fig. D.1. Okubo-Weiss parameter computed at the particle position as a function of the inverse Rossby number, for various Stokes numbers and $\stackrel{\sim }{ϵ}=0.1$
In the text

	Fig. F.1. Variation in dust radial drift velocity due to turbulence for $\stackrel{\sim }{ϵ}=1$ as a function of the inverse Rossby number and the Stokes number.
In the text