Dynamics of small heavy particles in homogeneous turbulence: a Lagrangian experimental study

Abstract We investigate the behaviour of microscopic heavy particles settling in homogeneous air turbulence. The regimes are relevant to the airborne transport of dust and droplets: the Taylor-microscale Reynolds number is $Re_\lambda = 289\text {--}462$, the Kolmogorov-scale Stokes number is $St = 1.2\text {--}13$ and the Kolmogorov acceleration is comparable to the gravitational acceleration (i.e. the Froude number $Fr = O(1)$). We use high-speed laser imaging to track the particles and simultaneously characterize the air velocity field, resolving all relevant spatio-temporal scales. The role of the flow sampled by the particles is spotlighted. In the present range of parameters, the particle settling velocity is enhanced proportionally to the velocity scale of the turbulence. Both gravity and inertia reduce the velocity fluctuations of the particles compared to the fluid; while they have competing effect on the particle acceleration, through the crossing trajectories and inertial filtering mechanisms, respectively. The preferential sampling of high-strain/low-vorticity regions is measurable, but its impact on the global statistics is moderate. The inertial particles have large relative velocity at small separations, which increases their pair dispersion; however, gravity offsets this effect by causing them to experience fluid velocities that decorrelate faster in time compared to tracers. Based on the observations, we derive an analytical model to predict the particle velocity and acceleration variances for arbitrary $St$, $Fr$ and $Re_\lambda$. This agrees well with the present observations and previous simulations and captures the respective effects of inertia and gravity, both of which play crucial roles in the transport.


Objective
The objective of the present study is to explore the effects of inertia and gravity on the motion of heavy particles in homogeneous turbulence. Given its relevance to countless aspects of natural, industrial and medical settings, the topic has been studied in depth and a large body of literature is summarized by excellent reviews (Balachandar & Eaton 2010;Gustavsson & Mehlig 2016). Reaching a predictive understanding of the transport of small airborne particles is even more crucial in the current pandemic (Mittal, Ni & Seo 2020). Still, solid evidence concerning important aspects of the particle transport (e.g. fall speed, velocity fluctuations, acceleration and dispersion) has remained elusive. The reasons include the difficulty of carrying out detailed measurements resolving all important spatio-temporal scales, the stringent hypotheses of theoretical and numerical models and the scarcity of one-to-one comparisons between computations and experiments. In particular, it is noteworthy that the majority of numerical studies have focused on the zero-gravity case. If this allows us to isolate the effect of inertia, it also impedes the direct validation in the laboratory. Indeed, several behaviours of heavy particles in turbulence were theorized and simulated for decades, e.g. the oversampling of regions of high strain and downward velocity fluctuations (Maxey 1987;Squires & Eaton 1991b;Wang & Maxey 1993); nevertheless, they were only recently demonstrated and quantified by experiments (Petersen, Baker & Coletti 2019).
Here, we consider the case of turbulent air laden with solid particles much smaller than the Kolmogorov scale η. We focus on the range of Stokes number (based on the Kolmogorov time scale) St ≈ 1-13, while the acceleration scale of the turbulence is of the same order as the gravitational acceleration. These conditions are especially relevant to the transport of dust and droplets in the atmosphere. We perform time-resolved imaging of both the fluid and the particle motion, resolving virtually all scales at play. The measurements of particle velocity, acceleration and relative motion highlight the influence of the turbulence and the competing effects of inertia and gravity, demonstrating how in the present regime both effects are crucial for the fate and transport of the dispersed phase. The analysis culminates with an analytical model that builds on the classic framework put forward by Csanady (1963) to estimate the motion of the particles from the properties of the fluid they sample. The model is shown to agree well with the present observations, as well as with data in the literature. The rest of the paper is organized as follows: in § 1.2 we summarize the theoretical background and key existing results relevant to the present study; in § 2 we describe the experimental apparatus, and the measurement and processing procedure; in § 3 we present and discuss the data; in § 4 we introduce the analytical model and compare it with the observations, before drawing conclusions in § 5.

Background
We briefly review some fundamental relations and concepts that will help interpret and model the behaviour of the dispersed phase. We consider particles much denser than the fluid and sufficiently small compared to any flow scale. We indicate properties of the particles, fluid and fluid at the particle location with subscripts 'p', 'f ' and 'fp', respectively. If the particle Reynolds number Re p (based on the particle diameter d p and a slip velocity from the fluid u s ) is sufficiently small, drag and gravity are the only forces usually retained in the equation of motion (Maxey & Riley 1983) (1.1) Here u p is the particle velocity, a p = du p /dt is the particle acceleration, u fp is the fluid velocity at the particle location (which we will refer to as the sampled-fluid velocity), τ p is the particle response time and g is the gravitational acceleration along the unit vectorê y . We normalize using the Kolmogorov scales for time, velocity and acceleration: τ η = (ν/ε) 1/2 , u η = (νε) 1/4 and a η = u η /τ η = (ε 3 /ν) 1/4 , respectively, where ν is the kinematic viscosity and ε is the dissipation rate. This yields The Stokes number St = τ p /τ η and the Froude number Fr = a η /g can be combined in the settling parameter Sv = St/Fr = (τ p g)/u η . We define the slip velocity as u s = u p − u fp , i.e. the particle velocity relative to the surrounding flow. (The opposite sign convention is also used in the literature, defining the slip velocity as the fluid velocity seen by the particle.) Denoting averaged quantities with angled brackets and fluctuating ones with a prime, the Reynolds decomposition reads a p = a p + a p = τ p −1 u fp + u fp − u p − u p − gê y . In equilibrium conditions ( a p = 0) and still fluid ( u fp = 0) the particles settle at a terminal velocity u p,y = u s,y = −τ p g. Turbulence can either increase or decrease the fall speed through different mechanisms (Nielsen 1993;Good, Gerashchenko & Warhaft 2012). Perhaps the best known among those is preferential sweeping, by which inertial particles oversample downward sides of turbulent eddies, leading to a net increase in settling velocity, especially up to St = O(1) (Maxey 1987;Wang & Maxey 1993). Consistently with this view, Good et al. (2014) probed the parameter space and found that settling enhancement was maximum for St = O(1) and Sv = O(1). From (1.4), the modification of the mean vertical velocity of the particles can be expressed as u y = u p,y + τ p g = u p,y − u s,y = u fp,y (negative for settling enhancement). That is, the settling rate modification is determined by the vertical fluid velocity sampled by the particles. It is debated which velocity scale governs the phenomenon, with various studies indicating settling enhancement proportional with the root-mean-square (r.m.s.) of the fluid velocity fluctuations u rms ≡ (u f ) 2 1/2 : typically u y ≈ −0.2u rms (Aliseda et al. 2002;Yang & Shy 2005;Huck et al. 2018). This enhancement was recently associated with a version of the sweep-stick mechanism in which particles oversample regions of small Lagrangian acceleration (Falkinhoff et al. 2020). Recently Momenifar & Bragg (2020) showed by theoretical arguments and numerical simulations that the enhancement is in fact governed by a range of velocity scales, whose width increases with St. These effects lead to a settling enhancement scaling with u rms for St = O(1) and Sv = O(1), but dependent on particle inertia outside this range as observed by Rosa & Pozorski (2017) and Tom & Bragg (2019). We note that, in the present study, we investigate the range in which settling enhancement scales with u rms , as evidenced in § 3.1.
Averaging over the square of (1.5) leads to (1.6) which relates the particle velocity and acceleration variance to the sampled-fluid velocity and slip velocity variance. Of particular relevance to such a relation is the framework developed by Csanady (1963) (building on earlier work by Tchen 1947;Hinze 1975), who proposed a link between the Lagrangian spectrum of the particle velocity (E p ) and the one of the sampled-fluid velocity (E fp ) through a response function H 2 where ω is the angular frequency in the Lagrangian frame of reference. The velocity and acceleration variances can be obtained from the spectra as The response function can be modelled from (1.1), by taking the Fourier transform of the particle velocity and sampled-fluid velocity (Csanady 1963) (1.10) Zhang, Legendre & Zamansky (2019) derived a more complete form of the response function including unsteady forces, which, however, are expected to be small for microscopic heavy particles. The dependence of the response function on the particle response time illustrates the particle inability to respond to fluctuations with frequencies greater than O(τ −1 p ). This behaviour, often termed inertial filtering, was clearly demonstrated in situations where gravity is absent or negligible (Ayyalasomayajula et al. 2006;Bec et al. 2006). Deutsch & Simonin (1991) made the important distinction between the spectra of the flow and the spectra of the flow sampled by the particles. In presence of gravity, already Yudine (1959) realized the importance of the drift through turbulent eddies, with particles experiencing fast-changing flow conditions: this crossing-trajectories effect leads to faster decorrelation of the particle motion (Squires & Eaton 1991a;Elghobashi & Truesdell 1992;Wang & Stock 1993). Consequently, gravitational drift enhances particle acceleration compared to zero-gravity conditions, counteracting the effect of inertial filtering (Ireland, Bragg & Collins 2016a,b). For weakly inertial particles (St < 1), preferential sampling of high-strain, low-vorticity regions may also contribute to increasing particle acceleration, as the fluid acceleration is higher in the strain-dominated regions (Bec et al. 2006).

Experimental apparatus
Experiments are performed in a chamber where a region of homogeneous anisotropic air turbulence is formed by two facing jet arrays. The facility was introduced and qualified in detail in Carter et al. (2016) and Carter & Coletti (2017; here, we only give a brief description for completeness. The chamber measures 2.4 × 2 × 1.1 m 3 in the x, y and z μ on (s) u rms (m s −1 ) u rms,x /u rms,y L (m) T E (s) T L (s) η (mm) τ η (ms) Re λ 0. directions (where x is aligned with the jet axis and y is in vertical direction) and has acrylic walls for optical access. Each jet array consists of 128 quasi-synthetic jets, individually operated according to a sequence proposed by Variano & Cowen (2008). The region of homogeneous turbulence (with negligible mean flow and shear) measures approximately 0.5 × 0.7 × 0.4 m 3 . The Reynolds number can be tuned by adjusting the average firing time of the jets (μ on ). In the present study we operate the jets in two modes, and the main turbulence properties for each are reported in table 1. For both cases, the region of homogeneity is substantially larger than the integral length scale of the turbulence, allowing for a natural inter-scale energy cascade without major influence of the boundary conditions. The materials and procedure followed to investigate the inertial particle transport are similar to those in Petersen et al. (2019). The particles are fed into the chamber via a 3 m vertical chute connected to the top of the chamber, being released at a steady rate using an AccuRate dry material feeder. We use three sizes of soda-lime glass beads (density ρ p = 2500 kg m −3 ), with mean diameters d p = 32 ± 7 μm, 52 ± 6 μm and 96 ± 11 μm. Particle response times are evaluated using the Schiller & Naumann correction (Clift, Grace & Weber 2005) where μ is the air dynamic viscosity and Re p,0 = d p τ p g/ν is the particle Reynolds number based on the still-air settling velocity τ p g. Iterative evaluation of (2.1) yields response times of τ p = 7.4, 17 and 47 ms respectively. The particles are expected to approach terminal velocity over a distance of the order of τ 2 p g. For the largest particles this is 2 cm, an order of magnitude smaller than the extent of the homogeneous region traversed before reaching the measurement field of view. Table 2 reports the main non-dimensional parameters for the five experimental cases obtained combining the different particle types and turbulence forcing. We focus on St = O(1) and Sv = O(1), with one case of larger St and Sv. For the particle volume fraction and mass fraction in this study (of order 10 −4 and 10 −7 , respectively), the flow properties are not expected to be modified by the loading (Petersen et al. 2019).

Measurement approach
Imaging is performed in the x-y symmetry plane at the centre of the region of homogeneous turbulence. The flow is seeded with 1-2 μm DEHS (di-ethyl-hexyl-sebacat) droplets which faithfully follow the flow. The flow is illuminated using a Nd:YLF single-pulse laser (Photonics, 30 mJ pulse −1 ) synchronized with a VEO640 camera mounting a 200 mm Nikon lens. An aperture number f # = 4 gives a thickness of the focal plane of 1.5 mm. The active portion of the camera sensor, and therefore the size of the  Table 3. Acquisition parameters for the two turbulence forcing cases; f aq and T aq represent the acquisition frequency and length of acquisition, respectively.
field of view (FOV), depends on the acquisition frequency. The latter is optimized to give a time separation between consecutive images of at most 0.1τ η , yielding the resolutions reported in table 3. For both Re λ , the FOV is much smaller than the region of homogeneous turbulence. For each case we record 10 separate runs, with total duration of the recordings of around 40 integral time scales. Data are processed using the procedure detailed in Petersen et al. (2019). Raw images are separated in particle-only images and tracer-only images, distinguishing DEHS droplets and glass beads based on brightness and size. Particle image velocimetry (PIV) is performed on the tracer-only images. We use an initial interrogation window of 64times64 pixels, refined to 32 × 32 pixels with 50 % overlap, for a vector spacing of 0.9 mm (2.9η and 3.7η for the Re λ = 289 and 462 cases respectively) that resolves the fine scales of turbulence (Worth, Nickels & Swaminathan 2010). Particle tracking velocimetry (PTV) is performed on the particle-only images, following the cross-correlation approach. The fluid velocity is evaluated at the particle locations using weighted linear interpolation of the four neighbouring velocity vectors. A comparison with cubic and spline interpolation shows no significant difference, but linear interpolation substantially decreases the computation time. Particle velocities and accelerations are determined from the particle position using convolution with the first and second derivatives of a Gaussian kernel, respectively. The width of the latter is chosen as 0.5τ η and 0.4τ η for the Re λ = 289 and 462 cases, respectively, following the procedure established for tracers in Voth et al. (2002) and Mordant, Lévêque & Pinton (2004) and applied to inertial particles in Gerashchenko et al. (2008), Nemes et al. (2017) and Ebrahimian, Sanders & Ghaemi (2019b,a).
The number of samples used to calculate particle statistics ranges between 0.8 × 10 6 and 5.6 × 10 6 for the different cases. Uncertainty in the statistics is affected by both random uncertainty, due to the finite sample size, and bias uncertainty, due to systematic errors in estimating the particle centroid and the local fluid velocity. Baker & Coletti (2021), who followed a similar time-resolved PIV/PTV approach for particle-laden turbulent boundary layers, showed that the bias uncertainty was negligible compared to the random uncertainty. Compared to their study, the present particles are much smaller with respect to the flow scales, and therefore the error associated with interpolating the fluid velocity at the particle location is correspondingly smaller. The centroid location uncertainty  Table 4. Random errors of the velocity and acceleration statistics based on the standard deviation of the last 20 % of data for the case with the lowest number of samples. (investigated for these same particles in Petersen et al. 2019) is an order of magnitude smaller than the typical particle displacement as in Baker & Coletti (2021), and therefore its impact is also deemed negligible compared to the random uncertainty. The latter is estimated based on the standard deviation of the last 20 % of data (Ebrahimian et al. 2019a), and is reported in table 4 for the main quantities.

Velocity
We begin by considering the statistics of the sampled fluid, focusing on the vertical component u fp,y which is most relevant to the settling process. Figure 1(a) shows the mean vertical sampled-fluid velocity for all particles, as well as the velocity conditioned on upward-/downward-moving particles (i.e. particles with positive/negative u fp,y ). The results are approximately independent of Sv and consistent with the scaling u fp,y /u η = C, where the constant C is the mean for the respective sets of particles. Given the expected dependence of the vertical velocities with St and Sv, here and in the following figure we shade the area outside of the investigated range to emphasize that the trends should not be extrapolated. According to previous studies mentioned in § 1, we expect u fp,y ≈ −0.2u rms , which over the present range of Re λ implies u fp,y ≈ −2u η . Indeed, for the unconditional average over all particles we find C ≈ −2. The results for the downward-moving and upward-moving subsets are consistent with C ≈ −5 and C ≈ 7, respectively. Thus, the upward-/downward-moving particles sample fluid regions with relatively large upward/downward velocity fluctuations, and these are of similar magnitude for both subsets. Therefore, the net settling enhancement stems from the downward particles being more numerous, not from their association with stronger downward events. Figure 1(b) shows the normalized slip velocity, the dashed line indicating u s,y /u η = −Sv, or u s,y = −τ p g, as theorized by Wang & Maxey (1993). This corresponds well to the measured data, except for the upward-moving particles of largest Sv, which have a significantly smaller slip velocity. This is likely related to the assumption a p = 0. The latter is strictly valid only for the ensemble of all particles, and not necessarily for specific subsets. In absence of this assumption, (1.4) predicts a smaller slip velocity for mean downward particle acceleration. We then consider the particle vertical mean velocity u p,y . Since in the considered regimes we have approximately u fp,y /u η = C and u s,y /u η = −Sv, we expect u p,y /u η = C − Sv. Figure 2(a) supports this scaling for the ensemble of all particles. Similarly, normalizing by the still-air terminal velocity leads to u p,y /(τ p g) = C/Sv − 1. This is confirmed in figure 2(b), where the scaling is shown to hold also for upward and downward-moving particles (with the respective values of the constant). We again stress that this cannot be extrapolated ad libitum: in the limit of both vanishing and infinite particle inertia, we expect u p,y /(τ p g) = −1.
The probability density function (p.d.f.) of the normalized vertical particle velocity, u p,y /u η , is shown in figure 3(a) for three selected cases, the other cases sharing the same trends. For increasing Sv, the distributions shifts to more negative values, as expected, and the distributions become positively skewed. This is in contrast with the findings of Baker et al. (2017), who used point-particle simulations. The disagreement for the heavier particles is not surprising, as the assumptions behind such simulations become questionable with increasing Re p . Figure 3(b) reports the p.d.f. of the absolute value of the particle velocity fluctuations, |u p,y |/u η , distinguishing between upward-and downward-moving particles. We observe no difference between both subsets, and therefore in the following we do not consider them separately. All Sv cases collapse well on a Gaussian distribution when normalized by the Kolmogorov velocity. This confirms the dominant role of the fluid fluctuations in determining the particle velocity fluctuations, for a wide range of response times and fall speeds.
We then investigate the scaling of the particle velocity fluctuations by considering the variance of u p = u fp + u s (3.1) Figure 4 displays the vertical components of the four terms in (3.1) for the various cases, normalized by Kolmogorov scaling. The horizontal components, not shown, behave similarly. The particle velocity variance (u p,y ) 2 is smaller than but comparable to the variance of the sampled-fluid velocity (u fp,y ) 2 (figure 4a), as also reported by Ireland et al. (2016a) in zero-gravity simulations. This confirms that, in the present range of parameters, the fluctuating energy of the particles is driven by the turbulent kinetic energy. Figure 4(b) indicates that the normalized slip velocity variance (u s,y ) 2 /u 2 η varies linearly with St. This is consistent with the scaling u s /u η ∝ St 1/2 in Balachandar (2009), derived for the present range (τ η < τ p < T L ) but in the absence of gravity. The covariance u fp,y u s,y also varies linearly with St and approximately equals − (u s,y ) 2 . Therefore, from (3.1) we have (u p,y ) 2 ≈ (u fp,y ) 2 − (u s,y ) 2 . As (u s,y ) 2 grows with St, we retrieve the influence of inertial filtering: heavier particles exhibit weaker velocity fluctuations with respect to the sampled-fluid fluctuations. We conclude this section by comparing the sampled-fluid and particle velocity variances against the fluid velocity variance, again focusing on the vertical components (figure 5). This allows us to quantitatively compare the fluctuating energy of the dispersed and carrier phase. The zero-gravity simulations of Ireland et al. (2016a) indicated that the particle velocity fluctuations can exceed the fluid fluctuations for St < 1, due to preferential sampling of energetic flow regions. In the present case, on the other hand, both inertia and gravity concur to reduce the particle fluctuating energy below the turbulent kinetic energy of the fluid, in agreement with the algebraic model of Wang & Stock (1993) and the results of Good et al. (2014) in similar ranges of St and Sv. In particular, figure 5 shows that the fluctuating energy of the fluid sampled by the particles is somewhat lower than (although comparable to) the unconditional turbulent kinetic energy. This implies that, as far as this observable is concerned, the effect of preferential sampling is relatively weak. The latter is indeed offset by gravitational drift, which reduces the particle ability to follow energetic fluid structures, hence (u fp,y ) 2 < (u f ,y ) 2 . This picture will be confirmed later, when analysing the fast decorrelation of the sampled-fluid velocity (see § 4.2). Inertial filtering further reduces (u p,y ) 2 with respect to (u fp,y ) 2 , and therefore the former is 10-30 % lower than (u f ,y ) 2 .

Acceleration
We present results for the vertical components of the particle acceleration, a p,y , the horizontal components behaving similarly. In addition, we consider the temporal derivative of the sampled-fluid velocity, du fp,y /dt. The p.d.f.s of the latter are presented in figure 6(a), showing similar distributions for the cases with the same Re λ . Indeed, for vanishing inertia du fp,y /dt equals the Eulerian acceleration which, in homogeneous turbulence, is a function of Re λ only (Hill 2002;Sawford et al. 2003). For a given Re λ , the distributions of du fp,y /dt become wider with Sv, in agreement with the simulations of Ireland et al. (2016b). This is a manifestation of the crossing-trajectories effect: for higher settling rate, the particles experience rapid changes of the sampled-fluid environment. The distributions of a p,y (figure 6b) are much narrower than the corresponding distributions of du fp,y /dt, due to inertial filtering. This effect becomes stronger with larger particle inertia, which in the present case implies an increase of both St and Sv.
The variances of du fp,y /dt and a p,y are quantified in figure 7. As observed above, the variance of du fp,y /dt is dominated by Re λ and increases with Sv. For tracers in homogeneous isotropic turbulence, the normalized acceleration variance of the fluid can be approximated as , hence for 917 A47-10 Re λ = 289 and 462 we expect a 0 = 3.6 and 4.0, respectively. Ireland et al. (2016b) showed a monotonic increase of (du fp,y /dt) 2 /u 2 η roughly with Sv 3/2 , independent of Re λ . In the present study, the variance of du fp,y /dt is larger than for tracers, but not as large as in Ireland et al. (2016b) and still dependent on Re λ . The difference with Ireland et al. (2016b) is reflected in the Froude number: in their study Fr = 0.052, while here Fr = 0.8 and 1.9, for Re λ = 289 and 462 respectively. For Fr = O(1), gravitational and Kolmogorov acceleration are of the same order of magnitude, i.e. fluid turbulence and gravity are expected to have comparable influences. The variance of a p,y (figure 7b) decreases approximately linearly with Sv in the present range. It is significantly smaller than the variance of du fp,y /dt and also smaller than a 0 , due to inertial filtering.

Preferential sampling
To quantify the extent and influence of preferential sampling, we discriminate between rotation-dominated and strain-dominated fluid regions using the second invariant of the velocity gradient tensor (Hunt, Wray & Moin 1988 where ω 2 = 2tr(Ω 2 ) is the enstrophy, s 2 = tr(S 2 ) is the squared strain rate and S and Ω are the symmetric and anti-symmetric part of the velocity gradient tensor, respectively. Due to the planar nature of the measurements, we can only consider the in-plane components of the velocity gradient tensor, which are not sufficient to fully describe the flow topology (Perry & Chong 1994). Still, especially in homogeneous turbulence, the in-plane part of the tensor provides important physical insight into the properties of high-enstrophy and high-strain structures (Cardesa et al. 2013) and captures the fundamental small-scale features (Fiscaletti, Ganapathisubramani & Elsinga 2015;Carter & Coletti 2018). Figure 8 shows the p.d.f. of Q evaluated at the particle locations, compared to the unconditioned fluid. For the smaller St considered, the particle-conditioned distributions display the expected under-sampling of rotation-dominated regions (Q > 0), which, however, becomes progressively weaker as St increases above unity. The effect of the small-scale features of the sampled fluid on the particle motion is depicted in figure 9, which plots the particle acceleration variance conditioned on strain rate (figure 9a) and enstrophy (figure 9b). In both cases, larger levels of small-scale turbulence activity correspond to stronger accelerations. Although the correlation with s 2 is somewhat stronger than with ω 2 , the similarity of both plots is consistent with the view that high-strain and high-enstrophy events are often concurrent (Worth & Nickels 2011;Yeung, Donzis & Sreenivasan 2012;Carter & Coletti 2018). Due to inertial filtering, the impact of the sampled-fluid topology decreases steeply with St. This is clearly shown in figure 9(c), displaying the variance of the particle acceleration conditioned on the sign of Q: particles in strain-dominated regions do display larger accelerations (Bec et al. 2006), but the effect becomes unmeasurable for St = 4.7 and larger.

Structure functions and pair dispersion
In this section we consider two-particle statistics, starting with the second-order Eulerian velocity structure function S 2 (r) = (u (x) − u (x + r)) 2 , where r is the separation vector. For fluid tracers in homogeneous isotropic turbulence, this scales as r 2 in the dissipative range, plateaus to the fluid velocity variance at large-scale separations, and follows the scaling predicted by Kolmogorov (1941) in the inertial range S 2 (r) = C 2 (εr) 2/3 or S 2 /u 2 η = C 2 (r/η) 2/3 , where the subscripts and ⊥ denote velocity components longitudinal and transverse to the separation vector, respectively, and C 2 ≈ 2 (Saddoughi & Veeravalli 1994). For non-tracer particles, inertia and gravity modify these trends. Simulations and experiments indicate that, when gravitational effects are negligible or absent, inertia leads to greater relative velocities between nearby particles, such that the structure function increasingly deviates from the r 2 scaling at small separations ( indicated a strong reduction of relative particle velocity at all separations. They attributed this to the decorrelation of the sampled-fluid velocity along the particle trajectories (which we shall confirm later), hindering the path-history effect and in turn decreasing the relative velocities. To our best knowledge, no previous experimental observation could verify this latter point. The longitudinal structure functions are presented in figure 10, the transverse components (not shown) showing analogous trends. Results are compared to the measured Eulerian structure function for tracers, which follows both r 2 scaling in the dissipative range and the scaling from (3.3) in the inertial range. Due to the finite laser sheet thickness (≈6η) we expect an overestimation of the relative velocity over the dissipative range (Dou et al. 2018), which, however, may not overwhelm the trend. At small separations (r 20η), the structure functions of the particles deviate from tracers with increasing St, confirming previous findings. In the inertial range, S 2 (r) roughly follows the r 2/3 scaling, but the values are significantly lower than for tracers. The gap persists at large scales, consistent with the fact that the inertial particle fluctuating energy (to which the structure function asymptotes for large r) is lower compared to the fluid, see figure 5. Although the competing effects of inertia and gravity cannot be separated here, these results appear to confirm the observation of Ireland et al. (2016b) that gravity reduces the relative particle velocities at all scales.  Kolmogorov (1941) for the inertial range, respectively. We then turn to particle pair separations as a function of time, r(t). In homogeneous isotropic turbulence, for tracer pairs with an initial separation r 0 in the inertial range, we expect the mean square separation to follow the ballistic scaling proposed by Batchelor (1950 where t B = (r 2 0 /ε) 1/3 is the characteristic time scale of an eddy of size r 0 . Recently, a more general time scale t 0 = S 2 (r 0 )/ε has been proposed by Bitane, Homann & Bec (2012). For t B t T L (or t 0 t T L ), the dispersion does not depend on r 0 and is expected to follow the Richardson-Obukhov scaling r(t) 2 = gεt 3 , where g ≈ 0.5 (Salazar & Collins 2012). For t in the dissipative range, particle inertia enhances pair dispersion at small times, due to the large relative velocities at small separation (Bec et al. 2010;Gibert, Xu & Bodenschatz 2010); while for larger t, the inertial filtering and path-history effects reduce pair dispersion compared to tracers (Bragg, Ireland & Collins 2016). To the best of our knowledge, the only previous investigations on the effect of gravity on pair dispersion are the numerical studies by Chang, Malec & Shaw (2015) and Dhariwal & Bragg (2019), who mostly focused on bi-dispersed particles sets.
The mean square separation for our inertial particles is presented in figure 11. Due to the nature of the measurements, the results are biased by the constraint that the trajectories cannot separate more than the laser sheet thickness in z. One way to account for this is to consider that the right-hand side in (3.5) is the geometric average of the structure function components at r 0 , i.e. (r(t) − r 0 ) 2 = (S 2 (r 0 ) + 2S 2⊥ (r 0 ))t 2 . Setting to zero the out-of-plane velocity, we can write the mean square separation for tracers as 7 3 C 2 (εr 0 ) 2/3 t 2 , which is plotted as a black solid line for reference. The effect of St (and indirectly of Sv) is represented in figure 11(a), for trajectories with initial separation 3 < r 0 /η < 4. With increasing particle inertia, the mean square separation grows. Still, all curves are generally below the expectation for tracers, illustrating the competing effect of inertia and gravity. The normalized form plotted in figure 11(b) illustrates the effect of the initial separation for the case St = 13: with increasing initial separations, the impact of the large relative velocities is weakened, and the curves tend to collapse on each other (which is also the case for smaller St, not shown). This confirms the strong influence of the uncorrelated motion of nearby particles, already highlighted by the structure functions. The plot also emphasizes how the particles follow the Batchelor regime up to t = O(τ η ), which corresponds as well to t = O(t B ) and t = O(t 0 ). Around t = O(τ η ) we observe a  Figure 11. Mean square separation of particle pairs with initial separation r 0 . Mean square separation of pairs with 3 < r 0 /η < 4 for all cases (a). Mean square separation of pairs with 3 < r 0 /η < 4 to 9 < r 0 /η < 10 (increasing with darker shades) for the St = 13 case (b). The solid thin lines indicate the scaling in (3.5). mild slowing down of the ballistic separation, also observed for tracers by Bourgoin (2015), after which the slope increases. However, given the bias due to the shape of the illuminated volume, caution should be exerted when interpreting the behaviour for relatively long times. Here, we just note that, for particles in zero gravity, the transition out of the ballistic regime is expected to ensue at times O(τ p ) (Bragg et al. 2016), i.e. much later than what we observe. The difference appears to be rooted in the effect of gravity and could be associated with the time scale of eddy crossing by the falling particles, similarly to recent findings for rising bubbles (Mathai et al. 2018). Further research is warranted on this point.

Modelling the velocity and acceleration variances
Here, we leverage and expand on the framework of Csanady (1963) to obtain expressions for the particle velocity and acceleration variances. To summarize the detailed explanation that follows, we integrate the Lagrangian particle velocity spectrum, which is obtained from the sampled-fluid velocity spectrum modulated by a response function. To extend the frequency range of the spectrum, we Fourier transform the Sawford (1991) expression of the velocity autocorrelation. As this was originally derived for the unconditional fluid, we substitute in it the Lagrangian time scales associated with the sampled fluid.

Lagrangian spectrum of particle velocity
For frequencies in the inertial range π/T L ω π/τ η , the Lagrangian spectrum for the fluid flow is expected to scale as E ∝ εω −2 (Tennekes & Lumley 1972;Yeung 2001). The spectrum can also be derived from the velocity autocorrelation as the two form a Fourier transform pair. In simple approaches considering only the inertial and large-scale dynamics, the auto-correlation function is well approximated by R(t) = (u ) 2 exp(−t/T L ) (Hinze 1975;Mordant et al. 2001;Zhang et al. 2019), valid for ω π/τ η (or τ η t) and corresponding to E = (u ) 2 T L /(1 + (ωT L ) 2 ). For small time separations, the autocorrelation deviates from the exponential, tending to a horizontal asymptote at t = 0 with a curvature proportional to the acceleration variance (Mordant et al. 2004). Sawford (1991) proposed an alternative formulation valid for all t, using a second time scale T 2 related to the finite acceleration variance Fourier transformation yields (4.2) Note that in this two-time-scale model T L no longer equals the correlation time. We use this expression to model the Lagrangian spectrum of the sampled-fluid velocity, by substituting the respective time scales T L,fp and T 2,fp . Using the response function (1.10) in combination with (1.8) and (1.9), we have . (4.4)

Time scales of the sampled-fluid velocity
Because evaluating equations (4.3) and (4.4) requires estimates for T L,fp and T 2,fp , we consider the issue of how those compare to T L and T 2 . For tracers, the Lagrangian integral time scale is related to the Eulerian integral time scale (T E ) as (Yeung 2001) where C 0 is the pre-factor in the expression of the Lagrangian velocity structure function. This depends on the turbulence Reynolds number, and based on a review of the literature Lien & D'Asaro (2002) (Ouellette et al. 2006). In addition, the numerical prefactor in (4.5) may depend on Reynolds number and large-scale properties of the turbulence (see, for example Sreenivasan 1998;Vassilicos 2015). Considering inertial particles, the Lagrangian integral time scale of the sampled fluid T L,fp is influenced by both inertia and gravity. For the range of St in the present study, however, the effect of St as derived empirically by Wang & Stock (1993) is negligible. (Alternatively, Jung, Yeo & Lee (2008) proposed an empirical expression of the St effect derived in zero gravity.) Gravity decreases T L,fp due to the crossing-trajectories effect, for which we adopt the expression proposed by Csanady (1963) (see also Pozorski & Minier 1998) where α = 1 and 4 for the time scales associated with the vertical and horizontal components, respectively. Alternatively, using (u f ) 2 /u 2 η = Re λ / √

(Hinze 1975) and
T L /τ η = 2(Re λ + 32)/( √ 15C 0 ) (Zaichik, Simonin & Alipchenkov 2003), this can be expressed as T L,fp τ η = 2(Re λ + 32) which is a function of Re λ and Sv only. For tracers, the time scale T 2 can be written as (Sawford 1991) A relation for T 2,fp is obtained from the variance of du fp /dt (4.9) Thus, we can empirically evaluate T 2,fp using the measured values for (du fp /dt) 2 and (u fp ) 2 . The values are presented in figure 12, along with the prediction for tracers in (4.8); T 2,fp is measurably smaller than T 2 , although no clear trend with St is discerned.
Note that (4.9) relates the two time scales by a Taylor time scale of the sampled flow, which can be defined as τ 2 fp,λ ≡ (u fp ) 2 / (du fp /dt) 2 (Tennekes & Lumley 1972;Pope 2000), such that τ 2 fp,λ = T L,fp T 2,fp . With the above estimates of T L,fp and T 2,fp , we use (4.1) to model the autocorrelations of the sampled-fluid velocity, R fp . We consider both horizontal and vertical components, and present the results as correlation coefficients, ρ fp , normalizing by the variance of the sampled-fluid velocity for the respective components. The modelled curves are plotted in figure 13, along with the corresponding measurements and the fluid velocity autocorrelation (which is also modelled via (4.1), using T L and T 2 ). For clarity we separate both Re λ cases, as these have different fluid time scales. The plots only extend to values of t for which the measured autocorrelations are based on at least 100 trajectories. This allows for limited time lags, still sufficient to highlight the trends. The sampled-fluid velocity decorrelates faster with Sv due to the crossing-trajectory effect. Despite quantitative differences with the measurements, (4.1) captures well this trend for short times. In figure 14 we compare the measured ρ fp against the particle velocity autocorrelation coefficients, ρ p ; the fluid velocity autocorrelation is also included for reference. The apparent trend is that ρ p decays more slowly for heavier and heavier particles, which is a consequence of inertia; while ρ fp decays more rapidly, which is a consequence of gravity.

Model predictions
We now compare the measured variances of the particle velocity and particle acceleration against the respective predictions from (4.3) and (4.4). In non-dimensional form , (4.10) .
(4.11) Red lines indicate the predictions from (4.10) and (4.11), coloured by Fr. Circles indicated measured data, coloured by Sv (see, for example figure 13 for colour coding). The '+' signs indicate case-specific predictions as discussed in the text for matching colours. Grey and black dashed lines in (a) indicate results from the model by Ayala et al. (2008) for Fr = 10 −3 and 10 2 respectively. The dashed line in (b) indicates the asymptote St −2 .
These are shown in figures 15(a) and 15(b), respectively, as a function of St. The expressions are essentially model forms of the particle response function, as they describe the variance of the particle velocity and acceleration with respect to the sampled-fluid flow. Here, and in the following, we again consider the vertical components only, the horizontal components leading to analogous conclusions. These 'generalized predictions' are compared with the measurements and with 'case-specific predictions'. The latter differ from the generalized predictions due to the difference in T 2,fp : for the generalized model we approximate T 2,fp ≈ T 2 (which is calculated for any St according to (4.8)), whereas for the case-specific predictions we use the empirically determined values for T 2,fp (figure 12). As can be seen, this simplification has a moderate impact and does not alter the trends. In the small St limit, the modelled particle velocity variance equals the variance of the sampled-fluid velocity, regardless of Fr. For finite St, the particle velocity variance (u p ) 2 / (u fp ) 2 decreases with both inertia and gravity (figure 15a). The dependence with St reflects inertial filtering, while the effect of Fr reflects the decrease of T L,fp with gravitational drift, which in turn damps (u p ) 2 (4.10). The effect of gravity becomes negligible for Fr 10. In the limits of either strong gravity (Fr 1) or large inertia (St 1), (u p ) 2 vanishes compared to (u fp ) 2 . Overall, the generalized model represents well the present data. Modelled results are qualitatively similar to results of the model presented in Ayala, Rosa & Wang (2008), the main quantitative difference resulting from the omission of a second time scale (T 2 ) in that model.
The normalized particle acceleration variance (a p ) 2 /( (u fp ) 2 τ −2 η ) (figure 15b) tends to a 0 /u 0 for small St, where u 0 ≡ a 2 f /u 2 η is a function of Re λ . For large St, it asymptotes to St −2 , or (a p ) 2 = (u fp ) 2 τ −2 p . For Fr 1 the particle acceleration variance displays a non-monotonic behaviour with St, due to the competing effects of inertia and gravity: an increase of inertia for a fixed Fr implies an increase of Sv and therefore of crossing-trajectory drift, augmenting du fp,y /dt and in turn also the particle acceleration. As St increases further, inertial filtering eventually dominates and the particle acceleration is dampened. As the sampled-fluid properties are not known a priori, (4.10) and (4.11) have limited predictive power. However, as shown in figure 5, the sampled-fluid velocity variance is marginally smaller than the fluid velocity variance, and one can approximate Hinze 1975) and substituting in (4.11) leads to where we assumed T 2,fp ≈ T 2 and T L,fp is obtained from (4.6). Because of the assumptions made, (4.12) is expected to underpredict somewhat the observations. The dependence on St, Fr and Re λ is illustrated in figure 16. For large St, the curves asymptote to Re λ St −2 / √ 15. At intermediate St and Fr 1, we retrieve the non-monotonic behaviour discussed above. This non-monotonicity was also reported by Ireland et al. (2016b) for Fr = 0.052 as well as by Parishani et al. (2015) for Fr ≈ 0.14. In general (4.12) agrees remarkably well with the trends reported by their simulations as well as with our experimental data. The present model demonstrates how the value of St that maximizes the particle acceleration increases with Fr, while the maximum shrinks and eventually disappears for Fr 0.1. This is the reason why the non-monotonic behaviour found by Ireland et al. (2016b) at Fr = 0.052 is not seen in our study at Fr = O(1), nor in previous studies at Fr ≈ 0.3 (Ayyalasomayajula et al. 2006) and Fr ≈ 30 (Volk et al. 2008). For increasing Re λ , the maximum of (a p ) 2 /a 2 η is found at larger St. This is consistent with the view that, with increasing Re λ , a wider range of scales may influence the particle motion, and therefore also particles with St > O(1) may respond strongly to the turbulence (Yoshimoto & Goto 2007).
We now use (4.11) to derive an expression for the slip velocity variance. Squaring and averaging the fluctuating part of the particle equation of motion, a p = −u s /τ p , we have (u s ) 2 = τ 2 p (a p ) 2 . Substituting in (4.11) gives where again we assumed T 2,fp ≈ T 2 and T L,fp is obtained from (4.6). Figure 17(a) verifies this prediction, comparing it with the case-specific predictions (using empirical estimates for T 2 ) and the measurements. The modelled slip velocity variance is in the ranges (u s,y ) 2 ≈ 0 in the small St limit, when particles faithfully follow the flow; to (u s,y ) 2 ≈ (u fp,y ) 2 at large St, when particles are ballistic and the slip velocity fluctuations effectively equal the velocity fluctuations of the sampled fluid. In between, the variance of the slip velocity increases with both inertia and gravity as both effects reduce the ability of the particles to follow the fluid fluctuations. The agreement with the measurements is remarkable. Also, the generalized model is consistent with the scaling derived by Balachandar (2009) for the different regimes: for the slip velocity variance, his arguments imply (u s,y ) 2 ∝ St 2 for τ p < τ η , (u s,y ) 2 ∝ St for τ η < τ p < T L , and (u s,y ) 2 ≈ constant for τ p > T L . Figure 17(b) shows that the covariance of the sampled-fluid velocity and slip velocity is approximately u fp u s ≈ − (u s ) 2 . This result is supported by the measurements (see figure 4), and can be derived from the assumptions of the model: comparing (4.13) and (4.10) implies (u p ) 2 = (u fp ) 2 − (u s ) 2 , and comparing the latter equality with the variance of the particle velocity (expressed as u p = u fp + u s ) implies u fp u s = − (u s ) 2 . Accordingly, the covariance ranges from approximately zero at small St, where there is no slip velocity, to u fp u s = (u fp ) 2 at large St, where the slip velocity equals the sampled-fluid velocity.
From the relationship u fp u s ≈ − (u s ) 2 (and substituting u s = τ p a p and/or u s = u p − u fp ) we can derive other covariances between the particle velocity, the sampled-fluid velocity, and the particle acceleration; those are especially useful to formulate stochastic models (Zamansky, Vinkovic & Gorokhovski 2013;Pope 2014). Normalizing by the respective r.m.s. values, we have expressions for the following correlation coefficients: , (4.14)  Figure 18. Modelled correlation coefficients of the vertical components of the particle velocity and sampled-fluid velocity (a) and of the particle acceleration and sampled-fluid velocity (b). Red lines indicate the predictions from (4.14) and (4.15), coloured by Fr. Circles indicated measured data, coloured by Sv (see, for example figure 13 for colour coding). The '+' signs indicate case-specific predictions as discussed in the text for matching colours.
The model predictions in (4.14) and (4.15) are presented in figure 18 for the vertical component. The correlation between the particle velocity and sampled-fluid velocity, ρ(u p u fp ), approximates unity in the small St limit as expected, and the model predicts the decorrelation due to inertia and gravity (figure 18a). The agreement with the measurements is satisfactory, although a comparison with data at lower Fr is needed to corroborate the prediction. In contrast, the particle acceleration does not correlate with the fluid velocity at small St (figure 18b), with ρ(a p u fp ) increasing with gravity and inertia. Indeed, in the ballistic limit, the slip velocity fluctuations are equal and opposite to the fluid velocity fluctuations, that is a p = −u s /τ p = u fp /τ p . The significant mismatch with the data (approximately 40 % for most of the cases) may partly be due to inherent uncertainty on the particle acceleration, but is also likely related to the simplistic assumption that drag and gravity are the only forces at play. Finally, the prediction that the particle acceleration is uncorrelated from the particle velocity (4.16) is confirmed by our measurements, from which we find ρ(a p u p ) < 0.05 for all cases.

Discussion and conclusions
We have experimentally investigated the transport of sub-Kolmogorov heavy particles in homogeneous turbulence. All relevant spatio temporal scales are resolved for the first time in a similar configuration, albeit through planar measurements. We consider ranges of St and Sv for which a rich particle-turbulence interaction is expected, including settling enhancement, inertial filtering and preferential sampling. The focus is on the respective roles of inertia and gravity, which have different and often competing influence on the particle motion. A unique feature of the present measurements is the access to the local properties of the turbulence experienced by the particles along their trajectories. The importance of the sampled-fluid properties is already clear from the mean vertical velocities. It is proposed that, in the present ranges of St and Sv, the particle settling velocity can be rationalized by assuming u fp,y /u η = C. Because, to first order, the settling enhancement u y = u fp,y , it is also proportional to the velocity scale of the turbulence, consistent with previously reported trends. The limited range of Re λ , however, cannot clarify whether this scale shall be u η , u rms , or a multi-scale quantity between those. The mean of the instantaneous slip velocity is well represented by Stokes drag with the Schiller & Neumann correction, at least for moderate St. Upward-and downward-moving particles sample fluid regions with large velocity fluctuations in those same directions, with similar magnitude in both cases. Therefore, the settling velocity is augmented due to the downward-moving particles being more numerous, not because the sampled downward fluid fluctuations are stronger than the upward ones.
The fluid fluctuations also play a dominant role in determining the particle velocity fluctuations. The variance of the particle velocity (u p ) 2 is comparable to but somewhat smaller than the sampled-fluid velocity variance (u fp ) 2 , due to inertial filtering; and the latter is slightly smaller than the fluid velocity variance (u f ) 2 , due to gravitational drift. While gravity and inertia have concurrent effects on the particle fluctuating energy, their influences on the particle accelerations are opposite to each other: the crossing-trajectories effect augments the temporal derivative of the sampled-fluid velocity, (du fp /dt) 2 , which act to enhance the particle acceleration variance; but this is offset (at least in the present range of parameters) by inertial filtering. The net result is that heavier particles display smaller r.m.s. accelerations and less intermittent acceleration p.d.f.s. The preferential sampling of high-strain/low-vorticity regions is measurable, but its global impact on the particle motion is weak.
The competing influences of inertia and gravity are on display also in the two-particle statistics. The uncorrelated component of the relative motion augments the particle velocity structure functions at small separations; while the reduced fluctuating energy of the particles (compared to tracers) has an opposite effect at inertial-range and large-scale separations. The large relative velocities of nearby particles, which increase with particle inertia, cause heavier particles to separate faster; still, the mean square separation is generally below the expectation for tracers. This is attributed to gravity causing the particles to experience fluid velocities that decorrelate faster in time, with respect to zero-gravity conditions. The inertial particles appear to transition out of the ballistic regime at earlier times compared to tracers, similarly to what recently shown for bubbles rising in homogeneous turbulence.
The planar nature of the measurements limits the full characterization of both the particle motion and the turbulent fluid flow. This is expected to affect especially the two-particle statistics in the form of a biased sampling of the trajectories, although this may not obscure the apparent trends. The three-dimensional tracking of the particles would overcome this limitation. Nevertheless, volumetric techniques are presently not capable of simultaneously capturing both phases at the required resolution: despite fast-paced advances in this area, the spatial resolution in volumetric PIV/PTV is still generally below what can be achieved by planar measurements (Discetti & Coletti 2018); nor these methods have systematically been adapted to multi-phase flows yet. The present study indeed highlights the central role of the sampled-fluid properties, and extending such an analysis to three-dimensional experiments remains a challenge.
Based on these experimental observations, we have derived an analytical model of particle velocity and acceleration inspired by the seminal work of Csanady (1963). This is based on applying a response function to the spectrum of fluid velocities experienced by the particles. To this end, we use the expression proposed by Sawford (1991) for the fluid velocity autocorrelation, in which we substitute estimates for the time scales T L and T 2 of the sampled fluid. We lack an analytical expression for the latter time scale, but we show that using the unconditional-fluid formulation (with classic estimates of the C 0 constant) has a small quantitative influence. In its basic form ((4.10) and (4.11)) the model provides the particle velocity and acceleration variances as a function of (u fp ) 2 . The model agrees generally well with the experimental observations, captures the respective effects of inertia and gravity over a wide range of the controlling parameters, and predicts correlations between particle and sampled-fluid velocities and accelerations. In particular, consistent with the arguments by Balachandar (2009), it predicts the variance of the slip velocity (u s ) 2 to scale as St 2 for τ p < τ η , as St for τ η < τ p < T L , and to plateau for τ p > T L .
Because we show that (u fp ) 2 is only slightly smaller than (u f ) 2 (at least in the present range of parameters), the model can be written in a weaker form with more predictive power by substituting (u fp ) 2 = u 2 η Re λ / √ 15 (Hinze 1975) in which T 2 /τ η is a function of Re λ only and T L,fp /τ η is a function of Re λ and Sv only, given by (4.8) and (4.7), respectively. When compared with the present experiments and recent simulations, in particular for the particle acceleration, this version of the model also agrees well with the observations and represents the complex dependency with inertia and gravity. In particular, it predicts the increase in r.m.s. particle acceleration with gravitational drift, and its non-monotonic dependence with St when Fr 1 as recently reported by Ireland et al. (2016b).
Taken together, the laboratory observations and the derived model indicate how, unless the turbulence acceleration is overwhelming (Fr 1), both inertia and gravity are key ingredients to understand the transport of heavy particles in homogeneous turbulence. This calls into question the practice of setting gravity to zero to isolate inertial effects. Such a consideration appears to be broadly applicable: recent studies on bubbles in homogeneous turbulence (Mathai et al. 2016) and heavy particles in turbulent boundary layers (Baker & Coletti 2021) reached a similar conclusion.
Funding. The present work was supported in part by the US Army Research Office, Division of Earth Materials and Processes (grant W911NF-17-1-0366), and Division of Fluid Dynamics (grant W911NF-18-1-0354).
Declaration of interests. The authors report no conflict of interest.