Why the Mittag-Leffler Function Can Be Considered the Queen Function of the Fractional Calculus?

In this survey we stress the importance of the higher transcendental Mittag-Leffler function in the framework of the Fractional Calculus. We first start with the analytical properties of the classical Mittag-Leffler function as derived from being the solution of the simplest fractional differential equation governing relaxation processes. Through the sections of the text we plan to address the reader in this pathway towards the main applications of the Mittag-Leffler function that has induced us in the past to define it as the Queen Function of the Fractional Calculus. These applications concern some noteworthy stochastic processes and the time fractional diffusion-wave equation We expect that in the next future this function will gain more credit in the science of complex systems. Finally, in an appendix we sketch some historical aspects related to the author’s acquaintance with this function.


Introduction
For a few decades, the special transcendental function known as the Mittag-Leffler function has attracted the increasing attention of researchers because of its key role in treating problems related to integral and differential equations of fractional order.
Since its introduction in 1903-1905 by the Swedish mathematician Mittag-Leffler at the beginning of the last century up to the 1990s, this function was seldom considered by mathematicians and applied scientists.
Before the 1990s, from a mathematical point of view, we recall the 1930 paper by Hille and Tamarkin [1] on the solutions of the Abel integral equation of the second kind, and the books by Davis [2], Sansone & Gerretsen [3], Dzherbashyan [4] (unfortunately in Russian), and finally Samko et al. [5]. Particular mention would be for the 1955 Handbook of High Transcendental Functions of the Bateman project [6], where this function was treated in Volume 3, in the chapter devoted to miscellaneous functions. For former applications we recall an interesting note by Davis [2] reporting previous research by Dr. Kenneth S. Cole in connection with nerve conduction, and the papers by Cole & Cole [7], Gross [8] and Caputo & Mainardi [9,10], where the Mittag-Leffler function was adopted to represent the responses in dielectric and viscoelastic media. More information are found in the Appendix of this survey.
In the 1960's the Mittag-Leffler function started to exit from the realm of miscellaneous functions because it was considered as a special case of the general class of Fox H functions, that can exhibit an arbitrary number of parameters in their integral Mellin-Barnes representation. However, in our opinion, this classification in a too general framework has, to some extent, obscured the relevance and the applicability of this function in applied sciences. In fact, most mathematical models are based on a small number of parameters, say 1 or 2 or 3, so that a general theory may be confusing whereas the adoption of a generalized Mittag-Leffler function with 2 or 3 indices may be sufficient.
Nowadays it is well recognized that the Mittag-Leffler function plays a fundamental role in Fractional Calculus even if with a single parameter (as originally introduced by Mittag-Leffler) just to be worth of being referred to as the Queen Function of Fractional Calculus, see Mainardi & Gorenflo [11]. We find some information on the Mittag-Leffler functions in any treatise on Fractional Calculus but for more details we refer the reader to the surveys of Haubold, Mathai and Saxena [12] and by Van Mieghem [13] and to the treatise by Gorenflo et al. [14], just devoted to Mittag-Leffler functions, related topics and applications.
The plan of this survey is the following. We start to give in Section 2 the main definitions and properties of the Mittag-Leffler function in one parameter with related Laplace transforms. Then in Section 3 we describe its use in the simplest fractional relaxation equation pointing out its compete monotonicity. The asymptotic properties are briefly discussed in Section 4. In Section 5 we briefly discuss the so called generalized Mittag-Leffler function, that is the 2-parameter Mittag-Leffler function. Of course further generalization to 3 and more parameter will be referred to specialized papers and books. Then in the following sections we discuss the application of the Mittag-Leffler function in some noteworthy stochastic processes. We start in Section 6 with the fractional Poisson process, and then in Section 7 with its application of the thinning of renewal processes. The main application are dealt in Section 8 where we discuss the continuous time random walks (CTRW) and then in Section 9 we point out the asymptotic universality. In Section 10 we discuss the time fractional diffusion-wave processes pointing out the role of the Mittag-Leffler functions in two parameters and their connection with the basic Wright functions. In Appendix A we find worthwhile to report the acquaintance of the author with the Mittag-Leffler functions started in the late 1960s and continued up to nowadays.
We recall that Sections 3-10 are taken from several papers by the author, published alone and with colleagues and former students. Furthermore we have not considered other applications of the Mittag-Leffler functions including, for example, anomalous diffusion theory in terms of fractional and generalized Langevin equations. On this respect we refer the readers to the articles of the author, see References [15,16], and to the recent book by Sandev and Tomovski [17] and references therein. For many items related to the Mittag-Leffler functions we refer again to the treatise by Gorenflo et al. [14].
We recognize that it is an entire function of order 1/α providing a simple generalization of the exponential function exp(z) to which it reduces for α = 1. For detailed information on the Mittag-Leffler-type functions and their Laplace transforms the reader may consult e.g., [6,18,19] and the recent treatise by Gorenflo et al. [14].
We also note that for the convergence of the power series in (1) the parameter α may be complex provided that (α) > 0. The most interesting properties of the Mittag-Leffler function are associated with its asymptotic expansions as z → ∞ in various sectors of the complex plane.
In this paper we mainly consider the Mittag-Leffler function of order α ∈ (0, 1) on the negative real semi-axis where is known to be completely monotone (CM) due a classical result by Pollard [20], see also Feller [21].
Let us recall that a function φ(t) with t ∈ IR + is called a completely monotone (CM) function if it is non-negative, of class C ∞ , and (−1) n φ (n) (t) ≥ 0 for all n ∈ IN. Then a function ψ(t) with t ∈ IR + is called a Bernstein function if it is non-negative, of class C ∞ , with a CM first derivative. These functions play fundamental roles in linear hereditary mechanics to represent relaxation and creep processes, see, for example, Mainardi [22]. For mathematical details we refer the interested reader to the survey paper by Miller and Samko [23] and to the treatise by Schilling et al. [24].
In particular we are interested in the function whose Laplace transform pair reads Here we have used the notation ÷ to denote the juxtaposition of a function of time f (t) with its Laplace transform The pair (3) can be proved by transforming term by term the power series representation of e α (t) in the R.H.S of (2). Similarly we can prove the following Laplace transform pair for its time derivative For this Laplace transform pair we can simply apply the usual rule for the Laplace transform for the first derivative of a function, that reads

The Mittag-Leffler Function in Fractional Relaxation Processes
For readers' convenience let us briefly outline the topic concerning the generalization via fractional calculus of the first-order differential equation governing the phenomenon of (exponential) relaxation. Recalling (in non-dimensional units) the initial value problem the following two alternatives with α ∈ (0, 1) are offered in the literature: where D 1−α t and * D α t denote the fractional derivative of order 1 − α in the Riemann-Liouville sense and the fractional derivative of order α in the Caputo sense, respectively.
For a generic order µ ∈ (0, 1) and for a sufficiently well-behaved function f (t) with t ∈ IR + the above derivatives are defined as follows, see for example, Gorenflo and Mainardi [18], Podlubny [19], Between the two derivatives we have the relationship Both derivatives in the limit µ → 1 − reduce to the standard first derivative but for In analogy to the standard problem (5), we solve the problems (7) and (8) with the Laplace transform technique, using the rules pertinent to the corresponding fractional derivatives, that we recall hereafter for a generic order µ ∈ (0, 1), We note that it is generally more cumbersome to use the Laplace transform pair for the Rieaman Liouville derivative (13) that for the Capute derivative (14). Indeed the rule (13) requires the initial value of the fractional integral of f (t) whereas the rule (14) simply requires the initial value of f (t). For this property the Caputo derivative is mostly used in physical problems where finite initial values are given.
Then we recognize that the problems (a) and (b) are equivalent since the Laplace transform of the solution in both cases comes out as that yields, in virtue of the Laplace transform pair (3) We thus recognize that the Mittag-Leffler function provides the solution to the fractional relaxation equation, as outlined, for example, by Gorenflo and Mainardi [18], Mainardi and Gorenflo [11], and Mainardi [22].
Furthermore, by anti-transforming the R.H.S of (3) by using the complex Bromwich formula, and taking into account for 0 < α < 1 the contribution from branch cut on the negative real semi-axis (the denominator s α + 1 does nowhere vanish in the cut plane −π ≤ arg s ≤ π), we get, see the survey by Gorenflo and Mainardi [18], where We note that this formula was obtained as a simple exercise of complex analysis without be aware of the Titchmarsh formula for inversion of Laplace transforms [25], revised by Gross and Levi [26] and by Gross [27]. This formula is rarely outlined in books on Laplace transforms so we refer the reader for example to Apelblat's book [28] for its presence.
Since K α (r) is non-negative for all r in the integral, the above formula proves that e α (t) is a CM function in view of the Bernstein theorem. This theorem provides a necessary and sufficient condition for a CM function as a real Laplace transform of a non-negative measure.
However, the CM property of e α (t) can also be seen as a consequence of the result by Pollard [20] because the transformation x = t α is a Bernstein function for α ∈ (0, 1). In fact it is known that a CM function can be obtained by composing a CM with a Bernstein function based on the following theorem: Let φ(t) be a CM function and let ψ(t) be a Bernstein function, then φ[ψ(t)] is a CM function.
As a matter of fact, K α (r) provides an interesting spectral representation of e α (t) in frequencies.
With the change of variable τ = 1/r we get the corresponding spectral representation in relaxation times, namely that can be interpreted as a continuous distributions of elementary (i.e., exponential) relaxation processes. As a consequence we get the identity between the two spectral distributions, that is a surprising fact pointed out in Linear Viscoelasticity by the author in his book [22]. This kind of universal/scaling property seems a peculiar one for our Mittag-Leffler function e α (t).
In Figure 1, we show K α (r) for some values of the parameter α. Of course for α = 1 the Mittag-Leffler function reduces to the exponential function exp(−t) and the corresponding spectral distribution is the Dirac delta generalized function centred at r = 1, namely δ(r − 1). In Figure 2, we show some plots of e α (t) for some values of the parameter α. It is worth to note the different rates of decay of e α (t) for small and large times. In fact the decay is very fast as t → 0 + and very slow as t → +∞. The Mittag-Leffler function turns out the basic function in relaxation processes of physical interested as occurring in viscoelastic and dielectric materials. We refer the readers for viscoelasticity, that is, to the contribution of the author including References [22,29,30] whereas for dielectric materials to the survey by Garrappa et al. [31]. For the pioneers who have pointed out the role of the Mittagf-Leffler function in mechanical and dielectric relaxation processes we refer to the recent survey by Mainardi and Consiglio [32].

Asymptotic Approximations to the Mittag-Lefler Function
We now report the two common asymptotic approximations of our Mittag-Lefflrer function. Indeed, it is common to point out that the function e α (t) matches for t → 0 + with a stretched exponential with an infinite negative derivative, whereas as t → ∞ with a negative power law. The short time approximation is derived from the convergent power series representation (2). In fact, The long time approximation is derived from the asymptotic power series representation of e α (t) that turns out to be, see [6] so that, at the first order, As a consequence the function e α (t) interpolates for intermediate time t between the stretched exponential and the negative power law. The stretched exponential models the very fast decay for small time t, whereas the asymptotic power law is due to the very slow decay for large time t. In fact, we have the two commonly stated asymptotic representations: The stretched exponential replaces the rapidly decreasing expression 1 − t α /Γ(1 + α) from (21). Of course, for sufficiently small and for sufficiently large values of t we have the inequality In Figures 3 and 4, we compare for α = 0.25, 0.5, 0.75, 0.90 in logarithmic scales the function e α (t) (continuous line) and its asymptotic representations, the stretched exponential e 0 α (t) valid for t → 0 (dashed line) and the power law e ∞ α (t) valid for t → ∞ (dotted line). We have chosen the time range 10 −5 ≤ t ≤ 10 +5 .
We note from Figures 3 and 4 that, whereas the plots of e 0 α (t) remain always under the corresponding ones of e α (t), the plots of e ∞ α (t) start above those of e α (t) but, at a certain point, an intersection may occur so changing the sign of the relative errors. The interested reader may consul the plots of the relative errors in the 2014 paper by the author [33] from which, in particular, Figures 1-4 have been extracted.

The Generalized Mittag-Leffler Function
In this survey we will devote our attention mainly to the classical Mittag-Leffler in one parameter α as introduced by Mittag-Leffler since 1903 defined by the power series in (1). We have just learned from the instructive E-print by Van Mieghem [13] that the series (1) was discussed by Hadamard in 1893, that is 10 years earlier than Mittag-Leffler himself.
As a matter of fact a straightforward generalization of the classical Mittag-Leffler function is obtained by replacing the additive constant 1 in the argument of the Gamma function in (1) by an arbitrary complex parameter β . It was formerly considered in 1905 by Reference [34] and soon later by Mittag-leffler himself, almost incidentally in one of his notes. Later, in the 1950's, such generalization was investigated by Humbert and Agarwal, with respect to the Laplace transformation, see References [35][36][37]. Usually, when dealing with Laplace transform pairs, the parameter β is required to be real and positive like α.
For this function we agree to use the notation Of course E α,1 (z) ≡ E α (z). The series is still convergent for all the complex plane C so the function (26) is still entire for (α) > 0 for any β ∈ C with order 1/ α so the additional parameter play any role on this respect. However the Laplace transform pairs concerning the Mittag-Leffler function (26) and its derivative are known to be with α, β > 0 and (s) > |λ| 1α , see, for example, Refs. [14,19,22], and We also note the following relation concerning the first derivative of the classical Mittag-Leffler function with the two-parameter Mittag-Leffler function usually overlooked by several authors but of easy proof: We report the plot of the function φ α (t) herewith in Figure 5. We note that Mittag-Leffler functions with more than two parameters were also dealt by several authors as pointed out in [14]. In particular, for the 3-parameter Mittag-Leffler function (known as Prabhakar function) and related operators we refer the reader to the recent survey by Giusti et al. [38] and references therein. Kiryakova has dealt in a number of papers the multi-index Mittag-Leffler functions, see for example [39].

The Fractional Poisson Process and the Mittag-Leffler Function
Hereafter we describe how the Mittag-Leffler function enters in the so-called fractional Poisson process. We are following the original approach by Mainardi et al. in [40] where the fractional Poisson process is referred to as the renewal process of the Mittag-Leffler type. However, an independent approach to the fractional Poisson process was given for example, by Laskin in [41].

Essentials of Renewal Theory
The concept of renewal process has been developed as a stochastic model for describing the class of counting processes for which the times between successive events are independent identically distributed (iid) non-negative random variables, obeying a given probability law. These times are referred to as waiting times or inter-arrival times. For more details see, for example, the classical treatises by Cox [42], Feller [21].
For a renewal process having waiting times T 1 , T 2 , . . . , let That is t 1 = T 1 is the time of the first renewal, t 2 = T 1 + T 2 is the time of the second renewal and so on. In general t k denotes the kth renewal.
The process is specified if we know the probability law for the waiting times. In this respect we introduce the probability density function (pd f ) φ(t) and the (cumulative) distribution function Φ(t) so defined: When the non-negative random variable represents the lifetime of technical systems, it is common to refer to Φ(t) as to the failure probability and to as to the survival probability, because Φ(t) and Ψ(t) are the respective probabilities that the system does or does not fail in (0, T]. A relevant quantity is the counting function N(t) defined as that represents the effective number of events before or at instant t. In particular we have Ψ(t) = P (N(t) = 0) . Continuing in the general theory we set thus F k (t) represents the probability that the sum of the first k waiting times is less or equal t and f k (t) its density. Then, for any fixed k ≥ 1 the normalization condition for F k (t) is fulfilled because In fact, the sum of k random variables each of which is finite with probability 1 is finite with probability 1 itself. By setting for consistency F 0 (t) ≡ 1 and f 0 (t) = δ(t), where for the Dirac delta generalized function in IR + we assume the formal representation we also note that for k ≥ 0 we have We now find it convenient to introduce the simplified * notation for the Laplace convolution between two causal well-behaved (generalized) functions f (t) and g(t) Being f k (t) the pd f of the sum of the k iid random variables T 1 , . . . , T k with pd f φ(t) , we easily recognize that f k (t) turns out to be the k-fold convolution of φ(t) with itself, so Equation (36) simply reads: Because of the presence of Laplace convolutions a renewal process is suited for the Laplace transform method. Throughout this paper we will denote by f (s) the Laplace transform of a sufficiently well-behaved (generalized) function f (t) according to and for δ(t) consistently we will have δ(s) ≡ 1 . Note that for our purposes we agree to take s real. We recognize that (38) reads in the Laplace domain where, using (32),

The Classical Poisson Process as a Renewal Process
The most celebrated renewal process is the Poisson process characterized by a waiting time pd f of exponential type, The process has no memory. Its moments turn out to be and the survival probability is We know that the probability that k events occur in the interval of length t is The probability distribution related to the sum of k iid exponential random variables is known to be the so-called Erlang distribution (of order k). The corresponding density (the Erlang pd f ) is thus so that the Erlang distribution function of order k turns out to be In the limiting case The results (44)- (46) can easily obtained by using the technique of the Laplace transform sketched in the previous section noting that for the Poisson process we have: and for the Erlang distribution: We also recall that the survival probability for the Poisson renewal process obeys the ordinary differential equation (of relaxation type)

The Renewal Process of Mittag-Leffler Type
A "fractional" generalization of the Poisson renewal process is simply obtained by generalizing the differential Equation (49) replacing there the first derivative with the integro-differential operator * D β t that is interpreted as the fractional derivative of order β in Caputo's sense, see Section 2. We write, taking for simplicity We also allow the limiting case β = 1 where all the results of the previous section (with λ = 1) are expected to be recovered. For our purpose we need to recall the Mittag-Leffler function as the natural "fractional" generalization of the exponential function, that characterizes the Poisson process. We again recall that the Mittag-Leffler function of parameter β is defined in the complex plane by the power series as stated in Section 2 where the parameter was denoted by α.
The solution of Equation (50) is known to be, see Section 3 so Then, the corresponding Laplace transforms read Hereafter, we find it convenient to summarize the most relevant features of the functions Ψ(t) and φ(t) when 0 < β < 1 . We begin to quote their series expansions convergent in all of IR suitable for t → 0 + and their asymptotic representations for t → ∞, and In contrast to the Poissonian case β = 1, in the case 0 < β < 1 for large t the functions Ψ(t) and φ(t) no longer decay exponentially but algebraically. As a consequence of the power-law asymptotics the process turns be no longer Markovian but of long-memory type. However, we recognize that for 0 < β < 1 both functions Ψ(t), φ(t) keep the "completely monotonic" character of the Poissonian case. as can be simply derived from Section 2. We recall that complete monotonicity of our functions Ψ(t) and φ(t) means (−1) n d n dt n Ψ(t) ≥ 0 , (−1) n d n dt n φ(t) ≥ 0 , n = 0, 1, 2, . . . , t ≥ 0 , (57) or equivalently, their representability as real Laplace transforms of non-negative generalized functions (or measures). For the generalizations of Equations (44)- (46), characteristic of the Poisson and Erlang distributions respectively, we must point out the Laplace transform pair with E  [19]. Then, by using the Laplace transform pairs (25) and Equations (52), (53), (58) in Equations (37) and (38), we have the generalized Poisson distribution, and the generalized Erlang pd f 's (of order k ≥ 1), The generalized Erlang distribution functions turn out to be

The Gnedenko-Kovalenko Theory of Thinning and the Mittag-Leffler Function
The thinning theory for a renewal process has been considered in detail by Gnedenko and Kovalenko [43] in the first edition of their book on Queue theory of 1968. However, the connection with the Laplace transform of the Mittag-Leffler function outlined at the end of this section in Equations (71) and (72), see also [44] and [45], is surprisingly not present in the second edition of the book by Gnedenko & Kovalenko in 1989.
We must note that other authors, like Szántai [46,47] speak of rarefaction in place of thinning. Let us sketch here the essentials of this theory: in the interest of transparency and easy readability we avoid the possible decoration of the relevant power law by multiplying it with a slowly varying function.
Denoting by t n , n = 1, 2, 3, . . . the time instants of events of a renewal process, assuming 0 = t 0 < t 1 < t 2 < t 3 < . . . , with i.i.d. waiting times T 1 = t 1 , T k = t k − t k−1 for k ≥ 2, (generically denoted by T), thinning (or rarefaction) means that for each positive index k a decision is made: the event happening in the instant t k is deleted with probability p or it is maintained with probability q = 1 − p, 0 < q < 1. This procedure produces a thinned or rarefied renewal process with fewer events (very few events if q is near zero, the case of particular interest) in a moderate span of time.
To compensate for this loss we change the unit of time so that we still have not very few but still a moderate number of events in a moderate span of time. Such change of the unit of time is equivalent to rescaling the waiting time, multiplying it with a positive factor τ so that we have waiting times τT 1 , τT 2 , τT 3 , . . . , and instants τt 1 , τt 2 , τt 3 , . . . , in the rescaled process. Our intention is, vaguely speaking, to dispose on τ in relation to the rarefaction parameter q in such a way that for q near zero in some sense the "average" number of events per unit of time remains unchanged. In an asymptotic sense we will make these considerations precise.
Denoting by F(t) = P(T ≤ t) the probability distribution function of the (original) waiting time T, by f (t) its density ( f (t) is a generalized function generating a probability measure) so that F(t) = t 0 f (t ) dt , and analogously by F k (t) and f k (t) the distribution and density, respectively, of the sum of k waiting times, we have recursively Observing that after a maintained event the next one of the original process is kept with probability q but dropped in favour of the second-next with probability p q and, generally, n − 1 events are dropped in favour of the n-th-next with probability p n−1 q, we get for the waiting time density of the thinned process the formula With the modified waiting time τ T we have hence the density f (t/τ)/τ, and analogously for the density of the sum of n waiting times f n (t/τ)/τ. The density of the waiting time of the rescaled (and thinned) process now turns out as In the Laplace domain we have f n (s) = f (s) n , hence (using p = 1 − q) from which by Laplace inversion we can, in principle, construct the waiting time density of the thinned process. By rescaling we get Being interested in stronger and stronger thinning (infinite thinning) let us now consider a scale of processes with the parameters τ (of rescaling) and q (of thinning), with q tending to zero under a scaling relation q = q(τ) yet to be specified.
We have essentially two cases for the waiting time distribution: its expectation value is finite or infinite. In the first case we put In the second case we assume a queue of power law type (dispensing with a possible decoration by a function slowly varying at infinity) Then, by the Karamata theory (see References [21,48]) the above conditions mean in the Laplace domain with a positive coefficient λ and 0 < β ≤ 1. The case β = 1 obviously corresponds to the situation with finite first moment (2.6a), whereas the case 0 < β < 1 is related to a power law queue with c = λ Γ(β + 1) sin(βπ)/π . Now, passing to the limit of q → 0 of infinite thinning under the scaling relation between the positive parameters q and τ, the Laplace transform of the rescaled density g q,τ (s) in (66) of the thinned process tends for fixed s to which corresponds to the Mittag-Leffler density Let us remark that Gnedenko and Kovalenko obtained (71) as the Laplace transform of the limiting density but did not identify it as the Laplace transform of a Mittag-Leffler type function. Observe that in the special case λ < ∞ we have β = 1, hence as the limiting process the Poisson process, as formerly shown in 1956 by Rényi [49].

The Continuous Time Random Walk (Ctrw) and the Mittag-Leffler Function
The name continuous time random walk (CTRW) became popular in physics after Montroll and Weiss (just to cite the pioneers) published a celebrated series of papers on random walks for modelling diffusion processes on lattices, see, for example, Reference [50], and the book by Weiss [51] with references therein. CTRWs are rather good and general phenomenological models for diffusion, including processes of anomalous transport, that can be understood in the framework of the classical renewal theory. In fact a CTRW can be considered as a compound renewal process (a simple renewal process with reward) or a random walk subordinated to a simple renewal process. Hereafter we will mainly follow the approach by Gorenflo & Mainardi, see, for example, Reference [52].
A spatially one-dimensional CTRW is generated by a sequence of independent identically distributed (iid) positive random waiting times T 1 , T 2 , T 3 , . . . , each having the same probability density function φ(t) , t > 0 , and a sequence of iid random jumps X 1 , X 2 , X 3 , . . . , in IR , each having the same probability density w(x) , x ∈ IR .
Let us remark that, for ease of language, we use the word density also for generalized functions in the sense of Gel'fand & Shilov [53], that can be interpreted as probability measures. Usually the probability density functions are abbreviated by pd f . We recall that φ(t) ≥ 0 with ∞ 0 φ(t) dt = 1 and w(x) ≥ 0 with +∞ −∞ w(x) dx = 1. Setting t 0 = 0 , t n = T 1 + T 2 + . . . T n for n ∈ IN , the wandering particle makes a jump of length X n in instant t n , so that its position is x 0 = 0 for 0 ≤ t < T 1 = t 1 , and x n = X 1 + X 2 + . . . X n , for t n ≤ t < t n+1 . We require the distribution of the waiting times and that of the jumps to be independent of each other. So, we have a compound renewal process (a renewal process with reward), compare Reference [42].
By natural probabilistic arguments we arrive at the integral equation for the probability density p(x, t) (a density with respect to the variable x) of the particle being in point x at instant t , in which δ(x) denotes the Dirac generalized function, and the survival function denotes the probability that at instant t the particle is still sitting in its starting position x = 0 . Clearly, Equation (73) satisfies the initial condition p(x, 0 + ) = δ(x).
Note that the special choice gives the pure renewal process, with position x(t) = N(t), denoting the counting function, and with jumps all of length 1 in positive direction happening at the renewal instants.
For many purposes the integral Equation (73) of CTRW can be easily treated by using the Laplace and Fourier transforms. Writing these as

then in the Laplace-Fourier domain Equation (73) reads
Introducing formally in the Laplace domain the auxiliary function and assuming that its Laplace inverse H(t) exists, we get, following Mainardi et al. [54], in the Laplace-Fourier domain the equation and in the space-time domain the generalized Kolmogorov-Feller equation with p(x, 0) = δ(x), where H(t) acts as a memory function.
If the Laplace inverse H(t) of the formally introduced function H(s) does not exist, we can formally set K(s) = 1/ H(s) and multiply (78) with K(s). Then, if K(t) exists, we get in place of (79) the alternative form of the generalized Kolmogorov-Feller equation with p(x, 0) = δ(x). where K(t) acts as a memory function Special choices of the memory function H(t) are (i) and (ii), see Equations (81) and (85): giving the exponential waiting time with In this case we obtain in the Fourier-Laplace domain and in the space-time domain the classical Kolmogorov-Feller equation giving the Mittag-Leffler waiting time with In this case we obtain in the Fourier-Laplace domain where * D β t denotes the fractional derivative of of order β in the Caputo sense, see Section 3. The time fractional Kolmogorov-Feller equation can be also expressed via the Riemann-Liouville fractional derivative D 1−β t , see again Section 3, that is with p(x, 0 + ) = δ(x). The equivalence of the two forms (88) and (89)  We note that the choice (i) may be considered as a limit of the choice (ii) as β = 1. In fact, in this limit we find H(s) ≡ 1 so H(t) = t −1 /Γ(0) ≡ δ(t) so that Equations (78)-(79) reduce to Equations (83)-(84), respectively. In this case the order of the Caputo derivative reduces to 1 and that of the R-L derivative to 0, whereas the Mittag-Leffler waiting time law reduces to the exponential. In the sequel we will formally unite the choices (i) and (ii) by defining what we call the Mittag-Leffler memory function whose Laplace transform is Thus we will consider the whole range 0 < β ≤ 1 by extending the Mittag-Leffler waiting time law in (86) to include the exponential law (82). (79) clearly may be supplemented by an arbitrary initial probability density p(x, 0) = f (x). The corresponding replacement of δ(x) by f (x) in (73) then requires in (76) multiplication of the term (1 − φ(s))/s by f (κ) and in (78) replacement of the LHS by H(s) s p(κ, s) − f (κ) . With p(x, 0) = δ(x) we obtain in p(x, t) the fundamental solution of Equation (79).

Note:
The probability density function for the waiting time distribution in terms of the Mittag-Leffler function was formerly given since 1995 by Hilfer [55][56][57]. In these papers the waiting time density was given with the Mittag-Leffler function in two parameters without noting the relation with the first derivative of the classical Mittag-Leffler function as stated in Equation (29). We also note that 10 years earlier Balakrishnan [58] had derived a similar expression without recognizing the Mittag-Leffler function. Like in the case of the thinning process dealt by Gnedenko-Kowalenko (see Section 7) once again the Mitag-Leffler function was unknown for the authors.

Manipulations: Rescaling and Respeeding
We now consider two types of manipulations on the CTRW by acting on its governing Equation (79)  (A) means change of the unit of time (measurement). We replace the random waiting time T by a waiting time τT, with the positive rescaling factor τ. Our idea is to take 0 < τ 1 in order to bring into near sight the distant future. In a moderate span of time we will so have a large number of jump events. For τ > 0 we get the rescaled waiting time density By decorating also the density p with an index τ we obtain the rescaled integral equation of the CTRW in the Laplace-Fourier domain as where, in analogy to (77), (B) means multiplying the quantity representing ∂ ∂t p(x, t) by a factor 1/a, where a > 0 is the respeeding factor: a > 1 means acceleration, 0 < a < 1 means deceleration. In the Laplace-Fourier representation this means multiplying the RHS of Equation (78) by the factor a since the expression s p(κ, s) − 1 corresponds to ∂ ∂t p(x, t).
We now chose to consider the procedures of rescaling and respeeding in their combination so that the equation in the transformed domain of the rescaled and respeeded process has the form Clearly, the two manipulations can be discussed separately: the choice {τ > 0, a = 1} means pure rescaling, the choice {τ = 1, a > 0} means pure respeeding of the original process. In the special case τ = 1 we only respeed the original system; if 0 < τ 1 we can counteract the compression effected by rescaling to again obtain a moderate number of events in a moderate span of time by respeeding (decelerating) with 0 < a 1. These vague notions will become clear as soon as we consider power law waiting times.
Defining now we finally get, in analogy to (78), the equation What is the combined effect of rescaling and respeeding on the waiting time density? In analogy to (77) and taking account of (96) we find φ τ,a (s) = 1 and so, for the deformation of the waiting time density, the essential formula Remark 2. The formula (99) has the same structure as the thinning formula (66) in Section 5 (just devoted to the thinning theory) by identification of a with q. In both problems we have a rescaled process defined by a time scale τ, and we send the relevant factors τ, a and q to zero under a proper relationship. However in the thinning theory the relevant independent parameter going to 0 is that of thinning (actually respeeding) whereas in the present problem it is the rescaling parameter τ.

Power Laws and Asymptotic Universality of the Mittag-Leffler Waiting Time Density
We have essentially two different situations for the waiting time distribution according to its first moment (the expectation value) being finite or infinite. In other words we assume for the waiting time pd f φ(t) either or For convenience we have dispensed in (101) with decorating by a slowly varying function at infinity the asymptotic power law. Then, by the standard Tauberian theory (see References [21,48]) the above conditions (100)-(101) mean in the Laplace domain the (comprehensive) asymptotic form where we have Then, fixing s as required by the continuity theorem of probability theory for Laplace transforms, taking and sending τ to zero, we obtain in the limit the Mittag-Leffler waiting time law. In fact, Equations (99) and (102) imply as τ → 0 with 0 < β ≤ 1, the Laplace transform of φ ML (t). This formula expresses the asymptotic universality of the Mittag-Leffler waiting time law that includes the exponential law for β = 1. It can easily be generalized to the case of power laws decorated with slowly varying functions, thereby using the Tauberian theory by Karamata (see again References [21,48]).

Comment:
The formula (105) says that our general power law waiting time density is gradually deformed into the Mittag-Leffler waiting time density as τ tends to zero.

Remark 3.
Let us stress here the distinguished character of the Mittag-Leffler waiting time density we can easily prove the identity Note that Equation (107) states the self-similarity of the combined operation rescaling-respeeding for the Mittag-Leffler waiting time density. In fact, (107) implies φ ML τ,a (t) = φ ML (t/c)/c with c = τ/a 1/β , which means replacing the random waiting time T ML by c T ML . As a consequences, choosing a = τ β we have Hence the Mittag-Leffler waiting time density is invariant against combined rescaling with τ and respeeding with a = τ β .
Observing (105) we can say that φ ML (t) is a τ → 0 attractor for any power law waiting time (101) under simultaneous rescaling with τ and respeeding with a = λτ β . In other words, this attraction property of the Mittag-Leffler probability distribution with respect to power law waiting times (with 0 < β ≤ 1) is a kind of analogy to the attraction of sums of power law jump distributions by stable distributions.

The Mittag-Leffler Functions W.R.T the Time Fractional Diffusion-Wave Equations and the Wright Functions
In this section we show the relations of the Mittag-Leffler function with the Wright function via Laplace and Fourier transformations, in order to provide other arguments to outline the role of the Mittag-Leffler in the Fractional Calculus. For this purpose, because of the necessity to work with two independent parameters we first recall the proper definitions of the Mittag-Leffler and the Wright function. Then we will consider the time fractional diffusion-wave equation with its fundamental solutions to the basic boundary value problem that urn out to be expressed in terms of some special cases of the Wright functions, the so called F and M functions. Finally we pay attention to some noteworthy formulas for the M-Wright function, including its connections with the Mittag-Leffler function.

Definitions and Main Properties of the Wright Functions
The classical Wright function, that we denote by W λ,µ (z), is defined by the series representation convergent in the whole complex plane, As a consequence W λ,µ (z) is an entire function for all λ ∈ (−1, +∞). Originally Wright assumed λ ≥ 0 in connection with his investigations on the asymptotic theory of partition [59,60] and only in 1940 he considered −1 < λ < 0, [61]. We note that in the Vol 3, Chapter 18 of the handbook of the Bateman Project [6], presumably for a misprint, the parameter λ is restricted to be non-negative, whereas the Wright functions remained practically ignored in other handbooks. In 1993 the present author, 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, as outlined in his papers published from 1994 to 1997, [62][63][64][65][66]. For other earlier treatments of this function we refer to the 1999 paper by Gorenflo, Luchko and Mainardi [67]).
In view of the asymptotic representation in the complex domain and of the Laplace transform the Wright functions were distinguished by the author in first kind (λ ≥ 0) and second kind (−1 < λ < 0) as outlined e.g., in the Appendix F of his book [22].
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 /Γ(µ).
Following the profs in the Appendix F in Reference [22] we get the following Laplace transform pairs of the Wright functions in terms of the Mittag-Leffler functions in two parameters, where r can be the time variable t > 0 or the space variable x > 0) for the first kind (λ ≥ 0) for the second kind (λ = −ν, 0 < ν < 1) The Wright functions of the first kind are useful to find the solutions of some (linear and non-linear) differential equations of fractional order as recently shown by Garra and Mainardi,[68].
Since the pioneering works in 1990's by the author, noteworthy cases of Wright functions of the second kind, known as auxiliary functions F and M play fundamental roles in solving the Signalling problem and the Cauchy value problem, respectively for the time fractional diffusion-wave equation.
We first recall hereafter these auxiliary functions in terms of the Wright functions of the second kind, following their power series representations. They read and The series representations of our auxiliary functions are derived from those of W λ,µ (z) in (109). We have: and where we have used the well-known reflection formula for the Gamma function,

The Time-Fractional Diffusion-Wave Equation and the Related Green Functions
For the reader's convenience let us recall the main formulas for the time fractional diffusion equations and their fundamental solutions (also referred to as the Green functions) for the Cauchy and Signalling problems. For more details we refer to References [69,70].
Denoting as usual x, t the space and time variables, and r = r(x, t) the response variable, the family of these evolution equations reads where the time derivative of order β is intended in the Caputo sense, namely is the operator * D β t , introduced in Section 3, but for order less than 1, see Equation (10), and a is a positive constant of dimension L 2 T −β . Thus we must distinguish the cases 0 < β ≤ 1 and 1 < β ≤ 2. We have It should be noted that in both cases 0 < β ≤ 1, 1 < β ≤ 2, the time fractional derivative in the L.H.S. of Equation (117) can be removed by a suitable fractional integration. leading to alternative forms where the necessary initial conditions at t = 0 + explicitly appear.
For this purpose we apply to Equation (117) the fractional integral operator of order β, namely For β ∈ (0, 1] we have: For β ∈ (1, 2] we have: Then, as a matter fact, we get the integro-differential equations: if 0 < β ≤ 1 : if 1 < β ≤ 2 : Denoting by f (x) , x ∈ IR and h(t) , t ∈ IR + sufficiently well-behaved functions, the basic boundary-value problems are thus formulated as following, assuming 0 < β ≤ 1, (a) Cauchy problem If 1 < β < 2 , we must add into (122) and (123) the initial values of the first time derivative of the field variable, r t (x, 0 + ) , since in this case the corresponding fractional derivative is expressed in terms of the second order time derivative. To ensure the continuous dependence of our solution with respect to the parameter β also in the transition from β = 1 − to β = 1 + , we agree to assume as it turns out from the integral forms (120)-(121).
In view of our subsequent analysis we find it convenient to set ν := β/2 , so 0 < ν ≤ 1/2 , ⇐⇒ 0 < β ≤ 1 , and from now on to add the parameter ν to the independent space-time variables x , t in the solutions, writing r = r(x, t; ν). For the Cauchy and Signalling problems we introduce the so-called Green functions G c (x, t; ν) and G s (x, t; ν), which represent the respective fundamental solutions, obtained when f (x) = δ(x) and h(t) = δ(t) . As a consequence, the solutions of the two basic problems are obtained by a space or time convolution according to It should be noted that in (126) G c (x, t; ν) = G c (|x|, t; ν) because the Green function of the Cauchy problem turns out to be an even function of x. According to a usual convention, in (127) the limits of integration are extended to take into account for the possibility of impulse functions centred at the extremes. Now we recall the results obtained in 1990's by the author that allow us to express the two Green functions in terms of the auxiliary functions F ν (ξ) and M ν (ξ) where, for x > 0, t > 0 acts as similarity variable. Then we obtain the Green functions in the space-time domain in the form We also recognize the following reciprocity relation for the original Green functions, Now F ν (ξ), M ν (ξ) are the auxiliary functions for the general case 0 < ν ≤ 1, which generalize those well known for the standard (Fourier) diffusion equation and for the standard (D'alembert) wave equation derived for ν = 1/2 and for ν = 1, respectively.

Some Noteworthy Results for the M ν Wright Function
In this survey we find worthwhile to concentrate our attention on a single auxiliary function, the M-Wright function, sometimes referred to as the Mainardi function. Indeed this function is indeed referred with this name in n the 1999 book by Podlubny [19], that is one of the most cited treatises on fractional calculus. Then this name is found in several successive papers and books related to fractional diffusion and wave processes, see for example, the relevant 2015 paper by Sandev et al. [71].
Let us now recall some interesting analytic results related to the so-called Mainardi function. One reason for the major attention is due to its straightforward generalization of the Gaussian probability density obtained for ν = 1/2, that is the fundamental solution of the Cauchy problem for the standard diffusion equation. Furthermore it allows an impressive visualization of the evolution with the order ν ∈ (0, 1) of the Green function of the Cauchy problem of the fractional diffusion wave Equation (129) as shown in the next figures with a = 1 and taking t = 1.
The readers are invited to look the YouTube video by my former student Armando Consiglio whose title is "Simulation of the M−Wright function", in which the author has shown the evolution of this function as the parameter ν changes between 0 and 0.85 in the interval (−5 < x < +5) of IR centered in x = 0 represented herewith in Figures 6 and 7 at fixed time t = 1.  The readers interested to have more details on the classical Wright functions would consult the recent survey by Luchko [72] and references therein.
In view of time-fractional diffusion processes related to time-fractional diffusion equations it is worthwhile to introduce the function in two variables which defines a spatial probability density in x evolving in time t with self-similarity exponent H = ν.
Of course for x ∈ IR we have to consider the symmetric version of the M-Wright function. obtained from (132) multiplying by 1/2 and replacing x by |x|.
Hereafter we provide a list of the main properties of this function, which can be derived from the Laplace and Fourier transforms for the corresponding Wright M-function in one variable presented in papers by Mainardi and recalled in the Appendix F of Reference [22].
For the Laplace transform of M ν (x, t) with respect to t > 0 and x > 0 we get respectively: For the Fourier transforms with respect to the spatial variable x we have for M ν (x, t) with x ∈ IR + , so that for the symmetric function M ν (|x|, t) we get Restricting our attention at the known analytic expressions of the M ν functions versus x at fixed time t = 1 we recall the following results for some special rational values of the parameter ν: ν = 1/3 (see Reference [22]) M 1/3 (x) = 3 2/3 Ai(x/3 1/3 ) , (137) ν = 1/2 (see Reference [22]) In the above equations Ai and Ai denote the Airy function and its first derivative.
Funding: This research received no external funding.
for me to devote my research interests to the application of fractional calculus in relaxation, oscillation phenomena governed by fractional ODEs and diffusion, wave phenomena governed by fractional PDEs. Once again I understood the relevance of the Mittag-Leffler functions but also that of the Wright functions, both of them classified as miscellaneous functions in the handbook of Bateeman project. I must note that, as far as I known, the Bateman handbook was the only one published in English to deal with these special function, and therefore accessible to me.
The year 1994 was the golden year for me as far as my acquaintance with fractional calculus and related special functions is concerned. Indeed I took profit by the acquaintance in three different conferences with the late Prof Gorenflo and Prof. Nigmatullin (in Bordeaux, France), with Prof. Podlubny and Prof. Caputo (in Astlanta, USA), and with Prof Virginia Kiryakova and the late Prof. Stankovic (in Sofia, Bulgaria), among other authorities of the fractional calculus. But it was with Prof Gorenflo that I started a collaboration for more than 20 years (1995-2015) motivated by our common interest towards the potential of the Mittag-Leffler functions in the applications of the fractional calculus.
Then, since 1997, I was interested in the emerging science of Econophysics thanks mainly to my younger colleague Enrico Scalas. With Gorenflo, Scalas and his student Raberto we published some papers on the advent of fractional Calculus in Econophysics, see e.g., [54] and my historical survey in Mathematics [77]. In 2007, on the occasion of the 80-birthday of Prof. Caputo, I published with Gorenflo a survey in Fractional Calculus and Applied Analysis [11] where I took the liberty to propose for the Mittag-Leffler function the (successful) title of the Queen Function of the Fractional Calculus. Some years earlier, Gorenflo had contacted the American Mathematical Society to give a specific number to the Mittag-Leffler function, that is 33E12, in the MSC classification. Gorenflo   Unfortunately, I lost Gorenflo's guidance and collaboration in 2015 when he suffered strong health troubles that led him to his death on 20 October 2017 at 87 years. He was Emeritus Professor of Mathematics at the Free University of Berlin since his retirement in 1998.