Statistical Lyapunov Theory Based on Bifurcation Analysis of Energy Cascade in Isotropic Homogeneous Turbulence: A Physical–Mathematical Review

This work presents a review of previous articles dealing with an original turbulence theory proposed by the author and provides new theoretical insights into some related issues. The new theoretical procedures and methodological approaches confirm and corroborate the previous results. These articles study the regime of homogeneous isotropic turbulence for incompressible fluids and propose theoretical approaches based on a specific Lyapunov theory for determining the closures of the von Kármán–Howarth and Corrsin equations and the statistics of velocity and temperature difference. While numerous works are present in the literature which concern the closures of the autocorrelation equations in the Fourier domain (i.e., Lin equation closure), few articles deal with the closures of the autocorrelation equations in the physical space. These latter, being based on the eddy–viscosity concept, describe diffusive closure models. On the other hand, the proposed Lyapunov theory leads to nondiffusive closures based on the property that, in turbulence, contiguous fluid particles trajectories continuously diverge. Therefore, the main motivation of this review is to present a theoretical formulation which does not adopt the eddy–viscosity paradigm and summarizes the results of the previous works. Next, this analysis assumes that the current fluid placements, together with velocity and temperature fields, are fluid state variables. This leads to the closures of the autocorrelation equations and helps to interpret the mechanism of energy cascade as due to the continuous divergence of the contiguous trajectories. Furthermore, novel theoretical issues are here presented among which we can mention the following ones. The bifurcation rate of the velocity gradient, calculated along fluid particles trajectories, is shown to be much larger than the corresponding maximal Lyapunov exponent. On that basis, an interpretation of the energy cascade phenomenon is given and the statistics of finite time Lyapunov exponent of the velocity gradient is shown to be represented by normal distribution functions. Next, the self–similarity produced by the proposed closures is analyzed and a proper bifurcation analysis of the closed von Kármán–Howarth equation is performed. This latter investigates the route from developed turbulence toward the non–chaotic regimes, leading to an estimate of the critical Taylor scale Reynolds number. A proper statistical decomposition based on extended distribution functions and on the Navier–Stokes equations is presented, which leads to the statistics of velocity and temperature difference.


Introduction
This article presents a review of previous works of the author regarding an original Lyapunov analysis of the developed turbulence which leads to the closures of the von Kármán-Howarth and Corrsin equations and to the statistics of both velocity and temperature difference [1][2][3][4][5][6][7]. This theory studies the fully developed homogeneous isotropic turbulence through the bifurcations of the incompressible Navier-Stokes equations using a specific statistical Lyapunov analysis of the fluid kinematic field. In addition, now it is introduced the energy cascade interpretation and explained some of the mathematical properties of the proposed closures. This work is organized into two parts. One is the reasoned review of previous results but with new demonstrations and theoretical procedures. The other one, presented in sections marked with asterisk symbol "*", concerns new theoretical issues of the proposed turbulence theory.
Although numerous articles were written which concern the closures of the Lin equation in the Fourier domain [8][9][10][11][12][13][14][15][16], few works address the closures of the autocorrelation equations in the physical space. These last ones, being based on the eddy-viscosity concept, describe diffusive closure models. Unlike the latter, the proposed Lyapunov theory provides nondiffusive closures in the physical space based on the property that, in developed turbulence, contiguous fluid particles trajectories continuously diverge. Thus, the main purpose of this review is to summarize the results of the previous works based on a theory which does not use the eddy-viscosity paradigm and to give new theoretical insights into some related issues.
The homogeneous isotropic turbulence is an ideal flow regime characterized by the energy cascade phenomenon where the diverse parts of fluid exhibit the same statistics and isotropy. On the other hand, the turbulent flows occurring in nature and in the various fields of engineering are generally much more complex than homogeneous isotropic turbulence. In such flows, spatial variations of average velocity and of other statistical flow properties can happen causing very complex simultaneous effects that add to the turbulent energy cascade and interact with the latter in a nontrivial fashion. Hence, the study of the energy cascade separately from the other phenomena requires the analysis of isotropic homogenous turbulence.
The von Kármán-Howarth and Corrsin equations are the evolution equations of longitudinal velocity and temperature correlations in homogeneous isotropic turbulence, respectively. Both the equations, being unclosed, need the adoption of proper closures [17][18][19][20]. In detail, the von Kármán-Howarth equation includes K, the term due to the inertia forces and directly related to the longitudinal triple velocity correlation k, which has to be properly modelled. The modeling of such term must take into account that, due to the inertia forces, K does not modify the kinetic energy and satisfies the detailed conservation of energy [18]. This latter states that the exchange of energy between wave-numbers is only linked to the amplitudes of such wave-numbers and of their difference [21]. Different works propose for the von Kármán-Howarth equation the diffusion approximation [22][23][24] where r and D = D(r) are separation distance and turbulent diffusion parameter, respectively and u 2 = u i u i /3 corresponds to the longitudinal velocity standard deviation. Following Equation (1), the turbulence can be viewed as a diffusivity phenomenon depending upon r, where K will include a term proportional to ∂ 2 f /∂r 2 . In the framework of Equation (1), Hasselmann [22] proposed, in 1958, a closure suggesting a link between k and f which expresses k in function of the momentum convected through a spherical surface. His model, which incorporates a free parameter, expresses D(r) by means of a complex expression. Thereafter, Millionshtchikov developed a closure of the form D(r) = k 1 ur, where k 1 represents an empirical constant [23]. Although both the models describe two possible mechanisms of energy cascade, in general, do not satisfy some physical conditions. For instance, the Hasselmann model does not verify the continuity equation for all the initial conditions, whereas the Millionshtchikov equation gives, following Equation (1), values of velocity difference skewness in contrast with experiments and energy cascade [18]. More recently, Oberlack and Peters [24] suggested a closure where D(r) = k 2 ru 1 − f , being k 2 a constant parameter. The authors show that such closure reproduces the energy cascade and, for a proper choice of k 2 , provides results in agreement with the experiments [24].
For what concerns the Corrsin equation, this exhibits G, the term responsible for the thermal energy cascade. This quantity, directly related to the triple velocity-temperature correlation m * , also needs adequate modellation. As G depends also on the velocity correlation, the Corrsin equation requires the knowledge of f , thus it must be solved together to the von Kármán-Howarth equation. Different works can be found in the literature which deal with the closure of Corrsin equation. Some of them study the self-similarity of the temperature correlation in order to analyze properties and possible expressions for G. Such studies are supported by the idea that the simultaneous effect of energy cascade, conductivity and viscosity, makes the temperature correlation similar in the time. This question was theoretically addressed by George (see [25,26] and references therein) which showed that the decaying isotropic turbulence reaches the self-similarity, while the temperature correlation is scaled by the Taylor microscale whose current value depends on the initial condition. More recently, Antonia et al. [27] studied the temperature structure functions in decaying homogeneous isotropic turbulence and found that the standard deviation of the temperature, as well as the turbulent kinetic energy, follows approximately the similarity over a wide range of length scales. There, the authors used this approximate similarity to estimate the third-order correlations and found satisfactory agreement between measured and calculated functions. On the other hand, the temperature correlation can be obtained using proper closures of von Kármán-Howarth and Corrsin equations suitable for the energy cascade phenomenon. On this argument, several articles has been written. For instance, Baev and Chernykh [28] (and references therein) analyzed velocity and temperature correlations by means of a closure model based on the gradient hypothesis which relates pair longitudinal second and third order correlations, by means of empirical coefficients.
Although other works regarding the von Kármán-Howarth equation were written [29][30][31][32][33], to the author's knowledge a physical-mathematical analysis based on basic principles which provides analytical closures of von Kármán-Howarth and Corrsin equations has not received due attention. Therefore, the aim of the this work is to present a review of the Lyapunov analysis presented in [1][2][3][4][5][6][7] and new theoretical insights into some related issues.
In the present formulation, based on the Navier-Stokes bifurcations, the current fluid placements, together with velocity and temperature fields, are considered to be fluid state variables. This leads to the closures of the autocorrelation equations and helps to interpret the mechanism of energy cascade as due to the continuous divergence of the contiguous trajectories.
In line with Ref. [3], the present work first addresses the problem for defining the bifurcations for incompressible Navier-Stokes equations, considering that these latter can be reduced to an opportune symbolic form of operators for which the classical bifurcation theory of differential equations can be applied [34]. In such framework, this analysis remarks that a single Navier-Stokes bifurcation will generate a doubling of the velocity field and of all its several properties, with particular reference to the characteristic length scales. If on one side the lengths are doubled due to bifurcations, on the other hand the characteristic scale for homogeneous flows in infinite domains is not defined. Hence, the problem to define the characteristic length-and therefore the flow Reynolds number-in such situation is also discussed. Such characteristic scale is here defined in terms of spatial variations of initial or current velocity field in such a way that, in fully developed homogeneous isotropic turbulence, this length coincides with the Taylor microscale. As far as the characteristic velocity is concerned, this is also defined in terms of velocity field so that, in developed turbulence, identifies the velocity standard deviation.
The trajectories bifurcations in the phase space of the velocity field are here formally dealt with using a proper Volterra integral formulation of the Navier-Stokes equations, whereas the turbulence transition is qualitatively analyzed through general properties of the bifurcations and of the route toward the fully developed chaos. This background, regarding the general bifurcations properties and the route toward the chaos, will be useful for this analysis.
The adopted statistical Lyapunov theory shows how the fluid relative kinematics can be much more rapid than velocity and temperature fields in developed turbulence, so that fluid strain and velocity fields are statistically independent with each other. Moreover, in addition to References [1][2][3][4][5][6][7], this analysis introduces the bifurcation rate of the velocity gradient, a quantity providing the frequency at which the velocity gradient determinant vanishes along fluid particles trajectories. The bifurcation rate, in fully developed turbulence, is shown to be much greater than the maximal Lyapunov exponent of the velocity gradient. This explains the energy cascade through the relation between material vorticity, Lyapunov vectors and bifurcation rate using the Lyapunov theory. In detail, the energy cascade can be viewed as a continuous and intensive stretching and folding process of fluid particles which involves smaller and smaller length scales during the fluid motion, where the folding frequency equals the bifurcation rate.
Next, the statistics of the Lyapunov exponents is reviewed. In agreement with Reference [6], we show that the local Lyapunov exponents are uniformely unsymmetrically distributed in their interval of variation. Unlike Reference [6] which uses the criterion of maximum entropy associated with the fluid particles placements, the isotropy and homogeneity hypotheses are here adopted. A further result with respect to the previous issues pertains the finite time Lyapunov exponents statistics: through the bifurcation analysis and the central limit theorem, we show that the finite time Lyapunov exponent tends to a fluctuating variable distributed following a normal distribution function.
Thereafter, the closure formulas of von Kármán-Howarth and Corrsin equations are derived through the Liouville equation and finite scale Lyapunov exponent statistics. These closures do not correspond to a diffusive model, being the result of the trajectories' divergence in the continuum fluid. Such formulas coincide with those just obtained in References [1,4,5] where it is shown that such closures adequately describe the energy cascade phenomenon, reproducing, negative skewness of velocity difference, the Kolmogorov law and temperature spectra in line with the theoretical argumentation of Kolmogorov, Obukhov-Corrsin and Batchelor [35][36][37], with experimental results [38,39] and with numerical data [40,41]. These closures are here achieved by using different mathematical procedures with respect to the other articles [1,4,5]. While the previous works derive such closures studying the local fluid act of motion in the finite scale Lyapunov basis [1,4] and adopting maximum and average finite scale Lyapunov exponents [5], here these closures are obtained by means of the local finite scale Lyapunov exponents PDF, showing that the assumptions of References [1,4,5] agree with this analysis, corroborating the previous results. Some of the properties of the proposed closures are then studied, with particular reference to the evolution times of the developed correlations and their self-similarity. In detail, as new result with respect the previous articles, this analysis shows that the proposed closures generate correlations self-similarity in proper ranges of separation distance, which is directly linked to the particles trajectories divergence.
Furthermore, a novel bifurcation analysis of the closed von Kármán-Howarth equation is proposed, which considers the route starting from the fully developed turbulence toward the non-chaotic regimes. This extends the discussion of the previous works and represents an alternative point of view for studying the turbulent transition. According to this analysis, the closed von Kármán-Howarth equation is decomposed in several ordinary differential equations through the Taylor series expansion of the longitudinal velocity correlation. This procedure, which also accounts for the aforementioned self-similarity, leads to estimating the Taylor scale Reynolds number at the transition. This latter is found to be 10, a value in good agreement with several experiments which give values around 10, and in particular with the bifurcations analysis of the energy cascade of Reference [3], which provides a critical Reynolds number of 10.13 if the route toward the turbulence follows the Feigenbaum scenario [42,43].
Finally, the statistics of velocity and temperature difference, of paramount importance for estimating the energy cascade, is reviewed. While References [1,2,4,7] determine such statistics through a concise heuristic method, this analysis uses a specific statistical decomposition of velocity and temperature which adopts appropriate stochastic variables related to the Navier-Stokes bifurcations. The novelty of the present approach with respect to the previous articles is that the random variables of such decomposition are opportunely chosen to reproduce the Navier-Stokes bifurcation effects and the isotropy: these are highly nonsymmetrically distributed stochastic variables following opportune extended distribution functions which can assume negative values. Such decomposition, able to reproduce negative skewness of longitudinal velocity difference, provides a statistics of both velocity and temperature difference in agreement with theoretical and experimental data known from the literature [44][45][46][47][48]. Here, in addition to References [1,2,4,7], a detailed mathematical analysis is presented which concerns the statistical properties of the aforementioned extended distribution functions in relation to the Navier-Stokes bifurcations.
In brief, the original contributions of the present work can be summarized as: (i) The bifurcation rate associated with the velocity gradient is shown to be much larger than the maximal Lyapunov exponent of the velocity gradient. (ii) As the consequence of (i), the energy cascade can be viewed as a succession of stretching and folding of fluid particles which involves smaller and smaller length scales, where the particle folding happens at the frequency of the bifurcation rate. (iii) As the consequence of (i), the central limit theorem provides reasonable argumentation that the finite time Lyapunov exponent is distributed following a gaussian distribution function. (iv) The proposed closures generate correlations self-similarity in proper ranges of variation of the separation distance which is directly caused by the continuous fluid particles trajectories divergence. (v) A specific bifurcation analysis of the closed von Kármán-Howarth equation is proposed which allows to estimate the critical Taylor scale Reynolds number in isotropic turbulence. (vi) A statistical decomposition of velocity and temperature is presented which is based on stochastic variables distributed following extended distribution functions. Such decomposition leads to the statistics of velocity and temperature difference, where the intermittency of these latter increases as Reynolds number and Péclet number rise.

Background
In the framework of the link between bifurcations and turbulence, this section deals with some of the fundamental elements of the Navier-Stokes equations and heat equation, useful for the present analysis. In particular, we will address the problem of defining an adequate bifurcation analysis for the Navier-Stokes equations and will analyze the meaning of the characteristic length scales when a homogeneous flow is in an infinite domain. All the considerations regarding the fluid temperature can also be applied to any passive scalar that exhibits diffusivity. A statistically homogeneous and isotropic flow with null average velocity is considered.
In order to formulate the bifurcation analysis, we start from the Navier-Stokes equations and the temperature equation where u = u(t, x), p = p(t, x) and ϑ = ϑ(t, x) are velocity, pressure and temperature fields, ν and χ = kρ/C p are fluid kinematic viscosity and thermal diffusivity, being ρ = const, k and C p density, fluid thermal conductivity and specific heat at constant pressure, respectively. In this study ν and χ are supposed to be independent from the temperature, thus Equation (2) is autonomous with respect to Equation (3), whereas Equation (3) will depend on Equation (2). To define the bifurcations of Equations (2) and (3), such equations are first expressed in the symbolic form of operators. To this end, in the momentum Navier-Stokes equations, the pressure field is eliminated by means of the continuity equation, thus Equations (2) and (3) are formally written aṡ in which N is a nonlinear quadratic operator incorporating −u · ∇ x u, −ν∇ 2 x u and the integral nonlinear operator which expresses the pressure gradient as a functional of the velocity field, being Therefore, p provides nonlocal effects of the velocity field [49] and the Navier-Stokes equations are reduced to be an integro-differential equation formally expressed by Equation (4). For what concerns Equation (5), it is the evolution equation of ϑ, where M is a linear operator of ϑ. Accordingly, transition and turbulence are caused by the bifurcations of Equation (4), where ν −1 plays the role of the control parameter. At this stage of the analysis, it is worth to remark the following two items: (a) there is no explicit methods of bifurcation analysis for integro-differential equations such as Equation (4). (b) since the flow is statistically homogeneous in an infinite domain, characteristic scales of the problem are not defined.
The item (a) can be solved according to the analysis method proposed by Ruelle and Takens in Reference [34]: it is supposed that the infinite dimensional space of velocity field {u} can be replaced by a finite-dimensional manifold, then Equation (4) can be reduced to be the equation of the kind studied by Ruelle and Takens in Reference [34]. Therefore, the classical bifurcation theory of ordinary differential equations [34,43,50] can be formally applied to Equation (4) and the present analysis can be considered valid within the limits of the formulation proposed in Reference [34].
For what concerns the characteristic length, a homogeneous flow in infinite domain is free from boundary conditions, thus the characteristic scale, being not defined, is here chosen in function of the spatial variations of the current velocity field. Thus, for all flow regimes in infinite regions, (i.e., non-chaotic, turbulent and transition flows), characteristic length and velocity, L and U respectively, are here chosen in terms of volume integrals of u in the following manner where V is the fluid domain volume, ":" denotes the Frobenius inner product and c = O(1) is a dimensionless constant which will be properly chosen. The flow Reynolds number is then defined in terms of U and L as Equation (8) provides an extension of the Taylor scale Reynolds number which applies for every flow regime. In particular, such definition holds also for non turbulent flows, where U and L, although not velocity standard deviation and statistical correlation scale, provide a generalization of the latter. In fully developed homogeneous turbulence, the volume integrals appearing in Equation (7) equal statistical averages calculated over the velocity field ensemble, such as velocity standard deviation and dissipation rate. Accordingly, in isotropic homogeneous turbulence, L and U identify, respectively, the Taylor scale λ T and standard deviation u of one of the velocity components and Re = U L/ν coincides with the Taylor scale Reynolds number R T . Such definitions (7) extend the concept of velocity variance and statistical correlation scale and will be used for the bifurcation analysis proposed in this work.

Navier-Stokes Bifurcations
Before introducing the bifurcations analysis of the Navier-Stokes equations in the operatorial form (4), it is worth remarking that a given point in the space of velocity fields setū ∈ {u}-or temperature fieldθ ∈ {ϑ}-corresponds to a spatial distribution including all its characteristics, in particular the length scales associated withū.
The bifurcations of Equation (4) happen when the Jacobian ∇ u N exhibits at least an eigenvalue with zero real part (NS-bifurcations), and this occurs when Such bifurcations are responsible for multiple velocity fieldsû which provides the same fieldu. In fact, during the fluid motion, multiple solutionsû andθ can be determined, at each instant, through inversion of Equation (4)u = N(u; ν) In the framework of the trajectories bifurcations in the phase space, the fluid motion can be expressed by means of Equation (4) and initial conditions u(0) and ϑ(0), using the following Volterra integral formulation where N and M are proper operators such that (12) Specifically, N is a nonlinear operator of u, where the image N ({u}) has the same structure of {u}, whereas M is linear with respect to ϑ and the image M({u} × {ϑ}) is isomorphic with {ϑ}. Thus, N and M admit in general the following jacobians According to Equation (11), a trajectory bifurcation happens when ∇ u N is singular, that is when det (∇ u N ) = 0 (14) and the multiple solutions of Equation (11), sayû andθ, are given in terms of u and ϑ through using the implicit functions theorem. Therefore, if velocity and temperature fields are supposed to be known for ν=ν 0 , the fields calculated for ν =ν 0 are formally expressed aŝ

Qualitative Analysis of the Route Toward the Chaos
With reference to Equation (10) or (16), when ν −1 is relatively small, N and N behave like linear operators and Equation (15) returnsû ≡ u(t, x) as unique solution. Increasing ν −1 , the Navier-Stokes equations encounter the first bifurcation at ν = ν 1 , the jacobian ∇ u N is singular there, and thereafter Equation (15) determines different velocity fieldsû with the corresponding length scales. A single bifurcation causes a doubling of u, that is, a doubling of the velocity values and of the length scales. Although the route toward the chaos can be of different kinds [34,42,43,51], one common element of these latter is that the number of encountered bifurcations at the onset of the chaotic regimes is about greater than three. Hence, if ν −1 is quite small, the velocity field can be represented by its Fourier series of a given basic scale. The first bifurcation introduces new solutionsû whose Fourier characteristic lengths are independent from the previous one. Thereafter, each bifurcation adds new independent scales, and, after the third bifurcation (ν −1 = ν −1 * ), the transition occurs, the several characteristic lengths and the velocity values appear to be continuously distributed and thus the velocity field is represented by the Fourier transform there. In such situations, a huge number of such solutions are unstable; u(t, x) tends to sweep the entire velocity field set and the motion is expected to be chaotic with a high level of mixing. As forθ, M and M are both linear operators of ϑ, thusθ follows the variations ofû.
If ν −1 does not exceed its critical value, say ν −1 * , the velocity fields satisfying Equation (15) are limited in number and this corresponds to the intermediate stages of the route toward the chaos. On the contrary, when ν −1 > ν −1 * , the region of developed turbulence where λ NS > 0 is observed, being λ NS the average maximal Lyapunov exponent of the Navier-Stokes equations, is formally calculated as and y is the Lyapunov vector associated with the Navier-Stokes equations. Then, ν −1 * depends on u, and Re * , calculated with Equation (7), can be roughly estimated as the minimum value of Re for which λ NS ≥ 0. Figure 1 qualitatively shows the route from non-chaotic regimes toward the developed turbulence. Specifically, Figure 1a,b report two bifurcation maps at a given instant, providing the velocity component u 1 in a point of the space and one characteristic scale of the velocity field in function of ν −1 . Figure 1c-e symbolically represent, for assigned values of ν, the velocity field set (points inside the dashed circle), three different solutions of the Navier-Stokes equations-say P, Q and R-and the several subsets σ 1 , σ 2 ,... which correspond to islands that are not swept during the fluid motion. The figure also depicts L=L(ν −1 ) and U=U(ν −1 ) (Figure 1f,g), formally calculated with Equation (7). Following Equation (16), these maps are not universal, as u 1 = u 1 (ν −1 ), = (ν −1 ), L = L(ν −1 ) and U = U(ν −1 ) do not represent universal laws and their order of magnitude will depend on velocity field at ν −1 0 . When ν −1 > ν −1 3 , the number of solutions diverges and the bifurcation tree of u 1 and drastically changes its structure showing tongue geometries that develop from the different bifurcations. As long as ν −1 does not exceed much ν −1 3 , the extension of such tongues is relatively bounded, whereas the measure of the islands σ k is quite large. This means that, although u 1 and exhibit chaotic behavior there, these do not sweep completely their variation interval, thus Equation (4) do not behave like an ergodic dynamic system there. This corresponds to Figure 1c, where the velocity fields P, Q and R, being differently placed with respect to σ k , k = 1, 2, .. will exhibit different values of average kinetic energy and dissipation rate in V. As ν −1 rises, these tongues gradually increase their extension whereas the measures of σ k diminish (see Figure 1d) until reaching a situation in which the bifurcation tongues overlap with each other and the islands σ k vanish ( Figure 1e). Such developed overlapping corresponds to the chaotic behavior of u 1 and , where these latter almost entirely describe their variation interval: Equation (4) behave like an ergodic dynamic system there, whereas all the velocity fields, in particular P, Q, and R, although different to each other, give the same values of average kinetic energy and dissipation rate in V. This is the onset of the fully developed turbulence.
As far as L and U are concerned, these are both functionals of u following Equation (7), accordingly their variations in terms of ν −1 are peculiar, with quite different results with respect to u 1 and . In particular, the structure of the first three bifurcations do not show important differences with respect to u 1 and , whereas, after the third bifurcation (ν −1 > ν −1 3 ), the chaotic regime begins and the bifurcation tree of U and L exhibits a completely different shape to the corresponding zone of u 1 and . In detail, the chaotic region extension of U and L appears to be more limited than that of u 1 and until to collaps in the lines A-B when ν −1 > ν −1 A . This is because the several bifurcations in correspond to a large number of solutions that show different levels of average kinetic energy and dissipations rate in V which are in some way comparable to each other, respectively. Hence, although the chaotic regime is characterized by myriad of values of u 1 and which widely sweep the corresponding ranges, L and U, being related to average kinetic energy and dissipation rate, will exhibit smaller variations. For relatively high values of ν −1 , when the velocity fluctuations behavior is ergodic, the averages calculated on phase trajectory tends to the spatial averages. The region of chaotic regime collaps into the line A-B there. Along such lines, for assigned ν, all the solutions-in particular P, Q and R-will exhibit the same level of kinetic energy and dissipation and this represents the regime of fully developed turbulence. The Reynolds number Re = ν −1 U L is shown in terms of ν −1 in Figure 2. Also this map is non universal as it depends on ν −1 0 . Nevertheless, such representation allows to identify the critical Reynolds number Re * = R T * = ν −1 * U * L * , the minimum value of R T for which the flow maintains statistically homogeneous and isotropic compatible with λ NS ≥ 0. Hence, a critical Reynolds number Re * = R T * will assume a unique value, represented by the point A of Figures 1 and 2, which plays the role of an universal limit in homogeneous isotropic turbulence. Then, ν * ≡ ν A , L * ≡ L A , U * = U A and the lines A-B represent regimes of fully developed homogeneous isotropic turbulence where We conclude this section by remarking that the characteristic length of the problem is an undefined quantity in infinite domain. Therefore, the length scales of u are used for determining the flow Reynolds number the critical value of which, Re * = ν −1 * U * L * has to be properly estimated. Accordingly, L * ≡ λ T * and U * ≡ u * , linked with each other, will depend on R T * and ν.
Such qualitative analysis is here used as background to formulate a specific bifurcation analysis of the velocity correlation equation and to determine an estimate of the critical Reynolds number R T * .

Kinematic Bifurcations. Bifurcation Rate
The Navier-Stokes bifurcations have significant implications for what concerns the relative kinematics of velocity field. This kinematics is described by the separation vector ξ (finite scale Lyapunov vector), which satisfies the following equationṡ being x(t) and y(t) = x(t) + ξ(t) two fluid particles trajectories. In the case of contiguous trajectories, |ξ| → 0, and Equation (19) where dx and ∇ x u(t, x) are, respectively, elemental separation vector and velocity gradient. One point of the physical space is of bifurcation for the velocity field (kinematic bifurcation) if ∇ x u(t, x) has at least an eigenvalue with zero real part and this happens when its determinant vanishes, that is, As seen, when R T > R T * , due to Navier-Stokes bifurcations, the velocity field evolution will be characterized by continuous distributions of length scales and velocity values. Therefore, for t > 0, the velocity gradient field will exhibit nonsmooth spatial variations where ∇ x u(t, x) = 0, and its determinant, det (∇ x u(t, x)), is expected to frequently vanish along fluid particles trajectories.
To justify this, one could search a link between such property and the statistics of the eigenvalues of ∇ x u which directly arises from the fluid incompressibility [52]. In this regard, observe that an arbitrary particle trajectory l t : x(t) belongs to the surface Σ 1 and identically satisfies the equation Thanks to Navier-Stokes bifurcations and fully developed turbulence hypothesis, for t > 0, Σ 1 and l t will show abrupt variations in their local placement, orientation and curvatures and will tend to sweep the entire physical space. On the other hand, the vanishing condition of velocity gradient determinant defines the surface Σ 2 = Σ 1 . Thus, the points which satisfy both the conditions (22) and (24) belong to the line l b = Σ 1 ∩ Σ 2 , and represent all the possible kinematic bifurcations which could happen along l t . Because of fully developed turbulence, l b will also show nonsmooth spatial variations and will tend to describe the entire physical space. Therefore, the kinematic bifurcations that occur along l t are obtained as l t ∩ l b , being l t , l b ∈ Σ 1 . As l t and l b are two different curves of the same surface Σ 1 that exhibit chaotic behaviors, their intersections are expected to be very frequent, forming a highly numerous set of points on Σ 1 according to the qualitative scheme of Figure 3 wherein l t and l b are represented by solid and dashed lines. Specifically, for R T > R T * , t > 0, the Navier-Stokes bifurcations produce the regime of fully developed turbulence, where length scales and velocity values are continuously doubled and this causes situations where the number of the intersections between l t and l b (kinematic bifurcations) diverges. To show this, the kinematic bifurcation rate is now introduced. This quantity, calculated along a fluid particle trajectory, is defined as follows: The rate S b can be much greater than the eigenvalues modulus of ∇ x u and than its maximal Lyapunov exponent. In fact, due to the Navier-Stokes bifurcations and to the hypothesis of fully developed chaos, the characteristic scales of u are continuously doubled, thus D ≡ det (∇ x u) is expected to be a function of the kind det (∇ x u) = D (y 1 , y 2 , ..., y n ) , where, due to bifurcations, n tends to diverge and For one assigned velocity field, from Equations (25) and (26), the simultaneous values of u and ∇ x (det(∇ x u)) can cause very frequent kinematic bifurcations whose rate can be significantly greater than the maximal Lyapunov exponent of Equation (19) or (20). In fact, following Equations (25) and (26), the order of magnitude of S b identifies the ratio (large scale velocity)-(small scale length) where the small scale n represents the minimum distance between to successive kinematic bifurcations encountered along fluid particle trajectory. This means that the changing rate of ∇ x u along l t can be much more rapid than the rate of divergence of two contiguous trajectories. At this stage of the present study, S b is assumed to be much greater than the maximal Lyapunov exponent of Equation (19) and its estimation will be performed in the following as soon as n is identified by means of this analysis.

Lyapunov Kinematic Analysis
The aim of this section is to discuss how, in fully developed turbulence, the fluctuations of fluid particles displacements and local strain can be much more rapid and statistically independent with respect to the time variations of velocity field. To analyze this, consider that, in fully developed turbulence, the Navier-Stokes bifurcations cause non smooth spatial variations of u(t, x) which in turn deternine very frequent kinematic bifurcations. Due to the fluid incompressibility, two fluid particles will describe chaotic trajectories, x(t) and y(t) = x(t) + ξ(t), which diverge with each other with a local rate of divergence quantified by the local Lyapunov exponent of finite scale ξ According to such definition ofλ, around to a given instant, t 0 , ξ andξ can be expressed as as long as |ξ| ≈ |ξ(t 0 )| = r, where Q is an orthogonal matrix giving the orientation of ξ with respect to the inertial frame R and ω E is the angular velocity of ξ with respect to R whose determination is carried out by means of a proper orthogonalization procedure of the Lyapunov vectors described in Reference [6]. The classical local Lyapunov exponent is obtained for |ξ| → 0,λ → Λ, that is On the other hand, dx can be expressed through Equation (20) as follows where the exponential denotes the series expansion of operators Although in developed turbulence the Navier-Stokes bifurcations cause abrupt spatial variations of velocity and temperature, with λ NS > 0, due to fluid dissipation, u and ϑ are in any case functions of slow growth of t ∈ (0, ∞), whereas ξ and dx, being not bounded by the dissipation effects, are functions of exponential growth of t. Therefore, in line with the analysis of Reference [3], and taking into account that S b >> sup λ , that ξ and dx are much more rapid than u(t, x) being sup λ >> λ NS , it follows that ξ and dx will exhibit power spectra in frequency intervals which are completely separated with respect to those of the power spectum of u. To study this, consider now the Taylor series expansion of u with respect to t of the trajectories equations, that is, x = u(0, x(t)) + ..., ξ = u(0, x(t) + ξ(t)) − u(0, x(t)) + ..., for finite scale |ξ|, dẋ = ∇ x u(0, x(t))dx + ..., for contiguous trajectories (34) The first terms (terms of 0 order) of such Taylor series do not correspond to time variations in velocity field, thus these do not modify the fluid kinetic energy. Furthermore, as sup λ >> λ NS (fully developed turbulence), such terms reproduce the particles trajectories as long as Following Equation (35), the fluctuations of ξ and dx are statistically independent with respect to the time variations of the velocity field. Next, sup λ >> λ NS , thus the number of kinematic bifurcation, which happen for 0 < t < O(1/λ NS ), is expected to be quite high and can be considered to be significative from the statistical point of view. Now, according to the mathematical analysis of the continuum media [53], the following map is considered which expresses the placement of material elements at the current time t in function of their referential position, say x 0 = x(0) [53]. From Equation (32), the local fluid strain ∂x(t)/∂x 0 is then an exponential growth function of t which, thanks to the above mentioned property of independence of dx from u(t, x), results to be independent and much faster with respect to the time variations of the velocity field. In fact, from the Lyapunov theory of kinematic field, such strain reads as where G is a proper fluctuating matrix whose elements G ij = O(1) are functions of of slow growth of t.
As long as t ∈ (0, a) we have that is ∂x(t)/∂x 0 is independent of the time variations of the velocity field. In brief, as sup λ >> λ NS , two time scales are here considered: one associated with the velocity field and the other one related to the relative fluid kinematics. Thus, ξ, ∂x(t)/∂x 0 andλ are statistically independent of u. Furthermore, due to very frequent kinematic bifurcations in (t, t + 1/λ NS ), ξ, local strain andλ are expected to be continuously distributed in their variation ranges. This conclusion is supported by the arguments in References [54,55] (and references therein), where the author remarks among other things that the fields u(t, x), (and therefore also u(t, x + ξ) − u(t, x)) produce chaotic trajectories also for relatively simple mathematical structure of u(t, x) (also for steady fields!).

*Turbulent Energy Cascade, Material Vorticity and Link with Classical Kinematic Lyapunov Analysis
By means of theoretical considerations based on the classical Lyapunov theory and on the property that the kinematic bifurcation rate is much larger than the maximal Lyapunov exponent of the velocity gradient, an interpretation of the kinetic energy cascade phenomenon is given which shows that η ≡ dx is much more rapid and statistically independent with respect to u. Following such considerations, the vorticity equation of a material element (material vorticity)-directly obtained making the curl of the incompressible Navier-Stokes equations-is compared with the evolution equation of η which follows the classical Lyapunov theory. These equations read as From such relations, it is apparent that, for inviscid fluids (ν = 0), the time variations of η and of ω along a fluid particle trajectory x=x(t) follow the same equation, thus ω identifies those particular Lyapunov vectors such that η ∝ ∇ x × u at the initial time t 0 . On the other hand, regardless of the initial condition η(0), η(t) tends to align with the direction of the maximum rising rate of the trajectories distance [56]. If ω(t 0 ) = k η(t 0 ), then ω(t) = k η(t), ∀t > t 0 (von Helmholtz), where k does not depend on t, while η is a fast growth function of t. Thus, following the Lyapunov theory, for inviscid fluids, |ω|, calculated along x=x(t), tends to exponentially rise with t. More in general, for inviscid fluids, ω and η are both fast growth (exponential) functions of the time, where ω tends to align to the direction of maximum growth rate of |η| [56]. A nonzero viscosity influences the time variations of the material vorticity making this latter a slow growth function of t ∈ (t 0 , ∞), whereas η and ξ remain in any case exponential growth functions of t. This implies that, for ν = 0, the characteristic time scales of u (and ϑ) and η are different and that after the time t 0 + a, the fluctuations of ξ result in being statistically independent from u. This holds also when ν → 0 for properly small length scales, except for ν = 0.
Based on the previous observations, the combined effect of very frequent bifurcations and stretching term ∇ x u ω produces the kinetic energy cascade. This phenomenon regards each fluid particle, where ∇ x u ω acts on the material vorticity in the same way in which ∇ x u η influences η. In fact, according to Equation (39), as long as |∇ x u ω| >> ν|∇ 2 x ω|, arbitrary material lines η-thus arbitrary material volumes built on different Lyapunov vectors η, that is, η 1 × η 2 · η 3 -moving along x(t), experience the material vorticity growth and deform according to the Lyapunov theory. According to the analysis of the previous section, such growth phenomenon, due to ∇ x u ω, preserves the average kinetic energy and corresponds to the continuous kinetic energy transfer from large to small scales that is, the kinetic energy cascade phenomenon. Due to the arbitrary choice of x(t), this pertains to all the fluid particles. For what concerns the thermal energy cascade, ϑ is a passive scalar, the temperature follows the velocity fluctuations according to Equation (15), thus the cascade of thermal energy is direct consequence of the mechanism of kinetic energy cascade.
In brief, the energy cascade can be linked to the material vorticity tendency to be proportional to the classical Lyapunov vectors whose modulus changes according to the Lyapunov theory. Specifically, according to Equation (39) and taking into account that S b >> sup λ , η is much faster and statistically independent with respect to the velocity field, while the energy cascade can be viewed as a continuous and intensive stretching and folding process of fluid particles which involves smaller and smaller length scales during their motion and where the particle folding process happens with a frequency given by the bifurcation rate.

Distribution Functions of u, ϑ, x, ξ andλ
Following the present formulation, u, ϑ, x and ξ are the fluid state variables. Therefore, the distribution function of u, ϑ, x and ξ, say P, varies according to the Liouville theorem associated with (4), (5) and (19) [57] ∂P ∂t where, according to the notation of Equations (4) and (5), δ/δu and δ/δϑ are functional partial derivatives with respect to u and ϑ, respectively and (∂/∂•) · stands for the divergence with respect to •. In line with the previous analysis and with References [5,6], P can be factorized as follows being F and P ξ the distribution functions of (u, ϑ) and of (x, ξ), respectively. It is worth remarking that Equation (41) represents the crucial point of this analysis, being the hypothesis of fully developed turbulence following the present formulation. The evolution equations of F and P ξ are formally obtained from Equation (40) and taking into account the aforementioned statistical independence (41). This allows to split the Liouville Equation (40) in the two following equations ∂F ∂t where the boundary conditions of P ξ read as In case of homogeneous and isotropic turbulence, P ξ does not depend on x and can be expressed in function of the finite scale r as follows where δ denotes the Dirac's delta and r k are uniformly distributed points on a sphere S of radius r due to isotropy hypothesis, being k a generic index indicating the several points on S. This leads to Alsoλ and ω E are statistically independent of the velocity field and are continuously distributed in their ranges of variation. In particular, the PDF ofλ, say P λ , can be calculated by means of P ξ with the Frobenius-Perron equation Now, in isotropic turbulence, the longitudinal component of the velocity differenceξ · ξ/r is uniformely distributed in its variation range as ξ sweeps S, while, due to the fluid incompressibility, λ is expected to vary in the interval (−λ S /2, λ S ), where λ S = sup λ . Therefore, substituting Equation (44) in Equation (46), we found thatλ uniformely sweeps (−λ S /2, λ S ), according to Observe that Equations (45) and (47) agree with the results of Reference [6], where the author shows that ξ andλ are both uniformely distributed in their ranges by means of the condition H = max compatible with certain constraints, being H the entropy associated with the kinematic state (x, ξ). This is because the isotropic homogeneous turbulence hypotheses, here expressed through Equations (44) and (45), correspond to the maximum of H. The causes of the nonsymmetric distribution ofλ with respect to the origin, also analyzed in Reference [6], are fluid incompressibility and alignment property of ξ with respect to the maximum rising rate direction. Following such property, regardless of the initial condition ξ(0), ξ(t) tends to align with the direction of the maximum rising rate of the trajectories distance [56]. Therefore, such a distribution function provides positive average Lyapunov exponents and gives the link between average and square mean values of the finite scale Lyapunov exponent according to where • ξ indicates the average of • calculated, through P ξ or P λ .

*Finite Time Lyapunov Exponents and Their Distribution in Fully Developed Turbulence
Altough the local Lyapunov exponentλ quantifies the local trajectories divergence in a point of space, in practice, the trajectory stability is evaluated by observing the particle motion in a finite time interval, say (t 0 , t 0 + τ). For this reason, it is useful to define the finite time Lyapunov exponent as the average ofλ in such time interval, that is This exponent trivially satisfies lim τ→0λ τ =λ.
If τ is properly high, a statistically significant number of kinematic bifurcations n can occur for t ∈ (t 0 , t 0 + τ), thusλ τ is in general a fluctuating variable which exhibits variations whose amplitude diminishes as τ increases. Accordingly,λ τ will be distributed following a Gaussian PDF in fully developed turbulence. In fact, due to the bifurcations encountered in (t 0 , t 0 + τ),λ τ can be written as sum of several terms, each of them related to the effects of a single bifurcation, that is, where ln( k / k−1 ) gives the contribution of the kth bifurcation starting from t 0 , being k−1 and k the Lyapunov vectors moduli calculated immediately before and after the kth bifurcation. On the other hand, due to fully developed chaos, each of such terms is expected to be statistically independent of all other ones and if τ → ∞, the number of encountered bifurcations n diverges. Hence, a proper variant of the central limit theorem can be applied and this would guarantee thatλ τ tends to a Gaussian stochastic variable [58]. The novelty of the present section consists in the implication that the property S b >> λ τ has on Equation (51). Such property should ensure that λ τ can be approximated to a gaussian stochastic variable also for certain finite values of τ. In fact, if τ ≈ 1/λ τ or τ 1/λ τ , the time interval (t 0 , t 0 + τ) should include a statistically significant number of kinematic bifurcations, thus the distribution function of λ τ is expected to be a Gaussian PDF, expecially for relatively high values of the Taylor scale Reynolds number.

Closure of von Kármán-Howarth and Corrsin Equations
Starting from the property of statistical independence (41) and adopting the Liouville theorem, the closure formulas of von Kármán-Howarth and Corrsin equations are here determined and the effects of the chaotic trajectories divergence on these closures are discussed.
In fully developed isotropic homogeneous turbulence, the pair correlation functions of longitudinal velocity components and of temperature, defined as satisfy the von Kármán-Howarth equation [17] and Corrsin equation [19,20], respectively, where von Kármán-Howarth and Corrsin equations are properly obtained from the Navier-Stokes and heat equations written in two points of space, say x and x + r. These correlation equations read as follows The boundary conditions associated with such equations are being u ≡ u 2 r , θ ≡ ϑ 2 , where λ T ≡ −1/ f (0) and λ θ ≡ −2/ f θ (0) are Taylor and Corrsin microscales, respectively. The quantities K and G, arising from inertia forces and convective terms, give the energy cascade and are expressed as [17,19,20] where the repeated index denotes the summation convention. Following the theory [17,19,20], K and G are linked to the longitudinal triple velocity correlation function k and to the triple correlation between u r and ϑ, according to As well known from the literature [17,19,20], without particular hypotheses about the statistics of u and ϑ, K and G are unknown quantities which can not be expressed in terms of f and f θ , thus at this stage of this analysis, both the correlations Equation (54) are not closed.
In order to obtain analytical forms of K and G, observe that these latter, representing the energy flow between length scales in the fluid, do not modify the total amount of kinetic and thermal energies [18,19]. Indeed, convective term, inertia and pressure forces determine interactions between Fourier components of velocity and temperature fields providing the transfer of kinetic and thermal energy between volume elements in the wavenumber space, whereas the global effect of such these interactions leaves u 2 and θ 2 unaltered [18,19]. On the other hand, the proposed statistical independence property (41) allows to write the time derivative of P as sum of two terms the first one of which, being related to ∂F/∂t, provides the time variations of velocity and temperature fields. The second one, linked to ∂P ξ /∂t, not producing a change of u 2 and θ 2 , identifies the energy cascade effect. Therefore, K and G arise from the second term of (58) and can be expressed, by means of the Liouville theorem (40) and Equation (42), in terms of material displacements ξ, taking into account flow homogeneity and fluid incompressibility. Specifically, from Equations (40)-(42), K and G, directly arising from −F∂(P ξξ )/∂ξ, are calculated as follows where U = {u} × {ϑ}, Ξ = {ξ} and dU and dΞ are the corresponding elemental volumes, and Integrating Equation (59) with respect to U , we obtain Again, integrating by parts Equation (61) with respect to Ξ, taking into account the boundary conditions (43) (P ξ ≡ 0, ∀ξ ∈ ∂Ξ) and the isotropy hypothesis, K and G are written as Now, the Lyapunov theory providesξ =λξ + ω E × ξ, and in isotropic homogeneous turbulence P ξ = δ(|ξ| − r)/4πr 2 , thus K and G are Furthermore, the finite scale Lyapunov theory also gives the relationship between velocity correlation and Lyapunov exponents according to where λ ξ and λ 2 ξ are linked with each other through Equation (48), therefore the closure formulas of K and G are in terms of autocorrelations and of their gradients These closure formulas do not include second order derivatives of autocorrelations, thus Equation (65) do not correspond to a diffusive model. The energy cascade expressed by Equation (65) is not based on the eddy viscosity concept, being the result of the trajectories divergence in the continuum fluid. This cascade phenomenon and Equation (65) are here interpreted as follows: (1) In fully developed chaos, the Navier-Stokes bifurcations determine a continuous distribution of velocity, temperature and of length scales, where one single bifurcation causes doubling of velocity, temperature, length scale and of all the properties associated with the velocity and temperature fields according to Equations (15) and (16). This leads to nonsmooth spatial variations of velocity field and very frequent kinematic bifurcations. The main asset of Equation (65) with respect to the other models is that Equation (65) are not based on phenomenological assumptions, such as, for instance, the eddy viscosity paradigm [22][23][24]28,29,33] but are obtained through theoretical considerations concerning the statistical independence of ξ from u and the Liouville theorem.

Remark 1.
At this stage of the present analysis, it is worth remarking on the importance of the hypothesis of the statistical independence of u and ξ expressed by Equation (41). This latter, expressing the hypothesis of fully developed turbulence following this study, leads to the analytical expressions of K and G separating the effects of the trajectories divergence in the physical space from those of the velocity field fluctuations in the Navier-Stokes phase space. Without such hypothesis, the energy cascade effect can not be expressed through the term −F∂(P ξξ )/∂ξ and using Equation (59), thus the proposed closures (65) cannot be determined.
Thanks to their theoretical foundation, Equation (65) do not exhibit free model parameters or empirical constants which have to be identified. These closure formulas coincide with those just obtained by the author in the previous works [1,4,5]. While References [1,4] derive such closures expressing the local fluid act of motion in the finite scale Lyapunov basis and using the frame invariance property of K and G, Reference [5] achieves the same formulas adopting maximum and average finite scale Lyapunov exponents, properly defined and the statistical independence of ξ and u.
Here, unlike References [1,4,5], Equation (65) are determined exploiting the unsymmetric distribution function ofλ just studied in Reference [6], showing that the assumptions of References [1,4,5] are congruent with the present analysis, corroborating the results of the previous work.
References [1,4] show that these closures adequately describe the energy cascade phenomenon and the energy spectra. In detail, K reproduces the kinetic energy cascade mechanism following the Kolmogorov law and G gives the thermal energy cascade in line with the theoretical argumentation of Kolmogorov, Obukhov-Corrsin and Batchelor [35][36][37], with experimental results [38,39] and with numerical data [40,41]. Moreover, Equation (65) allows the calculation of the skewness of ∆u r and ∂u r /∂r which is directly linked to the energy cascade intensity. This is [18] Then, substituting Equation (65) in Equation (66), the skewness of ∂u r /∂r is  Table 1 recalls the comparison, presented in Reference [5,6], between the value of the skewness H 3 (0) of this analysis and those achieved by the aforementioned works. The results were that the maximum absolute difference between the proposed value and the other results were less than 10%. Therefore, the proposed hypotheses, leading to the distribution function (47) and to the closures (65), seem to be adequate assumptions for estimating turbulent energy cascade and spectra.

Properties of the Proposed Closures
Here, some of the properties of the proposed closures (65) are renewed, with particular reference to the evolution times of the developed velocity and temperature autocorrelations. In detail, we will show that these correlations reach their developed shape in finite times which depend on the initial condition and that, after this period, the hypothesis of statistical independence could be not more verified. This result is given in Reference [5], where the author adopts a specific Lyapunov analysis using two exponents properly defined. Unlike in Reference [5], such a result is here achieved through the previously obtained local finite scale Lyapunov exponent distribution (47). To analyze this, the evolution equations of u, θ, λ T and λ θ are first obtained taking the coefficients of order r 0 and r 2 of Equation (54) arising from the Taylor series expansion of even powers of f and f θ [17,19,20] This leads to the following equations While Equation (69) do not depend on the particular adopted closures [17,19,20], Equation (70) are obtained using the proposed closures (65). On the other hand, it is useful to consider the fluctuations of the classical Lyapunov exponent, defined as which are related to f through Equations (64) and (68) in such a way that being Λ the root mean square ofΛ. Following Equation (70), the time variations of λ T , λ θ and Λ are now discussed. The first terms at the R.H.S. of Equation (70) provide the turbulent energy cascade, whereas the other ones arise from the fluid diffusivities. While these latter contribute to increasing both the correlation lengths, the energy cascade mechanism tends to reduce these scales and if such a mechanism is sufficiently stronger than diffusivities, then dλ T /dt < 0 and dλ θ /dt < 0.
For sake of our convenience, the condition ν = 0, χ = 0 is first studied. In this case, u and θ are both constants, whereas λ T , λ θ and Λ vary with t. In detail, λ T and λ θ are proportional to each other and vary linearly with time according to while Λ monotonically rises and goes to infinity in a finite time, being τ the dimensionless time.
When ν = χ = 0, the energy cascade provides that both the microscales decrease until to τ → 2, where both the correlations are considered to be fully developed, λ T → 0, λ θ → 0 and Λ → ∞ (see solid lines of Figure 4). Thus, the two correlations will exhibit developed shapes in finite times whose values depend on the initial condition Λ(0). The meaning that both the microscales are decreasing functions of τ is that kinetic and thermal energies are continuously transferred from large to small scales following the previous scheme. Next, as τ → 2, Λ → +∞ and this means that the velocity gradient diverges in a finite time depending on Λ(0) and that contiguous particles trajectories diverge with a growth rate infinitely faster than velocity and temperature fields. For ν > 0, χ > 0, then du/dt < 0 and dθ/dt < 0 in any case and f and f θ are here supposed to be fully developed as soon as dλ T /dt = 0 and dλ θ /dt = 0, respectively. These situations are qualitatively shown in the figure by the dashed lines for different values of R T and Pe, where R T = λ T u/ν and Pe = Pr R T are, respectively, Reynolds number and Péclet number, both referred to the Taylor microscale, being Pr = ν/χ the Prandtl number. When the initial microscales are relatively large, the diffusivities effects are quite smaller than the convective terms, the energy cascade is initially stronger than the diffusivities effects and both the microscales exhibit about the same trend just discussed for ν = χ = 0. According to Equations (69) and (70), the interval where τ ranges can be splitted in two subregions for both f and f θ . The first ones correspond to values of τ ∈ (0, 2) such that dλ T /dt < 0 and dλ θ /dt < 0, which are upper bounded by the endpoints τ 1 < 2, τ 2 < 2 where dλ T /dt(τ 1 ) = 0 and dλ θ /dt(τ 2 ) = 0 (dashed lines), respectively, being in general τ 1 = τ 2 . There, the kinetic and thermal energy cascades are momentarily balanced by viscosity and thermal diffusivity, respectively and both the autocorrelations can be considered fully developed. For both the correlations, such momentary balance happens in finite times τ < 2 which depend on the initial condition. As far as Λ is concerned, this initially coincides about with that obtained for ν = 0, then reaches its maximum for τ 2 and thereafter diminishes due to viscosity. When Λ achieves its maximum, dΛ/dt = 0, chaos and mixing reach their maximum levels, the correlations are about fully developed, thus relative kinematics and fluid strain change much more rapidly than velocity field. Thereafter, we observe regions where dΛ/dt < 0. There, due to the relatively smaller values of the microscales, the dissipation is stronger than the energy cascade and both the correlation lengths tend to rise according to Equation (70). Such a region, which occurs immediately after the condition dΛ/dt = 0, corresponds to the regime of decaying turbulence.
Observe that the proposed closures (65) are expected to be verified where dΛ/dt > 0, in which the Navier-Stokes bifurcations generate the regime of fully developed turbulence. On the contrary, in regime of decaying turbulence -dΛ/dt < 0-, after a certain time, say τ + > τ 1 ≈ 2, it results Λ/Λ(0) < 1. In such situations, the relative kinematic and fluid strain could be not faster than velocity field, thus the statistical independence hypothesis (41) could be not satisfied and Equation (65) will be not defined. Therefore, the condition τ ≈ 2 or Λ/Λ(0) < 1 provides a further limit of validity for the proposed closure formulas.

*Self-Similarity and Developed Correlations of the Proposed Closures
This section analyzes self-similarity and developed shape of f and f θ produced by the proposed closures. The new result with respect to the previous works consists in to remark that the proposed closures generate correlations self-similarity in proper ranges of r, which is directly related to the fluid trajectories divergence. To study this question, observe that a given function of t and r, say ψ = ψ(t, r), which completely exhibits self-similarity with respect to r as t changes, is a function of the kind ψ(t, r) = ψ r L(t) (74) and exactly satisfies the equation whereinL is the characteristic length associated with the specific problem. From such equation, the self-similarity of ψ is linked to the variation rate d lnL(t)/dt. Now, thanks to the mathematical structures of the proposed closures (65), and taking into account that f and f θ are both even functions of r which near the origin behave like Equation (68), K and G can be expressed through even power series of f as follows thus, the evolution equations of both the autocorrelations can be written in the following way Comparing Equations (75) and (77), it follows that the proposed closures (65) generate self-similarity in a range of variation of r where Λ/2r∂ f /∂r and Λ/2r∂ f θ /∂r are dominant with respect to the other terms. As the result, such self-similarity is directly caused by the continuous fluid trajectory divergence-quantified by Λ-which happens thank to very frequent kinematic bifurcations. In such these intervals, the correlations will exhibit self-similarity during their time evolution, thus f and f θ can be expressed there as follows In such regions, the energy cascade is intensive and much stronger than the diffusivities effects, thus following Equation (70), λ θ (t) is expected to be proportional to λ T (t) Next, as ϑ is a passive scalar, energy cascade and fluid diffusivities act on u and θ in such a way that their increments are proportional with each other. Therefore, far from the initial condition, we expect that Now, Equation (79) provides a link between the correlation scales and Pr. In fact, substituting Equation (79) in Equation (69), we obtain Furthermore, from Equation (70), also f IV (0) and f IV θ (0) are related to the Prandtl number Hence, the developed autocorrelations can be estimated searching for the solutions of the closed von Kármán-Howarth and Corrsin equations in the self-similar form (78) when dλ T /dt = dλ θ /dt = 0. This leads to the following ordinary differential equations system Several solutions of these equations were numerically obtained in [2,4], where the author shows that velocity and temperature correlations agree with the Kolmogorov law, with the theoretical arguments of Obukhov-Corrsin and Batchelor and with the numerical simulations and experiments known from the literature [19,[35][36][37][38][39][40][41].
For sake of reader convenience, Figures 5 and 6 report the velocity correlations and the corresponding spectra E(κ), T(κ) numerically calculated with the first equation of Equation (83) for R T = 100, 200, 300, 400, 500, 600, being where all these cases correspond to the same level of average kinetic energy. The integral correlation scale of f results to be a rising function of R T , while the triple longitudinal velocity correlation k maintains negative with a minimum of about −0.04 whose value is achieved for values of r/λ T which rise with the Reynolds number. For what concerns the spectra, observe that increasing κ, the kinetic energy spectra behave like E(κ) ≈ κ 4 near the origin, then exhibit a maximum and thereafter are about parallel to the dashed line κ −5/3 in a given interval of the wave-numbers. The size of this latter, which defines the inertial range of Kolmogorov, rises as R T increases. For higher values of κ, which correspond to scales less than the Kolmogorov length, E(κ) decreases more rapidly than in the inertial range. As K does not modify the kinetic energy, the proposed closure gives ∞ 0 T(κ)dκ ≡ 0.  From these solutions, the Kolmogorov constant C, here calculated as is shown in Table 2 in function of the Reynolds number, where ε = −3/2 du 2 /dt. The obtained values of C ≈ 2 are in good agreement with the corresponding values known from the literature. Next, Figure 7 shows the temperature spectra Θ(κ) and the temperature transfer function Γ(κ) calculated as follows [65]    in such a way that The variations of Θ(κ) with R T and Pr are quite peculiar and consistent with previous studies according to which there are regions where Θ(κ) exhibits different scaling laws Θ(κ) ≈ κ n . Following the proposed closures, n 2 when κ → 0 in any case. For Pr = 0.001, when R T ranges from 50 to 300, the temperature spectrum essentially exhibits two regions: one in proximity of the origin where n 2 and the other one, at higher values of κ, where −17/3 < n < −11/3, (value very close to −13/3). The value of n ≈ −13/3, here obtained in an interval around tor ≈1, is in between the exponent proposed by [36] (−17/3) and the value determined by [40] (−11/3) by means of numerical simulations. Increasing κ, n significantly diminishes and Θ(κ) does not show scaling law. When Pr = 0.01, an interval nearr ≈ 1 where −17/3 < n < −13/3 appears and this is in agreement with [36]. Next, for Pr = 0.1, the previous scaling law vanishes, whereas for R T = 50 and 100, n changes with κ and Θ(κ) does not show clear scaling laws. When R = 300, the birth of a small region is observed, where n ≈ −5/3 has an inflection point. For Pr = 0.7 and 1, with R T = 300, the width of this region is increased, whereas at Pr = 10 and R = 300, we observe two regions: one interval where n has a local minimum with n −5/3 and the other one where n exhibits a relative maximum, with n −1. For larger κ, n diminishes and the scaling laws disappear. The presence of the scaling law n −5/3 agrees with the theoretical arguments of [20,37] (see also [39,41] and references therein). Figure 7 also reports (on the bottom) the spectra Γ(κ) (solid lines) and T(κ) (dashed lines) which describe the energy cascade mechanism.

*Bifurcation Analysis of Closed von Kármán-Howarth Equation: From Fully Developed Turbulence Toward Non-Chaotic Regimes
Starting from non-chaotic regimes, the transition towards the fully developed turbulence happens through intermediate stages [34,42,43,51] which correspond to bifurcations where the relative Reynolds numbers show the same order of magnitude. This section presents a specific bifurcation analysis, which, unlike the classical route toward the chaos [34,42,43,51], analyzes the inverse route: the starting condition is represented by the fully developed homogeneous isotropic turbulence and the route followed is that towards the non-chaotic regime. Such route corresponds to the path B → A of Figures 1f,g and 2. Along the line B → A, R T gradually diminishes and the bifurcations of the closed von Kármán-Howarth equation, properly defined, will be here studied. This analysis estimates R * T through the closures (65) and their previously seen properties, where R * T defines the minimum value of R T for which the turbulence maintains fully developed, homogeneous and isotropic. This provides the order of maginitude of Re at the transition, indicating a further limit of the proposed closures.
In order to formulate a bifurcation analysis for the velocity correlation equation, consider now the various coefficients of the closed von Kármán-Howarth equation which arise from the even Taylor series expansion of f (t, r) = ∑ k f (k) 0 r k /k!. Each of such these coefficients corresponds to one of the following equations Such equations can be written by introducing the infinite dimensional state vector which represents the state of the longitudinal velocity correlation. Therefore, Equation (88), formally written asẎ are equivalent to the closed von Kármán-Howarth equation. Equation (90) defines a bifurcation problem where ν plays the role of control parameter. Thus, this bifurcation analysis studies the variations of Y caused by ν according to For ν > ν 0 , Y is formally calculated through the implicit functions inversion theorem where ∇ Y F is the jacobian ∂F/∂y. A bifurcation of Equation (90) happens when this jacobian is singular, that is, If ν 0 is quite small (R T properly large), the energy cascade is dominant with respect to the viscosity effects and ∇ Y F is expected to be nonsingular. Increasing ν, Y smoothly varies according to Equation (92) and thereafter the dissipation gradually becomes stronger than the energy cascade until reaching the first bifurcation where condition (93) occurs. With reference to Figure 2, this corresponds to the path B → A until to reach A. There, a hard loss of stability is expected for the fully developed turbulence toward non-chaotic regimes [66]. Therefore, R * T is calculated as that value of R T at bifurcation which gives the maximum of the largest real part of the eigenvalues of ∇ Y F [66,67] compatible with the current value of the average kinetic energy u 2 , that is, where l k , k=1, 2,... are the eigenvalues of ∇ Y F.
On the other hand, as previously seen, far from the initial condition, the energy cascade acts keeping f similar in the time in a given interval of variation of r. There, the evolution of f is expected to be described-at least in first approximation-by Equation (78) and this suggests that-under such approximation-the knowledge of u and λ T can be considered to be sufficient to describe the evolution of f . Hence, only the first two components of the state vector Y are taken which correspond to the coefficients of the order of r 0 and r 2 of Equation (88). Thus, thanks to the self-similarity, the infinite dimensional space where Y lies is replaced by a finite dimensional manifold and the state vector is reduced to f IV 0 plays the role of a parameter which characterizes the velocity correlation and the jacobian ∇ Y F reads as From Equation (97), as long as ν > 0 is properly small, det (∇ Y F) > 0. In order that a bifurcation happen, det (∇ Y F) must vanish for a certain value of ν and this implies that f IV 0 λ 4 T > −10/7. Thus, increasing ν, det (∇ Y F) /ν diminishes and there exists a value of ν where this jacobian determinant vanishes. To determine R * T , f IV 0 is eliminated through the bifurcation condition (det (∇ Y F) = 0) and Equation (97), that is, Therefore, the singular jacobian is and admits the following eigenvalues and eigenvectors l 1 , l 2 , y 1 and y 2 , respectively The eigenvalue l 2 ∈ R maintains positive for R T > 5 and reaches its maximum l 2max = 5ν/λ 2 T for R T = 10. Accordingly, R * T is estimated as which corresponds to f IV 0 = 0. Another characteristic value of R T is obtained in the case where both the eigenvalues vanish. This is R T = 5 and is expected to represent the onset of the decaying turbulence regime. In fact, in such situation, it is reasonable that f and λ T are Hence, f IV 0 λ 4 T 3 and R T 4, in agreement with the previous estimation.

Remark 2.
It is worth remarking that R * T provides the minimum of R T in fully developed isotropic homogeneous turbulence, thus this gives the order of magnitude of R T at the transition. Of course, the transition toward the chaos consists in intermediate stages (bifurcations of Navier-Stokes equations) where the turbulence is not developed and the velocity statistics does not exhibit, in general, isotropy and homogeneity. Hence, the obtained results provide the order of magnitude of R T at the transition. On the basis of this analysis, during the transition, R T ranges as The obtained value of R * T = 10 is in very good agreement with the bifurcations analysis of the turbulent energy cascade [3], where the author shows that, in the transition toward the developed turbulence, if the bifurcations cascade follows the Feigenbaum scenario [42,43], the critical Taylor scale Reynolds number is about 10.13 and occurs after three bifurcations.
We conclude this section by remarking the limits under which R * T is estimated. Such limits derive from the local self-similarity produced by the closures (65) which allow to consider only the first two equations of (88).

Velocity and Temperature Fluctuations
The purpose of this section is to obtain, by means of the previous Lyapunov analysis, formal expressions of velocity and temperature fluctuations which will be useful for estimating the statistics of these latter. For sake of our convenience, Navier-Stokes and thermal energy equations are now written in the following dimensionless divergence form where T and q denote, respectively, dimensionless stress tensor and heat flux, according to the Navier-Fourier laws being I the identity tensor, T v the viscous stress tensor and the pressure p is given according to Equation (6).
In order to obtain the analytical forms of velocity and temperature fluctuations, Equation (104) are first expressed in terms of referential coordinate x 0 where the repeated index denotes the summation convention. The adoption of the referential coordinates allows to factorize of ∂u/∂t and ∂ϑ/∂t as a product of two statistically uncorrelated matrices: one depending on velocity and temperature fields and the other representing the local fluid deformation. Velocity and temperature fluctuations are here obtained integrating Equation (106) in the set (t, a). Due to the alignment property of the Lyapunov vectors [56], exp(−Λt) rapidly goes to zero as t → ∞ in any case, whereas ∂T ij /∂x 0k and ∂q j /∂x 0k are functions of slow growth of t. Hence, velocity and temperature fluctuations are formally calculated integrating Equation (106) in the set (t, ∞) where ∂T ij /∂x 0k and ∂q j /∂x 0k are considered to be constant and equal to the corresponding values at the current time. Such fluctuations are then expressed in function of current velocity and temperature fields according to where |W jk | < ∞ as G −1 jk is represented by slow growth functions of t. It is worth to remark that Equation (107) are, in general, rough approximations of velocity and temperature fluctuations. Nevertheless, in fully developed turbulence, dx(t) is considered to be much more rapid than u(t, x), thus Equation (107) provide one accurate way to express velocity and temperature in terms of referential coordinates by means of the Lyapunov theory.

*Statistics of Velocity and Temperature Difference
In developed turbulence, longitudinal velocity and temperature difference, ∆u r = (u(t, x ) − u(t, x)) · r/r and ∆ϑ = ϑ(t, x ) − ϑ(t, x), r = x − x, play a role of paramount importance as these quantities describe energy cascade, intermittency and are linked to dissipation. This section analyzes the statistics of such quantities in fully developed homogeneous isotropic turbulence through the previously seen kinematic Lyapunov analysis and using a proper statistical decomposition of velocity and temperature. In order to determine this statistics, the Navier-Stokes bifurcations effect on ∆u r and ∆ϑ is first analyzed. To this purpose, ∆u r and ∆ϑ are expressed in function of current velocity and temperature through Equation (107) The several bifurcations happening during the fluid motion determine a continuous doubling of u in several functions, sayv k , k = 1, 2, ..., in the sense that each encountered bifurcation introduces new functionsv k whose characteristics are independent of the velocity field at previous time. Then, due to bifurcations, u is of the form It is worth remarking that, while u(t, x) is solution of the Navier-Stokes equations, the functionŝ v k are not. Therefore, the functionsv k are the result of the mathematical segregation due to bifurcations of a fluid state variable which physically only exist in combination, thus each of them is not directly observable. This implies that u will be distributed, in line with the Liouville theorem, according to a classical definite positive distribution function. On the contrary, each single functionv k , representing mathematical segregation of the fluid state, will be distributed following extended distribution functions which can exhibit negative values [68][69][70] compatible with conditions linked to the specific problem. These conditions mainly arise from (a) the Navier-Stokes equations and from (b) the isotropic hypothesis. For what concerns (a), in order that pressure and inertia forces can cause sizable variations of velocity autocorrelation, each termv k ≡ (v 1 ,v 2 ,v 3 ) will be distributed following highly nonsymmetric extended distribution function, for which As for (b), due to isotropic hypothesis, u would be distributed following a gaussian PDF [18], thus, according to the Navier-Stokes equations, pressure and inertia forces will not give contribution to the time derivative of the third statistical moment of u. Accordingly, the absolute value of odd statistical moments of order n ofv k is expected to be very high in comparison with the even statistical moments of order n + 1, that is, ki (n+1)/2 , n = 3, 5, 7, ..., i = 1, 2, 3.
This suggests that ∆u and u can be expressed through a specific statistical decomposition [71], as a linear combination of opportune stochastic variables ξ k that reproduce the doubling bifurcations effect and whose extended distribution functions satisfy Equations (111) and (112). Furthermore, as ϑ is a passive scalar, its fluctuations are the result of u and of thermal diffusivity, thus also ϑ is written by means of the same decomposition where U k and Θ k (k = 1, 2, ...) are coordinate functions of t and x, being ∇ x · U k = 0, ∀k and ξ k (k = 1, 2, ...) are dimensionless independent centered stochastic variables such that where q, providing the skewness of ξ k k = 1, 2..., satisfies to Therefore, the distribution functions of ξ k can assume negative values compatible with Equations (114) and (115).
Through the decomposition (113), we will show that the negative value of H where k ξ k are the contributions of inertia and pressure forces and of the fluid viscosity, respectively, whereas ∑ j ∑ k B jk ξ j ξ k and 1/Pe ∑ k b k ξ k arise from the convective term and fluid conduction. Because of turbulent isotropy, it is reasonable that u i and ϑ are both Gaussian stochastic variables [18,71,72], thus the various terms of Equation (116) satisfy the Lindeberg condition, a very general, necessary and sufficient condition for satisfying the central limit theorem [71,72]. Such theorem does not apply to ∆u i and ∆ϑ as these latter are the difference between two correlated Gaussian variables, thus their PDF are expected to be very different with respect to Gaussian distributions. To study the statistics of ∆u r and ∆ϑ, the fluctuations of these latter are first expressed in terms of ξ k ∆u r (r) = ∑ j ∑ k ∆A jk ξ j ξ k + 1 being In Equation (118), the matrices ∆A jk and ∆B jk are decomposed following their symmetric and antisymmetric parts, respectively S ujk , S θjk and Ω ujk , Ω θjk . These last ones give null contribution in Equation (117), whereas the terms arising from S ujk and S θjk are expressed as in which the first term of Equation (119) is decomposed in the following manner being and where ∑ + and ∑ − denote summations for (S Xjj > 0, S Xkk > 0) and (S Xjj ≤ 0, S Xkk ≤ 0) and n + X and n − X are the corresponding numbers of terms of such summations, whereas ∑ + j =k and ∑ − j =k indicate the sums of addends calculated for j = k corresponding to S Xjj > 0, S Xkk > 0 and S Xjj < 0, S Xkk < 0, respectively. The decomposition (119) and (120) and the definitions (121) lead to the following expression of velocity and temperature difference fluctuations Now, we show that ξ X , η X and ζ X , X = u, θ tend to uncorrelated gaussian variables. In fact, from Equation (121), η X and ζ X , X = u, θ are sums of random terms belonging to two different sets of uncorrelated stochastic variables (i.e., the sets for which S Xii < 0 and S Xii > 0), therefore η X and ζ X , are two uncorrelated stochastic variables such that η X = ζ X = 0, X = u, θ. Furthermore, as ξ k are statistically independent with each other, the central limit theorem applied to Equation (121) guarantees that both η X and ζ X tend to two uncorrelated centered gaussian random variables. As for ξ X , X = u, θ, the following should be considered: due to the analytical structure of Equation (122), each term of ξ X is a centered variable, thus ξ X = 0. Next, in Equation (122), the following terms i are mutually uncorrelated, as each of these is sum of random variables belonging to two different uncorrelated sets. Moreover, ∑ i =j ξ i ξ j includes several weakly correlated terms, whereas ∑ k g Xk ξ k is the sum of independent variables. On the other hand, due to hypothesis of fully developed chaos, the energy cascade, here represented by Equations (114), (115) and (117), will generate a strong mixing on the several terms of Equation (117), thus a proper variant of the central limit theorem can be applied to ξ X whose several terms are weakly dependent with each other [72]. As the result, ξ X , X = u, θ will tend to centered gaussian variables statistically independent of η X and ζ X .
Hence, the statistics of ∆u r and ∆ϑ is represented by the following structure functions of the independent centered gaussian stochastic variables ξ X , η X and ζ X for which ξ 2 X = η 2 X = ζ 2 X = 1.
where L u and L θ are now introduced to take into account that ξ X , η X and ζ X have standard deviation equal to unity. Thus and L X , S − X and S + X are parameters depending upon r which have to be determined. To this regard, it worth remarking that, in regime of fully developed isotropic turbulence in infinite domain, the numbers of parameters necessary to describe the statistics of ∆u r and ∆ϑ should be minimum compatible with assigned quantities which define the current state of fluid motion, such as average kinetic energy, temperature standard deviation and correlation functions. On the other hand, the evolution equation of f [17] requires the knowledge of the correlations of the third order k to be solved. Therefore, in fully developed homogeneous isotropic turbulence, the sole knowledge of f and k is here considered to be the necessary and sufficient information for determining the statistics of ∆u r . This implies that S + u is proportional to S − u through a proper quantity which does not depend on r, that is, where χ < 1 is a function of R T giving the skewness of ∆u r , which has to be identified. Accordingly, S u and L u will be determined in function of f and k as soon as χ = χ(Re) is known. For what concerns the temperature difference, observe that, due to turbulence isotropy, the skewness of ∆ϑ should be equal to zero and this gives Therefore, the structure functions of ∆u r and ∆ϑ read as Furthermore, again following the parameters minimum number, the ratio Ψ θ (r) ≡ S θ /L θ would be proportional to Ψ u (r) ≡ S u /L u through a proper coefficient depending upon the Prandtl number alone, that is where σ is a function of the Prandtl number which has to be determined. At this stage of the present analysis, we show that, in fully developed turbulence, L u and L θ are, respectively, functions of R T and Pe, resulting in L u ∝ R −1/2 T and L θ ∝ Pe −1/2 . In fact, from Equation (125) we obtain As | ξ 3 k | >>> 1, ξ i ξ j ξ k ξ l , first and third addend of Equation (130) are negligible with respect to second one, thus L u and L θ tend to functions of the kind where F u (r) and F θ (r) are functions of r which do not directly depend on R T and Pe. Hence, the dimensionless ∆u r and ∆ϑ, normalized with respect to the corresponding standard deviations, are expressed in function of R T and Pe and this identifies σ = √ Pr. Equation (132) provide peculiar structure functions giving the statistics of ∆u r and ∆ϑ. Now, if χ = χ(R T ) is considered to be known, L u and S u can be expressed in function of ∆u 2 r and ∆u 3 r , where this latter is calculated adopting the proposed closure (65). In fact, L u and S u are related to ∆u 2 r and ∆u 3 r through Equation (128) thus, L u , S u and Φ are expressed in function of f (r) and k(r) as S u (r) = 3/4 In the expression of L u (r) of Equations (134), the argument of the square root must be greater than zero and this leads to the following implicit condition for χ where the proposed closure (65) is taken into account. Inequality (135), solved with respect to χ, gives the upper limit for χ As far as the temperature difference is concerned, we have thus Equation (137) allows to calculate L θ in terms of the other quantities In Equations (134) and (138), the function χ = χ(R T ) has to be identified, and Φ(r) depends on the specific shape of f (r), where, due to the constancy of H The distribution functions of ∆u r and ∆ϑ are formally calculated through the Frobenius-Perron equation [57], taking into account that ξ X , η X and ζ X are independent identically distributed centered gaussian variables such that ξ 2 where δ is the Dirac delta, P(ξ, η, ζ) is the 3D gaussian PDF and ∆u r (ξ, η, ζ)) and ϑ(ξ, η, ζ)) are determined by Equation (132). In other words, the statistics of ∆u r and ∆ϑ can be inferred looking at the proposed statistical decomposition (113) which includes the bifurcations effects in isotropic turbulence. This is a non-Gaussian statistics, where the absolute value of the dimensionless statistical moments increases with R T and Pe. In detail, the dimensionless statistical moments of ∆u r and ∆ϑ are easily calculated in function of χ, Ψ u and Ψ θ where Φ(0) and χ = χ(R T ) have to be identified. To this end, we first analyze the statistics of ∂u r /∂r which, following the proposed Lyapunov analysis, exhibits a constant skewness H and H Accordingly, χ = χ(R T ) is implicitly expressed in function of Φ(0) √ R T . From Equation (143), u (0) = −3/7, admits limit resulting in χ(R T ) < 0 for properly small values of R T . On the other hand, in fully developed turbulence, the PDF of ∂u r /∂r exhibits non gaussian behavior (i.e., non gaussian tails) for ∂u r /∂r → ±∞, accordingly χ must be positive. Hence, the limit condition χ = 0 is supposed to be achieved for R T = R * T = 10 which represents the minimum value of R T for which the turbulence is homogeneous isotropic. This allows to identify Φ(0) by means of Equation (143) Thus, Equation (143) gives, in the implicit form, the variation law χ = χ(R T ) which is depicted in Figure 8. We conclude this section with the following considerations regarding the proposed analysis, and summarizing some of the results just obtained in the previous works.
For non-isotropic turbulence or in more complex situations with boundary conditions or walls, the velocity will be not distributed following a normal PDF, thus Equation (112) will be not verified, and Equation (132) will change its analytical structure incorporating stronger intermittent terms [72] giving the deviation with respect to the isotropic turbulence. Hence, the absolute statistical moments of ∆u r will be greater than those calculated through Equation (141), indicating that, in more complex cases than the isotropic turbulence, the intermittency of ∆u r can be significantly stronger.
Next, Ψ u and Ψ θ represent the ratios (large scale velocity)-(small scale velocity) and (large scale temperature)-(small scale temperature), respectively. In particular, Ψ u ∝ u/u s ≈ (u 2 /λ T )/(u 2 s /l s ) being l s and u s the characteristic small scale and the corresponding velocity. This means that u/u s ≈ λ T /l s ≈ √ R T and that the Reynolds number relative to u s and l s is u s l s /ν ≈ 1, that is l s and u s identify the Kolmogorov scale and the corresponding velocity. For what concerns Ψ θ , ϑ is a passive scalar, thus Ψ θ reads as Ψ θ ∝ θ/θ s ≈ θ/θ s (u/λ T )/(u s /l s )and this leads to u s l s /ν ≈ 1.
At this stage of the present analysis, we can show that the kinematic bifurcation rate S b , defined by Equation (25), is much larger than the kinematic Lyapunov exponents. In fact, S b can be also estimated as the ratio (large scale velocity)-(small scale length), where large scale velocity and small scale length are given by u and by the Kolmogorov scale, respectively. Taking into account the Kolmogorov scale definition and Equation (69), we obtain confirming the assumption made in the relative section. In fully developed turbulence, S b >> Λ and is a rising function of R T . As shown in Reference [1], the statistics given by Equations (139) and (141) agree with the experimental data presented in References [47,48]. There, in experiments using low temperature helium gas between two counter-rotating cylinders (closed cell), the PDF of ∂u r /∂r and its statistical moments are measured. Although the experiments regard wall-bounded flows, the measured PDF of velocity difference are comparable with the present results ( Equations (139) and (141)). Apart from a lightly non-monotonic evolution of H (4) u (0) and H (6) u (0) in [47,48], the dimensionless statistical moments of ∂u r /∂r exhibit same trend and same order of magnitude of the corresponding quantities calculated with Equation (141). In particular, the PDFs of ∂u r /∂r obtained with the present analysis show non gaussian tails which coincide with those measured in [47,48].
In Figure 9, the normalized PDFs of ∂u r /∂r, calculated with Equations (139) and (141) Figure 9c represents the enlarged region of Figure 9b, showing the tails of PDF for 5 < s < 8. According to Equations (139) and (141), the tails of the PDF rise with the Reynolds number in the interval 10 < R T < 700, whereas for R T > 700, smaller variations are observed. On the right-bottom, the results of [47] for R T = 255, 416, 514, 1035 and 1553 are shown. Despite the aforementioned non-monotonic trend (see Figure 9 (Right-bottom)), Figure 9c gives values of the PDFs and of the corresponding average slopes which agree with those obtained in [47], expecially for 5 < s < 8. To this regards, it is worth to remark that, in certain conditions, the flow obtained in the experiments of [47] could be quite far from the isotropy hypothesis, as such experiments pertain wall-bounded flows, where the walls could significantly influence the fluid velocity in proximity of the probe. In References [1,2,4] the scaling exponents ζ V (n) associated with the several moments of ∆u r (∆u r ) n ≈ A n r ζ V (n) , are calculated with Equation (132) through the following best fitting procedure. The statistical moments of ∆u r are first calculated in function of r using Equation (141) (see Figure 10 (Left)). Then, the scaling exponents ζ V (n) are identified through a minimum square method which, for each statistical moment, is applied to the following optimization problem J n (ζ V (n), A n )≡ r 2 r 1 ( (∆u r ) n − A n r ζ V (n) ) 2 dr = min, n = 1, 2, ...
where ( (∆u r ) n ) are calculated with Equation (141),r 1 is assumed to be equal to 0.1, whereasr 2 is taken in such a way that ζ V (3) = 1. The so obtained scaling exponents are shown in Figure 10 (Right-side) (solid symbols) where these are compared with those given by the Kolmogorov theories K41 [44] (dashed line) and K62 [45] (dotted line) and with the exponents calculated by She-Leveque [46] (continuous curve). For n < 4, ζ V (n) ≈ n/3 and for higher values of n, due to the nonlinear terms of Equation (132), ζ V (n) shows multiscaling behavior. The values of ζ V (n) here calculated are in good agreement with the She-Leveque data, and result to be lightly greater than those obtained in [46] for n > 8.   [44]. Dotted line is for Kolmogorov K62 data [45]. Continuous line is for She-Leveque data [46].
As far as the temperature difference statistics is concerned, Figure 11 (Left) shows the distribution function of ∂ϑ/∂r in terms of dimensionless abscissa s = ∂ϑ/∂r statistical moments of ∂ϑ/∂r and of ϕ and to the approximation of assuming the components of ∇ϑ to be statistically uncorrelated. Finally, observe that the experimental data of [47,74] allow to identify Φ(0). Table 3 reports a comparison between the value of Φ(0) calculated with the present theory and those obtained through elaboration of the experimental data of [47,74]. Form this comparison, the value of Φ(0) calculated with Equation (145) is in very good agreement with those obtained through the elaboration the data of [47,74]. Table 3. Identification of Φ(0) through elaboration of experimental data of [47,74] and comparison with the present analysis.

Conclusions
A review of previous theoretical results concerning an original turbulence theory is presented. The theoretical approaches here adopted, different with respect to the other articles, confirm and corroborate the results of the previous works.
In separate sections, novel issues regarding the proposed turbulence theory are presented, and are here summarized. - The bifurcation rate of velocity gradient, calculated along fluid particles trajectories is shown to be much larger than the maximal Lyapunov exponent of the kinematic field. -On the basis of the previous item, the energy cascade is viewed as a stretching and folding succession of fluid particles which gradually involves smaller and smaller scales. - The central limit theorem, in the framework of the bifurcation analysis, provides reasonable argumentation that the finite time Lyapunov exponent can be approximated by a gaussian random variable if τ ≈ 1/Λ. - The closures of von Kármán-Howarth and Corrsin equations given by this theory determine velocity and temperature correlations which exhibit local self-similarity directly linked to the continuous particles trajectories divergence. - The proposed bifurcation analysis of the closed von Kármán-Howarth equation studies the route from developed turbulence toward non-chaotic regimes and leads to an estimation of the critical Taylor scale Reynolds number in isotropic turbulence in agreement with the various experiments. -Finally, a specific statistical decomposition of velocity and temperature is presented. This decomposition, adopting random variables distributed following extended distribution functions, leads to the statistics of velocity and temperature difference which agrees with the data of experiments.
Funding: This research received no external funding