Optimal FPE for non-linear 1d-SDE. I: Additive Gaussian colored noise

Many complex phenomena occurring in physics, chemistry, finance, etc. can be reduced, by some projection process, to a 1-d SDE for the variable of interest. Typically, this SDE results both non linear and non Markovian, thus an exact equivalent Fokker Planck equation (FPE), for the probability density function (PDF), is not generally obtainable. However, the FPE is desirable because it is the main tool to obtain important analytical statistical information as the stationary PDF and the First Passage Time. Several techniques have been developed to deal with the finite correlation time $\tau$ of the noise in nonlinear SDE, with the aim of obtaining an effective FPE. The main results are the"best"FPE (BFPE) of Lopez et al. and the FPE obtained by using the"local linearization assumption"(LLA) introduced by Grigolini and Fox. In principle the BFPE is the best FPE obtainable by using a perturbation approach, where the noise is weak, but the correlation time can be large. However, when compared with numerical simulations of the SDE, the LLA FPE usually performs better than the BFPE. Moreover, the BFPE gives often"unphysical"results that reveal some flaws, problems that do not affect the LLA FPE. The common step of the perturbation approaches that lead to the BFPE is the interaction picture. In this work, that is presented in two companion papers, we prove that this issue affecting the BFPE is due to a non correct use of the interaction picture, a consequence of the pitfalls of non-linear dissipative systems. We will show how to cure this problem, so as to arrive to the real best FPE obtainable from a erturbation approach. However, the LLA FPE for 1d-SDE continue to be preferable for different reasons. In this first paper we shall consider nonlinear systems of interest perturbed by additive Gaussian colored noises.

In the present work we are interested in non-linear 1-d SDEs that can be written as: X = −C(X) + I(X)ξ(t). (1) where X is the variable of interest, −C(X) is the unperturbed velocity field, I(X) is the perturbation function, ξ(t) is the stochastic perturbation with zero average and autocorrelation function ϕ(t) = ξξ(t) / ξ 2 , the parameter controls the intensity of the perturbation.
Countless are the cases in which some physical (and not only) process can be described by a SDE like that in Eq. (1). To avoid making confusing the reader, in this paper we shall focus our attention on the additive case, namely when I(X) = 1: The more general case of multiplicative perturbation presents some specific critical issues that shall be addressed in the companion paper.
It is a standard result in Statistics that when the stochastic forcing is a white noise: , the SDE of Eq. (2) is completely equivalent to the following Fokker Planck Equation (FPE) for the Probability Density Function (PDF) P (X, t) of the X variable (we shall use the shorthand ∂ X := ∂/∂X): ∂ t P (X, t) = ∂ X C(X) + 2 ξ 2 ∂ 2 X P (X, t).
From the FPE of Eq. (3) the stationary PDF is given by P W,eq (X) = 1 Z e − X C(y) 2 ξ 2 dy (4) in which Z is the normalization constant. However, white noise is usually a too extreme oversimplification of real driven forcing of phenomena of interest. The importance of systems driven by colored noise has been recognized in a number of very different situations, e.g., statistical properties of dye lasers [7][8][9][10] chemical reaction rate [11][12][13][14], optical bistability [15,16], large scale Ocean/Atmosphere dynamics [17,18] an many others. Here we shall assume that the Gaussian stochastic process ξ(t) is characterized by a "finite" correlation time τ and by the unitary intensity ξ 2 τ = 1. It is well known that both in the present one-dimensional case and in the more general multidimensional one, if the unperturbed velocity vector field is linear, the Gaussian property of the (generally colored) noise ξ(t) is "linearly" transferred to the system of interest, thus the FPE structure does not break (see, e.g., [14,19]). On the contrary, in the case of non linear SDE and/or non Gaussian noise, for not vanishing τ the FPE breaks down. This is the case we are interested in this work. Several techniques have been developed to deal with the correlation time of the noise in nonlinear SDE, with the aim of obtaining an effective FPE that, with a good approximation, describes the evolution and the stationary properties of P (X, t).
They can be summarized in two main strands that correspond to two general techniques: the cumulant expansion technique [20,21] and the projection-perturbation methods (e.g., [22,23]). Each of these methods leads to a formally exact evolution equation for the PDF of the driven process. At this level the different descriptions are therefore equivalent. The exact formal results do not lend themselves to calculations nor give an FPE structure, therefore they require that approximations be made. The approximations made within these various formalisms involve truncations and/or partial resummations of infinite series in the and τ parameters. Not surprisingly, it has been argued and shown [2] that the effective FPE obtained from these different techniques are identical at the same level of approximation (time scale separation, weak perturbation, Gaussian noise etc.). In some recent works we have shown as it is possible to obtain an analytical result for a FPE for weak noises but without restrictions on the τ parameter [24][25][26][27][28]. In the next section we will shortly summarize this result.

THE FPE
From Eq. (2) it follows that, for any realization of the process ξ(u), with 0 ≤ u ≤ t, the time-evolution of the PDF of the whole system, that we indicate by P ξ (X; t), satisfies the following PDE: in which the unperturbed Liouville operator L a is L a := ∂ X C(X) (6) and the Liouville perturbation operator is A standard step of the perturbation method is to introduce the interaction representation, by which Eq. (5) becomes andL  (10) is also called the Lie evolution of the operator L I along the liouvillian L a , for a time −t. For a further use, we note that the Lie evolution of a product of operators is the product of the Lie evolution of the single operators: Assuming that at the initial time t = 0, the PDF P (x; 0) does not depend on the possible values of the process ξ (or that we wait enough to makes ineffective the initial conditions) from Eq. (8), to the first non vanishing power in the parameter , we get (see Appendix ??): from which, getting rid of the interaction picture, we have This is the standard result we obtain with any perturbation approach (for example, the Zwanzing projection method [24]). The next step is to rewrite, if possible, Eq. (12) as a FPE: Then, given the state dependent diffusion coefficient D(X), the stationary PDF of the FPE is easily obtained To pass from Eq. (12) to Eq. (13), the crucial term is the operator e L × a u [∂ X ]. In most of the literature concerning the Zwanzing projection method (e.g., [22]), the explicit FPE is obtained from Eq. (12) by assuming that τ , the decay time of the correlation function ϕ(t), is much smaller than the unperturbed dynamics driven by the liouvillian L a . In fact, in this case it is save to substitute, in Eq. (12), the power expansion (note the shorthand that leads to a FPE with a state dependent diffusion coefficient, given by a series of "moments" of the time u, weighted with the correlation function ϕ(u). However, such a series, as it is apparent from Eq. (15), contains secular terms and is not (generally) absolutely convergent. This is clearly shown in the example considered in Fig. 1. A way to avoid this problem is solving, without approximations, the Lie evolution of the differential operator ∂ X along the Liovillian L. In [28] this has been done for the general case of multidimensional sys-tems and multiplicative forcing. In the present simpler one-dimensional case, recalling that L a = ∂ X C(X), the Lie evolution of ∂ X , without approximations, can be directly obtained as follow: where • the Lie evolution along a deterministic (first order partial differential operator) Liouvillian of a regular function, here C(X), is just the back-time evolution of the same function along the flux generated by the same Liouvillian, Inserting Eq. (16) in the FPE of Eq. (12) we get, in a clear and straight way, the so called Best Fokker Planck Equation (BFPE) of Lopez, West and Lindenberg [2] (actually, our is a generalization of that, because we don't make the assumption that ϕ(t) = exp(−t/τ ) and we leave finite the time integration t [29]): For "enough" weak noise ξ(t), the BFPE should be the best possible approximation we can obtain from a perturbation approach to the SDE of Eq. (2). Actually we shall see in the following that some more prescriptions must be added to the formal expression of Eq. (18).
If we are interested in short times statistical features of the system, we have to leave the time t as the upper limit of integration in Eq. (18). However, if we are interested in large times (compared to the time scales of both the unperturbed system of interest and the correlation of the perturbation), or to stationary statistics, in order to avoid the divergence of the same integrals, we must check, case by case, the decaying property of the integrand. Actually, assuming that τ is not a small expansion parameter, we should expect that, as for any result from a time perturbation procedure, the BFPE describes well the dynamics of the PDF for times much smaller than some timet, proportional to the inverse of the expansion parameter 2 . For t ≥t, in principle, the BFPE could completely fail. Unfortunately, the relaxation process emerging after the coupling of the system of interest with the perturbation, takes times that generally are larger thant (order of ( 2 τ ) −1 ). Therefore, the sole perturbation approach cannot ensure, for example, that the stationary PDF (= lim t→∞ P (X; t)) of the FPE of Eq. (18) is a good approximation of the "true" one. But the cumulant approach can do that. In different perspectives and contexts this fact has been largely discussed and proved by van Kampen [30,31], Fox [32], Terwiel [33], and Roerdink [34] many years ago, working directly (and heavily) with the analytical expressions of the generalized cumulants defined by the t-ordered exponential. These results have been generalized in [6], where it is shown that for the SDE of Eq. (2), the projection perturbation approach leads to the BFPE of Eq. (18) plus terms that destroy the FPE structure that are at least order O( 4 τ 3 ).
Concerning the LLA FPE, West at al. have shown [2] that it can be formally derived from the BFPE of Eq. (18) in the following way: 1. by assuming that there is enough time-scale separation between the unperturbed dynamics and the decay time of the correlation function ϕ(t), so that the unperturbed dynamics X 0 (t − u) can be considered close to the initial position X; 2. given point 1 here above, instead of expanding directly 1 C(X 0 (t−u)) in powers of u (that would give rise to the same secular terms of the expansion given in Eq. (15)), West at al. expand its logarithm: truncating the series at the first order.
In fact, by using point 2 in Eq. (18), we are lead to the LLA FPE (generalized to finite times and to general correlation function of the noise): It is important to stress that in the case of linear systems of interest, namely for C(X) = γX, the series expansion of the r.h.s. of Eq. (19) stops exactly at the first order in u, while this does not happen expanding directly the term 1/C(X 0 (t − u)). Therefore, instead of using the West and al. pathway (represented by the points 1-2 here above) to go from the BFPE to the LLA FPE, the latter one can be though as the result of the approximation given by substituting the function C(X)/C(X 0 (t−u)) with and exponential decay function with state dependent decay coefficient: we get the following result for the state dependent diffusion coefficient of the FPE: that, for large times becomes in which the "cap" means Laplace transform. From Eq. (22) it results that while D LLA (X, ∞) exists and it is positive under quite general and clear conditions, the situation is quite more complex for D BF P E (X, ∞). A simple example may serve for illustration. Le us consider the case in which C(X) = α sinh(kX) and ϕ(t) = exp(−t/τ ). A straightforward calculation leads to C(X)/C(X 0 (t − u)) = cosh(αku) − cosh(kX) sinh(αku), that inserted in Eq. (22), with the convergence constraint αkτ < 1. For the D LLA of Eq. (22) we easily get where now the only constraint is that the dissipative flux is not divergent (namely, α > 0).  23) and (24) and exploiting Eq. (14), with a little algebra we arrive to the following stationary PDFs:  and From Fig. 3 we also see that the LLA stationary PDF given by Eq. (26) that, for t > 2X 2 is a complex number: for large times it is not defined! On the other hand, in this case the LLA diffusion coefficient is simply given by that, for large times, has the following simple limit: Hereafter we shall prove that the above mentioned flaws of the BFPE are due to a non correct implementation of the perturbation procedure, and we shall remedy to this situation.
For this purpose, we note that the negative values of the D BF P E of Eq. (22) is related to the fact that the kernel of the integral can be negative for some X value. For example, considering again the case of C(X) = α sinh(kX), we see from Fig. 4, solid lines, that, after a given timeū(X) that depends on X, the function C(X)/C(X 0 (t − u)) becomes negative. We also see that the larger the X value, the shorter the timeū(X). Thus, whatever the correlation decay time τ ∈ (0, 1/αk), there will be always someX value such that D BF P E (X, ∞) of Eq. (the case for C(X) = X 3 is detailed in Fig. 5). For "preceding" times −u with u >ū(X) there aren't points in the state-space that are connected to X by the flux generated by the velocity field −C(X). This is obviously due to the strong irreversible nature of the flux, that shrinks the state-space. At the end, this implies that for such strong dissipative fluxes, the back time evolution must be limited to times u <ū(X), namely we have to multiply by the Heaviside function Θ(ū(X)−u) any function of X 0 (t−u).
Therefore, the BFPE state dependent diffusion coefficient of Eqs. (21)-(22) must be corrected as in the following: Concerning the stationary PDF, the corrected BFPE result is obtained by inserting in Eq. (14) the here above corrected BFPE state dependent diffusion coefficient.
Thus, for the case C(X) = α sinh(kX), in place of the state diffusion coefficient (23), from Eq. (30) we get : By substituting in Fig. (2) the curves relative to D BF P E (X, ∞) with the corresponding ones relative to D cor BF P E (X, ∞) of Eq. (31), we obtain Fig. (6) where we can see that the state dependent diffusion coefficient D cor BF P E (X, ∞) is always positive and it is also not so much different from the D LLA (X, ∞). The stationary PDF for this case is given by using Eq. (31) in Eq. (14). Because of the integral in the exponent in Eq. (14), we cannot obtain an analytical expression, however, the numerical integration is easily done and the results, for different τ and values, are shown in Fig. (7). We can see that the stationary PDFs of the corrected BFPE are now quite close (within the limits of the perturbation approach) to those from the numerical simulations, also for large τ values and relatively large . However, the LLA FPE continues to perform better.
In the case of a cubic velocity field, namely for C(X) = X 3 , D cor BF P E (X, ∞) of Eq. (30) is now real (it is not a complex number as for D BF P E (X, ∞) ): where F (x) := e −x 2 x 0 e y 2 dy = e −x 2 √ π 2 erfi(x) is the Dawson function. As for the previous case, the corrected BFPE diffusion coefficient of Eq. (32) is positive and close to the corresponding LLA result (see Fig. 8). Also in this case, inserting in Eq. (14) the expression for the diffusion coefficient, here in Eq. (32), we don't obtain an analytical result for the stationary PDF, thus we have to resort to the numerical integration. The result is shown in Fig. 9. As we can see, for small τ values the corrected BFPE leads to a stationary PDF that is close to the LLA one, and both are good approximations of the PDF obtained by the numerical simulation of the SDE of Eq. (2). However, in this case, for τ > 1, both the corrected BFPE and the LLA FPE give a stationary PDF that is no longer so close to the one obtained directly from the SDE. This is not so surprising because the pure cubic velocity field is an extreme case of non-linear system. Despite this, in the tail of the PDF, the LLA result is close to the numerical simulations, while the corrected BFPE PDF remains far (see Fig. 10). This fact is worth stressing because the tails of the PDF affect important statistical observables as the first passage time or/and the waiting time distribution.
Finally, we take into account the diffusion in a periodic potential. In particular we choose C(X) = α sin(kX). In this case, the analytical results for both for the diffusion coefficient and for the stationary PDF are almost identical to the case of hyperbolic velocity field C(X) = α sinh(kX) here above already analyzed, provided we substitute the hyperbolic functions of X with the corresponding trigonometric ones. The main difference is that now the divergence of the flux is not asymptotically growing, thus the back time evolution X 0 (t− u) has not asymptotes (it does not diverges at all). In this case the function C(X)/C(X 0 (t− u)) is always positive, and simply increases with u as e kαu (see Fig. 11). Therefore for C(X) = α sin(kX) the "standard" BFPE formula of Eq. (22) for the diffusion coefficient can be used without need of corrections. It is clear that the periodic nature of the velocity flux do not allow any (not vanishing) stationary PDF. However, if the intensity of the stochastic perturbation is weak, starting from an ensemble located in one well, the system is subjected to a first relatively fast relaxation toward a metastable state, where the PDF is confined in the same well, then, for large times, an infinite diffusion process toward the other wells takes place. In Fig. 12 we see that for τ = 0.95 and = 0.5 there is a large difference between the BFPE and the LLA results. In particular the semi-log plot (see the inset) highlights the good agreement in the tails of the PDF, between the LLA PDF and the PDF obtained from the numerical simulation, while the BFPE result is not so good.

CONCLUSIONS
By definition, the BFPE is the best FPE we can get from a perturbation approach, starting from the SDE of Eq. (2), in which is the small parameter. The BFPE, although restricted to red noise perturbations, was obtained many years ago by Lopez, West and Lindenberg [2], but it present some flaws. In this work we have shown that a correct use of the interaction picture cures these problems of the BFPE. So revisited, the BFPE results close to the FPE obtained by the LLA of Grigolini [3,4] and Fox [5], in particular for weak noise perturbations. However, increasing the intensity of the perturbation, the differences become relevant: the LLA FPE performs better than the "cured" BFPE when compared with numerical simulations of the SDE of Eq. (2). Actually the LLA FPE seems to work well also for values of far exceeding the limit imposed by the perturbation approach. A clear explanation for this is beyond the scope of the present paper, and will be detailed in a forthcoming one. Here we just mention two facts: 1. the LLA approach is based on the assumption that, for any X value, we can safely substitute the unperturbed evolution of the function f (X, u) := C(X)/C(X 0 (t − u)) with an exponential time decays, with X dependent exponent: For one-dimensional dissipative systems, the exponential behavior of such a back time evolution is typical (see Fig. 13); 2. it is possible to rigorously demonstrate that when ξ(t) is a Gaussian stochastic process, the LLA makes vanishing all the terms of the projection perturbation approach that destroy the FPE structure: there aren't O( n ), n > 2 corrections to be considered.
Therefore, although we have obtained the correct BFPE for the SDE of Eq. (2), the sim- . It is apparent that at the time valueū where the back time evolution X 0 (t−u) diverges, the function C(X)/C(X 0 (t − u)) vanishes, while for larger times it is a purely imaginary number.  PDF obtained from Eq. 14) and by using D(X) = D cor BF P E (X, ∞) of Eq. (31). The dashed light gray lines correspond to the previous P BF P E curves, namely those from Eq. (23). Note that the corrected BFPE diffusion coefficient is always positive and behaves "well" also for αkτ > 1.    [36] G. P. Tsironis and P. Grigolini, Phys. Rev. Lett. 61, 7 (1988).