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 (Reference Homann, Dreher and Grauer2007). 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 Reference Frisch1995). 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 (Reference Frisch1995), § 8.6.5). Despite this, such an approximation is still useful for estimating moderate-order statistics, well beyond the central-limit approximation. Using Kolmogorov (Reference Kolmogorov1962) refined similarity hypothesis, and recent confirmations by Lawson et al. (Reference Lawson, Bodenschatz, Knutsen, Dawson and Worth2019), 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 (Reference Sreenivasan and Kailasnath1993). 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. (Reference Saw, Debue, Kuzzay, Daviaud and Dubrulle2018) 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 (Reference Goto and Vassilicos2008) (see also Coleman & Vassilicos (Reference Coleman and Vassilicos2009)). 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 (Reference Oka and Goto2021) 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.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 (Reference Toschi and Bodenschatz2009), 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 (Reference Hill2002b) 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 Reference Borgas1993; Sawford et al. Reference Sawford, Yeung, Borgas, Vedula, La Porta, Crawford and Bodenschatz2003; Biferale et al. Reference Biferale, Boffetta, Celani, Devenish, Lanotte and Toschi2004). 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 (Reference Gotoh and Fukayama2001), Bec et al. (Reference Bec, Biferale, Boffetta, Celani, Cencini, Lanotte, Musacchio and Toschi2006) and Yeung et al. (Reference Yeung, Pope, Lamorgese and Donzis2006), 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 Reference Obukhov and Yaglom1951; Hill & Wilczak Reference Hill and Wilczak1995)

(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 (Reference Obukhov and Yaglom1951). 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. (Reference Xu, Ouellette, Vincenzi and Bodenschatz2007) and display a power law with an even steeper exponent close to $-1$.

This behaviour extends beyond the transition scale introduced by Hill (Reference Hill2002a), 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 (Reference Nelkin1990), 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 (Reference Hill2002a). 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. (Reference Bec, Biferale, Boffetta, Celani, Cencini, Lanotte, Musacchio and Toschi2006) 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 Reference Wilkinson, Mehlig and Bezuglyy2006), a phenomenon also coined the sling effect (Falkovich et al. Reference Falkovich, Fouxon and Stepanov2002). 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 Reference Yeung, Sreenivasan and Pope2018), consistent with the behaviour of $\varDelta _a$. However, Buaria et al. (Reference Buaria, Pumir, Bodenschatz and Yeung2019) 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. Reference Bec, Biferale, Boffetta, Celani, Cencini, Lanotte, Musacchio and Toschi2006). Gorokhovski & Zamansky (Reference Gorokhovski and Zamansky2018) 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 (Reference Gorokhovski and Zamansky2018), 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. (Reference Fevrier, Simonin and Squires2005). 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. (Reference Fevrier, Simonin and Squires2005) 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.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 (Reference Bec and Chétrite2007). 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 (Reference Reade and Collins2000) that seems confirmed by the experimental measurements of Petersen, Baker & Coletti (Reference Petersen, Baker and Coletti2019) (see also Brandt & Coletti (Reference Brandt and Coletti2022)). The scaling that we observe also significantly deviates from the prediction of Balkovsky et al. (Reference Balkovsky, Falkovich and Fouxon2001) (see also Falkovich, Fouxon & Stepanov (Reference Falkovich, Fouxon and Stepanov2003)), 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. (Reference Bec, Cencini and Hillerbrand2007b) 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. (Reference Bec, Biferale, Cencini, Lanotte, Musacchio and Toschi2007a) 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. (Reference Bragg, Ireland and Collins2015) and Ariki et al. (Reference Ariki, Yoshida, Matsuda and Yoshimatsu2018) 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 (Reference Bec and Chétrite2007). 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 Reference Gaite2005), 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. (Reference Boffetta, De Lillo and Gamba2004) 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 (Reference Goto and Vassilicos2006) 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 (Reference Yoshimoto and Goto2007) 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. (Reference Sumbekova, Cartellier, Aliseda and Bourgoin2017) 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 (Reference Yoshimoto and Goto2007) and Sumbekova et al. (Reference Sumbekova, Cartellier, Aliseda and Bourgoin2017).

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 Reference Cristelli, Batty and Pietronero2012). 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. Reference De Marzo, Gabrielli, Zaccaria and Pietronero2021).