A Tutorial on the Basic Special Functions of Fractional Calculus
AA Tutorial on the Basic Special Functions of Fractional Calculus
FRANCESCO MAINARDIUniversity Bologna and INFNDepartment of Physics & [email protected]
Abstract:
In this tutorial survey we recall the basic properties of the special function of the Mittag-Leffler andWright type that are known to be relevant in processes dealt with the fractional calculus. We outline the majorapplications of these functions. For the Mittag-Leffler functions we analyze the Abel integral equation of the sec-ond kind and the fractional relaxation and oscillation phenomena. For the Wright functions we distinguish themin two kinds. We mainly stress the relevance of the Wright functions of the second kind in probability theorywith particular regard to the so-called M -Wright functions that generalizes the Gaussian and is related with thetime-fractional diffusion equation. Key–Words: Mittag-Leffler functions, Wright functions, Fractional Calculus, Laplace, Fourier and Mellin trans-forms, Probability theory, Stable distributions.
Paper published in WSEAS TRANSACTIONS on MATHEMATICS, Vol 19 (2020), pp. 74–98.DOI: 10.37394/23206.2020.19.8
The special functions of the Mittag-Leffer and Wrighttype in general play a very important role in the theoryof the fractional differential and integral equations.The purpose of this tutorial survey is to outline the rel-evant properties of the these functions outlining theirapplications.This work is organized as follows. In Section 2, werecall the essentials of the fractional calculus that pro-vide the necessary notions for the applications.In section 3, we start to define the Mittag-Leffler func-tions. For this purpose we introduce the Gamma func-tion and the classical Mittag-Leffler functions of oneand two parameters. Then we deal with the auxiliaryfunctions of the Mittag-Leffler type to be used in thenext sections.In Section 4 we apply the above auxiliary functionsof the Mittag-Leffler type to solve the Abel integralequations of the second kind, that are noteworthycases of Volterra integral equations.In Section 5 we finally consider the most famous ap-plications of the auxiliary functions of the Mittag-Leffler tyewpe, that is the solutions of the time frac-tional differential equations governing the phenomenaof fractional relaxation and fractional oscillationsIn Section 6 we start to define the Wright functions.For this purpose we distinguish two kinds of thesefunctions. Particular attention is devoted to two spe- cial cases of the Wright function of the second kindintroduced by Mainardi in the 1990’s in virtue of theirimportance in probability theory and for the time-fractional diffusion equations. Nowadays in the FCliterature they are referred to as the Mainardi func-tions. In contrast to the general case of the Wrightfunction, they depend just on one parameter ν ∈ [0 , . One of the Mainardi functions, known as the M -Wright function, generalizes the Gaussian func-tion and degenerates to the delta function in the limit-ing case ν = 1 .Then, in Section 7 we recall how the Mainardi func-tions are related to an important class of the proba-bility density functions (pdf’s) known as the extremalL´evy stable densities. This emphasizes the relevanceof the Mainardi functions in the probability theory in-dependently on the framework of the fractional dif-fusion equations. We present some plots of the sym-metric M -Wright function on IR for several parametervalues ν ∈ [0 , / and ν ∈ [1 / , .Finally concluding remarks are carried out in Section8 and two tutorial appendices on stable distributionsand the time fractional diffusion equation are addedfor readers’ convenience.The paper is competed with historical and biblio-graphical concerning the past approach of the authortowards the Wright functions. a r X i v : . [ m a t h . G M ] O c t The essentials of fractionalcalculus
This section is mainly based on the 1997 CISM surveyby Gorenflo and Mainardi [19].The
Riemann-Liouville fractional integral of order µ > is defined as t J µ f ( t ) := 1Γ( µ ) (cid:90) t ( t − τ ) µ − f ( τ ) dτ , (2 . where Γ( µ ) := (cid:90) ∞ e − u u µ − du , Γ( n + 1) = n ! is the Gamma function .By convention t J = I (Identity operator). We canprove semi-group property t J µ t J ν = t J ν t J µ = t J µ + ν , µ, ν ≥ . (2 . Furthermore we have for t > t J µ t γ = Γ( γ + 1)Γ( γ + 1 + µ ) t γ + µ , µ ≥ , γ > − . (2 . The fractional derivative of order µ > in the Riemann-Liouville sense is defined as the operator t D µ t D µ t J µ = I , µ > . (2 . If we take m − < µ ≤ m , with m ∈ IN we recognizefrom Eqs. (2.2) and (2.4) t D µ f ( t ) := t D m t J m − µ f ( t ) , (2 . hence, for m − < µ < m , t D µ f ( t ) = d m dt m (cid:20) m − µ ) (cid:90) t f ( τ ) dτ ( t − τ ) µ +1 − m (cid:21) , (2 . a ) and, for µ = m , t D µ f ( t ) = d m dt m f ( t ) , . (2 . b ) For completion t D = I . The semi-group property isno longer valid but for t > t D µ t γ = Γ( γ + 1)Γ( γ + 1 − µ ) t γ − µ , µ ≥ , γ > − . (2 . However, the property t D µ = t J − µ is not generallyvalid! The fractional derivative of order µ ∈ ( m − , m ] ( m ∈ IN ) in the Caputo sense is defined as the opera-tor t D µ ∗ such that, t D µ ∗ f ( t ) := t J m − µ t D m f ( t ) , (2 . hence, for f m − < µ < m , t D µ ∗ f ( t ) = 1Γ( m − µ ) (cid:90) t f ( m ) ( τ ) dτ ( t − τ ) µ +1 − m , (2 . a ) and for µ = m t D m ∗ f ( t ) =: d m dt m f ( t ) . (2 . b ) Thus, when the order is not integer the two fractionalderivatives differ in that the derivative of order m doesnot generally commute with the fractional integral.We point out that the Caputo fractional derivative sat-isfies the relevant property of being zero when appliedto a constant, and, in general, to any power function ofnon-negative integer degree less than m , if its order µ is such that m − < µ ≤ m . Gorenflo and Mainardi (1997) [19] have shown theessential relationships between the two fractionalderivatives (when both of them exist), for m − <µ < m , t D µ ∗ f ( t ) = t D µ (cid:34) f ( t ) − m − (cid:88) k =0 f ( k ) (0 + ) t k k ! (cid:35) , t D µ f ( t ) − m − (cid:88) k =0 f ( k ) (0 + ) t k − µ Γ( k − µ + 1) . (2 . In particular, if m = 1 so that < µ < , we have t D µ ∗ f ( t ) = t D µ (cid:2) f ( t ) − f (0 + ) (cid:3) , t D µ f ( t ) − f (0 + ) t − µ Γ(1 − µ ) . (2 . The
Caputo fractional derivative , represents a sortof regularization in the time origin for the
Riemann-Liouville fractional derivative .We note that for its existence all the limiting values f ( k ) (0 + ) := lim t → + f ( k ) ( t ) are required to be finite for k = 0 , , . . . . m − .We observe the different behaviour of the two frac-tional derivatives at the end points of the interval ( m − , m ) namely when the order is any positiveinteger: whereas t D µ is, with respect to its order µ , an operator continuous at any positive integer, t D µ ∗ isn operator left-continuous since lim µ → ( m − + t D µ ∗ f ( t ) = f ( m − ( t ) − f ( m − (0 + ) , lim µ → m − t D µ ∗ f ( t ) = f ( m ) ( t ) . (2 . We also note for m − < µ ≤ m , t D µ f ( t ) = t D µ g ( t ) ⇐⇒ f ( t ) = g ( t )+ m (cid:88) j =1 c j t µ − j , (2 . t D µ ∗ f ( t ) = t D µ ∗ g ( t ) ⇐⇒ f ( t ) = g ( t )+ m (cid:88) j =1 c j t m − j . (2 . In these formulae the coefficients c j are arbitrary con-stants.We point out the major utility of the Caputo fractionalderivative in treating initial-value problems for phys-ical and engineering applications where initial condi-tions are usually expressed in terms of integer-orderderivatives. This can be easily seen using the Laplacetransformation .Writing the Laplace transform of a sufficiently well-behaved function f ( t ) ( t ≥ ) as L { f ( t ); s } = (cid:101) f ( s ) := (cid:90) ∞ e − st f ( t ) dt , the known rule for the ordinary derivative of integerorder m ∈ IN is L { t D m f ( t ); s } = s m (cid:101) f ( s ) − m − (cid:88) k =0 s m − − k f ( k ) (0 + ) , where f ( k ) (0 + ) := lim t → + t D k f ( t ) . For the
Caputo derivative of order µ ∈ ( m − , m ] ( m ∈ IN ) we have L { t D µ ∗ f ( t ); s } = s µ (cid:101) f ( s ) − m − (cid:88) k =0 s µ − − k f ( k ) (0 + ) ,f ( k ) (0 + ) := lim t → + t D k f ( t ) . (2 . The corresponding rule for the
Riemann-Liouvilederivative of order µ is L { t D µt f ( t ); s } = s µ (cid:101) f ( s ) − m − (cid:88) k =0 s m − − k g ( k ) (0 + ) ,g ( k ) (0 + ) := lim t → + t D k g ( t ) , g ( t ) := t J m − µ f ( t ) . (2 . Thus the rule (2.14) is more cumbersome to be usedthan (2.13) since it requires initial values concerningan extra function g ( t ) related to the given f ( t ) througha fractional integral.However, when all the limiting values f ( k ) (0 + ) are fi-nite and the order is not integer, we can prove by thatall g ( k ) (0 + ) vanish so that the formula (2.14) simpli-fies into L { t D µ f ( t ); s } = s µ (cid:101) f ( s ) , m − < µ < m . (2 . For this proof it is sufficient to apply the Laplacetransform to Eq. (2.8), by recalling that
L { t ν ; s } = Γ( ν + 1) /s ν +1 , ν > − , (2 . and then to compare (2.13) with (2.14). We note that the Mittag-Leffler functions are presentin the Mathematics Subject Classification since theyear 2000 under the number 33E12 under recommen-dation of Prof. Gorenflo.A description of the most important properties ofthese functions (with relevant references up to thefifties) can be found in the third volume of the Bate-man Project edited by Erdelyi et al. (1955) in thechapter
XV III on Miscellaneous Functions [10].The treatises where great attention is devoted to thefunctions of the Mittag-Leffler type is that by Djr-bashian (1966) [9], unfortunately in Russian.We also recommend the classical treatise on complexfunctions by Sansone & Gerretsen (1960) [55].Nowadays the Mittag-Leffler functions are widelyused in the framework of integral and differentialequations of fractional order, as shown in all treatiseson fractional calculus.In view of its several applications in Fractional Calcu-lus the Mittag-Leffler function was referred to as theQueen Function of Fractional Calculus by Mainardi& Gorenflo (2007) [34].Finally, the functions of the Mittag-Leffler type havefound an exhaustive treatment in the treatise byGorenflo, Kilbas, Mainardi & Rogosin (2014) [15]As pioneering works of mathematical nature in thefield of fractional integral and differential equations,we like to quote Hille & Tamarkin (1930) [22] whohave provided the solution of the Abel integral equa-tion of the second kind in terms of a Mittag-Lefflerunction, and Barret (1954) [1] who has expressedthe general solution of the linear fractional differentialequation with constant coefficients in terms of Mittag-Leffler functions.As former applications in physics we like to quote thecontributions by Cole (1933) [5] in connection withnerve conduction, see also Davis (1936) [7],and byGross (1947) [21] in connection with mechanical re-laxation.Subsequently, Caputo & Mainardi (1971a), (1971b)[3, 4] have shown that Mittag-Leffler functions arepresent whenever derivatives of fractional order areintroduced in the constitutive equations of a linear vis-coelastic body. Since then, several other authors havepointed out the relevance of these functions for frac-tional viscoelastic models, see e.g.
Mainardi’s survey(1997) [31] and his (2010) book [32]. Γ( z ) The
Gamma function Γ( z ) is the most widelyused of all the special functions: it is usually dis-cussed first because it appears in almost every integralor series representation of other advanced mathemati-cal functions. The first occurrence of the gamma func-tion happens in 1729 in a correspondence between Eu-ler and Goldbach. We take as its definition the integralformula due to Legendre (1809) Γ( z ) := (cid:90) ∞ uz − − u du , R e ( z ) > . (3 . This integral representation is the most common for Γ( z ) , even if it is valid only in the right half-plane of C .The analytic continuation to the left half-plane is pos-sible in different ways. As will be shown hereafter ,the domain of analyticity D Γ of Γ( z ) turns out to be D Γ = C − { , − , − , . . . } . The most common continuation is carried out by the mixed representation due to Mittag-Leffler: and readsfor z ∈ D Γ Γ( z ) = ∞ (cid:88) n =0 ( − n n !( z + n ) + (cid:90) ∞ e − u u z − du . (3 . This representation can be obtained from the so-called
Prym’s decomposition , namely by splitting the inte-gral in (3.1) into 2 integrals, the one over the interval ≤ u ≤ which is then developed in a series, the other over the interval ≤ u ≤ ∞ , which, being uni-formly convergent inside C , provides an entire func-tion. The terms of the series (uniformly convergentinside D Γ ) provide the principal parts of Γ( z ) at thecorresponding poles z n = − n . So we recognize that Γ( z ) is analytic in the entire complex plane except atthe points z n = − n ( n = 0 , , . . . ), which turn out tobe simple poles with residues R n = ( − n /n ! . Thepoint at infinity, being an accumulation point of poles,is an essential non-isolated singularity. Thus Γ( z ) is atranscendental meromorphic function.The reciprocal of the Gamma function turns out to bean entire function. its integral representation in thecomplex plane was due to Hankel (1864) and reads z ) = 12 πi (cid:90) Ha e u u z du , z ∈ C , where Ha denotes the Hankel path defined as a con-tour that begins at u = −∞ − ia ( a > ), encirclesthe branch cut that lies along the negative real axis,and ends up at u = −∞ + ib ( b > ). Of course, thebranch cut is present when z is non-integer because u − z is a multivalued function; when z is an integer,the contour can be taken to be simply a circle aroundthe origin, described in the counterclockwise direc-tion. The Mittag-Leffler functions, that we denote by E α ( z ) , E α,β ( z ) are so named in honour of G¨ostaMittag-Leffler, the eminent Swedish mathematician,who introduced and investigated these functions in aseries of notes starting from 1903 in the frameworkof the theory of entire functions [45, 46, 47, 48].The functions are defined by the series representa-tions, convergent in the whole complex plane C forRe ( α ) > } E α ( z ) := ∞ (cid:88) n =0 z n Γ( αn + 1) ,E α,β ( z ) := ∞ (cid:88) n =0 z n Γ( αn + β ) , (3 . with β ∈ C . .Originally Mittag-Leffler assumed only the parameter α and assumed it as positive, but soon later the gen-eralization with two complex parameters was consid-ered by Wiman. [60]. In both cases the Mittag-Lefflerfunctions are entire of order / Re ( α ) . The integralrepresentation for z ∈ C introduced by Mittag-Leffleran be written as E α ( z ) = 12 πi (cid:90) Ha ζ α − e ζ ζ α − z dζ, α > . (3 . Using series representations of the Mittag-Lefflerfunctions it is easy to recognize E , ( z ) = E ( z ) = e z , E , ( z ) = e z − z ,E , ( z ) = cosh( z ) , E , ( − z ) = cos( z ) ,E , ( z ) = sinh( z ) z , E , ( − z ) = sin( z ) z , (3 . and more generally E α,β ( z ) + E α,β ( − z ) = 2 E α,β ( z ) ,E α,β ( z ) − E α,β ( − z ) = 2 z E α,α + β ( z ) . (3 . Furthermore, for α = 1 / , E / ( ± z / ) = e z (cid:104) ± z / ) (cid:105) = e z erfc ( ∓ z / ) , (3 . where erf (erfc) denotes the (complementary) errorfunction defined for z ∈ C as erf ( z ) := 2 √ π (cid:90) z e − u du , erfc ( z ) := 1 − erf ( z ) . In view of applications we introduce the followingcausal functions in time domain e α ( t ; λ ) := E α ( − λ t α ) ÷ s α − s α + λ , (3 . e α,β ( t ; λ ) := t β − E α,β ( − λ t α ) ÷ s α − β s α + λ , (3 . e α,α ( t ; λ ) := t α − E α,α ( − λ t α )= ddt e α ( − λ t α ) ÷ − λs α + λ . (3 . A function f ( t ) defined in IR + is completely monotone (CM) if ( − n f n ( t ) ≥ . The function e − t is theprototype of a CM function.For a Bernstein theorem a generic CM function reads f ( t ) = (cid:90) ∞ e − rt K ( r ) d r , K ( r ) ≥ . (3 . We have for λ > e α,β ( t ; λ ) := t β − E α,β ( − λt α ) CM iff < α ≤ β ≤ . (3 . Using the Laplace transform we can prove, followingGorenflo and Mainardi (1997) [19] that for < α < (with λ = 1 ) E α ( − t α ) (cid:39) − t α Γ( α + 1) · · · t → + ,t − α Γ(1 − α ) · · · t → + ∞ , (3 . and E α ( − t α ) = (cid:90) ∞ e − rt K α ( r ) d r (3 . with K α ( r ) = 1 π r α − sin( απ ) r α + 2 r α cos( απ ) + 1 > . (3 . In the following sections we will outline the key roleof the auxiliary functions in the treatment of integraland differential equations of fractional order, includ-ing the Abel integral equation of the second kind andthe differential equations for fractional relaxation andoscillation.Before closing this section we find it convenient toprovide the plots of the functions ψ α ( t ) = e α ( t ) := E α ( − t α ) , (3 . and φ α ( t ) = t − (1 − α ) E α,α ( − t α ) := − ddt E α ( − t α ) , (3 . for t ≥ and for some rational values of α ∈ (0 , .For the sake of visibility, for both functions we haveadopted linear and logarithmic scales. Logarithmicscales have been adopted in order to better point outtheir asymptotic behaviour for large times.It is worth noting the algebraic decay of ψ α ( t ) and φ α ( t ) ψ α ( t ) ∼ sin( απ ) π Γ( α ) t α ,φ α ( t ) ∼ sin( απ ) π Γ( α + 1) t ( α +1) , t → + ∞ . (3 . igure 1: Plots of ψ α ( t ) with α = 1 / , / , / , top: linear scales; bottom: logarithmic scales.Figure 2: Plots of φ α ( t ) with α = 1 / , / , / , top: linear scales; bottom: logarithmic scales. Let us now consider the Abel equation of the secondkind with α > , λ ∈ C : u ( t ) + λ Γ( α ) (cid:90) t u ( τ )( t − τ ) − α dτ = f ( t ) , (4 . In terms of the fractional integral operator such equa-tion reads (1 + λ J α ) u ( t ) = f ( t ) , (4 . and consequently can be formally solved as follows: u ( t ) = (1 + λJ α ) − f ( t )= (cid:32) ∞ (cid:88) n =1 ( − λ ) n J αn (cid:33) f ( t ) . (4 . Recalling the definition of the fractional integral theformal solution reads u ( t ) = f ( t ) + (cid:32) ∞ (cid:88) n =1 ( − λ ) n t αn − Γ( αn ) (cid:33) ∗ f ( t ) . (4 . Recalling the definition of the function, e α ( t ; λ ) := E α ( − λ t α ) = ∞ (cid:88) n =0 ( − λ t α ) n Γ( αn + 1) , (4 . where E α denotes the Mittag-Leffler function of or-der α , we note that for t > : ∞ (cid:88) n =1 ( − λ ) n t αn − Γ( αn ) = ddt E α ( − λt α ) = e (cid:48) α ( t ; λ ) , (4 . Finally, the solution reads u ( t ) = f ( t ) + e (cid:48) α ( t ; λ ) ∗ f ( t ) . (4 . Of course the above formal proof can be made rigor-ous. Simply observe that because of the rapid growthof the gamma function the infinite series in (4.4) and(4.6) are uniformly convergent in every bounded in-terval of the variable t so that term-wise integrationsand differentiations are allowed.However, we prefer to use the alternative technique ofLaplace transforms, which will allow us to obtain thesolution in different forms, including the result (4.7).Applying the Laplace transform to (4.1) we obtain (cid:20) λs α (cid:21) ¯ u ( s ) = ¯ f ( s ) = ⇒ ¯ u ( s ) = s α s α + λ ¯ f ( s ) . (4 .. ow, let us proceed to obtain the inverse Laplacetransform of (4.8) using the following Laplace trans-form pair e α ( t ; λ ) := E α ( − λ t α ) ÷ s α − s α + λ . (4 . We can choose three different ways to get the inverseLaplace transforms from (4.8), according to the stan-dard rules. Writing (4.8) as ¯ u ( s ) = s (cid:34) s α − s α + λ ¯ f ( s ) (cid:35) , (4 . a ) we obtain u ( t ) = ddt (cid:90) t f ( t − τ ) e α ( τ ; λ ) dτ . (4 . a ) If we write (4.8) as ¯ u ( s ) = s α − s α + λ [ s ¯ f ( s ) − f (0 + )] + f (0 + ) s α − s α + λ , (4 . b ) we obtain u ( t ) = (cid:90) t f (cid:48) ( t − τ ) e α ( τ ; λ ) dτ + f (0 + ) e α ( t ; λ ) . (4 . b ) We also note that, e α ( t ; λ ) being a function differen-tiable with respect to t with e α (0 + ; λ ) = E α (0 + ) = 1 , there exists another possibility to re-write (4.8),namely ¯ u ( s ) = (cid:34) s s α − s α + λ − (cid:35) ¯ f ( s ) + ¯ f ( s ) . (4 . c ) Then we obtain u ( t ) = (cid:90) t f ( t − τ ) e (cid:48) α ( τ ; λ ) dτ + f ( t ) , (4 . c ) in agreement with (4.7). We see that the way b ) is more restrictive than the ways a ) and c ) since it requires that f ( t ) be differentiable with L -transformable derivative. Generally speaking, we consider the following differ-ential equation of fractional order α > , for t ≥ : D α ∗ u ( t ) = D α (cid:32) u ( t ) − m − (cid:88) k =0 t k k ! u ( k ) (0 + ) (cid:33) = − u ( t ) + q ( t ) , (5 . where u = u ( t ) is the field variable and q ( t ) is a givenfunction, continuous for t ≥ . Here m is a positiveinteger uniquely defined by m − < α ≤ m , whichprovides the number of the prescribed initial values u ( k ) (0 + ) = c k , k = 0 , , , . . . , m − . In particular, we consider in detail the cases(a) fractional relaxation < α ≤ , (b) fractional oscillation < α ≤ . The application of the Laplace transform yields (cid:101) u ( s ) = m − (cid:88) k =0 c k s α − k − s α + 1 + 1 s α + 1 (cid:101) q ( s ) . (5 . Then, putting for k = 0 , , . . . , m − ,u k ( t ) := J k e α ( t ) ÷ s α − k − s α + 1 ,e α ( t ) := E α ( − t α ) ÷ s α − s α + 1 , (5 . and using u (0 + ) = 1 , we find, u ( t ) = m − (cid:88) k =0 c k u k ( t ) − (cid:90) t q ( t − τ ) u (cid:48) ( τ ) dτ . (5 . In particular, the formula (5.4) encompasses the solu-tions for α = 1 , , since α = 1 , u ( t ) = e ( t ) = exp( − t ) ,α = 2 , u ( t ) = e ( t ) = cos t, u ( t ) = J e ( t ) = sin t. When α is not integer, namely for m − < α < m , we note that m − represents the integer part of α (usually denoted by [ α ] ) and m the number of ini-tial conditions necessary and sufficient to ensure theuniqueness of the solution u ( t ) . Thus the m functions u k ( t ) = J k e α ( t ) with k = 0 , , . . . , m − representthose particular solutions of the homogeneous equa-tion which satisfy the initial conditions u ( h ) k (0 + ) = δ k h , h, k = 0 , , . . . , m − , and therefore theyrepresent the fundamental solutions of the fractionalequation (5.1), in analogy with the case α = m . Fur-thermore, the function u δ ( t ) = − u (cid:48) ( t ) = − e (cid:48) α ( t ) rep-resents the impulse-response solution .Now we derive the relevant properties of the basicfunctions e α ( t ) directly from their Laplace represen-tation for < α ≤ , e α ( t ) = 12 πi (cid:90) Br e st s α − s α + 1 ds , (5 . without detouring on the general theory of Mittag-Leffler functions in the complex plane. Here Br de-notes a Bromwich path, i.e. a line Re ( s ) = σ > andIm ( s ) running from −∞ to + ∞ .or transparency reasons, we separately discuss thecases (a) < α < and (b) < α < , recall-ing that in the limiting cases α = 1 , , we know e α ( t ) as elementary function, namely e ( t ) = e − t and e ( t ) = cos t . For α not integer the power function s α is uniquelydefined as sα = | s | α e i arg s , with − π < arg s < π , that is in the complex s -plane cut along the negativereal axis.The essential step consists in decomposing e α ( t ) intotwo parts according to e α ( t ) = f α ( t ) + g α ( t ) , as indi-cated below. In case (a) the function f α ( t ) , in case (b)the function − f α ( t ) is completely monotone ; in bothcases f α ( t ) tends to zero as t tends to infinity, fromabove in case (a), from below in case (b). The otherpart, g α ( t ) , is identically vanishing in case (a), butof oscillatory character with exponentially decreasingamplitude in case (b).For the oscillatory part we obtain via the residue the-orem of complex analysis, when < α < : g α ( t ) = 2 α e t cos ( π/α ) cos (cid:20) t sin (cid:18) πα (cid:19)(cid:21) . (5 . We note that this function exhibits oscillations withcircular frequency ω ( α ) = sin ( π/α ) and with an exponentially decaying amplitude withrate λ ( α ) = | cos ( π/α ) | = − cos ( π/α ) . For the monotonic part we obtain f α ( t ) := (cid:90) ∞ e − rt K α ( r ) dr , (5 . with K α ( r ) = − π Im (cid:32) s α − s α + 1 (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) s = r e iπ (cid:33) = 1 π r α − sin ( απ ) r α + 2 r α cos ( απ ) + 1 . (5 . This function K α ( r ) vanishes identically if α is aninteger, it is positive for all r if < α < , negativefor all r if < α < . In fact in (5.8) the denominatoris, for α not integer, always positive being > ( r α − ≥ . Hence f α ( t ) has the aforementioned monotonicityproperties, decreasing towards zero in case (a), in-creasing towards zero in case (b). We note that, in order to satisfy the initial condition e α (0 + ) = 1 , we find (cid:90) ∞ K α ( r ) dr = 1 if < α ≤ , (cid:90) ∞ K α ( r ) dr = 1 − /α if < α ≤ . In Figs. 3 and 4 we display the plots of K α ( r ) , that wedenote as the basic spectral function , for some valuesof α in the intervals (a) < α < , (b) < α < . Figure 3: Plots of the basic spectral function K α ( r ) for < α < : α = 0 . , . , . , . .. Figure 4: Plots of the basic spectral function − K α ( r ) for < α < : α = 1 . , . , . , . . . In addition to the basic fundamental solutions, u ( t ) = e α ( t ) , we need to compute the impulse-response solutions u δ ( t ) = − D e α ( t ) for cases (a)and (b) and, only in case (b), the second fundamentalsolution u ( t ) = J e α ( t ) . or this purpose we note that in general it turns outthat J k f α ( t ) = (cid:90) ∞ e − rt K kα ( r ) dr , (5 . with K kα ( r ) := ( − k r − k K α ( r )= ( − k π r α − − k sin ( απ ) r α + 2 r α cos ( απ ) + 1 , (5 . where K α ( r ) = K α ( r ) , and J k g α ( t ) = 2 α e t cos ( π/α ) cos (cid:20) t sin (cid:18) πα (cid:19) − k πα (cid:21) . (5 . For the impulse-response solution we note that the ef-fect of the differential operator D is the same as thatof the virtual operator J − . Hence the solutions for the fractional relaxation are: (a) 0 < α < ,u ( t ) = c u ( t ) + (cid:90) t q ( t − τ ) u δ ( τ ) dτ , (5 . a ) where u ( t ) = (cid:82) ∞ e − rt K α ( r ) dr ,u δ ( t ) = − (cid:82) ∞ e − rt K − α ( r ) dr , (5 . a ) with u (0 + ) = 1 , u δ (0 + ) = ∞ , and for t → ∞ u ( t ) ∼ t − α Γ(1 − α ) , u ( t ) ∼ t − α Γ(2 − α ) . (5 . a ) Hence the solutions for the fractional oscillation are: (b) 1 < α < ,u ( t ) = c u ( t ) + c u ( t ) + (cid:90) t q ( t − τ ) u δ ( τ ) dτ , (5 . b ) u ( t ) = (cid:90) ∞ e − rt K α ( r ) dr + 2 α e t cos ( π/α ) cos (cid:20) t sin (cid:18) πα (cid:19)(cid:21) ,u ( t ) = (cid:90) ∞ e − rt K α ( r ) dr + 2 α e t cos ( π/α ) cos (cid:20) t sin (cid:18) πα (cid:19) − πα (cid:21) ,u δ ( t ) = − (cid:90) ∞ e − rt K − α ( r ) dr − α e t cos ( π/α ) cos (cid:20) t sin (cid:18) πα (cid:19) + πα (cid:21) , (5 . b ) with u (0 + ) = 1 , u (cid:48) (0 + ) = 0 ,u (0 + ) = 0 , u (cid:48) (0 + ) = 1 ,u δ (0 + ) = 0 , u (cid:48) δ (0 + ) = + ∞ , and for t → ∞ u ( t ) ∼ t − α Γ(1 − α ) ,u ( t ) ∼ t − α Γ(2 − α ) ,u δ ( t ) ∼ − t − α − Γ( − α ) ; (5 . b ) In Figs. 2a and 2b we display the plots of the basicfundamental solution for the following cases, respec-tively :(a) α = 0 . , . , . , , (b) α = 1 . , . , . , , obtained from the first formula in (5.13a) and (5.13b),respectively.We now want to point out that in both the cases (a) and(b) (in which α is just not integer) i.e. for fractionalrelaxation and fractional oscillation , all the funda-mental and impulse-response solutions exhibit an al-gebraic decay as t → ∞ , as discussed above.This algebraic decay is the most important effect ofthe non-integer derivative in our equations, which dra-matically differs from the exponential decay presentin the ordinary relaxation and damped-oscillation phe-nomena.Figure 5: Plots of the basic fundamental solution u ( t ) = e α ( t ) with α = 0 . , . , . , . igure 6: Plots of the basic fundamental solution u ( t ) = e α ( t ) with α = 1 . , . , . , . We would like to remark the difference between frac-tional relaxation governed by the Mittag-Leffler typefunction E α ( − at α ) and stretched relaxation governedby a stretched exponential function exp( − bt α ) with α , a , b > for t ≥ . A common behaviour isachieved only in a restricted range ≤ t (cid:28) where E α ( − at α ) (cid:39) − a Γ( α + 1) t α = 1 − b t α (cid:39) e − b t α , b = a Γ( α + 1) . In Figs. 3a, 3b, 3c for α = 0 . , . , . we havecompared E α ( − t α ) ( full line) with its asymptotic ap-proximations exp [ − t α / Γ(1 + α )] ( dashed line) validfor short times, and t − α / Γ(1 − α ) ( dotted line) validfor long times.We have adopted log-log plots in order to betterachieve such a comparison and the transition from thestretched exponential to the inverse power-law decay.In Figs. 4a, 4b, 4c we have shown some plots of the basic fundamental solution u ( t ) = e α ( t ) for α =1 . , . , . , respectively.Here the algebraic decay of the fractional oscillationcan be recognized and compared with the two con-tributions provided by f α (monotonic behaviour, dot-ted line) and g α ( t ) (exponentially damped oscillation,dashed line) The zeros of the solutions of the fractionaloscillation
Now we find it interesting to carry out some investi-gations about the zeros of the basic fundamental so- Figure 7: Log-log plot of E α ( − t α ) for α =0 . , . , . .igure 8: Decay of the basic fundamental solution u ( t ) = e α ( t ) for α = 1 . , . , . ; full line = e α ( t ) , dashed line = g α ( t ) , dotted line = f α ( t ) . lution u ( t ) = e α ( t ) in the case (b) of fractional os-cillations. For the second fundamental solution andthe impulse-response solution the analysis of the ze-ros can be easily carried out analogously.Recalling the first equation in (5.13b), the required ze-ros of e α ( t ) are the solutions of the equation e α ( t ) = f α ( t )+ 2 α e t cos ( π/α ) cos (cid:20) t sin (cid:18) πα (cid:19)(cid:21) = 0 . (5 . We first note that the function e α ( t ) exhibits an odd number of zeros, in that e α (0) = 1 , and, for suffi-ciently large t , e α ( t ) turns out to be permanently neg-ative, as shown in (5.14b) by the sign of Γ(1 − α ) . The smallest zero lies in the first positivity intervalof cos [ t sin ( π/α )] , hence in the interval < t <π/ [2 sin ( π/α )] ; all other zeros can only lie in thesucceeding positivity intervals of cos [ t sin ( π/α )] , ineach of these two zeros are present as long as α e t cos ( π/α ) ≥ | f α ( t ) | . (5 . When t is sufficiently large the zeros are expected tobe found approximately from the equation α e t cos ( π/α ) ≈ t − α | Γ(1 − α ) | , (5 . obtained from (5.15) by ignoring the oscillation fac-tor of g α ( t ) and taking the first term in the asymp-totic expansion of f α ( t ) . As shown in the report [18]such approximate equation turns out to be useful when α → + and α → − . For α → + , only one zero is present, which is ex-pected to be very far from the origin in view of thelarge period of the function cos [ t sin ( π/α )] . In fact,since there is no zero for α = 1 , and by increasing α more and more zeros arise, we are sure that onlyone zero exists for α sufficiently close to 1. Putting α = 1 + (cid:15) the asymptotic position T ∗ of this zero canbe found from the relation (5.17) in the limit (cid:15) → + . Assuming in this limit the first-order approximation,we get T ∗ ∼ log (cid:18) (cid:15) (cid:19) , (5 . which shows that T ∗ tends to infinity slower than /(cid:15) , as (cid:15) → . For details see again the 1995 report byGorenflo & Mainardi [18].For α → − , there is an increasing number of zerosup to infinity since e ( t ) = cos t has infinitely manyzeros [in t ∗ n = ( n + 1 / π , n = 0 , , . . . ]. Puttingnow α = 2 − δ the asymptotic position T ∗ for theargest zero can be found again from (5.17) in the limit δ → + . Assuming in this limit the first-order ap-proximation, we get T ∗ ∼ π δ log (cid:18) δ (cid:19) . (5 . Now, for δ → + the length of the positivity intervalsof g α ( t ) tends to π and, as long as t ≤ T ∗ , thereare two zeros in each positivity interval. Hence, in thelimit δ → + , there is in average one zero per intervalof length π , so we expect that N ∗ ∼ T ∗ /π . Remark : For the above considerations we got inspira-tion from an interesting paper by Wiman (1905) [61],who at the beginning of the XX-th century, after hav-ing treated the Mittag-Leffler function in the complexplane, considered the position of the zeros of the func-tion on the negative real axis (without providing anydetail). The expressions of T ∗ are in disagreementwith those by Wiman for numerical factors; however,the results of our numerical studies carried out in the1995 report [18] confirm and illustrate the validity ofthe present analysis.Here, we find it interesting to analyse the phenomenonof the transition of the (odd) number of zeros as . ≤ α ≤ . . For this purpose, in Table I we report theintervals of amplitude ∆ α = 0 . where these tran-sitions occur, and the location T ∗ (evaluated within arelative error of . ) of the largest zeros found atthe two extreme values of the above intervals.We recognize that the transition from 1 to 3 zeros oc-curs as . ≤ α ≤ . , that one from 3 to 5 ze-ros occurs as . ≤ α ≤ . , and so on. The lasttransition is from 15 to 17 zeros, and it just occurs as . ≤ α ≤ . .N ∗ α T ∗ ÷ . ÷ .
41 1 . ÷ . ÷ . ÷ .
57 8 . ÷ . ÷ . ÷ .
65 14 . ÷ . ÷ . ÷ .
70 20 . ÷ . ÷
11 1 . ÷ .
73 27 . ÷ . ÷
13 1 . ÷ .
76 33 . ÷ . ÷
15 1 . ÷ .
79 39 . ÷ . ÷
17 1 . ÷ .
80 45 . ÷ . Table I N ∗ = number of zeros, α = fractional order T ∗ location of the largest zero. The classical
Wright function , that we denote by W λ,µ ( z ) , is defined by the series representation con-vergent on the whole complex plane C , W λ,µ ( z ) = ∞ (cid:88) n =0 z n n !Γ( λn + µ ) , λ > − , µ ∈ C . (6 . One of its integral representations for λ > − , µ ∈ C reads as: W λ,µ ( z ) = 12 πi (cid:90) Ha e σ + zσ − λ dσσ µ , (6 . where, as usual, Ha denotes the Hankel path. Then, W λ,µ ( z ) is an entire function for all λ ∈ ( − , + ∞ ) .Originally, in 1930’s Wright assumed λ ≥ in con-nection with his investigations on the asymptotic the-ory of partitions [64, 65], and only in 1940 [66] heconsidered − < λ < .We note that in the Vol 3, Chapter 18 of the handbookof the Bateman Project [10], presumably for a mis-print, the parameter λ is restricted to be non-negative,whereas the Wright functions remained practically ig-nored in other handbooks. In 1990’s Mainardi, beingaware only of the Bateman handbook, proved that theWright function is entire also for − < λ < in hisapproaches to the time fractional diffusion equation,see [28, 29, 30].In view of the asymptotic representation in the com-plex domain and of the Laplace transform for positiveargument z = r > ( r can denote the time variable t or the positive space variable x ) the Wright functionsare distinguished in first kind ( λ ≥ ) and second kind ( − < λ < ) as outlined in the Appendix F of thebook by Mainardi [32].It is possible to prove that the Wright function is en-tire of order / (1 + λ ) , hence of exponential typeif λ ≥ . , that is only for the Wright functionsof the first kind. The case λ = 0 is trivial since W ,µ ( z ) = e z / Γ( µ ) . Recurrence relations
Some of the properties, that the Wright functionsshare with the most popular Bessel functions, wereenumerated by Wright himself.Hereafter, we quote some relevant relations from thehandbook of Bateman Project Handbook [10]: λz W λ,λ + µ ( z ) = W λ,µ − ( z )+(1 − µ ) W λ,µ ( z ) , (6 . ddz W λ,µ ( z ) = W λ,λ + µ ( z ) . (6 . e note that these relations can easily be derived fromthe series or integral representations, (6.1) or (6.2). Generalization of the Bessel functions.
For λ = 1 and µ = ν +1 ≥ the Wright functions (ofthe first kind) turn out to be related to the well knownBessel functions J ν and I ν by the identities: J ν ( z ) = (cid:18) z (cid:19) ν W ,ν +1 (cid:32) − z (cid:33) ,I ν ( z ) = (cid:18) z (cid:19) ν W ,ν +1 (cid:32) z (cid:33) . (6 . In view of this property some authors refer to theWright function as the
Wright generalized Besselfunction (misnamed also as the
Bessel-Maitland func-tion ) and introduce the notation J ( λ ) ν ( z ) := (cid:18) z (cid:19) ν W λ,ν +1 (cid:32) − z (cid:33) = (cid:18) z (cid:19) ν ∞ (cid:88) n =0 ( − n ( z/ n n ! Γ( λn + ν + 1) , (6 . with λ > and ν > − . In particular J (1) ν ( z ) := J ν ( z ) . As a matter of fact, the Wright functions (ofthe first kind) appear as the natural generalization ofthe entire functions known as Bessel - Clifford func-tions , see e.g.
Kiryakova [23], and referred by Tri-comi [59] as the uniform Bessel functions , see alsoGatteschi [13].Similarly we can properly define I ( λ ) ν ( z ) . We note that two particular Wright functions of thesecond kind, were introduced by Mainardi in 1990’s[28, 29, 30] named F ν ( z ) and M ν ( z ) ( < ν < ),called auxiliary functions in virtue of their role in thetime fractional diffusion equations. These functionsare indeed special cases of the Wright function of thesecond kind W λ,µ ( z ) by setting, respectively, λ = − ν and µ = 0 or µ = 1 − ν . Hence we have: F ν ( z ) := W − ν, ( − z ) , < ν < , (6 . and M ν ( z ) := W − ν, − ν ( − z ) , < ν < , ((6 . These functions are interrelated through the followingrelation: F ν ( z ) = νzM ν ( z ) . (6 . The series and integral representations of the auxiliary functions are derived from those of the general Wrightfunctions. Then for z ∈ C and < ν < we have: F ν ( z ) = ∞ (cid:88) n =1 ( − z ) n n !Γ( − νn )= − π ∞ (cid:88) n =1 ( − z ) n n ! Γ( νn + 1) sin ( πνn ) , (6 . and M ν ( z )= ∞ (cid:88) n =0 ( − z ) n n !Γ[ − νn + (1 − ν )]= 1 π ∞ (cid:88) n =1 ( − z ) n − ( n − νn ) sin ( πνn ) , (6 . The second series representations in (6.10)-(6.11)have been obtained by using the well-known reflec-tion formula for the Gamma function Γ( ζ ) Γ(1 − ζ ) = π/ sin πζ . For the integral representation we have F ν ( z ) := 12 πi (cid:90) Ha e σ − zσ ν dσ, (6 . and M ν ( z ) := 12 πi (cid:90) Ha e σ − zσ ν dσσ − ν . (6 . As usual, the equivalence of the series and integralrepresentations is easily proved using the Hankel for-mula for the Gamma function and performing a term-by-term integration.Explicit expressions of F ν ( z ) and M ν ( z ) in terms ofknown functions are expected for some particular val-ues of ν as shown and recalled in [28, 29, 30], thatis M / ( z ) = 1 √ π e − z / , (6 . M / ( z ) = 3 / Ai ( z/ / ) , . (6 . Liemert and Klenie [24] have added the following ex-pression for ν = 2 / M / ( z ) = 3 − / e − z / (cid:104) / z Ai (cid:16) z / / (cid:17) − (cid:48) (cid:16) z / / (cid:17)(cid:105) (6 . where Ai and Ai (cid:48) denote the Airy function and itsfirst derivative. Furthermore they have suggested in In the RHS of Eq (6.10) we have corrected an error alreadypresent since the Appendix F of my 2010 book [32]. he positive real field IR + the following remarkablyintegral representation M ν ( x ) = 1 π x ν/ (1 − ν ) − ν · (cid:90) π C ν ( φ ) exp ( − C ν ( φ )) x / (1 − ν ) dφ, (6 . where C ν ( φ ) = sin(1 − ν )sin φ (cid:18) sin νφ sin φ (cid:19) ν/ (1 − ν ) , (6 . corresponding to equation (7) of the article written bySaa and Venegeroles [54] .Furthermore, it can be proved, see [42] that M /q ( z ) satisfies the differential equation of order q − d q − dz q − M /q ( z ) + ( − q q z M /q ( z ) = 0 , (6 . subjected to the q − initial conditions at z = 0 , de-rived from (6.15), M ( h )1 /q (0) = ( − h π Γ[( h + 1) /q ] sin[ π ( h + 1) /q ] , (6 . with h = 0 , , . . . q − . We note that, for q ≥ , Eq.(6.18) is akin to the hyper-Airy differential equationof order q − , see e.g. [Bender & Orszag 1987].We find it convenient to show the plots of the M -Wright functions on a space symmetric interval of IRin Figs 1, 2, corresponding to the cases ≤ ν ≤ / and / ≤ ν ≤ , respectively. We recognizethe non-negativity of the M -Wright function on IRfor / ≤ ν ≤ consistently with the analysison distribution of zeros and asymptotics of Wrightfunctions carried out by Luchko, see [25], [26].Figure 9: Plots of the M -Wright function as a functionof the x variable, for ≤ ν ≤ / . Figure 10: Plots of the M -Wright function as a func-tion of the x variable, for / ≤ ν ≤ . Let us consider the Laplace transform of the Wrightfunction using the usual notation W λ,µ ( ± r ) ÷ (cid:90) ∞ e − s r W λ,µ ( ± r ) dr , where r denotes a non negative real variable and s isthe Laplace complex parameter.When λ > the series representation of the Wrightfunction can be transformed term-by-term. In fact, fora known theorem of the theory of the Laplace trans-forms, see e.g. Doetsch (194) [8], the Laplace trans-form of an entire function of exponential type canbe obtained by transforming term-by-term the Taylorexpansion of the original function around the origin.In this case the resulting Laplace transform turns outto be analytic and vanishing at infinity. As a con-sequence, we obtain the Laplace transform pair for | s | > W λ,µ ( ± r ) ÷ s E λ,µ (cid:18) ± s (cid:19) , λ > , (6 . where E λ,µ denotes the Mittag-Leffler function in twoparameters. The proof is straightforward noting that ∞ (cid:88) n =0 ( ± r ) n n ! Γ( λn + µ ) ÷ s ∞ (cid:88) n =0 ( ± /s ) n Γ( λn + µ ) , and recalling the series representation of the Mittag-Leffler function, E α,β ( z ) := ∞ (cid:88) n =0 z n Γ( αn + β ) , α > , β ∈ C . or λ → + Eq. (6.20) provides the Laplace trans-form pair for | s | > , W + ,µ ( ± r ) = e ± r Γ( µ ) ÷ µ ) 1 s ∓ s E ,µ (cid:18) ± s (cid:19) , (6 . where, to remain in agreement with (6.20), we haveformally put, E ,µ ( z ) := ∞ (cid:88) n =0 z n Γ( µ ) := 1Γ( µ ) E ( z ) := 1Γ( µ ) 11 − z . We recognize that in this special case the Laplacetransform exhibits a simple pole at s = ± while for λ > it exhibits an essential singularity at s = 0 . For − < λ < the Wright function turns out to bean entire function of order greater than 1, so that careis required in establishing the existence of its Laplacetransform, which necessarily must tend to zero as s →∞ in its half-plane of convergence.For the sake of convenience we limit ourselves toderive the Laplace transform for the special case of M ν ( r ) ; the exponential decay as r → ∞ of the orig-inal function provided by (6.20) ensures the existenceof the image function. From the integral representa-tion (6.13) of the M ν function we obtain M ν ( r ) ÷ πi (cid:90) ∞ e − s r (cid:20)(cid:90) Ha e σ − rσ ν dσσ − ν (cid:21) dr = 12 πi (cid:90) Ha e σ σ ν − (cid:20)(cid:90) ∞ e − r ( s + σ ν ) dr (cid:21) dσ = 12 πi (cid:90) Ha e σ σ ν − σ ν + s dσ . Then, by recalling the integral representation of theMittag-Leffler function (3.4), E α ( z ) = 12 πi (cid:90) Ha ζ α − e ζ ζ α − z dζ , α > , z ∈ C , we obtain the Laplace transform pair M ν ( r ) := W − ν, − ν ( − r ) ÷ E ν ( − s ) , < ν < , . (6 . In this case, transforming term-by-term the Taylor se-ries of M ν ( r ) yields a series of negative powers of s , that represents the asymptotic expansion of E ν ( − s ) as s → ∞ in a sector around the positive real axis. We note that (6.22) contains the well-known Laplacetransform pair, see e.g. Doetsch [8] and Eq. (3.7): M / ( r ) := 1 √ π exp (cid:16) − r / (cid:17) ÷ E / ( − s ) = exp (cid:0) s (cid:1) erfc ( s ) , (6 . valid ∀ s ∈ C Analogously, using the more general integral repre-sentation (6.2) of the standard Wright function, wecan prove that in the case λ = − ν ∈ ( − , andRe ( µ ) > , we get W − ν,µ ( − r ) ÷ E ν,µ + ν ( − s ) , < ν < . (6 . In the limit as ν → + (thus λ → − ) we formallyobtain the Laplace transform pair W − ,µ ( − r ) := e − r Γ( µ ) ÷ µ ) 1 s + 1 := E ,µ ( − s ) (6 . Therefore, as λ → ± , and µ = 1 we note a sort ofcontinuity in the results (6.21) and (6.25) since W , ( − r ) := e − r ÷ s + 1) (6 . with s + 1) = (cid:40) (1 /s ) E ( − /s ) , | s | > E ( − s ) , | s | < . (6 . We here point out the relevant Laplace transform pairsrelated to the auxiliary functions of argument r − ν with < ν < , see for details the cited author’spapers r F ν (1 /r ν ) = νr ν +1 M ν (1 /r ν ) ÷ e − s ν , (6 . ν F ν (1 /r ν ) = 1 r ν M ν (1 /r ν ) ÷ e − s ν s − ν . (6 . We recall that the Laplace transform pairs in (6.28)were formerly considered by Pollard (1946) [51],Later Mikusinski (1959 )[44] got a similar resultbased on his theory of operational calculus, and fi-nally, albeit unaware of the previous results, Buchen& Mainardi (1975) [2] derived the result in a formalway. We note, however, that all these Authors werenot informed about the Wright functions. Aware ofthe Wright functions was Stankovic [58] who in 1970gave a rigorous proof of the Laplace transform pairsnvolving the Wright functions with first negative pa-rameter, here referred of the second kind,Hereafter we like to provide two independent proofsof (6.28) carrying out the inversion of exp( − s ν ) , ei-ther by the complex Bromwich integral formula or bythe formal series method. Similarly we can act for theLaplace transform pair (6.29).For the complex integral approach we deform theBromwich path Br into the Hankel path Ha , that isequivalent to the original path, and we set σ = sr .Recalling (6.13)-(6.14), we get L − [exp ( − s ν )] = 12 πi (cid:90) Br e sr − s ν ds = 12 πi r (cid:90) Ha e σ − ( σ/r ) ν dσ = 1 r F ν (1 /r ν ) = νr ν +1 M ν (1 /r ν ) . Expanding in power series the Laplace transform andinverting term by term we formally get, after recalling(6.12)-(6.13): L − [exp ( − s ν )] = ∞ (cid:88) n =0 ( − n n ! L − [ s νn ]= ∞ (cid:88) n =1 ( − n n ! r − νn − Γ( − νn )= 1 r F ν (1 /r ν ) = νr ν +1 M ν (1 /r ν ) . We note the relevance of Laplace transforms (6.24)and (6.28) in pointing out the non-negativity of theWright function M ν ( x ) for x > and the completemonotonicity of the Mittag-leffler functions E ν ( − x ) for x > and < ν < . In fact, since exp ( − s ν ) denotes the Laplace transform of a probability density(precisely, the extremal L´evy stable density of index ν , see [Feller (1971)]), the L.H.S. of (6.28) must benon-negative, and so also the L.H.S of F(24). As amatter of fact the Laplace transform pair (6.24) shows,replacing s by x , that the spectral representation of theMittag-Leffler function E ν ( − x ) is expressed in termsof the M -Wright function M ν ( r ) , that is for x ≥ E ν ( − x ) = (cid:90) ∞ e − rx M ν ( r ) dr, < ν < . (6 . We now recognize that Eq. (6.30) is consistent with aresult derived by Pollard (1948) [52]. It is instructive to compare the spectral representationof E ν ( − x ) with that of the function E ν ( − t ν ) . Werecall for t ≥ , E ν ( − t ν ) = (cid:90) ∞ e − rt K ν ( r ) dr, < ν < , (6 . with spectral function K ν ( r ) = 1 π r ν − sin( νπ ) r ν + 2 r ν cos ( νπ ) + 1= 1 π r sin( νπ ) r ν + r − ν + 2 cos ( νπ ) . (6 . The relationship between M ν ( r ) and K ν ( r ) is worthto be explored. Both functions are non-negative, inte-grable and normalized in IR + , so they can be adoptedin probability theory as density functions.The transition K ν ( r ) → δ ( r − for ν → is easy tobe detected numerically in view of the explicit repre-sentation (6.32). On the contrary, the analogous tran-sition M ν ( r ) → δ ( r − is quite a difficult matter inview of its series and integral representations. In thisrespect see the figure hereafter carried out in the 1997paper by Mainardi and Tomirotti [43].Figure 11: Plots of M ν ( r ) with ν = 1 − (cid:15) around themaximum r ≈ .Here we have compared the cases (a) (cid:15) = 0 . , (b) (cid:15) = 0 . , obtained by an asymptotic method origi-nally due to Pipkin (continuous line), 100 terms-series(dashed line) and the standard saddle-point method(dashed-dotted line).In the following Section we deal the asymptotic rep-resentations of the Wright functions for parameter λ = − ν not close to the singular case ν = 1 . For the asymptotic analysis in the whole complexplane for the Wright functions, the interested reader isreferred to Wong and Zhao (1999a),(1999b)[62, 63],who have considered the two cases λ ≥ and − < < separately, including a description of Stokes’discontinuity and its smoothing.For the Wright functions of the second kind, where λ = − ν ∈ ( − , , we recall the asymptotic ex-pansion originally obtained by Wright himself, that isvalid in a suitable sector about the negative real axisas | z | → ∞ , W − ν,µ ( z ) = Y / − µ e − Y × (cid:34) M − (cid:88) m =0 A m Y − m + O ( | Y | − M ) (cid:35) , (6 . with Y = Y ( z ) = (1 − ν ) ( − ν ν z ) / (1 − ν ) , (6 . where the A m are certain real numbers.Let us first point out the asymptotic behaviour of thefunction M ν ( r ) as r → ∞ . Choosing as a variable r/ν rather than r , the computation of the requestedasymptotic representation by the saddle-point approx-imation yields, see Mainardi & Tomirotti (1994) [42], M ν ( r/ν ) ∼ a ( ν ) r ( ν − / / (1 − ν ) × exp (cid:20) − b ( ν ) r / (1 − ν ) (cid:21) , (6 . where a ( ν ) and b ( ν ) are positive coefficients a ( ν ) = 1 (cid:112) π (1 − ν ) , b ( ν ) = 1 − νν > . (6 . The above evaluation is consistent with the first termin Wright’s asymptotic expansion (6.33) after havingused the definition (6.36).We point out that in the limit ν → − the function M ν ( r ) tends to the Dirac function δ ( r − , but in anon-symmetric way as shown in the two plots in figure11 of the previous subsection. Using the known completely monotone functions, thetechnique of the Laplace transform, and the Bern-stein theorem, one can prove non-negativity of someWright functions. Say, the function p ν,µ ( r ) = Γ( µ ) W − ν,µ − ν ( − r ) (7 . can be interpreted as a one-sided probability densityfunction (pdf) for < ν ≤ , ν ≤ µ (see [27]). To show this, we use the Laplace transform pair (6.24)that we rewrite in the form W − ν,µ − ν ( − r ) ÷ E ν,µ ( − s ) , < ν ≤ , (7 . and the fact that the Mittag-Leffler function E ν,µ ( − s ) is completely monotone for < ν ≤ , ν ≤ µ . Ac-cording to the Bernstein theorem, the function p ν,µ ( r ) is non-negative.To calculate the integral of p ν,µ ( r ) over IR + let usmention that it can be interpreted as the Laplace trans-form of p ν,µ at the point s = 0 or the Mellin transformat s = 1 . Using the Mellin integral transform of theWright function as in [26] leads now to the followingchain of equalities: (cid:90) ∞ p ν,µ ( r ) dr = (cid:90) ∞ Γ( µ ) W − ν,µ − ν ( − r ) dr = Γ( µ )Γ( s )Γ( µ − ν + νs ) (cid:12)(cid:12)(cid:12)(cid:12) s =1 = Γ( µ )Γ( µ ) = 1 . The Mellin transform technique allows us to calculatealso all moments of order s > of the pdf p ν,µ ( r ) onIR + : (cid:90) ∞ p ν,µ ( r ) r s dr = (cid:90) ∞ Γ( µ ) W − ν,µ − ν ( − r ) r s +1 − dr = Γ( µ )Γ( s + 1)Γ( µ + νs ) . (7 . For µ = 1 , the pdf p ν,µ ( r ) can be expressed in termsof the M -Wright function M ν ( r ) , < ν < definedby Eq. (6.8). As it is well known (see, e.g., [32]), M ν ( r ) can be interpreted as a one-sided pdf on IR + with the moments given by the formula (cid:90) ∞ M ν ( r ) r s dr = Γ( s + 1)Γ(1 + νs ) , s > . (7 . We find it worthwhile to recall the relations be-tween the Mainardi auxilary functions and the ex-tremal L´evy stable densities as proven in the 1997 pa-per by Mainardi and Tomirotti [43]. For an essentialaccount of the general theory of L´evy stable distribu-tions in probability In the present paper the interestedreader may consult the Appendix A in the present pa-per. More details can be found in the 1997 E-print byMainardi-Gorenflo-Paradisi [41] and in the 2001 pa-per by Mainardi-Luchko-Pagnini. [36], recalled in theAppendix F in [32].Then, from a comparison between the series expan-sions of stable densities according to the Fekller-Takayasu canonic form with index of stability α ∈ , and skewness θ ( | θ | ≤ min { α, − α } ), andthose of the Mainardi auxiliary functions in Eqs. (6.7)- (6.8), we recognize, see also [33], that the auxiliaryfunctions are related to the extremal stable densitiesas follows L − αα ( x ) = 1 x F α ( x − α ) = αx α +1 M α ( x − α ) , < α < , x ≥ . (7 . L α − α ( x ) = 1 x F /α ( x ) = 1 α M /α ( x ) , < α ≤ , −∞ < x < + ∞ . (7 . In the above equations, for α = 1 , the skewness pa-rameter turns out to be θ = − , so we get the singularlimit L − ( x ) = M ( x ) = δ ( x − . (7 . Hereafter we show the plots the extremal stable densi-ties according to their expressions in terms of the M -Wright functions, see Eq. (7.1), Eq. (7.1) for α = 1 / and α = 3 / , respectively.Figure 12: Plot of the unilateral extremal stable pdffor α = 1 / Figure 13: Plot of the bilateral extremal stable pdf for α = 3 / We point out that the most relevant applications of ourauxiliary functions, are when the variable is real. Inparticular we consider the case of the symmetric M -Wright function as a function of the variable | x | forall IR with varying its parameter ν ∈ [0 , becauserelated to the fundamental solution of the Cauchyproblem of the time fractional diffusion-wave equa-tion dealt in Appendix B In the following Figs. 14and 15 we compare the plots of the M ν ( | x | ) -Wrightfunctions in | x | ≤ for some rational values in theranges ν ∈ [0 , / and ν ∈ [1 / , , respectively.To gain more insight of the effect of the variationof the parameter ν we will adopt bot linear and log-arithmic scales for the ordinate. Thus in Fig. 14we see the transition from exp( −| x | ) for ν = 0 to / √ π exp( − x ) for ν = 1 / , whereas in Fig. 15 wesee the transition from / √ π exp( − x ) to the deltafunctions δ ( x ± for ν = 1 . In plotting M ν ( | x | ) atfixed ν for sufficiently large | x | the asymptotic repre-sentation (6.34)-(6.35) is useful since, as | x | increases,the numerical convergence of the series in (6.11) be-comes poor and poor up to being completely ineffi-cient: henceforth, the matching between the series andthe asymptotic representation is relevant. However, as ν → − , the plotting remains a very difficult task be-cause of the high peak arising around x = ± .Figure 14: Plots of M ν ( | x | ) with ν =0 , / , / , / , / for | x | ≤ ; top: linearscale, bottom: logarithmic scale.igure 15: Plots of M ν ( | x | ) with ν =1 / , / , / , for | x | ≤ : top: linear scale;bottom: logarithmic scale) The Fourier transform of the M -Wright function. The Fourier transform of the symmetric (and normal-ized) M -Wright function provides its characteristicfunction useful in Probability theory. F (cid:104) M ν ( | x | ) (cid:105) ≡ (cid:100) M ν ( | x | ):= 12 (cid:90) + ∞−∞ e iκx M ν ( | x | ) dx = (cid:90) ∞ cos( κx ) M ν ( x ) dx = E ν ( − κ ) . (7 . For this prove it is sufficient to develop in series thecosine function and use the formula for the absolutemoments of the M -Wright function in IR + . (cid:90) ∞ cos( κx ) M ν ( x ) dx = ∞ (cid:88) n =0 ( − n κ n (2 n )! (cid:90) ∞ x n M ν ( x ) dx = ∞ (cid:88) n =0 ( − n κ n Γ(2 νn + 1) = E ν, ( − κ ) . We also have (cid:90) ∞ sin( κx ) M ν ( x ) dx = ∞ (cid:88) n =0 ( − n κ n +1 (2 n + 1)! (cid:90) ∞ x n +1 M ν ( x ) dx = ∞ (cid:88) n =0 ( − n κ n +1 Γ(2 νn + 1 + ν ) = κE ν, ν ( − κ ) . We now consider M -Wright functions as spatial prob-ability densities evolving in time with self-similarity,that is M ν ( x, t ) := t − ν M ν ( xt − ν ) , x, t ≥ . (7 . These M -Wright functions are relevant for their com-position rules proved by Mainardi et al. in [36], andmore generally in [39] by using the Mellin Trans-forms.The main statement can be summarized as: Let M λ ( x ; t ) , M µ ( x ; t ) and M ν ( x ; t ) be M -Wrightfunctions of orders λ, µ, ν ∈ (0 , respectively,then the following composition formula holds for any x, t ≥ : M ν ( x, t ) = (cid:90) ∞ M λ ( x ; τ ) M µ ( τ ; t ) dτ, ν = λ µ. (7 . The above equation is also intended as a subordina-tion formula because it can be used to define subordi-nation among self-similar stochastic processes (withindependent increments), that properly generalize themost popular Gaussian processes, to which they re-duce for ν = 1 / .These more general processes are governed by time-fractional diffusion equations, as shown in papers ofour research group, see the survey paper by Mainardi,Mura and Pagnini [37] and references therein. Thesegeneral processes are referred to as Generalized greyBrownian Motions , that include both
Gaussian Pro-cesses (standard Brownian motion, fractional Brown-ian motion) and non-Gaussian Processes (Schneider’sgrey Brownian motion), to which the interested readeris referred for details.
In this survey we have outlined the basic proper-ties of the Mittag-Leffler and Wright functions. Wehave also considered a number of applications tak-ing into account special functions of these families.We have stressed their relations with fractional cal-culus. In particular, we have added a number of tu-torial appendices to enlarge the fields of applicabilityof the Wright functions, nowadays less known thanthe Mittag-Leffler functions to which they are relatedthrough integral transforms. cknowledgments
The research activity of the author is carried out in theframework of the activities of the National Group ofMathematical Physics (GNFM, INdAM), Italy.
Appendix A: L´evy stable distributions
The term stable has been assigned by the French math-ematician Paul L´evy, who in the 1920’s years starteda systematic research in order to generalize the cele-brated
Central Limit Theorem to probability distribu-tions with infinite variance. For stable distributionswe can assume the followingD
EFINITION : If two independent real random vari-ables with the same shape or type of distribution arecombined linearly and the distribution of the result-ing random variable has the same shape, the commondistribution (or its type, more precisely) is said to bestable .The restrictive condition of stability enabled L´evy(and then other authors) to derive the canonical form for the Fourier transform of the densities of these dis-tributions. Such transform in probability theory isknown as characteristic function .Here we follow the parameterization adopted in Feller(1971) [11] revisited in 1998 by Gorenflo & Mainardi[20] and popularized in 2001 by Mainardi, Luchko &Pagnini [36].Denoting by L θα ( x ) a (strictly) stable density in IR,where α is the index of stability and an θ the asymme-try parameter, improperly called skewness , its charac-teristic function reads: L θα ( x ) ÷ (cid:98) L θα ( κ ) = exp (cid:104) − ψ θα ( κ ) (cid:105) ,ψ θα ( κ ) = | κ | α e i ( sign κ ) θπ/ , ( A. with < α ≤ , | θ | ≤ min { α, − α } . ( A. . We note that the allowed region for the parameters α and θ turns out to be a diamond in the plane { α, θ } with vertices in the points (0 , , (1 , , (1 , − , (2 , , that we call the Feller-Takayasu diamond , seeFig.16. For values of θ on the border of the diamond(that is θ = ± α if < α < , and θ = ± (2 − α ) if < α < ) we obtain the so-called extremal stabledensities . Figure 16: The Feller-Takayasu diamond for L´evy sta-ble densities.We note the symmetry relation L θα ( − x ) = L − θα ( x ) , sothat a stable density with θ = 0 is symmetricStable distributions have noteworthy properties ofwhich the interested reader can be informed from theexisting literature. Here-after we recall some peculiar Properties :- The class of stable distributions possesses its owndomain of attraction , see e.g.
Feller (1971)[11].-
Any stable density is unimodal and indeed bell-shaped , i.e. its n -th derivative has exactly n zeros inIR, see Gawronski (1984) [14] and Simon (2015) [57].- The stable distributions are self-similar and infinitelydivisible . These properties derive from the canonicalform (A.1)-(A.2) through the scaling property of theFourier transform.
Self-similarity means L θα ( x, t ) ÷ exp (cid:104) − tψ θα ( κ ) (cid:105) ⇐⇒ L θα ( x, t ) = t − /α L θα ( x/t /α )] , ( A. where t is a positive parameter. If t is time, then L θα ( x, t ) is a spatial density evolving on time withself-similarity. Infinite divisibility means that for every positive inte-ger n , the characteristic function can be expressed asthe n th power of some characteristic function, so thatany stable distribution can be expressed as the n -foldconvolution of a stable distribution of the same type.Indeed, taking in (A.3) θ = 0 , without loss of gener-lity, we have e − t | κ | α = ‘ (cid:104) e − ( t/n ) | κ | α (cid:105) n ⇐⇒ L α ( x, t ) = (cid:2) L α ( x, t/n ) (cid:3) ∗ n , ( A. where (cid:104) L α ( x, t/n ) (cid:105) ∗ n := L α ( x, t/n ) ∗ . . . ∗ L α ( x, t/n ) is the multiple Fourier convolution in IR with n iden-tical terms.Only in special cases we get well-known probabilitydistributions.For α = 2 (so θ = 0 ), we recover the Gaussian pdf ,that turns out to be the only stable density with finitevariance, and more generally with finite moments ofany order δ ≥ . In fact L ( x ) = 12 √ π e − x / . ( A. All the other stable densities have finite absolute mo-ments of order δ ∈ [ − , α ) .For α = 1 and | θ | < , we get L θ ( x ) = 1 π cos( θπ/ x + sin( θπ/ + [cos( θπ/ , ( A. which for θ = 0 includes the Cauchy-Lorentz pdf , L ( x ) = 1 π
11 + x . ( A. In the limiting cases θ = ± for α = 1 we obtain the singular Dirac pdf ’s L ± ( x ) = δ ( x ± . ( A. In general we must recall the power series expansionsprovided by Feller (1971) [11]. We restrict our atten-tion to x > since the evaluations for x < can beobtained using the symmetry relation.The convergent expansions of L θα ( x ) ( x > ) turn outto befor < α < , | θ | ≤ α : L θα ( x ) = 1 π x ∞ (cid:88) n =1 ( − x − α ) n Γ(1 + nα ) n ! sin (cid:20) nπ θ − α ) (cid:21) , ( A. for < α ≤ , | θ | ≤ − α : L θα ( x ) = 1 π x ∞ (cid:88) n =1 ( − x ) n Γ(1 + n/α ) n ! sin (cid:20) nπ α ( θ − α ) (cid:21) . ( A. From the series (A.9) and the symmetry relation wenote that the extremal stable densities for < α < are unilateral, precisely vanishing for x > if θ = α ,vanishing for x < if θ = − α . In particular theunilateral extremal densities L − αα ( x ) with < α < have as Laplace transform exp( − s α ) .From a comparison between the series expansionsin (A.9)-(A.10) and in (6.10)-(6.11) for the auxiliaryfunctions F α ( x ) , M α ( x ) we recognize that for x > the auxiliary functions of the Wright type are relatedto the extremal stable densities as in Eqs. (6.28)-(6.29).More generally, all (regular) stable densities, given inEqs. (6.38)-(6.39), were recognized to belong to theclass of Fox H -functions, as formerly shown in 1986by Schneider [56], see also Mainardi-Pagnini-Saxena(2003) [40]. Appendix B: The time-fractionaldiffusion equation
There exist three equivalent forms of the time-fractional diffusion equation of a single order, twowith fractional derivative and one with fractional in-tegral, provided we refer to the standard initial condi-tion u ( x,
0) = u ( x ) .Taking a real number β ∈ (0 , , the time-fractionaldiffusion equation of order β in the Riemann-Liouville sense reads ∂u∂t = K β D − βt ∂ u∂x , ( B. in the Caputo sense reads ∗ D βt u = K β ∂ u∂x , ( B. and in integral form u ( x, t ) = u ( x )+ K β β ) (cid:90) t ( t − τ ) β − ∂ u ( x, τ ) ∂x dτ . ( B. where K β is a sort of fractional diffusion coefficientof dimensions [ K β ] = [ L ] [ T ] − β = cm /sec β .The fundamental solution (or Green function ) G β ( x, t ) for the equivalent Eqs. (B.1) - (B.3), that is the solu-tion corresponding to the initial condition G β ( x, + ) = u ( x ) = δ ( x ) ( B. an be expressed in terms of the M -Wright function G β ( x, t ) = 12 1 (cid:112) K β t β/ M β/ (cid:32) | x | (cid:112) K β t β/ (cid:33) . ( B. The corresponding variance can be promptly obtained σ β ( t ) := (cid:90) + ∞−∞ x G β ( x, t ) dx = 2Γ( β + 1) K β t β . ( B. As a consequence, for < β < the variance is con-sistent with a process of slow diffusion with similarityexponent H = β/ .The fundamental solution G β ( x, t ) for the time-fractional diffusion equation can be obtained by ap-plying in sequence the Fourier and Laplace transformsto any form chosen among Eqs. (B.1)-(B.3). Let usdevote our attention to the integral form (B.3) usingnon-dimensional variables by setting K β = 1 andadopting the notation J βt for the fractional integral.Then, our Cauchy problem reads G β ( x, t ) = δ ( x ) + J βt ∂ G β ∂x ( x, t ) . ( B. In the Fourier-Laplace domain, after applying formulafor the Laplace transform of the fractional integral andobserving (cid:98) δ ( κ ) ≡ , we get (cid:99)(cid:102) G β ( κ, s ) = 1 s − κ s β (cid:99)(cid:102) G β ( κ, s ) , from which for Re ( s ) > , κ ∈ IR (cid:99)(cid:102) G β ( κ, s ) = s β − s β + κ , < β ≤ .. ( B. Strategy (S1):
Recalling the Fourier transform pair ab + κ ÷ a b / e −| x | b / , a, b > , ( B. and setting a = s β − , b = s β , we get (cid:102) G β ( x, s ) = 12 s β/ − e −| x | s β/ . ( B. Strategy (S2):
Recalling the Laplace transform pair s β − s β + c L ↔ E β ( − ct β ) , c > , ( B. and setting c = κ , we have (cid:99) G β ( κ, t ) = E β ( − κ t β ) . ( B. Both strategies lead to the result G β ( x, t ) = 12 M β/ ( | x | , t ) = 12 t − β/ M β/ (cid:18) | x | t β/ (cid:19) , ( B. consistent with Eq. (B.5). The time-fractional drift equation
Let us finally note that the M -Wright function doesappear also in the fundamental solution of the time-fractional drift equation. Writing this equation in non-dimensional form and adopting the Caputo derivativewe have ∗ D βt u ( x, t ) = − ∂∂x u ( x, t ) , −∞ < x < + ∞ , t ≥ , ( B. where < β < and u ( x, + ) = u ( x ) . When u ( x ) = δ ( x ) we obtain the fundamental solution(Green function) that we denote by G ∗ β ( x, t ) . Follow-ing the usual approach we show that G ∗ β ( x, t ) = t − β M β (cid:18) xt β (cid:19) , x > , , x < , ( B. that for β = 1 reduces to the right running pulse δ ( x − t ) for x > .In the Fourier-Laplace domain, after applying the for-mula for the Laplace transform of the Caputo frac-tional derivative and observing (cid:98) δ ( κ ) ≡ , we get s β (cid:99)(cid:102) G ∗ β ( κ, s ) − s β − = + iκ (cid:99)(cid:102) G ∗ β ( κ, s ) , from which, for Re ( s ) > , κ ∈ IR (cid:99)(cid:102) G ∗ β ( κ, s ) = s β − s β − iκ , < β ≤ , . ( B. To determine the Green function G ∗ β ( x, t ) in the space-time domain we can follow two alternative strategiesrelated to the order in carrying out the inversions in(B.16).(S1) : invert the Fourier transform getting (cid:102) G β ( x, s ) and then invert the remaining Laplace transform;(S2) : invert the Laplace transform getting (cid:99) G ∗ β ( κ, t ) and then invert the remaining Fourier transform. Strategy (S1):
Recalling the Fourier transform pair ab − iκ F ↔ ab e − xb , a, b > , x > , ( B. and setting a = s β − , b = s β , we get (cid:102) G ∗ β ( x, s ) = s β − e − xs β . ( B. trategy (S2): Recalling the Laplace transform pair s β − s β + c L ↔ E β ( − ct β ) , c > , ( B. and setting c = − iκ , we have (cid:99) G ∗ β ( κ, t ) = E β ( iκt β ) . ( B. Both strategies lead to the result (B.15). In view ofEq. (7.5) we also recall that the M -Wright functionis related to the unilateral extremal stable density ofindex β . Then, using our notation for stable densities,we write our Green function (B.15) as G ∗ β ( x, t ) = tβ x − − /β L − ββ (cid:16) tx − /β (cid:17) . ( B. Appendix C: Historical andbibliographic notes
In the early nineties, precisely in his 1993 former anal-ysis of fractional equations describing slow diffusionand interpolating diffusion and wave-propagation, thepresent author [28], introduced the so called auxiliaryfunctions of the Wright type F ν ( z ) := W − ν, ( − z ) and M ν ( z ) := W − ν, − ν ( − z ) with < ν < , in or-der to characterize the fundamental solutions for typ-ical boundary value problems, as it is shown in theprevious sections. Being then only aware of the Hand-book of the Bateman project, where the parameter λ of the Wright function W λ,µ ( z ) was erroneously re-stricted to non-negative values, the author thought tohave originally extended the analyticity property ofthe original Wright function by taking ν = − λ with ν ∈ (0 , . Then the function M ν was referred toas the Mainardi function in the 1999 treatise by Pod-lubny [50] and in some later papers and books.It was Professor B. Stankovi´c during the presentationof the paper by Mainardi & Tomirotti [42] at the Con-ference
Transform Methods and Special Functions,Sofia 1994 , who informed the author that this exten-sion for − < λ < had been already made byWright himself in 1940 (following his previous papersin the thirties), see [66]. In a 2005 paper published inFCAA [35], devoted to the 80th birthday of Profes-sor Stankovi´c, the author used the occasion to renewhis personal gratitude to Professor Stankovi´c for thisearlier information that led him to study the originalpapers by Wright and to work (also in collaboration)on the functions of the Wright type for further applica-tions, see e.g. [16, 17] and [38]. For the above reasonsthe author preferred to distinguish the Wright func-tions into the first kind ( λ ≥ ) and the second kind( − < λ < ). It should be noted that in the book by Pr¨uss (1993)[53] we find a figure quite similar to our Fig 15-top re-porting the M -Wright function in linear scale, namelythe Wright function of the second kind. It was de-rived from inverting the Fourier transform expressedin terms of the Mittag-Leffler function following theapproach by Fujita [12] for the fundamental solutionof the Cauchy problem for the diffusion-wave equa-tion, fractional in time. However, our plot must beconsidered independent by that of Pr¨uss because theauthor (Mainardi) used the Laplace transform in hisformer paper presented at WASCOM, Bologna, Oc-tober 1993 [28] (and published later in a number ofpapers and in his 2010 book) so he was aware of thebook by Pr¨uss only later. References: [1] J.H. Barret (1954). Differential equations of non-integer order,
Canad. J. Math , 529–541.[2] P.W. Buchen and F. Mainardi (1975). Asymptotic ex-pansions for transient viscoelastic waves, Journal deM´ecanique , 597–608.[3] M. Caputo and F. Mainardi (1971a). A new dissi-pation model based on memory mechanism, Pureand Appl. Geophys. (PAGEOPH) , 134–147.[Reprinted in Fract. Calc. Appl. Anal. No 3, 309–324 (2007)][4] M. Caputo and F. Mainardi (1971b). Linear modelsof dissipation in anelastic solids,
Riv. Nuovo Cimento (Ser. II) , 161–198.[5] K.S. Cole (1933). Electrical conductance of biolog-ical systems, Electrical excitation in nerves, in Pro-ceedings Symposium on Quantitative Biololgy , ColdSpring Harbor, New York, Vol. 1, pp. 107–116.[6] A. Consiglio & F. Mainardi (2019). On the evo-lution of fractional diffusive waves.
Ricerche diMatematica , published on line 06 Dec. 2019.DOI: 10.1007/s11587-019-00476-6 [E-print: arxiv.org/abs/1910.1259 ][7] H.T. Davis (1936).
The Theory of Linear Operators ,The Principia Press, Bloomington, Indiana.[8] G. Doetsch (1974).
Introduction to the Theory andApplication of the Laplace Transformation , SpringerVerlag, Berlin.[9] M.M. Dzherbashyan (1966).
Integral Transforms andRepresentations of Functions in the Complex Plane ,Nauka, Moscow [in Russian]. There is also thetransliteration as Djrbashyan.10] A. Erd´elyi, W. Magnus, F. Oberhettinger and F. Tri-comi, (1955).
Higher Transcendental Functions , 3-rd Volume, McGraw-Hill, New York . [BatemanProject].[11] W. Feller (1971).
An Introduction to Probability The-ory and its Applications . Wiley, New York, Vol. II,Second Edition.[12] Y. Fujita (1990). Integro-differential equation whichinterpolates the heat equation and the wave equation,I, II,
Osaka J. Math. , 309–321, 797–804.[13] L.G. Gatteschi (1973). Funzioni Speciali , UTET,Torino.[14] W. Gawronski (1984). On the bell-shape of stabledistributions,
Annals of Probability , 230–242.[15] R. Gorenflo, A.A. Kilbas, F. Mainardi & S. Ro-gosin (2014). Mittag-Leffler Functions. Related Top-ics and Applications , Springer, Berlin. Second Edi-tion in preparation.[16] R. Gorenflo, Yu. Luchko and F. Mainardi (1999).Analytical properties and applications of the Wrightfunction.
Fract. Calc. Appl. Anal. , 383–414.[17] R. Gorenflo, Yu. Luchko and F. Mainardi (2000).Wright functions as scale-invariant solutions of thediffusion-wave equation, Journal of Computationaland Applied Mathematics , 175–191.[18] R. Gorenflo and F. Mainardi (1996). Frac-tional oscillations and Mittag-Leffler func-tions, Preprint No A-96-14, FachbereichMathematik und Informatik, Freie Universit¨atBerlin, Serie Mathematik, pp. 22. [E-print: ][19] R. Gorenflo and F. Mainardi (1997). Fractional calcu-lus: integral and differential equations of fractionalorder, in: A. Carpinteri and F. Mainardi (Editors):
Fractals and Fractional Calculus in Continuum Me-chanics , Springer Verlag, Wien, pp. 223–276. [E-print: http://arxiv.org/abs/0805.3823 ][20] R. Gorenflo and F. Mainardi (1998). Randomwalk models for space-fractional diffusion processes,
Fract. Calc. Appl. Anal. No 2, 167–190.[21] B. Gross (1947). On creep and relaxation,
J. Appl.Phys. , 212–221.[22] E. Hille and J.D. Tamarkin (1930). On the theory oflinear integral equations, Ann. Math. , 479–528.[23] V. Kiryakova (1994). Generalized Fractional Calcu-lus and Applications , Longman & J. Wiley, Harlow -New York. [Pitman Research Notes in Mathematics,Vol. 301] [24] A. Liemert and A. Klenie (2015). Fundamental Solu-tion of the Tempered Fractional Diffusion Equation.
J. Math. Phys , 113504.[25] Yu. Luchko (2000). On the asymptotics of zeros ofthe Wright function. Zeitschrift f¨ur Analysis und ihreAnwendungen , , 597–622.[26] Yu. Luchko (2019). The Wright function and its ap-plications, in A. Kochubei, Yu.Luchko (Eds.), Hand-book of Fractional Calculus with Applications
Vol.1: Basic Theory, pp. 241- 268. De Gruyter GmbH,2019, Berlin/Boston. Series edited by J. A.TenreiroMachado.[27] Yu. Luchko & V. Kiryakova (2013). The Mellin in-tegral transform in fractional calculus,
Fract. Calc.Appl. Anal. , 405–430.[28] F. Mainardi (1994). On the initial value problem forthe fractional diffusion-wave equation, in: Rionero,S. and Ruggeri, T. (Editors), Waves and Stability inContinuous Media . World Scientific, Singapore, pp.246–251. [Proc. VII-th WASCOM, Int. Conf. ”Wavesand Stability in Continuous Media”, Bologna, Italy,4-7 October 1993][29] F. Mainardi (1996a). The fundamental solutionsfor the fractional diffusion-wave equation.
AppliedMathematics Letters , No 6, 23–28.[30] F. Mainardi (1996b). Fractional relaxation-oscillation and fractional diffusion-wave phenomena.
Chaos, Solitons & Fractals , 1461–1477.[31] F. Mainardi (1997). Fractional calculus: some ba-sic problems in continuum and statistical mechanics,in: A. Carpinteri and F. Mainardi (Editors): Frac-tals and Fractional Calculus in Continuum Mechan-ics , Springer Verlag, Wien, pp. 291–348. [E-print: http://arxiv.org/abs/1201.0863 ][32] F. Mainardi (2010).
Fractional Calculus and Wavesin Linear Viscoelasticity . Imperial College Press,London and World Scientific, Singapore.[33] F. Mainardi and A. Consiglio (2020). The Wrightfunctions of the second kind in Mathemat-ical Physics,
Mathematics No 6, 884/1–26 .[DOI: 10.3390/MATH8060884] [E-print arXiv:2007.02098 ][34] F. Mainardi and R. Gorenflo (2007). Time-fractionalderivatives in relaxation processes: a tutorial sur-vey,
Fract. Calc. Appl. Anal. , 269–308. [E-print: http://arxiv.org/abs/0801.4914 ][35] F. Mainardi, R. Gorenflo and A. Vivoli (2005).Renewal processes of Mittag-Leffler and Wrighttype, Fract. Calc. Appl. Anal. , 7–38. [E-print http://arxiv.org/abs/math/0701455 ]36] F. Mainardi, Yu. Luchko and G. Pagnini(2001). The fundamental solution of the space-time fractional diffusion equation, Fract.Calc. Appl. Anal. , 153–192. [E-print: arxiv.org/abs/cond-mat/0702419 ][37] F. Mainardi, A. Mura and G. Pagnini (2010), The M Wright function in time-fractional diffusion pro-cesses: a tutorial survey,
International Journalof Differential Equations
Vol. 2010, Article ID104505, 29 pages. [DOI:10.1155/2010/104505] [E-print http://arxiv.org/abs/1004.2950 ][38] F. Mainardi and G. Pagnini (2003). The Wrightfunctions as solutions of the time-fractional diffu-sion equation.
Applied Mathematics and Computa-tion
No 1, 51–62.[39] F. Mainardi, G. Pagnini and R. Gorenflo(2003). Mellin transform and subordinationlaws in fractional diffusion processes,
Fract.Calc. Appl. Anal. No 4, 441–459. [E-print: http://arxiv.org/abs/math/0702133 ][40] F. Mainardi, G. Pagnini and R.K. Saxena (2005). Fox H functions in fractional diffusion, J. Comp. Appl.Math. , 321–331.[41] F. Mainardi, P. Paradisi and R. Gorenflo (1998),Probability distributions generated by fractionaldiffusion equations, Invited Lecture, Work-shop on Econophysics, Budapest 21-27 July1997. LaTeX Pre-print, Dept. of Physics,Bologna, January 1998, pp. ii +39. [E-print http://arxiv.org/abs/0704.0320 ] N.B.It would have appeared in J. Kertesz and I. Kondor(Editors),
Econophysics: an Emerging Science,Kluwer, Dordrecht , book NOT published![42] F. Mainardi and M. Tomirotti (1995). On a specialfunction arising in the time fractional diffusion-waveequation, in: P. Rusev, I. Dimovski and V. Kiryakova,(Editors),
Transform Methods and Special Functions,Sofia 1994 , Science Culture Technology, Singapore,pp. 171–183.[43] F. Mainardi and M. Tomirotti (1997). Seismic PulsePropagation with Constant Q and Stable ProbabilityDistributions. Annali di Geofisica
40, 1311–1328. [E-print: arxiv.org/abs/1008.1341 ][44] J. Mikusi´nski (1959). On the function whose Laplacetransform is exp ( − s α ) , Studia Math. , 191–198.[45] G.M. Mittag-Leffler (1903a). Une g´en´eralisation del’int´egrale de Laplace-Abel, C.R. Acad. Sci. Paris (Ser. II) , 537–539.[46] G.M. Mittag-Leffler (1903b). Sur la nouvelle fonc-tion E α ( x ) , C.R. Acad. Sci. Paris (Ser. II) , 554–558. [47] G.M. Mittag-Leffler (1904). Sopra la funzione E α ( x ) , Rendiconti R. Accademia Lincei (Ser. V) ,3–5.[48] G.M. Mittag-Leffler (1905). Sur la repr´esentationanalytique d’une branche uniforme d’une fonctionmonog`ene, Acta Math. , 101–181.[49] R.B. Paris (2019). Asymptotics of the special func-tions of fractional calculus, in: A. Kochubei,Yu.Luchko (Editors), Handbook of Fractional Calcu-lus with Applications , Vol. 1: Basic Theory, pp. 297-325. De Gruyter GmbH, 2019 Berlin/Boston. Seriesedited by J. A.Tenreiro Machado.[50] I. Podlubny (1999).
Fractional Differential Equa-tions , Academic Press, San Diego.[51] H. Pollard (1946). The representation of exp ( − x λ ) as a Laplace integral, Bull. Amer. Math. Soc. , 908–910.[52] H. Pollard, (1948). The completely monotonic char-acter of the Mittag-Leffler function E α ( − x ) , Bull.Amer. Math. Soc. , 1115–1116[53] J. Pr¨uss, (1993). Evolutionary Integral Equations andApplications . Birkhauser Verlag, Basel.[54] A. Saa and R. Venegeroles (2011). Alternative nu-merical computation of one-sided L´evy and Mittag-Leffler distributions,
Phys. Rev. E
84, 026702.[55] G. Sansone and J. Gerretsen. (1960).
Lectures on theTheory of Functions of a Complex Variable , Vol. I.
Holomorphic Functions , Nordhoff, Groningen.[56] W.R. Schneider (1986). Stable distributions: Foxfunction representation and generalization, in S. Al-beverio, G. Casati and D. Merlini, D. (Editors),
Stochastic Processes in Classical and Quantum Sys-tems , Springer Verlag, Berlin, pp. 497–511. [LectureNotes in Physics, Vol. 262][57] T. Simon (2015). Positive stable densities and thebell-shape,
Proc. Amer. Math. Soc.
No 2, 885–895.[58] B. Stankovi ˇ c (1970). On the function of E.M. Wright. Publ. de l’Institut Math´ematique, Beograd, NouvelleS´er. , 113–124.[59] F.G. Tricomi (1959). Funzioni Speciali , Gheroni,Torino.[60] A. Wiman (1905a). ¨Uber den Fundamentalsatz derTheorie der Funkntionen E α ( x ) . Acta Math. ,191–201.[61] A. Wiman (1905b). ¨Uber die Nullstellen der Funkn-tionen E α ( x ) . Acta Math. , 217–234.62] R. Wong and Y.-Q. Zhao (1999a). Smoothing ofStokes’ discontinuity for the generalized Bessel func-tion. Proc. R. Soc. London A
Proc. R. Soc. London A
Journal LondonMath. Soc. , 8, 71–79.[65] E.M. Wright (1935). The asymptotic expansion ofthe generalized Bessel function.
Proc. London Math.Soc. (Ser. II)
38, 257–270.[66] E.M. Wright (1940). The generalized Bessel functionof order greater than one.