Exact and explicit probability densities for one-sided L´evy stable distributions K. A. Penson1,∗ and K. G´orska1,2,† 1Laboratoire de Physique Th´eorique de la Mati`ere Condens´ee (LPTMC), Universit´e Pierre et Marie Curie, CNRS UMR 7600, Tour 13 - 5i`eme ´et., B.C. 121, 4 pl. Jussieu, F 75252 Paris Cedex 05, France 2Nicolaus Copernicus University, Institute of Physics, ul. Grudzia¸dzka 5/7, PL 87-100 Torun´, Poland Westudyfunctionsg (x)which are one-sided,heavy-tailedL´evy stableprobability distributions 1 α of index α, 0 < α < 1, of fundamental importance in random systems, for anomalous diffusion 1 0 and fractional kinetics. We furnish exact and explicit expressions for gα(x), 0 ≤ x < ∞, for all 2 α = l/k < 1, with k and l positive integers. We reproduce all the known results given by k ≤ 4 and present many new exact solutions for k > 4, all expressed in terms of known functions. This n will allow a ’fine-tuning’of α in order to adapt g (x) to a given experimental situation. α a J PACSnumbers: 05.40.Fb,05.40.-a,02.50.Ng 5 ] Theoreticaldescriptionofmanycollectivephysicalsys- been hampered for subjective and objective reasons [20, h c temswhichincludespecialsortofdisorderorrandomness 21]. The subjective ones include a certain reticence to e often requires a radical departure from classical diffu- usedistributionswithbothmeanandvariancediverging. m sive behaviour. On the probabilistic level this signifies Themainobjectivereasonisalackofknowledgeofg (x) α - the appearance of distributions with non-conventional formostvaluesofα. Theexistinginterpolationformulae t a characteristics, like diverging mean and variance along [15] appear to be cumbersome to use. t with all integer moments different from zeroth one. It seems that obtaining explicit g (x) for arbitrary s α . In this context the prominent rˆole plays the discovery 0 < α < 1 constitutes a true challenge: only for a t a of particular distributions with these properties, called limited number of values of α, i.e. α = 1/2 (s. g (x) 1/2 m now L´evy stable laws [1], whose generic example is above), 1/4 [22], 1/3 [23], 2/3 [24], and 3/4 [23] the ex- d- g”1s/t2a(bxle)”=mea(2n√s πhexr3e/2t)h−a1tetxhpe(−pr1o/d4uxc)t, xof≥ch0a;rathcteerwisotridc parlibciittraforrymαs o[4f,g9α](xis)oanrelyknofowanl.imTihteedfoursmeaalssoitlurteiqonuirfoers n functions (cf) of two such laws is a cf of another law seriesorasymptoticexpansions,whichmaybecomeprob- o of the same type [1]. The general distribution of that lematic, especially for small α. c [ type gα(x) canbe showntopossessthe cforthe Laplace The objective of this work is to present an universal transform of the form [2–4]: formula for g (x), α = l/k, with k > l positive integers, 3 α v ∞ which is exact and explicit. It reproduces all the known 3 e−pxgα(x)dx = e−pα, p>0, 0<α<1, (1) cases enumerated above, and yields an infinity of new Z 9 0 solutions for k >4, of which we quote, for the first time, 1 several instances. which is the well-known Kohlrausch-Williams-Watts 0 function [5] or stretched exponential. Several indepen- Eq. (1) for α=l/k can be inverted giving . 7 dent proofs can be given that g (x) obeying Eq. (1) is α 0 positive [2, 3, 6]. √kl 1 ll 0 g (x) = Gk,0 ∆(l,0) , (2) 1 The functions gα(x) are ubiquitous in many fields of l/k (2π)(k−l)/2 x l,k (cid:18)kkxl (cid:12)∆(k,0)(cid:19) condensedandsoftmatter physics [7–9], geophysics[10], (cid:12) v: meteorology [11], economics [12], fractional kinetics [13, (cid:12) i valid for all x 0, where Gm,n z (ap) is the Meijer G X 14], etc. For instance, the value α = 1/4 is thought to ≥ p,q (cid:16) |(bq)(cid:17) describe mechanical and dielectric properties of glassy function [29, 30] and ∆(k,a) = a,a+1,...,a+k−1 is a r k k k a polymers[15]. Itisalsoconfirmedthatthesamevalueof special list of k elements. Eq. (2) is listed without proof αisrelevantforastatisticaldescriptionofsubrecoillaser as a special case for ν =0 and a=1 of formula 2.2.1.19, cooling[16,17]. Ingeneral,numerousphenomena falling in vol. 5 of [30]. The detailed demonstration of Eq. (2) in the class of subdiffusion [18] call for gα(x), α < 1, with a combined use of Laplace and Mellin transforms in their theoretical description. On theoretical side, the will be given elsewhere. It turns out that the r. h. s. of L´evy stable distributions are essential tools in the study Eq.(2)isafinitesumofk 1generalizedhypergeometric − of random maps and resulting combinatorial structures functions of type F (ap) z [30]: [19]. TheactualuseofL´evystabletypedistributionshas p q(cid:16)(bq)(cid:12) (cid:17) (cid:12) k−1b (k,l) ll g (x)= j F 1,∆(l,1+jl/k) ( 1)k−l , ∗ [email protected] l/k Xj=1 x1+j kl l+1 k(cid:18) ∆(k,j+1) (cid:12)(cid:12) − kkxl(cid:19) † [email protected] (cid:12) (3) 2 where b (k,l) are numerical coefficients given by of k and l. The application of Gauss-Legendre multipli- j cation formula to both of these gamma functions allows bj(k,l)= kjl(j2l/πk)(√k−kll)/2 hQij=−11Γ(cid:0)i−kl−j(cid:1)1iΓhQi ik=−j1j+1Γ(cid:0)i−kj(cid:1)i, tohfeTcohideeffienactdiievfinactnasttaibogjn(ekoo,ffl)opoFufrqE’ssqoi.lnu(t4Ei)oq.n.,(E3)qsa.n(d3)thaendex(t4ra)cotvioenr i=1 l − k Q (cid:0) (cid:1) (4) Eq. (7) is clearly seen in practice, in conjunction with where Γ(y) is Euler’s gamma function. Eq. (3) is the ex- the use of computer algebra systems [27]. Since in re- act implementation of the program outlined in the fun- cent versions of these systems the hypergeometric func- damental work of Scher and Montroll [23], in which it tions F as well as the Meijer G function are fully im- p q was conjectured that g (x) can be expressed in terms plemented, their use permits high-precisioncalculations. l/k of F ’s. However in [23] actually only one new instance Forreader’sconveniencewe givein [28]the Maple(cid:13)R syn- p q g (x)waswrittendown. OurformulaEq.(3),afterap- tax for g (x), see Eq. (2) above. Our experience indi- 3/4 l/k propriatereductions in F ’s,see below,givesallexactly cates that for small α our results for small x are more p q known cases mentioned above [20–23], with g (x) [22] practical to use than the x 0 asymptotics given in 1/4 → [3]. The reason is that there the region of applicability g (x)= b1(4,1) F − −1 + b2(4,1) (5) of Mikusin´ski’s asymptotic expansion [3] shrinks to ex- 1/4 x5/4 0 2(cid:18)1/2,3/4(cid:12)44x(cid:19) x3/2 ceedingly small values of x. For example for α = 1/20, (cid:12) g (0)=0,butahugepeaking (x)appearsalready (cid:12) 1/20 1/20 × 0F2(cid:18)3/−4,5/4(cid:12)4−41x(cid:19)+ b3x(74/,41)0F2(cid:18)5/−4,3/2(cid:12)4−41x(cid:19) axtinxt∼hi1s0r−eg14io.nI.nIcnotnhteraospt,poousirtefolrimmuitlafoerwαor>k1fintehefoHruamny- (cid:12) (cid:12) bertexpansionEq.(7)isslowlyconvergentforx<αbut (cid:12) (cid:12) and offers an unlimited number of new solutions for approximation [3] works well as then g (x) is very close α gl/k(x), k>4, e.g.: to zero in a considerable region near x=0 (e.g. already for α = 5/6 the function g (x) is practically equal to 5/6 gp/5(x)=Xj=41 xbj1(+5j,pp/5) p+1F5(cid:18)1∆,∆(5(,pj+,11+)jp/5)(cid:12)(cid:12)55pxpp(cid:19), (6) IszImeIraoolnluxpFictgao.n3xa.l≈so0b.3e5s)e.enSufochr αa=pr4a/ct5i,cacollmypflaartertehgeiocnurfvoer (cid:12) p = 1,...,4, see Table I for coefficients in Eqs. (5) and InFig.1wecomparethreedistributionsforl/k=1/2, (6), etc. 1/3 and 1/4. The salient feature for l/k = 1/4 is the appearanceofasharpmaximumforverysmallx sothat The symbol ∆(k,a) in Eq. (3) permits one to en- thesethreecurvescanbebarelyshownonthesamescale. code all the possible cases of k and l in a single for- Analogously, for l/k = 1/5 the maximum of g (x) ap- mula. However we draw attention to the fact that can- 1/5 pears at x (1/5) 0.0002 and the value g (x ) 25. cellations will appear there due to the obvious identity 0 1/5 0 ≈ ≈ F (ap),(αr) x = F (ap) x , where (α ) is an For x < x0(1/5) the values of g1/5(x) are very close to p+r q+r(cid:16)(bq),(αr)(cid:12) (cid:17) p q(cid:16)(bq)(cid:12) (cid:17) r zero. As already mentioned above, for smaller values arbitrary sequenc(cid:12)e of r parameter(cid:12)s not equal to zero or of l/k this type of behaviour is even more pronounced (cid:12) (cid:12) to negative integers. Thus Eq. (5) is a sum of three and it explains a posteriori the difficulties encountered 0F2 functions, and likewise g3/4(x), which is not spec- in devising approximations valid for small l/k and small ified here, will be a sum of three 2F2 functions, neatly x[15,24]. InFig.2we presentthe comparisonofseveral confirming Eq.(C10)of [23], etc. In this manner for any distributions for values l/k 1/2. Here the ’sharpen- l/k,aclosedformofgl/k(x)canbeobtainedfromEq.(3). ing’ of the distributions, as l/≈k goes from 1/2 to smaller However,only for k 3 it can be written down in terms values, is very clearlyvisible but is less dramatic than in ≤ of standard special functions [5, 20, 21, 23, 24]. Fig.1. WepresentinFig.3thenewdistributionsg (x) p/5 A heuristic indication of how Eq. (3) comes about can given by Eq. (6) for p=2,3 and 4. be obtained from the series representation for gα(x) de- All these probability distributions share the following rived by Humbert [25], discussed for example by Hughes features: (a) g (x) 0, for x 0, where they present α [26] and used in [19]. The series → →−2+α −α an essential singularity x2(1−α) exp[ A(α)x1−α], A(α) > 0 [3]; (b) g ∼(x) B(α)x−−(1+α), for 1 ∞ ( 1)j+1 α → g (x)= − Γ(1+αj) sin(παj) (7) x , B(α) > 0, indicating heavy-tailed asymp- α π Xj=1 j!x1+αj tot→ics ∞for large x; (c) all their fractional moments M (µ) = ∞xµg (x)dx = Γ( µ/α)/[αΓ( µ)], for α 0 α − − is a convergent expansion valid for all 0 < α < 1 and real µ, <R µ < α, including Mα(0) = 1, are finite, −∞ x>0. For α=l/k the decomposition of the summation and are infinite otherwise; (d) g (x) are unimodal with α index j modulo k yields an equivalent representation of the maximum at x (α), and x (α) 0 as α 0. 0 0 → → g (x) as a sum of k 1 infinite series. The structure The distributions g (x) constitute basic ingredients of l/k α − of the coefficients in each of these series involve gamma all theories of anomalous diffusion where they are em- functionsratiosoftypeΓ(li+θ)/Γ(ki+φ),whereiisthe ployed to produce solutions P (x,t) in space-time do- α new summation index and θ and φ are simple functions main of various forms of the Fokker-Planck equations 3 j 1 2 3 4 bj(4,1) 4Γ1(3) 4−√1π √21Γ6π(43) — 4 √5Γ(1) √5Γ(2) √5Γ(3) √5Γ(4) bj(5,1) 20πB5 −20πA5 40πA5 −120πB5 √522/5Γ(1) √524/5Γ(2) √521/5Γ(3) √523/5Γ(4) bj(5,2) 10√·πΓ(3)5B −10√·πΓ(1)5A −100√·πΓ(9)5A 100√·πΓ(75)B 10 10 10 10 3√531/10Γ(1) √537/10Γ(2) 3√533/10Γ(3) 7√539/10Γ(4) bj(5,3) 10Γ(·2)Γ(75)B 50Γ(·4)Γ(145)A −25Γ( 1· )Γ(11)5A 7−50Γ(·8)Γ(13)5B 15 15 15 15 15 15 15 15 421/105−1/2√πΓ(1) 627/105−3/2√πΓ(2) 1423/105−5/2√πΓ(3) 1129/105−7/2√πΓ(4) bj(5,4) Γ·(3)Γ(1)Γ(11)5B Γ·(1)Γ(7)Γ(17)5A Γ·(9)Γ(3)Γ(13)A5 Γ·(7)Γ(9)Γ(19)B5 10 20 20 10 20 20 10 20 20 10 20 20 TABLE I. Coefficients of Eqs. (5) and (6); A=sin(π/5) and B =sin(2π/5). along with their fractional generalizations [6, 22]. For We believe that the exact forms of g (x) obtained α instance in [22] P (x,t) is given as a convolution (called in this work, along with their asymptotics for x α there inverse L´evy transform) of d g (t/s1/α) with and exact values of fractional moments, constitu→te∞a ds − α P1(x,t) being a normalized solutio(cid:2)n of the ord(cid:3)inary solid basis to extract a value of α best suited for an Fokker-Planckequation, see Eq. (1) in [22]. The explicit experimental situation at hand. Once it has been done, forms of g (x) presented here will permit further devel- such description can be further ’fine-tuned’ by choosing α opment of this ambitious approach. values of k and l which would optimise the choice of α. Theavailabilityofg (x)makesitpossibletofullyde- Wehope thatthis approachwillproveusefulinpractical l/k scribe the long tail distributions of carrier transit times applications. in amorphous materials such as As Se and TNF-PVK. 2 3 In fact in classic work [23] the measured values of α for We thank Professor E. Barkai for kindly informing us these two materials were α = 0.45 and α = 0.8 respec- about his results obtained in [22]. tively, compare Fig. 6 of [23]. These values were for a The authors acknowledge support from Agence Na- longtime intractabletheoretically. Fromnowon,setting tionale de la Recherche (Paris, France) under Program α = 9/20 and α = 4/5 in our Eqs. (2) and (3) directly No. ANR-08-BLAN-0243-2. provides the sought for framework for interpretation of these data. The appropriate distributions are presented as the curve II in Fig. (2), (α = 9/20) and the curve III in Fig. (3), (α=4/5). [1] J.-P. Kahane, in L´evy Flights and Related Topics in (1995). Physics, (Lecture Notes in Physics, vol. 450), edited [13] I. M. Sokolov, J. Klafter, and A. Blumen, Phys. Today by M.F. Shlesinger, G.M. Zaslavsky, and U. Frisch, 55, 48 (2002). (Springer,Berlin, 1995). [14] A. V. Chechkin, V. Yu. Gonchar, R. Gorenflo, N. Kora- [2] H.Pollard, Bull. Amer. Math. Soc. 52, 908 (1946). bel, and I. M. Sokolov, Phys.Rev. E78, 021111 (2008). [3] J. Mikusin´ski, Studia Math. 18, 191 (1959). [15] J. T. Bendler, J. Stat. Phys.36, 625 (1984). [4] W. R. Schneider, in Stochastic Processes in Classical [16] F. Bardou, J.-P. Bouchaud, O. Emile, A. Aspect, and and Quantum Systems (Lecture Notes in Physics, vol. C. Cohen-Tannoudji, Phys. Rev.Lett. 72, 203 (1994). 262), edited by S. Albeverio, G. Casati, and D. Merlini [17] F. Bardou, J.-P. Bouchaud, A. Aspect, and C. Cohen- (Springer,Berlin, 1986). Tannoudji, L´evy Statistics and Laser Cooling (Cam- [5] R. S. Anderssen, S. A. Husain, and R. J. Loy, ANZIAM bridge University Press, 2002). J. 45, C800 (2004). [18] T. Koren, J. Klafter, and M. Magdziarz, Phys. Rev. E [6] I.M. Sokolov, Phys.Rev.E 63, 011104 (2000). 76, 031129 (2007). [7] R.Metzler and J. Klafter, J. Phys.A 37, R161 (2004). [19] P. Flajolet and R. Sedgewick, Analytic Combinatorics [8] P.-G. deGennes, Macromol. 35, 3785 (2002). (Cambridge University Press, 2009). [9] R.Hilfer, Phys. Rev.E 65, 061510 (2002). [20] V. V. Ucha˘ikin, Zh. E´ksp. Teor. Fiz. 115, 2113 (1999) [10] O.Sottolongo-Costa, J. C.Antoranz,A.Posadas, F.Vi- [Sov. Phys. JETP 88, 1155 (1999)]. dal, and A. Vazquez, Geophys. Res. Lett. 27, 1965 [21] V. V. Ucha˘ikin, Usp. Fiz. Nauk 173, 847 (2003) [Sov. (2000). Phys. Usp.46, 821 (2003)]. [11] M. Lagha and M. Bensebti, hal-00194153 (2007). [22] E. Barkai, Phys. Rev.E 63, 046118 (2001). [12] R. N. Mantegna and H. E. Stanley, Nature 376, 46 4 FIG.1. Comparison ofg (x): thecurvesI,IIandIIIcorre- FIG. 2. Comparison of g (x): the curves I, II, III, IV and l/k l/k spond to l/k=1/4 (s. Eq. (5)),1/3 and 1/2, respectively. V correspond to l/k =2/5,9/20,1/2,11/20 and 3/5, respec- tively. Calculations were performed using Eqs. (3) and (4). [23] H. Scher and E. W. Montroll, Phys. Rev. B 12, 2455 (1975). [24] E.W.MontrollandJ.T.Bendler,J.Stat.Phys.34,129 (1984). [25] P. Humbert,Bull. Soc. Math. Fr. 69, 121 (1945). [26] B. D. Hughes, Random Walks and Random Environ- ments, vol. 1 (Clarendon Press, Oxford,1995). [27] We havemade extensiveuseof Maple(cid:13)R in thiswork. [28] Here is the Maple(cid:13)R procedure LevyDist(k, l, x) used FIG. 3. Comparison of g (x): the curves I, II and III cor- l/k respondtol/k=2/5,3/5and 4/5, respectively. Calculations were performed using Eq. (6). to calculate Eq. (2): LevyDist := proc(k, l, x) simplify(convert( sqrt(k * l) * MeijerG([[],[seq(j1/l, j1 = 0..l-1)]] , [[seq(j2/k, j2 = 0..k-1)],[]], l^l /(k^k * x^l)) /(x*(2*Pi)^((k-l)/2)), StandardFunctions)); end; .