Fluctuations of Lyapunov Exponents in homogeneous and isotropic turbulence
aa r X i v : . [ phy s i c s . f l u - dyn ] S e p Fluctuations of Lyapunov exponents in homogeneous and isotropic turbulence
Richard D. J. G. Ho, ∗ Andr´es Arm´ua, † and Arjun Berera ‡ School of Physics and Astronomy, University of Edinburgh, Edinburgh EH9 3FD, United Kingdom (Dated: September 4, 2019)In the context of the analysis of the chaotic properties of homogeneous and isotropic turbulence,direct numerical simulations are used to study the fluctuations of the finite time Lyapunov exponent(FTLE) and its relation to Reynolds number, lattice size and the choice of the steptime usedto compute the Lyapunov exponents. The results show that using the FTLE method producesLyapunov exponents that are remarkably stable under the variation of the steptime and lattice size.Furthermore, it reaches such stability faster than other characteristic quantities such as energy anddissipation rate. These results remain even if the steptime is made arbitrarily small. A discrepancyis also resolved between previous measurements of the dependence on the Reynolds number of theLyapunov exponent. The signal produced by different variables in the steady state is analyzedand the self decorrelation time is used to determine the run time needed in the simulations toobtain proper statistics for each variable. Finally, a brief analysis on MHD flows is also presented,which shows that the Lyapunov exponent is still a robust measure in the simulations, although theLyapunov exponent scaling with Reynolds number is significantly different from that of magneticallyneutral hydrodynamic fluids. ∗ [email protected] † [email protected] ‡ [email protected] I. INTRODUCTION
It is well-known that turbulence implies chaotic behaviour [1, 2]. The description of chaos in turbulent flows hasa wide range of applications ranging from geophysical to engineering flows [3–7]. In particular, the chaotic study ofturbulence is strongly associated with the problem of predictability [3, 4, 8–10]. Chaos theory has also been applied toplasmas and astrophysical flows [11–13]. In the past few decades, there has been an increasing interest in numericallydetermining the degree of chaos in turbulent flows.In recent times, high performance computing has made possible to study chaotic properties of 3D flows usingdirect numerical simulations (DNS) [14–20]. This consists in evolving the Navier-Stokes equations (NSE) using nomodelling. Almost two decades ago, DNS in 2D flows were used to study predictability [21, 22]. Chaotic properties of3D turbulence have been also studied using approximate models such as EDQNM closure approximation [23] or shellmodels [10, 24, 25]. All these works have studied chaos and predictability using the Eulerian description of fluids. Adifferent approach commonly found in literature uses the Lagrangian description of turbulent flows to describe chaos[26–29].Whereas the Lagrangian approach deals with dispersion of pairs of fluid particles, the Eulerian approach measuresthe evolution of two different flows with similar initial conditions. The study of Lagrangian chaos is widely appliedin geophysical [4, 30–36] and astrophysical applications [37–39] as well as the study of engineering flows [40, 41] andplasmas [42–44]. In these contexts, chaos is used to characterize mixing and transport properties as well as determiningLagrangian coherent structures (LCS) [45]. Unlike Eulerian chaos, several statistical analyses of Lagrangian chaoshave been done in the past [26, 27, 46].Eulerian chaos measures the evolution in state space of two initially close realizations of the velocity field u ≡ u ( x , t )and u + δ u , where δ u is infinitesimal. Given the chaotic nature, the difference between both realizations growsexponentially. This can be expressed as | δ u ( t ) | ∼ | δ u (0) | exp ( λt ) . (1)Here λ is the maximal Lyapunov exponent. This exponent is a measure of the level of chaos in a turbulent flow [47].In other words, it sets the timescale in which the evolution of the flow can be predicted.Lagrangian chaos statistics have been widely studied compared to the Eulerian case. For the Eulerian case, themethod used in [15, 16] produces several chaotic realizations that produces different Lyapunov exponents as the flowevolves, this allows to perform a statistical study of Lyapunov exponents. In [16], it is noted that the time evolution ofsuch Lyapunov exponent suffers large and fast fluctuations which differ from the fluctuations usually observed in otherquantities such as energy or Reynolds number. Hence, a deeper study of the fluctuations of the Lyapunov exponentsand their statistical properties is needed. Within a turbulent system there can be two types of fluctuations. One is ofstochastic nature, it occurs due to the interaction with noisy environments such as molecular collisions in the smallscale. The second type is given by the chaotic nature of the evolution in the NSE due to the non-linear term. These arecalled deterministic fluctuations. In this work we focus on describing deterministic fluctuations. Unlike experiments,numerical simulation do not introduce unwanted stochastic effects unless deliberately implemented. Thus, in thiswork, every fluctuation is deterministic.There are a couple of reasons why we study the fluctuations in the Lyapunov exponents in the Eulerian descriptionof a turbulent flow. One is the magnitude and timescales of such fluctuations are significantly different from othercharacteristic quantities such as energy, Reynolds number and dissipation. This shows the need to understand theproperties of the flow that are particular to chaos and their relation to other parameters in the flow. The other ischaracterizing the distributions of flow quantities such as Reynolds number, dissipation rate or Lyapunov exponent isthe key in any statistical description of turbulence. Applications may be possible since the stability of any quantitymay be useful to characterize or calibrate systems.In 1979, Ruelle predicted that this exponent should be proportional to the inverse of the smallest timescale inthe system, that is the Kolmogorov time τ = ( ν/ε ) / , where ε is the dissipation rate and ν the kinematic viscosity[48]. Ruelle’s prediction uses simple arguments based on dimensional analysis. First, Ruelle states that the maximalLyapunov exponent should be associated with the smallest timescale of the system, τ . Subsequent work related λ to the Reynolds number, Re = U L/ν , where U is the rms velocity and L the integral length scale. Here L =(3 π/ E ) R E ( k ) /k dk , E being the turbulent kinetic energy and E ( k ) the energy spectrum. Using Kolmogorov (K41)theory [49], it can be deduced that λ ∝ (1 /T )Re / , where T = L/U is the large eddy turnover time.The relation λ ∼ /τ proposed by Ruelle is based on very simple considerations. Nevertheless, in practice it isobserved that the product λτ is not constant, but it tends to increase slightly with Reynolds number [14, 15]. Usinga multifractal model, and considering intermittency, the following scaling relation is derived [24] λ = DT Re α , (2)where α = (1 − h ) / (1 + h ), and h is the H¨older exponent obtained from the structure function h| u ( x + r ) − u ( x ) |i ∼ V l h , where l = r/L . In Kolmogorov theory, h = 1 / α = 1 / α < .
5, whereas numerical results show values of α > . II. METHOD
In this work, several simulations evolving homogeneous and isotropic flows (HIT) were run using a pseudospectralmethod and a fully-dealiased DNS code in a periodic box with N = 32 , , , and 512 collocation points.Each flow is evolved using the incompressible NSE, ∂ t u = − ( u · ∇ ) u + ν ∇ u + f , (3) ∇ · u = 0 , (4)where u ( x , t ) is the velocity field and f represents an external force per unit volume. A negative damping forcingscheme is implemented. It consists in using the lowest modes in the velocity field ˆ u ( k , t ) (the Fourier transform of u ( x , t )) as a forcing function, ˆ f ( k , t ) = ( ε E f ˆ u ( k , t ) if 0 < | k | ≤ k f , , (5)where E f is the energy of the velocity field contained in the forcing bandwidth. This well tested forcing function[51, 52] has the advantage that the dissipation rate ε is a parameter that can be set a priori in each simulation. Inthis case, ε is set to 0 . k f = 2 .
5, so only the large scales are forced. The field is initializedwith a Gaussian distribution for components of the velocity field including a random seed. A full description of thecode and the forcing can be found in [53].Once the fluid evolution reaches a statistically steady state, the averages of
L, T , U, ε and the finite time Lyapunovexponent λ are computed over multiple large eddy turnovers. However, a certain level of fluctuation is always observed.This is possibly associated with the forcing mechanism together with the irregular spatio-temporal behaviour given bythe highly non-equilibrium nature of the turbulent regime. A complete understanding and characterizations of thesefluctuations are not the aim of this work, but it would be an interesting step to take in the future. Instead, the maingoal of this study is to understand the relations between these variables and also with other simulation parameterssuch as the lattice size N . The finite time Lyapunov exponents are computed for each simulation using the methoddescribed in the following section.As well as looking at the evolution of the NSE, we perform a similar analysis for the case of homogeneous andisotropic MHD flows. This analysis consists in evolving the velocity field u and the magnetic field b using the followingincompressible MHD equations ∂ t u = − ( u · ∇ ) u + ν ∇ u + ( ∇ × b ) × b + f , (6) ∂ t b = ( b · ∇ ) u − ( u · ∇ ) b + η ∇ b , (7)together with the incompressibility condition given by Eq.(4) and ∇ · b = 0 , (8)where η is the magnetic diffusion coefficient. For MHD, the Lyapunov exponents measure the divergence of a combi-nation of both the velocity and magnetic fields as described in the following section. A. Finite Time Lyapunov Exponents method
Finite time Lyapunov exponents are widely studied in many applications within the Lagrangian description [4, 26,27, 45, 54]. According to the theory, the Maximal Lyapunov exponent do not depend on the initial infinitesimalperturbation [55], but a slight dependence is usually present for finite perturbations in finite times [56]. Hence, theMaximal Lyapunov exponent is estimated by averaging over an ensemble of different realizations of finite perturbationsin finite times.In this section, we present the method used in this work to compute the FTLE, which is also used in [15]. Inthe Eulerian description, the procedure used to measure Lyapunov exponents consists in evolving two fields that areinitially similar. Given the original field u , another field u = u + δ u is created. Then, we evolve the system overa finite steptime ∆ t , measure the field difference, and then introduce a successive perturbation, starting the processagain each time. The field difference is defined as δ t = p E ud ( t ) , (9) E ud = 12 Z d k | ˆ u ( k , t ) − ˆ u ( k , t ) | . (10)The perturbation introduced every steptime ∆ t is given by u ( x , ∆ t ) = u ( x , ∆ t ) + δ δ ∆ t δ u ( x , ∆ t ) , (11)where δ is set to 10 − and ∆ t = 0 . λ is measured as ˜ λ = 1∆ t ln (cid:18) δ ∆ t δ (cid:19) , (12)and this procedure is repeated, obtaining a set of values for ˜ λ . The mean Lyapunov exponent and its variance werecalculated in the usual way λ = h ˜ λ i , (13) σ λ = 1 n X h ˜ λ i − h ˜ λ i , (14)where n is the number of FTLEs measured.In the case of MHD, we use the same procedure. The only difference is that δ t is taken as a combination of theincrease in both the velocity and magnetic field difference δ t = p E ud ( t ) + 2 E bd ( t ) , (15)where E bd = 12 Z d k | ˆ b ( k , t ) − ˆ b ( k , t ) | . (16)The advantage of the FTLE method is that it produces a sample of Lyapunov exponents that allows us to performa statistical analysis of the measurements. Other methods exist in literature, for instance, the direct method consistsin introducing a small perturbation once the system reaches a steady state and then letting the two fields evolveuntil they become uncorrelated [14, 17]. This method reveals interesting features because it displays the long timeevolution of the spectrum of the field difference | δ u | (in which a transient and a saturation stage are found [14]).However, this procedure takes a relatively large computational time to get just one Lyapunov exponent. In contrast,the FTLE method is useful to obtain a large number of Lyapunov exponents which can be used to analyze thestatistical distribution for a given steady state. This procedure shows that fluctuations in the Lyapunov exponentsare quite significant and this fact was not previously noted using the direct method. Hence, it is interesting to analyzethe source of such fluctuations and see if they can be explained by other fluctuating quantities in the system. III. HYDRODYNAMIC RESULTS
The values of λ measured using the FTLE method produces a time signal with large fluctuations. These fluctuationsare briefly mentioned in [15] and some aspects are studied in [16, 18]. In our work, we further analyze the propertiesof these signals. We look at the mean and standard deviation of λ and study their dependence with different flowparameters to see if these variations are of physical origin or if they are caused by the choice of parameters in thenumerical simulations.First, we look at the probability distributions of Re and the FTLEs. Second, we look at how the fluctuations in λ depend on the lattice size and Reynolds number. Third, we find a resolution to the discrepancy in results betweenthose of [14] which quoted α = 0 . ± .
03 and [15], which had a higher value of α = 0 . ± .
05. Fourth, we observehow the Lyapunov exponent fluctuations depend on the choice of the steptime parameter ∆ t of the FTLE method.Fifth, we study the chaotic behaviour of decaying turbulence and use the high sensitivity of the system with respectto the Lyapunov exponent to test which is the proper large timescale to describe Eulerian chaos. A. Distributions of λ and Re When the fluid reaches a steady state, its statistical quantities such as Re and total energy, fluctuate around amean value. We first measure the fluctuations in Re, and then observe if these are consistent with the fluctuationsmeasured for λ .In Figure 1 we show the typical distributions obtained for Re in the steady state, measured for roughly 10-20 largeeddy turnover times. In Figure 2 we show the typical distributions obtained for ˜ λ in the steady state, measured forthe same run time as the equivalent distribution for the Reynolds number.Surprisingly, most of the distributions of FTLEs are well approximated by a Gaussian whereas the Reynolds numbershistograms were poorly approximated by a Gaussian distribution and even far from a bell shaped distribution in mostcases. For a much longer measurement time, the distributions for Re and other large scale statistics approximate aGaussian distribution more than those in Figure 1. The distribution of Figure 3 shows the distribution for a singlerun with a measurement time of approximately 500 large eddy turnover times, which corroborated that a Gaussiandistribution is approximated for longer run times. However, even after 10-20 large eddy turnover times, the Restatistics approximate a Gaussian very poorly. This result may be due to the fact that successive values of Re in asimulation are not random but strongly correlated, whereas for the Lyapunov exponent the correlation is weaker.It is remarkable that given Equation (2), the distribution of Lyapunov exponents approximates to a Gaussian formwhen the Reynolds numbers does not even come close to this. Apparently some source of randomization is introducedinto the measurement, which could well be a consequence of the perturbation method, which occurs at a relativelyhigh rate. This is one of the most important observations in this paper.In Figure 4, we present the distribution of λ for a simulation with a run time of 200 T , obtaining n = 11863 valuesfor ˜ λ . From this Figure, we can see that, even though FTLE distributions approximate to a Gaussian, it does nottend to an exact normal distribution for large n and some noise is still observed. The distribution in this Figure hasskewness -0.12 and flatness 2.66, where skewness is the third standardized moment of the probability distributionand flatness is the fourth standardized moment. For a Gaussian distribution we would expect skewness to be 0 andflatness to be 3. However, for this simulation, the skewness of total energy, Re, and dissipation is always positive andin the range 0.42-0.46. Flatness is slightly above 3 for these statistics.We may measure the FTLE immediately after initializing the simulation, or wait until after the system reaches astatistically steady state. If the FTLE is measured immediately after initialization, there is a initial transient until itreaches a statistically steady state, but this transient is shorter than for the other large scale statistics such as totalenergy and Reynolds number. This is also provides evidence that the correlation timescale for the Lyapunov exponentis much shorter than for Reynolds number, energy and dissipation rate. . . . . . . . . Re . . . . . . . P r o b a b ili t y d e n s i t y (a)
340 360 380 400 420 440 Re . . . . . . P r o b a b ili t y d e n s i t y (b)
850 900 950 1000 1050 1100 1150 Re . . . . . P r o b a b ili t y d e n s i t y (c) FIG. 1. Normalized distributions of Re for three simulations: (a) N = 128 and ν = 0 .
01 , (b) N = 256 and ν = 0 . N = 512 and ν = 0 . B. Dependence of λ and Re on lattice size We look at the dependence of the Lyapunov exponent λ and its fluctuations σ λ on the lattice size. The Lyapunovexponent, Reynolds number and their fluctuations are measured for different runs keeping viscosity and dissipationrate fixed ( ν = 0 .
01 and ε ∼ . N = 32 , , , , ∼ − λ on lattice size N for different initial seeds. As can be seen, the Figureshows no correlation or dependence of λ on lattice size. The mean values λ varies little with N compared to thefluctuations represented in the error bars. In Figure 6, we show the dependence of the fluctuations of the Lyapunovexponent σ λ on the lattice size. A slight trend is observed in the Figure, although the variation is not significantcompared to that given by the variation in Re in the following section and the uncertainty in the value of σ λ is notrelevant given that this affects only the second significant figure in every measurement we perform. C. Dependence of λ on Re After determining the variation of the fluctuations with the lattice size, we observe how λ and its fluctuationsvary depending on the Reynolds number. For this, eleven simulations were run varying viscosity in order to obtaindifferent values of Re ranging between 80 and 1200. Table I shows parameter values for all simulations. Once thevalues for λ, σ λ ,Re and σ Re are obtained, a weighed least squared fit is performed as shown in Figure 7 to Equation(2), obtaining a value for α = 0 . ± .
02 and D = 0 . ± . T and λ are studied. Then, using the relation in Eq.(2), the fluctuations of the quantities .
25 0 .
30 0 .
35 0 .
40 0 .
45 0 . ˜ λ P r o b a b ili t y d e n s i t y (a) . . . . . ˜ λ . . . . . . . . . P r o b a b ili t y d e n s i t y (b) . . . . . . . . ˜ λ . . . . . P r o b a b ili t y d e n s i t y (c) FIG. 2. Normalized distributions of ˜ λ for three simulations: (a) N = 128 and ν = 0 .
01, (b) N = 256 and ν = 0 . N = 512 and ν = 0 . on the r.h.s are propagated to estimate the expected level of fluctuations of the Lyapunov exponent on the l.h.s, thuscorroborating if the level of fluctuations is consistent on both sides. It is observed that the level of fluctuations in λ is larger than that expected from Eq.(2). This might indicate that there are other effects on the system introducingthese variations, although it is important to recall that Eq.(2) is not to be taken as a fundamental relation, sinceit is derived from dimensional considerations. According to the theory, FTLEs have some dependence on the initialfinite perturbation and the finite time in which it is measured, so the fluctuations can have origin not only in thefluctuations of the steady state but also by the successive perturbations introduced by the FTLE method and thesteptime ∆ t .Figure 8 plots the fluctuations in Re, σ Re , against their mean values. It shows that there is a linear dependencein the observed range of Re. The following relation is found, σ Re = c Re, with c = 0 . ± .
01. Similarly for therelation between λ and σ λ , a linear dependence is observed and the following relation is measured, σ λ = c λ , with c = 0 . ± .
02. In our simulations, the values for T and σ T are approximately constant, with σ T /T ∼ . σ T and σ Re and usingthe measured values we compare the result to the σ λ obtained using the FTLE procedure. All fluctuations aredeterministic and are measured in the steady state. The following relation is obtained σ λ ≈ (cid:12)(cid:12)(cid:12)(cid:12) ∂λ∂Re (cid:12)(cid:12)(cid:12)(cid:12) σ Re + (cid:12)(cid:12)(cid:12)(cid:12) ∂λ∂T (cid:12)(cid:12)(cid:12)(cid:12) σ T (17)Here, we consider D and α as universal constants with no fluctuations associated to them. Their uncertainties do notstem from deterministic fluctuations as for Re and T , instead they have origin in the linear regression performed toestimate their values. This distinction is highly important for this analysis. To compare fluctuations, it is advantageousto consider the standard deviation of λ relative to its mean value, so we expect σ λ λ ≈ α σ Re Re + σ T T = 0 . .
65 70 75 80 85 90 95 100 Re . . . . . . . . . P r o b a b ili t y d e n s i t y FIG. 3. Histogram for the distribution of Re in a single run with measurement time of approximately 500 T . The solid blackline represents the Gaussian function with mean and variance corresponding to the data sample .
25 0 .
30 0 .
35 0 .
40 0 . ˜ λ P r o b a b ili t y d e n s i t y FIG. 4. Histogram for the distribution of ˜ λ in a single run with N = 128 and ν = 0 .
01, obtaining a large sample of n = 11863.The solid black line represents the Gaussian function with variance and mean given by the numerical data. N . . . . . . . . λ FIG. 5. Dependence of λ on lattice size N , different shapes represent different initial seeds. N . . . . . σ λ T FIG. 6. Dependence of Lyapunov exponent fluctuations σ λ T vs. N , different shapes represent different seeds for the initialrandom field for each lattice size N ν ε Re λ σ λ T T E τ k max η ν is the kinematic viscosity, ǫ is the dissipation rate, Re is the Reynolds number, λ is the maximal Lyapunov exponent and σ λ its standard deviation, T = L/U and T E = E/ε are two different definitions for the large eddy turnover time, τ is the Kolmogorov microscale time, k max ≈ N/ − η = ( ν /ε ) / is the Kolmogorov microscale length and k max η is theresolution. Comparing this expected value to the measured value σ λ λ = 0 .
2, we find that the level of fluctuation we expect fromEq.(2) is inconsistent with the one measured. Such fluctuations appear to have origin not only in the fluctuationsof T and Re. There are many possible explanations for the high level of fluctuations in λ . The perturbation in theFTLE procedure or the forcing could affect the value of σ λ . Another possible cause is that even for infinitesimalperturbations, turbulent flows are not characterized by a unique maximal Lyapunov exponent, contrary to what thetheory indicates. The reason might well be a combination of all those previously mentioned. The important fact isthat σ λ /λ fluctuations are significant and they cannot be simply estimated by measuring the values of Re, σ Re , T and σ T .The relationship between fluctuations of the Lyapunov exponent and the mean Reynolds number is also measured.Figure 9 shows a plot of σ λ against Re with a fit to the scaling relation σ λ T ∝ Re γ . The value found for oursimulations is γ = 0 . ± . D. Effect of definition of T on α We note that previous results which have analyzed the relationship in Equation (2) have found disagreement inthe value of α which they produce. For example, previous results using the direct method for DNS have given an α = 0 . ± .
03 [14], which agrees with the finding using an FTLE method as done here. However, other FTLE DNS0 Re × − × × × λ T FIG. 7. (Color online) λT vs. Re obtained for eleven simulations where T = L/U . The linear fit is represented by a solidgray (red) line and Ruelle’s prediction is represented by the dotted gray line.
200 400 600 800 1000 1200 Re σ R e DataLinear fit
FIG. 8. σ Re vs. Re for statistically steady state of the flow and linear fit for slope c = 0 . ± . Re × − × − × − × − σ λ T FIG. 9. Data and linear fit to determine power law σ λ T ∼ Re γ results have found 0 . ± .
05 [15]. This section offers a resolution to this discrepancy in the results.1One of the main features of turbulence is that it involves many characteristic scales. Furthermore, in homogeneousand isotropic turbulence, the absence of boundary conditions sets an extra complication in determining the large scalequantities. Using different definitions of the large length and time scales can result in different analyses. The keydifference between the prior works was that these used different choices for the definition of L and T . In [14], L isthe integral length scale defined in section I and T = L/U , whereas in [15], this value was defined differently, being T E = E/ε and then L E = U T E . The forcing used in [14] was the same used in this work and the same as used in[15].The difference between these definitions deserves some attention. The definition of the integral length scale L stemsfrom the correlation function R ij ( r ) = h u i ( x ) u j ( x + r ) i and can be associated with the size of the largest eddies.Hence, the first definition, given by the integral length scale, can be understood as the average time that it takes fora large eddy initially occupying a given region is space to move to a different region which is completely uncorrelatedwith the initial one. The second definition T E takes the total energy E and the dissipation rate ε , which is associatedwith the characteristic time in which the largest eddies transfer energy to smaller eddies in the energy cascade.Taylor in 1935 proposed the following relation for a HIT flow in a steady state, using dimensional analysis [57], ε = C ε U L , (18)where C ε is known as the dimensionless dissipation rate. Evidence has shown that this coefficient tends to a constantvalue C ε, ∞ ≈ . T and T E in the large Reynolds numbers limitwhere C ε ≈ . E = 3 U /
2, these two definitions differ only by a constant factor, being T E = 3 T . Insuch case, there should be no discrepancy when using these different definitions to test Ruelle’s relation. Nevertheless,if we consider that the dimensionless dissipation has some Reynolds number dependence, then both definitions willhave an impact in the scaling relation between λ and Re given by Eq.(2).We can improve the accuracy of this analysis by considering a useful expression for the dimensionless dissipationrate derived by McComb et al [50]. The following expression is derived using the von K´arm´an-Howarth equationand performing an asymptotic expansion in the second and third-order structure functions in powers of the inverseReynolds number, obtaining C ε ( Re ) = C ε, ∞ + CRe + O (cid:18) Re (cid:19) , (19)where C is a constant that does not depend on Re. Also in [50], these constants are computed using DNS. Thenumerical values obtained are C ε, ∞ = 0 . ± .
006 and C = 18 . ± .
3, which are in agreement with other values inliterature [59]. Our simulations show that the relation of C ε (Re) as a function of Re is in strong agreement with theresults of [50].Considering this result, we see that the relation between T E and T carries some extra dependence on the Reynoldsnumber, T E = Eε = 3 L C ε U = 32 C ε ( Re ) T . (20)The scaling relation in Equation (2) is re-tested using the definition of T E . We expect that when we use the definition T E instead of T , a different exponent α E is found instead of α . So λT E = D E Re α E . (21)Figure 10 shows the data for λ as a function of Re using the definition T E . Using a linear fit, we obtain the valuesfor Equation (21) of α E = 0 . ± .
006 and D E = 0 . ± . T and T E . Starting from Eq.(21), the following relation is obtained λ T E T C ε = D E Re α E ,λT = f ( Re ) ≡ D E Re α E C ε ( Re ) . (22)2 Re λ T E FIG. 10. (Color online) λ vs Re using the definition of T E = E/ε . The linear fit is indicated by a solid gray (red) line whereasRuelle’s prediction is denoted by the dotted line. Re × − × × × λ T FIG. 11. (Color online) Data obtained for λ vs T , the straight solid gray (red) line represent the linear fit, the straight greydotted line represents Ruelle’s prediction and the curved dashed light gray (cyan) line represents the prediction given by f (Re). This relation sets an extra dependence on the Reynolds number that is not present in the original Ruelle’s relationin Eq.(2). Figure 11 shows the comparison between the data obtained, Ruelle’s prediction, and the prediction givenby f (Re) in Eq.(22) by using the two different definitions T and T E . It can be seen that the slight curved deviationfrom the line is well predicted by f (Re).We can also predict the expected value of α when considering the extra Re dependence in Eq.(19) and given themeasured value for α E and see if it is consistent with the one found in the previous section. To compare bothexponents, we consider equations (2) and (22). Then we take the following derivative to determine the scaling relation α = d ln( λT ) d ln Re = α E + Re d ln ( C ( Re )) dRe ,α = α ( Re ) = α E −
11 + C ε, ∞ C Re . (23)As expected from Figure (11), the scaling is not a constant since it has a slight dependence on Re. To compare withthe expected value of α previously obtained, we first evaluate α (Re) using the same values of Re measured in thesimulations. Then we take the average, obtaining h α i exp = α E − *
11 + C ε, ∞ C Re + ≈ . , (24)3which is within one standard deviation of the measured value of α . This shows that both procedures are consistent.Although the data in Figure 10 is more linearly correlated. This resolves the apparent disagreement between resultsin [14] and [15], and is one of the main results in this paper.We also measure the scaling relation σ λ T E ∝ Re γ E and compare it to that found in [15]. In that work, the authorsclaim that the uncertainty for this value is large but obtain γ ∗ ∼ .
2. In our simulations, we find a smaller value γ E = 0 . ± .
05, which is closer to the value of 0 . ± .
05 for γ found previously using the integral length scale todefine T . E. FTLE steptime dependence of λ When the FTLE method is used, it is necessary to set certain parameters such as the magnitude of the perturbationor the steptime ∆ t between successive perturbations. It is reasonable to suppose that for large values of the steptimethe results are more stable since the maximal Lyapunov exponent is defined for t → ∞ . However, in the case ofFTLEs some dependence on the initial perturbation is expected, hence, the drawback is that for finite times, i.e.∆ t << T E , the growth might not reach stability and thus the estimate of λ and σ λ might be inaccurate for reasonswe give in Section IV. On the other hand, having shorter steptime would be an advantage in the sense that it allowsone to obtain a larger sample of Lyapunov exponents in a significantly shorter simulation time. .
00 0 .
05 0 .
10 0 .
15 0 .
20 0 .
25 0 . . . λ .
00 0 .
05 0 .
10 0 .
15 0 .
20 0 .
25 0 . ∆ t . . . σ λ FIG. 12. Dependence of λ and σ λ on the steptime ∆ t . Figure 12 shows the result of varying the steptime ∆ t on the Lyapunov exponent and its fluctuations. The simulationhas Re = 50. As can be seen in the Figure, the values of λ and σ λ are remarkably stable even for very short steptimesusing the FTLE method. Indeed, at the extreme lowest value, the steptime is so short that it is only 5-10 simulationtimesteps. At the other extreme, the direct method, previously described, is akin to having a very high steptimesuch that, within one steptime, the system becomes decorrelated. Even using the direct method, the mean valuesof λ are similar to the ones found using the FTLE method.These results suggest that the steptime is of very littleimportance to the measured value of λ , and so it can be measured with a very low steptime which helps to reduce thecomputational cost. This finding can be useful in certain applications that use the Lyapunov exponent to characterizeturbulent flows [18], especially if we find a way to relate the Lyapunov exponent to other quantities that characterizethe flow such as Reynolds number, energy or dissipation rate. This is the main finding in this paper and togetherwith the fast stability described in Section IV. This means that the Lyapunov exponent is a robust measure in thesetype of simulations. F. Decaying turbulence
We analyze the chaotic behaviour in a decaying turbulent flow. The procedure consists in evolving an initial fielduntil it enters a power law decay phase. After that, a copy of the field is created and perturbed infinitesimally. Then,we measure the evolution of the energy difference in time E d ( t ) and compare it with predictions. Using the relationin Eq.(2) for steady states, we are able to predict the evolution of E d ( t ) by adapting it to a decaying flow. For thiswe consider the difference field at a given simulation time t i as | δ u ( t i ) | ≡ δu i . For each simulation time, we assumethat the difference grows exponentially. Thus we have4 δu i +1 ≈ δu i exp ( λ i +1 δt ) | u i +1 || u i | , (25)where δt ≡ t i +1 − t i is the timestep of the simulation. We multiply by the factor | u i +1 | / | u i | to account for the decayin the magnitude of | u | , that we assume uniform. It is important to remark that this approximation also assumesthat the Lyapunov exponent vary slowly in time, thus ( λ i +1 + λ i ) / ≈ λ i +1 .Continuing for i + 2, we obtain δu i +2 = δu i +1 exp ( λ i +1 δt + λ i +2 δt ) | u i +2 || u i | , (26)and so, in general δu f = δu i exp Z fi dtλ ( t ) ! | u f || u i | , (27)where we replace the infinite sum for the integral over time, and λ ( t ) is a function of time.We use this prediction to discriminate between the different functional dependencies of λ given in Eqs. (2) and (21).This is a convenient comparison method because the predicted value of δu f in Eq.(27) is very sensitive to changes in λ ( t ). Using these laws and the statistical properties of the unperturbed field, we can numerically perform the integralin Eq.(27) without having to measure Lyapunov exponents in the decaying regime. Thus, we predict the behaviourof E d ( t ) given its initial value E d (10). We use a simulation with N = 1024 and we perturb the flow at wavenumber k = 40.Figure 13 shows the influence of the different functional forms of λ using Eq.(27) to predict E d ( t ) as well as themeasured value of E d ( t ) for the above mentioned simulation. As can be seen, the prediction given by (21), that uses T E as the large time scale is the best fit to the measured data. This fit can be improved by using slightly highervalues of D E and α E , which would still be within one standard deviation of the fit for the forced data. This isfurther evidence that the functional form of λ predicted by Eq.(21) is the best fit to data. This has implications forthe physical origin of the Eulerian chaos in turbulence, since we cannot derive Eq.(21) from the prediction λ ∼ /τ consistently with the Kolmogorov theory. As such, chaos in HIT must have some different physical basis. − − − −
10 100 E d ( t ) t FIG. 13. Different predictions of the functional form for λ . The red dotted line corresponds to the prediction in (21) the bluedotted line to the prediction in (2). The solid black line corresponds to the measured data. IV. TIMESCALES IN HIT
Quantities such as energy, Reynolds number, Lyapunov exponents and dissipation rate fluctuate around a meanvalue during a steady state. We look at the timescales of these time signals produced in our simulations. Especially welook at the fluctuations in the Lyapunov exponent and we compare it with those of energy and dissipation rate. One of5the interesting features noted in [16] is that the fluctuations in the Lyapunov exponents are faster than those of energyor dissipation rate. This can occur due to the perturbation we introduce every steptime ∆ t . Here a more systematicstudy of this timescale is presented. We observed the autocorrelation of the signal to measure its self-decorrelationtime, which we will simply call decorrelation time. Given the above findings that histograms for some magnitudesreach a stable distribution faster than for others, we may wonder what is the run time necessary to get a stablemeasurement of different variables depending on the timescale inherent in each time signal. In order to determinesuch time, it is useful to look at the decorrelation time.Previous results show that E ud has a limit to its growth rate which is roughly equal to ε [14]. We also know that thevalue of E ud saturates at 2 E As such, the minimum time that it takes for two fields to become completely decorrelatedis 2
E/ǫ = 2 T E . This result means that the timescale that is relevant for the energy and other large scale propertiesare dependent on T E .For this purpose, we use the autocorrelation in time ρ X, ∆ t for some random variable X . The autocorrelation isrelated to the covariance ρ X,Y , defined as ρ X,Y = E [( X − µ X )( Y − µ Y )] σ X σ Y , (28)where E [ ... ] is the expectation, Y is some random variable, µ is the mean, and σ the standard deviation. If we thendefine Y as the value of X after some time ∆ t (in this section ∆ t is not related in anyway to the steptime ∆ t used inprevious sections), then we can define ρ X, ∆ t as equivalent to ρ X,Y here. In this case, X and Y should have the samestatistics, and so the equation simplifies to ρ X, ∆ t = E [ XX ∆ t ] − µ X V ar [ X ] , (29)where X ∆ t is the value of X after time ∆ t , and V ar [ ... ] is the variance.We perform simulations of forced HIT using a value of ν = 0 .
006 on a lattice of size 128 with ε = 0 .
1. Thesimulations had T E = 5 .
43 in simulation time units, with Re = 130. In this instance, we run it until simulation time1000, which is over 450 large eddy turnover times. − ρ E , ∆ t ∆ t FIG. 14. Solid line shows complementary autocorrelation of energy 1 − ρ E, ∆ t for Re = 130. Dashed line shows scaling with(∆ t ) , dotted line shows constant 1. Inset is same as main plot but with log-log axes. Figure 14 shows (1 − ρ E, ∆ t ) for the simulation described in the previous paragraph. For small ∆ t we can approximate1 − ρ E, ∆ t = (∆ t/T d ) , where we define T d as the characteristic decorrelation time. In this case, we find that T d = 5 . T E . Similarly, we find that the point at which ρ E, ∆ t = 0 is at ∆ t = 11 . T E . Weexpect that the important timescale for measurement may depend on the length scale which is most important forthose statistics. For instance, whilst total energy will be dominated by large scale structures, the dissipation shouldbe dominated by the smallest lengthscales. Figure 15 shows 1 − ρ ε, ∆ t for the same simulation. Approximating1 − ρ ε, ∆ t = (∆ t/T d ) gives T d = 2 .
58, which is approximately half the value of T E . Similarly, the time at which ρ ε, ∆ t < T E .We calculate the autocorrelation for E ( k ) in order to look at the dependence on the length scale explicitly. Theplot for ρ E (1) , ∆ t looks very similar to those for ρ E, ∆ t . However, the plots for k >
1, even as low as k = 2, look more6 − ρ ε , ∆ t ∆ t FIG. 15. Solid line shows complementary autocorrelation of dissipation 1 − ρ ε, ∆ t for Re = 130. Dashed line shows scaling with(∆ t ) , dotted line shows constant 1. Inset is same as main plot but with log-log axes. like the plots for ρ ε, ∆ t , despite the fact that k = 2 is within the forcing range. All of these plots reach ρ E ( k ) , ∆ t ≤ T E when k >
1. For moderately high k , there is a (∆ t/T d ) dependence at low ∆ t . But, when welook at the autocorrelation for very high wave number, something unexpected happens, which is that a new scalingappears. − ρ E ( ) , ∆ t ∆ t FIG. 16. Solid line shows complementary autocorrelation of E (26), 1 − ρ E (26) , ∆ t for Re = 130. Dashed line shows scaling with(∆ t ) , dotted line shows scaling with (∆ t ) / . Figure 16 shows 1 − ρ E (26) , ∆ t . The value k = 26 is chosen because it is roughly the inverse of the Kolmogorov mi-croscale for this simulation. The new scaling is approximately (∆ t ) / . Even with this new scaling, the autocorrelationstill becomes decorrelated before T E .We now look at the autocorrelation for λ itself. Figure 17 shows 1 − ρ λ, ∆ t . These statistics are measured on asimulation with a much shorter run time and with FTLE steptime ∆ t = 0 .
05. We find that T d ≈ . T /
4. This is further evidence that the FTLE should reach stable statistics faster thanfor instance, the Reynolds number or total energy. This is one of the main results in this paper. We may expect thislonger correlation time for the statistics like total energy, Re, and dissipation compared to the Lyapunov exponentjust from looking at their evolution in time. The graphs for Lyapunov exponent fluctuate more wildly. The correlationbetween Lyapunov exponents has a much shorter timescale which is why the statistics of the Lyapunov exponent aremore quickly approximated by a Gaussian. However, the Lyapunov exponent becomes decorrelated even faster thanthe energy at the Kolmogorov microscale, even though we would expect those scales to be evolving the fastest in thewhole system. This might indicate that the perturbations introduced in the FTLE method have a relevant effect onthe Lyapunov exponent self-correlation time. However we cannot consider successive values of λ to be uncorrelatedsince the decorrelation timescale T d is much larger than the steptime.Given the decorrelation time T d , we need to determine the run time T r necessary to obtain proper statistics for7 − ρ λ , ∆ t ∆ t FIG. 17. Solid line shows complementary autocorrelation of λ , 1 − ρ λ, ∆ t for Re = 130. Dashed line shows scaling with (∆ t ) ,dotted line shows constant 1. Inset is same as main plot but with log-log axes. each quantity. It is important to mention that T r is not the actual run time of the simulation but the time in whichthe simulation is in the steady state, so initial transients are not considered.For this purpose, a useful quantity to look at is the standard error of the mean σ X . We may consider a set ofmeasurements of a quantity X , with its own mean X and standard deviation σ . If instead, we consider a subsetof n measurements, we can measure a mean value of X ( X sub ) that will usually differ from the real mean value X .Repeating this procedure many times for different subsets, we can get a distribution of X sub . The standard deviationof such a distribution is defined as the standard error of the mean σ X . It is known that for Gaussian distributions,the relation σ X /σ = 1 / √ n holds, where σ is the standard deviation of the entire set, so as n grows, the estimate ofthe mean improves.To associate the value of n to the time signals obtained in our simulations, we divided the entire run time in n segments of length given by the decorrelation time T d , so n = T r /T d , then, the relation between the standard error ofthe mean and the standard deviation is now σ X σ = r T d T r , (30)where σ is now the standard deviation of the signal in a simulation with T r → ∞ . As a rule of thumb, the standarderror of the mean should be no greater than 10% of the standard deviation, so a value of T r ≈ T d will ensure thiscondition. That is the reason why in the case of measuring Lyapunov exponents, for which T d ≈ . ≈ T E /
10, arun time of T r = 50 ≈ T E is enough to obtain proper statistics, whereas for the energy, that same run time isnot enough. Note that the above criteria to determine the run time is uniquely dependent on the decorrelation time,since it is not necessary to measure σ and σ X once we observe that the distribution is approximately Gaussian.On the other hand, the proper run time needed to obtain proper statistics for the energy is approximately 100 T E ,which is ten times greater than that needed in the case of Lyapunov exponents. The measurement of Lyapunovexponents are the most robust in these simulations. This suggests that the Lyapunov exponents are not only stableto variations in the steptime and lattice size, but also the run time needed to compute them is much faster thanfor other flow quantities. Thus, using relations such as Ruelle’s, that connect the Lyapunov exponent to Reynoldsnumber, is useful to attempt to characterize the flow in a quick and robust way. V. MHD
Twelve simulations were run using the incompressible MHD equations described in Equations (6 - 7). Thesesimulations keep magnetic Prandtl number, Pm = ν/η , set to 1 for all simulations. The forcing is set so that the sumof magnetic and kinetic dissipation is approximately 0.1 and magnetic helicity is kept low. The FTLE method wasused to obtain the Lyapunov exponents. The values of Re ranged between 40 and 1300. Figure 18 shows histograms forthe distribution of λ and Re obtained from one of the simulations. Table II shows parameter values for all simulations.Similar to the hydrodynamic case, whilst the λ are well approximated by a Gaussian distribution, the value for Re isnot. As is known, analysis of the chaotic properties of MHD simulations presents noisier results than the NSE case[17]. As well, we find that characteristic quantities and their fluctuations are less correlated. MHD involves more8 N ν ε Re λ σ λ T τ k max η ν is the kinematic viscosity, ǫ is the kinematicdissipation rate, Re is the Reynolds number, λ is the maximal Lyapunov exponent and σ λ its standard deviation, T is the largeeddy turnover time, τ is the Kolmogorov microscale time, k max ≈ N/ − η = ( ν /ε ) / is the Kolmogorov microscale length and k max η is the resolution. .
15 0 .
20 0 .
25 0 .
30 0 .
35 0 . ˜ λ P r o b a b ili t y d e n s i t y (a)
75 80 85 90 95 Re . . . . . . . . p r o b a b ili t y d e n s i t y (b) FIG. 18. Distributions of λ and Re for the MHD case for one simulation. time and length scales than magnetically neutral flows, given that for each velocity field length and time scale, thereis a corresponding time or length scale for the magnetic field. As well, there are new relations between these scalesthat are possible given the presence of magnetic diffusion. Hence, it is to be expected that some dependencies thatare obtained from dimensional considerations in hydrodynamics are lost in MHD.The first relationship to test for MHD simulations is the Ruelle relation λ ∝ /τ . Here, τ = ( ν/ǫ ) / is theKolmogorov time, and the dissipation is the kinetic dissipation. Figure 19 shows λ against 1 /τ . A linear fit λ = A + Bτ − is done. The measured values are A = 0 . ± .
02 and B = 0 . ± . τ − , showing that for the values of τ − we explore, the linearapproximation is not accurate whereas for those values that are not explored here the linear approximation becomesreasonable. Exploring a wider range of values to evaluate this relation is not the main goal of this project and thatwould require computational cost which is not pursued for this work.Figure 20 shows the relationship between λ and Re. This plot allows us to test the scaling relation in Equation (2).From the Figure, we measure a value of α = 0 . ± .
02, which is quite far from the value of α ≈ . σ λ and Re is very weak, contrary to the hydrodynamic case, where there was a strongcorrelation as is shown in Figure 9.9 τ − . . . . . . λ FIG. 19. λ vs τ − (the inverse of Kolmogorov time τ = ( ν/ǫ ) / ) and linear fit represented by the solid red line. Re × − λ T FIG. 20. λT vs. Re using the FTLE method in a MHD fluid We now compare the dependence of σ λ on λ and σ Re on Re, as was done for hydrodynamics previously. Figure 21shows these relations for our MHD simulations. There are similarities and differences with the hydrodynamic results.Unlike for NSE, in MHD the fluctuations in the Lyapunov exponents are not strongly related to λ . However, similarto the NSE case, in MHD the fluctuations of Re are strongly related to Re itself. From the Figure, a linear fit for theMHD simulations of σ Re = c Re results in a measured value for the slope of c = 0 . ± . c = 0 . ± . λ and 1 /τ is strong, the strong relationship between τ and Re only holds in hydrodynamics.Since these relations are not uniquely determined in MHD, the behavior of λ is not as strongly related to Re. Thelack of a uniquely determined relation may be the reason why the linear relation between σ λ and λ is not present inMHD.Similar to hydrodynamics, we test the sensitivity of the MHD FTLE results to changes in steptime. Figure 22shows the dependence of λ and σ λ on the steptime ∆ t of the FTLE method for one MHD simulation. As in thehydrodynamic case, both λ and σ λ seem to be stable even though a very wide range of ∆ t is explored. VI. DISCUSSION AND CONCLUSIONS
We studied Eulerian chaos in homogeneous and isotropic flows evolved using the incompressible Navier-Stokesequations. We observed the finite time Lyapunov exponents and the dependence of fluctuations on other parameters0 Re σ R e FIG. 21. Distributions of λ and Re for the MHD case. No correlation was shown between λ and σ λ . The relation between Reand its fluctuations presented a linear relation with a fit given by the solid line. . . . . . . . . λ . . . . . . ∆ t . . . σ λ FIG. 22. Dependence of λ and σ λ on the steptime ∆ t of the FTLE method for MHD flows. such as Reynolds number, lattice size and steptime in the FTLE procedure. There are four main results which wemay take from the analysis done in this paper.The first is that the FTLEs are approximately Gaussian distributed. We do not expect that they are fully Gaussianhowever, because the statistics on which it depends, and turbulent statistics in general, are known to be non-Gaussian.We may try to find out the true distribution of the FTLEs by relating it to previous results found by analysingHIT with dynamical systems theory. Previously, it has been found that when the NSE are modelled with a negativedamping force, at low Re a relaminarization process occurs akin to that found in parallel wall-bounded shear flows[61]. The lifetime of this process was found to depend super-exponentially on the Reynolds number. We can associatethe parts of the distribution of the FTLE where λ < E/ε rather than
L/U . The study of the chaotic behaviour in decayingturbulence shows that the evolution of the divergence is quite sensitive to the different functional forms proposed.Numerical data supports the idea that
E/ε is the one that best fit the data. We also see that the lattice size is notreally a factor, and under-resolved simulations also capture the chaotic behaviour. As such, the smallest length scalesmay not be as relevant to the chaos and as shown in previous works [14, 15, 19], these are the scales in which the twofield initially close become decorrelated first. Conversely, large scale properties take a longer time to decorrelate andthus are more predictable.Third, from the analysis of Figure 12, we see that both the Lyapunov exponent and its fluctuations are very weaklydependent on the steptime used in the FTLE method. Practically, this result means that the time needed to get ameasure of this chaotic property can be achieved at low computational cost.Fourth, motivated by understanding the chaotic properties of the system, we analyzed the timescale of the signalsof the energy, dissipation rate and Lyapunov exponent in the steady state. We look at the time autocorrelation andwe see that the energy only becomes decorrelated after 2 T E . However, the small scale properties become decorrelatedfaster than T E and the Eulerian maximal Lyapunov exponent becomes decorrelated faster than any of these otherstatistics, suggesting that it is the most robust to measurement in simulations, reaching a faster stability. From thisanalysis we obtain a useful rule that determines the run time T r needed to obtain proper statistics of a quantity byconsidering the decorrelation time in its signal T d , being T r ≈ T d . For instance, we see that for the total energy,averages need to be taken for run times T r ≈ T E , whereas for other properties such as Lyapunov exponents ordissipation rate, shorter run times can be used. This run time rule we give in relation to the self correlation timeis useful not only when looking at chaotic properties of the flow but also for any other study that needs to measurequantities in the steady state.The robustness of the results when measuring the maximal Lyapunov exponent suggest that it is a stable statisticalproperty of the system. We also see that it reaches stability faster than other statistical properties of the system andour analysis cements it as an extra parameter of HIT which is important in understanding the interesting inter-scaledynamics inherent in turbulent flows. This could have further applications in studies that use the Lyapunov exponentto characterize the turbulent flow as in [18]. Furthermore, if we learn the relation between the Lyapunov exponentand other flow quantities, we could use such an exponent to determine these other quantities in a faster and morestable way. For instance, Ruelle’s relation relates λ and Re and this can be used to estimate Re in a shorter run time.Finally, in this paper some of the FTLEs’ properties measured for hydrodynamic flows were measured for the case ofisotropic and incompressible MHD flows. The distribution of FTLEs is approximately Gaussian as in the case of NSE.The remarkable stability of FTLEs against variations in the steptime is also present in MHD flows. Nevertheless, thedependence on Reynolds number is not the same as in the case of NSE. Ruelle’s relation is derived from assumingthat the Kolmogorov time is the smallest timescale in the system. Whereas this is true for incompressible NSE, othercharacteristic timescales are present in the smallest scales in MHD that may break this relation. ACKNOWLEDGEMENTS
This work used the Cirrus UK National Tier-2 HPC Service at EPCC funded by the University of Edinburghand EPSRC (EP/P020267/1). R.D.J.G.H. was supported by the U.K. Engineering and Physical Sciences ResearchCouncil (EP/M506515/1) and A.A. was supported by the University of Edinburgh. A.B. acknowledges funding fromthe U.K. Science and Technology Facilities Council. [1] D. Ruelle and F. Takens, Comm. Math. Phys. , 167 (1971).[2] R. Deissler, Phys. Fluids , 1453 (1986).[3] E. N. Lorenz, J. Atmospheric Sci. , 130 (1963).[4] S. Yoden and M. Nomura, J. Atmospheric Sci. , 1531 (1993).[5] T. Masuoka, Y. Takatsu, and T. Inoue, Microscale Thermophys. Eng. , 347 (2003).[6] A. Musker, AIAA J. , 655 (1979).[7] D. C. Poirel and S. J. Price, AIAA J. , 1960 (2001).[8] C. Leith, J. Atmospheric Sci. , 145 (1971).[9] C. Leith and R. Kraichnan, J. Atmospheric Sci. , 1041 (1972).[10] E. Aurell, G. Boffetta, A. Crisanti, G. Paladin, and A. Vulpiani, J. Phys. A: Math. Gen. , 1 (1997).[11] J. Kurths and H. Herzel, Sol. Phys. , 39 (1986).[12] R. Grappin, J. Leorat, and A. Pouquet, J. de Physique , 1127 (1986).[13] W. Huang, W. Ding, D. Feng, and C. Yu, Phys. Rev. E , 1062 (1994).[14] A. Berera and R. D. J. G. Ho, Phys. Rev. Lett. , 024101 (2018).[15] G. Boffetta and S. Musacchio, Phys. Rev. Lett. , 054102 (2017). [16] P. Mohan, N. Fitzsimmons, and R. D. Moser, Phys. Rev. Fluids , 114606 (2017).[17] R. D. Ho, A. Berera, and D. Clark, Phys. Plasmas , 042303 (2019).[18] G. Nastac, J. W. Labahn, L. Magri, and M. Ihme, Phys. Rev. Fluids , 094606 (2017).[19] S. Mukherjee, J. Schalkwijk, and H. J. Jonker, J. Atmospheric Sci. , 2715 (2016).[20] Y. C. Li, R. D. Ho, A. Berera, and Z. Feng, arXiv preprint arXiv:1908.04838 (2019).[21] G. Boffetta, A. Celani, A. Crisanti, and A. Vulpiani, Phys. Fluids , 724 (1997).[22] G. Boffetta and S. Musacchio, Phys. Fluids , 1060 (2001).[23] O. M´etais and M. Lesieur, J. Atmospheric Sci. , 857 (1986).[24] A. Crisanti, M. Jensen, A. Vulpiani, and G. Paladin, Phys. Rev. Lett. , 166 (1993).[25] M. Yamada and Y. Saiki, Nonlinear Process. Geophys. , 631 (2007).[26] N. de Divitiis, Adv.Math. Phys. (2018).[27] G. Lapeyre, Chaos: An Interdiscip. J. NonlinearSci. , 688 (2002).[28] L. Biferale, G. Boffetta, A. Celani, B. Devenish, A. Lanotte, and F. Toschi, Phys. Fluids , 115101 (2005).[29] L. Biferale, G. Boffetta, A. Celani, B. Devenish, A. Lanotte, and F. Toschi, Phys. Fluids , 111701 (2005).[30] S. Vannitsem, Chaos: An Interdiscip. J. NonlinearSci. , 032101 (2017).[31] H. Guo, W. He, T. Peterka, H.-W. Shen, S. M. Collis, and J. J. Helmus, IEEE Trans. Vis. Comput. Graph. , 1672(2016).[32] D. Garaboa-Paz, J. Eiras-Barca, and V. P´erez-Mu˜nuzuri, Earth Syst. Dyn. , 865 (2017).[33] F. d’Ovidio, V. Fern´andez, E. Hern´andez-Garc´ıa, and C. L´opez, Geophys. Res. Lett. (2004).[34] G. Haller and T. Sapsis, Chaos: An Interdiscip. J. NonlinearSci. , 023115 (2011).[35] G. Haller, Annu. Rev. Fluid Mech. , 137 (2015).[36] E. R. Abraham and M. M. Bowen, Chaos: An Interdiscip. J. NonlinearSci. , 373 (2002).[37] A. C.-L. Chian, S. S. Silva, E. L. Rempel, M. Goˇsi´c, L. R. B. Rubio, R. A. Miranda, and I. S. Requerey, arXiv:1904.08472(2019).[38] E. Rempel, A.-L. Chian, F. J. Beron-Vera, S. Szanyi, and G. Haller, Mon. Notices Royal Astron. Soc. , L108 (2016).[39] A. C.-L. Chian, E. L. Rempel, G. Aulanier, B. Schmieder, S. C. Shadden, B. T. Welsch, and A. R. Yeates, The Astrophys.J , 51 (2014).[40] J.-N. Tang, C.-C. Tseng, and N.-F. Wang, Acta Mech. Sin. , 612 (2012).[41] T. S. Pedersen, P. K. Michelsen, and J. J. Rasmussen, Phys. Plasmas , 2939 (1996).[42] K. Padberg, T. Hauff, F. Jenko, and O. Junge, New J. Phys. , 400 (2007).[43] M. Falessi, F. Pegoraro, and T. Schep, J. Plasma Phys. (2015).[44] J. Misguich, R. Balescu, H. Pecsell, T. Mikkelsen, S. Larsen, and Q. Xiaoming, Plasma Phys. Control. Fusion , 825(1987).[45] S. C. Shadden, F. Lekien, and J. E. Marsden, Phys. D: Nonlinear Phenom. , 271 (2005).[46] P. L. Johnson and C. Meneveau, Phys. Fluids , 085110 (2015).[47] I. Shimada and T. Nagashima, Prog. Theor. Exp. Phys. , 1605 (1979).[48] D. Ruelle, Phys. Lett. A , 81 (1979).[49] A. N. Kolmogorov, Proc. Royal Soc. Lond. A , 15 (1991).[50] W. D. Mccomb, A. Berera, S. R. Yoffe, and M. F. Linkmann, Phys. Rev. E (2015).[51] L. Machiels, Phys. Rev. Lett. , 3411 (1997).[52] T. Ishihara and Y. Kaneda, in Statistical Theories and Computational Approaches to Turbulence (Springer, 2003) pp.177–188.[53] S. R. Yoffe, arXiv:1306.3408 (2013).[54] S. L. Brunton and C. W. Rowley, Chaos: An Interdiscip. J. NonlinearSci. , 017503 (2010).[55] V. I. Oseledets, Trans. Moscow Mater. Soc. , 179 (1968).[56] I. Goldhirsch, P.-L. Sulem, and S. A. Orszag, Phys. D: Nonlinear Phenom. , 311 (1987).[57] G. I. Taylor, Proc. Royal Soc. Lond. A , 421 (1935).[58] C. R. Doering and P. Constantin, Phys. Rev. E , 4087 (1994).[59] K. R. Sreenivasan, Phys. Fluids , 528 (1998).[60] T. Gotoh, D. Fukayama, and T. Nakano, Phys. Fluids , 1065 (2002).[61] M. F. Linkmann and A. Morozov, Phys. Rev. Lett.115