# Characterizing inertial particle clustering with 1D signals

Daniel Odens Mora, Alberto Aliseda, Alain H. Cartellier, Martin Obligado

► To cite this version:

Daniel Odens Mora, Alberto Aliseda, Alain H. Cartellier, Martin Obligado. Characterizing inertial particle clustering with 1D signals. Measurement Science and Technology, 2023, 35 (2), pp.025302. 10.1088/1361-6501/ad031f. hal-04805752

HAL Id: hal-04805752

https://hal.science/hal-04805752v1

Submitted on 25 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Characterizing 1D Inertial Particle Clustering

Daniel Odens Mora, Alain Cartellier, and Martin Obligado*

Univ Grenoble Alpes, CNRS, Grenoble-INP,
LEGI, F-38000, Grenoble, France †

Alberto Aliseda

Department of Mechanical Engineering,
University of Washington, Seattle, Washington 98195-2600, USA

(Dated: June 3, 2020)
Abstract

Preferential concentration is a common phenomenon found in turbulent flows seeded with inertial
particles. Although it has been studied extensively, there are still many open questions about its
fundamental physics. These gaps hinder the reconciliation of different experimental and numerical
studies into a single coherent quantitative view, which is needed to enable accurate high resolution
modeling.

This work examines the influence of the measuring technique, and in particular of the dimen-
sionality of the measurement (1D line, 2D planes or 3D volumes) on the characterization of cluster
properties (a consequence of preferential concentration), and proposes an approach to disentangle
the cluster-characterizing results from random contributions that could contaminate the cluster
statistical analysis.

First, we studied this effect by projecting 2D and 3D data snapshots containing clustering onto a 1D axis. The objective was to simulate 1D sensors (widely used experimentally) with different sensing lengths.

These projected records were analysed by unidimensional Voronoï tessellations. Our results
demonstrate that average mean clustering properties could be retrieved, if the measuring window
is equal or larger than the Kolmogorov length scale (η), and smaller than about ten percent of the
integral length scale of the turbulence L. These observations are consistent with 2D and 3D data
taken under similar experimental conditions. This agreement validates our projection approach as
it adequately captures the behavior of a 1D probe immersed in a 2D, or a 3D flow.

Additionally, we found that in 1D the raw probability density function (PDF) of Voronoï cells
does not provide error-free information on the clusters size or local concentration. We propose
a methodology to correct for this measurement bias, based on the histograms of the number of
particles within a cluster. The analysis of these histograms helps to explain the power-law behavior
previously observed in the clusters size PDFs in 2D and 3D data. Moreover, the histogram analysis
shows that such power law is not very robust in 1D, as it expands less than one decade.

Finally, we show that to condition the statistics with the number of particles inside each cluster
also allows to discern between turbulence-driven clustering and particle concentration fluctuations
due to randomness. Our test then complements the classical cluster identification algorithm.

<footnote>* Martin.Obligado@univ-grenoble-alpes.fr</footnote>
# I. INTRODUCTION

Turbulent flows laden with inertial particles constitute an active research area within multiphase fluid mechanics due to their potential applications in fields such as: planetary formation, pollutant dispersion, and cloud formation [1, 2]. Hence, experiments are not only useful to validate the multiple numerical approaches, but also to disentangle their underlying physics, and finally reach comprehensive understanding.

Several methods are available to characterize particle-clusters, with Voronoï tessellations [3–6] becoming increasingly popular in both experimental studies employing visualization techniques, e.g., particle tracking velocimetry (PTV) [4, 7] and numerical simulations [8, 9].

Despite its current popularity, it was early recognized [10] that there could be pitfalls in the Voronoï method that could affect its results. The open questions are: to what extent physics from a 3D particle spatial distribution measured by 2D or 1D techniques can be captured by the Voronoï method, and how the measurement configuration (1D, 2D, line diameter, plane thickness, etc.) could impact its results. The study of Monchaux [10] was the first to survey these possible biases in context of particle laden-flows. By projecting a 3D numerical dataset onto a 2D plane, they studied the effect of the laser sheet thickness ($L_{th}$ in their notation) on the resulting 2D Voronoï statistics. The study concludes that 2D Voronoï analysis is reliable for the most common values of laser sheet thicknesses ($L_{th} \in [2\eta - 6\eta]$) used in 2D imaging experiments, and very robust under random thinning, which simulates the effect of missed particles. However, care is needed when comparing 2D Voronoï statistics from different data sets obtained under different experimental or optical conditions, as sub-sampling results heavily depend on the particle seeding and the turbulence scales.

On the other hand, a recent study [11] reports that a 1D Voronoï analysis (referred also as 1DVOA from now on) performed on a record taken by an optical probe via phase detection [12], is unable to capture preferential concentration in experimental conditions under which clustering is observed using 2D Voronoï analysis [7]. Thus, 1D and 2D Voronoï analyses could yield contradictory results regarding the presence of preferential concentration within the flow.

The study of Mora et al.[11] suggests that the probe's measuring volume (a region in space where the instrument can detect particles transiting) has an important impact on this

<footnote>† Also at Univ Grenoble Alpes, CNRS, Grenoble-INP, LEGI, F-38000, Grenoble, France</footnote>
problem. This study reveals that below a certain threshold of the measuring volume, the 1D
Voronoï method is unable to capture the spatial correlations between the recorded particles.
This observation, however, is strongly affected by the particle concentration, as it was proven
in 2D projections [10], and therefore, it cannot be attributed to insufficient resolution of the
measuring instrument. Indeed, it is found that while an optical probe (with a measuring
length 10 times smaller than the Kolmogorov length scale of the flow) was unable to retrieve
clustering, it did recover the local liquid fraction, the particle size, and the particle velocity
distributions consistent with previous measurements.

Biases involving 1D measurements have also been found in other methods to determine
the existence of preferential concentration. For instance, when using 1D pair correlation
functions, it has been found that excluded volume effects (due to the finite size of the
droplets) could accumulate and had a combined effect across a range of scales, biasing the
results obtained by this method [13, 14]. There is, however, a fundamental difference between
1D pair correlation functions and 1D Voronoï tessellations. The former method subdivides
the dataset domain into segments of a certain scale, whereas, the latter method conducts
analysis that is considered to be scale free. Excluded volume effects do not accumulate in
1D Voronoï tessellations given that 1D Voronoï cells are computed from their immediate
neighbours. Conversely, the finite particle size impacts the value of the criterion used to
assess the presence of clustering; clustering is present if the standard deviation of the Voronoi
cell size σ<sub>v</sub> coming from the particle dataset is larger than the respective one coming from
the equivalent random Poisson distribution (RPP). Uhlmann has reported [15–17] via 3D
Voronoï tessellations that, for spherical particles with diameters up to 5η, the respective
magnitude of RPP standard deviation is reduced by 10%. These observations imply that
the nature of the pitfalls found in pair correlation functions, and 1D Voronoï tessellations,
although related, are distinct.

Given the available technology, and experimental apparatuses, one could question the
use of a ‘1D’ technique to characterize a 3D phenomenon. Still, these complementary ex-
perimental techniques, for example, phase detection optical probes [12], or phase doppler
interferometry (PDI) [18], have been proven useful, as they provide additional information
not so easily available from 2D/3D measurements. Moreover, when identifying the actual
droplet spatial distribution within a cloud [1, 2, 19] or probing pollutant dispersion within
cities [20, 21], in situ measurements are usually acquired by 1D eulerian instruments.
Hence, it is worth examining whether the 1D Voronoï statistics have any extra inherent biases when used to quantify preferential concentration, i.e., when the approach to analyze 2D or 3D data is extended to 1D records. The latter approach usually consists on building the Voronoï sizes probability density functions (PDFs), and then draw conclusions based on the particle local concentration [4, 7, 8, 15, 16, 22], which is easily available via Voronoï tessellations. The presence of these biases is not restricted to 1D probes. Analogue problems could arise in 2D and 3D measurements. The 2D Voronoi analysis, however, has been shown to be very robust for different laser sheet thicknesses and concentrations [10].

In this work, we study the impact of the instrument measuring window size (MWS) on 1DVOA statistics conducted on ‘projected’ records from numerical and experimental data sets. Then, different MWSs ‘emulate’ eulerian measurements taken with probes of different sizes in the same particle-laden flow. The later analysis allowed us to provide a range of MWSs (or probe sizes) for which evidence of preferential concentration could be consistently retrieved, under similar experimental values of the Taylor-based Reynolds number ($Re_λ$), and global liquid fraction ($φ_v$).

According to our analyses, this range’s lower bound (labeled $MWS_†$) is of the same order of the Kolmogorov lengthscale ($η$). In other words, our analyses suggest that to retrieve clustering for similar values of $Re_λ$, and $φ_v$ one should measure with a probe with a window size of order $η$, i.e., $MWS_† = O(η)$. On the other hand, we observe that for the analyzed experimental data and DNS (Direct Numerical Simulation) dataset this range’s upper bound (labeled $MWS_⋆$), after which evidence of clustering becomes weaker, seems to be located at measuring window size close to a tenth of the integral length scale. $MWS_⋆ \sim L/10$ where $L$ is the integral length scale of the carrier phase.

Furthermore, to characterize the clusters properties, we looked for biases in the probability density functions (PDFs) of cluster local concentration, and linear size. Our data strongly suggest that the method developed to analyze preferential concentration in 2D [23] cannot be (in general) directly extended to 1D Voronoï tessellations analyses. Such direct extension, under some experimental conditions, is not conclusive regarding the presence of clustering within the flow, and undermines the capacity of 1DVOA to characterize the phenomenon in 1D unambiguosly. The origin of this limitation could be attributed to the loss of information intrinsic to the 1D measuring technique, and it is supported by the weaker correlation seen among the recorded particles positions.
To further explore these biases, we developed a theoretical model to compute the PDFs of cluster sizes coming from particle positions spatially distributed according to a random Poisson process (RPP). The model underlying principles suggest that conditioning the clusters PDF by the number of particles within clusters could aid to properly find and characterize ‘turbulence’ driven clusters. The method, therefore, allows to disentangle randomness (as high density regions are also present within a RPP distribution due to random fluctuations of the local concentration) from turbulence (where clusters are controlled by the topology of the turbulent flow) in a 1D signal.

## II. THE EXPERIMENTAL SETUP AND METHODS

### A. Experimental setup and numerical database

The experiment was conducted in the close-circuit wind tunnel ‘Lespinard’ that has been extensively used to study particle clustering under turbulent conditions [7, 23]. A sketch of the experimental setup is depicted in figure 1a.

Turbulence was produced by means of an active grid [24] operated in triple random or open mode [25]. At the measuring station (region #3 in fig 1a), the unladen turbulence has been experimentally found to be very close to a statistically homogeneous isotropic state (figure 1b) under similar conditions of $Re_λ$, and $\eta$ [26] to those found in the present work.

Downstream of the grids, a rack of 36 spray nozzles generated inertial water droplets with a polydisperse diameter distribution. This polydispersity (see figure 1c) has been previously quantified via phase doppler interferometry [7].

At the measuring station, a high speed camera collected 4500 images of the light scattered by the droplets from a 1-mm thick laser plane. The images, with an area of (120×100 mm²) and collected to be statistical independent realizations, were post-processed to identify the location of the droplet centers (more details about the experimental setup can be found in [27]).

Along with experimental data, 30 numerical snapshots from a DNS data base [28] (https://turbase.cineca.it/init/routes/#/logging/view_dataset/3/tabmet) were surveyed. Each file contained the trajectories of 1280 particles positions integrated over 3300 time steps.
<figure><img src="image_2.png" /><figcaption>FIG. 1: a) Sketch of the wind tunnel experimental setup. 1, 2, and 3 refer to the locations of the active grid, the injection rack, and the measurement region, respectively. The measuring region downstream distance was taken from the beginning of the injector rack.</figcaption></figure>

The shaded region illustrates the extend of the laser sheet. The transverse square cross-section has dimensions of 750 × 750 mm². b) Velocity power spectral density for the active grid (AG) Re_λ ≈ 250, and for the open grid (OG) Re_λ ≈ 30. Both spectra were obtained via hot-wire anemometry. c) Droplet Diameter D_p distribution. The symbol (o) refers to data from [7], and the (—) line refers to a log fit (parameters shown in the plot legend).
<table><caption>TABLE I: Main turbulence parameters for experiments and DNS data used on this work. $λ = u'\sqrt{15ν/ε}$ is the Taylor length scale, and $ν \sim 1.5 \times 10^{-5} [m^2s^{-1}]$ is the air viscosity. $u'$ is the rms of the streamwise velocity fluctuations. The turbulent Reynolds number is defined as $Re_λ = u'λ/ν$ and $D_p$ is the value of the most probable diameter which was used to compute the Stokes number $St_η = \frac{(D_p/η)^2}{36}(1 + 2ρ_p/ρ_f)$, see [7]. $ρ_p/ρ_f$ is the density ratio between the particles, and the carrier phase. $ε$, and $L$ are the carrier dissipation, and integral length scale, respectively. AG/OG stands for the random or open mode of the active grid. $φ_v$ is volume fraction for the experimental data. As numerical data is a '1-way' coupling point particle simulation, we assume it is very diluted $φ_v \ll 10^{-6}$ [31]. More details on how the experimental unladen flow parameters have been calculated can be found on [25].</caption><thead><tr><th>Dataset</th><th>Grid Mode</th><th>Re<sub>λ</sub></th><th>St<sub>η</sub></th><th>εL<sup>4</sup>/ν<sup>3</sup></th><th>D<sub>p</sub>/η</th><th>L/η</th><th>λ/η</th><th>φ<sub>v</sub></th><th>ρ<sub>p</sub>/ρ<sub>f</sub></th></tr></thead><tbody><tr><td>EXP-2D-AG-A</td><td>Random</td><td>250</td><td>0.9</td><td>4.3 × 10<sup>8</sup></td><td>0.125</td><td>175</td><td>35</td><td>1.2 × 10<sup>-5</sup></td><td>800</td></tr><tr><td>EXP-2D-AG-B</td><td>Random</td><td>250</td><td>0.9</td><td>4.3 × 10<sup>8</sup></td><td>0.125</td><td>175</td><td>35</td><td>2.3 × 10<sup>-5</sup></td><td>800</td></tr><tr><td>EXP-2D-OG-A</td><td>Open</td><td>30</td><td>0.05</td><td>9.0 × 10<sup>4</sup></td><td>0.032</td><td>20</td><td>7</td><td>1.2 × 10<sup>-5</sup></td><td>800</td></tr><tr><td>EXP-2D-OG-B</td><td>Open</td><td>30</td><td>0.05</td><td>9.0 × 10<sup>4</sup></td><td>0.032</td><td>20</td><td>7</td><td>2.3 × 10<sup>-5</sup></td><td>800</td></tr><tr><td>DNS [29, 30]</td><td>-</td><td>185</td><td>1</td><td>1.1 × 10<sup>10</sup></td><td>-</td><td>314</td><td>26</td><td>-</td><td>≈ 1</td></tr></tbody></table>

These snapshots contained inertial point particles with a Stokes number close to unity; $St_η = τ_p/τ_η \approx 1$ where $τ_p = ρ_p D_p^2/18μ$ is the particle relaxation time, $τ_η = \sqrt{ν/ε}$ is the Kolmogorov timescale, $D_p$ is the particle diameter, $ε$ the turbulent dissipation rate, $μ$ the dynamic viscosity of the carrier phase, and $ρ_f$, and $ρ_p$ are the carrier and particle density, respectively.

The snapshots files can be found in the database labeled as: RM-2008-LIGHT-512.St6.XX.h5, which have particles with $τ_p = 0.048282$, and $β = 3ρ_f/(ρ_f + 2ρ_p) = 0$, i.e., $ρ_f \ll ρ_p$. More details of the numerical setup ($Re_λ \approx 185$) can be found in [29, 30]. The experimental and simulation parameters are summarized in tables I, and II.
<table><caption>TABLE II: Datasets summary. $N_{snap}$ is the number of snapshots, $N_p$ is the number of particles, $\langle N_p/N_{snap} \rangle$ is the average number of particles per snapshot, and $\sigma_{\langle N_p/N_{snap} \rangle}$ is the standard deviation of the mentioned average.</caption><thead><tr><th>Dataset</th><th>N<sub>snap</sub></th><th>⟨N<sub>p</sub>/N<sub>snap</sub>⟩</th><th>σ<sub>⟨N<sub>p</sub>/N<sub>snap</sub>⟩</sub></th></tr></thead><tbody><tr><td>EXP-2D-AG-A</td><td>4500</td><td>1 × 10<sup>3</sup></td><td>280</td></tr><tr><td>EXP-2D-AG-B</td><td>4500</td><td>5 × 10<sup>3</sup></td><td>1200</td></tr><tr><td>EXP-2D-OG-A</td><td>4500</td><td>3 × 10<sup>3</sup></td><td>330</td></tr><tr><td>EXP-2D-OG-B</td><td>4500</td><td>1 × 10<sup>4</sup></td><td>800</td></tr><tr><td>DNS</td><td>30</td><td>4 × 10<sup>4</sup></td><td>0</td></tr><tr><td>RPP</td><td>1000</td><td>10<sup>3</sup>, 10<sup>4</sup>, 10<sup>5</sup></td><td>0</td></tr></tbody></table>

## B. Voronoï Tesselations

Originally, Monchaux et al. [23] proposed the use of Voronoï tessellations to quantify preferential concentration in turbulent flows. In 2D, the Voronoï tessellation yields a collection of areas $A$ from particles spatial positions [3]. Statistics can be then computed for this collection, e.g., the average area $\langle A \rangle$, or the area standard deviation $\sigma_A$. The Voronoï analysis is usually conducted by considering the area collection normalized by its mean $\mathcal{V} = A/\langle A \rangle$, giving $\langle \mathcal{V} \rangle = 1$.

Monchaux et al. [23] suggested that Voronoï area PDFs can be characterized by their standard deviation ($\sigma_v$), a conclusion they substantiated by the log-normal behavior present in these PDFs, and the collapse seen under the log-normal transformation for different cases studied (with different values of $Re_\lambda$, or concentration $\phi_v$). These observations led them to conclude that $\sigma_v$ estimates the ‘intensity’ of clustering.

Then, to quantify ‘clustering’, $\sigma_v$ is compared with the respective standard deviation coming from particles spatially distributed following a RPP distribution ($\sigma_{RPP}$), which by definition has no spatial correlations at any scale. Thus, preferential concentration is present if $\sigma_v > \sigma_{RPP}$. It is clear that under this criterion, the voids (or area outliers of the area collection) contribute the most to the numerical value of $\sigma_v$ [7].

Likewise, the unidimensional Voronoï tessellation (1DVOA) generates a collection of
<figure><img src="image_3.png" /><figcaption>FIG. 2: Undimensional Voronoï tessellation (1DVOA). For a given particle position Z with left, and right neighbour particle $Z_L$, and $Z_R$ respectively, the length of the Voronoï cell is given by $L = |Z_R - Z_L|/2$.</figcaption></figure>

<table><caption>TABLE III: Summary of the Voronoï PDF expressions for 1,2,3 dimensions, and their respective properties as reported by Ferenc and Nedá [3] (1D and 2D) and Tanemura [32] (3D).</caption><thead><tr><th>Dimension</th><th>PDF expression</th><th>σ<sub>RPP</sub></th><th>Model</th></tr></thead><tbody><tr><td>1</td><td>4 × V<sup>1.0</sup>e<sup>-2.0V</sup></td><td>√1/2</td><td>Analytical</td></tr><tr><td>2</td><td>343/15√7/2π × V<sup>2.5</sup>e<sup>-3.5V</sup></td><td>≈ 0.53</td><td>Fit</td></tr><tr><td>3</td><td>345/7 × V<sup>3.8</sup>e<sup>-4.0V</sup><sup>1.17</sup></td><td>≈ 0.42</td><td>Fit</td></tr></tbody></table>

lengths $L$, not to be confused with the turbulence integral length scale $\mathcal{L}$ (figure 2). However, Mora et al.[11] reported that clustering could be present but not retrieved via 1D measurements [1, 13], and therefore, the criterion $\sigma_V > \sigma_{RPP}$ cannot be directly extended to 1D measurements, and it may not be entirely conclusive for 1DVOA.

For clarity, we will define some acronyms below, but will often recur to the acronym's verbatim. First, angle brackets $\langle...\rangle$ represent an ensemble average. Second: our short notation for the Voronoï tessellations variables goes as; $V = M_{voro}/\langle M_{voro} \rangle$ where $M_{voro}$ represents the length $L$, area $A$, or volume $V$ of a Voronoï cell. Thirdly, a 1D/2D/3D Voronoï tessellation analysis will be referred as 1DVOA, 2DVOA, or 3DVOA, respectively. Finally, the respective standard deviation of a RPP distribution ($\sigma_{RPP}$), which has no correlations at any scale, can be analytically computed: $\sigma_{RPP} \approx 0.71/0.53/0.42$ for 1DVOA/2DVOA/3DVOA, respectively (see table III, and references [3, 32]). The 2D and 3D Voronoï diagrams were computed using the python library *freud* [33].
### C. 1D ‘Projected’ Voronoï Analysis

Following the method proposed by Mora et al. [11], we projected particle positions coming from previous 2D and 3D data, which exhibit preferential concentration, orthogonally onto a line.

This approach aims to emulate the basic features of an Eulerian measurement: a fixed probe in space (e.g. downstream from the grid of a wind tunnel). Thus, the particle coordinates orthogonal to the bulk velocity ($U_∞$) were ‘frozen as the particles were advected by the flow figure 3a).

For the experimental 2D data sets, this projection process will be denoted as $2D_{EXP} \rightarrow 1D_{\perp}$, where $2D_{EXP}$ stands for 2D experimental images, and $1D_{\perp}$ is a uni-dimensional orthogonal projection into 1D, i.e., onto the streamwise axis $\gamma$ (figure 3a). For the numerical 3D data, the notation is ($3D_{DNS} \rightarrow 1D_{\perp}$), where $3D_{DNS}$ is the 3D dimensional DNS data set (figure 3b). We applied this procedure for each data snapshot (table II).

Thus, our ‘projection’ algorithm [11] comprised the following steps:

1. A random vertical $\hat{y}$ coordinate was generated. This coordinate set the position of the axis $\gamma$ onto which the particles were projected.

2. A symmetric measuring volume size (MWS) window was defined intended to quantify the effect of the instrument finite spatial resolution. At all times, the $\hat{y}$ coordinate choice was restricted so that the MWS width was inside the computational domain.

3. All the particle *centers* that lay within this measuring strip ($\gamma$ and width MWS) were projected onto the axis $\gamma$, i.e., their horizontal coordinate was recorded.

4. 1DVOA was performed over the particle location line projections.

5. Three additional elements were considered for the 3D data:

(a) A $\hat{z}$ transverse coordinate was generated to position a cylinder axis.

(b) The generated cylinder diameter was equal to the measuring window.

(c) All particles within this cylindrical volume were subsequently projected on the cylinder axis $\gamma$ (see figure 3b).
<figure><img src="image_4.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(c)</figcaption></figure>

FIG. 3: a) Sketch of the $2D_{EXP} \rightarrow 1D_{\perp}$ particles centers projection for an arbitrary image. MWS is the measuring window size, $\hat{y}$ is the randomly generated vertical coordinate of the axis $\gamma$ over which the points are orthogonally projected. b) Sketch of the $3D_{DNS} \rightarrow 1D_{\perp}$ particles centers projection for an arbitrary DNS snapshot. MWS is the measuring window size equal to the cylinder diameter, $\hat{y}$ and $\hat{z}$ are the randomly generated coordinates of the axis $\gamma$ onto which the points are orthogonally projected. c) Average number of projected samples per snapshot from 1D sampling ($2D_{EXP}$ or $3D_{DNS} \rightarrow 1D_{\perp}$), and its dependency with the measuring window size. At very small MWS with respect to $\mathcal{L}$ the average number of samples captured is small, which is directly linked to lack of clustering recently reported [11].
The previous algorithm seems to capture the basic features of an equivalent experimental eulerian measurement, and is suitable to examine the effect of different measuring window sizes have on 1D Voronoï statistics. For instance, most quasi-unidimensional instruments [12, 18, 34] are particle ‘counters’, which usually yield a list of events with their respective eulerian ‘arrival’ time. One would then expect that the larger the window size (the region in space where the instrument detects the particle transiting), the larger the number of droplets is ‘detected’, a behavior this ‘projection’ approach does indeed capture (see figure 3c).

### III. RESULTS

#### A. PDFs and Voronoï cell standard deviation

We first checked if preferential concentration was present in the 2D and 3D data. Using the same 2D datasets, Mora et al.[11] and Obligado et al.[27] have found strong evidence of preferential concentration by means of planar Voronoï tessellations. Next, we verified the existence of clustering in the 3D numerical database [29, 30] by means of tridimensional Voronoï Tessellations [3].

Visual inspection of the probability density function (PDF) (see figure 4a) of the normalized Voronoï cell volume ($V = V/\langle V \rangle$) suggests the presence of preferential concentration within the DNS data. This is further confirmed by the larger value of the standard deviation of $V$ (following the criterion proposed by Monchaux et al. [23]) with respect to the RPP distribution, which has no correlation at any scale (see section II B), i.e., $\sigma_V \approx 0.62 > \sigma_{3D_{RPP}} \approx 0.42$, (see table III).

Having confirmed the presence of preferential concentration within the 3D data, we applied the projection algorithm to it (see section II C). Then, we conducted a unidimensional Voronoï tessellation analysis to these $3D_{DNS} \rightarrow 1D_{\perp}$ projections.

The probability density functions of the Voronoï cells sizes $V = L/\langle L \rangle$ computed for several measuring window sizes (MWS), reveal a clear trend on these PDFs with varying window size (figure 4b): the larger the window size was, the closer the PDF shape approached the respective RPP distribution.

The right tail of the Voronoï PDFs ($V \gg 1$) exhibits large changes with varying measuring
<figure><img src="image_7.png" /><figcaption>FIG. 4: a) Probability density function plot of 3DVOA for the DNS data [29, 30].
Following the criterion of Monchaux et al. [23], it is clear that the DNS data contains
clustering, as $\sigma_V \approx 0.62 > \sigma_{3D_{RPP}} \approx 0.42$ is larger than its equivalent one for a 3D RPP
distribution. b) PDF plot of $(3D_{DNS} \rightarrow 1D_{\perp})$ 1DVOA for several MWS.</figcaption></figure>

window size. The importance of this cannot be overstated given that the criterion proposed
by Monchaux et al.[23] to determine the existence of clustering relies on the numerical value
of standard deviation $\sigma_v$, which in turn heavily depends [7] on the cell values much larger
than the mean $\langle L \rangle$, i.e., $\sigma_v = \int_0^\infty (V - \langle V \rangle)^2 PDF(V) dV$.

Furthermore, when computing $\sigma_v$ for our 2D, and 3D projections (figure 5a), these results show that the presence of ‘clustering’ via this criterion can only be recovered above a certain window size ($MWS_\perp$), below which the evidence of clustering is not conclusive, i.e., $\sigma_v/\sigma_{RPP} \approx 1$. In fact, for some data, $\sigma_v/\sigma_{RPP} < 1$ could lead to the wrong conclusion of a more ordered sub-poissonian distribution [1, 11, 13].

Given the smaller number of samples available with the projection method at smaller
window sizes (see figure 3c), one could attribute the result $\sigma_v/\sigma_{RPP} < 1$ to insufficient
statistical convergence. On the contrary, having statistical convergence does not necessarily
guarantee that the spatial correlations within the data are captured [13].

The latter can be illustrated by applying the same projection algorithm to a three-
dimensional random set (labeled $3D_{RPP} \rightarrow 1D_{\perp}$ in table II). Although this data set is,
by definition, random, it shows a similar a transition region where $\sigma_v/\sigma_{RPP} < 1$ (figure 5b).
These data sets attain their ‘theoretical’ value of $\sigma_v$ after a certain size which depends on to
<figure><img src="image_8.png" /><figcaption>FIG. 5: a) 1DVOA standard deviation evolution for the datasets used. The larger the concentration, the higher $\sigma_V$ for fixed measuring window size. The peak location follows the relation $\text{MWS}_{\star} \approx \mathcal{L}/10$, where $\mathcal{L}$ the integral length scale of the flow (see table I). However, its value depends on $Re_{\lambda}$, as these and previous studies have shown [7]. b) 1DVOA standard deviation evolution of projections coming from synthetic random 3D data ($3D_{RPP} \rightarrow 1D_{\perp}$). $N_P$ stands for the number of points inside the 3D domain for 1000 synthetic snapshots.</figcaption></figure>

the number of particles in the domain (concentration). Hence, it is the instrument and its capacity to capture particle spatial correlations which ultimately determines to which extent the 1DVOA is successful. Accurately identifying these correlations yields the proposed picture of clusters and voids [23], as there cannot be voids without clusters. Moreover, if evidence of clustering ($\sigma_V/\sigma_{RPP} > 1$) is retrieved by means of 1DVOA, it is reliable, as a strong signature of clustering is present even after randomly removing 70% of the data points (figures 6a and 6b).

The previous observations are very relevant from an experimental point of view, as they aid to explain why two distinct measuring instruments with dissimilar window sizes could yield opposite results; as seen in the experiments of Sumbekova [35] and Mora et al. [11], ran under almost identical conditions but with different 1D measuring instruments. For the former experiment, a PDI device ($\text{MWS} \approx \eta$) was used, while for the latter, optical probes were used ($\text{MWS} \approx 0.1\eta$). This order of magnitude difference in the instrument window size is key to understand the origin of this contradiction (figure 5a).
<figure><img src="image_9.png" /><figcaption>FIG. 6: 1D PDI data from Sumbekova [35] for varying degrees of thinning (randomly removing samples from their records). The larger the percentage the larger the number of original particles removed. a) Thinned Voronoï PDFs b) Their corresponding normalised standard deviation ($\sigma_V/\sigma_{RPP}$).</figcaption></figure>

These results, therefore, indicate that a sensitivity analysis has to be performed if a uni-
dimensional Voronoï analysis retrieves no clustering where higher dimensional techniques
do (2DVOA or 3DVOA). However, considering the observed strong dependency of the mini-
mum window size (MWS$_{†}$) on the concentration values $\phi_v$ for which $\sigma_V/\sigma_{RPP} > 1$ (figures 5a
and 5b), and the strong dependency of $\sigma_V/\sigma_{RPP}$ on the carrier phase conditions Re$_{\lambda}$, such
sensitivity analysis is not straightforward, and it should be tuned to the specific conditions
of the experiment.

Interestingly with our projection algorithm, we also recovered that the degree of clustering
increases with $Re_\lambda$ in agreement with previous studies [7, 22]. Our results, and the data of
Sumbekova [35] suggest that an instrument with a measuring volume equal or larger than
the Kolmogorov length scale ($\eta$) should be used if evidence of preferential concentration is
to be retrieved by means of 1D Voronoï tessellation analysis. Likewise by means of RDF
functions, Saw et al. [36] have recovered evidence of preferential concentration using a probe
with a volume of order $\eta$.

Although the choice of $\eta$ as MWS$_{†}$ could be seen as arbitrary, it seems to be adequate for
the current range of $Re_\lambda$ attainable in several experimental facilities (and in our data sets)
under similar experimental conditions. More importantly, this criterion does not contradict
previously reported clustering at scales smaller than η [37]. If these correlations were to be
present at smaller scales than η, they would be effectively captured (figure 5a).

Despite being only intended to capture the basics of an eulerian measurement, our ap-
proach also seems to capture interesting information related to preferential concentration
physics, giving extra credence to its underlying hypotheses.

On the one hand, it is remarkable that the peak of the ratio σ<sub>V</sub>/σ<sub>RPP</sub> (maximum clustering
degree in figure 5a) occurs at a measuring window size close to a tenth of the integral length
scale, i.e., MWS<sub>⋆</sub> ≈ L/10. This peak is representative of the maximum spatial correlation
present at the data set, as our projection procedure does not increase the actual correlation in
the data (see 3D RPP projections results in figure 5b). Hence, this window size, or scale, can
be connected to the multi-scale mechanisms proposed to explain preferential concentration
[38, 39].

For instance, the sweep-stick mechanism [38, 40, 41] considers that the scaling behavior
for voids (the counterpart of clusters) in the inertial range is not described by a single
scale, but instead follows a self-similar behavior within a range of scales (window size),
ℓ<sub>max</sub>/ℓ<sub>min</sub> = O(10) (where ℓ denotes scales from the inertial range of the turbulent flow).
We can therefore interpret MWS<sub>⋆</sub> as representative of the maximum interaction among all
scales between the carrier turbulent flow and the particles. Also, the study of Bragg et al.
[42] reports that their proposed mechanism responsible of clustering attains its peak at a
scale close to r = 200η, and with L = 800η, these observations are in the same order of
magnitude of MWS<sub>⋆</sub> ≈ L/10.

B. Average clusters size and measuring window size

We will now examine the effects that the observed biases have on the cluster linear size $L_C$. We started by computing $L_C$ following the most widely accepted cluster identification algorithm [23]: first we selected all the normalized Voronoï cell sizes $\mathcal{V}$ that were below a threshold value, $\mathcal{V} < \mathcal{V}_{th}$. This threshold is defined as the closest crossing (to the left of the RPP peak) between the RPP PDF, and our data PDF (see figure 7a), i.e., $\mathcal{V}_{th} = \mathcal{V} = \mathcal{V}|_{RPP} < 1$ [10, 23, 26].

The second step in the cluster algorithm is to find (within the collection $\mathcal{V} < \mathcal{V}_{th}$) cells
which share at least one edge. Then, two or more neighbouring cells were considered to be
a cluster ($N_{PC} \ge 2$, where $N_{PC}$ is the number of particles inside the cluster).

We took the same Voronoï cell threshold ($V_{th} \approx 0.55$) for all 2D and 3D projections. It is important to note that this cluster identification algorithm is applied individually to each snapshot (see section IIC and table II). Otherwise, spurious results in the different metrics, not shown here, could arise.

The results (see figure 7a & 7b) for the different data sets, and constrained to the MWSs for which $\sigma_V/\sigma_{RPP} > 1$, show two distinct regions: one where the average cluster size has a power-law decay behavior with increasing measuring window, and another where $\langle L_C \rangle/\eta$ evolves slowly with MWS.

The former region occurs at larger window sizes than the minimum for which $\sigma_V/\sigma_{RPP} > 1$ (see in figure 5a, the region where MWS$\gg$MWS$_{\dagger}$, the minimum window size). The power law behavior is a direct result of the projection method, and can be easily explained if one takes into account that: first, the average Voronoï cell value $\langle L \rangle$, which is inversely proportional to the particle number ($N_p$), i.e., $\langle L \rangle \propto N_p^{-1}$, also exhibits an analogous behavior with the window size (see figure 3c). Second, it has to be noted that the Voronoï cells below the threshold ($V < V_{th}$) are the ones that contribute the most to $\langle L \rangle$ [7]. Hence, it is not surprising that clusters become smaller (and denser) as the concentration is locally increased by projecting more particles onto a fixed domain; a leap in concentration yields a proportional reduction in $\langle L_C \rangle$, a potential setback on the applicability of 1DVOA for cluster characterization.

Given these considerations, we proceed to analyze the transition region for which $\langle L_C \rangle$ varies slowly with MWSs. The extent of this region seems heavily dependent on the liquid fraction $\phi_v$ and $Re_\lambda$. More importantly, the estimated average cluster lengths ($\langle L_C \rangle/\eta$) are again in agreement with previous 2D and 3D studies. Although early research shows [1, 15, 21, 39] that average cluster size is of order $10\eta$ for $Re_\lambda = O(100)$, recent research has shown [7, 45, 46] that $L_C$ may grow for larger values of $Re_\lambda$, for instance, for $Re_\lambda = O(300)$ [7, 22, 46] report $L_C = O(20 - 100\eta)$ in agreement with our observations.

In fact, the data (figures 5a and 7b) suggest that the maximum clustering window size $MWS_\star \approx \langle L_C \rangle \approx 0.1\mathcal{L}$. This conjecture seems to be supported by previously published data, for instance; Monchaux and Dejoan [9] reported $\langle L_C \rangle/\eta \approx 2-4$, with $\mathcal{L}/\eta \approx 30$, and $Re_\lambda \approx 30$, which is in the same order of magnitude of $\langle L_C \rangle/\eta \approx 0.1\mathcal{L}/\eta$; Obligado et al. [22] and Sumbekova et al. [7] outlined $\langle L_C \rangle/\eta \approx 10-20$ with $\mathcal{L}/\eta \approx 500$ and $Re_\lambda \in [200, 400]$;
and recent experiments from Petersen et al. [46] reported $\langle L_C \rangle / \eta \approx 10-40$ with $\mathcal{L}/\eta \approx 500$ and $Re_\lambda \in [200, 500]$, which are again in rough agreement with MWS\* $\approx \mathcal{L}/10$. The 3D DNS projected data also follows a similar trend, i.e., $MWS\* \approx \langle L_C \rangle \approx 40\eta$, which is close to the value found from 3D Voronoï tessellation analysis.

Thus, the expression $\langle L_C \rangle \approx \mathcal{L}/10$ could constitute a scaling for the cluster characteristic size, under similar experimental conditions. This latter consequence is remarkable, and consistent with the view of the strong role of turbulence in controlling the spatial correlation in the particle concentration field [40, 41]. For instance, despite having on average the same number of projected particles (see figures 3c and 7b), the two 2D data sets with the highest concentrations, $\phi_v = 2.3 \times 10^{-5}$, have characteristic sizes differing by almost an order of magnitude in agreement with previous 2D studies [7, 22].

Once again, these results give credence to our conservative recommendation for a minimum window size close to $\eta$, as such probe will not only be able to recover evidence of preferential concentration, but also will retrieve similar values and trends of $\langle L_C \rangle / \eta$, under similar experimental conditions of $Re_\lambda$, and $\phi_v$.

## IV. CLUSTER SIZE PDF

Having checked the effects of the measuring window size on the average cluster size $\langle L_C \rangle$, we proceed to study the sensitivity of the PDF of clusters length ($L_C$) to varying MWSs.

The analyses for experimental active grid data ($Re_\lambda \approx 250, 10^5\phi_v \approx 2.3$, the data from the other experimental conditions exhibited a similar behavior) yield two outstanding observations (see figures 8a and 8b). First, the lack of a conclusive power law on the right tail of the $L_C$ PDF. Second, there is a rough collapse for the PDFs of all window size with previous experimental data.

The absence of a power law behavior is surprising, as it was previously found in multiple 2D and 3D Voronoï studies [4, 8, 9, 16, 22, 46], and has been thought to be representative of the ‘turbulence-driven’ clusters, as the arguments of Goto and Vassilicos [38] suggest that the PDF of the voids areas should scale as $f_{AV} \sim A^{-5/3}$ based on Kolmogorov scalings, and a ‘-5/3’ powerlaw in the velocity power spectral density. Remarkably, such power law is also found for the PDF of clusters area $A_C$, or cluster volume $V_C$ in both 2D and 3D. For our unidimensional records, such power-law behaviour barely exists (if at all), as it can only be
<figure><img src="image_10.png" /><figcaption>FIG. 7: a) Probability density functions (PDFs) of $(2D_{EXP} \rightarrow 1D_{\perp})$ 1DVOA (EXP-2D-AG-B) for three different measuring window sizes (MWS). The ○ (SKA) markers correspond to the data of Sumbekova [35] for $10^5 \phi_v = 2$. b) Average linear cluster size $\langle L_C \rangle$ vs different measuring window sizes (MWS) for the data analyzed. The filled markers represent the axis on the left where the cluster average size is normalized by $\eta$. The open symbols correspond to the axis on the right, where the cluster average size is normalized by $\mathcal{L}$ (the integral length scale).</figcaption></figure>

identified for less than a decade for different datasets (figure 8b).

The second observation is somehow baffling, as one could not unambiguously anticipate that these projected data PDFs (numerical and experimental) would collapse at different MWSs. The existence of relative good agreement between the projected data and previous experiments is also intriguing [35, 44] (figure 8b).

This observation has also been reported for the PDF of cluster volumes ($V_C$) in numerical simulations. For instance, Uhlmann and collaborators [6, 15, 16] have also found by means of 3D Voronoï tessellations that if the same clustering algorithm (and with the same threshold values $V_{th}$, see section II B) is applied to data coming from a 3D random distribution, and to clustering containing data, both volume cluster PDFs attain (to some extent) a close resemblance. This pitfall is not minor, as it is critical to distinguish between the clustering coming from the turbulence and from the random fluctuations in the data. Otherwise, it is very difficult to correctly assess the influence of the turbulence on particle laden related variables such as particle settling velocity [4].
Hence, we conducted a similar analysis as the one of Uhlmann’s group to explore the origin of both observations: the lack of power law behavior, and the close PDFs resemblance in 1D. Our approach was to apply the 1DVOA to synthetic data from a uniform random distribution with a uniform distribution of probability, and applying the same threshold value to those Voronoï cells with sizes smaller than 0.55, i.e., $V < V_{th} = 0.55$. These cells were then processed using the cluster identification algorithm described in section III B. As random datasets have no actual clusters, the algorithm then detects the high concentration regions present on the RPP. This analysis will then be useful to compare these random fluctuations with the clusters generated by the interaction between inertial particles and a background turbulent flow.

The results of this analysis evince that the PDF of cluster sizes from 1D synthetic random data collapse remarkably well with the clusters PDFs generated from previous 1D experimental PDI particle records [35, 44] under turbulent conditions (see figure 8b). As this covers most of the available clusters data, our next objective is to check whether these observations are a result of the Voronoï tessellation algorithm in 1D, or if they are a consequence of the cluster detection method [23] in any dimension.

## A. Clusters PDF Model

We developed an analytical model for cluster size probability distribution $f(L_C)$ coming from synthetic 1D random data (a detailed explanation can be found in the A). Our aim with this model was to gain insight into: a very weak (if existent at all) power law, and a rather universal shape for the cluster size PDF for any input 1D data fed into the 1DVOA.

We used a mixture PDF model [47], a model based on PDFs linear superposition, and as such each PDF $f_i$ is multiplied by a weight $\alpha_i$, i.e., $f_{mix} = \sum_{i}^{N} \alpha_i f_i$, where N is the number of PDFs to combine. We then superimposed the PDFs of clusters of two, three, four, and up to N particles to construct a PDF which represents the magnitudes of $L_C$.

The question was, however, how to choose the coefficients $\alpha_i$. A sensible approach is to propose that these weights are proportional to the probability of clusters having $i$ particles within the cluster collection. In other words, $\alpha_{N_{PC}}$ = counts of clusters of size $N_{PC}$ / the total number of clusters. We then obtained these weights by computing histograms ($S_N$) conditioned on the number of particles in a cluster; $N_{PC}$ (see figure 9a, and equations A10
<figure><img src="image_11.png" /><figcaption>FIG. 8: In the figures, the markers (o) and (Δ) correspond to the experimental data from Sumbekova [26], and Bateson et al.[44] (taken in a different facility), respectively. a) PDFs of 1DVOA of linear cluster size $L_C/\langle L_C \rangle$ for several window sizes (MWSs). b) PDFs of 1DVOA of linear cluster size $L_C/\langle L_C \rangle$ for a random uniform distribution (RPP [3]), experimental data, and the model proposed. For all clusters having between 2 and 25 cluster cells ($N_{PC}$). The model proposed here represents well the right tail of the RPP data when $N_{PC} \ge 2$. The condition $V \le V_{th} = 0.55$ was employed for clustering computation.</figcaption></figure>

- A12).

Our PDF mixture model based on these weights (see open symbols in figure 9a), has good agreement with the clusters PDF found by applying the cluster identification algorithm to synthetic 1D RPP data (see figure 8b), and therefore, it strongly supports the absence of a power law decay in the cluster size PDF coming from a 1D RPP distribution.

Moreover, the construction of the histograms $S_N$ for turbulence-induced data reveals that $S_N$ histograms are power law distributed (see filled symbols in figure 9a) up to certain extent for the data containing preferential concentration both in 1D and 2D, whereas they decay exponentially for randomly generated RPP data both in 1D and 2D. On the contrary, 3D randomly generated data still shows a power-law similar to the one present for the turbulent DNS dataset.

It therefore appears that the existence of a power law behavior in the clusters PDF is correlated to the functional dependency of the weights distribution ($S_N$) in any dimension (see figure 9a). This observation then suggests that the individual PDFs ($f_i$) should attain
a particular shape at increasing values of $N_{PC}$. Given that our model estimates the PDF of $L_C|_{N_{PC}}$ assuming a sum of independent random variables, i.e., $L_C|_{N_{PC}} = X_1 + X_2 + X_3 + \dots + X_{N_{PC}}$ (see in the appendix figure 15c, and equations A1 and A9), we found that at increasing $N_{PC}$ this PDF tends to a Gaussian distribution, as stated by the central limit theorem [48].

Although this theorem is valid for sums of independently distributed random variables, it has been rigorously proven [49] (invoking strong-mixing conditions [50, 51]) that some consequences of it can be extended to sums of correlated variables having limited normal or gamma distributions; $f_Γ(\mathcal{V}, p, k) = k^p \mathcal{V}^{p-1} e^{-k\mathcal{V}}/\Gamma(k)$. Taking into account that; the Voronoï cells $(\mathcal{V})$ PDF (see table III) are special cases, or numerical fits of the gamma distribution (for instance, in 1D $k=p=2$ [3, 52]), and that sums of random correlated variables are implied when computing the clusters PDF, a central limit theorem [49] should hold for the cluster algorithm of Monchaux et al. [23]. We can then conclude that the histograms of $S_N$ indeed control shape of the cluster size PDF, and thereby, they provide a robust criterion to identify turbulence-driven clusters.

From these previous results, three conclusions can be drawn. First, the weak power law regime seen in the right tail of the cluster size PDF in 1D is due to the limited extend of such behavior in their respective 1D $S_N$ histograms, as shown by the $S_N$ behavior computed from the experimental records of Sumbekova [35], and Bateson et al.[44]. This behavior in the $S_N$ histograms is really robust and independent of the number of samples (for the PDI). For instance, the data of Sumbekova had $10^6$ samples, whereas the data of Bateson et al. had $10^4$ samples. Hence, this weaker power law behavior is more likely to come from the loss of correlation inherent to the 1D measurements, and not due to insufficient statistical convergence.

Second, the small ‘compact’ clusters sizes, which are on the left tail of the PDF, have also a behavior close to a power law, i.e., $f_C(\mathcal{M}_C/\langle\mathcal{M}_C\rangle < 1) \approx (\mathcal{M}_C/\langle\mathcal{M}_C\rangle)^a$, where $\mathcal{M}$ stands for length $L$, area $A$ or volume $V$, and $a$ is twice the absolute value of the algebraic exponent found in the base RPP Voronoï cell PDFs (see table III), i.e., $a \approx 2/5/8$ in 1D/2D/3D, respectively.

Third, the stronger power law behavior previously reported for the 2D/3D area/volume clusters PDF [8, 22] can be alternatively explained by the much wider extent (compared to 1D) of such powerlaw behavior in the histograms of $S_N$ (figure 9a). Our analysis is
<figure><img src="image_12.png" /><figcaption>FIG. 9: a) Histogram for number of particles within a cluster the different datasets. Filled symbols represent random data, whereas open symbols represent the turbulent datasets.</figcaption></figure>

b) Probability distribution function of clusters volume $V_C/\langle V_C \rangle$ for 3D DNS data, synthetic data, and the PDF mixture model using the weights found in figure 9a. The condition $\mathcal{V} \le \mathcal{V}_{th} = 0.62$ was employed for clustering computation.

also in agreement with the findings of Uhlmann and collaborators [6, 15, 16] who report the presence of a power law for the clusters volume PDF (3DVOA) when the clustering algorithm is applied to a random Poisson distribution with no correlations at any scale. The latter result further supports the strong role that the $S_N$ histograms have on the clusters PDF despite the increased uncertainty of our mixture PDF model in higher dimensions due to the approximations involved in the base Voronoï PDFs in 2D and 3D (for details see A 1, and table III).

Thus, it seems justified to use this model to analyze some properties, and possible biases found in 2D, and 3D Voronoï analyses, as our model for the PDFs of cluster areas/volumes have the same leading order behavior: with the absence/presence of power laws if the weights $S_N$ have exponential/power law distributions (see figures 9b, 10a and 10b).

For instance, in the context of 2D Voronoï tessellations, an additional step has been recently suggested to suppress spurious edge effects [46, 53]. According to these studies, these edge effects cause ripples in the cluster areas PDFs, $A_C/\langle A_C \rangle$, and therefore, one has to require that all cells surrounding the detected clusters have a cell size below the threshold $\mathcal{V}_{th}$.
<figure><img src="image_13.png" /><figcaption>FIG. 10: a) Probability distribution function of cluster area $A_C/\langle A_C \rangle$ for a random uniform distribution, and the model proposed using the weights found in figure 9a, and using weights from a ZIPF distribution $\alpha_i \propto N_{PC}^{-1}$. Indeed using the latter weights yields a strong power behavior. b) PDF of clusters area normalized by average Voronoï cell area $\langle A \rangle$ for the synthetic data, our model. The presence of ripples within the analytical PDF seems to challenge the argument of Zamansky et al. [53] to cater for the 2DVOA edge effects.</figcaption></figure>

Although we recovered ripples in the 2D cluster area PDFs (figures 10a and 10b), it is questionable that we could attribute the existence of these ripples to edge effects, as our analytical model (convolutions via Fourier transforms) should not be affected by them. From our model construction, it appears that these ripples occur at the boundary that mixes the individual PDFs, e.g., where the PDF for clusters containing 2 particles merges with the PDF for clusters containing 3 particles. Given that having clusters on these boundaries is less likely, it could explain such oscillations in the PDF. These ripples, however, disappear or become less ‘intense’ for larger cluster sizes containing many particles. From our analysis, filtering out these ripples is equivalent to leaving out small, power law dependent clusters.

## B. An approach to disentangle turbulence from random fluctuations

Considering the results of the previous section, a complementary methodology to the classical algorithm [23] is needed to discern turbulence-driven clusters from concentration fluctuations of random uncorrelated particles.
<figure><img src="image_14.png" /><figcaption>FIG. 11: PDFs of 3DVOA (DNS data of Baker et al. [8] for $St = 11$, and $S_v = \tau_p g/u' = 0$) for normalized cluster volume $V_C/\langle V_C \rangle$. Baker et al. suggest that turbulence driven clusters are the ones for which for $V_C > \eta^3$ or $V_C/\langle V_C \rangle > 1$. The figure shows that their criterion might fail (see figure 11b) for almost the first decade of clusters (several points). b) CDF of 3DVOA (DNS data of Baker et al. [8] for $St = 11$, and $S_v = \tau_p g/u' = 0$) for normalized cluster volume $V_C/\langle V_C \rangle$. The vertical dashed line (---) line represents up to the range where RPP-3D generated data follows a distinguishable power law (see figure 11a).</figcaption></figure>

With this goal, we briefly recount proposed fixes found in the literature, and analyze their fitness to solve the mentioned problem. Baker et al. [8] have proposed an amendment in the context of 3D Voronoï tessellations. Their underlying principle is that only the clusters lying on the right tail of $V_C/\langle V_C \rangle$ (see figure 11a) should be considered as turbulence-driven ‘coherent’ clusters. However, our data, and histogram analyses (figures 9a, and 9b) challenge this argument, since the clustering algorithm applied to synthetic 3D RPP data still yields a self-similar behavior for the clusters volume distribution ($V_C/\langle V_C \rangle$). This behavior lasts for at least a decade, before transitioning into an exponential decaying behavior in agreement with the $S_N$ analysis. Thus, by following strictly this criterion, a large number of ‘random’ (see figure 11b) clusters could remain after the additional step proposed, and be counted as turbulence-induced clusters (see figures 11a and 11b), defeating the purpose of such amendment.

An alternative test to distinguish randomness from turbulence-driven clusters has been proposed by Uhlmann and collaborators [6, 15]. It involves comparing the 3D clusters mean
aspect ratio when clustering is present, as turbulence-driven clusters would have larger
aspect ratio than randomly generated cells, due to gravity and enhanced particle settling.
This criterion, however, could become less accurate as the clusters become less columnar.

On the other hand, the analyses presented in the previous section suggest a methodology
based on segregated cluster probability $S_N$ (figure 9a). According to these analyses, one ex-
pects that the right tail of the turbulence-driven cluster size distribution (in any dimension),
conditioned on the number of particles in the cluster ($N_{PC}$) would preserve the powerlaw
behavior at increasing thresholds of $N_{PC}$; whereas for random data this behavior would be
eventually lost at increasing $N_{PC}$.

For the 1D case, this segregation approach seems promising, as despite having a very weak
power-law dependency (perhaps due to the limited amount of data from the records from
Sumbekova [35]), the clustering containing data seems to preserve its power law behavior
up to cluster of four or more points ($N_{PC} \le 4$) (see figure 12a).

In 2D, the segregation approach is very effective (figure 12b) as the right tail PDF power
law behavior (figure 9a) is conserved at increasing the values of $N_{PC}$. This is in agree-
ment with the $S_N$ analysis (figure 9a), and consistent with the mentioned conjecture that
turbulence-driven clustering will preserve its power law behavior at increasing values of
increasing $N_{PC}$.

In 3D, the clustering containing results apply to a larger extent their power law de-
pendency, and the respective 3D random data eventually decays exponentially at around
$N_{PC} > 20$ (see figure 12c), consistent with the transition found in $S_N$ (figure 9a).

Thus, an analysis of the clusters histograms conditioned on the number of particles $N_{PC}$,
is a promising complementary step to detect the presence of turbulence-driven clusters, as
it seems capable of detecting and avoiding the possible biases present in cluster identifying
algorithms [8, 23] in any dimension.

# V. CLUSTERS CONCENTRATION 1D

## A. Average cluster concentration

Knowledge of the particle local concentration is of the foremost importance in the study
of particle-laden flows [31], and its potential applications. The easy access to local concen-
tration maps is what gives Voronoï tessellations the upper hand with respect to different
available methods [10, 23].

However, taking into account the results from sections III A-III B, the probe measur-
ing volume will have an impact on the estimated local concentration values, when these
tessellations are applied to 1D records. In this section, we conduct a sensitivity analysis
on the projected data (analogous to the one found in section III B for the average cluster
size), and for the same data found in table I. We then define the cluster concentration as
$C_C = N_{PC}/L_C$, where $N_{PC}$ is the number of particles inside the cluster, and $L_C$ is its length.
The average concentration of the entire record is defined as $C_0 = N_P/L_R$, where $N_P$ is the
total number of particles detected over the total length recorded $L_R$, which is close to the
inverse of the average Voronoï length [3] $C_0 \sim 1/\langle L \rangle = f(\text{MWS})$.

The analysis by varying windows shows that the average cluster concentration value
$\langle C_C/C_0 \rangle$ has a non-monotonic behavior which depends on the measuring window size, and
the bulk liquid fraction $\phi_v$ (although only one dataset is shown in figure 13 the remaining
data sets -not shown here- exhibited the same behavior). There is also a transition region
similar to the one found in the normalized average cluster size $\langle L_C \rangle/\eta$ (see figure 7b), and
after this region the concentration tends to a value close to 3.

For the transition region, the data points with a measuring size of order $\eta$ are in good
agreement with previous 1D and 2D experimental data, i.e., $\langle C_C/C_0 \rangle \in [2.0 - 5.0]$ [23,
27, 35]. Although a comparison between the 1D and the 2D average cluster concentration
values $\langle C_C/C_0 \rangle$ was not conclusive due to the uncertainty in the lower to higher dimension
extrapolation [54], it seems that our approach also captures the leading order evolution of
the average cluster concentration.

The latter supports again that under similar experimental conditions an instrument with
with measuring window size of order $\eta$, is able to recover global average clustering parameters
which are in agreement with previous data captured with higher dimensional techniques.

On the other hand, at increasing MWSs values, the cluster concentration approaches
asymptotically to $\langle C_C/C_0 \rangle \approx 3.0$; cluster density varies proportionally with the MWS
in agreement with arguments of section III B. It is, however, important to notice that
$\langle C_C/C_0 \rangle \approx 3.0$ is very close to the value attained when the same algorithm is applied
to the denser regions of RPP distributions in any dimension. We can estimate this con-
centration analytically for the RPP distribution, if we take account that $C_0/C \approx \langle L \rangle|_{0}^{V_{th}} =$
$K \int_0^{V_{th}} V^2 e^{-2V} dV$ (K is a normalization constant see equation A2) for $V_{th} = 0.55$, $C_0/C \approx \langle L \rangle|_{0}^{V_{th}} \approx 0.332$, or $C/C_0 \approx 3$. Given the similarity between the RPP, and turbulence-driven values, care has to be taken when drawing conclusions based on average cluster concentration.

## B. Concentration PDFs

We now explore the concentration PDFs, and their sensitivity to variations with the measuring window size. Similar to the clusters PDF, the 1D concentration PDFs of synthetic random and clustering containing data have a close resemblance (figure 14a).

The concentration PDF also follows a -4 exponent, as reported by Sumbekova via PDI measurements [26]. This exponent, however, is a trivial result of the unidimensional testellation, as it can be analytically estimated; given that $C_C \sim L_C^{-1}$, so that $PDF_{C_C}(C_C) \sim PDF_{L_C}(1/C_C) \times C_C^{-2}$ by the chain rule [55]. Then, for the region $L_C/\langle L \rangle \ll 1$, $C_C \gg 1$, where the cluster size PDF$_{L_C}$ exhibits a power law with exponent of 2 (figure 8b), it follows that $PDF_{C_C}(C_C) \sim L_C^4 \sim C_C^{-4}$, and thereby, the -4 exponent.

Thus, to gain more insight into the effects of turbulence on $C_C$, we conditioned the clusters based on their number of particles, as proposed in section IV B. Although mild differences are found between clustering, and random generated data at very small $N_{PC}$, the preferential concentration data shows that it is almost an order of magnitude more likely to find very dense regions $C_C/C_0 > 4$ than in the RPP (see figure 14b) when $N_{PC} > 8$. Previous experiments have shown that particles belonging to ‘compact’ (very dense) clusters exhibit peculiar properties, for instance, enhanced settling velocity due to collective effects [4, 43]. Given these results, it is questionable to neglect those ‘compact’ clusters found at the left region of the clusters PDFs, as proposed by the amendment of Baker et al. [8] also discussed in section IV B.

As a final observation, it seems advisable that future experiments should acquire very long records in order to disentangle, and correctly study statistics conditioned on local concentration as 1D measurements are unable to register enough turbulence-driven clusters events. For instance, the PDI data of Sumbekova [35], contained in its entirety $10^6$ samples, which yielded $10^4$ clusters with $N_{PC} \ge 2$, and $10^3$ clusters for $N_{PC} \ge 10$, which undeniably impacts the convergence of any conditioned statistic.
VI. CONCLUSIONS

Unidimensional Voronoï tessellations to analyse preferential concentration provide results consistent with the mean values of cluster size and concentration reported in previous studies under the same experimental conditions. However, quantifying preferential concentration by means of this analysis has some biases that need to be considered.

It is clear that if unidimensional measurements via 1DVOA recover evidence of particle clustering, i.e., $\sigma_V/\sigma_{RPP} > 1$ following the criteria of Monchaux et al.[23], this evidence is reliable as random distributions cannot create spurious traces of clustering given adequate statistical convergence.

On the contrary, if 1DVOA yields $\sigma_V/\sigma_{RPP} \approx 1$, it cannot be directly concluded that preferential concentration is absent within the flow.

The origin of this effect has been tracked to the measuring window size (MWS): if it is too small or too wide, it might not be able to detect preferential concentration present within the turbulent flow by means of the unidimensional Voronoï Analysis (1DVOA). Thus, MWS plays an important role on preferential concentration measurements that needs to be understood and taken into account.

To capture preferential concentration with 1DVOA, MWS should be a fraction of the expected cluster length, $\langle L_C \rangle$. In the experimental conditions reported here, our results suggest the range of MWS that retrieve preferential concentration may have as lower bound the Kolmogorov length scale ($\eta$).

This guidelines are justified by our sensitivity analysis, which shows that the maximum degree of clustering (the magnitude of $\sigma_V$) is found at $MWS_{\star} \sim L/10$, where $L$ is the integral length scale of the carrier phase. In general, this expression seems to follow $MWS_{\star} \sim \langle L_C \rangle$, which combined with $MWS_{\star} \sim L/10$, yields that this window should be of order $\eta$ consistent with previously reported values of $L_C = O(10-20)\eta$ under similar conditions.

Hence, a MWS within the latter range may be able to recover meaningful values of the average cluster linear size $\langle L_C/\eta \rangle$, as well as the average cluster concentration $\langle C_C/C_0 \rangle$.

However, it seems that these thresholds (in general) could depend on specific experimental conditions, or DNS simulation conditions, and therefore an iterative procedure for varying the MWS should be put in place if evidence of preferential concentration is not recovered at the first try by unidimensional analysis. This conclusion is supported by
the 1DVOA performed on the numerical projections of 3D numerical, and 2D experimen-
tal data ($3D_{DNS} \rightarrow 1D_{\perp}$, $2D_{EXP} \rightarrow 1D_{\perp}$, respectively), and previously reported quasi-
unidimensional PDI data.

A pitfall of the 1DVOA is that the ‘raw’ cluster linear size $L_c/\langle L_c \rangle$ and cluster concentra-
tion $C_c/C_0$ probability density functions (PDFs) might not be insightful for characterizing
preferential concentration, as the loss of information inherent to the 1D projections in quasi-
unidimensional experimental methods, e.g., PDI, or optical probes, weakens the correlations
present within data, and yields as consequence the disappearance of a power law behavior in
the right tail of the cluster size PDF.

We develop a simple theoretical model for cluster PDFs from random distributions, which ruled out the existence of a power law within the resulting PDF. However, given the urgent need of characterizing these PDFs, we proposed an alternative filtering approach to disen-tangle randomness from turbulence clustering. This approach conditions the cluster size PDF by the number of particles in a cluster, $N_{PC}$, and is based on the premise that if pref-erential concentration is present in the data, the right tail of the clusters PDF will conserve its power law dependence, as in respective clusters histogram $S_N$. This could have deep implications when analyzing the impact of collective effects on particles settling velocity via quasi-unidimensional measuring techniques.

Our data suggests that the 2D and 3D Voronoi tessellations (2DVOA/3DVOA) are much
more robust to characterize the cluster size PDF than its 1D counterpart. In fact, 2DVOA
is very robust and its biases are minimal, as previously demonstrated [56].

3DVOA, however, presents biases that a recently proposed cutoff criterion [8] does not consider. The results presented for the classical cluster detection algorithm [23] applied to 3D random synthetic data, which does not contain correlations at any scale, yielded the presence a power law on the right tail of the cluster size PDF, as reported elsewhere [15, 16]. In 3D, our filtering approach seems promising, as it allows to clearly determine whether the clusters data is coming from a random, or turbulent origin. In the latter case our results confirm that a powerlaw is preserved for a larger extent when compared to the former case.
## VII. ACKNOWLEDGMENTS

The authors would like also to thank A2PS for all the material, and support provided during the study of the optical probe. We thank Markus Uhlmann and Agatha Chouippe from KIT for fruitful discussions, and Laure Vignal for her help and insight regarding the experimental setup. The numerical (DNS) Lagrangian data set was provided by the TurBase [28] repository under the project EuHIT: a European High performance Infrastructures in Turbulence (Grant agreement # 312778). Our work has been partially supported by the LabEx Tec21 (Investissements dAvenir Grant Agreement No. ANR-11-LABX-0030), and by the ANR Project No. ANR-15-IDEX-02.

## Appendix A: Clusters PDF Model RPP Distribution

In order to model the normalized cluster size PDF for $L_C/\langle L_C \rangle$ coming from randomly generated data (RPP distribution, see table III). The following assumptions were made:

* The normalized cluster size PDF for a cluster population with an arbitrary number of cluster points $L_C/\langle L \rangle|_{N_{PC}}$, where $N_{PC} \ge 2$ is the number of points inside the cluster, is equal to the sum $N_{PC}$ of the independent and identically distributed [57] variables $X_i$, with $i = 1, 2, 3, \dots, N_{PC}$;

$$ \frac{L_C|_{N_{PC}}}{\langle L \rangle} = X_1 + X_2 + X_3 + \dots + X_{N_{PC}} \quad (A1) $$

The random variable $X_i$ is distributed as:

$$ f(\mathcal{V}) = K_f \begin{cases} \mathcal{V}e^{-2\mathcal{V}} & \text{if } 0 < \mathcal{V} < \mathcal{V}_{th} \\ 0 & \text{otherwise} \end{cases} \quad (A2) $$

$$ \langle \mathcal{V} \rangle|_{0<\mathcal{V}<\mathcal{V}_{th}} = \frac{e^{2\mathcal{V}_{th}} - 2\mathcal{V}_{th}^2 - 2\mathcal{V}_{th} - 1}{e^{2\mathcal{V}_{th}} - 2\mathcal{V}_{th} - 1} \quad (A3) $$

$$ \frac{\langle L_C|_{N_{PC}} \rangle}{\langle L \rangle} = N_{PC} \langle \mathcal{V} \rangle|_{0<\mathcal{V}<\mathcal{V}_{th}} \quad (A4) $$

where $\mathcal{V}$ is the normalized cell size, $K_f = 0.25(1 - (2\mathcal{V}_{th} + 1)\exp(-2\mathcal{V}_{th}))$ a constant that accounts for the normalization of the PDF, namely,
$\int_{0}^{\mathcal{V}_{th}} f(\mathcal{V})d\mathcal{V} = 1$. $f$ is the theoretical model PDF for $\mathcal{V}$ proposed by Ferenc et Néda [3] but with its domain bounded by $\mathcal{V}_{th}$, which is the threshold to compute the clusters.

The sum of $N_{PC}$ independent variables has a PDF equal to the convolution of their respective individual PDFs [55], e.g., $N_{PC} = 2$ (the simplest case):

$$f_2(Z) = \int_{-\infty}^{\infty} f_X(X)f_Y(Z-X)dX \qquad (A5)$$

with $Z = X + Y$ (these are dummy random variables, it is equally valid $Z = X_1 + X_2$), considering the support of components distribution, a traditional technique consist in dividing the range of the new random variable in two:

$$f_2(Z) = K_2 \begin{cases} \int_0^Z X e^{-2X} (Z-X) e^{2(X-Z)} dX & 0 < Z \le \mathcal{V}_{th} \\ \int_{Z-\mathcal{V}_{th}}^{\mathcal{V}_{th}} X e^{-2X} (Z-X) e^{2(X-Z)} dX & \mathcal{V}_{th} < Z \le 2\mathcal{V}_{th} \end{cases} \qquad (A6)$$

Where $K_2$ is a normalizing constant. As the original variables had a support from $0 < Z \le \mathcal{V}_{th}$ is straight forward to see that the new support is $0 < Z \le 2\mathcal{V}_{th}$, further simplifying:

$$f_2(Z) = K_2 e^{-2Z} \begin{cases} \int_0^Z XZ - X^2 dX & 0 < Z \le \mathcal{V}_{th} \\ \int_{Z-\mathcal{V}_{th}}^{\mathcal{V}_{th}} XZ - X^2 dX & \mathcal{V}_{th} < Z \le 2\mathcal{V}_{th} \end{cases} \qquad (A7)$$

After carrying on the integration:

$$f_2(Z) = K_2 e^{-2Z} \begin{cases} Z^3 & 0 < Z \le \mathcal{V}_{th} \\ 6Z\mathcal{V}_{th}^2 - Z^3 - 4\mathcal{V}_{th}^3 & \mathcal{V}_{th} < Z \le 2\mathcal{V}_{th} \end{cases} \qquad (A8)$$

The figure 15a illustrates that there is good agreement with the right tail of the distribution between numerically generated data, and the model here proposed. The disagreement in the left tail comes from the nature of this assumption, as we ignore the influence of the neighbours in the construction of the Voronoï cells. However, it seems that this approximation is satisfactory for the right tail of the distribution.

For larger values of $N_{PC}$, the convolution in 'physical' space becomes cumbersome, and thereby, the duality between the convolution and the Fourier transform, i.e., $\mathcal{F}\{F*$
$G\} = \hat{F} \cdot \hat{G}$ is going to be employed to compute the PDF of $L_C/\langle L \rangle|_{N_{PC}}$, for $N_{PC} \ge 3$, whence:

$$f_{N_{PC}} = \mathcal{F}^{-1}\{\prod_{k=1}^{N_{PC}} \hat{F}\} \qquad (A9)$$

being $\hat{F}$ the Fourier transform of the PDF found in equation A2. This approach (see figures 15a-15b) has good agreement regarding the PDF right tail despite the increase in statistical uncertainty due to the decrease in the number of samples for an increasing $N_{PC}$ (see figure 9a)

It is expected that as $N_{PC} \to \infty$ the PDF collapses into a normal/gaussian distribution due to the central limit theorem [55], as seen in figure 15c. The average of this distribution is given by equation A4, e.g., for $N_{PC} = 30$, $\langle L_C \rangle/\langle L \rangle|_{N_{PC}=30} \approx 10$.

* The second assumption deals with the composition of the PDF for clusters having between $2 \le N_{PC} \le N_{PC}^*$, with $N_{PC}^* = 3, 4, \dots, N$ points.

The approach taken was to make a mixture model [47], i.e., the composite PDFs will be a weighted sum of the **normalized PDFs** for each $N_{PC}$. This is written as:

$$f_{2 \le N_{PC} \le N} = \sum_{k=2}^{N} \alpha_{k-2} f_k l \qquad (A10)$$

$$\sum_{k=2}^{N} \alpha_{k-2} = 1 \qquad (A11)$$

$$\int_{0}^{nV_{th}} f_n(\mathcal{Z})d\mathcal{Z} = 1 \qquad (A12)$$

Being able to compute $f_k$, the weights $\alpha_{k-2}$ were modeled following a numerical experiment (see figure 9a), which revealed that in 1D there is a predominantly presence of 2-point clusters, and that in turn this presence halves for 3-points clusters, and halves again for 4-point clusters, and so on. This can also be seen as the probability of finding a n-point cluster. The result is plotted in Figure 8b, and it is repeated here for clusters between $2 \le N_{PC} \le 25$ (see figure 15d), the latter being the larger cluster detected in
the numerical experiment. To change from $L_C/\langle L \rangle$ to $L_C/\langle L_C \rangle$, it is straight forward following the chain rule [55];

$$f_Y(y) = f_X(v(y))|v'(y)| \qquad (A13)$$

which completes the explanation of the model here proposed.

Given these previous results and the central limit theorem, one conjectures that for $L_C/\langle L \rangle \gg 1$ the right tail of the clusters PDF will be within an envelope determined only by the weights $\alpha_i$ distribution, regardless of the presence of turbulence. Hence, it validates our model, as our interest is to examine right tail of such PDF.

The left tail however (as explained in the main text) behave as $f(M_C/\langle M \rangle < 1) \approx (M_C/\langle M \rangle)^\alpha$ behaves as twice this exponent, e.g., $\alpha \approx 2/5/8$ in 1D/2D/3D. The lower exponent than the respective on from our model $\alpha = 3$ is due to our assumptions, as we independently $N_{PC}$ cells PDFs, whereas in reality, due to the tessellation construction the distance between the cells ‘centers’ is **correlated**, and such correlation seems to be particularly important for small, ‘compact’ clusters, as for large clusters the joint PDF can successfully modeled as a product of marginal PDFs, as shown.

# 1. Higher Dimensions

The extension to the previous arguments to higher dimensions has the pitfall of not having analytical PDFs available to do the computations, but rather fits proposed by Ferenc and Néda [3], and therefore, a mismatch is expected. The weights also have a different behavior than their counter part in 1D (see figure 9a). For instance, in 2DVOA the figures 16a-16b reveal that there is an acceptable agreement again with the right tail of the numerical experiment. Proceeding into the compound PDF (equation A10), figure 10a shows that the model captures the changes in the behavior of the numerical PDF (inflection points), however the slope is not completely well captured, and effect that might be due to not enough numerical samples, or the inaccuracy of the fits proposed by Ferenc and Néda [3] (error propagation).
[1] R. A. Shaw, Particle-turbulence interactions in atmospheric clouds, Annual Review of Fluid Mechanics **35**, 183 (2003).

[2] M. L. Larsen, R. A. Shaw, A. B. Kostinski, and S. Glienke, Fine-scale droplet clustering in atmospheric clouds: 3d radial distribution function from airborne digital holography, Physical review letters **121**, 204501 (2018).

[3] J. S. Ferenc and Z. Néda, On the size distribution of Poisson Voronoi cells, Physica A: Statistical Mechanics and its Applications **385**, 518 (2007), arXiv:0406116 [cond-mat].

[4] P. Huck, C. Bateson, R. Volk, A. Cartellier, M. Bourgoin, and A. Aliseda, The role of collective effects on settling velocity enhancement for inertial particles in turbulence, Journal of Fluid Mechanics **846**, 1059 (2018).

[5] W. Yuan, L. Zhao, H. I. Andersson, and J. Deng, Three-dimensional voronoï analysis of preferential concentration of spheroidal particles in wall turbulence, Physics of Fluids **30**, 063304 (2018).

[6] A. Chouippe and M. Uhlmann, On the influence of forced homogeneous-isotropic turbulence on the settling and clustering of finite-size particles, Acta Mechanica, **1** (2019).

[7] S. Sumbekova, A. Cartellier, A. Aliseda, and M. Bourgoin, Preferential concentration of inertial sub-Kolmogorov particles: The roles of mass loading of particles, Stokes numbers, and Reynolds numbers, Physical Review Fluids **2**, 24302 (2017), arXiv:1607.01256.

[8] L. Baker, A. Frankel, A. Mani, and F. Coletti, Coherent clusters of inertial particles in homogeneous turbulence, Journal of Fluid Mechanics **833**, 364 (2017).

[9] R. Monchaux and A. Dejoan, Settling velocity and preferential concentration of heavy particles under two-way coupling effects in homogeneous turbulence, Physical Review Fluids **2**, 104302 (2017).

[10] R. Monchaux, M. Bourgoin, and A. Cartellier, Analyzing preferential concentration and clustering of inertial particles in turbulence, International Journal of Multiphase Flow **40**, 1 (2012).

[11] D. O. Mora, A. Aliseda, A. Cartellier, and M. Obligado, Pitfalls measuring 1d inertial particle clustering, in *iTi Conference on Turbulence* (Springer, 2018) pp. 221–226.
[12] M. Hong, A. Cartellier, and E. J. Hopfinger, Characterization of phase detection optical probes for the measurement of the dispersed phase parameters in sprays, *International Journal of multiphase flow* **30**, 615 (2004).

[13] A. B. Kostinski and R. A. Shaw, Scale-dependent droplet clustering in turbulent clouds, *Journal of fluid mechanics* **434**, 389 (2001).

[14] R. A. Shaw, A. B. Kostinski, and M. L. Larsen, Towards quantifying droplet clustering in clouds, *Quarterly Journal of the Royal Meteorological Society* **128**, 1043 (2002).

[15] M. Uhlmann and T. Doychev, Sedimentation of a dilute suspension of rigid spheres at intermediate galileo numbers: the effect of clustering upon the particle motion, *Journal of Fluid Mechanics* **752**, 310 (2014).

[16] M. Uhlmann and A. Chouippe, Clustering and preferential concentration of finite-size particles in forced homogeneous-isotropic turbulence, *Journal of Fluid Mechanics* **812**, 991 (2017), arXiv:1612.00318.

[17] M. Uhlmann, Voronoï tesselation analysis of sets of randomly placed finite-size spheres, *Physica A: Statistical Mechanics and its Applications*, 124618 (2020).

[18] W. Bachalo and M. Houser, Phase/doppler spray analyzer for simultaneous measurements of drop size and velocity distributions, *Optical Engineering* **23**, 235583 (1984).

[19] A. Jaczewski and S. P. Malinowski, Spatial distribution of cloud droplets in a turbulent cloud-chamber flow, *Quarterly Journal of the Royal Meteorological Society* **131**, 2047 (2005).

[20] R. E. Britter and S. R. Hanna, Flow and Dispersion in Urban Areas, *Annual Review of Fluid Mechanics* **35**, 469 (2003).

[21] S. D. Sabatino, R. Buccolieri, and P. Salizzoni, Recent advancements in numerical modelling of flow and dispersion in urban areas: a short review, *International Journal of Environment and Pollution* **52**, 172 (2013).

[22] M. Obligado, T. Teitelbaum, A. Cartellier, P. Mininni, and M. Bourgoin, Preferential concentration of heavy particles in turbulence, *Journal of Turbulence* **15**, 293 (2014).

[23] R. Monchaux, M. Bourgoin, and A. Cartellier, Preferential concentration of heavy particles: A Voronoï analysis, *Physics of Fluids* **22**, 10.1063/1.3489987 (2010).

[24] L. Mydlarski, A turbulent quarter century of active grids: from Makita (1991) to the present, *Fluid Dynamics Research* **49**, 061401 (2017).
[25] D. O. Mora, E. Muñiz Pladellorens, P. Riera Turró, M. Lagauzere, and M. Obligado, Energy cascades in active-grid-generated turbulent flows, Phys. Rev. Fluids 4, 104601 (2019).

[26] S. Sumbekova, Concentration préférentielle de particules inertielles : la structure et la dynamique de clusters, Ph.D. thesis, Universite Grenoble Alpes (2016).

[27] M. Obligado, A. Cartellier, A. Aliseda, T. Calmant, and N. de Palma, Study on preferential concentration of inertial particles in homogeneous isotropic turbulence via big-data techniques, Physical Review Fluids 5, 024303 (2020).

[28] R. Benzi, L. Biferale, F. Bonaccorso, H. Clercx, A. Corbetta, W. Möbius, F. Toschi, F. Salvadore, C. Cacciari, and G. Erbacci, Turbase: a software platform for research in experimental and numerical fluid dynamics, in *2017 International Conference on High Performance Computing & Simulation (HPCS)* (IEEE, 2017) pp. 51–57.

[29] J. Bec, L. Biferale, A. Lanotte, A. Scagliarini, and F. Toschi, Turbulent pair dispersion of inertial particles, Journal of Fluid Mechanics 645, 497 (2010).

[30] J. Bec, L. Biferale, M. Cencini, A. Lanotte, and F. Toschi, Intermittency in the velocity distribution of heavy particles in turbulence, Journal of Fluid Mechanics 646, 527 (2010).

[31] S. Balachandar and J. K. Eaton, Turbulent Dispersed Multiphase Flow, Annual Review of Fluid Mechanics 42, 111 (2010).

[32] M. Tanemura, Statistical distributions of poisson voronoi cells in two and three dimensions, FORMA-TOKYO- 18, 221 (2003).

[33] V. Ramasubramani, B. D. Dice, E. S. Harper, M. P. Spellings, J. A. Anderson, and S. C. Glotzer, freud: A software suite for high throughput analysis of particle simulation data (2019), arXiv:1906.06317.

[34] A. Cartellier and E. Barrau, Monofiber optical probes for gas detection and gas velocity measurements: conical probes, International Journal of Multiphase Flow 24, 1265 (1998).

[35] S. Sumbekova, A. Aliseda, A. Cartellier, and M. Bourgoin, Clustering and Settling of Inertial Particles in Turbulence, in *Proceedings of the 5th International Conference on Jets, Wakes and Separated Flows (ICJWSF2015)*, Vol. 185 (2016) pp. 475–482.

[36] E. W. Saw, R. A. Shaw, J. P. Salazar, and L. R. Collins, Spatial clustering of poly-disperse inertial particles in turbulence: II. Comparing simulation with experiment, New Journal of Physics 14, 105031 (2012).
[37] F. Toschi and E. Bodenschatz, Lagrangian Properties of Particles in Turbulence, Annual Review of Fluid Mechanics **41**, 375 (2009).

[38] S. Goto and J. C. Vassilicos, Self-similar clustering of inertial particles and zero-acceleration points in fully developed two-dimensional turbulence, Physics of Fluids **18**, 1 (2006).

[39] P. J. Ireland, A. D. Bragg, and L. R. Collins, The effect of reynolds number on inertial particle dynamics in isotropic turbulence. part 1. simulations without gravitational effects, Journal of Fluid Mechanics **796**, 617 (2016).

[40] S. Goto and J. C. Vassilicos, Sweep-stick mechanism of heavy particle clustering in fluid turbulence, Physical Review Letters **100**, 1 (2008).

[41] S. W. Coleman and J. C. Vassilicos, A unified sweep-stick mechanism to explain particle clustering in two- and three-dimensional homogeneous, isotropic turbulence, Physics of Fluids **21**, 1 (2009).

[42] A. D. Bragg, P. J. Ireland, and L. R. Collins, Mechanisms for the clustering of inertial particles in the inertial range of isotropic turbulence, Physical Review E **92**, 023029 (2015).

[43] A. Aliseda, A. Cartellier, F. Hainaux, and J. C. Lasheras, Effect of preferential concentration on the settling velocity of heavy particles in homogeneous isotropic turbulence, Journal of Fluid Mechanics **468**, 77 (2002).

[44] C. P. Bateson and A. Aliseda, Wind tunnel measurements of the preferential concentration of inertial droplets in homogeneous isotropic turbulence, Experiments in fluids **52**, 1373 (2012).

[45] T. Wittemeier and J. S. Shrimpton, Explanation of differences in experimental and computational results for the preferential concentration of inertial particles, Computers & Fluids **173**, 37 (2018).

[46] A. J. Petersen, L. Baker, and F. Coletti, Experimental study of inertial particles clustering and settling in homogeneous turbulence, Journal of Fluid Mechanics **864**, 925 (2019).

[47] S. Frühwirth-Schnatter, *Finite mixture and Markov switching models* (Springer Science & Business Media, 2006).

[48] S. M. Ross, S. M. Ross, S. M. Ross, S. M. Ross, and E.-U. Mathématicien, *A first course in probability* (Prentice Hall Upper Saddle River, NJ, 1998).

[49] D. DeRiggi, A central limit theorem for correlated variables with limited normal or gamma distributions, Communications in Statistics-Theory and Methods **48**, 5213 (2019).
[50] I. A. Ibragimov, A note on the central limit theorems for dependent random variables, Theory of Probability & Its Applications **20**, 135 (1975).

[51] R. C. Bradley Jr, Central limit theorems under weak dependence, Journal of Multivariate Analysis **11**, 1 (1981).

[52] T. Kiang, Random fragmentation in two and three dimensions, Zeitschrift fur Astrophysik **64**, 433 (1966).

[53] R. Zamansky, F. Coletti, M. Massot, and A. Mani, Turbulent thermal convection driven by heated inertial particles, Journal of Fluid Mechanics **809**, 390 (2016).

[54] B. Spinewine, H. Capart, M. Larcher, and Y. Zech, Three-dimensional Voronoï imaging methods for the measurement of near-wall particulate flows, Experiments in Fluids **34**, 227 (2003).

[55] R. V. Hogg, J. McKean, and A. T. Craig, *Introduction to mathematical statistics* (Pearson Education, 2005).

[56] R. Monchaux, Measuring concentration with Voronoi diagrams: The study of possible biases, New Journal of Physics **14**, 095013 (2012).

[57] N. L. Johnson, S. Kotz, and N. Balakrishnan, *Continuous univariate distributions, vol. 2* (John Wiley\& Sons, 1994).
<figure><img src="image_15.png" /><figcaption>FIG. 12: PDFs of characteristic cluster sizes compensated by the respective voids exponents given by the work of Goto and Vassilicos [38] (7/3 and 5/3 for 1D and 2D in the order given), by Uhlmann and collaborators [15, 16] (16/9 for 3D). Filled symbols come from the data in table I, whereas open symbols are data coming from their respective RPP. a) PDFs of $L_C/\langle L_C \rangle$. Filled symbols are experimental data PDI from Sumbekova [35] and the 1D-RPP. b) PDFs of $A_C/\langle A_C \rangle$ for our 2D experimental data base, and the 2D-RPP. c) PDFs of normalized cluster volume $V_C/\langle V_C \rangle$ for the DNS data, and the 3D-RPP.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>FIG. 13: Average cluster concentration $\langle C_C/C_0 \rangle$ vs the measuring window scale for numerical and experimental data. Filled symbols refer to MWS for which $\sigma_V/\sigma_{RPP} > 1$, whereas blank symbols refer to $\sigma_V/\sigma_{RPP} < 1$ (figure 5a). The X refers to the PDI data from Sumbekova[26] $\langle C_C/C_0 \rangle \sim 4.56$. The larger the number of points projected, the closer the cluster concentration to the value of 3.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>FIG. 14: a) PDF of the cluster concentration over global concentration $C_C/C_0$. For experimental sampled data (EXP-2D-AG-B) ($2D_{EXP} \rightarrow 1D_{\perp}$), and different measuring window sizes (MWS). b) PDF of the cluster concentration over global concentration $C_C/C_0$. Solid symbols are experimental PDI data [26], whereas open symbols are from RPP-1D data. The vertical (---) line represents the asymptotic value $C_C/C_0 \approx 3.0$ found in figure 13. It can be seen that large cluster concentrations have higher probability in clustering containing data than in random data.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>FIG. 15: a) PDF of the cluster size $L_C/\langle L \rangle$, where $\langle L \rangle$ is the average Voronoï cell size. The plot reveals that the model (equation A8) has good agreement with the right tail of numerical generated data. The plot also shows the Fourier computed PDF. b) PDF of the cluster size $L_C/\langle L \rangle$. The plot reveals that the model (equation A8) has good agreement with the right tail of numerical generated data. The plot also shows the Fourier computed PDF. c) PDF of the cluster size $L_C/\langle L \rangle$ for $N_{PC} = 30$. This figure illustrates that as $N_{PC}$ increases the PDF tends to a normal distribution as expected by the central limit theorem [55]. d) PDF of 1DVOA for normalized linear cluster size $L_C/\langle L_C \rangle$ for a random uniform distribution, experimental data, and the model proposed. For all clusters having between 2 and 25 points ($N_{PC}$). The model here proposed represents well the right tail of the RPP data, and the marker (□) corresponds to data from [26], which computed the clusters with the condition $N_{PC} \ge 1$. The condition $V \le V_{th} = 0.55$ was employed for clustering computation.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>FIG. 16: a) PDF of the cluster size $A_C(A)$. The plot reveals that the model (equation A8) has good agreement with the right tail of numerical generated data. The plot also shows the Fourier computed PDF. b) PDF of the cluster size $A_C$. The plot reveals that the model (equation A8) has good agreement with the right tail of numerical generated data. The plot also shows the Fourier computed PDF.</figcaption></figure>
open markers (measuring window size) / η

<figure><img src="image_20.png" /></figure>

closed markers (measuring window size) / $\mathcal{L}$