A Tutorial on the Basic Special Functions of Fractional Calculus
Francesco Mainardi
math.GM 26A33, 33E12, 34A08, 34C26

References:

12.5pt










A Tutorial  on the Basic  Special Functions of Fractional Calculus
     
[t]2.3in 12.5pt
     FRANCESCO MAINARDI
       University Bologna and INFN
    Department of Physics & Astronomy


    francesco.mainardi@bo.infn.it
 0in







 
 
   
[b]6.9in 12.5pt Abstract:


In this tutorial survey we recall  the basic properties of the 
special function of the Mittag-Leffler and Wright type that are known  to be relevant
in processes dealt with the fractional calculus.
We outline the major applications of these functions.
For the Mittag-Leffler functions 
we analyze  the Abel integral equation of the second kind and 
the fractional relaxation and oscillation phenomena.
For the Wright functions we distinguish them in two kinds.
We mainly stress the relevance of the Wright functions of the second kind  in probability theory with particular regard to the so-called 
M-Wright functions that generalizes the Gaussian  
and is related  with the time-fractional diffusion equation. 



 [4mm] Key–Words:  
Mittag-Leffler functions,
Wright functions, Fractional Calculus,
Laplace, Fourier  and Mellin transforms,  
 Probability theory,   Stable distributions.
 



    
   Paper published in
WSEAS TRANSACTIONS on MATHEMATICS, Vol 19 (2020), pp. 74–98.

 DOI: 10.37394/23206.2020.19.8
========================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================



empty empty



  
  
   
  

§ INTRODUCTION
 

The special functions of the Mittag-Leffer and Wright type in general play a very important role in the theory of the fractional  differential and integral  equations. 
The purpose of this tutorial survey  is to outline the relevant properties of
 the these functions outlining their applications.


   This work is organized as follows. 
 In Section 2, we recall  the essentials  of the fractional calculus
 that provide the necessary notions for the applications.


In section 3, we start to define the Mittag-Leffler functions.
For this purpose we  introduce
the Gamma function and the classical Mittag-Leffler functions of one and two parameters. Then we deal with the auxiliary functions of the Mittag-Leffler type
 to be used in the next sections.




 

 In Section 4 we apply the above auxiliary functions of the Mittag-Leffler type  to solve the Abel integral equations of the second kind, that are noteworthy cases of Volterra integral
  equations.
  
 
  In Section 5 we finally consider the most famous applications of the auxiliary functions of the Mittag-Leffler tyewpe, that is the solutions of the time fractional differential equations
  governing the phenomena of fractional relaxation and fractional oscillations
  
  
  In Section 6  we start to define the Wright  functions.
  For this purpose    we distinguish two kinds of these functions.
  Particular attention is devoted   to two special cases of the Wright function of the second kind introduced by Mainardi in the 1990's 
 in virtue of their importance in probability theory and for the time-fractional diffusion equations.
  Nowadays in the FC literature they are referred to as the Mainardi functions. 
In contrast to the general case of the Wright function, they depend just on one parameter ν∈[0,1).
One of the Mainardi functions, known as the M-Wright function,  generalizes the
 Gaussian function  and degenerates to the delta function in the limiting case ν=1.
 

Then, in Section 7  we recall how the Mainardi functions are related to
 an important class of the probability density functions (pdf's) known  as the extremal Lévy stable densities. This  emphasizes the relevance of the Mainardi functions 
 in the probability theory   independently on the framework of the fractional diffusion equations.
  We present some plots of the symmetric M-Wright function  on
  for several parameter values ν∈ [0, 1/2] and ν∈ [1/2,1]. 
 
 
 Finally concluding remarks are carried out in Section 8 and two 
 tutorial appendices on stable distributions and the time fractional diffusion equation are added for readers' convenience. 
 

 The paper is competed with  historical and bibliographical concerning the past approach of  the author towards the Wright functions. 

 


§ THE ESSENTIALS OF  FRACTIONAL 
 CALCULUS
 

This section is mainly  based on the 1997 CISM survey by  Gorenflo and Mainardi <cit.>.
 The Riemann-Liouville fractional integral of order μ >0  is defined as

    _t J^μ  f(t) :=    Γ(μ ) ∫_0^t  (t-τ)^μ-1  f(τ )  dτ , 
    (2.1)

where 

    Γ(μ):= ∫_0^∞^-uu^μ-1  du ,  Γ(n+1)= n!

is theGamma function.


By convention _tJ^0 = I (Identity operator). 
We can prove semi-group property

    _tJ^μ  _tJ^ν =  
       _tJ^ν  _tJ^μ =   _tJ^μ +ν ,  μ ,ν≥ 0 .  
     (2.2)

Furthermore we have for t>0




    _t J^μ  t^γ=
       Γ(γ +1)Γ(γ +1+ μ  ) 
         t^γ+μ ,
      μ≥ 0 ,
       γ >-1 .
     (2.3)


The fractional derivative of order μ >0
in the Riemann-Liouville sense  is defined  
 as the operator _tD^μ 
 




    _tD^μ  _tJ^μ  = I ,  μ >0 . (2.4)

If we take 
  m-1 <μ≤ m ,  with m∈ we recognize from Eqs. (2.2) and (2.4)

    _tD^μ f(t) :=    _tD^m  _tJ^m-μ f(t) , (2.5)

hence, for  m-1 <μ  < m,

    _tD^μ  f(t) =
      d^m dt^m[
      Γ(m-μ )∫_0^t
        f(τ) dτ (t-τ )^μ  +1-m], 
    (2.5a)
 
 and, for μ=m,
 
    _tD^μ  f(t) =
    d^m  dt^m f(t) ,.
    (2.5b)

For completion 
_tD^0 = I.
The semi-group property is no longer valid but  
for t>0  

    _t D^μ  t^γ=
       Γ(γ +1)Γ(γ +1-μ ) 
         t^γ-μ,
      μ≥ 0,
       γ >-1.(2.6)

However, the property _t D^μ =  _t J^-μ  is not generally valid!

The fractional derivative of order μ∈ (m-1,m] (m∈) in the
Caputo sense  is defined as the operator
_tD_*^μ  such that,

    _tD_*^μ f(t) :=    _tJ^m-μ  _tD^m  f(t) ,
    (2.7)

hence, for f  m-1 <μ  < m,

    _tD_*^μ f(t) =  
     Γ(m-μ ) ∫_0^t
     f^(m)(τ)  dτ (t-τ )^μ  +1-m ,
    (2.7a)

and for μ=m

    _tD_*^m f(t)=:
    d^m  dt^m f(t) .
    (2.7b)

Thus, when the order is not integer the two fractional derivatives
 differ in that  the derivative of order m
does not generally commute with the fractional integral.
We point out that the  Caputo fractional derivative
  satisfies the  relevant property
of being zero when applied to a constant, and, in general,
to any power function  of non-negative integer degree less than m ,
if its order μ is such that m-1<μ≤ m .

Gorenflo and Mainardi (1997) <cit.>
have shown the essential  relationships between the two fractional
derivatives 
(when both of them exist),
for  m-1 <μ <m,

    _tD_*^μ f(t)   =
    {[         _tD^μ [ f(t) -
      ∑_k=0^m-1 f^(k)(0^+) t^k k!] ,; _tD^μ  f(t) -
        ∑_k=0^m-1f^(k)(0^+)  
    t^k-μΓ(k-μ+1) . ].
    (2.8)

In particular, if m=1 so that  0<μ <1,
we have

    _tD_*^μ f(t) =  {[     _tD^μ [ f(t) -
       f(0^+) ]  ,; _tD^μ f(t) - f(0^+)  t^-μΓ(1-μ) . ].
    (2.9)

The Caputo fractional derivative,
represents a sort of 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→ 0^+ f^(k)(t)
are required to be finite for k=0,1, 2. … m-1.






We observe the different behaviour  of the two fractional derivatives
at the end points of the interval (m-1,m)
namely when the order is any positive integer:
whereas  _tD^μ is, with respect to its order μ ,
 an operator continuous
 at any positive integer,
_tD_*^μ   is an operator left-continuous
since

    [ lim_μ→ (m-1)^+ _tD_*^μ f(t) 
     = 
       f^(m-1)(t) - f^(m-1) (0^+) ,;                        lim_μ→ m^- 
     _tD_*^μ f(t)  =   f^(m)(t) . ](2.10)






We also note for m-1 < μ≤ m ,

    _t D^μ  f(t)  =  _tD^μ  g(t)
        ⟺ 
      f(t) = g(t) + ∑_j=1^m c_j  t^μ -j ,  (2.11)



    _tD_* ^μ  f(t)  =   _tD_*^μ  g(t)
        ⟺ 
      f(t) = g(t) +  ∑_j=1^m c_j  t^m-j . (2.12)

In these formulae the coefficients c_j are arbitrary constants.



We point out the major utility
of the Caputo fractional derivative
in treating initial-value problems for physical and engineering
applications where initial conditions are usually expressed in terms of
integer-order derivatives. This can be easily seen
using the Laplace transformation.
Writing the Laplace transform of a sufficiently well-behaved function
f(t) (t≥ 0) as

    ℒ{ f(t);s}=  f(s)
     := ∫_0^∞^  -st  f(t)  dt ,

 the known rule for
 the ordinary derivative of integer order m ∈ is 

    ℒ{ _tD^m  f(t) ;s} =
          s^m  f(s)
       -∑_k=0^m-1    s^m  -1-k  f^(k)(0^+)  ,

where 

    f^(k)(0^+) := lim_t→ 0^+ _tD^k f(t) .
 

For the Caputo derivative of order μ∈ (m-1,m] (m ∈) 
we have

    [                                                                                ℒ{ _tD_*^μ f(t) ;s}  =
          s^μ f(s)
       -∑_k=0^m-1    s^μ  -1-k  f^(k)(0^+) ,;                                                                                                                         f^(k)(0^+) := lim_t→ 0^+  _tD^kf(t) . ](2.13)


The corresponding rule for the 
Riemann-Liouvile  derivative of order μ   is


    [                                                                           ℒ{ _tD_t^μ  f(t);s}  =
         s^μ f(s)
       -∑_k=0^m-1 s^m -1-k  g^(k)(0^+) ,;                                                                                      g^(k)(0^+) := lim_t→ 0^+  _tD^kg(t) ,  g(t):=  _tJ^m-μ f(t)  . ](2.14)

 Thus the rule  (2.14) is     more cumbersome to be used
than (2.13) since it requires initial values  concerning
an extra function g(t) related to the given f(t)  
through a  fractional integral.






However, when all the limiting   values f^(k)(0^+)
are finite and the order is not integer, we can prove
by  that all g^(k)(0^+)
 vanish  so that  
the formula (2.14)  simplifies into

    ℒ{ _tD^μ  f(t);s} =
          s^μ f(s)  ,
    	     m-1 <μ< m .
    	  (2.15)

	  For this proof it is sufficient to
	  apply the Laplace transform to Eq. (2.8),
	  by recalling that  
	  
    ℒ{ t^ν ;s} =
         Γ(ν+1)/ s^ν+1 , ν >-1 ,(2.16)

	 and then to compare  (2.13) with (2.14). 




§ IHE FUNCTION OF MITTAG-LEFFLER TYPE

We note that the   Mittag-Leffler functions are present
in the Mathematics Subject Classification since the year 2000
under the number 33E12 under recommendation of Prof. Gorenflo.
A description of the most important properties of these
functions (with relevant references up to the fifties)
can be found in the third volume of the  Bateman Project
edited by Erdelyi et al. (1955)
in the chapter XVIII on Miscellaneous Functions
<cit.>.
The  treatises	where great attention is
devoted to the functions of the Mittag-Leffler type is that
 by Djrbashian (1966)
<cit.>, unfortunately in Russian.
We also recommend
the classical treatise on complex functions by
Sansone & Gerretsen (1960) <cit.>.







Nowadays the Mittag-Leffler functions are widely used in the framework
of integral and differential equations of fractional order, as shown  
in all treatises on fractional calculus.


In view of its several applications in Fractional Calculus the Mittag-Leffler function 
was referred to as the Queen Function of Fractional Calculus by Mainardi & Gorenflo (2007) <cit.>.
Finally, the functions of the Mittag-Leffler type have  found an exhaustive treatment in the treatise by
 Gorenflo, Kilbas, Mainardi & Rogosin (2014)
 <cit.>













As pioneering works of mathematical nature in the field of fractional
integral and differential equations, we like
to quote 

Hille & Tamarkin (1930)
<cit.>
who  have provided the solution
of the Abel integral equation of the second kind
in terms of a Mittag-Leffler  function,
and  Barret (1954) <cit.>
 who has expressed the   general
solution of the linear	fractional differential equation with constant
coefficients in terms of Mittag-Leffler functions.

As former applications in physics we like to quote the contributions
by  Cole (1933) <cit.> 
in connection with nerve conduction, see
  also  Davis (1936) <cit.>,
and by	Gross (1947) <cit.>
in connection with  mechanical relaxation.
Subsequently, Caputo & Mainardi (1971a), (1971b)
 <cit.>  
have shown that  Mittag-Leffler functions are present whenever
derivatives of fractional order are
introduced in the constitutive equations of a linear viscoelastic body.
Since then, several other authors
have pointed out  the relevance of these functions for
fractional viscoelastic models, see Mainardi's  survey (1997)
<cit.>
and his (2010) book
<cit.>.


















 §.§ 
The Gamma function: Γ(z)



Re  (z)   Im  (z)
arg  z
Re  (s)
Re  (p)     Re  (q)
Re   Im 
B
Γ(z)
D_Γ
Ai 
Erfc 
ŭ
u 
u 



Γ(z)    Γ(α)
Γ(α ,  z)
γ(α ,  z)
D_Γ
e E
Ai 
Erfc 
Ei 
Ein 
loglog 

The   Gamma function  Γ(z) is the most widely
used of all the special

functions: it is usually discussed first because it appears in almost every
integral or series representation of  other advanced mathematical
functions.
The first occurrence of the gamma function happens in 1729 
in a correspondence between Euler and Goldbach.
We  take as its  definition the integral formula due to Legendre (1809)

    Γ(z) :=
      ∫_0^∞ u^ z-1 ^ -u  du  ,  >0 .
    (3.1)

This integral representation is the most common for , even if
it is valid only in the right half-plane of .
The  analytic continuation   to the left half-plane is possible in different ways.
As  will be shown hereafter ,  the domain of analyticity D_Γ
of     turns out to be

    D_Γ =  - {0, -1, -2, …} .

        
        
        The most common continuation is carried out by the mixed representation due to Mittag-Leffler: 
   and reads for  z ∈ 
   
    Γ (z) =  ∑_n=0^∞(-1)^n/n!(z+n)  +  ∫_1^∞ e^-u u^z-1  du  .(3.2)

This representation can be obtained  from the so-called
Prym's decomposition, namely by
splitting the
integral in (3.1) into 2 integrals, the one over the interval
0 ≤ u ≤ 1 which is then developed in a series,
the other over the interval
1 ≤ u ≤∞, which, being uniformly convergent inside ,
 provides an entire function.
The terms of the series (uniformly convergent inside )
provide the principal parts of 
at the corresponding  poles z_n= -n .
So we recognize
that  is analytic
in the entire complex plane except at the points z_n=-n (n=0,1, …),
  which turn out to be  simple poles with residues  R_n = (-1)^n/n!.
The point at  infinity, being an accumulation point of poles,
is an essential non-isolated singularity.
Thus    is a transcendental meromorphic function.     


The reciprocal of the Gamma function turns out to be an entire function. its integral representation in the complex plane was due to Hankel (1864) and reads

    Γ (z) =
     2π i ∫_Ha^u/u^z  du  ,   z ∈ ,
  
where Ha denotes the Hankel path
 defined as a contour that
begins at u =-∞ - ia (a>0), encircles the  branch
cut that lies along the negative real axis, and ends up at
u= - ∞ + ib (b>0). 
Of course, the branch 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 around the origin, described in the counterclockwise
direction.







 §.§ 
The classical Mittag-Leffler functions


The Mittag-Leffler functions, that we denote by E_α(z), E_α, β(z) are so named  
in honour of Gösta Mittag-Leffler, the eminent Swedish mathematician, who introduced and 
investigated these functions
 in a series of notes starting from 1903 in the framework of the theory of entire 
 functions 
 <cit.>.  
The functions are defined by the series representations, convergent in the whole complex plane  for Re (α) >0}

    {[         E_α (z) := ∑_n=0^∞z^n/Γ (α n+1),; E_α, β (z ) := ∑_n=0^∞z^n/Γ(α n+ β   ) , ].
       (3.3)

   with β∈..
   Originally Mittag-Leffler assumed only the parameter α and assumed it as positive,
    but soon later the generalization 
   with two complex parameters was considered
   by Wiman. <cit.>.
   In both cases the Mittag-Leffler functions are entire of order 
   1/Re(α). 
   
 
The integral representation
 for z∈
 introduced by Mittag-Leffler can be written as
   
    E_α (z) =
      2π i ∫_Haζ ^α -1 ^ ζ/ζ ^α -z 
            dζ ,   α  >0.
     (3.4)

 
 Using  series representations of the Mittag-Leffler functions  it is easy to recognize

    -0.2truecm
    {-0.2truecm
    [     E_1,1(z)= E_1(z)  = ^z,   E_1,2(z)= ^z - 1/z,;                                                  ;       E_2,1(z^2)= cosh(z),   E_2,1(-z^2)= cos(z),;                                                  ; E_2,2(z^2)= sinh(z)/z,
         E_2,2(-z^2)=sin(z)/z, ] . (3.5)

and more generally

    {[     -0.2truecm
    E_α,β(z) + E_α,β(-z)= 2  E_2α,β(z^2),;                                                     ; -0.2truecm
    E_α,β(z) - E_α,β(-z)= 2z  E_2α,α+ β(z^2). ] . (3.6)

Furthermore, for α=1/2,

    [ E_1/2 (±  z^1/2) 
    -0.2truecm =
         ^ z [ 1+  (± z^1/2) ];                   -0.2truecm          = ^ z   (∓ z^1/2) , ](3.7)

 where erf (erfc) denotes the (complementary) error function
defined for z∈ as

    (z) := 2√(π) ∫_0^z ^ -u^2 du  ,
       (z) :=	1 -   (z) .


































 

 §.§ The auxiliary functions of the Mittag-Leffler type

 
 
 
 
In view of applications we introduce the following causal functions in time domain

    e_α (t;λ ) := E_α(-λ  t^α)
       ÷s^α -1/ s^α +λ ,
      (3.8)

  
 
    e_α,β  (t;λ):=
       t^β   -1   E_α,β(-λ  t^α)
       ÷s^α -β/ s^α +λ , (3.9)

	
    [         e_α,α  (t;λ) 
     -0.2truecm   :=
        t^α  -1   E_α,α(-λ  t^α);                        -0.2truecm = d/dt e_α(-λ  t^α)÷
     -λ/ s^α +λ. ](3.10)
	   






 A  function f(t) defined in ^+ is 
completely monotone (CM)  if (-1)^n  f^n(t) ≥ 0.
The function  ^-t is the prototype of a CM function.
For a Bernstein theorem a generic CM function reads



    f(t)= ∫_0^∞^-rt  K(r) r  ,   K(r) ≥ 0 .(3.11)


We have  for λ >0 

    [ e_α,β(t;λ):=
      t^β-1  E_α,β( - λ t^α);                                      ;                  CM iff   0<α≤β≤  1 . ](3.12)

  
  
  
Using the Laplace transform we can prove, following Gorenflo and Mainardi (1997) <cit.>
    that  for 0<α <1  (with λ=1)

    -0.1truecm
     E_α(-  t^α) ≃{-0.2truecm
    [ 1-  t^α/Γ(α +1)
     ⋯-0.2truecm                      t→ 0^+,;     t^-α/Γ(1-α)
     ⋯-0.2truecm                      t → +∞, ]. (3.13)

and 

    E_α(-  t^α)=
    ∫_0^∞^-rt  K_α(r)  r (3.14)

 with
 
    K_α (r)
    
    
     =π  r^α -1 sin (απ)/r^2α + 2  r^α cos  (απ) +1
    	 > 0 . (3.15)

   
In the following sections we will outline the key role of the auxiliary functions
in the treatment of integral and differential equations of fractional order, including the Abel 
	integral equation of the second kind and the differential equations for
	fractional relaxation and oscillation.   







Before closing this section we find it convenient to provide the plots of the functions

    ψ_α(t) = e_α(t) := E_α (-t^α) , (3.16)

and

    ϕ _α (t) =  t^-(1-α  )  E_α  ,α(-  t^α):=  - d/dt E_α (- t^α) , (3.17)

for t≥ 0 and for some rational values of α∈ (0,1].

For the sake of visibility, for both functions we have adopted linear and logarithmic scales.
 Logarithmic scales have been adopted in order to better point out their
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.18)











§ ABEL INTEGRAL EQUATION 
 OF THE SECOND KIND


Let us now consider the Abel equation of the  second kind
with  α >0 , λ∈:

    u(t)  + λΓ(α ) ∫_0^t u(τ) (t-τ )^1-α  dτ = f(t) ,
     	(4.1)

In terms of the fractional integral operator such equation reads

    ( 1 +λ  J^α)  u(t)=f(t) ,(4.2)

and consequently can be  formally solved as follows:

    [                                    u(t)                  = (1 + λ J^α)^-1  f(t);                                         = (1+ ∑_n=1^∞ (-λ) ^n   J^α n)   f(t) . ](4.3)




Recalling the definition of the fractional integral the formal solution reads

    u(t) = f(t) +  ( ∑_n=1^∞ (-λ) ^n  t_+^α n-1Γ(α n)) *  f(t) .
      (4.4)

  
  
  
  Recalling the definition of the 
function,
  
    e_α (t;λ ) := E_α (-λ  t^α ) =
    ∑_n=0^∞(-λ t^α)^n
     Γ (α n+1) ,
         (4.5)

where
E_α
denotes the  Mittag-Leffler function of order α ,
 we note that for t>0:

    ∑_n=1^∞ (-λ) ^n  t_+^α n-1Γ (α n)
     = d  dt  E_α  (-λ t^α )
     = e_α ^'(t;λ )   , (4.6)


Finally, the solution  reads

    u(t)= f(t)+ e_α^' (t;λ) *  f(t)  . (4.7)

 
Of course the above formal proof


can be made  rigorous.
Simply observe that because of the rapid growth of the gamma function the
infinite series in (4.4) and (4.6) are uniformly convergent in every
bounded interval of the variable t so that term-wise integrations and
differentiations are allowed.



However,
we prefer to use the alternative technique of Laplace transforms,
which will allow us  to obtain
the solution in different forms, including  the result (4.7).
Applying the Laplace transform to (4.1) we
obtain

    [1 +λ s^α] u̅(s) =  f̅(s)
      ⟹ u̅(s) = s^α s^α +λ f̅(s)
     . (4..8)

Now,  let us proceed to obtain the inverse Laplace transform
of (4.8) using    the following
Laplace transform pair	


   
    e_α (t;λ) := E_α (-λ  t^α)
      ÷ s^α -1 s^α +λ .
     (4.9)


We can choose  three different
ways to get the
inverse Laplace transforms   from (4.8), according to the standard rules.

Writing  (4.8) as

    u̅(s) = s  [
     s^α -1 s^α+ λ f̅(s)] ,
    (4.10a)
 we obtain

    u(t) =d dt ∫_0^t  f(t-τ )   e_α (τ;λ)
        dτ .      (4.11a)

If we write  (4.8) as

    u̅(s) = s^α -1 s^α +λ 
      [ s f̅(s) - f(0^+)] + f(0^+) s^α -1 s^α +λ ,	  (4.10b)

we obtain

    u(t) =  ∫_0^t
     f^'(t-τ)       e_α (τ;λ) 	dτ
       +	f(0^+)   e_α(t;λ ) .
       (4.11b)



We also note that,  e_α (t;λ) being a
function differentiable with respect to t
with 

    e_α (0^+;λ)= E_α	(0^+)=1,

there exists another possibility to re-write (4.8), namely

    u̅(s) = [ s  s^α -1 s^α+ λ -1]  f̅(s) +
      f̅(s) .
    (4.10c)

Then  we obtain

    u(t) =
       ∫_0^t   f(t-τ)  e_α^' (τ;λ)  dτ
      + f(t)  ,	  (4.11c)

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
 ℒ-transformable derivative.


 


§ FRACTIONAL RELAXATION
 AND  OSCILLATIONS

Generally speaking, we consider
the following differential equation of fractional order
α >0 , for t≥ 0:

    [                        D_*^α  u(t)
    -0.2truecm =
    D^α ( u(t) - ∑_k=0^m-1t^k  k! 
    u^(k) (0^+));                                    -0.2truecm                            = - u(t) + q(t)  , ](5.1)

where u=u(t) is the field variable and q(t) is a given function,
continuous for t≥ 0 .
Here m is a positive integer uniquely defined
by m-1 <α≤ m ,
which provides the number of the prescribed initial values
u^(k)(0^+) = c_k ,   k=0,1,2, … , m-1 .




In particular, we consider in detail the cases


(a) fractional relaxation  0<α≤ 1 ,


(b) fractional oscillation  1<α≤ 2 .	



The application of the Laplace transform

yields

    u(s) = ∑_k=0^m-1 c_k  s^α -k -1
    s^α +1	    + s^α+1  q(s) .
     (5.2)

Then,  putting  for	k=0,1,…,m-1 ,

    [ u_k(t) :=  J^k e_α (t)  ÷s^α -k -1 s^α +1 ,;       e_α (t)  :=E_α(-t^α)  ÷s^α-1 s^α +1 , ](5.3)

and using u_0(0^+) =1 ,
we find, 

    u(t) = ∑_k=0^m-1 c_k   u_k (t)
       - ∫_0^t  q(t-τ )  u_0 ^'(τ)  dτ .
        (5.4)

In particular, the  formula (5.4)
encompasses the solutions
for α = 1 , 2 ,  since

    α =1  ,   u_0(t) =e_1(t) = exp (-t) ,


    α=2, 
       u_0(t) =e_2(t)  =cos t,  
     u_1(t) = J^1 e_2(t) =sin  t.



When α is not integer, namely for m-1 <α <m ,
we note that m-1 represents
the integer part of α (usually denoted by [α])
and m  the number  of initial  conditions
necessary and sufficient to ensure the uniqueness of the
solution u(t). Thus the m functions
u_k(t) = J^k e_α (t) with k =0, 1, …, m-1
represent those particular solutions

of the homogeneous equation which satisfy
the initial conditions
u_k^(h) (0^+)  = δ _k h ,
h,k =0, 1, …, m-1   ,
and therefore they represent the fundamental solutions
of the fractional equation (5.1), in analogy with the case
α =m .  
Furthermore, the function
u_δ (t) = -u_0^'(t) =  -e_α ^' (t)
represents the impulse-response solution.
Now we derive the relevant properties of the basic
functions e_α (t)  directly   from their
Laplace representation for 0<α≤ 2,

    e_α (t) = 2π i ∫_Br^ st  s^α -1 s^α +1  ds
       ,(5.5)

 without detouring on the general theory of Mittag-Leffler functions
in the complex plane.
Here Br denotes a Bromwich path, a line
Re (s) = σ >0 
and Im (s) running from -∞ to +∞.



For transparency reasons, we  separately discuss the cases
(a) 0<α <1  and
    (b) 1<α <2 ,
recalling that in the limiting cases
α =1 ,  2 ,
we  know e_α (t) as elementary function,
 namely
   e_1(t) = ^ -t
and e_2(t) = cos  t .
For α not integer the power function  s^α is uniquely
defined as
 s^α = |s|^α ^  i arg s ,
with -π < arg  s <π ,
that is in the complex
s-plane cut along the negative real axis.
The essential step  consists in decomposing e_α (t) into two parts
according to
e_α(t)  = f_α(t) + g_α(t)  ,
as indicated below.
In case (a) the function f_α(t) ,
in case (b) the function -f_α(t)
is completely monotone; in both cases  f_α (t)
tends to zero as t tends to infinity,
from above in case (a), from below in case (b).
The other part, g_α (t) ,
is identically vanishing in case (a), but of
oscillatory character with exponentially	decreasing amplitude
in case (b).


For the oscillatory part we obtain
via the residue theorem of complex analysis,
when  1<α<2: 

    g_α (t) = 2α ^  t cos (π /α ) cos [t sin (πα) ].  (5.6)

We note that this function exhibits oscillations with
circular frequency
 
    ω  (α ) = sin (π/α)

and with an  exponentially decaying amplitude with  rate

    λ (α ) = |cos (π /α)| = - cos (π /α)
     .



For the  monotonic part we obtain

    f_α (t) := ∫_0^∞^  -rt 
         K_α (r)  dr , (5.7)

with

    [                                               K_α (r)                =
      - π Im ( .s^α -1 s^α +1|_ s=r ^iπ);                                                       = π  r^α -1 sin (απ)
        r^2α + 2  r^α cos  (απ) +1 . ](5.8)

This function K_α(r) vanishes identically if α is
an integer, it is positive for all r if 0<α <1 ,
negative for all r if 1<α <2 .
In fact in (5.8) the denominator  is, for α not integer, always
positive being	> (r^α -1)^2 ≥ 0 .


Hence f_α (t) has the aforementioned
monotonicity properties, decreasing towards zero in case (a),
increasing towards zero in case (b).
We  note that, in order to satisfy the initial condition
e_α (0^+) =1,  we find

    ∫_0^∞ K_α (r)  dr =1  if  0<α≤ 1 ,


    ∫_0^∞ K_α (r)  dr = 1- 2/α if  1<α≤ 2 .

In Figs. 3 and 4 we display the plots of K_α (r),
that we denote as the
basic spectral function,
 for some values of α in the intervals
(a) 0<α <1 ,
(b) 1<α <2  .





















 



In addition to the   basic fundamental solutions, u_0(t)= e_α (t) ,
we need to compute the impulse-response solutions
u_δ (t) = - D^1   e_α (t) for cases (a) and (b)
and,  only in case (b), the second fundamental solution
u_1(t) = J^1   e_α (t) .

For this purpose we note that in general  it turns out that

    J^k   f_α (t) =
        ∫ _0^∞^  -r t  K_α^k(r)  dr  ,
      (5.9)

with
  
    [                                                             K_α^k(r)                                             := (-1)^k   r^-k  K_α(r);                                                                      =
         	(-1)^k π  r^α -1-k sin (απ)
        r^2α + 2  r^α cos  (απ) +1 , ](5.10)

where K_α (r) = K_α^0(r) ,  and

    J^k g_α (t) =
     2α ^  t cos (π /α ) cos ł[t sin ł(πα)̊  - k πα]̊.
        (5.11)




 For the impulse-response solution we note that the effect of the differential operator D^1
 is the same as that of the virtual operator J^-1 .


Hence  the solutions  for the
fractional relaxation  are:


  (a) 0<α <1 ,

    u(t) = c_0   u_0(t)
    + ∫_0^t q(t-τ)  u_δ(τ)  dτ ,(5.12a)

where

    [     u_0(t) = ∫_0^∞^  -rt K_α^0(r) dr ,; u_δ (t) =
     -∫_0^∞^  -rt K_α^-1(r) dr , ](5.13a)

  with

    u_0(0^+) =1 ,   u_δ(0^+) =∞ ,


and for t →∞

    u_0(t) ∼t^-αΓ(1-α ) ,  
     u_1(t) ∼t^1-αΓ(2-α ) .(5.14a)



Hence  the solutions  for the
fractional oscillation  are:

 (b)  1<α <2 ,

    u(t) = c_0   u_0(t) + c_1   u_1(t)
        + ∫_0^t q(t-τ)  u_δ(τ)  dτ ,
    (5.12b)



    [                                         u_0(t)                     = ∫_0^∞^  -rt K_α^0(r)  dr;                                                   +
       2α ^ t cos (π /α ) cos[t sin ł(πα)]	 ,;                                         u_1(t)                       = ∫_0^∞^ -rt K_α^1(r) dr;                                                +
        2α ^  t cos (π /α ) cos[t sin (πα) -πα],;                                         u_δ(t)                   =  -∫_0^∞^  -rt K_α^-1(r) dr;                                                  -  2α ^   t cos (π /α ) cos[t sin (πα) +πα], ](5.13b)

with

    u_0(0^+) =1,   u_0^'(0^+)=0,


    u_1(0^+) =0,   u_1^'(0^+)=1,

 
    u_δ(0^+) =0,   u_δ^'(0^+) = +∞,

  and for t →∞
  
    {[                                  u_0(t) ∼t^-αΓ(1-α ),;                                 u_1(t) ∼t^1-αΓ(2-α ),;                            u_δ (t) ∼ - t^-α-1 Γ(-α ); ] .
     (5.14b)




In Figs. 2a and 2b we display the plots of the basic fundamental solution  for the
following cases, respectively :

 (a) α = 0.25 ,   0.50 ,  0.75 ,   1 ,


(b) α = 1.25 ,   1.50 ,  1.75 ,   2 ,


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)
for  fractional relaxation
and fractional oscillation,
all the fundamental and impulse-response solutions
exhibit an algebraic decay   as t →∞ ,
as discussed above.
This algebraic decay is the most important effect of
the non-integer derivative in our equations, which
dramatically differs from the
exponential decay present in the ordinary relaxation
 and damped-oscillation phenomena.
 




















 


We would like to remark  the difference between fractional relaxation
governed by the Mittag-Leffler type function
E_α (-a t^α) and stretched relaxation governed by a
stretched    exponential function exp (-b t^α )
with α  ,  a  , b >0 for t ≥ 0 .
A common behaviour is achieved only in a restricted range
0 ≤ t ≪  1 where

    [                         E_α (-a t^α) ≃ 1 - a Γ(α +1)  t^α
     =  1 - b   t^α;                                             ≃^-b  t^α,
        b = a Γ(α +1). ]




-0.15truecm
























0.5truecm
In Figs. 3a, 3b, 3c  for
α =0.25, 0.50, 0.75 
we have compared 
 E_α (-t^α) (full line)
with its asymptotic approximations
exp [- t^α/Γ(1+α)]
(dashed line)
valid for short times,  and t^-α /Γ(1-α)
(dotted line)
valid for long times.
We have adopted log-log  plots
  
in order to better achieve such a comparison

and 
the transition  from the stretched exponential
to the inverse power-law decay.
In Figs. 4a, 4b, 4c we have shown
  some plots of the basic fundamental
solution u_0(t) = e_α (t)
for   α = 1.25 ,   1.50 ,  1.75, respectively.
Here the algebraic decay of the fractional oscillation
can be recognized and compared with the two contributions provided by
f_α (monotonic behaviour, dotted line)
and g_α (t) (exponentially
damped oscillation, dashed line)



The zeros of the solutions of the fractional 
 oscillation
Now we find it interesting to carry out
some investigations
about the zeros of the basic fundamental solution
u_0(t) = e_α (t)   in the case (b) of fractional oscillations.
For the second fundamental solution
and the  impulse-response solution the analysis of the zeros
can be easily carried out analogously.


Recalling the first equation in (5.13b),
the required zeros of e_α (t) are the solutions of
the  equation
  
    e_α (t) =  f_α (t)  +
    2α ^  t cos  (π/α ) cos [t sin ł(πα)̊]
    =0 . (5.15)

We first note  that
the function e_α (t)  exhibits an odd number
of zeros, in that  e_α (0)= 1 ,
and,
for sufficiently large t,   e_α (t)
turns out to be permanently negative, as shown in (5.14b)
by the sign of Γ(1-α ) .


The smallest zero lies in the first positivity
interval of cos  [t sin (π/α)] ,
hence in the interval 0<t< π/[2  sin (π/α)] ;
all other zeros can only lie in the succeeding
positivity intervals   of cos  [t sin (π/α)] ,
in each of these two zeros are present as long as

    2α ^  t  cos  (π/α )≥ |f_α (t)| . (5.16)





When t is sufficiently large the zeros are expected
to be found approximately from the equation

    2α ^  t cos  (π/α )≈t^-α |Γ(1-α)| ,(5.17)

obtained from (5.15)
by ignoring the oscillation factor of  g_α (t)
and taking the first  term in the asymptotic  expansion
of  f_α (t).
As shown in the report <cit.>
 such
approximate equation turns out to be
useful when α→ 1^+ and α→ 2^- .



For α→ 1^+ ,  only one zero is present,
which is expected to be very far from the origin in view of the
large 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 only one zero exists  for α sufficiently close to 1.
Putting α = 1+ϵ
the asymptotic position T_* of this zero
can be found from the relation (5.17) in the limit
ϵ→ 0^+ .
Assuming in this  limit 
the first-order approximation,	we get

    T_*∼logł(2ϵ)̊ , (5.18)

which shows that T_* tends to infinity slower than
1/ϵ ,  as ϵ→ 0 .
For details see again the  1995 report
by Gorenflo & Mainardi <cit.>.



For α→ 2^-, there is an increasing number of zeros
up to infinity
since e_2(t) = cos  t has infinitely many zeros
[in t^*_n = (n+1/2)π ,  n = 0,1,…].
Putting now α =2-δ the asymptotic position T_* for
the largest zero
can be found again   from (5.17)
in the limit δ→ 0^+ .
Assuming in this  limit 
the first-order approximation,	we get

    T_* ∼12π δ log( 1 δ) .(5.19)


Now, for δ→ 0^+ the  length of the positivity
intervals of g_α (t) tends to π
and, as long as t ≤ T_* , there are two zeros in
each positivity interval. Hence, in the limit δ→ 0^+ ,
there is in average one zero per interval of  length π ,
so we expect that
N_*∼   T_*/π .



Remark : For the above considerations we got inspiration
from an interesting paper by Wiman (1905) <cit.>, who

at the beginning of the XX-th  century,
 after having treated
the Mittag-Leffler function  in the complex plane,
considered the	position of the zeros
of the function  on the negative real axis
(without providing  any detail).
The expressions  of  T_*  are in disagreement with those
by Wiman for  numerical factors; however,
the results of our numerical studies carried out in the 1995 report 
<cit.>
confirm and  illustrate the validity of the present   analysis.
Here, 
we find it interesting to analyse
the phenomenon of the transition
of the	(odd) number of  zeros as
1.4 ≤α≤ 1.8 .
For this purpose, in Table I we report	the intervals
of amplitude Δα  = 0.01  where
these transitions occur, and the
location T_* (evaluated  within a relative
error of 0.1%)  of the  largest zeros found at the
two extreme  values  of the above intervals.



 We recognize that the transition from 1 to 3 zeros occurs as
1.40 ≤α≤ 1.41,   that one from  3 to 5 zeros occurs
as 1.56 ≤α≤ 1.57, and so on. The last transition
 is from 15 to 17 zeros, and it just occurs
as 1.79 ≤α≤  1.80 .


    #
         
    #
      
          height 2 pt                      N_*      α      T_*   
          height 2 pt                   
       height 2 pt                      1÷ 3      1.40 ÷ 1.41      1.730÷ 5.726   
          height 2 pt                      3÷ 5      1.56÷ 1.57      8.366÷ 13.48   
          height 2 pt                   
       height 2 pt                      5÷ 7      1.64 ÷ 1.65      14.61÷ 20.00   
          height 2 pt                   
       height 2 pt                      7÷ 9      1.69 ÷ 1.70      20.80÷ 26.33   
          height 2 pt                   
       height 2 pt                      9÷ 11      1.72 ÷ 1.73      27.03÷ 32.83   
          height 2 pt                   
       height 2 pt                      11÷ 13      1.75 ÷ 1.76      33.11÷ 38.81   
          height 2 pt                   
       height 2 pt                      13÷ 15      1.78 ÷ 1.79      39.49÷ 45.51   
          height 2 pt                   
       height 2 pt                      15÷ 17      1.79 ÷ 1.80      45.51÷ 51.46   
          height 2 pt                    7  Table I



N_* = number of zeros, α = fractional order
 T_* location of the largest zero.




§ THE FUNCTIONS OF THE  WRIGHT TYPE
  

The classical Wright function, that we denote by W_λ , μ(z), is defined by the series representation convergent on the whole complex plane
,

    -0.2truecm
    W_λ , μ(z) =∑_n=0^∞z^n/n!Γ (λ n + μ),  λ > -1,  μ∈.
    (6.1)

One of its integral representations 
for λ > -1,  μ∈
reads as: 

    W_λ , μ(z) = 1/2π i∫_Ha ^σ + zσ ^-λdσ/σ ^μ, 
    (6.2)

where, as usual, Ha denotes the Hankel path. 

Then, 
W_λ , μ(z) is  an entire function for all
λ∈ (-1, +∞).
 Originally, in 1930's  Wright assumed λ≥ 0 in connection with his investigations on the asymptotic theory of partitions
<cit.>, 
  and only in 1940 <cit.>   he considered
  -1 < λ < 0.
  

We note that in the Vol 3, Chapter 18 of the handbook of the Bateman Project
<cit.>, presumably for a misprint,
 the parameter λ  is restricted to be non-negative,
whereas  the Wright functions remained  practically ignored   in other  handbooks.  
In 1990's   Mainardi,   being aware only  of  the Bateman handbook,
proved that the Wright function is entire also for -1<λ<0  in his approaches to the time fractional diffusion equation, 
see <cit.>. 


In view of the asymptotic representation in the complex domain 
and of the Laplace transform for positive argument  z=r>0
(r can denote the time variable t or the positive space variable x)
  the Wright functions are distinguished in
   first kind (λ≥ 0) and second kind
 (-1< λ < 0) 
 as outlined in the Appendix F of the book by Mainardi
 <cit.>.
It is possible to prove that the Wright function is entire of order
1/(1+λ) , hence
of exponential type if λ≥ 0 .,
that is only for the Wright functions of the first kind. 
The case λ =0 is trivial
since
W_0, μ(z) = ^  z/ Γ(μ ) .

Recurrence relations
Some of the properties,  that the Wright functions  share with
the most popular Bessel functions,
were enumerated by Wright himself.
Hereafter, we quote some relevant relations from
the handbook of Bateman Project Handbook  <cit.>: 


    -0.1truecm
    λ z   W_λ ,λ + μ(z) 
    = W_λ ,μ -1(z) + (1-μ )  W_λ ,μ(z) ,
       (6.3)


    d /  dz  W_λ ,μ(z) =
        W_λ ,λ +μ(z)  . (6.4)


We note that these relations can easily be derived from the 
series or integral representations,  (6.1) or (6.2).
 
Generalization of the Bessel functions.


 For λ =1 and μ =ν +1≥ 0 the Wright functions
(of the first kind) turn out to be related to the well known
Bessel functions J_ν  and I_ν  by the  identities:

    [ J_ν(z) = (z/2 )^ν 
          W_1, ν + 1(- z^2/4) ,;  I_ν(z) = (z/2 )^ν 
          W_1, ν + 1( z^2/4) . ](6.5)

In view of this property  some authors refer to the Wright function as
the  Wright  generalized  Bessel function
(misnamed also as the Bessel-Maitland function)
and introduce the notation

    [                                              J_ν ^(λ) (z)                    :=
    (z/2 )^ν 
          W_λ, ν + 1(- z^2/4);                                                           = (z/  2  )^ν∑_n=0^∞(-1)^n (z/2)^2n/  n! Γ(λ  n + ν +1) , ](6.6)

	with  λ >0 and ν >-1.
In particular	J_ν ^(1) (z) := J_ν  (z).

As a matter of fact, the Wright functions (of the first kind)
appear as the natural generalization of the  entire functions
 known as Bessel - Clifford functions,  see 
 Kiryakova <cit.>,
 and referred  by Tricomi <cit.>  
  
 as  the uniform Bessel functions, see also
Gatteschi <cit.>.
 

  Similarly we can properly define I_ν ^(λ) (z). 



 

 §.§ The Mainardi auxiliary functions
 
 
 
  
  
  
   
  
We note that  two particular Wright functions of the second kind, 
were  introduced by Mainardi in 1990's 
<cit.>
named F_ν (z) and M_ν (z) (0<ν<1),
called  
auxiliary functions in virtue of their role in the time fractional  diffusion equations.
 These functions
 
 are  indeed special cases of the Wright function of the second kind
  W_λ , μ(z) by setting, respectively,
  λ = -ν and μ = 0 or μ = 1-ν. 
  Hence we have:

    F_ν(z) := W_-ν , 0(-z),    0 < ν < 1,
    (6.7)

and

    M_ν(z) := W_-ν , 1- ν(-z),   0 < ν < 1,
    ((6.8)

 
These functions are interrelated through the following relation: 

    F_ν(z) = ν z M_ν(z).
    (6.9)




The series and integral representations of the auxiliary functions are derived from those of the general Wright functions. Then
for z∈ and 0<ν<1 we have:

    [                            F_ν(z)-0.2truecm                  =  ∑_n=1^∞(-z)^n/n!Γ(-ν n);                                  -0.2truecm = 
    -1/π∑_n=1^∞(-z)^n/n!Γ(ν n + 1)sin(πν n), ](6.10)



and 

    [                            M_ν(z)-0.2truecm          = ∑_n=0^∞(-z)^n/n!Γ[-ν n + (1-ν )];                                  -0.2truecm = 1/π∑_n=1^∞(-z)^n-1/(n-1)!Γ(ν n)sin(πν n), ](6.11)



The second series  representations in (6.10)-(6.11) 
 have been obtained by using the well-known reflection formula for the Gamma function[In the RHS of Eq (6.10) we have corrected an error already present since the  Appendix F of my 2010 book  <cit.>.]  
 
    Γ(ζ) Γ(1-ζ)  =π /sin πζ .

For the integral representation we have


    F_ν(z) := 1/2π i∫_Ha^σ - zσ ^νdσ, 
    (6.12)


and 

    M_ν(z) := 1/2π i∫_Ha^σ - zσ ^νdσ/σ ^1-ν. 
    (6.13)


As usual, the equivalence of the series and integral representations
is easily proved   using the Hankel formula for the Gamma
function and performing a term-by-term integration.	


Explicit expressions of F_ν(z) and M_ν(z) in terms of known functions are expected for some particular values of ν
as shown and recalled in 
<cit.>, that is

    M_1/2(z) = 1/√(π)^-z^2 /4,
    (6.14)




    M_1/3(z) = 3^2/3(z/3^1/3),. 
    (6.15)


Liemert and Klenie <cit.> have added the 
following expression for ν=2/3

    [                     M_2/3(z) 
    -0.2truecm                       = 3^-2/3 ^-2z^3/27;                               -0.2truecm [3^1/3 z (z^2/3^4/3) -3^' ( z^2/3^4/3) ] ](6.16)


where   and ^' denote the Airy function
and its first derivative.

Furthermore they  have 
suggested in the positive real field  ^+ the following remarkably integral representation

    [                         M_ν(x) 
    -0.2truecm                        =
    1/π x^ν/(1-ν)/1-ν;                                 -0.2truecm ·∫_0^π C_ν(ϕ) exp(-C_ν(ϕ)) 
    x^1/(1-ν)  dϕ, ](6.17)


where

    C_ν(ϕ) = sin(1-ν)/sinϕ ( sinνϕ/sinϕ)^ν/(1-ν) ,
    (6.18)
 


corresponding to equation (7) of the article written by Saa and Venegeroles
<cit.> .
Furthermore, it can be  proved, see <cit.>
 that
M_1/q(z)  satisfies the   differential equation 
of order q-1

    d^q-1/  dz^q-1  M_1/q(z) +
      (-1)^q/  q  z  M_1/q(z) =0 , (6.18)

subjected to the q-1 initial conditions at z=0,
derived from (6.15),

    M_1/q^(h)(0) =
       (-1)^h/π Γ[(h + 1)/q] sin [π  (h+1)/q]
       ,   
       (6.19)

  with h = 0, 1, …  q-2.
We note that, for q≥ 4 , Eq. (6.18) is  akin to the
hyper-Airy  differential equation  of order q-1 ,
see [Bender & Orszag 1987].






We find it convenient to show the plots of the M-Wright functions 
on a space symmetric interval of   in Figs 1, 2, corresponding to the cases 
0≤ν≤ 1/2 and 
1/2 ≤ν≤ 1, respectively.
   We recognize 
the non-negativity of the M-Wright function on  for
1/2 ≤ν≤ 1  consistently with the  analysis
 on distribution of zeros and asymptotics of Wright functions
 carried out by Luchko, see <cit.>, <cit.>.





 

 §.§ Laplace transform pairs related to the Wright function

Let us consider the Laplace transform of the Wright function
using the usual notation

    W_λ ,μ (± r)     ÷ ∫_0^∞^  -s r  W_λ ,μ (± r)  dr
        ,

where r denotes a non negative real variable

and s is the Laplace complex parameter.



When  λ > 0  the series representation of the Wright function can
be transformed term-by-term. In fact, for a known theorem
of the theory of the Laplace transforms, see Doetsch (194)
<cit.>,
the  Laplace transform of an entire function
of exponential type can be obtained  by transforming term-by-term the
Taylor expansion of the
original function around the origin.
In this case the resulting Laplace transform turns out
to be   analytic and vanishing at infinity.



As a consequence, we obtain  the Laplace transform pair
for |s| >  0

    W_λ ,μ (± r)  ÷ s  E_λ ,μł(±s)̊ ,
       λ > 0 ,  (6.20)

where E_λ ,μ denotes  the Mittag-Leffler function
in two parameters.  
The proof is  straightforward noting that

    ∑_n=0^∞(± r)^n/  n! Γ (λ n+μ ) ÷ s ∑_n=0^∞(± 1/s)^n/Γ (λ n+μ ) ,

and recalling the series representation
of the  Mittag-Leffler function,

    E_α, β (z ) := ∑_n=0^∞z^n/Γ(α n+ β   ) ,
     α >0 ,  β∈ .





For λ→ 0^+ Eq. (6.20) provides the Laplace transform pair
for |s| >  0,

    [    W_0^+,μ (± r) =
    ^ ± r/Γ(μ ); ÷Γ(μ)  s∓ 1
      =  s  E_0 ,μ(±s) ] ,
         (6.21)

where, to remain in agreement with (6.20), we have formally put,

    E_0,μ (z) := ∑_n=0^∞z^n /Γ(μ  )
    := Γ(μ)  E_0(z) :=
      Γ(μ) 1-z .

We recognize that in this special case the

Laplace transform exhibits a simple pole at
s = ± 1 while for λ >0 it exhibits an essential
singularity at s = 0 .



For -1<λ <0 the Wright function  turns out to be
an entire function of order greater than 1, so that
care is  required in establishing the existence
of its  Laplace transform, which necessarily must tend to
zero as s →∞ in its half-plane of convergence.
For the sake of convenience we limit ourselves to derive the
Laplace transform  for the special case of M_ν (r) ;
the exponential decay   as r →∞  of the original
function  provided by (6.20) ensures the existence of
the image function.
From the integral representation (6.13) of the M_ν function we obtain

    M_ν (r)  ÷  2π i ∫_0^∞^  -s  r [ ∫_Ha^σ -rσ ^ν dσ/σ^1-ν]  dr


    =  2π i ∫_Ha^σ σ^ν-1 [
     ∫_0^∞^  - r(s+σ ^ν )  dr ]  dσ


    = 2π i ∫_Ha^σ σ^ν-1/σ ^ν +s  dσ .

 
 

Then, by recalling the integral representation  of the Mittag-Leffler
function (3.4),
  
    E_α (z) =
      2π i ∫_Haζ ^α -1 ^ ζ/ζ ^α -z 
            dζ ,   α  >0 ,
       z∈ ,

we obtain the Laplace transform pair

    [ M_ν (r):= W_-ν, 1-ν(-r);  ÷  E_ν (-s),   0<ν<1,. ](6.22)

In this case, transforming term-by-term
 the Taylor series 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 Laplace
transform pair, see Doetsch <cit.>  and 
Eq. (3.7):

    -0.2truecm
    [                M_1/2(r)
     -0.2truecm              := 
    √(π) exp(- r^2/ 4);                          -0.2truecm ÷  E_1/2 (-s) = exp(s^2)  erfc( s), ](6.23)

  valid ∀ s ∈

 


Analogously, using the more general integral representation  (6.2) of the standard Wright function, 
we can prove that  in the case λ =-ν∈ (-1,0) and Re (μ) >0, we get

    W_-ν ,μ(-r)  ÷ 
       E_ν  , μ +ν(-s) ,   0<ν  <1 . (6.24)


In the limit as ν→ 0^+ (thus λ→ 0^-)
we formally obtain the Laplace transform pair

    [  W_0^-, μ(-r)
     -0.2truecm           :=
     ^  -r/Γ(μ );                -0.2truecm ÷Γ(μ ) s+1 :=  E_0 ,μ(-s) ](6.25)

Therefore, as λ→  0^± , and μ =1 we note a sort of continuity
in the  results (6.21)  and (6.25)    since
 
    W_0,1(-r) :=^-r÷(s+1)(6.26)
 
 with 
 
    (s+1)={[                   -0.2truecm (1/s)  E_0(-1/s) ,    |s|>1;;                   -0.2truecm         E_0(-s) ,    |s|<1 . ].
    (6.27)




We here point out the relevant  Laplace transform pairs related to the
 auxiliary functions of argument r^-ν
 with  0<ν <1,    see for details the cited author's papers  

    r  F_ν( 1/r^ν) =
       ν/  r^ν +1   M_ν( 1/r^ν) ÷ ^ -s^ν , (6.28)



    ν  F_ν( 1/r^ν) =
       1/  r^ν   M_ν( 1/r^ν) ÷ ^  -s^ν/s^1-ν . (6.29)

 We recall that the Laplace transform pairs  in  (6.28) were formerly considered
 by Pollard (1946) <cit.>,



Later  Mikusinski (1959 )<cit.> got a similar result  based on his  theory of operational calculus, 
and finally, albeit unaware of  the previous results,
  Buchen & Mainardi (1975)
<cit.>  
   derived the result in a formal way.
We note, however, that all these Authors   were not informed about the Wright functions.
Aware of the Wright functions was Stankovic <cit.> 
who in 1970 gave a rigorous proof of the Laplace transform pairs involving the Wright functions with first negative parameter, here referred of the second kind,







Hereafter we like to provide  two independent
proofs of (6.28) carrying out
the inversion of exp (-s^ν) , either by  the complex Bromwich integral
formula or  by  the formal series method. Similarly we can act for the Laplace transform
pair (6.29).
  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 σ = s r.
Recalling  (6.13)-(6.14), we get  

    L^-1 [  exp( -s^ν)] =
        2π i ∫_Br^  sr - s^ν    ds

 
    =
        2π i  r  ∫_Ha^ σ -(σ/r) ^ν    dσ


    =   r  F_ν( 1/r^ν)
      =  ν/  r^ν +1   M_ν( 1/r^ν)
      .

 Expanding in power series the Laplace transform and inverting  term by term we formally get, after recalling
 (6.12)-(6.13):  

    L^-1 [ exp( -s^ν)] =
       ∑_n=0^∞(-1)^n  /  n! L^-1 [ s^ν n]

  
    =  ∑_n=1^∞(-1)^n  /  n!  r^-ν n-1/Γ(-ν n)


    = 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 the Wright function M_ν(x) 
 for x>0 and  the complete monotonicity
 of the Mittag-leffler functions E_ν(-x) for x>0 and 0<ν<1. 
 In fact, since exp( -s^ν)
 denotes the Laplace transform of a probability density 
 (precisely, the extremal Lévy stable density of index ν, 
 see [Feller (1971)]),
 the L.H.S. of (6.28) must be non-negative, and so also the L.H.S of F(24). 
 As a matter of fact the Laplace transform pair (6.24) shows, replacing s by x,
  that the spectral representation of the Mittag-Leffler function E_ν(-x) 
  is expressed in terms of the M-Wright function M_ν(r), that is
  for x≥ 0
 
    -0.1truecm 
     E_ν(-x) = ∫_0^∞^  -rx  M_ν(r)  dr , 
     0<ν<1.  (6.30)

 We now recognize that Eq. (6.30) is consistent  with a result
  derived by Pollard (1948) <cit.>.
  
  
  
  It is instructive to compare the spectral representation of E_ν(-x) with that
 of the function E_ν(-t^ν).
 We recall  for t≥ 0,
 
    -0.1truecm
     E_ν(-t^ν)= ∫_0^∞^  -rt  K_ν(r)  dr, 
     0<ν<1, 
      (6.31)

 with spectral function

    [                          K_ν (r) 
     -0.2truecm = 
    π r^ν -1 sin (νπ)/r^2ν + 2 r^ν cos (νπ)+ 1;                                    -0.2truecm    = 
      π r sin (νπ)/r^ν + r^-ν+ 2 cos (νπ) . ](6.32)

The relationship between M_ν(r) and K_ν(r) is worth to be explored. 
Both  functions are non-negative,  integrable and normalized in ^+,
 so they  can be adopted in probability theory as  density functions.      




The transition K_ν(r)→δ(r-1)  for ν→ 1 
is easy to be detected numerically in view of the explicit  representation  (6.32). 
On the contrary, the analogous transition 
M_ν(r)→δ(r-1)  
is quite a difficult matter in view of its series and integral representations.
In this respect see the figure hereafter carried out
in the 1997 paper by Mainardi and Tomirotti <cit.>.











Here we have compared  the cases (a) ϵ =0.01 , (b) ϵ =0.001 ,
 obtained by an asymptotic method originally 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 representations of the Wright functions for parameter λ =-ν not close to the singular case 
ν=1.
 

 §.§ The Asymptotic representations

For the  asymptotic analysis in the whole complex plane
for the Wright functions, the interested reader is referred to 
Wong and Zhao (1999a),(1999b)<cit.>,
who
 have considered  the two   cases λ≥ 0 and -1<λ <0
separately,   including a description of   Stokes' discontinuity and its smoothing.
 For the Wright functions of the second kind,  where 
 λ =-ν∈ (-1,0) ,   we recall the asymptotic expansion  
originally obtained by Wright himself, 
that is  valid in a suitable
sector about the negative real axis as |z| →∞,

    -0.2truecm
    [                    W_-ν ,μ(z)
     -0.2truecm                    =
        Y^  1/2-μ ^  - Y;                                -0.2truecm ×[
      ∑_m=0^M-1 A_m    Y^-m + O(|Y|^-M)] , ](6.33)

with 

    Y = Y(z) =(1-ν )  (-ν ^ν  z)^1/(1-ν ) , 
     (6.34)

where the A_m are certain real numbers.



Let us first  point out    the asymptotic behaviour 
of the function M_ν(r) as r →∞. 



Choosing as
a variable r/ν rather than r, the computation of the requested 
asymptotic representation  by the saddle-point approximation
yields,  see Mainardi & Tomirotti (1994) <cit.>,

    [      M_ν (r/ν )
      -0.2truecm ∼
       a(ν )  r^(ν -1/2)/(1-ν);                   -0.2truecm      ×exp[- b(ν) r^1/(1-ν)], ](6.35)

where a(ν) and b(ν) are positive coefficients

    a(ν) = √(2π (1-ν)),  
      b(ν) = 1-ν/ν    >0. (6.36)

The above evaluation is consistent with the first term in
Wright's asymptotic expansion (6.33)  after having used the definition (6.36).
 We  point out that in the  limit ν→ 1^- the function M_ν(r) tends to the
   Dirac  function δ(r-1), but in a non-symmetric way as shown in 
   the two plots in   figure 11 of the previous subsection.
 


  

§ THE WRIGHT FUNCTION IN PROBABILITY THEORY

  Using the known completely monotone functions, the technique of the Laplace transform, and the Bernstein theorem, one can prove non-negativity of some Wright  functions. Say, the function

    p_ν,μ(r) = Γ(μ)  W_-ν,μ-ν(-r)
    (7.1)

can be interpreted as a one-sided probability density function (pdf)
 for 0<ν≤ 1, ν≤μ 
(see <cit.>). To show this, we use the Laplace transform pair
 (6.24)  that we rewrite in the form

    W_-ν,μ-ν(-r)
     ÷ E_ν,μ (-s),   0<ν≤ 1, (7.2)
 
and the fact that the Mittag-Leffler function E_ν,μ (-s) is completely monotone for 0<ν≤ 1, ν≤μ. According to the Bernstein theorem, the function p_ν,μ(r) is non-negative. 


To calculate the integral of p_ν,μ(r) over ^+ let us mention that it can be interpreted as the Laplace   transform of p_ν,μ at the point s=0 or the Mellin transform at s=1. 
Using the Mellin integral transform of the Wright function as in 
<cit.> leads now to the following chain of equalities:

    [                         ∫_0^∞ p_ν,μ(r)   dr-0.2truecm                        = 
    ∫_0^∞Γ(μ)  W_-ν,μ-ν(-r)  dr;                                            -0.2truecm =  
     . Γ(μ) Γ(s)/Γ(μ-ν + ν s)|_s=1 
    = Γ(μ) /Γ(μ) = 1. ]

The Mellin transform technique allows us to calculate also all
moments of order s>0 of the pdf p_ν,μ(r)  on 
^+:

    [ -0.5truecm
    ∫_0^∞  p_ν,μ(r)  r^s   dr-0.2truecm         = 
    ∫_0^∞Γ(μ) W_-ν,μ-ν(-r)  r^s+1-1  dr;                                     -0.2truecm                   =   
    Γ(μ) Γ(s+1)/Γ(μ + ν s). ](7.3)

For μ = 1, the pdf p_ν,μ(r)  can be expressed in terms of the 
M-Wright function M_ν(r), 0<ν<1 defined by 
Eq. (6.8). As it is well known (see, e.g.,
<cit.>), 
 M_ν(r) can be interpreted as a one-sided pdf on
 ^+ with the moments given by the formula

    ∫_0^∞ M_ν(r)  r^s  dr = Γ(s+1) Γ(1 +ν s),
       s>0.
    (7.4)






 §.§ The Mainardi auxiliary functions 

 as extremal stable densities

We find it 
worthwhile to recall   the relations between the Mainardi auxilary functions   and the extremal  Lévy stable densities as proven
 in the 1997 paper by Mainardi and Tomirotti <cit.>.
For an essential account of the general theory of 
 Lévy stable distributions in probability
In the present  paper the interested reader may consult the Appendix A in the present paper. More details can be found
in  the 1997  E-print by Mainardi-Gorenflo-Paradisi
<cit.> and in the 2001 paper by Mainardi-Luchko-Pagnini. <cit.>, 
recalled in the Appendix F in <cit.>.
Then, from a comparison between  the series expansions of stable densities 
 according to the Fekller-Takayasu canonic form 
with index of stability α∈ (0,2] and 
skewness θ  (|θ| ≤min{α, 2-α}),
  and  those of the Mainardi  auxiliary functions in 
 Eqs. (6.7) - (6.8), 
we recognize, see also <cit.>,  
that the 
auxiliary functions  are related to the extremal stable densities
as follows

    [            -0.1truecm
    L_α ^-α (x)
     -0.2truecm =  x   F_α  (x^-α) =
     α x^α  +1   M_α (x^-α),;                                    -0.2truecm                           0<α  <1,     x ≥ 0. ](7.5)
 

    -0.1truecm
    [           L_α ^α  -2(x)
    -0.2truecm = x   F_1/α(x) =
     1 α   M_1/α(x) ,;                         -0.2truecm               1<α≤ 2 ,    -∞<x<+∞. ](7.6)

 In the above equations, for α=1, the skewness parameter turns out to be
  θ =-1, 
 so  we get the singular limit
 
    L_1^-1(x)= M_1(x)= δ(x-1) .
     (7.7)


Hereafter we show the plots the extremal stable densities according to their expressions in terms of the M-Wright functions, see   
Eq. (<ref>), Eq. (<ref>)
for α=1/2 and
 α=3/2, respectively.
 
	 
	 
	
	
	
	





 §.§ The plots and the Fourier transform 

of the symmetric M-Wright  function
  

 We point out  that the most relevant applications of   our auxiliary functions, 
   are when the variable is real.
   In particular we consider the case of the symmetric M-Wright function
   as a function of the variable |x| for all  with varying its  parameter
   ν∈ [0,1] because related to the fundamental solution of the Cauchy 
   problem of the time fractional diffusion-wave equation dealt in Appendix B
   
	 In the following Figs. 14  and 15   
we compare the plots of the  M_ν (|x|)-Wright functions
in  |x| ≤ 5 for some rational values  in the ranges  ν∈ [0,1/2]
and ν∈ [1/2, 1], respectively.
To gain more insight of the effect of the variation of the parameter ν
we will adopt bot linear and logarithmic scales for the ordinate.
Thus in Fig. 14 
 we see the transition from exp (-|x|) for ν=0
to 1/√(π) exp (-x^2) for ν=1/2, whereas
in Fig. 15 we 
see the transition from   1/√(π) exp (-x^2) 
to the delta functions δ(x± 1) for ν=1.

In plotting M_ν (|x|) at fixed ν for sufficiently large |x|
the asymptotic representation (6.34)-(6.35) is  useful
 since, as |x| increases,
the numerical convergence of the series in (6.11)
becomes poor and poor up to being completely inefficient:
henceforth, the matching between the series and the asymptotic
representation  is relevant.  
However, as ν→ 1^-,
the plotting remains a very difficult task because
of the high peak arising around x= ± 1.












  



The Fourier transform of the  M-Wright function.
The Fourier  transform of  the symmetric (and normalized) M-Wright 
function provides its characteristic function useful in Probability theory.  

    [                                                                    ℱ[ 2 M_ν(|x|)]
    ≡2 M_ν(|x|);                                                               := 
    2∫_-∞^+∞^iκ x  M_ν(|x|)  dx;                                                = ∫_0^∞cos (κ x)   M_ν(x)   dx = E_2ν (-κ^2) . ](7.9)

 For this prove it is sufficient to develop in series the cosine function and use the formula  for  the absolute moments of the  M-Wright function in ^+. 

    [                                                                          ∫_0^∞cos (κ x)  M_ν (x)   dx;                                                     =
    ∑_n=0^∞ (-1)^n κ^2n/(2n)! ∫_0^∞x^2n  M_ν(x)  dx;                                                    = 
    ∑_n=0^∞ (-1)^n κ^2n/Γ(2ν n +1) = E_2ν,1(-κ^2) . ]

 We also have 
 
    [                                                                                            ∫_0^∞sin (κ x)  M_ν (x)   dx;                                                                  =∑_n=0^∞ (-1)^n κ^2n+1/(2n+1)! ∫_0^∞x^2n+1  M_ν(x)  dx;                                                             = 
    ∑_n=0^∞ (-1)^n κ^2n+1/Γ(2ν n +1+ν) = κ E_2ν, 1+ν(-κ^2) . ]




 §.§ Subordination formulas

We now consider  M-Wright functions as spatial probability densities evolving in time 
with self-similarity, that is

    M_ν(x,t):=t^-ν M_ν(x t^-ν)  ,   x,t≥ 0  . (7.10)

These M-Wright functions are relevant for their composition rules proved  
 by Mainardi et al. in <cit.>,  and more generally 
 in  <cit.>
 by using the Mellin Transforms. 
The main statement can be summarized as: 
Let M_λ(x;t), M_μ(x;t)
and M_ν(x;t) be M-Wright functions
of  orders
λ, μ, ν∈ (0,1)
 respectively,
then the following composition formula 
holds for any x, t ≥ 0:

    -0.2truecm
    M_ν(x,t) =∫_0^∞
    M_λ(x;τ)  M_μ(τ;t)
      dτ,  ν = λ μ.
    (7.11)




The above equation  is also intended as a   subordination formula because it can be used  
to define subordination among   self-similar stochastic processes 
(with independent increments), that properly generalize the most popular Gaussian processes,
to which they reduce for ν =1/2. 
These more general processes are governed by time-fractional diffusion equations,
as shown in   papers of our research group, 
see the survey paper by Mainardi, Mura and Pagnini <cit.> and references therein.

 

These general processes are referred to as Generalized grey Brownian Motions, that include
both Gaussian Processes  (standard Brownian motion, fractional Brownian motion)
and non-Gaussian Processes (Schneider's  grey Brownian motion),
to which 
the interested reader is referred for details.







§ CONCLUSIONS

 

In this  survey we have outlined  the basic properties of the Mittag-Leffler and  Wright functions. We have also  considered  a number of applications 
taking into account  special functions of these families.
We have stressed their relations with fractional calculus.
In particular, we have added a number of tutorial appendices to enlarge the fields
of applicability of  the Wright functions, nowadays less known than the Mittag-Leffler functions to which they are related through integral transforms.  
 



 

§ ACKNOWLEDGMENTS
 
The research activity of the author
is carried out in the framework of the activities of the National Group of Mathematical Physics (GNFM, INdAM), Italy.




§ APPENDIX A: LÉVY STABLE DISTRIBUTIONS





The term stable has been assigned by the French  mathematician Paul Lévy,
 who in the 1920's years started a systematic research 
in order to generalize the celebrated Central Limit Theorem to 
probability distributions  with infinite variance.  

For stable distributions we can assume the following 


 Definition:
If two independent real random variables
with the same shape or type of distribution are combined linearly and
the distribution of the resulting random variable has  the same shape,
the common distribution (or its type, more precisely) is said to be
stable.


The restrictive condition of stability enabled Lévy (and then other authors) to derive
the canonical form for the Fourier transform of the densities of these distributions. 
Such transform in probability theory is known  as characteristic function. 



Here we follow the parameterization adopted in  Feller (1971)
<cit.>
 revisited in 1998 by  Gorenflo & Mainardi 
<cit.>  and  popularized in 2001  
 by Mainardi, Luchko & Pagnini
 <cit.>.
 Denoting by  L_α^θ(x) a (strictly) stable density  in , 
 where α is the index of stability and  
 anθ the asymmetry parameter, improperly called skewness, 
 its characteristic function reads: 

    [                     L_α^θ(x)-0.2truecm     ÷L_α ^θ(κ)
    =   exp[- ψ_α ^θ(κ )] ,;                                       ;                             -0.2truecm ψ_α  ^θ(κ ) =
        |κ|^α ^  i (κ)θπ/2 , ] (A.1)

 with 

    0<α≤ 2 ,  
     |θ| ≤ min {α  ,2-α} .
     (A.2).


 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,0), (1,1), (1,-1), (2,0),
that we call the Feller-Takayasu diamond,
see Fig.16.
For values of θ on the border of the diamond
(that is θ = ±α if 0<α  < 1, and θ = ± (2-α ) if 1<α  <2)
we obtain the so-called extremal stable densities.








  We note the symmetry relation
L_α ^θ (-x)=	L_α ^-θ (x), so that a stable density with θ=0 is symmetric
Stable distributions have noteworthy properties of which the interested reader 
  can be informed from the existing literature. 
 Here-after we recall some  peculiar 
 Properties:
 - The class of  stable distributions possesses its own  domain of attraction, 
see Feller (1971)<cit.>.
- Any stable  density  is unimodal and indeed bell-shaped, its n-th derivative has exactly n zeros in , see 
Gawronski (1984) <cit.> and  
Simon (2015) <cit.>.
 - The stable distributions are  self-similar and  infinitely divisible.
These properties derive from the canonical form (A.1)-(A.2) through the scaling property of the Fourier transform.



Self-similarity means

    [                                                          L_α ^θ (x,t) ÷exp [-t ψ _α  ^θ(κ)];                                               ⟺ L_α ^θ (x,t)  = t^-1/α L_α ^θ (x/t^1/α ) ], ] (A.3)

where t is a positive parameter.
If t is time, then  L_α ^θ (x,t) is a spatial density evolving on time with self-similarity. 

Infinite divisibility means that  for every positive integer
n,  the characteristic function  can be expressed as the nth power of some characteristic function, 
so that

any stable distribution can be expressed as the	n-fold convolution of a
stable distribution of the same type. 
Indeed, taking in (A.3) θ=0, without loss of generality, we have

    [                                  ^-t|κ|^α                      = 
    `[^-(t/n)|κ|^α]^n;                                           ⟺  L_α ^0 (x,t)  = 
    [L_α ^0 (x,t/n)]^*n , ](A.4)
  
where   

    [L_α ^0 (x,t/n)]^*n :=
    L_α ^0 (x,t/n) *  … * L_α ^0 (x,t/n)

is the multiple Fourier convolution in  with n 
identical terms. 



Only in special cases 
we get well-known probability distributions.
For α  =2 (so θ =0), we recover the  Gaussian pdf, that turns out to be the
only stable density with finite variance, and more generally with finite  moments of any order
δ≥ 0. In fact

    L_2^0(x) = 1/2√(π)^  -x^2/4
     .
    (A.5)

All the other stable densities have finite absolute moments  of order δ∈ [-1, α ).

For α  =1  and  |θ| <1, we get

    -0.2truecm L_1^θ (x) =
     1/π cos (θπ/2)/[x+   sin (θπ/2)]^2 +[ cos (θπ/2)]^2 ,
     (A.6)

which  for θ=0 includes the  Cauchy-Lorentz pdf,

    L_1^0(x) = 1/π1/1+x^2
     .
    (A.7)

In the limiting cases   θ = ± 1 for α =1 we obtain
 the singular Dirac pdf's
 
    L_1^± 1(x)=δ(x ± 1) .
     (A.8)




In general we must recall the power series expansions 
provided by  Feller (1971) <cit.>.   
We restrict our attention to x>0
since the evaluations for x<0  can be obtained using the symmetry relation.

The convergent expansions of L_α ^θ (x)  (x>0) turn out to be
 for 0<α   <1 ,  |θ| ≤α :
 
    L_α  ^θ (x) =1π x ∑_n=1^∞
       (-x^-α)^n  Γ (1+ nα  ) n! sin[ nπ 2(θ -α  )],
        (A.9)

for 1<α≤ 2 ,   |θ| ≤ 2-α :

    L_α  ^θ (x) =1π x ∑_n=1^∞
       (-x)^n Γ (1+ n/α  ) n! sin[ nπ 2α(θ -α  )].
     (A.10)

From the series  (A.9) and the  symmetry relation
we note that the
extremal stable densities for 0<α   <1 are
unilateral, precisely vanishing for x>0 if θ =α,
vanishing for x<0 if θ =-α.
In particular the unilateral extremal densities L_α ^-α(x) with 0<α <1 
have as Laplace transform
exp (-s^α ). 



 From a comparison between  the series expansions in (A.9)-(A.10) and 
in (6.10)-(6.11) for the auxiliary functions F_α(x), M_α(x) 
we recognize that for x>0  
the auxiliary functions of the Wright type are related to the extremal stable densities as in Eqs. (6.28)-(6.29).








 



More generally,  all (regular) stable densities, given in Eqs. (6.38)-(6.39),  were recognized to belong 
to the class of Fox H-functions,
as formerly shown in 1986 by  Schneider 
<cit.>, 
see also Mainardi-Pagnini-Saxena (2003)
<cit.>.



  
  
 


§ APPENDIX B: THE TIME-FRACTIONAL 
 DIFFUSION EQUATION

There exist three equivalent forms of the time-fractional diffusion equation of a single order,
two with fractional derivative and one with fractional integral,

provided we refer to the standard initial condition u(x,0)= u_0(x). 
Taking a real number β∈ (0,1), the time-fractional diffusion equation 
of order β in the Riemann-Liouville sense reads  
     
    ∂ u/∂ t= K_β D^1-β_t ∂^2 u/∂ x^2 ,
    (B.1)

 in the Caputo sense reads

    _*D^β_t u = K_β ∂^2 u/∂ x^2 ,
    (B.2)

and in integral form

    u(x,t)=u_0(x)+K_β 1/Γ(β)∫_0^t 
    (t-τ)^β-1 ∂^2 u(x,τ)/∂ x^2  dτ .
    (B.3)

where K_β is a sort of fractional diffusion coefficient of dimensions 
[K_β]= [L]^2  [T]^-β= cm^2/sec^β.



The fundamental solution 
(or Green function) 𝒢_β(x,t) for  the  equivalent Eqs. (B.1) - (B.3),
that is the solution corresponding to the initial condition

    𝒢_β(x,0^+)= u_0(x) = δ(x) (B.4)

can be expressed in terms of the M-Wright function  

    𝒢_β(x,t)= 1/2√(K_β)  t^β/2  M_β/2(|x|/√(K_β)  t^β/2)  .
    (B.5)

The corresponding variance can be  promptly obtained 
 
    σ _β ^2(t) := ∫_-∞^+∞ x^2  𝒢_β(x,t)  dx =  
       2/Γ(β+1)   K_β t^β . (B.6)

 As a consequence, for 0<β <1 the variance is consistent 
 with a  process of slow diffusion with similarity exponent H=β/2.
 
 
 
 
 
 
 







The  fundamental solution _β(x,t) for the time-fractional diffusion equation can be obtained
by applying in sequence the Fourier and Laplace transforms to  any form  chosen
among Eqs. (B.1)-(B.3).
Let us devote our attention to the integral form (B.3)
using non-dimensional variables by setting K_β =1 
and adopting the notation J_t^β for the fractional integral.
Then,  our  Cauchy problem  reads

    _β(x,t)= δ(x) + J_t^β ∂^2 _β/∂ x^2(x,t) 
     .(B.7)

In  the Fourier-Laplace domain, after applying   formula for the
Laplace transform   of the fractional integral 
and observing δ(κ ) ≡ 1,
we get
 
    _β(κ ,s) = s
         -κ^2/s^β _β(κ ,s)  ,

from which     for Re (s) > 0 , κ∈

    _β(κ ,s)
       =   s^β -1/s^β + κ ^2 ,   0<β≤ 1 .
     . (B.8)

	
	
	












Strategy (S1): Recalling the Fourier transform pair

    a/b +κ^2 ÷ a/2 b^1/2 ^  -|x|b^1/2 ,   a,b>0 ,
    (B.9)

and setting a=s^β-1, b=s^β, we get

    _β (x,s)= 2s^β/2-1 ^  -|x|s^β/2 .(B.10)

 Strategy (S2): Recalling the Laplace transform pair
 
    s^β-1/s^β +c   E_β(-c t^β) ,   c>0 ,
     (B.11)

and setting c=κ^2, we  have

    _β (κ ,t) = E_β(-κ^2 t^β) .
    (B.12)

Both strategies lead to the result

    _β(x,t) = 2M_β/2(|x|,t) = 2  t^-β/2  
    M_β/2(|x|/t^β/2) , (B.13)

consistent with Eq. (B.5). 



The time-fractional drift equation
Let us finally note that the M-Wright function does appear also in the fundamental solution of 
the time-fractional drift equation. 
Writing this equation in non-dimensional form 
and adopting  the Caputo derivative we have 

    _*D^β_t u(x,t) = -  ∂/∂ x u(x,t)  , -∞<x< +∞ ,   t≥ 0 , (B.14)

where 0<β<1 and u(x,0^+) = u_0(x).
When u_0(x)= δ(x) we  obtain the fundamental solution (Green function) that we denote by
_β^*(x,t).
Following the usual approach  we show that 

    _β^*(x,t) = 
     {[ t^-β  M_β(x/t^β) ,              x>0 ,;                0 ,             x<0  , ].
      (B.15)

that for β=1 reduces to the right running pulse δ(x-t) for x>0.


In  the Fourier-Laplace domain, after applying the  formula for the
Laplace transform   of the Caputo fractional derivative 
and observing δ(κ ) ≡ 1,
we get
 
    s^β _β^*(κ ,s)- s^β-1 = +iκ _β^*(κ ,s)  ,

from which,    for Re (s) > 0,  κ∈  

    _β^*(κ ,s)
       =   s^β -1/s^β -iκ ,   0<β≤ 1 , 
     . (B.16)

To determine the  Green function _β^*(x,t)

in the space-time domain we can follow two
alternative  strategies related to the
order in carrying out the inversions in (B.16).


(S1) : invert  the Fourier transform
getting _β (x,s)
   and   then invert the remaining  Laplace transform;


(S2) : invert  the Laplace transform getting _β^* (κ ,t)
and then invert the remaining  Fourier transform.



Strategy (S1): Recalling the Fourier transform pair

    a/b -iκ  a/ b ^  -xb ,   a,b>0 ,  x>0 ,
    (B.17)

and setting a=s^β-1, b=s^β, we get

    _β^* (x,s)= s^β-1 ^  -x s^β .(B.18)

 Strategy (S2): Recalling the Laplace transform pair
 
    s^β-1/s^β +c   E_β(-c t^β) ,  c>0 ,(B.19)

and setting c=-iκ, we  have

    _β^* (κ ,t) = 
    E_β(iκ t^β) .(B.20)

Both strategies lead to the result (B.15).

In view of Eq. (7.5) we also recall that the M-Wright function is related to the unilateral 
extremal stable density of index β. Then,
using our notation  for stable densities,
 we write our Green function (B.15) as 

    _β^*(x,t) = t/β x^-1-1/β  L_β^-β(t x^-1/β) . 
    (B.21)





§ APPENDIX C:  HISTORICAL AND 
 BIBLIOGRAPHIC NOTES
  
In the early nineties, precisely in his 1993 former analysis of  fractional equations
describing  slow diffusion and interpolating diffusion and wave-propagation,
the present author <cit.>,
  introduced the so called auxiliary functions of the Wright type

F_ν (z) :=   W_-ν , 0(-z) and
M_ν (z) :=  W _-ν , 1-ν(-z) with 0<ν <1,

in order to characterize the fundamental solutions
for typical boundary value problems, as it is shown in the previous sections.   
 
Being then only aware of the Handbook of the  Bateman project,
where the parameter λ of the Wright function
W_λ,μ(z)
was erroneously restricted to  non-negative values,

the author thought to have originally extended the  analyticity  property  of the original Wright  function
by taking  ν= -λ  with ν∈ (0,1).


Then  the function M_ν was  referred to as the 
Mainardi functionin the 1999 treatise by  Podlubny
<cit.>
 and in some later papers and books.







 


 It was  Professor B. Stanković
during the presentation of the paper
by Mainardi & Tomirotti
<cit.> 
at the Conference Transform Methods and Special Functions, Sofia 1994,
who informed the author that this extension for -1<λ < 0
had  been already made  by Wright himself in 1940 (following
his previous papers in the thirties), see <cit.>.

In  a 2005  paper published in FCAA <cit.>,  devoted to the 80th birthday of Professor Stanković,
 the author  used the occasion
to renew  his personal gratitude to Professor
Stanković  for this earlier information that led
 him to study the original papers by Wright and to work  (also in collaboration)
on the functions of the Wright type for further applications,
 see <cit.> and <cit.>.

For the above reasons the author preferred to distinguish the Wright functions into the first kind 
(λ≥ 0)   and the second kind (-1<λ <0).


It should be noted that in the book by Prüss (1993) 
<cit.> 
we find a figure quite similar to our Fig  15-top reporting the M-Wright function in linear scale, namely the Wright function of the second kind. It was  derived from inverting the Fourier transform expressed in terms of the Mittag-Leffler function following the approach 
by Fujita <cit.>
for the fundamental  solution of the Cauchy problem for the diffusion-wave equation, fractional in time.  

However, our plot must be considered independent by that of Prüss because the author (Mainardi)  used the Laplace transform in his former paper 
 presented at WASCOM, Bologna, October  1993
<cit.>
(and published later in a number of papers and in his 2010 book)  so he  was  aware of the book by Prüss only later. 






99







Barret 54
J.H. Barret (1954).
Differential equations of non-integer order,
Canad. J. Math 6, 529–541.

Buchen-Mainardi 75
P.W.  Buchen  and  F. Mainardi  (1975).
  Asymptotic expansions for transient  viscoelastic waves,
  Journal de Mécanique 14, 597–608.

Caputo-Mainardi PAG71
 M. Caputo  and F.  Mainardi  (1971a).
   A new dissipation model based on memory mechanism,
  Pure and Appl. Geophys. (PAGEOPH) 91, 134–147.
  [Reprinted in Fract. Calc.  Appl. Anal.
  10 No 3, 309–324 (2007)]

 Caputo-Mainardi RNC71
 M. Caputo and  F. Mainardi  (1971b).
  Linear models of dissipation in  anelastic solids,
  Riv. Nuovo Cimento (Ser. II) 1, 161–198. 
  

Cole 33
 K.S. Cole  (1933).
 Electrical conductance of biological systems,
 Electrical excitation in nerves, in
Proceedings Symposium on  Quantitative Biololgy, Cold Spring
Harbor,  New York,  Vol. 1, pp. 107–116.






Consiglio-Mainardi RIM2019
A. Consiglio &   F. Mainardi (2019).
On the evolution of fractional diffusive waves.
Ricerche di Matematica, published on line 06 Dec. 2019.
 DOI: 10.1007/s11587-019-00476-6
[E-print: arxiv.org/abs/1910.1259]
  
Davis BOOK36
H.T.  Davis  (1936).
The Theory of Linear Operators,
  The Principia Press, Bloomington, Indiana. 
  
  Doetsch BOOK74
  G. Doetsch  (1974).
  Introduction to the Theory and Application of the Laplace
  Transformation,
  Springer Verlag, Berlin.

 
 Dzherbashyan BOOK66
 M.M. Dzherbashyan  (1966).
 Integral Transforms and
  Representations of Functions in the Complex Plane, 
  Nauka, Moscow  [in Russian]. 
  There is also the transliteration as Djrbashyan.

 Erdelyi BATEMAN
 
 A. Erdélyi, W. Magnus, F. Oberhettinger  and F. Tricomi, (1955).
Higher Transcendental Functions, 3-rd Volume,
 McGraw-Hill, New York . [Bateman Project]. 








 
 
 
 
 
 


Feller BOOK71
 W. Feller (1971).
An Introduction to Probability Theory and its Applications.
  Wiley, New York, Vol. II, Second Edition.
   


Fujita 90
 Y. Fujita  (1990).
  Integro-differential equation which interpolates the heat equation
 and the wave equation, I, II,
  Osaka J. Math. 27, 309–321, 797–804. 

Gatteschi BOOK73
L.G. Gatteschi (1973).
Funzioni Speciali, 
  UTET, Torino.

Gawronski 84
W.  Gawronski (1984).  
  On the bell-shape of stable distributions,
  Annals of Probability 12, 230–242. 

GKMS BOOK14
R. Gorenflo,  A.A. Kilbas, F. Mainardi & S. Rogosin (2014).
Mittag-Leffler Functions. Related Topics and Applications, Springer, Berlin. 
Second Edition in preparation.


 
 GoLuMa 99
 R. Gorenflo, Yu.  Luchko and F. Mainardi (1999).
Analytical properties and applications of the Wright function.
Fract. Calc.  Appl. Anal.   2, 383–414.


GoLuMa 00
R. Gorenflo,  Yu. Luchko   and  F.  Mainardi  (2000).
Wright functions as scale-invariant solutions of the diffusion-wave
 equation,
Journal of  Computational and   Applied Mathematics
     118, 175–191.

Gorenflo-Mainardi BERLIN95
R. Gorenflo and F. Mainardi (1996).
 Fractional oscillations and Mittag-Leffler functions,
  Preprint No A-96-14, 
  Fachbereich Mathematik und Informatik,
  Freie Universität Berlin, Serie    Mathematik, pp. 22.
  [E-print:   www.math.fu-berlin.de/publ/index.html] 

Gorenflo-Mainardi CISM97
R. Gorenflo and F. Mainardi (1997).
  Fractional calculus:
integral and differential equations of fractional order,
   in:  A. Carpinteri and  F. Mainardi  (Editors):
Fractals and Fractional Calculus in Continuum Mechanics,
 Springer Verlag, Wien, 
 pp. 223–276.
[E-print: http://arxiv.org/abs/0805.3823]


Gorenflo-Mainardi FCAA98
R. Gorenflo and F. Mainardi (1998).
Random walk models for space-fractional diffusion processes,
Fract. Calc. Appl.  Anal. 1 No 2, 167–190. 

Gross  JAP47
B.  Gross   (1947).
On creep and relaxation,
J. Appl. Phys. 18, 212–221.

Hille-Tamarkin 30
E. Hille  and J.D. Tamarkin (1930).
  On the theory of linear integral equations,
  Ann. Math.  31, 479–528.

Kiryakova BOOK94
V. Kiryakova  (1994).
Generalized Fractional Calculus and Applications,
Longman & J. Wiley, Harlow - New York. 
[Pitman Research Notes in Mathematics, Vol. 301]

Liemert-Kleine JMP2015 
 A. Liemert  and A.   Klenie  (2015).
 Fundamental Solution of the Tempered Fractional Diffusion Equation.
 J. Math. Phys  56, 113504. 
 
  Luchko 2000
  Yu. Luchko (2000). 
 On the asymptotics of zeros of the Wright function. 
 Zeitschrift für Analysis und ihre Anwendungen,
 19, 597–622.
 
 Luchko HFCA
 Yu.  Luchko (2019).
The Wright function and its applications,
in A. Kochubei, Yu.Luchko (Eds.),
Handbook of Fractional Calculus with Applications
Vol. 1: Basic Theory, pp. 241- 268. 
De Gruyter GmbH, 2019, Berlin/Boston.
 Series edited  by J. A.Tenreiro Machado.
 
 LK 2013
Yu. Luchko & V. Kiryakova (2013).
 The Mellin integral transform in fractional calculus, 
Fract. Calc. Appl. Anal. 16, 405–430.

Mainardi WASCOM93
F. Mainardi  (1994).
 On the initial value problem for the fractional diffusion-wave equation,
in:  Rionero, S. and  Ruggeri, T. (Editors),
Waves and Stability in Continuous Media.
World Scientific, Singapore, pp. 246–251.
 [Proc. VII-th  WASCOM, Int. Conf. "Waves and Stability in Continuous
   Media",  Bologna, Italy, 4-7 October 1993]



Mainardi AML96
 F. Mainardi  (1996a).
The fundamental solutions for the fractional diffusion-wave equation.
 Applied Mathematics Letters,9  No 6,  23–28.


Mainardi CSF96
 F. Mainardi (1996b).
  Fractional relaxation-oscillation and fractional
  diffusion-wave phenomena.
  Chaos, Solitons & Fractals 7, 1461–1477.

Mainardi CISM97
F. Mainardi  (1997).
  Fractional calculus:
some basic problems in continuum and statistical mechanics,
   in:   A. Carpinteri  and  F. Mainardi  (Editors):
Fractals and Fractional Calculus in Continuum Mechanics,
 Springer Verlag, Wien, 
  pp. 291–348.
  [E-print: http://arxiv.org/abs/1201.0863]
  
 






 

   Mainardi BOOK10 
F. Mainardi  (2010).
 Fractional  Calculus and Waves in Linear Viscoelasticity.
  Imperial College Press, London and World Scientific, Singapore.

 Mainardi-Consiglio MATHEMATICS2020
F. Mainardi and A. Consiglio (2020).
The Wright functions of the second kind in Mathematical Physics,
Mathematics 8 No 6, 884/1–26
.[DOI: 10.3390/MATH8060884]
[E-print arXiv:2007.02098]


Mainardi-Gorenflo FCAA07
F. Mainardi  and  R. Gorenflo  (2007).
Time-fractional derivatives in relaxation processes: a tutorial survey,
Fract. Calc.  Appl. Anal. 10,   269–308. 
[E-print: http://arxiv.org/abs/0801.4914]




Mainardi-Gorenflo-Vivoli FCAA05
 F. Mainardi, R. Gorenflo and A. Vivoli  (2005).
 Renewal processes of Mittag-Leffler and Wright type,
 Fract. Calc.  Appl. Anal. 8,  7–38.  
[E-print http://arxiv.org/abs/math/0701455]




Mainardi-Luchko-Pagnini FCAA01
F.  Mainardi, Yu.  Luchko and G.  Pagnini (2001).
 The fundamental solution of the space-time fractional diffusion equation,
   Fract. Calc.  Appl. Anal.  4, 153–192.

 [E-print: arxiv.org/abs/cond-mat/0702419]


Mainardi-Mura-Pagnini IJDE10
F. Mainardi, A. Mura and G. Pagnini (2010),
The M Wright function in time-fractional diffusion processes: a tutorial survey, International Journal of Differential Equations Vol. 2010, Article ID 104505, 29 pages. 
[DOI:10.1155/2010/104505] [E-print http://arxiv.org/abs/1004.2950]


Mainardi-Pagnini AMC03 
F. Mainardi and G. Pagnini (2003). 
The Wright functions as solutions of the time-fractional diffusion equation. 
Applied Mathematics and Computation 141 No 1, 51–62.

Mainardi-Pagnini-Gorenflo  FCAA03
F. Mainardi,  G. Pagnini  and  R. Gorenflo (2003).
Mellin transform and subordination laws  in fractional diffusion processes,
Fract. Calc.  Appl. Anal. 6 No 4, 441–459. 
[E-print: http://arxiv.org/abs/math/0702133]

Mainardi-Pagnini-Saxena JCAM03
F. Mainardi,  G. Pagnini  and R.K.  Saxena (2005).
Fox H functions in fractional diffusion,
   J. Comp. Appl. Math. 178, 321–331.

Mainardi-Paradisi-Gorenflo BUDAPEST97
F. Mainardi, P. Paradisi and R. Gorenflo (1998), 
Probability distributions generated by fractional diffusion equations,
Invited Lecture, Workshop on Econophysics, Budapest 21-27 July 1997.
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!
Mainardi-Tomirotti TMSF94
F.  Mainardi  and M.  Tomirotti (1995).
On a special function
  arising  in  the time fractional diffusion-wave equation,
  in: P. Rusev,  I. Dimovski   and  V. Kiryakova,  (Editors),
Transform Methods and Special Functions, Sofia 1994,
    Science Culture Technology, Singapore, pp. 171–183.


Mainardi-Tomirotti GEO97
 F. Mainardi and M.  Tomirotti (1997).
Seismic Pulse Propagation with Constant Q and Stable Probability
 Distributions.
Annali di Geofisica 40, 1311–1328.

[E-print: arxiv.org/abs/1008.1341]

Mikusinski 59
  J. Mikusiński  (1959).
    On the function whose Laplace transform
    is exp  (- s^α ) ,
    Studia Math. 18, 191–198.

Mittag-Leffler 03a
G.M. Mittag-Leffler (1903a).
Une généralisation de l'intégrale de Laplace-Abel,
C.R. Acad. Sci. Paris (Ser. II) 137, 537–539.

Mittag-Leffler 03b
G.M. Mittag-Leffler  (1903b).
Sur la nouvelle fonction E_α (x),
C.R. Acad. Sci. Paris (Ser. II) 137, 554–558.

Mittag-Leffler 04
G.M. Mittag-Leffler  (1904).
Sopra la  funzione E_α (x),
 Rendiconti R. Accademia Lincei (Ser. V) 13, 3–5.


 Mittag-Leffler 05
G.M. Mittag-Leffler  (1905).
Sur la représentation analytique d'une branche uniforme
  d'une fonction monogène,
   Acta Math. 29, 101–181.




Paris HFCA
R.B. Paris (2019).
Asymptotics of the special functions of
fractional calculus,
in: A. Kochubei, Yu.Luchko (Editors),
Handbook of Fractional Calculus with Applications,
Vol. 1: Basic Theory, pp. 297-325.
De Gruyter GmbH, 2019 Berlin/Boston.
 Series edited by J. A.Tenreiro Machado.

Podlubny BOOK99
I. Podlubny (1999).
Fractional Differential Equations,
 Academic Press, San Diego.


Pollard 46
      H.  Pollard  (1946).
 The representation of exp ( -x^λ)  as a Laplace
 integral,  Bull. Amer. Math. Soc. 52, 908–910.


Pollard 48
       H.  Pollard,  (1948). 
    The completely monotonic character of the Mittag-Leffler function 
    E_α (-x),
   Bull. Amer. Math. Soc. 54,  1115–1116

Pruss BOOK93
J. Prüss, (1993).
Evolutionary Integral Equations and Applications.
 Birkhauser Verlag, Basel.
 
Saa-Venegeroles PRE2011 
 A. Saa and  R. Venegeroles (2011).
 Alternative numerical computation of one-sided Lévy and Mittag-Leffler distributions,
 Phys. Rev. E  84, 026702.
 
 Sansone-Gerretsen 60
  G. Sansone  and J. Gerretsen. (1960).
Lectures on the Theory of Functions of a Complex Variable,
 Vol. I. Holomorphic Functions,
 Nordhoff, Groningen.  

Schneider LNP86
W.R. Schneider  (1986).
 Stable distributions: Fox function representation and generalization,
 in  S. Albeverio, G.  Casati   and  D. Merlini, D. (Editors),
 Stochastic Processes in Classical and Quantum Systems,
Springer Verlag, Berlin,  pp.  497–511. 
 [Lecture Notes in Physics, Vol. 262]

Simon 15
T. Simon (2015).
Positive stable densities and the bell-shape,
Proc. Amer. Math. Soc. 143 No 2, 885–895.

Stankovic BEOGRAD1970 
B. Stankovič (1970).
 On the function of E.M. Wright.
 Publ. de l’Institut Mathématique, Beograd,
Nouvelle Sér. 10, 113–124.

Tricomi BOOK59
F.G. Tricomi (1959).
Funzioni Speciali, 
  Gheroni, Torino.
  

Wiman 1905a
A. Wiman (1905a).
Über den Fundamentalsatz der Theorie der Funkntionen
E_α(x).
Acta Math. 29, 191–201.

Wiman 1905b
A. Wiman (1905b).
Über die Nullstellen der Funkntionen E_α(x).
 Acta Math. 29, 217–234.

Wong 99a
 R. Wong and Y.-Q.  Zhao   (1999a).
Smoothing of Stokes' discontinuity for the generalized Bessel function.
Proc. R. Soc. London A455,  1381–1400.

Wong 99b
R.   Wong and  Y.-Q.  Zhao (1999b).
Smoothing of Stokes' discontinuity for the generalized Bessel function II,
Proc. R. Soc. London A 455,  3065–3084.

Wright 33
 E.M. Wright  (1933).
On the coefficients of power series having exponential singularities.
 Journal London Math. Soc.,  8, 71–79.

Wright 35
E.M. Wright (1935).
The asymptotic expansion of the generalized Bessel function.
Proc. London Math. Soc. (Ser. II)  38,  257–270.






Wright 40
E.M. Wright (1940).
 The generalized Bessel function of order greater than one.
Quart. J. Math., Oxford Ser. 11, 36–48.






