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 Reference Wilczek and Meneveau2014), 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 Reference Yakhot and Donzis2017; Khurshid et al. Reference Khurshid, Donzis and Sreenivasan2023) 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 (Reference Yakhot and Donzis2017).

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 Reference Wilczek and Meneveau2014),

(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 Reference Meneveau2011). 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 Reference Itskov2015). Statistical isotropy then allows us to employ tensor function representation theory (e.g. Rivlin & Ericksen Reference Rivlin and Ericksen1955; Itskov Reference Itskov2015) to express the unclosed conditional averages in (2.8) as isotropic tensor functions of the velocity gradient (e.g. Pope Reference Pope1975; Novikov Reference Novikov1993). 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 Reference Pennisi and Trovato1987), 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 Reference Lund and Novikov1992) 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 Reference Tian, Livescu and Chertkov2021; Buaria & Sreenivasan Reference Buaria and Sreenivasan2023). 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 Reference Pennisi and Trovato1987). 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. Reference Meurer2017). 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 (Reference Carbone and Wilczek2022).

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.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 Reference Wilczek and Meneveau2014), 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 Reference Wilczek and Meneveau2014)

(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 (Reference Wilczek and Meneveau2014) 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 Reference Wilczek and Meneveau2014; Johnson & Meneveau Reference Johnson and Meneveau2016), 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 Reference Carbone, Iovieno and Bragg2020). 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 Reference Leppin and Wilczek2020), 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 Reference Betchov1956), the perturbation theory at a small Reynolds number for the full Navier–Stokes equations (Wyld Reference Wyld1961), together with the constant dissipation rate imposed by the stochastic forcing (Novikov Reference Novikov1965). 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 Reference Carbone and Wilczek2022). 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 Reference Girimaji and Pope1990; Johnson & Meneveau Reference Johnson and Meneveau2016; Leppin & Wilczek Reference Leppin and Wilczek2020), 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 Reference Wyld1961). 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 Reference Wilczek and Meneveau2014)

(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 (Reference Leppin and Wilczek2020). 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.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. Reference Livan, Novaes and Vivo2018), 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 Reference Wigner1955)

(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.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 Reference Dresselhaus and Tabor1992; Tom, Carbone & Bragg Reference Tom, Carbone and Bragg2021).

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 Reference Betchov1956; Davidson Reference Davidson2015)

(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 Reference Yakhot and Orszag1986), while it weakly increases with the Reynolds number with an exponent of order $0.1$ in numerical simulations (Ishihara, Gotoh & Kaneda Reference Ishihara, Gotoh and Kaneda2009). 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. Reference Ashurst, Kerstein, Kerr and Gibson1987). 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 Reference Vieillefosse1982), 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 Reference Novikov1990; Cantwell Reference Cantwell1992).

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. Reference Ashurst, Kerstein, Kerr and Gibson1987). 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 Reference Yakhot and Donzis2017; Gotoh & Yang Reference Gotoh and Yang2022; Khurshid et al. Reference Khurshid, Donzis and Sreenivasan2023). 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 Reference Yakhot and Donzis2017). 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. Reference Khurshid, Donzis and Sreenivasan2023), 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 Reference Yakhot and Donzis2017), 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 Reference Betchov1956), we analyse the p.d.f. of the strain-rate tensor using elements of the theory of random matrices (Livan et al. Reference Livan, Novaes and Vivo2018). 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 Reference Wigner1955; Dyson Reference Dyson1962), 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 Reference Wigner1955) has been more recently derived in this field (Shtilman, Spector & Tsinober Reference Shtilman, Spector and Tsinober1993; Gotoh & Yang Reference Gotoh and Yang2022).

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. Reference Livan, Novaes and Vivo2018)

(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 Reference James1977), 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 (Reference Betchov1956) 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 Reference Lund and Novikov1992).

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. Reference Carbone, Iovieno and Bragg2020), 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. Reference Ashurst, Kerstein, Kerr and Gibson1987).

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 Reference Meneveau2011). 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 Reference Vieillefosse1982) 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 Reference Vieillefosse1982), 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.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 Reference Leppin and Wilczek2020; Carbone & Wilczek Reference Carbone and Wilczek2021)

(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 Reference Wilczek and Meneveau2014). 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 (Reference Leppin and Wilczek2020), 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. Reference Chevillard, Lévêque, Taddia, Meneveau, Yu and Rosales2011). 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 Reference Ghira, Elsinga and da Silva2022) 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 Reference Leppin and Wilczek2020).