2.1. Evolution of geometric invariants

Equation (1.1) describes a nine-dimensional dynamical system in three spatial dimensions. However, for an isotropic system only geometric invariants, i.e. combinations of the coefficients that are invariant under rotation of the coordinate system, carry physical information (Lumley Reference Lumley2007, appendix A2.6). These invariants amount to the coefficients $c_k$ of the characteristic polynomial (2.2) of $\textsf{$\boldsymbol{\sigma}$}$ (cf.  Lang Reference Lang1993; Garibaldi Reference Garibaldi2004; Lawson & Dawson Reference Lawson and Dawson2015),

(2.2)\begin{equation} P(\lambda) = \det( \lambda \boldsymbol{\mathsf{I}}_3 - \textsf{$\boldsymbol{\sigma}$} ) = \lambda^3 - X \lambda^2 + Y \lambda - Z \lambda^0, \end{equation}

where $\boldsymbol{\mathsf{I}}_3$ is the $3\times 3$ identity matrix, and

(2.3a–c)\begin{equation} X = \text{Tr}(\textsf{$\boldsymbol{\sigma}$}) , \quad Y = \frac{1}{2} \left(\text{Tr}(\textsf{$\boldsymbol{\sigma}$})^2 - \text{Tr}(\textsf{$\boldsymbol{\sigma}$}^2) \right) , \quad Z = \det(\textsf{$\boldsymbol{\sigma}$}) . \end{equation}

The time evolution of $X$ is found by taking the trace of (1.1). Due to the incompressibility of the fluid the noise term $\boldsymbol{\mathsf{s}}$ is trace-less and vanishes. Moreover, $\text {Tr}(\textsf{$\boldsymbol{\sigma}$} ^2)$ can be expressed in terms of $X$ and $Y$, yielding

(2.4a)\begin{equation} \dot X ={-} X (1+X) + 2 Y . \end{equation}

The time derivative of $Y$ amounts to $\dot Y = X \dot X - ({1}/{2}) \, ({{{\mathrm {d}}}}/{{{\mathrm {d}}} t}) \text {Tr}(\textsf{$\boldsymbol{\sigma}$} ^2)$, and the time derivative of $\text {Tr}(\textsf{$\boldsymbol{\sigma}$} ^2)$ can be worked out by multiplying (1.1) by $\textsf{$\boldsymbol{\sigma}$}$ and taking the trace

(2.4b)\begin{align} \dot Y &= X \dot X + \left( \text{Tr}( \textsf{$\boldsymbol{\sigma}$}^2 ) + \text{Tr}( \textsf{$\boldsymbol{\sigma}$}^3 ) - \text{Tr}(\boldsymbol{\mathsf{s}} \textsf{$\boldsymbol{\sigma}$}) \right) \nonumber\\ &={-}X^2 (1+X) + 2 X Y + X^2 - 2 Y + X (X^2-2Y) - X Y + 3 Z - \text{Tr}(\boldsymbol{\mathsf{s}} \textsf{$\boldsymbol{\sigma}$}) \nonumber\\ &={-}Y (2+X) + 3 Z - \text{Tr}(\boldsymbol{\mathsf{s}} \textsf{$\boldsymbol{\sigma}$}) . \end{align}

In the second step we used the Cayley–Hamilton theorem, $P(\textsf{$\boldsymbol{\sigma}$} ) = 0$, to express $\text {Tr}( \textsf{$\boldsymbol{\sigma}$} ^3 )$ in terms of $X$, $Y$ and $Z$.

To derive the temporal evolution of $Z$ we multiply (1.1) with $\textsf{$\boldsymbol{\sigma}$} ^{-1}$, take the trace and use Jacobi's formula:

(2.4c)\begin{equation} \text{Tr}( \log(\textsf{$\boldsymbol{\sigma}$})) = \log(\det(\textsf{$\boldsymbol{\sigma}$})). \end{equation}

Hence, we find

(2.4d)\begin{equation} \dot Z = Z \frac{{{\mathrm{d}}}}{{{\mathrm{d}}} t} \log Z = Z \text{Tr}\!\left( \textsf{$\boldsymbol{\sigma}$}^{{-}1} \frac{{{\mathrm{d}}}}{{{\mathrm{d}}} t} \textsf{$\boldsymbol{\sigma}$} \right) ={-}Z \left(3+X\right) + \text{Tr}\left( \boldsymbol{\mathsf{s}} \textsf{$\boldsymbol{\sigma}$}^* \right) , \end{equation}

where $\textsf{$\boldsymbol{\sigma}$} ^*= Z \textsf{$\boldsymbol{\sigma}$} ^{-1}$ is the conjugate matrix of $\textsf{$\boldsymbol{\sigma}$}$.

As a final remark we note that two-dimensional turbulence may be considered as flow in a plane where all perpendicular velocity components vanish. In that case also their gradients vanish, and consequently also the determinant $Z$ of the gradient matrix. Moreover, by a straightforward calculation, one verifies that $Y$ is the determinant of the non-vanishing $2\times 2$ sub-matrix of $\textsf{$\boldsymbol{\sigma}$}$. Therefore, (2.4a) and (2.4b) with $Z=0$ provide the description of inertial particles in two-dimensional turbulence. In the one-dimensional case (2.4a) with $Y = 0$ and $Z = 0$ reduces to a well-studied one-dimensional model to study caustics (Gustavsson & Mehlig Reference Gustavsson and Mehlig2016).

The dynamical system (2.4) constitutes a nonlinear stochastic differential equation in a Lagrangian point of view when $\text {Tr}( \boldsymbol{\mathsf{s}} \textsf{$\boldsymbol{\sigma}$} )$ and $\text {Tr}( \boldsymbol{\mathsf{s}} \textsf{$\boldsymbol{\sigma}$} ^* )$ are interpreted as noise terms acting on the $Y$ and $Z$ dynamics, respectively. These coupling terms model the inertial effects of the turbulent dynamics, inducing the deviation of the particle from the fluid dynamics.

3.1. The flow upon approaching a caustic

The deterministic dynamical system (3.1) has fixed points at

(3.6a–d)\begin{equation} P_0 = (0,0, 0) , \quad P_1 = ({-}1,0, 0) , \quad P_2 = ({-}2,1, 0) , \quad P_3 = ({-}3,3,-1) . \end{equation}

There is a remarkable correspondence between triples of these points and the constants $C_0, \ldots, C_3$:

$C_0 = 0$ determines trajectories in the plane $Z_0 = 0$ spanned by $P_0$, $P_1$, $P_2$;

$C_1 = 0$ determines trajectories in the plane $Y_0 = -3 Z_0$ spanned by $P_0$, $P_1$, $P_3$;

$C_2 = 0$ determines trajectories in the plane $X_0 = - 2 Y_0 -3 Z_0$ spanned by $P_0$, $P_2$, $P_3$;

$C_3 = 0$ determines trajectories in the plane $Y_0 = -1 -X_0 -Z_0$ spanned by $P_1$, $P_2$, $P_3$.

The explicit solution (3.2) implies that these are invariant planes of the dynamics. Trajectories that start in one of the planes will stay in the plane, and trajectories that start in the intersection of any pair of these planes have a dynamics on a straight line defined by the intersection of the planes. Therefore, also the straight lines through any two of the fixed points are invariants of motion. Linear stability analysis reveals that the flow along these lines is oriented such that $P_i$ has $i$ unstable directions. In particular, $P_0$ is linearly stable.

For $C_3 \neq 0$ and $t \to \infty$ the denominators of (3.2) increase faster than the numerators. Hence, $P_0$ is a global attractor of the dynamics which renders the particle velocity field incompressible i.e. where particles follow the surrounding fluid. Almost all trajectories will approach the origin. Trajectories in the half space that contains the origin can approach zero without divergence. Trajectories in the plane $C_3 = 0$ will either approach $P_1$ or $P_2$. Trajectories on the opposite side of that plane will encounter a root of the denominator at some finite time $t_c$, that leads to a flip of the signs of $X$, $Y$ and $Z$. This finite-time singularity amounts to a steepening of the slope. The velocity gradients diverge and change sign. Similar events happen in wave breaking, and in the formation of caustics in optical reflections (Berry Reference Berry1976; Beven Reference Beven2019). An initial condition $(X_0, Y_0, Z_0)$ will undergo $0 \leq \nu \leq 3$ caustics when the polynomial has $\nu$ real roots with ${E_c} = \exp (t_c) > 1$, and it has undergone $0\leq \nu _p \leq 3-\nu$ caustics when the polynomial has $\nu _p$ roots with $0<{E_c}<1$.

There are regions where (1)-caustics occur. Furthermore, there is a phase-space region where (2)-caustics occur and a region where (3)-caustics occur. (3)-caustics are thereby defined by phase-space trajectories that diverge three times before they approach the global attractor $P_0$.

The topological constraints on the flow stipulate that the domains of different types of trajectories are separated by the four invariant planes spanned by triples of three fixed points, and by the manifolds where cubic polynomial has double and triple roots.

The solutions, (3.2), of the deterministic dynamics, (3.1), are ratios of polynomials in the variable ${{\mathrm {e}}}^t$. In the generic case the roots of the numerator and denominator are distinct. In that case there are numbers $\alpha$ and $\beta$ such that the approach to the singularity is always linear in leading order,

(3.7)\begin{equation} \lvert X(t) \rvert \gg 1 \quad \Rightarrow \quad Y(t) \sim \alpha X(t) \quad \land \quad Z(t) \sim \beta X(t) . \end{equation}

Far away from the origin the evolution of the trajectory towards a caustic follows a straight line to infinity with $X(t) \to -\infty$ with some values $\alpha, \beta \in \mathbb {R}$. Interestingly, there are infinite many possibilities of a caustic to occur because $\alpha$ and $\beta$ are real parameters. Then it returns from the exactly opposite position in the three-dimensional phase-space. It still follows (3.7) with the same numbers $\alpha$ and $\beta$, and $X(t) \gg 1$ and decaying. Subsequently, the phase-space trajectory either goes through another caustic, where $\alpha$ and $\beta$ take other values, or it settles down to the origin. A trajectory can undergo at most three caustics, that correspond to the three roots of the denominator of (3.2).

The deterministic dynamics (3.1) always ends up asymptotically in the global attractor $P_0$ and is therefore not chaotic.

3.2. Caustics in the two-dimensional model

It is instructive to note that the $(X,Y)$-plane is an invariant manifold of the three-dimensional dynamics (3.1), and that this manifold is attractive as long as $X > -3$. Basic features of the dynamics can therefore be understood already in the two-dimensional dynamics where it is much easier to visualise the structure of the flow. For $Z \equiv 0$ the two-dimensional deterministic dynamics

(3.8a)\begin{gather} \dot{X}(t) ={-}X(t) \left( 1 + X(t) \right) + 2 \, Y(t), \end{gather}

(3.8b)\begin{gather}\dot{Y}(t) ={-}Y(t) \left( 2 + X(t) \right) \end{gather}

has solutions

(3.9a)\begin{gather} X(t)= \frac{-2c_0- c_1 \, {\rm e}^t }{ c_0 + c_1 \, {\rm e}^t + c_2 \, {\rm e}^{2t} }, \end{gather}

(3.9b)\begin{gather}Y(t) = \frac{c_0}{ c_0 + c_1 \, {\rm e}^t + c_2 \, {\rm e}^{2t} } \end{gather}

with $c_0$, $c_1$, $c_2$ determined by the initial condition $(X_0,Y_0)$ of a phase-space trajectory at time $t=0$,

(3.10)\begin{equation} \left. \begin{aligned} c_0 & = Y_0,\\ c_1 & ={-} (X_0+2 Y_0)\\ c_2 & = 1+X_0+Y_0 . \end{aligned} \right\} \end{equation}

This dynamics has fixed points at

(3.11a–c)\begin{equation} Q_0 = (0,0) , \quad Q_1 = ({-}1,0), \quad Q_2 = ({-}2,1) , \end{equation}

i.e. the fixed points $P_0$, $P_1$, and $P_2$ of (3.6a–d) with $Z=0$. Accordingly, the invariant manifolds of the two-dimensional system are straight lines through every pair of these fixed points, and the flow along these lines is orientated such that $Q_i$ has $i$ unstable directions.

During a caustic the phase-space trajectories are to leading order linear functions with slope $\gamma$ such that

(3.12)\begin{equation} Y(t) \sim \gamma X(t) , \end{equation}

and there are periodic boundary conditions in the variables $X$ and $Y$ during a caustic at infinity.

The flow field, the invariant lines and the fixed points are shown in figure 1. This figure nicely shows that the origin is an attractor of the flow. Small finite excitations from zero decay right away, as long as they amount to initial conditions to the upper right of the manifold spanned by $Q_1$ and $Q_2$,

(3.13)\begin{equation} Y ={-}1-X \end{equation}

that is marked by a solid red line. However, the trajectory makes a large excursion when their initial conditions are located below the red line (3.13). In particular, it must then go through a caustic before it can return to the origin. Different trajectories with initial conditions below the red line (3.13) are shown in figure 2. One can see that they have to perform a caustic event until the relaxation towards the global attractor sets in. At the global attractor the particle velocity fields is incompressible and no further caustic events can occur.

Figure 1.

The flow field, the invariant lines and the fixed points. Arrows show the direction of the flow field of (3.8) in the $(X,Y)$-plane. The speed of the flow is provided by a colour coding with a label provided by the bar to the right. The three fixed points $Q_0$, $Q_1$ and $Q_2$ are marked by a red $\boldsymbol {+}$, a green $\blacksquare$ and a blue $\bullet$, respectively. The invariant lines that proceed through all pairs of fixed points are marked by solid lines where the colours serve for a better visualisation.

Figure 2.

Different types of caustics in the two-dimensional model. The dashed, dash-dotted and dotted line show trajectories featuring a ($1^-$)-caustic, a ($1^+$)-caustic and a 2-caustic, respectively.

Similar to the three-dimensional system one can classify the trajectories based on the roots of the common denominator $Y_0 - (x_0 + 2 Y_0) \exp (t) + (1 + X_0 + Y_0) (\exp (t))^2$ of the solutions (3.9). The denominator is a quadratic function in $\exp (t)$ with roots ${E_c}=\exp (t_c)$. The polynomial has a double root at $Y_0 = X_0^2/2$ (yellow line in figure 2). Trajectories undergo $0 \leq \nu \leq 2$ caustics if the denominator has $\nu$ roots with ${E_c} > 1$, and they experienced $0 \leq \nu _p \leq 2-\nu$ caustics in their past when the denominator has $\nu _p$ roots with $0<{E_c}<1$. A straightforward, even though tedious calculation provides a classification that can also be inferred for the inspection of the flow (figure 1).

1. (i) For initial condition with $X_0 < -1$ and $Y_0 = 0$ there are roots with ${E_c} = 0$ and ${E_c}>1$, respectively. These initial conditions lie on trajectories that started close to $Q_1$, and they undergo a single caustic where they diverge towards the left. They return from the right, and converge towards $Q_0$.

2. (ii) For initial condition with $X_0 < -1$ and $0 < Y_0 < -1-X_0$ there are roots with ${E_c} < 0$ and ${E_c}>1$, respectively. These initial conditions lie on trajectories that started close to $Q_2$, and they undergo a single caustic where they diverge towards the upper left. They return from the lower right, and converge towards $Q_0$. These trajectories are said to exhibit ($1^+$)-caustics. An example is provided by the dash-dotted line in figure 2.

3. (iii) Initial conditions with $X_0 > -1$ and $0 \geq Y_0 > -1-X_0$ amount to post-caustic locations on the trajectories described in (i) for $Y_0=0$ and (ii) for $Y_0>0$, respectively.

4. (iv) Initial conditions with $X_0 < -2$ and $Y_0 = -1-X_0$ proceed along the red line, originating at $P_2$, diverging to the upper right and then approaching $P_1$ along the red line.

5. (v) For initial condition with $-1-X_0 < Y_0 < X_0^2/2$ there are two roots with ${E_c}>1$. These initial conditions lie on trajectories that started close to $Q_2$, undergo a caustic where they diverge towards the upper left. They return below the red line from the lower right, undergo another caustic towards the lower left and then they converge towards $Q_0$. These trajectories are said to exhibit (2)-caustics. An example is provided by the dotted line in figure 2.

6. (vi) Initial conditions with $Y_0<0$ and $Y_0 < -1-X_0$ amount to points in between the two caustics the former trajectories. Trajectories initialised at these points perform a ($1^-$)-caustic (dashed line in figure 2). The vast majority of noise-induced caustics is of this type.

7. (vii) For $X_0<-2$ and $Y_0 = X_0^2/2$ the trajectories are initiated on the line where the denominator has a double root. They evolve towards the upper left along the parabola indicated as solid orange line in figure 2, and subsequently they approach $Q_2$ along the positive branch of the parabola. This special trajectory separates the region where initial conditions perform two caustics, and where they return to the origin without undergoing a caustic.

8. (viii) All other initial conditions approach directly $Q_0$ without undergoing caustics.

3.3. Caustics in the one-dimensional model

The explicit solutions of the deterministic evolution imply that caustics involve a divergence of the divergence, $X=\boldsymbol {\nabla } \boldsymbol {\cdot } {{\boldsymbol {v}}}$, of the particle velocity field, ${{\boldsymbol {v}}}$ (Wilkinson & Mehlig Reference Wilkinson and Mehlig2005; Wilkinson et al. Reference Wilkinson, Mehlig and Bezuglyy2006). A first impression on the aerosol dynamics can then be obtained by considering only the one-dimensional dynamics along the invariant manifold $Y_0=Z_0=0$ (Gustavsson et al. Reference Gustavsson, Meneguz, Reeks and Mehlig2012; Gustavsson & Mehlig Reference Gustavsson and Mehlig2013,  Reference Gustavsson and Mehlig2016)

(3.14)\begin{equation} \left. \begin{gathered} \dot X(t) ={-} X(t) \left( 1 + X(t) \right) ={-} \frac{{{\mathrm{d}}}}{{{\mathrm{d}}} X} \varPhi( X(t) )\\ \text{with}\ \varPhi(X) = \frac{X^2}{6}( 3 + 2 X ) . \end{gathered} \right\} \end{equation}

For the initial condition $X(0) = X_0$, (3.14) has the solution

(3.15)\begin{equation} X(t) = \frac{ X_0 }{ (1+X_0) \, {{\mathrm{e}}}^t - X_0 } \end{equation}

and the system undergoes a caustic when ${E_c} = \exp (t_c) = X_0 / (1+X_0) > 1$, i.e. for initial condition $X_0 < -1$ where $X=0$ corresponds to the unstable point of the one-dimensional dynamical system. For $X_0>0$ the trajectory underwent a caustic in the past, and for $-1< X_0<0$ the trajectories relax directly towards the stable fixed point at the origin.

At this point we provide another property for the analysis of caustics, the caustic duration $T_c$. As the formulas are quite easy for the one-dimensional system we provide it first here. A similar and straightforward calculation holds true also in the two- and three-dimensional cases. An estimation for the duration of caustic events is especially interesting for the next section where we analyse the effects of applying noise to the deterministic dynamics, i.e. (2.4). The duration of a caustic will be defined as the time $T_c$ where the particle proceeds in a region $X \not \in [-X_D, X_D]$ where the forcing due to the noise is negligible as compared with the deterministic force due to the steep decay of the potential $\varPhi (X)$. From (3.15) one infers that

(3.16)\begin{equation} X_D = \frac{ -X_D }{ (1-X_D) \, {{\mathrm{e}}}^{\langle T_c \rangle} + X_D } \quad \Rightarrow \quad \langle T_c \rangle = \ln \frac{ -1 -X_D }{ 1 - X_D } . \end{equation}

Here, $X_D$ is a user-defined parameter and typically values are $X_D > 1$. In a case with noise this parameter will depend on the noise strength of the turbulent flow field. It can be determined by analysing when the phase space trajectories approach straight lines in the mean, i.e. when trajectories during a caustic stay in a cylinder of some fixed diameter $\delta$. As a remark, in the two- and three-dimensional cases one finds curves $Y_D(X_D)$ and planes $Z_D(X_D,Y_D)$, respectively.

4.1. Frequency of caustics in the one-dimensional model

The one-dimensional model reduces to the white-noise assumption from earlier works, i.e. Wilkinson & Mehlig (Reference Wilkinson and Mehlig2003) and Gustavsson & Mehlig (Reference Gustavsson and Mehlig2016). They analysed

(4.2)\begin{equation} \dot{X} ={-} X (1+X) + \sqrt{2 D} \epsilon(t) \end{equation}

for delta-correlated white noise $\epsilon (t)$ with noise amplitude $D$. In their work they analysed the one-dimensional case with great detail up to finding exact solutions for the Lyapunov exponents which describe the noised dynamics. However, for higher dimensions the studied dynamics in the white noise limit have not yet been exactly solved. The difficulty is that the multi-dimensional deterministic drift part of the resulting set of Langevin equations is neither potential nor solenoidal. Therefore, general mathematical frameworks cannot be applied.

As the rate $j_c$ at which caustics occur is of importance for our work, we present a quick derivation of its asymptotic for large noise. We start by transforming the one-dimensional dynamics (4.2) into a steady-state Fokker–Planck equation (Gustavsson & Mehlig Reference Gustavsson and Mehlig2016),

(4.3)\begin{equation} -j_c(D) = X \left(1+X\right) \rho_s(X) + D \partial_X \rho_s(X). \end{equation}

Here, $\rho _s$ corresponds to the probability density in the steady state, and $j_c$ determines the steady-state probability current which corresponds to the rate at which caustics occur. Multiplying (4.3) with the integrating factor $\exp ( {X^3}/{3D} + {X^2}/{2D} )$ provides the formal solution

(4.4)\begin{equation} \rho_s(X) ={-} \frac{j_c}{D} \exp\left({-\frac{X^3}{3D}-\frac{X^2}{2D}}\right) \int_{-\infty}^{X} \exp\left({\frac{\hat{X}^3}{3D}+\frac{\hat{X}^2}{2D}}\right) \, {{\mathrm{d}}} \hat{X} . \end{equation}

The density, $\rho _s(X)$, is normalised over the real numbers, $X \in \mathbb {R}$. Hence, one can integrate (4.4) and derive an equation for $j_c$

(4.5)\begin{equation} -j_c(D) = \frac{D}{\int_{-\infty}^{\infty} \exp\left({-\dfrac{X^3}{3D}-\dfrac{X^2}{2D}}\right) \int_{-\infty}^{X} \exp\left({\dfrac{\hat{X}^3}{3D}+\dfrac{\hat{X}^2}{2D}}\right) \, {{\mathrm{d}}} \hat{X} \, {{\mathrm{d}}} X }. \end{equation}

In the large-noise case one can solve the integrals in (4.5) by noting that the $X^2$-term in the exponent are sub-dominant for $D\gg 1$, and adopting the substitution

(4.6a,b)\begin{equation} y = \frac{X}{D^{1/3}}, \quad \hat{y} = \frac{\hat{X}}{D^{1/3}}. \end{equation}

This turns (4.5) into

(4.7)\begin{align} -j_c(D) &= \dfrac{D^{1/3}}{ \int_{-\infty}^{\infty} \exp\left({-\dfrac{y^3}{3} -\dfrac{y^2}{2D^{1/3}}}\right) \int_{-\infty}^{y} \exp\left({\dfrac{\hat{y}^3}{3} +\dfrac{\hat{y}^2}{2D^{1/3}}}\right) \, {{\mathrm{d}}}\hat{y} \, {{\mathrm{d}}} y }\nonumber\\ &\sim \dfrac{D^{1/3}}{ \int_{-\infty}^{\infty} \int_{-\infty}^{y} \exp\left({\dfrac{\hat{y}^3-y^3}{3} }\right) \, {{\mathrm{d}}}\hat{y} \, {{\mathrm{d}}} y} \quad \text{for } D \gg 1. \end{align}

Next we observe that $\hat {y}^3-y^3=(\hat {y}-y)(y^2+y\hat {y}+\hat {y}^2)$ and introduce the substitution $z = \hat {y}-y$ and $\hat {z} = y+z/2$. After some arithmetic one finds that

(4.8)\begin{equation} -j_c(D) \sim \dfrac{D^{1/3}}{ \int_{-\infty}^{0} {{\mathrm{e}}}^{{z^3}/{12}} \ \int_{-\infty}^{\infty} \,{{\mathrm{e}}}^{ z \hat{z}^2 } \, {{\mathrm{d}}}\hat{z} \, {{\mathrm{d}}} z } \sim \dfrac{D^{1/3}}{ \sqrt{\rm \pi} \int_{-\infty}^{0} \, {{\mathrm{e}}}^{{z^3}/{12}} \sqrt{-z} \, {{\mathrm{d}}} z } \quad \text{for } D \gg 1 . \end{equation}

Substituting $x=-z^3/12$ turns the integral in the denominator into a Gamma function and provides the final result

(4.9)\begin{equation} -j_c(D) \sim \dfrac{3^{5/6}}{2^{1/3} {\rm \pi}^{1/2} \varGamma\left(\dfrac{1}{6}\right)} D^{1/3} \quad \text{for } D \gg 1 . \end{equation}

This result is compared with numerical results in § 4.3.

4.2. Caustics in the two- and three-dimensional models

The stochastic system (2.4) constitutes a noise-driven excitable system (Pikovsky & Kurths Reference Pikovsky and Kurths1997; Lindner & Schimansky-Geier Reference Lindner and Schimansky-Geier1999; Anishchenko et al. Reference Anishchenko, Astakhov, Neiman, Vadivasova and Schimansky-Geier2003) An exact solution of this dynamics lies beyond the scope of the present paper. Hence, we discuss caustics in the two- and three-dimensional models based only on numerical work, and we start with the two-dimensional sub-system

(4.10a)\begin{gather} \dot{X} ={-}X(1+X) + 2Y, \end{gather}

(4.10b)\begin{gather}\dot{Y}={-}Y(2+X) + \sqrt{2 D_1} \ {\epsilon}_1(t). \end{gather}

The position of the critical manifold (red line in figure 2) suggests that the closest point to cross the manifold lies to the lower left of the origin. Inspection of the flow field in figure 1 revels that this is also the region that is most easily reached by trajectories. Trajectories that start in the vicinity of this region will perform a ($1^-$)-caustic. Hence, we expect that these should most frequently be observed. This is in line with the findings of numerical integration of the system. The ($1^-$)-caustics occur in over $95$% of the cases. A typical caustic of this type is shown by the dashed line in figure 3. However, the other types of caustics can also be observed even if only rarely. The dotted line in figure 3 shows a $(2)$-caustic.

Figure 3.

Two caustics observed in a numerical simulation of (4.10) for noise strength $D=10$. Approximately $95\,\%$ of the caustics observed in this case are ($1^-$)-caustics (black dashed line). Most of the others are ($1^+$)-caustics and rarely one observes a (2)-caustic (black dotted line). The grey band demonstrates the region where the trajectories evolve in a tube of width $\delta$ during a caustic. Here we choose $\delta = 1$.

For the three-dimensional model (2.4), the critical manifold is spanned by the fixed points $P_1, P_2, P_3$, and the closest point to cross this manifold is situated in the $Z=0$ plane right at the position identified in the two-dimensional model. As long as $D_2 < D_1$ the dynamics stays close to the $Z=0$ plane and one observes a very similar dynamics as in the two-dimensional system; with the rare exception of a (3)-caustic. When $D_2 \geq D_1$ more and more trajectories exhibit a non-trivial dynamics in the $Z$ coordinate.

4.3. Rate of occurrence

The rate at which caustics occur is numerically estimated as the number of observed caustics $N$ divided by the total simulation time, $T_{{sim}}$,

(4.11)\begin{equation} - j_c(D) \simeq \frac{ N }{ T_{{sim}} } . \end{equation}

In all simulations we limit $N$ to 5000 caustics. To reach this number we perform simulations for noise strengths $D \geq 10$. The results of the simulations are shown in figure 4. For the one-dimensional model ($\blacksquare$) the simulation results lie right on our analytical prediction (4.9), i.e. a power law with exponent $1/3$. For the two-dimensional model ($\times )$ and the three-dimensional model with $D_1 \geq D_2$ ($\bigcirc$, $\blacklozenge$) we observe the same power-law dependence, but now with a smaller exponent of approximately $0.2$. When $D_1 < D_2$ in the three-dimensional model, the growth is even slower ($\bigstar$). In line with the assessment of the qualitative features of the trajectories of § 4.2 the two-dimensional sub-system (4.10) provides a good approximation of the rate of occurrence of caustics in the full three-dimensional dynamics (2.4) as long as $D_1 \gtrsim D_2$.

Figure 4.

The rate at which caustics occur in simulations. From top to bottom the data points indicate the different models. Here $\blacksquare$ shows data for the one-dimensional model, $\times$ shows data for the two-dimensional model, $\bigcirc$ shows data for the three-dimensional model with $D = D_1 = D_2$, $\blacklozenge$ shows data for the three-dimensional model with increasing $D_1=D$ and $D_2$ fixed at $10$ and $\bigstar$ shows data for the three-dimensional model with increasing $D_2=D$ and $D_1$ fixed at $10$. The blue line shows the prediction (4.9), which fits perfectly with the numeric data of the one-dimensional model. The red line corresponds to the power law $0.2 \ D^{1/5}$. This line provides a very good fit of the data of the two- and three-dimensional models where the noise on the $Z$-dynamics is sub-dominant.

The different exponents of the power laws imply that the noise strength of the effective noise in the one-dimensional model must not be interpreted as a turbulent noise strength.

4.4. Probability to observe caustics

The analysis of the previous subsection reveals that caustics occur much less frequently than expected based on the analysis of the one-dimensional model. However, the linear phase-space dependence established in (3.7) provides an opportunity to effectively identify runaway trajectories, even while the divergence of the velocity is not yet or no longer close to its divergence. In this novel approach the relative probability $P_{{obs}}$ to identify a caustic when observing a particle at random amounts to the sum over all caustic times $T_c$ of the $N$ caustics observed in a simulation divided by the total time of the simulation, $T_{{sim}}$. As suggested in § 3.3 the duration of the caustic $T_c$ is provided here by the time at which the phase-space trajectory is locked at a straight line within an $\delta$-tube so that (3.7) holds to a good approximation, and the noise is sub-dominant in the dynamics. This is demonstrated by the grey bands in figure 3 that indicate the $\delta =1$ tubes for the displayed caustics. Based on this algorithm we define

(4.12)\begin{equation} P_{{obs}} = \frac{1}{T_{{sim}}} \sum_{i=1}^{N} T_{c,i} = \frac{N}{T_{{sim}}} \frac{1}{N} \sum_{i=1}^{N} \, T_{c,i} \simeq j_c \langle T_c \rangle. \end{equation}

In the last step of this equation we used that $N/T_{{sim}}$ converges to the rate, $j_c$, of the occurrence of caustics, and that the total duration of all caustics divided by the number of caustics amounts to the average duration $\langle T_c \rangle$ of a caustic.

The simulation results are reported in figure 5. Overall, the probability to observe a caustic is located in the single-digit percentage range, and it increases in all cases as the noise amplitudes become larger. For the one-dimensional model the numerical data ($\blacksquare$) lie right on the prediction (solid blue line) that is obtained as the product of (4.9) and (3.16). The probability to observe a caustic increases roughly logarithmically with the noise strength $D$. The data for the two-dimensional dynamics ($\times$) follow the same trend. However, in this case the probability is smaller by an ample factor of two. The data for the three-dimensional dynamics with $D=D_1>D_2$ ($\blacklozenge$) takes the same values as for the two-dimensional dynamics. For $D = D_1 = D_2$ the probabilities ($\bigcirc$) are smaller, and the data for $D=D_2 > D_1$ take the smallest probabilities ($\bigstar$). This shows that a faithful discussion of the probability to observe caustics must be based on a careful analysis of the noise terms.

Figure 5.

The fraction of time, $P_{obs}$, that a particle spends in a caustic, i.e. it approaches infinity in a path that fits into a tube of diameter $\delta$. See figure 4 for the description of symbols. The solid blue line shows the theory curve of the observation probability in the one-dimensional case, which fits the data nicely.