On some properties of the Mittag-Leffler function $E_\alpha(-t^\alpha)$, completely monotone for $t>0$ with $0<\alpha<1$

We analyse some peculiar properties of the function of the Mittag-Leffler (M-L) type, $e_\alpha(t):= E_\alpha(-t^\alpha)$ for $0<\alpha<1$ and $t>0$, which is known to be completely monotone (CM) with a non negative spectrum of frequencies and times, suitable to model fractional relaxation processes. We first note that these two spectra coincide so providing a universal scaling property of this function. Furthermore, we consider the problem of approximating our M-L function with simpler CM functions for small and large times. We provide two different sets of elementary CM functions that are asymptotically equivalent to $e_\alpha(t)$ as $t \to 0$ and $t \to \infty$.


Introduction
Since a few decades the special transcendental function known as Mittag-Leffler function has attracted an increasing attention of researchers because of its key role in treating problems related to integral and differential equations of fractional order.
Since its introduction by the Swedish mathematician Mittag-Leffler at the beginning of the last century up to the 1990's, this function was seldom considered by mathematicians and applied scientists.
Before the 1990's, from a mathematical point of view ,we recall the 1930 paper by Hille and Tamarkin [18] on the solutions of the Abel integral equation of the second kind, and the books by Davis [8], Sansone & Gerretsen [40], Dzherbashyan [10] (unfortunately in Russian), and finally Samko et al. [38]. Particular mention would be for the 1955 Handbook of High Transcendental Functions of the Bateman project [11], where this function was treated in the chapter devoted to miscellaneous functions. For former applications we recall an interesting note by Davis [8] reporting a previous research by Dr. Kenneth S. Cole in connection with nerve conduction, and the papers by Cole & Cole [7], Gross [16] and Caputo & Mainardi [5,6], where the Mittag-Leffler function was adopted to represent the responses in dielectric and viscoelastic media.
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, see e.g.the books by Kiryakova [23], Kilbas and Saigo [21], Marichev [30], Mathai & Saxena [32], Mathai et al. [33], Srivastava et.al. [43].
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 science. 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, see e.g. Beghin & Orsingher [3], Capelas et al. [4], Sandev et al. [39], Tomovski et al. [46]. Multi-index Mittag-Leffler functions have been introduced as well, see e.g. Kilbas et al. [19], Kilbas & Saigo [20], Kiryakova [24], Kiryakova & Luchko [25], but their extensive use has not yet been pointed out in applied sciences up to now.

The Mittag-Leffler function in fractional relaxation processes
The Mittag-Leffler function is defined by the following power series, convergent in the whole complex plane, We recognize that it is an entire function providing a simple generalization of the exponential function to which it reduces for α = 1. We also note that for the convergence of the power series in (2.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. For detailed asymptotic analysis, which includes the smooth transition across the Stokes lines, the interested reader is referred to Wong and Zhao [48].
In this paper we limit ourselves to 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 [37], see also Feller [12].
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 e.g Mainardi [28]. For mathematical details we refer the interested reader to the survey paper by Miller and Samko [34] and to the recent book by Schilling et al. [41].
In particular we are interested to the function that provides the solution to the fractional relaxation equation, see Gorenflo and Mainardi [15], Mainardi and Gorenflo [29], Mainardi [28].
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 whose solution is 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 e.g. Gorenflo and Mainardi [15], Podlubny [35], 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 (2.3), we solve the problems (2.5a) and (2.5b) with the Laplace transform technique, using the rules pertinent to the corresponding fractional derivatives. The problems (a) and (b) are equivalent since the Laplace transform of the solution in both cases comes out as that yields our function The Laplace transform pair can be proved by transforming term by term the power series representation of e α (t) in the R.H.S of (2.2). Furthermore, by anti-transforming the R.H.S of (2.11) 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 also Gorenflo and Mainardi [15], Since K α (r) is non-negative for all r in the integral, the above formula proves that e α (t) is 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 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 distribution 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 [28]. This kind of universal/scaling property seems a peculiar one for our Mittag-Leffler function e α (t). In Fig 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).

Asymptotic approximations to the Mittag-Lefler function
In Fig 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 two common asymptotic approximations
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.2). In fact, The long time approximation is derived from the asymptotic power series representation of e α (t) that turns out to be, see Erdélyi (1955), 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 (3.1). Of course, for sufficiently small and for sufficiently large values of t we have the inequality e 0 α (t) ≤ e ∞ α (t) , 0 < α < 1 .
pointing out a continuous line at an error 1% under which the approximations can be considered reliable.     We note from Figs 3-7 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 steep cusps occurring to the left of t = 10 0 in Figs 5, 6, and 7 indicates that there the relative error is falling down to zero.

The two rational asymptotic approximations
We now propose a new set of CM functions approximating e α (t): {f α (t), g α (t)}, alternative to {e 0 α (t), e ∞ α (t)}, obtained as the first Padè approximants [0/1] to the power series in t α (2.2) and (3.2), respectively. For more details on the theory of Padè Approximants we refer e.g. to Baker [1]. We thus obtain the following rational functions in t α : Now we prove the inequality as a straightforward consequence of the reflection formula of the gamma function. In fact, recalling the definitions (3.7)-(3.8) we have for ∀t ≥ 0 and α ∈ (0, 1): In Figs 8-12 LEFT we compare for α = 0.25, 0.5, 0.75, 0.90, 0.99 in logarithmic scales the function e α (t) (continuous line) and its rational asymptotic representations, f α (t) valid for t → 0 (dashed line) and g α (t) valid for t → ∞ (dotted line). We have chosen the time range 10 −5 ≤ t ≤ 10 +5 . In the RIGHT we have shown the plots of the relative errors (no longer in absolute values) pointing out a continuous line at an error 1% under which the approximation can be considered reliable.     Fig.9 Approximations f α (t) (dashed line) and g α (t) (dotted line) to e α (t) (LEFT) and the corresponding relative errors (RIGHT) in 10 −5 ≤ t ≤ 10 +5 for α = 0.50.    As a matter of fact, from the plots in Figs 8-12, we recognize, for the time range and for α ∈ (0, 1) considered by us, the relevant inequality that is g α (t) and f α (t) provide lower and upper bounds to our Mittag-Leffler function e α (t) . This of course can be seen as a conjecture that we leave as an open problem to be proved (or disproved) by specialists of CM functions.
We also see that whereas the short time approximation f α (t) turns out to be good only for small times, the long time approximation g α (t) is good (surprisingly) also for short times, falling down only in an intermediate time range. In fact from the RIGHT of Figs 8-12 we can estimate the ranges of validity when the relative error is less than 1%. We could further discuss on them.
Concerning Padè Approximants (PA) for the Mitag-Leffler functions E α (−x) for x > 0 and 0 < α < 1 we like to cite Freed et al. [13]for an appendix devoted to the table of PA, Starovotov and Starovotova [44] for theorems on the best uniform rational approximations. Quite recently Zeng & Chen [49] have found a global Padè approximation with degree 2 of the generalized Mittag-Leffler function E α,β (x) with x > 0, 0 < α < 1 and β ≥ α. In this interesting pre-print the uniform approximation can account for both the Taylor series for small arguments and asymptotic series for large arguments, but no upper and lower bounds are obtained as in our case. We point out the fact that our first PA's [0/1] are constructed from two different series (in positive and negative powers of t α ) both resulting CM functions whereas the successive PA's of higher order can exhibit oscillations for t > 0.

Conclusions
We have discussed some noteworthy properties of the Mittag-Leffler function E α (−t α ) with 0 < α < 1 in the range t > 0.
Being a completely monotone (CM) function, because of the Bermstein theorem this function admits non-negative frequency and time spectra. We have pointed out that these two spectra are equal, so providing a universal scaling property.
Furthermore, in view of a numerical approximation we have compared two different sets of approximating CM functions, asymptotically equivalent to E α (−t α ) for small and large times: the former commonly used in the literature, the latter probably new. The last set is constituted by two simple rational functions that provide upper and lower bounds, at least in our numerical examples on a large time range.
We are allowed to the conjecture that this bounding property is always valid for any t > 0 and for any α ∈ (0, 1): an open problem offered to specialists of CM functions.
Last but not the least, we note that our rational function approximating E α (−t α ) for large times provides (surprisingly) a good numerical approximation also for short times, so failing only in an intermediate range.

Final remark
Since the appearance of the first version (May 2013), this paper was submitted to the attention of some colleagues and to any interested reader of arXiv in order to get the proof of our conjecture. Only recently (October 2013) a proof was provided by Thomas Simon (University of Lille, France) based on probability arguments, see [42]. We note, however, that in July 2013 Renato Spigler (University of Rome 3, Italy) proved, in part, the conjecture made in eq. (3.11), for x := t α in some right neighborhood of x = 0. This set is the most critical in establishing such conjecture, see