Stretched-Gaussian asymptotics of the truncated Lévy flights for the diffusion in nonhomogeneous media
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] J un Stretched-Gaussian asymptotics of thetruncated L´evy flights for the diffusion innonhomogeneous media
Tomasz Srokowski
Institute of Nuclear Physics, Polish Academy of Sciences, PL – 31-342 Krak´ow,Poland
Abstract
The L´evy, jumping process, defined in terms of the jumping size distributionand the waiting time distribution, is considered. The jumping rate depends onthe process value. The fractional diffusion equation, which contains the variablediffusion coefficient, is solved in the diffusion limit. That solution resolves itselfto the stretched Gaussian when the order parameter µ →
2. The truncation ofthe L´evy flights, in the exponential and power-law form, is introduced and thecorresponding random walk process is simulated by the Monte Carlo method. Thestretched Gaussian tails are found in both cases. The time which is needed to reachthe limiting distribution strongly depends on the jumping rate parameter. Whenthe cutoff function falls slowly, the tail of the distribution appears to be algebraic.
Key words:
Diffusion; Fractional equation; Truncated L´evy flightsPACS 02.50.Ey, 05.40.Fb, 05.60.-k
Transport in physical systems can be described in terms of a jumping processwhich is completely defined by two probability distributions. They determinespatial and temporal characteristics of the system and are called the jumpingsize, λ ( x ), and waiting time, w ( t ), probability distributions, respectively. Thestochastic trajectory consists of infinitely fast jumps which are separated byintervals when the Brownian particle is at rest (e.g., due to traps). Usually,the distributions λ and w are regarded as mutually independent and processeswhich they generate are studied in the framework of the decoupled continuoustime random walk (CTRW). If the distribution w ( t ) is Poissonian, the jumpstake place at random, uniformly distributed times and the corresponding pro-cess is Markovian. The algebraic form of w ( t ), which possesses long tails, is Preprint submitted to Elsevier Science 3 November 2018 f particular interest; it leads to long rests and then to a weak, subdiffusivetransport [1]. The jumping size distribution λ can assume – in order to be sta-ble – either the Gaussian form, or obey the general L´evy distribution whichcorresponds to processes with long tails, frequently encountered in many areasof physics, sociology, medicine [2] and many others. In the diffusion limit ofsmall wave numbers, the master equation for the decoupled CTRW, and forthe case of the L´evy distributed jumping size, resolves itself to the fractionaldiffusion equation with a constant diffusion coefficient.However, many physical phenomena require taking into account that the spacehas some structure and that the diffusion coefficient must actually be variable.It is the case when one considers the transport in porous, inhomogeneous me-dia and in plasmas, as well as for the L´evy flights in systems with an externalpotential [3]. In the Hamiltonian systems, the speed of transport depends onregular structures in the phase space [4]. Any description of the diffusion onfractals must involve the variable diffusion coefficient, and that one in thepower-law form [5,6]. Similarly, the diffusion on multifractal structures can beregarded as a superposition of the solutions which correspond to the individ-ual fractals [7]. The power-law form of the diffusion coefficient has been alsoused to describe e.g. the transport of fast electrons in a hot plasma [8] andthe turbulent two-particle diffusion [9]. The presence of long jumps indicates ahigh complexity and the existence of long-range correlations. One can expectthat especially in such systems the waiting time depends on the position. Thejumping process, stationary and Markovian, which takes into account thatdependence, has been proposed in Ref.[10]: the distribution w ( t ) is Poissonianwith a x -dependent jumping rate ν ( x ). In this paper, we consider that processfor the L´evy distributed λ ( x ) and solve the corresponding fractional equa-tion. The solution consistently takes into account that the symmetric L´evydistribution, defined by the characteristic function exp( −| k | µ ), 0 < µ ≤ ∼ | x | − µ − , for 0 < µ < µ = 2.The L´evy process is characterized by very long jumps since the tails of theL´evy distribution fall slowly and the second moment is divergent. In a concrete,realistic problem, however, the available space is finite and the asymptoticsdiffers from the L´evy tail. One can take that into account by introducing acutoff to the stochastic equations by multiplying the jumping size distribution– in the L´evy form – by a function which falls not slower than 1 / | x | , usuallyby the exponential or the algebraic function. Since the resulting probabilitydistribution has the finite variance, in the diffusion approximation the processis equivalent to the Gaussian limit µ → µ → µ →
2. Sec.III is devoted to thetruncated L´evy flights: we present the probability distributions, obtained fromsimulations of random walk trajectories, for both exponential and power-lawform of the cutoff. The results are summarized in Sec.IV.
We consider the random walk process defined by the waiting time probabilitydistribution w ( t ) and the jump-size distribution λ ( x ). They are of the form w ( t ) = ν ( x )e − ν ( x ) t , (1)where ν ( x ) = | x | − θ ( θ > −
1) and λ ( x ) = q /π ∞ Z exp( − K µ k µ ) cos( kx ) dk, (2)respectively. The latter expression corresponds to the symmetric L´evy distri-bution. The master equation for the above process is of the form [10] ∂∂t p ( x, t ) = − ν ( x ) p ( x, t ) + Z ν ( x ′ ) λ ( | x − x ′ | ) p ( x ′ , t ) dx ′ . (3)In the diffusion limit of small wave numbers, the Eq.(3) can be reduced to thefractional equation. Indeed, the Fourier transform of jump-size distribution, e λ ( k ) = exp( − K µ | k | µ ), can be expanded e λ ( k ) = 1 − K µ | k | µ + . . . and themaster equation in the Fourier space becomes the following equation [15] ∂ e p ( k, t ) ∂t = − K µ | k | µ e F [ ν ( x ) p ( x, t )] . (4)3he above approximation means that the summation over jumps is substi-tuted by the integral and it agrees with the exact result if the jumps aresmall, compared to the entire trajectory. One can demonstrate by estimatingthe neglected terms in the Euler-Mclaurin summation formula that the ap-proximation fails near the origin and for µ ≤ µ > ∂p ( x, t ) ∂t = K µ ∂ µ [ ν ( x ) p ( x, t )] ∂ | x | µ , (5)where ∂ µ /∂ | x | µ is the Riesz-Weyl derivative, defined for 1 < µ < ∂ µ ∂ | x | µ f ( x, t ) = −
12 cos( πµ/ − µ ) ∂ ∂x ∞ Z −∞ f ( x ′ , t ) | x − x ′ | µ − dx ′ . (6)Eq.(5) for the constant diffusion coefficient – often generalized to take into ac-count non-Markovian features of the trapping mechanism in the framework ofCTRW by substituting the simple time derivative by the fractional Riemann-Liouville derivative [1] – is frequently considered and solved by means of avariety of methods. The solution resolves itself to the L´evy stable distributionwith the asymptotic power-law x -dependence and divergent second moment.The method of solution which is especially interesting for our considerationsinvolves the Fox functions. The well-known result of Schneider [18] states thatany L´evy distribution, both symmetric and asymmetric, can be expressed as H , , . In order to solve the Eq.(5) for the case of the variable jumping rate, ν ( x ) = | x | − θ , we conjecture that the solution also belongs to the class of func-tions H , , and it has the scaling property. Then the ansatz is the following p ( x, t ) = N a ( t ) H , , a ( t ) | x | (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ( a , A ) , ( a , A )( b , B ) , ( b , B ) , (7)where N is the normalization constant. The Fox function is defined in thefollowing way [19,20] 4 mnpq z (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ( a p , A p )( b q , B q ) = 12 πi Z L χ ( s ) z s ds, (8)where χ ( s ) = Q m Γ( b j − B j s ) Q n Γ(1 − a j + A j s ) Q qm +1 Γ(1 − b j + B j s ) Q pn +1 Γ( a j − A j s ) . (9)The contour L separates the poles belonging to the two groups of the gammafunction in the Eq.(9). Evaluation of the residues leads to the well-knownseries expansion of the Fox function: H mnpq z (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ( a p , A p )( b q , B q ) = m X h =1 ∞ X ν =0 Q mj =1 ,j = h Γ( b j − B j b h + νB h ) Q nj =1 Γ(1 − a j + A j b h + νB h ) Q qj = m +1 Γ(1 − b j + B j b h + νB h ) Q pj = n +1 Γ( a j − A j b h + νB h ) ( − ν z ( b h + ν ) /B h ν ! B h , (10)We will try to solve the fractional equation (5) by inserting the function (7).This procedure, if successful, would allow us to find conditions for the coeffi-cients and the function a ( t ). The idea to assume the solution in the form (7)is motivated by the following property of the Fox function z σ H mnpq z (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ( a p , A p )( b q , B q ) = H mnpq z (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ( a p + σA p , A p )( b q + σB q , B q ) (11)which eliminates the algebraic factor by means of a simple shift of the coef-ficients. One can demonstrate that Eq.(5) cannot be satisfied, in general, bythe function (7) for any choice of the parameters. However, as long as we areinterested only in the diffusion limit of small wave numbers | k | (large | x | ), thehigher terms in the characteristic function expansion can be neglected. There-fore we require that the Eq.(5) should be satisfied by a function which agreeswith the exact solution only up the terms of the order | k | µ in the Fourier space.Note that this condition does not introduce any additional idealization sincethe Eq.(5) itself has been constructed on the same assumption: the higherterms in the | k | expansion of λ ( x ) have been also neglected.First, we need the Fourier transform of the Fox function. Since the process issymmetric, we can utilize the formula for the cosine transform which yields5lso a Fox function but of the enhanced order: ∞ Z H mnpq x (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ( a p , A p )( b q , B q ) cos( kx ) dx = πk H n +1 ,mq +1 ,p +2 k (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (1 − b q , B q ) , (1 , / , , (1 − a p , A p ) , (1 , / . (12)The function e p ( k, t ) is then of the order H , , . The Fourier transform of thefunction p θ = x − θ p ( x, t ) can be obtained in the same way. Next, we insert theappropriate functions into the Eq.(4) and expand both sides of the equation byusing the formula (10). Let us denote the expansion coefficients of the functions e p ( k, t ) and e p θ ( k, t ) by h σ,ν and h ( θ ) σ,ν , respectively; σ assumes the values 1 and2. Applying the Eq.(10) yields e p ( k, t ) = h , + h , | k | + h , | k | (1 − a ) /A − + h , | k | (2 − a ) /A − + h , k + . . . (13)and e p θ ( k, t ) = h ( θ )1 , + h ( θ )1 , | k | + h ( θ )2 , | k | (1 − a + θA ) /A − + . . . . (14)After inserting the above expressions to the Eq.(4), we find some simple rela-tions among the coefficients of the Fox function by comparison of the expo-nents. To get the term | k | µ on lhs, which corresponds to the term k on rhs, weneed the condition (2 − a ) /A − µ . Moreover, we attach the third term onrhs to the first one by putting (1 − a + θA ) /A − a = 1 − (1 − θ ) / ( µ + θ ) and A = 1 / ( µ + θ ).The coefficients h , and h ( θ )1 , vanish identically since the gamma function inthe denominators has its argument equal 0. The only remaining term can beeliminated by assuming the condition 1 − b − B (1 − θ ) = 0; then h , = 0since that term contains the function Γ(1 − b − B (1 − θ )) in the denominator.Therefore, such choice of the coefficients reduces the function (13) to the twoterms and it makes it identical with the probability distribution for the L´evyprocess in the diffusion limit. Finally, the Eq.(4) becomes a simple differentialequation which determines the function a ( t ):˙ ξ = K µ h ( θ )0 h , ξ − θ/µ , (15)6here h ( θ )0 = h ( θ )1 , + h ( θ )2 , and ξ ( t ) = a − µ . The solution a ( t ) = K µ h ( θ )0 h , θµ ! t − / ( µ + θ ) (16)corresponds to the initial condition p ( x,
0) = δ ( x ). The coefficient h , can bedetermined directly from Eq.(10), whereas h ( θ )0 = (2 π ) − a − θ e p (0 , t ) = π − a − θ R ∞ p ( x, t ) dx can be expressed in terms of the Mellin transform from the Fox function, χ ( s ),and then easily evaluated.In the limit µ →
2, the asymptotic behaviour of the fractional equationchanges qualitatively. It is no longer algebraic; the tails of the distribution,and consequently the tails of the Fox function, have to assume the exponen-tial form for µ = 2. It is possible only if the algebraic contribution from theresidues vanishes. If n = 0, all poles are outside the contour L and H m, p,q isgiven by the integral over a vertical straight line: H m, p,q = 12 πi w + i ∞ Z w − i ∞ χ ( s ) z s ds, (17)where w < Re ( b i /B i ). The above integral is usually neglected if n > b and B which allows us todetermine these coefficients: b = 1 − (1 − θ ) / (2 + θ ) and B = 1 / (2 + θ ).Note that the above choice of ( b , B ) yields h , = 0 in Eq.(13) and then thesolution of Eq.(5) is exact up to the order k for any µ . The most generalsolution of the fractional equation in the form of the function H , , , whichinvolves the required conditions, is the following p ( x, t ) = N a ( t ) H , , a ( t ) | x | (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (1 − − θµ + θ , µ + θ ) , ( a , A )( b , B ) , (1 − − θ θ , θ ) . (18)The coefficients in the main diagonal have a simple interpretation. ApplyingEq.(10) to (18) reveals that p ( x, t ) behaves as | x | b /B for | x | → b and B are responsible for the shape of the probabilitydistribution near the origin. On the other hand, the parameters a and A H mnpq z (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ( a p , A p )( b q , B q ) = H mnpq z (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (1 − b q , B q )(1 − a p , A p ) (19)and by expansion according to Eq.(10): the leading term is of the form p ( x, t ) ∼| x | (2 − a ) /A = | x | − − µ ( | x | → ∞ ). Now it becomes clear why our methodof solving the Eq.(5) did not determine the parameters ( b , B ). Since weneglected higher terms in the k -expansion, the region of small | x | remainedbeyond the scope of the approximation. However, we can supplement that so-lution by referring directly to the master equation (3) which reveals a simplebehaviour near the origin. That equation is satisfied by the stationary solu-tion 1 /ν ( x ) = | x | θ , for any normalizable λ ( x ). Obviously, such p ( x, t ) cannotbe interpreted as the probability density distribution since the normalizationintegral diverges in infinity but it properly reproduces that distribution forsmall | x | . Therefore we obtain the additional condition b = θB which im-proves the agreement of our solution with the solution of the master equation.The probability distribution for the case θ = 0 can be easily found by the di-rect solution of the fractional equation which is exact and yields the followingvalues of the parameters: b = 0, B = 1, a = 1 /
2, and A = 1 / µ = 2, the main diagonal in Eq.(18) can be eliminated and thesolution of the fractional equation (5) in the limit µ → p ( x, t ) = N a ( t ) H , , a ( t ) | x | (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ( a , A )( b , B ) . (20)The function a ( t ) is given by Eq.(16) in the following form a ( t ) = " K (2 + θ ) h h t − / (2+ θ ) , (21)where h = Γ( b + B (1 − θ )) / Γ( a + A (1 − θ )) and h = Γ( b +3 B ) / Γ( a +3 A )are the coefficients of the k -expansion of the functions | x | − θ p ( x, t ) and p ( x, t ),respectively. The asymptotic expression for p ( x, t ) can be obtained from theestimation of the integral (17) by the method of steepest descents [21]. The8esult reads H , , ( z ) ≈ e iπ ( α ′ − / E ( z e iπα ) , (22)where E ( z ) = 12 πiα ∞ X j =0 C j ( βα α z ) (1 − α ′ − j ) /α exp( βα α z ) /α (23)and α = B − A , α ′ = a − b + 1 / β = A A /B B , z = a | x | . The finalresult appears to be a stretched Gaussian, modified by an algebraic factor anda series which converges to a constant for z → ∞ : H , , ( z ) ≈ πα z (1 − α ′ ) /α exp( − β /α αz /α ) β (1 − α ′ ) /α α − α ′ ∞ X j =0 ( − j C j β − j/α α − j z − j/α . (24)The above expression has been obtained by Wyss [22] as the expansion of H , , ( z ) which satisfies a generalized diffusion equation. It is the integral equa-tion in respect to the time variable (non-Markovian) which resolves itself tothe fractional diffusion equation; that equation is commonly used to handlethe subdiffusive processes in the framework of the CTRW [12,23,1]. The func-tion which contains the stretched Gaussian, modified by the power-law factor,is used as the asymptotic form of the propagator for diffusion on fractals, e.g.on the Sierpi´nski gasket [24].The coefficients C j are defined by means of the following expressionΓ(1 − a + A s )Γ(1 − b + B s ) ( βα α ) − s ≡ G = ∞ X j =0 C j Γ( αs + α ′ + j ) . (25)They can be explicitly evaluated by a subtraction of the consecutive termsand by taking the limit s → ∞ . More precisely, C j are given by the followingrecurrence formula C j = lim s →∞ " G Γ( αs + α ′ ) − C − C αs + α ′ − . . . − C j ( αs + α ′ ) . . . ( αs + α ′ + j − ×× ( αs + α ′ )( αs + α ′ + 1) . . . ( αs + α ′ + j −
1) (26)where C = √ πα α ′ − / A / − a B b − / (27)9as been obtained from the expansion of the gamma functions by means of theStirling formula. The exponent of the stretched Gaussian, 1 /α , is connectedwith higher moments of the distribution p ( x, t ) and it cannot be determinedin the framework of the diffusion approximation.All moments of the distribution p ( x, t ) are convergent. In particular, the vari-ance, which determines the diffusion properties of the system, is given by theexpansion coefficient h in a simple way: h x i = − ∂ ∂k e p (0 , t ) = h a − ∼ t θ . (28)For θ = 0 the diffusion coefficient D = lim t →∞ h x i ( t ) / t assumes a finite valueand the diffusion is normal. The case θ = 0 means the anomalous diffusion:either the enhanced one for θ <
0, or the subdiffusion for θ >
0; the diffusioncoefficients are then ∞ or 0, respectively. The kind of diffusion depends onlyon θ and it is not sensitive on free parameters. The anomalous diffusion isfrequently encountered in physical phenomena, in particular in complex anddisordered systems [25], as well as in dynamical systems [4].On the other hand, the fractional equation (5) for µ = 2 assumes a form ofthe simple diffusion equation ∂p ( x, t ) ∂t = K ∂ [ | x | − θ p ( x, t )] ∂x (29)which can be solved exactly just by assuming the scaling form of the solution p ( x, t ) = a ( t ) f ( a ( t ) x ). The functions a ( t ) and f ( ax ) are derived by insertingto the Eq.(29) and by separation of the variables [26]. That procedure finallyyields p ( x, t ) = N ( K t ) − θ θ | x | θ exp − | x | θ K (2 + θ ) t ! . (30)Eq.(29) follows from the master equation (3) with the Gaussian λ ( x ), when oneneglects all terms higher than the second one in the Kramers-Moyal expansion.That procedure is justified if jumps are small and ν ( x ) is a smooth function[27]. We can expect that the solution of the master equation converges withtime to Eq.(30) and this convergence is fast for small | θ | . Convergence of thetails must be slow, because the contribution from large jumps is substantialfor large | x | . Indeed, we will demonstrate in the following that for large | θ | , inparticular for θ close to −
1, a very long time is required.The result (30) is useful for further improvement of the solution (18). Thecomparison of Eqs. (24) and (30) yields the conditions for the Fox function10oefficients which ensure the proper limit µ →
2: lim µ → ( B − A ) = 1 / (2 + θ )and lim µ → (1 − α ′ ) /α = θ . By inserting the coefficients which follow from thoselimiting values to the Eq.(18) and by assuming, in addition, that B = 1, weobtain a particular solution of the fractional equation (5) in the form p ( x, t ) = N a ( t ) H , , a ( t ) | x | (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (1 − − θµ + θ , µ + θ ) , ( + θ (1+ θ )2+ θ , − θ )( θ, , (1 − − θ θ , θ ) . (31)Therefore, the above solution of the fractional equation takes into account theKramers-Moyal result (30) in the limit µ → θ → The L´evy flights are characterized by the power-law tail with the exponentsmaller than 3 which implies the infinite variance. However, in physical sys-tems, which are limited in space, the variance must always be finite. Thefiniteness of the available space must be taken into account in any attemptto simulate the random walk in lattices, for example, in a model of turbu-lence which describes a transport in a Boltzmann lattice gas [28]. One canconstrain the jumping size either by taking into account that the Brownianparticle actually possesses a finite velocity, i.e. to substitute the L´evy flightsby the L´evy walks, or by introducing some cutoff in the jumping size prob-ability distribution. As a result, the variance becomes finite and, in the caseof the mutually independent jumps, the random walk probability distributionconverges to the Gaussian, according to the central limit theorem. The trun-cation of the L´evy flights can be accomplished either by a simple removingof the tail [29] or by multiplying the tail by some fast falling function. Anobvious choice in this context is the exponential exp( − γ | x | ) and the algebraicfunction | x | − β , where β ≥ − µ . The former case was considered by Koponenin Ref.[30], where an analytic expression for the characteristic function wasderived. The L´evy flights with exponential truncation serve as a model for phe-nomena in many fields, e.g. in turbulence [31], solar systems, economy. Thedistribution function of velocity and magnetic-field vector differences withinsolar wind can be reasonably fitted in this way [32]. In the framework of theeconomic research, the L´evy process is a natural model of the financial assetsflow and the fractional equations are applied to characterize the dynamicsof stock prices, where rare, non-Gaussian events are frequently encountered,in particular if the market exhibits high volatility. One can improve in this11 ,0 0,2 0,4 0,6 0,8 1,010 -3 -2 -1 p ( x ,t ) / a a|x| Fig. 1. The distributions p ( x, t ) obtained from trajectory simulations for the pow-er-law cutoff (35) with µ ′ = 5 . t = 5, t = 10, t = 50 and t = 200 (from top tobottom). The distribution predicted by Eq.(30) is also presented and it coincideswith the curve for t = 200. Other parameters are the following: θ = 1, σ = 0 . a = t − / (2+ θ ) . way the Black-Scholes model, commonly used to price the options, which isrestricted to the Gaussian distributions. In order to incorporate the finite-ness of the financial system to the fractional equations formalism by meansof the exponential truncation of the L´evy tails, models known as CGMY andKoBoL have been devised [33]. The non-Markovian fractional equation, whichresults from the CTRW model with the L´evy distributed and exponentiallytruncated jump size, has been considered in Ref.[34]. The solutions do notexhibit a typical scaling at small time but they converge, asymptotically, tothe stretched Gaussian which is predicted by the subdiffusive case of CTRWwith the Gaussian step-size distribution.In Ref. [30], the following jumping size distribution, in the form of the L´evytail multiplied by the exponential factor, has been introduced: λ ( x ) = N e − γ | x | | x | − µ − (32)and the characteristic function for the process has been evaluated. In theabove formula γ ≥ N is the normalization constant. For the symmetricprocess, the normalized Fourier transform from λ ( x ) is given by [30,34] e λ ( k ) = 4 π µ Γ( µ ) tan( πµ/ γ µ − ( k + γ ) µ/ cos( µ arctan( k/γ )] + 1 . (33)We keep only the terms of the lowest order in | k | . The expansion of the ex-12 T Fig. 2. The time T needed to reach the distribution (30), as a function of θ , for thecase of the algebraic cutoff (35) with the parameters µ ′ = 5 . σ = 0 . pression (33) produces the following result e λ ( k ) ≈ π µ Γ( µ ) γ µ − tan( πµ/ µ − µ ) k ≡ − K E k . (34)The diffusion process is then described by Eq.(29) with K = K E .On the other hand, one can apply the power-law cutoff to the L´evy tail. InRef.[35] the fractional equation of the distributed order, with the constantdiffusion coefficient, was introduced; it implies the cutoff in the form of thepower-law function with the exponent 5 − µ . Since the equation involves boththe fractional and the diffusion component, there is no simple scaling. In thispaper, we assume the jumping size distribution in the form of the modifiedL´evy tail: λ ( x ) = | x | ≤ σµ ′ σ µ ′ | x | − µ ′ − for | x | > σ, (35)where µ ′ = µ + β > e λ ( k ) = 1 − µ ′ µ ′ −
2) ( kσ ) − [Γ(1 − µ ′ ) cos( πµ ′ / | k | σ ) µ ′ + . . . . (36)In this case we have K = µ ′ σ / µ ′ −
100 200 300 400 5000,00,51,01,52,02,53,0 t Fig. 3. The shape parameter δ of the distribution tail exp( − const | ax | δ ) as a functionof time for the exponential cutoff (32). The curves correspond to the cases: θ = 1, θ = 0 and θ = − . µ = 1 . γ = 1. In the following, we evaluate the probability density distribution p ( x, t ) forboth forms of the L´evy flight cutoff, exponential and algebraic, by means of therandom walk trajectory simulations. The waiting time is sampled from Eq.(1)and the jumping size from either (32) or (35). Fig.1 presents the power-law casefor θ = 1, which corresponds to the subdiffusion. The results are comparedwith the Kramers-Moyal limiting distribution (30), which is expected to bereached at large time. We observe a rapid convergence for small and medium | x | -values, whereas the tails reach the form (30) at about t = 200. Nevertheless,the shape of the tails is always stretched exponential ∼ exp( − const | ax | δ ) andthe index δ rises with time. The speed of convergence to the distribution (30)strongly depends on the parameter θ . In Fig.2 we present the convergencetime T as a function of θ . It is relatively short only for small | θ | ; for thiscase the kernel in the master equation (3) changes weakly with | x | and thehigher terms in the Kramers-Moyal expansion soon become negligible. T risesrapidly for the negative θ : the estimation presented in the figure suggests thatthe dependence T ( θ ) is exponential, ∼ exp( − θ ), which yields T ∼ when θ approaches -1. The rapid growth of the time needed to reach convergenceof the tails for the negative θ may also be related to a specific shape of thedistribution. The tails become flat for θ <
0, and the asymptotics emergesfirst for very large | x | . On the other hand, the probability density to stay inthe origin, p (0 , t ), is then infinite: we have p ( x, t ) ∼ t − (1+ θ ) / (2+ θ ) | x | θ ( | x | ≪ θ = 0 we obtain the usual resultfor diffusion: p (0 , t ) ∼ / √ t .Application of the exponential cutoff to the L´evy tail produces similar prob-14 ,6 0,8 1,0 1,2 1,410 -3 -2 -1 p ( x ,t ) / a a|x| Fig. 4. The distributions p ( x, t ) at t = 100 and t = 200 obtained from trajectorysimulations for the exponential cutoff with µ = 1 . γ = 2. These curves coincideand assume the numerically estimated shape ∼ a | x | exp( − . a | x | ) . . The distribu-tion predicted by Eq.(30) is also presented. The parameter θ = 1 and a = t − / (2+ θ ) . ability density distributions to those presented in Fig.1 and they are alsocharacterized by the stretched Gaussian tails. However, the convergence rateof the index δ to the value 2 + θ , predicted by the solution (30), is smallerthan for the power-law truncation because the kernel in Eq.(3) is steeper inthis case. In Fig.3 we present the dependence δ ( t ) for all kinds of the diffu-sion. Initially, the exponent rises fast but then it begins to stabilize and thecurves approach the asymptotic values very slowly, especially for the superdif-fusive case of the negative θ . Discrepancies from Eq.(30) are more pronouncedfor γ = 2. This case is presented in Fig.4: though the rescaled distributionseems to be stabilized already for t = 100, its shape differs from (30) also forintermediate values of a | x | .The exponential asymptotics of the probability density distributions is guar-antied by the existence of the finite second moment. However, the time neededto reach that asymptotics can be so large [29] that it cannot be observed in anypractical realizations of the process. In fact, the convergence rate to the normaldistribution is governed by the third moment, according to the Berry-Ess´eentheorem [37] which refers to the sum of mutually independent variables, sam-pled from the same distribution. It states that if the third moment ρ is finite,then the deviation of a given distribution from the normal one is less than33 ρ/ (4 σ √ n ), where n the number of steps and σ is the standard deviation. If ρ = ∞ , the convergence to the Gaussian may be problematic in practice. Todemonstrate that the case of divergent third moment is exceptional also forour process, let us consider the power-law cutoff of the L´evy tail, Eq.(35), withthe parameter µ ′ = 2 .
5. The tail of the resulting random walk distribution,presented in Fig.5, is no longer exponential but it assumes the power-law form15 -6 -5 -4 -3 -2 -1 t=5 t=10 t=20 t=200 p ( x ,t ) / a a|x| Fig. 5. The distributions p ( x, t ), plotted in the double logarithmic scale, obtainedfrom trajectory simulations for the power-law cutoff with µ ′ = 2 . θ = 1 and σ = 0 . | x | − δ which persists to the largest times, numerically accessible. The param-eter δ is constant and well determined for small time: δ = 3 .
4. Moreover, itseems to be independent of θ . The presence of algebraic tail is not restrictedto the case µ ′ ≤
3, for which the third moment is divergent; also for slightlylarger values of µ ′ it is clearly visible. In the case µ ′ = 3 . δ = 4 . δ = 4 . µ ′ = 3 . < δ < . < δ < δ = 3 . We have solved the fractional equation which follows from the master equationfor the jumping process in the diffusion approximation of small wave numbers.The jumping size distribution has the L´evy form. The jumping rate depends onposition and then the diffusion coefficient in the fractional equation is variable.We have considered the jumping rate in the algebraic form ν ( x ) = | x | − θ , wellsuited for the diffusion on self-similar structures. The generalization to otherdependences ν ( x ) is possible [7]. It has been demonstrated that the equation16s satisfied by the scaling formula which can be expressed in terms of theFox function H , , , when one neglects higher terms in the k -expansion. In thelimit µ →
2, the solution reduces itself to the Fox function of the lower orderand it exhibits the stretched-Gaussian asymptotic form (24). The exponentof that function, 1 /α , is related to the higher moments and it cannot beuniquely determined in the diffusion approximation. The solution predicts allkinds of diffusion, both normal and anomalous, which are distinguished by theparameter θ . The diffusion equation can be solved exactly for the case µ = 2and that solution constitutes the Kramers-Moyal approximation of the masterequation. The requirement that the fractional equation solution should agreewith that result in the limits µ → t → ∞ allows us to find additionalconditions for the Fox function coefficients.In the approximation of small wave numbers, the problem of truncated L´evyflights coincides with that of the Gaussian jump sizes. We have applied twoforms of the cutoff: the exponential and the algebraic ones to study the ran-dom walk process by the Monte Carlo method. In most cases, the probabilitydensity distributions converge with time to the Kramers-Moyal result whichpredicts α = 1 / (2 + θ ). However, that convergence appears very slow if θ isfar away from 0, especially for θ <
0. If the truncation function is steep, thedistribution seems to assume the stabilized asymptotic shape which differsslightly from Eq.(30) (see Fig.4) and this conclusion may indicate a limit ofapplicability of the Kramers-Moyal approximation. In other cases, form (30)is actually reached after a long time. For smaller time, the scaled distributionis time-dependent, in disagreement with the ansatz (7). However, the distribu-tion tail is always power-law and the exponent δ depends on time very weakly(Fig.3), compared to the function a ( t ). One can expect that in an experimen-tal situation, when the observation time is finite, the distributions are unableto converge to (30) and they may reveal the values of δ smaller than 2 + θ .The convergence of the distribution to (30) becomes problematic not only forvery sharp cutoffs but, conversely, for the functions which possess a large, inparticular infinite, third moment. Numerical simulations predict in this casethe power-law asymptotics and no trace of the exponential tail could be found.The value of the exponent of power-law tail agrees with observations, e.g. forthe financial data. References [1] R. Metzler, J. Klafter, Phys. Rep. 339 (2000) 1.[2] B. J. West, W. Deering, Phys. Rep. 246 (1994) 1.[3] D. Brockmann, T. Geisel, Phys. Rev. Lett. 90 (2003) 170601.
4] G. M. Zaslavsky, Phys. Rep. 371 (2002) 461.[5] B. O’Shaughnessy, I. Procaccia, Phys. Rev. Lett. 54 (1985) 455.[6] R. Metzler, T. F. Nonnenmacher, J. Phys. A 30 (1997) 1089.[7] T. Srokowski, Phys. Rev. E 78 (2008) 031135.[8] A. A. Vedenov, Rev. Plasma Phys. 3 (1967) 229.[9] H. Fujisaka, S. Grossmann, S. Thomae, Z. Naturforsch. Teil A 40 (1985) 867.[10] A. Kami´nska, T. Srokowski, Phys. Rev. E 69 (2004) 062103.[11] T. Geisel, J. Nierwetberg, A. Zacherl, Phys. Rev. Lett. 54 (1985) 616.[12] G. Zumofen, J. Klafter, Phys. Rev. E 47 (1993) 851.[13] J. Klafter, A. Blumen, M. F. Shlesinger, Phys. Rev. A 35 (1987) 3081.[14] M. F. Shlesinger, B. J. West, J. Klafter, Phys. Rev. Lett. 58 (1987) 1100.[15] T. Srokowski, A. Kami´nska, Phys. Rev. E 74 (2006) 021103.[16] E. Barkai, Chem. Phys. 284 (2002) 13.[17] A. V. Chechkin, V. Y. Gonchar, J. Klafter, R. Metzler, L. V. Tanatarov, J.Stat. Phys. 115 (2004) 1505.[18] W. R. Schneider, in: S. Albeverio, G. Casati, D. Merlini (Eds.), StochasticProcesses in Classical and Quantum Systems, Lecture Notes in Physics, Vol.262, Springer, Berlin, 1986.[19] C. Fox, Trans. Am. Math. Soc. 98 (1961) 395.[20] A. M. Mathai, R. K. Saxena,
The H-function with Applications in Statisticsand Other Disciplines , (Wiley Eastern Ltd., New Delhi, 1978).[21] B. L. S. Braaksma, Compos. Math. 15 (1964) 239.[22] W. Wyss, J. Math. Phys. 27 (1986) 2782.[23] J. Klafter, G. Zumofen, J. Phys. Chem. 98 (1994) 7366.[24] H. E. Roman, Phys. Rev. E 51 (1995) 5422.[25] J.-P. Bouchaud, A. Georges, Phys. Rep. 195 (1990) 12.[26] Kwok Sau Fa, E. K. Lenzi, Phys. Rev. E 67 (2003) 061105.[27] N. G. van Kampen,
Stochastic Processes in Physics and Chemistry (North-Holland, Amsterdam, 1981).[28] F. Hayot, L. Wagner, Phys. Rev E 49 (1994) 470.[29] R. Mantegna, H. E. Stanley, Phys. Rev. Lett. 73 (1994) 2946.
30] I. Koponen, Phys. Rev. E 52 (1995) 1197.[31] B. Dubrulle, J.-Ph. Laval, Eur. Phys. J. B 4 (1998) 143.[32] R. Bruno, L. Sorriso-Valvo, V. Carbone, B. Bavassano, Europhys. Lett. 66(2004) 146.[33] A. Cartea, D. del-Castillo-Negrete, Physica A 374 (2007) 749.[34] A. Cartea, D. del-Castillo-Negrete, Phys. Rev. E 76 (2007) 041105.[35] I. M. Sokolov, A. V. Chechkin, J. Klafter, Physica A 336 (2004) 245.[36] M. Marseguerra, A. Zoia, Physica A 377 (2007) 1.[37] W. Feller,
An introduction to probability theory and its applications (John Wileyand Sons, New York, 1966), Vol.II.[38] X. Gabaix, P. Gopikrishnan, V. Plerou, H. E. Stanley, Nature 423 (2003) 267.[39] P. Gopikrishnan, M. Meyer, L. A. N. Amaral, H. E. Stanley, Eur. Phys. J. B 3(1998) 139.[40] V. Plerou, P. Gopikrishnan, L. A. Nunes Amaral, M. Meyer, H. E. Stanley,Phys. Rev. E 60 (1999) 6519.[41] H. E. Stanley, Physica A 318 (2003) 279.(John Wileyand Sons, New York, 1966), Vol.II.[38] X. Gabaix, P. Gopikrishnan, V. Plerou, H. E. Stanley, Nature 423 (2003) 267.[39] P. Gopikrishnan, M. Meyer, L. A. N. Amaral, H. E. Stanley, Eur. Phys. J. B 3(1998) 139.[40] V. Plerou, P. Gopikrishnan, L. A. Nunes Amaral, M. Meyer, H. E. Stanley,Phys. Rev. E 60 (1999) 6519.[41] H. E. Stanley, Physica A 318 (2003) 279.