Uncovering the multifractality of Lagrangian pair dispersion in shock-dominated turbulence

Lagrangian pair dispersion provides insights into mixing in turbulent flows. By direct numerical simulations (DNS) we show that the statistics of pair dispersion in the randomly forced two-dimensional Burgers equation, which is a typical model of shock-dominated turbulence, is very different from its incompressible counterpart because Lagrangian particles get trapped in shocks. We develop a heuristic theoretical framework that accounts for this -- a generalization of the multifractal model -- whose prediction of the scaling of Lagrangian exit times agrees well with our DNS.

The statistics of the relative pair dispersion of Lagrangian particles (tracers) in turbulent flows has numerous physical applications [1][2][3][4][5][6][7][8][9].Richardson's law, a pioneering work in this field, states that, in a turbulent fluid, the mean-squared displacement between a pair of tracers, R 2 (t) , at time t, scales as R 2 (t) ∼ t 3 [10].In other words, there is a dynamic exponent z = 3/2.Richardson's law can be viewed as a consequence of Kolmogorov's 1941 (K41) scaling theory of turbulence [see, e.g., Ref. 9].It is now well established that the K41 theory is incomplete because it does not account for intermittency, which arises predominantly because of the most dissipative structures in the flow and leads to multifractality [11].In brief, intermittency and multifractality in homogeneous and isotropic turbulence is characterised by the non-trivial scaling properties of moments of velocity differences, ⟨δv(ℓ) p ⟩, over length scales, ℓ, i.e., ⟨δv(ℓ) p ⟩ ∼ ℓ ζp , with the multiscaling exponent ζ p a nonlinear function of p (K41 yields simple scaling with ζ p = p/3) [12].Therefore, it is natural to hypothesize that Richardson's law must have multiscaling corrections too.In particular, we expect dynamic multiscaling [13][14][15][16][17], characterised by not one but an infinity of dynamic exponents; specifically, different moments of R(t) should scale with different powers of t.However, even the most recent experimental [18][19][20][21][22] and numerical [20,[23][24][25][26][27] measurements of R 2 provide only limited support to Richardson's law, because the power-law behavior is observed over a range of scales that covers a decade or so.Hence, there seems to be little hope of extracting the dynamic-multiscaling behavior from experimental or numerical measurements of moments of R(t).
It turns out that the Lagrangian exit time, conventionally defined as the time taken for the separation between a pair of tracers (Lagrangian interval) to exceed a given threshold (e.g., the doubling time) [28][29][30][31], provides a robust measure for dynamic multiscaling.From the statistics of these exit times it is possible to extract dynamic multiscaling exponents that match the predictions from the multifractal model of turbulence [32][33][34].
We have, so far, outlined Lagrangian pair dispersion in incompressible turbulence.We turn now to shockdominated turbulence, which refers to highly compressible turbulence wherein the total energy in the irrotational modes is significantly larger than that in the solenoidal modes [35]; this poses formidable challenges that we must confront because such flows are widely prevalent in many astrophysical systems [36].Over the last two decades, several groups have studied the multifractal properties of compressible turbulence, ranging from flows that are weakly compressible to those that are shock-dominated [37][38][39][40][41][42][43][44].Additionally, the dynamics of tracers and heavy inertial particles in such systems have also been investigated [45][46][47][48][49].However, the possible multifractal generalization of Richardson's law and exit time statistics to compressible turbulence has not been established yet.
We develop the theoretical framework that is required for this generalization to shock-dominated turbulence.It rests on the crucial observation that tracers strongly cluster at the shocks, which comprise the most dissipative structures, making Lagrangian pair dispersion in such flows qualitatively different from that in incompressible turbulence.Consequently, we define two exit times -(i) the doubling time and the (ii) halving time -as the times taken for a Lagrangian interval to (a) go beyond twice and (b) shrink below half its initial length, respectively.We carry out explicit direct numerical simulations (DNSs) of a simplified model of shock-dominated turbulence, namely, the randomly-forced two-dimensional (2D) Burgers equation, which is not only rich enough to display the complexities of such turbulence [37,50,51], but also simple enough for us to develop a heuristic theory for our DNS results.We find that the statistics of the doubling times agree with those that follow from the conventional application of the multifractal model, but those of halving times do not.We generalize this framework to account for the clustering of tracers on shocks and obtain therefrom the scaling exponents for the moments of the distribution of halving times.The results from our heuristic theory are in excellent agreement with those from our DNS.In a similar vein, we define and investigate the scaling properties of doubling and halving frequencies of Lagrangian intervals, and try to understand them on the basis of the underlying flow intermittency.The randomly-forced 2D Burgers equation is: u(x, t) is the Eulerian velocity at position x and time t, and ν is the coefficient of viscosity.f (x, t) is a zero-mean, irrotational, Gaussian, white-in-time random force whose Fourier components f (k, t) [with k the wave vector] obey the constraint given in (1).The equal-time and dynamicscaling properties of the Burgers equation, which can be mapped to the KPZ equation [52], have been studied extensively, but, most often, in one dimension (1D) through DNSs and renormalization-group methods [53][54][55][56][57][58][59][60][61][62][63][64][65] We perform high-resolution pseudo-spectral DNSs of (1) on a bi-periodic square domain with side L = 2π.After the DNS reaches a non-equilibrium-statistically-stationary state (NESS), we calculate the energy spectrum E(k), which shows a power-law regime that is consistent with where δu ∥ (r) Here the symbol ⟨•⟩ denotes averaging over the NESS.Equation (1) exhibits bi-scaling, given the type of forcing we use [60,66], i.e., the scaling exponents Our DNS results largely agree with this [67].In the NESS, we seed the flow uniformly with tracers and track For a Lagrangian interval of initial length R 0 , within the inertial range, a naive application of the multifractal model of turbulence, as done previously [28][29][30], yields the following scaling exponents of the order-p moments of the doubling times, τ D : where The second of these equations gives the bridge relation between the dynamic exponent, κ D p , and the equal-time exponent, ζ −p .We probe the validity of this bridge relation only for p < 1, because the structure functions of negative order, S −p (r), exist only in this regime [68,69].We calculate κ D p and ζ −p , from our DNS and plot them versus p in Fig. 2. Clearly (4b) holds within error bars.By construction, the multifractal model does not distinguish between the scaling behaviors of doubling times and halving times, τ H .This may lead us to expect, for compressible turbulence, that the moments of τ H should have the same scaling as the moments of τ D , i.e., ⟨τ p Our DNS, see Fig. (2), shows that this naive expectation is false, because κ D p ̸ = κ H p .We must, therefore, generalize the multifractal model in order to construct a theory that yields the pdependence of κ H p shown in Fig. (2).We first outline the standard argument that is used to understand Richardson's law [9].A Lagrangian interval of length, R(t), un- Dynamic scaling exponents, of order p, as a function of p for the doubling times κ D p (triangles) and halving times κ H p (squares) of Lagrangian intervals.Clearly, contrary to the naive expectation based on the conventional multifractal model, κ D p ̸ = κ H p .The region shaded in light pink is the multifractal prediction (4), where the equal-time scaling exponents ζ−p and their error bars are calculated from our DNS.The dashed line is our theoretical prediction (8), which agrees very well with the numerical results.
dergoes a Brownian motion (we restrict ourselves to 2D) with a diffusivity, K(R), that depends on R itself, i.e., the probability distribution function (PDF) of R, namely, W (R, t), satisfies The R-dependence of K(R) is deduced from the following dimensional arguments: K ∼ (δ R V ) 2 t cor , where δ R V is the typical velocity difference across the scale, R, and t cor is the typical correlation time.If we now use the K41 forms, δ R V ∼ R 1/3 and t cor ∼ R/δ R V ∼ R 2/3 , we get K ∼ R 4/3 which, when substituted in (5), yields Richardson's law [70].
To carry out a similar calculation for 2D Burgers turbulence, we need an appropriate scaling form of K(R).The first, straightforward, generalization is to replace the where h is the scaling exponent of the velocity field.The second key idea is to recognize that, when we consider a Lagrangian interval of length R here, we must distinguish between the following three possibilities for this interval: the typical correlation time is determined by sweeping as both ends of the interval are trapped in the same shock, so t cor ∼ R, whence K ∼ R 2h+1 .

Case (3) [interval is away from shocks]:
The interval can decrease in only one of the following two ways: (A) Case (3A): Both ends of the interval get trapped, at a later time, in the same shock, so the arguments used in case (1) apply.Therefore, at late times, K ∼ R 2h+1 .(B) Case (3B): The two ends of the interval get trapped in two different shocks, which then approach each other.Hence, δ R V is then the velocity difference between the two shocks, which does not depend on R. Since t cor ∼ R, this yields K ∼ R at late times.
In all of these cases, the calculation of the PDF of halving times is a first-passage problem for (5), with two boundaries, one at R → ∞ and the other at R = R 0 /2, where R 0 ≡ R(t = 0).Given the forms of K(R) in Cases ( 1)-( 3), it is sufficient to calculate the halving time PDFs, P 1 and P 2 , for K ∼ R 2h+1 and K ∼ R, respectively.We obtain [71] (6a) where A 1 and A 2 are numerical constants [see Appendix B for details].By collecting all the cases together, we find the overall PDF of τ H , for large τ H : where w 1 , w 2 , w 3A , and w 3B are the relative weights for the Cases discussed above.As the shocks are onedimensional structures, w 1 is negligible, but the other weights are not; and in the limit R 0 → 0, w 3A P 1 is the dominant contribution.Thereby, we obtain the following result for the moments of the PDF of τ H : In the final step we substitute h = 1/3 [60,66].Equation 8 agrees well with our numerical results [see the dashed line in Fig. (2)].
Several important comments are now in order: (a) We emphasize that, although the halving-time dynamic exponent κ H p depends linearly on the order p, the statistics of halving times is not Gaussian -the tail of the PDF of τ H is exponential [see ( 6) and (7)].Furthermore, (8) does not imply simple dynamic scaling because we have shown explicitly that differently defined time-scales have different scaling exponents.
(b) Although we discuss halving and doubling times, explicitly, our results apply to exit times in general.
(c) Even though our DNS study is limited to the randomly-forced 2D Burgers equation 1, our theoretical arguments can be generalized to other, richer varieties of shock-dominated turbulence in a straightforward manner, as long as the equal-time scaling exponents, ζ p , are known.A critical problem with the multifractal description of the scaling of exit times is that the bridge relation 4, requires the calculation of negative moments of δ R V .One way to avoid the calculation of such moments is to consider the statistics of inverse exit times.In particular, we define the halving and doubling frequencies to be ω H ≡ 1/τ H and ω D ≡ 1/τ D , respectively.Thereby we define two new sets of dynamic scaling exponents, As we have discussed already, the multifractal model for incompressible turbulence does not distinguish between them [28][29][30] clearly, χ H p ̸ = χ D p .Our extension of the multifractal model predicts the correct exponents in the following manner.The scaling of the moments of ω D and ω H , for p > 1, is determined primarily by the small-τ D and smallτ H behaviors of their respective PDFs.The short-time growth of R 0 can occur only if the interval is away from the shocks [Case (3), but at short times].Therefore, for all such intervals, K ∼ R 4/3 and consequently χ D p = 2p/3 for all p [dashed line in Fig. (3)], in agreement with the results of our DNS.As for χ H p , at short times, even for incompressible turbulence, intervals can decrease, hence the the prediction of the multifractal model holds.
In an earlier paper [65] we have investigated dynamic multiscaling in the 1D stochastically forced Burgers equation by using a slightly different formulation than we have used here.Instead of halving times, we have considered interval-collapse times, i.e., the time it takes for a Lagrangian interval to shrink to a point.We have shown that, in 1D, the interval-collapse-time PDF has powerlaw tails.In contrast, we now find that, in 2D, the PDF of halving times has exponential tails.In both 1D and 2D cases, the shocks play a key role in determining the multiscaling properties.
In summary, we have shown how to generalize the multifractal model for incompressible-fluid turbulence to the stochastically forced 2D Burgers equation, which is rich enough to display the complexities of shock-dominated turbulence.Our DNS demonstrates clearly that the statistics of halving times deviate starkly from those of doubling times.By generalizing the standard argument that is used to derive Richardson's law from K41 theory, we have developed a theory that leads to a natural way of understanding the statistics of halving and doubling times.The key idea in this theory is that, when we consider a Lagrangian interval of length R, we must distinguish between Cases (1)-( 3).Our theoretical arguments can be generalized to shock-dominated compressible turbulence, if the equal-time exponents ζ p are known.
We expect shock-dominated turbulence in astrophysical systems at both high Mach and Reynolds numbers, e.g., in the interstellar media and molecular clouds, where the turbulence is driven by supernovae explosions.The dynamics of Lagrangian intervals in such flows provides useful insights into transport and mixing in these systems, which influence chemical kinetics and the rates of formation of stars and planetesimals [36].We expect our generalization of Richardson's law to apply to such systems.For compressible turbulence, where the irrotational and solenoidal components of the flow have similar energies, our theory may require further generalization.Our work brings out the importance of calculating the statistics of both halving and doubling times of Lagrangian intervals from DNSs of compressible turbulent flows, for trans-sonic, super-sonic, and hyper-sonic cases.2. Case (3B): at late times the two ends of the interval get trapped in two different shocks which then approach each other.The velocity difference δ R V is the velocity difference between two shocks; this does not depend on the distance between them, i.e., δ R V is independent of R; furthermore, In all of these cases, the calculation of the PDF of halving times is a first-passage time problem for the following 2D Fokker-Planck equation for the the PDF W (R, t) of R: with absorbing boundaries at R → ∞ and R = R 0 /2 ≡ R(t = 0)/2.It is sufficient to do such a calculation for the two cases with K ∼ R 2h+1 and K ∼ R. We make the change of variable ξ = 1/R, which maps the radial domain to the interval [0, 2/R 0 ].The Fokker-Planck equation in ξ then reads, We first consider the case K ∼ R 1+2h , i.e., K = K 0 ξ −1−2h , where K 0 is an R−independent constant.We define y = ξ (2h−1)/2 and substitute it in (B3) to obtain We solve (B4) with the boundary condition W (y = Y, t) = 0, where Y = (2/R 0 ) (2h−1)/2 , by separation of variables.We obtain where J q is the order-q Bessel function of the first kind, where q = (1 + 2h) (1 − 2h), c n is the n−th zero of J q , and A n is the normalization constant.The initial distribution is W (ξ, t = 0) = W 0 (ξ) = δ(ξ − ξ * )/2πξ, where ξ * = 1/R 0 .In the variable y we have where σ = 1/2 (2h−1)/2 .We set t = 0, multiply both sides of (B5) by y q+1 J q (c m y/Y ), integrate from y = 0 to y = Y σ, and by using the orthogonality properties of Bessel functions, we get for any positive integer m.Thus, The survival probability is given by In the new variable y, the area element ξdξdθ = 2 1−2h y (2h−5)/(1−2h) dydθ, so For large t, the n = 1 term dominates over the others.Hence, by retaining only its contribution and substituting z = y/Y in the integrand, we get By replacing Y by (2/R 0 ) (2h−1)/2 , we can evaluate the halving-time distribution P 1 (τ H , R 0 ) as follows: whence we obtain the following asymptotic form: where A 1 is a numerical constant.
Similarly, for the intervals with K ∼ R, we obtain the asymptotic form of their halving-time distribution: For Cases (3A) and (3B), let τ s be the time at which both ends of the interval get trapped on the same or different shocks, respectively.Then, for τ H ≫ τ s , the halving-time distributions for Cases (3A) and (3B) reduce to P 1 and P 2 , respectively.Therefore, the overall halving-time PDF P(τ H , R 0 ), for an interval of initial size R 0 , is a weighted sum of all these cases [see Eq. ( 7) in the main text].Thus, by retaining only the leading-order contribution in the limit R 0 ≪ 1, we get which yields Eq. ( 8) in the main text.
To validate (B15), we evaluate the halving-time cumulative PDFs Q(τ H ) from our DNS, for different initial interval sizes R 0 , by using the rank-order method, which ensures that the CPDfs are free of binning errors.According to (B15) this CPDF scales as ; this confirms our theoretical result.Furthermore, to verify the assumption that τ H ≫ τ s for intervals belonging to Cases (3A) and (3B), we evaluate the following numerically: 1.The halving-time PDFs P Λ,H , of intervals whose lengths initially exceed a fixed threshold, Λ = 1.5R 0 , before contracting.
2. The PDFs P Λ,s , of τ s for these intervals.
The intervals in Cases (3A) and (3B) lie in the smooth regions of the flow at t = 0, so, at short times, they grow; and, as they do so, their ends fall on shocks which eventually lead to their contraction.Therefore, P Λ,H and P Λ,s are representative PDFs of τ H and τ s , respectively, for these intervals.We define τ * to be the time at which such an interval size exceeds Λ and plot P Λ,H and P Λ,s as functions of τ ′ H = (τ H − τ * )/τ * and τ s ′ = (τ s − τ * )/τ * , respectively, for different values of R, in Fig. (7a).We observe that the tails of P Λ,H are much broader than those of P Λ,s for every value of R 0 .This provides compelling evidence in favor of our assumption that, generally τ H ≫ τ s , for intervals belonging to Cases (3A) and (3B).For the purpose of visualization, we plot some representative time series of sizes, R(t), of these intervals, coloured by the respective means of the instantaneous velocity divergences at their ends in Fig. (7b).
In Fig. 8 we present plots of the dependence of the moments of halving and doubling times on the initial interval size; we also give plots of negative-order equal-time structure functions.

FIG. 1 .
FIG. 1.Typical snapshots of a part of our domain in the non-equilibrium statistically steady state of turbulence.(a) The divergence of the velocity field, ∇ • u; shocks are the dark filament-like structures with large negative divergences; instantaneous tracer positions are marked by red dots.The tracers cluster on the shocks.Surface plots of (b) the x−component and (c) the y− component of the velocity field u(x, t).The shocks are the large negative jumps.
their subsequent motion.We present a pseudo-gray-scale plot of ∇ • u in Fig. (1a), in which shocks are visible as dark filamentary structures.Tracers, shown in red, accumulate on these shocks.In Fig. (1b) and Fig. (1c) we give, respectively, illustrative surface plots of the x− and y− components of u(x, t) at the same instant of time as for Fig. (1a).In these surface plots, the shocks appear as large negative jumps.Henceforth, we non-dimensionalize all quantities by using L and T L ≡ L I /u rms , the largeeddy-turnover time, with u rms the root-mean-square velocity and L I the integral length scale [see Appendix A for details].

5 FIG. 3 .
FIG. 3. Dynamic scaling exponents, of order p as a function of p for the doubling frequencies χ D p (squares) and halving frequencies χ H p (triangles) of Lagrangian intervals.Clearly, contrary to the naive expectation based on the conventional multifractal model, χ D p ̸ = χ H p for all p.The region shaded in light pink is the predicted bridge relation χ H p = p − ζp, where the equal-time scaling exponents ζp and their error bars are calculated from our DNS.The dashed line is our prediction for χ D p .
so the naive expectation is χ H p = χ D p = p − ζ p .On the contrary, our DNS shows (within the error bars of Fig. (3)], χ H p = p − ζ p and χ D p = 2p/3;

FIG. 4 . 5 × 4 E 3 ϵ 1 / 4 the
FIG. 4. (a)Log-log plot versus the wave number k of the energy spectrum E(k) (averaged over shells of radius k in Fourier space and over time in the nonequilibrium statistically steady turbulent state); the inertial range extends over nearly a decade and a half in k, with E(k) ∼ k −5/3 (black line).Inset: Net energy flux Π(k) through the wave-number scale k, displaying a direct cascade of energy; the red straight line indicates Π(k) ∼ A log k, where A is a constant, as suggested on theoretical grounds for the randomly-forced 1D Burgers and 3D randomly forced Navier-Stokes turbulence[64,76].(b) Equal-time structure function exponents, ζp, versus p; the black lines represent bifractal scaling.Inset: Log-log plots of the Eulerian equal-time structure functions, Sp(r), versus r for p = 1 (circles), p = 2 (squares), p = 3 (triangles), p = 4 (triangles) and p = 6 (stars); the shaded region denotes the regime of best power-law fit within the inertial regime; ζp are calculated via local slope analyses of these curves.

FIG. 5 .
FIG. 5. Schematic diagram of a typical case of a Lagrangian interval of length R lying along a shock (solid line).A and B marks the two ends of the interval.

)
In Fig.(6), we plot Q(τ H ), for different values of R 0 in the inertial range, on a semilog graph.All the plots are consistent with exponential tails, which collapse onto a single line (shown dashed in the figure) on rescaling τ H by R 1−2h 0

FIG. 6 .
FIG. 6. Cumulative PDF of halving times The inset shows the semi-log plots of the cumulative PDFs (CPDFs) Q(τH) of halving times.The main figure shows the collapse of the PDFs, after τH is scaled by R 1−2h 0 (see text for our theoretical result).The dashed straight line is drawn for reference.The main figure and the inset use the same colors and symbols.

FIG. 7 .
FIG. 7. (a) Semilog plots of PΛ,s(τ ′ s ) and PΛ,H(τ ′ H ) for different values of R0.This indicates that τH is generally much greater than τs for a typical interval belonging to Cases (3A) or (3B), for any R0.(b) Representative time series of interval sizes, R(t), for R0 = 0.01; colors denote the mean of the values of ∇ • u at the two ends of the interval R(t).

FIG. 8 .
FIG. 8. Log-log plots of some of the (a) positive fractional moments (whose scaling exponents are shown in Fig. 2 in the main text) and (b) negative moments (whose scaling exponents are shown in Fig. 3 in the main text) of τH and τD as functions of the initial interval size R0.(c) Negative-order equal-time functions for different values of p; the scaling exponents ζ−p are also shown in Fig.2in the main text.In all these figures, the vertical lines indicate the boundaries of the regimes of local-slope analyses, which yield the various scaling exponents and their error bars.For all cases, the scaling regimes extend over more than a decade in the independent variable.