Asymptotic predictions on the velocity gradient statistics in low-Reynolds-number random flows: onset of skewness, intermittency and alignments

Abstract

Stirring a fluid through a Gaussian forcing at a vanishingly small Reynolds number produces a Gaussian random field, while flows at higher Reynolds numbers exhibit non-Gaussianity, cascades, anomalous scaling and preferential alignments. Recent works (Yakhot & Donzis, Phys. Rev. Lett., vol. 119, 2017, 044501; Khurshid et al., Phys. Rev. E, vol. 107, 2023, 045102) investigated the onset of these turbulent hallmarks in low-Reynolds-number flows by focusing on the scaling of the velocity gradients and velocity increments. They showed that the onset of power-law scalings in the velocity gradient statistics occurs at low Reynolds numbers, with the scaling exponents being surprisingly similar to those in the inertial range of fully developed turbulence. In this work we address the onset of turbulent signatures in low-Reynolds-number flows from the viewpoint of the velocity gradient dynamics, giving insights into its rich statistical geometry. We combine a perturbation theory of the full Navier–Stokes equations with velocity gradient modelling. This procedure results in a stochastic model for the velocity gradient in which the model coefficients follow directly from the Navier–Stokes equations and statistical homogeneity constraints. The Fokker–Planck equation associated with our stochastic model admits an analytic solution that shows the onset of turbulent hallmarks at low Reynolds numbers: skewness, intermittency and preferential alignments arise in the velocity gradient statistics through a smooth transition as the Reynolds number increases. The model predictions are in excellent agreement with direct numerical simulations of low-Reynolds-number flows.

1. Introduction

The Reynolds number characterizes the range of active scales involved in the flow of a Newtonian fluid, determining the transition from a laminar motion to a disordered and turbulent state (Reynolds 1883). At a vanishingly small Reynolds number, the flow obeys linear equations, and stirring the fluid through a Gaussian random forcing produces a Gaussian random velocity field. On the other hand, flows at higher Reynolds numbers undergo a nonlinear evolution and exhibit highly non-Gaussian turbulent features. These distinguishing marks of the turbulent dynamics include cascades (e.g. Alexakis & Biferale 2018; Ballouz & Ouellette 2020; Vela-Martín & Jiménez 2021), anomalous scaling of the velocity increments (e.g. Benzi et al. 1995; Chen et al. 2005), extreme intermittency and preferential alignments of the velocity gradients (e.g. Ashurst et al. 1987; Lund & Rogers 1994; Buaria et al. 2019; Buaria & Pumir 2022). Recent works (Yakhot & Donzis 2017; Sreenivasan & Yakhot 2021; Gotoh & Yang 2022; Khurshid, Donzis & Sreenivasan 2023) characterized the onset of these turbulent features and highlighted a striking similarity between the scalings in low-Reynolds-number flows and fully developed turbulence, which motivates further investigations of low-Reynolds-number random flows.

One approach to address the onset of turbulent motion consists of considering the three-dimensional, incompressible Navier–Stokes equations for a statistically isotropic flow without boundaries, driven by a large-scale Gaussian forcing. This set-up excludes the effect of any boundary condition and partially overcomes the lack of universality of low-Reynolds-number flows (Gotoh & Yang 2022). Such idealized flows can be investigated analytically at a small Reynolds number by formulating a perturbation theory of the Navier–Stokes equations (Wyld 1961). This direct approach provides complete insight into the full velocity field at low Reynolds numbers but involves several technical complications. For example, the terms in the series expansion soon become excessively complicated and violations of Galilean invariance occur (Yakhot & Donzis 2018).

A recent insightful approach focused on the scaling of the velocity gradient moments and structure functions at a low Reynolds number (Yakhot & Donzis 2017; Sreenivasan & Yakhot 2021) starting from the Hopf equation for the characteristic functional of the velocity field (Hopf 1952). This scaling analysis showed that, surprisingly, low-Reynolds-number flows without an inertial range, not detected even with the aid of the extended self-similarity (Benzi et al. 1993), have qualitatively the same scalings as fully developed turbulence (Schumacher, Sreenivasan & Yakhot 2007). In particular, the scaling exponents of the structure functions, with the Reynolds number and spatial separation, observed in low-Reynolds-number flows match well the scaling exponents predicted by a theory relying on very-high-Reynolds-number hypotheses (Yakhot & Sreenivasan 2004). Those scalings are observed at spatial separations $r\gtrapprox \eta$, where $\eta$ is the Kolmogorov length scale, and at a Reynolds number based on the Taylor microscale $\textit {Re}_\lambda \gtrapprox 9$, whereas any anomalous scaling is negligible at $\textit {Re}_\lambda \leq 3$ (Yakhot & Donzis 2017). While this scaling analysis fully characterizes the velocity increment statistics across the scales, it does not shed light on the rich statistical geometry of the flow. For example, the alignments and interplay between the strain rate and the vorticity cannot be inferred.

Here, we address the onset of non-Gaussianity in flows driven by a random forcing from the viewpoint of the velocity gradients. The velocity gradient encodes many distinguishing features of the turbulent state (Meneveau 2011), and it comprehensively characterizes the geometry of the vorticity and strain rate. As the Reynolds number increases, the velocity gradient transitions from a Gaussian random matrix state (Livan, Novaes & Vivo 2018) to a turbulent state, featuring skewness of the longitudinal components, associated with the cascade of kinetic energy (Eyink 2006; Carbone & Bragg 2020; Johnson 2020), and preferential configurations of the strain-rate eigenvalues (Betchov 1956; Lund & Rogers 1994; Davidson 2015), as well as intermittency and preferential alignments between the vorticity and the strain-rate eigenvectors (Ashurst et al. 1987; Buaria, Bodenschatz & Pumir 2020).

To analytically capture the onset of non-Gaussianity, we construct a model for the velocity gradients that is directly derived from the randomly forced Navier–Stokes equations at low Reynolds numbers. The main challenge in formulating such a model for the gradient dynamics stems from the non-locality of turbulence. The drastic reduction of degrees of freedom in going from the full Navier–Stokes equations to a small system of ordinary differential equations governing the gradient dynamics comes at the cost of introducing unclosed terms (Meneveau 2011). Those unclosed terms consist of the traceless/anisotropic pressure Hessian and the viscous Laplacian of the velocity gradient, which require modelling. Recent phenomenological models for the velocity gradient have proven effective in reproducing the small-scale turbulence statistics at moderately large Reynolds numbers (e.g. Girimaji & Pope 1990; Chertkov, Pumir & Shraiman 1999; Chevillard & Meneveau 2006; Wilczek & Meneveau 2014; Johnson & Meneveau 2016; Pereira, Moriconi & Chevillard 2018; Leppin & Wilczek 2020).

In contrast to phenomenological models at high Reynolds numbers, our low-Reynolds-number model for the velocity gradients can be derived directly from the Navier–Stokes equations. This direct derivation reduces the number of modelling hypotheses and free parameters. Indeed, at order zero in the Reynolds number, the velocity field is Gaussian, allowing for the exact computation of the non-local/unclosed terms in the equations governing the velocity gradient (Wilczek & Meneveau 2014). We use those exact expressions of the unclosed terms at zero Reynolds number to construct an expansion in the Reynolds number of the velocity gradient dynamics. Then, we close the model by using the asymptotic weak-coupling expansion of the full Navier–Stokes equations at small Reynolds number (Wyld 1961), combined with the two statistical homogeneity constraints on the incompressible velocity gradient (Betchov 1956; Carbone & Wilczek 2022). The resulting model for the single-time statistics of the velocity gradient does not feature adjustable parameters and does not require any input from simulations or experiments. Furthermore, we extend the model to predict the full temporal dynamics by using two adjustable model parameters fitted from direct numerical simulation (DNS) results. This extended model captures both the velocity gradient single-time statistics and the time correlations.

The presented model is associated with a Fokker–Planck equation (FPE) for the velocity gradient probability density function (p.d.f.), in which the Reynolds number and forcing parameters are in one-to-one correspondence with the same parameters featured in the forced Navier–Stokes equations. This FPE admits asymptotic analytic solutions, thus yielding an analytic approximation to the velocity gradient p.d.f. We also provide extensive comparisons of the analytical results to DNS data. Since the analytical results are only asymptotically valid, at a sufficiently small Reynolds number, the DNS results help to determine up to which Reynolds number our analytical results hold. We also included some DNS results at higher Reynolds numbers to illustrate which low-Reynolds-number features persist in turbulent flows.

Attempts to analytically predict the high-dimensional velocity gradient p.d.f. so far built upon phenomenological models for the small-scale turbulent dynamics (Chertkov et al. 1999; Moriconi, Pereira & Grigorio 2014; Apolinário, Moriconi & Pereira 2019). In general, the analytical expression for the p.d.f. follows from field-theoretical methods employed to solve the nonlinear Langevin equation governing the flow dynamics (Martin, Siggia & Rose 1973), consisting of the renormalized action method with a one-loop correction (Kleinert & Schulte-Frohlinde 2001; Cavagna et al. 2021). The resulting form of the p.d.f. of the gradient is typically not integrable analytically, thus preventing the direct computation of marginal distributions and moments, computed instead via Monte-Carlo sampling (Moriconi et al. 2014). Differently, our prediction of the p.d.f. of the velocity gradient is analytically integrable and it allows us to derive expressions for the marginal distributions and the moments of the velocity gradients. Those analytic expressions explicitly relate the quantities of interest, e.g. skewness, kurtosis and preferential alignments of the gradients, to the Reynolds number and forcing parameters, thus rationalizing the onset of non-Gaussianity in low-Reynolds-number random flows.

The paper is organized as follows. Section 2 presents the derivation of a FPE for the velocity gradient p.d.f. from the Navier–Stokes equations. Section 3 specializes the FPE for low-Reynolds-number flows, while § 4 presents its analytic solution. The comparison between the model/analytic predictions and low-Reynolds-number DNS is presented in §§ 5 and 6, while the conclusions and outlook are discussed in § 7. Appendices A and B describe the set-up of the numerical simulations and the low-Reynolds-number expansion of the Navier–Stokes equations used to determine the model coefficients.

2. Fokker–Planck equation for the velocity gradient probability density in statistically isotropic flows

In this section we obtain the general FPE governing the single-time/single-point statistics of the velocity gradient. We then specialize it for flows at low Reynolds numbers.

2.1. Governing equations and reference scales

We begin with the three-dimensional incompressible Navier–Stokes equations, driven by an external Gaussian stochastic forcing $\boldsymbol {F}$,

(2.1a)$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} = 0, \end{gather}$$

(2.1b)$$\begin{gather}\partial_t \boldsymbol{u} + \textit{Re}_\gamma\left[(\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla}) \boldsymbol{u} + \boldsymbol{\nabla} P \right] = \nabla^2 \boldsymbol{u} + \sigma \boldsymbol{F}, \end{gather}$$

where $\boldsymbol {u}(\boldsymbol {x},t)$ is the velocity field, $P(\boldsymbol {x},t)$ is the pressure field and $\textit {Re}_\gamma$ is the Reynolds number. The equation governing the velocity gradient dynamics is obtained by taking the gradient of (2.1),

(2.2a)$$\begin{gather} \textrm{Tr}(\boldsymbol{\mathsf{A}}) = 0, \end{gather}$$

(2.2b)$$\begin{gather}\partial_t \boldsymbol{\mathsf{A}} + \textit{Re}_\gamma\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{\mathsf{A}} ={-} \textit{Re}_\gamma\left(\boldsymbol{\mathsf{AA}} + \boldsymbol{\mathsf{H}} \right) + \nabla^2 \boldsymbol{\mathsf{A}} + \sigma \boldsymbol{\nabla} \boldsymbol{F}, \end{gather}$$

where $A_{ij} = \boldsymbol {\nabla }_j u_i$ is the velocity gradient, $H_{ij} = \boldsymbol {\nabla }_i\boldsymbol {\nabla }_j P$ is the pressure Hessian and standard matrix product is implied. The Gaussian tensorial noise $\boldsymbol {\nabla } \boldsymbol {F}$ in (2.2) has zero mean, is statistically isotropic and white in time. Its single-point statistics are fully specified by its single-point correlation, which in Cartesian component notation reads

(2.3)\begin{equation} \left\langle \boldsymbol{\nabla}_j F_i (\boldsymbol{x},t) \boldsymbol{\nabla}_q F_p(\boldsymbol{x},t') \right\rangle = \delta(t-t')\left( 4\delta_{ip}\delta_{jq}-\delta_{iq}\delta_{jp}-\delta_{ij}\delta_{pq} \right), \end{equation}

where $\delta _{ij}$ denotes the Kronecker delta. The two-point correlation of the stochastic forcing is detailed in Appendix A.

Equations (2.1) and (2.2) are written in non-dimensional variables suited for the upcoming low-Reynolds-number expansion. We employ non-dimensional variables throughout the paper, and denote the corresponding dimensional variables with a bar when needed. To construct reference scales for a suitable non-dimensionalization of the variables, we introduce a reference length $\bar {\gamma }_0$, related to the damping role of the viscous Laplacian at very low Reynolds numbers. When the velocity field is multipoint Gaussian, the conditional Laplacian of the velocity gradient reduces to a linear damping (Wilczek & Meneveau 2014), and the corresponding damping length scale is (in dimensional variables)

(2.4)\begin{equation} \bar{\gamma}_0 = \sqrt{ -\frac{ \left\langle \bar{A}_{ij}\bar{A}_{ij} \right\rangle}{\left\langle \bar{A}_{ij}\bar{\nabla}^2\bar{A}_{ij} \right\rangle}}. \end{equation}

The length scale $\bar {\gamma }_0$ depends on the spatial correlation of the forcing, and we determine it in Appendix B. Using $\bar {\gamma }_0$, together with the kinematic viscosity of the fluid $\bar {\nu }$ and the Kolmogorov time scale $\bar {\tau }_\eta =1/\sqrt {\langle \|\bar {\boldsymbol{\mathsf{A}}}\|^2\rangle }$ (angle brackets indicating ensemble average), we define the non-dimensional variables employed in (2.1) as

(2.5a–d) \begin{equation} t = (\bar{\nu}/\bar{\gamma}_0^2)\bar{t}, \quad \boldsymbol{x} = \bar{\boldsymbol{x}}/\bar{\gamma}_0, \quad \boldsymbol{u} = (\bar{\tau}_\eta/\bar{\gamma}_0)\bar{\boldsymbol{u}}, \quad \sigma^2 = (\bar{\gamma}_0^2\bar{\tau}_\eta^2/\bar{\nu})\bar{\sigma}^2. \end{equation}

Based on these quantities, we can define a Reynolds number according to

(2.6)\begin{equation} \textit{Re}_\gamma = \frac{\bar{\gamma}_0^2}{\bar{\nu}\bar{\tau}_\eta}, \end{equation}

which expresses the ratio between the velocity gradient magnitude $\bar {\tau }_\eta ^{-1}$ and the viscous damping $\bar {\nu }/\bar {\gamma }_0^2$. The Reynolds number (2.6) weighs the nonlinearities in (2.1), thus allowing us to take the small-Reynolds-number limit properly. The zeroth-order solution of the non-dimensional Navier–Stokes equation (2.1) consists of a Gaussian random velocity field resulting from a stochastically driven diffusion equation. By gradually increasing the Reynolds number, the nonlinear terms come into play leading to non-Gaussian statistics.

In the following, we analyse this ‘onset of non-Gaussianity’ at low Reynolds numbers as $\textit {Re}_\gamma$ increases. In previous works (e.g. Yakhot & Donzis 2017; Khurshid et al. 2023) the onset of non-Gaussianity has been investigated by means of the more common Reynolds number based on the Taylor microscale, $\textit {Re}_\lambda$. In Appendix B we show that at low Reynolds numbers, $\textit {Re}_\lambda$ and $\textit {Re}_\gamma$ are proportional, namely

(2.7)\begin{equation} \textit{Re}_\lambda \simeq 1.5 \textit{Re}_\gamma, \end{equation}

with the proportionality factor being weakly dependent on the correlation of the stochastic forcing (the reported value refers to our simulation set-up). The moderate dependence of the proportionality factor in (2.7) on the forcing correlation gives some robustness to characterizing the onset of non-Gaussianity by means of either $\textit {Re}_\lambda$ and $\textit {Re}_\gamma$, and we will see that our estimate for the critical $\textit {Re}_\lambda$ is consistent with that of Yakhot & Donzis (2017).

2.2. Fokker–Planck equation for the velocity gradient invariants

Equation (2.2) is associated with a FPE for the p.d.f. of the velocity gradient $f(\boldsymbol{\mathsf{A}}; t)$. In Cartesian components, the FPE for an ensemble of fluid particles sharing the same instantaneous configuration of the velocity gradient $\boldsymbol{\mathsf{A}}$ reads (Wilczek & Meneveau 2014),

(2.8)\begin{align} \frac{\partial f}{\partial t} &= \frac{\partial}{\partial A_{ij}} \left[\vphantom{\left(4\frac{\partial f}{\partial A_{ij}} - \frac{\partial f}{\partial A_{ji}} \right)} \textit{Re}_\gamma\left(A_{ik}A_{kj} + \left\langle H_{ij}(\boldsymbol{x},t)\middle|\boldsymbol{\mathsf{A}}\right\rangle\right)f - \big\langle\nabla^2 A_{ij}(\boldsymbol{x},t)\big|\boldsymbol{\mathsf{A}}\big\rangle\, f \right.\nonumber\\ &\quad + \left.\frac{1}{2}\sigma^2\left(4\frac{\partial f}{\partial A_{ij}} - \frac{\partial f}{\partial A_{ji}} \right)\right], \end{align}

where $\langle {\cdot } |\boldsymbol{\mathsf{A}} \rangle$ denotes the ensemble average conditional on the velocity gradient configuration $\boldsymbol{\mathsf{A}}$. The conditional averages in (2.8), namely the anisotropic part of the conditional pressure Hessian and the viscous Laplacian of the gradient, are unclosed terms that cannot be computed based only on the gradient at a single point (Meneveau 2011). Thus, the simplifications of going from the equation for the whole velocity gradient (2.2) (encoding the full space–time complexity of the velocity gradient field realizations) to an equation for the single-time/single-point statistics of the gradients (2.8), come with the cost of introducing unclosed terms.

In statistically isotropic flows, the velocity gradient p.d.f. $f$ is rotationally invariant, namely a function of only the five independent velocity gradient invariants (Itskov 2015). Statistical isotropy then allows us to employ tensor function representation theory (e.g. Rivlin & Ericksen 1955; Itskov 2015) to express the unclosed conditional averages in (2.8) as isotropic tensor functions of the velocity gradient (e.g. Pope 1975; Novikov 1993). The conditional averages are represented as combinations of basis tensors $\boldsymbol{\mathsf{B}}^n$ with coefficients $\gamma _n$ that depend upon the velocity gradient invariants $\mathcal {I}_k$. More specifically, the basis tensors are formed through the strain-rate tensor $\boldsymbol{\mathsf{S}}=(\boldsymbol{\mathsf{A}} + \boldsymbol{\mathsf{A}}^\top )/2$ and rotation-rate tensor $\boldsymbol{\mathsf{W}}=(\boldsymbol{\mathsf{A}} - \boldsymbol{\mathsf{A}}^\top )/2$, and they read

(2.9)\begin{equation} \left.\begin{gathered} \boldsymbol{\mathsf{B}}^1 = \boldsymbol{\mathsf{S}}, \quad \boldsymbol{\mathsf{B}}^3 = \widetilde{\boldsymbol{\mathsf{SS}}}, \quad \boldsymbol{\mathsf{B}}^5 = \boldsymbol{\mathsf{S}} \boldsymbol{\mathsf{W}} + \boldsymbol{\mathsf{W}} \boldsymbol{\mathsf{S}}, \quad \boldsymbol{\mathsf{B}}^7 = \widetilde{\boldsymbol{\mathsf{S}} \boldsymbol{\mathsf{WW}}} + \widetilde{\boldsymbol{\mathsf{WW}} \boldsymbol{\mathsf{S}}}, \\ \boldsymbol{\mathsf{B}}^2 = \boldsymbol{\mathsf{W}}, \quad \boldsymbol{\mathsf{B}}^4 = \boldsymbol{\mathsf{S}} \boldsymbol{\mathsf{W}} - \boldsymbol{\mathsf{W}} \boldsymbol{\mathsf{S}}, \quad \boldsymbol{\mathsf{B}}^6 = \widetilde{\boldsymbol{\mathsf{WW}}}, \quad \boldsymbol{\mathsf{B}}^8 = \widetilde{\boldsymbol{\mathsf{SS}} \boldsymbol{\mathsf{W}}} + \widetilde{\boldsymbol{\mathsf{W}} \boldsymbol{\mathsf{SS}}}, \end{gathered}\right\} \end{equation}

where the tilde indicates the traceless/anisotropic part of the tensor, and the standard matrix product is implied. It would be necessary to include two additional basis tensors in (2.9) to fix a possible degeneracy of the basis (Pennisi & Trovato 1987), but we ignore those zero-measure configurations. The independent invariants $\mathcal {I}_k$ formed through the velocity gradient read

(2.10)\begin{equation} \left.\begin{gathered} \mathcal{I}_1 = \textrm{Tr}\left(\boldsymbol{\mathsf{S}} \boldsymbol{\mathsf{S}}\right), \quad \mathcal{I}_3 = \textrm{Tr}\left(\boldsymbol{\mathsf{S}} \boldsymbol{\mathsf{S}} \boldsymbol{\mathsf{S}}\right), \quad \mathcal{I}_5 = \textrm{Tr}\left(\boldsymbol{\mathsf{S}} \boldsymbol{\mathsf{S}} \boldsymbol{\mathsf{W}} \boldsymbol{\mathsf{W}}\right),\\ \mathcal{I}_2 = \textrm{Tr}\left(\boldsymbol{\mathsf{W}} \boldsymbol{\mathsf{W}}\right), \quad \mathcal{I}_4 = \textrm{Tr}\left(\boldsymbol{\mathsf{S}} \boldsymbol{\mathsf{W}} \boldsymbol{\mathsf{W}}\right). \end{gathered}\right\} \end{equation}

A sixth invariant would be necessary to fix the handedness of the strain-rate eigenvector basis with respect to the vorticity. However, this sixth invariant is uniquely determined in terms of the other five only up to a sign (Lund & Novikov 1992) and we do not consider it as an independent variable.

With the above definitions of the basis tensors (2.9) and invariants (2.10), the drift term in the FPE (2.8) can be compactly written as

(2.11) \begin{equation} \textit{Re}_\gamma\left(\boldsymbol{\mathsf{AA}} + \left\langle \boldsymbol{\mathsf{H}}(\boldsymbol{x},t)\middle|\boldsymbol{\mathsf{A}}\right\rangle\right) - \big\langle\nabla^2 \boldsymbol{\mathsf{A}}(\boldsymbol{x},t)\big|\boldsymbol{\mathsf{A}}\big\rangle = \sum_{n=1}^{8} \gamma_n(\mathcal{I}) \boldsymbol{\mathsf{B}}^n. \end{equation}

The number of basis tensors (2.9) that are necessary to form a basis and to represent the most general drift (2.11), depends on the functional form of the polynomial coefficients $\gamma _n(\mathcal {I})$. If the coefficients are polynomials of the invariants, sixteen basis tensors are necessary to form a basis (Tian, Livescu & Chertkov 2021; Buaria & Sreenivasan 2023). However, if we relax the constraints on the coefficients, thus allowing them to be generic functions of the velocity gradient invariants (2.10), then we need just ten basis tensors to form a complete basis to represent any traceless tensor function (Pennisi & Trovato 1987). Furthermore, two of those ten basis tensors are linear combinations of the others, except when the strain-rate eigenvalues coincide and/or the vorticity is an eigenvector of the strain rate. Those configurations occur with zero probability, since during the dynamical evolution of the velocity gradients two of the eigenvalues may be very close at a certain time, but they will be pushed far apart by the dynamics at some later time. Equivalently, if we draw our random matrix $\boldsymbol{\mathsf{A}}$ from a realistic turbulent distribution, the probability of realizations featuring coincident strain-rate eigenvalues and/or vorticity parallel to some of the strain-rate eigenvectors is zero. We finally have the eight basis tensors (2.9), the same in number as the independent Cartesian components of a traceless matrix.

The general expression (2.11) allows us to formulate the FPE (2.8) in terms of the invariants by carrying out the contractions using symbolic calculus (Meurer et al. 2017). For notation simplicity, we denote the p.d.f. of the velocity gradient with $f$ regardless of its argument, $f(\mathcal {I}(\boldsymbol{\mathsf{A}})) \equiv f(\boldsymbol{\mathsf{A}})$ since $f$ is a probability density with respect to the traceless tensor $\boldsymbol{\mathsf{A}}$ and the invariants $\mathcal {I}_k$ are employed only for a simpler (isotropic) parametrization. Finally, the steady-state FPE (2.8) reduces to

(2.12)\begin{align} \gamma_n \phi^n f + M_{kl}Z^{ln} \frac{\partial \left(\gamma_n f\right)}{\partial \mathcal{I}_k} - \frac{\sigma^2}{2} M_{kl}\left( 4\phi^l-\phi'^l \right)\frac{\partial f}{\partial \mathcal{I}_k} = \frac{\sigma^2}{2} M_{kl}M_{pq}\left( 4Z^{lq}-Z'^{lq} \right)\frac{\partial^2 f}{\partial \mathcal{I}_k \partial \mathcal{I}_{p}} , \end{align}

with sum over repeated indices implied. The symmetric matrix $\boldsymbol{\mathsf{Z}}$ in (2.12) is the metric tensor $Z^{ln}(\mathcal {I}) = B^l_{ij}B^n_{ij}$, which in matrix form reads

(2.13) \begin{align} \boldsymbol{\mathsf{Z}} = {\left[\begin{matrix}\mathcal{I}_1 & 0 & \mathcal{I}_3 & 0 & 0 & \mathcal{I}_4 & 2 \mathcal{I}_5 & 0\\0 & - \mathcal{I}_2 & 0 & 0 & - 2 \mathcal{I}_4 & 0 & 0 & - 2 \mathcal{I}_5\\\mathcal{I}_3 & 0 & \frac{\mathcal{I}_1^{2}}{6} & 0 & 0 & - \frac{\mathcal{I}_1 \mathcal{I}_2}{3} + \mathcal{I}_5 & \frac{\mathcal{I}_1 \mathcal{I}_4}{3} + \frac{2 \mathcal{I}_2 \mathcal{I}_3}{3} & 0\\0 & 0 & 0 & \mathcal{I}_1 \mathcal{I}_2 - 6 \mathcal{I}_5 & 0 & 0 & 0 & 0\\0 & - 2 \mathcal{I}_4 & 0 & 0 & - \mathcal{I}_1 \mathcal{I}_2 + 2 \mathcal{I}_5 & 0 & 0 & - \mathcal{I}_1 \mathcal{I}_4 - \frac{\mathcal{I}_2 \mathcal{I}_3}{3}\\\mathcal{I}_4 & 0 & - \frac{\mathcal{I}_1 \mathcal{I}_2}{3} + \mathcal{I}_5 & 0 & 0 & \frac{\mathcal{I}_2^{2}}{6} & \frac{\mathcal{I}_2 \mathcal{I}_4}{3} & 0\\2 \mathcal{I}_5 & 0 & \frac{\mathcal{I}_1 \mathcal{I}_4}{3} + \frac{2 \mathcal{I}_2 \mathcal{I}_3}{3} & 0 & 0 & \frac{\mathcal{I}_2 \mathcal{I}_4}{3} & - \frac{\mathcal{I}_1 \mathcal{I}_2^{2}}{2} + 3 \mathcal{I}_2 \mathcal{I}_5 + \frac{2 \mathcal{I}_4^{2}}{3} & 0\\0 & - 2 \mathcal{I}_5 & 0 & 0 & - \mathcal{I}_1 \mathcal{I}_4 - \frac{\mathcal{I}_2 \mathcal{I}_3}{3} & 0 & 0 & \frac{\mathcal{I}_1^{2} \mathcal{I}_2}{2} - 3 \mathcal{I}_1 \mathcal{I}_5 + \frac{2 \mathcal{I}_3 \mathcal{I}_4}{3}\end{matrix}\right]}, \end{align}

while $Z'^{ln} = B^l_{ij}B^n_{ji}$ denotes a modified metric tensor. The symbols $\phi ^n(\mathcal {I}) = \partial B^n_{ij}/\partial A_{ij}$ indicate the divergence of the basis tensors

(2.14)\begin{equation} \boldsymbol{\phi} = \left[5 \quad 3 \quad 0 \quad 0 \quad 0 \quad 0 \quad \frac{10 \mathcal{I}_2}{3} \quad 2 \mathcal{I}_1\right], \end{equation}

while the divergence of the transposed basis tensors is $\phi '^n = \partial B^n_{ji}/\partial A_{ij}$. The components of the derivatives of the invariants are collected in the matrix $M_{kl}=(\partial \mathcal {I}_k/\partial A_{ij})B^n_{ij}Z^{-1}_{nl}$ (with $\boldsymbol{\mathsf{Z}}^{-1}$ denoting the inverse of the metric tensor (2.13)), which in matrix form reads

(2.15)\begin{equation} \boldsymbol{\mathsf{M}} = \left[\begin{matrix}2 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & -2 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 3 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & -1 & 1 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 1 & -1\end{matrix}\right]. \end{equation}

The quantities (2.13), (2.14), (2.15) characterizing the tensor basis (2.9) can all be derived from the Christoffel symbols computed in Carbone & Wilczek (2022).

Equation (2.12) is a second-order partial differential equation for the function $f(\mathcal {I})$ of the five variables $\mathcal {I}_k$. The model coefficients $\gamma _n(\mathcal {I})$ determine the properties and complexity of the FPE, and we will now compute the model coefficients $\gamma _n$ for low-Reynolds-number flows.

3. Determining the model coefficients

We introduce the two hypotheses necessary to obtain the presented model: the coefficients $\gamma _n$ in (2.11) are constant and the expansion (2.11) is truncated to basis tensors up to degree two in the velocity gradient. While these assumptions are exact for the conditional pressure Hessian at first order in Reynolds number, they introduce modelling approximations for the viscous stresses. We will test the validity of the assumptions a posteriori, by comparing the model coefficients $\gamma _n$ with the same coefficients computed directly from DNS of low-Reynolds-number flows (see figure 7).

3.1. Zeroth-order conditional velocity gradient Laplacian

At order zero in $\textit {Re}_\gamma$, the equation governing the velocity gradient dynamics (2.2) is linear, and the resulting flow has Gaussian statistics. For a multipoint Gaussian random field, the conditional Laplacian reduces to a linear damping (Wilczek & Meneveau 2014), and in the non-dimensional variables (2.5a–d) we have

(3.1)\begin{equation} \left\langle \nabla^2 \boldsymbol{\mathsf{A}}( \boldsymbol{x},t)\middle|\boldsymbol{\mathsf{A}}\right\rangle ={-}\boldsymbol{\mathsf{A}} + \textit{O}\left(\textit{Re}_\gamma\right). \end{equation}

The conditional velocity gradient Laplacian for a Gaussian random field contributes to the model coefficients (2.11) at order zero in Reynolds number, while the higher-order viscous corrections require modelling.

3.2. Zeroth-order conditional anisotropic pressure Hessian

For a Gaussian random field, that is at order zero in Reynolds number, the conditional traceless/anisotropic pressure Hessian takes the simple form (Wilczek & Meneveau 2014)

(3.2)\begin{equation} \left\langle \tilde{\boldsymbol{\mathsf{H}}}(\boldsymbol{x},t)\middle|\boldsymbol{\mathsf{A}}\right\rangle ={-}\tfrac{2}{7}\widetilde{\boldsymbol{\mathsf{SS}}} - \tfrac{2}{5}\widetilde{\boldsymbol{\mathsf{WW}}} + h_4 \left(\boldsymbol{\mathsf{SW}} - \boldsymbol{\mathsf{WS}}\right) +\textit{O}(\textit{Re}_\gamma), \end{equation}

while the local/isotropic part of the Hessian is specified by the incompressibility condition $\textrm {Tr}(\boldsymbol{\mathsf{H}})=-\textrm {Tr}(\boldsymbol{\mathsf{AA}})$. The conditional anisotropic pressure Hessian for a Gaussian random field contributes to the model coefficients (2.11) at order one in Reynolds number. Therefore, we know from Wilczek & Meneveau (2014) the exact first-order correction to the gradient dynamics at small $\textit {Re}_\gamma$ due to the pressure Hessian.

The coefficient $h_4$ weighing $\boldsymbol{\mathsf{B}}^4 = \boldsymbol{\mathsf{SW}} - \boldsymbol{\mathsf{WS}}$ in (3.2) depends on the structure of the Gaussian flow through the correlation function (Wilczek & Meneveau 2014; Johnson & Meneveau 2016), but previous works have shown that it does not contribute to the single-point statistics of the velocity gradient. This can be seen geometrically since $\boldsymbol{\mathsf{B}}^4$ rotates the strain-rate eigenframe while leaving unchanged the vorticity orientation with respect to the eigenframe itself (Carbone, Iovieno & Bragg 2020). It can also be seen from the FPE for the gradient p.d.f. (2.12) since $\boldsymbol{\mathsf{B}}^4$ contracts to zero with all the other basis tensors (Leppin & Wilczek 2020), as from the fourth row/column of the metric tensor (2.13). Therefore, we can ignore the contributions from $\boldsymbol{\mathsf{B}}^4$ as long as we are concerned with single-time/single-point statistics.

3.3. Tensor representation of the conditional averages and resulting velocity gradient model

We model the higher-order contributions from the unclosed terms by employing the general expression (2.11), truncated at degree two in the velocity gradient. We consider the basis tensors (2.9) from $\boldsymbol{\mathsf{B}}^1$ to $\boldsymbol{\mathsf{B}}^6$, by keeping in mind that $\boldsymbol{\mathsf{B}}^4$ can be ignored as long as we focus on single-point statistics. Also, we assume that the coefficients $\gamma _n$ in (2.11) are constant. These hypotheses, together with the exact expression of the zeroth-order conditional averages (3.1), (3.2), yield the following representation of the drift term (2.11):

(3.3)\begin{align} &\textit{Re}_\gamma\left(\boldsymbol{\mathsf{A}}\boldsymbol{\mathsf{A}} + \left\langle \boldsymbol{\mathsf{H}}(\boldsymbol{x},t)\middle|\boldsymbol{\mathsf{A}}\right\rangle \right) -\left\langle\nabla^2 \boldsymbol{\mathsf{A}}(\boldsymbol{x},t)\middle|\boldsymbol{\mathsf{A}}\right\rangle= \boldsymbol{\mathsf{A}} -\textit{Re}_\gamma^2\left[\delta_1\boldsymbol{\mathsf{S}} + \delta_2\boldsymbol{\mathsf{W}}\right] \nonumber\\ &\quad +\textit{Re}_\gamma\left[\left(\tfrac{5}{7}-\delta_3\right)\widetilde{\boldsymbol{\mathsf{SS}}} + \gamma_4\left(\boldsymbol{\mathsf{SW}}-\boldsymbol{\mathsf{WS}}\right) + \left(1-\delta_5\right)\left(\boldsymbol{\mathsf{SW}}+\boldsymbol{\mathsf{WS}}\right) + \left(\tfrac{3}{5}-\delta_6\right)\widetilde{\boldsymbol{\mathsf{WW}}}\right]. \end{align}

Here the constant coefficients $\delta _i$ are of order one and are to be determined. We include second-order terms in $\textit {Re}_\gamma$ to keep the variance of the gradients constant for all Reynolds numbers, as is the case for the stochastically driven Navier–Stokes equations (2.1). The powers in the Reynolds numbers in (3.3) have been chosen so that the model equations remain unchanged under the transformation $\bar {t}\to -\bar {t}$ and $\bar {\nu }\to -\bar {\nu }$, as is the case for the Navier–Stokes equations (2.1). We will test a posteriori the trend of the model coefficients in terms of the Reynolds number against the DNS data (see figure 7).

The constant model parameters $\delta _1,\delta _2,\delta _3,\delta _5,\delta _6$ remain to be determined. To do this, we use the two Betchov homogeneity constraints (Betchov 1956), the perturbation theory at a small Reynolds number for the full Navier–Stokes equations (Wyld 1961), together with the constant dissipation rate imposed by the stochastic forcing (Novikov 1965). This results in constraints on the average of the velocity gradient invariants

(3.4a)$$\begin{gather} \left\langle \mathcal{I}_1 \right\rangle = \tfrac{1}{2}, \quad \left\langle \mathcal{I}_1+\mathcal{I}_2 \right\rangle = 0, \quad \left\langle \mathcal{I}_3+3\mathcal{I}_4 \right\rangle = 0, \end{gather}$$

(3.4b)$$\begin{gather}\left\langle \mathcal{I}_3 \right\rangle = S_3\textit{Re}_\gamma, \quad \left\langle \mathcal{I}_5 \right\rangle ={-}\tfrac{1}{12} + X_5\textit{Re}_\gamma^2, \end{gather}$$

where $S_3$ and $X_5$ are constant parameters. The first relation in (3.4a) follows from the constant variance of the gradients imposed by the stochastic forcing. It implies that the Kolmogorov time scale $\tau _\eta =1/\sqrt {2\left \langle \mathcal {I}_1 \right \rangle }$ is one in the non-dimensional variables (2.5a–d). The other relations in (3.4a) are the two independent Betchov homogeneity constraints that can be formulated using solely the velocity gradient (Carbone & Wilczek 2022). The homogeneity relations (3.4a) have already been employed to reduce the number of parameters in numerical simulations of velocity gradient models at high Reynolds numbers (Girimaji & Pope 1990; Johnson & Meneveau 2016; Leppin & Wilczek 2020), and here we can impose those constraints analytically, as we will see below. Relations (3.4b) follow from a low-Reynolds-number expansion of the Navier–Stokes equations (Wyld 1961). The quantities $S_3$ and $X_5$ represent, respectively, the rate of change of the third- and fourth-order moments of the velocity gradient with the Reynolds number, starting from a Gaussian zeroth-order configuration. These coefficients depend on the forcing correlation, and we determine them in Appendix B.

In the constraints (3.4) the brackets indicate the ensemble average, that is, for a generic function of the velocity gradient $\varphi (\boldsymbol{\mathsf{A}})$,

(3.5)\begin{equation} \left\langle \varphi(\boldsymbol{\mathsf{A}}) \right\rangle = \int\textrm{d}\boldsymbol{\mathsf{A}} \ f\left(\mathcal{I}(\boldsymbol{\mathsf{A}});\delta_i\right)\varphi(\boldsymbol{\mathsf{A}}), \end{equation}

where $f(\mathcal {I};\delta _i)$ is the p.d.f. of the velocity gradient governed by the FPE (2.12), with the drift term (2.11). The FPE (2.12) admits perturbative analytic solutions at small Reynolds numbers, as described in further detail in the next section. This analytic solution $f(\mathcal {I};\delta _i)$ depends parametrically on the model coefficients $\delta _i$, and it allows us to compute the ensemble averages in (3.4) analytically. Therefore, the five constraints (3.4) constitute a linear system of five equations for the five model parameters $\delta _1,\delta _2,\delta _3,\delta _5,\delta _6$. The solution of this linear system is

(3.6)\begin{equation} \left.\begin{gathered} \delta_1 ={-} \frac{1824 S_{3}^{2}}{35} - \frac{144 X_{5} }{7},\quad \delta_2 = \frac{96 S_{3}^{2} }{7} + \frac{240 X_{5}}{7}, \\ \delta_3 = \frac{5}{7} + \frac{120 S_{3} }{7}, \quad \delta_4 = \gamma_4 -h_4, \\ \delta_5 = 1 - \zeta_5 - \frac{72 S_{3} }{7} - \frac{180 X_{5} }{7 S_{3}}, \quad \delta_6 = \frac{3}{5} - \frac{6 \zeta_5}{5} - \frac{936 S_{3} }{35} - \frac{216 X_{5} }{7 S_{3}}. \end{gathered}\right\}\end{equation}

Additionally, the noise variance $\sigma ^2=1/15$ is fixed by the constraint $\tau _\eta =1$ at $\textit {Re}_\gamma =0$ (at which the FPE (2.12) with the drift (3.3) is an Ornstein–Uhlenbeck process). We assume that the noise variance $\sigma ^2$ is independent of the Reynolds number.

We have now determined all the single-time/single-point model coefficients, $\delta _i$ and $\sigma$. The model features the exact first-order contribution in $\textit {Re}_\gamma$ from the pressure Hessian, while the modelling hypotheses concern the representation of the higher-order corrections. The resulting model for the velocity gradient single-time statistics does not feature any free parameter, not requiring any parameter scan to match the DNS results. Still, there are two free gauge parameters, namely $\gamma _4$ and $\zeta _5$, which do not affect single-time statistics, but only multi-time correlations. This gauge stems from the fact that the Gaussian and isotropic p.d.f. with unity Kolmogorov time scale (e.g. Wilczek & Meneveau 2014)

(3.7)\begin{equation} f_0 = \frac{225 \sqrt{5}}{{\rm \pi}^{4}} \exp\left[- \left( 4\delta_{ip}\delta_{jq}+ \delta_{iq}\delta_{jp}\right)A_{ij}A_{pq}\right] \end{equation}

solves the nonlinear steady-state FPE

(3.8)\begin{equation} \frac{\partial}{\partial A_{ij}} \left[ \left(A_{ij} + \textit{Re}_\gamma \zeta_5\left( B^5_{ij}+\frac{6}{5}B^6_{ij}\right) + \textit{Re}_\gamma \gamma_4 B^4_{ij}\right)f_0 + \frac{1}{30}\left(4\frac{\partial f_0}{\partial A_{ij}} - \frac{\partial f_0}{\partial A_{ji}} \right)\right] = 0 \end{equation}

for all $\gamma _4$ and $\zeta _5$. A particular case of this gauge for the zeroth-order Gaussian solution has been observed by Leppin & Wilczek (2020). We will determine the gauge parameters $\gamma _4$ and $\zeta _5$ in § 6, where we focus on multi-time statistics, with the aid of DNS data.

4. Analytic approximation of the velocity gradient p.d.f. at small Reynolds numbers

The FPE (2.12) with the drift term specified by (3.3) admits perturbative solutions in the Reynolds number. The solution is expanded up to second order,

(4.1)\begin{equation} f(\mathcal{I}) \sim f_0(\mathcal{I}) + \textit{Re}_\gamma f_1(\mathcal{I}) + \textit{Re}_\gamma^2 f_2(\mathcal{I}), \end{equation}

and we solve for $f_i$ at all orders, in the form of a polynomial of the invariants times the zeroth-order Gaussian solution. In particular, plugging the expansion (4.1) into the FPE (2.12), comparing terms of the same order in Reynolds number, and imposing the average constraints (3.4) yields the following asymptotic solution of the FPE at small $\textit {Re}_\gamma$:

(4.2a)\begin{gather} f_0 = \frac{225 \sqrt{5}}{{\rm \pi}^{4}} \exp({- 5 \mathcal{I}_1 + 3 \mathcal{I}_2}),\end{gather}

(4.2b)\begin{gather} f_1 = \frac{3600 \sqrt{5} S_{3} \left(25 \mathcal{I}_3 - 21 \mathcal{I}_4\right)}{7 {\rm \pi}^{4}} \exp({- 5 \mathcal{I}_1 + 3 \mathcal{I}_2}), \end{gather}

(4.2c)\begin{align} f_2 &= \frac{720 \sqrt{5}}{49 {\rm \pi}^{4}} \left(- 16\,320 S_{3}^{2} \mathcal{I}_1 \mathcal{I}_2 - 6860 S_{3}^{2} \mathcal{I}_1 - 1344 S_{3}^{2} \mathcal{I}_2^{2} - 140 S_{3}^{2} \mathcal{I}_2 + 50\,000 S_{3}^{2} \mathcal{I}_3^{2} \right. \nonumber\\ &\quad + 84\,000 S_{3}^{2} \mathcal{I}_3 \mathcal{I}_4 + 35\,280 S_{3}^{2} \mathcal{I}_4^{2} + 42\,240 S_{3}^{2} \mathcal{I}_5 + 2240 S_{3}^{2} - 22\,950 X_{5} \mathcal{I}_1 \mathcal{I}_2 - 1575 X_{5} \mathcal{I}_1 \nonumber\\ &\quad \left .+ 1890 X_{5} \mathcal{I}_2^{2} - 1575 X_{5} \mathcal{I}_2 + 59400 X_{5} \mathcal{I}_5\right) \exp({- 5 \mathcal{I}_1 + 3 \mathcal{I}_2}). \end{align}

We remark that while the solution (4.1) is only asymptotic in $\textit {Re}_\gamma$, the terms $f_i$ in (4.2) solve exactly, for all $\mathcal {I}_k$, the expanded partial differential equation (2.12) at each order in the Reynolds number. Therefore, while we deal with an asymptotic expansion at small Reynolds number, there is no explicit assumption on the magnitude of the velocity gradients. To assess the asymptotic solution (4.2) by symbolic computation, we rewrote the terms in (4.2) as functions of the Cartesian components of the velocity gradient, inserted the resulting expression into the Cartesian FPE (2.8), and checked that the remainder is zero up to second order in Reynolds number.

The main advantage of the simple expansion (4.2) is that it gives full analytic access to the relevant moments of the p.d.f. This allows us to investigate analytically the onset of non-Gaussianity at small Reynolds number. A shortcoming of this expansion is that the approximate p.d.f. (4.2) may not be positive for all $\mathcal {I}_k$ and $\textit {Re}_\gamma$. The approximation to the velocity gradient p.d.f. (4.2) changes sign when the polynomial prefactor (which multiplies the Gaussian exponential part) vanishes. Setting the polynomial prefactor to zero is equivalent to solving a quadratic equation for $\textit {Re}_\gamma$. The roots of that algebraic equation, as functions of the velocity gradient invariants, show that the p.d.f. can become negative at low $\textit {Re}_\gamma$ only for large values of the invariants. At those large values of the velocity gradient invariants the Gaussian exponential part $\exp (- 5 \mathcal {I}_1 + 3 \mathcal {I}_2)$ has already strongly decayed. The comparison with the numerical results in § 5 will show that this positivity issue is indeed negligible in the explored range of Reynolds numbers.

An alternative way to obtain an asymptotic solution of the FPE (2.12) consists of the effective action method, proposed in Martin et al. (1973) for stochastic differential equations. This method yields solutions of the form $\exp (-\mathcal {S}(\mathcal {I}))$, where $\mathcal {S}$ is the effective action, featuring higher-order polynomials in $\mathcal {I}_k$ and $\textit {Re}_\gamma$. Such a p.d.f. does not give analytic access to the moments, which are usually computed numerically via Monte-Carlo sampling (Moriconi et al. 2014). The effective action involves renormalized noise variance and model coefficients (Apolinário et al. 2019) aiming to improve the accuracy and range of validity of the analytic predictions. Our model targets $\textit {Re}_\gamma$ very small, and the equations themselves are only valid for low Reynolds numbers. Therefore, in this set-up, the simple expansion yielding the solution (4.2) seems appropriate, and we leave more advanced approaches for future work.

4.1. Volume elements and the moments of the velocity gradient

The solution (4.2) allows us to compute ensemble averages analytically. The ensemble average (3.5) is conveniently computed in the strain-rate eigenframe. We express all the invariants $\mathcal {I}_k$ in terms of the strain-rate eigenvalues $\lambda _i$ and vorticity principal components $\omega _i=\boldsymbol {v}_i\boldsymbol {\cdot }\boldsymbol {\omega }$, where $\boldsymbol {v}_i$ are the strain-rate eigenvectors. Additionally, we use coordinates in the strain-rate eigenframe based on the vorticity magnitude $\omega =\|\boldsymbol {\omega }\|$ and the alignments between the vorticity vector and strain-rate eigenvectors, $\hat {\omega }_i\equiv \omega _i/\omega$. This procedure transforms the integral (3.5) into

(4.3)\begin{equation} \left\langle \varphi(\mathcal{I}) \right\rangle =\int\textrm{d}\boldsymbol{\mathsf{S}}\int\textrm{d}\omega_1\,\textrm{d}\omega_2\,\textrm{d}\omega_3 \ f(\mathcal{I}) \varphi(\mathcal{I}), \end{equation}

with the integrations in the strain-rate eigenframe specified by

(4.4a)$$\begin{gather} \int\textrm{d}\boldsymbol{\mathsf{S}} = \int_{0}^{\infty}\textrm{d}\lambda_2\int_{\lambda_2}^{+\infty}\textrm{d}\lambda_1 \,J_S(\lambda) + \int_{-\infty}^{0}\textrm{d}\lambda_2\int_{{-}2\lambda_2}^{+\infty}\textrm{d}\lambda_1 \,J_S(\lambda), \end{gather}$$

(4.4b)$$\begin{gather}\int\textrm{d}\omega_1\,\textrm{d}\omega_2\,\textrm{d}\omega_3 = \int_{0}^{+\infty}\textrm{d}\omega^2\int_0^1\textrm{d}\hat{\omega}_1^2\int_0^1\textrm{d} \hat{\omega}_2^2 \,J_\omega(\omega,\hat{\omega}), \end{gather}$$

and where the invariants are also expressed as functions of the strain-rate variables,

(4.5)\begin{equation} \left.\begin{gathered} \mathcal{I}_1 = \sum_i \lambda_i^2, \quad \mathcal{I}_2 ={-}\frac{\omega^2}{2}, \\ \mathcal{I}_3 = \sum_i \lambda_i^3, \quad \mathcal{I}_4 = \frac{\omega^2}{4}\sum_i \lambda_i\hat{\omega}_i^2, \quad \mathcal{I}_5 = \frac{\omega^2}{4}\sum_i \lambda_i^2\left(\hat{\omega}_i^2-1\right). \end{gathered}\right\}\end{equation}

Due to the statistical isotropy, rotations of the strain-rate eigenframe can be integrated out. This reduces the dimensionality of the ensemble average integral from eight (i.e. the independent components of the traceless $\boldsymbol{\mathsf{A}}$ in (3.5)) to five (i.e. the independent strain-rate eigenvalues and the vorticity principal components in (4.4)). However, integrating out rotations also comes at the cost of introducing volume elements (Livan et al. 2018), namely $J_S$ and $J_\omega$ in (4.4a). The volume element $J_S$ associated with the transformation from a standard Cartesian reference frame to the strain-rate eigenframe consists of the Wigner repulsion term (Wigner 1955)

(4.6)\begin{equation} J_S = 2{\rm \pi}^2\left|(\lambda_1-\lambda_2)(\lambda_2-\lambda_3)(\lambda_1-\lambda_3)\right| = \sqrt{2}{\rm \pi}^2\sqrt{\mathcal{I}_1^3-6\mathcal{I}_3^2}. \end{equation}

The absolute value in (4.6) drops when employing the ordered strain-rate eigenvalues $\lambda _1>\lambda _2>\lambda _3$ as integration variables. Due to incompressibility, $\lambda _1+\lambda _2+\lambda _3=0$, and we have two possible configurations, $\lambda _1>\lambda _2>0$ or $\lambda _1>-2\lambda _2>0$, which specify the integration bounds in (4.4a). Finally, going from Cartesian coordinates to the vorticity magnitude/orientation in the strain-rate eigenframe introduces the volume element

(4.7)\begin{equation} J_\omega = \frac{1}{8}\sqrt{\frac{\omega^2}{\hat{\omega}^2_1\hat{\omega}^2_2\hat{\omega}^2_3}} , \end{equation}

where $\hat {\omega }_3^2 = 1-\hat {\omega }_1^2-\hat {\omega }_2^2$, and we can consider only positive values of the vorticity principal components $\hat {\omega }_i$ since the invariants $\mathcal {I}_k$ (4.5) are even functions of the vorticity.

5. Comparison of single-point/single-time statistics

We compare the single-time/single-point velocity gradient statistics resulting from our model and from DNS at low Reynolds numbers. We point out several qualitative changes in the dynamics and statistical geometry of the gradient as the Reynolds number increases until a transition to turbulence. The DNS set-up is detailed in Appendix A, while all the theoretical predictions follow by integrating out variables from the asymptotic solution (4.1), making use of the expressions for the ensemble average (4.4). The DNS and model parameters are in one-to-one correspondence. In particular, the zeroth-order conditional viscous damping and the correlations of the forcing are the same in the model and in the DNS, as described in Appendices A and B.

5.1. Moments of the velocity gradient invariants

Increasing the Reynolds number starting from zero leads to the onset of the skewness, intermittency and preferential alignments in the gradient statistics. We analyse these three features separately by looking at the strain-rate and vorticity statistics from the strain-rate eigenframe viewpoint (Dresselhaus & Tabor 1992; Tom, Carbone & Bragg 2021).

Figure 1(a) shows the normalized moments of the strain-rate longitudinal components as a function of the Reynolds number $\textit {Re}_\gamma$. The moments of the invariants are related to the moment of the longitudinal strain-rate component through (Betchov 1956; Davidson 2015)

(5.1a–c)\begin{equation} \left\langle A_{11}^2 \right\rangle = \tfrac{2}{15}\left\langle \mathcal{I}_1 \right\rangle, \quad \left\langle A_{11}^3 \right\rangle = \tfrac{8}{105}\left\langle \mathcal{I}_3 \right\rangle, \quad \left\langle A_{11}^4 \right\rangle = \tfrac{4}{105}\left\langle \mathcal{I}_1^2 \right\rangle. \end{equation}

The skewness of the strain rate is quantified via the third-order moment of the strain rate $\left \langle \mathcal {I}_3 \right \rangle$, while the even and higher-order moments, e.g. $\left \langle \mathcal {I}_1^2 \right \rangle$, quantify the intermittency. After the rapid initial increase, the skewness of the longitudinal strain-rate component approaches its typical high-Reynolds-number value of $-0.5/-0.6$. At large Reynolds number, this skewness is predicted to be constant by a renormalization approach (Yakhot & Orszag 1986), while it weakly increases with the Reynolds number with an exponent of order $0.1$ in numerical simulations (Ishihara, Gotoh & Kaneda 2009). At $\textit {Re}_\gamma =\textit {O}(1)$, the strain-rate statistics already display a remarkable skewness, while the intermittency is negligible. This is because the third-order moments, e.g. $\left \langle \mathcal {I}_3 \right \rangle$, grow linearly with the Reynolds number, while even-order moments, e.g. $\left \langle \mathcal {I}_1^2 \right \rangle$, grow quadratically. From this we infer that the skewness is a dominant feature in low-Reynolds-number flows, while the intermittency is negligible.

Figure 1. Onset of skewness, intermittency and alignments at low Reynolds number, in terms of the expansion parameter $\textit {Re}_\gamma$. (a) Normalized moments of the longitudinal velocity gradient component $A_{11}$. Note that we plot $\langle (-A_{11})^n\rangle$ since the odd-order moments are negative. (b) Alignments between the vorticity and the strain-rate eigenvectors. Solid coloured lines are from DNS, while black dashed lines are the analytic predictions from our low-Reynolds-number model. The inset shows the deviation of the alignments from the Gaussian configuration, $|3\left \langle \hat {\omega }_i^2 \right \rangle -1|$, as a function of $\textit {Re}_\gamma$, in a log-log scale.

Figure 1(b) shows a striking qualitative change in the statistical geometry of the velocity gradient at low Reynolds numbers: the switching in the preferential alignments between the strain rate and the vorticity. The alignments are characterized through the normalized vorticity components in the strain-rate eigenframe $\hat {\omega }_i=\boldsymbol {v}_i\boldsymbol {\cdot }\boldsymbol {\omega }/\|\boldsymbol {\omega }\|$ (where $\boldsymbol {v}_i$ are the strain-rate eigenvectors), ordered based on the corresponding strain-rate eigenvalue $\lambda _i$, with $\lambda _1\geq \lambda _2\geq \lambda _3$. The vorticity aligns with the most extensional strain-rate eigenvector at low Reynolds number, and only at $\textit {Re}_\gamma \gtrapprox 5$ it aligns with the intermediate eigenvector. At $\textit {Re}_\gamma <5$, the alignments are as if the vorticity was a material line in the fluid flow subject to a persistent strain. Then, at $\textit {Re}_\gamma \simeq 5$, the transition to turbulence begins: the somewhat counter-intuitive alignment between the intermediate strain-rate eigenvector and the vorticity establishes, as in fully developed turbulence (Ashurst et al. 1987). The onset of this peculiar alignment at a higher Reynolds number is consistent with the local nonlinearities in the velocity gradient dynamics becoming more relevant. For example, the restricted Euler model (Vieillefosse 1982), in which the gradient is driven only by the local/anisotropic part of the nonlinear term, also displays alignment between the vorticity and the intermediate strain-rate eigenvector while approaching a finite-time singularity (Novikov 1990; Cantwell 1992).

By integrating out the strain-rate eigenvalues and vorticity magnitude from the asymptotic solution of the FPE (4.1), with the ensemble average expressed in the form of (4.4), we can get an analytic approximation of the vorticity principal components as a function of the Reynolds number

(5.2)\begin{equation} \left.\begin{gathered} \left\langle \hat{\omega}_1^2 \right\rangle = \frac{1}{3} - \frac{6 \sqrt{30} S_{3} \textit{Re}_\gamma}{25 \sqrt{\rm \pi}} + \frac{4224 S_{3}^{2} \textit{Re}_\gamma^{2}}{1225} + \frac{1188 X_{5} \textit{Re}_\gamma^{2}}{245}, \\ \left\langle \hat{\omega}_2^2 \right\rangle = \frac{1}{3} - \frac{8448 S_{3}^{2} \textit{Re}_\gamma^{2}}{1225} - \frac{2376 X_{5} \textit{Re}_\gamma^{2}}{245}, \\ \left\langle \hat{\omega}_3^2 \right\rangle = \frac{1}{3} + \frac{6 \sqrt{30} S_{3} \textit{Re}_\gamma}{25 \sqrt{\rm \pi}} + \frac{4224 S_{3}^{2} \textit{Re}_\gamma^{2}}{1225} + \frac{1188 X_{5} \textit{Re}_\gamma^{2}}{245}. \end{gathered}\right\}\end{equation}

Both $\left \langle \hat {\omega }_1^2 \right \rangle$ and $\left \langle \hat {\omega }_2^2 \right \rangle$ start from a random uniform value $\left \langle \hat {\omega }_i^2 \right \rangle =1/3$ at zero Reynolds number. Then $\left \langle \hat {\omega }_1^2 \right \rangle$ grows linearly with $\textit {Re}_\gamma$ while $\left \langle \hat {\omega }_2^2 \right \rangle$ grows only quadratically. The preferential alignments are closely related to the onset of the strain-rate skewness. Indeed, the fact that $S_3<0$ (analytically derived through the Wyld expansion in Appendix B) not only implies an average direct energy cascade, but it also implies that the vorticity preferentially aligns with the most extensional direction at low Reynolds number, as seen from (5.2). However, $\left \langle \hat {\omega }_1^2 \right \rangle$ has a strong negative second-order contribution in $\textit {Re}_\gamma$ and eventually, $\left \langle \hat {\omega }_2^2 \right \rangle$, which has a positive growth rate, takes over. The analytic approximation (5.2) suggests that the growth of $\left \langle \hat {\omega }_2^2 \right \rangle$ with $\textit {Re}_\gamma$ can only take place if $X_5$ is negative, and its magnitude is large enough compared with $S_3^2$. Furthermore, it is known that the vorticity avoids the compressional direction in turbulence (Ashurst et al. 1987). This lack of alignment establishes already at very low Reynolds numbers, and it strongly depends on the fact that $S_3$ and $X_5$ are both negative. Our model captures this lack of alignment since $\left \langle \hat {\omega }_3^2 \right \rangle$ in (5.2) has negative first- and second-order contributions, so it quickly decreases with the Reynolds number.

The comparison between the DNS and model predictions for the strain-rate moments and strain rate–vorticity alignments in figure 1 shows that our low-Reynolds-number model is quantitatively accurate up to $\textit {Re}_\gamma \simeq 1$, in qualitative agreement with the DNS until $\textit {Re}_\gamma \simeq 5$ and then breaks down, as expected for a perturbative weak-coupling approach. Interestingly, the model can predict the switching of the strain rate–vorticity alignments, that is essentially encoded in the first- and second-order corrections in $\textit {Re}_\gamma$ to the statistical alignments. However, the analytically predicted Reynolds number at which the flipping of the preferential alignments takes place is slightly lower than the actual one. This is because the perturbation parameter is already relatively large in the transition regime, $\textit {Re}_\gamma \simeq 5$.

The results highlight many qualitative changes as the Reynolds number increases, consistent with the observations in recent works (Yakhot & Donzis 2017; Gotoh & Yang 2022; Khurshid et al. 2023). These works have identified a transition to turbulence based on the sudden growth of the even velocity gradient moments and the onset of a power law as their growth begins, around $\textit {Re}_\lambda \simeq 9$. The transition is usually identified by means of the Reynolds number based on the Taylor microscale $\textit {Re}_\lambda$, while we localize the transition in terms of $\textit {Re}_\gamma$, that is the expansion parameter in our analysis. The two Reynolds numbers are proportional close to the Gaussian state, $\textit {Re}_\lambda \simeq 1.5\textit {Re}_\gamma$ with moderate dependence on the forcing details (see Appendix B), while they are non-trivially related at higher Reynolds numbers, their relation depending on the spatial correlations of the forcing. According to the estimate (2.7), the transition Reynolds number $\textit {Re}_\gamma \simeq 5$ corresponds to $\textit {Re}_\lambda \simeq 7.5$, that is consistent with the transition Reynolds number estimated via the analysis of the velocity gradient scaling exponents (Yakhot & Donzis 2017). Summarizing, when $\textit {Re}_\lambda ={O}(10)$ (more quantitatively $\textit {Re}_\lambda$ between 5 and 10), the scaling exponents of the velocity gradient moments as functions of the Reynolds number change and, as an independent indicator, the alignments between the strain rate and vorticity flip towards a turbulence-like configuration.

The transition from a Gaussian state to a turbulence-like configuration is smooth, as shown by our model and DNS results. This is in agreement with the results in the literature (e.g. Khurshid et al. 2023), which do not show sharp changes of the velocity gradient scaling exponents in the vicinity of the critical Reynolds number. Analogously, by taking the strain rate–vorticity alignments as an indicator of the velocity gradient state, we observe the flipping of the preferential alignments at a critical $\textit {Re}_\gamma$, but the transition towards a turbulent configuration is smooth and smeared in the vicinity of that critical Reynolds number. Based on our results, we interpret the transition to turbulence reported in the literature (e.g. Yakhot & Donzis 2017), not as a sharp change in the scaling exponents of the velocity gradient moments in terms of $\textit {Re}_\lambda$, but rather as an approximation to a smooth transition between two power-law trends with different exponents.

5.2. Cartesian velocity gradient components

The changes in the moments of the velocity gradient invariants reflect into the onset of non-Gaussianity in the statistics of the velocity gradient Cartesian components. In particular, the onset of skewed strain-rate statistics already at a very low Reynolds number results in the non-Gaussianity of the velocity gradient longitudinal component p.d.f., shown in figure 2(a–c). As the Reynolds number increases, the left tail develops, while the right tail becomes slightly sub-Gaussian near the p.d.f. core. Our model captures the slight increase of the left tail up to $\textit {Re}_\gamma =\textit {O}(1)$, even though the non-Gaussianity at such low Reynolds numbers is very mild compared with the non-Gaussianity in high-Reynolds-number turbulence (figure 2d). The model can capture the p.d.f. for approximately six decades for $\textit {Re}_\gamma \leq 1$, but the quantitative agreement deteriorates with increasing $\textit {Re}_\gamma$, and we can approximate the p.d.f. of $A_{11}$ for only three decades at $\textit {Re}_\gamma =2.5$ due to the onset of an unphysical right tail. While the skewness of the velocity gradient statistics is evident already at low $\textit {Re}_\gamma$, the signatures of intermittency show up only at larger Reynolds numbers. This reflects in slight changes of the velocity gradient off-diagonal component p.d.f., as shown in figure 2(e–g). Some intermittency is noticeable only at $\textit {Re}_\gamma >1$, and again very mild compared with the intermittency levels of high-Reynolds-number turbulence (figure 2h). Our analytic prediction can capture the slight increase of the tails of the p.d.f. and respects the symmetry (lack of skewness) of the transverse component statistics.

Figure 2. The p.d.f. of the longitudinal (a–d) and transverse (e–h) components of the velocity gradient at various Reynolds numbers. The solid lines are from DNS, while the dashed black lines indicate the analytic prediction from our low-Reynolds-number model. The thinner dashed line indicates the zeroth-order Gaussian p.d.f.

5.3. Strain-rate statistics

In order to better understand the origin of the strain-rate skewness and to get insight into the preferential configuration of the strain-rate eigenvalues observed in fully developed turbulence (Betchov 1956), we analyse the p.d.f. of the strain-rate tensor using elements of the theory of random matrices (Livan et al. 2018). At zero Reynolds number, the velocity field is multipoint Gaussian, and the strain rate is an orthogonal Gaussian random matrix. Gaussian random matrices are well-established tools in many research areas (Wigner 1955; Dyson 1962), but less frequently employed in fluid dynamics. For example, the distribution of the dissipation rate corresponding to a Gaussian strain rate (i.e. the p.d.f. of the norm of a Gaussian symmetric matrix Wigner 1955) has been more recently derived in this field (Shtilman, Spector & Tsinober 1993; Gotoh & Yang 2022).

Due to incompressiblity and statistical isotropy, the strain-rate p.d.f. is parameterized through two quantities, namely the two independent (unordered) strain-rate eigenvalues $\lambda _{1}$ and $\lambda _{2}$, or the two principal invariants $\mathcal {I}_1 = 2(\lambda _1^2+\lambda _2^2+\lambda _1\lambda _2)$ and $\mathcal {I}_3 = 3\lambda _1\lambda _2(\lambda _1+\lambda _2)$. Therefore, changing coordinates from the standard Cartesian basis to the strain-rate eigenframe brings a remarkable dimensionality reduction. This change of coordinates implies that the p.d.f. of the strain-rate p.d.f. $f_S(\lambda (\boldsymbol{\mathsf{S}}))$ is related to the p.d.f. of its eigenvalues $\rho _{\lambda } (\lambda )$ through (Livan et al. 2018)

(5.3)\begin{equation} f_S(\lambda) = \frac{\rho_\lambda(\lambda)}{J_S(\lambda)}, \end{equation}

where $J_S$ is the volume element of the Stiefel manifold (James 1977), defined in (4.6). Since the volume element $J_S$ has a convoluted form, the p.d.f. of the strain rate parameterized through the eigenvalues, $f_S(\lambda )$, is visually much simpler than the p.d.f. of the eigenvalues themselves $\rho (\lambda )$. We will therefore plot and discuss $f_S(\lambda )$, shown in figure 3 for various Reynolds numbers. Furthermore, we employ the rescaled sum and difference of the eigenvalues, so that the contours of the strain-rate p.d.f. at vanishingly small Reynolds number are circles. The contours of the p.d.f. in figure 3 are logarithmically equispaced, the coloured solid lines are from DNS and the black dashed lines represent the analytic prediction obtained by integrating out the vorticity from (4.1) in the fashion of (4.4)

(5.4)\begin{align} f_S(\lambda) &= \left( \frac{50 \sqrt{30}}{{\rm \pi}^{{5}/{2}}} -\textit{Re}_\gamma \frac{60\,000 \sqrt{30} S_{3}}{7 {\rm \pi}^{{5}/{2}}} \lambda_{1} \lambda_{2} \left(\lambda_{1} + \lambda_{2}\right) + \textit{Re}_\gamma^2 \frac{20\,000 \sqrt{30} S_{3}^{2}}{49 {\rm \pi}^{{5}/{2}}} \right. \nonumber\\ &\quad\times \left.\left( 1800 \lambda_{1}^{4} \lambda_{2}^{2} + 3600 \lambda_{1}^{3} \lambda_{2}^{3} + 1800 \lambda_{1}^{2} \lambda_{2}^{4} - 42 \lambda_{1}^{2} - 42 \lambda_{1} \lambda_{2} - 42 \lambda_{2}^{2} + 7\right)\vphantom{\frac{60\,000 \sqrt{30} S_{3}}{7 {\rm \pi}^{{5}/{2}}}}\right)\nonumber\\ &\quad\times\exp({- 10 \lambda_{1}^{2} - 10 \lambda_{1} \lambda_{2} - 10 \lambda_{2}^{2}}). \end{align}

The analytic prediction (5.4) captures the strain-rate p.d.f. up to $\textit {Re}_\gamma =\textit {O}(1)$. Discrepancies appear in the zero-measure regions in which two of the strain-rate eigenvalues coincide. This could be due to the choice of the basis tensors (2.9), which are not linearly independent when the strain rate is in a degenerate configuration.

Figure 3. The p.d.f. of the strain rate $f_S(\lambda )$ parameterized through its unordered and rescaled eigenvalues $\lambda _i$, at various Reynolds numbers. The colourmap and coloured solid contours are from DNS, while the dashed black lines indicate the corresponding analytic prediction from our low-Reynolds-number model. The thin dashed lines in the high-Reynolds-number plot indicate the empirical parameterization (5.5). The colourmap is $\log _{10}$ scale, and the contours are equispaced in $\log _{10}$ scale with unit increments.

As the Reynolds number increases, the contours of the strain-rate p.d.f. in figure 3 transition from a circular to a triangular shape. The contours elongate towards large and positive $\lambda _1+\lambda _2$, that is, large and negative $\lambda _3$, while shrinking along $\lambda _1-\lambda _2$, especially when $\lambda _3$ is large. This shape indicates a preferential state with a negative and large eigenvalue, with the other two being smaller, positive and close to each other. This configuration has been observed by Betchov (1956) from a phenomenological analysis of the strain-rate dynamics in high-Reynolds-number turbulence, and (5.4) analytically shows the onset of this preferential configuration. The skewness in the p.d.f. (5.4) is due to the first-order contribution in $\textit {Re}_\gamma$, which amplifies the probability of regions in which the product $S_3\lambda _1\lambda _2\lambda _3$ is positive. This corresponds to a preferentially positive intermediate strain-rate eigenvalue $\lambda _2$ (Lund & Novikov 1992).

Figure 3 also shows that the shape of the strain-rate p.d.f.s at low and moderate Reynolds numbers qualitatively differ. The edges displayed by the strain-rate p.d.f. at moderately large Reynolds numbers are not present at lower Reynolds numbers and constitute a high-Reynolds-number feature. The symmetry of the p.d.f. with respect to the axes $\lambda _1=\lambda _2$, $\lambda _1=-\lambda _2/2$ and $\lambda _2=-\lambda _1/2$ remains at all Reynolds numbers, simply because the eigenvalues can be arbitrarily interchanged. Even at large Reynolds numbers however, the contours of the strain-rate p.d.f. look surprisingly simple. In fact, they can be parameterized through a combination of the principal invariants, namely

(5.5)\begin{equation} \alpha_1\left(f_S\right)\mathcal{I}_1^3 + \alpha_2\left(f_S\right)\mathcal{I}_3\mathcal{I}_1^{3/2}+ \alpha_3\left(f_S\right)\mathcal{I}_3^2=1, \end{equation}

where the quantities $\alpha _i$ depend on the contour level $f_S$ itself. This parameterization is shown in figure 3 for the large-Reynolds-number strain-rate p.d.f.

5.4. Vorticity statistics

We now analyse the distribution of the vorticity components in the strain-rate eigenframe and explore the non-monotonic trend of the strain rate–vorticity preferential alignments with increasing Reynolds number observed in figure 1(b). As done for the strain rate in the previous section, we compensate the kinematic effects due to the volume element $J_\omega$ (4.7) that arises from employing magnitude/orientation coordinates in the strain-rate eigenframe. More specifically, we integrate out the strain rate and the vorticity magnitude from the velocity gradient p.d.f., according to (4.4), to obtain the marginal p.d.f. of the normalized vorticity principal components

(5.6)\begin{equation} \rho_{\hat{\omega}}(\hat{\omega}^2) = \int\textrm{d}\omega \, J_\omega \int\textrm{d}\boldsymbol{\mathsf{S}}\,f(\mathcal{I}\left(\boldsymbol{\mathsf{S}},\boldsymbol{\omega})\right), \end{equation}

where the integration with respect to the strain rate is defined in (4.4a). We then weigh the p.d.f. of the alignments $\rho _{\hat {\omega }}$ by the denominator of the volume element $J_\omega$ (4.7), which does not depend on the vorticity magnitude and can be pulled out of the integral in (5.6), thus yielding

(5.7)\begin{equation} f_{\hat{\omega}}(\hat{\omega}^2) = \mathcal{N}\sqrt{\hat{\omega}_1^2\hat{\omega}_2^2\hat{\omega}_3^2} \rho_{\hat{\omega}}(\hat{\omega}^2) , \end{equation}

where $\mathcal {N}$ is a normalization factor such that $\int \textrm {d}\hat {\omega }_1^2\,\textrm {d}\hat {\omega }_2^2\ f_{\hat {\omega }}=1$. The advantage of such reweighting of the p.d.f. is that it compensates for the complicated features of the volume element $J_\omega$, resulting in visually simple trends of the alignment p.d.f., as shown in figure 4.

Figure 4. The p.d.f. of the squared normalized vorticity principal components weighted by the volume element, as in (5.7), at various Reynolds numbers. Here $\hat {\omega }_1$ and $\hat {\omega }_2$ are the normalized vorticity components along the most extensional and intermediate strain-rate eigendirections, respectively. The colourmap and coloured solid contours refer to the DNS results, while the black dashed lines refer to the corresponding analytic prediction (5.8). The numbers on the contours indicate the value of the p.d.f. on that contour level.

In figure 4 we compare the reweighted p.d.f. of the alignments $f_{\hat {\omega }}$ from DNS (coloured lines) with our low-Reynolds-number model prediction (black dashed lines). The asymptotic prediction of $f_{\hat {\omega }}$ follows from integrating out the strain rate and vorticity magnitude from the asymptotic solution (4.1) and it reads

(5.8)\begin{align} f_{\hat{\omega}}(\hat{\omega}^2) &= 2 \,{-}\, \textit{Re}_\gamma\frac{18 \sqrt{30} S_{3} }{5 \sqrt{\rm \pi}}(2 \hat{\omega}_{1}^2 + \hat{\omega}_{2}^2 - 1) +\textit{Re}_\gamma^2\left(\vphantom{\frac{59856 S_{3}^{2} \hat{\omega}_{2}^2}{245}} 228 S_{3}^{2} \hat{\omega}_{1}^4 + 228 S_{3}^{2} \hat{\omega}_{1}^2 \hat{\omega}_{2}^2 \right.\nonumber\\ &\quad \left.+\, 228 S_{3}^{2} \hat{\omega}_{1}^2 + 66 S_{3}^{2} \hat{\omega}_{2}^4 - \frac{59\,856 S_{3}^{2} \hat{\omega}_{2}^2}{245} + \frac{21\,912 S_{3}^{2}}{245} - \frac{10\,692 X_{5} \hat{\omega}_{2}^2}{49} + \frac{3564 X_{5}}{49}\right), \end{align}

where, the vorticity components are ordered according to the corresponding strain-rate eigenvalue, with $\lambda _1\geq \lambda _2\geq \lambda _3$, and two of the normalized vorticity principal components suffice to parameterize the alignment p.d.f. since $\sum _i\hat {\omega }_i^2=1$.

At very small $\textit {Re}_\gamma$, the alignment distribution in figure 4 is close to random uniform, corresponding to an almost constant $f_{\hat {\omega }}$ and to the lack of preferential alignments. At low Reynolds numbers, the probability density reweighted by the volume element, $f_{\hat {\omega }}$, displays almost straight contours, which have a slope such that the p.d.f. takes larger values where $\hat {\omega }_1^2$ is large. The analytic solution (5.8) shows that the contours at first order in the Reynolds number consist indeed of straight lines with slope $-2$ in the $\hat {\omega }_1^2$–$\hat {\omega }_2^2$ plane. This results in the preferential alignment between the vorticity and the most extensional strain-rate eigendirection, as previously observed in figure 1. Moreover, the first-order correction enhances the probability of regions in which $\hat {\omega }_1^2+\hat {\omega }_2^2$ is large since $S_3<0$. This shows a connection between the negative skewness of the strain-rate statistics and the lack of alignment between the vorticity and the most compressional strain-rate direction. As the Reynolds number increases, the contours of the alignment p.d.f. 4 tilt toward preferentially large $\hat {\omega }_2^2$, that is a slope larger than $-1$ in the $\hat {\omega }_1^2$–$\hat {\omega }_2^2$ plane. The contours remain approximately straight at low Reynolds numbers, with mild variations of the reweighted p.d.f. magnitude on its support. The analytic solution (5.8) quantitatively captures the p.d.f. and the tilting of the contours with increasing $\textit {Re}_\gamma$, the tilting stemming from second-order terms in $\textit {Re}_\gamma$. At large Reynolds numbers, the p.d.f. varies more rapidly on its support, and its contours visually deviate from straight lines. However, the low- and high-Reynolds-number reweighted p.d.f. of the alignments still look remarkably similar, up to a tilting of the contours.

The tilting of the contours of the normalized vorticity principal components p.d.f. toward preferentially large $\hat {\omega }_2^2$ reflects in changes of the marginal p.d.f.s of the cosines $\hat {\omega }_i=|\boldsymbol {v}_i\boldsymbol {\cdot }\boldsymbol {\omega }|$, where $\boldsymbol {v}_i$ are the ordered strain-rate eigenvectors. In this case, we could not compute the analytic approximation to the marginal p.d.f. of the alignments due to the technical difficulty introduced by the volume element $1/|\hat {\omega }_1\hat {\omega }_2\hat {\omega }_3|$ featured in the integrals (4.4). Instead, we compute the p.d.f. of the alignments by numerically solving the Langevin equation associated with the FPE (2.8) for an ensemble of velocity gradients

(5.9)\begin{align} \frac{\textrm{d} \boldsymbol{\mathsf{A}}}{\textrm{d} t} &={-}\left[ \left(1 - \textit{Re}_\gamma^2\delta_1\right)\boldsymbol{\mathsf{S}} + \left(1 - \textit{Re}_\gamma^2\delta_2\right)\boldsymbol{\mathsf{W}} + \textit{Re}_\gamma\left(\frac{5}{7}-\delta_3\right)\widetilde{\boldsymbol{\mathsf{SS}}} + \textit{Re}_\gamma\gamma_4\left(\boldsymbol{\mathsf{SW}}-\boldsymbol{\mathsf{WS}}\right)\right. \nonumber\\ &\quad \left.+\,\textit{Re}_\gamma\left(1-\delta_5\right)\left(\boldsymbol{\mathsf{SW}}+\boldsymbol{\mathsf{WS}}\right) + \textit{Re}_\gamma\left(\frac{3}{5}-\delta_6\right)\widetilde{\boldsymbol{\mathsf{WW}}}\right] + \sigma \boldsymbol{\varGamma} , \end{align}

where $\boldsymbol {\varGamma }$ is a white-in-time tensorial noise that has the same (single-point/single-time) statistics of the forcing $\boldsymbol {\nabla } \boldsymbol {F}$ (2.3).

Figure 5 shows the p.d.f. of the normalized vorticity principal components $\hat {\omega }_i$, obtained from DNS and from numerical integration of our low-Reynolds-number model (5.9). At very small $\textit {Re}_\gamma$, there are no preferential alignments, and the distribution of the normalized principal components is random uniform. As $\textit {Re}_\gamma$ increases, the vorticity preferentially aligns first with the most extensional strain-rate principal direction while avoiding alignment with the compressional direction. Since the spinning of the fluid element causes compression along the vorticity direction (Carbone et al. 2020), the lack of alignment between vorticity and the contracting direction hinders the growth of the most compressional velocity gradients. The transient alignment between the vorticity and the extensional strain-rate eigenvector lasts up to $\textit {Re}_\gamma \simeq 5$, and then the well-known preferential alignment with the intermediate eigenvector settles (Ashurst et al. 1987).

Figure 5. Vorticity-strain rate alignments quantified by the p.d.f. of the ordered vorticity principal components at various Reynolds numbers. The solid coloured lines refer to the DNS, while the black dashed lines are from the numerical solution of our low-Reynolds-number model (5.9).

5.5. Velocity gradient principal invariants

Finally, we focus on the joint p.d.f. of the velocity gradient principal invariants, namely

(5.10a,b)\begin{equation} \mathscr{Q} ={-}\tfrac{1}{2}\left(\mathcal{I}_1+\mathcal{I}_2\right), \quad \mathscr{R} ={-}\tfrac{1}{3}\left(\mathcal{I}_3+3\mathcal{I}_4\right), \end{equation}

that is a hallmark in the study of the turbulent velocity gradient dynamics (Meneveau 2011). The volume element characterizing the $\mathscr {R}$–$\mathscr {Q}$ space prevented us from analytic integration of the full p.d.f. (4.1) to obtain the marginal p.d.f. of the principal invariants. As for the strain rate–vorticity alignments presented above, we obtain the marginal $\mathscr {R}$–$\mathscr {Q}$ p.d.f. by numerically solving the Langevin equation (5.9).

The $\mathscr {R}$–$\mathscr {Q}$ p.d.f. entangles all the information on the strain rate and vorticity in a somewhat complicated way. Indeed, integrating out the eigenvectors of $\boldsymbol{\mathsf{A}}$, in order to get the marginal $\mathscr {R}$–$\mathscr {Q}$ p.d.f., does not immediately relate to integrating out rotations of the reference frame, as it does instead for the strain rate. Also, the $\mathscr {R}$–$\mathscr {Q}$ p.d.f. features two kinematically different regions, one below the Vieillefosse line (Vieillefosse 1982) corresponding to real velocity gradient eigenvalues and one above it corresponding to complex eigenvalues. This topological difference causes edges in the p.d.f. of the principal invariants, even though the gradient statistics are smooth.

Nonetheless, the $\mathscr {R}$–$\mathscr {Q}$ p.d.f. displays a distinguishing mark of the velocity gradient statistics: the classical teardrop shape, elongated towards the right Vieillefosse tail (Vieillefosse 1982), as shown in figure 6. At zero Reynolds number, the p.d.f. takes the well-known Gaussian shape, symmetric along the $\mathscr {R}$ axis due to invariance of (2.2) under the sign flipping $\boldsymbol{\mathsf{A}}\rightarrow -\boldsymbol{\mathsf{A}}$. As $\textit {Re}_\gamma$ increases, the characteristic teardrop shape arises. Our model quantitatively captures the onset of the skewness along the right Vieillefosse tail, which persists at higher Reynolds numbers. However, the similarity between the low- and high-Reynolds-number p.d.f.s is only qualitative, the large-Reynolds-number p.d.f. displaying significantly more elongated tails. This is also because the principal invariants are powers of the velocity gradient of degrees two and three, and tiny tails in the strain-rate eigenvalues and vorticity components distributions result in pronounced tails of the $\mathscr {R}$–$\mathscr {Q}$ p.d.f.

Figure 6. Probability density of the velocity gradient principal invariants (5.10a,b) in $\log _{10}$ scale. The colourmap and coloured solid contours refer to the DNS, while the black dashed lines are from the numerical solution of our low-Reynolds-number model (5.9).

6. Comparison of multi-time statistics and velocity gradient dynamics

In the previous section we showed how the asymptotic solution (4.1) quantitatively captures the single-time/single-point velocity gradient statistics obtained from DNS up to $\textit {Re}_\gamma \simeq 1$, with a qualitative agreement up to $\textit {Re}_\gamma \simeq 5$. Now, we investigate whether our model can also reproduce the full Lagrangian dynamics of the velocity gradients, as characterized by velocity gradient sample realizations and time correlations.

6.1. Model coefficients from the DNS

To compare the velocity gradient time series generated by the DNS with our low-Reynolds-number model predictions, we first need to fully specify all the model coefficients (3.6). Indeed, in § 2 we have identified two gauge coefficients, $\gamma _4$ and $\zeta _5$, which only affect multi-time statistics while leaving the single-time p.d.f.s unchanged. We compute the model coefficients (3.6) directly from the DNS and compare them with their analytic predictions. This allows us to fit the two gauge coefficients, $\gamma _4$ and $\zeta _5$, from DNS data, thus enabling our low-Reynolds-number model to capture the velocity gradient temporal dynamics.

From the DNS of low-Reynolds-number random flows, we have access to space–time realizations of the pressure Hessian and viscous Laplacian, i.e. the unclosed terms in our single-time/single-point modelling approach. Therefore, we can compute from DNS the averages of the unclosed terms conditional on the local velocity gradient via tensor function representation theory (Leppin & Wilczek 2020; Carbone & Wilczek 2021)

(6.1a)$$\begin{gather} \left.\left\langle \tilde{\boldsymbol{\mathsf{H}}}(\boldsymbol{x},t) \middle| \boldsymbol{\mathsf{A}} \right\rangle \right. = \sum_{n=1}^8 h_n \boldsymbol{\mathsf{B}}^n, \end{gather}$$

(6.1b)$$\begin{gather}\left.\left\langle \nabla^2\boldsymbol{\mathsf{A}}(\boldsymbol{x},t) \middle| \boldsymbol{\mathsf{A}} \right\rangle \right. + \boldsymbol{\mathsf{A}} = \textit{Re}_\gamma^2\sum_{n=1}^2 \delta_n \boldsymbol{\mathsf{B}}^n + \textit{Re}_\gamma\sum_{n=3}^8 \delta_n \boldsymbol{\mathsf{B}}^n. \end{gather}$$

Here $\boldsymbol{\mathsf{B}}^n$ are the basis tensors (2.9), $h_n$ in (6.1a) are the conditional pressure Hessian components and $\delta _n$ in (6.1b) are the components of the viscous corrections. To be consistent with our low-Reynolds-number model formulation, here we use the same representation of the pressure Hessian and viscous Laplacian that we employed to derive the model coefficients (3.6). From the viscous Laplacian (6.1b), we subtract the linear damping part, $-\boldsymbol{\mathsf{A}}$, and then we introduce second-order corrections in the Reynolds number proportional to the first two coefficients $\delta _1$ and $\delta _2$, and first-order corrections proportional to the other basis tensors. To extract the coefficients $h_n$ and $\delta _i$ from the DNS data, we use the following property of any conditional average (we illustrate this only for the anisotropic pressure Hessian):

(6.2)\begin{equation} \left\langle \left.\left\langle \tilde{H}_{ij}(\boldsymbol{x},t) \middle| \boldsymbol{\mathsf{A}} \right\rangle \right. T_{ij}(\boldsymbol{\mathsf{A}}) \right\rangle = \left\langle \tilde{H}_{ij}(\boldsymbol{x},t) T_{ij}\left(\boldsymbol{\mathsf{A}}(\boldsymbol{x},t)\right) \right\rangle. \end{equation}

Here $\boldsymbol{\mathsf{T}}(\boldsymbol{\mathsf{A}})$ is a tensor function of the gradient only. At most, the model coefficients depend on all the five invariants (2.10), while our main modelling hypothesis is that we limit ourselves to constant coefficients. This hypothesis is exact for the conditional pressure Hessian and the zeroth-order viscous Laplacian at very low Reynolds number, while it constitutes an approximation for higher-order corrections.

Thanks to the property (6.2), we obtain a linear system for the constant model coefficients by taking the double contraction between the tensor representations (6.1) and the basis tensors (2.9) and then averaging

(6.3a)$$\begin{gather} \sum_{n=1}^8 \left\langle Z^{mn} \right\rangle h_n = \left\langle \tilde{H}_{ij}(\boldsymbol{x},t) B^m_{ij}\left(\boldsymbol{\mathsf{A}}(\boldsymbol{x},t)\right) \right\rangle, \end{gather}$$

(6.3b)$$\begin{gather}\textit{Re}_\gamma^2 \sum_{n=1}^2 \left\langle Z^{mn} \right\rangle\delta_n + \textit{Re}_\gamma \sum_{n=3}^8 \left\langle Z^{mn} \right\rangle\delta_n = \left\langle \left(\nabla^2 A_{ij}(\boldsymbol{x},t) + \boldsymbol{\mathsf{A}}(\boldsymbol{x},t)\right) B^m_{ij}\left(\boldsymbol{\mathsf{A}}(\boldsymbol{x},t)\right) \right\rangle, \end{gather}$$

where $Z^{mn}$ is the metric tensor (2.13). Solving the linear system (6.3) yields the coefficients $h_n$ and $\delta _n$, plotted in figure 7 as functions of the Reynolds number $\textit {Re}_\gamma$.

Figure 7. Components of the conditional anisotropic pressure Hessian $h_n$ (a,b) components of the conditional viscous Laplacian $\delta _n$ (see (6.1) for their definition), as functions of the Reynolds number $\textit {Re}_\gamma$. The coloured points indicate the coefficients computed from DNS using (6.3), with the point types differentiating the order of the corresponding basis tensor (2.9). The grey transparent lines are from the analytic prediction (3.6), after fitting the gauge parameters $\gamma _4$ and $\zeta _5$. The model coefficients $\delta _7$ and $\delta _8$ are set to zero and not shown here.

Figure 7(a) shows that the pressure Hessian coefficients $h_3$ and $h_6$ (relative to the basis tensors $\widetilde {\boldsymbol{\mathsf{SS}}}$ and $\widetilde {\boldsymbol{\mathsf{WW}}}$, respectively) constitute the dominant contribution at low Reynolds numbers, with $h_3\to -2/7$ and $h_6\to -2/5$ at vanishingly small $\textit {Re}_\gamma$ (Wilczek & Meneveau 2014). The coefficient $h_7$ is largely compensated by the viscous part, which is beneficial for modelling since we hypothesized that the third-order basis tensors do not contribute to the gradient dynamics. At $\textit {Re}_\gamma \simeq 5$, the coefficient $h_3$, which mitigates the strain-rate self-amplification, becomes larger in magnitude than the coefficient $h_7$, which instead hinders the centrifugal forces due to the rotation of the fluid element. At the same Reynolds threshold, the conditional pressure Hessian component along $\boldsymbol{\mathsf{B}}^1\equiv \boldsymbol{\mathsf{S}}$ starts increasing.

Figure 7(b) shows that the corrections to the viscous linear damping encoded in $\delta _1$ and $\delta _2$ are moderate, and the model qualitatively captures the growth in magnitude of $\delta _1$ and $\delta _2$, which slowly become more negative as the Reynolds number increases. The fact that the coefficients $\textit {Re}_\gamma \delta _1$ and $\textit {Re}_\gamma \delta _2$ remain small at small Reynolds numbers justifies the assumption of quadratic (or, at least, higher order than linear) corrections to the viscous damping. As for the pressure Hessian, also for the viscous stress, the dominant contributions come from the second-order basis tensors, $\boldsymbol{\mathsf{B}}^3$ to $\boldsymbol{\mathsf{B}}^6$. The cubic terms are subleading at small Reynolds number and become relevant only at $\textit {Re}_\gamma \gtrapprox 1$. While the contributions to the symmetric part $\boldsymbol{\mathsf{B}}^7$ from the pressure Hessian and viscous stress compensate each other, the viscous anti-symmetric part proportional to $\boldsymbol{\mathsf{B}}^8$ becomes relevant when $\textit {Re}_\gamma \gtrapprox 1$. This indicates that velocity gradient models at higher Reynolds numbers would need to take into account those third-order basis tensors.

Remarkably, the coefficient $\delta _3$, relative to $\boldsymbol{\mathsf{B}}^3\equiv \widetilde {\boldsymbol{\mathsf{SS}}}$, is uniquely determined by the skewness growth rate at a small Reynolds number $S_3$, as shown by the asymptotic prediction (3.6). The strong relationship between the onset of skewness in the gradient statistics and the coefficient of $\boldsymbol{\mathsf{B}}^3$ has been noticed before in the numerical experiments of Leppin & Wilczek (2020), and here it is shown analytically. Furthermore, both $S_3$ and $\delta _3$ simultaneously match their DNS values (cf. figures 1(a) and 7(b)), which have been obtained in two independent ways ($S_3$ by looking at $\left \langle \mathcal {I}_3 \right \rangle$ as a function of $\textit {Re}_\gamma$, and $\delta _3$ by computing the conditional viscous Laplacian). This supports the consistency of our modelling approach.

Finally, the analytic expressions of the coefficients $\delta _4$, $\delta _5$ and $\delta _6$ in (3.6) feature the gauge terms $\gamma _4$ and $\zeta _5$. While the single-time statistics are independent of those two terms, the time correlations are affected. We fit the values of these gauge terms from DNS data, and with

(6.4a,b)\begin{equation} \gamma_4\simeq{-}0.022, \quad \zeta_5\simeq 0.013, \end{equation}

the analytic prediction (3.6) matches well the coefficients $\delta _n$ independently computed from DNS, as in figure 7(b).

6.2. Lagrangian trajectories

Now that we have fully determined the model coefficients, we analyse the velocity gradient along fluid–particle trajectories generated by our model and by DNS. The model coefficients necessary to reproduce the single-time velocity gradient statistics have been determined analytically. However, fitting DNS data is necessary to set the gauge terms in (3.6), affecting the velocity gradient time correlations. The Langevin equation (5.9) and the associated FPE (2.8) are designed to match single-time/single-point statistics, and the results can be interpreted in either an Eulerian or Lagrangian sense. In the following, we assume a Lagrangian viewpoint.

Figure 8 shows a longitudinal component of the velocity gradient tensor along fluid–particle trajectories at various Reynolds numbers, from the model and from DNS. For a vanishingly small Reynolds number, the noise dominates and the realizations follow a tensorial Ornstein–Uhlenbeck process (Chevillard et al. 2011). As the Reynolds number increases, the effect of the random forcing weakens, time correlations establish and larger negative gradients persist. The symmetry about $A_{11}=0$ breaks as the Reynolds number increases, with large negative gradients being more likely, as evident especially at a large Reynolds number. The larger-Reynolds-number trajectories observed on the time scale of a few Kolmogorov times appear smooth, indicating that the details of the stochastic forcing are concealed by the rich small-scale turbulent structure of the flow. The velocity gradient realizations from our low-Reynolds-number model are visually similar to the realizations from DNS up to $\textit {Re}_\gamma =\textit {O}(1)$, and we are now going to quantify this similarity by comparing their time correlations.

Figure 8. Time evolution of the longitudinal velocity gradient component at low Reynolds number, from DNS (a,c) and from our low-Reynolds-number model (2.2) (b,d). The bottom panel (e) shows trajectories at a moderately large Reynolds number. Curves at different transparency refer to various samples in the same simulation. The plots share the same horizontal axis, with the dimensional time normalized by the Kolmogorov time scale, $\bar {t}/\bar {\tau }_\eta =\textit {Re}_\gamma t$.

6.3. Time correlations from DNS and from the low-Reynolds-number model

In figure 9 we compare the normalized time correlations of the strain and rotation rates

(6.5a,b)\begin{equation} C_{\boldsymbol{\mathsf{S}}}(t) = \frac{\left\langle S_{ij}(t)S_{ij}(0) \right\rangle}{\left\langle S_{ij}(0)S_{ij}(0) \right\rangle}, \quad C_{\boldsymbol{\mathsf{W}}}(t) = \frac{\left\langle W_{ij}(t)W_{ij}(0) \right\rangle}{\left\langle W_{ij}(0)W_{ij}(0) \right\rangle}, \end{equation}

as obtained from the DNS and from our low-Reynolds-number model (5.9). At very small $\textit {Re}_\gamma$, the velocity gradient follows a linear Langevin equation so that the correlations decay exponentially in time, with characteristic time scale $\bar {\gamma }_0^2/\bar {\nu }$. Therefore, in the non-dimensional time $t$ (2.5a–d), the correlations collapse at small Reynolds numbers, as evident from the insets in figure 9. However, as the Reynolds number increases, the correlations decay with a time scale shorter than $\bar {\gamma }_0^2/\nu$, indicating that the nonlinearities reduce the normalized time correlations expressed in the non-dimensional variables (2.5a–d).

Figure 9. Normalized time correlations of the strain rate (a) and rotation rate (b), for various Reynolds numbers, as functions of the time lag normalized by the Kolmogorov time scale, $\bar {t}/\bar {\tau }_\eta =\textit {Re}_\gamma t$. Solid lines are from DNS, while the symbols refer to our low-Reynolds-number model (2.2). The insets show the correlations in a semi-logarithmic scale, as a function of the non-dimensional time lag, $t=(\bar {\nu }/\bar {\gamma }_0^2)\bar {t}$. The black dashed lines in the insets indicate the expected zero-Reynolds-number correlation, $C_{\boldsymbol {A}}= \exp (-t)$.

Furthermore, the strain rate is correlated over shorter time scales, as in figure 9(a), and the bulk of the correlations already establishes at $\textit {Re}_\gamma = \textit {O}(10)$, that is, just after the transition to turbulent configuration begins. Instead, the vorticity correlates over longer time scales and those long-lasting correlations only establish at relatively large Reynolds numbers. The extreme vorticity intermittency and long-lasting correlations then constitute the main missing ingredients in low-Reynolds-number random flows, as compared with turbulent flows. This is consistent with recent observations (e.g. Ghira, Elsinga & da Silva 2022) that the most intense vortical structures fully develop only at very large Reynolds numbers.

Our model can quantitatively predict the velocity gradient time correlations up to $\textit {Re}_\gamma =\textit {O}(1)$, and to this end, the gauge terms discussed in the previous sections are crucial. Moreover, the model performs better on the vorticity time correlations than on the strain-rate correlations. The difficulty in predicting the strain-rate time correlations has been recently observed also in reduced-order models for the velocity gradient at high Reynolds number (Leppin & Wilczek 2020).

7. Conclusions

We have derived a model for the velocity gradient in low-Reynolds-number flows governed by the stochastically forced Navier–Stokes equations. The model parameters follow from the homogeneity constraints on the velocity gradient moments and a weak-coupling expansion of the Navier–Stokes equations. The expansion parameter is the Reynolds number $\textit {Re}_\gamma$ (2.6), comparing the magnitudes of the velocity gradient fluctuations and the viscous damping. The model features the exact first-order corrections to the velocity gradient dynamics due to the pressure Hessian at small Reynolds numbers. The Betchov homogeneity constraints $\left \langle \mathcal {I}_1+\mathcal {I}_2 \right \rangle =\left \langle \mathcal {I}_3+3\mathcal {I}_4 \right \rangle =0$, the growth of the strain-rate skewness, $\left \langle \mathcal {I}_3 \right \rangle = S_3\textit {Re}_\gamma$, and fourth-order moment, $\left \langle \mathcal {I}_5 \right \rangle = -1/12 + X_5\textit {Re}_\gamma ^2$, are also imposed exactly (see (2.10) for the definitions of the invariants $\mathcal {I}_k$). Modelling is necessary to devise a closure for the viscous Laplacian contributions and second-order pressure Hessian contributions such that the trend of the average invariants at small $\textit {Re}_\gamma$ hold and the Betchov homogeneity constraints are fulfilled. The resulting model can quantitatively predict the onset of non-Gaussianity in random low-Reynolds-number flows up to Reynolds numbers of order one, and qualitative agreement persists until a transition to the velocity gradient configuration reminiscent of turbulence. Beyond this transition, the model predictions break down, as expected for a perturbative weak-coupling approach.

The asymptotic solution to the FPE associated with our stochastic model is integrable analytically, thus granting access to the moments of the velocity gradients at low Reynolds numbers through closed asymptotic expressions. This allowed us to quantify the onset of hallmark turbulent features, like the alignments between the vorticity and the strain rate, the preferential configuration of the strain-rate eigenvalues and the characteristic teardrop shape of the velocity gradient principal invariants p.d.f. We summarize the main observations below.

At low Reynolds numbers, the results highlight a direct link between a forward energy cascade and preferential alignment between the vorticity and the most extensional direction, both because $S_3<0$. At higher Reynolds numbers, corrections to the fourth-order velocity gradient moments, related to intermittency, lead to the preferential alignment between the vorticity and the intermediate strain-rate direction. The model also shows how the skewness growth rate $S_3<0$ affects the shape of the strain-rate p.d.f., breaking the symmetry of the Gaussian configuration and producing a preferentially large and negative strain-rate eigenvalue. This skewness is also associated with the onset of the teardrop shape of the velocity gradient principal invariants p.d.f. While the first-order corrections to the skewness of the velocity gradient statistics are evident already at $\textit {Re}_\gamma <1$, the onset of the intermittency is slower, quadratic in $\textit {Re}_\gamma$, and evident only at $\textit {Re}_\gamma >1$. The model for the single-time statistics is closed, not requiring any input from the DNS. On the other hand, fitting two gauge parameters from the DNS data is necessary to capture the multi-time statistics of the velocity gradient. After that fitting, our model can produce velocity gradient time series visually and statistically similar to those from DNS at low Reynolds numbers, and we can explore the evolution of the velocity gradient time correlations in terms of $\textit {Re}_\gamma$. The results show that the strain-rate time correlations approach those at high Reynolds numbers already at small Reynolds numbers $\textit {Re}_\gamma =\textit {O}(10)$. Larger differences between low and high Reynolds numbers are observed for the time correlations of vorticity, showing that the long-lasting time correlations of the vorticity constitute a high-Reynolds-number feature.

At moderate Reynolds numbers, $\textit {Re}_\gamma \simeq 5$, we observe the flipping of the strain rate–vorticity alignments towards a configuration reminiscent of turbulence, in which the vorticity preferentially aligns with the intermediate strain-rate eigenvector. Around the same critical $\textit {Re}_\gamma$ intermittency becomes relevant. The transition value $\textit {Re}_\gamma \simeq 5$ corresponds to a Reynolds number based on the Taylor microscale $\textit {Re}_\lambda \simeq 7.5$, consistent with the critical Reynolds number identified via the analysis of the velocity increments/gradients scaling exponents (Yakhot & Donzis 2017; Gotoh & Yang 2022; Khurshid et al. 2023). At very low Reynolds numbers, we have $\textit {Re}_\lambda \simeq 1.5 \textit {Re}_\gamma$, with the proportionality coefficient moderately dependent on the forcing correlations, thus justifying the characterization of the onset of non-Gaussianity via either of the Reynolds numbers $\textit {Re}_\gamma$ or $\textit {Re}_\lambda$.

Finally, both the model and the DNS results show that the transition from a Gaussian state to a turbulence-like configuration is smooth, smeared around a critical Reynolds number. This finding is consistent with previous results (Yakhot & Donzis 2017; Gotoh & Yang 2022; Khurshid et al. 2023) that do not show sharp changes of the velocity gradient statistics in the vicinity of the critical Reynolds number. Based on this underlying smoothness, we interpret the transition to turbulence, as quantified by the change in the scaling exponents of the velocity gradient moments in terms of the Reynolds number (e.g. Yakhot & Donzis 2017), as an approximation to a smooth transition between two power-law trends with different exponents.

As an outlook, the methodology developed here may also be applicable to model the p.d.f. of the velocity gradient starting from the Navier–Stokes equations at higher Reynolds numbers, for example, by employing renormalization techniques (Yakhot & Orszag 1986; Apolinário et al. 2019; Verma 2023).