Detailed characterization of extreme clustering at near-contact scales in isotropic turbulence

Abstract Recent measurements of inertial particles in isotropic turbulence (Hammond & Meng, J. Fluid Mech., vol. 921, 2021, A16) revealed surprising extreme clustering of particles at near-contact separations $(r)$, whereby the radial distribution function, $g(r)$, grows from $O(10)$ to $O(10^3)$ with a $(r/a)^{-6}$ scaling (where $a$ is the particle radius), and a surprising upturn of the mean inward particle-pair relative velocity (MIRV). Hydrodynamic interactions (HIs) were proposed to explain the extreme clustering, but despite predicting the correct scaling $(r/a)^{-6}$, the HI theory underpredicted $g(r)$ by at least two orders of magnitude (Bragg et al., J. Fluid Mech., vol. 933, 2022, A31). To further understand the extreme clustering phenomenon and the relevance of HI, we characterize $g(r)$ and particle-pair kinematics for Stokes numbers $0.07 \leq St \leq 3.68$ in a homogeneous isotropic turbulence chamber using three-dimensional (3-D) particle tracking resolved to near–contact. A drift–diffusion equation governing $g(r)$ is presented to investigate the kinematic mechanisms of particle pairs. Measurements in all 24 conditions show that when $r/a\lessapprox 20$, extreme clustering consistently occurs, scaling as $g(r) \sim (r/a)^{-k}$ with $4.5 \leq k \leq 7.6$, which increases with $St$. Here $g(r)$ varies with $St$, particle size, density and polydispersity in ways that HI cannot explain. The extreme clustering region features an inward drift contributed by particle-pair turbophoresis and an inward radial relative acceleration. The latter indicates an interparticle attractive force at these separations that HI also cannot explain. The MIRV turns upward when approaching the extreme clustering region, opposite to direct numerical simulation predictions. These observations further support our previous assessment that extreme clustering arises from particle–particle interactions, but HI is not the main mechanism.


Introduction
Understanding how particles collide in isotropic turbulence is critical for modelling droplet coalescence in cloud turbulence (Shaw 2003;Grabowski & Wang 2013), particle coagulation in the chemical and process industries (McMillan et al. 2013) and mass entrainment within combustion engines (Wei et al. 2011).The rate of collisions in an isotropic flow can be described by the collision kernel (Sundaram & Collins 1997;Wang, Wexler & Zhou 2000), which depends on two statistics: the radial distribution function (RDF), g(r), and the mean of the inward particle-pair radial relative velocity (RV) w r − , abbreviated here as 'mean inward relative velocity (MIRV)', evaluated when the centre-to-centre particle separation r equals the contact value (e.g.2a for monodisperse particles).
It is well known that turbulence causes particle clustering through the interaction of particle inertia and turbulent eddies (Maxey 1987;Chun et al. 2005;Bec et al. 2007;Gustavsson & Mehlig 2011;Bragg & Collins 2014a;Bragg, Ireland & Collins 2015a,b;Gustavsson & Mehlig 2016).In the past, experimental measurements of the RDF (Salazar et al. 2008) and RV (Dou et al. 2018), resolving particle separations (r) down to the Kolmogorov scale (η), have confirmed findings from direct numerical simulations (DNS) based on one-way coupling of turbulence laden with small and heavy inertial particles: the RDF increases to values of O(10 0 -10 1 ) as r decreases down to O(η).However, recent advances in laser particle tracking techniques have enabled measurements at sub-Kolmogorov separations down to near-contact scales (Yavuz et al. 2018;Hammond & Meng 2021;Bragg et al. 2022) and revealed surprising extreme clustering of inertial particles never seen before.As r decreases below O(η), g(r) grows from O(1) to O(10 3 ), scaling as r −k where k ≈ 6 (Hammond & Meng 2021;Bragg et al. 2022).Hammond & Meng (2021) suggest that this extreme clustering, completely unexpected from inertial particle-turbulence interactions, must be driven by particle-particle interactions that become appreciable at small separations.
A frequently considered particle-particle interaction mechanism is hydrodynamic interaction (HI), i.e. the effect on a particle from the altered flow field due to the presence of another particle (Wang et al. 2005;Dhanasekaran, Roy & Koch 2021).Yavuz et al. (2018) attempted to explain their surprising extreme clustering data by extending the drift-diffusion model of Chun et al. (2005) to describe the clustering of weakly inertial particles with St 1 subject to HI, where the Stokes number St is defined as the particle response time τ p divided by the Kolmogorov time scale τ η .However, Bragg et al. (2022) found that the theory of Yavuz et al. (2018) contained several errors.Once these errors are corrected, the HI theory predicts that the RDF should behave as (Bragg et al. 2015b) g(r) ∼ (r/a) −St 2 μ 4 exp(μ 1 (r/a) −6 + (Stμ 2 + St 2 μ 3 )(r/a) −1 ), (1.1) where r is the particle-pair separation, a is the particle radius and μ i are coefficients that depend on the flow properties (the detailed definitions of μ i are discussed in § 5.2 and can be found in more detail in Bragg et al. (2022)).The leading HI term exp(μ 1 (r/a) −6 ) is the far-field form of the result derived in Brunk, Koch & Lion (1997) describing the clustering due to HI independent of St.In the limit as St → 0, one obtains the far-field relation g(r) − 1 ∼ (r/a) −6 as the leading contribution of HI to clustering.Both Hammond & Meng (2021) and Bragg et al. (2022) presented experimental measurements of g(r) that display an r −6 scaling in the extreme clustering regime, which suggest that HI could be the cause of the observed extreme clustering.However, the values of g(r) they measured were orders of magnitudes higher than the HI theory prediction for weakly inertial particles, even at values of r/a where the far-field HI assumption of the theory apply.Bragg et al. (2022) considered the assumptions made by the HI theory to determine if they could explain the discrepancy between the theory and the experiments.All assumptions were found to be valid for the experiments and could not explain the discrepancy.Thus, Bragg et al. (2022) speculated that the agreement of the r −6 scaling with experiments could be just a coincidence.They went on to suggest that 'the mechanism for the extreme clustering observed here, in Hammond & Meng (2021), and Yavuz et al. (2018) remains something of a mystery' and 'the particle equation of motion invoked in the theory is clearly missing some vital effect, which future work must seek to uncover' (Bragg et al. 2022).
To further understand the extreme clustering phenomenon and the relevance of HI, in this study we characterize extreme clustering in detail under a broad range of St conditions.
We aim to answer the following questions.
(i) Does extreme clustering consistently occur, at what separation, and to what magnitude?(ii) Do all RDF measurements follow the r −6 scaling predicted by HI theory?(iii) How does St affect the RDF, and is this consistent with the HI theory?(iv) What other flow/particle properties affect extreme clustering, and what does this reveal about the possible mechanism?
To that end, this study expands on our recent RDF (Hammond & Meng 2021;Bragg et al. 2022) and RV (Hammond & Meng 2021) measurements in an enclosed homogeneous isotropic turbulent airflow chamber using a high-resolution 3-D particle tracking technique designed for small-separation measurements.We varied the particle radius, particle shell thickness (which will vary the particle density) and chamber fan speed to achieve a broad range of particle inertia values with 0.07 ≤ St ≤ 3.68.For each case, we measured g(r), particle RV and relative acceleration (RA) statistics.To gain further insights into the potential driving forces of extreme clustering, we evaluate a purely kinematic relation governing g(r) derived exactly from an underlying probability density function (p.d.f.) master equation.This allows us to analyse the experimental data in terms of different contributions to the drift and diffusion mechanisms related to clustering (Bragg & Collins 2014a), whose behaviour could provide insights into at least the qualitative behaviour of the forces and mechanisms underlying the extreme clustering.We investigated the effects of St on the measured statistical quantities and compared the results with the HI theory.

Kinematic theory
Since the force(s) responsible for generating the extreme particle clustering at near-contact scales observed experimentally are yet to be understood, we cannot write a dynamical equation that governs the particle-pair relative motion as the particles approach each other.We can, however, derive a purely kinematic relation that governs g(r) to gain useful insight into the processes governing its behaviour and perhaps the elusive mechanism(s) responsible for the extreme particle clustering at near-contact scales observed experimentally in Yavuz et al. (2018), Hammond & Meng (2021) and Bragg et al. (2022).
Let r(t) and w(t) denote the particle-pair separation (based on the particle centres) and RV vectors, and w r (t) ≡ r(t) −1 r(t) • w(t) the radial (or longitudinal) component of the RV.Kinematically, we have ṙ(t) ≡ w r (t), where r(t) ≡ r(t) .By taking moments of the equation governing the joint p.d.f. of r(t) and w r (t), a set of equations governing the p.d.f.(r, t) ≡ δ(r(t) − r) may be derived (Bragg et al. 2015b where D/Dt ≡ ∂ t + w r (t) r ∇ r , r is a time-independent coordinate, • r denotes an ensemble average conditioned on r(t) = r and S 2 ≡ [w r (t) − w r (t) r ] 2 r is the second-order particle RV structure function.
The time-dependent RDF g(r, t) is related to via a constant β as = βg (Chun et al. 2005;Bragg et al. 2015b), and substituting this into (2.1) and rearranging gives an exact drift-diffusion equation governing the RDF: For a statistically stationary flow, (1.1) yields where Φ is the constant particle-pair mass flux (PPMF).In a system where there is no particle agglomeration (such as our experimental system; see Bragg et al. (2022)) and no injection of particle pairs at the large scales of the flow, a zero-flux state is expected with Φ = 0.In our experiments g > 0 (regions where g = 0 do not occur), which then implies w r (t) r = 0 from (2.4), and so (2.3) reduces to 0 = −S 2 ∇ r g + g( ẇr (t) r − ∇ r S 2 ). (2.5) Since in a turbulent flow S 2 > 0, the sign of ∇ r g is governed entirely by the sign of the drift term ẇr (t) r − ∇ r S 2 , which is not sign definite.In regimes where g(r) grows as r reduces, i.e. ∇ r g < 0, we must have which simply means that for the RDF to increase as r decreases, the total drift acting on the particle pair must be negative, i.e. inward.As such, it is an inward drift that produces particle clustering in the flow.The first term, ẇr (t) r , is the mean RA between the particle pair at separation r.We abbreviate this term to RA.The behaviour of this term will depend directly on the forces that act upon the particles.When RA is negative, ẇr (t) r < 0, the forces on particle pairs on average act in the negative r direction (meaning inward), thus promoting clustering.
The second term, −∇ r S 2 , is the negative gradient of the second-order structure function of the particle-pair RV with respect to r.We designate this as the particle-pair turbophoresis (PT).When the PT is negative, −∇ r S 2 < 0, it contributes to an inward drift of particles and thus promotes clustering.This condition corresponds to ∇ r S 2 > 0, or a positive gradient of S 2 , under which S 2 decreases with decreasing r.The particle-pair inward drift by this mechanism means that particles move in the direction of decreasing S 2 , or decreasing fluctuations of their radial RV.This turbophoretic drift phenomenon is analogous to the turbophoretic drift mechanism that plays a crucial role in governing particle concentrations in turbulent boundary layers (Reeks 1983;Johnson et al. 2020;Bragg, Richter & Wang 2021), where particles drift in the direction of the negative gradients of the wall-normal particle velocity variance (e.g.towards the wall).But in the case of PT, the particle-pair inward drift arises from the gradient of inertial particle-pair radial RV fluctuations.In summary, (2.6) shows that there are two distinct kinematic quantities associated with clustering, which can be examined in experimental measurement data.Exploration of the behaviour of ( ẇr )(t) r and −∇ r S 2 could help us gain insight into the physical processes generating the extreme clustering phenomenon.

Flow and particles
The 3-D particle tracking experiments were performed in an enclosed homogeneous isotropic turbulence (HIT) chamber using a LaVision particle tracking velocimetry system and LaVision's Four-Pulse Shake-The-Box (STB) method (Sellappan, Alvi & Cattafesta 2020) described by Hammond & Meng (2021).The flow chamber (pictured in figure 1) was a 1 m diameter truncated icosahedron ('soccer ball'-shaped) turbulent airflow facility.It was equipped with 20 symmetrically placed fans blowing to the centre, producing HIT in a 48 mm region at the centre of the chamber (Dou et al. 2016).The chamber was electrically grounded and the net electric charge on particles measured in situ was found to be negligible (Bragg et al. 2022).The turbulence strength was varied by changing the fan rotation speed.The Froude number Fr, defined as the ratio of the Kolmogorov scale for fluid accelerations to the acceleration due to gravity, of all flow conditions was in the range 4.4 < Fr < 24.2.For each case, we calculated the normalized settling velocity, Sv = St/Fr, to determine the possible impact of gravity on the particle dynamics.The largest value was Sv = 0.28, with the majority having Sv < 0.1.In all cases, Sv is far too small for gravity to be having a leading-order impact on the particle dynamics.A full characterization of the base turbulence in the HIT chamber was described in Dou et al. (2016).Optical glass was mounted on five facets of the 'soccer ball', allowing laser illumination, imaging and optical diagnostics.Particles to be studied were pneumatically injected from the side.
Table 1 shows all the flow conditions and particle properties for the 25 sets of experiments.We studied five different sets, for which one particle property condition (rows under 'Particle properties' in table 1) was examined in the HIT chamber under five independent flow conditions.This led to a total of 25 different experimental conditions in the range 246 < Re λ < 357 and 0.07 < St < 3.68.Among these 25 conditions (with distinct St values given), the a = 14.25 μm and St = 0.74 condition (dark grey) was previously presented in Hammond & Meng (2021), and an additional 12 cases (light grey) were studied in Bragg et al. (2022).The current study adds 12 new cases (white) to complete the matrix of conditions.Unfortunately, one set of data (a = 8.75 μm, Re λ = 324, St = 0.49) was accidentally deleted and unable to be reproduced, since the LaVision system necessary for reproducing the experiments was a university-shared instrumentation and had been returned.
The particles used in the experiments were three types of hollow glass spheres (3M glass bubbles, series K25, S60 and iM16K), each type with a fixed shell thickness, and thus the density varied as particle radius varied for a given type.The particle types and radius conditions were chosen to specifically achieve the variation and range of St examined.At purchase, the particles had a broad diameter distribution from less than 5 μm to 105 μm.To produce near-monodisperse particle samples, we sieved the particles into narrower size ranges using a Gilson GA-8 sonic sieve with standard test sieves following the procedure described by Dou et al. (2018).From the sieved particles we utilized five particle samples in four particle diameter ranges (5-10 μm, 15-20 μm, 25-32 μm and 38-45 μm).The average density of each particle sample was measured using a Micromeritics Accu-Pyc II 1340 gas-displacement pycnometer, documented in table 1.We used sieve mesh size increments of either 5 μm or 7 μm with variance in radius being either ±1.25 μm or ±1.75 μm, respectively.As a decreases, this variance became increasingly significant compared with the mean particle radius, thereby increasing particle polydispersity.The polydispersity (φ = 2σ a /a) is defined here as the range of particle radii 2σ a in a final sample divided by the average particle radius a in that sample.

Optical diagnostics
The Four-Pulse STB particle tracking velocimetry set-up and method were described by Hammond & Meng (2021).As shown in figure 1, the set-up consists of two Photonics Nd-YLF lasers and four Phantom Veo 640L cameras arranged in a cross configuration, 20 • from the centre plane in each cardinal direction.The laser beams travel through a series of optics, resulting in a 10 mm × 30 mm laser sheet entering the chamber.The actual measurement volume utilized is 10 mm × 30 mm × 50 mm in the centre of the chamber.In addition, the temperature inside the flow chamber was monitored during operation and kept within ±1 • C.

Particle tracking and processing
After the injected particles were allowed to equilibrate over approximately 100 large-eddy turnover times (≈30 s), particle tracking was performed on three independent sets of 3093 realizations of the particle-laden flow for the a = 3.75 μm, a = 8.75 μm and a = 20.75 μm particles, and five independent sets of 3093 realizations for the a = 14.25 μm particles.Using LaVision DaVis software, we defined one particle track as four consecutive images of the same particle from four laser pulses, temporally separated by t 1 , t 2 and t 3 as shown in figure 1.The setting of these time intervals varied with Re λ as shown in table 1.The tracks were then processed using the STB function in DaVis.All t values and STB input parameters are listed in table 1.The particle volume fraction was kept at approximately 2.2 × 10 −5 (equivalent to 0.002 particles per pixel), which makes our system a dilute suspension.(Please note that figure 6 in Bragg et al. (2022) is a photograph of particles sampled from the chamber using a microscope slide with two-sided tape, for the purpose of analysing the shape and agglomeration (possible fusion) of particles; this image is not an experimental image of particles dispersed in the turbulent flow.)For all particles detected in each realization of the turbulence, the relative position, RV and RA for every possible particle-pair combination were calculated.These kinematic terms were interpolated at the midpoint of the four-pulse track, between pulses two and three, using the position and time data from all four pulse points.These values were then used to calculate the two critical particle-pair statistics in the collision kernel, g(r) and w r (t) − r , as well as two components of the drift mechanism, ẇr (t) r and −∇ r S 2 .We accounted for the physical boundary of the experiment in calculating the RDF as described in Hammond & Meng (2021).To calculate the gradient of S 2 , we used a second-order-accurate discrete gradient operation.Since discrete gradients amplify random scatter in the data, we first used a first-order Savitzky-Golay low-pass filter with a window length of 3 on S 2 before calculating the gradient.We provide a detailed analysis of the measurement uncertainty in Appendix A.

Measurement uncertainties
Measurement uncertainty was calculated in much the same way as Hammond & Meng (2021).It is described in detail in Appendix A, which includes a comprehensive table of uncertainty values and supplemental figures with additional error bars.We have quantified statistical convergence through confirmation of minimal standard error across the measurements of the RDF.Of the possible uncertainties inherent in our measurement equipment, the interpolation uncertainty contributes most significantly to our position measurement uncertainty.The interpolation uncertainty arises from the finite amount of time between pulses two and three during our four-pulse tracking method.Though t 2 is quite small, it is still finite and causes uncertainty when calculating particle position between these two pulses.The interpolation uncertainty is negligible at large separations but becomes more significant near contact.Even accounting for this uncertainty, the trends in our data are still visible, as they extend outside of the uncertainty ranges.This is made clear in the figures in Appendix A, which include visualizations of the position uncertainty for select cases as horizontal error bars.Particle position uncertainty is significantly smaller than the interpolation uncertainty but is reported here for completeness.In addition, the overall value of the reported statistics could be influenced by the input variables for the STB algorithm.Though there are many user-defined input variables, it has been found that the allowable triangulation error, , is the most important (Novara et al. 2019).By varying by ±10 % for our 50 Hz cases, we produce the vertical error bars present in the figures throughout the paper.

The RDF results
We present the comprehensive RDF results on a log-log scale in figure 2, with g(r) as a function of normalized separations r/η in figure 2(a,c,e,g,i) and of r/a in figure 2(b,d, f,h,j).Figure 2(a,b), 2(c,d), 2(e, f ), 2(g,h) and 2(i,j) correspond to the five rows of experimental conditions in table 1, with St generally decreasing from top to bottom.Within each row, i.e. for each particle type, curves of different colours represent g(r) measurements at increasing fan speeds and thus increasing St.For all tests, the smallest separation we were able to measure accurately was r/a = 2.07 where a is the particle radius, which translates to different r/η values for different types of particles and fan speeds.Note that r/a = 2.00 is contact for monodisperse particles.
All 24 sets of g(r) data consistently show extreme clustering, with g(r) reaching O( 102 )-O(10 4 ).They demonstrate similar trends as r decreases from infinity, exhibiting first a moderate growth and then a rapid transition to explosive growth (extreme clustering) before finally levelling off.We denote the scaling exponent of the RDF by k, such that g(r) − 1 ∼ (r/a) −k .Evidently, k varies significantly in these different regions.For convenience of discussion, we identify four distinct regions: (I) inertial clustering region; (II) transition region; (III) extreme clustering region; (IV) decorrelation region.These regions are marked in figure 2 (b,d, f,h,j), where the insets provide a zoom-in of the transition region (II).Although our focus is on understanding the extreme clustering (region III), for completeness we will discuss the results from r/a = ∞ to r/a = 2.07 near contact, beginning in region I, the inertial clustering region.

Inertial clustering region (I)
From r/η = 100 down to roughly 1 or 2, g(r) gradually increases but remains at O(1), exhibiting the familiar clustering due to inertial particle interaction with turbulence, which 982 A21-8  has been extensively studied in theory (Chun et al. 2005;Zaichik & Alipchenkov 2007;Bragg & Collins 2014a), DNS (Saw et al. 2012a;Rosa et al. 2013;Ireland, Bragg & Collins 2016;Dhanasekaran & Koch 2022) and experiments (Salazar et al. 2008).Therefore, we refer to this region as the inertial clustering region (I).In this region, the scaling exponent k is small, ranging from approximately 0.07 to 0.22, and the data for g(r) is consistent with the power-law scaling (r/η) −c 1 that is well-documented in the literature, where the scaling exponent c 1 is specific to the inertial clustering region (I) (Reade & Collins 2000;Chun et al. 2005;Saw et al. 2012a,b).Note that we use k to designate the scaling exponent for every region of g(r), while c 1 is designated only for the inertial clustering region.The end of this region is in the range 0.8 < r/η < 2.7 depending on particle type and Re λ conditions.
To examine the effect of St on the inertial clustering, in figure 3 we plot the average scaling exponent k from the RDF measured in region I (denoted by the coloured shapes) as a function of St.To compare our data with the literature, we also plot the c 1 values averaged over Re λ from DNS of non-interacting small and heavy inertial particles in isotropic turbulence (Ireland et al. 2016), denoted by the black plus signs.Each experimental particle type (unique combination of radius, a, density, ρ, and polydispersity, φ) is presented as a group with a specific colour-symbol designation.The DNS predicts that as St increases from 0, c 1 increases from 0 quickly, reaching a maximum value of c 1 ≈ 0.  2021) is represented here in figure 3 by the filled fuchsia square.In their paper they reported the scaling exponent in the inertial clustering region (I) for this case to be k = 0.39, but this was a result of averaging the inertial clustering region (I) and the transition region (II), since they did not separate out the latter.Based on our current region definition, we recalculated the scaling exponent for their case to be k = 0.14 in the inertial clustering region (I) as shown in figure 3. Hammond & Meng (2021) also commented that their measured k value was much lower than c 1 = 0.69 predicted by DNS for monodisperse particles at a similar St value (St = 0.7) and attributed the discrepancy to the polydispersity of their particle sample.The presence of polydispersity has been known to diminish the value of c 1 for the overall particle sample based on the least clustered particle population (Saw et al. 2012a,b).However, in figure 3, we see no clear trend for k as the polydispersity, φ, increases for our samples.In fact, the particles that matched DNS (the red circles) had the highest polydispersity.No conclusions can be drawn, however, since in our experiments polydispersity was not an independent parameter but arose as a side effect of sieving particles of different sizes (smaller particles had larger relative size range and thus polydispersity).
The results clearly show that the scaling exponent k (or c 1 ) under matching St conditions has different values based on different particle properties.For example, the red circle case (a = 3.75 μm) and blue diamond case (a = 8.75 μm) under the same St (0.23) show an order of magnitude difference in k (0.22 and 0.07, respectively), as do the yellow star and green triangle cases (both a = 20.75 μm, with ρ = 0.5 g cm −3 and 0.3 g cm −3 , respectively) at similar St conditions (St = 1.96, k = 0.15 and St = 1.91, k = 0.08, respectively).The k value is more similar for the fuchsia square (a = 14.25 μm) and green triangle (a = 20.75 μm) cases at the same St (St = 0.74, with k = 0.14 and k = 0.09, respectively), and even more so for the blue diamond (a = 8.75 μm) and fuchsia square (a = 14.25 μm) cases at similar St conditions (St = 0.37, k = 0.12 and St = 0.36, k = 0.14, respectively).This inconsistency in k at similar St values is indicative that besides St, other parameters related to the particles are also influencing the value of k.
By examining the trend of k versus St in conjunction with particle properties, we see some interesting trends.The scaling exponent, k, increases as St increases for the smaller particles (a = 3.75 μm and 8.75 μm), but decreases as St increases for the larger particles (a = 14.25 μm and 20.75 μm).Besides particle size, some other particle properties could also be playing a role in the inertial clustering.The two a = 20.75 μm cases (green triangle and yellow star) present drastically different k magnitudes, where the denser particle type (ρ = 0.5 g cm −3 ) exhibited larger k values under the same St.However, it is premature to pinpoint particle density as the cause for the different k values; the particles in the green and yellow cases are two types of hollow glass spheres (K25 and S60, respectively) of the same radius endowed with different shell thicknesses.Besides different densities, they could come with different surface properties which we are unaware of.We therefore cautiously conjecture, based on our experimental observations, that besides particle inertia (represented by St), particle radius and particle density or a related particle property may also play a role in particle clustering, even at large separations in the inertial clustering region (I).This is consistent with the DNS results in Daitche (2015), where particles governed by the full Maxey & Riley (1983) equation were simulated in isotropic turbulence.Daitche (2015) showed that the clustering of the particles can depend strongly on the particle-to-fluid density ratio ≡ ρ p /ρ f as well as on a/η (within the a/η 1 regime).However, for the values of and a/η in our experiments, the DNS results of Daitche (2015) would indicate that the effects of and a/η are weak (e.g.see figure 8 in Daitche (2015)).Hence the dependence we see of k on and a/η, in addition to St, cannot be explained within the context of the additional forces in the Maxey & Riley (1983) equation that are important when is not sufficiently large or a/η sufficiently small.

Transition region (II)
We define the transition region (II) as that where the scaling exponent k increases to O(10 0 ).Details of g(r) in this region are shown in the insets in figure 2. In this region, the scaling exponent k starts to rise rapidly as the separation r decreases.The boundaries of this region are marked by the dotted lines in figure 2(b,d, f,h,j).As the particle size decreases, the onset of the transition region (II) occurs at increasingly large normalized separation distances r/a and the width of the region increases.
In the transition region, g(r) does not grow as a power law, and therefore in this region k is a function of r.To further investigate how the transition region (II) is affected by St and the particle properties, we plot the average gradient of k (equivalent to the second-order gradient of g(r)) in the transition region (II) (denoted by ∇ r k ) versus St in figure 4. , φ = 0.66 a = 8.75 µm, ρ p = 0.74 g cm -3 , φ = 0.28 a = 14.25 µm, ρ p = 0.31 g cm -3 , φ = 0.24 a = 20.75 µm, ρ p = 0.3 g cm -3 (K25), φ = 0.17 a = 20.75 µm, ρ p = 0.5 g cm -3 (S60), φ = 0.17 Figure 4 shows a weak increase of ∇ r k as St increases for the three sets of particles at St < 1 (red circles, blue diamonds, fuchsia squares) and a relative independence of St for the two sets of particles with St > 1, the green triangle (a = 20.75 μm, ρ = 0.3 g cm −3 ) and yellow star (a = 20.75 μm, ρ = 0.5 g cm −3 ) cases.Within each particle sample, an increase in St (through an increase in fan speed) corresponds to a slightly swifter transition.However, different particle samples with similar St do not have similar ∇ r k .This indicates that particle properties may also be contributing to the swiftness of transition between the inertial clustering (I) and extreme clustering (III) regions.

Extreme clustering region (III)
The extreme clustering region (III), as shaded in figure 2(b,d, f,h,j), is characterized by a drastic increase of g(r) with a constant scaling exponent k.Its onset occurs around r/η ≈ 1; thus, the extreme clustering is predominantly in the sub-Kolmogorov region.The fact that this dramatic change in the behaviour of g(r) occurs at sub-Kolmogorov scales suggests that it is likely driven by particle-particle interactions and not by particle-turbulence interactions (Hammond & Meng 2021).For curves in each row (corresponding to different fan speeds and thus St values for the same particle sample), the onset of the extreme clustering region (III) occurs at different r/η but collapses to the same r/a, which ranges from 7 to 25.From top to bottom in figure 2(b,d, f,h,j), i.e. as the particle radius a decreases and simultaneously the polydispersity φ increases, the onset occurs sooner and the region gets generally narrower, giving way to an increasingly earlier onset of the levelling off, i.e. the onset of the decorrelation region (IV).
In figure 2(b,d, f,h,j) we plot a dashed line to represent the average slope of g(r) in the extreme clustering region (III) for each particle type on the log-log plots.It is seen that as particles decrease in size, the scaling exponent k decreases.In figure 5   behaviour III to be where g(r) scales specifically as r −6 .This is accurate for their single case (represented as the filled fuchsia square in figure 5) but cannot be generalized for all the St and particle radius conditions in our expanded experiments.Also note that Bragg et al. (2022) did not differentiate between the k values and reported g(r) ∼ r −6 for all their cases.
The general trend that k increases with St is especially obvious within a given particle type (data of the same symbols in figure 5).However, the degree to which k increases as St increases is not consistent among different particle types and has no appreciable trend in terms of change in particle radius.Comparing the red circle (a = 3.75 μm) and the blue diamond (a = 8.75 μm) cases, as the radius increases, k decreases in general.But comparing the blue diamond and the fuchsia square (a = 14.25 μm) cases, k increases in general, contradictory to the first comparison.Further comparing the fuchsia square case with both the green triangle and the yellow star cases (both a = 20.75 μm), k decreases for the green triangle case but increases for the yellow star case.It appears that an increase in density for the same particle type, as seen in the green triangle and yellow star cases (ρ = 0.3 g cm −3 and 0.5 g cm −3 , respectively), causes a general increase in k in this region.Other particle properties seem to be playing a role in influencing the value of k.At the extremes, the red circle case (a = 3.75 μm) and blue diamond case It should be kept in mind that in our experiment, St was varied with a combination of particle properties and the fan speed in the turbulence chamber.For the same particle sample, the increase of St was achieved by increasing the fan speed alone.Moreover, just as we have seen in the inertial clustering region (I) and the transition region (II), here in the extreme clustering region (III), particle properties, separate from St, are influencing k and thus the degree of extreme clustering.

Decorrelation region (IV)
When the particle-pair separation further decreases past the extreme clustering region (III), g(r) begins to flatten, as expected from polydispersity in the samples.This is the decorrelation region (IV), so named because polydispersity decorrelates particle-pair relative motions and causes the RDF to flatten (Chun et al. 2005;Saw et al. 2012b;Dhariwal & Bragg 2018).In brief, samples that have a relatively wide variation (from polydispersity, φ) from the average St will have particles with differing responses to the fluid accelerations, independent of particle-pair separation.At a critical separation, the relative velocities become independent of separation and the drift mechanism driving the clustering vanishes, leading to the onset of the plateau of the RDF.A larger polydispersity will result in a larger critical separation, which is clearly observed in our data.Further, the critical separation is constant for a given particle type regardless of the fan speed condition (changing average St but not the variation from the average St).These two conditions reinforce that the behaviour of the decorrelation region (IV) is a predictable function of the polydispersity, φ, of the sample and unrelated to any particle-turbulence interactions.

Particle-pair kinematics results
As previously discussed in § 2, (2.6) shows that there two distinct kinematic effects that are exclusively associated with particle clustering and the growth of the RDF, ẇr (t) r and −∇ r S 2 .For g(r) to increase as the separation decreases, the total drift term in (2.6) must be negative, i.e. ẇr (t) r − ∇ r S 2 < 0. In this section we will examine the total drift and the two contributions to it, namely the particle-pair radial RA ( ẇr (t) r ) and PT (−∇ r S 2 ), along with the behaviour of g(r) predicted by (2.6).A deeper understanding of these kinematic quantities can lead to insights into the forces and mechanisms that dominate in the extreme clustering region (III) and near-contact region, e.g. by suggesting how the underlying force responsible for the extreme clustering must behave as a function of r.

Total drift
In figure 6, we present the total drift ẇr (t) r − ∇ r S 2 normalized by the Kolmogorov acceleration, a η (figure 6b,d, f,h,j), for each particle type and the corresponding RDF (figure 6a,c,e,g,i), both in terms of r/a.As before, the five rows of plots correspond to the five rows of particle conditions in table 1, and curves of different colours represent measurements at varying fan speeds and thus varying St.The boundaries of the regions are marked with dotted lines, and the extreme clustering region (III) is denoted by the shaded region.The total drift plots have a horizontal red dashed line at drift equal to zero to visually separate when the term is positive or negative.
Figure 6 shows that for all particle types the total drift is extremely negative in the extreme clustering region (III), attaining values that are up to two orders of magnitude larger than a η .Moreover, although difficult to discern, the total drift is also slightly negative in the inertial clustering region (I).Such negative drift is consistent with particle clustering based on the kinematic theory discussed in § 2. As polydispersity increases and the decorrelation region (IV) widens, the total drift term begins to fluctuate wildly, oscillating between positive and negative values.This wild fluctuation is likely due to the statistics not being fully converged and a relic of calculating the turbophoresis term (as discussed in § 4.2.3), which requires calculating the gradient of the second-order structure function.
It is peculiar that the total drift becomes slightly positive in the transition region (II) for all cases.A positive drift term should be inhibiting clustering according to (2.6), but the RDF is still increasing as separations decrease in this region.This anomaly could be due to one or both of the assumptions made in our kinematic theory being violated at these separations, namely, statistical stationarity (leading to (2.4)) and zero PPMF (leading to (2.5)).To understand this anomaly, we plot the normalized PPMF in figure 7, defined as the RDF multiplied by the particle-pair radial RV normalized by the Kolmogorov velocity: These plots include error bars determined via propagation using the uncertainty for the RDF and radial RV (Moffat 1988).For all cases, the PPMF is small through the inertial clustering region (I) and transition region (II).From this, we can say that the zero-PPMF assumption holds and thus cannot explain the anomaly of positive drift in the transition region (II).
Interestingly, the PPMF deviates from zero at smaller separations.For all cases except the a = 20.75 μm and ρ = 0.5 g cm −3 case, the PPMF deviation is positive and occurs in the decorrelation region (IV).Within the extreme clustering region (III), however, the deviations from zero are within measurement uncertainties.The a = 20.75 μm and ρ = 0.5 g cm −3 case is unique in that the PPMF deviation is positive, at least an order of magnitude larger than in the other cases, and occurs in the extreme clustering region (III).
According to (2.1), any deviation of PPMF from a constant value indicates a lack of statistical stationarity of the particle statistics.The statistical stationarity of our experimental set-up was investigated and documented in Dou et al. (2016); however, that study did not evaluate the PPMF.It is possible that this quantity is more sensitive to non-stationary effects than the quantities considered in Dou et al. (2016) and that the particles in our experiments may in fact not have fully attained a stationary state.This would be surprising given the time scale of the experiments, but could be caused by particle-particle interactions.Further investigation is required to understand this.4.2.2.The RA Figure 8(a,c,e,g,i) presents the mean particle-pair radial RA, ẇr (t) r , normalized by the Kolmogorov acceleration, a η , as a function of r/a.Recall that a negative RA indicates that particle pairs are accelerating towards each other, which means they are experiencing a net attractive force.In the inertial clustering region (I), the RA remains constant and very slightly negative at O(−0.1), which is consistent with the relatively weak clustering indicated by the RDF in this region.Just before the onset of the transition region (II), the RA remains negative and the magnitude increases significantly, indicating the onset of a significant attractive force that drives particles towards each other, which would contribute to the extreme clustering region (III).In region III, the RA remains negative, with its magnitude increasing to a maximum and then decreasing.The decrease in magnitude continues into the decorrelation region (IV) as particle pairs approach contact.Samples with a higher polydispersity, φ, experience lower RA magnitudes and an earlier onset of a plateau in RA near contact, mirroring the plateau observed in the RDF.4.2.3.The PT Figure 8(b,d, f,h,j) presents the negative gradient of the second-order structure function of particle-pair RV, −∇ r S 2 , normalized by the Kolmogorov acceleration, a η , as a function of r/a.We have termed this component PT.It describes the tendency of particles to move in the direction of decreasing spatial fluctuations of particle-pair RV.Recall that when PT is negative, it contributes to clustering.In the inertial clustering region (I), PT remains very slightly negative at O(−0.1).In the transition region (II), as the separation decreases, PT turns positive, reaching a peak magnitude just before the extreme clustering region (III), then decreases and becomes negative in the extreme clustering region (III).The PT remains negative and reaches a maximum magnitude much larger than RA in region III, indicating that PT is dominating the total drift term in region III.As discussed in § 4.2.1, the PT behaviour in region II cannot easily be explained when compared with the behaviour of g(r), as g(r) is increasing as separations decrease but the total drift (dominated by PT in this region) is positive.However, as separations decrease in region III, PT becomes negative, realigning with the behaviour of g(r), which continues to increase as separations decrease.In the decorrelation region (IV), PT is mostly positive but exhibits some extreme fluctuations, which drastically intensify as the particle radius decreases.Such fluctuations are likely amplifications of measurement uncertainties and postprocessing calculations, as briefly discussed in § 4.2.1.

The MIRV results
In order to evaluate the collision kernel, from our particle tracking data we extract the particle-pair MIRV, which is defined by w r (t) − r = 0 −∞ w r p(w r |r) dw r where w r is the particle-pair radial RV and p(w r |r) is the p.d.f. of radial RV conditioned on r.The magnitude of the measured MIRV normalized by the Kolmogorov velocity, u η , is presented in figure 9, both as a function of r/η (figure 9a,c,e,g,i) and as a function of r/a (figure 9b,d, f,h,j) for all the experimental conditions.For three of the cases (a = 20.75 μm and ρ = 0.5 g cm −3 , a = 14.25 μm and a = 3.75 μm) at Re λ = 324 we also plot data from available previous experiments (Dou et al. 2018) and DNS (Ireland et al. 2016), represented by the purple plus signs and purple dashed lines, respectively.
For the three cases with corresponding previous measurements and DNS data in figure 9(a,c,e,g,i), the agreement with previous experiments and DNS is excellent for separations larger than r/η = 20.The results all show an overlapping monotonic decreasing MIRV as the particle-pair separation decreases.This trend is predicted by DNS of one-way coupled point-particles subject to drag forces in turbulence (Ayala et al. 2008;Ireland et al. 2016;Dou et al. 2018;Hammond & Meng 2021).As the separation decreases below r/η = 20, all MIRV data continue to decrease, but the rate of decrease starts to diverge: our new data the slowest and DNS the fastest.At some smaller r/η values, our new MIRV starts to turn upward, which is not seen in previous measurements or DNS.The  turning point collapses onto a single r/a for each particle type, ranging from r/a ≈ 10 for the largest particles to r/a ≈ 50 for the smallest particles.

Collision kernel results
Using g(r) and w r (t) − r at r/a = 2.07 to approximate contact (r = 2a), we estimate the collision kernel for all the experimental conditions.Following Hammond & Meng (2021), the non-dimensional collision kernel can be calculated as The resulting non-dimensional collision kernel as a function of St is shown in figure 10, along with the values from DNS of one-way coupled, particle-laden turbulence (Ireland et al. 2016).Comparing the experimental collision kernel results with those of DNS, we find that our experimental estimates are O(10 3 -10 5 ) larger than DNS predictions.This is expected, since our values for the RDF at near contact are three to four orders of magnitude larger than predicted from DNS, and our values of the MIRV at contact are up to an order of magnitude larger than those from DNS as well.
The data point shown by a filled fuchsia square is at the identical condition to that previously reported in Hammond & Meng (2021).However, they originally made an error, accidentally reporting K * C (2a)/u η , which was subsequently corrected in a corrigendum (Hammond & Meng 2023).Here we show the corrected K * C (2a) = 3.6 × 10 4 for the St = 0.74 and a = 14.25 μm particles in the current study as well as in the corrigendum (Hammond & Meng 2023).

Discussion
In Hammond & Meng (2021) and Bragg et al. (2022), the authors argue that the surprising extreme clustering found in their experiments must be a result of particle-particle interactions.Bragg et al. (2022) started by hypothesizing that the particle-particle interactions responsible for extreme clustering were HIs, as did Yavuz et al. (2018), although the latter had several errors in their theory.Once the theory was corrected in Bragg et al. (2022), significant disagreement was found between the near-contact RDF predicted by the HI theory and all experimental results presented in the above three papers.
Hammond & Meng (2021) and Bragg et al. (2022) found that their measured RDF did in fact behave as g(r) − 1 ∝ r −6 as predicted by HI theory in the limit as St → 0. However, the magnitude of their experimental RDF at contact was at least two orders of magnitude larger than the theory prediction.Based mainly on the unexplainable large difference in the magnitude of g(r) predicted by the theory and that observed in the experiments, Bragg et al. (2022) concluded that the 'particle equation of motion invoked in the [HI] theory is clearly missing some vital effect, which future work must seek to uncover' and 'the mechanism for the extreme clustering observed here, in Hammond &Meng (2021), andYavuz et al. (2018) remains something of a mystery'.
To gain a better understanding of the extreme clustering phenomenon and the relevance of HI therein, we expanded the experimental conditions of Hammond & Meng (2021) and Bragg et al. (2022) by sweeping a broad range of St values.

The RDF and extreme clustering
Our RDF measurement data confirmed that the extreme clustering occurs consistently for all conditions, initiating around the start of the sub-Kolmogorov regime.For any given particle sample, extreme clustering starts at the same r/a for different fan speeds (and thus different Re λ and St).This indicates that extreme clustering is more directly related to particle properties than to turbulence.For any given particle sample, the magnitude of the RDF at contact increases monotonically with increasing St (i.e.fan speed), but across different particle samples (with varying sizes and densities) the influence of St is more complicated: the same St may correspond to different RDF values.
Among particles of different sizes and densities, the qualitative behaviour of the RDF across the entire measured r/a range was rather similar, showing four regions with distinct characteristics -inertial clustering (I), transition (II), extreme clustering (III) and decorrelation (IV) -as r/a decreases from large values to near contact.In each region, the RDF was influenced by particle properties in addition to St, with different particle samples having similar St behaving differently.

Comparison with HI theory
For convenience of discussion, recall (1.1), which is the RDF equation from the HI theory as derived in Bragg et al. (2022).This equation identifies four terms contributing to the RDF, associated with their μ i coefficients.While this equation was derived using the far-field (r/a 1) asymptotic forms of the mobility coefficients describing the HI, these are known to be valid down to r/a ≈ 2.05 and are very accurate for r/a ≥ 3 (Brunk et al. 1997).Therefore, any errors in the theory associated with the use of the far-field mobility coefficients will be negligible in the region of r/a ≥ 3, which includes the region where the RDF exhibits extreme clustering.The μ 1 term is the far-field form of the leading HI contribution that is independent of particle inertia and therefore not a function of St.The μ 2 and μ 3 terms represent the HI contributions that are dependent on particle inertia, i.e. the St effect on HI.The μ 4 term describes the clustering due to particle inertia interacting with turbulence, which happens even in the absence of HI (Chun et al. 2005).Since both μ 2 and μ 3 are negative (Bragg et al. 2022), the theory indicates that particle inertia should reduce the HI-induced mechanism for clustering (although as discussed in Bragg et al. (2022), the overall effect of increasing St is to increase the RDF due to the contribution of the non-HI term St 2 μ 4 ).
In the limit as St → 0, the HI-induced clustering mechanism in the theory of Bragg et al. (2022) is the same as that in Brunk et al. (1997).The basic mechanism is that the particles disturb the flow field surrounding other nearby particles in the flow, and when the strain rate associated with the disturbance flow is negative along the particle-pair separation direction, the particle velocity field will be compressed and the particles will cluster (Brunk et al. 1997;Bragg et al. 2022).This mechanism leads to the contribution g(r) − 1 ∼ (r/a) −6 in (1.1).The theory of Bragg et al. (2022) only accounts for the far-field HI; the near-field behaviour that dominates near contact includes the lubrication effect that inhibits actual particle contact when the fluid is modelled using the continuum hypothesis (Lambert, Weynans & Bergmann 2018;Ababaei et al. 2021).These near-field effects were included in Brunk et al. (1997), but the effects of particle inertia were not.In the limit as St → 0, (1.1) becomes and since the theory assumes that r/a 1, a Taylor expansion of the exponential function leads to . (5.2) Applying a Taylor expansion to (1.1) for r/a 1, we obtain the following for finite St: . (5.3) If we assume that HI explains the extreme clustering in our experiments when St 1, then we would expect that the measured g(r) primarily behaves according to the μ 1 (r/a) −(St 2 μ 4 +6) term.From Bragg et al. (2022), μ 4 is positive, which means that the scaling exponent in the extreme clustering range should be k = St 2 μ 4 + 6 ≥ 6 .However, figure 5 shows that our measured scaling exponent k is considerably smaller than 6 when St 1, reaching values as low as around 4.5.(Note that even though Bragg et al. (2022) reported a scaling exponent k = 6 with no St-dependent modulation, we have reprocessed their data in § 4 and found that k in their data is in the range 4.5 < k 6.6, with the single case presented in Hammond & Meng (2021) happening to be k = 6.) As previously discussed in Bragg et al. (2022), the discrepancy between the theory predictions and their experimental results cannot be explained by any of the assumptions made in the theory, as they are all valid for the experiments.Specifically, the HI theory assumes the following: (i) gravitational settling is negligible; (ii) the particles are point-particles; (iii) the particle Reynolds number, Re p , is small, leading to Stokes flow around the particles; (iv) other forces (such as added mass, pressure gradient etc.) are negligible; (v) particle-particle interactions are not many-body.These assumptions also hold for experiments in the current study.For (i), as mentioned in § 3.1, our particles' settling velocities are much too small for gravity to have a leading-order impact on the particle dynamics, and certainly far too small to explain the extreme clustering.For (ii), our particles have diameters from 7.5 μm to 41.5 μm, at least an order of magnitude smaller than the Kolmogorov lengths (101-179 μm).For (iii), all our experimental conditions have Re p = a 2 /τ n ν 1.For (iv), our particle-to-fluid density ratio is O(100), so we can consider other forces such as added mass, pressure gradient etc. to be small compared with the drag force on the particle (see Daitche 2015).Finally, for (v), as shown in § 5 of Bragg et al. (2022), more than half of our experimental conditions had an average of only one satellite particle around a primary particle in the extreme clustering region, indicating that these conditions did not exhibit many-body interactions.Our current experiments report all of the same experiments presented in Bragg et al. (2022) with the addition of more conditions, and all of our cases exhibited the extreme clustering.Therefore, the two-body interaction assumption cannot explain the discrepancy.Related to this, it is important to clarify that large values of the RDF associated with the extreme clustering do not imply that in the flow there are many particles clustered together in close proximity.High values of the RDF mean only that the probability of finding two particles at a given separation is much larger than it would be if the particles were uniformly distributed.It says nothing regarding the probability that there will be three or more particles in close proximity.
There is therefore overwhelming evidence that the HI theory completely fails to explain the extreme clustering observed throughout our experiments.Hence, we agree with Bragg et al. (2022) that 'new, yet-to-be-identified physical mechanisms are at play, requiring further investigation and new theories'.

Kinematics
Without knowledge of all the forces at play in our experiments, we cannot write a definitive dynamic equation for the particle pairs.Therefore, we instead sought to understand the kinematic effects governing g(r), and we identified two mechanisms contributing to the total drift of particle pairs, which must be negative to generate clustering.Our RDF data in both the inertial clustering (I) and the extreme clustering (III) regions were consistent with these drift mechanisms.The total drift in the extreme clustering region (III) is significantly negative, supporting that g(r) rapidly increases as r decreases.The PT term dominates in the total drift in region III, suggesting that the particle-pair motion underlying the extreme clustering is greatly influenced by the spatial gradient of the particle-pair RV fluctuations: particles prefer to move in directions of decreasing fluctuations.It is clear that further understanding the mechanisms of the PT term will provide meaningful insight into the particles' behaviour leading to the extreme clustering; however, interpreting this term depends on knowing the forces acting on the particles.
Additionally, we have observed that the RA is negative within the extreme clustering region (III), indicating that an inward, attractive force is acting to pull particles together.This force cannot be explained by HI, since our experiments lie outside the St → 0 limit, and modifications to HI due to particle inertia (from any non-zero value of St) should inhibit clustering, as discussed in § 5.2.
In the context of monodisperse inertial point-particles governed by Stokes drag forces with zero PPMF, the results in Bragg & Collins (2014a,b) show that the RA is positive, so that the RA term acts to inhibit particle clustering.This means that the Stokes drag force could not be responsible for the extreme clustering seen in our experiments.In the same context, PT corresponds to a drift in the direction of decreasing fluid velocity fluctuations, leading to a growth of the RDF as the particles approach.The underlying physical process associated with −∇ r S 2 can be understood in terms of a path-history, symmetry-breaking effect (Bragg & Collins 2014a;Bragg et al. 2015a,b).However, this physical interpretation of the effect captured by −∇ r S 2 does not necessarily apply to situations where the force acting on the particles is not a drag force.Indeed, as we saw when considering our experimental results, in the transition region (II), S 2 increases with decreasing r (i.e.−∇ r S 2 > 0), highlighting that the physical processes governing −∇ r S 2 in this region are not the same as those in the well-studied regime of point-particles subject to Stokes drag forces.
982 A21-23 5.4.The MIRV The MIRV between particle pairs is not only required for calculating the collision kernel but also important for understanding the behaviour of particles approaching each other.Our experimental findings on the MIRV in the sub-Kolmogorov regime are strikingly different from those previously reported in both experiments (Dou et al. 2018) and DNS (Ireland et al. 2016).In an earlier study of the same 'soccer ball' HIT chamber as in the current study, we performed planar four-frame particle tracking velocimetry measurements of particle RV using some of the same particles (Dou et al. 2018), resolving data down to separations of r/η = 1.There, higher MIRV than DNS predictions was also observed in the separations of r/η 10.Dou et al. (2018) conjectured that the deviation of experimental measurements from DNS may have been caused by either measurement uncertainty due to their thick two-dimensional laser sheet or the polydispersity in their particle samples.However, we suggest that neither is the case.
Our data show the same trends in the range 1 ≤ r/η ≤ 20 as observed in Dou et al. (2018), and in fact our results deviate even further from DNS.Our measurement system allows for much more accurate and 3-D measurements than available to Dou et al. (2018), so we can say that measurement uncertainty is not responsible for the deviations from DNS.As we have shown in the previous sections, the effects of polydispersity are most prominently seen in the decorrelation region (IV), but the deviations from DNS observed by Dou et al. (2018) occurred in the inertial clustering region (I), and our deviations occur throughout regions I-III.Additionally, the deviations occur for all our experimental conditions regardless of polydispersity level.More importantly, the particle samples with the lowest polydispersity deviate from DNS more significantly.So, polydispersity cannot be responsible for the observed deviations from DNS either.
Moreover, the MIRV in our measurements not only trended upward compared with DNS when 1 ≤ r/η ≤ 20 (where the data in Dou et al. (2018) reached their resolution limit) but continued to trend upward in the sub-Kolmogorov regime before reaching a maximum between 0.5 ≤ r/η ≤ 1 and trending downward again approaching contact.We suggest that the deviations from DNS seen in the current experiment and partly in the previous experiment by Dou et al. (2018) are the results of particle-particle interactions not accounted for in DNS.This speaks again for the proposition that another significant force is at play in experiments that has been ignored in theory and simulations.The consequence of ignoring this critical effect is clear when comparing our collision kernel with those of DNS, where our collision kernel reached values O(10 3 ) to O(10 5 ) larger than those calculated from DNS.

Conclusion
We have presented the most comprehensive 3-D particle tracking measurements to date of inertial particle clustering in isotropic turbulent flow at near-contact separations under 24 experimental conditions by broadly sweeping the Stokes number for inertial particles (0.07 ≤ St ≤ 3.68).We have examined the RDF, higher-order kinematics and MIRV, and consistently observed extreme clustering under all conditions.For all cases, the measured RDF reaches between O(10 2 ) and O(10 4 ) near contact, which is two to four orders of magnitude higher than HI theory predictions, thus confirming Bragg et al. (2022)'s observation.Although the measured RDF in the extreme clustering region in all cases roughly follows the (r/a) −6 trend predicted by the HI theory, the scaling exponent ranges from 4.5 to 7.6 in a St-dependent manner that the HI theory cannot explain.To gain insight into the kinematic mechanisms behind extreme clustering, we have developed a purely kinematic equation governing the RDF and examined contributions to the inward drift responsible for clustering.Analyses of the inward drift from experimental data consistently reveal that in the extreme clustering region (III), particle pairs experience an inward relative acceleration, which indicates the existence of an attractive force that cannot be attributed to the Stokes drag force or HIs.Moreover, the MIRV values at sub-Kolmogorov separations in all the cases are significantly elevated compared with those from previous DNS and experiments.All RDF and particle-pair kinematic quantities as functions of r/a reveal a dependence on particle properties beyond those captured by St, while the HI theory assumes, because of the equation of motion from which it is constructed, that in addition to r/a, only St should matter.Therefore, this study further strengthens the notion that there is a vital near-contact effect that greatly enhances particle clustering which existing DNS and theory do not capture.trends are still real, since the confidence windows for even the smallest separations do not overlap in MIRV and RA, and the plateaus and slope of the RDF are still significant.

A.3. Particle position uncertainty
The recordings of positions themselves have uncertainties which will impact the precision of the RDF and RV values.Through the use of the vibration isolation and volume self-calibration, we expect the position uncertainties captured to be within 0.15 pixels (Novara et al. 2019), which, based on our camera's pixel size of 21 μm, corresponds to a position uncertainty of δx = 3.2 μm and an uncertainty of δr ≈ 4.5 μm estimated via propagation (Moffat 1988).These uncertainties are small and overshadowed by the interpolation uncertainty.

A.4. Tracking input sensitivity and uncertainty
There are many user-defined parameters when running the Four-Pulse STB particle tracking algorithm, the most significant of which is the allowable triangulation error, .According to Novara et al. (2019), it is the most consequential parameter affecting the output.We used = 1.5 voxel (where 1 voxel ≈ 21 μm), since this value produced the most tracks when compared with surrounding values.To quantify the uncertainty in the RDF and RV calculations, was varied by ±10 % for all of the 50 Hz cases (40 Hz for a = 8.75 μm).We then took twice the standard deviation of these results as the vertical error bars shown in figures 2, 6(a,c,e,g,i), 8(a,c,e,g,i) and 9.The results show that the RV was not significantly impacted by triangulation error at small separations, but the RDF shows a large possibility for variations.These allowed variations do not change the order of magnitude of the results for the RDF in the 20.75 μm and 14.25 μm diameter cases but seem to cast some doubt in the 8.75 μm and 3.75 μm radius cases in the decorrelation region (IV).We still believe the overall trend of the RDF to be true for these cases, considering that the allowable variation within the extreme clustering region (III) still indicates the upward trend.

Figure 1 .
Figure 1.Diagram of the Four-Pulse STB system including HIT chamber, cameras, laser orientations (b), photograph of the HIT chamber and camera set-up (a) and timing scheme of the laser pulses and camera shutters (c).

Figure 2 .
Figure 2. Comprehensive RDF measurement results for five types of particles, presented as g(r) versus r/η (a,c,e,g,i) and r/a (b,d, f,h,j), consistently show extreme clustering.All axes are in log scale.In (b,d, f,h,j) the normalized particle separation distance r/a is divided into four regions (labelled I-IV) separated by vertical dotted lines: I, inertial clustering region; II, transition region; III, extreme clustering region (shaded); V, decorrelation region.Insets show detailed views of the transition region (II).The dashed line in the extreme clustering region (III) of each row represents the average fit in the region for all cases, labelled with r −k , where k is the scaling exponent averaged over all cases in that row.The error bars are from Re λ = 324 for all cases except a = 8.75 μm, which has error bars for the Re λ = 277 case.
Figure 3. Average value of the scaling exponent k in g(r) ∼ r −k over the separation r in the inertial clustering region (I) as a function of St with uncertainty error bars.The case presented in Hammond & Meng (2021) is represented as a filled fuchsia square.Linear fits are plotted for each particle type.The black plus signs are values from DNS (Ireland et al. 2016) averaged over Re λ .
7 for 0.5 < St < 1, and then slowly decreases as St continues to increase.The experimentally measured k value for our smallest particles (a = 3.75 μm, red circles) matches well with the DNS-predicted c 1 values under the same St, but for other particles the experimental k values are much smaller than the c 1 values predicted by DNS at the same St, and different particles have different k under the same St.However, we can discern a possible peak around St = 0.6 in the experiments, mirroring the trend of the DNS curve.The single case (St = 0.74) studied by Hammond & Meng (

Figure 4 .
Figure 4. Mean value of the gradient of k (second-order gradient of g(r)) in the transition region (II), which represents the swiftness of transition, versus St.The case presented in Hammond & Meng (2021) is represented as a filled fuchsia square.Linear fit lines are plotted for each particle type.

Figure 5 .
Figure 5. Average value of the scaling exponent k in g(r) ∼ r −k over the separation r in the extreme clustering region (III) as a function of St including vertical error bars as the uncertainty in k.The case presented in Hammond & Meng (2021) is represented as a filled fuchsia square.Linear fit lines of k versus St are plotted as solid lines for each particle type.

Figure 6 .
Figure 6.The RDF (a,c,e,g,i) and total drift ẇr (t) r − ∇ r S 2 normalized by the Kolmogorov acceleration (b,d, f,h,j) presented as functions of r/a.The vertical dotted lines represent the boundaries of the regions described in § 4.1, the shaded region is the extreme clustering region (III), and the horizontal dashed red lines in (b,d, f,h,j) indicate where the total drift is equal to zero.Panels (a,b), (c,d), (e, f ), (g,h) and (i,j) are for different particle types, and the different-coloured curves represent different St conditions.The error bars in (a,c,e,g,i) are from Re λ = 324 for all cases except a = 8.75 μm, which has error bars for the Re λ = 277 case.

Figure 7 .
Figure 7.The normalized PPMF presented as a function of r/a.The vertical dotted lines represent the boundaries of the regions described in § 4.1, and the shaded region is the extreme clustering region (III).Each plot is for a different particle type and the different-coloured curves represent different St conditions.The error bars are from Re λ = 324 for all cases except a = 8.75 μm, which has error bars for the Re λ = 277 case.Note that for only the a = 20.75 μm case, the y-axis ranges from −5000 to 5000 instead of −1000 to 1000.

Figure 8 .
Figure 8.The two contributions in the total drift term normalized by the Kolmogorov acceleration presented as a function of r/a: (a,c,e,g,i) mean particle-pair radial RA ẇr (t) r ; (b,d, f,h,j); PT −∇ r S 2 .The vertical dotted lines represent the boundaries of the regions described in § 4.1, the shaded region is the extreme clustering region (III) and the horizontal dashed red lines indicate where RA and PT are equal to zero.Each row is for a different particle type and the different-coloured curves represent different St conditions.The error bars in (a,c,e,g,i) are from Re λ = 324 for all cases except a = 8.75 μm, which has error bars for the Re λ = 277 case.

Figure 9 .
Figure 9. Comprehensive MIRV, presented as St decreases in general from (a,b) to (i,j), versus r/η (a,c,e,g,i) and r/a (b,d, f,h,j).Vertical dotted lines in (b,d, f,h,j) show division of regions.The shaded region in (b,d, f,h,j) is the extreme clustering region (III).The error bars are from Re λ = 324 for all cases except a = 8.75 μm, which has error bars for the Re λ = 277 case.Panels (a,b), (e,f ) and (i,j) include DNS from Ireland et al. (2016) and previous experimental results from Dou et al. (2016).

DNSFigure 10 .
Figure 10.Normalized collision kernel as a function of St.The case presented in Hammond & Meng (2021) is represented as a filled fuchsia square.Linear fit lines are plotted for each particle type.Here DNS values from Ireland et al. (2016) are shown averaged over Re λ .

Figure 11 .
Figure 11.Measurement uncertainty for (a,d,g,j,m) the RDF, (b,e,h,k,n) the MIRV and (c, f,i,l,o) the RA.Blue horizontal error bars represent the separation uncertainty, and the black vertical error bars represent the uncertainty in the measured values which have been shown previously in the respective figures.Only one case per particle type is shown for clarity.Error bars from Re λ = 324 are shown for all cases except a = 8.75 μm, for which Re λ = 277 is shown.Dotted vertical lines indicate the separation of the regions I-IV.The RA plots have a horizontal dashed red line at zero for reference.

Table 1 .
Hammond & Meng (2021)ns.STB inputs were used for all conditions.Flow conditions and particle properties are combined to calculate the St for each condition.The condition presented inHammond & Meng (2021)is highlighted in dark grey, and the cases presented inBragg et al. (2022)are highlighted in light grey and include the case presented inHammond & Meng (2021).New cases are left white.*Themeasurement data for the a = 8.75 μm and Re λ = 324 condition, corresponding to St = 0.49, was inadvertently deleted and unable to be recovered.

Table 2 .
Experimental statistics: for each experimental condition, the number of turbulence realizations (images) captured and the average number of particles per frame in those images, the maximum RSE for each case's RDF calculation and in which bin (r/a) it occurred, and the maximum uncertainty in the separation calculation δr in , normalized by a.