the Creative Commons Attribution 4.0 License.

the Creative Commons Attribution 4.0 License.

# Intricate relations among particle collision, relative motion and clustering in turbulent clouds: computational observation and theory

### Ewe-Wei Saw

### Xiaohui Meng

Considering turbulent clouds containing small inertial particles, we investigate the effect of particle collision, in particular collision–coagulation, on particle clustering and particle relative motion. We perform direct numerical simulation (DNS) of coagulating particles in isotropic turbulent flow in the regime of small Stokes number (*St*=0.001–0.54) and find that, due to collision–coagulation, the radial distribution functions (RDFs) fall off dramatically at scales *r*∼*d* (where *d* is the particle diameter) to small but finite values, while the mean radial component of the particle relative velocity (MRV) increases sharply in magnitude. Based on a previously proposed Fokker–Planck (drift-diffusion) framework, we derive a theoretical account of the relationship among particle collision–coagulation rate, RDF and MRV. The theory includes contributions from turbulent fluctuations absent in earlier mean-field theories. We show numerically that the theory accurately accounts for the DNS results (i.e., given an accurate RDF, the theory could produce an accurate MRV). Separately, we also propose a phenomenological model that could directly predict MRV and find that it is accurate when calibrated using fourth moments of the fluid velocities. We use the model to derive a general solution of RDF. We uncover a paradox: the past empirical success of the differential version of the theory is theoretically unjustified. We see a further shape-preserving reduction of the RDF (and MRV) when the gravitational settling parameter (*S*_{g}) is of order *O*(1). Our results demonstrate strong coupling between RDF and MRV and imply that earlier isolated studies on either RDF or MRV have limited relevance for predicting particle collision rate.

- Article
(1939 KB) -
Supplement
(791 KB) - BibTeX
- EndNote

The motion and interactions of small particles in turbulence have fundamental implications for atmospheric clouds; specifically, they are relevant for the timescale of rain formation, particularly in warm clouds (Falkovich et al., 2002; Wilkinson et al., 2006; Grabowski and Wang, 2013) (a similar problem also applies to planet formation in astrophysics; Johansen et al., 2007). They are also important for engineers who are designing future, greener combustion engines, as this is a scenario they wish to understand and control in order to increase fuel efficiency (Karnik and Shrimpton, 2012). Cloud particles or droplets, due to their inertia, are known to be ejected from turbulent vortices and thus form clusters – regions of enhanced particle density (Wood et al., 2005; Bec et al., 2007; Saw et al., 2008; Karpińska et al., 2019); this together with droplet collision is of direct relevance for the mentioned applications. Due to the technical difficulty of obtaining extensive and systematic experimental or field data on particle/droplet collision in turbulent cloud, many of the recent studies rely on direct numerical simulation (DNS), examples of which can be found in, e.g., Onishi and Seifert (2016) and Wang et al. (2008), and reference therein. Up until now, we have not had definitive answers to basic questions such as how to calculate the particle collision rate from basic turbulence-particle parameters and what the exact relation between collision and particle clustering and/or motions is, for, as we shall see, our work reveals that collision–coagulation causes profound changes in particle relative velocity statistics and particle clustering, questioning earlier understanding of the problem. The difficulty of this problem is in part related to the fact that turbulence is, even by itself, theoretically virtually intractable due to its nonlinear and complex nature.

The quest for a theory of particle collision in turbulence started in 1956 when Saffman and Turner (1956) derived a mean-field formula for the collision rate of finite-size, inertia-less particles. In another landmark work (Sundaram and Collins, 1997), a general relation among the collision rate (*R*_{c}), particle clustering and mean particle relative radial velocity was presented: ${R}_{c}/\left({n}_{\mathrm{1}}{n}_{\mathrm{2}}V\right)=\mathrm{4}\mathit{\pi}{d}^{\mathrm{2}}g\left(d\right)\langle {w}_{r}\left(d\right){\rangle}_{\ast}$, where *g*(*r*) is the particle radial distribution function (RDF), *w*_{r} is the radial component of relative velocity between two particles, $\langle \cdot {\rangle}_{\ast}$ denotes averaging over particle pairs, 〈*w*_{r}(*d*)〉_{∗} is the mean radial component of relative particle velocity (MRV), *n*_{i} denotes the global averages of particle number density, *V* is the spatial volume of the domain and *d* is the particle diameter. The remarkable simplicity of this finding inspired a “separation paradigm”, which is the idea that one could study the RDF or MRV separately (which is technically easier), and the independent results from the two may be combined to accurately predicts *R*_{c} (an idea that we subsequently challenge). Another work of special interest here is the drift-diffusion model by Chun et al. (2005) (hereafter CK theory) (note that there are other similar theories: Balkovsky et al., 2001; Zaichik and Alipchenkov, 2003). The CK theory, derived for non-colliding particles in the limit of a vanishing particle Stokes number, *St* (a quantity that reflects the importance of the particle's inertia in dictating its motion in turbulence), correctly predicted the power-law form of the RDF (Reade and Collins, 2000; Saw et al., 2008) and has seen remarkable successes over the years, including the accurate account of the modified RDF of particles interacting electrically (Lu et al., 2010) and hydrodynamically (Yavuz et al., 2018).

Here, we first present results on RDF and MRV for particles undergoing collision–coagulation^{1}. The data are obtained via direct numerical simulation (DNS), which is the gold-standard computational method in terms of accuracy and completeness for solving the most challenging fluid dynamics problem i.e., turbulent flows. It is worth noting that the focus of our work is on the fundamental relationship between collision, RDF and MRV and on highlighting differences from the case with non-colliding particles (Chun et al., 2005). To that end, we have designed the DNS to have an idealized setup similar to what was done in Chun et al. (2005), which would allow us to identify the effects of particle collision–coagulation without doubt. As a result, this limits the direct applicability to real systems (these limitations are detailed in Sect. 4.5).

Analysis of the DNS results is followed by a theoretical account of the relations between collision rate, RDF and MRV, which includes mean-field contributions (Saffman and Turner, 1956; Sundaram and Collins, 1997) and contributions from turbulent fluctuations (absent from earlier theories, Saffman and Turner, 1956; Sundaram and Collins, 1997). The theory is derived from the Fokker–Planck (drift-diffusion) framework first introduced in the CK theory (Chun et al., 2005).
We shall see that the main effect of collision–coagulation is the enhanced asymmetry in the particle relative velocity distribution^{2} and that this leads to nontrivial outcomes.

To observe how particle collision–coagulation affects RDF and MRV, we performed direct numerical simulation (DNS) of steady-state isotropic turbulence embedded with particles of finite but sub-Kolmogorov size. We solve the incompressible Navier–Stokes equations (Eq. 1) using the standard pseudo-spectral method (Rogallo, 1981; Pope, 2000; Mortensen and Langtangen, 2016) inside a triply periodic cubic box.

where *ρ*, *p*, *ν* and **f** are the fluid mass density, pressure, kinematic viscosity and imposed forcing respectively. The velocity field is discretized on a 256^{3} grid.
Aliasing resulting from Fourier transform of truncated series is removed via a $\mathrm{2}/\mathrm{3}$ dealiasing rule (Rogallo, 1981).
A statistically stationary and isotropic turbulent flow is achieved by continuously applying random forcing to the lowest wave numbers until the flow's energy spectrum is in steady state (Eswaran and Pope, 1988).
Second-order Runge–Kutta time stepping was employed. Further details of such a standard turbulence simulator can be found in, e.g., Pope (2000), Rogallo (1981) and Mortensen and Langtangen (2016).
The accuracy of DNS for turbulent flows has been experimentally validated for decades (see, e.g., the compilation of results in Pope, 2000).

Particles in the simulations are advected via a viscous Stokes drag force (Maxey and Riley, 1983):

where **u** and **v** are the local fluid and particle velocity respectively, and *τ*_{p} is the particle inertia response time, defined as ${\mathit{\tau}}_{p}=\frac{\mathrm{1}}{\mathrm{18}}({\mathit{\rho}}_{p}/\mathit{\rho}-\mathrm{1})({d}^{\mathrm{2}}/\mathit{\nu})$, where *ρ*_{p} is the particle mass density, and *d* is the particle diameter.
As mentioned, this work focuses on the fundamental relationship between collision–coagulation, RDF and MRV, as well as on addressing the validity of the theory (to be described). It is thus, beneficial to keep the DNS setting idealized (and in the regime relevant for the theory) for the sake of clarity when interpreting results. To that end, the DNS does not include inter-particle hydrodynamic interactions (HDIs) and gravitational settling, nor does it consider the effects of temperature and humidity variation and phase transitions. Such practice is not uncommon in studies designed to isolate and address fundamental issues related to particles dynamics in turbulence; examples that are closely related to the current setup and/or problem include Sundaram and Collins (1997), Chun et al. (2005), Bec et al. (2007), Salazar et al. (2008), Wang et al. (2008), Woittiez et al. (2009) and Voßkuhle et al. (2013). However, such an approach certainly limits the direct applicability of our results to some realistic problems in the atmosphere; these limitations will be detailed in Sect. 4.5, where a discussion of the effects of gravity and HDI is also given.

In this context, the particle Stokes number, defined as ${\mathit{\tau}}_{p}/{\mathit{\tau}}_{\mathit{\eta}}$, where *τ*_{η} is the Kolmogorov timescale, can be expressed as $\mathit{\text{St}}=\frac{\mathrm{1}}{\mathrm{18}}({\mathit{\rho}}_{p}/\mathit{\rho}-\mathrm{1})(d/\mathit{\eta}{)}^{\mathrm{2}}$, where *η* is the Kolmogorov length scale. Time stepping of the particle motion is done using a second-order modified Runge–Kutta method with an “exponential integrator” that is accurate, even for *τ*_{p} much smaller than the fluid's time step (Ireland et al., 2013). The particles introduced into the simulation are spherical and are of the same size, the initial number of particles is 10^{7} and they are randomly distributed in space. Particles collide when their volumes overlap and a new particle is formed conserving volume and momentum (Bec et al., 2016).
We continuously, randomly, inject new particles so that the system is in a steady state after some time. Statistical analysis is done at steady state on monodisperse particles (i.e., particles with the same *St*). Experimental validation of the accuracy of such a particle simulation scheme in DNS can be found in Salazar et al. (2008), Saw et al. (2012, 2014) and Dou et al. (2018).

Values of key parameters of the DNS are given in Table 1. Values of other parameters and further details can be found in the Supplement.

As described in Chun et al. (2005), in the limit of *St*≪1, particle motions are closely tied to the fluid velocity and, to leading order, completely specified by the particle position and fluid velocity gradients. We consider the Fokker–Planck equation which is closed and deterministic (see, e.g., Appendix J in Pope, 2000):

where $P\equiv P({r}_{i},t|{\mathrm{\Gamma}}_{ij}\left(t\right))$ is the (per volume) probability density function (PDF) for a secondary particle
to be at vector position *r*_{i} relative to a primary particle at time *t*, conditioned on a fixed and known history of the velocity gradient tensor along the primary particle's trajectory Γ_{ij}(*t*), and *W*_{i} is the mean velocity of secondary particles relative to the primary, under the same condition. Note that *W*_{i} is a conditional average, while *w*_{i} denotes a realization of relative velocity between two particles.

From this, one could derive an equation for 〈*P*〉:

where 〈.〉 implies ensemble averaging over primary particle histories (note that 〈*W*_{r}〉≡ unconditional mean of *w*_{r}, averaged over all particle pairs, i.e., the MRV). This equation, however, is not closed due to the correlation between the fluctuating terms *W*_{i} and ${P}^{\prime}\equiv P-\langle P\rangle $.
The correlation 〈*W*_{i}*P*^{′}〉 can be written in terms of a drift flux and diffusive flux (detailed derivation is well described in Chun et al., 2005), such that we have

where the drift flux is

and the diffusive flux is

where **r**^{′} satisfies a characteristic equation:
$\partial {r}_{i}^{\prime}/\partial {t}^{\prime}={W}_{i}({\mathbf{r}}^{\prime},{t}^{\prime})$, with boundary condition ${r}_{i}^{\prime}={r}_{i}$ at ${t}^{\prime}=t$.

Finally we note that since particles are allowed to collide–coagulate in our theory, we use the conventional definition of MRV: $\langle {W}_{r}\rangle \equiv \langle {w}_{r}{\rangle}_{\ast}$. In some works that consider non-colliding (ghost) particles, the conditional mean $\langle {w}_{r}|{w}_{r}\le \mathrm{0}{\rangle}_{\ast}$ must be used for the purpose of calculating the mean collision rate, since $\langle {w}_{r}{\rangle}_{\ast}=\mathrm{0}$ in isotropic turbulence (Chun et al., 2005).

We compute the RDF via
$g\left(r\right)={N}_{\text{pp}}\left(r\right)/\left[\frac{\mathrm{1}}{\mathrm{2}}N\right(N-\mathrm{1})\mathit{\delta}{V}_{r}/V]$,
where *N*_{pp}(*r*) is the number of particle pairs found to be separated by distance *r*, and *δ**V*_{r} is the volume of a spherical shell of radius *r* and infinitesimal thickness *δ**r*.

Figure 1 shows the RDFs obtained for monodisperse particles of various Stokes numbers and sizes. Two cases (*St*=0.22 and 0.54) are shown in panel a, and two more (*St*=0.054 and 0.001) are shown in panel b. In this work, we focus on the smaller values of *St* since
the theory which we shall consider is also only applicable in the *St*≪1 regime. However, we have included the *St*=0.54 case to demonstrate that the observations to be described also extends to finite *St*. In all cases, except one, the particles are of the same size $d={d}_{\ast}$, where *d*_{∗} represents the specific value of ${d}_{\ast}=\mathrm{9.49}\times {\mathrm{10}}^{-\mathrm{4}}$ dm, chosen so that the particle sizes are about *O*(0.1) times the Kolmogorov scale (*η*), thus allowing us to still observe a regime (3*d**≲**r**≲*30*η*) of power-law RDFs. To shows the effect of changing particle size, panel a also includes a case of $\mathit{\text{St}}=\mathrm{0.54},d=\mathrm{2}{d}_{\ast}$ for comparison. Looking at panel a, apart from the apparent power-law behavior of the RDFs at intermediate values of *r*, the most striking feature of these RDFs for colliding–coagulating particles is that they fall off dramatically in the *r*∼*d* regime. This is very different from what was seen in earlier studies of non-colliding particles, where *g*(*r*) are simple power laws (Chun et al., 2005; Saw et al., 2008).
We also see that as *r* approaches *d*, the steepness of the curve (see, e.g., the blue circles) increases as *g*(*r*) drops off; this and the fact that the abscissa is logarithmic implies that $\partial g/\partial r$ is increasing exponentially in the process. As a consequence, it is difficult to discern from these plots if the limit of *g*(*r*) at particle contact (*r*→*d*) is still nonzero. This is an important question as ${lim}_{r\to d}\left[g\right(r\left)\right]=\mathrm{0}$ implies that the mean-field formula of Sundaram and Collins (1997) has zero contribution towards *R*_{c}; i.e., the collision rate is solely due to turbulent fluctuations. It is only by re-plotting *g*(*r*) versus *r*−*d* (see insets in Fig. 1) and using a remarkable resolution that is 10^{3} finer than *d* that we see a convincing trend supporting a finite *g*(*r*→*d*).
Also clear in panel a is the observation that with changing particle size (*d*), the location of the sharp fall-off merely shifts to where the new value of *d* is.

The strong effect of particle collision on the RDF (also on MRV as we shall see later) challenges the validity of the separation paradigm. We note that similar fall-off of RDF was previously observed (Sundaram and Collins, 1997), but a complete analysis and theoretical understanding were lacking. Also, a study on multiple collisions (Voßkuhle et al., 2013) had hinted at the potential problem with the separation paradigm.

Another observation is that in the power-law regime (3*d**≲**r**≲*30*η*), the RDFs appear (as expected) as straight lines with slopes (i.e., power-law exponents) that increase with *St* and are numerically consistent with those found for non-colliding particles (see, e.g., Saw et al., 2012).

## 4.1 Theoretical account via drift-diffusion theory

To theoretically account for the new findings, we make some derivations that are partially similar to the ones in Chun et al. (2005) but under a new constraint due to coagulations: at contact (*r*=*d*), the radial component of the particle relative velocities can not be positive^{3},
while with increasing *r* the constraint is gradually relaxed. The
first consequence of this is that the distribution of the radial
component of the relative particle velocity (*W*_{r}) is highly
asymmetric at *r*≈*d*; i.e., the PDFs of positive *W*_{r} values are very
small (this constitutes the “enhanced asymmetry” mentioned
earlier). Thus for *r*≈*d*, 〈*W*_{r}〉 must be negative.
In Sect. 3, we showed that in the *St*≪1 limit, one could derive a master equation (Eq. 4), reproduced here for clarity:

where ${q}_{i}^{d}$ is the drift flux (of probability due to turbulent fluctuation) defined in Eq. (5), and ${q}_{i}^{D}$ is the diffusive flux defined in Eq. (6).

We then expand *W*_{i}, $\partial {W}_{l}/\partial {r}_{l}$ and
(consequentially) the fluxes as perturbation series with *St*
as the small parameter (details in the Supplement or Chun et al., 2005). The coagulation constraint affects the values of the coefficients of these series. For the drift flux, the leading-order terms (in powers of *St*) are

with ${a}_{ik}^{\left(\mathrm{1}\right)}={\mathit{\tau}}_{\mathit{\eta}}\langle {\mathrm{\Gamma}}_{ik}\left(t\right){\mathrm{\Gamma}}_{lm}\left({t}^{\prime}\right){\mathrm{\Gamma}}_{ml}\left({t}^{\prime}\right)\rangle $ and
${a}_{ki}^{\left(\mathrm{2}\right)}={\mathit{\tau}}_{\mathit{\eta}}^{\mathrm{2}}\langle {\mathrm{\Gamma}}_{ij}\left(t\right){\mathrm{\Gamma}}_{jk}\left(t\right){\mathrm{\Gamma}}_{lm}\left({t}^{\prime}\right){\mathrm{\Gamma}}_{ml}\left({t}^{\prime}\right)\rangle $, where Γ_{ij} is the *i**j*th component of the fluid's velocity gradient tensor at the particle position (the *a*_{ik} values are thus related to two time correlations of moments of velocity gradients; Chun et al. (2005) showed that ${a}_{ik}^{\left(\mathrm{2}\right)}\propto \stackrel{\mathrm{\u203e}}{{S}^{\mathrm{2}}}-\stackrel{\mathrm{\u203e}}{{R}^{\mathrm{2}}}$, where $(\stackrel{\mathrm{\u203e}}{{S}^{\mathrm{2}}},\stackrel{\mathrm{\u203e}}{{R}^{\mathrm{2}}})$ denotes the average fluid (strain rate tensor, rotation rate tensor) squared at particle positions). As explained earlier, coagulation constraint causes the PDF of relative particle velocities to become highly asymmetric for *r*∼*d*; thus ${a}_{ik}^{\left(\mathrm{1}\right)}$ is nonzero at this scale. This is very different from the case of non-colliding particles (Chun et al., 2005), where ${a}_{ik}^{\left(\mathrm{1}\right)}$ is always zero due to statistical isotropy. Under the constraint, DNS gives ${\int}_{-\mathrm{\infty}}^{t}{a}_{ik}^{\left(\mathrm{1}\right)}\mathrm{d}{t}^{\prime}\approx -\mathrm{0.18}$ s^{−1}
and ${\int}_{-\mathrm{\infty}}^{t}{a}_{ki}^{\left(\mathrm{2}\right)}\mathrm{d}{t}^{\prime}\approx \mathrm{2.45}$ s^{−1}
(more in the Supplement). Thus for *r*∼*d*, the drift flux is negative for large *St* but becomes positive^{4}
when *St* is below the value of ≈0.07; and in the limit of *St*→0, it is dominated by the first term in Eq. (7).

The term ${q}_{i}^{D}$ is a “nonlocal” diffusion caused by fluctuations and can be estimated using a model that assumes the particle relative motions are due to a series of random uniaxial straining flows (Chun et al., 2005). Chun et al. (2005) showed that, generally, ${q}_{i}^{D}$ has an integral form (due to non-locality), and only in the special case where *g*(*r*) is a power law may it be cast into a differential form (similar to a local diffusion). In view of the nontrivial *g*(*r*) observed here, we must proceed with the integral form:

where ${R}_{\mathrm{0}}\equiv {r}_{\mathrm{0}}/r$, with *r*_{0} as the initial separation
distance of a particle pair before a straining event. *F* is the
probability density function for the duration of each event, *f*_{I} is
determined by relative prevalence of extensional versus compressional
strain events (more details in the Supplement or Chun et al., 2005) and Ω is the solid angle for the axis of the straining flow. Note that due to coagulation, the *R*_{0} integration starts from $d/r$. We differ crucially from the CK theory via the introduction of the (positive) factor *c*_{st}, which can be shown to be equal to $\left|{c}_{\mathrm{1}}\right|$, where *c*_{1} is the power-law exponent of the RDF the particles would have assuming they are non-colliding (details in the Supplement).

By definition, *g*(*r*)≡*α*〈*P*〉. Periodic boundaries in our DNS imply that *α*=*V*
(more in the Supplement). Using this and the fact that the problem only has radial (*r*) dependence, we rewrite Eq. (4) as

where the content inside [.] gives the total flux. For a system in
steady state, the first term in Eq. (8) is zero, and upon integrating with limits [*d*,*r*], we have

where we have identified the total flux at contact (*r*=*d*) as the negative of the (always positive) normalized collision rate ${R}_{c}^{*}\equiv {R}_{c}/\left(\mathrm{4}\mathit{\pi}\left[N\right(N-\mathrm{1})/\mathrm{2}]/V\right)$,
and comparing with Eq. (7), we see that

with the specific values of the *t*^{′} integrals already given above. For clarity, we reiterate that on the left side of Eq. (9), we have the diffusive flux (${q}_{r}^{D}$), mean-field flux (*r*^{2}*g*(*r*)〈*W*_{r}〉) and drift flux (${q}_{r}^{d}$), while on the right, the total flux is given in terms of the normalized collision rate (${R}_{c}^{*}$). We note that this equation embodies the full relationship among RDF, MRV and collision rate.

## 4.2 Ansatz and accuracy of the theory

Simple analytical solutions
to Eq. (9) may be elusive due to its integral nature (a consequence of the non-local diffusive flux). However, one could gain insights into it and test its accuracy via numerical solutions. To that end, we begin with a simple ansatz for *g*(*r*), then we curve-fit the ansatz to the DNS-produced RDF (*g*_{DNS}(*r*)). This enables us to, firstly, verify that the ansatz could accurately represent *g*_{DNS}(*r*) and, secondly, to obtain a “calibrated” ansatz that is a numerically accurate representation of *g*_{DNS}(*r*). We then show that Eq. (9), supplied with the calibrated ansatz, could numerically predict 〈*W*_{r}〉(*r*) (i.e., MRV) that agrees well with the DNS-produced MRV. In short, we will show that given a “correct” *g*(*r*), Eq. (9) produces the “correct” 〈*W*_{r}〉(*r*).

The ansatz has the form *g*_{a}(*r*)=*g*_{0}(*r*)*g*_{s}(*r*), with ${g}_{s}\left(r\right)={c}_{\mathrm{0}}{r}^{-{c}_{\mathrm{1}}}$, i.e., the RDF form for non-colliding particles (Chun et al., 2005) under the same conditions. As a first-order analysis, we let *g*_{0}, which embodies the effects of collision, take the simplest form that could still capture the main features of the RDFs seen in Fig. 1. Specifically, we let ${g}_{\mathrm{0}}\left(r\right)={c}_{\mathrm{00}}(r-{d}_{\mathrm{1}}{)}^{{c}_{\mathrm{10}}}$, where *c*_{00}(*r*) and *c*_{10}(*r*) are each piecewise constant quantities that switch from their small *r* to large *r* values at a crossover-scale *r*_{c} (of the order of *d*); i.e., *g*_{0} is a two-piece power law of *r*−*d*_{1}. (Note that our earlier finding of $g(r\to d)>\mathrm{0}$ implies that *d*_{1}<*d*.)

From a given DNS-produced RDF (*g*_{DNS}(*r*)), we first obtain a calibrated *g*_{s} by fitting ${c}_{\mathrm{0}}{r}^{-{c}_{\mathrm{1}}}$ to *g*_{DNS}(*r*) in the power-law regime *d*≪*r**≲*10*η* (see the red line in Fig. 1b). Next, we compute the DNS estimate of *g*_{0} via ${g}_{\mathrm{0}}^{\text{DNS}}={g}_{\text{DNS}}\left(r\right)/{g}_{s}\left(r\right)$, which is essentially a compensated RDF (see the ◦ plot in Fig. 1b). To get a calibrated *g*_{0}, we then fit the general form of *g*_{0} given above to ${g}_{\mathrm{0}}^{\text{DNS}}$ (see the cyan line in Fig. 1b; note that each time, two pieces of power laws are fitted to one ${g}_{\mathrm{0}}^{\text{DNS}}$, and *r*_{c} results naturally from the intersection of the two). Figure 1b shows the calibrated ansatz for the case of *St*=0.054 and verifies its accuracy (the red line is *g*_{s}(*r*), the cyan line is *g*_{0}(*r*) and the dashed black line (*g*_{0}⋅*g*_{s}) accurately reproduces *g*_{DNS}(*r*)). The inset in Fig. 1b shows that ${g}_{\mathrm{0}}^{\text{DNS}}\left(r\right)$ (i.e., ${g}_{\text{DNS}}\left(r\right)/{g}_{s}\left(r\right)$) is roughly *St*-independent for *St*≪1.

Next, we numerically evaluate the integral in the first term of Eq. (9).
The *St*≪1 assumption allows us to approximate *g*(*r*,*St*) inside the integral by its zero-*St* cousin $g(r,\mathit{\text{St}}\to \mathrm{0})$ (Chun et al., 2005).
In practice, we replace *g*(*r*,*St*) with the ansatz fitted to the DNS result of $g(r,\mathit{\text{St}}=\mathrm{0.001})$. Next, we use the DNS data to estimate *A*_{τ} and compute ${R}_{c}^{*}$ and *c*_{st} (for this case, DNS gives ${R}_{c}^{*}=\mathrm{9.69}\times {\mathrm{10}}^{-\mathrm{10}}\phantom{\rule{0.125em}{0ex}}{\mathrm{dm}}^{\mathrm{3}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$; ${c}_{\mathit{\text{st}}}=\left|{c}_{\mathrm{1}}\right|$ as mentioned earlier). Finally we use Eq. (9) to predict 〈*W*_{r}〉(*r*).

Comparison of the predicted 〈*W*_{r}〉(*r*) with the values obtained directly from the DNS is shown in Fig. 2. The prediction shown was made for the case of *St*=0.054, to be compared with its DNS counterpart (the ◦ symbols). (We also show the DNS result for *St*=0.001 and 0.11 to highlight an observation that 〈*W*_{r}〉(*r*) is almost *St*-independent in this small *St* regime.)
We showed earlier that for *r*∼*d*, *A*_{τ} is given by Eq. (10). However, as stated earlier, as *r* increases, the (statistical) asymmetry induced by collision–coagulation gradually becomes subdominant to the isotropy of turbulent fluctuation. Statistical isotropy implies ${a}_{ik}^{\left(\mathrm{1}\right)}=\mathrm{0}$ (Chun et al., 2005), a fact our DNS data confirm. Thus, for *r*≫*d*, *A*_{τ} equals the order *St*^{2} term in Eq. (10), exactly the same as the results of Chun et al. (2005) for non-colliding particles.
For this reason, we show two versions of the prediction: $\langle {W}_{r}{\rangle}_{r\sim d}^{\text{theory}}$ and $\langle {W}_{r}{\rangle}_{r\gg d}^{\text{theory}}$, which are respectively obtained by setting *A*_{τ} to its small *r* and large *r* limits ($-\mathrm{2.6}\times {\mathrm{10}}^{-\mathrm{3}}$ s^{−1}, $\mathrm{7.1}\times {\mathrm{10}}^{-\mathrm{3}}$ s^{−1}) respectively. The agreement between DNS and the predictions is noteworthy, especially for small *r*. At *r*≈2*d*, the DNS result shows a weak tendency to first follow the upward trend of $\langle {W}_{r}{\rangle}_{r\sim d}^{\text{theory}}$ and then drops off significantly at *r**≳*2.5*d*. The latter is consistent with the fact that $\langle {W}_{r}{\rangle}_{r\gg d}^{\text{theory}}$ is below $\langle {W}_{r}{\rangle}_{r\sim d}^{\text{theory}}$, but the drop is sharper than predicted.

## 4.3 Phenomenological model of MRV

Alternatively, Eq. (9) may be solved for the correct *g*(*r*), if 〈*W*_{r}〉 is given. As we are assuming *St*≪1, particle velocity statistics may be approximated by their fluid counterparts (Chun et al., 2005); i.e., we may replace 〈*W*_{r}〉 with 〈*W*_{r}〉_{St=0}, the latter being the MRV of fluid particles. Hence, if 〈*W*_{r}〉_{St=0} is known, it may be used, together with Eq. (9),
to predict RDF of any finite but small *St*. Figure 2a shows that 〈*W*_{r}〉_{St>0} from the DNS do not change significantly for $\mathit{\text{St}}\in [\mathrm{0.001},\mathrm{0.1}]$, supporting this approach^{5}.

Here we provide a simple, first-order model for 〈*W*_{r}〉_{St=0}. We limit ourselves to the regime of small particles (*d*≪*η*) and anticipate that 〈*W*_{r}〉 is non-trivial (nonzero) only for *r*∼*d*, a fact observable in Fig. 2a. We also assume that the relative trajectories of particles are rectilinear at such small scales. The coagulation constraint then implies that: in the rest frame of a particle (call it P1), a second particle nearby must move in such a way that the angle (*θ*) between its relative velocity and relative position (seen by P1)
must satisfy ${\mathrm{sin}}^{-\mathrm{1}}(d/r)\le \mathit{\theta}\le \mathit{\pi}$, under the convention of ${\mathrm{sin}}^{-\mathrm{1}}\left(x\right)\in [-\frac{\mathit{\pi}}{\mathrm{2}},\frac{\mathit{\pi}}{\mathrm{2}}]$
(more in the Supplement). We can thus write (by treating negative and positive *w*_{r} separately, applying the K41 theory (Kolmogorov, 1941) and the bounds on *θ*, details in the Supplement), for *St*≪1, that

where $\langle .{\rangle}_{\ast}$ denotes averaging over particle pairs, *p*_{+} (*p*_{−}) is the probability of a realization of *w*_{r} being positive (negative), and ${P}_{\mathit{\theta}}^{+}$ is a conditional PDF such that ${P}_{\mathit{\theta}}^{+}\equiv P\left(\mathit{\theta}|{w}_{r}\ge \mathrm{0}\right)\equiv P\left(\mathit{\theta}|\mathit{\theta}\in [\mathrm{0},\frac{\mathit{\pi}}{\mathrm{2}}]\right)$, *θ*_{m} is the lower bound of *θ* described above.
For a first-order account, we neglect skewness in the distribution of particle relative velocities and set ${p}_{\pm}=\mathrm{0.5}$. Following Kolmogorov (1941), we have set $\langle {w}_{r}|{w}_{r}<\mathrm{0}{\rangle}_{\ast}={\mathit{\xi}}_{\text{\u2013}}r$, where ${\mathit{\xi}}_{\pm}={C}_{s}\sqrt{\mathit{\epsilon}/\left(\mathrm{15}\mathit{\nu}\right)}$, (*C*_{s} is a Kolmogorov constant, we found *C*_{s}=0.76 by matching *ξ*_{–}*r* to the first-order fluid velocity structure-function from the DNS).

A simple phenomenological model for *P*(*θ*) may be constructed using the (statistical) central-limit theorem by assuming that the angle of approach *θ* at any time is the sum of many random-incremental rotations in the past, thus we write

where *N*exp […] is the circular normal distribution, i.e., analog of Gaussian distribution for angular data; sin (*θ*) results from integration over azimuthal angles (*ϕ*). We set ${\mathit{\mu}}_{\mathit{\theta}}=\frac{\mathit{\pi}}{\mathrm{2}}$ (neglect skewness in fluid's relative velocity PDF) and obtain *K* by matching the transverse to longitudinal ratio of structure functions (TLR) of the particle relative velocities with the ones via the DNS data; *N* is determined via normalization of *P*(*θ*). Figure 2b shows the 〈*w*_{r}〉_{∗} derived via Eqs. (11) and (12), using *K* calibrated with TLR of second-, fourth-, and sixth-order structure functions respectively. The results have a correct qualitative trend of vanishing values at large *r* that increase sharply as *r* approaches *d*, with the fourth-order's result giving the best agreement with DNS. Currently we do not have a satisfactory rationale to single out the fourth order. The TLR of different orders giving differing results may imply that our first-order model may be incomplete, possibly due to over-simplification in Eq. (12) or to the inaccuracy of the rectilinear assumption ($d/\mathit{\eta}$ in the DNS may be insufficiently small).

## 4.4 Differential version of the theory, its validity and its solution

We now discuss an important but precarious theoretical issue. Chun et al. (2005) clearly showed that the non-local diffusion (${q}_{r}^{D}$) may only be converted from its general integral form into a differential version when the underlying RDF is a simple power law. However, Lu et al. (2010) and Yavuz et al. (2018), working in two very different scenarios, found that their predictions using the differential form of the theory agree well with experiments, even when the RDFs involved were clearly not power laws. We shall attempt to remedy this apparent paradox in future work. To examine how well this albeit unjustified method works here, we recast Eq. (9) into its differential form (Chun et al., 2005):

where *B*_{nl}=0.0397 (this value is computed from our DNS), and *B*_{nl} is expected to depend on flow characteristics, e.g., *R*_{λ} and *τ*_{η} (more in the Supplement).
Using Eq. (13), the same *g*_{s}*g*_{0} ansatz,
we make another prediction for 〈*W*_{r}〉_{St=0.054}, which is plotted in Fig. 2a (dashed green line).
This prediction is far from the DNS at *r*∼*d* but performs as well as the integral version at *r*≫*d* (the jump in the curve is just an artifact from the kink in the ansatz).

One advantage of Eq. (13) is that it allows for a general solution, which we now give, assuming 〈*W*_{r}〉 is given by Eqs. (11) and (12):

with $q\left(r\right)={R}_{c}^{*}{\mathit{\tau}}_{\mathit{\eta}}/\left({B}_{nl}{r}^{\mathrm{4}}\right)$, $\mathit{\beta}\left(r\right)=\mathrm{exp}\left[\int p\left(r\right)\mathrm{d}r\right]$ and $p\left(r\right)=\left[{A}_{\mathit{\tau}}r-\langle {w}_{r}{\rangle}_{\ast}\right]{\mathit{\tau}}_{\mathit{\eta}}/\left({B}_{nl}{r}^{\mathrm{2}}\right)$ (more in the Supplement).

## 4.5 Effects of gravity and other limitations

Thus far, we have not considered the effects of gravity on the particles. Here we provide a glimpse into the role of gravity (a detailed analysis is beyond the scope of the present study). In keeping with the scope of current work, we restrict ourselves to the case of monodisperse particles only. For this, we rerun the DNS cases of *St*=0.054 and 0.54 with gravity (to be compared with the zero-gravity case). The new particle advection equation is $\mathrm{d}\mathbf{v}/\mathrm{d}t=(\mathbf{u}-\mathbf{v})/{\mathit{\tau}}_{p}+\mathit{g}$ (all other details of the DNS remain unchanged). We choose to have the particle settling parameter ${S}_{g}\equiv {\mathit{\tau}}_{p}g/{u}_{\mathit{\eta}}$ (where *u*_{η} is the Kolmogorov velocity scale) be in the range *O*(0.1)−*O*(1)
(this is achieved by letting $\left|\mathit{g}\right|=\mathrm{10}$ dm s^{−2}). As a result, the range of *S*_{g} and *St* explored here is well aligned with measured values in natural clouds (Siebert et al., 2010). For the case of *St*=0.054 (*S*_{g}=0.49), we find no discernible difference for both RDF and MRV between the “with gravity” and zero-gravity results (corresponding figures in the Supplement). For the *St*=0.54 (*S*_{g}=4.9) case, Fig. 3 shows the effects of gravity on the RDF and MRV. We see that the slope (exponent of *g*(*r*) in the range $d\ll r<\mathrm{20}\mathit{\eta}$) of the RDF in the gravitational case is reduced by about 15 % compared to the zero-gravity case (*S*_{g}=0).
However, the shape of the RDF in the collision regime (*r*∼*d*) is approximately preserved, suggesting that a construct of the form *g*_{collision}×*g*_{gravity} may be a good first-order model for the full RDF (close examination of the compensated RDFs gives substantial support for this idea, details in the Supplement). These observations imply that as *S*_{g} increases from *O*(0.1) to *O*(1), the effects of gravity on RDF grow from negligible to significant but not dominant; the main effect is the reduction of the exponent while the collision-related “modulation” (*g*_{collision}) remains largely intact. The inset of Fig. 3 shows that the MRV is also weakened by gravity, though the statistical noise limits the strength of this conclusion. Lastly, it is worth noting that in the complementary DNS by Woittiez et al. (2009) that included gravity but not actual collisions, a much stronger gravitational effect was found on the statistics of bidisperse particles relative to the monodisperse case.

As mentioned, the fundamental focus of our work precludes the DNS and theory from considering a number of complexities relevant to some applications. As a result, this limits the direct quantitative applicability of our results to some realistic problems (e.g., in clouds). Besides gravity, another neglected factor is the hydrodynamic inter-particle force (HDI). Recent works, e.g., Yavuz et al. (2018); Bragg et al. (2022) found that the HDI also has a strong impact on RDF for *r*∼*d*. For monodisperse particles with small to moderate *St*, the HDI is expected to be more important than gravity. While we expect that the HDI should not alter the qualitative trend that *g*(*r*) should fall towards a small value at *r*→*d* (the same applies to the observed trend of MRV), it is likely that the HDI and collision would affect RDF and MRV in a coupled manner.

Also neglected is the influence of temperature, humidity and the vapor–liquid phase transition, which are important in the atmospheric clouds. These factors have a substantial impact on the polydispersity of small droplets (see, e.g., Kumar et al., 2012, 2014). However, for monodisperse statistics considered here, they are likely to play a minor role (they will be more important when future works consider the full polydisperse problem).

One limitation of the theory stems from the assumption of *St*≪1 and its corollary that particle velocity statistics in this regime are *St*-independent (Chun et al., 2005), which limits the theory's applicability to real systems. This implies that MRV should be *St*-independent in this regime. Our DNS results (spanning 2 orders of magnitude in *St*) shown in Fig. 2 give some support to the latter. However, unlike the theoretical prediction for MRV of case *St*=0.054 (Fig. 2), we have found that the prediction for *St*=0.11 is discernibly below the DNS result (figure in the Supplement). This could be due to the finite *St* effect not captured by the theory or could be for other reasons (details in the Supplement). Hence, a finite *St* extension of the theory is desirable to improve its applicability to real systems.

To conclude, we observed that collision strongly affects the RDF and MRV and imposes strong coupling between them^{6}. This challenges the efficacy of a “separation paradigm” and suggests that results from any studies that preclude particle collision have limited relevance for predicting collision statistics.
We have presented a theory for particle collision–coagulation in turbulence (based on a Fokker–Planck framework) that explains the above observations
and verified its accuracy by showing that 〈*W*_{r}〉 could be accurately predicted using a sufficiently accurate RDF. The theory accounts for the full collision–coagulation rate, which includes contributions from mean field and fluctuations, and as such,
our work complements and completes earlier mean-field theories (Saffman and Turner, 1956; Sundaram and Collins, 1997). We showed that a simple model of particle approach angles could capture the main features of 〈*W*_{r}〉 and used it to derive a general solution for RDF from the differential version of the theory. We uncovered a possible paradox regarding the past empirical successes of the differential drift-diffusion equation (see Sect. 4.4).
Further shape-preserving reduction of the RDF and MRV was observed when the gravitational settling parameter (*S*_{g}) is of order *O*(1).
Our findings provide new perspectives of particle collision and its relation with clustering and relative motion, which have implications for atmospheric clouds or generally for systems involving colliding particles in unsteady flows.

The data (model codes) may be obtained from the corresponding author upon request.

The supplement related to this article is available online at: https://doi.org/10.5194/acp-22-3779-2022-supplement.

EWS oversaw the conception and execution of the project. EWS did the theoretical derivations in collaboration with XM. XM and EWS conducted the numerical simulation and data analysis. EWS and XM wrote the article.

The contact author has declared that neither they nor their co-author has any competing interests.

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

This work was supported by the National Natural Science Foundation of China (grant 11872382) and by the Thousand Young Talents Program of China. We thank Jialei Song for their help. We thank Wai Chi Cheng, Jianhua Lv, Liubin Pan and Raymond A. Shaw for discussion and suggestions.

This research has been supported by the National Natural Science Foundation of China (grant no. 11872382).

This paper was edited by Graham Feingold and reviewed by two anonymous referees.

Balkovsky, E., Falkovich, G., and Fouxon, A.: Intermittent Distribution of Inertial Particles in Turbulent Flows, Phys. Rev. Lett., 86, 2790, 2001. a

Bec, J., Biferale, L., Cencini, M., Lanotte, A., Musacchio, S., and Toschi, F.: Heavy Particle Concentration in Turbulence at Dissipative and Inertial Scales, Phys. Rev. Lett., 98, 084502, 2007. a, b

Bec, J., Ray, S. S., Saw, E. W., and Homann, H.: Abrupt growth of large aggregates by correlated coalescences in turbulent flow, Phys. Rev. E, 93, 031102, 2016. a

Bragg, A. D., Hammond, A. L., Dhariwal, R., and Meng, H.: Hydrodynamic interactions and extreme particle clustering in turbulence, J. Fluid Mech., 933, A31, https://doi.org/10.1017/jfm.2021.1099, 2022. a

Chun, J., Koch, D. L., Rani, S. L., Ahluwalia, A., and Collins, L. R.: Clustering of aerosol particles in isotropic turbulence, J. Fluid Mech., 536, 219–251, 2005. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x

Dou, Z., Bragg, A. D., Hammond, A. L., Liang, Z., Collins, L. R., and Meng, H.: Effects of Reynolds number and Stokes number on particle-pair relative velocity in isotropic turbulence: a systematic experimental study, J. Fluid Mech., 839, 271–292, 2018. a

Eswaran, V. and Pope, S. B.: An examination of forcing in direct numerical simulations of turbulence, Comput. Fluids, 16, 257–278, 1988. a

Falkovich, G., Fouxon, A., and Stepanov, M. G.: Acceleration of rain initiation by cloud turbulence, Nature, 419, 151, 2002. a

Grabowski, W. W. and Wang, L.-P.: Growth of cloud droplets in a turbulent environment, Annu. Rev. Fluid Mech., 45, 293–324, 2013. a

Ireland, P. J., Vaithianathan, T., Sukheswalla, P. S., Ray, B., and Collins, L. R.: Highly parallel particle-laden flow solver for turbulence research, Comput. Fluids, 76, 170–177, 2013. a

Johansen, A., Oishi, J. S., Mac Low, M.-M., Klahr, H., Henning, T., and Youdin, A.: Rapid planetesimal formation in turbulent circumstellar disks, Nature, 448, 1022–1025, 2007. a

Karnik, A. U. and Shrimpton, J. S.: Mitigation of preferential concentration of small inertial particles in stationary isotropic turbulence using electrical and gravitational body forces, Phys. Fluids, 24, 073301, 2012. a

Karpińska, K., Bodenschatz, J. F. E., Malinowski, S. P., Nowak, J. L., Risius, S., Schmeissner, T., Shaw, R. A., Siebert, H., Xi, H., Xu, H., and Bodenschatz, E.: Turbulence-induced cloud voids: observation and interpretation, Atmos. Chem. Phys., 19, 4991–5003, https://doi.org/10.5194/acp-19-4991-2019, 2019. a

Kolmogorov, A. N.: The local structure of turbulence in incompressible viscous fluid for very large Reynolds numbers, Dokl. Akad. Nauk SSSR+, 30, 299–303, 1941. a, b

Kumar, B., Janetzko, F., Schumacher, J., and Shaw, R. A.: Extreme responses of a coupled scalar–particle system during turbulent mixing, New J. Phys., 14, 115020, 2012. a

Kumar, B., Schumacher, J., and Shaw, R. A.: Lagrangian mixing dynamics at the cloudy–clear air interface, J. Atmos. Sci., 71, 2564–2580, 2014. a

Lu, J., Nordsiek, H., Saw, E. W., and Shaw, R. A.: Clustering of Charged Inertial Particles in Turbulence, Phys. Rev. Lett., 104, 184505, 2010. a, b

Maxey, M. R. and Riley, J. J.: Equation of motion for a small rigid sphere in a nonuniform flow, Phys. Fluids, 26, 883–889, 1983. a

Mortensen, M. and Langtangen, H. P.: High performance Python for direct numerical simulations of turbulent flows, Comput. Phys. Commun., 203, 53–65, 2016. a, b

Onishi, R. and Seifert, A.: Reynolds-number dependence of turbulence enhancement on collision growth, Atmos. Chem. Phys., 16, 12441–12455, https://doi.org/10.5194/acp-16-12441-2016, 2016. a

Pope, S. B.: Turbulent Flows, Cambridge Univ. Press, Cambridge, UK, https://doi.org/10.1017/CBO9780511840531, 2000. a, b, c, d, e

Reade, W. C. and Collins, L. R.: Effect of preferential concentration on turbulent collision rates, Phys. Fluids, 12, 2530, https://doi.org/10.1063/1.1288515, 2000. a

Rogallo, R. S.: Numerical experiments in homogeneous turbulence, vol. 81315, National Aeronautics and Space Administration, 1981. a, b, c

Saffman, P. and Turner, J.: On the collision of drops in turbulent clouds, J. Fluid Mech., 1, 16–30, 1956. a, b, c, d

Salazar, J. P. L. C., de Jong, J., Cao, L., Woodward, S. H., Meng, H., and Collins, L. R.: Experimental and numerical investigation of inertial particle clustering in isotropic turbulence, J. Fluid Mech., 600, 245–256, 2008. a, b

Saw, E. W., Shaw, R. A., Ayyalasomayajula, S., Chuang, P. Y., and Gylfason, A.: Inertial Clustering of Particles in High-Reynolds-Number Turbulence, Phys. Rev. Lett., 100, 214501, 2008. a, b, c

Saw, E.-W., Shaw, R. A., Salazar, J. P., and Collins, L. R.: Spatial clustering of polydisperse inertial particles in turbulence: II. Comparing simulation with experiment, New J. Phys., 14, 105031, 2012. a, b

Saw, E.-W., Bewley, G. P., Bodenschatz, E., Sankar Ray, S., and Bec, J.: Extreme fluctuations of the relative velocities between droplets in turbulent airflow, Phys. Fluids, 26, 111702, 2014. a

Siebert, H., Gerashchenko, S., Gylfason, A., Lehmann, K., Collins, L., Shaw, R., and Warhaft, Z.: Towards understanding the role of turbulence on droplets in clouds: in situ and laboratory measurements, Atmos. Res., 97, 426–437, 2010. a

Sundaram, S. and Collins, L.: Collision statistics in an isotropic particle-laden turbulent suspension. Part 1. Direct numerical simulations, J. Fluid Mech., 335, 75–109, 1997. a, b, c, d, e, f, g

Voßkuhle, M., Lévêque, E., Wilkinson, M., and Pumir, A.: Multiple collisions in turbulent flows, Phys. Rev. E, 88, 063008, 2013. a, b

Wang, L.-P., Ayala, O., Rosa, B., and Grabowski, W. W.: Turbulent collision efficiency of heavy particles relevant to cloud droplets, New J. Phys., 10, 075013, https://doi.org/10.1088/1367-2630/10/7/075013, 2008. a, b

Wilkinson, M., Mehlig, B., and Bezuglyy, V.: Caustic activation of rain showers, Phys. Rev. Lett., 97, 48501, 2006. a

Woittiez, E. J., Jonker, H. J., and Portela, L. M.: On the combined effects of turbulence and gravity on droplet collisions in clouds: a numerical study, J. Atmos. Sci., 66, 1926–1943, 2009. a, b

Wood, A. M., Hwang, W., and Eaton, J. K.: Preferential concentration of particles in homogeneous and isotropic turbulence, Int. J. Multiphas. Flow, 31, 1220, 2005. a

Yavuz, M., Kunnen, R., Van Heijst, G., and Clercx, H.: Extreme small-scale clustering of droplets in turbulence driven by hydrodynamic interactions, Phys. Rev. Lett., 120, 244504, 2018. a, b, c

Zaichik, L. I. and Alipchenkov, V. M.: Pair dispersion and preferential concentration of particles in isotropic turbulence, Phys. Fluids, 15, 1776, 2003. a

^{1}

Coagulation is, in a sense, the simplest outcome of collision. In the companion paper we shall argue that the major qualitative conclusions of our work also apply to cases with other collisional outcomes.

^{2}

In the collision-less case, the asymmetry is much weaker and is related to viscous dissipation of energy in turbulence (Pope, 2000).

^{3}

In other words, particles may approach each other (and collide), but they can not be created at contact and then separate.

^{4}

Here a positive ${q}_{r}^{d}$ merely reflects a deficit in the inward flux of neighboring particles since we find that ${q}_{r}^{d}+{q}_{r}^{D}$ is always negative.

^{5}

This is true in the relatively idealized system simulated but may not apply to the general problem that includes other effects

^{6}

This statement also holds for other types of collisional outcomes (not only for collision–coagulation), but the details of the specific outcomes should be different from the current case.