Fokker-Planck equation with memory: the crossover from ballistic to diffusive processes in many-particle systems and incompressible media

The unified description of diffusion processes that cross over from a ballistic behavior at short times to normal or anomalous diffusion (sub- or superdiffusion) at longer times is constructed on the basis of a non-Markovian generalization of the Fokker-Planck equation. The necessary non-Markovian kinetic coefficients are determined by the observable quantities (mean- and mean square displacements). Solutions of the non-Markovian equation describing diffusive processes in the physical space are obtained. For long times, these solutions agree with the predictions of the continuous random walk theory; they are, however, much superior at shorter times when the effect of the ballistic behavior is crucial.


Introduction
Although particle diffusion processes have been studied for about two centuries [1], there are still some subtle issues that need to be faced at the present time. In this paper, we point out these subtle difficulties and offer ways to overcome them.
The quantitative theory of diffusion processes begins in 1855 with the phenomenological solution of the diffusion problem by Fick. Fick employed an empirical definition of the particle flux through the surface of a subvolume (Fick's first law) and the continuity equation which reflects the conservation of particles. This combination results in the diffusion equation (the Fick's second law) which defines the time evolution of the probability distribution function (PDF) of particle concentration. The variance of this PDF grows in time according to where D is the diffusion coefficient which in general depends on the particles and medium in which they diffuse. The equation solved by the PDF is the classical diffusion equation The solution of this equation with the initial condition f (r, t = 0) = δ(r ), where δ(x) is the Dirac delta-function, reads f (r, t ) = exp −r 2 /(4D t ) πD t .
where α 1. We expect α 1 when the diffusion steps are correlated, with persistence for α > 1 and anti-persistence for α < 1 [3]. This law is usually valid only at long times and for t ≪ t c we may again find ballistic behavior. All the issues discussed above for classical diffusion reappear in the context of anomalous diffusion, both in 1 and in higher dimensions. The bulk of this paper will deal with establishing the methods to achieve a consistent theory for the PDF that is valid for all times. As a preparation for more complex situations, in section 2 we review the telegraph equation that regularizes all the issues raised above in 1 dimension. In section 3 we generalize the methodology of the telegraph equation to Fokker-Planck equations with arbitrary non-local kernels in time. In the same section we show how to relate this kernel to the observable mean square displacement of a diffusive process. In the next section we turn to the Langevin equation for a guidance how to compute the mean square displacement to uniquely determine the associated Fokker-Planck equation. In the next Section we generalize the methods used in this paper to the 3-dimensional case taking into account hydrodynamic interactions. The last section offers a summary and conclusions.

Review: the telegraph equation in 1 dimension.
In order to develop a diffusion model with possible finite speed of propagation, it is necessary to simultaneously take into account the mean particle velocity c and the time length t c of the mean free path [12]. Then, the diffusion process is defined by the following equation where the diffusion coefficient D = c 2 t c .
In the limit t c → 0, the time evolution of the PDF is possible only if D 0, e.g., the mean particle velocity c → ∞. In this case, equation (7) is reduced to the classical diffusion equation (2) and the solution of this equation with the initial condition (3) is equation (4). This solution is normalized to unity and its variance agrees with equation (1).
The delta-function initial condition can be considered as a limit of a sequence of continues functions, for example Comparing with equation (2) we conclude that normal diffusion is an inverted process to the limit defined by equation (8), the initial delta-function spreads in space when time increases, keeping the center of mass at rest at the origin.
In the limit t c → ∞, equation (7) transforms to the wave equation Its solution with the initial condition defined by equation (3) reads This solution looks as though the initial conditions propagate with the velocity c without changing the shape of the distribution.
For an arbitrary value of t c , the solution of the telegraph equation is given by three contributions (see, e.g., [15] and the references cited therein) Here, the two dimensionless variables are the normalized time τ = t /t c and the normalized length ξ = r / D t c . Θ(x) is the Heaviside function. The variance of this distribution reads . The delta-function is graphically represented by narrow Gaussians [see equation (8)]. The functions f I ,I I (ξ,τ)Θ(τ − ξ) [see equations (13) and (14)] are shown in panels (b) and (c).

13004-3
This form clearly crosses over from ballistic to normal diffusion as τ increases. The interpretation of the three contributions is as follows: the first term corresponds to the damped ballistic propagation of the kind defined by equation (10). The next two terms exist only inside the compact region that expands with the velocity c. The functions f I (ξ, τ) and f I I (ξ, τ) are given by and where I ν (z) are the modified Bessel functions of the first kind. The asymptotic behavior of these functions is defined by I ν (z) z→∞ ∼ exp(z)/ 2πz, therefore, in the long time limit τ ≫ 1, equation (11) is reduced to the solution of the classical diffusion equation given by equation (4). Various terms in the solution of the telegraph equation, equation (11), are shown in figure 1. Physically, the first term describes the part of the initial mass of diffusers at the origin that shoots out as a propagating solution which decreases exponentially in time, thus providing more and more mass to the diffusive part of the transport process. The function f I (ξ, τ)Θ(τ − ξ) represents the spreading of diffusers whose mass peaks at the origin without participating in the propagation of the first function. Note that this solution is zero at t = 0, but it is finite for any positive time. The third term f I I (ξ, τ)Θ(τ − ξ) is seen to first increase in weight and in amplitude and later on to decrease. It stems from the existence of a finitely rapid front that does not allow the unphysical infinite speed of propagation that characterizes the solution (4). Of course, the sum of the three contributions is normalized.

From the telegraph equation to the time-nonlocal Fokker-Planck equation
In order to motivate the use of the Fokker-Planck equation with time non-local kernels, we rederive the telegraph equation in the following way: let us begin with the continuity equation in the following form: Here, Γ(r, t ) is the particle flux. In its turn, this flux satisfies the equation [12] ∂Γ(r, t ) where Γ 0 (r, t ) is the same as the first Fick's law The solution of equation (16) is given by [17,18] Substitution of this equation into equation (15) yields the time nonlocal diffusion integro-differential equation After time differentiation, this equation is reduced to the telegraph equation defined by equation (7).

13004-4
where W (t ) is the kernel responsible for the non-Fickian behavior of the diffusion process. The Laplace transform of the solution of this equation is given by [19] f (r, s) = 1 where the function P (r, s) is the Laplace transform of the solution of the auxiliary equation with the same initial condition which is identical with equation (2). The Laplace transform of the solution corresponding to the initial condition defined by equation (3) follows from equation (4) where K 1/2 (z) is the modified Bessel function of the third kind. Substitution of equation (23) into equation (21) yields the solution of equation (20) f (r, s) = 2 π r W (s) The Laplace transform of an even moment of the distribution is defined by Substitution of equation (24) into equation (25) yields For the moment of zero order, it is necessary to take into account that lim m→0 mΓ(2m) = 1/2, therefore, 〈r 0 〉 s = 1/s, i.e., the normalization condition.
For the variance (m = 1), equation (26) reads The kernel of the time-nonlocal Fokker-Planck equation is defined by the mean square displacement and, therefore, [see equation (20)] the PDF is also defined by this quantity.

Determining the kernel from the Langevin equation
The conclusion of the last section is that in order to obtain the appropriate Fokker-Planck equation for a given process, we need to determine the time dependence of the variance 〈r 2 〉 t . One way to do so is to measure this moment from experimental data. On the other hand, if this data are not available, or if one wants to derive this information from the physics of the problem, another starting point can be the generalized Langevin equation [6,7].
The standard Langevin equation is Newton's second law applied to a Brownian particle where the random force acting on a particle is taken into account. As was shown in [20][21][22], the generalized Langevin equation is written in terms of a time non-local friction force: where R(t ) is the random force component which is uncorrelated with the velocity that has a zero mean 〈R(t )〉 = 0. The autocorrelation of the random force is related to the kernel in equation (28) by the fluctuation-dissipation theorem [21,22]: where m is the mass of a particle. For the kernel of the kind γ(t ) = γ 0 δ(t ), equation (28) is transformed to the standard Langevin equation. In the general case, the Fourier transform of the random force correlation function is colored, for example cf. [23]).
In dimensionless variables defined in the section 2 with the diffusion coefficient defined by D = k B T t c /m and 1/t c = ∞ 0 γ(t )dt , the Laplace transform of the mean square displacement follows from equation (28) (see, e.g., [24]) and reads where the dimensionless kernel satisfies in the time domain ∞ 0 γ(τ)dτ = 1. One can see from this equation that if lim s→∞γ (s)/s → 0, the function 〈ξ 2 〉 s ∼ 1/s 3 and the mean square displacement exhibit the ballistic behavior at short times 〈ξ 2 〉 t ∼ t 2 . The shape of the functionγ(s) at small s is responsible for asymptotic behavior of the mean square displacement in the time domain.
Substitution of equation (30) into equation (27) yields the kernel This equation settles the relation between the memory kernel of the time-nonlocal Fokker-Planck equation and the memory kernel of the Langevin equation.
The solution of equation (20) given by equation (24) with the memory kernel defined by equation (31) Equation (32) consists of two terms, where The Taylor expansion of equation (34) is defined by For the modified Bessel function of the third kind ∂ n ∂z n Equation (33) It is suitable to rewrite equation (33) in the following form Taking into account that equation (37) can be written as Let lim s→∞γ (s) = γ 0 . Under this assumption, the sum in equation (40) asymptotically tends to zero and the function defined by equation (40) reads and the function f I I (ξ, s) from equation (38) is given by The inverse Laplace transform of equation (42) yields proceeding to limit lim τ→0 f I I (ξ, τ) = 0 shows that this contribution to the PDF has nothing to do with the initial condition and is responsible for diffusion of injecting particles during the transport process.
The nonvanishing term in the function f I (τ, s) in the limit of large s is defined by Its inverse Laplace transform yields Therefore, the part of the PDF under consideration which is responsible for the initial condition and the further impulse propagation is contained in the first summand in equation (38). Now we can isolate from the PDF defined by equation (32) the part corresponding to the diffusion process of the particles which are lost by the propagating impulse where the function Φ(r, s) is defined by equation (32)

13004-7
The inverse Laplace transform of equation (47) yields the diffusive part of the PDF in the time domain The importance of this result is that the explicit Heaviside function takes upon itself the discontinuity in the function f diff (ξ, τ). The function Φ(ξ, τ − ξ) in the time domain is a continuous function. If it does not have an analytical representation in a closed form, it can be evaluated numerically, for example using the direct integration method [25]. Summing together the results (45) and (48) in the time domain we get a general solution of the non-Markovian problem with a short-time behavior, in the form From this solution one can see that the diffusion repartition of the PDF occurs inside the spatial diffusion domain which increases in a deterministic way. The first term in equation (49) corresponds to the propagating delta-function which is inherited from the initial conditions, and it keeps decreasing in time at the edge of the ballistically expanding domain.
In order to estimate the PDF in the long time limit, it is necessary to suggest the form of the memory kernel of the Langevin equation at small s. The reasonable choice is given by One can see from equation (30) This result coincides with the PDF from [27]. Below we demonstrate the application of the time-nonlocal approach with explicit examples.

Standard asymptotic diffusion
For γ(s) = 1 in equation (31) In this case, equation (46) reads where u = s + 1/2, therefore, the inverse Laplace transform is given by The function Φ(ξ, u) for the particular case under consideration follows from the definition given by where From the initial value theorem, it follows that Therefore, the inverse Laplace transform of equation (56) is given by The inverse Laplace transform of equation (57) is given by [28] Ψ(ξ, t ) = I 0 where I 0 (z) is the modified Bessel function.

Substitution of equation (60) into equation (59) yields
Substitution of this equation into equation (49) reduces the solution to the result for the telegraph equation given by equation (11). The diffusion part of this solution [the sum of graphics shown in panels (b) and (c) in figure 1] is displayed in figure 2. As was discussed above in the long time limit, this part approaches the solution of the classical diffusion problem given by equation (2). In order to measure how fast the convergence takes place, it is convenient to estimate the time dependence of the kurtosis of the distribution defined by κ = 〈ξ 4 〉 〈ξ 2 〉 2 − 3.
(62)  For the Gauss distribution κ = 0, for the pure ballistic propagation κ = −2, moments are defined by equation (26) and for the kernel from equation (53) the second moment is given by equation (12). Calculation of the fourth moment yields The kurtosis calculated with these moments is shown in figure 3. At short times, the propagation is ballistic with the following transition to the standard diffusion during a large time interval.

Anomalous diffusion
In order to join the short time ballistic and long time anomalous regimes [see equation (6) and equation (50)], the Laplace transform of the memory kernel of the Langevin equation can be defined by It follows from this equation that and, therefore, the singular part of the PDF is given by The diffusion part of the PDF was estimated numerically and is presented in figure 4. The corresponding kurtosis is shown in figure 5.  In a different way, the kernel of the time-nonlocal Fokker-Planck equation can be defined by the mean square displacement given in the time domain. In [14] this quantity was proposed in a form interpolating the short time ballistic and long time anomalous behavior where 0 α 2 and t 0 is the crossover characteristic time, at t ≪ t 0 , the law (67) describes the ballistic regime and at t ≫ t 0 the fractional diffusion.  Introduce now dimensionless variables 〈ξ 2 〉 τ = 〈∆r 2 〉 t /(2D α t α 0 ) and τ = t /t 0 . With these variables, the last equation reads

13004-11
The Laplace transform of equation (68) is given by where Γ(a, s) is the incomplete gamma function. At large s equation can be written as Substitution of equation (69)  The time evolution of the kurtosises of these distributions is displayed in figure (7). The reader should appreciate the tremendous role of memory. The ballistic part which is represented by the advancing and reducing delta-function sends backwards the probability that it loses due to the exponential decay seen in figure 1, panel (a). Thus, the diffusive part of the PDF is replenished with time and asymptotically approaches the Markovian distribution. The speed of convergence and the evolution of the PDF shape depends on the interim behavior of the mean square displacement (or, which is the same, of the memory kernel of the Langevin equation).

Memory kernel
The friction force acting upon a spherical particle has been derived in [29,30] (the modern discussion can be found in [31]) and is given by where F fr (t ) is the friction force, R is the particle radius, and ρ and η are the density and viscosity of the solvent. Substitution of this force into the Langevin equation yields the fractional equation, and the Laplace transform of the memory kernel can be found in a straightforward way [32]. For simplicity instead of this approach we will use the results obtained [16] for the mean square displacement.
A system with hydrodynamic memory has two characteristic times where M s is the mass of a solvent particle, S = 6πηR is the Stokes friction coefficient and One can see from equation (73) that τ F > 0. Therefore, it is reasonable to introduce dimensionless variables τ = t /τ F and β = τ R /τ F . Then, the projection of the mean square displacement onto one of the axis reads 〈∆x 2 (τ)〉 = 2Dτ F τ − 2 βτ/π 1/2 + (β − 1) where D is the diffusion coefficient, erfc(z) is the complimentary error function and coefficients α 1 and α 2 are given by The dependence of the coefficients α 1 and α 2 on the parameter β is shown in figure 8.  At short times τ ≪ 1, the mean square displacement is given by the Taylor expansion of equation (74) One can see that the mean square displacement defined by equation (76) corresponds to the ballistic motion.
In the long time limit, the mean square displacement is given by 4. β = 9 (M = 0). In this case, the real parameters α 1,2 = (3 ∓ 5)/2. is the value of the parameter β, the slower is the convergence to the asymptotic limit.
It is necessary to note that the predictions of equation (74) are in a good agreement with experimental data [9] The Laplace transform of the mean square displacement given by equation (74) reads 〈 ∆ξ 2 (s)〉 = 2 therefore, the memory kernel of the Langevin equation is given by Substitution of equation (80) to equation (31) yields the kernel of the time-nonlocal Fokker-Planck equation with due regard for hydrodynamics effects W (s) = 1 s + 1 + βs . (81)

Probability density function
Solutions of the telegraph equation in higher than one dimension were considered in [33] and it was shown that the result is unphysical, in some regions the PDF becomes negative. Nevertheless, a reasonable PDF in explicit form was obtained for two-dimensional systems within the framework of a persistent random walk model [34]. Therefore, the possibility to describe the high dimension diffusion process with finite speed of propagation using the time-nonlocal approach is an open question.
In the general case, equation (20) reads in dimensionless units where d is the space dimension and the Laplace operator is given by For pure ballistic propagation, the kernel in equation (82) In the limit τ → 0, the function defined by equation (95) tends to three-dimensional delta-function, i.e., the initial condition of the problem under consideration. In infinitesimal time, the delta-function due to hydrodynamic interactions becomes smooth and the PDF is free from the singular contribution. Therefore, there is no need to pull out the ballistic contribution from the full PDF, and the solution of the problem is given by where the inverse Laplace transform of the continuous function is This can be estimated numerically. The results of calculations for different values of the parameter β are shown in figure 10. Note that a similar ballistic peak of the PDF was observed in simulations of the test particle transported in pseudoturbulent fields [36]. The kurtosis of the PDF for different values of the parameter β are compared in figure 11.

Conclusion
The strategy followed in this paper is to construct a time-nonlocal Fokker-Planck equation which reproduces the time dependence of the mean square displacement of an underlying process throughout the time domain. It should be stressed that the same mean square displacement may correspond, in general, to different models for the time-dependence of the PDF. Thus, the predictions of the time nonlocal Fokker-Planck approach should be compared to experimental data for the PDF to assess its scope of validity and the quality of the approximation. The advantage of the proposed model is that it can deal with diffusion processes that cross-over from a ballistic to a fractional behavior when time increases from short to long values, respectively.
The general one-dimensional solution (49) demonstrates the effect of the temporal memory in the form of a partition of the probability distribution function inside a growing spatial domain which increases in a deterministic way. The approach provides a solution that exists at all times, and, in particular, is free from the instantaneous action puzzle. An extension of the employed approach to higher spatial dimensions is used to study the implications of hydrodynamic interactions on the shape of the PDF. It is shown that singular ballistic contribution to the PDF is smoothed out during the propagation. The expansion given by equation (37) splits the Laplace transform of the PDF into two parts. One part is the ballistic part which is the solution of an inhomogeneous differential equation in spatial variable with the initial condition of the problem. The second one is the diffusion part, which solved the homogeneous equation and is zero at the initial time. In its turn, the diffusion part is divided again into a few terms which correspond, in general, to different auxiliary equations. Therefore, each diffusion term in s-space can be multiplied by a coefficient C (s) so that the normalization condition and the mean square displacement are preserved. Moreover, there is a freedom of modelling here, and this freedom has not been exhausted in this paper. Further analysis and comparison with experimental systems are necessary to nail down the various options of modeling different physical realizations. This research will be the subject of the coming publications.