3.1.1. Clustering fraction and length scale

The fraction of particles forming clusters, denoted as $\chi _{cl}$ , is a key metric for understanding how the tendency to aggregating evolves across different flow regimes. This is defined as the number of particles contained in clusters, normalised by the total particle count. Figure 4(a) depicts $\chi _{cl}$ as a function of $Ca$ at various $\phi$ . As already discussed in Shin & Coletti (Reference Shin and Coletti2024), $\chi _{cl}$ increases with decreasing $Ca$ (as capillarity prevails over drag) and increasing $\phi$ (as lubrication prevails over drag). Moreover, for increasing $Ca$ at a given $\phi$ , $\chi _{cl}$ approaches the level obtained in the synthetic data by RSA (dotted line). This implies that, when the fluid turbulence is sufficiently high, the statistical proximity of non-overlapping particles is largely responsible for the observed trend.

Figure 4. (a) Clustered particle fraction ( $\chi _{cl}$ ) as a function of areal fraction ( $\phi$ ) across different capillary numbers ( $Ca$ ). (b) Normalised cluster size ( $\langle 2R_g \rangle _w/{d}_p$ ) as a function of $\phi$ at various $Ca$ levels. In each case, the dotted line represents the clustering fractions from random configurations for the respective $\phi$ values, established through RSA up to $\phi = 0.54$ .

Figure 4(b) shows how the mean cluster length scale depends on $Ca$ and $\phi$ . To avoid biases induced by the large number of small clusters, we define the weighted-average radius as

(3.1) \begin{equation} \langle R_{g} \rangle _{w} \equiv \frac {\sum _{i} n_{p,i} R_{g,i}}{\sum _{i}n_{p,i}}, \end{equation}

where $R_{g,i}$ is the radius of gyration of the $i$ th cluster and $n_{p,i}$ is the number of particles forming it. At low $Ca$ , particles are heavily influenced by strong attractive capillary forces exceeding the disruptive effect of drag, leading to extended clusters. These can grow larger than a typical eddy size $L_F$ , comprising hundreds of particles. With increasing $Ca$ , these objects have higher chance of being exposed to extensional strain exceeding the capillarity, leading to more frequent break-up and thus smaller aggregates. Again, with increasing $Ca$ the radius obtained by RSA is approached. Overall, both $\chi _{cl}$ and $\langle R_{g} \rangle _{w}$ indicate that, when the turbulent forcing is intense and $Ca \gg 1$ , the clustering is increasingly governed by the stochastic adjacency of particles, rather than inter-particle interactions.

Figure 5. Probability density functions (PDFs) of the cluster size, $R_g$ , normalised by the exponential decay length, ( $l_{exp }$ ): (a) RSA-generated random configurations at $\phi = 0.27$ , 0.44 and 0.53. The corresponding values of $l_{exp }/{d}_p$ are 0.20, 0.33 and 0.87, respectively. (b) The PDFs at $\phi = 0.27$ for various experimental capillary numbers: $Ca= 0.16$ , 0.23, 0.29, 0.87, 1.22, 1.51, 1.82, 2.10, 2.49 and 2.69, with the corresponding values of $l_{exp }/{d}_p= 5.66$ , 4.16, 3.76, 0.41, 0.37, 0.33, 0.29, 0.32, 0.25 and 0.21, respectively. (c) The PDFs at $\phi =0.53$ for $Ca = 0.59$ , 0.87, 1.22, 1.51, 1.82, 2.10, 2.49 and 2.69 with the corresponding values of $l_{exp }= 1.17$ , 1.00, 1.00, 0.95, 0.88, 0.83, 0.79 and 0.74. The black dashed line in each panel represents an exponential decay, $\exp (-R_g/l_{exp })$ .

The PDFs of $R_{g}$ offer further insight into how the clustering dynamics changes across regimes. Figure 5(a) shows that clusters formed by RSA exhibit an exponential distribution, independently of $\phi$ . An exponential fit of the probability $p(R_g) \sim \exp (-R_g / l_{exp })$ , applied to the RSA data as well as to experimental data that exhibit exponential decay, defines a characteristic length $l_{exp }$ used to normalise the distributions in the different scenarios. Moderately dense systems (figure 5 b) exhibit a transition with varying $Ca$ : in the capillary-driven clustering regime ( $Ca \lt 1$ ) the cluster radius PDFs possess long tails, with significant probability of large aggregates spanning the entire domain; while in the drag-driven break-up regime ( $Ca \gt 1, \phi \lt 0.40$ ) drag disrupts the large clusters and the distribution gradually approaches the RSA behaviour. This observation supports the view that the clustering dynamics in moderately dilute systems is predominantly influenced by the random proximity of particles.

The scenario is different for denser systems (figure 5 c). Here, both in the capillarity-driven regime ( $Ca \lt 1$ ) and in the lubrication-driven clustering regime ( $Ca \gt 1, \phi \gt 0.40$ ), the decay of the PDFs is much slower compared with RSA. This implies that the formation of very large clusters cannot be explained by excluded volume effect at high concentration, and are instead sustained by lubrication forces. At even higher concentrations, few large clusters percolate through the whole system. This scenario is discussed in the next section.

3.1.2. Cluster percolation and self-similarity

Dense particle systems are known to exhibit sharp transitions to a percolating state when their concentration increases above a threshold. In such a state, exceptionally large clusters are formed that span the entire domain. The threshold is system-dependent: for example, Shen, O’Hern & Shattuck (Reference Shen, O’Hern and Shattuck2012) found a critical area fraction $\phi ^{*}=0.549$ by simulating periodic 2-D layers of frictionless particles under compression, while the experiments of Alicke, Stricker & Vermant (Reference Alicke, Stricker and Vermant2023) found $\phi ^{*}=0.3-0.36$ for particle layers under shear and compression.

The present system also exhibits a transition to a percolating state. This is illustrated in figures 6(a) and 6(b), where instantaneous realisations are shown for two representative cases at $Ca = 0.33$ . While both have comparable concentrations $\phi = 0.44$ and $0.53$ , only the latter displays percolating clusters. The transition is quantified in figure 6(c), plotting the percolation fraction $\chi _p$ (number of imaging frames featuring at least one cluster percolating through the whole field of view, divided by the total number of acquired frames) as a function of $\phi$ , for a representative value of $Ca=0.33$ . While we do not vary $\phi$ finely enough to identify a precise threshold $\phi ^{*}$ , the transition appears relatively sharp. Visual inspection also confirms that the formation of a percolation cluster is hindered by increasing $Ca$ , which favours more homogeneous patterns (figure 6 d).

The trend across the parameter space is summarised in figure 6(e). For $Ca \ll 1$ , we estimate $\phi ^{*} \sim 0.4 - 0.5$ , comparable to classic configurations of compressed/sheared particle layers (Alicke et al. Reference Alicke, Stricker and Vermant2023). On the other hand, for $Ca \geqslant 1$ , the threshold is significantly higher, $\phi ^{*} \sim 0.55 - 0.6$ . This is consistent with the notion that percolation naturally emerges in self-similar systems, where the influence of a preferential scale characterising the pattern is weak or absent (Aharony & Stauffer Reference Aharony and Stauffer2003). At large $Ca$ , the flow increasingly imposes the forcing length scale $L_F$ that disrupts system-spanning clusters, and percolation can only be attained at the largest concentrations. In particular, as we showed in Shin & Coletti (Reference Shin and Coletti2024), for $Ca\gt 1$ , the typical cluster size exceeds $L_{F}$ only for areal fractions above $\phi \sim 0.55$ (which we now identify as the approximate percolation threshold in this regime).

The dominant role of the forcing scale, as opposed to the larger scales where energy accumulates in Q2-D turbulence, was already stressed in studies of Lagrangian dispersion (Xia et al. Reference Xia, Francois, Punzmann and Shats2013, Reference Xia, Francois, Punzmann and Shats2019). For both particle clustering and dispersion, such a role is likely rooted in the persistent nature of eddies of size $L_F$ (see discussion in Xia et al. Reference Xia, Francois, Punzmann and Shats2013). On the other hand, when inter-particle attraction becomes dominant, cohesive clusters much larger than $L_F$ appear. One may hypothesise that those might be related to commensurately large-scale fluid motions, either due to the inverse cascade or engendered by the influence of the clusters themselves on the flow. To investigate this possibility, we define the Eulerian velocity autocorrelation $\Gamma _E(r)=\langle \boldsymbol{u}(\boldsymbol{x}_0,t) \cdot \boldsymbol{u}(\boldsymbol{x}_0 + \boldsymbol{r},t) \rangle _{\boldsymbol{x}_0,t} / \langle \boldsymbol{u}(\boldsymbol{x}_0,t) \cdot \boldsymbol{u}(\boldsymbol{x}_0,t) \rangle _{\boldsymbol{x}_0,t}$ and compare it for both particles and tracers. This is shown in figure 7 for the representative case $Ca=2.69$ and $\phi =0.62$ . It highlights how, even when particles form large percolating clusters, the correlation scale of their motion matches that of the fluid flow (approximately $0.5L_F$ , see Shin et al. Reference Shin, Coletti and Conlin2023). We deduce that the emergence of system-spanning clusters is not related to fluid motions at scales larger than $L_F$ .

Figure 6. Snapshots with identified clusters at the same $Ca = 0.33$ but different areal fractions: (a) $\phi = 0.44$ , and (b) $\phi = 0.53$ . Particles within the same cluster are indicated with the same colour. (c) Fraction of frames $\chi _p$ where a system-spanning cluster is present as a function of $\phi$ at fixed $Ca = 0.33$ . (d) Examples of cluster morphology for systems with $\phi \approx 0.44$ at different capillary numbers: $Ca = 0.29$ , 0.59 and 1.82 (from left to right). At this $\phi$ , $Ca \ll 1$ entails the formation of a percolating structure, while higher $Ca$ results in non-connected structures. (e) Value of $\chi _p$ in the $Ca-\phi$ space. Panels (a) and (b) correspond to the points inside the red solid rectangle, and panel (d) corresponds to cases inside the red dotted rectangle.

Figure 7. Examples of the Eulerian velocity autocorrelation for the DL configuration. The dashed line represents tracers ( $Re=1047$ ), while the red line corresponds to a dense suspension of particles at the same level of forcing ( $Ca=2.69$ and $\phi =0.62$ ), where they form a large, system-spanning cluster.

In order to further characterise the present system within the percolation framework, we analyse the cluster size distributions. Unlike in sub-§ 3.1.1 where we focused on a geometric scale, here, we focus on the size, intended as the number of particles $n_p$ contained in each cluster. In particular, we evaluate the cluster number density $f(n_p)$ , namely the number density of clusters of size $n_p$ per unit area. In 2-D self-similar particle systems (e.g. patterns of droplets condensing on surfaces (Stricker et al. Reference Stricker, Grillo, Marquez, Panzarasa, Smith-Mannschott and Vollmer2022), we expect $f(n_p) \sim (n_p/n_{p,max})^{-\xi }(n_{p,max} A_p)^{-\theta }$ , where $n_{p,max}$ is the size of the largest cluster in the distribution, $\xi$ is the so-called polydispersity exponent characterising the scaling range, $A_p=\pi {d}_p^{2}/4$ is the projected area of a particle and $\theta$ is a trivial exponent depending on the dimensionality of the system, which has here the value of 1 to guarantee dimensional consistency.

Figure 8(a) reports $f(n_p)$ for cases of comparable areal fraction, $\phi \sim 0.44$ , displaying self-similar scaling ranges with a polydispersity exponent that varies with $Ca$ (considering comparable or higher concentrations, not shown for brevity, leads to analogous conclusions). To extract the polydispersity exponents, we rescale $f(n_p)$ in figure 8 to have $f(n_p)(n_{p,max}A_p)^{\theta }$ , which is expected to scale as $(n_p/n_{p,max})^{-\xi }$ .

Figure 8. (a) Cluster size distributions for a particle area fraction $\phi \approx 0.44$ at different capillary numbers $Ca$ . Here, $n_p$ is the number of particles within a cluster, characterising the cluster size and $f(n_p)$ is the areal number density of clusters of size $n_p$ . Dashed lines denote experiments in the SL configuration, while solid lines are from the DL2 set-up. (b) Normalised form of the same data, where $n_{p,max}$ is the size of the largest cluster, $A_p = \pi {d}_p^{2}/4$ is the projected area of a particle and $\theta = 1$ is a trivial exponent from dimensional analysis. Two limiting values of the polydispersity exponents, $\xi =2.05$ (dash-dot) and $\xi =5/2$ (dotted), which correspond to the predicted Fisher’s exponent (Fisher Reference Fisher1967), are indicated.

For $Ca \lt 1$ , the scaling range is broad and $\xi$ is somewhat smaller than but comparable to the value of 2.05 characterising 2-D percolation, as observed in randomly generated clusters on regular lattices (Stauffer Reference Stauffer1979). The discrepancy may be attributed to finite size effects.

This indicates that, in the capillary-driven clustering regime, the present system shares similarities (at least topologically) with classic percolating particle systems. For large capillary numbers $Ca \gt 1$ , on the other hand, the scaling range is narrower as the forcing hinders the formation of self-similar structures. Over approximately a decade, however, the size distribution is consistent with a polydispersity exponent $\xi = 5/2$ . This corresponds to the well-known prediction derived by Fisher (Reference Fisher1967) for a coagulation process under the assumption that each agglomerate has equal probability of interacting with any other. Such a condition is typically not satisfied in standard percolation processes, leading to shallower exponents (Hoshent et al. Reference Hoshent, Stauffer, Bishop, Harrison and Quinn1979; Stauffer Reference Stauffer1979). We therefore hypothesise that, in the present configuration, the forcing by the underlying turbulence contributes to randomising the system, effectively restoring the validity of Fisher’s assumption. We remark that this is the first observation of a percolation transition in turbulent particle suspensions, and its detailed properties deserve further investigation.

Figure 9. Fractal dimension $d_{frac}$ obtained via the box-counting method. (a) An example of normalised box count as a function of box side length for a cluster (inset), (b) $d_{frac}$ mapped in the $Ca-\phi$ space and (c) $d_{frac}$ as a function of the mean cluster size. The red dotted line indicates the mean $d_{frac}$ from varying $\phi$ values obtained via RSA.

3.1.3. Clusters’ internal structure

The different mechanisms at play in the $Ca-\phi$ space can affect not only the clusters’ size but also their internal structure. Even though particle patterns may not be strictly self-similar, an operational definition of fractal dimension $d_{frac}$ can be useful to quantify the compactness and complexity of their geometry. To estimate it, we select in each frame the cluster containing the largest number of particles, and we reconstruct its binary image based on the centroids and diameter of its particles. We then apply to it the box-counting method (Monchaux, Bourgoin & Cartellier Reference Monchaux, Bourgoin and Cartellier2012; Baker et al. Reference Baker, Frankel, Mani and Coletti2017), with $L$ the box size (figure 9 a). The data are normalised by the maximum count $N_{b,max}$ obtained from the smallest box size ( $L=O({d}_p)$ ) yielding a plot of $N_b / N_{b,max}$ as a function of $L$ , where $N_b$ is the number of boxes that cover the binarised cluster. The normalised curves from all analysed clusters are averaged to extract $d_{frac}$ by linear fitting in logarithmic scale.

Figure 9(b) maps the fractal dimension in the $Ca-\phi$ space. Similarly to $\chi _{cl}$ and $\langle R_g \rangle _w$ previously discussed, $d_{frac}$ drops sharply in moderately dense suspensions ( $\phi \lt 0.40$ ) within a narrow range of $Ca\lt 1$ . Further increasing $Ca$ leads to $d_{frac}$ flattening around 1.5, close to the baseline value of $1.52 \pm 0.2$ from RSA which exhibits negligible $\phi$ dependence. When $Ca \gt 1$ , the contour lines become largely horizontal, signalling a major role of $\phi$ . The lubrication-driven clusters are indeed compact, the particles filling the open spaces in their internal structure.

We now consider how $d_{frac}$ varies with $\langle R_{g} \rangle _w$ (figure 9 c). While small clusters (up to $\langle R_g \rangle _w / {d}_p = O(10)$ ) exhibit values of $d_{frac}$ similar to those generated by RSA, for larger ones $d_{frac}$ increases with size, tending to (although never reaching) the value of 2 characteristic of round (non-fractal) objects. For large $\langle R_g \rangle _w / {d}_p = O(10^{2})$ , $d_{frac}$ is significantly lower for $Ca \lt 1$ . This indicates a convoluted structure, which can be the result of hit-and-stick collisions, frequent when the capillary forces are too strong to be overcome by drag. In contrast, clusters at $Ca\gt 1$ have a more dynamic composition, with particles continuously joining and leaving the aggregate due to relatively weaker cohesive forces. Thus, we can draw an analogy between the cluster formation at low $Ca$ and the diffusion-limited cluster aggregation (DLCA) which is known to result in more open (less compact) clusters (Meakin Reference Meakin1983; Lazzari et al. Reference Lazzari, Nicoud, Jaquet, Lattuada and Morbidelli2016). Notably, unlike DLCA of colloids driven by Brownian diffusion, our system employs non-Brownian particles driven by the underlying turbulent flow. In systems exhibiting aggregation of particles in turbulence, the fractal dimension of the aggregates has been evaluated both experimentally (e.g. Maggi, Mietta & Winterwerp Reference Maggi, Mietta and Winterwerp2007) and numerically (e.g. Zhao et al. Reference Zhao, Pomes, Vowinckel, Hsu, Bai and Meiburg2021). In particular, Zhao et al. (Reference Zhao, Pomes, Vowinckel, Hsu, Bai and Meiburg2021) provided empirical scaling for $d_{frac}$ as a function of the key non-dimensional parameters. The present results qualitatively agree with their model, e.g. concerning the increase of $d_{frac}$ with larger cluster size. A quantitative comparison, however, is hindered by the fact that the present configuration is quasi-two-dimensional.

A natural metric to assess the degree of ordering/packing inside clusters is the hexatic order parameter $\psi _{6}$ , which measures the local ordering based on the orientation of the particle’s six nearest neighbours. For the particle $j$ , it is defined as

(3.2) \begin{equation} \psi _{6,j} = \frac {1}{6}\left |\sum _{k=1}^{6}\exp (i6\theta _{jk})\right |, \end{equation}

where $k$ is an index for the six nearest neighbouring particles, $\theta _{jk}$ is the angle between the reference direction and the line connecting the particle $j$ and its neighbouring particle $k$ .

We assess $\psi _6$ for all individual clustered particles, ranging between 0 (dislocations and disclinations) and 1 (hexagonal crystal). Figure 10(a) compares the PDFs of $\psi _6$ for the same $\phi =0.44$ and two different $Ca$ . It is clear that lower $Ca$ yields higher probabilities of large $\psi _6$ values. Increasing $Ca$ pushes the PDF towards lower $\psi _6$ values, approaching the distribution obtained by RSA at the corresponding area fraction.

We further consider the average hexatic order parameter over all clustered particles, $\langle \psi _6 \rangle$ , in the $Ca-\phi$ space (figure 10 b). This displays a similar trend to previously discussed metrics: it diminishes as the system transitions from the capillary-driven clustering to drag-driven regime, and increases transitioning to lubrication-driven clustering. Plotting $\langle \psi _6 \rangle$ as a function of $\phi$ (figure 10 c) indicates how the high- $Ca$ cases follow closely the values obtained by RSA. This indicates again that clusters formed under strong turbulence forcing tend to resemble those formed by a random process.

Figure 10. Hexatic order parameter ( $\psi _6$ ) of the clustered particles. (a) The PDFs of $\psi _6$ at $\phi = 0.44$ but different $Ca$ values: 0.09 (red) and 2.69 (blue). The black dashed line represents values from RSA for comparison. (b) Value of $\langle \psi _6 \rangle$ mapped in the $Ca-\phi$ space. (c) Value of $\langle \psi _6 \rangle$ as a function of $\phi$ across varying $Ca$ . The black dashed line represents values obtained from RSA.

3.1.4. Cluster lifetime

The observation that clusters formed at low $Ca$ are more compact than those at high $Ca$ suggests that the former may also be longer lived. Verifying this hypothesis, and in general understanding the temporal dynamics of particle clusters within turbulent flows, requires estimating a cluster’s lifetime. The challenge, however, is represented by the dynamic nature of such entities where particles continuously join and escape.

Following Liu et al. (Reference Liu, Shen, Zamansky and Coletti2020), we regard clusters in successive time steps as the same entity if they share the majority of their particles. This method mirrors the philosophical inquiry of the Ship of Theseus (see, e.g. Rea Reference Rea1995). For example, consider two clusters A and B in two successive time steps. The fraction of forward-in-time connection is defined as the number of particles shared by A and B divided by the total number of particles in A; while the fraction of backward-in-time connection is the number of shared particles divided by the total number of particles in B. Clusters A and B are considered continuous manifestations of the same cluster if their fraction of forward- and backward-in-time connections both exceed 0.5, and their lifetime is the time during which such a condition is continuously verified. We then estimate the weighted-averaged cluster lifetime $\langle \tau _{cl} \rangle _{w}$ as

(3.3) \begin{equation} \langle \tau _{cl} \rangle _{w} \equiv \frac {\sum _{i} n_{p,i} \tau _{cl,i}}{\sum _{i}n_{p,i}}, \end{equation}

where $\tau _{cl,i}$ is the lifetime of the $i$ th cluster.

To compare among different cases, the cluster lifetime is to be normalised by a time scale that reflects the particle–fluid dynamics. A natural choice is the eddy turnover time of the underlying flow, $L_F/u_{rms,f}$ . This is verified by analysing the dispersion and stretching of the group of particles that initially belong to a cluster. We use singular value decomposition (Baker et al. Reference Baker, Frankel, Mani and Coletti2017; Petersen, Baker & Coletti Reference Petersen, Baker and Coletti2019), which identifies the first principal axis along the main direction of spread, the second axis being orthogonal to it. The corresponding singular values, $s_1$ and $s_2$ , are proportional to the spread of particle centroids along these axes. Figure 11(a) shows a representative example of temporal evolution of the normalised length scale $\langle R_g (t) / R_{g,0} \rangle$ and anisotropy $[s_2 (t) / s_1 (t)]/[s_{2,0}/s_{1,0}]$ , where the subscript 0 indicates the time when the cluster is first identified. The trend confirms that sizeable changes to the spatial organisation occur within a time $O(L_F / u_{rms,f})$ .

Figure 11. (a) Temporal evolution of the cluster size $\langle R_g(t)/R_{g,0} \rangle$ and cluster anisotropy $\langle [s_2(t)/s_1(t)]/[s_{2,0}/s_{1,0}] \rangle$ for an example case with $Ca = 2.49$ and $\phi = 0.44$ . (b) Contour map of the weighted average cluster lifetime $\langle \tau _{{cl}} \rangle _w$ normalised by the flow turnover time $L_F/u_{{rms},f}$ across the $Ca-\phi$ space.

For a given eddy turnover time, however, the cluster lifetime can be significantly larger or smaller depending on the balance of the forces at play. Figure 11(b) displays a contour map of the normalised cluster lifetime in the $Ca-\phi$ space. As expected, clusters have longer lifetimes for lower $Ca$ , where cohesive capillary forces dominate. As $Ca$ increases, particularly in the moderately dense regime ( $\phi \lt 0.4$ ), clusters can live for just a fraction of the eddy turnover time before getting disrupted by the high strain rate of the intensely turbulent flow. On the other hand, in the denser regimes ( $\phi \gt 0.4$ ) the cluster lifetimes can far exceed $L_F/u_{rms,f}$ irrespective of $Ca$ . This is due to large clusters being brought about by excluded volume effects and kept together by lubrication. The particle dynamics within these large clusters, however, does vary substantially depending on $Ca$ , as we show in the next section.

3.2.1. Single-particle dispersion

In Shin & Coletti (Reference Shin and Coletti2024), we demonstrated how Eulerian properties of the particle motion, in particular their fluctuating kinetic energy, vary within the $Ca-\phi$ space. Here we focus on Lagrangian transport, which we first characterise in terms of the single-particle dispersion. It is worth remarking that the peculiar properties of 2-D turbulent flows impact the Lagrangian dispersion. In particular, as shown by Xia et al. (Reference Xia, Francois, Punzmann and Shats2013) and confirmed by Shin et al. (Reference Shin, Coletti and Conlin2023), the length scale that determines the dispersion process is the forcing scale $L_F$ . This is at odds with the behaviour of 3-D turbulence, where dispersion is driven by the range of scales in which most energy reside. This characteristic was confirmed also in the transport of finite-size particles in Q2-D turbulence (Xia et al. Reference Xia, Francois, Punzmann and Shats2019).

In the following, we consider the displacement $\boldsymbol{r}_i$ of the $i$ th particle from an initial position at time $t_0$ , and calculate the mean squared displacement (MSD) by averaging over all particles and initial times

(3.4) \begin{equation} \big\langle r^{2}(t) \big\rangle = \big\langle |\boldsymbol{x}_{i}(t_{0}+t) - \boldsymbol{x}_{i}(t_{0})|^{2} \big\rangle _{i,t_{0}}, \end{equation}

where $\boldsymbol{x}_i$ represents the centroid position of the $i$ th particle. As the Lagrangian velocity autocorrelation $\Gamma _{L}(t) = \langle {\boldsymbol{u}_{i}(t_{0})\cdot \boldsymbol{u}_{i}(t_{0}+t)}\rangle _{i,t_{0}} / \langle {\boldsymbol{u}_{i}(t_{0})\cdot \boldsymbol{u}_{i}(t_{0})}\rangle _{i,t_{0}}$ approximates an exponential decay of characteristic time $t_L$ , one expects the classic scaling $\langle r^{2}(t) \rangle \approx u_{rms}^{2} t^{2}$ for $t \ll t_L$ , and $\langle r^{2}(t) \rangle \approx 2u_{rms}^{2} t_{L} t$ for $t \gg t_L$ , where $u_{rms}$ is the root-mean-squared velocity. The approximately exponential decay of $\Gamma _L$ for the turbulent flow was demonstrated in Shin et al. (Reference Shin, Coletti and Conlin2023) over a Lagrangian time scale of the flow that we term $t_{L,f}$ . For the particles, $\Gamma _L$ is shown by figure 12(a) in the example case $Ca = 0.87$ and $\phi = 0.26$ . The characteristic time of the Lagrangian particle motion, $t_{L,p}$ , is taken as the e-fold decay time of $\Gamma _L$ and found to be comparable to the turnover time of the particle motion, $L_F/u_{rms,p}$ . Figure 12(b) plots the MSD for the particles, confirming the expected short-time and long-time behaviours with a ballistic-to-diffusive transition around $t \sim t_{L,p}$ .

Figure 12. An example of (a) the Lagrangian velocity autocorrelation function (VACF) and (b) the MSD of particles for the case $Ca = 0.87$ and $\phi = 0.27$ . The blue dotted line and red dash-dot line represent the short-term ballistic and long-term diffusive predictions, respectively, based on an exponentially decaying VACF with characteristic time scale $t_{L,p}$ . The time is normalised by $t_{L,p}$ , and the MSD is normalised by $L_p^{2} = (u_{rms,p} t_{L,p})^{2}$ .

The diffusive behaviour allows us to evaluate the long-term diffusivity associated with the particle motion as (Taylor Reference Taylor1922; Davidson Reference Davidson2015)

(3.5) \begin{equation} K_p=\frac {1}{2} \frac {d \big\langle r^{2}(t) \big\rangle }{dt} \approx u_{rms,p} L_{p} = u_{rms,p}^{2} t_{L,p}, \end{equation}

where $L_p = u_{rms,p} t_{L,p}$ is a characteristic length scale of the particle motion. This is compared with the turbulent diffusivity of the fluid flow $K_f$ , defining the normalised particle diffusivity

(3.6) \begin{equation} \frac {K_p}{K_f} = \frac {u_{rms,p}^{2} t_{L,p}}{u_{rms,f}^{2} t_{L,f}} = \frac {\langle E_{k,p} \rangle }{\langle E_{k,f} \rangle } \frac {t_{L,p}}{t_{L,f}}, \end{equation}

where $\langle E_{k,p} \rangle$ and $\langle E_{k,f} \rangle$ are the mean kinetic energies of particles and fluid, respectively. Figures 13(a) and 13(b) map in the $Ca-\phi$ space the two factors: the particle-to-fluid ratio of fluctuating energy $\langle E_{k,p} \rangle /\langle E_{k,f} \rangle$ , and the particle-to-fluid ratio of Lagrangian time scales $t_{L,p}/t_{L,f}$ , respectively. Figure 13(a) reflects the observations made by Shin & Coletti (Reference Shin and Coletti2024). In the capillary-driven clustering regime ( $Ca \lt 1$ ) clusters can grow larger than the typical eddy size $L_F$ , thus their motion is hindered as they filter out the energy of the smaller flow scales. In the drag-driven break-up regime ( $Ca \gt 1, \phi \lt 0.4$ ), where the tendency to clustering is reduced, individual particles and small clusters have little inertia and thus possess a fluctuating energy comparable to the tracers. In the lubrication-driven clustering regime ( $Ca \gt 1, \phi \gt 0.4$ ), the particle kinetic energy is largely dissipated by the interaction with adjacent particles. Figure 13(b) shows how, if a particle belongs to a large compact cluster, it exhibits a long-time-correlated motion as it moves collectively with other members of the cluster. This results in $t_{L,p}/t_{L,f}$ exceeding unity for dense concentrations, in particular in the lubrication-driven clustering regime. The particle diffusivity in figure 13(c) is the result of both competing trends: reduced energy and longer time correlation for large/compact clusters, and vice versa. Quantitatively, the trend of the fluctuating energy prevails, and the map of $K_p/K_f$ in the $Ca-\phi$ space closely resembles that of $\langle E_{k,p} \rangle /\langle E_{k,f} \rangle$ . In conclusion, the clustering fostered by capillarity and lubrication reduces the particles’ ability to diffuse through the system.

Figure 13. Contour maps in the $Ca-\phi$ space of (a) the ratio of Lagrangian time scales $t_{L,p} / t_{L,f}$ , (b) the ratio of particle to fluid kinetic energy $\langle E_{k,p} \rangle / \langle E_{k,f} \rangle$ and (c) the normalised particle turbulence diffusivity $K_p / K_f$ .

3.2.2. Velocity structure functions

Pair statistics are a classic tool to reveal how particles move relatively to each other. In order to highlight the dynamics of the aggregates, we restrict our attention to pairs of particles belonging to the same clusters, although this does not change the trend with respect to the unconditioned data analysis. We begin by considering the longitudinal and transverse second-order velocity structure functions, defined respectively as

(3.7) \begin{equation} S_2^L(r) = \left \langle \left ( \delta \boldsymbol{u}_{ij} \cdot \frac {\boldsymbol{r}_{ij}}{r} \right )^{2} \right \rangle , \end{equation}

(3.8) \begin{equation} S_2^T(r) = \left \langle \left ( \delta \boldsymbol{u}_{ij} - \delta \boldsymbol{u}_{ij} \cdot \frac {\boldsymbol{r_{ij}}}{r} \right )^{2} \right \rangle . \end{equation}

Here, $\delta \boldsymbol{u}_{ij} = \boldsymbol{u}_j - \boldsymbol{u}_i$ is the velocity difference between the $i$ th and $j$ th particles whose separation vector is $\boldsymbol{r}_{ij} = \boldsymbol{x}_j - \boldsymbol{x}_i$ , with $|\boldsymbol{r}_{ij}| = r$ . The longitudinal component corresponds to $\delta \boldsymbol{u}_{ij}$ projected along the direction of $\boldsymbol{r}_{ij}$ , which is responsible for the particles’ separation, and the transverse component accounts for the element of $\delta \boldsymbol{u}_{ij}$ perpendicular to $\boldsymbol{r}_{ij}$ , leading to rotation of the separation vector.

Figures 14(a–c) compares structure functions for both particles and tracers for the same level of forcing, normalised by $u_{rms,p}^{2}$ and $u_{rms,f}^{2}$ , respectively. The three identified regimes display distinctly different behaviours. In the capillary-driven clustering regime (figure 14 a), both $S_2^L$ and $S_2^T$ are lower than those of tracers. Particles are tightly bound by capillarity within compact clusters, thus they translate collectively with relatively small longitudinal velocity differences. The transverse component of the relative velocity is inhibited to a degree that depends on the cluster size, but in general much less so than its longitudinal counterpart. In the drag-driven break-up regime (figure 14 b), clusters cannot be sustained in the region of high strain rate, thus they tend to reside in large vortex cores where they spin as quasi-rigid bodies. As those compact clusters spin, the particles in them that are separated by relatively large distances attain a rotational motion (with respect to each other) even faster than fluid parcels at equivalent separation. In those same clusters, on the other hand, inter-particle interactions limit the longitudinal separation rate. In the lubrication-driven clustering regime (figure 14 c), particle adjacency imposed by excluded volume leads to the formation of large clusters spanning multiple eddies. While their lifetime is relatively long (see § 3.1.4), they are continuously deformed by the underlying turbulent flow, and their normalised $S_2^L$ and $S_2^T$ approach those of tracers.

Figure 14. Examples of the second-order velocity structure functions ( $S_2$ ) of particles from different experiments: (a) $Ca = 0.16$ and $\phi = 0.14$ , (b) $Ca = 1.82$ and $\phi = 0.27$ and (c) $Ca = 2.69$ and $\phi = 0.62$ . In each panel, $S_2^L$ (black squares) and $S_2^T$ (red circles) are normalised by $u_{rms,p}^{2}$ . The dashed and dash-dot lines represent $S_2^L$ and $S_2^T$ of tracers at the same level of forcing, respectively, normalised by $u_{rms,f}^{2}$ .

3.2.3. Particle pair dispersion

Finally, we characterise the particle relative dispersion in a Lagrangian framework. To this end, we measure the mean squared separation (MSS) by tracking the evolution of the separation vector $\boldsymbol{r}_{ij}$ between particle pairs as a function of time

(3.9) \begin{equation} MSS(t) = \big\langle |\boldsymbol{r}_{ij}(t_0 + t) - \boldsymbol{r}_{ij}(t_0)|^{2} \big\rangle . \end{equation}

We focus on particle pairs that are initially adjacent, i.e. for which, at the initial time, $r_{ij,0} \sim {d}_p$ .

To define the normalisation time scale, we first resort to the theory of 2-D turbulence. In the enstrophy cascade subrange, classical arguments for the separation law of particle pairs in smooth flows predict an exponential growth of the MSS at short times, with a characteristic time scale $t_{\zeta } = \zeta ^{-1/3}$ , where $\zeta \equiv \nu \langle |\nabla \omega |^{2} \rangle$ is the enstrophy dissipation rate and $\omega$ is the vorticity (Jullien Reference Jullien2003). We estimate $\zeta$ by modelling the flow field as an array of Taylor–Green vortices (Shin et al. Reference Shin, Coletti and Conlin2023)

(3.10) \begin{equation} \zeta /\nu = \big\langle |\nabla \omega |^{2} \big\rangle = \frac {4 \pi ^{4} u_{rms,f}^{2}}{L_{F}^{4}} .\end{equation}

This simplification yields quantitative agreement with the PIV measurements (figure 15 a).

Figure 15. (a) Relationship between the enstrophy dissipation rate divided by the kinematic viscosity ( $\zeta / \nu$ ) and $u_{rms,f}^{2} / L_F^{4}$ . Markers are obtained from PIV analysis of tracer experiments, and the red dotted line represents the estimation based on Taylor–Green vortices with a slope of $4\pi ^{4}$ . (b) Examples of the MSS of tracers in the DL configuration at different levels of forcing, with an initial separation of 1.84 mm. The dashed and dash-dot lines indicate $t^{2}$ and $t^3$ scalings, respectively. (c) Compensated MSS, where the MSS is divided by $(t / t_\zeta )^3$ , highlighting the $t^3$ scaling as a plateau.

Based on this estimate, we approximate the characteristic time scale as $t_{\zeta } \approx (4 \pi ^{4} \nu u_{rms,f}^{2} / L_F^{4})^{-1/3}$ . For short times ( $t \ll t_{\zeta }$ ), the MSS is primarily influenced by the initial velocity difference, $\langle \delta u_{r}^{2}(t) \rangle \approx \langle \delta u_{r_{0}}^{2} \rangle$ . We therefore rescale the MSS as $\langle |\boldsymbol{r}_{ij}(t_0 + t) - \boldsymbol{r}_{ij}(t_0)|^{2} \rangle /(t_{\zeta }^{2} \langle \delta u_{r_{0}}^{2} \rangle )$ , similar to Li et al. (Reference Li, Wang, Qi and Coletti2024) who studied pair-dispersion in free-surface turbulence, with $t_{\zeta }$ replacing their transition time $t_{D}$ between ballistic and diffusive behaviours of the relative velocity $\delta u$ .

Figure 15(b) presents normalised MSS curves from tracers in the DL configuration, with an initial separation $r_{ij,0} \sim {d}_p$ . This is close to the range of initial separation of $2-4\eta$ that Tan & Ni (Reference Tan and Ni2022) proposed as an upper limit for observing the super-diffusive scaling $\langle |\boldsymbol{r}_{ij}(t_0 + t) - \boldsymbol{r}_{ij}(t_0)|^{2} \rangle / \sim t^3$ in single-phase homogeneous turbulence. This scaling was originally proposed by Richardson (Reference Richardson1926) based on a scale-dependent diffusivity, and later justified by Obukhov (Reference Obukhov1941) as consistent with Kolmogorov’s (Reference Kolmogorov1941) theory. According to the latter, the dissipation rate of the turbulent kinetic energy is the only governing parameter in the inertial range, which dimensionally leads to the super-diffusive scaling. Here, the MSS curves collapse well on the chosen normalisation. For $t/t_{\zeta } \ll 1$ , we note that experimental limitations such as resolution constraints affect the precise measurement of very small relative displacements of nearby tracers. This is likely the cause of the weak decline in separation observed at low Reynolds numbers. For $t/t_{\zeta } \gt 1$ , the MSS curves display the super-diffusive scaling, and this persists until the separation grows beyond $L_{F}$ , where the dispersion becomes purely diffusive.

The $t^3$ -scaling range (highlighted by the compensated plots in figure 15 c) is therefore limited by the scales of the system, and possibly by factors such as cross-scale contamination, which have made the super-diffusive regime particularly elusive (Bourgoin et al. Reference Bourgoin, Ouellette, Xu, Berg and Bodenschatz2006; Salazar & Collins Reference Salazar and Collins2009; Bitane, Homann & Bec Reference Bitane, Homann and Bec2012; Elsinga, Ishihara & Hunt Reference Elsinga, Ishihara and Hunt2022; Tan & Ni Reference Tan and Ni2022). The limited range over which the super-diffusive scaling is observed, and the fact that such an observation typically requires small initial separations (below the inertial range prescribed by Richardson–Obukhov theory) has created significant debate on the nature of the process in both 2-D and 3-D turbulence (Kellay & Goldburg Reference Kellay and Goldburg2002; Elsinga et al. Reference Elsinga, Ishihara and Hunt2022). This issue is, however, outside the scope of the present study: our interest is rather on how the particle dispersion compares with the fluid tracer dispersion, as we now discuss.

The observable width of the super-diffusive regime depends on the scale separation between the forcing scale $L_F$ and the initial separation. For separations much larger than $L_F$ , the behaviour ultimately becomes diffusive, with a corresponding time scale roughly given by $t_{L,p} \sim L_F/u_{rms,p}$ . In our analysis, the temporal scale separation between the inertial (ballistic or super-diffusive) regime and the diffusive regime scales approximately as $t_{L,p}/t_{\zeta } \sim L_{F}^{-1/3}$ , modulated by the velocity ratio $u_{rms,f}/u_{rms,p}$ at the given flow $Re$ . Thus, increasing $L_F$ would, in a purely temporal sense, slightly reduce the window over which super-diffusion is observed. Conversely, from a spatial perspective, it is crucial that $L_F$ is significantly larger than the typical initial separation (which is lower bounded by the particle diameter ${d}_p$ ); indeed, our results indicate that, if the initial separation exceeds approximately $0.25L_F$ , the $t^3$ -scaling is not apparent and replaced by ballistic scaling ( $\sim t^{2}$ ). Moreover, when the particle velocity $u_{rms,p}$ is much lower than the fluid velocity $u_{rms,f}$ – as observed at high areal fractions – the effective temporal separation is enhanced, thereby extending the super-diffusive regime.

Figure 16. Examples of the MSS of particle pairs with an initial separation of ${d}_p = 1.84$ mm: (a) $\phi = 0.27$ , and (b) $\phi = 0.71$ . The dashed line corresponds to the MSS of tracers from the most turbulent flow in the experiments. (c) The MSS curves for high $\phi = 0.62$ across varying $Ca$ values, with time normalised by the collision time scale $t_{col} = h / u_{rel}$ and MSS by $h^{2}$ . The dotted line indicates the MSS prediction by (3.12). (d) The MSS curves at $Ca = 2.49$ across different areal fractions, with lines indicating the MSS predictions by (3.12) at $\phi = 0.09$ , 0.27, 0.44 and 0.71 (dash-dot), respectively.

The relative dispersion of the particles is illustrated in figure 16 for different values of $Ca$ and $\phi$ , again for $r_0 \sim {d}_p$ . A dashed line showing the MSS of tracers for a representative case (weakly dependent on the forcing, as shown in figure 15) is included for reference. In the moderately dense case ( $\phi =0.27$ , figure 16 a), we observe slower initial dispersion at $t/t_{\zeta } \lt 1$ , due to capillarity and lubrication inhibiting separation at small distances. Once particle pairs grow far enough apart that cohesive interactions become negligible, the curves at different $Ca$ collapse and exhibit a super-diffusive scaling for over a decade, until the separation grows beyond $L_F$ .

At $\phi = 0.71$ (figure 16 b), the effect of $Ca$ becomes more prominent. At $Ca \ll 1$ , as the particle velocity is significantly slower than the fluid velocity (see figure 12 a), it takes much longer than $t_{\zeta }$ for a pair to separate enough to break away from their mutual attraction. Conversely, at $Ca \gt 1$ particle separation occurs more rapidly and the super-diffusive regime is retrieved. This scaling of relative dispersion in such densely packed cases is striking, as most particles remain contained within a massively large, percolating cluster. This highlights the dynamic nature of clusters formed at high $Ca$ , in contrast to the quasi-rigid clusters formed at low $Ca$ . This change in behaviour is not captured by the analysis of cluster lifetime, which shows no noticeable trend with $Ca$ at $\phi = 0.71$ (see figure 11 b).

When using $t_{\zeta }$ for normalisation, the transition to the super-diffusive regime at dense concentrations occurs at different times. This is not surprising, as the fluid enstrophy is unlikely to be the governing factor at such high particle concentrations. We hypothesise that, for dense suspensions at high $Ca$ (i.e. in the lubrication-driven clustering regime), the frequent inter-particle interactions via lubrication, rather than the fluid enstrophy, become the dominant mechanism governing the energy dissipation and in turn the relative dispersion. For each particle, the energy dissipation rate due to lubrication is $\varepsilon _{lub} \sim F_{lub} u_{rel}$ , where $F_{lub}$ is the lubrication force experienced by a particle moving with a velocity $u_{rel}$ relative to its neighbour. Considering the particle number density per area $n = 4\phi /(\pi {d}_p^{2})$ and the area density $\rho _A = \rho _f \delta$ , where $\delta$ is the fluid layer thickness, the relevant short-range energy dissipation rate per unit mass can be approximated as

(3.11) \begin{equation} \varepsilon _{lub} \approx F_{lub} u_{rel} \frac {n}{\rho _A} \approx \frac {3\pi \mu {d}_p^{2} u_{rel}^{2} }{2h} \frac {4\phi /(\pi {d}_p^{2}) }{\rho _f \delta } = \frac {6 \nu _f u_{rel}^{2}}{\delta {d}_p} \frac {\phi }{[(\phi _{max}/\phi )^{0.5} - 1]}, \end{equation}

where $h = {d}_p[(\phi _{max}/\phi )^{0.5}-1]$ is the mean surface-to-surface separation between neighbouring particles, with $\phi _{max} \approx 0.9064$ being the maximum areal fraction from the 2-D hexagonal close packing.

Under the assumption that $\varepsilon _{lub}$ is the main governing parameter, dimensional arguments yield immediately $MSS \sim \varepsilon _{lub} t^3$ , and therefore

(3.12) \begin{equation} MSS(t) \sim \varepsilon _{lub} t^3 \sim \frac{6\nu _f u_{rel}^2}{\delta {d}_p} \frac {\phi }{[(\phi _{max}/\phi )^{0.5} - 1]}t^{3} .\end{equation}

Here, the relevant length and time scale are associated with inter-particle collisions (meant in a general sense as direct interaction/contact), which we take as $h$ and $t_{col} = h/u_{rel}$ , respectively. Since the longitudinal second-order structure function $S_{2,p}^L/u_{rms,p}^{2}$ of particles in the lubrication-driven clustering regime shows excellent quantitative agreement with $S_{2,f}^L/u_{rms,f}^{2}$ of tracers (figure 14 c), we approximate $u_{rel}(r) \sim u_{rms,p} (S_{2,f}^L(r)/u_{rms,f}^{2})^{0.5}$ , evaluated at $r=h$ . These length and time scales correspond to the distance and time over which a particle in a pair travels to encounter another neighbour whose centre is $h+{d}_p$ away from its own. For $\phi \gt 0.4$ , the particles in densely packed suspensions are expected to forget their initial relative velocity through mutual interactions, with lubrication governing the decorrelation of their relative velocity.

This expectation is confirmed in figure 16(c,d), which plots the MSS curves for $\phi =0.62$ with varying $Ca$ values, and for $Ca=2.49$ with varying $\phi$ values, respectively, normalised by $t_{col}$ and $h^{2}$ . Figure 16(c) shows that this normalisation produces a tight collapse of the curves, which display a super-diffusive regime for $t \gt t_{col}$ . Moreover, the super-diffusive MSS estimated for the highest $Ca$ using (3.12), based on the assumption that lubrication is the dominant energy dissipation mechanism in dense cases (indicated with a dotted line), shows good agreement with the experimental curve. We stress that such a scaling argument is expected to apply within a pre-factor of order unity. The agreement between (3.12) and the experimental MSS implies that $\varepsilon _{lub}$ takes on the role of the turbulent dissipation rate in single-phase flows.

To further inspect the role of short-range interactions in MSS across different concentrations, we plot the normalised MSS curves at $Ca = 2.49$ across varying $\phi$ values in figure 16(d). Dashed lines calculated using (3.12) indicate how the agreement progressively improves for increasing concentrations, as one transitions from the drag-driven break-up to the lubrication-driven clustering regime.