Tsallis Extended Thermodynamics Applied to 2-d Turbulence: Lévy Statistics and q-Fractional Generalized Kraichnanian Energy and Enstrophy Spectra

The extended thermodynamics of Tsallis is reviewed in detail and applied to turbulence. It is based on a generalization of the exponential and logarithmic functions with a parameter q. By applying this nonequilibrium thermodynamics, the Boltzmann-Gibbs thermodynamic approach of Kraichnan to 2-d turbulence is generalized. This physical modeling implies fractional calculus methods, obeying anomalous diffusion, described by Lévy statistics with q < 5/3 (sub diffusion), q = 5/3 (normal or Brownian diffusion) and q > 5/3 (super diffusion). The generalized energy spectrum of Kraichnan, occurring at small wave numbers k, now reveals the more general and precise result k−q. This corresponds well for q = 5/3 with the Kolmogorov-Oboukov energy spectrum and for q > 5/3 to turbulence with intermittency. The enstrophy spectrum, occurring at large wave numbers k, leads to a k−3q power law, suggesting that large wave-number eddies are in thermodynamic equilibrium, which is characterized by q = 1, finally resulting in Kraichnan’s correct k−3 enstrophy spectrum. The theory reveals in a natural manner a generalized temperature of turbulence, which in the non-equilibrium energy transfer domain decreases with wave number and shows an energy equipartition law with a constant generalized temperature in the equilibrium enstrophy transfer domain. The article contains numerous new results; some are stated in form of eight new (proven) propositions.


Introduction
There are primarily two limiting cases to describe turbulent flows, a microscopic and a macroscopic description. The first is based on the entire ensemble of atoms and molecules in the flow field. The characteristic amount of a considered entity, respectively the number of degrees of freedom of such a large-scale system, e.g., an atmospheric turbulent flow, is of the order of the value of the Avogadro number (6.022 × 10 23 ), which is very large. In 1996 Castaign [1] pointed out that controlling each of these single degrees is impossible and a futile program. The related problem is so large that it can at best be handled by molecular dynamics and supercomputing. On the macroscopic side, it is the description of fluid behavior by the Navier-Stokes Equations (NSE), which reduce the problem to only four macroscopic variables (however, be aware that these variables are functions of space and time), namely three velocity components and the pressure. Naively, one could now assume e.g., dilute gases, Boson gases, photon radiation, phonons in solids, fermions and ferromagnetic magnons at low temperatures, etc. (e.g., [61,63], etc.).
The systems, studied in this section, exhibit order/disorder phenomena. High order is described by small entropy and, vice versa. We support the viewpoint that applying the right entropic functional is essential for physical modelling (see [59]). The entropy associated to a BG physical system was proposed by Boltzmann [34] and refined by Gibbs [35] for general systems by the entropy: where the p i denotes the microscopic probabilities of the system occupying a partial region Ω i (1≤ i ≤ N) of the total phase space Ω. The quantity k is a constant. In case of standard Boltzmann-Gibbs thermodynamics this is the Boltzmann constant k = k B . Here, we denote it by k in reference to the generalizations to nonlinear dynamical systems, where it takes other values. These probabilities obey the normalization: For equal probabilities p i = 1/N one derives: When two probabilistic independent subsystems A and B, with numbers of states N A and N B , are put in contact, the joint probabilities obey the special relation: In this case the entropy functional (compare with Equation (1)) is: By applying Equation (2) to the two parentheses in Equation (5), this simplifies to: With the help of Equation (1), it follows that the Boltzmann-Gibbs entropy of a joint system A + B with independent subsystems A and B is additive, viz.: Furthermore, the BG entropy is at its maximum at equal probabilities, shows expansibility to new states and concavity (for proofs of these features see e.g., Tsallis [59,64]). Consider a system in thermal equilibrium with phase space Ω. For a microscopic or at least a very small part of this system, with a phase space element Ω i , with some simplifying assumptions, the probabilities are given by: (see for example [7]). E i is the energy related to the domain Ω i of the system. This expression is called Boltzmann factor. Furthermore, the cited derivation yields: with the Boltzmann constant k, the Kelvin temperature T and the partition function of Boltzmann-Gibbs type, given by:

Kraichnan's BG-Equilibrium Thermodynamics of 2-d and 3-d Turbulent Flow Fields
Now, the question is: can these thermodynamic laws describe turbulence? Kraichnan saw a kind of analogy between a weakly coupled Boson gas below the Bose-Einstein condensation temperature threshold in thermal equilibrium and turbulence. He stated [65]: There is a fairly close dynamical analogy in which the number density and the kinetic energy of the Bosons play the roles of kinetic energy and squared vorticity, which is the specific enstrophy of a two-dimensional turbulent field, see [66,67]. We state that this early introduced example is one of the most instructive ones to show how thermodynamics principally applies to turbulence.
From now on, we assume the fluid to be incompressible. Let us start with the divergence of the momentum equation for such a fluid (in the following we mainly follow [68]): This transforms to: which, owing to incompressibility, div → u = 0, simplifies to the relation: For a divergence-free specific force field, div → f = 0, it follows that: Furthermore, it is evident that: Thus, Equation (14) becomes: 1 Because the Laplace operator of the pressure is expressible as Γ ∆pdV = Γ div (∇p) dV, the integral over a domain Γ, with boundary ∂Γ and for smooth pressure p, can be transformed by the divergence theorem to: in which → n is the unit normal vector on ∂Γ pointing to the outside of the domain Γ. For a Neumann problem, for which: is solved in Γ subject to the flux boundary condition: prescribed on ∂Γ, we regard the right-hand side of Equation (18) as prescribed in Γ and know that this so-called Neumann problem for the Poisson equation is mathematically well posed and solvable. The view point given by Equations (18) and (19) allows us to infer that the pressure p is a quadratic functional of the fluid velocity field → u . Therefore, we have learnt that in an incompressible fluid the pressure p at each instant of time t is fully determined by the total velocity field → u in the (entire) fluid domain and its initial pressure distribution p 0 at time t 0 , where the velocity field takes the form → u 0 , and so: p in which ψ is a differentiable function. We, thus, can define the quadratic operator: so that the momentum equation (see the expression in square brackets of Equation (11)) can be written as: Foias et al. [68] write: This equation is not the usual form of the Navier-Stokes equation, but it puts in evidence the fact that those equations may be regarded as functional evolution equations for the velocity → u . Kraichnan [69] (see also [3,65,[70][71][72]) studied the force-free Euler equations in this sense as: However, he additionally applied a spectral Ritz-Galerkin truncation with a filtering operator that in the following shall be outlined in detail.
We consider the NSE with periodic boundary conditions, e.g., the solutions of the equation in square brackets of Equation (11) with boundary conditions that are spatially periodic and for simplicity Entropy 2018, 20, 109 7 of 41 assume that they have the maximum period L in each of the three directions. Then, a Fourier series expansion of the velocity is: with the dimensionless integer wave number vector → n = (n 1 , n 2 , n 3 ) and the physical wave number vector → k = κ (n 1 , n 2 , n 3 ) in which κ has dimension 1 length . In particular, the following conditions hold: The lowest component of the wave number is: If we now cut the spectrum also at a high wave number: then it follows for the two-fold truncated velocity vector that: For more details on the truncation procedure, cutting away eddies of higher wave numbers than that of a predefined threshold value κ" (n i > N, i = 1, 2, 3), respectively setting their Fourier components equal to zero, see also [70,73,74]. This is tantamount to eliminating all Fourier components with wave numbers larger than κ", a fact that is reminiscent of introducing the Kolmogorov dissipation length l K , if we assume that κ" is the wave number closest or identical to this inverse dissipation length: Galerkin-Ritz spectral procedures project a vector from the entire vector space onto a space spanned by the N first eigenvectors (of the Stokes operator, for a definition using duality see [68]), which is orthonormal: Any velocity vector → u may be written by means of the vectors of an orthonormal basis and the corresponding Fourier componentsû n : As described above, a Ritz Galerkin procedure cuts away the highest modes. In a modification we demand that in the summation (31) it removes both, lowest and largest modes (see Figure 1). A projection is described by: As described above, a Ritz Galerkin procedure cuts away the highest modes. In a modification we demand that in the summation (31) it removes both, lowest and largest modes (see Figure 1). A projection is described by:  Three domains are distinguished: a virtual range that belongs to virtual eddies of diameter larger than the fluid basin (range to the left), the so-called inertial range with the full ensemble of eddies with a diameter L down to l K and the dissipation range with eddies that would have an eddy diameter smaller than the Kolmogorov dissipation scale if they would exist. The circular objects characterize eddies of typical sizes of their domains. The standard Ritz Galerkin truncation neglects wave numbers larger than the bounding wave number value κ".
A Ritz Galerkin truncation can be interpreted as a projection method to convert a continuous operator problem, e.g., an infinite set of differential equations into a discrete problem, in our case a finite set of evolution equations describing a dynamical system (see [75]).
Then, a quadratic series development is introduced (compare similarity with the definition in [70]): where here the terms x i and u i represent independent Cartesian replicas of the space and velocity fields and the terms A imn are spectral coupling coefficients. In this projection operation the constraint: remains valid, because the Galerkin truncation operator P n commutes with derivatives. This is a modified interpretation of Kraichnan's method to account for the fluid incompressibility. It follows that Kraichnan's 2-d and 3-d systems conserve each a subspace, namely in the two-dimensional problem defined by the 2-d kinetic energy per unit of surface mass ρ (kg m −2 ): and the enstrophy: whereas, in the three-dimensional case, it is the 3-d kinetic energy: Entropy 2018, 20, 109 9 of 41 and the helicity: For brevity reasons, we introduced a different notation, namely the substitutions given by Ω ε = ε, It is clear that if κ" tends to infinity that a transformation from the discrete to the continuous case, can be performed by replacing sums by corresponding integrals.

Remark 1.
The enstrophy definition which in Equation (36) has been introduced is a special case where incompressibility is assumed to hold. The more general definition of enstrophy is: where |. . .| denotes the Frobenius norm, which is the sum of the absolute values of all the Jacobian matrix elements J i,j = ∂u i /dx j , i ∈ {1, 2, 3}. Now, following Kraichnan and going back to the discrete description, one distinguishes: (i) Energy conservation: with d = 2 for χ = ε and d = 3 for χ = H (we wrote "∝" to neglect the constant ρ/2 in this and the following equations). With Equations (30) and (32) one derives: Note the difference in the summation indices of Equations (41)! (ii) Enstrophy conservation: Advantageous is that the operations "taking the derivatives and applying the Ritz-Galerkin projection" commute.
We start with Equation (24) and in the following auxiliary calculation we omit for brevity reasons the tilde superscript in the velocity components: Furthermore, denoting conjugate complex numbers by overbars, it follows (using Einstein's summation convention) that: With the discrete version of Equations (39) and (44), it is concluded that: (45) or: Next, we again apply the Ritz-Galerkin truncation and Equations (44)- (46) and obtain: With the analogous procedure as described by Equations (40) and (41), we can now derive for the enstrophy the final representation: where k = → k denotes the modulus of the wave number → k .
(iii) Helicity conservation: Because in this article we are mainly concerned with 2-d turbulence, we shall not prove representations (49). The 3-d case shall be treated in future work elsewhere.
Notice that all these conserved quantities are quadratic forms of the velocity vector. However, only for the Euler equation they are strictly conserved. Then, for the discretized version, Kraichnan introduced a set of probability densities by: with the partition function Z BGχ . The generalized energy E χ , and its average quantity < E χ >, are now defined as: Because in Equations (40), (46), (48) and (49), E χ and Ω χ are quadratic functions of u i orû n , we know that this system of evolution equations, constituting a conservative dynamical system, is obeying Gaussian statistics.
A Fourier transform of E ε for 2-d turbulence leads to the following generalized energy-enstrophy spectrum (see Refs. [3,70] and for its derivation the Appendix A, Equation (A34)): where α and β are constants. On the other hand, a Fourier transform of E H for 3-d turbulence leads to the energy spectrum (see also [3,70]): in whichˆ E H (k) is the Fourier transform of E H (k). Kraichnan calls expressions in Equations (52) and (53) the 'mean modal intensity spectra'. Then, the 'isotropic energy and enstrophy spectrum' of 2-d turbulence follows by Equations (54). The second formula is only valid, if one is dealing with power laws: In the 2-d situation Equation (52) implies: By applying Equation (54), the energy and enstrophy spectra for the low and high wave number regimes are: Whereas the energy spectrum shows the useless value zero, the enstrophy spectrum is in excellent agreement with computer experiments (see Section 8 below). This suggests that at high wave numbers the occurring small eddies are close to or even in thermodynamic equilibrium. However, the experimentators state that their simulation results may be slightly erroneous because of aliasing and slow convergence of their calculations. This is motivation for us to also q-generalize the results for the large wave number regime. In future, if even higher precision experiments would be available, the power law exponent of this spectrum could be adjusted with a small non-equilibrium contribution. In any case q values very close to "1", respectively power law exponents close to "−3" are expected to occur. Kraichnan called the high wave number regime 'inverse cascade' of two-dimensional turbulence ( [65,76,77]).
We ask ourselves: What is the essential conclusion of this instructive spectral example, derived by Kraichnan fifty years ago? The answer to this question is given by: Proposition 1. By a Ritz-Galerkin truncation of the 2-d and 3-d Euler equations, these equations can be transformed to a finite system of evolution equations with two conserved quadratic forms that obey Gaussian statistics.
Furthermore, we have corroborated and now state: Truncating the momentum equation and applying Boltzmann-Gibbs equilibrium statistical mechanics (thermodynamics) leads to a loss of essential features of low-wave number turbulence. It is a too bold approach of modeling turbulent fields, like the application of a linear and local closure scheme (see Table 1). Today it is known that the entropy functional and the partition function, respectively, define 'universality classes' showing identical 'critical exponents' (see [31]). This means that a system close to its criticality exhibits a kind of universal behavior that is not determined by the specific microscopic interactions of neighboring particles. It was shown, for instance, that a liquid-gas transition and a ferromagnetic system (lattice gas and Ising model), in the region just below criticality, possess with high precision the same scaling exponents, because they belong to the same universality class (Refs. [29,30]). This agrees with the statement of Goldenfeld [27] that phenomena with the same set of critical exponents are said to form a universality class. Therefore, one ought to distinguish between exact equivalences, where two models are analogous to one another (with a unique mapping from one to the other) and approximate equivalences, in which the partition functions may be (slightly) different, but nevertheless comprise two systems that show the same or similar behavior near criticality. According to Goldenfeld [27] members of the same universality class have the following three properties in common: (1) The symmetry group of the Hamiltonian, (2) The dimensionality of the physical problem, (3) The range of the forces (short or long).
With these facts, we assert that the Boltzmann-Gibbs thermodynamics is inadequate to describe turbulence sufficiently accurately. In nonlinear dynamics (see [53]) different approaches have been worked out, which are also presented in Tsallis [59]. They involve different entropy functionals, as e.g., the Normalized entropy, Shannon entropy, Rényi entropy, Kolmogorov-Sinai entropy, Escort entropy, etc. Up-to-present, not all questions have been answered concerning what kind of generalized entropy functional is optimal and which one best applies to a specific complex system.

An Introduction to the Extended Thermodynamics of Tsallis
Below we demonstrate that yet a further entropy functional, called Tsallis entropy, is particularly optimal to describe turbulence. This entropy can be directly related to fractal geometry and Lévy flight statistics (see [78,79]). Furthermore, the above class of entropies have many features in common and the Tsallis entropy enjoys simple relations with other entropies listed e.g., in [59].
Let us introduce the Tsallis entropy in a metaphoric manner (see [59]). To this end, we aim for the derivation of Equation (3) by considering the linear differential equation: in which, for simplicity, we concentrate on the case of a single probability, i.e., N equals "1"; for this case: Writing (57) as: its integration yields: in which log e C is the constant of integration. It follows that: With Equation (60) the constant C is determined by expressing the requirement of certainty (which is equivalent to S BG (1) = 0) as: implying C = 1, so that, owing to Equation (60): which is also the entropy of a system with elements having equal probabilities that was already anticipated by Equation (3). With Equations (58) and (63) this is generalized to become: (compare with Equation (1) for a single state p) for a system of N elements with the single element probability p. Furthermore, Equation (61) reduces to: At this stage the differential equation with its initial condition (57), is suggestive for a generalization of the Boltzmann-Gibbs entropy functional to nonlinearity. Following this recipe, Tsallis proposed a nonlinear scaling with value q ∈ R of the linear differential Equation (57), which is: Writing (66) alternatively as: a straightforward integration yields the result: where C is a constant of integration, which is determined to be C = 1, since for N = 1, S q = 0. Thus, the general integral of (66) is: Now, we shall define a generalization of the exponential function, called 'q-exponential function' (suggested by Tsallis [59]) and given by Equation (69): This function is displayed in Figure 2 for different values of q. In agreement with relation (70) and n = 1/(1 − q) the well-known formula, defining the usual exponential function, emerges if the limit q → 1 of Equation (70) is taken, viz.: and n = 1/(1 − q) the well-known formula, defining the usual exponential function, emerges if the limit 1 → q of Equation (70) is taken, viz.: Kolmogorv-Oboukov function Figure 2. The q-exponential function x q e for typical values of q. After Tsallis [59], reproduced with additions, e.g, the Kolmogorov-Oboukov function (see Section 5).
From Equations (69) and (71), we draw the following inference: For q = 1, the BG-entropy S BG and the Tsallis entropy S q are the same.
On the other hand, the inverse function of (69) is: This now suggests the definition of the 'q-logarithmic function': where here, log e x = ln x is used to avoid a confusion between the q exponent and the basis b of the logarithmic function log b. The definition (74) also implies: Moreover, in agreement with Equation (74) and the assignment n = 1/(1 − q) the formula defining the usual logarithmic function is deduced from the following chain: From the above presented metaphor and Equations (73) and (74), Tsallis, for equal probabilities, postulated the following generalization of Equation (3): In Figure 3 this function is displayed for k = 1 and different q-values.
From the above presented metaphor and Equations (73) and (74), Tsallis, for equal probabilities, postulated the following generalization of Equation (3): (77) In Figure 3 this function is displayed for k = 1 and different q-values.  With p = 1/N, q S in Equation (77) can be expressed as a function of p: With p = 1/N, S q in Equation (77) can be expressed as a function of p: Notice that it must be carefully checked which rules of the usual logarithm apply to its generalized form. The relation ln q (1/p) is not identical to −ln q (p) unless q = 1; thus, we cannot write this equation in the q-generalized form of Equation (1) for a single probability. The next step is to demand: In this generalization process, arbitrary probabilities must be introduced as: In agreement with Equation (76), we conclude that, according to (80) and (1) for N = 1: and, therefore, as already stated after Equation (72), equilibrium thermodynamics of the Boltzmann-Gibbs type is a special case of the, thus, proposed Tsallis non-equilibrium thermodynamics. Inserting definition (75) of the q-logarithmic function into (80) leads to:  (2), is identical to: This q-entropy or Tsallis entropy generalizes the entropy of an equilibrium to that of a non-equilibrium thermodynamic system, respectively from systems with linear to those with nonlinear and complex behavior; the latter is suitable for turbulence. Notice that this entropy initially was introduced by ideas of fractal geometry. It follows that: Next, a new set of probabilities is introduced by defining, with index ES denoting "Escort": The complete set {P i } is called 'escort probabilities'. They obviously also fulfill the normalization condition: The condition q < 1 makes the escort probabilities larger than the p i s, and q > 1 makes them smaller. Therefore, q < 1 enhances the occurrence of rare or large-scale events and q > 1 enhances those of frequent or small-scale events. We experience here a connection with Lévy statistics in which the case q < 5/3 corresponds to sub diffusivity, q = 5/3 describes the standard Brown or normal diffusion, obeying Gaussian statistics, and the case q > 5/3 belongs to super diffusivity (see Tables 2 and 3).
We will identify more quantitative relations between Lévy statistics, anomalous diffusion and extensive thermodynamics in Section 5. When q < 1, it is important that summations must exclude all probabilities of zero value (because of singular behavior), whereas for q > 1 the summation may include all occurring probabilities.
Next, let us investigate the compositional property of the Tsallis entropy. Then, from Equation (74), we conclude that: This is identical to the following extended mathematical expression: which with use of Equation (87) is now rewritten to take the form: With the help of Equations (73) and (74) these expressions are transformed back to Tsallis entropies. This yields: which demonstrates non-additivity of S q for q = 1. Depending on q, it follows that: Therefore, for q < 0, S q (p) is a convex function, whereas for equal probabilities, p 1 = p = 1/2 and p 2 = 1 − p = 1/2, it takes a minimum (Figure 4, left). On the other hand, for q > 0 it is concave with a maximum for equal probabilities in this case (Figure 4, right). Slowly, it becomes transparent that non-extensive thermodynamics is an appropriate framework to tackle anomalous diffusion and, thus, also turbulence (see [80]).  We stress that all q entropies are strictly convex or concave, a feature that not all generalized entropies (e.g., Renyi, Escort, etc.) fulfil. After Tsallis [59] reproduced with additions. e.g, the Kolmogorov-Oboukov function (see Section 5).

Relation between Lévy Statistics and Tsallis Extended Thermodynamics
To make the entire situation more transparent, we now outline the connections between Lévy walk and Lévy flight dynamics (see [80][81][82]), respectively, and the Tsallis non-extensive and nonadditive thermodynamics (see [59,64]).
The mean square displacement of diffusive problems is generalized in [80]: We stress that all q entropies are strictly convex or concave, a feature that not all generalized entropies (e.g., Renyi, Escort, etc.) fulfil. After Tsallis [59] reproduced with additions. e.g., the Kolmogorov-Oboukov function (see Section 5).

Relation between Lévy Statistics and Tsallis Extended Thermodynamics
To make the entire situation more transparent, we now outline the connections between Lévy walk and Lévy flight dynamics (see [80][81][82]), respectively, and the Tsallis non-extensive and non-additive thermodynamics (see [59,64]).
The mean square displacement of diffusive problems is generalized in [80]: with time t and exponent α = 2H = 1, describing Brownian diffusion. H denotes the so-called 'Hurst exponent'. After Alemany and Zanette [81] the jump probability distribution can be obtained from Equation (1) by maximizing the entropy S BG , in the continuous case written as an integral: subject to the two auxiliary constraints, namely: and: the quantity σ is the standard deviation. For simplicity, we have chosen the 1-d representation.
For higher dimensions the considerations are analogous (see [81]). The problem of maximizing the entropy (93), subject to the constraints (94) and (95) is technically performed by maximizing the Lagrange function: in which λ 1 and λ 2 are the Lagrange parameters. The solution of this extremum principle determines these Lagrange parameters as well as the jump probability distribution, which takes the form: where for details it is suggested to consult e.g., the book on fluctuation theory by Montroll and Lebowitz [82]. The proof of Equation (97) is given as a special case of the proof of Equation (103) (paragraph below and [81]). A Fourier transform of Equation (97) for small wave number k (large wave length λ) leads to the characteristic function (see [81]) of the following form: With innovative imagination, this is now generalized as:  (97) and (98)). An inverse Fourier transform of this equation reveals the asymptotic Lévy distribution appropriate for large distances, with the scaling property: This is called the long tail of the Lévy probability distribution. To be precise, the approximate large distance (small k) Fourier transform of (99) is (100), a result of the linear equilibrium theory, which is then made nonlinear in (99), and this approximation in Fourier space is assumed to be representative for the long-tail behaviour of the Lévy probability distribution expressed as Equation (100).
The same idea of applying an extremum principle and by taking into account generalizations of the constraints (94) and (95), replacing the probability p by the escort probability P, viz.: with the quantity σ q playing the role of a q-generalized standard deviation and with the q-entropy S q (see Equation (83)), viz.: leads to the result (see [81]), with the generalized Boltzmann constant k: and the two Lagrange multipliers λ 1 and λ 2 . The applied method is outlined in more detail in the next section and the proof of Equation (103) is given in [79]). Furthermore, more information on the Lagrange multipliers is given in Section 10. For large distances x, setting the exponents of the formulas Equations (100) and (103) equal to one another, yields: These two relations directly connect Lévy flight statistics with the thermodynamic formalism of Tsallis. The corresponding characteristic parameters are presented in Table 3.

Escort Probability Distribution and Expectation Values
The entropy optimization method, as applied by Kraichnan in 1967 (see [70]) for BG thermodynamics, shall now be analogously applied to the Tsallis thermodynamic theory, because it is a prerequisite to perform thermodynamics with generalized entropic forms. With some small deviations, we follow mainly this reference.
Let us now, in analogy to Equation (85), introduce the continuous escort probability distribution: in which the integral is introduced as a Lebesque integral with continuous measure. The distribution is also normalized, viz.: Furthermore, by defining the q-average value of a function f as: it follows for f = x, that: which is the q-average value of x.
To optimize the entropy distribution S q [P], in analogy to Equation (96), we define the Lagrange function Φ, containing the two constraints Equations (106) and (107) for the energy function given by f (x) = E(x): Substituting the q entropy, Equation (102), yields: in which E q is the average energy defined below (see Equation (115)) and where the denominator q − 1 can be taken as a quantitative measure of "Non-BG behavior". Maximization is obtained by performing the optimization dΦ[P]/dP = 0; it leads to the optimum result P max , where, upon dropping the index, one obtains (see [59]): The optimization process determines the generalized Boltzmann factor P [E(x)], which is called here Tsallis factor. One can show that the Lagrange multiplier λ 1 , due to Equation (106), factorizes out (see [59]). Now, a q-generalized partition function is defined by: When applying the q-exponential function (70) and definition (112), the Tsallis factor (111) is found to be given by (see also in [79]): Next, the average energy is calculated by:

Fractional Calculus: A Promising Method to Describe Turbulence
In analogy to Kraichnan's calculus of the energy-enstrophy spectrum, to obtain the q-generalized versions of Equations (52), (55) and (56), we require fractional calculus. To be able to generalize the steps between these equations, the Fourier transform of a derivative and its square must be calculated. Therefore, we deal in a general manner with fractional derivatives. The reason is that the common Riemannian derivative, in our case, insufficiently describes turbulent phenomena. The adequate new tools for our problem are fractional derivatives and integrals. Such were, for example, also applied with success by the present authors to turbulence modeling of turbulent shear flows (see e.g., [26]).
In Fourier calculus the following main principle holds: where F denotes the Fourier transform. Now, we study the special case n = 0, which yields: We observe that by a Fourier transform a derivative of order m of a function f(x) is transformed to a power law of the wave number k, with exponent m, times the Fourier transformed function. This is the basic rule to generalize the usual derivative to the Fourier fractional derivative of order m = q: The rule for usual integer-order derivatives and power laws with integer exponents, is generalized in fractional calculus to real and complex numbers (see e.g., Refs. [83,84]); in our context real-order derivatives are sufficient. By applying the Fourier fractional derivative to a power law, one finds (see e.g., [84]): Now, we set q = 1, which yields: With the formula (see [85]): which is the usual Riemannian differentiation rule for a power law with exponent α. In fractional calculus there is no consensus on the differentiation of a constant. There are essentially two definitions: where definition (123) is due to Caputo's fractional derivative [84] and Equation (123) is the solution for a Riemannian fractional derivative [84] with a constant c being a function of Gamma functions and is motivated by Equation (119) with α = 0. As the first definition leads also to the useless result (56, left), derived with BG thermodynamics, the second one is very promising and, therefore, leads us in a slightly pragmatic way to give it preference (see Section 8).

Fractional Generalization of Kraichnan's Spectra and Their Validation by Numerical Experiments
With such knowledge, Equation (54) is now q-generalized in the following manner: Now, a Fourier transform of E ε for 2-d turbulence leads to the following generalized energy-enstrophy spectrum (its proof is outlined in the Appendix A, see especially Equation (A33)): Two wave number regions emerge from this formula. The first is the small wave number regime, for which, from Equation (125), it follows that: Now, by applying Equation (123), one calculates the fractional derivative of this constant: which is different from Equation (56). Notice that in Equation (127), with q = 1, a reciprocal energy spectrum does not result, because in this case c = 0 will guarantee agreement with the usual differrentiation, so that (56) also occurs as special case of this new treatment. This would be a wave number spectrum related to the thermodynamic equilibrium of a turbulent flow field at small wave numbers, which does not exist! However, as we can see in the third line of Table 4, there is another special case, namely the one corresponding to Brownian diffusion, which is characterized by q = 5/3, leading to the well verified Kolmogorov-Oboukov energy spectrum: This result is in best agreement with computer experiments for 2-d turbulence performed and reported by Lilly [86] and reviewed by Kraichnan [70]. Let us summarize: Proposition 3. Tsallis' extended thermodynamics generalizes Kraichnan's constant turbulent energy by a power-law energy spectrum with exponent −q, which relates this generalized spectrum directly with the Kolmogorov-Oboukov −5/3 law and its fractal generalization, e.g., given by the fractal β-model (see e.g., Frisch [3]), describing turbulent flow fields with intermittency (see below). With Kolmogorov's assumption that the energy spectrum depends only on the wave number modulus k and the dissipation rate ε, the exponent −5/3 power law is obtained for two as well as for three dimensional turbulence (see [87] and Figure 5). Notice that this spectrum is valid only at low wave numbers up to an intersecting value k I after which a k −3 power law (see also Figure 5) is observed that dominates the behaviour up to Kolmogorov's dissipation wave number ( [70,87]). This result is in best agreement with computer experiments for 2-d turbulence performed and reported by Lilly [86] and reviewed by Kraichnan [70]. Let us summarize: Proposition 3. Tsallis' extended thermodynamics generalizes Kraichnan's constant turbulent energy by a power-law energy spectrum with exponent −q, which relates this generalized spectrum directly with the Kolmogorov-Oboukov −5/3 law and its fractal generalization, e.g., given by the fractal β-model (see e.g., Frisch [3]), describing turbulent flow fields with intermittency (see below). With Kolmogorov's assumption that the energy spectrum depends only on the wave number modulus k and the dissipation rate ε, the exponent −5/3 power law is obtained for two as well as for three dimensional turbulence (see [87] and Figure 5). Notice that this spectrum is valid only at low wave numbers up to an intersecting value I k after which a 3 − k power law (see also Figure 5) is observed that dominates the behaviour up to Kolmogorov's dissipation wave number ( [70,87]).

Figure 5.
For comparisons of the power laws describing the spectra of the inertial energy and the enstrophy range early computer simulations of Lilly [86] have already been cited and interpreted by Kraichnan [70]. These numerical experiments show a convergence toward the predicted values of the power law exponents, which for the highest computation time (2200 time steps) are at low wave number close to -5/3 and at high wave number close to the equilibrium value -3. From Lilly [86] with changes.
Furthermore, notice also that a Kolmogorov-Oboukov turbulent flow is thermodynamically equally a non-equilibrium process that, on the other hand, shows simplicity by its lack of intermittency. For ∞ → Re , turbulent fields converge toward a state obeying the Kolmogorov- Figure 5. For comparisons of the power laws describing the spectra of the inertial energy and the enstrophy range early computer simulations of Lilly [86] have already been cited and interpreted by Kraichnan [70]. These numerical experiments show a convergence toward the predicted values of the power law exponents, which for the highest computation time (2200 time steps) are at low wave number close to −5/3 and at high wave number close to the equilibrium value −3. From Lilly [86] with changes.
Furthermore, notice also that a Kolmogorov-Oboukov turbulent flow is thermodynamically equally a non-equilibrium process that, on the other hand, shows simplicity by its lack of intermittency. For Re → ∞, turbulent fields converge toward a state obeying the Kolmogorov-Oboukov energy spectrum. Tsallis [59] wrote that chaotic systems with strong chaos (showing at least one positive Lyapunov exponent) approach thermodynamic equilibrium.
Weakly chaotic and turbulent systems show complexer behavior than infinite Reynolds number turbulent flows. These systems are closer to criticality. Table 4, in its last line, shows that turbulent systems usually are characterized by q values between 5/3 and 3. Therefore, with (127) also the next refined energy spectrum follows: The intermittency exponent µ is closely related to the structure function ς p of order p of the fluid velocity (see [3]). From the q-value range in the last line of Table 4, we also deduce that: and conclude that Tsallis' q factor copes directly with basic quantities of research on turbulence, which are known for many decades. Note that a power law exponent with a magnitude larger than 5/3 was a hypothesis of Kolmogorov [88] and Oboukov [89] for an intermittency correction (see also in [46]), which by the present analysis is more than just qualitatively confirmed. A comparison of Equation (130) with a corresponding formula in [22] reveals Equation (130). In this equation d denotes the spatial dimension and D the Hausdorff-Besicovitch dimension, which is a fractal dimension. The difference d-D is called codimension. All this is in order as long as one works with Lévy flight statistics (see e.g., [90][91][92]), where a cascade of Lévy flight sizes directly relates to corresponding diameters of eddies in a way that basic kinematic laws of turning eddies with different angular velocities are fulfilled (see [22]). This is also in full agreement with results of the fractal β-model [3]. On the other hand, Gotoh and Kraichnan [90] critically comment on the application of Tsallis' non-equilibrium thermodynamics (with a constant q-factor) to turbulence. Arimitsu and Arimitsu [93] make the link of Tsallis thermodynamics to the multifractal model of turbulence. In this model the intermittency and, thereby, its exponent µ are scale dependent. Notice that intermittency increases with decreasing scale. These considerations are in perfect agreement with the direct relation (130) between µ and q, where a scale dependent intermittency exponent also calls for a scale-dependent q factor. In the language of Tsallis thermodynamics, different wave number sub systems show different q values.
The second is the large wave number regime, for which, from Equation (125), it follows that: Now, by applying Equation (119), the fractional derivative of this function is derived, which yields: which for q = 1 is identical with Kraichnan's former result (56), derived with BG equilibrium thermodynamics, showing an 'exponent minus three power law'. Also this result is in good agreement with the numerical experiments of Lilly [86], presented in Figure 5. The good coincidence between the equilibrium thermodynamic theory of Kraichnan, being the Tsallis thermodynamics special case with q = 1, and the experiments suggests that small eddies, which have smaller turnover and life times compared to large eddies (see Section 9) are close to or even exactly in thermal equilibrium. If in future experimentally determined small deviations to the power law exponent "−3" would be observed, then, with high probability they would be due to small nonequilibrium effects. In this context, note that Kraichnan proposed a small logarithmic correction to the enstrophy range of the spectrum [70].

Justification of the Quadratic Form of the Energy as a Functional of Real Space Coordinates
In this section we give a motivation of the quadratic form of the energy functional, which Kraichnan discovered. From the constancy of the real space and velocity variables (see below), viz.: also it follows the constancy of its sum (for the dimension of → x see explanation following Equation (160)). This explains the statement of Kraichnan [70] that the constant energy surfaces in phase space are hyperspheres. To describe the eddy fluctuation motion, we start with a 2-d coupled dynamical system .
where for simplicity we neglected all the wavy overbars. Differentiation of (Equation (134), left) leads to ..
By differentiating (Equation (134), right) and substituting (Equation (134), left) a second second-order differential equation follows, so that we now face a system of two ordinary differential equations .. .. (137) These are differential equations of harmonic oscillators for the two coordinates x 1 and x 2 with the solutions: The squared velocity is: Substituting Equations (134), it follows that: By a substitution of solution (138) into (139) and the application of Pythagora's law for the trigonometric functions the same result is obtained. The energy of this rotating eddy is: In a usual harmonic oscillator, the total energy E is the sum of the kinetic energy T and the potential energy U, which is constant. In the oscillation process these two exchange their energies under the restriction: On the other hand, in a turbulent eddy there is an exchange between two kinetic energies of the two coupled oscillators: We now rewrite this in the notation of the fractal beta model (see [3]) by renaming r by l/2 and assigning the index 1 to several quantities: Notice that in turbulence research the largest eddy is denoted by the index 0. To cope with Kraichnan's notation, we start the summations with index 1. Kraichnan stated that the total energy E is constant: Now, we take advantage of the fractal-β model (see [3]). Thereby, we follow the presentation in [22]. In this model eddies are assumed to occur in classes of different sizes l n (see Figure 6). The largest eddies show a diameter l 1 . The ratio of eddy sizes of two neighboring classes is: where, following an assumption by Kolmogorov, in Figure 6 we have chosen a decay of eddies to such of half size (b = 2). Notice that b is an important parameter of a Lévy flight distribution, which in a special fractal model is connected with the eddy size ratio (see [22]). Now, it follows that: In the fractal β-model eddies do not fully occupy the available space. In each step from one eddy class to the next a reduction by the factor β (0 < β < 1) occurs, where this β gave the model its name. Then, the active space fraction is: Here, we do not need to specify the spatial dimension. However, in 2-d turbulence one could imagine that the excessively large eddies are perfectly two dimensional, whereas finer structures (smaller eddies) fluctuate in all three spatial directions. It is assumed that at some intermediate wave number a small transition region occurs. A first assumption is that the eddies of the first class fully fill the available space. From Equation (148) it follows that: Then an increasing intermittency to smaller scales is observed. The kinetic fluctuation energy in an eddy of class n is: From this it follows that the velocity of eddies of class n is given by: The time in which the energy in an eddy of class n is transferred to one of class n + 1 is assumed to be the time of a single turn of an eddy of class n; consequently, this time is called 'turnover time' of an eddy: For further remarks on t n see Ref. [22]. The locality hypothesis of Kolmogorov [87] states that the transfer rate of energy from an eddy of one to its next class is constant and given by: The range in which this is the case is called the inertial range. Now, it follows that: where Equation (152) was substituted into (154). From this it follows that: For the largest eddies one has: where by writing in Equation (156) with intention twice l 1 with their specific exponents each, simplifies the following calculus. Now, by substituting Equation (156) into (155) yields: A comparison of this exponent with Equation (130) reveals that: with the Tsallis parameter q. Now, Equation (157) is rewritten: With Equation (147) it follows from (159) that: To obtain the total energy of the eddies of all classes, this equation is inserted into the sum of Equation (145): The final result is: This deserves the formal: Proposition 4. By scaling laws the total turbulent kinetic energy can be expressed as a function of the energy contained in the eddies of the largest size E 1 only. The turbulent kinetic energy of the eddies of the remaining classes E 2 . . . E N is described by the multiplicative factor κ − 1 times E 1 . This is a direct consequence of the self-similarity of the eddy cascade. The special case N = 1 leads to the correct solution of a total energy belonging to that of the single class containing only largest eddies: The special case q = 1 yields a division "0/0". Therefore, to Equation (162) the rule of 'Bernoullide l'Hôpital' is applied by differentiating the nominator and denominator with respect to b, yielding: which leads to the energy equipartition law between the eddies of different classes: (165) For emphasis, the results (160) and (165) suggest:

Proposition 5. Equilibrium BG-turbulent flows show energy equipartition in the eddies of all classes and, consequently, for infinite Reynolds number flow its turbulent kinetic energy diverges.
Next, with Equations (144) and (162) takes the form: By modifying the coordinates by the stretching transformation that this space coordinate has the dimension of a velocity-one finds: which is the result that also confirms Equation (A8) in the Appendix A below, q.e.d. The same results, with a different constant κ, could also be derived by taking in Equation (162) the occupation probability into account. This type of modeling includes the birth rates and life times of eddies of different classes (for details see [22]).  The special case N = 1 leads to the correct solution of a total energy belonging to that of the single class containing only largest eddies: The special case q = 1 yields a division "0/0". Therefore, to Equation (162) the rule of 'Bernoulli-de l'Hôpital' is applied by differentiating the nominator and denominator with respect to b, yielding: which leads to the energy equipartition law between the eddies of different classes: For emphasis, the results (160) and (165) suggest: Proposition 5. Equilibrium BG-turbulent flows show energy equipartition in the eddies of all classes and, consequently, for infinite Reynolds number flow its turbulent kinetic energy diverges.
Next, with Equations (144) and (162) takes the form: By modifying the coordinates by the stretching transformation , 2}-note that this space coordinate has the dimension of a velocity-one finds: which is the result that also confirms Equation (A8) in the Appendix A below, q.e.d. The same results, with a different constant κ, could also be derived by taking in Equation (162) the occupation probability into account. This type of modeling includes the birth rates and life times of eddies of different classes (for details see [22]).

Proposition 6.
By applying a Ritz-Galerkin truncation the continuous spectrum of turbulent flows is discretized. It belongs to a system of 2N coupled oscillators with N pairs each with its own fluctuation frequency. The large number of frequencies ω n , n ∈ {1, . . . N} scale in a self-similar manner.
A description by periodic or quasi-periodic motion finally rests on a somewhat crude model of turbulence. However, this does not shed any remarkable shadow on this kind of valuable description of isotropic turbulence. Rather, Kraichnan's model and its generalization in this article offers considerable deep insights into the physics of 2-d turbulent flows, as it will also become clear in the next section.

The Lévy Flight Probability Distribution and the Weakly q-Deformed Gaussian Distribution
In Section 5 we were working with the Non-Gaussian Lévy flight probability distribution (103). On the other hand, in Section 8 the proof of Equation (125) was performed with a weakly q-deformed Gaussian probability distribution (see insert in Equation (A3) below). Therefore, the question could be raised, whether this might be an inconsistency in our modeling. However, in this section profound explanations are given on these two distributions, and it is shown that they are so closely related that in the present physical modeling no claim is justified that the procedure is deprived by rationality.
We start by clarifying what we define as a weakly q-deformed Gaussian distribution and what was already earlier defined to be a strongly or just a q-deformed Gaussian distribution.
Firstly, a weakly q-deformed Gaussian distribution is proposed by writing a usual Gaussian probability distribution as an escort probability distribution by combining Equation (97) with (105): Then, it follows that: As demanded, the case q = 1 leads to the original classical Gaussian probability distribution (97). Definition (169) was applied in Equation (A4), and the motivation for this choice is given below.
Secondly, according to e.g., [94,95], a q-deformed Gaussian distribution is defined by applying Tsallis' generalized exponential function (70): where also here, based on considerations in Section 4, for q = 1 the classical Gaussian distribution (97) follows. The generalized standard deviation σ q is defined by Equation (101). It is Kraichnan and his co-author Gotoh [90] who give us confirmation that the weakly deformed q-Gaussian probability distribution is the right choice in the context of our work. They state that with the Ritz-Galerkin description of the truncated dynamical system and the three constraints: in phase space, an energy form exists that is a hyper sphere. Then also when q = 1, the energy shell shows a homogenous probability density, and if additionally, N → ∞, the probability distribution converges toward the Gaussian distribution: Note, that we assume that our finite dimensional dynamical system contains a very high number (2N) of binary oscillators, so that in equations following (171), also for finite N, we will write instead of the symbol "≈" the sign "=". Furthermore, we know that in turbulence dissipation occurs, an effect which destroys the Liouville property c). On the other hand, in a delicate equilibrium between a random forcing and dissipation, which e.g., can be described by a (fractional) Langevin equation, Tsallis thermodynamics can be extended to weakly q-deformed Gaussian distributions as applied in our proof in the Appendix A, which is motivated by Gotoh and Kraichnan's arguments above.
Next, with the statement b) above, we generalize with help of (105) distribution (171) as: in which tildes have been dropped and where we have written x 2 for → x 2 . Our approach is based on the escort generalized Boltzmann probability distribution and yields: By inserting the energy (51) (see also Equation (A9) below), it follows that: Comparing the prefactors of Equations (172) and (174) yields: and by comparing the arguments of the exponential functions, with u 2 = 2E, yields: which, with Z ESε = 2β q / q 2 N , is consistent with Equation (125). Notice that the constant in Equation (176) depends on the space dilatation or contraction given in front of Equation (167)! Next, Equation (172) is rewritten as: Assuming that the argument in the exponential function in (177) is an order of magnitude smalller than unity, a linear Taylor series expansion is applied to yield: In the next step the prefactor on the left-hand side of (178) is incorporated into the squared brackets; this yields: . (179) The first term in the parenthesis is compared with the first one of Equation (103). This process yields: The same procedure for the second terms of Equations (103) and (179) and solving the emerging equation for λ 2 leads to: Inserting (180) results in: Next, we consider the BG special case. Inserting q = 1 into Equation (180) delivers: which corresponds to a diverging first Lagrange parameter. The limit q → 1 applied to Equation (182) reveals: and when we apply Equation (165): Next, we introduce a 'generalized temperature' of turbulent flows. We assume that eddies of class n with energy E n are characterized by their generalized temperature T n . Thermodynamics suggests that the relation between energy and generalized temperature is: We may call T n the generalized turbulent Kelvin temperature of an eddy of class n. Applying this relation with n = 1 to Equation (185) yields: This leads us to: Proposition 7. In turbulent BG flows the inverse second Lagrange parameter λ 2 defines the generalized temperature of the eddies being members of the first class.
Next, we apply the energy equipartition law of equilibrium turbulent flows (see Equation (160)), viz.: and conclude with help of (160) and (186) for the entire energy-enstrophy spectrum that: By applying this new result to 2-d turbulence, we now create: These findings are consistent with Kraichnan's statement that in the enstrophy range, where eddies belonging to different wave number intervals are in thermal equilibrium, there is no energy transfer (see [69,70,86]).
Finally, it is stated that the probability distribution (103) and the escort generalized weak Gaussian probability distribution (172), applied in our model, do not lead to an inconsistency, because Equation (103) can be interpreted to be a first-order Taylor expansion of Equation (172). Equally satisfying is that comparisons between the two probability distributions lead to important correct statements on the equilibrium and non-equilibrium thermodynamics of eddies of different classes.

Discussion of Results, Conclusions and Outlook
This review is mainly based on work of Kraichnan [69,70] and Tsallis [59,64,78] and our extensions. A further important contribution, namely the relation between Tsallis thermodynamics and Lévy statistics is due to Alemany and Zanette [81]. Our paper incorporates and combines these works, but also provides a majority of new results; some remained unanswered and were hanging for decades as open ends. In this memoir a unique modern theory of 2-d turbulence was created by e.g., implementing new ideas of fractional calculus, so that our new theory now presents itself with completeness and shows perfect agreement with experiments. We have also found a physically clearly defined generalized temperature of turbulent flows. Its variation acts as a measure of deviation of a turbulent system from equilibrium.
Kraichnan applied equilibrium BG thermodynamics to 2-d and 3-d turbulence. For 2-d turbulence two partial regimes are observed, a low wave number regime, where energy is transferred down the cascade, and a large wave number regime in which the transported physical quantity is the enstrophy [96], flowing toward small wave numbers. The generalization of the 3-d turbulence energy spectrum with Tsallis nonextensive formalism and fractional calculus seems promising, but has not yet been derived and is proposed for further investigation.
On the other hand, Tsallis has been a forerunner of generalizing thermodynamics to non-equilibrium processes. Turbulence at moderate Reynolds number is a highly non-equilibrium physical phenomenon and demands new nonlocal concepts that are still in their infancy. A large selection of different entropic forms was proposed, which luckily have some mutual connections, so that it is often more a matter of taste than correctness to apply one or another. The Tsallis entropic form works with a new formalism that generalizes the significant classical statistical based terms by a q-factor and, thus, leaves well-established formulas of equilibrium thermodynamics preserved and in a generalization to non-equilibrium thermodynamics practically unaltered, if the basic functions are q-generalized (e.g., by the q-exponential function, the q-logarithmic function or the q-deformed fractional derivative [54]). We introduced the weakly q-deformed Escort-Gaussian probability distribution, which is the right choice for a Ritz-Galerkin truncated system of the NSE that develops to a high-degree dynamical coupled system of 2N oscillators or N eddies, respectively. For a large number N of eddy classes, the Gaussian distribution is a good approximation and in the limit as N → ∞ it is even the correct one. We found that the Non-Gaussian Lévy probability distribution can be sought to be a Taylor expansion to this limiting case. However, it is well-known that most complex behavior occurs above criticality and at medium Reynolds numbers and simplifies in the infinite Reynolds number limit.
According to a study of Alemany and Zanette (see [81]) Tsallis' extended thermodynamics is directly related to Lévy walks and flights. And Lévy statistics had been applied to successfully describe turbulence (see e.g., [3,22])! A direct conclusion is that Tsallis' contributions to statistical physics yield a favorite tool to build a solid foundation of thermodynamics of turbulence. However, Gotoh and Kraichnan [90] offered critical remarks concerning the application of Tsallis thermodynamics with a single q-value to turbulence. They also discussed the application of a wave number dependent q value to overcome occurring problems. Such emerge naturally in our article in the context of modeling turbulent fields with intermittency (compare with Equation (130)).
Furthermore, it is known that Lévy random walks and flights can be related to fractional Langevin and Fokker-Planck equations that are based on fractional calculus. It is evident that the dynamics belonging to Richardson's and Mandelbrot's fractal geometry is fractional dynamics and calculus. Egolf and Hutter [26] introduced fractional calculus to turbulence modeling. Therefore, it is a direct consequence that Tsallis' thermodynamics can be formulated with the help of fractional calculus. Without this tool it is impossible to generate by a differentiation of a constant energy the Kolmogorov-Oboukov k −5/3 energy intensity spectrum. Moreover, a fractional generalization of the enstrophy is of importance to generalize Kraichnan's energy and enstrophy intensity spectra.
With all these new developments, it was tempting for us to apply Tsallis' thermodynamics to generalize Kraichnan's important work on the Ritz-Galerkin truncated turbulent energy and enstrophy spectra of 2-d turbulence. The new findings, related to the energy and enstrophy spectra, are consistent with three well-known types of descriptions of turbulence with an increasing level of complexity: (1) the BG thermodynamic results of Kraichnan; (2) the spectrum of Richardson, Oboukov and Kolmogorov; and (3) newer results of Lévy and Frisch, etc., describing turbulence including intermittency effects.
In the context of the proof of the generalized 2-d turbulence spectrum of Kraichnan, it has become clear that the usual definition of enstrophy with a square of the first derivative of the velocity relates intimately to Boltzmann-Gibbs (q = 1) equilibrium thermodynamics and requires a generalization with e.g., the nonextensive (q = 1) Tsallis thermodynamics including fractional derivatives. In future, a generalization of the definition of enstrophy in articles and standard text books may be demanded.
All these new studies show that non-equilibrium thermodynamics and fractional calculus, involving nonlocality, are likely to play a crucial role in modern turbulence research and computation. Whereas such theoretical approaches have started to actively bloom, time seems now also ready for fluid dynamic computational experts to introduce into their numerical algorithms and codes new fractional concepts. Modern zero-equation turbulence models, as e.g., the DQTM [23], are more than just a compensation of the eliminated wave number terms [96] e.g., by a Ritz-Galerkin truncation method. They contain scaling laws and couple, in a statistical sense, small with large wave number eddies in a nonlocal manner. Therefore, we conjecture and predict that in future low-order fractional turbulence modeling will be superior and more accurate compared to the past and present results of higher-order conventional calculus computations, based on linear and local concepts, and will decrease the demanded computation power (CPU time) likely by at least an order of magnitude. Energy Cascade where, different than in 3-d turbulence, in the 2-d case enstrophy is transported from large to small wave number eddies (up the turbulent energy cascade). In later years he was also involved in important developments on anomalous scaling and intermittency and in this context developed the exactly solvable Kraichnanian Turbulence Model. Intermittency occurs as a result of the presence of two phases in a turbulent flow, namely laminar streaks and swirling coherent structures and is directly linked with the phase change character of turbulence.

In memory of Robert Harry Kraichnan (1928-2008):
He received his PhD from MIT in 1949 and became a member of the Institute for Advanced Studies in Princeton. After working as theoretician at Columbia University and the Courant Institute of Mathematical Sciences (New York University) he became a free-lance researcher, supported e.g., by MIT, Los Alamos National Laboratory, NASA, etc. He received the 2003 Dirac Medal and numerous other prestigious prizes. Furthermore, he was the last assistant of Albert Einstein. Being a master also of (quantum) field theory and statistical physics, he attributed his excellent knowledge in a forty years' passionate activity to turbulence, which is one of the most important still unsolved problems of classical physics. A highlight of his activity in this field was the development of the Direct Interaction Approximation. It is a generalization of the mixing-length concept, introduced by Ludwig Prandtl. After 1967 he devoted strong efforts to two-dimensional turbulence, as it occurs, for example, in the atmosphere and oceans, and applied equilibrium thermodynamic methods, in an analogy to the theory of Boson gases, to fluid dynamics (see in this article). He discovered the Inverse Energy Cascade where, different than in 3-d turbulence, in the 2-d case enstrophy is transported from large to small wave number eddies (up the turbulent energy cascade). In later years he was also involved in important developments on anomalous scaling and intermittency and in this context developed the exactly solvable Kraichnanian Turbulence Model.

In memory of Robert Harry Kraichnan (1928-2008):
He received his PhD from MIT in 1949 and became a member of the Institute for Advanced Studies in Princeton. After working as theoretician at Columbia University and the Courant Institute of Mathematical Sciences (New York University) he became a freelance researcher, supported e.g., by MIT, Los Alamos National Laboratory, NASA, etc. He received the 2003 Dirac Medal and numerous other prestigious prizes. Furthermore, he was the last assistant of Albert Einstein. Being a master also of (quantum) field theory and statistical physics, he attributed his excellent knowledge in a forty years' passionate activity to turbulence, which is one of the most important still unsolved problems of classical physics. A highlight of his activity in this field was the development of the Direct Interaction Approximation. It is a generalization of the mixing-length concept, introduced by Ludwig Prandtl. After 1967 he devoted strong efforts to two-dimensional turbulence, as it occurs, for example, in the atmosphere and oceans, and applied equilibrium thermodynamic methods, in an analogy to the theory of Boson gases, to fluid dynamics (see in this article). He discovered the Inverse Energy Cascade where, different than in 3-d turbulence, in the 2-d case enstrophy is transported from large to small wave number eddies (up the turbulent energy cascade). In later years he was also involved in important developments on anomalous scaling and intermittency and in this context developed the exactly solvable Kraichnanian Turbulence Model.
Author Contributions: No experiments have been performed. Both authors have contributed essentially to the scientific content and writing of this article.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Proof of the Generalization of Kraichnan's Energy-Enstrophy Spectrum
To determine the energy-enstrophy spectrum of 2-d isotropic turbulence, we apply a Fourier transform to the generalized energy E ε ( x ): This is identical to: respectively:ˆ where the escort probability distribution function (105) was inserted. Then, this yields: It is clear that in the context of this proof also the enstrophy, Equation (42), requires a generalization by fractional calculus methods. Therefore, we introduce the modified enstrophy: with the fractional derivatives given by their two most common notations. Notice that, to perform this proof, we will not discuss which of the various fractional derivatives is best suited for this special application; this will be left for future clarification. The employed property of transformation of a fractional derivative of a certain order to a power law exponent with a value of exactly this order is valid for practically all versions of fractional derivatives, e.g., the Fourier fractional derivative [84], Riemann-Liouville fractional derivative [83], etc. It is straightforward to see that, with the help of Equations (41) and (48), this generalization leads to the following modification of Equation (51): where a possibly occurring constant in front of k 2q is absorbed in β. Now, an important question is, whether in the fractal case the subspace conserving energy and enstrophy still remains a quadratic form. It is expected that in phase space this could be a fractal subset (strange attractor). However, in Section 9 we prove that even if the total energy presents itself in Euclidian form that it is based upon fractal features with self-similarity. Therefore, in agreement with Kraichnan we write: u 2 = a ik x i x k = a 11 x 2 1 + a 12 x 1 x 2 + a 21 x 2 x 1 + a 22 x 2 2 , a 12 = a 21 .
Then a rotational main axis transformation ( x i → x i ) is applied and the coordinates are compressed or stretched in such a manner that relation (A7) simplifies to: where in the following for simplicity the convex curved overheads will be dropped. In Section 9 we have already outlined a motivation and proof of this relation. Notice that this does not imply a decoupling of Fourier modes. It is evident that the sum of the replicas for each spatial component inserted into the square quantities of Equation (A8) describes self and interaction mode coupling. Now, it follows that: This is inserted into (A4) and yields: 2 e −ik 1 x 1 e −ik 2 x 2 e −qβ q (α+β k 2q ) x 2 1 e −q(α+β k 2q ) x 2 2 dx 1 dx 2 . (A10) An advantage is that it is possible to separate the double-integrals in x 1 and x 2 :