Extreme velocity gradients in turbulent flows

Fully turbulent flows are characterized by intermittent formation of very localized and intense velocity gradients. These gradients can be orders of magnitude larger than their typical value and lead to many unique properties of turbulence. Using direct numerical simulations of the Navier-Stokes equations with unprecedented small-scale resolution, we characterize such extreme events over a significant range of turbulence intensities, parameterized by the Taylor-scale Reynolds number ($R_\lambda$). Remarkably, we find the strongest velocity gradients to empirically scale as $\tau_K^{-1} R_\lambda^{\beta}$, with $\beta \approx 0.775 \pm 0.025$, where $\tau_K$ is the Kolmogorov time scale (with its inverse, $\tau_K^{-1}$, being the {r.m.s.} of velocity gradient fluctuations). Additionally, we observe velocity increments across very small distances $r \le \eta$, where $\eta$ is the Kolmogorov length scale, to be as large as the {r.m.s.} of the velocity fluctuations. Both observations suggest that the smallest length scale in the flow behaves as $\eta R_\lambda^{-\alpha}$, with $\alpha = \beta - \frac{1}{2}$, which is at odds with predictions from existing phenomenological theories. We find that extreme gradients are arranged in vortex tubes, such that strain conditioned on vorticity grows on average slower than vorticity, approximately as a power law with an exponent $\gamma<1$, which weakly increases with $R_\lambda$. Using scaling arguments, we get $\beta=(2-\gamma)^{-1}$, which suggests that $\beta$ would also slowly increase with $R_\lambda$. We conjecture that approaching the limit of infinite $R_\lambda$, the flow is overall smooth, with intense velocity gradients over scale $ \eta R_\lambda^{-1/2}$, corresponding to $\beta = 1$.


I. INTRODUCTION
Quantitative studies of turbulence in incompressible flows reveal that the averaged dissipation rate of turbulent kinetic energy, ǫ , is independent of kinematic viscosity, ν, when ν → 0 or equivalently when the turbulence intensity, i.e., the Reynolds number, is very high [1,2]. This empirical result, also known as the zeroth law of turbulence, implies that the amplitude of velocity gradients grows on average as ( ǫ /ν) 1/2 . However, the fluctuations of velocity gradients are orders of magnitude larger than this average value, a phenomenon referred to as small-scale intermittency [3,4]. Such extreme events play a crucial role in numerous physical processes in both nature and engineering, e.g. turbulent dispersion [5], cloud physics [6], turbulent combustion in jet engines [7,8], and are also conjectured to be connected to regularity and smoothness of fluid equations [9,10]. Hence, understanding their formation and statistical properties is of central importance in developing a complete theory of turbulence [4]. The complexity of the problem is apparent in Fig. 1, which shows the structure the velocity gradients. The strongly intermittent nature of turbulence is clearly visible by the highly inhomogeneous distribution of the regions of very intense gradients (see also [2]). * dhawal.buaria@ds.mpg.de Fluid turbulence involves a wide range of spatial scales, from approximately the system size all the way down to the very finest scale, corresponding to the largest gradients. In this respect, it can be viewed as an emblematic example for other complex dynamical systems, where such extreme events are also observed [11,12], including the climate system [13], with its far-reaching implications.
Ever since Kolmogorov formulated and refined his seminal hypotheses [14], intermittency in turbulence has been the subject of many studies [4]. In particular, detailed investigations demonstrate that the very large fluctuations in velocity gradients become more extreme with increasing Reynolds number [15,16]. While there have been theoretical proposals to describe quantitatively the Reynolds number dependence of velocity gradient fluctuations [17][18][19][20], they have remained difficult to verify due to lack of reliable data. In fact, directly measuring the most intense fluctuations, experimentally or numerically, over a reasonable range of Reynolds number is a very challenging endeavor, as very high spatial and temporal resolution is required to accurately resolve such fluctuations. As recently pointed out [21], this demand can be even stricter in numerical simulations than previously expected. Consequently, such high resolution investigations have been so far restricted to low Reynolds numbers [16,23].
In this work, we characterize the dependence of the extreme velocity gradients on the Reynolds number, and illuminate the underlying physical processes. To this end, we use high resolution direct numerical simulations (DNS) of isotropic turbulence, based on highly accurate Fourier pseudo-spectral methods. To accurately resolve the extreme gradients, all our simulations were carried out with a small-scale resolution at least 3-4 times higher than typical turbulence simulations, along with appropriate temporal resolution [21]. Going up to grids of 8192 3 points, we have obtained results at Taylor-scale Reynolds number (R λ ) ranging from 140 to 650. In order to characterize the gradients, we first consider the probability density functions (PDFs) of square of the norm of strain and vorticity, which represent the symmetric and skew-symmetric components of the velocity gradient tensor respectively. They are analogous to dissipation rate and enstrophy and have the same mean values (within a prefactor) in isotropic turbulence, given by 1/τ 2 K , where τ K = (ν/ ǫ ) 1/2 is the Kolmogorov time scale [16]. Consistent with previous works [15,22], we observe that the PDFs of these quantities, when normalized by τ K exhibit tails that become broader with increasing R λ . By further characterizing these PDFs, we demon-strate that their tails can be collapsed very well over the range of R λ covered here, when instead normalized by the time scale: which implies that the strongest gradients in the flow correspond to a time scale τ ext which increasingly decreases with respect to τ K as R λ increases (and hence the strongest gradients in the flow grow as τ −1 K R λ β ). Numerically, we find that β ≈ 0.775 ± 0.025. The tails of the PDFs of velocity increments δu r , normalized by the Kolmogorov velocity scale u K (= (ν ǫ ) 1/4 ), over distances r ≤ η (where η = (ν 3 / ǫ ) 1/4 is the Kolmogorov length scale), also become broader when R λ increases. On the contrary, when normalized by the r.m.s. of velocity fluctuations u ′ , the tails grow very slowly. With the understanding that the most intense gradients in the flow occur with velocity increments of order u ′ over a scale η ext we conclude that η ext ∼ ηR λ −α , with α = β − 1 /2, represents the smallest scale in the flow. The collapse of the tails of PDFs of δu r , when normalized by either u K R λ β or u ′ R λ α , supports these findings. Comparisons with existing theoretical predictions [17,20], point to difficulties in explaining our data. However, these theories utilize the phenomenological definition that the smallest scales in the flow correspond to a local Reynolds number of unity [3], which, contrary to the numerical results of [24,25] and also our own, does not appear to be satisfied at the location of intense gradients, as utilized in current work to characterize the smallest scales.
Consistent with earlier works [24,26,27], we find that the structures corresponding to the largest velocity gradients to be vortex tubes. We do not find extreme events in strain and vorticity to be colocated [16,28]. Conditional averaging shows that intense strain is always likely to be accompanied by equally intense vorticity. However, intense vorticity is found to be accompanied by relatively less intense strain, with an approximate power law dependence corresponding to exponent γ < 1, which very slowly increases with R λ . With the interpretation that η ext is the radius of most intense vortex tubes, we use simple scaling arguments to relate it to the conditional strain, and thereby relate γ to β. This suggests that β would also slowly increase with R λ . We conjecture that the limit β = 1 (and α = 0.5), as predicted by some intermittency theories, would only be realized for R λ → ∞.
The rest of the manuscript is organized as follows. In Section II, we describe our numerical methods. Our numerical results concerning the scaling of extreme velocity gradients are presented in Section III. The structure of regions of very intense velocity gradients is investigated in Section IV. Section V contains a discussion, comparing our results with existing theories, and then providing an an alternative description connected to flow structure examined in Section IV. We briefly discuss the implications of our results on future DNS and experiments in Section VI. Finally, we present our conclusions in Section VII.

II. NUMERICAL APPROACH AND DATABASE
The present work is based on DNS of the incompressible Navier-Stokes equations where u is the velocity field (satisfying ∇ · u = 0), p is pressure, and f is the forcing term used to maintain a stationary state [29,30]. The equations are solved utilizing a massively parallel implementation of Rogallo's pseudo-spectral algorithm [31], whereby the aliasing errors are controlled by a combination of truncation and phase-shifting [32]. We use explicit second-order Runge-Kutta scheme for time integration, with the time step ∆t subject to a constraint for numerical stability expressed in terms of the Courant number, C = ∆t ∆x (||u|| 1 ) max , where || · || 1 represents the L 1 -norm and the maximum is taken over all (N 3 ) grid points. The flow simulated is homogeneous and isotropic with periodic boundary conditions, on a cubic domain of (2π) 3 for all cases.
As stressed earlier, appropriate numerical resolution of the small scales is crucial to our study of extreme velocity gradients. Spatial resolution in pseudo-spectral DNS is typically measured by the parameter k max η, where k max = √ 2N/3 is the largest wavenumber resolved and η is the Kolmogorov length scale. Equivalently, one can use the ratio ∆x/η (≈ 2.96/k max η), where ∆x = 2π/N is the grid spacing. Most turbulence simulations, aimed at reaching high Reynolds number, are in the range 1 ≤ k max η ≤ 2 [28,33]. However, resolution studies have shown that such a resolution is inadequate for studying extreme events in velocity gradients [16,21,23]. Hence, we have consistently used k max η ≈ 6 in all the runs shown here. Additionally, we have also used a Courant number of 0.3, instead of 0.6 in previous studies e.g. [16,22,28], as it was recently found that the latter led to spurious over-prediction of the gradients [21]. Resolution studies presented in [21] and our own tests confirm that the resolution used here is adequate to address the questions asked in this work. The database used here and the corresponding simulation parameters are listed in Table I. The Taylorscale Reynolds numbers (R λ ) considered here are similar to those in some previous works [16,22], but with a much higher small-scale resolution as emphasized earlier. These high resolution simulations were recently used in [21]. In the present work, we simply restarted these runs (which were already in a stationary state) and extended them to substantially longer times to greatly improve statistical convergence. We list the length of the current simulation T in terms of the large-eddy turnover time T E . The statistical results shown here were obtained by analyzing N s instantaneous snapshots for each run. Whereas a long simulation is typically desirable for sampling accuracy, finite resources have limited the value of T for the highest resolution runs. Nevertheless, since our focus is on highly intermittent velocity gradients, one can improve sampling by simply analyzing more snapshots for a given simulation length as the Reynolds number increases. This is justified, both by the increase of the ratio of time scales T E /τ K with R λ , and also by the increasingly smaller time scales associated with the extreme events, as discussed in the manuscript.

A. PDFs of vorticity and strain
In order to study small-scale intermittency, in this subsection, we characterize the velocity gradient tensor by its two quadratic invariants [4], namely Ω = ω i ω i , where ω = ∇ × u is the vorticity, and Σ = 2s ij s ij , where s ij is the strain rate tensor defined as s ij = (∂u i /∂x j +∂u j /∂x i )/2. The former is the enstrophy and the latter is simply the dissipation divided by viscosity i.e. Σ = ǫ/ν. In isotropic turbulence, as considered here, Ω = Σ = 1/τ 2 K , where τ K is the Kolmogorov time scale. Investigating the extreme events in Ω or Σ amounts to focusing on the outmost parts of the (wide) tails of their PDFs. Issues of statistical convergence makes a precise determination of these quantities extremely difficult. With the available data, we estimated the statistical error in each bin, by looking at the sample to sample fluctuations across various snapshots used to determine the statistics. We kept only bins with an error less than 20% compared to the mean. These PDFs, converged with respect to both the small-scale resolution [21] and statistical sampling, allowed us to determine the properties of the extreme events, as presented next. Fig. 2 shows the PDF of (a) Ω and (b) Σ, normalized by their mean value, 1/τ 2 K , at various Reynolds numbers. We primarily observe that as the Reynolds number increases, the tails of these PDFs get wider and extend to much higher values. Equivalently, the likelihood of finding a value of Ωτ 2 K or Στ 2 K larger than a given large value increases with the Reynolds number. This is expected and consistent with previous studies. One notices, however, that the part of the PDFs, corresponding to events smaller than about 10 times the mean, appear to approximately collapse for different Reynolds numbers. This can be seen in the insets of Fig. 2, which show a zoomed in version.
The existence of increasingly large fluctuations, as shown in Fig. 2, leads us to ask how large are the extreme gradients and how quickly do they grow with increasing Reynolds number. We propose to answer this by rescaling the PDFs. Namely, we use a different time (1)), to rescale the extreme values. Denoting f Ω (Ω e ) and f Σ (Σ e ) as PDFs of Ω e = Ωτ 2 ext and Σ e = Στ 2 ext respectively, Fig. 3 shows R λ δ f Ω (Ω e ), and R λ δ f Σ (Σ e ). The factor R λ δ provides a measure of how rare the largest fluctuations of Ω e or Σ e are, when R λ increases. As shown in Fig. 3, using β ≈ 0.775 and δ ≈ 4.0, the wide tails of rescaled PDFs are almost perfectly collapsed. This indicates that while the average events in Ω and Σ scale as τ −2 K , the most extreme events behave like τ −2 K R λ 2β . The exponents β and δ, used in Fig. 3 to collapse the large tails of the PDFs, can be also empirically determined by utilizing a functional form of the tails of PDFs of Ω and Σ. While theories have proposed several functional forms for the entire range of PDFs [4,[34][35][36], the stretched exponential function is known to empirically fit the tails of the PDFs very accurately [15,16,21,37,38]. Since the tails of the PDFs in Fig. 3 collapse, we use the following stretched exponential functional form to mathematically verify the value of β: where x = Ωτ 2 K or Στ 2 K (or alternatively Ω/ Ω and ǫ/ ǫ respectively in the notation of [16,21]) and a, b, c are the fitting parameters. Applying a change of variable The collapse shown in Fig. 3 implies that such that the constants b ′ and a ′ are independent of R λ . Thus, the dependence of b 1/c as a function of R λ , provides a direct access to β.
To determine the coefficients, we simply fit the logarithm of the PDF to the functional form (log a − bx c ) of Eq. (3). We choose the fitting window to be x ≥ 50, which sufficiently excludes the region around the mean value, where the PDFs appear to collapse for various R λ (as shown in insets of Fig. 2). We also explicitly checked by extending the fitting range to smaller values, but found that the results remained virtually unchanged. The determination of the three parameters a, b and c then leads to a non-linear regression. However, since non-linear regression can be very sensitive to the initial guess values for the fitting parameters -especially the value of the exponent c in this particular case -determining the parameters directly in such a manner can result in significant error [39]. On the other hand, if the value of c is known beforehand, then a very robust fit can be obtained, since the fitting procedure reduces to a linear regression to determine only a and b. The exact values of a and b would obviously substantially differ for different values of c, but this would not matter if they all provide the same value of β (which as shown next, is the case).
The values of the exponent c in previous numerical studies [16,21] were found to be close to the range 0.23 − 0.25, with a possible scatter within 0.19 − 0.29 and no clear dependence on Reynolds number (e.g. see Table 4 of [16]). Keeping this in mind, we therefore fit the PDF by assuming fixed values of c, ranging from 0.19 to 0.29 in increments of 0.02, and determine the parameters b and a for PDFs Ω and Σ for all available R λ . Note that a wider range for c may be considered, but this chosen range falls within the error obtained from a naive nonlinear regression and hence for values outside the chosen range, the quality of fit starts deteriorating. To provide a measure, for the chosen values of c, the coefficient of determination (R 2 ) was greater than 0.995 for every fit. Additionally for each c, the resulting values of a and b are always obtained with greater than 95% confidence, resulting in negligible error bars. In fact, these values are even found to be quite insensitive to minor variations in the fitting window, e.g., our fits compare extremely well with those of [16], who considered a fitting window of 5 ≤ x ≤ 100 for R λ ≤ 240. In this regard, we make a note that the results of [16] can only be trusted for R λ ≤ 240, since the higher R λ runs were affected by resolution issues, as reported in [21]. Nevertheless, the excellent quality of fit is evident in Fig. 3 (also see Fig.9 of [16] which is also in near perfect agreement with our fits). The dependence of the coefficient b 1/c on R λ at various values of c is shown in Fig. 4. The data points correspond to curve fits for both Ω and Σ (thus giving two sets of points for each c). For the sake of clarity, the values of b 1/c are divided by their corresponding values at R λ = 650, which imposes that all the curves shown in Fig. 4 pass through 1 at R λ = 650 (since higher R λ provides a larger fitting range and hence can be expected to be most robust). This also allows us to directly compare the data points for every c value considered. We find that all sets of points superpose reasonably well, and remarkably point to a similar power law in R λ . This collapse demonstrates that the determination of the exponent β is not very sensitive to the precise value of c, at least with the available data. However, weak deviations from scaling cannot be ruled out, especially if an even larger range of R λ is considered in future. In fact, we will later (in Section V) present arguments supporting a very weak growth of β with R λ .
We would like to further clarify that ideally the best curve fits to Eq. (3) may lead to a dependence of c on the Reynolds number R λ . However, choosing a fixed c substantially improves the quality of curve fit and also minimizes the sensitivity to the fitting range. In fact, a fixed value of c also helps in determining the scaling without ambiguity, since if c is a function of R λ , the constant b ′ in Eq. (5) will also become a function of R λ , which would recursively require additional non-linear fits to obtain β. The minor deviations for different c values at various R λ , can also be possibly explained by this. Nevertheless, the good collapse seen for such a wide range of c values provides clear evidence that the scaling proposed provides a compelling description of our data, at least over the available range of R λ . We will see that this is also further supported by results shown in section III B. Finally, by fitting a power law through the obtained data points (marked by dashed line), we obtain β = 0.775 ± 0.025 (where the error bar takes into account the variation across different c values), which was used in scaling the PDFs in Fig. 3. The same procedure for the parameter a (not shown) gives δ ≈ 4.0, with deviations of approximately 5-10%. Thus, systematically characterizing the PDFs of Ω and Σ, we are able to mathematically determine that the strongest gradients in the flow grow as τ −1 K R λ β , with β = 0.775 ± 0.025. We again emphasize that obtaining such a result required statistically well converged PDFs (in turn requiring adequate spatial and temporal resolutions) over a wide enough range of Reynolds numbers.

B. PDFs of velocity increments
In order to further validate the scaling obtained in Section III A, we next investigate the PDFs of velocity increments. In simplified notation, velocity increments are given as δu r = u(x + r) − u(x), where the separation distance r can be either in the direction of u (longitudinal increments) or perpendicular to u (transverse increments). Over very small distances, the velocity differences essentially reduces to the velocity gradients (within a constant factor), aside from systematic but small errors introduced by use of finite differencing. Hence, we can expect velocity increments over small distances to show the same scaling as derived earlier. However, to extract information about the gradients, one still needs to ensure that r is sufficiently small. We note in this respect that the high resolution of our runs, k max η ≈ 6, effectively allows us to calculate increments over a very small distance (r ≈ η/2). A benefit of using velocity increments is that their 1D surrogates can be also obtained and verified using experiments [40]. Fig. 5a shows the PDF of δu r , normalized by the Kolmogorov velocity scale u K , for r/η ≈ 0.5, at various Reynolds numbers. Both the longitudinal and transverse components are shown (in dashed and solid lines respectively). Since we are interested only in the magnitude of the increments, we take the absolute value of δu r . Consistent with the results shown in Fig. 2, the initial part of the PDFs, corresponding to moderate events, superpose very well (see inset of Fig. 5a) and as Reynolds number increases, the tails start growing. Both the longitudinal and transverse increments show the same behavior, with the transverse component being expectedly larger [3].
To further characterize the velocity increments, we next consider the PDFs of δu r /u ′ , where u ′ is the r.m.s. of velocity fluctuations. Fig. 5b shows the PDFs of δu r /u ′ corresponding to r = η/2, at various Reynolds numbers and for both longitudinal and transverse components. The corresponding PDFs for r = η are shown in the inset. Even at such a small separations, we observe that the velocity differences can be as high as u ′ , for the entire range of Reynolds number considered. This appears to be consistent with the observation of [24,25] at R λ 170. Additionally, it appears that while the probability density for a given δu r /u ′ decreases with increasing R λ , the extent of δu r /u ′ itself slowly increases.
Following similar ideas as in Fig. 3b, we determine the PDFs of δu r /u K , rescaled by R λ −β , shifted by a factor R λ δ . The result is shown in Fig. 5c corresponding to r/η ≈ 0.5 for both longitudinal and transverse components. Note for the PDFs of δu r /u ′ this corresponds to rescaling by R λ −α , with α = β−0.5, since u ′ /u K ∼ R λ 1/2 [3] (the importance of α is discussed later in Section V). We find that the rescaled PDFs collapse very well for different Reynolds numbers. The scatter towards the very end of the tails (especially for the transverse component) can be attributed to lack of statistical convergence for the endmost bins. In the inset of Fig. 5c, we repeat the exercise, but now for PDFs corresponding to r/η ≈ 1. The superposition, although comparatively worse with respect to r/η = 0.5, still remains very good. An alternative approach to investigate velocity increments is to consider the quantity δu r /r, which for sufficiently small r, is a proxy for the velocity gradient and thus also independent of r. Hence, it is tempting to use the same scaling based on τ ext as before, to collapse the tails of various PDFs of δu r /r (as done in Fig. 5c). Fig. 6 shows such rescaled PDFs for several values of r/η, at R λ = 140 and 650. Once again, we find that the tails of the curves for r/η ≈ 0.5 collapse very well for both R λ . However at r/η ≈ 1, the curve for R λ = 650 starts  to deviate from this collapse. At r/η ≈ 2, the curve for R λ = 650 significantly deviates from that for R λ = 140.
This provides yet further evidence that the resolution of ∆x/η ≈ 0.5 is adequate for both R λ . On the other hand, ∆x/η ≈ 1, while adequate for R λ = 140, is insufficient for R λ = 650. Additionally, we also see that as r/η grows, the deviations of the curves at R λ = 650 increase faster than those at R λ = 140. This also provides a hint that the smallest length scale in the flow is actually smaller than η and additionally decreasing with increasing R λ . This observation will be further analyzed and discussed in Section V.

IV. STRUCTURE OF REGIONS OF INTENSE VORTICITY AND STRAIN
In order to gain some understanding on the structure of regions of intense gradients, we first use flow visualization. Vorticity arranged in tube-like structures has been repeatedly seen in DNS over a very large range of Reynolds numbers -from very low at R λ ≈ 45 [26] all the way to R λ ≈ 1100 [27] (although the small-scale resolution in these simulations was limited). Whether such tubes carry the most intense regions of velocity gradients in the flow, however, has been questioned by a recent study [28], which suggested that the largest values of Ω and Σ appear colocated and without any coherent structure. It is important to note that these observations may have been affected by the numerical artefacts documented in [21]. One of the motivations of the present work is to revisit the issue. We again stress that the present visualizations are based on DNS at much higher spatial resolution than previously available. Fig. 7 shows a collection of instantaneous snapshots from the R λ = 650 run, focusing on the region of most intense gradients. Since vorticity is in general larger than strain, the domains are chosen such their centers correspond to the maximum value of Ω. The various panels show different contour thresholds (indicated as C in subfigure caption) in cyan for Ωτ 2 K and in red for Στ 2 K . In the first panel (Fig. 7a), a domain of 301 3 grid points or (150η) 3 is shown with the contour threshold of 50 for both vorticity and strain. Structures consisting of clusters of vortex tubes, qualitatively similar to e.g. [27], are readily seen. The spacing between neighboring tubes widely varies: some vortices are relatively isolated, others seem to be more strongly interacting with their surrounding. Large values of strain are mostly located around large vortices, a phenomenon noticed many times (see [28] and references therein).
Panel (b) zooms into the region of most intense gradients, showing a domain of (50η) 3 with a contour thresh-old of 300. The structure is composed of two closely interacting vortex tubes, wrapped around by intense strain. In panels (c) and (d), contour levels are successively increased to 500 and 1500 respectively and we also further zoom in to show a domain of (25η) 3 in (d). The vortex tube structure becomes very distinct, whereas the region occupied by strain reduces substantially in (c) and completely disappears in (d). Note, the largest value of Ωτ 2 K is equal to about 3000 (at the center of the domain in each panel). In comparison, the largest value of Στ 2 K is about 1800, located in top left corner of domain shown in (c) -and hence no coherent strain region is visible in (d). Although not explicitly shown here, we also confirmed that the velocity increments around the center of each panel Fig. 7 correspond to far tails of PDF of δu r as shown in Fig. 5, i.e., δu r ≃ u ′ .
We analyzed many such flow fields corresponding to different snapshots and virtually all of them show a qualitatively similar behavior, i.e., as the contour thresholds are increased tube-like vorticity structures become prominent and high-strain regions shrink and disappear at lower values than the high-vorticity regions. While not directly evident in Fig. 7, we also find that the locations of maximum values of vorticity and strain are typically separated by at least 10 − 20η and never coincident, e.g. in Fig. 7c, These observations confirm that the largest values of Σ are much smaller than the large values of Ω and the regions for large values of Σ and Ω are not co-located. Hence, we conclude that visualizations in [16,28] were also affected by resolution issues reported in [21]. We remark in this respect that new independent tests at R λ = 1300 and k max η = 3, with a time step twice smaller than in [28] -although not shown here -confirm that the qualitative aspect of the regions of extreme vorticity/strain are similar to that shown in Fig. 7.
In order to quantify the relation between strain and vorticity, we next consider their conditional expectations with respect to each other -shown in Fig. 8 for various R λ . For low values of Ω or Σ, the conditional dependencies are very weak, i.e., strain and vorticity appear to be decorrelated. However, for conditional values greater than unity, i.e., the mean value, the conditional expectations clearly increase, seemingly showing a power law. Comparison with a dashed line of slope 1 (on log-log coordinates), suggests that Ω|Σ ∼ Σ 1 . In contrast This implies that intense events in strain are always likely to be accompanied by equally strong events in vorticity, whereas the strain is comparatively weaker in very intense vortices. This appears to be consistent with the earlier observations of vorticity being more intermittent than strain [16,26,41] and ultimately concerns with the inter-relationship of vorticity and strain, which is still an open question in turbulence. Note that [16] shows a similar plot as Fig. 8, however their curve for Σ|Ω spuriously approaches a slope of 1 (for large Ω) because of resolution issues [21].
Interestingly, Fig. 8b also suggests that the exponent γ slowly increases with R λ . By fitting approximate power laws, we find that γ varies from 0.60 − 0.72 over the range of R λ considered here (see inset of Fig. 8b), although the variation appears to get weaker as R λ increases. This naturally leads to the question of what the limit of R λ → ∞ entails, which our data is unable to answer conclusively. Theoretical considerations suggest that γ = 1 in the large R λ limit [4,42,43]. Given the very slow increasing trend of γ, it is evident that extremely high R λ would be necessary to realize γ = 1, if at all possible (a simple sigmoidal or power law extrapolation suggests γ = 0.99 would be realized for R λ 20000). Therefore, the differences between strain and vorticity are expected to persist, even at the highest turbulence levels on earth. A fundamental understanding of Eq. (6) from first principles, i.e. a determination of the strain acting on a given vortex, resulting from the tangle of vortices as shown in Fig. 7 is still an open question in turbulence. As we will suggest in Section V B, the power law dependence on the strain conditioned on vorticity in fact provides a way to understand the scaling exponent β.

V. THEORETICAL CONSIDERATIONS
In this section, we discuss the observation that the extreme velocity gradient fluctuations scale as τ −1 ext = τ −1 K × R λ β . We first compare our result with existing theories, especially the multifractal model, and thereafter provide a new description for the observed scaling which relates to the structure of the flow discussed earlier in Section IV.
In simplest terms, extreme velocity gradients result from large velocity differences over a very small length scale; the largest velocity gradient in the flow can be written as proportional to δu max /η ext , where δu max is largest velocity difference over the smallest length scale η ext [19], which given the flow structure, can be physically interpreted as the radius of the smallest vortex tube [24]. Our notation that the largest velocity gradients scale as τ −1 ext therefore implies τ ext ∼ η ext /δu max .
Thus, the question how large the gradients can grow entails answering how large δu can become over the smallest length scale η ext . Based on earlier resolution studies [16,21] and also on the results presented in Section III B, the smallest scale η ext can be defined by the resolution at which the PDFs of the gradients have converged -which for the present range of R λ , gives η/2 ≤ η ext η. Notice that one could formally define scales smaller than η ext , but given a smooth velocity field, the velocity increments at such scales will simply decrease linearly with the scale size, with respect to those at η ext (as also demonstrated by the velocity increment PDFs for r/η = 0.5 and 1 at R λ = 140 in Fig. 6). Thus, any length scale smaller than η ext would show the same scaling as η ext itself and would be immaterial for the purpose of present study.

A. Comparisons with existing theories
It is natural to interpret our results using existing theories. To this end, we begin by reviewing the multifractal model, which provides explicit predictions concerning the smallest scales in the flow. In the multifractal model, as well as in some other phenomenological approaches, a recurring concept is that of the fluctuating local viscous cutoff scale, say η x , defined such that the velocity increment over this scale has a local Reynolds number of unity [17,20]: which essentially results from equating the viscous time scale η 2 x /ν to the convective time scale η x /δu. In the multifractal framework, the velocity increment over a distance r is given as δu r /u ′ ∼ (r/L) h , where L is the energy-injection scale, and h is the local Hölder exponent within an interval [h min , h max ] such that a fractal set D(h) can be determined for every h. Thereafter, following the derivation of [17], the smallest scale in the flow can readily be obtained corresponding to the minimum Hölder exponent where η is the Kolmogorov length scale. It also follows which implies β = 2α. Earlier works have suggested that h min = 0 [17,19], which gives α = 0.5, β = 2α = 1 and hence δu max ∼ u ′ . The value of β ≈ 0.775 derived earlier can be obtained by using h min ≈ 0.06, and would additionally imply α = β/2 ≈ 0.39. However, h min ≈ 0.06, or rather any non-zero positive value of h min , implies α < 0.5, and hence suggests that δu max /u ′ would decrease with increasing R λ . In fact, h min > 0 also suggests that the range of δu r /u ′ decreases with R λ for a fixed r/η. Our observation in Fig. 5b, which demonstrates that the range of δu r /u ′ does not show any sign of decreasing with R λ at r/η 1 -and rather appears to be slowly increasing, does not unambiguously support the decay of δu max /u ′ implied by the theory. At the same time, since η/2 < η ext η over the present range of R λ , the extent of PDFs in Fig. 5b implies δu max u ′ for the strongest gradients. We notice that the probability density for δu u ′ appears to decrease slowly with R λ , however, the excellent superposition of the PDFs in Fig. 5c indicates that the decay is at best algebraic (decreasing as R λ −δ ), and therefore, the probability should remain finite even as R λ → ∞. The above suggests h min = 0 to allow δu max ∼ u ′ , leading to β = 1 (for α = 0.5) as suggested by [17,19], but at odds with the observed value of β = 0.775 (and the corresponding h min ≈ 0.06).
In our view, the inconsistency above is a result of the assumptions built into the extension of the multifractal theory, originally developed to describe the inertial scales, to far dissipative scales. In particular, the definition in Eq. (8) obtained by equating the convective and dissipative time scales, while reasonable for inertial range (where the rate of energy transfer across scales can be assumed to be constant), does not appear justified at smallest scales, where dissipation dominates. This is readily observed in Fig. 5a, where the local Reynolds number from the tails (corresponding to strongest gradients residing in vortex tubes) increases steadily with the R λ (much strongly than η ext decreases with R λ ). In fact, earlier works based on DNS at R λ 170, have already have shown that the local Reynolds numbers corresponding to the vortex tubes where the intense gradients are localized are much larger than unity, and appear to scale differently from the multifractal prediction [24,25]. Our numerical results at significantly higher R λ (and also higher small-scale resolution) further strengthen this conclusion and puts into question the (phenomenological) criterion that the smallest scale in the flow can be determined by a local Reynolds number of order unity and hence, also the relation α = β/2. However, an additional remark is necessary in this context. While we associate the smallest scales of motion with the extreme gradients (which appear to reside in vortex tubes), the theoretical constructs resulting from multifractal considerations are based on the local scaling, where the smallest scales simply correspond to the minimum Hölder exponent h min , without any explicit connections to the flow structure. As a result, there is no guarantee that the scales resulting from h min actually correspond to the structures observed in Fig. 7, calling for some caution when comparing our results with the multifractal theory.
An alternative description, which also utilizes Eq. (8), is that of Yakhot and Sreenivasan [20]. In their approach, the even moments of the velocity increment (δu) 2n (or the structure functions), each correspond to a unique dissipative length scale η n , such that the smallest possible scale in the flow corresponds to n → ∞. Thereafter, again utilizing Eq. (8), η n can be related to the anomalous inertial range scaling exponents of structure functions and n → ∞ results in a similar prediction as that of multifractal theory with h min = 0, i.e., α = 0.5. However, once again, such an approach does not appear as justified given the lack of evidence for Eq. (8) to define the smallest scales. In fact, numerical results of [16,23,44] all suggest that the smallest scales in fact grow weaker than the prediction of Yakhot-Sreenivasan theory (and hence also that of multifractal theory for h min = 0). This observation is once again reinforced by the results presented in current work.

B. Alternative description in light of strain-vorticity dynamics
In view of the apparent shortcomings of intermittency theories reviewed in the previous subsection, we propose, in order to reconcile our observation concerning the very large velocity differences and the exponent β < 1, a different description that directly relates to the flow structure explored in Section IV.
Based on Fig. 5b, we propose that the strongest gradients correspond to δu max ∼ u ′ over the smallest scale η ext . This is in line with observations of [24,25] and also h min = 0 as postulated by [17]. Note, this assumption also essentially implies that the PDFs of velocity (and hence velocity increments) are bounded [9]. Thereafter, substituting δu max ∼ u ′ and τ ext = τ K R λ −β into Eq. (7) gives where we have used u ′ /u K ∼ R λ 1/2 from classical scaling estimate [3]. The value β = 0.775 ± 0.025 found numerically leads to α = 0.275 ± 0.025, which is also the value used in Fig. 5c to collapse the PDFs of δu r /u ′ for r = η/2 < η ext . Note, the convergence of δu r /r to the velocity gradient ensures that the PDFs have a well-defined limit for r ≤ η ext (and hence the PDFs at r < η ext can be simply obtained by linearly rescaling the PDF at r = η ext by a factor r/η ext , which essentially is the same as R λ α for r = η/2. This value of α ≈ 0.275 can be further verified by considering the PDFs shown in Fig. 6. While for r ≤ η ext , δu r /r converges to the velocity gradient, systematic deviations arise for r > η ext . These deviations from the gradient can be accordingly quantified by analyzing the higher order terms in a Taylor series expansion of δu r /r and can be shown to be approximately proportional to r/η ext [16,20]. Since η ext /η decreases when R λ increases, at a given value of r/η, the deviations from the PDFs, especially in the tail, from their limiting form at r/η ext ≪ 1 also increases, as clearly seen in Fig. 6. In addition, we find that the deviations of the PDF tails at fixed values of r/η ext to be independent of R λ . For α ≈ 0.275, η ext /η decreases by approximately 1.53 between R λ = 140 and 650. In contrast, taking the value of α predicted by the multifractal theory, α = β/2 ≈ 0.39, leads to a variation of 1.82 in the ratio η ext /η. The r/η values at R λ = 140 shown in Fig. 6 are chosen, as close as possible, within these factors (of 1.53 and 1.82), compared to the r/η values for R λ = 650 (the curves corresponding to 1.53 are shown in dashed-dotted lines, whereas curves corresponding to 1.82 are shown in dashed lines). As visible, the PDFs corresponding to the factor of 1.53 between the two R λ cases collapse remarkably well (especially as r/η increases), hence providing an alternative means to verify α ≈ 0.275.
Interestingly, alongside some scaling arguments to evaluate α, Eq. (12) leads to two different limits (which incidentally also correspond to previously reported cases in literature). The first limit corresponds to simply assuming that the smallest length scale in the flow is the Kolmogorov length scale, i.e., η ext = η. Using this, we get α = 0 and β = 0.5, and hence This result was derived in [24,25], based on the analysis of DNS data at relatively low Reynolds numbers (R λ 170), which the present work greatly improves upon. The second limit consists in taking into account the extreme fluctuations of the velocity gradients. Physically, the smallest scale in a flow can be thought to result from a balance between viscosity ν and strain Σ (= 2s ij s ij , as defined earlier), which leads to the expression of the length scale: η ext ≃ (ν 2 /Σ) 1/4 , familiar in a number of contexts [45]. Assuming η ext = η, which leads to Eq. (13), amounts to a mean field approximation, consisting in replacing the strain by its averaged value (as η is calculated from the mean dissipation). Taking into account the large fluctuations of Σ results in η ext being smaller than η [18]. In this regard, the second limit can be simply derived by evaluating η ext using the maximum value of strain (Σ max ), i.e., η ext = (ν 2 /Σ max ) 1/4 . Thereafter, using Σ max ∼ τ −2 ext and δu max ∼ u ′ based on earlier results, it follows from Eq. (7) where we have used ν/u 2 K = τ K and u ′ /u K ∼ R λ 1/2 . This implies β = 1 and α = 0.5 from Eq. (12), as also predicted by intermittency models discussed earlier [17,19,20]. However, this is not completely surprising, as defining η ext based on Σ max with Σ max ∼ ντ −2 ext also leads to β = 2α, which is essentially the multifractal prediction, and in conjunction with β = α + 0.5 derived in Eq. (12) gives β = 1 (and α = 0.5). Additionally for this scenario, the local Reynolds number, which can be written as η 2 ext ν −1 τ −1 ext using Eq. (7), comes out to be constant as inherently assumed in intermittency theories discussed earlier.
The numerically observed value of β ≈ 0.775 lies between β = 0.5 and 1, which suggests that η ext results from a strain, intermediate between the two limits considered before. In fact, this is precisely what we observed in Section IV. As noted earlier, Eq. (6) (and Fig. 8b) suggests that the strain acting on a very intense vortex tube is significantly weaker than naively expected by postulating Σ ∝ Ω. A simplified estimate consists in substituting Σ max in the argument leading to Eq. (14) by τ −2 K (τ 2 K Ω max ) γ , as suggested by Eq. (6). Thereafter, we get The limits of β = 0.5 and 1 correspond to γ = 0 and 1 respectively. In view of the weak dependence of γ shown in the inset of Fig. 8b, Eq. (15) suggests a dependence of β on R λ . The values of γ observed over the range of R λ studied here, 0.60 γ 0.72, implies a variation of β in the range: 0.72 β 0.78, which is quantitatively consistent with β ≈ 0.775 determined empirically in Section III. The weak variation of β implied by Eq. (15) may also explain the slight deviations from scaling seen for R λ = 140 in Fig. 4. In fact in Fig. 4, considering only data points at R λ = 140 and 240, the slope corresponds to β ≈ 0.73, which appears to be remarkably consistent with that obtained from γ for these R λ . On the other hand, the interesting possibility that γ → 1 when R λ → ∞ would then suggest, in view of Eq. (15), that β → 1, as originally expected by some theories (albeit corresponding to a constant local Reynolds number much larger than unity). However, the very slow variation of γ shown the inset of Fig. 8b would indicate that β = 1 would be attained at extremely large values of R λ , likely larger than practically relevant. In this regard, the scope of existing predictions in understanding finite Reynolds number scaling appears to be severely limited. Whereas the prediction of the exponent γ and its dependence on R λ is a very challenging task, we briefly note that the cascade model of She and Leveque [46] presents a similar idea, though with shortcomings. The model postulates that the locally averaged dissipation field ǫ r at a scale r and the corresponding moment ratios: ǫ (p) r = ǫ p+1 r / ǫ p r , are related to the hierarchy of complex structures in the flow. The most singular structures correspond to ǫ (∞) r , which in the phenomenology of [46] obeys the power law dependence: ǫ (∞) r ≃ ǫ (L/r) µ , with µ = 2/3. While their original arguments were postulated for inertial scales, if one were to extend the cascading process down to smallest scale, i.e., r = η ext = ν 3 /ǫ , it leads to α = (3µ)/(8 − 2µ). Using µ = 2/3 as proposed by She-Leveque then gives α = 0.3, which is close to our current prediction of α = 0.275 ± 0.025. However, they also suggest h min = 1/9 within the multifractal formalism, which using Eqs. (9)-(11), gives α = 0.3, β = 0.6 and δu max ∼ u ′ R λ −0.2 , which are clearly inconsistent with our data. Ultimately, phenomenological descriptions (as those of [46]) are at best weakly connected to flow structures, and typically assume a constant value of the exponents such as β and α, thus ignoring any possible dependence on R λ , as suggested from Eq. (15) and Fig. 8. Hence, it appears that the closeness of α between the She-Leveque model and our current result is only fortuitous.
In conclusion, our analysis of the most intense vortex structures observed in the flow relates the exponent β with the properties of the strain acting on vortices, and in particular with the exponent γ defined by Eq. (6). The weak variation of γ with R λ , see Fig. 8, implies that β should increase with R λ . A very natural conjecture is that the symmetry between strain and vorticity, clearly broken at finite R λ , will be restored as R λ → ∞, and that the exponents γ and β both tend to 1, corresponding to earlier predictions [17,20]. Understanding the R λdependence of the strain acting on intense vortex tubes appears as an essential question in this regard, that deserves renewed theoretical attention.

VI. IMPLICATIONS FOR SIMULATIONS AND EXPERIMENTS
The identification of the smallest scale η ext , characteristic of the largest velocity gradients in the flow, which decreases faster than η when R λ increases, has some obvious consequences for the resolution constraints required in both DNS and experiments. In DNS, it is typical for most studies based to be performed with a k max η or ∆x/η held constant across the range of R λ simulated (e.g. see [33,47]). On the other hand, in experiments, the resolution, determined by the probe size or the data acquisition frequency, often gets worse as R λ increases [40]. The present results, however, show that in studies focused on intermittency, one must continuously improve ∆x/η as R λ is increased to adequately resolve the smallest scales, i.e., ∆x/η ext should be held constant across various simulations. In fact, this suggestion was also put forward by [20], though their criterion was stricter than the present numerical results suggest.
The simulations presented here suggest, based on PDFs of various components of the velocity gradient tensor, that a resolution of ∆x/η ≈ 1 is sufficient for R λ = 240, but barely insufficient for R λ = 390 (this is also evident from Fig. 6). An earlier resolution study at R λ ≤ 240 [16], also supports this. In fact, in [16], the authors also explored a practical approach to determine the necessary resolution based on calculating the error between p-th order structure function and its analytic behavior for small distances (obtained from Taylor series expansion), such that the result is a function of p and R λ . However, such expressions are limited to small values of p, since the derivation retained only a small number of terms in the Taylor expansion and hence also cannot be generalized to estimate η ext .
Nevertheless, based on previous and current results, it follows empirically that ∆x/η to accurately resolve the velocity gradients should be where R λ * ≈ 300 is the reference Taylor-scale Reynolds number, at which η ext ≈ η. The above relation provides a practical resolution criteria for future simulations at even larger problem sizes than considered here. While α ≈ 0.275 for the current range of R λ , we anticipate newer simulations at higher R λ would progressively update α and also quantify its dependence on R λ (though given the slow growth of α with R λ , a very substantial range of R λ might be required). The constraint provided by Eq. (16) should also apply to experimental investigations, which are currently capable of providing data at much higher R λ compared to DNS [40,48]. A simple estimate suggests that η/η ext 3 corresponding to these laboratory experiments at R λ ≈ 6000 − 10000. While this correction is unlikely to affect the dynamics in wind tunnel experiments [40], it might enhance the quantum effects in liquid-He experiments at [48]. However, more quantitative studies, even at relatively lower R λ , would be useful, since currently resolving even η in such high R λ experiments is an outstanding technical challenge.

VII. CONCLUSIONS
Using very well-resolved DNS of isotropic turbulence, both in space and time, at Taylor-scale Reynolds number R λ ranging from 140 to 650, we have characterized the extreme fluctuations of the velocity gradients. In particular, we focused on the square of vorticity, Ω = ω i ω i , and strain, Σ = 2s ij s ij (synonymous with enstrophy and dissipation), which have the same mean value (equal to 1/τ 2 K ). Whereas the PDFs of Ωτ 2 K and Στ 2 K superpose well around their mean values, the extents of their tails strongly grows as R λ increases. We find that these tails can be empirically collapsed by using a smaller time scale τ ext , defined as τ ext = τ K × R λ −β , implying that the extreme velocity gradients in the flow grow as τ −1 ext . The numerical results indicate that β ≈ 0.775 ± 0.0025.
The above result is further validated by analyzing the PDFs of velocity increments δu r at distances r equal to or less than the Kolmogorov length scale η. Our results show that δu r can be as large as the velocity r.m.s. u ′ , and slowly increases with R λ . The excellent superposition of the rescaled PDFs of R λ −α δu r /u ′ , with α = β− 1 /2 over the range of R λ covered in this study, suggests that the largest velocity gradients consist of velocity increments of approximately u ′ over a size η ext ∼ η R λ −α . The existence of scales smaller than η is consistent with previous ideas, although our results do not quantitatively support the existing phenomenological theories. In particular, the assumption that extreme events correspond to a local Reynolds number of unity does not appear to be justified in vortex tubes, where extreme gradients are found to reside -as revealed by flow visualizations in Figs. 1 and 7 and also consistent with previous studies [24,27].
Further analysis of flow structures around intense gradients in Fig. 7 reveals that vorticity and strain are generally not spatially colocated, as suggested in some earlier studies [16,28]. Conditional averaging shows that strain acting on intense vorticity, is on average weaker than the vorticity, and shows an approximate power law behavior given by Eq. (6). It is important to note that the structures shown in Fig. 7 taken in isolation do not lead to a strong vorticity amplification. This suggests that the stretching, necessary to create the very intense velocity gradients, in a representation of the velocity field in terms of a Biot-Savart equation [49], could originate from a non-local mechanism. The observation that the strain acting on intense vortices is significantly weaker than corresponding vorticity (as reflected in exponent γ < 1 in Eq. (6)) can be viewed as a consequence of this non-locality. Using scaling analysis, we are able to quantitatively relate β with the exponent γ. The weak increase in γ with R λ suggests the same for β. This leaves open the possibility that β could asymptote to 1 (and α to 0.5), in the limit of R λ → ∞ -a simple extrapolation of our data suggests that β ≥ 0.99 would require R λ 20, 000. However, such high Reynolds numbers might not be feasible experimentally or numerically. To conclude, explaining the scalings discussed here, especially in the light of γ, remains an outstanding theoretical challenge. Much remains to be learned by analyzing wellresolved data from even higher Reynolds numbers than considered here, from both DNS and experiments.