Langevin equation in complex media and anomalous diffusion

The problem of biological motion is a very intriguing and topical issue. Many efforts are being focused on the development of novel modelling approaches for the description of anomalous diffusion in biological systems, such as the very complex and heterogeneous cell environment. Nevertheless, many questions are still open, such as the joint manifestation of statistical features in agreement with different models that can also be somewhat alternative to each other, e.g. continuous time random walk and fractional Brownian motion. To overcome these limitations, we propose a stochastic diffusion model with additive noise and linear friction force (linear Langevin equation), thus involving the explicit modelling of velocity dynamics. The complexity of the medium is parametrized via a population of intensity parameters (relaxation time and diffusivity of velocity), thus introducing an additional randomness, in addition to white noise, in the particle's dynamics. We prove that, for proper distributions of these parameters, we can get both Gaussian anomalous diffusion, fractional diffusion and its generalizations.


General condition for the emergence of anomalous diffusion
Diffusion is described through the following simple, but general stochastic equation: being V t a stochastic process describing a generic random fluctuating signal.
Here X t and V t are the position and velocity of a particle moving in a random medium, respectively. For a generic, nonstationary process, the two-time Probability Density Function (PDF) p(V 1 , t 1 ; V 2 , t 2 ) depends on both times t 1 and t 2 . Similarly, the correlation function is, in general, a function of the times t 1 and t 2 1 . 5 Now, by integrating in time the above kinematic equation (1), making the square and the ensemble average, we get the Mean Square Displacement (MSD): where, in order to get : X t = X 0 , we assumed a uniform initial position X 0 . In the stationary case, the two-time statistics, including the correlation function, depends only on the time lag t = |t 1 − t 2 |, and the above formula reduces to: or, equivalently: where R(t) = V t1+t · V t1 = V t · V 0 is the stationary correlation function.
Notice that these expressions have very general validity, independently of the particular statistical features of V t . These expressions were firstly published by Taylor in 1921 [1], which implicitly formulated the following: 10 Theorem (Taylor 1921) Given the stationary correlation function R(t), let us define the correlation time scale: Then, if the following condition occurs: normal diffusion always emerges in the long-time regime: thus defining the long-time spatial diffusivity D X : independently from the details of the microdynamics driving the fluctuating velocity V t .
Taylor's theorem gives in Eq. (6) the general conditions to get normal dif- 15 fusion, i.e., a linear scaling in the variance: X 2 ∼ t). This result has a very general validity, independently from the statistical features of the stochastic process V t . The theorem also establishes the regime of validity of normal diffusion, given by the asymptotic condition t τ c . As a consequence, the emergence of anomalous diffusion is strictly connected to the failure of the assumption (6). 20 In particular, we get two different cases: • Superdiffusion: (10) • Subdiffusion: In order to get τ c = 0 and, thus, subdiffusion, velocity anti-correlations must emerge. This means that there exist time lags t such that R(t) < 0 (e.g., the anti-persistent Fractional Brownian Motion, with H < 0.5). Being R(0) = V 2 st > 0, in subdiffusion the correlation function is surely positive in the short-  (6) is the main guiding principle exploited here to derive stochastic models for anomalous diffusion.

30
The Fractional Brownian Motion (FBM) B H (t) was introduced by Mandelbrot and Van Ness in their famous 1968's paper [3]. Since then, thousands of papers have been devoted to both theoretical investigations and applications of FBM (see, e.g., [4] for a review). FBM is a Gaussian process with self-similar stationary increments and long-range correlations. In formulas, FBM has the following 35 properties: 3 Interestingly, this relation is here derived in a very general framework, i.e., for a generic fluctuating signal Vt, with the only assumption of the existence of a stationary regime in the long-time limit. As known, the stationary condition usually emerges in correspondence of motion reaching an equilibrium state. However, the stationary condition is more general with respect to equilibrium and, for this reason, we prefer to leave the notation "st" for "stationary" instead of "eq" for "equilibrium". 4 A correlation time scale, different from the above definition of τc can be sometimes introduced for subdiffusion (e.g., the time period in a harmonic correlation function), but it does not have the meaning of discriminating a long-time regime with normal diffusion from a short-time regime.
• B H (t) has stationary increments; • B H (t) has a Gaussian distribution for t > 0;

40
• the correlation function is given by: The FBM increments are given by: The process V δt (s) is also called fractional Gaussian noise 5 . Both B H (t) and V δt (s) are self-similar stocastic processes but, at variance with B H (t), the increments V δt (s) are also stationary, i.e., their statistical features do not depend on s, but only on δt. V δt (s) is a Gaussian process and is uniquely defined by the mean, variance and correlation function, which are derived from the above listed properties of FBM: Then, we can say that FBM is a Gaussian process with stationary and selfsimilar increments V s (δt), while FBM is Gaussian, self-similar but not stationary. Eq. (14) also shows that, with the exception of the standard Brownian motion (H = 1/2), increments V δt (s) are not independent each other. Fractional Gaussian noise and FBM are exactly self-similar, i.e., they satisfy the 45 relationship: X(at) = a H X(t), the increment V 1 (s) with δt = 1 is usually considered in both theoretical and experimental studies, as a generic δt can be obtained by simply rescaling the process with the self-similarity relationship. In Fig. 1 the increment correlation functions of a persistent (H > 0.5) and of an antipersistent (H < 0.5) FBM are compared. It is evident that antipersistent 50 FBM is associated with anticorrelations, and this is the reason why subdiffusion emerges in this case. The asymptotics of the correlation function are easily obtained by rewriting it in the following way (see [4], pages 6-7): being, for x = δt/t < 1: The limit t → ∞ corresponds to x → 0 and the Taylor expansion of h H (x) gives: so that [3,4]: Regarding the correlation time τ c defined in Eq. (5), we can exploit the same asymptotic expansion used for R(t). Firstly, we apply Eq. (5) to a finite time t: so that: τ c = lim t→∞ τ c (t). Then, for the fractional Gaussian noise we get: Analogously to R(t), this can be written as: and, for x < 1, h H+1/2 (x) is again given by Eq. (16), but with H + 1/2 instead of H. Then, an asymptotic formula similar to Eq. (17) can be derived: and, finally: Clearly, the mathematical limit t → ∞ corresponds to the physical regime t δt. Exploiting the asymptotic behavior of τ c (t) given in Eq. (23), we can now derive the values of the correlation time scale τ c = τ c (∞): The three cases correspond to persistent (superdiffusive) FBM, normal Brownian motion and antipersistent (subdiffusive) FBM, respectively.

Box 2. Properties of R(t) and g(τ )
The use of Laplace transform, defined by the expression: gives important information about the normalization and moments of distributions. The stationary correlation function R(t) and the distribution g(τ ) are related by Eq. (8,MT). For any choice of the distribution g(τ ), the correlation function R(t) and g(τ ) must satisfy the following properties: (i) The distribution g(τ ) must be a PDF normalized to 1: which determines a constrain on the behavior of the first derivative of the correlation function: (ii) The MSD is a power-law of time with superdiffusive scaling 1 < φ < 2 in the asymptotic long-time limit: where C1 and C2 are proper constants and the second asymptotic limit follows from the Tauberian theorem [5]. From Eq. (3) or Eq. (4) it results: we get equivalently the following expression for the stationary correlation function: with η = 2 − φ, 0 < η < 1. Note that the above limits can be equivalently written as asymptotic behaviors, e.g.: R(t) ∼ t φ−2 for t → ∞, which means that the function R(t) is approximated by C3t φ−2 in the long time range.
(iii) The MSD at time zero is zero: (iv) Furthermore being 0 < R(0) < ∞, from Eq.(9,MT) the distribution g(τ ) must have non-zero, finite mean: The properties that must be satisfied by the stationary correlation function R(t) and by the PDF g(τ ) are listed in the above Box 1.
We now prove the following where τ * must be introduced to get an adimensional parameter as argument of L −η η . The mean correlation time is given by: so that we have: The normalization constant C can be obtained by imposing L[g(τ )](0) = 1. Exploiting the relationship ∞ s exp(−ξτ )dξ = exp(−ξτ )/τ and making the change of variables τ = ( τ /C)τ , we get: and: Substituting this relationship into Eq. (28) we finally get Eq. (14,MT), which is a properly normalized PDF.
(ii) superdiffusive scaling: We now prove that R(t) ∼ t −η , with 0 < η < 1, a condition leading to the superdiffusive scaling for the position variance: This can be proven thanks to the integral representation of the extremal Lévy density: Hence, we have: where τ * = τ Γ(1/η)/η. It is useful to rewrite the expression as: which can be solved through the residues theorem considering the poles s/η+1 = −n or s = n, with n = 0, 1, 2..∞.

70
In the first case we have: where each term of the series is obtained by the limit: When t → ∞ only the first term survives and we find: Substituting τ * = τ Γ(1/η)/η, we finally get Eq. (15,MT), from which we obtain the superdiffusive scaling of the position variance σ 2 Considering the poles in the other semi-plane, s = n with n = 0, 1, 2..∞, we find that: converges to R(0) = ν τ , as already shown before.
(iii) MSD at time zero is zero:

Example:
In the special case η = 1/2, the extremal Lévy function corresponds to the Lévy-Smirnov distribution, the whole exercise can be solved analitycally and we may consider for simplicity τ Γ(1/η) η = 1 : Solving the integral the analytical form of the correlation function turns to be: which leads to the following exact formula for the position variance: which, apart from the conditional statistics, is essentially the same as Eq. (9). For a standard OU process with fixed ν and τ , V 2 |V 0 , τ, ν st = V 2 eq and Eq. (41) relates the diffusion (D X ) and relaxation (τ ) properties through the equilibrium distribution ( V 2 eq ). In his 1905 paper [2], Einstein studied the Brownian motion in a gas at equilibrium, where velocity distribution is given by the Maxwell-Boltzmann law. In this case, the Einstein-Smoluchowsky relation becomes: being T , m and k the gas temperature, the Brownian particle mass and the Boltzmann constant, respectively.

Numerical scheme for the Langevin equation
In order to avoid stability problems, the numerical algorithm for the simulation of Eqs. (1) and 4,MT) was implemented using an implicit scheme with order of strong convergence 1.5 [6]. This is given by the following expression: being V n = V (n∆t), ∆t the time step, ∆W n = W (t n + ∆t) − W (t n ) the increments of the Wiener process, a(V ) = −V /τ and b = √ 2ν the drift and noise terms, respectively. Further, we have: being u 1 (n) and u 2 (n) two independent random numbers with uniform distributions in [0, 1]. A suitable time step ∆t, also depending on the time scale τ , is necessary to maintain the accuracy of the numerical scheme. To take into account both the ensemble variability of the relaxation time τ , which is different for different trajectories, and the time variability of drift and noise terms along the same trajectory, we applied a variable time step according to the scheme given in Ref. [7]: This adaptive time step allows to avoid any problem of convergence and accuracy in the numerical scheme, Eqs. (43) and (44). At the same time, in the range of short τ , this algorithm can give very short time steps, thus determining very long simulation times for a consistent number of trajectories. To overcome this 85 problem we note that the short time regime τ τ of the PDF g(τ ) does not significantly affect the anomalous scaling of diffusion, which mostly depends on the asymptotic tail of the distribution g(τ ). A cut-off was then introduced in the short-time regime. By comparing the numerical simulations with theoretical results we chose the cut-off value τ min = 0.004, much smaller that τ , which is 90 always of the order 0.5 − 1 for all sampled sets of τ .

Numerical algorithm for the random generator of τ
Here we describe a method to generate random variables τ distributed according to the law of Eq. (14,MT), where A(η) is the normalization coefficient, and τ is already dimensionless. For this, we use a well-known inverse transform sampling method (see, e.g. [8]), so the procedure is straightforward.

95
First, we generate a set of extremal Lévy density random numbers L −η η (τ ) by using the generator described in Refs. [9,10], see Eq. (3.2) of the latter paper, and extract its histogram. Since the beginning of the histogram has much statistical noise (red curve in Fig. 2a), it is a good solution to replace these values with analytical asymptote at small arguments [11] (blue curve in Fig. 2a). Moreover, we also expand the histogram with another asymptote, at large τ s (green curve in Fig. 2a): where Then, we divide the obtained histogram by argument and find the normalization coefficient numerically in order that the resulting PDF is normalized to unity. Finally, we calculate the semi-analytical cumulative distribution function (CDF) (see Fig. 2): where δτ i is the i th histogram's bin width, i = 0, 1, 2..n. Now, we draw a random variable τ obeying the target pdf (46) with where u ∈ [0, 1) is a uniformly distributed random variable: F −1 is a numerically (or if u > F (x n ), semi-analytically) inverted CDF. Let us take out a verification and compare the original PDF g(τ ) used for the simulations and the histogram of the generated 10 7 random numbers with 100 this algorithm g sim (τ ). The result is shown in Fig. 3. At intermediate values of τ the inaccuracy is about 1%, increasing due to statistical error at very small and large τ s (where g(τ ) is small).
The software for the numerical simulations were written in C++ language (Debian gcc 4.9) and Python 2.7 and can be downloaded at the following web-site:  110 We here provide an intuitive presentation of the Schneider grey noise, the grey Brownian motion and the generalized grey Brownian motion. More rigouros details can be found in [12,13,14,15,16,17,18,19].

Schneider grey noise, gBM and ggBM
The grey noise is a generalization on the basis of the Mittag-Leffler function of the white noise. The Mittag-Leffler function E β (z) is defined as and it is a generalization of the exponential function that is recovered as special case when β = 1, i.e., E 1 (−z) = e −z . As well as the exponential function, when 0 < β < 1, the Mittag-Leffler function is a completely monotonic function. A useful formula for what follows is For any characteristic functional Φ(z) there exists a unique probability measure µ such that and if Φ(z) = E β (−z 2 ), 0 < β < 1, the probability measure µ is the so-called Schneider grey noise [12, 13,17]. When β = 1 we have E 1 (−z 2 ) = e −z 2 , and 115 the Gaussian white noise follows. Let us introduce the stochastic process X(t) driven by the noise µ and we look for its probability density function. The characteristic function is where function ϕ α (t) takes into account what remains of parameter t after the integration, and it is related to the scaling in time of X(t). By the inversion of (58) we have the probability density function of X(t) as follows where M β/2 is the M-Wrigth/Mainardi function. By using (58) and (56), we have that the variance of X(t) is In the same spirit, the correlation function of the process X(t) can be computed. In fact from (58) it holds and by applying again formula (56) the correlation function results to be Now we discuss how to establish function ϕ α (t). Let 1 [a,b] be the indicator function such that it is equal to 1 when a < t < b and to 0 elsewhere. In analogy with the Wiener process where the Brownian motion is B(t) = t 0 dW (τ ), we write the process X(t) as where X 0 is a random variable equivalent in distribution to X(t) but independent of t, i.e., the probability density function of X 0 is p 0 (x) = p(x, t = 1). From (63) we have that and from comparison with (60) and (62), we obtain that ϕ α (t) is established through the stochastic process 1 [a,b] that meets Finally we observe that, by setting ϕ 2 α (t) = t α , X(t) is the Brownian motion when α = β = 1, and we refer to it as the grey Brownian motion and the generalized grey Brownian motion when 0 < α = β < 1 and 0 < α < 2, 0 < β < 1, respectively. Moreover, in order to have a process with stationary increments we assume ϕ 2 α (t, s) = |t − s| α , the correlation function results to be The corresponding stochastic process is obtained with a randomly-scaled Gaussian process, i.e., a Gaussian process multiplied for a non-negative independent randon variable not dependent on time.
From integral representation formulae of the M function [20], we have that X 0 has the same density of X(1) if, for example, we state X 0 = √ Λ B(1) where Λ is a non-negative random variable distributed according to M β and B(1) is a Gaussian variable. Finally, we obtain that Looking at (60) and (62), the process Finally, by setting H = α/2, the trajectories of the process X(t) can be generated by Since the fBm X H (t) is fully characterized by the variance and the correlation 120 functio, the process X(t) is also fully characterized by the variance and the correlation function. With a somewhat forced terminology, the term ggBM can be thought to include any randomly scaled Gaussian process, i.e., any processes defined by the product of a Gaussian process with an independent and constant non-negative random 125 variable.

Mainardi distribution and Lévy densities
Fractional diffusion processes are a generalization of classical Gaussian diffusion, mainly in the direction of the time-fractional diffusion, i.e., by replacing the first derivative in time with a time-fractional derivative, and in the direction of the 130 space-fractional diffusion, i.e., by replacing the second derivative in space with a space-fractional derivative. In the case of time-fractional diffusion the Gaussian particle density is generalized by the so-called M -Wright/Mainardi functions [22,23], and in the case of the space-fractional diffusion the particle density is generalized by the so-called Lévy stable densities [11].
The M-Wright/Mainardi function M ν (r), r ≥ 0, 0 < ν < 1, is defined by the series: and it provides a generalization of the Gaussian and Airy functions: Moreover, the following limit holds: The M density function is related to the Mittag-Leffler function through the Laplace transform: and it has an exponential decay for r → ∞, i.e.: which allows for finite moments that can be computed through the formula: ∞ 0 r q M ν (r)dr = Γ(q + 1) Γ(νq + 1) , q > −1 .
The asymptotic behaviour for |z| → ∞ is the power-law and, for extremal densities, the following exponential decay holds for z → 0: Moreover, the following limit holds: A remarkable formula of the Lévy density is the following integral representation for z ≥ 0, 0 < β < 1: that, in the special case α q = 2, θ q = 0, provides the following link with the Gaussian density [20,24]: The M ν (r) function, r ≥ 0, 0 < ν < 1, and the extremal Lévy density L −ν ν (r) are related by the formula: In the present paper we consider such special densities in order to highlight the relation of the proposed formulation with the fractional diffusion. However, the asympototic behaviour of the modeled diffusion can be achieved by using the asymptotic behaviour of the involved densities. This means, by using exponential and power-law functions rather than special functions. with: The nonlocal operators t D β * and x D α 0 are the Caputo fractional time derivative and the Riesz-Feller space derivative, respectively (see [11] for the definition of these operators). This is the same equation discussed in Refs. [11,25], but with a generalized fractional diffusivity A α different from 1.