Caustic Frequency in 2D Stochastic Flows Modeling Turbulence

Stochastic flows mimicking 2D turbulence in compressible media are considered. Particles driven by such flows can collide and we study the collision (caustic) frequency. Caustics occur when the Jacobian of a flow vanishes. First, a system of nonlinear stochastic differential equations involving the Jacobian is derived and reduced to a smaller number of unknowns. Then, for special cases of the stochastic forcing, upper and lower bounds are found for the mean number of caustics as a function of Stokes number. The bounds yield an exact asymptotic for small Stokes numbers. The efficiency of the bounds is verified numerically. As auxiliary results we give rigorous proofs of the well known expressions for the caustic frequency and Lyapunov exponent in the one-dimensional model. Our findings may also be used for estimating the mean time when a 2D Riemann type partial differential equation with a stochastic forcing loses uniqueness of solutions.


Introduction
We address a stochastic flow [1] covered by the following Ito stochastic differential equations dr = vdt, dv = −(v/τ)dt + dw(t, r), where phase variables r = (x, y) and v = (u, v) are interpreted as the position and velocity of fluid particles moving in a turbulent flow. Here the random forcing w(t, r) is assumed to be a Gaussian white noise in time and a statistically homogeneous random field in space, i.e., Ew(t 1 , r 1 )w(t 2 , r 2 ) T = min(t 1 , t 2 )B(r 1 − r 2 ), where B(r) is its space covariance matrix. A constant non-random parameter τ is called the Lagrangian correlation time.
In particular, the motion of a single particle in the framework (1) is covered by a Langevin equation for the velocity where w(t) is a standard Brownian motion in 2D, thereby the Lagrangian velocity of the particle is simply a well known (two-dimensional) Ornstein-Uhlenback process and its position is an integrated Ornstein-Uhlenback process.
In general the motion of any number n of the particles is described by a diffusion process in R 2n for the vector of positions/velocities (r 1 , v 1 , . . ., r n , v n ) with a drift and diffusivity matrix expressed in terms of τ and B [2].
To introduce the problem we focus on, let us supply (1) with initial conditions r| t=0 = a, v| t=0 = u 0 (a) where u 0 (a) is a given deterministic function (the initial Eulerian velocity field), i.e., each particle starts from a certain position with a certain velocity completely determined by its position, and introduce the Jacobian where r(t, a) is the solution of (1) satisfying the above initial conditions or, in simple words, the position of the particle at moment t starting (t = 0) at position a.
An important physical meaning of the Jacobian can be seen from the following. Let the initial position a be a random variable independent of the flow with probability density π a (x), then for the conditional density π r(t,a) (x) conditioned on the flow is given by π r(t,a) (x) = π a (x)/J(t, x) It is well known [3] that for incompressible flows J(t, a) ≡ 1 , hence the density of the particles does not change at any point in the fluid. A remarkable and one of the most important feature of model (1) is that it yields a compressible flow for any value of τ and any function u 0 (a), e.g., [4]. Moreover, J(t, a) may take zero values leading to an infinite density at certain points of the fluid called caustics.
Let ν be the mean frequency (intensity) of caustic occurrence and the top Lyapunov exponent (LE) for (1). Despite numerous studies of this important quantity in physical literature (see recent comprehensive reviews in [4,5]) almost all remarkable analytical results were obtained in the one dimensional case, among them explicit expressions [6,7] should be mentioned first where z = s −4/3 /4 , s is Stokes number (the ratio of Lagrangian and Eulerian time scales, e.g., [2]), M is the modulus of the Airy function. Furthermore, to our best knowledge no rigorous proofs of formulas (3) were presented so far. Furthermore, worth noting expressions for ν and λ in a quite different setup when the noise is a telegraph process rather then the Gaussian white noise [6]. Even though that results were proven rigorously it still was a one dimensional case that hard to consider as a turbulence model. As for two dimensions, we found no analytical results whatsoever concerning with the dependence of ν on s, while some asymptotics for LE were obtained in [2,6,7]. Moreover, exact expressions for LE were derived as well in very special degenerate cases [7].
Advances in this work could be summarized as follows.
(ii) In the framework of general 2D model (1) a system of stochastic differential equations was derived in four unknowns, involving the Jacobian.
(iii) The system was efficiently investigated for two cases of a special forcing in (1) yielding an exact asymptotic of ν(s) as s → 0 and reasonable numerical boundaries for ν(s) in a wide range of s. That results may be used for estimating the mean time when the corresponding SPDE for the Eulerian velocity field loses its uniqueness.
(iv) In one of that cases it was found that LE has exactly same expression as given in (3).
(v) For an isotropic turbulence the system was elaborated via introducing longitudinal and normal Stokes numbers.
In physical literature the model (1) most often was used to describe so called inertial particles such as water drops in atmospheric clouds, particulate matter or living organisms in the turbulent upper layer of oceans, and many other phenomena. Undoubtedly, real turbulent flows in the ocean and atmosphere are essentially 3D phenomena. However, the introduction of three dimensional perturbations does not destroy the main features of the cascade picture according to [8], implying that 2D turbulence phenomenology establishes a realistic picture of turbulent fluid flows.
In Eulerian terms Equation (1) cover Lagrangian motion in the Eulerian velocity field u(t, x) satisfying the following equation By the method of characteristics one can find that for some (maybe short) interval (0, t) the solution of this equation exists and unique since at the initial moment J(0, a) = 1. However, after some time the uniqueness is lost (the Jacobian vanishes) due to a very weak dissipation modeled by the last term on the left hand side (LHS). Thus, estimates of the mean of the first moment when J(t, a) hits zero would provide us with an idea when (4) loses uniqueness. In terms of the non-linear wave theory, one can treat that as the first moment of wave breaking.
Worth noting that if the dissipation term is replaced by a classical friction k∆u, then we arrive at the well known stochastic Burger's equation which has a unique solution for all t thereby no caustics may occur.
Some of our results can be predicted by simply proceeding to dimensionless variables in (4). Let U and L be some typical velocity and length scales, respectively. Changing u, w, x, t to u/U, w/U, x/L, t/τ one gets from (4) where the dimensionless parameter is called Stokes number. Thus, if s → 0 then the underlying Eulerian velocity field satisfies a linear equation and hence the intensity of caustics should tend to zero. In the opposite case s → ∞ the non-linear term dominates and one may expect that the number of caustics per time grows indefinitely.
The paper is organized as follows. A similar one-dimensional problem is exactly solved in Section 1. The solution is essentially used in investigating 2D phenomena. A deterministic case (w = 0) is briefly discussed in Section 2. In Section 3, we derive a system of equations containing J in the general case of an arbitrary homogeneous forcing. For particular forms of the forcing rigorous estimates for the mean number of caustics ν(s) are found in Section 4. Furthermore, ν(s) is computed by solving a simple 1D parabolic equation. The case of an isotropic forcing (unsolved yet) is mentioned in Section 5. Conclusions are gathered in Section 6 and some details are brought to Appendix A.

One-Dimensional Case
A part of results from this section is well known from physical literature [4,7,9], presented sometimes in a quite vague fashion. We formulate and prove them rigorously and give more details. Consider with Ew(t, x) = 0 and For J(t, a) = ∂x ∂a one can get by direct differentiation (6) in a Introduce the following dimensionless quantity p = τ J 1 /J and apply the Ito formula to the last system. As a result we have dp = 1 Then we proceed to a dimensionless time t/τ −→ t and obtain where w is a standard Wiener process and Worth noting that if the velocity and length scales of the flow are chosen as U = b(0)τ and L = −2b(0)/b (0), respectively, then (8) coincides with the earlier introduced Stokes number (5).
Assume that zeros of J(t) are prime that is typical for stationary processes, then they can be identified with moments of explosion of p(t).

Proposition 1.
Process p is explosive(e.g., [10]), more exactly where the explosion time S is finite P(S < ∞) = 1 and its expectation can be explicitly computed where Ai(z), Bi(z) are the Airy functions and M(z) their modulus.
Notice that (9) implies the expression (2) for the intensity of caustics since ν = 1/E(S). Explosiveness of p(t) follows from the Feller's criteria [10]. To prove (9) we introduce µ(x) as the mean time to explosion under condition p(0) = x, then solve Then we set E(S) = µ(∞).
and arrive at (9). Details of computations are the same as in [11]. For the purpose of studying 2D models the knowledge of the mean explosion time is not enough.
and let u(t, x) be the same probability under condition p(0) = x, then In [12], it was shown that u(t, x) satisfies the initial value problem To ensure uniqueness of solution of (10) we add the natural boundary conditions assuming that the limit u(t, ∞) exists similarly to m(∞).
A simplest Euler scheme was used for solving (10) numerically. Some details and validation of the choice of the scheme parameters are given in Appendix A. Notice that our goal is not to evaluate and minimize the error of the numerical computations, but rather to illustrate that (10) can be solved quite accurately and efficiently with very simple tools.  To proceed to investigating LE we notice that after explosion process p(t) can be continued by starting over, i.e., by solving same Equation (7) with the initial condition lim t→S+0 p(t) = ∞.
The positive sign follows from relation J 1 = dJ/dt and the fact that zeros of J(t) are prime. Thus, p(t) takes different signs on different sides of a zero of J(t). (7) is expressible in terms of the ergodic mean of p

LE for flow
as follows from (2) and definition of p(t), see [2] for details. It also can be exactly found Proposition 2. The ergodic mean of p is given bȳ The statement is proven in Appendix A and here we just make some comments.
Notice that the density of the stationary distribution π(p) of p(t) does exist and is given by −∞ e y 2 /2s 2 +y 3 /3s 2 dy (12) where C is a normalized constant. One can find the following asymptotic Thus, the integral for the invariant mean is formally divergent, however if the integral is meant as a Cauchy principal value then it coincides with the RHS of (11).

Deterministic 2D Case
In this section, we plug in (1) w = 0. The resulting linear equations are easy to solve and get the following expression for the Jacobian and (a, b), (u 0 , v 0 ) are the initial position and velocity, respectively, related by and S be the first time when J(S) = 0. Easy to find that In particular if Q = 0 and P < 0 do not depend on (a, b) at all, then G = R 2 and S , defined for all s > 0, does not depend on the initial point as well. Thus, the equation loses uniqueness exactly at moment t = S. If one defines u(t, x) after t > S as a solution of (13) with the initial condition u(S + 0, x) = u 0 (x) then ν = 1/S could be interpreted as the intensity of caustics in time. In Figure 2, we show a few curves ν(s) for different s 0 in order to compare them visually with a stochastic case addressed in the next sections. If S depends on the initial point, then its interpretation in terms of Eulerian setup (13) is not that simple and will not be discussed here.

System of Equations for Jacobian in General 2D Case
Now we return to the stochastic model (1). In a few works it was pointed out that there is a closed equation for the matrix Namely, e.g., [2,4] An analysis of that equation led to some general important conclusions, but it is of little help for our purposes because in general it cannot be reduced to efficiently handled scalar equations.
Let us first rewrite (1) in the coordinate-wise form with r = (x, y), a = (a, b) where It is not possible to obtain a closed equation for J(t), but it can be included in a system of four equation as it is shown in Appendix A. Namely, for dimensionless time t = t/τ denoted by the same letter, Jacobian can be represented as where dimensionless random function m 1 (t) is included into the following system s is Stokes number defined similarly to (8) and an exact expression for it is given in Appendix A.
To define processes η i (t) we nondimensionalize w (i) (t) and set where the subs mean derivatives. Thus, η's are dependent Wiener (non-standard) processes with a covariance matrix Edη (i) (t, x, y)dη (j) (t, x, y) = σ ij η dt given by where b ij is a dimensionless version of b ij and partial derivatives are taken at x = 0, y = 0.

Two Special Cases
Under certain conditions imposed on the forcing in (14) the first moment S of J(t) to hit zero turns out to be the minimum from the first explosion moments for two independent 1D processes described by (7).
Next due to the zero initial conditions m 3 = m 4 ≡ 0. Introduce By adding and subtracting first two equations in (15) we get two separated equations dp = (−p − p 2 )dt + sdw 1 , dq = (−q − q 2 )dt + sdw 2 (20) for independent identically distributed processes p(t) and q(t) Model 2. In this model we assume with independent U and V such that where the prime means the derivative in the space coordinate. From (16, 17, 21, 22) it follows that are independent identically distributed Wiener processes since η 1 = −2(U + V ) and Thus, m 2 = m 4 ≡ 0 and introducing arrive at the same equations (20). Below we show simulated velocity fields at a particular time moment for both models (Figure 3) where periodic in space forcings were used.
Let S 1 and S 2 be the first explosion moments for p(t) and q(t), respectively, then Certainly the knowledge of ES i is not enough to find E(S), but the latter can be expressed in terms ofF(t) = P(S 1 > t) as Manipulating with this formula it is not difficult to present an example where the expected value of the minimum of two independent identically distributed random variables X, Y is finite while the expectation of each variable is infinite due to heavy tails of X, Y distribution. So, theoretically speaking ES i and E(S) may differ by an order of magnitude. Fortunately it is not the case here since the tail of distribution of S 1 is exponential and can be evaluated [13]. Namely from Theorem 1.1 in [13] it follows that the exponential moment is finite for all γ = γ(s) satisfying and g(s) = max x : x −∞ exp(y 3 /3s 2 + y 2 /2s 2 )dy ∞ x exp(−y 3 /3s 2 − y 2 /2s 2 )dy .
Hence, it is reasonable to assume that with γ satisfying (24). In view of this assumption On the other side obviously that Thus, for the mean number of caustics ν(s) = 1/E(S) one gets from (9) For small s asymptotics of the lower bound and upper bound coincide and we get ν(s) ≈ 1 4πs 3/2 e − 1 6s 2 , while for large s the corresponding asymptotics differ just by a constant 0.4s 2/3 < ν(s) < 1.25s 2/3 . Approximation (25) is quite speculative and indeed it greatly overestimates ν(s). It can be seen by comparing the bounds with exact curve ν(s) obtained from solving the corresponding PDE forF(t) (Section 1). The upper and lower bounds for ν(s) are shown in Figure 4 as well as its numerical version.  In Model 2 LE cannot be found by such simple tools because the equations for x and y components do not split.
Finally, notice that LE in Model 1 is given by the same expression (11) as in 1D case. Indeed, in this case each component of the separation process ∆r = r(t, a) − r(t, 0) is the separation process for one dimensional flow (6). Hence |∆r| ≈ |a| exp (λt) for small |a| and large t, where λ is LE for (6). Then our claim follows from Definition (2).

Isotropic Forcing
Assume that in the original equations the forcing is isotropic, then its covariances are given by [3] where r = x 2 + y 2 . Assume smooth longitudinal and normal correlation functions Then from (17) it can be derived that (15) takes form here w j , j = 1, 2, 3, 4 are independent standard Wiener processes.
No essential progress is made in analyzing this system yet because none of variables can be eliminated. The only reason for presenting it here is a great importance of the isotropic case for applications.

Conclusions and Discussion
While one dimensional stochastic models for inertial particles have been comprehensively addressed, refs. [4,6,7,9] no essential progress in analytical studies for 2D turbulence have been reported yet. In this work for the first time the intensity of the caustics was rigorously treated in the framework of two dimensional stochastic flows modeling 2D turbulence in compressible media.
Namely, we revealed two cases, where the mean time to explosion in the particle density , µ(s), can be analytically estimated for the full range of Stokes number s and can be accurately evaluated by solving an initial/boundary problem for a one dimensional parabolic equation. Worth noting that both boundaries provide an exact asymptotic as s → 0. For both models the reciprocal ν(s) = 1/µ(s), interpreted as the intensity of caustics, increases monotonically from zero to infinity that is in full agreement with physics of phenomena in question.
The reported advance in the proposed models is due to reduction in the number of unknowns in the system (15) from four to two. Alas, for the most interesting isotropic case such a reduction is not possible, but a hope for finding asymptotics of ν(s) still remain. This is the main direction of our future studies.
Another interesting area of further consideration is the Equation (4) describing the Eulerian velocity field generating Lagrangian motion covered by model (1). It should be recognized that the interpretation of µ and ν in terms of solutions of (4) is not clear enough except the deterministic case with the initial velocity field u 0 (x) for which |∂u 0 /∂x| and ∇ · u 0 do not depend on x at all. A natural and important question is how the mean time until the uniqueness loss depends on the initial velocity field. The proof is based on the following statement Lemma A1. If there exist real m and smooth bounded function V(x) such that and then with probability onep = m Proof. From (A2) and Ito formula it follows that V(p(t)) satisfies dV = (p − m)dt + sV dw.
By integrating both sides we get The second term on the right hand side goes to zero as T → ∞ due to the large numbers law. The difference on the left hand side (LHS) of (A4) also converges to zero because of the boundness of V and condition (A1). Lemma is proven.
Notice importance of (A1). If V(∞) = V(−∞) then the jumps of V at explosion moments can accumulate leading to a non-zero limit of LHS. Now we directly construct such a function V(x) and determine m. First, we take a bounded solution of (A2) for which V(−∞) = 0. To ensure (A1) one should set The integral for M is meant as Cauchy principal value. That allows for changing the order of integration after substitution t = u − y, y = ys 2/3 Integration in y leads to To complete the proof one should account for [11], This relation also can be derived from the fact that the both sides satisfy the same differential equation and same initial conditions [14].