1. Introduction
Understanding how particles collide in isotropic turbulence is critical for modelling droplet coalescence in cloud turbulence (Shaw Reference Shaw2003; Grabowski & Wang Reference Grabowski and Wang2013), particle coagulation in the chemical and process industries (McMillan et al. Reference McMillan, Shaffer, Gopalan, Chew, Hrenya, Hays, Karri and Cocco2013) and mass entrainment within combustion engines (Wei et al. Reference Wei, Yu, Grout, Chen and Ma2011). The rate of collisions in an isotropic flow can be described by the collision kernel (Sundaram & Collins Reference Sundaram and Collins1997; Wang, Wexler & Zhou Reference Wang, Wexler and Zhou2000), which depends on two statistics: the radial distribution function (RDF), 
$g(r)$, and the mean of the inward particle-pair radial relative velocity (RV) 
$\langle w_r \rangle ^-$, 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 Reference Maxey1987; Chun et al. Reference Chun, Koch, Rani, Ahluwalia and Collins2005; Bec et al. Reference Bec, Biferale, Cencini, Lanotte, Musacchio and Toschi2007; Gustavsson & Mehlig Reference Gustavsson and Mehlig2011; Bragg & Collins Reference Bragg and Collins2014a; Bragg, Ireland & Collins Reference Bragg, Ireland and Collins2015a,Reference Bragg, Ireland and Collinsb; Gustavsson & Mehlig Reference Gustavsson and Mehlig2016). In the past, experimental measurements of the RDF (Salazar et al. Reference Salazar, De Jong, Cao, Woodward, Meng and Collins2008) and RV (Dou et al. Reference Dou, Ireland, Bragg, Liang, Collins and Meng2018), resolving particle separations (
$r$) down to the Kolmogorov scale (
$\eta$), 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\unicode{x2013}10^1)$ as 
$r$ decreases down to 
$O(\eta )$. However, recent advances in laser particle tracking techniques have enabled measurements at sub-Kolmogorov separations down to near-contact scales (Yavuz et al. Reference Yavuz, Kunnen, Van Heijst and Clercx2018; Hammond & Meng Reference Hammond and Meng2021; Bragg et al. Reference Bragg, Hammond, Dhariwal and Meng2022) and revealed surprising extreme clustering of inertial particles never seen before. As 
$r$ decreases below 
$O(\eta )$, 
$g(r)$ grows from 
$O(1)$ to 
$O(10^3)$, scaling as 
$r^{-k}$ where 
$k\approx 6$ (Hammond & Meng Reference Hammond and Meng2021; Bragg et al. Reference Bragg, Hammond, Dhariwal and Meng2022). Hammond & Meng (Reference Hammond and Meng2021) 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. Reference Wang, Ayala, Kasprzak and Grabowski2005; Dhanasekaran, Roy & Koch Reference Dhanasekaran, Roy and Koch2021). Yavuz et al. (Reference Yavuz, Kunnen, Van Heijst and Clercx2018) attempted to explain their surprising extreme clustering data by extending the drift–diffusion model of Chun et al. (Reference Chun, Koch, Rani, Ahluwalia and Collins2005) to describe the clustering of weakly inertial particles with 
$St\ll 1$ subject to HI, where the Stokes number 
$St$ is defined as the particle response time 
$\tau _p$ divided by the Kolmogorov time scale 
$\tau _\eta$. However, Bragg et al. (Reference Bragg, Hammond, Dhariwal and Meng2022) found that the theory of Yavuz et al. (Reference Yavuz, Kunnen, Van Heijst and Clercx2018) contained several errors. Once these errors are corrected, the HI theory predicts that the RDF should behave as (Bragg et al. Reference Bragg, Ireland and Collins2015b)
\begin{equation} g(r) \sim (r/a)^{{-}St^2 \mu_4} \exp(\mu_1 (r/a)^{{-}6} + (St \mu_2 + St^2 \mu_3 ) (r/a)^{{-}1} ), \end{equation}
where 
$r$ is the particle-pair separation, 
$a$ is the particle radius and 
$\mu _i$ are coefficients that depend on the flow properties (the detailed definitions of 
$\mu _i$ are discussed in § 5.2 and can be found in more detail in Bragg et al. (Reference Bragg, Hammond, Dhariwal and Meng2022)). The leading HI term 
$\exp (\mu _1 (r/a)^{-6})$ is the far-field form of the result derived in Brunk, Koch & Lion (Reference Brunk, Koch and Lion1997) describing the clustering due to HI independent of 
$St$. In the limit as 
$St \rightarrow 0$, one obtains the far-field relation 
$g(r)-1 \sim (r/a)^{-6}$ as the leading contribution of HI to clustering.
Both Hammond & Meng (Reference Hammond and Meng2021) and Bragg et al. (Reference Bragg, Hammond, Dhariwal and Meng2022) 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. (Reference Bragg, Hammond, Dhariwal and Meng2022) 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. (Reference Bragg, Hammond, Dhariwal and Meng2022) 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 (Reference Hammond and Meng2021), and Yavuz et al. (Reference Yavuz, Kunnen, Van Heijst and Clercx2018) 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. Reference Bragg, Hammond, Dhariwal and Meng2022).
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 Reference Hammond and Meng2021; Bragg et al. Reference Bragg, Hammond, Dhariwal and Meng2022) and RV (Hammond & Meng Reference Hammond and Meng2021) 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\leq St \leq 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 Reference Bragg and Collins2014a), 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.
2. 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. (Reference Yavuz, Kunnen, Van Heijst and Clercx2018), Hammond & Meng (Reference Hammond and Meng2021) and Bragg et al. (Reference Bragg, Hammond, Dhariwal and Meng2022).
Let 
${\boldsymbol {r}}(t) \text{ and } {\boldsymbol {w}}(t)$ denote the particle-pair separation (based on the particle centres) and RV vectors, and 
$w_r (t)\equiv \Vert\boldsymbol {r}(t)\Vert^{-1} \boldsymbol {r}(t)\boldsymbol {\cdot }{\boldsymbol {w}}(t)$ the radial (or longitudinal) component of the RV. Kinematically, we have 
$\dot {r}(t)\equiv w_r (t)$, where 
$r(t)\equiv \Vert\boldsymbol {r}(t)\Vert$. 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. 
$\varrho (r,t)\equiv \langle \delta (r(t)-r)\rangle$ may be derived (Bragg et al. Reference Bragg, Ireland and Collins2015b; Johnson, Bassenne & Moin Reference Johnson, Bassenne and Moin2020),
\begin{gather} \frac{{\rm D}}{{\rm D}t} \varrho={-}\varrho\boldsymbol{\nabla}_r \langle w_r(t)\rangle_r, \end{gather}
\begin{gather}\varrho \frac{{\rm D}}{{\rm D}t}\langle w_r(t)\rangle_r={-}\boldsymbol{\nabla}_r (\varrho S_2 )+\varrho \langle \dot{w}_r (t)\rangle _r, \end{gather}
where 
${{\rm D}}/{{\rm D}t}\equiv \partial _t+\langle w_r (t)\rangle _r \boldsymbol {\nabla }_r$, 
$r$ is a time-independent coordinate, 
$\langle \cdot \rangle _r$ denotes an ensemble average conditioned on 
$r(t)=r$ and 
$S_2\equiv \langle [w_r (t)-\langle w_r (t)\rangle _r ]^2\rangle _r$ is the second-order particle RV structure function.
The time-dependent RDF 
$g(r,t)$ is related to 
$\varrho$ via a constant 
$\beta$ as 
$\varrho =\beta g$ (Chun et al. Reference Chun, Koch, Rani, Ahluwalia and Collins2005; Bragg et al. Reference Bragg, Ireland and Collins2015b), and substituting this into (2.1) and rearranging gives an exact drift–diffusion equation governing the RDF:
\begin{equation} 0={-}S_2 \boldsymbol{\nabla}_r g+g\left(\langle \dot{w}_r (t)\rangle_r-\boldsymbol{\nabla}_r S_2-\frac{{\rm D}}{{\rm D}t} \langle w_r (t)\rangle_r\right). \end{equation}For a statistically stationary flow, (1.1) yields
\begin{equation} g\langle w_r (t)\rangle_r=\varPhi, \end{equation}
where 
$\varPhi$ 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. (Reference Bragg, Hammond, Dhariwal and Meng2022)) and no injection of particle pairs at the large scales of the flow, a zero-flux state is expected with 
$\varPhi =0$. In our experiments 
$g>0$ (regions where 
$g=0$ do not occur), which then implies 
$\langle w_r (t)\rangle _r=0$ from (2.4), and so (2.3) reduces to
\begin{equation} 0={-} S_2 \boldsymbol{\nabla}_r g+g(\langle \dot{w}_r (t)\rangle_r-\boldsymbol{\nabla}_r S_2 ). \end{equation}
Since in a turbulent flow 
$S_2 > 0$, the sign of 
$\boldsymbol {\nabla }_r g$ is governed entirely by the sign of the drift term 
$\langle \dot {w}_r (t)\rangle _r-\boldsymbol {\nabla }_r S_2$, which is not sign definite. In regimes where 
$g(r)$ grows as 
$r$ reduces, i.e. 
$\boldsymbol {\nabla }_r g<0$, we must have
\begin{equation} \langle \dot{w}_r (t)\rangle_r-\boldsymbol{\nabla}_r S_2<0, \end{equation}
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, 
$\langle \dot {w}_r (t)\rangle _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, 
$\langle \dot {w}_r (t)\rangle _r<0$, the forces on particle pairs on average act in the negative 
$r$ direction (meaning inward), thus promoting clustering.
The second term, 
$-\boldsymbol {\nabla }_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, 
$-\boldsymbol {\nabla }_r S_2<0$, it contributes to an inward drift of particles and thus promotes clustering. This condition corresponds to 
$\boldsymbol {\nabla }_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 Reference Reeks1983; Johnson et al. Reference Johnson, Bassenne and Moin2020; Bragg, Richter & Wang Reference Bragg, Richter and Wang2021), 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 
$\langle (\dot {w}_r )(t)\rangle _r$ and 
$-\boldsymbol {\nabla }_r S_2$ could help us gain insight into the physical processes generating the extreme clustering phenomenon.
3. Experiments
3.1. 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 Reference Sellappan, Alvi and Cattafesta2020) described by Hammond & Meng (Reference Hammond and Meng2021). 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. Reference Dou, Pecenak, Cao, Woodward, Liang and Meng2016). The chamber was electrically grounded and the net electric charge on particles measured in situ was found to be negligible (Bragg et al. Reference Bragg, Hammond, Dhariwal and Meng2022). 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. (Reference Dou, Pecenak, Cao, Woodward, Liang and Meng2016). 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.
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).
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_\lambda < 357$ and 
$0.07 < St < 3.68$. Among these 25 conditions (with distinct 
$St$ values given), the 
$a=14.25\ \mathrm {\mu } \mathrm {m} \text{ and } St = 0.74$ condition (dark grey) was previously presented in Hammond & Meng (Reference Hammond and Meng2021), and an additional 12 cases (light grey) were studied in Bragg et al. (Reference Bragg, Hammond, Dhariwal and Meng2022). The current study adds 12 new cases (white) to complete the matrix of conditions. Unfortunately, one set of data 
$(a=8.75\ \mathrm {\mu } \mathrm {m} , Re_\lambda =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.
Table 1. Experimental conditions. STB inputs were used for all conditions. Flow conditions and particle properties are combined to calculate the 
$St$ for each condition. The condition presented in Hammond & Meng (Reference Hammond and Meng2021) is highlighted in dark grey, and the cases presented in Bragg et al. (Reference Bragg, Hammond, Dhariwal and Meng2022) are highlighted in light grey and include the case presented in Hammond & Meng (Reference Hammond and Meng2021). New cases are left white.
*The measurement data for the 
$a=8.75\ \mathrm {\mu } \mathrm {m}$ and 
$Re_{\lambda }=324$ condition, corresponding to 
$St=0.49$, was inadvertently deleted and unable to be recovered.
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\ \mathrm {\mu } \mathrm {m}$ to 
$105\ \mathrm {\mu } \mathrm {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. (Reference Dou, Ireland, Bragg, Liang, Collins and Meng2018). From the sieved particles we utilized five particle samples in four particle diameter ranges (
$5\unicode{x2013}10\ \mathrm {\mu } \mathrm {m}, 15\unicode{x2013}20\ \mathrm {\mu } \mathrm {m}, 25\unicode{x2013}32\ \mathrm {\mu } \mathrm {m}$ and 
$38\unicode{x2013}45\ \mathrm {\mu } \mathrm {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\ \mathrm {\mu } \mathrm {m}$ or 
$7\ \mathrm {\mu } \mathrm {m}$ with variance in radius being either 
$\pm 1.25 \ \mathrm {\mu } \mathrm {m}$ or 
$\pm 1.75\ \mathrm {\mu } \mathrm {m}$, respectively. As 
$a$ decreases, this variance became increasingly significant compared with the mean particle radius, thereby increasing particle polydispersity. The polydispersity 
$(\phi =2\sigma _a/a)$ is defined here as the range of particle radii 
$2\sigma _a$ in a final sample divided by the average particle radius a in that sample.
3.2. Optical diagnostics
The Four-Pulse STB particle tracking velocimetry set-up and method were described by Hammond & Meng (Reference Hammond and Meng2021). 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^{\circ }$ from the centre plane in each cardinal direction. The laser beams travel through a series of optics, resulting in a 
$10\ {\rm mm} \times 30\ {\rm mm}$ laser sheet entering the chamber. The actual measurement volume utilized is 
$10\ {\rm mm} \times 30\ {\rm mm} \times 50\ {\rm mm}$ in the centre of the chamber. In addition, the temperature inside the flow chamber was monitored during operation and kept within 
$\pm 1\,^\circ {\rm C}$.
3.3. Particle tracking and processing
After the injected particles were allowed to equilibrate over approximately 100 large-eddy turnover times (
$\approx$30 s), particle tracking was performed on three independent sets of 3093 realizations of the particle-laden flow for the 
$a=3.75\ \mathrm {\mu } \mathrm {m} , a=8.75\ \mathrm {\mu } \mathrm {m}$ and 
$a=20.75\ \mathrm{\mu } \mathrm {m}$ particles, and five independent sets of 3093 realizations for the 
$a=14.25\ \mathrm {\mu } \mathrm {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 
$\Delta t_1, \Delta t_2$ and 
$\Delta t_3$ as shown in figure 1. The setting of these time intervals varied with 
$Re_\lambda$ as shown in table 1. The tracks were then processed using the STB function in DaVis. All 
$\Delta t$ values and STB input parameters are listed in table 1.
The particle volume fraction was kept at approximately 
$2.2 \times 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. (Reference Bragg, Hammond, Dhariwal and Meng2022) 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 
$\langle w_r (t)\rangle _r^{-}$, as well as two components of the drift mechanism, 
$\langle \dot {w}_r (t)\rangle _r$ and 
$-\boldsymbol {\nabla }_r S_2$. We accounted for the physical boundary of the experiment in calculating the RDF as described in Hammond & Meng (Reference Hammond and Meng2021). 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.
3.4. Measurement uncertainties
Measurement uncertainty was calculated in much the same way as Hammond & Meng (Reference Hammond and Meng2021). 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 
$\Delta 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, 
$\epsilon$, is the most important (Novara et al. Reference Novara, Schanz, Geisler, Gesemann, Voss and Schröder2019). By varying 
$\epsilon$ by 
$\pm 10\,\%$ for our 
$50\ {\rm Hz}$ cases, we produce the vertical error bars present in the figures throughout the paper.
4. Results
4.1. 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/\eta$ 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/\eta$ values for different types of particles and fan speeds. Note that 
$r/a = 2.00$ is contact for monodisperse particles.
Figure 2. Comprehensive RDF measurement results for five types of particles, presented as 
$g(r)$ versus 
$r/\eta$ (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_{\lambda } = 324$ for all cases except 
$a = 8.75\ \mathrm {\mu } {\rm m}$, which has error bars for the 
$Re_{\lambda } = 277$ case.
All 24 sets of 
$g(r)$ data consistently show extreme clustering, with 
$g(r)$ reaching 
$O(10^2)\unicode{x2013}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 \sim (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 =\infty$ to 
$r/a=2.07$ near contact, beginning in region I, the inertial clustering region.
4.1.1. Inertial clustering region (I)
From 
$r/\eta =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 has been extensively studied in theory (Chun et al. Reference Chun, Koch, Rani, Ahluwalia and Collins2005; Zaichik & Alipchenkov Reference Zaichik and Alipchenkov2007; Bragg & Collins Reference Bragg and Collins2014a), DNS (Saw et al. Reference Saw, Salazar, Collins and Shaw2012a; Rosa et al. Reference Rosa, Parishani, Ayala, Grabowski and Wang2013; Ireland, Bragg & Collins Reference Ireland, Bragg and Collins2016; Dhanasekaran & Koch Reference Dhanasekaran and Koch2022) and experiments (Salazar et al. Reference Salazar, De Jong, Cao, Woodward, Meng and Collins2008). 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 \text{ to } 0.22$, and the data for 
$g(r)$ is consistent with the power-law scaling 
$(r/\eta )^{-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 Reference Reade and Collins2000; Chun et al. Reference Chun, Koch, Rani, Ahluwalia and Collins2005; Saw et al. Reference Saw, Salazar, Collins and Shaw2012a,Reference Saw, Shaw, Salazar and Collinsb). 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/\eta < 2.7$ depending on particle type and 
$Re_\lambda$ 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_\lambda$ from DNS of non-interacting small and heavy inertial particles in isotropic turbulence (Ireland et al. Reference Ireland, Bragg and Collins2016), denoted by the black plus signs. Each experimental particle type (unique combination of radius, 
$a$, density, 
$\rho$, and polydispersity, 
$\phi$) 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 \approx 0.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\ \mathrm {\mu } {\rm 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.
Figure 3. Average value of the scaling exponent 
$k$ in 
$g(r) \sim 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 (Reference Hammond and Meng2021) 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. Reference Ireland, Bragg and Collins2016) averaged over 
$Re_{\lambda }$.
The single case 
$(St = 0.74)$ studied by Hammond & Meng (Reference Hammond and Meng2021) 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 (Reference Hammond and Meng2021) 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. Reference Saw, Salazar, Collins and Shaw2012a,Reference Saw, Shaw, Salazar and Collinsb). However, in figure 3, we see no clear trend for 
$k$ as the polydispersity, 
$\phi$, 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\ \mathrm {\mu } \mathrm {m}$) and blue diamond case (
$a = 8.75\ \mathrm {\mu } \mathrm {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\ \mathrm {\mu } \mathrm {m}$, with 
$\rho = 0.5\ \mathrm {g}\ {\rm cm}^{-3}$ and 
$0.3\ \mathrm {g}\ {\rm 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\ \mathrm {\mu } \mathrm {m}$) and green triangle (
$a = 20.75\ \mathrm {\mu } \mathrm {m})$ cases at the same 
$St$ (
$St = 0.74,\text{ with } k = 0.14$ and 
$k = 0.09$, respectively), and even more so for the blue diamond (
$a = 8.75\ \mathrm {\mu } \mathrm {m}$) and fuchsia square (
$a = 14.25\ \mathrm {\mu } \mathrm {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\ \mathrm {\mu }\mathrm {m}$ and 
$8.75\ \mathrm {\mu }\mathrm {m}$), but decreases as 
$St$ increases for the larger particles (
$a = 14.25\ \mathrm {\mu }\mathrm {m}$ and 
$20.75\ \mathrm {\mu }\mathrm {m}$). Besides particle size, some other particle properties could also be playing a role in the inertial clustering. The two 
$a=20.75\ \mathrm {\mu }\mathrm {m}$ cases (green triangle and yellow star) present drastically different 
$k$ magnitudes, where the denser particle type (
$\rho =0.5\ {\mathrm {g}}\ \mathrm {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 (Reference Daitche2015), where particles governed by the full Maxey & Riley (Reference Maxey and Riley1983) equation were simulated in isotropic turbulence. Daitche (Reference Daitche2015) showed that the clustering of the particles can depend strongly on the particle-to-fluid density ratio 
$\varrho \equiv \rho _p/\rho _f$ as well as on 
$a/\eta$ (within the 
$a/\eta \ll 1$ regime). However, for the values of 
$\varrho$ and 
$a/\eta$ in our experiments, the DNS results of Daitche (Reference Daitche2015) would indicate that the effects of 
$\varrho$ and 
$a/\eta$ are weak (e.g. see figure 8 in Daitche (Reference Daitche2015)). Hence the dependence we see of 
$k$ on 
$\varrho$ and 
$a/\eta$, in addition to 
$St$, cannot be explained within the context of the additional forces in the Maxey & Riley (Reference Maxey and Riley1983) equation that are important when 
$\varrho$ is not sufficiently large or 
$a/\eta$ sufficiently small.
4.1.2. 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 
$\langle \boldsymbol {\nabla }_r k\rangle$) versus 
$St$ in figure 4. Here 
$\langle \boldsymbol {\nabla }_r k\rangle$ represents the swiftness of transition, showing how quickly a given case experiences the transition from the inertial clustering region (I) to the extreme clustering region (III) since it represents how quickly 
$k$ is changing, or the change in slope of the RDF. A large value of 
$\langle \boldsymbol {\nabla }_r k\rangle$ represents a swift and abrupt transition (
$k$ is changing quickly as a function of 
$r$), while a low value of 
$\langle \boldsymbol {\nabla }_r k\rangle$ is indicative of a slow, smooth transition (
$k$ is changing slowly as a function of 
$r$). Understanding how swiftly a given curve transitions can ultimately provide insight into how competing forces exchange dominance when transitioning between regions I and III, where the clustering is potentially driven by fundamentally different physics.
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 (Reference Hammond and Meng2021) is represented as a filled fuchsia square. Linear fit lines are plotted for each particle type.
Figure 4 shows a weak increase of 
$\langle \boldsymbol {\nabla }_r k\rangle$ 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\ \mathrm {\mu }\mathrm {m} ,\rho = 0.3\ \mathrm {g}\ \mathrm {cm}^{-3}$) and yellow star (
$a = 20.75\ \mathrm {\mu }\mathrm {m} ,\rho = 0.5\ \mathrm {g}\ \mathrm {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 
$\langle \boldsymbol {\nabla }_r k\rangle$. This indicates that particle properties may also be contributing to the swiftness of transition between the inertial clustering (I) and extreme clustering (III) regions.
4.1.3. 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/\eta \approx 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 Reference Hammond and Meng2021). 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/\eta$ 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 
$\phi$ 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 we plot the 
$k$ value for each individual case. The full range of 
$k$ was approximately 
$4.6 < k < 7.6$ in this region, increasing as 
$St$ increases. Note that Hammond & Meng (Reference Hammond and Meng2021) define the 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. (Reference Bragg, Hammond, Dhariwal and Meng2022) did not differentiate between the 
$k$ values and reported 
$g(r) \sim r^{-6}$ for all their cases.
Figure 5. Average value of the scaling exponent 
$k$ in 
$g(r) \sim 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 (Reference Hammond and Meng2021) is represented as a filled fuchsia square. Linear fit lines of 
$k$ versus 
$St$ are plotted as solid lines for each particle type.
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\ \mathrm {\mu }\mathrm {m}$) and the blue diamond (
$a = 8.75\ \mathrm {\mu }\mathrm {m}$) cases, as the radius increases, 
$k$ decreases in general. But comparing the blue diamond and the fuchsia square (
$a=14.25\ \mathrm {\mu }\mathrm {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\ \mathrm {\mu }\mathrm {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 (
$\rho = 0.3\ \mathrm {g}\ \mathrm {cm}^{-3}$ and 
$0.5\ \mathrm {g}\ \mathrm {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\ \mathrm {\mu }\mathrm {m}$) and blue diamond case (
$a = 8.75\ \mathrm {\mu }\mathrm {m}$) under the same 
$St$ (0.23) present quite different 
$k$ values (
$k = 5.69$ and 4.96, respectively), while the fuchsia square (
$a= 14.25\ \mathrm {\mu }\mathrm {m}$) and green triangle (
$a = 20.75\ \mathrm {\mu }\mathrm {m}$) cases have much closer 
$k$ values at the same 
$St$ (0.74) (
$k = 6.07$ and 
$k = 5.90$, respectively).
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.
4.1.4. 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. Reference Chun, Koch, Rani, Ahluwalia and Collins2005; Saw et al. Reference Saw, Shaw, Salazar and Collins2012b; Dhariwal & Bragg Reference Dhariwal and Bragg2018). In brief, samples that have a relatively wide variation (from polydispersity, 
$\phi$) 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, 
$\phi$, of the sample and unrelated to any particle–turbulence interactions.
4.2. Particle-pair kinematics results
As previously discussed in § 2, (2.6) shows that there are two distinct kinematic effects that are exclusively associated with particle clustering and the growth of the RDF, 
$\langle \dot {w}_r (t)\rangle _r$ and 
$-\boldsymbol {\nabla }_r S_2$. For 
$g(r)$ to increase as the separation decreases, the total drift term in (2.6) must be negative, i.e. 
$\langle \dot {w}_r (t)\rangle _r-\boldsymbol {\nabla }_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 (
$\langle \dot {w}_r (t)\rangle _r$) and PT (
$-\boldsymbol {\nabla }_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$.
4.2.1. Total drift
In figure 6, we present the total drift 
$\langle \dot {w}_r (t)\rangle _r-\boldsymbol {\nabla }_r S_2$ normalized by the Kolmogorov acceleration, 
$a_\eta$ (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. The RDF (a,c,e,g,i) and total drift 
$\langle \dot {w}_r (t) \rangle _r - \boldsymbol {\nabla }_rS_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_{\lambda } = 324$ for all cases except 
$a = 8.75\ \mathrm {\mu }\mathrm {m}$, which has error bars for the 
$Re_{\lambda } = 277$ case.
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_\eta$. 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:
\begin{equation} PPMF= g(r)(\langle w_r(t)\rangle_r/u_\eta). \end{equation}These plots include error bars determined via propagation using the uncertainty for the RDF and radial RV (Moffat Reference Moffat1988). 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).
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_{\lambda } = 324$ for all cases except 
$a = 8.75\ \mathrm {\mu }\mathrm {m}$, which has error bars for the 
$Re_{\lambda } = 277$ case. Note that for only the 
$a = 20.75\ \mathrm {\mu }\mathrm {m}$ case, the 
$y$-axis ranges from 
$-5000$ to 5000 instead of 
$-1000$ to 1000.
Interestingly, the PPMF deviates from zero at smaller separations. For all cases except the 
$a=20.75\ \mathrm {\mu }\mathrm {m} \text{ and } \rho =0.5\ \mathrm {g}\ \mathrm {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\ \mathrm {\mu }\mathrm {m} \text{ and } \rho =0.5\ \mathrm {g}\ \mathrm {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. (Reference Dou, Pecenak, Cao, Woodward, Liang and Meng2016); 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. (Reference Dou, Pecenak, Cao, Woodward, Liang and Meng2016) 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, 
$\langle \dot {w}_r (t) \rangle _r$, normalized by the Kolmogorov acceleration, 
$a_\eta$, 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, 
$\phi$, experience lower RA magnitudes and an earlier onset of a plateau in RA near contact, mirroring the plateau observed in the RDF.
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 
$\langle \dot {w}_r (t) \rangle _r$; (b,d,f,h,j); PT 
$-\boldsymbol {\nabla }_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_{\lambda } = 324$ for all cases except 
$a = 8.75\ \mathrm {\mu }\mathrm {m}$, which has error bars for the 
$Re_{\lambda } = 277$ case.
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, 
$-\boldsymbol {\nabla }_r S_2$, normalized by the Kolmogorov acceleration, 
$a_\eta$, 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.
4.3. 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 
$\langle w_r (t)\rangle _r^-=\int _{-\infty }^0 w_r p(w_r|r)\,{\rm d}w_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_\eta$, is presented in figure 9, both as a function of 
$r/\eta$ (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\ \mathrm {\mu }\mathrm {m} \text{ and } \rho = 0.5\ \mathrm {g}\ \mathrm {cm}^{-3} , a = 14.25\ \mathrm {\mu }\mathrm {m}$ and 
$a =3.75\ \mathrm {\mu }\mathrm {m}$) at 
$Re_\lambda = 324$ we also plot data from available previous experiments (Dou et al. Reference Dou, Ireland, Bragg, Liang, Collins and Meng2018) and DNS (Ireland et al. Reference Ireland, Bragg and Collins2016), represented by the purple plus signs and purple dashed lines, respectively.
Figure 9. Comprehensive MIRV, presented as 
$St$ decreases in general from (a,b) to (i,j), versus 
$r/\eta$ (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_{\lambda } = 324$ for all cases except 
$a = 8.75\ \mathrm {\mu }\mathrm {m}$, which has error bars for the 
$Re_{\lambda } = 277$ case. Panels (a,b), (e,f) and (i,j) include DNS from Ireland et al. (Reference Ireland, Bragg and Collins2016) and previous experimental results from Dou et al. (Reference Dou, Pecenak, Cao, Woodward, Liang and Meng2016).
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/\eta = 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. Reference Ayala, Rosa, Wang and Grabowski2008; Ireland et al. Reference Ireland, Bragg and Collins2016; Dou et al. Reference Dou, Ireland, Bragg, Liang, Collins and Meng2018; Hammond & Meng Reference Hammond and Meng2021). As the separation decreases below 
$r/\eta =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/\eta$ 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\approx 10$ for the largest particles to 
$r/a\approx 50$ for the smallest particles.
The surprising upturn of MIRV happens consistently for all our experimental conditions, reaching a local maximum in the extreme clustering region (III) and then decreasing again as separations approach contact. Like the RDF and RA, samples with higher polydispersity experience a plateau in the decorrelation region (IV). For all cases, the upturn occurs near 
$r/\eta =1$, suggesting that the cause of this upturn is less related to particle–turbulence interactions and is more likely caused by particle–particle interactions.
4.4. Collision kernel results
Using 
$g(r)$ and 
$\langle w_r (t)\rangle _r^-$ at 
$r/a=2.07$ to approximate contact (
$r=2a$), we estimate the collision kernel for all the experimental conditions. Following Hammond & Meng (Reference Hammond and Meng2021), the non-dimensional collision kernel can be calculated as
\begin{equation} K_C^*(2a)=\frac{K_C(2a)}{(2a)^2 u_\eta }=\frac{4{\rm \pi}}{u_\eta} g(r=2a) \langle w_r (t)\rangle_{r=2a}^-. \end{equation}
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. Reference Ireland, Bragg and Collins2016). Comparing the experimental collision kernel results with those of DNS, we find that our experimental estimates are 
$O(10^3\unicode{x2013}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.
Figure 10. Normalized collision kernel as a function of 
$St$. The case presented in Hammond & Meng (Reference Hammond and Meng2021) is represented as a filled fuchsia square. Linear fit lines are plotted for each particle type. Here DNS values from Ireland et al. (Reference Ireland, Bragg and Collins2016) are shown averaged over 
$Re_{\lambda }$.
The data point shown by a filled fuchsia square is at the identical condition to that previously reported in Hammond & Meng (Reference Hammond and Meng2021). However, they originally made an error, accidentally reporting 
$K_C^*(2a)/u_\eta$, which was subsequently corrected in a corrigendum (Hammond & Meng Reference Hammond and Meng2023). Here we show the corrected 
$K_C^*(2a)=3.6\times 10^4$ for the 
$St=0.74 \text{ and } a=14.25\ \mathrm {\mu }\mathrm {m}$ particles in the current study as well as in the corrigendum (Hammond & Meng Reference Hammond and Meng2023).
5. Discussion
In Hammond & Meng (Reference Hammond and Meng2021) and Bragg et al. (Reference Bragg, Hammond, Dhariwal and Meng2022), the authors argue that the surprising extreme clustering found in their experiments must be a result of particle–particle interactions. Bragg et al. (Reference Bragg, Hammond, Dhariwal and Meng2022) started by hypothesizing that the particle–particle interactions responsible for extreme clustering were HIs, as did Yavuz et al. (Reference Yavuz, Kunnen, Van Heijst and Clercx2018), although the latter had several errors in their theory. Once the theory was corrected in Bragg et al. (Reference Bragg, Hammond, Dhariwal and Meng2022), 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 (Reference Hammond and Meng2021) and Bragg et al. (Reference Bragg, Hammond, Dhariwal and Meng2022) found that their measured RDF did in fact behave as 
$g(r)-1\propto r^{-6}$ as predicted by HI theory in the limit as 
$St\rightarrow 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. (Reference Bragg, Hammond, Dhariwal and Meng2022) 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 (Reference Hammond and Meng2021), and Yavuz et al. (Reference Yavuz, Kunnen, Van Heijst and Clercx2018) 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 (Reference Hammond and Meng2021) and Bragg et al. (Reference Bragg, Hammond, Dhariwal and Meng2022) by sweeping a broad range of 
$St$ values.
5.1. 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_\lambda$ 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.
5.2. 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. (Reference Bragg, Hammond, Dhariwal and Meng2022). This equation identifies four terms contributing to the RDF, associated with their 
$\mu _i$ coefficients. While this equation was derived using the far-field (
$r/a\gg 1$) asymptotic forms of the mobility coefficients describing the HI, these are known to be valid down to 
$r/a \approx 2.05$ and are very accurate for 
$r/a \geq 3$ (Brunk et al. Reference Brunk, Koch and Lion1997). 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 \geq 3$, which includes the region where the RDF exhibits extreme clustering. The 
$\mu _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 
$\mu _2$ and 
$\mu _3$ terms represent the HI contributions that are dependent on particle inertia, i.e. the 
$St$ effect on HI. The 
$\mu _4$ term describes the clustering due to particle inertia interacting with turbulence, which happens even in the absence of HI (Chun et al. Reference Chun, Koch, Rani, Ahluwalia and Collins2005). Since both 
$\mu _2$ and 
$\mu _3$ are negative (Bragg et al. Reference Bragg, Hammond, Dhariwal and Meng2022), the theory indicates that particle inertia should reduce the HI-induced mechanism for clustering (although as discussed in Bragg et al. (Reference Bragg, Hammond, Dhariwal and Meng2022), the overall effect of increasing 
$St$ is to increase the RDF due to the contribution of the non-HI term 
$St^2\mu _4$).
In the limit as 
$St\to 0$, the HI-induced clustering mechanism in the theory of Bragg et al. (Reference Bragg, Hammond, Dhariwal and Meng2022) is the same as that in Brunk et al. (Reference Brunk, Koch and Lion1997). 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. Reference Brunk, Koch and Lion1997; Bragg et al. Reference Bragg, Hammond, Dhariwal and Meng2022). This mechanism leads to the contribution 
$g(r)-1 \sim (r/a)^{-6}$ in (1.1). The theory of Bragg et al. (Reference Bragg, Hammond, Dhariwal and Meng2022) 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 Reference Lambert, Weynans and Bergmann2018; Ababaei et al. Reference Ababaei, Rosa, Pozorski and Wang2021). These near-field effects were included in Brunk et al. (Reference Brunk, Koch and Lion1997), but the effects of particle inertia were not. In the limit as 
$St\rightarrow 0$, (1.1) becomes
\begin{equation} g(r)\sim \exp \left(\mu _1 \left(\frac{r}{a}\right)^{{-}6}\right), \end{equation}
and since the theory assumes that 
$r/a\gg 1$, a Taylor expansion of the exponential function leads to
\begin{equation} g(r)\sim 1+\mu _1 \left(\frac{r}{a}\right)^{{-}6}. \end{equation}
Applying a Taylor expansion to (1.1) for 
$r/a\gg 1$, we obtain the following for finite 
$St$:
\begin{align} g(r)&\sim\left(\frac{r}{a}\right)^{{-}St^2\mu_4}\left(1+\left(\mu_1\left(\frac{r}{a}\right)^{{-}6}+(St\mu_2+St^2\mu_3)\left(\frac{r}{a}\right)^{{-}1}\right)\right) \nonumber\\ &=\left(\frac{r}{a}\right)^{{-}St^2\mu_4}+\mu_1\left(\frac{r}{a}\right)^{-(St^2\mu_4+6)}+(St\mu_2+St^2\mu_3)\left(\frac{r}{a}\right)^{-(St^2\mu_4+1)}. \end{align} If we assume that HI explains the extreme clustering in our experiments when 
$St\ll 1$, then we would expect that the measured 
$g(r)$ primarily behaves according to the 
$\mu _1 (r/a)^{-(St^2 \mu _4+6)}$ term. From Bragg et al. (Reference Bragg, Hammond, Dhariwal and Meng2022), 
$\mu _4$ is positive, which means that the scaling exponent in the extreme clustering range should be 
$k=St^2 \mu _4+6\geq 6$ . However, figure 5 shows that our measured scaling exponent 
$k$ is considerably smaller than 
$6$ when 
$St\ll 1$, reaching values as low as around 
$4.5$. (Note that even though Bragg et al. (Reference Bragg, Hammond, Dhariwal and Meng2022) 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 \lesssim 6.6,$ with the single case presented in Hammond & Meng (Reference Hammond and Meng2021) happening to be 
$k=6$.)
As previously discussed in Bragg et al. (Reference Bragg, Hammond, Dhariwal and Meng2022), 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\ \mathrm {\mu } {\rm m}$ to 
$41.5\ \mathrm {\mu } {\rm m}$, at least an order of magnitude smaller than the Kolmogorov lengths (
$101\unicode{x2013}179\ \mathrm {\mu } {\rm m}$). For (iii), all our experimental conditions have 
$Re_p = a^2/\tau _n \nu \ll 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 Reference Daitche2015). Finally, for (v), as shown in § 5 of Bragg et al. (Reference Bragg, Hammond, Dhariwal and Meng2022), 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. (Reference Bragg, Hammond, Dhariwal and Meng2022) 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. (Reference Bragg, Hammond, Dhariwal and Meng2022) that ‘new, yet-to-be-identified physical mechanisms are at play, requiring further investigation and new theories’.
5.3. 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\rightarrow 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 (Reference Bragg and Collins2014a,Reference Bragg and Collinsb) 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 
$-\boldsymbol {\nabla }_r S_2$ can be understood in terms of a path-history, symmetry-breaking effect (Bragg & Collins Reference Bragg and Collins2014a; Bragg et al. Reference Bragg, Ireland and Collins2015a,Reference Bragg, Ireland and Collinsb). However, this physical interpretation of the effect captured by 
$-\boldsymbol {\nabla }_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. 
$-\boldsymbol {\nabla }_r S_2>0$), highlighting that the physical processes governing 
$-\boldsymbol {\nabla }_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.
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. Reference Dou, Ireland, Bragg, Liang, Collins and Meng2018) and DNS (Ireland et al. Reference Ireland, Bragg and Collins2016). 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. Reference Dou, Ireland, Bragg, Liang, Collins and Meng2018), resolving data down to separations of 
$r/\eta = 1$. There, higher MIRV than DNS predictions was also observed in the separations of 
$r/\eta \leq 10$. Dou et al. (Reference Dou, Ireland, Bragg, Liang, Collins and Meng2018) 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\leq r/\eta \leq 20$ as observed in Dou et al. (Reference Dou, Ireland, Bragg, Liang, Collins and Meng2018), 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. (Reference Dou, Ireland, Bragg, Liang, Collins and Meng2018), 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. (Reference Dou, Ireland, Bragg, Liang, Collins and Meng2018) 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\leq r/\eta \leq 20$ (where the data in Dou et al. (Reference Dou, Ireland, Bragg, Liang, Collins and Meng2018) reached their resolution limit) but continued to trend upward in the sub-Kolmogorov regime before reaching a maximum between 
$0.5\leq r/\eta \leq 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. (Reference Dou, Ireland, Bragg, Liang, Collins and Meng2018) 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.
6. 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 \leq St \leq 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. (Reference Bragg, Hammond, Dhariwal and Meng2022)'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.
Acknowledgements
We thank the National Science Foundation for support through the Major Research Instrumentation Program (Award 1828544, Program Manager H. Chelliah).
Funding
This research received no specific grant from any funding agency in the commercial or not-for-profit sectors.
Declaration of interests
The authors report no conflict of interest.
Appendix A
A.1. Sample size and statistical convergence
Table 2 includes the number of realizations and average particles per frame for each realization. For analysis of each case, the particle-pair separations were binned into 100 bins, the edges of which were based on appropriate scaling relative to the Kolmogorov length of the case. To ensure that the data are statistically significant, we aim for the RDF statistics in each bin to converge with minimal relative standard error (RSE). Additionally, we report the maximum RSE for each case and the bin (in terms of 
$r/a$) for which the maximum RSE occurred in table 2.
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 
$\delta r_{in}$, normalized by 
$a$.
A.2. Interpolation uncertainty
Our interpolation technique gathers four position measurements of a given particle and then interpolates a position, velocity and acceleration at the track midpoint, between pulses two and three (see figure 1). The accuracy of this interpolation is limited by the 
$\Delta t_2$ set for the experiments. Despite this time being very small, it is still finite and thus can cause uncertainty in the interpolation. To quantify this uncertainty, we calculate the root mean square radial distance travelled by particles between the second and third pulse at each bin as 
$\delta r_{in}=\Delta t_2\sqrt {\langle {w_r}^2 (r)\rangle }$ where 
$\langle {w_r}^2 (r)\rangle$ is the variance of the particle-pair radial RV p.d.f. This maximum value for each case is also included in , and figure 11 shows the 
$50\ {\rm Hz}$ test for each particle type with blue horizontal error bars indicating the confidence interval for 
$r\pm \delta r_{in}$ for the RDF, MIRV and RA. The interpolation uncertainty is negligible at larger 
$r$ but widens as 
$r$ decreases. It is clear that the RDF, MIRV and RA 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.
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_{\lambda } = 324$ are shown for all cases except 
$a = 8.75\ \mathrm {\mu }\mathrm {m}$, for which 
$Re_{\lambda } = 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.
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. Reference Novara, Schanz, Geisler, Gesemann, Voss and Schröder2019), which, based on our camera's pixel size of 
$21\ \mathrm {\mu }\mathrm {m}$, corresponds to a position uncertainty of 
$\delta x=3.2\ \mathrm {\mu }\mathrm {m}$ and an uncertainty of 
$\delta r \approx 4.5\ \mathrm {\mu }\mathrm {m}$ estimated via propagation (Moffat Reference Moffat1988). 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, 
$\epsilon$. According to Novara et al. (Reference Novara, Schanz, Geisler, Gesemann, Voss and Schröder2019), it is the most consequential parameter affecting the output. We used 
$\epsilon =1.5$ voxel (where 
$1~\mathrm {voxel}\approx 21~\mathrm {\mu }\mathrm {m}$), since this value produced the most tracks when compared with surrounding values. To quantify the uncertainty in the RDF and RV calculations, 
$\epsilon$ was varied by 
$\pm 10\,\%$ for all of the 
$50\ {\rm Hz}$ cases (
$40\ {\rm Hz}$ for 
$a=8.75\ \mathrm {\mu }\mathrm {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\ \mathrm {\mu }\mathrm {m}$ and 
$14.25\ \mathrm {\mu }\mathrm {m}$ diameter cases but seem to cast some doubt in the 
$8.75\ \mathrm {\mu }\mathrm {m}$ and 
$3.75\ \mathrm {\mu }\mathrm {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.
Open access
