The Wright functions of the second kind in Mathematical Physics

In this review paper we stress the importance of the higher transcendental Wright functions of the second kind in the framework of Mathematical Physics.We first start with the analytical properties of the classical Wright functions of which we distinguish two kinds. We then justify the relevance of the Wright functions of the second kind as fundamental solutions of the time-fractional diffusion-wave equations. Indeed, we think that this approach is the most accessible point of view for describing non-Gaussian stochastic processes and the transition from sub-diffusion processes to wave propagation. Through the sections of the text and suitable appendices we plan to address the reader in this pathway towards the applications of the Wright functions of the second kind. Keywords: Fractional Calculus, Wright Functions, Green's Functions, Diffusion-Wave Equation,


Introduction
The special functions play a fundamental role in all fields of Applied Mathematics and Mathematical Physics because any analytical results are expressed in terms of some of these functions. Even if the topic of special functions can appear boring and their properties mainly treated in handbooks, we would promote the relevance of some of them not yet so well known. We devote our attention to the Wright functions, in particular with the class of the second kind. These functions, as we will see hereafter, are fundamental to deal with some non-standard deterministic and stochastic processes. Indeed the Gaussian function (known as the normal probability distribution) must be generalized in a suitable way in the framework of partial differential equations of non-integer order. This work is organized as follows. In Section 2 we introduce the Wright functions, entire in the complex plane that we distinguish in two kinds in relation on the value-range of the two parameters on which they depend. In particular we devote our attention on two Wright functions of the second kind introduced by Mainardi with the term of auxiliary functions. One of them, known as M-Wright function generalizes the Gaussian function so it is expected to play a fundamental role in non-Gaussian stochastic processes. Indeed In Section 3 we show how the Wright functions of the second kind are relevant in the analysis of time-fractional diffusion and diffusion-wave equations being related to their fundamental solutions. This analysis leads to generalize the known results r of the standard diffusion equation in the one-dimensional case, that is recalled in Appendix A by means of auxiliary functions as particular cases of the Wright functions of the second kind known as M-Wright or Mainardi functions. For readers' convenience, in Appendix B we will also provide a introduction to the time-derivative of fractional order in the Caputo sense We remind that nowadays, as usual, by fractional order we mean a non-integer order,so that the term "fractional" is a misnomer kept only for historical reasons. In Section 4 we consider again the Mainsrdi auxiliary functions functions for their role in probability theory and in particular in the framework of Lévy stable distributions whose general theory is recalled in Appendix C. In Section 5 we show how the auxiliary functions turn out to be included in a class that we denote the four sister functions. On their turn these four functions depending on a real parameter ν ∈ (0, 1) are the natural generalization of the three sisters functions introduced in Appendix A devoted to the standard diffusion equation. The attribute of sisters was put by one of us (F. M.) because of their inter-relations, in his lecture notes on Mathematical Physics, so it has only a personal reason that we hope to be shared by the readers. Finally, in Section 6, we provide some concluding remarks paying attention to work to be done in the next future. We point out that we have equipped our theoretical analysis with several plots hoping they will be considered illuminating for the interested readers. We also note that we have limited our review to the simplest boundary values problems of equations in one space dimension referring the readers to suitable references for more general treatments in Section 3.1.

The Wright functions of the second kind and the Mainardi auxiliary functions
The classical Wright function, that we denote by W λ,µ (z), is defined by the series representation convergent in the whole complex plane, The integral representation reads as: where Ha − denotes the Hankel path: this one is a loop which starts from −∞ along the lower side of negative real axis, encircles with a small circle the axes origin and ends at −∞ along the upper side of the negative real axis. W λ,µ (z) is then an entire function for all λ ∈ (−1, +∞). Originally Wright assumed λ ≥ 0 in connection with his investigations on the asymptotic theory of partition [48,49] and only in 1940 he considered −1 < λ < 0, [50]. We note that in the Vol 3, Chapter 18 of the handbook of the Bateman Project [10], presumably for a misprint, the parameter λ is restricted to be non-negative, whereas the Wright functions remained practically ignored in other handbooks. In 1993 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, that will be dealt in a next Section. In view of the asymptotic representation in the complex domain and of the Laplace transform for positive argument z = r > 0 (r can be the time variable t or the 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 [35]. In particular, for the asymptotic behaviour, we refer the interested reader to the two papers by Wong and Zhao [46,47], and to the surveys by Luchko and by Paris in the Handbook of Fractional Calculus and Applications, see respectively [25], [41], and references therein. We note that the Wright functions are entire of order 1/(1 + λ) hence only the first kind functions (λ ≥ 0) are of exponential order whereas the second kind functions (−1 < λ < 0) are not of exponential order. The case λ = 0 is trivial since W 0,µ (z) = e z /Γ(µ). As a consequence of the distinction in the kinds, we must point out the different Laplace transforms proved e.g. in [15], [35], see also the recent survey on Wright functions by Luchko [25]. We have: for the first kind, when λ ≥ 0 for the second kind, when −1 < λ < 0 and putting for convenience ν = −λ so 0 < ν < 1 Above we have introduced the Mittag-Leffler function in two parameters α > 0, β ∈ C defined as its convergent series for all z ∈ C .
For more details on the special functions of the Mittag-leffler type we refer the interested readers to the treatise by Gorenflo et al [14], where in the forthcoming 2-nd edition also the Wright functions are treated in some detail.
In particular, two Wright functions of the second kind, originally introduced by Mainardi and named F ν (z) and M ν (z) (0 < ν < 1), are called auxiliary functions in virtue of their role in the time fractional diffusion equations considered in the next section. These functions, F ν (z) and M ν (z), are indeed special cases of the Wright function of the second kind W λ,µ (z) by setting, respectively, λ = −ν and µ = 0 or µ = 1 − ν.
Hence we have: Those functions are interrelated through the following relation: which reminds us the second relation in 7, seen for the standard diffusion equation.
The series representations of the auxiliary functions are derived from those of W λ,µ (z). Then: where it has been used in both cases the reflection formula for the Gamma function (Eq. 11) among the first and the second step of Eqs. (9) and (10), Also the integral representations of the auxiliary functions are derived from those of W λ,µ (z). Then: Explicit expressions of F ν (z) and M ν (z) in terms of known functions are expected for some particular values of ν as shown and recalled by Mainardi in the first 1990's in a series of papers [29,30,31,32], that is Liemert and Klenie [21] have added the following expression for ν = 2/3 where Ai and Ai denote the Airy function and its first derivative. Furthermore they have suggested in the positive real field IR + the following remarkably integral representation , corresponding to equation (7) of the article written by Saa and Venegeroles [42] . Let us point out the asymptotic behaviour of the function M ν (x) as x → +∞. Choosing as a variable x/ν rather than x, the computation of the asymptotic representation by the saddle-point approximation carried out by Mainardi and Tomirotti yields, see [39] and [35], The above evaluation is consistent with the first term in the original asymptotic expansion by Wright in [49,50] after having used the definition of .M -Wright function Now we find it convenient to show the plots of the M -Wright functions on a space symmetric interval of IR in Figs 1, 2, corresponding to the cases 0 ≤ ν ≤ 1/2 and 1/2 ≤ ν ≤ 1, respectively. We recognize the nonnegativity of the M -Wright function on IR for 1/2 ≤ ν ≤ 1 consistently with the analysis on distribution of zeros and asymptotics of Wright functions carried out by Luchko, see [22], [25].
where D is a positive constant whose dimensions are L 2 T −β and u = u(x, t; β) is the field variable, which is assumed again to be a causal function of time. The Caputo fractional derivative is recalled in the Appendix B so that in explicit form the TFDWE (21) splits in the following integro-differential equations: In view of our analysis we find convenient to put: We can then formulate the basic problems for the Time Fractional Diffusion-Wave Equation using a correspondence with the two problems for the standard diffusion equation. Denoting by f (x) and g(t) two given, sufficiently well-behaved functions, we define: a) Cauchy problem (26) u(x, 0 + ; ν) = 0, 0 ≤ x < +∞; u(0 + , t; ν) = g(t), u(+∞, t; ν) = 0, t > 0 If 1/2 < ν ≤ 1 corresponding to 1 < β ≤ 2 we must consider also the initial value of the first time derivative of the field variable u t (x, 0 + ; ν), since in this case Eq. (21) turns out to be akin to the wave equation and consequently two linear independent solutions are to be determined. However, to ensure the continuous dependence of the solutions to our basic problems on the parameter ν in the transition from ν = (1/2) − to ν = (1/2) + , we agree to assume u t (x, 0 + ; ν) = 0.
For the Cauchy and Signalling problems, following the approaches by Mainardi, see e.g. [29] and related papers, we introduce now the Green functions G c (x, t; ν) and G s (x, t; ν) that for both problems can be determined by the LT technique, so extending the results known from the ordinary diffusion equation. We recall that the Green functions are also referred to as the fundamental solutions, corresponding respectively to f (x) = δ(x) and g(t) = δ(t) with δ(·) is the Dirac delta generalized function The expressions for the Laplace Transforms of the two Green's functions are: Now we can easily recognize the following relation: which implies for the original Green functions the following reciprocity relation for x > 0'and t > 0 and 0 < ν < 1: where z is the similarity variable and F ν (z) and M ν (z) are the Mainardi auxilary functions introduced in the previous section. Indeed Eq. (30) is the generalization of Eq. (A.8) that we have seen for the standard diffusion equation due to the introduction of the time fractional derivative of order ν Then, the two Green functions of the Cauchy and Signalling problems turn out to be expressed in terms of the two auxiliary functions as follows.
For the Cauchy problem we have For the Signalling problem we have: that generalizes Eq. (A.7).

4.1.
Complements to the time-fractional diffusion-wave equations. The boundary value problems dealt previously can be considered with a source data function f x) and g(t) different from the Dirac generalized functions, in particular with box-type functions as it has been carried out recently by us, see [8].
The TFDWE can be generalized in 2D and 3D space dimensions. so consequently the Wright functions play again a fundamental role. However, we prefer to refer the interested reader to the literature, in particular to the papers by Luchko and collaborators [22,23,24,25], [27,28], [1], by Hanyga [18] and to the recent analysis by Kemppainen [19]. All of them are originated in some way from the seminal paper by Schneider & Wyss [43]. In some of these papers the authors have considered also fractional differentiation both in time and in space, so that they have generalized to more than one dimension the former analysis by Mainardi, Luchko & Pagnini [37] on the space-time fractional diffusion-wave equations.

The M −Wright functions in probability theory and the stable distributions
We recognize that the Wright M -function with support in IR + can be interpreted as probability density function (pdf ) because it is non negative and also it satisfies the normalization condition: We now provide more details on these densities in the framework of the theory of probability. Fundamental quantities about the Wright M −function are the absolute moments of order δ > −1 in R + , that are finite and turn out to be: The result is based on the integral representation of the M −Wright function: The exchange between two integrals and the following identity contributed to the final result for Eq. (35): In particular, for δ = n ∈ IN, the above formula provides the moments of integer order. Indeed recalling the Mittag-Leffler function introduced in Eq. (5) with α = ν and β = 1: the moments of integer order can also be computed from the Laplace transform pair proved in the Appendix F of [35] as follows:

The auxiliary functions versus 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 [40]. For readers' convenience we refer to Appendix C for an essential account of the general Lévy stable distributions in probability. Indeed, from a comparison between the series expansions of stable densities in (C.8)-(C.9) and of the auxiliary functions in Eqs. (9) -(10), we recognize that the auxiliary functions are related to the extremal stable densities as follows In the above equations, for α = 1, the skewness parameter turns out to be θ = −1, so we get the singular limit Hereafter we show the plots the extremal stable densities according to their expressions in terms of the M -Wright functions, see Eq. (40)., Eq. (41) for α = 1/2 and α = 3/2, respectively. We recognize that

5.2.
The symmetric M − Wright function. We easily recognize that extending the function M ν (x) in a symmetric way to all of IR (that is putting x = |x|) and dividing by 2 we have a symmetric pdf with support in all of IR. As the parameter ν changes between 0 and 1, the pdf goes from the Laplace pdf to two half discrete delta pdfs passing for ν = 1/2 through the Gaussian pdf.
To develop a visual intuition, also in view of the subsequent applications, we show the plots of the symmetric M −Wright function on the real axis at t = 1 for some rational values of the parameter ν ∈ [0, 1]  The readers are invited to look the YouTube video by Consiglio whose title is "Simulation of the M −Wright function", in which the author shows the evolution of this function as the parameter ν changes between 0 and 0.85 in a finite interval of IR centered in x = 0. Finally we compute the characteristic function for the symmetric Eq. 43 is obtained by developing in series the cosine function and using Eq. 35. In particular: 5.3. The Wright M-function in two variables. In view of timefractional diffusion processes related to time-fractional diffusion equations it is worthwhile to introduce the function in two variables From Eq. (55) we derive the Fourier transform of M ν (|x|, t) with respect to x ∈ IR, Using the Mellin transforms, Mainardi et al. [38] derived the following interesting integral formula of composition, Special cases of the Wright M-function are simply derived for ν = 1/2 and ν = 1/3 from the corresponding ones in the complex domain, see Eqs. (28)- (29). We devote particular attention to the case ν = 1/2 for which we get the Gaussian density in IR, For the limiting case ν = 1 we obtain

The four sisters
In this section we show how some Wright functions of the second kind can provide an interesting generalization of the three sisters discussed in Appendix A. The starting point is a (not well-known) paper published in 1970 by Stankovic [45], where (in our notation) the following Laplace transform pair is proved rigorously: where x and t are positive. We note that the Stankovic formula can be derived in a formal way by developing the exponential function in positive power of s and inverting term by term as described in the Appendix F of the book by Mainardi [35]. We recognize that the Laplace Transforms of the Three Sisters functions φ(x, s), ψ(x, s) and χ(x, s) are particular cases of the Eq. (52) for ν = 1/2, that is of according to the following scheme: -φ(x, s) with µ = 1, -ψ(x, s) with µ = 0, -χ(x, s) with µ = 1/2. If ν is no longer restricted to ν = 1/2 we define Four Sisters functions as follows Here we show some plots of these functions, both in the t and in the x domain for some values of ν (ν = 1/4, 1/2, 3/4). Note that for ν = 1/2 we only find three functions, (the Three Sisters functions) of Appendix A

Conclusions
In our survey on the Wright functions we have distinguished two kinds, pointing out the particular class of the second kind. Indeed these functions have been shown to play key roles in several processes governed by non Gaussian processes, including sub-diffusion, transition to wave propagation, Lévy stable distributions. Furthermore, we have devoted our attention to four functions of this class that we agree to called the Four Sisters functions. All these items justify the relevance of the Wright functions of the second kind in Mathematical Physics.

Appendix A: The standard diffusion equation and the three sisters
In this Appendix let us recall the Diffusion Equation in the onedimensional case where the constant D > 0 is the diffusion coefficient coefficient, whose dimensions are L 2 T −1 and x , t denote the space and time coordinates, respectively. Two basic problems for Eq. (7) are the Cauchy and Signalling ones introduced hereafter In these problems some initial values and boundary conditions are set; specify the values attained by the field variable and/or by some of its derivatives on the boundary of the space-time domain is an essential step to guarantee the existence, the uniqueness and the determination of a solution of physical interest to the problem, not only for the Diffusion Equation.
Two data functions f (x) and g(t) are then introduced to write formally these conditions; some regularities are required to be satisfied by f (x) and g(t), and in particular f (x) must admit the Fourier transform or the Fourier series expansion if the support is finite, while h(t) must admit the Laplace Transform. We also require without loss of generality that the field variable u(x, t) is vanishing for t < 0 for every x in the spatial domain. Given these premises, we can specify the two aforementioned problems.
In the Cauchy problem the medium is supposed to be unlimited (−∞ < x < +∞) and to be subjected at t = 0 to a known disturbance provided by the data function f (x). Formally: This is a pure initial-value problem (IVP) as the values are specified along the boundary t = 0. In the Signalling problem the medium is supposed to be semi-infinite (0 ≤ x < +∞) and initially undisturbed. At x = 0 (the accessible end) and for t > 0 the medium is then subjected to a known disturbance provided by the causal function g(t). Formally: This problem is referred to as an initial boundary value problem (IBVP) in the quadrant {x, t} > 0. For each problem the solutions turn out to be expressed by a proper convolution between the data functions and the Green functions G, that are the fundamental solutions of the problems. For the Cauchy problem we have: For the Signalling problem we have: Following the lecture notes in Mathematical Physics by Mainardi [34], we note that the following relevant property is valid for {x, t} > 0: According to Mainardi' s notations, Eq. (7) is known as reciprocity relation, F (z) and M (z) are called auxiliary functions and z is the similarity variable.
A particular case of the Signalling problem is obtained when g(t) = H(t) (the Heaviside unit step function) and the solution u(x, t) turns out to be expressed in terms of the complementary error function: where the sign ÷ is used for the juxtaposition of a function with its Laplace transform. We easily note that Eq. (7) is related to the Step-Response problem, Eq. (7) is related to the Signalling problem and Eq. (7) is related to the Cauchy problem. Following the lecture notes by Mainardi [34] we agree to call the above functions the three sisters functions for their role in the standard diffusion equation. They will be discussed with details hereafter. Everything that we have said above will be found again as a special case of the Time Fractional Diffusion Equation where the time derivative of the first order is replaced by a suitable time derivative of non-integer order.
It is easy to demonstrate that each of them can be expressed as a function of one of the 2 others three sisters (table 1). Table 1. Relations among the three sisters in the Laplace domain.
The three sisters in the t domain may be all directly calculated by making use of the Bromwich formula taking account of the contribution of the branch cut of √ s and of the pole of 1/s. Wee obtain: Then, through the substitution ρ = √ r, we arrive at the Gaussian integral and, consequently, we find the previous explicit expressions of the three sisters, that is: reminding the definition of the complementary error function. Alternatively, we can compute the three sisters in t domain by using the relations among the three sisters in the Laplace domain listed in table 1. But in this case one of the three sisters in t domain must be already known. Assuming to know φ(a, t) from Eq. (A.11), we get: -ψ(a, t) from ψ(a, s) = s φ(a, s). Indeed, noting since φ(a, 0 + ) = 0 we can obtain (A.12), namely where a is seen as a parameter,.
Indeed it immediately follows (A.13), namely For more details we refer the reader again to [34].

Appendix B: Essentials of Fractional Calculus
Fractional calculus is the field of mathematical analysis which deals with the investigation and applications of integrals and derivatives of arbitrary order. The term fractional is a misnomer, but it is retained for historical reasons, following the prevailing use. This appendix is based on the 1997 surveys by Gorenflo and Mainardi [16] and by Mainardi [33]. For more details on the classical treatment of fractional calculus the reader is referred to the nice and rigorous book by Diethelm [9] published in 2010 by Springer in the series Lecture Notes in Mathematics. According to the Riemann-Liouville approach to fractional calculus, the notion of fractional integral of order α (α > 0) is a natural consequence of the well known formula (usually attributed to Cauchy), that reduces the calculation of the n−fold primitive of a function f (t) to a single integral of convolution type. In our notation the Cauchy formula reads where IN is the set of positive integers. From this definition we note that f n (t) vanishes at t = 0 with its derivatives of order 1, 2, . . . , n − 1 .
For convention we require that f (t) and henceforth f n (t) be a causal function, i.e. identically vanishing for t < 0 .
In a natural way one is led to extend the above formula from positive integer values of the index to any positive real values by using the Gamma function. Indeed, noting that (n − 1)! = Γ(n) , and introducing the arbitrary positive real number α , one defines the Fractional Integral of order α > 0 : where IR + is the set of positive real numbers. For complementation we define J 0 := I (Identity operator), i.e. we mean J 0 f (t) = f (t) . Furthermore, by J α f (0 + ) we mean the limit (if it exists) of J α f (t) for t → 0 + ; this limit may be infinite.
We note the semigroup property J α J β = J α+β , α , β ≥ 0 , which implies the commutative property J β J α = J α J β , and the effect of our operators J α on the power functions These properties are of course a natural generalization of those known when the order is a positive integer.
Introducing the Laplace transform by the notation L {f (t)} := ∞ 0 e −st f (t) dt = f (s) , s ∈ C , and using the sign ÷ to denote a Laplace transform pair, i.e. f (t) ÷ f (s) , we note the following rule for the Laplace transform of the fractional integral, which is the generalization of the case with an n-fold repeated integral.
After the notion of fractional integral, that of fractional derivative of order α (α > 0) becomes a natural requirement and one is attempted to substitute α with −α in the above formulas. However, this generalization needs some care in order to guarantee the convergence of the integrals and preserve the well known properties of the ordinary derivative of integer order.
Denoting by D n with n ∈ IN , the operator of the derivative of order n , we first note that D n J n = I , J n D n = I , n ∈ IN , i.e. D n is leftinverse (and not right-inverse) to the corresponding integral operator J n . In fact we easily recognize from (B.1) that As a consequence we expect that D α is defined as left-inverse to J α . For this purpose, introducing the positive integer m such that m − 1 < α ≤ m , one defines the Fractional Derivative of order α > 0 as (B.6) Defining for complementation D 0 = J 0 = I , then we easily recognize that D α J α = I , α ≥ 0 , and Of course, these properties are a natural generalization of those known when the order is a positive integer. Note the remarkable fact that the fractional derivative D α f is not zero for the constant function f (t) ≡ 1 if α ∈ IN . In fact, (B.7) with γ = 0 teaches us that This, of course, is ≡ 0 for α ∈ IN, due to the poles of the gamma function in the points 0, −1, −2, . . . . We now observe that an alternative definition of fractional derivative was introduced by Caputo in 1967 [3] in a geophysical journal and in 1969 [4] in a book in Italian. Then the Caputo definition was adopted in 1971 by Caputo and Mainardi [5,6] in the framework of the theory of Linear Viscoelasticity. Nowadays it is usually referred to as the Caputo fractional derivative and reads (B.9) We note that there are a number of discussions on the priority of this definition that surely was formerly considered by Liouville as sated by Butzer and Westphal [2]. However Liouville did not recognize the relevance of this representation derived by a trivial integration by part whereas Caputo, even if unaware of the Riemann-Liouville representation, promoted his definition in several papers over all for the applications where the Laplace transform plays a fundamental role. We agree to denote (B.9) as the Caputo fractional derivative to distinguish it from the standard Riemann-Liouville fractional derivative (B.6). The Caputo definition (B.9) is of course more restrictive than the Riemann-Liouville definition (B.6), in that requires the absolute integrability of the derivative of order m. Whenever we use the operator D α * we (tacitly) assume that this condition is met. We easily recognize that in general unless the function f (t) along with its first m − 1 derivatives vanishes at t = 0 + . In fact, assuming that the passage of the m-derivative under the integral is legitimate, one recognizes that, for m − 1 < α < m and t > 0 , and therefore, recalling the fractional derivative of the power functions (B.7) , The alternative definition (B.9) for the fractional derivative thus incorporates the initial values of the function and of its integer derivatives of lower order. The subtraction of the Taylor polynomial of degree m−1 at t = 0 + from f (t) means a sort of regularization of the Riemann-Liouville fractional derivative. In particular for 0 < α < 1 we get . According to the Caputo definition, the relevant property for which the fractional derivative of a constant is still zero can be easily recognized, We now explore the most relevant differences between the two fractional derivatives (B.6) and (B.9). We observe, again by looking at (B.7), that D α t α−1 ≡ 0 , α > 0 , t > 0 . From above we thus recognize the following statements about functions which for t > 0 admit the same fractional derivative of order α , In these formulas the coefficients c j are arbitrary constants.
For the two definitions we also note a difference with respect to the formal limit as α → (m − 1) + . From (B.6) and (B.9) we obtain respectively, We now consider the Laplace transform of the two fractional derivatives. For the standard fractional derivative D α the Laplace transform, assumed to exist, requires the knowledge of the (bounded) initial values of the fractional integral J m−α and of its integer derivatives of order k = 1, 2, . . . , m − 1 . The corresponding rule reads, in our notation, The Caputo fractional derivative appears more suitable to be treated by the Laplace transform technique in that it requires the knowledge of the (bounded) initial values of the function and of its integer derivatives of order k = 1, 2, . . . , m − 1 , in analogy with the case when α = m . In fact, by using (B.4) and noting that (B.19) we easily prove the following rule for the Laplace transform, Indeed, the result (B.20), first stated by Caputo by using the Fubini-Tonelli theorem, appears as the most "natural" generalization of the corresponding result well known for α = m .
In particular Gorenflo and Mainardi have pointed out the major utility of the Caputo fractional derivative in the treatment of differential equations of fractional order for physical applications. In fact, in physical problems, the initial conditions are usually expressed in terms of a given number of bounded values assumed by the field variable and its derivatives of integer order, no matter if the governing evolution equation may be a generic integro-differential equation and therefore, in particular, a fractional differential equation.

Appendix C: The Lévy stable distributions
We now introduce the so-called Lévy stable distributions.. The term stable has been assigned by the French mathematician Paul Lévy, who, in the tuenties of the last century, 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 canonic form for the characteristic function of the densities of these distributions. Here we follow the parameterization by Feller [11,12] revisited by Gorenflo & Mainardi in [17], see also [37]. Denoting by L θ α (x) a generic stable density in IR, where α is the index of stability and and θ the asymmetry parameter, improperly called skewness, its characteristic function reads: 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 Figure 10. 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 relevant existing literature. Here-after we recall some peculiar Properties: -The class of stable distributions possesses its own domain of attraction, see e.g. [12].
-Any stable density is unimodal and indeed bell-shaped, i.e. its n-th derivative has exactly n zeros in IR, see Gawronski [13], Simon [44] and Kwaśnicki [20]. -The stable distributions are self-similar and infinitely divisible. These properties derive from the canonic form (C.1) through the scaling property of the Fourier transform. Self-similarity means 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 (C.1) θ = 0, without loss of generality, we have where L 0 α (x, t/n) * n := L 0 α (x, t/n) * L 0 α (x, t/n) * · · · * L 0 α (x, t/n) is the multiple Fourier convolution in IR with n identical terms.
Only for a few particular cases, the inversion of the Fourier transform in (C.1) can be carried out using standard tables, and well-known probability distributions are obtained. 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 (C. 6) In the limiting cases θ = ±1 for α = 1 we obtain the singular Dirac pdf 's L ±1 1 (x) = δ(x ± 1) . (C.7) In general, we must recall the power series expansions provided in [12]. 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 , |θ| ≤ α : From the series in (C.8) 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 support in IR + and Laplace transform exp(−s α ). For α = 1/2 we obtain the so-called Lévy-Smirnov pdf : As a consequence of the convergence of the series in (C.8)-(C.9) and of the symmetry relation we recognize that the stable pdf 's with 1 < α ≤ 2 are entire functions, whereas with 0 < α < 1 have the form where Φ 1 (z) and Φ 2 (z) are distinct entire functions.The case α = 1 (|θ| < 1) must be considered in the limit for α → 1 of (C.8)-(C.9), because the corresponding series reduce to power series akin with geometric series in 1/x and x, respectively, with a finite radius of convergence. The corresponding stable pdf 's are no longer represented by entire functions, as can be noted directly from their explicit expressions (C.5)-(C.6). We omit to provide the asymptotic representations of the stable densities referring the interested reader to Mainardi et al (2001) [37]. However, based on asymptotic representations, we can state as follows; for 0 < α < 2 the stable pdf 's exhibit fat tails in such a way that their absolute moment of order δ is finite only if −1 < δ < α. More precisely, one can show that for non-Gaussian, not extremal, stable densities the asymptotic decay of the tails is For the extremal densities with α = 1 this is valid only for one tail (as |x| → ∞), the other (as |x| → ∞) being of exponential order. For 1 < α < 2 the extremal pdf 's are two-sided and exhibit an exponential left tail (as x → −∞) if θ = +(2 − α) , or an exponential right tail (as x → +∞) if θ = −(2 − α) . Consequently, the Gaussian pdf is the unique stable density with finite variance. Furthermore, when 0 < α ≤ 1, the first absolute moment is infinite so we should use the median instead of the non-existent expected value in order to characterize the corresponding pdf .