Paths to caustic formation in turbulent aerosols

The dynamics of small, yet heavy, identical particles in turbulence exhibits singularities, called caustics, that lead to large fluctuations in the spatial particle-number density, and in collision velocities. For large particle, inertia the fluid velocity at the particle position is essentially a white-noise signal and caustic formation is analogous to Kramers escape. Here we show that caustic formation at small particle inertia is different. Caustics tend to form in the vicinity of particle trajectories that experience a specific history of fluid-velocity gradients, characterised by low vorticity and a violent strain exceeding a large threshold. We develop a theory that explains our findings in terms of an optimal path to caustic formation that is approached in the small inertia limit.

The dynamics of small, yet heavy, identical particles in turbulence exhibits singularities, called caustics, that lead to large fluctuations in the spatial particle-number density, and in collision velocities. For large particle, inertia the fluid velocity at the particle position is essentially a white-noise signal and caustic formation is analogous to Kramers escape. Here we show that caustic formation at small particle inertia is different. Caustics tend to form in the vicinity of particle trajectories that experience a specific history of fluid-velocity gradients, characterised by low vorticity and a violent strain exceeding a large threshold. We develop a theory that explains our findings in terms of an optimal path to caustic formation that is approached in the small inertia limit.
Ensembles of heavy particles in turbulence, such as water droplets in turbulent clouds [1] or dust grains in the turbulent gas of protoplanetary disks [2,3], may exhibit large fluctuations of the particle-number density and of their relative velocities [4][5][6][7]. These fluctuations are enhanced by the formation of caustics [8][9][10], i.e., folds of the particle distribution over configuration space. Caustic formation is an effect of particle inertia, driven by the fluid-velocity gradients, that gives rise to a multivalued particle-velocity field. Due to this multivaluedness, often called the 'sling effect' [11,12], particles may approach each other, and possibly collide, at large relative velocities. Accordingly, caustics have an important impact on the distribution of relative velocities [13][14][15][16], and are a crucial ingredient to theories for collision rates and collision outcomes [17][18][19]. In effect, caustic formation may increase the variance of the particle size distribution in turbulent aerosols because, on the one hand, caustics facilitate particle growth by enhancing collision rates [14,16]. Increased collision velocities may, on the other hand, lead to fragmentation, and thus to reduced particle sizes [3].
Caustics have been extensively studied in direct numerical simulations (DNSs) of particles in turbulence [12] and model flows [9,[20][21][22][23][24][25]. Recent numerical studies [26,27] found that high-velocity collisions tend to occur where the turbulent strain is large, but this cannot be explained in terms of the white-noise models usually used to study caustic formation [8][9][10]. A precise understanding of how caustics form at small particle inertia, including the local flow conditions that lead to their formation, is crucial for the identification of caustics in experiments [28] and for sampling them efficiently in DNS [12].
In this Letter, we describe a significant step towards a detailed understanding of how caustics form in turbulence. Using a DNS of two-dimensional turbulence, we show that whether a caustic forms or not depends on the history of the fluid-velocity gradients experienced by closeby particles, not just upon instantaneous correlations between particle positions and flow gradients (pref-FIG. 1. Illustration of caustic formation in two spatial dimensions. Shown are three nearby particles. A caustic is formed when faster particles overtake slower ones. At a caustic, the areaV of the parallelepiped spanned by separation vectors between the three particles, shown in grey, vanishes. erential concentration [7,29,30]). When particle inertia is small, we find a most likely history, i.e., an "optimal path" to caustic formation. To determine this path is an optimal-fluctuation problem, similar in nature to localisation due to optimal potential fluctuations in disordered conductors [31], population extinction due to environmental and population-size fluctuations [32,33], and shock formation in Burgers turbulence [34,35]. Based on this observation, we develop a theory that explains how the strain and vorticity change along the optimal path to caustic formation: The fluid strain performs a time-localised, violent fluctuation that exceeds a large threshold, while vorticity remains small. Our results explain qualitatively why DNSs of particles in turbulence show increased collision rates in straining regions [26,27]. Even at finite inertia, the optimal path leaves a clear mark in the data, providing criteria for the identification and the efficient sampling of caustics in experiments and in DNSs.
In a dilute suspension of small, heavy, spherical particles, the dynamics of a single particle is approximately given by Stokes' law [7], Here, x and v denote particle position and velocity; τ p = 2a 2 ρ p /(9ρ f ν) is the particle-relaxation time which depends on the particle size a, the kinematic viscosity ν of the fluid, and the particle and fluid densities, ρ p and ρ f , respectively. The turbulent fluid-velocity field, evaluated FIG. 2. Path to caustic formation for St = 0.31. Numerical results using a DNS of two-dimensional turbulence. The path density is colour coded, normalised to unity along each slice of constant Z to improve visibility. (a) Path density in the Z-A 2 plane. The dashed line shows the approximation, Z ∼ A 2 . The solid line shows the evolution of the stable fixed point according to Eq. (6). The dot marks the bifurcation to instability (see main text). (b) Individual contributions from strain (Z vs σ 2 , positive axis) and vorticity (Z vs −ω 2 , negative axis). (c) Path densities for σ 2 (positive axis) and −ω 2 (negative axis) as a function of time before caustic formation at t = 0.
at the particle position, is denoted by u[x(t), t].
To describe caustic formation, we consider the parallellepiped spanned by d + 1 nearby particles in d spatial dimensions. How the spatial volumeV (t) of this object contracts or expands under the nonlinear dynamics (1) is determined by the spatial Jacobian Since the dynamics (1) takes place in 2d-dimensional phase space, spatial subvolumesV may collapse in finite time,V → 0, when a caustic forms [7]. Figure 1 shows a typical particle configuration that leads to a caustic in two spatial dimensions. As is well known, caustic formation is closely related to the dynamics of the particle-velocity gradients, which reads, in dimensionless form [7,11], with initial condition Z(t 0 ) = A(t 0 ). Here, the Stokes number St = τ p /τ K is a dimensionless measure of particle inertia; Z ij = τ p ∂v i (t)/∂x j (t) and A ij = τ p ∂u i (t)/∂x j (t) are the dimensionless matrices of particle-velocity gradients and fluid-velocity gradients, respectively. In Eq. (2), time is dedimensionalised by the Kolmogorov time, τ 2 K = 2 i,j=1 (∂u i /∂x j ) 2 s = τ 2 p Tr(AA T ) s , where · s denotes a steady-state ensemble average. Here and in the following, we use the abbreviations Z(t) = Z[x(t), t] and A(t) = A[x(t), t]. Using (2), one finds [7] where Z = TrZ is the divergence of the field of particlevelocity gradients. Hence, a necessary condition for caustic formation is that Z escapes to negative infinity. Apart from the Reynolds number Re that specifies the turbulence intensity, the particle dynamics is determined by the Stokes number St. For small St, particle detachment is characterised by A − Z, which is typically of the order of St, and thus small. Caustic formation requires the activation of the nonlinear term in Eq. (2) that drives the particle-velocity gradients into a caustic. This, in turn, requires rare and violent fluctuations of the fluidvelocity gradients A of the order of unity.
We determine the dominant events that drive caustic formation by measuring the statistics of paths in the joint space of Z and the fluid velocity gradients A. In isotropic turbulence, the properties of any statistical quantity must be invariant under rotations. In addition to Z, we therefore map out the paths of the invariants obtained from the symmetric (S) and antisymmetric (O) parts of the fluid-velocity gradient matrix A: ω 2 = TrOO T and σ 2 = TrSS T . Figure 2 shows the paths to caustic formation obtained by numerical simulation, using a DNS of two-dimensional incompressible turbulence. Our simulations are performed in a periodic box forced at the large scales, in the regime of direct cascade of enstrophy. A drag-friction term ensures steady-state turbulence; see the Supplemental Material (SM) [36] for more details. The path density in Fig. 2 is colour coded, with the highest densities shown in yellow. Figure 2(a) shows paths to caustic formation in the [37,38] that discerns hyperbolic from elliptic regions in the flow. We see that most paths (yellow regions) that reach a caustic at Z = −∞ pass a large fluid-gradient threshold A 2 ≈ 0.2. The solid and dashed lines are explained in our analysis below. Figure 2(b) shows that typical paths to caustic formation correspond to large strain σ 2 . Vorticity ω 2 , by contrast, remains small for the majority of paths. Figure 2(c) shows the time evolution of σ 2 and ω 2 prior to caustic formation at t = 0. We observe that while ω 2 remains small, σ 2 increases sharply, reaches a large value, and then decreases again. The majority of the large strain, however, persists until the caustic is formed, suggesting that caustics preferentially form in regions of large strain. Appealing to optimal-fluctuation theory, our numerical results point towards an optimal path that underlies caustic formation, characterised by small vorticity and a violent strain. Although the spread in our data is quite large at the value of the Stokes number we used, St = 0.31, the optimal path leaves a strong mark in our data, reflected by the yellow streaks in Fig. 2.
We explain our observations using an optimal fluctuation approach. The first step is to analyse the fixed-point structure and the bifurcations of Eq. (2). To this end, we expand the equation of motion (2) for the 2 × 2 matrix Z in a basis of matrices generated by the identity matrix I and This basis is orthogonal with respect to the inner product defined for two matrices M and N by M, N = 1 2 Tr[MN], so that e i , e j = g ij = g ij = diag(−1, 1, 1) ij and e T i = e i = g ij e j . Here and in the following, we use the Einstein sum convention. We denote the three-vectors corresponding to Z and A by z respectively. This formulation in terms of the Lorentzian metric g ij [39] is convenient because it disentangles the strain and vorticity parts of the fluid-velocity gradients A. We have O = ωe 1 , so that A 1 = ω describes the vorticity. The other components, A 2 and A 3 , describe the strain, S = A 2 e 2 + A 3 e 3 . Similarity transformations of Z,Z = PZP −1 leave Z invariant, but transform z i by means of a proper Lorentz transformation,z i = Λ i j z j . The same holds for transformations of A. The matrix Λ which transforms A i and z i has the properties Λ T gΛ = g, and det Λ = 1.
Expanding Eq. (2) in the basis (4), we obtain As the time derivatives on the left-hand side of Eq. (5) are multiplied by St 1, we expect the dynamics of Z and z i to take place in the vicinity of its stable fixed points, if they exist. For A i = 0, we find three fixed points, Z = z i = 0, Z = −2, z i = 0, and Z = −1, z i z i = 1/4, whose stability is determined by the eigenvalues of the stability matrix of (5). The fixed point Z = z i = 0 is stable for where the fixed point disappears. We conclude that when A 2 < 1/8, the dynamics (5) takes place in the vicinity of the stable fixed point obtained from the implicit equation When A 2 = 2A i A i exceeds 1/8, the fixed point ceases to exist, and the nonlinear dynamics (5) drives Z to negative infinity, forming a caustic. The evolution (6) of the stable fixed point as a function of A 2 and σ 2 is shown as the solid lines in Figs. 2(a) and 2(b). The fixed points become unstable at the bifurcation point, Hence, for A 2 < 1/8, particle neighbourhoods are stable and are continuously deformed by the fluid-velocity gradients, according to Eq. (6). For A 2 > 1/8, however, the neighbourhoods become unstable and collapse after a short time. Expanding Eq. (6) for small Z, one obtains the approximation Z ∼ −A 2 [dashed lines in Figs. 2(a) and 2(b)] used by Maxey [29] to explain the preferential concentration of heavy particles in incompressible turbulence [7,30]. This approximation fails to describe caustic formation because it predicts that Z remains finite, and thus leads to the incorrect conclusion that particle neighbourhoods are always stable.
Our stability analysis of Eq. (5) explains the qualitative shape of the paths in Fig. 2(a). However, it misses some of the important results of our DNS. In particular, the stability analysis does not explain why only the strain contributes to caustic formation and vorticity remains small [ Fig. 2(b)], and it has no bearing on the time evolution of the large gradient fluctuations shown in Fig. 2(c). Finally, we observed in Fig. 2(a) that the threshold reached by most paths is actually slightly larger than 1/8 = 0.125, the value predicted by our stability analysis.
In order to explain these parts of our observations we need to go beyond the stability analysis, and consider how the fluid-velocity gradients reach the large threshold required to render particle neighbourhoods unstable. We do this in the following by computing their optimal fluctuation.
The steady-state correlation functions of S and O, evaluated along particle trajectories in isotropic and homogeneous turbulence, have the general form and S ik (t)O jl (t ) s = 0. In two spatial dimensions, the tensors in Eqs. (7) are given by C S ijkl = 1/4(δ ij δ kl + δ il δ jk − δ ik δ jl ) and C O ijkl = 1/2(δ ij δ kl − δ il δ jk ). We express Eqs. (7) in terms of the basis (4) to obtain the steady-state correlations of A i , All other correlations of A i are zero in the steady state. The right-hand sides of Eqs. (8) are parametrised as with the non-dimensional correlation times s and o of S and O, respectively. For tracer particles with St = 0 one has C S (0) = C O (0) = 1/2. Inertial particles with St > 0 tend to avoid vortical regions due to preferential concentration [29,30], so that C O (St) < 1/2. The amplitude C S , on the other hand, remains approximately equal to 1/2 [26]. In our two-dimensional numerics, we also find C S (St) ≈ 1/2 for Stokes numbers between 0.21 and 0.51. . The colour-coded path density is the same as in Fig. 2(b). (b) σ 2 opt (t) from theory (dashed line). The colour-coded path density is the same as in Fig. 2(c).
To describe how the fluid-velocity gradients reach the required threshold values, we model A i (t) as independent, stationary Gaussian processes with zero mean. For this class of processes, the most probable (optimal) fluctuation A i opt (t) to reach a given threshold can be obtained by optimal-fluctuation methods, as we show in the SM [36]. By minimising the action associated with the path probability, we find that the optimal fluctuation of the fluid-velocity gradients is free of vorticity, ω 2 opt = 0, in agreement with our DNSs in Figs. 2(b) and (c). This result is intuitive: Vorticity contributes to A 2 with a negative sign, so that any fluctuation of A 2 = σ 2 − ω 2 that reaches the large threshold A 2 th with finite vorticity requires an even larger strain contribution, to make up for vorticity. The optimal way to reach the threshold value is thus through paths that are vorticity free, whereas the probabilities of paths with finite vorticity are exponentially suppressed. The optimal path for the strain, by contrast, is found to be a time-localised fluctuation, σ 2 opt = A 2 th e −2(t−t th ) 2 /s 2 , given by a Gaussian function peaked at time t th < 0 [36].
Using the optimal gradient fluctuation (σ 2 opt , ω 2 opt ), we now obtain the explicit form for the optimal path (Z opt , σ 2 opt ) as a function of time. For St 1, the lefthand sides of Eqs. (5) are small most of the time. To evaluate Z opt , we therefore use A i opt as an input into Eq. (5). Making use of the fact that the vorticity is zero along the optimal path, ω opt = 0, we find that A opt = A i opt e i and Z opt = 1/2Z opt +z i opt e i can be brought into diagonal form by a Lorentz transformation [36]. The equations for the diagonal entries (eigenvalues) λ ± opt of Z opt decouple into two equations, with initial conditions λ ± (t 0 ) = ±A th e −(t0−t th ) 2 /s 2 / √ 2. The uncoupled Eqs. (10) are solved numerically, which yields the optimal path Z opt using Z opt = λ + opt + λ − opt . We note that for finite St, the threshold value A 2 th , determined numerically [36], exceeds the value 1/8 obtained from the stability analysis, in agreement with our DNS. The reason is that the optimal strain fluctuation σ 2 opt (t) decreases for t > t th . In order for Z opt to reach negative infinity, σ 2 opt (t) must exceed 1/8 for a finite time so that Z opt can become large and negative. The time for which the threshold must exceed 1/8 decreases as St becomes smaller, and we recover A 2 → 1/8 in the limit St → 0.
In Fig. 3, we compare our theory and DNSs. The dashed line in Fig. 3(a) shows (Z opt , σ 2 opt ) obtained from Eqs. (10), with St = 0.31 and dimensionless correlation time s = 2.1 determined numerically. We observe qualitative agreement between the theoretically obtained optimal path and the (yellow) regions of high path density. However, Eqs. (10) slightly overestimate the threshold value A 2 th . The likely reason is that the Stokes number in our DNS is too large to closely follow our analytical results, valid for St 1. Figure 3(b) shows the time dependence of σ 2 opt (t) obtained from theory with St = 0.31 (dashed line) and the corresponding path density from our DNS. We observe that the theory correctly predicts the localised, violent fluctuation of the strain prior to the formation of the caustic, which occurs at t = 0. A considerable fraction of the large strain required to initialise the caustic persists at the time of caustic formation, which explains why caustics form preferentially in regions of large strain [26,27]. For other Stokes numbers between 0.24 and 0.51, our DNS results lead to the same conclusion (data not shown).
In conclusion, we explained caustic formation in turbulent aerosols at small particle inertia as an optimalfluctuation problem. In order for caustics to form, the fluid-velocity gradients must follow an optimal path, characterised by small vorticity and a violent strain that exceeds a large threshold. The remnants of the optimal fluctuation at the time of caustic formation result in a strong instantaneous correlation between large strains and caustic events.
Since caustics give rise to a multivalued particlevelocity field, and thus to high relative particle velocities, our results provide an explanation for the recently observed, instantaneous correlation between particle collisions and intense strain [26,27]. The characteristic shape of the optimal path to caustic formation will allow one to identify caustics in experiments, and the strong instantaneous correlation of caustics and strain makes it possible to efficiently sample caustics in simulations.
The stability analysis described in this Letter can be generalised to three dimensions, where it reveals that particle neighbourhoods become unstable when the two invariants Q = −TrA 2 /2 and R = −TrA 3 /3 reach large thresholds in the Q-R plane. We therefore speculate that optimal-fluctuation methods also explain caustic formation in three-dimensional turbulence at small Stokes numbers. In this supplemental material we provide additional details regarding the analytical calculations summarised in the main text, and for the numerical simulations used to obtain the results shown in Figs. 2 and 3 in the main text.

I. INSTANTON CALCULATION FOR THE FLUID-VELOCITY GRADIENTS
In this section we explain in detail the instanton calculation outlined in the main text. We consider the problem of the gradient process A i reaching a threshold value given by In particular, we explain how to derive the time dependence of the optimal paths for the gradient processes A i opt , ω 2 opt and σ 2 opt . We model the individual components A i of the fluid-velocity gradient matrix A by independent, stationary Gaussian processes with zero mean and steady-state correlation functions [cf. Eqs. (8) and (9) in the main text], (S1) Figure 1 shows our numerical data for the correlation functions. Figure 1(a) shows that TrO(t)O(0) s is well approximated by a Gaussian function, f (t) = exp(−t 2 ). The correlation function of O(t), shown in Fig. 1(b), is close to an exponential function, so that we choose g(t) = exp(−t). In Eqs. (S1), o and s are the dedimensionalised correlation times defined in the main text. From the gray line in Fig. 1  To obtain the optimal fluctuation to reach a threshold, we consider the transition probability of the gradients A i from A i (t 0 ) = 0 at time t 0 , to a large threshold value A i (t th ) = a i at time t th . Because the processes A i are independent, we compute this transition probability for a single process A(t) with arbitrary correlation function. For the Ornstein-Uhlenbeck process, it is well known [1] how to compute the optimal path to reach a threshold. Here we use a similar method, but generalised to any stationary Gaussian process. The processes we consider need not be Markovian as the Ornstein-Uhlenbeck process, and may not have a representation in terms of a stochastic differential equation.
The transition probability of the zero-mean process A(t) from zero at time t 0 , to a at time t th can be written in terms of a path integral where G[z(t)] denotes the generating functional of A(t) and N is a normalisation factor. For Gaussian processes, Here c t0 is the conditional correlation function of A, defined as The conditioning implies c t0 (t 0 , t) = c t0 (t, t 0 ) = 0 for all t ≥ t 0 . As t 0 → −∞, A becomes stationary, so that Here c s (t − t ) denotes the steady-state correlation function of A(t) which also defines the variance ε 2 of the process, c s (0) = ε 2 . Using the above expressions, we write Eq. (S2) as with action S t0 [z, A] given by We consider weak-inertia case St 1 which corresponds to ε 1. In this case, the path integration in Eq. (S7) is well-known to concentrate around an 'optimal path' (z opt , A opt ) [3,4]. This path is obtained by a variation of the action S t0 [z, A] over z(t) that vanishes at the optimal path: This determines the optimal path A opt for A(t): The transition probability (S2), and thus the action S t0 [z, A], depends on A(t) only through the final point, A(t th ) = a. As a consequence, z opt (t) has the general from The constant C is evaluated by using the boundary condition A opt (t th ) = a. We find , and thus, A opt (t) = a c t0 (t, t th ) c t0 (t th , t th ) . (S11) Upon substituting Eqs. (S11) and (S10) into Eq. (S7), the action S t0 [z opt , A opt ] evaluated along the optimal path (z opt , A opt ) reads . (S12) Minimising S t0 [z opt , A opt ] over initial times t 0 [1] gives t 0 → −∞. This yields for t ≤ t th , (S13) For t > t th , the process is unconstrained and relaxes from a to zero as t → ∞. Unconstrained relaxation does not contribute to the action [5]. Furthermore, for time-reversal invariant processes A(t), the optimal path for the relaxation from the threshold a must be identical to the time-reversed optimal trajectory A opt for reaching the threshold. This argument gives A opt for all t, (S14) In other words, for a time-reversible Gaussian, stationary process with zero mean, the optimal path for reaching a threshold a is, up to a constant, given by the steady-state correlation function. The action S depends only on the variance ε 2 , and is otherwise independent of the details of the process A(t). For the Ornstein-Uhlenbeck process, our results agree with those in Ref. [1] based on the Onsager-Machlup action [6]. We now be apply Eqs. (S13) and (S14) to the specific, independent processes A i with steady-state correlation functions (S1) and f (t) = exp(−t 2 ) and g(t) = exp(−t). Comparing Eq. (S14) with (S1) we have The parameters a i are obtained by minimising the action S(a i ) in Eq. (S16), From this result, we obtain as stated in the main text.

II. SOLUTION TO THE EIGENVALUE EQUATIONS FOR Z
We now show how we obtain our solutions for the optimal path for Z = TrZ = λ + + λ − for small St 1 [see main text below Eq. (10)]. We first explain how to derive Eq. (10) in the main text. We then briefly discuss our numerical method to determine the threshold value A 2 th , and to produce Fig. (3) in the main text. To this end, we write At time t = t 0 we have Z(t 0 ) = A(t 0 ), so that Z opt (t 0 ) = 0. We conclude that z 1 opt (t) = 0 for all times. Now, we perform a Lorentz transform, which in this case is just a simple rotation, of the remaining components The leaves the general form of Eqs. (S20) invariant, so that we obtain after the transform FromZ(t 0 ) =Ã(t 0 ) we again concludez 2 opt (t) = 0 for all times. Using the initial decomposition we obtainZ opt = Z opt /2I +z 3 e 3 , where e 3 is diagonal with elements e 3 = diag(1, −1). Hence, we find that after the transformation, Z opt is in diagonal form with eigenvalues λ + opt = Z opt /2 +z 3 opt and λ − opt = Z opt /2 −z 3 opt . Taking time-derivatives of these eigenvalues, and using Eqs. (S21), we obtain the uncoupled equations for the eigenvalues λ + opt and λ − opt , presented in Eq. (10) in the main text. As the trace of a matrix is the sum of its eigenvalues, we have Z opt = TrZ opt = λ + opt + λ − opt . These equations need to be solved numerically with the initial conditions λ + opt (t 0 ) = A th √ 2 e −(t0−t th ) 2 /s 2 and λ − opt (t 0 ) = − A th √ 2 e −(t0−t th ) 2 /s 2 . To determine the threshold value A th , we start at A th = 1/8 and increase A th incrementally until λ − opt escapes to negative infinity. By returning to the previous A th , and reducing the step size, we obtain an iterative approximation of the smallest value A th for which λ − opt → −∞. When computing the optimal trajectories shown in Fig. (3) in the main text, the threshold value A th must be exceeded slightly, in order to enable λ − to escape to negative infinity in finite time. The amount by which the threshold is exceeded is a free parameter in our model. Motivated by the fact that Gaussian fluctuations around the optimal path A i opt are of the order of the square-root of the variance (St 2 C s (St)/4) 1/2 , determined by the action (S17), we choose this parameter equal to St/(2 √ 2). We have checked that our results are insensitive to small changes of this parameter value.

III. DIRECT NUMERICAL SIMULATIONS
In this section, we provide some details of the numerical method we used to solve the two-dimensional Navier-Stokes equations and the particle dynamics, resulting in Figs. 2 and 3 in the letter.