Statistics of finite scale local Lyapunov exponents in fully developed homogeneous isotropic turbulence

The present work analyzes the statistics of finite scale local Lyapunov exponents of pairs of fluid particles trajectories in fully developed incompressible homogeneous isotropic turbulence. According to the hypothesis of fully developed chaos, this statistics is here analyzed assuming that the entropy associated to the fluid kinematic state is maximum. The distribution of the local Lyapunov exponents results to be an unsymmetrical uniform function in a proper interval of variation. From this PDF, we determine the relationship between average and maximum Lyapunov exponents, and the longitudinal velocity correlation function. This link, which in turn leads to the closure of von K\`arm\`an-Howarth and Corrsin equations, agrees with results of previous works, supporting the proposed PDF calculation, at least for the purposes of the energy cascade main effect estimation. Furthermore, through the property that the Lyapunov vectors tend to align the direction of the maximum growth rate of trajectories distance, we obtain the link between maximum and average Lyapunov exponents in line with the previous results. To validate the proposed theoretical results, we present different numerical simulations whose results justify the hypotheses of the present analysis.


Introduction
In fully developed turbulence, the finite scale Lyapunov exponents of the fluid kinematic field are of paramount importance because i) describe the turbulent energy cascade phenomenon, and ii) give the fluid viscous dissipation when the length scale goes to zero. One of the characteristics of these exponents in turbulence is that their statistics is related to the instantaneous velocity field, whereas does not depend directly on the time variations of this latter. Such characteristic, which represents the crucial point of this work, is consequence of the fact that the times of variations of the velocity field, which changes according to the Navier-Stokes equations, are much greater than those of the fluid displacements which follow the fluid kinematics [1,2]. Specifically, whereas the velocity field is a function of slow growth of the time, the distance between particles and the local fluid deformation are unbounded quantities which rise exponentially with the time [1]. Accordingly, the fluid particles displacements statistics, being not directly linked to the time variations of the velocity field, can be considered to be the result of the current fluid act of motion.
It is worth to remark that, the adoption of the finite scale Lyapunov exponents in place of the classical ones is justified by the fact that the turbulence is a complex phenomenon involving numerous length scales, each of them is characterized by different properties. The several perturbations of finite size will vary following nonlinear differential equations out of tangent space, thus the sole use of classical Lyapunov exponents is not adequate to describe the perturbations behavior associate to the different length scales and therefore the energy cascade.
Next, as a result of item i), the knowledge of the Lyapunov exponents statistics would lead to the closure of the equations of velocity and temperature correlations and therefore to the determination of kinetic energy and temperature spectra. Therefore, the distribution of such exponents is very important for quantifying the effects of the turbulent energy cascade.
Several works dealing with the closure of the correlations equations are present in the literature, as for instance [3,4,5,6,7,8,9]. Nevertheless, to the author knowledge, the influence of the finite scale Lyapunov exponents statistics on turbulent energy cascade and on the closure of the correlations equations has not received due attention. Hence, the objective of the present work is to develop a theoretical analysis based on the aforementioned properties which leads to determine the statistics of the local finite scale Lyapunov exponent in fully developed homogeneous isotropic turbulence for incompressible fluids. This local exponent, defined as provides the instantaneous growth rate of the distance ρ = √ ξ · ξ between two fluid particles trajectories x(t) and y(t) = x(t) + ξ(t), where ξ is the separation vector (finite scale Lyapunov vector). This exponent is linked to the so called finite size Lyapunov exponent, FSLE, defined in Ref. [10] in a more general framework of dynamic systems. Specifically, FSLE can be obtained in terms ofλ as a proper average ofλ over a given interval, say is an assigned threshold slightly greater than the unity to avoid interferences between the scales, whereas ξ(t k ), k=1, 2,... are rescaled along the direction ξ(t k )/|ξ(t k )|.
Because of non-smooth spatial variations of the velocity field,λ can exhibit fluctuations of sizable amplitude with respect to its average value, thus λ plays the role of a stochastic variable and will be distributed according to a certain PDF. To define the magnitude of these oscillations, average and maximum finite scale Lyapunov exponents,λ and λ + respectively, are so definedλ where • and • ξ ·ξ≥0 denote the average over the entire ensemble of ξ, and the average calculated on the ensemble whereξ · ξ ≥ 0.
As the consequence of the aforementioned properties, the present analysis assumes that the kinematics of a pair of fluid particles, characterized by ξ, is much faster and statistically independent with respect to the time variations of velocity field. This property, just discussed in [1,2] for what concerns the closure of von Kármán-Howarth and Corrsin equations [11] [12,13], was previously supported by the arguments presented in Ref. [14] (and references therein), where the author observes that: i) the velocity fields u(t, x) produce chaotic trajectories also for relatively simple mathematical structure of u(t, x). ii) the flows given by u(t, x) stretch and fold continuously and rapidly causing an effective mixing of the particles trajectories.
Through the hypotheses of fully developed chaos and fluid incompressibility, we first estimate the interval of variation ofλ, and thereafter determine the distribution function ofλ by maximizing the entropy associated to the fluid kinematic state. The maximization of such entropy is justified by the fact that the regime of fully developed turbulence corresponds to a situation of maximum chaos where the bifurcations cause a total loss of the initial condition data of x and ξ. As the consequence, we show thatλ results to be uniformely distributed in such interval, in particular, we determine the relationship betweenλ and λ + , resulting λ + = 2λ.
A further confirmation of such link is obtained by exploiting the alignment property of ξ, following which the Lyapunov vectors tend to align the direction of the maximum growth rate of ρ [15].
In order to compare the proposed statistics with the average energy cascade effects,λ and λ + are expressed in terms of the longitudinal velocity correlation function f = u r u ′ r /u 2 where u 2 = 1/3 u · u , u r = u(t, x) · r/r and u ′ r = u(t, x + r) · r/r. This relationship leads to the closure formulas of von Kármán-Howarth and Corrsin equations, and coincides with that just presented in Ref. [1] where the author adopts only λ + andλ without considering the distribution of the local exponents. There, the proposed closure formula leads to a value of the skewness of ∂u r /∂r equal to -3/7, in agreement with those obtained by the several authors with direct numerical simulation of the Navier-Stokes equations (DNS) [16,17,18], and Large-eddy simulations (LES) [19,20,21]. Therefore, the here proposed statistics should be adequate for describing the distribution ofλ, at least for what concerns the main properties of the energy cascade phenomenon. The novelty of this work with respect to the literature, and in particular with respect to Ref. [1], is represented by the statistics ofλ and its distribution function. This latter provides much more detailed statistical informations about dissipation and energy cascade than the analysis of Ref. [1] which, adopting average exponents, describes only the mean effects of energy cascade and dissipation. The present analysis provides, by means of such PDF, the estimation of all the Lyapunov exponent statistical moments. Accordingly, this study recovers, among the other things, the results of Ref. [1] and, in addition, gives all the statistical properties arising from this specific PDF, in particular, the deviations of such effects of dissipation and energy cascade with respect to their average values. Furthermore, to justify the hypotheses of the proposed statistics, the theoretical results of this latter are compared with numerical simulations of a proper differential system representing the incompressible fluid kinematics.

Interval of variation ofλ in incompressible turbulence
This section proposes an analysis for estimating the range of variations ofλ based on the hypotheses of fluid incompressibility and fully developed chaos.
The present analysis starts from the consideration that the turbulent energy cascade is related to the fluid particles trajectories divergence, therefore the relative fluid kinematics plays an important role in the estimation of the properties of such energy cascade [1,2]. The relative kinematics is expressed by finite scale Lyapunov vector ξ which satisfieṡ where u = u(t, x) varies according to the Navier-Stokes equations, and x(t) and y(t) = x(t)+ξ(t) are two fluid particles trajectories. The local divergence between x(t) and y(t) is quantified byλ which, due to the bifurcations of the kinematic field (kinematic bifurcations), exhibits oscillations of sizable amplitude with respect to its average value. These bifurcations continuously happen in those points of the space where det (∇u(t, x)) = 0 (4) Observe that the kinematic bifurcations defined by Eq. (4) are not the Navier-Stokes equations bifurcations (dynamic bifurcation), but arise from these latter [2]. In fact, the Navier-Stokes bifurcations frequentely occur determining continuously non-smooth spatial variations of u(t, x) which in turn lead to the condition (4) in the several points of the space. The definition ofλ given by Eq. (1) implies that ξ can be locally expressed as whereλ plays the role of the stochastic variable and Q(t) is a proper orthogonal matrix providing the orientation of ξ with respect to the inertial frame R. Accordingly,ξ isξ in which ω defines the angular velocity of ξ with respect to R.
For an assigned velocity field, finite scale Lyapunov exponents and vectors are formally calculated with Eq. (1) through the following orthogonalization procedure: 1) the maximal local Lyapunov exponent, sayλ 1 , is first obtained by choosing the direction e 1 (t) ≡ ξ 1 /|ξ 1 | which maximizesλ in Eq. (1), b) the second exponentλ 2 ≤λ 1 is calculated by selecting the direction e 2 (t) ≡ ξ 2 /|ξ 2 | in the subspace (plane) orthogonal to ξ 1 which maximizesλ, c) finally, the third oneλ 3 ≤λ 2 ≤λ 1 corresponds to the direction e 3 ≡ ξ 3 /|ξ 3 | normal to both e 1 and e 2 , where |ξ 1 |=|ξ 2 |=|ξ 3 |. The so obtained vectors system E ≡ (e 1 , e 2 , e 3 ) defines a rigid space which moves with respect to R with a given angular velocity ω E depending on the local fluid motion. Therefore, ξ k andλ k , k = 1, 2, 3, are locally expressed by The classical local Lyapunov exponentsΛ k are defined for ξ k → 0, k = 1, 2, 3. Now, in order to estimate the set of variations ofλ, observe that, due to fluid incompressibility, ∇ · u ≡ 0, and the classical local exponents obey to the following conditionΛ In general, Eq.(8) does not hold for finite scale Lyapunov vectors. Nevertheless, Eq. (8) is valid for those finite scale vectors for which the volume ξ 1 × ξ 2 · ξ 3 is locally preserved, therefore, without lack of generality, such these exponents can be written in the form where ε and λ m are variables depending on the current act of motion, and Now, following the hypothesis of fully developed chaos, it is reasonable that λ ranges in the set (λ 0 , λ S ), where λ S > 0 assumes its maximum value compatible with Eq. (9), and λ 0 is consequentely calculated. This implies It is worth to remark that λ S depends on the instantaneous velocity field which in turn is the result of time evolution of the fluid motion starting from the initial condition following the Navier-Stokes equations. Hence, at the current time, λ S assumes values related to viscous dissipation and kinetic energy both consequence of the fluid motion.

Incompressible fully developed turbulence
This section studies the distribution functions P = P (t, x, ξ) and P λ = P λ (t,λ) in incompressible fully developed turbulence, where P (t, x, ξ) and P λ (t,λ), are, respectively, the distribution functions of fluid particles position and Lyapunov vector, and of local Lyapunov exponent. In particular, we will show that the proposed statistics leads to the following relations and to the link between λ + and the longitudinal velocity correlation function. According to Eq. (1), P λ depends on P , and is expressed in terms of this latter through the Frobenius Perron equation where δ is the Dirac's delta. Hence, the distribution function P is first studied. This PDF changes with the time according to the Liouville theorem [22] associated to Eqs. (1). This theorem, arising from the following relation and from Eqs. (1), provides the evolution equation of P [22] ∂P ∂t where ∇ x · (•) and ∇ ξ · (•) denote the divergence of (•) defined in the spaces {x} and {ξ} respectively, and dx and dξ are the elemental volumes in the corresponding spaces.
Taking into account Eq. (15), and that the homogeneous isotropic turbulence is defined for unbounded fluid domains, P will satisfy the following boundary condition Accordingly, the statistical average of an integrable function of x and ξ, say ζ, is calculated in terms of P In particular, average and maximum finite scale Lyapunov exponents are Remark: It is very important to observe that the Liouville theorem in the form of Eq. (16) holds also when the variable velocity field is replaced with the same field frozen at the current time. This is due to the hypothesis that the times of variations of the velocity field, which changes according to the NavierStokes equations, are much larger than those associated with x and ξ whose times of variations are of the order of 1/λ. Following such hypothesis, the statistics of x and ξ does not depend on the time variations of u and can be considered to be the result of the instantaneous fluid act of motion.

Distribution function of x and ξ
For sake of our convenience, we introduce the quantity y ≡ (x, ξ) ∈ A ⊂ {y} which defines the relative kinematics, being {y} ≡ {x} × {ξ}, thus the Liouville theorem reads as Now, the entropy H associated to y is defined as Due to fluid incompressibility, ∇ y ·f ≡ 0, therefore from the Liouville theorem the entropy rate identically vanishes In order to estimate the steady distribution of y, observe that the fully developed turbulence corresponds to a situation of maximum chaos where the kinematic bifurcations cause a total loss of the initial condition data of y. Accordingly, it is reasonable that H assumes its maximum value compatible with Eqs. (15) and (21) with ∂P/∂t = 0. This corresponds to the following variational problem in which L(P, ∇ y P ) = −P ln P + ηP + χ∇ y · (P f) is the lagrangian of the problem, and η and χ = χ(y) are the Lagrange multipliers associated to the conditions (15) and (21), respectively. The maximum of J is then obtained as steady condition for J (δJ = 0) through the variational calculus and this leads to the Euler-Lagrange equation where A is the normalization constant related to η whose value is calculated with Eq. (15), whereas χ(y) is obtained through the Liouville equation. This latter gives in which (• ⊗ •) denotes the dyadic product between vectors. Integrating Eq. (29), we obtain where Φ is an arbitrary solenoidal vector field. Now, the determinant of (f ⊗ f) identically vanishes, thus Eq. (30) admits solutions only for Φ ≡ 0.
Moreover, as (f ⊗ f) exhibits minimum rank, Eq. (30) has only the trivial solution i.e. y is uniformely distributed on A ⊂ {y}, being in which m {•} indicates the measure of the set •.

Lyapunov exponents distribution function
In order to estimate P λ , the following should be noted: Ifλ is given, Eq. (2) corresponds to a hypersurface Σλ of the space {y} whose equation reads as Whenλ varies, Σλ and its points y * ≡ (x * , ξ * ) will change according to Eq.
where n ≡ ∇ y G/|∇ y G| and τ are, respectively, local normal and tangent unit vectors to Σλ both calculated in y * , being βτ not determined. On the other hand, dy * /dλ can be determined considering that x * does really not depend onλ, whereas the variations of ξ * are related to those ofλ by means of Eqs. (5) and (6). Thus, dy * /dλ is locally expressed as Asξ is solenoidal, the surface integrals of dy * /dλ · n over arbitrary hypersurfaces Σλ ( flow of dy * /dλ through Σλ ) assume the same value which does not depend onλ, i.e.
therefore, from Eq. (34), we obtain Now, the distribution function ofλ is calculated substituting Eq. (32) in the Frobenius-Perron equation (14) and considering that y is uniformely distributed on A P λ (λ) = y P (y) δ(G(y;λ))dy = 1 m(A) A δ(G(y;λ))dy where the surface integral over Σ λ ⊂ A, described by G(y;λ) = 0, is formally calculated according to the Minkowski measure theory [23]. Taking into account Eq. (37) and thatλ ∈ (−λ S /2, λ S ),λ results to be uniformely distributed in (−λ S /2, λ S ), being Hence,λ and λ + are calculated in terms of P λ with Eqs. (19) and (20) λ ≡ λ = λ S 4 > 0, Next, it is usefull to calculate the mean square λ 2 According to Eq. (41), the mean square ofλ equals the square of the average ofλ calculated forξ · ξ ≥ 0, and the standard deviation is proportional to λ It is worth to remark that the obtained distribution (39) is the result of three effects, of which two are in competition with each other. The first one of these, due to the Lyapunov vectors tendency to align the maximum growth rate direction of ρ, is responsible for variation intervals in whichλ >0, producing trajectories instability. The second one, related to the fluid incompressibility, acts in opposite sense preserving the volume, determining regions whereλ <0. The third element is the chaotic regime, here given by imposing H =max, which causes a continuous distribution ofλ ∈ (−λ S /2, λ S ).

Analysis through alignment property of ξ
One reasonable confirmation of the previous results is here given by exploiting the alignment property of ξ, following which ξ tends to align to the direction of the maximum growth rate of ρ [15], and the fluid incompressibility. This leads to an alternative way to achieve Eq. (12). Accordingly, λ + is now calculated adopting a proper PDF P + = P + (t, x, ξ) obtained as projection of P (t, x, ξ) at the time t+ τ , where x is considered to be constant and τ is the Lyapunov time defined by After the time τ , the alignment tendency of ξ provides that mostly all the Lyapunov vectors calculated at x =const will be such thatξ · ξ ≥ 0. The vectors lying in subspaces of {ξ} orthogonal to the maximum rising rate direction of ρ which do not follow such alignment form a null measure set in {y}, thus P + is a distribution function which satisfies and is expressed in function of P and τ where P = P (t, x, ξ) and ∇ ξ P = ∇ ξ P (t, x, ξ). Neglecting the higher order terms, P + is calculated as This PDF identically satisfies Eqs. (15). In fact, the integral over {x} × {ξ} of the first term at the R.H.S. of Eq. (46) is equal to one, whereas the integral of the second one can be reduced to a proper surface integral of P over ∂ {ξ} where P ≡ 0 through Green's identity, thus this identically vanishes. Furthermore P + exhibits the same entropy of P at least of higher order terms, which in turn does not vary with the time due to the fluid incompressibility The second addend at the R.H.S. of Eq. (47) identically vanishes as it can be reduced to be a surface integral of P over ∂ {ξ} where P ≡ 0. Therefore At least of higher order terms, P + maintains the same level of informations of P , thus it is adequate to estimate λ + . Accordingly, this latter is calculated as Substituting Eq. (46) and (43) into Eq. (49), we have Integrating by parts the second addend and taking into account the fluid incompressibility and the boundary conditions (17), we obtain in agreement with previous results.

Lyapunov exponents in terms of longitudinal velocity correlation
The results of the previous analysis allow to achieve the link between λ + and the longitudinal velocity correlation function. In fact, the standard deviation of longitudinal velocity difference directly depends on u 2 and f according to On the other hand, the Lyapunov theory gives the longitudinal velocity difference in terms ofλ Taking into account the isotropy and Eq. (41), (u ′ r − u r ) 2 reads as thus, λ + = λ + (r),λ =λ(r) are expressed in function of the finite scale r by means of the longitudinal velocity correlation f Equation (56) coincides with that proposed in Ref. [1] which leads to the closure formulas of the von Kármán and Corrsin equations. Unlike Ref.

Reference
Simulation R T H 3 (0) Present analysis ---3/7 = -0.428... [16] DNS 202 -0.44 [17] DNS 45 -0.47 [18] DNS 64 -0.40 [19] LES < 71 -0.40 [20] LES ∞ -0.40 [21] LES 720 -0.42 calculated with the proposed analysis and those obtained by the aforementioned works. It results that the maximum absolute difference between the proposed value and the other ones results to be less than 10 %. Furthermore, other studies [24,25,26] have shown that the closure formulas referable to this analysis provide kinetic energy and temperature spectra which exhibit scaling laws in agreement to the theoretical arguments of Kolmogorov, Obukov-Corrsin and Batchelor [27,28,29]. Therefore, the adopted hypothesis H=max and the consequent distribution (39) seem to be adequate assumptions at least for what concerns the estimation of turbulent energy cascade main effects.

Results and Discussions
In order to justify the plausibility of the previous hypotheses, this section presents one statistical analysis of numerical simulations of a simple differential system representing incompressible fluid kinematics. This system is properly chosen in such a way that it can exhibit simple mathematical structure, and chaotic behavior corresponding to an expected high value of entropy H. To achieve this latter condition, we adopt an adequate differential system which shows a "weak" or "reduced" link between velocity and spatial coordinates and an expected huge number of bifurcations per unit time (bifurcations rate). For this reason, each velocity component is assumed to be depending only on one single spatial coordinate with opportune scaling factors. Thus, the chosen differential system is given by the following where q, giving different scaling factors along the coordinate directions, will be properly selected to study the system behavior and to obtain a high bifurcations rate. The velocity field is periodic and C ≡ (0, 2π/q) × (0, 2π/q 2 ) × (0, 2π) represents the smallest regions of periodicity. The velocity gradient is then and its determinant, det(∇ x u) = cos qx cos q 2 y cos z, vanishes in those points where at least one of qx, q 2 y and z assumes values π(1/2 ± k), k = 0, 1, 2, .... Accordingly, large values of bifurcations rate are expected when q is opportunely high. Observe that, due to its peculiar analytical structure, the system generates trajectories which can differ from isotropic turbulence, while it is not sure that H assumes its maximim value. Nevertheless, as its Jacobian determinant vanishes in numerous points, it is reasonable that the proposed system shows a high entropy and a behavior which can be in some way compared to the fully developed chaos. Hence, although the regime H= max and Eq. (57) can correspond to different PDFs, say P λ and P λE , respectively, it is reasonable that the two distributions show elements in common, expecially for what concerns the interval whereλ ranges. In Table 2: Statistical parameters for |ξ| = 10 −7 . qλ where I 1 , I 2 ,λ/λ + and σ/λ are statistical parameters related to the PDF shape. Following Eqs. (59), it is plausible that the occurrencesλ > 0 are about two time thoseλ < 0, and that most of the occurrences ofλ happen in (−λ S /2, λ S ), where λ S is here obtained in function of λ  Tables 2 and 3 report the mentioned statistic parameters, and Figs. 1 and 2 show the distribution functions for different values of q. The set of differential equations is given by Eqs. (57) and by the corresponding evolution equations of ξ following Eqs. (3). During the process of integration, the Lyapunov vectors are continuously rescaled (at each time step) in order to maintain the initial scale |ξ(0)| at which ξ and λ are both referred. The scheme of integration is a fourth-order Runge-Kutta method with automatic adaptive step size, which returns equidistant samples x(t k ), ξ(t k ), t k = k∆t, k=1, 2, ..., N, being ∆t properly computed. Specifically, the time of integration is assumed T = 30000/λ 0max , beingλ 0max the maximum Lyapunov exponent at t = 0, and N= 2. 10 6 . Thereafter, the distribution functions ofλ are numerically computed through statistical elaboration of the simulations data. In order to obtain an expected high value of the bifurcations rate, the assumed values of q are such that q ∈ (2,25). For what concerns the initial conditions, all the simulations are computed for x 0 = 0.1, y 0 = 0.2, z 0 = 0.3, with an initial orientation of ξ corresponding to the minimum rising rate of ρ.
In line with the chaotic trajectories behaviors, we selected the simulations where x describes almost completely at least C, or at least the union of regions of space equivalent to C. The first set of simulations, obtained for |ξ|= 10 −7 , is reported in Fig. 1, and the other one, depicted in Fig. 2, is computed for |ξ| =10 −3 . Solid and dotted lines represent, respectively, P λE and the estimates of P λ obtained from Eqs. (39) and (60). It is apparent that, although the two PDFs differ with each other, the computed values ofλ are unsymmetrically distributed in such a way that I 1 , I 2 , standard deviation and exponents ratio, satisfy the relations (59). In particular, σ results to be proportional toλ in any case with a proportionality constant around to √ 3 (0.8 ÷ 2.5). We conclude this section by observing that variations of q and of initial conditions in terms of position x(0) and/or orientation of ξ(0), can produce changing in the shape of P λE and of the aforementioned parameters, depending on how the system describes its phase space. In particular, if x sweeps almost completely at least C -or at least the union of regions of {x} equivalent to C-λ is found to be unsymmetrically distributed, and the values of I 1 , I 2 , σ/λ andλ/λ + agree with those given in Eqs. (59). On the contrary, in the cases where x only partially sweeps C -or the equivalent union of parts of space-sizable differences can be observed with respect to the present results in terms of PDF shape and statistical parameters.

Conclusions
The distribution function of the finite scale local Lyapunov exponent of the kinematic field was studied in homogeneous isotropic turbulence. Based on reasonable assumptions regarding the fully developed chaos and the fluid incompressibility, the shape of such distribution and the range of variations ofλ are determined. This distribution results to be an uniform function in a proper non-symmetric interval of variations. The results arising from such PDF, in particular the link between λ + andλ and the closure for the longitudinal velocity correlation equation, agree with those presented in Ref. [1], and this should support the hypothesis H = max. An alternative way to determine the link between such exponents is also presented, which is based on the alignment property of the Lyapunov vectors. Direct simulations of a very simple differential system representing incompressible fluid kinematics give results which corroborate the hypotheses of the present analysis.
Thus, the von Kármán-Howarth equation gives the relationship between the statistical moments (∆u r ) 2 and (∆u r ) 3 in function of r.

Competing Interests
The author declares that there is no conflict of interests regarding the publication of this article.

Acknowledgments
This work was partially supported by the Italian Ministry for the Universities and Scientific and Technological Research (MIUR), and received no specific grant from any funding agency in the public, commercial or not-forprofit sectors.