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.


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 2 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;Meneveau 2011;Buaria et al. 2019;Buaria & Pumir 2022).Recent works (Yakhot & Donzis 2017;Sreenivasan & Yakhot 2021;Gotoh & Yang 2022;Khurshid et al. 2022) 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 threedimensional, incompressible Navier-Stokes equations for a statistically isotropic flow without boundaries, driven by a large-scale Gaussian forcing.This setup 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 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 et al. 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  , where  is the Kolmogorov length scale, and at a Reynolds number based on the Taylor micro-scale Re  9, whereas any anomalous scaling is negligible at Re  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 et al. 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 & Pumir 2021).
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 modeling.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 et al. 1999;Chevillard & Meneveau 2006;Wilczek & Meneveau 2014;Johnson & Meneveau 2016;Pereira et al. 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 weakcoupling 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 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 for the velocity gradient probability density function (PDF), 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 Fokker-Planck equation admits asymptotic analytic solutions, thus yielding an analytic approximation to the velocity gradient PDF.
Attempts to analytically predict the high-dimensional velocity gradient PDF are limited and, so far, built upon phenomenological models for the small-scale turbulent dynamics (Chertkov et al. 1999;Moriconi et al. 2014;Apolinário et al. 2019).In general, the analytical expression for the PDF follows from field-theoretical methods employed to solve the nonlinear Langevin equation governing the flow dynamics (Martin et al. 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 PDF 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 PDF 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 random flows.
The paper is organized as follows.Section 2 presents the derivation of a Fokker-Planck equation (FPE) for the velocity gradient PDF from the Navier-Stokes equations.Section 3 specializes the FPE for low-Reynolds flows, while section 4 presents its analytic solution.
The comparison between the model/analytic predictions and low-Reynolds DNS is presented in section 5, while the conclusions and outlook are discussed in section 7. The appendices A and B describe the setup of the numerical simulations and the low-Reynolds number expansion of the Navier-Stokes equations used to determine the model coefficients.

Fokker-Planck equation for the velocity gradient probability density in statistically isotropic flows
In this section, we obtain the general Fokker-Planck equation (FPE) governing the singletime/single-point statistics of the velocity gradient.We then specialize it for flows at low Reynolds numbers.
2.1.Equation governing the velocity gradient dynamics We begin with the three-dimensional incompressible Navier-Stokes equations, driven by an external Gaussian stochastic forcing , (2.1a) where (, ) is the velocity field, (, ) is the pressure field and Re  is the Reynolds number.Equations (2.1) are written in non-dimensional variables suited for the upcoming low-Reynolds number expansion.We will 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 γ0 , related to the damping role of the viscous Laplacian at low Reynolds numbers (see (3.1) for its definition).With γ0 , together with the kinematic viscosity of the fluid ν and the Kolmogorov time scale τ = 1/ √︁ ∇ ū 2 (angle brackets indicating ensemble average), we define the non-dimensional variables employed in equation (2.1) as Based on these quantities, we can define a Reynolds number according to which expresses the ratio between the velocity gradient magnitude τ−1  and the viscous damping ν/ γ2 0 .The Reynolds number (2.3) weighs the nonlinearities in equation (2.1), thus allowing 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 analyze this "onset of non-Gaussianity" at low Reynolds numbers.
The stochastically-driven Navier-Stokes equations (2.1) are associated with a Fokker-Planck equation for the velocity gradient single-time/single-point probability density (Wilczek & Meneveau 2014).To obtain it, we start with the Langevin equation governing the velocity gradient dynamics, obtained from taking the gradient of equation (2.1b), (2.4a) where    = ∇    is the velocity gradient,    = ∇  ∇   is the pressure Hessian and standard matrix product is implied.The Gaussian tensorial noise ∇ in (2.4) has zero mean, is statistically isotropic and white in time.It is fully specified by its correlation, which in Cartesian component notation reads

Focus on Fluids articles must not exceed this page length
where    denotes the Kronecker delta.The spatial correlation of the stochastic forcing is specified in Appendix A.
2.2.Fokker-Planck equation for the velocity gradient invariants Equation (2.4) is associated with a Fokker-Planck equation (FPE) for the probability density function (PDF) of the velocity gradient  (A; ).In Cartesian components, the FPE for an ensemble of fluid particles sharing the same instantaneous configuration of the velocity gradient A reads (Wilczek & Meneveau 2014),

𝜕 𝑓 𝜕𝑡
where •|A denotes the ensemble average conditional on the velocity gradient configuration A. The conditional averages in (2.6), namely the conditional anisotropic pressure Hessian and 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.4) (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.6), come with the cost of introducing unclosed terms.
In statistically isotropic flows, the velocity gradient probability density function  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.6) as isotropic tensor functions of the velocity gradient.The conditional averages are represented as combinations of basis tensors B  with coefficients   that depend upon the velocity gradient invariants I  .More specifically, the basis tensors are formed through the strain-rate tensor S = (A + A )/2 and rotation-rate tensor W = (A − A )/2, and they read where the tilde indicates the traceless/anisotropic part of the tensor, and standard matrix product is implied.It would be necessary to include two additional basis tensors in (2.7) to fix a possible degeneracy of the basis (Pennisi & Trovato 1987).Indeed, some of the tensors (2.7) become linearly dependent when two of the strain-rate eigenvalues coincide or when the vorticity is an eigenvector of the strain rate.We neglect those zero-measure configurations.The independent invariants I  formed through the velocity gradient read (2.8) A sixth invariant would be necessary to fix the handedness of the strain-rate eigenvector basis with respect to the vorticity.

Determining the model coefficients
We introduce the two hypotheses necessary to obtain the presented model: the coefficients   in (2.9) are constant and the expansion (2.9) 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   with the same coefficients computed directly from direct numerical simulations of low-Reynolds number flows (see figure 7).
3.1.Zeroth-order conditional velocity gradient Laplacian At order zero in Re  , the equation governing the velocity gradient dynamics (2.4) is linear, and the resulting flow has Gaussian statistics.For a multi-point Gaussian random field, the conditional Laplacian reduces to a linear damping (Wilczek & Meneveau 2014).Therefore the conditional Laplacian in dimensional variables is of the form and in the non-dimensional variables (2.2), we just have ∇ 2 A|A ∼ −A.The conditional velocity gradient Laplacian for a Gaussian random field contributes to the model coefficients (2.9) at order zero in Reynolds number, while the higher-order viscous corrections require modelling.The reference length γ0 depends on the spatial correlation of the forcing, and we determine it in Appendix B.
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 a simple form (Wilczek & Meneveau 2014) while the local/isotropic part of the Hessian is specified by the incompressibility condition Tr(H) = −Tr(A 2 ).The conditional anisotropic pressure Hessian for a Gaussian random field contributes to the model coefficients (2.9) at order one in Reynolds number.Therefore, we know from Wilczek & Meneveau (2014) the exact first-order correction to the gradient dynamics at small Re  due to the pressure Hessian.
The coefficient ℎ 4 weighing B 4 = SW − WS in (3.2) depends on the structure of the Gaussian flow through the correlation function, 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 B 4 rotates the strain-rate eigenframe while leaving unchanged the vorticity orientation with respect to the eigenframe itself (Carbone et al. 2020).It can also be seen from the Fokker-Planck equation for the gradient PDF (2.10) since 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.11).Therefore, we can ignore the contributions from B 4 as long as we are concerned with single-time/single-point statistics.

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.9), truncated at degree two in the velocity gradient.We consider the basis tensors (2.7) from B 1 to B 6 , by keeping in mind that B 4 can be ignored as long as we focus on single-point statistics.Also, we assume that the coefficients   in (2.9) 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.9) where the constant coefficients   are of order one and are to be determined.We include second-order terms in Re  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 equation (3.3) have been chosen so that the model equations remain unchanged under the transformation t → −t and ν → − ν, 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  1 ,  2 ,  3 ,  5 ,  6 remain to be determined.To do this, we use the two Betchov homogeneity constraints (Betchov 1956), the perturbation theory at small Reynolds 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 where  3 and  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   = 1/ √︁ 2 I 1 is unitary in the non-dimensional variables (2.2).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 been already employed to reduce the number of parameters in numerical simulations of velocity gradient models at high Reynolds numbers (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  3 and  5 represent, respectively, the rate of change of the third-and fourthorder 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 (A), where  (I;   ) is the PDF of the velocity gradient governed by the FPE (2.10), with the drift term (2.9).The FPE (2.10) admits perturbative exact solutions at small Reynolds numbers, as described in further detail in the next section.This analytic solution  (I;   ) depends parametrically on the model coefficients   , 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  1 ,  2 ,  3 ,  5 ,  6 .The solution of this linear system is Additionally, the noise variance  2 = 1/15 is fixed by the unitary Kolmogorov time constraint at Re  = 0 (at which the FPE (2.10) with the drift (3.3) is an Ornstein-Uhlenbeck process).
We assume that the noise variance  2 is independent of the Reynolds number.
We have now determined all the single-time/single-point model coefficients,   and .The model features the exact first-order contribution in Re  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  4 and  5 , which do not affect single-time statistics, but only multi-time correlations.This gauge stems from the fact that the Gaussian and isotropic PDF with unitary Kolmogorov time scale (e.g.Wilczek & Meneveau 2014) for all  4 and  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  4 and  5 in section 6, where we focus on multi-time statistics, with the aid of DNS data.

Analytic approximation of the velocity gradient PDF at small Reynolds numbers
The Fokker-Planck equation (2.10) with the drift term specified by (3.3) admits perturbative solutions in the Reynolds number.The solution is expanded up to second order and we solve for   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.10), 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 Re We remark that while the solution (4.2) is only asymptotic in Reynolds number, the terms   in (4.2) solve exactly, for all I  , the expanded partial differential equation (2.10) 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.2b) by symbolic computation, we rewrote the terms in (4.2c) as functions of the Cartesian components of the velocity gradient, inserted the resulting expression into the Cartesian FPE (2.6), 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 PDF.This allows us to investigate analytically the onset of non-Gaussianity at small Reynolds number.Instead, a shortcoming of this expansion is that the solution is not positive for all I  , as it should be for a proper PDF.However, this positivity breaking happens when the polynomial pre-factors in (4.1) vanish, that is, when the gradients are of the order of Re −1/3  .At such large gradients, the Gaussian factor exp(−5I 1 + 3I 2 ) has already strongly decayed and the consequent error is very small.The comparison with the numerical results in Section 5 will confirm that this positivity issue is indeed negligible at low Reynolds numbers.
An alternative way to obtain an asymptotic solution of the FPE (2.10) consists of the effective action method, proposed in Martin et al. (1973) for stochastic differential equations.This method yields solutions of the form exp(−S(I)) where S is the effective action, featuring higher-order polynomials in I  and Re  .Such a PDF 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 Re  very small, and the equations themselves are only valid for low Reynolds numbers.Therefore, in this setup, the simple expansion yielding the solution (4.2) seems appropriate, and we leave more advanced approaches for future work.

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 I  in terms of the strain-rate eigenvalues   and vorticity principal components   =   • , where   are the strain-rate eigenvectors.Additionally, we use coordinates in the strain-rate eigenframe based on the vorticity magnitude  =  and the alignments between the vorticity vector and strain-rate eigenvectors, ω ≡   /.This procedure transforms the integral (3.5) into with the integration operators in the strain-rate eigenframe specified by and where the invariants are also expressed as functions of the strain-rate variables, 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 A in (3.5)) to five (i.e., the independent strainrate eigenvalues and the vorticity principal components in (4.4)).However, integrating out rotations also comes at the cost of introducing volume elements The absolute value in (4.6) drops when employing the ordered strain-rate eigenvalues  1 >  2 >  3 as integration variables.Due to incompressibility,  1 +  2 +  3 = 0, and we have two possible configurations,  1 >  2 > 0 or  1 > −2 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 and we can consider only positive values of the vorticity principal components ω since the invariants I  (4.5) are even functions of the vorticity.

Comparison of single-point/single-time statistics
We compare the single-time/single-point velocity gradient statistics resulting from our model and from direct numerical simulations at low Reynolds numbers.We point out several qualitative changes in the dynamics and statistical geometry of the gradient as the Reynolds number increases till a transition to turbulence.The DNS setup 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 the Appendices A and B.

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 analyze those three features separately by looking at the strain-rate and vorticity statistics from the strain-rate eigenframe viewpoint (Dresselhaus & Tabor 1992;Tom et al. 2021).
Figure 1 (5.1) The skewness of the strain rate is quantified via the third-order moment of the strain-rate I 3 , while the even and higher-order moments, e.g.I 2 1 , quantify the intermittency.After the rapid initial increase, the skewness of the strain-rate longitudinal component saturates towards 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 et al. 2009).At Re  = 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.I 3 , grow linearly with the Reynolds number, while even-order moments, e.g.I 2 1 , grow quadratically.Recent works (Yakhot & Donzis 2017; Gotoh & Yang 2022), 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.The DNS data indicate that the transition to turbulent scalings starts at Re  ≈ 9 (Yakhot & Donzis 2017), corresponding to Re  ≈ 5 (see figure 10).Our velocity gradients model encodes by construction the different scaling of odd and even moments at low Reynolds number.Indeed, those scalings can be analytically derived through a small-Reynolds-number expansion of the Navier-Stokes equations, as outlined in Appendix B.
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 ω =   • /  (where   are the strain-rate eigenvectors), ordered based on the corresponding strain-rate eigenvalue   , with  1 2  3 .The vorticity aligns with the most extensional strain-rate eigenvector at low Reynolds number, and only at Re  5 it aligns with the intermediate eigenvector.At Re  < 5, the alignments are as if the vorticity was a material line in the fluid flow subject to a persistent strain.Then, at Re  ≈ 5, the transition to turbulence begins: the somewhat counter-intuitive alignment between the intermediate strain-rate eigenvector and the vorticity establishes, as observed 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 simple 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.
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 (5.2) Both ω2 1 and ω2 2 start from a random uniform value ω2  = 1/3 at zero Reynolds number.Then ω2 1 grows linearly with Re  while ω2 2 grows only quadratically.The preferential alignments are closely related to the onset of the strain-rate skewness.Indeed, the fact that  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 equation (5.2).However, ω2 1 has a strong negative second-order contribution in Reynolds and eventually, ω2 2 , which has a positive growth rate, takes over.The analytic approximation (5.2) suggests that the growth of ω2 2 with Re  can only take place if  5 is negative, and its magnitude is large enough compared to  2 in qualitative agreement with the DNS until Re  ≈ 5 and then breaks down, as expected for a weak-coupling perturbative approach.The results highlight many qualitative changes as the Reynolds number increases, consistent with the observations in recent works (Yakhot & Donzis 2017; Gotoh & Yang 2022).However, the Taylor Reynolds number is usually employed to investigate the transition to turbulence in the literature, while we localize the transition at Re  ≈ 5.In our small-Reynolds number analysis, the actual expansion parameter is Re  , which is non-trivially related to the Taylor Reynolds number, their relation depending on the spatial correlations of the forcing.

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 analyze the PDF of the strain-rate tensor using elements of the theory of random matrices (Livan et al. 2018).At infinitesimally small Reynolds number, the velocity field is multi-point 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 has been only recently re-derived in this field (Gotoh & Yang 2022).
Due to incompressiblity and statistical isotropy, the strain-rate PDF is parameterized through two quantities, namely the two independent (unordered) strain-rate eigenvalues  1 and  2 , or the two principal invariants 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 PDF of the strain-rate PDF   ( (S)) is related to the PDF of its eigenvalues   () through (Livan et al. 2018) where   is the volume element of the Stiefel manifold (James 1977), defined in (4.6).Since the volume element   has a convoluted form, the PDF of the strain rate parameterized through the eigenvalues,   (), is visually much simpler than the PDF of the eigenvalues themselves ().We will therefore plot and discuss   (), shown in figure 2 for various Reynolds numbers.Furthermore, we employ rescaled sum and difference of the eigenvalues, so that the contours of the strain-rate PDF at vanishingly small Reynolds number are circles.
The contours of the PDF in figure 2 are logarithmically equispaced, the coloured solid lines are from the DNS, and the dotted black lines represent the analytic prediction obtained by integrating out the vorticity from (4.1) in the fashion of (4.4) (5.4) As the Reynolds number increases, the contours of the strain-rate PDF in figure 2 transition from a circular to a triangular shape.The contours elongate towards large and positive  1 +  2 , that is large and negative  3 , while shrinking along  1 −  2 , especially when  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 equation (5.4) analytically shows the onset of this preferential configuration.The skewness in the PDF (5.4) is due to the first-order contribution in Re  , which amplifies the probability of regions in which the product  3  1  2  3 is negative.This corresponds to a preferentially positive intermediate strain-rate eigenvalue  2 (Lund & Novikov 1992).The analytic prediction (5.4) captures the strain-rate PDF up to Re  = 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.7), which are not linearly independent when the strain rate is in a degenerate configuration.
Figure 2 also shows that the shape of the strain-rate PDFs at low and moderate Reynolds numbers qualitatively differ.The edges displayed by the strain-rate PDF at moderately large Reynolds numbers are not present at lower Reynolds numbers and constitute a high-Reynolds feature.The symmetry of the PDF with respect to the axes  1 =  2 ,  1 = − 2 /2 and  2 = − 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 PDF look surprisingly simple.In fact, they can be parameterized through a combination of the principal invariants, namely where the quantities   depend on the contour level   itself.This parameterization is shown in figure 2 for the large-Reynolds number strain-rate PDF.
The onset of skewed and elongated tails of the strain-rate PDF results in the non-Gaussianity of the strain-rate longitudinal component PDF, shown in figure 3.As the Reynolds number increases, a left tail develops, while the right branch becomes slightly sub-Gaussian near the PDF core.Our model captures the slight increase of the left tail up to Re  = O(1), even though the non-Gaussianity at such low Reynolds numbers is very mild compared to the non-Gaussianity in high-Reynolds number turbulence.

Vorticity statistics
We now analyze 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 remove the kinematic effects due to the volume element   (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 PDF, according to equation (4.4), to obtain the marginal PDF of the normalized vorticity principal components where the integration with respect to the strain rate is defined in (4.4a).We then weigh the PDF of the alignments  ω by the denominator of the volume element   (4.7), which does not depend on the vorticity magnitude and can be pulled out of the integral in (5.6), thus yielding where N is a normalization factor such that The advantage of such reweighting of the PDF is that it removes the complicated features of the volume element   , resulting in visually simple trends of the alignment PDF, as shown in figure 4.
In figure 4 we compare the re-weighted PDF of the alignments  ω from DNS (coloured lines) with our low-Reynolds number model prediction (black dotted lines).The asymptotic prediction of  ω follows from integrating out the strain rate and vorticity magnitude from the asymptotic solution (4.1) and it reads where, the vorticity components are ordered according to the corresponding strain-rate eigenvalue, with  1 >  2 >  3 , and two of the normalized vorticity principal components suffice to parameterize the alignment PDF since  ω2  = 1.At very small Re  , the alignment distribution in figure 4 is close to random uniform, corresponding to an almost constant  ω and to the lack of preferential alignments.At low Reynolds numbers, the probability density re-weighted by the volume element,  ω , displays almost straight contours, which have a slope such that the PDF takes larger values where ω2 1 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 ω2 1 -ω2 2 plane.This results in the preferential alignment between the vorticity and the most extensional strainrate eigendirection, as previously observed in figure 1.Moreover, the first-order correction enhances the probability of regions in which ω2 1 + ω2 2 is large since  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 PDF 4 tilt toward preferentially large ω2 2 , that is a slope larger than −1 in the ω2 1 -ω2 2 plane.The contours remain approximately straight at low Reynolds numbers, with mild variations of the re-weighted PDF magnitude on its support.The analytic solution (5.8) quantitatively captures the PDF and the tilting of the contours with increasing Re  , the tilting stemming from second-order terms in Re  .At large Reynolds numbers, the PDF varies more rapidly on its support, and its contours visually deviate from straight lines.However, the low-and high-Reynolds number re-weighted PDF 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 PDF toward preferentially large ω2 2 reflects in changes of the marginal PDFs of the cosines ω = |  • |, where   are the ordered strain-rate eigenvectors.In this case, we could not compute the analytic approximation to the marginal PDF of the alignments due to the technical difficulty introduced by the volume element 1/| ω1 ω2 ω3 | featured in the integrals (4.4).Instead, we compute the PDF of the alignments by numerically solving the Langevin equation associated with the FPE (2.6) for an ensemble of velocity gradients where  is a white-in-time tensorial noise which has the same (single-point/single-time) statistics of the forcing ∇ (2.5).
Figure 5 shows the PDF of the normalized vorticity principal components ω , obtained from DNS and from numerical integration of our low-Reynolds number model (5.9).At very small Re  , there are no preferential alignments, and the distribution of the normalized principal components is random uniform.As Re  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 till Re  ≈ 5, and then the well-known preferential alignment with the intermediate eigenvector settles (Ashurst et al. 1987).

Velocity gradient principal invariants
Finally, we focus on the joint PDF of the velocity gradient principal invariants, namely that is a hallmark in the study of the turbulent velocity gradient dynamics (Meneveau 2011).
The volume element characterizing the R-Q space prevented us from analytic integration of the full PDF (4.1) to obtain the marginal PDF of the principal invariants.As for the strain rate-vorticity alignments presented above, we obtain the marginal R-Q PDF by numerically solving the Langevin equation (5.9).The R-Q PDF entangles all the information on the strain rate and vorticity in a somewhat complicated way.Indeed, integrating out the eigenvectors of A, in order to get the marginal R-Q PDF, does not immediately relate to integrating out rotations of the reference frame, as it is instead for the strain rate.Also, the R-Q PDF 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 PDF of the principal invariants, even though the gradient statistics are smooth.
Nonetheless, the R-Q PDF 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 infinitesimally small Reynolds number, the PDF takes the well-known Gaussian shape, symmetric along the R axis due to invariance of equations (2.4) under the sign flipping A → −A (Meneveau 2011).As Re  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-Reynolds and high-Reynolds number PDFs is only qualitative, the large Reynolds number PDF displaying significantly more elongated tails.This is also because the principal invariants are moments 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 R-Q PDF.

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 direct numerical simulations up to Re  1, with a qualitative agreement till Re  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.

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 section 2, we have identified two gauge coefficients,  4 and  5 , which only affect multi-time statistics while leaving the single-time PDFs unchanged.We compute the model coefficients (3.6) directly from the DNS data and compare them with their analytic predictions.This allows us fitting the two gauge coefficients,  4 and  5 , from DNS data, thus enabling our low-Reynolds number model to capture the velocity gradient temporal dynamics.
From the direct numerical simulation of low-Reynolds 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) Here B  are the basis tensors (2.7), ℎ  in (6.1a) are the conditional pressure Hessian components, and   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, −A, and then we introduce second-order corrections in the Reynolds number proportional to the first two coefficients  1 and  2 , and first-order corrections proportional to the other basis tensors.
To extract the coefficients ℎ  and   from the DNS data, we use the following property of any conditional average (we illustrate this only for the anisotropic pressure Hessian) where T (A) is a tensor function of the gradient only.At most, the model coefficients depend on all the five invariants (2.8), 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.7) and then averaging 8 ∑︁ =1   ℎ  =    (, )    (A(, )) , (6.3a) where   is the metric tensor (2.11).Solving the linear system (6.3)yields the coefficients ℎ  and   , plotted in figure 7 as functions of the Reynolds number Re  .Figure 7(a) shows that the pressure Hessian coefficients ℎ 3 and ℎ 6 (relative to the basis tensors SS and WW respectively) constitute the dominant contribution at low Reynolds numbers, with ℎ 3 → −2/7 and ℎ 6 → −2/5 at vanishingly small Re  (Wilczek & Meneveau 2014).The coefficient ℎ 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 Re  5, the coefficient ℎ 3 , which mitigates the strain-rate selfamplification, becomes larger in magnitude than the coefficient ℎ 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 B 1 ≡ S starts increasing.
Figure 7(b) shows that the corrections to the viscous linear damping encoded in  1 and  2 are moderate, and the model qualitatively captures the growth in magnitude of  1 and  2 , which slowly become more negative as the Reynolds number increases.The fact that the coefficients Re   1 and Re   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, B 3 to B 6 .The cubic terms are subleading at small Reynolds number and become relevant only at Re  1.While the contributions to the symmetric part B 7 from the pressure Hessian and viscous stress compensate, the viscous anti-symmetric part proportional to B 8 becomes relevant when Re  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  3 , relative to B 3 ≡ SS, is uniquely determined by the skewness growth rate at small Reynolds  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 B 3 has been noticed before in the numerical experiments of Leppin & Wilczek (2020), and here it is shown analytically.Furthermore, both  3 and  3 simultaneously match their DNS values (cf.figure 1(a) and 7(b)), which have been obtained in two independent ways ( 3 by looking at I 3 as a function of Re  , and  3 by computing the conditional viscous Laplacian).This supports the consistency of our modelling approach.
Finally, the analytic expressions of the coefficients  4 ,  5 and  6 in (3.6), feature the gauge terms  4 and  5 .While the single-time statistics are independent of those two terms, the time correlations are affected.We fit the values of those gauge terms from DNS data, and with

Lagrangian trajectories
Now that we have fully determined the model coefficients, we analyze 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 Fokker-Planck equation (2.6) 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 the longitudinal component of the strain rate along fluid-particle trajectories at various Reynolds numbers, from the model and from the DNS.For a vanishingly small Reynolds number, the noise dominates and the realizations follow a tensorial Ornstein-Uhlenbeck process.As the Reynolds number increases, the effect of the random forcing weakens, time correlations establish and larger negative gradients persist.The symmetry about  11 = 0 breaks as Reynolds 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 Re  = O(1), and we are now going to quantify this similarity by comparing their time correlations.as obtained from the DNS and from our low-Reynolds number model (5.9).At very small Re  , the velocity gradient follows a linear Langevin equation so that the correlations decay exponentially in time, with characteristic time scale γ2 0 /.Therefore, in the non-dimensional time  (2.2), 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 γ2 0 /, indicating that the convective nonlinearities reduce time correlations.Furthermore, the strain rate is correlated over shorter time scales, as in figure 9(a), and the bulk of the correlations already establishes at Re  = 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 to turbulent flows.This is in agreement with recent observations (e.g.Ghira et al. 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 Re  = 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).

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 follows directly from a weakcoupling expansion of the Navier-Stokes equations, combined with homogeneity constraints on the velocity gradient moments, and it is asymptotically exact at small Reynolds numbers.It features the exact first-order corrections to the velocity gradient dynamics due to the pressure Hessian, while the viscous contributions are subject to modelling hypotheses.
We have explored the rich statistical geometry of the velocity gradient tensor, by comparing the results from direct numerical simulation and our model.The model can quantitatively predict the onset of non-Gaussianity in random low-Reynolds number flows up to unitary Reynolds number, and qualitative agreement persists till a transition to the velocity gradient configuration reminiscent of turbulence.Beyond this transition, the model predictions break down, as expected for a weak-coupling perturbative approach.Single-time statistics follow from the model without inputs from the DNS.On the other hand, capturing the time correlations requires fitting two gauge parameters from the DNS data.
The asymptotic solution to the model Fokker-Planck equation associated with our stochastic model is integrable analytically, thus characterizing the onset of skewness, alignments and intermittency at low Reynolds number 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 PDF.
The methodology developed here may also be applicable to model the PDF of the velocity gradient starting from the Navier-Stokes equations at larger Reynolds numbers, for example by employing renormalization techniques (Yakhot & Orszag 1986;Apolinário et al. 2019).
where F indicates the spatial Fourier transform,    ( ) =    −     / 2 is the projection tensor on the plane orthogonal to the wavevector k, k = k is the wavevector norm, i is the imaginary unit, and  is the transformed velocity field, The domain length is L = 2 in code units, and the integer wavevectors components range in the interval −/2 < k < /2, where  is the number of resolved Fourier modes in each direction (Canuto et al. 2006).The grid spacing in Fourier space in each direction is Δ k = 2/ L (and Δ k = 1 in code units).The low-Reynolds number simulations resolve 64 3 Fourier modes, with spatial resolutions ranging from kmax η 40 to kmax η 5 with increasing Reynolds number Re  .Here η = ( ν3 / ε) In all the simulations, the complex Gaussian random forcing in (A 1) is limited to low wavenumbers, and it has the form where ā and b are real, vectorial, white-in-time Gaussian random processes with zero mean and unitary variance, ā The sum in (A 3) is extended to the forced wavenumbers k , 1 k < K, with K = √ 7 (in code units).The noise amplitude in (A 1) is proportional to the wavenumber, so that the kinetic energy spectrum   is proportional to  2 at low wavenumbers.For Re  → 0, equation (A 1) tends to an Ornstein-Uhlenbeck process in which the Fourier modes of the velocity evolve independently from each other and have variance  2 ∝ σ2 0 , independent of the viscosity and the wavevector.Therefore, the three-dimensional velocity spectrum scales as   ∼  2 0  2 , compatible with thermal equilibrium of the lowest Fourier modes.
By employing the pre-factor (Δ k) 2 in the definition (A 3), σ0 has the physical units of inverse time and it determines the Kolmogorov time scale of the flow, independently of the domain size L. Furthermore, the forcing (A 3) is in good approximation statistically isotropic, and the correlations of its gradients approximate the isotropic correlation of the forcing (2.5) employed in our low-Reynolds number model.By substituting the definition (A 3) into equation (2.5), it follows that the dimensional forcing amplitude σ0 used in the DNS, and the non-dimensional amplitude  = 1/ √ 15 employed in the Fokker-Planck equation (2.6), are related through where γ0 is the characteristic scale of the damping (B 12), and K is the maximum forced wavenumber.We compute the DNS parameter σ0 that gives a unitary Kolmogorov time scale in the following Appendix B. The simulation parameters are summarized in table 1.
To get an idea of the parameter range and flow configuration, we report in figure 10 the velocity spectra together with a few characteristic Reynolds numbers and scale ratios.Figure 10(a) shows the kinetic energy spectra from the simulations at various Reynolds numbers.As the Reynolds number increases, energy cascades towards the small scales and the high-Plugging the series expansion of the velocity field (B 3) into the Navier-Stokes equations (B 4) gives the following relations for the series coefficients  (  ) , where , ,  are non-negative integers and the propagators read The external Gaussian forcing is delta-correlated in time, with zero mean and correlation where  = γ0 k, the wavevectors k have integer components, and the forced wavevectors belong to the interval {1 k 2 < 7} (cf.DNS setup in Appendix A).The discrete forcing (B 8) is the Fourier transform with respect to time of the forcing employed in the DNS (A 3), in which the wavevectors are discrete by construction.This discrete forcing F is only approximately isotropic, and the spatial correlation of the forcing employed in the Wyld expansion and in the DNS (B 8) corresponds only approximately to the isotropic form (2.5) used in our low-Reynolds number model.However, for the chosen forced wavenumbers range, the discrepancies between the actual forcing correlation and its isotropic limit are tiny.

B.2. Moments of the velocity gradient at order zero in the Reynolds number
We are interested in the velocity gradients moments resulting from the expanded equations (B 6), which can be derived from the Fourier-space velocity field as follows.The Fourier transform of the strain rate is a linear operator on the velocity field, and at any order in the Reynolds number we have The zeroth-order strain-rate magnitude follows from equation (B 6a) combined with the definition of the forcing correlation (B 8), where  =   3 /(2) 4 is the convolution normalization factor.The result (B 10) is valid not only at order zero, but for all Reynolds numbers since the random forcing F produces a constant variance of the gradients (Furutsu 1963;Novikov 1965), independent of the weight of the nonlinearities in (2.1).Imposing a unitary Kolmogorov time scale, I 1 = 1/2, yields σ0 0.027, (B 11) that is the parameter employed in our DNS (cf.table 1).Finally, we need to compute the damping parameter γ0 from the dimensional velocity gradient field since it constitutes the reference length in our low-Reynolds number model.This reference length follows from the conditional Laplacian (B 12) in a Gaussian random where the wavevector norm k and the coefficient γ0 are dimensional (in code units).
The expansion parameter occurring in the asymptotic solution of the FPE (4.2) Re  , is proportional to γ2 0 , as in (2.3).The relation between Re  and the Reynolds number based on the Taylor microscale Re  , that is usually employed in the literature to investigate the onset of non-Gaussianity (Gotoh & Yang 2022), depends on the details of the forcing.
We tackle the frequency-wavenumbers integrals (B 17) by means of a combination of symbolic calculus and discrete integration.The integration with respect to the frequency is carried out analytically by using the residue theorem (Kleinert & Schulte-Frohlinde 2001).On the other hand, the integration over the wavevectors reduces to a discrete summation over the discrete forced wavevectors, due to the form of the forcing correlation (B 8), consisting of a superposition of Dirac delta functions.Evaluating the integral (B 17), and the corresponding one for  5 , we finally get  3 −0.027, 5 −0.0010, (B 18) matching well the trend of the moments from our DNS at low Re  (cf.figure (1)).
For the parameters  5 the procedure is the same as sketched for  3 , but the expansion goes up to second order in the Reynolds number and there is one additional integration variable.Also, the average of the sixth order moments of  (0) now splits into fifteen pairs by the Wick theorem (differently from equation (B 16), in which we got just three pairs).Therefore, the direct computation of integrals like (B 17) presented here is computationally expensive, and the (equivalent) numerical integration of the expanded Navier-Stokes equations (B 4) is preferable.This integration yields the fields ũ( ) (B 3), from which one can compute the moments of the gradient in physical space through Fast Fourier Transform (the FFT requiring just 13 3 grid points for these low-order moments and the forcing (B 8)).Systematic study of expressions like (B 17), aiming to an efficient symbolic/numeric evaluation of these integrals, is the subject of ongoing work.
(a)  shows the normalized moments of the strain-rate longitudinal components as a function of the Reynolds number Re  .The moments of the invariants are related to the moment of the longitudinal strain-rate component through(Betchov 1956; Davidson 2015)

Figure 1 :
Figure 1: Onset of skewness, intermittency and alignments at low Reynolds number, in terms of the expansion parameter Re  .(a) Normalized moments of the strain-rate longitudinal component  11 .Note that we plot (− 11 )  since the odd-order moments are negative.(b) Alignments between the vorticity and the strain-rate eigenvectors.Solid coloured lines are from the DNS, while black dotted lines are the analytic predictions from our low-Reynolds number model.

Figure 2 :
Figure 2: PDF of the strain rate   () parameterized through its unordered and rescaled eigenvalues   , at various Reynolds numbers.The colormap and coloured solid contours are from the DNS, while the dotted black lines indicate the corresponding analytic prediction from our low-Reynolds number model.The thin dotted lines in the high-Reynolds number plot indicate the empirical parameterization (5.5).The colormap is log 10 scale, and the contours are equispaced in log 10 scale with unitary increments.

Figure 3 :
Figure 3: PDF of the strain-rate longitudinal component at various Reynolds numbers.The solid lines are from the DNS, while the dotted black lines indicate the analytic prediction from our low-Reynolds number model.The thinner dotted line indicates the zeroth-order Gaussian PDF.

Figure 4 :
Figure 4: PDF of the squared vorticity principal components weighted by the volume element, as in (5.7), at various Reynolds numbers.ω1 and ω2 are the normalized vorticity components along the most extensional and intermediate strain-rate eigendirections, respectively.The colormap and coloured solid contours refer to the DNS results, while the black dotted lines refer to the corresponding analytic prediction (5.8).The numbers on the contours indicate the value of the PDF on that contour level.

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

Figure 6 :
Figure 6: Probability density of the velocity gradient principal invariants (5.10) in log 10 scale.The colormap and coloured solid contours refer to the DNS, while the black dotted lines are from the numerical solution of our low-Reynolds number model (5.9).

Figure 7 :
Figure 7: Components of the conditional anisotropic pressure Hessian ℎ  (a) and (b)components of the conditional viscous Laplacian   (see equation (6.1) for their definition), as functions of the Reynolds number Re  .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.7).The gray transparent lines are from the analytic prediction (3.6), after fitting the gauge parameters  4 and  5 .The model coefficients  7 and  8 are set to zero and not shown here.

Figure 8 :
Figure 8: Time evolution of the velocity gradient longitudinal component at low Reynolds number, from DNS (a,c) and from our low-Reynolds number model (2.4) (b,d).The bottom panel (e) shows trajectories at 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, t/ τ = Re  .

Figure 9 :
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, t/ τ = Re  .Solid lines are from the DNS, while the symbols refer to our low-Reynolds number model (2.4).The insets show the correlations in a semi-logarithmic scale, as a function of the non-dimensional time lag,  = ( ν/ γ2 0 ) t.The dashed black dotted lines in the insets indicate the expected zero-Reynolds number correlation,   = exp(−).