Transforming Gaussian diffusion into fractional a generalized law of large numbers approach

14 Pages·0.28 MB·English

by  E. Barkai

Transforming Gaussian diffusion into fractional, a generalized law of large numbers
approach
E. Barkai
Department of Chemistry and Center for Materials Science and Engineering. Massachusetts Institute of Technology. 77
Massachusetts Avenue Cambridge, MA 02139.
0
(February 1, 2008)
0
0
2
tion. TheFFPEisanasymptoticequationwhichextends
n
a theCTRWtoincludetheeffectofanexternalforcefield.
The fractional Fokker–Planck equation (FFPE) [R. Met-
J The CTRW itself was not designedfor a description of a
zler, E. Barkai, J. Klafter Phys. Rev. Lett. 82, 3563 (1999)]
7 particle in external force field thus the FFPE describes
describes an anomalous sub diffusive behavior of a particle
new type of behavior not explored in depth yet.
in an external force field. In this paper we present the solu-
]
h tionoftheFFPEintermsofanintegraltransformation. The The main result in this paper concerns a transforma-
c transformation maps the solution of ordinary Fokker–Planck tion of ordinary Gaussian diffusion into fractional diffu-
e equation onto the solution of the FFPE. We investigate in sion [5,31–33]. Let Pα(r,t) be the propagator (Green’s
m
detail the force free particle and the particle in uniform and function)offractionalFokker–Planckequation(aspecial
- harmonic fields. The meaning of the transformation is ex- case is the fractional diffusion equation), here 0<α<1
t
a plained based on the asymptotic solution of the continuous isthe fractionalexponentandα=1isthe standardcase
st timerandomwalk(CTRW).Wealsofindanexactsolutionof describedbyordinaryFokker-Planckequation. Thenthe
t. the CTRW and compare the CTRW result with the integral solution, Pα(r,t), is found based on the transformation
a solution of theFFPE for theforce free case.
m ∞ dN(s,t)
02.50.-a,r05.40.Fb,05.30.Pr Pα(r,t)= P1(r,s)ds, (1)
- ds
d Z0 (cid:20) (cid:21)
n and
o
t
c I. INTRODUCTION
N(s,t)=1 L (2)
[ − α s1/α
(cid:18) (cid:19)

1 The continuous time random walk (CTRW), intro- denotes the “inverse” one sided L´evy stable distribution
v duced by Montroll and Weiss [1,2], describes different (i.e., L (x) is the one sided L´evy stable distribution).
1 α
typesofdiffusionprocessesincludingthestandardGaus- For example considerthe forcefree case. The solutionof
8
0 sian diffusion, sub diffusion [2–5] and L´evy walks [6–12]. the FFPE Pα(r,t), with initial conditions concentrated
1 The CTRW has been a useful toolin many diversefields on the origin, is found by transforming P (r,s) and the
1
0 and over the last three decades e.g., it was used to de- transformed function is the well known Gaussian solu-
0 scribe transportin disorder medium [4,13,14]and in low tion of the integer diffusion equation with initial con-
0
dimensional chaotic systems [8,15–17]. More recently ditions concentrated on the origin. The transformation
/
t fractionalkineticequationswereinvestigatedasatoolde- Eq. (1)besidesits practicalvalue forfindingsolutionsof
a
m scribingphenomenologicallyanomalousdiffusion[18–24]. the FFPE also explains its meaning and relation to the
In these equations fractional derivatives [25] replace or- CTRW (see details below).
-
d dinary integer derivatives in the standardinteger kinetic We list previous work on the integral transformation
n equation. It is believed that fractional calculus can be Eq. (1)
o
used to model non integer diffusion phenomena. Schnei- 1BouchaudandGeorges[5]fortheasymptoticlongtime
c
: der and Wyss have formulated the fractional diffusion limit of the CTRW, then P1(r,s) is Gaussian,
v equation (see Sec. III) [26] and under certain condi- 2 Zumofen and Klafter [31]for CTRW on a lattice, then
i
X tions [27–29] the fractional diffusion equation describes P (r,s) is solution of an ordinary random walk on a lat-
1
r theasymptoticlargetimebehaviorofthedecoupledsub- tice,
a diffusive CTRW (with r2 tδ; δ < 1). Thus some 3SaichevandZaslavsky[32]foronedimensionalsolution
h i ∼
whatlikeordinaryrandomwalkswhichcanbe described offractionaldiffusionequation,theyhavealsoconsidered
asymptoticlybythediffusionequation,socantheCTRW an extension which includes L´evy flights,
be described by fractional diffusion equation. 4 Barkai and Silbey [33] for the fractional Ornstein-
Recently, Metzler, Barkai and Klafter [30] have inves- Uhlenbeck process.
tigatedfractionalFokker–Planckequation(FFPE)tode- We generalize these results for the dynamics described
scribe an anomalous sub-diffusion motion in an external bythe FFPE,andinvestigateingreaterdetailEq. (1)in
non–linearforcefield. Intheabsenceoftheexternalforce the context of CTRW and fractional diffusion equation
field the FFPE reduces to the fractional diffusion equa- in dimensions d=1,2,3.
1
This paper is organized as follows. In Sec. (II) we is the Laplace transform of P(r,t).
follow [5] and derive the transformation Eq. (1) based W (r,s) willgenerallydepend onf(r), howeverwe are
on the CTRW. We find an exact solution to the CTRW interestedonlyinthelargetimebehaviorofP(r,t)mean-
valid for short and long times and compare this solution ing that only contributions from large s are important.
with the result obtained from Eq. (1). In Sec. (III) we From standard central limit theorem we know
considerthefractionaldiffusionequation. Weinvestigate
1 r2
its integral solution Eq. (1) as well as the Fox function W (r,s) G(r,s)= exp . (6)
solution found in [26]. In Sec. (IV) we show how to →s→∞ (4πs)d/2 (cid:18)−4s(cid:19)
obtain solutions of the FFPE using Eq. (1). We give
We have used the assumption that the system is unbi-
as examples the motion of a particle in a uniform and
asedanduse convenientunits. For mostsystems andfor
harmonic force fields. In Sec. (V) a brief summary is
large times N (s,t) is concentrated on s = t/τ and
given. CT av
∞
τ = tψ(t)dt is the averaged pausing time. In this
av 0
case P(r,t) will become Gaussian when t . When
R → ∞
N (s,t) is broad a non-Gaussian behavior is found.
II. CONTINUOUS TIME RANDOM WALK CT
This is the case when τ = or in other words when
av
∞
the Laplace transform of ψ(t); ψ(u) behaves as
In the decoupled version of the continuous time ran-
dom walk (CTRW), a random walker hops from site to ψ(u) 1 uα+ 0<α 1, (7)
site and at each site it is trapped for a random time [2]. ∼ − ··· ≤
For this well known model two independent probability andwehaveusedconvenientunits. Shlesinger[3]showed
densitiesdescribetherandomwalk. Thefirstisψ(t),the that in this case an anomalous sub-diffusive behavior is
probability density function (PDF) of the independent found, r2 tα. Becausethe steps areindependent and
identically distributed (IID) pausing times between suc- based ohn cio∼nvolution theorem of Laplace transform
cessive steps. The second is the PDF f(r) for the IID
displacements of the random walker at each step. Thus N (s,u)=
CT
the CTRW describes a process for which the particle is
trapped on the originfor time τ1, it then jumps to r1, it 1 ψ(u)
is trapped on r1 for time τ2 and then it jumps to a new − exp sln[ψ(u)]
u { }∼
location, the process is then renewed. In what follows
we assume f(r) has a finite variance and a zero mean. uα−1exp( suα), (8)
Theasymptoticbehaviorofthe decoupledCTRWiswell −
investigated [2,3,31,34–36] and we now summarize some and N (s,u) is the Laplace transform of N (s,t).
CT CT
results fromthe CTRW literature which are of relevance Following[5]wereplacethesummationinEq. (5)with
to our work. integration and use Eq. (6) to find
Let P(r,t) be the PDF of finding the CTRW particle
∞
at r at time t. Let N (s,t) be the probability that s
CT P(r,u) n(s,u)G(r,s)ds (9)
steps are made in the time interval (0,t) so clearly ∼
Z0
∞ and according to Eq. (8)
N (s,t)=1 (3)
CT
s=0 n(s,u) uα−1exp( suα). (10)
X ≡ −
and the subscript CT denoted the CTRW. Because the Using the inverse Laplace transform we find
model is decoupled
∞
∞ P (r,t) n(s,t)G(r,s)ds (11)
∼
P(r,t)= NCT(s,t)W (r,s), (4) Z0
Xs=0 and n(s,t) is the inverse Laplace transform of n(s,u)
[5,32,33]
andW (r,s)istheprobabilitydensitythataparticlehas
reached r after s steps. In what follows we shall use the
1 t t
Fourier (r→k) and Laplace (t→u) transforms, we use n(s,t)= αs1+1/αlα s1/α (12)
the conventionthatthe argumentof a function indicates (cid:18) (cid:19)
in which space the function is defined, e.g., l (z) in Eq. (12) is a one sided L´evy stable probability
α
density whose Laplace transform is
∞
P (r,u)= N (s,u)W (r,s), (5)
CT ∞
Xs=0 lα(u)= e−uxlα(x)dx=e−uα. (13)
Z0
2
According to Eq. (11) the large time behavior of the A. CTRW Solution
CTRW solution is reached using an integral transforma-
tion of the Gaussian solution of ordinary diffusion pro- Let us consider an example of the (decoupled) CTRW
cesses [i.e., of G(r,s)]. As we shall see similar transfor- process. The solution of the CTRW in k,u space is [2]
mationscanbeusedtosolvethe FFPEandinparticular
the fractional diffusion equation. 1 ψ(u) 1
P (k,u)= − . (17)
The kernel n(s,t) is a non negative PDF normalized u 1 ψ(u)f(k)
−
according to
Usually CTRW solutions are found based on numerical
∞
inverse Fourier–Laplace transform of Eq. (17). Here we
n(s,t)ds=1, (14)
find an exact solution of the CTRW process for a spe-
Z0
cial choice of f(r) and ψ(t). Our solution is an infinite
itreplacestheCTRWprobabilityNCT(s,t). Noticethat sum of well known functions. We assume the PDF of
the PDFn(s,t), like NCT(s,t), is independent of the di- jumptimesψ(t)tobeonesidedL´evystabledensitywith
mensionalityoftheproblemd. Somepropertiesofn(s,t) ψ(u)=exp( uα) [40]. Displacements are assumedto be
are given in [32]. All moments of n(s,t) are finite and Gaussian an−d then
are given later in Eq. (87). It is easy to show that
n(s,t) = [dN(s,t)/ds] with N(s,t) defined in Eq. (2). 1 r2
W (r,s)= exp . (18)
n(s,t)iscalledthe “inverse”onesidedL´evystableprob- (4πs)d/2 −4s
(cid:18) (cid:19)
ability density [37].
Whenα=1,wefindn(s,t)=δ(s t)andhenceP(r,t) is exact, not only asymptotic. For this choice of PDFs
−
is Gaussian. This is expected since when α=1 the first the solution of the CTRW can be found explicitly. We
moment of pausing times τ is finite, therefore the law use
av
of large numbers is valid, and hence we expect that the
1 exp( uα)
number of steps s inthe randomwalkscheme will follow N (s,u)= − − exp( suα) (19)
CT
u −
s t/τ . When α < 1 the law of large numbers is
av
∼
not valid and instead the random number of steps s is andtheconvolutiontheoremofLaplacetransformtofind
described by n(s,t). Thus the transformation, Eq. (11),
has a meaning of a generalized law of large numbers. N0(t)=1 Lα(t)
−
Eq. (12)givesthekerneln(s,t)intermsofaonesided
t t
L´evy stable density. Schneider [39] has expressed L´evy N (s,t)=L L (20)
stable densities in terms of a Fox H function [41,42] CT α s1/α − α (s+1)1/α
(cid:18) (cid:19) (cid:18) (cid:19)
and
1 ( 1,1)
lα(z)= αz2H11,,10 z−1| ( 1−/α,1/α) . (15) t
(cid:20) − (cid:21)
L (t)= l (t)dt (21)
α α
Asymptotic behaviors of the one sided L´evy density can Z0
be found in Feller’s book [38] or based on the asymp-
is the one sided L´evy stable distribution. Inserting Eqs.
totic behaviors of the H Fox function [39,41,42]. For the
(18), (20) in Eq. (4) we find
casesα=1/3,1/2,2/3,closedformequationsintermsof
known functions, may be found in [39] (notice that [39] P(r,t)=[1 L (t)]δ(r)+
α
−
pointsoutrelevanterrorsintheliteratureonL´evystable
densities).
∞
Important for our purposes is the result obtained by t t
L L
Tunaley[34]alreadyin1974whichexpressedtheasymp- s=1( α(cid:20)s1/α(cid:21)− α"(s+1)1/α#)×
totic behaviorofthe CTRW solutionP (r,t), in terms of X
its Fourier–Laplace transform, as shown in [34] 1 r2
exp . (22)
uα−1 (4πs)d/2 (cid:18)−4s(cid:19)
P(k,u) . (16)
∼ uα+k2
The first term on the right hand side describes random
walks for which the particle did not leave the origin
In Appendix A we verify that Eq. (16) is indeed the
within the observation time t, the other terms describe
Fourier–Laplace transform of Eq. (11). The inversion of
random walks where the number of steps is s. In Fig. 1
equation (16) was accomplished by Tunaley [34] in one
we show the solution of the CTRW process in a scaling
dimensionandbySchneiderandWyss[26]indimensions
form. Weconsiderathreedimensionalcase,α=1/2and
two and three (see more details below).
use
3
1 tαm
L (t)= Γ[1/2,1/(4t)]. (23) C Γ(1+m) . (27)
1/2 √π m,d Γ(1+αm)
Here Γ(, ) denotes the incomplete Gamma function. To derive Eq. (27) we used the small u expansion of the
· ·
The figure shows r3P(r,t) versus the scaling variable Laplace transform of Eq. (25) and Tauberian theorem.
ξ =r2/tα. Forlargetimestthesolutionconvergestothe Later we shall show that the moments in Eq. (27) are
asymptoticsolutionfoundbasedontheintegraltransfor- identicaltothemomentsobtaineddirectlyfromtheinte-
mation Eq. (11). gral transformation Eq. (11). Hence our interpretation
of Eq. (11) as the asymptotic solution of the CTRW is
0.05 justified,[howeverourderivationofEq. (27),isbasedon
a specific choice of ψ(t) and f(r)].
0.045
0.04
III. FRACTIONAL DIFFUSION EQUATION
0.035
0.03
The fractionaldiffusionequationdescribes the asymp-
P(r,t)0.025 totic behavior of the CTRW Eq. (11). Balakrishnan
3r [43]hasderivedafractionaldiffusionprocessbasedupon
0.02
a generalization of Brownian motion in one dimension.
0.015 Schneider and Wyss [26,44] have formulated the follow-
ing fractional diffusion equation describing this process
0.01
0.005 ∂P (r,t)
= D1−α 2P (r,t), (28)
∂t 0 t ∇
0
0 10 20 30 40 50 60 70 80 90
ξ
FIG.1. Weshowr3P(r,t)versusξ=r2/tα fortheCTRW and 0 < α < 1. The fractional Riemann Lioville deriva-
tive D1−α in Eq. (28) is defined [25]
processandforα=1/2,d=3. Thecurvesinthefigure,t=2 0 t
(dotdashed),t=20(dotted),t=200 (dashed)andt=2000
(stars) are converging in the limit t → ∞ to the asymptotic D1−αZ(t)= 1 ∂ tdt′ Z(t′) . (29)
mastercurve(solid). Thismastercurvewascomputednumer- 0 t Γ(α)∂tZ0 (t−t′)1−α
ically basedonEq. (11). Itisdifficulttodistinguishbetween
The Fourier–Laplace solution of the fractional diffusion
the results for t = 200,2000 and the asymptotic result. Not
equation, with initial conditions P(r,t=0)=δ(r),
shown is thedelta function contribution at r=0.
uα−1
Let us calculate the Cartesian moments P (k,u)= (30)
uα+k2
M(2m , ,2m )=
1 ··· d is identical to the right hand side of Eq. (16). Therefore
∞ ∞ the integral solution of fractional diffusion equation is
dx dx x2m1 x2mdP(x , ,x ,t) (24)
1··· d 1 ··· d 1 ··· d ∞
Z−∞ Z−∞ P(r,t)= n(s,t)G(r,s)ds. (31)
with non negative integers m1,m2, . Clearly the Z0
···
odd moments are equal zero and from normalization
ThussolutionEq. (11)whichwasonlyasymptoticinthe
M(0, ,0)=1. In Appendix B we find
CTRWframeworkisexactwithinthefractionaldiffusion
···
equation approach. The integral solution Eq. (31) for
M(2m , ,2m )=
1 ··· d d = 1 was investigated by Saichev and Zaslavsky, we
∞ investigatealsothecasesd=2,3whichexhibitbehaviors
t t
C L L sm (25) different than the one dimensional case.
m,ds=1( α(cid:20)s1/α(cid:21)− α"(s+1)1/α#) As mentioned, Schneider and Wyss [26] have inverted
X
Eq. (30) finding the solution to the fractional diffusion
with
equation in terms of rather formal Fox function. Later
22m 1 we shall show that moments of P(r,t) calculated based
C = Πd Γ m + (26)
m,d πd/2 i=1 i 2 ontheintegralsolutionEq. (31)indimensionsd=1,2,3
(cid:18) (cid:19) are identical to the moments calculated based upon the
andm= dm >0. IntheAppendix wealsoshowthat Foxfunctionsolution. Inthissenseweshowthatthetwo
i i
for t solutions are identical.
→∞P At this stage we have only shown that the two ap-
M(2m , ,2m ) proaches, asymptotic CTRW and fractional diffusion
1 d
··· ∼
4
equation are identical. The reader might be wondering Γ n+d 2nΓ(1+n/2)
rn(t) =Ω 2 tαn/2. (35)
why should we bother with the fractionalequationif the h i d Γ(1+αn/2)
(cid:0) (cid:1)
CTRW approach gives identical results. The situation
is some what similar to the standard diffusion equation For n=2
whichpredictsaGaussianevolution. Thediffusionequa- 2dtα
tion α = 1 is an asymptotic equation describing much r2(t) = . (36)
h i Γ(1+α)
more general random walks. The main advantage of the
diffusion equation over a random walk approach is its When α = 1 the moments in Eq. (35) are identical to
simplicity. Solutions with special boundary conditions the Gaussian moments in Eq. (33).
(reflecting/absorbing) are relatively simple and usually In Appendix C we calculate the Cartesian moments
capture the essence of the more complex random walks. defined in Eq. (24). We find
Another extension of the diffusion equation is the diffu-
sioninexternalfieldasdescribedbyFokker–Planckequa- Γ(1+m)
M(2m , ,2m )=C tαm (37)
tions,suchanextensionforrandomwalksiscumbersome. 1 ··· d m,dΓ(1+αm)
The same is true for the fractional diffusion equation. It
can serve as a phenomenological tool describing anoma- whichisidenticaltotherighthandsideofequation(27).
lous diffusion. As we shall show in the next section we Thus the moments of the integral solution of fractional
can include the effects of an external field. On the other diffusionequationareidenticalto the asymptotic behav-
hand CTRW by definition is not built to consider an ex- ior of the moments of the CTRW found in the previous
ternal field. section.
Eq. (37) was also found in [26] based on Fox function
solution of fractional diffusion equation (see details next
A. Integral Solution subsection). Since the integral solution Eq. (31) gives
identical moments to those found using the Fox function
We first notice that the integral solution Eq. (31) solutiononecanassumethatthetwosolutionsareidenti-
shows that P (r,t) is normalized and non-negative. The cal. Itwouldbe interestingtoshowinamoredirectway
normalization is easily seen with the help of Eq. (14) that the Fox function solution is identicalto the integral
andthe nonnegativityofP (r,t)isevidentbecauseboth solution. In Appendix D we use a theorem on Fox func-
n(s,t) and G(r,s) are non negative. tions and attempt to give such a proof. Unfortunately
Calculation of spherical moments rn(t) is now con- we fail. The interestedreaderis referredto Appendix D.
h i
sidered. The calculation follows two steps, the first is to
calculate rn(s) forGaussiandiffusion(orfindtheGaus-
h i B. Fox Function Solution
sianmomentsinatextbook)thenusingthetransforma-
tion defined in Eq. (31) we find the moments rn(t) for
the fractional case. More precisely in Laplacehu spacie In [26], the Mellin transform was used to find the so-
lutionofthefractionaldiffusionequationintermsofFox
∞
H function
rn(u) = dsn(s,u) G(r,s)rndr , (32)
h i
Z0 (cid:20)Z (cid:21) P (r,t)=α−1π−d/2r−d
×
the Gaussian spherical moments are
(1,1)
H20 2−2/αr2/αt−1 (38)
G(r,s)rndr=ΩdΓ n+d 2nsn/2 (33) 12(cid:18) |(d/2,1/α),(1,1/α)(cid:19)
2
Z (cid:18) (cid:19)
The asymptotic expression for this solution is [26]
with Ω =1/√π, Ω =1andΩ =2/√π. The moments
1 2 3
for the fractional case are found using Eq. (32) P (r,t) καr−dξ2(2d−α) exp λ1ξ2−1α , (39)
∼ −
rn(u) = (cid:16) (cid:17)
h i where ξ r2/tα is the scaling variable,
≡
n+d ∞ κα =π−d/22−2−dα (2 α)−1/2α[α(d+1)/2−1]/(2−α)] (40)
Ω 2nΓ dssn/2n(s,u)= −
d
2
(cid:18) (cid:19)Z0 and
n+d
Ω 2nΓ Γ(1+n/2)u−1−αn/2. (34) λ =(2 α)αα/(2−α)2−2/(2−α). (41)
d 1
2 −
(cid:18) (cid:19)
Eq. (39) is valid for ξ >>1. For a brief introduction to
Using the inverse Laplace transformation we find
Fox functions and its application see [46].
5
Thebehaviorforξ <<1isfoundinAppendixEbased solution are not tabulated and hence from a practical
on the asymptotic expansion of the H function [41,42]. point of view one has to use numerical methods to find
For dimensions d=1 we find the solution. Also the convergence of the asymptotic so-
lution for ξ << 1 is not expected to be fast and hence
1 ∞ ( 1)nξn/2
P (r,t)= − (42) the integral solution is of special practical use.
2tα/2 n!Γ[1 α(n+1)/2] We consider the case α=1/2 and d=3. Using
n=0 −
X
and for d=3 1
l (z)= z−3/2exp[ 1/(4z)] (47)
1 ∞ ( 1)nξn/2 1/2 2√π −
P (r,t)= − . (43)
4πt3α/2ξ1/2 n!Γ[1 α(1+n/2)] withz >0. WefindtheintegralsolutionP(r,t)Eq. (31)
n=0 −
X
using numericalintegration. We haveused Mathematica
The leading terms in these expansions are for d=1
whichgaveallthenumericalresultswithoutdifficulty. In
1 1 Fig. 2weshowP(r,t)versusronasemilogplotandfor
P (r,t)= (44)
2Γ(1 α/2)tα/2 ··· different times t. Close to the origin r = 0 we observe a
− sharp increase of P(r,t) as predicted in Eq. (45).
and for d=3 More detailed behavior of P(r,t) is presented in Fig.
1 r−1 3 where we show rP(r,t) versus r. In Fig. 3 we also
P(r,t)= + . (45) exhibit linear curves based on the asymptotic expansion
4πΓ(1 α) tα ···
− Eq. (43) which predicts
Weseethatford=1andα=1,P(r =0,t)=1/(2√πt),
rP(r,t) C (t) rC (t) (48)
as expected for this normal case. We also see that for ∼ 1 − 2
d = 3, α = 1 and when r 0 the solution diverges like
6 → where
P(r,t) 1/r. ThisbehaviorisnotunphysicalandP(r,t)
∼
is normalized according to 4π P(r,t)r2dr =1. 1
C (t)= (49)
Ford=2theasymptoticexpansionoftheFoxfunction 1 4π3/2t1/2
R
isnotknown(seemoredetailsinAppendixE).UsingEq.
and C (t) = 1/[4πΓ(1/4)t3/4]. Eq. (48) is valid when
(31) one can show that for d=2, and ξ <<1 2
ξ = r2/t1/2 << 1 and as expected we see in Fig. 3
P(r,t) 1 ln tα/2/r . (46) that forthis case the numericalintegralsolutionandthe
∼ πΓ(1 α)tα asymptotic solution (48) agree well.
− h i
Eqs. (44-46)werederivedindependentlybyA.I.Saichev
102
[45].
The CTRW behavior on the origin is different than
101
what we have found for the fractional diffusion equation
Eqs. (44-46). Within the CTRW on the decay on the 100
origin is described with the first term on the right hand
side of Eq. (22), [1 L (t)]δ(r). This term describes 10−1
α
−
random walks for which the particle is trapped on its
initial location during the time of observation t. This P(r)10−2
CTRW term decays like t−αδ(r) for long times and no
10−3
such singular term is found in the solution of the frac-
tional. In dimensions d = 2,3 the decay of the CTRW 10−4
singular term is as slow as the decay of P(0,t) found
from the fractional diffusion equation. Hence on the ori- 10−5
gin the fractionaldiffusion approximationdoes not work
well. In contrast, for normal random walks, the singu- 10−6
0 0.5 1 1.5 2 2.5 3 3.5
lar term decays exponentially with time and then the r
FIG.2. P(r,t)versusr onasemilogplotandfordifferent
diffusion approximation works well already after an ex- times t = 0.02 (dotted), t = 2 (dash dotted) and t = 200
ponentially short time. (dashed). Note that P(r = 0,t) = ∞ hence the solution in
figureis cutoff close to r→0.
C. Example
Here we use the integral solution of the fractional dif-
fusion equationfor a specific example. The Fox function
6
0.5 0.12 FIG.4. Scaling behavior of P(r,t) for α=1/2 and d=3.
0.4 t=0.02 0.1 t=0.2 We show r3P(r,t) versus ξ = r2/tα for times t = 0.02 (dot-
0.3 ted), t=2 (dash dotted) and t=200 (dashed), due to exact
0.08
r P(r,t) 00..12 0.06 stcharleiengcuorfvetshienstohluetfiiognuroen.eAclsaonnshotowdnist(isnogliudischurbveetsw)eaerne tthhee
0.04
0 two asymptotic behaviors valid for small and large values of
−0.1 0.02 the scaling variable ξ. For ξ >>1 we used Eq. (39) and for
−0.2 0
0 0.5 1 1.5 0 0.5 1 1.5 ξ<<1 we used the leading term in the expansion Eq. (43).
r r
x 10−3
0.04
0.035 t=2 3.2 t=200 IV. FRACTIONAL FOKKER–PLANCK
EQUATION
P(r,t) 0.03 3
r 0.025
Let us consider the fractional Fokker–Planckequation
2.8
0.02 (FFPE) [30] describing the stochastic evolution of a test
0.015 2.6 particle under combined influence of external force field
0 0.5 1 1.5 0 0.5 1 1.5
r r F(x) and a thermal heat bath. The equation reads
FIG. 3. Behavior of P(r,t) close to the origin r = 0 for
α = 1/2 and d = 3. We show rP(r,t) versus r for different
P˙ (x,t)=K D1−αLˆ P(x,t) (50)
times t = 0.02,0.2,2,200 (dashed curves). Because P(r,t) α 0 t fp
diverges when r → 0 like P(r,t) ∼ 1/r, the figures shows
that rP(r,t) → C1(t). Also shown (linear curves) the ap- and
proximation Eq. (48), as expected when ξ =r2/tα <<1 the ∂ F(x) ∂2
approximation is in good agreement with the integral solu- Lˆfp =−∂x k T + ∂x2 (51)
tion. b
is the Fokker–Planck operator. K is a generalized dif-
α
FinallyinFig. 4weshowourresultsinascalingform, fusion constant and T is temperature. The fractional
similar to the way we presented the CTRW results in derivative D1−α was defined in Eq. (29). We consider
0 t
Fig. 1 which shows the CTRW solution for both short theonedimensionalcaseandextensionstohigherdimen-
andlongtimes. Wepresentr3P(r,t)versusξfordifferent sions are straight forward. When F(x) = 0 the FFPE
choices of time t. We observecollapse of all curves,both Eq. (50)reducestothefractionaldiffusionequation(28).
for shorter and longer times, onto one master curve. We Whenα=1werecovertheordinaryFokker-Planckequa-
also show the asymptotic behaviors ξ >> 1 and ξ << tion. The stationary solution of the FFPE is the Boltz-
1, Eq. (39) and Eq. (43) respectively. The numerical mann distribution and the equation is compatible with
integralsolutionEq. (31)agreeswellwiththeasymptotic linear response theory [30].
behaviors in the appropriate regimes. Comparing Fig. InthissectionwefindanintegralsolutionoftheFFPE.
(1)andFig. (4)weseethatthe solutionofthefractional We show that the solution is normalized and non neg-
diffusionequationapproximatestheexactsolutionofthe ative, an issue not discussed in [30]. Without loss of
CTRW very well when t and r=0. generality we set K = 1. We consider initial condition
→∞ 6 α
P(x,t=0)=δ(x x ).
0
−
0.05 Firstlet usshowthatthe solutionis normalized. Inte-
grating Eq. (50) with respect to x and using the bound-
0.045
ary conditions
0.04
∂P(x,t)
0.035 P(x,t) = =0 (52)
x=±∞ x=±∞
| ∂x |
0.03
P(r,t)0.025 we find
3r ∞
0.02 P˙(x,t)dx=0, (53)
Z−∞
0.015
meaning that normalization is conserved as it should.
0.01
We write the solution of the FFPE Eq. (50) in terms
0.005 of an integral of a product of two functions
00 10 20 30 40 50 60 70 80 90 ∞
ξ P(x,t)= dsn˜(s,t)P (x,s) (54)
1
Z0
where
7
∂P1(x,s) =Lˆ P (x,s). (55) Thusn˜(s,u)=n(s,u)foundinthecontextofthesolution
FP 1
∂s of CTRW Eq. (10).
We see that the integral solution of the FFPE has a
P (x,s) is a normalized solution of the ordinary Fokker-
1 similar structure as the solutions of the fractional diffu-
Planck equation with initial conditions P (x,s = 0) =
1 sion equation or equivalently of the asymptotic CTRW.
δ(x x ). Methods of solution of Eq. (55) are given in
− 0 The solution is Eq. (54) with n(s,t) defined in Eq. (12)
Risken’s book [47]. In what follows we shall prove that
andP (x,s)beingthesolutionofcorrespondingordinary
n˜(s,t)=n(s,t), defined in Eq. (12). 1
Fokker-Planckequation. Itiseasytounderstandthatthe
Weusethe LaplacetransformofEq. (54)andnormal-
solution, P(x,t) is normalized and non negative because
ization condition of P(x,t) to show
both P (x,t) and n(s,t) are normalized PDFs.
1
∞ In [30] a different method of solution of the fractional
n˜(s,u)ds=1/u. (56) Fokker–Planck equation was investigated. It can be
Z0 shownthat P(x,t) canbe expanded inan eigenfunction
Hencen˜(s,t)isnormalizedaccordingto ∞n˜(s,t)ds=1. expansion,whichissimilartothestandardeigenfunction
0
The Laplace transform of Eq. (50) reads expansion of solution of ordinary Fokker–Planck equa-
R tion, however eigen modes decay according to Mittag–
uP(x,u) δ(x x )=u1−αLˆ P (x,u), (57) Lefflerrelaxationinsteadofthestandardexponentialde-
0 fp
− −
cay of the modes in the integer Fokker–Planckequation.
inserting Eq. (54) in Eq. (57) we find It can be shown that the integral solution investigated
here is identical to the eigen function expansion in [30].
∞
u n˜(s,u)P (x,s)ds δ(x x )=
1 0
− −
Z0
A. Example 1, Biased Fractional Wiener Process
∞
u1−α n˜(s,u)Lˆ P (x,s)ds, (58)
fp 1
Z0 Considerthebiasedfractionaldiffusionprocessdefined
withageneralizeddiffusioncoefficientK andauniform
integrating by parts using Eq. (55) we find α
force field F(x) = F > 0. For this case the mean dis-
∞ placement grows slower than linear with time according
u n˜(s,u)P (x,s)ds δ(x x )=
1 0 to
− −
Z0
FK tα
α
x(t) = . (64)
h i k TΓ(1+α)
u1−α[n˜( ,u)P (x,s= ) n˜(0,u)P (x,0)] b
1 1
∞ ∞ −
ThewellknownsolutionoftheordinaryFokker-Planck
u1−α ∞ ∂ n˜(s,u) P (x,s)ds. (59) equation is
1
− ∂s
Z0 (cid:20) (cid:21) x F˜s 2
1
FromEq. (56)n˜(∞,u)=0andsinceP1(x,0)=δ(x−x0) P1(x,s)= 4πsexp−(cid:16) −4s (cid:17) . (65)
we may rewrite Eq. (59)
 
∞ un˜(s,u)+u1−α ∂ n˜(s,u) P (x,s)ds= Eq. (65) describes a biased Wiener process and F˜ =
1
Z0 (cid:26) (cid:20)∂s (cid:21)(cid:27) F/(kbT). The solution P(x,t) for the fractional case is
found using the transformation Eq. (54). In Laplace u
1 u1−αn˜(0,u) δ(x x ). (60)
0 space the solution is
− −
(cid:2) (cid:3)
Eq. (60) is solved once both of its sides are equal zero; P(x,u)=
therefore two conditions must be satisfied, the first
F˜uα−1τα F˜(x 1+4(uτ)α x)
exp − | | (66)
n˜(0,u)=uα−1, (61) 1+4(uτ)α " p 2 #
and the second andpτα = 1/(F˜K )2. For (τu)α << 1, corresponding to
α
the long time behavior of the solution P(x,t), we find
∂
n˜(s,u)= uαn˜(s,u). (62)
∂s − F˜uα−1ταexp F˜ταuαx x>0
−
The solution of Eq. (62) with initial condition Eq. (61) P(x,u)  (cid:16) (cid:17)
is ∼F˜uα−1ταexp F˜ x F˜ xταuα x 0.
− | |− | | ≤
n˜(s,u)=uα−1exp(−uαs). (63)  (cid:16) (cid:17) (67)
8
Since 0 P(x,u) uα−1andhenceaccordingtoTaube- 0.3
−∞ ∼
rian thReorem −0∞P(x,t)∼1/tα we have P(x,t)=0 for 0.02.52 t=2 0.15 t=10
x<0 whent . Therefore,using the inverseLaplace
transformofE→Rq.∞(67)the asymptoticbehaviorof P(x,t) P(x,t)0.15 0.1
0.1
is 0.05
0.05
αA11/αx1+t1/αlα A1/αtx1/α 0<x −05 0 5 10 −05 0 5 10 15
lim P(x,t)= x x
t→∞  (cid:0) (cid:1)
0 0>x,
 0.04 0.02
 (68) 0.03 t=250 0.015 t=103
andA=F˜τα. IntegratingEq. (68)we find the distribu- P(x,t)0.02 0.01
tion function
0.01 0.005
x t
tl→im∞Z−∞P(x,t)dx=1−Lα(cid:18)A1/αx1/α(cid:19), (69) FI0G.05. Th20e dxyn40amic6s0of P(x,t)0 fo0r the5f0raxctio10n0al d1i5ff0u-
valid for x > 0. Eq. (69) was derived also in [37] based sion process in external field F with α = 1/2. We show
the solution for four different times. The dotted curve is
on the biased CTRW, thus as expected the solution of
the asymptotic solution valid for large t. The maximum of
the fractional Fokker-Planck equation converges to the
P(x,t) is for all times on theinitial condition x=0. We use
solution of the CTRW in the limit of large t.
On the origin one can use Tauberian theorem to show F =kbT =K1/2 =1.
A
P(0,t) t−α (70) B. Example 2, the Fractional Ornstein-Uhlenbeck
∼ Γ(1 α)
− Process
valid for long times. For the case F = 0 we have found
P(0,t) t−α/2 so as expected the decay on the origin is We consider as a second example the fractional
fasterfo∼rthebiasedcasesinceparticlesaredriftingaway Ornstein–Uhlenbeck (OU) process, namely the motion
from the origin. ofthe testparticle inharmonicoscillator. This casecan-
In Fig. 5 we presentthe solutionfor the case α=1/2, not be analyzed using the CTRW in a direct way. The
then CTRW formalism considers only uniformly biased ran-
dom walks and the fractional processes in non-uniform
P(x,t)= fields are not uniformly biased. We consider α = 1/2,
F(x)/k T = x and use the well known solution of the
b
1 ∞ 1 s2 (x F˜s)2 ordinary OU−process [47,48]
ds exp −
tK12/2π Z0 √4πs "−4K12/2t − 4s # 1 (x x e−s)2
0
q (71) P1(x,s)= 2π(1 e−2s)exp"− 2(−1 e−2s) #. (73)
− −
p
which is evaluated numerically. For large times we have The solution of the fractional OU process is then found
using numerical integration of Eq. (54) using Eqs. (47)
A A2x2
and (73). Our results are presented in Fig. 6. We have
P(x,t) exp (72)
∼ √πt − 4t considered an initial condition x = 1/2 and we observe
(cid:20) (cid:21) 0
a strong dependence of the solution on the initial con-
for0<x. Asseeninthefiguretheexactresultexhibitsa
dition. A cusp on x = x is observed for all times t,
0
strong sensitivity on initial condition and the maximum
thus the initial condition has a strong influence on the
ofP(x,t)islocatedonx=0. Thisisdifferentthanordi-
solution. The solution approaches the stationary Gaus-
nary diffusion process in which the maximum of P(x,t)
sian shape slowly in a power law way and the solution
is on x(t) . The curves in Fig. 5 are similar to those
h i deviates from Gaussian for any finite time. Unlike the
observedby Scher and Montroll [4] based on lattice sim-
ordinary Gaussian OU process, the maximum of P(x,t)
ulation of CTRW and also by Weissman et al [36] who
is not on the average x(t) but rather the maximum is
investigatedbiasedCTRWusingananalyticalapproach. h i
for short times located on the initial condition.
The FFPE solution presented here is much simpler than
the CTRW solution, still it captures all the important
features of the more complex CTRW result.
9
1.4 1.4 anditcanserveas apracticaltoolfor finding solutionof
1.2 1.2 certain fractional kinetic equations.
t=0.02 t=0.2
1 1
P(x,t)00..68 00..68
VI. ACKNOWLEDGMENTS
0.4 0.4
0.2 0.2
0 0 EB thanks A. I. Saichev and G. M. Zaslavsky for cor-
−2 0 2 −2 0 2
x x respondence and J. Klafter, R. Metzler and G. Zumofen
for discussions.
1.4 1.4
1.2 1.2
t=2 t=20
1 1
VII. APPENDIX A
P(x,t)00..68 00..68
0.4 0.4 We rewrite Eq. (16)
0.2 0.2
∞
0 −2 0 2 0 −2 0 2 P (k,u) uα−1 e−s(uα+k2)ds, (74)
x x ∼
FIG. 6. The dynamics of P(x,t) for the fractional Orn- Z0
stein–Uhlenbeck process and α=1/2. We show the solution the inverse Fourier transform of Eq. (74) is
for four different times. The dotted curve is the stationary
solution, i.e. the Boltzmann distribution which for the Har- P (r,u) uα−1 −1 ∞e−s(uα+k2)ds , (75)
monic potential is Gaussian. Notice the cusp on the initial ∼ F
condition x=x0=1/2. (cid:26)Z0 (cid:27)
and −1 is the inverse Fourier transform. Changing the
F
Like the ordinary OU process the fractional OU pro- order of integration over parameter s and the operation
cess has a special role. The ordinary OU process de- −1 (this is later justified by the identity of moments of
F
scribestwotypesofbehaviors,thefirstisanoverdamped theintegralsolutionandthatfoundin[26])wefindusing
motion of a particle in harmonic potential, the second Eq. (8)
is the velocity of a Brownian particle modeled by the
∞ 1
Langevinequation,thelatteristhebasisfortheKramers P(r,u) n(s,u) exp r2/4s ds. (76)
equation. InasimilarwaythefractionalOUdescribesan ∼Z0 (4πs)d/2 −
(cid:0) (cid:1)
overdampedandanomalousmotionin Harmonicpoten-
Applying the inverseLaplacetransformto this equation,
tial considered in this section and in [30], it can also be
changing the order of integrations over s and the inverse
usedtomodelthevelocityofaparticleexhibitingaL´evy
Laplace operation, we find Eq. (11).
walk type of motion [33]. The fractional OU process is
the basis of the fractional Kramers equation introduced
recentlybyBarkaiandSilbey[33],thisequationdescribes VIII. APPENDIX B
superdiffusionwhileinthisworkwehaveconsideredsub
diffusion.
We investigate the Cartesian moments of the CTRW.
These moments are
V. SUMMARY M(2m , ,2m )=
1 d
···
Fractional diffusion equation is an asymptotic equa- ∞ ∞
tion which predicts the behavior of the decoupled con- dx1··· dxdx21m1xm2 2···x2dmd×
tinuous time random walk in the sub diffusive regime. Z−∞ Z−∞
The fractional Fokker–Planck equation considers such a
∞
sub-diffusive motion in an external force field and close t t
L L
to thermal equilibrium. In this work we have considered s=1( α(cid:20)s1/α(cid:21)− α"(s+1)1/α#)×
anintegraltransformationwhichgivesthesolutionofthe X
fractionalFokker–Planckequationintermsofsolutionof 1 x2+ x2
exp 1 ··· d . (77)
ordinary Fokker–Planck equation. Solution of ordinary (4πs)d/2 − 4s
(cid:18) (cid:19)
Fokker–Planckequation can be found based on different
Changing order of integration and summation, we use
analytical and numerical methods [47–51]. The integral
the identity
transformation describes also the long time behavior of
the CTRW in dimension d = 1,2,3. Thus the transfor- ∞ ∞ 1 x2+ x2
mationmapsGaussiandiffusionontofractionaldiffusion, Z−∞dx1···Z−∞dxd(4πs)d/2 exp(cid:18)− 1 4·s·· d(cid:19)×
10