Density of Inertial Particles: Exactly Solvable 2D Models

Leonid I. Piterbarg

Department of Mathematics, University of Southern California

Kaprielian Hall, Room 108, 3620 Vermont Avenue, Los Angeles, CA 90089-2532

Abstract

Inertial particles in 2D driven by a Gaussian white noise forcing are
considered. For two examples of the forcing (compressible and incompressible)
upper and lower bounds are found for the mean number of caustics as a function
of Stokes number. Efficiency of the bounds is verified by numerical methods.

Keywords: Inertial particles, Caustics, Lagrangian stochastic models

Introduction

In compressible fluid flows the density of Lagrangian particles finite at an
initial moment can go to infinity in finite time as a flow is evolving. Such
density explosions are related to occurence of caustics in the underlying
velocity field. One of the most important such phenomena is motion of inertial
particles in turbulence ( a comprehensive review can be found in [1]).

In this work we address the mean frequency (intensity) of caustic occurence
$\nu$ in the framework of a stochastic flow modeling Lagrangian motion in a
turbulent velocity field. The stochastic flow concept allows one for
describing joint statistics of several particles or even a particle continuum
at any moment.

The main goal is to introduce two particular models in two dimensions, very
simple but not trivial, for which the dependence of $\nu$ on Stokes number $s$
(the ratio of Lagrangian and Eulerian time scales, [2]) can be analytically
investigated to a certain extent. Namely, first, lower and upper bounds for
$\nu(s)$ are provided yielding an exact asymptotic as $s\to 0$ and, second,
$\nu(s)$ is accurately evaluated for a wide range of $s$ by solving a $1D$
parabolic equation numerically.

A general version of a stochastic flow addressed here has been introduced in
[3] and assumes finite particle velocities and a white noise type of
accelerations. Thus, it is not a Kraicnan model of turbulence widely used in
physics of fluids [2]. More exactly, assume that the position
$\mathbf{r}(t,\mathbf{a})$ of a particle at time $t$ with the initial position
$\mathbf{a}$ and its velocity $\mathbf{v}(t,\mathbf{a})$ are described by the
following Ito equations

$d\mathbf{r}=\mathbf{v}dt,\ \
d\mathbf{v}=-(\mathbf{v}/\tau)dt+d\mathbf{w}(t,\mathbf{r}),\ \ \
\mathbf{r}(0,\mathbf{a})=\mathbf{a},\ \
\mathbf{v}(0,\mathbf{a})=\mathbf{u}_{0}(\mathbf{a})$ $None$

where $\mathbf{u}_{0}(\mathbf{x})$ is an initial (Eulerian) velocity field,
$\mathbf{w}(t,\mathbf{r}$) is a random forcing with zero mean which is a
Brownian motion in time and a statistically homogeneous random field in space
, i.e.

$E\mathbf{w}(t_{1},\mathbf{r}_{1})d\mathbf{w}(t_{2},\mathbf{r}_{2})^{T}=\min(t_{1},t_{2})\mathbf{B}(\mathbf{r}_{1}-\mathbf{r}_{2})$

where $\mathbf{B}(\mathbf{r})$ is its space covariance matrix, all the vectors
are two-dimensional, and finally $\tau$ is the Lagrangian correlation time
scale.

In particular, the motion of a single particle in the framework (1) is covered
by a Langevin equation for the velocity

$d\mathbf{v}=-(\mathbf{v}/\tau)dt+\mathbf{B}(\mathbf{0})^{1/2}d\mathbf{w}(t)$

where $\mathbf{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
$(\mathbf{r}_{1},\mathbf{v}_{1},...,\mathbf{r}_{n},\mathbf{v}_{n})$ with a
drift and diffusivity matrix expressed in terms of $\tau$ and $\mathbf{B}$,
[3].

In this work we mostly address the Jacobian

$J(t,\mathbf{a})=\displaystyle\det\left(\frac{\partial\mathbf{r}}{\partial\mathbf{a}}\right)$
$None$

and partially the top Lyapunov exponent (LE)

$\lambda=\displaystyle\lim_{t\to\infty}\frac{1}{t}\lim_{|\mathbf{a}|\to
0}\ln\left(\frac{|\mathbf{r}(t,\mathbf{a})-\mathbf{r}(t,\mathbf{0})|}{|\mathbf{a}|}\right)$
$None$

which depends on the direction of $\mathbf{a}$ for anisotropic flows.

A consideration of two-particle motion suffices for studying both
characteristics .

An important physical meaning of the Jacobian can be seen from the following.
Let the initial position $\mathbf{a}$ be a random variable independent of the
flow with pdf $\pi_{0}(\mathbf{a})$, then for the conditional pdf
$\pi(t,\mathbf{r})$ conditioned on the flow is given by

$\pi(t,\mathbf{r}(t,\mathbf{a}))=\pi_{0}(\mathbf{r}(t,\mathbf{a}))/J(t,\mathbf{a})$

Hence, if for some $t_{0}$ the Jacobian takes zero value, the density becomes
infinite at the point $\mathbf{r}(t_{0},\mathbf{a})$ what we will call an
occurrence of a caustic at this point.

We show that for a certain class of forcings caustics arise with probability
one, they form a stationary process in time and sensible bounds for the mean
number of caustics are given.

In Eulerian terms equations (1) cover Lagrangian motion in the Eulerian
velocity field $\mathbf{u}(t,\mathbf{x})$ satisfying the following equation

$\frac{\partial\mathbf{u}}{\partial
t}+\mathbf{u}\cdot\nabla\mathbf{u}+\mathbf{u}/\tau=\frac{\partial\mathbf{w}}{\partial
t},\ \ \ \ \mathbf{u}|_{t=0}=\mathbf{u}_{0}(\mathbf{x})$ $None$

By the method of characteristics one can find that for some short interval
$(0,t)$ the solution of this equation exists and unique since at the initial
moment $J(0,\mathbf{a})=1$. However after some time the uniqueness is lost
(the Jacobian takes a zero value) 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,\mathbf{a})$ hits zero would provide us with an idea
when (4) looses uniqueness. In terms of the non-linear wave theory, one can
treat that as the first moment of wave breaking.

Some of our results can be predicted by simply procceding to dimensionless
variables in (4). Let $U$ and $L$ be some typical velocity and length scales
respectively. Changing $\mathbf{u},\mathbf{w},\mathbf{x},t$ to
$\mathbf{u}/U,\mathbf{w}/U,\mathbf{x}/L,t/\tau$ one gets from (4)

$\frac{\partial\mathbf{u}}{\partial
t}+s\mathbf{u}\cdot\nabla\mathbf{u}+\mathbf{u}=\frac{\partial\mathbf{w}}{\partial
t}$

where the dimensionless parameter

$s=\displaystyle\frac{U\tau}{L}$ $None$

is called Stokes number.

Thus, if $s\to 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\to\infty$ 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 ($\mathbf{w}=\mathbf{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 $\nu(s)$ are found in Section 4. Also $\nu(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.

1\. One-dimensional case

A part of results from this section is well known from physical literature,
[1]. We formulate them in rigorous form and give more details.

Consider

$dx=udt,\ \ du=-(u/\tau)dt+dw(t,x),\ \ \ x(0)=a,\ \ u(0)=u_{0}(a)$ $None$

with $Ew(t,x)=0$ and

$Ew(t_{1},x_{1})dw(t_{2},x_{2})=min(t_{1},t_{2}))b(x_{1}-x_{2})$

For

$J(t,a)=\frac{\partial x}{\partial a}$

one can get by direct differentiation (6) in $a$

$dJ=J_{1}dt,\ \ dJ_{1}=-(J_{1}/\tau)dt+Jd\tilde{w}(t)$

where

$J_{1}(t,a)=\frac{\partial u}{\partial a},\ \ \ \ E\tilde{w}(t)=0,\ \
E\tilde{w}(t)^{2}=-b^{\prime\prime}(0)t$

By setting $p=\tau J_{1}/J$, introducing dimensionless time $t=t/\tau$, and
using Ito formula obtain

$dp=(-p-p^{2})dt+\sqrt{2}sdw$ $None$

where $w$ is a standard Wiener process and

$s=\sqrt{-b^{\prime\prime}(0)\tau^{3}}$ $None$

. Worth noting that if the velocity and length scales of the flow are chosen
as $U=\sqrt{b(0)\tau}$ and $L=\sqrt{-b(0)/b^{\prime\prime}(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. [4]), more exactly

$\quad\lim_{t\rightarrow S-0}p(t)=-\infty,$

where the explosion time $S$ is finite

$\quad P(S<\infty)=1$

and its expectation can be explicitly computed

$E(S)=2\pi^{2}\sqrt{z}M(z)^{2},\ \ z=s^{-4/3}/4$ $None$

where

$M(z)=\sqrt{Ai(z)^{2}+Bi(z)^{2}},\quad$

$Ai(z),Bi(z)$ are the Airy functions and $M(z)$ their modulus.

Explosiveness of $p(t)$ follows from the Feller’s criteria, [4].

To prove (9) introduce $\mu(x)$ as the mean time to explosion under condition
$p(0)=x$, then solve

$\quad L\mu(x)=-1,\quad\mu(-\infty)=0$

and set

$E(S)=\mu(\infty).$

where

$L=(-x-x^{2})\frac{\partial}{\partial x}+s^{2}\frac{\partial^{2}}{\partial
x^{2}}$

is the generator of (7). Details of computations can be found in [5].

For the purpose of studying $2D$ models the knowledge of the mean explosion
time is not enough. Introduce

$\bar{F}(t)=P(S>t)$

and let $u(t,x)$ be the same probability under condition $p(0)=x$, then

$\bar{F}(t)=u(t,\infty)$

In [6] it is shown that $u(t,x)$ satisfies the initial value problem

$\displaystyle\frac{\partial u(t,x)}{\partial t}=L_{x}u(t,x),\ \ u(0,x)=1$
$None$

To ensure uniqueness of solution of (10) we add the natural boundary
conditions assuming that the limit $u(t,\infty)$ exists similarly to
$m(\infty)$.

$u(t,-\infty)=0,\ \ u_{x}(t,\infty)=0$

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.
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.

Graphs of $\bar{F}(t)$ for few values of $s$ are shown in Fig.1

Figure 1. $P(S>t)$ versus $t$ for $s=0.5,1,2$

After explosion process $p(t)$ can be continued by starting over, i.e. by
solving same eq. (7) with the initial condition

$\lim_{t\rightarrow S+0}p(t)=\infty$

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)$.

LE for flow (7) is expressible in terms of the ergodic mean of $p$

$\lambda=\bar{p}/\tau,\ \ \
\bar{p}=\displaystyle\lim_{t\to\infty}\frac{1}{t}\int_{0}^{t}p(s)ds$

as follows from (3) and definition of $p(t)$, see [3] for details. It also can
be exactly found

Proposition 2.

The ergodic mean of $p$ is given by

$\bar{p}=-1+\frac{M^{\prime}(z)}{2\sqrt{z}M(z)}$ $None$

The statement is proven in Appendix and here we just make some comments.
Notice that the density of the stationary distribution $\pi(p)$ of $p(t)$ does
exist and is given by

$\pi(p)=Ce^{-p^{2}/2s^{2}-p^{3}/3s^{2}}\int_{-\infty}^{p}e^{y^{2}/2s^{2}+y^{3}/3s^{2}}dy$
$None$

where $C$ is a normalized constant.

One can find the following asymptotic

$\pi(p)\sim|p|^{-2},\quad p\rightarrow\infty.$

Thus, the integral for the invariant mean

$Ep=\int_{-\infty}^{\infty}p\pi(p)dp$

is formally divergent, however if the integral is meant as a Cauchy principal
value then it coincides with the RHS of (11).

For bounded functions of $p$ the ergodicity holds true with a conventional
interpretation of the integral over the invariant density.

Proposition 3

Process $p(t)$ is ergodic , i.e. for any bounded function $f(p)$

$\lim_{T\rightarrow\infty}\frac{1}{T}\int_{0}^{T}f(p(t))dt=\int_{-\infty}^{\infty}f(p)\pi(p)dp$

2\. Deterministic 2D Case

Let us assume $\mathbf{w}=0$ in (1) and proceed to dimensionless variables as
described in Introduction before eq. (5). As a result the expression for the
Jacobian (2) is found in the explicit form

$J(t)=Ps^{2}\left(1-e^{-t}\right)^{2}+Qs\left(1-e^{-t}\right)+1$

where $s$ is Stokes number defined in (5)

$\displaystyle P=\frac{\partial(u_{0},v_{0})}{\partial(a,b)},\ \
Q=\frac{\partial u_{0}}{\partial a}+\frac{\partial v_{0}}{\partial b}$

and $(a,b)$, / / ($u_{0},v_{0})$ are the initial position and velocity
respectively related by $u_{0}=u_{0}(a,b),\ \ v_{0}=v_{0}(a,b)$.

Let

$G=\\{(P,Q)|\ Q^{2}>4P\\}\bigcap\\{(P,Q)|\ P<0\ \rm{or}\ Q<0\\}$

and $S$ be the first time when $J(S)=0$. Easy to find that

$\displaystyle
S=\left\\{\begin{array}[]{lll}\displaystyle-\ln\left(1-\frac{s_{0}}{s}\right)&{\rm
if}\ \ (P,Q)\in G\ \ {\rm and}\ \ s>s_{0}\\\ \\\ \infty&\ {\rm
otherwise}\end{array}\right.$

where

$\displaystyle s_{0}=\frac{-Q-\sqrt{Q^{2}-4P}}{2P}$

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

$\frac{\partial\mathbf{u}}{\partial
t}+s\mathbf{u}\cdot\nabla\mathbf{u}+\mathbf{u}=\mathbf{0},\ \ \
\mathbf{u}(0,\mathbf{x})=\mathbf{u}_{0}(\mathbf{x})$ $None$

looses uniqueness exactly at moment $t=S$. If one defines $u(t,x)$ after $t>S$
as a solution of (13) with the initial condition

$\mathbf{u}(S+,\mathbf{x})=\mathbf{u}_{0}(\mathbf{x})$

then $\nu=1/S$ could be interpreted as the intensity of caustics in time. In
Figure 2 we show few curves $\nu(s)$ for different $s_{0}$ in order to compare
them visually with a stochastic case addressed in the next sections.

Figure 2. Dependence of the caustic rate $\nu(s)$ on Stokes number $s$ in
absence of noise.

If $S$ depends on the initial point, then its interpretation in terms of
Eulerian set up (13) is not that simple and will not be discussed here.

3\. System of Equations for Jacobian in general 2D Case

Now we return to the stochastic model (1). In few works it was pointed out
that there is a closed equation for the matrix

$\mathbf{Z}=\frac{\partial\mathbf{v}}{\partial\mathbf{a}}\left(\frac{\partial\mathbf{r}}{\partial\mathbf{a}}\right)^{-1}$

Namely, e.g. [3]

$\displaystyle
d\mathbf{Z}=-\left(\mathbf{Z}/\tau+\mathbf{Z}^{2}\right)dt+d\left(\frac{\partial\mathbf{w}}{\partial\mathbf{a}}\right)$

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
$\mathbf{r}=(x,y),\mathbf{a}=(a,b)$

$\begin{array}[]{lll}du=-(u/\tau)dt+dw^{(1)}(t,x,y),\quad
dv=-(v/\tau)dt+dw^{(2)}(t,x,y),\quad dx=udt,\quad dy=vdt\\\ \\\ x(0)=a,\ \
y(0)=b,\ \ \ \ u(0)=u_{0}(a,b),\ \ v(0)=v_{0}(a,b)\end{array}$ $None$

where

$Edw^{(i)}(t,0,0)dw^{(j)}(t,x,y)=b^{ij}(x,y)dt$

and $b^{ij}$ are entries of $\mathbf{B}$ The goal is to investigate time
behavior of the Jacobian

$J(t)\equiv\frac{\partial x}{\partial a}\frac{\partial y}{\partial
b}-\frac{\partial x}{\partial b}\frac{\partial y}{\partial a}$

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. Namely, for
dimensionless time $t=t/\tau$ denoted by the same letter, Jacobian can be
represented as

$J(t)=\exp\left(\int_{0}^{t}m_{1}(t)dt\right),$

where dimensionless random function $m_{1}(t)$ is included into the following
system

$\begin{array}[]{lll}dm_{1}&=&(-m_{1}+(m_{4}^{2}-m_{1}^{2}-m_{2}^{2}-m_{3}^{2})/2)dt+sd\eta_{1},\\\
dm_{2}&=&(-m_{1}m_{2}-m_{2})dt+sd\eta_{2},\\\
dm_{3}&=&(-m_{1}m_{3}-m_{3})dt+sd\eta_{3},\\\
dm_{4}&=&(-m_{1}m_{4}-m_{4})dt+sd\eta_{4}\end{array}$ $None$

$s$ is Stokes number defined similarly to (8) and an exact expression for it
is given in Appendix.

To define processes $\eta_{i}(t)$ we nondimensionalize $w^{(i)}(t)$ and set

$\eta_{1}=w^{(2)}_{y}-w_{x}^{(1)},\ \ \eta_{2}=w^{(2)}_{y}+w_{x}^{(1)},\ \
\eta_{3}=w^{(2)}_{x}-w_{y}^{1},\ \ \eta_{4}=w^{(2)}_{x}+w_{y}^{(1)},\ $ $None$

where the subs mean derivatives. Thus $\eta$’s are dependent Wiener (non-
standard) processes with a covariance matrix
$Ed\eta^{(i)}(t,x,y)d\eta^{(j)}(t,x,y)=\sigma^{ij}_{\eta}dt$ given by

$(\sigma^{ij}_{\eta})=\left(\begin{tabular}[]{llll}$-b^{11}_{xx}-b^{22}_{yy}+2b^{12}_{xy}$&$-b^{22}_{yy}+b^{11}_{xx}$&$b^{12}_{xx}+b^{12}_{yy}$&$b^{12}_{xx}-b^{12}_{yy}$\\\
\\\
$-b^{22}_{yy}+b^{11}_{xx}$&$-b^{11}_{xx}-b^{22}_{yy}-2b^{12}_{xy}$&$-b^{12}_{xx}+b^{12}_{yy}$&$-b^{12}_{xx}-b^{12}_{yy}$\\\
\\\
$b^{12}_{xx}+b^{12}_{yy}$&$-b^{12}_{xx}+b^{12}_{yy}$&$-b^{22}_{xx}-b^{11}_{yy}+2b^{12}_{xy}$&$-b^{22}_{xx}+b^{11}_{yy}$\\\
\\\
$b^{12}_{xx}-b^{12}_{yy}$&$-b^{12}_{xx}-b^{12}_{yy}$&$-b^{22}_{xx}+b^{11}_{yy}$&$-b^{22}_{xx}-b^{11}_{yy}-2b^{12}_{xy}$\\\
\\\ \end{tabular}\right)$ $None$

where $b^{ij}$ is a dimensionless version of $b^{ij}$ and partial derivatives
are taken at $x=0,y=0$.

4\. 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).

In this section we assume a zero initial velocity field

$\mathbf{u}_{0}(\mathbf{x})\equiv\mathbf{0}$

that implies zero initial conditions $m_{k}(0)=0,\ \ k>1$ in (15)

Model 1

Assume

$w^{(1)}(t,x,y)=w^{(1)}(t,x),\ \ w^{(2)}(t,x,y)=w^{(2)}(t,y)$ $None$

where $w^{(1)}(t,x),w^{(2)}(t,y)$ are independent and

$E\left(w^{(1)}_{x}\right)^{2}=E\left(w^{(2)}_{y}\right)^{2}$ $None$

From (16,18) it follows that

$\eta_{3}=0,\ \ \eta_{4}=0$

Then from (16,17,19) $E\eta_{1}^{2}=E\eta_{2}^{2},\ \ \ E\eta_{1}\eta_{2}=0$
and hence

$w_{1}=(\eta_{1}+\eta_{2})/2,\ \ w_{2}=(\eta_{1}-\eta_{2})/2$

are independent identically distributed Wiener processes.

Next due to the zero initial conditions $m_{3}=m_{4}\equiv 0$. Introduce

$p=\frac{m_{1}+m_{2}}{2},\ \ q=\frac{m_{1}-m_{2}}{2}$

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}$ $None$

for independent identically distributed processes $p(t)$ and $q(t)$

Model 2

In this model we assume

$w^{(1)}(t,x,y)=U(t,x+y)+V(t,x-y),\ \ w^{(2)}(t,x,y)=-U(t,x+y)+V(t,x-y)$
$None$

with independent $U$ and $V$ such that

$E\left(U^{\prime}\right)^{2}=E\left(V^{\prime}\right)^{2}$ $None$

From (16, 17, 21, 22) it follows that

$\eta_{2}=0,\ \ \eta_{4}=0$

and

$w_{1}=(\eta_{1}+\eta_{3})/2,\ \ w_{2}=(\eta_{1}-\eta_{3})/2$

are independent identically distributed Wiener processes since
$\eta_{1}=-2(U^{\prime}+V^{\prime})$ and $\eta_{3}=-2(U^{\prime}-V^{\prime})$

Thus $m_{2}=m_{4}\equiv 0$ and introducing

$p=\frac{m_{1}+m_{3}}{2},\ \ q=\frac{m_{1}-m_{3}}{2}$

arrive at the same equations (20).

Below we show simulated velocity fields at a particular time moment for both
models (Fig.3) where periodic in space forcings were used.

Figure 3. Snapshots of velocity fields for Model 1 (left) and Model 2 (right)

Let $S_{1}$ and $S_{2}$ be the first explosion moments for $p(t)$ and $q(t)$
respectively, then

$S=\min(S_{1},S_{2})$

Certainly the knowledge of $E{S_{i}}$ is not enough to find $E(S)$, but the
latter can be expressed in terms of $\bar{F}(t)=P\left(S_{1}>t\right)$ as

$ES=\int_{0}^{\infty}\bar{F}(t)^{2}dt$

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
$E{S_{i}}$ and $E(S)$ may differ in an order. Fortunately it is not the case
here since the tail of distribution of $S_{1}$ is exponential and can be
evaluated [5].

Figure 4. The theoretical upper and lower bounds for the caustic rate $\nu(s)$
as functions of Stokes number $s$ and its numerical version

Namely from Theorem 1.1 in [7] it follows that the exponential moment

$E\left(e^{\gamma S_{1}}\right)<\infty$ $None$

is finite for all $\gamma=\gamma(s)$ satisfying

$\gamma(s)<\frac{1}{g(s)}$ $None$

and

$g(s)=\max\left\\{x:\left(\int_{-\infty}^{x}\exp(y^{3}/3s^{2}+y^{2}/2s^{2})dy\right)\left(\int_{x}^{\infty}\exp(-y^{3}/3s^{2}-y^{2}/2s^{2})dy\right)\right\\}$

Hence it is reasonable to assume that

$\bar{F}(t)\geq e^{-\gamma t}$ $None$

with $\gamma$ satisfying (24). In view of this assumption

$E(S)>\frac{1}{2g(s)}$

On the other side obviously that

$E(S)<E(S_{1})$

Thus for the mean number of caustics $\nu(s)=1/E(S)$ one gets from (9)

$\displaystyle\left(2\pi^{2}\sqrt{z}M(z)^{2}\right)^{-1}<\nu(s)<2g(s)$

For small $s$ asymptotics of the lower bound and upper bound coincide and we
get

$\nu(s)\approx\displaystyle\frac{1}{4\pi s^{3/2}}e^{-\frac{1}{6s^{2}}}$

while for large $s$ the corresponding asymptotics differ just by a constant

$0.4s^{2/3}<\nu(s)<1.25s^{2/3}$

Approximation (25) is quite speculative and indeed it greatly overestimates
$\nu(s)$. It can be seen by comparing the bounds with exact curve $\nu(s)$
obtained from solving the corresponding PDE for $\bar{F}(t)$ (Section 1). The
upper and lower bounds for $\nu(s)$ are shown in Figure 4 as well as its
numerical version.

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
$\Delta\mathbf{r}=\mathbf{r}(t,\mathbf{a})-\mathbf{r}(t,\mathbf{0})$ is the
separation process for one dimensional flow (6). Hence
$|\Delta\mathbf{r}|\approx|\mathbf{a}\exp{(\lambda t)}$ for small
$|\mathbf{a}|$ and large $t$, where $\lambda$ is LE for (6). Then our claim
follows from definition (3).

In Model 2 LE cannot be found by such simple tools because the equations for
$x$ and $y$ components are not split.

5\. Isotropic forcing

Assume that in the original equations the forcing is isotropic, then its
covariances are given by [2]

$\displaystyle b^{11}(x,y)$ $\displaystyle=$ $\displaystyle
b_{N}(r)+\frac{x^{2}}{r^{2}}(b_{L}(r)-b_{N}(r)),\quad
b^{22}(x,y)=b_{N}(r)+\frac{y^{2}}{r^{2}}(b_{L}(r)-b_{N}(r)),$ $\displaystyle
b^{12}(x,y)$ $\displaystyle=$ $\displaystyle
b^{21}(x,y)=\frac{xy}{r^{2}}(b_{L}(r)-b_{N}(r)),$

where $r=\sqrt{x^{2}+y^{2}}$. Assume smooth longitudinal and normal
correlation functions

$\displaystyle b_{L}(r)$ $\displaystyle=$ $\displaystyle
b_{0}-\frac{1}{2}\beta_{L}r^{2}+O(r^{4}),$ $\displaystyle b_{N}(r)$
$\displaystyle=$ $\displaystyle b_{0}-\frac{1}{2}\beta_{N}r^{2}+O(r^{4}),$

with $b_{0},\beta_{L},\beta_{N}>0$. Introduce
$s_{L}=\sqrt{2\beta_{L}\tau^{3}},\ s_{N}=\sqrt{2\beta_{N}\tau^{3}}$

Then from (17) it can be derived that (15) takes form

$\begin{array}[]{lll}dm_{1}&=&(-m_{1}+(m_{4}^{2}-m_{1}^{2}-m_{2}^{2}-m_{3}^{2})/2)dt+s_{L}dw_{1},\\\
dm_{2}&=&(-m_{1}m_{2}-m_{2})dt+s_{N}dw_{2},\\\
dm_{3}&=&(-m_{1}m_{3}-m_{3})dt+s_{L}dw_{3},\\\
dm_{4}&=&(-m_{1}m_{4}-m_{4})dt+s_{N}dw_{4},\end{array}$

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.

6\. Conclusion and Discussion

While one dimensional stochastic models for inertial particles have been
comprehensively studied, [1,8] there was a lack of $2D$ examples with a
significant analytical advances. Here we suggested two such examples, where
the mean time to explosion in the particle density , $\mu(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. 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 $\nu(s)$ still remain. For
both models the reciprocal $\nu(s)=1/\mu(s)$, interpreted as the intensity of
caustics, increases monotonically from zero to infinity that is not a surprise
at all.

Regarding to a Riemann type equation for the underlying velocity field (4),
$\mu(s)$ can be viewed as a time scale on which a solution remains unique. It
should be recognized that the interpretation of $\mu$ and $\nu$ in terms of
the Eulerian velocity field is not clear enough except the deterministic case
with the initial velocity field $\mathbf{u}_{0}(\mathbf{x})$ for which
$|\partial\mathbf{u}_{0}\partial\mathbf{x}/|$ and $\nabla\cdot\mathbf{u}_{0}$
do not depend on $\mathbf{x}$ at all.

Appendix

A1. Proof of Proposition 2

The proof is based on the following statement

Lemma . If there exist real $m$ and smooth bounded function $V(x)$ such that

$\quad V(\infty)=V(-\infty)=0$ $None$

and

$\quad LV(x)=x-m,$ $None$

then with probability one

$\bar{p}=m$

Proof . From (A2) and Ito formula it follows that $V(p(t))$ satisfies

$\quad dV=(p-m)dt+sV^{\prime}dw.$

By integrating both sides we get

$\frac{1}{T}(V(p(T))-V(p(0))=\frac{1}{T}\int_{0}^{T}p(t)dt-m+\frac{s}{T}\int_{0}^{T}V^{\prime}(p(t))dw$
$None$

The second term on the right hand side goes to zero as $T\to\infty$ 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 $\quad V(\infty)\neq V(-\infty)$ 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)

$V(x)=\displaystyle\frac{1}{s^{2}}\int_{-\infty}^{x}\exp\left(\frac{y^{3}/3+y^{2}/2}{s^{2}}\right)\int_{y}^{\infty}\exp\left(-\frac{u^{3}/3+u^{2}/2}{s^{2}}\right)(u-m)dudy.$

for which $V(-\infty)=0$. To ensure (A1) one should set

$m=\frac{M}{N}$

where

$\begin{array}[]{lll}M=\int_{-\infty}^{\infty}\exp\left(\frac{y^{3}/3+y^{2}/2}{s^{2}}\right)\int_{y}^{\infty}\exp\left(-\frac{u^{3}/3+u^{2}/2}{s^{2}}\right)ududy\\\
\\\
N=\int_{-\infty}^{\infty}\exp\left(\frac{y^{3}/3+y^{2}/2}{s^{2}}\right)\int_{y}^{\infty}\exp\left(-\frac{u^{3}/3+u^{2}/2}{s^{2}}\right)dudy\end{array}$
$None$

The integral for $M$ is meant as Cauchy principal value. That allows for
changing the order of integration after substitution

$t=u-y,\ \ y^{\prime}=ys^{2/3}$

Integration in $y^{\prime}$ leads to

$\quad m=\frac{1}{2}\left(-1+\frac{g^{\prime}(z)}{2\sqrt{z}g(z)}\right)$

where

$\quad g(z)=\int_{0}^{\infty}\exp(zt-t^{3}/12)\frac{dt}{\sqrt{t}}.$

To complete the proof one should account for, [5],

$\quad g(x)=\pi^{-1}M(x)^{2}.$

This relation also can be derived from the fact that the both sides satisfy
the same differential equation and same initial conditions, [9].

Finally, notice that the expression for $m$ after changing the order of
integration in (A4) turns to the mean of $p(t)$ with respect to the invariant
measure $\pi(p)$ given in (12)

Proposition 2 is proven. Proposition 3 can be proven by the same arguments.

A2. Details of the computational algorithm for solving (10)

Let $h$ and $\Delta t$ be space and time respectively, $X>0$ a big enough
number, $x_{i}=-X+(i-1)h,\ \ n=(n-1)\Delta t$ the corresponding space/time
grid, $i=1,...,M,\ \ n=1,...,N$ where $h=2X/M,\ \ \Delta t=t/N$. Set
$b_{i}=\Delta t(x_{i}+x_{i}^{2}+s^{2}/h),\ \ a=\Delta ts^{2}/h$. Then a
standard Eulerian scheme for (23) is written as

$\begin{array}[]{lll}u_{i+1}(n)=u_{i}(n)+hy_{i}(n)\\\ \\\
u_{i+1}(n+1)-u_{i+1}(n)=ay_{i+1}(n)-b_{i}y_{i}(n)\end{array}$ $None$

Notice that (A5) is well posed because $u_{i}(1)=1,\ i=1,...,M,\ \ \
u_{1}(n)=0,\ \ y_{M}(n)=0,\ \ n=1,...,N$

The choice of parameters $X,h,\Delta t$ was dictated by standard conditions on
ratio $\Delta t/h$ and constraints

$\sum_{1}^{N}\bar{F}(n)\Delta t\approx\mu(s)$

where $\mu(s)=E(S)$ is the exact value of the expectation of $S$ found from
(9). In particular $\mu(1/\sqrt{2}$ =8.7735 while under our choice of
$X=7.5,\Delta t=0.001,M=61$ it was obtained that

$\sum_{1}^{N}\bar{F}(n)\Delta t\approx 8.7785$

A similar comparison for other $s$ can be seen from Figure 4

A3. Derivation of (15)

For the purpose of nondimensionalizing let us redenote the forcing on RHS of
(14) by capital letters $W^{(i)}(t,x,y),\ \ i=1,2$. Then by differentiating
(14) in $a$ and $b$ we get

$\displaystyle du_{a}=-(u_{a}/\tau)dt+x_{a}d\xi_{1}+y_{a}d\xi_{2},\quad
dv_{a}=-(v_{a}/\tau)dt+x_{a}d\xi_{1}+y_{a}d\xi_{4},\quad dx_{a}=u_{a}dt,\quad
dy_{a}=v_{a}dt$ $\displaystyle
du_{b}=-(u_{b}/\tau)dt+x_{b}d\xi_{1}+y_{b}d\xi_{2},\quad
dv_{b}=-(v_{b}/\tau)dt+x_{b}d\xi_{1}+y_{b}d\xi_{4},\quad dx_{b}=u_{a}dt,\quad
dy_{b}=v_{b}dt$

where

$\xi_{1}=\partial W^{(1)}/\partial x,\quad\xi_{2}=\partial W^{(1)}/\partial
y,\quad\xi_{3}=\partial W^{(2)}/\partial x,\quad\xi_{4}=\partial
W^{(2)}/\partial y$

and the covariance matrix of $\xi$’s is given by

$\left(\begin{tabular}[]{llll}$-B^{11}_{xx}$&$0$&$-B^{12}_{xx}$&$-B^{12}_{xy}$\\\
\\\ $0$&$-B^{11}_{yy}$&$-B^{12}_{xy}$&$-B^{12}_{yy}$\\\ \\\
$-B^{12}_{xx}$&$-B^{12}_{xy}$&$-B^{22}_{xx}$&$0$\\\ \\\
$-B^{12}_{xy}$&$-B^{12}_{yy}$&$0$&$-B^{22}_{yy}$\\\ \\\ \end{tabular}\right)$

where $B^{ij}$ is the space covariance matrix of
$\left(W^{(1)},W^{(2)}\right)$

Introduce

$\displaystyle J_{1}(t)$ $\displaystyle\equiv$ $\displaystyle\frac{\partial
x}{\partial a}\frac{\partial v}{\partial b}-\frac{\partial x}{\partial
b}\frac{\partial v}{\partial a},$ $\displaystyle J_{2}(t)$
$\displaystyle\equiv$ $\displaystyle\frac{\partial y}{\partial
a}\frac{\partial u}{\partial b}-\frac{\partial y}{\partial b}\frac{\partial
u}{\partial a},$ $\displaystyle J_{3}(t)$ $\displaystyle\equiv$
$\displaystyle\frac{\partial x}{\partial a}\frac{\partial u}{\partial
b}-\frac{\partial x}{\partial b}\frac{\partial u}{\partial a},$ $\displaystyle
J_{4}(t)$ $\displaystyle\equiv$ $\displaystyle\frac{\partial y}{\partial
a}\frac{\partial v}{\partial b}-\frac{\partial y}{\partial b}\frac{\partial
v}{\partial a},$ $\displaystyle J_{5}(t)$ $\displaystyle\equiv$
$\displaystyle\frac{\partial u}{\partial a}\frac{\partial v}{\partial
b}-\frac{\partial u}{\partial b}\frac{\partial v}{\partial a},$

Applying Ito formula we obtain

$\displaystyle dJ$ $\displaystyle=$ $\displaystyle(J_{1}-J_{2})dt,\quad~{}$
$\displaystyle dJ_{1}$ $\displaystyle=$
$\displaystyle(-J_{1}/\tau+J_{5})dt+Jd\xi_{4},$ $\displaystyle dJ_{2}$
$\displaystyle=$ $\displaystyle(-J_{2}/\tau-J{5})dt+Jd\xi_{1},$ $\displaystyle
dJ_{3}$ $\displaystyle=$ $\displaystyle(-J_{3}/\tau)dt+Jd\xi_{2},$
$\displaystyle dJ_{4}$ $\displaystyle=$
$\displaystyle(-J_{4}/\tau)dt+Jd\xi_{3},$ $\displaystyle dJ_{5}$
$\displaystyle=$
$\displaystyle(-J_{5}/\tau)dt+J_{1}d\xi_{1}+J_{2}d\xi_{2}-J_{3}d\xi_{3}-J_{4}d\xi_{4},$

To proceed to dimensionless variables assume

$W^{(i)}(t,x,y)=\sigma\tau^{1/2}w^{(i)}(t/\tau,x/L,y/L),\ \ i=1,2$

where $w^{(i)}$ are dimensionless Wiener processes, thereby in dimensionless
variables

$Edw^{(i)}(t,0,0)dw^{(j)}(t,x,y)=b^{ij}(x,y)dt$

with

$B^{ij}(x,y)=\sigma^{2}b^{ij}(x/L,y/L)$

Introduce a Stokes number similarly to (8)

$s=\sigma\tau^{3/2}/L$

The result is

$\begin{array}[]{lll}dJ&=&I_{1}dt,{}\\\
dI_{1}&=&(-I_{1}+2J_{5})dt+sJd\eta_{1},\\\ dI_{2}&=&-I_{2}dt+sJd\eta_{2},\\\
dI_{3}&=&-I_{3}dt+sJd\eta_{3},\\\ dI_{4}&=&-I_{4}dt+sJd\eta_{4}\end{array}$
$I_{1}=\tau(J_{1}-J_{2}),\quad I_{2}=\tau(J_{1}+J_{2}),\quad
I_{3}=\tau(J_{4}-J_{3}),\quad I_{4}=\tau(J_{3}+J_{4})$

where

$\eta_{1}=(\xi_{4}-\xi_{1})/h,\eta_{2}=(\xi_{4}+\xi_{1})/h,\eta_{3}=(\xi_{3}-\xi_{2})/h,\eta_{4}=(\xi_{3}+\xi_{2})/h,\
\ \ h=\sigma\tau^{1/2}$

We do not need a stochastic equation for $J_{5}$ since it can be found from
the following equation

$I_{2}^{2}-I_{1}^{2}-I_{4}^{2}+I_{3}^{2}=-4JJ_{5},$ $None$

which follows from

$J_{1}J_{2}+JJ_{5}=J_{3}J_{4}.$

We can reduce the number of freedom degrees by introducing
$m_{k}=I_{k}/J,\quad k=1,...,5$ where $m_{5}$ is eliminated by means of (A6).
The result is system (15)

References

1\. K. Gustavsson and B. Mehlig, Statistical models for spatial patterns of
inertial particles in turbulence, _arXiv:1412.4374v1 [physics.flu-dyn]_ ,
2014.

2\. A. S. Monin and A.M. Yaglom, ”Statistical Fluid Mechanics: Mechanics of
Turbulence”, MIT Press, Cambridge, MA, 1975.

3\. L. I. Piterbarg, The top Lyapunov exponent for a stochastic flow modeling
the upper ocean turbulence. _SIAM J. Appl. Math._ , 62, 777-800, 2001.

4\. I. Karatzas and S.E. Shreve, _Brownian Motion and Stochastic Calculus_ ,
Springer-Verlag., NY-Berlin-Heidelberg, 1991.

5\. B. I. Halperin, Green’s functions for a particle in a one-dimensional
random potential, _Phys. Rev.,_ 139, A104-117, 1965.

6\. I. Karatzas I. and S.E. Shreve, Distribution of the time to explosion for
one-dimensional diffusions _Probab. Theory Relat. Fields_ , 164, 1027-1069,
2016.

7\. O. Loukianov, D. Loukianova, and S. Song, Spectral gaps and exponential
integrability of hitting times for linear diffusions, _Ann. Inst. H. Poincaré
Probab. Statist. Volume 47, Number 3 , 679-698.,_ , 2011.

8\. G. Falkovich , S. Musacchio , L. Piterbarg , and M. Vucelja (2007)
Inertial particles driven by a telegraph noise, _Physical Review E_ , 76, 2007

9\. M. Abramowitz and I. A. Stegun, _Handbook of Mathematical Functions_ Dover
Publications, INC., New York, 1970