Inertial-particle dynamics in turbulent flows: caustics, concentration fluctuations, and random uncorrelated motion

We discuss the relation between three recent approaches of describing the dynamics and the spatial distribution of particles suspended in turbulent flows: phase-space singularities in the inertial particle dynamics (caustics), real-space singularities of the deformation tensor, and random uncorrelated motion. We discuss how the phase- and real-space singularities are related. Their formation is well understood in terms of a local theory. We discuss implications for random uncorrelated motion. Our results are supported by results of direct numerical simulations of inertial particles in model flows.


Abstract.
We have performed numerical simulations of inertial particles in random model flows in the white-noise limit (at zero Kubo number, Ku = 0) and at finite Kubo numbers. Our results for the moments of relative inertial particle velocities are in good agreement with recent theoretical results (Gustavsson & Mehlig 2011a) based on the formation of phase-space singularities in the inertial particle dynamics (caustics). We discuss the relation between three recent approaches describing the dynamics and spatial distribution of inertial particles suspended in turbulent flows: caustic formation, real-space singularities of the deformation tensor, and random uncorrelated motion. We discuss how the phase-and real-space singularities are related. Their formation is well understood in terms of a local theory. We discuss implications for random uncorrelated motion.

Introduction
The dynamics of particles suspended in randomly mixing or turbulent flows ('turbulent aerosols') has been studied intensively for several decades. Recently, substantial progress in understanding the dynamics of turbulent aerosols has been achieved (see the papers published in this special issue and the references cited therein).
The phenomenon of spatial clustering of independent point particles subject to Stokes drag in turbulent flows is now well understood: below the dissipative length scale (where the fluid flow is smooth) the particles eventually cluster onto a fractal set in configuration space. The corresponding fractal dimension has been determined by means of direct numerical simulations (Bec 2003) as well as theoretical approaches , Gustavsson & Mehlig 2011b. Different mechanisms ['preferential concentration' (Maxey 1987) and 'multiplicative amplification' , Gustavsson & Mehlig 2011b)] contribute to spatial clustering. A third mechanism giving rise to particle clustering was recently studied by following the deformation of an infinitesimally small volume of particles transported along a particle trajectory ['full Lagrangian method' (Ijzermans et al. 2010)]. The small volume may vanish at isolated singular points in time, giving rise to instantaneous singularities in the particle-concentration field. Using this approach the statistical properties of these singularities were analysed by Meneguz & Reeks (2011).
One important reason for studying spatial clustering of inertial particles is that this phenomenon is argued to enhance the rate at which collisions occur in turbulent aerosols at small values of the 'Stokes number'. This dimensionless parameter, St = (γτ ) −1 , is given in terms of the particle damping rate γ and the relevant correlation time τ of the flow. Both are defined below.
Arguably spatial clustering has an effect upon the collision rate at small Stokes numbers. But there is a second mechanism that leads to a significant enhancement of the collision rate as the Stokes number increases: direct numerical simulations of particles suspended in turbulent flows (Sundaram & Collins 1997, Wang et al. 2000 show that relative particle velocities at small separations increase substantially as the Stokes number is varied beyond a threshold of order unity. In (Falkovich et al. 2002, Wilkinson et al. 2006) this behaviour was explained by the occurence of singularities in the particle dynamics, causing large relative velocities at small separations. These singularities occur as the phase-space manifold folds, as illustrated in Fig. 1. As a consequence, particle velocities at a given point in space become multi-valued, causing large velocity differences between nearby particles. The boundaries of the folding region are referred to as 'caustics' , Crisanti et al. 1992. It was shown that the rate of caustic formation is an activated process , Gustavsson & Mehlig 2012. This explains the sensitive dependence of the rate of caustic formation upon the Stokes number observed in direct numerical simulations of particles in turbulence (Pumir & Falkovich 2007).
An alternative way of characterising relative velocities of inertial particles was suggested in (Fevrier et al. 2005, Simonin et al. 2006. The authors of these papers decomposed inertial particle velocities into two contributions: a spatially correlated, smoothly varying 'filtered' velocity field, and a random, spatially and temporally uncorrelated contribution, commonly referred to as 'random uncorrelated motion' (Reeks et al. 2006, Masi et al. 2011. The aim of this paper is twofold. First we summarise results of numerical simulations of particles suspended in model flows . Our numerical results for the moments of relative velocities of inertial particles are in quantitative agreement with recent analytical results based on the notion of caustic formation (Gustavsson & Mehlig 2011a). Second we demonstrate that caustic formation not only provides an understanding of relative velocities at small separations, caustic formation also explains spatial clustering due to singularities in the local deformation tensor, and the existence and properties of random uncorrelated motion.
We conclude the introduction by summarising our results in more detail. In this paper we show that recent predictions by Wilkinson et al. (2006) and Gustavson & Mehlig (2011a) based on the notion of caustic formation describe many aspects of the fluctuations of relative velocities at small separations. We compare formulae for the moments of relative velocities (Eqs. (18) and (20) below) to new results of numerical simulations of one-and two-dimensional models for inertial particles suspended in whitenoise flows, and for a three-dimensional kinematic simulation of particles suspended in an incompressible flow field with an energy spectrum typical of the small scales of Figure 2. (Left) multi-valued velocities of particles suspended in a two-dimensional random flow with finite Ku and St as described in Section 2.1. The base of each red arrow corresponds to a particle position (taken to be on a regular grid in the x-yplane). The orientation of the velocity is that of the arrow. All arrows have the same length, the magnitudes of the velocities are not shown. The blue line delineates the position of the caustics in the x-y-plane. The region of multi-valued velocities ends in a cusp that is only approximately resolved. In Section 4 it is explained how multivalued velocities between caustics give rise to so-called random uncorrelated motion. Parameters: Ku = 1, St = 10. (Right) particle-density in the x-y-plane, showing significantly enhanced particle-number density in the vicinity of the caustic line. Same parameters as above. Black corresponds to high density, white to low density. turbulence. We find good agreement. This demonstrates that Eqs. (18) and (20) which were derived in the white-noise limit, are valid more generally.
Further we examine the prediction by (Fevrier et al. 2005, Simonin et al. 2006) that the so-called longitudinal second-order structure function of relative velocities tends to a finite value at vanishing separations in the presence of random uncorrelated motion. The analytical theory, Eqs. (18) and (20) below, shows that this is true for sufficiently large Stokes numbers [the case examined numerically by (Simonin et al. 2006)]. But at Stokes numbers smaller than a critical value, the structure function tends to zero, despite the fact that there may still be a substantial singular (multi-valued) contribution to relative velocities due to the formation of caustics.
We discuss in detail that the singularities of the deformation tensor are in fact caustic singularities, as pointed out by ). We study the dynamics of the deformation tensor J, and the matrix Z of particle-velocity gradients. We show that as det J approaches zero, TrZ → −∞. We briefly remark upon the statistical properties of the singularities (Meneguz & Reeks 2011).
In summary, we demonstrate that the notion of random uncorrelated motion, and the occurrence of zeroes in the local deformation tensor can both be explained in terms of caustic formation, both qualitatively and in many ways quantitatively. Last but not least our results indicate that the white-noise approximation successfully describes many aspects of turbulent aerosols.
The remainder of this paper is organised as follows. In Section 2 we introduce the models analysed in this paper: inertial particles suspended in a two-dimensional incompressible random flow in the white-noise limit, and a kinematic simulation of inertial particle dynamics. Section 3 summarises what is known about the rate of caustic formation and discusses consequences for the fluctuations of relative particle velocities. We compare the analytical theory to results of numerical simulations of the models described in Section 2. In Section 4 we briefly review the notion of random uncorrelated motion, and compare the conclusions of (Fevrier et al. 2005, Simonin et al. 2006) to our analytical and numerical results. In Section 5 we describe the dynamics of the local deformation tensor and its correspondence to the dynamics of the matrix of particlevelocity gradients. Finally, Section 6 contains our conclusions.

Model
The motion of small, non-interacting spherical particles suspended in a flow is commonly approximated bẏ Here r and v are the position and velocity of a particle, u(r, t) is the velocity field evaluated at the particle position, γ is the viscous damping rate, and dots denote time derivatives. The components of the vector r are denoted by r j , j = 1, . . . , d in d dimensions. The components of u and v are referred to in an analogous way. Sometimes it is more convenient to denote the components of r by (x, y, z) instead of (r 1 , r 2 , r 3 ). We use the two notations interchangeably. For Eq.
(1) to be valid, it is assumed that the particle Reynolds number is small, that Brownian diffusion of the particles can be neglected, and that the particle density is much larger than that of the fluid. We also assume that the velocity field u varies smoothly on small spatial and temporal scales with smallest length-and time scales η and τ (the Kolmogorov scales for turbulent flows). The typical magnitude of the velocity field is denoted by u 0 .
In dimensionless units (t = t ′ /γ, r = ηr ′ , v = γηv ′ , u = γηu ′ and dropping the primes), the equation of motion becomeṡ The Stokes number does not appear explicitly in this equation, but the fluctuations of the dimensionless velocity u depend upon St (see Eq. (4) below). In addition to the Stokes number, the dynamics is characterised by a second dimensionless number, the 'Kubo number' Ku ≡ u 0 τ /η. We note that turbulent flows have Ku ∼ 1. In the remainder of this paper we frequently refer to these two dimensionless numbers. For a discussion of further dimensionless parameters see . The numerical results shown in the following were obtained for two different models. These models are introduced in the following two subsections.

Random-flow model
Following (Wilkinson & Mehlig 2003) we approximate the incompressible velocity field u(r, t) in Eq.
(2) by a Gaussian random function that varies smoothly in space and time. We discuss results for one-and two-dimensional versions of the random-flow model. The one-dimensional case is most easily analysed, the two-dimensional incompressible case is important (since one-dimensional flows are special, they are always compressible which give rise to a path-coalescence transition (Wilkinson & Mehlig 2003)). A twodimensional incompressible velocity field can be written in terms of a stream function ψ(r, t): Here e 3 is the unit vector ⊥ to the x-y-plane. We assume that ψ(r, t) is a Gaussian random function with ψ = 0 and correlation function in dimensionless variables.
In this paper we also refer to results of a one-dimensional random-flow model. This is defined in an analogous fashion in terms of a Gaussian random flow velocity u(x, t) with zero mean and correlation function We note the one-dimensional flow is compressible. The numerical data shown in Figs. 1 and 2 are obtained by computer simulations of the models described above. We simplify the model by linearising Eq.
(2). This yields the following equation for the dynamics of a small separation R = r 1 − r 2 and velocity difference V = v 1 − v 2 between two particles: Here A is the matrix of fluid velocity gradients, with elements A ij = ∂u i /∂r j .
To simplify further, we take the white-noise limit of this model. This limit corresponds to Here ǫ is a dimensionless measure of the particle inertia introduced by (Mehlig & Wilkinson 2004) [see also ]. We take c 1 = 1 for one-dimensional flows [this is consistent with the convention used in (Gustavsson & Mehlig 2011a)].
For incompressible two-dimensional flows we take c 2 = 1/2, as in (Gustavsson & Mehlig 2011b). In the white-noise limit, the instantaneous value of the velocity gradient A in (6) becomes independent of the particle position. In two spatial dimensions, we denote the independent random increments of the elements A 11 , A 12 and A 21 of A in a small time step δt by δa 1 , δa 2 , and δa 3 . Note that A 22 = −A 11 since the flow is incompressible. One finds: The results shown in Figs. 3, 4, and 5 are obtained by computer simulations of this model, approximating the time-dependence of A(r(t), t) as a white-noise signal.

Kinematic simulation
As an alternative to the single-scale white-noise model introduced in the previous subsection, we simulate a turbulent incompressible velocity field in a three-dimensional periodic box by a large number of Fourier modes varying randomly in space and time.
The modes are chosen in such a way that the associated energy spectrum approximates a prescribed form, namely that originally used by (Kraichnan 1970). The model is identical to that used by (Ijzermans et al. 2010, Meneguz & Reeks 2011. For convenience we briefly summarise its relevant features below. For details, we refer the reader to (Ijzermans et al. 2010, Meneguz & Reeks 2011. In dimensionless form, the incompressible velocity field u(r, t) is represented as a Fourier series of N modes (N = 200 in our simulations): with random coefficients a (n) and b (n) , random wave numbers k (n) , and random frequencies ω (n) . In order to guarantee the periodicity of the flow in a cube of dimensions L × L × L, the allowed wave number components k (n) 2π is the integral length scale of the flow. The integer numbers m (n) i are chosen randomly in such a way that the lengths k (n) = k (n) · k (n) are approximately equal to the ideal wave number k (n) id . The latter is determined by the energy spectrum as follows: As mentioned above, the energy spectrum E(k) is taken to be (Kraichnan 1970): This spectrum is representative for low-Reynolds-number turbulence (Spelt & Biesheuvel 1997). The maximum of E(k) is located at k = 1 and the total kinetic energy ∞ 0 E(k)dk = 3/2. This corresponds to 3u 2 0 /2 in dimensional form. The use of the Kraichnan energy spectrum results in a relatively small separation of scales; in our simulations, the smallest wave number k (1) ≃ 0.25 and the largest wave number k (N ) ≃ 2.14. The frequencies ω (n) are chosen randomly from a Gaussian distribution with zero mean and a variance proportional to k (n) . This implies that the Kubo number is of order unity. Following (Spelt & Biesheuvel 1997), we take the variance to be 0.4 k (n) . Finally, the coefficients a (n) and b (n) are determined by choosing a random direction in Cartesian space, and by picking a length randomly from a Gaussian distribution with zero mean and a variance 9/(2N). By doing so, the mean kinetic energy at a given position in spacē is approximately equal to 3/2 for all values of r.

One spatial dimension
As illustrated in Fig. 1, caustics form when the phase-space manifold folds over. In one spatial dimension this happens when the slope of the manifold becomes infinite, that is when z = ∂v/∂x → −∞. The rate at which this occurs is determined by the equation of motion for z (Wilkinson & Mehlig 2003): Here A = ∂u/∂x represents the random driving by the fluid-velocity gradients. In the case of independent particles (which we consider here), z goes through infinity in a symmetrical fashion. At large values of |z|, the random driving can be neglected, so thatż ≈ −z − z 2 . The corresponding deterministic probability distribution of z reads ρ(z) = C/[z(1 + z)], and is valid in the tails of z. In the white-noise limit, Eq. (15) is equivalent to a Fokker-Planck equation for the distribution of z. In (Wilkinson & Mehlig 2003) this equation was solved in one spatial dimension. The resulting rate of caustic formation [called 'rate of crossing caustics' by Wilkinson & Mehlig (2003)] can be written as (Gustavsson & Mehlig 2012): where ǫ 2 = Ku 2 St (see Section 2.1). In Eq. (16), Ai(y) is the Airy function. In the limit of small values of ǫ, this expression exhibits the asymptotic behaviour Eq. (16) shows that the number of caustics increases rapidly as ǫ 2 passes through 1/6 , Wilkinson et al. 2006. This sensitive dependence is commonly referred to as an 'activated law', in analogy with the sensitive temperature depdence of chemcial reaction rates in Arrhenius' law. Gustavsson & Mehlig (2012) computed the one-dimensional rate of caustic formation at small but finite Kubo numbers and found it to sensitively depend on the Stokes number: in this case too, the St-dependence exhibits 'activated form': where S is an St-dependent 'action'. In the white-noise limit, S = 1/(6St), consistent with Eq. (17). As Fig. 1 shows, particle-velocities become multi-valued between two caustics in the wake of a singularity, giving rise to large relative velocities between nearby particles. While the rate of caustic formation is determined by the rate at which the local quantity z = ∂v/∂x tends to −∞, the distribution of relative velocities at small particle separations is determined by the solution of the full non-local equations (6) for particle separations and relative velocities (Gustavsson & Mehlig 2011a).
A consequence of large relative velocities at small separations is that between caustics, particles collide frequently with large relative velocities (c.f. Fig. 1), giving rise to a large collision rate (we note, however, that in this paper it is assumed that the particles are independent point particles that do not actually collide).
By contrast, in the absence of caustics, particles may still approach each other due to fluctuations of the underlying flow-velocity field. At small separations the flow is smooth, and in this regime relative velocities between particles are expected to tend to zero as the particles in question approach each other.
Which one of these two mechanisms of bringing particles together makes the dominant contribution to the collision rate depends upon the value of St and on the particle size a (separation 2a at the point of contact). Relative velocities of particles thrown at each other due to the formation of caustics are expected to make the dominant contribution to the collision rate if St is large and/or when the particles are sufficiently small. Particles slowly approaching each other ('logarithmic diffusion') dominate otherwise.
In the white-noise limit, and in one spatial dimension, an asymptotic approximation for the moments of relative velocities at small separations was derived in (Gustavsson & Mehlig 2011a): Here X = x 1 − x 2 and V = v 1 − v 2 are the separation and the relative velocity of a pair of particles, and ρ(X, V ) is their distribution function. It is assumed that |X| ≪ 1 and p > −1. Further, D 2 is the correlation dimension of the phase-space The second term in Eq. (18), C p , is due to multi-valued velocities between caustics. This contribution, in one spatial dimension, does not depend upon |X| for small values of |X|. In other words: it remains finite as |X| → 0. This is a consequence of the fact that as the manifold in Fig. 1 folds over, particles initially far apart are thrown at each other quickly.
The first term in Eq. (18), B p |X| p+D 2 −1 , vanishes as |X| → 0. It constitutes the main contribution to m p (X) in the absence of caustics and is affected by spatial clustering: for a given value of p, the exponent is smallest (and thus the contribution largest) when D 2 attains its minimum as a function of Stokes number.
In Eq. (18) the case p = 1 is of particular importance, since m 1 (X) is closely related (yet not identical) to the collision rate between particles at small separations X = 2a. It is expected that the coefficient of the caustic contribution in Eq. (18), C 1 , is proportional to the caustic formation rate J caustic (Wilkinson et al. 2006).

Two and three spatial dimensions
In two and three spatial dimensions the caustic rate can be found in a way similar to the one-dimensional case , Gustavsson & Mehlig 2011b: the matrix Z with elements Z ij = ∂v i /∂r j obeys the equation: Here A is the matrix of fluid-velocity gradients introduced in Section 2, with elements A ij = ∂u i /∂r j . In analogy with the one-dimensional case, tr(Z) → −∞ as caustics are formed. In the white-noise limit, we expect that the rate of caustic formation is again given by (17). In ) numerical factors in Eq. (17) slightly different from 1/6 were quoted in two and three spatial dimensions. More recent numerical results (not shown) show that the asymptote (17) is approached very slowly as ǫ becomes small. Our best estimates at the smallest values of ǫ indicate that the factor in the argument of the exponential in (17) is asymptotically the same (equal to 1/6) in one, two, and three spatial dimensions.

Moments of relative velocities in two and three spatial dimensions obey laws analogous to (18). At small separations
Here R = |R| and v R ≡ V ·ê R is the radial projection of the relative velocity between two particles at separation R. Further, d 2 is the spatial correlation dimension, it is assumed that the Stokes number is small enough so that d 2 ≤ d. As in the one-dimensional result, Eq. (18), there are two contributions to the moments of relative velocities [compare the parameterisation of the St-dependence of the collision rate suggested by Wilkinson et al. (2006)]. The second term in Eq. (20) is due to multi-valued velocities between caustics. But note that in two and three spatial dimensions, not all particle pairs thrown together give rise to close approaches. The reason is that in addition to having one relative coordinate pass zero at finite relative velocity (so that a caustic occurs), the other coordinates must be small, i.e. only particles heading sufficiently towards each other as the caustic occurs end up at small enough separations to contribute to the small R velocity moments. This explains the geometrical factor R d−1 in (20). It is absent in one spatial dimension, d = 1. Fig. 3 shows comparisons of Eq. (20) with results of numerical simulations of the random-flow model described in Section 2.1. The parameter d 2 in Eq. (20) is determined as follows. Setting p = 0 in (20) and taking the limit R → 0 defines the spatial correlation dimension d 2 . The latter is found numerically by fitting m 0 (R) to the power law R d 2 −1 .
We now describe how the fits in Fig. 3 were obtained. The parameter d 2 was taken from Fig. 4. The coefficients B p and C p in Eq. (20) were fitted to the numerical results for different parameter values. The fitting region (the range of R over which Eq. (20) is fitted) lies between the dashed lines in Fig. 3. We observe good agreement between the numerical results and fits to Eq. (20). In particular, the results clearly show that the moments m p scale as R d−1 for small values of R, independently of p. Fig. 5 shows the coefficient C 1 of the caustic contribution obtained in this way as a function of ǫ −2 . Since this contribution requires the formation of caustics, we expect C 1 to exhibit an ǫ-dependence of the form (17). Fig. 5 shows that this is indeed the case. Fig. 6 shows results for m 0 (R) and m 1 (R) obtained by kinematic simulations of the random-flow model described in Section 2.2 for two values of the Stokes number, St = 0.4 and St = 0.7. As expected, m 0 (R) is of power-law form, reflecting spatial clustering. The corresponding correlation dimensions are shown, as a function of St, in Fig. 7. The correlation dimension exhibits the expected minimum (here at St ≈ 0.4). Corresponding results for direct numerical simulations of particles in turbulent flows have been obtained by a number of authors (Wilkinson et al. 2010, Bec et al. 2010, Chun et al. 2005. The green squares in Fig. 6 correspond to numerical results for m 1 (R) as a function of R. Consider first the left panel (St = 0.4). At small separations R we expect that m 1 (R) should scale as R d−1 = R 2 , while it should scale as R d 2 ≈ R 2.4 at large values of R. Despite the fact that the two powers are rather similar, the two scalings can be distinguished in Fig. 6. In the right panel (St = 0.7), the caustic contribution R d−1  dominates. Given the data available from the kinematic simulations, it is more difficult to reliably determine the St-dependence of C 1 by fitting (solid green line in Fig. 6). Our best estimates are shown in Fig. 8. The fits and the corresponding error bars were obtained by a non-linear least-squared fit using MATLAB 2011. We find that C 1 depends very sensitively on St, as expected because the formation of caustics is an activated process. We expect (Gustavsson & Mehlig 2011a, Gustavsson & Mehlig 2012 that the St-dependence of C 1 follows the law J caustic /γ ∼ exp[−S(St)/Ku 2 ]. However, the range of Stokes numbers for which C 1 can reliably be estimated is too small to determine the form of the function S(St).

Random uncorrelated motion
Singularities in the inertial-particle dynamics (corresponding to the formation of caustics) give rise to multi-valued particle velocities at locations in space bounded by caustics: any identical particles that are very close may move at substantially different velocities. Figs. 1 and 2 illustrate this fact in one and two spatial dimensions. This implies in particular that the relative motion of inertial particles cannot be captured in terms of a 'hydrodynamic' approximation describing the particle velocities in terms of a smooth velocity field. In particular Fevrier et al. (2005) infer from their DNS calculations of inertial-particle motion in a homogeneous isotropic and stationary turbulent flow field that the velocity of a particle at a position r(t) at time t in a single realisation of the carrier flow field u(r, t) is given by the sum of two components v(r(t) = r, t) = v(r, t) + δv(r(t) = r, t) .
Here v(r, t) is a smoothly varying filtered velocity field which Fevrier et al. (2005) refer to as the 'mesoscopic Eulerian particle velocity field'. Values for the smooth component in any realisation are found by dividing the spatial domain into cells and calculating the average velocity associated with the number of particles in each individual cell (the number of particles in each cell being sufficiently large to form a statistically stationary average). The residual component δv is termed the 'quasi Brownian velocity distribution component' by Fevrier et al. (2005). It is now commonly referred to as 'random uncorrelated motion', or RUM for short (Reeks et al. 2006, Masi et al. 2011). This residual RUM part is assumed to be uncorrelated with the smooth part and with itself at infinitesimally small separations in space and time.
The existence of multi-valued velocities between caustics is consistent with a singular contribution to the particle velocities, of the form of Eq. (21). We infer that the extent of random uncorrelated motion [its relative contribution compared to the smooth part in Eq. (21)] must depend sensitively on the value of St, since the rate of caustic formation exhibits this sensitive dependence on the Stokes number.
Let us consider the implications of Eq. (21) and the accompanying assumptions for the second moment of the relative radial velocity between two particles, v This result is of the same form as Eq. (20). The two right-most terms in (22) correspond to the caustic contribution in (20). In other words, Eq. (20) provides a quantitative prediction for the contribution of random uncorrelated motion to the moments of relative radial velocities. Consider for example the form of the so-called 'longitudinal structure functions' for relative velocities of the suspended particles. Simonin et al. (2006) argue that the second-order structure function remains finite as the spatial separation R between particle velocities tends to zero. In the notation of the previous section, the second-order structure function is given by The limiting behaviour of s (2) (R) can be deduced from Eq. (20). From this equation we see that m 0 (R) ∼ R min{d 2 ,d}−1 . The correlation dimension d 2 saturates to d at a critical Stokes number, St c (c.f. Fig. 4 where d 2 = d for ǫ 2 > ǫ 2 c ≈ 1). For St > St c the suspended particles are uniformly distributed in space [see also (Bec et al. 2010, Salazar & Collins 2012]. Let us consider this case. As R → 0, the caustic contribution C 2 R d−1 to m 2 (R) dominates in Eq. (20). This implies that as argued in (Simonin et al. 2006). For St < St c , by contrast, we find More precisely, s (2) (R) tends to zero as g −1 (R) when R → 0 (the pair correlation function g(R) is given by g(R) = m 0 (R)/R d−1 ). We emphasise that the behaviour (24) of the structure function in two and three spatial dimensions must be distinguished from the fact that the moments m p (X) of relative velocities in one spatial dimension always approach a positive constant as |X| → 0 when St > 0. Indeed, we have shown that s (2) (R) may approach zero as R → 0, yet multi-valued particle velocities do still give rise to a substantial singular contribution to the moments of relative velocities, as a consequence of singularities giving rise to caustics.
Let us compare these findings to the results shown in Fig. 3(a) in Simonin et al. (2006). The data shown in this figure (except perhaps the data set labeled '1') imply that the structure function approaches a positive constant as R → 0. We conclude that the data sets shown (possibly with the exception of '1') correspond to Stokes numbers larger than St c . It should be noted that Simonin et al. (2006) define their Stokes number St L in terms of the integral time scale of the turbulent flow. Here and in a large part of the literature on inertial particles in turbulent flows the Stokes number St is defined in terms of the Kolmogorov time τ . Since usually St ≫ St L it is plausible that most data sets in Fig. 3(a) in Simonin et al. (2006) We conclude by noting that it has been shown (see ) and references cited therein) that the maximal Lyapunov exponent describing the dynamics of inertial particles suspended in incompressible flows is positive. This implies that the inertial particle dynamics is chaotic. In the limit of very large Stokes numbers, inertial particle dynamics is thus similar to the random motion of molecules in a gas [gas-kinetic limit, see (Abrahamson 1975)]. This justifies the view that there is a random uncorrelated component to the inertial particle dynamics. It is a consequence of the formation of caustics.

Singularities in particle concentration
Changes to the local concentration of inertial particles suspended in mixing flows can be described by the deformation tensor J with elements J ij = ∂r i /∂r j (0) evaluated along a particle trajectory r(t) with initial position r(0). The matrix J describes the relative motion of infinitesimally close particles. In particular, the volume spanned by the separation vectors between d+1 infinitesimally close particles in d spatial dimensions is given by δV = |J|δV 0 , where J ≡ det(J) and δV 0 is the initial volume, see Fig. 9.
Nothing prevents J from occasionally changing sign. This implies that the volume δV may shrink to zero, giving rise to a singularity in the local particle concentration ∝ δV −1 , Ijzermans et al. 2010. The singularities influence the tails of the distribution of local particle concentration, making particle clustering highly non-Gaussian and intermittent (Meneguz & Reeks 2011). The zeroes of J correspond to the formation of caustics ). This fact is illustrated in Fig. 1: as J → 0 we see that z → −∞. In the following we discuss the dynamics of z and J in one spatial dimension, and then the dynamics of Z and J in two and three spatial dimensions. Singularities in the local particle density due to caustics occur also in a collisionless medium of weakly interacting particles. As a model for the early structure of the universe, the corresponding linear equation of motion r(t) = r 0 + t v(r 0 ) has been analysed by Zeldovich and collaborators. For a review and a discussion of the connection between this problem and Burgers' equation see (Shandarin & Zeldovich 1989).

One spatial dimension
In one spatial dimension we analyse the joint dynamics of z = ∂v/∂x and J = ∂x/∂x 0 , where x 0 = x(0) is the initial particle position. Noting thatJ = ∂v/∂x 0 we see that Figure 9. Illustrates how an infinitesimal area element δA(t) in two spatial dimensions spanned by the separation vectors δR 1 and δR 2 between three initially close particles is transported along the particle trajectories.
z =J/J. The dynamics of z is governed by Eq. (15) which in turn yields an equation for the dynamics of J: with A = ∂u/∂x. This is the one-dimensional analogue of Eq. (2.20) in (Ijzermans et al. 2010). The singularities z → −∞ and J → 0 occur simultaneously. This can be seen in the deterministic limits of Eqs. (15) and (26): assume that z is large. Then (15) can be approximated byż = −z − z 2 . When J is small then (26) is approximatelyJ = −J . These two equations are solved by Consider an initial condition z 0 < −1. In this case, singularities in z and J occur as t passes through t 0 = ln(z 0 /(1 + z 0 )) for both solutions (27). Thus, the rate at which J passes 0 is identical to the rate at which z tends to −∞.

Two and three spatial dimensions
In two and three spatial dimensions the situation is analogous. The matrices Z and J are related by Z =JJ −1 . Eq. (19) gives the motion of Z and the corresponding equation This equation is identical to Eq. (2.20) in (Ijzermans et al. 2010). In analogy with the one-dimensional case, the deterministic solution is found to be: where Z =JJ −1 is obtained from (29). Singularities occur when the determinant J ≡ det(J) vanishes, or equivalently when TrZ =J /J diverges. The determinant of J is obtained from (29) in two and three spatial dimensions J d=3 = J 0 [1 + T 1 + T 2 + Z 0 − e −t (T 1 + 2T 2 + 3Z 0 ) where the invariants J 0 ≡ det J 0 , Z 0 ≡ det Z 0 , T 1 ≡ TrZ 0 and T 2 ≡ [(TrZ 0 ) 2 − Tr(Z 2 0 )]/2 were defined. Depending on the initial condition Z 0 =J 0 J −1 0 , J may pass zero at a finite time t 0 . NowJ and J cannot pass zero simultaneously (assuming that J(t) is a regular function, thenJ(t 0 ) = 0 implies that J(t) has a double root at t 0 ). It follows that TrZ is singular at t 0 . We have explicitly checked in two spatial dimensions that this is the case.

Conclusions
In this paper we have compared three recent approaches to describing inertial particle dynamics: caustic formation giving rise to multi-valued particle velocities, the notion of random uncorrelated motion, and spatial clustering as a consequence of singularities in the local deformation tensor J.
We have shown that clustering due to singularities of J can be explained in terms of caustic formation. Furthermore we have compared the consequences of the hypothesis of random uncorrelated motion with predictions for the fluctuations of relative velocities in random-flow models. The hypothesis of random uncorrelated motion leads to an expression for the moments of relative velocities that consists of two terms: a smooth part, and a contribution due to random uncorrelated motion. This expression corresponds precisely to Eqs. (18) and (20) for the moments of relative velocities obtained in (Gustavsson & Mehlig 2011a). These theoretical results, describing the effect of caustics upon the fluctuations of relative velocities, make it possible to quantify the degree of random uncorrelated motion, commonly measured in terms of the longitudinal structure function s (2) (R): for Stokes numbers below a critical value, s (2) (R) tends to zero as the separation R → 0.
We have performed numerical simulations of one-and two-dimensional random-flow models in the white-noise limit as well as kinematic simulations at finite Kubo numbers. We have found that results of these simulations are consistent with Eqs. (18) and (20).
Recently, two comprehensive studies of inertial particle dynamics using direct numerical simulations of particles suspended in turbulent flows were published (Bec et al. 2010, Salazar & Collins 2012. A detailed comparison between the analytical theory and the results of these direct numerical simulations for the distribution and the moments of relative velocities will be published elsewhere (Cencini et al. 2012).
Last but not least, we remark that the phenomenon of clustering and relative particle dynamics in turbulent flows analysed here has much in common with the way particles are transported and deposited in turbulent boundary layers (Young & Leeming 1997): enhanced particle concentrations are observed near the wall, corresponding to the clustering of inertial particles in turbulent flows. Moreover, as in the case of particles suspended in turbulent flows, particle inertia gives rise to large impact velocities [referred to as 'free flight to the wall' (Brooke et al. 1994)].