Spatial clustering of polydisperse inertial particles in turbulence: I. Comparing simulation with theory

Particles that are heavy compared to the fluid in which they are embedded (inertial particles) tend to cluster in turbulent flow, with the degree of clustering depending on the particle Stokes number. The phenomenon is relevant to a variety of systems, including atmospheric clouds; in most realistic systems particles have a continuous distribution of sizes and therefore the clustering of ‘polydisperse’ particle populations is of special relevance. In this work a theoretical expression for the radial distribution function (RDF) for mono- and bidisperse inertial particles in the low Stokes number limit (Chun et al 2005 J. Fluid Mech. 536 219–51) is compared with the results of a direct numerical simulation of particle-laden turbulence. The results confirm the power-law form of the RDF for monodisperse particles with St ≲ 0.3. The clustering signature occurs at scales ≲10–30 times the Kolmogorov scale, consistent with a dissipation-scale mechanism. The theory correctly predicts the decorrelation scale below which bidisperse particles cease to cluster because of their distinct inertial response. A ‘saturation’ effect was observed, however, in which the power-law exponent is limited by the least clustered particle population. An expression is presented with which a polydisperse RDF can be obtained from the mono- and bidisperse RDFs and the particle size distribution. The DNS data clearly show that the effect of polydispersity is to diminish clustering, and place a bound on the level of polydispersity required to approximate a monodisperse system; this result is of relevance to experimental studies and realistic systems.

2 scale below which bidisperse particles cease to cluster because of their distinct inertial response. A 'saturation' effect was observed, however, in which the power-law exponent is limited by the least clustered particle population. An expression is presented with which a polydisperse RDF can be obtained from the mono-and bidisperse RDFs and the particle size distribution. The DNS data

Introduction
It is a common observation, for instance when stirring milk in a cup of tea, that turbulence produced by the stirring action causes the initially inhomogeneous field of milk to become homogeneously mixed. In this mixing process, the randomly rotating and stretching fluid motion in turbulence causes the initially simple geometry of the milk field to be distorted and stretched, thus increasing the surface area bordering the 'milky' and clear fluid. When this bordering area grows sufficiently, the rate of molecular diffusion across it dominates the other transport mechanisms, and the mixture becomes homogeneous down to molecular scales; see e.g., [1]. Similar mixing processes are ubiquitous in nature, such as the mixing of smokestack plumes and nutrients in rivers, among many other examples.
One might then expect, for example, that something analogous occurs at the boundary of atmospheric clouds. Clouds mix with the surrounding clear atmosphere, to be sure, but the 'marker' in this case consists of discrete particles rather than a continuous scalar, and this results in a range of interesting possibilities. Let us restrict ourselves to the simplest case of liquid clouds (turbulent air containing water droplets). Observations suggest that the distribution of water droplets is, in fact, not always uniform at small scales ( 10 −2 m) (e.g. [2,3]). The key to understanding this clustering phenomenon is the inertia of the particles. Unlike the 'milk' discussed earlier, which is made of particles with size comparable to molecular scales, water droplets have finite size and have a much higher mass density than the surrounding fluid. Thus while the former is passively advected by the fluid, and its motion at any time mimics that of the 3 local fluid motion, the latter (inertial particle) may decorrelate from the advecting fluid motion due to its inertia and gravitational settling. This inertial de-correlation can lead to the clustering of inertial particles in turbulent flow at fine scales (comparable to the smallest structures of the turbulent motion).
This phenomenon, referred to in the literature as either inertial clustering or preferential concentration, is of course not restricted to clouds, but applies to many particle laden, unsteady flows such as fuel droplets in combustion engines, planetesimal formation [4] and so on. Inertial clustering may be important for understanding the dynamics of such systems. For instance, it has been proposed that inertial clustering may lead to increased collision rate of water droplets in clouds, which may help explain the short time scale of rain formation [5][6][7].
This paper focuses on the fundamental investigation of inertial clustering, its multi-scale nature and how it depends on dynamical parameters of the fluid and particles (St, Re, S g , to be defined later on). Specifically, our attention is drawn not only to the clustering of like-sized particles, but also the clustering of particles of different sizes. In most realistic applications including clouds, particles have a continuous distribution of sizes and therefore these cross correlations are needed to understand how inertial clustering influences the evolution of the size distribution. Our purpose is to present a unified comparison of the theoretical, numerical and laboratory results relevant to the multi-size clustering phenomenon.
Owing to the significant amount and distinct nature of the background material specific to the simulation, on the one hand, and the experiment, on the other hand, this work is presented in two companion papers. Part I (this paper) will focus on the results of direct numerical simulation (DNS) and its comparison with diffusion-drift theory. Part II [8] will focus on our experimental work and a comparison of the experimental result with the same simulation results used here. Taken together, Parts I and II allow us to make an indirect comparison of the results between theory and experiment. These general observations and conclusions will be presented, for the sake of logical consistency, in Part II.
In Part I we proceed in section 2 with a short survey of the relevant theories. In section 3, we describe the simulation methodology. In section 4, we present the simulation results and discuss its comparison with theories. A summary of the work presented in Part I and the conclusions that can be drawn from it independent of the experimental data are presented in section 5.

Theories and quantification of inertial clustering
It has been proposed that inertial clustering of particles in turbulence is the result of particles being centrifuged out of regions of high fluid vorticity (highly rotating) as a result of their inertia and thus preferentially concentrating in the regions of high strain. In fact, evidence of inertial particles preferentially concentrating in regions of low vorticity and high strain is abundant (see, e.g., [9][10][11][12][13][14][15][16][17][18][19][20]). The detailed dependence of clustering on particle inertia and spatial scale, however, is still a subject of investigation. Spatial correlations between particles of different inertia, despite their importance in applications, are even less well understood.
The scale dependence of spatial particle correlations can be quantified through the radial distribution function (RDF), which is the ratio of the probability of finding a particle pair separated by a distance r normalized by the same probability for a randomly distributed mixture: where ψ(r ) is the sum over the number of particles found at distance r from each 'trial particle' (for simplicity the trial particles are taken as the total population of N particles). δV r /V is the ratio of an infinitesimal volume at distance r (i.e. a spherical shell of thickness δr ) to the total sample volume. Other metrics exist for quantification of clustering [21] but the RDF has the advantage of being directly related to the droplet collision rate [22], so we use this here. When a population consists of particles of two types i and j (e.g., two different diameters) the RDF can be generalized to the bidisperse form where now ψ i j (r ) is the number of particles of type j encountered at distance r from the N i trial particles of type i. A derivation of equation (2), along with a discussion of some subtleties associated with finite volume and edge effects, is given in section 3.2 of [23]. We expect the degree of clustering to depend upon the particle Stokes number, which is a dimensionless measure of particle inertia defined as where τ p ≡ ρ p d 2 /(18ν) is the particle stopping time, τ k ≡ (ν/ε) 1/2 is the Kolmogorov time scale, r k ≡ ν 3/4 /ε 1/4 is the Kolmogorov length scale, d is the particle diameter, ε is the turbulent kinetic energy dissipation rate, ρ p and ρ are the particle and fluid densities, respectively, and ν is the kinematic viscosity of the fluid. Realistic systems of particles, such as those occurring in atmospheric clouds or various engineering suspensions, possess a range of particle sizes (and possibly particle densities as well). Once the mono-and bidisperse RDFs are defined for all particle sizes, it is possible to write a size-averaged or polydisperse RDF as where p(St)dSt is the fraction of particles with Stokes numbers between St and St + dSt. The equivalent expression in discrete form is where N i is the number of particles of type i, N ≡ N c i=1 N i is the total number of particles and N c is the number of size categories. It is thus sufficient to analyze the bidisperse RDF in order to predict the RDF for the entire particle population without regard to the specific particle types.
The validity of equations (4) and (5) can be verified by expanding equation (5) using equation (2), yielding The double summation of ψ i j over the two indices is equivalent to not distinguishing particles of different St in the evaluation of the RDF. This is by definition the polydisperse RDF.
In the monodisperse case, studies reveal that clustering is strong even at scales much smaller than the scales of the smallest turbulent eddies (Kolmogorov length scale), growing with the inverse of spatial scales as a power law (e.g. [11]): We shall call c 0 the power-law pre-factor (or simply pre-factor) and c 1 the power-law exponent (or clustering exponent). This observation called into question the completeness of the vortexejection picture of inertial clustering since it is not obvious how eddies with length scales at or above the Kolmogorov length scale could cause clustering at scales far below the Kolmogorov length scale. The power-law scaling of the monodisperse RDF can be explained theoretically as a diffusion-drift process [22,24,25]. In essence, the theories depend on a perturbative expansion about a fluid particle trajectory with the particle Stokes number serving as the small parameter. In the limit St 1, the perturbation result for a monodisperse population of particles predicts that c 1 ∝ St 2 . The theories all assume Stokes flow around each sphere and Re p < 1, where Re p is the Reynolds number of the local flow around the sphere. Other assumptions that go into the theory include: particle size is much smaller than flow structures (d/r k 1); the effect of gravitational settling of particle is negligible, which can be stated in terms of the small gravitational settling parameter, S g ≡ v g /v k , where v g = τ p g is the particle terminal falling speed and v k ≡ (εν) 1/4 is the Kolmogorov velocity scale; and negligible flow modification by the particles, which implies dilute particle loading (see, e.g., [26]). According to [22], sub-Kolmogorov clustering arises because inertial particles on average have a drift relative velocity that tends to bring particles closer together. The steady form of the RDF results from the balance between this inward drift velocity and the dispersing effect of the chaotic fluid motion (turbulent diffusion).
The theoretical description of Chun et al [22] also provides a functional form for the bidisperse RDF g i j : where r c = 5.0|St i − St j |r k is a 'saturation' scale below which the power-law increase of the RDF ceases and the RDF approaches a constant. The parameters c 0 and c 1 have the same interpretation as in the monodisperse case, and in fact it should be noted that equation (8) reduces to equation (7) when r c = 0. Equation (8) is therefore the generalized form of the RDF. Besides the Stokes number dependence of g(r ) via c 1 , as described by equations (7) and (8) with the monodisperse prediction of c 1 ∝ St 2 , several questions have remained largely unexplored. The pre-factor c 0 is important in experiments and applications, for example, because it defines the range of scales at which the onset of inertial clustering occurs. A further complication in most realistic flows is the fact that the clustering resulting from inertial effects appears simultaneously with clustering due to mixing across large-scale inhomogeneities [16,27]. Whether or how the two can be separated has been considered in a previous, related experiment [16], but remains an open question. Finally, the possible Reynolds number dependence of inertial clustering has been a longstanding question (e.g. [28]), and is quite relevant for understanding the role of inertial clustering in geophysical flows (e.g. [29]). This set of papers provides insight into each of these topics.
The work presented here is significantly influenced by the RDF approach and the theoretical interpretation of Chun et al [22]. However, we must point out that alternative theoretical approaches to this problem exist (e.g. [30][31][32][33]). Of potential relevance to the present work is the fractal interpretation of inertial clustering. Following [33], inertial clustering is understood as resulting from the dissipative nature of the dynamic equations that leads to a contraction of the phase space of the system. Since the dissipation becomes strong in the sub-Kolmogorov scales (where fluid strain is strongest), particles will converge to a dynamically evolving attractor (a set of geometrical points). The cluster field will become fractal due to scale-invariant dynamics at these scales. A statistical measure of clustering can be provided via the (fractal) correlation dimension D 2 (see [33]). The correlation dimension and the clustering exponent are related via In this sense our results, cast in the language of c 1 , are also relevant to the fractal approach.

Direct numerical simulation (DNS)
DNS of turbulence containing particles has been the primary tool for the scientific investigation of the dynamics of inertial particles in turbulence (e.g. [11,[33][34][35][36]). A direct comparison of the experimental findings with the available theoretical and simulation studies to date is complicated by the fact that most numerical studies have primarily focused on monodisperse (see e.g., [22,24,37]) and bidisperse (see, e.g., [12,22]) particles, whereas experiments and most practical flows typically involve a polydisperse population of particles. There are indications that the RDF signature of inertial clustering is a strong function of polydispersity [15,38,39], and we further investigate this here.
In order to make more meaningful comparisons, we perform DNS of a polydisperse distribution of particles. This allows us to compare our experimental findings directly with the DNS in two ways: (i) a single mode of particle Stokes number, as occurs in an experiment attempting to achieve a single-Stokes-number condition; and (ii) two finite-width modes of particle Stokes number, as occurs in an experiment attempting to achieve bidisperse conditions. We will refer to the former as 'quasi-monodisperse' and the latter as 'quasi-bidisperse'. We first compare DNS with theory for the mono-and bidisperse settings. The DNS serves as a bridge for an indirect comparison between theory and experiments that will be the focus of Part II of this paper [8].

Simulating the fluid turbulence
We simulate homogeneous and isotropic turbulence inside a three-dimensional cubic domain with periodic boundary conditions. The fluid flow is obtained as a solution to the incompressible Navier-Stokes equations and the corresponding continuity equation: where u(x, t) is the fluid velocity field, P(x, t) is the fluid pressure field and F(x, t) is a fluid forcing that injects kinetic energy to maintain the flow. The equations are solved on a 256 3 grid using the standard pseudospectral technique (involving fast Fourier transforms; for details see [40]). Energy is continuously injected into the first two wavenumbers to maintain a Table 1. Turbulence parameters in arbitrary units (except for the first and the last, which are dimensionless) in the DNS. R λ is the Taylor microscale Reynolds number, ε is the turbulent energy dissipation rate, u is the turbulent root-meansquare velocity, ν is the fluid kinematic viscosity, L = π/(2u 2 ) 15ν/ε is the Taylor microscale, r k is the Kolmogorov length scale, T is the large eddy turnover time, τ k is the Kolmogorov time scale, t is the fluid time step and k max is the maximum resolved wavenumber of the simulation.
turbulent flow that is statistically stationary, homogeneous and isotropic [41]. Table 1 lists the major parameters of the flow. Further details of the numerical method can be found elsewhere (e.g. [15,28]).

Simulating the particles
Particles with 250 discrete Stokes numbers in the range 0.01 St 1.2 with uniform increments of δSt = 0.005 between each category, and with each Stokes number having 8000 particles, are introduced randomly into the flow and advanced in time. Hereafter, we refer to each simulated Stokes number as a Stokes number 'line' (in short, St-line). This unique setup makes it possible to study the effects of polydispersity and also provides the flexibility to match an arbitrary distribution, including the measured distribution from the experiment. It is important to note here that we are making the assumption that δSt = 0.005 is small enough to adequately discretize the continuous particle size distribution. This assumption should be reasonable when the finite range of particle Stokes numbers considered is large compared to this spacing, which is valid for the distributions considered in this study. The particles are initially uniformly distributed throughout the domain. The DNS was run for 18T , where T ≡ L/u is the large eddy turnover time, to allow the particles to equilibrate with the fluid, and then particle statistics were gathered until 60T . The choice of 18T as the equilibrium time is a very conservative one since RDFs are apparently steady at significantly earlier times. Other studies have used as low as 6T under similar conditions (see e.g., [28]).
The governing equations for the particles are where x i and v i are the position and velocity of the ith particle, respectively, while u(x i ) is the undisturbed fluid velocity at the particle position, x i . These equations are a simplified version of the Maxey-Riley equations [42] under the assumption of heavy (ρ p /ρ f 1) and small (d/r k 1) particles. Specifically, under these conditions, one may neglect the added mass, Basset history and Faxen corrections. These assumptions are well satisfied in the context of the present work, since ρ p /ρ f ∼ 1000 for water droplets in air and d/r k ∼ 0.1 for the range of droplet sizes used in the experiments. Regarding the second condition, it is generally accepted that r k under-estimates the smallest scales of fluid velocity gradient by a factor of ∼10 (see, e.g., [16]), which makes this condition less stringent. The simulation also assumes that particles are point-like and neglects reverse coupling of the particles on the flow (both conditions are well approximated by the experiments discussed in Part II). Finally, we neglect gravitational settling.
The DNS work of [34] found no appreciable effect of gravity on the particle concentration statistics for S g ≡ v g /u k < 3, where v g = τ p g is the gravitational settling velocity and u k is the Kolmogorov velocity scale. This is consistent with the experimental findings of Saw et al [16], where S g ranges from O(10 −2 ) to O(1). Further technical details include the use of a two-stage second-order Runge-Kutta method (Heun's method) in integrating equation (11) to obtain particle tracks. Fluid velocities at particle centers are obtained via an eighth-order Lagrangian interpolation scheme similar to that described in [43]. The smallest particles were advanced multiple time steps within each fluid time step in order to achieve the specified error tolerance.

Results and discussions
We begin by showing snapshots of particle fields from the DNS in figure 1. Both panels correspond to the same instant in time and the same spatial sub-domain (a 'thin slice' 3r k thick and 2π × 2π across), with the right panel showing particles of higher Stokes number. Clearly, the particles appear clustered with clustering more pronounced at higher Stokes numbers. Closer observation of the right panel reveals some geometrical structure of the clusters. Some of the clusters appear as thin lines, which suggests that they are two-dimensional sheets in three-dimensional space. Further geometric analysis shall be the topic of future work, but we note that since we can see from equation (9) that the clustering exponent (c 1 ) is related to the fractal dimension of the particle field (D 2 ), the presence of sheet-like clusters implies that c 1 will have values approaching unity (D 2 will have values close to 2, the value for a perfect two-dimensional object (e.g. [44])). We shall see in section 4.1 that this is indeed the case (cf figure 3). Throughout our studies, c 1 has never been found to exceed 0.8; it is interesting to consider whether this tentative upper bound is related theoretically to the sheet-like clusters that are observed. Figure 2 shows g(r ) for quasi-monodisperse particles, calculated using three-dimensional particle positions, for various Stokes numbers. Inertial clustering is found to occur for r/r k 20 and increases with St but starts to show saturation at St ∼ 1. These results are 'quasi'monodisperse in the sense that particles from more that one St line were used in the calculation of each g(r ) in order to achieve sufficient statistical convergence. Specifically, for each St level studied we included particles of five St-lines centered on the St of interest; thus in effect we are considering particles with St ± 0.01 (recall that the line spacing is δSt = 0.005 and there are 8000 particles per line). The detailed study (see section 4.2) confirms that the effect of this degree of St averaging results in only a few per cent reduction in the power-law exponent of g(r ) at the length scale resolved here (r/r k 0.1). This reduction is negligible compared with the level of statistical fluctuation in these results. We are thus confident that these results are relevant for monodisperse particles, although strictly speaking they represent lower bounds.

Comparing DNS results with theory: monodisperse
Next, we compare the power-law exponents (cf equation (7)) of these RDFs with the theoretical prediction of [22]. This is shown in figure 3, where the solid parabolic curve is from where ε(t) is the instantaneous kinetic energy dissipation rate, ζ (t) is the instantaneous enstrophy (square of vorticity) times the kinematic viscosity, σ X is the standard deviation of the variable X normalized by its mean X , ρ X Y is the cross-correlation coefficient and T X Y is the correlation time normalized by the Kolmogorov time scale (for details see [22]). In the same work, using semi-empirical inputs from DNS performed at R λ = 47.1, a value of 6.6 was obtained as the coefficient in front of this St 2 -law. The figure clearly shows good numerical agreement between the two results for St 0.3 even though it is difficult to confirm the St power. For example, the precision of the data is not sufficient to rule out a linear scaling. We partially address this point by performing a linear fit (method: error accounting least-square fitting using modelŷ = mx + y o ) on the DNS data in the range St 0.4 and found that y o = −0.09 ± 0.04. However, since a negative intercept is unphysical we conclude that our result supports positive curvature for small St, which is consistent with the St 2 scaling predicted by Chun et al [22].
The theory for finite St is still open [25], but the peaking and subsequent decrease of c 1 with increasing St, starting around St = 0.5, are likely related to the physical argument that heavier particles, with inertial response times larger than the correlation time scale of the smallest (dissipative) eddies (typical ∼τ k ), do not respond well to the centrifugal effects of these eddies, thus resulting in weaker inertial clustering. Other works (see, e.g., [11,33]) have shown that particles become less clustered at larger St.
One point worth addressing here is the fact that the theoretical results rely semi-empirically on a DNS that was run at R λ = 47.1, lower than in the present simulation (R λ = 143). The good  [45,46]. This simulation was performed at R λ = 83).
numerical agreement suggests that inertial clustering may not be a strong function of Reynolds number. This matter will be addressed more closely in Part II [8]. On the other hand, if one considers equation (12), it is not obvious that the complicated coefficient would have only a weak Reynolds dependence. Currently, there is no clear explanation for this [28].
Alternatively, Chun et al [22] also provide a more general form for c 1 , still in the limit of St 1: where S 2 ≡ S l j S l j and R 2 ≡ R l j R l j (summation over repeated indices implied) are the second invariant of the rate of strain and rate of rotation tensors, respectively (i.e. , and · p implies averaging over an ensemble of inertial particles. Using this formulation and their DNS results for S 2 p and R 2 p , the resulting function c 1 (St) is still a curve but the region of St 2 scaling was limited to St 0.05. At higher St this curve was found to be below the theoretical curve of 6.6St 2 . The question of Reynolds dependence raised above is equally valid for this formulation. Our findings here imply that if the theory were correct, these statistics should be weak functions of Reynolds number (an observation further strengthened in Part II [8] with the help of experimental data).
Finally, figure 4 shows a comparison of independent DNS results on the trend of the clustering exponent, c 1 , as a function of Stokes number, St. The agreement between these results while not surprising is by no means trivial, since different researchers use slightly different simulation implementations (turbulence forcing, fluid velocity interpolation, data reduction). The close agreement seen here increases the confidence in the data.
We conclude that there is good agreement between theory and DNS for monodisperse particles at small Stokes numbers, up to St 0.3. The disagreement at larger St is not surprising since the perturbation theories are only valid in the limit of St 1.  Figure 5 shows how the RDF changes from the case of (quasi) monodisperse to polydisperse. In particular, the slope (power-law exponent) decreases with increasing polydispersity, and a plateau region appears at small r , which grows with increasing polydispersity. The power-law model for the RDF is modified by polydispersity. The appearance of the plateau at small scales in the RDF signifies that one will not see clustering at or below these scales. A physical interpretation of this saturation is that, on average, particles of different St tend to cluster at slightly offset locations. This is hardly surprising and reflects the fact that particle trajectories are deterministic functions of their inertia. Saturation of the RDF with decreasing scale can also be interpreted as resulting from the differential response of the two particles to fluid accelerations [22]. The diffusion resulting from this differential response is independent of particle separation distance, and therefore leads to the form of equation (8), with r c being the scale at which this acceleration diffusion is of the same order as the inward drift.

Effects of polydispersity on the radial distribution function (RDF)
In many situations, one may wish to know how polydispersity affects c 1 . Figure 6 shows c 1 as a function of the width of the Stokes number distribution, St. We note that the clustering exponent c 1 decreases linearly with increasing Stokes number width St, up to approximately 0.1. This results in a simple rule of thumb that the Stokes number width at which c 1 is diminished by 10%, St 10% , is roughly St mid /5 in the range of Stokes numbers studied.
This monotonic reduction of the power-law exponent with increasing St further motivates the need to consider polydisperse particle distributions in realistic applications. Equation (8) relates the polydisperse RDF, to the mono-and bidisperse RDFs, and in the next section the bidisperse theory is compared to the DNS results. Figure 7 illustrates the qualitative behavior of the quasi-bidisperse RDF. Here the value of St 1 is fixed while St 2 is varied for each RDF in order to show the effect of bidispersity (RDFs are invariant with respect to the exchange of indices, i.e. g 12 (r ) = g 21 (r )). The general trend reveals a power-law tapering off to a plateau, with the transition shifting to larger scales as the difference between St 1 and St 2 is increased. This behavior is similar to that observed for polydisperse particles in section 4.2, which is not surprising since the two cases are connected through equation (4). However, as can be seen from comparing the two panels of figure 7, the behavior of g 12 is asymmetric depending on whether St 2 is decreased or increased relative to the fixed St 1 . When St 2 is decreased with respect to St 1 the slope of g 12 (at the inflection point) diminishes and a plateau starts to appear. On the other hand, when St 2 is increased from St 1 , the slope of g 12 remains nearly constant. The magnitude of g 12 at r above the scale corresponding to the inflection point even becomes slightly larger than the monodisperse case. We will return to this point later.

Comparing DNS with theory: bidisperse RDF
The general shape of each of the RDFs in figure 7 is consistent with the saturated power-law model described by equation (8) [22]. As noted earlier, in the theory r c is the scale below which flattening of g 12 occurs, while c 0 and c 1 have the same interpretation as in the monodisperse case. Hence, equation (8) reduces to the monodisperse RDF (equation (7)) when r c = 0 and is therefore a generalization that encompasses both cases. The theory also predicts (in the limit of St 1, and with semi-empirical input from DNS at R λ = 47) that This provides a theoretical basis for interpretation and comparison with the DNS results reported here: a nonlinear regression analysis of the RDFs to equation (8) is performed in the range of r/r k 0.2 -6, resulting in empirical values for the saturation scale r c and the power-law slope c 1 . We note that by inspection, in all cases, equation (8) fits the data reasonably well. Figure 8 shows a comparison of theory (equation (14)) and DNS for the trend of the flattening scale, r c . The agreement is surprisingly good even at higher Stokes numbers (e.g. St 0.4), especially considering that the theory was developed in the limit of St 1. It is clear that r c behaves symmetrically and linearly with respect to the origin of St 2 − St 1 (hereafter denoted as St 21 ). This implies that r c only depends on the absolute value of St 21 . This trend of g 12 can be understood as follows. Particles of different St cluster with 'sharpness' characteristic of their St. The 'clustering sites' of particles with different Stokes numbers are still in the vicinity of each other (e.g. in the regions of high low vorticity); however, they are slightly offset and r c reflects an isotropic average of these offsets. These offsets presumably have isotropic statistics over the particle field, at least when there is no effect of gravity or other non-isotropic factors.
In order to compare the theoretical prediction represented by equation (15) with the DNS, one needs to determine the power-law exponent, c 1 , which is only recovered at sufficiently small scales (r < r k ). In reality, r c grows quickly with | St 21 | and its flattening effect approaches r k for | St 21 | as small as 0.03, preventing the clear observation of a power-law exponent c 1 in the bidisperse RDFs for | St 21 | 0.03 (see equation (14) and figure 8). Since much of our data fall into this category, we are unable to make a direct comparison with the theory of equation (15). However, as is evident in figure 7 and the top-left panel of figure 9, there exists a range of scales (r/r k 1-6) in between r k and the large-scale flattening of the RDF, where a power-lawlike region is still observed. We present here an analysis of the RDFs in this range. In order to avoid confusion, we define a new power-law exponent, c 1 , which corresponds to the slope of the RDF (in log-log axes) in this intermediate range. Defined this way, c 1 is empirically obtained when the RDFs are fitted by equation (8) in the range of r/r k 0.2-6 (we note that in the monodisperse cases, the magnitudes of c 1 are nominally not less than 80% of their asymptotic counterparts, c 1 , which are obtained by fitting over the full range of resolved scales).
We now consider how well equation (15) describes this intermediate range. This range can be thought of as the early part of the transition from dissipation scales to inertial scales. Figure 9 shows that there is, in fact, a rather marked departure from the theoretical expectation. In the top-left panel it can be seen that with one Stokes number fixed (St 1 ) and the other varied (St 2 ), the slope of the RDF does not increase monotonically over the full range as predicted by equation (15). Instead, the power-law exponent increases for St 2 < St 1 and then saturates for St 2 St 1 (we emphasize that g 12 (r ) ≡ g 21 (r ) by definition). This trend is observed for all Stokes numbers considered in this study. The top-right panel of figure 9 further illustrates this point. Here the power-law exponent of the RDF (c 1 ) is plotted versus St 2 , for various conditions on St 1 . The black plot, as the control, shows the results from the monodisperse case (St 1 = St 2 as St 2 varies). The red, purple and cyan plots are for fixed St 1 with the corresponding values of 0.14, 0.4, 0.7, respectively. In each of the cases, the bidisperse c 1 follows closely the trend of the monodisperse results for St 2 < St 1 , but is almost constant when St 2 > St 1 . As can be anticipated from these results already, compensating for the bidisperse RDFs by equation (15) leads to a poor collapse (not shown). To what extent this saturation behavior is a shortcoming of the theory in its valid dissipation-range application, as opposed to a deviation relevant for the transition to the inertialrange scaling, remains an open question. Generally, however, the unexpected behavior of the bidisperse c 1 can be understood as a 'saturation' effect, where c 1 is always strongly limited by the least clustered set of particles. We speculate that while particles of types 1 and 2 are each correlated to a degree that is controlled by their respective Stokes numbers, their cross correlation, within this intermediate range of scales, is controlled by the least correlated of the two particles. To make a better prediction that reflects the observed saturation effect, c 1 can be approximated by the smaller of the two monodisperse c 1 , each corresponding to one of St 1 and St 2 : where c 1 (St i ) refers to the power-law exponents of the RDF of monodisperse particles with the corresponding Stokes number. Indeed, reasonable collapse of the RDFs within the range r/r k 1-6 was obtained when g 12 (r ) is compensated by equation (16), as shown in the bottom panels of figure 9. We conclude that, in the regime of finite Stokes number difference we have studied (| St| of the order of 0.01-0.1), the bidisperse theory can be extended to capture the behavior of the clustering saturation scale, r c , but does not describe the full complexity of the behavior of the intermediate power-law exponent c 1 . To account for the observed saturation effect, c 1 can be approximated by the minimum value of monodisperse c 1 , for the two Stokes numbers considered. It should be understood that equation (16) leaves open the exact dependence of c 1 on St. It is therefore adequate to extend the Chun et al [22] theory, allowing it to make reasonably accurate predictions of the polydisperse RDF, even with finite Stokes number differences (cf equation (4)).

Summary and conclusions
Part I of this study has focused on comparing the theory of Chun et al [22] with DNS for monoand bidisperse particle sizes. Through the formula shown in equation (8), these results can be used to calculate the RDF for an arbitrary particle size distribution (polydispersity).
Our results confirm that the RDF arising from inertial clustering is power-law-like and extends to the limit of very small length scales (r r k ) for monodisperse particles. Quantitative agreement was found between DNS and the monodisperse theory for St 0.3. For non-zero St, the RDF becomes larger than unity at separation distances r < 10 -30 r k . This finding supports the understanding that inertial clustering is driven by the dissipation-scale fluid motions and is consistent with the general consensus that the turbulent enstrophy reaches a maximum around 10-20 Kolmogorov lengths (r k ). In a practical sense, this allows for some relaxation of the resolution requirements in an experiment designed to investigate the RDF of inertial particles.
For the case of bidisperse particles we found that the theory can be extended to the regime of finite Stokes number differences (section 4.3). Specifically, when comparing the bidisperse RDFs in this regime, we found that the theory was able to predict the trend of the flattening scale r c very well, but fails to account for the 'saturation' effect found in the trend of c 1 (the experimental data in Part II confirm this effect). The observed RDFs suggest that the bidisperse c 1 depends only on the most weakly correlated particles (see equation (16)).
Besides this, we have also studied the effects of polydispersity on inertial clustering, in view of its importance for many realistic applications (e.g. clouds) and interpretation of experimental data. In section 4.2 we showed, using the DNS data, that the effect of polydispersity is to diminish clustering. Specifically, when polydispersity is increased, the RDF slope is lowered and a plateau region appears below a small length scale that grows with polydispersity. We further found that the reduction of the slope (c 1 ) becomes significant ( 10%) when St St mid /5. We also provided a mathematical relation that allows one to predict the polydisperse RDF given the mono-and bidisperse RDFs and the particle size distribution (equation (4)).