# Quality Modeling Under A Relaxed Natural Scene Statistics Model

Abhinau K. Venkataramanan  Department of Electrical and Computer Engineering
The University of Texas at Austin
Austin, USA
<EMAIL_ADDRESS>Alan C. Bovik  Department of Electrical and Computer
Engineering
The University of Texas at Austin
Austin, USA
<EMAIL_ADDRESS>
###### Abstract

Information-theoretic image quality assessment (IQA) models such as Visual
Information Fidelity (VIF) and Spatio-temporal Reduced Reference Entropic
Differences (ST-RRED) have enjoyed great success by seamlessly integrating
natural scene statistics (NSS) with information theory. The Gaussian Scale
Mixture (GSM) model that governs the wavelet subband coefficients of natural
images forms the foundation for these algorithms. However, the explosion of
user-generated content on social media, which is typically distorted by one or
more of many possible unknown impairments, has revealed the limitations of
NSS-based IQA models that rely on the simple GSM model. Here, we seek to
elaborate the VIF index by deriving useful properties of the Multivariate
Generalized Gaussian Distribution (MGGD), and using them to study the behavior
of VIF under a Generalized GSM (GGSM) model.

###### Index Terms:

Visual Information Fidelity, Generalized Gaussian Scale Mixture, Differential
Entropy, Kurtosis.

##  I Introduction

The field of Natural Scene Statistics (NSS) modeling of images originates from
the seminal work of Ruderman [1], later formalized in the form of the Gaussian
Scale Mixture (GSM) model that governs the subband coefficients of natural
images [2]. The GSM model forms the backbone of successful Full-Reference (FR)
Image Quality Assessment (IQA) models such as Information Fidelity Criterion
(IFC) [3], Visual Information Fidelity (VIF) [4], and VMAF [5], as well as
Reduced Reference (RR) Video Quality Assessment (VQA) models such as ST-RRED
[6] and SpEED-QA [7]. The GSM model also underlies the perceptually
significant divisive normalization processess used in popular No-Reference
(NR) IQA models like DIIVINE [8], BRISQUE [9], and NIQE [10].

However, images may undergo distortions in multiple stages (e.g. acquisition
and compression) or may coincidentally occur, commingling to form aggregate
distortions that are hard to model or even characterize. Distortions such as
these often arise in User Generated Content (UGC) [11] [12], where combined
degradations are hard to capture using the GSM model. A comprehensive
evaluation of NR IQA models on UGC can be found in [13].

Here, we provide a theoretical analysis of VIF under a Generalized GSM (GGSM)
model of NSS [14]. We achieve this by deriving novel results on the
Multivariate Generalized Gaussian Distribution (MGGD). Further, we also derive
results on the multivariate kurtosis, which we use to approximate the true
value of VIF, assuming the images being analyzed follow a GGSM model.

##  II Visual Information Fidelity

We begin by reviewing VIF, which is an information-theoretic FR IQA model. Let
the coefficients obtained after a bandpass transform on the reference image,
such as a steerable pyramid decomposition [15], be
$\mathbf{C_{i}}\in\mathbb{R}^{M}$, $i=1\dots N$. Assuming a signal attenuation
and additive noise distortion model, where
$\mathbf{D_{i}}=g_{i}\mathbf{C_{i}}+\mathbf{V_{i}}$

are distorted bandpass coefficients, $g_{i}$ are deterministic scalars, and
$\mathbf{V_{i}}\sim\mathcal{N}(0,\sigma_{v}^{2}\mathbf{I})$.

VIF assumes that the $\mathbf{C_{i}}$ are distributed as a Gaussian Scale
Mixture (GSM). That is, each $\mathbf{C_{i}}$ can be expressed as
$\mathbf{C_{i}}=Z_{i}\mathbf{U_{i}}$, where $Z_{i}$ are non-negative scalar
random variables (RVs) and $\mathbf{U_{i}}\sim\mathcal{N}(0,\mathbf{C_{U}})$,
$\mathbf{C_{U}}\in\mathbb{R}^{M\times M}$, independent of $Z_{i}$.

In VIF, uncertainty of perception, such as neural noise is modelled as
additive white Gaussian noise (AWGN), leading to observed reference and
distorted bandpass coefficients $\mathbf{E_{i}}$ and $\mathbf{F_{i}}$, given
by

$\mathbf{E_{i}}=Z_{i}\mathbf{U_{i}}+\mathbf{N_{i}},$ (1)
$\mathbf{F_{i}}=g_{i}Z_{i}\mathbf{U_{i}}+\mathbf{N^{\prime}_{i}},$ (2)

where
$\mathbf{N_{i}},\mathbf{N^{\prime}_{i}}\sim\mathcal{N}(0,\sigma_{n}^{2}\mathbf{I})$
are independent and identically distributed (iid). Then, assuming a
$K$-subband decomposition, the VIF index is defined as

$VIF=\frac{\sum\limits_{k=1}^{K}\sum\limits_{i=1}^{N}I(\mathbf{C_{i}^{k}};\mathbf{F_{i}^{k}}|Z_{i}^{k}=z_{i}^{k})}{\sum\limits_{k=1}^{K}\sum\limits_{i=1}^{N}I(\mathbf{C_{i}^{k}};\mathbf{E_{i}^{k}}|Z_{i}^{k}=z_{i}^{k})}.$
(3)

##  III Multivariate Generalized Gaussian Distribution

The probability density function (pdf) $f_{\mathbf{U}}$ of an $M$-dimensional
MGGD RV $\mathbf{U}\sim\mathcal{MGGD}(0,\alpha,\mathbf{C_{U}})$ is given by

$f_{U}(\mathbf{u})=\frac{\alpha\Gamma(\frac{M}{2})\det(\mathbf{C_{U}})^{-\frac{1}{2}}}{\pi^{\frac{M}{2}}\Gamma(\frac{M}{2\alpha})2^{\frac{M}{2\alpha}}}\exp\left(-\frac{1}{2}\left(\mathbf{u}^{T}\mathbf{C_{U}}^{-1}\mathbf{u}\right)^{\alpha}\right).$
(4)

If $\mathbf{U_{i}}\sim\mathcal{MGGD}(0,\alpha,\mathbf{C_{U}})$ in (1) - (2),
where $\alpha>0$, the resulting model is called a Generalized GSM (GGSM). When
$\alpha=1$, the GGSM and GSM models are identical.

The first result we derive in this section is an expression for the
differential entropy of an MGGD, which is hitherto known only for the scalar
case ($M=1$).

######  Lemma 1.

The differential entropy $h(\mathbf{U})$ of a RV
$\mathbf{U}\sim\mathcal{MGGD}(0,\alpha,\mathbf{C_{U}})$ is given by

$h(\mathbf{U})=\frac{M}{2\alpha}-\log\frac{\alpha\Gamma(\frac{M}{2})}{\pi^{\frac{M}{2}}\Gamma(\frac{M}{2\alpha})2^{\frac{M}{2\alpha}}}+\frac{1}{2}\log\det(\mathbf{C_{U}}).$
(5)

We will omit the proof of Lemma 1, since the proof of Lemma 2 contains all the
techniques used here.

The second result we derive is the expression for the Fisher information
Matrix (FIM), under translation, of an MGGD RV.

######  Definition 1.

If $\mathbf{X}$ has a probability density function $f_{X}$, its Fisher
information matrix is defined as

$\mathbf{J(X)}=E\left[\left(\frac{\nabla
f_{X}(\mathbf{X})}{f_{X}(\mathbf{X})}\right)\left(\frac{\nabla
f_{X}(\mathbf{X})}{f_{X}(\mathbf{X})}\right)^{T}\right],$

where $\nabla$ denotes the gradient operator.

######  Lemma 2.

The FIM $\mathbf{J(U)}$ of an MGGD is finite iff
$\alpha>\frac{1}{2}-\frac{M}{4}$, and is given by

$\mathbf{J(U)}=\frac{2^{2-\frac{1}{\alpha}}\alpha^{2}}{M\Gamma(\frac{M}{2\alpha})}\Gamma\left(2+\frac{M-2}{2\alpha}\right)\mathbf{C_{U}}^{-1}.$
(6)

This quantity is easy to derive for an MGGD when $M=1$, but no such an
expression for the multivariate case is available in the literature. This
result will be used to apply Lemma 3 in Section IV.

###### Proof.

Recall that under an invertible linear transformation
$\mathbf{Y}=\mathbf{AX}$, we know that $\mathbf{Y}$ has the pdf

$f_{Y}(\mathbf{y})=\frac{1}{\det{\mathbf{A}}}f_{X}(\mathbf{A}^{-1}\mathbf{y}).$
(7)

Differentiating with respect to $\mathbf{y}$, and applying the chain rule for
gradients, it follows that the FIM of $\mathbf{Y}$ is given by

$\mathbf{J(Y)}=(\mathbf{A}^{-1})^{T}\mathbf{J(X)}\mathbf{A}^{-1}.$ (8)

Now, consider transformation
$\mathbf{Y}=\mathbf{C_{U}}^{-\frac{1}{2}}\mathbf{U}$. Then, using (7), the
$(p,q)$ th term of $\mathbf{J(Y)}$ may be expressed

$[\mathbf{J(Y)}]_{pq}=\alpha^{2}\int\limits_{\mathbb{R}^{d}}y_{p}y_{q}\left\lVert\mathbf{y}\right\rVert_{2}^{4(\alpha-1)}f_{Y}(\mathbf{y})d\mathbf{y}.$
(9)

But, when $p\neq q$, (9) is the integral of an odd function, hence,

$[\mathbf{J(Y)}]_{pq}=0,\quad\text{if }p\neq q.$ (10)

When $p=q$, the integral may be evaluated using $M$-spherical coordinates,
which we omit the details of for brevity. Then, using (8) yields the result in
(6).

Note that the condition on $\alpha$ is necessary and sufficient for the Gamma
function in the numerator of (6) to be defined. ∎

##  IV Bounds on VIF under the GGSM Model

###  IV-A Bounds on Differential Entropy

Next, we derive bounds on VIF under the GGSM model. For simplicity, consider
one subband of the bandpass decomposition, yielding random fields
$\mathbf{C}$, $\mathbf{D}$, $\mathbf{E}$, and $\mathbf{F}$ defined as in
Section II.

To derive bounds on VIF, we will use bounds on the differential entropies of
sums of RVs. Let $\mathbf{X}$ and $\mathbf{Y}$ be two independent RVs. Then,
using Shannon’s Entropy Power Inequality (EPI) [16], we obtain.

$h(\mathbf{X}+\mathbf{Y})\geq\frac{M}{2}\log\left(\exp\left(\frac{2}{M}h(\mathbf{X})\right)+\exp\left(\frac{2}{M}h(\mathbf{Y})\right)\right).$
(11)

To obtain an upper bound, we use the following result [17].

######  Lemma 3.

Let
$\mathcal{F}=\\{\mathbf{X}\in\mathbb{R}^{M}|E[\log(1+\left\lVert\mathbf{X}\right\rVert_{2})]<\infty\\}$,
i.e., the set of RVs having finite logarithmic moments. If
$\mathbf{X}\in\mathcal{F}$ has finite $h(\mathbf{X})$ and $\mathbf{J(X)}$, and
$\mathbf{Z}\sim\mathcal{N}(\mu,\sigma^{2}\mathbf{I})$ is independent of
$\mathbf{X}$, then

$h(\mathbf{X}+\mathbf{Z})\leq
h(\mathbf{X})+\frac{M}{2}\log\left(1+\frac{\sigma^{2}}{M}\operatorname{tr}\mathbf{J(X)}\right).$

Since MGGD is a bounded pdf with a finite covariance, its logarithmic moment
exists. Hence, $\mathbf{U}\in\mathcal{F}$. From Section III, we know that
$h(\mathbf{U})$ is always finite and $\mathbf{J(U)}$ has finite entries iff
$\alpha>\frac{1}{2}-\frac{M}{4}$. Under this mild constraint, we can apply
Lemma 3 to the entropy terms in VIF.

###  IV-B Bounds on VIF

In (3), the mutual information term corresponding to the reference image is

$I(\mathbf{C_{i}};\mathbf{E_{i}}|z_{i})=h(\mathbf{C_{i}}+\mathbf{N_{i}}|z_{i})-h(\mathbf{N_{i}}|z_{i}).$
(12)

Since $\mathbf{C_{i}}=Z_{i}U$, then
$\mathbf{C_{i}}|_{Z_{i}=z_{i}}\sim\mathcal{MGGD}(0,z_{i}^{2}\mathbf{C_{U}},\alpha)$.
Using this relationship and the expression for the differential entropy of a
Gaussian RV, we find the following lower and upper bounds on (12),
respectively:

$I_{l}(\mathbf{C_{i}};\mathbf{E_{i}}|z_{i})\leq
I(\mathbf{C_{i}};\mathbf{E_{i}}|z_{i})\leq
I_{u}(\mathbf{C_{i}};\mathbf{E_{i}}|z_{i})$ (13)

where

$I_{l}(\mathbf{C_{i}};\mathbf{E_{i}}|z_{i})=\frac{M}{2}\log\left(1+\frac{z_{i}^{2}}{2\pi
e\sigma_{n}^{2}}\exp\left(\frac{2h(\mathbf{U})}{M}\right)\right),$ (14)

and

$I_{u}(\mathbf{C_{i}};\mathbf{E_{i}}|z_{i})=h(\mathbf{U})+\frac{M}{2}\log\left(\frac{z_{i}^{2}}{2\pi
e\sigma_{n}^{2}}+\frac{\operatorname{tr}\mathbf{J(U)}}{2\pi eM}\right).$ (15)

Similarly, the mutual information term in 3 corresponding to the distorted
image is

$I(\mathbf{C_{i}};\mathbf{F_{i}}|z_{i})=\sum\limits_{i=1}^{N}h(g_{i}\mathbf{C_{i}}+\mathbf{V_{i}}+\mathbf{N^{\prime}_{i}}|z_{i})-h(\mathbf{V_{i}}+\mathbf{N^{\prime}_{i}}|z_{i}).$
(16)

As above,
$g_{i}\mathbf{C_{i}}|_{Z_{i}=z_{i}}\sim\mathcal{MGGD}(0,g_{i}^{2}z_{i}^{2}\mathbf{C_{U}},\alpha)$.
Using this relationship and the expression for the differential entropy of a
Gaussian RV, we derive the following lower and upper bounds on (16),
respectively:

$I_{l}(\mathbf{C_{i}};\mathbf{F_{i}}|z_{i})\leq
I(\mathbf{C_{i}};\mathbf{F_{i}}|z_{i})\leq
I_{u}(\mathbf{C_{i}};\mathbf{F_{i}}|z_{i})$ (17)

where

$\begin{split}I_{l}(\mathbf{C_{i}};\mathbf{F_{i}}|z_{i})=\frac{M}{2}&\log\left(1+\frac{(g_{i}z_{i})^{2}}{2\pi
e\left(\sigma_{n}^{2}+\sigma_{v}^{2}\right)}\exp\left(\frac{2h(\mathbf{U})}{M}\right)\right),\end{split}$
(18)

and

$I_{u}(\mathbf{C_{i}};\mathbf{F_{i}}|z_{i})=h(\mathbf{U})+\frac{M}{2}\log\left(\frac{(g_{i}z_{i})^{2}}{2\pi
e(\sigma_{v}^{2}+\sigma_{n}^{2})}+\frac{\operatorname{tr}\mathbf{J(U)}}{2\pi
eM}\right).$ (19)

Repeating this analysis for all $K$ subbands, we obtain the following bound on
VIF under the GGSM model.

$\begin{split}\frac{\sum\limits_{k=1}^{K}\sum\limits_{i=1}^{N}I_{l}(\mathbf{C_{i}^{k}};\mathbf{F_{i}^{k}}|z_{i}^{k})}{\sum\limits_{k=1}^{K}\sum\limits_{i=1}^{N}I_{u}(\mathbf{C_{i}^{k}};\mathbf{E_{i}^{k}}|z_{i}^{k})}&\leq\text{VIF}\leq\frac{\sum\limits_{k=1}^{K}\sum\limits_{i=1}^{N}I_{u}(\mathbf{C_{i}^{k}};\mathbf{F_{i}^{k}}|z_{i}^{k})}{\sum\limits_{k=1}^{K}\sum\limits_{i=1}^{N}I_{l}(\mathbf{C_{i}^{k}};\mathbf{E_{i}^{k}}|z_{i}^{k})}.\end{split}$
(20)

##  V Approximating VIF In Practical Applications

When applying this model to practical applications, like quality assessment of
UGC pictures, finding bounds on the values of VIF supplies insufficient
accuracy. However, obtaining an expression for VIF assuming a more general
framework is difficult, since the distribution of the sum of MGGD and Gaussian
RVs is not easily characterized, even in the scalar case [18] [19].

To overcome this limitation, [18] and [19] approximated the distribution of
the sum of independent scalar GGD RVs by a GGD RV, using a moment-matching
method. Specifically, if $U,V$ are independent scalar GGD RVs, and $W=U+V$,
then $W$ is approximated as a GGD RV $\tilde{W}$, such that $W$ and
$\tilde{W}$ share the same means, variances and kurtoses. In the multivariate
case, MGGDs are modelled using the mean, covariance matrix and Mardia’s
multivariate kurtosis [20] [21].

######  Definition 2 (Mardia’s Kurtosis).

Given an $M$-dimensional RV $\mathbf{X}$ having mean $\boldsymbol{\mu}$ and
covariance $\boldsymbol{\Sigma}$, Mardia’s Kurtosis of $X$ is given by

$\gamma_{2}(\mathbf{X})=E[\left((\mathbf{X}-\boldsymbol{\mu})^{T}\boldsymbol{\Sigma}^{-1}(\mathbf{X}-\boldsymbol{\mu})\right)^{2}]-M(M+2).$

The work in [18] and [19] is restricted to the scalar case because an
expression for Mardia’s kurtosis of sums of RVs is hitherto unknown. We bridge
this gap by deriving the following results. In the following lemma, let
$\mathbb{S}^{M}_{++}$ denote the set of $M\times M$ symmetric positive
definite matrices.

######  Lemma 4.

Let $\mathbf{X}$ and $\mathbf{Y}$ be independent $M$-dimensional RVs having
zero means, without loss of generality, and covariance matrices
$\boldsymbol{\Sigma}_{\mathbf{X}},\boldsymbol{\Sigma}_{\mathbf{Y}}$
respectively. Let $\mathbf{Z}=\mathbf{X}+\mathbf{Y}$, having covariance
$\boldsymbol{\Sigma_{Z}}=\boldsymbol{\Sigma_{X}}+\boldsymbol{\Sigma_{Y}}$.
Further, for any RVs $\mathbf{A},\mathbf{B}$ having covariances
$\boldsymbol{\Sigma_{A}},\boldsymbol{\Sigma_{B}}\in\mathbb{S}^{M}_{++}$, let
$\boldsymbol{\Delta_{AB}}=\boldsymbol{\Sigma_{A}}^{\frac{1}{2}}\boldsymbol{\Sigma_{B}}^{-1}\boldsymbol{\Sigma_{A}}^{\frac{1}{2}}$,
and
$\rho_{\mathbf{A}\mathbf{B}}=2\left\lVert\boldsymbol{\Delta_{AB}}\right\rVert_{F}^{2}+\operatorname{tr}\left(\mathbf{\Delta_{AB}}\right)^{2}.$
Then, Mardia’s Kurtosis of $\mathbf{Z}$ sastisfies the relationship

$\begin{split}\gamma_{2}(\mathbf{Z})=&E\left[\left(\mathbf{X}^{T}\boldsymbol{\Sigma_{Z}}^{-1}\mathbf{X}\right)^{2}\right]-\rho_{\mathbf{X},\mathbf{Z}}+\\\
&E\left[\left(\mathbf{Y}^{T}\boldsymbol{\Sigma_{Z}}^{-1}\mathbf{Y}\right)^{2}\right]-\rho_{\mathbf{Y},\mathbf{Z}}.\end{split}$
(21)

The proof of Lemma 4 follows from the definitions and properties of
multivariate moments and cumulants defined in [22], and evaluating
expectations of Gaussian RVs. We omit the details of this proof for brevity.

######  Definition 3.

A RV $\mathbf{X}$ is said to Elliptically Distributed with mean
$\boldsymbol{\mu}$, covariance $\boldsymbol{\Sigma}$ and a “generating
function” $g$, i.e.,
$\mathbf{X}\sim\mathcal{ED}(\boldsymbol{\mu},\boldsymbol{\Sigma},g)$ if the
pdf of $\mathbf{X}$ is of the form

$f_{\mathbf{X}}(\mathbf{x})=C\det(\Sigma)^{\frac{1}{2}}g((\mathbf{x}-\boldsymbol{\mu})^{T}\boldsymbol{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu})).$

######  Theorem 1.

Let
$\mathbf{X}\sim\mathcal{ED}\left(0,\boldsymbol{\Sigma_{X}},g_{\mathbf{X}}\right)$,
$\mathbf{Y}\sim\mathcal{ED}\left(0,\boldsymbol{\Sigma_{Y}},g_{\mathbf{Y}}\right)$
be independent $M$-dimensional RVs, and $\mathbf{Z}=\mathbf{X}+\mathbf{Y}$.
For a generating function $g$, let
$\lambda_{g}^{(M)}=E[V_{1}^{4}]/E[V_{1}^{2}V_{2}^{2}]-3$, where $V_{i}$ are
components of the $M$-dimensional RV $\mathbf{V}\sim\mathcal{ED}(0,I,g)$, and
let ∘2 denote elementwise squaring. Then,

$\displaystyle\gamma_{2}(\mathbf{Z})=$
$\displaystyle\frac{\left(\rho_{\mathbf{X}\mathbf{Z}}+\lambda_{g_{\mathbf{X}}}^{(M)}\operatorname{tr}\left(\boldsymbol{\Delta}_{\mathbf{X}\mathbf{Z}}^{\circ
2}\right)\right)\left(\gamma_{2}\left(\mathbf{X}\right)+M\left(M+2\right)\right)}{M\left(M+2+\lambda_{g_{\mathbf{X}}}^{(M)}\right)}+$
$\displaystyle\frac{\left(\rho_{\mathbf{Y}\mathbf{Z}}+\lambda_{g_{\mathbf{Y}}}^{(M)}\operatorname{tr}\left(\boldsymbol{\Delta}_{\mathbf{Y}\mathbf{Z}}^{\circ
2}\right)\right)\left(\gamma_{2}\left(\mathbf{Y}\right)+M\left(M+2\right)\right)}{M\left(M+2+\lambda_{g_{\mathbf{Y}}}^{(M)}\right)}-$
$\displaystyle\rho_{\mathbf{X}\mathbf{Z}}-\rho_{\mathbf{Y}\mathbf{Z}}.$ (22)

###### Proof.

Let $\mathbf{V}\sim\mathcal{ED}(0,I,g)$ be an $M$-dimensional RV. Observing
that the distribution of $\mathbf{V}$ is invariant to permutations of its
components, it follows that $E[Y_{i}^{4}]=E[Y_{1}^{4}]$, and
$E[Y_{i}^{2}Y_{j}^{2}]=E[Y_{1}^{2}Y_{2}^{2}]$ for all $i\neq j$. Using this
symmetry and the definition of $\lambda_{g_{\mathbf{X}}}^{(M)}$, it follows
that

$\displaystyle\gamma_{2}(\mathbf{V})$
$\displaystyle=E[\left\lVert\mathbf{V}\right\rVert^{4}]-M(M+2)$
$\displaystyle=ME[V_{1}^{4}]+M(M-1)E[V_{1}^{2}V_{2}^{2}]-M(M+2)$
$\displaystyle=ME[V_{1}^{2}V_{2}^{2}](\lambda_{g}^{(M)}+M+2)-M(M+2)$ (23)

Consider the term
$E\left[\left(\mathbf{X}^{T}\boldsymbol{\Sigma}_{\mathbf{Z}}^{-1}\mathbf{X}\right)^{2}\right]$.
Using the transformation
$\mathbf{\tilde{X}}=\boldsymbol{\Sigma_{X}}^{-\frac{1}{2}}\mathbf{X}$ and the
definition of $\boldsymbol{\Delta_{AB}}$, it follows that

$E\left[\left(\mathbf{X}^{T}\boldsymbol{\Sigma}_{\mathbf{Z}}^{-1}\mathbf{X}\right)^{2}\right]=E\left[\left(\mathbf{\tilde{X}}^{T}\boldsymbol{\Delta_{XZ}}\mathbf{\tilde{X}}^{T}\right)^{2}\right].$
(24)

Note that $\mathbf{\tilde{X}}\sim\mathcal{ED}(0,I,g_{\mathbf{X}})$, and since
kurtosis is invariant to linear transformations,
$\gamma_{2}(\mathbf{\tilde{X}})=\gamma_{2}(\mathbf{X})$. Expanding the
polynomial in (24), we observe that the expectation of terms having odd powers
is 0. Further, by the same symmetry argument used in (V), it follows that

$\displaystyle
E\left[\left(\mathbf{X}^{T}\boldsymbol{\Sigma}_{\mathbf{Z}}^{-1}\mathbf{X}\right)^{2}\right]=$
$\displaystyle E[\tilde{X}_{1}^{2}\tilde{X}_{2}^{2}]\sum_{i\neq
j}\left(2\left(\Delta_{\mathbf{XZ}}^{ij}\right)^{2}+\Delta_{\mathbf{XZ}}^{ii}\Delta_{\mathbf{XZ}}^{jj}\right)$
$\displaystyle+E[\tilde{X}_{1}^{4}]\sum_{i}\left(\Delta_{\mathbf{XZ}}^{ii}\right)^{2}.$
(25)

Using (V), it follows that (V) evaluates to the first term in (1). Repeating
the analysis for $\mathbf{Y}$ completes the proof. ∎

For an MGGD, which is an elliptical distribution, it can be shown that
$\lambda_{g}^{(M)}=0$. In addition, it can be shown that Mardia’s Kurtosis of
Gaussian RVs is zero. Using these properties, we obtain the following
corollary of Theorem 1.

######  Corollary 1.1.

Let $\mathbf{X}$ be an MGGD RV having zero mean, shape parameter $\alpha$ and
covariance matrix $\boldsymbol{\Sigma}$,
$\mathbf{Y}\sim\mathcal{N}(0,\sigma^{2}\mathbf{I})$ independent of
$\mathbf{X}$ and $\mathbf{Z}=\mathbf{X}+\mathbf{Y}$. Then,

$\gamma_{2}(\mathbf{Z})=\gamma_{2}(\mathbf{X})\frac{\rho_{\mathbf{XZ}}}{M(M+2)}.$

Using Corollary 1.1, we can utilize moment-matching to approximate the
distributions of $\mathbf{E_{i}}$ and $\mathbf{F_{i}}$ by MGGD RVs. As a
result, we are able to approximate the mutual information terms in (12) and
(16), using the expression for differential entropy of MGGD from Lemma 1.

To illustrate this, consider $I(\mathbf{C_{i}};\mathbf{E_{i}}|z_{i})$ in (12).
Note that $h(\mathbf{N_{i}})$ is known, since $\mathbf{N_{i}}$ is a Gaussian
RV. Let $\boldsymbol{\Sigma_{U}}$ be the covariance matrix of $\mathbf{U}$.
Then, given $Z_{i}=z_{i}$, the covariance of $\mathbf{C_{i}}$ is
$z_{i}^{2}\boldsymbol{\Sigma_{U}}$. Using Corollary 1.1, the distribution of
$\mathbf{E_{i}}=\mathbf{C_{i}}+\mathbf{N_{i}}$ is approximated as
$\mathbf{\tilde{E_{i}}}\sim\mathcal{MGGD}(0,\mathbf{C_{\tilde{E_{i}}}},\beta_{i})$
by moment-matching, such that

$\text{Cov}(\mathbf{\tilde{E_{i}}})=\text{Cov}(\mathbf{E_{i}})=z_{i}^{2}\boldsymbol{\Sigma_{U}}+\sigma_{n}^{2}\mathbf{I},$
(26)
$\gamma_{2}(\mathbf{\tilde{E_{i}}})=\gamma_{2}(\mathbf{E_{i}})=\gamma_{2}(\mathbf{U})\frac{\rho_{\mathbf{C_{i}}\mathbf{E_{i}}}}{M(M+2)}.$
(27)

Knowing the covariance and kurtosis, the parameters
$\mathbf{C_{\tilde{E_{i}}}}$ and $\beta_{i}$ may be estimated, following the
standard procedure for MGGDs [21]. Using these estimates in (5), we
approximate $h(\mathbf{C_{i}}+\mathbf{N_{i}})\approx
h(\tilde{\mathbf{E_{i}}})$,

using which $I(\mathbf{C_{i}};\mathbf{E_{i}}|z_{i})$ may be approximated. A
similar method may be used to approximate
$I(\mathbf{C_{i}};\mathbf{F_{i}}|z_{i})$ in (16). Approximate values of mutual
information obtained in this manner may be used in (3) to approximate VIF.

##  VI Conclusion and Discussion

In this work, we derived novel expressions for properties of MGGD RVs, using
which we analyzed the VIF model under a GGSM natural image model. In addition,
we proposed a method to approximate VIF using a novel result regarding
Mardia’s kurtosis of the sum of independent RVs. This framework extends easily
to other information-theoretic quality assessment models, such as ST-RRED and
SpEED-QA. To the best of our knowledge, this work is the first theoretical
extension of an IQA model to a relaxed model of NSS.

The only obstacle in the way of testing this work is the lack of a
multivariate GGSM modeling algorithm. The authors of [23] proposed iterative
maximum likelihood and MCMC-based Bayesian methods to estimate the parameters
of a scalar GGSM, under some additional assumptions on the mixing distribution
$Z_{i}$. An iterative modeling algorithm was proposed in [14], but after one
iteration of that method, the divisively normalized vectors lie on the surface
of an ellipsoid. As a result, no proper continuous distribution can model
these vectors, making this method unstable and inaccurate in practice.
Developing a sound GGSM modeling algorithm will pave the way for quality
assessment models that are resilient to deviations from ideal NSS.

## References

  * [1] Daniel L Ruderman,  “The statistics of natural images,”  Network: Computation in Neural Systems, vol. 5, no. 4, pp. 517–548, 1994.
  * [2] M.J. Wainwright and E.P. Simoncelli,  “Scale mixtures of gaussians and the statistics of natural images,”  in International Conference on Neural Information Processing Systems. 1999, p. 855–861, MIT Press.
  * [3] H.R. Sheikh, A.C. Bovik, and G. De Veciana,  “An information fidelity criterion for image quality assessment using natural scene statistics,”  IEEE Transactions on Image Processing, vol. 14, no. 12, pp. 2117–2128, 2005.
  * [4] H.R. Sheikh and A.C. Bovik,  “Image information and visual quality,”  IEEE Transactions on Image Processing, vol. 15, no. 2, pp. 430–444, 2006.
  * [5] Z. Li, A. Aaron, I. Katsavounidis, A.K. Moorthy, and M. Manohara,  “Toward a practical perceptual video quality metric,”  The Netflix Tech Blog, vol. 6, 2016.
  * [6] R. Soundararajan and A. C. Bovik,  “Video quality assessment by reduced reference spatio-temporal entropic differencing,”  IEEE Transactions on Circuits and Systems for Video Technology, vol. 23, no. 4, pp. 684–694, 2013.
  * [7] C. G. Bampis, P. Gupta, R. Soundararajan, and A. C. Bovik,  “SpEED-QA: Spatial efficient entropic differencing for image and video quality,”  IEEE Signal Processing Letters, vol. 24, no. 9, pp. 1333–1337, 2017.
  * [8] A. K. Moorthy and A. C. Bovik,  “Blind image quality assessment: From natural scene statistics to perceptual quality,”  IEEE Transactions on Image Processing, vol. 20, no. 12, pp. 3350–3364, 2011.
  * [9] A. Mittal, A. K. Moorthy, and A. C. Bovik,  “No-reference image quality assessment in the spatial domain,”  IEEE Transactions on Image Processing, vol. 21, no. 12, pp. 4695–4708, 2012.
  * [10] A. Mittal, R. Soundararajan, and A. C. Bovik,  “Making a “completely blind” image quality analyzer,”  IEEE Signal Processing Letters, vol. 20, no. 3, pp. 209–212, 2013.
  * [11] Z. Sinno and A. C. Bovik,  “Large-scale study of perceptual video quality,”  IEEE Transactions on Image Processing, vol. 28, no. 2, pp. 612–627, 2019.
  * [12] Y. Li, S. Meng, X. Zhang, S. Wang, Y. Wang, and S. Ma,  “Ugc-video: Perceptual quality assessment of user-generated videos,”  in IEEE Conference on Multimedia Information Processing and Retrieval, 2020, pp. 35–38.
  * [13] Z. Tu, Y. Wang, N. Birkbeck, B. Adsumilli, and A.C. Bovik,  “UGC-VQA: Benchmarking blind video quality assessment for user generated content,”  arXiv preprint arXiv:2005.14354, 2020.
  * [14] P. Gupta, A.K. Moorthy, R. Soundararajan, and A.C. Bovik,  “Generalized Gaussian scale mixtures: A model for wavelet coefficients of natural images,”  Signal Processing: Image Communication, vol. 66, pp. 87–94, 2018.
  * [15] E.P. Simoncelli and W.T. Freeman,  “The steerable pyramid: A flexible architecture for multi-scale derivative computation,”  in International Conference on Image Processing. IEEE, 1995, vol. 3, pp. 444–447.
  * [16] C.E. Shannon,  “A mathematical theory of communication,”  The Bell System Technical Journal, vol. 27, no. 3, pp. 379–423, 1948.
  * [17] J. Fahs and I. Abou-Faycal,  “A new tight upper bound on the entropy of sums,”  Entropy, vol. 17, no. 12, pp. 8312–8324, 2015.
  * [18] Q. Zhao, H. Li, and Y. Shen,  “On the sum of generalized gaussian random signals,”  in International Conference on Signal Processing, 2004, vol. 1, pp. 50–53 vol.1.
  * [19] H. Soury and M. Alouini,  “New results on the sum of two generalized gaussian random variables,”  in IEEE Global Conference on Signal and Information Processing, 2015, pp. 1017–1021.
  * [20] K. V. Mardia,  “Measures of multivariate skewness and kurtosis with applications,”  Biometrika, vol. 57, no. 3, pp. 519–530, 1970.
  * [21] E. Gómez, M.A. Gomez-Viilegas, and J.M. Marín,  “A multivariate generalization of the power exponential family of distributions,”  Communications in Statistics - Theory and Methods, vol. 27, no. 3, pp. 589–600, 1998.
  * [22] T. Kollo,  “Multivariate skewness and kurtosis measures with an application in ica,”  Journal of Multivariate Analysis, vol. 99, no. 10, pp. 2328–2338, 2008.
  * [23] Z. Rahnamaei, N. Nematollahi, and R. Farnoosh,  “The location-scale mixture exponential power distribution: A Bayesian and maximum likelihood approach,”  Journal of Applied Mathematics, vol. 2012, pp. 978282, Oct 2012.


# The waltz of tiny droplets and the flow they live in

S. Ravichandran<EMAIL_ADDRESS>Nordic Institute for Theoretical
Physics, KTH Royal Institute of Technology and Stockholm University, Stockholm
SE 10691  Interdisciplinary Programme in Climate Studies, Indian Institute of
Technology Bombay, Mumbai 400076  Rama Govindarajan<EMAIL_ADDRESS>International Centre for Theoretical Sciences, Tata Institute of Fundamental
Research, Bengaluru 560089, India.

###### Abstract

This article describes the dynamics of small inertial particles centrifuging
out of a single vortex. It shows the importance of caustics formation in the
vicinity of a single vortex: both for particle collisions and void formation.
From these single-vortex studies we provide estimates of the role of caustics
in high Reynolds number turbulence, and in the case of clouds, estimate how
they may help in rain initiation by bridging the droplet-growth bottleneck. We
briefly describe how the Basset-Boussinesq history force may be calculated by
a method which does not involve huge memory costs, and provide arguments for
its possible importance for droplets in turbulence. We discuss how phase
change could render cloud turbulence fundamentally different from turbulence
in other situations.

##  I Introduction

###  I.1 The flow situations

The dynamics of finite-sized particles or droplets in turbulent flow is
relevant in a wide range of natural and industrial settings. In the dispersal
of pollutants due to a volcanic eruption or from industrial flue gases, in a
sandstorm, or a cloud, or when a river carrying sediment disgorges itself into
the ocean, we encounter such interaction. We are interested here in spherical
particles, which are so small that to a very good approximation they are in
Stokes flow relative to the surrounding fluid.

Why is the dynamics, in flow, of individual Stokesian particles interesting?
First, such dynamics is often not passive, in that the particle could be
inertial, with finite acceleration relative to the flow. Such particles
therefore may not follow fluid streamlines, and this is responsible for a
variety of complex behaviour. In background turbulent flow, particle inertia
leads to clustering [1], caustics formation [2, 3], and, in some cases,
coalescence and growth [4, 5, 6, 7, 8] (the general features of such flows are
reviewed in [9]; see also the reviews by [10, 11] of properties relevant to
clouds). Further, the particles mentioned are often of density different from
the fluid, and sediment or rise under gravity [12, 13, 14, 15]. Notice that we
use the term ‘particles’ to stand for droplets as well, except when growth by
coalescence or phase change is relevant, when we will refer to droplets. This
interchangeable notation is applicable when droplets are sufficiently small
for surface tension to ensure that they remain spherical in shape. When
particles are not in dilute suspension, their interactions lead to new
physics, including instabilities, see e.g. [16, 17]. But here we discuss
dilute suspensions, where inter-particle interactions can be ignored.

We are particularly interested in droplets in atmospheric clouds. Clouds are a
crucial part of the Earth’s energy balance. Current estimates suggest that
clouds are responsible for a significant fraction of both the albedo (i.e.,
reflecting away of incoming solar ‘short-wave’ radiation) of the planet, as
well as the greenhouse effect (i.e., trapping of outgoing ‘long-wave’
radiation) [see, e.g. 18]. The respective magnitudes of these contributions
from different patches of cloud depend on their altitude, extent and
thickness, as well as the droplet size distribution within. The size of each
droplet in a cloud is a changing function of time, since droplets can whittle
away by evaporation or grow by condensation, depending on whether their
immediate neighbourhood is subsaturated or supersaturated in water vapour. The
energy released during these transformations plays a crucial role in the
dynamics of clouds, and in processes such as turbulent entrainment and mixing.

The dynamics of inertial droplets are also thought to be responsible for the
rapid onset of rain in warm (i.e. ice-free) clouds ([1, 19, 2, 20]). The
millimetre-size droplets required for rain initiation cannot be formed either
by vapour diffusion which is too slow for droplets larger than $O(10)$ micron,
or by gravitational collisions-coalescence which is not sufficiently rapid for
droplets smaller than $O(50)$ micron. The range $O(10-50)$ micron is called
the ‘droplet growth bottleneck’, and is thought to be ‘bridged’ by the
dynamics of inertial droplets in turbulent flow. It is known, for instance
that the presence of a few [one-in-a-million, see e.g. 19] large droplets
leads to a significantly more rapid generation of rain-sized droplets. We
posit here that these large droplets may arise due to turbulence-induced
collisions of $10$ micron droplets due to the action of small but powerful
single vortices.

In the rest of this review, in §I.2 we introduce the basic equations governing
the dynamics of inertial droplets and then discuss the formation of caustics
due to particle inertia, neglecting higher order effects, in §II. We describe
the effects of a nominally higher-order correction, namely the Basset history
term, in §III, showing that its effects are negligible for raindrop growth due
to gravitational settling. We then discuss the effects of condensation of
water vapour onto droplets on the dynamics of clouds in §IV, including some of
our own contributions to the study of the interactions of droplet inertia,
turbulence and phase change.

###  I.2 The Maxey-Riley equation

The dynamics of a heavy sphere in unsteady Stokes flow is described by the
famous Maxey-Riley (MR) equation [21]

$\displaystyle\dot{\mathbf{y}}$ $\displaystyle=$
$\displaystyle{\mathbf{v}}(t)$ (1) $\displaystyle\dot{\mathbf{v}}$
$\displaystyle=$
$\displaystyle-\frac{1}{St}({\mathbf{v}}-{\mathbf{u}})-\frac{3}{\sqrt{2\beta
St}}\left[\frac{1}{\sqrt{\pi
t}}({\mathbf{v}}(0)-{\mathbf{u}}(0))+\frac{1}{\sqrt{\pi}}\int_{0}^{t}\frac{\dot{\mathbf{v}}(s)-\dot{\mathbf{u}}(s)}{\sqrt{t-s}}ds\right]+\frac{1}{Fr^{2}}\mathbf{e}_{g}$
(2)

where ${\mathbf{y}}$ and ${\mathbf{v}}$ respectively are the position and
velocity of the particle at time $t$, and the dots refer to differentiation in
time along the particle trajectory. We prescribe that the density ratio
$\beta\equiv\rho_{p}/\rho_{f}\gg 1$, where $\rho_{p}$ and $\rho_{f}$
respectively are the densities of the particle and the fluid. In this limit,
the Froude number $Fr^{-2}\equiv gT_{f}^{2}/L_{f}$, where $g$ is the
acceleration due to gravity and $T_{f}$ and $L_{f}$ are characteristic flow
time and length scales respectively. The Stokes number $St\equiv\tau/T_{f}$,
where

$\tau\equiv\frac{2\beta a^{2}}{9\nu},$ (3)

is the particle time scale, $a$ being the radius of the particle and $\nu$ the
kinematic viscosity of the fluid. Equation (2) is in general nonlinear, since
${\mathbf{u}}={\mathbf{u}}({\mathbf{y}}(t))$, the velocity of the fluid at the
particle location, can depend nonlinearly on ${\mathbf{y}}$. The Faxen
correction terms due to the curvature of the velocity field on the scale of
the particle have been neglected in this representation of the MR equation.
This is reasonable because Stokesian particles in turbulence tend to be far
smaller than the structures in the flow. Most often [22, 23, 3] in the absence
of gravity, this equation is further truncated, after the first term on the
right hand side, reducing it to the simplified Maxey-Riley (SMR) equation:

$\dot{\mathbf{v}}=-\frac{1}{St}({\mathbf{v}}-{\mathbf{u}}).$ (4)

The Basset-Boussinesq history term (‘Basset-history’ or ‘BBH’ for short),
i.e., the second term on the right hand side of equation (2), has been dropped
in this approximation. At first sight, this seems to be reasonable for
particles of small Stokes number, since we have retained all terms up to
$O(St)$ and neglected terms of $O(St^{3/2})$ and higher. We shall proceed with
this assumption in the next section and examine it thereafter.

We ask first how turbulence, in the absence of gravity, determines particle
dynamics. This question may be answered in some degree by obtaining, in
detail, the dynamics near a single Lamb-Oseen vortex, as we do in Section II
where we discuss the role and importance of single-vortex caustics. We then
investigate the relevance of our findings in a generic high Reynolds number
turbulence.

##  II Dynamics near a Lamb-Oseen vortex, and the caustics question

We begin this discussion by defining caustics as situations where two or more
particles arrive simultaneously at the same place with different velocities.
It is believed that caustics are important for enhanced particle collisions
[e.g. 2, 4, 24] which in turn are a crucial ingredient, in the case of
droplets, for growth by coalescence. Collisions are an important reason why we
have rain at all, and the background turbulence in clouds, with its associated
vorticity, is considered to be a lead player in this. A natural way to divide
a turbulent flow in the context of particle dynamics is into regions of
vorticity and regions of strain, since inertial particles are known to
centrifuge out of the former and collect in the latter [25, 26]. It is well
known that caustics form when particles centrifuging out of different vortices
collide with each other [2]. But are single-vortex caustics possible? This was
answered in the affirmative by [27] and is discussed in the following. In this
study we limit ourselves to particles which are far heavier than the fluid
they are embedded in, but the analysis can easily be modified to include a
finite density ratio, which in itself can make the dynamics more interesting.
In the following, when we use the term ‘caustics’, we are referring to single-
vortex caustics.

Patches of vorticity occur in a range of shapes and sizes in turbulence, but a
given patch of vorticity, if far away from other patches, will axisymmetrise
into a Lamb-Oseen vortex in whose cross-section vorticity varies as a
gaussian:

$\omega=\omega_{v}\exp\left[-\frac{r_{d}^{2}}{r_{v}^{2}}\right],$ (5)

where $\omega_{v}=\Gamma/(\pi r_{v}^{2})$ is the vorticity at the vortex
centre, $r_{d}$ is the dimensional radial distance from the centre of the
vortex and $r_{v}$ is the characteristic size of the vortex at a given time.
In fact, as was pointed out by an anonymous referee of another paper, a Lamb-
Oseen vortex can be termed the “drosophila” of turbulence, a building block
from whose behaviour we may glean general principles much as one does in
biology using that organism. In two-dimensional turbulence, it was verified
[28] that vortex patches tend to be circular and display a Gaussian cross-
section of vorticity [see also 29]. There is increasing evidence in high
Reynolds number simulations that intermittent vortices tend to be tubelike.

In two dimensions, in the neighbourhood of a Lamb-Oseen vortex, the SMR
equation in cylindrical-polar ($r-\theta$) coordinates reads [27]

$\displaystyle\ddot{r}+\dot{r}=\frac{\zeta^{2}}{r^{3}},$ (6)
$\displaystyle\dot{\zeta}+\zeta=1-e^{-St_{v}r^{2}}$ (7)

where we have nondimensionalised the equations with the special length scale

$L\equiv\left(\frac{\Gamma\tau}{2\pi}\right)^{1/2},$ (8)

the particle time scale $\tau$, and denoted $\zeta\equiv r^{2}\dot{\theta}$ as
the non-dimensional angular momentum of the particle. Note that $L$
encompasses properties of the flow and the particle. The ratio of $L$ and the
vortex size $r_{v}$ is also a ratio of time scales of particles and vortex,
and appears in the form of a vortex time-scale based Stokes number

$St_{v}\equiv\frac{\omega_{v}\tau}{2}=\frac{\Gamma\tau}{2\pi r_{v}^{2}}.$ (9)

Now, the Stokes number in turbulent flow is normally defined in terms of the
large time scale $T_{f}$ in the flow, which is much bigger than the time scale
of the small vortices. Particles could thus be of small Stokes number $St$ in
terms of the large scale and large Stokes number $St_{v}$ in terms of a single
small-scale vortex. We are interested in vortices in this regime.

###  II.1 Point vortex, inner solution: the formation of caustics

We note that for a point vortex, i.e., when $r_{v}=0$, equations (6) and (7)
become parameter free. The solution in this limit describes the dynamics of
Stokesian particles of any Stokes number near any point vortex irrespective of
its strength. Equation (7) is now easily solvable, and yields

$\zeta=1-(1-\zeta_{0})e^{-t}.$ (10)

Here and in the following, a subscript $0$ indicates the value of the quantity
at the initial time. Note that $\zeta\sim O(1)$, and for the analysis below we
will prescribe $\zeta_{0}\neq 0$. Equation (6) supports a boundary-layer
structure, and, for a point vortex, we may solve exactly for the dynamics
close to, and very far away from, the vortex centre. The inner solution, for
small distances from the vortex centre, is only relevant at short times, since
this region is empty of particles once they get centrifuged out. Following the
standard procedure for singular perturbation problems, we define inner
variables

$R_{i}\equiv\frac{r}{\delta_{i}}\quad{\rm and}\quad
T_{i}\equiv\frac{t}{\epsilon_{i}},$ (11)

which are $O(1)$, with $\delta_{i},\epsilon_{i}\ll 1$. Substituting these, and
equation (10) in equation (6), we find that the lowest order balance requires
$\epsilon_{i}\sim\delta_{i}^{2}$. Since these are arbitrary constants defined
by us, we may set $\epsilon_{i}=\delta_{i}^{2}$. At this order, there is no
further information on what exactly $\epsilon_{i}$ and $\delta_{i}$
individually are, but this information is not necessary either. Equation (6)
reduces at the lowest order to

$R_{i}^{\prime\prime}=\frac{\zeta_{0}^{2}}{R_{i}^{3}}$ (12)

where primes refer to differentiation with respect to $T_{i}$. Equation (12)
is an autonomous nonlinear ordinary differential equation, whose solution is

$R_{i}^{2}=\left[\frac{\zeta_{0}^{2}}{R_{i0}^{2}}+{R^{\prime}}_{i0}^{2}\right]T_{i}^{2}+2R_{i0}R^{\prime}_{i0}T_{i}+R_{i0}^{2}.$
(13)

Now caustics are formed when two sets of particles starting out at different
$R_{0}$ arrive at the same radial location $R_{c}$ at a particular time
$T_{c}$, with the subscript $c$ denoting caustics. The caustics time for a
given set of initial conditions may be obtained from equation (13). For
example, if two particles, initially at $R_{01}$ and $R_{02}$ have the same
$\zeta_{0}$ and start out with zero radial velocity, we have

$T_{c}=\frac{R_{01}R_{02}}{\zeta_{0}}.$ (14)

We arrive at the conclusion that the inner solution supports the formation of
caustics for a point vortex. Notably, all terms in equation (14) are $O(1)$,
which means that the caustics time $t_{c}=\epsilon_{i}T_{c}$ is very short.

###  II.2 Outer solution: no caustics

Figure 1: Fermat’s spiral describing the centrifugation of a heavy particle in
the outer region of a vortex. As time goes on, the radial distance covered per
cycle reduces and the time taken to cover each cycle increases.

At large time, all particles evacuate the vicinity of the vortex, and occupy
the region at large $r$. Moreover equation (10) gives $\zeta=1$. The outer
equations may be obtained in a similar manner as above by defining outer
variables, this time as $R_{o}=\delta_{o}r$ and $T_{o}=\epsilon_{o}t$ where
$\delta_{o},\epsilon_{o}\ll 1$ to get, at the lowest order,

$R_{o}^{\prime}=\frac{1}{R_{o}^{3}},$ (15)

the primes now referring to differentiation with respect to $T_{o}$. Note that
the subscript $o$ stands for the outer equation, and is not to be confused
with the subscript $0$ for initial conditions. The balance this time yields
$\epsilon_{o}=\delta_{o}^{4}$. Equation (15) may be solved immediately to
yield $R_{o}^{4}-R_{o0}^{4}=4T_{o}$. At large times $R_{o}\gg R_{o0}$, so we
may neglect the latter in comparison to the former. When put back in the
physical variables, we have $r^{4}=4t$, with $r^{2}\dot{\theta}=1$. In other
words $\theta/r^{2}=1/2$, i.e., a particle far away from the vortex will
execute a Fermat’s spiral, shown in schematic in figure 1, thus centrifuging
ever more slowly out of the vicinity.

It is evident that from equation (15) that a given radial distance away from
the vortex is associated with a unique particle radial velocity $\dot{R}_{o}$.
This means that particle dynamics may be described in the outer region in
terms of a field for the particle velocity, namely
$\dot{R}_{o}=\dot{R}_{o}(R_{o},\theta,T_{o})$. This in turn means there is no
possibility of caustics in this region.

###  II.3 The complete solution for single-vortex caustics

The singular perturbation study above is carried out just to be instructive.
The complete SMR for a particle near a Lamb-Oseen vortex, equations (6) and
(7), are trivial to solve numerically and we have done so. From the lowest-
order singular perturbation theory, we saw for a point vortex that caustics
happen at early times and $r\ll 1$ and never at late times and $r\gg 1$. We
may thus expect that the limit for caustics will happen at a caustics radius
$r_{c}\sim O(1)$. In fact, for a range of initial conditions [27] we find that
the caustics radius is given by

$r_{cd}^{2}\sim\Gamma\tau,\quad{\rm or}\quad r_{c}\sim 1,$ (16)

where the subscript $d$ denotes a dimensional quantity.

Figure 2:  The range of initial locations (the region between the red dashed
lines, where $r_{0}=r_{ci}$ and the blue dashed lines, where $r_{0}=r_{c}$)
broadens as the Stokes number $St_{v}$ (equation 9) increases. For each curve,
caustics occur in locations bounded by these vertical asymptotes. Thus, for
$St_{v}\lesssim 2.3$, no caustics occur. These curves were obtained using
inertial finite-time Lyapunov exponents (iFTLEs, see [30] in the radial
direction, the values of which are plotted as the ordinate. Figure 3:
Schematic of single-vortex caustics. The vortex is shown by the filled blue
circle. Caustics can happen only inside $r_{c}$, and an example of an inner
light green ring of particles just about to overtake the outer dark green ring
is shown. Particles within the caustics radius $r_{c}$ get centrifuged out
very rapidly, thus evacuating this region. Just beyond, an intensely clustered
ring of particles collects and continues centrifuging out at a far slower
rate. This is a singular perturbation problem with distinguished limits at
radii well within and well outside the caustics radius.

For a Lamb-Oseen vortex of finite size too, the complete solution of the SMR
equation shows caustics, see figure 2. The difference is that while for a
point vortex, particles starting from any $r<r_{c}$ participated in caustics,
for the Lamb-Oseen vortex we have caustics occurring for particles starting
within $r_{ci}\leq r\leq r_{c}$, where $r_{ci}$ is an inner limit for caustics
formation, which decreases as the Stokes number $St_{v}$ increases, as seen
from figure 2. Further, all caustics happen at times of $O(\tau)$, and the
region within $r_{c}$ is devoid of particles at later times, irrespective of
the value of $r_{ci}$. In other words, the region
$r_{d}\lesssim(\Gamma\tau)^{1/2}$ where caustics form is also the region which
gets voided extremely quickly (on the particle timescale) of particles due to
rapid centrifugation. This will be important for our calculations below of
void volumes. Another finding is that in the case of single-vortex caustics,
there is a direct and strong correlation between collisions and the formation
of caustics. Collisions are overwhelmingly the result of caustics formation,
and so caustics are crucial for droplet growth. Figure 2 indicates in
schematic that there is a sharp increase in the number density of particles
just outside $r_{c}$, which is one reason for enhanced collisions. This highly
clustered band then centrifuges out slowly on the Fermat spiral time scale.
Incidentally we do have rare occurrences of collisions without caustics, when
the centres of two particles come closer together than the sum of their radii.
More details are available in [27, 7].

The calculation of the caustics time in equation (14) was made for two heavy
spherical particles of identical size and therefore equal time constant
$\tau$. In the case of two particles with different $\tau$, equation (13) can
be rewritten in dimensional form for each particle, and the caustics radius
and time calculated. In Ref. [7, figure 3 therein], we showed that even very
small differences in the Stokes numbers can vastly change the “polydisperse
caustics radius”, denoted here by $r_{cp}$ and quantified in figure 4. The
particle starting at a smaller distance from the vortex centre is denoted by
the subscript $1$ and the other one by $2$. If the outer particle is larger,
the initial conditions determine whether the smaller one can overtake it. The
differences in figures 4 between the two initial conditions therefore merits
further study. When the inner particle is larger, $r_{cp}$ is larger by two
orders of magnitude than $r_{c}$. At first sight, it appears that our findings
in the earlier section are thus rendered irrelevant. This is emphatically not
the case: our conclusions hold for a distribution of particle sizes, unless
this distribution is extremely wide. The relevant caustics radii are still the
individual $r_{c}$’s. Well beyond this, each particle follows its outer
solution of executing an ever-more slowly widening Fermat spiral, and approach
velocities between particles, and therefore collision rates, are extremely
low. Nevertheless, $r_{cp}$ is technically the caustics radius for
polydisperse particles. A tangible effect of polydispersity is that the bigger
droplet’s radial velocity being larger than the smaller one’s, it can collide
with and collect more droplets rather easily, and so polydisperse caustics
afford a “rich get richer” scenario, by which droplets which are already
larger are strongly favoured to grow more easily.

Figure 4:  Nondimensional polydisperse caustics radii as a function of the
ratio of time scales $\xi=\tau_{2}/\tau_{1}$ of outer to inner particle. Small
differences in size lead to large variations. Left: particles initially
following circular streamlines; right: particles initially stationary. See
text for a discussion of relevance.

We shall use these findings on dynamics near one vortex in the next subsection
to estimate particle collision due to caustics in high Reynolds number
turbulence. We shall see that vortices of smaller and smaller scale contribute
more and more significantly to the total number of caustics events.

###  II.4 An estimate of single-vortex caustics in high Reynolds number
turbulence

A vast amount of new knowledge is being accrued about high Reynolds number
turbulence under homogeneous isotropic conditions, though these do not include
the thermodynamics of phase change. As we shall see in section IV, due to the
heat released during condensation, small-scale vorticity is likely to be more
prevalent in a cloud than in the common homogeneous isotropic turbulence, but
our analysis here based on ‘normal’ turbulence will provide a possibly
conservative estimate of droplet growth due to caustics near a single vortex.
For this we rely heavily on the excellent study by Buaria and Pumir [31],
showing how extreme vorticity follows universal behaviour at high Reynolds
number. They have studied Taylor microscale-based Reynolds numbers
$R_{\lambda}$ up to $1300$, but their demonstration of universality allows for
extrapolation to higher $R_{\lambda}$. Such extrapolation is of course subject
to validation against future simulations. The vorticity field is known to be
far more intermittent than the strain field, and we shall find this to be very
important. Following [31], we define $\Omega\equiv\omega_{i}\omega_{i}$. At
the small scales, we have [31]

$\frac{r_{v}}{\eta_{K}}\simeq(\Omega
t_{K}^{2})^{-\gamma/4}\quad\simeq\left(\frac{t_{K}\Gamma}{r_{v}^{2}}\right)^{-\gamma/2}$
(17)

where $\eta_{K}$ and $t_{K}$ are the Kolmogorov length and time scales,
$u_{v}$ is a typical velocity of a vortex of size $r_{v}$, and $\gamma$ is a
parameter obtained from direct numerical simulations. Here $\Gamma\sim
r_{v}u_{v}\sim\Omega^{1/2}r_{v}^{2}$, where $u_{v}$ is a velocity scale of the
vortex. As $R_{\lambda}\to\infty$, it is expected that $\gamma\to 1$, but in
our range of interest, we may take $\gamma\simeq 0.8$.

What follows is an order of magnitude calculation for vorticity at small
scales. For caustics to form, our single-vortex calculations have shown that
we must have the caustics radius $r_{c}$ to be larger than the vortex radius
$r_{v}$. Using equations (17) and (16), we have the simple requirement for
caustics, that

$\Omega\geq\frac{1}{\tau^{2}},\quad{\rm or}\ \Omega_{min}=\frac{1}{\tau^{2}},$
(18)

where $\Omega_{min}$ is the squared vorticity amplitude of the largest sized
vortex that can give rise to caustics. In other words, vortices of $St_{v}\geq
1$, i.e., of time scale equal to or faster than the particle time scale
participate in caustics. This estimate is exactly in line with the
observations of [32], in direct numerical simulation of particulate channel
flow, that significant voids due to centrifugation are formed for vortices
whose turnover time is comparable to the particle timescale. From equation 9,
this is equivalent to

$\frac{2\pi}{\omega_{v}}=\tau_{p},\text{or }St_{v}=\pi,$

which is comparable to the critical Stokes number of $2.3$ for caustics to
occur (see figure 2). Our contention is that what [32] see in their
simulations is a consequence of the caustics radius, and of (18).

As the Reynolds number increases, so does the intermittency in $\Omega$ [31],
and the fastest, or most extreme, time scale $t_{ext}$ in the flow reduces,
according to

$t_{ext}\sim t_{K}R_{\lambda}^{-\xi},$ (19)

where $\xi$ is an $O(1)$ number, changing very slowly with $R_{\lambda}$.
Again, for our purposes, $\xi\simeq 0.8$. Buaria & Pumir show that when
$\Omega$ is scaled by $t_{ext}^{-2}$, the tails of its probability density
over a range of $R_{\lambda}$ from $140$ to $1300$ collapse into a single
curve.

To simplify our estimate, we assume all particles are of uniform size. The
void area $a_{void}$ around a vortex, within which there are no particles, is
comparable to the area of the caustics region, i.e.,

$a_{void}\sim\Gamma\tau\sim u_{v}r_{v}\tau\sim\Omega^{1/2}r_{v}^{2}\tau\sim
x^{(1-\gamma)/2}\eta_{K}^{2}\frac{\tau}{t_{K}},$ (20)

where we have used equation (17), and defined $x\equiv\Omega t_{K}^{2}$.
Knowledge of the probability density $p(x)$ of $x$ enables us to obtain the
void fraction due to particle evacuation in the flow. Due to homogeneity, we
can obtain the void area fraction on any plane, and this would be the same
fraction in other planes, and therefore in the entire volume. Assuming
stationarity, $p(x)\Delta x|_{x1}$ is the fraction of the volume of the flow
which is occupied by vortices with $x$ in the range $x_{1}$ to $x_{1}+\Delta
x$. Similarly on any planar cross section in the flow the fraction of area
$A(x)$ occupied by vorticity in the corresponding range would be
$p(x){\color[rgb]{0,0,0}\Delta x_{x1}}$. This area is split into an average of
$A(x)/r_{v}^{2}$ vortices, where $r_{v}^{2}$ is the area of a vortex in order
of magnitude. On a plane of unit area, since

$\int_{0}^{\infty}p(x)dx=1,\quad{\rm we\ have}\quad\int_{0}^{\infty}A(x)dx=1,$
(21)

and the void fraction may be obtained, using equation (20), as

$f_{void}=\int_{x_{min}}^{\infty}\frac{p(x)}{r_{v}^{2}}a_{void}(x)\
dx\quad=\frac{\tau}{t_{K}}\int_{x_{min}}^{\infty}p(x)x^{1/2}\ dx.$ (22)

Note the lower limit of the integral, $x_{min}=\Omega_{min}t_{K}^{2}$, the
minimum $x$ for caustics, below which no void is created. This void fraction
is attained on very fast time scales, and may be taken as the instantaneous
fraction for a given vorticity field.

The idea is to find the squared amplitude of vorticity $\Omega_{min}$ beyond
which vortices will participate in caustics, and then use equations (20) and
(22) to obtain the fraction of volume in the turbulent flow occupied by voids.
A small fraction of the droplets which have evacuated these voids would
participate in caustics events, giving opportunity to coalesce and make bigger
drops. It is an oft quoted estimate of [19] that all we need for a runaway
growth of droplets into raindrops is a very small minority of larger drops, of
the order of one in a million drops. This is consistent with our finding of
“rich getting richer” by single vortex caustics. Using this, we may estimate
the generation of caustics by intermittent vortices in turbulent flows. This
exercise is carried out in the following section to give us an estimate of the
minimum turbulence levels which can produce rain.

###  II.5 Cloud estimates

In cloud physics we often encounter the term ‘droplet growth bottleneck’. This
is outlined as follows. A nascent cloud consists of large numbers of cloud
condensation nuclei (CCN), typically aerosol particles of the order of 1
micron in size, in a background slightly supersaturated in water vapour. Water
vapour at these concentrations can only condense into liquid water on CCNs,
and will not undergo spontaneous condensation in their absence. It is known
that diffusion of water vapour and condensation onto CCNs results, in a very
short time, in cloud droplets of size $\sim 10$ microns. But beyond this size,
their growth by condensation is too slow to create raindrops, as seen in
figure 5(a). This estimate is made accounting for radial diffusion of water
vapour onto to a droplet in a background which is 5 percent supersaturated,
allowing for the reduction in supersaturation near the droplet caused both by
condensation and the resulting increase in temperature. We mention that this
level of supersaturation in higher than seen in most clouds. A falling droplet
of instantaneous radius $a$, which is larger than the average radius $a_{0}$,
has another growth mechanism accessible to it: being larger, it will sediment
faster than the smaller droplets, since its dimensional terminal velocity

${\mathbf{v}}_{td}=g\tau_{p}{\mathbf{e}}_{g},\quad{\rm or\ in\ nondimensional\
terms,}\quad{\mathbf{v}}_{t}=-\frac{St}{Fr^{2}}.$ (23)

From this, and equation (3), we see that the bigger drop sediments faster than
the average by a factor $(a/a_{0})^{2}$, and can therefore participate in what
we would call “gravity induced caustics” events, i.e., catch up with and
overtake smaller droplets. The justification for assuming that a droplet
always falls at its terminal velocity is given as follows. If we drop the
history term in (2) but retain the gravity term, this equation may be readily
solved, as is often done in undergraduate fluid mechanics, to give

${\mathbf{v}}={\mathbf{v}_{t}}\left[1-e^{-t/St}\right]+{\mathbf{v}}_{0}e^{-t/St}.$
(24)

In a short time, of $O(St)$, the particle attains its terminal velocity, and
we may assume this process to be instantaneous.

Figure 5: Droplet growth (a) by condensation in a quiescent environment
supersaturated by $5\%$ in water vapour, and (b) by gravitational settling and
the associated collisions and coalescence. The initial radius of the larger
drop is $2^{1/3}a_{s}$ for the blue solid line and $1.01a_{s}$ for the red
dashed line. With $a_{s}=10$ microns, $N=10^{9}/m^{3}$ and $F=1$, and air-
water properties, the time shown is in seconds.

Consider a large three-dimensional box of quiescent fluid, with $N$ identical
small droplets per unit volume of size $a_{s}$ uniformly distributed in space.
The subscript $s$ stands for ‘small’. The suspension is dilute, so the
droplets do not interact with each other. Equation (24) tells us that all
these droplets are sedimenting steadily at their terminal velocity
${\mathbf{v}}_{ts}$. Now consider one larger droplet of radius $a$, slightly
bigger than $a_{s}$. This droplet falls faster than the smaller ones, and will
thus collide with the smaller droplets below it, which are within a cylinder
of radius $a+a_{s}$ from the centre of the large drop. We take it that a
fraction $F$ of collisions result in coalescence, allowing the bigger drop to
keep growing. Since the suspension is dilute, we may take it that the mean
collision time is far greater than the time scale of the big droplet. Using
equation (23), the growth of the larger droplet is then given by (from the
collision rate [20, in, e.g., equation 2 thereof])

$\frac{dy}{dt}=\frac{\pi
FNg\tau_{s}a_{s}^{2}}{3}\frac{(y-1)(y+1)^{3}}{y^{2}},$ (25)

where $y=a/a_{s}$. Figure 5(b) shows that the growth due to gravity by the
above equation, assuming every collision results in coalescence, i.e. $F=1$.
The smaller droplets are taken to be $10$ microns in radius, with a high
number density of $N=10^{9}/m^{3}$, which corresponds to cumulonimbus clouds.
The droplet size is seen to diverge at a finite time, but we note that that
when the drop is bigger than $\approx 50$ micron [e.g. 33], it will no longer
be in Stokes flow, and equation (25) would cease to be valid. If two smaller
drops coalesce to form the initial large drop, its radius would be
$2^{1/3}a_{s}$. Such a drop would grow to $10a_{s}$ in about 4 minutes,
whereas a drop initially at $1.01a_{s}$ would take twice that long. With a
smaller $N$, and with $F<1$, the growth rate would be proportionately slower.
Since rain is estimated to happen within about ten minutes of the initiation
of condensation, and unlike assumed above, collision efficiencies for these
small droplets are $F\ll 1$ [20, quoting Pruppacher and Klett], such growth is
too slow to create raindrops. We caution the reader again that the analysis
above is constrained by the applicability of the SMR equations 2. Large
droplets have larger drag forces acting on them with the drag increasing as
$|{\mathbf{v}}|^{2}$ for sufficiently large droplets. By accounting for finite
Reynolds number effects, it has been shown by others that droplets which are
already of size 50-100 microns can grow rapidly by gravitational collisions
and coalescence to millimetre size.

Figure 6:  Schematic of the droplet-growth bottleneck in clouds. Droplets
smaller than about 10 microns grow by diffusion and condensation, and those
bigger than 50-100 microns by gravitational settling and the associated
collisions and coalescence. Turbulence is said to bridge this gap. While
growth in the higher end of this size range could be enhanced by preferential
concentration, we propose single-vortex caustics as being effective near
intermittent vortices for the smaller end of this range.

We are thus confronted with the ‘droplet growth bottleneck’ question, which
asks how droplets grow rapidly from about 10 to 50 microns. Preferential
concentration in strain-dominated regions in turbulence is often presented as
the answer. But in simulations of turbulence, such preferential concentration
has only been seen to take place for particles of Stokes number of order
unity, based on the Kolmogorov timescale. Cloud droplets of $10$ micron radius
typically have Stokes numbers much smaller than unity in this scale. Figure 6
summarises this problem. Below, we make a case for single-vortex caustics to
be effective in the early part of the growth process beyond the diffusion
range.

We make estimates for droplet collisions due to caustics in a cumulonimbus
cloud, assuming that cloud turbulence is unchanged by phase change, and
assuming that the collapse of the tails of the vorticity probability density
seen by [31] may be extended to higher $R_{\lambda}$. In a cumulonimbus cloud
we have $Re\simeq 10^{8}$, i.e., $R_{\lambda}\simeq 10^{4}$, and turbulent
dissipation $\epsilon\simeq 0.1m^{2}/s^{3}$ [11]. This yields the large length
and velocity scales to be $L_{f}\simeq 300m$ and $U_{f}\simeq 3m/s$, and
Kolmogorov scales $\eta_{K}\simeq 3\times 10^{-4}m$, $u_{K}\simeq 0.03m/s$,
and therefore $t_{K}\simeq 0.01s$. We shall see below that this mechanism is
likely to be active in clouds of lower turbulence levels as well.

$R_{\lambda}$ | $t_{K}$ | $\eta_{K}$ | $x_{min}=\Omega_{min}t_{K}^{2}$ | $f_{void}$
---|---|---|---|---
$650$ | $0.1s$ |  $10^{-3}$ m | $2500$ | $\sim 10^{-10}$
$1300$ | $0.027s$ | $5\times 10^{-4}m$ | $200$ | $\sim 10^{-4}$
$10^{4}$ | $10^{-2}s$ | $3\times 10^{-4}m$ | $25$ | $O(1)$
Table 1:  Estimates of squared vorticity amplitude of the vortex of the
largest time scale that can participate in caustics, and the volume fraction
devoid of droplets as functions of the Taylor microscale Reynolds number.
Turbulence data provided in Fig.2a of [31] is used, along with extrapolation
to large $R_{\lambda}$ as suggested in that study. Calculations are based on a
particle time scale of 2 milliseconds, corresponding to a water droplet of
$10^{-6}m$ radius in air.

The probability density function for squared vorticity amplitude is provided
in Figure 2a of [31]. We use this to generate the estimates shown in Table 1,
for three values of $R_{\lambda}$. Consistent with estimates for turbulent
dissipation $\epsilon$ provided for clouds by [11], we prescribe
$\epsilon=0.01$ and $0.001m^{2}/s^{3}$ for $R_{\lambda}=1300$ and $650$
respectively. With these, and assuming constant dissipation for large scales,
(down to the Kolmogorov scale), we obtain the Kolmogorov time and length
scales shown in the table. Now a $10$ micron droplet in a cloud has a time-
scale $\tau\simeq 2$ milliseconds. Using equation (18), we obtain the minimum
squared vorticity amplitude, above which caustics can form, and display these
in the table. We note that all vortices which can form single-vortex caustics
in a cloud are far smaller than the Kolmogorov scale, i.e., they are intense
intermittent vortices. This points to the need for simulations which resolve
well below the Kolmogorov scale, up to $\eta_{ext}$. The highest $R_{\lambda}$
for which such well-resolved direct numerical simulations have been performed
in homogeneous isotropic turbulence is $1300$. Already at this Reynolds number
a large fraction of all vorticity is contained in the intermittent range
$\eta_{K}\geq\eta>\eta_{ext}$. From fig. 2a of [31], and using equations (20)
and (22), we can obtain the fraction of the flow devoid of particles for
$R_{\lambda}=650$ and $1300$. The case of $R_{\lambda}=10^{4}$, consistent
with deep convective clouds, is treated differently. Since direct numerical
simulations are not available, we can get at best a crude estimate. We assume
here that the collapse of the probability density curves seen in figure 3a of
[31] extends up to this Reynolds number. Under this assumption, we have
$p(x\times t_{K}^{2}/t_{ext}^{2})|_{10^{4}}=p(x^{*}\times
t_{K}^{2}/t_{ext}^{2})|_{1300}$, yielding $x^{*}/x=0.005$. This provides us
the transformation to use the probability density function of
$R_{\lambda}=1300$ to obtain an estimate for $R_{\lambda}=10^{4}$. But there
is a source of error here. The collapse of data is valid only for the tails of
the distribution, whereas $x^{*}$ in this case works out to be $0.131$.
Secondly, the data for low $x^{*}$ are not easy to read off in the data of
[31]. To make allowances for this we obtain only an order of magnitude
estimate for void fraction at $R_{\lambda}=10^{4}$, as being of $O(1)$. Our
scaling leads us to expect that a significant fraction of all droplets at this
Reynolds number will be participating in caustics events, leading to rapid
rainfall initiation. But indeed nothing short of actual cloud simulations
which are resolved well past the Kolmogorov scale at this high Reynolds number
will provide us the real estimate. In particular, the collapse with
$R_{\lambda}$ seen up to $1300$ may not persist up to another order of
magnitude increase.

The presence of voids in the spatial distribution of inertial particles is
well-known [34, 35]. In particular, Bec & Chétrite [35] suggest a model in
which vortices eject particles of timescales comparable to their turnover
times, and show that this model is able to reproduce observed statistics of
spatial distributions of particles. Our ideas are similar, except that our
calculations are based on an exact model for the flow and its coupling with
particle inertia.

Clouds are not always this turbulent and an $R_{\lambda}$ of $1300$ may not be
unrealistic in some cloud situations. So we restrict our further discussion to
lower Reynolds numbers. In a volume of a cubic meter within a cumulonimbus
cloud, there are $O(10^{9})$ droplets of 10 micron radius, so by our estimate,
$10^{5}$ drops per cubic metre may be termed as caustics droplets. In other
clouds this may be an order of magnitude lower, since $N$ is lower. We have
found [27, 7] that caustics droplets that do undergo collisions do so on time
scales comparable to the inertial particle timescale $\tau$, and within a
distance comparable to the caustics radius $\sqrt{\Gamma\tau}$ of the vortex.
Allowing for a collision efficiency $F$ of $0.01-0.1$ we expect about 1 in
every hundred thousand or a million droplets to coalesce with another, and
become a factor of $2^{1/3}$ larger in radius than average. These larger drops
have a correspondingly higher Stokes number, and larger caustics radius, so
they can be thrown out further and attain higher momentum than the average
droplet, with enhanced probability for further collisions and growth due to
both turbulence and gravity. The argument, following Kostinski & Shaw [19], is
that it is the initial collisions that are crucial in determining the ‘luck’
of the droplets that ultimately become rain, since the later collisions become
progressively more likely. Therefore, a mechanism such as ours that explain
the initial collisions can be a crucial part of the story. We may conclude
that, by this mechanism, single-vortex caustics are extremely unlikely to
initiate rainfall at $R_{\lambda}$ of $650$, but at $1300$ and beyond, they
can initiate rain, with increasing rapidity for higher Reynolds numbers.

##  III Basset-Boussinesq history

In the above we had worked with the SMR equation, whereas particle dynamics
are described by the MR equation (2). We notice that the next term in the
power series of $St$ is the Basset-history force, being a half-power of Stokes
smaller than the Stokes drag term. This term is usually neglected, not because
we are convinced about its smallness of magnitude, but simply because it is
too hard to compute, if one were to follow the naive approach of performing
the integral at every instant of time for every particle [36]. It is clear
that this approach is not feasible for any large flow. For example a cloud
contains a billion droplets per cubic metre, and performing this integral for
every single droplet, each time from $t=0$ to the present time, involves not
only linearly increasing computation time for each time interval, but also
impossibly large storage, since one would need to store, for every particle,
its acceleration and also that of the fluid at the location of the particle at
all past instants. We proposed a novel approach [37] to solving this problem.
The main advance provided by this approach is that memory costs are
drastically reduced, and the memory burden is constant in time. It also
involves significantly lower time costs. The reader is referred to [37] for
more details, but the basic idea is explained here. For several simple flow
situations, our approach affords explicit analytical solutions to the MR
equation which were unknown hitherto, and a numerical method for the general
case. This was made possible by noticing that the term within the square
brackets in the MR equation (2), is just the negative of the half (fractional)
derivative in time of ${\mathbf{q}}\equiv{\mathbf{v}}-{\mathbf{u}}$. Following
the idea of Vishal Vasan, we now write a heat equation in a fictional
coordinate $x$, as

${\mathbf{q}}_{t}\equiv{\mathbf{q}}_{xx},$ (26)

so the MR equation (2) may be rewritten simply as

$\dot{\mathbf{v}}=-\frac{{\mathbf{q}}}{St}+\frac{3}{(2\beta St)^{1/2}}\
{\mathbf{q}}_{x}+\frac{1}{Fr^{2}}\mathbf{e}_{g}.$ (27)

Note that equation (27) is valid only on the $x=0$ boundary of equation (26),
and forms its boundary condition. An excellent case for study here is that of
a particle sedimenting under gravity in a quiescent ambient. Equation (24)
showed how a particle attains its terminal velocity, given by (23), within a
few particle time-scales. In [37] on the other hand, we show how the BBH force
alters the particle’s behaviour in a fundamental way. By our approach, we can
write down the complete analytical solution of (27), which includes the
effects of the BBH, for this problem. For long times, the solution may be
expressed as

${\mathbf{v}}={\mathbf{v}_{t}}-3{\mathbf{v}}_{t}\sqrt{\frac{St}{2\beta\pi
t}}+O(t^{-3/2}).$ (28)

Accounting for the BBH term does not change the terminal velocity, but the
particle approaches it algebraically, i.e., slowly, rather than by the
exponentially vanishing transients we obtained without history.

In particular, in [37] we solved the complete MR equation (27) to obtain the
effect of BBH on the growth by collision and coalescence of a larger droplet
sedimenting through a sea of identical smaller droplets. The particle density
ratio $\beta$ chosen there, $O(10)$, was relevant to other settings of
droplet-laden flows, but not directly to clouds, where $\beta=10^{3}$. For a
dilute suspension, the inter-collision time is large enough that we may assume
equation (28) to be valid. The result for $\beta=10^{3}$ is plotted in figure
7 and compared against the estimate given by equation (25).

Figure 7:  Growth in time of an initially larger droplet as it sediments
through a uniform distribution of smaller droplets, accreting them as it does
so. Black solid line: with Basset history. Red dashed line: without Basset
history.

Thus, while neglecting the Basset history term may in general lead to
incorrect results, Fig. 7 shows that for the small droplet number densities
and large density ratios typical of clouds, the BBH term is indeed negligible,
and the algebraic growth of settling velocity (Eq. 28) does not affect the
result.

##  IV Inertia–Thermodynamics coupling

Thus far, we have concerned ourselves with the influence of the flow on the
droplets, and neglected the influence of the droplets on the flow. Typically,
neglecting the momentum coupling due to the droplets on the flow is
justifiable when the volume fractions are sufficiently small. In clouds
however, although the volume fraction is only $O(10^{-6})$, droplets can have
a profound effect on the flow. This is because droplets in clouds also act as
nuclei for condensation, and the accompanying heat release and the consequent
buoyancy modification cannot be neglected.

Droplets, as we have seen, are evacuated from the vicinity of strong vortices.
Since droplets also act as condensation nuclei, due to their evacuation, we
expect practically no condensation to occur in the void regions surrounding
strong vortices in a supersaturated parcel (such as a rising cumulus cloud).
What happens in a supersaturated environment when most CCN evacuate the
vicinity of strong vortices is that condensation takes place selectively
outside these voided spaces, and latent heat is released outside, and the
regions occupied by droplets becomes warmer. So strong vortices are now
situated within cold patches which are about the size of the caustics radius.
If the temperature differences within and outside the patches are large
enough, the flow changes qualitatively due to the presence of droplets. This
is because buoyancy effects kick in, and baroclinic torque is created. As a
result, such flows with phase change display a larger amount of kinetic energy
in the small scales than turbulent flows without. The increase in kinetic
energy is due to the conversion by baroclinic torque of the increased
potential energy due to temperature inhomogeneity, which in turn is due to
selective condensation outside void regions. A sample result is shown in
figure 8, which is taken from [38]. There, we showed that this effect can
occur if the vortex strengths and droplet sizes are such that the product of
the inertial Stokes number $St_{v}$ and the condensation Stokes number
$St_{s}$ is sufficiently large. In a model simulation consisting of just two
vortices, we found that if $St_{v}St_{s}>\sim 10^{2}$, the flow with phase
change is qualitatively different from the flow without. The reader is
referred to [38] for a scaling argument for why this product is important.
Intuitively, we may argue that if the condensation time scale is too small,
condensation will be complete before droplets have a chance to centrifuge out
of vortices. On the other hand, if the droplet time scale is too small, their
inertia will be too weak to create void spaces. We may now make an estimate,
using data from very large simulations, whether this mechanism can operate in
typical conditions in cumulus clouds in the Earth’s atmosphere.

The typical inertial timescale in cumulus clouds, where the typical droplet
has a radius of $O(10)$ micron, is $O(10^{-3})$s. The timescale for phase
change depends on the concentration of droplets and weakly on the ambient
temperature, and can be estimated to be $O(10)$s, using

$\tau_{s}=\frac{C\rho_{s}^{0}}{4\pi Na},$

where $C\approx 10^{7}m\cdot skg^{-1}$ is a thermodynamic constant, $N$ is the
number density of droplets of radius $a$, and $\rho_{s}^{0}$ is the saturation
vapour density at some base temperature $T_{0}$. Thus, the product of these
timescales is $O(10^{-2})$s2. The Kolmogorov timescale, for a dissipation rate
of $\epsilon=10^{-2}m^{2}s^{-3}$, typical of cumulus and cumulonimbus clouds
[11], is $\tau_{k}=O(10^{-3/2})$s, giving the Kolmogorov ‘vorticity scale’
$\omega_{k}=10^{3/2}s^{-1}$. The product of the inertial and phase-change
Stokes numbers is

$St_{p}\times St_{s}=\omega_{k}^{2}\tau_{p}\tau_{s}=O(10).$

Thus, the mechanism in [38], of flow dynamics being changed due to droplets,
is not expected to operate for typical Kolmogorov vortices. However,
intermittent vortices, we have seen, have vorticities orders of magnitude
larger, and thus the product above will be sufficiently large. And since the
void fraction becomes significant only at high Reynolds numbers, this is an
additional requirement. Our study points to the creation of a larger fraction
of smaller and stronger vortices (increased intermittency) if this mechanism
were to be in operation.

Figure 8:  The vorticity at $t=50$, similar to figures 16 and 17 of [38]. The
domains are $20\times 20$ in the $x$ and $y$ directions, and vorticity
contours are plotted between the values $(-1,1)$. The nondimensional Reynolds
number in the simulations is $Re=2\times 10^{4}$, the phase-change timescale
is $St_{s}=|\omega|\tau_{s}=200$ and the droplet Stokes number is
$St=|\omega|\tau=1$, where $|\omega|$ is the vorticity scale of the flow, and
$\tau_{s}$ and $\tau$ are the timescales of condensation and droplet inertia
respectively. The simulations were performed using a pseudospectral method in
the $\omega-\psi$ formulation in two dimensions, with $2/3$rds de-aliasing and
a second order exponential time-differencing (ETD-2) scheme. The details may
be found in [38].

Since, as we have argued, the primary requirement for the mechanism to operate
is the strong vortices associated with intermittent behaviour, this
modification of the turbulence is self-sustaining. Explicit evidence for this
is the subject of an ongoing study.

Before we end, we mention how phase change can bring about two beautiful cloud
formations: mammatus and asperitas, described in [14] and [39] respectively.
As small droplets sediment below the clouds, evaporation can create cooling
and a denser layer of air can form just below the clouds. Once this layer
becomes pronounced enough, Rayleigh-Taylor instability can take place, giving
rise to startling lobe formations, which are called mammatus clouds. When
shear is added to this mix, we get a more disordered cloud formation with
sharper structures below the cloud base, with some alignment in the shear
direction. These are the rare and recently notified asperitas clouds.

##  V Summary, conclusions, and discussion about future work

In this article we have examined the consequences of inertial particles (or
droplets) centrifuging out of vortical structures in turbulent flow. We have
restricted our attention to small spherical droplets or particles in the
unsteady Stokes limit in dilute suspension, whose dynamics may be described by
the Maxey-Riley equation where Faxen correction terms have been neglected.
Given the increasing evidence that turbulence at the small scale contains a
significant fraction of tubelike vortices we have treated the Lamb-Oseen
vortex as the building block of turbulence, and written down, in (6) and (7),
the MR equations for a particle in this flow. When taken to the point vortex
limit, these equations can be shown to support a boundary-layer structure.
Single-vortex caustics can form only when
$r_{d}<r_{cd}\simeq(\Gamma\tau)^{1/2}$. Thus, caustics formation is only
relevant for vortices whose size $r_{v}<r_{c}$, which we have shown here to
consist of vortices smaller than the Kolmogorov scale. Our arguments do not
change significantly by vortex stretching, as in the case of Burgers vortices
where we have shown that greater stretching rates lead to faster caustic
formation [8]. This is important, given the prevalence of vortex stretching in
turbulent flows [40, 41].

The connection we find between caustics formation and vortex evacuation
enables us to calculate void fractions in the flow, made up of circular
(cylindrical in three dimensions) patches around strong vortices which are
devoid of particles. Evacuating particles have a much bigger chance of
colliding with other particles at short times and, in the case of droplets,
coalesce and grow. Using turbulence data at very high Reynolds numbers in
conjunction with these ideas, we make the audacious propositions that single-
vortex caustics could be important in rain initiation, and that rain
initiation needs $R_{\lambda}>10^{3}$.

We then showed, extending the computation of [37] to a density ratio
$\beta=10^{3}$ representative of clouds, that the omission of the Basset-
history term is justifiable. While neglecting the BBH term alters the manner
in which terminal velocity is attained by a falling droplet is algebraic
rather than exponential, the effect of this for dilute suspensions with cloud-
like density ratios is small. Finally, we describe how cloud turbulence, due
to phase change, can be different from turbulence in other situations.

Several leads for future work emerge from this discussion. First of all,
larger and larger Reynolds number direct numerical simulations, which solve
the Navier-Stokes equations along with droplet dynamics and the thermodynamics
of phase change, are urgently needed, to understand rain initiation in clouds
and many other basic physical processes in clouds. Since it is unimaginable
that such simulations can include the dynamics of each droplet, a clever
superdrop method which includes the physics of caustics is crucial to develop.
We note that there are many other flow situations where BBH will be important,
and our conclusion above is only applicable to particles which are far denser
than the fluid, and moreover in dilute suspension. Indeed, the effects of the
history term are striking when the density ratio is smaller. For example,
carbonaceous matter sinking through the ocean has density comparable to the
background fluid, and correct estimates of carbon sequestration by this
process, and its part in mitigating climate change, could be important.

The formation of different cloud shapes contains, apart from breathtaking
beauty, much physics, including hitherto undescribed instabilities, and again
this area is just opening up for studies. Our findings are also relevant in
sprays and industrial flows, and the effects of zones devoid of particles in
such flows need to be studied.

###### Acknowledgements.

Many thanks are due to Alain Pumir for sharing some data and clarifying some
aspects of their recent study, which forms a crucial input for some of our
arguments. RG gratefully acknowledges support of the Department of Atomic
Energy, Government of India, under project no. RTI4001. SR was supported at
Nordita under the Swedish Research Council grant no. 638-2013-9243. Nordita is
partially funded by Nordforsk.

## References

  * Shaw  _et al._ [1998] R. A. Shaw, W. C. Reade, L. R. Collins, and J. Verlinde, Preferential concentration of cloud droplets by turbulence: Effects on the early evolution of cumulus cloud droplet spectra, Journal of the Atmospheric Sciences 55, 1965 (1998).
  * Falkovich  _et al._ [2002] G. Falkovich, A. Fouxon, and M. G. Stepanov, Acceleration of rain initiation by cloud turbulence, Nature 419, 151 (2002).
  * Wilkinson and Mehlig [2005] M. Wilkinson and B. Mehlig, Caustics in turbulent aerosols, Europhysics Letters (EPL) 71, 186 (2005).
  * Wilkinson  _et al._ [2006] M. Wilkinson, B. Mehlig, and V. Bezuglyy, Caustic activation of rain showers, Physical Review Letters 97, 048501 (2006).
  * Mehlig  _et al._ [2007] B. Mehlig, V. Uski, and M. Wilkinson, Colliding particles in highly turbulent flows, Physics of Fluids 19, 098107 (2007).
  * Vosskuhle  _et al._ [2014] M. Vosskuhle, A. Pumir, E. Lévêque, and M. Wilkinson, Prevalence of the sling effect for enhancing collision rates in turbulent suspensions, Journal of Fluid Mechanics 749, 841 (2014).
  * Deepu  _et al._ [2017] P. Deepu, S. Ravichandran, and R. Govindarajan, Caustics-induced coalescence of small droplets near a vortex, Physical Review Fluids 2, 024305 (2017).
  * Agasthya  _et al._ [2019] L. Agasthya, J. R. Picardo, S. Ravichandran, R. Govindarajan, and S. S. Ray, Understanding droplet collisions through a model flow: Insights from a burgers vortex, Physical Review E 99, 063107 (2019).
  * Balachandar and Eaton [2010] S. Balachandar and J. K. Eaton, Turbulent dispersed multiphase flow, Annu. Rev. Fluid Mech 42, 111 (2010).
  * Shaw [2003] R. A. Shaw, Particle-turbulence interactions in atmospheric clouds, Annual Review of Fluid Mechanics 35, 183 (2003).
  * Grabowski and Wang [2013] W. W. Grabowski and L.-P. Wang, Growth of cloud droplets in a turbulent environment, Annual Review of Fluid Mechanics 45, 293 (2013).
  * Bec  _et al._ [2014] J. Bec, H. Homann, and S. S. Ray, Gravity-driven enhancement of heavy particle clustering in turbulent flow, Physical Review Letters 112, 1 (2014).
  * Ireland  _et al._ [2016] P. J. Ireland, A. D. Bragg, and L. R. Collins, The effect of reynolds number on inertial particle dynamics in isotropic turbulence, Part II: simulations with gravitational effects, Journal of Fluid Mechanics 796, 659 (2016).
  * Ravichandran  _et al._ [2020] S. Ravichandran, E. Meiburg, and R. Govindarajan, Mammatus cloud formation by settling and evaporation, J. Fluid Mech. 899, A27 (2020).
  * Omand  _et al._ [2020] M. M. Omand, R. Govindarajan, J. He, and A. Mahadevan, Sinking flux of particulate organic matter in the oceans: Sensitivity to particle characteristics, Scientific Reports 2020 10:1 10, 1 (2020).
  * Crowley [1971] J. M. Crowley, Viscosity-induced instability of a one-dimensional lattice of falling spheres, Journal of Fluid Mechanics 45, 151 (1971).
  * Chajwa  _et al._ [2020] R. Chajwa, N. Menon, S. Ramaswamy, and R. Govindarajan, Waves, algebraic growth, and clumping in sedimenting disk arrays, Physical Review X 10, 041016 (2020).
  * Lohmann and Feichter [2005] U. Lohmann and J. Feichter, Global indirect aerosol effects: A review, Atmospheric Chemistry and Physics 5, 715 (2005).
  * Kostinski and Shaw [2005] A. B. Kostinski and R. A. Shaw, Fluctuations and luck in droplet growth by coalescence, Bulletin of the American Meteorological Society 86, 235 (2005).
  * Wilkinson [2016] M. Wilkinson, Large deviation analysis of rapid onset of rain showers, Physical Review Letters 116, 018501 (2016).
  * Maxey and Riley [1983] M. R. Maxey and J. J. Riley, Equation of motion for a small rigid sphere in a nonuniform flow, Physics of Fluids 26, 883 (1983).
  * Bec [2003] J. Bec, Fractal clustering of inertial particles in random flows, Physics of Fluids 15, L81 (2003).
  * Bec  _et al._ [2005] J. Bec, A. Celani, M. Cencini, and S. Musacchio, Clustering and collisions of heavy particles in random smooth flows, Physics of Fluids 17, 073301 (2005).
  * Falkovich and Pumir [2007] G. Falkovich and A. Pumir, Sling effect in collisions of water droplets in turbulent clouds, Journal of the Atmospheric Sciences 64, 4497 (2007).
  * Maxey [1987] M. R. Maxey, The motion of small spherical particles in a cellular flow field, Physics of Fluids 30 (1987).
  * Squires and Eaton [1991] K. D. Squires and J. K. Eaton, Preferential concentration of particles by turbulence, Physics of Fluids A: Fluid Dynamics 3, 1169 (1991).
  * Ravichandran and Govindarajan [2015] S. Ravichandran and R. Govindarajan, Caustics and clustering in the vicinity of a vortex, Phys. Fluids 27 (2015).
  * Ramadugu  _et al._ [2022] R. Ramadugu, P. Perlekar, and R. Govindarajan, Surface tension as the destabiliser of a vortical interface, Journal of Fluid Mechanics 936, 10.1017/JFM.2022.97 (2022).
  * Moriconi  _et al._ [2022] L. Moriconi, R. M. Pereira, and V. J. Valadão, Circulation statistics and the mutually excluding behavior of turbulent vortex structures,  (2022).
  * Sudharsan  _et al._ [2016] M. Sudharsan, S. L. Brunton, and J. J. Riley, Lagrangian coherent structures and inertial particle dynamics, Physical Review E 93, 033108 (2016).
  * Buaria and Sreenivasan [2022] D. Buaria and K. R. Sreenivasan, Intermittency of turbulent velocity and scalar fields using three-dimensional local averaging, PHYSICAL REVIEW FLUIDS 7, 72601 (2022).
  * Motoori  _et al._ [2022] Y. Motoori, C. Wong, and S. Goto, Role of the hierarchy of coherent structures in the transport of heavy small particles in turbulent channel flow, Journal of Fluid Mechanics 942, A3 (2022).
  * Good  _et al._ [2014] G. Good, P. Ireland, G. Bewley, E. Bodenschatz, L. Collins, and Z. Warhaft, Settling regimes of inertial particles in isotropic turbulence, Journal of Fluid Mechanics 759, R3 (2014).
  * Bec  _et al._ [2007] J. Bec, L. Biferale, M. Cencini, A. Lanotte, S. Musacchio, and F. Toschi, Heavy particle concentration in turbulence at dissipative and inertial scales, Physical Review Letters 98, 084502 (2007).
  * Bec and Chétrite [2007] J. Bec and R. Chétrite, Toward a phenomenological approach to the clustering of heavy particles in turbulent flows, New Journal of Physics 9, 77 (2007).
  * Haller [2019] G. Haller, Solving the inertial particle equation with memory, Journal of Fluid Mechanics 874, 1 (2019).
  * Prasath  _et al._ [2019] S. G. Prasath, V. Vasan, and R. Govindarajan, Accurate solution method for the Maxey-Riley equation, and the effects of Basset history, Journal of Fluid Mechanics 868, 428 (2019).
  * Ravichandran and Govindarajan [2017] S. Ravichandran and R. Govindarajan, Vortex-dipole collapse induced by droplet inertia and phase change, Journal of Fluid Mechanics 832, 745 (2017).
  * Ravichandran and Govindarajan [2022] S. Ravichandran and R. Govindarajan, Instability driven by settling and evaporation in a shear flow: A model for asperitas clouds, Physical Review Fluids 7, 010501 (2022).
  * Ohkitani [2002] K. Ohkitani, Numerical study of comparison of vorticity and passive vectors in turbulence and inviscid flows, Physical Review E 65, 046304 (2002).
  * Hamlington  _et al._ [2008] P. E. Hamlington, J. Schumacher, and W. J. A. Dahm, Local and nonlocal strain rate fields and vorticity alignment in turbulent flows, Physical Review E - Statistical, Nonlinear, and Soft Matter Physics 77, 026303 (2008).