Clustering of settling charged particles in turbulence: theory and experiments

Atmospheric clouds, electrosprays and protoplanetary nebula (dusty plasma) contain electrically charged particles embedded in turbulent flows, often under the influence of an externally imposed, approximately uniform gravitational or electric force. We have developed a theoretical description of the dynamics of such systems of charged, sedimenting particles in turbulence, allowing radial distribution functions (RDFs) to be predicted for both monodisperse and bidisperse particle size distributions. The governing parameters are the particle Stokes number (particle inertial time scale relative to turbulence dissipation time scale), the Coulomb-turbulence parameter (ratio of Coulomb ‘terminal’ speed to the turbulence dissipation velocity scale) and the settling parameter (the ratio of the gravitational terminal speed to the turbulence dissipation velocity scale). The theory is compared to measured RDFs for water particles in homogeneous, isotropic air turbulence. The RDFs are obtained from particle positions measured in three dimensions using digital holography. The measurements verify the general theoretical expression, consisting of a power law increase in particle clustering due to particle response to dissipative turbulent eddies, modulated by an exponential electrostatic interaction term. Both terms are modified as a result of the gravitational diffusion-like term, and the role of ‘gravity’ is explored by imposing a macroscopic uniform electric field to create an enhanced, effective gravity.


Introduction
Collections of macroscopic electrically charged particles are found in atmospheric clouds [1,2], colloids [3,4], dusty plasmas [5,6], protoplanetary nebula [7] and electrosprays [8], and the electrostatic interaction forces can strongly influence the dynamics and spatial correlations between particles [3]- [6], [9,10]. In many of these systems, the relative role of 'settling' due to a gravitational force (or other constant, externally imposed forces) is also of importance to the particle dynamics and spatial distribution [11,12]. We ask the question: How does the presence of unipolar charge (one sign) and gravitational settling affect the spatial arrangement or clustering of particles in turbulence? The behaviour of such particles is relevant particularly to the atmosphere, where cloud particles are both charged [2] and gravitationally settled within a turbulent medium. The spatial arrangement of charged, settling particles is of possible relevance because spatial clustering is one aspect of the collision-coalescence process responsible for rain initiation [11,13,14].
Some comments about what is already known about the spatial distribution of particles in turbulence will help us to place the present work in the proper context. In the absence of charge and gravitational forces, inertial particles have been observed to form spatial 3 clusters by avoiding regions in the fluid with high vorticity and getting accumulated in regions of high strain [15,16]. Alternative interpretations for the clustering are the sweep-stick mechanism in which particles tend to congregate and are advected with lowacceleration points in the fluid [17], and a fractal geometry interpretation based on the dissipative particle dynamics [18]. This clustering can be quantified by the radial distribution function (RDF), which can be defined for particles of the same size (monodisperse) or for particles of two different sizes (bidisperse). The monodisperse RDF for electrically neutral particles in the absence of gravity increases as a power law with decreasing spatial scale, with a power-law exponent that is a function of the particle inertia [13,19,20]. The bidisperse RDF between two particle sizes (bidisperse RDF), again electrically neutral in the absence of gravity, has been shown to increase first as a power law and then plateau to a constant value with decreasing spatial scale [21,22]. The addition of gravitational settling has been found in computational studies to cause the bidisperse RDF for neutral particles to plateau at a larger length scale than it would without settling [12]. Finally, the monodisperse RDF for charged particles has been experimentally shown to behave first as an increasing power law, as with the uncharged case, and then to reach a peak followed by a rapid decay with decreasing spatial scale [10].
The increase in the RDF with decreasing scale can be described theoretically as the result of a time-averaged drift between particles as a result of their finite inertia [13,21]. For example, the drift-diffusion (Fokker-Planck) model of Chun et al [21] yields functional forms of the RDF consistent with the behaviour of monodisperse and bidisperse uncharged particles [21], and of monodisperse charged particles [10]. In this work, we further extend the model of Chun et al [21] to include both gravitational settling and electrostatic interactions for a bidisperse distribution of particles. We then evaluate the theory by comparing its predictions to data obtained from an experiment in which inertial particles in turbulence are charged and placed in several 'effective' gravitational fields (combination of the earth's gravity and an externally imposed electric field).
We now formally pose the problem of the behaviour of a pair of charged particles settling in a turbulent flow. Consider spherical particles of two different sizes (superscripts in square brackets or subscripts without brackets of 1 and 2 will be used to distinguish them) that have unipolar electrical charge q 1 and q 2 , and that are in the presence of a uniform gravitational field of magnitude g. The particles reside in a turbulent flow with Kolmogorov length and timescales η and τ η , with fluid density ρ f . Let d be the particle diameter, ρ p be the particle density, m be the particle mass, q be the particle charge and τ p be the particle inertial response time. In the Stokes drag regime, τ p = d 2 ρ p /18νρ f , where ν is the kinematic viscosity of the fluid. The relative strengths of particle inertia, gravitational settling and electrostatic interactions compared to turbulence on the motions of the particles can be quantified by the Stokes number, settling parameter and a new Coulomb-turbulence number, respectively, where u g = τ p g is the gravitational terminal settling speed, u q is the product of τ p and the Coulomb acceleration summed for two charged particles separated by distance η, β ≡ τ p /m is particle mobility [23] and u η is the Kolmogorov velocity scale u η = η/τ η . The Stokes number St is the ratio between the particle's inertial response time and the Kolmogorov time scale. The settling parameter Sg is the ratio of the particle's terminal settling velocity (it is assumed that ρ p ρ f ) to the Kolmogorov velocity. The Coulomb-turbulence number is the ratio between the 'terminal' speed resulting from the Coulomb force between two charged particles separated by the Kolmogorov length scale, and the Kolmogorov velocity. Use of the word 'terminal' in this context is intended to draw attention to the analogy with gravitational settling since both come from time-independent forces. The fact that the Coulomb force is spatially dependent of course means that the analogy is imperfect, but for small τ p we can expect that the fluid-particle relative speed can be described by the 'terminal' or steady state value to within reasonable approximation. The general problem as defined above is highly complex and we therefore choose limits relevant to practical applications, which allow for significant simplification. The assumptions are (i) isotropic homogeneous turbulence in an incompressible Newtonian fluid, (ii) d η, (iii) ρ p ρ f , (iv) Re particle 1 (Stokes flow), (v) dilute limit (λ > η, where λ is the mean particle separation, and negligible multi-particle electrostatic interactions), (vi) modest gravity (Sg 1), (vii) modest electrostatic force over most of the dissipation range (Ct 12 1), (viii) weak particle inertia (St 1).
Several comments can be made about these limits. The settling parameter can be written as where a η is the Kolmogorov acceleration scale [24], and it follows therefore that limits (vi) and (viii) are related. For example, limit (vi) can be written as St a η /g, so under common conditions such as for water droplets in turbulent air, condition (viii) is automatically satisfied if (vi) is satisfied. In limits (vi) and (vii), the use of 1 implies that the dimensionless quantity can safely reach the order of unity, but the theory will break down in the limit 1. The primary contributions made here, beyond already published work [10], are: the drift-diffusion theory is stated more completely in its application to charged particles, with a clearer discussion of relevant assumptions; the theory is extended to include bidisperse particle size distributions; for bidisperse particles the influence of gravitational settling is taken into account; further experimental results are presented with better resolved RDFs (a larger number of particle pairs sampled); and the variation of effective gravity through imposition of a largescale electric field. This paper is organized as follows. The equation of motion for the relative velocity between particles is developed in section 2. The drift-diffusion model for the RDF is introduced and the contributions of inertia, settling and electrostatic interactions to the drift and diffusivity are developed in section 3. Then, the functional forms of the monodisperse and bidisperse uncharged and monodisperse charged RDFs are reviewed and the bidisperse charged RDF derived in section 4. In section 5, the experiment designed to provide monodisperse and bidisperse RDFs of settling charged particles in turbulence is described. The derived functional 5 forms of the RDF (from section 4) are then compared to the experimental RDFs in section 6. The main results and remaining questions for future investigation are presented in section 7.

Equation of relative motion for a pair of particles
Consider a pair of particles numbered 1 and 2 that might be electrically charged in a turbulent flow and are in a gravitational field. Take a reference frame centred on particle 1 and let the position of particle 2 relative to the first be r. Then, using the assumptions of Stokes drag (Stokes flow), the dilute particle limit (no three-body electrostatic interactions) and no Debye shielding, the equations of motion for the two particles are where x [ j] and v [ j] are the position and velocity of the jth particle, u is the fluid velocity, and k = 1/4π 0 , with 0 being the permittivity of free space. The relative velocity w ≡ v [2] − v [1] = dr/dt of the particles is then The gravitational term cancels, but settling is implicitly contained in v [1] i . We make this explicit by substituting in the equation of motion for particle 1 from equation (4),

Simplifications and assumptions
If r is small compared to the dissipation eddy scale of ∼10 η to ∼20 η [25], the flow is smooth, allowing for Taylor expansions in the velocity field. Then, to first order, the first term on the right side of equation (6) becomes ( i j r j − w i )/τ [2] p , where i j ≡ ∂u i /∂ x j is the fluid velocity gradient tensor evaluated at x [1] . If the particle Stokes numbers St are small and there are no other very strong forces on the particles, the particles nearly follow the flow as tracers.

6
These conditions of St 1, and modest gravitational and electrostatic forces, Sg 1 and Ct 1, imply that there is no significant coupling with the fluid forces, and therefore the resulting particle velocities can be treated independently. Specifically, gravitational settling results to first order in the addition of a terminal speed gτ p to the relevant particle velocity component. Likewise, electrostatic interaction results in the addition of a Coulomb terminal speed τ p (kq 1 q 2 /mr 2 )(r i /r ) to the particle velocity. The measurements of Lu et al [10] were consistent with the assumption of independent treatment of turbulent and electrostatic forces.
The precise parameter range (St, Sg, Ct) over which these assumptions are valid is left to future investigations. However, for St that are not small, very large settling terminal velocities, and/or r small enough that the Coulomb terminal velocity is very large, these assumptions would clearly be invalid. We now turn to the problem of using this differential equation for the relative velocity between two particles to obtain the RDF.

The radial distribution function (RDF) balance equation
Spatial correlations among particle positions are quantified here via the RDF g 12 (r ), essentially the probability that any two particles of type 1 and 2 are separated by length r , relative to the probability of a perfectly random arrangement. Particle types 1 and 2 here denote two different τ p ; note that the monodisperse RDF is obtained when 1 and 2 are identical. Using a Fokker-Planck theoretical framework, Chun et al [21] showed that in steady state the g 12 (r ) results from a balance between a mean, inward relative drift velocity between particles and turbulent diffusion that tends to even out the spatial distribution. The balance equation is 0 = − w p g 12 + D 12 dg 12 dr , where · denotes averages over particle paths and D 12 is a particle turbulent diffusivity resulting from turbulent dispersion and relative motion due to different particle accelerations.

Drift velocity
If the Stokes numbers St are small, the particles nearly follow the flow as tracers, allowing r i and w i to be expanded in powers of St 2 . This results in a mean drift velocity due solely to the first term (inertia term) of equation (6), in the absence of settling and electrostatic interactions [21]. This inertial drift is where S and R are the dimensionless second invariants of the strain and rotation rate tensors, which are related to the velocity gradient tensor as τ 2 η l j jl = S 2 − R 2 . Since inertial particles tend to favour strain over rotation, the drift is negative, meaning that particles tend to drift towards each other [21].
This expression can be extended to include the influence of gravitational and electrostatic forces. For particles of different τ p , gravitational settling causes particles to approach or recede from each other depending on whether the particle with the larger τ p is above or below the other. Both geometrical arrangements are equally likely, so gravitational settling does not contribute to 7 the mean drift velocity. Assuming unipolar charge, the electrostatic interaction adds an outward, relative Coulomb terminal velocity of magnitude (β 1 + β 2 )(kq 1 q 2 /r 2 ). Hence, the electrostatic interaction contributes to the mean drift velocity, which becomes [10]

Diffusivity
Chun et al [21] argued that the particle diffusivity in equation (7) can be written as where the first term is a 'nonlocal' turbulent diffusivity with a dimensionless Batchelortype dispersion constant B nl . Modelling turbulence as a sequence of random extensional and compressional flows along one axis (rotation assumed to have a negligible impact), Chun et al [21] obtained B nl = 0.0926 to leading order. We will use this value in the analysis presented in this paper. The second term D is a diffusivity that comes from temporal correlation of the relative motion of the pair of particles as a result of particle accelerations. More precisely, it comes from an autocorrelation integral of w, or the radial-radial part of the diffusivity tensor where i and j are tensor indices as opposed to denoting particles 1 and 2. For small St without settling or electrostatic interactions, the relative velocity between the two particles is approximately where a i is the fluid acceleration. The mean drift velocity (given by equation (9)) makes a negligible contribution. The second equality follows from the first-order approximation a i (x [1] ) = a i (x [2] ) = a i in the dissipation range [21]. It has been demonstrated that when Sg 1, there is strong coupling between the gravitational settling velocity and the turbulence-induced velocity fluctuations [12,26]. By analogy, we expect that the limit Ct 12 1 will result in strong coupling between the Coulomb terminal velocity and the turbulence-induced velocity fluctuations. The Coulomb-turbulence number is defined relative to the Kolmogorov length scale, but one can also consider a similarly defined function of inter-particle distance r : given the inverse-square dependence of the Coulomb force and the smooth velocity field within the dissipation range, such that u(r ) ≈ u η (r/η), it follows that Ct 12 (r ) ∝ r −3 . We may expect, therefore, that there will be a separation distance sufficiently small that nonlinear coupling of the turbulence and Coulomb velocities will be present. In the absence of other forces, however, we expect that very few particles are able to penetrate this Coulomb-dominated region. The implication is that the condition Ct 12 1, interpreted as Coulomb-turbulence number less than or of the order of unity, is automatically satisfied. In the context of this work then, we will assume that r is sufficiently large that the Coulomb terminal velocity is of the order of the inertial drift velocity and therefore negligible in the consideration of the diffusivity. This assumption has been checked empirically already [10] and will be further evaluated in the experimental section here, but it nevertheless remains an area of interest for future investigation.
Gravity does not fluctuate, however, and affects each particle by giving it, to first order, an average velocity in the downward direction equal to its terminal settling velocity. For particles of different sizes, this adds to the relative velocity resulting from differential response to fluid acceleration. In the limit of dilute particle loading, and in a still fluid, the particle terminal speed is gτ p and the difference in terminal speeds of two particles is (τ p )g. Due to the effects of turbulence and mass loading [27], the terminal speed of a particle becomes Cgτ p , where C presumably is a function of the turbulence parameters, the particle number density and possibly τ p . This makes the difference in terminal speeds (Cτ p )g. If C is significantly different from unity, one can utilize an effective g e ≡ g (Cτ p )/ (τ p ) that makes the difference in terminal speeds (τ p )g e . Including the difference in terminal velocities, the total relative velocity between the two particles where the sign of the gravitational term reflects which particle is above the other. The diffusivity then becomes noting that cross terms are zero due to isotropy (i.e. the fluid acceleration is not correlated with the direction of g i ). Furthermore, because g i is fixed, we can write where a is any component of the fluid acceleration and the factor of 1/3 comes from a surface integral over all polar and azimuthal angles of the squared radial component of the gravitational acceleration vector. The first term D a is the contribution to diffusivity due to different responses to the fluid acceleration, which Chun et al [21] showed can be approximated as where a 0 is a constant associated with the acceleration variance a 2 = a 0 a 2 η , a η is the Kolmogorov acceleration η/τ 2 η , and τ a is the acceleration correlation time scale. They found that a 0 has a slight Reynolds number dependence, but has a value of a 0 ≈ 1.5 for direct numerical simulations with Re λ of 47.1 and 57.3 [21]. For the experiments described in this paper in section 5, which have an Re λ of 80, we will use the same value. From direct numerical simulation results, Chun et al [21] obtained a fluid acceleration correlation time of τ a = 1.5τ η .
Why should the gravitational acceleration, which is constant, appear as a diffusion term, as in D g ? Clearly the integral, as written with a lower limit of −∞, diverges, and in fact we would expect gravity in isolation to lead to ballistic rather than diffusive particle separation. An important observation, however, is that our system consists of particles embedded in a turbulent flow containing a continuous and large range of eddy sizes and time scales. It is significant that we have confined ourselves to the limit Sg 1, with Sg defined with respect to the Kolmogorov velocity u η . As particle separation exceeds one Kolmogorov eddy size η, the larger eddies with velocity scales much greater than gravitational terminal speed, will cause explosive particle dispersion and resultant randomization of particle positions, as well as randomization of particle relative velocity. In short, one could define an Sg r for any eddy in the inertial subrange with characteristic size r η, resulting in Sg r 1, meaning negligible gravitational settling effects at that scale. The implication, therefore, is that the terminal fall speeds are much smaller than most of the eddy velocities in the inertial range. The droplets are still being randomized, therefore, such that on a short (i.e. Kolmogorov) time scale, droplets can be thought of as having a correlated relative drift, but on longer time scales they may be exchanged, with gravity tending to cause the droplets to approach each other or to recede from each other with equal probability.
Given this physical picture, it is clearly unrealistic to integrate from −∞ to t in D g since the relative motion of particles of 1 and 2 do not have an infinite correlation time. Instead, there is a finite correlation time related to the particle dispersion caused by eddies with r η (i.e. distribution of w increases in width with increasing r ). The contribution to the diffusivity by gravity can be approximated as where τ g is the gravitational correlation time, and finally It is important to note that D is not a function of r , which will be important when the balance equation for the RDF (equation (7)) is solved in section 4. Because we are working in the range Sg 1, we can expect relatively weak decoupling of particles from the fluid due to gravity. For small Sg, it is reasonable to assume that τ g = τ a , since decorrelation is due primarily to fluid motions rather than gravitational settling. As the settling parameter difference |Sg 2 − Sg 1 | increases, we might expect τ g to become smaller than τ a , because the relative settling velocity tends to decrease the time required for two particles initially in the same turbulent eddy to end up in two different eddies, thereby losing their correlation. The functional form of this dependence is unknown, but because we limit this work to small values of the settling parameter, we make the simplest assumption, τ g = τ a . We will revisit this assumption when comparing the theory with experimental data.

The general case
We now have the necessary equations to obtain expressions for the RDF that include the effects of gravity and electrostatic repulsion. Substituting the equations for w p and D 12 (equations (9) and (10)) into the steady state balance equation for the RDF (equation (7)), we obtain This equation can be simplified by separating variables and using parallel notation as that of Chun et al [21] in defining c 1 and r c , The parameter c 1 is defined as and can be interpreted as the ratio of the inward inertial drift of the particles to an outward turbulent diffusive 'drift' (D 12 − D )/r . Note that settling, even in the monodisperse case, will likely change the relative sampling of strain and rotation experienced by particles, so we may expect some change in c 1 . We might expect that this would be a reduction because of increased randomization caused by gravitationally induced 'decoupling' of the particle from the fluid, and indeed the experimental results presented later seem to support that conclusion. The parameter r c can be defined as the scale at which turbulent dispersion is comparable to the accelerationinduced diffusivity, B nl r 2 c /τ η = D . Physically, it will be seen that r c manifests itself as a length scale below which the diffusivity due to particles 1 and 2 having different sizes causes the RDF to flatten. Using the expression for D derived in the last section and the assumption τ g = τ a , we can write This is the same r c as that from Chun et al [21], with the exception of the additional settling term when particles are in a gravitational field. Note that r c = 0 for a monodisperse distribution of particles, even if they are settling. We note that equation (22) derived here is similar in functional form to an expression stated by Ayala et al [12], apparently obtained empirically from their direct numerical simulation results. Their expression differs by a factor of π/8 instead of 1/3, and the presence of 404.61/Re λ instead of τ a /B nl τ η . The latter can be restated as τ g = τ a = 37.467τ η /Re λ , which is consistent with our earlier assumption of τ g = τ a , but suggests a weak Reynolds number dependence in contrast to our assumption that the time scales are constant. Significantly, their expression for τ g is independent of Sg, also consistent with our earlier assumption.
Finally, we obtain the RDF by integrating and then exponentiating equation (20), where A is an integration constant. At this point, the specifics of whether the distribution of particles is monodisperse versus bidisperse and uncharged versus charged determine the final functional form of the RDF in equation (23). We note, however, that this expression already depends explicitly on the three dimensionless parameters introduced in section 1, St, Sg and Ct. The dependence on the settling parameters can be made more explicit by rewriting equation (22) as In the following subsections, specific expressions for g(r ) under various conditions are obtained from equation (23). The first three are previous results but we show them here for completeness and to verify that the new expression simplifies to previously derived limits. Once specific expressions are derived, the theory will be compared to experiments in which the three dimensionless parameters are varied.

Monodisperse uncharged particles
For monodisperse (RDF denoted as g as opposed to g 12 ) uncharged particles, both r c and Ct are zero. Then letting A = c 0 η c 1 , equation (23) simplifies to the form given by Chun et al [21], where c 0 is an integration constant that is linked to the scale at which the RDF transitions from a power law to its inertial range form [22].

Bidisperse uncharged, settling particles
For bidisperse uncharged particles, Ct 12 is zero. Then letting A = c 0 (η 2 + r 2 c ) c 1 /2 , equation (23) has the form derived by Chun et al [21], but with the difference that r c now includes effects of gravitational settling when g = 0. Again, c 0 is an integration constant that is linked to the RDF transition scale.

Monodisperse charged particles
For monodisperse charged particles, r c is zero and Ct 12 (equation (3)) reduces to Then letting A = c 0 η c 1 , equation (23) simplifies to the expression of Lu et al [10], with the same physical meaning for c 0 and where Equation (28) is simply the monodisperse uncharged RDF (equation (25)) modulated by an exponential that goes to unity for large r and decays rapidly towards zero for small r .

Bidisperse charged, settling particles
12 with the same physical meaning for c 0 . The exponential term must go to unity as r goes to infinity, as with the monodisperse charged case. This requires that B = (π/2) Ct 12 η 3 /B nl r 3 c , giving a final expression for the bidisperse charged RDF, where the trigonometric identity arctan(r c /r ) = (π/2) − arctan(r/r c ) was used and As with the monodisperse charged RDF, equation (31) is essentially the bidisperse uncharged RDF (equation (26)) modulated by an exponential that approaches unity for large r and decays rapidly towards zero for small r . The exponential term, however, now depends on the transition scale r c in a nontrivial way. The bracketed r c /r − arctan(r c /r ) term in the exponent yields the monodisperse expression (equation (28)) in the limit r c /r 1, using the expansion arctan(z) ≈ z − z 3 /3. In the opposite limit, r c /r 1, the bracketed term becomes r c /r − π/2 ≈ r c /r .

Turbulence chamber, optics and digital holography
In order to evaluate the diffusion-drift theory for charged particles in the presence of a gravitational field, we have carried out an experimental study of particles in turbulence within the relevant parameter range (St, Sg, Ct). In this section, we describe the experimental methodology and then in section 6 we present the results. The turbulence chamber, optics and digital holographic method are described elsewhere [10], [28]- [30], but we review them briefly here. A diagram of the optical configuration and a picture of the chamber are shown in figure 1.
The turbulence chamber is an air-filled cube 50 cm on a side with a speaker-driven jet at each vertex, inspired by the chamber of Hwang and Eaton [31]. The speakers are driven by random uncorrelated voltages, and the combined effect generates approximately homogeneous and isotropic turbulence as measured by laser-Doppler velocimetry and three-dimensional (3D) particle tracking. Particle positions and velocities are obtained using digital inline holography with a two-camera system, subsequent digital reconstruction of the holograms, particle matching between the two cameras and 3D nearest-neighbour particle tracking [28].
The optical system consists of a frequency doubled Nd:YLF laser double-pulsed for velocimetry (20 ns pulse duration, 800 µs separating a double-pulse pair and 3 s between double-pulses). The beam is passed through a beam expander and spatial filter, and is then split into two beams and guided onto two CCD cameras such that the beams intersect orthogonally at the centre of the turbulence chamber. The purpose of the two-camera system is to minimize the chance of spurious particle detection and to improve the precision of particle spacing. The holograms are recorded by two Illunis XMV-11000 cameras with 4008 × 2672, 9 µm pixels and the resulting measurement volume over which the views of the two cameras overlap is approximately 36.1 × 24.0 × 36.1 mm 3 .
The energy dissipation rate ε was determined from second-order longitudinal structure functions obtained via 3D particle tracking [28,32]. The root-mean-square velocity fluctuation The turbulence chamber is a cube with a speaker at each vertex. The optical system consists of a double-pulsed laser, a beam expander with spatial filter, and a beamsplitter and mirrors to guide the expanded beam through the centre of the chamber and onto the two CCD cameras. u rms was also found from the velocities of the tracked particles. From these, the Taylormicroscale Reynolds number Re λ , the Kolmogorov length scale η = (ν 3 /ε) 1/4 , the Kolmogorov time scale τ η = √ ν/ε, the Kolmogorov velocity scale u η = η/τ η and the Kolmogorov acceleration scale a η = η/τ 2 η were calculated. These turbulence parameters are listed in table 1.

Particles
The turbulent flow is seeded with water droplets generated by a vibrating-orifice aerosol generator (TSI model 3450). The droplet generator was placed at the top of the turbulence chamber and the droplets were allowed to settle into the chamber. The droplets are typically generated with a diameter of approximately 35 µm (varies a little from experiment to experiment) with a dispersion of about ±3 µm. After being generated, many of the droplets collide and coalesce after leaving the orifice to form 'doublets' and 'triplets' with diameters of approximately 44 and 50 µm, although the triplets are too few in number to get well-resolved RDFs and therefore will not be used in this paper. In order to calculate dimensionless parameters, it is necessary to know the droplet diameters as accurately as possible. The droplet size distribution is initially characterized using a threshold (i.e. pixel counting) method in the reconstructed holograms. The relative droplet sizes can be determined quite precisely using this method [28] but variations in laser intensity due to spatial filtering and other optical adjustments, as well as variations in the laser itself, lead to a constant offset for all particle diameters that change from experiment to experiment. The offset is found using the known diameter ratios that exist between singlets, doublets and triplets and comparing those to the well-defined peaks in the size distribution obtained from the holography. As such, the diameter of the singlets for a given experiment is determined with an uncertainty of about ±3 µm. The obtained singlet diameter is also in agreement with the diameter calculated by applying the Bernoulli equation to flow at the droplet generator exit: the particle generation frequency, the diameter of the orifice and the water pressure of the water feed in the droplet generator are known, and a resultant diameter can be calculated. Once the particle diameters are obtained from the holography, the size distribution can be divided into bins corresponding to singlets and doublets, thereby allowing monodisperse and bidisperse statistics to be calculated.
Electrical charging of the particles is an inherent part of the generation process itself, which involves breakup of a water jet passing through an orifice. The magnitude of the charging (charge-to-mass ratio) is inversely related to the orifice size [8]. The distribution of electric charges per singlet drop was checked by sending the stream of particles into the uniform electric field existing between two oppositely charged, parallel conducting plates. The particles are allowed to reach their Coulomb terminal speed, which is then measured with a laser-Doppler velocimeter. The width of the distribution for each diameter mode (singlet, doublet, etc) was found to be accounted for by the width of the diameter distribution obtained from the holography, and doublets and triplets were confirmed to have double and triple the charge as singlets. Thus, the charge-to-mass ratio for the droplets is relatively uniform, with all singlets having approximately the same charge, doublets having double that charge and so on.
For the individual experiments presented in this paper, the droplet charges were measured by sending the droplet stream from the vibrating orifice into a metal collector and measuring the current between the collector and the ground with an electrometer (Keithley 642 electrometer) when at steady state (net charge of the cup is at equilibrium). A flow of air was drawn out of the bottom of the collector through a metal screen by a pump to prevent droplets from escaping from the top of the collector. The charge of the droplets is determined using the measured current, the production rate of singlet drops set on the vibrating-orifice aerosol generator, and the narrowness of the droplet charge distribution found previously with laser-Doppler velocimetry, resulting in an uncertainty of approximately 10%.

Controlling effective gravity
In order to investigate the effect of gravity on spatial clustering over a wider range of Sg without changing parameters such as St and Ct, a vertical approximately uniform electric field is externally imposed inside the turbulence chamber. The resulting downward gravitational and electrostatic forces experienced by the particles are collectively an 'effective gravity' g e that can be independently controlled over a limited range. Since the charge distribution for singlet droplets is narrow and the charge-to-mass ratio is constant for singlets, doublets, etc, the increased g e is approximately the same for the different particle sizes. To generate the electric field, the top and bottom of the turbulence chamber were lined with metal to make the chamber into a parallel plate capacitor (there is a hole in the top so that particles from the generator can enter the chamber). The plates are connected to an adjustable high-voltage power supply (up to 10 kV) to apply a controllable vertical electric field.
The imposed electric field is not completely vertical and uniform throughout the entire chamber because of the cubic aspect ratio and the fact that the top plate is approximately 50 × 50 cm 2 and the bottom plate is 35 × 35 cm 2 . We can be confident, however, that by symmetry the field is vertical and approximately uniform along the line between the centre of the top and bottom plates and the horizontal plane through the middle of the chamber. In these experiments, we set the turbulence to be relatively weak, with ε = 0.02 m 2 s −3 , so that the particles settle from the top where they enter the chamber downward through the centre of the chamber (measurement volume of the two-camera system) without exploring a large portion of the chamber. This means that they spend most of their time in the regions of the chamber where the electric field is vertical and uniform, so it is reasonable to characterize the experiments with a single g e . The nonuniformity in electric field also implies that the effective gravity cannot simply be found by measuring the potential difference between the two plates, however. Instead, the effective gravity was determined empirically by measuring the difference in terminal velocity between singlets and droplets, and then using well-established drag laws along with the known particle diameters [33]. This gives the difference in terminal velocities as g e (τ p ), and as a result, deviations from the ideal particle terminal velocity due to turbulence, mass loading, etc are all included in this effective gravity. (It should be noted, however, that the mass loading and volume fractions are sufficiently small that deviations due to collective effects should be negligible [33].)

Applying theory to experimental results
Three experiments were conducted (referred to as runs 1, 2 and 3) with close to the same singlet droplet diameter (differences are smaller than the uncertainty) but different effective gravity. All had enough singlets and doublets for well-resolved monodisperse and bidisperse RDFs but not enough triplets. The effective gravity, the approximate mean interparticle spacing (n −1/3 ∼ λ where n is the particle number density) and the particle parameters for runs 1-3 are given in table 2. Since n −1/3 is about 7-10 mm for all runs and 10η = 6.4 mm, the particle distribution is dilute on the dissipation scale.
All RDFs were calculated by the methods outlined by Saw [22]. This consists in estimating g(r ) and g 12 (r ) for discrete spherical shells of finite thickness followed by normalization utilizing the large-scale RDF (RDF from a particle field consisting of all detected particles from all laser pulse pairs). This normalization serves to remove apparent clustering caused by cut-off of spherical shells by the finite measurement volume, and apparent spatial nonuniformity due to variations in particle detection probability over the view volume (e.g. due to laser intensity variations). The normalization is valid when the apparent clustering is independent of and scale separated from the clustering of interest, both of which are valid in this experiment [22].
The calculated RDFs can be compared to the theory (equations (28) and (31)) in a phenomenological way by fitting the 'constants' that either are not determined by the theory itself, such as c 0 , or have theoretical or empirical values from computational work, which can Table 2. Effective gravity (g e ), the approximate interparticle spacing (n −1/3 ∼ λ, where n is the particle number density) and particle parameters for the individual experiments for the best estimate of the singlet particle diameter. The particle parameters are the particle diameters, charges, Stokes numbers St, settling parameters Sg, and monodisperse and bidisperse Ct ratios. A numerical subscript of 1 or 2 indicates that the parameter is for singlets or doublets, respectively. Ct with one subscript is monodisperse for that particle size and that with two subscripts is bidisperse between the particle sizes indicated by the two subscripts. Run  then be compared to the fitted values, such as c 1 and c 2,mono for monodisperse particles or c 1 , c 2,bi and r c for bidisperse particles. All fits are done for r between 400 µm, where the size of the radial shell gives only a small number of pair counts and the finite thickness of the shell begins to cause a significant bias, and 15η so as to only include the dissipation range. The fits were accomplished by weighted (based on counting uncertainty) nonlinear least squares on a semilog scale (log r ). In order to ascertain the effect of the uncertainty in the singlet diameter, the fits are performed using the best estimate singlet diameter, and also that diameter ±3 µm, with the effective gravity recalculated for each one (except for run 1, for which g e = g exactly).

Clustering of monodisperse electrically charged particles
Monodisperse RDFs for runs 1-3 with fits (equation (28)) are shown in figure 2. The fitted parameters along with the ratio α of the fitted c 2,mono to the theoretical value (equation (29)) are listed in table 3. The fit curves are independent of the assumed diameter of the singlet particles but the fitted values of c 2,mono are not since Ct is a function of the droplet diameter. Hence, the presented uncertainty in the fitted value of c 2,mono reflects both the uncertainty from the nonlinear fitting and the uncertainty in the singlet particle diameter. Qualitatively, equation (28) fits the experimental RDFs in figure 2 well. Since c 0 is not determined by the theory presented in this paper and depends on the transition to the inertial range [22], it is not of primary interest here. The fitted values of c 1 , the inertial clustering power law exponent, are in the range of or less than (for runs 2 and 3 with higher g e ) the values obtained for monodisperse uncharged particle distributions in direct numerical simulations, which is consistent with size dispersion and settling [22].
The reduction in c 1 is most pronounced for the doublets in runs 2 and 3. These are the cases for which settling is the strongest, as quantified by Sg in table 2. The inertial clustering exponent c 1 is defined in terms of the difference in the amount of strain and rotation sampled by the particles, as described by equation (21). The observed pattern in c 1 may indicate that gravitational settling causes inertial particles to more evenly sample regions of fluid experiencing strain and rotation, thereby reducing the strength of inertial clustering. There is not enough information to conclude anything about the effect of electrostatic repulsion on the Table 3. Fitted parameters for fitting equation (28) to the monodisperse charged RDFs in figure 2. α is the ratio of the fitted value of c 2,mono to the theoretical value (equation (29)). Uncertainties reflect the uncertainty from the fits and the uncertainty in the singlet particle diameter.  [10] also found that α was approximately 2 for the two experiments presented, which had St of 0.3 and 0.6 with values of Ct of 1.3 and 0.85, respectively, and an effective gravity of g. We find this consistent offset to be compelling, but at this point we can only speculate as to possible explanations. Firstly, the values of Ct tested here may be sufficiently large so as to invalidate the assumption that the Coulomb terminal velocity is simply added to each particle's velocity in finding the mean 'drift' velocity (section 3.2). That could lead to an increase in the coefficient in the exponential term c 2,mono , although it would seem remarkable that there would be a consistent bias of 2 in multiple experiments. Secondly, the droplet charge could have been underestimated for all particles, perhaps due to the presence of ambient ions in the air that reduced the charge measurement. That such a process would yield exactly a factor of 2 under all conditions also seems remarkable. Thirdly, the results may indicate that B nl is approximately half the size determined by the model utilized by Chun et al [21]. This could be a result of either limitations or inaccuracies in that model for turbulence in general, or limitations of the model when applied to particles with the parameter ranges investigated here. Further investigation is needed to understand the α ≈ 2 observation.
Finally, we revisit the assumptions made in section 3.3 of weak coupling between turbulence-induced fluctuations and gravitational or Coulomb velocities. These were expressed as Sg 1 and Ct 12 1, and the implication was that if either dimensionless parameter is much greater than unity, or in other words if the gravitational or Coulomb velocity is much greater than the Kolmogorov velocity, then strong coupling would occur. The sedimentation parameter is of the order of 1 for all experimental conditions, but may be on the borderline of acceptability for doublet particles in runs 2 and 3 (Sg = 4.7 and 5.7, respectively). Computational work shows significant coupling for Sg 1, but the transition bears further study [11,26]. For now, we take the agreement with theory, given reasonable fitted parameters over a range of effective gravity, as suggestive that the assumption is valid. The assumption involving Ct 12 is more troublesome because of the inherent r −2 dependence of the Coulomb force. The Coulombturbulence number is less than 1 for all conditions except for doublet particles in runs 2 and 3, where it is of the order of 1 (Ct = 1.07 for both). The Coulomb-turbulence number is based on droplet separation of η, which is 640 µm, and given the r −3 dependence of the velocity ratio u q (r )/u(r ) (essentially a scale-dependent Ct 12 ; see section 3.3), we would expect nearly a factor of 10 increase at η/2. Below that scale, the assumption of negligible coupling may be violated: the theory is confined, therefore, to scales within the dissipation range (r 20η), but not so small that Coulomb forces become strongly dominant. As anticipated, however, the number of particle pairs separated by less than such length scales is very small: figure 2 shows g(r ) 0.2 for all r < η/2. Again, the agreement with theory suggests that the assumption is reasonable. We leave further consideration of the transition to significant coupling to future work.

Clustering of bidisperse electrically charged particles
Bidisperse RDFs for runs 1-3 (with several fits of equation (31), to be discussed later) are shown in figure 3. Comparing g 12 of run 1 for g e = g to that of the other experiments, the effect of increasing the difference in settling parameters by increasing effective gravity is to flatten or reduce the decay of the RDF with decreasing scale. This matches the prediction from equation (24) that r c is an increasing function of the difference in settling parameters and the fact that increasing r c in the derived bidisperse charged RDF (equation (31)) flattens and reduces the decay that comes with the decreasing scale. Qualitatively, this can be interpreted as settling droplets filling the Coulomb 'exclusion' region around particles as gravitational forces increase in strength relative to electrostatic forces.
In principle, the fit of equation (31) can be carried out with four free parameters (c 0 , c 1 , c 2,bi and r c ), but we wish to evaluate the theoretical expressions for c 2,bi and r c obtained in section 4, so we plot fits using three combinations: constraining r c to be its theoretical value (equation (22) with τ g = τ a = 1.5 τ η , black curve in figure 3); constraining r c as before and c 2,bi to its theoretical value (equation (32), red thin curves in figure 3); and constraining r c as before and c 2,bi to twice its theoretical value, as suggested from the experimental monodisperse charged RDFs in section 6.1 (green dashed curves in figure 3). As in the monodisperse analysis, c 0 and c 1 were not constrained because there is little theory for c 0 , and c 1 is expected to deviate somewhat due to the influence of gravitational settling on the 'sampling' of the flow by particles (cf equation (21)). For the two fits with c 2,bi constrained, curves are shown for the best estimate of singlet droplet diameter and for the nominal singlet diameter ±3 µm to span the uncertainty range.
Qualitatively, the different fits match the experimental radial distribution functions in figure 3 and several of the constraint combinations cannot be obviously preferred given the experiment data uncertainties. Even with both r c and c 2,bi constrained, the theory captures the tendency for gravitational 'diffusion' to overcome electrostatic repulsion as g e increases. It is apparent, however, that of the two, the curve fit with c 2,bi constrained to be twice its theoretical value (green dashed curves) is closer to the experimental RDFs for all three experiments. This qualitative pattern for c 2,bi can also be seen in the last column of table 4, which lists the values of α (the ratio of the fitted c 2,bi to the theoretical value in equation (32)) coming from fits for which only r c is constrained to be its theoretical value (equation (22)). This is the same pattern as that found for α with the monodisperse charged RDFs (section 6.1), and the same possible explanations apply here as well. It should be noted, however, that constraining r c to equation (22) with τ g = τ a = 1.5τ η does not take into account our speculation that τ g may be a function of the difference in settling parameters and therefore not have a simple, constant value.
In order to evaluate the representation of gravitational settling in the theory, we are particularly interested in the behaviour of r c as effective gravity and therefore the relative  settling speed between singlet and doublet particles is varied. We consider the situation with c 2,bi constrained to be the theoretical value and twice the theoretical value. Those two situations result in the r c values given in table 4, and are seen to increase in a similar manner to the  (22)) along with fitted values (g e /a η ) 2 after dividing both sides by the square of the difference in Stokes numbers (St 2 − St 1 ) 2 . Hence, the data should follow a linear relationship between r 2 c /η 2 (St 2 − St 1 ) 2 and (g e /a η ) 2 insofar as this value of τ g is correct. Using the fits to the bidisperse RDFs in which c 2,bi is constrained to be its theoretical value and two times its theoretical value, this quantity is plotted in figure 4. The theoretical linear dependence when τ g = τ a = 1.5τ η is shown for comparison. Indeed, for a given constrained value of c 2,bi in figure 4, r 2 c /η 2 (St 2 − St 1 ) 2 is an increasing function of (g e /a η ) 2 , as expected. Forcing c 2,bi to be twice its theoretical value leads to significantly improved agreement with the theory line, indicating that τ g = τ a = 1.5τ η is valid for the range of Sg tested if c 2,bi is actually double its theoretical value. This provides another reason to further investigate c 2,bi , just as in the case of c 2,mono .

Conclusions
Building on the drift-diffusion model of Chun et al [21], the contributions of electrostatic interactions between particles and gravitational settling in homogeneous isotropic turbulence to particle clustering have been developed. The theory is argued to hold under the following dimensionless conditions. The particle Stokes number, representing the particle inertial time scale relative to turbulence dissipation time scale, is small, St 1; the Coulomb-turbulence parameter, representing the ratio of Coulomb 'terminal' speed to turbulence dissipation velocity scale, is of the order of 1 or smaller, Ct 1; and the settling parameter, representing the ratio of the gravitational terminal speed to turbulence dissipation velocity scale, is of the order of 1 or smaller, Sg 1. These conditions allow for the simplification that settling and electrostatic interactions act independently, leading to a superposition of their respective terminal velocities to the velocities particles would experience without settling and electrostatic interactions. With these contributions, the drift-diffusion model (equation (7)) was used to obtain the functional form of the RDF to quantify particle clustering for a bidisperse distribution of electrically charged particles (both settling and nonsettling) and re-obtain the form of the RDF derived by Lu et al [10] for a monodisperse distribution of charged particles. The derived RDFs turn out to be the RDFs for uncharged particles multiplied by an exponential that depends on electrostatic and turbulence quantities (bidisperse and monodisperse forms are slightly different). The functional form decays with decreasing scale and goes to unity in the limit r → ∞.
In order to test these extensions to the drift-diffusion model and the resulting functional forms of the RDFs, experiments were performed on electrically charged water droplets (narrow distribution of charge per droplet) in a speaker-driven turbulence chamber (approximately homogeneous and isotropic turbulence), with 3D particle positions obtained using dual-camera digital inline holography. Two different sizes of water droplets were used, with the larger size being the result of coalescence of two droplets of the smallest size. The effective gravity experienced by the water droplets was controlled by applying an approximately vertical and uniform electric field between conducting plates lining the top and bottom of the turbulence chamber.
The derived monodisperse charged RDF fits well to the RDFs obtained in the experiments, which was seen previously by Lu et al [10]. As in that previous work, the coefficient inside the exponential term of the RDF (c 2,mono ) appears to be twice as large as the theory presented in this paper predicts. This apparent discrepancy is as yet unexplained. Possible reasons include underestimation of the electrical charge on each droplet or a flaw in one or more of the assumptions in the theory, although it would be surprising for such problems to yield a consistent factor of 2 bias; a further possible explanation, but providing a consistent bias, would be a difference in the Batchelor-type dispersion constant from that predicted by the model used by Chun et al [21]. Either way, this warrants further investigation. The inertial clustering exponent c 1 was observed to decrease slightly, but monotonically with increasing g e , and this also is an area in which further theory is needed.
The derived bidisperse charged RDF also fits the experimental data, even when the coefficient in the exponential term (c 2,bi ) and the clustering saturation scale (r c ) that leads to flattening of the RDF are both constrained to the theoretical values. The theory better matches the observations when based on twice the value for c 2,bi derived in this paper. As in the monodisperse case, the peculiarity of this factor of two in the Coulomb coefficient remains an open question. Finally, the dependence of r c on g e is reasonable over the range studied, which spans 1-2.3g, but further studies are needed for investigating the relative roles and magnitudes of the Lagrangian acceleration correlation time τ a and the gravitational diffusivity correlation time τ g .
Future experimental work aimed at addressing these uncertainties will require more precisely resolved RDFs, with the particle size, particle charge and effective gravity more accurately determined. Even in the current state, however, the theory is suitably confirmed by the experimental measurements to justify further investigation of the role of charge and settling in various settings, from clouds to dusty plasmas.