Vortex stretching and enstrophy production in high Reynolds number turbulence

An essential ingredient of turbulent flows is the vortex stretching mechanism, which emanates from the non-linear interaction of vorticity and strain-rate tensor and leads to formation of extreme events. We analyze the statistical correlations between vorticity and strain rate by using a massive database generated from very well resolved direct numerical simulations of forced isotropic turbulence in periodic domains. The grid resolution is up to $12288^3$, and the Taylor-scale Reynolds number is in the range $140-1300$. In order to understand the formation and structure of extreme vorticity fluctuations, we obtain statistics conditioned on enstrophy (vorticity-squared). The magnitude of strain, as well as its eigenvalues, is approximately constant when conditioned on weak enstrophy; whereas they grow approximately as power laws for strong enstrophy, which become steeper with increasing $R_\lambda$. We find that the well-known preferential alignment between vorticity and the intermediate eigenvector of strain tensor is even stronger for large enstrophy, whereas vorticity shows a tendency to be weakly orthogonal to the most extensive eigenvector (for large enstrophy). Yet the dominant contribution to the production of large enstrophy events arises from the most extensive eigendirection, the more so as $R_\lambda$ increases. Nevertheless, the stretching in intense vorticity regions is significantly depleted, consistent with the kinematic properties of weakly-curved tubes in which they are organized. Further analysis reveals that intense enstrophy is primarily depleted via viscous diffusion, though viscous dissipation is also significant. Implications for modeling are nominally addressed as appropriate.

An essential ingredient of turbulent flows is the vortex stretching mechanism, which emanates from the non-linear interaction of vorticity and strain-rate tensor and leads to formation of extreme events. We analyze the statistical correlations between vorticity and strain rate by using a massive database generated from very well resolved direct numerical simulations of forced isotropic turbulence in periodic domains. The grid resolution is up to 12288 3 , and the Taylor-scale Reynolds number is in the range 140 − 1300. In order to understand the formation and structure of extreme vorticity fluctuations, we obtain statistics conditioned on enstrophy (vorticity-squared). The magnitude of strain, as well as its eigenvalues, is approximately constant when conditioned on weak enstrophy; whereas they grow approximately as power laws for strong enstrophy, which become steeper with increasing R λ . We find that the well-known preferential alignment between vorticity and the intermediate eigenvector of strain tensor is even stronger for large enstrophy, whereas vorticity shows a tendency to be weakly orthogonal to the most extensive eigenvector (for large enstrophy). Yet the dominant contribution to the production of large enstrophy events arises from the most extensive eigendirection, the more so as R λ increases. Nevertheless, the stretching in intense vorticity regions is significantly depleted, consistent with the kinematic properties of weakly-curved tubes in which they are organized. Further analysis reveals that intense enstrophy is primarily depleted via viscous diffusion, though viscous dissipation is also significant. Implications for modeling are nominally addressed as appropriate.

I. INTRODUCTION
Small-scale intermittency, a hallmark of fluid turbulence, refers to the occurrence of sudden and intense fluctuations of velocity gradients [1,2], as routinely reflected in long tails of their strongly non-Gaussian probability distributions [3][4][5][6]. Understanding such extreme events is of paramount practical interest. For example, strong strain rates can dramatically enhance mixing of scalars or adversely affect flame propagation in reacting flows [7,8], whereas strong vortical motions engender clustering of particles, facilitating cloud and rain formation [9,10]. From a fundamental standpoint, amplification of gradients is an essential component of the energy cascading process, leading to generation of small scales in turbulent flows [11][12][13]. Thus, characterizing intermittency and the associated generative mechanisms is at the heart of turbulence theory and modeling [14,15].
A key mechanism responsible for the formation of such extreme events and the small-scales is the process of vortex stretching [11], which results from non-linear coupling between vorticity and strain-rate tensor, respectively defined as, ω = ∇ × u, and s ij = 1 2 (∂u i /∂x j + ∂u j /∂x i ), where u  1. (a) Probability density functions (PDFs) of alignments between the vorticity unit vectorω, and the eigenvectors e i , corresponding to the eigenvalues λ i of the strain-rate tensor (with λ 1 ≥ λ 2 ≥ λ 3 ). (b) PDF of β, defined in Eq. (6), which measures the relative strength of the intermediate eigenvalue with respect to the overall strain amplitude. Solid lines for Taylor-scale Reynolds number R λ = 1300 and dashed lines for R λ = 140. There seems to be virtually no effect of increasing R λ on the PDFs. is the velocity field. This is evident from the vorticity transport equation: where D/Dt = ∂ t + u · ∇ is the material derivative, ν is kinematic viscosity and the vector ω j s ij gives the aforementioned vortex stretching term. In all turbulent flows, the vortex stretching term directly results in vorticity amplification and hence net production of enstrophy (vorticitysquared) [16]. Numerous studies over last few decades, both from experiments and direct numerical simulations (DNS), have revealed some robust features of vortex stretching, which seemingly appear to be universal across different turbulent flows. Notably, the vorticity vector preferentially aligns with the eigenvector corresponding to the intermediate eigenvalue of the strain tensor, which in turn is positive (on average). [14][15][16][17][18]. Fig. 1 summarizes these findings utilizing data from DNS of isotropic turbulence (the details of which are discussed in Section II). Quantitatively the results can be expected to be slightly different for other turbulent flows, but qualitatively they behave similarly. Remarkably, the results also appear virtually independent of the Reynolds number. The persistent correlations between strain and vorticity demonstrated in Fig. 1, highlight some important aspects of vortex stretching in turbulence. Particularly, the alignment between vorticity and the intermediate strain eigenvector suggests a depletion of stretching, compared to what can be expected if vorticity were to align with the most extensive strain eigenvector. However, the statistics shown in Fig. 1 are obtained from a uniform sampling of the flow, and do not distinguish quiescent regions from those where extreme events reside. It is well known that the extreme events in the flow are organized in localized vortex tubes, where enstrophy is orders of magnitude larger than its global average [4,6,19], whereas quiescent regions are often believed to be structureless. Thus, the mechanisms of vortex stretching may differ significantly between such regions. In this regard, statistics conditioned on the strength of vorticity or strain offer the unique prospect of understanding how small-scale structures are produced by the flow, by specifically isolating the extreme events from the moderate and the weak events. This can be especially useful in understanding the dynamics of vorticity amplification [14]. In addition to providing fundamental insight, such conditional statistics are also very useful in turbulence modeling, especially in PDF methods or similar approaches [20][21][22][23].
In this paper, our objective is to present a detailed fundamental investigation of turbulence small-scale structure in light of vortex stretching and resulting enstrophy production. To this end, we utilize pseudo-spectral DNS of isotropic turbulence in a periodic domain, which is the most efficient numerical tool to study the small-scale properties of turbulence. One important purpose of the current study is also to understand the effect of increasing Reynolds number. With that in mind, we utilize a massive DNS database with Taylor-scale Reynolds number R λ ranging from 140 to 1300 on grid sizes 1024 3 to 12288 3 respectively -with particular attention on having very good small-scale resolution to resolve the extreme events accurately [24]. We compute statistics related to vortex stretching, conditioned on magnitude of vorticity. We find that the alignment of vorticity with intermediate strain eigenvector is strongly enhanced as vorticity becomes more intense. This is qualitatively consistent with the notion that the most intense velocity gradient structure have a quasi two-dimensional (2D) structure, for which the role of the largest strain eigenvalue in vortex stretching is inhibited [25,26]. However, contrary to this, we find that the dominant contribution to enstrophy production still comes from the most extensive eigendirection; though the contribution from intermediate direction is also comparable. Nevertheless, the stretching in the most intense vorticity regions is significantly depleted. Further analysis shows that viscous diffusion of enstrophy is primarily responsible for arresting extreme events (and dominates over viscous dissipation of enstrophy). We also propose improved functional forms for modeling the statistics reported in this work.
The rest of the manuscript is organized as follows. In Section II, we present the database used in this work. Correlations between strain and vorticity, in particular their dependence on the Reynolds number, are discussed in Section III. Using conditional averages, the dependence of these correlations on the vorticity are presented in Sections IV Finally, we summarize our results in Section V. Some details on numerical resolution, crucial in the study of small scale turbulent statistics, are presented in the Appendix A.

II. NUMERICAL APPROACH AND DATABASE
The data used in the present work was generated via DNS of the incompressible Navier-Stokes equations where u is the divergence free velocity field (∇ · u = 0), p is the pressure, ρ is the fluid density, ν is the kinematic viscosity and f is the forcing term at large scales to achieve a statistically stationary state [27]. We utilize a massively parallel implementation of Rogallo's pseudo-spectral algorithm [28] in a periodic (2π) 3 cubic domain, with explicit second-order Runge-Kutta scheme for time integration. The aliasing errors are controlled by a combination of truncation and phaseshifting [29], resulting in highest resolved wavenumber given as k max = √ 2N/3, where N is the number of grid points in each direction.
The DNS database used in the current work is summarized in Table I, along with the main simulation parameters. The runs with Taylor-scale Reynolds numbers, R λ , in the range 140 ≤ R λ ≤ 650 were also utilized in our recent work [6] and all have a very high spatial resolution, k max η ≈ 6, or alternatively ∆x/η ≈ 0.5, where ∆x = 2π/N is the grid spacing and η is the Kolmogorov length scale. 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 η ≤ 2 (or equivalently, 1.5 ≤ ∆x/η ≤ 3) [5,30]. As established in our recent work [6], such a resolution is more than adequate to accurately resolve and study the small-scales at the R λ R λ N 3 k max η T E /τ K T sim N s 140 1024 3  ), the number of grid points (N 3 ), spatial resolution (k max η), ratio of large-eddy turnover time (T E ) to Kolmogorov time scale (τ K ), length of simulation (T sim ) in statistically stationary state and the number of instantaneous snapshots (N s ) used for each run to obtain the statistics. considered, although the resolution constraints appear to be more severe as the Reynolds number increases.
In addition to our previous runs, we consider a new run at R λ = 1300 with a resolution of k max η ≈ 3 (or ∆x/η ≈ 1) [31]. While this resolution may not be sufficient at R λ = 1300 to accurately determine some of the quantities discussed in [6], such as the probability distributions of strain or vorticity (a resolution of k max η ≈ 4.5 would be more appropriate, according to the arguments in [6]), we note that our focus in the current work are the conditional statistics, which are known to be less sensitive to resolution [32]. (We also provide additional tests in the Appendix A). Additionally, the length of this new run is comparatively short. However, the characteristic time scale of the most intense events becomes much smaller than the Kolmogorov time scale τ K as R λ increases [6], allowing for adequate statistical convergence by simply sampling more frequently in time [24].

III. VORTICITY AND STRAIN CORRELATIONS: UNCONDITIONAL STATISTICS
We begin by analyzing unconditional results on vortex stretching, which will be useful for a better understanding of the conditional statistics described in the next section. The results presented here are generally consistent with previous numerical studies at lower Reynolds numbers. By improving the numerical resolution, and by extending the range of values of R λ considered, we provide additional insight, especially into the Reynolds number dependence.
In order to quantify the intensity of vorticity we simply utilize the enstrophy Ω = ω i ω i . From the vorticity transport equation in Eq. (1), the following transport equation for enstrophy can be derived The viscous terms have been decomposed into a diffusion contribution ν∇ 2 Ω, and a purely dissipative contribution, ǫ ω = 2ν||∇ω|| 2 . In statistically stationary isotropic turbulence, as considered here, applying averaging to above equation gives which gives the balance between production and dissipation of enstrophy. The above relation is in fact also approximately valid for other turbulent flows at sufficiently high Reynolds numbers [11]. Our analysis is primarily focused around this enstrophy production term, which is evidently affected by both the local alignment between strain and vorticity and their absolute magnitudes. This interaction between strain and vorticity is most readily described in the eigenframe of the strain tensor [17,33]. The strain tensor can be diagonalised into an orthonormal basis, with its eigenvalues given by λ i for i = 1, 2, 3 (with λ 1 ≥ λ 2 ≥ λ 3 ) and the corresponding eigenvectors denoted by e i . Incompressibility imposes λ 1 + λ 2 + λ 3 = 0, which in turn renders λ 1 to be always positive (stretching) and λ 3 to be always negative (compressive). Thus the enstrophy production term can be rewritten as which isolates individual contribution from each eigendirection.   Various individual contributions in each eigendirection are computed and listed in Table II and also plotted in Fig. 2. All (dimensional) quantities are non-dimensionalized appropriately by the Kolmogorov time scale τ K , defined as Ω = ǫ /ν = 1/τ 2 K , where ǫ = 2νs ij s ij is the energy dissipation rate. We first consider the mean of each eigenvalue as shown first in the table and plotted in Fig. 2a. Note that they sum up to zero due to incompressibility. Interestingly, it appears the contribution in each direction slowly decreases with increasing R λ , however, the ratio λ 1 / λ 2 is approximately constant (at about 4.3). Alternatively, the quantity β defined as which measures the relative amplitude of λ 2 , is also independent of R λ . This is indeed to be expected, based on Fig. 1, where we saw that the PDF of β itself does not change with R λ .
Although not shown, we also observe that the PDF of λ 2 /λ 1 is also independent of R λ (and qualitatively similar to that of β). Additionally, their PDFs suggest that the most probable value of λ 1 : λ 2 : λ 3 is approximately 2.85 : 1 : −3.85, which is in good agreement with earlier works [17,33].
On the other hand, the mean-square values of individual eigenvalues (the next quantity in Table II), once scaled with τ 2 K , are essentially independent of R λ . Note by definition, they sum up to 1/2, since λ i λ i = s ij s ij . Similarly, the mean-square of alignment cosine between vorticity and eigenvectors of strain -defined as (e i ·ω) 2 , whereω is the unit vector along vorticity -also do not reveal any dependence on R λ . This is again consistent with Fig. 1a, where the PDFs of these cosines were seen to be practically identical at the lowest and highest R λ in this work. Note that the square of cosines is considered because they have the nice property of summing up to unity for all three directions (in addition to being bounded between 0 and 1 for each individual direction for the case of perfect orthogonal and parallel alignment respectively).
Next in Table II, we consider the individual contribution to net enstrophy production through the terms λ i (e i ·ω) 2 . It follows that the contribution from the first direction is always positive and that from the third direction is always negative. We find that λ 1 (e 1 · ω) 2 is approximately twice as large as − λ 3 (e 3 · ω) 2 . Additionally, the intermediate direction is strongly positive resulting in net positive enstrophy production. Once non-dimensionalized, we find that the absolute values of the three components λ i (e i · ω) 2 τ 3 K all increase with R λ . Importantly, despite the strong alignment with intermediate eigenvector, the net contribution to enstrophy production from the first eigendirection is always larger than the second. In fact, the difference between the contributions of the first and second strain eigendirections increases with R λ .
To further clarify, the individual contribution in each eigendirection as a fraction of total production is shown in Fig. 2c (absolute value is shown for comparison). Interestingly, we find that the relative contribution from the intermediate direction is virtually constant with increasing R λ , whereas the contributions from the first and third directions both increase in the same manner (though they are opposite in signs). However, it appears that this increase is prominent only at smaller R λ and an asymptotic state is possibly reached at high R λ , such that the relative contribution from each eigendirection is constant, and approximately in the ratio 0.78 : 0.62 : −0.40.
Whereas our results summarized in Table II and Fig. 2 agree qualitatively with those of earlier studies, we notice a significant quantitative difference compared to the results of [14,33,34], especially concerning the individual eigenvalues and their contribution to net enstrophy production. On the other hand, results of [35] at R λ ∼ 100 based on simulations of grid turbulence, appear to be in reasonable agreement with our results at R λ = 140. This suggests that the disagreement with the results of [14,33,34] is likely due to practical difficulties in experimentally measuring the velocity gradient tensor with sufficient accuracy [14,36].
Finally, the increase in net enstrophy production with R λ , documented in Fig. 2b, can be readily explained utilizing the well-known relation: where S is the skewness of the longitudinal velocity gradients [16,37]. It is also well known that due a forward cascade of energy in (three-dimensional) turbulent flows, S is negative and its magnitude slowly increases with R λ , approximately as a power law [19]. The final column of Table II lists the values of −S for different R λ , which are in excellent agreement with Eq. (7). On a similar note, given the far improved small-scale resolution utilized in current work, our values of S suggest a minor correction to the power law fit reported in [19], with the new fit given as This approximate power law also applies for the net enstrophy production (and hence net enstrophy dissipation through Eq. (4)), albeit with a different prefactor.

IV. STATISTICS CONDITIONED ON VORTICITY/ENSTROPHY
The unconditional correlations of strain and vorticity reviewed in the previous section revealed some very robust features. In this section, to specifically understand the formation of extreme events, we condition the correlations on the magnitude of vorticity. More precisely, we use Ωτ 2 K as the conditioning variable, which is simply the enstrophy divided by its mean value, i.e., Ω/ Ω . Since the distribution of Ωτ 2 K varies over a very large range of values, the conditioning is performed using logarithmically spaced bins, with 8 bins per decade (same as in [6]). To facilitate comparisons between runs at different R λ , the investigated conditional averages are also non-dimensionalized by Kolmogorov scales (as necessary).

A. Alignment
We first investigate the conditional alignment properties between vorticity and the eigenvectors of strain tensor by considering the conditional PDFs of the directional cosines e i ·ω. Similar to Fig. 1a, the absolute magnitude of the cosine is used, since its sign is immaterial for enstrophy production. The conditional PDFs corresponding to the intermediate eigenvector, i.e., i = 2, are shown in Fig. 3a, for the highest R λ (= 1300) in this work and for the range 10 −2 ≤ Ωτ 2 K ≤ 10 3 ; for comparison the marginal PDF (from Fig. 1a) is shown as a black dashed line. For small values of Ωτ 2 K , the preferential alignment between e 2 andω is reduced, and PDF of e 2 ·ω approaches a uniform distribution. On the other hand, for large values of Ωτ 2 K , the alignment between e 2 and ω becomes much stronger, with the PDF of e 2 ·ω. becoming increasingly concentrated around 1. The alignment appears to be almost independent of Ωτ 2 K for Ωτ 2 K 100, suggesting an asymptotic state for the most extreme events.
The conditional PDFs for i = 1 and 3 are shown in Fig. 3b and c respectively. The marginal PDFs shown in black dashed lines (again from Fig. 1a), did not reveal any preferential alignment between vorticity and e 1 , with its distribution close to uniform; whereas vorticity was preferentially orthogonal to e 3 . The conditional PDFs ofω · e 1 on vorticity reveals a very weak preferential alignment between e 1 and vorticity at small values of Ωτ 2 K . However, this alignment disappears when Ωτ 2 K 1. For Ωτ 2 K 10,ω and e 1 appear to be preferentially orthogonal to each other (albeit very weakly). On the other hand, the conditional PDFs ofω · e 3 , suggest that vorticity is always preferentially orthogonal to e 3 , except at small conditioning values, where the PDF approaches an uniform distribution (similar to those in other directions). To summarize, as the value of Ωτ 2 K becomes larger, both e 1 and e 3 appear to become increasingly orthogonal to vorticity, with the effect being strongly pronounced for e 3 . Simultaneously, the alignment of vorticity with e 2 becomes extremely strong at large Ωτ 2 K . The above result can be concisely illustrated by considering the second moment of the conditional PDFs, i.e., (e i ·ω) 2 |Ω and plotting it as a function of Ωτ 2 K . Since the directional cosine is bounded between 0 and 1, the conditional expectation is also bounded between 0 and 1. (with a value of 1/3 for a uniformly distributed PDF). In addition, the second moment has the nice property that the contributions in three directions adds up to unity, i.e., 3 i=1 (e i ·ω) 2 |Ω = 1, for any value of Ω. Such a representation also allows us to systematically investigate the R λ dependence, as evident from Fig. 3d, which shows the conditional expectation (e i ·ω) 2 |Ω for i = 1, 2, 3 and various R λ . The behavior of the conditional expectations in Fig. 3d, is completely consistent with that of conditional PDFs in Fig. 3a-c, and clearly demonstrates the incipient plateau-like behavior for intense vorticity. Interestingly, from this figure, it is also evident that the curves show a very weak dependence on R λ .
Overall, the above results point to very different statistical properties of the strain field in regions of intense vorticity, compared to those in regions of moderate or weak enstrophy. Manifestly, the weakest events are 'structure-less', with all alignments being equally likely. On the other hand, intense enstrophy regions appear to have a very specific structure, with the conditional alignments (as seen from Fig. 3d) approximately in the ratio 0.2 : 0.7 : 0.1. Since the most intense vorticity events are typically found to be arranged in tubes with weak curvatures, the above result appears consistent with an effectively 2D structure -where the most extensive and compressive eigenvectors lie in the equatorial plane, and the intermediate eigenvector is along the tube axis. In fact, such an alignment effect was also established in [25], and shown to result from kinematic constraints. A similar behavior is also observed in the interaction between vortex structures in the Euler equations, where the formation of very intense structure, shaped as 2D vortex sheets is often observed, see e.g. [26]. For such structures, the largest and weakest strain eigenvectors are perpendicular to vorticity, and amplification is due to the weak intermediate eigenvector. This results in an effective reduction of the nonlinearity [14].

B. Strain and eigenvalues
In addition to the alignments of vorticity and strain eigenvectors, the magnitude of strain and its eigenvalues are of obvious importance in determining net enstrophy production. We characterize the magnitude of the strain tensor by Σ = 2s ij s ij , which is simply the dissipation rate divided by viscosity, i.e., Σ = ǫ/ν. In homogeneous turbulence, as considered here, we have Σ = Ω = 1/τ 2 K , allowing us to compare the extreme values of both Σ and Ω with respect to their same mean value.
The conditional expectation Σ|Ω is shown in Fig. 4a. Similar results were presented in [6], albeit without the data at R λ = 1300. As noted in [6], the conditional expectation for smaller values of Ω is approximately constant, suggesting vorticity and strain are uncorrelated. Thereafter, as Ω increases, the expectation approximately increases as a power law, Σ|Ω ∼ Ω γ , with γ < 1implying that strain acting on intense vortices is comparatively weaker than the vorticity. The exponent γ, shown in the inset of Fig. 4a, increase very weakly with R λ , suggesting that γ = 1 would be realized as R λ → ∞. Assuming a simple sigmoid or power law dependence of γ on R λ suggested that γ ≈ 1 would only be realized for extremely large R λ (which are practically unrealizable in both experiments and numerical simulations). Overall the new data at R λ = 1300 appear to be consistent with conclusions drawn in [6] (which utilized data at R λ ≤ 650). In this regard, we make a note that the power law representation Σ|Ω ∼ Ω γ , is an empirical observation, but nonetheless very useful from a modeling perspective. This power law behavior appears to be very robust, as indicated in Fig. 4b, which shows the conditional expectations compensated by Ω γ -with curves for all R λ showing a robust plateau (with variations of less than 5%).
The conditional relation between strain and vorticity is further analyzed by considering the conditional expectations of eigenvalues, i.e., λ i |Ω . Figure 5a shows the conditional expectations for i = 1 and 2. The properties of the third eigenvalue can be easily inferred from the relation λ 3 = −(λ 1 + λ 2 ). As evident from Fig. 5a, we find that the conditional expectation of λ 2 is always positive, throughout the entire range of Ω and smaller than that of λ 1 , approximately by a constant factor ≈ 8 (note that λ 1 is by construction always positive). Furthermore, both eigenvalues show the same qualitative behavior, matching that of Σ|Ω . This is not surprising since Σ can be written as Σ = 2(λ 2 1 + λ 2 2 + λ 2 3 ). Evidently, the contribution to Σ is dominated by λ 1 and λ 3 since their magnitudes are significantly larger than that of λ 2 . Consequently, one can expect λ 1 (and -λ 3 ) to scale as Σ 1/2 (even when conditioned on Ω). Although not explicitly shown, we indeed find this to be the case, with λ 1 |Ω and −λ 3 |Ω approximately scaling as Ω γ/2 .
Using similar arguments as above, one could possibly expect λ 2 |Ω to also scale as Ω γ/2 . However, this is not the case and the approximate power law we find corresponds to an exponent smaller than γ/2. This can be readily observed by considering either the ratio λ 2 /λ 1 or the parameter β defined in Eq. (6), which is bounded between −1 and 1. The conditional expectation β|Ω is shown in Fig. 5b. For Ωτ 2 K 1, β does not vary much, but for larger Ω, β appears to slowly decrease, almost independently of R λ . Thus, it appears reasonable to assume λ 2 |Ω ∼ Ω γ/2−δ . Our data suggest δ ≈ 0.05 (with the coefficient of determination exceeding 99%).
To further investigate the behavior of the intermediate eigenvalue, we show the conditional PDFs of β in Fig. 5c for R λ = 1300 (with the marginal PDF in black dashed line). Consistent with β|Ω , the PDFs show very minor variations for small Ωτ 2 K . For larger Ωτ 2 K , a slightly stronger variation can be seen, with the PDFs peaking at slightly smaller positive values of β. Interestingly, for all conditioning values, there is virtually no change in PDFs for negative value of β. This suggests that the likelihood of λ 2 to be negative (or positive) stays the same, with just positive values becoming comparatively smaller in magnitude at large Ω. This suggests that for intense vorticity events, there is a somewhat stronger cancellation between λ 1 and λ 3 . This is consistent with the notion that the regions where vorticity is most intense are essentially two-dimensional [25], with the largest and the smallest eigenvalues of strain being almost canceling each other, the value of the intermediate eigenvalue being much smaller than the two other ones. This is corroborated by the observation that the exponent of the power law for λ 2 |Ω is smaller than for λ 1 |Ω (and −λ 3 |Ω ).

C. Enstrophy production
We now investigate the net enstrophy production resulting from the correlations between strain and vorticity discussed earlier. Figure 6a shows the conditional average ω iωj s ij |Ω (made dimensionless by multiplying by τ K ), which simply equals ω i ω j s ij |Ω /Ω, since Ω is the conditioning variable. The conditional average ω iωj s ij |Ω can also be interpreted as the effective strain acting on vorticity, which ultimately leads to production of enstrophy. The dashed line shown in Fig. 6a represents a Ω 1/2 dependence, which would correspond to strain as strong as the local vorticity. Figure 6a shows that for very small values of Ωτ 2 K , ω iωj s ij |Ω grows faster than ∝ Ω 1/2 -despite earlier result in Fig. 4 showing that the magnitude of strain conditioned on vorticity is mostly constant for small Ω. This can be understood by realizing first that for Ω → 0, the conditional production term also goes to zero, since the strain and vorticity are completely decorrelated. However, as Ω increases, while the strain and vorticity magnitude are still decorrelated, there is some weak preferential alignment between vorticity and first two eigenvectors of strain (as noticed in Fig. 3). This incipient alignment is in precise agreement with the strong initial growth of production term.
On the other hand, for large Ω, the alignments are approximately in an asymptotic state, whereas Σ|Ω is weaker than Ω 1 -suggesting a similar behavior ω iωj s ij |Ω . The results shown in Fig. 6a are consistent with this expectation. To better show the dependence on R λ , Fig. 6b zooms on the domain Ωτ 2 K 1. The R λ dependence appears to be qualitatively similar to that also observed for strain and its eigenvalues in Section IV B. However, given the different power law behaviors of λ 1 |Ω and λ 2 |Ω , it is reasonable to model the curves in Fig. 6b as a sum of these two power laws: Note, such a functional form assumes that the conditional alignments are approximately constant with Ω (for Ωτ 2 K 1), which is justified from the observation in Fig. 3d. In earlier works [22,23], this dependence was modeled as a single power law, leading to some inconsistencies, until an ad hoc tuning of the model was performed. Our data suggest c 1 ≈ 0.088 and c 2 ≈ c 1 /20 (with the coefficient of determination exceeding 99%).
Further insight into enstrophy production is provided by Fig. 6c, which shows the fractional contribution in each eigendirection to the total production (based on Eq. (5)). Since the contribution corresponding to λ 3 is always negative, we plot its absolute magnitude, for explicit comparison with the contributions from other eigenvalues on log-log scales. Two very different behaviors corresponding to very weak and very large events respectively are observed. For weak events, the individual relative contributions are all significantly larger than unity, implying a strong cancellations to produce the net production of enstrophy (we recall that the individual contributions sum up to unity). This is also consistent with the observation above that the lack of alignment between vorticity and the eigenvalues of strain at small values of Ωτ 2 K results in a sum of the contributions of the three eigenvectors ≈ 0 i.e. much weaker than each individual term. A consequence of the very weak alignment of vorticity with any of the three strain eigenvectors (see Fig. 3d) is that these individual contributions can be estimated to be in the ratio of the magnitude of eigenvalues, which indeed appears to be the case. However, as the value of Ωτ 2 K increases, vorticity aligns preferentially with some particular eigenvectors of strain, and the individual contributions to the net productions change accordingly.
To better understand the R λ dependence, Fig. 6d shows a zoomed in version of Fig. 6c, focusing only on the first and second eigendirections. The relative contribution to enstrophy production from the first eigendirection start decreasing in magnitude when Ωτ 2 K increases, whereas the positive contribution from the intermediate eigendirection starts increasing. For extreme events, the contributions from the first and the second eigendirections appear to be comparable (whereas that from the third direction is significantly smaller in magnitude). It turns out that the contribution from the second direction is slightly larger than that from first for moderate events. However, this difference clearly becomes smaller as R λ increases. On the other hand, for most intense events, the dominant contribution surprisingly always comes from the first direction (and also appears to increase with the conditioning value, whereas the contribution for the second direction appears to very slowly decrease).

D. Viscous destruction of enstrophy
While vortex stretching on an average leads to enstrophy production, the viscous terms lead to enstrophy dissipation, with their effects balancing each other in the (statistically) stationary state. Of the two viscous terms in Eq. (3), the term ǫ w = 2ν ∂ω i ∂x j ∂ω i ∂x j contributes to enstrophy dissipation, whereas the other term ν∇ 2 Ω is on average zero (because of homogeneity) and only leads to spatial redistribution of enstrophy. However, when conditioning on Ω, the second term may contribute to the enstrophy budget differently, depending on the strength of vorticity [38]. Figures 7a and b show the conditional averages of ǫ w and −ν∇ 2 Ω respectively, appropriately non-dimensionalized by τ 3 K . While ǫ ω is always positive, the requirement ∇ 2 Ω = 0 imposes that the quantity −ν ∇ 2 Ω|Ω cannot be positive definite. For this reason and also since Fig. 7b utilizes log scales, positive values of −ν ∇ 2 Ω|Ω , which correspond to destruction of enstrophy, are represented by solid lines, whereas negative values corresponding to production of enstrophy, are shown in dashed lines.
As performed in previous subsections, we can demarcate the results based on the strength of vorticity. As evident in Fig. 7, both terms are roughly constant at small Ω conditioning values and of virtually identical magnitude, but of opposite signs. This suggests that the overall contribution of viscous terms at small Ω values is negligible. For Ωτ 2 K 1, both terms are positive and strongly increase with Ω. The latter is consistent with the notion that the viscous diffusion term overall acts to deplete enstrophy in regions of high enstrophy and thereafter diffuse it to regions of weaker enstrophy (resulting in no net contribution to average enstrophy budget). Similar to Fig. 4, the curves for ǫ ω |Ω appear to get steeper with increasing R λ . On the other hand, the curves for −ν ∇ 2 Ω|Ω remarkably do not appear to change slope with increasing R λ (but are slightly shifted, with the zero-crossing point weakly dependent on R λ ).
Since both the viscous terms are dimensionally equivalent to the enstrophy production term, which in turn is equivalent to third moment of vorticity, we compare the curves in Fig. 7a and b with dashed line representing Ω 3/2 . Similar to enstrophy production term (in Fig. 6a and b), the growth of ǫ w |Ω is slower than Ω 3/2 , though similar trends with R λ are observed in both cases with the curves becoming steeper and slightly closer to Ω 3/2 . On the other hand, for −ν ∇ 2 Ω|Ω , all the curves appear to conform to Ω 3/2 scaling. Interestingly, we observe that while both conditional viscous terms are of comparable magnitude, the diffusion term dominates at large Ω values, suggesting that destruction of enstrophy is dominated by the diffusion term for the most extreme events.
To interpret the results of Fig. 7, we recall that ν∇ 2 Ω describes a diffusion process, transferring enstrophy from regions where it is large to regions where it is smaller. Figure 7 allows us to distinguish two ranges, for Ωτ 2 K 1, and for Ωτ 2 K 1. The ostensibly simple dependence of ν∇ 2 Ω|Ω ∝ Ω 3/2 for Ωτ 2 K 1 can possibly be explained as follows. Using the observation that the largest vorticity fluctuations correspond to narrow weakly curved vortex tubes, we can approximately write ν∇ 2 Ω ∼ νΩ/ℓ 2 , where ℓ is the characteristic radius of the vortex tube. Thereafter, taking νΩ/ℓ 2 ∼ Ω 3/2 implies that Ω 1/2 ℓ 2 /ν ∼ 1. Thus, the vortices that contribute to viscous diffusion have a size such that the local Reynolds number is of order unity. Interestingly, a similar argument is often utilized in defining a local fluctuating dissipation length scale, akin to a multifractal description [1]. However, note the current argument stems from considering the dissipation of vorticity, and is conceptually different than the other one, which relies on a local energy balance (and does not adequately describe the intense vortices [4,6]).

V. CONCLUSIONS
The vortex stretching mechanism, resulting from non-linear coupling of vorticity ω i and strainrate tensor s ij , plays a central role in the formation of extreme events and small-scales in turbulent flows. In this paper, we have systematically investigated the statistical correlations underlying the vortex stretching mechanism in direct numerical simulations of stationary isotropic turbulence across a wide range of Reynolds numbers (140 ≤ R λ ≤ 1300). As commonly done, the correlations are studied in the eigenframe of the strain-rate tensor -for which, by definition, the first and third eigenvalues are positive and negative respectively. In this context, it is well known that vorticity preferentially aligns with the eigenvector of strain corresponding to the intermediate eigenvalue, λ 2 , which is on the average positive and thus leads to net production of enstrophy (vorticity-squared) This behavior is very robust across all R λ , as shown in Fig. 1. However, as suggested before [14], this does not necessarily imply that the dominant contribution to production of enstrophy also arises from the intermediate eigendirection. We confirm this and show that the dominant contribution arises from the first eigendirection, since the corresponding eigenvalue is significantly larger than the intermediate eigenvalue. However, the contribution from the second (intermediate) eigendirection is also comparable.
In order to understand the formation and structure of extreme events [6], we investigated the coupling between vorticity and strain tensor by extracting statistics conditioned on the enstrophy Ω = ω i ω i . The alignment between vorticity and the eigenvector corresponding to the intermediate eigenvalue of strain is found to be much stronger than the averaged alignment shown in Fig. 1a when Ωτ 2 K ≫ 1, where τ K is the Kolmogorov time scale defined as τ 2 K = 1/ Ω . This trend is consistent with visualizations of the most intense vorticity structures in turbulent flows, which turn out to be almost straight vortex tubes; and also with many attempts to follow numerically the formation of very strong velocity gradients in the Euler equations, which also revealed a strong tendency to form quasi two-dimensional structures, where vorticity can only be aligned with the intermediate strain eigenvector. Due to this enhanced alignment, we find that the contribution of second eigendirection is larger than that from first for moderately strong vorticity events, however the difference becomes smaller as R λ increases. In contrast, for the most extreme events (Ωτ 2 K 100), the contribution from first eigendirection is surprisingly larger than that from second.
Interestingly, we find that the conditional averages of strain magnitude (given by Σ = 2s ij s ij ) and the eigenvalues of strain, are mostly constant for Ωτ 2 K < 1, and approximately vary power laws for Ωτ 2 K > 1. The exponent γ, which describes the power law for Σ, i.e., Σ|Ω τ 2 K ∼ (Ωτ 2 K ) γ , also describes the dependence of the first and third eigenvalue, where γ is a function of R λ , in accordance with [6]. However, a minor correction, independent of R λ , is necessary for the intermediate eigenvalue. Combining the results, the enstrophy production can also be expressed by a simple functional form, involving the sum of two power laws, see Eq. (9). In contrast, an analysis of the viscous terms responsible for destruction of enstrophy reveals that the primary contribution for large enstrophy comes from viscous diffusion of enstrophy (ν∇ 2 Ω), rather than viscous dissipation (ν||∇ω|| 2 ). Whereas for small enstrophy the diffusion term actually leads to enstrophy production (since ultimately this term does not contribute to the overall enstrophy budget). We highlight the importance of these results in light of turbulence modeling (albeit an explicit attempt is left for future work).
In discussing our results, we have used several times the remark that the most intense vortex structures are shaped as weakly curved vortex tubes [4,6,19], which, as noticed by [25], may provide a kinematic explanation for the observed alignment between vorticity and the eigenvectors of strain. Our results, obtained by conditioning the statistics on enstrophy, appear to corroborate this picture. However, some of the other results uncovered in this work, such as the conditional expectation of β, cannot be simply explained by assuming that the regions with most intense vorticity and strain are due to weakly curved Burgers vortices. This suggests existence of a more complicated structure of the vortex tubes, involving also axial velocity [39]. In this light, it is also important to note that the analysis performed here purely relies on local single-point correlations. However, vorticity and strain are non-locally related to each other through the Biot-Savart relation [40]. Thus, a complete picture of vortex stretching can only emerge by also characterizing the nonlocal effects. Such an analysis is ongoing and will be reported separately.
We conclude by noticing that, although we have observed over the range of Reynolds numbers considered (140 ≤ R λ ≤ 1300) a systematic dependence on R λ of the statistical properties of the velocity gradient tensor conditioned on enstrophy or strain, the intriguing question remains as to whether the quantities studied here tend to a limiting form when R λ → ∞.