Homogeneous turbophoresis of heavy inertial particles in turbulent flow

Abstract

Heavy particles suspended in turbulent flow possess inertia and are ejected from violent vortical structures by centrifugal forces. Once piled up along particle paths, this small-scale mechanism leads to an effective large-scale drift. This phenomenon, known as ‘turbophoresis’, causes particles to leave highly turbulent regions and migrate towards calmer regions, explaining why particles transported by non-homogeneous flows tend to concentrate near the minima of turbulent kinetic energy. It is demonstrated here that turbophoretic effects are just as crucial in statistically homogeneous flows. Although the average turbulent activity is uniform, instantaneous spatial fluctuations are responsible for inertial-range inhomogeneities in the particle distribution. Direct numerical simulations are used to probe particle accelerations, specifically how they correlate to local turbulent activity, yielding an effective coarse-grained dynamics that accounts for particle detachment from the fluid and ejection from excited regions through a space- and time-dependent non-Fickian diffusion. This leads to cast fluctuations in particle distributions in terms of a scale-dependent Péclet number ${\textit {Pe}}_\ell$, which measures the importance of turbulent advection compared with inertial turbophoresis at a given scale $\ell$. Multifractal statistics of energy dissipation indicate that $ {\textit {Pe}}_\ell \sim \ell ^\delta /\tau _{p}$ with $\delta \approx 0.84$. Numerical simulations support this behaviour and emphasise the relevance of the turbophoretic Péclet number in characterising how particle distributions, including their radial distribution function, depends on $\ell$. This approach also explains the presence of voids with inertial-range sizes, and the fact that their volumes have a non-trivial distribution with a power-law tail $p(\mathcal {V}) \propto \mathcal {V}^{-\alpha }$, with an exponent $\alpha$ that tends to 2 as ${\textit {Pe}}_\ell \to 0$.

1. Introduction

The transport of small, heavy particles by a developed turbulent flow is a common occurrence in nature and industry. Whether they are droplets in air, dust in gas, or sediments in water, these particles are often smaller than the smallest active scale of the fluid and have a larger mass density. They thus possess inertia, resulting in their detachment from the carrier fluid and uneven spatial distributions, a phenomenon known as preferential concentration. This is important in determining the interactions between these particles, such as collisions and aggregation. It also alters the transfers of momentum, kinetic energy and heat in the particle-laden fluid. One notable example of inertial particles is water droplets in atmospheric clouds. As stressed by Jonas (1996), turbulence triggers variability in droplet sizes that can explain why the time scales for rain initiation are much shorter than those predicted by mean-field arguments. Pinsky & Khain (1997) (see also Shaw (2003)) demonstrated that the preferential concentration of droplets affects their growth by condensation and coalescence. Heterogeneities have been observed in situ (see e.g. Kostinski & Shaw 2001) and their small-scale effects have been quantified to improve droplet collision rates (see Reade & Collins 2000; Falkovich, Fouxon & Stepanov 2002). Still, many challenging questions raised in clouds involve interactions over a huge range of scales and thus, cannot be addressed without having recourse to large-eddy simulations (LES). Such approaches need ad hoc parameterisations of particle dynamics and their microphysical interactions, as discussed for instance in Morrison et al. (2020). Planet formation by dust aggregation in the early Solar system is another important natural instance of inertial particles, which raises similar issues. Local fluctuations in the particle concentration trigger gravitational collapse and thus the formation of larger objects. Because of rotation around the star, dust particles migrate in large-scale anticyclonic Keplerian vortices (Gerosa, Méheut & Bec 2023) or in pressure bumps (Johansen et al. 2007). It is probably in these regions that primary accretions occur, but the effect of turbulence is still unclear (Johansen et al. 2015). A better understanding requires developing models to quantify dust clustering in the inertial range of turbulence (see e.g. Hartlep & Cuzzi 2020) and designing LES tools that cope with astrophysical specificities. Other natural situations where inertial particles occur include plankton ecology in the ocean (Seuront, Schmitt & Lagadeuc 2001) and seed dispersion above plant canopies (Pan, Chamecki & Isard 2014). In all cases, a precise description of large-scale fluctuations in particle density is crucial.

Equivalent questions arise in engineering. When optimising droplet vaporisation in injection sprays (Sahu, Hardalupas & Taylor 2018) or monitoring particulate fouling (Henry, Minier & Lefèvre 2012), it is important to understand how inertial particles distribute over scales comparable to the larger scales of the carrier turbulent flow. The complexity of flow geometries and inhomogeneities in industrial applications give a critical role to the spatial variations of the time-averaged particle density. Much effort has thus been dedicated to derive effective transport equations for the average concentration field. In this context, Caporaloni et al. (1975) unveiled a fundamental mechanism in which turbulence inhomogeneities drive particles out of the most excited regions of the flow and concentrate them in quieter zones. They dubbed this phenomenon turbophoresis (see also Reeks (1983)), in analogy to thermophoresis, where temperature gradients cause a motion of diffusive particles towards colder regions of space. Reeks (1983, 1992) proposed closures of the kinetic equations for the particle phase-space distribution to derive effective diffusion equations for the average spatial concentration. This leads to the particle fluxes due to inertia being described by a Fick law, where the coefficient of diffusion is related to the local Lagrangian correlation of the fluid velocity. Such arguments have been successfully employed to explain why particles in turbulent channel flows tend to migrate towards the walls (see e.g. Marchioli & Soldati 2002; Kuerten & Vreman 2005; Sardina et al. 2012; Fouxon et al. 2018; Brandt & Coletti 2022). However, the dependence of the diffusion coefficient on the particle Stokes number is not yet fully understood. Belan, Fouxon & Falkovich (2014) (see also Belan (2016)) showed that particles with sufficient inertia escape from low-kinetic-energy regions, leading to a localisation–delocalisation phase transition. De Lillo et al. (2016) examined the case of turbulent flows with an inhomogeneous forcing and found that turbophoretic effects are more pronounced at intermediate particle inertia. Mitra, Haugen & Rogachevskii (2018) interpreted this behaviour as a balance between turbophoretic and turbulent diffusions.

The applicability of turbophoresis to particle transport in flows with average inhomogeneities raises questions about its relevance in homogeneous situations. In homogeneous isotropic turbulence, instantaneous snapshots reveal spatial fluctuations of kinetic energy throughout the inertial range. Meanwhile, particle distributions display heterogeneous concentrations characterised by large-scale quasiuniform regions, localised voids and sheet-like clusters, as observed for instance by Eaton & Fessler (1994). To quantify inertial-range particle distributions, different observables are needed compared with those used for the dissipative range. At small scales, particle distributions exhibit multifractal scaling properties (see Hogan, Cuzzi & Dobrovolskis 1999; Bec et al. 2011; Schmidt, Fouxon & Holzner 2017; Bec, Gustavsson & Mehlig 2024) and are fully characterised by a dimension spectrum that depends solely on the Stokes number. The unified picture of the joint dependence on length scale and response time arises from the fact that dissipative-range dynamics involve a unique time scale determined by the typical amplitude of velocity gradients. This is in contrast to the hierarchy of time scales involved in inertial-range physics. In the two-dimensional inverse cascade, Boffetta, De Lillo & Gamba (2004) found that particles concentrate quasiuniformly on thin filamentary structures separated by voids whose distribution follows a universal scaling law. However, in the random, white-in-time, self-similar flows considered by Bec, Cencini & Hillerbrand (2007b), such scaling is absent, and particle distributions are characterised by local fractal dimensions determined by the scale-dependent Stokes number $\textit {St}_{\ell } \propto \tau _{p}/\ell ^{2/3}$ (Balkovsky, Falkovich & Fouxon 2001), defined by non-dimensionalising the particle response time by the turnover time at the observation scale $\ell$. Both of these scenarios coexist in three-dimensional turbulence, as pointed up by Bec et al. (2007a), Yoshimoto & Goto (2007) or inferred from the sweep-stick mechanism of Goto & Vassilicos (2008). The intricate spatial correlations of the pressure gradient, or equivalently of the fluid acceleration, play a key role. By using Voronoï tessellations, Monchaux, Bourgoin & Cartellier (2010, 2012) introduced a definition of particle clusters and found that their size distribution follows a universal scaling law independent of the Stokes number. This was confirmed by Baker et al. (2017), who showed that clusters preferentially sample regions of the flow with higher strain and lower vorticity. However, Bragg, Ireland & Collins (2015) found that this statistical bias depends on inertia and is actually quantified by the scale-dependent Stokes number $\textit {St}_{\ell }$. These arguments led them to predict scale invariance for two-particle statistics when $\textit {St}_{\ell }\ll 1$, which was confirmed by Hartlep, Cuzzi & Weston (2017) using a cascade multiplier approach. Ariki et al. (2018) further argued that the pair correlation function follows a universal power-law $\propto \textit {St}_{\ell }^2$ using a Lagrangian renormalisation closure. The wavelet analysis conducted by Matsuda, Schneider & Yoshimatsu (2021) shows intermittent particle densities, with a stronger contribution from voids observed at smaller spatial scales. However, the question of whether scale invariance holds in the inertial-range distributions of particles and, if so, which mechanisms are involved, remains ambiguous.

To shed new light on these issues, an appropriate effective model for inertial-range particle dynamics is expected to be useful. While there are various simplified approaches to dilute particle suspensions, reviewed for instance by Balachandar & Eaton (2010), the Eulerian field representations of the particle phase proposed by Ferry & Balachandar (2001) provide promising tools. In this approach, the particle velocity is tied to the carrier phase, with the effect of inertia being regarded as a compressible correction proportional to the fluid velocity acceleration. While this approximation has often shown its relevance, it remains limited to the asymptotics of small particle inertia, and it combines very different time scales, as acceleration is influenced by dissipative-range physics. Fevrier, Simonin & Squires (2005) extended these considerations to large Stokes numbers by assuming that the particle motion can be seen as the sum of a mesoscopic velocity and a random component. The latter term corresponds to a diffusive motion, is uncorrelated in space, and has been found to properly reproduce particle properties when their response time is much larger than the turbulent large-eddy turnover time. This contribution, dubbed random uncorrelated motion by Reeks, Fabbro & Soldati (2006), was used by Gustavsson et al. (2012) in synthetic random flows and shown to suitably describe the effect of fold caustics on the particle kinetics. This approach relies on the idea that turbulence has only a cumulative effect along particle paths, as long as the latter have a sufficiently long correlation time. However, fluctuations do not need to be averaged over times prescribed by the particles’ lag, but this procedure can rather stem from a spatial or temporal coarse-graining of the turbulent field, thus incorporating the effect of instantaneous spatial inhomogeneities.

We aim here to introduce a model capable of effectively describing and quantifying particle dynamics within the inertial range of fully developed turbulent flows. We argue that small-scale detachments from the fluid can be elucidated by examining the accelerations experienced by the particles. Rather than simply filtering out fluctuations in a time-reversible manner, these detachments cause particles to carry forward past fluctuations. Building on the phenomenology introduced by Bec & Chétrite (2007), we argue that this mechanism cumulates over time, leading to an ejection process that causes non-Fickian particle fluxes. Our proposed model utilises an Itô, rather than Stratonovich, diffusion process with a diffusion coefficient that varies based on the local flow activity. This statistically homogeneous turbophoresis can be used to quantify inhomogeneities in the particle distribution and, to some extent, reconcile the various viewpoints discussed above. Our analysis is grounded in the results of direct numerical simulations conducted at large Reynolds numbers and relies on a comprehensive evaluation of particle accelerations.

The paper is structured as follows. In § 2 we introduce our settings and discuss the relevant observables for our analysis. We also provide a general appreciation of the correlations between particle concentrations and instantaneous inhomogeneities in turbulent activity. In § 3 we develop a Lagrangian perspective and conduct a detailed statistical analysis of particle acceleration to quantify particle detachment from the fluid across varying response times and observation scales. In § 4 we shift our focus to the Eulerian frame and use acceleration results to derive an effective equation for the particle coarse-grained density. From this model, we draw properties of the inertial-range distribution and discuss specifically the implications of this approach to the distribution of voids. Finally, in § 5 we summarise our findings and discuss possible perspectives.

2. Models, simulations and spatial coarse-graining

2.1. Homogeneous isotropic turbulence and energy dissipation

We investigate the behaviour of particles passively suspended in a three-dimensional fluid flow. The velocity field of the fluid, denoted by ${\boldsymbol {u}}(\boldsymbol {x},t)$, satisfies the incompressible Navier–Stokes equations

(2.1)\begin{equation} \partial_t {\boldsymbol{u}} + {\boldsymbol{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol{u}} ={-} (1/\rho_{f})\,\boldsymbol{\nabla} p + \nu \nabla^2 {\boldsymbol{u}} +\boldsymbol{f} \quad \text{with} \ \boldsymbol{\nabla}\boldsymbol{\cdot} {\boldsymbol{u}} = 0, \end{equation}

where $p$ represents the pressure, $\rho _{f}$ is the mass density of the fluid, $\nu$ is its kinematic viscosity and $\boldsymbol {f}$ is an external volume force. The force is prescribed with homogeneous and isotropic statistics and is correlated on large scales in both space and time. The force injects kinetic energy into the flow at an average rate of $\varepsilon = \langle \boldsymbol {f}\boldsymbol {\cdot }{\boldsymbol {u}}\rangle$. We perform direct numerical simulations of (2.1) using the pseudospectral code LaTu on the triply periodic box $[0,2 {\rm \pi}]^3$ and employ third-order Runge–Kutta time marching. The details of the code can be found in Homann, Dreher & Grauer (2007). Two sets of simulations are carried out with different resolutions. Corresponding numerical and turbulent parameters are presented in table 1.

Table 1. Parameters of the numerical simulations: $N^3$, number of collocation points; $\nu$, fluid kinematic viscosity; $\Delta t$, time step; $\varepsilon$, average kinetic-energy dissipation rate; $\eta =\nu ^{3/4}/\varepsilon ^{1/4}$, Kolmogorov dissipative scale; $u_{rms}$, root mean square velocity; $L=u_{rms}^3/\varepsilon$, large scale; $R_\lambda =\sqrt {15}\,u_{rms}^2/(\varepsilon \nu )^{1/2}$, Taylor-scale Reynolds number; $N_{p}$ number of particles for each value of the Stokes number.

After a certain period of time, the fluid velocity field ${\boldsymbol {u}}$ reaches a statistical steady state characterised by multifractal statistics of the local dissipation rate $\varepsilon _{loc}({\boldsymbol {x}}) = (\nu /2)\,\mathrm {tr}(\boldsymbol {\nabla }{\boldsymbol {u}}({\boldsymbol {x}})+\boldsymbol {\nabla }{\boldsymbol {u}}^{\mathsf {T}}({\boldsymbol {x}}))^2$ (see e.g. Frisch 1995). This is evidenced from the scale-dependent statistics of the coarse-grained dissipation $\varepsilon _\ell$ obtained by averaging the local dissipation over the ball $\mathcal {B}_{\ell }({\boldsymbol {x}})$ of centre ${\boldsymbol {x}}$ and diameter $\ell$, that is $\varepsilon _\ell ({\boldsymbol {x}}) \equiv (1/|\mathcal {B}_{\ell }|)\int _{\mathcal {B}_{\ell }({\boldsymbol {x}})}\varepsilon _{loc}({\boldsymbol {x}}')\,{\rm d}^3x'$. For $\eta \ll \ell \ll L$, the probability distribution of $\varepsilon _\ell$ takes the form

(2.2)\begin{equation} p(\varepsilon_\ell)\,\mathrm{d}\varepsilon_\ell = \left(\ell/L \right)^{3-\mathcal{F}(\alpha)}\mathrm{d}\mu(\alpha) \quad \mbox{with}\ \varepsilon_\ell = \varepsilon(\ell/L)^{\alpha-1}, \end{equation}

where $\mathcal {F}(\alpha )$ is the multifractal spectrum, which can be interpreted as the dimension of the fractal set on which the scale-averaged dissipation is $\propto \ell ^{\alpha -1}$ when $\ell /L\to 0$, and $\mathrm {d}\mu (\alpha )$ corresponds to the weight associated with each singularity exponent $\alpha$. In Kolmogorov 1941 phenomenology, there are no fluctuations of $\varepsilon _\ell$ and $\mathcal {F}(\alpha )=-\infty$ except for $\alpha =1$, for which $\mathcal {F}(1) = 3$. Figure 1(a) shows the multifractal spectrum obtained from numerical measurements of $\varepsilon _\ell$. We find that log-normal statistics, for which the dimension spectrum is a parabola $\mathcal {F}(\alpha ) = 3-(\alpha -1-\mu /2)^2/(2\mu )$, provide a good approximation for the lowest values of $\alpha$. It is, however, known that log-normal distributions have several shortcomings due to the non-conservative nature of the cascade models on which they are based (see discussion in Frisch (1995), § 8.6.5). Despite this, such an approximation is still useful for estimating moderate-order statistics, well beyond the central-limit approximation. Using Kolmogorov (1962) refined similarity hypothesis, and recent confirmations by Lawson et al. (2019), the statistics of the fluid velocity can be related to fluctuations in $\varepsilon _\ell$. For the longitudinal structure functions $S_n^\parallel (\ell ) = \langle [\hat {\boldsymbol {\ell }}\boldsymbol {\cdot } ({\boldsymbol {u}}({\boldsymbol {x}}+\boldsymbol {\ell })-{\boldsymbol {u}}({\boldsymbol {x}}))]^n\rangle$, the log-normal approximation with intermittency parameter $\mu$ predicts a scaling behaviour $S_n^\parallel (\ell ) \sim \ell ^{\zeta _n}$, where $\zeta _n = n/3+(\mu /18)(3n-n^2)$. In our simulations, we observe $\mu \approx 0.26$, consistent with the seminal work of Sreenivasan & Kailasnath (1993). This results in $\zeta _2\approx 0.696$, $\zeta _4\approx 1.276$, $\zeta _6\approx 1.740$, which are in good agreement with experimentally measured values (see Saw et al. (2018) for a recent review).

Figure 1. Statistics of the scale-averaged dissipation rate $\varepsilon _\ell$ for $R_\lambda = 460$. (a) Measured dimension spectrum $\mathcal {F}(\alpha )$ as a function of $\alpha = 1+\log (\varepsilon _\ell /\varepsilon ) / \log (\ell /L)$ at various scales $\ell$. The log-normal approximation for $\mu =0.26$ is shown as a dashed line. (b) Moments of order $p=2$, $3$, $4$ of the average dissipation rate $\varepsilon _\ell$ conditioned on its value $\varepsilon _{\ell '}$ in a larger box of size $\ell ^\prime =128\eta$, as a function of $\ell /\ell ^\prime$ and different conditionings as indicated in the colour bar. Log-normal behaviours are displayed as dashed lines.

Multifractal statistics are often interpreted phenomenologically as resulting from the random multiplicative cascade experienced by the coarse-grained dissipation. This scenario suggests that the probability distribution (2.2) should also apply to the fluctuations of $\varepsilon _{\ell }({\boldsymbol {x}})$ conditioned on the observed value of $\varepsilon _{\ell ^\prime }({\boldsymbol {x}})$ at the same location but over a larger scale $\ell ^\prime >\ell$. As shown in figure 1(b) for $\ell ^\prime = 128\eta$, numerical simulations confirm this feature, revealing a scaling regime with an exponent that is closely approximated by the log-normal prediction. These multiscale statistics play a crucial role in investigating the coarse-grained dynamics of transported particles, as we will discuss in more detail later on.

2.2. Particles, preferential sampling and concentrations

After the fluid flow reaches a statistical steady state, we introduce heavy, inertial, point-like particles that are homogeneously seeded with velocities equal to that of the fluid at their positions. The trajectories $\boldsymbol {x}_{p}(t)$ of these particles follow

(2.3a,b)\begin{equation} \frac{\mathrm{d}\kern0.7pt\boldsymbol{x}_{p}}{\mathrm{d}t} = \boldsymbol{v}_{p},\quad \frac{\mathrm{d}\boldsymbol{v}_{p}}{\mathrm{d}t} = \boldsymbol{a}_{p} ={-} \frac{1}{\tau_{p}}\left[\boldsymbol{v}_{p} - {\boldsymbol{u}}(\boldsymbol{x}_{p},t)\right]. \end{equation}

Particles are assumed much smaller than the Kolmogorov dissipative scale $\eta$, and sufficiently massive to neglect so added-mass, Magnus, and history effects. The viscous drag intensity is given by the response time $\tau _{p} = \rho _{p}\,d_{p}^2 / (18\nu \rho _{f})$, where $\rho _{p}$ is the particle mass density and $d_{p}$ its diameter. This time is used to define the Stokes number $\textit {St} = \tau _{p} / \tau _{\eta }$, with $\tau _{\eta } = (\nu / \epsilon )^{1/2}$ denoting a Kolmogorov dissipative time scale. The Stokes number measures particle inertia. When $\textit {St}\ll 1$, the particles almost follow the flow and behave as tracers. When $\textit {St}\gg 1$, they detach from the flow and behave ballistically. We adopt a Lagrangian approach in our simulations, where particles’ trajectories are tracked by integrating (2.3a,b) with the fluid velocity at their location obtained by linear interpolation from the grid. We use 10 different values of the Stokes number in the range $\textit {St} \in [0.1, 6.5]$ and, for each value of $\textit {St}$, a number $N_{p}$ of particles that roughly corresponds to one particle per box of size $(9\eta )^3$.

Upon reaching a statistically stationary state, the particle distributions exhibit highly non-uniform patterns and strongly correlate with the turbulent structures of the flow, as depicted in figure 2(a). The spatial arrangement of particles shows voids in the most active regions of the flow, where dissipation is high, sheet-like clusters that encapsulate these voids, and quasiuniform distributions in regions with lower turbulent intensity. These concentration fluctuations are attributed to the inertial-range motions of particles, as the sizes of the regions are much larger than the dissipative scale $\eta$. To filter out dissipative-range effects, we introduce the coarse-grained particle density $\langle \rho _{p} \rangle _{\ell }$. It is obtained by counting the number of particles in small boxes of size $\ell$, which define a partition of the spatial domain. Figure 2(b) shows $\langle \rho _{p} \rangle _{\ell }$ obtained with $\ell = 32\eta$. The spatial variations of particle dynamics also serve as a marker for the different regions of the flow. In figure 2(c), we show the coarse-grained root mean square acceleration obtained by averaging the squared modulus of acceleration for all particles located in given boxes of size $\ell$. Particle voids clearly correspond to high accelerations, indicating that concentration fluctuations are caused by detachment from the fluid and expulsion from active regions. It is worth noting that this mechanism differs somewhat from the conventional picture of inertial ejection by centrifugal forces from isolated vortices, as the thickness of vortex filaments is several times smaller than the coarse-graining scale $\ell$. This suggests that the observed particle ejections results from the collective effect of multiple turbulent structures.

Figure 2. (a) Two-dimensional slice of the instantaneous energy dissipation field $\varepsilon _{loc}({\boldsymbol {x}})$ for $R_\lambda =460$, together with particle positions for $\textit {St}=1$ shown as black dots. Colours stand for $\log _{10}(\varepsilon _{loc} /\varepsilon )$. (b) Particle coarse-grained density $\langle \rho _{p} \rangle _{\ell }$ in the upper-half of the same slice obtained for $\ell = 32\eta$. Colours are again on a logarithmic scale. (c) Root mean square particle acceleration, coarse-grained over the same grid, here normalised by $(\varepsilon ^3/\nu )^{1/4}$.

The observed correlations between particle concentrations and Lagrangian accelerations prompt a discussion on whether these findings can be explained by the sweep-stick mechanism, originally proposed by Goto & Vassilicos (2008) (see also Coleman & Vassilicos (2009)). In this scenario, particles ‘stick’ along the manifolds where fluid acceleration is orthogonal to the directions associated with the largest eigenvalues of its gradient. These particle-laden manifolds are then ‘swept’ by the fluid flow. The sweep-stick mechanism is most effective in describing particles’ concentrations for $\textit {St}\approx 1$ and on observation scales within the dissipative range of turbulence. Oka & Goto (2021) extended this approach by proposing that, for inertial-range response times ($\tau _\eta \ll \tau _{p} \ll \tau _L$), particles cluster near manifolds derived from a coarse-grained acceleration field at a scale $\ell$ chosen such that the associated turnover time $\tau _\ell = \varepsilon ^{-1/3}\ell ^{2/3}$ matches $\tau _{p}$. However, by focusing on a single resonant scale, this approach may overlook the multiscale nature of particle clustering, which requires accounting for turbulent fluctuations across a broader range of scales to fully capture the complex dynamics involved.

The importance of turbulent fluctuations and the tendency of particles to concentrate in regions of low turbulent activity can be further quantified by measuring their preferential sampling of energy dissipation within the flow. The mean value of $\varepsilon _\ell$ computed along the paths of particles with different Stokes numbers is shown as a function of the coarse-graining scale $\ell$ in figure 3(a). Particles sample preferentially regions where $\varepsilon _\ell$ is lower than the average dissipation $\varepsilon$, even when their inertia is weak, see figure 3(b). This bias persists in the inertial range, indicating that it stems from agitation accumulated along particle paths rather than instantaneous ejection from the flow small-scale structures. Measurements of the multifractal spectrum evaluated at particle positions confirm this tendency, as shown in figure 3(c) for $\ell =32\eta$. The dependence on $\textit {St}$ is weak and visible only at negative values of the singularity exponent corresponding to the most violent events. At $\alpha _{p}>1$, the dimension spectra associated with different Stokes numbers are almost undistinguishable. This suggests that preferential sampling results from the expulsion of particles from the most singular regions rather than convergence towards calmer ones.

Figure 3. (a) Mean coarse-grained dissipation $\varepsilon _\ell$ at particle location as a function of the averaging scale $\ell$ for various Stokes numbers and $R_\lambda = 460$. (b) Discrepancy shown this time as a function of $\textit {St}$ for $\ell =16$, $32$ and $64\eta$. Solid lines are fits $\propto \exp (-c(\ell )/\textit {St}^{1/2})$. (c) Dimension spectrum of the singularity exponent $\alpha _{p} =1+ \log (\varepsilon _\ell (\boldsymbol {x}_{p})/\langle \varepsilon _\ell (\boldsymbol {x}_{p}) \rangle ) / \log (\ell /L)$ at particle position for $\ell = 32\eta$, various $\textit {St}$, tracers, and the log-normal approximation.

The observed correlations between the dynamical and concentration properties of particles and the instantaneous inertial-range inhomogeneities of the turbulent flow suggest that the underlying mechanisms are akin to turbophoresis in non-homogeneous flows, at least qualitatively. Specifically, particles tend to move away from regions with high turbulent activity, forming voids and follow the fluid in calmer zones. To provide quantitative support for these ideas, we aim to develop effective equations for an averaged particle density. In the study of turbophoresis in non-homogeneous flows, these equations are obtained by averaging over either the realisations of turbulence or time in statistically stationary and ergodic situations. However, such classical averages are not applicable to instantaneous particle distribution in homogeneous turbulence. Nevertheless, we expect that a similar effective dynamics can be derived from a low-pass-filtered viewpoint, where the coarse-grained average $\langle {\cdot } \rangle _\ell$ plays a central role.

3. Non-homogeneous diffusion of Lagrangian trajectories

We revisit here the classical approach used to develop stochastic Langevin models for turbulent transport (see e.g. Minier (2016), for a review). The approach is based on the assumption that while Lagrangian velocities are correlated over time scales of the order of the integral time scale, acceleration becomes uncorrelated much faster, justifying an approximation of trajectories as diffusive processes. We begin in § 3.1 by providing effective approximations for the second-order statistics of fluid acceleration. We then extend these approximations to inertial particles in § 3.2, specifically to describe their spatially averaged acceleration. Finally, we use these results in § 3.3 to approximate particle dynamics as a diffusion process with a space- and time-dependent diffusion coefficient.

3.1. Fluid acceleration

Turbulent accelerations of fluid particles are among the most striking signatures of intermittency. At the turn of the century, significant advances in direct numerical simulations and in particle-tracking experimental techniques have enabled detailed investigations into acceleration statistics (see Toschi & Bodenschatz (2009), for a review). These studies revealed that the variance of acceleration deviates from its dimensional estimate and exhibits a notable dependence on the Reynolds number. Specifically, it can be expressed as $\langle |\boldsymbol {a}|^2 \rangle = A_2(\textit {Re})\,\varepsilon ^{3/2}/\nu ^{1/2}$, where $A_2$ accounts for this dependence, with $\textit {Re} = R_\lambda ^2/15 = \varepsilon ^{1/3}L^{4/3}/\nu = u_{rms}L/\nu$ denoting the large-scale Reynolds number. Hill (2002b) found that at moderate values of the Reynolds number, Taylor's scaling suggests $A_2 \propto \textit {Re}^{1/2}$, assuming that acceleration is dominated by pressure gradients. At large $\textit {Re}$, intermittency prevails and $A_2 \propto \textit {Re}^\gamma$, where $\gamma$ can be estimated using multifractal approaches (see e.g. Borgas 1993; Sawford et al. 2003; Biferale et al. 2004). Indeed, assuming that $|\boldsymbol {a}|^2 \sim \varepsilon _{\eta _\alpha }^{3/2}/\nu ^{1/2}$ with $\eta _\alpha$ such that $Re_{\eta _\alpha } = \varepsilon _{\eta _\alpha }^{1/3}\eta _\alpha ^{4/3}/\nu = 1$ and $\varepsilon _{\eta _\alpha } = \varepsilon (\eta _\alpha /L)^{\alpha -1}$, one can relate the fluctuation of acceleration to that of the singularity exponent $\alpha$. A saddle-point argument then yields $\gamma = \sup _\alpha [3(\mathcal {F}(\alpha )+3)/(\alpha +3)] -9/2$. Using the log-normal approximation of § 2.1, one gets $\gamma = -3[2 \mu +\sqrt {(\mu -32) \mu +64}-8]/(2 \mu ) \approx 0.078$ for $\mu = 0.26$. To match the two behaviours expected at moderate and large Reynolds numbers, we introduce the ad hoc approximation

(3.1)\begin{equation} A_2(\textit{Re}) \approx \frac{a\,\textit{Re}^{\gamma}}{[1+(R_\star{/}\textit{Re})^{1/2}]^{1-2\gamma}}. \end{equation}

Figure 4(a) compares this fit with numerical measurements by Gotoh & Fukayama (2001), Bec et al. (2006) and Yeung et al. (2006), together with current simulations. The approximation (3.1) with $\gamma =0.078$, $a = 6.2$ and $R_\star =80$ provides a reasonably good agreement.

Figure 4. (a) Normalised variance of the fluid acceleration as a function of the Reynolds number $\textit {Re} = R_\lambda ^2/15$. Data from several numerical studies are shown as symbols. The dashed line is a behaviour $\propto \textit {Re}^\gamma$ with $\gamma = 0.078$. The solid line corresponds to the fitting formula (3.1) with $a = 6.2$ and $R_\star =80$. (b) Spatial correlation of the fluid acceleration $C(r) = \langle \boldsymbol {a}({\boldsymbol {r}},t)\boldsymbol {\cdot }\boldsymbol {a}(0,t)\rangle$ with $r = |{\boldsymbol {r}}|$, shown for the two Reynolds numbers of our dataset. The dashed line is $\propto r^{\zeta _4-2}$ with $\zeta _4 = 1.276$. The dotted curve is the parabolic approximation $C(r) \approx C(0)\,[1-(r/\lambda _1)^2/2]$ with $\lambda _1 = [3\,C(0)/Q(0)]^{1/2} \approx 5.3\eta$. The solid line is the approximation (3.4) that displays a behaviour ${\simeq }C(0)\,\lambda _1/r$ at large separations.

We now examine the spatial correlations of acceleration, which will be important in approximating particle displacement later. In an isotropic flow, the correlation tensor components in the longitudinal and transverse directions to a given separation ${\boldsymbol {r}}$ are interrelated. In a homogeneous and isotropic flow, $C(r) = \langle \boldsymbol {a}({\boldsymbol {r}},t)\boldsymbol {\cdot }\boldsymbol {a}(0,t)\rangle$ depends only on the distance $r=|{\boldsymbol {r}}|$ and satisfies (Obukhov & Yaglom 1951; Hill & Wilczak 1995)

(3.2)\begin{equation} \frac{1}{r^2}\frac{\mathrm{d}}{\mathrm{d} r}\left(r^2 \frac{\mathrm{d}C}{\mathrm{d} r} \right) \approx{-}Q(r), \quad \mbox{with} \ Q(r) \equiv \left\langle \partial_ju_i({\boldsymbol{r}},t)\, \partial_iu_j({\boldsymbol{r}},t) \, \partial_lu_k(0,t) \, \partial_ku_l(0,t) \right\rangle, \end{equation}

where summation is assumed over repeated indices. This relation assumes that pressure gradients dominate acceleration and uses the Poisson equation to express them in terms of velocity gradients. In homogeneous isotropic flow, the right-hand side of (3.2) can be expressed as $Q(r) = (1/6)\,\partial _{ijkl} D_{ijkl}({\boldsymbol {r}})$, where $D_{ijkl}({\boldsymbol {r}}) = \langle [u_i({\boldsymbol {r}})-u_i(0)] [u_j({\boldsymbol {r}})-u_j(0)] [u_k({\boldsymbol {r}})-u_k(0)] [u_l({\boldsymbol {r}})-u_l(0)]\rangle$ is the fourth-order structure function. Since $D_{ijkl} \propto r^{\zeta _4}$, we expect acceleration correlations to decrease as $C(r) \propto r^{\zeta _4-2}$ when $\eta \ll r \ll L$. Therefore, we get $\propto r^{-0.724}$, which is steeper than the K41 prediction $\propto r^{-2/3}$ proposed by Obukhov & Yaglom (1951). Figure 4(b) shows the spatial correlations of acceleration for our two numerical simulations. Our data are consistent with the experimental measurements of Xu et al. (2007) and display a power law with an even steeper exponent close to $-1$.

This behaviour extends beyond the transition scale introduced by Hill (2002a), which is derived from the Taylor expansion of correlations at small separations. Equation (3.2) yields

(3.3)\begin{equation} C(r) = \frac{1}{r}\int_0^r r^{\prime2}\,Q(r^\prime)\,\mathrm{d}r^\prime + \int_r^\infty r^\prime\,Q(r^\prime)\,\mathrm{d}r^\prime. \end{equation}

Hence, the correlation function $C(r)$ can be approximated to leading order as $C(r) \approx C(0) [1- (1/2)\,(r/\lambda _1)^2]$ as $r$ approaches zero. Here, $\lambda _1 = [3\,C(0)/Q(0)]^{1/2}$ with $Q(0) = \langle [\mathrm {tr} (\boldsymbol {\nabla }{\boldsymbol {u}})^2 ]^2\rangle >0$ is the length scale characterising the parabolic decay of the acceleration correlations, analogous to the Taylor microscale for velocity correlations.

It is worth noting that both $C(0)$ and $Q(0)$ exhibit an intermittent dependence on the Reynolds number. While $C(0)/(\varepsilon ^{3/2}/\nu ^{1/2})$ scales as $\textit {Re}^{\gamma }$, the Reynolds-number dependence of $Q(0)$ relates to that of the fourth-order moment of velocity gradients. As found by Nelkin (1990), the $p$th order moments of fluid velocity derivatives can be expressed using multifractal formalism, similar to the variance of acceleration. Specifically, assuming that $|\boldsymbol {\nabla }{\boldsymbol {u}}| \sim (\varepsilon _{\eta _\alpha }/\nu )^{1/2}$, with the same definition as before of the fluctuating dissipative scale $\eta _\alpha$, one obtains $\langle |\boldsymbol {\nabla }{\boldsymbol {u}}|^p\rangle \sim (\varepsilon /\nu )^{p/2}\,\textit {Re}^{\chi _p}$, where $\chi _p = \sup _\alpha [3(p(1-\alpha )/2+3-\mathcal {F}(\alpha ))/(\alpha +3)]$. This leads to $Q(0)/(\varepsilon /\nu )^{2} \sim \textit {Re}^{\chi _4}$, where the log-normal approximation gives $\chi _4 = -3(3 \mu +\sqrt {(\mu -48) \mu +64}-8)/(2 \mu )$. In summary, we have $\lambda _1 \propto \eta \,\textit {Re}^{(\gamma -\chi _4)/2}$. Using $\mu =0.26$ for the intermittency parameter, we obtain $\gamma \approx 0.078$ and $\chi _4\approx 0.217$, which yields $\lambda _1 \propto \eta \,\textit {Re}^{-0.069}$, consistent with the prediction of Hill (2002a). Our numerical simulations reveal that $\lambda _1/\eta = 5.51$ for $R_\lambda =290$ and $\lambda _1/\eta = 5.32$ for $R_\lambda =460$, confirming a weak dependence on Reynolds number. Nevertheless, as shown in figure 4(b), deviations from the predicted inertial-range scaling persist for scales much larger than $\lambda _1$.

We interpret the observed behaviour as an extended contribution from small scales to the integral relation (3.3). For $r>\lambda _1$, the first term always gives a contribution ${\simeq }Q(0)\,\lambda _1^3/(3\,r) = C(0)\,\lambda _1/r$, obtained by evaluating the integral over the interval $0< r^\prime <\lambda _1$. Furthermore, separations $r^\prime$ in the inertial range contribute to both integrals a term $\propto \varepsilon ^{4/3} r^{-2/3} (r/L)^{\zeta _4-4/3}$, with a universal constant determined by the fourth-order structure function, independent of Reynolds number. Balancing these two terms, we find that the first contribution is dominant as long as $C(0)\,\lambda _1/\eta \gg (\varepsilon ^{3/2}/\nu ^{1/2})\,(r/\eta )^{\zeta _4-1} (\eta /L)^{\zeta _4-4/3}$, which is satisfied for $r \ll \lambda _2 =\eta \, [(L/\eta )^{4/3-\zeta _4}\,C(0)/(\varepsilon ^{3/2}/\nu ^{1/2})\,\lambda _1/\eta ]^{1/(\zeta _4-1)} \sim \eta \,\textit {Re}^{\alpha }$ where $\alpha = [1-3\,\zeta _4/4+3\,\gamma /2-\chi _4/2]/(\zeta _4-1)$. The log-normal approximation gives $\alpha \approx 0.189$, which is smaller than $3/4$, consistently ensuring that $\lambda _2\ll L$. This second crossover scale is much larger than $\lambda _1$, hence ensuring the existence of a range of separations $\lambda _1\ll r \ll \lambda _2$ over which the correlations of acceleration behave as $C(r) \simeq C(0)\,\lambda _1/r$ and this range increases with $\textit {Re}$. The numerical data of figure 4(b) confirm this picture. The scaling observed at $r\gtrsim 10\,\eta$ extends further in the inertial range as $\textit {Re}$ increases, and corresponds to $C(r) \propto 1/r$ with a constant that depends weakly on $\textit {Re}$. Both this regime and the small-scale parabolic approximation of the correlation can be matched by the following ad hoc formula

(3.4)\begin{equation} C(r) \approx \frac{A_2(\textit{Re})\,\varepsilon^{3/2}}{\nu^{1/2}[1+(r/\lambda_1)^2]^{1/2}}. \end{equation}

This approximation, shown as a solid line in figure 4(b), is in good agreement with numerical data. In the following, we will use this formula to coarse-grain the particle dynamics.

3.2. Particle accelerations

We focus here on the statistical properties of the acceleration $\boldsymbol {a}_{p} = \mathrm {d}\boldsymbol {v}_{p}/\mathrm {d}t$ of inertial particles. Figure 5(a) shows its variance as a function of the Stokes number. Our measurements agree with those of Bec et al. (2006) and, as they span larger values of the Reynolds numbers, they allow us to substantiate and extend several observations made in that work.

Figure 5. (a) Variance of particle acceleration for the two Reynolds numbers as a function of the Stokes number. The two solid curves are corresponding fits of the form (3.5) with $b=0.42$ and $c = 0.17$. (b) Relative discrepancy in acceleration variance in log–log coordinates. The solid curve is $\varDelta _a = \exp (-b/\textit {St}^{1/2})$.(c) Time autocorrelations of particle acceleration for $R_\lambda = 290$ and various $\textit {St}$. (d) Integral correlation time of acceleration. The solid line corresponds to (3.6) with $\tau _{I}(0) = 2.15\tau _\eta$, $b=0.42$ and $d = 0.2$.

First, we observe that our data, corresponding to two different Reynolds numbers, collapse reasonably well on the top of each other when plotted as a function of $\textit {St}$ and rescaled by the acceleration variance of tracers. This can be seen in figure 5(b), which shows the relative discrepancy in acceleration variance $\varDelta _a \equiv [\langle |\boldsymbol {a}|^2\rangle - \langle |\boldsymbol {a}_{p}|^2\rangle ] / \langle |\boldsymbol {a}|^2\rangle$. Although a weak Reynolds-number dependence is noticeable at very small Stokes numbers, one difficulty distinguishes deviations from possible statistical or numerical errors. Therefore, most effects of intermittency are accounted for by the factor $A_2(\textit {Re})$ introduced in § 3.1. This suggests that acceleration variance can be approximated as $\langle |\boldsymbol {a}_{p}|^2\rangle \approx \langle |\boldsymbol {a}|^2\rangle [1-\varDelta _a(\textit {St})]$, where $\varDelta _a$ is a non-dimensional function of the Stokes number with no significant dependence on $\textit {Re}$.

The second observation is an abrupt reduction in the acceleration variance at small but finite values of $\textit {St}$. There is a drop of over 25 % from $St=0$ to $St=0.1$, which we interpret as a consequence of preferential sampling, specifically of particle ejection from violent small-scale vortical structures. Our data suggest that the relative discrepancy $\varDelta _a$ increases faster than any a power law of $\textit {St}$. This is evidenced by its convexity when plotted in log–log coordinates in figure 5(b), indicating that the acceleration variance may have an essential singularity at $\textit {St}=0$. Such a dependence on Stokes number has been observed previously for the rate at which fold caustics occur (Wilkinson, Mehlig & Bezuglyy 2006), a phenomenon also coined the sling effect (Falkovich et al. 2002). These same events drive the abrupt depletion observed for energy dissipation in § 2.2 and here for acceleration variance. Figures 3(b) and 5(b) show that both discrepancies are well-fitted by a curve $\propto \exp (-c/\textit {St}^{1/2})$, where $c$ depends weakly on $\textit {Re}$. This can be interpreted as a contribution from the probability that the local Stokes number $\tau _{p}|\boldsymbol {\nabla }\boldsymbol {u}|$ is sufficiently large for the particle to detach from the flow, and thus that $\tau _\eta |\boldsymbol {\nabla }\boldsymbol {u}|\gtrsim \textit {St}^{-1}$. At high $\textit {Re}$, the distribution of turbulent velocity gradients is known to display stretched-exponential tails with an exponent ${\approx }1/2$ (see Yeung, Sreenivasan & Pope 2018), consistent with the behaviour of $\varDelta _a$. However, Buaria et al. (2019) found that the constant in the exponential has a significant dependence on the Reynolds number. Therefore, to further refine our discussion, it will be necessary to better understand this dependence in future studies.

Deviation to this singular behaviour occurs at $\textit {St}\gtrsim 1$. Preferential sampling becomes less important, and acceleration statistics are dominated by the particle delay on the flow: their velocity is given by low-pass filtering the fluid velocity over time scales smaller than $\tau _{p}$ (see Bec et al. 2006). Gorokhovski & Zamansky (2018) used such considerations to estimate $\langle |{\boldsymbol {u}}-\boldsymbol {v}_{p}|^2\rangle \simeq \langle |{\boldsymbol {u}}(t)-{\boldsymbol {u}}(t-\tau _{p})|^2\rangle \propto \varepsilon \tau _{p}$, where the last relation assumes $\tau _{p}\gg \tau _\eta$ and uses the inertial-range scaling of the second-order Lagrangian structure function. Consequently, the variance of acceleration approaches a power-law $\langle |\boldsymbol {a}_{p}|^2\rangle \propto St^{-1}$ when $\textit {St}\gg 1$.

The two asymptotics $\textit {St}\ll 1$ and $\textit {St}\gg 1$ can be matched by the ad hoc formula

(3.5)\begin{equation} \langle |\boldsymbol{a}_{p}|^2 \rangle \approx \frac{A_2(\textit{Re})\,\varepsilon^{3/2}}{\nu^{1/2}} \,\frac{1-\exp({-}b/\textit{St}^{1/2})}{(1+c\,\textit{St}^2)^{1/4}}, \end{equation}

which, as seen from figure 5(a), gives a fairly good approximation of particle acceleration variance up to $\textit {St}\approx 7$. Note that this fitting formula differs from other proposals, such as the one suggested by Gorokhovski & Zamansky (2018), which aimed to also capture response times larger than the large-eddy turnover time $\tau _{L} = u_{rms}^2/\varepsilon$. As the response times of our particles lie below $\tau _L$ (i.e. are such that $\textit {St}\ll \textit {Re}^{1/2}$), we use hereafter (3.5).

We now turn to two-time statistics of the particle acceleration, focusing on the autocorrelation $\mathcal {A}(t) = \langle \boldsymbol {a}_{p}(t)\boldsymbol {\cdot }\boldsymbol {a}_{p}(0) \rangle /\langle |\boldsymbol {a}_{p}|^2 \rangle$. The results are shown in figure 5(c). Figure 5(d) shows the integral time $\tau _{I} = \int |\mathcal {A}(t)|\,\mathrm {d}t$ as a function of $\textit {St}$. At small Stokes numbers, it approaches the value for tracers, $\tau _{I}(0) \approx 2.15\tau _\eta$. Deviations occur due again to preferential sampling. Ejection from small-scale vortical structures leads to particles concentrating in regions where the local dissipative time scale is larger than its average. Dimensionally, we expect $\tau _\eta ^{\!@{part}}/\tau _\eta \simeq [\langle |\boldsymbol {a}|^2\rangle / \langle |\boldsymbol {a}_{p}|^2\rangle ]^{1/3}$ and, assuming that $\tau _{I} \approx \tau _\eta ^{\!@{part}}$, we get $\tau _{I}(\textit {St}) \simeq \tau _{I}(0) [ 1-\exp (-b/\textit {St}^{1/2})]^{-1/3}$ for $\textit {St}\ll 1$. On the other hand, for large Stokes numbers, the particle response time effectively filters out all the flow time scales below it, resulting in $\tau _{I}(\textit {St})\propto \tau _{p}$. These two regimes can be matched using the fitting formula

(3.6)\begin{equation} \tau_{I}(\textit{St}) = \frac{\tau_{I}(0) [ 1+d\,\textit{St}]^{5/6}}{[ 1-\exp({-}b/\textit{St}^{1/2})]^{1/3}}, \end{equation}

where $\tau _{I}(0) = 2.15\tau _\eta$ is the value measured from tracers, $b=0.42$ is obtained from the acceleration variance, and $d = 0.2$ provides a good agreement with the data of figure 5(d).

To complete this survey, we examine the spatially averaged particle acceleration $\langle \boldsymbol {a}_{p} \rangle _\ell ({\boldsymbol {x}},t)$. This quantity is defined as the instantaneous average acceleration of all particles that, at time $t$, are located within a ball $\mathcal {B}_{\ell }$ of diameter $\ell$ centred at position ${\boldsymbol {x}}$. We are particularly interested in the statistical properties of $\langle \boldsymbol {a}_{p} \rangle _\ell$ when conditioned on the local value of the spatially averaged dissipation rate $\varepsilon _\ell$, which is calculated in the same region $\mathcal {B}_{\ell }$ at the same time. Our goal is to account for the intermittency and variability of turbulence using Kolmogorov's refined similarity hypothesis. This hypothesis suggests that the statistical properties of turbulent quantities at a scale $\ell$ should be expressed in terms of the local dissipation $\varepsilon _\ell$, rather than its global average $\varepsilon$. This approach allows us to relate local fluctuations in small-scale quantities, such as acceleration, to the inertial-range fluctuations of the dissipation field. According to dimensional analysis, the conditional statistics of the coarse-grained particle acceleration $\langle \boldsymbol {a}_{p} \rangle _\ell$, once normalised by $(\varepsilon _\ell ^3/\nu )^{1/4}$, should depend only on two parameters: the local Stokes number $\textit {St}_{\ell } = \tau _{p}/(\nu /\varepsilon _\ell )^{1/2}$, which is obtained by non-dimensionalising $\tau _{p}$ with the local Kolmogorov time $(\nu /\varepsilon _\ell )^{1/2}$; and the local Reynolds number $\textit {Re}_{\ell } = u_\ell \ell /\nu = \varepsilon _\ell ^{1/3}\ell ^{4/3}/\nu$, which characterises the instantaneous turbulence intensity in the region of size $\ell$.

We can express the conditional mean-squared coarse-grained acceleration of particles as

(3.7)\begin{equation} \left\langle |\langle \boldsymbol{a}_{p} \rangle_\ell|^2 \,\middle | \,\varepsilon_\ell \right\rangle \propto \frac{1}{\ell^3} \int_{0}^{\ell/2} \left\langle \boldsymbol{a}_{p}({\boldsymbol{r}},t)\boldsymbol{\cdot}\boldsymbol{a}_{p}(0,t) \, \middle | \, \varepsilon_\ell \right\rangle r^2\,\mathrm{d}r. \end{equation}

Based on our earlier analysis of the spatial correlations of the fluid acceleration and the approximation (3.4), we can assume that for $\textit {Re}_{\ell }\gg 1, \langle \boldsymbol {a}_{p}({\boldsymbol {r}},t)\boldsymbol {\cdot }\boldsymbol {a}_{p}(0,t) \,|\, \varepsilon _\ell \rangle \approx \langle |\boldsymbol {a}_{p}|^2 \,|\, \varepsilon _\ell \rangle / [1+(r/\lambda _1(\varepsilon _\ell ))^2]^{1/2}$, where $\lambda _1(\varepsilon _\ell ) \propto \eta (\varepsilon _\ell )\,\textit {Re}_{\ell }^{(\gamma -\chi _4)/2}$ is the cut-off scale of acceleration spatial correlations associated with the local conditioning dissipation $\varepsilon _\ell$. The conditional acceleration variance $\langle |\boldsymbol {a}_{p}|^2 \,| \, \varepsilon _\ell \rangle$ is obtained from (3.5) by replacing $\varepsilon$, $\textit {Re}$ and $\textit {St}$ with their local values $\varepsilon _\ell$, $\textit {Re}_{\ell }$ and $\textit {St}_{\ell }$ given by the instantaneous spatially averaged dissipation. Using $\lambda _1(\varepsilon _\ell ) \propto \ell \,\textit {Re}_{\ell }^{-\beta }$ with $\beta = 3/4 + (\chi _4-\gamma )/2 \approx 0.819$, we obtain

(3.8)\begin{equation} \left\langle |\langle \boldsymbol{a}_{p} \rangle_\ell|^2 \,\middle | \,\varepsilon_\ell \right\rangle \approx \frac{e}{\textit{Re}_{\ell}^{\beta}} \frac{A_2(\textit{Re}_{\ell})\,\varepsilon_\ell^{3/2}}{\nu^{1/2}}\,\frac{1-\exp\left({-}b/\textit{St}_{\ell}^{1/2}\right)}{\left(1+c\,\textit{St}_{\ell}^2\right)^{1/4}}, \end{equation}

with $e>0$. Local fluctuations about this average are described by the coarse-grained variance of acceleration, which we can obtain by replacing the dissipation rate with the conditioning value in (3.5). We get

(3.9)\begin{equation} \left\langle \langle |\boldsymbol{a}_{p}|^2 \rangle_\ell-|\langle \boldsymbol{a}_{p} \rangle_\ell|^2 \,\middle | \,\varepsilon_\ell \right\rangle\approx \frac{A_2(\textit{Re}_{\ell})\,\varepsilon_\ell^{3/2}}{\nu^{1/2}}\,\frac{1-\exp\left({-}b/\textit{St}_{\ell}^{1/2}\right)}{\left(1+c\,\textit{St}_{\ell}^2\right)^{1/4}}\left(1-\frac{e}{\textit{Re}_{\ell}^\beta}\right) . \end{equation}

These approximations are valid under the conditions $Re_\ell \gg 1$ and $\ell \ll \lambda _2$. The assumption $Re_\ell \gg 1$ ensures that turbulent scaling and refined similarity hypotheses apply, while $\ell \ll \lambda _2$ is required to match the observed scaling $C(r)\propto r^{-1}$ of acceleration spatial correlations.

Figure 6 shows scatter plots of the conditional mean-squared coarse-grained acceleration and the coarse-grained variance of acceleration obtained from numerical simulations. The solid lines in the figure represent the predictions (3.8) and (3.9), which are based on the approximations made in the preceding text and fitted parameters. The close agreement between the numerical data and the predictions supports the validity of our approximations.

Figure 6. Coarse-grained statistics of particle acceleration conditioned on the local dissipation $\varepsilon _\ell$ for $R_\lambda = 460$. (a) Mean-square coarse-grained acceleration as a function of the local Stokes number $\textit {St}_{\ell } = \tau _{p}/(\nu /\varepsilon _\ell )^{1/2}$ for various $\textit {Re}_{\ell } = \varepsilon _\ell ^{1/3}\ell ^{4/3}/\nu$ labelled by different colours. Solid lines show prediction (3.8), while symbols correspond to numerical measurements for different $\textit {St}$. (b) Same representation for the coarse-grained variance of particle acceleration. Solid lines are to the fit (3.9) with $e=2$.

3.3. An effective diffusion process

To derive effective equations for the particle coarse-grained dynamics, we combine all ingredients from previous analyses. Using (2.3a,b), we can write the particle velocity as $\boldsymbol {v}_{p}(t) = {\boldsymbol {u}}(\boldsymbol {x}_{p}(t),t) -\tau _{p}\,\boldsymbol {a}_{p}(t)$, allowing us to express its displacement over a time $\delta t$ as

(3.10)\begin{equation} \delta\boldsymbol{x}_{p}(t) \equiv \boldsymbol{x}_{p} (t+\delta t)-\boldsymbol{x}_{p}(t) = \int_t^{t+\delta t}{\boldsymbol{u}}(\boldsymbol{x}_{p}(s),s)\,\mathrm{d}s -\tau_{p}\, \int_t^{t+\delta t}\boldsymbol{a}_{p}(s)\,\mathrm{d}s. \end{equation}

We choose $\delta t$ to be much smaller than the Lagrangian correlation time $\tau _{Lag}$ of ${\boldsymbol {u}}$ to ensure that the fluid velocity along particle path does not vary significantly in $[t,t+\delta t]$. Thus, the first integral on the right-hand side of (3.10) can be approximated as ${\boldsymbol {u}}(\boldsymbol {x}_{p}(t),t)\, \delta t + O(\delta t/\tau _{Lag})^2$. All fluctuations and dependences on particle inertia are entailed in the second integral. Additionally, if we assume that $\delta t$ is much longer than the correlation time $\tau _{I}(\textit {St})$ of the particle acceleration, we can apply the central-limit theorem and write

(3.11)\begin{equation} \int_t^{t+\delta t}\boldsymbol{a}_{p}(s)\,\mathrm{d}s \stackrel{\mathrm{law}}{\sim} \mathcal{N}\!\left(\overline{\boldsymbol{a}_{p}}\, \delta t, [\overline{\boldsymbol{a}_{p}\otimes\boldsymbol{a}_{p}}-\overline{\boldsymbol{a}_{p}}\otimes \overline{\boldsymbol{a}_{p}}]\, \tau_{I}\,\delta t\right) + O(\delta t/\tau_{I})^{3/2}, \end{equation}

where $\otimes$ is the outer product and $\mathcal {N}\!(\boldsymbol {m}, \boldsymbol{\mathsf{C}})$ denotes a multivariate normal random variable with mean $\boldsymbol {m}$ and covariance matrix $\boldsymbol{\mathsf{C}}$. The Lagrangian time average $\overline {({\cdot })}$ introduced here is obtained by time integration along particle paths over the interval $[t,t+T]$, assuming the limit $T/\tau _{I}\to \infty$. It contains information about the turbulent state in which the particle is at the initial time $t$ and is crucial to account for inertial-range fluctuations. To estimate this time average, we use an Eulerian spatial average over a coarse-graining scale $\ell$, so that $\overline {({\cdot })} \simeq \langle {\cdot } \rangle _\ell (\boldsymbol {x}_{p}(t),t)$. This estimate assumes that $T$ is chosen of the order of the turnover time $\tau _\ell = \varepsilon ^{-1/3}\ell ^{2/3}$ associated with $\ell$, and hence that $\tau _\ell \gg \tau _{I}(\textit {St})$. Preferential sampling by particles, which naturally arises from the Lagrangian average, is now accounted for by evaluating the Eulerian average at the current particle position $\boldsymbol {x}_{p}$.

Under these assumptions, we can now express the particle displacement as

(3.12)\begin{equation} \delta\boldsymbol{x}_{p}(t) \approx \left[{\boldsymbol{u}}(\boldsymbol{x}_{p}(t),t)-\tau_{p}\,\langle \boldsymbol{a}_{p}(t) \rangle_\ell\right] \delta t + \boldsymbol{\sigma}_\ell(\boldsymbol{x}_{p}(t), t)\,\delta \boldsymbol{W}(t). \end{equation}

Here, $\delta \boldsymbol {W}$ denotes the increment of the three-dimensional Wiener process, and $\boldsymbol {\sigma }_\ell$ is a tensorial diffusion coefficient that satisfies

(3.13)\begin{equation} \tfrac{1}{2}\boldsymbol{\sigma}_\ell\,\boldsymbol{\sigma}_\ell^{\mathsf{T}} = \boldsymbol{\mathsf{D}}_\ell \quad\mbox{with}\ \boldsymbol{\mathsf{D}}_\ell= \tfrac{1}{2}\tau_{p}^2 \tau_{I}\, [\langle \boldsymbol{a}_{p}\otimes\boldsymbol{a}_{p} \rangle_\ell-\langle \boldsymbol{a}_{p} \rangle_\ell\otimes \langle \boldsymbol{a}_{p} \rangle_\ell]. \end{equation}

This diffusion coefficient not only depends on the particle response time and the coarse-graining scale, but also fluctuates in space and time. Taking the limit $\delta t \to 0$ while keeping $\tau _{I}\ll \delta t\ll \tau _\ell = \varepsilon _\ell ^{1/2}\ell ^{2/3}$, we can write the effective displacement (3.12) as the stochastic differential equation

(3.14)\begin{equation} \mathrm{d}\kern0.7pt\boldsymbol{x}_{p}(t) \approx \left[{\boldsymbol{u}}(\boldsymbol{x}_{p}(t),t)-\tau_{p}\,\langle \boldsymbol{a}_{p}(t) \rangle_\ell\right] \mathrm{d} t +\boldsymbol{\sigma}_\ell(\boldsymbol{x}_{p},t)\,\mathrm{d} \boldsymbol{W}(t), \end{equation}

where $\boldsymbol {\sigma }_\ell$ is given by (3.13). The diffusion appears here as a multiplicative noise, which we define using the Itô convention. This is imposed by the requirement that in the statistical steady state, the average particle velocity should vanish, i.e. $\langle \mathrm {d}{\kern0.9pt}\boldsymbol {x}_{p}(t)/\mathrm {d}t \rangle = 0$. Since $\langle {\boldsymbol {u}}(\boldsymbol {x}_{p}(t),t) \rangle = 0$ and $\langle \langle \boldsymbol {a}_{p} \rangle _\ell \rangle = \langle \langle \boldsymbol {a}_{p} \rangle \rangle _\ell = 0$, the contribution of noise should vanish as well.

The proposed model (3.14) for particle dynamics share some similarities with the model introduced by Fevrier et al. (2005). In both cases, the drift term, the ‘mesoscopic Eulerian particle velocity’ in their work, is the sum of the fluid velocity and a residual one. In our model, this residual velocity is proportional to the filtered particle acceleration. Both models also include a noise term. However, while the ‘quasi-Brownian velocity’ of Fevrier et al. (2005) satisfies a molecular chaos assumption and is uncorrelated in space, we identify it in our model as a diffusion with a space–time dependent coefficient that fluctuates due to turbulent agitation. As a result, this contribution is correlated over inertial-range separations.

Particle inertia affect both drift and diffusion in the stochastic equation (3.14). These two contributions have different weights at different scales. They balance each other at a scale $\ell _{diff}$, estimated as $\ell _{diff} = \langle \boldsymbol{\mathsf{D}}_\ell ^{ii} \rangle /[ \tau _{p} \langle |\langle \boldsymbol {a}_{p} \rangle _\ell |^2 \rangle ^{1/2} ]$. Diffusion dominates at scales smaller than $\ell _{diff}$ and is negligible at larger scales. Thus, the diffusion term is relevant only when $\ell _{diff}$ is larger than the coarse-graining scale $\ell$. Using considerations from the previous subsection, and in particular (3.8) and (3.9) for the statistics of the coarse-grained acceleration, we can approximate for $\ell \gg \eta$:

(3.15)\begin{equation} \frac{\ell_{diff}}{\ell} = \frac{\tau_{p}}{2\,\ell} \frac{\left\langle\tau_{I} \left[ \langle |\boldsymbol{a}_{p}|^2 \rangle_\ell-|\langle \boldsymbol{a}_{p} \rangle_\ell|^2\right]\right\rangle}{\langle |\langle \boldsymbol{a}_{p} \rangle_\ell|^2 \rangle^{1/2}} \simeq \varPsi(\textit{St})\left( \frac{\ell}{\eta}\right)^{{-}1+(2/3)(\gamma+\beta)}. \end{equation}

The exponent is negative, indicating that the diffusive scale becomes very small when the coarse-graining scale $\ell$ is far inside the inertial range. Numerical measurements of $\ell _{diff}$, reported in figure 7(a), obtained from the coarse-grained statistics of particle acceleration, confirm the power-law behaviour (3.15) at $\ell \gg \eta$. We also observe that, for the moderate values of the Stokes number considered, $\ell _{diff}$ is always smaller than $\ell$. Based on previous acceleration correlation measurements, we expect that the constant $\varPsi$ behaves as

(3.16)\begin{equation} \varPsi(\textit{St}) \propto \frac{\textit{St}(1+d\,\textit{St})^{5/6}[1-\exp({-}b/\textit{St}^{1/2})]^{1/6}}{(1+c\,\textit{St}^2)^{1/8}}. \end{equation}

This prediction, shown as a solid curve in figure 7(b), compares well with the numerical measurements shown as circles. Extrapolating this behaviour to higher Stokes numbers, we obtain that $\varPsi \sim \textit {St}^{3/2}$, implying that neglecting diffusion requires choosing a coarse-graining scale such that $\ell /\eta \gg \textit {St}^{(3/2)/[1-(2/3)(\gamma +\beta )]}\approx \textit {St}^{3.73}$. Note that this condition applies to particle response times in the inertial range but still smaller than the fluid velocity Lagrangian correlation time $\tau _{Lag}$. This condition is more restrictive than the classical idea that the particle response time should be smaller than the eddy turnover time $\tau _\ell = \varepsilon ^{-1/3}\ell ^{2/3}$ associated with the coarse-graining scale, which would instead lead to $\ell /\eta \gg \textit {St}^{3/2}$.

Figure 7. (a) Ratio of the diffusive scale $\ell _{diff}$ to the coarse-graining scale $\ell$ for various Stokes numbers at $R_\lambda = 460$. The scale $\ell _{diff}$ is the ‘Batchelor’ scale above which the diffusive term in (3.14) becomes negligible. The dashed line has a slope $-1+(2/3)(\gamma +\beta ) \approx -0.402$, as predicted from the log-normal approximation. (b) Coefficient $\varPsi$ of the power law (3.15). Symbols are numerical measurements and the solid curve is the prediction (3.16).

In this section, we have shown that the coarse-grained dynamics of inertial particles, when averaged over time, can be effectively modelled using a stochastic equation that includes both drift and diffusion terms. The contributions from particle inertia, which distinguish heavy particles from fluid elements, are governed by the coarse-grained particle acceleration. This acceleration, which fluctuates both spatially and temporally, serves as a clear indicator of turbulent activity. Our analysis shows that for sufficiently large temporal coarse-graining scales, or equivalently small Stokes numbers, the diffusive effects become negligible, leading to dynamics that are predominantly governed by drift. In the following section, we will focus on this asymptotic regime and develop a further level of modelling that allows us to derive an effective dynamics for the Eulerian coarse-grained density of particles, shifting the emphasis from temporal to spatiotemporal averaging.

4. Particle transport as an Eulerian ejection process

4.1. Model dynamics for the particle density

In the previous section, we introduced an effective velocity field $\boldsymbol {v}_{p}^{eff} = {\boldsymbol {u}} - \tau _{p},\langle \boldsymbol {a}_{p} \rangle _{\ell }$, which describes the Lagrangian dynamics of particles at large temporal averaging scales (corresponding to the limit of weak inertia) in terms of the coarse-grained particle acceleration. The effective dynamics, characterised by both drift and diffusion terms, simplify as diffusion becomes negligible at sufficiently large scales or small Stokes numbers. While this approach provides insights into how inertia influences particle dynamics over time, to fully capture their spatial distribution in turbulent flows, we shift to an Eulerian framework. This involves reformulating the problem by tracking the evolution of the coarse-grained particle density $\langle \rho _{p} \rangle _{\ell }$ within a volume $\mathcal {B}_{\ell }$ of size $\ell$ – figure 8(a), combining spatial and temporal averages. This transition is not merely a change in perspective but a complementary methodology, enabling us to model the macroscopic behaviour of particles as a continuum, thus better quantifying how they cluster or disperse across different scales of turbulence.

Figure 8. Sketch of the ejection process. (a) A quasi-Lagrangian viewpoint is adopted to emphasise discrepancies due to inertia. Outgoing fluxes correspond to particles that acquired large-enough accelerations within $\mathcal {B}_{\ell }$. (b) The reference cell $(i,j,k)$ ejects particles with a rate $\varPhi _t(i,j,k)$ to all its neighbours and receives their individual contributions.

We adopt a quasi-Lagrangian approach and follow the control volume in its motion with the fluid velocity ${\boldsymbol {u}}$, while considering its exchanges with its Eulerian neighbours. To evaluate the fluxes due to particles’ inertia at the boundary $\partial \mathcal {B}_{\ell }$ of the control volume, we distinguish between outgoing and incoming fluxes. Some particles leave the volume because they have acquired a large-enough acceleration inside $\mathcal {B}_{\ell }$, and the outgoing flux should thus be controlled by the coarse-grained acceleration $\langle \boldsymbol {a}_{p} \rangle _{\ell }$ computed inside the reference volume. This flux can be expressed as

(4.1)\begin{equation} \varPhi_t({\boldsymbol{x}}) \approx \int_{\partial\mathcal{B}_{\ell}({\boldsymbol{x}})} \left\langle\left(-\tau_{p}\,\boldsymbol{a}_{p} \boldsymbol{\cdot} \boldsymbol{n}\right)\theta(-\boldsymbol{a}_{p} \boldsymbol{\cdot} \boldsymbol{n})\, \rho_{p}\right\rangle_\ell \mathrm{d}S \approx 3\ell^2 \tau_{p}\langle \rho_{p} \rangle_{\ell}\,|\langle \boldsymbol{a}_{p} \rangle_{\ell}|. \end{equation}

Here, $\theta$ denotes the Heaviside function, $\boldsymbol {n}$ is the unit vector normal to the surface of $\mathcal {B}_{\ell }$, and the average is taken over accelerations satisfying $\boldsymbol {a}_{p}\boldsymbol {\cdot } \boldsymbol {n} < 0$ to account only for outgoing particles. Assuming isotropic distribution of the outgoing flux, this signed average can be approximated by $(1/2) |\langle \boldsymbol {a}_{p} \rangle _{\ell }|$. The control volume $\mathcal {B}_{\ell }$ is chosen as a cube with edge length $\ell$, and the spatial domain is tiled by such cubes, see figure 8(b). The time evolution of the mass $\ell ^3\langle \rho _{p} \rangle _{\ell }$ of particles contained in the cell $(i,j,k)$ is then given by

(4.2)\begin{align} \frac{\mathrm{D}}{\mathrm{D}t} \left[\ell^3 \langle \rho_{p} \rangle_{\ell}\right] &={-}\varPhi_t(i,j,k) +\frac{1}{6}\left[\varPhi_t(i-1,j,k) +\varPhi_t(i+1,j,k)+ \varPhi_t(i,j-1,k)\right. \nonumber\\ &\left.\quad +\varPhi_t(i,j+1,k) +\varPhi_t(i,j,k-1) +\varPhi_t(i,j,k+1)\right]. \end{align}

Here $\mathrm {D}/\mathrm {D}t = \partial _t + {\boldsymbol {u}}\boldsymbol {\cdot }\boldsymbol {\nabla }$ denotes the material derivative along the trajectories of fluid elements. Mass is lost from the outgoing flux $\varPhi _t(i,j,k)$ in the reference cell and gained from the outgoing flux coming from its six neighbours on the cubic tiling. The right-hand side of (4.2) corresponds to the discrete Laplacian of the outgoing flux $\varPhi _t$. By considering the mass evolution on scales much larger than the coarse-graining scale $\ell$, we can write a continuous limit which reads

(4.3)\begin{equation} \partial_t \langle \rho_{p} \rangle_{\ell} + {\boldsymbol{u}} ({\boldsymbol{x}},t) \boldsymbol{\cdot}\boldsymbol{\nabla}\langle \rho_{p} \rangle_{\ell} \approx \nabla^2 \left[ \kappa_{\ell}({\boldsymbol{x}},t)\,\langle \rho_{p} \rangle_{\ell} \right ],\quad\mbox{with}\ \kappa_\ell = \tau_{p}\ell\,|\langle \boldsymbol{a}_{p} \rangle_{\ell}|/2. \end{equation}

The position- and time-dependent coarse-grained diffusion coefficient, $\kappa _{\ell }$, appears inside the Laplacian as expected for an ejection process. This model for particle transport provides a quantitative extension of the phenomenological ideas proposed in Bec & Chétrite (2007). As we will discuss later, the underlying ejection process gives rise to specific features in the probability distribution of the spatially averaged density, $\langle \rho _{p} \rangle _{\ell }$.

The diffusive term in (4.3) can be expressed as the divergence of the flux vector $\boldsymbol {\varphi }_t = -\kappa _{\ell }\,\boldsymbol {\nabla } \langle \rho _{p} \rangle _{\ell } -\langle \rho _{p} \rangle _{\ell }\,\boldsymbol {\nabla } \kappa _{\ell }$, which consists of two distinct contributions. The first corresponds to osmotic forces, resulting in classical Fickian diffusion that enhances mixing alongside fluid advection. The second arises from turbophoretic forces due to convection by the velocity $\boldsymbol {\nabla } \kappa _{\ell }$, which drive particles from regions with high $\kappa _{\ell }$, characterised by strong particle accelerations and high turbulent activity, to regions with low $\kappa _{\ell }$. The turbophoretic contribution is responsible for the preferential sampling of particles in the inertial range, as qualitatively discussed in § 2.2. To determine whether turbophoretic forces are strong enough to induce significant concentration fluctuations and inertial-range voids, we need to compare the magnitudes of the terms in $\boldsymbol {\varphi }_t$. In particular, turbophoretic forces dominate when $|\boldsymbol {\nabla }\kappa _\ell |/\kappa _{\ell }>|\boldsymbol {\nabla }\langle \rho _{p} \rangle _{\ell }| /\langle \rho _{p} \rangle _{\ell }$, which implies that the scale of variation of the diffusion coefficient (and thus of the particle acceleration) should be smaller than that of density.

The balance between fluid flow convection and turbophoretic diffusion can be characterised at a given coarse-graining scale $\ell$ by a dimensionless Péclet number that we define as

(4.4)\begin{equation} \textit{Pr}_{\ell} = \frac{\ell\delta u(\ell)}{\langle \kappa_\ell^2 \rangle^{1/2}} = \frac{2\sqrt{3}[S_2^{{\parallel}}(\ell)]^{1/2}}{\tau_{p}\left\langle|\langle \boldsymbol{a}_{p} \rangle_{\ell}|^2\right\rangle^{1/2}}. \end{equation}

Here, $\delta u(\ell )$ is the typical fluid velocity fluctuation at scale $\ell$, which is estimated by the square root of the second-order longitudinal structure function $S_2^{\parallel }$. In the inertial-range, $S_2^{\parallel } \sim \ell ^{\zeta _2}$ with $\zeta _2\approx 0.696$. Equation (3.8) shows that for the spatially averaged acceleration, $\langle |\langle \boldsymbol {a}_{p} \rangle _{\ell }|^2\rangle \propto A_2(\textit {Re}_{\ell })\,\textit {Re}_{\ell }^{-\beta }\sim \ell ^{(4/3)(\gamma -\beta )}$, which leads to the scaling behaviour $\textit {Pr}_{\ell } \sim \ell ^{\delta }$ with $\delta = \zeta _2/2 -(2/3)(\gamma -\beta )\approx 0.84$ for coarse-graining scales $\ell$ in the inertial range. Figure 9(a) confirms this power-law dependence. Regarding the dependence on the Stokes number, we have $\textit {Pr}_{\ell } \sim \textit {St}^{-1}$ when $\textit {St}\ll 1$, as shown in figure 9(b). The numerical data indicate that the Péclet number can reach values larger than 1 for both $\ell \gg \eta$ and $\textit {St}\ll 1$. The scaling laws for small Stokes number and large averaging scale give $\textit {Pr}_{\ell }\sim (\ell /\eta )^\delta /\textit {St}$, which results in a Péclet number much higher than unity when $\ell /\eta \gg \textit {St}^{1/\delta } \sim \textit {St}^{1.19}$. This scaling is distinct from those discussed in the previous section based on Lagrangian considerations.

Figure 9. Scale-dependent Péclet number $\textit {Pr}_{\ell }$ defined from (4.4) shown in (a) as a function of the coarse-graining scale $\ell$ for various Stokes numbers and in (b) as a function of the Stokes number for various coarse-graining scales $\ell$. The dashed line in (a) shows a behaviour $\textit {Pr}_{\ell } \propto (\ell /\eta )^{\delta }$ with $\delta = \zeta _2/2 + (2/3)\,(\gamma -\beta ) \approx 0.84$ that is expected for $\ell \gg \eta$. The dashed line in (b) corresponds to $\textit {Pr}_{\ell }\propto \textit {St}^{-1}$ that prevails at small values of the Stokes number.

4.2. Distribution of the coarse-grained density

Based on our previous arguments, we anticipate that for sufficiently large scales, the Péclet number $\textit {Pr}_{\ell }$ defined in (4.4) captures alone dependences upon both the Stokes number $\textit {St}$ and the coarse-graining scale $\ell$. This asymptotic regime corresponds to the range of parameter values where the approximation (4.3) accurately describes particle dynamics, and we expect that their clustering behaviour will primarily depend on $\textit {Pr}_{\ell }$.

We start with examining the radial distribution function, or pair distribution function $g(\ell )$, which describes the probability of finding two particles at a distance $\ell$, normalised by the probability for a uniform distribution. It can be expressed through the second-order moment of the coarse-grained density $\langle \rho _{p} \rangle _{\ell }$, namely $g(\ell ) = \langle \langle \rho _{p} \rangle _{\ell }^2\rangle / \langle \langle \rho _{p} \rangle _{\ell }\rangle ^2$. For a uniform distribution, we have $\langle \rho _{p} \rangle _{\ell } \equiv \rho _0 = \langle \langle \rho _{p} \rangle _{\ell }\rangle$, so $g(\ell )=1$. Deviations from uniformity as a function of the scale-dependent Péclet number are shown in figure 10(a). Data associated with different values of the Stokes number collapse onto a unique master curve, when the coarse-graining scale $\ell$ is chosen far enough in the inertial range. This curve shows two distinct scaling regimes: one at moderate Péclet numbers and another at large values.

Figure 10. (a) Deviation from uniformity of the radial distribution function, $g(\ell )-1$, as a function of the scale-dependent Péclet number $\textit {Pr}_{\ell }$ for $R_\lambda =460$ and various Stokes numbers. The dashed lines show the expected behaviours at moderate and large values of $\textit {Pr}_{\ell }$. The solid curve is an approximation $\propto \textit {Pr}_{\ell }^{\zeta _2/\delta -2} / [1+\textit {Pr}_{\ell }/P]^{\zeta _2/\delta }$ with $P = 200$. (b) Probability density function (p.d.f.) of the particle coarse-grained density $\langle \rho _{p} \rangle _{\ell }$ for six different values of the Péclet number (different symbols), each corresponding to two Stokes numbers (different colours) and two coarse-graining scales combining to the same $\textit {Pr}_{\ell }$.

For large values of $\textit {Pr}_{\ell }$, we can express the coarse-grained density as $\langle \rho _{p} \rangle _{\ell } = \rho _0 +\delta \rho$ with $\delta \rho \ll \rho _0$. To leading order, the perturbation satisfies $\partial _t \delta \rho + {\boldsymbol {u}}\boldsymbol {{\cdot }\nabla } \delta \rho = \rho _0 \nabla ^2 \kappa _{\ell }$. For statistically stationary deviations to uniformity, we get $\delta u(\ell )\, \delta \rho \sim \rho _0\kappa _{\ell } / \ell$, which implies that $\delta \rho \sim \textit {Pr}_{\ell }^{-1}$ and the variance scales as $g(\ell )-1 \sim \textit {Pr}_{\ell }^{-2}$. At lower Péclet numbers and higher Stokes numbers, when deviations to uniformity are still small, the velocity contribution is dominated by the large-scale advection, so we have $u_{rms}\delta \rho \sim \kappa _\ell /\ell$. This means that deviations from uniformity depend on $\kappa _\ell$, but not on fluid velocity fluctuations at the scale $\ell$. Thus, we have $\delta \rho \sim \kappa _\ell /\ell \sim \delta u(\ell ) /\textit {Pr}_{\ell }$. Using $\textit {Pr}_{\ell } \sim \ell ^{\delta }$, we obtain the second scaling regime $g(\ell )-1 \sim \textit {Pr}_{\ell }^{\zeta _2/\delta -2}$. A solid curve in figure 10(a) shows an ad hoc approximation matching these two asymptotic laws. It provides a reasonable fit to the numerical measurements.

Let us contextualise our results with respect to previous findings on how particles recover a uniform distribution at large scales. When $\ell$ becomes large, $g(\ell )$ tends to unity and in our approach, we find that $\log g(\ell ) \approx g(\ell )-1 \sim \textit {Pr}_{\ell }^{-2} \sim \tau _{p}^2\ell ^{-1.68}$. Such an algebraic dependence differs from the exponential decay proposed by Reade & Collins (2000) that seems confirmed by the experimental measurements of Petersen, Baker & Coletti (2019) (see also Brandt & Coletti (2022)). The scaling that we observe also significantly deviates from the prediction of Balkovsky et al. (2001) (see also Falkovich, Fouxon & Stepanov (2003)), who proposed that the radial distribution depends primarily on the scale-dependent Stokes number, with $\log g(\ell ) \propto \textit {St}_{\ell }^2 \sim \tau _{p}^2\ell ^{-4/3}$ when $\textit {St}_{\ell }\ll 1$. However, the relevance of $\textit {St}_{\ell }$ has so far been demonstrated only in models assuming the fluid velocity is a white noise process. For instance, for velocities in the Kraichnan ensemble, it has been shown by Bec et al. (2007b) that $\log g(\ell ) \propto (\mathcal {D}_2(\textit {St}_\ell )-3)\log \ell \sim \textit {St}_\ell ^2\,\log \ell$, which is not a pure function of $\textit {St}_{\ell }$. Moreover, the direct numerical simulations of Bec et al. (2007a) at moderate Reynolds numbers suggest that particle distributions primarily depend on a scale-dependent contraction rate $\propto \tau _{p}\ell ^{-5/3}$, without clear evidence of scaling for the radial distribution function in the inertial range. More recent numerics by Bragg et al. (2015) and Ariki et al. (2018) indicate a scaling $\log g(\ell ) \sim \ell ^{-4/3}$, albeit with uncertainties on the Stokes number dependence. Our approach reveals a second power-law regime, persisting up to $\textit {Pr}_{\ell }\approx 100$, where $\log g(\ell ) \sim \textit {Pr}_{\ell }^{\zeta _2/\delta -2}\sim \tau ^{1.17}\ell ^{-1.39}$, potentially masquerading a behaviour $\propto \ell ^{-4/3}$. It would be of interest to reassess previous measurements of the radial distribution function in the light of present findings.

To complement our analysis, we turn to the p.d.f. $p(\langle \rho _{p} \rangle _{\ell })$ of the coarse-grained density. Figure 10(b) displays numerical measurements for six different high values of the Péclet number. Remarkably, data obtained from various combinations of the particle response time $\tau _{p}$ and the coarse-graining scale $\ell$, resulting in the same $\textit {Pr}_{\ell }$, exhibit a reasonable collapse, within the range of statistical errors. This confirms the significance of the scale-dependent Péclet number in characterising density fluctuations. The observed probability distributions manifest distinctive features. Both tails, associated with small and large values of $\langle \rho _{p} \rangle _{\ell }$, are broader than those expected for a Poisson distribution corresponding to a uniform particle density. These deviations can be explained by the ejection process framework developed in Bec & Chétrite (2007). Specifically, we find that large densities occur more frequently than the quasi-Gaussian tail of the Poisson distribution. The p.d.f.s exhibit a subexponential behaviour $p(\langle \rho _{p} \rangle _{\ell })\propto \exp (-C\langle \rho _{p} \rangle _{\ell }\log \langle \rho _{p} \rangle _{\ell })$, which is clearly captured by our data, as evident in figure 11(a).

Figure 11. (a) The p.d.f. of the coarse-grained density $p(\langle \rho _{p} \rangle _{\ell })$ shown for $\ell =64\eta$ represented as a function of $\langle \rho _{p} \rangle _{\ell }\log \langle \rho _{p} \rangle _{\ell }$ to evidence the subexponential tail at large masses. (b) Cumulative distribution function (CDF) of the coarse-grained density for different Stokes numbers, $\ell =64\eta$ and $R_\lambda =460$ (same labels as panel a), showing a power-law behaviour $\propto \langle \rho _{p} \rangle _{\ell }^{\xi }$ at small values. Inset (c) Exponent $\xi$ as a function of the scale-dependent Péclet number, measured for various $\textit {St}$ and coarse-graining scales $\ell$.

Regarding the left-hand tail, quasiempty regions occur also more frequently than in a simple Poisson process. Density distributions follow there a power-law $p(\langle \rho _{p} \rangle _{\ell })\propto \langle \rho _{p} \rangle _{\ell }^{\xi -1}$, as depicted in figure 11(b). This behaviour is again a characteristic feature of ejection processes. To provide a heuristic explanation, we consider the approximation (4.3) of the dynamics, where the evolution of the coarse-grained density along a particle trajectory is given by

(4.5)\begin{equation} \tilde{\rho}_{p} \equiv \langle \rho_{p} \rangle_{\ell}(\boldsymbol{x}_{p}(t),t) \approx \rho_0 \exp \int_{-\infty}^t \nabla^2\kappa_\ell(\boldsymbol{x}_{p}(s),s)\,\mathrm{d}s. \end{equation}

We can thus write the cumulative distribution of the Lagrangian density $\tilde {\rho }_{p}$ for $\tilde {\rho }_{p}\ll \rho _0$ as

(4.6)\begin{align} \mathrm{Pr}(\tilde{\rho}_{p}<\rho) &\approx \mathrm{Pr}\left(\int_{-\infty}^t \nabla^2\kappa_\ell(\boldsymbol{x}_{p}(s),s)\,\mathrm{d}s<\log\frac{\rho}{\rho_0}\right) \nonumber\\ &\approx \max_N \left[\mathrm{Pr}\left(\tau_{\ell} \nabla^2\kappa_\ell < \frac{1}{N} \log\frac{\rho}{\rho_0}\right)\right]^N. \end{align}

Here, we have decomposed the Lagrangian integral of the turbophoretic term into a sum of $N$ equally distributed independent random variables $\tau _{\ell } \nabla ^2\kappa _\ell$, where $\tau _{\ell }$ represents the correlation time of the ejection rate along particle paths. The asymptotic $\rho \ll \rho _0$ behaviour is then obtained by optimising $N$, which represents the number of times mass must be ejected to create a void. If this number is of the order of unity, the above formula samples the (negative) tail of the distribution of $\nabla ^2\kappa _\ell$. A power-law behaviour arises because it is more favourable to choose a value of $N$ of the order of $|\log (\rho /\rho _0)|$, indicating that empty regions are more likely to results from persistent ejections rather than rare, violent events leading to instantaneous voids. Thus, by writing the optimum as $N = -n\log (\rho /\rho _0)$ with $n=O(1)$, the CDF becomes

(4.7)\begin{equation} \mathrm{Pr}(\tilde{\rho}_{p}<\rho) \propto \left(\rho/\rho_0\right)^{\tilde{\xi}}, \quad\mbox{with } \tilde{\xi} ={-}n\log \mathrm{Pr}\left(\tau_{\ell} \nabla^2\kappa_\ell(\boldsymbol{x}_{p}(t),t) <{-}1/n\right). \end{equation}

When the Péclet number is large enough, advection dominates, resulting in the correlation time $\tau _\ell$ being given by the eddy-turnover time at scale $\tau _{\ell }\simeq \ell /\delta u(\ell )$. Consequently, $\mathrm {Pr}(\tau _{\ell } \nabla ^2\kappa _\ell <-1/n) \simeq \mathrm {Pr}(|\ell \delta u(\ell )/\kappa _\ell | < n) \propto \textit {Pr}_{\ell }^{-1}$ for $\textit {Pr}_{\ell }\gg 1$. Furthermore, the exponent $\tilde {\xi }$ is bounded from below by zero in order for the probability distribution of $\tilde {\rho }_{p}$ to be normalisable. Eulerian statistics are then obtained by accounting for the additional factor of $\rho /\rho _0$ involved in the Lagrangian average, because it is itself weighted by the particle density. This finally leads us to write the probability distribution of the Eulerian coarse-grained density as

(4.8)\begin{equation} \mathrm{Pr}(\langle \rho_{p} \rangle_{\ell}<\rho) \propto \left(\rho/\rho_0\right)^{\xi},\quad \mbox{with } \xi = 1+\tilde{\xi} \approx 1+f\max\left[0,\log(\textit{Pr}_{\ell}/{Pe}_\star)\right], \end{equation}

where $f$ is a positive constant. Figure 11(c) displays the measured exponent $\xi$ as a function of the scale-dependent Péclet number. The exponent saturates at $\xi =1$ for $\textit {Pr}_{\ell }<{Pe}_\star \approx 5.5$ and is larger than 1 above that threshold, increasing as $f\log (\textit {Pr}_{\ell }/{Pe}_\star )$ with $f\approx 0.75$ for larger values, confirming the prediction given by (4.8).

4.3. Distribution of voids

We now shift our attention to the large voids that prominently emerge in the spatial distribution of particles. As we observed in § 2.2, the sizes of these empty regions span the entire inertial range, even at moderate Stokes numbers. Our goal here is to investigate to what extent the statistics of these voids can be explained by the effective diffusion (4.3) introduced in § 4.1.

To detect these voids numerically, we rely on the spatially averaged density. They are defined as connected sets of empty cells, identified by a label-propagation algorithm. The volume $\mathcal {V}$ of each void is determined by counting the number of cubes with a volume $\ell ^3$ that it encompasses. While alternative techniques for void detection, such as Delaunay tessellations (see e.g. Gaite 2005), may offer better algorithmic efficiency and the ability to define voids in a parameter-free manner, they yield the same results as presented below. Therefore, we have chosen to continue working with the spatially averaged density, which is central to the model for particle coarse-grained dynamics proposed in § 4.1. Figure 12 displays two-dimensional slices of the three-dimensional distribution of voids for a coarse-graining scale of $\ell =16\eta$ and for three different values of the particle response time. These distributions are shown at the same instant of time and in the same slice as the local kinetic energy dissipation rate in figure 2(a). The comparison of these two figures clearly identifies voids as regions with high turbulent activity. Furthermore, there are evident correlations between the empty regions associated with different Stokes numbers. One such correlation can be observed for the circled greenish structure, where the intensity of voids increases with $\textit {St}$. Conversely, in other cases, exemplified by the lower orangish structure circled with dots, a void that exists at small $\textit {St}$ can be filled by particles with a larger inertia.

Figure 12. Two-dimensional cuts of the instantaneous distribution of voids for $R_\lambda =460$ and (a) $\textit {St}=0.4$, (b) $\textit {St}=1$ and (c) $\textit {St}=2.5$. The time and the position of this slice are the same as in figure 2. Voids are obtained as connected empty cubes of size $\ell = 16\,\eta$.

Figure 13(a) presents the complementary cumulative probability distributions of void volumes obtained for different Stokes numbers and an elementary coarse-graining scale $\ell =16\eta$. These distributions exhibit broad tails at large inertial-range volumes, displaying a distinctive power-law behaviour $\mathrm {Pr}(\mathcal {V}>v) \propto v^{-\zeta }$, where $\zeta \approx 1$ is particularly evident for the highest Stokes number values. These statistics remain robust when using alternative definitions of voids, or when changing either the total number of particles $N_{p}$ or the coarse-graining scale $\ell$. Similar power-law dependencies have been previously observed in the probability distribution of void sizes. In the two-dimensional inverse cascade, Boffetta et al. (2004) found an intermediate regime where the p.d.f. of void areas behaves as $p(a)\propto a^{-1.8}$, independent of the Stokes number, with an exponential cutoff at larger sizes. Goto & Vassilicos (2006) proposed a self-similar distribution of void areas with $p(a)\propto a^{-5/3}$, arising from sweep-stick mechanisms where particles preferentially trace fluid zero-acceleration points. Extending these arguments to three dimensions, Yoshimoto & Goto (2007) predicted a power-law exponent $\zeta = 7/9\approx 0.778$ for the cumulative distribution of void volumes, with reasonable numerical support at moderate values of the Reynolds number. Figure 13(a) showcases this behaviour for comparison. Additional evidence supporting this shallow trend comes from grid-turbulence experiments by Sumbekova et al. (2017) and analyses employing Voronoï diagrams, where they found $p(a)\propto a^{-1.8\pm 0.1}$ for void areas in two-dimensional cross-sections of the three-dimensional particle distribution. Assuming a relationship of the form $p(v) \sim v^{-1/3} p(a)$ with $a\sim v^{2/3}$, these observations suggest $\zeta \approx 0.53\pm 0.07$. However, our data clearly show a steeper slope, even for Stokes numbers exceeding those considered in both Yoshimoto & Goto (2007) and Sumbekova et al. (2017).

Figure 13. (a) Complementary CDF $\mathrm {Pr}(\mathcal {V}>v)$ of the void volumes $\mathcal {V}$ shown for $R_\lambda =460$ and various Stokes numbers, as labelled. The numerical data is represented by symbols, while the solid lines correspond to approximations of the form (4.10) – see text for parameter values. (b) Same data, rescaled to emphasise the log-normal behaviour for $v\gg v_\star$, shown as a solid line, which is independent of the Stokes number. Inset is the measured power-law exponent $\zeta$ as a function of the Stokes number. The solid line corresponds to $\zeta = 1-0.17\log \textit {St}$.

We revisit here void statistics in light of the ejection process that we introduced to model particle dynamics in the inertial range. The probability that the volume of a void exceeds the value $v$ can be estimated as the probability of finding very few particles in a cube of size $\ell \sim v^{1/3}$. This implies that the coarse-grained density is there of the order of or smaller than $v^{-1}$. Thus, we can write $\mathrm {Pr}(\mathcal {V}>v) \sim \mathrm {Pr}(\langle \rho _{p} \rangle _{\ell }\lesssim v^{-1})$. Using the asymptotic behaviour (4.8) for the distribution of $\langle \rho _{p} \rangle _{\ell }$ at small values, we obtain

(4.9) \begin{equation} \mathrm{Pr}(\mathcal{V}>v) \sim (v/\eta^3)^{-\xi({\scriptsize\textit{Pr}_{\ell}})} \sim (v/\eta^3)^{{-}1}\exp({-\log(v/\eta^3)\,f\max[0,\log({\scriptsize\textit{Pr}_{\ell}}/{\scriptsize{\textit Pe}}_\star)]}). \end{equation}

Choosing $\ell$ to be of the order of $v^{1/3}$, we have $\textit {Pr}_{\ell }\sim (\ell /\eta )^\delta /\textit {St}\propto (v/\eta ^3)^{\delta /3}/\textit {St}$, resulting in

(4.10)\begin{equation} \mathrm{Pr}(\mathcal{V}>v) \sim \begin{cases} \textit{St}^{{-}3/\delta}(v/v_\star)^{{-}1} & \mbox{if } v< v_\star, \\ \textit{St}^{{-}3/\delta}(v/v_\star)^{-\zeta(\textit{St})}\exp({-g[\log(v/v_\star)]^2}) & \mbox{if } v\geqslant v_\star. \end{cases} \end{equation}

Here, $g$ is a positive constant, $v_\star \propto \eta ^3\textit {St}^{3/\delta }$ and the exponent $\zeta$ has a logarithmic dependence on the Stokes number of the form $\zeta \approx 1-h\log (\textit {St}/\textit {St}_\star )$, where $h>0$. Figure 13(a) shows such predictions for the distribution of void volumes along with numerical data. Reasonable agreement is obtained by choosing for fitting parameters $v_\star /\eta ^3=4000\,\textit {St}^{3/\delta }$, $g=0.0085$, $h=0.17$ and $\textit {St}_\star =1$. The measurements shown in figure 13(b) corroborate these values. Figure 13(b) represents the rescaled complementary cumulative distribution of void volumes as a function of $v/v\star$. Despite statistical noise, data associated with various Stokes numbers (symbols) seem to collapse for $v>v_\star$ onto the log-normal master curve $\exp [-0.0085[\log (v/v\star )]^2]$ (solid line). The measured exponent $\zeta$ is represented in the inset in figure 13(b). It follows $\zeta \approx 1-0.17\log \textit {St}$ for $\textit {St}\lesssim 1$ and saturates to $\zeta \approx 1$ for larger $\textit {St}$.

It is worth noting that the intermediate power-law behaviour that we observe in the distribution of void sizes can be interpreted in terms of Zipf's law (see e.g. Cristelli, Batty & Pietronero 2012). Samples following this law exhibit coherence and adhere to certain dynamical constraints, which are satisfied when the size dynamics of the objects under consideration can be described as a multiplicative process. In the context of turbophoresis, interpreted as an ejection process, this framework naturally emerges, as the mass of particles ejected from a given cell is proportional to its volume. For such a coherent process, the exponent $\zeta =1$ represents a classical case. It arises when large voids are formed through the merging of smaller, independent voids with uncorrelated histories, as may occur at large Stokes numbers. The growth rate of a large void becomes proportional to the probability of intersecting other empty regions, which, in turn, is proportional to its volume. This process, known as ‘preferential attachment’, leads to an exponent of $\zeta =1$ (see De Marzo et al. 2021).

5. Concluding remarks

In this paper, we have presented convincing evidence that the phenomenon of turbophoresis, previously thought to occur only in turbulent flows containing inhomogeneities, also manifests in statistically homogeneous situations. This effect arises from the instantaneous non-uniformities intrinsic to turbulent flows, spanning the whole inertial range. Our direct numerical simulations clearly illustrate the ejection of inertial particles from highly active regions of the flow, leading to their concentration in calmer regions. Remarkably, this behaviour persists in spatially coarse-grained representations of both the flow and the particles, resulting in strong correlations between the spatially averaged particle concentration and the fluctuations in turbulent kinetic energy dissipation within the inertial range.

The fluctuations in particle acceleration play a crucial role in the turbophoresis process. When particles experience pure Stokes drag, these acceleration fluctuations govern their deviations from fluid motion. Through analytical and phenomenological arguments, as well as a detailed analysis of numerical simulations, we have gained insights into the statistics of particle acceleration. This includes understanding spatial and temporal correlations, as well as the influence of fluid flow intermittency on second-order statistics. Building upon these insights, we have introduced approximations for the inertial-range dynamics of particles in terms of effective diffusion equations with a diffusivity that varies in both space and time. The diffusion coefficient is expressed in terms of the coarse-grained particle acceleration, which, in turn, is determined by local turbulent activity. These approximations hold when spatial averaging scales are sufficiently large or particle inertia is sufficiently small, ensuring that higher-order corrections to this dynamics remain negligible. Our study integrates two complementary perspectives: the Lagrangian view, which involves temporal averaging of particle dynamics; and the Eulerian view, which incorporates both temporal and spatial averaging. The Lagrangian model captures time-dependent effects of inertia, while the Eulerian model provides a macroscopic description of particle distributions as a continuum. Although neither model is fully closed due to the need for detailed statistical input, the preliminary studies that we present offer a foundation for linking particle dynamics with coarse-grained turbulent statistics, particularly under Kolmogorov's (1962) refined similarity hypothesis. An important finding is that inertial-range particle dynamics depend solely on a local Péclet number that quantifies the relative importance of advection by the fluid flow compared with inertia-induced diffusion at a given coarse-graining scale $\ell$. Notably this Péclet number exhibits a non-trivial power-law dependence on the observation scale, $\textit {Pr}_{\ell }\sim (\ell /\eta )^{0.84}/\textit {St}$, where the exponent is prescribed by the intermittent statistics of the fluid velocity and deviates significantly from the value $2/3$ that would be obtained by dimensional analysis.

The diffusive models we have developed provide means to infer of the distribution of particles in the inertial range. Similarly to other situations where diffusiophoresis is at play (see e.g. Raynal et al. 2018; Basset et al. 2022), the effective diffusivity of particles can be interpreted as an effective coarse-grained compressibility that is responsible for particle clustering. Specifically, we demonstrate in this work that the statistics of the coarse-grained particle density $\langle \rho _{p} \rangle _{\ell }$ at a given inertial-range scale $\ell \gg \eta$ depend solely on the scale-dependent Péclet number $\textit {Pr}_{\ell }$. Furthermore, these diffusive models predict that the p.d.f.s of $\langle \rho _{p} \rangle _{\ell }$ exhibit algebraic tails at small values and allow for the characterisation of the associated exponent as a function of $\textit {Pr}_{\ell }$. For large masses, the models predict a superexponential behaviour that is also well reproduced by our direct numerical simulations. However, statistics that span different scales, such as the distribution of voids, display more intricate dependencies. Nonetheless, we find that the probability distribution of void volumes follows a power law with exponent steeper than $-2$ at intermediate values, transitioning to a log-normal tail at larger values. Our direct numerical simulations demonstrate a reasonably good agreement with this prediction, emphasising the need to revisit previous work on void statistics in the light of these potential behaviours.

The introduction of space-dependent diffusions in this study presents a novel framework for incorporating inertial particles into models or LES of turbulent flows. The coarse-grained particle density can be effectively approximated using diffusion equations derived from spatial averaging, with a fluctuating diffusion coefficient determined by the local turbulent dissipation rate. To test, calibrate and validate this approach, further numerical simulations that integrate the effective advection-diffusion equations at various coarse-graining resolutions are necessary. Although beyond the scope of this work, this perspective holds promise for future work.