Scaling of acceleration statistics in high Reynolds number turbulence

The scaling of acceleration statistics in turbulence is examined by combining data from the literature with new data from well-resolved direct numerical simulations of isotropic turbulence, significantly extending the Reynolds number range. The acceleration variance at higher Reynolds numbers departs from previous predictions based on multifractal models, which characterize Lagrangian intermittency as a naive extension of Eulerian intermittency. The disagreement is even more prominent for higher-order moments of the acceleration. Instead, starting from a known exact relation, we relate the scaling of acceleration variance to that of Eulerian fourth-order velocity gradient and velocity increment statistics. This prediction is in excellent agreement with the variance data. Our work highlights the need for models that consider Lagrangian intermittency independent of the Eulerian counterpart.


Introduction:
The acceleration of a fluid element in a turbulent flow, given by the Lagrangian derivative of the velocity, resulting from the balance of forces acting on it, is arguably the simplest descriptor of its motion. This is directly reflected in the Navier-Stokes equations: where, u is the divergence-free velocity (∇ · u = 0), p the kinematic pressure, ν is the kinematic viscosity and f is a forcing-term. Besides its fundamental role in the study of turbulence [1][2][3][4], understanding the statistics of acceleration is of paramount importance for diverse range of applications constructed around stochastic modeling of transport phenomena in turbulence [5][6][7][8]. The application of Kolmogorov's 1941 phenomenology implies that the variance (and higher-order moments) of any acceleration component a can be solely described by the mean-dissipation rate ǫ and ν [9][10][11]: where a 0 is thought to be a universal constant. However, extensive numerical and experimental work has shown that a 0 increases with Reynolds number [12][13][14][15][16][17][18][19][20]. Thus, obtaining data on a 0 and modeling its R λvariation has been a topic of considerable interest. While several theoretical works have focused on acceleration statistics [21][22][23][24], the most notable procedure -but whose validity should not be taken for granted -stems from the multifractal model [25][26][27], which quantifies acceleration intermittency (and, in general, the intermittency of other Lagrangian quantities) by adapting to the Lagrangian viewpoint the well-known Eulerian framework, based either on the energy dissipation rate [28] or velocity increments [29]. A key result from this consideration is that a 0 ∼ R χ λ , χ ≈ 0.135, where R λ is the Taylorscale Reynolds number. While data from direct numerical simulations (DNS) and experiments do not directly Normalized acceleration variance a0 = a 2 /( ǫ 3/2 ν −1/2 ) as a function of R λ . The data can be prescribed by a simple R 0.25 λ power-law at high R λ , in contrast to the previously proposed empirical fit given by Eq. (3). display this power-law, it was nevertheless presumed to to be asymptotically correct at very large R λ , and an empirical interpolation formula [16,19], with χ = 0.135, c 1 = 1.9, c 2 = 85, was suggested to fit the data, showing reasonable success [16]. An alternative scaling: a 0 ∼ R 0.25 λ was proposed by Hill [21], which was indistinguishable from Eq.(3) at low-R λ [16,20]; we discuss the veracity of this proposal later.
In this Letter, we revisit the scaling of acceleration variance (and higher-order moments) by presenting new DNS data at higher R λ . The new variance data agrees with previous lower R λ data, where the R λ range over-laps, but increasing deviations from Eq. (3) occur at higher R λ . Results for high-order moments show even stronger deviations from previous predictions. Further analysis shows that the extension of Eulerian multifractal models to the Lagrangian viewpoint is the source of this discrepancy. We develop a statistical model which shows excellent agreement with variance data at high R λ , and also provide an updated interpolation fit to include low R λ data.
Direct Numerical Simulations: The DNS data utilized here correspond to the canonical setup of forced stationary isotropic turbulence in a periodic domain [30], allowing the use of highly accurate Fourier pseudo-spectral methods [31]. The novelty is that we have simultaneously achieved very high Reynolds number and the necessary grid resolution to accurately resolve the small-scales [32,33]. The data correspond to the same Taylor-scale Reynolds number R λ range of 140 − 1300 attained in recent studies [34][35][36][37], which have adequately established convergence with respect to resolution and statistical sampling. The grid resolution is as high as k max η K ≈ 6 which is substantially higher than k max η K ≈ 1 − 2, utilized in previous acceleration studies [16,19,20]; k max = √ 2N/3 is the maximum resolved wavenumber on a N 3 grid and η K = ν 3/4 ǫ −1/4 is the Kolmogorov length scale. This improved small-scale resolution is especially necessary for capturing higher-order statistics of acceleration, since acceleration is even more intermittent than spatial velocity gradients [12,19].
Acceleration variance: Figure 1 shows the compilation of data from various sources including data from both DNS [19,20] [38] and bias-corrected experiments [39]. We have also included DNS data obtained directly from Lagrangian trajectories of fluid particles [40][41][42], which give identical results for acceleration variance [43]. As evident, while Eq. (3) works for the previous range of R λ , it does not fit the new data. In fact, a R 0.25 λ scaling is more appropriate at higher R λ , and as discussed later, the failure of multifractal models in fitting higher-order moments is even more conspicuous. To gain clarity on this point, it is useful to discuss the multifractal models first.
Acceleration scaling from multifractals: The key idea in multifractal approaches is to quantify the intermittency of acceleration in terms of the intermittency of Eulerian velocity gradients or dissipation rate. Assuming a simple phenomenological equivalence between temporal and spatial derivatives, acceleration can be written in terms of dissipation rate and viscosity as a ∼ ǫ 3/4 ν −1/4 . Thus, the moments of acceleration are obtained as: Alternatively, where a K = ǫ 3/4 ν −1/4 , i.e., acceleration based on Kolmogorov variables. Since Eulerian intermittency dictates that ǫ q = ǫ q for any q = 1 [28,29], the key assumption in its extension to Lagrangian intermittency is that the p-th moment of acceleration scales as the (3p/4)-th moment of ǫ [25,27]. The scaling of ǫ q / ǫ q can be obtained by several approaches, all leading to similar results. We briefly summarize a few approaches below, with additional details in the Supplementary Material [44].
The above framework leads to the result: An approximation for F (α), such as the p-model [45,47], can be used to obtain τ q . The p-th moment of acceleration can then be simply obtained as [48] Instead of dissipation, one can also start by taking the velocity increment δu r over scale r to be Hölder continuous: δu r /u ′ ∼ (r/L) h , where h is the local Hölder exponent and D(h) is the corresponding multifractal spectrum. A scale-dependent dissipation rate ǫ r can then be defined as ǫ r ∼ (δu r ) 3 /r, which reduces to the true dissipation for the viscous-cutoff defined by the condition δu r r/ν ≃ 1. This framework leads to the same result as in Eq. (6), corresponding to α = 3h and F (α) = D(h). A well-known approximation for D(h) is given by the She-Leveque model [49]. Finally, we can also use the  log-normal model [50], which gives τ q = 3µq(q − 1)/4, even though it is untenable for very large q [29]. Here, µ is the intermittency exponent, with experiments and DNS suggesting µ ≈ 0.25 [51,52].
The scaling of acceleration moments obtained from these three approaches and also from DNS data are listed in Table I, up to sixth-order. All approaches give essentially the same result for the acceleration variance, with the exponent of about 0.135 used in Eq. (3). However, the high-R λ DNS data clearly do not conform to any of the power laws shown in Table I. The results for normalized fourth and sixth order moments, also plotted in Fig. 2, clearly show that the power-laws increasingly differ from multifractal predictions.
As noted earlier, the use of multifractals is primarily motivated by Eq. (4). To get better insight, in Fig.3a, we plot a 0 and ǫ 3/2 / ǫ 3/2 versus R λ . While the latter shows a clear R 0.14 λ scaling as anticipated from multifractals (and also the log-normal model), the former shows a steeper scaling of R 0. 25 λ . An even more general and direct test is presented in Fig.3b, by checking the validity of Eq. (4) for different p values. The data clearly suggest that the acceleration intermittency, being stronger, cannot be described by extending the Eulerian intermittency of the dissipation-rate. In fact, a similar observation has been made for Lagrangian velocity structure functions, where extensions of the p-model and the She-Leveque model severely underpredicts their intermittency (i.e., overpredicts the inertial-range exponents) [53].
It is worth considering if one might describe the scaling of acceleration moments in terms of enstrophy Ω = |ω| 2 (ω being the vorticity), instead of dissipation. This change addresses the likelihood that acceleration is influenced more by transverse velocity gradients than by longitudinal ones [54,55]. In isotropic turbulence Ω = ǫ /ν, but the higher moments differ; enstrophy being more intermittent [32,56]. The resulting modification to Eq. (4) is: a p ∼ Ω 3p/4 ν 2p/4 . However, as tested in Fig. 3a,b, the differences arise only for large p; even then, it is not sufficient to explain the stronger intermittency of acceleration (also see Supplementary [44]).
Acceleration variance from fourth-order structure function: A statistical model for acceleration variance is now obtained using a methodology similar to that proposed by [21], but differing in some crucial aspects. From Eq. (1), acceleration variance can be obtained directly as [57] The viscous contribution is known to be small and can be ignored [13]. An exact relation for variance of pressure-gradient is also known [58,59]: where the Ds are the fourth order longitudinal, transverse and mixed structure functions, in order. The above results can be rewritten as [21]: where H χ is defined by Eqs. (8)- (9). At sufficiently high R λ ( 200), DNS data [13,20] confirm that H χ ≈ 0.65 (also see Supplementary [44]). We can normalize both sides by Kolmogorov-scales to write Assuming standard scaling regimes [29], we can write where F is the flatness of ∂u/∂x, ξ 4 is the inertial-range exponent, and C 4 , C are constants which depend on R λ ; ℓ is a crossover scale between the viscous and inertial range and is determined by matching the two regimes as Now, taking we have Finally, from piecewise integration of Eq. (11), it can be shown that (see Supplementary Material [44] for intermediate steps): Substituting the R λ -dependencies, we get The values of α, β and ξ 4 are in principle obtainable from Eulerian intermittency models. The exponent α simply corresponds to τ 2 in Eq. (6), since F ∼ ǫ 2 / ǫ 2 . Multifractal and log-normal models predict α = τ 2 ≈ 0.38. The DNS data for F are shown in Fig. 4a, giving α ≈ 0.387, in excellent agreement with the prediction, and also with previous experimental and DNS results in literature [18,20].
On the other hand, intermittency models predict ξ 4 ≈ 1.28 [49]. Our DNS shows ξ 4 ≈ 1.3, which is well within statistical error bounds. Finally, the prediction for β from multifractal model is β ≈ (4−3ξ 4 )/2, which reduces to β ≈ µ/3 for log-normal model; both predictions give β ≈ 0.08 (also see Supplementary [44]). Figure 4b shows the normalized fourth-order structure function from our DNS data, using ξ 4 ≈ 1.3. Note, as expected, the inertialrange increases with R λ . The inset of the bottom panel shows C 4 , giving β ≈ 0.2. This observed β is substantially larger than 0.08 anticipated from multifractal and log-normal models.
The use of α = 0.387, β = 0.2 and ξ 4 = 1.3 in Eq. (17) leads to which is in excellent agreement with the high-R λ data shown in Figs. 1 and 3a. The exponent 0.25 is virtu-ally insensitive to a small variation in ξ 4 , but is significantly impacted by the choice of β = 0.2 (instead of 0.08). Moreover, the use of β ≈ 0.08 in Eq. (17) gives a 0 ∼ R 0.15 λ , which is essentially the same as the exponent 0.14 obtained earlier in Table I. This shows the robustness of piecewise integration leading to the result in Eq. (17) and also suggests that the discrepancy from multifractal prediction is due to the exponent β (and hence the proportionality constant C 4 ). In this regard, the role of β needs to be further explored, especially in relation to the inadequacy of Eq. (4).
We note that the exponent 0.25 was also suggested by Hill [21]. However, Hill arrived at this result by deriving that a 0 ∼ F 0.79 and F ∼ R 0.31 λ based on [60]; evidently, the current data do not agree with both of these results. It appears that the two errors fortuitously cancelled out each other to give the 0.25 exponent. Finally, we point out that the exponent 0.25 describes the data for R λ 200. To describe the data at lower R λ , an empirical interpolation formula in the spirit of Eq. (3) can be devised with χ = 0.25. Least-square fit gives c 1 ≈ 0.89, c 2 ≈ 40 (also see Supplementary [44]).
Conclusions: The moments of Lagrangian acceleration are known to deviate from classical K41 phenomenology due to intermittency. Attempts were made to quantify these deviations by extending the Eulerian multifractal models to the Lagrangian viewpoint and devising an ad-hoc interpolation formula to agree with available data from DNS and experiments. The first contribution of this article is to present new very well resolved DNS data on Lagrangian acceleration at higher R λ , and show that they disagree with the results from multifractal models, and the interpolation formula. The disagreement gets increasingly stronger with the moment order. As a second contribution, the article devises a statistical model that is able to correctly capture the scaling of acceleration variance. While this framework does not seem amenable for generalization to higher-order moments, our results show that the intermittency of Lagrangian quantities remains an open problem, even more compellingly than before.
Acknowledgments: We thank P.K. Yeung and Luca Biferale for providing helpful comments on an earlier draft of the manuscript. We gratefully acknowledge the Gauss Centre for Supercomputing e.V. for providing computing time on the supercomputers JUQUEEN and JUWELS at Jülich Supercomputing Centre (JSC), where the simulations reported in this paper were performed. * dhawal.buaria@nyu.edu [1]  The underlying idea behind all approaches is that acceleration a can be dimensionally written in terms of dissipation rate ǫ and viscosity ν as Thus, the moments of acceleration are obtained as Normalizing both side by Kolmogorov scales, we have where a K = ǫ 3/4 ν −1/4 . Thus, our goal is to obtain the normalized moments of dissipation. We present here the results from various intermittency models, providing extra details compared to the main text.

A. Multifractal approach based on velocity increments
Within the multifractal framework, the velocity increment δu r over scale r is assumed to be Hölder continuous: where h is the local Hölder exponent and D(h) is the corresponding multifractal spectrum. Within this paradigm, a scale-dependent dissipation ǫ r can be defined as which can be rewritten as where we have used ǫ ∼ u ′3 /L based on dissipation anomaly. We can obtain the moments of ǫ r [1] as: Now ǫ r reduces to true dissipation when r corresponds to the viscous cutoff, often defined by the condition which can be rewritten as where Re = u ′ L/ν. Using the steepest gradient argument, and Re ∼ R 2 λ , the scaling of the dissipation moments at large Re can be obtained as Thus, it readily follows that The exponents can be obtained for any standard D(h). A well known approximation is given by the She-Leveque model [2]: where d * = (1 − 3h * )/(1 − γ), h * = 1/9 and γ = 2/3.
Longitudinal versus transverse directions : In principle, the multifractal model does not differentiate between longitudinal and tranverse directions, i.e., it predicts that velocity increments and gradients in both longitudinal and transverse directions scale similarly. However, it is now well established that this is not the case,with tranverse increments/gradients being somewhat more intermittent [3][4][5][6]. Based on this, an ad-hoc modification to D(h) in Eq. (12) was proposed in ref. [4], with h * = 1/9 but γ = 1/2, to empirically fit the scaling exponents of transverse structure functions. Thereafter, this result was later utilized to predict the scaling of acceleration variance in [7]. However, we note that this extension is unjustified, since this modified D(h) (originally intended to describe transverse structure functions) already fails to describe the scaling of transverse velocity gradients. For instance, this modified D(h) predicts the following for enstrophy: Ω 3/2 / Ω 3/2 ∼ R 0.23 λ (as opposed to the 0.14 predicted for scaling dissipation). However, this prediction is at clear odds with DNS data in Fig. 3a, which shows that both dissipation and enstrophy moments are close to the 0.14 prediction. Essentially, this ad-hoc modification to D(h) grossly overpredicts the scaling of transverse velocity gradients and hence cannot be utilized to describe the scaling of acceleration also.

B. Multifractal approach based on dissipation
The multifractal approach based on dissipation directly deals with the dissipation field, with the starting point being [8] where α is the local Hölder exponent, with the multifractal spectrum given by some F (α). In the literature [8], the 1D spectrum f (α) is more commonly used, which is simply given by f (α) = F (α) − 2. As is evident, this expression is essentially equivalent to that in Eq. (6), with α = 3h. It also follows that the moments of ǫ r are given as Now, the viscous cutoff at which ǫ r reduces to true dissipation is defined by the scale which, by using Eq. (13), leads to Here, we have also used ǫ ∼ u ′3 /L and Re = u ′ L/ν. Using the steepest gradient argument, and Re ∼ R 2 λ , the scaling of the dissipation moments at large Re can be obtained as In the literature, the result is often stated in terms of the Renyi dimension D q , which is the Legendre transform of F (α). Stated in terms of D q the previous result can be rewritten as [9] τ q = 6(q − q) , where the unique functionq(q) is determined from An approximation for D q is provided by the p-model [10]: for p 1 = 0.7 and p 2 = 1 − p 1 .

C. Kolmogorov's lognormal model
From the  lognormal model [11], it is well-known that where µ is the intermittency exponent, with standard estimates suggesting µ = 0.25 [12]. This result also corresponds to D q = 1 − µq/2.

II. INTEGRATION OF FOURTH ORDER STRUCTURE FUNCTION
As noted in the main text, acceleration variance can be obtained by piecewise integration of the fourth-order structure function D 1111 (r): where H χ is defined by the above relation and is assumed to be constant at sufficiently high R λ . This is confirmed by previous [3] and current DNS data, both plotted in Fig. 1, showing H χ ≈ 0.65 for R λ 200. The right hand side in Eq. (22) can be obtained via piecewise integration, utilizing: where F is the flatness of ∂u/∂x, ξ 4 is the corresponding inertial-range exponent. Note that both C 4 and C are functions of R λ . The scale ℓ corresponds to the crossover between dissipation and inertial range, and is obtained by invoking continuity of D 1111 (r) at r = ℓ, i.e., (F/225)(ℓ/η K ) 4 = C 4 (ℓ/η K ) ξ4 , which leads to To simplify notation, we momentarily take r/η K → r, ℓ/η K → ℓ, L/η K → L and evaluate the integral in Eq. (22) as As ξ 4 < 2, we can ignore the last two terms with L, since they will both decrease with increasing R λ (note C ∼ R λ 2 and L ∼ R λ 3/2 [1]). Thus, we get Note that we have F ℓ 4 ∼ C 4 ℓ ξ4 from Eq. (24), which implies F ℓ 2 ∼ C 4 ℓ ξ4−2 . Thus, both terms in Eq. (28) have the same scaling, and we can simply write where we have replaced ℓ → ℓ/η K . The R λ scaling of F and C 4 are now taken to be which, upon substitution in Eq. (24), also gives Finally, substituting these relations in Eq. (29), we get a 0 ∼ R λ (2α−αq+2β)/(4−ξ4) , which is reported in the main text. Using α ≈ 0.387, β ≈ 0.2 and ξ 4 ≈ 1.30, we get a 0 ∼ R λ 0.25 .
III. UPDATED INTERPOLATION FORMULA FOR a0 VS. R λ As evident from Fig. 1 (and Fig. 3a) of the main text, the result a 0 ∼ R λ 0.25 describes the data excellently for R λ 200. To approximately model the data for lower R λ , it might be worthwhile to consider an interpolation formula similar to Eq. (3) of the main text, with χ = 0.25. A least-square fit gives c 1 ≈ 0.89 and c 2 ≈ 40. The fit is shown in the figure below. We note that this result is only meant to be an empirical approximation. In principle, other interpolation formulae can also be devised to obtain a more robust fit at low R λ .

IV. SCALING OF STRUCTURE FUNCTION PROPORTIONALITY CONSTANTS
The inertial-range scaling of p-th order structure functions can be written as: where C p are the proportionality constants and ζ p is the inertial range. As per K41, ζ p = p/3 and C p are universal constants independent of R λ . However, due to intermittency ζ p is a non-trivial function of p and C p depends on R λ (except for p = 3). To derive the R λ -scaling of C p , we start again from Eq. (4). Standard multifractal arguments [1] give: Normalizing the above equation by Kolmogorov variables gives: Thereafter, using u ′ /u K ∼ R λ 1/2 and L/η K ∼ R λ 3/2 [1] leads to: Comparing with Eq. (34) gives: For p = 4, we have C 4 ∼ R λ β . Thus, the multifractal prediction gives β = (4 − 3ζ 4 )/2. For K62 log-normal, we have ζ p = 4/3 − 2µ/9, giving β = µ/3. Both multifractals and K62 log-normal predict ζ 4 ≈ 1.28 [1], leading to β ≈ 0.08, which is substantially lower than the β ≈ 0.2 suggested by DNS results, as reported in the main text.