Hydrodynamic interactions and extreme particle clustering in turbulence

From new detailed experimental data, we found that the Radial Distribution Function (RDF) of inertial particles in turbulence grows explosively with $r^{-6}$ scaling as the collision radius is approached. We corrected a theory by Yavuz et al. (Phys. Rev. Lett. 120, 244504 (2018)) based on hydrodynamic interactions between pairs of weakly inertial particles, and demonstrate that even this corrected theory cannot explain the observed RDF behavior. We explore several alternative mechanisms for the discrepancy that were not included in the theory and show that none of them are likely the explanation, suggesting new, yet to be identified physical mechanisms are at play.

Small, inertial particles can spontaneously cluster in incompressible turbulent flows, an effect considered important for droplet collision rates in atmospheric clouds [1,2] and planetesimal formation in turbulent circumstellar disks [3]. However, even in the absence of particle inertia, hydrodynamic interactions (HI) between pairs of particles can also lead to particle clustering [4]. This behavior has been explored theoretically for inertia-free particles (St = 0, where Stokes number St defined as the particle response time τ p divided by the Kolmogorov timescale τ η ) in low Reynolds number flows [5] as well as turbulent flows [4]. These analyses show that the Radial Distribution Function (RDF) g(r), which quantifies clustering, scales according to separation r as g(r) ∼ r −6 , for r exceeding a few multiples of the particle diameter (i.e. the "far-field" regime).
Because HI occur over scales on the order of the particle size [4], which is often on the order of microns, it is challenging to experimentally observe the effects of HI on particle clustering in turbulent flows. Indeed, until recently such experimental inquiry was prohibited due to spatiotemporal resolution and perspective overlap preventing the identification of particle pairs with very small separations [6,7]. Yavuz et al [8] presented the first experimental evidence of extreme inertial particle clustering as r approaches the collision radius. By extending the analysis of the inertia-free theory by Brunk et al [4] to the case of weakly inertial particles (St 1) using the driftdiffusion model of Chun et al [9], Yavuz et al claimed that this RDF enhancement was due to the combined effect of particle inertia and particle-pair HI [8]. Their resulting theory predicted that g(r) would scale as r −6 for a regime of r; however, they did not observe this scaling in their experiments. They speculated that the r −6 regime must occur at a scale smaller than they were able to reliably measure. Their theory did, however, predict a * andrew.bragg@duke.edu † The first two authors have equal contribution. ‡ huimeng@buffalo.edu contribution to the RDF arising from particle inertia that appears at slightly larger separations, and they claimed that this contribution explained the extreme clustering observed in their experiments. However, as we demonstrate, the theory by Yavuz et al [8] contains a number of errors. When these errors are corrected, the theory actually indicates that the inertial contribution to HI hinders clustering instead of enhancing it. Therefore, this theory by Yavuz et al [8] cannot explain the extreme clustering they observed. In addition, their measurement data exhibited significant scatter at the small-r region where the RDF increased explosively, which does not inspire confidence in the new data they offer.
Recently, Hammond and Meng [6] reported a new experimental RDF measurement of inertial particles (St = 0.74, and a = 14.25µm, where a is particle radius) in isotropic turbulence at r down to near-contact (r/a = 2.07). Using a novel particle tracking approach based on four-pulse Shake-the-Box (4P-STB), they obtained RDF and relative velocity results at unprecedentedly high resolution. When r went below r/η = O(1) (η: Kolmogorov length) corresponding to r/a = O(10), their measured g(r) grew explosively due to particle-particle interactions. The order of magnitude of g(r) at near-contact matched that of Yavuz et al [8]; however, their improved measurement resolution clearly showed that g(r) scaled as r −6 in this regime. This scaling is reminiscent of the prediction for the RDF of inertia-free particles subject to HI [4], hinting that HI may have driven the observed explosive clustering. In order to understand the mechanics behind these observed extreme clustering data, more theoretical analysis is required. Moreover, to validate the theory, more experimental data is desired.
In this letter, we present additional experimental RDF measurements over a range of St (0.07 to 1.06) and particle radius a (3.75 to 20.75µm), closely examine the theoretical analysis by Yavuz et al [8], correct errors contained therein, and compare the scaling exponents predicted by the corrected theory with those of the new experimental dataset. We found enormous discrepancies between the theoretical predictions and experimental results. Therefore, we also investigated potential explana-tions for these discrepancies.
Theory -We now derive a solution for g(r) while providing a physical explanation for how HI leads to particle clustering in turbulence. As shown in the Supplemental Material, the steady-state RDF g(r) may be expressed exactly as (with t → ∞) [10] where W is the relative velocity between two particles, ∂ s ξ ≡ W, and · r denotes an ensemble average conditioned on the particles having the separation ξ(t) = r. For fluid particles, ∇ · W = 0 in an incompressible flow and so g(r) = 1, i.e. they do not cluster. However, if ∇ · W is finite, clustering may occur with g(r) > 1. If we consider monodisperse particle-pairs with radius a that experience HI, then for St → 0 we have ∇ · W = λS [4], where λ ≥ 0 is a non-dimensional, nonlinear function of r/a that characterizes the HI, and S is the fluid strain-rate parallel to the particle-pair separation vector. Since λ ≥ 0, then the particle field will be compressed in regions where S < 0, and dilated in regions where S > 0. That ∇ · W = 0 is due to the disturbance fields in the flow produced by displacement of the fluid around the two particles, which in turn generates forces on the particles. This force either causes the particles to be attracted or repelled from each other, and vanishes for fluid particles (a = 0) since they do not disturb the flow.
Using ∇ · W = λS in (1) we see that g(r) > 1 is associated with a preference for trajectories with t 0 λ(ξ(s))S (s) ds < 0, that arises precisely because the particles are compressed into regions where S < 0. This phenomenon is similar to the case of inertial particles with St 1 (without HI) whose clustering is driven by preferential sampling of weak-vorticity, high-strain regions of the flow [9,11,12], that arises due to the particles being centrifuged out of vortical regions of the flow [13].
The HI effect on clustering is dependent on St. Since HI only occur when r is sufficiently small, we define a as the lengthscale of the hydrodynamic disturbance, below which HI become appreciable. At r > a , HI are not important, and the clustering arises solely due to how inertia modifies the particle interaction with the turbulence [11,12]. For r < a and St 1, the physical mechanism leading to RDF enhancement comes from particles being compressed into regions where S < 0 as discussed above, with sub-leading corrections to the trajectories due to inertia. For r < a and St ≥ O(1), the mechanism generating g(r ≤ a ) > 1 will be strongly affected by the non-local dependence of W(ξ(s), s) upon the turbulence the particles have experienced along their path-history at times s < s [11,12].
While (1) is useful for understanding how particles cluster, it is not straightforward to derive from this a closed expression for g(r). Yavuz et al [8] developed a theoretical model for g(r) in the regime St 1, based on the drift-diffusion models of [4,9]. The model assumes monodisperse particle-pairs that experience HI and Stokes drag forces, suspended in a turbulent flow. Unfortunately, we have found several errors in their analysis. In the Supplemental Material we derive in detail the correct version of the theory, confining attention to the far-field asymptotic behavior (r a, although in practice the far-field asymptotics are valid down to r/a ≈ 2.05 [4]), and retaining terms up to order St 2 , the same regime considered in [8]. This gives rise to where explicit forms of the coefficients µ 1 , µ 2 , µ 3 , µ 4 are given in the Supplemental Material. Two of the four terms in (2) have been characterized in the literature previously. In the absence of HI, µ 1 = µ 2 = µ 3 = 0, and g(r) ∼ (r/a) −St 2 µ4 describes the clustering due solely to particle inertia [9]. The leading HI contribution exp(µ 1 a 6 /r 6 ) is the far-field form of the result derived in Brunk et al [4], which is independent of St and describes the clustering due to HI that can occur even for St = 0.
The O(St 2 ) inertial contribution to the clustering arising from HI, exp(St 2 µ 3 a/r) was first derived in [8], where they determined µ 3 by fitting exp(St 2 µ 3 a/r) to their experimental data, obtaining µ 3 > 0. However, we show in the Supplemental Material that the theory specifies µ 3 ≤ 0, meaning this term inhibits clustering. Therefore, the scaling of the explosive RDF observation of Yavuz. et al. [8] cannot be justifiably associated with exp(St 2 µ 3 a/r).
Finally, the leading order contribution in St is exp(Stµ 2 a/r). In Yavuz et al [8] this contribution is absent since they argued that the third-order correlation on which µ 2 depends is zero for isotropic turbulence. As a result of this they concluded that the leading-order effects of particle inertia occurs at O(St 2 ). However, as we demonstrate in the Supplemental Material, this correlation cannot be zero. It is associated with the average amplification of the strain-rate and vorticity fields in turbulence, which are associated with the energy cascade [14] and are necessarily finite in three-dimensional turbulence [15]. Not only does this mean that for the corrected theory, inertia affects g(r) at O(St), instead of O(St 2 ), but also µ 2 ≤ 0, so that the leading order effect of inertia in the presence of HI is to suppress the clustering, not enhance it.
Experiments -The experimental dataset presented in this paper was acquired by following the new smallr particle tracking methodology developed in Hammond and Meng [6]. We measured g(r) at 10 different St and 4 different particle sizes, listed in Table I. The result at St = 0.74, a = 14.25µm is identical to those of [6].
The experiments were performed in the enclosed, fan-driven, Homogeneous Isotropic Turbulence chamber (HIT chamber). The complete turbulence characteristics of this chamber are detailed in [16]. The particles   [16].
were hollow spheres by 3M (3M Glass Bubbles, types K25, S60, and IM16K), which allowed particle size control through sieving and inertia control through choice of particle type [17]. We sieved the originally widely polydisperse particles to acquire narrower, albeit still polydisperse, size distributions for each of the four particle samples. Particle density was measured with a Micromeritics accu-Pyc II 1340 gas pycnometer. Uncertainties in r and g(r) were calculated following the method of [6] and are presented in the Supplemental Material. These two uncertainties had similar magnitudes across all conditions. Convergence of the RDF was achieved and the standard error was < 2%. For the a = 14.25µm particles, 15465 realizations were acquired for statistical convergence as in [6]. For the remaining three particle types, 9279 realizations were acquired.
Results -The experimental results for g(r) all 13 flow and particle combinations are shown in Fig. 1. At larger r, we observe the behavior g(r) ∼ r −St 2 µ4 from (2). The RDF in this regime is consistent with previous experiments [18]. When r/a decreases to r/a ≈ 30 in Fig.  1(a) for a = 3.75µm and r/a ≈ 12 in Fig. 1(b) for a = 14.25µm, g(r) grows explosively for all St, attaining values that are two orders of magnitude larger than those observed in previous simulations of inertial particles in turbulence without HI [19]. In this explosive regime, g(r) − 1 ∝ (r/a) −6 . This is consistent with the far-field form of (2) in the limit St → 0.
When r/a further decreases to below r/a ≈ 10 in Fig.  1(a) and r/a ≈ 3.5 in Fig. 1(b), g(r) flattens out. This is most likely due to particle polydispersity, which is known to cause g(r) to asymptote to a constant value at r < r c , where r c is a cut-off scale that increases with increasing polydispersity in the system [9,[20][21][22][23][24]. In our experiments, as particle size decreased, polydispersity increased (quantified by the ratio of the standard deviation to the mean of the particle size distribution), and correspondingly for smaller a the flattened region in g(r) broadened. 25µm. Also shown is the behavior g(r) ∝ r −6 . Plot (c) compares results with same/similar St but different a and Re λ (see Table I).
sive regime, g(r) increases with increasing St, with the scaling g(r) − 1 ∝ (r/a) −6 preserved. This would indicate that increasing St weakly but consistently enhances the RDF, within the uncertainty of the RDF measurement. However, even for the cases where St 1 such that the theory applies, this is fundamentally inconsistent with (2), according to which particle inertia does not affect the r −6 scaling, and for the HI contributions, reduces rather than enhances the clustering.
To investigate if the RDF collapses on r/a as predicted by the theory, in Fig. 1 (c) we plot RDF at three different St, each obtained for two different particle radii a. It can be seen clearly that for St = 0.36 and St = 0.37 (red curves) the RDF generally collapses over decreasing r, up until r/a ≈ 6. For the St = 0.74 (blue) case, the results might collapse down to r/a ≈ 4, as this is within the uncertainty of r measurement. For the St = 0.23 (green) case, the results clearly do not collapse, which is inconsistent with (2).
Using data from direct numerical simulations (DNS) to prescribe the fluid statistics on which µ 1 depends (see Supplemental Material), we find µ 1 ≈ 31.98. By contrast, for St = 0.07 and a = 3.75µm we find the proportionality coefficient for the relation g(r) − 1 ∝ (r/a) −6 to be 3 × 10 8 , and for St = 0.23 and a = 8.75µm it is 1.5 × 10 6 . This is an enormous discrepancy between the theory and experiments.
Discussion -Yavuz et al [8] claimed that their data does not allow them to observe g(r) ∼ exp(µ 1 a 6 /r 6 ), but that they do observe g(r) ∼ exp(St 2 µ 3 a/r), to which they fit their data to indirectly obtain µ 3 . As discussed earlier, this claim is highly problematic because while their fit yields µ 3 > 0, the theory requires µ 3 ≤ 0 (and using our DNS we estimate µ 3 ≈ −40.15). As such, their observations cannot be justifiably associated with g(r) ∼ exp(St 2 µ 3 a/r). Moreover, for some of the cases, their fit to g(r) ∼ exp(St 2 µ 3 a/r) is not that strong. Indeed, their case with a = 10µm, St = 0.19 is quite well described by g(r) − 1 ∝ (r/a) −6 over the range 6 < r/a < 11, the same scaling we observe. However, just as found in our data, their data implies a proportionality constant orders of magnitude larger than the theoretical µ 1 .
All of these findings show that the extreme clustering observed here and in [8] cannot be correctly described by a theory based on the HI of weakly inertial particlepairs, contrary to the claims of [8]. We have therefore sought to understand which assumptions in the theory may be responsible for its catastrophic failure to predict g(r) (noting that although the theory correctly predicts the scaling g(r)−1 ∝ (r/a) −6 , it probably does so for the wrong reasons given the enormous quantitative errors). To this end, we investigated four potential error sources.
First, the theory assumes HI between a particle pair; however, many-body HI could occur if three or more particles are found within each other's hydrodynamic disturbance field. To test if many-body HI occurred in our experiments, we calculated the average number of particles in a sphere of radius R around a test particle of radius a, conditioned on there being at least one satellite particle around the test particle, denoted by N (R|θ). For particle-pairs, N (R|θ) = 1, while N (R|θ) ≥ 2 indicates more than two particles in the sphere. Figure 2 shows the results for N (R|θ), where sphere size R/a can be likened to separation r/a. We found that in the range of r/a where g(r) grows explosively the 14.25µm and 8.75µm particles have N (R|θ) close to 1, while the 3.75µm and 20.75µm particles have N (R|θ) up to 3. This variation is due to different particle number densities. Although this could mean that many-body HI is playing a role for the 3.75µm and 20.75µm particles, many-body HI definitely do not for the 14.25µm and 8.75µm particles. Since extreme clustering is observed among all our experiments, many-body HI cannot be the fundamental cause of the discrepancy with the theory. Second, particles were assumed smooth spheres in the theory. To test this assumption, we sampled particles from the flow facility during operation. Microscope images show that 87% of particles were smooth spheres with no agglomeration (see Supplemental Material for images). Therefore this is not a cause of the discrepancy between theory and experiment.
Third, the neglect of other physically relevant forces on the particle-pair motion in the experiments, such as electrostatic and/or van der Waals forces etc. It is straightforward to show, however, that these forces would lead to behavior that is very different than g(r) − 1 ∝ (r/a) −6 (see, e.g. [25]), and therefore cannot be the explanation. Moreover, we have measured the charge level in the flow facility, determining it as < O(10 −16 )C (see Supplemental Material). Therefore, Coulomb forces on the particles were negligible.
Fourth, the particle Reynolds number Re p was assumed to be small in the theory [4,8,9] such that Stokes flow around the particles is assumed. If we use the expression Re p = a 2 /τ η ν [4] then for our experiments, Re p 1. Therefore, this assumption holds.
Conclusions -More experimental evidence of extreme clustering of inertial particles at small separations in a turbulent flow corroborates earlier observations [6,8] and allows for a clearer look into the scaling of g(r) and the influence of St and a. Our data confirms g(r) − 1 ∝ r −6 in the explosive scaling regime, contrary to Yavuz et al [8]. We demonstrate that the corrected theory based on weakly inertial particle-pair HI cannot explain the extreme clustering, since the theory predicts an inhibition rather than enhancement of g(r) by the inertial contribution to HI, while in experiments St weakly increases the extreme clustering. Moreover, the theoretical predictions for the RDF are in error by orders of magnitude. As such, the extreme clustering observed here and in [8] remains something of a mystery. The particle equation of motion invoked in the theory is clearly missing some vital effect, which future work must seek to uncover.

Supplemental Material for "Hydrodynamic Interactions and
Extreme Particle Clustering in Turbulence" 1 Theory

Expression for the RDF
Let r p (t) denote the separation of a particle-pair, w p (t) ≡ṙ p (t) their relative velocity, and r a timeindependent coordinate field. The exact solution for the probability density function ( where we have assumed the particles are initially uniformly distributed over the domain V , · r denotes an ensemble average conditioned on ξ(t) = r, where ξ is defined by ∂ s ξ ≡ W. The formal definition of the field W is where the operator · r p (0),w p (0) r,u denotes an ensemble average over all initial separations and initial relative velocities, conditioned on a given realization of the flow u, and conditioned on r p (t) = r. This field reduces to W(r, t) = w p (t|r) in the limit where there is a unique initial condition that generates a trajectory satisfying r p (t) = r, such as is the case for fluid particles.
For a statistically stationary, isotropic system, the dependence upon r reduces to a dependence upon r ≡ r , and the PDF ρ(r) is related to the radial distribution function (RDF) as where N is the total number of particles in the flow, and n ≡ N/V . In the thermodynamic limit, g(r) = V ρ(r) and we therefore finally obtain which is the result quoted in the paper.

Drift-diffusion model for the RDF
In [2], the authors developed a theoretical model for how weakly inertial particles cluster in the presence of HI and turbulence. The basis for their result was the drift-diffusion models of [3,4]. However in going through their analysis we have found several significant errors. Therefore, in the following, we re-derive the correct form of the theory and discuss the differences with the result obtained in [2]. Following [4], the steady-state form of ρ is determined by the transport equation (it is to be understood that t → ∞) where q D and q d are the diffusion and drift vectors, respectively. The transport equation above provides an integro-differential equation for the PDF ρ(r). In the limit where the correlation timescale of the flow is very small, the equation becomes differential, since over these short correlation times ξ(s) ≈ r, and this is referred to as the "diffusion approximation" [3,4]. As discussed in [4], the diffusion approximation is not valid in real turbulence since the correlation time of the flow is of the same order as that on which ξ(s) evolves. Nevertheless, in [4] it is demonstrated that in real turbulence, the diffusion approximation gives the correct functional forms in the model, and the quantitative error associated with the diffusion approximation can be corrected for using a "non-local" correction coefficient, for which a model was constructed in [4]. In view of this, we will also adopt the diffusion approximation (with the non-local correction coefficient to be explicitly included later) which is the form of the drift-diffusion model used in [2]. To construct solutions for ρ(r) in the regime St 1, we must specify the correlations in q D and q d , which are constructed using solutions to the particle equation of motion.
Just as in [2], we consider monodisperse pairs of small, heavy, inertial particles, subject to HI, whose equation of motion is where Υ p is the total angular velocity of the two particles (sum of their angular velocities), J p = J(x p (t), r p (t), t), f (x, r, t) ≡ Γ · r − 2a D rr and x refers to the arguments of the velocity gradient Γ(x, t) ≡ ∇u and the strain-rate S(x, t) ≡ (∇u + ∇u )/2. The terms A, B, C, D, E are nondimensional functions of r ≡ r only, and we will discuss them in more detail momentarily. Note that in the expression for f above, it has been assumed that r is sufficiently small so that the background flow (i.e. the flow field in the absence of HI) is locally linear. In a turbulent flow this formally requires assuming r η, where η is the Kolmogorov length scale. However, it should be appreciated that in practice this locally linear flow assumption is known to be valid up to r = O(10η) (e.g. see figure  12 of [5]). This is due to the fact that the Kolmogorov theory underestimates the scale at which viscous and inertial forces balance, which may be due in part to the phenomena of nonlinear depletion [6].
By simply inserting the expression for w p (t) given in (10) into the termẇ p that appears in (10), then to order O(St 2 ) we obtain Then, based on the definition of the field W(r, t) given in (2) we obtain whereΘ Note that owing to the definition of W(r, t), in (14) the terms J and f are explicitly J(x p (t), r, t) and f (x p (t), r, t), i.e. measured at fixed separation r, but along the time-dependent trajectory x p (t).
For clarity, we now switch to index notation, and write the result for W in the form of an expansion in St so that to order St 2 we have i (r, t) ∂ ∂r l W [1] l (r, s) ds .
The required far-field forms are then Here we have thrown away the terms involvingΘ, since as noted in [2], this term is sub-leading in the far-field under the assumptions St 1 and that the local flow around the particles is steady Stokes flow (already assumed in the particle equation of motion). The results above match exactly with those in [2] (although our notation differs), and retain the leading order contributions from the HI effects.
The first contribution in (19) is l (r, s) ds ∼ 75 where τ S is the correlation timescale of the strain-rate, and S 2 ≡ S p ab (t)S p ab (t) . This result is the same as the far-field version of the drift velocity derived in [3].
In [2], it is argued that the order St contributions to q d disappear due to isotropy of the flow. We now investigate this claim. The second contribution in (19) involves Invoking isotropy, incompressibility and the Betchov relations [9] we have and In [2], it is assumed that this quantity is zero, but we will now show that it is not. Just as in [2,4] we may write this as where Σ denotes the autocorrelation for the quantity. By introducing the strain-rate and vorticity, the singletime invariant can be written as This invariant is not zero; the first term on the rhs is the strain-rate production term, and the second is the enstrophy production term (i.e. vortex stretching) [10]. These invariants are always finite in three dimensional turbulence, and are directly connected to the energy cascade process [11]. Thus, remarkably, the processes governing the energy cascade also contribute to inertial particle clustering in the presence of HI. Note, however, that these terms are zero for two dimensional turbulence.
Putting this together we obtain where τ Σ is the timescale associated with the autocorrelation Σ, and Γ 3 ≡ Γ p bc (t)Γ p bd (t)Γ p dc (t) . Note that since Γ 3 < 0 in three-dimensional turbulence [10], then this contribution actually opposes the clustering of the particles.
The third contribution in (19) Similar to before where Σ denotes the autocorrelation for the quantity, and we have used the Betchov relations to obtain Putting this altogether we obtain where τ Σ is the timescale associated with the autocorrelation Σ . The fourth contribution in (19) t 0 W [1] i (r, t) ∂ ∂r l W [1] l (r, s) ds These integrals involve the quantity (with differing index labels) and this may be re-written using isotropy and an autocorrelation function as where Γ 4 1 ≡ A aabb (t, t), Γ 4 2 ≡ A abab (t, t), Γ 4 3 ≡ A abba (t, t). Using this we then obtain for the fourth contribution in (19) t 0 W [1] i (r p (t), t) ∂ ∂r l W [1] l (r p (s), s) ds We now gather together all the contributions to (19), and at each order in St only retain up to the leading order contribution from HI, yielding Concerning the diffusion term, in the far-field this is to leading order [2] where B nl is the non-local correction coefficient discussed earlier. Using (42) and (41) in (5), together with (3), we obtain the following solution for the RDF where Note that we have used ∼ in the solution for g(r) to reflect the fact that for St 1 there is an O(1) boundary condition at r = η [4,5].
This result in (43) differs from that of [2] in a number of crucial ways. First, in their analysis µ 2 is absent since they incorrectly assumed that Γ 3 = 0. The fact that µ 2 = 0 means that the leading order HI contribution associated with particle inertia occurs at O(St), not at O(St 2 ) as claimed in [2]. Moreover, since Γ 3 < 0, then µ 2 < 0, meaning that the leading order effect of particle inertia is to reduce the clustering associated with HI. Second, in the study of [2] they do not directly compute or estimate µ 4 , but rather infer its behavior indirectly by fitting their RDF expression to their experimental data. With this procedure they find that µ 3 > 0. However, the theory shows that µ 3 < 0, if we assume Γ 4 1 , Γ 4 2 , Γ 4 3 are all positive, as seems reasonable (see also below). As discussed in the paper, we suggest that this is because the behavior they observed in their experiment is not described by the theory, and hence the coefficient inferred from their data does not in fact correspond to µ 3 .
Finally, we note that we could have split up the correlations involving Γ p ij into separation contributions from the strain-rate and vorticity, as was done in [4]. This allows the theory to capture the influence of the different Lagrangian timescales associated with the strain-rate and vorticity fields [5]. However, this would have further complicated out result and would not change the conclusions drawn from our analysis.
The value µ 2 ≈ −0.02 is very small, and as a consequence the O(St) contribution in (43) is only the leading order inertial correction to the HI induced clustering if St < µ 2 /µ 3 ≈ 6.14 × 10 −4 , i.e. extremely small. Yavuz et al claimed that Γ 3 = 0 for an isotropic flow, and therefore µ 2 = 0. We argued that this is incorrect, and our DNS supports this showing Γ 3 = −0.15τ −3 η . Instead, the reason why µ 2 is small is because according to our DNS, τ Σ = 0.0034τ η . This occurs because the associated autocorrelation function passes through zero at time lag ≈ 2.4τ η and then exhibits a significant negative loop, leading to a small value for the correlation timescale τ Σ .
The DNS data also confirms our expectation that µ 3 is negative, demonstrating that the positive exponent measured in [2] by fitting the RDF expression to their experimental data cannot be associated with µ 3 .
The predictions from the theory (assuming boundary condition g(r = 100a) = 1) are shown in figure 1 for the far-field regime r 3a for which they are supposed to be valid. Figure 1 (a) clearly shows that while inertia enhances the RDF, it does not enhance it by orders of magnitude. Furthermore, figure 1 (b) shows the predictions when using µ 4 = 0 in order to consider the role of the contributions that arise solely due to HI. From this it is clear that the inertial modification of the HI effects act to reduce the clustering compared to the St = 0 case, rather than enhance it. Therefore, the enhancements due to inertia observed in figure 1 (a) are solely due to the term involving µ 4 which describes the clustering that arises even in the absence of HI, and is associated with the centrifuge effect and the associated preferential sampling of strain over rotation along the particle trajectory [4].

Error bars and uncertainty
The uncertainties have been discussed in detail in [13]. For the uncertainty in r, the primary source of error arises from the use of track interpolation to identify particle position and thus particle-pair separation. This uncertainty appears since the non-zero radial relative velocity may change the instantaneous value of r in the time between laser pulses. If fluctuations in r occur that are over shorter timescales than the time between frames ∆t 2 , these fluctuations will not be recorded in the track. As such, we take the product of the standard deviation of radial relative velocity and ∆t 2 to estimate the range of separations which may contribute to the recorded data at the given datapoint. For the uncertainty in RDF, we found that the random error (quantified by the standard error) was extremely small (< 2%), since the data was well-converged. The potential for variation in RDF instead arose instead in the selection of inputs for STB. As such, we varied the most important, consequential input parameter in STB, the maximum allowable triangulation error by ±10%, and took twice the standard deviation of the resulting RDFs as the vertical error bar in RDF.
In Fig. 2, we plot the interpolation uncertainty (in r) as a shaded region for the middle-most St (St = 0.16 in Fig. 2a and St = 0.74 in Fig 2b) in both plots of the experimental results. Similarly, we plot the RDF uncertainty by ensemble forecast as vertical error bars for the middle-most St. As a result, we find that the variations among the different St results across r are within the uncertainty limits, but there are clear, weak positive trends in RDF with St. Additionally, the comparison of the uncertainty in r of the St = 0.74 case in Fig. 2(b) with the results of the main letter shows that for across the four particle types, the results of the a = 8.75µm, 14.25µm, and 20.75µm particles generally collapse over decreasing r in the explosive r −6 regime.

Particle photographs to test for agglomeration
To test if some of the particles inside the flow facility had agglomerated, we sampled particles from the HIT chamber, and took images using a microscope. The particles were sampled by applying an adhesive to a glass microscope slide, and inserting it into the flow facility during fan operation while seeded with the a = 21µm particles. A crop of this microscope image is shown in Fig. 3.
Of 1684 uniformly sampled particles from the photograph, 87% were regular spheres, while 7% were double spheres (two spheres connected in the manner two soap bubbles touch, with varying individual radii) and 6% were cracked or broken. As stated in the Letter, due to the enormous difference between the theoretical prediction and the experimental results, it is not expected that the 13% of the particles that deviated from a regular spherical shape solely led to the observed difference between theory and experiments.

Electric charge effect on the RDF experiments
It was argued in the main text that electric charge does not play an important role in the particle clustering we observe. Here, we support this argument with experimental measurements of charge in our chargeminimized flow facility, where we find the Coulomb force is negligible even at the smallest observable separation r ≈ 2.07a, and thus cannot be responsible for the observed extreme clustering.
Particle charging can potentially occur in our facility by friction charging due to particle contact with the fans and inside walls of the flow chamber, a phenomenon known as the triboelectric effect. The magnitude of the charge generated by this process depends on the difference in the work functions of the two materials [14]. To prevent particle charging, we coated the flow facility fans and walls in carbon conductive shielding paint (Stewart Macdonald), with work function ∼ 5 eV [15,16], to match the work function of silicon dioxide, the primary compound of the particles [15]. The flow facility was grounded, such that the conductive coating would mitigate residual or pre-existing particle charge.
To test if this material combination produced minimal charge, we measured particle charge in a smallscale turbulent flow facility [17]. The experiment was performed on the a = 14.25µm glass bubble particles inside of a turbulent impinging flow tube with identical fans and surface coatings to the HIT chamber. We found that the resulting charge distribution was bipolar, with a mean of 3.5 × 10 −17 C (C: Coulomb) and standard deviation 4.0 × 10 −16 C.
To verify that the level of charge was similar in our full-scale HIT chamber, we also measured the electric charge on particles in situ. We used our charge measurement technique described in [18], modified to observe deflection of particle trajectories in a strong electric field using Shake-the-Box particle tracking. By comparing the particle terminal velocity before and after the introduction of an electric field, we estimated the electric drift velocity v x for charge measurements [18]. We found that even when initially detectably charged particles were injected into the flow facility (≈ 1 × 10 −14 C), after a time duration shorter than the startup time of the experiments (∼ 30s), particle charge was undetectable. We have calculated the resolution limit of this in-situ charge measurement system as 1.2 × 10 −15 C, based on the particle position resolution of 0.5 pixels. The electric charge is expected to be below this resolution limit.
To test if particle charge at the level of 10 −15 C could explain the extreme clustering observed, we compared the measured relative inward velocities of the particles with that which would be expected for oppositely-sign charged particles under idealized Coulomb attraction with charge q = 10 −15 C. In particular, we consider Coulomb attraction between two low-Reynolds-number particles (subject to Stokes drag) with charge magnitude q and opposite charge sign, for which the magnitude of their inward relative velocity would be w(r) = k e q 1 q 2 6πµa where k e is Coulomb's constant and µ is the dynamic viscosity of the fluid. In [13], we have reported inward relative velocity statistics as a function of r for the St = 0.74, a = 28.5µm particles used in this study. If q ≈ 10 −15 C, then (48) predicts relative velocities that are smaller by two orders of magnitude compared to the measurements. In the regime where g(r) − 1 ∝ (r/a) −6 , this gap widens to three orders of magnitude. Since the predicted relative velocity due to the Coulomb force between two particles with the largest possible electric charge is 3 orders of magnitude smaller than the average observed relative velocities, this offers additional evidence that the Coulomb force is weak at the observed spatial scales, and not an explanation for the extreme clustering observed in the experiments.