{'fraction_non_alphanumeric': 0.04426578845408756, 'fraction_numerical': 0.013996512979465324, 'mean_word_length': 3.9398325358851674, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 26, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Universal long-time properties of Lagrangian statistics in the Batchelor regime and their application to the passive scalar problem

We consider transport of dynamically passive quantities in the Batchelor regime of smooth in space velocity field. For the case of arbitrary temporal correlations of the velocity we formulate the statistics of relevant characteristics of Lagrangian motion. This allows to generalize many results obtained previously for the delta-correlated in time strain, thus answering the question of universality of these results.


INTRODUCTION
The problem of passive scalar transport by turbulent flows has received much attention lately. The progress achieved has been made possible mainly by the introduction of the Kraichnan model [1]. Within the model the turbulent velocity statistics is believed to be Gaussian, scale invariant in space and δ correlated in time, which allows one to write closed equations on the correlation functions of the scalar. Such a velocity has only a few rough features in common with realistic flows, which are intermittent and have a finite scale-dependent correlation time, contrary to what is assumed in the model. Nevertheless, it seems that many interesting properties of the statistics are inherent in the dynamics, rather than due to the intermittency of the velocity statistics itself. Unusually for the turbulence theory numerous results have been obtained analytically using the Kraichnan model.
Having reached an understanding of this model, it is then natural to generalize its results, passing to more realistic flows. However, due to the complicated interplay between spatial and temporal properties of the velocity, one encounters various difficulties in introducing a meaningful velocity field with a finite correlation time. The only case where this was easily done is the so called Batchelor regime [2], where the spatial structure of the velocity is rather simple, and therefore one can separate space and time dependencies. It appears in the limit of large Prandtl numbers, which is the ratio of the fluid viscosity to the diffusivity of the transported quantity. In studying advection below the viscous length, the correlation functions of the velocity are smooth functions of space, which allows one to introduce an effective description with v α = σ α β (t)r β [2]. In this way time and space become completely separated.
Our aim here is to investigate the degree of universality of the passive scalar statistics for arbitrary temporal correlations of the velocity. We utilize the close rela-tion between the statistics of Lagrangian trajectories in a turbulent flow and the statistics of the passive scalar. Therefore, it seems reasonable to separate the problems, first investigating the Lagrangian motion and next applying the results to particular problems. For the Batchelor regime only a few degrees of freedom characterize the Lagrangian dynamics, which makes the problem solvable.
The plan of this paper is as follows. First we pass to the comoving reference frame in the equation for a passive scalar, which allows us to consider the Lagrangian mapping as an affine transformation, characterized by a random matrix. After its probability distribution function (PDF) is found, we consider several particular examples of the scalar statistics both for the decaying and forced turbulence. We show that the statistics of the scalar can be found by integration of the distribution function with a kernel, depending on the problem in question.

I. GENERAL RELATIONS
Advection of a passive scalar ϑ by incompressible velocity field v is described by the equation where κ is the molecular diffusivity. We shall be interested in the limit of small but finite κ. In the case of continuous injection of the scalar, one should add a source φ(t, r) into the right-hand side of Eq. (1.1). Let us consider a blob of the scalar having a size L much smaller than the viscous length of the velocity. The variation of the velocity on the scale of the blob is much smaller than the large homogeneous velocity transferring the blob as a whole. To account for a slow variation of the form of the blob due to relative motion of the particles, it is natural to pass to the reference frame moving with the velocity of a particle within the blob [14,15]. Since the velocity is a smooth function on the scale of the blob, it can be expanded in a Taylor series, thus leading to the equation Here σ αβ (t) is the matrix of the velocity derivatives taken at the chosen Lagrangian point. Incompressibility implies σ αα = 0. For turbulent flows σ should be regarded as a random matrix, having a finite correlation time τ , which is the Lagrangian correlation time of the velocity. The complete information about the Lagrangian flow, defined by ∂ t R = σR, is contained in the matrix W , satisfying This generates an affine transformation of space points, so that a vector R transforms as R(t) = W (t)R(0). The volume conservation guarantees det W = 1. The motion of the particles of the scalar differs from that of the space points due to the nonzero diffusivity. To investigate this motion we introduce the "inertia tensor" of the blob [14], where N = dr ϑ(t, r) is the number of particles of the scalar. It is easy to check that N is conserved by the full equation (1.2). It turns out that I contains all the necessary information, and will appear in the following sections as the result of formal calculations. The tensor I satisfies the closed dynamical equation The initial condition depends on the form of the initial blob, generally I αβ ∼ L 2 . We shall see that for problems with spatial isotropy it is enough to consider I αβ (0) = L 2 δ αβ . One can check that I can be expressed via W in a way that is nonlocal in time. Since I is symmetric and incorporates the diffusion, instead of working with W it will be more convenient for us to work directly with Eq. (1.4). The dynamics of the symmetric matrix I can be separated into the nontrivial essential dynamics of its eigenvalues and the trivial dynamics of the angular degrees of freedom. It is thus natural to reformulate the dynamics for the eigenvalues directly, excluding irrelevant angular degrees of freedom. For the case that is δ correlated in time, this can be done exactly, resulting in the Calogero-Sutherland model [5]. We shall show that for a finite correlation time of σ the angular degrees of freedom can also be effectively excluded. The reason for this is that only the large-time dynamics of the eigenvalues is important for our purposes so that in many respects (but not all) the matrix σ appears to be δ correlated in time.
Before we proceed with the derivation, it is useful to understand qualitatively the typical dynamics of a blob. If the amplitude of the velocity fluctuations is large enough (the precise condition will be formulated below), the term κ on the right-hand side of Eq. (1.4) can be disregarded during the initial stage of the evolution. Then one can make sure that I coincides with L 2 W W T . According to the Oseledets theorem [16], at large enough times the logarithms of the eigenvalues of the latter matrix are asymptotically equal to 2λ i t. The Lyapunov exponents λ 1 , . . ., λ d do not depend on a particular realization of σ; hence they are important characteristics of the system.
The above implies that the directions corresponding to positive and negative Lyapunov exponents will grow or decrease correspondingly, and that the blob will become an ellipsoid with the length of its main axis changing as exp(λ i t). The orientation of the ellipsoid can be arbitrary. The smallest dimension will decrease exponentially with the rate λ d , until at t ≈ |λ d | −1 ln L 2 |λ d |/κ it reaches a scale r dif = κ/|λ d |, where the diffusive spreading of particles makes further contraction impossible. Later, the smallest direction will fluctuate around r dif . This will not affect other directions, that will continue to change according to their Lyapunov exponents. In order to have a wide separation of the scales L and r dif , one should require a large value of the Peclet number This ensures that the time needed to reach the diffusion scale is large, so that the above arguments are valid. Apart from the typical event described here, we shall also need the distribution of all outcomes. This is the aim of Sec. II. Although it is not difficult to work with the arbitrary dimensionality of space, we shall consider the physical dimensionalities d = 2 and 3 only.

II. STATISTICS OF I
To separate the angular degrees of freedom from the radial ones, it is natural to represent I as follows: Here R is an orthogonal matrix composed of the eigenvectors of I, and Λ is a diagonal matrix with the eigenvalues e 2ρ1 , . . . e 2ρ d along the diagonal (we believe that the eigenvalues are ordered, so that We do not assume a summation over the repeating indices in Eqs. (2.2) and (2.3). This system of equations is not very useful for analyzing the general case. However, one can notice that if during the evolution the eigenvalues become widely separated, that is ρ 1 ≫ · · · ≫ ρ d , the system is greatly simplified. In this case the antisymmetric matrix Ω becomes and due to Eqs. (2.2,2.3) the dynamics of the angular degrees of freedom is independent on the eigenvalues. Therefore, Eq. (2.2) can be resolved: Here ρ 0i are some constants of the order of unity that should be determined by matching with the initial period of separation of the eigenvalues. Integrals in Eq. (2.5) determine the dynamics of ρ i . We consider times much larger thanτ , the correlation time ofσ, which is generally less than or of the order of τ . The form of the probability distribution function of ρ is different if we considerτ ≪ t < ∼ |λ d | −1 ln Pe or t > ∼ |λ d | −1 ln Pe, depending on whether the diffusion has started to be relevant or not. If the former case is considered, one can disregard the second term in Eq. (2.5), thus obtaining ρ i ≈ t 0 dt ′σ ii (t ′ ). We recognize the case of the central limit theorem. However, the Gaussian distribution describes only the bulk of the most probable events leaving rare events out of the domain of its validity. We shall need a more general expression [17] which can be derived from the following considerations. The integrals can be considered as sums of a large number n ≈ t/τ of independent identically distributed random variables. Thus we investigate the distribution of X given by Without loss of the generality we can assume that x i = 0. If the generating function exp[iyx i ] of each x is exp[−s(y)], then X has the generating function exp[−ns(y)]. To find the distribution function of X, one should make the inverse Fourier transform At large n this integral can be calculated in the saddlepoint approximation. Writing the extremum condition, we see that y extr is a function of the argument X/n, which implies P(X) ∝ exp[−nS(X/n)]. For X ≪ n one can expand S in the Taylor series and obtain P ∝ Here ∆ is the variation of x i . This is nothing but the central limit theorem. On the other hand, if we increase n, keeping the ratio X/n a constant of the order of unity, we can assert that ln P ∝ −n. This has a simple interpretation. If X is of the order of n, only realizations where most of x i are of the same sign contribute. Therefore we can model the situation by the binomial random process, which gives just the above result.
If we replace sum (2.6) by the integral t 0 dt ′ x(t ′ ) of a random function x over time t much larger than the correlation time, we should only note that the characteristic function of X is proportional to exp[−ts(y)], and then proceed as above. We used the fact that the characteristic function is an exponent of the cumulant generating function. The derivation is easily generalized for several quantities. Thus the distribution functions in d = 2 and 3 are given by the formulas Here S 2 (x 1 ) and S 3 (x 1 , x 2 ) are some functions depending on the details of the statistics of σ. In the δ correlated case one can find the explicit expression for S 2,3 [6,13] (see also Appendix A). The constants λ i are nothing but the Lyapunov exponents, which are expressed via the statistics of σ in the following way (cf. Ref. [16]) To have a self-consistent picture, we should assume that the spectrum of the Lyapunov exponents is nondegenerate, that is λ i > λ i+1 . Physically nondegeneracy of the Lyapunov exponents means that a blob is unstable with respect to the fluctuations, leading to a separation of the lengths of its sides. Noticing that the Lagrangian point r = 0 is a saddle point for the incompressible flow, one can easily verify that the strain directions corresponding to the further elongation of the blob prevail. Therefore, during a time τ of approximately constant strain, the blob will be on average further elongated. At t ≫ (λ i − λ i+1 ) −1 we can disregard effects originating from the boundary ρ i = ρ i+1 . Equations (2.7) and (2.8) are not valid in a narrow region near the boundary which has a width of the order of unity. Since it is much smaller than λ i t, we can use the step function θ to model the form of the PDF near the boundary.
Due to the incompressibility condition the exponents satisfy d i=1 λ i = 0. Then, in order to have a spectrum that is nondegenerate in d = 2, one should only require that λ 1 > 0, which is the same saying that trajectories diverge exponentially. In d = 3, it is necessary to supply some information about λ 2 . If the statistics of σ is symmetric with respect to time reversion, then λ 2 = 0 [13]. However, if this is not the case, it is generally nonzero. In Appendix B we find the expression for λ 2 if the correlation time of σ is small, which shows that its sign is generally arbitrary.
The form of the functions S 2,3 depends on particular details of the statistics of σ. However, it is possible to make two general statements about these functions. First, one can assert that at small x the expansion is valid, reproducing the central limit theorem. The constants C ij are defined as Here . . . stands for irreducible correlation function. The integrals should be calculated over an interval, much larger then the correlation time ofσ. Note that the condition of incompressibility ensures that d j=1 C ij = 0. When x i are of the order of unity, the functions S 2,3 have no singularities and change smoothly. The quadratic expansion of S 2,3 is valid as long as whereτ is the correlation time ofσ. In the δ-correlated case it holds everywhere (Appendix A). The normalization of P is determined by the quadratic part of S 2,3 , since most of the probability is concen- Now consider t > ∼ |λ d | −1 ln Pe. The diffusion is irrelevant for ρ i having a non-negative Lyapunov exponent. However, there is a finite probability, increasing with t, that ρ d reaches the diffusion scale. This requires an account of the last term on the right-hand side of Eq. (2.2) or Eq. (2.5). The diffusion will not allow ρ d to decrease much below ln(|λ d |/κ). On the contrary, negative λ d will prevent it from increasing. As a result, the corresponding ρ i will be distributed stationarily around the value ln(κ/|λ d |). Relaxation times associated with this distribution are diffusion independent and thus are much less then t. On the other hand, ρ having non-negative Lyapunov exponents are the integrals over the whole evolution time t, so that their values at time t are not sensitive to the last period of evolution with duration of the order of the relaxation time of ρ d . This means that fixing their values at time t ≫ |λ d | −1 ln Pe will not affect the distribution of ρ d , and the whole probability distribution function P is factorized (cf. Refs. [3,8]). In d = 2 we can write Here P st is the stationary distribution of ρ 2 . In d = 3 the situation is more complicated. While λ 3 is always negative, λ 2 can be both positive and negative. The form of the PDF will be different for these two cases. If λ 2 ≥ 0, Since ρ d is independent of the rest of ρ i the functions S 2,3 are the same as in Eqs. (2.7) and (2.8).
If λ 2 is negative, then at t ≫ |λ 2 | −1 ln(|λ 2 |/κ) the distribution over ρ 2 will also become steady and concentrated near ln(|λ 2 |/κ) which by order of magnitude is equal to ln(|λ 3 |/κ). Therefore, our assumption that ρ 2 ≫ ρ 3 is incorrect. Still ρ 1 ≫ ρ 2,3 , and the equation for ρ 1 is separated from the other variables. Then the distribution function is equal to Finally, let us note that since the configuration space SO(d) of the rotation matrix R [see Eq. (2.1] is finite, and since there is no preferred direction in space, at large t the matrix is distributed uniformly over the sphere.
The basic result obtained above is the special scaling form of the probability density functions. It is this universal form which lies in the origin of the results derived below.

III. DECAYING TURBULENCE
As a first application, let us consider decay of a passive scalar ϑ. The problem is posed as follows: given a random distribution of the scalar density ϑ 0 at t = 0, find its statistics at t > 0. In the framework of the Kraichnan model the single-point statistics was considered by Son [6], who obtained the following long-time asymptotic behavior: where γ α in d = 3 is equal to α(6 − α)D/4 for 0 ≤ α < 3 and 9D/4 otherwise (D is a parameter characterizing the strength of the fluctuations of σ). The same decay law has been claimed for the gradients of the scalar. Here we consider the problem for arbitrary correlation time of σ. Our consideration shows that due to the abovementioned special form of Eqs. (2.11-2.13) the law (3.1) is valid for an arbitrary statistics of σ both for the singlepoint value of the scalar and its gradient. In the δcorrelated limit we obtain a result for γ α different from that of Ref. [6]. The results also show that the basic assumption of Ref. [18] is incorrect. The following qualitative picture, supported by the calculations presented below, explains the decay. First, consider a single blob initially having a characteristic size L and containing N particles of the scalar. As velocity stretches the blob, the number of particles does not change, contrary to the volume of the blob. At t > ∼ |λ d | −1 ln(|λ d |L 2 /κ) the dimensions of the blob with negative Lyapunov exponents are frozen at r dif , while the rest keep growing exponentially, resulting in an exponential growth of the total volume of the blob. Since the volume is proportional to √ det I = exp ( ρ i ), one has |ϑ| α ∝ exp (−α ρ i ) , where the averaging should be done with the help of the PDF discussed above, that is |ϑ| α ∝ dρ P(ρ) exp (−α ρ i ). The result is determined by a compromise between two competing factors. While the averaged quantity exp (−α ρ i ) favors smaller values of ρ i , the maximum of the probability is attained when each growing ρ i is equal to λ i t. Obviously, for larger α the volume acquires more importance, so that the main contribution is made by smaller blobs, which are less probable but have larger concentration of the scalar. So, at small α, the deviation from the average growth λt is small, and γ α is determined by the Gaussian part of the PDF. This gives a parabolic dependence on α. On the other hand, if α is large enough, the main contribution is due to a blob having a minimal possible volume which is of the order of L d . The decay exponent is fully determined by the probability to have such a blob and hence is α independent [3,6]. Note that this picture implies that the exponential decay holds at t > ∼ |λ d | −1 ln Pe.
If instead of a single blob one takes a spatially homogeneous problem, this consideration should be slightly modified. At large t, initially uncorrelated blobs are brought close to each other because of the contraction along a certain direction, and then they overlap diffusively. Since the number of overlapping blobs is large, due to the central limit theorem it is rather ϑ 2 which is inversely proportional to volume. Therefore, |ϑ| α ∝ exp (−α ρ i /2) .
Formally, one should solve Eq. (1.2) with the initial condition ϑ(0, r) = ϑ 0 (r). The solution is From the qualitative arguments it is clear that the longtime asymptotic should be independent of the particular form of the distribution. We will take the simplest statistics, which is Gaussian with the pair correlation function This particular form is chosen for further convenience.
In what follows we set L = 1. It is possible to generalize the calculation for arbitrary χ and show that the results are independent of its form.
To proceed, we introduce the generating function of ϑ: Here we assume averaging over the statistics of σ and the initial distribution of the scalar. The simplest part is to perform averaging over the Gaussian field ϑ 0 . To do this, we substitute Eq. (3.2) into Eq. (3.5) and using the expression for the characteristic function of a Gaussian random variable [17], obtain Substituting χ(k) = (8π) d/2 χ 0 exp(−2k 2 ) and integrating over k we obtain We used the fact that I = W W T + Q, which can be verified by writing equation on W W T + Q and comparing it with Eq. (1.4). Using Eq. (3.6), one can find Here C α is a numerical constant. Equation (3.7) reduces the problem to an averaging of powers of det I = exp(2ρ i ), where P is the probability density function of ρ discussed above. In the large time limit this integral can be calculated in the saddle-point approximation. The calculation is slightly different for d = 2 and 3.

A. Two-dimensional case
In d = 2 integral (3.8) should be calculated only over ρ 1 , since the distribution over ρ 2 is stationary. The saddle-point equation is It is clear that ρ 1 ∝ t. As long as one can use the quadratic expansion, which valid at least at small α, the solution of this equation is ρ 1 = (λ 1 − αC 11 /2)t, hence At α > α cr = −2S ′ 2 (−λ 1 ) the value of ρ 1 becomes much smaller than λ 1 t; the integral (3.8) is determined by the boundary of the integration region, and therefore γ α = S 2 (−λ 1 ), independent of α.
The domain of validity of Eq. (3.9) depends on the value of the parameter λτ . If it is much smaller than unity, we can use the quadratic approximation to S 2 everywhere. This case effectively corresponds to the Kraichnan limit [6].
In the opposite limit λτ > ∼ 1, the quadratic expansion of S 2 cannot be used for α > ∼ 1/(τ C 11 ), and Eq. (3.9) is valid only at α ≪ 1/(τ C 11 ). The form of the intermediate region is not universal, and depends on the particular form of S 2 and hence on details of the statistics of σ.

B. Three-dimensional case
In d = 3 the result is similar to that of d = 2, though the consideration is slightly more complicated, due to the presence of an additional degree of freedom. There are two cases to be considered. If λ 2 < 0, then P is given by Eq. (2.13) and the calculation is the same as for d = 2.
If λ 2 ≥ 0, both degrees of freedom ρ 1 and ρ 2 are active and one should use the PDF (2.12). The saddle-point equations are where x 1 = (ρ 1 − λ 1 t)/t and x 2 = (ρ 2 − λ 2 t)/t. Again, the beginning of the curve is determined by the Gaussian part of P, and γ α is parabolic: with ρ 1 = (λ 1 + αC 13 /2)t and ρ 2 = (λ 2 + αC 23 /2)t. At α larger than α cr calculated below, we have the αindependent behavior Depending on the parameters, two different types of behavior can occur at α < α cr . First, it is possible that as α increases, ρ 2 will grow more slowly with t and at certain α will become much smaller than λ 2 t. At larger α the integration over ρ 2 will be determined by the region ρ 2 ≪ λ 2 t and one should replace system (3.10) by the single equation The other possibility is that at certain α the difference ρ 1 − ρ 2 becomes much smaller than ρ 1,2 . Because of the constraint ρ 1 > ρ 2 , for larger α integral (3.8) is determined by the boundary ρ 1 = ρ 2 of the domain of integration and the saddle-point equation becomes (3.14) Then α cr = − [∂ 1 S(x 1 , x 2 ) + ∂ 2 S(x 1 , x 2 )] at x 1 = −λ 1 and x 2 = −λ 2 . Geometrically, the first case corresponds to elongated ellipsoids, and the second one to sheets having two largest dimensions of the same order. If the parameter λτ is small enough, these changes of the regime occur within the Gaussian part of S 3 . Then one can perform a more detailed investigation. The first regime is realized if λ 1 > λ 2 C 13 /C 23 and C 23 < 0. Then, at α > −2λ 2 /C 23 Eq. (3.11) should be replaced by If λ 2 = 0, Eq. (3.11) has no region of validity, and at α > 0 one should use Eq. (3.15) which becomes In particular, within the Kraichnan model, Eq. (3.16) is valid for 0 ≤ α ≤ α cr . Substituting λ 1 and C ij (see Appendix A) one finds α cr = 4, for α < α cr and γ α = 3D for α > α cr . Our result is different from the one obtained in Ref. [6], which coincides rather with Eq. (3.11). An exact solution for α = 2 (see Appendix C and Ref. [9]) supports Eq. (3.17). The reason for the discrepancy is the following. Despite the fact that ρ 2 ∝ ln(κ/D), it is impossible to ignore it completely. If this were done, the anticorrelation between ρ 1 and ρ 2 , existing due to the incompressibility condition would lead to the growth of ρ 2 , thus making the calculation inconsistent.
The second regime takes place if C 23 > C 13 and C 23 > C 13 λ 2 /λ 1 . Then, starting from α = 2(λ 1 − λ 2 )/(C 23 − C 13 ) Eq. (3.11) should be replaced by another formula. Although the dependence on α is still parabolic, the coefficients are rather cumbersome, so we do not write this here.

C. Gradients of the decaying scalar
In the same manner one can consider the decay of the gradients of the scalar. In analogy, we can look for the correlation functions |ω| α , where ω = ∇ϑ(t, r)| r=0 . As in the case of single-point scalar statistics, these correlation functions decay exponentially in time. It was claimed in Ref. [6] that the decay law of the scalar and its gradient is the same within the Kraichnan model. Here we show that this is actually the case for arbitrary correlated strain. Qualitatively it follows from the estimate that |∇ϑ| ≈ ϑ/l, where l = exp(ρ d ) is the smallest dimension of the blob. As explained above, ϑ and l can be considered as independent, while the statistics of l is stationary. Thus the decay of the gradient is solely due to the change of the density of the scalar. More formally, one has Introducing the function Z(y) = exp [i (y, ω)] , and averaging it over the initial distribution (3.4) we obtain a formula similar to Eq. (3.6) Then, making a Fourier transform over y, we obtain the PDF of ω: Considering this expression in the eigenbasis of the matrix I, we observe that |ω| α ∼ |ω d | α , since ρ d is smaller than the rest of the ρ is. Recalling that the distribution over ρ d is stationary, we immediately obtain that which, due to Eq. (3.7), gives the same law of decay.

A. Single-point distribution of ϑ
In this section we shall investigate the steady state distribution of a passive scalar which occurs in the presence of a stationary source. For this purpose we introduce a random function φ(t, r) on the right-hand side of Eq. (1.2), injecting blobs of the scalar with the characteristic size L. Due to the linearity of the problem, the scalar field at the moment t = 0 is given by a superposition of the scalar injected at earlier instants of time. Each realization of σ can be characterized by a parameter t * (cf. Ref. [3]), such that the smallest dimension of blobs injected at t ≈ −t * approaches r dif at t = 0. The ambiguity in the definition of t * is of the order of |λ d | −1 which is much smaller than the typical stretching time |λ d | −1 ln Pe. Considering the motion of the scalar injected at −t * < ∼ t < 0 one may neglect diffusion, and the scalar is simply transfered along Lagrangian trajectories. On the other hand, as discussed in Sec. III, the contribution of the scalar injected at t < ∼ −t * is exponentially small. Thus t * separates diffusive and diffusionless regimes. One can write the following approximate formula: If the correlation time of the source is much smaller than t * , for a fixed realization of σ integral (4.1) can be considered as a Gaussian variable with zero average and the dispersion proportional to t * , so that after averaging over φ, for the single-point PDF [3] one obtains where χ 0 = dt φ(t, 0)φ(0, 0) . The effective Gaussianity of the pumping has its limitations due to finite correlation time of φ [analogously to the discussion of Eqs. (2.7) and (2.8)]. Since we work in the comoving reference frame, this correlation time is very small, and hence only the tail will be affected. At the end of the section we discuss the implications of this. To proceed with formal calculations, we should specify statistics of φ. Here we shall take φ as a Gaussian field with the pair-correlation function where χ(r) is the same as in sec. III. The expression for the generating function Z = exp(iyϑ) follows Here I(t, t ′ ) is a matrix, satisfying Eq. (1.4) with respect to t and the initial condition I(t ′ , t ′ ) = 1. Naturally, the integration is performed over the initial time, summing up the blobs injected at different times. The integral in Eq. (4.3) gives the formal definition of t * entering Eq. (4.2): Introducing the distribution function p(t * ), we rewrite Eq. (4.2) as Since t * is a functional of the whole trajectory ρ i (t, t ′ ), one needs more information than contained in the simultaneous distribution function of ρ. However, the following approximation, becoming exact at ln Pe → ∞, reduces the problem to single-time statistics. We shall neglect the configurations for which the smallest dimension of the blob once reached r dif starts to grow. Then the realizations for which t * is larger than some T , and those for which the blob injected at −T has ρ 3 (0, −T ) > r dif , are the same, leading us to the following formulas: Here one should substitute PDF's (2.7) and (2.8), since t * is determined by the diffusionless regime. These equations define nothing but the flux of the probability out of the region ρ 3 > ln(κ/|λ d |), once returns are disregarded.
Investigation of the above integrals shows that p(t * ) has the following properties. Its main body is concentrated in the vicinity of t * = |λ d | −1 ln Pe, and has a width of the order of √ ln Pe. On the other hand, its tail t * ≫ |λ d | −1 ln Pe decays exponentially, where c is equal to S 2 (−λ 1 ) in d = 2 and S 3 (−λ 1 , −λ 3 ) in d = 3. The intermediate region is not universal and depends on the details of σ. Note that since Eq. (4.7) gives the probability that at large time t * a blob has not yet decayed, therefore c is equal to the limiting value of γ α (see sec. III). This information allows one to calculate the probability distribution function P(ϑ). If ϑ ≪ ln(|λ d /κ), it is the central peak of p(t * ) that determines the scalar PDF: As we increase ϑ, at ϑ > ∼ ln Pe the details of the distribution of t * become important. The Gaussian regime [Eq. (4.8)] will turn into some nonuniversal asymptotic. Nonetheless, at ϑ ≫ ln Pe, due to Eq. (4.7) the universality is restored: For a δ-correlated σ one can find the complete function P(ϑ). Then S 2,3 is Gaussian (Appendix A) and the result can be found in the saddle-point approximation. In d = 2 the result coincides with that of Refs. [4,10,11]. In d = 3 we obtain the formula (cf. Refs. [10,11]) ln P ∝ −3 ln 2 Pe + Dϑ 2 2χ 0 − ln Pe for |ϑ| < 4 χ 0 /D ln Pe, and ln P ∝ 6 ln Pe − 4 3 ln 2 Pe + 3Dϑ 2 8χ 0 otherwise. The change of the regime is related to the fact that the two dimensions of the contributing blobs start to be equal to r dif . This result is different from the one presented in Ref. [10]. The difference can be qualitatively explained as follows. In our case the structures of the scalar making the main contribution to the PDF are columns, with the two smallest dimensions of the same order. They appear because of the anticorrelation originating from the incompressibility condition: for t * larger than a mean value, ρ 3 should decrease slower than λ 3 t, which by virtue of the anticorrelation leads to a decrease of ρ 2 faster than λ 2 t. Staring from a certain value of ϑ, both ρ 2 and ρ 3 decrease at the same rate (an analogous phenomenon is described in subsection III B). This structure is different from the ansatz proposed in Ref. [10]. Equation (4.2) should be modified if φ has a finite correlation time τ φ [see the discussion leading to Eqs. (2.7] and (2.8) ). That is where f (x) deviates from x 2 /(2χ 0 ) at x > ∼ 1/τ φ . This may affect only the tail of P(ϑ). If the parameter τ φ √ χ 0 c is much smaller than unity the tail is determined completely by the region where f ∝ x 2 , and one obtains the asymptotic result (4.9). Conversely, at τ φ √ χ 0 c > ∼ 1 the form of f should be accounted. Nevertheless, one can easily check that the exponential tail survives with a decrement depending on the form of f .

B. Gradients
Here we briefly consider the statistics of the scalar gradients. Within the Kraichnan model the problem was solved in Ref. [7]. From the qualitative picture presented there, one can conclude that the PDF is determined by the short-time fluctuations of σ, and hence is nonuniversal. The following considerations support the conclusion.
In a way, similar to the one leading to Eqs. Analogously to the case of the scalar density, the gradient field at t = 0 is given by a superposition of contributions of blobs injected at earlier moments of time. We observe that the contribution of each blob into ω 2 ≡ (∇ϑ) 2 is determined by two factors: the value of the scalar density χ 0 (det I) −1/2 and the inverse size of the blob contained in I −1 αβ . Not all the blobs make a contribution to Eq. (4.10). Indeed, the size of the blobs injected at |t ′ | ≪ |t * |, (where t * is defined in subsection IV A) is much larger than r dif and, therefore the value of the gradient will be small. On the other hand, the scalar injected at |t ′ | ≫ |t * | has an exponentially small density and hence does not contribute. Thus the distribution of ω is determined by the blobs injected at t ′ ≈ −t * which have the minimum possible size provided the diffusion is still ineffective at t = 0. Each realization of σ can thus be roughly characterized by two relevant parameters. The first one is the lateral dimension l of the thinnest blobs, for which the diffusion can still be neglected at t < 0. This is related to the very last stage of the evolution, when blobs of the smallest size of the order of r dif may undergo a strong rapid contraction, increasing the gradient without dissipating the scalar. Let us stress that the fluctuation should be short lived in order to suppress the diffusive spreading of the particles.
The other parameter is the duration T 0 of the injection stage for these blobs, showing how many blobs approach r dif at t ≈ 0. There is no average strain during this period, so that blobs injected at −t * − T 0 < ∼ t ′ < ∼ −t * all have a size of the order L at t ≈ −t * . Since at t > −t * the blobs move in the same velocity field, they all have approximately the same size at t = 0. Formally, the number of relevant blobs is expressed by the formula [c.f. Eq.
On the average T 0 ∼ |λ d | −1 and l ∼ r dif , so that . Nevertheless, studying the tail of the gradients PDF, it is necessary to take into account the large deviations of these parameters. The probability of a large value of T 0 is related to configurations of small strain (see Secs. III and IV A), and decays as exp(−cT 0 ). Writing P(ω) ∼ exp[−ω 2 l 2 /(2χ 0 T 0 )] T0,l , one can average over T 0 and find To average over l one notes that the tail l ≪ r dif of the probability distribution function of l is related to the tail of P st (ρ d ) [see Eqs. (2.11-2.13)] via l = exp(ρ d ). Indeed, both are determined by the probability of a strong and rapid contraction from r dif to l ≪ r dif . Hence Within the Kraichnan model ln P st ∝ −l −2 (see Appendix A) and the result ln P(ω) ∝ −|ω| 2/3 of Ref. [7] easily follows. In general the fluctuations of the smallest dimension take place at times of order τ near t = 0, and therefore are related to the single-time distribution ofσ, which is nonuniversal. We conclude that the gradient statistics is nonuniversal, and cannot be predicted unless some specific information is supplied. For example, if the distribution of l falls off very fast at small l, the distribution of ω will have an exponential tail at very large ω. In particular, this can explain the results of numerical simulations [3,19], where a cutoff can be due either to the grid step or imposed by hand [3]. If the tail of P st behaves according to ln[P st (l)] ∝ −l −α , the tail of P(ω) has the stretched exponential form: ln[P(ω)] ∝ −|ω| α/(α+1) .

V. CONCLUSION
We considered a passive scalar advected by a random large-scale velocity field. Our purpose was to establish the degree of universality of the scalar statistics for arbitrary an correlated velocity. The investigation can be reduced to the statistics of different Lagrangian characteristics of the smooth flow. In the limit of large a Peclet number, part of the relevant information is contained in the long-time asymptotic properties of the Lagrangian statistics, which is shown to possess a universal form. The scalar quantities related to long-time evolution thus manifest universal statistical features. Generally, these are the central part and the tail of the corresponding PDF. We considered several particular examples of such quantities: the decay of the scalar density and its gradient, and the scalar density in the forced case. Conversely, the statistics of the gradients in the forced case requires information about short-time fluctuations of the velocity, and is thus sensitive to its details.
The application of the Lagrangian statistics established here is not only restricted to the above examples. One can slightly modify the procedure to consider manypoint correlation functions of the scalar, say the PDF of the scalar difference at two points [10], correlation functions out of the convective interval [12], and other problems. A similar scheme could be applied to other passive quantities, like vectors [20] and tensors [21].
The tensor structure is fixed by incompressibility condition. Zero correlation time allows to write the Fokker-Planck equation for the probability distribution of R and Λ. Integrating out the angular degrees of freedom, one can see that the equation obtained is equivalent to the following Langevin dynamics [5] where ξ i are random Gaussian delta-correlated processes with the following correlation functions Let us now consider a typical evolution of the eigenvalues. At t = 0 all the eigenvalues ρ i are equal to zero. Then, during short initial period of time, all ρ start to differ. We can always arrange ρ so that ρ 1 > ρ 2 > .. > ρ d . We observe that then the ballistic terms coth(2ρ i − 2ρ j ) are arranged in the same order so that the eigenvalues will continue to separate and at t ≫ D −1 the following inequalities will hold ρ 1 ≫ ρ 2 ≫ . . . ≫ ρ d . If this is the case, we can substitute the hyperbolic cotangents by ±1 and get the following equations This simplified dynamics can be easily turned into the probability density functions of ρ [13,6] P(t, ρ 1 , As explained in the main text these formulae are valid at times 1/D ≪ t ≪ 1/D ln(DL 2 /κ) and the form of the PDF near the boundary can be modeled by the step function. At times t ≫ 1/D ln(DL 2 /κ) one should also calculate the stationary PDF of ρ d . They are readily found from the one-dimensional Fokker-Planck equation In this section we assume that the correlation time τ of σ is small, namely Dτ ≪ 1, where D = tr dtσ T (0)σ(t) characterizes the amplitude of the fluctuations of σ. We investigate the effect of finite correlation time on the Lyapunov spectrum. In d = 2 there are no essential changes with respect to the δ-correlated case since both λ 1 and λ 2 = −λ 1 get small corrections in τ leaving all qualitative features unchanged. However in d = 3, one can ask whether λ 2 will shift from its zero value at τ = 0, and whether the correction is positive or negative at finite τ . We demonstrate that already the first order correction in τ leads to a generally non-zero value of λ 2 in d = 3, which can be both positive and negative.
Here it is more convenient to parameterize the angular degrees of freedom of W W T by the eigenvectors e i instead of matrix R (see Eqs. (2.2, 2.3)), given by e α i = R iα . If the eigenvalues are separated, that is ρ 1 ≫ ρ 2 ≫ ρ 3 , the equations for ρ 1 and corresponding to it eigenvector e 1 decouple ∂ t ρ 1 = (e 1 , σe 1 ) , ∂ t e 1 = σe 1 − e 1 (e 1 , σe 1 ) .
In the calculation it is convenient to deal with symmetric matrices, which is achieved by decomposing σ into symmetric and antisymmetric parts σ = s + ω and introducing e 1 = M n with ∂ t M = ωM , M (0) = 1 so that λ 1 = (n,sn) , ∂ t n =sn − n(n,sn) withs = M T sM . To find the first order correction to λ 1 we integrate the above differential equation from 0 to t, and then iterate the obtained expression once. After averaging over the directions of n(0) we find which is generally non-zero and has no definite sign.