\section{Introduction and statement of the main results}

\label{intro}

The classical Fisher-Tippett-Gnedenko theorem states that the limit distributions arising from $a_n (\max(X_1, \ldots, X_n) - b_n)$ with $a_n > 0, b_n \in{\mathbb{R}}$ and $(X_1, \ldots, X_n)$ an i.i.d real sample, can be classified up to positive affine transformation into three families:
$$\left\{\begin{array}{ll}
{\bf W}_\rho = {\bf L}^{1/\rho} & \quad \mbox{(Weibull distribution)}\\
{\bf F}_\rho = {\bf L}^{-1/\rho} & \quad \mbox{(Fr\'echet distribution)}\\
{\bf G} = \log {\bf L} & \quad \mbox{(Gumbel distribution)}
\end{array}\right.$$
where ${\bf L}$ is the unit exponential random variable, $\rho$ is a positive parameter and, with an abuse of notation which we will make throughout the paper, we have identified a random variable with its law. From the distribution function viewpoint, the three above extreme laws can also be obtained as the unique solution to a certain ordinary differential equation. More precisely, if $F(x)$ stands for a distribution function on ${\mathbb{R}}$ and ${\bar F} (x) = 1- F(x)$ denotes its associated survival function, the following equations
$$\left\{\begin{array}{lll}
F'(x) =  \rho x^{\rho -1} {\bar F}(x), & x >0, & F(0) = 0\\
{\bar F}'(x) = -\rho x^{-\rho -1} F(x), & x >0, & F(0) = 0\\
F'(x) = e^x {\bar F}(x), & x\in{\mathbb{R}} &
\end{array}\right.$$
have each a unique solution which is respectively given by the Weibull distribution ${\bf W}_\rho,$ the Fr\'echet distribution ${\bf F}_\rho$ and the Gumbel distribution ${\bf G}.$

Notice that those three equations involve a logarithmic derivative and that they are solved via the exponential function. In this paper, we will consider some extensions of these equations in the context of fractional calculus. Throughout, we shall refer to the Appendix for all definitions and notations on the fractional integrals and derivatives that we will consider. In fractional calculus, a fundamental role is played by the classical Mittag-Leffler function
$$E_\alpha(z) \; =\; \sum_{n\ge 0} \frac{z^n}{\Ga(1+n\alpha)}, \qquad \alpha > 0,\, z\in\CC,$$which can be viewed as a generalization of the exponential function. We refer to Chapter 3 in \cite{GKMR} for a modern account on this function, and also to Chapter 2 therein for an interesting historical overview.

Let us first discuss an example. It is well-known from the general results of Barrett \cite{B1954} - see also Lemma 3.24 and the inversion formula (E.1.10) in \cite{GKMR} - that for every $\alpha,\lambda > 0$ the function $x\mapsto E_\a(-\lambda x^\alpha)$ solves on ${\mathbb{R}}^+$ the following fractional differential equation
\begin{equation}
\label{Barr}
\Dap (1 -f)(x)\; =\;  \lambda f(x),
\end{equation}
where $\Dap$ is the progressive Liouville fractional derivative on the half-axis. Besides, it follows from the works of Pillai \cite{P1990} that for every $\alpha\in(0,1]$ the function $x\mapsto E_\a(-x^\alpha)$ is the survival function of a distribution on ${\mathbb{R}}^+.$ More precisely, one has
\begin{equation}
\label{MLU}
E_\a(-x^\alpha)\; =\; \esp\left[ e^{-x^\alpha {\bf Z}_\alpha^{-\alpha}}\right]\; =\; {\mathbb{P}}\left[ {\Z_\a} \times {\bf L}^{\frac{1}{\alpha}} > x\right]
\end{equation}
where ${\Z_\a}$ has a standard positive $\alpha-$stable distribution with the normalization $\esp[e^{-x{\Z_\a}}] = e^{-x^\alpha}$ and, here and throughout, the product is assumed to be independent. This shows that the distribution function of the random variable ${\Z_\a} \times {\bf L}^{\frac{1}{\alpha}}$ solves the fractional differential equation
\begin{equation}
\label{Barre}
\Dap F (x)\; =\;  {\bar F}(x)
\end{equation}
on $(0,\infty),$ with the initial condition $F(0) = 0.$

The above fact can be used to display another, dual example involving the regressive Liouville fractional derivative on the half-axis $\Dm$ and the generalized Mittag-Leffler function
$$E_{\alpha,\beta}(z) \; =\; \sum_{n\ge 0} \frac{z^n}{\Ga(\beta +n\alpha)}, \qquad \alpha,\beta > 0,\; z\in\CC.$$
On the one hand, it follows indeed from the formulas (4.10.13) and (E.2.6) in \cite{GKMR} that for every $\alpha,\lambda > 0$ the function $x\mapsto \Ga(\alpha)\, E_{\alpha,\alpha}(-\lambda x^{-\alpha})$ solves on $(0,\infty)$ the fractional differential equation
$$\Dm (1 -f)(x)\; =\;  \lambda x^{-2\alpha} f(x).$$
On the other hand, the above Pillai result shows that for every $\alpha\in(0,1],$ the function
$$\Ga(\alpha)\, E_{\alpha,\alpha}(-x^{-\alpha})\; =\; \Ga(1+\alpha) \,E_\a'(-x^{-\alpha})\; =\; {\mathbb{P}}\left[  ({\bf Z}_\alpha^{-1})^{(\alpha)}\times{\bf L}^{-\frac{1}{\alpha}}\le x\right]$$
is a distribution function on ${\mathbb{R}}^+,$ where the second equality follows from an elementary transformation of (\ref{MLU}) involving the size-bias $({\bf Z}_\alpha^{-1})^{(\alpha)}$ of order $\alpha$ of the inverse positive $\alpha-$stable random variable (recall that for $t\in{\mathbb{R}}$ and a positive random variable $X$ such that $\esp [X^t] < \infty$ the size-bias of order $t$ of $X$ is the random variable $X^{(t)}$ whose law is defined by
$$\esp [f(X^{(t)})]\; =\; \frac{\esp [ X^t f(X)]}{\esp[X^t]}$$
for all $f :{\mathbb{R}}^+\to{\mathbb{R}}$ bounded continuous, and that $\esp[{\bf Z}_\alpha^{-\alpha}] = 1/\Ga(1+\alpha)$). All of this shows that the distribution function of the random variable $({\bf Z}_\alpha^{-1})^{(\alpha)} \times {\bf L}^{-1/\alpha}$ solves the fractional differential equation
\begin{equation}
\label{Barres}
\Dm {\bar F} (x)\; =\;  x^{-2\alpha} F(x)
\end{equation}
on $(0,\infty),$ with the initial condition $F(0) = 0.$\\

In this paper, we wish to study more general fractional equations than (\ref{Barre}) and (\ref{Barres}), which are natural extensions of the above differential equations characterizing the classical extreme distributions.  Our findings involve the $\alpha-$stable subordinator $\{\sigma^{(\a)}_t\, t\ge 0\},$ which is the real L\'evy process starting from zero such that $\sigma^{(\a)}_1 \stackrel{d}{=} {\Z_\a}.$ For every $\alpha\in [0,1],$ its Laplace transform is given by
$$\esp[e^{-\lambda \sigma^{(\a)}_t}]\; =\; e^{-t\lambda^\alpha},\qquad \lambda, t\ge 0.$$
Observe that $\sigma^{(\a)}$ is a pure drift for $\alpha = 1$ that is $\sigma^{(1)}_t =t,$ and a pure killing at an exponential time ${\bf L}$ for $\alpha =0$ that is $\sigma^{(0)}_t = \infty$ on $\{t \ge {\bf L}\}$ and $\sigma^{(0)}_t = 0$ on $\{t < {\bf L}\}.$ Our first main result gives a fractional extension of the Weibull distribution.

\begin{theo}
\label{Weib}
For every $\lambda, \rho > 0$ and $\alpha\in [0,1],$ there exists a unique distribution function solving the fractional differential equation
\begin{equation}
\label{Weib+}
\Dap F (x)\; =\; \lambda\, x^{\rho -\alpha} {\bar F} (x)
\end{equation}
on $(0,\infty),$ with the initial condition $F(0) = 0.$ The corresponding random variable is
$${\W}_{\a,\lbd,\rho} \; \stackrel{d}{=}\; {\bf W}_\rho \; \times\; \left( \lambda \int_0^\infty \left( (1- \sigma^{(\a)}_t)_+^{}\right)^{\rho -\alpha} \, dt\right)^{-\frac{1}{\rho}}.$$
\end{theo}

In the above statement, we have used the standard notation $x_+ = \max (x,0)$ for $x\in{\mathbb{R}}.$ Observe that the integral on the right-hand side is finite a.s. for every $\rho > 0$ and $\alpha\in [0,1]:$ this is clear for $\alpha =1$ since $\rho > 0,$ and when $\alpha < 1$ this follows from the fact that $\sigma^{(\a)}$ is a non-decreasing c\`adl\`ag process which crosses the level 1 a.s. by a jump - see e.g. Theorem III.4 in \cite{B1996}. The above result shows that the fractional index $\alpha\in [0,1]$ of the derivative $\Dap$ gives rise to a non-trivial multiplicative perturbation of the Weibull random variable ${\bf W}_\rho$ given by the power of a certain Riemannian integral of the stable subordinator, whereas the parameter $\lambda$ is simply a scaling constant with ${\W}_{\a,\lbd,\rho} \stackrel{d}{=} \lambda^{-1/\rho}\, {\bf W}_{\alpha, 1,\rho}.$ One has also the identities\begin{equation}
\label{WeibIDs}
{\bf W}_{1, \rho,\rho}\; \stackrel{d}{=}\; {\bf W}_\rho\qquad\quad\mbox{and}\qquad\quad {\bf W}_{0, 1,\rho}\; \stackrel{d}{=}\; \frac{{\bf W}_\rho}{{\bf W}_\rho}
\end{equation}
where, here and throughout, the quotient is assumed to be independent. The random variable on the right has a Pareto distribution of type III - see \cite{ARY} for a study of the latter distribution, and the mapping in law $\alpha\mapsto {\bf W}_{\alpha,\rho^\alpha\!, \rho}$ can be viewed as a parametrized arc connecting this Pareto III distribution and the Weibull distribution, the parameter being the index of the underlying stable subordinator. \\

Our second main result gives a fractional extension of the Fr\'echet distribution.

\begin{theo}
\label{Fresh}
For every $\lambda, \rho > 0$ and $\alpha\in [0,1],$ there exists a unique distribution function solving the fractional differential equation
\begin{equation}
\label{Fresh+}
\Dm {\bar F}(x)\; =\; \lambda\, x^{-\rho -\alpha} F(x)
\end{equation}
on $(0,\infty),$ with the initial condition $F(0) = 0.$ The corresponding random variable is
$${\F}_{\a,\lambda,\rho} \; \stackrel{d}{=}\; {\bf F}_\rho \; \times\; \left( \lambda \int_0^\infty \left( 1+ \sigma^{(\a)}_t\right)^{-\rho -\alpha} \, dt\right)^{\frac{1}{\rho}}.$$
\end{theo}

In the above statement the integral on the right-hand side is finite a.s. by the law of the iterated logarithm at infinity for $\sigma^{(\a)}$ - see e.g. Theorem III.11 in \cite{B1996}. As above, the index $\alpha$ of the derivative $\Dm$ produces a multiplicative perturbation of the Fr\'echet random variable ${\bf F}_\rho$ via a Riemannian integral of $\sigma^{(\a)}$, whereas the parameter $\lambda$ is a scaling constant with ${\F}_{\a,\lambda,\rho} \stackrel{d}{=} \lambda^{1/\rho}\, {\bf F}_{\alpha, 1,\rho}.$ One has also the identities
\begin{equation}
\label{FreshIDs}
{\bf F}_{1, \rho,\rho}\; \stackrel{d}{=}\; {\bf F}_\rho\qquad\quad\mbox{and}\qquad\quad {\bf F}_{0, 1,\rho}\; \stackrel{d}{=}\; \frac{{\bf F}_\rho}{{\bf F}_\rho}\; \stackrel{d}{=}\; \frac{{\bf W}_\rho}{{\bf W}_\rho}\cdot
\end{equation}
It is interesting to observe from the two above theorems that the two mappings in law
$$\alpha\;\mapsto\; {\bf L}^{\pm\frac{1}{\rho}} \; \times\;\left( \rho^\alpha  \int_0^\infty \left( 1\mp \sigma^{(\a)}_t\right)^{\pm\rho -\alpha}_+ \, dt\right)^{\mp\frac{1}{\rho}},$$
with a sign switch at $\alpha = 0$ corresponding to the trivial equation $F = x^{\rho} {\bar F},$ produce a parametrized arc connecting the classical extreme random variables ${\bf W}_\rho$ and ${\bf F}_\rho.$ The traditional role of the $\alpha-$stable subordinator is to define a fractional Laplacian via the underlying subordinated semi-group whose marginals are the symmetric $\beta-$stable distributions with $\beta = 2\alpha$, the densities of the latter being up to some multiplicative constant the solutions to the fractional Cauchy problem
$$\frac{\partial f}{\partial t} \; =\; \frac{\partial^{\beta} f}{\partial x^{\beta}} $$
on $(0,\infty)\times{\mathbb{R}}^+$ - see \cite{GM1998} and the references therein for more on this standard subject. The above results show that the $\alpha-$stable subordinator is also involved, by means of its Riemannian integrals, in the solution to some fractional differential equations naturally associated to the Weibull or the Fr\'echet distribution. \\

The classical Gumbel distribution is the limit in law of either $\rho({\bf W}_\rho -1)\stackrel{d}{\longrightarrow}{\bf G}$ or $\rho(1-{\bf F}_\rho)\stackrel{d}{\longrightarrow}{\bf G}$ as $\rho\to\infty.$ In order to define a fractional Gumbel distribution, it is natural from the above to introduce the random variable
$${\bf G}_\alpha \; =\; \log\left(\int_0^\infty e^{-\sigma^{(\a)}_t}\, dt\right)$$
for every $\alpha\in [0,1],$ the a.s. convergence of the integral being a well-known consequence of
$$\esp\left[\int_0^\infty e^{-\sigma^{(\a)}_t}\, dt\right]\; =\;\int_0^\infty \esp\left[ e^{-\sigma^{(\a)}_t}\right]\, dt\; =\;\int_0^\infty e^{-t}\, dt\; =\;1. $$
We then have the following result involving the progressive Liouville fractional derivative on the line $\Dp$, which can be guessed at the formal limit $\rho\to\infty$ after a logarithmic change of variable in Theorem \ref{Weib} or Theorem \ref{Fresh}, and whose derivation is actually rigorous.

\begin{theo}
\label{Gumb}
For every $\lambda > 0$ and $\alpha\in [0,1],$ there exists a unique distribution function solving the fractional differential equation
\begin{equation}
\label{Gumb+}
\Dp F (x)\; =\; \lambda^\alpha e^{\lambda x}\, {\bar F}(x)
\end{equation}
on ${\mathbb{R}}.$ The corresponding random variable is
$$\Gal \; \stackrel{d}{=}\;\lambda^{-1}\left( {\bf G} \, -\, {\bf G}_{\alpha}\right).$$
\end{theo}

Unlike the fractional Weibull or Fr\'echet distributions, the perturbation on the standard Gumbel distribution induced by the parameter $\alpha$ of the progressive derivative $\Dp$ is linear and not multiplicative. Again, the parameter $\lambda$ is a scaling constant with $\Gal \stackrel{d}{=} \lambda^{-1} {\bf G}_{\alpha,1}.$ One has also the identities
$${\bf G}_{1, 1}\; \stackrel{d}{=}\; {\bf G}\qquad\quad\mbox{and}\qquad\quad {\bf G}_{0, 1}\; \stackrel{d}{=}\; {\bf G}\, -\,{\bf G}.$$
The random variable on the right has a standard logistic distribution - see \cite{St} for an account on the latter distribution, and the mapping in law $\alpha\mapsto {\bf G}_{\alpha,1}$ can be viewed as a parametrized arc connecting the logistic distribution and the Gumbel distribution. \\

The proof of the three above theorems is divided in two parts. In Section 2 we prove the uniqueness of the solutions to a more general class of equations. These results have an independent interest, because the uniqueness problem is not always addressed in the literature on fractional calculus. It is classical in analysis to show that the solution of a differential equation must be the fixed point of an integral equation, and we use the same method, in the framework of fractional calculus. In Section 3 we show the existence of a probability law solving the above equations, and we establish the explicit multiplicative, respectively additive, factorizations. This is done via a one-to-one correspondence with a Kilbas-Saigo function, respectively a Le Roy function, which leads to a family of positive random variables characterized by their entire moments and previously studied in \cite{LS}, in a more general context. \\

The Kilbas-Saigo functions are three-parameter generalizations of the classical Mittag-Leffler functions $E_\a$ and $E_{\alpha,\beta}$ defined by the convergent series representation
$$E_{\alpha,m,l} (z)\; =\; \sum_{n\ge 0} \left(\prod_{k=1}^n \frac{\Ga(1+ \alpha((k-1)m +l))}{\Ga(1 + \alpha((k-1)m + l +1))}\right) z^n, \qquad z\in\CC,$$
for $\alpha,m >0$ and $l > -1/\alpha,$ with the convention made here and throughout that an empty product always equals 1. Note that $E_\alpha = E_{\alpha,1,0}$ and that $\Ga(\beta) E_{\alpha,\beta} = E_{\alpha,1,\frac{\beta-1}{\alpha}}.$ We refer to Chapter 5.2 in \cite{GKMR} for an account, including an extension to complex values of the parameter $l.$ The Le Roy functions are simple generalizations of the exponential function defined by
$$\cL_\alpha(z)\; =\; \sum_{n\ge 0} \frac{z^n}{(n!)^\alpha}, \qquad z\in\CC,$$
for $\alpha > 0.$ Introduced in \cite{L1899} in the context of analytic continuation, these functions are much less studied than Mittag-Leffler's. See however the recent paper \cite{GMR} and the references therein. In this paper, we can hinge upon the fractional extreme distributions to deduce some analytical features of these two interesting classes of special functions, in analogy to some known properties of the classical Mittag-Leffler functions. More precisely, we characterize their complete monotonicity on the negative half-line, we prove certain monotonicity properties with respect to the parameters, we derive their exact asymptotic behaviour at $-\infty$, and we establish uniform and optimal hyperbolic bounds. In particular, we prove the complete monotonicity of the function $x\mapsto E_{\alpha,m,m-1}(-x)$ for every $\alpha\in (0,1]$ and $m >0,$ solving an open question stated in \cite{DMV}. In a less complete way, we also study the infinite divisibility of the fractional extreme distributions. All these analytical results are to be found in Section 4.

\section{Some uniqueness results on fractional hazard rates}

In this section we prove the uniqueness of distribution functions solving fractional equations of the type (\ref{Weib+}), (\ref{Fresh+}) or (\ref{Gumb+}), where the power function is replaced by a more general hazard rate.  We repeat that all definitions and notations on the fractional operators that we will consider here can be found in Appendix A.1.

\subsection{The Weibull case} We consider the equation
\begin{equation}
\label{pb1b+}
\Dap F\; =\; h {\bar F}, \qquad F(0) =0,
\end{equation}
where $F$ is a distribution function and $h : (0,\infty)\to{\mathbb{R}}^+$ is measurable and locally bounded. In the case $\alpha = 1,$ there exists a solution to (\ref{pb1b+}) if and only if
$$\int_0^\infty h(t)\, dt \; =\; \infty,$$
with a unique solution given by
$${\bar F}(x)\; =\; \exp\Big\{-\int_0^xh(t)dt\Big\}\cdot$$
Recall that the function $h$ is called either the reliability function or the hazard rate of the underlying positive random variable. In the case $\alpha =0,$ there exists a solution to (\ref{pb1b+}) if and only if $h$ is non-decreasing, $h(0) = 0$ and $h(x)\to\infty$ as $x\to\infty,$ with a unique solution given by
$${\bar F}(x)\; =\; \frac{1}{1+h(x)}\cdot$$
In order to state our result in the fractional case $\alpha\in (0,1),$ we introduce the following linear operator
$$A^{\alpha,h}_{0+} : f\;\mapsto\; {\rm I}^{\alpha}_{0+}(hf)$$
which is well-defined on measurable functions from $(0,\infty)$ to ${\mathbb{R}}^+.$

\begin{theo}\label{t1+}
Assume that there exists $\rho > 0$ such that $h(x) = O(x^{\rho-\alpha})$ as $x\to 0.$ Then, if it exists, the distribution function satisfying \eqref{pb1b+} is uniquely defined  by the convergent series
\begin{equation}
\label{rpf6+}
{\bar F} (x) \; =\; \sum_{n\ge 0} (-1)^n (A^{\alpha,h}_{0+})^n{\bf 1}\, (x), \qquad x > 0.
\end{equation}
\end{theo}

\proof

We begin by transforming \eqref{pb1b+} into an integral equation. Since $F$ is a distribution function on ${\mathbb{R}}^+,$ there exists a probability measure $\mu$ on ${\mathbb{R}}^+$ such that
\begin{equation}
\label{ItFr}
F(x)\; =\; \int_0^x \mu(dt)\; =\; {\rm I}^{\alpha}_{0+}(F_{\alpha,\mu})(x),\qquad x > 0,
\end{equation}
where we have set
$$F_{\alpha,\mu}(x) \; =\; \frac{1}{\Ga(1-\alpha)}\int_0^x (x-t)^{-\alpha} \,\mu(dt), \qquad x >0,$$
and the second equality in (\ref{ItFr}) is a direct consequence of Fubini's theorem and of the evaluation of the Beta integral of the first kind. Moreover, it is easy to see that the function $F_{\alpha,\mu},$ which might take infinite values, is nevertheless locally integrable at zero since $\mu$ is a probability. Hence, applying ${\rm I}^\alpha_{0+}$ on both sides of (\ref{pb1b+}), we can use the inversion formula (\ref{Fr0+}) and get
$$F\; =\; {\rm I}^{\alpha}_{0+}( h{\bar F})\; =\; A^{\alpha,h}_{0+}({\bar F})$$
on $(0,\infty).$ This leads to the fixed point equation ${\bar F} = {\bf 1} - A^{\alpha,h}_{0+}({\bar F})$ and, by the linearity of $A^{\alpha,h}_{0+}$, to
$${\bar F} (x) \; =\; \sum_{k = 0}^{n-1} (-1)^k (A^{\alpha,h}_{0+})^k{\bf 1}\, (x)\; +\; (-1)^{n} (A^{\alpha,h}_{0+})^{n}({\bar F})\, (x), \qquad x > 0,$$
for every $n\ge 1.$ Fixing now $x >0,$ the assumption made on $h$ implies that there exists a constant $c >0$ such that $h(t) \le c t^{\rho -\alpha}$ for every $t\in [0,x].$ Since moreover ${\bar F}(t)\le 1$ for every $t\in[0,x],$ an immediate induction based on the Beta integral of the first kind implies
$$0\;\le\;(A^{\alpha,h}_{0+})^{n}({\bar F})\, (x)\; \le\; \left(\prod_{k=1}^n \frac{\Ga(1-\alpha +k\rho)}{\Ga(1+k\rho)}\right) c^n x^{\rho n}\;\to\; 0\quad\mbox{as $n\to\infty,$}$$
where the convergence towards zero follows e.g. from (1.1.5) in \cite{AAR1999}. This completes the proof.

\endproof

\begin{rem}
\label{rt1+}
{\em (a) It can be proved that the series defined in \eqref{rpf6+} converges uniformly on every compact set and that it tends to 1 as $x\to 0.$ On the other hand, because of the alternate signs it is difficult to guess whether it is non-increasing and converges to zero as $x\to \infty.$ In the particular case when $h$ is a power function, Theorem \ref{Weib} gives a positive answer with the help of a Kilbas-Saigo function and the $\alpha-$stable subordinator. It would be interesting to know if there are other hazard rate functions $h$ such that the series in \eqref{rpf6+} is indeed a survival function on $(0,\infty).$

\medskip

(b)  In a different direction, sharing a certain analogy with the previous item,  the authors have introduced in \cite{TaVa2016, TaVa2017} generalized fractional distributions which are not conventional and classical distributions with fractional hazard rates. In \cite{TaVa2018}, the stochastic approximation of fractional probability distribution have been studied.

\medskip

(c) The above proof shows the more general fact that under the same assumption on $h$, for every $\alpha\in (0,1)$ there exists a unique bounded solution to
$$\Dap (f(0) -f)\; =\;  h f$$
on $(0,\infty),$ which is given by
$$f\; =\; f(0)\;\times\;\sum_{n\ge 0} (-1)^n (A^{\alpha,h}_{0+})^n{\bf 1}.$$
This can be viewed as an extension of (\ref{Barr}) which handles the case when $h$ is a positive constant.}

\end{rem}

\subsection{The Fr\'echet case} We consider the equation
\begin{equation}
\label{pb1b-}
\Dm {\bar F}\; =\; h F, \qquad F(0) =0,
\end{equation}
given on distribution functions, where $h : (0,\infty)\to{\mathbb{R}}^+$ is measurable and locally bounded. In the case $\alpha = 1,$ there exists a solution to (\ref{pb1b-}) if and only if $h$ is locally integrable on $(0, \infty]$ and such that
$$\int_0^\infty h(t)\, dt \; =\; \infty,$$
with a unique solution given by
$$F(x)\; =\; \exp\Big\{-\int_x^\infty h(t)dt\Big\}\cdot$$
 In the case $\alpha =0,$ there exists a solution to (\ref{pb1b-}) if and only if $h$ is non-increasing, $h(0+) =\infty$ and $h(x)\to 0$ as $x\to \infty,$ with a unique solution given by
$$F(x)\; =\; \frac{1}{1+h(x)}\cdot$$
In order to state our result in the fractional case $\alpha\in (0,1),$ we introduce the following linear operator
$$A^{\alpha,h}_{-} : f\;\mapsto\; {\rm I}^{\alpha}_{-}(hf),$$
which is well-defined on measurable functions from $(0,\infty)$ to ${\mathbb{R}}^+.$

\begin{theo}
\label{t1-}
Assume that there exists $\rho > 0$ such that $h(x) = O(x^{-\rho-\alpha})$ as $x\to \infty.$ Then, if it exists, the distribution function satisfying \eqref{pb1b-} is uniquely defined by the convergent series
\begin{equation}
\label{rpf6-}
F (x) \; =\; \sum_{n\ge 0} (-1)^n (A^{\alpha,h}_{-})^n{\bf 1}\, (x), \qquad x > 0.
\end{equation}
\end{theo}

\proof

The proof is analogous to that of \eqref{t1-}, except that we deal with survival functions. Since ${\bar F}$ is a survival function on ${\mathbb{R}}^+,$ there exists a probability measure $\mu$ on ${\mathbb{R}}^+$ such that
$${\bar F}(x)\; =\; \int_x^\infty \mu(dt)\; =\; {\rm I}^{\alpha}_{-}({\bar F}_{\alpha,\mu})(x),\qquad x > 0,$$
where we have set
$${\bar F}_{\alpha,\mu}(x) \; =\; \frac{1}{\Ga(1-\alpha)}\int_x^\infty (t-x)^{-\alpha} \,\mu(dt), \qquad x >0,$$
and used Fubini's theorem together with the evaluation of the Beta integral of the second kind. The function ${\bar F}_{\alpha,\mu}$ is locally integrable at infinity since for every $y >0$ one has
$$\int_y^\infty {\bar F}_{\alpha,\mu}(x) \, dx \; =\; \frac{1}{\Ga(2-\alpha)}\int_y^\infty (t-y)^{1-\alpha} \,\mu(dt) \; =\; \frac{1}{\Ga(2-\alpha)}\int_y^\infty (s-y)^{-\alpha} {\bar F}(s) \,ds\; <\; \infty,$$
the finiteness of the third integral following from the fact that $\Dm{\bar F}$ must be finite on $(0,\infty).$ Hence, we can apply the inversion formula (\ref{Fr0-}) and get
${\bar F} = {\rm I}^{\alpha}_{-}( h F) = A^{\alpha,h}_{-}(F)$ on $(0,\infty),$ which leads to $F = {\bf 1} - A^{\alpha,h}_{-}(F)$ and then to
$$F (x) \; =\; \sum_{k = 0}^{n-1} (-1)^k (A^{\alpha,h}_{-})^k{\bf 1}\, (x)\; +\; (-1)^{n} (A^{\alpha,h}_{-})^{n}(F)\, (x), \qquad x > 0,$$
for every $n\ge 1.$ Fixing now $x >0,$ the assumption made on $h$ implies that there exists a constant $c >0$ such that $h(t) \le c t^{-\rho -\alpha}$ for every $t > x.$ Since moreover $F(t)\le 1$ for every $t >x,$ an induction based on the Beta integral of the second kind implies
$$0\;\le\;(A^{\alpha,h}_{-})^{n}(F)\, (x)\; \le\; \left(\prod_{k=1}^n \frac{\Ga(k\rho)}{\Ga(\alpha+k\rho)}\right) c^n x^{-\rho n}\;\to\; 0\quad\mbox{as $n\to\infty,$}$$
which completes the proof.

\endproof

\begin{rem}
\label{rt1-}
{\em (a) It can be proved that the series defined in \eqref{rpf6-} converges uniformly on every compact set of $(0,\infty)$ and that it tends to 1 as $x\to \infty.$ As above, it is not easy to guess from the alternate signs whether the series is non-decreasing and converges to zero as $x\to 0.$ The case when $h$ is a power function gives a positive answer in Theorem \ref{Fresh}, with the help of another Kilbas-Saigo function. It would be interesting to know if there are other functions $h$ such that the series in \eqref{rpf6-} is indeed a distribution function on $(0,\infty).$

\medskip

(b) The above proof shows that under the same assumption on $h$, for every $\alpha\in (0,1)$ there exists a unique bounded function having a limit $\ell$ at infinity and solving
$$\Dm (\ell -f)\; =\;  h f$$
on $(0,\infty),$ which is given by
$$f\; =\; \ell \;\times\;\sum_{n\ge 0} (-1)^n (A^{\alpha,h}_{-})^n\,{\bf 1}.$$
Observe that this solution is zero if $\ell =0.$}
\end{rem}

\subsection{The Gumbel case} We consider the equation
\begin{equation}
\label{pb1b}
\Dp F\; =\; h {\bar F}, \qquad F(x) > 0\;\; \mbox{on ${\mathbb{R}}$}
\end{equation}
given on distribution functions, where $h : {\mathbb{R}}\mapsto{\mathbb{R}}^+$ is measurable and locally bounded. In the case $\alpha = 1,$ there exists a solution to (\ref{pb1b}) if and only if $h$ is locally integrable at $-\infty$ and such that
$$\int_{\mathbb{R}} h(t)\, dt \; =\; \infty,$$
with a unique solution given by
$${\bar F}(x)\; =\; \exp\Big\{-\int_{-\infty}^xh(t)dt\Big\}\cdot$$
In the case $\alpha =0,$ there exists a solution to (\ref{pb1b+}) if and only if $h$ is non-decreasing, $h(x) \to 0$ as $x\to -\infty$ and $h(x)\to\infty$ as $x\to\infty,$ with a unique solution given by
$${\bar F}(x)\; =\; \frac{1}{1+h(x)}\cdot$$
In order to state our result in the fractional case $\alpha\in (0,1),$ we introduce the following linear operator
$$A^{\alpha,h}_{+} : f\;\mapsto\; {\rm I}^{\alpha}_{+}(hf),$$
which is well-defined on measurable functions from ${\mathbb{R}}$ to ${\mathbb{R}}^+.$ The following result is a simple variation on Theorem \ref{t1-}.

\begin{theo}
\label{t1G}
Assume that there exists $\lambda > 0$ such that $h(-x) = O(e^{-\lambda x})$ as $x\to \infty.$ Then, if it exists, the distribution function satisfying \eqref{pb1b} is uniquely defined by the convergent series
\begin{equation}
\label{rpf6}
{\bar F} (x) \; =\; \sum_{n\ge 0} (-1)^n (A^{\alpha,h}_{+})^n{\bf 1}\, (x), \qquad x\in{\mathbb{R}}.
\end{equation}
\end{theo}

\proof
Setting $G (x) = {\bar F}(-x),$ we see that $G$ is a distribution function on ${\mathbb{R}}$ with ${\bar G}(x) = F(-x).$ Changing the variable transforms \eqref{pb1b} into
$$\Dm {\bar G}\; =\; g G, \qquad G(x) > 0\;\; \mbox{on ${\mathbb{R}}$},$$
with $g(x) = h(-x).$ An argument similar to the proof of Theorem \ref{t1-}, where the evaluation of the Beta integral of the second kind is replaced by that of the Gamma integral, gives then the unique solution
$$G (x) \; =\; \sum_{n\ge 0} (-1)^n (A^{\alpha,g}_{-})^n{\bf 1}\, (x), \qquad x \in{\mathbb{R}}.$$
Changing the variable backwards, we obtain the required \eqref{rpf6}.

\endproof

\section{Proof of the main theorems}

In this section we show the existence of the real random variables associated to the fractional differential equations (\ref{Weib+}), (\ref{Fresh+}) and (\ref{Gumb+}), and we express them in terms of the unit exponential random variable and an integral transform of an independent $\alpha-$stable subordinator. The main ingredient in the proof is the following infinite independent product
$${\bf T}(a,b,c)\; =\; \prod_{n\ge 0} \left(\frac{a+nb +c}{a+nb}\right) {\bf B}_{a+ nb, c}$$
where, here and throughout, ${\bf B}_{a,b}$ denotes a standard Beta random variable. We refer to Section 2.1 in \cite{LS} for more details on this product, including the fact that it is a.s. convergent for every $a,b,c > 0.$ We also mention from Proposition 2 in \cite{LS} that its Mellin transform is
$$\esp[{\bf T}(a,b,c)^s]\;=\; \left( \frac{\Ga(ab^{-1})}{\Ga ((a+c)b^{-1})}\right)^s\times\; \frac{[a+c; b]_s}{[a; b]_s}$$
for every $s > -a,$ where $[z;\delta]_s$ stands for the generalized Pochhammer symbol which is defined in (\ref{Poch}) below. In general, we shall refer to the Appendix A.2 for all the properties of Barnes' double Gamma function and its associated Pochhammer symbol that we will need.

\subsection{Proof of Theorem \ref{Weib}}

We first consider the case $\alpha\in (0,1).$ The uniqueness is a direct consequence of Theorem \ref{t1+} with $h(x) = \lambda x^{\rho -\alpha}.$ Moreover, we know by (\ref{rpf6+}) that a distribution solving (\ref{Weib+}), if it exists, has survival function
$${\bar F} (x) \; =\; \sum_{n\ge 0} (-1)^n (A^{\alpha,h}_{0+})^n{\bf 1}\, (x), \qquad x > 0.$$
Since
$$(\Iap t^\beta)(x) \; =\; \frac{\Ga(\beta + 1)}{\Ga(\alpha +\beta +1)} \, x^{\alpha +\beta +1}$$
for every $\beta > -1,$ an induction implies
$${\bar F}(x) \; = \; \sum_{n\ge 0}\left(\prod_{k=1}^{n}\frac{\Gamma(k\rho+1-\alpha)}{\Gamma(k\rho+1)}\right)(-\lambda  x^{\rho})^n\; = \; E_{\alpha,\frac{\rho}{\alpha},\frac{\rho}{\alpha}-1}(-\lambda x^\rho)$$
for every $x\ge 0.$ Observe that alternatively, the fact that $E_{\alpha,\frac{\rho}{\alpha},\frac{\rho}{\alpha}-1}(-\lambda x^\rho)$ is a solution to (\ref{Weib+}) follows from Theorem 5.27 and the inversion formula (E.1.10) in \cite{GKMR}.\\

It thus remains to prove that $x\mapsto E_{\alpha,\frac{\rho}{\alpha},\frac{\rho}{\alpha}-1}(-\lambda x^\rho)$ is a survival function on ${\mathbb{R}}^+$ and to identify the underlying positive random variable. For every $z\in\CC,$ one has
$$E_{\alpha,\frac{\rho}{\alpha},\frac{\rho}{\alpha}-1}(z)\; =\; \sum_{n\ge 0} a_n(\alpha,\rho) \,\frac{z^n}{n!}$$
with
$$a_n(\alpha,\rho)\;= \;\rho^{-n} \,\prod_{k=1}^n\frac{\Gamma(k\rho +1 -\alpha)}{\Gamma(k\rho)}\cdot$$
Let us now consider the positive random variable
$$\frac{\Ga(\rho +1-\alpha)}{\Ga(\rho +1)}\; {\bf T}(1,\rho^{-1}, (1-\alpha)\rho^{-1}).$$
By the aforementioned Proposition 2 in \cite{LS} and the concatenation formula (\ref{Conca1}), the positive entire moments of the latter random variable are given by
$$\delta^{n} \left( \frac{[1+(1-\alpha)\delta;\delta]_n}{[1;\delta]_n}\right)\; =\; a_n(\alpha,\rho)$$
where $\delta = \rho^{-1}.$ The Stirling formula (\ref{Stirling}) implies
$$a_n(\alpha,\rho)^{-\frac{1}{2n}}\; \sim\; \kappa n^{\frac{\alpha-1}{2}}\qquad
\mbox{as $n\to \infty$}$$
for some positive constant $\kappa,$ so that Carleman's criterion
$$\sum_{n\ge 1} a_n(\alpha,\rho)^{-\frac{1}{2n}}\; =\; \infty$$
is fulfilled, and the law of the latter random variable is determined by its positive entire moments. Finally, it follows from Theorem (b) (i) in \cite{LS} with $q =\rho -\alpha$ that
$$\frac{\Ga(\rho +1-\alpha)}{\Ga(\rho +1)}\; {\bf T}(1,\rho^{-1}, (1-\alpha)\rho^{-1})\; \stackrel{d}{=}\; \cA(\alpha, 0,\rho -\alpha)\; \stackrel{d}{=}\; \int_0^\infty \left( 1- \sigma^{(\a)}_t\right)_+^{\rho -\alpha} \, dt.$$
All in all, we have shown that
\begin{equation}
\label{KS1}
E_{\alpha,\frac{\rho}{\alpha},\frac{\rho}{\alpha}-1}(z)\; =\; \esp\left[\exp \left\{ z\, \int_0^\infty \left( 1- \sigma^{(\a)}_t\right)_+^{\rho -\alpha} \, dt\right\}\right], \qquad z\in\CC.
\end{equation}
This implies
\begin{eqnarray*}
E_{\alpha,\frac{\rho}{\alpha},\frac{\rho}{\alpha}-1}(-\lambda x^\rho) & = & \esp\left[\exp\left\{ -\lambda x^{\rho}\, \int_0^\infty \left( 1- \sigma^{(\a)}_t\right)_+^{\rho -\alpha} \, dt\right\}\right]\\
& = & {\mathbb{P}}\left[\, {\bf L}\, >\, x^{\rho}\,\left( \lambda\int_0^\infty \left( 1- \sigma^{(\a)}_t\right)_+^{\rho -\alpha} \, dt\right)\right]\; = \; {\mathbb{P}}\left[ {\W}_{\a,\lbd,\rho} > x\right]
\end{eqnarray*}
for every $x\ge 0,$ which completes the proof for $\alpha\in(0,1).$ The case $\alpha =1$ is that of the classical Weibull distribution and was already discussed in the introduction. Finally, the case $\alpha = 0$ amounts to solving $F = \lambda x^\rho (1-F),$ which yields ${\bar F}(x) = 1/(1+ \lambda x^\rho)$ and
$${\bf W}_{0,\lambda,\rho}\;\stackrel{d}{=}\; \left( \frac{{\bf L}}{\lambda{\bf L}}\right)^{\frac{1}{\rho}}\;\stackrel{d}{=}\; {\bf W}_\rho\,\times\,\left( \lambda \int_0^\infty \left( 1- \sigma^{(0)}_t\right)_+^\rho \, dt\right)^{-\frac{1}{\rho}}.$$
\hfill$\square$

\begin{rem}
\label{rW1}
{\em (a) The main result of \cite{JSW2018} implies the identification
$$\left(\frac{\Ga(\rho +1-\alpha)}{\Ga(\rho)}\; {\bf T}(1,\rho^{-1}, (1-\alpha)\rho^{-1})\right)^{-\frac{1}{\rho}}\stackrel{d}{=}\;\,\cG(\rho +1-\alpha, 1-\alpha)$$
where $\cG(m,a)$ is the generalized stable random variable with parameters $m > a > 0,$ whose density is up to normalization the unique positive solution to
\begin{equation}
\label{JSW1}
 {\rm I}^a_{0+} f \; =\; x^m f
\end{equation}
on $(0,\infty).$ This yields the further identity in law
\begin{equation}
\label{JSW2}
{\W}_{\a,\lbd,\rho}\; \stackrel{d}{=}\; \left(\rho\lambda^{-1}\right)^{\frac{1}{\rho}}{\bf W}_\rho\,\times\, \cG(\rho +1-\alpha, 1-\alpha).
\end{equation}
In particular - see the introduction in \cite{JSW2018} for the third identity, one has
$${\bf W}_{\alpha,1,\alpha}\; \stackrel{d}{=}\; \alpha^{1/\alpha}{\bf W}_{\alpha,\alpha,\alpha}\; \stackrel{d}{=}\;  (\alpha^{1/\alpha} \cG(1, 1-\alpha))\,\times\,{\bf W}_\alpha\; \stackrel{d}{=}\; {\Z_\a}\,\times\,{\bf L}^{1/\alpha},$$
in accordance with the Pillai distribution mentioned in the introduction since
$${\mathbb{P}}\left[ {\bf W}_{\alpha,1,\alpha} > x\right]\; =\; E_{\alpha,1,0}(-x^\alpha)\; =\; E_{\alpha}(-x^\alpha), \qquad x\ge 0.$$
Notice that the integro-differential equation (\ref{JSW1}) shares some formal similarities with (\ref{Weib+}). Nevertheless it is essentially different because it deals with densities whereas (\ref{Weib+}) deals with distribution functions.

\medskip

(b) There exist unique solutions to fractional differential equations of the type (\ref{Weib+}) without the restriction to distribution functions. The main result of \cite{BDP2008} states that for every $\alpha\in (0,1],$ there is a unique solution to
$$(\alpha+1)\Dap f\; = \; -x f$$
satisfying the boundary condition
$$f(x)\; \sim\; \frac{x^{\alpha-1}}{\Ga(\alpha)\Ga(1/(\alpha+1))}$$
at zero, which is the density function of the running maximum of a spectrally positive $(\alpha+1)-$stable L\'evy process starting from zero.

\medskip

(c) With the above notation, one has
$$a_n(\alpha,\rho)\; =\; \frac{n!}{\Phi_{\alpha,\rho}(1)\cdots\Phi_{\alpha,\rho}(n)}$$
with
\begin{eqnarray*}
\Phi_{\alpha,\rho}(x) \; = \; \frac{\Ga(1+\rho x)}{\Ga(1-\alpha +\rho x)}& = & \frac{1}{\Ga(1-\alpha)}\; +\; \frac{\alpha}{\rho\Ga(1-\alpha)}\int_0^\infty (1- e^{-x u}) \frac{e^{-u\rho^{-1}}}{(1-e^{-u\rho^{-1}})^{\alpha +1}}\, du
\end{eqnarray*}
a Bernstein function. By Proposition 3.3 in \cite{CPY}, we deduce
\begin{equation}
\label{Txi}
\frac{\Ga(\rho +1-\alpha)}{\Ga(\rho+1)}\; {\bf T}(1,\rho^{-1}, (1-\alpha)\rho^{-1})\;\stackrel{d}{=}\;\int_0^\infty e^{-\xi^{(\alpha,\rho)}_t}\, dt
\end{equation}
where $\left\{ \xi^{(\alpha,\rho)}_t\! , \; t\ge 0\right\}$ is the subordinator having Laplace exponent
$$\esp\left[ e^{-\lambda \xi^{(\alpha,\rho)}_t}\right]\; = \; e^{-t \Phi_{\alpha,\rho}(\lambda)}.$$
This leads to
$${\W}_{\a,\lbd,\rho}\; \stackrel{d}{=}\; {\bf W}_\rho\;\times\; \left( \lambda \int_0^\infty e^{-\xi^{(\alpha,\rho)}_t}\, dt \right)^{-\frac{1}{\rho}}.$$
Let us notice that the identification (\ref{Txi}) can also be deduced from Corollary 5 in \cite{LS} is the case $q = \rho -\alpha > -\alpha$ and ${\hat \rho} =1,$ with the notation therein. Observe finally that this is consistent with the limiting case $\alpha = 1$ with $\xi^{(1,\rho)}_t = \rho\, t$ and $\alpha = 0$ with $\xi^{(0,\rho)}_t = \sigma^{(0)}_t.$

\medskip

(d) The above proof shows that the function
$$x\;\mapsto\;E_{\alpha,m,m-1}(-x)$$
is completely monotone (CM) for every $\alpha\in (0,1]$ and $m >0.$ This can be viewed as a generalization of the classic result by Pollard that $E_\alpha(-x) = E_{\alpha,1,0}(-x)$ is CM for every $\alpha\in (0,1]$ - see e.g. Proposition 3.23 in \cite{GKMR}. As already mentioned in the introduction, this also solves a conjecture stated in \cite{DMV} - see Section 4 and equations (10) and (11) therein. Notice that the formula (\ref{KS1}) implies the Bernstein representation
$$E_{\alpha,m,m-1}(-x)\; =\; \esp\left[ \exp -x \left\{ \int_0^\infty \left( 1- \sigma^{(\a)}_t\right)_+^{\alpha (m-1)} \, dt\right\}\right], \qquad x\ge 0.$$
For $m = 1,$ setting $T_\alpha = \inf\{ t > 0,\; \sigma^{(\a)}_t > 1\}\stackrel{d}{=} {\bf Z}_\alpha^{-\alpha}$ we obtain
$$E_\alpha(-x)\; =\; E_{\alpha,1,0}(-x)\; =\; \esp [e^{-x T_\alpha}]\; =\; \esp [e^{-x {\bf Z}_\alpha^{-\alpha}}],$$
a well-known fact following from our discussion prior to the statement of Theorem \ref{Weib}. See \cite{TS2015} for other CM functions related to $E_\alpha.$ In Section \ref{FurthA} below, we will generalize this fact and show some further analytical properties of the Kilbas-Saigo function $E_{\alpha,m,m-1}.$}

\end{rem}

We end this section with a convergent series representation, in the non-explicit case $\alpha\in(0,1),$ for the density $f^{\bf W}_{\alpha,\lambda,\rho}$ of ${\W}_{\a,\lbd,\rho}.$ This is an immediate consequence of a term-by-term differentiation of the survival function
$${\mathbb{P}}[{\W}_{\a,\lbd,\rho} > x]\; =\; E_{\alpha,\frac{\rho}{\alpha},\frac{\rho}{\alpha} -1}(-\lambda x^\rho)$$
which was obtained during the proof of Theorem \ref{Weib}.

\begin{coro}
\label{WeibD}
For every $\alpha\in (0,1),$ the density of ${\W}_{\a,\lbd,\rho}$ has the following convergent series representation on $(0,\infty)$
$$f^{\bf W}_{\alpha,\lambda,\rho}(x)\; =\; \lambda \, x^{\rho -1}\,\sum_{n\ge 0}\left( \prod_{j=0}^n \frac{\Ga(j\rho + \rho + 1 -\alpha)}{\Ga(j\rho +\rho)}\right) \frac{(-\lambda x^\rho)^n}{\rho^n\,n!}\cdot$$
\end{coro}

\subsection{Proof of Theorem \ref{Fresh}}

The proof is analogous to that of Theorem \ref{Weib}, except that we will deal with distribution functions instead of survival functions. We begin with the case $\alpha\in (0,1)$ and the uniqueness follows from Theorem \ref{t1-}. Besides, we know by (\ref{rpf6-}) that a distribution solving (\ref{Fresh+}), if it exists, must have distribution function
$$F (x) \; =\; \sum_{n\ge 0} (-1)^n (A^{\alpha,h}_{-})^n{\bf 1}\, (x), \qquad x > 0,$$
with $h(x) = \lambda x^{-\rho -\alpha}.$  Since
$$(\Im t^{-\beta})(x) \; =\; \frac{\Ga(\beta -\alpha)}{\Ga(\beta)} \, x^{\alpha -\beta}$$
for every $\beta > \alpha,$ an induction implies
$$F(x) \; = \; \sum_{n\ge 0}\left(\prod_{k=1}^{n}\frac{\Gamma(k\rho)}{\Gamma(k\rho+\alpha)}\right)(-\lambda  x^{-\rho})^n\; = \; E_{\alpha,\frac{\rho}{\alpha},\frac{\rho-1}{\alpha}}(-\lambda x^{-\rho})$$
for every $x > 0.$ Alternatively, the fact that $E_{\alpha,\frac{\rho}{\alpha},\frac{\rho-1}{\alpha}}(-\lambda x^{-\rho})$ is a solution to (\ref{Fresh+}) follows from Theorem 5.30 and the inversion formula (E.1.10) in \cite{GKMR}. It remains to prove that $x\mapsto E_{\alpha,\frac{\rho}{\alpha},\frac{\rho-1}{\alpha}}(-\lambda x^{-\rho})$ is a distribution function on $(0,\infty)$ and to identify the underlying positive random variable. For every $z\in\CC,$ one has
$$E_{\alpha,\frac{\rho}{\alpha},\frac{\rho-1}{\alpha}}(z)\; =\; \sum_{n\ge 0} b_n(\alpha,\rho) \,\frac{z^n}{n!}$$
with
$$b_n(\alpha,\rho)\;= \;\rho^{-n} \,\prod_{k=1}^n\frac{\Gamma(k\rho +1)}{\Gamma(k\rho+\alpha)}\cdot$$
Reasoning exactly as above implies that $\{b_n(\alpha,\rho), \; n\ge 0\}$ is the determinate integer moment sequence of the positive random variable
$$\frac{\Ga(\rho)}{\Ga(\rho +\alpha)}\; {\bf T}(1 +\alpha\rho^{-1},\rho^{-1}, (1-\alpha)\rho^{-1})\; \stackrel{d}{=}\;\int_0^\infty \left( 1+ \sigma^{(\a)}_t\right)^{-\rho -\alpha} \, dt,$$
where the identity in law follows from Corollary 3 in \cite{LS}. We have hence shown that
\begin{equation}
\label{KS2}
E_{\alpha,\frac{\rho}{\alpha},\frac{\rho-1}{\alpha}}(z)\; =\; \esp\left[\exp \left\{ z\, \int_0^\infty \left( 1+ \sigma^{(\a)}_t\right)^{-\rho -\alpha} \, dt\right\}\right], \qquad z\in\CC.
\end{equation}
This implies
\begin{eqnarray*}
E_{\alpha,\frac{\rho}{\alpha},\frac{\rho-1}{\alpha}}(-\lambda x^{-\rho}) & = & \esp\left[\exp\left\{ -\lambda x^{-\rho}\, \int_0^\infty \left( 1+ \sigma^{(\a)}_t\right)^{-\rho -\alpha} \, dt\right\}\right]\\
& = & {\mathbb{P}}\left[\, {\bf L}\, \ge\, x^{-\rho}\,\left( \lambda\int_0^\infty \left( 1+ \sigma^{(\a)}_t\right)^{-\rho -\alpha} \, dt\right)\right]\; = \; {\mathbb{P}}\left[ {\F}_{\a,\lambda,\rho} \le x\right]
\end{eqnarray*}
for every $x > 0,$ which completes the proof for $\alpha\in(0,1).$ As for Theorem \ref{Fresh} the remaining cases $\alpha =0$ and $\alpha =1$ are elementary and we leave the details to the reader.

\hfill$\square$

\begin{rem}\label{rF1} {\em (a) Contrary to ${\W}_{\a,\lbd,\rho},$ the factor
$$\int_0^\infty \left( 1+ \sigma^{(\a)}_t\right)^{-\rho -\alpha} \, dt$$
appearing in the decomposition of ${\F}_{\a,\lambda,\rho}$ cannot be expressed as a generalized stable law. On the other hand, this factor can also be viewed as the perpetuity of some subordinator: rewriting
$$b_n(\alpha,\rho)\; =\; \frac{n!}{\Psi_{\alpha,\rho}(1)\cdots\Psi_{\alpha,\rho}(n)}$$
with the Bernstein function
\begin{eqnarray*}
\Psi_{\alpha,\rho}(x) \; = \; \frac{\Ga(\alpha +\rho x)}{\Ga(\rho x)} & = & \frac{\alpha}{\rho\Ga(1-\alpha)}\int_0^\infty (1- e^{-x u}) \frac{e^{-u\rho^{-1}}}{(1-e^{-u\rho^{-1}})^{\alpha +1}}\, du,
\end{eqnarray*}
we obtain as above
\begin{equation}
\label{Tzeta}
\frac{\Ga(\rho)}{\Ga(\rho+\alpha)}\; {\bf T}(1+\alpha\rho^{-1},\rho^{-1}, (1-\alpha)\rho^{-1})\;\stackrel{d}{=}\;\int_0^\infty e^{-\zeta^{(\alpha,\rho)}_t}\, dt
\end{equation}
where $\left\{ \zeta^{(\alpha,\rho)}_t\! , \; t\ge 0\right\}$ is the subordinator having Laplace exponent
$$\esp\left[ e^{-\lambda \zeta^{(\alpha,\rho)}_t}\right]\; = \; e^{-t \Psi_{\alpha,\rho}(\lambda)}.$$
This leads to
$${\F}_{\a,\lambda,\rho}\; \stackrel{d}{=}\; {\bf F}_\rho\;\times\; \left( \lambda \int_0^\infty e^{-\zeta^{(\alpha,\rho)}_t}\, dt \right)^{\frac{1}{\rho}}.$$
Let us again notice that the identification (\ref{Tzeta}) can be deduced from Corollary 5 in \cite{LS} is the case $q = -\rho -\alpha < -\alpha$ and ${\hat \rho} =0$ - see also Remark 10 therein. Observe also that the limiting cases $\alpha = 0$ and $\alpha =1$ are consistent, with respectively $\zeta^{(0,\rho)}_t = \sigma^{(0)}_t$ and $\zeta^{(1,\rho)}_t = \rho\, t.$

\medskip

(b) Since $\Psi_{1-\alpha,\rho}(x)\Phi_{\alpha,\rho}(x) = \rho x,$ we have
$$ n!\; =\;\rho^n  \, \times\, b_n(1-\alpha,\rho)\,\times\, a_n(\alpha,\rho)$$
for all $\alpha\in[0,1],\rho > 0$ and $n\ge 1.$ By moment determinacy, this implies the following factorization of the unit exponential law
\begin{eqnarray}
\label{LRho}
{\bf L} & \stackrel{d}{=} & \left( \int_0^\infty \rho \left( 1+ \sigma_t^{(1-\alpha)}\right)^{-\rho -1+\alpha} \, dt\right) \times \left(\int_0^\infty \left( (1- \sigma^{(\a)}_t)_+\right)^{\rho -\alpha} \, dt\right)\nonumber\\
& \stackrel{d}{=} & \left( \int_0^\infty \left( 1+ \rho^{-1}\sigma_t^{(1-\alpha)}\right)^{-\rho -1+\alpha} \, dt\right) \times \left(\int_0^\infty \left( (1- \rho^{-1}\sigma^{(\a)}_t)_+\right)^{\rho -\alpha} \, dt\right)
\end{eqnarray}
which is valid for all $\alpha\in[0,1]$ and $\rho > 0.$ For $\rho = \alpha,$ this factorization reads
$${\bf L}\;\stackrel{d}{=}\; \left( \int_0^\infty \frac{\alpha\, dt}{1+ \sigma_t^{(1-\alpha)}}\right) \times\; {\bf Z}_\alpha^{-\alpha}\; \stackrel{d}{=}\; {\bf L}^\alpha\,\times\,{\bf Z}_\alpha^{-\alpha},$$
where the first identity follows from Remark \ref{rW1} (d) and the second one from (3.4) in \cite{LS}. The simple identity ${\bf L}\stackrel{d}{=} {\bf L}^\alpha\times{\bf Z}_\alpha^{-\alpha}$ is well-known as Shanbhag-Sreehari's identity. It has been thoroughly discussed in Section 3 of \cite{BY2001} from the point of view of perpetuities of subordinators, and their associated remainders. Observe also that changing the variable and letting $\rho\to\infty$ in (\ref{LRho}) leads to
$${\bf L}\;\stackrel{d}{=}\; \left( \int_0^\infty e^{-\sigma_t^{(1-\alpha)}} dt \right) \times\left( \int_0^\infty e^{-\sigma_t^{(\alpha)}} dt \right),$$
another classic identity obtained in \cite{CPY} - see Example E therein. Last, it is interesting to mention the following identity, which follows at once from (\ref{LRho}), Theorem \ref{Weib} and Theorem \ref{Fresh}:
$${\bf W}_{\alpha,1,\rho}\,\times\,{\bf F}_{1-\alpha,\rho,\rho}^{-1}\;\stackrel{d}{=}\;\left( \frac{{\bf L}\times{\bf L}}{{\bf L}} \right)^{\frac{1}{\rho}}.$$

(c) The above proof shows that the function
$$x\;\mapsto\;E_{\alpha,m,m-\frac{1}{\alpha}}(-x)$$
is CM for every $\alpha\in (0,1]$ and $m >0.$ This is a generalization of the fact that $ E_{\alpha,1,1-\frac{1}{\alpha}}(-x)\; =\; \Ga(\alpha)E_{\alpha,\alpha}(-x)$ is CM for every $\alpha\in (0,1],$ which is itself a direct consequence of the aforementioned Pollard theorem because $\alpha E_\alpha'(-x) = E_{\alpha,\alpha} (-x).$ The formula (\ref{KS2}) implies the Bernstein representation
$$E_{\alpha,m,m-\frac{1}{\alpha}}(-x)\; =\; \esp\left[ \exp -x \left\{ \int_0^\infty \left( 1+ \sigma^{(\a)}_t\right)^{-\alpha (m+1)} \, dt\right\}\right], \qquad x\ge 0.$$
For $m = 1,$ with the notation of Remark \ref{rW1} (d) we obtain
$$E_{\alpha,1,1-\frac{1}{\alpha}}(-x)\; =\; \Ga(\alpha)E_{\alpha,\alpha}(-x) \; =\; \Ga(1+\alpha)E'_{\alpha}(-x)\; =\; \Ga(1+\alpha)\esp \left[ T_\alpha\, e^{-x T_\alpha}\right]\; =\; \esp \left[ e^{-x T_\alpha^{(1)}}\right]$$
where $T^{(1)}_\alpha$ is the size-bias of order 1 of $T_\alpha.$ This implies the curious identity
$$T_\alpha^{(1)} \; \stackrel{d}{=}\; \int_0^\infty \left( 1+ \sigma^{(\a)}_t\right)^{-2\alpha} \, dt.$$
In a different direction, is is worth recalling that for every $\beta > \alpha$ and $\alpha\in (0,1]$ the function
$$E_{\alpha,1,\frac{\beta-1}{\alpha}}(-x)\; =\; \Ga(\beta)E_{\alpha,\beta}(-x)\; =\;\esp \left[ e^{-x \,{\bf B}_{\alpha,\beta -\alpha}^\alpha \times\, T_\alpha^{(1)}}\right]$$
is also CM, where the Bernstein representation involving the $\alpha-$power of a standard Beta distribution follows directly from Lemma 4.26 in \cite{GKMR} - see also the references therein. In Section \ref{FurthA} below, we will come back to this example together with further analytical properties of the Kilbas-Saigo functions $E_{\alpha,m,m-\frac{1}{\alpha}}.$}

\end{rem}

We end this section with a convergent series representation in the non-explicit case $\alpha\in(0,1)$ for the density $f^{\bf F}_{\alpha,\lambda,\rho}$ of ${\F}_{\a,\lambda,\rho}.$ This is a consequence of a term-by-term differentiation of the distribution function
$${\mathbb{P}}[{\F}_{\a,\lambda,\rho} \le x]\; =\; E_{\alpha,\frac{\rho}{\alpha},\frac{\rho-1}{\alpha}}(-\lambda x^{-\rho})$$
which was obtained during the proof of Theorem \ref{Fresh}.

\begin{coro}
\label{FreshD}
For every $\alpha\in (0,1),$ the density of ${\F}_{\a,\lambda,\rho}$ has the following convergent series representation on $(0,\infty)$
$$f^{\bf F}_{\alpha,\lambda,\rho}(x)\; =\; \lambda \, x^{-\rho -1}\,\sum_{n\ge 0}\left( \prod_{j=0}^n \frac{\Ga(j\rho + \rho + 1)}{\Ga(j\rho +\rho +\alpha)}\right) \frac{(-\lambda x^{-\rho})^n}{\rho^n\,n!}\cdot$$
\end{coro}

\subsection{Proof of Theorem \ref{Gumb}}

The argument is shorter than for Theorems \ref{Weib} and \ref{Fresh}. We first consider the case $\alpha\in (0,1).$ The uniqueness is a direct consequence of Theorem \ref{t1G}. We next compute, by Fubini's theorem, the survival function
$${\mathbb{P}}[\Gal > x]\; =\; {\mathbb{P}}[{\bf L} > e^{\lambda x +{\bf G}_\alpha}] \; = \; \esp\left[ e^{-e^{\lambda x +{\bf G}_\alpha}}\right]\; =\; \sum_{n\ge 0} \frac{(-1)^n e^{\lambda n x}}{n!}\, \esp\left[ e^{n{\bf G}_\alpha}\right]$$
for every $x\in{\mathbb{R}}.$ On the other hand, since
$$e^{{\bf G}_\alpha}\;\stackrel{d}{=}\;\int_0^\infty e^{-\sigma^{(\a)}_t}\, dt,$$
we know from Proposition 3.3 in \cite{CPY} that $\esp\left[ e^{n{\bf G}_\alpha}\right] = (n!)^{1-\alpha}$ for all $n\ge 0.$
This implies
\begin{equation}
\label{LREM}
{\bar F} (x)\; =\; {\mathbb{P}}[\Gal > x]\; =\;\sum_{n\ge 0} \frac{(-1)^n e^{\lambda n x}}{(n!)^\alpha}\; =\; \cL_\alpha(-e^{\lambda x}).
\end{equation}
A direct integration based on the Gamma integral and Fubini's theorem shows finally that
$$\Dp F (x)\; =\; \lambda^\alpha e^{\lambda x} {\bar F}(x)$$
for every $x\in{\mathbb{R}}$ as required. The case $\alpha = 1$ was already discussed in the introduction with a unique solution $F(x) = {\mathbb{P}}[{\bf G} \le \lambda x] = {\mathbb{P}}[{\bf G}_{1,\lambda} \le x]$, whereas the unique solution in the case $\alpha = 0$ is obviously $F(x) = 1/(1+ e^{-\lambda x}),$ which is the distribution function of $\lambda^{-1} ({\bf G} -{\bf G}) = {\bf G}_{0,\lambda}.$

\hfill$\square$

\begin{rem}
\label{rG1}
{\em (a) The above proof also shows that the unique distribution function solving the fractional differential equation
$$\Dm {\bar F} (x)\; =\; \lambda^\alpha e^{-\lambda x} F(x), \qquad F(x) > 0\;\mbox{on ${\mathbb{R}},$}$$
is $F(x) = \cL_\alpha(-e^{-\lambda x}) = {\mathbb{P}}[-\Gal \le x].$

\medskip

(b) It is easy to deduce from the representations of the fractional extreme distributions in terms of integrals of the stable subordinator the following convergences in law
$$\rho\lambda^{-1}\left({\bf W}_{\alpha, \rho^\alpha\!, \rho} - 1\right)\;\stackrel{d}{\longrightarrow}\;\Gal\qquad\mbox{and}\qquad \rho\lambda^{-1}\left( 1- {\bf F}_{\alpha,\rho^\alpha\! ,\rho}\right)\;\stackrel{d}{\longrightarrow}\;\Gal$$
as $\rho\to\infty,$ for every $\alpha\in [0,1]$ and $\lambda > 0.$ Observe that the case $\alpha=\lambda =1$ amounts to the aforementioned convergences in law $\rho({\bf W}_\rho -1)\stackrel{d}{\longrightarrow}{\bf G}$ and $\rho(1-{\bf F}_\rho)\stackrel{d}{\longrightarrow}{\bf G}.$}

\end{rem}

As above, we finish this paragraph with a convergent series representation in the non-explicit case $\alpha\in(0,1)$ for the density $f^{\bf G}_{\alpha,\lambda}$ of $\Gal,$ which is a consequence of a term-by-term differentiation of the survival function ${\mathbb{P}}[\Gal > x] = \cL_{\alpha}(-e^{\lambda x}).$

\begin{coro}
\label{GumbD}
For every $\alpha\in (0,1],$ the density of $\Gal$ has the following convergent series representation on ${\mathbb{R}}$
$$f^{\bf G}_{\alpha,\lambda}(x)\; =\; \sum_{n\ge 1} \frac{(-1)^{n-1} \lambda n\, e^{\lambda n x}}{(n!)^\alpha}\cdot$$
\end{coro}

\section{Further properties}

\label{FurthA}

\subsection{On the complete monotonicity of the Kilbas-Saigo function} In this paragraph, motivated by the previous examples arising from the fractional Weibull and Fr\'echet distributions, we wish to characterize the complete monotonicity of all functions $x\mapsto E_{\alpha,m,l}(-x)$ on $(0,\infty).$ We begin with the following result on generalized Pochhammer symbols, which is reminiscent of Proposition 5.1 and Theorem 6.2 in \cite{D2010} and has an independent interest.

\begin{lem}
\label{Gtype}
Let $a,b,c,d$ and $\delta$ be positive parameters. There exists a positive random variable $Z ={\bf Z}[a,c;b,d;\delta]$ such that
\begin{equation}
\label{Gdelta}
\esp[Z^s]\; =\; \frac{[a;\delta]_s[c;\delta]_s}{[b;\delta]_s[d;\delta]_s}
\end{equation}
for every $s > 0,$ if and only if $b+d \le a+c$ and $\inf\{b,d\} \le \inf\{a, c\}.$ This random variable is absolutely continuous on $(0,\infty)$, except in the degenerate case $a=b=c=d.$ Its support is $[0,1]$ if $b+d = a+c$ and $[0,\infty)$ if $b+d < a+c.$
\end{lem}

\proof

We discard the degenerate case $a=b=c=d,$ which is obvious with $Z=1.$ By (\ref{Malm}) and some rearrangements - see also (2.15) in \cite{LS}, we first rewrite
$$\log\left( \frac{[a;\delta]_s[c;\delta]_s}{[b;\delta]_s[d;\delta]_s}\right) \; =\; \kappa\, s\; + \; \int_{-\infty}^0 (e^{sx} - 1 - sx) \left( \frac{e^{- b\vert x\vert} + e^{- d\vert x\vert} - e^{- a\vert x\vert} - e^{- c\vert x\vert}}{\vert x\vert (1-e^{-\vert x\vert})(1-e^{-\delta \vert x\vert})}\right) dx$$
for every $s > 0,$ where $\kappa$ is some real constant. By convexity, it is easy to see that if $b+d \le a+c$ and $\inf\{b,d\} \le \inf\{a, c\},$ then the function $z\mapsto z^b + z^d - z^a - z^c$ is positive on $(0,1).$ This implies that the function
$$x\;\mapsto\;  \frac{e^{- b\vert x\vert} + e^{- d\vert x\vert} - e^{- a\vert x\vert} - e^{- c\vert x\vert}}{\vert x\vert (1-e^{-\vert x\vert})(1-e^{-\delta \vert x\vert})}$$
is positive on $(-\infty, 0)$ and that it can be viewed as the density of some L\'evy measure on $(-\infty,0),$ since it integrates $1\wedge x^2.$ By the L\'evy-Khintchine formula, there exists a real infinitely divisible random variable $Y$ such that
$$\esp[e^{sY}] \; =\;\frac{[a;\delta]_s[c;\delta]_s}{[b;\delta]_s[d;\delta]_s}$$
for every $s > 0,$ and the positive random variable $Z = e^Y$ satisfies (\ref{Gdelta}). Since we have excluded the degenerate case, the L\'evy measure of $Y$ is clearly infinite and it follows from Theorem 27.7 in \cite{S1999} that $Y$ has a density and the same is true for $Z.$

Assuming first $b+d = a+c,$ a Taylor expansion at zero shows that the density of the L\'evy measure of $Y$ integrates $1\wedge \vert x\vert$ and we deduce from (\ref{Malm}) the simpler formula
$$\log\esp[e^{sY}]\; =\; \log\left( \frac{[a;\delta]_s[c;\delta]_s}{[b;\delta]_s[d;\delta]_s}\right) \; =\; -\int^{\infty}_0 (1 - e^{-sx}) \left( \frac{e^{- bx} + e^{- dx} - e^{- ax} - e^{- cx}}{x(1-e^{-x})(1-e^{-\delta x})}\right) dx.$$
By the L\'evy-Khintchine formula, this shows that the ID random variable $Y$ is negative. Moreover, its support is $(-\infty,0]$ since its L\'evy measure has full support and its drift coefficient is zero - see Theorem 24.10 (iii) in \cite{S1999}, so that the support of $Z$ is $[0,1].$

Assuming second $b+d < a+c,$ the same Taylor expansion as above shows that the density of the L\'evy measure of $Y$ does not integrate $1\wedge \vert x\vert$ and the real L\'evy process associated to $Y$ is hence of type C with the terminology of \cite{S1999} - see Definition 11.9 therein. By Theorem 24.10 (i) in \cite{S1999}, this implies that $Y$ has full support on ${\mathbb{R}},$ and so does $Z$ on ${\mathbb{R}}^+\!.$ \\

It remains to prove the only if part of the Lemma. Assuming $a \le d$ and $b\le c$ without loss of generality, we first observe that if $a < b$ then the function
$$s\;\mapsto\; \frac{[a;\delta]_s[c;\delta]_s}{[b;\delta]_s[d;\delta]_s}$$
is real-analytic on $(-b,\infty)$ and vanishes at $s = - a > -b,$ an impossible property for the Mellin transform of a positive random variable. The necessity of $b+d\le a+c$ is slightly more subtle and hinges again upon infinite divisibility. First, setting $\varphi(z) = z^b + z^d - z^a - z^c$ and $z_* =\inf\{ z > 0, \; \varphi(z) < 0\},$ it is easy to see by convexity and a Taylor expansion at 1 that if $b+d > a+c,$ then $z_* < 1$ and $\varphi(z) < 0$ on $(z_*,1)$ with $\varphi(z)\sim (b+d-a-c)(z-1)$ as $z\to 1.$ Introducing next the ID random variable $V$ with Laplace exponent
$$\log\esp[e^{sV}]\; =\; - \kappa\, s\; + \; \int_{\log z_*}^0 (e^{sx} - 1 - sx) \left( \frac{e^{- a\vert x\vert} + e^{- c\vert x\vert} - e^{- b\vert x\vert} - e^{- d\vert x\vert}}{\vert x\vert (1-e^{-\vert x\vert})(1-e^{-\delta \vert x\vert})}\right) dx,$$
we obtain the decomposition
$$\log\left( \frac{[a;\delta]_s[c;\delta]_s}{[b;\delta]_s[d;\delta]_s}\right)\; +\; \log\esp[e^{sV}] \; =\;\int_{-\infty}^{\log z_*} (e^{sx} - 1 - sx) \left( \frac{e^{- b\vert x\vert} + e^{- d\vert x\vert} - e^{- a\vert x\vert} - e^{- c\vert x\vert}}{\vert x\vert (1-e^{-\vert x\vert})(1-e^{-\delta \vert x\vert})}\right) dx,$$
whose right-hand side is the Laplace exponent of some ID random variable $U$ having an atom because its L\'evy measure, whose support is bounded away from zero, is finite - see Theorem 27.4 in \cite{S1999}. On the other hand, the random variable $V$ has an absolutely continous and infinite L\'evy measure and hence it has also a density. If there existed $Z$ such that (\ref{Gdelta}) holds, then the independent decomposition $U\stackrel{d}{=} V  + \log Z$ would imply by convolution that $U$ has a density as well. This contradiction finishes the proof of the Lemma.

\endproof

\begin{rem} {\em (a) By the Mellin inversion formula, the density of ${\bf Z}[a,c;b,d;\delta]$ is expressed as
$$f(x) \; =\; \frac{1}{2\ii\pi x} \int_{s_0 -\ii \infty}^{s_0 + \ii\infty} x^{-s} \left( \frac{[a;\delta]_s[c;\delta]_s}{[b;\delta]_s[d;\delta]_s} \right) \, ds$$
over $(0,\infty)$ for any $s_0 > -\inf\{b,d\}.$ From this expression, it is possible to prove that this density is real-analytic over the interior of the support. We omit details. Let us also mention by Remark 28.8 in \cite{S1999} that this density is positive over the interior of its support.

\medskip

(b) With the standard notation for Pochhammer symbols, the aforementioned Proposition 5.1 and Theorem 6.2 in \cite{D2010} show that
$$s\;\mapsto\;\frac{(a)_s(c)_s}{(b)_s(d)_s}$$
is the Mellin transform of a positive random variable if and only if $b+d \ge a+c$ and $\inf\{b,d\} \ge \inf\{a, c\}.$ This fact can be proved exactly as above, in writing
$$\log\left( \frac{(a)_s(c)_s}{(b)_s(d)_s}\right) \; =\; - \int^{\infty}_0 (1-e^{-sx}) \left( \frac{e^{- ax} + e^{- c x} - e^{- b x} - e^{- d x}}{x(1-e^{-x})}\right) dx.$$
This expression also shows that the underlying random variable has support $[0,1]$ and that it is absolutely continuous, save for $a+c = b+d$ where it has an atom at zero. We refer to \cite{D2010} for an exact expression of the density on $(0,1)$ in terms of the classical hypergeometric function.}

\end{rem}

We can now state the main result of this paragraph, which characterizes the CM property for $E_{\alpha,m,l}(-x)$ on $(0,\infty).$

\begin{prop}
\label{KSCM}
Let $\alpha, m > 0$ and $l > -1/\alpha.$ The Kilbas-Saigo function
$$x\;\mapsto\; E_{\alpha,m,l}(-x)$$
is {\em CM} on $(0,\infty)$ if and only if $\alpha\le 1$ and $l\ge m-1/\alpha.$ Its Bernstein representation is
\begin{equation}
\label{BKS}
E_{\alpha,m,l}(-x)\; =\; \esp\left[ \exp -x \left\{ {\bf X}_{\alpha,m,l}\;\times\int_0^\infty \left( 1+ \sigma^{(\a)}_t\right)^{-\alpha (m+1)} \, dt\right\}\right]
\end{equation}
with $\delta =1/\alpha m$ and ${\bf X}_{\alpha,m,l} = {\bf Z}[1+1/m, (\alpha l +1)\delta;1, 1/m + (\alpha l +1)\delta;\delta].$
\end{prop}

\proof
Assume first $l\ge m-1/\alpha$ and let
$${\bf Y}_{\alpha,m,l}\;=\; {\bf X}_{\alpha,m,l}\;\times\int_0^\infty \left( 1+ \sigma^{(\a)}_t\right)^{-\alpha (m+1)}.$$
By Lemma \ref{Gtype} and Proposition 2 in \cite{LS}, its Mellin transform is
\begin{eqnarray*}
\esp[({\bf Y}_{\alpha,m,l})^s] & = &  \delta^{s}\, \frac{[1+\delta;\delta]_s[(\alpha l +1)\delta;\delta]_s}{[1;\delta]_s[1/m +(\alpha l +1)\delta;\delta]_s}\\
& = & \Ga(1+s)\,\times\frac{[(\alpha l +1)\delta;\delta]_s}{[1/m +(\alpha l+1)\delta;\delta]_s}
\end{eqnarray*}
where in the second equality we have used (\ref{Conca4}). By Fubini's theorem, the moment generating function of ${\bf Y}_{\alpha,m,l}$ reads
\begin{eqnarray*}
\esp[e^{z {\bf Y}_{\alpha,m,l}}] & = & \sum_{n\ge 0}\; \esp[({\bf Y}_{\delta,m,\eta})^n]\,\frac{z^n}{n!}\\
& = & \sum_{n\ge 0} \left(\frac{[(\alpha l +1)\delta;\delta]_n}{[1/m +(\alpha l +1)\delta;\delta]_n}\right) z^n\\
& = &  \sum_{n\ge 0} \left(\prod_{j=0}^{n-1}\frac{\Ga(\alpha(jm +l)+1)}{\Ga(\alpha(jm +l +1)+1)}\right) z^n\; =\; E_{\alpha,m,l}(z)
\end{eqnarray*}
for every $z\ge 0,$ where in the third equality we have used (\ref{Conca1}) repeatedly. The latter identity is extended analytically to the whole complex plane and we get, in particular,
$$ E_{\alpha,m,l}(-x) \; =\; \esp[e^{-x {\bf Y}_{\alpha,m,l}}],\qquad x\ge 0.$$
This shows that $E_{\alpha,m,l}(-x)$ is CM with the required Bernstein representation.\\

We now prove the only if part. If $E_{\alpha,m,l}(-x)$ is CM, then we see by analytic continuation that $E_{\alpha,m,l}(z)$ is the moment generating function on $\CC$ of the underlying random variable $X,$ whose positive integer moments read
$$\esp[X^n]\; =\; n!\,\times\left(\prod_{j=0}^{n-1}\frac{\Ga(\alpha(jm +l)+1)}{\Ga(\alpha(jm +l +1)+1)}\right),\quad n\ge 0.$$
If $\alpha > 1,$ Stirling's formula implies $\esp[X^n]^{\frac{1}{n}}\to 0$ as $n\to\infty$ so that $X\equiv 0,$ a contradiction because $E_{\alpha,m,l}$ is not a constant. If $\alpha = 1$ and $l+1 < m,$ then
$$\esp[X^n]\; = \;\frac{n!}{(c)_n m^n}\; \sim\; \frac{n^{1-c}}{m^n}\quad \mbox{as $n\to\infty,$}$$
with $c = (l+1)/m \in (0,1).$ In particular, the Mellin transform $s\mapsto\esp[X^s]$ is analytic on $\{\Re(s)\ge 0\},$ bounded on $\{\Re(s) = 0\},$ and has at most exponential growth on $\{\Re(s) > 0\}$ because
$$\vert \esp[X^s]\vert\; \le\; \esp\left[ X^{\Re(s)}\right]\; =\; \left( \esp\left[ X^{[\Re(s)] +1}\right]\right)^{\frac{\Re(s)}{[\Re(s)] +1}}$$
by H\"older's inequality. On the other hand, the Stirling type formula (\ref{Stirling}) implies, after some simplifications,
$$\delta^{-s} \frac{[1+\delta;\delta]_s[c;\delta]_s}{[1;\delta]_s[c+\delta;\delta]_s}\; = \; \delta^{-s} s^{1-c} (1+ o(1))\quad \mbox{as $\vert s\vert \to\infty$ with $\vert \arg s\vert < \pi$}$$
and this shows that the function on the left-hand side, which is analytic on $\{\Re(s)\ge 0\},$ has at most linear growth on $\{\Re(s) = 0\}$ and at most exponential growth on $\{\Re(s) > 0\}.$ Moreover, the above analysis clearly shows that
$$\esp[X^n] \; = \; \delta^{-n} \frac{[1+\delta;\delta]_n[c;\delta]_n}{[1;\delta]_n[c+\delta;\delta]_n}$$
for all $n\ge 0$ and by Carlson's theorem - see e.g. Section 5.81 in \cite{T1939}, we must have
$$\esp[X^s] \; = \; \delta^{-s} \frac{[1+\delta;\delta]_s[c;\delta]_s}{[1;\delta]_s[c+\delta;\delta]_s}$$
for every $s > 0,$ a contradiction since Lemma \ref{Gtype} shows that the right-hand side cannot be the Mellin transform of a positive random variable if $\eta < 0.$ The case $\alpha < 1$ and  $l+1/\alpha < m$ is analogous. It consists in identifying the bounded sequence
$$ \frac{1}{n!}\,\times\left(\prod_{j=0}^{n-1}\frac{\Ga(\alpha(jm +l +1)+1)}{\Ga(\alpha(jm +l)+1)}\right)$$
as the values at non-negative integer points of the function
$$\delta^{-s}\times \frac{[1;\delta]_s[1/m + (\alpha l +1)\delta;\delta]_s}{[1+\delta;\delta]_s[(\alpha l +1)\delta;\delta]_s}\; =\; \delta^{-s} e^{-(1-\alpha) s \ln (s) + \kappa s + O(1)}\quad\mbox{as $\vert s\vert \to \infty$ with $\vert \arg s\vert < \pi,$}$$
where the purposeless constant $\kappa$ can be evaluated from (\ref{Stirling}). On $\{\Re(s) \ge 0\},$ we see that this function has growth at most $e^{\pi (1-\alpha) \vert s\vert /2}$ and we can again apply Carlson's theorem. We omit details.

\endproof

\begin{rem}

\label{RCM}

{\em (a) When $m =1,$ applying (\ref{Conca1}) we see that the random variable ${\bf X}_{\alpha, 1,l}$ has Mellin transform
$$\esp[({\bf X}_{\alpha,1,l})^s]\; =\; \frac{[2;\delta]_s[l+1/\alpha;\delta]_s}{[1;\delta]_s[1 + l + 1/\alpha;\delta]_s}\; =\; \frac{(\alpha)_{\alpha s}}{(\beta)_{\alpha s}}$$
with $\beta = 1+\alpha l\ge \alpha.$ This shows ${\bf X}_{\alpha, 1,l}\stackrel{d}{=} {\bf B}_{\alpha,\beta -\alpha}^\alpha,$ and we recover the Bernstein representation of the CM function $\Ga(\beta) E_{\alpha,\beta}(-x)$ which was discussed in Remark \ref{rF1} (c). Notice also the very simple expression for the Mellin transform
$$\esp[({\bf Y}_{\alpha,1,l})^s] \; =\; \frac{\Ga(1+\alpha l)\Ga(1+s)}{\Ga(1 +\alpha (l+s))}\cdot$$

(b) Another simplification occurs when $l+1/\alpha = k m$ for some integer $k\ge 1.$ One finds
$$\esp[({\bf X}_{\alpha,m,km -1/\alpha})^s]\; =\; \frac{[k;\delta]_s[1 +1/m;\delta]_s}{[1;\delta]_s[k +1/m;\delta]_s}\; =\; \prod_{j=1}^{k-1} \left( \frac{(\alpha j m)_{u}}{(\alpha(jm+1))_{u}}\right)$$
for $u =\alpha m s \ge 0,$ which implies
$${\bf X}_{\alpha, m, km - 1/\alpha}\; \stackrel{d}{=}\; \left( {\bf B}_{\alpha m, \alpha}\times\cdots\times {\bf B}_{\alpha m (k-1), \alpha}\right)^{\alpha m}.$$
In general, the law of the absolutely continuous random variable ${\bf X}_{\alpha,m,l}$ valued in $[0,1]$ seems to have a complicated expression.

\medskip

(c) As seen during the proof, the random variable ${\bf Y}_{\alpha,m,l}$ defined by the Bersntein representation
$$E_{\alpha, m, l}(-x) \; =\; \esp [e^{-x {\bf Y}_{\alpha, m, l}}]$$
has Mellin transform
\begin{equation}
\label{MellY}
\esp[({\bf Y}_{\alpha,m,l})^s]\; =\; \Ga(1+s)\, \times\frac{[(\alpha l +1)\delta;\delta]_s}{[1/m +(\alpha l+1)\delta;\delta]_s}
\end{equation}
with $\delta = 1/\alpha m,$ for every $s > -1.$ By Fubini's theorem, this implies
$$\int_0^\infty E_{\alpha, m, l}(-x)\, x^{s-1}\, dx \; = \;  \Ga(s)\,\esp [{\bf Y}_{\alpha,m,l}^{-s}]\; =\; \Ga(s)\Ga(1-s)\, \times\frac{[(\alpha l +1)\delta;\delta]_{-s}}{[1/m +(\alpha l+1)\delta;\delta]_{-s}}$$
for every $s\in (0,1).$ Notice that this formula, which seems unnoticed in the literature on the Kilbas-Saigo function, remains valid for $l \in(-1/\alpha, m-1/\alpha)$ by the analyticity of the map $l\mapsto E_{\alpha, m, l} (-x).$ For $m = 1,$ we recover from (\ref{Conca1}) the formula
$$\int_0^\infty E_{\alpha, \beta}(-x)\, x^{s-1}\, dx \; = \; \frac{1}{\Ga(\beta)} \int_0^\infty E_{\alpha,1, \frac{\beta-1}{\alpha}}(-x)\, x^{s-1}\, dx \; = \;\frac{\Ga(s)\Ga(1-s)}{\Ga(\beta - \alpha s)}$$
which is given in (4.10.3) of \cite{GKMR}, as a consequence of the Mellin-Barnes representation of $E_{\alpha,\beta}(z).$ Notice that there is no such Mellin-Barnes representation for $E_{\alpha,m,l}(z)$ in general.}

\end{rem}

\subsection{Asymptotic behaviour of the densities}

\label{ABTD}

In this paragraph we study the behaviour of the density functions of the fractional Weibull and Fr\'echet distributions at both ends of their support. To this end, we evaluate the Mellin transforms of ${\W}_{\a,\lbd,\rho}$ and ${\F}_{\a,\lambda,\rho}.$ The case of the fractional Gumbel distribution requires different arguments and will be handled in Paragraph \ref{SCLR}.

\subsubsection{The Weibull case} As a consequence of Corollary \ref{WeibD}, we first obtain
$$f^{\bf W}_{\alpha,\lambda, \rho}(x) \; \sim\; \left(\frac{\lambda\,\Ga(\rho +1-\alpha)}{\Ga(\rho)}\right) x^{\rho -1}\quad\mbox{as $x\to 0.$}$$
The behaviour of the density at infinity is less immediate and we will need the exact expression of the Mellin transform of ${\W}_{\a,\lbd,\rho},$ which has an interest in itself.

\begin{prop}
\label{MASW}
The Mellin transform of ${\W}_{\a,\lbd,\rho}$ is
$$\esp\left[{\W}_{\a,\lbd,\rho}^s\right]\; =\; \left(\frac{\rho^\alpha}{\lambda}\right)^{\frac{s}{\rho}}\Ga(1+ s\rho^{-1})\,\times\,\frac{[\rho +(1-\alpha);\rho]_{-s}}{[\rho;\rho]_{-s}}$$
for every $s\in (-\rho,\rho).$ As a consequence, one has
$$f^{\bf W}_{\alpha,\lambda, \rho}(x)\; \sim\; \left(\frac{\rho }{\lambda \Ga(1-\alpha)}\right) x^{-\rho -1} \quad \mbox{as $x\to\infty.$}$$
\end{prop}

\proof
We start with a more concise expression of (\ref{MellY}) for $l = m-1,$ which is a direct consequence of (\ref{Conca4}):
$$\esp[({\bf Y}_{\alpha,\frac{\rho}{\alpha},\frac{\rho}{\alpha} -1})^s]\; =\; \rho^{-s} \times\frac{[1+(1-\alpha)\rho^{-1};\rho^{-1}]_s}{[1;\rho^{-1}]_s}\cdot$$
By Theorem \ref{Weib}, we deduce
\begin{eqnarray*}
\esp\left[{\W}_{\a,\lbd,\rho}^s\right] & = &\esp\left[\left( \frac{{\bf L}}{\lambda {\bf Y}_{\alpha,\frac{\rho}{\alpha},\frac{\rho}{\alpha} -1}}\right)^{\frac{s}{\rho}}\right]\\
& = & \left(\frac{\rho}{\lambda}\right)^{\frac{s}{\rho}}\Ga(1+ s\rho^{-1})\,\times\, \frac{[1+(1-\alpha)\rho^{-1};\rho^{-1}]_{-s\rho^{-1}}}{[1;\rho^{-1}]_{-s\rho^{-1}}}\\
& = & \left(\frac{\rho^\alpha}{\lambda}\right)^{\frac{s}{\rho}}\Ga(1+ s\rho^{-1})\,\times\,\frac{[\rho +(1-\alpha);\rho]_{-s}}{[\rho;\rho]_{-s}}
\end{eqnarray*}
for every $s\in (-\rho,\rho)$ as required, where the third equality comes from (\ref{Conca3}). The asymptotic behaviour of the density at infinity is then a standard consequence of Mellin inversion. First, we observe from the above formula and (\ref{zeroes}) that the first positive pole of $s\mapsto \esp\left[{\W}_{\a,\lbd,\rho}^s\right]$ is simple and isolated in the complex plane at $s = \rho,$ with
\begin{eqnarray*}
\esp\left[{\bf W}_{\alpha,\lambda,\rho}^s\right] & \sim & \left(\frac{\rho^\alpha}{\lambda}\right)\,\times\,\frac{[\rho +(1-\alpha);\rho]_{-\rho}}{[\rho;\rho]_{-s}} \\
& \sim & \left(\frac{\rho^{\rho +\alpha}}{\lambda}\right)\,\times\,\frac{[\rho +(1-\alpha);\rho]_{-\rho}}{[2\rho;\rho]_{-\rho}}\,\times\,(\rho)_{-s}\\
& =& \left(\frac{\rho}{\lambda \Ga(1-\alpha)}\right) \times\, \Gamma (\rho-s)\;\sim\; -\left(\frac{\rho}{\lambda \Ga(1-\alpha)}\right)\times\, \frac{1}{s-\rho}
\end{eqnarray*}
as $s\to \rho,$ where the second asymptotics comes from (\ref{Conca4}) and the third equality from (\ref{Conca2}). Therefore, applying Theorem 4 (ii) in \cite{FGD1995} - beware the correction $(\log x)^k \to (\log x)^{k-1}$ to be made in the expansion of $f(x)$ therein, we obtain
$$f^{\bf W}_{\alpha,\lambda,\rho}(x)\; \sim\; \left(\frac{\rho }{\lambda \Ga(1-\alpha)}\right) x^{-\rho -1} \quad \mbox{as $x\to\infty$}$$
as required.
\endproof

\begin{rem}
\label{rW2}
{\em (a) Another proof of the asymptotic behaviour at infinity can be obtained from (\ref{JSW2}). By multiplicative convolution the latter implies, setting $f^\cG_{\alpha,\rho}$ for the density of $\cG(\rho +1-\alpha,1-\alpha),$
\begin{eqnarray*}
f^{\bf W}_{\alpha,\lambda,\rho}(x) & = & \lambda \, x^{\rho-1} \int_0^\infty f^\cG_{\alpha,\rho}(y)\,y^{-\rho}\, e^{-\frac{\lambda}{\rho}(\frac{x}{y})^\rho} dy \\
& = & \left(\frac{\lambda}{\rho}\right)^{\frac{1}{\rho}}  \int_0^\infty f^\cG_{\alpha,\rho}\left( x(\rho \lambda^{-1} t)^{-\frac{1}{\rho}}\right)\,t^{-\frac{1}{\rho}}\,e^{-t}\, dt\\
& \sim & \left(\frac{\rho }{\lambda \Ga(1-\alpha)}\int_0^\infty t \,e^{-t}\, dt \right) x^{-\rho -1}\; =\; \left(\frac{\rho }{\lambda \Ga(1-\alpha)}\right) x^{-\rho -1}
\end{eqnarray*}
as $x\to\infty,$ where for the asymptotics we have used the Proposition in \cite{JSW2018} and a direct integration. This argument does not make use of Mellin inversion and is overall simpler than the above, but it does not convey to the Fr\'echet case.

\medskip

(b) The Mellin transform simplifies for $\alpha = 0$ and $\alpha = 1\! :$ we recover
$$\esp[{\bf W}_{0,\lambda,\rho}^s] \; = \; \lambda^{-\frac{s}{\rho}}\Ga(1+s\rho^{-1})\Ga(1-s\rho^{-1})\qquad\mbox{and}\qquad \esp[{\bf W}_{1,\lambda,\rho}^s] \; = \; \left(\frac{\rho}{\lambda}\right)^{\frac{s}{\rho}}\Ga(1+s\rho^{-1})$$
in accordance with the scaling property ${\W}_{\a,\lbd,\rho} \stackrel{d}{=} \lambda^{-1/\rho}{\bf W}_{\alpha,1,\rho}$ and (\ref{WeibIDs}), where the first equality follows from (\ref{Conca1}) and (\ref{Rhorho}). The Mellin transform takes a simpler form in two other situations.

\begin{itemize}

\item For $\rho=\alpha,$ we obtain from (\ref{MellY}), (\ref{Conca1}) and (\ref{Conca2})
$$\esp[({\bf Y}_{\alpha,1,0})^s]\; =\; \frac{\Ga(1+s)}{\Ga(1+\alpha s)}\; =\; \esp[{\bf Z}_\alpha^{-\alpha s}],$$
in accordance with Remark \ref{rW1} (d). This yields
$${\bf W}_{\alpha,\lambda,\alpha}\;\stackrel{d}{=}\; \left( \frac{{\bf L}}{\lambda {\bf Y}_{\alpha,1,0}}\right)^{\frac{1}{\alpha}}\;\stackrel{d}{=}\; \lambda^{-\frac{1}{\alpha}}\,{\bf Z}_\alpha\,\times\,{\bf L}^{\frac{1}{\alpha}},$$
an identity which was already discussed for $\lambda = 1$ in the introduction as the solution to (\ref{Barre}). The Mellin transform reads
$$\esp[{\bf W}_{\alpha,\lambda,\alpha}^s] \; = \; \lambda^{-\frac{s}{\alpha}}\,\frac{\Ga(1+s\alpha^{-1})\Ga(1-s\alpha^{-1})}{\Ga(1-s)}\cdot$$

\item For $\rho=1-\alpha,$  where we obtain from (\ref{Conca2})
$$\esp[{\bf W}_{1-\rho,\lambda,\rho}^s] \; = \; \left(\frac{\rho}{\lambda}\right)^{\frac{s}{\rho}}\,\frac{\Ga(1+s\rho^{-1})\Ga(\rho -s)}{\Ga(\rho)}\qquad\mbox{and}\qquad {\bf W}_{1-\rho,\lambda,\rho}\; \stackrel{d}{=}\; \left(\frac{\rho}{\lambda}\right)^{\frac{1}{\rho}}\,{\bf L}^{\frac{1}{\rho}}\,\times\,{\bf \Gamma}^{-1}_{\rho}$$
having denoted by ${\bf \Gamma}_t,$ here and throughout, the standard Gamma random variable with parameter $t.$

\end{itemize}

\medskip

(c) Integrating the density and using ${\mathbb{P}}[{\W}_{\a,\lbd,\rho} > x] = E_{\alpha,\frac{\rho}{\alpha}, \frac{\rho}{\alpha}-1}(-\lambda x^\rho),$ we obtain at once the following asymptotic behaviour at infinity for any $\alpha\in (0,1]$ and $m > 0:$
$$E_{\alpha,m, m-1}(-x)\;\sim\; \frac{1}{\Ga(1-\alpha)\, x}\qquad\mbox{as $x\to\infty.$}$$
This behaviour seems unnoticed in the literature on the Kilbas-Saigo function, and turns out to be the same as that of the classical Mittag-Leffler function $E_\alpha(-x)$ - see e.g. (3.4.15) in \cite{GKMR}. It is actually possible to get the behaviour of $E_{\alpha,m, l}(-x)$ at infinity for all $l>-1/\alpha$ with the help of the Mellin transform computed in Remark \ref{RCM} (c). Notice however that the first positive pole might not be simple, for example when $m > 1$ and $l = m-1 -1/\alpha.$ We shall not discuss this behaviour here, save for $l = m-1/\alpha$ in the framework of the Fr\'echet distribution - see Remark \ref{rF2} (c).

\medskip

(d) The four examples discussed in (b) above have Mellin transform expressed in terms of the quotient of a finite number of Gamma functions, making it possible to use a Mellin-Barnes representation of the density in order to get a full asymptotic expansion at infinity. For example, using the standard notation of Definition C.1.1 in \cite{AAR1999}, one obtains from (1.8.28) in \cite{KST} the expansion
$$f^{\bf W}_{\alpha,\lambda,\alpha}(x)\; =\; \lambda\, x^{\alpha-1} \,E_{\alpha,\alpha}(-\lambda x^\alpha)\; \sim\; \sum_{n\ge 1} \frac{n\alpha\, x^{-1-n\alpha}}{\lambda^n\Ga(1-n\alpha)}$$
which is everywhere divergent. Unfortunately, the Mellin transform of ${\W}_{\a,\lbd,\rho}$ might have poles of variable order and there does not seem to exist any general formula for the full asymptotic expansion at infinity of the density of ${\W}_{\a,\lbd,\rho}$.}
\end{rem}

\subsubsection{The Fr\'echet case} As a consequence of Corollary \ref{FreshD}, we first have
$$f^{\bf F}_{\alpha,\lambda, \rho}(x) \; \sim\; \left(\frac{\lambda\,\Ga(\rho +1)}{\Ga(\rho+\alpha)}\right) x^{-\rho -1}\quad\mbox{as $x\to \infty.$}$$
The behaviour of the density at zero is less immediate and we will need, as above, the exact expression of the Mellin transform of ${\F}_{\a,\lambda,\rho},$ whose strip of analyticity is larger than for ${\W}_{\a,\lbd,\rho}.$

\begin{prop}
\label{MASF}
The Mellin transform of ${\F}_{\a,\lambda,\rho}$ is
$$\esp\left[{\F}_{\a,\lambda,\rho}^s\right]\; =\; \left(\frac{\rho^\alpha}{\lambda}\right)^{-\frac{s}{\rho}}\Ga(1- s\rho^{-1})\,\times\,\frac{[\rho +1;\rho]_{s}}{[\rho +\alpha;\rho]_{s}}$$
for every $s\in (-\rho-\alpha,\rho).$ As a consequence, one has
$$f^{\bf F}_{\alpha,\lambda, \rho}(x)\; \sim\; \left(\frac{\rho^{\frac{\alpha^2}{\rho}}(\rho +\alpha) }{\lambda^{1 +\frac{\alpha}{\rho}}}\,\Ga(1+\alpha)\, G(1-\alpha; \rho)\, G(1+\alpha; \rho)\right) x^{\rho +\alpha -1} \quad \mbox{as $x\to 0.$}$$
\end{prop}

\proof The evaluation of the Mellin transform is done as for the fractional Weibull distribution, starting from the more concise expression
$$\esp[({\bf Y}_{\alpha,\frac{\rho}{\alpha},\frac{\rho -1}{\alpha}})^s]\; =\; \rho^{-s} \times\frac{[1+\rho^{-1};\rho^{-1}]_s}{[1 +\alpha\rho^{-1};\rho^{-1}]_s}$$
and writing
$$\esp\left[{\F}_{\a,\lambda,\rho}^s\right] \; = \; \esp\left[\left( \frac{{\bf L}}{\lambda {\bf Y}_{\alpha,\frac{\rho}{\alpha},\frac{\rho-1}{\alpha}}}\right)^{-\frac{s}{\rho}}\right]\; =\;\left(\frac{\rho^\alpha}{\lambda}\right)^{-\frac{s}{\rho}}\Ga(1- s\rho^{-1})\,\times\,\frac{[\rho +1;\rho]_{s}}{[\rho +\alpha;\rho]_{s}}\cdot$$
The asymptotic behaviour of $f^{\bf F}_{\alpha,\lambda, \rho}(x)$ at zero follows then as for that of $f^{\bf W}_{\alpha,\lambda, \rho}(x)$ at infinity, in considering the residue at the first negative pole $s = -(\rho+\alpha)$ which is simple and isolated in the complex plane, applying Theorem 4 (i) in \cite{FGD1995} - with the same correction as above, and making various simplifications. We omit details.

\endproof

\begin{rem}
\label{rF2}{\em (a) Comparing the Mellin transforms, Propositions \ref{MASW} and \ref{MASF} imply the interesting factorization
\begin{equation}
\label{FacWF}
{\W}_{\a,\lbd,\rho}^{-1}\;\stackrel{d}{=}\; {\F}_{\a,\lambda,\rho}\;\times\;{\bf Z}(\rho +1-\alpha,\rho+\alpha;\rho,\rho+1;\rho).
\end{equation}
In general, it follows from Proposition \ref{KSCM} that for every $\alpha, m, \lambda > 0$ and $l > m- 1/\alpha,$ there exists a positive random variable having distribution function $E_{\alpha,m, l}(-\lambda x^{-\alpha m}),$ and which is given by (\ref{MellY}), (\ref{BKS}) and Theorem \ref{Fresh} as the independent product
$${\bf F}_{\alpha,\lambda, \alpha m}\;\times\;({\bf X}_{\alpha,m,l})^{\frac{1}{\alpha m}} \;\stackrel{d}{=}\; {\bf F}_{\alpha,\lambda, \alpha m}\;\times\;{\bf Z}(\alpha l +1,\alpha (m+1);\alpha m,\alpha(l+1) +1;\alpha m),$$
where the identity in law follows from (\ref{Conca3}). In this respect, the fractional Fr\'echet distributions can be viewed as the ``ground state'' distributions associated to the Kilbas-Saigo functions $E_{\alpha,m, l},$ in the limiting case $l = m-1/\alpha.$

\medskip

(b) As above, the Mellin transform simplifies for $\alpha = 0,1\! :$ in accordance with (\ref{FreshIDs}), we get
$$\esp[{\bf F}_{0,\lambda,\rho}^s] \; = \; \lambda^{\frac{s}{\rho}}\,\Ga(1+s\rho^{-1})\Ga(1-s\rho^{-1})\qquad\mbox{and}\qquad \esp[{\bf F}_{1,\lambda,\rho}^s] \; = \; \left(\frac{\lambda}{\rho}\right)^{\frac{s}{\rho}}\Ga(1-s\rho^{-1}).$$
The Mellin transform also takes a simpler form in the same other situations as above.

\begin{itemize}

\item For $\rho=\alpha,$ with
$$\esp[{\bf F}_{\alpha,\lambda,\alpha}^s] \; = \; \lambda^{\frac{s}{\alpha}}\,\frac{\Ga(\alpha)\Ga(1+s\alpha^{-1})\Ga(1-s\alpha^{-1})}{\Ga(\alpha+s)}\cdot$$
This yields the identity ${\bf F}_{\alpha,\lambda,\alpha}\,\stackrel{d}{=}\, \lambda^{\frac{1}{\alpha}}\,({\bf Z}^{-1}_\alpha)^{(\alpha)}\,\times\,{\bf L}^{-\frac{1}{\alpha}},$ which was discussed in the introduction for $\lambda =1$ as the solution to (\ref{Barres}). This is also in accordance with Remark \ref{rF1} (c), since
$$(T_\alpha^{(1)})^{\frac{1}{\alpha}}\;\stackrel{d}{=}\; (({\bf Z}^{-\alpha}_\alpha)^{(1)})^{\frac{1}{\alpha}}\;\stackrel{d}{=}\;({\bf Z}^{-1}_\alpha)^{(\alpha)}.$$
Notice that the constant appearing in the asymptotic behaviour of the density at zero is also simpler: one finds
\begin{equation}
\label{Fala}
f^{\bf F}_{\alpha,\lambda, \alpha}(x)\; \sim\; \left(\frac{2\alpha\,\Ga(1+\alpha)}{\lambda^2\,\Ga(1-\alpha)}\right) x^{2\alpha -1} \quad \mbox{as $x\to 0.$}
\end{equation}

\item For $\rho=1-\alpha,$ with
$$\esp[{\bf F}_{1-\rho,\lambda,\rho}^s] \; = \; \left(\frac{\lambda}{\rho}\right)^{\frac{s}{\rho}}\,\Ga(1-s \rho^{-1})\Ga(1+s)\quad\mbox{and}\quad{\bf F}_{1-\rho,\lambda,\rho}\; \stackrel{d}{=}\; \left(\frac{\lambda}{\rho}\right)^{\frac{1}{\rho}}\,{\bf L}^{-\frac{1}{\rho}}\,\times\,{\bf L}.$$
Here, the density converges at zero to a simple constant: one finds
$$f^{\bf F}_{1-\rho,\lambda, \rho}(x)\; \to\; \left(\frac{\rho}{\lambda}\right)^{\frac{1}{\rho}} \Ga(1+\rho^{-1})\qquad \mbox{as $x\to 0.$}$$
\end{itemize}

\medskip

(c) Integrating the density and using ${\mathbb{P}}[{\F}_{\a,\lambda,\rho} \le x] = E_{\alpha,\frac{\rho}{\alpha}, \frac{\rho -1}{\alpha}}(-\lambda x^{-\rho}),$ we obtain the following asymptotic behaviour at infinity for any $\alpha\in (0,1]$ and $m > 0:$
$$E_{\alpha,m, m-\frac{1}{\alpha}}(-x)\;\sim\; (\alpha m)^{\frac{\alpha}{m}}\Ga(1+\alpha)\, G(1-\alpha; \alpha m)\, G(1+\alpha; \alpha m) \, x^{-1-\frac{1}{m}}\qquad\mbox{as $x\to\infty.$}$$
For $m=1,$ this behaviour matches the first term in the full asymptotic expansion
$$E_{\alpha,1, 1-\frac{1}{\alpha}}(-x)\; =\; \Ga(\alpha) E_{\alpha,\alpha}(-x)\;\sim\; \Ga(\alpha)\sum_{n\ge 1} \frac{(-1)^n}{\Ga(-\alpha n) x^{n+1}}\cdot$$
As for $E_{\alpha,m, m-1}(-x),$ such a full asymptotic expansion seems difficult to obtain for all values of $m.$}
\end{rem}

\subsection{Optimal bounds for the distribution functions}

In Theorem 4 of \cite{TS2014}, the following uniform hyperbolic bounds are obtained for the classical Mittag-Leffler function:
\begin{equation}
\label{MLB}
\frac{1}{1 +\Ga(1-\alpha) x}\; \le\; E_\alpha(-x)\; \le\; \frac{1}{1 +\frac{1}{\Ga(1+\alpha)}\, x}
\end{equation}
for every $\alpha\in[0,1]$ and $x\ge 0.$ The constants in these inequalities are optimal because of the asymptotic behaviours
$$E_\alpha(-x)\; \sim\;\frac{1}{\Ga(1-\alpha) x}\quad\mbox{as $x\to\infty$}\qquad\mbox{and}\qquad 1- E_\alpha(-x)\; \sim\; \frac{x}{\Ga(1+\alpha)} \quad\mbox{as $x\to 0.$}$$
In this paragraph, we shall obtain analogous bounds for the Kilbas-Saigo functions $E_{\alpha,m,m-1}(-x)$ and $E_{\alpha,m,m-\frac{1}{\alpha}}(-x),$ which are associated to the fractional Weibull and Fr\'echet distributions. We begin with the following monotonicity properties, of independent interest.

\begin{prop}
\label{Mono}
Fix $\alpha\in (0,1]$ and $x\in{\mathbb{R}}.$ The functions
$$m\; \mapsto\; E_{\alpha,m,m-1}(x)\qquad\mbox{and}\qquad m\; \mapsto\; E_{\alpha,m,m-\frac{1}{\alpha}}(x)$$
are decreasing on $(0,\infty)$ if $x > 0$ and increasing on $(0,\infty)$ if $x < 0.$
\end{prop}

\proof

 This is a direct consequence of (\ref{KS1}), (\ref{KS2}), and the fact that $\sigma^{(\a)}_t > 0$ for every $t > 0.$

\endproof

\begin{rem}
{\em It would be interesting to know if the same property holds for $m \mapsto E_{\alpha,m,m - l}(x)$ and any $l \le 1/\alpha.$ However, only the cases $l = 1$ and $l =1/\alpha$ seem to involve the $\alpha-$stable subordinator in a direct way.}
\end{rem}

As in \cite{TS2014}, our analysis to obtain the uniform bounds will use some notions of stochastic ordering. Recall that if $X,Y$ are real random variables such that $\esp[\varphi(X)] \le \esp[\varphi(Y)]$ for every $\varphi :{\mathbb{R}} \to{\mathbb{R}}$ convex, then $Y$ is said to dominate $X$ for the convex order, a property which we denote by $X\prec_{cx} Y.$ The following result on convex orderings for infinite Beta products has an independent interest.

\begin{lem}
\label{Tcx}
For every $a,b,c > 0$ and $d \ge c,$ one has
$${\bf T}(a,b,c)\;\prec_{cx}\; {\bf T}(a,b,d).$$
\end{lem}

\proof

By the definition of ${\bf T}(a,b,c)$ and the stability of the convex order by mixtures - see Corollary 3.A.22 in \cite{SSh}, it is enough to show
$$(a+b) {\bf B}_{a,b}\;\prec_{cx}\;(a+c) {\bf B}_{a,c}$$
for every $a,b > 0$ and $c \ge b.$ Using again Corollary 3.A.22 in \cite{SSh} and the standard identity ${\bf B}_{a,c}\,\stackrel{d}{=}\,{\bf B}_{a,b}\times{\bf B}_{a+b,c-b},$ we are reduced to
$$\left(\frac{a+b}{a+c}\right)\; =\; \esp[{\bf B}_{a+b,c-b}]\;\prec_{cx}\;{\bf B}_{a+b,c-b}$$
which is a direct consequence of Jensen's inequality.

\endproof

The following result is a generalization of (\ref{MLB}), which deals with the case $m=1$ only, to all Kilbas-Saigo functions $E_{\alpha,m,m-1}.$ The argument is much simpler than in the original proof of (\ref{MLB}).

\begin{prop}
\label{KSO1}
For every $\alpha\in[0,1], m >0$ and $x\ge 0,$ one has
$$\frac{1}{1 +\Ga(1-\alpha) x}\; \le\; E_{\alpha,m,m-1}(-x)\; \le\; \frac{1}{1 +\frac{\Ga(1+\alpha(m-1))}{\Ga(1+\alpha m)}\, x}\cdot$$
\end{prop}

\proof

The first inequality is a direct consequence of Proposition \ref{Mono}, which implies
\begin{eqnarray*}
E_{\alpha,m,m-1}(-x) & \ge & \esp\left[\exp \left\{ -x\, \int_0^\infty \left( 1- \sigma^{(\a)}_t\right)_+^{-\alpha} \, dt\right\}\right]\\
& = & \esp\left[ e^{-x\,\Ga(1-\alpha)\,{\bf L}}\right]\; = \; \frac{1}{1 +\Ga(1-\alpha) x}
\end{eqnarray*}
for $x\ge 0,$ where the first equality follows from Theorem (b) (ii) in \cite{LS}. For the second equality, we come back to the infinite product representation
\begin{eqnarray*}
\int_0^\infty \left( 1- \sigma^{(\a)}_t\right)_+^{\rho -\alpha} \, dt & \stackrel{d}{=} & \frac{\Ga(\rho +1-\alpha)}{\Ga(\rho +1)}\; {\bf T}(1,\rho^{-1}, (1-\alpha)\rho^{-1})
\end{eqnarray*}
which follows from Theorem (b) (i) in \cite{LS}, as in the proof of Theorem \ref{Weib}. Lemma \ref{Tcx} implies then
\begin{eqnarray*}
\int_0^\infty \left( 1- \sigma^{(\a)}_t\right)_+^{\rho -\alpha} \, dt & \prec_{cx} & \frac{\Ga(\rho +1-\alpha)}{\Ga(\rho +1)}\; {\bf T}(1,\rho^{-1}, \rho^{-1})\; \stackrel{d}{=}\; \frac{\Ga(\rho +1-\alpha)}{\Ga(\rho +1)}\;{\bf L}
\end{eqnarray*}
where the identity in law follows from (2.7) in \cite{LS}. Using (\ref{KS1}) with $\rho = \alpha m$ and the convexity of $t\mapsto e^{-xt}$, we obtain the required
$$E_{\alpha,m,m-1}(-x)\; \le\; \frac{1}{1 +\frac{\Ga(1+\alpha(m-1))}{\Ga(1+\alpha m)}\, x}\cdot$$
\endproof

\begin{rem}
\label{rKS1}
{\em (a) As for the classical case $m=1$, these bounds are optimal because of the aforementioned asymptotic behaviours
$$E_{\alpha, m, m-1}(-x)\; \sim\;\frac{1}{\Ga(1-\alpha) x}\;\mbox{as $x\to\infty$}\quad\mbox{and}\quad 1\,-\, E_{\alpha, m, m-1}(-x)\; \sim\; \frac{\Ga(1+\alpha (m-1))}{\Ga(1+\alpha m)}\, x \;\mbox{as $x\to 0.$}$$

\medskip

(b) It is easy to check that the above proof also implies
$$E_{\alpha,m,m-1}(x)\; \le\; \frac{1}{(1 -\Ga(1-\alpha) x)_+}$$
for every $\alpha\in[0,1], m > 0$ and $x\ge 0.$ This seems unnoticed even in the classical case $m=1.$
}
\end{rem}

Our next result is a uniform hyperbolic upper bound for the Kilbas-Saigo function $ E_{\alpha, m, m-\frac{1}{\alpha}},$ with an optimal power exponent by Remark \ref{rF2} (c) and an optimal constant since
$$1\,-\, E_{\alpha, m, m-\frac{1}{\alpha}}(-x)\; \sim\; \frac{\Ga(\alpha m)\, x}{\Ga(\alpha(m+1))}\quad \mbox{as $x\to 0.$}$$
\begin{prop}
\label{KSO2}
For every $\alpha\in(0,1], m >0$ and $x\ge 0,$ one has
$$E_{\alpha,m,m-\frac{1}{\alpha}}(-x)\; \le\; \frac{1}{\left( 1 + \frac{\Ga(1+\alpha m)}{\Ga(1+\alpha (m+1))}\, x\right)^{1+\frac{1}{m}}}\cdot$$
\end{prop}

\proof

The inequality is derived by convex ordering as in Proposition \ref{KSO1}: one has
\begin{eqnarray*}
\int_0^\infty \left( 1+ \sigma^{(\a)}_t\right)_+^{-\rho -\alpha} \, dt & \stackrel{d}{=} & \frac{\Ga(\rho)}{\Ga(\rho +\alpha)}\; {\bf T}(1 +\alpha\rho^{-1},\rho^{-1}, (1-\alpha)\rho^{-1})\\
& \prec_{cx} & \frac{\Ga(\rho)}{\Ga(\rho +\alpha)}\; {\bf T}(1 +\alpha\rho^{-1},\rho^{-1}, \rho^{-1})\; \stackrel{d}{=}\; \frac{\Ga(\rho +1)}{\Ga(\rho +1+\alpha)}\;{\bf \Gamma}_{1+\frac{\alpha}{\rho}}
\end{eqnarray*}
where the first identity follows from Corollary 3 in \cite{LS} as in the proof of Theorem \ref{Weib}, the convex ordering from Lemma \ref{Tcx} and the second identity from (2.7) in \cite{LS}. Then, using (\ref{KS2}) with $\rho = \alpha m,$ we get the required inequality.
\endproof

As in Proposition \ref{KSO1} we believe that there is also a uniform lower bound, with a more complicated optimal constant which can be read off from the asymptotic behaviour of the density at zero obtained in Proposition \ref{MASF}:
\begin{equation}
\label{ConjLow}
E_{\alpha,m,m-\frac{1}{\alpha}}(-x) \;\ge\; \frac{1}{(1 +(\alpha m)^{-\frac{\alpha}{m+1}}(\Ga(1+\alpha)\, G(1-\alpha; \alpha m)\, G(1+\alpha; \alpha m))^{-\frac{m}{m+1}}\, x)^{1+\frac{1}{m}}}\cdot
\end{equation}
Unfortunately, the proof of this general inequality still eludes us. The monotonicity property observed in Proposition \ref{Mono} does not help here, giving only the trivial lower bound zero. The discrete factorizations which are used in \cite{TS2014} are also more difficult to handle in this context, because the Mellin transform underlying $E_{\alpha,m,m-\frac{1}{\alpha}}$ is expressed in terms of generalized Pochhammer symbols. In the case $m=1,$ we could however get a proof of (\ref{ConjLow}).

\begin{prop}
\label{KSO3}
For every $\alpha\in(0,1)$ and $x\ge 0,$ one has
$$E_{\alpha,1,1-\frac{1}{\alpha}}(-x)\; \ge\; \frac{1}{\left( 1 +\sqrt{\frac{\Ga(1-\alpha)}{\Ga(1+\alpha)}}\, x\right)^{\! 2}}\cdot$$
\end{prop}

\proof By Remark \ref{rF1} (c) we have $E_{\alpha,1,1-\frac{1}{\alpha}}(-x) = \esp \left[ e^{-x T_\alpha^{(1)}}\right]$ and since
$$\esp \left[ e^{-x \,{\bf \Gamma}_2}\right]\; =\; \frac{1}{\left( 1 + x\right)^2}$$
for every $x \ge 0,$ as in the proof of Theorem 4 in \cite{TS2014} it is enough to show that
\begin{equation}
\label{StoO}
T_\alpha^{(1)}\,\prec_{st}\, \sqrt{\frac{\Ga(1-\alpha)}{\Ga(1+\alpha)}}\, {\bf \Gamma}_2,
\end{equation}
where $\prec_{st}$ stands for the usual stochastic order between two real random variables. Recall that $X\prec_{st} Y$ means ${\mathbb{P}}[X\ge x]\le{\mathbb{P}}[Y\ge x]$ for every $x\in{\mathbb{R}}.$ Since $T_{1/2}\stackrel{d}{=} 2\sqrt{{\bf \Gamma}_{1/2}},$ the case $\alpha =1/2$ is explicit and the stochastic ordering can be obtained directly. More precisely, the densities of both random variables in (\ref{StoO}) are respectively given by
$$\frac{x}{2}\; e^{-x^2/4}\qquad\mbox{and}\qquad \frac{x}{2}\; e^{-x/\sqrt{2}}$$
on $(0,\infty),$ where they cross only once at $x= 2\sqrt{2}.$ It is a well-known and easy result that this single intersection property yields (\ref{StoO}) - see Theorem 1.A.12 in \cite{SSh}.

The argument for the case $\alpha\neq 1/2$ is somehow analogous, but the details are more elaborate because the density of $T_\alpha^{(1)}$ is not explicit anymore. We proceed as in Theorem C of \cite{TS2014} and first consider the case where $\alpha$ is rational. Setting $\alpha = p/n$ with $n > p$ positive integers and $X_\alpha = T_\alpha^{(1)} = (Z_\alpha^{-\alpha})^{(1)}$ we have, on the one hand,
\begin{eqnarray*}
\esp[(X_{\alpha})^{n s}] & = & \frac{\esp[(T_\alpha)^{1+ns}]}{\esp[T_\alpha]} \\
& = & \frac{\Ga(2+ns) \Ga(1+pn^{-1})}{\Ga(1+pn^{-1} +ps)}\\
& = & \frac{n^{ns}}{p^{ps}}\,\times \,\esp\left[ \left({\bf B}_{\frac{2}{n}, \frac{1}{p} -\frac{1}{n}}\right)^s\right]\,\times\,\frac{\prod_{i=3}^{n+1} (in^{-1})_s}{\prod_{j=2}^{p} (jp^{-1} +n^{-1})_s}
\end{eqnarray*}
for every $s > -2n^{-1},$ where in the third equality we have used repeatedly the Legendre-Gauss multiplication formula for the Gamma function - see e.g. Theorem 1.5.2 in \cite{AAR1999}. The same formula implies, on the other hand,
\begin{eqnarray*}
\esp\left[\left(\sqrt{\frac{\Ga(1-\alpha)}{\Ga(1+\alpha)}}\; {\bf \Gamma}_2\right)^{\!\! ns}\, \right] & = & \frac{n^{ns}\, \kappa_\alpha^s}{p^{ps}}\,\times \,\esp\left[ \left({\bf \Gamma}_{\frac{2}{n}}\right)^s\right]\,\times\,\left( \prod_{i=3}^{n+1} (in^{-1})_s\right)\\
& = & \frac{n^{ns}}{p^{ps}}\,\times\,\esp\left[ \left(\kappa_\alpha\,\times \,{\bf \Gamma}_{\frac{2}{n}}\,\times\, \prod_{j=2}^{p}\; {\bf \Gamma}_{\frac{j}{p} +\frac{1}{n}}\right)^{\!s}\,\right]\,\times\, \frac{\prod_{i=3}^{n+1} (in^{-1})_s}{\prod_{j=2}^{p} (jp^{-1} +n^{-1})_s}
\end{eqnarray*}
for every $s > -2n^{-1},$ with the notation
$$\kappa_\alpha\; =\; \left(\prod_{i=1}^p\frac{\Ga(ip^{-1} - n^{-1})}{\Ga(ip^{-1} + n^{-1})}\right)^{\frac{n}{2}}.$$ Since
$$\frac{\prod_{i=3}^{n+1} (in^{-1})_s}{\prod_{j=2}^{p} (jp^{-1} +n^{-1})_s}\; =\;\esp\left[ \left(\prod_{i=2}^{p}\; {\bf B}_{\frac{i+1}{n}, \frac{i}{p} - \frac{i}{n}}\,\times\! \prod_{j=p+1}^{n} {\bf \Gamma}_{\frac{j+1}{n}}\right)^{\!s}\;\right]$$
for every $s > -3n^{-1},$ by factorization and Theorem 1.A.3(d) in \cite{TS2014} we are finally reduced to show
$${\bf B}_{\frac{2}{n}, \frac{1}{p} -\frac{1}{n}} \prec_{st}  \left(\prod_{i=1}^p\;\frac{\Ga(ip^{-1} - n^{-1})}{\Ga(ip^{-1} + n^{-1})}\right)^{\frac{n}{2}} \!\times \,{\bf \Gamma}_{\frac{2}{n}}\,\times\, \prod_{j=2}^{p}\; {\bf \Gamma}_{\frac{j}{p} +\frac{1}{n}}$$
for every $n >p$ positive integers. The latter is equivalent to
$$({\bf B}_{\frac{2}{n}, \frac{1}{p} -\frac{1}{n}})^{\frac{2}{n}} \prec_{st}  \left(\prod_{i=2}^p\;\frac{\Ga(ip^{-1} - n^{-1})}{\Ga(ip^{-1} + n^{-1})} \right)\times \left({\bf \Gamma}_{\frac{2}{n}}\,\times\, \prod_{j=2}^{p}\; {\bf \Gamma}_{\frac{j}{p} +\frac{1}{n}}\right)^{\frac{2}{n}}$$
and this is proved via the single intersection property exactly as for (5.1) in \cite{TS2014}: the random variable on the left-hand side has an increasing density on $(0,1),$ whereas the random variable on the right-hand side has a decreasing density on $(0,\infty),$ both densities having the same positive finite value at zero. We omit details. This completes the proof of (\ref{StoO}) when $\alpha$ is rational. The case when $\alpha$ is irrational follows then by a density argument.

\endproof

\begin{rem}
\label{rKS3}
{\em (a) It is easy to check from (\ref{Conca2}) and (\ref{Rhorho}) that
$$\frac{\Ga(1+\alpha)}{\Ga(1-\alpha)}\; =\;\alpha^\alpha\,\Ga(1+\alpha)\, G(1-\alpha; \alpha )\, G(1+\alpha; \alpha),$$
so that Proposition \ref{KSO3} leads to (\ref{ConjLow}) for $m=1,$ in accordance with the estimate (\ref{Fala}). In general, the absence of a tractable complement formula for the product $G(1-\alpha; \delta )\, G(1+\alpha; \delta)$ makes however the constant in (\ref{ConjLow}) more difficult to handle.

\medskip

(b) Combining Propositions \ref{KSO3} and \ref{KSO2} implies the following optimal bounds on the generalized Mittag-Leffler function $E_{\alpha,\alpha}(-x)$ for every $\alpha\in (0,1)$ and $x\ge 0,$ to be compared with (\ref{MLB}):
$$\frac{1}{\left( 1 +\sqrt{\frac{\Ga(1-\alpha)}{\Ga(1+\alpha)}}\, x\right)^2}\;\le\; \Ga(\alpha)\,E_{\alpha,\alpha}(-x)\; \le\; \frac{1}{\left( 1 +\frac{\Ga(1+\alpha)}{\Ga(1+ 2\alpha)}\, x\right)^2}\cdot$$
Notice that letting $\alpha\to 1$ leads to the trivial bound $0\le e^{-x}\le (2/(2+x))^2.$}
\end{rem}

Our last result in this paragraph gives analogous bounds for the generalized Mittag-Leffler functions $E_{\alpha,\beta}(-x)$ with $\alpha\neq\beta$ whenever they are completely monotone, that is for $\beta > \alpha$ - see Remark \ref{rF1} (c). Although there is no direct connection to fractional extreme laws, we include this result here because of its independent interest as a generalization of (\ref{MLB}).

\begin{prop}
\label{KSO4}
For every $\alpha\in(0,1], \beta > \alpha$ and $x\ge 0,$ one has the optimal bounds
$$\frac{1}{1 +\frac{\Ga(\beta-\alpha)}{\Ga(\beta)}\, x}\;\le\; \Ga(\beta)\,E_{\alpha,\beta}(-x)\; \le\; \frac{1}{1 +\frac{\Ga(\beta)}{\Ga(\beta+\alpha)}\, x}\cdot$$
\end{prop}

\proof By the last equality in Remark \ref{rF1} (c) we have
$$\Ga(\beta)\,E_{\alpha,\beta}(-x)\; =\; \esp\left[ e^{-x \, {\bf Y}_{\alpha,1,l}}\right]$$
with $l = (\beta - 1)/\alpha > 1-1/\alpha$ and ${\bf Y}_{\alpha, 1, l} \stackrel{d}{=} {\bf B}_{\alpha,\beta -\alpha}^\alpha \times\, T_\alpha^{(1)}.$ From Remark \ref{RCM} (a), one obtains
\begin{equation}
\label{MomXab}
\esp \left[ ({\bf Y}_{\alpha,1,l})^s\right]\; =\; \frac{\Ga(1+s) \Ga(\beta)}{\Ga(\beta + \alpha s)}
\end{equation}
for every $s > -1,$ which implies the factorization ${\bf L} \stackrel{d}{=} {\bf Y}_{\alpha,1,l}\times ({\bf \Gamma}_\beta)^{\alpha}.$ Since, by Jensen's inequality,
$$\frac{\Ga(\beta +\alpha)}{\Ga(\beta)}\; =\; \esp\left[ ({\bf \Gamma}_\beta)^\alpha\right] \prec_{cx} ({\bf \Gamma}_\beta)^{\alpha},$$
we deduce from Corollary 3.A.22 in \cite{SSh} the convex ordering
$${\bf Y}_{\alpha,1,l}\prec_{cx} \frac{\Ga(\beta)}{\Ga(\beta +\alpha)}\;{\bf L}$$
which, as above, implies
$$\Ga(\beta)\,E_{\alpha,\beta}(-x)\; \le\; \frac{1}{1 +\frac{\Ga(\beta)}{\Ga(\beta+\alpha)}\, x}$$
for every $x\ge 0.$

The argument for the other inequality is analogous to that of Proposition \ref{KSO3}. By density, we only need to consider the case $\alpha = p/n$ and $\beta = (p+q)/n$ with $p<n$ and $q$ positive integers. By (\ref{MomXab}) and the Legendre-Gauss multiplication formula, we obtain
$$\esp \left[ ({\bf Y}_{\alpha,1,l})^{ns}\right]\; =\;\frac{n^{ns}}{p^{ps}}\,\times \,\esp\left[ \left({\bf B}_{\frac{1}{n},\frac{q}{np}}\right)^s\right]\,\times\,\frac{\prod_{i=2}^{n} (in^{-1})_s}{\prod_{j=1}^{p-1} (jp^{-1} + (p+q)(np)^{-1})_s}$$
for every $s > -n^{-1}.$ On the other hand, one has
$$\esp\left[ \left( \frac{\Ga(\beta-\alpha)}{\Ga(\beta)}\;{\bf L}\right)^{\! ns}\,\right] \; = \; \frac{n^{ns}}{p^{ps}}\;\esp\left[ \left(\kappa_{\alpha,\beta}\,\times \,{\bf \Gamma}_{\frac{1}{n}}\,\times\, \prod_{j=1}^{p-1}\; {\bf \Gamma}_{\frac{j}{p} +\frac{p+q}{np}}\right)^{\!\! s}\,\right]\,\times\,\frac{\prod_{i=2}^{n} (in^{-1})_s}{\prod_{j=1}^{p-1} (jp^{-1} + (p+q)(np)^{-1})_s}$$
with
$$\kappa_{\alpha,\beta}\; =\; p^p \left(\frac{\Ga(qn^{-1})}{\Ga((p+q)n^{-1})} \right)^n.$$
Comparing these two formulas, we are reduced to show
$$({\bf B}_{\frac{1}{n},\frac{q}{np}})^{\frac{1}{n}}\,\prec_{st}\, p^{\frac{p}{n}} \left(\frac{\Ga(qn^{-1})}{\Ga((p+q)n^{-1})} \right)\times\left({\bf \Gamma}_{\frac{1}{n}}\,\times\, \prod_{j=1}^{p-1}\; {\bf \Gamma}_{\frac{j}{p} +\frac{p+q}{np}}\right)^{\frac{1}{n}}$$
for every $p< n$ and $q$ positive integers. This is obtained in the same way as above via the single intersection property. We leave the details to the reader.

\endproof

\subsection{Some properties related to infinite divisibility}

\label{SPID}

In this paragraph we derive some infinite divisibility properties for the fractional extreme distributions, in the spirit of the Corollary in \cite{JSW2018}. Recall that the law of a positive random variable $X$ is called a generalized Gamma convolution ($X\in\cG$ for short) if there exists a suitably integrable deterministic function $a : {\mathbb{R}}^+\to{\mathbb{R}}^+$ such that
$$X\; \stackrel{d}{=}\; \int_0^\infty a(t)\, d\Gamma_t$$
where $\{\Ga_t, \, t\ge 0\}$ is the Gamma subordinator. Equivalently, one has $X\in\cG$ iff its log-Laplace exponent reads
$$-\log \esp[e^{-\lambda X}]\; =\; a\lambda\, +\, \int_0^\infty (1- e^{-\lambda x})\, k(x) \frac{dx}{x}$$
with $a \ge 0$ and $k(x)$ a CM function. This representation shows that a random variable $X\in\cG$ is also infinitely divisible ($X\in{\mathcal{I}}$ for short). An important subclass of $\cG$ is that of hyperbolically completely monotone random variables, which we will denote by $\cH.$ By definition, one has $X\in\cH$ iff $X$ has a positive density $f_X$ on $(0,\infty)$ such that $f_X(uv) f_X(u v^{-1})$ is CM in the variable $v+v^{-1}$ for all $u >0.$ We refer to \cite{Bond} for a classic account on the classes $\cG$ and $\cH$, including the above facts and much more. See also \cite{JRY2008} for a more recent survey. In Chapter 7 of \cite{Bond}, the class $\cG$ is extended to distributions on the real line, under the denomination EGGC. More precisely, an infinitely divisible distribution on ${\mathbb{R}}$ is called an EGGC if its L\'evy measure has a density on ${\mathbb{R}}^*$ of the type $\vert x\vert^{-1} k(x)$ with $x\mapsto k(x)$ and $x\mapsto k(-x)$ being CM functions on $(0,\infty).$ In the following, we will say that such L\'evy measures belong to the Thorin class. We will also use the same notation $X\in\cG$ to denote EGGC distributions, since there shall be no ambiguity on the support of $X.$

Our analysis is based on the following two lemmas on infinite Beta products, which have an independent interest. The first one is a precision made on (2.5) in \cite{LS}, whereas the second one is an extension of the main argument for proof of the Corollary in \cite{JSW2018}. This extension was already discussed in Remark 1 therein but we give some detail for the sake of completeness, and for the independent interest of the logarithmic estimate (\ref{Tabc}).

\begin{lem}
\label{Pow}
For every $a,b,c > 0$ one has
$$\left( \frac{\Ga(a+c)}{\Ga(a)\,b^c}\; {\bf T}(ab^{-1},b^{-1},cb^{-1})\right)^{\!\frac{1}{b}}\,\stackrel{d}{=}\,\;\frac{\Ga((a+c)b^{-1})}{\Ga(ab^{-1})}\; {\bf T}(a,b,c)$$
\end{lem}

\proof
By (2.5) in \cite{LS}, it is enough to show that the random variables on both sides have the same expectation. By Proposition 2 in \cite{LS}, this amounts to
$$\frac{G((a+c+1)b^{-1}; b^{-1})\, G(ab^{-1}; b^{-1})}{G((a+c)b^{-1}; b^{-1})\, G((a+1)b^{-1}; b^{-1})}\; =\; \frac{\Ga((a+c)b^{-1})\,b^{\frac{c}{b}}}{\Ga(ab^{-1})}$$
which is a consequence of (\ref{Conca2}). Alternatively, the identity in law can be obtained from (\ref{Conca3}).

\endproof

\begin{lem}
\label{HCMT}
For every $a,b,c > 0$ one has ${\bf T}(a,b,c)^{-1}\,\in\,\cG$ and
$${\bf T}(a,b,c)\,\in\,\cH\;\Leftrightarrow\;{\bf T}(a,b,c)\,\in\,{\mathcal{I}}\;\Leftrightarrow\; c\,\ge\, b.$$
\end{lem}

\proof The fact that ${\bf T}(a,b,c)^{-1}\in\cG$ is derived as in Corollary 8 in \cite{LS}. For the if part of the equivalence, we first notice that the case $c=b$ is obvious since  by (2.7) in \cite{LS} one has ${\bf T}(a,b,b) \stackrel{d}{=} a^{-1}{\bf \Gamma}_a,$ whose density is HCM. If $c > b,$ we decompose
$${\bf T}(a,b,c)\; \stackrel{d}{=}\; {\bf T}(a,b,b)\,\times\,{\bf T}(a+b, b, c-b)\; \stackrel{d}{=}\; a^{-1}{\bf \Gamma}_a\,\times\, {\bf T}(a+b, b, c-b),$$
where the first identity in law follows from (2.4) in \cite{LS}, and we can conclude exactly as in the proof of the Corollary in \cite{JSW2018}, since $a+b >a.$ For the only if part of the equivalence, we need to show $c < b\Rightarrow {\bf T}(a,b,c)\not\in{\mathcal{I}}.$ To do so, we first deduce from (\ref{Stirling}) and some simplifications the limit behaviour
$$\frac{1}{s} \left( \esp\left[ \left( {\bf T}(a,b,c)\right)^{\frac{bs}{c}}\right]\right)^{\frac{1}{s}}\;\to \; \frac{1}{c \, {\rm e}}\,\left( \frac{\Ga(ab^{-1})}{\Ga((a+c)b^{-1})}\right)^{\frac{b}{c}}\qquad \mbox{as $s\to\infty.$}$$
Applying, as for the Corollary in \cite{JSW2018}, Lemma 3.2 in \cite{CSY1999}, implies
\begin{equation}
\label{Tabc}
\log {\mathbb{P}} \left[ {\bf T}(a,b,c) > x\right]\;\sim\; -\,c\, \left( \frac{\Ga((a+c)b^{-1})}{\Ga(ab^{-1})}\, x\right)^{\frac{b}{c}}\qquad \mbox{as $x\to\infty.$}
\end{equation}
This shows that ${\bf T}(a,b,c)\not\in{\mathcal{I}}$ whenever $b >c,$ since its upper tail probabilities are then superexponentially small - see e.g. Theorem 26.8 in \cite{S1999}.

\endproof

\begin{rem}
\label{RTabc}
{\em For $a=b,$ it follows after some simplifications from the Theorem and the Proposition in \cite{JSW2018}, with $\alpha = c$ and $m = a+c$ therein, that the density of ${\bf T}(a,a,c)$ is equivalent to
$$\left(\frac{(2\pi)^{\frac{a+c}{2} -1}\, a^{\frac{c}{a}}\,(\Ga(1+ca^{-1}))^{\frac{(a+c)(1+c)}{2c} -1}}{\sqrt{ac}\,G(a+c,a)}\right)\, x^{\frac{(a+c)(1+c)}{2c} -2}\, e^{-c\,\left(\Ga(1+ca^{-1})\, x\right)^{\frac{a}{c}}}\qquad \mbox{as $x\to\infty.$}$$
This is an improvement at the natural scale of the logarithmic estimate (\ref{Tabc}). Notice that in the case $a=b=c,$ this behaviour matches the exact formula $\frac{a^{a} x^{a-1}}{\Ga(a)}\, e^{-ax}$ for the density of ${\bf T}(a,a,a)\stackrel{d}{=} a^{-1}{\bf \Gamma}_a.$ The exact behaviour of the density of ${\bf T}(a,b,c)$ at infinity for all $a,b,c > 0$ is an interesting open question.}
\end{rem}

We can now state the main result of this paragraph.

\begin{prop}
\label{HCMWFG}
For every $\lambda, \rho > 0$ and $\alpha\in[0,1]$ one has

\begin{itemize}

\item ${\W}_{\a,\lbd,\rho}\in\cG$ if $\rho\le 1$ and ${\W}_{\a,\lbd,\rho}\in\cH$ if $\rho\le 1-\alpha,$

\item ${\F}_{\a,\lambda,\rho}\in\cH$ if $\rho\le 1-\alpha,$

\item $\Gal\in\cG.$
\end{itemize}

\end{prop}

\proof

It follows from Theorem \ref{Weib}, Theorem \ref{Fresh}, Lemma \ref{Pow} and Formula (2.7) in \cite{LS} that
$${\W}_{\a,\lbd,\rho}\;\stackrel{d}{=}\;\left(\frac{\rho^\alpha}{\lambda}\right)^{\frac{1}{\rho}} \frac{\Ga(1+\rho^{-1})}{\Ga(1+(1-\alpha)\rho^{-1})}\,\times\, \frac{{\bf T}(\rho,\rho,1)}{{\bf T}(\rho,\rho,1-\alpha)}$$and
$${\F}_{\a,\lambda,\rho}\;\stackrel{d}{=}\;\left(\frac{\lambda}{\rho^\alpha}\right)^{\frac{1}{\rho}} \frac{1}{\Ga(1+\alpha\rho^{-1})}\,\times\, \frac{{\bf T}(\rho+\alpha,\rho,1-\alpha)}{{\bf T}(\rho,\rho,1)}\cdot$$
If $\rho\le 1-\alpha,$ by the second statement in Lemma \ref{HCMT} all random variables involved on the right-hand side belong to $\cH$, and Theorem 5.1.1 in \cite{Bond} implies that ${\W}_{\a,\lbd,\rho}$ and ${\F}_{\a,\lambda,\rho}$ belong to $\cH$ as well. The fact that ${\W}_{\a,\lbd,\rho}\in\cG$ for $\rho\le 1$ follows from the first statement in Lemma \ref{HCMT} and Theorem 6.2.1 in \cite{Bond}. Last, a consequence of (\ref{Malm1}) is
\begin{eqnarray*}
\log \esp[e^{ s\Gal)}] & = & \log \Ga(1+ s\lambda^{-1})\; -\; (1-\alpha)\, \log \Ga(1+ s\lambda^{-1}) \\
& = & -\alpha\gamma\lambda^{-1} s\, +\, \int_{\mathbb{R}} (e^{sx} - 1 - sx)\,({\bf 1}_{\{x <0\}} \, +\, (1-\alpha){\bf 1}_{\{ x > 0\}})\, \frac{dx}{\vert x\vert (e^{\lambda \vert x\vert} -1)}
\end{eqnarray*}
for every $s >-\lambda.$ By the L\'evy-Khintchine formula, this shows that $\Gal$ is infinitely divisible and that its L\'evy measure belongs to the Thorin class, in other words $\Gal\in\cG.$

\endproof

\begin{rem}
\label{RHCM}
{\em (a) The same proof shows that ${\bf G}_\alpha\in\cG,$ since
$$\log \esp[e^{ s{\bf G}_\alpha)}] \; = \; (1-\alpha)\, \log \Ga(1+ s) \; = \; (1-\alpha)\left( -\gamma s\, +\, \int_{-\infty}^0 (e^{sx} - 1 - sx)\,\frac{dx}{\vert x\vert (e^{\vert x\vert} -1)}\right)$$
for every $s > -1,$ so that ${\bf G}_\alpha\in{\mathcal{I}}$ with a L\'evy measure in the Thorin class. On the other hand, Theorem 2.1 in \cite{BL2015} shows a superexponential behaviour for the density of ${\bf L}_\alpha$ at infinity, and this implies as above that ${\bf L}_\alpha\not\in{\mathcal{I}}$ except for $\alpha = 0$ or $\alpha = 1.$

\medskip

(b) The above proof also shows that in general, ${\W}_{\a,\lbd,\rho}$ and ${\F}_{\a,\lambda,\rho}$ can be expressed as the independent quotient of two random variables in $\cG.$ Unfortunately, this does not allow one to infer further infinite divisibility properties. We believe however that ${\F}_{\a,\lambda,\rho}\in{\mathcal{I}}$ for all $\rho > 0.$ See \cite{BS2013} for a panorama of results related to the infinite divisibility of the classical extreme distributions.}
\end{rem}

\subsection{Some complements on the Le Roy function}

\label{SCLR}

In this paragraph we derive some miscellaneous results on the Le Roy function
$$\cL_\alpha(x) \; =\; \sum_{n\ge 0} \frac{x^n}{(n!)^\alpha}\cdot$$
This function played a role in the proof of Theorem \ref{Gumb} and can be viewed as another generalization of the exponential function, which it is interesting to compare to the classical Mittag-Leffler function $E_\alpha(x).$ Here and throughout we discard the explicit cases $\cL_0(x) = E_0(x) = 1/(1-x)$ and $\cL_1(x) = E_1(x) = e^x.$\\

We begin with the asymptotic behaviour at infinity. Le Roy's original result - see \cite{L1899} p. 263 - reads
$$\cL_\alpha (x)\; \sim\; \frac{(2\pi)^{\frac{1-\alpha}{2}}}{\sqrt{\alpha}}\, x^{\frac{1-\alpha}{2}}\, e^{\alpha x^{\frac{1}{\alpha}}}\qquad \mbox{as $x\to\infty,$}$$
and is obtained by a variation on Laplace's method. The latter method can be used to solve Exercise 8.8.4 in \cite{Olver}, which states
\begin{equation}
\label{Olv1}
\cL_\alpha (-x)\; = \; \frac{2 (2\pi)^{\frac{1-\alpha}{2}}}{\sqrt{\alpha}}\,  x^{\frac{1-\alpha}{2\alpha}}\, e^{\alpha \cos(\pi/\alpha) x^{\frac{1}{\alpha}}}\, \left( \sin\left( \pi (2\alpha)^{-1} + \alpha \sin\left(\pi\alpha^{-1}\right) x^{\frac{1}{\alpha}}\right)\, +\, O(x^{-\frac{1}{\alpha}})\right)
\end{equation}
for $\alpha\ge 2$ and
\begin{equation}
\label{Olv2}
\cL_\alpha (-x)\; \sim \; \frac{1}{\alpha^\alpha\, \Ga(1-\alpha) \, x\, (\log x)^\alpha}
\end{equation}
for $\alpha \in (1,2),$ as $x\to \infty.$ The following estimate, which seems to have passed unnoticed in the literature, completes the picture.

\begin{prop}
\label{Olv3}
For every $\alpha\in (0,1),$ one has
$$\cL_\alpha (-x)\; \sim \; \frac{1}{\Ga(1-\alpha) \, x\, (\log x)^\alpha}\qquad \mbox{as $x\to\infty.$}$$
\end{prop}

\proof

By (\ref{LREM}), we have
$$\cL_\alpha (-x)\; =\; {\mathbb{P}}\left[ {\bf L} > x {\bf L}_\alpha\right]\; =\; \int_0^\infty e^{-xt} \, f_\alpha(t)\, dt$$
with ${\bf L}_\alpha = e^{{\bf G}_\alpha}$ having density $f_\alpha$ on $(0,\infty).$ On the one hand, recalling
$$\esp \left[ {\bf L}_\alpha^s\right]\; =\; \Ga(1+s)^{1-\alpha}$$
for every $s > -1,$ we have $f_\alpha = e_{1-\alpha}$ with the notation of \cite{BL2015} and we can apply Theorem 2.4 therein to obtain
\begin{equation}
\label{Asfa}
f_\alpha(x)\;\sim\; \frac{1}{\Ga(1-\alpha)\, (-\log x)^\alpha}\qquad \mbox{as $x\to 0.$}
\end{equation}
Plugging this estimate into the above expression for $\cL_\alpha (-x),$ we conclude the proof by a direct integration.
\endproof

\begin{rem}
\label{Rolv}
{\em (a) The estimate (\ref{Asfa}) also gives the asymptotic behaviour of the density of $\Gal$ at the right end of the support. Indeed, by multiplicative convolution the density of $e^{\lambda \Gal}$ on $(0,\infty)$ writes
$$\int_0^\infty e^{-xy} \; y f_\alpha(y)\, dy \; \sim\; \frac{1}{\Ga(1-\alpha)\, x^2\, (\log x)^\alpha}\qquad\mbox{as $x\to\infty,$}$$
where the estimate follows from (\ref{Asfa}) as in the above proof. A change of variable implies then
$$f^{\bf G}_{\alpha,\lambda}(x)\;\sim\;\left(\frac{\lambda^{1-\alpha}}{\Ga(1-\alpha)}\right)\, x^{-\alpha}\, e^{-\lambda x}\qquad\mbox{as $x\to\infty.$}$$
The asymptotic behaviour of the density at the left end can be obtained as in Paragraph \ref{ABTD} via the moment generating function
$$\esp\left[ e^{s\Gal}\right]\; =\; \Ga(1+ s\lambda^{-1})\,\Ga(1-s\lambda^{-1})^{1-\alpha}, \qquad \vert s\vert <\lambda.$$
Reasoning as in Proposition \ref{MASW} via the converse mapping theorem leads to
$$f^{\bf G}_{\alpha,\lambda}(-x)\;\sim\;\lambda\, e^{-\lambda x}\qquad\mbox{as $x\to\infty,$}$$
in accordance with the first term in the expansion given in Corollary \ref{GumbD}. Observe that this converse mapping argument does not work directly for estimating $f^{\bf G}_{\alpha,\lambda}(x)$ at the right end, because of the fractional singularity in the moment generating function.

\medskip

(b) In the case $\alpha =2,$ one has $\cL_2(x) = I_0(2\sqrt{x})$ and $\cL_2(-x) = J_0(2\sqrt{x})$ for all $x\ge 0,$ where $I_0$ and $J_0$ are the classical, modified or not, Bessel functions with index 0. In particular, a full asymptotic expansion for $\cL_2$ at both ends of the support is available, to be deduced e.g. from (4.8.5) and (4.12.7) in \cite{AAR1999}. These expansions also exist when $\alpha$ is an integer since $\cL_\alpha$ is then a generalized Wright function - see Chapter F.2.3 in \cite{GKMR} and the original articles by Wright quoted therein. The case when $\alpha$ is not an integer is open, and might be technical in the absence of a true Mellin-Barnes representation.}
\end{rem}

Our next result characterizes the connection between the entire function $\cL_\alpha(z)$ and random variables. Recall that a function $f : \CC\to\CC$ which is holomorphic in a neighbourhood $\Omega$ of the origin is a moment generating function (MGF) if there exists a real random variable $X$ such that
$$f(z) \; =\; \esp\left[ e^{z X}\right],\qquad z\in\Omega.$$
In particular, $\cL_0$ is the MGF of the exponential law ${\bf L}$ and $\cL_1$ is that of the constant variable ${\bf 1}.$

\begin{prop}
\label{LRCM}
The function $\cL_\alpha(z)$ is the {\em MGF} of a real random variable if and only if $\alpha\le 1.$ In this case, one has
$$\cL_\alpha(z)\; =\; \esp\left[ e^{z {\bf L}_\alpha}\right], \qquad z\in\CC.$$
\end{prop}

\proof

The if part is a consequence of (\ref{LREM}) as in the proof of Proposition \ref{Olv3}. For the only if part, the estimates (\ref{Olv1}) and (\ref{Olv2}) show that $\cL_\alpha(z)$ takes negative values on ${\mathbb{R}}^-$ when $\alpha  > 1,$ so that it cannot be the moment generating function of a real random variable.

\endproof

Observe that since ${\bf L}_\alpha$ is non-negative, the above result also shows $\cL_\alpha(-x)$ is CM on $(0,\infty)$ if and only if $\alpha\le 1,$ echoing Pollard's classical result for the Mittag-Leffler $E_\alpha(-x)$ - see Proposition 3.23 in \cite{GKMR}. One can ask if there are further complete monotonicity properties for $\cL_\alpha,$ as in \cite{TS2015} for $E_\alpha.$ In a different direction, the following result gives a monotonicity property which is akin to Proposition \ref{Mono}.

\begin{prop}
\label{LRMono}
The mapping $\alpha\mapsto\cL_\alpha (x)$ decreases on $[0,1]$ for every $x\in{\mathbb{R}}.$
\end{prop}

\proof The fact that $\alpha\mapsto\cL_\alpha (x)$ decreases on ${\mathbb{R}}^+$ is obvious for $x\ge 0,$ by the definition of $\cL_\alpha.$ To show the property on $[0,1]$ for $x <0,$ we will use a convex ordering argument. More precisely, the Malmsten formula (\ref{Malm1}) and the L\'evy-Khintchine formula show that for every $t \in [0,1],$ the random variable ${\bf G}_{1-t} = \log {\bf L}_{1-t}$ is the marginal at time $t$ of a real L\'evy process, since $\esp[e^{\ii z {\bf G}_{1-t}}] =  \Ga(1+ \ii z)^t = e^{t \psi(z)}$ for every $z\in{\mathbb{R}},$ with
$$\psi(z)\; =\; -\gamma\ii z\, +\, \int^0_{-\infty} (e^{\ii zx} - 1 - \ii zx)\, \frac{dx}{\vert x\vert (e^{\vert x\vert} -1)}\cdot$$
This is actually well-known - see again Example E in \cite{CPY}. By independence and stationarity of the increments of a L\'evy process, we deduce that there exists a multiplicative martingale $\{M_t, t\in [0,1]\}$ such that $M_t \stackrel{d}{=} {\bf L}_{1-t}$ for every $t \in [0,1]$ and Jensen's inequality implies
$${\bf L}_\beta\;\prec_{cx}\; {\bf L}_\alpha$$
for every $0\le\alpha\le\beta\le 1.$ Applying the definition of convex ordering to the function $\varphi(x) = e^{-x},$ we get $\cL_\beta(-x)\le\cL_\alpha(-x)$ for every $x > 0$ and $0\le\alpha\le\beta\le 1,$ as required.
\endproof

\begin{rem}
\label{Picq}

{\em (a) In the terminology of \cite{HPRY}, the family $\{{\bf L}_{1-\alpha}, \; \alpha\in[0,1]\}$ is a peacock, whose associated multiplicative martingale is here completely explicit. We refer to \cite{HPRY} for numerous examples of explicit peacocks related to exponential functionals of L\'evy processes. Observe from Lemma \ref{Tcx} that the family $\{{\bf T}(a,b,t), \; t >0\}$ is also a peacock.

\medskip

(b) It is easily seen from Corollary 3 in \cite{LS} that ${\bf T}(1+\alpha q^{-1},q^{-1}, (1-\alpha) q^{-1})\stackrel{d}{\longrightarrow} {\bf L}_\alpha$ as $q\to\infty.$ On the other hand, we could not prove that the family $\{{\bf T}(1+(1-\alpha) q^{-1},q^{-1}, \alpha q^{-1}), \; \alpha\in[0,1]\}$ is a peacock for any fixed $q >0.$}
\end{rem}

Letting $\alpha\to 0,1$  in the above proposition leads to the bounds
$$e^{-x}\;\le\;\cL_\beta(-x)\;\le\;\cL_\alpha(-x)\;\le\; \frac{1}{1+x}$$
for every $x \ge 0$ and  $0\le\alpha\le\beta\le 1,$ to be compared with the less complete bounds (6.9) in \cite{TS2014} for $E_\alpha(-x).$ The hyperbolic upper bound is optimal as in Propositions \ref{KSO1} and \ref{KSO2}, since $1- \cL_\alpha(-x)\sim x$ as $x\to 0.$ The exponential lower bound is thinner than the order given in the estimate (\ref{Asfa}). On the other hand, it does not seem that stochastic ordering arguments can help to get a uniform estimate involving a logarithmic term. The following last proposition gives alternative bounds on $\cL_\alpha(-x)$ in terms of Kilbas-Saigo functions. It is a direct consequence of the Bernstein representations given in Remark \ref{rW1} (d), Remark \ref{rF1} (c) and Proposition \ref{LRCM}, and of the monotone character of the function $s\mapsto (1-xs^{-1})_+^s$ on $(0,\infty)$ for every $x\in{\mathbb{R}}.$ We omit details.

\begin{prop}
\label{LRKS}
For every $\alpha\in [0,1]$ and $x, m > 0,$ one has
$$\left\{ \begin{array}{l}
E_{\alpha, m,m-\frac{1}{\alpha}}(-(\alpha (m+1))^\alpha x)\;\le\;\cL_\alpha(-x)\;\le\; E_{\alpha, m+1,m}(-(\alpha m)^\alpha x),\\
E_{\alpha, m+1,m}((\alpha m)^\alpha x)\;\le\;\cL_\alpha(x)\;\le\; E_{\alpha, m,m-\frac{1}{\alpha}}((\alpha (m+1))^\alpha x).\end{array}\right.$$
Besides, $E_{\alpha, m,m-\frac{1}{\alpha}}((\alpha (m+1))^\alpha x)$ and $E_{\alpha, m+1,m}((\alpha m)^\alpha x)\to\cL_\alpha(x)$ as $m\to\infty,$ for every $x\in{\mathbb{R}}.$
\end{prop}

\medskip

\section*{Appendix}

\subsection*{A.1. Fractional integrals and derivatives} In this paragraph, we fix the notation on the fractional operators which are used throughout the paper. This is an excerpt from the beginning of Chapter 2 in \cite{KST}. We will consider only three kinds of such operators which are the most familiar ones, and our fractional parameter $\alpha$ will always be supposed in $[0,1].$ There are certainly many other fractional operators with a larger family of fractional parameters, and we refer to the whole Chapter 2 in \cite{KST} for an account.

\subsubsection*{{\em A.1.1.} Progressive Liouville operators on the half-axis} For every $\alpha\in (0,1),$ the operator $f\mapsto \Iap(f)$ with
$$\Iap (f)(x)\; =\; \frac{1}{\Ga (\alpha)}  \int_{0}^{x} (x-u)^{\alpha-1}f(u) \, du,\quad  x > 0,$$
is well-defined, taking possibly infinite values, on measurable functions $f : (0,\infty)\to {\mathbb{R}}^+.$ It is easy to see that if $f$ is integrable at zero, then so is $\Iap(f)$ - see Lemma 2.1 in \cite{KST} for a more general result. The corresponding fractional derivative $f\mapsto \Dap f$ is such that
$$\Dap (f) (x)\; =\; \frac{\mathrm{d}}{\mathrm{d} x}\left( {\rm I}^{1-\alpha}_{0+} (f)\right) (x)$$
and is well-defined almost everywhere as soon as $f = \Iap(g)$ for some $g$ integrable at zero. Moreover, for such functions there is an inversion formula
 \begin{equation}
\label{Fr0+}
\Iap\left(\Dap(f)\right)\; =\; \Iap(g)\; =\; f,
\end{equation}
which is valid almost everywhere - see Lemma 2.5 in \cite{KST}. These operators are extended to the boundary cases $\alpha = 0$ with ${\rm I}^0_{0+} = {\rm D}^0_{0+} = {\rm Id}$ and $\alpha = 1$ with  ${\rm I}^1_{0+}$ and ${\rm D}^1_{0+}$ being respectively the usual running integral and derivative - see (2.1.7) in \cite{KST}.

\subsubsection*{{\em A.1.2.} Regressive Liouville operators on the half-axis}
 For every $\alpha\in (0,1),$ the operator $f\mapsto \Im(f)$ with
$$\Im (f)(x)\; =\; \frac{1}{\Ga (\alpha)}  \int_{x}^\infty (u-x)^{\alpha-1}f(u) \, du,\quad  x > 0,$$
is well-defined on measurable functions $f : (0,\infty)\to {\mathbb{R}}^+.$ It is easy to see that if $f$ is integrable at infinity, then so is $\Im(f)$ - see again Lemma 2.1 in \cite{KST}. The corresponding fractional derivative $f\mapsto \Dm f$ is such that
$$\Dm (f) (x)\; =\; -\frac{\mathrm{d}}{\mathrm{d} x}\left( {\rm I}^{1-\alpha}_{-} (f)\right) (x)$$
and is well-defined as soon as $f = \Im(g)$ for some $g$ integrable at infinity. Moreover, for such functions there is an inversion formula
 \begin{equation}
\label{Fr0-}
\Im\left(\Dm(f)\right)\; =\; \Im(g)\; =\; f,
\end{equation}
which is valid almost everywhere. These operators are extended to the boundary cases $\alpha = 0$ with ${\rm I}^0_{-} = {\rm D}^0_{-} = {\rm Id}$ and $\alpha = 1$ with  ${\rm I}^1_{-}$ and ${\rm D}^1_{-}$ being respectively the usual running integral and the opposite of the usual derivative - see again (2.1.7) in \cite{KST}.

\subsubsection*{{\em A.1.3.} Progressive Liouville operators on the real axis}  For every $\alpha\in (0,1),$ the operator $f\mapsto \Ip(f)$ with
$$\Ip (f)(x)\; =\; \frac{1}{\Ga (\alpha)}  \int_{-\infty}^x (x-u)^{\alpha-1} f(u) \, du,\quad  x \in{\mathbb{R}},$$
is well-defined on measurable functions $f : {\mathbb{R}} \to {\mathbb{R}}^+.$ Observe that $\Ip(f) (-x) = \Im (g) (x)$ with $g(x) = f(-x)$ so that we can transfer to $\Ip$ the properties on $\Im.$ In particular, if $f$ is integrable at $-\infty$, then so is $\Ip(f).$ The corresponding fractional derivative $f\mapsto \Dp f$ is such that
$$\Dp (f) (x)\; =\; -\frac{\mathrm{d}}{\mathrm{d} x}\left( {\rm I}^{1-\alpha}_{+} (f)\right) (x)$$
and is well-defined as soon as $f = \Ip(g)$ for some $g$ integrable at $-\infty$. Moreover, for such functions there is an inversion formula $\Ip\left(\Dp(f)\right) = \Ip(g) = f,$ which is valid almost everywhere. These operators are extended to the boundary cases $\alpha = 0$ with ${\rm I}^0_{+} = {\rm D}^0_{+} = {\rm Id}$ and $\alpha = 1$ with  ${\rm I}^1_{+}$ and ${\rm D}^1_{+}$ being respectively the usual running integral and derivative.

\subsection*{A.2. Barnes' double Gamma function} In this paragraph, mostly taken from \cite{BK1997} and \cite{K2011} to which we refer for further results, we gather some useful facts about Barnes' double Gamma function $G(z;\delta).$ For every $\delta > 0,$ this function is defined as the unique solution to the functional equation
\begin{equation}
\label{Conca1}
G(z+1;\delta)\; =\; \Gamma(z\delta^{-1}) G(z;\delta)
\end{equation}
with normalization $G(1;\delta) =1.$ This function is holomorphic on $\CC$ and admits the following Malmsten type representation
\begin{equation}
\label{Malm}
G(z;\delta)\; =\; \exp\int_0^\infty \left( \frac{1-e^{-zx}}{(1-e^{-x})(1-e^{-\delta x})} \, -\, \frac{ze^{-\delta x}}{1-e^{-\delta x}}\, +\, (z-1)(\frac{z}{2\delta} -1)e^{-\delta x}\, - \, 1\right)\frac{dx}{x}
\end{equation}
which is valid for $\Re(z) >0$ - see (5.1) in \cite{BK1997}. Putting (\ref{Conca1}) and (\ref{Malm}) together and making some simplifications, we recover the standard Malmsten formula for the Gamma function
\begin{equation}
\label{Malm1}
\Ga(1+z)\; =\; \exp\left\{ -\gamma z\, +\, \int^0_{-\infty} (e^{zx} - 1 - zx)\, \frac{dx}{\vert x\vert (e^{\vert x\vert} -1)}\right\}
\end{equation}
for every $z >-1,$ where $\gamma$ is Euler's constant. The following Stirling type asymptotic behaviour
\begin{equation}
\label{Stirling}
\log G(z;\delta)\; -\; \frac{1}{2 \delta} \left( z^2 \log z\, -\, (\frac{3}{2} + \log \delta)\, z^2 \, -\,(1+\delta)\, z\log z\right)\; -\; A \, z\; -\; B \log z\; \to\; C
\end{equation}
is valid for $\vert z\vert\to\infty$ with $\vert \arg (z) \vert < \pi,$ for some real constants $A, B$ and $C$ which are given in (4.5) of \cite{BK1997}. There is a second concatenation formula
\begin{equation}
\label{Conca2}
G(z+\delta;\delta)\; =\; (2\pi)^{(\delta-1)/2} \delta^{1/2 -z}\Gamma(z) G(z;\delta)
\end{equation}
which is valid for all $z\in\CC,$ the right-hand side being understood as an anxiolytic extension when $z$ is a non-positive integer - see (4.6) in \cite{K2011} and the references therein. Observe that (\ref{Conca1}) and (\ref{Conca2}) lead readily to the closed formula
\begin{equation}
\label{Rhorho}
G(\delta;\delta)\; =\; G(1+\delta;\delta)\; =\; (2\pi)^{(\delta-1)/2} \delta^{-1/2}.
\end{equation}
In this paper we make an extensive use of the following Pochhammer type symbol
\begin{equation}
\label{Poch}
[a;\delta]_s\; =\; \frac{G(a +s;\delta)}{G(a;\delta)}
\end{equation}
which is well-defined for every $a,\delta > 0$ and $s > - a.$ The following formula
\begin{equation}
\label{Conca3}
[a\delta^{-1};\delta^{-1}]_{s\delta^{-1}} \; =\; (2\pi)^{s(1/\delta-1)/2}\, \delta^{s^2/2\delta - s(1+(1-2a)/\delta)/2}\,[a;\delta]_s
\end{equation}
can be deduced from (4.10) in \cite{K2011} - beware the different normalization for $G(1;\delta)$ therein which becomes irrelevant when considering the Pochhammer type symbol. Notice also that (\ref{Conca2}) yields
\begin{equation}
\label{Conca4}
\delta^s\,[a+\delta;\delta]_s\; =\; (a)_s \,
[a;\delta]_s
\end{equation}
with the standard notation
$$(a)_s \; =\; \frac{\Ga(a+s)}{\Ga(a)}$$
for the usual Pochhammer symbol. Finally, we observe from the double product representation of $G(z,\delta)$ - see e.g. (4.4) in \cite{K2011}, that for every $a, \delta > 0$ one has
\begin{equation}
\label{zeroes}
\inf\{s > 0, \;  [a;\delta]_{-s} = 0\}\; =\; a
\end{equation}
and that this zero is simple and isolated on the complex plane.