The Path Integral Formulation of Climate Dynamics

The chaotic nature of the atmospheric dynamics has stimulated the applications of methods and ideas derived from statistical dynamics. For instance, ensemble systems are used to make weather predictions recently extensive, which are designed to sample the phase space around the initial condition. Such an approach has been shown to improve substantially the usefulness of the forecasts since it allows forecasters to issue probabilistic forecasts. These works have modified the dominant paradigm of the interpretation of the evolution of atmospheric flows (and oceanic motions to some extent) attributing more importance to the probability distribution of the variables of interest rather than to a single representation. The ensemble experiments can be considered as crude attempts to estimate the evolution of the probability distribution of the climate variables, which turn out to be the only physical quantity relevant to practice. However, little work has been done on a direct modeling of the probability evolution itself. In this paper it is shown that it is possible to write the evolution of the probability distribution as a functional integral of the same kind introduced by Feynman in quantum mechanics, using some of the methods and results developed in statistical physics. The approach allows obtaining a formal solution to the Fokker-Planck equation corresponding to the Langevin-like equation of motion with noise. The method is very general and provides a framework generalizable to red noise, as well as to delaying differential equations, and even field equations, i.e., partial differential equations with noise, for example, general circulation models with noise. These concepts will be applied to an example taken from a simple ENSO model.


Introduction
The equations that govern the evolution of the atmosphere and the ocean have been known for a long time and have been extensively investigated. To investigate them, several numerical methods that exploit the first order time derivatives to obtain the time evolution, have been intensely developed. The equations showed a strong sensitivity to small perturbations, both in the initial conditions as well as in the parameters defining them, giving rise to the entire field of dynamical chaos [1].
The chaotic nature of the dynamics stimulated the application of methods and ideas derived from statistics and statistical dynamics. For instance, ensemble systems are used to make weather predictions which are designed to sample the phase space around the initial condition. Such an approach has been shown to substantially improve the usefulness of the forecasts since it allows forecasters to issue probabilistic forecasts. The implicit assumption is that the presence of various sources of errors, coupled with the intrinsic sensitivity of the evolution equations to small errors [1], makes a single forecast not so useful [2,3].
The concept has gained a large consensus because it has been shown to be relevant to various dynamical problems. Numerical experiments driven by external forcing, such as those used with prescribed SST (Sea Surface Temperature) or even prescribed concentrations of greenhouses gases in climate change experiments, have shown that the response to external forcing is still sensitive to errors, either because of uncertainties in the initial condition or in the model formulation. Ensemble experiments are now commonly used in these cases [4][5][6].
These works shifted the dominant paradigm of interpreting the evolution of atmospheric flows (and the ocean to some extent, see [7]) attributing an increasing importance to the probability distribution of the variables of interest rather than to a single representation. The ensemble experiments can be considered as crude attempts to estimate the evolution of the probability distribution of the climate variables, which is the relevant quantity for practice. Other interesting quantities, as variance and correlation functions, can be obtained from the Probability Distribution Function (PDF). The ensemble mean of temperature, for instance, cannot be considered simply as the average of the available ensemble members, but as the simplest estimation of the expectation value.
Finding an equation for the evolution of the PDF is far from trivial. Hasselmann [8] has shown that a stochastic component is consistent with the basic principles of the atmospheric/ocean dynamics and whereas other investigators [9][10][11][12][13][14] have shown that some aspects of the atmosphere dynamics can be described by simple models with a stochastic component. It is also possible to estimate the stochastic component from observations [15,16].
The addition of stochastic noise to the evolution equation results in a multidimensional Langevin-like equation that can be shown to support a Fokker-Planck equation for the evolution of the probability distribution of the state vector. This result is very interesting since the Fokker-Planck equation is linear, even if the corresponding evolution equation may be non-linear. However, the Fokker-Planck equation is obtained in a phase space with the dimensions corresponding to the number of degrees of freedom of the original equations. Even a very simple general circulation model can easily have hundreds of degrees of freedom and a numerical approach is not feasible. This paper shows that it is possible to write the evolution of the probability distribution as a functional integral of the same kind introduced by Feynman [17] in quantum mechanics, using some of the methods and results developed in statistical physics [18,19]. The approach allows obtaining a formal solution to the Fokker-Planck equation corresponding to the Langevin-like equation of motion with noise. The method is very general and it provides a framework easily generalizable to red noise, as well as to delay differential equations, and even field equations, i.e. partial differential equations with noise. The approach has been proved useful in fields other than physics, such as polymer theory, chemistry and even financial markets [20][21][22]. There are also applications to other relevant problems in geosciences : turbulence fluids [23][24][25], Lyapunov exponents [26], data assimilation [27], or wave propagation in random media [28,29]. The first quantum field theory formalism describing additive noise was developed by Martin, Siggia and Rose [30], by using a different kind of approach, a method similar to the canonical quantization. The path integral technique, however, is relatively less known in the field of Climatology.
In this paper, the authors attempt to solve stochastic differential equations with the Path Integral technique. This method is applied to solve a linear simple model and a non-linear one, relevant to climatological problems, to demonstrate the power of this tool. Although the technique seems involuted, it could be very easily generalized and could also be the basis for applications to field equations arising in a field theory. This method has only been used with simple linear and non-linear ENSO models, which contain only time depending variables. The aim of this paper is to stimulate interest in the path integral technique for application in the investigation of the Global Climate System. The authors' hope is to use the formalism of the field variables to face, with this technique, more complicated models, by applying this method to study general circulation models with noise.
The remainder of this paper is organized as follows. Section 0 introduces and summarizes the general theoretical foundation and Section discusses the calculation of the integrals. Section introduces the concept of Green's matrix and functions. Section introduces a discussion of perturbation expansion applied to nonlinear cases. In Section, these concepts are applied to an example taken from a simple ENSO model and Section concludes.

The Path Integral Formulation
Langevin equation and probability. The systems describing the atmosphere or the ocean can be written as coupled Langevin equations: where q(t)~(q 1 (t), . . . ,q d (t)) represents a trajectory in R d and f m (q) represents a differentiable function of q. It is assumed that there are d degrees of freedom, and in what follows it will be considered a Gaussian white noise E m (t). This kind of noise is characterized by its 1-point correlation functions, the averages, that are equal to zero, and by the 2-points correlation functions: In the equation above, d mn is the Kronecker delta and Q measures the strength of the correlation. For simplicity, Q is taken as a constant, and the variances of different E m (t) noise terms are equal. The equations above are not the most general stochastic first order differential equations. Time translation invariance has been explicitly assumed, and the same variance has been used for different variables, but those restrictions are not really limiting and it has been assumed for simplicity [31].
The Langevin equations (1) generate a time-dependent probability density function for a stochastic vector q(t), given the value of this vector at initial time, which can be written formally as: in which q 0 and t 0 are the initial conditions, and d is the Dirac delta. This probability is just the ensemble average over the solutions of the Langevin equations (1); STE denotes an average with respect to the probability distribution of the realizations of the stochastic variables E m (t). P(q,Tjq 0 ,t 0 ) is the conditional probability to find the system in q at time T starting from the point q 0 at the time t 0 . (q m (T){q m ) is the difference between a point of the trajectory obtained with the Langevin equation (1) at the time T, and a fixed point in the configurations space. The trajectory depends on the initial condition q 0 , at time t 0 . Although q 0 doesn't appear in the right-hand side of the equation, that expression implicitly depends on it by means q m (T).
Using the Gaussian nature of the noise, starting from the equation above, it is possible to write a Fokker-Planck equation for P(q,Tjq 0 ,t 0 ), see for example [31]: The formal solution of this equation can be written as a path integral [19] P(q,Tjq 0 ,t 0 )~ð where ½Dq(t) means that the integration is done over all paths q(t) that go from t 0 to T. The functional S(q) is the continuous Onsager-Machlup action which in the white noise case for the last equation summed over repeated index it is used. The extra divergence term in the action is associated with the difficulty of defining the derivative of a stochastic process. These expressions are symbolic and, they have to be defined by a discretization rule. In fact, a functional integral is well-defined only if it is assigned a formal continuos expression and a discretization rule. The process paths, which are solutions of the Langevin equation, are continuous as Dt?0, but they are not differentiable, and the ordinary rules of calculus must be modified to come up with a consistent definition. In the case of a simple additive noise, the pathologies do not show up, but if there is multiplicative noise, it is absolutely necessary to choose an interpretation. In the following, the Stratonovich interpretation will be used as the discretization rule, which allows treating the fields as differentiable, and therefore to use them in the ordinary rules of calculus. In the case of weak additive noise, the divergence term drops simplifying the action: Expectation values for a generic quantity F (q(T)) can be obtained by where S T here is just a time average, and P(q 0 ,t 0 ) is the distribution that describes the system at the initial time t 0 .
Integrating over the initial conditions using P(q 0 ,t 0 ), the average depends only on the point q(T). The correlation can be obtained by using a polynomial expressions of the q components on the functional F . Stochastic equations and path integrals have a mathematical meaning only if it is a discretization is associated to them. One can apply a discretization, for instance denoting the initial and final times by t 0 and T, respectively, with n~1, . . . ,N. The probability distribution of the discretized noise is given by : If the Langevin equation (1) is integrated in an infinitesimal time interval Dt, the discretized equation becomes The conditional probability, that the system will be in the state q nz1 at time t nz1 given that it was in q n at time t n , could be defined with the following symbol, where d is the Dirac delta. On the right-hand side of the equation above, the only variables which appear are q nz1 , q n , Dt, therefore it is necessary to always use a notation for which the transition probability depends explicitly on time, for instance t nz1 and t n . This will be more explicit, as it will soon be shown, when the summation in the action is transformed into an integral with extremes depending on the initial and final time of the transition, see Eq. (6). This means that time is a variable, not an index, and it is coherent with the fact that the PDF, which satisfies the Fokker-Planck equation, is time-dependent. In order to obtain the unconditional probability p(q nz2 ,t nz2 ), one would have to use the Kolmogorov-Chapman equation, the probability for the entire path can be obtained where it has been defined The S N functional plays the role of the action as in classical mechanics and it is also known as the Onsager-Machlup functional. Probability cannot be exactly analytically computed for a non-linear f, but with a linear f the integral is Gaussian and can be computed. From the Eq. (11) it is possible to see that (2pQDt) N=2 , these quantities always have to be considered in the limit approximation. There are N{1 integrations over the possible intermediate values of the path, and the end points q 0 , q N are fixed. Note that there are N factors in the denominator of the Eq. (11), 1 (2pQDt) N=2 , and so presumably a normalization factor will have to be introduced later, since they can be divergent when N??. The choice of the discretization is important because the term is ill-defined when the small time step limit is studied, and it must be treated carefully. It turns out that Feynman's original choice of symmetrizing the term [17] as is equivalent to choosing the Stratonovich interpretation. Different continuous formal expressions exist for the functional integral, which, with the appropriate discretization rule, define the same stochastic process. The Stratonovich mean point formulation is useful to analytically treat problems. In particular it is connected to the possibility of using the usual techniques of integral calculus. With this kind of discretization, it is possible to define all the terms that in the limits become the action seen before Eq. (6). The discretization, beyond giving meaning to the expressions above, gives a recipe to explicitly compute those quantities. The propagator. The probability of reaching q N at T from any point q 0 at t 0 , obeying the initial distribution P(q 0 ,t 0 ), is then given by which describes the evolution of the probability distribution from time t 0 to time T. It is the solution to the Fokker-Planck equation.
The final integration over q 0 resolves the normalization issues previously mentioned and a final result is obtained. It is also possible to write Eq. (13) as where a symbol for the kernel has been introduced that propagates the solution from time t 0 to time t N~T ; this expression is also known as the propagator. This equation is the analogous of Eq. (5) discretized. The concept of the path integrals recurring in these formulas is illustrated in Fig. 1. The probability of reaching q N starting at q 0 is composed by the sum of all paths that may take all possible intermediate values at intermediate times. Their contribution must be integrated for all possible values. For further details about the path integral, refer to [22].
Considering that: the expression for the probability in the continuous case is given by and continuous time propagator from time t 0 to T is Eq. (16) is the probability of finding the system in the state q at time T, given the initial distribution P(q 0 ,t 0 ) at time t 0 .

Calculating the Path Integral
Practically, analytically computable path integrals are rare, and they are essentially limited to Gaussian integrals, which, as previously noticed, are obtained when f(q) is linear. They can be analytically calculated from the discretization previously introduced only in particular cases [22]. It is possible to consider an approximate method for the computation derived from the steepest descent method (or saddle point method) [31]. The path integral is dominated by the minima of the action, which are the trajectories that minimize the action functional. It can be approximated by a series of Gaussian integral, one for each minimum of the action, considering fluctuations around these The Path Integral Formulation of Climate Dynamics PLOS ONE | www.plosone.org trajectories and computing the approximate integral. In this way the path integral can be separated into two simpler factors; the first one contains the stationary conditions, and the second contains a term that can be transformed, using the projection on eigenfunctions, in Gaussian integrals. If f is considered as non-linear, this method results useful because it lets us use a simple perturbation expansion technique.
Let the function f(q n ) be a linear operator A. In this case the action can be written as and the path integrals become Following the steepest descent method, trajectories that minimize the action must be found. However, there is a problem associated with the fact that, for a system of the present form, there are two solutions to the first order of variations, which correspond to the equation of motion without noise. The solutions correspond to the choice q~r, so that _ r r~+Ar r(0)~q 0 , the unperturbed trajectory corresponds to the plus sign. Obviously it would be desirable to be able to investigate the perturbation around this solution, but this is complicated because the particular value of the action in this case is zero, making a traditional expansion impossible. However, as pointed out by [32], there is a method that allows the expansion along the correct solution and also satisfies both boundary conditions for the integration in the action. It is necessary to introduce a change of variables quantity q~rzg, so that the action (18) can be written as because r satisfies the equation of motion without noise. The boundary conditions on this expression are given by The measure of the integral does not change, since it is a linear transformation, and ½Dq(t) is transformed in ½Dg(t) without the adjoint of a new factor in front of it. One can now substitute around an unperturbed trajectory g c (t) so that deviations of order ffiffiffiffi Q p are introduced obeying the boundary conditions y(0)~y(T)~0 Substituting Eq. (21) in the action (20), it is and integrating by parts the various terms and using the boundary conditions, it is obtained Therefore, if a g c is chosen, which satisfies the equation with the given boundary conditions the action can be divided into two parts: the explicit terms depending on the boundary conditions and implicitly on the unperturbed solution r, and a term that depends only on the fluctuations y, The term S 1 does not depend on the varying path y(t) and therefore can be taken out from the integration in Eq. (19), whereas the term S 2 will depend only on time T, which is often called the prefactor. The propagator (19) can then be written as with boundary conditions y(0)~y(T)~0. The remaining calculation can be finished by observing that the action in the paths y is then equivalent to a Sturm-Liouville boundary problem for the differential operator L The operator L is self-adjoint and therefore has a complete orthonormal set of eigenfunctions w l i with real eigenvalues m l i , l~1, . . . ,?, i~1, . . . ,d. The eigenfunction and eigenvalues are d-multiple infinities as a consequence of the dimensionality d of the operator. The variables y can be expanded in a series of the complete orthonormal eigenfunctions and Using this approach, the path integral Eq. (17) can be written as an infinite set of Gaussian integrals over the coefficients of the expansion. A change of variables from the q i 's to the c i 's will allow the execution of the integral. The functional path integral becomes an integral for the coefficients c 2 li , because in varying them, all possible paths are obtained. Since L is self-adjoint, it can be diagonalized by a unitary transformation with a unit Jacobian for the change of variables, therefore the path integral measure remains the same, and the boundary conditions are satisfied by the eigenfunctions. The integral is then formed by an infinite number of Gaussian integrals, and it can be obtained that, or K(q,T; q 0 ,t 0 ) lim L??
The product is reduced to the inverse root of the determinant of L.
This determinant and the constant, which contain the temporal step, are usually regularized considering the ratio between this propagator and the propagator for a free evolution.

Generating Functions
The calculation of the n-points correlation functions, that will be used to compute the correlations in the following examples, is complicated, but it can be simplified by introducing the moment generating functional The functional derivative of the expression above provides the expectation value for the mean. Remember that the functional derivative is defined as follows where d on the right-hand side is a Dirac delta, while the notation on the left is the usual notation for the functional derivation. The higher order correlations can be obtained by repeating the process: and for a generic functional F , it is possible to prove that The formalism of the derivation operator, appearing within the scope of F , means that one has to substitute the functional derivatives in place of the usual variables on which the operator F is defined as in the Eq. (32), where q m (t)q n (t) are substituted with The following paragraph shows how it is possible to compute Z½J for a general case with non-linear f in the Langevin equations.

Perturbation Expansions
Feynman diagrams. The path integral formulation adapts itself very naturally to the definition of perturbation corrections of various kinds, for example, it can be used to compute corrections to the probability distribution and to the correlation functions. In fact, because of the general complexity of the action, it will be difficult to know the exact distribution computing the integrals. Although the technique seems involuted, it can be very easily generalized and can be the basis for applications to field equations arising in a field theory.
Consider the propagator for a non-linear evolution, _ q q{Aq{mf(q)~0, where m is a parameter that measures the strength of the non-linear terms, that is an extension of (19). The same coordinate transformation described in Sect. (), q~rzg, can be introduced so that the action can be written as an extension of Eq. (20) Clearly the new measure ½Dg(t) also has to be considered. The quadratic nature of the action creates a potential problem because the expansion of the terms, according to powers of the coupling constant m, generate terms of the form _ g gf(g) that couples state variables with derivatives. It is possible to overcome this problem using the Hubbard-Stratonovich transformation [33,34] extended to the multidimensional case, that is a generalization of the identity for the functional integrals. If the propagator is considered in its discretized form Eq.(15), and for each integral that appears the identity above is used when the continuos limit is restored, the propagator becomes K(g,T; g 0 ,0)~ð D½y(t) The auxiliary functions y(t) are defined over the entire time axis. This transformation introduces new integrations that can be summarized as D½y(t). The field f(t)~{iy(t) can be introduced and the trace of the linear part can be taken from the functional integrals, as it does not depend on the paths, yielding or The subscript V has been added to underscore the dependence of this propagator on the non-linear terms in the second exponential exp (V ), whereas the quadratic terms are contained in S 0 . The term V (t) contains higher order terms in q(t) (hence in g(t) and w(t)) that reflect the impact of the non-linear interactions. The propagator corresponding to the quadratic part describes the evolution of the system without interaction and therefore can be described as the free evolution of the system. Usually it can be computed exactly: K 0 (g,T; g 0 ,0)~ð D½q whereas, in the presence of interactions, it is In other words, the propagator for the problem is the expected value of the interaction with respect to the probability distribution of the unperturbed, usually linear, problem. In the presence of a small coupling constant m, the exponential for the interaction can be expanded in series, yielding successive corrections to the free propagator These expectation values can be computed using the generating functional Eq. (33).
Perturbation expansion for the correlation functions. It is useful, for the following computations, to define a scalar product as (x,y)~x Ã: y, where the asterisk indicates Hermitian conjugation, (: Ã )~ ( (:) T . The generating function can also be written for the non-linear case using the transformed action (37). It is convenient to write it using the real vector J~(j,k)~(j 1 ,j 2 ,k 3 ,k 4 ) as the source term, so that where exp ( : For a small coupling constant m, the exponential in a Taylor series can be expanded to obtain When the function of the path V is a polynomial, every term is the expectation value of the terms of the series expansion of the exponential, and each one can be obtained by differentiating the generating function of the free evolution. The series can be formally exponentiated and written for the generating function of the non-linear case analogously to Eq. (33), that must be normalized by Z(0). The expression for the quadratic generating function can be written as where a zero subscript has been added to indicate that it is the generating function for the linear evolution. Introducing the vector u~(g,w), it is possible to write where D {1 is the Hermitian operator It is possible to obtain an explicit form for Z 0 ½J by inserting u~u c zw, with which the numerator becomes: We can find u c so that D {1 u c {J~0, and then The remaining path integral over w(t) is eliminated by the normalization, therefore the generating function is given by The solution u c can be expressed in terms of the Green's function of the operator D {1 , and the final form of the generating function is This is a general expression; in fact, in the linear case a relation formally identical to the one above is obtained.

Results and Discussion
The Case of the ENSO A simple model of the ENSO system based on the recharge theory was proposed years ago [35]. Following this model, ENSO can be described by a simple linear system The matrix in the equation above is indicated with L.
The action for this system is given by  K 0 (z 1T ,z 2T ,Tjz 10 ,z 20 ,0)~b e bT 2pQ sinh (bT) With the choice of parameters proposed in [35], c~1,c~0:75,r~0:25,a~0:125,b 0~2 :5,m~2=3, the system undergoes stable oscillations, and the entries of the corresponding matrix L are b~0 and w~ffi ffiffiffiffiffiffiffiffiffi 3=32 p . The corresponding propagator can be written as Fig . 2 shows the probability distribution, obtained for a propagator for an initial probability distribution, that is, a delta function at the origin. It is a Gaussian (the figure shows only the section for z 2T~0 ), whose standard deviation increases with time. The system is analogous to a Brownian motion with the particle diffusing in the entire space. The period of the oscillation is close to 20 months and the separate members of the ensemble deviate rapidly as the system evolves. Fig. 3a shows the evolution of the individual members of the ensemble as the oscillation gains larger amplitude. The basic linear oscillation is neutral, so the stochastic fluctuations create the amplification effect, which later will result in the flattening of the probability distribution. For values of m smaller than the critical value, the oscillation is damped, but the stochastic forcing can counterbalance it, permitting a statistical equilibrium. Fig. 3b shows the time evolution for the damped case, and it is possible to see how the divergence is considerably slowed down. Depending on the magnitude of the stochastic force Q, a different value of m is necessary for equilibrium.
The probability distribution is correctly estimated by the propagator as it can be seen in Fig. 4. The zeroth order generating function can be obtained from the Green's function as in Eq. (55). The 2-point correlation function is given by the second functional derivative of Z 0 (J), Considering more derivatives, one might also investigate higher order statistics such as the skewness. The Green's function G 11 for the ENSO model in the transformed coordinates is given by where H here is the sign function. In this way the standard deviation is given by equal time correlations (t~t) Considering the evolution for a semi-infinite domain, when T becomes very large, it will be obtained It is interesting to note that the same time correlation does not depend on the oscillating part of the solution and the frequency w does not appear anywhere. The autocorrelation for positive lags t 0~t {t is given by The cross-correlations in these coordinates are identically zero, but, going back to the (h,h) coordinates, they will recover the correlations shown in [35].
In the same paper [35], a non-linear extension of the standard model is proposed. The non-linear terms represent the negative feedback of the thermocline, and involve the strength of the coupling between the wind stress and the SST; they are cubic in h and h. The extra term appears only in the equation for the temperature as

{E(bhzh) 3
This expression can be used to get the non-linear terms in the action (36) to obtain the perturbation expansion in power of the interaction coefficient E, which corrects the free (linear) propagator in the presence of non-linear terms. The expansion is rather   tedious and, to illustrate the point, the system will be somewhat simplified reducing the non-linear term to a simple form, obtaining a simplified version of the cubic non-linear term in the system (53) that will result in d dt where e, function of b, measures the strength of the non-linearity. The action for this system is given by Eq. (37), where z plays the role of g. The relevant terms in the action are those deriving from w T f(zzr) which, in this case, reduce to the interaction terms between w 2 and z 2 , {ebw 2 g 3 2 . There are also terms deriving from the divergence in the action. The interaction terms are therefore given by: The generating function for these terms is then given by Eq. (44) that can be expanded in power of e, where for convenience the numbering j~(j 1 ,j 2 ) and k~(k 3 ,k 4 ) have been introduced. As one can see from Eq. (33), the functional derivatives have to be evaluated at the same time point, t, and they correspond to the powers of the dynamical variables.
As an example, the correction of the temporal covariance of z 1 will be computed to demonstrate the approach. This covariance is given by the 2-point correlation function, as in Sect. (4), The basic rules of the functional derivation are given by and therefore the two derivatives in Z V (J) will eliminate all terms with less than two j,k, whereas the terms with a larger number of (j,k) will be eliminated by the evaluation at J~(j,k)~0. Due to these mechanisms, the derivative only selects quadratic terms in the expansion of Z V (J). The other term in the first order expansion will be obtained by taking four derivatives, three with respect to j 2 , and one with respect to k 4 . There are two terms of this kind The denominator is given by the following expression (G 24 (t,t)zG 42 (t,t))G 22 ,(t,t)dt : Propagators and interactions can be graphically seen in Fig. 5 and Fig. 6. The numerator is more complicated because now there are two more derivatives. The same arguments used before now lead to the conclusion that only the terms with three Green's functions will survive. The problem is combinatorial and is well known in quantum field theory. It is essentially the same as finding all possible combinations of six points in time: the "external" points, t 1 ,t 2 , and the "internal" points t that are going to be integrated over. Depending on which of the six j or k the derivatives will operate, different kinds of integrals will be generated. The zero order in e is simply G 22 (t 1 ,t 2 ), but for the first order we need to count the contribution from V I . The quadratic term in z 2 will result in The combinatorial analysis indicates that in all there are 4 4! terms given by the four time points; (t 1 ,t 2 ,t,t) are treated, organized in such a way that M 1~1 6 and M 2~8 . More complicated expressions are obtained from the quartic terms. In this case there are three Green's functions involved: G 22 , G 24 and G 42 . Firstly considering the combination with G 24 , it can be seen that there are 5! Ã 3~360 terms,  simplifications can be obtained because the numerator can be factored to the first order in e so that the normalization can be completely canceled at the denominator. The G 22 (t 1 ,t 2 ) can be collected to obtain for the numerator, or at the first order in e The first parenthesis cancels with the numerator and the final expression for the variance is obtained vz 1 (t 1 )z 1 (t 2 )w~G 22 (t 1 ,t 2 )zother terms in e : This is the unperturbed variance corrected by the non-linear terms.
The terms in the perturbation expansion can be expressed with a graphical representation via Feynman diagrams like those in Fig.  5. In this problem there are three kinds of propagators, corresponding to the matrix entries of the Green's matrix. The diagonal entries generate the propagator of the state variable z, and the off diagonal terms, which turn out to be symmetric, generating the propagator connecting the state variable to the auxiliary variables w. The Green's function G 22 (t 1 ,t 2 ) can be graphically expressed with a straight line. On the other hand, the G 24 propagator can be seen as a dashed-continuos line. The points t 1 and t 2 are the external lines of the graph, the time point t is recurring twice and is therefore special, because it has two lines that must be connected with the other point.
The quadratic terms Eq. (60) can be graphically written as in Fig. 7. The (b) graph in the figure represents the integral where the G 22 (t 1 ,t 2 ) propagator can be factored out. It is an example of the fact that these kinds of terms show up graphically since they are made up of separate parts. The so-called ''disconnected'' graph, in this example it is the product of G 22 (t 1 ,t 2 ) and Ð T 0 G 22 (t,t)dt. The terms corresponding to z 3 2 w 2 are more complicated. The internal vertex is of order four and has four lines, which must be connected with two external points. A four line vertex corresponds to the product of two Green's functions, in this case a G 22 and a G 24 , because there are only two external lines. The other two lines must be closed on themselves. The graphs are shown in Fig. 8, without showing all the possible symmetries and exchanges that produce all the 720 terms.
The disconnected graphs are the product of the component graphs, therefore the final correction to the variance or 2-point correlation can be written in the form vz 1 (t 1 )z 1 (t 2 )w~G 22 The results are shown in Fig. 9. The figure shows the time evolution of the variance at equal times t 1~t2 of an ensemble of 2000 numerical simulations. The solid line for the linear case concurs with the theoretical value at equilibrium, Q=2b~8, within the errors. The first order estimate of the non-linear equilibration gives 7.35 and 6.50 for e~0:1 and e~0:3 which are also in concordance with the results.

Conclusions
This paper has shown that the path integral formulation and functional methods can be used for stochastic equations derived from the type of equation of motion that are used to describe the atmosphere and the ocean. These equations pose special complications because the evolution equations are first order in time causing an action that introduces coupling terms between the velocity terms and the forcing function.
This problem prevents a straightforward application of the method as in quantum physics, however, it can be treated by a careful consideration of the boundary conditions. Complications in higher than one dimensions can be treated using the Stratonovich-Hubbard transformation. A perturbation expansion can then be designed for non-linear cases based on the calculation of the generating function for the n-points correlation functions and Feynman diagrams can be introduced.
In this paper the path integral technique is applied to solve a linear simple model and a non-linear one, related to the Climate System, to demonstrate of the power of this tool. Although the technique seems involuted, it could be very easily generalized and could also be the basis for applications to field equations arising in a field theory. This method has only been used with linear and non-linear simple ENSO models, which contain only depending on time variables. The aim of this paper is to stimulate interest in the path integral technique to study the Global Climate System. The authors' hope is to use the formalism of the field variables to face, with this technique, more complicated models, such as applying this method to study general circulation models with noise.