Large Scale Emerging Properties from Non Hamiltonian Complex Systems

The concept of “large scale” depends obviously on the phenomenon we are interested in. For example, in the field of foundation of Thermodynamics from microscopic dynamics, the spatial and time large scales are order of fraction of millimetres and microseconds, respectively, or lesser, and are defined in relation to the spatial and time scales of the microscopic systems. In large scale oceanography or global climate dynamics problems the time scales of interest are order of thousands of kilometres, for space, and many years for time, and are compared to the local and daily/monthly times scales of atmosphere and ocean dynamics. In all the cases a Zwanzig projection approach is, at least in principle, an effective tool to obtain class of universal smooth “large scale” dynamics for few degrees of freedom of interest, starting from the complex dynamics of the whole (usually many degrees of freedom) system. The projection approach leads to a very complex calculus with differential operators, that is drastically simplified when the basic dynamics of the system of interest is Hamiltonian, as it happens in Foundation of Thermodynamics problems. However, in geophysical Fluid Dynamics, Biology, and in most of the physical problems the building block fundamental equations of motions have a non Hamiltonian structure. Thus, to continue to apply the useful projection approach also in these cases, we exploit the generalization of the Hamiltonian formalism given by the Lie algebra of dissipative differential operators. In this way, we are able to analytically deal with the series of the differential operators stemming from the projection approach applied to these general cases. Then we shall apply this formalism to obtain some relevant results concerning the statistical properties of the El Niño Southern Oscillation (ENSO).


Introduction
One of the main present challenges of modern science research is to give a formal and well grounded justification to the emergence of new general properties and laws of motion in complex systems when we change the temporal and/or spatial scales of observation, from the local one to the larger and collective one.How Thermodynamics, Onsager regression laws, Arrhenius law for chemical reactions can be directly derived from microscopic dynamics [1][2][3] are examples concerning this main issue, where for all of them the local time and space scales are order of nanoseconds and nanometers, respectively while the concept of "macrodynamics" applies to times larger just as some micro seconds and sizes from few micrometers.In many other contexts however, the local dynamics can take some days and involve tens or hundreds of kilometers in space, while the large scale emerging phenomena have time scales of some years and can have world wide effects, as happens, for example in the case of the El Niño/La Niña oscillations.
Actually, some examples where to apply the famous statement more is different are listed in the celebrated work of Anderson [4], where we can also find a first tentative of defining a hierarchy of different sciences disciplines, where the next one emerges as a larger scale view of the previous one.
A very interesting and particular discipline, surprisingly not included in the cited list of Anderson, is geophysical fluid dynamics, where different kind of dynamics and specific features continuously emerge, changing the space scale.Just a rough classification leads to start from the microscopic one, where friction and thermal dissipation dominate, to the mesoscales (order of tens-hundreds of kilometers) where eddy/vortex dynamics and locally driven currents are the major players, till to very large scale geostrophyic dynamics, where the inertia of large mass of water, the Coriolis force, the latitude dependence of the sun irradiation and the shape of the Ocean Basins play the main role.
It is clear that the emergence of the regular and universal large scale (let us say, in general, "macroscopic") properties from a complex/chaotic, irregular and system-dependent local ("microscopic") dynamics, implies a partition of the set of macroscopic systems into large equivalence classes of the corresponding microscopic ones.What are the criteria behind this classification is still a debated issue.
It is in other words to understand what are the dynamic and equilibrium properties, or the parameters of the system really relevant for the emergence of certain general macroscopic features.For example, it is well known that the chemical reaction rate (a macroscopic quantity) of a dilute solute in a solvent, depends only on some dynamical properties of the solvent (the fast "thermal bath"), as its correlation function and its susceptibility to external perturbations, besides to the reaction potential of the reactive particles [5].
Concerning the foundation of Thermodynamics, one way to address this issue is to take into account that Statistical mechanics and Thermodynamics are strictly related to the "standard" (i.e., with constant friction and diffusion coefficients) Fokker Planck Equation (FPE) for the Probability Density Function (PDF) of the system of interest.In fact, the thermodynamic properties of the system is included in the equilibrium solution of the standard FPE, i.e., the Canonical PDF (or Boltzman distribution) ∼ exp(−E/k B T), where E is the energy of the system of interest (a closed system that weakly exchange energy with its surrounding) and T is the temperature (k B is the Boltzman constant).Thus the universality of the validity of Thermodynamics and of the linear regression laws are led back to the ubiquitous emergence of the FPE and its Canonical equilibrium PDF [1][2][3].The transport coefficients of the FPE determine the macroscopic statistical dynamics and stationary features of the system of interest.We can say that the FPE is a direct viable root to Thermodynamics.Therefore the problem of finding the physical and formal ground to Thermodynamics starting from microscopic dynamics is shifted to the problem of obtaining the FPE for some large scale variables of interest hiding the details of the dynamics of the deterministic Hamiltonian microscopic system.
It is worthwhile stressing that, on the contrary, the FPE is usually obtained (and was originally derived) in the field of Markovian stochastic processes as the only possible truncation of the Kramers Moyal expansion of the Master Equation [6][7][8].In other words, in the stochastic picture, the FPE is equivalent to some (multivariate in general) Brownian Motion where the equilibrium PDF and the temperature are introduced "by hand" [9].
Here we shall show how it is possible to obtain a FPE for the PDF of some variables of interest, representing the large scale observation point of view of a general complex deterministic dynamical system, without introducing ad hock statistical assumptions.The variables of interest are the "observables" of the system, of which the measure is accessible only as smooth average quantities and superimposed fast fluctuations.Thus the FPE for the PDF of the variables of interest shall emerge just from the projection procedure.This means that the formal analytic expressions of the transport coefficients of the FPE, as functions of the dynamical properties of the microscopic system, shall give automatically the criteria for the already mentioned partition into large equivalence classes of deterministic systems for which similar large scale behaviors emerge.Actually, this is just what was already done in some previous papers concerning the foundation of Thermodynamics [1][2][3], by using a Zwanzig-like projection formalism.However, in these cases the "fundamental" approach requires to stay in a Hamiltonian framework.The Hamiltonian nature of the system simplify a lot the projection formal procedure, and in particular allows to skip the very difficult problem of managing the infinite series of differential operators that correspond to the diffusion parts of the FPE (see next Section).In fact for Hamiltonian systems a fundamental relationship exists between the diffusion (second order partial derivatives) and the dissipation (first order partial derivatives) terms of the FPE.This relationship, often cited as the Fluctuation Dissipation theorem, allow us to deal only with the deterministic parts of the FPE, i.e., the first order partial derivatives, and to obtain the diffusion coefficients just by using the Fluctuation Dissipation relationship (see Reference [10] for a clear example of such an approach).From the physicist point of view, the fluctuation and dissipation processes come out from the perturbation and the corresponding reaction forces, respectively, of the fast microscopic degrees of freedom (often called thermal bath) to the slow system of interest.Thus fluctuation and dissipation are the way the exchange of energy between the macroscopic and the microscopic variables is balanced at equilibrium, from a statistical point of view.
On the other hand in many other contexts of interest, as climate dynamics problems [11] or some geophysical fluid dynamics phenomena like El Niño/La Niña episodes [12], the system of interest of the corresponding models is intrinsically dissipative and, more in general, non Hamiltonian (for example, the interaction between the atmosphere and the Ocean is usually described by a non Hamiltonian formalism [12][13][14][15][16]).In non Hamiltonian cases fluctuation and dissipation could not balance to each other, they could be of different "strength" and origin.Moreover, the equilibrium PDF (if any) is not known "a priori" and often it is not possible to obtain an analytic expression for it.Therefore in these non Hamiltonian cases where we apply the projection approach to hide the variables that we do not (or we cannot) observe, we cannot avoid working with power series of differential operators.Here we propose the use of the Lie algebra formalism to manage these cumbersome formal expressions.
In some recent works [17,18] we started this approach considering two dimensional quasi-Hamiltonian systems of interest, where the broken of the Hamiltonian structure is only due to a standard friction force and to the non Hamiltonian nature of the interaction with the rest of the system.In these particular cases we were able to re-sum the power series of operators stemming from the projection approach.
In the present contribute we retrace most of the steps of the paper [18], but we shall release the limit in the number of degrees of freedom and we shall stay in a more general (also non Hamiltonian) framework, taking advantage of a more deeper use of the Lie algebra formalism.However we don't want to get lost the reader (and ourselves) with too general and formal issues, thus we shall introduce and use some results of the Lie algebra just for what we strictly need for our task, without pretending to be rigorous from a mathematical point of view.For a little bit more formal and rigorous approach of the Lie evolution of differential operators we refer the reader to another work just submitted elsewhere for review.

Defining the Class of Systems
According to what we have stated in the Introduction, our aim is to obtain the large time scale statistical behavior of some variables of interest, interacting with a given deterministic system of which the detailed dynamics is not known.This means that we assume that the whole system can be divided, from a dynamical and formal point of view, in two parts.The first one, called system of interest (or part a), represents the large scale observation point of view.Namely, the variables of interest or the "observables" of the system, of which the measure is accessible and where the path of the single trajectory looks quite random, but averaging the trajectory paths we have, for large times, a smooth dynamics.The second one, that we name booster (or part b) [5,[19][20][21], represents the "microscopic" dynamics of the system, of which the detailed trajectories are not accessible to the measure.We suppose that the interaction between these two parts is weak enough to justify a perturbation approach and we are interested to observe the system for large times respect to the typical time scale of the booster.Notice that in principle we do not make any assumption concerning the time scale separation between the dynamics of the unperturbed booster and that of the unperturbed system of interest, however, we must be aware of the fact that any resonance condition between the two systems could severely affect the perturbation approach.
From a formal point of view, the above picture can be made explicit in the following way (the dot stays for d/dt): where the N variables of interest are represented by the vector x (for example, the position and the velocity of a reactant particle [22] or the equatorial Pacific thermocline depth and the sea surface temperature anomalies for the El Niño-La Niña phenomena [23]), while (ξ, π) are the variables of the states of the booster.Note that the (usually collective) booster variable ξ forces the system of interest via the term ∆ I(x)ξ and, in turn, the booster dynamics is affected by the smooth function ∆ h(x), representing the "reaction force" of the system of interest.
For the sake of simplicity here the system of interest is directly perturbed by only one booster variable, and that is why we need to represent it with a specific symbol (ξ).Actually it is straightforward (but a little cumbersome) to generalize the present approach to the case where all the booster variables enter in the first line of the system of Equation ( 1), instead of only ξ.This would involve the cross correlation matrix and the vector of the response functions of the booster variables instead of the autocorrelation function (Equation (10)) and the response function (Equation ( 23)) of the sole ξ variable, respectively.
The ∆ parameter controls the interaction strength between the system of interest and the booster.For vanishing ∆ we say that the two systems are unperturbed.The interaction between the two parts is assumed to be weak, thus for large time scales the statistics of the new dynamics of the system of interest, emerging from the weak interaction with the booster, takes times order of 1/(∆ 2 ).We assume also that the dynamics of the unperturbed system of interest, guided by the velocity vector −C(x) is known.For example, −C(x) could be a smooth electro-magnetic field, or a linear forcing as in the case of the recharge oscillator model mimicking the El Niño Southern Oscillations (ENSO) [12][13][14][15], a pressure field, the gradient of a potential with a reaction threshold for chemical reactions etc.
For what the dynamics of the booster is concerned, we want to leave it as general as it is possible, thus we don't assume to know the explicit expressions of the vector velocity field (F, Q) but we shall make just some weak assumptions that we need to apply the projection approach.Before that, however, we must introduce the Liouville picture of the description of the system.By assumption, the macroscopic dynamics cannot depend too much on the specific state of the whole system (the "microscopic state"), thus, as already stressed in the Introduction, our effective macroscopic state is defined, at any given time, by a PDF on the whole states space.
In other words, we take advantage of the Gibbs idea introduced in the statistical mechanics context: Gibbs considered an infinite number of copies of the original system, where each copy differs from the other from the microscopic point of view, but not from the macroscopic one, he gave this construct the name ensemble, and studied the evolution of the ensemble density in the whole phase space.Thus, following Gibbs we study the behavior of large scale "observables" of interest of a general class of dynamical systems, where the only possible measurable quantities of the system are averages in the states space of the system, using a PDF as weight.According with that, we assume that we operate in a (at least Banach) space where a linear temporal evolution operator L of the PDF ρ(x, ξ, π; t) is defined, namely: We call L the Liouville operator or the Liouvillian.Thus we introduce also the associated Liouville adjoint operator L + in the following way: where the function f (x, ξ, π) belongs to the Banach space, and the symbol . . .means the average with the PDF ρ(x, ξ, π; t) as weight.
From a formal point of view, the Liouvillian of the system of Equation ( 1) can be rewritten as where L a is the Liouville operator of the unperturbed (i.e., for ∆ = 0) system of interest that, using the continuity condition, has the following explicit expression (repeated indexes means summation from 1 to N, where N is the number of variables of interest): while is the Liouvillian that perturbs the system of interest and, finally, M(ξ, π, ∆ h(x)) is the formal Liouville operator concerning the booster system.As for the related velocity field (F, Q), we leave completely unspecified this operator (it can be a differential operator, a master equation, etc.), but we assume that for small ∆ we can use the following expansion in power series of ∆: where the operators L b and W are the zero-th (unperturbed) and first-order contribution, respectively, of this power series expansion.In this way the Liouvillian of Equation ( 4) can be rewritten as For further use, we introduce the left and the right representations of a differential operator of the vector field with the following definition: in general a Liouvillian is written in the left representation, as for L a and L I in Equations ( 5) and ( 6) respectively, i.e., in the form ∂ j f j (x), with the components of the basis of the vector field on the left respect to the velocity field.This is because it evolves the PDF in the expectation calculus of the observable of interest.However in standard (pushforward) differential calculus a vector field is written in the right representation g j (x)∂ j i.e., with the components of the basis of the vector field on the right.Of course we have ∂ j f j (x) = f j (x)∂ j + (∂ j f j (x)), thus the left and the right representations coincides for solenoidal vector fields (as it is in the Hamiltonian cases).Notice that in the Liouville framework, to any Liouvillian L = ∂ j f j (x) (left representation) is associated the adjoint one defined as L + = − f j (x)∂ j (right representation), that, in the Lagrangianian point of view evolves directly the observable of interest inside the expectation calculus, instead of evolving the PDF.
Now, let us better specify the booster properties.From Equation (1), we see that the booster is defined by the following general equations of motion: where we have used K(t) as a generic time dependent external perturbation in place of the function ∆h(x).Thus, for K = 0 the booster is unperturbed.The main assumptions concerning the booster are: (i) the autocorrelation function ϕ(t) of the booster variable ξ decays in time so quickly that it has finite moments (finite correlation time), (ii) the booster responds linearly to the external perturbation K(t) (linear response).

The definition of the (normalized) autocorrelation time ϕ(t) is
where the symbol . . .0 means an average over the unperturbed equilibrium PDF.The last equality in Equation ( 10) follows from the definition of the Liouville operator of the unperturbed booster.Thus we have which is, at the same time, the definition of the booster correlation time τ.The (ii) hypothesis of linear response must be considered from a statistical point of view.Namely, let us assume that for t < 0 the booster is unperturbed (K = 0) and that the corresponding mean value of ξ vanishes: ξ 0 = 0 (it is always possible to redefine the booster variable to satisfy this condition).If a perturbation K(t) is applied at time t = 0, for linear response we mean that ξ(t) K , i.e., the average value of ξ, for t > 0, can be approximated by: where S(t) is the booster response function.In the special case where the external perturbation is constant, Equation ( 12) yields where the susceptibility χ(t) is defined by Thus χ(0) = 0 and, due to the property (i), we assume that it reaches in times of the order of the correlation time τ, the stationary asymptotic value χ(∞) ≡ χ.We introduce now another dynamical function, c(t), defined by with the property c(0) = 1 and c(∞) = 0.The function c(t) is related to the response function of Equation ( 12) by the equation The decay time of the susceptibility is defined as In the Liouville picture, the linear response hypothesis has been implicitly in part used with Equation (7), that, with the substitution ∆h(x) → K(t), implies; where ℘(ξ, π; t) is the reduced PDF for the sole booster variables.Moreover-although arguable from a rigorous point of view-we make the ansatz that from the linear response hypothesis it follows that it is possible to expand in power series of K the perturbed PDF of the booster around its equilibrium unperturbed state: with ℘ 1 (ξ, π; 0) = 0. Inserting Equation (19) into Equation ( 18), we get This first-order PDF can be used to write the response of the variable ξ of the booster as Interchanging the order of the temporal and phase space integrations (see below for a comment about that) we have Comparing Equation (22) with Equation ( 12) we see that the response function of the booster variable ξ to external perturbations is given by It must be pointed out that interchanging the order of integration between the states space and the time, as done to pass from Equation (21) to Equation (22), is not allowed in the case where Fubini's theorem does not hold, as for example when we have logarithmic divergences due to the presence of pure gravitational or Coulombian forces.Thus we exclude such situations or we assume that a cut-off function can be introduced to avoid such divergences without modifying the physics processes.Moreover, it is worth mentioning that while we have clearly stated that the booster is usually a deterministic system (it should be intended as the representation of the fundamental deterministic dynamics behind the observed large scale phenomenon), the present approach applies also to stochastic booster.Actually the correlation function and the response function can be often more easily obtained from stochastic systems rather than from deterministic ones.Moreover, we could also include Lèvy and more in general nonlocal diffusion processes.Namely, we could take into account the many cases that can be treated by the fractional Laplace formalism [24].However, for the sake of simplicity we do not consider here these particular, but not peculiar, cases.
We think worthwhile stressing again that in the Hamiltonian case, where the perturbation Liouvillian W is a Poisson Brackets, if we assume that the equilibrium PDF is the canonical one, i.e., ℘ eq,0 (ξ, π) ∝ exp(−H b (ξ, π)/k B T), where H(ξ, π) is the unperturbed booster Hamiltonian, from Equations ( 12) and ( 23) we get the Fluctuation Dissipation Theorem: c(t) = ϕ(t), as in standard linear response theories such as the celebrated one of Kubo [8], from which ϑ = τ.However in general the booster is not Hamiltonian and we don't know the exact expressions of the equations of motion, thus we do not assume that the equation c(t) = ϕ(t) holds true.
Note that the request that the booster responds linearly to external perturbation is quite strong.In fact the hypothesis of linear response, thought reasonable for many cases of interest, it is not so easy to demonstrate in general (see [1,2,25,26] for a deep discussion about the linear response of Hamiltonian and non Hamiltonian systems).
Actually, the result of Equation ( 23) is not new and is one of the steps of both [3,18] that we retrace in the present paper for the reader convenience, as declared at the end of the Introduction.
As already stressed, starting from the Liouville equation of Equation ( 8) and using a Zwanzig-like formal projection approach in the perturbation version of References [3,10,27], the equation of motion for the PDF of the part of interest can be obtained.This approach hides the dynamics and the initial state of the "microscopic" part (ξ, π) by an integration procedure, by using the assumption (i) for the booster, and considering observation times "t" much larger than the booster correlation decay times.Consequently, the evolution of the PDF of the part of interest is not anymore deterministic and, to the lowest non vanishing order on ∆, the time evolution of the reduced PDF results governed by the following integro-differential equation: where ϕ(u) and S(u) are just the booster normalized auto-correlation function and the booster response function as defined in Equations ( 10) and ( 23), respectively.Notice that, although to get rid of the non homogeneous term (initial state) of Equation ( 24) we used both the assumption (i) for the booster, and the fact that we are interested to large time scale dynamics (i.e., t >> τ, see [3] for details about this point) we are not allowed to substitute with ∞ the upper limit of the integrals.In fact we do not know if they are convergent or divergent for t >> τ.Actually, even if the autocorrelation function and the susceptibility vanishes for times much larger than τ, we still do not know what happen for the rest of the expression under the integrals, that are evolved by the unperturbed Luiovillian L a .
As it has been shown in [18], the last term in the r.h.s. of Equation ( 24), stemming from the reaction of the system of interest on the booster, is a first order partial differential operator, i.e., it induces an additional deterministic drift that, in certain cases, can be identified as a dissipation force.We shall call this operator the feedback term of the effective Liouvillian of the system of interest.In the present paper the main focus shall be on the second term in the r.h.s. of Equation (24).In particular we shall try to obtain an analytic expression for the operator that we find in this term and that is represented by the expression Expanding the exponentials in the above expression, we get, in general, an infinite power series of differential operators, containing all the derivative orders.This power series of differential operators is usually drastically simplified by using some approximations, involving time scales separation (see for example Reference [27]), or specific system arguments [10].In Reference [18], for the specific case we dealt with there (two dimensional Hamiltonian system of interest with a linear dissipative force), we have shown that Equation (25) gives rise exactly to a first order partial differential operator and we were able to obtain the analytic expression of this first order partial differential operator.
Here we want to generalize these results, i.e., we shall show how and when it is possible in general to get a formal expression of our central operator of Equation (25), in terms of the basis (∂/∂x 1 , ..., ∂/∂x N ) of the vector field of the fluxes on the states space of the system of interest.To simplify the notation, let us define ∂ j ≡ ∂/∂x j , j = 1, ..., N, and drop out the argument x in the unperturbed (∆ = 0) velocities vector field C i of Equation (5).
We think worthwhile to point out that the Lie algebra approach we shall introduce to work with the series of operators given in Equation ( 25), leads to general results that, as we shall see in the next Section, allows us to include in a single formal framework also the feedback term (the last term in the r.h.s. of Equation ( 24)) already successfully treated in [18].In fact, this term contains the series of operators exp(L a u)h(x) exp(−L a u) that can be considered as a special case of the one of Equation (25), where the operator L I is just a function.

The Lie Derivative
To manage the power series of operators in the second and third terms of the r.h.s of Equation (24) we shall introduce now a Lie derivative definition that can be considered, in the way we use it, as a generalization of the ordinary back derivative, in the sense that when it is applied to functions it coincides with minus the ordinary time derivative along the unperturbed trajectory, but it is also applicable to operators.
The Lie derivative, also called the adjoint-Lie operator, A × [. . .] associated with a given operator A, applied to a given operator (or function) B, is (here) defined by the standard commutation rule: where C is the generic set of functions/operators following the Lie derivative and its argument.More precisely, the definition in Equation ( 26) corresponds to the Lie derivative "along A" of the generic operator B usually introduced in differential geometry.For the Lie derivative the usual Leibniz rules of the standard derivative holds (very easy to demonstrate): and, likewise: With this "operative" definition we introduce also the exponential adjoint Lie operator as that using the Hadamard Lemma (see Appendix A) can be written in the following way from which it follows the trivial but very important property: Actually the Lie derivative here introduced is formally similar to the evolution operator in the Heisenberg picture of Quantum Mechanics.However we do not use the Hamiltonian to evolve functions and operators, but the Liovillian.The former is a second order differential operator, while the latter is a first order one, a very important difference.Although we are working with Classical Mechanics, we do not transform our Liouvllian in the ordinary Poisson Brackets for two reasons: the first one is that we are working with non Hamiltonian systems, the second one is that we want to apply our Lie derivative and adjoint-Lie time evolution also to differential operators, and not only to functions.
Exploiting the above definitions and the Hadamard Lemma, and making explicit the formal expressions of the Liouvillians L a and L I given in Equations ( 5) and ( 6), respectively, we can rewrite Equation ( 24) as (we use also Comparing the above equation with the previous equivalent Equation ( 24) the advantage should be clear: in the previous one we have some exponentials of the unperturbed Liovillian L a applied generically to all the expressions on the right of the same exponentials, including the PDF, without the possibility of using the useful algebraic distributive property e L a u BD • • • = (e L a u B)(e L a u D) . . ., because B and D are operators; in Equation (32), on the contrary, this property can be exploited.
It is noticeable, though not so surprising and already used in [18], that as it regards the Lie evolution of functions, the Liouvillian L a = ∂ i C i and the opposite of its adjoint C i ∂ i (from Equation (3) we have L + a = −C i ∂ i ) generate the same flux, corresponding to the simple backward evolution of the function itself (i.e., the left and right representations of the Liouvillian are, in this sense, equivalent, despite the fact that the velocity vector field is no solenoidal).In some way, this equivalence generalizes the case of symplectic or co-symplectic fluxes, where we have the identity and Here above we have used the following notation for the unperturbed backward evolution of x: Notice that in the r.h.s. of the above equation the time "t" does not enter, thus it should not be present also in the l.h.s. of the same equation.Actually, for the l.h.s. of Equation (35) we should use the more rigorous notation x 0 (x; u), where the first argument of the function x 0 (a; u) is the initial condition, i.e., x 0 (a; 0) = a.However, we think that the expression in Equation ( 35) makes the notation less heavy, and it can be directly connected to the usual Lagrange (or Langevin) picture, where x(t) = x.
The demonstration of Equation (33) follows by making explicit the commutators and using the Leibeniz rules for the derivative operator.Then, it is straightforward (though cumbersome) to check that the principle of mathematical induction can be used to shift this property to all the terms of the series that defines the exponential adjoint-Lie operator of Equation (34).
Equation (34) allows us to substitute in Equation (32) the Lie evolution of the functions I j (x), (∂ k I k (x)) and h(x) with the ordinary advective back time unperturbed evolution, i.e., I j (x 0 (t − u)), (∂ k I k (x)) x=x 0 (t−u) and h(x 0 (t − u)), respectively: As we have already stressed, the last line of Equation ( 36) corresponds to the deterministic drift (for example a friction) induced by the feedback of the system of interest on the booster.The second-to-last line is like the "noise-induced drift" that appears in Stratonovich systems because the correlation time of the "noise" (here the deterministic booster ξ variable) is not vanishing.Now we remain with the problem of evaluating the previous line (the second one), and in particular to solve the Lie evolution along the unperturbed Liouvillian L a of the basis of the vector field ∂ j .Unfortunately, unlike the functions, for the Lie evolution along L a of differential operators we have not, in general, the equivalence between L a = ∂ i C i and −L + a = C i ∂ i .However this equivalence is desirable because it makes possible a direct connection with some classical results on Lie algebra (that here we do not exploit), where the differential operators are usually written in the right form instead of the left one.As we show in Appendix B, if and only if one of the following equivalent conditions are satisfied: , the Liouvillian commutes with its adjoint, 2. ∀j | 1 ≤ j ≤ N, (∂ j ∂ i C i ) = 0, i.e., we have constant divergence of the velocity vector field C i , then for the Lie evolution of vector fields, L a is actually equivalent to −L + a , i.e., for any element α j (x)∂ j of the vector space, we have One example of class of systems that satisfy the above condition are the set of the co-symplectic Hamiltonian systems subjected also to linear non conservative forces (e.g., friction), such as Volterra gyrostats [28][29][30], spin systems interacting with magnetic fields, the inviscid or the dissipative Euler equation or some low order model approximations of the Navier-Stokes equations [28][29][30][31], including, for example, the famous Lorenz system [32][33][34].Another very important and large class of dynamical systems for which Equation (37) holds true is that of the standard symplectic Hamiltonian systems subjected to some linear friction or explosive force, as, just to quote a few examples among so many, the "ubiquitous" Duffing oscillator, the celebrated Recharge Oscillator mimicking the El Ninõ Southern Oscillation (ENSO) [12][13][14]23] or the usual picture of a chemical reaction process in a solvent, where a particle (the reactant) reacts escaping from a potential well, by jumping over a barrier of height E b [5,[35][36][37].Actually, for the pure symplectic/co-symplectic systems it is convenient to use a specific procedure to resum the series of operator given by the Lie evolution of vector fields, taking advantage of the energy conservation property of the flux (see [3,18,23] for details).
Here to stay in a more general context, we do not assume that the constant divergence condition holds true, however, the specific application examples that we shall consider later ahead, shall share the validity of this property.In Appendix C we show that the Lie evolution along a deterministic Liouvillian (i.e., that contains only first order partial derivatives as it is for L a ), of the basis ∂ j of the vector space is a first order partial derivative operator too, of the same left kind (1 ≤ j, k ≤ N): where β j,k (x; u) is a enough smooth function of x.The above results means that e L × a u [. . .] is a Lie Group on the vector field of the Manifold of the states of the system (or, better, on the tangent bundle of the Manifold).This is an almost trivial but very important result because it states that the second line of the master equation of Equation (36) contains second order (and not higher ones) derivatives, thus it gives rise to a diffusion process (plus other Stratonovich-like drift terms for the average evolution of the system of interest variables, of course).
From Equation (38) it is easy to obtain the β j,k (x; u) functions, in terms of the unperturbed trajectories of the system of interest (that we always suppose to know).In fact, first off all we note that by applying the l.h.s. of Equation (38) to x h (i.e., the h-th component of the vector x) and exploiting (34) we have: where is the unperturbed forward time evolution of x h .Then, by applying the r.h.s. of Equation ( 38) to x h we get: Finally, comparing Equation (39) with Equation ( 41), we get At the end we can summarize the main points of this Section in the following way: i.e., the Lie evolution of a product of operators/functions is equal to the product of the Lie evolution of the operators/functions; 2. e L × a u [ f (x)] = f (x 0 (t − u)), i.e., the Lie evolution along the unperturbed Liouvillian L a of an analytic function f (x) is not an operator, but just the advective backward time evolution of the same function; 3. e L × a u ∂ j = ∂ k β j,k (x; u), namely, the Lie evolution along the unperturbed Liouvillian L a of the first order differential operators is a Lie Group on the vector field of the Manifold of the states of the system of interest (or, better, on the tangent bundle of the Manifold).4. Once the unperturbed time evolution of the system of interest is known, from Equation (42) we have the coefficients of the Lie-evolved vector field.
With some algebra it is possible to directly verify that when the divergence of the velocity vector field C i is constant, i.e., when Equation (37) holds true, then we have ∂ k β j,k (x; u) = 0, thus both the left and the right representation of e L × a u [.

The Fokker Planck Equation
Now we are in the position to write the FPE for the PDF of the part of interest of our general system of Equation (1).Exploiting Equation (38) in Equation (36) we get (note: where β j,k (x; u) is given in Equation (42).
Notice that when both the divergence of the unperturbed velocity field is constant and the perturbation forcing is divergence free (as in the Hamiltonian cases), the third row in the r.h.s. of the above equation vanishes.
To make stationary the FPE coefficients, we should replace with ∞ the upper limit t of integration.This can be done taking advantage of the fact that we are interested in large time scale dynamics (compared with the decay time of the booster correlation functions).However, in general we cannot be sure that for t → ∞ the integrals in the above equations converge.For example, if I i ∝ x and the system of interest is a linear dissipative oscillator with friction coefficient λ, we have that x 0 (t − u) diverges as about exp(λt/2).This could make divergent also the integrals.In this respect, it is apparent that to avoid the divergence of the last integral of Equation ( 43), at least one of the following conditions must hold true: 2. the relaxation time of the unperturbed system of interest, that we indicate with 1/λ, is much larger than that of the unperturbed booster, i.e., 1/λ << τ (making the reasonable assumption that the susceptibility decays in a time order of τ).
Of course a similar analysis should be carried out for the second and the third integrals of the same equation.Taking into account that the convergence of the integrals must be checked on a case by case basis, in Equation ( 43) we replace with ∞ the upper limit t of integration and we get the following generalized FPE: in which the diffusion matrix D and the drift vector G, induced by the interaction with the external variable ξ, are given by and respectively, where, again, β j,k (x; u) is given in Equation (42).Equations ( 45) and ( 46) are important results because they give the formal transport coefficients of the generalized FPE in Equation (44), stemming from the perturbation projection approach applied to the very large class of dynamical systems given in Equation (1).At the same time, from Equations ( 45) and ( 46) it is possible to establish which are the conditions for standard statistics (Gaussian, Canonical, Lèvy, Thermodynamics, linear Onsager regressions [1][2][3], etc.) to emerge for the system of interest (see next Section for details about that).Thus Equations ( 45) and ( 46) can be used to identify the partition of the set of systems into equivalence classes for the emergence of different kind of regular and universal statistical properties.
Looking at these expressions for the transport coefficients, we see that they are expressed in terms of convolutions between the autocorrelation function, for the diffusion coefficients, and the response function, for the drift coefficients, of the variable ξ, with the unperturbed back-time evolution of functions of the variables of interest.Now the question is: we have got the FPE of Equation ( 44), but are the cumbersome expressions in Equations ( 45) and ( 46) really useful to obtain usable information about the transport coefficients of this FPE?
Usually the information about the booster dynamics is limited to experimental, observational or model data.The first case refers to laboratory experiments, thus the correlation function and the response function can be obtained as a fit of the data from targeted experiments.The second case concerns, for example, climatological or geophysical phenomena, as El Ninõ/La Ninã [23], that are large scale oceanic events induced mainly by the interaction with the fast atmosphere.Here the correlation function can still be obtained from data observation, but the response function of the fast atmosphere (for example) can be only obtained by simplified analytic or numerical models.
The last case refers to relatively simple low order models, usually derived starting from complex fundamental differential equations (for example fluid dynamical "building blocks" equations) of complex processes [38], reduced to a simplified set of ordinary differential equations for few variables (for instance the Lorenz models of atmospheric circulation [39][40][41][42][43]), through a series of approximations and hypotheses.In this case the autocorrelation function and the response function are obtained analytically or by numerical simulations.
It is clear that the expressions in Equations ( 45) and ( 46) can be managed only on a case by case basis.Here below we present how we can work with them in some representative cases.

Two Dimensional Hamiltonian and Quasi-Hamiltonian Systems of Interest: Recovering Previous Results and Thermodynamics
In [3,18] it has been shown that when the system of interest is both two dimensional and quasi-Hamiltonian, i.e., as already defined, the broken of the Hamiltonian structure is only due to one linear dissipation force and/or to the non Hamiltonian nature of the interaction with the booster, it is possible to directly re-sum the power series of operator exp(L × a u)[∂ i ] stemming from the projection approach.In this way a generalized FPE for the PDF of the system of interest is obtained.In particular, when the system of interest is purely Hamiltonian, Thermodynamics is recovered.Now we shall show that in these cases from the central Equation ( 42) we obtain the same results.
For a two dimensional oscillator, we have is the unperturbed potential (from now on, the prime symbol on a function means the partial derivative respect to the x variable of the same function), and Γv is a linear dissipative force.Moreover, for the sake of simplicity, we also assume that the interaction between the system of interest and the booster is Hamiltonian, i.e., we have an interaction potential is h(x) ξ, from which I x = 0 and and Equations ( 44)-( 46) become: and Using Equation ( 16), and integrating by part, the last equation can be written as: To recover the standard results concerning the derivation of the generalized Brownian motion from microscopic dynamics, we start searching for something looking like a friction and a potential term in the drift term of Equation (50).It is clear that if the system of interest is a linear oscillator then the unperturbed back-evolution of the velocity, v 0 (t − u), appearing in G v , is proportional to the "initial" values v and x: v 0 (t − u) = a(u)x + b(u)v, where a(u) and b(u) are functions with initial conditions a(0) = 0, b(0) = 1.Thus, if also the interaction potential h(x) is linear, both the friction and the renormalized potential are directly obtained from Equation (50).The same holds in the case of very large time scale separation between the system of interest and the booster dynamics.In fact, in this case we can use the approximation v 0 (t − u) ∼ v + V (x)u.In general, to get from the drift coefficient G v of Equation ( 50) a term proportional to the velocity v and one proportional to V (x) we use the identity v 0 (t and we obtain: In fact, in the above expression, we can associate the last term in the r.h.s., that is proportional to the velocity v, with a generalized friction stemming from the interaction with the booster, i.e., we introduce the following state-dependent friction coefficient: The previous equation shows that, in line with previous results (see for example [1][2][3]), the "reaction" of the booster to the interaction with the system of interest gives rise to a feedback term that is seen as an additive friction force for the same system of interest.According with the Linear Response Theory, the relationships given in Equations ( 16) and ( 23) allow us to relate the response of the booster (and then the friction coefficient) to some equilibrium correlation function of the booster.We stress again, however, that in general this correlation function is not just the autocorrelation function of the variable ξ perturbing the system of interest: it depends on the specific expressions of the booster equilibrium PDF and of the perturbation operator W of Equation (23).Thus the Fourier transform of the susceptibility is not (necessarily) directly related to the power spectrum of the forcing ξ(t).
Concerning the first and second terms in the r.h.s. of Equation ( 51), they can be considered as parts of a generalized force, that depends on both x and v, defined by Now, let us consider the diffusion coefficients of Equation ( 48), where we have to evaluate, for this case, the therms β v,x (x, v; x) and β v,v (x, v; x) of which the general expressions are given in the central result of Equation (42).In Appendix D we show that int the present case the following identities hold true: that, by the definition in Equation (38), means: The above result coincides with that in Equation ( 53) of [18].Using Equation (54) in Equation ( 48) we get and all the results of [18] follow directly.However, to show how Thermodynamics is recovered, we think worthwhile to delve into a little bit more detail about the pure Hamiltonian case, i.e., when we put Γ = 0 in Equations ( 47)- (56).First, we start observing that, with a little algebra, for Γ = 0, using Equations ( 51)-( 56) the generalized FPE of Equation ( 44) can be written in the following way: where the local "temperature" k B Θ is defined by the effective liouvillian for the system of interest is and, finally, the effective (renormalized) potential U(x, v) is defined by the following equation: The above results allow us to classify the systems from which we recover Thermodynamics as a large time scale emergent property.This happens when the equilibrium PDF of the FPE of Equation ( 57) is canonical: ρ eq,0 (x, v) ∝ exp(−H e f f /k B Θ) with H e f f = U(x) + v 2 /2, in fact it is well known that from the Canonical PDF Thermodynamics follows.On the other hand it is clear that the canonical PDF is the equilibrium solution of Equation ( 57) if and only if the temperature k B Θ is a constant (as it must be for a true temperature definition) and the renormalized potential U depends only on the x variable (as a true potential must be).It is easy to verify that this happens in the following cases 1. when, concerning the booster, the standard Fluctuation Dissipation Theorem holds true, i.e., the booster susceptibility decays as the auto-correlation function: c(u) = ϕ(u); 2. when both the system of interest and the interaction with the booster are linear: x 0 (t − u) = a(u)x + b(u)v, h(x) ∝ x, where a(u) and b(u) are functions of the time u; 3. for large (but finite) time scale separation between the booster and the system of interest dynamics: We refer to [3,18] for details, and we quickly focus our attention only to the first case.We think noticeable that in this case, in spite of the fact that the friction coefficient γ of Equation ( 52) and the diffusion coefficients D v,v and D v,x of Equation ( 56) depend on both the x and v coordinates, the temperature-like of Equation ( 58) is a constant and the renormalized potential is slightly different from the bare potential V(x): U(x) = V(x)∆ 2 χ h 2 (x)/2.Notice also that the condition c(u) = ϕ(u) can be realized in many different cases, and this fact contribute, in our opinion, to explain why Thermodynamics is so universal.For example it is not surprising that we have c(u) = ϕ(u) when the booster is a thermal bath, i.e., a thermalized system with an infinite number of degrees of freedom as in standard Linear Response Theory [8], or when we have a microcanonical booster with a large number of degrees of freedom [2].However, there are less obvious, but not less common cases where the standard Fluctuation Dissipation Theorem holds true, as for example happens for a general class of simple deterministic maps [20] or in geophysical fluid dynamics [25,44,45].

A Linear System of Interest with a Nonlinear Interaction with the Booster
Le us assume that the unperturbed velocity field of the system of interest can be approximated as being linear: C i = C ik x k , where 1 ≤ i, k ≤ N and C ik are constants, components of a N × N matrix that we name C. In this case the unperturbed time evolution of the system of interest is a linear function of the initial position: where B kl (u) are the components of the matrix B(u) given by from which Just for the sake of simplicity, we also assume that the reaction function h(x) is a linear function of the system of interest variables: h(x) = h k x k (it is straightforward to generalize the results releasing this assumption).Using Equations ( 61)- (63) in Equation (44), we obtain a more standard FPE: with From Equation ( 62), i.e., from the fact that the equation of motion of the system of interest are a system of linear ODE, it follows that in the above equations the functions B lk (u) can be expressed as linear combinations, with constant coefficients, of the real and imaginary parts of exp(λ i u), where λ i , i = 1, ..., N are the eigenvalues of the matrix C (see Equations ( 71) and (72) for the case N = 2).
The FPE of Equation ( 64) is written in the conservative form where the current j is given by j i ≡ (R i (x) + G ik x k + D ik ∂ k ) ρ(x; t).If the perturbation vector field I does not depend on x, it is easy to study the conditions for which the divergence of j goes to 0 for t → ∞.
If these conditions hold true, in particular, defining D and G as the N × N matrix with elements D ik and G ik , respectively, if the N × N matrix A defined by (superscript "T" denotes transposition) is positive-definite, then, we obtain the following Gaussian function as the equilibrium DF: as it can be directly verified substituting the above expression in Equation ( 64) and using Equation (69).
To show the importance of the result in Equations ( 64)-( 66), we consider a specific application: the Recharge Oscillator Model (ROM) [13,14,[46][47][48][49][50][51][52] mimicking the El Niño Southern Oscillations (ENSO): where the ζ = x 1 and T = x 2 variables refers to the depth of the thermocline of the West Equatorial Pacific Ocean, and to the Sea Surface Temperature (SST) of the East Equatorial Pacific Ocean, respectively.In this case the eigenvalues of the matrix C are λ ± = −Γ/2 ± iΩ where Γ ≡ γ ζ + γ T and Note that if Γ 2 /4 > ω 2 + γ ζ γ T , in the above equations we must substitute the trigonometric functions with the corresponding hyperbolic ones.Both the thermocline depth and the SST are subjected to a fast forcing by the interaction with the atmosphere, mainly due to the Madden Julian Oscillations (MJO) and the Westerly Wind Burst (WWB) [38,[53][54][55][56].These can be represented as multiplicative perturbations to the ROM [23,57,58].Considering the dynamics of the atmosphere as the generic "rest of the system" identified by the booster variable ξ, and the other "irrelevant" variables π as in Equation ( 1), we can write: where the interaction vector field (I ζ (T), I T (T)) depends on the anomalous temperature T as in the emulation of an enhanced air-sea coupling as SST anomalies increase (see for example [15]).
Comparing the above equation with that of Equation ( 1) we see that here the reaction term h(x) is vanishing.The goal is to describe the statistics of the ENSO, represented by the variables (ζ, T) perturbed by the atmosphere as in the above equation.Usually, for the sake of simplicity, the perturbation ∆ I ζ (ζ, T)ξ of the thermocline depth is not considered [23,57], but here we shall take into account also this effectively important contribution.From Equations (73) and (74), we have where Through the Zwanzig projection approach we arrive to the result of Equations ( 64)- (67) in which where From the above expressions we deduce that the diffusion terms of the FPE are second order polynomials in ζ, T with constant coefficients (the subscripts i, j stand for the subscripts ζ or T): Using then in Equations ( 78) and (79) the explicit results for the B ij coefficients given in Equation ( 72), we get that the terms D (k) ij are linear combination of the real and of the imaginary parts of the Laplace transform of the auto-correlation function ϕ(t) evaluated at the point, in the complex plane, corresponding to the eigenvalues λ ± = −Γ/2 ± iΩ of the unperturbed ROM.It is clear that in general, with second order polynomials as diffusion coefficients of the FPE, the stationary DF shall not be Gaussian, but (if it exists) typically it will be skewed with some kurtosis values and power law tails.Kurtosis and skewedness of the histogram of the frequency of the SST anomalies during El Niño/La Niña phenomena are actually obtained from the observation data [23].From the fact that the diffusion coefficients are second order polynomials of the system of interest variables, we can get closed first order ODE for the moments.In fact, it is straightforward to verify that from the FPE of Equations ( 76)-(79), we get a system of first order ODE for the n−th moments, that does not involves moments of higher order.For example, the equation of motion for the average of ζ and T is given by a couple of first order linear ODE, i.e., it corresponds to a forced dumped oscillator.This fact can be used to get easily all the relevant statistical information of the ENSO, but this is beyond the scope of the present paper.
There are of course many other cases of systems that can be treated in a similar way, for example connected RLC electric circuits, mechanical rotors supported by magnets [59], just to cite a couple of them.

Discussion
In this paper we show how and in which cases we obtain a large scale emergent dynamics for some variables of interest, weakly interacting with a general deterministic system.Because we assume that the single trajectories of the whole deterministic system are not accessible for observation, the information we can get are of statistical nature.Thus the knowledge of the PDF is the goal of the present approach.Actually the main result we obtain here is the generalized FPE for the PDF for the variable of interest.The diffusion-like coefficients of this generalized FPE are given in terms of the analytic expressions in Equation (42).The transport coefficients of the FPE determine the large time scale dynamics of the statistics of the system of interest, that takes times order of 1/∆ 2 , where ∆ is the small parameter that controls the weak coupling between the variables of interest and the general deterministic booster.
Thus, with this work we shed some light on the problem of the emergence of regular and smooth regression laws and standard Statistical Mechanics on Nature.Actually, the explicit formal results for the transport coefficients of this generalized FPE, allow us to establish which are the conditions for standard statistics (Gaussian, Canonical, Levy, Thermodynamics, linear Onsager regressions [1][2][3], etc.) to emerge for the system of interest.
Notice that, while usually the projection procedure is applied to Hamiltonian "microscopic" systems, where dissipation and diffusion are related by the Fluctuation Dissipation theorem, and stem from hiding the many and fast degrees of freedom of the "thermal bath" [10,[60][61][62][63][64][65][66], the results of the present paper allow us to generalize the procedure to the much broader class of non Hamiltonian systems.Thus the applications of our results are not limited to the field of foundation of Statistical Mechanics and Thermodynamics.Some examples are classic or quantum dissipative rotators, spin systems interacting with magnetic fields [67,68], as Landau-Lifshitz systems in micromagnetism [69], and many others from which an important class belongs to the field of geophysical fluid dynamics.In fact, while the Navier-Stokes (NS) equations in the simplest inviscid approximation become the Euler equation of motion, i.e., a pure conservative co-symplectic dynamical system, more in general, the NS equations contain dissipative terms, thus cannot have a pure Hamiltonian representation.Actually, they can be approximated by a set of coupled nonlinear ordinary differential equations, called low-order models (LOM), for example in the form of coupled Volterra gyrostats [28][29][30].
An interesting large scale geophysical phenomenon that can be approached in this way is the El Niño-Southern Oscillation (ENSO) [23], a naturally occurring event in the tropical Pacific that has global impacts of great relevance to society.The ENSO stems form the weak interaction between the dissipative Pacific Equatorial Ocean and the forcing atmosphere.The formal application of the results of the present work to the ENSO case has been introduced in Section 5.2 and, from a more quantitative point of view, this is an issue on which we are currently working.Appendix B Here we demonstrate the following proposition, introduced in Section 3: Proposition A1.If L a = ∂ i C i , the following conditions are equivalent: 1. L × a [L + a ] = (L + a ) × [L + a ] = 0, i.e., the Liouvillian commutes with its adjoint, 2. ∀j | 1 ≤ j ≤ N, (∂ j ∂ i C i ) = 0, i.e., we have constant divergence of the velocity vector field C i , 3. for the Lie evolution of vector fields of the tangent bundle, L a is equivalent to −L + a , i.e., for any element α j (x)∂ j of the vector space, we have e L × a u α j (x)∂ j = e (−L + a ) × u α j (x)∂ j , (A6) from which, taking the first order in u: Proof.We start the demonstration showing that condition (2) implies condition (1).This is very easy because , that is vanishing if condition (2) holds.The reversal is not so obvious, but it comes out considering L a as the generator of the flux in the Manifold of x: dt → the divergence of the velocity vector field C i is constant.Now we demonstrate that from condition (3) condition (1) follows.Actually this is easy, in fact it is enough to choose α j = C j in Equation (A7) to have directly L × a [L + a ] = (−L + 0 ) × [L + a ] = 0.For the reverse we use the equivalence between conditions (1) and (2) and we demonstrate that (2) implies (3).As in the cases of the previous propositions, to demonstrate Equation (A6) assuming condition (2) to hold, we start demonstrating the equation for the first order term, corresponding to Equation (A7), then it is easy, by induction, to shift the same property to all the terms of the power series that defines the adjoint-Lie evolution in Equation (A6).For this purpose we first use the Leibniz rule: L × a α j (x)∂ j = L × a α j (x) ∂ j + α j (x) L × a ∂ j .In the first term of the r.h.s. of this equation we already know that we can substitute L a with −L + a (see Section 3).Thus it remains to show that the same holds for the second term, but that is easy to do: 2) is satisfied, ends the demonstration.
Notice that when one of the three equivalent conditions of Proposition A1 are satisfied, using the Lie time evolution we can substitute L a = ∂ j C j with (L a − L + a )/2 and we recover the antisymmetric property of the Liouvillian operator also for non Hamiltonian (symplectic or co-symplectic) fluxes.This fact could be used for an eigenvalue/eigenvector approach to the Lie evolution of operators, but it is beyond the scope of the present work.

Appendix C
In this appendix we demonstrate Equation (38), i.e., that the Lie evolution along the unperturbed Liouvillian L a of the basis of the vector field is a first order differential operator.Actually we shall demonstrate something more general: Proposition A2.If A and B are generic partial differential operators with derivatives up to order m and l, respectively (including the case m = 0 or/and l = 0 when A or/and B are just functions of the variables of the system), then the Lie evolution of B along A, defined as e A × u [B], is generally a partial differential operator with