Statistical properties of the coarse-grained velocity gradient tensor in turbulence: Monte-Carlo simulations of the tetrad model

A proper description of the velocity gradient tensor is crucial for understanding the dynamics of turbulent flows, in particular the energy transfer from large to small scales. Insight into the statistical properties of the velocity gradient tensor and into its coarse-grained generalization can be obtained with the help of a stochastic ‘tetrad model’ that describes the coarse-grained velocity gradient tensor based on the evolution of four points. Although the solution of the stochastic model can be formally expressed in terms of path integrals, its numerical determination in terms of the Monte-Carlo method is very challenging, as very few configurations contribute effectively to the statistical weight. Here, we discuss a strategy that allows us to solve the tetrad model numerically. The algorithm is based on the importance sampling method, which consists here of identifying and sampling preferentially the configurations that are likely to correspond to a large statistical weight, and selectively rejecting configurations with a small statistical weight. The algorithm leads to an efficient numerical determination of the solutions of the model and allows us to determine their qualitative behavior as a function of scale. We find that the moments of order n⩽4 of the solutions of the model scale with the coarse-graining scale and that the scaling exponents are very close to the predictions of the Kolmogorov theory. The model qualitatively reproduces quite well the statistics concerning the local structure of the flow. However, we find that the model generally tends to predict an excess of strain compared to vorticity. Thus, our results show that while some physical aspects are not fully captured by the model, our approach leads to a very good description of several important qualitative properties of real turbulent flows.

3 the inertial range, where the influences of viscosity and of the flow forcing are not significant. Our approach follows earlier theoretical works, which showed that significant insight can be gained for understanding the statistical properties of a passive scalar advected by a turbulent flow by following a set of Lagrangian particles in the flow [19]- [21]. In this spirit, the construction introduced in [1] is based on following four Lagrangian points in the flow. This allows us to construct effectively a finite difference version, M, of the velocity derivative tensor, which can be considered as a surrogate coarse-grained velocity gradient tensor 1 . Physically the evolution of M is coupled to the dynamics of the geometry of the set of points, which is described by three vectors, ρ i a , where the subscript a refers to the direction in space, and the superscript i labels the three vectors. Based on a set of reasonable assumptions, a stochastic model, denoted in the following as the 'tetrad model', has been proposed to describe the evolution of M and ρ [1]. Many of the observables introduced in [1] can be measured, in direct numerical simulations (DNS), and now also in well-controlled laboratory flows, thanks to the spectacular progress in particle tracking velocimetry methods. In this context, obtaining reliable solutions of the tetrad model is expected to provide new insight when compared with high-Reynolds-number data in DNS or laboratory flow.
The probability distribution function (PDF) of the two invariants Q = −tr(M 2 )/2 and R = −tr(M 3 )/3 leads to a characterization of the physical properties of the matrix M, in particular concerning the topology of the flow [22] and its dependence on scale. Additionally, the alignment properties between vorticity and the eigenvectors of the rate of the strain matrix provide a useful characterization of the flow properties [23]. The invariants involving the vorticity and the strain rate correlation [24] are sensitive to these alignment properties and a priori to the increased intermittency as the scale r decreases and thus also provide important information on the flow structure as a function of scale.
The stationary solution of the tetrad model, describing the (Eulerian) statistics of the velocity gradient tensor, is formally expressible in terms of a path integral. The properties of M at a scale r can be determined by generating random paths, corresponding to tetrads with a radius of gyration growing up to the integral scale, where the velocity field is assumed to be Gaussian. This algorithm leads to a Monte-Carlo method for determining the solution of the model. Although straightforward in principle, the corresponding algorithm turns out to be immensely inefficient, as an overwhelming majority of the drawn configurations contribute with a minuscule probability. The problem becomes more severe as the scale r becomes smaller. This situation suggests that very few paths actually count.
In these circumstances, a possible strategy consists in identifying the optimal path, that is, to look for semiclassical solutions of the problem [25]. This method leads to solutions that reproduce qualitatively a number of features observed in DNS [25]. A fundamental limitation of the semiclassical method, however, is that it is difficult to estimate quantitatively the error made and to systematically improve the solution. An alternative strategy is to devise a way to explicitly sample appropriately the parts of phase space that contribute most, instead of exploring uniformly the whole space. This approach, known as importance sampling, has been followed here. The main difficulty in implementing it is to find the appropriate way to sample the configuration space. 4 The determination of the PDF of R and Q requires the integration over all possible configurations with a given set of values of (R, Q), and over all the paths connecting the tetrad at scale r to a tetrad at scale L. While straightforward Monte-Carlo implementations proved to be vastly inefficient, an improvement in the spirit of the semiclassical method was proposed in [26]. In practice, however, the resulting algorithm did not allow us to determine reliably the solution for values of r L/4. The failure of the calculations in [26] indicated that a strategy for sampling efficiently the space of configurations had to be introduced. The aim of this paper is to show how to sample efficiently the paths that matter statistically.
The proper idea to solve the problem comes from importance sampling techniques [27]. A major source of inspiration for the implementation here comes from the Pruned-Enriched algorithm developed by Grassberger [28] for various problems of counting chain configurations in polymer physics and other statistical mechanics problems. The main idea consists in generating more configurations ('enriching') around a trajectory that has a strong statistical weight, and removing configurations ('pruning') when their statistical weight becomes too low. In fact, we show here that it is possible to formulate our problem in such a way that the strategy developed in [28] is applicable.
The resulting algorithm allows us to obtain numerical solutions of the model, with reasonable statistical accuracy over a range of scale L/256 r L/2 with a limited amount of computational resources: our runs typically took ∼10 days on a personal computer.
In its present form, the model of [1] involves three dimensionless parameters: α, the effective reduction of nonlinearity, as well as the coefficients C u and C η characterizing the amplitude of the stochastic terms acting on the ρ and on the M equations, respectively. One of the major conclusions obtained with the semiclassical approximation [25], with a slightly different form of the model, was that the value of the coefficient α had the most significant impact on the qualitative behavior of the solution: for α 0.5, the nonlinearity is so weak that the energy transfer has the wrong sign, whereas for α 0.5, the scale dependence of the strain component of M is incorrect. Thus, the semiclassical approximation suggested that the physically consistent value of α was constrained to be around 0.5. The full Monte-Carlo calculation, using the original model [1], suggests a more subtle picture.
The presentation of our results is organized as follows. In section 2, we discuss the formulation of the model. Section 3 presents the numerical Monte Carlo used here, and in particular, the implementation of the importance sampling method. The various moments of order n 4 are introduced in section 4. The scaling properties of the solutions of the model are discussed in section 5; in particular, we show that the scalings of the moments are very close to the predictions of the Kolmogorov K41 theory [29,30] for all the values of the parameters. We then discuss in section 6 the structure of the PDF of the (R, Q) invariants and their dependence on scale and on the various parameters of the model. Section 7 contains our results concerning the densities of invariants (vorticity, strain and vortex stretching) predicted by the model. Finally, we summarize our results and conclude in section 8.

Formulation
The main idea of developing the present model comes from the observation that the motion of a few Lagrangian points in a turbulent flow can lead to simple phenomenological approximations, which can in turn lead to original predictions concerning the Eulerian correlation functions of the flow [19]- [21]. The model, postulated in [1], is formulated in terms of two coupled stochastic differential equations, modeling the evolution of M, the coarse-grained velocity gradient tensor, along with the three vectors ρ i (i = 1, 2, 3), which describe the shape of the tetrad with respect to the center of mass, dM ab dt where the matrix k is the inverse of ρ. In the notation used here, the lower indices a, b, etc refer to the space dimension and the upper indices i, j, etc specify one among the three vectors describing the shape of the tetrad. The deterministic terms on the left-hand side of equations (1) and (2) represent the effect of the dynamics induced by scales of order |ρ|. The matrix , defined by equation (3), is symmetric with a trace 1; it provides a coupling of the geometry with the dynamics of M. Lastly, the stochastic terms in equations (1) and (2),η andû, represent the random effect of the small scales of motion. They are represented by Gaussian, white in time noise terms, with a scale dependence prescribed by dimensional considerations, consistent with the Kolmogorov scaling [29], where ρ 2 = i,a ρ i a 2 . The noise term acting on ρ,û, is assumed here to act only in the direction transverse to the nine-dimensional vector ρ, which is not expected to be a significant restriction.
It is convenient to reformulate slightly the model, by scaling time by (ρ 2 / ) 1/3 . Defining the equations of the problem become dm ab dτ − where the variance of the noise terms simply reads In this form, the equations become scale independent, in the sense that they are invariant under the transformation ρ → λρ.

Determination of statistical properties
The model is used here to compute the statistics of the velocity gradient tensor as a function of scale. To this end, one expresses the steady-state solution of the Fokker-Planck equation for the PDF P(m, ρ, τ ) associated with equations (7)- (10) in terms of the PDF at integral scale, where N is the normalization constant. By using Green's functions, one obtains the expression for the solution as a path integral, where S represents the action: The simple interpretation of the path integral representation equation (12) is that the PDF P(M 0 , ρ 0 ) is obtained by summing over all possible solutions of equations (7) and (8), integrated backwards in time, and starting at t = 0 with the fixed boundary condition M 0 and ρ 0 , up to the time −T at which the value of |ρ| reaches L, the integral scale. The action term defined by equation (13) expresses that the noise terms have the proper variance, as given by equations (9) and (10). This representation invites a Monte-Carlo treatment. In principle, by simply integrating equations (7), (8) and (3) backwards in time, up to the time −T where the norm of ρ reaches L, and by attributing the statistical weight given by equation (11) to the configuration, one can effectively determine numerically the integral path, equation (12).
The algorithm is thus formulated as a 'first passage problem': the solution is recorded whenever the size |ρ| reaches a prescribed value. The property already noticed that equations (7)-(10) do not depend explicitly on scale then leads to the following simplification. To compute the statistics at scale r , one may simply generate the solutions of equations (7)-(10), with |ρ(τ = 0)| = 1, and integrate until the value of |ρ(τ = −T )| = L/r , and compute the statistical weight given by equation (11).

Numerical solution of the tetrad model
The effort to determine reliable solutions of the tetrad model has been inhibited by a serious difficulty; the number of configurations that contribute significantly is very small, the more so as the ratio L/r increases. We analyze here the origin of this difficulty, and propose an algorithm that solves the problem.

Monte-Carlo algorithm: origin of the difficulty
The reason why the Monte-Carlo algorithm is vastly inefficient can be understood with the help of the following reasoning, based on figure 1. Figure 1(a) shows the evolution of the norm |ρ|(−τ )| as the system of equations (7)-(10) is integrated backwards in time, with |ρ| = 1 initially. The signs show the time of first passage when the value |ρ| = 2 and |ρ| = 4, which give the time at which to compute the statistical weight, using equation (11), corresponding to scales r = L/2 and r = L/4. Figure 1(b) shows the trajectory shown in figure 1(a) (full line), as well as two other trajectories (dashed line and dotted line), in the (log 2 (|ρ|), log 2 (|m|)) plane. As expected naively, based, e.g., on Kolmogorov considerations [1,25], the norm of |m| grows as |ρ| grows, plausibly as |ρ| 2/3 . This implies that the statistical weight computed with equation (11) decreases as |ρ| increases. The most important feature seen in figure 1(b) is the large level of fluctuations in the values of |m| from trajectory to trajectory at a given level of |ρ|. The solutions shown contribute, at a ratio L/r = 4 (vertical dashed line in figure 1(b)), to a variation by a factor of ∼4 in the magnitude of |m|. This very large dispersion results in variations of the Gaussian weight, equation (11), over many orders of magnitude. The problem becomes even more severe as the value of |ρ|, hence L/r , increases, because of the overall growth of |m| with |ρ|, which makes the Gaussian weight even more selective.
The trajectories contributing significantly to the statistics are clearly those with the lowest possible values of |m| at a given value of |ρ|. To obtain reliable results, one may ignore trajectories with high values of |m| and focus instead on configurations with low values of |m|. We explain in turn the strategy followed here to achieve this goal.

Focusing on important configurations
One may first note that selecting initial configurations with the lowest possible value of the norm |m|(τ = 0) is very likely to lead to configurations with a larger statistical weight. The importance sampling strategy [27] has been devised to address this sort of difficulty and has been applied in this spirit to compute the PDF in the (R, Q) plane in [26]. The resulting algorithm certainly improves the quality of the statistics at moderate values of |ρ|. It does not allow us, however, to go to scale ratios larger than L/r 8, because the algorithm is eventually plagued by the growing dispersion of the values of |m| as |ρ| increases.
The main idea of this paper is inspired by the Pruned-Enriched algorithm of Grassberger [28]. Our implementation of the original idea goes as follows. Consider, for the sake of simplicity, a set of trajectories with a fixed initial value of m(τ = 0) = m 0 and a fixed value of ρ 0 , corresponding to a regular tetrad with size |ρ 0 | = 1. The norm of m, as well as the dispersion in the values of |m|, increases as the solution is integrated (backwards) in time, hence when |ρ| grows. Among a set of trajectories integrated up to a time −τ ρ 0 , where |ρ|(−τ ρ 0 ) = ρ 0 , it is likely that those with high values of |m(−τ ρ 0 )|, which therefore do not contribute very much at scale L/r = ρ 0 , do not contribute very much either at larger values of |ρ|. In comparison, there is a much greater chance that a trajectory with a low value of |m|(−τ ρ 0 ) may contribute significantly at larger values of |ρ|. Thus, the lower the value of |m|(−τ ρ 0 ), the more the likelihood that the trajectory contributes to the statistics at higher values of |ρ|. In these circumstances, we choose to generate several copies of trajectories that have a low value of |m|(−τ ρ 0 ) and to selectively remove trajectories with a high value of |m|(−τ ρ 0 ) for a set of fixed values of ρ 0 .
The calculation of the averaged statistical weight W m 0 ,ρ thus goes as follows. We begin by selecting a discrete subset of values ofρ n ; in practice, we took typically up to a value ρ = 256 (which corresponds to p = 8 in equation (14)). We generate N c0 configurations, and estimate, using the original algorithm, a (crude) estimate of the statistical weight W m 0 ,ρ n , by averaging over these N c0 configurations. We then generate N c1 configurations, to which we assign at τ = 0 the weight w c = 1. When the solution reaches the next valuē ρ n , we compare the product p of the weight w c of the individual configuration and of the Gaussian term given by equation (11) to the averaged statistical weight at that scale, W m 0 ,ρ n . The Pruned-Enriched algorithm proceeds with the following rules.
(i) If the product p is smaller than a fraction (typically a third) of W m 0 ,ρ n , keep the configuration with probability 1/2, and double the weight w c of the remaining configurations. Update the value of W m 0 ,ρ n with p c and proceed with the integration when appropriate. (ii) If p is larger than 1/3 of W m 0 ,ρ n , but smaller than 3 times W m 0 ,ρ n , simply update the value of W m 0 ,ρ n with p; then proceed with the integration with an unchanged value of w c . (iii) If p is larger than 3 times W m 0 ,ρ n , generate n r new replicated configurations, each of them with the same value of ρ, and m, but each with a weight equal to w c /n r . Update the value of W m 0 ,ρ n with p c , and proceed with the integration with the new n r configuration.
The values 1/3 and 3 used to separate cases (i), (ii) and (iii) were as suggested in [28]. We checked that these values were also appropriate to our problem. The number of new configurations generated in case (iii) was chosen in the range 4-8 in our calculations.
As demonstrated below, this allowed us to solve, to a large extent, the difficulty of spending too much time on irrelevant configurations, thus allowing us to focus only on very useful trajectories.

Numerical methods and convergence
The integration of equations (7)-(10) has been carried out by using a straightforward Euler scheme, formally of the order of 1/2 [31]. The variance of the noise term u depends on ρ, which suggests that the noise in equation (8) is multiplicative. However, the noise u is transverse to ρ, considered as a vector in the (3 × 3)-dimensional space: it acts only on the angular part of the vector ρ, not on its norm. Adding a longitudinal component to the noise term, of the form , would make the noise term in the ρ equation explicitly |ρ| dependent, which would require more sophisticated methods to integrate equations (7)-(10); see [32]. We have checked that adding a longitudinal contribution to the noise in ρ does not change the conclusions of this work significantly.
As figure 1(b) reveals, the value of M grows as |ρ| increases. The term m × d ln(|ρ| 2/3 )/dt in equation (7) suggests a growth of the norm of |m| like |ρ| 2/3 , consistent with the usual Kolmogorov estimates [1,25]. To maintain a proper accuracy in the numerical integration, we gradually decrease the integration time step as |ρ| increases. More precisely, each time the value of ρ reaches a new value ofρ n , the value of the time step is changed to dt n = dt 0 × (ρ 0 ρ n ) 2/3 . The statistical properties we observed did not show any significant variation when dt 0 was varied in the range 10 −3 dt 0 10 −2 . For most of the runs discussed here, we have used a value of dt 0 = 5 × 10 −3 .
Our runs were organized in the following way. In order to compute the PDF in the (R, Q) plane, we first laid out a grid in the (R 0 , Q 0 ) plane, typically of size 91 2 , with a width equal to R = Q = 0.1, symmetric with respect to 0. We then drew a set of random values of m 0 (the initial condition), with a Gaussian probability P i , corresponding to a product of individual Gaussian distributions for the nine elements of the matrix m 0 , with a weight equal to σ 2 m 0 = 1.5, with the constraint that tr(m 0 ) = 0. The generated configuration with m 0 is then weighted with the inverse of the probability of m 0 , P i (m 0 ) [26,27]. The equations of motion equations (7)-(9) are then integrated backwards in time. The corresponding configuration contributes to the probability, conditioned on the values R = −tr(m 3 0 )/3 × (ρ 2 n / ) and Q = −tr(m 2 0 )/2 × (ρ 2 n / ) 2/3 at scaleρ n . As explained in the previous subsection, we begin by generating N c0 configurations, allowing us to generate a coarse estimate of the PDF P(R, Q,ρ n ). We also estimate the conditional average of quantities such as the strain rate, the enstrophy and other invariants of the velocity derivative tensor. Then, we generate a large number N c1 of configurations, using the Pruned-Enriched algorithm, defined above.
Convergence clearly depends on the number N c1 of configurations drawn, as well as on the number n r of replicated configurations of high probability at step (iii). The runs discussed here have been carried out with typically over N c1 = 5 × 10 7 configurations. The number of replicated configurations n r was varied, so the total number of configurations emanating from a single one did not exceed a few thousand. To this end, we varied n r between 4 and 8. Each run typically lasted for a couple of weeks on a work station. Comparing several values of n r for a set of values of α, C u and C η convinced us of the good quality of our results.

Invariants of the velocity gradient tensor
Intermittency in turbulent flows is often discussed in terms of the velocity structure functions, i.e. of the moments of the velocity differences between two points separated by a distance r , and by the scaling properties of the resulting quantities. It is convenient here to use the classification proposed in [24] to discuss in a statistically homogeneous and isotropic flow the moments of order n 4 involving various components of the velocity gradient tensor, A, . . . For n = 2, the φ 2 moments are expressible in terms of the two invariants φ 2 tr(s 2 ) = s 2 and φ 2 ω 2 = ω 2 , where s is the rate of the strain matrix (s = 1 2 (A + A T )) and ω the vorticity (ω a = abc A bc ). Similarly, the third moments are expressible in terms of φ 3 s 3 = tr(s 3 ) and φ 3 ωsω = ω.s.ω . Lastly, the fourth moments φ 4 can be expressed in terms of φ 4 s 4 = tr(s 2 ) 2 , φ 4 ω 2 s 2 = ω 2 tr(s 2 ) , φ 4 ωssω = ω.s.s.ω and φ 4 ω 4 = (ω 2 ) 2 . The arguments of [24] involve only isotropy and incompressibility. The statistical properties of M are expected to be invariant under rotation. Applying the arguments of [24] to the matrix M leads to the conclusion that the moments are expressible, for n = 2, in terms of 2 S 2 = tr(S 2 ) and for n = 3 in terms of and for n = 4 in terms of S.S. and where We recall here that the invariants Q and R are simply expressed in terms of the other second and third invariants, according to The invariants Q and R are simply the coefficients of the characteristic polynomial of M, which makes them particularly useful for discussing the structure of the velocity gradient tensor [22]. In the case of the moments of the velocity derivative A, homogeneity imposes further that φ 2 ω 2 = 2φ 2 s 2 and φ 3 s 3 = −3φ 3 ωsω . This feature does not extend to the coarse-grained velocity gradient tensor, as defined here. In fact, no particular relationship exists between 2 S 2 and 2 2 or between 3 S 3 and 3 .S. . The fact that the coarse-grained velocity gradient tensor defined here does not possess the same homogeneity properties as its small-scale counterpart is a consequence of the definition in [1], as pointed out in [16]. Switching to a definition of the coarse-grained velocity gradient tensor that would satisfy automatically the constraint 2 2 = 2 2 S 2 could conceivably alleviate some of the difficulties concerning the excess of strain generated by the model; see section 5.3.
It is worth noting that each matrix M is defined, up to a rotation, by five numbers (three vorticity and two strain components); hence, among the nine quantities introduced above, only five are independent for a given M. Clearly, tr(S 2 ) 2 , ( 2 ) 2 and tr(S 2 ) 2 can be simply expressed in terms of tr(S 2 ) and 2 . Still the ensemble averages over the flow of ( 2 ) 2 are not obviously related to 2 . The study of all the invariants listed here would provide insights into the fluctuations of the tensor and into its intermittency properties [24].
The various invariants introduced in equations (16)- (18) were determined numerically and are discussed in the following sections.

Scaling laws
We have followed the solution down to a scale ∼L/256, which allows us to discuss the scalings of the solutions of the model, at least in this limited range of scales.

Elementary scaling considerations
Before we present the numerical solutions of the model, we discuss the expected scaling solutions, in the spirit of the Kolmogorov K41 theory [29,30]. According to the K41 theory, the velocity difference between two points u(|ρ|) behaves as a function of distance |ρ| between the two points like u ∝ ( |ρ|) 1/3 . As a result, the coarse-grained velocity gradient tensor M is expected to behave as M ∝ u/|ρ| ∝ ( /|ρ| 2 ) 1/3 .
It is appropriate to compare the scaling of the solutions of the model with the Kolmogorov K41 scaling, as explained here. According to the K41 phenomenology, the moments of the coarse-grained velocity gradient tensor, equations (16) and (17), should behave as These scaling laws have been shown to provide a good approximation for the second and third moments in low-Reynolds-number simulations [25]. To show possible deviations from the Kolmogorov scaling laws, we have plotted the nth moments of M, divided by the predictions of equation (21): n × (|ρ| 2 / ) n/3 . The dependence of these ratios as a function of |ρ| provides an indication of how well the Kolmogorov scaling holds. The K41 theory effectively assumes self-similarity of the velocity field as a function of scale, in the sense that the moments of the velocity differences at scale |ρ|, once expressed in units of the standard deviation at scale |ρ|, are all independent of |ρ|. A consequence of this property is that the structure of the flow, as characterized in the (R, Q) plane, should be independent of |ρ|.
The tetrad model, as defined by equations (1)- (3), formally admits solutions that scale according to the Kolmogorov scaling. Yet, the existence of 'naive scaling' solutions does not necessarily prevent the existence of anomalous scaling [13,21], one of the most significant lessons learned while solving the Kraichnan model of passive scalar advection. Whether the solutions of the tetrad model actually follow the 'naive' scaling law is therefore not a priori obvious. We demonstrate in the following subsections that the deviations from the Kolmogorov scaling and from self-similarity predicted by the model, if any, are relatively small.

Scaling law in the absence of nonlinearity (α = 1)
We begin to discuss here the conceptually simplest case, where the nonlinearity in equation (7) is absent, that is, when α = 1. As we find numerically, and justify analytically, the scaling laws in this case follow the Kolmogorov K41 phenomenology.
It is a simple matter to see that when the noise term in equation (7) is turned off (C η = 0), the quantity m/|ρ| 2/3 is invariant, which corresponds to M = constant in equation (1). As a consequence, this conservation law shows that all moments of M are independent of |ρ| when C η = 0.
The noise term in equation (7) leads, as expected, to a qualitatively very different behavior, as illustrated in figure 2, that shows the scale dependence of one of the (compensated) second ( 2 ) and fourth ( 4 ) order moments of M, defined by equations (16) and (18), respectively. Figure 2 demonstrates that, to a very good accuracy, the various moments exhibit a power law: Similar observations were made for all the second and fourth moments, allowing us to define 2 S 2 , 4 S 4 , 4 2 S 2 and 4 SS in the same spirit as in equations (16), (18), (22) and (23). The dependence of the coefficients 2 2 and 4 4 shown in figure 2(c) is consistent with simple power laws: 2 ∼ C 2/3 η and 4 ∼ C 4/3 η . Note that the odd moments of M all vanish when α = 1. The observed dependence can be understood with the help of the following simplified argument. To analyze the 'first passage' problem that determines the statistical weight of each configuration, we return to equations (1) and (2). In the absence of nonlinearity and of noise (α = 1, C η = 0), the value of M(τ = 0), with |ρ| = ρ 0 is fixed. The time it takes for the value of ρ to reach |ρ| = L is T L ∼1/|M(τ = 0)|ln(L/ρ 0 ). In the presence of noise, C η = 0, M evolves only as a result of the noise term: The right-hand side of equation (24) is a Gaussian random noise term, with a zero mean, and a variance given by  The dependence on C η of the moments, see panel (c), suggests that the second moments behave as C 2/3 η and the fourth moments as C 4/3 η . Our numerical results for the second and fourth moments listed in section 4 indicate that the scalings are all consistent with the scaling relations documented here for vorticity.
The largest contribution to the noise is expected to come from early times (when |ρ| ∼ρ 0 is small), and the time T L has been simply estimated from the zero-noise case. As This estimate is consistent with the numerical observation that the moments of the order of 2, 2 , and 4, 4 , scale as ∼C 2/3 η ρ −4/3 , ∼C 4/3 η ρ −8/3 , completely consistent with the results of figure 2.
The argument above provides a simplified description of the growth of M, which does not fully take into account the coupling between ρ and M.
Our numerical results indicate that in the absence of the nonlinear term (α = 1), the third moments of M vanish, consistent with an underlying M → −M symmetry. This observation is consistent with the expectation that the lack of M → −M symmetry in the solution of the model is associated with the nonlinear term of the evolution equation. However, establishing directly the symmetry of the solutions from the equation of motion equations (1) and (2) or from the associated Fokker-Planck equation is not completely obvious. The evolution of the geometric variables ρ changes completely when the sign of M(0) in equation (2) is flipped. The symmetry is restored markedly when M(0) is sufficiently small, or more precisely, when ln(L/|ρ|)/|M(0)| is sufficiently large, so the noise term in equation (1) dominates. Still, our numerical results demonstrate that in the absence of nonlinearity (α = 1), the symmetry R → −R is very well observed.
We generally found that the value of C u does not affect the moments measured here very much.

Scaling law when α = 1
The nonlinear term in equation (7) is multiplied by (1 − α). Thus, varying α away from 1 effectively changes the strength of the nonlinearity in the problem.
The main result demonstrated here is that the nonlinear term does not modify the scaling dependence of the moments very much compared to the α = 1 case; see figure 3 in the case C u = 0.15 and C η = 0.5. The clear trend seen in figure 3 is a decrease of the moments when α increases. This trend has been observed for all the values of C u and C η we have investigated numerically (table 1).
Vortex stretching, i.e. the spontaneous amplification of vorticity by the flow, is an important feature of turbulent flows. In terms of the velocity gradient tensor, this property is expressed by the fact that the third invariant ω.s.ω is positive. In the model of [1], the quantity |ρ| 2 × tr(M 2 M T ) is directly related to the transfer of energy in the flow [1,33]: energy transfer from large to small scales is reflected by a positive value of − tr(M 2 M T ) . Note that an energy transfer independent of scale is equivalent to − tr(M 2 M T ) ∝ |ρ| 2 . For simplicity, we refer to − tr(M 2 M T ) in the following as the energy transfer term. Elementary manipulations show that In the case of the velocity gradient tensor, a consequence of homogeneity is that [34] tr(m 3 ) = 0, which implies that φ 3 s 3 = −3φ 3 ωsω . Thus, the production of vorticity (φ 3 ωsω > 0) is related to the fact that φ 3 s 3 is negative. The solutions of the model determined numerically here are such that 3 Another consequence of the fact that stretching is positive, and that the (random) matrices M are traceless and isotropic, is that the intermediate eigenvalue of the rate of strain should be mostly positive [34]. We have checked directly this property and compared the distribution of √ 6λ 2 / S 2 1/2 with the corresponding distribution in real flow. The distribution from the model  becomes more asymmetric when α diminishes. The predictions of the model are in qualitative agreement with DNS and experimental measurements [18].
As revealed by figure 3, the moments recorded here all follow a scaling law very close to the Kolmogorov K41 scaling: M p ∝ ( /|ρ| 2 ) p/3+δ p , with a value of δ p consistent with δ p = 0. Judging from our data, the potential departure from the Kolmogorov scaling law is obtained at α ≈ 0.3-0.5. The departure of the exponents from the Kolmogorov K41 laws did not exceed |δ 2 | 0.035, |δ 3 | 0.075 and |δ 4 | 0.10, which implies that the difference between the Kolmogorov K41 exponent and the measured exponent is less than ∼7.5%. We note that the values of the corrections to the K41 exponents, measured from the velocity structure functions, are also very weak (δ 2 ≈ 0.05 and δ 4 ≈ −0.05; [35]), compared to the values found here. A systematic investigation of the deviations from the K41 scaling behavior would require the investigation of higher-order moments of M, and over a larger range of scales than the range covered in this paper.
On general grounds, it is normal to expect that increasing the nonlinearity by decreasing α should lead to increased discrepancies with respect to the 'normal' scaling behavior, obtained when α = 1. It is a bit surprising to observe here that the deviations with respect to the K41  law are strongest for α = 0.3-0.5 and tend to decrease when α is further reduced towards 0. This effect is possibly related to an extra-conservation law, obtained when α = 0. In the α = 0 case, and in the purely deterministic case (C u = C η = 0), the quantity ρ i a · M ab · ρ i c · M cb (with summation over all repeated indices) is conserved. This invariant is the sum of the kinetic energies of all particles, with respect to the center of mass. This invariance is related to the fact that the energy transfer due to the nonlinearity in the model vanishes when α = 0 [1]. This conservation law seems to influence the solution in all sorts of ways, as demonstrated below.
Ignoring (weak) potential corrections to the normal scalings, the dependence of the moments is fully characterized by the numerical prefactors. The variation of this prefactor as a function of α is shown in figure 4 for the second (a), fourth (b) and third ((c) and (d)) moments (C u = 0.15 and C η = 0.5). Figure 4 reveals that the even moments are approximately constant at small values of α and decrease for α 0.3. We observe that the prefactor associated with the decay of the third moment vanishes when α → 1. In fact, figure 4(d) suggests that the third moment is proportional to (1 − α), for α 0.1. This observation is consistent with the notion that a nonvanishing value of the third order moments 3 is induced by the nonlinear term in equation (1), which is proportional to (1 − α). A full analytic understanding of the generation of skewness when (1 − α) > 0 and of the precise dependence of the third moments as a function of the parameter, in particular of (1 − α), requires a more complete analysis of the equation of motion.
In addition to a decrease of the moments as α increases, figure 4 reveals a precise ordering between the various moments. For small values of α (α 0.9), the strain variance 2 S 2 , defined by equation (16), is found to be larger than the vorticity variance 2 2 . Figure 4(a) indicates that vorticity prevails (slightly) when α → 1, for α 0.9. In the case of the fourth moments, our numerical results indicate that the fourth moment of strain 4 S 4 is larger than the fourth moment of vorticity 4 4 , at least for values of α 0.7, the trend being reversed at larger values of α. In fact, the relative decay of the strain as α approaches the value 1 is such that for α very close to 1, the fourth moment of the strain becomes comparable to the fourth-order mixed strain-vorticity moment 4 S 2 2 . The moment 4 SS is always found to be significantly smaller than the other moments. The precise ordering between the various moments, shown in figure 4 for the particular set of values C u = 0.15 and C η = 0.5, is in fact completely general. For all the values of the parameters we have studied (see table 1) the second moment tr(S 2 ) exceeds 2 for α 0.9, and the fourth moment tr(S 4 ) exceeds ( 2 ) 2 for α 0.7.
It is interesting to compare the ratios between the second-and fourth-order moments (see figures 4(a) and (b)) with the ratios corresponding to a Gaussian distribution of traceless matrices, or with the numbers determined in DNS. In the case of second moments, it is straightforward to estimate the ratio between strain and vorticity, as defined here, for the Gaussian ensemble of matrices: the quantities 2 S 2 and 2 2 are in the ratio 5 : 6. The model predicts that the vorticity invariant 2 2 is as expected in the Gaussian ensemble for α = 1 (in the absence of nonlinearity), but that the strain invariant 2 S 2 is comparatively larger as soon as α 0.9. Similarly, the fourth-order invariants 4 S 4 , 4 S 2 2 , 4 SS and 4 4 are, for a Gaussian ensemble of matrices, in the ratio 2.1 : 3 : 1 : 10 [24]. We find that for α = 1, in the absence of nonlinearity, the ratio between 4 S 4 and 4 SS is close to 3. This is consistent with the almost complete lack of alignment between vorticity and the eigenvectors of the strain matrix in the solutions of the model. However, 4 4 is smaller than 10 × 4 SS . On the other hand, 4 S 4 is larger than 2 × 4 S 2 2 . At smaller values of α, the fluctuations of tr(S 2 ) tend to become much more dominant compared to the fluctuations of vorticity 2 , consistent with the results for the second moments.
The main conclusion here is that the solutions of the model generally underpredict the importance of vorticity and overpredict the strain component. The fact that the two fourth-order moments 4 .S.S. and 4 S 2 2 are in the same ratio as in the case of Gaussian random matrices also indicates that the correlations between vorticity and the rate of strain eigenvectors are not properly captured by the model. These features are at odds with the results found in DNS for the velocity gradient, obtained at very high Reynolds numbers (R λ 10 3 ) [36], at least under the reasonable assumption that the trends observed at small scale in the inertial range properly represent the behavior obtained in the limit |ρ| → 0. This possibly suggests a weakness of the model.
The results obtained here should be contrasted with the semiclassical results. The parameter α was found in [25] to actually modify significantly the scaling exponent of tr(S 2 ) (which deviates significantly from −4/3 when α 0.35) and the sign of the energy transfer (which becomes negative when α 0.5). The full Monte-Carlo solutions do not confirm this picture: the signs of tr(M 2 M T ) always correspond to an energy transfer from large to small scales. The difference is possibly due to fluctuations not properly taken into account in the semiclassical approximation. It should also be stressed that the exact form of the model used in [25] is slightly different in the way the noise term in equation (2), responsible for relaxation towards isotropy of the spatial variables ρ, is implemented. It would be of interest to understand the origin of the difference between the two solutions.

Probability distribution function (PDF) of the (R, Q) invariants
This section is devoted to a discussion of the PDF of the invariants R and Q. We briefly recall here that the characteristic polynomial of the matrix M is As a result, the eigenvalues of the matrix M are fully determined by R and Q. The (R, Q) plane is separated by the zero-discriminant line 4Q 3 + 27R 2 = 0: above (respectively below) this line, M has two complex conjugate eigenvalues and one real one (respectively three real eigenvalues). The former case corresponds to a flow locally spiralling around an axis, whereas the latter situation corresponds to a purely straining flow. Thus, the PDF in the (R, Q) plane provides essential insight into the local structure of the flow [22]. The solutions all show a relatively mild change in shape compared to the symmetric, shell-shaped PDF obtained for a Gaussian velocity field [1]. In comparison with the solutions obtained in the semiclassical approximation [25], the shapes obtained with the full model show qualitatively a much closer similarity to the solutions obtained in experiments or in DNS [18].

Dependence of PDF on α
The solution for α = 0.9, corresponding to a very weak nonlinear term, is close to symmetric with respect to the vertical axis. The corresponding symmetry R → −R is observed when α = 1; it is the result of the underlying symmetry M → −M discussed above. As the value of α is decreased from α = 0.9 to α = 0.5, the nonlinearity increases, and the lack of symmetry R → −R becomes increasingly visible. One observes, in particular, a growth of the probability along the separatrix 4Q 3 + 27R 2 = 0 along the R 0 side. The comparatively large probability in this region of the (R, Q) plane has been reported many times in DNS, as well as in the semiclassical solution of the model [25]. In the latter case, this increased probability had been associated with the growth of the nonlinear term. In fact, the simplest (isotropic) closure approximation for the pressure Hessian in the evolution equation for A leads to a divergence in a finite time of the matrix A, with R → ∞, along the zero-discriminant line [3,22]. The increase in probability along this part of the (R, Q) plane is generally thought to be a dynamic consequence of this effect [15].
In this sense, it is a bit surprising to observe that when α decreases further (α 0.3), the tendency to accumulate high probability regions along the zero-discriminant line, in the R 0 part of the (R, Q) plane, saturates. Our numerical observation, however, is generally consistent with the results of figure 4, which also indicate that the moments stop to grow when α diminishes past ∼0.3. Figure 4 may even suggest a slight decrease of the moments when α decreases from 0.1 to 0, which is difficult to establish with accuracy in view of the large numerical error bars.
This observation is plausibly related to the peculiarities of the limit α → 0, induced by the conservation law of kinetic energy of the particles alluded to above when α = 0 and in the absence of noise terms. The invariants R and Q have been rescaled by tr(S 2 ) 3/2 and tr(S 2 ) , respectively. The evolution of the PDF is particularly visible in the development of the probability distribution along the zero-discriminant line 4Q 3 + 27R 2 = 0 along the positive R > 0 side. The isoprobability lines shown on the graphs correspond to probability levels 10 −n , where n is an integer (n 6).

Dependence of PDF on the amplitude of noise
From a qualitative point of view, the solutions of the model result from a competition between the purely stochastic terms (the right-hand side of equations (1) and (2)) and the nonlinear term. We consider here the influence of the noise term in the M-equation, equation (1), and measured directly by the parameter C η . The effect of noise on the structure of the PDF in the (R, Q) plane is shown in figure 6 at scale |ρ| = L/4, and for α = 0.5 and C u = 0.15. Increasing the parameter C η results in a strong reduction of the lack of symmetry R → −R. Indeed, the solution corresponding to C η = 0.25 exhibits a much more pronounced excess of probability along the zero-discriminant line on the R 0 side than the solution with C η = 2.
This observation can be qualitatively understood as follows. Increasing the stochastic term, while keeping all parameters fixed, tends to diminish the relative effect of the nonlinearity, responsible for the lack of M → −M symmetry. It is therefore reasonable to expect that increasing C η tends to make the solution more symmetric with respect to the vertical axis in the (R, Q) plane. . The lack of symmetry decreases more weakly when the noise amplitude is increased. The isoprobability contours are as in figure 5. Figure 7 shows the evolution through scale of the PDF in the (R, Q) plane for the solution of the model corresponding to α = 0.7, C u = 0.15 and C η = 1. The quality of the numerically determined PDF degrades slightly as the scale |ρ| decreases. This is a consequence of the statistical difficulties in determining precise solutions of the model. The quality of the solution is, however, sufficient to show that the evolution of the PDF as a function of scale is very limited. The same feature has been observed for all the values of the parameters that we have considered.

Scale dependence of PDF
The very weak dependence of the shape of the PDF in the (R, Q) plane as a function of scale |ρ| suggests that the flow properties are almost self-similar in the inertial range. This is generally consistent with the results found in real flows, which show that the deviations from the self-similar K41 picture are relatively weak. In this sense, the solutions of the model properly reproduce this feature of turbulent flows (see also [18]).

Density of strain, enstrophy, enstrophy production and energy transfer in the (R, Q) plane
Further information on the local structure of the flow in the (R, Q) plane can be obtained by representing the local density of dynamic quantities, such as strain, tr(S 2 ), enstrophy 2 , enstrophy production .S. and energy transfer −tr(M 2 M T ). For each of the above observables, the represented density is equal to the conditional PDF as a function of R and Q, multiplied by the PDF of (R, Q). (respectively negative, 0) contour lines are indicated by dashed (respectively dashed-dotted, full) lines. Enstrophy production is positive (respectively negative) for negative (respectively positive) values of R, the positive contribution being larger than the negative one, resulting globally in net enstrophy production, as expected in turbulent flows. In the same spirit, the energy dissipation term, −tr(M 2 M T ) (d), consists of a strong positive contribution along the separatrix in the (R, Q) plane, on the R 0 side, and a weaker negative contribution along the separatrix, on the R 0 side. As it was the case for the enstrophy production term, the positive contribution prevails over the negative one, thus resulting in a net energy transfer towards small scales, again, as expected in turbulent flows. The densities shown in figure 8 are qualitatively very similar to the densities obtained either from DNS or in real experiments [18]. This qualitative similarity, to be documented elsewhere, is encouraging as it demonstrates that the model captures important aspects of the statistics. Figure 8 corresponds to a tetrad size of |ρ| = L/4, and for a set of values α = 0.5, C u = 0.15 and C η = 0.5. The general aspect of the density plots does not change much when the size |ρ| decreases. As was the case for the PDF in the (R, Q) plane, changing α or C η results in quantitative changes of the density. Figure 9 illustrates the effect of changing systematically α in the case of the enstrophy production density. When α = 0.1 (figure 9(a)) the region in the (R, Q) plane where enstrophy production is positive has a much larger extension than the region corresponding to negative enstrophy production. In contrast, when α = 0.7 (figure 9(c)) the regions of positive and negative enstrophy production tend to become comparable in extension in the (R, Q) plane. Consistent with the observations of figures 4(c) and (d) and 5, diminishing the value of α results in destroying the R → −R symmetry, inherited from the M → −M symmetry when α = 1. Similarly, increasing the value of the noise term C η also results in restoring the symmetry R → −R (not shown).

Summary and discussions
The obtainment of reliable solutions of the tetrad model [1] is the main achievement of this paper. The ability to sample the very rare configurations that really contribute to the solution is the key ingredient, enabling the generation of reliable solutions for the model. We have used methods based on importance sampling strategies and, in particular, on the pruned-and-enlarged methods [28]. The results obtained here are a significant improvement on solutions determined with the help of semiclassical methods. Our numerical solutions allow us to compare quantitatively the predictions of the model with the available experimental and DNS results. We discuss in turn the success and the shortcomings of the model. The solutions obtained here generally show qualitative agreement with several of the main features observed in the Navier-Stokes turbulence. The low order moments of the solutions behave, for tetrad scale in the inertial range, as power laws, with scaling exponents very close to the Kolmogorov K41 range. Weak deviations with respect to the K41 scaling, consistent with the experimentally observed values for the second and fourth moments, are possibly obtained for values of the parameter α ≈ 0.3 − 0.5. More importantly, the model correctly reproduces the signs of all the terms associated with small-scale generation and energy transfer in turbulent flows. The predicted structures of the PDFs in the (R, Q) plane are also much closer to the shapes observed numerically and experimentally [18] than the semiclassical solutions. Lastly, the predicted density of enstrophy, strain, enstrophy production and energy transfer in the (R, Q) plane are in very good qualitative agreement with what is observed experimentally [18].
While the model properly reproduces the important properties of turbulent flows as listed above, it fails to capture some other aspects of turbulent flows. One of the main features of the model is that it tends to predict an excess of strain, compared to vorticity. The available numerical result [25], at moderate Reynolds numbers, suggests that vorticity is larger than strain. This trend is confirmed by recent experimental results [18]. The growth of the strain with respect to vorticity in our model is clearly induced by the nonlinear term: figure 4 demonstrates that even a weak nonlinearity (α 0.9) is sufficient to make strain prevalent compared to vorticity. Interestingly, this suggests that, while the nonlinearity in the model correctly produces the lack of R → −R symmetry also observed in DNS or experiments, this nonlinearity introduces too strong a bias towards strong strain.
Related to the observation that the model does not predict the proper balance between strain and vorticity, we also observed that the ratios between the fourth order moments of the coarsegrained velocity gradient tensor obtained from the model seem to be at odds with the properties obtained for the velocity gradient tensor obtained in DNS at very high Reynolds numbers. Additionally, the model does not predict the relationship between vorticity and eigenvectors of strain observed numerically and in experiments: the tendency to align with the intermediate eigenvector of S is much weaker in the model than in reality [18,37].
Thus, while the solutions of the model correctly reproduce many important features of turbulent flows, the model misses several other properties. This situation points to some limitations of the model, as formulated here, that need to be identified and, if possible, fixed. In this situation, it is worth stressing that the possibility of formulating a model of turbulence in terms of stochastic differential equations and solving it accurately with the help of limited computer resources (a workstation) provides a way of investigating systematically the essential ingredients leading to the generation of small scales in turbulent flows, the understanding of which is crucial for many applications of turbulent flows.