Reduction of dust radial drift by turbulence in protoplanetary disks

Context. Dust particles in protoplanetary disks rotate at velocities exceeding those of the surrounding gas due to a lack of pressure support. Consequently, they experience a headwind from the gas that drives them toward the central star. Radial drift occurs on timescales much shorter than those inferred from disk observations or those required for dust to aggregate and form planets. Additionally, turbulence is often assumed to amplify the radial drift of dust in planet-forming disks when modeled through an e ﬀ ective viscous transport. However, the local interactions between turbulent eddies and particles are known to be signiﬁcantly more intricate than in a viscous ﬂuid. Aims. Our objective is to elucidate and characterize the dynamic e ﬀ ects of Keplerian turbulence on the mean radial and azimuthal velocities of dust particles. Methods. We employed 2D shearing-box incompressible simulations of the gas, which is maintained in a developed turbulent state while rotating at a sub-Keplerian speed. Dust is modeled as Lagrangian particles set at a Keplerian velocity, therefore experiencing a radial force toward the star through drag. Results. Turbulent eddies are found to reduce the radial drift, while simultaneously enhancing the azimuthal velocities of small particles. This dynamic behavior arises from the modiﬁcation of dust


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 slowing 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 evidences suggest that disk turbulence exhibits lower amplitude than previously believed (Flaherty et al. 2015(Flaherty et al. , 2018;;Villenave et al. 2022).This is consistent with the saturated non-linear 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 non-negligible compared to the strength of turbulence, resulting in a reshape of the tur-bulent 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 clustering pf 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 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 particles radial drift, while concurrently enhancing their azimuthal drift.

Gas
We assume 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 consider the gas flow to be two-dimensional.We consider a small region located at distance r 0 from the star, rotating at a rate Ω along with the gas, and we employ the periodic shearing box approximation.If (r, θ) are the polar coordinates centered onto the star, we define the local coordinates x = r−r 0 and y = r 0 θ.The total gas velocity solves: 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 u = v + 3 2 Ω x e y is assumed periodic in the frame distorted by the mean shear.We conduct numerical simulations with 256 × 1024 collocation points using the open-source pseudo-spectral solver Snoopy (see Lesur & Longaretti 2005).
The turbulent flow is maintained in a statistically steady state by introducing both a forcing term, f , and a large-scale dissipation D = γ (v + 3 2 Ω x e y ).The prescribed force is random, Gaussian, homogeneous, isotropic in the sheared frame, with a zero mean, white noise in time, and with spatial correlations concentrated at large scales.This implementation allows us to investigate a generalized form of 2D turbulence.Although the resulting flow exhibits 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 τ ω = ⟨ω 2 ⟩ −1/2 , where ω = ∂ x u y − ∂ y u 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 have selected values of Ro in the range of [0.1, 10] (see Appendix A for details on the simulation parameters).Notice that at small Ro, rotation becomes predominant over turbulence.The effects of such weak turbulence on particle dynamics have been understudied until now.

Dust
The position and velocity (R p , U p ) of a particle in the rotating and sheared frame solve 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 ρ p a/(ρ g c 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 ϵ = (Ω 2 K − Ω 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 prescribing the level of gas turbulence, we can define the dimensionless drift parameter ε = ϵτ 2 f /ℓ f .We adopt the fiducial value ε = 0.1, to examine the effects of weak but non-negligible 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 −2 , 10].We neglect both particle-particle interactions and their feedback onto the gas.

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, 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 NT p,x = −τ p ϵ/(1+τ 2 p Ω 2 ) (as detailed in Appendix E).We therefore define the relative velocity ∆U p,x /U NT p,x = (⟨U p,x ⟩ − U NT p,x )/U NT p,x , where ⟨•⟩ denotes the average over all particles and times.This quantity measures the variation of 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 −1 and small Stokes numbers.This indicates that, for these parameters, turbulence substantially slows down dust radial drift.It is noteworthy that, for Ro −1 = 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 here needed to achieve sufficiently high statistics.

Radial drift
The mean radial drift velocity of particles is presented in Fig. velocities are smaller than in the laminar case.For Ro −1 > 0.7, dust lies longer in anticyclones, where it gets locked to the gas radial velocity 1 .Turbulent decrease of dust drift reaches up to an order of magnitude in some cases (i.e., circles for which ∆U p,x /U NT p,x = −90% in Fig. 1).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 −1 = 0.49 particles mostly reside between eddies, yet their radial drift is still slowed down by turbulence.For Ro −1 = 0.10, instead, small particles drift faster toward the star, as usually assumed for strong turbulence.We will 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 NT p,x , indicating large solids weakly coupled to the turbulent flow.

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 NT p,y = τ 2 p Ωϵ/2(1 + τ 2 p Ω 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 −1 ).It is important to notice that the magnitude of the strongest azimuthal drift velocity represents only a small fraction of the Keplerian velocity (⟨U p,y ⟩ ≈ 10 −4 r 0 Ω for Ro −1 = 0.10 and St = 0.1).Therefore, it would not be currently 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.

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.In Fig. 2, particularly for Ro −1 = 0.10, one can notice that 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 ones 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 (Fig. 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.

Reduction of radial dust drift
Figure 2 shows that, for Ro −1 > 0.10, the radial velocity of dust particles is decreased by turbulence.Two different regimes should be considered here.For Ro −1 > 0.7 particles mostly remain trapped inside anticyclones.Therefore, the reduction of 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 −1 = 0.49, where particles are instead ejected from eddies, counter-intuitive 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 of this radial drift reduction lies in the Coriolis force, which induces an interchange of azimuthal and radial velocities: The decrease of radial drift velocities by turbulence can therefore be easily explained.While preferential sweeping would normally enhance the radial velocity of dust, it also increase 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 of dust radial drift velocities.Moreover, as rotation rate augments, shear renders turbulence increasingly non-isotropic (Gerosa et al. 2023).In particular, eddies become strongly stretched azimuthally, as can be seen when comparing Fig. 4b (Ro −1 = 0.10) and Fig. 4d (Ro −1 = 0.49).This, in turn, affects the dynamics of dust.To better understand this mechanism, in Fig. 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.

Estimation of the turbulent parameter α
Turbulent intensity is classically defined as the ratio between the amplitude of the turbulent fluctuations u of the gas and 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: (5) The modifications of both dust velocities, due to preferential sweeping, and particle concentrations can be estimated through this value.Considering H = 10 l f , the α parameter can be as large as 10 −2 for Ro −1 = 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 −1 = 6.19 we obtain α = 10 −5 .At these low turbulence levels, not extensively studied in the literature, turbulent anticyclones are concentrating dust and slowing down drift.Indeed, 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 −1 = ∞), 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| 2 ⟩).
For this parameter we find values of two orders of magnitude lower than 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.

Dependence on the drift parameter ε
The shearing box approximation holds true only when the box is located at a distance r 0 from the star much larger than the largest scales of the simulation.This requires r 0 ≫ l f .Consequently, we obtain the condition ϵ ≫ (Ω 2 K − Ω 2 ) l f .The acceptable values of (Ω 2 K − Ω 2 ) are determined from: where n is the index of the power law describing pressure as a function of radius.In the literature, typical disk properties yield 0.996 < v 2 y /v 2 K < 0.999.For our analysis, we make the arbi-trary choice v 2 y /v 2 K = 0.998.Therefore, for slow rotation rates, ε can be as small as 10 −4 , while faster rotations necessitate considering larger drift parameter (e.g., εmin ≈ 0.2 for Ro −1 = 10).In Fig. 5 the black hatched region shows the inaccessible values of ε for the shearing box approximation as a function of the inverse Rossby number.Choosing a larger (smaller) value for v 2 y /v 2 K shifts the inaccessible region to smaller (larger) values of the drift parameter ε.In our exploration of the parameter space, we investigated various values of ε within the accessible range.Figure 5 presents the results on dust concentration and radial velocity variation as a function of ε and Ro −1 for St = 0.1.It demonstrates that, for ε ≤ 1, dust is ejected from eddies, and its velocity is enhanced at small Ro −1 .Conversely, dust tends to cluster more readily and experiences a decrease in radial drift for Ro −1 > 0.7.The results at ε = 1 further corroborate the findings presented in this paper, as they exhibit the same behavior as those at ε = 0.1, while being much further away from the inaccessible region (see Appendix F).However, at ε = 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 of dust radial drift is quantitatively smaller for the same Ro and St at larger ε.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.

Conclusions
We addressed in this Letter 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 in overcoming the radial drift barrier and 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, as numerous studies have indicated a twodimensionalization 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, we did not consider in this study 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.

Fig. 1 .
Fig. 1.Variation of dust radial drift velocity due to turbulence for ε = 0.1 as a function of the inverse Rossby number and the Stokes number.

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

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

Fig. 4 .
Fig. 4. Vorticity ω, normalized by τ ω , and dust particle position R p (black dots) for St ≈ 0.1.The snapshots in a) and c) have been obtained for a simple cellular flow, while b) and d) from our direct numerical simulations at ε = 0.1 and Ro −1 = 0.10 for b), Ro −1 = 0.49 for d).

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