Self-attenuation of extreme events in Navier-Stokes turbulence

Turbulent fluid flows are ubiquitous in nature and technology, and are mathematically described by the incompressible Navier-Stokes equations (INSE). A hallmark of turbulence is spontaneous generation of intense whirls, resulting from amplification of the fluid rotation-rate (vorticity) by its deformation-rate (strain). This interaction, encoded in the non-linearity of INSE, is non-local, i.e., depends on the entire state of the flow, constituting a serious hindrance in turbulence theory and in establishing regularity of INSE. Here, we unveil a novel aspect of this interaction, by separating strain into local and non-local contributions utilizing the Biot-Savart integral of vorticity in a sphere of radius R. Analyzing highly-resolved numerical turbulent solutions to INSE, we find that when vorticity becomes very large, the local strain over small R surprisingly counteracts further amplification. This uncovered self-attenuation mechanism is further shown to be connected to local Beltramization of the flow, and could provide a direction in establishing the regularity of INSE.

A parcel of fluid moving at velocity u(x, t) in a flow, where x ∈ R 3 is the spatial location and t is time, simultaneously undergoes rotation and shape deformation, respectively characterized by the vorticity vector ω = ∇×u and the strain rate tensor S ij = (∂ j u i + ∂ i u j )/2. Its evolution in time can thereby be described by the incompressible Navier-Stokes equations (INSE) written as the vorticity equation [1]: where D t = ∂ t + u j ∂ j is the material derivative and ν is the kinematic viscosity of the fluid. This equation simply expresses that along a parcel trajectory vorticity is nonlinearly stretched by the strain rate, and also subjected to viscous damping. An essential aspect of this stretching term is that it leads to amplification of vorticity, i.e. generation of enstrophy Ω = ω i ω i , via the production term P Ω = ω i ω j S ij , as readily seen by taking the dot-product of Eq. (1) with ω i [2]. The rate at which enstrophy is amplified, and whether it can overcome viscous damping to blow-up in finite time, remains one of the outstanding unsolved Clay Millennium Prize problems [3,4]. It is known that for a finite-time blow-up, P Ω must grow unbounded. In addition, it has also been proven that this unbounded growth can possibly only occur when the viscosity ν is sufficiently small [4], which would correspond to turbulent solutions of the INSE. In fact, it is well known that Ω is highly intermittent in turbulent flows, attaining values hundreds or thousands times its mean, becoming even more extreme as the relative strength of viscosity is decreased [5][6][7][8][9][10]. However, these extreme events are typically found to be arranged in tube-like structures [5][6][7][8][9][10], with geometrical properties * dhawal.buaria@ds.mpg.de deterring maximum possible amplification [2,11,12]. Nevertheless, the question remains open whether the non-linear amplification could overcome viscous damping when the flow is sufficiently turbulent.
A fundamental difficulty in analyzing Eq. (1) arises from the non-local coupling between vorticity and strain rate; which implies that strain acting on vorticity at a point, as in Eq. (1), is in fact coupled to the entire state of the flow. Specifically, this non-locality can be quantified by expressing the strain tensor as a Biot-Savart integral of the vorticity field over the entire 3D spatial domain: where r = x − x ′ (with r = |r|) and ǫ ijk is the alternating Levi-Civita symbol. Thus, the amplification of vorticity can be entirely written in terms of vorticity itself, but the above integral poses a serious mathematical challenge in understanding the mechanisms encoded in the nonlinearity. In the current work, by evaluating the above integral numerically, we provide evidence that as vorticity is amplified to large values, the strain induced locally will ultimately act to attenuate its further amplification. In order to extract the local strain induced from vorticity amplification, we consider the following decomposition, by splitting the integration domain into a spherical neighborhood of radius R and the remaining domain [13,14]: where [· · ·] denotes the integrand in Eq.(2). The first term S N L ij is the non-local or background strain acting on the vorticity to stretch it, whereas S L ij is the local strain, induced by the vorticity in its neighborhood in response to the stretching. Thereafter, the production term can also be decomposed as as P Ω = P L Ω + P N L Ω , where P L,N L Ω = ω i ω j S L,N L ij . For such a decomposition, explicit bounds on P N L Ω can be established in terms of the total kinetic energy of the flow [15]. Thus, an unbounded growth of P Ω is only possible through P L Ω . However, our results will demonstrate, that when R is small enough, the term P L Ω remarkably acts to attenuate extreme vorticity fluctuations. Further analysis reveals that this attenuation is also connected to local Beltramisation of the flow, i.e., preferrential alignment of vorticity with velocity, which is expected to deplete the growth of nonlinearity [16].
To analyze the complex interaction between strain and vorticity, we utilize our unique database generated through direct numerical simulations (DNS) of the INSE. The simulations correspond to canonical setup of forced homogeneous and isotropic turbulence in a periodic domain [9], and are performed using the well-known Fourier pseudo-spectral methods, thus allowing us to obtain any quantity of interest with highest accuracy practicable [17]. It is instructive to note that that the mathematical results typically obtained in R 3 can be readily generalized to our simulation in the T 3 torus. Using the largest grid sizes currently feasible in turbulence simulations, of up to 12288 3 points [18,19], the Taylor-scale Reynolds number R λ , which quantifies the turbulence intensity, is varied from 140 to 1300 in our simulations (corresponding to fully developed turbulence). Special attention is given to faithfully resolve the small-scales and hence the extreme events [10], keeping the grid spacing smaller than the Kolmogorov length scale, η = (ν 3 / ǫ ) 1/4 , based on the mean dissipation rate of kinetic energy ǫ , where the average · is taken over the 3D spatial domain and also multiple realizations. Note that the mean enstrophy Ω , is equal to ǫ /ν, due to underlying homogeneity [1]. Additional details about our DNS and database are provided in the Methods section.

RESULTS
Efficient determination of the local and nonlocal strain: While the vorticity and strain fields can be easily obtained from DNS, we have devised an efficient method to compute the local and non-local strain fields, without directly evaluating the prohibitively expensive Biot-Savart integral over the entire domain. As shown in [13], using a Taylor-series expansion of vorticity over a distance R, the non-local strain S N L (x, R) can be expressed in terms of the total strain as follows: Starting from the above expression and transforming it to Fourier space (where the differential operator ∇ 2 reduces to a simple multiplication by −k 2 ), leads to the relation where(·) denotes the Fourier transform, k is the wavenumber vector with k = |k| and f (kR) is an infinite series. In practice, truncating f (kR) to a finite number of terms can at best provide approximate results [13]. However, as derived in the Supplementary, one can show that f (kR) converges to the following expression: This allows us to evaluate the Biot-Savart integral in Eq.2 by applying a simple transfer function to the total strain rate in Fourier space, and thus to obtain S L,N L ij (and P L,N L Ω ) very accurately for any value of R. Interestingly, it is worth noting that f (kR) in Eq. (6) corresponds to the sinc function in 3D, which also happens to be the Fourier transform of a box or top-hat filter (of radius R), commonly utilized in other disciplines, e.g. large-eddy simulation (LES), signal processing. Thus, evaluating the non-local strain essentially reduces to a filtering operation on the total strain.
Visualization of extreme events: Figure 1 illustrates our main result, namely that the local contribution to stretching, P L Ω , is in fact negative in the neighborhood of extreme vorticity events. The visualizations shown in Fig. 1 focus on a small domain of size (50η) 3 around one of the extreme vorticity events in the flow (with the most intense vorticity at the center). Figure 1a and b show isosurfaces of enstrophy, respectively at 100 and 1000 times the mean value corresponding to moderate and intense events, and illustrate the characteristic vortex-tube structure [8][9][10]. The cut through the midplane of the domain is shown in Fig. 1c, and demonstrates the sharp variation of enstrophy across the cross section of the tubes. Figure 1d-f show the total production P Ω for the same field. In Fig. 1d, isosurfaces are shown for levels ±400 (with cyan and red corresponding to positive and negative values respectively), which approximately correspond to moderate enstrophy (as shown in Fig. 1a). Whereas in Fig. 1e, isosurfaces are shown for ±1000, which correspond to intense enstrophy (as shown in Fig. 1b). In Fig. 1f, the 2D contour field at the mid-plane is shown. The main observation is that P Ω is overwhelmingly positive, which is anticipated given large enstrophy in these tubes, and also from dynamical constraints of turbulence [20].
Finally, Fig. 1g-i shows the contribution P L Ω from local strain for R = 2η. In Fig. 1g and h, isosurfaces are shown for levels ±50 and ±200 respectively, once again corresponding to moderate and intense enstrophy events respectively. Unlike P Ω which is always positive on average Middle row: enstrophy production based on total strain, suitably non-dimensionalized by mean of enstrophy, at thresholds of (d) ±400, and (e) ±1000, which approximately correspond to moderate and intense enstrophy, shown in (a) and (b) respectively. (f) 2D contours at the mid-plane. The production terms based on total strain is overwhelmingly positive. Bottom row: enstrophy production based on local strain, once again suitably non-dimensionalized by mean enstrophy, at thresholds of (h) ±50, and (g) ±200, again corresponding to moderate and intense enstrophy shown in (a) and (b) respectively. (i) 2D contours at the mid-plane, revealing that the production term based on local strain is strongly negative in the regions of intense vorticity. For each row, the thresholds shown in first two isosurfaces plots are marked by dashed and solid lines respectively in last 2D-contour field plot.
[20], the mean of P L Ω has no such constraints. For moderate values shown in Fig. 1g, we find that the volumes occupied by positive and negative values are comparable. However, for intense value shown in Fig. 1h, negative stretching rate is more prevalent especially around the center where vorticity is maximum. This is corroborated by Fig. 1i, which shows the 2D contour level of P L Ω at the mid plane and reveals that both negative and positive values occur in the outer regions of the tubes where vorticity is not very intense; whereas large negative values occur inside the tubes, where vorticity is most intense.
Let us briefly mention that the flow structure presented in Fig. 1 represents one generic scenario of how the regions of intense vorticity look like. Needless to say, we inspected many such regions, and note that all of them qualitatively behave in the same manner, and essentially lead to the same conclusion. We have included another such example in the Supplementary.
Conditional statistics: To establish the quantitative significance of the observations in Fig. 1, Fig. 2a shows the average of P L Ω /Ω conditioned on Ω, for R = η and 2η, and various Reynolds numbers. Note P L Ω /Ω = ω iωj S L ij (where(·) is the corresponding unit vector) and provides the measure of effective strain engendering enstrophy production, irrespective of the strength of vorticity [2,11]. The Taylor expansion in Eq. (4) implies that for small R, S L ij can be written as which suggests that the local strain is in fact proportional to the Laplacian of the total strain. Hence for comparison, we have also shown the conditional expectation ω iωj η 2 ∇ 2 S ij |Ω in Fig. 2a, and P L Ω is accordingly multiplied by 10η 2 /R 2 . The conditional production term is virtually zero for small to moderate values of Ω -consistent with strong cancellation between negative and positive values seen in Fig. 1g. However, as Ω gets larger, the expectation P L Ω |Ω becomes negative for all Reynolds numbers and strongly increases in magnitude with Ω. We note that the values of P L Ω are overwhelmingly negative for large Ω, as corroborated by the observation (not shown in figure) that conditional expectations of |P L Ω | and | − P L Ω | are virtually equal. In addition, in Fig. 2b, we have show the conditional root-mean-square (rms) of the fluctuations of the P L Ω , normalized in the same manner as Fig. 2a. Once again, we have included the corresponding curve for η 2 ∇ 2 S ij for comparison. Remarkably, we observed the exact behavior as seen in panel Fig. 2a (except the curves are all on the positive side, because the rms is always positive by definition). At the same time, we note that the curves in both Fig. 2a and b, have comparable values, i.e., the mean and rms are comparable (especially for large Ω). This reaffirms that P L Ω is predominantly negative when conditioned on large values of Ω, and thus consolidates the observed self-attenuation mechanism. Finally, it is Negative contribution of local strain to production of enstrophy. (a) Averaged enstrophy production due to the local strain, P L Ω , conditioned on enstrophy normalized by its mean value. The curves shown correspond for R/η = 1 and 2, at Taylor-scale Reynolds numbers R λ = 390 − 1300. For comparison, we also show the contribution based on η 2 ∇ 2 Sij , which is the limiting value of local strain for small R as noted in Eq. (7). (Accordingly the curves for R/η = 1 and 2 are also adjusted by a factor of 10η 2 /R 2 ). (b) The conditional root-mean-square σ X|Ω of the local enstrophy production term (X =ωiωjS L ij ), defined as σ 2 X|Ω = X 2 |Ω − X|Ω 2 . Similar normalization as panel (a) is used.
worth noting that as R/η becomes smaller the curves for a given Reynolds number expectedly approach the analytical limit given by Eq. 7. The result for R/η = 0.5 (not shown), was found to be virtually indistinguishable from the corresponding curve showing the analytical limit.
The observation from Figs.1 and 2 that extreme vorticity fluctuations are accompanied by negative values of P L Ω indicates that the strain induced locally acts to prevent further growth of enstrophy. It is important to realize that this mechanism is separate from viscous diffusion or dissipation of enstrophy [21], but still acts in conjunction with it. Additionally, this self-attenuating mechanism is far stronger than a mere reduction (depletion) of nonlinearity [16,22]. Depletion of non-linearity essentially refers to weakening of vortex stretching (compared to its maximum possible amplitude) [23], which is evidently reflected in alignment of vorticity with intermediate eigenvector of strain tensor and hence the weak curvature of vortex tubes [11,24] -as also seen in Fig. 1a-b. However, the presence of self-attenuation suggests that the non-linearity itself could be capable of preventing a runaway blowup, even as viscosity gets small (as suggested by Fig. 2a, where the increase of P L Ω is merely shifted to larger values of Ω as viscosity decreases). A careful mathematical analysis of this mechanism and determining mathematical bounds on P L Ω could possibly reveal a path in establishing global regularity of INSE.
Connection to helicity: The presence of negative local stretching accompanying intense vorticity raises additional questions about the local flow structure. Given that intense vorticity is arranged in tubes with weak curvature, additional insight could be obtained by a simple kinematic analysis of stretching generated by such structures. To this end, we consider a simple axisymmetric vortex tube with a radius of curvature R c [25][26][27]. Utilizing a curvilinear polar coordinate system: (r,θ,ŝ), which respectively correspond to unit vectors in the radial direction, the azimuthal direction and the direction tangent along the (curved) axis of the tube, we assume that the vorticity is of the form ω = ω s (r, s)ŝ+ω θ (r, s)θ. The component ω s corresponds to azimuthal velocity in the tube similar to a two-dimensional Burgers vortex [28], whereas the component ω θ comes from axial velocity along the tube. Thereafter, utilizing Eq. (7), one can derive (as shown in the Supplementary): which gives the local stretching induced by the vortex tube as sum of two terms, involving F and G, which are functions of ω s and ω θ and their derivatives.
The term with a cos θ/R c dependence results from the curvature of the tube and produces a dipolar structure, with positive and negative contributions depending on the sign of cos θ [29] -consistent with the structure seen in Fig. 1g and i. In contrast, the term independent of cos θ acts as a monopole. Based on the results shown in Figs. 1 and 2, the sign of F must be positive, and would result in attenuation of intense vorticity by the local strain. Interestingly, F is identically zero if the component ω θ vanishes, i.e. there is no axial flow velocity. This suggests that some local alignment between vorticity and velocity must occur when vorticity is large. Interestingly, a similar conclusion can also be reached by realizing that the non-linear terms in INSE, in Eq. (1), can be rewritten as ∇ × (u × ω). Thus, local Beltramization, i.e., alignment of u and ω in regions of large enstrophy would essentially act to restrict the non-linear amplification [16].
The above prediction is consistent with earlier results at low R λ [30], as well as with our own results at significantly higher R λ in Fig. 3, which shows the conditional average of the cosine between velocity and vorticity, conditioned on enstrophy. The average is taken over the absolute value, since the sign of the cosine is immaterial to measure the degree of Beltramization (note that the dot product of velocity and vorticity is not sign-definite). For small values of Ω, the average stays constant at 0.5, consistent with a uniform distribution of the cosine. However, the conditional average increases at large Ω, in good correlation with the increase of the magnitude of P L Ω seen in Fig. 2a. Thus, in fully developed turbulence, the intense whirling motions (vortex tubes), emblematic of the small-scale structures, are innately three-dimensional and helical.

DISCUSSION
In conclusion, we have utilized very well resolved numerical simulations of fully developed turbulence to investigate extreme fluctuations of vorticity, which can be considered as signatures of potential singularities of INSE. Our results show that when vorticity is strongly amplified, the non-linearity in its local neighborhood remarkably counteracts further amplification, instead of enhancing it. In addition, this effect gets stronger as vorticity gets stronger and also as Reynolds number increases (or viscosity decreases). Thus, our results suggest that the non-linearity -which is responsible for amplification in the first place -also encodes a mechanism (in conjunction with viscosity), which can prevent a finitetime singularity from occurring. A deeper understanding of this self-attenuation mechanism based on a clear physical argument could help to set stronger mathematical bounds on the stretching of vorticity [15,31], and could be an essential ingredient to prove global regularity of the INSE [3].
Another important observation in this regard is the local Beltramization of the flow in regions of large enstrophy -highlighting the helical nature of the small-scales of turbulence (which are structurally arranged in vortex tubes). While it was anticipated that reduction (depletion) of non-linearity would lead to such helicity [16], the uncovered self-attenuation mechanism shows that the effect is in fact much stronger, and directly counteracts vorticity amplification. A promising direction in this regard could be to extend the ideas based on helical decomposition to further analyze this local Beltramization [32,33]. Indeed, a recent work has established global regularity for a decimated version of INSE which enforces helicity to be sign-definite [34]. A possible extension to full INSE, in light of the uncovered self-attenuation mechanism, presents an important challenge for future work.
Finally, it is worth noting that our numerical simulations of stationary isotropic turbulence do not address specific initial value problems, such as those involving collisions between two or more vortex tubes [29,[35][36][37]. Such special flow configurations are routinely studied to investigate the development of a possible finite-time singularity, mostly in the context of inviscid flows (ν = 0), i.e., the Euler equations. However, a conclusive demonstration of a blowup or lack thereof still remains elusive [38]. While complicated interactions between vortex tubes already occur in our simulations, it remains to be understood how the ideas developed here would apply to these special configurations.

Direct numerical simulations:
The data utilized in the current work are generated through direct numerical simulations (DNS) of the incompressible Navier-Stokes equations (INSE) where u is the divergence free velocity field (∇ · u = 0), P is the pressure, ρ is the fluid density, ν is the kinematic viscosity, and f corresponds to large scale forcing used to maintain a statistically stationary state [39]. The equations are solved using a massively parallelized version of the well-known Fourier pseudo-spectral algorithm of Rogallo (1981) [40]. The aliasing errors resulting from the convolution sums are controlled by grid shifting and spherical truncation [41]. Our DNS corresponds to the canonical setup of homogeneous and isotropic turbulence with periodic boundary conditions on a cubic domain of side length L 0 = 2π, which is ideal for studying small scales and hence extreme events at highest Reynolds numbers possible [9]. The domain is discretized using N 3 grid points, with uniform grid spacing ∆x = L 0 /N in each direction. We utilize explicit second-order Runge-Kutta for time integration, where the time step ∆t is subject to the Courant number (C) constraint for numerical stability: ∆t = C∆x/||u|| ∞ (where || · || ∞ is the L ∞ norm). The DNS database used in the current work is summarized in Table I, along with the main simulation parameters. An important consideration in studying extreme events is that of spatial resolution, which is measured in pseudo-spectral DNS by the parameter k max η, where k max = √ 2N/3 is the maximum resolved wavenumber on a N 3 grid and η is the Kolmogorov length scale. Equivalently, one can use the ratio ∆x/η which is approximately equal to 3/k max η. The runs with Taylor-scale Reynolds numbers, R λ , in the range 140 ≤ R λ ≤ 650 were also utilized in our recent work [10] and all have a very high spatial resolution, k max η ≈ 6 (or ∆x/η ≈ 0.5). This resolution should be compared to the one used in comparable numerical investigations of turbulence at high Reynolds numbers, which are mostly in the range 1 ≤ k max η ≤ 1.5 [9,42] -which do not resolve the extreme events adequately. In addition to our previous runs, we have performed a new run at significantly higher R λ of 1300, on a larger 12288 3 grid with a small-scale resolution of k max η = 3 (or ∆x/η ≈ 1). This is one of the largest DNS reported to date -comparable with [18] which also reported results from 12288 3 run at R λ = 2300, but with k max η ≈ 1 (where the small-scales were not properly resolved).
We have also listed the simulation length T sim used for generating independent ensembles, in terms of the large-eddy turnover time (T E ) or the Kolmogorov time scale (τ K ). The statistical results are obtained by averaging over N s independent ensembles, which are uniformly spread out over the simulation length. Note, the range of time scales is typically given by the ratio T E /τ K , which scales linearly with R λ [1]. However, the time scale of extreme events which we consider here is smaller than τ K , getting even smaller as R λ increases [10].

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author on request.

CODE AVAILABILITY
The simulation and post-processing codes that have been used to produce the results of this study are available from the corresponding author on request.