Lyapunov exponents for particles advected in compressible random velocity fields at small and large Kubo numbers

We calculate the Lyapunov exponents describing spatial clustering of particles advected in one- and two-dimensional random velocity fields at finite Kubo numbers Ku (a dimensionless parameter characterising the correlation time of the velocity field). In one dimension we obtain accurate results up to Ku ~ 1 by resummation of a perturbation expansion in Ku. At large Kubo numbers we compute the Lyapunov exponent by taking into account the fact that the particles follow the minima of the potential function corresponding to the velocity field. The Lyapunov exponent is always negative. In two spatial dimensions the sign of the maximal Lyapunov exponent \lambda_1 may change, depending upon the degree of compressibility of the flow and the Kubo number. For small Kubo numbers we compute the first four non-vanishing terms in the small-Ku expansion of the Lyapunov exponents. By resumming these expansions we obtain a precise estimate of the location of the path-coalescence transition (where \lambda_1 changes sign) for Kubo numbers up to approximately Ku = 0.5. For large Kubo numbers we estimate the Lyapunov exponents for a partially compressible velocity field by assuming that the particles sample those stagnation points of the velocity field that have a negative real part of the maximal eigenvalue of the matrix of flow-velocity gradients.


Introduction
Consider many small tracer particles advected in a random or chaotic compressible flow. An initially uniform scatter of particles cannot remain uniform because particles advected in smooth compressible flows cluster together. An example of this effect is discussed by Sommerer and Ott [1] who describe experiments following fluorescent tracers floating on the surface of an unsteady flow. Since the particles are constrained to the surface of the flow, they experience local up-and down-welling regions as sources and sinks, rendering the surface flow compressible. As a consequence the particles form fractal patterns. The authors of [1] interpret these patterns in terms of random dynamical maps and estimate the Lyapunov fractal dimension. This dimension is computed from the Lyapunov exponents by means of the Kaplan-Yorke formula [2]. The Lyapunov exponents λ 1 > . . . > λ d (here d is the spatial dimension) describe the long-term evolution of the patterns formed by the particles. The maximal Lyapunov exponent λ 1 describes the dynamics of an initially infinitesimal separation between two particles. When λ 1 < 0 separations between nearby particles must decrease on the long run, clustering is strong. This regime was referred to as 'path-coalescence phase' in [3]. The path-coalescence transition occurs at λ 1 = 0: when λ 1 > 0 separations typically grow, but clustering can nevertheless be substantial [4]. The sum λ 1 +λ 2 describes the evolution of a small area element spanned by the separation vectors of three nearby particles, and so forth.
Refs. [5,6] summarise results of direct numerical simulations of tracers floating on the surface of turbulent flows, and characterise the resulting fractal patterns in terms of their Lyapunov dimensions.
Another example is that of inertial particles in turbulent flows. Finite inertia allows particles to detach from the flow. This effect gives rise to fractal patterns of inertial particles suspended in incompressible flows [7,8]. When the particle inertia is small it is commonly argued that the resulting fractal patterns can be understood in terms of a model that describes particles advected in a slightly compressible particle-velocity field [9]. The small correction term that renders the particle-velocity field compressible at small particle inertia was first derived by Maxey [10].
This approach is frequently used in the literature to explain spatial clustering (so-called 'preferential concentration') of inertial particles suspended in turbulent flows. We remark that this approach must fail when inertial effects become stronger. In this case the particle-velocity field develops singularities (so-called 'caustics') that preclude the existence of a smooth particle-velocity field (see e.g. [11,12,13,14]). The singularities give rise to large relative velocities between nearby particles [15,16,17].
Many authors have studied the Lyapunov exponents of particles advected in turbulent, random, and chaotic velocity fields numerically. Analytical results could only be derived in certain limiting cases though. A limit that allows analytical progress in terms of diffusion approximations is Ku → 0. The 'Kubo number' Ku = u 0 τ /η is a dimensionless measure of the correlation time τ of the fluctuations of the underlying velocity field, u 0 is the typical speed of the flow and η its correlation length. In the limit of Ku → 0 3 the problem of calculating the Lyapunov exponents simplifies considerably. In this limit the flow causes many weakly correlated small displacements of the particles and diffusion approximations can be used to compute the exponents for random Gaussian flows [18,8] and in the Kraichnan model [19,20,21]. The exponents describe the fluctuations of small separations between particles (much smaller than the correlation length or the Kolmogorov length η), inertial-range fluctuations are not relevant in this limit, and thus the results for smooth random velocity fields and for the Kraichnan model are equivalent when Ku → 0. Lyapunov exponents for inertial particles in the limit of Ku → 0 were computed in Refs. [3,22,23] in one, two, and three spatial dimensions respectively.
Less is known about clustering at finite Kubo numbers where the particles have sufficient time to preferentially sample the sinks of the underlying velocity field. This effect is important in the examples mentioned above, but it is not captured by theories formulated in terms of diffusion approximations in the limit of Ku → 0. At large Kubo numbers the spatial patterns formed by the particles must depend on the details of the fluctuations of the underlying flow, but it is not known how to analytically compute the Lyapunov exponents of particles advected in compressible velocity fields at finite Kubo numbers. We note however that an expression for the maximal Lyapunov exponent in incompressible two-dimensional flows at finite Kubo numbers was obtained by Chertkov et al. [24]. Approximating the fluctuations of the flowvelocity gradient by telegraph noise with a finite correlation time, Falkovich et al. computed Lyapunov exponents in one-dimensional and incompressible two-dimensional models for advected and inertial particles [25,26]. Also, Dhanagare et al. [27] recently investigated the spatial clustering of particles advected in compressible random renovating flows.
In this paper we compute the Lyapunov exponents for particles advected in one-and two-dimensional compressible Gaussian random velocity fields with finite Kubo numbers (the model is defined in Section 2). We use an approach recently developed to describe incompressible turbulent aerosols at finite Kubo numbers [28], generalising a method used by Wilkinson [29] to compute the Lyapunov exponent for particles advected in a one-dimensional random velocity field to lowest order in Ku. This approach expresses the fluctuations of the flow-velocity gradient along the particle trajectories at finite Kubo numbers in terms of correlation functions of the flow velocity and its derivatives at fixed positions in space. A perturbation expansion in Ku is obtained by iteratively refining an approximation for the paths taken by the particles [28].
In Section 3 we develop perturbation series to order Ku 12 in one spatial dimension. By comparison with computer simulations we show that a Padé-Borel resummation of the series yields accurate results up to Ku ∼ 1. For Ku ≫ 1 the resummation fails. In this case the particles are predominantly found near stagnation points of the flow (where the flow velocity vanishes) with negative velocity gradients, that is near the minima of the corresponding potential function. In this regime the Lyapunov exponent is determined by the flow-gradient fluctuations near these points, and we show how to compute Table 1 Conversion table comparing different parametrisations of compressibility in d-dimensional random velocity fields. The parameters β and Γ were introduced in Ref. [22] in two and in Ref. [32] in three spatial dimensions. The parameter ∆C in Eq. (6.21) in Ref. [33], and ℘ in Eq. (57) in Ref. [20] (see also Ref. [21]).
the exponent using the Kac-Rice formula [30,31] for counting singular points of random functions. Section 4 summarises the corresponding results for two-dimensional velocity fields. For small Kubo numbers we compute the first four non-vanishing terms in a perturbation expansion (to order Ku 8 ). We compare the results of a Padé-Borel resummation of this series with results of numerical simulations. We find that resummation of the perturbation series provides an accurate estimate of the location of the path-coalescence transition (the degree of compressibility where the maximal Lyapunov exponent λ 1 changes sign) for Kubo numbers up to ∼ 0.5.
For much larger Kubo numbers, particles in a purely compressible velocity field (that can be written as the gradient of a potential function) spend most of their time near the minima of the potential function, that is near stagnation points of the velocity field with negative real part of the maximal eigenvalue of the matrix of flow-velocity gradients. As in the one-dimensional case the Lyapunov exponents can be computed using the Kac-Rice formula. Our results agree well with those of numerical simulations of particles in velocity fields with a compressible component, at large but finite Kubo numbers.
Section 5 summarises our conclusions.

Model
We study particles advected in one-and two-dimensional Gaussian random velocity fields. In one dimension the equation of motion iṡ Here x t denotes the position of the particle at time t, the dot denotes a time derivative, and u(x, t) is a Gaussian random velocity field. We write u = u 0 ∂ψ/∂x where u 0 is the typical speed of the flow, and ψ(x, t) is a Gaussian random function with zero mean values and correlation function The correlation length is denoted by η, and the correlation time is denoted by τ . In two spatial dimensions we writė The velocity field is defined as [22]: whereê z is the unit vector in the z-direction and where ψ and φ are independent Gaussian random functions with zero means and correlation functions The first term in Eq. (4) is an incompressible (or 'solenoidal') contribution. The second term is a compressible (or 'potential') contribution. As in one dimension the speed-, length-, and time scales of the flow are denoted by u 0 , η and τ . The Kubo number is given by Ku = u 0 τ /η. In the following we adopt dimensionless units t = τ t ′ , x = ηx ′ , u = u 0 u ′ and we drop the primes. A second dimensionless parameter of the problem is the degree of compressibility. Ref. [22] introduced the parameter It ranges from 1/3 (β → ∞, compressible) to Γ = 3 (β = 0, incompressible).
In the limit of Ku → 0, the maximal Lyapunov exponent is negative for Γ ≤ 1 (β ≥ 1) and positive otherwise. Other authors parametrise the degree of compressibility in other ways. Table 1 compares different definitions.

One spatial dimension
The Lyapunov exponent describes the long-term growth (or decline) of the separation δx t between two initially infinitesimally close particles. It is computed by linearising the equation of motion (1) to find the dynamics of a small separation δx t between two neighbouring particles. Using the dimensionless variables introduced in Section 2 we have: Here the flow-velocity gradient at position x at time t is denoted by A(x, t). It follows that the Lyapunov exponent is given by by the average flow-velocity gradient evaluated at the particle position x t : 6 where · · · denotes an average over flow realisations. The Lyapunov exponent is computed by expanding the implicit solution of (1). In dimensionless units it is given by: Here x 0 is the initial particle position, and ξ t = x t − x 0 is the difference between the trajectory of a particle and its initial position (to be distinguished from the separation δx t between two neighbouring particles at time t). Since ξ t is proportional to Ku it can be considered small provided that Ku is sufficiently small. In this case we expand u(x t , t) in powers of ξ t : Inserting ξ t = x t − x 0 from Eq. (10) into Eq. (11) and iterating Eq. (11) yields a perturbation series of u(x t , t) in terms of powers of Ku. In the same way an expansion of A(x t , t) is found. To third order in Ku we obtain for example: where , and so forth. Averaging yields an expression for the Lyapunov exponent in terms of time integrals of Eulerian correlation functions of the velocity field and its derivatives. Evaluating these correlation functions requires computing averages of products of u(t) and its spatial derivatives ∂ k x u(t) ≡ ∂ k x u(x, t)| x=x0 (for k = 1, 2, . . .) evaluated at different times. For a Gaussian random velocity field we use Wick's theorem which states that the average of a product of n Gaussian variables z 1 , . . . z n is equal to the sum of all ways of decomposing the product into a products of covariances. The averaged product z 1 , . . . z n is calculated using the known covariances z i z j . Averages of products of an odd number of factors vanish when z i = 0. In one spatial dimension, the covariances of the velocity and its spatial derivatives are determined by Eq. (2): with m, n = 0, 1, 2, . . . . In this way we obtain an expansion of the Lyapunov exponent in powers of Ku. The final result to order Ku 12 is: This series expansion is asymptotically divergent: it diverges for any fixed value of Ku but every partial sum of the series approaches λ as Ku → 0. At large orders k the coefficients c k in the series (14) are of the form [34] c with S close to 1/6, a ≈ 0.23, and b ≈ 1. The series (14) can be resummed by Padé-Borel resummation [34]. The result is expressed as the Laplace transform of the so-called 'Borel sum' (assumed to have a finite radius of convergence due to the extra factor of 1/k!): The Lyapunov exponent is estimated by The integration path C is taken to be a ray in the upper right quadrant of the complex plane. In order to compute the integral the Borel sum must be analytically continued outside its radius of convergence. This can be achieved by 'Padé approximants' [35]. We know 6 non-zero coefficients in the sum, this allows us to compute the Padé approximant of third orders in Ku 2 in numerator and denominator: Eqs. (17) and (18) together determine an approximation for the Lyapunov exponent. The corresponding result is shown in Fig. 1, compared with results of numerical simulations of the model. We observe good agreement for Kubo numbers up to order unity. If more coefficients in the perturbation series were known, higher-order Padé approximants could be computed to improve the accuracy of the Padé-Borel resummation.
For large values of Ku the resummation fails. We now show how to approximate the Lyapunov exponent at large but finite values of Ku (the limit t → ∞ in Eq. (7) is taken at a finite value of Ku). At large Kubo numbers the particles spend most of their time in the vicinity of the minima of the 'potential' V (x, t) = −ψ(x, t). An approximate expression for the Lyapunov exponent can be obtained by averaging the gradient A at the minima, that is at the stagnation points of the flow velocity u with A < 0. The distribution of A at u = 0 can be estimated using the Kac-Rice formula [30,31]: Here P [u, A] is the joint distribution of the velocity field u and its gradient A. It is determined by the correlation function of ψ given in Section 2. We find: Since p is symmetric in A, the distribution of negative values of A at u = 0 is given by 2 p(A). The Lyapunov exponent is thus given by The limiting behaviour (21) is shown in Fig. 1. It is in good agreement with the results of numerical simulations. 9 4 Two spatial dimensions

Small-Ku limit
Consider first the case of small Kubo numbers. Now there are two Lyapunov exponents to compute, describing the time evolution of the distance |δx t | between two neighbouring particles and of the infinitesimal area element δA t spanned by the separation vectors between three neighbouring particles: As in the one-dimensional case, these Lyapunov exponents are computed by linearising the equation of motion, Eq. (3) in this case. The dynamics at small separations δx between two neighbouring particles is Here A is the flow-gradient matrix with elements A ij ≡ ∂u i /∂x j . The Lyapunov exponents (23) are calculated from These relations are analogous to the one-dimensional Eq. (9). In Eq. (25), the vector n t ≡ δx t /|δx t | is a time-dependent unit vector aligned with the separation vector between the two particles. Its dynamics follows from Eq. (24): The expressions (24) - (27) can be expanded analogously to the one-dimensional case described in the previous section. For the model flow given by Eqs. (4) and (5) we find to order Ku 8 : (1 + β 2 ) 2 + Ku 6 423 + 972β 2 + 597β 4 + 20β 6 6(1 + β 2 ) 3 − Ku 8 164136 + 521517β 2 + 591081β 4 + 255803β 6 + 19927β 8 144(1 + β 2 ) 4 (28) and  Table 1 The lowest order, Ku 2 , is consistent with the results quoted in [18,19,20,21,8]: the maximal Lyapunov exponent changes sign at β c = 1. For β = 0, Eqs. (28) and (29) yield the Lyapunov exponents for particles advected in a two-dimensional incompressible Gaussian random velocity field with finite Kubo number: λ 1 + λ 2 = 0 (there cannot be clustering of particles advected in an incompressible flow), and Eq. (28) corresponds to the advective limit of Eq. (8) in Ref. [28] (the limit of zero Stokes number, St → 0, must be taken in this equation describing the the maximal Lyapunov exponent of inertial particles). The series (28) and (29) are expected to be asymptotically divergent. Just as the one-dimensional result (14), Eqs. (28) and (29) are expected to fail at large Kubo numbers.
The series (28) and (29) can be resummed as in the one-dimensional case. The results are seen in Fig. 2. Panel a shows results for the maximal Lyapunov exponent for β close to β c = 1 (the location of the path-coalescence transition in the limit of Ku → 0). We see that at finite Kubo numbers the path-coalescence transition occurs at β c (Ku) < 1. Comparing the first two terms in Eq. (28) shows that to order Ku 2 : For very small values of Ku this agrees with the numerical results in Fig. 2b. This panel shows the location of the path-coalescence transition in the Ku-β 2 plane. At large values of Ku, Padé-Borel resummations of the perturbation series (28) substantially improve the result. We observe good agreement between the numerical results and those of the resummation for values of Ku up to approximately 0.5. The results summarised in Fig. 2 show that at larger Kubo numbers less compressibility is needed to turn the maximal Lyapunov exponent negative. This is consistent with the behaviour observed at very large Kubo numbers: in the following section we show that the particles preferentially sample the attracting stagnation points of the velocity field in this limit. The contribution of these points increases as the Kubo number becomes larger. The observation that the effect of the compressible part of the velocity field is amplified at large Kubo numbers is consistent with numerical results in random renovating flows [27]. The dynamics of particles advected on the surface of a turbulent flow, by contrast, show a different behaviour [5]. In this case it is observed that the path-coalescence transition occurs at larger values of the compressibility for larger Kubo numbers. It would be of interest to determine which particular property of the turbulent flow gives rise to this effect.

Large-Ku limit
The large-Ku limit in two-dimensional compressible flows is solved as in one spatial dimension. The required distribution p(A) of the flow-gradient matrix A at u = 0 is found to be: where a = (A 11 , A 12 , A 21 , A 22 ) and The above expressions correspond to flows with both potential and solenoidal components because we expect that the expressions for the Lyapunov exponents derived below are not only valid for purely potential flows, but also yield good estimates for flows with a small solenoidal component. We change coordinates s ± = (A 11 ± A 22 )/2 and u ± = (A 12 ± A 21 )/2 to obtain The Lyapunov exponents are given by the eigenvalues of the strain matrix σ ± = s + ± s 2 − + u 2 + − u 2 − at the zeroes of u subject to the constraint Re σ ± < 0. This condition is equivalent to the condition TrA < 0 and det A > 0 and can be expressed as We have:

12
The factor 4 is a normalisation factor due to the fact that we only consider matrices A with Re σ ± < 0. Further Θ(z) takes the value Θ(z) = 1 when z > 0, and zero otherwise. To evaluate the integral we change variables according to with 0 ≤ w < ∞, −∞ < ξ < ∞ and 0 ≤ φ < 2π. Upon integrating over ϕ and ξ we obtain: Performing the remaining integrals we find an approximation for the Lyapunov exponents: The sum of the Lyapunov exponents is given by Eqs. (39) and (40) give estimates for the Lyapunov exponents for large but finite values of Ku, as opposed to Eqs. (28)-(29) that give the corresponding expressions for small values of Ku. Let us consider the purely compressible limit β → ∞ in (39) and (40): Both the maximal exponent λ 1 and the sum λ 1 +λ 2 are negative in this limit, as expected: the particles converge to the minima of V (x, t) = −ψ(x, t) at large Kubo numbers. Eqs. (41) and (42) (28) and (29) show. Fig. 3a shows the asymptotic result (41) for λ 1 in comparison with results of numerical simulations for particles suspended in a two-dimensional compressible (purely potential) velocity field. We observe good agreement.
We expect Eqs. (39) and (40) to give reliable estimates when the particles are typically found very close to the minima of the potential. Fig. 3b shows numerical results for the Lyapunov exponents at Ku = 100 in partially compressible flows as a function of β 2 , compared with Eqs. (39) and (40). We observe that Eqs. (39) and (40) yield reasonable estimates even for small degrees of compressibility. We note that the theory must fail in the incompressible limit (β = 0), and is only an approximation for finite values of β. The results show, however, that the dynamics is dominated by the attracting stagnation points of the velocity field at large Kubo numbers.

Conclusions
In this paper we have computed the Lyapunov exponents of small tracer particles advected in one-and two-dimensional compressible random velocity fields at finite Kubo numbers. For small Kubo numbers we have obtained results by Padé-Borel resummation of perturbation expansions in Ku. For large Kubo numbers we have computed the Lyapunov exponents using the Kac-Rice formula. At finite Kubo numbers the resulting Lyapunov exponents are determined by the details the velocity-field fluctuations (at small values of Ku by the Eulerian n-point functions of the velocity field and its derivatives, and at large values of Ku by the statistics of its stagnation points). Our results generalise earlier results [18,19,20,21,8] for compressible flows 14 obtained in the limit Ku → 0 to finite Kubo numbers. We find that λ ∼ Ku 2 as Ku → 0 and λ ∼ Ku at large (but finite) values of Ku. We have demonstrated that the analytical results are in good agreement with results of numerical simulations, and provide accurate estimates of the location of the path-coalescence transition for particles advected in compressible flows with finite Kubo numbers.
The limit β → 0 at large Kubo numbers remains to be analysed. For small values of β the stagnation points determining the Lyapunov exponents attract only weakly because the corresponding potential minima are shallow: Ku must be very large for the minima not to disappear before particles are attracted. In this limit a fraction of particles spends an appreciable amount of time on close-to closed orbits. This contribution is not accounted for in the derivation of Eqs. (39) and (40). For this reason these results are approximate, unless β is infinity. At β = 0 and Ku = ∞ the two-dimensional dynamics (3) corresponds to a one-dimensional Hamiltonian system. In this case the Lyapunov exponents must vanish. It would be interesting (but outside the scope of this paper) to consider, if possible, an expansion around this steady case.
We conclude by commenting on two further implications of our results. First, as mentioned in the introduction, spatial clustering of weakly inertial particles in incompressible velocity fields is often described in terms of a model where the particles are advected in a 'synthetic' velocity field with a small compressible component (due to the particle inertia). It is known that this approach must fail when the inertia is large (because of the formation of caustics). But there is also a problem in the small-inertia limit. Consider the sum of the Lyapunov exponents for inertial particles in a Gaussian random flow at finite Kubo numbers, Eq. (9) in Ref. [28]: Here St is the 'Stokes number' characterising the importance of particle inertia. The limit St = 0 corresponds to advective dynamics. Comparing Eq. (43) with the leading-order term of Eq. (29) at small St and β would lead us to conclude that weakly inertial particles are described by an advective model with 'effective compressibility' β = 3/2 KuSt. But this does not give the correct result for λ 1 (c.f. Eq. (8) in Ref. [28]), neither does this correspondence yield consistent results for terms of higher order in Ku (determined by higher-order correlation functions of the velocity field). This shows that care is required when approximating the dynamics of weakly inertial particles by advection in a weakly compressible velocity field: in general the statistics obtained by sampling along particle trajectories with actual inertial velocities, and in the effective compressible velocity field are different. Second, a related Ku-expansion was recently used to compute the tumbling rate of small axisymmetric particles in three-dimensional random velocity fields at finite Kubo numbers [36]. It turns out that the resummation works well also for the series expansion of the tumbling rate.