The statistical geometry of material loops in turbulence

Material elements – which are lines, surfaces, or volumes behaving as passive, non-diffusive markers – provide an inherently geometric window into the intricate dynamics of chaotic flows. Their stretching and folding dynamics has immediate implications for mixing in the oceans or the atmosphere, as well as the emergence of self-sustained dynamos in astrophysical settings. Here, we uncover robust statistical properties of an ensemble of material loops in a turbulent environment. Our approach combines high-resolution direct numerical simulations of Navier-Stokes turbulence, stochastic models, and dynamical systems techniques to reveal predictable, universal features of these complex objects. We show that the loop curvature statistics become stationary through a dynamical formation process of high-curvature folds, leading to distributions with power-law tails whose exponents are determined by the large-deviations statistics of finite-time Lyapunov exponents of the flow. This prediction applies to advected material lines in a broad range of chaotic flows. To complement this dynamical picture, we confirm our theory in the analytically tractable Kraichnan model with an exact Fokker-Planck approach.

We carried out our analysis for two additional direct numerical simulations (DNS) of the Navier-Stokes equation, one at Taylor-scale Reynolds number R λ ≈ 147 and one at R λ ≈ 334. The simulation details are summarized in Table S1. The curvature statistics of an ensemble of material loops were determined in each simulation as described in Methods. The curvature distributions of the supplementary simulations are shown in Fig. S1. Remarkably, the large-curvature power-law exponents, determined by best fits, are almost the same across the simulations, indicating that there may be no significant Reynolds-number dependence in this range. This becomes even more apparent in Fig. S2a, where the curvature distributions are shown to almost perfectly collapse for the latest point in time in the three simulations when nondimensionalized by η. As a function of the integral length, the distributions shift toward larger κ with increasing Reynolds number (Fig. S2b). Also the curvature peak PDFs in Figure S3 show no measurable Reynolds-number dependence. Consistent with our theory, their high-curvature exponent differs from the curvature PDF exponent by 1, mostly within error bars.
The curvature peak number above different thresholds as a function of time is shown in Figure S4 for the two Supplementary Table S1. Main DNS Parameters. Our simulations are run on three-dimensional periodic domains of side length 2π discretized on a real-space grid with N 3 points. The Kolmogorov length and time scales, η and τη, respectively, are computed from the mean kinetic energy dissipation ε and the kinematic viscosity ν. Based on the largest wavenumber kmax resolved by our code, we compute the resolution criterion kmaxη. Using the root-mean-squared velocity component u 2 1/2 and the energy spectrum E(k), we define the integral length L = π 2 u 2 dk k E(k). The integral time scale is computed as T = L u 2 −1/2 . Although loops and FTLE simulations are initialized with identical fields and parameters, the flows eventually diverge due to numerical rounding errors and chaos. The loops simulations contain nL material loops, initially sampled by 300 tracer particles per loop whose number increases roughly exponentially due to refinement in intervals of τ ref . For the FTLE simulations, n tracer trajectories are integrated along with the flow field.
Supplementary Figure S1. Curvature PDFs in supplementary simulations. Both in the low-Reynolds-number (a) and the high-Reynolds-number (b) supplementary simulation, they strongly resemble the distribution of the main simulation. Their power-law exponents, determined by best fits, are almost identical across all simulations.   Figure S2. Comparison of curvature PDFs at the latest simulation times, scaled by the Kolmogorov length (a) and the integral length (b). In Kolmogorov units, hardly any trend is visible whereas the peak of the PDF shifts toward larger curvature values with increasing Reynolds number when nondimensionalized in integral units.
Supplementary Figure S3. Curvature peak statistics in the low-Reynolds-number (a) and the high-Reynolds-number (b) supplementary simulation look the same as in the main simulation. Like for the curvature PDF, the large-curvature power-law exponents, determined by best fits, are almost identical across all simulations.
Supplementary Figure S4. Mean curvature peak number as a function of time in the low-Reynolds-number (a) and the high-Reynolds-number (b) supplementary simulation. The mean number of peaks above different thresholds grows exponentially, proportional to the mean length of the loops. Lines are vertically shifted to compare their growth rate, which is why they are not necessarily ordered as a function of the threshold condition. Best fits to the last third of the total peak number curve yield β = (0.19580 ± 0.00028)/τη (a) and β = (0.26207 ± 0.00014)/τη (b). As in the main simulation, the standard error from the fit is so small that we neglect it in the following. Insets: Curvature peak distribution at the latest simulation time indicating the different thresholds.
supplementary simulations. Notably, our observation that the curvature peak number grows proportionally to the mean length of the loops carries over to these Reynolds numbers. The corresponding growth rate β in units of the Kolmogorov time appears to increase as a function of Reynolds number. Finally, we also determine FTLE statistics by integrating the deformation tensor along trajectories of randomly distributed particles. This is done for 25 million tracer particles in an additional simulation at the smaller Reynolds number and for one billion tracer particles in an additional simulation at the larger Reynolds number that use the same initial condition as the corresponding loops simulations. In Figure S5, we determine the steepest-descent minima needed for our theoretical prediction. Notice that the minimum is taken at values of ρ p that increase with Reynolds number and therefore reach further into the tail of the FTLE distribution. Hence in the large-Reynoldsnumber simulation, despite the enormous number of tracer particles tracked, the minimum can only be resolved up to t ≈ 25τ η . As in the main simulation, we determine the asymptotic value of the minimum by fitting the algebraic decay function (24) to those data points with t ≥ t min ( Figure S5, insets). The time t min is chosen based on the weighted mean squared error, as explained in Figure S6, in order to filter out a transient regime of the decay. The resulting exponents are α = A = 0.61 ± 0.08 for the low-Reynolds-number simulation, α = A = 0.54 ± 0.11 for the main simulation and α = A = 0.55 ± 0.07 for the high-Reynolds-number simulation, all of them consistent with the measured curvature PDF power-law exponents.

SUPPLEMENTARY NOTE 2. NUMERICAL ANALYSIS OF GENERALIZED LYAPUNOV EXPONENTS
As a complementary approach to computing Cramér functions, we may also use generalized Lyapunov exponents (GLE) in order to determine the exponent α based on the implicit equation (47), in short: is the first standard GLE as opposed to the curvature-peak GLE. For their numerical computation, we adopt the method from ref. [1]. We first compute the cumulant-generating function of ρ 1 (t)t, given by log exp(qρ 1 (t)t) , as a function of q and t. In order to estimate Supplementary Figure S6. Different fit choices for extrapolating the steepest-descent minimum. In order to compute the steepest-descent minimum (15), in principle we need a fully resolved Cramér function. In practice, however, we only have finite-time estimates of the Cramér function, which yield finite-time estimates of the minimum (our "data points" in the insets of Figs. 5 and S5). We use the simple decay function (24) to capture and extrapolate the evolution of these data points as a function of the time at which the finite-time Cramér function is computed. In order to obtain a satisfactory fit, it is helpful to leave out a transient regime of data points t < tmin such that only the asymptotic behavior is captured by the fit. Here we show the resulting values of α = A (top) and the decay exponent C (middle) for different choices of tmin in all of the simulations (left to right). The time scale B is not shown. We justify the choice of tmin by computing the weighted mean squared error (MSE, bottom) given by N i=1 (δi/σi) 2 /(N − 3). Here, N is the number of data points included, δi is the deviation of the fit from the i-th data point, and σi is the error of the i-th data point, given by the maximum of the two-sided error computed from the error envelopes in Figs. 5 and S5. The quantity N i=1 (δi/σi) 2 is minimized by the fit, which we divide by the number of degrees of freedom N − 3 (number of data points minus number of fit parameters) for comparability. We choose tmin to be the start of the first plateau of the weighted MSE (red lines). L 1 (q), we perform an affine fit of the cumulant-generating function in the range t ∈ [t max /2, t max ], as exemplified in Fig. S7a for t max = 40τ η . The slope of each fit including its standard error becomes our estimate of L 1 (q). The same procedure is applied to L p (q).
For the main simulation, the results are shown in Fig. S7b. We can first read off the value of β by evaluating L 1 (q) at q = 1. Indeed, the different estimates for L 1 (1) appear to converge toward the value of β previously estimated by other means. In order to estimate α, we need to read off the intersection of the β-line with L p (q). For this curvature-peak GLE, we observe stronger fluctuations as a function of t max . For small t max , we expect the estimates of the cumulant-generating function to be accurate. However, if t max is too small, we have not yet reached convergence of the t → ∞ limit in the GLE. For larger t max , we improve on the convergence of the GLE, but we also rely more heavily on extreme values of ρ p (t) (especially for large q), which are limited by our sample size. Hence we expect the best estimate to be found at intermediate t max . For the main simulation, we indeed find those intermediate curves to come closest to the value of α estimated from the curvature PDF. For the supplementary simulations (Fig. S8), the same argumentation holds and the estimates fluctuate as a function of t max to a certain extent. If we simply read off α from the intersection of lines, the GLE method slightly overestimates α for all simulations, possibly due to both sampling and time-convergence limitations.
Supplementary Figure S7. Generalized Lyapunov exponents in the main simulation. a Plotting the cumulant-generating function of ρ1(t)t for fixed argument q as a function of time (solid lines), the GLE L1(q) can be estimated as the asymptotic slope of the curve by an affine fit (dashed lines) on the interval t ∈ [tmax/2, tmax] where in this case tmax = 40τη. b Generalized Lyapunov exponents can be used to estimate β and α. The first standard GLE L1(q) (solid lines) is shown for tmax ranging from 12.08τη (yellow) to 39.68τη (violet). The curvature-peak GLE Lp(q) is shown for the same times (dashed lines). For comparison, we also show the line q = 1 (solid, black), the value of β estimated from curvature peak number (dotted, grey) and the value of α estimated from the curvature PDF (solid, red).
Supplementary Figure S8. Generalized Lyapunov exponents in the supplementary simulations, at low Reynolds number (a) and at high Reynolds number (b). As in the main simulation, the first standard GLE L1(1) converges toward the previously estimated value of β, but the convergence is slower in the high-Reynolds-number simulation. Similarly, the value of α is slightly overestimated in both simulations. In a, tmax ranges from 12.10τη (yellow) to 39.94τη (violet). In b, tmax ranges from 16.27τη (yellow) to 35.26τη (violet). Note that the value of β as a function of the Kolmogorov time scale τη differs from the one in Fig. S4 because τη is slightly different in the FTLE simulation.

SUPPLEMENTARY NOTE 3. MATERIAL LINE BUNDLES AND FLUX CANCELLATIONS
In the discussion, we explain how our results may help to shed light on the problem of flux cancellations in magnetohydrodynamics (MHD). Flux cancellations occur when magnetic field lines with opposite orientation are brought closely together by the flow and thus cancel each other in the integration of magnetic flux. In our simulations, we observe that the material lines are brought into such a configuration quite frequently (see Fig. S9). They even tend to form bundles of lines where about half of the lines has opposite orientation. Remarkably, these bundles appear to be strongly related to the sharp folds that lead to curvature peaks since the folds are typically observed in the middle or at the end of such bundles. Therefore, it may be worthwhile studying folds of the magnetic field in MHD (possibly detected by curvature peaks) and their relation to flux cancellations.
Supplementary Figure S9. Loops form pairs and bundles of almost parallel material lines, which are closely associated with sharp folds. The loop snapshot is identical to Fig. 1 of the manuscript, taken at t = 27τη in the main simulation at R λ ≈ 216.

SUPPLEMENTARY NOTE 4. SPATIAL AND TEMPORAL RESOLUTION OF LOOP TRACKING
For tracking material loops in our simulations, we treated the sample points of the loops as Lagrangian tracers. Over time, new particles are inserted to sustain the necessary loop resolution. Since using higher-order time-stepping methods for the particles would require histories that are not available for newly inserted particles, we resorted to first-order Euler time stepping. In order to ensure the quality of our results, we tested the code with different temporal resolutions of flow and particle time stepping and different spatial resolution conditions of the loops. The results for simulations at R λ ≈ 146 are shown in Fig. S10 where we compare curvature statistics for the different resolutions. We observe that neither an improved temporal resolution nor an improved spatial resolution of the loops leads to noticeable changes in the curvature distribution.  Figure S10. Curvature distributions for different temporal and spatial resolutions of the loops, taken from 1000 loops after 7.49τη in simulations at R λ ≈ 146. ∆t is the time-step size of the simulation, including field and particle time stepping. ∆s denotes the maximum distance of sampling points of the loops enforced by the refinement. We adjusted the additional refinement condition based on curvature accordingly, effectively doubling the density of sample points uniformly along the loops when moving from ∆s = 0.1η to ∆s = 0.05η and further to ∆s = 0.025η. Other parameters are identical with the low-Reynolds-number simulation shown in Supplementary Note 1, which uses ∆t = 0.018τη and ∆s = 0.1η. Note that the curved shape of the distribution can be attributed to the early time (t = 7.49τη) in the loops' evolution.

SUPPLEMENTARY NOTE 5. GEOMETRIC EVOLUTION EQUATIONS
Here, we derive the evolution equations (54)-(57). Implicitly, all quantities considered in these equations are evaluated along a material line element L = L(φ, t), whose evolution is given by the tracer equation (1). Therefore which is (54). The evolution equation (55) for the tangent vector of the Frenet-Serret frame [2], can then be directly computed The normal vector of the Frenet-Serret frame is defined aŝ where s denotes an arc-length parameterization of the line, i.e. |∂ s L| = 1. Since the transform from s to φ is timedependent, evaluating ∂ t at constant φ and at constant s are different operations. Here, we always want to take ∂ t at constant φ (i.e. at the same tracer particle). Then, ∂ s and ∂ t do not commute. Having this in mind, we compute (summation over repeated indices implied) In order to swap the t-and s-derivatives, we notice that ∂ s = dφ We then insert the evolution equation (S5) for the tangent vectort, where in the last step we used (S6) and ∂ s u k =t m ∂ m u k . Since |∂ st | =κ and δ ij =t itj +n inj +b ibj , we have where we also used orthonormality of the Frenet-Serret frame. This is (56). Finally, in order to derive the curvature evolution equation, we use a simple definition as a function of the tangent vector with arc-length parameterization,κ which is equivalent to definition (2). Using the previous results, we obtain This is (57), which can also be found in ref. [3].

SUPPLEMENTARY NOTE 6. FOKKER-PLANCK EQUATION IN THE KRAICHNAN MODEL
In Methods, we laid out the terms that need to be calculated in order to arrive at the Fokker-Planck equation in the Kraichnan model. In order to proceed, we combine (53) with the evolution equations derived in the previous section and get We want to evaluate these averages using the Gaussian integration by parts formula [4][5][6] combined with the correlation tensor (51). This works analogously for all of them. So let us focus on one of the averages. By introducing delta functions, we can consider the velocity field at the Eulerian coordinate x and take the derivative out of the average (the other quantities are still evaluated at L(φ, t)), Then Gaussian integration by parts yields (S16) The product rule for the functional derivative yields five different terms, which can all be treated in the same way. Let us again focus on a single one of them, namely . (S17) In order to determine the response function δκ δu k (z,t) , we formally integrate the curvature evolution equation (S12), . (S18) By causality, the initial condition will not depend on u k (z, t) for t > 0. The integrand will not depend on u k (z, t) for all t < t either, and the only contribution to the functional derivative can come from the time t = t. Sinceκ,n, andt are integrated quantities of the delta-correlated field u, we expect them to be continuous in time, just like a Wiener process is continuous while its differential is not. Hence their response functions will only be finite, thus not contribute to the integral. Using that the response function becomes where the factor 1 2 comes from the fact that only half of the delta function δ(t−t ) is contained in the integration range [0, t]. Using integration by parts and the sifting property of the delta function, our term M can thus be simplified to Although the spatial correlation tensor R ik (r) can be freely chosen, we can restrict its functional form by assuming isotropy and incompressibility of the Kraichnan field. By isotropy the correlation tensor must be even, hence odd derivatives vanish at zero, e.g. ∂ j ∂ n ∂ o R ik (0) = 0. For the second derivatives, we know that they must have the general form of an isotropic rank-4 tensor [7], By definition, this tensor must be symmetric under exchange of j and n, which implies B = C. Finally, incompressibility implies δ ij Q ik jn = 0 so that leading to the general form [8] Analogous arguments can be used to deduce the general form (59) of the fourth-order derivatives of the correlation tensor. Using this expression, our term M can be evaluated, by orthonormality of the Frenet-Serret frame.
In the following, we list all the terms that need to be computed along with the results of their evaluation. By Gaussian integration by parts and the product rule for functional derivatives, the averages of (S14) split into and Evaluating each of these terms as explained previously yields where Q and P are defined through (58) and (59), respectively. Inserting these results into (S14) yields the Fokker-Planck equation (20).

SUPPLEMENTARY NOTE 7. NUMERICAL RESULTS IN THE KRAICHNAN MODEL
Here, we present a numerical analysis of curvature statistics in the Kraichnan model. To this end, we interpret the tracer equation (1) as a Langevin equation. Since Itô and Stratonovich interpretations coincide for this equation, we may use the Euler-Maruyama scheme [9] to integrate particle trajectories. In every time step, the Gaussian flow field is computed on 1024 3 grid points with a model energy spectrum as described by Pope [10, p. 232], with η a viscous length scale, L ≈ 946η an integral length scale and the functions which determine the large-and small-scale behavior. The spectrum (S44) integrates to the total energy E = 1.05 (code units). The temporal resolution is ∆t = 7.15 × 10 −7 (code units) and the spatial resolution can be quantified by k max η ≈ 2.0, where k max is the maximum resolved wavenumber. Furthermore, we choose β = 5.2, c L = 6.03 and c η = 0.40. For particle time stepping, the field is interpolated using spline interpolation with continuous derivatives up to and including third order computed over a kernel of 12 3 grid points. The loops are adaptively refined as described for the Navier-Stokes simulations in Methods. Figure S11 shows a visualization of an initially circular loop deformed by the Kraichnan field for 1.2Q −1 . Visually, it shares many features of material loops in Navier-Stokes turbulence but appears slightly more compact (compare Fig. 1). The geometric similarities also manifest in the curvature statistics ( Figure S12a), which display the same type of unimodal distribution with power-law tails. Over time, the PDF converges to the stationary solution (21) of the Fokker-Planck equation, featuring the power-law tails κ 1 and κ −18/7 .
The Kraichnan model also forms curvature peaks, whose distribution ( Figure S12b) qualitatively resembles the one in Navier-Stokes turbulence. Our theory predicts the high-curvature tail to scale as a power law with exponent −11/7, which is confirmed by the simulation. In order to form our theory, we made the empirical observation that the curvature peak number grows proportional to the mean line length. This is also what we observe in the Kraichnan model ( Figure S13), where the mean line length can be computed analytically to be proportional to e 4Qt [11].  Figure S13. Mean number of curvature peaks in the Kraichnan model above different thresholds over time, vertically shifted for comparison. Consistent with our observation in Navier-Stokes turbulence, the curves become asymptotically proportional to the mean arc length of loops, which grows as e 4Qt as predicted by theory. Inset: Curvature peak distribution at t = 1.2Q −1 indicating the different thresholds.

SUPPLEMENTARY NOTE 8. LYAPUNOV EXPONENTS BY QR-DECOMPOSITION
In Methods, we have defined finite-time Lyapunov exponents as the growth rate of singular values of the deformation tensor (32). In practice, however, we instead compute the QR-decomposition of the deformation tensor for numerical stability, where Q(t) is orthogonal and R(t) is upper triangular. Note that the matrix Q(t) is unrelated to the constant Q introduced previously, which quantifies velocity gradient fluctuations in the Kraichnan model. The growth rate of the diagonal elements of R(t) can then be interpreted as an alternative definition of FTLEs [1], Note that this definition of FTLEs depends on the choice of the coordinate system. Complementing the Methods section, we here show that peak curvature dynamics of a parabola are exactly captured by this alternative definition of FTLEs. We start from a parabolic material line as defined in (25). Given that the flow is statistically isotropic, FTLE statistics should not depend on the choice of the coordinate system. We can therefore assume that the line is aligned with the coordinate axes, i.e. k(0) = e 1 and l(0) = e 2 . In this case we have where we have used the fact that Q is orthogonal and R is upper triangular. Furthermore, we get Hence, the peak curvature given by (30) can be written as Inserting the alternative definition of FTLEs (S48) yields κ p (t) = e [ρ 1 (t)−2ρ 2 (t)]t κ p (0), an exact analog to the approximate equation (37). This means that the FTLEs defined by the QR-decomposition precisely capture the curvature growth of parabolic line elements. In that sense, the QR-definition of FTLEs is very suitable for our purposes.