Caustic formation in a non-Gaussian model for turbulent aerosols

Caustics in the dynamics of heavy particles in turbulence accelerate particle collisions. The rate $\mathscr{J}$ at which these singularities form depends sensitively on the Stokes number St, the non-dimensional inertia parameter. Exact results for this sensitive dependence have been obtained using Gaussian statistical models for turbulent aerosols. However, direct numerical simulations of heavy particles in turbulence yield much larger caustic-formation rates than predicted by the Gaussian theory. In order to understand possible mechanisms explaining this difference, we analyse a non-Gaussian statistical model for caustic formation in the limit of small St. We show that at small St, $\mathscr{J}$ depends sensitively on the tails of the distribution of Lagrangian fluid-velocity gradients. This explains why different authors obtained different St-dependencies of $\mathscr{J}$ in numerical-simulation studies. The most-likely gradient fluctuation that induces caustics at small St, by contrast, is the same in the non-Gaussian and Gaussian models. Direct-numerical simulation results for particles in turbulence show that the optimal fluctuation is similar, but not identical, to that obtained by the model calculations.


I. INTRODUCTION
Particle inertia causes heavy particles suspended in an incompressible turbulent flow to cluster and form fractal spatial patterns [1], because it allows the particles to detach from the flow.This happens at caustics, singular points of the particle dynamics, where the volume of the infinitesimal particle neighbourhood collapses to zero [2][3][4][5][6].Between caustics, the particle velocities are multi-valued, giving rise to anomalously large relative particle-velocities at small separations [3,[7][8][9].This effect is also known as the "sling effect" [3,10,11].It causes spatial continuum descriptions of the inertial particle velocity field to fail.To some degree, the nonsmooth contribution from multivaluedness can be modeled using spatially uncorrelated fluctuations [12], or more accurate approaches [13][14][15].
Theoretical analysis of Gaussian statistical models resolving individual particle trajectories shows that, at small Stokes numbers, caustics form by an optimal fluctuation that involves large fluid strain and zero fluid vorticity [16].In the plane spanned by the invariants Q = − 1  2 TrA 2 and R = −detA of the traceless matrix A of fluid-velocity gradients, the optimal fluctuation that induces a caustic follows the right branch of the Vieillefosse line 27  4 R 2 + Q 3 = 0 [17].However, turbulent fluid-velocity gradients are not Gaussian distributed, which raises the question to which extent non-Gaussian fluctuations change these results.
To address this question, we formulate a statistical model that exhibits non-Gaussian fluid-velocity statistics and analyse caustic formation in the small-St limit [16][17][18][19][20] using optimal-fluctuation theory.We find that the optimal fluctuation agrees with that in the Gaussian models at small St: caustics form along the right branch of the Vieillefosse line [Fig.1(a)].The rate of caustic formation J τ K , by contrast, depends strongly on the shape of the non-Gaussian tails of the distribution of A. At small Stokes numbers, its St-dependence is determined by a non-universal action that is sensitive to the tails of the distribution of fluid-velocity gradients [Fig.1(b)].As St → 0, the rate of caustic formation obeys an asymptotic law in the form in d spatial dimensions.Here, F is a scaling function that we determine explicitly and M ≥ 1 is a model parameter described in Section III. Figure 1(c) shows the collapse of the data onto the scaling function F (a) for small enough St.For St √ M ≫ 1, the non-Gaussian model coincides with the Gaussian model and we find inversely quadratic scaling in St, −ln(J τ K ) ∝ St −2 , in agreement with earlier models [3,18,21].For St √ M ≪ 1, by contrast, the St scaling is inversely linear, −ln(J τ K ) ∝ St −1 .The different St-scalings arise because the rate of caustic formation is determined by the tails of the probability distribution of the fluid-velocity gradients, which depend on the relative magnitudes of M and St.The incomplete collapse of the data in Fig. 1(c) is a finite-St effect discussed in Section VI, as indicated by the correction term O(St −1/3 ) in Eq. (1).Including the leading-order correction in St (Section VI) improves the collapse, see inset of Fig. 1(c).
In contrast to the sensitive dependence of J on the parameters, the most likely history of fluid-velocity gradients that induces a caustic -the optimal fluctuation -is the same for all parameters M of the non-Gaussian model, and identical to that of the Gaussian model [16,17] for small St.This indicates that the optimal fluctuation that induces a caustic is more robust than the rate of its occurrence, which is strongly model dependent.To emphasise this point, we compare our results with direct numerical simulations (DNS) of particles in turbulence and find qualitative agreement for the optimal fluctuation, providing further evidence for its robustness.
Our results explain why studies aimed at comparing the rate of caustic formation J from DNS with that from idealised (and often Gaussian) models fail: For small St, J depends of the tails of the distribution of fluid-velocity gradients which is < l a t e x i t s h a 1 _ b a s e 6 4 = " x j P j 8 l / G 9 L N 0 w 8 p 5 T U A j j 5 x d 6 J a / l < / l a t e x i t > ⇠ 0.6a < l a t e x i t s h a 1 _ b a s e 6 4 = " V A 6 W C W W A A w / L L Y C n q t z k k i q C R f 8 = " > A A A C 0 H i c b V H B T t w w E P W m t K W h L d A e e 7 G 6 Q u p h t U o q B D 3 0 s F U v X J A o s I C 0 S Z H t T H a t O H F k T 1 B X 0 a r q l U O v 7 Z f w L / x N 7 W W L W G A k y 0 9 v 5 n m e Z 3 i t p M U o u u 4 E T 1 a e P n u + + i J c e / n q 9 f r G 5 p s T q x s j Y C i 0 0 u a M M w t K  < l a t e x i t s h a 1 _ b a s e 6 4 = " r 8 K Z E s d L A 3 / D k k U j G 8 h B S E n y r 5 g = " > A A A C 0 H i c b V H B T t w w E P W m t K W h L d A e e 7 G 6 Q u p h t U o q B D 3 0 s F U v X J A o s I C 0 S Z H t T H a t O H F k T 1 B X 0 a r q l U O v 7 Z f w L / x N 7 W W L W G A k y 0 9 v 5 n m e Z 3 i t p M U o u u 4 E T 1 a e P n u + + i J c e / n q 9 f r G 5 p s T q x s j Y C i 0 0 u a M M w t K      19) (thick black line), using the St → 0 limit of the scaling variable a in Eq. ( 1).The inset shows the same, but using the finite Stokes correction for a in Eq. ( 21) with δ λ th from the theory in Table I.
non-Gaussian and unknown in general.

II. CAUSTIC FORMATION
In the Stokes approximation, the equation of motion for a small, yet heavy spherical particle reads [6] Here x x x and v v v are particle position and velocity, u u u(x x x,t) is the turbulent velocity at the particle position at time t, and dots denote time derivatives.Gravitational settling is neglected (the effect of settling on caustic formation is discussed in Refs.[22,23]).The Stokes approximation applies for small enough particles that are much heavier than the fluid.In this limit, particle inertia can nevertheless be important, even when the particles are small.Caustics occur when particle-velocity field folds over configuration space, so that the particle-velocity gradients ∂ v i /∂ x j tend to −∞ [5,16,17].To identify caustics, one therefore follows the particle-velocity gradients and determines when they diverge.To this end, one integrates [3] alongside Eq. ( 2).The matrices A and Z contain the fluid-velocity gradients A i j = τ p ∂ u i /∂ x j and particle-velocity gradients Z i j = τ p ∂ v i /∂ x j at the particle position, non-dimensionalised with the particle-response time τ p .Equations ( 2) and ( 3) are used to identify caustic locations.Approaches that aim to resolve the particle number density close to a caustic, by contrast, require also higher-order derivatives of the particle velocity [15,24].We consider caustic formation at weak particle inertia, St ≪ 1.In this case, the dynamics of Z occurs on a time scale of order τ p , much shorter than τ K , the time scale of changes of A. In this "persistent limit" [16][17][18][19][20], changes of A can be treated as approximately constant.In this approximation, caustics occur whenever Eq. ( 3) has no stable fixed points, because then Z escapes to negative infinity in finite time.In three spatial dimensions, the fixed points of Eq. ( 3) vanish whenever the fluidvelocity gradients A exceed a threshold in the Q-R plane, parameterised by the line [17] In the chosen non-dimensionalisation, the elements of A are typically of the order of St, and thus small when St ≪ 1.Since the threshold (4) is of order unity, this means that rare large fluctuations of A are needed to reach the threshold line.
An equivalent condition can be formulated in terms of the eigenvalues of A [16,17,25].It requires that the most negative eigenvalue of A be real and drop below a negative threshold that tends to −1/4 as St → 0. In contrast to Eq. ( 4), this condition is independent of the spatial dimension d.
Inertialess tracers move along the stream lines of the flow and sample configuration space homogeneously [26].As a consequence, the Lagrangian statistics of A agrees with the Eulerian one.Inertial particles, by contrast, preferentially sample straining regions [27] and distribute inhomogeneously over configuration space.This leads to different magnitudes of the Lagrangian correlation functions of of strain S = (A + A T )/2 and vorticity O = (A − A T )/2 for finite St [28][29][30][31], Here, f S (x) and f O (x) are normalised correlation functions with f S (0) = f O (0) = 1 and C S (St) and C O (St) quantify the relative magnitude of strain and vorticity along inertial particle trajectories.Preferential concentration [32,33] leads to C S (St) > C O (St), i.e., strain dominates at finite St [27].However, preferential concentration is absent for very small St, so that For Gaussian-distributed fluid-velocity gradients, Eqs.(5) imply that the steady-state probability P s (A) for the occurrence of a realisation A of fluid-velocity gradients is given by where we dropped the normalisation factor because it is subleading in the limit of small St.To leading order in St, the logarithmic probability for A to reach the threshold A th determines the rate of caustic formation, −ln(J τ K ) ∼ −lnP(A th ).For Gaussian fluid-velocity gradients (6), this leads to the inverse quadratic scaling in St found in Refs.[16,17], −ln(J τ K ) ∝ St −2 as St → 0.

III. NON-GAUSSIAN MODEL
Motivated by the fact that turbulent fluid-velocity gradients have non-Gaussian tails, we explore the effects of such tails on caustic formation.To this end, we formulate a non-Gaussian model for the turbulent fluid-velocity gradients A driving Eq. ( 3).It is based on an ensemble of M identical, time-independent random velocity fields u u u m (x x x), all with identical smooth spatial correlation functions.We superpose these u u u m (x x x) with random, time-dependent coefficients c m (t) to obtain a non-Gaussian model for turbulent fluctuations, Here c m (t) are independent, identical Gaussian processes with zero mean and covariance where f (0) = 1 and with the correlation time τ c of u u u(x x x,t) defined as τ c = ∫ ∞ 0 dt f (t).For √ M ≫ |A i j |, the central-limit theorem implies that A is Gaussian distributed and P s (A) is given by Eq. ( 6), so that the results of [16,17] are recovered.For large fluid-velocity gradients and finite M, however, the central-limit theorem does not apply.In this case, the distribution of A has non-Gaussian tails, whose form we discuss in the next section.

IV. OPTIMAL-FLUCTUATION THEORY
For finite M and large fluid-velocity gradients, the central limit theorem cannot be not applied to Eq. (7).In this case, we find by a saddle-point analysis of the non-Gaussian model The derivation of Eq. ( 9) is summarised in Appendix A. We see that the argument of the square root in Eq. ( 9) is just 2M times the Gaussian action on the right-hand side of Eq. ( 6), which simplifies the analysis that follows.
As mentioned earlier, at small St ≪ 1, caustics form whenever the most negative eigenvalue of A is real and drops below a negative threshold that tends to −1/4 for St → 0 [16,17,25].In the persistent limit, the non-dimensional rate of caustic formation J τ K is equal to the probability for this event, P s (A th ), to leading order in St.In order to compute this probability, we first express the right-hand side of Eqs. ( 6) and ( 9) in terms of the eigenvalues λ λ λ = (λ 1 ,...,λ d ) T of S and in terms of the non-zero singular values µ µ µ = (µ 1 ,..., where Transforming the joint probability density of A into the probability of λ λ λ and µ µ µ in this way involves a Jacobian determinant.This so-called Vandermode determinant, however, is independent of St and does therefore not contribute at this order.We marginalise the distribution (10) over λ 1 ,...,λ d−1 and µ 1 ,..., µ ⌊d/2⌋ using the contraction principle [34] and compute by constrained minimisation.The constraint ∑ d n=1 λ n = 0 arises because A must be traceless in incompressible flow.We enforce the constraint by adding a Lagrange multiplier.Evaluating the stationary point of the resulting Lagrangian gives λ * k = −λ d /(d −1) and µ * j = 0, independently of H .This implies that the optimal gradient configuration that leads to a caustic has vanishing vorticity, O * ≡ 0, as shown in Refs.[16,17], while one of the eigenvalues of S * , λ d here, becomes large and negative.All other eigenvalues λ k are most likely equal and positive, and of 1/(d − 1) times the magnitude of λ d , in keeping with the incompressibility constraint.Hence after diagonalisation and ordering of the eigenvalues, the most likely threshold configuration to induce a caustic has the form where λ th is the threshold value for the most-negative eigenvalue of A * th and e is the diagonal matrix For d = 3, the fact that all subleading eigenvalues λ * k are identical, λ * k = λ th /(d − 1), implies that caustics are most likely induced by fluid-velocity gradients that lie on the Vieillefosse line, not only in the Gaussian [17] but also in the non-Gaussian model.Figure 1(a) shows the excursion of fluid-velocity gradients in the Q-R plane at a time of order τ K before caustic formation, obtained from simulations of the non-Gaussian model with M = 1.We see that the fluid-velocity gradients prior to caustic formation lie beyond the threshold (solid line) and close to the Vieillefosse line (thick dashed line), as predicted by our analysis.Note that most gradients in Fig. 1(a) lie considerably below the threshold line.The reason is that at non-zero St [= 0.1 in Fig. 1(a)], the threshold λ th is larger than for St → 0, because the strict separation of time scales breaks down and caustic formation takes a finite time for St > 0. This means that the smallest eigenvalue of the matrix of fluid-velocity gradients must fall below −1/4 for a finite time in order to enable caustic formation [16,17,25].Consequently, at finite St, the threshold value λ th in Eq. ( 12) must be adjusted, i.e., λ th = 1/4 + δ λ th , where δ λ th is a correction that vanishes as St → 0. From the numerical solution of a model system, we show in Section V that δ λ th scales as δ λ th ∼ St 2/3 for St ≪ 1.

V. SHAPE OF OPTIMAL FLUCTUATION
We now consider the shape of the optimal fluctuation A * (t) of fluid-velocity gradients as a function of time that causes the most negative eigenvalue to drop below −λ th at time t th .To this end, we note that the matrix e in Eq. ( 12) can be chosen as an element of an orthonormal basis of traceless matrices with respect to the inner product ⟨X,Y⟩ = d −1 (d − 1)Tr(X T Y) defined by the trace Tr of the product of two matrices X and Y.The normalisation is chosen here such that ⟨e,e⟩ = 1.In such a basis, the most likely way to reach the threshold A * th is by a large excursion of the amplitude A * (t) associated with the basis element e, while the amplitudes of all other other basis elements remain small.In this way, the amplitude A * (t) alone determines the shape of the optimal fluctuation A * (t) for small St, reflecting the effecitively one-dimensional nature of caustic formation at small St [16][17][18]25].The eigenvalue threshold −λ th in Eq. ( 12) is then reached when A * (t) reaches A * (t th ) = −λ th at time t th .
In the previous Section, we discussed that the breakdown of time-scale separation increases the threshold value λ th (St) for finite St. Using the one-dimensional parametrisation (14), we now make this statement more precise by considering the combined dynamics of the optimal fluctuations of fluid-velocity gradients and particle-velocity gradients.To this end, we model the optimal fluctuation A * (t) as the optimal fluctuation of an Ornstein-Uhlenbeck process.For this particular Gaussian process, the joint optimal fluctuations of the fluid-velocity gradient A * (t) and the escaping part of the particle-velocity gradient Z * (t) obey a < l a t e x i t s h a 1 _ b a s e 6 4 = " i c E x t e 5 3 z I u c L 3 i M + 4 0 t q / T A k  closed set of differential equations, see e.g.[21] for a derivation.These equations read where 0 ≤ t ≤ t c and p z and p A are conjugate momenta that enable non-trivial optimal fluctuations A * and Z * .The parameter and the non-dimensionalised Lagrangian strain correlation time κ S = τ −1 K ∫ ∞ 0 dt f S (t) reflect the properties of fluid-velocity gradients along Lagrangian trajectories [28,31] in dimensions d > 1.In order to generate a caustic, Z * in Eqs.(15) must escape from the trivial fixed point (Z * ,A * , p z , p A ) = (0,0,0,0) at t = 0 and arrive at Z * → −∞ at time t = t c , the time at which a caustic forms.The condition that the fluctuation be optimal leads to the additional boundary conditions p z (t c ) = p A (t c ) = 0.
To find the optimal fluctuation that satisfies Eqs.(15) with the prescribed boundary conditions, we employ a numerical shooting method.Taking initial conditions on the unstable subspace of the trivial fixed point, we "shoot" trajectories obeying Eqs.(15), to "hit" the target boundary condition at t = t th .Since the trivial fixed point has two unstable directions, we can parametrise the initial conditions by a single angle θ and write (Z * ,A * , p z , p A ) = ε(cosθ e e e 1 + sinθ e e e 2 ), with ε ≪ 1 at t = 0. Here, the vectors e e e 1 and e e e 2 furnish an orthonormal basis of the unstable subspace of the trivial fixed point.By shooting trajectories that obey Eqs.(15) with initial conditions parametrised by θ , we optimise θ to hit the target boundary conditions Z * (t c ) → −∞, and p z (t c ) = p A (t c ) = 0.The trajectory A * (t) obtained in this way gives the amplitude A * (t) of the optimal fluctuation A * (t) in Eq. ( 14) for small but finite St when the fluid-velocity gradients are Ornstein-Uhlenbeck processes.
In order to understand how the optimal fluctuation obtained from the shooting method compares with the optimal fluctuation of the non-Gaussian model, we first consider a simplified, one-dimensional model with σ 2 S = κ S = 1. Figure 2(a) shows the resulting amplitude A * (t) for different St plotted vs. t − t th .We observe that the optimal fluctuation approaches the shape of the correlation function of the Ornstein-Uhlenbeck process with κ S = 1 as St → 0 [Fig.2(a)].This has been shown to be the case for all Gaussian processes in this limit [16,17].In Appendix B we show that this also holds for the non-Gaussian process studied here, which indicates that the optimal fluctuation of the non-Gaussian model approaches the correlation function in a similar way.
To verify that the optimal fluctuations of the non-Gaussian model for small St follow the same trend as observed in Fig. 2(a) for the Ornstein-Uhlenbeck process, we compare A * (t) obtained from the shooting method with results of numerical simulations of the non-Gaussian model.To simplify the simulations and to thus allow for smaller St, we consider again the one-dimensional version of the non-Gaussian model, obtained from by Eqs. ( 7) and (3) with d = 1.In the one-dimensional model, the spatial dependence of the fluid-velocity gradients is neglected and the c m (t) are taken as independent Ornstein-Uhlenbeck processes with σ 2 S = κ S = 1.To obtain the optimal fluctuation, we integrate Eqs. ( 7) and (3) until a caustic forms and backtrace the history of the fluid-velocity gradient A(t) prior to caustic formation at t = t c .From the resulting gradient trajectories we compute for M = 1 the mean, shown in Fig. 2(b) for different values of St, and the trajectory density, shown in Fig. 2(c) for one small value of St.We observe that the trajectory density in Fig. 2(c) sharply focuses in a narrow region of high density.This is the optimal fluctuation, which is well approximated by the mean value of the fluid-gradient at each time step (red line).The mean values in Fig. 2(b) show a similar trend as the Ornstein-Uhlenbeck process in Fig. 2(a), approaching the shape of the correlation function as St becomes smaller.Moreover, Fig. 2(b) shows that for decreasing St the threshold time t th where A * (t) has its minimum approaches the caustic formation time t c .This implies that after a caustic has been initiated at t = t th , it takes shorter and shorter time for it to form at t = t c as St decreases.Thus, as St → 0 caustics form instantaneously whenever the amplitude A * (t) reaches −1/4 [16,17].Note also that the smaller St, the larger the magnitude of A * (t) at t = t c .Since the optimal fluctuation A * (t) is a pure strain in our model, this implies that caustic formation events at small St are strongly correlated with straining regions of the flow.This provides an explanation for increased particle collisions in straining regions observed in Refs.[35][36][37].
The mean value for small St in Fig. 2(c) agrees well with the optimal fluctuation obtained from the shooting method (dashed line), even though M = 1, while the shooting method applies strictly only in the limit M → ∞.This suggests that the optimal fluctuation is insensitive to changes of the model when St is small.Further evidence for the robustness of the optimal fluctuation is presented in Fig. 2(d).The figure shows δ λ th , the magnitude by which the optimal fluctuation falls below the value −1/4 as a function of St, obtained from numerical simulations of the one-dimensional model and from the shooting method.We observe a mild dependence of δ λ th on M at finite St, where M = ∞ corresponds to the Gaussian model.In the limit St → 0 all threshold values approach −1/4 with scaling δ λ th ∝ St 2/3 extracted numerically, albeit with M-dependent prefactors.The same St 2/3 scaling can also be obtained from a bound derived in Bätge et al. [25], by using the optimal fluctuation to establish a relation between the duration of exceeding the threshold and δ λ th , D and ∆ in their Eq.( 8), respectively.We conclude that the optimal fluctuation of fluid-velocity gradients is generally robust against changes of the model.The rate of caustic formation, by contrast, is sensitive to the tails of the distribution, as shown in Fig. 1(b).
The observed robustness of the optimal fluctuation with respect to changes of M justifies using A * (t) obtained from the shooting method for finite St as the amplitude for the d-dimensional optimal fluctuation in Eq. ( 14).With the appropriate parameters σ 2 S = 2C S d −1 (d + 2) −1 ∼ 1/15 and κ S obtained from simulations of the non-Gaussian model in three spatial dimensions, d = 3, we numerically extract the small-St scaling of δ λ th from the shooting method.The results are given in Table I.We observe that κ S depends on M, which introduces a weak M-dependence of the scaling prefactor of δ λ th extracted from the shooting method.

VI. DATA COLLAPSE ONTO SCALING FUNCTION
With the one-dimensional picture (14) of caustic formation developed in the previous section, we now derive the scaling function [black line in Fig. 1(c)] onto which the caustic formation data collapses for small St.To this end, we again diagonalise A(t) at the time of caustic formation.We then obtain from Eq. ( 7) the following expression for the component of the fluid velocity gradient that reaches the threshold A th : where A m = ⟨A m ,e⟩ is the amplitude of A m along the basis matrix e in Eq. ( 13).Due to homogeneity and stationarity, the steadystate distribution P s (A) of A, is independent of x x x and t, for any single point in space and time.Furthermore, since c m and A m are independent, A follows a large-deviation principle [34] for M ≫ 1, with rate function F (a).For Gaussian c m and A m , the function F (a) can be obtained explicitly by first computing the generating function of A and applying a Legendre transform [34].We find data for turbulent fluid-velocity gradients scatters broadly along and slightly below the Vieillefosse line.For the non-Gaussian model [Fig.1(a)], by contrast, the scatter appears more centred on the Vieillefosse line.
The elongated scatter of turbulent fluid-velocity gradients is likely due the fact that the non-linear terms in the Navier-Stokes equations self amplify fluid-velocity gradients along this branch [41,42].This generates large gradient excursions and thus leads to anomalously large fluctuations along the right branch of the Vieillefosse line.That the turbulent gradients scatter below the Vieillefosse line, instead of centering on the line as in the non-Gaussian model, can be understood by decomposing Q = −trS 2 /2 + trOO T /2 and R = −trS 3 /3 − trSO 2 .If O = 0, no pair of (R,Q) may lie above the Vieillefosse line by its definition.In this case, the distribution is concentrated just below the Vieillefosse line, for both the Gaussian model and the DNS (not shown).The effect of non-zero O is to first shift Q upwards, because trOO T ≥ 0. Second, a negative value of trSO 2 shifts R to the right, and vice versa.This explains why the DNS data lies below the Vieillefosse line in the Q − R plane.The reason is that trSO 2 is unlikely to take negative values.This is apparent from Figure 3(b), which shows that the probability to have negative trSO 2  is very small at the onset of caustic formation (the onset time is indicated by the vertical dashed line).This peculiar coupling between strain and vorticity as the strain grows large is specific to turbulence.In the statistical model, the distribution of trSO 2 is almost symmetric, see Fig. 3(c), explaining why the data in Fig. 1(a) is scattered both above and below the Vieillefosse line.
The strong interrelation between strain and vorticity in turbulence discussed above also results in vorticity giving a finite, correlated contribution to the optimal fluctuation of Q (not shown) [43].This is in contrast to the Gaussian and non-Gaussian models considered here, where the optimal vorticity contribution vanishes [16,17].Despite these differences, the shape of the optimal fluctuation in the DNS is still quite well approximated by the correlation function, as shown in Fig. 3(d), where the trajectory density of the eigenvalues of the turbulent strain S, are depicted prior to caustic formation.The optimal fluctuation theory described above assumes that the mean values ⟨λ i ⟩ vanish.This is consistent, because they are of higher order than δ λ ∼ St 2/3 for small St. Since the Stokes number in Fig. 3(d) is not very small, we subtracted the mean values from λ i .The solid blue lines show the theoretical predictions for the optimal fluctuations, given by the correlation function f S of S normalised to a fitted threshold value.We observe qualitative agreement.
Finally, Fig. 3(e) shows the St-dependence of rate of caustic formation J τ K for particles in turbulence from Refs.[10,38], using DNS with different Reynolds numbers Re λ .In the literature, different authors have fitted different exponents, ranging between minus one and minus two, to the St-dependence of −lnJ τ K [10,38], motivated by results for the rate of caustic formation in Gaussian statistical models [5].The data in Fig. 3(e) appears to be consistent with inverse linear scaling St −1 for some of the cases, but the scaling exponent depends strongly on Re λ .This is consistent with the observations in Ref. [25] and supported by the analysis of the non-Gaussian statistical model: the small-St scaling of −lnJ τ K is strongly affected by the extreme-value statistics of fluid-velocity gradients in turbulence, not captured by Gaussian models.

VIII. CONCLUSIONS
We formulated a non-Gaussian statistical model for turbulent aerosols, with non-Gaussian tails of the distribution of fluidvelocity gradients.We analysed how caustics form in this model at small Stokes numbers, in the persistent limit.Using large-deviation theory we computed the rate of caustic formation for the non-Gaussian model and demonstrated that caustic formation is sensitive to the tails of the distribution of fluid-velocity gradients, because they must overcome a threshold for a caustic to form: The smaller the Stokes number, the larger the fluid-velocity gradients must be to initiate a caustic.For St ⪅ 0.1 the required fluid-velocity gradients lie far in the tails of the distribution.
While the rate of caustic formation depends sensitively on the tails of the fluid-velocity gradient distribution seen by the particles, the form of the optimal fluctuation is robust against changes of the parameters in the non-Gaussian model.In particular, the optimal fluctuation of the non-Gaussian model has small vorticity, just like the Gaussian one, and it follows the right branch of the Vieillefosse line at small Stokes numbers.
What are the implications of these results for caustics in three-dimensional homogeneous isotropic turbulent flow?First, our analysis of the non-Gaussian model shows that the rate of caustic formation is sensitive to the tails of the fluid-velocity gradient distribution.As a consequence, the rate of caustic formation is not expected to have a simple scaling form such as St −1 or St −2 .This explains why numerical studies of the rate of caustic formation using direct numerical simulation of turbulence at different Reynolds numbers yielded inconclusive results [10,38], simply because the St-dependence of −ln(J τ K ) is not a power law.These conclusions are consistent with the results of Ref. [25], who analysed the Reynolds-number dependence of the rate of caustic formation using DNS of turbulence.Second, the optimal fluctuation of fluid-velocity gradients leading to caustic formation in turbulence is similar to the predictions of our non-Gaussian model.In particular, the fluid-velocity gradients that initiate caustics in turbulence lie close to the right branch of the Vieillefosse line.This robustness explains why particle collisions in turbulence tend to occur close to the Vieillefosse line [37].Third, the vorticity contribution to the optimal fluctuation in turbulence, however, is not simply scattered around zero but performs a small, yet correlated excursion.This suggests an interdependence of vorticity and strain in the tails of the joint distribution of turbulent strain and vorticity [25].We remark that one can produce an optimal fluctuation with non-zero vorticity by adding third-and fourth-order terms to the (quadratic) Gaussian action in Eq. ( 6).Whether or not the finite optimal vorticity in turbulence has a similar origin is an open question.
< l a t e x i t s h a 1 _ b a s e 6 4 = " O h e q l P N k R m T B c c D / z + x z s B x H d 3 8 = " > A A A C Q H i c b V H L S s N A F J 3 4 r P X V 6 t L N a B E q S E m k q M u C G 5 c V + p I 2 l M l k 0 g 6 Z S c L M p F B C v s K t / o x / 4 R + 4 E 7 e u n L R Z m N Y L A 4 d z X + e e c S J G p T L N D 2 N j c 2 t 7 Z 7 e 0 V 9 4 / O D w 6 r l R P e j K M B S Z d H L J Q D B w k C a M B 6 S q q G B l E g i D u M N J 3 / I c s 3 5 8 R I W k Y d N Q 8 I j Z H k 4 B 6 F C O l q e d 6 M n I 8 i N K r c a V m N s x F w H V g 5 a A G 8 m i P q 8 b 5 yA 1 x z E m g M E N S D i 0 z U n a C h K K Y k b Q 8 i i W J E P b R h A w 1 D B A n 0 k 4 W i l N 4 q R k X e q H Q L 1 B w w f 7 t S B C X c s 4 d X c m R m s r V X E b + l x v G y r u 3 E x p E s S I B X i 7 y Y g Z V C L P z o U s F w Y r N N U B Y U K 0 V 4 i k S C C t t U m F L N l t I T x Y v 6 V h 2 k k n O h h f K / S i j 5 b X + D J + I G W J u W t a m W q s Wr o P e T c O 6 b T S f m r W W m d t b A m f g A t S B B e 5 A C z y C N u g C D D h 4 A a / g z X g 3 P o 0 v 4 3 t Z u m H k P a e g E M b P L 3 J q r + E = < / l a t e x i t > (a) < l a t e x i t s h a 1 _ b a s e 6 4 = " U e S 3 5 K x T TG O F O Z f u Y U r f d 3 O 6 f K Y = " > A A A C O H i c b V H L S s N A F J 3 4 r P X V 6 t L N a B F c S E m k q C s p u H H Z S l / Q h j K Z 3 L R D J w 9 m J o U S + g V u 9 W f 8 E 3 f u x K 1 f 4 K T N w r R e G D i c + z r 3 j B N x J p V p f h g b m 1 v b O 7 u F v e L + w e H R c a l 8 0 p F h L C i 0 a c h D 0 X O I B M 4 C a C u m O P Q i A c R 3 O H S d y W O a 7 0 5 B S B Y G L T W L w P b J K G A e o 0 R p q v k 8 L F X M q r k I v A 6 s D F R Q F o 1 h 2 T g f u C G N f Q g U 5 U T K v m V G y k 6 I U I x y m B c H s Y S I 0 A k Z Q V / D g P g g 7 W S h d I 4 v N e N i L x T 6 B Q o v 2 L 8 d C f G l n P m O r v S J G s v V X E r + l + v H y r u 3 E x Z E s Y K A L h d 5 M c c q x O n Z 2 G U C q O I z D Q g V T G v F d E w E o U q b k 9 u S zh b S k / l L W p a d p J L T 4 b n y S Z T S 8 l p / w g T E l H B 3 X t S m W q s W r o P O T d W 6 r d a a t U r 9 I b O 3 g M 7 Q B b p C F r p D d f S E G q i N K A L 0 g l 7 R m / F u f B p f x v e y d M P I e k 5 R L o y f X 3 t Y r Q M = < / l a t e x i t > R < l a t e x i t s h a 1 _ b a s e 6 4 = " e 9 m 4 Z 1 2 2 g + N W s t M 7 e 3 B M 7 A B a g D C 9 y B F n g E b d A F G H D w A l 7 B m / F u f B p f x v e y d M P I e 0 5 B I Y y f X 3 R A r + I = < / l a t e x i t > (b) < l a t e x i t s h a 1 _ b a s e 6 4 = " / u W e C I 4 u h d L k I c / x o u f G c E 3 P Y k s 4 X 0 Z P G S j m U n m e R s e K H c j z J a X u m / 8 I m Y I e a m F W 2 q t W r h O u h d N 6 y b R v O p W W / d 5 / a W w S k 4 B 5 f A A r e g B R 5 B G 3 Q B B g y 8 g F f w Z r w b n 8 a X 8 b 0 s L R l 5 z w k o h P H z C 6 D V s B E = < / l a t e x i t > St < l a t e x i t s h a 1 _ b a s e 6 4 = " / / H u j 8 A K e h g x + Q p 8 X 9 u / j 5 b X + D J + I G W J u W t a m W q s W r o P e T c O 6 b T S f m r W W m d t b A m f g A t S B B e 5 A C z y C N u g C D D h 4 A a / g z X g 3 P o 0 v 4 3 t Z u m H k P a e g E M b P L 3 Y W r + M = < / l a t e x i t > (c) < l a t e x i t s h a 1 _ b a s e 6 4 = " C Z b / M W H x F B B + O U g m S 9 8 o 6 m I G N s U 7 G a y 8 j / c s N E e n d 2 S s M 4 k S T E y 0 V e w q C M Y O Y A d C k n W L K Z A g h z q r R C P E E c Y a l 8 K m z J Z n P h i e I l H d N O M 8 n Z 8 E K 5 H 2 e 0 u F L / 4 R M + R c y d V 5 S p 5 q q F 6 6 D X 0 M 0 b v f n U r L f u c 3 v L 4 B S c g 0 t g g l v Q A o + g D b o A g x i 8 g F f w p r 1 r n 9 q X 9 r 0 s L W l 5 z w k o h P b z C 8 y M s I Y = < / l a t e x i t > ⇠ 0.3a 2 < l a t e x i t s h a 1 _ b a s e 6 4 = " j 8 S y t w 4 9 T H m u F 1 e l w e P u V r O S H Y w 6 2 6 D 4 z h t v X n / + l K m L W r P 2 5 5 b S Q H m g q l s 9 u i D v f + f d D e 3 U 1 d U s g K Y W x 8 6 E 2 G S Q Z 7 4 F G + 9 e u B l Y 6 O b m i b 7 g y F N n P T L / h G / l W m c w N I Q R 5 y P D U C h 0 n l L z r W a h a H b d H x / r w / B y c d + v N P f / r b d H X x e 7 H y V v C P v y Q c S k 1 0 y I H v k g A y J I B P y m / w h f 4 P D 4 E f w M / h 1 U x p 0 F p q 3 Z C m C y 3 8 2 S e U n < / l a t e x i t > 10 0 < l a t e x i t s h a 1 _ b a s e 6 4 = " r 8 K Z E s d L A 3 / D k k U j G 8 h B S E n y r 5 g = " > A A A C 0 H i c b V H B T t w w E P W m t K W h L d A e e 7 G 6 Q u p h t U o q B D 3 0 s F U v X J A o s I C 0 S Z H t T H a t O H F k T 1 B X 0 a r q l U O v 7 Z f w L / x N 7 W W L W G A k y 0 9 v 5 n m e Z 3 i t p M U o u u 4 E T 1 a e P n u + + i J c e / n q 9 f r G 5 p s T q x s j Y w 6 2 6 D 4 z h t v X n / + l K m L W r P 2 5 5 b S Q H m g q l s 9 u i D v f + f d D e 3 U 1 d U s g K Y W x 8 6 E 2 G S Q Z 7 4 F G + 9 e u B l Y 6 O b m i b 7 g y F N n P T L / h G / l W m c w N I Q R 5 y P D U C h 0 n l L z r W a h a H b d H x / r w / B y c d + v N P f / r b d H X x e 7 H y V v C P v y Q c S k 1 0 y I H v k g A y J I B P y m / w h f 4 P D 4 E f w M / h 1 U x p 0 F p q 3 Z C m C y 3 8 4 s u U o < / l a t e x i t > 10 1

1 <
8 f 6 8 P w e m H Y f x x u P 9 t v z / 6 v N z 5 B t p F b 9 E 7 F K N P a I S + o m O U I I Z K d I W u 0 e 8 g C W z w M / h 1 U x r 0 l p o 3 a C W C y 3 9 w u O Z r < / l a t e x i t > 10 l a t e x i t s h a 1 _ b a s e 6 4 = " L A e D w m r N j O H F z H f m Y 4 S 2 P s 6 m A d

1 <
8 f 6 8 P w e m H Y f x x u P 9 t v z / 6 v N z 5 B t p F b 9 E 7 F K N P a I S + o m O U I I Z K d I W u 0 e 8 g C W z w M / h 1 U x r 0 l p o 3 a C W C y 3 9 w u O Z r < / l a t e x i t > 10 l a t e x i t s h a 1 _ b a s e 6 4 = " z s 7 9 A i d t F q b 1 w s D h 3 N e 5 Z 5 y I M 6 l M 8 w M X t r Z 3 d v e K + 6 W D w 6 P j k 3 L l t C v D W F D o 0 J C H o u 8 Q C Z w F 0 F F M c e h H A o j v c O g 5 0 6 c 0 3 5 u B k C w M 2 m o e g e 2 T c c A 8 R o n S Z H r u Z T 8 L z e I l f d g J y y I Y g U B X S 3 y Y m 6 o 0 E j P N l w m g C o + 1 4 B Q w b R W g 0 6 I I F R p c 3 J b 0 t l C e j J / S d u y k 1 R y O j x X P o 1 S W t 7 o T 5 i C m B H u L k r a V G v d w k 3 Q v a 1 Z d 7 V 6 q 1 5 t P G b 2 F t E 5 u k T X y E L 3 q I G e U R N 1 E E W A X t A r e s P v + B N / 4 e 9 V a Q F n P W c o F / j n F 5 b E r R I = < / l a t e x i t > a < l a t e x i t s h a 1 _ b a s e 6 4 = " U 6 2 h Y c h N 9 k c u g S k G U i I 9 c k b 3 t K FIG. 1. Caustic formation in a non-Gaussian model for a turbulent aerosol.(a) Probability of Q and R at the onset of caustic formation (colour coded).Parameters: St = 0.1 and M = 1.The thick dashed line is the right branch of the Vieillefosse line (see text).The solid line is the threshold line for caustic formation, Eq. (4).(b) Rate of caustic formation for the non-Gaussian model for different values of M. Also shown are the limiting St scalings, dashed lines.(c) Collapse of the small-St data from panel (b) (markers) onto the scaling function F (a) given in Eq. (19) (thick black line), using the St → 0 limit of the scaling variable a in Eq. (1).The inset shows the same, but using the finite Stokes correction for a in Eq. (21) with δ λ th from the theory in TableI.

<
l a t e x i t s h a 1 _ b a s e 6 4 = " O h e q l P Nk R m T B c c D / z + x z s B x H d 3 8 = " > A A A C Q H i c b V H L S s N A F J 3 4 r P X V 6 t L N a B E q S E m k q M u C G 5 c V + p I 2 l M l k 0 g 6 Z S c L M p F B C v s K t /o x / 4 R + 4 E 7 e u n L R Z m N Y L A 4 d z X + e e c S J G p T L N D 2 N j c 2 t 7 Z 7 e 0 V 9 4 / O D w 6 r l R P e j K M B S Z d H L J Q D B w k C a M B 6 S q q G B l E g i D u M N J 3 / I c s 3 5 8 R I W k Y d N Q 8 I j Z H k 4 B 6 F C O l q e d 6 M n I 8 i N K r c a V m N s x F w H V g 5 a A G 8 m i P q 8 b 5 y j 5 b X + D J + I G W J u W t a m W q s W r o P e T c O 6 b T S f m r W W m d t b A m f g A t S B B e 5 A C z y C N u g C D D h 4 A a / g z X g 3 P o 0 v 4 3 t Z u m H k P a e g E M b P L 3 J q r + E = < / l a t e x i t > (a) < l a t e x i t s h a 1 _ b a s e 6 4 = " x j P j 8 H V q H d a m w l f C + p w Q 9 1 B H r f k = " > A A A C Q H i c b V H L S s N A F J 3 4 r P X V 6 t L N a B E q S E m k q M u C G 5 c V + p I 2 l M l k 0 g 6 Z S c L M p F B C v s K t / o x / 4 R + 4 E 7 e u n L R Z m N Y L A 4 d z X + e e c S J G p T L N D 2 N j c 2 t 7 Z 7 e 0 V 9 4 / O D w 6 r l R P e j K M B S Z d H L J Q D B w k C a M B 6 S q q G B l E g i D u M N J 3 / I c s 3 5 8 R I W k Y d N Q 8 I j Z H k 4 B 6 F C O l q e d 6 M n I 8 i N O r c a V m N s x F w H V g 5 a A G 8 m i P q 8 b 5 y j 5 b X + D J + I G W J u W t a m W q s W r o P e T c O 6 b T S f m r W W m d t b A m f g A t S B B e 5 A C z y C N u g C D D h 4 A a / g z X g 3 P o 0 v 4 3 t Z u m H k P a e g E M b P L 3 Y W r + M = < / l a t e x i t > (c) < l a t e x i t s h a 1 _ b a s e 6 4 = " Y T k l m m 4 h b S k / l L W p a d p J L T 4 b n y S Z T S 8 l p / w g T E l H B 3 X t S m W q s W r o P O T d W 6 r d a a t U r 9 I b O 3 g M 7 Q B b p C F r p D d f S M G q i N K A L 0 g l 7 R m / F u f B p f x v e y d M P I e k 5 R L o y f X 1 x E r P I = < / l a t e x i t > A < l a t e x i t s h a 1 _ b a s e 6 4 = " B f T X h l 3 o x / 4 R + 4 E 7 e u n L R Z m N Y L A 4 d z X + e e c S J G p T L N D 2 N j c 2 t 7 Z 7 e 0 V 9 4 / O D w 6 r l R P e j K M B S Z d H L J Q D B w k C a M B 6 S q q G B l E g i D u M N J 3 / I c s 3 5 8 R I W k Y d N Q 8 I j Z H k 4 B 6 F C O l q e d 6 M n I 8 6 K R X 4 0 r N b e 9 m 4 Z 1 2 2 g + N W s t M 7 e 3 B M 7 A B a g D C 9 y B F n g E b d A F G H D w A l 7 B m / F u f B p f x v e y d M P I e 0 5 B I Y y f X 3 R A r + I = < / l a t e x i t > (b) < l a t e x i t s h a 1 _ b a s e 6 4 = " Y T k l m m 4 h b S k / l L W p a d p J L T 4 b n y S Z T S 8 l p / w g T E l H B 3 X t S m W q s W r o P O T d W 6 r d a a t U r 9 I b O 3 g M 7 Q B b p C F r p D d f S M G q i N K A L 0 g l 7 R m / F u f B p f x v e y d M P I e k 5 R L o y f X 1 x E r P I = < / l a t e x i t > A < l a t e x i t s h a 1 _ b a s e 6 4 = " / u W e C I 4 u h d L k I c / x o u f G c E 3 P Y z I = " > A A A C P 3 i c z / A L 3 I l b d 0 7 a L E z r h Y H D u a 9 z z z g R o 1 K Z 5 o d R 2 t j c 2 t 4 p 7 1 b 2 9 g 8 O j 6 q 1 4 5 4M Y 4 F J F 4 c s F A M H S c J o Q L q K K k Y G k S C I O 4 z 0 H f 8 h y / d n R E g a B h 0 1 j 4 j N 0 S S g H s V I a W q Q j A S H z y o d V + t m w 1 w E X A d W D u o g j / a 4 Z p y N 3 B D H n A Q K M y T l 0 D I j Z S d I K I o Z S S u j W J I I Y R 9 N y F D D A H E i 7 W Q h O I U X m n G h F w r 9 A g U X 7 N + O B H E p 5 9 z R l R y p q V z N Z e R / u W G s v D s 7 o U E U K x L g 5 S I v Z l C F M L s e u l Q Q r N h c A 4 Q F 1 V o h n i K Bs N I e F b Z k s 4 X 0 Z P G S j m U n m e R s e K H c j z J a X u m / 8 I m Y I e a m F W 2 q t W r h O u h d N 6 y b R v O p W W / d 5 / a W w S k 4 B 5 f A A r e g B R 5 B G 3 Q B B g y 8 g F f w Z r w b n 8 a X 8 b 0 s L R l 5 z w k o h P H z C 6 D V s B E = < / l a t e x i t > St < l a t e x i t s h a 1 _ b a s e 6 4 = " Y T k l m m 4 h b S k / l L W p a d p J L T 4 b n y S Z T S 8 l p / w g T E l H B 3 X t S m W q s W r o P O T d W 6 r d a a t U r 9 I b O 3 g M 7 Q B b p C F r p D d f S M G q i N K A L 0 g l 7 R m / F u f B p f x v e y d M P I e k 5 R L o y f X 1 x E r P I = < / l a t e x i t > A < l a t e x i t s h a 1 _ b a s e 6 4 = " 9 z R 1 J z a i L G E D i W o w n L X 1 a I b i S j A = " > A A A C U 3 i c b V F N S w J B G J 7 d r M y y r G 5 1 m Z L A o G w 3 p D q F 0 C X o Y q A W 6 C K z 4 6 w N O / v B z L u C L A v 9 m q 7 1 Z z r 0 W 7 o 0 q x 7 S e m H g e Z / 3 + x k 3 F l y B Z X 0 Z 5 k p h d W 2 9 u F H a 3 C p v 7 1 R 2 9 7 o q S i R l H R q J S D 6 7 R D H B Q 9 Y B D o I 9 x 5 K R w B X s y f X v 8 v j T m E n F o 7 A N k 5 g 5 A R m F 3 O O U g K Y G l Y M a n M M g 7 c s A 0 + z 0 o g 8 k m X k P 2 a B S t e r W 1 P B f Y M 9 B F c 2 t N d g 1 j v r D i C Y B C 4 E K o l T P t m J w U i K B U 8 G y U j 9 R L C b U J y P W 0 z A k A V N O O j 0 i w y e a G W I o 1 3 S t P x P 0 X z o 0 q x 7 S e m H g e Z / 3 + x k 3 E l y B Z X 0 Z Z m 5 j c 2 s 7 v 1 P Y 3 d s v H p T K h 1 0 V x p K y D g 1 F K J 9 d o p j g A e s A B 8 G e I 8 m I 7 w r W c y f 3 W b w 3 o 1 3 S t P x P 0 X z o 0 q x 7 S e m H g e Z / 3 + x k 3 E l y B Z X 0 Z Z m 5 j c 2 s 7 v 1 P Y 3 d s v H p T K h 1 0 V x p K y D g 1 F K J 9 d o p j g A e s A B 8 G e I 8 m I 7 w r W c y f 3 W b w 3 o b j a 1 m n x L m x Q U / w p K m e Y F Y S r + D Y o L x S n j t R s 8 k h o F q a U F I L S 0 u 3 I x B w 2 C r G d r U + r e 2 s R m / Z J r L y j r l e v m a + m L v K b N T / s 3 C 9 S 3 o K K q Y 0 3 1 N i 1 8 D W 6 O + 9 5 p / + T q p H d x 3 t i 7 w w 7 Z N / a D e e y M X b D f b M C G T L A / 7 I 7 d s w f n 0 X l 2 / r a a 1 J b T g H 2 2 F q 3 2 C 5 S Z t G c = < / l a t e x i t > th < l a t e x i t s h a 1 _ b a s e 6 4 = " z s FIG. 2. One-dimensional model for caustic formation.(a) Optimal fluctuations of fluid-velocity gradients, A * (t), obtained from the shooting method for the Gaussian case (see text) for different St vs. t −t th .The thick horizontal line shows the threshold −λ th = −1/4.(b) Same as (a) but for simulations of the non-Gaussian model with M = 1.The lines show the mean over realisations at each time step.(c) Optimal fluctuation as a function of t −tc using M = 1 and St = 0.04, where tc is the time at which the caustic forms.The color map shows the probability distribution of A(t), normalised at each time step by the maximal probability.The solid line is the mean obtained in panel (b) and the dashed line shows the result from the shooting method in panel (a).(d) The depth δ λ th of the optimal fluctuation from simulations of the one-dimensional models with M = 1, 2 and ∞.Also shown are the results of the shooting method (see text), as a dashed line.

TABLE I .
Numerically obtained non-dimensional Lagrangian strain correlation times κ S and scaling prefactors δ λ th St −2/3 for different parameters M in the non-Gaussian model.