Paradoxical diffusion: Discriminating between normal and anomalous random walks
Commonly, normal diffusive behavior is characterized by a linear dependence of the second central moment on time, , while anomalous behavior is expected to show a different time dependence, with for subdiffusive and for superdiffusive motions. Here we demonstrate that this kind of qualification, if applied straightforwardly, may be misleading: There are anomalous transport motions revealing perfectly “normal” diffusive character (), yet being non-Markov and non-Gaussian in nature. We use recently developed framework (magdziarz2007b, , Phys. Rev. E 75, 056702 (2007)) of Monte Carlo simulations which incorporates anomalous diffusion statistics in time and space and creates trajectories of such an extended random walk. For special choice of stability indices describing statistics of waiting times and jump lengths, the ensemble analysis of paradoxical diffusion is shown to hide temporal memory effects which can be properly detected only by examination of formal criteria of Markovianity (fulfillment of the Chapman-Kolmogorov equation).
pacs:05.40.Fb, 05.10.Gg, 02.50.-r, 02.50.Ey,
Usually various types of diffusion processes are classified by analysis of the spread of the distance traveled by a random walker. If the mean square displacement grows like with the motion is called subdiffusive, in contrast to normal () or superdiffusive () situations. For a free Brownian particle moving in one dimension, a stochastic random force entering its equation of motion is assumed to be composed of a large number of independent identical pulses. If they posses a finite variance, then by virtue of the standard Central Limit Theorem (CLT) the distribution of their sum follows the Gaussian statistics. However, as it has been proved by Lévy and Khintchine, the CLT can be generalized for independent, identically distributed (i.i.d) variables characterized by non-finite variance or even non-finite mean value. With a Lévy forcing characterized by a stability index independent increments of the particle position sum up yielding bouchaud1990 (), see below. Such enhanced, fast superdiffusive motion is observed in various real situations when a test particle is able to perform unusually large jumps dubkov2008 (); shlesinger1995 (); metzler2004 (). Lévy flights have been documented to describe motion of fluorescent probes in living polymers, tracer particles in rotating flows and cooled atoms in laser fields. They serve also as a paradigm of efficient searching strategies Viswanathan1996 (); sims2008 (); reynolds2009 () in social and environmental problems with some level of controversy Edwards2007 ().
In contrast, transport in porous, fractal-like media or relaxation kinetics in inhomogeneous materials are usually ultraslow, i.e. subdiffusive klafter1987 (); metzler2000 (); metzler2004 (). The most intriguing situations take place however, when both effects – occurrence of long jumps and long waiting times for the next step – are incorporated in the same scenario zumofen1995 (). The approach to this kind of anomalous motion is provided by continuous time random walks (CTRW) which assume that the steps of the walker occur at random times generated by a renewal process. In particular, a mathematical idealization of a free Brownian motion (Wiener process ) can be then derived as a limit (in distribution) of i.i.d random (Gaussian) jumps taken at infinitesimally short time intervals of non-random length. Other generalizations are also possible, e.g. can be defined as a limit of random Gaussian jumps performed at random Poissonian times. The characteristic feature of the Gaussian Wiener process is the continuity of its sample paths. In other words, realizations (trajectories) of the Wiener process are continuous (although nowhere differentiable) doob1942 (). The process is also self-similar (scale invariant) which means that by rescaling and another Wiener process with the same properties is obtained. Among scale invariant stable processes, the Wiener process is the only one which possesses finite variance janicki1994 (); saichev1997 (); doob1942 (); dubkov2005b (). Moreover, since the correlation function of increments depends only on time difference and increments of non-overlapping times are statistically independent, the formal differentiation of yields a white, memoryless Gaussian process vankampen1981 ():
commonly used as a source of idealized environmental noises within the Langevin description
Here stands for the drift term which in the case of a one-dimensional overdamped motion is directly related to the potential , i.e. .
In more general terms the CTRW concept may asymptotically lead to non-Markov, space-time fractional noise , and in effect, to space-time fractional diffusion. For example, let us define where the number of summands is statistically independent from and governed by a renewal process with . Let us assume further that , belong to the domain of attraction of stable distributions, and , whose corresponding characteristic functions , with the density , are given by
for . Here the parameter denotes the stability index, yielding the asymptotic long tail power law for the -distribution, which for is of the type. The parameter () characterizes the scale whereas () defines an asymmetry (skewness) of the distribution.
Note, that for , , the stable variable is defined on positive semi-axis. Within the above formulation the counting process satisfies
where denotes the integer part of the number and stands for the stable distribution of random variable , i.e. . Moreover, since
asymptotically one gets
where and are constants. The resulting (in general non-Markov) process becomes self-similar Lévy random walk tunaley1974 (); saichev1997 (); shlesinger1995 (); kotulski1995 (); nielsen2001 (); barkai2002 (); uchaikin2003 (); srokowski2009 (), i.e.
The asymptotic form given by Eq. (8) can be easily derived scalas2006 (); saichev1997 (); barkai2002 (); germano2009 () for decoupled CTRW by applying Tauberian theorems to the Montroll-Weiss montroll1965 () expression
for the Laplace-Fourier transform of . The latter satisfies the integral (master) equation
Here stands for a waiting time probability density function (PDF) whereas denotes a jump length PDF. With a suitable change of time and space variables, in the limit of , , the Laplace-Fourier transform can be written as
The inverse Laplace transform of can be expressed in a series form saichev1997 ():
where is the Mittag-Leffler function of order
Further application of the inverse Fourier transform in yields saichev1997 () a final series representation of :
where . The above series are divergent for . However, for , the summation has been shown saichev1997 () to produce the closed analytical formula
where (as previously) .
The exact solution of the decoupled CTRW, as given by the infinite sum (7) of stable probability densities, has been studied by Barkai barkai2002 () based on a special choice of PDFs and . In particular, in barkai2002 () the extremely slow convergence of certain CTRW solutions to the (fractional) diffusion approximation has been discussed. The rigorous proof of equivalence between some classes of CTRWs and fractional diffusion has been given by Hilfer and Anton hilfer1995 ()
In this article we investigate CTRW scenarios which, in an asymptotic limit, yield paradoxical diffusion, i.e. the non-Markovian superdiffusive process taking place under sublinear operational time. The combination of long flights and long breaks between them is responsible for the characteristic shape of the PDF and scaling properties of moments. In particular, for , the paradoxical diffusion process exhibits the same scaling as an ordinary Brownian motion despite its PDF is significantly different from Gaussian.
In a forthcoming Section, the relation between fractional calculus and CTRW approach is briefly reminded and the experimental results based on numerical PDF estimators are presented. Section III is devoted to the discussion of scaling properties of moments. In Section IV detection and analysis of memory effects in empirical series of the CTRW-type realizations are proposed and critically tested.
Ii Relation between CTRW and fractional calculus
The theory of stochastic integration of a corresponding Ito-Langevin equation with respect to a general CTRW “measure” has been developed in a series of papers scalas2006 (); meerschaert2006 (); magdziarz2007 (); germano2009 (). Here, we study statistical properties of such a motion constrained to the initial position . To achieve the goals, we adhere to the scheme of stochastic subordination magdziarz2007b (); magdziarz2008 (); magdziarz2007 (), i.e. we obtain the process of primary interest as a function by randomizing the time clock of the process using a different clock . In this approach stands for -stable subordinator, , where denotes a strictly increasing -stable process whose distribution yields a Laplace transform . The parent process is composed of increments of symmetric -stable motion described in an operational time
and in every jump moment the relation is fulfilled. The (inverse-time) subordinator is (in general) non-Markovian hence, as it will be shown, the diffusion process possesses also some degree of memory. The above setup has been recently proved magdziarz2007b (); magdziarz2008 (); magdziarz2007 () to give a proper stochastic realization of the random process described otherwise by a fractional diffusion equation:
with the initial condition . In the above equation denotes the Riemannn-Liouville fractional derivative defined by the relation
and stands for the Riesz fractional derivative with the Fourier transform . Eq. (18) has been otherwise derived from a generalized Master equation metzler1999 (). The formal solution to Eq. (18) can be written metzler1999 () as:
For processes with survival function (cf. Eq. (11)) given by the Mittag-Leffler function Eq. (14), this solution takes an explicit form saichev1997 (); jespersen1999 (); metzler1999 (); scalas2006 (); germano2009 ()
In this paper, instead of investigating properties of an analytical solution to Eq. (18), we switch to a Monte Carlo method magdziarz2007b (); magdziarz2008 (); magdziarz2007 (); gorenflo2002 (); meerschaert2004 () which allows generating trajectories of the subordinated process with the parent process in the potential free case, i.e. for . The assumed algorithm provides means to examine the competition between subdiffusion (controlled by a -parameter) and Lévy flights characterized by a stability index . From the ensemble of simulated trajectories the estimator of the density is reconstructed and statistical qualifiers (like quantiles) are derived and analyzed.
As mentioned, the studied process is self-similar (cf. Eq. (9)). We further focus on examination of a special case for which . As an exemplary values of model parameters we choose (Markovian Brownian diffusion) and (subordination of non-Markovian sub-diffusion with Lévy flights). Additionally we use and as Markovian and non-Markovian counterparts of main cases analyzed. Fig. 1 compares trajectories for all exemplary values of and . Straight horizontal lines (for ) correspond to particle trapping while straight vertical lines (for ) correspond to Lévy flights. The straight lines manifest anomalous character of diffusive process.
To further verify correctness of the implemented version of the subordination algorithm magdziarz2007b (), we have performed extensive numerical tests. In particular, some of the estimated probability densities have been compared with their analytical representations and the perfect match between numerical data and analytical results have been found. Fig. 2 displays numerical estimators of PDFs and analytical results for with (Gaussian case, left top panel), with (Cauchy case, right top panel), with (left bottom panel) and with (right bottom panel). For those last two cases, the expressions for has been derived saichev1997 (), starting from the series representation given by Eq. (15). For the appropriate formula reads
while for the probability density is
and are the integral exponential function and the Airy function respectively. We have also compared results of simulations and Eq. (16) for other sets of parameters , . Also there, the excellent agreement has been detected (results not shown).
Figure 3 and 4 display time-dependent probability densities and corresponding cumulative distribution functions () for “short” and for, approximately, an order of magnitude “longer” times. The persistent cusp sokolov2002 () located at is a finger-print of the initial condition and is typically recorded for subdiffusion induced by the subordinator with . For Markov Lévy-Wiener process dybiec2006 (); denisov2008 () for which the characteristic exponent , the cusp disappears and PDFs of the process become smooth at . In particular, for the Markovian Gaussian case (, ) corresponding to a standard Wiener diffusion, PDF perfectly coincides with the analytical normal density .
The presence of Lévy flights is also well visible in the power-law asymptotic of CDF, see bottom panels of Figs. 3 and 4. Indeed, for independently of the actual value of the subdiffusion parameter and at arbitrary time, for . Furthermore, all PDFs are symmetric with median and modal values located at the origin.
Iii Scaling properties of moments
The self-similar character of the process (cf. Eq. (9)) is an outcome of allowed long flights and long breaks between successive steps. In consequence, the whole distribution scales as a function of with the width of the distribution growing superdiffusively for and subdiffusively for . This scaling is also clearly observable in the behavior of the standard deviation and quantiles , defined via the relation , see Figs. 5, 6. For random walks subject to superdiffusive, long-ranging trajectories (), the asymptotic scaling is observed for sufficiently long times, cf. Fig. 5. On the other hand, normal (Gaussian) distribution of jumps superimposed on subdiffusive motion of trapped particles () clearly shows rapid convergence to the law. Notably, both sets and lead to the same scaling , although in the case , the process is non-Markov, in contrast to a standard Gaussian diffusion obtained for . Thus, the competition between subdiffusion and Lévy flights questions standard ways of discrimination between normal (Markov, ) and anomalous (generally, non-Markov ) diffusion processes.
Indeed, for , the process is not only self-similar but it is also memoryless (i.e. Markovian). In such a case, the asymptotic PDF is -stable janicki1994 (); weron1995 (); weron1996 (); dybiec2006 () with the scale parameter growing with time like , cf. Eq. (3). This is no longer true for subordination with when the underlying process becomes non-Markovian and the spread of the distribution follows the -scaling (cf. Fig. 6, right panels).
Some additional care should be taken when discussing the scaling character of moments of bouchaud1990 (); fogedby1998 (); jespersen1999 (). Clearly, Lévy distribution (with ) of jump lengths leads to infinite second moment (see Eqs. (3) and (4))
irrespectively of time . Moreover, the mean value of stable variables is finite for only ( for symmetric case under investigation). Those observations seem to contradict demonstration of the scaling visible in Fig. 5 where standard deviations derived from ensembles of trajectories are finite and grow in time according to a power law. A nice explanation of this behavior can be given following argumentation of Bouchaud and Georges bouchaud1990 (): Every finite but otherwise arbitrarily long trajectory of a Lévy flight, i.e. the stochastic process underlying Eq. (18) with , is a sum of finite number of independent stable random variables. Among all summed stable random variables there is the largest one, let say . The asymptotic form of a stable densities
together with the estimate for allow one to estimate how standard deviations grows with a number of jumps . In fact, the largest value can be estimated from the condition
which locates most of the “probability mass” in events not exceeding the step length (otherwise, the relation states that occurred at most once in trials, bouchaud1990 ()). Alternatively, can be estimated as a value which maximizes probability that the largest number chosen in trials is
Due to finite, but otherwise arbitrarily large number of trials , the effective distributions becomes restricted to the finite domain which size is controlled by Eq. (28). Using the estimated threshold, see Eq. (28) and asymptotic form of stable densities, see Eq. (25), it is possible to derive an estimate of
Finally, after jumps
Consequently, for Lévy flights standard deviation grows like a power law with the number of jumps . In our CTRW scenario incorporating competition between long rests and long jumps, the number of jumps grows sublinearly in time, , leading effectively to with and . Since in any experimental realization tails of the Lévy distributions are almost inaccessible and therefore effectively truncated, analyzed sample trajectories follow the pattern of the scaling, which is well documented in Fig. 5.
Iv Discriminating memory effects
Clearly, by construction, for the limiting “anomalous diffusion” process is non-Markov. This feature is however non-transparent when discussing statistical properties of the process by analyzing ensemble-averaged mean-square displacement for the parameters set , (e.g. ) when contrary to what might be expected, , similarly to the standard, Markov, Gaussian case. This observation implies a different problem to be brought about: Given an experimental set of data, say time series representative for a given process, how can one interpret its statistical properties and conclude about anomalous (subdiffusive) character of underlying kinetics? The similar question has been carried out in a series of papers discussing use of transport coefficients in systems exhibiting weak ergodicity breaking (see he2008 () and references therein).
To further elucidate the nature of simulated data sets for , we have adhered to and tested formal criteria vankampen1981 (); fulinski1998 (); dybiec2009b (); dybiec-sparre () defining the Markov process. The standard formalism of space- and time- continuous Markov processes requires fulfillment of the Chapman-Kolmogorov equation ()
along with the constraint for conditional probabilities which for the “memoryless” process should not depend on its history. In particular, for a hierarchy of times , the following relation has to be satisfied
Eqs. (31) and (32) have been used to directly verify whether the process under consideration is of the Markovian or non-Markovian type. From Eq. (31) squared cumulative deviation between LHS and RHS of the Chapman-Kolmogorov relation summed over initial () and final () states has been calculated fulinski1998 ()
The same procedure can be applied to Eq. (32) leading to
Figure 7 presents evaluation of (top panel) and (bottom panel) for and as a function of the intermediate time . It is seen that deviations from the Chapman-Kolmogorov identity are well registered for processes with long rests when subdiffusion wins competition with Lévy flights at the level of sample paths. The tests based on (see Eq. (33)) and (see Eq. (34)) have comparative character. The deviations and are about three order of magnitudes higher for the parameter sets and than and values for the Markovian counterparts with and , , respectively. Performed analysis clearly demonstrates non-Markovian character of the limiting diffusion process for and the findings indicate that scaling of PDF, or, in consequence, scaling of the variance and interquantile distances (see Fig. 6) do not discriminate satisfactory between “normal” and “anomalous” diffusive motions zumofen1995 (). In fact, linear in time spread of the second moment does not necessarily characterize normal diffusion process. Instead, it can be an outcome of a special interplay between subdiffusion and Lévy flights combined in the subordination . The competition between both processes is better displayed in analyzed sample trajectories where combination of long jumps and long trapping times can be detected, see Fig. 1.
Summarizing, by implementing Monte Carlo simulations which allow visualization of stochastic trajectories subjected to subdiffusion (via time-subordination) and superdiffusive Lévy flights (via extremely long jumps in space), we have demonstrated that the standard measure used to discriminate between anomalous and normal behavior cannot be applied straightforwardly. The mean square displacement alone, as derived from the (finite) set of time-series data does not provide full information about the underlying system dynamics. In order to get proper insight into the character of the motion, it is necessary to perform analysis of individual trajectories. Subordination which describes a transformation between a physical time and an operational time of the system sokolov2000 (); magdziarz2007b () is responsible for unusual statistical properties of waiting times between subsequent steps of the motion. In turn, Lévy flights are registered in instantaneous long jumps performed by a walker. Super- or sub- linear character of the motion in physical time is dictated by a coarse-graining procedure, in which fractional time derivative with the index combines with a fractional spatial derivative with the index . Such situations may occur in motion on random potential surfaces where the presence of vacancies and other defects introduces both – spatial and temporal disorder sancho2004 (). We believe that the issue of the interplay of super- and sub-diffusion with a crossover resulting in a pseudo-normal paradoxical diffusion may be of special interest in the context of e.g. facilitated target location of proteins on folding heteropolymers belik2007 () or in analysis of single particle tracking experiments lubelski2008 (); he2008 (); metzler2009 (); bronstein2009 (), where the hidden subdiffusion process can be masked and appear as a normal diffusion.
Acknowledgements.This project has been supported by the Marie Curie TOK COCOS grant (6th EU Framework Program under Contract No. MTKD-CT-2004-517186) and (BD) by the Foundation for Polish Science. The authors acknowledge many fruitfull discussions with Andrzej Fuliński, Marcin Magdziarz and Aleksander Weron.
- (1) M. Magdziarz and A. Weron, Phys. Rev. E 75, 056702 (2007).
- (2) J. P. Bouchaud and A. Georges, Phys. Rep. 195, 127 (1990).
- (3) A. A. Dubkov, B. Spagnolo, and V. V. Uchaikin, Int. J. Bifurcation Chaos. Appl. Sci. Eng. 18, 2649 (2008).
- (4) Lévy Flights and Related Topics in Physics, edited by M. F. Shlesinger, G. M. Zaslavsky, and J. Frisch (Springer Verlag, Berlin, 1995).
- (5) R. Metzler and J. Klafter, J. Phys. A: Math. Gen. 37, R161 (2004).
- (6) G. M. Viswanathan et al., Nature (London) 381, 413 (1996).
- (7) D. W. Sims et al., Nature (London) 451, 1098 (2008).
- (8) A. M. Reynolds and C. J. Rhodes, Ecology 90, 877 (2009).
- (9) A. M. Edwards et al., Nature (London) 449, 1044 (2007).
- (10) J. Klafter, A. Blumen, and M. F. Shlesinger, Phys. Rev. A 35, 3081 (1987).
- (11) R. Metzler and J. Klafter, Phys. Rep. 339, 1 (2000).
- (12) G. Zumofen and J. Klafter, Phys. Rev. E 51, 1818 (1995).
- (13) J. L. Doob, Annals of Mathematics 43, 351 (1942).
- (14) A. Janicki and A. Weron, Simulation and Chaotic Behavior of -Stable Stochastic Processes (Marcel Dekker, New York, 1994).
- (15) A. I. Saichev and G. M. Zaslavsky, Chaos 7, 753 (1997).
- (16) A. A. Dubkov and B. Spagnolo, Fluct. Noise Lett. 5, L267 (2005).
- (17) N. G. van Kampen, Stochastic Processes in Physics and Chemistry (North–Holland, Amsterdam, 1981).
- (18) J. Tunaley, J. Stat. Phys. 11, 397 (1974).
- (19) M. Kotulski, J. Stat. Phys. 81, 777 (1995).
- (20) Lévy Processes: Theory and Applications, edited by O. E. Barndorff-Nielsen, T. Mikosch, and S. I. Resnick (Birkhäuser, Boston, 2001).
- (21) E. Barkai, Chem. Phys. 284, 13 (2002).
- (22) V. V. Uchaikin, Physics Uspekhi 46, 821 (2003).
- (23) T. Srokowski, Physica A 388, 1057 (2009).
- (24) E. Scalas, Physica A 362, 225 (2006).
- (25) G. Germano, M. Politi, E. Scalas, and R. L. Schilling, Phys. Rev. E. 79, 066102 (2009).
- (26) E. W. Montroll and G. H. Weiss, J. Math. Phys. 6, 167 (1965).
- (27) R. Hilfer and L. Anton, Phys. Rev. E 51, R848 (1995).
- (28) M. M. Meerschaert and E. Scalas, Physica A 370, 114 (2006).
- (29) M. Magdziarz, A. Weron, and K. Weron, Phys. Rev. E 75, 016708 (2007).
- (30) M. Magdziarz, A. Weron, and J. Klafter, Phys. Rev. Lett. 101, 210601 (2008).
- (31) R. Metzler, E. Barkai, and J. Klafter, Europhys. Lett. 46, 431 (1999).
- (32) S. Jespersen, R. Metzler, and H. C. Fogedby, Phys. Rev. E 59, 2736 (1999).
- (33) E. Scalas, R. Gorenflo, and F. Mainardi, Phys. Rev. E 69, 011107 (2004).
- (34) R. Gorenflo et al., Chem. Phys. 284, 521 (2002).
- (35) M. M. Meerschaert and C. Tadjeran, J. Comput. Appl. Math. 172, 65 (2004).
- (36) I. M. Sokolov, J. Klafter, and A. Blumen, Phys. Today 55, 48 (2002).
- (37) B. Dybiec, E. Gudowska-Nowak, and P. Hänggi, Phys. Rev. E 73, 046104 (2006).
- (38) S. I. Denisov, W. Horsthemke, and P. Hänggi, Phys. Rev. E 77, 061112 (2008).
- (39) A. Weron and R. Weron, Lect. Not. Phys. 457, 379 (1995).
- (40) R. Weron, Statist. Probab. Lett. 28, 165 (1996).
- (41) H. C. Fogedby, Phys. Rev. E 58, 1690 (1998).
- (42) Y. He, S. Burov, R. Metzler, and E. Barkai, Phys. Rev. Lett. 101, 058101 (2008).
- (43) A. Fuliński et al., Phys. Rev. E 58, 919 (1998).
- (44) B. Dybiec, J. Stat. Mech. P08025 (2009).
- (45) B. Dybiec and E. Gudowska-Nowak, Europhys. Lett. 88, 10003 (2009).
- (46) I. M. Sokolov, Phys. Rev. E 63, 011104 (2000).
- (47) J. M. Sancho et al., Phys. Rev. Lett. 92, 250601 (2004).
- (48) V. V. Belik and D. Brockmann, New J. Phys. 9, 54 (2007).
- (49) A. Lubelski and J. Klafter, J. Phys. Chem. B 112, 12740 (2008).
- (50) R. Metzler et al., Acta Phys. Pol. B 40, 1315 (2009).
- (51) I. Bronstein et al., Phys. Rev. Lett. 103, 018102 (2009).