On the effect of finite-time correlations on the turbulent mixing in smooth chaotic compressible velocity fields

Most theoretical results about turbulent mixing have been obtained for ideal flows that are delta-correlated in time. As is often believed, those ideal flows are, with regard to mixing, very similar to real flows with a finite correlation time. However, recent results show that these two cases may differ considerably. In this article we study the effects of finite correlation time in a chaotic smooth statistically isotropic two-dimensional velocity field. As mixing is predominantly determined by the statistics of the stretching of material elements (e.g. lines “painted” onto a liquid), in this article we focus on the characteristics describing such stretching: finite-time Lyapunov exponents and the Lyapunov dimension. For these quantities, we derive analytical expressions as functions of the correlation time and the compressibility of the velocity field, and we investigate these expressions numerically. The results agree well with numerical results of other authors, and are useful for understanding several physical phenomena, e.g. patchiness of pollution spreading on an ocean and kinematic magnetic dynamos.


INTRODUCTION
Turbulent mixing of passive fields affects our daily life in various ways and plays an important role under a wide variety of physical conditions.Classical examples of scalar fields that can be considered in a good approximation to be passive and that are often subject to turbulent mixing include the temperature of fluids, the concentration of pollutants, nutrients, oxygen, etc. Therefore it is not surprising that turbulent advection has been studied extensively over several decades, theoretically, numerically, and experimentally.For more examples and details, cf.reviews [1][2][3][4][5][6][7].
An example of a passive vector field that is subject to a turbulent mixing is provided by the linear stage of the action of magnetic dynamos, for a review see [8].All the cosmic magnetic fields are due to the action of magnetic dynamos: when large masses of cosmic plasma are moving, the magnetic Reynolds number is very large and resistivity can be neglected: the magnetic flux is conserved through material loops.In three dimensions, due to the approximate incompressibility of the plasma and the exponential stretching of material lines, the cross-sectional area of magnetic flux tubes decreases exponentially, which gives rise to the amplification of the magnetic field.If there is no macroscopic field initially, molecular fluctuating fields will be amplified; at the initial linear stage, the magnetic pressure is negligible as compared with the thermal pressure, hence the magnetic field is passive.Note that for a working magnetic dynamo, in addition to the stretching of material lines, the so-called twist-fold action is also required (cf.[8]).In the context of the current paper it is worth noting that the theoretical results regarding the linear magnetic dynamos at sub-Kolmogorov scale rely upon the Lyapunov exponents of chaotic flows (cf.[9]).
The theoretical studies have been based almost exclusively on the idealized model of the Kraichnan flows [10], which are delta-correlated in time (i.e. with a zero correlation time).However, during the last decade, numerical evidence has been gathered showing that the match is far from perfect: in the case of compressible flows such as the flows at the free-slip surface of a turbulent liquid, the theoretical predictions based on the Kraichnan flows are incorrect, cf.[11,12].Recent theoretical studies [13] have not been able to explain the semi-numerical findings of [11,12].
In the case of chaotic flows with smooth streamlines, the properties of stretching statistics are most conveniently described in terms of finite-time Lyapunov exponents (FTLEs), which describe, roughly speaking, how the logarithms of the length, width, and height of a material parallelepiped will evolve.It is well known that the probability density function of FTLEs of smooth statistically isotropic Kraichnan flows evolves according to the convective-diffusive equation with constant coefficients; the values of these constants can be expressed in terms of the dimensionality and compressibility of the flow [1,14].It is the ratio of the coefficients of the convective and diffusive terms of this equation that defines the qualitative features of the mixing process; the coefficients themselves are just the largest FTLE and its variance.In the case of a sufficiently strong compressibility, the largest FTLE can become negative, which gives rise to a fractal clustering of dye particles [15][16][17].
In what follows, we derive analytical expressions for the FTLEs and the Lyapunov dimension as functions of compressibility and show that the results are in a good agreement with earlier simulation results.

THE MODEL AND METHOD
In the case of smooth chaotic velocity fields, the largest FTLE can be found as the average increment of the length of material vectors during the correlation time T of the velocity field in logarithmic scale divided by T (cf.[18][19][20]).This increment is defined by the strain tensor of a pair of material vectors, i.e. the time integral of the velocity gradient tensor, during the correlation time.Since what is important is only the time integral and not the detailed behaviour as a function of time, we shall assume for the sake of simplicity that the velocity gradient tensor is piece-wise constant: we assume that it is constant during the correlation time T , at the end of which it takes a new random and uncorrelated value.
In the case of compressible fields, the overall surface area of the basin occupied by the convective field is still constant, therefore we need to have both divergent and convergent regions present simultaneously.In order to make analytical progress, we adopt the simplest possible geometry and assume that the basin is divided into two regions, one with a positive velocity divergence d 1 and the other with a negative divergence d 2 = −d 1 .In order to keep velocities bounded, we add another constraint, namely that the determinants of the velocity gradients add up to zero.
Material vectors are stretched by the velocity gradients ∇v 1 and ∇v 2 .If a material vector x spends a time t in a velocity field v, then it transforms into the vector Dx, where the matrix D = e t∇v .As the real symmetric matrix D T D is positive semidefinite, its eigenvalues are real and non-negative.Thus, it is diagonalizable via a spatial rotation into diag(p 2 , q 2 ), where p 2 and q 2 are the eigenvalues of D T D. Then the squared length of the unit vector x = ( cos α sin α ) transforms into 2 = x T D T Dx = p 2 cos 2 α + q 2 sin 2 α.
The finite-time characteristic exponent (FTCE) for the length of a material vector spending time t in a constant velocity gradient, which corresponds to p i and q i , is given by the following average over directions: (The values of p i and q i also depend implicitly on t. The integral is given in [21, eq.2.6.38.4] and explained briefly in the appendix of this article.)The concept of FTCE for the length was introduced in [18].While asymptotically it becomes equal to the largest Lyapunov exponent, for finite times there is a small difference between the largest FTLE and the FTCE for the length (FTCE is a weighted average of all the FTLEs where the largest FTLE dominates).If the material vector spends time t 1 in environment 1 and time t 2 in environment 2, then its FTCE is the weighted average Let us define d i = ∇ • v i = Tr ∇v i and order the two environments so that d 1 ≥ 0 (region of expansion) and d 2 ≤ 0 (region of contraction).Particles from region 1 are eventually expelled into region 2; any particle that has entered region 2 stays there forever.Next, we consider the evolution of N 0 particles starting from region 1, which has an area A, and we wish to find the number N 1 (t) of them remaining in region 1 after time t.By Gauss's theorem, the flux of the particles out of region 1 is The average FTCE for material vectors that start from region 1 and end up in region 2, calculated for total time T and averaged over such material vectors whose number is N 12 , is where N 12 = N 0 − N 1 (T ).Adding the N 1 (T ) material vectors that spend all the time T in region 1 and the other N 0 vectors that spend all the time in region 2, we obtain the average FTCE for the length of a material vector: This λ 1 is also the FTLE for the longer semiaxis of a material ellipse.To obtain the FTLE for its shorter semiaxis, λ 2 , we first find the FTCE for its area, λ A .Formula (5) can be inverted to say that a material area in region 1 grows proportionally to e d 1 t , therefore the FTCE for a material area staying in region 1 is d 1 .For an area staying in region 2, the corresponding value is Analogously to Eq. ( 3), the FTCE for a material area that spends time t 1 in region 1 and time t 2 in region 2 is the weighted average The derivation of Eq. ( 7) readily generalizes to material areas.Substituting Λ 1 by d 1 , Λ 2 by −d 1 , Λ 12 by Λ A12 and taking the integral, we obtain the FTCE for areas, From this, we can calculate the FTLE for the shorter semiaxis of a material ellipse and the Lyapunov dimension: For every point in our numerical calculations, we generated ∇v 1 as a random matrix whose every entry is sampled from a uniform distribution between −5 and 5, and flipped its sign if its trace turned out to be negative (to ensure that region 1 is expelling).Then we calculated the compressibility ( If no such real ∇v 2 existed, we discarded the whole data point.Constraining ∇v 2 into being a sum of a diagonal and an antisymmetric matrix did not cause any loss of generality, because any real-valued tensor can be rotated into such a form (its symmetric part is diagonalizable by an orthogonal transformation), and rotation does not influence the FTLEs.
It is useful to express the correlation time through a dimensionless quantity called the Kubo number.For a given Kubo number K, the correlation time of our velocity field is because our construction implies that ∇v 1 = ∇v 2 .
(Equivalently, we could have normed the velocity gradients by ∇v 1 and equated the Kubo number and the correlation time.)This relation follows the definition used in article [13]: there K = u 0 T /η with u 0 being the typical speed of the flow and η being its correlation length.By using Eq. ( 2) and expressing p and q in terms of the matrix elements of the velocity gradient, it can be shown that for a general velocity gradient where ∆ = 4bc + (a − d) 2 is the discriminant of the eigenvalue equation for the velocity gradient.We used it in Eq. ( 7) to compute the larger FTLE λ 1 ; as the remaining integral in that equation could be evaluated analytically into a complicated expression involving hypergeometric functions, we preferred to calculate the integral numerically.Finally, we used Eqs ( 9)-( 11) to obtain λ A , λ 2 , and D L .

RESULTS AND SPECIAL CASES
The numerical results thus obtained are depicted in Fig. 1 together with the exact solution for deltacorrelated flows (for K = 0), D L = 2/(1 + 2℘).(This prediction of the Kraichnan model is given, e.g. in [11, eq. (3)].)Just as in the simulations in [11], when the correlation time grows, the Lyapunov dimension is reduced for small compressibilities and increased for large compressibilities.However, contrary to the cited article, we observe that D L is not completely determined by ℘.Instead, there is a range of possible values, with lower Lyapunov dimensions occurring more often in our statistical samples.If we averaged our results for a fixed compressibility, we would get a graph that would be very similar to the cited one, but such an averaging needs motivation.At the compressibility ℘ = 1, the point clouds condense into a single point whose location depends on the Kubo number.The reason is that conditions (13) constrain the velocity gradients into the form or rotations thereof.For such a velocity field, the discriminants ∆ 1 = ∆ 2 = 1, and Eq. ( 15) gives where the "+" of the "±" is for Λ 1 and the "−" is for Λ 2 .Substituting those into Eq.(3) then gives after substituting this into (6), we obtain ) (This integral can be quite simply proved by showing that its derivative with respect to K is zero.)This allows us to find, for the case ℘ = 1, analytical expressions for the FTLEs and the Lyapunov dimension, using Eqs ( 7) and ( 9)-( 11): where S = K + 2 ln cosh K 2 .For small Kubo numbers, D L ∼ 2 3 + 2 27 K; for large Kubo numbers, D L tends to unity as D L ∼ 1 − ln 2 2K−2−ln 2 .The fact that there is a range of possible D L for every ℘ suggests that there may be another parameter that, together with ℘, determines the Lyapunov dimension uniquely.The parameter must be dimensionless, and it must not change if the velocity fields are rotated.As it appears, (det ∇v 1 )/ ∇v 1 2 is a good choice.(To generalize this extra parameter to arbitrary velocity fields with non-constant gradients, an appropriate choice would be its root-mean-square.)Indeed, when it is added as the third dimension to the plots, it can be seen that the simulated points form a smooth surface.The shape of the surface is depicted in Fig. 2. As is evident from the contour plots, when the correlation time T becomes large (as does the Kubo number K), the contour lines tend to straighten and fan out from the point where ℘ = det ∇v 1 = 0.In the lower left corner of the plots there is a triangular region where the Lyapunov dimension practically vanishes, and there is a characteristic line above which it grows monotonically with the slope.All such behaviour can be characterized analytically.
Namely, after a long correlation time, essentially all the particles have entered the region of contraction (region 2) and stayed there for a long time.This means that all their characteristic dynamics is governed by the properties of the environment in the region of contraction, and the contributions of the region of expansion can be neglected.Therefore, for a large correlation time, λ 1 ∼ Λ 2 (T ).Looking again at Eq. ( 15), it can be seen that the qualitative behaviour of Λ 2 (T ) depends on the sign of the discriminant of the eigenvalue equation for ∇v 2 .Let the discriminant be ∆ 2 .For ∆ 2 > 0, the hyperbolic cosine grows exponentially, and √ ∆ 2 where d 2 is the divergence of v 2 , as previously.On the other hand, for ∆ 2 ≤ 0, only the first term remains: in this case, λ 1 ∼ d 2 .In both cases, the FTCE for areas, λ A ∼ d 2 .
The sign of ∆ 2 is uniquely determined by the parameter or equivalently by its counterpart δ 1 , calculated from v 1 , that by conditions (13) equals −δ 2 .The discriminant ∆ 2 is positive if and only if δ 2 < 1/4.The value of the Lyapunov dimension at the limit of long correlation times is also determined by δ 2 : Indeed, as det ∇v 2 = − det ∇v 1 , the isolines of D L are the (straight) lines where the ratio δ 1 -the slope of a contour line in Fig. 2b -is constant.The region of D L ∼ 0 is bordered by the slope δ 1 = −1/4.

CONCLUSION
We have derived analytically the finite-time Lyapunov exponents in smooth chaotic statistically isotropic velocity fields and used these results to express the Lyapunov dimension as a function of the compressibility of the velocity field.These results show that there is a clear departure between the ideal Kraichnan flows and the real time-correlated flows.Furthermore, unlike in the case of Kraichnan flows, the Lyapunov exponent is not defined purely by the compressibility, but it depends also on other statistical features of the velocity gradient tensor (the root-mean-square of its determinant can be used as a parameter).
Our results are in a good agreement with earlier numerical studies [11,12].We foresee use of our findings in the context of marine applications: first of all, for the studies of the pollution transport by sea currents on the free water surface, cf.[22].
As a further development, our approach will be extended to the calculation of the variance of the finitetime Lyapunov exponents and to the non-smooth velocity fields using the chaotic triplet-map model, cf.[23][24][25].

ACKNOWLEDGEMENT
The research was supported by the European Union Regional Development Fund (Centre of Excellence TK124: "Centre for Nonlinear Studies").
being the Frobenius norm), and chose ∇v 2 , if it existed, as one of the four unique matrices of the form a b −b d that give the same compressibility ℘, but whose determinant and trace fulfill the conditions det ∇v 2 = − det ∇v 1 , Tr ∇v 2 = − Tr ∇v 1 .

Fig. 1 .
Fig. 1.Lyapunov dimension as a function of compressibility.Grey dots depict the observed pairs of the Lyapunov dimension D L and compressibility ℘ for different Kubo numbers K. Points under D L = 1 are in the "strongly compressible" regime.Continuous lines plot the Kraichnan relationship D L = 2/(1 + 2℘).

Fig. 2 .
Fig. 2. Contour lines of the Lyapunov dimension D L .The darkest shading corresponds to D L ≈ 0.8 on (a) and D L ≈ 0 on (b); at the lightest shading, D L ≈ 2. Outside the shaded triangular area it is impossible to fulfill conditions (13).