Structure and role of the pressure Hessian in regions of strong vorticity in turbulence

Abstract Amplification of velocity gradients, a key feature of turbulent flows, is affected by the non-local character of the incompressible fluid equations expressed by the second derivative (Hessian) of the pressure field. By analysing the structure of the flow in regions where the vorticity is the highest, we propose an approximate expression for the pressure Hessian in terms of the local vorticity, consistent with the existence of intense vortex tubes. Contrary to the often used simplification of an isotropic form for the pressure Hessian, which in effect inhibits vortex stretching, the proposed approximate form of the pressure Hessian enables much stronger vortex stretching. The prediction of the approximation proposed here is validated with results of direct numerical simulations of turbulent flows.


Introduction
The main difficulties to understand three-dimensional incompressible turbulent flows arise from two intrinsic properties of fluid mechanics, namely nonlinearity and non-locality.Nonlinearity leads to the generation of small-scale structures in turbulent flows and to the intermittent formation of highly localized and intense velocity gradients (Taylor 1938).
It is well established that strong vorticity regions emerge, forming tube-like structures (Siggia 1981;She, Jackson & Orszag 1990;Jackson, She & Orszag 1991;Jiménez et al. 1993;Ishihara, Gotoh & Kaneda 2009;Buaria et al. 2019), surrounded by sheet-like structures of strong strain regions (Buaria et al. 2019;Buaria, Pumir & Bodenschatz 2021).Interestingly, the structural difference between regions of intense strain and intense vorticity results in a strong asymmetry: vorticity is large in regions of intense strain, but in comparison strain is not as strong in regions of intense vorticity (Buaria & Pumir 2022).On the other hand, the non-locality in turbulent flows is captured by the pressure field p(x), which satisfies the Poisson equation: where m, the velocity gradient tensor, is defined by m ij ≡ ∂u i /∂x j , u being the fluid velocity.The classical solution of (1.1) leads to an expression of the pressure at any point as an integral of the velocity gradient over the whole space.The pressure field in turn enters the evolution equations for the velocity gradient tensor: In (1.2), ν is the kinematic viscosity, and h p ij ≡ ∂ 2 p/∂x i ∂x j the pressure Hessian tensor.A detailed understanding of the pressure Hessian is therefore a prerequisite to describe the dynamics of the velocity gradient tensor.Earlier studies have investigated the role of pressure Hessian in (1.2) (Borue & Orszag 1998;Tsinober, Ortenberg & Shtilman 1999;Carbone, Iovieno & Bragg 2020;Zhou & Yang 2023), or in the equation for the coarse-grained velocity gradient tensor (Chertkov, Pumir & Shraiman 1999;Yang, Pumir & Xu 2020;Tom, Carbone & Bragg 2021).An important issue consists in understanding the role of velocity gradient structure on pressure (Pumir 1994;Vlaykov & Wilczek 2019), and more generally, in developing closure models for the pressure Hessian in the Lagrangian dynamics of velocity gradient (Vieillefosse 1982(Vieillefosse , 1984;;Cantwell 1992;Martın, Dopazo & Valiño 1998a;Martın et al. 1998b;Chevillard et al. 2008;Meneveau 2011;Wilczek & Meneveau 2014;Lawson & Dawson 2015;Johnson & Meneveau 2016;Tian, Livescu & Chertkov 2021).Here, our general objective is to propose an approximate expression for the pressure Hessian, using a closed expression in terms of the velocity gradient, with the aim of providing insight on the dynamics of the velocity gradient tensor itself.
The simplifying assumption, known as restricted Euler (RE), consists in considering an isotropic functional form for h p , i.e. h p ij = 1 3 ∇ 2 pδ ij , with δ ij being the Kronecker symbol and ∇ 2 p given by (1.1).This approximate form expresses h p locally in terms of strain and vorticity, and is formally obtained by neglecting non-local contributions when expressing the solution of (1.1) (Ohkitani & Kishiba 1995).Remarkably, this leads to a simple dynamical system which can be integrated with elementary means (Vieillefosse 1982(Vieillefosse , 1984;;Cantwell 1992), therefore providing a qualitative dynamics, in a statistical sense, in the (R, Q) plane, where R ≡ − 1 3 tr(m 3 ) and Q ≡ − 1 2 tr(m 2 ) are two invariants of m.In particular, the simplified dynamics leads to a singularity along the separatrix 4Q 3 + 27R 2 = 0 and R > 0, which provides an explanation for the excess of probability of the invariants along this separatrix (Cantwell 1993).However, the RE assumption predicts that in vorticity-dominated regions (Q > 0), the flow moves quickly away from regions with positive vortex stretching (R < 0) towards regions where vortex stretching is negative.In comparison, the results of direct numerical simulation (DNS) show that the approximate form of the pressure Hessian obtained from the RE approximation leads to an excessively strong inhibition of vortex stretching (Tsinober 2009;Buaria & Pumir 2023).This clearly implies that the deviatoric part of h p ij , denoted H p ij in the following, and defined as the difference between the pressure Hessian and the isotropic component, , plays an essential role.The difficulty, however, is that H p ij is determined by an integral expression (Ohkitani & Kishiba 1995).
Here, we focus on regions with strong vorticity, and propose that the pressure Hessian tensor is determined primarily by local vorticity values due to the asymmetric distribution of vorticity and strain.We will present our theory and propose that where ω i ≡ ijk m jk is the vorticity ( ijk is the Levi-Civita tensor) and Ω ≡ ω i ω i is twice the enstrophy, and τ K ≡ 1/ Ω 1/2 denotes the Kolmogorov time scale.Throughout this text, • denotes the ensemble average of a fluctuating quantity.We then discuss the implications of this result on the Lagrangian dynamics of the velocity gradient.Our analysis shows that the contributions from the deviatoric parts of pressure Hessian are not negligible and counteract the contribution of nonlinear terms to the leading order, and leads to an enhancement of vortex stretching in moderately strong vorticity regions.

The DNS data
We first introduce the two DNS datasets of homogeneous and isotropic turbulence used in this work.Dataset A is downloaded from snapshot 5 of the 'Forced isotropic turbulence dataset on 8192 3 Grid' dataset (Yeung,  Then the statistic for Ωτ 2 K = 25 would be 2.3 × 10 6 , and 2.2 × 10 6 for Ωτ 2 K = 100.For comparison, we also downloaded 512 3 equally spaced data points of the whole snapshot, which would provide information for statistics of unconditional and lower values of vorticity. Dataset B was obtained by simulating the Navier-Stokes equations using the spectral code described in Pumir (1994) with 512 3 grid points and R λ = 210, with a lower spatial resolution compared with Dataset A: k max η K ≈ 1.6.While this resolution is insufficient to reliably capture the statistics of the largest velocity gradients in the flow, we explicitly checked that the probability density function (p.d.f.) of enstrophy, Ω, coincides with the data shown in figure 2 of Buaria et al. (2019) at R λ = 240, for Ωτ 2 K 50.Other tests convinced us that resolution is not an issue over the range of velocity gradient intensity considered here.We saved 21 configurations of the full velocity fields, over a time of the order of ∼10 eddy turnover times, which we analysed to produce the data shown here.

Analysis of the pressure Hessian in regions of very large vorticity
We begin by focusing on the antisymmetric part of the velocity gradient tensor in (1.2), which involves vorticity ω and thus the enstrophy Ω: where is the rate of strain tensor.Equation (3.1) implies that enstrophy amplification results from a competition between vortex stretching ω i s ij ω j and the viscous diffusion.Notoriously, the pressure Hessian does not explicitly appear in the equations for ω i and Ω.On the other hand, the pressure Hessian ∂ 2 p/∂x i ∂x j plays an explicit role in the dynamic equation of strain s ij : and, therefore, pressure Hessian will affect the vortex stretching dynamics by acting on strain (Buaria & Pumir 2023).

Heuristic derivation
We start by conditioning the equation of the dynamics of strain, (3.2), on enstrophy Ω and taking the conditional average: Earlier work has shown that vorticity dominates the strain in extreme vorticity regions, in the sense that Σ|Ω ∼ Ω γ , with Σ ≡ 2s ij s ji , and with γ < 1 when Ωτ 2 K 1 (Buaria et al. 2019;Buaria, Bodenschatz & Pumir 2020;Buaria & Pumir 2022).We use the corresponding feature, namely Σ Ω, to analyse the magnitudes of the various terms appearing in (3.3).This is at the origin of the simple functional form for h p ij that we propose in this work.
We first notice that both the operators D/Dt and ν∇ 2 have the dimension of the inverse of a time, 1/T.Given that the fastest time scale in extreme vorticity regions is 1/Ω 1/2 , we estimate that Besides, the second and third terms on the left-hand side of (3.3) scale as Thus the third term dominates the first, second and fifth terms on the left-hand side of (3.3), as a result, we should have h 1.This argument suggests that the pressure Hessian is determined by the local value of vorticity in extreme vorticity regions.For simplicity, we denote and, subtracting out the trace, we obtain its deviatoric part as The approximations leading to (3.6) suggest that the difference between H p ij and H p ij is mostly due to the strain component, which is explicitly neglected in our heurstic derivation.This is consistent with our own numerical results, see e.g.figure 1, which indicate that the average of tr[(H p − H p ) 2 ] conditioned on Ω grows as ∼Ω 2γ , whereas the average of tr[(H p ) 2 ] conditioned on Ω grows as ∼ Ω 2 .This implies, in particular, that the p.d.f. of the difference between H p and its asymptotic form is becoming narrower, in units of tr[(H p ) 2 ]|Ω when Ω is very large.

Structure of the pressure Hessian in regions of large vorticity
Before exploring the implications of (3.6) on the dynamics of the velocity gradient tensor, we first discuss the properties of the tensor H p , and compare them with the results of DNS.
It is easy to see that the tensor has a zero eigenvalue in the eigendirection ω, along with a degenerate eigenvalue 1 4 ω 2 , in the plane perpendicular to ω.This is generally consistent with DNS results.Nomura & Post (1998), Lawson & Dawson (2015) and Buaria & Pumir (2023)  4 ω 2 and 0; the intermediate one reaching a value slightly smaller than 1 4 ω 2 .To investigate further the structure of the eigenvalues of h p , we study the joint p.d.f. of its invariants.Following Tom et al. (2021), we define where b is the deviatoric part of the pressure Hessian, H  We now demonstrate that the structure of the p.d.f.s of (ζ, χ ) conditioned on Ωτ 2 K 1 is fully consistent with the presence of intense vortex tube-like structure (Siggia 1981;She et al. 1990;Jiménez et al. 1993;Ishihara et al. 2009;Buaria et al. 2019).To this end, we consider a Burgers vortex, characterized by the standard definition of the Reynolds number, Re = Γ /ν, where Γ is the circulation of the vortex.At the centre of the vortex (Burgers 1948) with Re 1, the eigenvalues of the pressure Hessian are also 1 4 ω 2 , 1 4 ω 2 and 0; and the eigendirection corresponding to the 0 eigenvalue aligns with ω (Andreotti 1997;Horiuti 2001;Horiuti & Takagi 2005).We now compare the structure of the pressure Hessian in a turbulent flow with that in a Burgers vortex, whose properties depend only on the distance r to the axis.Following Andreotti (1997), we use the dimensionless radius r * = r a/ν, where a is the external strain rate acting on the vortex.We note that the Reynolds number Re amounts to the vorticity at the centre divided by the external strain a.As already stated, the maximum value of the probability in all panels of figure 2, indicated by the green squares, corresponds to the axisymmetric configuration on the axis of the vortex (r * = 0).Furthermore, the dependence of the invariants ζ and χ on r * up to r * = 1 in a Burgers vortex is shown as the red line with crosses in figure 2(d), see also (4) of Andreotti (1997).As r * increases, the solution moves up in the (ζ, χ ) plane, first towards the left border and then back towards the centre.We note that its variation with respect to Reynolds number Re is negligible.The joint p.d.f. of (ζ, χ ), particularly the ridges, closely follows the Burgers vortex solution.We also observe that when conditioning on a given value of Ωτ 2 K , as done in figure 2, the points at a given value of (ζ, χ ) away from the maximum, or ζ = χ = √ 3/9, correspond to points at some distances away from the centre, or r * > 0, of vortices with increasing intensity.These vortices are therefore less probable, which explains that the maximum probability in figure 2 corresponds to ζ = χ = √ 3/9.Additionally, the distribution becomes more concentrated along the red curve when Ωτ 2 K increases.Figure 2 therefore demonstrates that the Burgers vortex model provides a compelling qualitative description of the pressure Hessian.We notice, however, that any axisymmetric vortex model would lead to qualitatively similar observations near the centre of the vortex, consistent with (3.6) for the deviatoric part of the pressure Hessian tensor.

Dynamics in the (R, Q) plane
To illustrate the dynamical implications of (3.6), we project the equations of motion for m on the plane (R, Q), where the invariants R and Q are defined by In (4.1) and (4.2), and throughout this text, we assume summation of repeated indices.Physically, Q expresses the difference between enstrophy and the square of the strain rate, so Q = 1 2 Ω − tr(s 2 ) > 0 corresponds to vorticity-dominated regions and Q < 0 corresponds to strain-dominated regions.Similarly, R = − 1 3 tr(s 2 ) − 1 4 ω • s • ω < 0 corresponds to the vortex-stretching-dominated region.The exact Lagrangian evolution equations for the invariants R and Q can be readily derived from (1.2): where H p ij is the deviatoric part of the pressure Hessian, and We stress that these two terms in (4.3) and (4.4) are unclosed.The RE approximation introduced above provides a closure consisting in neglecting the deviatoric contributions to the pressure Hessian: H p = 0.As mentioned already (Buaria & Pumir 2023), the RE assumption fails to predict the qualitative behaviour of the dynamics in the (R, Q) plane when Q is positive and large.Specifically, the RE assumption leads to dR/dt = 2 3 Q 2 > 0, but the results of DNS show that dR/dt < 0 in some region of the second quadrant of the (R, Q) plane, see, for example, figures 3(a) and 3(d) of Wilczek & Meneveau (2014) and figure 5 of Yang et al. (2023).The observation that pressure leads to a 'depression of nonlinearity' (Borue & Orszag 1998) suggests the functional where 0 < α < 1 is a model parameter (Chertkov et al. 1999).This model predicts that the pressure terms induced by H p ij in (4.3) and (4.4) anti-align with the RE 983 R2-7 terms in the (R, Q) plane, in quantitative contradiction with numerical results for Q > 0, and particularly for R < 0. Attempts to improve the RE approximation have assumed mostly some local expressions for the pressure Hessian, with an explicit relation between H p and the velocity gradient tensor m (Wilczek & Meneveau 2014).In the following, we will discuss theoretically the role of pressure Hessian in the vorticity-dominated regions, and derive a local analytical expression for the contributions of H p ij in Eqs.(4.3) and (4.4), to leading order when Ωτ 2 K 1.We now use the approximate form of the deviatoric part of the pressure Hessian, (3.6), to analyse the dynamics in the (R, We recall that the expressions for dQ/dt and dR/dt obtained by using the RE approximation, provided by (4.3) and (4.4) are −3R and 2 3 Q 2 , respectively.In the spirit of the approach in § 3.1, based on the quantitative difference between Σ and Ω, namely with Σ Ω when Ωτ 2 K 1, we further notice that when Ωτ 2 K 1, R 1 3 s ij s jk s ki , so the H p term cancels, to the leading order in Σ/Ω, approximately 1/3 of the RE term in the Q equation.More significantly, in the limit The lowest-order contribution in a formal development in (Σ/Ω) involves terms such as − 1 12 ω 2 s ij s ji − 1 4 ω i s ik s kj ω j , which is negative.We observe that this implies dR/dt < 0, which is generally consistent with the DNS results (see e.g.figure 5 of Yang et al. 2023).Earlier work aimed at describing the pressure Hessian (Chevillard et al. 2008;Wilczek & Meneveau 2014) did not particularly focus on high-enstrophy regions.It would be interesting to compare the corresponding predictions in light of the approximation.We note that the tetrad model (Chertkov et al. 1999;Yang et al. 2020Yang et al. , 2023) ) involves quadratic terms in the vorticity introduced through the full tensor m, insufficient to compensate for the strong vortex-stretching reduction due to the RE approximation.
Figure 3(a) shows the DNS results for dQ/dt and dR/dt on the (R, Q) plane (blue arrows) and compares them with the theoretical predictions (magenta arrows).Figure 3(a) demonstrates that when Qτ 2 K 5 and |R τ 3 K | is not too large, the magenta arrows, which represent our theoretical predictions (rightmost terms of (4.5), (4.6)), agree well with the blue arrows, which represent the 'exact' DNS values.Figure 3 also directly compares the ratios between the numerically determined values tr(H p • m) (b) and tr(H p • m 2 ) (c) and the values found by replacing H p by H p .The agreement for the tr(H p • m 2 ) generally improves when Q increases, as expected when Ωτ 2 K increases.The effect is weaker for the term tr(H p • m).We note that the quality of the agreement at R ≈ 0 degrades when Qτ 2 K increases, which is due to the small values of both numerator and denominator.Figure 4 provides further insight on the dynamics, by comparing the contributions of H p to dQ/dt and dR/dt in the (R, Q) plane, with the leading-order terms in (4.5) and (4.6), namely R and − 2 3 Q 2 .Figure 4 reveals that when Qτ 2 K > 5 and |Rτ 3 K | is not too large, the green arrows, which represent our theoretical predictions to leading order are reassuringly close to the blue arrows, which represent the values obtained directly from DNS.
Remarkably, we notice that using − 2 3 Q 2 for the contribution due to the deviatoric part of the pressure Hessian leads to dR/dt = 0 in the inviscid case.Whether the blue arrows in figure 4(a) are longer or shorter than the green ones in fact determines the sign of dR/dt.We observe that dR/dt > 0 when Qτ 2 K 20 for Dataset B with R λ = 210.At lower values  of Qτ 2 K , and for Rτ 3 K sufficiently negative, on the contrary, dR/dt < 0. For these values, the dynamics tends to keep the values (R, Q) in the (R < 0, Q > 0) quadrant, where stretching is largest.This contrasts sharply with the predictions of the RE approximation, which tends to inhibit vortex stretching (Buaria & Pumir 2023).A more precise comparison between tr(H p • m 2 ) and − 2 3 Q 2 is shown in figure 4(b).At relatively low value of Qτ 2 K , and for Rτ 3 K −3, tr(H p • m 2 ) is more negative than − 2 3 Q 2 , consistent with the observation that dR/dt < 0. As Qτ 2 K increases, however, the region where dR/dt < 0 is restricted to much larger values of Rτ 3 K < 0, whereas at small, intermediate values of |Rτ 3 K |, dR/dt > 0.

Conclusion
In this work, we have proposed a very simple approximate expression for the deviatoric part of the pressure Hessian in terms of the local field, valid in regions of very large vorticity in turbulent flows.Namely, we find that ∂ 2 p/∂x i ∂x j ∼ − 1 4 (ω i ω j − ω 2 δ ij ) when ω 2 τ 2 K 1.This simple expression comes from an elementary analysis of the flow, and in particular, on the structure of the solutions (vortex tubes) in these flow regions.We stress that the deviatoric term of the pressure Hessian is essentially non-local, which implies that (3.6) can be only approximately valid.Despite these limitations, the observations from DNS support our theoretical discussions and provide further evidence for our conclusions.The numerical results in homogeneous isotropic flows presented here therefore demonstrate that the considerations developed in this work provide potentially valuable insight into the closure models of pressure Hessian in the Lagrangian dynamics of velocity gradient.Although one expects on general grounds that our conclusions should extend to other turbulent flows, e.g. in the presence of mean shear and/or walls at large Reynolds numbers, it will be interesting to check the validity of our approximation for other classes of turbulent flows.
The expression (3.6) reduces to the simple explicit (local) relation between H p ij and w ij ≡ 1 2 (m ij − m ji ), the antisymmetric part of the velocity gradient tensor: H p ij = w ik w kj − 1 3 tr(w 2 )δ ij .We note in this respect that recent machine-learning approaches (Tian et al. 2021;Buaria & Sreenivasan 2023) have looked systematically for an expansion of H p in terms of all the invariants of the velocity gradient tensor, m.In this sense, (3.6) can be compared with that from machine learning approaches.Alternatively, the expression provides a point of comparison which may be used in a physics-informed machine-learning approach.It is our expectation that the physical considerations initiated in this work can be further developed, and can lead to clearer understanding of the structure and modelling of turbulent flows.
observed that when Ωτ 2 K 1, ω strongly aligns with eigendirection e ph 3 of the pressure Hessian h p ij ; and the largest and the smallest eigenvalues of h p ij conditioned on large Ω tend to approach 1

Figure 2 .
Figure 2. Joint p.d.f. of the dimensionless invariants ζ and χ of the pressure Hessian conditioned on (a) Ωτ 2 K = 1, (b) Ωτ 2 K = 4, (c) Ωτ 2 K = 25 and (d) Ωτ 2 K = 100.The red line with cross symbols in (d) corresponds to a Burgers vortex.The white dashed-dotted lines show the iso-probability contours corresponding to the value 1, and the dashed lines to 2 n .

Figure 2
Figure 2 shows the joint p.d.f. of ζ and χ for Ωτ 2 K = 1, 4, 25 and 100, respectively.The eigenvalues of the pressure Hessian exhibit a high-probability region near the configuration corresponding to hp ij , where λ p 1 = λ p 2 = h p ii /2 and λ p 3 = 0 so χ = ζ = √ 3/9 ≈ 0.1925.This configuration is indicated by the green squares in figure 2.We now demonstrate that the structure of the p.d.f.s of (ζ, χ ) conditioned on Ωτ 2

Figure 3 .
Figure 3. (a) Results for dQ/dt and dR/dt in the (R, Q) plane.The red arrows show the RE term (−3R, 2 3 Q 2 ), the blue arrows the deviatoric pressure Hessian term (−m ij H p ji , −m ij m jk H p ki ) and the magenta arrows our theoretical predictions for large Ω, i.e. the rightmost terms of (4.5), (4.6).Panels (b,c) show the ratio between the contributions of the deviatoric part of the pressure Hessian, H p , and those of the approximate form, H p , (3.6), to dQ/dt (b) and dR/dt (c).The ratios, plotted as a function of Rτ 3 K , increase towards 1 as Qτ 2 K increases.The large value for dQ/dt close to Rτ 3 K ≈ 0 corresponds to very small values of both the terms tr(H p • m) and tr( H p • m).Dataset B was used to construct the figure.
Figure 4. (a) Pressure contribution to dQ/dt and dR/dt in the (R, Q) plane.The red arrows correspond to the RE term (−3R, 2 3 Q 2 ), which corresponds to the isotropic (local) contribution of the pressure Hessian, the blue arrows to the anisotropic pressure Hessian term (−m ij H p ji , −m ij m jk H p ki ) and the green arrows show the simplified theoretical predictions at large Ω, (R, − 2 3 Q 2 ).Panel (b) shows the ratio between the contribution of H p to dR/dt, compared with the simplified form −2Q 2 /3.Dataset B (R λ = 210) was used to construct the figure.
of the deviatoric part of the pressure Hessian tensor, H p (blue line) and of the difference, (H p − H p ) (red line) conditioned on Ω.In agreement with the approximation proposed in the text, the second moment of H p grows as Ω 2 , whereas the second moment of (H p − H p ) grows as Ω 2γ , with γ < 1. Dataset B was used to construct the figure.
Funding.P.-F.Y., H.X. and G.W.H. are grateful to the Natural Science Foundation of China (NSFC) Basic Science Center Program 'Multiscale Problems in Nonlinear Mechanics' (grant no.11988102).P.-F.Y. and H.X. also acknowledge financial support from the NSFC grants 12202452, 11672157 and 91852104.A.P. received support from Agence Nationale de la Recherche (ANR) under grant number Contract TILT No. ANR-20-CE30-0035.Declaration of interests.The authors report no conflict of interest.