format: text

text:
# The clustering instability of inertial particles spatial distribution in turbulent flows

## I Introduction

Formation and evolution of aerosols and droplets inhomogeneities (clusters) are of fundamental significance in many areas of environmental sciences, physics of the atmosphere and meteorology (e.g., smog and fog formation, rain formation), transport and mixing in industrial turbulent flows (like spray drying, pulverized-coal-fired furnaces, cyclone dust separation, abrasive water-jet cutting) and in turbulent combustion (see, e.g., CST98 ; T77 ; C80 ; S86 ; EKPR97 ; SR98 ; EKRI98 ; VY00 ). The reason is that the direct, hydrodynamic, diffusional and thermal interactions of particles in dense clusters strongly affect the character of the involved phenomena. Thus, e.g., enhanced binary collisions between cloud droplets in dense clusters can cause fast broadening of droplet size spectrum and rain formation (see, e.g., VY00 ). Another example is combustion of pulverized coal or sprays whereby reaction rate of a single particle or a droplet differs considerably from a reaction rate of a coal particle or a droplet in a cluster (see, e.g., AR92 ; ARD94 ).

Analysis of experimental data shows that spatial distributions of droplets in clouds are strongly inhomogeneous (see, e.g., PB89 ; KM93 ; B99 ; KS01 ). Small-scale inhomogeneities in particles distribution were observed also in laboratory turbulent flows SE91 ; EF94 ; FKE94 ; HL00 .

It is well-known that turbulence results in a relaxation of inhomogeneities of concentration due to turbulent diffusion, whereas the opposite process, e.g., a preferential concentration (clustering) of droplets and particles in turbulent fluid flow still remains poorly understood.

In this paper we suggest a theory of clustering of particles and droplets in turbulent flows. The clusters of particles are formed due to an instability of their spatial distribution suggested in Ref. EKRI96  and caused by a combined effect of a particle inertia and a finite velocity correlation time. Particles inside turbulent eddies are carried out to the boundary regions between them by inertial forces. This mechanism of the preferential concentration acts in all scales of turbulence, increasing toward small scales. An opposite process, a relaxation of clusters is caused by a scale-dependent turbulent diffusion. The turbulent diffusion decreases towards to smaller scales. Therefore, the clustering instability dominates in the Kolmogorov inner scale $`\eta `$, which separates inertial and viscous scales. Exponential growth of the number of particles in the clusters is saturated by their collisions.

In our previous study EKRI96  we suggested and analyzed qualitatively an idea that inertia of particles may lead to their clustering. Later this idea was questioned by our quantitative analysis EKRS00 ; EKRS99  of the Kraichnan model of turbulent advection of particles by the delta-correlated in time random velocity field. It was proved that the clustering of inertial particles does not occur in the Kraichnan model. The latter result may be considered as counterexample.

The main quantitative result of the theory of clustering instability of inertial particles, suggested in this paper, is the existence of this instability under some conditions that we determined. We showed that the inertia of the particles is only one of the necessary conditions for particles clustering in turbulent flow. In the present study we found a second necessary condition for the clustering instability: a finite correlation time of the fluid velocity field which in the suggested theory results in a nonzero divergence of the field of Lagrangian trajectories. This time is equal zero in the above mentioned Kraichnan model (see K68 ) which was the reason for the disappearance of the instability in this particular model.

In this study we used a model of the turbulent velocity field with a finite correlation time which drastically changes the dynamics of inertial particles. In the framework of this model of the velocity field we rigorously derived the sufficient conditions for the clustering instability. We demonstrated the existence of the new phenomena of strong and weak clustering of inertial particles in a turbulent flow. These two types of the clustering instabilities have different physical meaning and different physical consequences in various phenomena. We computed also the instability thresholds which are different for the strong and weak clustering instabilities.

## II Qualitative analysis of strong and weak clustering

### II.1 Basic equations in the continuous media approximation

In this study we used the equation for the number density $`n(t,𝒓)`$ of particles advected by a turbulent velocity field $`𝒖(t,𝒓)`$:

$`{\displaystyle \frac{n(t,𝒓)}{t}}+\mathbf{}[n(t,𝒓)𝒗(t,𝒓)]=D\mathrm{\Delta }n(t,𝒓),`$ (1)
where $`D=kT/6\pi \nu \rho a`$ is the coefficient of molecular (Brownian) diffusion, $`\nu `$ is the fluid kinematic viscosity, $`\rho `$ and $`T`$ are the fluid density and temperature, respectively, $`a`$ is the radius of a particle and $`k`$ is the Boltzmann constant. Due to inertia of particles their velocity $`𝒗(t,𝒓)𝒖(t,𝒓),`$ e.g., the field $`𝒗(t,𝒓)`$ is not divergence-free even for div $`𝒖=0`$ (see EKRI96 ). Equation (1) implies conservation of the total number of particles in a closed volume. Consider

$`\mathrm{\Theta }(t,𝒓)=n(t,𝒓)\overline{n},`$ (2)
the deviation of $`n(t,𝒓)`$ from the uniform mean number density of particles $`\overline{n}`$. Equation for $`\mathrm{\Theta }(t,𝒓)`$ follows from Eq. (1):

$`{\displaystyle \frac{\mathrm{\Theta }(t,𝒓)}{t}}`$ $`+`$ $`[𝒗(t,𝒓)\mathbf{}]\mathrm{\Theta }(t,𝒓)`$
$`=`$ $`\mathrm{\Theta }(t,𝒓)\text{div}𝒗(t,𝒓)+D\mathrm{\Delta }\mathrm{\Theta }(t,𝒓).`$
Here we assumed that the mean particles velocity is zero. We also neglected the term $`\overline{n}\text{div}𝒗`$ describing an effect of an external source of fluctuations. This term does not affect the growth rate of the instability. In the present study we investigate only the effect of self-excitation of the clustering instability, and we do not consider an effect of the source term on the dynamics of fluctuations. The source term $`\overline{n}\text{div}𝒗`$ causes another type of fluctuations of particle number density which are not related with an instability and are localized in the maximum scale of turbulent motions. A mechanism of these fluctuations is related with perturbations of the mean number density of particles by a random divergent velocity field. The magnitude of these fluctuations is much lower than that of fluctuations which are caused by the clustering instability.

In our qualitative analysis of the problem we use Eq. (II.1) written in a co-moving with a cluster reference frame. Formally, this may be done using the Belinicher-L’vov (BL) representation (for details, see LP2 ; LP1 ). Let $`𝝃_\text{L}(t_0,𝒓|t)`$ will be Lagrangian trajectory of the reference point \[in the particle velocity field $`𝒗(t,𝒓)`$\] located at $`𝒓`$ at time $`t_0`$ and $`𝝆_\text{L}(t_0,𝒓|t)`$ be an increment of the trajectory:

$`𝝆_\text{L}(t_0,𝒓|t)`$ $`=`$ $`{\displaystyle _{t_0}^t}𝒗[\tau ,𝝃_\text{L}(t_0,𝒓|\tau )]\text{d}\tau ,`$ (4)
$`𝝃_\text{L}(t_0,𝒓|t)`$ $``$ $`𝒓+𝝆_\text{L}(t_0,𝒓|t).`$
By definition $`𝝆_\text{L}(t_0,𝒓|t_0)=0`$ and $`𝝃(t_0,𝒓|t_0)=𝒓`$. Define as $`𝒓_0`$ a position of a center of a cluster at the “initial” time $`t_0=0`$ (for the brevity of notations hereafter we skip the label $`t_0`$ ) and consider a “co-moving” reference frame with the position of the origin at $`𝜻_0(t)𝜻(𝒓_0|t)`$. Then BL velocity field $`\stackrel{~}{𝒗}(𝒓_0|t,𝒓)`$ and BL velocity difference $`𝑾(𝒓_0|t,𝒓)`$ are defined as

$`\stackrel{~}{𝒗}(𝒓_0|t,𝒓)`$ $``$ $`𝒗[t,𝒓+𝝆_\text{L}(𝒓_0|t)],`$ (5)
$`𝑾(𝒓_0|t,𝒓)`$ $``$ $`\stackrel{~}{𝒗}(𝒓_0|t,𝒓)\stackrel{~}{𝒗}(𝒓_0|t,𝒓_0).`$ (6)
Actually the BL representation is very similar to the Lagrangian description of the velocity field. The difference between the two representations is that in the Lagrangian representation one follows the trajectory of every fluid particle $`𝒓+𝝆_\text{L}(𝒓,t_0|t)`$ (located at $`𝒓`$ at time $`t=t_0`$), whereas in the BL-representation there is a special initial point $`𝒓_0`$ (in our case the initial position of the center of the cluster) whose trajectory determines the new coordinate system (see LP2 ; LP1 ). With time the BL-field $`\stackrel{~}{𝒗}(𝒓_0|t,𝒓)`$ becomes very different from the Lagrangian velocity field. It must be noted that the simultaneous correlators of both, the Lagrangian and the BL-velocity fields, are identical to the simultaneous correlators of the Eulerian velocity $`𝒗(𝒓,t)`$. The reason is that for stationary statistics the simultaneous correlators do not depend on $`t`$, and in particular one can assume $`t=t_0`$.

Similar to Eq. (5) let us introduce BL representation for $`\mathrm{\Theta }(t,𝒓)`$:

$`\stackrel{~}{\mathrm{\Theta }}(𝒓_0|t,𝒓)\mathrm{\Theta }[t,𝒓+𝝆_\text{L}(𝒓_0|t)].`$ (7)
In BL variables defined by Eqs. (5)-(7), Eq. (II.1) reads:

$`{\displaystyle \frac{\stackrel{~}{\mathrm{\Theta }}(𝒓_0|t,𝒓)}{t}}+[𝑾(𝒓_0|t,𝒓)\mathbf{}]\stackrel{~}{\mathrm{\Theta }}(𝒓_0|t,𝒓)`$
$`=`$ $`\stackrel{~}{\mathrm{\Theta }}(𝒓_0|t,𝒓)\text{div}𝑾(𝒓_0|t,𝒓)+D\mathrm{\Delta }\stackrel{~}{\mathrm{\Theta }}(𝒓_0|t,𝒓).`$
The difference between Eqs. (II.1) and (II.1) is that Eq. (II.1) involves only velocity difference (6) in which the velocity $`\stackrel{~}{𝒗}(𝒓_0|t,𝒓_0)`$ of the cluster center is subtracted.

### II.2 Rigid-cluster Approximation

Consider qualitatively a time evolution of different statistical moments

of the deviation $`\mathrm{\Theta }(t,𝒓)`$ defined by

$`_q(t)|\mathrm{\Theta }(t,𝒓)|^q_v,`$ (9)
assuming that at the initial time, $`t=0`$, the spatial distribution of particles is almost homogeneous, all moments $`_q(0)`$ are small, where $`_v`$ denotes the ensemble averaging over random velocity field $`𝒗`$. In order to eliminate the kinematic effect of sweeping of the cluster as a whole we consider Eq. (II.1) in the BL-representation, Eq. (II.1). Since the simultaneous moments of any field variables in the Eulerian and in the BL-representations coincide, the moments $`_q(t)`$ can be written as

$`_q(t)=|\stackrel{~}{\mathrm{\Theta }}(𝒓_0|t,𝒓)|^q_v.`$ (10)
Our conjecture is that on a qualitative level we can consider the role of each term in the Eq. (II.1) separately, assuming some reasonable, time independent, frozen shape $`\theta (x)`$ of a distribution $`\stackrel{~}{\mathrm{\Theta }}(𝒓_0|t,𝒓)`$ inside a cluster:

$`\stackrel{~}{\mathrm{\Theta }}(𝒓_0|t,𝒓)=A(t)\theta \left({\displaystyle \frac{|𝒓𝒓_0|}{\mathrm{}_{\text{cl}}}}\right).`$ (11)
Here $`A(t)`$ is time-dependent amplitude of a cluster, $`\mathrm{}_{\text{cl}}`$ is the characteristic width of the cluster. Shape function $`\theta (x)`$ may be chosen with a maximum equal one at $`x=0`$ and unit width. Real shapes of various clusters in the turbulent ensemble are determined by a competition of different terms in the evolution equation (II.1). However, we believe that particular shapes affect only numerical factors in the expression for the growth rate of clusters and do not effect their functional dependence on the parameters of the problem which is considered in this subsection.

#### II.2.1 Effect of turbulent diffusion

The advective term in the LHS of Eq. (II.1) results in turbulent diffusion inside the cluster. This effect may be modelled by renormalization of the molecular diffusion coefficient $`D`$ in the right hand side (RHS) of Eq. (II.1) by the effective turbulent diffusion coefficient $`D_\text{T}`$ with a usual estimate of $`D_\text{T}`$:

$`DD+D_\text{T},D_\text{T}\mathrm{}_{\text{cl}}v_{\text{cl}}/3.`$ (12)
Hereafter $`v_{\text{cl}}`$ is the mean square velocity of particles at the scale $`\mathrm{}_{\text{cl}}`$. Instead of the full Eq. (II.1) consider now a model equation

$`{\displaystyle \frac{\stackrel{~}{\mathrm{\Theta }}(𝒓_0|t,𝒓)}{t}}=D_\text{T}\mathrm{\Delta }\stackrel{~}{\mathrm{\Theta }}(𝒓_0|t,𝒓),`$ (13)
which accounts only for turbulent diffusion. Equations (10) and (13) yield:

$`_q(t)/tq|\stackrel{~}{\mathrm{\Theta }}(𝒓_0|t,𝒓)|^{q1}D_\text{T}\mathrm{\Delta }|\stackrel{~}{\mathrm{\Theta }}(𝒓_0|t,𝒓)|_w.`$ (14)
Substituting distribution (11) we estimate Laplacian in Eq. (14) as $`1/\mathrm{}_{\text{cl}}^2`$. Equations (13) and (14) imply that

$`_q(t)/t=qD_\text{T}_q(t)/\mathrm{}_{\text{cl}}^2.`$ (15)
The solution of Eq. (15) reads:

$`_q(t)`$ $`=`$ $`_q(0)\mathrm{exp}[\gamma _{\text{dif}}(q)t],`$
$`\gamma _{\mathrm{dif}}(q)`$ $``$ $`qD_\text{T}/\mathrm{}_{\text{cl}}^2,`$ (16)
where $`\gamma _{\text{dif}}(q)`$ denotes a contribution to the damping rate of $`_q(t)`$ caused by turbulent diffusion.

#### II.2.2 Effect of particles inertia

In this subsection we show that the term $`\stackrel{~}{\mathrm{\Theta }}\text{div}𝑾`$ in the RHS of Eq. (II.1) can result in an exponential growth of $`_q(t)\mathrm{exp}[\gamma _{\text{in}}(q)t],`$ i.e., in the instability. We denoted here the contribution to the growth rate of $`_q(t)`$, caused by the inertia of particles, by $`\gamma _{\text{in}}(q)`$. In order to evaluate $`\gamma _{\text{in}}(q)`$ we neglect now in Eq. (II.1) both, the convective term in the LHS of this equation (i.e. the the turbulent velocity difference inside the cluster) and the molecular diffusion term. The resulting equation reads:

$`{\displaystyle \frac{\stackrel{~}{\mathrm{\Theta }}(𝒓_0|t,𝒓)}{t}}=\stackrel{~}{\mathrm{\Theta }}(𝒓_0|t,𝒓)\text{div}𝑾(𝒓_0|t,𝒓).`$ (17)
The main contribution to the BL-velocity difference $`𝑾(𝒓_0|t,𝒓)`$ in the RHS of this equation is due to the eddies with size $`\mathrm{}_{\text{cl}}`$, the characteristic size of the cluster. Denote by $`v_{\text{cl}}`$ the characteristic velocity of these eddies and by $`\tau _v\mathrm{}_{\text{cl}}/v_{\text{cl}}`$ the corresponding correlation time. In our qualitative analysis we neglect the $`𝒓`$ dependence of $`\text{div}𝑾(𝒓_0|t,𝒓)`$ inside the cluster and consider the divergence in Eq. (17) as a random process $`b(t)`$ with a correlation time $`\tau _v`$:

$`\text{div}𝑾(𝒓_0|t,𝒓)b(t).`$ (18)
Together with the decomposition (11) this yields the following equation for the cluster amplitude $`A(t)`$:

$`{\displaystyle \frac{A(t)}{t}}=A(t)b(t).`$ (19)
The solution of Eq. (19) reads:

$`A(t)=A_0\mathrm{exp}[I(t)],I(t){\displaystyle _0^t}b(\tau )𝑑\tau .`$ (20)
Integral $`I(t)`$ in Eq. (20) can be rewritten as a sum of integrals $`I_n`$ over small time intervals $`\tau _v`$:

$$I(t)=\underset{n=1}{\overset{t/\tau _v}{}}I_n,I_n(t)\underset{(n1)\tau _v}{\overset{n\tau _v}{}}b(\tau )𝑑\tau .$$

In our qualitative analysis integrals $`I_n`$ may be considered as independent random variables. Using the central limit theorem we estimate the total integral

$$I(t)\sqrt{I_n^2_v}\sqrt{N}\zeta ,I_n^2_v=b^2_v\tau _v^2,$$

where $`\mathrm{}_v`$ denotes averaging over turbulent velocity ensemble, $`\zeta `$ is a Gaussian random variable with zero mean and unit variance, $`N=t/\tau _v.`$ Now we calculate

$$_q(t)=\mathrm{\Theta }^qP(\zeta )𝑑\zeta ,P(\zeta )=(1/\sqrt{2\pi })\mathrm{exp}(\zeta ^2/2).$$

Therefore, $`_q(t)=J_1\mathrm{exp}(q^2S^2N/2),`$ where

$$J_1=(1/\sqrt{2\pi })\mathrm{exp}[(\zeta qS\sqrt{N})^2/2]𝑑\zeta 1.$$

Since the main contribution to the integral $`J_1`$ arises from $`\zeta qS\sqrt{N}`$, the parameter $`q`$ cannot be large. In this approximation the $`q`$moments

$$_q(t)=_q(0)\mathrm{exp}[\gamma _{\text{in}}(q)t]$$

with $`\gamma _{\text{in}}(q)`$ being the growth rate of the $`q`$-th moment due to particles inertia which is given by

$`\gamma _{\text{in}}(q){\displaystyle \frac{1}{2}}\tau _v[\text{div}𝑾(𝒓_0|t,𝒓)]^2_vq^2.`$ (21)
#### II.2.3 Qualitative picture of the clustering instability

In previous subsections we evaluated the contributions to the growth rate of $`_q(t)`$ due to the turbulent diffusion $`\gamma _{\text{dif}}(q)`$, Eq. (16), and due to the particles inertia $`\gamma _{\text{in}}(q)`$, Eq. (21). The total growth rate may be evaluated as a sum of these contributions:

$`_q(t)`$ $`=`$ $`_q(0)\mathrm{exp}(\gamma _qt),`$ (22)
$`\gamma _q`$ $``$ $`\gamma _{\text{dif}}(q)+\gamma _{\text{in}}(q),`$
$`\gamma _q`$ $``$ $`qD_\text{T}/\mathrm{}_{\text{cl}}+{\displaystyle \frac{1}{2}}\tau _v[\text{div}𝑾(𝒓_0|t,𝒓)]^2_vq^2.`$
Clearly, the instability is caused by a nonzero value of $`\tau _v(\text{div}𝑾)^2,`$ i.e., by a compressibility of the particle velocity field $`𝒗(t,𝒓)`$.

Compressibility of fluid velocity itself $`𝒖(t,𝒓)`$ (including atmospheric turbulence) is often negligible, i.e., $`\text{div}𝒖0`$. However, due to the effect of particles inertia their velocity $`𝒗(t,𝒓)`$ does not coincide with $`𝒖(t,𝒓)`$ (see, e.g., M87 ; MC86 ; WM93 ; MCW96 ), and a degree of compressibility, $`\sigma _v`$, of the field $`𝒗(t,𝒓)`$, may be of the order of unityEKRI96 ; EKR96 ; EKRS00 .

Parameter $`\sigma _v`$ is defined as

$`\sigma _v[\text{div}𝒗]^2/|\mathbf{}\times 𝒗|^2.`$ (23)
Note that parameter $`\sigma _v`$ is independent of the scale of the turbulent velocity field. It characterizes a compressible part of the velocity field as a whole. The main contribution to this parameter comes from the scales which are of the order of the Kolmogorov scale $`\eta `$.

For inertial particles $`\text{div}𝒗\tau _\text{p}\mathrm{\Delta }P/\rho `$, where $`\tau _\text{p}`$ is the particle response time,

$`\tau _\text{p}=m_\text{p}/6\pi \rho \nu a=2\rho _\text{p}a^2/9\rho \nu ,`$ (24)
where $`m_\text{p}`$ and $`\rho _\text{p}`$ are the mass and material density of particles, respectively. The fluid flow parameters are: pressure $`P`$, Reynolds number $`e=Lu_\text{T}/\nu `$, the dissipative scale of turbulence $`\eta =Le^{3/4}`$, the maximum scale of turbulent motions $`L`$ and the turbulent velocity $`u_\text{T}`$ in the scale $`L`$. Now we can estimate $`\sigma _v`$ as

$`\sigma _v(\rho _\text{p}/\rho )^2(a/\eta )^4(a/a_{})^4`$ (25)
(see EKRS00 ), where $`a_{}`$ is a characteristic radius of particles. For $`a>a_{}`$ it is plausible to correct this estimate as follows:

$`\sigma _v{\displaystyle \frac{a^4}{a^4+a_{}^4}}.`$ (26)
For water droplets in the atmosphere $`\rho _\text{p}/\rho 10^3`$ and $`a_{}\eta /30`$. For the typical value of $`\eta 1`$mm it yields $`a_{}30\mu `$m. On windy days when $`\eta `$ decreases, the value of $`a_{}`$ correspondingly becomes smaller.

Then we estimate $`\tau _v[\text{div}𝒗]^2`$ as $`2\sigma _v/\tau _v`$, because $`[\text{rot}𝒗]^22\tau _v^2`$. Assuming that the cluster size $`\mathrm{}_{\text{cl}}`$ is of the order of the inner scale of turbulence, $`\eta `$, we have to identify $`\tau _v`$ with a turnover time of eddies in the inner scale $`\eta `$, $`\tau _v\tau _\eta \eta /v_\eta =(L/u_\text{T})e^{1/2}`$. Thus, the growth rate $`\gamma _q`$ of the $`q`$-th moment in Eq. (22) may be evaluated as

$`\gamma _q\gamma _{\text{cl}}q(qq_{\text{cr}}),\gamma _{\text{cl}}{\displaystyle \frac{\sigma _v}{\tau _\eta }},q_{\text{cr}}{\displaystyle \frac{1}{3\sigma _v}}.`$ (27)
Clearly, the moments with $`q>q_{\text{cr}}`$ are unstable. Equations (25) and (27) imply that it happens when $`a>a_{q,\mathrm{cr}}`$ where $`a_{q,\mathrm{cr}}=a_{1,\mathrm{cr}}/q^{1/4}`$ is the value of $`a`$ at which $`q_{\text{cr}}=q`$. The largest value of $`a_{q,\mathrm{cr}}`$ corresponds to the instability of the first moment, $`|\mathrm{\Theta }|`$: $`a_{1,\mathrm{cr}}0.8a_{}`$, $`a_{2,\mathrm{cr}}0.84a_{1,\mathrm{cr}}`$, $`a_{3,\mathrm{cr}}0.76a_{1,\mathrm{cr}}`$, $`a_{4\mathrm{c}\mathrm{r}}0.71a_{1,\mathrm{cr}}`$, etc.

Note that if $`|\mathrm{\Theta }|`$ grows in time then almost all particles can be accumulated inside the clusters (if we neglect a nonlinear saturation of such growth). We define this case as a strong clustering. On the other hand, if $`q_{\text{cr}}>1`$ the first moment $`|\mathrm{\Theta }|`$ does not grow, and the clusters contain a small fraction of the total number of particles. This does not mean that the instability of higher moments is not important. Thus, e.g., the rate of binary particles collisions is proportional to the square of their number density $`n^2=(\overline{n})^2+|\mathrm{\Theta }|^2`$. Therefore, the growth of the 2nd moment, $`|\mathrm{\Theta }|^2`$, (which we define as a weak clustering) results in that binary collisions occur mainly between particles inside the cluster. The latter can be important in coagulation of droplets in atmospheric clouds whereby the collisions between droplets play a crucial role in a rain formation. The growth of the $`q`$-th moment, $`|\mathrm{\Theta }|^q`$, results in that $`q`$-particles collisions occur mainly between particles inside the cluster. The growth of the negative moments of particles number density (possibly associated with formation of voids and cellular structures) was discussed in BF01  (see also SZ89 ; KS97 ).

In the above qualitative analysis whereby we considered only one-point correlation functions of the number density of particles, we missed an important effect of an effective drift velocity which decreases a growth rate of the clustering instability. For the one-point correlation functions of the number density of particles the effective drift velocity is zero for homogeneous and isotropic turbulence. However, in the equations for two-point and multi-point correlation functions of the number density of particles the effective drift velocity is not zero and as we will see in the next section it increases a threshold for the clustering instability.

## III The clustering instability of the 2nd moment

### III.1 Basic equations

In the previous section we estimated the growth rates of all moments $`|\mathrm{\Theta }|^q`$. Here we present the results of a rigorous analysis of the evolution of the two-point 2nd moment

$`\mathrm{\Phi }(t,𝑹)\mathrm{\Theta }(t,𝒓)\mathrm{\Theta }(t,𝒓+𝑹).`$ (28)
In this analysis we used stochastic calculus \[e.g., Wiener path integral representation of the solution of the Cauchy problem for Eq. (1), Feynman-Kac formula and Cameron-Martin-Girsanov theorem\]. The comprehensive description of this approach can be found in EKRS99 ; ZRS90 ; ZMR88 ; EKR00 .

We showed that a finite correlation time of a turbulent velocity plays a crucial role for the clustering instability. Notably, an equation for the second moment $`\mathrm{\Phi }(t,𝑹)`$ of the number density of inertial particles comprises spatial derivatives of high orders due to a non-local nature of turbulent transport of inertial particles in a random velocity field with a finite correlation time (see Appendix A and EKRS00 ). However, we found that equation for $`\mathrm{\Phi }(t,𝑹)`$ is a second-order partial differential equation at least for two models of a random velocity field:

Model I. The random velocity with Gaussian statistics of the integrals $`_0^t𝒗(t^{},𝝃)𝑑t^{}`$ and $`_0^tb(t^{},𝝃)𝑑t^{}`$, see Appendix B.

Model II. The Gaussian velocity field with a small yet finite correlation time, see Appendix C.

In both models equation for $`\mathrm{\Phi }(t,𝑹)`$ has the same form:

$`\mathrm{\Phi }/t=\widehat{}\mathrm{\Phi }(t,𝑹),`$ (29)
$`\widehat{}=B(𝑹)+2𝑼(𝑹)\mathbf{}+\widehat{D}_{\alpha \beta }(𝑹)_\alpha _\beta ,`$
but with different expressions for its coefficients. The meaning of the coefficients $`B(𝑹)`$, $`𝑼(𝑹)`$ and $`\widehat{D}_{\alpha \beta }(𝑹)`$ is as follows:

– Function $`B(𝑹)`$ is determined only by a compressibility of the velocity field and it causes generation of fluctuations of the number density of inertial particles.

– The vector $`𝑼(𝑹)`$ determines a scale-dependent drift velocity which describes a transfer of fluctuations of the number density of inertial particles over the spectrum. Note that $`𝐔(𝑹=0)=0`$ whereas $`B(𝑹=0)0.`$ For incompressible velocity field $`𝐔(𝑹)=0,`$ $`B(𝑹)=0`$.

– The scale-dependent tensor of turbulent diffusion $`\widehat{D}_{\alpha \beta }(𝑹)`$ is also affected by the compressibility.

In very small scales this tensor is equal to the tensor of the molecular (Brownian) diffusion, while in the vicinity of the maximum scale of turbulent motions this tensor coincides with the usual tensor of turbulent diffusion. Tensor $`\widehat{D}_{\alpha \beta }(𝑹)`$ may be written as

$`\widehat{D}_{\alpha \beta }(𝑹)`$ $`=`$ $`2D\delta _{\alpha \beta }+D_{\alpha \beta }^\text{T}(𝑹),`$ (30)
$`D_{\alpha \beta }^\text{T}(𝑹)`$ $`=`$ $`\stackrel{~}{D}_{\alpha \beta }^\text{T}(0)\stackrel{~}{D}_{\alpha \beta }^\text{T}(𝑹).`$
In Appendix B we found that for Model I:

$`B(𝑹)`$ $``$ $`2{\displaystyle _0^{\mathrm{}}}b[0,𝝃(𝒓_1|0)]b[\tau ,𝝃(𝒓_2|\tau )]𝑑\tau ,`$ (31)
$`𝑼(𝑹)`$ $``$ $`2{\displaystyle _0^{\mathrm{}}}𝒗[0,𝝃(𝒓_1|0)]b[\tau ,𝝃(𝒓_2|\tau )]𝑑\tau ,`$
$`\stackrel{~}{D}_{\alpha \beta }^\text{T}(𝑹)`$ $``$ $`2{\displaystyle _0^{\mathrm{}}}v_\alpha [0,𝝃(𝒓_1|0)]v_\beta [\tau ,𝝃(𝒓_2|\tau )]𝑑\tau .`$
For the $`\delta `$-correlated in time random Gaussian compressible velocity field the operator $`\widehat{}`$ is replaced by $`\widehat{}_0`$ in the equation for the second moment $`\mathrm{\Phi }(t,𝑹)`$, where

$`\widehat{}_0B_0(𝑹)+2𝑼_0(𝑹)\mathbf{}+\widehat{D}_{\alpha \beta }(𝑹)_\alpha _\beta ,`$
$`B_0(𝑹)=_\alpha _\beta \widehat{D}_{\alpha \beta }(𝑹),`$ (32)
$`U_{0,\alpha }(𝑹)=_\beta \widehat{D}_{\alpha \beta }(𝑹)`$
(for details see EKRS00 ; EKRS99 ). In the $`\delta `$-correlated in time velocity field the second moment $`\mathrm{\Phi }(t,𝑹)`$ can only decay in spite of the compressibility of the velocity field. The reason is that the differential operator $`\widehat{}_0_\alpha _\beta \widehat{D}_{\alpha \beta }(𝑹)`$ is adjoint to the operator $`\widehat{}_0^{}\widehat{D}_{\alpha \beta }(𝑹)_\alpha _\beta `$ and their eigenvalues are equal. The damping rate for the equation

$`\mathrm{\Phi }/t=\widehat{}_0^{}\mathrm{\Phi }(t,𝑹)`$ (33)
has been found in Ref. EKR95  for a compressible isotropic homogeneous turbulence in a dissipative range:

$`\gamma _2={\displaystyle \frac{(3\sigma _\text{T})^2}{6\tau _\eta (1+\sigma _\text{T})(1+3\sigma _\text{T})}}.`$ (34)
Here $`\sigma _\text{T}`$ is the degree of compressibility of the tensor $`D_{\alpha \beta }^\text{T}(𝑹)`$. For the $`\delta `$-correlated in time incompressible velocity field $`(\sigma _\text{T}=0)`$ Eq. (33) was derived in Ref. K68 . Thus, for the Kraichnan model of turbulent advection (with a delta correlated in time velocity field) the clustering instability of the 2nd moment does not occur.

A general form of the turbulent diffusion tensor in a dissipative range is given by

$`D_{\alpha \beta }^^\mathrm{T}(𝑹)=(C_1R^2\delta _{\alpha \beta }+C_2R_\alpha R_\beta )/\tau _\eta ,`$ (35)
$`C_1=2(2+\sigma _\text{T})/3(1+\sigma _\text{T}),`$
$`C_2=2(2\sigma _\text{T}1)/3(1+\sigma _\text{T}).`$
The parameter $`\sigma _\text{T}`$ is defined by analogy with Eq. (23):

$`\sigma _\text{T}{\displaystyle \frac{\mathbf{}𝑫_\text{T}\mathbf{}}{\mathbf{}\times 𝑫_\text{T}\times \mathbf{}}}={\displaystyle \frac{_\alpha _\beta D_{\alpha \beta }^^\mathrm{T}(𝑹)}{_\alpha _\beta D_{\alpha ^{}\beta ^{}}^^\mathrm{T}(𝑹)ϵ_{\alpha \alpha ^{}\gamma }ϵ_{\beta \beta ^{}\gamma }}},`$ (36)
where $`ϵ_{\alpha \beta \gamma }`$ is the fully antisymmetric unit tensor. Equations (23) and (36) imply that $`\sigma _\text{T}=\sigma _v`$ in the case of $`\delta `$-correlated in time compressible velocity field. Equations (31) show that for a finite correlation time identities (32) are violated and

$$B(𝑹)B_0(𝑹),𝑼(𝑹)𝑼_0(𝑹).$$

For a random incompressible velocity field with a finite correlation time the tensor of turbulent diffusion $`D_{\alpha \beta }(𝑹)=\tau ^1\xi _\alpha (𝒓_1)\xi _\beta (𝒓_2)`$ \[see Eq. (78)\] and the degree of compressibility of this tensor is

$`\sigma _\text{T}={\displaystyle \frac{(\mathbf{}𝝃)^2}{(\mathbf{}\times 𝝃)^2}},`$ (37)
where $`𝝃(𝒓_1|t)`$ is the Lagrangian displacement of a particle trajectory which paths through point $`𝒓_1`$ at $`t=0.`$ Remarkably that Taylor T21  obtained the coefficient of turbulent diffusion for the mean field in the form $`D_\text{T}(𝑹=0)=\tau ^1\xi _\alpha (𝒓_1|t)\xi _\alpha (𝒓_1|t)`$.

### III.2 Clustering instability in Model I

Let us study the clustering instability for the model of the random velocity with Gaussian statistics of the integrals

$$_0^t𝒗(t^{},𝝃)𝑑t^{},_0^tb(t^{},𝝃)𝑑t^{},$$

see Appendix B. In this model Eq. (29) in a non-dimensional form reads:

$`{\displaystyle \frac{\mathrm{\Phi }}{\stackrel{~}{t}}}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Phi }^{\prime \prime }}{m(r)}}+\left[{\displaystyle \frac{1}{m(r)}}+(UC_2)r^2\right]{\displaystyle \frac{2\mathrm{\Phi }^{}}{r}}+B\mathrm{\Phi },`$
$`1/m(r)`$ $``$ $`(C_1+C_2)r^2+2/\mathrm{Sc},`$ (38)
where $`𝑼U𝑹`$ and Sc$`=\nu /D`$ is the Schmidt number. For small inertial particles advected by air flow Sc$`1`$. The non-dimensional variables in Eq. (38) are $`rR/\eta `$ and $`\stackrel{~}{t}=t/\tau _\eta `$, $`B`$ and $`U`$ are measured in the units $`\tau _\eta ^1`$. Consider a solution of Eq. (38) in two spatial regions.

a. Molecular diffusion region of scales. In this region $`r\text{Sc}^{1/2}`$, and all terms $`r^2`$ (with $`C_1,C_2`$ and $`U`$) may be neglected. Then the solution of Eq. (38) is given by

$`\mathrm{\Phi }(r)=(1\alpha r^2)\mathrm{exp}(\gamma _2t),`$
where

$`\alpha `$ $`=`$ $`\text{Sc}(B\gamma _2\tau _\eta )/12,B>\gamma _2\tau _\eta .`$
b. Turbulent diffusion region of scales. In this region $`\text{Sc}^{1/2}r1`$, the molecular diffusion term $`1/\text{Sc}`$ is negligible. Thus, the solution of (38) in this region is

$`\mathrm{\Phi }(r)=A_1r^\lambda \mathrm{exp}(\gamma _2t),`$
where

$`\lambda `$ $`=`$ $`(C_1C_2+2U\pm iC_3)/2(C_1+C_2),`$
$`C_3^2`$ $`=`$ $`4(B\gamma _2\tau _\eta )(C_1+C_2)(C_1C_2+2U)^2.`$
Since the total number of particles in a closed volume is conserved:

$$_0^{\mathrm{}}r^2\mathrm{\Phi }(r)𝑑r=0.$$

This implies that $`C_3^2>0`$, and therefore $`\lambda `$ is a complex number. Since the correlation function $`\mathrm{\Phi }(r)`$ has a global maximum at $`r=0,`$ $`C_1>C_22U`$. The latter condition for very small $`U`$ yields $`\sigma _T3`$. For $`r1`$ the solution for $`\mathrm{\Phi }(r)`$ decays sharply with $`r`$. The growth rate $`\gamma _2`$ of the second moment of particles number density can be obtained by matching the correlation function $`\mathrm{\Phi }(r)`$ and its first derivative $`\mathrm{\Phi }^{}(r)`$ at the boundaries of the above regions, i.e., at the points $`r=\text{Sc}^{1/2}`$ and $`r=1`$. The matching yields $`C_3/2(C_1+C_2)2\pi /\mathrm{ln}\text{Sc}`$. Thus,

$`\gamma _2`$ $`=`$ $`{\displaystyle \frac{1}{\tau _\eta (1+3\sigma _\text{T})}}[{\displaystyle \frac{200\sigma _\text{U}(\sigma _\text{T}\sigma _\text{U})}{3(1+\sigma _\text{U})}}{\displaystyle \frac{(3\sigma _\text{T})^2}{6(1+\sigma _\text{T})}}`$ (39)
$`{\displaystyle \frac{3\pi ^2(1+3\sigma _\text{T})^2}{(1+\sigma _\text{T})\mathrm{ln}^2\text{Sc}}}]+{\displaystyle \frac{20(\sigma _\text{B}\sigma _\text{U})}{\tau _\eta (1+\sigma _\text{B})(1+\sigma _\text{U})}},`$
where we introduced parameters $`\sigma _\text{B}`$ and $`\sigma _\text{U}`$ defined by

$$B=20\sigma _\text{B}/(1+\sigma _\text{B}),U=20\sigma _\text{U}/3(1+\sigma _\text{U}).$$
(40)
Note that the parameters $`\sigma _\text{B}\sigma _\text{U}\sigma _v`$. For the $`\delta `$-correlated in time random compressible velocity field $`\sigma _\text{B}=\sigma _\text{U}=\sigma _\text{T}=\sigma _v`$. Figure 1 shows the range of parameters $`(\sigma _v,\sigma _\text{T})`$ for $`\sigma _\text{B}=\sigma _\text{U}=\sigma _v`$ in the case of $`\text{Sc}=10^3`$ (curve c), $`\text{Sc}=10^5`$ (curve b) and $`\text{Sc}\mathrm{}`$ (curve a). The dashed line $`\sigma _v=\sigma _\text{T}`$ corresponds to the $`\delta `$-correlated in time random compressible velocity field. This is a limiting line for the curve ”a”. Figure 1 demonstrates that even a very small deviations from the $`\delta `$-correlated in time random compressible velocity field results in the instability of the second moment of the number density of inertial particles. The minimum value of $`\sigma _\text{T}`$ required for the clustering instability is $`\sigma _\text{T}0.26`$ and a corresponding value of $`\sigma _v0.12`$ (see Fig. 1). For smaller value of $`\sigma _v`$ the clustering instability can occur, but it requires larger values of $`\sigma _\text{T}.`$

Notably, in Model II of a random velocity field (i.e., the Gaussian velocity field with a small yet finite correlation time) the clustering instability occurs when $`\sigma _v>0.2`$ (see Appendix C). Indeed, the growth rate $`\gamma _2`$ of the second moment of particles number density is determined by equation:

$`\mathrm{\Gamma }`$ $`=`$ $`\stackrel{~}{B}(\sigma _v)\mathrm{St}^2{\displaystyle \frac{(3\sigma _v)^2}{6(1+\sigma _v)(1+3\sigma _v)}}`$
$`{\displaystyle \frac{8(1+3\sigma _v)}{3(1+\sigma _v)}}\left({\displaystyle \frac{\pi }{\mathrm{ln}\mathrm{Sc}}}\right)^2,`$
$`\stackrel{~}{B}(\sigma _v)`$ $`=`$ $`12\left(b_2+{\displaystyle \frac{b_3a_1}{4a_2^2}}{\displaystyle \frac{b_1}{2a_2}}\right),`$
where St$`=\overline{\tau }_{\text{ren}}/\tau _\eta `$ is the Strouhal number, $`\mathrm{\Gamma }=\gamma _2(1+\overline{\tau }_{\text{ren}}\gamma _2)^2,`$ and

$`a_1`$ $`=`$ $`{\displaystyle \frac{2(19\sigma _v+3)}{3(1+\sigma _v)}},a_2={\displaystyle \frac{2(3\sigma _v+1)}{3(1+\sigma _v)}},`$
$`b_1`$ $`=`$ $`{\displaystyle \frac{1}{27(1+\sigma _v)^2}}(121278\sigma _v3067\sigma _v^2),`$
$`b_2`$ $`=`$ $`{\displaystyle \frac{850}{9}}\left({\displaystyle \frac{\sigma _v}{1+\sigma _v}}\right)^2,`$
$`b_3`$ $`=`$ $`{\displaystyle \frac{1}{27(1+\sigma _v)^2}}(36+466\sigma _v+2499\sigma _v^2).`$
Figure 2 shows the range of parameters $`(\text{Sc},\sigma _v)`$ for $`\text{Sc}=10^3`$ (curve c), $`\text{Sc}=10^5`$ (curve b) and $`\text{Sc}\mathrm{}`$ (curve a).

It is seen from FIG. 2 that for $`\sigma _v>0.2`$ the second-order correlation function of the number density of inertial particles can grow in time exponentially (i.e., $`\gamma _2>0)`$ even for very small Strouhal numbers. For example, in the vicinity of $`\sigma _v=3,`$ the growth rate $`\gamma _2`$ of the clustering instability of the second-order correlation function is given by

$$\gamma _2=4\times 10^3\left(\frac{\tau u_\eta }{\eta }\right)^2\frac{(3\sigma _v)^2}{240}\frac{20}{3}\left(\frac{\pi }{\mathrm{ln}\mathrm{Sc}}\right)^2.$$
(41)
The sufficient condition for the exponential growth of the second moment of a number density of inertial particles is $`\mathrm{Sc}>\mathrm{Sc}^{(\mathrm{cr})},`$ where the critical Schmidt number $`\mathrm{Sc}^{(\mathrm{cr})}`$ is given by $`\mathrm{Sc}^{(\mathrm{cr})}=\mathrm{Sc}(\gamma _2=0).`$ The clustering instability occurs when the degree of compressibility of particles velocity $`\sigma _v>0.2,`$ i.e., for particles and droplets with the radius $`a_{}>25.4\mu `$m . Equation (27) also yields a similar value $`\sigma _{\mathrm{cr}}1/6`$ for the threshold of the instability of the 2nd moment (at $`q_{\mathrm{cr}}=2`$). Note that Eq. (34) is written for $`\text{Sc}\mathrm{}`$.

## IV Nonlinear effects

The compressibility of the turbulent velocity field with a finite correlation time can cause the exponential growth of the moments of particles number density. This small-scale instability results in formation of strong inhomogeneities (clusters) in particles spatial distributions. The linear analysis does not allow to determine a mechanism of saturation of the clustering instability. As can be seen from Eq. (39) molecular diffusion only depletes the growth rates of the clustering instability at the linear stage (contrary to the instability discussed in Ref. BF01 ). The clustering instability is saturated by nonlinear effects.

Now let us discuss a mechanism of the nonlinear saturation of the clustering instability using on the example of atmospheric turbulence with characteristic parameters: $`\eta 1`$mm, $`\tau _\eta (0.10.01)`$s. A momentum coupling of particles and turbulent fluid is essential when $`m_\mathrm{p}n_{\text{cl}}\rho `$, i.e., the mass loading parameter $`\varphi =m_\mathrm{p}n_{\text{cl}}/\rho `$ is of the order of unity (see, e.g., CST98 ). This condition implies that the kinetic energy of fluid $`\rho 𝒖^2`$ is of the order of the particles kinetic energy $`m_\mathrm{p}n_{\text{cl}}𝒗^2,`$ where $`|𝒖||𝒗|`$. This yields:

$`n_{\text{cl}}a^3(\rho /3\rho _\mathrm{p}).`$ (42)
For water droplets $`\rho _\mathrm{p}/\rho 10^3`$. Thus, for $`a=a_{}30\mu `$m we obtain $`n_{\text{cl}}10^4`$ cm<sup>-3</sup> and the total number of particles in the cluster of size $`\eta `$, $`N_{\mathrm{cl}}\eta ^3n_{\text{cl}}10`$. This values may be considered as a lower estimate for the “two-way coupling” when the effect of fluid on particles has to be considered together with the feed-back effect of the particles on the carrier fluid. However, it was found in 02LP  that turbulence modification by particles is governed by the ratio of the particle energy and  the total energy of the suspension (rather then the energy of the carrier fluid) and thus by parameter $`\varphi (1+\varphi )`$ (rather then by $`\varphi `$ itself). Thus we expect that the two-way coupling can only mitigate but not stop the clustering instability.

An actual mechanism of the nonlinear saturation of the clustering instability is “four way coupling” when the particle-particle interaction is also important. In this situation the particles collisions result in effective particle pressure which prevents further grows of concentration. Particles collisions play essential role when during the life-time of a cluster the total number of collisions is of the order of number of particles in the cluster. The rate of collisions $`Jn_{\text{cl}}/\tau _\eta `$ can be estimated as $`J4\pi a^2n_{\text{cl}}^2|𝒗_{\mathrm{rel}}|`$. The relative velocity $`𝒗_{\mathrm{rel}}`$ of colliding particles with different but comparable sizes can be estimated as $`|𝒗_{\mathrm{rel}}|\tau _\mathrm{p}|(𝒖\mathbf{})𝒖|\tau _\mathrm{p}u_\eta ^2/\eta `$. Thus the collisions in clusters may be essential for

$`n_{\text{cl}}a^3(\eta /a)(\rho /3\rho _\mathrm{p}),\mathrm{}_\mathrm{s}a(3a\rho _\mathrm{p}/\eta \rho )^{1/3},`$ (43)
where $`\mathrm{}_\mathrm{s}`$ is a mean separation of particles in the cluster. For the above parameters ($`a=30\mu `$m) $`n_{\text{cl}}3\times 10^5`$cm<sup>-3</sup>, $`\mathrm{}_\mathrm{s}5a150\mu `$m and $`N_{\text{cl}}300`$. Note that the mean number density of droplets in clouds $`\overline{n}`$ is about $`10^210^3`$cm<sup>-3</sup>. Therefore the clustering instability of droplets in the clouds increases their concentrations in the clusters by the orders of magnitude.

In all our analysis we have neglected the effect of sedimentation of particles in gravity field which is essential for particles of the radius $`a>100\mu `$m. Taking $`\mathrm{}_{\text{cl}}\eta `$ we assumed implicitly that $`\tau _\mathrm{p}<\tau _\eta `$. This is valid (for the atmospheric conditions) if $`a60\mu `$m. Otherwise the cluster size can be estimated as $`\mathrm{}_{\text{cl}}\eta (\tau _\mathrm{p}/\tau _\eta )^{3/2}`$.

Our estimates support the conjecture that the clustering instability serves as a preliminary stage for a coagulation of water droplets in clouds leading to a rain formation.

## V Discussion

In this study we investigated the clustering instability of the spatial distribution of inertial particles advected by a turbulent velocity field. The instability results in formation of clusters, i.e., small-scale inhomogeneities of aerosols and droplets. The clustering instability is caused by a combined effect of the particle inertia and finite correlation time of the velocity field. The finite correlation time of the turbulent velocity field causes the compressibility of the field of Lagrangian trajectories. The latter implies that the number of particles flowing into a small control volume in a Lagrangian frame does not equal to the number of particles flowing out from this control volume during a correlation time. This can result in the depletion of turbulent diffusion.

The role of the compressibility of the velocity field is as follows. Divergence of the velocity field of the inertial particles $`\text{div}𝒗=\tau _p\mathrm{\Delta }P/\rho .`$ The inertia of particles results in that particles inside the turbulent eddies are carried out to the boundary regions between the eddies by inertial forces (i.e., regions with low vorticity and high strain rate). For a small molecular diffusivity $`\text{div}𝒗dn/dt`$ \[see Eq. (1)\]. Therefore, $`dn/dt\tau _p\mathrm{\Delta }P/\rho .`$ Thus there is accumulation of inertial particles (i.e., $`dn/dt>0)`$ in regions with $`\mathrm{\Delta }P<0.`$ Similarly, there is an outflow of inertial particles from the regions with $`\mathrm{\Delta }P>0.`$ This mechanism acts in a wide range of scales of a turbulent fluid flow. Turbulent diffusion results in relaxation of fluctuations of particles concentration in large scales. However, in small scales where turbulent diffusion is small, the relaxation of fluctuations of particle concentration is very weak. Therefore the fluctuations of particle concentration are localized in the small scales.

This phenomenon is considered for the case when density of fluid is much less than the material density $`\rho _p`$ of particles $`(\rho \rho _p).`$ When $`\rho \rho _p`$ the results coincide with those obtained for the case $`\rho \rho _p`$ except for the transformation $`\tau _p\beta _{}\tau _p,`$ where

$`\beta _{}=2\left(1+{\displaystyle \frac{\rho }{\rho _p}}\right)\left({\displaystyle \frac{\rho _p\rho }{2\rho _p+\rho }}\right).`$
For $`\rho \rho _p`$ the value $`dn/dt\beta _{}\tau _p\mathrm{\Delta }P/\rho .`$ Thus there is accumulation of inertial particles (i.e., $`dn/dt>0)`$ in regions with the minimum pressure of a turbulent fluid since $`\beta _{}<0.`$ In the case $`\rho \rho _p`$ we used the equation of motion of particles in fluid flow which takes into account contributions due to the pressure gradient in the fluid surrounding the particle (caused by acceleration of the fluid) and the virtual (”added”) mass of the particles relative to the ambient fluid MR83 .

The exponential growth of the second moment of a number density of inertial particles due to the small-scale instability can be saturated by the nonlinear effects (see Section IV). The excitation of the second moment of a number density of particles requires two kinds of compressibilities: compressibility of the velocity field and compressibility of the field of Lagrangian trajectories, which is caused by a finite correlation time of a random velocity field. Remarkably, the compressibility of the field of Lagrangian trajectories determines the coefficient of turbulent diffusion (i.e. the coefficient $`D_{\alpha \beta }`$ near the second-order spatial derivative of the second moment of a number density of inertial particles in Eq. (29). The compressibility of the field of Lagrangian trajectories causes depletion of turbulent diffusion in small scales even for $`\sigma _v=0`$. On the other hand, the compressibility of the velocity field determines a coefficient B(r) near the second moment of a number density of inertial particles in Eq. (29). This term is responsible for the exponential growth of the second moment of a number density of particles.

Summary:

$``$ We showed that the physical reason for the clustering instability in spatial distribution of particles in turbulent flows is a combined effect of the inertia of particles leading to a compressibility of the particle velocity field $`𝒗(t,𝒓)`$ and a finite velocity correlation time.

$``$ The clustering instability can result in a strong clustering whereby a finite fraction of particles is accumulated in the clusters and a weak clustering when a finite fraction of particle collisions occurs in the clusters.

$``$ The crucial parameter for the clustering instability is a radius of the particles $`a`$. The instability criterion is $`a>a_{\mathrm{cr}}a_{}`$ for which $`(\text{div}𝒗)^2=|\text{rot}𝒗|^2`$. For the droplets in the atmosphere $`a_{}30\mu `$m. The growth rate of the clustering instability $`\gamma _{\text{cl}}\tau _\eta ^1(a/a_{})^4`$, where $`\tau _\eta `$ is the turnover time in the viscous scales of turbulence.

$``$ We introduced a new concept of compressibility of the turbulent diffusion tensor caused by a finite correlation time of an incompressible velocity field. For this model of the velocity field, the field of Lagrangian trajectories is not divergence-free.

$``$ We suggested a mechanism of saturation of the clustering instability - particle collisions in the clusters. An evaluated nonlinear level of the saturation of the droplets number density in clouds exceeds by the orders of magnitude their mean number density.

###### Acknowledgements.

We have benefited from discussions with I. Procaccia. This work was partially supported by the German-Israeli Project Cooperation (DIP) administrated by the Federal Ministry of Education and Research (BMBF), by the Israel Science Foundation and by INTAS (Grant 00-0309). DS is grateful to a special fund for visiting senior scientists of the Faculty of Engineering of the Ben-Gurion University of the Negev and to the Russian Foundation for Basic Research (RFBR) for financial support under grant 01-02-16158.

## Appendix A Basic equations in the model with a random renewal time

In this Appendix we derive Eq. (63) for the simultaneous second-order correlation function $`\mathrm{\Phi }(t,𝒓)`$ which serves as a basis for further analysis in Appendixes B and C under some simplifying model assumptions about the statistics of the velocity field.

SSmethodExact solution of dynamical equations for a given velocity field

#### A.0.1 Simple case: no molecular diffusion

Consider first Eq. (1) for the number density of particles $`n(t,𝒓)`$ in the case $`D=0`$:

$`{\displaystyle \frac{n(t,𝒓)}{t}}+\mathbf{}[n(t,𝒓)𝒗(t,𝒓)]=0,`$ (44)
when all particles are transported only by advection. Solution of Eq. (44) with the initial condition $`n(s,𝒓)`$ is given by

$`n(t,𝒓)`$ $`=`$ $`G(t,𝒓)n[s,𝝃_\text{L}(t,𝒓|s)],`$ (45)
where $`𝝃_\text{L}(t,𝒓|s)`$ is the Lagrangian trajectory of the particle which is located at coordinate $`𝒓`$ at time $`t`$. Here we label the particles at present moment of time $`t`$ and consider a current time $`s<t`$ as moments in the past. This differs from a usual approach, see Eqs. LABEL:rL, when particles are labelled at the initial time $`t_0`$, and a current time $`t>t_0`$. Therefore in the equations below it is more convenient to redefine Lagrangian displacement $`𝝆_\text{L}(t,r|s)\stackrel{~}{𝝆}_\text{L}(t,r|s)=𝝆_\text{L}(t,r|s)`$. Now Eqs. LABEL:rL can be written as

$`\stackrel{~}{𝝆}_\text{L}(t,𝒓|s)`$ $`=`$ $`{\displaystyle _s^t}𝒗[\tau ,𝝃_\text{L}(t,𝒓|\tau )]\text{d}\tau ,`$ (46)
$`𝝃_\text{L}(t,𝒓|s)`$ $``$ $`𝒓\stackrel{~}{𝝆}_\text{L}(t,𝒓|s).`$ (47)
The Green function is the functional of $`𝝃_\text{L}(t,𝒓|s)`$:

$`G(t,𝒓,s)`$ $`=`$ $`\mathrm{exp}\left\{{\displaystyle _s^t}b[\tau ,𝝃_\text{L}(t,𝒓|\tau )]𝑑\tau \right\},`$
$`b(t,𝒓)`$ $``$ $`\mathbf{}𝒗(t,𝒓),`$ (48)
Introduce the shift operator

$`\mathrm{exp}[\stackrel{~}{𝝆}_\text{L}\mathbf{}]=1\stackrel{~}{𝝆}_\text{L}\mathbf{}+{\displaystyle \frac{1}{2!}}[\stackrel{~}{𝝆}_\text{L}\mathbf{}]^2\mathrm{}`$ (49)
which acts as follows:

$`\mathrm{exp}[\stackrel{~}{𝝆}_\text{L}\mathbf{}]n(t,𝒓)=n(t,𝒓\stackrel{~}{𝝆}_\text{L}).`$ (50)
One can validate relation LABEL:shift1 by Taylor series expansion of the function $`n(t,𝒓\stackrel{~}{𝝆}_\text{L})`$. Now Eq. (45) can be rewritten as follows:

$`n(t,𝒓)=G(t,𝒓,s)\mathrm{exp}[\stackrel{~}{𝝆}_\text{L}(t,𝒓|s)\mathbf{}]n(s,𝒓).`$ (51)
#### A.0.2 Molecular diffusion as a Wiener process

Consider now the full Eq. (1) with $`D0`$ whereby particles are transported by both, fluid advection and molecular diffusion. It was found by Wiener (see, e.g., ZRS90 ) that Brownian motion (molecular diffusion) can be described by the  Wiener random process $`𝒘(t)`$ with the following properties:

$`𝒘(t)_𝒘=0,w_i(t+\tau )w_j(t)_𝒘=\tau \delta _{ij}.`$ (52)
Here $`\mathrm{}_𝒘`$ denotes the mathematical expectation over the statistics of the Wiener process. Introduce the Wiener trajectory $`𝝃_\text{W}(t,𝒓|s)`$ (which usually is called the Wiener path) and the Wiener displacement $`𝝆_\text{W}(t,𝒓|s)`$ as follows:

$`𝝃_𝒘(t,𝒓|s)`$ $``$ $`𝒓𝝆_\text{W}(t,𝒓|s),`$ (53)
$`𝝆_\text{W}(t,𝒓|s)`$ $`=`$ $`{\displaystyle _s^t}𝒗[\tau ,𝝃_𝒘(t,𝒓|\tau )]\text{d}\tau +\sqrt{2D}𝒘(ts).`$
Comparison of this formula with Eqs. LABEL:rL1 shows that in the limit $`D0`$, $`𝝃_\text{W}(t,𝒓|s)𝝃_\text{L}(t,𝒓|s)`$ and $`𝝆_\text{W}(t,𝒓|s)𝝆_\text{L}(t,𝒓|s)`$.

In Refs. \[EKR00, \] it was shown that solution of Eq. (1) (with $`D0`$) can be written as solution LABEL:A4 of Eq. (44) (with $`D=0`$) by replacement $`\stackrel{~}{𝝆}_\text{L}(t,𝒓|s)\stackrel{~}{𝝆}_\text{W}(t,𝒓|s)`$ and then averaging over the statistics of the Wiener processes LABEL:W:

$`n(t,𝒓)=G(t,𝒓,s)\mathrm{exp}[𝝆_\text{W}(t,𝒓|s)\mathbf{}]n(s,𝒓)_𝒘.`$
(54)
SSvelocityTwo-step averaging over velocity statistics

#### A.0.3 Model of a random velocity field

Note that Eq. (54) is a solution of Eq. (1) at a given realization of the random velocity field. Our next goal is to determine the simultaneous correlation functions

$`\overline{n}(t)`$ $`=`$ $`n(t,𝒓)_𝒗,`$ (55)
$`\mathrm{\Phi }(t,𝒓_2𝒓_1)`$ $`=`$ $`n(t,𝒓_1)n(t,𝒓_2)_𝒗\overline{n}^2(t),`$
averaged over the stationary, space homogeneous statistics of turbulent velocity field, where $`\mathrm{}_𝒗`$ denotes this averaging. Since the initial distribution $`n(t_0,𝒓)`$ is assumed to be homogeneous in space, $`\overline{n}(t)`$ is independent of spatial coordinate, and $`\mathrm{\Phi }(t,𝒓_2𝒓_1)`$ depends only on the difference $`𝒓_2𝒓_1`$.

In order to simplify the averaging procedure LABEL:def-cor we consider a model of random velocity field which fully looses memory at some instants of renewal $`\tau _j`$. For $`t_1`$ and $`t_2`$ inside a renewal interval \[$`\tau _j<t_1,t_2<\tau _{j+1}`$\] the velocity pair correlation function is defined as

$`^{\alpha \beta }(t_2t_1,𝒓_2𝒓_1)v_\alpha (t_1,𝒓_1)v_\beta (t_2,𝒓_2)_𝒗,`$ (56)
where $`\mathrm{}_𝒗`$ denotes averaging over “intrinsic statistics” of the velocity field. In our model the velocity fields before and after renewals are statistically independent. The interval between the renewal instants $`\tau _j`$ may be the same or randomly distributed, say with the Poisson statistics. In the latter case the full averaging $`\mathrm{}_𝒗`$ may be considered as a two-stage process. First one calculates $`\mathrm{}_𝒗`$ and then averages over the statistics of the renewal time $`\tau _{\text{ren}}`$, which is denoted as $`\mathrm{}_{\text{ren}}`$:

$`\mathrm{}_𝒗\mathrm{}_𝒗_{\text{ren}}.`$ (57)
For the Poisson statistics of $`\tau _j`$

$`F^{\alpha \beta }(t_2t_1,𝒓_2𝒓_1)v_\alpha (t_1,𝒓_1)v_\beta (t_2,𝒓_2)_𝒗`$ (58)
$`=`$ $`^{\alpha \beta }(t_2t_1,𝒓_2𝒓_1)\mathrm{exp}(|t_2t_1|/\overline{\tau }_{\text{ren}}),`$
where $`\overline{\tau }_{\text{ren}}`$ is a mean renewal time. It would be useful to define the correlation time of the function $`^{\alpha \beta }`$ as follows

$`\tau _v(R)={\displaystyle ^{\alpha \beta }(\tau ,𝑹)𝑑\tau /^{\alpha \beta }(0,𝑹)}.`$ (59)
Certainly this model of the random velocity field cannot be considered as universal. However, it reproduces important features of some flows (see, e.g., Ref. \[LST00, \]).

#### A.0.4 Averaging procedure

Our model involves three random processes:

1. The Wiener random process which describes Brownian (molecular) diffusion.
2. Poisson process for a random renewal time.
3. The random velocity field between the renewals.

Eq. (54) presents $`n(t,𝒓)`$ after the first step, i.e., it describes the number density at a given realization of a velocity field. Using Eq. (54) we obtain

$`n(t,𝒓_1)n(t,𝒓_2)=G(𝒓_1)G(𝒓_2)\mathrm{exp}[𝝃^{}(𝒓_1)\mathbf{}_1`$
$`+𝝃^{}(𝒓_2)\mathbf{}_2]n(s,𝒓_1)n(s,𝒓_2)_{𝒘𝒘},`$ (60)
where $`\mathbf{}_1=/𝒓_1`$ and $`\mathbf{}_2=/𝒓_2`$ and $`_{𝒘𝒘}`$ denotes averaging over two independent Wiener processes determining two Wiener paths. Hereafter for simplicity we use the following notations: $`G(𝒓)G(t,𝒓,s)`$ and $`𝝃^{}(𝒓)𝝃^{}(t,𝒓|s).`$

Now we average Eq. (60) over a random velocity field for a given realization of a Poisson process:

$`\stackrel{~}{\mathrm{\Phi }}(t,𝒓_2𝒓_1)=n(t,𝒓_1)n(t,𝒓_2)_𝒗(\overline{n})^2`$
$`=`$ $`G(𝒓_1)G(𝒓_2)\mathrm{exp}[𝝃^{}(𝒓_1)\mathbf{}_1+𝝃^{}(𝒓_2)\mathbf{}_2]_{𝒘𝒘}_𝒗`$
$`\times \stackrel{~}{\mathrm{\Phi }}(t_0,𝒓_1𝒓_2)\}.`$
Here the time $`t_0`$ is the last renewal time before time $`t`$ and $`t^{}=tt_0`$ is a random variable. Thus, averaging of the functions

$$G(𝒓_1)G(𝒓_2)\mathrm{exp}[𝝃^{}(𝒓_1)\mathbf{}_1+𝝃^{}(𝒓_2)\mathbf{}_2],\mathrm{\Phi }(t_0,𝒓_1𝒓_2)$$

is decoupled into two time intervals because the first function is determined by the velocity field after the renewal while the second function $`\mathrm{\Phi }(t_0,𝒓_1𝒓_2)`$ is determined by the velocity field before renewal. Now we take into account that for the Poisson process any instant can be chosen as the initial instant. We average Eq. (A.0.4) over the random renewal time. The probability density $`p(t^{})`$ for a random renewal time is given by

$`p(t^{})=\overline{\tau }_{\text{ren}}^1\mathrm{exp}(t^{}/\overline{\tau }_{\text{ren}}).`$ (62)
Thus the resulting averaged equation for “fully” averaged correlation function $`\mathrm{\Phi }(t,𝑹)=\stackrel{~}{\mathrm{\Phi }}(t,𝑹)_{\text{ren}}`$, defined by Eq. (55), assumes the following form:

$`\mathrm{\Phi }(t,𝑹)=\overline{\tau }_{\text{ren}}^1{\displaystyle _0^t}\widehat{P}(\tau ,𝑹)\mathrm{\Phi }(t\tau ,𝑹)\mathrm{exp}(\tau /\overline{\tau }_{\text{ren}})𝑑\tau `$
$`+\mathrm{exp}(t/\overline{\tau }_{\text{ren}})\widehat{P}(t,𝑹)\mathrm{\Phi }_0(𝑹).`$ (63)
The first term in Eq. (63) describes the case when there is at least one renewal of the velocity field during the time $`t`$ (i.e., the Poisson event), whereas the second term describes the case when there is no renewal during the time $`t`$. Here $`\mathrm{\Phi }_0(𝑹)=\mathrm{\Phi }(t=0,𝑹)`$ and

$`\widehat{P}(t,𝑹)=G(𝒓_1)G(𝒓_2)`$ (64)
$`\times \mathrm{exp}[𝝃^{}(𝒓_1)\mathbf{}_1+𝝃^{}(𝒓_2)\mathbf{}_2]_{𝒘𝒘}_𝒗`$
$`=\mathrm{exp}g(𝒓_1)+g(𝒓_2)+𝝃^{}(𝒓_1)\mathbf{}_1+𝝃^{}(𝒓_2)\mathbf{}_2_{𝒘𝒘}_𝒗,`$
where $`G(𝒓)=\mathrm{exp}[g(𝒓)]`$. Equation LABEL:C14 is simplified in Appendixes B and C under the additional assumptions about the velocity field statistics.

## Appendix B velocity field with Gaussian Lagrangian trajectories

Consider a model of a random velocity field where Lagrangian trajectories, i.e., the integrals $`𝐯(\mu ,𝛏)𝑑\mu `$ and $`b(\mu ,𝛏)𝑑\mu `$ have Gaussian statistics. Using an identity $`\mathrm{exp}(a\eta )_\eta =\mathrm{exp}(\frac{1}{2}a^2)`$ in Eq. (64) we obtain

$`\widehat{P}(\mu ,𝑹)`$ $`=`$ $`\mathrm{exp}[\mu \widehat{}],`$ (65)
where

$`\widehat{}`$ $`=`$ $`B(𝑹)+2U_\alpha (𝑹)_\alpha +\widehat{D}_{\alpha \beta }(𝑹)_\alpha _\beta ,`$
$`\mu B(𝑹)`$ $`=`$ $`g(𝒓_1)g(𝒓_2)_{𝒘𝒘}_𝒗,`$ (66)
$`\mu U_\alpha (𝑹)`$ $`=`$ $`\xi _\alpha ^{}(𝒓_1)g(𝒓_2)_{𝒘𝒘}_𝒗,`$
$`\widehat{D}_{\alpha \beta }(𝑹)`$ $`=`$ $`D_{\alpha \beta }(0)D_{\alpha \beta }(𝑹),`$
$`\mu D_{\alpha \beta }(𝑹)`$ $`=`$ $`\xi _\alpha ^{}(𝒓_1)\xi _\beta ^{}(𝒓_2)_{𝒘𝒘}_𝒗.`$
Here $`\eta `$ is a Gaussian random variable with zero mean value and unit variance and $`G(𝒓)_𝒘_𝒗=1.`$ The latter yields

$$g_𝒘_𝒗=\frac{1}{2}\stackrel{~}{g}^2_𝒘_𝒗,g=g_𝒘_𝒗+\stackrel{~}{g},$$

where $`\stackrel{~}{g}_𝒘_𝒗=0.`$ When correlation time $`\tau _v(R)`$, Eq. (59), is much less then the current time $`t`$ and $`\overline{\tau }_{\text{ren}}`$, these correlation functions are given by

$`B(𝑹)`$ $`=`$ $`2{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}b[0,𝝃(𝒓_1)]b[\mu ^{},𝝃(𝒓_2)]_{𝒘𝒘}_𝒗𝑑\mu ^{},`$ (67)
$`U_\alpha (𝑹)`$ $`=`$ $`2{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}v_\alpha [0,𝝃(𝒓_1)]b[\mu ^{},𝝃(𝒓_2)]_{𝒘𝒘}_𝒗𝑑\mu ^{},`$
$`D_{\alpha \beta }(𝑹)`$ $`=`$ $`2{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}v_\alpha [0,𝝃(𝒓_1)]v_\beta [\mu ^{},𝝃(𝒓_2)]_{𝒘𝒘}_𝒗𝑑\mu ^{},`$
where we used an identity

$`{\displaystyle \underset{0}{\overset{\mu }{}}}a_\alpha (\mu ^{},𝒓_1)𝑑\mu ^{}{\displaystyle \underset{0}{\overset{\mu }{}}}c_\beta (\mu ^{\prime \prime },𝒓_2)𝑑\mu ^{\prime \prime }_𝒘_𝒗`$
$`2\mu {\displaystyle \underset{0}{\overset{\mathrm{}}{}}}a_\alpha (0,𝒓_1)c_\beta (\mu ^{},𝒓_2)_𝒘_𝒗𝑑\mu ^{}.`$
Eq. (65) allows to rewrite Eq. (63) as

$`\mathrm{\Phi }(t,𝑹)`$ $`=`$ $`{\displaystyle \frac{1}{\overline{\tau }_{\text{ren}}}}\left[{\displaystyle _0^t}\mathrm{exp}(\mu \widehat{}_1)𝑑\mu \right]\mathrm{\Phi }(t,𝑹)`$ (68)
$`+\mathrm{exp}(t\widehat{}_1)\mathrm{\Phi }(t,𝑹),`$
where

$`\widehat{}_1=\widehat{}{\displaystyle \frac{}{t}}{\displaystyle \frac{1}{\overline{\tau }_{\text{ren}}}}.`$ (69)
To derive Eq. (68) we used the following identity

$`\mathrm{\Phi }(t\mu ,𝑹)=\mathrm{exp}\left(\mu {\displaystyle \frac{}{t}}\right)\mathrm{\Phi }(t,𝑹),`$ (70)
which follows from the Taylor expansion

$`f(t+\tau )={\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}\left(\tau {\displaystyle \frac{}{t}}\right)^m{\displaystyle \frac{f(t)}{m!}}=\mathrm{exp}\left(\tau {\displaystyle \frac{}{t}}\right)f(t).`$ (71)
In particular,

$$\mathrm{\Phi }_0(𝑹)=\mathrm{\Phi }(tt,𝑹)=\mathrm{exp}\left(t\frac{}{t}\right)\mathrm{\Phi }(t,𝑹).$$

Evaluating the integral in Eq. (68) we obtain

$`\left[\mathrm{exp}(t\widehat{}_1)1\right]\left(\widehat{}_1+\overline{\tau }_{\text{ren}}^1\right)\mathrm{\Phi }(t,𝑹)=0.`$ (72)
Here we used the commutativity relation

$$\widehat{}_1\mathrm{exp}(t\widehat{}_1)=\mathrm{exp}(t\widehat{}_1)\widehat{}_1.$$

Thus, finally

$`{\displaystyle \frac{\mathrm{\Phi }}{t}}=\left[B(𝑹)+2𝐔(𝑹)\mathbf{}+\widehat{D}_{\alpha \beta }(𝑹)_\alpha _\beta \right]\mathrm{\Phi }(t,𝑹),`$
(73)
Note that in the limit $`\overline{\tau }_{\text{ren}}\mathrm{}`$, Eq. (73) describes the evolution of $`\mathrm{\Phi }(t,𝑹)`$ in the model of the random velocity field without renewals.

## Appendix C Gaussian velocity field with a small yet finite correlation time

Here we consider a random Gaussian velocity field with a small $`\overline{\tau }_{\text{ren}}`$. Using Eq. (70) we rewrite Eq. (63) in the form

$$\left\{\frac{1}{\overline{\tau }_{\text{ren}}}\underset{0}{\overset{t}{}}\widehat{P}(\tau ,𝑹)\mathrm{exp}\left(\frac{\tau }{\overline{\tau }_{\text{ren}}}\widehat{M}\right)𝑑\tau 1\right\}\mathrm{\Phi }(t,𝑹)=0,$$
(74)
where $`\widehat{M}=1+\overline{\tau }_{\text{ren}}(/t)`$ and we neglected the last term in Eq. (63) for small $`\tau `$. Expanding the function $`\widehat{P}(\tau ,𝑹)`$ in Taylor series in the vicinity of $`\tau =0`$ we obtain

$`\left\{{\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}\overline{\tau }_{\text{ren}}^k\left[{\displaystyle \frac{^k\widehat{P}(\tau ,𝑹)}{\tau ^k}}\right]_{\tau =0}\widehat{M}^{(k+1)}1\right\}\mathrm{\Phi }(t,𝑹)=0,`$
(75)
where we used that

$$\underset{0}{\overset{t}{}}\tau ^k\mathrm{exp}\left(\frac{\tau }{\overline{\tau }_{\text{ren}}}\widehat{M}\right)𝑑\tau =k!\overline{\tau }_{\text{ren}}^{k+1}\widehat{M}^{(k+1)}.$$

Neglecting the terms $`O(\overline{\tau }_{\text{ren}}^5)`$ in Eq. (75) we obtain

$`\widehat{M}^2{\displaystyle \frac{\mathrm{\Phi }(t,𝑹)}{t}}`$ $`=\overline{\tau }_{\text{ren}}[\left({\displaystyle \frac{^2\widehat{P}(\tau ,𝑹)}{\tau ^2}}\right)_{\tau =0}`$
$`+\overline{\tau }_{\text{ren}}^2\left({\displaystyle \frac{^4\widehat{P}(\tau ,𝑹)}{\tau ^4}}\right)_{\tau =0}]\mathrm{\Phi }(t,𝑹),`$
where we used that the expansion of the operator $`\widehat{P}(\tau ,𝑹)`$ into Taylor series (for small $`\tau )`$ for a random Gaussian velocity field has only even powers of $`\tau .`$ Thus, the equation for the correlation function $`\mathrm{\Phi }(t,𝑹)`$ is given by

$`\widehat{M}^2{\displaystyle \frac{\mathrm{\Phi }(t,𝑹)}{t}}`$ $`=`$ $`[B(𝑹)+2𝑼(𝑹)\mathbf{}+\widehat{D}_{\alpha \beta }(𝑹)_\alpha _\beta ]\mathrm{\Phi },`$ (77)
where

$`\widehat{D}_{\alpha \beta }(𝑹)`$ $`=`$ $`{\displaystyle \frac{\overline{\tau }_{\text{ren}}}{2}}\stackrel{~}{\xi }_\alpha \stackrel{~}{\xi }_\beta G(𝒓_1)G(𝒓_2)_{𝒘𝒘}_𝒗,`$ (78)
$`U_\alpha (𝑹)`$ $`=`$ $`{\displaystyle \frac{1}{\overline{\tau }_{\text{ren}}}}g(𝒓_2)\xi _\alpha ^{}(𝒓_1)_{𝒘𝒘}_𝒗+{\displaystyle \frac{1}{2\overline{\tau }_{\text{ren}}}}g(𝒓_1)g(𝒓_2)\stackrel{~}{\xi }_\alpha _{𝒘𝒘}_𝒗,`$ (79)
$`B(𝑹)`$ $`=`$ $`{\displaystyle \frac{1}{\overline{\tau }_{\text{ren}}}}g(𝒓_1)g(𝒓_1)_{𝒘𝒘}_𝒗,`$ (80)
Here for the homogeneous turbulent velocity field EKRS99 :

$$\stackrel{~}{𝝃}=𝝃^{}(𝒓_2)𝝃^{}(𝒓_1),\mathbf{}=/𝑹,G=\overline{G}+g,g_{𝒘𝒘}_𝒗=0,\overline{G}=G_{𝒘𝒘}_𝒗=1.$$

Using the expansion of $`𝝃(\overline{\tau }_{\text{ren}},𝒓)`$ and $`g[\overline{\tau }_{\text{ren}},𝝃(𝒓)]`$ into Taylor series of a small time $`\overline{\tau }_{\text{ren}}`$ after the lengthly algebra we obtain

$`\widehat{D}_{\alpha \beta }(𝑹)`$ $`=`$ $`2D\delta _{\alpha \beta }+2\overline{\tau }_{\text{ren}}[\stackrel{~}{f}_{\alpha \beta }(𝑹)+\mathrm{St}^2Q_{\alpha \beta }(𝑹)],`$ (81)
$`Q_{\alpha \beta }(𝑹)`$ $`=`$ $`3[(_\nu f_{\mu \beta })(_\mu f_{\alpha \nu })\stackrel{~}{f}_{\mu \nu }_\nu _\mu f_{\alpha \beta }]+24A_\alpha A_\beta +12(A_\mu _\mu f_{\alpha \beta }\stackrel{~}{f}_{\alpha \beta }_\mu A_\mu )20\stackrel{~}{f}_{\alpha \mu }_\mu A_\beta ,`$ (82)
$`U_\alpha (𝑹)`$ $`=`$ $`2\overline{\tau }_{\text{ren}}\{A_\alpha \mathrm{St}^2[(_\nu A_\mu )(_\mu f_{\alpha \nu })+10A_\mu _\mu A_\alpha +12A_\alpha _\mu A_\mu ]\},`$ (83)
$`B(𝑹)`$ $`=`$ $`2\overline{\tau }_{\text{ren}}\{_\mu A_\mu +\mathrm{St}^2[(_\nu A_\mu )(_\mu A_\nu )6(_\mu A_\mu )^2]\},`$ (84)
$`A_\alpha `$ $`=`$ $`_\beta f_{\alpha \beta },\stackrel{~}{f}_{\alpha \beta }=f_{\alpha \beta }(0)f_{\alpha \beta }(𝑹),f_{\alpha \beta }(𝑹)=v_\alpha (𝒓_1)v_\beta (𝒓_2)_𝒗,`$
and St$`=\overline{\tau }_{\text{ren}}/\tau _\eta `$ is the Strouhal number. In these calculations we neglected small terms $`O(\mathrm{St}^2R^3^3).`$ Our analysis showed that the neglected small terms do not affect the growth rate of the clustering instability. In Eqs. (82)-(84) we assumed that the correlation function $`f_{\alpha \beta }`$ for homogeneous, isotropic and compressible velocity field is given by

$$f_{\alpha \beta }(𝑹)=\frac{u_\eta ^2}{3}\left[(F+F_c)\delta _{\alpha \beta }+\frac{RF^{}}{2}P_{\alpha \beta }+RF_c^{}R_{\alpha \beta }\right],$$
(85)
(see EKR95 ), and in scales $`0<R1`$ incompressible $`F(R)`$ and compressible $`F_c(R)`$ components of the random velocity field are given by

$$F(R)=(1R^2)/(1+\sigma _v),F_c(R)=\sigma _vF(R),$$

in scales $`R1`$ the functions $`F=F_c=0.`$ Here $`R`$ is measured in the units of $`\eta ,`$ $`P_{\alpha \beta }(R)=\delta _{\alpha \beta }R_{\alpha \beta },`$ $`R_{\alpha \beta }=R_\alpha R_\beta /R^2`$ and $`F^{}=dF/dR.`$ Turbulent diffusion tensor $`D_{\alpha \beta }(R)`$ is determined by the field of Lagrangian trajectories $`𝝃`$ \[see Eq. (78)\]. Due to a finite correlation time of a random velocity the field of Lagrangian trajectories $`𝝃`$ is compressible even if the velocity field is incompressible $`(\sigma _v=0).`$ Indeed, for $`\sigma _v=0`$ we obtain

$$(\mathbf{}𝝃)^2_𝒘_𝒗=\frac{20}{3}\mathrm{St}^4.$$

Using Eqs. (81)-(85) we calculate the functions $`\widehat{D}_{\alpha \beta }(𝑹),`$ $`U_\alpha (𝑹)`$ and $`B(𝑹):`$

$`\widehat{D}_{\alpha \beta }(𝑹)`$ $`=`$ $`[2D+R^2(a_3+\mathrm{St}^2b_6)]\delta _{\alpha \beta }`$ (86)
$`+R^2(a_4+\mathrm{St}^2b_4)R_{\alpha \beta },`$
$`U_\alpha (𝑹)`$ $`=`$ $`R_\alpha (a_5+\mathrm{St}^2b_5),`$ (87)
$`B`$ $`=`$ $`a_6+\mathrm{St}^2b_2,`$ (88)
where $`b_2=\frac{51}{47}b_5`$ and

$`a_5`$ $`=`$ $`{\displaystyle \frac{20\sigma _v}{3(1+\sigma _v)}}={\displaystyle \frac{a_2}{3}},a_3={\displaystyle \frac{2\sigma _v+4}{3(1+\sigma _v)}},`$
$`a_4`$ $`=`$ $`{\displaystyle \frac{4\sigma _v2}{3(1+\sigma _v)}},b_5={\displaystyle \frac{2350}{27}}\left({\displaystyle \frac{\sigma _v}{1+\sigma _v}}\right)^2,`$
$`b_6`$ $`=`$ $`{\displaystyle \frac{12+872\sigma _v+433\sigma _v^2}{27(1+\sigma _v)^2}},`$
$`b_4`$ $`=`$ $`{\displaystyle \frac{2(12203\sigma _v+1033\sigma _v^2)}{27(1+\sigma _v)^2}}.`$
We will show here that the combined effect of particles inertia $`(\sigma _v0)`$ and finite correlation time $`(\mathrm{St}0)`$ results in the excitation of the clustering instability whereby under certain conditions there is a self-excitation of the second moment of a number density of inertial particles. This instability causes formation of small-scale inhomogeneities in spatial distribution of inertial particles.

The equation for the second-order correlation function for the number density of inertial particles reads

$$\widehat{M}^2\frac{\mathrm{\Phi }(t,R)}{t}=\frac{\mathrm{\Phi }^{\prime \prime }}{m(R)}+\stackrel{~}{\lambda }(R)\mathrm{\Phi }^{}+B\mathrm{\Phi }$$
(89)
\[see Eqs. (88)\], where the time $`t`$ is measured in units of $`t_\eta ,`$ and

$`\mathrm{\Phi }^{}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Phi }}{R}},\mathrm{\Phi }^{\prime \prime }={\displaystyle \frac{^2\mathrm{\Phi }}{R^2}},{\displaystyle \frac{1}{m}}={\displaystyle \frac{2(1+X^2)}{\text{Sc}}},`$
$`\stackrel{~}{\lambda }`$ $`=`$ $`{\displaystyle \frac{2[2+X^2(1+2C)]}{R\mathrm{Sc}}},C={\displaystyle \frac{a_1+\mathrm{St}^2b_1}{4\beta }},`$
$`\beta `$ $`=`$ $`{\displaystyle \frac{a_2+\mathrm{St}^2b_3}{2}},X(R)=\sqrt{\mathrm{Sc}\beta }R,R=|𝒓_2𝒓_1|,`$
$`a_1`$ $`=`$ $`{\displaystyle \frac{2(19\sigma _v+3)}{3(1+\sigma _v)}},a_2={\displaystyle \frac{2(3\sigma _v+1)}{3(1+\sigma _v)}},`$
$`b_1`$ $`=`$ $`{\displaystyle \frac{1}{27(1+\sigma _v)^2}}(121278\sigma _v3067\sigma _v^2),`$
$`b_3`$ $`=`$ $`{\displaystyle \frac{1}{27(1+\sigma _v)^2}}(36+466\sigma _v+2499\sigma _v^2).`$
In order to obtain a solution of Eq. (89) we use a separation of variables, i.e., we seek for a solution in the following form:

$$\mathrm{\Phi }(t,R)=\widehat{\mathrm{\Phi }}(R)\mathrm{exp}(\gamma _2t),$$

whereby $`\gamma _2`$ is a free parameter which is determined using the boundary conditions

$$\widehat{\mathrm{\Phi }}(R=0)=1,\widehat{\mathrm{\Phi }}(R\mathrm{})=0.$$

Here $`\gamma _2`$ is measured in units of $`1/t_\eta `$. Since the function $`\mathrm{\Phi }(t,R)`$ is the two-point correlation function, it has a global maximum at $`R=0`$ and therefore it satisfies the conditions:

$`\widehat{\mathrm{\Phi }}^{}(R=0)`$ $`=`$ $`0,\widehat{\mathrm{\Phi }}^{\prime \prime }(R=0)<0,`$
$`\widehat{\mathrm{\Phi }}(R=0)`$ $`>`$ $`|\widehat{\mathrm{\Phi }}(R>0)|.`$
Then Eq. (89) yields

$`\mathrm{\Gamma }\widehat{\mathrm{\Phi }}(R)={\displaystyle \frac{1}{m(R)}}\widehat{\mathrm{\Phi }}^{\prime \prime }+\stackrel{~}{\lambda }(R)\widehat{\mathrm{\Phi }}^{}+B\widehat{\mathrm{\Phi }},`$ (90)
where $`\mathrm{\Gamma }=\gamma _2(1+\overline{\tau }_{\text{ren}}\gamma _2)^2`$ . Equation (90) has an exact solution for $`0R<1:`$

$`\widehat{\mathrm{\Phi }}(X)`$ $`=`$ $`S(X)X(1+X^2)^{\mu /2},`$ (91)
$`S(X)`$ $`=`$ $`\text{Re}\{A_1P_\zeta ^\mu (iX)+A_2Q_\zeta ^\mu (iX)\},`$
and $`P_\zeta ^\mu (Z)`$ and $`Q_\zeta ^\mu (Z)`$ are the Legendre functions with imaginary argument

$$Z=iX,\mu =C\frac{3}{2},\zeta =\frac{1}{2}\pm \sqrt{C^2\kappa },\kappa =\frac{B\mathrm{\Gamma }}{2\beta }.$$

Solution of Eq. (89) can be analyzed using asymptotics of the exact solution (91). This asymptotic analysis is based on the separation of scales (see, e.g., ZMR88 ; EKR95 ). In particular, the solution of Eq. (89) has different regions where the form of the functions $`m(R)`$ and $`\stackrel{~}{\lambda }(R)`$ are different. The functions $`\widehat{\mathrm{\Phi }}(R)`$ and $`\widehat{\mathrm{\Phi }}^{}(R)`$ in these different regions are matched at their boundaries in order to obtain continuous solution for the correlation function. Note that the most important part of the solution is localized in small scales (i.e., $`R1).`$ Using the asymptotic analysis of the exact solution for $`X1`$ allowed us to obtain the necessary conditions of a small-scale instability of the second moment of a number density of inertial particles. The results obtained by this asymptotic analysis are presented below.

The solution (91) has the following asymptotics: for $`X1`$ (i.e., in the scales $`0R1/\sqrt{\mathrm{Sc}})`$ the solution for the second moment $`\widehat{\mathrm{\Phi }}`$ is given by

$`\widehat{\mathrm{\Phi }}(X)=\{1(\kappa /6)[X^2+O(X^4)]\}.`$ (92)
For $`X1`$ (i.e., in the scales $`1/\sqrt{\mathrm{Sc}}R<1)`$ the function $`\widehat{\mathrm{\Phi }}`$ is given by

$$\widehat{\mathrm{\Phi }}(X)=\text{Re}\{AX^{C\pm \sqrt{C^2\kappa }}\}.$$
(93)
When $`C^2\kappa <0`$ the second-order correlation function for a number density of inertial particles $`\widehat{\mathrm{\Phi }}`$ is given by

$$\widehat{\mathrm{\Phi }}(R)=A_3R^C\mathrm{cos}(\nu _I\mathrm{ln}R+\phi ),\nu _I=\sqrt{\kappa C^2},$$

where $`C>0`$ and $`\phi `$ is the argument of the complex constant $`A.`$ For $`R1`$ the second-order correlation function for the number density of inertial particles is given by

$`\widehat{\mathrm{\Phi }}(R)=(A_4/R)\mathrm{exp}(R\sqrt{3\mathrm{\Gamma }/2}),`$ (94)
where $`\mathrm{\Gamma }>0`$. Since the total number of particles in a closed volume is conserved, i.e., particles can only be redistributed in the volume,

$$_0^{\mathrm{}}R^2\widehat{\mathrm{\Phi }}(R)𝑑R=\widehat{\mathrm{\Phi }}(k=0)=0.$$

The latter yields $`\phi =\pi /2`$ for $`\mathrm{ln}\mathrm{Sc}1`$ and $`\mathrm{\Gamma }1.`$ When $`C^2\kappa >0`$, the solution (93) cannot be matched with solutions (92) and (94). Thus, the condition $`C^2\kappa <0`$ is the necessary condition for the existence of the solution for the correlation function. The condition $`C>0`$ provides the existence of the global maximum of the correlation function at $`R=0`$.

Matching functions $`\widehat{\mathrm{\Phi }}`$ and $`\widehat{\mathrm{\Phi }}^{}`$ at the boundaries of the above-mentioned regions yields coefficients $`A_k`$ and $`\mathrm{\Gamma }.`$ In particular, the eigenvalue $`\mathrm{\Gamma }`$ is given by Eq. (III.2).


text:
# Evidence of fast rotation in dwarf elliptical galaxies

## 1 Introduction

Although dwarf elliptical galaxies (dEs) constitute the dominant galaxy population in nearby clusters, their origin and nature are still unclear. One of the main questions under debate is whether dEs represent the extension of the “classical” ellipticals (Es) to lower luminosities, or, on the contrary, whether they are the result of distinct formation and evolution processes (for a review see, e.g., Ferguson & Binggeli 1994). While early observational works suggested a dichotomy in the structural properties of both galaxy families (e.g. Kormendy 1985), Binggeli & Jerjen (1998) (see also Jerjen & Binggeli 1997) have shown that, if a Sérsic law is used to fit the observed luminosity profiles, dEs appear as the true low-luminosity extension of the classical Es. Previous workers have also studied the possible dichotomy between dEs and Es by comparing the flattening distributions (e.g. Ryden & Terndrup 1994), clustering properties (Conselice, Gallagher & Wyse 2001) and stellar populations (Gorgas et al. 1997, hereafter G97; Drinkwater et al. 2001) of dEs in nearby clusters. Another key issue that should provide important clues to the open question of the possible dichotomy is the analysis of rotational properties.

From the work of Davies et al. (1983) (see also Halliday et al. 2001), it is well known that low-luminosity E galaxies, in contrast with giant Es, are supported by rotation. If dEs were the extension to even lower luminosities of the classical Es, one should expect to find relatively high rotational velocities along the major axis of these objects. The observational data in this subject require high signal-to-noise (S/N) spectra at faint levels of surface brightnesses and it is, therefore, very scarce. Until recently, rotation curves had only been measured for 6 dEs: two dwarfs in Virgo (VCC 351 and IC 794; Bender & Nieto 1990, hereafter BN90), the 3 dwarfs companions of M31 (e.g. Bender, Paquet & Nieto 1991) and the Fornax dE (Mateo et al. 1991). In all these cases, dEs were found to rotate too slowly to be consistent with an oblate, isotropic body flattened by rotation. Very recently, new data have been added to this sparse statistics. On the one hand, Geha, Guhathakurta & van der Marel (2001, hereafter GGV01) have presented kinematical data for a sample of 4 Virgo dEs, finding that, in agreement with the previous results, they are slow rotators. On the other hand, De Ricke et al. (2001, hereafter DR01) have obtained deep spectroscopic data of the dE FS76 (in the NGC 5044 group), showing that the galaxy is rotating as fast as predicted by the isotropic models. In this Paper we revisit the question of rotation in dEs by practically doubling the observational sample so far available. We will show that most of the galaxies in the new sample tend to be rotationally supported and, therefore, fast rotation in dwarf elliptical galaxies is not a rare event, as previously thought.

## 2 Galaxy sample, observations and data reduction

The galaxies analysed in this work correspond to those already presented in G97. This sample of Virgo galaxies, which is listed in Table 1, comprises 4 dwarf ellipticals and two dwarf lenticular galaxies. Although NGC 4489 is classified as S0 (non dwarf) in Binggeli, Sandage & Tammann (1985), it exhibits a mean surface brightness within the effective radius of $`SB_e=22.24`$ mag arcsec-2 (Faber et al. 1989), which, in the $`SB_e`$–absolute magnitude plane, places it closer to the locus of dwarf ellipticals than to the low-luminosity early-type galaxies111At the absolute magnitude of NGC 4489, classical ellipticals and S0’s attain typical values of $`SB_e`$ between 18 and 21 mag arcsec-2, whereas for the dwarf galaxies of this sample $`22.5<SB_e<23`$ mag arcsec-2.. It must be noted also that the present sample is biased toward the brightest dwarf early-type galaxies of the Virgo cluster.

Long-slit spectroscopic observations along the major axes of the galaxies were carried out with the INT at La Palma. Although the relevant observational parameters were already described in G97, we want to highlight that high-quality spectra, covering the 4700–6100 Å range with a spectral resolution of $`\sigma =60`$ km s-1 ($`\mathrm{FWHM}2.5`$ Å) were obtained out to, typically, the effective radii of the galaxies.

The data reduction and analysis were performed using the REDucm E package (Cardiel et al. 1998; Cardiel 1999), specially developed for the parallel treatment of data and errors in the reduction of spectroscopic observations. In this sense, each fully processed spectrum is accompanied by an associated error spectrum, which contains the propagation of the initial random errors (due to photon statistics and read-out noise) throughout the cascade of arithmetic manipulations in the reduction procedure. Full details concerning the data reduction will be presented in Pedraz (2002).

We complemented the observations of the targets with a large set (39) of template stars from the Lick/IDS stellar library (Gorgas et al. 1993; Worthey et al. 1994). The template list was selected to provide an appropriate coverage in stellar atmospheric parameters ($`T_{\mathrm{eff}}`$ from 3247 to 6575 K, $`\mathrm{log}g`$ from 1.11 to 4.84, and \[Fe/H\] from $`2.15`$ to 0.35), suitable for the modelling of old stellar populations.

## 3 Measurement of kinematic data

In order to estimate radial velocities and velocity dispersions for each galaxy spectrum, we wrote a dedicated program within REDucm E which follows the MOVEL and OPTEMA algorithms described by González (1993). The MOVEL algorithm, an improvement of the classic Fourier quotient method by Sargent et al. (1977), is an iterative procedure in which a galaxy model is processed in parallel to the galaxy spectrum. In this way, a comparison between the input and recovered broadening functions for the model allows to correct the galaxy power spectrum from any imperfections of the data handling in Fourier space. The main improvement of the procedure is introduced through the OPTEMA algorithm, which is able to overcome the typical template mismatch problem. First, the 39 individual stellar spectra of the observing run were binned to create a set of 8 template spectra corresponding to different spectral types and luminosity classes (namely F0V, G0V, K0V, M0V, G5III, G8III, K0III, K5III). For each galaxy central spectrum, these 8 template spectra are scaled, shifted and broadened according to a first guess of the $`\gamma `$ (mean line strength), $`V`$ (radial velocity)and $`\sigma `$ (velocity dispersion) parameters. The next step is to find the linear combination of the template spectra that best matches the observed galaxy spectrum. This provides a first composite template which is fed to the MOVEL algorithm. The output kinematic parameters are then used to create an improved composite template and the process is iterated until it converges. This iterative approach then provides an optimal template while simultaneously computing the radial velocity and velocity dispersion of the galaxy spectrum. In Figure 1 we show a typical fit between the central spectrum of a galaxy of the sample and its corresponding optimal template corrected with the derived kinematic parameters. It must be noted that, once an optimal template is derived for the central spectrum of a given galaxy, this same template is used to derive the kinematics of the different spectra along the galaxy radius. This is a fair approximation since these galaxies exhibit relatively flat stellar population gradients along the radius (G97).

For each galaxy spectrum, random errors in the derived kinematic parameters were computed by numerical simulations. In each simulation, a bootstrapped galaxy spectrum, computed from the error spectra provided by the reduction with REDucm E by assuming Gaussian errors, is fed to the algorithms described above (note that a different optimal template is computed in each simulation). Errors in the parameters are then calculated as the unbiased standard deviations of the different solutions. These final errors are expected to be quite realistic, since they incorporate all the uncertainties of the whole reduction process, from the first reduction steps (e.g. flatfielding) to the final measurements of the kinematics.

In order to assess the sensitivity of the above fitting procedure to measure velocity dispersions of only a fraction of the spectral resolution ($`\sigma _{\mathrm{instr}}`$), we carried out several Monte Carlo simulations in which we added Gaussian noise to artificially broadened stars. To summarize, we conclude that, in order to measure a velocity dispersion of $`\sigma _{\mathrm{instr}}/2`$ (30 km/s for our instrumental setup) with a relative error below 50 % ($`<15`$ km/s) we needed galaxy spectra with a minimum S/N (per Å) of 20–25.

Central velocity dispersions $`\sigma _0`$, together with their corresponding random errors, for the galaxy sample are listed in Table 1. We must note that S/N for the central spectra ranged from 38 to 98. The central velocity dispersion of NGC 4489 has been already measured by several authors. Our result ($`\sigma _0=51.1\pm 1.8`$) agrees with some of the previous measurements, like those of Davies et al. (1983; $`\sigma _0=53\pm 15`$), Faber et al. (1989; $`\sigma _0=49\pm 7`$) and González (1993; $`\sigma _0=48\pm 3`$), but it is somewhat smaller than the values given by Simien & Prugniel (1997; $`\sigma _0=63\pm 13`$) and, especially, Smith et al. (2000; $`\sigma _0=72\pm 7`$). We must note that the spectral resolution of these two last works ($`\sigma _{\mathrm{inst}}100`$ km s-1) is dangerously high to derive accurate velocity dispersions at these levels. IC 794 has also been included in previous kinematics works. Our central velocity dispersion ($`\sigma _0=41.1\pm 3.2`$) is marginally smaller than the value obtained by BN90 ($`\sigma _0=52\pm 6`$) but it perfectly agrees with the recent result of GGV01 ($`\sigma _0=41.6\pm 0.9)`$222This value has been computed, after the digitization of Fig. 1 in GGV01, as the error weighted mean of all the measurements with $`r<2^{\prime \prime }`$.. Central velocity dispersions for NGC 4431, IC 3393 and UGC 7436 have already been listed in the compilation by Bender et al. (1992). We must note that these values (68, 55 and 45 km/s respectively) are significantly larger than our measurements (42, 25 and 35 km/s). The original source of those values is unpublished data from R. Bender and, therefore, we do not have information enough to untangle the origin of this discrepancy. In any case, given the high S/N ratios of our spectra and the comparison with the previous published works, we are confident that our velocity dispersions are free of significant biases.

Folded rotation velocity curves along the major axes of the galaxies are presented in Figure 2. For the outer regions, we co-added a sufficient number of spectra in the spatial direction to guarantee a minimum signal-to-noise per Å of 20. In the case of NGC 4415, NGC 4431 and UGC 7436 the profiles extend out to the effective radius of the galaxies, whereas, for the other three galaxies, the data reach roughly to half the effective radius. The quality of the data and reliability of the error bars is apparent by the typical agreement between symmetrical points at both galaxy sides. An exception is the rotation curve for NGC 4489. In this case, a displacement of 1.3 arcsec in the kinematic centre can explain the asymmetry.

## 4 Discussion

In order to quantify the rotational support of the objects we follow the usual procedure of comparing the maximum rotational velocity ($`V_{\mathrm{max}}`$) to the mean velocity dispersion ($`\overline{\sigma }`$). A conservative estimate of the parameter $`V_{\mathrm{max}}`$ was computed as the error-weighted mean of the two data pairs with the highest rotational velocities. In some galaxies, like NGC 4431 and UGC 7436, the rotation curves increase to the outermost points observed, so the derived $`V_{\mathrm{max}}`$ should be considered as lower limits to the actual values. To compute the mean velocity dispersion we coadded all the individual spectra with radii between $`1^{\prime \prime }`$ and the effective radius. Prior to this, we shifted all the spectra to the same wavelength scale by using the rotation curves displayed in Fig 2. The resulting spectra exhibit S/N ratios, per Å, ranging from 62 to 119. Final values of $`V_{\mathrm{max}}`$ and $`\overline{\sigma }`$, together with their corresponding errors, are listed in Table 1.

In Figure 3 we have plotted the ratio $`V_{\mathrm{max}}/\overline{\sigma }`$ against mean ellipticity (compiled from the most recent available estimates in the literature – see caption to Table 1). The previous measurements from BN90 (for VCC 351), GGV01 (for VCC 917, VCC 1254 and VCC 1876) and DR01 (for FS 76) are also included. For IC 794, also in the samples of BN90 and GGV01, we only plot the ratio computed in this work. However, it is important to note that our estimate for this galaxy ($`V_{\mathrm{max}}/\overline{\sigma }=0.08\pm 0.04`$) is in perfect agreement with the upper limits given by the previous works ($`<0.20`$ in BN90, and $`<0.12`$ in GGV01). This gives further support to the data presented in this paper. We must also note that, to be fully consistent, we have recalculated the $`V_{\mathrm{max}}/\overline{\sigma }`$ ratio for the galaxy in DR01 using the same procedure than for our sample (DR01 give a ratio of peak velocity to central velocity dispersion of $`0.33\pm 0.15`$). After digitizing the measurements in their Fig. 2, we obtain a ratio of $`0.30\pm 0.04`$.

The curve in Figure 3 shows the expected relationship for an oblate spheroid with isotropic velocity distribution, thus rotationally flattened, from Binney (1978). As it is apparent from this figure, while the dwarf galaxies included in the previous works (with the exception of DR01) were found to rotate too slowly to be consistent with the isotropic model, 5 out of the 6 dwarf early-type galaxies studied in this paper (and 3 out of the 4 dEs) are compatible with being rotationally supported. We recall that the rotational velocities for the two galaxies at $`\epsilon 0.4`$ could be in fact lower limits of their actual values.

To investigate this issue in more detail, in Figure 4 we present a plot of the anisotropy parameter (ratio of the observed value of $`V_{\mathrm{max}}/\overline{\sigma }`$ to the one predicted by the isotropic model) against absolute magnitude. In order to compare the properties of dEs with those of giant and low-luminosity Es, we have included the galaxies from the compilation of Bender, Burstein & Faber (1992). Concerning the three galaxies in common with their sample (NGC 4431, IC 3393 and UGC 7436; see above), we must note that their $`V_{\mathrm{max}}/\overline{\sigma }`$ values (1.00, 0.50 and 0.44, respectively) are not incompatible with ours ($`0.71\pm 0.16`$, $`0.65\pm 0.14`$, $`0.47\pm 0.11`$).

In their study of the fundamental plane, Bender et al. (1992) (see also BN90) noted that galaxies with anisotropic velocity distributions were found at both high and low values of mass. Apart from confirming these findings, the new results presented in Figure 4 show that the brightest early-type dwarfs of the sample tend to be rotationally supported, extending the behaviour of the low-luminosity ellipticals towards the bright end of the dwarf family. Furthermore, these data indicate that, for fainter luminosities ($`M_B17`$), and as it was the case for the bright ellipticals, dwarf galaxies exhibit a range of rotational properties (i.e. not all dEs are supported by anisotropy).

In spite of these results, the present sample is still too small to derive firm conclusions about the possible relation between rotational support and other parameters. For instance, there are nucleated and non-nucleated dwarfs among both the rotating and non-rotating subsamples. Note that the search for correlations between rotational support and other properties is essential to test the possible existence of a dichotomy within the dE family (as suggested by Ryden et al. 1999) similar to that found for the Es (in isophote shapes, rotation, and core properties).

It is also specially interesting to compare the rotational properties with the stellar population of dEs. We must highlight that, with the exception of NGC 4489 (a possible transition object between dS0 and S0), the only galaxy of our sample flattened by anisotropic velocity dispersion (IC 794) is clearly the most metal-rich and youngest object in G97 (see their Figure 1.b). If confirmed by additional stellar populations measurements of other non-rotating dEs, this still speculative result would give support to a scenario in which recent star formation episodes are linked to a decrease in rotational support. In this sense, this would agree with the suggestion by BN90 that supernova-driven winds, triggered by star formation, could cause anisotropic supported objects. Also, both low rotation and young stellar populations are simultaneously predicted by the harassment models of Moore, Lake and Katz (1998), in which cluster dE are the result of morphological transformation of accreted spiral galaxies. This tentative correlation between rotation and age is in the opposite sense to that found for classical ellipticals, in which fast-rotating Es tend to have disky isophotes and younger ages (Bender 1988; de Jong & Davies 1997).

In any case, and independently of the possible link between rotation and star formation, the present results are compatible with the diversity of rotational properties predicted by the harassment scenario (compare our Fig. 4 with Fig. 7 of Moore et al. 1998). More data is clearly needed to understand the driving parameters governing the dynamical properties of dEs, as well as the possible correlation of rotation with luminosity profiles, isophotal shapes, core properties (nucleated versus non-nucleated) and stellar populations.

## Acknowledgments

The INT is operated on the island of La Palma by the Royal Greenwich Observatory at the Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias. This work was supported in part by the Spanish Programa Nacional de Astronomía y Astrofísica under grant AYA2000-977.