# Intermittency and collisions of fast sedimenting droplets in turbulence

Itzhak Fouxon1<EMAIL_ADDRESS>Seulgi Lee2<EMAIL_ADDRESS>Changhoon
Lee1,2<EMAIL_ADDRESS>1 School of Mathematics and Computing, Yonsei
University, Seoul 03722, South Korea  2 Department of Mechanical Engineering,
Yonsei University, Seoul 03722, South Korea

###### Abstract

We study theoretically and numerically spatial distribution and collision rate
of droplets that sediment in homogeneous isotropic Navier-Stokes turbulence.
It is assumed that, as it often happens in clouds, typical turbulent
accelerations of fluid particles are much smaller than gravity. This was shown
to imply that the particles interact weakly with individual vortices and, as a
result, form a smooth flow in most of the space. In weakly intermittent
turbulence with moderate Reynolds number $Re_{\lambda}$, rare regions where
the flow breaks down can be neglected in the calculation of space averaged
rate of droplet collisions. However, increase of $Re_{\lambda}$ increases
probability of rare, large quiescent vortices whose long coherent interaction
with the particles destroys the flow. Thus at higher $Re_{\lambda}$, that
apparently include those in the clouds, the space averaged collision rate
forms in rare regions where the assumption of smooth flow breaks down. This
intermittency of collisions implies that rain initiation could be a strongly
non-uniform process. We describe the transition between the regimes and
provide collision kernel in the case of moderate $Re_{\lambda}$ describable by
the flow. The distribution of pairwise distances (radial distribution function
or RDF) is shown to obey a separable dependence on the magnitude and the polar
angle of the separation vector. Magnitude dependence obeys a power-law with a
negative exponent, manifesting multifractality of the droplets’ attractor in
space. We provide the so far missing numerical confirmation of a relation
between this exponent and the Lyapunov exponents and demonstrate that it holds
beyond the theoretical range. The angular dependence of the RDF exhibits a
maximum at small angles quantifying particles’ formation of spatial columns.
We provide typical dimensions of the columns, which belong in the inertial
range. We derive the droplets’ collision kernel using that in the considered
limit the gradients of droplets’ flow are Gaussian. We demonstrate that as
$Re_{\lambda}$ increases the columns’ aspect ratio decreases, eventually
becoming one when the isotropy is restored. We propose how the theory could be
constructed at higher $Re_{\lambda}$ of clouds by using the example of the
RDF.

##  I Introduction

Turbulence of air in warm clouds accelerates collisions of water droplets and
thus must be included in studies of precipitation [1, 2, 3, 4, 5, 6]. This
inclusion is of high interest since it could help to resolve the bottleneck
problem in rain formation. The bottleneck is caused by the narrowness of the
size distribution of droplets created by condensation of vapor. The size
proximity implies that the difference in settling velocities of the droplets
is so small that their gravitational collisions would take long times, that
are incompatible with the observations. Hence, the passage to the later stages
of formation of large precipitating drops, where the growth occurs by
gravitational collisions of different size droplets, demands something beyond
gravity. Turbulence, that creates size dispersion by introducing mechanism for
collisions of equal-size droplets, can be the missing factor. However, the
relevance of turbulence has not been quantified so far [7] and assessment of
how much turbulence influences the rain formation is an open problem. A main
problem is the extremely high Reynolds number that holds in the clouds. Thus
numerical studies can be performed only at much smaller Taylor microscale
Reynolds numbers $Re_{\lambda}$ than in the clouds. This makes rigorous
theoretical studies of droplets’ collisions in the Navier-Stokes turbulence
(NST), which is the focus of this paper, specially significant. However
theoretical studies of droplets’ behavior in high-$Re_{\lambda}$ NST are
obstructed by having to deal with flow whose statistics is unknown. The
problem is aggravated because intermittency produces significant probability
of rare events that may locally accelerate the collision rates by a large
factor in comparison with estimates using typical events.

The lack of knowledge of statistics of turbulence by itself does not prohibit
quantitative predictions for the turbulent transport. In fact, the theory is
able to provide accurate quantitative predictions for the NST [4]. This is
thanks to independence of certain properties of turbulent transport of the
details of the statistics (universality). As it is often the case in
statistical physics, the universality holds due to the appearance of sum of a
large number of independent random variables in the analysis. Here we study
the case where the particle sediments so quickly through an individual vortex
that the vortex perturbs its motion only weakly. Thus the droplet’s velocity
is determined by effects of many independent vortices accumulated during the
particle velocity relaxation time $\tau$. The issuing universality demands
that intermittency is not too strong so that the probability of large
quiescent vortices is not too high (intermittency both increases the
probability of strong bursts and of large quiescent regions [8]). Otherwise
individual vortex is so large that the particle never leaves it during the
interaction time $\tau$. The motion would then be determined by interaction
with a single vortex destroying the universality. This limits the study to
moderate $Re_{\lambda}$ well below those in clouds where the statistics is
very intermittent. Still we demonstrate that the results provide a useful
reference point and some of them do generalize to high $Re_{\lambda}$.

As far as the case of weak intermittency is concerned, the present work
continues the study [9] which provided successful quantitative predictions for
correlation dimension and Lyapunov exponents of particles in the Navier-Stokes
turbulence. The predictions were confirmed by direct numerical simulations at
$Re_{\lambda}=70$ (claimed in [10] contradiction with some data is eliminated
below). Here we provide the theory of the collision kernel and perform
simulations at the same $Re_{\lambda}=70$, that demonstrate full agreement
with the theory. We derive angle-dependent radial distribution function (RDF)
that gives probability of finding a pair of particles at a given vector
separation. The RDF determines the collision kernel by providing probability
of a pair at collision distance. The angle dependence of the RDF allows us to
obtain vertical and horizontal dimensions of particles’ columns in space whose
formation is the signature of the studied fast sedimentation limit, see Fig.
1. In the remainder of the Introduction we provide outlook (currently lacking
in the literature) at the existing rigorous theories and describe the paper’s
organization.

###  I.1 Theory at weak gravity and small inertia and its breakdown in clouds

Until [9], a complete rigorous theory for inertial particles in the NST
existed only at negligible gravity and small but finite inertia [4, 11]. The
magnitude of inertia is measured by the dimensionless Stokes number $St$ which
is the ratio of $\tau$ and the typical turnover-time of the viscous scale
vortices, the Kolmogorov time $\tau_{K}$, see [8]. The inertia is usually
considered to be small at $St\ll 1$ where the particle velocity would relax
quickly to the local velocity of the flow and the particles would nearly trace
the flow [12]. However we demonstrate in this subsection that intermittency of
turbulence quite certainly makes the effect of inertia in clouds large, even
at $St\ll 1$.

The work [4] introduced a picture of motion of particles in space that holds
at $St\ll 1$ independently of the Reynolds number and, thus, would hold also
in the clouds. This picture can be understood by using the concept of the
local Stokes number, similar to that of the local Reynolds number [8]. This
number is necessary in order to describe strong contrasts in the strength of
particle-vortex interaction throughout a high Reynolds number flow. The
contrasts can be described via the local energy dissipation rate $\epsilon$,
defined as the product of the kinematic viscosity $\nu$ and the square of
velocity gradients. The rate undergoes strong intermittent fluctuations in
clouds, see [1] and references therein. We observe that the Kolmogorov time,
defining $St$, is given by $\tau_{K}\equiv\sqrt{\nu/\epsilon_{0}}$ where
$\epsilon_{0}$ is the mean energy dissipation. Correspondingly the local
relevance of inertia can be characterised by the local Stokes number
$\tau\sqrt{\epsilon/\nu}$. If the local Stokes number is not small, then
strong interaction of particles with the local vortices occurs in that region.
Due to the intermittency of turbulence, the magnitude of fluctuations of the
local Stokes number depends on $Re_{\lambda}$, implying dependence of
particles’ statistics on both $St$ and $Re_{\lambda}$.

At $St\ll 1$, for the more probable events turbulent vortices are slow,
$\tau\sqrt{\epsilon/\nu}\ll 1$. The particles’ motion in these vortices is
smooth and ordered in space so that there is no crossing of trajectories of
different particles. A flow of particles $\bm{v}(t,\bm{x})$ can be introduced
[12], providing the changes of their coordinates $\bm{x}(t)$ as
$\dot{\bm{x}}=\bm{v}(t,\bm{x}(t))$. However for fast vortices with time-scale
$\lesssim\tau$, characterized by non-small local Stokes number, the turbulent
driving creates jets of particles that separate from the turbulent flow, as in
a sling, and cause the trajectories’ crossing [4]. This “sling effect” was
observed experimentally [17] and at any $Re_{\lambda}$, including those in
clouds, it is confined to well isolated small regions of fast vortices.
Increase of $Re_{\lambda}$, at least in the high-$Re_{\lambda}$ limit, would
make regions of the sling effect even more rare in space because the regions
of quiescent turbulence increase in size with the Reynolds number [8]. This
creates significant difficulties in experimental and numerical measurements
which must have resolution that increases with $Re_{\lambda}$.

The above picture implies that space averages, such as the average rate of
droplets’ collisions, can be found as sums of contributions of flow and sling
regions (termed below ”flow contribution” and ”sling contribution”). The
contributions of these regions into the collision rate were calculated in [4].
It was demonstrated that, despite that the regions of the sling effect are
rare in space, their contribution into the collision kernel can be significant
because they create optimal conditions for collisions [4, 18].

In contrast with the above qualitative picture, that is valid at any Reynolds
number, quantitative predictions of [4, 11] break down at increasing
$Re_{\lambda}$. The predictions for the RDF rely on the assumption that at
$St\ll 1$ the RDF is determined by vortices whose turnover time, similarly to
the most probable turnover time $\tau_{K}$, is much smaller than $\tau$.
Making that assumption, the theory implies that distribution of inertial
particles in the turbulent flow is multifractal. Multifractality manifests
itself in the RDF that obeys at small distances a power-law with an exponent
$-\alpha$. Here $\alpha$ is positive \- the RDF diverges at zero separation
because the particles’ concentration is singular. An explicit formula for
$\alpha$ via a high order moment of the turbulent velocity gradients is
provided in [4, 11]. Application of the standard phenomenology of turbulence
to that formula gives that at large $Re_{\lambda}$ the magnitude of the
exponent is proportional to the product of $St^{2}$ and a positive power of
Reynolds number. This implies that in the limit of large $Re_{\lambda}$,
however small $St$ is, the exponent, calculated under the described
assumptions, is larger than $3$. This however implies that the RDF becomes
non-integrable at the origin which contradicts finiteness of the total number
of particles (the number’s second moment is proportional to the integral of
the RDF over distance).

The reason for the described contradiction is that the assumption of the
theory that the RDF is determined by vortices with turnover time smaller than
$\tau$ becomes inconsistent at large $Re_{\lambda}$. In fact, the theory
predicts its own breakdown via the formula for $\alpha$. Intermittency of
turbulence implies appreciable presence in the flow of vortices whose turnover
time is much smaller than $\tau_{K}$. More precisely the probability that the
turnover time is smaller than $\tau_{K}$ by a power of the Reynolds number is
appreciable [8], cf. the recent [13]. As a result, in the limit of large
Reynolds numbers the gradients that determine $\alpha$ grow with
$Re_{\lambda}$ according to a non-trivial power-law [16]. This implies that
the time-scale of the vortices that determine $\alpha$ decreases with
$Re_{\lambda}$ as $\tau_{K}Re_{\lambda}^{-q}$ with $q>0$ (according to the
standard phenomenology time scales as the inverse gradient). Thus, considering
increase of $Re_{\lambda}$ at a fixed small $St$, at large enough
$Re_{\lambda}$ the time-scale of the relevant vortices becomes of order
$\tau$. We designate the corresponding threshold Reynolds number by
$Re_{\lambda}^{*}$. This number is not sharp and is determined by order of
magnitude only via the asymptotic equality
$\tau_{K}\left(Re_{\lambda}^{*}\right)^{-q}\sim\tau$. Thus at
$Re_{\lambda}\gtrsim Re_{\lambda}^{*}$ the assumption that the relevant
vortices have time-scale larger than $\tau$ is inconsistent. We conclude that
the RDF is determined at $Re_{\lambda}\gtrsim Re_{\lambda}^{*}$ by vortices
whose characteristic turnover time is of order $\tau$ or less. Anomalously
strong, short-lived vortices become relevant to the RDF at higher
$Re_{\lambda}$ because intermittency makes their probability appreciable. We
also conclude that an approach other than that of [4, 11] must be devised at
$Re_{\lambda}\gtrsim Re_{\lambda}^{*}$.

The above considerations demonstrate that the actual validity condition of the
small inertia theory of [4, 11] is $St\ll Re^{-q}$ and not $St\ll 1$, cf. [14,
15]. The value of $q$ can hardly be predicted and must be fixed numerically.
Thus the theory of [4, 11] is a so-called asymptotic theory at small $St$ i.e.
it holds at $St\to 0$ when $Re_{\lambda}$ is fixed.

The above demonstration of limitations of the theory of [4] due to
intermittency is performed by using the laws holding in the limit of large
$Re_{\lambda}$. This does not signify that intermittency becomes relevant only
at $Re_{\lambda}\to\infty$. Simulations of [14] of the particles’ motion in
the NST demonstrated that for $St\sim 0.1$ the small inertia theory of the RDF
of [4] applies accurately at the moderate $Re_{\lambda}=21$. However,
considering higher $Re_{\lambda}$, already at $Re_{\lambda}=47$, that theory
fails. We see that smallness of $St$ does not guarantee the smallness of
inertial effects already at rather moderate $Re_{\lambda}$. Despite that the
criterion of validity of the theory, $St\ll Re^{-q}$, derived at
$Re_{\lambda}\to\infty$, cannot be used at $Re_{\lambda}$ of [14], we can use
it to get very rough idea of magnitude of $q$. We find that $q\sim 0.5$ i.e.
it is a number of order one. Using $q=0.5$ we find that for $Re_{\lambda}\sim
10^{4}$, that can hold in the clouds, the effects of inertia on the RDF are
appreciable at $St$ as small as $0.01$ (which covers the whole range of
relevant droplet sizes, see Sec. II). More precise considerations demand
extensive future numerical work.

The above conclusion transfers to the collision kernel which, under the
assumption of small inertia, is given by a sum of two terms [4]. One term
describes the contribution of vortices whose turnover time is much smaller
than $\tau$. That term is proportional to the RDF and, as we saw, it becomes
determined by vortices with turnover time of order $\tau$ at higher
$Re_{\lambda}$ . The other term in the collision kernel describes the
contribution of the sling effect which by definition is due to vortices with
time-scale $\tau$ or smaller. Thus in the limit of high $Re_{\lambda}$ the
collision kernel is due to vortices with turnover time $\tau$ or smaller.

We summarize how increase of $Re_{\lambda}$ results in the breakdown of the
theory of [4]. The theory demonstrates that the clustering rate, holding in
the regions of the smooth motion of the droplets, is determined by vortices
whose characteristic time-scale decreases with $Re_{\lambda}$ as a power-law.
The theory holds at moderate $Re_{\lambda}$ where this time-scale is much
larger than $\tau$ as necessary for the self-consistency of the assumption of
smooth motion. When $Re_{\lambda}$ increases, faster vortices with decreasing
turnover time become relevant, until their time-scale becomes comparable with
$\tau$ at some $Re^{*}_{\lambda}$, cf. [14]. When this happens the separation
of the total rate of collisions into the contributions of the regions of the
smooth motion and of the sling effect, that was used in [4] for the
calculations, breaks down. Both contributions are determined by ”resonant”
vortices with lifetime of order $\tau$. The corresponding changes in the
theory will be published elsewhere.

###  I.2 Fast sedimentation theory

Applicability of the above theory to clouds is limited by three factors. The
first factor is the Stokes number which is not necessarily small for the
droplets taking part in the rain formation. The second factor is gravity,
which in clouds is strong and not weak. The magnitude of gravity can be
measured by the Froude number $Fr$ given by the ratio of the typical
(Lagrangian) acceleration of the fluid particles
$\epsilon_{0}^{3/4}/\nu^{1/4}$ and gravitational acceleration $g$, i.e.
$Fr\equiv\epsilon_{0}^{3/4}/(g\nu^{1/4})$. Thus, [4] applies at $Fr\to\infty$
when gravity is negligible. In contrast, in the clouds even rather strong
cloud turbulence with $\epsilon_{0}=2000cm^{2}/s^{3}$ gives $Fr=0.5$. For
typical turbulence with smaller $\epsilon_{0}$, the Froude number is yet
smaller. The last limiting factor is the high $Re_{\lambda}$ of the clouds, as
described above. The two former limitations were overcome in [9] who
constructed the theory that holds at strong gravity, $Fr\to 0$, and any $St$
including $St\gg 1$. However this theory, similarly to that considered in the
previous subsection, is also asymptotic and has the validity condition $Fr\ll
Re^{-q^{\prime}}$ with $q^{\prime}>0$, as will be discussed in this work. Thus
how small $Fr$ must be, for the theory to be valid at a given $Re_{\lambda}$,
is not known. It was found in numerical simulations of [9] that at
$Re_{\lambda}=70$ the validity of the theory demands $Fr\leq 0.033$ which
implies quite weak typical accelerations $\epsilon_{0}^{3/4}/\nu^{1/4}$ as
compared with the gravity. The theory was also shown to apply reasonably well
at $Fr=0.05$.

The main limitation of [9] is again the Reynolds number - the intermittency
increases the fluctuations of the gradients and destroys the approximations
made. Nevertheless this theory, in comparison with that of [4], incorporates
the effects of the gravitational sedimentation of the droplets, that are
strong in precipitating clouds, giving us qualitative insight and the
possibility of consistent interpolation. It seems reasonable, given that
clouds are indeed characterized by small $Fr$, to make the small $Fr$ theory a
starting point for approaching collisions in clouds. It is to this theory that
this paper is devoted. Here we provide more detailed predictions than in [9],
confirm them numerically and describe how the theory breaks down at higher
Reynolds number.

Before continuing the development of the approach of [9], a controversy must
be faced. The theoretical formula of [9] for the correlation codimension
$\alpha$, studied in detail below, was tested numerically in [10] at
$Re_{\lambda}=398$. The comparison was done for $Fr=0.052$ and a range of
$St$. They observed that while [9] describe correctly the independence of
$\alpha$ of $St$ for $St\geq 1$, quantitatively the predictions are wrong by
about $50$ per cent. It is unfortunate that the performed comparison contained
two mistakes. The equation $(4.21)$ of [10] which the authors considered as
the prediction of [9] for $\alpha$, is factor of $2$ smaller than the actual
prediction made in [9]. If the correct formula is used then the discrepancy is
about $30$ and not $50$ per cent. This seems to be as much as one can hope
for, because [10] use the predictions outside the domain of validity of the
theory, which is $\alpha\ll 1$, cf. [11]. In fact, the result of [10]
completely agrees with the simulations of [9] that demonstrated that at
$Fr=0.05$ there are significant deviations from the theory. The inequalities
$\alpha\ll 1$ and $Fr\ll 1$ differ much because of the large numerical factor:
we have $\alpha\approx 13Fr$, see below and Sec. IV.

We stress the asymptotic character of the described theories in order to avoid
future misunderstandings. Both theories of [4] and [9] hold rigorously in the
limits of $St\to 0$ and $Fr\to 0$, respectively, when the Reynolds number is
held fixed. Thus, if it is found that the predictions of [9] are invalid, then
it tells that $Fr$ is too large and by decreasing $Fr$ the theory will be made
to hold true. Similar situation holds for $St\to 0$ theory of [4]. Thus the
observation of [14] of breakdown of the theory of [4] at $Re_{\lambda}\geq 47$
only tells that at this $Re_{\lambda}$ the theory applies at $St$ smaller than
those considered in [14]. These $St$ are determined by the condition that
their corresponding $\tau$ are much smaller than the lifetime of the vortices
that determine the sum of the Lyapunov exponents in [4]. As an example of the
use of the asymptotic theory beyond its region of validity we will demonstrate
below that the numbers obtained from $Fr\to 0$ theory at $Re_{\lambda}=70$ in
[9] can be used for predicting the values observed at $Re_{\lambda}=398$ in
[10].

###  I.3 Review of developments before the work [9]

We review the developments leading to the small $Fr$ theory of [9], both to
give credit and address concerns of an anonymous referee. This theory was
aimed to explain the observations of [19] of patterns of particles sedimenting
in turbulence. Flow description of the particle motion was employed, which is
rigorously valid at $Fr\to 0$ at other parameters, including the Reynolds
number, fixed, cf. above. The earliest clue to the possibility of the flow
description seemingly was made in [15] who observed that increasing gravity at
a fixed, not necessarily small $St$, damps the sling effect. The reason is
that faster sedimentation of the particle shortens the interaction time during
which an individual vortex swings the particle before the shooting. Thus the
sling events, that destroy single-valuedness of the flow [4], become more rare
at increasing gravity and at a sufficiently strong gravity the flow may become
single-valued despite a possibly large $St$.

The value of gravity at which a smooth single-valued flow of droplets exists
in the NST was provided in a Master thesis [20]. The study was done
disregarding the effects of intermittency and predicted that at $Fr\ll 1$ the
flow is well-defined for typical vortices, independently of $St$. In
application to clouds this implies that all droplets with size within the
range from $10$ to $60$ microns, which is the range where turbulence has
relevance, move in most of the space according to a size-dependent smooth
flow. The condition, that the regions where the smooth flow description breaks
down are rare, happens to imply that the flow’s compressibility is necessarily
weak (this coincidence occurs because the same property of inertia that causes
the sling effect also causes the compressibility of the particles’ flow). This
allows to derive the droplets’ distribution from the general solution for
tracers’ distributions in a weakly compressible random flow [11, 4].

Later, the observation that the sling effect is deactivated at large gravity
was done in a model two-dimensional flow studied in the regime where the
turbulent flow constitutes a small perturbation of the particle’s trajectory
in [21], see also [22, 19] and cf. similar expansion in [12]. All these works
but [22] used for the study of the sling effect the blowup equation introduced
in [4]. The work [22] provided a different outlook by observing numerically
that at $Fr\to 0$ average velocity difference of nearby particles scales
linearly with the distance between these particles. Thus [22] concluded that
in this limit the particles become tracers in an effective flow. Care is
needed though. Indeed, if there is a flow, then the linear scaling holds for
each realization and hence also statistically. However linear scaling of the
average velocity difference would also hold for a velocity field with
(effective) discontinuities as in Burgers turbulence [23]. Thus the observed
linear scaling can be used as an indication of the existence of the flow
however not as its proof. Finally [9], who worked independently of [21, 22]
(the first arxiv version of [9] appeared in the same year with statement of
independent work), gave a rigorous theory of the droplets’ flow able to
provide quantitative predictions for the particles’ behavior in the NST.
Direct numerical simulations (DNS) of the NST were done for weakly
intermittent turbulence with $Re_{\lambda}=70$ and confirmed the theory. The
work demonstrated how complete calculations can be performed, despite that
there is no explicit formula for the droplets’ flow in terms of the underlying
fluid flow (the relationship between these flows is non-local both in space
and in time).

The limit of large gravity studied in [9] assumes that the droplet’s settling
velocity is larger than the typical turbulent velocity at the Kolmogorov scale
[8]. However the settling velocity must still be smaller than the integral
scale velocity, to fit the applications in clouds. Thus turbulence is not a
small perturbation of the particles’ trajectories, as in [21] or some
qualitative considerations of [22], because the particle velocity coincides
with the local flow velocity in the leading order. The work [9] differed from
previous works by aiming at a complete, rigorous $Fr\to 0$ theory for
particles in the NST without modelling assumptions. This comes as the next
effort to get a rigorous theory for the NST, after the $St\to 0$ theory of
[4]. For the considered $Re_{\lambda}=70$ it was observed numerically that the
theory is accurate for $Fr\leq 0.03$, with some theoretical predictions
holding up to $Fr=0.1$, cf. above.

Certain aspects of the theory of [9] were observed previously. Due to fast
sedimentation the separation of particle pairs is driven effectively by white
noise, as it was observed for the NST in [24]. The authors originally
considered the case of zero gravity and large inertia, $St\gg 1$. The large
inertia causes the particle to drift fast through the fluid so that the flow
”looks to it” as a white noise. However, [24] observe that, using the
gravitational drift instead of inertial one, gives the answer for the case
with gravity. The value of the first Lyapunov exponent for the case with
gravity was provided. It can also be seen from this work that pair separation
is effectively horizontal, the result which was significantly developed
further in [22]. Horizontality of separation implies that the particles spend
more time when located one above the other which will be seen as columns in
space observed in [22, 19]. Independently, the applicability of the white
noise model, and preferential alignment of the vector inter-pair distance with
the vertical, were observed in the model flow of [21], who also derived
numerically the Lyapunov exponents of their model. The observation that
gravity enhances preferential concentration at $St\gtrsim 1$ and can result in
strong particle clustering was done in [21, 22, 19] independently. In
contrast, at $St\ll 1$, gravity decreases the clustering where three different
asymptotic regimes exist [9].

###  I.4 Organization of the paper

All studies in this paper are performed for the Navier-Stokes equations of
incompressible flow and no model of turbulence is used. Since the text is
rather long then for reader’s convenience we describe organization of the
material. In the next section we introduce the equation of motion, discussing
its applicability in the high $Re_{\lambda}$ case of clouds, and the numerical
scheme. The main message of this section for future studies is that in high
$Re_{\lambda}$ turbulence we need to use a more fundamental description of the
particle-vortex interaction than in the usually used equation of particle
motion, in order to be certain that we do not miss significant contributions
in the collision kernel. This is because the flow fluctuations with scales
smaller than the particle size might be relevant due to intermittency.

In section III we describe how increase of $Re_{\lambda}$ invalidates small
$Fr$ theories in the simplest context. We derive the remarkably simple
structure of separation of close particles at $Fr\to 0$: the vertical
component of the separation is conserved and the horizontal component evolves
as separation of two tracers in a white-noise velocity known as the Kraichnan
model [36]. This structure is the reason for columns’ formation. We derive the
spectrum of the Lyapunov exponents extending the results of [9]. We observe
that horizontal motions are driven by vortices much larger than the Kolmogorov
scale whose characteristic size increases with $Re_{\lambda}$. This leads to
breakdown of the theory above certain $Re_{\lambda}$, where the statistics
becomes isotropic and particles would not form columns in space.

The next section is devoted to reviewing the developments that preceded the
theory of [9] and description of the main predictions of that theory. It is
demonstrated, in a significantly more detailed way than in [9], that the DNS
of the NST at $Re_{\lambda}=70$, indicate unequivocally that [9] provides us
with a completely valid theory that can be used in the domain of its validity.
The contradiction with simulations of [10] at $Re_{\lambda}=398$ is due to the
application of the theory outside the domain of validity (this is besides that
[10] use a wrong numerical factor in studying the prediction of [9]). Section
V is devoted to explanation that the theory of [9] breaks down at increasing
$Re_{\lambda}$. The reason is growing intermittency that implies both
increasing regions of calm turbulence, allowing long coherent particle-eddy
interactions, and increasing relevance of strong bursts of velocity gradients.
It must be stressed again that the theory’d breakdown does not mean that there
is some critical $Re_{\lambda}$ where it stops to work. Rather, it says that
the range of small $Fr$ where [9] holds, shrinks to zero with increasing
$Re_{\lambda}$ in a power-law fashion.

The rest of the paper is devoted to the case of $Fr$ so small that the flow
description of [9] holds. We provide complete theory of collisions in this
limit and its numerical confirmation. Section VI derives the RDF at not too
small angles and same size particles. The main difference from [9] is the
recognition of the fact that the smoothness scale of the droplets’ flow is
larger than the Kolmogorov scale by order of magnitude. We provide numerical
data for the RDF (which was not done in [9]) and demonstrate that they confirm
the theory. We also derive in this section a sum rule. That provides the
probability density function (PDF) of the inter-pair distance irrespective of
the pair’s orientation in space, which equals the angle-averaged RDF. We
demonstrate that small angles, that correspond to preferential vertical
orientation of the pairs, can be neglected in the PDF of the distances. This
imposes a constraint on the magnitude of preferential orientation.

The complete angle dependence of the RDF of equal size particles is derived in
section VII. Section VIII extends the calculations to different size particles
by providing the bidisperse RDF. Section IX describes reduction of collision
kernel of different size particles to the RDF. We apply Yaglom-type relation
for rewriting the kernel and calculate the average velocity of approach of
colliding droplets. The difference from the classical work [27] is that
compressibility demands somewhat different approach in treating the kernel.
The approach velocity can be calculated completely due to Gaussianity of
gradients of the droplets’ flow. The next Section collects the information to
complete the calculation of the collision kernel. Section XI studies the
possibilities for interpolation of the results to the Reynolds numbers
characteristic of clouds.

The paper is quite long and for convenience of the reader we provide a rather
detailed summary and outlook in Section XII. Appendices provide technical
details of calculations. The main results of this work are: demonstration of
clash of $Fr\to 0$ and $Re_{\lambda}\to\infty$ limits; theory and numerical
confirmation of angle-dependent RDF at $Fr\to 0$; description of dimensions of
particles’ columns in space (see Fig. 1) and conjecture on the RDF in the high
Reynolds number turbulence in clouds.

##  II Equation of motion and its applicability in clouds

In this section, we consider the Newton equations of motion of the droplets.
The equations are characterized by three dimensionless parameters - the Froude
number $Fr=\epsilon_{0}^{3/4}/[g\nu^{1/4}]$, the Stokes number
$St=\tau/\tau_{K}$ and the Taylor microscale Reynolds number $Re_{\lambda}$,
see the Introduction. Our definition of $Fr$ is one of a number of possible
definitions. It is chosen here because of its independence of the parameters
of the droplets: $Fr$ characterizes turbulence with respect to gravity and is
a property of the flow only, cf. [14] and below. All the three parameters play
significant role in defining the nature of the particle trajectories. We
describe the parameter ranges that are relevant for the rain formation
problem. We stress the possibility that the usually used equation of motion
could actually not apply in the clouds because of the strong intermittency.

Typical values—We consider spherical droplets with radii $a$ from $10$ up to
$60$ microns, which is the size range where turbulence is most relevant in the
formation of larger droplets [1]. For a given droplet radius $a$, the Stokes
and Froude numbers are not independent. We consider the Stokes number
$St_{s}(a)\equiv\tau_{s}(a)\sqrt{\epsilon_{0}/\nu}$ defined with the help of
the Stokes time $\tau_{s}(a)=(2/9)(\rho/\rho_{a})(a^{2}/\nu)$ where $\rho_{a}$
is the density of air and $\rho$ is the droplet density (this is not $St$ in
the rest of the paper, see below). We have
$St_{s}(a)=0.02(a/15)^{2}\sqrt{\epsilon_{0}/10}$ where $a$ is measured in
microns, $\epsilon_{0}$ is measured in $cm^{2}/s^{3}$ and we use the numerical
values of $\nu=0.15\ cm^{2}/s$ and water-to-air density ratio of $814$. This
formula can be used for studying the size dependence of $St_{s}$ at different
levels of turbulence, as characterized by $\epsilon_{0}$. For instance, using
$a$ of $10$ microns we find $St_{s}\gtrsim 0.01$ for $\epsilon_{0}\gtrsim 10$.
We conclude from considerations of subsection I.1 that in the clouds, in the
whole range of relevant droplet sizes, the RDF and the collision kernel are
determined by vortices with time-scale of order $\tau$.

We also observe that $St_{s}(a)=0.5(a/15)^{2}Fr^{2/3}$. We find that at
$Fr=0.03$ (corresponding to $\epsilon_{0}=48\ cm^{2}/s^{3}$ by
$\epsilon_{0}/10=517Fr^{4/3}$), the range of $St_{s}(a)\gtrsim 1$ corresponds
to droplets’ radii larger than forty microns. When $Fr=0.05$ (which is
$\epsilon_{0}$ about $100cm^{2}/s^{3}$), the range of sizes with
$St_{s}(a)\gtrsim 1$ is somewhat larger than thirty microns. For not weak
turbulence with $Fr=0.29$ (with $\epsilon\approx 1000\ cm^{2}/s^{3}$), we have
$St_{s}(a)\gtrsim 1$ at size of 15 - 20 microns. Thus allows to see the range
of sizes with non-small inertia for different strengths of turbulence.

Equation—We use the effective linear drag description within which the
droplet’s coordinate $\bm{x}(t)$ and velocity $\bm{v}(t)$ obey

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\frac{d\bm{x}}{dt}=\bm{v},\
\ \frac{d\bm{v}}{dt}=-\frac{\bm{v}-\bm{u}(t,\bm{x}(t))}{\tau}+\bm{g},$ (1)

where $\bm{u}(t,\bm{x})$ is the incompressible homogeneous turbulent flow
(obeying the Navier-Stokes equation), $\bm{g}$ is the vector of gravitational
acceleration, and $\tau$ is the effective relaxation time, described below. We
assume that the droplets form a dilute gas; therefore, the droplet’s motion
can be considered independently of other droplets. The equation’s validity
demands that the particles are much denser than the fluid, which is true for
liquid droplets in air. The condition of validity of the linear friction force
is that the flow changes weakly at the scale of the particle, i.e. $a$ is much
smaller than the local viscous scale [8], and the Reynolds number of the flow
perturbation due to the particle $Re_{p}$, given by the product of the drift
velocity $|\bm{v}-\bm{u}(t,\bm{x}(t))|$ and $2a/\nu$, is small,
$Re_{p}=2|\bm{v}-\bm{u}(t,\bm{x}(t))|a/\nu\ll 1$. We consider these conditions
with the account of intermittency.

Role of fluctuations of viscous scale—Due to intermittency, the local viscous
scale has strong fluctuations where with appreciable probability its values
may be smaller than the Kolmogorov scale
$\eta\equiv\nu^{3/4}/\epsilon_{0}^{1/4}$ by a power of $Re_{\lambda}$, see
e.g. [8] for theoretical references and [13] for a recent numerical
observation. In clouds, the Kolmogorov scale is one or two orders of
magnitudes larger than relevant droplets’ sizes $a$, see numerical estimates
below. Since the viscous scale associated with extreme events is smaller than
the Kolmogorov scale by $Re_{\lambda}^{\delta}$, where $\delta$ is about $1/2$
according to [13], then it seems quite certain that Eqs. (1) do not apply to
extreme events in clouds. However, this is not necessarily of concern since
these events might be irrelevant for the space-averaged rate of droplets’
collisions. What is more relevant is whether the viscous scale associated with
those vortices that are relevant is larger than $a$. The present work
indicates that the kernel is determined by events for which the local Froude
number $\epsilon^{3/4}/[g\nu^{1/4}]$ is of order one. For these events
velocity gradients are larger than the typical value of
$\sqrt{\epsilon_{0}/\nu}$ by a factor of $Fr^{-2/3}$ (with $Fr$ always
standing for $\epsilon_{0}^{3/4}/[g\nu^{1/4}]$). The corresponding local
viscous scale is smaller than the Kolmogorov scale by factor of $Fr^{-1/3}$
where we assume that the local viscous scale is of order of
$\nu^{3/4}\epsilon^{-1/4}$. Since in practice $Fr^{1/3}\sim 1$ then the
assumption that the local viscous scale of cloud vortices relevant for
collision kernel is much larger than $a$ appears self-consistent.

Some reservations need to be made. We used in the above estimates the standard
phenomenology of turbulence [8] which could be invalid [13]. This might demand
a reconsideration. Moreover self-consistency of the assumption that the local
viscous scale of relevant vortices is much larger than $a$ may not yet
guarantee that vortices with scale smaller than $a$ can be neglected in the
collision kernel. For instance, self-consistency argument fails for the sling
effect: the calculation of the collision kernel that neglects the sling effect
is self-consistent in weakly intermittent turbulence (see the Introduction),
yet the sling effect’s contribution may be appreciable [4]. Therefore it might
be necessary to perform a separate study of the effect of vortices with
spatial scale $\lesssim a$. These vortices could cause a kind of ”sling
effect” of their own because they could generate large velocity difference of
nearby particles. This would demand separate account of these vortices in the
collision kernel. It seems that the only way to study this possibility is by
performing numerical simulations at high $Re_{\lambda}$. These simulations
must use the description of motion that is more fundamental than Eq. (1) and
applies also to vortices with size $\lesssim a$. This is left for future work.

Remaining condition of $Re_{p}\ll 1$—For the study of this condition it is
useful to integrate the velocity equation

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\bm{v}(t)=\bm{v}(0)\exp\left(-\frac{t}{\tau}\right)+\int_{0}^{t}dt^{\prime}\exp\left[-\frac{t-t^{\prime}}{\tau}\right]\left(\frac{\bm{u}[t^{\prime},\bm{x}(t^{\prime})]+\bm{g}\tau}{\tau}\right).$
(2)

We find that, after the particle spent in the flow time larger than $\tau$, we
have

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\bm{v}-\bm{u}(t,\bm{x}(t))=\bm{g}\tau+\int_{-\infty}^{t}dt^{\prime}\exp\left[-\frac{t-t^{\prime}}{\tau}\right]\left(\frac{\bm{u}[t^{\prime},\bm{x}(t^{\prime})]-\bm{u}[t,\bm{x}(t)]}{\tau}\right),$
(3)

where we rearranged the terms so that the RHS provides the velocity of the
drift with respect to the local flow. The drift has a contribution due to the
sedimentation velocity $\bm{g}\tau$ and the inertial lag behind the flow. For
the study of the RHS we start with instructive case of negligible gravity.

Figure 1: DNS snapshots of particles’ distributions in space at fixed $St=1$,
$Re_{\lambda}=70$ and increasing $Fr$ of (a) $Fr=0.0125$, (b) $Fr=0.0167$, (c)
$Fr=0.0333$, and (d) $Fr=0.05$. The distributions are multifractal. The
positive difference of $3$ and the fractal dimension increases linearly with
$Fr$ in the considered range. The fractal structure is anisotropic. Pronounced
columns of particles predominantly point in the direction of gravity, that was
slightly slanted from the $y$-axis to avoid numerical artifact due to
periodicity. We demonstrate in the text that characteristic dimensions of the
columns are $g\tau^{2}$ in the direction of gravity and $l_{c}\simeq 10\eta$
in the perpendicular direction. This is confirmed by the Figure. Indeed $St=1$
and Kolmogorov scale $\eta=0.024$, imply $g\tau^{2}=\eta/Fr$. Thus the theory
predicts that characteristic longest dimension of columns is $0.024/Fr$ and
the characteristic shortest dimension is $0.24$. This agrees with the
snapshots where columns become more convoluted at increasing $Fr$. There is no
apparent correlation with the vorticity component in the perpendicular
direction, which would hold for the centrifugal mechanism of clustering.

###  II.1 Drift velocity at negligible gravity

Here we consider the self-consistency of the assumption $Re_{p}\ll 1$ in the
case of negligible gravity, $Fr\to\infty$.

Case of $St\ll 1$—Provided that the local Stokes number
$\tau\sqrt{\epsilon/\nu}$ is much smaller than one, we have the Maxey formula
[12]

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\bm{v}-\bm{u}(t,\bm{x}(t))=-\tau\bm{a}(t,\bm{x}(t))+O\left(\frac{\tau^{2}}{\tau_{\nu}^{2}}\right),\
\ \bm{a}\equiv\partial_{t}\bm{u}+(\bm{u}\cdot\nabla)\bm{u},\ \
\frac{1}{\tau_{\nu}}\equiv\left|\frac{1}{a}\frac{da}{dt}\right|,$ (4)

where we introduced the field of Lagrangian accelerations of the fluid
particles $\bm{a}(t,\bm{x})$ and local time of variations of viscous scale
eddies $\tau_{\nu}$. The above formula is readily confirmed by using Taylor
expansion of velocity difference in the integrand in Eq. (3) at
$t^{\prime}=t$. The equation breaks down for fast vortices with
$\tau_{\nu}\lesssim\tau$. However we demonstrated in the Introduction that, in
the limit of large $Re_{\lambda}$, both the RDF and the collision kernel are
determined by vortices with $\tau_{\nu}\sim\tau$. We conclude that for clouds
Eq. (4), despite that it is true in most of the space, see the Introduction,
is irrelevant even for $St\ll 1$ (the equation could still be used for
calculating the collision kernel of droplets with extremely small $St$,
however this would not have practical relevance).

We see from the above that we cannot estimate the drift velocity with the help
of Eq. (4). The velocity is estimated by observing that the relevant vortices
with $\tau_{\nu}\sim\tau$ cause the particles to move with respect to the flow
at velocity of order $\sqrt{\nu/\tau}$. This is the typical velocity of the
(local) viscous scale eddies with turnover time $\tau$ (this assumes the usual
phenomenology of turbulence [8] which is not obviously true at any
$Re_{\lambda}$, cf. above and [13]). We conclude that the Reynolds number
$Re_{p}=2|\bm{v}-\bm{u}(t,\bm{x}(t))|a/\nu$ of relevant flow perturbations due
to the particle is of order $\sqrt{a^{2}/\nu\tau}$. Using the Stokes formula
$\tau_{s}(a)=2a^{2}\rho/(9\nu\rho_{a})$ and $\tau\sim\tau_{s}(a)$ we find that
the condition $Re_{p}\ll 1$ gives $\sqrt{\rho_{a}/\rho}\ll 1$. Here $\rho_{a}$
is the density of air and $\rho$ is the droplet density. This condition is
independent of the particle radius $a$ and it is obeyed by liquid droplets in
air. We remark that the obtained condition $\sqrt{a^{2}/\nu\tau}\ll 1$
coincides with the condition that $a$ is much smaller than the local viscous
scale $\sqrt{\nu\tau}$ of vortices with time-scale $\tau$.

We conclude that Eq. (1) is valid at $St\ll 1$ for all relevant vortices.
Moreover we can asymptotically continue this conclusion to $St\sim 1$ since at
$St\sim 1$ the statistics is still determined by vortices with
$\tau_{\nu}\sim\tau$. Thus at zero gravity and $St\lesssim 1$ the usage of the
equation of motion seems valid, up to reservations described above.

Case of $St\gg 1$—It remains that we consider the case of $St\gg 1$ and
negligible gravity where the time-scale $\tau$ belongs in the inertial range
of turbulence (there could also be the case of $\tau$ comparable or larger
than the eddy turnover time of the integral scale turbulence. However this
case does not seem to have applications in the rain formation problem and will
not be considered). If $Re_{\lambda}$ is moderate so that intermittency is
negligible, then we can use the dimensional Kolmogorov-type estimate
$\sqrt{\epsilon_{0}\tau}$ for the drift velocity [24]. This gives
$Re_{p}=\sqrt{\epsilon_{0}\tau}(a/\nu)=\sqrt{St}(a/\eta)$ which might not be
small for large $St$ particles. At higher $Re_{\lambda}$, where intermittency
is relevant, we can use Landau-type argument [8]. Within it, the velocity is
estimated by changing $\epsilon_{0}$ with the local energy dissipation rate
$\epsilon$ i.e. is given by $\sqrt{\epsilon\tau}$. The resulting changes
depend on which flow fluctuations are relevant in the collision kernel at
$St\gg 1$. Theoretical study is yet to be done.

Our conclusion is that $Re_{p}\ll 1$ is self-consistent for particles with not
too large $St$. However for certain $St>1$, whose value depends on the
intermittency, we have $Re_{p}\gtrsim 1$ for relevant fluctuations of
turbulence, the flow perturbation due to the particle is non-linear and Eqs.
(1) break down.

###  II.2 Reynolds number of flow perturbation $Re_{p}$ at $Fr\ll 1$ and
definition of $St$

We consider the full Eqs. (3) in the case of non-negligible gravity. We assume
$Fr\ll 1$ aiming at asymptotic description of the clouds, as explained in the
Introduction. In this case the drift velocity is given in the leading order by
the sedimentation velocity in the still air,
$\bm{v}-\bm{u}[t,\bm{x}(t)]\approx\bm{g}\tau$; see [9]. This can be readily
demonstrated at $St\ll 1$ where the particles’ acceleration $\dot{\bm{v}}$ is
of the order of the acceleration of the fluid particles. This, by definition
of $Fr\ll 1$, is much smaller than $\bm{g}$, so that the LHS of the velocity
equation in Eq. (1) is negligible, giving
$\bm{v}-\bm{u}[t,\bm{x}(t)]\approx\bm{g}\tau$. At $St\gtrsim 1$ we use that
Eq. (3) implies
$\left|\bm{v}-\bm{u}(t,\bm{x}(t))-\bm{g}\tau\right|\sim\left|\bm{u}(t-\tau,\bm{x}(t-\tau))-\bm{u}(t,\bm{x}(t))\right|$
where we estimate the integral as the velocity difference in the particle
frame at time lag $\tau$. This difference is either due to turbulence’s
changes in time or in space. The former are given by the typical velocity
$\sqrt{\epsilon_{0}\tau}$ of eddies with time-scale $\tau$, see above. In
turn, the variations in space are caused by the particle crossing during time
$\tau$ of the distance $g\tau^{2}$ due to sedimentation. At moderate
$Re_{\lambda}$ with negligible intermittency, this contributes to
$\left|\bm{u}(t-\tau,\bm{x}(t-\tau))-\bm{u}(t,\bm{x}(t))\right|$ the
Kolmogorov velocity difference $(\epsilon_{0}g\tau^{2})^{1/3}$ between spatial
points separated by $g\tau^{2}$. We find the estimate

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\left|\bm{v}-\bm{u}(t,\bm{x}(t))-\bm{g}\tau\right|\sim
max\left[\sqrt{\epsilon_{0}\tau},(\epsilon_{0}g\tau^{2})^{1/3}\right]\sim\sqrt{\epsilon_{0}\tau}max\left[1,St^{1/6}Fr^{-1/3}\right].$
(5)

We conclude from this formula by using $Fr\ll 1$ and $St\gtrsim 1$ that the
drift velocity is $\bm{g}\tau$ by using

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\frac{\left|\bm{v}-\bm{u}(t,\bm{x}(t))-\bm{g}\tau\right|}{g\tau}\sim
max\left[FrSt^{-1/2},Fr^{2/3}St^{-1/3}\right]\ll 1,$ (6)

which proves $\bm{v}-\bm{u}(t,\bm{x}(t))\approx\bm{g}\tau$. At higher Reynolds
numbers where intermittency is relevant, significant changes might be
necessary. Their implementation demands the currently missing knowledge of
which type of intermittent fluctuations, whose values by themselves cover a
wide range of orders of magnitude, are relevant.

We find that the Reynolds number $Re_{p}$ of the flow perturbation by the
particle, $Re_{p}=2|\bm{v}-\bm{u}(t,\bm{x}(t))|a/\nu$ equals $2ag\tau/\nu$.
This number is small at sizes smaller than thirty microns where $\tau$ in Eq.
(1) is the Stokes time $\tau_{s}(a)=(2/9)(\rho/\rho_{a})(a^{2}/\nu)$, given by
the droplet’s mass $4\pi\rho a^{3}/3$ divided by the coefficient of the Stokes
force $6\pi\rho_{a}\nu a$. In contrast, at radii from $30$ to $60$ microns we
have $Re_{p}\gtrsim 1$. There an effective function $\tau(a)$, different from
$\tau_{s}(a)$ must be used as $\tau$ in Eqs. (1); see [55, 1]. It is with the
help of this function that we define the Stokes number
$St(a)=\tau(a)\sqrt{\epsilon_{0}/\nu}$, which determines the strength of the
droplets’ inertia, and not $\tau_{s}(a)$, cf. [55, 14, 9]. The times $\tau(a)$
and $\tau_{s}(a)$ are of the same order in the size range of interest; the use
of $\tau(a)$ in equation of motion, as compared with the use of $\tau_{s}(a)$,
was found to decrease the collision kernel of larger droplets by up to $26$
per cent [55].

###  II.3 Range of considered parameters

In the rest of the paper, unless told otherwise, we assume that $St\gtrsim 1$,
since there is already a well-developed theory for droplets whose size obeys
$St(a)\ll 1$. The traditional small $St$ theory relies on Eq. (4), see [12,
4]. However in some cases Eq. (4) does not hold at $St\ll 1$ because of the
gravity; see the study of all possibilities in [9]. We remark that other sets
of dimensionless parameters, different from our $St$, $Fr$ and $Re_{\lambda}$
are in use, cf. e.g. the studies [55, 1] who use the parameter
$Sv=g\tau(a)/(\epsilon\nu)^{1/4}=St/Fr$, which is the ratio of the
gravitational settling velocity of the particle and the velocity of turbulent
eddies at the Kolmogorov scale. This parameter mixes characteristics of
particles, gravity, and turbulence and is much larger than one in the range of
$Fr\ll 1$ and $St\gtrsim 1$ that we study.

In clouds the velocity at the integral scale of turbulence $L$ is much larger
than the sedimentation velocity. Thus we assume that the first term in
$\bm{v}(t)\approx\bm{u}(t,\bm{x}(t))+\bm{g}\tau$ dominates the sum, that is
the typical value of the integral scale velocity $(\epsilon_{0}L)^{1/3}$ is
much larger than $g\tau$ for relevant $\tau$ (for the limit where $\bm{v}(t)$
is dominated by sedimentation and $g\tau\gg(\epsilon_{0}L)^{1/3}$ see [12] and
also [22, 21]). We must keep $\bm{g}\tau$ in $\bm{v}(t)$, despite that it is
much smaller than $\bm{u}(t,\bm{x}(t))$ since it dominates at small scales
where the clustering occurs, exceeding the Kolmogorov scale velocity. The
approximation $\bm{v}(t)\approx\bm{u}(t,\bm{x}(t))+\bm{g}\tau$ is equivalent
to stating that the droplets are transported by the flow
$\bm{u}(t,\bm{x})+\bm{g}\tau$. This flow is incompressible and would bring no
clustering. Inclusion of further corrections to this approximation is
necessary for the description of inhomogeneous spatial distribution of the
droplets in the steady state.

###  II.4 Direct numerical simulations

In order to test below the theoretical predictions we performed DNS of the
particle-laden isotropic turbulence. The Navier-Stokes equation was
numerically solved on $128^{3}$ grids using a spectral method on the periodic
cubic domain to describe the homogeneous isotropic turbulence at
$Re_{\lambda}=70$. The equation of particle motion (Eq. 1) was solved by
taking into account the linear Stokes drag and gravity together. The initial
positions and velocities of the particles are random and local fluid
velocities, respectively. Information of fluid quantities at the particle
position was obtained using the fourth-order Hermite interpolation scheme [65,
66]. The details of the numerics can be found in Refs. [67, 68, 69, 70].
Typical snapshots of the particles’ distribution in the steady state are shown
in Fig. 1. The angle dependence of the RDFs was computed in a statistically
steady condition with five different populations of $N_{p}=9,971,200$
particles for each Froude number, $Fr=0.0125,0.0167,0.025,0.033,0.05$. The
Stokes number is fixed at $St=1$. Here, $N_{p}$ was determined by
$N_{\eta}=2.33$, which is the average number of particles within a sphere with
radius $\eta$. For the estimation of the correlation codimension $\alpha$
using Eq. (44), the Lyapunov exponents were computed by releasing many pairs
of particles. The initial distance between the particles is set to
$10^{-4}\eta$, and the change in distance between two particles, the area
between three particles, and the volume constructed from four particles were
measured on the basis of Gram-Schmidt renormalization for a period of
$80\tau_{\eta}$ after the transient period due to the arbitrary initial
condition. There are 10,000 sets of pairs released in one flow field, and data
are collected over a total of 12 flow fields. The same method was used in our
earlier works on droplets and bubbles [9, 71]. In our simulations the only
parameter which we changed was the gravity acceleration that changed $Fr$.

##  III Lyapunov exponents, columns and their disappearance as $Re_{\lambda}$
grows

In this section we study exponential separation of two particles in the
viscous range. This is probably the shortest path to deriving the small-scale
columnar structure formed by the particles. We provide a concise formula for
the Lyapunov exponent $\lambda_{1}$ of the droplets via the energy spectrum of
turbulence $E(k)$. These results were obtained in [9] where also further
references are provided. Here we use a simpler approach, similar to [24], cf.
[22]. This allows us to address the breakdown of the theory at increasing
Reynolds number. The breakdown implies that the columns that were observed in
the direct numerical simulation might disappear if $Re_{\lambda}$ is increased
at fixed $Fr$ and $St$. We also provide the Lyapunov exponent in a range of
non-small $Fr$, that was not derived previously.

It is a direct consequence of the equation of motion (1), that for $r$ deep
inside the viscous range, the separation vector of the equal size droplets
$\bm{r}$ obeys

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\dot{\bm{r}}=\bm{V},\ \
\dot{\bm{V}}=-\frac{\bm{V}-s\bm{r}}{\tau};\ \
s_{ik}(t)=\nabla_{k}u_{i}(t,\bm{x}(t)),$ (7)

where $s_{ik}$ is the matrix of gradients of turbulent flow in the frame of
one of the particles whose trajectory is designated by $\bm{x}(t)$. We will
use this equation for the calculation of $\lambda_{1}$ which will be seen to
be determined by vortices whose size $L_{c}$ is in the inertial range and much
larger than the Kolmogorov scale $\eta$. Therefore the result can be applied
to the evolution of $r(t)$ as long as $r(t)\ll L_{c}$ and not under the more
stringent $r\ll\eta$.

We first consider Eq. (7) at moderate Reynolds numbers with weak intermittency
where rigorous study of $Fr\to 0$ limit is possible. The impact of
intermittency on our considerations will be considered later (it must be
observed however that $L_{c}\gg\eta$ is a manifestation of intermittency. Our
detailed assumption will be seen below).

###  III.1 Derivation at moderate $Re_{\lambda}$

We study the parameters’ range of $Fr\ll 1$ and $St\gtrsim 1$ assuming that
$Re_{\lambda}$ is moderate so that several quantities in the study below obey
the Kolmogorov-type estimates [8]. The particle moves with respect to the
local flow at the velocity $\bm{g}\tau$, as explained above. This implies [4]
that the correlation time of the matrix $s_{ik}(t)$, providing velocity
gradients in the frame of the droplet, is the minimum of the Kolmogorov time
$\tau_{K}$ and the time $\tau_{g}\equiv\eta/(g\tau)$. Here $\tau_{g}$ is the
time during which the droplet moving through the flow at velocity $\bm{g}\tau$
crosses the correlation length of the turbulent velocity gradients, the
Kolmogorov scale $\eta$. Observing that $\tau_{g}/\tau_{K}=Fr/St$ is small in
the considered range of parameters ($Fr\ll 1$ and $St\gtrsim 1$) we conclude
that in our case $s_{ik}(t)$ varies in time at the time-scale $\tau_{g}$.
During this time-scale, which is much less than the turnover time of the
viscous scale eddies $\tau_{K}$, the field of the velocity gradients does not
change in time appreciably if considered in the frame of the fluid particles.
Therefore the temporal correlation function
$\left\langle\nabla_{k}u_{i}(t,\bm{x}(t))\nabla_{r}u_{p}(t^{\prime},\bm{x}(t^{\prime}))\right\rangle$
equals the spatial correlation function
$\left\langle\nabla_{k}u_{i}(t,\bm{x}(t))\nabla_{r}u_{p}(t,\bm{x}(t)+\bm{g}\tau(t^{\prime}-t))\right\rangle$.
Moreover the product of the typical value of $s_{ik}(t)$, which is
$\tau_{K}^{-1}$, and of the correlation time $\tau_{g}$ of $s(t)$ is small.
Considering this product as a small parameter, the leading order approximation
is the limit of zero correlation time where $s(t)$ in Eq. (7) can be replaced
by white noise [24].

We demonstrate in detail how the white noise description of the effect of $s$
on the evolution of $\bm{r}$ and $\bm{V}$ in Eq. (7) arises. We observe that
Eq. (7) gives

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\bm{V}(t)=\bm{V}(0)\exp\left(-\frac{t-t^{\prime}}{\tau}\right)+\frac{1}{\tau}\int_{0}^{t}\exp\left(-\frac{t-t^{\prime}}{\tau}\right)s(t^{\prime})\bm{r}(t^{\prime})dt^{\prime}.$
(8)

We concentrate on the phenomena where the particles stay inside the viscous
scale for times much larger than $\tau$. There the above equation holds at
$t\gg\tau$ and the first term can be neglected. There are physical phenomena
that occur on a time-scale $\tau$ or smaller so that for them $t\gg\tau$ is
not the range of interest. An example is the sling effect [4] where the
particles detach from the flow at the time of the beginning of the sling
effect, taken to be $t=0$, and then move ballistically toward each other so
that the first term in Eq. (8) dominates $\bm{V}(t)$. In this case, despite
that $r\ll\eta$ and the Taylor expansion of the velocity difference made in
Eq. (7) holds, what changes the distance between the particles is not the
difference of the turbulent flow velocities at the positions of the particles,
which is given by $s\bm{r}$, but their inertia. The sling effect lasts for
times of order $\tau$ and is not describable by the limit of $t\gg\tau$.

Phenomena that happen at times larger than $\tau$ include separation of two
infinitesimally close trajectories in the six-dimensional phase space. This
defines the first Lyapunov exponent $\lambda_{1}$ that describes the growth of
the separation vector $(\bm{r},\bm{v})$ between the trajectories via

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\lambda_{1}=\lim_{t\to\infty}\frac{1}{2t}\ln\left(\frac{r^{2}(t)+\tau^{2}V^{2}(t)}{r^{2}(0)+\tau^{2}V^{2}(0)}\right),$
(9)

where the dimensional factor $\tau$ is irrelevant and used only for having a
dimensionally uniform expression. Here, besides that $r(0)\ll\eta$, also the
initial velocity difference is small so that $r(t)$ stays much smaller than
$\eta$ during times much larger than $\tau$. The limit above exists and is
given by the same constant for almost all trajectories (i.e. with possible
exception of initial positions with zero volume in the phase space) [78].
Therefore disregarding the initial period of evolution of duration of order
$\tau$ (or simply setting $\bm{V}(0)=0$ using that the limit is independent of
the vector $(\bm{r}(0),\bm{V}(0))$ as long as this vector is non-zero) we can
use instead of Eq. (8) the simplified equation

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\bm{V}(t)=\frac{1}{\tau}\int_{-\infty}^{t}\exp\left(-\frac{t-t^{\prime}}{\tau}\right)s(t^{\prime})\bm{r}(t^{\prime})dt^{\prime}.$
(10)

We study the regime determined implicitly by the condition
$\lambda_{1}\tau_{g}\ll 1$, whose explicit form will be provided later by
obtaining $\lambda_{1}$. We introduce a separation time $\Delta t$ that obeys
$\tau_{g}\ll\Delta t\ll min[\lambda_{1}^{-1},\tau]$. This $\Delta t$ exists
because we assumed $\lambda_{1}\tau_{g}\ll 1$ and because
$\tau_{g}/\tau=Fr/St^{2}$ is small by $Fr\ll 1$ and $St\gtrsim 1$. Then we can
write

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\bm{V}(t)\approx\frac{1}{\tau}\int_{-\infty}^{t}\exp\left(-\frac{t-t^{\prime}}{\tau}\right)s_{\Delta
t}(t^{\prime})\bm{r}(t^{\prime})dt^{\prime},\ \ s_{\Delta
t}(t)\equiv\int_{t}^{t+\Delta t}s(t^{\prime})\frac{dt^{\prime}}{\Delta t},$
(11)

where $s_{\Delta t}$ is $s(t)$ coarse-grained over time-scale $\Delta t$. We
used that $\Delta t$ is smaller than both the characteristic time of
variations of $\bm{r}(t)$, which is $\lambda_{1}^{-1}$ and the characteristic
time of variations of the exponent $\tau$.

We observe that $s_{\Delta t}$ is proportional to the integral of the random
process $s(t)$ over a time-interval which is much larger than the correlation
time of this process, $\tau_{g}$. This implies by the central limit theorem
(CLT) that the statistics of $s_{\Delta t}$ is Gaussian. Indeed we can write
the integral in the definition of $s_{\Delta t}$ as sum of contributions of
many intervals whose duration is the correlaton time $\tau_{g}$, i.e.
$s_{\Delta t}(t)=\left(\Delta
t\right)^{-1}\sum_{i=1}^{N}\int_{t+(i-1)\tau_{g}}^{t+i\tau_{g}}s(t^{\prime})dt^{\prime}$
where $N=\Delta t/\tau_{g}$ is much larger than one (we can use in these
considerations $\Delta t$ that is an integer multiple of $\tau_{g}$). Then,
since the contributions of different intervals are independent, we find that
$s_{\Delta t}$ is sum of a large number of independent random variables and
the CLT applies. A more rigorous proof can be constructed by applying the
cumulant expansion theorem [57] to the characteristic function of $s_{\Delta
t}$. We conclude that the distribution of $s_{\Delta t}(t)$ is fully
determined by the mean and the dispersion that fix a Gaussian distribution
uniquely.

It can be checked that the mean of $s$ is negligible, see [9] and Appendix A
(we remark that the usual argument which uses that $\langle s_{ik}\rangle$ is
proportional to $\delta_{ik}$ due to isotropy of the small-scale turbulence
fails in this case. Gravity makes the statistics anisotropic). Therefore the
statistics of $s_{\Delta t}(t)$ coincides with the statistics of
$\int_{t}^{t+\Delta t}s^{K}(t^{\prime})dt^{\prime}/(\Delta t)$ where
$s^{K}(t)$ is a Gaussian matrix process with zero mean and zero correlation
time which dispersion is picked so that dispersions of $\int_{t}^{t+\Delta
t}s^{K}(t^{\prime})dt^{\prime}/(\Delta t)$ and $s_{\Delta t}(t)$ agree. The
white noise matrix process $s^{K}(t)$ arises in the Kraichnan model of the
turbulent transport as indicated by the superscipt [36]. It is readily seen
that our condition on the dispersion of $s^{K}(t)$ demands that
$\int_{-\infty}^{\infty}\left\langle s_{ik}(t)s_{pr}(t^{\prime})\right\rangle
dt^{\prime}$ equals $\int_{-\infty}^{\infty}\left\langle
s^{K}_{ik}(t)s^{K}_{pr}(t^{\prime})\right\rangle dt^{\prime}$. Finally using
our previous consideration of correlations of $s_{ik}(t)$ we find

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\left\langle
s^{K}_{ik}(t)s^{K}_{pr}(t^{\prime})\right\rangle=\delta(t-t^{\prime})\kappa_{ikpr},\
\
\kappa_{ikpr}\equiv\int_{-\infty}^{\infty}\left\langle\nabla_{k}u_{i}(0)\nabla_{r}u_{p}(\bm{g}\tau
t)\right\rangle dt,$ (12)

where the correlation function in the integrand is the equal-time correlation
function of turbulent velocity gradients. Calculation of $\kappa_{ikpr}$ made
in [9] reveals that in the leading order $s^{K}_{ik}$ is a $2\times 2$ random
matrix since all the entries of the matrix that contain index $z$ vanish. This
$2\times 2$ matrix obeys the usual statistics of two-dimensional Kraichnan
model determined by restriction of $\kappa_{ikpr}$ to indices different from
$z$

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\kappa_{\alpha\beta\gamma\delta}=D\left(3\delta_{\alpha\gamma}\delta_{\beta\delta}-\delta_{\alpha\beta}\delta_{\gamma\delta}-\delta_{\alpha\delta}\delta_{\gamma\beta}\right),\
\ D\tau=\frac{\pi\int_{0}^{\infty}E(k)kdk}{8g}\sim Fr,$ (13)

where $E(k)$ is the energy spectrum of turbulence, cf. [24]. Here and below,
the Greek indices take values of $1$ or $2$ and summation over repeated
indices is implied (no confusion can be caused between $\alpha$ as index and
$\alpha$ as exponent). Two-dimensional Kraichan model
$\dot{\bm{v}}_{\perp}=-\left(\bm{v}_{\perp}-s^{K}\bm{r}_{\perp}\right)/\tau$,
where the subscript stands for horizontal components of the vectors, is well
studied [77, 30]. It is seen from the formula above that $D\tau\sim Fr\ll 1$.
This corresponds to the limit of small inertia where the approximation
$\dot{\bm{r}}_{\perp}\approx s^{K}\bm{r}_{\perp}$ holds. Thus the problem
reduces to solved evolution of distance between the two fluid particles in the
Kraichnan model [36]. This gives immediately the spectrum of the Lyapunov
exponents $\lambda_{i}$ defined by asymptotic growth rates of hypersurfaces
composed of particles. Thus $\lambda_{1}$, $\lambda_{1}+\lambda_{2}$, and
$\sum_{i=1}^{3}\lambda_{i}$ provide the asymptotic logarithmic growth rates of
the infinitesimal line, surface, and volume elements at large times,
respectively [36]. We have $\lambda_{1}=-\lambda_{3}$, that describes behavior
of $\bm{r}_{\perp}$ and $\lambda_{2}=0$ that describes conservation of the
vertical component of $\bm{r}$. The value of $\lambda_{1}$ is [36, 9]

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\lambda_{1}=D=\frac{\pi\int_{0}^{\infty}E(k)kdk}{8g\tau}.$
(14)

Self-consistency of the derivation demands that the obtained $\lambda_{1}$
obeys $\lambda_{1}\tau_{g}\ll 1$. This gives the condition $D\tau\ll 1$ which
is equivalent to $Fr^{2}/St^{2}\ll 1$. The inequality holds in the studied
range of parameters proving the self-consistency. The above prediction was
confirmed in [9] numerically by direct numerical simulations of the motion of
inertial particles in the Navier-Stokes turbulence with $Re_{\lambda}=70$ .
The result was found to be quantitatively accurate for $Fr\leq 0.05$, where
linear dependence on $Fr$ holds, and reasonably good for $Fr\leq 0.1$, see
Fig. 2. For higher $Fr$ the equation’s validity depends on $St$. The higher
the $St$ is, the more accurate the above formula is.

The general validity condition of the effective white-noise description above
is $Fr/St\ll 1$. This condition guarantees that correlation time of velocity
gradients in the particle’s frame is $\tau_{g}$. It also ensures that both
$\lambda_{1}\tau_{g}\ll 1$ and $\tau_{g}/\tau=Fr/St^{2}\ll 1$ hold, allowing
to pass from Eq. (10) to Eq. (11).

We observe that $Fr/St\ll 1$ can hold also for non-small $D\tau\propto Fr$ if
$St\gg 1$. Since the Froude number in clouds can be of order one then the case
of $St\gg Fr\sim 1$ might have practical applications. Therefore we provide
the results in this case for possible future use. We assume non-small Froude
number, $Fr\gtrsim 1$, and $St\gg Fr$, so that the white noise description
holds. If $D\tau\geq 1$ then $\lambda_{1}$ of the Kraichan model can be
written as
$\lambda_{1}=D^{1/3}\tau^{-2/3}{\tilde{\lambda}}_{1}\left[(D\tau)^{-1/3}\right]$.
Here ${\tilde{\lambda}}_{1}$ is a function that takes values of order one at
$D\tau\geq 1$ slowly changing with its argument from $0.5$ at $D\tau=1$ to $2$
at zero argument corresponding to $D\tau=\infty$, see [24] and references
therein. We find then that the self-consistency condition
$\lambda_{1}\tau_{g}\ll 1$ holds by $\lambda_{1}\tau_{g}\sim Fr^{4/3}/St^{2}$.
We conclude that

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\lambda_{1}=\frac{1}{\tau}\left(\frac{\pi\int_{0}^{\infty}E(k)kdk}{8g}\right)^{1/3}{\tilde{\lambda}}_{1}\left[\left(\frac{\pi\int_{0}^{\infty}E(k)kdk}{8g}\right)^{-1/3}\right],\
\ 1\lesssim Fr\ll St.$ (15)

Testing this prediction in the DNS is left for future work.

The above results imply that the droplets will form columns in space. Indeed
in the leading order the droplets separate horizontally keeping their vertical
separation conserved. Once the random trajectories of the particles bring one
of the particles on the top of the other, this pair configuration is preserved
in time indefinitely (more precisely for long time due to instability). The
particles form meta-stable bound states and after some time the space will be
divided into columns of particles that would coalesce and form a single
column. In reality the next order corrections cause gradual dissolution of the
vertical pair configuration so that columns are only the more probable
configuration of the particles. This phenomenon is quantitatively described in
later sections by studying the angular dependence of the radial distribution
function.

###  III.2 Implications of dissipation range statistics and breakdown at large
$Re_{\lambda}$

Here we assume that $Fr$ is fixed at a value much smaller than $1$ and study
the impact of increasing intermittency at growing $Re_{\lambda}$ on the
considerations above. The increase of $Re_{\lambda}$ can have two-fold effect
on $\lambda_{1}$. It can influence the estimate $\lambda_{1}\tau\sim Fr$,
implied by Eqs. (14)-(15), and it can also invalidate the assumptions made in
the derivation. We start from considering the former effect.

Comparison with Lyapunov exponent of tracers—It is useful to make comparison
with first Lyapunov exponent of tracer particles in the Navier-Stokes
turbulence $\lambda_{1}^{t}$. The dimensionless product
$\lambda_{1}^{t}\tau_{K}$, decays with $Re_{\lambda}$. This was predicted by
using the multifractal model in [44] and confirmed numerically in [45], see
also [40]. The reason for the decay is that increasing intermittency of
turbulence increases size of regions with quiescent quasilaminar turbulence
where the separating pair of particles consequently stays longer [8]. This
depletes chaos as measured by the dimensionless Lyapunov exponent [44]. We
consider if there is a similar dependence of $\lambda_{1}$ in Eq. (14) on
$Re_{\lambda}$.

Influence of non-trivial structure of the dissipation range—The dimensionless
Lyapunov exponent $\lambda_{1}\tau$ obeys

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\lambda_{1}\tau=\frac{\pi\int_{0}^{\infty}E(k)kdk}{8g}=c_{0}Fr;\
\ \
c_{0}\equiv\frac{\pi\nu^{1/4}\int_{0}^{\infty}\\!\\!E(k)kdk}{8\epsilon_{0}^{3/4}},$
(16)

where we defined dimensionless constant $c_{0}$. The integral in $c_{0}$ is
determined by the form of the turbulence energy spectrum $E(k)$ in the
dissipative range. Therefore it can be assumed to be independent of the large-
scale forcing so that $c_{0}$ is a dimensionless function of $Re_{\lambda}$.
It can be seen from the numerical data of [9] and Fig. 2 that $c_{0}\approx
1.5$ at $Re_{\lambda}=70$. The existing data on the spectrum can be used to
study the dependence of $c_{0}$ on $Re_{\lambda}$. It was found in [83, 53]
that the dissipation range spectrum observed in their DNS was well described
with

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!E(k)=C(k\eta)^{{\tilde{\alpha}}}\exp\left[-\beta
k\eta\right],$ (17)

where $C$, ${\tilde{\alpha}}$, and $\beta$ are functions of the Reynolds
number $Re_{\lambda}$. We find assuming that the contribution of smaller
wavenumbers can neglected and using the equation above that

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!c_{0}=\frac{\pi\int_{0}^{\infty}\\!\\!E(k)kdk}{8\eta\int_{0}^{\infty}\\!\\!E(k)k^{2}dk}=\frac{\pi\int_{0}^{\infty}\\!\\!k^{1+{\tilde{\alpha}}}\exp\left[-\beta
k\eta\right]dk}{8\eta\int_{0}^{\infty}\\!\\!k^{2+{\tilde{\alpha}}}\exp\left[-\beta
k\eta\right]dk}=\frac{\pi\beta}{8(2+{\tilde{\alpha}})},$ (18)

where we used $\int_{0}^{\infty}\\!\\!E(k)k^{2}dk=\epsilon_{0}/\nu$ and
$\eta=\nu^{3/4}/\epsilon_{0}^{1/4}$. The calculation is self-consistent
provided that ${\tilde{\alpha}}>-2$ since otherwise
$\int_{0}^{\infty}\\!\\!k^{1+{\tilde{\alpha}}}\exp\left[-\beta k\eta\right]dk$
diverges at small $k$ where Eq. (17) does not apply. Moreover
$2+{\tilde{\alpha}}$ must be not too small since otherwise the contribution of
small wavenumbers would still be appreciable. The measurement of $\alpha$
provided in [83] and described by the fit $\alpha=7.3Re_{\lambda}^{-0.47}-2.9$
gives that the condition ${\tilde{\alpha}}>-2$ is obeyed by $Re_{\lambda}$
smaller than about $85$. For $Re_{\lambda}=70$ we have $2+\alpha=0.09$ which
is too small for the calculation to be self-consistent so the fits for
$\alpha$ and $\beta$ provided in [83, 53] cannot be used for evaluating
$c_{0}$. However these fits give unequivocal indication that $c_{0}$ has
appreciable growth with $Re_{\lambda}$. For further details on the spectrum
see [53] for the DNS aspects and [64] for the theory.

Correlation length of gradients of particles’ flow is much larger than the
Kolmogorov scale— Further insight into $\lambda_{1}$ is reached by rewriting
the above prediction for the Lyapunov exponent as

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\lambda_{1}=D=\frac{\kappa_{\alpha\beta\alpha\beta}}{8}=\frac{1}{8}\int_{-\infty}^{\infty}\left\langle\nabla_{\beta}u_{\alpha}(0)\nabla_{\beta}u_{\alpha}(\bm{g}\tau
t)\right\rangle
dt=\frac{1}{4g\tau}\int_{0}^{\infty}\left\langle\nabla_{\beta}u_{\alpha}(0)\nabla_{\beta}u_{\alpha}(z{\hat{z}})\right\rangle
dz,$ (19)

where ${\hat{z}}$ is directed upwards. The RHS is similar to the dispersion of
the finite-time Lyapunov exponent studied in [44], which behaves as time
integral of the different time Lagrangian correlation function of velocity
gradients
$\int_{-\infty}^{\infty}\left\langle\nabla_{k}u_{i}(0)\nabla_{k}u_{i}(t)\right\rangle
dt$. The last integral can be written as $2t_{c}/\tau_{K}^{2}$ where
$t_{c}\equiv\int_{0}^{\infty}\left\langle\nabla_{k}u_{i}(0)\nabla_{k}u_{i}(t)\right\rangle
dt/\left\langle(\nabla\bm{u})^{2}\right\rangle$ is the effective correlation
time $t_{c}$. It was found that $t_{c}/\tau_{K}$ obeys rather strong increase
with the Reynolds number given by $Re_{\lambda}^{\kappa}$ with $\kappa\sim
0.1$ which qualitatively agrees with the numerical observations in [45]. The
reason is that at increasing $Re_{\lambda}$, due to intermittency, the regions
of moderate velocity gradients of order $\tau_{K}^{-1}$ become larger both in
space and in time [8]. This results in power-law increase of $t_{c}$ with
$Re_{\lambda}$. This makes it highly probable that also the spatial
correlation length $L_{c}$ defined by

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!L_{c}\equiv\frac{\int_{0}^{\infty}\left\langle\nabla_{\beta}u_{\alpha}(0)\nabla_{\beta}u_{\alpha}(z{\hat{z}})\right\rangle
dz}{\left\langle\nabla_{\beta}u_{\alpha}\nabla_{\beta}u_{\alpha}\right\rangle}=\frac{5\nu\int_{0}^{\infty}\left\langle\nabla_{\beta}u_{\alpha}(0)\nabla_{\beta}u_{\alpha}(z{\hat{z}})\right\rangle
dz}{2\epsilon_{0}},$ (20)

increases with $Re_{\lambda}$ so that $L_{c}/\eta\propto
Re_{\lambda}^{\Delta_{\lambda}}$ with $\Delta_{\lambda}\sim 0.1$ (the power-
law dependence is probably valid quantitatively at large $Re_{\lambda}$ and at
moderate $Re_{\lambda}$ is valid qualitatively only). Here we used that small-
scale isotropy, incompressibility and spatial homogeneity imply that single-
point statistics of turbulent velocity gradients obeys

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\left\langle\nabla_{k}u_{i}\nabla_{r}u_{p}\right\rangle=\frac{\epsilon_{0}}{30\nu}\left(4\delta_{ip}\delta_{kr}-\delta_{ik}\delta_{pr}-\delta_{ir}\delta_{pk}\right),\
\
\left\langle\nabla_{\beta}u_{\alpha}\nabla_{\beta}u_{\alpha}\right\rangle=\frac{2\epsilon_{0}}{5\nu}.$
(21)

(This formula is found by differentiation of the velocity pair correlation
function in [82].) We find from Eq. (19) that

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\lambda_{1}\tau=\frac{\epsilon_{0}L_{c}}{10g\tau\nu}=Fr\frac{L_{c}}{10\eta},$
(22)

which would give $\lambda_{1}\tau\sim FrRe_{\lambda}^{\Delta_{\lambda}}$ on
assuming the power-law dependence of $L_{c}/\eta$ on $Re_{\lambda}$. Using
$\lambda_{1}\tau=1.5Fr$, which was observed in the simulations of [9] at
$Fr\leq 0.033$, see Fig. 2, we find that $L_{c}=15\eta$ at $Re_{\lambda}=70$.
The large numerical factor demonstrates failure of dimensional estimates. The
reason for this factor is the non-trivial structure of the energy spectrum in
the dissipation range. Indeed, $L_{c}$ is proportional to $c_{0}$ considered
previously.

Theory breakdown at large $Re_{\lambda}$—The observation that the scale of
relevant flow configurations $L_{c}$ is much larger than $\eta$ and grows with
$Re_{\lambda}$ as $L_{c}/\eta\propto Re_{\lambda}^{\Delta_{\lambda}}$ implies
breakdown of the assumptions of the theory at large $Re_{\lambda}$. The
derivation of $\lambda_{1}$ assumed that the correlation time of velocity
gradients in the frame of the particle is the time during which the droplet
crosses the spatial correlation length $\eta$. However the correlation length
of relevant gradients is $L_{c}$ and not $\eta$, which is much larger than
$\eta$ already at moderate $Re_{\lambda}=70$. Using in the previous
considerations $L_{c}/(g\tau)$ as the correlation time of $s(t)$, we find that
the condition that $s_{\Delta t}$ is a sum of large number of independent
random variables demands that $L_{c}/(g\tau)\ll\tau$ and
$L_{c}/(g\tau)\ll\lambda_{1}^{-1}$. We observe that the former condition is
stronger because $\lambda_{1}\tau$ is never large, $\lambda_{1}\tau\lesssim
1$, at Stokes numbers of order one considered here. Indeed, inertia causes the
Lyapunov exponent of the particles, that tend to move ballistically, to be
smaller than the Lyapunov exponent of tracer particles $\lambda_{1}^{t}$.
Therefore $\lambda_{1}^{t}\tau\sim St$ implies $\lambda_{1}\tau\lesssim 1$. We
conclude that the condition of the theory applicability boils down to
$L_{c}/(g\tau)\ll\tau$. We find assuming the power-law dependence of
$L_{c}/\eta$ on $Re$ that our derivation of $\lambda_{1}$ holds provided that

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!L_{c}\ll g\tau^{2},\ \ \
\ \
15\left(\frac{Re_{\lambda}}{70}\right)^{\Delta_{\lambda}}\ll\frac{St^{2}}{Fr}.$
(23)

We observe that there is large numerical factor in Eq. (23). This limits the
theory’s applicability to rather small $Fr$: at $Re_{\lambda}=70$ deviations
from the prediction for $\lambda_{1}$ are observed at $Fr$ as small as $0.05$,
see [9]. The increase of $Re_{\lambda}$ will further decrease $Fr$ at which
the theory applies.

We saw that increase of $Re_{\lambda}$ makes our calculation of $\lambda_{1}$
inconsistent starting from $Re_{\lambda}$ for which $L_{c}\sim g\tau^{2}$.
However Eqs. (14)-(15) could also become invalid because of breakdown of the
assumed Gaussianity of $s_{\Delta t}(t)$ in Eq. (11), despite $\Delta t$ being
much larger than the correlation time of $s(t)$. This could happen because
cumulants of order higher than two, which the Gaussian approximation neglects
[57], involve higher order moments of the velocity gradients. These moments
due to intermittency would contain higher powers of $Re_{\lambda}$. This would
invalidate their discarding at large enough $Re_{\lambda}$. The resulting
criterion is similar to Eq. (23) and brings the same conclusion that the
derivation breaks down at increasing $Re_{\lambda}$.

We describe how the predictions of the effective white-noise description are
used in practice assuming that the validity conditions hold and $St\sim 1$.
The first step is to derive $D$ from the spectrum of turbulence by using Eq.
(14) and check if $D\tau\ll 1$. If yes, then $D$ is the predicted value of the
Lyapunov exponent. If not, then the white-noise description fails and other
treatment is necessary. Thus if we consider $\lambda_{1}\tau$ as a function of
$Re_{\lambda}$ at fixed $Fr\ll 1$ and $St\sim 1$, then $\lambda_{1}\tau$
increases with $Re_{\lambda}$ until $\lambda_{1}\tau$ becomes of order one. It
is seen that this happens when the correlation time
$L_{c}/(g\tau)\sim\lambda_{1}\tau_{K}^{2}$ of $s(t)$ is of order $\tau_{K}$.
It seems that at further increase of $Re_{\lambda}$ the correlation time,
which cannot be larger than $\tau_{K}$, gets fixed at the Lagrangian
correlation time $\tau_{K}$ and the statistics of $s(t)$ becomes similar to
the isotropic statistics of turbulent velocity gradients. This would lift
anisotropy of particles’ distribution in space and make columnar structures
disappear.

We conclude that the above considerations indicate that increase of
$Re_{\lambda}$ would smear the columnar structure and make it to disappear
altogether at the high $Re_{\lambda}$ holding in clouds. In fact we
demonstrate below that $L_{c}$ and $g\tau^{2}$ provide horizontal and vertical
dimensions of the particles’ columns respectively. Here $L_{c}$ must be
considered as the correlation length of the velocity gradients of the
particles’ flow which need not and is not the same as its counterpart for
turbulence. Starting from $Re_{\lambda}$ for which $L_{c}\sim g\tau^{2}$,
which is where the calculation of $\lambda_{1}$ becomes invalid, vertical and
horizontal dimensions of the columns become similar, the particles’ structures
are isotropic and columns are no longer preferential. Numerical studies of
dependence of $L_{c}$ on $Re_{\lambda}$ would provide further insight into the
discussion of this section and are left for future work.

##  IV Predictions for $Fr\to 0$ limit and confirmation

In this section we revisit predictions of [9] preparing the ground for the
study of how intermittency affects validity of [9] performed in the next
Section. The main observation of [9] is that smallness of Lagrangian
acceleration of the fluid particles in comparison with $g$ results in the
smooth spatial motion of the droplets. The particles’ velocities after
transients are uniquely determined by their spatial positions on which they
depend in a differentiable manner. This conclusion is reached by estimating
the accelerations with the typical value $\epsilon_{0}^{3/4}/\nu^{1/4}$. The
smoothness holds irrespective of the Stokes number, which is significant since
droplets in the clouds often have $St\gtrsim 1$ where without gravity
particles’ motion would not be spatially smooth.

The limitations of the above observation are in the usage of typical
accelerations which might not be relevant in view of intermittency of small-
scale turbulence. Thus at $Re_{\lambda}=690$, which is much smaller than
$Re_{\lambda}$ of the clouds, the accelerations’ flatness of $55$ was observed
in [32]. It is then not obvious at all if the estimation of the role of
gravity by using typical turbulent accelerations is adequate in the clouds.

The above issue can be considered similarly to the case of negligible gravity
and $St\ll 1$. We can decompose the flow domain into the major part, where
turbulent accelerations are much smaller than gravity, and rare regions of
vigorous turbulence where turbulent accelerations are larger or comparable
with $g$. For qualitative considerations we can estimate local turbulent
acceleration of the fluid particles as $\epsilon^{3/4}/\nu^{1/4}$, cf. [33].
The local Froude number $\epsilon^{3/4}/(g\nu^{1/4})$ is small in most of the
space. In this calm major part of the space the droplets form a smooth spatial
flow as explained in [9]. In contrast, in regions with local Froude number of
order one, sling effect with particles’ jets holds. The rate of collisions is
then given by the sum of contributions of the large calm region of turbulence,
where smooth flow holds, and widely spaced regions of jets. This decomposition
is identical to that at weak gravity and small inertia, see the Introduction
and [4].

For small enough $Fr$, however large $Re_{\lambda}$ is, the collision kernel
is due to the major volume fraction with smooth flow and the theory of [9]
applies. In contrast, if we fix $Fr$ and increase $Re_{\lambda}$ then the
characteristic acceleration of relevant vortices of the smooth portion of the
flow increases until it becomes of order $g$ (the sling effect would become
non-negligible already at smaller $Re_{\lambda}$, cf. [4]). At higher
$Re_{\lambda}$ the collisions due to both smooth flow and sling effect, occur
predominantly in the rare regions of ”resonant” vortices with acceleration of
order $g$ and the theory of [9] does not apply. We provide details on this
non-commutativity of $Fr\to 0$ and $Re_{\lambda}\to\infty$ limits below.

###  IV.1 Theory of [9]: predictions, confirmation, contradiction

Figure 2: Lyapunov exponents: theory vs. measurement, (a) linear and (b) log-
log plots. The small $Fr$ theory of [9] is valid at smallest $Fr$, however
breaks down at increasing $Fr$. The linear growth $|\lambda_{3}|\tau=1.5Fr$
predicted in [9] is observed. However $\lambda_{1}$, predicted to be equal to
$-|\lambda_{3}|$, obeys a different behavior, describable by the fit
$\lambda_{1}\tau=1.5Fr^{0.9}$ (we used in the text that
$\lambda_{1}\tau=1.5Fr$ is valid at $Fr\leq 0.033$). The ratio
$(|\lambda_{3}|-\lambda_{1})/\lambda_{1}$, that vanishes in the leading order
of the theory, grows roughly linearly with $Fr$ remaining reasonably small at
$Fr=0.033$ however its value $\simeq 0.5$ at $Fr=0.05$ indicates the theory
breakdown. Similarly the compressibility ratio
$C_{KY}\equiv|\sum_{i=1}^{3}\lambda_{i}/\lambda_{3}|$ deviates from the
predicted linear dependence on $Fr$ at $Fr=0.05$, remaining, as predicted,
quite small. The exponent $\lambda_{2}$ corresponds to vertical separation of
the particles and it remains small at all $Fr$ with
$max[\lambda_{2}/\lambda_{1}]=0.05$ and
$max[\lambda_{2}/|\sum_{i=1}^{3}\lambda_{i}|]=0.11$, both attained at
$Fr=0.05$. Thus $\lambda_{2}$ can be neglected in $\sum_{i=1}^{3}\lambda_{i}$.
Compressibility of the flow grows with $Fr$ because the horizontal flow,
described by $\lambda_{1}$ and $\lambda_{3}$, becomes more compressible.
Figure 3: Shown are the DNS values of the correlation codimension $\alpha$,
obtained from numerical observations of the RDF, twice the Kaplan-Yorke
(Lyapunov) codimension $C_{KY}$, fixed from the numerical observations of the
Lyapunov exponents, and the low-$Fr$ theory prediction
$3\pi\int_{0}^{\infty}E(k)kdk/(2g)$ that gives $13Fr$. All figures in the
paper use $St=1$ and $Re_{\lambda}=70$.

We describe the results for the spectrum of the Lyapunov exponents
$\lambda_{i}$ of spatial motion of the particles in the light of more detailed
numerical simulations performed for this work, see Fig. 2. It is immediate
consequence of the white noise description introduced in the previous section
that the spectrum has time-reversal symmetry [36]; that is
$\lambda_{1}=-\lambda_{3}$ and $\lambda_{2}=0$. Here the vanishing of
$\lambda_{2}$ describes conservation of the vertical component of the
separation in the leading order, see above. This approximation provides the
linear order term of the dependence of $\lambda_{i}$ in $Fr$. The first
Lyapunov exponent vanishes at $Fr=0$ where the particles’ trajectories in
space stop to be chaotic. This is because they become effectively the
trajectories of particles sedimenting in still air. This is not completely
trivial because the particles’ velocity in the leading order is still the
local turbulent flow $\bm{u}(\bm{x}(t),t)$.

To linear order in $Fr$ the sum of the Lyapunov exponents
$\sum_{i=1}^{3}\lambda_{i}$, which is a main measure of clustering, is zero
and there is no clustering. The leading order term in
$\sum_{i=1}^{3}\lambda_{i}$ is quadratic in $Fr$ and it was derived in [9] via
the energy spectrum of the turbulent flow $E(k)$. The theoretical predictions
are confirmed in Fig. 2. The Figure also shows that
$\left|\sum_{i=1}^{3}\lambda_{i}\right|$ grows with $Fr$ not because of the
growth of $\lambda_{2}$ but rather because of growing compressibility of the
horizontal flow, $max(\lambda_{2}/|\sum_{i=1}^{3}\lambda_{i}|)\ll 1$.

The main application of the Lyapunov exponents to the distribution of the
particles is to the calculation of the Kaplan-Yorke codimension $C_{KY}$, i.e.
the difference of the space dimension and the Kaplan-Yorke dimension. It can
be seen by using $\lambda_{1}+\lambda_{2}>0$ and $\sum_{i=1}^{3}\lambda_{i}<0$
in the dimension’s definition in [35] that in our case
$C_{KY}=\left|\left(\sum_{i=1}^{3}\lambda_{i}\right)/\lambda_{3}\right|$. The
codimension in this case has interpretation of compressibility ratio and it is
small because of the flow’s weak compressibility. The results of [11] for
attractors of weakly compressible flows give that the steady state
fluctuations of the particles’ concentration are lognormal and the rest of
fractal dimensions can be obtained from the Kaplan-Yorke dimension. For
instance the so-called correlation codimension $\alpha$, that is defined as
the scaling exponent in the power-law $r^{-\alpha}$ for the RDF (probability
to find a pair of droplets at distance $r$; see the Introduction) obeys

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\alpha=2C_{KY}=\frac{2\left|\sum_{i=1}^{3}\lambda_{i}\right|}{\left|\lambda_{3}\right|},$
(24)

cf. [4]. It is significant for the sequel that this is a universal relation
that holds for any weakly compressible flow [11]. This explains the central
role of the effective spatial flow of the droplets in the theory.

We determined the RDF numerically, which was not done in [9], in order to test
Eq. (24), see below and Fig. 3. The usage in the above formula of the
(corroborated) formulas of [9] that provide $\lambda_{i}$ via $E(k)$ gives,

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\alpha=\frac{3\pi\int_{0}^{\infty}E(k)kdk}{2g},\
\ Fr\to 0.$ (25)

We find by comparison with Eqs. (14) and (22) that

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\alpha=12\lambda_{1}\tau=Fr\frac{6L_{c}}{5\eta}.$
(26)

These are remarkably simple formulas that must however be taken with a grain
of salt: it is an asymptotic result and at high $Re_{\lambda}$ its validity
may demand unpractically small $Fr$ as in the discussion of $\lambda_{1}$
above, see also below.

We observe that $\alpha$ predicted by Eq. (25) is independent of the Stokes
number: it is a property of turbulence and not of particles. The reason for
the particle size independence is that the stretching rate and compression
rate of infinitesimal volumes of particles, which ratio determines the fractal
dimensions, have identical dependence on $St$. Thus statistics of attractors
of particles with different sizes is identical. This does not tell that these
attractors coincide in space. Spatial locations of these attractors differ and
only their average properties agree.

Comparison with other works—The above size-independence explains the
observations of [22] in the relevant range; see also [10]. It was found in
[22] that at $St\gtrsim 1$ the correlation codimension is $St-$independent,
which was explained by further developing the formalism of [24]. The authors
managed to demonstrate that $\alpha=CFr$ where $C$ is a constant prefactor
which is independent of $St$. The calculation of $C$ was not provided and it
could be plausibly thought that $C\sim 1$. However [9] demonstrated that Eq.
(25) gives $C\simeq 13$ at $Re_{\lambda}=70$. This is a large value that would
be difficult to obtain without calculation. This large numerical factor
implies that clustering, whose strength is measured by $\alpha$, can be strong
at very small $Fr$.

Highest $Fr$ describable by [9]—The predictions of the $Fr\to 0$ theory allow
detailed comparison with the DNS data. It was observed in [9] that, as already
mentioned, at $Re_{\lambda}=70$, the predictions for the Lyapunov exponents
hold at $Fr\leq 0.033$. However some of them break down already at $Fr=0.05$,
see Fig. 2 where data, more detailed than in [9], are presented. At the same
time, for some quantities, the theory was found to apply at larger $Fr$ also.
This is because deviations of different quantities may compensate each other
as happens to be the case of $C_{KY}/(\lambda_{1}\tau)$ predicted in [9] to be
equal to six, independently of the details of turbulence, of $Fr$ and of $St$.
This quantity is accidentally close to the theoretical value at all $Fr$
within fifteen per cent discrepancy. The prediction $\alpha=2C_{KY}\propto Fr$
holds within four per cent discrepancy up to $Fr=0.05$, see Fig. 3. Moreover,
the data of [9] demonstrate that the theoretical predictions for $C_{KY}$ and
$\lambda_{1}$ hold at $St=1$ or $St=2$ up to $Fr\leq 0.1$ and break down only
at $Fr=0.2$. In contrast, at $St=0.5$ the prediction for $\lambda_{1}$ fails
at $Fr=0.1$. These observations are in agreement with [22] who observed
different behavior at $St\geq 1$ and $St<1$ for $Fr=0.05$.

Examination of the claim that [9] is invalid—The success of the theory of [9]
is evident from the above. However [10] observed that Eq. (25) does not work
at $St=1$ and $Fr=0.052$ and concluded that the equation is wrong. This is
despite that the equation was confirmed in [9] (strictly speaking [9]
confirmed $C_{KY}=3\pi\int_{0}^{\infty}E(k)kdk/(4g)$ which together with
$\alpha=2C_{KY}$, confirmed in the present work, validates Eq. (25) at
$Re_{\lambda}=70$). The discrepancy however is the indication of incorrect use
of Eq. (25), rather than its invalidity. The calculations of [9] are rigorous
asymptotic calculations at $Fr\to 0$ and can be checked. Their application
depends on the Reynolds number, since the theory holds non-uniformly in
$Re_{\lambda}$ as we stress in this work. In the case of [10] the discrepancy
is both due to factor of $2$ mistake in the formula and testing of the theory
at $Fr$ which was already claimed to be beyond the theory in [9]. The simple
rule of the thumb is that the theoretical prediction for $\alpha$ applies
provided that $\alpha\ll 1$, cf. below. Since $\alpha$ studied in [10] does
not obey the inequality then the observed deviation is more than reasonable,
cf. above.

##  V Theory’s limitation: clash of $Fr\to 0$ and $Re_{\lambda}\to\infty$
limits

In this section we provide detailed exposition of the theory of [9] in order
to demonstrate that increase of the Reynolds number at fixed $Fr$ invalidates
this theory. The main observation of [9] is that in the $Fr\to 0$ limit, one
can define the field $\bm{v}_{a}(t,\bm{x})$, which provides the velocity of
the droplet of radius $a$ located at time $t$ at point $\bm{x}(t)$:

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\dot{\bm{x}}=\bm{v}_{a}(t,\bm{x}(t)),\
\
\partial_{t}\bm{v}_{a}+(\bm{v}_{a}\cdot\nabla)\bm{v}_{a}=-\frac{\bm{v}_{a}-\bm{u}}{\tau(a)}+\bm{g},$
(27)

where the partial differential equation (PDE) on $\bm{v}_{a}(t,\bm{x})$ is
implied by the equation of motion on differentiating
$\dot{\bm{x}}=\bm{v}_{a}(t,\bm{x}(t))$ over time [4]. Below we omit the
subscript of $\bm{v}_{a}(t,\bm{x})$ unless the radius to which the flow
pertains needs to be referred. The flow is introduced implicitly as the
solution of the above PDE. An explicit formula for $\bm{v}(t,\bm{x})$ via
$\bm{u}(t,\bm{x})$, similar to the small Stokes numbers’ Eq. (4), is
unavailable. The solution can only become a well-defined single-valued field
after transients, lasting for times of order $\tau(a)$, during which the
initial condition is forgotten. Indeed, one can devise initial conditions for
Eq. (27) producing multi-valued flow after a short time.

It is readily seen that at finite $Fr$ the assumption of single-valued field
is inconsistent due to the blowup (sling) events at which $\nabla\bm{v}$
explodes. Indeed, spatial differentiation of Eq. (27) , and the passage in the
resulting equation to the frame moving with the particle, gives a closed
matrix ODE for the gradients
$\sigma_{ik}(t)\equiv\nabla_{k}v_{i}(t,\bm{x}(t))$ [4],

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\dot{\sigma}\\!+\\!\sigma^{2}\\!=\\!-\frac{\sigma\\!-\\!s}{\tau},$
(28)

where $s_{ik}(t)\equiv\nabla_{k}u_{i}(t,\bm{x}(t))$. If the quadratic term in
the LHS is not small then it would cause finite-time explosion of the
gradients. That explosion would signal the breakdown of the flow description
due to the flow becoming multi-valued [4, 14]. In contrast, if $\sigma^{2}$ in
Eq. (28) is small then, in the leading order, the gradients of the droplets’
flow in the particle frame are given by the finite expression,

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\sigma\approx\sigma_{l},\
\
\sigma_{l}\equiv\int_{-\infty}^{t}\exp\left(-\frac{t-t^{\prime}}{\tau}\right)s(t^{\prime})\frac{dt^{\prime}}{\tau},$
(29)

where the subscript stands for linearization of Eq. (28).

We first disregard intermittency as in our study of the Lyapunov exponent. As
in that study, the correlation time of $s_{ik}(t)$ is given by the smallest of
the Kolmogorov time $\tau_{K}$ and the sedimentation time
$\tau_{g}=\eta/g\tau$. We find as previously that the time is $\tau_{g}$ due
to $\tau_{g}/\tau_{\eta}=Fr/St\ll 1$. Moreover, we have
$\tau_{g}/\tau=Fr/St^{2}\ll 1$; therefore, the effective integration interval
$\tau$ in Eq. (29) is much larger than the correlation time of $s(t)$. Thus
$\sigma$ is effectively a sum of large number of independent identically
distributed random variables and is Gaussian, cf. the study of $\lambda_{1}$.

The Gaussianity implies that the condition that the probability of sling
events is small is tantamount to the condition that the dispersion
$\langle\sigma_{l}^{2}\rangle$ is much smaller than $1/\tau^{2}$ since the
bulk of the probability is determined by the dispersion (in writing matrix as
a scalar, as $\sigma_{l}^{2}$, the characteristic value is implied). Here we
disregard the non-zero average of $s_{ik}$ that exists due to the combination
of preferential concentration and anisotropy. This average can be excluded by
considering $\sigma_{l}-\langle\sigma_{l}\rangle$ instead of $\sigma_{l}$ and
observing that smallness of $\langle\sigma_{l}\rangle$ makes it irrelevant in
the estimates here and below, see Appendix A. Averaging the square of Eq.
(29), it is found that $\langle\sigma_{l}^{2}\rangle\tau^{2}\sim Fr$, see [9]
and Eq. (40) below. Therefore, in the leading order in $Fr\ll 1$, the
nonlinear term in Eq. (28) can be self-consistently neglected. The flow is
single-valued and well-defined in most of the space. This does not guarantee
though that the rare regions of slings cannot provide appreciable contribution
to some quantities, see the Introduction and [4].

We observe that the effective domain of integration in Eq. (29) is
$(t-\tau,t)$ so that
$\sigma_{l}(t)\tau\\!\sim\int_{t-\tau}^{t}\\!s(t^{\prime})dt^{\prime}$. Then
the considerations of the last paragraph imply that $\sigma_{l}$ is
statistically similar to instantaneous turbulent flow gradient $\nabla\bm{u}$
averaged over spatial interval of order $g\tau^{2}$ passed by sedimenting
particle in time $\tau$. This object resembles a similar average of local
energy dissipation rate, which is a standardly measured quantity [47]. The
measurements of [48] indicate that this quantity must undergo strong
fluctuations in clouds. The fluctuations, that are weak at $Re_{\lambda}=70$
considered in [9], destroy the particles’ flow at larger Reynolds numbers.

###  V.1 Flow breakdown at increasing $Re_{\lambda}$

The Kolmogorov-type estimates used after Eq. (29) are changed profoundly by
intermittency. Below we perform order of magnitude calculations only, not
writing the matrix indices.

Condition for well-defined flow of the particles—It is immediate consequence
of our study in Sec. III that the relevant spatial scale in the calculation of
$\langle\sigma_{l}^{2}\rangle$ is $L_{c}$ and not $\eta$. We have assuming
that $L_{c}/(g\tau)\lesssim\tau$, cf. Eq. (23)

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\langle\sigma_{l}^{2}\rangle\tau^{2}\sim\int_{t-\tau}^{t}dt_{1}dt_{2}\langle
s(t_{1})s(t_{2})\rangle\sim\frac{\tau}{\tau_{K}^{2}}\frac{L_{c}}{g\tau}=Fr\frac{L_{c}}{\eta}\sim
15Fr\left(\frac{Re_{\lambda}}{70}\right)^{\Delta_{\lambda}},$ (30)

where we observed that $\int\langle s(0)s(t)\rangle dt$ is estimated as the
product of characteristic value $\tau_{K}^{-2}$ of $s^{2}$ and the correlation
time $L_{c}/(g\tau)$ and we used notations of Sec. III. We observe that for
$St\sim 1$ the condition $\langle\sigma_{l}^{2}\rangle\tau^{2}\ll 1$ coincides
with the validity condition of the calculation of $\lambda_{1}$ given by Eq.
(23). Since the last condition is $L_{c}/(g\tau)\ll\tau$ then our assumption
$L_{c}/(g\tau)\lesssim\tau$ made in the calculation for Eq. (30) is self-
consistent. If $\langle\sigma_{l}^{2}\rangle\tau^{2}\ll 1$ holds then, by
Chebyshev’s inequality, the probability of the sling effect,
$|\sigma_{l}|\tau\sim 1$, is small and motion in most of the space is smooth,
cf. the discussion after Eq. (23). In contrast, at high $Re_{\lambda}$ obeying
$FrRe_{\lambda}^{\Delta_{\lambda}}\gtrsim 1$, the sling effect is typical. Due
to $St\sim 1$ the particle strongly interacts with vortex correlated over the
correlation length $L_{c}$ and the flow description breaks down. We consider
it highly plausible that this transition happens at $Re$ lower than those in
the clouds.

###  V.2 Strong anisotropy of the particles’ flow at small scales

We have introduced in subsection III.2 the scale $L_{c}$ which is the spatial
scale of turbulent flow gradients that determine the Lyapunov exponent. We
indicated that this scale must be understood as the correlation length of the
particles’ flow gradients. Here we refine this statement by demonstrating that
this is the correlation length only for generically oriented distances between
the points where correlations are studied. For vertical distances the
correlations are of longer range. As a result gradients of the particles’ flow
are correlated over spatial region similar to a vertically oriented column.
This anisotropy leaves direct mark on the multifractal distribution of
particles, considered later, where the multifractal structures extend over
similar columns.

The flow of the particles integrates non-trivially the turbulent flow via Eq.
(27). This implies that spatio-temporal properties of $\bm{v}(t,\bm{x})$ and
$\bm{u}(t,\bm{x})$ are quite different. We consider pair correlations of
$\nabla_{i}v_{k}(\bm{x})$ at $t=0$. We have according to Eq. (29) that

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\nabla_{i}v_{k}(\bm{x})\approx\int_{-\infty}^{0}\exp\left(\frac{t}{\tau}\right)s_{ik}(t)\frac{dt}{\tau},\
\ s_{ik}(t)\equiv\nabla_{k}u_{i}(t,\bm{q}(t,\bm{x})),$ (31)

where we introduced the Lagrangian trajectories of the particles’ flow
$\bm{q}(t,\bm{x})$. These are particles’ trajectories that pass $\bm{x}$ at
$t=0$. The flow’s existence implies that these trajectories are unique and
solve

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\partial_{t}\bm{q}(t,\bm{x})=\bm{v}(t,\bm{q}(t,\bm{x})),\
\ \bm{q}(t=0,\bm{x})=\bm{x}.$ (32)

We have

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\left\langle\nabla_{i}v_{k}(0)\nabla_{l}v_{m}(\bm{x})\right\rangle\approx\int_{-\infty}^{0}\exp\left(\frac{t_{1}+t_{2}}{\tau}\right)\left\langle\nabla_{k}u_{i}(t_{1},\bm{q}(t_{1},0))\nabla_{l}u_{m}(t_{2},\bm{q}(t_{2},\bm{x}))\right\rangle\frac{dt_{1}dt_{2}}{\tau^{2}}.$
(33)

The exponential decay factor in above integral effectively cuts the
integration domain off at times $|t_{i}|$ smaller or of order $\tau$. However
in the studied regime the particles’ separation during times of order $\tau$
is negligible due to $\lambda_{1}\tau\sim Fr\ll 1$, see Section III. Thus we
have $\bm{q}(t,\bm{x})-\bm{q}(t,0)\approx x$ for $|t|\lesssim\tau$. The
trajectories’ positions with respect to each other are determined by the
gravitational sweeping, similarly to the previous sections, and we may write

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\left\langle\nabla_{i}v_{k}(0)\nabla_{l}v_{m}(\bm{x})\right\rangle\approx\int_{-\infty}^{0}\exp\left(\frac{t_{1}+t_{2}}{\tau}\right)\left\langle\nabla_{k}u_{i}(\bm{g}\tau
t_{1})\nabla_{l}u_{m}(\bm{x}+\bm{g}\tau
t_{2})\right\rangle\frac{dt_{1}dt_{2}}{\tau^{2}},$ (34)

where the angular brackets here stand for usual spatial averaging. We see
that, provided that $|\bm{x}|$ is much smaller the correlation length $L_{c}$
of typical turbulent velocity gradients that determine the above average, we
can omit $\bm{x}$ in the above equation. We find that
$\left\langle\nabla_{i}v_{k}(0)\nabla_{l}v_{m}(\bm{x})\right\rangle\approx\left\langle\nabla_{i}v_{k}(0)\nabla_{l}v_{m}(0)\right\rangle$
for $x\ll L_{c}$. In contrast, if $x\gg L_{c}$ then, for a generic direction
of $\bm{x}$, the arguments of $\nabla\bm{u}$ in Eq. (34) are separated at all
relevant $t_{i}$ by distance much larger than $L_{c}$ and the correlation
function $\left\langle\nabla_{i}v_{k}(0)\nabla_{l}v_{m}(\bm{x})\right\rangle$
is much smaller than
$\left\langle\nabla_{i}v_{k}(0)\nabla_{l}v_{m}(0)\right\rangle$. There is
however an exception. If the ”initial” separation $\bm{x}$ is vertical and
much smaller than $g\tau^{2}$ then within time much smaller than $\tau$ the
argument $\bm{x}+\bm{g}\tau t_{2}$ becomes similar to $\bm{g}\tau t_{2}$ and
we still have
$\left\langle\nabla_{i}v_{k}(0)\nabla_{l}v_{m}(\bm{x})\right\rangle\approx\left\langle\nabla_{i}v_{k}(0)\nabla_{l}v_{m}(0)\right\rangle$.
It is readily seen that these considerations extend to those $\bm{x}$ that,
besides the vertical component smaller than $g\tau^{2}$, have a non-zero
horizontal component that is smaller than $L_{c}$. We conclude that velocity
gradients are correlated over cylinders (columns) with height of order
$g\tau^{2}$ and radius of the base of order $L_{c}$ (of course the shape is
not sharply defined). The cylinder is elongated, with height larger than
diameter by a large factor studied in detail later.

The significant difference in spatial correlations of flows of the fluid and
of the particles indicates that most probable gradients of droplets’ flow are
not created by the most probable gradients of turbulence, because the
connection between the gradients is not local in time. Rather gradients of the
particles’ flow are created by eddies with local viscous scale of order
$L_{c}$.

###  V.3 Breakdown of Gaussianity due to intermittency

The other effect of the intermittency is that it destroys Gaussianity of
$\sigma_{l}$. This can occur when the probability of the sling effect, as
measured by $\langle\sigma_{l}^{2}\rangle\tau^{2}$, is still small,
$\langle\sigma_{l}^{2}\rangle\tau^{2}\ll 1$. In that case the flow description
applies however the flow gradients are non-Gaussian. We demonstrate breakdown
of Gaussianity by considering the fourth moment. We have using Eq. (29)

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\langle\sigma_{l}^{4}\rangle\tau^{4}\sim
3\int_{t-\tau}^{t}\\!\langle s(t_{1})s(t_{2})\rangle\langle
s(t_{3})s(t_{4})\rangle dt_{1}dt_{2}dt_{3}dt_{4}+\int_{t-\tau}^{t}\\!\langle
s(t_{1})s(t_{2})s(t_{3})s(t_{4})\rangle_{c}dt_{1}dt_{2}dt_{3}dt_{4},$ (35)

where we decomposed the fourth-order correlation function of $s(t)$ into the
sum of the reducible part, given by the first term in the RHS, and the
irreducible part or the cumulant [57], designated by the subscript $c$ and
given by the last term. We have using properties of the different-time
correlations of $s(t)$, described previously, the estimate

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\int_{t-\tau}^{t}\langle
s(t_{1})s(t_{2})\rangle\langle s(t_{3})s(t_{4})\rangle
dt_{1}dt_{2}dt_{3}dt_{4}\sim Fr^{2}Re_{\lambda}^{2\Delta_{\lambda}}.$ (36)

If intermittency is negligible (so that $Re_{\lambda}^{\Delta_{\lambda}}\sim
1$ and the last term in Eq. (35) obeys $Re_{\lambda}-$independent estimate
$\tau(\eta/g\tau)^{3}/\tau_{K}^{4}\sim Fr^{3}/St^{2}$) then the last term in
Eq. (35) would be smaller than the term coming from the reducible part of the
correlation function by the factor of $\eta/(g\tau^{2})\sim Fr/St^{2}$.
Indeed, the irreducible correlation function $\langle
s(t_{1})s(t_{2})s(t_{3})s(t_{4})\rangle_{c}$ decays quickly when the
separation of any pair of $t_{i}$ becomes larger than the correlation time
$\eta/(g\tau)$. Thus the effective domain of integration over $t_{i}$ is
smaller than for the reducible contribution by $Fr/St^{2}$, the fact that lies
at the origin of the Gaussianity considered previously.

The intermittency causes the probability of $s$ much larger than the rms value
$\tau_{K}^{-1}$ to be non-negligible. Thus equal-time fourth order moment of
$s$ is proportional to $\tau_{K}^{-4}$ times a positive power of
$Re_{\lambda}$. As a result, at large $Re_{\lambda}$, the fourth order-
cumulant is larger than the reducible contribution $3\langle
s^{2}\rangle^{2}\sim\tau_{K}^{-4}$ by a power of the Reynolds number. Since
correlation time also depends on $Re_{\lambda}$ then it is not completely
obvious, however very plausible, that the last term of the LHS of Eq. (35)
contains a higher power of $Re_{\lambda}$ than
$Re_{\lambda}^{2\Delta_{\lambda}}$ of the reducible contribution. We obtain

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\frac{\int_{t-\tau}^{t}\langle
s(t_{1})s(t_{2})s(t_{3})s(t_{4})\rangle_{c}dt_{1}dt_{2}dt_{3}dt_{4}}{\int_{t-\tau}^{t}\langle
s(t_{1})s(t_{2})\rangle\langle s(t_{3})s(t_{4})\rangle
dt_{1}dt_{2}dt_{3}dt_{4}}\sim\frac{FrRe_{\lambda}^{\Delta}}{St^{2}},$ (37)

where the exponent $\Delta$ is positive and absorbs all dependencies on
$Re_{\lambda}$. Then the increase of $Re_{\lambda}$ at fixed $Fr$ and $St$
causes the irreducible contribution in $\langle\sigma_{l}^{4}\rangle$ to
become dominant. This happens at $Re_{\lambda}^{\Delta}\sim St^{2}/Fr$ which
at $St\sim 1$ occurs when the flow is still well-defined provided that
$\Delta>\Delta_{\lambda}$, cf. the comment after Eq. (30). Thus in this case
we will have a flow with non-Gaussian statistics of the gradients similarly to
the turbulence itself. In contrast, if $\Delta<\Delta_{\lambda}$ then, as long
as the flow is well-defined at all, the gradients’ statistics is Gaussian.
Similar observations hold for higher order moments
$\langle\sigma_{l}^{2n}\rangle$ with $n>2$. The converse of the above
statements is that decrease of $Fr$, at fixed $St$ and however large
$Re_{\lambda}$, causes the statistics of $\sigma_{l}$ to become Gaussian and
makes the flow description valid.

##  VI Radial distribution function at $Fr\to 0$: spatial correlations and sum
rule

We have considered above the general properties of the flow of the particles
and how increase in the intermittency due to growth of the Reynolds number can
change its properties. The rest of the paper is devoted to the case of $Fr$ so
small that the flow description is valid. This range depends on $St$ and
$Re_{\lambda}$ as explained above.

In this section we consider the radial distribution function (RDF). Since the
flow is weakly compressible then, generally, the RDF obeys scaling below the
correlation length of the gradients of the droplets’ flow [11]. We saw in
subsection V.2 that the gradients are correlated over columns whose height is
of order $g\tau^{2}$ and diameter of order of $L_{c}$. Here we demonstrate
numerically that, indeed, the scaling of the RDF holds over these columns, all
of whose dimensions are much larger than the Kolmogorov scale $\eta$. This
effect is solely due to gravity: multifractality at negligible gravity and
small inertia holds below the same correlation length $\eta$ as the turbulent
velocity gradients have [4]. We also derive here a sum rule demonstrating that
despite that the RDF, due to strong anisotropy, has a non-trivial angular
dependence, still the angle-averaged RDF is rather simple.

###  VI.1 Small-scale flow is two-dimensional

Main properties of the droplets’ flow at $Fr\to 0$ are weak compressibility
and small scale two-dimensionality. The weak compressibility can be seen by
taking the trace of Eq. (29) and using $trs=0$. This gives that in the leading
order $tr\sigma=0$ and the droplets’ flow is incompressible. Clustering occurs
due to small compressible component of the flow. The leading order term for
$tr\sigma$ can be found by taking the trace of Eq. (28). The solution of the
obtained equation gives (this equation in [9] contains a typo with a wrong
sign that is irrelevant for the calculations there)

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!tr\sigma(t)=-\int_{-\infty}^{t}dt^{\prime}\exp\left(-\frac{t-t^{\prime}}{\tau}\right)tr\sigma^{2}(t^{\prime})dt^{\prime}s\approx-\int_{-\infty}^{t}dt^{\prime}\exp\left(-\frac{t-t^{\prime}}{\tau}\right)tr\sigma_{l}^{2}(t^{\prime})dt^{\prime},$
(38)

where the last term is estimated as $\tau\sigma_{l}^{2}$. We find that
compressibility ratio $tr\sigma/\sigma_{l}$ is of order $\tau\sigma_{l}\ll 1$.

Droplets’ flow at small scales is described by the Taylor expansion of the
velocity difference below the smoothness scale $\eta_{p}(\theta)$ of the
particles’ flow ($\theta$ is the polar angle of $\bm{r}$, with $z-$axis
directed upwards)

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\bm{v}(\bm{x}_{2})-\bm{v}(\bm{x}_{1})\approx(\bm{r}\cdot\nabla)\bm{v}(x_{1}),\
\ r\ll\eta_{p}(\theta),\ \ \bm{r}\equiv\bm{x}_{2}\bm{x}_{1}.$ (39)

Here we introduced the angle-dependent correlation length $\eta_{p}(\theta)$.
The study that we performed in subsection V.2 implies that
$\eta_{p}(\theta)\sim L_{c}$ for $\theta\gg L_{c}/(g\tau^{2})$ and
$\eta_{p}(\theta)\sim g\tau^{2}$ for $\theta\lesssim L_{c}/(g\tau^{2})$. Here
$L_{c}/(g\tau^{2})\ll 1$, see below, and $0\leq\theta\leq\pi/2$ (the range of
$\pi/2\leq\theta\leq\pi$ is obtained by interchange of $\bm{x}_{i}$ in Eq.
(39)).

The matrix $\nabla_{i}v_{k}$ is a Gaussian $2\times 2$ matrix in the leading
order in $Fr$. This observation was made in [9] for the Lagrangian statistics
in the particle frame, see also Sec. III. We demonstrate in Appendix B that
the same proof applies for the Eulerian statistics at a fixed spatial point.
One finds that $\nabla\bm{v}$ does not have vertical components or derivatives
of horizontal components in the vertical direction in the leading order in
$Fr$. The statistics of planar gradients is Gaussian so that they are
described fully by zero mean and the dispersion,

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\langle\nabla_{k}v_{i}\nabla_{r}v_{p}\rangle=\frac{c_{0}Fr\left(3\delta_{ip}\delta_{kr}-\delta_{ik}\delta_{pr}-\delta_{ir}\delta_{kp}\right)}{2\tau^{2}},$
(40)

where $c_{0}$ is defined in Eq. (16). This formula confirms the order of
magnitude estimate $\langle(\nabla\bm{v})^{2}\rangle\tau^{2}\sim Fr$ made
above. Taking its trace over $x$ and $y$ gives
$\langle\left(\nabla_{x}v_{x}+\nabla_{y}v_{y}\right)^{2}\rangle=0$ so that the
planar flow is incompressible. There are only three random variables that
describe the gradients given by the components of the zero trace $2\times 2$
matrix. We stress that these properties describe small scales only: the
absolute motion of the droplet is three-dimensional, resembling that of the
tracer, by our assumption that the droplet velocity is dominated by large-
scale turbulence.

###  VI.2 RDF

The RDF, designated by $g_{12}(\bm{r})$, determines the probability density
function (PDF) $P_{12}(\bm{r})$ of finding a droplet of radius $a_{2}$ at a
distance $\bm{r}$ from a droplet with radius $a_{1}$. The PDF is proportional
to the mean concentration $\langle n_{2}\rangle$ so that we can write
$P_{12}(\bm{r})=\langle n_{2}\rangle g_{12}(\bm{r})$ where $g_{12}(\bm{r})$ is
a property of the two-particle motion. Here and below $n_{i}(\bm{x},t)$ is the
concentration of droplets with radius $a_{i}$. The RDF also provides the pair-
correlation function of concentrations, $\langle
n_{1}(0)n_{2}(\bm{r})\rangle=\langle n_{1}\rangle\langle n_{2}\rangle
g_{12}(\bm{r})$, see Appendix C.

Pair correlations of concentration $\langle n_{1}(0)n_{2}(\bm{r})\rangle$ are
determined quite similarly to the pair correlation function of passive scalar
in incompressible turbulence in the presence of the scalar source [36]. The
reason is that the continuity equation, obeyed by the droplets’ concentration
$n$, can be written as,

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\partial_{t}\ln
n+(\bm{v}\cdot\nabla)\ln n=-\nabla\cdot\bm{v},$ (41)

We see that the flow’s divergence, whose correlation length is $\eta_{p}$, can
be considered as the source of concentration fluctuations. In the case of
passive scalar the pair correlation is determined by the time that the
particles spent below the correlation length of the source before the
observation, when the correlations accumulate [36]. Similarly here $\langle
n_{1}(0)n_{2}(\bm{r})\rangle$ is determined by the mean duration
$t^{*}_{12}(\bm{r})$ of the time interval that the particles spent below the
correlation length of the source $\eta_{p}$ before they approach by distance
$r$ where $r\ll\eta_{p}$ (here $\eta_{p}$ is much smaller than the vessel size
which is of order of integral scale of turbulence or larger).

The relevance of $t^{*}_{12}(\bm{r})$ is seen by considering motion of only
two droplets of radii $a_{1}$ and $a_{2}$ in a finite-size vessel over long
time. The droplets’ separation is most of the time of order of the size of the
vessel and is much larger than $\eta_{p}$. At these scales the flow is
incompressible so that there is no build-up of correlations - the particles
approach or separate with the same probability without creating large-scale
inhomogeneities of the concentration. However occasionally the particles
approach by a small distance $\bm{r}$. This happens after they approach by
distance $\eta_{p}$. The mean time between the last two events is
$t^{*}_{12}(\bm{r})$ and during this time-interval the particles have
effective attraction because they experience the same local gradient of the
flow whose divergence, on average, is negative, which leads to correlations’
growth [38].

The explicit representation of $g_{12}(\bm{r})$ via $t^{*}_{12}(\bm{r})$ is
found by starting from straightforward generalization of the formula of [11]
that uses cumulant expansion theorem and weakness of compressibility

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\ln
g_{12}(\bm{r})=\int_{-\infty}^{0}dt_{1}dt_{2}\langle\nabla\cdot\bm{v}_{1}(t_{1},\bm{q}_{1}(t_{1},0))\nabla\cdot\bm{v}_{2}(t_{2},\bm{q}_{2}(t_{2},\bm{r}))\rangle,$
(42)

where $\bm{v}_{i}(t,\bm{x})$ is short-hand notation for
$\bm{v}_{a_{i}}(t,\bm{x})$, see Eq. (27), and $q_{i}(t,\bm{x})$ are the
Lagrangian trajectories of the flow $\bm{v}_{i}(t,\bm{x})$, see Eq. (32). Here
in the leading order in the small compressibility we can use for $\bm{q}_{i}$
the Lagrangian trajectories of the solenoidal component of $\bm{v}_{i}$ which
guarantees that
$\left\langle\nabla\cdot\bm{v}_{i}(t,\bm{q}_{i}(t,\bm{r}))\right\rangle=0$
(otherwise dispersion must be used in the above equation, see [11]). We
consider first the case of $a_{1}=a_{2}$. This is the only case where
$g_{12}(0)$ is infinite because the trajectories $q_{i}(t,0)$ coincide causing
divergence of the integral in Eq. (42) at $r=0$. We find from Eq. (42) by
using the formula of [61] for the sum of Lyapunov exponents $\sum\lambda_{i}$
at weak compressibility that [11]

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\ln
g(\bm{r})=2\left|\sum\lambda_{i}\right|t^{*}(\bm{r}),\ \
2\left|\sum\lambda_{i}\right|=\int_{-\infty}^{\infty}\langle\nabla\cdot\bm{v}(t,\bm{q}(0,0))\nabla\cdot\bm{v}(t,\bm{q}(t,0))\rangle
dt,$ (43)

where $g(\bm{r})=g_{11}(\bm{r})$ and
$\bm{q}(t,\bm{r})\equiv\bm{q}_{1}(t,\bm{r})$. Here the time
$t^{*}(\bm{r})\equiv t^{*}_{11}(\bm{r})$ can be defined as the time of
backward in time divergence of trajectories starting at distance $r$ to a
distance $\eta_{p}$, which is equivalent to the definition above. Beyond
$t^{*}(\bm{r})$ velocity divergence factors in Eq. (42) decorrelate and the
integral converges rapidly. We use that the PDF of $t^{*}(\bm{r})$ is strongly
peaked at the average value so the fluctuations of $t^{*}(\bm{r})$ can be
neglected, see [11] and below. We also observe that since the product of
concentration and infinitesimal volume is the conserved mass then
$|\sum\lambda_{i}|$ is also asymptotic growth rate of the concentration. This
can be used to show that $g(\bm{r})$ equals squared density accumulated by
compression of initially unit density from volume of size $\eta_{p}$ during
time $t^{*}(\bm{r})$, see [4, 9].

Since the divergence of the trajectories backward in time is horizontal, then
here $\eta_{p}$ is an appropriately defined correlation length $l_{c}$ of
$\nabla\cdot\bm{v}$ in the horizontal direction. We saw that $l_{c}$ must be
of order $L_{c}$ and we keep it as a free parameter for fitting numerical
results. The time $t^{*}(\bm{r})$ diverges logarithmically at small $r$ since
the divergence of the trajectories is exponential. Therefore, the divergence
of $g(r)$ at small $r$ is a power-law, see Eq. (43). For finding this law we
use that the fluctuations of time $t^{*}(\bm{r})$ defined as solution of
$|\bm{q}(t,\bm{r})-\bm{q}(t,0)|=l_{c}$ are weak. Indeed, similarly to the
limit for the Lyapunov exponent in a forward in time flow, we also have a
realization-independent limit for the first Lyapunov exponent of the time-
reversed flow $\lim_{t\to-\infty}\lim_{r\to
0}|t|^{-1}\ln(|\bm{q}(t,\bm{r})-\bm{q}(t,0)|/r)$. This logarithmic rate of
backward in time separation can be obtained by considering time-reversed
evolution of small volumes of particles, and using the compressibility
smallness. It is found that the rate equals the same $|\lambda_{3}|$ that
characterizes the rate of particles approach forward in time, see details in
[11, 36, 9, 38, 36]. Thus at large finite $|t|$ we have, with probability
close to one,
$|t|^{-1}\ln(|\bm{q}(t,\bm{r})-\bm{q}(t,0)|/r)\approx|\lambda_{3}|$ (we saw
previously that $|\lambda_{3}|\approx\lambda_{1}$). Moreover, in the leading
order the separation is two-dimensional and is determined by the horizontal
component $r\sin\theta$ of the initial separation (the vertical component
$r\cos\theta$ is conserved). We find at $t<0$ the separation law
$|t|^{-1}\ln(|\bm{q}(t,\bm{r})-\bm{q}(t,0)|/r\sin\theta)\approx|\lambda_{3}|$.
This gives $t^{*}(\bm{r})\approx|\lambda_{3}|^{-1}\ln(l_{c}/r\sin\theta)$. We
find from Eq. (43)

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!g(\bm{r})=\left(\frac{l_{c}}{r\sin\theta}\right)^{\alpha},\
\ \theta\gg\theta^{*},\ \
\alpha=2\left|\frac{\sum_{i=1}^{3}\lambda_{i}}{\lambda_{3}}\right|=\frac{3\pi\int_{0}^{\infty}E(k)kdk}{2g}=12c_{0}Fr.$
(44)

Here we introduced a lower threshold for the above considerations, a small
angle $\theta^{*}$. If the backward in time separation of $\bm{q}(t,\bm{r})$
and $\bm{q}(t,0)$ was rigorously two-dimensional, the form of $g(\bm{r})$ in
Eq. (44) would hold at any $\theta$, producing divergence at $\theta=0$. In
fact the next order corrections make the separation weakly three-dimensional.
Thus particles that are initially separated vertically disperse due to the
components of $\sigma$ neglected in Eq. (40). This regularizes the divergence
and introduces the lower limit $\theta^{*}$ of applicability in Eq. (44). The
small angle $\theta^{*}$, studied in section VII, is defined by the condition
that the separation velocities due do horizontal and vertical components of
the separation vector are of the same order when the particle separation forms
angle $\theta^{*}$ with the vertical. Eq. (44) was derived in [9] using the
smoothness scale of turbulence, the viscous scale $\eta$, instead of the
smoothness scale of the droplets’ flow $l_{c}$. The difference is
insignificant at $Fr\to 0$ where $\alpha$ is so small that
$(l_{c}/\eta)^{\alpha}\approx 1$ and the formulas coincide. However when
larger $Fr$ are considered, the difference becomes relevant, see below. We
stress that the formula for $\alpha$ is independent of the particle size and
holds for droplets driven by the Navier-Stokes turbulence without
approximations [9]. In this work we provide the numerical confirmation of the
prediction (which was not done in [9]), see Fig. 4. We also observe the large
numerical factor of $12$ above.

Figure 4:  Angle dependence of the RDF at fixed distance $r$ between the
particles deep in the smoothness range, (a) normalized RDF and (b)
$\sin\theta^{*}(Fr)$. The smoothness range is about $10\eta$, and the
considered range of $r$ was an order of magnitude smaller. The calculations
were performed separately for $r$ in small windows centered around different
values of $r$, and the results for different $r$ produced identical angle
dependence. Numerical observations for five different $Fr$ at fixed $St=1$ are
presented by different colors and compared with the $\sin^{-\alpha}\theta$ law
(solid line), that holds in the range of angles near $\pi/2$ as predicted
asymptotically in [9]. The indefinite growth of the RDF toward smaller
$\theta$ is stopped at angles where separation rates dictated by the vertical
and horizontal components are of the same order. This causes flattening of the
RDF at the characteristic angle $\theta^{*}$ determined from
$g(r,0)/g(r,\pi/2)=\sin^{-\alpha}\theta^{*}$ studied in section VII.

###  VI.3 Correlation dimension and angle-averaged RDF: sum rule

The angle-averaged (equal-size) RDF can be derived from the result of [9] that
the concentration $n_{r}$ coarse-grained over scale $r$, which is much smaller
than $l_{c}$, obeys

$\displaystyle\frac{\langle n_{r}^{2}\rangle}{\langle
n\rangle^{2}}=\left(\frac{l_{c}}{r}\right)^{\alpha},\ \ r\ll l_{c};\ \
n_{r}(\bm{x})=\int_{|\bm{x}^{\prime}-\bm{x}|<r}n(\bm{x}^{\prime})\frac{d\bm{x}^{\prime}}{4\pi
r^{3}/3}.$ (45)

This result holds provided that $\alpha\ll 1$ so that concentration
fluctuations at $r=l_{c}$ are weak [11]. It derives from the observation that
mass $\int_{|\bm{x}^{\prime}-\bm{x}|<r}n(\bm{x}^{\prime})d\bm{x}^{\prime}$
observed in the ball at $t=0$ is compressed by the flow effectively from an