Stochastic solution of space-time fractional diffusion equations

PHYSICAL REVIEW E, VOLUME 65, 041103 Stochastic solution of space-time fractional diffusion equations Mark M. Meerschaert* Department of Mathematics, University of Nevada, Reno, Nevada 89557-0084 David A. Benson Desert Research Institute, 2215 Raggio Parkway, Reno, Nevada 89506-0220 Hans-Peter Scheffler Fachbereich Mathematik, University of Dortmund, 44221 Dortmund, Germany Boris Baeumer Department of Mathematics and Statistics, University of Otago, Dunedin, New Zealand ~Received 16 July 2001; published 28 March 2002! Classical and anomalous diffusion equations employ integer derivatives, fractional derivatives, and other pseudodifferential operators in space. In this paper we show that replacing the integer time derivative by a fractional derivative subordinates the original stochastic solution to an inverse stable subordinator process whose probability distributions are Mittag-Leffler type. This leads to explicit solutions for space-time fractional diffusion equations with multiscaling space-fractional derivatives, and additional insight into the meaning of these equations. DOI: 10.1103/PhysRevE.65.041103 PACS number~s!: 05.40.Fb, 02.50.2r, 05.10.Gg, 47.55.Mh I. INTRODUCTION Space-fractional diffusion equations @1– 4# have been useful as models of anomalous transport in many diverse disciplines, including finance, semiconductor research, biology, and hydrogeology @5–7#. In the context of the flow in porous media, fractional space derivatives model large motions through highly conductive layers or fractures, while fractional time derivatives describe particles that remain motionless for extended periods of time. Dissolved solutes may sorb to solid material @8# or diffuse into immobile-water zones of various sizes @9#. The scalar space-fractional diffusion equation governs Lévy motion, and the tail parameter a of the Lévy motion equals the order of the fractional derivative. Solutions to the vector space-fractional diffusion equation are operator Lévy motions @10# that may scale at different rates in different directions. The matrix exponent of the fractional derivative is related to the scaling rates in a similar manner @11,12#. A more general diffusion equation governs any Lévy process X(t) @13,14#. The probability density p(x,t) of any such process solves a diffusion-type equation ] p ~ x,t ! 5L p ~ x,t ! ; ]t p ~ x,0 ! 5 d ~ x! , ~1! where L is the generator of the Feller semigroup S t f (x) 5 * f (x2y) p(y,t)dy @15–17#. In this case, we say that X(t) is the stochastic solution to Eq. ~1!. The generator L f (x) 5limt↓0 t 21 @ S t f (x)2 f (x) # . If X(t) is an a -stable Lévy motion without drift, then L is a fractional derivative operator of order a . *Electronic address: mcubed@unr.edu 1063-651X/2002/65~4!/041103~4!/$20.00 Fractional time derivatives are important in reactive transport, since solutes may interact with the immobile porous medium in highly nonlinear ways. There is evidence that solutes may sorb for random amounts of time that have a power law distribution @8#, or move into irregularly sized blocks of relatively immobile water, producing similar behavior @9#. If the first moment of these time delays diverges, then a fractional time derivative applies @6#. The fractional time derivative ] g g(t)/ ] t g for 0, g ,1 is the inverse Laplace transform of s g g(s), where g(s)5L@ g(t) # is the usual Laplace transform. In this paper, we find the stochastic solution to the space-time fractional diffusion equation ] g q ~ x,t ! ]tg 5Lq ~ x,t ! 1 d ~ x! t 2g . G ~ 12 g ! ~2! We show that if X(t) is the stochastic solution to Eq. ~1! then X(V t ) is the corresponding solution to Eq. ~2!, where V t is the inverse Lévy process @18# for the stable subordinator with index g . The fractional time derivative subordinates X(t) to the inverse stable subordinator V t . The space-time fractional diffusion equation is also connected with scaling limits of continuous time random walks ~CTRW, see @6#!. The spatial operator L depends on the jump size distribution @11,12#. A fractional time derivative of order 0, g ,1 pertains when the random waiting time T between jumps satisfies P(T.t)'t 2 g so that E(T)5`. The infinite mean waiting time CTRW limit is the finite mean waiting time CTRW limit, subordinated to the inverse stable subordinator V t . The random variable V t has a Mittag-Leffler distribution @19# previously noted in connection with fractional time derivatives @4,20# and relaxation @21#. 65 041103-1 ©2002 The American Physical Society MEERSCHAERT, BENSON, SCHEFFLER, AND BAEUMER PHYSICAL REVIEW E 65 041103 c 2 g N [ct] ⇒V t II. CTRW SCALING LIMITS CTRW were introduced @22,23# to study random walks on a lattice. They are now used in physics to model a wide variety of phenomena connected with anomalous diffusion @23–25#. With finite mean waiting times, the jump process is asymptotically linear, and the CTRW behaves in a manner similar to the original random walk for large time @20,26#. For a scalar process, finite variance jumps lead to Brownian motion in the scaling limit. Infinite variance jumps with power law tails lead to Lévy motion. Vector jumps with finite second moments lead to multivariable Brownian motion. Vector jumps with power law tails lead to multivariable Lévy motion, or operator Lévy motion if the power law behavior varies with the direction of motion @11,12#. The speed of convergence to the CTRW scaling limit, and the implications for fractional diffusion modeling, are discussed in a recent paper of Barkai @27#. Many physical applications involve infinite mean waiting times @21,28#. Introducing infinite mean waiting times has the effect of subordinating the CTRW scaling limit to the inverse process of a stable subordinator whose index g is the same as the power law tail index of the waiting times. Essentially, this is because the counting process for particle jumps is inverse to the jump time process. The jump time process is asymptotically the stable subordinator, so the counting process for particle jumps is asymptotically the inverse stable subordinator. A rigorous mathematical proof appears in @29#. We recount the basic ideas here to emphasize the physical applications. Given iid positive random variables J i let T n n 5 ( i51 J i denote the time of the nth particle jump. The pon sition of the particle after the nth jump is W(n)5 ( i51 Yi where Yi are iid and assumed independent of J i . Then N t 5max$n:Tn<t% counts the number of particle jumps by time t.0 and the CTRW variable W(N t ) gives the position of the particle at time t.0. If Y has zero mean and finite second moments, the simple random walk of particle jumps c 21/2W~@ ct # ! ⇒X~ t ! as c→`, ~3! where the scaling limit X(t) is a Brownian motion. Shrinking the spatial coordinates by c 1/2 compensates expanding the time scale by c according to the central limit theorem. If P(J.t)'t 2 g for some 0, g ,1 then c 21/g T [ct] ⇒B t as c→` ~4! according to the extended central limit theorem @15# where the scaling limit B t is the stable subordinator process @13#. The g -stable random variable B t is totally positively skewed, hence this Lévy process is strictly increasing. The inverse process V t 5inf$ t:B t . t % is also called the hitting time or first passage time process. Using the fact that T n ,N t are inverse, so that $ N t >x % 5 $ T ⌈x⌉ <t % , along with Eq. ~4! yields as c→`. Hence N [ct] 'c g V t , and together with Eq. ~3! this yields c 2 g /2W~ N [ct] ! ' ~ c g ! 21/2W~@ c g V t # ! ⇒X~ V t ! as c→`, so that the Brownian motion X(t) is subordinated to the inverse stable subordinator V t . The inverse processes have inverse distributional scaling B ct 5c 1/g B t and V ct 5c g V t , and together with the classical scaling for Brownian motion X(ct)5c 1/2X(t) this shows that the CTRW limit is subdiffusive X~ V ct ! 5X~ c g V t ! 5c g /2X~ V t ! with Hurst index H5 g /2,1/2. Since P(V t <t)5 P(B t > t ) 5 P(t 1/g B 1 > t )5 P„(B 1 / t ) 2 g <t… the random variable V t has the same density function as ( t /B 1 ) g . The density g g of the stable random variable B 1 has Laplace transform g L@ g g (t) # 5e 2s . Computing moments of (t/B 1 ) g shows that V t has a Mittag-Leffler distribution @19#. If p(x,t) is the density of X(t) then a conditioning argument along with a simple change of variable shows that X(V t ) has density q ~ x,t ! 5 E 5 t g ` 0 p„x, ~ t/s ! g …g g ~ s ! ds E ` 0 p ~ x,u ! g g ~ tu 21/g ! u 21/g 21 du. ~5! Analytical estimates in @29# show that q(k,t)>C i ki 2b for large i ki , so X(V t ) does not have a normal density and hence cannot be a fractional Brownian motion @30#. If P( i Yi .r)'r 2 a for some 0, a ,2 then X(t) is an a -stable Lévy motion and the CTRW limit X(V t ) has Hurst index H5 g / a . If the tail index varies with the spatial coordinate, operator norming applies @12#. Then X(ct)5c E X(t) leads to X(V ct )5c g E X(t) so that the Hurst index H5 g E is a matrix. For a diagonal exponent E5diag(1/a 1 , . . . ,1/a d ) the ith coordinate X i (t) is an a i -stable Lévy motion and X i (V t ) is self-similar with Hurst index g / a i . Diagonalizable matrix exponents introduce a change of coordinates. Repeated eigenvalues thicken probability tails by a logarithmic factor, and complex exponents introduce rotations, leading to discrete scale invariance @31#. In every case, the scaling limit X(t) of the simple random walk is subordinated by the inverse stable subordinator V t and the density changes from p(x,t) to q(x,t) via Eq. ~5! when infinite mean waiting times are introduced. Next, we show that this change corresponds to a fractional time derivative in the diffusion equation. III. TIME-FRACTIONAL DIFFUSION EQUATIONS Wyss @32# and Schneider and Wyss @33# studied a timefractional diffusion equation. Zaslavsky @34# introduced the space-time fractional kinetic equation ~2! for Hamiltonian chaos. When L52 v ] / ] x1D ] a / ] u x u a Saichev and Zaslavsky @35# show that if p(x,t) solves Eq. ~1! then the function q(x,t) given by Eq. ~5! solves Eq. ~2!. When a 52 they call the stochastic solution to Eq. ~2! a ‘‘fractal 041103-2 STOCHASTIC SOLUTION OF SPACE-TIME . . . PHYSICAL REVIEW E 65 041103 Brownian motion.’’ We prefer the term ‘‘time-fractional diffusion’’ to avoid confusing this process with the well-known fractional Brownian motion. In fact, if X(t) is a Brownian motion, the stochastic solution to Eq. ~1! in this case, then the stochastic solution to Eq. ~2! is X(V t ) where V t is the inverse g -stable subordinator. This interesting stochastic process is self-similar with Hurst index g /2 so it is subdiffusive. However X(V t ) does not have a Gaussian distribution and it does not have stationary increments @29#, so it is not fractional Brownian motion. Barkai, Metzler, and Klafter @20# introduce a fractional Fokker-Planck equation equivalent to Eq. ~2! with tional Cauchy problem ~2! via the inverse Lévy transform ~5!. We summarize the essentials here in order to clarify the argument. Use s g 21 5L@ t 2 g /G(12 g ) # and take Laplace-Fourier transforms (x°k,t°s) in Eq. ~2! to get s g q(k,s) 5 c (k)q(k,s)1s g 21 , where c (k) is the Fourier symbol of L, so that F @ L f (x) # 5 c (k) f (k). Then q ~ k,s ! 5 5 2 L52 ] V 8~ x ! ] 1K 1 2 . ]x mh1 ]x Barkai @4# applies Eq. ~5!, which he calls the inverse Lévy transform of p(x,t), to the solution of Eq. ~1! in order to solve this fractional Fokker-Planck equation. Scalar solutions to Eq. ~1! with L52 v S a a ] 12 b ] 11 b ] 1D 1 a ]x 2 ] ~ 2x ! 2 ]xa s g 21 s g 2 c ~ k! E ` are multivariable stable densities @11#, where erator with Fourier symbol E i u i 51 exp$ 2 @ s g 2 c ~ k!# u % du ~6! 0 g e 2s u 5e 2(su E ` 0 1/g ) g 5 E ` e 2su 1/g v 0 g g~ v ! d v e 2st g g ~ u 21/g t ! u 21/g dt. ~7! Then compute g s g 21 e 2s u 5 5 is the op- 21 d g u ds 1 gu E ` 0 SE ` 0 e 2st g g ~ u 21/g t ! u 21/g dt D te 2st g g ~ u 21/g t ! u 21/g dt ~8! and combine with Eq. ~6! to write q(k,s) as E Sg E E SE ` ~ ik• u ! a m ~ u ! d u . 0 1 u ` 0 ` e 2st 0 0 D t p ~ k,u ! g g ~ u 21/g t ! u 21/g 21 du dt. g Now invert the Laplace transform to obtain q ~ k,t ! 5 where the generalized fractional derivative D te 2st g g ~ u 21/g t ! u 21/g dt p ~ k,u ! du ` 5 L52 v •“1 12 “•A“1F, E 0 g 5 If a 52, this integral reduces to (ik)A(ik) 8 where the matrix A has i j component * u i u j m( u )d u , and solutions are vector Brownian motion. Operator stable densities, where the stable index depends on the coordinate, solve Eq. ~1! with Ff ~ x ! 5 ` s g 21 e 2s u p ~ k,u ! du, a L52 v •“1D“ m a “m E using * `0 e 2au du5a 21 and p(k,t)5e c (k)t , which follows g g from Eq. ~1!. Use d(e 2s u )/ds52 g s g 21 ue 2s u to get g g g s g 21 e 2s u 52( g u) 21 d(e 2s u )/ds. Recall that e 2s 5L@ g g (t) # and write D are a -stable densities @36,37#, purely symmetric when the skewness b 50 @2,38# and maximally skewed when b 51 @3,34#. When a 52 the skewness b is irrelevant, and the solutions are normal densities. Vector solutions for 5s g 21 t g E ` 0 p ~ k,u ! g g ~ tu 21/g ! u 21/g 21 du and invert the Fourier transform to get Eq. ~5!. IV. CONCLUSIONS @ f ~ x2y ! 2 f ~ x ! 1y“ f ~ x !# d f ~ y ! and d f (r E u )5r 22 dr m( u )d u is an operator stable Lévy measure @10,12,39#. These are all abstract Cauchy problems @16,17# whose solution p(x,t) is the family of densities for a Lévy process, a stationary independent increment process that includes Brownian motion and ~operator! Lévy motion as special cases. Baeumer and Meerschaert @40# give a rigorous mathematical proof that any solution to the abstract Cauchy problem ~1! is transformed to a solution of the frac- Infinite mean waiting times subordinate CTRW scaling limits to an inverse stable subordinator, equivalent to applying an inverse Lévy transform ~5! to the solution density. Since the solutions to time-fractional diffusion equations are also obtained via the inverse Lévy transform, a g -fractional time derivative in a diffusion equation has the effect of subordinating the stochastic solution to the inverse process of a g -stable subordinator. When applied to the classical diffusion equation, this procedure produces the ‘‘fractal Brownian mo- 041103-3 MEERSCHAERT, BENSON, SCHEFFLER, AND BAEUMER PHYSICAL REVIEW E 65 041103 tion’’ of Saichev and Zaslavsky as the solution to the timefractional diffusion equation, a model for subdiffusion. This interesting stochastic process is not the same as fractional Brownian motion, but rather a completely different stochastic process. For CTRW models with coupled memory, in which particle jumps Yi and waiting times J i are dependent, the effect of infinite mean waiting times is more complicated @26#. In a forthcoming paper we will show that infinite mean waiting times in a coupled CTRW model also induce subor- @1# @2# @3# @4# @5# @6# @7# @8# @9# @10# @11# @12# @13# @14# @15# @16# @17# S. Bochner, Proc. Natl. Acad. Sci. U.S.A. 35, 368 ~1949!. A. Compte, Phys. Rev. E 53, 4191 ~1996!. A. Compte, Phys. Rev. E 55, 6821 ~1997!. E. Barkai, Phys. Rev. E 63, 046118 ~2001!. R. Hilfer, Applications of Fractional Calculus in Physics ~World Scientific, Singapore, 2000!. An excellent review is given by R. Metzler and J. Klafter, Phys. Rep. 339, 1 ~2000!. R. Lenormand and B. Wang, Transp. Porous Media 18, 263 ~1995!; D. Benson, S. Wheatcraft, and M. Meerschaert, Water Resour. Res. 36, 1403 ~2000!; D. Benson, et al., Transp. Porous Media 42, 211 ~2001!; B. Baeumer, et al., Water Resour. Res. 37, 1543 ~2001!. Y. Hatano and N. Hatano, Water Resour. Res. 34, 1027 ~1998! R. Haggerty, S. McKenna, and L. Meigs, Water Resour. Res. 36, 3467 ~2000!. M. Meerschaert and H.P. Scheffler, Limit Theorems for Sums of Independent Random Vectors ~Wiley, New York, 2001!. M.M. Meerschaert, D.A. Benson, and B. Baumer, Phys. Rev. E 59, 5026 ~1999!. M.M. Meerschaert, D.A. Benson, and B. Baeumer, Phys. Rev. E 63, 021112 ~2001!. J. Bertoin, Lévy Processes ~Cambridge University Press, Cambridge, 1996!. K.I. Sato, Lévy Processes and Infinitely Divisible Distributions ~Cambridge University Press, Cambridge, 1999!. W. Feller, An Introduction to Probability Theory and Its Applications, 2nd ed. ~Wiley, New York, 1971!, Vol. II. A. Pazy, Semigroups of Linear Operators and Applications to Partial Differential Equations, Applied Mathematical Sciences Vol. 44 ~Springer, New York, 1983!. W. Arendt et al., Vector-valued Laplace Transforms and Cauchy Problems, Monographs in Mathematics ~BirkhaeuserVerlag, Berlin, 2001!. dination by an inverse stable subordinator, but in that case the two processes are dependent. ACKNOWLEDGMENTS M.M.M. and B.B. were partially supported by NSF DES Grant No. 9980484, and D.A.B. was partially supported by NSF-DES Grant No. 9980489 and DOE-BES Grant No. DEFG03-98ER14885. @18# N. Bingham, Z. Wahrscheinlichkeitstheor. Verwandte. Geb. 17, 1 ~1971!. @19# L. Bondesson, G. Kristiansen, and F. Steutel, Statist. Probab. Lett. 28, 271 ~1996!. @20# E. Barkai, R. Metzler, and J. Klafter, Phys. Rev. E 61, 132 ~2000!. @21# M. Kotulski and K. Weron, Lecture Notes in Statistics ~Springer, New York, 1995!, Vol. 114, p. 379. @22# E.W. Montroll and G.H. Weiss, J. Math. Phys. 6, 167 ~1965!. @23# H. Scher and M. Lax, Phys. Rev. B 7, 4491 ~1973!. @24# J. Klafter, A. Blumen, and M.F. Shlesinger, Phys. Rev. A 35, 3081 ~1987!. @25# V.V. Uchaikin and V.M. Zolotarev, Chance and Stability. Stable Distributions and Their Applications ~VSP, Utrecht, 1999!. @26# M. Kotulski, J. Stat. Phys. 81, 777 ~1995!. @27# E. Barkai, e-print cond-mat/0108024. @28# M. Shlesinger, J. Klafter, and Y.M. Wong, J. Stat. Phys. 27, 499 ~1982!. @29# M. Meerschaert and H. P. Scheffler ~unpublished!. @30# B. Mandelbrot and J.W. Van Ness, SIAM Rev. 10, 422 ~1968!. @31# D. Sornette, Phys. Rep. 297, 5 ~1998!. @32# W. Wyss, J. Stat. Phys. 27, 2782 ~1986!. @33# W.R. Schneider and W. Wyss, J. Math. Phys. 30, 134 ~1989!. @34# G. Zaslavsky, Physica D 76, 110 ~1994!. @35# A. Saichev and G. Zaslavsky, Chaos 7, 753 ~1997!. @36# A. Chaves, Phys. Lett. A 239, 13 ~1998!. @37# D. Benson, S. Wheatcraft, and M. Meerschaert, Water Resour. Res. 36, 1413 ~2000!. @38# R. Metzler, J. Klafter, and I.M. Sokolov, Phys. Rev. E 58, 1621 ~1998!. @39# Z. Jurek and J.D. Mason, Operator-Limit Distributions in Probability Theory ~Wiley, New York, 1993!. @40# B. Baeumer and M. Meerschaert, Fractional Calculus and Appl. Analysis 4, 481 ~2001!. 041103-4