Extreme spatial clustering by fractal catastrophes

We analyse the spatial inhomogeneities ('spatial clustering') in the distribution of inertial particles accelerated by a space-time dependent random force. To quantify spatial clustering, the phase-space dynamics of the particles must be projected to configuration space, resulting in mathematical catastrophes ('caustics') whose quantitative contribution to spatial clustering is not understood. We show how to solve this problem by projecting the phase-space finite-time Lyapunov exponents (FTLEs) to configuration space. Applying our method in one spatial dimension, we find that in the long-time limit 'fractal catastrophes', caustics that arise from the projection of a dynamical fractal attractor, make a distinct and universal contribution to the distribution of spatial FTLEs. This contribution gives rise to an extreme form of spatial clustering, and breaks a fluctuation relation for white-in-time Gaussian force fields.

Particles in physical environments are often subject to forces that appear to fluctuate randomly. Examples are particles in turbulence, such as water droplets in clouds [1] or dust in the turbulent gas of protoplanetary disks [2], but also small floaters on the free surface of a fluid in motion [3], charged particles accelerated in the atmosphere [4][5][6], Chladni patterns of dust grains on randomly vibrating plates [7], or particles diffusing in spatially varying temperature fields [8]. When the particle momenta are damped by friction, the phase-space dynamics is dissipative, leading to the formation of fractal patterns in the particle distribution [3,9,10].
The resulting strong inhomogeneities in the spatial distribution of particles ('spatial clustering') have been studied in detail in experiments [11][12][13]. Spatial clustering is of key significance because it brings particles close together, and thus affects particle interactions through collisions [14], evaporation or condensation [15], or chemical reactions [16]. Fractal clustering is quantified by the large-deviation statistics [17][18][19] of finite-time Lyapunov exponents (FTLEs) [20][21][22][23][24][25][26], measuring the evolution of infinitesimal volumes spanned by nearby particle trajectories. The distribution of the FTLEs allows to compute the long-time statistical properties of the dynamically evolving fractal attractor to which the particles converge.
The inertial particle dynamics occurs in phase space, and the statistical properties of the phase-space attractor are determined by the phase-space FTLEs. To describe spatial clustering one must project the phase-space FTLEs to configuration space. Since it is not understood how this projection affects the distribution of FTLEs, a first-principles theory of spatial clustering is missing, despite the significance of the problem.
The source of the difficulties is well known [10,[27][28][29][30], and illustrated in Fig. 1: because the inertial phasespace dynamics generates folds [ Fig. 1(a)], the spatial projection becomes many-to-one, causing spatially projected, infinitesimal neighbourhoods of particles to collapse to single points in configuration space [ Fig. 1(b)]. These singular points are cusp or fold catastrophes, also called 'caustics' [32,33]. For smooth phase-space manifolds, caustics are known to lead to finite-time singu- FIG. 1. Fractal phase-space and configuration-space clustering (schematic). (a) Fractal attractor in phase-space (position x, momentum p) at given time t [31]. (b) Increased spatial particle density (x, t) in the vicinity of caustic folds. (c) The magnification illustrates fractal clustering. larities in the spatial particle density, suggesting that caustics may increase spatial clustering [33]. These considerations are, however, too imprecise to quantify how caustics affect spatial clustering. More importantly, they rest on the notion of a smooth manifold, and it is unclear to which extent they apply to fractal phasespace attractors [ Fig. 1(c)]. Most existing theories aim at quantifying spatial clustering by perturbation theory [10,23,24,[34][35][36][37][38]. As the effect of caustics is nonperturbative [34,39,40] these theories naturally fail, already at small perturbation parameters. Up to now, these complications have made a first-principles description of spatial clustering impossible.
In this Letter, we show how to compute the statistics of spatial FTLEs for inertial particles accelerated by a smooth, random force field, and demonstrate the method in one dimension. As the main result, we find that caustics lead to an extreme form of spatial clustering which manifests itself in an exponentially increased probability of large negative spatial FTLEs, independent of the details of the force field. We show the implications of these conceptual insights for white-in-time Gaussian force fields and discuss how they apply to higherdimensional systems.
Consider the dynamics of the position x t and momentum p t of a particle of mass m in d spatial dimensions, At time t and for δx t much smaller than the correlation length of f , we can linearise the force field around (x t , p t ) to obtain the linear dynamics with the force-gradient matrix F ij = ∂ j f i . The solution of Eq.
(2) is expressed in terms of the Green function J(t) , and T exp denotes the time-ordered exponential, evaluated along (x, p). Writing J = V R, we decompose J into a rotation R and a stretch tensor V with positive and real eigenvalues (e σ (1) t t , . . . , e σ (2d) t t ). The exponents σ (i) t are the phase-space FTLEs discussed above, conventionally arranged in non-increasing order [24]: The second equation reflects that tr A(t) = −γd is constant [10,32]. The sub-sums ∑ n j=1 σ (j) t , n ≤ 2d, are the rates of growth and contraction of n-dimensional phase-space volumes spanned by n+1 nearby trajectories. For ergodic dynamics, the FTLEs have definite limits, λ i = lim t→∞ σ (i) t , called Lyapunov exponents [41]. We define the left Cauchy-Green tensor B ≡ JJ T = V 2 with eigenvalues (e 2σ (1) t t , . . . , e 2σ (2d) t t ). Studying B instead of V is convenient since the former obeys the closed equationḂ = BA T + AB. As shown in [22] this allows to formulate evolution equations for the FTLEs and for the orthogonal matrix O that diagonalises B. For long enough times and given a non-degenerate spectrum of Lyapunov exponents, the dynamics of O decouples from that of the FTLEs [22]. Expressing O in terms of simple rotations , leads to a closed stochastic equation for α t . The phase-space FTLEs are then solutions of stochastic integrals over α t of the form σ ii . This method works in principle in any spatial dimension d and for any force gradient F(x t , t).
As the number of Euler angles grows quadratically with d, we consider d = 1, noting that higher-dimensional generalisations are conceptually straightforward but require much more computational effort. For d = 1, on the other hand, there is only one Euler angle α t ≡ α (1) t , and one independent phase-space FTLE, σ Evolution of a small, two-dimensional phase-space volume, grey area, around a reference trajectory (x, p), dashdotted line. The two perpendicular axes are stretched or contracted by ∼ e tσ (i) t , and rotated by αt. The dashed line shows the projection ∼ e tσ t of the stretching direction to configuration space. The curved arrow indicates the direction of the steady-state flux. (b) Steady-state probability density of αt for white-in-time force fields with different ε.
Appendix A for details): Here, F t ≡ F (x t , t) and σ (2) t is obtained from the second constraint in Eq. (3). Since equation (4a) is invariant under the shift α t → α t + kπ, k ∈ Z, we impose periodic boundary conditions on [−π 2, π 2). Eqs. (4) admit the following interpretation, shown schematically in Fig. 2(a): Consider a two-dimensional phase-space volume, say a circle of unit area at t = 0. At t > 0, the circle is deformed by e tσ (1) t along one eigendirection of B and by e tσ (2) t along the other one. The initial circle is thus squeezed into an ellipse with decreased phase-space volume V t ∼ e −γt , due to the dissipative nature of the dynamics. At the same time, the eigensystem of B rotates by α t .
The system (4) is more conveniently written in terms of Z t = γ tan α t = m −1 δp t δx t [27,28], which measures the local particle-velocity gradient along the reference trajectory. We obtain The particle-velocity gradient regularly escapes to −∞ and immediately reappears at +∞ [27,28] which corresponds to a local fold of a phase-space neighbourhood over coordinate space, a caustic [33]. Note, however, that σ (1) t remains finite at all times. In terms of α t , a caustic corresponds to α t transitioning from −π 2 to [π 2] + . Because d dt α −π 2 = −γ, this transition happens deterministically with speed −γ. If F t is statistically stationary, α t reaches a non-equilibrium steady state with density P st (α t = a), which has a finite flux of magnitude J, the caustic-formation rate [28,29]. This implies P st (α t = −π 2) = P st (α t = [π 2] + ) = J γ.
To derive the statistics of the spatial FTLEσ t , we must project (σ To this end, consider the spatial projection of a neighbourhood δx t , δpt mγ ↦ δx t around (x t , p t ), at time t. The projection is represented as the dashed line in Fig. 2(a). After an initial transient, the absolute value δx t of the projected phase-space separation scales as At large t, α t and σ (1) t evolve on different time scales. As discussed above, α t passes −π 2 with rate J. The large-deviation principle (6), on the other hand, implies that σ (1) t stabilises in the vicinity of the Lyapunov exponent λ 1 . Therefore, the joint distribution of α t and σ (1) t factorises, so that the cumulative distribution function P (σ t ≤ s) can be written as (see Appendix B) Using P st (cos α t ≤ x) ∼ 2Jx γ for x ≪ 1 and Eq. (6), we find thatσ t obeys a large-deviation principle with rate functionÎ(s) that is closely related to the phase-space rate function I(s) (see Appendix B for details): with s 0 = argmax s∈R {−s − I(s)} and Eqs. (9) and (10) are the main results of this Letter. They show that the projection of the phase-space FTLEs to configuration space, (σ t ) ↦σ t , generally preserves the large deviation form ofσ t , but alters the rate func-tionÎ(s). Thus, for finite caustic formation rate J > 0, I(s) scales linearly in s for s < s 0 . This linear scaling is a universal contribution from the spatial projection and independent of the details of the force field. Botĥ I(s) and I(s) are convex, so thatÎ(s) < I(s) for s < s 0 . Hence, the probability of large, negativeσ t is exponentially enhanced by a factor of ∼ e t[I(s)−Î(s)] compared to the probability of observing large, negative phase-space FTLEs. Negative spatial FTLEs correspond to spatially localised clusters of particles, so that caustics lead to an extreme increase in the probability of clustering events.
Eqs. (9) and (10) have important consequences, and allow for a quantitative description of the caustic contribution to the statistics of (x, t), as we show in the following. We perform a Legendre transform ofÎ(s) which yields the spatial SCGF: (7) and (5b),σ t can be expressed aŝ as the the local, spatial particle density around (x t , p t ) at time t. As Λ(−1) is a function of ⟨ t ⟩, the expression for k = −1 in Eq. (11) follows from mass conservation (see Appendix C). Furthermore, Λ(k) = ∞ for k < −1 implies that ⟨ k t ⟩ = ⟨ δx −k ⟩ grows, for k > 1, super-exponentially in time. The moments of the fixed-frame spatial density (x, t) [ Fig. 1(b)] are obtained from t (measured along trajectories), by weighting each trajectory by its phase-space volume V t ∼ e −γt [24]. Since V t is independent of the trajectory, extreme spatial clustering entails that also ⟨ (x, t) k ⟩ grows superexponentially in time for k > 1. By contrast, in the overdamped limit [24] and in perturbative schemes [23] spatial density moments grow at most exponentially.
Let us analyse the consequences of Eqs. (9)-(11) for fractal spatial clustering. The spatial and the phasespace correlation dimensions,D 2 and D 2 , satisfy the implicit equationsΛ(−D 2 ) = 0 [45,46] and Λ(−D 2 ) = 0 [37,38,47]. Using this, Eq. (11) immediately im-pliesD 2 = min{D 2 , 1}, which proves previous conjectures on the basis of numerical investigations in turbulent aerosols [10,32,36]. The saturation ofD 2 at unity is a consequence of the linear part inÎ(s) for s < s 0 , and thus a result of the caustics. We conclude that when caustics events are rare, the projection has no effect on the correlation dimension andD 2 = D 2 . Above a set of critical values of the model parameters, however, the singularities shown in Fig. 1(b) become space-filling, resulting in the saturation.
For the phase-space dynamics (1), which generate the dynamical fractal attractor illustrated in Fig. 1(a), the statistics of FTLEs determines not only D 2 but the whole spectrum of fractal phase-space dimensions, D q [20,21,25]. This spectrum characterises the inhomogeneity of the fractal measure [48][49][50]. We consider the moments of the probability (mass) M R contained in a small phase-space ball of radius R = R = δx 2 + (δp mγ) 2 around (x, p). These moments scale as ⟨M n R ⟩ ∼ R ξn , for R ≪ 1 [25,32]. The so-called singularity exponents ξ n are related to the Renyi dimensions D q [48][49][50] by ξ n = nD n+1 . Adapting the formalism described in Ref. [25] to our model, we find that, for a positive maximum Lyapunov exponent λ 1 > 0, the exponents ξ n are given by For the special case of white-in-time force fields, we now show how to compute Λ(k) and I(s) explicitly. This allows to study in more detail the impact of the spatial projection onÎ(s), and to evaluate Eqs. (12). When f (x, t) has Gaussian statistics with zero mean and vanishing correlation time, the gradient F t is a Gaussian white noise with zero mean, and correlation ⟨F t F t ′ ⟩ = 2D F δ(t − t ′ ), where D F is the diffusion constant. In this case, the model depends solely on the dimensionless parameter ε 2 = m −2 γ −3 D F , and the dynamics of δx t , δpt mγ decouples from that of (x t , p t ) [10]. From the Fokker-Planck equation corresponding to Eq. (4a) (interpreted in the Stratonovich sense) we obtain the steady-state density P st (α t = a) analytically [28,51]. Figure 2(b) shows P st (α t =a) as a function of a for different values of ε. For small ε, α t stays close to zero most of the time, but is more and more likely to approach −π 2 as ε increases. Note that P st (α t = −π 2) = P st α t = [π 2] + = J γ as anticipated above, where J is obtained from the normalisation condition of the density [10].
Since F t is white in time, Λ(k) can be calculated as the leading eigenvalue of the tilted generator L k = (z 2 γ −2 + 1) −k 2 L (z 2 γ −2 + 1) k 2 + kz [52] associated with the large-deviation statistics [19,42,43] of σ (1) t , where L = −z(z + γ) d dz + ε 2 γ 3 d 2 dz 2 is the infinitesimal generator of the process (5a). We assume the leading eigenvalue of L k and its adjoint L † k to be unique and real [53]. In this case, and under certain conditions on the right and left eigenfunctions r k and l k [42], we have Since the dynamics for α t smoothly transitions from −π 2 to [π 2] + we demand that all eigenfunctions r k and l k are symmetric for large z . With help of (13), we can formulate a fluctuation relation for Λ(k) and, using the Legendre transform, carry it over to I(s) (see Appendix D), so that Note that I(s) has a reflection point at s = −γ 2. As shown in [54], this is a general property of pairs of FTLEs in dissipative systems of the form (1). The fluctuation relation (14b) connects the probability of clustering events (σ (1) t < 0) to the probability of voids (σ (1) t > 0) and follows from the statistical invariance of the force field F t under time reversal [54,55] (see also Appendix E). By taking a derivative of (14b) with respect to s at s = 0, we obtain I ′ (−γ 2) = −γ, so that s 0 = −γ 2 andÎ 0 = Λ(−1). We thus find from Eqs. (9)- (11), that the fluctuation relations (14) are not fulfilled by the projected quantitiesÎ(s) andΛ(k). This implies that the projection to configuration space destroys the balance (14) by introducing additional clustering due to caustics.
We solve Eqs. (13) numerically by a shooting method similar to that described in [30,56] for general k. For even integer k we use an exact analytical scheme (see Appendix F). Figure 3(a) shows the resulting Λ(k). As expected, the SCGF is convex [18,19], yet not a simple parabola as obtained from perturbation theory [22,24].
The explicit calculation of Λ(k) and I(s) allows us to obtain the singularity spectrum ξ n from Eq. (12). Figure 3(b) shows ξ n for 0 ≤ n ≤ 2. We observe that ξ n increases as a function of n and levels off to ξ n = ξ ∞ for n > n crit . The values of ξ ∞ and n crit can be expressed in terms of ξ 1 and Λ(−1) (see Appendix G). That ξ n is a non-linear function of n implies anomalous scaling of the phase-space mass moments ⟨M n r ⟩ in n. In other words, the particle distribution has non-Gaussian tails, even though the driving force is Gaussian [25,32]. As a result, the phase-space attractor [ Fig. 1(a)] is multifractal, in accordance with the numerical observations in Ref. [32].
In summary, we have shown how to quantify the effects of fractal catastrophes (caustics from projecting a dynamical fractal attractor) upon spatial clustering of inertial particles in a random force field. For one spatial dimension, we showed that these caustics lead to an exponential increase of the probability to observe large negative spatial FTLE, a universal law of extreme spatial clustering. The universal linear scaling of the left FLTE tail implies that spatial density moments above a critical order grow super -exponentially in time, and that the spatial correlation dimensionD 2 obeys a projection formula.
For white-in-time Gaussian force fields we could calculate the distribution of phase-space FTLEs explicitly. It exhibits a fluctuation relation, associating the probabilities of phase-space regions with large positive and large negative FTLE, voids and clusters. In the distribution of the spatial FTLE we showed that this balance is destroyed, a consequence of the extreme clustering.
An important question is how the conceptual insights obtained in one dimension carry over to systems in higher dimensions. A complete analysis, as provided for d = 1 in this Letter, poses formidable challenges because it involves a much more complex angular dynamics. However, our results should extend to the growth rate t −1 logV t of an infinitesimal d-dimensional spatial volumeV t [23,24], quantifying the long-time probability of observing local particle-rich regions in configuration space. We speculate that its rate function has a universal linear negative tail, resulting in extreme spatial clustering as in the onedimensional case.

ACKNOWLEDGMENTS
We thank Stellan Östlund and Anshuman Dubey for their comments on the manuscript. JM, KG, and BM were supported by the grant Bottlenecks for particle growth in turbulent aerosols from the Knut and Alice Wallenberg Foundation, Dnr. KAW 2014.0048, and in part by VR grant no. 2017-3865. Computational resources were provided by C3SE and SNIC.
Appendix A: Derivation of the dynamics of αt and σ (1) t , Eqs. (4) We start with the long-time evolution equations derived for O and σ (i) t in Ref. [22]: where In d = 1, we can parametrise O using a single angle α t . The matrices O and A [Eq. (2) in the main text] take the forms: Using Eq. (A2) we first express Ω as Ω = T(OAO T ) 21 so that the second equation in Eq. (A1) reads TO d dt α t = TO(OAO T ) 21 . Using this, we can read off the equation of motion for α t , d dt α t = (OAO T ) 21 . In conclusion, we obtain from Eqs. (A1) and (A4) the dynamics d dt α = −γ sin αt(sin αt + cos αt) These equations are equivalent to Eqs. (4) in the main text. It is easy to check from Eq. (A4) that d dt [tσ In the main text, we argued that the joint distribution of α t and σ (1) t factorises in the long time limit. In this case, the cumulative distribution function of α t conditional on a certain value of σ (1) t simplifies to We first use this expression to derive Eq. (8) in the main text. We have where we have used Eq. (7) in the second step. Conditioning on the value of σ (1) t we now write In the second step we have used Eq. (B1) and the large deviation form of σ (1) t , Eq. (6) in the main text. Eq. (B3) is Eq. (8) in the main text. Note that ifσ t obeys a largedeviation principle with rate functionÎ(s), the cumulative distribution function is of the form From the expressions P st (α = −π 2) = P st (α = [π 2] + ) = J γ derived in the main text, it follows for x ≪ 1 as stated in the main text. In order to obtain P (σ t ≤ s), we distinguish two cases in Eq. (B3): s ≥ −γ 2 and s < −γ 2.

Case s < −γ 2
For s < −γ 2 we have s < s ′ and thus e t(s−s ′ ) ≪ 1 in the integral in Eq. (B3). So that we can write Eq. (B3) using Eq. (B5) as where we have kept the constants to show that the linear part in the rate function disappears in case J = 0. Note that N t comes from the sub-leading terms in the saddle point approximation of the integral.

Case s ≥ −γ 2
In this case, we need to split the integral into two parts, where s > s ′ and s < s ′ , so that e t(s−s ′ ) ≫ 1 and e t(s−s ′ ) ≪ 1, respectively. We assume the case s = s ′ has measure zero. We obtain, similarly to before Which one of the two terms is leading, depends on location of the infimum of I(s) + s. As in the main text, we call s 0 = argmin s ′ ∈R {I(s ′ ) + s ′ }. If s 0 < s, then inf s ′ ≥s {I(s ′ )+s ′ } = I(s)+s, and the first term in Eq. (B7) is leading. If, on the other hand, s 0 ≥ s, the second term is leading. We conclude, for general s and using Eq. (B6), where the constantÎ 0 depends on the location of s 0 relative to −γ 2. We havê These are Eqs. (9) and (10) in the main text. Note that I(s) is a continuous and convex function of s. However, if s 0 < −γ 2 then it is not differentiable at s = −γ 2.
Appendix C: Calculation ofΛ(k), Eq. (16) The spatial scaled cumulant-generating functionΛ(k) is determined by Legendre transform ofÎ(s): If s 0 ≥ −γ 2 thenÎ(s) is differentiable and we readily obtain Eq. (11) in the main text: The situation becomes a bit more complicated when s 0 < −γ 2. In this case,Î(s) is not continuously differentiable at s = −γ 2 and we must consider the left and right limit of the derivatives, which gives a range of k values for whichΛ(k) is linear.
Appendix D: Derivation of fluctuation relation (14) using the tilted generator The tilted generator and its adjoint obey the eigenvalue equations The idea is to bring the first equation into the same form as the second one by a suitable change of variables and then compare the corresponding largest eigenvalues Λ(k). We first apply the transformations and introduce the change of variables z → y = z + γ 2.
The functionsr k (y) andl k (y) then obey the eigenvalue equations Now, we transform y → −y and k → −k − 2 in Eq. (D5) to obtain Assuming non-degeneracy of the leading eigenvalue, we compare this equation to (D4) and obtain: The fluctuation relation [Eq. (14a) in the main text] follows directly from the shift k → k − 1: The fluctuation relation for I(s), Eq. (14b) in the main text, is obtained from that of Λ(k) by Legendre transform: where we choose a N = 1. Substituting (F1) into (D4), we obtain a recurrence relation for a n which terminates at N = k, for positive integer k. In order to satisfy the boundary conditions that r k (z) must to be symmetric for large argument, we need to restrict N to positive even integer k. The recurrence relation (F1) can be written as an eigenvalue problem for the vector a = (a 0 , a 1 , . . . , a k−1 , 1): with the (k + 1) × (k + 1)-dimensional matrix L given by the matrix: The matrix L is of the Metzler type, which means that all its off-diagonal entries are non-negative. For this kind of matrices one can prove that its largest eigenvalue (and thus Λ(k) for positive, even integer k), is strictly real [58]. Using the fluctuation relation given in Eq. (D7) and stated in the main text, we can extend this result to negative even integer k. We conclude that Λ(k) is real for all even integer k. For all finite k, the approach described here can be used to write Λ(k) as the dominant root of a polynomial of order k + 1.
Solving numerically for the largest root gives exact expression for Λ(k), for even integer k. Note that Λ(k) is dimensionless here. To reobtain the dimensional Λ(k) we need to substitute Λ(k) → Λ(k)γ −1 in the solutions of Eqs. (F4).