# TheM-Wright Function in Time-Fractional Diffusion Processes: A Tutorial Survey

Authors: Francesco Mainardi; Antonio Mura; Gianni Pagnini
Journal: International Journal of Differential Equations
 (2010)

---

## Abstract

In the present review we survey the properties of a transcendental function of the Wright type, nowadays known asM-Wright function, entering as a probability density in a
relevant class of self-similar stochastic processes that we generally refer to as time-fractional
diffusion processes. Indeed, the master equations governing these processes
generalize the standard diffusion equation by means of time-integral operators interpreted
as derivatives of fractional order. When these generalized diffusion processes are properly
characterized with stationary increments, the M-Wright function is shown to play the
same key role as the Gaussian density in the standard and fractional Brownian motions.
Furthermore, these processes provide stochastic models suitable for describing phenomena
of anomalous diffusion of both slow and fast types.


---

## Body

## 1. Introduction

By time-fractional diffusion processes, we mean certain diffusion-like phenomena governed by master equations containing fractional derivatives in time whose fundamental solution can be interpreted as a probability density function (pdf) in space evolving in time. It is well known that, for the most elementary diffusion process, the Brownian motion, the master equation, is the standard linear diffusion equation whose fundamental solution is the Gaussian density with a spatial variance growing linearly in time. In such case we speak about normal diffusion, reserving the term anomalous diffusion when the variance grows differently. A number of stochastic models for explaining anomalous diffusion have been introduced in literature; among them we like to quote the fractional Brownian motion; see, for example, [1, 2], the Continuous Time Random Walk; see, for example, [3–6], the Lévy flights; see, for example, [7], the Schneider grey Brownian motion; see [8, 9], and, more generally, random walk models based on evolution equations of single and distributed fractional order in time and/or space; see, for example, [10–18].In this survey paper we focus our attention on modifications of the standard diffusion equation, where the time can be stretched by a power law (t→tα, 0<α<2) and the first-order time derivative can be replaced by a derivative of noninteger order β (0<β≤1). In these cases of generalized diffusion processes the corresponding fundamental solution still keeps the meaning of a spatial pdf evolving in time and is expressed in terms of a special function of the Wright type that reduces to the Gaussian when β=1. This transcendental function, nowadays known as M-Wright function, will be shown to play a fundamental role for a general class of self-similar stochastic processes with stationary increments, which provide stochastic models for anomalous diffusion, as recently shown by Mura et al. [19–22].In Section2 we provide the reader with the essential notions and notations concerning the integral transforms and fractional calculus, which are necessary in the rest of the paper. In Section 3 we introduce in the complex plane ℂ the series and integral representations of the general Wright function denoted by Wλ,μ(z) and of the two related auxiliary functions Fν(z), Mν(z), which depend on a single parameter. In Section 4 we consider our auxiliary functions in real domain pointing out their main properties involving their integrals and their asymptotic representations. Mostly, we restrict our attention to the second auxiliary function, which we call M-Wright function, when its variable is in ℝ+ or in all of ℝ but extended in symmetric way. We derive a fundamental formula for the absolute moments of this function in ℝ+, which allows us to obtain its Laplace and Fourier transforms. In Section 5 we consider some types of generalized diffusion equations containing time partial derivatives of fractional order and we express their fundamental solutions in terms of the M-Wright functions evolving in time with a given self-similarity law. In Section 6 we stress how the M-Wright function emerges as a natural generalization of the Gaussian probability density for a class of self-similar stochastic processes with stationary increments, depending on two parameters (α,β). These processes are defined in a unique way by requiring the determination of any multipoint probability distribution and include the well-known standard and fractional Brownian motion. We refer to this class as the generalized grey Brownian motion (ggBm), because it generalizes the grey Brownian motion (gBm) introduced by Schneider [8, 9]. Finally, a short concluding discussion is drawn. In Appendix A we derive the fundamental solution of the time-fractional diffusion equation. In Appendix B we outline the relevance of the M-Wright function in time-fractional drift processes entering as subordinators in time-fractional diffusion.


## 2. Notions and Notations

### 2.1. Integral Transforms Pairs

In our analysis we will make extensive use of integral transforms of Laplace, Fourier, and Mellin types; so we first introduce our notation for the corresponding transform pairs. We do not point out the conditions of validity and the main rules, since they are given in any textbook on advanced mathematics.Let(2.1)f̃(s)=ℒ{f(r);r→s}=∫0∞e-srf(r)dr
be the Laplace transform of a sufficiently well-behaved function f(r) with r∈ℝ+, s∈ℂ, and let(2.2)f(r)=ℒ-1{f̃(s);s→r}=12πi∫Bre+srf̃(s)ds
be the inverse Laplace transform, where Br denotes the so-called Bromwich path, a straight line parallel to the imaginary axis in the complex s-plane. Denoting by↔ℒthe justaposition of the original function f(r) with its Laplace transform f̃(s), the Laplace transform pair reads(2.3)f(r)↔ℒf̃(s).Let(2.4)f̂(κ)=ℱ{f(x);x→κ}=∫-∞+∞e+iκxf(x)dx
be the Fourier transform of a sufficiently well-behaved function f(x) with x∈ℝ, κ∈ℝ, and let(2.5)f(x)=ℱ-1{f̂(κ);κ→x}=12π∫-∞+∞e-iκxf̂(κ)dκ
be the inverse Fourier transform. Denoting by↔ℱthe justaposition of the original function f(x) with its Fourier transform f̂(κ), the Fourier transform pair reads(2.6)f(x)↔ℱf̂(κ).Let(2.7)f*(s)=ℳ{f(r);r→s}=∫0∞rs-1f(r)dr
be the Mellin transform of a sufficiently well-behaved function f(r) with r∈ℝ+, s∈ℂ, and let(2.8)f(r)=ℳ-1{f*(s);s→r}=12πi∫Brr-sf*(s)ds
be the inverse Mellin transform. Denoting by↔ℳthe justaposition of the original function f(r) with its Mellin transform f*(s), the Mellin transform pair reads(2.9)f(r)↔ℳf*(s).


### 2.2. Essentials of Fractional Calculus with Support inℝ+

Fractional calculus is the branch of mathematical analysis that deals with pseudodifferential operators that extend the standard notions of integrals and derivatives to any positive noninteger order. The term fractional is kept only for historical reasons. Let us restrict our attention to sufficiently well-behaved functionsf(t) with support in ℝ+. Two main approaches exist in the literature of fractional calculus to define the operator of derivative of noninteger order for these functions, referred to Riemann-Liouville and to Caputo. Both approaches are related to the so-called Riemann-Liouville fractional integral defined for any order μ>0 as(2.10)Jtμf(t):=1Γ(μ)∫0t(t-τ)μ-1f(τ)dτ.
We note the convention Jt0=I (Identity) and the semigroup property(2.11)JtμJtν=JtνJtμ=Jtμ+ν,μ≥0,ν≥0.
The fractional derivative of order μ>0 in the Riemann-Liouville sense is defined as the operatorDtμ which is the left inverse of the Riemann-Liouville integral of order μ (in analogy with the ordinary derivative), that is,(2.12)DtμJtμ=I,μ>0.
If m denotes the positive integer such that m-1<μ≤m, we recognize, from (2.11) and (2.12),Dtμf(t):=DtmJtm-μf(t); hence,(2.13)Dtμf(t)={dmdtm[1Γ(m-μ)∫0tf(τ)dτ(t-τ)μ+1-m],m-1<μ<m,dmdtmf(t),μ=m.
For completeness we define Dt0=I.On the other hand, the fractional derivative of orderμ>0 in the Caputo sense is defined as the operator *Dtμ such that *Dtμf(t):=Jtm-μDtmf(t); hence,(2.14)*Dtμf(t)={1Γ(m-μ)∫0tf(m)(τ)dτ(t-τ)μ+1-m,m-1<μ<m,dmdtmf(t),μ=m.
We note the different behavior of the two derivatives in the limit μ→(m-1)+. In fact,(2.15)μ→(m-1)+{Dtμf(t)→DtmJt1f(t)=Dt(m-1)f(t),*Dtμf(t)→Jt1Dtmf(t)=Dt(m-1)f(t)-Dt(m-1)f(0+),
where the limit for t→0+ is taken after the operation of derivation.Furthermore, recalling the Riemann-Liouville fractional integral and derivative of the power law fort>0,(2.16)Jtμtγ=Γ(γ+1)Γ(γ+1+μ)tγ+μ,Dtμtγ=Γ(γ+1)Γ(γ+1-μ)tγ-μ,μ>0,γ>-1,
we find the relationship between the two types of fractional derivative(2.17)Dμ[f(t)-∑k=0m-1tkk!f(k)(0+)]=*Dtμf(t).
We note that the Caputo definition for the fractional derivative incorporates the initial values of the function and of its integer derivatives of lower order. The subtraction of the Taylor polynomial of degree m-1 at t=0+ from f(t) is a sort of regularization of the fractional derivative. In particular, according to this definition, the relevant property that the derivative of a constant is zero is preserved for the fractional derivative.Let us finally point out the rules for the Laplace transform with respect to the fractional integral and the two fractional derivatives. These rules are expected to properly generalize the well-known rules for standard integrals and derivatives.For the Riemann-Liouville fractional integral, we have(2.18)ℒ{Jtμf(t);t→s}=f̃(s)sμ,μ≥0.
For the Caputo fractional derivative, we consequently get(2.19)ℒ{*Dtμf(t);t→s}=sμf̃(s)-∑k=0m-1sμ-1-kf(k)(0+),m-1<μ≤m,
where f(k)(0+):=limt→0+f(k)(t). The corresponding rule for the Riemann-Liouville fractional derivative is more cumbersome and it reads(2.20)ℒ{Dtμf(t);t→s}=sμf̃(s)-∑k=0m-1[DtkJt(m-μ)]f(0+)sm-1-k,m-1<μ≤m,
where the limit for t→0+ is understood to be taken after the operations of fractional integration and derivation. As soon as all the limiting values f(k)(0+) are finite and m-1<μ<m, formula (2.20) for the Riemann-Liouville derivative is simplified into(2.21)ℒ{Dtμf(t);t→s}=sμf̃(s),m-1<μ<m.
In the special case f(k)(0+)=0 for k=0,1,m-1, we recover the identity between the two fractional derivatives. The Laplace transform rule (2.19) was practically the key result of Caputo [23, 24] in defining his generalized derivative in the late sixties. The two fractional derivatives have been well discussed in the 1997 survey paper by Gorenflo and Mainardi [25]; see also [26], and in the 1999 book by Podlubny [27]. In these references the authors have pointed out their preference for the Caputo derivative in physical applications where initial conditions are usually expressed in terms of finite derivatives of integer order.For further reading on the theory and applications of fractional calculus, we recommend the recent treatise by Kilbas et al. [28].


## 2.1. Integral Transforms Pairs

In our analysis we will make extensive use of integral transforms of Laplace, Fourier, and Mellin types; so we first introduce our notation for the corresponding transform pairs. We do not point out the conditions of validity and the main rules, since they are given in any textbook on advanced mathematics.Let(2.1)f̃(s)=ℒ{f(r);r→s}=∫0∞e-srf(r)dr
be the Laplace transform of a sufficiently well-behaved function f(r) with r∈ℝ+, s∈ℂ, and let(2.2)f(r)=ℒ-1{f̃(s);s→r}=12πi∫Bre+srf̃(s)ds
be the inverse Laplace transform, where Br denotes the so-called Bromwich path, a straight line parallel to the imaginary axis in the complex s-plane. Denoting by↔ℒthe justaposition of the original function f(r) with its Laplace transform f̃(s), the Laplace transform pair reads(2.3)f(r)↔ℒf̃(s).Let(2.4)f̂(κ)=ℱ{f(x);x→κ}=∫-∞+∞e+iκxf(x)dx
be the Fourier transform of a sufficiently well-behaved function f(x) with x∈ℝ, κ∈ℝ, and let(2.5)f(x)=ℱ-1{f̂(κ);κ→x}=12π∫-∞+∞e-iκxf̂(κ)dκ
be the inverse Fourier transform. Denoting by↔ℱthe justaposition of the original function f(x) with its Fourier transform f̂(κ), the Fourier transform pair reads(2.6)f(x)↔ℱf̂(κ).Let(2.7)f*(s)=ℳ{f(r);r→s}=∫0∞rs-1f(r)dr
be the Mellin transform of a sufficiently well-behaved function f(r) with r∈ℝ+, s∈ℂ, and let(2.8)f(r)=ℳ-1{f*(s);s→r}=12πi∫Brr-sf*(s)ds
be the inverse Mellin transform. Denoting by↔ℳthe justaposition of the original function f(r) with its Mellin transform f*(s), the Mellin transform pair reads(2.9)f(r)↔ℳf*(s).


## 2.2. Essentials of Fractional Calculus with Support inℝ+

Fractional calculus is the branch of mathematical analysis that deals with pseudodifferential operators that extend the standard notions of integrals and derivatives to any positive noninteger order. The term fractional is kept only for historical reasons. Let us restrict our attention to sufficiently well-behaved functionsf(t) with support in ℝ+. Two main approaches exist in the literature of fractional calculus to define the operator of derivative of noninteger order for these functions, referred to Riemann-Liouville and to Caputo. Both approaches are related to the so-called Riemann-Liouville fractional integral defined for any order μ>0 as(2.10)Jtμf(t):=1Γ(μ)∫0t(t-τ)μ-1f(τ)dτ.
We note the convention Jt0=I (Identity) and the semigroup property(2.11)JtμJtν=JtνJtμ=Jtμ+ν,μ≥0,ν≥0.
The fractional derivative of order μ>0 in the Riemann-Liouville sense is defined as the operatorDtμ which is the left inverse of the Riemann-Liouville integral of order μ (in analogy with the ordinary derivative), that is,(2.12)DtμJtμ=I,μ>0.
If m denotes the positive integer such that m-1<μ≤m, we recognize, from (2.11) and (2.12),Dtμf(t):=DtmJtm-μf(t); hence,(2.13)Dtμf(t)={dmdtm[1Γ(m-μ)∫0tf(τ)dτ(t-τ)μ+1-m],m-1<μ<m,dmdtmf(t),μ=m.
For completeness we define Dt0=I.On the other hand, the fractional derivative of orderμ>0 in the Caputo sense is defined as the operator *Dtμ such that *Dtμf(t):=Jtm-μDtmf(t); hence,(2.14)*Dtμf(t)={1Γ(m-μ)∫0tf(m)(τ)dτ(t-τ)μ+1-m,m-1<μ<m,dmdtmf(t),μ=m.
We note the different behavior of the two derivatives in the limit μ→(m-1)+. In fact,(2.15)μ→(m-1)+{Dtμf(t)→DtmJt1f(t)=Dt(m-1)f(t),*Dtμf(t)→Jt1Dtmf(t)=Dt(m-1)f(t)-Dt(m-1)f(0+),
where the limit for t→0+ is taken after the operation of derivation.Furthermore, recalling the Riemann-Liouville fractional integral and derivative of the power law fort>0,(2.16)Jtμtγ=Γ(γ+1)Γ(γ+1+μ)tγ+μ,Dtμtγ=Γ(γ+1)Γ(γ+1-μ)tγ-μ,μ>0,γ>-1,
we find the relationship between the two types of fractional derivative(2.17)Dμ[f(t)-∑k=0m-1tkk!f(k)(0+)]=*Dtμf(t).
We note that the Caputo definition for the fractional derivative incorporates the initial values of the function and of its integer derivatives of lower order. The subtraction of the Taylor polynomial of degree m-1 at t=0+ from f(t) is a sort of regularization of the fractional derivative. In particular, according to this definition, the relevant property that the derivative of a constant is zero is preserved for the fractional derivative.Let us finally point out the rules for the Laplace transform with respect to the fractional integral and the two fractional derivatives. These rules are expected to properly generalize the well-known rules for standard integrals and derivatives.For the Riemann-Liouville fractional integral, we have(2.18)ℒ{Jtμf(t);t→s}=f̃(s)sμ,μ≥0.
For the Caputo fractional derivative, we consequently get(2.19)ℒ{*Dtμf(t);t→s}=sμf̃(s)-∑k=0m-1sμ-1-kf(k)(0+),m-1<μ≤m,
where f(k)(0+):=limt→0+f(k)(t). The corresponding rule for the Riemann-Liouville fractional derivative is more cumbersome and it reads(2.20)ℒ{Dtμf(t);t→s}=sμf̃(s)-∑k=0m-1[DtkJt(m-μ)]f(0+)sm-1-k,m-1<μ≤m,
where the limit for t→0+ is understood to be taken after the operations of fractional integration and derivation. As soon as all the limiting values f(k)(0+) are finite and m-1<μ<m, formula (2.20) for the Riemann-Liouville derivative is simplified into(2.21)ℒ{Dtμf(t);t→s}=sμf̃(s),m-1<μ<m.
In the special case f(k)(0+)=0 for k=0,1,m-1, we recover the identity between the two fractional derivatives. The Laplace transform rule (2.19) was practically the key result of Caputo [23, 24] in defining his generalized derivative in the late sixties. The two fractional derivatives have been well discussed in the 1997 survey paper by Gorenflo and Mainardi [25]; see also [26], and in the 1999 book by Podlubny [27]. In these references the authors have pointed out their preference for the Caputo derivative in physical applications where initial conditions are usually expressed in terms of finite derivatives of integer order.For further reading on the theory and applications of fractional calculus, we recommend the recent treatise by Kilbas et al. [28].


## 3. The Functions of the Wright Type

### 3.1. The General Wright Function

The Wright function, which we denote byWλ,μ,(z), is so named in honor of E. Maitland Wright, the eminent British mathematician, who introduced and investigated this function in a series of notes starting from 1933 in the framework of the asymptotic theory of partitions; see [29–31]. The function is defined by the series representation, convergent in the whole z-complex plane(3.1)Wλ,μ(z):=∑n=0∞znn!Γ(λn+μ),λ>-1,μ∈ℂ.
Originally, Wright assumed that λ≥0, and, only in 1940 [32], he considered -1<λ<0. We note that in Chapter 18 of Vol. 3 of the handbook of the Bateman Project [33], devoted to Miscellaneous Functions, presumably for a misprint, the parameter λ of the Wright function is restricted to be nonnegative. When necessary, we propose to distinguish the Wright functions in two kinds according to λ≥0 (first kind) and -1<λ<0 (second kind).For more details on Wright functions the reader can consult, for example, [34–41] and references therein.Theintegral representation of the Wright function reads(3.2)Wλ,μ(z)=12πi∫Haeσ+zσ-λdσσμ,λ>-1,μ∈ℂ,
where Ha denotes the Hankel path. We remind that the Hankel path is a loop that starts from -∞ along the lower side of the negative real axis, encircles the circular area around the origin with radius ϵ→0 in the positive sense, and ends at -∞ along the upper side of the negative real axis. The equivalence of the series and integral representations is easily proved using Hankel formula for the Gamma function(3.3)1Γ(ζ)=∫Haeuu-ζdu,ζ∈ℂ
and performing a term-by-term integration. In fact,(3.4)Wλ,μ(z)=12πi∫Haeσ+zσ-λdσσμ=12πi∫Haeσ[∑n=0∞znn!σ-λn]dσσμ=∑n=0∞znn![12πi∫Haeσσ-λn-μdσ]=∑n=0∞znn!Γ[λn+μ].
It is possible to prove that the Wright function is entire of order 1/(1+λ); hence, it is of exponential type only if λ≥0 (which corresponds to Wright function of the first kind). The case λ=0 is trivial since W0,μ(z)=ez/Γ(μ), provided that μ≠0,-1,-2,….


### 3.2. The Auxiliary Functions of the Wright Type

Mainardi, in his first analysis of the time-fractional diffusion equation [42, 43], aware of the Bateman handbook [33], but not yet of the 1940 paper by Wright [32], introduced the two (Wright-type) entire auxiliary functions,(3.5)Fν(z):=W-ν,0(-z),0<ν<1,Mν(z):=W-ν,1-ν(-z),0<ν<1,
interrelated through(3.6)Fν(z)=νzMν(z).As a matter of fact, functionsFν(z) and Mν(z) are particular cases of the Wright function of the second kind Wλ,μ(z) by setting λ=-ν and μ=0 or μ=1, respectively.Hereafter, we provide the series and integral representations of the two auxiliary functions derived from the general formulas (3.1) and (3.2), respectively.Theseries representations for the auxiliary functions read(3.7)Fν(z):=∑n=1∞(-z)nn!Γ(-νn)=1π∑n=1∞(-z)n-1n!Γ(νn+1)sin(πνn),(3.8)Mν(z):=∑n=0∞(-z)nn!Γ[-νn+(1-ν)]=1π∑n=1∞(-z)n-1(n-1)!Γ(νn)sin(πνn).
The second series representations in (3.7)-(3.8) have been obtained by using the reflection formula for the Gamma functionΓ(ζ)Γ(1-ζ)=π/sinπζ.As an exercise, the reader can directly prove that the radius of convergence of the power series in (3.7)-(3.8) is infinite for 0<ν<1 without being aware of Wright's results, as it was shown independently by Mainardi [42]; see also [27].Furthermore, we haveFν(0)=0 and Mν(0)=1/Γ(1-ν). We note that relation (3.6) between the two auxiliary functions can be easily deduced from (3.7)-(3.8), by using the basic property of the Gamma function Γ(ζ+1)=ζΓ(ζ).Theintegral representations for the auxiliary functions read(3.9)Fν(z):=12πi∫Haeσ-zσνdσ,(3.10)Mν(z):=12πi∫Haeσ-zσνdσσ1-ν.
We note that relation (3.6) can be obtained also from (3.9)-(3.10) with an integration by parts. In fact,(3.11)Mν(z)=∫Haeσ-zσνdσσ1-ν=∫Haeσ(-1νzddσe-zσν)dσ=1νz∫Haeσ-zσνdσ=Fν(z)νz.The equivalence of the series and integral representations is easily proved by using the Hankel formula for the Gamma function and performing a term-by-term integration.


### 3.3. Special Cases

Explicit expressions ofFν(z) and Mν(z) in terms of known functions are expected for some particular values of ν. Mainardi and Tomirotti [43] have shown that for ν=1/q, where q≥2is a positive integer, the auxiliary functions can be expressed as a sum of simpler (q-1) entire functions. In the particular cases q=2 and q=3, we find(3.12)M1/2(z)=1π∑m=0∞(-1)m(12)mz2m(2m)!=1πexp(-z24),(3.13)M1/3(z)=1Γ(2/3)∑m=0∞(13)mz3m(3m)!-1Γ(1/3)∑m=0∞(23)mz3m+1(3m+1)!=32/3Ai(z31/3),
where Ai denotes the Airy function.Furthermore, it can be proved thatM1/q(z) satisfies the differential equation of order q-1(3.14)dq-1dzq-1M1/q(z)+(-1)qqzM1/q(z)=0,
subjected to the q-1 initial conditions at z=0, derived from (3.14), such that(3.15)M1/q(h)(0)=(-1)hΓ[1-(h+1)/q]=(-1)hπΓ[h+1q]sin[π(h+1)q],
with h=0,1,…,q-2. We note that, for q≥4, (3.14) is akin to the hyper-Airy differential equation of order q-1; see, for example, [44]. Consequently, the auxiliary function Mν(z) could be considered a sort of generalized hyper-Airy function. However, in view of further applications in stochastic processes, we prefer to consider it as a natural (fractional) generalization of the Gaussian function, similarly as the Mittag-Leffler function is known to be the natural (fractional) generalization of the exponential function.To stress the relevance of the auxiliary functionMν(z), it was also suggested the special name M-Wright function, a terminology that has been followed in literature to some extent.Some authors including Podlubny [27], Gorenflo et al. [34, 35], Hanyga [45], Balescu [46], Chechkin et al. [12], Germano et al. [47], and Kiryakova [48, 49] refer to the M-Wright function as the Mainardi function. It was Professor Stanković, during the presentation of the paper by Mainardi and Tomirotti [43] at the Conference of Transform Methods and Special Functions, Sofia 1994, who informed Mainardi, being aware only of the Bateman Handbook [33], that the extension for -1<λ<0 had been already made just by Wright himself in 1940 [32], following his previous papers published in the thirties. Mainardi, in the paper [50] devoted to the 80th birthday of Prof. Stanković, used the occasion to renew his personal gratitude to Prof. Stanković for this earlier information that led him to study the original papers by Wright and work also in collaboration on the functions of the Wright type for further applications; see, for example, [34, 35, 51].Moreover, the analysis of the limiting casesν=0 and ν=1 requires special attention. For ν=0 we easily recognize from the series representations (3.7)-(3.8)(3.16)F0(z)≡0,M0(z)=e-z.
The limiting case ν=1 is singular for both auxiliary functions as expected from the definition of the general Wright function when λ=-ν=-1. Later we will deal with this singular case for the M-Wright function when the variable is real and positive.


## 3.1. The General Wright Function

The Wright function, which we denote byWλ,μ,(z), is so named in honor of E. Maitland Wright, the eminent British mathematician, who introduced and investigated this function in a series of notes starting from 1933 in the framework of the asymptotic theory of partitions; see [29–31]. The function is defined by the series representation, convergent in the whole z-complex plane(3.1)Wλ,μ(z):=∑n=0∞znn!Γ(λn+μ),λ>-1,μ∈ℂ.
Originally, Wright assumed that λ≥0, and, only in 1940 [32], he considered -1<λ<0. We note that in Chapter 18 of Vol. 3 of the handbook of the Bateman Project [33], devoted to Miscellaneous Functions, presumably for a misprint, the parameter λ of the Wright function is restricted to be nonnegative. When necessary, we propose to distinguish the Wright functions in two kinds according to λ≥0 (first kind) and -1<λ<0 (second kind).For more details on Wright functions the reader can consult, for example, [34–41] and references therein.Theintegral representation of the Wright function reads(3.2)Wλ,μ(z)=12πi∫Haeσ+zσ-λdσσμ,λ>-1,μ∈ℂ,
where Ha denotes the Hankel path. We remind that the Hankel path is a loop that starts from -∞ along the lower side of the negative real axis, encircles the circular area around the origin with radius ϵ→0 in the positive sense, and ends at -∞ along the upper side of the negative real axis. The equivalence of the series and integral representations is easily proved using Hankel formula for the Gamma function(3.3)1Γ(ζ)=∫Haeuu-ζdu,ζ∈ℂ
and performing a term-by-term integration. In fact,(3.4)Wλ,μ(z)=12πi∫Haeσ+zσ-λdσσμ=12πi∫Haeσ[∑n=0∞znn!σ-λn]dσσμ=∑n=0∞znn![12πi∫Haeσσ-λn-μdσ]=∑n=0∞znn!Γ[λn+μ].
It is possible to prove that the Wright function is entire of order 1/(1+λ); hence, it is of exponential type only if λ≥0 (which corresponds to Wright function of the first kind). The case λ=0 is trivial since W0,μ(z)=ez/Γ(μ), provided that μ≠0,-1,-2,….


## 3.2. The Auxiliary Functions of the Wright Type

Mainardi, in his first analysis of the time-fractional diffusion equation [42, 43], aware of the Bateman handbook [33], but not yet of the 1940 paper by Wright [32], introduced the two (Wright-type) entire auxiliary functions,(3.5)Fν(z):=W-ν,0(-z),0<ν<1,Mν(z):=W-ν,1-ν(-z),0<ν<1,
interrelated through(3.6)Fν(z)=νzMν(z).As a matter of fact, functionsFν(z) and Mν(z) are particular cases of the Wright function of the second kind Wλ,μ(z) by setting λ=-ν and μ=0 or μ=1, respectively.Hereafter, we provide the series and integral representations of the two auxiliary functions derived from the general formulas (3.1) and (3.2), respectively.Theseries representations for the auxiliary functions read(3.7)Fν(z):=∑n=1∞(-z)nn!Γ(-νn)=1π∑n=1∞(-z)n-1n!Γ(νn+1)sin(πνn),(3.8)Mν(z):=∑n=0∞(-z)nn!Γ[-νn+(1-ν)]=1π∑n=1∞(-z)n-1(n-1)!Γ(νn)sin(πνn).
The second series representations in (3.7)-(3.8) have been obtained by using the reflection formula for the Gamma functionΓ(ζ)Γ(1-ζ)=π/sinπζ.As an exercise, the reader can directly prove that the radius of convergence of the power series in (3.7)-(3.8) is infinite for 0<ν<1 without being aware of Wright's results, as it was shown independently by Mainardi [42]; see also [27].Furthermore, we haveFν(0)=0 and Mν(0)=1/Γ(1-ν). We note that relation (3.6) between the two auxiliary functions can be easily deduced from (3.7)-(3.8), by using the basic property of the Gamma function Γ(ζ+1)=ζΓ(ζ).Theintegral representations for the auxiliary functions read(3.9)Fν(z):=12πi∫Haeσ-zσνdσ,(3.10)Mν(z):=12πi∫Haeσ-zσνdσσ1-ν.
We note that relation (3.6) can be obtained also from (3.9)-(3.10) with an integration by parts. In fact,(3.11)Mν(z)=∫Haeσ-zσνdσσ1-ν=∫Haeσ(-1νzddσe-zσν)dσ=1νz∫Haeσ-zσνdσ=Fν(z)νz.The equivalence of the series and integral representations is easily proved by using the Hankel formula for the Gamma function and performing a term-by-term integration.


## 3.3. Special Cases

Explicit expressions ofFν(z) and Mν(z) in terms of known functions are expected for some particular values of ν. Mainardi and Tomirotti [43] have shown that for ν=1/q, where q≥2is a positive integer, the auxiliary functions can be expressed as a sum of simpler (q-1) entire functions. In the particular cases q=2 and q=3, we find(3.12)M1/2(z)=1π∑m=0∞(-1)m(12)mz2m(2m)!=1πexp(-z24),(3.13)M1/3(z)=1Γ(2/3)∑m=0∞(13)mz3m(3m)!-1Γ(1/3)∑m=0∞(23)mz3m+1(3m+1)!=32/3Ai(z31/3),
where Ai denotes the Airy function.Furthermore, it can be proved thatM1/q(z) satisfies the differential equation of order q-1(3.14)dq-1dzq-1M1/q(z)+(-1)qqzM1/q(z)=0,
subjected to the q-1 initial conditions at z=0, derived from (3.14), such that(3.15)M1/q(h)(0)=(-1)hΓ[1-(h+1)/q]=(-1)hπΓ[h+1q]sin[π(h+1)q],
with h=0,1,…,q-2. We note that, for q≥4, (3.14) is akin to the hyper-Airy differential equation of order q-1; see, for example, [44]. Consequently, the auxiliary function Mν(z) could be considered a sort of generalized hyper-Airy function. However, in view of further applications in stochastic processes, we prefer to consider it as a natural (fractional) generalization of the Gaussian function, similarly as the Mittag-Leffler function is known to be the natural (fractional) generalization of the exponential function.To stress the relevance of the auxiliary functionMν(z), it was also suggested the special name M-Wright function, a terminology that has been followed in literature to some extent.Some authors including Podlubny [27], Gorenflo et al. [34, 35], Hanyga [45], Balescu [46], Chechkin et al. [12], Germano et al. [47], and Kiryakova [48, 49] refer to the M-Wright function as the Mainardi function. It was Professor Stanković, during the presentation of the paper by Mainardi and Tomirotti [43] at the Conference of Transform Methods and Special Functions, Sofia 1994, who informed Mainardi, being aware only of the Bateman Handbook [33], that the extension for -1<λ<0 had been already made just by Wright himself in 1940 [32], following his previous papers published in the thirties. Mainardi, in the paper [50] devoted to the 80th birthday of Prof. Stanković, used the occasion to renew his personal gratitude to Prof. Stanković for this earlier information that led him to study the original papers by Wright and work also in collaboration on the functions of the Wright type for further applications; see, for example, [34, 35, 51].Moreover, the analysis of the limiting casesν=0 and ν=1 requires special attention. For ν=0 we easily recognize from the series representations (3.7)-(3.8)(3.16)F0(z)≡0,M0(z)=e-z.
The limiting case ν=1 is singular for both auxiliary functions as expected from the definition of the general Wright function when λ=-ν=-1. Later we will deal with this singular case for the M-Wright function when the variable is real and positive.


## 4. Properties and Plots of the Auxiliary Wright Functions  in Real Domain

Let us state some relevant properties of the auxiliary Wright functions, with special attention to theMν function in view of its role in time-fractional diffusion processes.


### 4.1. Exponential Laplace Transforms

We start with the Laplace transform pairs involving exponentials in the Laplace domain. These were derived by Mainardi in his earlier analysis of the time-fractional diffusion equation; see, for example, [42, 52],(4.1)1rFν(1rν)=νrν+1Mν(1rν)↔ℒe-sν,0<ν<1,(4.2)1νFν(1rν)=1rνMν(1rν)↔ℒe-sνs1-ν,0<ν<1.
We note that the inversion of the Laplace transform of the exponential exp(-sν) is relevant since it yields for any ν∈(0,1) the unilateral extremal stable densities in probability theory, denoted by Lν-ν(r) in [53]. As a consequence, the nonnegativity of both auxiliary Wright functions when their argument is positive is proved by the Bernstein theorem. We refer to Feller's treatise [54, 55] for Laplace transforms, stable densities and Bernstein theorem. The Laplace transform pair in (4.1) has a long history starting from a formal result by Humbert [56] in 1945, of which Pollard [57] provided a rigorous proof one year later. Then, in 1959 Mikusiński [58] derived a similar result on the basis of his theory of operational calculus. In 1975, albeit unaware of the previous results, Buchen and Mainardi [59] derived the result in a formal way. We note that all the above authors were not informed about the Wright functions. To our actual knowledge the former author who derived the Laplace transforms pairs (4.1)-(4.2) in terms of Wright functions of the second kind was Stanković in 1970; see [39].Hereafter we would like to provide two independent proofs of (4.1) carrying out the inversion of exp(-sν), either by the complex Bromwich integral formula following [42], or by the formal series method following [59]. Similarly we can act for the Laplace transform pair (4.2). For the complex integral approach we deform the Bromwich path Br into the Hankel path Ha, that is equivalent to the original path, and we set σ=sr. Recalling the integral representation (3.9) for the Fν function and (3.6), we get(4.3)ℒ-1[exp(-sν);s→r]=12πi∫Bresr-sνds=12πir∫Haeσ-(σ/r)νdσ=1rFν(1rν)=νrν+1Mν(1rν).
Expanding in power series the Laplace transform and inverting term by term, we formally get(4.4)ℒ-1[exp(-sν)]=∑n=0∞(-1)nn!ℒ-1[sνn]=∑n=1∞(-1)nn!r-νn-1Γ(-νn)=1rFν(1rν)=νrν+1Mν(1rν),
where now we have used the series representation (3.7) for the function Fν along with the relationship formula (3.6).


### 4.2. Asymptotic Representation for Large Argument

Let us point out the asymptotic behavior of the functionMν(r) when r→∞. Choosing a variable r/ν rather than r, the computation of the desired asymptotic representation by the saddle-point approximation is straightforward. Mainardi and Tomirotti [43] have obtained(4.5)Mν(rν)~a(ν)r(ν-1/2)/(1-ν)exp[-b(ν)r1/(1-ν)],a(ν)=12π(1-ν)>0,b(ν)=1-νν>0.
The above evaluation is consistent with the first term in the asymptotic series expansion provided by Wright with a cumbersome and formal procedure for his general function Wλ,μ when -1<λ<0; see [32]. In 1999 Wong and Zhao have derived asymptotic expansions of the Wright functions of the first and second kind in the whole complex plane following a new method for smoothing Stokes' discontinuities; see [40, 41], respectively.We note that, forν=1/2 as (4.5) provides the exact result consistent with (3.12),(4.6)M1/2(2r)=1πe-r2⇔M1/2(r)=1πe-r2/4.
We also note that in the limit ν→1- the function Mν(r) tends to the Dirac generalized function δ(r-1), as can be recognized also from the Laplace transform pair (4.1).


### 4.3. Absolute Moments

From the above considerations we recognize that, for theM-Wright functions, the following rule for absolute moments in ℝ+ holds(4.7)∫0∞rδMν(r)dr=Γ(δ+1)Γ(νδ+1),δ>-1,0≤ν<1.
In order to derive this fundamental result, we proceed as follows on the basis of the integral representation (3.10):(4.8)∫0∞rδMν(r)dr=∫0∞rδ[12πi∫Haeσ-rσνdσσ1-ν]dr=12πi∫Haeσ[∫0∞e-rσνrδdr]dσσ1-ν=Γ(δ+1)2πi∫Haeσσνδ+1dσ=Γ(δ+1)Γ(νδ+1).
Above we have legitimized the exchange between integrals and used the identity(4.9)∫0∞e-rσνrδdr=Γ(δ+1)(σν)δ+1,
along with the Hankel formula of the Gamma function. Analogously, we can compute all the moments of Fν(r) in ℝ+.


### 4.4. The Laplace Transform of theM-Wright Function

Let the Mittag-Leffler function be defined in the complex plane for anyν≥0 by the following series and integral representation; see, for example, [33, 60]:(4.10)Eν(z)=∑n=0∞znΓ(νn+1)=12πi∫Haζν-1eζζν-zdζ,ν>0,z∈ℂ.
Such function is entire of order 1/α for α>0 and reduces to the function exp(z) for ν>0 and to 1/(1-z) for ν=0. We recall that the Mittag-Leffler function for ν>0 plays fundamental roles in applications of fractional calculus like fractional relaxation and fractional oscillation; see, for example, [25, 26, 61, 62], so that it could be referred as the Queen function of fractional calculus. Recently, numerical routines for functions of Mittag-Leffler type have been provided, e.g., by Freed et al. [63], Gorenflo et al. [64] (with MATHEMATICA), Podlubny [65] (with MATLAB), Seybold and Hilfer [66].We now point out that theM-Wright function is related to the Mittag-Leffler function through the following Laplace transform pair:(4.11)Mν(r)↔ℒEν(-s),0<ν<1.
For the reader's convenience, we provide a simple proof of (4.11) by using two different approaches. We assume that the exchanges between integrals and series are legitimate in view of the analyticity properties of the involved functions. In the first approach we use the integral representations of the two functions obtaining(4.12)∫0∞e-srMν(r)dr=12πi∫0∞e-sr[∫Haeσ-rσνdσσ1-ν]dr=12πi∫Haeσσν-1[∫0∞e-r(s+σν)dr]dσ=12πi∫Haeσσν-1σν+sdσ=Eν(-s).
In the second approach we develop in series the exponential kernel of the Laplace transform and we use the expression (4.7) for the absolute moments of the M-Wright function arriving to the following series representation of the Mittag-Leffler function:(4.13)∫0∞e-srMν(r)dr=∑n=0∞(-s)nn!∫0∞rnMν(r)dr=∑n=0∞(-s)nn!Γ(n+1)Γ(νn+1)=∑n=0∞(-s)nΓ(νn+1)=Eν(-s).
We note that the transformation term by term of the series expansion of the M-Wright function is not legitimate because the function is not of exponential order; see [67]. However, this procedure yields the formal asymptotic expansion of the Mittag-Leffler function Eν(-s) as s→∞ in a sector around the positive real axis; see, for example, [33, 60], that is,(4.14)∑n=0∞∫0∞e-sr(-r)ndrn!Γ(-νn+(1-ν))=∑n=0∞(-1)nΓ(-νn+1-ν)1sn+1=∑m=1∞(-1)m-1Γ(-νm+1)1sm~Eν(-s),s→∞.


### 4.5. The Fourier Transform of the SymmetricM-Wright Function

TheM-Wright function, extended on the negative real axis as an even function, is related to the Mittag-Leffler function through the following Fourier transform pair:(4.15)Mν(|x|)↔ℱ2E2ν(-κ2),0<ν<1.
Below, we prove the equivalent formula(4.16)∫0∞cos(κr)Mν(r)dr=E2ν(-κ2).
For the proof it is sufficient to develop in series the cosine function and use formula (4.7) for the absolute moments of the M-Wright function:(4.17)∫0∞cos(κr)Mν(r)dr=∑n=0∞(-1)nκ2n(2n)!∫0∞r2nMν(r)dr=∑n=0∞(-1)nκ2nΓ(2νn+1)=E2ν(-κ2).


### 4.6. The Mellin Transform of theM-Wright Function

It is straightforward to derive the Mellin transform of theM-Wright function using result (4.7) for the absolute moments of the M-Wright function. In fact, setting δ=s-1 in (4.7), by analytic continuation it follows:(4.18)Mν(r)↔ℳΓ(s)Γ(ν(s-1)+1),0<ν<1.


### 4.7. Plots of the SymmetricM-Wright Function

It is instructive to show the plots of the (symmetric)M-Wright function on the real axis for some rational values of the parameter ν. In order to have more insight of the effect of the parameter itself on the behavior close to and far from the origin, we adopt both linear and logarithmic scale for the ordinates.In Figures1 and 2 we compare the plots of the Mν(x)-Wright functions in-5≤x≤5 for some rational values of ν in the ranges ν∈[0,1/2] and ν∈[1/2,1], respectively. In Figure 1 we see the transition from exp(-|x|) for ν=0 to 1/πexp(-x2) for ν=1/2, whereas in Figure 2 we see the transition from 1/πexp(-x2) to the delta functions δ(x±1) for ν=1. Because of the two symmetrical humps for 1/2<ν≤1, the Mν function appears bimodal with the characteristic shape of the capital letter M.Plots of the symmetricMν-Wright function with ν=0,1/8,1/4,3/8,1/2 for-5≤x≤5: (a) linear scale, (b) logarithmic scale.
(a)(b)Plots of the symmetricM-Wright function with ν=1/2,5/8,3/4,1 for-5≤x≤5: (a) linear scale, (b) logarithmic scale.
(a)(b)In plottingMν(x) at fixed ν for sufficiently large x, the asymptotic representation (4.5)-(4.6) is useful since, as x increases, the numerical convergence of the series in (3.8) decreases up to being completely inefficient: henceforth, the matching between the series and the asymptotic representation is relevant and followed by Mainardi and associates; see, for example, [38, 51, 53, 68]. However, as ν→1-, the plotting remains a very difficult task because of the high peak arising around x=±1. For this we refer the reader to the 1997 paper by Mainardi and Tomirotti [69], where a variant of the saddle point method has been successfully used to properly depict the transition to the delta functions δ(x±1) as ν approaches 1. For the numerical point of view we like to highlight the recent paper by Luchko [70], where algorithms are provided for computation of the Wright function on the real axis with prescribed accuracy.


### 4.8. The𝕄-Wright Function in Two Variables

In view of the time-fractional diffusion processes that will be considered in the next sections, it is worthwhile to introduce the function in the two variables(4.19)𝕄ν(x,t):=t-νMν(xt-ν),0<ν<1,x,t∈ℝ+,
which defines a spatial probability density in x evolving in time t with self-similarity exponent H=ν. Of course for x∈ℝ we have to consider the symmetric version obtained from (4.19) multiplying by 1/2 and replacing x by |x|.Hereafter we provide a list of the main properties of this function, which can be derived from Laplace and Fourier transforms of the correspondingM-Wright function in one variable.From (4.2) we derive the Laplace transform of 𝕄ν(x,t) with respect to t∈ℝ+ as(4.20)ℒ{𝕄ν(x,t);t→s}=sν-1e-xsν.
From (4.10) we derive the Laplace transform of 𝕄ν(x,t) with respect to x∈ℝ+ as(4.21)ℒ{𝕄ν(x,t);x→s}=Eν(-stν).
From (4.15) we derive the Fourier transform of 𝕄ν(|x|,t) with respect to x∈ℝ as(4.22)ℱ{𝕄ν(|x|,t);x→κ}=2E2ν(-κ2tν).
Moreover, using the Mellin transform, Mainardi et al. [71] derived the following integral formula:(4.23)𝕄ν(x,t)=∫0∞𝕄λ(x,τ)𝕄μ(τ,t)dτ,ν=λμ.
Special cases of the 𝕄-Wright function are simply derived for ν=1/2 and ν=1/3 from the corresponding ones in the complex domain; see (3.12)-(3.13). We devote particular attention to the case ν=1/2 for which we get from (4.6) the Gaussian density in ℝ:(4.24)𝕄1/2(|x|,t)=12πt1/2e-x2/(4t).
For the limiting case ν=1 we obtain(4.25)𝕄1(|x|,t)=12[δ(x-t)+δ(x+t)].


## 4.1. Exponential Laplace Transforms

We start with the Laplace transform pairs involving exponentials in the Laplace domain. These were derived by Mainardi in his earlier analysis of the time-fractional diffusion equation; see, for example, [42, 52],(4.1)1rFν(1rν)=νrν+1Mν(1rν)↔ℒe-sν,0<ν<1,(4.2)1νFν(1rν)=1rνMν(1rν)↔ℒe-sνs1-ν,0<ν<1.
We note that the inversion of the Laplace transform of the exponential exp(-sν) is relevant since it yields for any ν∈(0,1) the unilateral extremal stable densities in probability theory, denoted by Lν-ν(r) in [53]. As a consequence, the nonnegativity of both auxiliary Wright functions when their argument is positive is proved by the Bernstein theorem. We refer to Feller's treatise [54, 55] for Laplace transforms, stable densities and Bernstein theorem. The Laplace transform pair in (4.1) has a long history starting from a formal result by Humbert [56] in 1945, of which Pollard [57] provided a rigorous proof one year later. Then, in 1959 Mikusiński [58] derived a similar result on the basis of his theory of operational calculus. In 1975, albeit unaware of the previous results, Buchen and Mainardi [59] derived the result in a formal way. We note that all the above authors were not informed about the Wright functions. To our actual knowledge the former author who derived the Laplace transforms pairs (4.1)-(4.2) in terms of Wright functions of the second kind was Stanković in 1970; see [39].Hereafter we would like to provide two independent proofs of (4.1) carrying out the inversion of exp(-sν), either by the complex Bromwich integral formula following [42], or by the formal series method following [59]. Similarly we can act for the Laplace transform pair (4.2). For the complex integral approach we deform the Bromwich path Br into the Hankel path Ha, that is equivalent to the original path, and we set σ=sr. Recalling the integral representation (3.9) for the Fν function and (3.6), we get(4.3)ℒ-1[exp(-sν);s→r]=12πi∫Bresr-sνds=12πir∫Haeσ-(σ/r)νdσ=1rFν(1rν)=νrν+1Mν(1rν).
Expanding in power series the Laplace transform and inverting term by term, we formally get(4.4)ℒ-1[exp(-sν)]=∑n=0∞(-1)nn!ℒ-1[sνn]=∑n=1∞(-1)nn!r-νn-1Γ(-νn)=1rFν(1rν)=νrν+1Mν(1rν),
where now we have used the series representation (3.7) for the function Fν along with the relationship formula (3.6).


## 4.2. Asymptotic Representation for Large Argument

Let us point out the asymptotic behavior of the functionMν(r) when r→∞. Choosing a variable r/ν rather than r, the computation of the desired asymptotic representation by the saddle-point approximation is straightforward. Mainardi and Tomirotti [43] have obtained(4.5)Mν(rν)~a(ν)r(ν-1/2)/(1-ν)exp[-b(ν)r1/(1-ν)],a(ν)=12π(1-ν)>0,b(ν)=1-νν>0.
The above evaluation is consistent with the first term in the asymptotic series expansion provided by Wright with a cumbersome and formal procedure for his general function Wλ,μ when -1<λ<0; see [32]. In 1999 Wong and Zhao have derived asymptotic expansions of the Wright functions of the first and second kind in the whole complex plane following a new method for smoothing Stokes' discontinuities; see [40, 41], respectively.We note that, forν=1/2 as (4.5) provides the exact result consistent with (3.12),(4.6)M1/2(2r)=1πe-r2⇔M1/2(r)=1πe-r2/4.
We also note that in the limit ν→1- the function Mν(r) tends to the Dirac generalized function δ(r-1), as can be recognized also from the Laplace transform pair (4.1).


## 4.3. Absolute Moments

From the above considerations we recognize that, for theM-Wright functions, the following rule for absolute moments in ℝ+ holds(4.7)∫0∞rδMν(r)dr=Γ(δ+1)Γ(νδ+1),δ>-1,0≤ν<1.
In order to derive this fundamental result, we proceed as follows on the basis of the integral representation (3.10):(4.8)∫0∞rδMν(r)dr=∫0∞rδ[12πi∫Haeσ-rσνdσσ1-ν]dr=12πi∫Haeσ[∫0∞e-rσνrδdr]dσσ1-ν=Γ(δ+1)2πi∫Haeσσνδ+1dσ=Γ(δ+1)Γ(νδ+1).
Above we have legitimized the exchange between integrals and used the identity(4.9)∫0∞e-rσνrδdr=Γ(δ+1)(σν)δ+1,
along with the Hankel formula of the Gamma function. Analogously, we can compute all the moments of Fν(r) in ℝ+.


## 4.4. The Laplace Transform of theM-Wright Function

Let the Mittag-Leffler function be defined in the complex plane for anyν≥0 by the following series and integral representation; see, for example, [33, 60]:(4.10)Eν(z)=∑n=0∞znΓ(νn+1)=12πi∫Haζν-1eζζν-zdζ,ν>0,z∈ℂ.
Such function is entire of order 1/α for α>0 and reduces to the function exp(z) for ν>0 and to 1/(1-z) for ν=0. We recall that the Mittag-Leffler function for ν>0 plays fundamental roles in applications of fractional calculus like fractional relaxation and fractional oscillation; see, for example, [25, 26, 61, 62], so that it could be referred as the Queen function of fractional calculus. Recently, numerical routines for functions of Mittag-Leffler type have been provided, e.g., by Freed et al. [63], Gorenflo et al. [64] (with MATHEMATICA), Podlubny [65] (with MATLAB), Seybold and Hilfer [66].We now point out that theM-Wright function is related to the Mittag-Leffler function through the following Laplace transform pair:(4.11)Mν(r)↔ℒEν(-s),0<ν<1.
For the reader's convenience, we provide a simple proof of (4.11) by using two different approaches. We assume that the exchanges between integrals and series are legitimate in view of the analyticity properties of the involved functions. In the first approach we use the integral representations of the two functions obtaining(4.12)∫0∞e-srMν(r)dr=12πi∫0∞e-sr[∫Haeσ-rσνdσσ1-ν]dr=12πi∫Haeσσν-1[∫0∞e-r(s+σν)dr]dσ=12πi∫Haeσσν-1σν+sdσ=Eν(-s).
In the second approach we develop in series the exponential kernel of the Laplace transform and we use the expression (4.7) for the absolute moments of the M-Wright function arriving to the following series representation of the Mittag-Leffler function:(4.13)∫0∞e-srMν(r)dr=∑n=0∞(-s)nn!∫0∞rnMν(r)dr=∑n=0∞(-s)nn!Γ(n+1)Γ(νn+1)=∑n=0∞(-s)nΓ(νn+1)=Eν(-s).
We note that the transformation term by term of the series expansion of the M-Wright function is not legitimate because the function is not of exponential order; see [67]. However, this procedure yields the formal asymptotic expansion of the Mittag-Leffler function Eν(-s) as s→∞ in a sector around the positive real axis; see, for example, [33, 60], that is,(4.14)∑n=0∞∫0∞e-sr(-r)ndrn!Γ(-νn+(1-ν))=∑n=0∞(-1)nΓ(-νn+1-ν)1sn+1=∑m=1∞(-1)m-1Γ(-νm+1)1sm~Eν(-s),s→∞.


## 4.5. The Fourier Transform of the SymmetricM-Wright Function

TheM-Wright function, extended on the negative real axis as an even function, is related to the Mittag-Leffler function through the following Fourier transform pair:(4.15)Mν(|x|)↔ℱ2E2ν(-κ2),0<ν<1.
Below, we prove the equivalent formula(4.16)∫0∞cos(κr)Mν(r)dr=E2ν(-κ2).
For the proof it is sufficient to develop in series the cosine function and use formula (4.7) for the absolute moments of the M-Wright function:(4.17)∫0∞cos(κr)Mν(r)dr=∑n=0∞(-1)nκ2n(2n)!∫0∞r2nMν(r)dr=∑n=0∞(-1)nκ2nΓ(2νn+1)=E2ν(-κ2).


## 4.6. The Mellin Transform of theM-Wright Function

It is straightforward to derive the Mellin transform of theM-Wright function using result (4.7) for the absolute moments of the M-Wright function. In fact, setting δ=s-1 in (4.7), by analytic continuation it follows:(4.18)Mν(r)↔ℳΓ(s)Γ(ν(s-1)+1),0<ν<1.


## 4.7. Plots of the SymmetricM-Wright Function

It is instructive to show the plots of the (symmetric)M-Wright function on the real axis for some rational values of the parameter ν. In order to have more insight of the effect of the parameter itself on the behavior close to and far from the origin, we adopt both linear and logarithmic scale for the ordinates.In Figures1 and 2 we compare the plots of the Mν(x)-Wright functions in-5≤x≤5 for some rational values of ν in the ranges ν∈[0,1/2] and ν∈[1/2,1], respectively. In Figure 1 we see the transition from exp(-|x|) for ν=0 to 1/πexp(-x2) for ν=1/2, whereas in Figure 2 we see the transition from 1/πexp(-x2) to the delta functions δ(x±1) for ν=1. Because of the two symmetrical humps for 1/2<ν≤1, the Mν function appears bimodal with the characteristic shape of the capital letter M.Plots of the symmetricMν-Wright function with ν=0,1/8,1/4,3/8,1/2 for-5≤x≤5: (a) linear scale, (b) logarithmic scale.
(a)(b)Plots of the symmetricM-Wright function with ν=1/2,5/8,3/4,1 for-5≤x≤5: (a) linear scale, (b) logarithmic scale.
(a)(b)In plottingMν(x) at fixed ν for sufficiently large x, the asymptotic representation (4.5)-(4.6) is useful since, as x increases, the numerical convergence of the series in (3.8) decreases up to being completely inefficient: henceforth, the matching between the series and the asymptotic representation is relevant and followed by Mainardi and associates; see, for example, [38, 51, 53, 68]. However, as ν→1-, the plotting remains a very difficult task because of the high peak arising around x=±1. For this we refer the reader to the 1997 paper by Mainardi and Tomirotti [69], where a variant of the saddle point method has been successfully used to properly depict the transition to the delta functions δ(x±1) as ν approaches 1. For the numerical point of view we like to highlight the recent paper by Luchko [70], where algorithms are provided for computation of the Wright function on the real axis with prescribed accuracy.


## 4.8. The𝕄-Wright Function in Two Variables

In view of the time-fractional diffusion processes that will be considered in the next sections, it is worthwhile to introduce the function in the two variables(4.19)𝕄ν(x,t):=t-νMν(xt-ν),0<ν<1,x,t∈ℝ+,
which defines a spatial probability density in x evolving in time t with self-similarity exponent H=ν. Of course for x∈ℝ we have to consider the symmetric version obtained from (4.19) multiplying by 1/2 and replacing x by |x|.Hereafter we provide a list of the main properties of this function, which can be derived from Laplace and Fourier transforms of the correspondingM-Wright function in one variable.From (4.2) we derive the Laplace transform of 𝕄ν(x,t) with respect to t∈ℝ+ as(4.20)ℒ{𝕄ν(x,t);t→s}=sν-1e-xsν.
From (4.10) we derive the Laplace transform of 𝕄ν(x,t) with respect to x∈ℝ+ as(4.21)ℒ{𝕄ν(x,t);x→s}=Eν(-stν).
From (4.15) we derive the Fourier transform of 𝕄ν(|x|,t) with respect to x∈ℝ as(4.22)ℱ{𝕄ν(|x|,t);x→κ}=2E2ν(-κ2tν).
Moreover, using the Mellin transform, Mainardi et al. [71] derived the following integral formula:(4.23)𝕄ν(x,t)=∫0∞𝕄λ(x,τ)𝕄μ(τ,t)dτ,ν=λμ.
Special cases of the 𝕄-Wright function are simply derived for ν=1/2 and ν=1/3 from the corresponding ones in the complex domain; see (3.12)-(3.13). We devote particular attention to the case ν=1/2 for which we get from (4.6) the Gaussian density in ℝ:(4.24)𝕄1/2(|x|,t)=12πt1/2e-x2/(4t).
For the limiting case ν=1 we obtain(4.25)𝕄1(|x|,t)=12[δ(x-t)+δ(x+t)].


## 5. Fractional Diffusion Equations

Let us now consider a variety of diffusion-like equations starting from the standard diffusion equation whose fundamental solutions are expressed in terms of theM-Wright function depending on space and time variables. The two variables, however, turn out to be related through a self-similarity condition.


### 5.1. The Standard Diffusion Equation

The standard diffusion equation for the fieldu(x,t) with initial condition u(x,0)=u0(x) is(5.1)∂u∂t=K1∂2u∂x2,-∞<x<∞,t≥0,
where K1 is a suitable diffusion coefficient of dimensions [K1]=[L]2[T]-1=cm2/sec. This initial-boundary value problem can be easily shown to be equivalent to the Volterra integral equation(5.2)u(x,t)=u0(x)+K1∫0t∂2u(x,τ)∂x2dτ.
It is well known that the fundamental solution (usually referred as the Green function), which is the solution corresponding to u0(x)=δ(x), is the Gaussian probability density evolving in time with variance (mean square displacement) proportional to time. In our notation we have(5.3)𝒢1(x,t)=12πK1t1/2e-x2/(4K1t),(5.4)σ12(t):=∫-∞+∞x2𝒢1(x,t)dx=2K1t.
This variance law characterizes the process of normal diffusion as it emerges from Einstein's approach to Brownian motion (Bm); see, for example, [72].In view of future developments, we rewrite the Green function in terms of theM-Wright function by recalling (3.12), that is,(5.5)𝒢1(x,t)=121K1t1/2M1/2(|x|K1t1/2).
From the self-similarity of the Green function in (5.3) or (5.5), we are led to write(5.6)𝒢1(x,t)=1K1tH𝒢1(|x|K1tH,1),
where H=1/2 is the similarity (or Hurst) exponent and ξ=|x|/(K1t1/2) acts as the similarity variable. We refer to the one-variable function 𝒢1(ξ) as the reduced Green function.


### 5.2. The Stretched-Time Standard Diffusion Equation

Let us now stretch the time variable in (5.1) by replacing t with tα where 0<α<2. We have(5.7)∂u∂(tα)=Kα∂2u∂x2,-∞<x<+∞,t≥0,
where Kα is a sort of stretched diffusion coefficient of dimensions [Kα]=[L]2[T]-α=cm2/secα. It is easy to recognize that such equation is akin to the standard diffusion equation but with a diffusion coefficient depending on time; K1(t)=αtα-1Kα. In fact, using the rule(5.8)∂∂tα=1αtα-1∂∂t,
we have(5.9)∂u∂t=αtα-1Kα∂2u∂x2,-∞<x<+∞,t≥0.
The integral form corresponding to (5.7)—(5.9) reads(5.10)u(x,t)=u0(x)+αKα∫0t∂2u(x,τ)∂x2τα-1dτ.
The corresponding fundamental solution is the stretched-time Gaussian(5.11)𝒢α(x,t)=12πKαtα/2e-x2/(4Kαtα)=121Kαtα/2M1/2(|x|Kαtα/2).
The corresponding variance(5.12)σα2(t):=∫-∞+∞x2𝒢α(x,t)dx=2Kαtα
is characteristic of a general process of anomalous diffusion, precisely of slow diffusion for 0<α<1, and of fast diffusion for 1<α<2.


### 5.3. The Time-Fractional Diffusion Equation

In literature there exist two forms of the time-fractional diffusion equation of a single order, one with Riemann-Liouvile derivative and one with Caputo derivative. These forms are equivalent if we refer to the standard initial conditionu(x,0)=u0(x), as shown in [73].Taking a real numberβ∈(0,1), the time-fractional diffusion equation of order β in the Riemann-Liouville sense reads(5.13)∂u∂t=KβDt1-β∂2u∂x2,
whereas in the Caputo sense reads(5.14)*Dtβu=Kβ∂2u∂x2,
where Kβ is a sort of fractional diffusion coefficient of dimensions [Kβ]=[L]2[T]-β=cm2/secβ. Like for diffusion equations of integer order (5.1) and (5.7)-(5.9), we consider the equivalent integral equation corresponding to our fractional diffusion equations (5.13)-(5.14) as(5.15)u(x,t)=u0(x)+Kβ1Γ(β)∫0t(t-τ)β-1∂2u(x,τ)∂x2dτ.
The Green function 𝒢β(x,t) for the equivalent (5.13)–(5.15) can be expressed, also in this case, in terms of the M-Wright function, as shown in Appendix A by adopting two different approaches, as follows:(5.16)𝒢β(x,t)=121Kβtβ/2Mβ/2(|x|Kβtβ/2).
The corresponding variance can be promptly obtained from the general formula (5.5) for the absolute moment of the M-Wright function. In fact, using (5.5) and (5.16) and after an obvious change of variable, we obtain(5.17)σβ2(t):=∫-∞+∞x2𝒢β(x,t)dx=2Γ(β+1)Kβtβ.
As a consequence, for 0<β<1 the variance is consistent with a process of slow diffusion with similarity exponent H=β/2. For further reading on time-fractional diffusion equations and their solutions the reader is referred, among others, to [38, 51, 53] and [74, 75].


### 5.4. The Stretched Time-Fractional Diffusion Equation

In the fractional diffusion equation (5.13), let us stretch the time variable by replacing t with tα/β where 0<α<2 and 0<β≤1. We have(5.18)∂u∂tα/β=KαβDtα/β1-β∂2u∂x2,
namely,(5.19)∂u∂t=αβtα/β-1KαβDtα/β1-β∂2u∂x2,
where Kαβ is a sort of stretched diffusion coefficient of dimensions [Kαβ]=[L]2[T]-α=cm2/secα that reduces to Kα if β=1 and to Kβ if α=β. Integration of (5.19) gives the corresponding integral equation [21](5.20)u(x,t)=u0(x)+Kαβ1Γ(β)αβ∫0tτα/β-1(tα/β-τα/β)β-1∂2u(x,τ)∂x2dτ,
whose Green function 𝒢αβ(x,t) is(5.21)𝒢αβ(x,t)=121Kαβtα/2Mβ/2(|x|Kαβtα/2),
with variance(5.22)σα,β2(t):=∫-∞+∞x2𝒢α,β(x,t)dx=2Γ(β+1)Kαβtα.
As a consequence, the resulting process turns out to be self-similar with Hurst exponent H=α/2 and a variance law consistent with both slow diffusion if 0<α<1 and fast diffusion if 1<α<2. We note that the parameter β does explicitly enter in the variance law (5.22) only as in the determination of the multiplicative constant.It is straightforward to note that the evolution equations of this process reduce to those for time-fractional diffusion ifα=β<1, for stretched diffusion if α≠1 and β=1, and finally for standard diffusion if α=β=1.


## 5.1. The Standard Diffusion Equation

The standard diffusion equation for the fieldu(x,t) with initial condition u(x,0)=u0(x) is(5.1)∂u∂t=K1∂2u∂x2,-∞<x<∞,t≥0,
where K1 is a suitable diffusion coefficient of dimensions [K1]=[L]2[T]-1=cm2/sec. This initial-boundary value problem can be easily shown to be equivalent to the Volterra integral equation(5.2)u(x,t)=u0(x)+K1∫0t∂2u(x,τ)∂x2dτ.
It is well known that the fundamental solution (usually referred as the Green function), which is the solution corresponding to u0(x)=δ(x), is the Gaussian probability density evolving in time with variance (mean square displacement) proportional to time. In our notation we have(5.3)𝒢1(x,t)=12πK1t1/2e-x2/(4K1t),(5.4)σ12(t):=∫-∞+∞x2𝒢1(x,t)dx=2K1t.
This variance law characterizes the process of normal diffusion as it emerges from Einstein's approach to Brownian motion (Bm); see, for example, [72].In view of future developments, we rewrite the Green function in terms of theM-Wright function by recalling (3.12), that is,(5.5)𝒢1(x,t)=121K1t1/2M1/2(|x|K1t1/2).
From the self-similarity of the Green function in (5.3) or (5.5), we are led to write(5.6)𝒢1(x,t)=1K1tH𝒢1(|x|K1tH,1),
where H=1/2 is the similarity (or Hurst) exponent and ξ=|x|/(K1t1/2) acts as the similarity variable. We refer to the one-variable function 𝒢1(ξ) as the reduced Green function.


## 5.2. The Stretched-Time Standard Diffusion Equation

Let us now stretch the time variable in (5.1) by replacing t with tα where 0<α<2. We have(5.7)∂u∂(tα)=Kα∂2u∂x2,-∞<x<+∞,t≥0,
where Kα is a sort of stretched diffusion coefficient of dimensions [Kα]=[L]2[T]-α=cm2/secα. It is easy to recognize that such equation is akin to the standard diffusion equation but with a diffusion coefficient depending on time; K1(t)=αtα-1Kα. In fact, using the rule(5.8)∂∂tα=1αtα-1∂∂t,
we have(5.9)∂u∂t=αtα-1Kα∂2u∂x2,-∞<x<+∞,t≥0.
The integral form corresponding to (5.7)—(5.9) reads(5.10)u(x,t)=u0(x)+αKα∫0t∂2u(x,τ)∂x2τα-1dτ.
The corresponding fundamental solution is the stretched-time Gaussian(5.11)𝒢α(x,t)=12πKαtα/2e-x2/(4Kαtα)=121Kαtα/2M1/2(|x|Kαtα/2).
The corresponding variance(5.12)σα2(t):=∫-∞+∞x2𝒢α(x,t)dx=2Kαtα
is characteristic of a general process of anomalous diffusion, precisely of slow diffusion for 0<α<1, and of fast diffusion for 1<α<2.


## 5.3. The Time-Fractional Diffusion Equation

In literature there exist two forms of the time-fractional diffusion equation of a single order, one with Riemann-Liouvile derivative and one with Caputo derivative. These forms are equivalent if we refer to the standard initial conditionu(x,0)=u0(x), as shown in [73].Taking a real numberβ∈(0,1), the time-fractional diffusion equation of order β in the Riemann-Liouville sense reads(5.13)∂u∂t=KβDt1-β∂2u∂x2,
whereas in the Caputo sense reads(5.14)*Dtβu=Kβ∂2u∂x2,
where Kβ is a sort of fractional diffusion coefficient of dimensions [Kβ]=[L]2[T]-β=cm2/secβ. Like for diffusion equations of integer order (5.1) and (5.7)-(5.9), we consider the equivalent integral equation corresponding to our fractional diffusion equations (5.13)-(5.14) as(5.15)u(x,t)=u0(x)+Kβ1Γ(β)∫0t(t-τ)β-1∂2u(x,τ)∂x2dτ.
The Green function 𝒢β(x,t) for the equivalent (5.13)–(5.15) can be expressed, also in this case, in terms of the M-Wright function, as shown in Appendix A by adopting two different approaches, as follows:(5.16)𝒢β(x,t)=121Kβtβ/2Mβ/2(|x|Kβtβ/2).
The corresponding variance can be promptly obtained from the general formula (5.5) for the absolute moment of the M-Wright function. In fact, using (5.5) and (5.16) and after an obvious change of variable, we obtain(5.17)σβ2(t):=∫-∞+∞x2𝒢β(x,t)dx=2Γ(β+1)Kβtβ.
As a consequence, for 0<β<1 the variance is consistent with a process of slow diffusion with similarity exponent H=β/2. For further reading on time-fractional diffusion equations and their solutions the reader is referred, among others, to [38, 51, 53] and [74, 75].


## 5.4. The Stretched Time-Fractional Diffusion Equation

In the fractional diffusion equation (5.13), let us stretch the time variable by replacing t with tα/β where 0<α<2 and 0<β≤1. We have(5.18)∂u∂tα/β=KαβDtα/β1-β∂2u∂x2,
namely,(5.19)∂u∂t=αβtα/β-1KαβDtα/β1-β∂2u∂x2,
where Kαβ is a sort of stretched diffusion coefficient of dimensions [Kαβ]=[L]2[T]-α=cm2/secα that reduces to Kα if β=1 and to Kβ if α=β. Integration of (5.19) gives the corresponding integral equation [21](5.20)u(x,t)=u0(x)+Kαβ1Γ(β)αβ∫0tτα/β-1(tα/β-τα/β)β-1∂2u(x,τ)∂x2dτ,
whose Green function 𝒢αβ(x,t) is(5.21)𝒢αβ(x,t)=121Kαβtα/2Mβ/2(|x|Kαβtα/2),
with variance(5.22)σα,β2(t):=∫-∞+∞x2𝒢α,β(x,t)dx=2Γ(β+1)Kαβtα.
As a consequence, the resulting process turns out to be self-similar with Hurst exponent H=α/2 and a variance law consistent with both slow diffusion if 0<α<1 and fast diffusion if 1<α<2. We note that the parameter β does explicitly enter in the variance law (5.22) only as in the determination of the multiplicative constant.It is straightforward to note that the evolution equations of this process reduce to those for time-fractional diffusion ifα=β<1, for stretched diffusion if α≠1 and β=1, and finally for standard diffusion if α=β=1.


## 6. Fractional Diffusion Processes with Stationary Increments

We have seen that any Green function associated to the diffusion-like equations considered in the previous section can be interpreted as the time-evolving one-pointpdf of certain self-similar stochastic processes. However, in general, it is not possible to define a unique (self-similar) stochastic process because the determination of any multipoint probability distribution is required; see, for example, [22].In other words, starting from a master equation which describes the dynamic evolution of a probability density functionf(x,t), it is always possible to define an equivalence class of stochastic processes with the same marginal density function f(x,t). All these processes provide suitable stochastic representations for the starting equation. It is clear that additional requirements may be stated in order to uniquely select the probabilistic model.For instance, considering (5.19), the additional requirement of stationary increments, as shown by Mura et al.; see [19–22], can lead to a class {Bα,β(t),t≥0}; called “generalized” grey Brownian motion (ggBm), which, by construction, is made up of self-similar processes with stationary increments and Hurst exponent H=α/2. Thus {Bα,β(t),t≥0} is a special class of H-sssi processes, which provide non-Markovian stochastic models for anomalous diffusion, of both slow type (0<α<1) and fast type (1<α<2). According to a common terminology, H-sssi stands for H-self-similar-stationary-increments, see for details [2].TheggBm includes some well known processes, so that it defines an interesting general theoretical framework. The fractional Brownian motion (fBm) appears for β=1 and is associated with (5.7); the grey Brownian motion (gBm), defined by Schneider [8, 9], corresponds to the choice α=β, with 0<β<1, and is associated to (5.13), (5.14), or (5.15); finally, the standard Brownian motion (Bm) is recovered by setting α=β=1 being associated to (5.1). We should note that only in the particular case of Bm the corresponding process is Markovian.In Figure3 we present a diagram that allows to identify the elements of the ggBm class. The top region 1<α<2 corresponds to the domain of fast diffusion with long-range dependence. We remind that a self-similar process with stationary increments is said to possess long-range dependence if the autocorrelation function of the increments tends to zero like a power function and such that it does not result integrable, see for details [2]. In this domain the increments of the process Bα,β(t) are positively correlated, so that the trajectories tend to be more regular (persistent). It should be noted that long-range dependence is associated to a non-Markovian process which exhibits long-memory properties. The horizontal line α=1 corresponds to processes with uncorrelated increments, which model various phenomena of normal diffusion. For α=β=1 we recover the Gaussian process of the standard Brownian motion. The Gaussian process of the fractional Brownian motion is identified by the vertical line β=1. The bottom region 0<α<1 corresponds to the domain of slow diffusion. The increments of the corresponding process Bα,β(t) turn out to be negatively correlated and this implies that the trajectories are strongly irregular (antipersistent motion); the increments form a stationary process which does not exhibit long-range dependence. Finally, the diagonal line (α=β) represents the Schneider grey Brownian motion (gBm).Figure 3
Parametric class of generalized grey Brownian motion.Here we want to define theggBm by making use of the Kolmogorov extension theorem and the properties of the M-Wright function. According to Mura and Pagnini [21], the generalized grey Brownian motion Bα,β(t) is a stochastic process defined in a certain probability space such that its finite-dimensional distributions are given by(6.1)fα,β(x1,x2,…,xn;γα,β)=(2π)-(n-1)/22Γ(1+β)ndetγα,β∫0∞1τn/2M1/2(ξτ1/2)Mβ(τ)dτ,
with(6.2)ξ=(2Γ(1+β)-1∑i,j=1nxiγα,β-1(ti,tj)xj)1/2
and covariance matrix(6.3)γα,β(ti,tj)=1Γ(1+β)(tiα+tjα-|ti-tj|α),i,j=1,…,n.
The covariance matrix (6.3) characterizes the typical dependence structure of a self-similar process with stationary increments and Hurst exponent H=α/2; see, for example, [2].Using (4.23), for n=1, (6.1) reduces to(6.4)fα,β(x,t)=14tα∫0∞𝕄1/2(|x|t-α/2,τ)𝕄β(τ,1)dτ=12t-α/2Mβ/2(|x|t-α/2).
This means that the marginal density function of the ggBm is indeed the fundamental solution (5.21) of (5.18)-(5.19) with Kαβ=1. Moreover, because M1(τ)=δ(τ-1), for β=1, putting γα,1≡γα, we have that (6.1) provides the Gaussian distribution of the fractional Brownian motion(6.5)fα,1(x1,x2,…,xn;γα,1)=(2π)-(n-1)/22detγαM1/2((2∑i,j=1nxiγα-1(ti,tj)xj)1/2),
which finally reduces to the standard Gaussian distribution of Brownian motion as α=1.By the definition used above, it is clear that, fixedβ, Bα,β(t) is characterized only by its covariance structure, as shown by Mura et al. [20, 21]. In other words, the ggBm, which is not Gaussian in general, is an example of a process defined only through its first and second moments, which indeed is a remarkable property of Gaussian processes. Consequently, the ggBm appears to be a direct generalization of Gaussian processes, in the same way as the M-Wright function is a generalization of the Gaussian function.


## 7. Concluding Discussion

In this review paper we have surveyed a quite general approach to derive models for anomalous diffusion based on a family of time-fractional diffusion equations depending on two parameters:α∈(0,2), β∈(0,1].The unifying topic of this analysis is the so-calledM-Wright function by which the fundamental solutions of these equations are expressed. Such function is shown to exhibit fundamental analytical properties that were properly used in recent papers for characterizing and simulating a general class of self-similar stochastic processes with stationary increments including fractional Brownian motion and grey Brownian motion.In this respect, theM-Wright function emerges to be a natural generalization of the Gaussian density to model diffusion processes, covering both slow and fast anomalous diffusion and including non-Markovian property. In particular, it turns out to be the main function for the special H-sssi class of stochastic processes (which are self-similar with stationary increments) governed by a master equation of fractional type.



---

*Source: 104505-2010-02-11.xml*