A Guide to Special Functions in Fractional Calculus

Abstract

Dedicated to the memory of Professor Richard Askey (1933–2019) and to pay tribute to the Bateman Project. Harry Bateman planned his “shoe-boxes” project (accomplished after his death as Higher Transcendental Functions, Vols. 1–3, 1953–1955, under the editorship by A. Erdélyi) as a “Guide to the Functions”. This inspired the author to use the modified title of the present survey. Most of the standard (classical) Special Functions are representable in terms of the Meijer G-function and, specially, of the generalized hypergeometric functions ${}_p{F}_q$. These appeared as solutions of differential equations in mathematical physics and other applied sciences that are of integer order, usually of second order. However, recently, mathematical models of fractional order are preferred because they reflect more adequately the nature and various social events, and these needs attracted attention to “new” classes of special functions as their solutions, the so-called Special Functions of Fractional Calculus (SF of FC). Generally, under this notion, we have in mind the Fox H-functions, their most widely used cases of the Wright generalized hypergeometric functions ${}_p{\Psi }_q$ and, in particular, the Mittag–Leffler type functions, among them the “Queen function of fractional calculus”, the Mittag–Leffler function. These fractional indices/parameters extensions of the classical special functions became an unavoidable tool when fractalized models of phenomena and events are treated. Here, we try to review some of the basic results on the theory of the SF of FC, obtained in the author’s works for more than 30 years, and support the wide spreading and important role of these functions by several examples.

1. Historical Introduction

Special functions are particular mathematical functions that have more or less established names and notations due to their importance in mathematical analysis, functional analysis, geometry, physics, astronomy, statistics or other applications (Wikipedia: Special Functions [1]). It might be Euler, who started to talk, since 1720, about lots of the standard special functions. He defined the Gamma-function as a continuation of the factorial, also the Bessel functions and looked after the elliptic functions. Several (theoretical and applied) scientists started to use such functions, introduced their notations and named them after famous contributors. Thus, the notions as the Bessel and cylindrical functions; the Gauss, Kummer, Tricomi, confluent and generalized hypergeometric functions; the classical orthogonal polynomials (as Laguerre, Jacobi, Gegenbauer, Legendre, Tchebisheff, Hermite, etc.); the incomplete Gamma- and Beta-functions; and the Error functions, the Airy, Whittaker, etc. functions appeared and a long list of handbooks on the so-called “Special Functions of Mathematical Physics” or “Named Functions” (we call them also “Classical Special Functions”) were published. We mention only some of them in this survey.

As Richard Askey (to whose memory we dedicate this survey) confessed in his lectures [2] on orthogonal polynomials and special functions: “Now, there are relatively large number of people who know a fair amount about this topic. Nevertheless, …most of the mathematicians are totally unaware of the power of the special functions. They react to a paper which contains a Bessel function or Legendre polynomial by turning immediately to the next paper”, and continued: “Hopefully these lectures will show … how useful hypergeometric functions can be. Very few facts about them are known, but these few facts can be very useful in many different contexts. So, my advice is to learn something about hypergeometric functions: or, if this seems too hard or dull a task, get to know someone who knows something about them. If you already know something about these functions, share your knowledge with a colleague or two, or a group of students. Every large university and research laboratory should have a person who not only finds things in the Bateman Project (i.e., [3]), but can fill in a few holes in this set of books … In any case, I hope my point has been made; special functions are useful and those who need them and those who know them should start to talk to each other … The mathematical community at large needs the education on the usefulness of special functions more than most other people who could use them …”.

As a participant in the NATO International Conference on Special Functions and Applications 2000 (Arizona State University), the author had the chance to witness the late night long discussions (mainly between Richard Askey and Oleg Marichev) for the merits and competition of the two great projects on Special Functions, on which the Computer Algebra Systems packages Mathematica and Maple are based, the Bateman Project [3] and the NIST Project [4] based on the Abramowitz–Stegun handbook [5] and on a more recent one, edited by Olver–Lozier–Boisvert–Clark [6].

The author of this survey was tempted to start paying attention to Special Functions by having the handbook [3] on her desk, while working on a M.Sc. thesis. We cite some texts from the Preface of this encyclopedia book, known as the Bateman Project: “…During his last years he (Professor Harry Bateman) had embarked upon a project whose successful completion, he believed, would prove of great value to scientists in all fields. He planned an extensive compilation of “special functions”—solutions of a wide class of mathematically and physically relevant functional equations. He intended to investigate and to tabulate properties of such functions, inter-relations between such functions, their representations in various forms, their macro- and micro-scopic behavior, and to construct tables of important definite integrals involving such functions …anyone who has been faced with the task of handling and discussing and understanding in detail the solution to an applied problem which is described by a differential equation is painfully familiar with the disproportionately large amount of scattered research on special functions one must wade through in the hope of extracting the desired information …”. In the time of Bateman’s death (1946) his notes amounted to a veritable mountain of paper. His card-catalogue alone filled several dozen cardboard boxes (the famous “shoe-boxes”). …“Bateman planned his Project as a ‘Guide to the Functions’ on a gigantic scale ….the great importance of such a work hardly needs emphasis …(this) would have made this book as a kind of ‘Greater Oxford Dictionary of Special Functions’ (from the Introduction to [3])”. This project resulted in publication of five important reference volumes ([3,7]), under the editorship of Arthur Erdélyi, in association with W. Magnus, F. Oberhettinger and F.G. Tricomi.

In 2007, the Askey–Bateman Project was announced by Mourad Ismail as a five- or six-volume encyclopedic book series on special functions, based on the works of both Harry Bateman and Richard Askey. Starting in 2020, Cambridge University Press began publishing Volumes 1 and 2 of this Encyclopedia of Special Functions with series editors Mourad Ismail and Walter Van Assche [8].

2. Preliminaries—Basic Definitions and Facts

We give here only a short background on the considered Special Functions of Fractional Calculus (SF of FC). As for the standard special functions and same for the SF of FC, most of them are entire functions of the complex variable z or analytic ones in disks in ℂ. We skip the details on defining single-valued branches of the considered functions, functional spaces and operators’ properties there (see our previous works, e.g., ref. [9] (§5.5.i)). In addition, we limit ourselves to the Fox H-functions of one complex variable, as enough general level to expose our approach and results.

Among the long list of handbooks and surveys dedicated not only to classical SF but also to the SF of FC, we mention only few of them of the few decades: Mathai–Saxena [10], 1973; Marichev [11], 1978; Srivastava–Gupta–Goyal [12], 1982; Srivastava–Kashyap [13], 1982; Prudnikov–Brychkov–Marichev [14], 1992; Kiryakova [9], 1994; Yakubovich–Luchko [15], 1994; Podlubny [16], 1999; Kilbas–Saigo [17], 2004; Kilbas–Srivastava–Trujillo [18], 2006; Mathai–Haubold [19], 2008; Mathai–Saxena–Haubold [20], 2010; Mainardi [21], 2010; Gorenflo–Kilbas–Mainardi–Rogosin [22], 2014–2020; the recent ones as Cohl–Ismail [23], 2020; Assche–Ismail [8], 2020; Mainardi [24], 2020; etc. (see more sources also the survey paper Machado-Kiryakova [25]). The basic classes of SF considered here are shortly discussed below.

Definitions of the Basic Special Functions

We refer to the survey by Mainardi–Pagnini [26] that points out the pioneering role of Salvatore Pincherle on developing the generalized hypergeometric functions (and, thus, later appearing G-functions) by means of Mellin–Barnes integrals, where a historical note from the Bateman Project [3] (Vol. 1, p. 49) is cited: “... Of all integrals which contain Gamma functions in their integrands the most important ones are the so-called Mellin-Barnes integrals. Such integrals were first introduced by S. Pincherle, in 1888; their theory has been developed in 1910 by H. Mellin … and they were used for a complete integration of the hypergeometric differential equation by E.W. Barnes, 1908.”

Definition 1.

(Ch. Fox [27], 1961, see books as [9,12,14,18], and other earlier and latest ones) The Fox H-function is a generalized hypergeometric function, defined by means of the Mellin–Barnes type contour integral

$$H_{p,q}^{m,n}\left[z\left|\begin{array}{c}{(a_i,A_i)}_1^p\\ {(b_j,B_j)}_1^q\end{array}\right.\right]=\frac{1}{2\pi i}\underset{ℐ}{\int }{\mathcal{H}}_{p,q}^{m,n}\left(s\right)z^{-s}ds,with{\displaystyle \underset{p,q}{\overset{m,n}{\mathcal{H}}}}\left(s\right)=\frac{{\displaystyle \prod _{j=1}^m}\Gamma (b_j+B_{j}s){\displaystyle \prod _{i=1}^n}\Gamma (1-a_i-A_{i}s)}{{\displaystyle \prod _{j=m+1}^q}\Gamma (1-b_j-B_{j}s){\displaystyle \prod _{i=n+1}^p}\Gamma (a_i+A_{i}s)},$$

(1)

with complex variable $z\ne 0$ and a contour $ℐ$ in the complex domain; the orders $(m,n,p,q)$ are non-negative integers so that $0\le m\le q$, $0\le n\le p$, the parameters $A_i>0,B_j>0$ are positive and $a_i,b_j$, $i=1,\dots ,p;j=1,\dots ,q$ are arbitrary complex such that $A_i(b_j+l)\ne B_j(a_i-l^{\prime }-1),l,l^{\prime }=0,1,2,\dots ;i=1,\dots ,n;j=1,\dots ,m$. Note that the integrand ${\mathcal{H}}_{p,q}^{m,n}\left(s\right)$ with $s\mapsto -s$ is the Mellin transform of the H-function (1).

The details on the properties of the Fox H-function and types of contour L can be found in many contemporary handbooks on SF as [12,14,18], where its behavior is described in terms of the following parameters:

$$\begin{array}{c}\rho ={\displaystyle \prod _{i=1}^p}{A}_i^{-A_i}{\displaystyle \prod _{j=1}^q}{B}_j^{B_j};\Delta ={\displaystyle \sum _{j=1}^q}{B}_j-{\displaystyle \sum _{i=1}^p}{A}_i;\gamma =\underset{s\to \infty ,s\in {ℐ}_{i\infty }}{lim}\text{Re}s,\\ \mu ={\displaystyle \sum _{j=1}^q}{b}_j-{\displaystyle \sum _{i=1}^p}{a}_i+\frac{p-q}{2};a^{*}={\displaystyle \sum _{i=1}^n}{A}_i-{\displaystyle \sum _{i=n+1}^p}{A}_i+{\displaystyle \sum _{j=1}^m}{B}_j-{\displaystyle \sum _{j=m+1}^q}{B}_j.\end{array}$$

(2)

Depending on the values in (2), the H-function is a function analytic of z in disks $\left|z\right|<\rho$ or outside them, in some sectors, or in the whole complex plane. In particular, the integral (1) converges (see [14] (§8.3)), if one of the following conditions is satisfied: (1) $ℐ={ℐ}_{i\infty }$: $a^{*}>0,\left|\text{arg}z\right|<a^{*}\pi /2$; (2) $ℐ={ℐ}_{i\infty }$: $a^{*}\ge 0$, $\left|\text{arg}z\right|=a^{*}\pi /2$, $\gamma \Delta <-1-\text{Re}\mu$; (3) $ℐ={ℐ}_{-i\infty }$: $\Delta >0,0<|z<\infty$, or $\Delta =0,0<\left|z\right|<\rho$, or $\Delta =0,a^{*}\ge 0$, $\left|z\right|=\rho ,\text{Re}\mu <0$; or (4) $ℐ={ℐ}_{+i\infty }$: $\Delta <0,0<|z<\infty$, or $\Delta =0,\left|z\right|>\rho$, or $\Delta =0,a^{*}\ge 0$, $\left|z\right|=\rho ,\text{Re}\mu <0$. The contour ${ℐ}_{-i\infty }$ (respectively, ${ℐ}_{+i\infty }$) is a left (respectively, right) loop in some horizontal strip that begins the point $-\infty +i{\phi }_1$ (respectively, $+\infty +i{\phi }_1$) keeping all poles of the functions $\Gamma (b_j+B_{j}s),j=1,2,\dots ,m$ on the left side, and those of $\Gamma (1-a_i-A_{i}s),i=1,2,\dots n$ on its right side, and ends at the point $-\infty +i{\phi }_2$ (respectively, $+\infty +i{\phi }_2$), where ${\phi }_1<{\phi }_2$. The contour ${ℐ}_{i\infty }$ starts at the point $\gamma -i\infty$ and ends at $\gamma +i\infty$ in a way to separate the mentioned poles, same as for ${ℐ}_{\pm i\infty }$.

For studies on the behavior of the H-function around the singular points, one can see also the work of Karp [28], commenting and revisiting the results of Braaksma [29].

If all $A_i=B_j=1,i=1,\dots ,p;j=1,\dots ,q$, the H-function $H_{p,q}^{m,n}\left[z\left|\begin{array}{c}{(a_i,1)}_1^p\\ {(b_j,1)}_1^q\end{array}\right.\right]$ reduces to the Meijer G-function (C.S. Meijer [30], 1936–1941; see details in [3] (Vol. 1) and all above-mentioned books):

$$G_{p,q}^{m,n}\left[z\left|\begin{array}{c}{\left(a_i\right)}_1^p\\ {\left(b_j\right)}_1^q\end{array}\right.\right]=\frac{1}{2\pi i}\underset{ℐ}{\int }{\mathcal{G}}_{p,q}^{m,n}\left(s\right)z^{-s}ds$$

$$=\frac{1}{2\pi i}\underset{ℐ}{\int }\frac{{\displaystyle \prod _{j=1}^m}\Gamma (b_j+s){\displaystyle \prod _{i=1}^n}\Gamma (1-a_i-s)}{{\displaystyle \prod _{j=m+1}^q}\Gamma (1-b_j-s){\displaystyle \prod _{i=n+1}^p}\Gamma (a_i+s)}{z}^{-s}ds,z\ne 0.$$

(3)

In this case, the behavior of the function (3) depends on conditions (2) with $\rho =1$, $\Delta =q-p$, $\delta =m+n-\frac{p+q}{2}$. Although simpler than (1), the G-function is yet enough general as it incorporates most of the Classical SF (known also as Named SF) and many elementary functions (see lists of examples in [3] (Vol. 1), [9] (Appendix C)).

The basic SF of FC that are Fox H-functions but do not reduce to Meijer G-functions in the general case (of irrational $A_j,B_k$) are the following generalized hypergeometric functions, extending the more popular ${}_p{F}_q$-functions.

Definition 2.

(see [9,14,22]) The Wright generalized hypergeometric function ${}_p{\Psi }_q\left(z\right)$, called also Fox–Wright function (abbrev. as F-W g.h.f. or Wright g.h.f.), is defined as:

$${}_p{\Psi }_q\left[\left.\begin{array}{c}(a_1,A_1),\dots ,(a_p,A_p)\\ (b_1,B_1),\dots ,(b_q,B_q)\end{array}\right|z\right]={\displaystyle \sum _{k=0}^{\infty }}{\displaystyle \frac{\Gamma (a_1+k{A}_1)\dots \Gamma (a_p+k{A}_p)}{\Gamma (b_1+k{B}_1)\dots \Gamma (b_q+k{B}_q)}}{\displaystyle \frac{z^k}{k!}}$$

(4)

$$=H_{p,q+1}^{1,p}\left[-z\left|\begin{array}{c}(1-a_1,A_1),\dots ,(1-a_p,A_p)\\ (0,1),(1-b_1,B_1),\dots ,(1-b_q,B_q)\end{array}\right.\right].$$

(5)

It was introduced and studied by Sir Edward Maitland (E.-M.) Wright in a series of his works (e.g., [31,32], pp. 1933–1940). In denotations for the parameters (2), the power series (4) defines an entire function of z if $\Delta >-1$; it is absolutely convergent in the disk $\left\{\right|z|<\rho \}$ for $\Delta =-1$; and it is the same for $\left|z\right|=\rho$ if $Re\left(\mu \right)>1/2$, (see details, for example, in [33]).

When all $A_1=\dots =A_p=1,B_1=\dots =B_q=1$, the Wright g.h.f. reduces to the generalized hypergeometric ${}_p{F}_q$-function which itself is a case of the G-function (3) (for early details, see [3] (Vol. 1)):

$${}_p{\Psi }_q\left[\left.\begin{array}{c}(a_1,1),\dots ,(a_p,1)\\ (b_1,1),\dots ,(b_q,1)\end{array}\right|z\right]=c{}_p{F}_q(a_1,\dots ,a_p;b_1,\dots ,b_q;z)={\displaystyle \sum _{k=0}^{\infty }}\frac{{\left(a_1\right)}_k\dots {\left(a_p\right)}_k}{{\left(b_1\right)}_k\dots {\left(b_q\right)}_k}\frac{z^k}{k!}$$

(6)

$$=G_{p,q+1}^{1,p}\left[-z\left|\begin{array}{c}1-a_1,\dots ,1-a_p\\ 0,1-b_1,\dots ,1-b_q\end{array}\right.\right];$$

where

$$c=\left[{\displaystyle \prod _{i=1}^p}\Gamma \left(a_i\right)/{\displaystyle \prod _{j=1}^q}\Gamma \left(b_j\right)\right],{\left(a\right)}_k:=\Gamma (a+k)/\Gamma \left(a\right).$$

In general (that is, except in certain integer values of parameters when the series terminates or fails to make sense), ${}_p{F}_q$ converges for all finite z if $p\le q$, converges for $\left|z\right|<1$ if $p=q+1$ and diverges for all $z\ne 0$ if $p>q+1$. The simplest particular cases are the Gauss hypergeometric function ${}_2{F}_1$, the Kummer (confluent hypergeometric) function ${}_1{F}_1$ and the Bessel type functions ${}_0{F}_1$.

A very special and important case of SF of FC, as a H-function and also as a Wright ${}_p{\Psi }_q$-function, is the “Queen”-function of FC (see [34]), namely the Mittag–Leffler (M-L) function, which has recently enjoyed many extensions (along with many basic elementary and known SF as its particular cases) and wide applications in solutions of fractional order models. This is the topic of Section 4 and Section 5.

3. On the Use of Some G- and H-Functions in Theory of Integral Transforms and Special Functions

The Meijer G-function includes most of elementary and special functions (the classical ones) as particular cases, one can find lists of these, say in [3] (Vol. 1), [9] (Appendix C), [11,14]. Naturally, the more general Fox H-function incorporates all cases of the G-functions, and much more the SF of FC. Here, we attract readers’ attention to the use of two basic classes of G- and H-functions with specific orders: $G_{0,m}^{m,0}$, respectively, $H_{0,m}^{m,0}$ with $m=q$, $n=p=0$; and $G_{m,m}^{m,0}$, respectively, $H_{m,m}^{m,0}$ with $m=p=q,n=0$.

3.1. Use of G- and H-Functions as Kernels of Laplace Type Integral Transforms

The Laplace transform

$$\mathcal{L}\{f\left(t\right);s\}={\displaystyle \underset{0}{\overset{\infty }{\int }}}exp(-st)f\left(t\right)dt,\text{Re}s>\mu ,$$

(7)

is usually considered for functions $f\left(t\right)$ of the form

$$\left\{f\left(t\right)=t^p\tilde{f}\left(t\right),p>-1,\tilde{f}\in C[0,\infty );f\left(t\right)=O(exp\mu t),t\to \infty ,\mu \in ℝ\right\}.$$

Definition 3.

The G- and H-transforms (see, for example [35], also [15,36]) of the form

$$\mathcal{G}\{f\left(t\right);s\}={\displaystyle \underset{0}{\overset{\infty }{\int }}}{G}_{p,q}^{m,n}\left[st\left|\begin{array}{c}{\left(a_j\right)}_1^p\\ {\left(b_k\right)}_1^q\end{array}\right.\right]f\left(t\right)dt,resp.\mathcal{H}\{f\left(t\right);s\}={\displaystyle \underset{0}{\overset{\infty }{\int }}}{H}_{p,q}^{m,n}\left[st\left|\begin{array}{c}{(a_j,A_j)}_1^p\\ {(b_k,B_k)}_1^q\end{array}\right.\right]f\left(t\right)dt,$$

are said to be generalized integral transforms of Laplace type when

$$\delta =m+n-\frac{p+q}{2}>0,resp.a^{*}={\displaystyle \sum _{j=1}^n}{A}_j-{\displaystyle \sum _{j=n+1}^p}{A}_j+{\displaystyle \sum _{k=1}^m}{B}_k-{\displaystyle \sum _{k=m+1}^q}{B}_k>0,$$

and are considered in suitable functional spaces of “transformable” functions.

In 1958, Obrechkoff [37] introduced a far reaching generalization of the Laplace and Meijer transforms, particular cases of which were studied by many authors years later, mainly for the purposes of operational calculi for different classes of differential operators. His aims were to extend the theorem of S. Bernstein for absolutely monotonic functions representable by means of Laplace–Stieltjes transforms, when the conditions for nth derivatives are replaced by similar ones with more general differential operators. The Obrechkoff transform was defined as

$$\mathcal{F}\left(s\right)={\displaystyle \underset{0}{\overset{\infty }{\int }}}\Phi \left(st\right)f\left(t\right)dt$$

with a kernel $\Phi \left(s\right)$ given by the integral representation

$$\Phi \left(s\right)={\displaystyle \underset{0}{\overset{\infty }{\int }}}\dots {\displaystyle \underset{0}{\overset{\infty }{\int }}}{u}_1^{{\beta }_1}\dots u^{{\beta }_p}exp\left(-u_1-\dots -u_p-\frac{s}{u_1\dots u_p}\right)d{u}_1\dots d{u}_p.$$

(8)

Later, in 1966, Dimovski [38] introduced a class of differential operators of Bessel type and of arbitrary integer order $m>1$, called by the author (as for example in [9]) as hyper-Bessel differential operators. They have the alternative representations

$$\begin{array}{ll}Bf\left(t\right)& =t^{{\alpha }_0}\frac{d}{dt}{t}^{{\alpha }_1}\frac{d}{dt}\dots t^{{\alpha }_{m-1}}\frac{d}{dt}{t}^{{\alpha }_m}f\left(t\right)\\ & =t^{-\beta }{P}_m\left(t\frac{d}{dt}\right)f\left(t\right)=t^{-\beta }{\displaystyle \prod _{k=1}^m}\left(t\frac{d}{dt}+\beta {\gamma }_k\right)f\left(t\right),t>0,\end{array}$$

(9)

with arbitrary parameters ${\alpha }_0,{\alpha }_1,\dots ,{\alpha }_m$, $\beta :=m-({\alpha }_0+{\alpha }_1+\dots +{\alpha }_m)>0$, ${\gamma }_k:=\frac{1}{\beta }({\alpha }_k+{\alpha }_{k+1}+\dots +{\alpha }_m)$, $k=1,\dots ,m$, $P_m$ a polynomial of degree m. Evidently for $m=\beta =2$, ${\gamma }_{1,2}=\pm \frac{\nu }{2}$, one has the second-order Bessel differential operator $B_{\nu }$ with the Bessel function $y\left(t\right)=J_{\nu }\left(t\right)$ satisfying $B_{\nu }y\left(t\right)=-y\left(t\right)$. For other choices of parameters, many other particular differential operators appear in equations of mathematical physics, operational calculus and applied analysis. To combine the Mikusinski type algebraic approach to operational calculus for (9) with a transform method, Dimovski used a modified Obrechkoff transform (we shortly call it also Obrechkoff transform), defined as

$$\mathcal{O}\left\{f\left(t\right);s\right\}=\beta {\displaystyle \underset{0}{\overset{\infty }{\int }}}{t}^{\beta ({\gamma }_m+1)-1}K\left[{\left(st\right)}^{\beta }\right]f\left(t\right)dt=\beta {\displaystyle \underset{0}{\overset{\infty }{\int }}}\lambda (t,s)f\left(t\right)dt,$$

with the kernel-function

$${\displaystyle K\left(s\right)={\displaystyle \underset{0}{\overset{\infty }{\int }}}\dots {\displaystyle \underset{0}{\overset{\infty }{\int }}}exp\left(-u_1-\dots -u_{m-1}-{\displaystyle \frac{s}{u_1\dots u_{m-1}}}\right){\displaystyle \prod _{k=1}^m}{u}_k^{{\gamma }_m-{\gamma }_k-1}d{u}_1\dots d{u}_{m-1}.}$$

(10)

In [9] (Ch.3), also in other works like [39], we proved that the kernel-functions (8) and (10) of the Obrechkoff transforms are representable as Meijer’s $G_{0,m}^{m,0}$-functions, namely (for a proof see, e.g., Lemma 1 of [39]):

$$\Phi \left(s\right)=G_{0,p+1}^{p+1,0}\left[s\left|\begin{array}{c}--\\ {({\beta }_k+1)}_1^p,0\end{array}\right.\right],\lambda (t,s)=s^{-\beta ({\gamma }_m+1)+1}{G}_{m,m}^{m,0}\left[{\left(st\right)}^{\beta }\left|\begin{array}{c}--\\ {({\gamma }_k-\frac{1}{\beta }+1)}_1^m\end{array}\right.\right].$$

(11)

Therefore, the Obrechkoff transform appears to be a G-transform of Laplace type (because $\delta =m/2>0$), and its theory has been further developed in whole details (convolution theorems, real and complex inversion formulas, images, examples, etc.) more easily by using the tools of the G-functions (see for example [9] (Ch.3), [39]).

Another not observed fact was that functions of the form of kernels (8) and (10) of the Obrehkoff transform were studied yet in 1937 by Erdélyi [40]. He might be the first who derived a relation between the ${}_0{F}_{m-1}$-functions (what we mention next as hyper-Bessel functions) and these kernel-functions (formula (7.4) in [40]). However, at the time of Erdélyi’s work [40], 1937, the next step in introducing the G-functions had not yet been done by Meijer [30], 1946. Obrechkoff himself made no attempts to identify the kernel-function $\Phi \left(s\right)$ with some known special functions and studied its properties “ad hoc”. Thus, the $G_{0,m}^{m,0}$-functions seemed to appear in use for the hyper-Bessel operators and related integral transforms in author’s works since 1980 (see [9] (Ch.3), [41]).

Next, the generalized Obrechkoff transform (a fractionalized analog) was introduced and studied by Kiryakova [9] (Ch.5), Al-Mussalam–Kiryakova–Tuan [42] and Yakubovich–Luchko [15], with the Fox $H_{m,m}^{m,0}$-function as kernel:

$$\mathcal{B}\left(s\right)={\mathcal{B}}_{\left({\rho }_i\right),\left({\mu }_i\right)}\left\{f\left(t\right);s\right\}={\displaystyle \underset{0}{\overset{\infty }{\int }}}{H}_{0,m}^{m,0}\left[st\left|\begin{array}{c}--\\ {({\mu }_i-\frac{1}{{\rho }_i},\frac{1}{{\rho }_i})}_1^m\end{array}\right.\right]f\left(t\right)dt.$$

(12)

We call it as multi-index Borel–Dzrbashjan transform, because for $m=1$ it is reduced to the Borel transform

$${\mathcal{B}}_{\left(\rho \right),\left(\mu \right)}\left\{f\left(t\right);s\right\}={\displaystyle \underset{0}{\overset{\infty }{\int }}}exp(-s^{\rho }{t}^{\rho })t^{\mu \rho -1}f\left(t\right)dt$$

(13)

whose kernel appears to be a $H_{0,1}^{1,0}$-function. This integral transform was shown by Dzrbashjan [43,44] to have inversion formula involving the Mittag–Leffler function $E_{1/\rho ,\mu }$. The generalized Obrechkoff transform (12) is a tool in operational calculus for fractional multi-order analogs of hyper-Bessel differential operators (9), formally of the kind

$$D_{\left({\rho }_i\right),\left({\mu }_i\right)}f\left(t\right)=t^{-1}{\displaystyle \prod _{i=1}^m}\left(t^{1+(1-{\mu }_i){\rho }_i}{D}_{t^{{\rho }_i}}^{1/{\rho }_i}{t}^{({\mu }_i-1){\rho }_i}\right)f\left(t\right),$$

(14)

in the same way as the Laplace transform, the Obrechkoff transform and its particular cases serve for the classical differentiation, respectively for the hyper-Bessel operators (9).

In the studies on these Laplace type G- and H-integral transforms, we used essentially the theory of the G- and H-functions, mainly of the cases of orders $(m,0;0,m)$. Note that, for example, $G_{0,m}^{m,0}\left(s\right)$ is an analytic function in the sector $|args|<(m/2)\pi$ (where in this case $\delta =m/2>0$). Some additional necessary results on these G- and H- kernel functions were derived by Kiryakova [9] (Appendix), as Lemmas B.1–B.4, Corollaries B.5–B.7, Formula (E.21), etc.

From the known representations of some elementary and special functions in G- and H-terms, one observes many particular cases of simpler Laplace type integral transforms. Namely, the Laplace and Borel–Dzrbashjan transforms (7) and (13) are Obrechkoff transform (10) and multi-index Borel–Dzrbashjan transform (12), respectively, for $m=1$, since

$$exp(-s)=G_{0,1}^{1,0}\left[s\left|\begin{array}{c}--\\ 0\end{array}\right.\right],exp(-s^{\rho }{t}^{\rho })=H_{0,1}^{1,0}\left[st\left|\begin{array}{c}--\\ (\mu -\frac{1}{\rho },\frac{1}{\rho })\end{array}\right.\right].$$

For $m=2$, we have the classical Meijer transform as a case of the Obrechkoff transform, related to the Bessel differential operator $B_{\nu }={\displaystyle \frac{d}{dt}}{t}^{1-\nu }{\displaystyle \frac{d}{dt}}{t}^{\nu }$:

$${\mathcal{K}}_{\nu }\left\{f\left(t\right);s\right\}={\displaystyle \underset{0}{\overset{\infty }{\int }}}\sqrt{st}{K}_{\nu }\left(t\right)f\left(t\right)dt,since{K}_{\nu }\left(s\right)=\frac{1}{2}{G}_{0,2}^{2,0}\left[\frac{s^2}{4}\left|\begin{array}{c}--\\ \frac{\nu }{2},\frac{-\nu }{2}\end{array}\right.\right],$$

(15)

the kernel Macdonald function $K_{\nu }\left(s\right)$ has such a G-function representation.

In a series of papers [45,46], Krätzel introduced a generalization of the Meijer transform (again with $m=2$), and further a more general one of the type of Obrechkoff transform for arbitrary integer $m>1$,

$${\mathcal{L}}_{\nu }^{\left(m\right)}\left\{f\left(t\right);s\right\}={\displaystyle \underset{0}{\overset{\infty }{\int }}}{\lambda }_{\nu }^{\left(m\right)}\left[m{\left(st\right)}^{1/m}\right]f\left(t\right)dt:={\displaystyle \underset{0}{\overset{\infty }{\int }}}\Lambda (s,t)f\left(t\right)dt.$$

(16)

He used the transformation (16) for operational calculus for the following (hyper-Bessel type) operator or order $m>1$:

$$B_{\nu }^{\left(m\right)}=\frac{d}{dt}{t}^{\frac{1}{m}-\nu }{\left(t^{1-\frac{1}{m}}\frac{d}{dt}\right)}^{m-1}{t}^{\nu +1-\frac{2}{m}}\mathrm{with}\beta =1,{\gamma }_1=0,{\gamma }_k=\nu +\frac{k-2}{m},k=2,\dots ,m.$$

(17)

As expected, we can represent the Krätzel kernel in terms of the G-function corresponding to (11):

$$\Lambda (s,t)={\displaystyle \underset{0}{\overset{\infty }{\int }}}\dots {\displaystyle \underset{0}{\overset{\infty }{\int }}}\left[{\displaystyle \prod _{k=1}^{m-1}}{u}_k^{\nu -1+\frac{k-1}{m}}\right]exp\left(-u_1-\dots -u_{m-1}-\frac{st}{u_1...u_{m-1}}\right)d{u}_1\dots d{u}_{m-1}$$

$$=s^{-\nu -1+\frac{2}{m}}{G}_{0,m}^{m,0}\left[st\left|\begin{array}{c}--\\ 0,(\nu +\frac{k-2}{m}{)}_2^m\end{array}\right.\right].$$

(18)

Krätzel started from the simple case $m=2$ with a kernel of the form ${\displaystyle \underset{0}{\overset{\infty }{\int }}}{u}^{\gamma -1}exp(-u-st/u)du$ (with some variations as $t\mapsto t^{\rho },\rho >0$), close to the Macdonald function (15), which is often called as the Krätzel function. Then, many other authors continued to study it and established its relations to hypergeometric functions. We can refer to such works by Kilbas–Saxena–Trujillo [47], Mathai–Haubold [48], etc. In a paper by Glaeske–Kilbas–Saigo [49], a fractionalized analog of the Krätzel transform (7) was introduced, where instead of integer $m>1$ in the transformation (16), they took a (fractional) parameter $\rho >0$. Then, naturally, its kernel is represented as a H-function (due to some variations in the definition, it appears as $H_{1,2}^{2,0}$ instead of $H_{0,2}^{2,0}$). Relations with operators of fractional calculus are considered, but one should mention that such an integral transform is analog of the generalized Obrechkoff transform (12) for a fractional order differential operator of the form (14). In all these mentioned cases, the kernel functions have the form of (8) as also studied earlier by Erdélyi [40]. We conclude here the list of cases of the Obrechkoff transform with emphasize on the works by Ditkin–Prudnikov (as [50]) on operational calculi for (hyper-Bessel) operators of the form

$$B_1=\frac{d}{dt}t\frac{d}{dt},\text{and more generally},B_m=\frac{d}{dt}t\frac{d}{dt}t\frac{d}{dt}\dots \frac{d}{dt}=t^{-1}{\left(t\frac{d}{dt}\right)}^m,m\ge 2.$$

(19)

For $m=2$, the corresponding integral transform is a variant of the Meijer transform (with $\nu =0$), and in the general case $m>1$, Ditkin and Prudnikov [50] made use of an integral transform of the form

$\mathcal{B}\left\{f\left(t\right);s\right\}=2{\int }_0^{\infty }{E}_{0m}\left(st\right)f\left(t\right)dt,$ where we can represent the kernel $E_{0m}$ as a $G_{0,m}^{m,0}$ − function.

For more details on the Obrechkoff type transforms with kernels $G_{0,m}^{m,0}$ and $H_{0,m}^{m,0}$, their properties, images and special cases, see Kiryakova [9] (Ch.3, Ch.5), [39].

3.2. Use of G- and H-Functions as Kernels in Generalized Fractional Calculus

For basic background on Fractional Calculus (FC) as theory of operators of integration and differentiation of arbitrary (fractional) order, and its closely related topics as special functions (SF) and integral transforms, we refer to the books by Samko–Kilbas–Marichev [51], Podlubny [16], Kilbas–Srivastava–Trujillo [18], and Yakubovich–Luchko [15], as well as one by the author [9], among many others. For wider list, see, for example, Machado–Kiryakova [25]. In our works, and mainly for the needs of the SF theory, we consider the Riemann–Liouville (R-L) type integrals and their corresponding derivatives of R-L and Caputo type, respectively, their generalizations involving G- and H-functions in the kernels. Note that we concentrate on the left-hand side variants and skip details (in most cases being similar) for the Weyl-type, right-hand sided operators.

The basic tools in our studies are the fractional integration operators of the form $\tilde{I}f\left(z\right)=z^{{\delta }_0}{I}_{\beta }^{\gamma ,\delta }f\left(z\right),{\delta }_0\ge 0$, to which we refer as “classical fractional integrals”, where

$$I_{\beta }^{\gamma ,\delta }f\left(z\right)=\frac{1}{\Gamma \left(\delta \right)}{\displaystyle \underset{0}{\overset{1}{\int }}}{(1-\sigma )}^{\delta -1}{\sigma }^{\gamma }f\left(z{\sigma }^{\frac{1}{\beta }}\right)d\sigma =\frac{z^{-\beta (\gamma +\delta )}}{\Gamma \left(\delta \right)}{\displaystyle \underset{0}{\overset{z}{\int }}}{(z^{\beta }-{\xi }^{\beta })}^{\delta -1}{\xi }^{\beta \gamma }f\left(\xi \right)d\left({\xi }^{\beta }\right),$$

(20)

is the Erdélyi-Kober operator (E-K) of integration of order $\delta \ge 0$, depending on two additional parameters $\gamma \in ℝ,\beta >0$. In this general form, it is introduced in Sneddon’s works [52] and considered in some books (for example, [9] (Ch.2), [15,18,51]). The earlier versions with $\beta =1$, $\beta =2$ are due to Kober and Erdélyi. The R-L operator of integration $R^{\delta }$ appears as a case with one parameter only, for $\gamma =0$, $\beta =1$, ${\delta }_0=\delta \ge 0$,

$$R_{0+,z}^{\delta }f\left(z\right):=R^{\delta }f\left(z\right)=z^{\delta }{I}_1^{0,\delta }f\left(z\right);\mathrm{conversely},I_1^{\gamma ,\delta }f\left(z\right)=z^{-\gamma -\delta }{R}^{\delta }{z}^{\gamma }f\left(z\right).$$

(21)

The E-K fractional derivative $D_{\beta }^{\gamma ,\delta }$, corresponding to (20), is defined explicitly almost simultaneously in the works of Kiryakova [9] (Ch.2) and Yakubovich–Luchko [15] (Ch.3). It serves as an interpretation of the formal inversion formula ${\left\{I_{\beta }^{\gamma ,\delta }\right\}}^{-1}=I_{\beta }^{\gamma +\delta ,-\delta }$, namely:

$$D_{\beta }^{\gamma ,\delta }f\left(z\right)=D_n{I}_{\beta }^{\gamma +\delta ,n-\delta }f\left(z\right)={\displaystyle \prod _{j=1}^n}\left(\frac{1}{\beta }z\frac{d}{dz}+\gamma +j\right)I_{\beta }^{\gamma +\delta ,n-\delta }f\left(z\right),n-1<\delta \le n,n\in \mathbb{N}.$$

(22)

Here, the simplest integer order derivative ${(d/dz)}^n$ in the definition of the R-L fractional derivative $D^{\delta }f\left(z\right):={\left({\displaystyle \frac{d}{dz}}\right)}^n{R}^{n-\delta }f\left(z\right)$, is replaced by an auxiliary differential operator $D_n$ of integer order, a polynomial of $(zd/dz)$. The Caputo-type R-L and E-K fractional derivatives are defined in the same way but with exchanged order of the nonnegative order integration and the integer order differentiation (see, e.g., [53]).

The notion for generalized operators of fractional integration was introduced by Kalla in his 1969–1979 works (see the survey [54]), who suggested their common form

$$If\left(z\right)={\displaystyle \underset{0}{\overset{1}{\int }}}\Phi \left(\sigma \right){\sigma }^{\gamma }f\left(z\sigma \right)d\sigma =z^{-\gamma -1}{\displaystyle \underset{0}{\overset{z}{\int }}}\Phi \left(\frac{\xi }{z}\right){\xi }^{\gamma }f\left(\xi \right)d\xi ,$$

where $\Phi \left(\sigma \right)$ can be an arbitrary continuous (analytical) function for which the integral makes sense. The idea of such generalized fractional calculus is to replace the elementary function in the kernel of R-L and E-K operators (20) (and, say, the logarithmic kernel in the Hadamard integral) by some special function. Variants with the Gauss-, Bessel-, Whittaker-, arbitrary G- and H-functions appeared in papers of several authors (see historical details and references in [54,55]). If such a special function $\Phi$ is taken to be too general or too specific, only some formal operational rules for the corresponding fractional calculus can be derived. The lucky hint in our studies was to choose suitably the kernel-functions Φ to be of the form of $G_{m,m}^{m,0}$- and $H_{m,m}^{m,0}$-functions. Then, the operators of the generalized fractional calculus happen to be also commutative products of classical operators of FC, namely of finite number of Erdélyi-Kober operators. Thus, the tools of the special functions and the wide use of the classical FC are combined into a Generalized Fractional Calculus (GFC) in Kiryakova [9], with developed full theory and many illustrated applications in different areas of analysis, differential equations, special functions and integral transforms. Below, we briefly review the basic definitions and few results on this GFC.

Definition 4.

(Kiryakova, [9] (Ch.5)) We define the multiple E-K integral (of multiplicity $m>1$), by means of the real parameters’ sets $({\delta }_1\ge 0,\dots ,{\delta }_m\ge 0)$ (multi-order of integration) $({\gamma }_1,\dots ,{\gamma }_m)$ (multi-weight) and $({\beta }_1>0,\dots ,{\beta }_m>0)$ (additional multi-parameter), as:

$$I_{\left({\beta }_k\right),m}^{\left({\gamma }_k\right),\left({\delta }_k\right)}f\left(z\right):={\displaystyle \underset{0}{\overset{1}{\int }}}{H}_{m,m}^{m,0}\left[\sigma \left|\begin{array}{c}{({\gamma }_k+{\delta }_k+1-\frac{1}{{\beta }_k},\frac{1}{{\beta }_k})}_1^m\\ {({\gamma }_k+1-\frac{1}{{\beta }_k},\frac{1}{{\beta }_k})}_1^m\end{array}\right.\right]f\left(z\sigma \right)d\sigma ,$$

(23)

if ${\displaystyle \sum _{k=1}^m}{\delta }_k>0$; and as the identity operator: $I_{\left({\beta }_k\right),m}^{\left({\gamma }_k\right),(0,...,0)}f\left(z\right)=f\left(z\right)$, if ${\delta }_1={\delta }_2=\dots ={\delta }_m=0$.

It is important to mention that, for the particular conditions (2), the above kernel $H_{m,m}^{m,0}$-function is analytic function in the unit disk and $H_{m,m}^{m,0}\left(\sigma \right)\equiv 0$ for $\left|\sigma \right|>1$ (Kiryakova, ref. [9]).

In the case of all equal βs: ${\beta }_1={\beta }_2=\dots ={\beta }_m=:\beta >0$, integral (23) has a simpler form with a Meijer $G_{m,m}^{m,0}$-function ([9] (Ch.1)), which is also analytic in unit disk and $G_{m,m}^{m,0}\left(\sigma \right)\equiv 0$ for $\left|\sigma \right|>1$,

$$I_{(\beta ,\dots ,\beta ),m}^{\left({\gamma }_k\right),\left({\delta }_k\right)}f\left(z\right):=I_{\beta ,m}^{\left({\gamma }_k\right),\left({\delta }_k\right)}f\left(z\right)={\displaystyle \underset{0}{\overset{1}{\int }}}{G}_{m,m}^{m,0}\left[\sigma \left|\begin{array}{c}{({\gamma }_k+{\delta }_k)}_1^m\\ {\left({\gamma }_k\right)}_1^m\end{array}\right.\right]f\left(z{\sigma }^{1/\beta }\right)d\sigma =\left[{\displaystyle \prod _{k=1}^m}{I}_{\beta }^{{\gamma }_k,{\delta }_k}\right]f\left(z\right).$$

(24)

In both cases of (23) and (24), the operators of the form

$$\tilde{I}f\left(z\right)=z^{{\delta }_0}{I}_{\left({\beta }_k\right),m}^{\left({\gamma }_k\right),\left({\delta }_k\right)}f\left(z\right),\tilde{I}f\left(z\right)=z^{{\delta }_0}{I}_{\beta ,m}^{\left({\gamma }_k\right),\left({\beta }_k\right)}f\left(z\right),with{\delta }_0\ge 0,$$

(25)

are called generalized fractional integrals of multi-order $({\delta }_1,\dots ,{\delta }_m)$.

The important decomposition property (for proof, see for example, [9] (Th.1.2.10, Th.5.2.1), says that the same GFC integrals (23) and (24) can be represented, instead of using the kernel H- and G-functions, by repeated integral representations for the commutative product of classical E-K operators (20):

$$I_{\left({\beta }_k\right),m}^{\left({\gamma }_k\right),\left({\delta }_k\right)}f\left(z\right):=\left[{\displaystyle \prod _{k=1}^m}{I}_{{\beta }_k}^{{\gamma }_k,{\delta }_k}\right]f\left(z\right)$$

$$={\displaystyle \underset{0}{\overset{1}{\int }}}\dots {\displaystyle \underset{0}{\overset{1}{\int }}}\left[{\displaystyle \prod _{k=1}^m}\frac{{(1-{\sigma }_k)}^{{\delta }_k-1}{\sigma }_k^{{\gamma }_k}}{\Gamma \left({\delta }_k\right)}\right]f\left(z{\sigma }_1^{1/{\beta }_1}\dots {\sigma }_m^{1/{\beta }_m}\right)d{\sigma }_1\dots d{\sigma }_m.$$

(26)

In the book [9] and subsequent papers, we provided a full set of operational properties of the operators (23) and (24) that justify their names as operators of GFC, as semigroup property, formal inversion formula, reduction to identity or to the conventional integration operators for special parameters’ choice.

Analogously to the R-L and E-K fractional derivatives, we define the corresponding generalized fractional derivatives. The auxiliary differential operator $D_{\eta }$ is chosen on the base of the specific differential relations for the kernel function, derived for the G-functions, and especially for $G_{m,m}^{m,0}$ by Kiryakova [9] (App., Lemmas B.3, B.4, Cor. B.6) and for $H_{m,m}^{m,0}$ by Kiryakova [9] (Ch.5, Lemma 5.1.7)and Kiryakova–Luchko [53] (Lemma 18).

Definition 5.

(Kiryakova [9]) Let $D_{\eta }$ be the following polynomial of $z\left({\displaystyle \frac{d}{dz}}\right)$ of degree ${\eta }_1+\dots +{\eta }_m$:

$$D_{\eta }=\left[{\displaystyle \prod _{r=1}^m}{\displaystyle \prod _{j=1}^{{\eta }_r}}\left(\frac{1}{{\beta }_r}z\frac{d}{dz}+{\gamma }_r+j\right)\right],with{\eta }_k:=\left\{\begin{array}{cc}\left[{\delta }_k\right]+1,\hfill & for noninteger{\delta }_k,\hfill \\ {\delta }_k,\hfill & for integer{\delta }_k,\hfill \end{array}\right.k=1,\dots ,m.$$

(27)

The multiple (m-tuple) Erdélyi–Kober fractional derivative of R-L type of multi-order $({\delta }_1\ge 0,\dots ,{\delta }_m\ge 0)$ is defined by means of the differ-integral operator:

$$D_{\left({\beta }_k\right),m}^{\left({\gamma }_k\right),\left({\delta }_k\right)}f\left(z\right):=D_{\eta }{I}_{\left({\beta }_k\right),m}^{({\gamma }_k+{\delta }_k),({\eta }_k-{\delta }_k)}f\left(z\right)=D_{\eta }{\displaystyle \underset{0}{\overset{1}{\int }}}{H}_{m,m}^{m,0}\left[\sigma \left|\begin{array}{c}{({\gamma }_k+{\eta }_k+1-\frac{1}{{\beta }_k},\frac{1}{{\beta }_k})}_1^m\\ {({\gamma }_k+1-\frac{1}{{\beta }_k},\frac{1}{{\beta }_k})}_1^m\end{array}\right.\right]f\left(z\sigma \right)d\sigma .$$

(28)

Similarly, the Caputo-type generalized fractional derivative was introduced by Kiryakova and Luchko [53], as

$${}^{*}{D}_{\left({\beta }_k\right),m}^{\left({\gamma }_k\right),\left({\delta }_k\right)}f\left(z\right)=I_{\left({\beta }_k\right),m}^{({\gamma }_k+{\delta }_k),({\eta }_k-{\delta }_k)}{D}_{\eta }f\left(z\right).$$

(29)

In the case ${\beta }_1=\dots ={\beta }_m:=\beta >0$, simpler representations involving the Meijer G-function hold for the R-L and Caputo-type “derivatives” which correspond to the generalized fractional integral (24):

$$\begin{array}{cc}\hfill & D_{\beta ,m}^{\left({\gamma }_k\right),\left({\delta }_k\right)}f\left(z\right)=D_{\eta }{I}_{\beta ,m}^{({\gamma }_k+{\delta }_k),({\eta }_k-{\delta }_k)}f\left(z\right)={\displaystyle \left[{\displaystyle \prod _{r=1}^m}{\displaystyle \prod _{j=1}^{{\eta }_r}}\left(\frac{1}{\beta }z\frac{d}{dz}+{\gamma }_r+j\right)\right]I_{\beta ,m}^{({\gamma }_k+{\delta }_k),({\eta }_k-{\delta }_k)}f\left(z\right)},\\ & {}^{*}{D}_{\beta ,m}^{\left({\gamma }_k\right),\left({\delta }_k\right)}f\left(z\right)=I_{\beta ,m}^{({\gamma }_k+{\delta }_k),({\eta }_k-{\delta }_k)}{D}_{\eta }f\left(z\right).\end{array}$$

(30)

More generally, the differ-integral/integro-differential operators of the form

$$\tilde{D}f\left(z\right)=D_{\left({\beta }_k\right),m}^{\left({\gamma }_k\right),\left({\delta }_k\right)}{z}^{-{\delta }_0}f\left(z\right)=z^{-{\delta }_0}{D}_{\left({\beta }_k\right),m}^{({\gamma }_k-\frac{{\delta }_0}{\beta }),\left({\delta }_k\right)}f\left(z\right),\text{and}$$

$$\widetilde{{}^{*}D}f\left(z\right)={}^{*}{D}_{\left({\beta }_k\right),m}^{\left({\gamma }_k\right),\left({\delta }_k\right)}{z}^{-{\delta }_0}f\left(z\right)\text{with}{\delta }_0\ge 0,$$

(31)

are all called generalized (multiple, multi-order) fractional derivatives (of R-L or Caputo type).

Next, in Section 8, we often use also the notion of (generalized) fractional differintegrals. We have in mind either (generalized) fractional integrals or derivatives or compositions of some E-K fractional integrals and some E-K fractional derivatives. These appear as meanings of operators (26) when part of the order’s δs are non-negative and the other parts are negative.

For the functional spaces (here, we mainly limit to weighted analytical functions of complex z), mapping properties, long list of operational properties, images, etc., we refer, for example, to the work of Kiryakova [9,53,56].

We use also a further extension of the generalized fractional integrals (23), based on the so-called Wright–Erdélyi–Kober (W-E-K) operator of fractional integration (see [57]), with parameters as in E-K integral: $\delta \ge 0$, $\gamma$ real, $\beta >0$ and additional parameter $\lambda >0$, where the Wright–Bessel (Bessel–Maitland) function of the form $J_{\nu }^{\mu }$ (see (57)) in Section 5) is used in the kernel:

$$W_{\beta ,\lambda }^{\gamma ,\delta }f(z):=I_{\beta ,\lambda ,1}^{\gamma ,\delta }f\left(z\right)=\lambda {\displaystyle \underset{0}{\overset{1}{\int }}}{\sigma }^{\lambda (\gamma +1)-1}{J}_{\gamma +\delta -\lambda (\gamma +1)/\beta }^{-\lambda /\beta }({\sigma }^{\lambda }f\left(z\sigma \right)d\sigma .$$

(32)

One can show that, for $\lambda =\beta$, the above kernel-function reduces to the kernel of the E-K operator, therefore the W-E-K integration becomes the E-K one. Using compositions of W-E-K operators (32), Kalla and Galue [57] tried to develop a next step in the generalized fractional calculus with $H_{m,m}^{m,0}$ kernel-functions that have the same structure but different parameters λ_ks and λ_ks in upper and low rows. Some revisions and properties of these operators were further provided by Kiryakova [58,59,60].

Definition 6.

For integer $m\ge 1$ and real parameters ${\delta }_k\ge 0$, ${\gamma }_k$, ${\beta }_k>0,{\lambda }_k>0,{\beta }_k\ge {\lambda }_k$, $k=1,\dots ,m$, we define the multiple Wright–Erdélyi–Kober (W-E-K) fractional integrals, as follows:

$$\tilde{I}f\left(z\right)=I_{\left({\beta }_k\right),\left({\lambda }_k\right),m}^{\left({\gamma }_k\right),\left({\delta }_k\right)}f\left(z\right):={\displaystyle \underset{0}{\overset{1}{\int }}}{H}_{m,m}^{m,0}\left[\sigma \left|\begin{array}{c}{({\gamma }_i+{\delta }_i+1-\frac{1}{{\beta }_i},\frac{1}{{\beta }_i})}_1^m\\ {({\gamma }_i+1-\frac{1}{{\lambda }_i},\frac{1}{{\lambda }_i})}_1^m\end{array}\right.\right]f\left(z\sigma \right)d\sigma =\left[{\displaystyle \prod _{k=0}^m}{W}_{{\beta }_k,{\lambda }_k}^{{\gamma }_k,{\delta }_k}\right]f\left(z\right),$$

(33)

if ${\displaystyle \sum _{i=1}^m}{\delta }_i>0$; and as the identity operator: $\tilde{I}f\left(z\right)=f\left(z\right)$, when ${\delta }_1={\delta }_2=\dots ={\delta }_m=0$ and ${\lambda }_k={\beta }_k,k=1,\dots ,m$. For ${\gamma }_k>-1,k=1,\dots ,m$ and the above-mentioned conditions on the other parameters, the operators (33) are shown to preserve the space of analytic functions in disks or in starlike complex domains.

If ${\beta }_k={\lambda }_k,k=1,\dots ,m$, the “new” operators of GFC (33) coincide with operators (23). The corresponding generalized fractional derivatives $D_{\left({\beta }_k\right),\left({\lambda }_k\right),m}^{\left({\gamma }_k\right),\left({\delta }_k\right)}$ are defined by means of differential-integral operators similar to those for (28).

Here, we mention some few of the numerous special cases of the above defined GFC operators, to emphasize the particular elementary and special functions appearing in their kernels, and thus as cases of the kernel $H_{m,m}^{m,0}$- and $G_{m,m}^{m,0}$-functions.

For $m=1$, we have the kernel-functions:

$$H_{1,1}^{1,0}\left[\sigma \left|\begin{array}{c}(\gamma +\delta ,1/\beta )\\ (\gamma ,1/\beta )\end{array}\right.\right]=\beta {\sigma }^{\beta -1}{G}_{1,1}^{1,0}\left[{\sigma }^{\beta }\left|\begin{array}{c}\gamma +\delta \\ \gamma \end{array}\right.\right]=\beta \frac{{\sigma }^{\beta \gamma +\beta -1}{(1-{\sigma }^{\beta })}^{\delta -1}}{\Gamma \left(\delta \right)},$$

(34)

thus the generalized fractional integrals and derivatives (23) and (28) reduce to the corresponding E-K (20) and (22) and R-L operators (21): $I_{\beta ,1}^{\gamma ,\delta }=I_{\beta }^{\gamma ,\delta }$, $D_{\beta ,1}^{\gamma ,\delta }=D_{\beta }^{\gamma ,\delta }$, $R^{\delta }$ and $D^{\delta }$. Many other integration and differentiation operators introduced and used by different authors appear as their special cases.

In the case $m=2$, the kernel functions $H_{2,2}^{2,0}$ and $G_{2,2}^{2,0}$ reduce to a Gauss hypergeometric function or its variations, for example:

$$H_{2,2}^{2,0}\left[\sigma \left|\begin{array}{c}({\gamma }_1+{\delta }_1+1-\frac{1}{\beta },\frac{1}{\beta }),({\gamma }_2+{\delta }_2+1-\frac{1}{\beta },\frac{1}{\beta })\\ ({\gamma }_1+1-\frac{1}{\beta },\frac{1}{\beta }),({\gamma }_2+1-\frac{1}{\beta },\frac{1}{\beta })\end{array}\right.\right]=G_{2,2}^{2,0}\left[{\sigma }^{\beta }\left|\begin{array}{c}{\gamma }_1+{\delta }_1,{\gamma }_2+{\delta }_2\\ {\gamma }_1,{\gamma }_2\end{array}\right.\right]$$

$$=\frac{{\sigma }^{\beta {\gamma }_2}{(1-{\sigma }^{\beta })}^{{\delta }_1+{\delta }_2-1}}{\Gamma ({\delta }_1+{\delta }_2)}{}_2{F}_1({\gamma }_2+{\delta }_2-{\gamma }_1,{\delta }_1;{\delta }_1+{\delta }_2;1-{\sigma }^{\beta }).$$

(35)

Therefore, the generalized fractional integrals in this case are known as hypergeometric fractional integrals; some of them were introduced and studied by, e.g. Love, Saxena, Saigo and Hohlov (see [54]).

For $m=3$, we have as special case the Marichev–Saigo–Maeda (M-S-M) operators of FC, the integration operators introduced and studied by Marichev (1974) and Saigo et al. (1996, 1998) (see [55]). This is because their kernel-function, the Appel $F_3$ function (Horn function)

$F_3\left(a,a^{\prime },b,b^{\prime },c,z,\xi \right)={\displaystyle \sum _{m,n=0}^{\infty }}\frac{{\left(a\right)}_m{\left(a^{\prime }\right)}_n{\left(b\right)}_m{\left(b^{\prime }\right)}_n}{{\left(c\right)}_{m+n}}\frac{z^m{\xi }^n}{m!n!},\left|z\right|<1,\left|\xi \right|<1$ (see, e.g., [3,14]) is a case of the GFC kernel-functions $H_{3,3}^{3,0}$ and $G_{3,3}^{3,0}$ (see, for example, [14], §8.4.51, (2)):

$$\frac{{(1-\sigma )}^{c-1}}{\Gamma \left(c\right)}{F}_3\left(a,a^{\prime },b,b^{\prime },c,1-\frac{1}{\sigma },1-\sigma \right)$$

$$=G_{3,3}^{3,0}\left[\sigma \left|\begin{array}{c}a+b,c-a^{\prime },c-b^{\prime }\\ a,b,c-a^{\prime }-b^{\prime }\end{array}\right.\right]=H_{3,3}^{3,0}\left[\sigma \left|\begin{array}{c}(a+b,1),(c-a^{\prime },1),(c-b^{\prime },1)\\ (a,1),(b,1),(c-a^{\prime }-b^{\prime },1)\end{array}\right.\right],\text{Re}c>0.$$

(36)

Let $m\ge 1$ be an arbitrary integer, but all δs be equal integers, say ${\delta }_1=\dots ={\delta }_m=1$. Then, from (24), we obtain the hyper-Bessel integral operators L (we denote below their kernel by $G_1$) that correspond to the hyper-Bessel differential operators (9) of arbitrary (higher) integer order $m>1$. In practice, these are operators of integer multi-orders $(1,1,\dots ,1)$, but their fractional powers ${ℐ}^{\lambda }$, $\lambda >0$ have been represented (Kiryakova [9,41]) as GFC integrals of multi-order $(\lambda ,\lambda ,\dots ,\lambda )$ with kernels $G_{\lambda }$, where

$$G_1\left(\sigma \right)=G_{m,m}^{m,0}\left[\sigma \left|\begin{array}{c}{({\gamma }_k+1)}_1^m\\ {\left({\gamma }_k\right)}_1^m\end{array}\right.\right],G_{\lambda }\left(\sigma \right)=G_{m,m}^{m,0}\left[\sigma \left|\begin{array}{c}{({\gamma }_k+\lambda )}_1^m\\ {\left({\gamma }_k\right)}_1^m\end{array}\right.\right].$$

The kernel of ${ℐ}^{\lambda }$ in the form $G_{\lambda }$ appeared also in the work of McBride [61]. These expressions gave us the hint how to introduce our GFC, replacing $(\lambda ,\lambda ,\dots ,\lambda )$ by arbitrary fractional multi-order $({\delta }_1,{\delta }_2,\dots ,{\delta }_m)$, explanations are in [41]. We can mention also the Gelfond–Leontiev [62] operator (47) generated by the multi-index M-L functions (see next section and the works by Kiryakova [63,64]), as a more general example of operators of fractional multi-order where the Fox $H_{m,m}^{m,0}$-functions serve as kernels.

The H-functions of the form $H_{p,q}^{q,0}$, of which the kernel functions of (23) are cases with $p=q=m$, were studied in series of papers by Karp. In [28], he revisited the Braaksma results [29] for the H-function’s behavior in the neighborhood of the singular points and its analytical continuation. There he commented also works on applications of H-functions not only in fractional calculus, but also widely in statistics, including the book by Mathai–Saxena–Haubold [20].

In relation to the use of the $G_{m,m}^{m,0}$-functions (the kernel-functions of GFC integrals (24)) in applications to statistics, it is interesting to note that, in 1958, Kabe [65] explored them in statistics, as density functions of a random variable. He also distinguished the cases $m=1$ and $m=2$ (mentioned above) related to the kernel-functions of the E-K and of the hypergeometric fractional integrals, respectively, (34) and (35). Studies on the closely related $G_{m+1,m+1}^{m,1}$-functions as R-L integrals of $G_{m,m}^{m,0}$ can be found in the work by Karp [66].

4. Mittag–Leffler Functions and Their Extensions

The Mittag–Leffler (M-L) function $E_{\alpha }\left(z\right)$ was introduced by G. Mittag–Leffler ([67], 1902–1905), extended to 2-parameters as $E_{\alpha ,\beta }\left(z\right)$ by A. Wiman [68] and studied later by P. Humbert and R.P. Agarwal [69]. It was presented in the Bateman Project [3], Vol. 3 (1954), in a chapter for “Miscellaneous Functions”. However, for long time, it was ignored in the other handbooks on special functions because the applied scientists suffered from the lack of tables for its Laplace transforms. Although arising from the studies of Mittag–Leffler on a problem not related to fractional calculus, but on analytical continuation of series to maximal starlike domain (Mittag–Leffler star), nowadays, the M-L function is the most popular and most exploited SF of FC. It was titled as the “Queen”-function of FC by Gorenflo and Mainardi in 1997 (see also the very recent survey by Mainardi [34]). The basic theory and more details, can be found, for example, in [22,43] (see also, e.g., [9,16,70,71]).

Definition 7.

The Mittag–Leffler (M-L) functions $E_{\alpha }$ and $E_{\alpha ,\beta }$, are entire functions of order $\rho =1/\alpha$ and type 1, defined by the power series

$$E_{\alpha }\left(z\right)={\displaystyle \sum _{k=0}^{\infty }}{\displaystyle \frac{z^k}{\Gamma (\alpha k+1)}},E_{\alpha ,\beta }\left(z\right)={\displaystyle \sum _{k=0}^{\infty }}{\displaystyle \frac{z^k}{\Gamma (\alpha k+\beta )}},\alpha >0,\beta >0.$$

(37)

As “fractional index” $(\alpha >0)$ analogs of the exponential and trigonometric functions that satisfy ODEs of first and second order ($\alpha =1,2$), the M-L functions serve as solutions of fractional order differential and integral equations. An example is the Rabotnov function, called also “fractional exponent”, $y\left(z\right)=z^{\alpha -1}{E}_{\alpha ,\alpha }\left(z^{\alpha }\right)$ that solves the simplest fractional order differential equation $D^{\alpha }y\left(z\right)=y\left(z\right)$. Let us refer also to the pioneering work by Hille–Tamarkin [72], where the solution of the Abel integral equation of the second kind was provided in terms of a M-L function. As far as the Laplace transform images are mentioned, one can find these for the M-L type functions and their kth derivatives in the work of Podlubny ([16] (S.1.2.2)):

$$\mathcal{L}\left\{z^{\alpha k+\beta -1}{E}_{\alpha ,\beta }^{\left(k\right)}(\pm \lambda z^{\alpha }));s\right\}=(\frac{k!s^{\alpha -\beta }}{s^{\alpha }\mp \lambda }{)}^{k+1},\text{Re}s>{\left|\lambda \right|}^{1/\alpha }.$$

A Mittag–Leffler type function with three indices, known as the Prabhakar function [73], is also often studied and used (for details, see [22,70,71,74,75] and other contemporary books and surveys on M-L type functions):

$$E_{\alpha ,\beta }^{\gamma }\left(z\right)={\displaystyle \sum _{k=0}^{\infty }}\frac{{\left(\gamma \right)}_k}{\Gamma (\alpha k+\beta )}\frac{z^k}{k!},\alpha ,\beta ,\gamma \in C,\text{Re}\alpha >0;$$

(38)

where ${\left(\gamma \right)}_0=1,{\left(\gamma \right)}_k=\Gamma (\gamma +k)/\Gamma \left(\gamma \right)$ denotes the Pochhammer symbol. Its Laplace transform has the form

$$\mathcal{L}\left\{E_{\alpha ,\beta }^{\gamma }\left(\lambda z^{\alpha }\right);s\right\}=\frac{s^{-\beta }}{{(1-\lambda s^{-\alpha })}^{\gamma }}.$$

For $\gamma =1$, we get the M-L function $E_{\alpha ,\beta }$, and, if additionally $\beta =1$, then it is $E_{\alpha }$.

These M-L type functions are simple cases of the Wright g.h.f. and of the H-function, namely:

$$E_{\alpha ,\beta }\left(z\right)={}_1{\Psi }_1\left[\left.\begin{array}{c}(1,1)\\ (\beta ,\alpha )\end{array}\right|z\right]=H_{1,2}^{1,1}\left[-z\left|\begin{array}{c}(0,1)\\ (0,1),(1-\beta ,\alpha )\end{array}\right.\right],$$

$$E_{\alpha ,\beta }^{\gamma }\left(z\right)=\frac{1}{\Gamma \left(\gamma \right)}{}_1{\Psi }_1\left[\left.\begin{array}{c}(\gamma ,1)\\ (\beta ,\alpha )\end{array}\right|z\right]=H_{1,2}^{1,1}\left[-z\left|\begin{array}{c}(1-\gamma ,1)\\ (0,1),(1-\beta ,\alpha )\end{array}\right.\right].$$

Another generalization of M-L function (37) with additional parameters, for example $l\in ℂ,\mu \in ℝ$, was considered by Gorenflo–Kilbas–Rogosin [76], and its relations to FC operators:

$$E_{\alpha ,\mu ,l}\left(z\right)={\displaystyle \sum _{k=0}^{\infty }}{c}_k{z}^k,\mathrm{with}{c}_k={\displaystyle \prod _{j=0}^{k-1}}\frac{\Gamma \left[\alpha \right(j\mu +l)+1]}{\Gamma \left(\alpha \right(j\mu +l+1)+1]}.$$

A vector index extension of (37) appeared in the works by Luchko et al. (e.g., [15,77,78]) on operational calculus’ methods for some fractional order PDE and multi-term FO differential equations. Under the name multi-index (multiple) M-L function, it was introduced by Kiryakova [63,79] using a different approach, as to be the generating function of Gelfond–Leontiev generalized integration and differentiation operators (47) (see Definition 9) and inspired from the paper by Dzrbashjan [44] on M-L type function with $2\times 2$ indices. Further, this class of functions were studied in details by Kiryakova [59,80], Kilbas–Koroleva–Rogosin [81], Paneva–Konovska [74] and many other followers. Luchko et al. also considered multivariate analogs of the so-called vector index M-L functions [78].

Definition 8. (Kiryakova [59,80]) Let $m>1$ be an integer, $({\alpha }_1>0,{\alpha }_2>0,\dots ,{\alpha }_m>0)$ and $({\beta }_1,{\beta }_2,\dots ,{\beta }_m)$ be arbitrary real parameters. By means of these sets of “multi-indices”, the multi-index Mittag–Leffler function (abbrev. as multi-M-L f.) is defined as:

$$E_{\left({\alpha }_i\right),\left({\beta }_i\right)}\left(z\right):=E_{\left({\alpha }_i\right),\left({\beta }_i\right)}^{\left(m\right)}\left(z\right)={\displaystyle \sum _{k=0}^{\infty }}\frac{z^k}{\Gamma ({\alpha }_{1}k+{\beta }_1)\dots \Gamma ({\alpha }_{m}k+{\beta }_m)}.$$

(39)

Under weakened restrictions on αs (or their real parts) not to be obligatory all non-negative, the study was extended by Kilbas et al. [81].

As a further extension of both Prabhakar function (38) and of the ($2m$) multi-index M-L functions (39), Paneva–Konovska [74,82] introduced and studied the so-called $\left(3m\right)$-parametric (multi-index) Mittag–Leffler functions, similar to (39) but with additional set of parameters $({\gamma }_1,\dots ,{\gamma }_m)$:

$$E_{\left({\alpha }_i\right),\left({\beta }_i\right)}^{\left({\gamma }_i\right),m}\left(z\right)={\displaystyle \sum _{k=0}^{\infty }}\frac{{\left({\gamma }_1\right)}_k\dots {\left({\gamma }_m\right)}_k}{\Gamma ({\alpha }_{1}k+{\beta }_1)\dots \Gamma ({\alpha }_{m}k+{\beta }_m)}\frac{z^k}{{(k!)}^m}.$$

(40)

For $m=1$, one has the Prabhakar function, and, for ${\gamma }_1=\dots ={\gamma }_m=1$, these are (39).

The Mellin transforms of (39), (40) and their particular cases can be found in [83].

The so-called Le Roy type function has been an object of several recent studies, e.g., by Gerhold [84], Garra–Polito [85], Garrappa–Rogosin–Mainardi [86], Garrappa–Orsingher–Polito [87], as a new special function

$$F_{\alpha ,\beta }^{\left(\gamma \right)}\left(z\right)={\displaystyle \sum _{k=0}^{\infty }}\frac{z^k}{{[\Gamma (\alpha k+\beta )]}^{\gamma }},$$

(41)

which is an entire function of $z\in ℂ$ for parameters $\text{Re}\left(\alpha \right)>0$, $\beta \in ℝ$ and $\gamma >0$. This resembles to the M-L function (for $\gamma =1$) and to the multi-index M-L function (39) (for integer $\gamma =m$, all ${\alpha }_i=\alpha ,{\beta }_i=\beta ,i=1,\dots ,m$). The function (41) appeared as extension of the function $R_{\gamma }\left(z\right)={\displaystyle \sum _{k=0}^{\infty }}{z}^k/{[(k+1)!]}^{\gamma }$, introduced by E. Le Roy [88] (1899), similarly to the purposes of G. Mittag–Leffler [67] (1903) to study analytical continuations of the sums of power series, and it seems they were working in competition on such ideas. Similar to the M-L type functions, (41) is involved in solutions of various problems, including a Convey–Maxwell–Poison distribution for different degrees of over- and under-dispersion.

Some Basic Properties of the Multi-Index Mittag-Leffler Functions

The basic properties and results for the functions (39) and long lists of their examples, all of them having wide applications in solutions of integer- and fractional-order models, are provided in our previous papers (e.g., [59,60,79,80]). Some of them are reminded here.

Theorem 1.

The multi-index M-L functions (39) are entire functions with the following order ρ and type σ:

$$\frac{1}{\rho }={\alpha }_1+\dots +{\alpha }_m,\frac{1}{\sigma }={\left(\rho {\alpha }_1\right)}^{\rho {\alpha }_1}\dots {\left(\rho {\alpha }_m\right)}^{\rho {\alpha }_m},$$

(42)

respectively with α_is replaced by $Re\left({\alpha }_i\right)$s. Note that the type $\sigma >1$ for $m>1$ and only for $m=1$ (classical case (37)): $\sigma =1$. The following asymptotic estimate holds:

$$|E_{\left({\alpha }_i\right),\left({\beta }_i\right)}{\left(z\right)|\le C|z|}^{\rho \left(\right(1/2)+\mu -(m/2\left)\right)}\text{exp}\left({\sigma \left|z\right|}^{\rho }\right),\mu :={\beta }_1+\dots +{\beta }_m,for\left|z\right|\to \infty .$$

The $\left(3m\right)$-parameters M-L type functions (40) are also entire functions with the same order and type as in (42), see [74,82].

Lemma 1.

The multi-index M-L functions (39) are important examples of the Wright generalized hypergeometric functions ${}_p{\Psi }_q$ and of the Fox H-functions:

$$E_{\left({\alpha }_i\right),\left({\beta }_i\right)}\left(z\right)=E_{\left({\alpha }_i\right),\left({\beta }_i\right)}^{\left(m\right)}\left(z\right)={}_1{\Psi }_m\left[\left.\begin{array}{c}(1,1)\\ {({\beta }_i,{\alpha }_i)}_1^m\end{array}\right|z\right]=H_{1,m+1}^{1,1}\left[-z\left|\begin{array}{c}(0,1)\\ (0,1),{(1-{\beta }_i,{\alpha }_i)}_1^m\end{array}\right.\right].$$

(43)

Thus, the following Mellin–Barnes type integral representation holds (cf. with (1)):

$$E_{\left({\alpha }_i\right),\left({\beta }_i\right)}\left(z\right)=\frac{1}{2\pi i}\underset{ℐ}{\int }\frac{\Gamma \left(s\right)\Gamma (1-s)}{{\displaystyle \prod _{i=1}^m}\Gamma ({\beta }_i-s{\alpha }_i)}{(-z)}^{-s}ds,z\ne 0,$$

based on the Mellin transform (see [59,83]; also [18] (p. 48)):

$$\mathcal{M}\left\{E_{\left({\alpha }_i\right),\left({\beta }_i\right)}(-z);s\right\}=\frac{\Gamma \left(s\right)\Gamma (1-s)}{{\displaystyle \prod _{i=1}^m}\Gamma ({\beta }_i-s{\alpha }_i)},0<Re\left(s\right)<1.$$

(44)

Additionally, as shown by Paneva–Konovska [74,82], the $\left(3m\right)$-parametric functions (40) can be represented as

$$E_{\left({\alpha }_i\right),\left({\beta }_i\right)}^{\left({\gamma }_i\right),m}\left(z\right)=A{}_m{\Psi }_{2m-1}\left[\left.\genfrac{}{}{0pt}{}{({\gamma }_1,1),...,({\gamma }_m,1)}{({\beta }_1,{\alpha }_1),...,({\beta }_m,{\alpha }_m),(1,1),...,(1,1)}\right|z\right]$$

$$=A{H}_{m,2m}^{1,m}\left[-z\left|\genfrac{}{}{0pt}{}{(1-{\gamma }_1,1),...,(1-{\gamma }_m,1)}{{[(0,1),(1-{\beta }_i,{\alpha }_i)]}_1^m}\right.\right],\mathrm{with}A={\left[{\displaystyle \prod _{i=1}^m}\Gamma \left({\gamma }_i\right)\right]}^{-1},$$

(45)

which is in agreement with (43) for ${\gamma }_1=\dots ={\gamma }_m=1$.

As an analog of the Laplace transform ($\mathcal{L}$), relationship between the classical M-L function (37) and the classical Wright function: $\mathcal{L}\left\{\varphi (\alpha ,\beta ;z);s\right\}={\displaystyle \frac{1}{s}}{E}_{\alpha ,\beta }({\displaystyle \frac{1}{s}})$ (see in the books [16,18]), we derive the following new relation.

Lemma 2.

$$\mathcal{L}\left\{{}_0{\Psi }_m\left[\begin{array}{c}-\\ ({\beta }_1,{\alpha }_1),...,({\beta }_m,{\alpha }_m)\end{array}|z\right];s\right\}=\frac{1}{s}{E}_{\left({\alpha }_i\right),\left({\beta }_i\right)}\left(\frac{1}{s}\right),Re\left(s\right)>0.$$

(46)

Note that we can consider the ${}_0{\Psi }_m$-functions on the left-hand side as “fractional indices” analogs of the ${}_0{F}_m$-functions, that is of the hyper-Bessel functions $J_{{\nu }_1,\dots ,{\nu }_m}^{\left(m\right)}$ of Delerue [89], related to the hyper-Bessel operators (9) as their eigenfunctions, and discussed further as special cases of (39). For details on these special functions, see Kiryakova [9] (Ch.3).

Various relations for the multi-M-L functions in terms of the operators of classical FC and GFC have been derived in our previous works (e.g., [59,80]). First, let us consider the so-called Gelfond–Leontiev (G-L) operators of generalized integration and differentiation, generated by the coefficients of an entire function $\phi \left(\sigma \right)$. For the theory of the G-L operators in general, see Gelfond and Leontiev’s paper [62]) of 1951, and for details in the case when the mentioned entire function is taken to be the M-L function or multi-index Mittag–Leffler function, we refer to Kiryakova [9] (Ch.1), [59,63,79]. Here, we only remind the definition of the G-L operators related to $\phi \left(\sigma \right)=E_{\left({\alpha }_i\right),\left({\beta }_i\right)}\left(\sigma \right):={\displaystyle \sum _{k=0}^{\infty }}{b}_k{z}^k$ whose coefficients $b_k=1/\left(\Gamma ({\alpha }_{1}k+{\beta }_1)\dots \Gamma ({\alpha }_{m}k+{\beta }_k)\right)$ are taken as multipliers’ sequences below.

Definition 9.

(Kiryakova [63,64]) For functions $f\left(z\right)={\displaystyle \sum _{k=0}^{\infty }}{a}_k{z}^k$ analytic in a disk $\left\{\right|z|<R\}$, we consider the operators

$$\tilde{D}f\left(z\right):=D_{\left({\alpha }_i\right),\left({\beta }_i\right)}f\left(z\right)={\displaystyle \sum _{k=1}^{\infty }}{a}_k\frac{b_{k-1}}{b_k}{z}^{k-1},\widetilde{ℐ}f\left(z\right):={ℐ}_{\left({\alpha }_i\right),\left({\beta }_i\right)}f\left(z\right)={\displaystyle \sum _{k=0}^{\infty }}{a}_k\frac{b_{k+1}}{b_k}{z}^{k+1},$$

(47)

and call them multiple (multi-index) Dzrbashjan–Gelfond–Leointiev (D-G-L) differentiations and integrations, respectively. These are generated by the multi-index M-L functions and the name of Dzrbashjan is used in addition to Gelfond–Leontiev to honor his contribution to one of the first deep studies on M-L type functions, the book [43].

Evidently, $D_{\left({\alpha }_i\right),\left({\beta }_i\right)}{ℐ}_{\left({\alpha }_i\right),\left({\beta }_i\right)}f\left(z\right)=f\left(z\right)$, and it is proven that the radii of convergence (and analyticity) of resulting analytical functions in (47) are the same R as for $f\left(z\right)$. According to Theorem 3 in [79], operators (47) can be analytically extended outside the disks to starlike domains and represented as operators of GFC, as follows:

$$\tilde{D}f\left(z\right)=z^{-1}{D}_{(1/{\alpha }_i),m}^{({\gamma }_i-1-{\alpha }_i),\left({\alpha }_i\right)}f\left(z\right)-\left[{\displaystyle \prod _{i=1}^m}\frac{\Gamma \left({\gamma }_i\right)}{\Gamma ({\gamma }_i-{\alpha }_i)}\right]\frac{f\left(0\right)}{z},\widetilde{ℐ}f\left(z\right)=z{I}_{(1/{\alpha }_i),m}^{({\gamma }_i-1),\left({\alpha }_i\right)}f\left(z\right).$$

(48)

To start with the classical FC operators for the multi-index M-L functions, we state the following

Lemma 3.

(Kiryakova [80] (Lemma 3.2)) For any fixed $l,1\le l\le m$ and integration order ${\delta }_l>0$, we have for the E-K fractional integral the relation

$$I_{1/{\alpha }_l}^{{\gamma }_l-1,{\delta }_l}{E}_{\left({\alpha }_i\right),({\gamma }_1,...,{\gamma }_l,...,{\gamma }_m)}\left(\lambda z\right)=E_{\left({\alpha }_i\right),({\gamma }_1,...,{\gamma }_l+{\delta }_l,...,{\gamma }_m)}\left(\lambda z\right),\lambda \ne 0,$$

that is, a fractional integration can transform a multi-M-L function into another one with same α_is and corresponding parameter γ_l increased by the order of integration to ${\gamma }_l+{\delta }_l$.

Applying E-K fractional integrals of the form $I_{1/{\alpha }_i}^{{\gamma }_i-1,{\delta }_i}$ successively m-times ($i=1,\dots ,m$) to (39) and using the composition (decomposition) property (26), we obtain for the generalized fractional integrals (23) the image:

$$I_{(1/{\alpha }_i),m}^{({\gamma }_i-1),\left({\delta }_i\right)}{E}_{\left({\alpha }_i\right),\left({\gamma }_i\right)}\left(\lambda z\right)=E_{\left({\alpha }_i\right),({\gamma }_i+{\delta }_i)}\left(\lambda z\right).$$

(49)

Then, for ${\delta }_i:={\alpha }_i,i=1,\dots ,m$, and applying the operational rules for the operators $I_{\left({\beta }_i\right),m}^{\left({\gamma }_i\right),\left({\delta }_i\right)}$ and $D_{\left({\beta }_i\right),m}^{\left({\gamma }_i\right),\left({\delta }_i\right)}$ of GFC, the following generalized fractional integration and differentiation relations follow:

$$\left(\lambda z\right)I_{(1/{\alpha }_i),m}^{({\gamma }_i-1),\left({\alpha }_i\right)}{E}_{\left({\alpha }_i\right),\left({\gamma }_i\right)}\left(\lambda z\right)=E_{\left({\alpha }_i\right),\left({\gamma }_i\right)}\left(\lambda z\right)-\frac{1}{\Gamma \left({\gamma }_1\right)\dots \Gamma \left({\gamma }_m\right)},$$

$$D_{(1/{\alpha }_i),m}^{({\gamma }_i-1-{\alpha }_i),\left({\alpha }_i\right)}{E}_{\left({\alpha }_i\right),\left({\gamma }_i\right)}\left(\lambda z\right)=\left(\lambda z\right)E_{\left({\alpha }_i\right),\left({\gamma }_i\right)}\left(\lambda z\right)+\frac{1}{\Gamma ({\gamma }_1-{\alpha }_1)\dots \Gamma ({\gamma }_m-{\alpha }_m)},$$

(50)

as analogs of the classical relation $z^{\alpha }{D}^{\alpha }{E}_{\alpha }\left(\lambda z\right)=\lambda z^{\alpha }{E}_{\alpha }\left(\lambda z\right)+{\displaystyle \frac{1}{\Gamma (1-\alpha )}}$ for the R-L derivative $D^{\alpha }=z^{-\alpha }{D}_1^{-\alpha ,\alpha }$.

It remains to combine the results (48) and (50) to verify the fact that the multi-index M-L functions that generate the G-L operators (47) appear as their eigenfunctions:

Theorem 2.

The multi-index Mittag–Leffer function (39) satisfies the differential equation of fractional multi-order $({\alpha }_1,\dots ,{\alpha }_m)$:

$$\tilde{D}{E}_{\left({\alpha }_i\right),\left({\beta }_i\right)}\left(\lambda z\right)=D_{{\alpha }_i,{\beta }_i}{E}_{\left({\alpha }_i\right),\left({\beta }_i\right)}\left(\lambda z\right)=\lambda E_{\left({\alpha }_i\right),\left({\beta }_i\right)}\left(\lambda z\right),\lambda \ne 0.$$

(51)

The classical Poisson integral formula, representing the Bessel function via the cosine-function ([3] (Vol. 2)), can be written in terms of an E-K fractional integral, as

$$J_{\nu }\left(z\right)=\frac{2}{\sqrt{\pi }\Gamma (\nu +1/2)}{\left(\frac{z}{2}\right)}^{\nu }{\displaystyle \underset{0}{\overset{1}{\int }}}{(1-t^2)}^{\nu -1/2}cos\left(zt\right)dt=\frac{1}{\sqrt{\pi }}{\left(\frac{z}{2}\right)}^{\nu }{I}_{1/2}^{-1/2,\nu +1/2}\{cosz\}.$$

(52)

This representation has been extended in our works [9] (Ch.4), [90] for the hyper-Bessel functions (58), $m\ge 2$, that is for the ${}_0{F}_{m-1}$-functions, via generalized fractional integrals (24) of the function ${cos}_m$. The details follow in Section 8. For the multi-index M-L functions, a Poisson type integral representation of the kind of (52) has to explore the more general fractional calculus operators from Definition 6. This is a part of the general results discussed in Section 8, but we expose it here as to close (at least partly) the topic with some properties of the multi-index Mittag–Leffler functions.

Theorem 3.

(Kiryakova [59]) Let ${\alpha }_k>1$, ${\beta }_k\ge \frac{k}{m}$, $k=1,\dots ,m$. Then, we have the following Poisson-type integral representation of the multi-index M-L functions my means of multiple W-E-K fractional integrals (33) of the cosine function (54) of order m (from the next section):

$$E_{\left({\alpha }_k\right),\left({\beta }_k\right)}(-z)=c^{*}{I}_{{(1/{\alpha }_k)}_1^m,{\left(1\right)}_1^m,m}^{{(\frac{k}{m}-1)}_1^m,{({\beta }_k-\frac{k}{m})}_1^m}\left\{{cos}_m\left(m{z}^{1/m}\right)\right\}$$

$$=c^{*}{\displaystyle \underset{0}{\overset{1}{\int }}}{H}_{m,m}^{m,0}\left[\sigma \left|\begin{array}{c}{({\beta }_k-{\alpha }_k,{\alpha }_k)}_1^m\\ {(k/m-1,1)}_1^m\end{array}\right.\right]{cos}_m\left(m{\left(z\sigma \right)}^{1/m}\right)d\sigma ,with{c}^{*}:=\sqrt{m/{\left(2\pi \right)}^{m-1}}.$$

(53)

Remark 1.

The above result is parallel with (52) for the Bessel functions. If we take ${\alpha }_k=1$, ${\beta }_k=\frac{k}{m}$, the above GFC operator, the multiple W-E-K fractional integral, has a multi-order $(0,\dots ,0)$ and since also ${\lambda }_k={\beta }_k$, it turns into identity. Then, the $E_{\left({\alpha }_k\right),\left({\beta }_k\right)}$-function reduces to the ${cos}_m\left(z\right)$-function. It is similar in the simplest case to the Bessel function with index $\nu =-1/2$: $J_{-1/2}\left(z\right)=\sqrt{\frac{2}{\pi z}}cosz$. More generally, it is also known that the Bessel functions of “semi-integer” indices (called also“spherical functions”for their use in theory of spherical waves) are reducible to trigonometric functions or to integer order operators of them: $J_{n-1/2}\left(z\right)={\displaystyle \frac{{\left(2z\right)}^{n+1/2}}{\sqrt{\pi }}}{\displaystyle \frac{d^n}{{\left(d{z}^2\right)}^n}}\left\{{\displaystyle \frac{cosz}{z}}\right\},n=0,1,2,\dots$. In the case of multi-index M-L functions (39), we can call multi-indices of the form ${\alpha }_k=1$, ${\beta }_k:={\nu }_k-\frac{k}{m}=0,1,2,\dots$; for $k=1,\dots ,m$, as“semi-integer multi-indices”. A corollary of Theorem 3 tells that for such multi-indices the functions $E_{\left({\alpha }_k\right),\left({\beta }_k\right)}$ reduce either directly to generalized trigonometric functions, or to integer order integral or differential operators of them.

The results for the images of the multi-index Mittag–Leffler functions (39) and (40) under GFC integrals and derivatives, or under their particular cases a R-L, E-K, Saigo, Marichev–Saigo–Maeda operators, etc. can be written from the general results in Section 7 according to definition via the Wright g.h.f. ${}_1{\Psi }_m$.

Series in systems of special functions, in the general cases of $2m$- and $3m$-parameters M-L functions and their particular case (mentioned in next section) as the M-L function, Parbhakar function, multi-index and fractional analogs of the Bessel- and hyper-Bessel functions, were studied recently in details by Paneva–Konovska in a series of papers and in the book [74], especially with respect to their convergence in complex domain, including Cauchy–Hadamard, Abel, Tauber type, Hardy–Littlewood and Ostrovski type theorems.

5. Examples of M-L Type and Multi-Index M-L Functions

5.1. For $m=1$, this is the classical M-L function $E_{\alpha ,\beta }\left(z\right)$ with all its special cases:

•$\alpha >0,\beta =1$: $E_{0,1}\left(z\right)={\displaystyle \frac{1}{1-z}}$; $E_{1,1}\left(z\right)=exp\left(z\right)$; $E_{2,1}\left(z^2\right)=\cosh \left(z\right)$, $E_{2,1}(-z^2)=\cos \left(z\right)$; $E_{1/2,1}\left(z^{1/2}\right)$ = $exp\left(z\right)\left[1+\mathrm{erf}\left(z^{1/2}\right)\right]$ = $exp\left(z\right)\mathrm{erfc}(-z^{1/2})$ = $exp\left(z\right)\left[1+{\displaystyle \frac{1}{\sqrt{\pi }}}\gamma ({\displaystyle \frac{1}{2}},z)\right]$ (the error functions, or incomplete gamma functions);

•$\beta \ne 1$: $E_{1,2}\left(z\right)={\displaystyle \frac{e^z-1}{z}}$; $E_{1/2,2}\left(z\right)={\displaystyle \frac{\mathrm{sh}\sqrt{z}}{z}}$; $E_{2,2}\left(z\right)={\displaystyle \frac{\mathrm{sh}\sqrt{z}}{\sqrt{z}}}$; the Miller-Ross function $z^{\nu }{E}_{1,\nu +1}\left(az\right)$; etc.;

•$\beta =\alpha$: the α-exponential (Rabotnov) function $y_{\alpha }\left(z\right)=z^{\alpha -1}{E}_{\alpha ,\alpha }\left(z^{\alpha }\right)$.

• The trigonometric functions of order m, and, respectively the hyperbolic functions of order m:

$${cos}_m\left(z\right)={\displaystyle \sum _{j=0}^{\infty }}\frac{{(-1)}^j{z}^{mj}}{\left(mj\right)!}=E_{m,1}(-z^m),$$

(54)

$$y\left(z\right)={cos}_m\left(z\right)\mathrm{is} \mathrm{the} \mathrm{solution} \mathrm{of} \mathrm{IVP}{y}^{\left(m\right)}\left(z\right)=-y\left(z\right),y\left(0\right)=1,y^{\left(j\right)}\left(0\right)=0,j=1,\dots ,m-1;$$

$$k_r(z,m)={\displaystyle \sum _{j=0}^{\infty }}\frac{{(-1)}^j{z}^{mj+r-1}}{(mj+r-1)!}=z^{r-1}{E}_{m,r}(-z^m),r=1,2,\dots ;k_1(z,m):={\cos }_m\left(z\right)=E_{m,1}(-z^m),$$

$$h_r(z,m)={\displaystyle \sum _{j=0}^{\infty }}\frac{z^{mj+r-1}}{(mj+r-1)!}=z^{r-1}{E}_{m,r}\left(z^m\right),r=1,2,\dots ;h_1(z,m):={\cosh }_m\left(z\right)=E_{m,1}\left(z^m\right),$$

can also be expressed in terms of the M-L function (see in [3] (Vol. 3) and [16] (Ch.1)); and the same for their fractionalized versions, as by Plotnikov [91] and Tseytlin [92]:

$$S{c}_{\alpha }\left(z\right)={\displaystyle \sum _{k=0}^{\infty }}\frac{{(-1)}^k{z}^{(2-\alpha )m+1}}{\Gamma \left(\right(2-\alpha )m+2)}=z{E}_{2-\alpha ,2}(-z^{2-\alpha }),$$

$$C{s}_{\alpha }\left(z\right)={\displaystyle \sum _{k=0}^{\infty }}\frac{{(-1)}^k{z}^{(2-\alpha )m}}{\Gamma \left(\right(2-\alpha )m+1)}=E_{2-\alpha ,1}(-z^{2-\alpha }),$$

and by Luchko–Srivastava [77]:

$${sin}_{\lambda ,\mu }\left(z\right)={\displaystyle \sum _{k=0}^{\infty }}\frac{{(-1)}^k{z}^{2k+1}}{\Gamma (2\mu k+2\mu -\lambda +1)}=z{E}_{2\mu ,2\mu -\lambda +1}(-z^2),$$

$${cos}_{\lambda ,\mu }\left(z\right)={\displaystyle \sum _{k=0}^{\infty }}\frac{{(-1)}^k{z}^{2k}}{\Gamma (2\mu k+\mu -\lambda +1)}=E_{2\mu ,\mu -\lambda +1}(-z^2),$$

(see details again in Podlubny [16] (Ch.1)).

• Here, we mention also the so-called Lorenzo–Hartley functions [93], the F-function and its generalization the R-function, shown to be solutions of some linear fractional differential equations. We can represent them in terms of M-L function, namely, for $z>0$, $c=0$, $q\ge 0$, $\nu \le q$:

$$F_q(a,z)={\displaystyle \sum _{k=0}^{\infty }}\frac{a^k{z}^{(k+1)q-1}}{\Gamma \left(\right(k+1\left)q\right)}=z^{q-1}{E}_{q,q}\left(az\right),$$

$$R_{q,\nu }(a,0,z)={\displaystyle \sum _{k=0}^{\infty }}\frac{a^k{z}^{(k+1)q-1-\nu }}{\Gamma \left(\right(k+1)q-\nu )}=z^{q-1}{E}_{q,q-\nu }\left(az\right).$$

5.2. For $m=2$: We start with the not enough popular M-L type function of Dzrbashjan [44], with $2\times 2$ indices, which he denoted alternatively by (we need to set $1/{\rho }_i:={\alpha }_i,{\mu }_i:={\beta }_i,i=1,2$):

$${\Phi }_{{\rho }_1,{\rho }_2}(z;{\mu }_1,{\mu }_2)={\displaystyle \sum _{k=0}^{\infty }}\frac{z^k}{\Gamma ({\mu }_1+\frac{k}{{\rho }_1})\Gamma ({\mu }_2+\frac{k}{{\rho }_2})}:=E_{(\frac{1}{{\rho }_1},\frac{1}{{\rho }_2}),({\mu }_1,{\mu }_2)}\left(z\right)=E_{({\alpha }_1,{\alpha }_2),({\beta }_1,{\beta }_2)}\left(z\right).$$

(55)

Dzrbashjan found the order and type of this entire function, claimed on few simple particular cases, and considered some integral relations between (55) and Mellin transforms on a set of axes. Then, he developed a theory of integral transforms in the class L₂, involving kernel close to functions (55) and, further, proposed approximations of entire functions in L₂ for an arbitrary finite system of axes in complex plane starting from the origin.

The $2\times 2$-indices M-L type functions (55) were also studied in detail by Luchko in the recent paper [94]. He allowed the parameters ${\rho }_1,{\rho }_2$ to be also negative or zero and called them “4-parameters Wright functions of second kind”, separating the cases when ${\rho }_1+{\rho }_2>0$, ${\rho }_1+{\rho }_2=0$ or ${\rho }_1+{\rho }_2<0$.

Some of the simple cases of (55), as mentioned and denoted in Dzrbashjan [44], are:

• the M-L function itself: $E_{\frac{1}{\rho },\mu }\left(z\right)=E_{(\frac{1}{\rho },0),(\mu ,1)}\left(z\right)={\Phi }_{\rho ,\infty }(z;\mu ,1)$; ${\displaystyle \frac{1}{1-z}}=E_{(0,0),(1,1)}\left(z\right)={\Phi }_{\infty ,\infty }(z;1,1)$; the Bessel function: $J_{\nu }\left(z\right)={\left(\frac{z}{2}\right)}^{\nu }{E}_{(1,1),(\nu +1,1}\left(-\frac{z^2}{4}\right)={\left(\frac{z}{2}\right)}^{\nu }{\Phi }_{1,1}\left(-\frac{z^2}{4};1,\nu +1\right)$; etc.

To these examples, we added (see, e.g., Kiraykova [59]) the following cases:

• The Struve and Lommel functions (see [3] (Vol. 2); and details in [9] (App.,(C.8)), [79,80]):

$$s_{\mu ,\nu }\left(z\right)=\frac{1}{4}{z}^{\mu +1}{E}_{(1,1),\left(\right(3-\nu +\mu )/2,(3+\nu +\mu )/2)}(-\frac{z^2}{4}),H_{\nu }\left(z\right)=\frac{1}{\pi 2^{\nu -1}{(1/2)}_{\nu }}{s}_{\nu ,\nu }\left(z\right).$$

• The “classical” Wright function that arose in the studies of Fox ([95], 1928), Wright ([31], 1933) and Humbert and Agarwal ([69], 1953) and was also referred to in Erdélyi et al. [3] (Vol. 3). Initially, Wright [31] defined this function only for $\alpha >0$, then extended its definition for $\alpha >-1$ [32]. Now, we see this is a case of multi-index M-L function with $m=2$:

$$\varphi (\alpha ,\beta ;z):=W_{\alpha ,\beta }\left(z\right)={\displaystyle \sum _{k=0}^{\infty }}\frac{1}{\Gamma (\alpha k+\beta )}\frac{z^k}{k!}={}_0{\Psi }_1\left[\left.\begin{array}{c}-\\ (\beta ,\alpha )\end{array}\right|z\right]=E_{(\alpha ,1),(\beta ,1)}^{\left(2\right)}\left(z\right),$$

(56)

which is entire function of order $1/(1+\alpha )$. The survey papers by Gorenflo–Luchko–Mainardi [96] and Mainardi–Consiglio [97] reflect in detail its analytical properties and applications, see also the book [22] as well as the related literature. In the case $\alpha \ge 0$, the Wright function is said to be of first kind, and for $-1<\alpha <0$ of second kind. The latter survey [97] concentrates on the Wright function of second kind. It is noted that the first kind Wright function is of exponential order, while the second kind is not of exponential order, and naturally they have different asymptotic behaviors, Laplace transforms, etc. (see also Luchko [94]). The function (56) plays an important role in the solutions of linear partial fractional differential equations as the fractional diffusion-wave equation studied by Nigmatullin (1984–1986, to describe the diffusion process in media with fractal geometry, $0<\alpha <1$) and by Mainardi et al. (since 1994, for propagation of mechanical diffusive waves in viscoelastic media, $1<\alpha <2$). In the form $M(z;\beta )=\varphi (-\beta ,1-\beta ;-z),\beta :=\alpha /2$, this function is recently called as the Mainardi function (see [16] (Ch.1)). In our denotations, it is a multi-index M-L function with $m=2$ and a Dzrbashjan function (55): $M(z;\beta )=E_{(-\beta ,1),(1-\beta )}^{\left(2\right)}(-z)$ and has its own particular cases, such as $M(z;1/2)=1/\sqrt{\pi }exp(-z^2/4)$ and the Airy function, $M(z;1/3)=3^{2/3}Ai(z/3^{1/3})$. Note also that, for $\alpha =0$, the Wright function (56) reduces to the exponent, since $\varphi (0,\beta ;z)={\displaystyle \sum _{k=0}^{\infty }}{z}^k/(k!\Gamma \left(\beta \right))=(1/\Gamma \left(\beta \right))exp\left(z\right)$.

In alternative form and denotation, the Wright function (56) is known as the Wright–Bessel function or is misnamed as the Bessel–Maitland function:

$$J_{\nu }^{\mu }\left(z\right)=\varphi (\mu ,\nu +1;-z)={}_0{\Psi }_1\left[\left.\begin{array}{c}-\\ (\nu +1,\mu )\end{array}\right|-z\right]={\displaystyle \sum _{k=0}^{\infty }}\frac{{(-z)}^k}{\Gamma (\nu +k\mu +1)k!}=E_{(1/\mu ,1),(\nu +1,1)}^{\left(2\right)}(-z),$$

(57)

again as an example of the Dzrbashjan function. It is an obvious (and was introduced as such by Sir E. Maitland Wright [32]) “fractional index” analog of the classical Bessel function $J_{\nu }\left(z\right)=c(z/2){}_0{F}_1(z^2/4)$, more exactly, of the Bessel–Clifford function $C_{\nu }\left(z\right)$.

Several further “fractional-indices” generalizations of $J_{\nu }\left(z\right)$ and $J_{\nu }^{\mu }\left(z\right)$ are found in the studies of other authors (details are in [59]), and we can represent all of them as multi-index M-L functions. One of them is the so-called generalized Wright–Bessel(–Lommel) functions, introduced by Pathak ([98], 1966),

$$J_{\nu ,\lambda }^{\mu }\left(z\right)={(z/2)}^{\nu +2\lambda }{\displaystyle \sum _{k=0}^{\infty }}\frac{{(-1)}^k{(z/2)}^{2k}}{\Gamma (\nu +k\mu +\lambda +1)\Gamma (\lambda +k+1)}$$

$$={(z/2)}^{\nu +2\lambda }{E}_{(1/\mu ,1),(\nu +\lambda +1,\lambda +1)}^{\left(2\right)}\left(-{(z/2)}^2\right),\mu >0.$$

For $\mu =1$, it includes the mentioned Lommel and Struve functions, e.g., $J_{\nu ,\lambda }^1\left(z\right)=const{S}_{2\lambda +\nu -1,\nu }\left(z\right)$. A next example is the generalized Lommel–Wright function with four indices, introduced by de Oteiza, Kalla and Conde ([99], 1986), with $r>0,n\in N,\nu ,\lambda \in C$:

$$J_{\nu ,\lambda }^{r,n}\left(z\right)={(z/2)}^{\nu +2\lambda }{\displaystyle \sum _{k=0}^{\infty }}\frac{{(-1)}^k{(z/2)}^k}{\Gamma (\nu +kr+\lambda +1)\Gamma {(\lambda +k+1)}^n}$$

$$={(z/2)}^{\nu +2\lambda }{E}_{(1/r,1,\dots ,1),(\nu +\lambda +1,\lambda +1,\dots ,\lambda +1)}^{(n+1)}\left(-{(z/2)}^2\right).$$

5.3. The above is an interesting example of a multi-M-L function with $m=n+1$.

Other particular cases of multi-index ($2m$-parameters) M-L functions with greater multiplicity $m\ge 2$ are:

• For arbitrary $m\ge 2$: let $\forall {\alpha }_i=0$ and $\forall {\beta }_i=1$, $i=1,\dots ,m$. Then, from definition (39), we get again the geometric series

$$E_{(0,0,\dots ,0),(1,1,\dots ,1)}^{\left(m\right)}\left(z\right)={\displaystyle \sum _{k=0}^{\infty }}{z}^k=\frac{1}{1-z}.$$

• Consider the case $m\ge 2$, $\forall {\alpha }_i=1$, $i=1,\dots ,m$. Then, the function

$$E_{(1,1,\dots ,1),({\beta }_1,\dots ,{\beta }_m)}^{\left(m\right)}\left(z\right)={}_1{\Psi }_m\left[\left.\begin{array}{c}(1,1)\\ {({\beta }_i,1)}_1^m\end{array}\right|z\right]={[{\displaystyle \prod _{i=1}^m}\Gamma \left({\mu }_i\right)]}^{-1}{}_1{F}_m\left(1;{\beta }_1,{\beta }_2,\dots ,{\beta }_m;z\right)$$

reduces to ${}_1{F}_m$-function and also to a Meijer’s $G_{1,m+1}^{1,1}$-function. Denote ${\beta }_i={\gamma }_i+1,i=1,\dots ,m$, and let additionally one of the ${\beta }_i$ be 1, e.g., ${\beta }_m=1$, i.e., ${\gamma }_m=0$. Then, the multi-index M-L function becomes a ${}_0{F}_{m-1}$-function, that is, a hyper-Bessel function in the sense of Delerue [89] (see also [9] (Ch.3)):

$$J_{{\gamma }_i,\dots ,{\gamma }_{m-1}}^{(m-1)}\left(z\right)={\left(\frac{z}{m}\right)}^{{\displaystyle \sum _{i=1}^{m-1}}{\gamma }_i}{E}_{(1,1,\dots ,1),({\gamma }_1+1,{\gamma }_2+1,\dots ,{\gamma }_{m-1}+1,1)}^{\left(m\right)}\left(-{\left(\frac{z}{m}\right)}^m\right)$$

(58)

$$={\left[{\displaystyle \prod _{i=1}^{m-1}}\Gamma ({\gamma }_i+1)\right]}^{-1}{\left(\frac{z}{m}\right)}^{{\displaystyle \sum _{i=1}^{m-1}}{\gamma }_i}{}_0{F}_{m-1}\left({\gamma }_1+1,{\gamma }_2+1,\dots ,{\gamma }_{m-1}+1;-{\left(\frac{z}{m}\right)}^m\right)$$

$$:={\left[{\displaystyle \prod _{i=1}^{m-1}}\Gamma ({\gamma }_i+1)\right]}^{-1}{\left(\frac{z}{m}\right)}^{{\displaystyle \sum _{i=1}^{m-1}}{\gamma }_i}{j}_{{\gamma }_1,..,{\gamma }_{m-1}}^{(m-1)}(-z),$$

(59)

where $j_{{\gamma }_1,..,{\gamma }_{m-1}}^{(m-1)}$ is called as normalized hyper-Bessel function.

This representation suggests that the multi-index M-L functions (39) with arbitrary $({\alpha }_1,\dots ,{\alpha }_m)$≠$(1,\dots ,1)$ can be interpreted as fractional-indices analogs of the hyper-Bessel functions (58) and (59), which themselves are multi-index (but integer) analogs of the Bessel function. Functions (58) and (59) are closely related to the hyper-Bessel differential operators (9) (see Section 3.1), and form a fundamental system of solutions of the differential equations of the form $By\left(z\right)=\lambda y\left(z\right)$; the details are found in Kiryakova [9] (Ch.3, Th.3.4.3). For example, if the hyper-Bessel operator (9) is with $\beta =m$, ${\gamma }_1<{\gamma }_2<\dots <{\gamma }_m=0<{\gamma }_1+1$, the solution of the Cauchy problem $By\left(z\right)=-y\left(z\right),y\left(0\right)=1,y^{\left(j\right)}\left(0\right)=0,j=1,\dots ,m-1,$ is given by the normalized hyper-Bessel function (59): $y\left(z\right)=j_{{\gamma }_1,..,{\gamma }_{m-1}}^{(m-1)}(-z)$. Closely related functions are also the Bessel–Clifford functions of order m:

$$C_{{\nu }_1,\dots ,{\nu }_m}\left(z\right)={\displaystyle \sum _{k=0}^{\infty }}\frac{{(-1)}^k{z}^k}{\Gamma ({\nu }_1+k+1)\dots \Gamma ({\nu }_m+k+1)k!}=E_{(1,\dots ,1),({\nu }_1+1,\dots ,{\nu }_m+1,1)}^{(m+1)}(-z).$$

Let us mention the special functions appearing in a very recent paper by Ricci [100]. He considered the so-called Laguerre derivative $D_{ℐ}={\displaystyle \frac{d}{dz}}z{\displaystyle \frac{d}{dz}}$ and its iterates $D_{mℐ}={\displaystyle \frac{d}{dz}}z{\displaystyle \frac{d}{dz}}z\dots {\displaystyle \frac{d}{dz}}z$, same as the particular hyper-Bessel differential operators (19) considered in operational calculus by Ditkin and Prudnikov [50], as mentioned in Section 3.1. Then, the L-exponentials $e_1\left(z\right),e_2\left(z\right),\dots ,e_m\left(z\right),\dots$, which are eigenfunctions of $D_{mL}$, that is, $D_{mL}{e}_m\left(\lambda z\right)=\lambda e_m\left(\lambda z\right)$, are shown in [100] to have the form

$$e_m\left(z\right)={\displaystyle \sum _{k=0}^{\infty }}\frac{z^k}{{(k!)}^{m+1}}={}_0{F}_m(-;1,1,\dots ,1;z)={}_1{\Psi }_{m+1}\left[\left.\begin{array}{c}(1,1)\\ (1,1),(1,1),...,(1,1)\end{array}\right|z\right].$$

(60)

Then, these are examples of the hyper-Bessel functions (58) and of the multi-index Mittag–Leffler functions $E_{(1,\dots ,1),(1,\dots ,1)}^{(m+1)}\left(z\right)$ as well. In [100], applications to population dynamics and in solutions of linear dynamical systems of these SF and of the related Laguerre-type Bell polynomials and Laguerre-type generalized hypergeometric functions are discussed.

• The Rabotnov function (the $\alpha$-exponential function), presented in 5.1., appeared in Rabotnov’s works on application of fractional order operators in mechanics of solids. It is interesting to consider its multi-index analog, that is the case with all ${\beta }_i={\alpha }_i=\alpha >0,i=1,\dots ,m$. This is the function

$$y_{\alpha }^{\left(m\right)}\left(z\right)=z^{\alpha -1}{E}_{(\alpha ,\dots ,\alpha ),(\alpha ,\dots ,\alpha )}^{\left(m\right)}\left(z^{\alpha }\right)=z^{\alpha -1}{\displaystyle \sum _{k=0}^{\infty }}\frac{z^{\alpha k}}{{[\Gamma (\alpha +\alpha k)]}^m}.$$

(61)

Observe that, for $\alpha =1$, we get the Ricci function (60), namely: $e_{m-1}\left(z\right)={\displaystyle \sum _{k=0}^{\infty }}{\displaystyle \frac{z^k}{{[k!]}^m}}$, and also a case of the original Le Roy function with $\gamma =m$.

• In general, for rational values of $\forall {\alpha }_i,i=1,\dots ,m$, the functions (39) are reducible to generalized hypergeometric functions ${}_1{F}_m$ and to Meijer’s G-functions $G_{1,m+1}^{1,1}$, that is, to classical special functions.

Remark 2.

Note that all the results we derived for the multi-index M-L functions can be applied for their particular cases mentioned above.

6. Other Special Cases of the Wright Generalized Hypergeometric Functions ${}_{\mathbf{p}}{\Psi }_{\mathbf{q}}$

6.1.Virchenko and Ricci generalized hypergeometric functions. In [101] and some other papers, Virchenko studied some generalized hypergeometric functions denoted by ${}_2{R}_1^{\tau }\left(z\right)$ and ${}_1{\Phi }_1^{\tau }\left(z\right)$, as well as their integral representations, relations and applications to the generalized Legendre functions $P_k^{m,m}\left(z\right),Q_k^{m,n}\left(z\right)$, gamma functions, Laguerre’s functions, etc.

•

$${}_2{R}_1^{\omega ,\mu }(a,b;c;z)=\frac{\Gamma \left(c\right)}{\Gamma \left(a\right)\Gamma \left(b\right)}{\displaystyle \sum _{k=0}^{\infty }}\frac{\Gamma (a+k)\Gamma (b+\frac{\omega }{\mu }k)}{\Gamma (c+\frac{\omega }{\mu }k)}·\frac{z^k}{k!}.$$

For ${\displaystyle \frac{\omega }{\mu }}:=\tau >0$, and $a,b,c$ - complex, $a+k\ne 0,-1,-2,\dots$; $b+\tau k\ne 0,-1,-2,\dots$, $k=0,1,2,\dots$; $\left|z\right|<1$, it is rewritten as

$${}_2{R}_1^{\tau }(a,b;c;z)=\frac{\Gamma \left(c\right)}{\Gamma \left(a\right)\Gamma \left(b\right)}{\displaystyle \sum _{k=0}^{\infty }}\frac{\Gamma (a+k)\Gamma (b+\tau k)}{\Gamma (c+\tau k)}·\frac{z^k}{k!},$$

which is nothing but the Wright g.h.f. ${\displaystyle \frac{\Gamma \left(c\right)}{\Gamma \left(a\right)\Gamma \left(b\right)}}{}_2{\Psi }_1\left[\left.\begin{array}{c}(a,1),(b,\tau )\\ (c,\tau )\end{array}\right|z\right]$. Virchenko also proposed some examples of elementary functions for these special functions, e.g., ${\left(ln(1+z)\right)}_{\tau }$ and ${\left(arcsinz\right)}_{\tau }$; some generalized incomplete B-function; the Gauss function ${}_2{F}_1$; etc.

•

$${}_1{\Phi }_1^{\tau }(a;c;z)=\frac{\Gamma \left(c\right)}{\Gamma \left(a\right)}{\displaystyle \sum _{k=0}^{\infty }}\frac{\Gamma (a+\tau k)}{\Gamma (c+\tau k)}·\frac{z^k}{k!},$$

and, in Virchenko [101], generalizations of the gamma function, incomplete gamma function, probability integrals and Laguerre’s functions are introduced by means of ${}_1{\Phi }_1^{\tau }\left(z\right)$, which is a Wright g.h.f. of the form ${\displaystyle \frac{\Gamma \left(c\right)}{\Gamma \left(a\right)}}{}_1{\Psi }_1\left[\left.\begin{array}{c}(a,\tau )\\ (c,\tau )\end{array}\right|z\right]$, and, according to our classifications in Section 8, a confluent type g.h.f.

• In 5.3., the recent paper by Ricci [100] is mentioned for the Laguerre-type derivatives and related special functions. Along with the functions (60), there he also considered the Laguerre-type (L-) Bessel functions, L-type Gauss hypergeometric functions and the Laguerre-type generalized hypergeometric functions ${}_{ℐp}{F}_q$. They can be shown to be representable by ${}_p{F}_{q+1}$, thus also as ${}_p{\Psi }_{q+1}$, namely:

$${}_{ℐp}{F}_q(a_1,..,a_p;b_1,\dots ,b_q;z)={\displaystyle \sum _{k=0}^{\infty }}\frac{a_1^{\left(k\right)}\dots a_p^{\left(k\right)}}{b_1^{\left(k\right)}\dots b_q^{\left(k\right)}}·\frac{z^k}{{(k!)}^2}$$

$$={\displaystyle \sum _{k=0}^{\infty }}\frac{a_1^{\left(k\right)}\dots a_p^{\left(k\right)}}{b_1^{\left(k\right)}\dots b_q^{\left(k\right)}{\left(1\right)}^{\left(k\right)}}·\frac{z^k}{k!}={}_p{F}_{q+1}(a_1,\dots ,a_p;b_1,\dots ,b_q,1;z).$$

(62)

6.2.Mainardi-Masina and Paris generalized exponential integrals. In [102], Mainardi and Masina introduced a generalized exponential integral ${\text{Ein}}_{\alpha }\left(z\right)$ by replacing the exponential function in the complementary exponential integral $\text{Ein}\left(z\right)$ by the Mittag–Leffler function $E_{\alpha }\left(z\right)$ and mentioned the physical applications for $0<\alpha <1$ in the studies of the creep features of linear viscoelastic models. In the recent paper [103], Paris made the next step to involve the two-parameter M-L function, namely to consider the generalized exponential integral

$${\text{Ein}}_{\alpha ,\beta }\left(z\right)=z{\displaystyle \sum _{k=0}^{\infty }}\frac{{(-1)}^k{z}^{\alpha k}}{(ak+1)\Gamma (\alpha k+\alpha +\beta )},\text{which for}\beta =1\text{gives}{\text{Ein}}_{\alpha }\left(z\right).$$

(63)

As observed, this function can be seen as a case of the Wright g.h.f. with $p=q=2$, namely

$${\text{Ein}}_{\alpha ,\beta }\left(z\right)=z{\displaystyle \sum _{k=0}^{\infty }}\frac{\Gamma (\alpha k+1)\Gamma (k+1)}{\Gamma (\alpha k+2)\Gamma (\alpha k+\alpha +\beta )}\frac{{(-z^{\alpha })}^k}{k!}=z{}_2{\Psi }_2\left[\left.\begin{array}{c}(1,\alpha ),(1,1)\\ (2,\alpha ),(\alpha +\beta ,\alpha )\end{array}\right|-z^{\alpha }\right].$$

Paris studied in details the asymptotic expansion of (63) for $\left|z\right|\to \infty$. In [102,103], generalized Sine and Cosine integrals are also considered (of the kind mentioned in 5.1.), for example ${\text{Sin}}_{\alpha ,\beta }\left(z\right)={\text{Ein}}_{2\alpha ,\beta -\alpha }\left(z\right)$, with their asymptotics and plots for different values of parameters.

6.3.The so-called k-analogs of special functions. Claims on inventing and studying “new” classes of special functions in several recent papers have been based on the extended notion of the k-Gamma function, $k>0$. However, in all such works, its representation in terms of the classical Gamma-function is explicitly written there, and then is ignored:

$${\Gamma }_k\left(s\right)={\displaystyle \underset{0}{\overset{\infty }{\int }}}exp(-\frac{t^k}{k})t^{s-1}dt=k^{\frac{s}{k}-1}\Gamma \left(\frac{s}{k}\right),s\in ℂ,\text{Re}\left(s\right)>0,$$

(64)

where $\Gamma (.)$ is the classical Gamma-function.

In addition, the k-Pochhammer symbol is used in the next denotations:

$${\left(\lambda \right)}_{\nu ,\kappa }:={\Gamma }_k(\lambda +\nu \kappa )/{\Gamma }_k\left(\lambda \right),\lambda \in ℂ\setminus \left\{0\right\},\nu \in ℂ,\text{with}{\Gamma }_k\text{as in}\left(64\right).$$

(65)

In [104], using the above two definitions, we showed that most of these “new” functions are in fact some known special functions, namely Wright g.h.f. and its cases. For the details of establishing the mentioned relations, see Kiryakova [104]. In addition, in the references lists of [104,105], one can find the particular authors/sources mentioned below.

• A generalized k-Bessel function was introduced by Gehlot ([106], 2014), and studied by Mondal ([107], 2016) and Shaktawat et al. ([108], 2017). It is defined by

$$W_{\nu ,c}^k\left(z\right)={\displaystyle \sum _{n=0}^{\infty }}\frac{{(-c)}^n}{{\Gamma }_k(nk+\nu +k)}·\frac{{(z/2)}^{2n+\frac{\nu }{k}}}{n!},z\in ℂ,k>0,\text{Re}\left(\nu \right)>-1,c\in ℂ.$$

(66)

However, after simple exercise, the function (66) can be represented as a Wright g.h.f. ${}_0{\Psi }_1$, and even as the simpler g.h.f. ${}_0{F}_1$ of the same type as the classical Bessel function:

$$W_{\nu ,c}^k\left(z\right)={(z/2)}^{\frac{\nu }{k}}{\displaystyle \sum _{n=0}^{\infty }}\frac{{[-c{\left(\frac{z}{2}\right)}^2]}^n}{k^{n+1+\left(\frac{\nu }{k}\right)}\Gamma (n+1+\left(\frac{\nu }{k}\right))\Gamma (n+1)}=\dots$$

$$=\frac{{\left(\frac{z}{2}\right)}^{\frac{\nu }{k}}}{k^{1+\left(\frac{\nu }{k}\right)}}{\displaystyle \sum _{n=0}^{\infty }}\frac{{[-\left(\frac{c}{k}\right){\left(\frac{z}{2}\right)}^2]}^n}{\Gamma (1+\left(\frac{\nu }{k}\right)+n.1)\Gamma (1+n.1)}$$

$$=\frac{{\left(\frac{z}{2}\right)}^{\frac{\nu }{k}}}{k^{1+\left(\frac{\nu }{k}\right)}}{}_1{\Psi }_2\left[\left.\begin{array}{c}(1,1)\\ (1+\frac{\nu }{k},1),(1,1)\end{array}\right|-\frac{c}{k}{\left(\frac{z}{2}\right)}^2\right]=\frac{{(z/2)}^{\nu /k}}{k^{1+(\nu /k)}}{}_0{\Psi }_1\left[\left.\begin{array}{c}--\\ (1+\frac{\nu }{k},1)\end{array}\right|-\frac{c}{k}{\left(\frac{z}{2}\right)}^2\right]$$

$$=\frac{{\left(\frac{z}{2}\right)}^{\frac{\nu }{k}}}{k^{1+\left(\frac{\nu }{k}\right)}\Gamma (1+\nu )}{}_0{F}_1\left(-;1+\frac{\nu }{k};-\frac{c}{k}\frac{z^2}{4}\right).$$

(67)

Indeed, if we take $k=1$ and $c=1$, this function reduces to the classical Bessel function: $W_{\nu ,1}^1\left(z\right)=\frac{{(z/2)}^{\nu }}{\Gamma (1+\nu )}{}_0{F}_1\left(-;1+\nu ;-\frac{z^2}{4}\right)$. For $k>0$ and $c=1$ Gehlot [106] used (66) as a solution of a k-Bessel differential equation. Mondal [107] studied some properties of (66) for arbitrary $c\in C$. Shaktawat et al. [108] evaluated the Marichev–Saigo–Maeda (M-S-M) operators of FC

$$I^{a,a^{\prime },b,b^{\prime },c}f\left(z\right)=z^{c-a-a^{\prime }}{\displaystyle \underset{0}{\overset{1}{\int }}}\frac{{(1-\sigma )}^{c-1}}{\Gamma \left(c\right)}{\sigma }^{-a^{\prime }}{F}_3(a,a^{\prime },b,b^{\prime };c;1-\sigma ,1-\frac{1}{\sigma })f\left(z\sigma \right)d\sigma$$

(68)

of this function. Since its kernel Appel $F_3$-function is a H-function (36) with $m=3$, in view of author’s result from Corollary 3 in Section 7, it is well expected that the result appears in terms of a ${}_3{\Psi }_4$-function (because the indices of ${}_0{\Psi }_1$ are increased by 3 under the 3-tuple FC integral).

• Generalized k-Mittag–Leffler function. It was studied by many authors, for example in its simplest case by Gupta and Parihar ([109], 2014) in the form

$$E_{k,\alpha ,\beta }\left(z\right)={\displaystyle \sum _{n=0}^{\infty }}\frac{z^n}{{\Gamma }_k(\alpha n+\beta )}.$$

This function has various further extensions, such as the generalized k-Mittag–Leffler function by Nisar–Eata–Dhaifalla–Choi ([110], 2016):

$$E_{\kappa ,\alpha ,\beta }^{\eta ,\delta ,p,q}\left(z\right)={\displaystyle \sum _{n=0}^{\infty }}\frac{{\left(\eta \right)}_{qn,\kappa }}{{\Gamma }_k(\alpha n+\beta ){\left(\delta \right)}_{pn,\kappa }}{z}^n,\text{with}\kappa ,p,q\in {ℝ}_{+};\alpha ,\beta ,\eta ,\delta \in ℂ,$$

(69)

$$\text{and}\text{min}\left\{\text{Re}\right(\alpha ),\text{Re}(\beta ),\text{Re}(\eta ),\text{Re}(\delta \left)\right\}>0;q\le \text{Re}\left(\alpha \right)+p.$$

Again, by using (64) and (65), it can be transformed into a Wright g.h.f. (see [104], Case 5.2), namely:

$$E_{\kappa ,\alpha ,\beta }^{\eta ,\delta ,p,q}\left(z\right)=k^{1-\frac{\beta }{k}}\frac{\Gamma (\delta /k)}{\Gamma (\eta /k)}{}_2{\Psi }_2\left[\left.\begin{array}{c}(\frac{\eta }{k},\frac{q\kappa }{k}),(1,1)\\ (\frac{\delta }{k},\frac{p\kappa }{k}),(\frac{\beta }{k},\frac{\alpha }{k})\end{array}\right|k^{\frac{(q-p)\kappa -\alpha }{k}}z\right].$$

Nisar–Eata–Dhaifalla-Choi [110] put efforts to evaluate FC operators’ images of (69) by the standard techniques, and as expected in view of the general results in next Section 7 Theorem 5, Corollarys 1–3) these appear in terms of ${}_5{\Psi }_5$-functions (for the M-S-M operators (68)), in particular, as ${}_4{\Psi }_4$-functions (for the Saigo operators (78)) and ${}_3{\Psi }_3$-functions (for the R-L and E-K operators). In addition, the pathway integrals (that are related to E-K integrals) are calculated there.

• The generalized multi-index Bessel function. In a series of papers, Nisar et al. ([111], 2017, 2019) introduced and studied the function

$$J_{{\left({\beta }_j\right)}_m,\kappa ,b}^{{\left({\alpha }_j\right)}_m,\gamma ,c}\left(z\right)={\displaystyle \sum _{k=0}^{\infty }}\frac{c^k{\left(\gamma \right)}_{\kappa k}}{{\displaystyle \prod _{j=1}^m}\Gamma ({\alpha }_{j}k+{\beta }_j+\frac{b+1}{2})}\frac{z^k}{k!},m=1,2,3,\dots ,$$

(70)

with the Pochhammer symbol denotation (65) for ${\left(\gamma \right)}_{\kappa k}$; and for ${\alpha }_j,{\beta }_j,\gamma ,b,c\in C$, $j=1,2,\dots ,m$; ${\displaystyle \sum _{j=1}^m}\text{Re}\left({\alpha }_j\right)>\text{max}\{0,\text{Re}\left(\kappa \right)-1\}$; $\kappa >0$, $\text{Re}\left({\beta }_j\right)>0$, $\text{Re}\left(\gamma \right)>0$. As shown in Kiryakova [104], this is only a very special case of the Wright generalized hypergeometric function ${}_1{\Psi }_m$, namely:

$$J_{{\left({\beta }_j\right)}_m,\kappa ,b}^{{\left(\alpha \right)}_m,\gamma ,c}\left(z\right)=\frac{1}{\Gamma \left(\gamma \right)}{\displaystyle \sum _{n=0}^{\infty }}\frac{\Gamma (\varkappa n+\gamma )}{{\displaystyle \prod _{j=1}^m}\Gamma \left({\alpha }_{j}n+({\beta }_j+\frac{b+1}{2})\right)}\frac{{\left(cz\right)}^n}{n!}$$

(71)

$$=\frac{1}{\Gamma \left(\gamma \right)}{}_1{\Psi }_m\left[\begin{array}{c}(\gamma ,\varkappa )\\ {({\beta }_j+\frac{b+1}{2},{\alpha }_j)}_{j=1}^m\end{array};cz\right]$$

$=\frac{1}{\Gamma \left(\gamma \right)}{H}_{1,m+1}^{1,1}\left[-cz\left|\begin{array}{c}(1-\gamma ,\varkappa )\\ (0,1),{(1-{\beta }_j-\frac{b+1}{2})}_{j=1}^m\end{array}\right.\right]$, that is, it is also a Fox H-function.

Then, the R-L fractional integral (21) can be evaluated as part of Kiryakova’s general results in next Section 7 (Theorem 5, in particular Corollary 1 for $m=1$, $\gamma =\beta =1$), or directly from Kilbas’ Theorem 2 in [33], which is a variant of Lemma 1 in Kiryakova [112]. Taking there $p=1,q=m$, $c_1=\gamma ,C_1=\kappa$, $d_j={\beta }_j+\frac{b+1}{2},D_j={\alpha }_j$ and $\mu =1$, one obtains the following R-L image for the multi-index Bessel function (70):

$$I^{\lambda }\left\{t^{\delta -1}{J}_{{\left({\beta }_j\right)}_m,\varkappa ,b}^{{\left({\alpha }_j\right)}_m,\gamma ,c}\left(z\right)\right\}=\frac{1}{\Gamma \left(\gamma \right)}{z}^{\delta +\lambda -1}{}_2{\Psi }_{m+1}\left[\begin{array}{c}(\gamma ,\varkappa ),(\delta ,1)\\ {({\beta }_j+\frac{b+1}{2},{\alpha }_j)}_1^m,(\lambda +\delta ,1)\end{array};cz\right].$$

This was to be the result in Theorem 1, Equation (2.4) in arXiv:1706.08039 [111], its v1: 2017, but was written wrongly there—similarly looking but involving a ${}_2{\Psi }_2$-function. The evident true result involves the Wright function ${}_2{\Psi }_{m+1}$ (see Kiryakova [104] (5.3.)), as later corrected in v2: 2019 of [111].

• A special case of (70) appears as a kind of generalized multi-index Mittag–Leffler function. It was introduced by Saxena and Nishimoto ([113], 2010). As mentioned by Agarwal–Rogosin–Trujillo ([114], 2015), it is representable also as a Wright g.h.f. ${}_1{\Psi }_m$, namely:

$$E_{{({\alpha }_j,{\beta }_j)}_m}^{(\gamma ,\kappa )}\left(z\right)={\displaystyle \sum _{n=0}^{\infty }}\frac{{\left(\gamma \right)}_{\kappa n}}{{\displaystyle \prod _{j=1}^m}\Gamma ({\beta }_j+{\alpha }_{j}n)}·\frac{z^n}{n!}=\frac{1}{\Gamma \left(\gamma \right)}{}_1{\Psi }_m\left[\left.\begin{array}{c}(\gamma ,\kappa )\\ {({\beta }_j,{\alpha }_j)}_1^m\end{array}\right|z\right].$$

(72)

Therefore, all the GFC operators (say the R-L, E-K, Saigo, M-S-M operators) of this special function can be evaluated by means of the general results in Section 7, Corollaries 1–3 there. For $m=1,b=-1$, this is the SF considered by Srivastava and Tomovski ([115], 2009): ${\displaystyle E_{\alpha ,\beta }^{\gamma ,\kappa }\left(z\right)={\displaystyle \sum _{n=0}^{\infty }}\frac{{\left(\gamma \right)}_{\kappa n}}{\Gamma (\alpha n+\beta )}·\frac{z^n}{n!}}$.

• Similar, but simpler, is the case of the generalized Lommel-Wright function from the paper by Agarwal–Jain–Agarwal–Baleanu ([116], 2018), which is commented in Kiryakova [117]. It has a representation as a Wright g.h.f. as follows:

$$J_{\omega ,\theta }^{\phi ,m}\left(z\right)={\left(\frac{z}{2}\right)}^{\omega +2\theta }{\displaystyle \sum _{k=0}^{\infty }}\frac{{(-1)}^k{\left(\frac{z}{2}\right)}^{2k}}{{\left(\Gamma (\theta +k+1)\right)}^m\Gamma (\omega +k\phi +\theta +1)}$$

(73)

$$={\left(\frac{z}{2}\right)}^{\omega +2\theta }{}_1{\Psi }_{m+1}\left[(1,1);(\theta +1,1),\dots ,(\theta +1,1),(\omega +\theta +1,\phi );-z^2/4\right],\phi >0.$$

Note, additionally, that (73) is an example of the multi-index Mittag–Leffler function (39), namely: $J_{\omega ,\theta }^{\phi ,m}\left(z\right)={\left(\frac{z}{2}\right)}^{\omega +2\theta }{\left(\frac{z}{2}\right)}^{\omega +2\theta }{E}_{(1,\dots ,1,\phi ),(\theta +1,\dots ,\theta +1,\omega +\theta +1)}^{(m+1)}\left(-{\left(\frac{z}{2}\right)}^2\right)$. Then, all the FC images of (73) evaluated in the commented paper follow at once by our general results (see details in [117]).

6.4.The S-function. It was introduced by Saxena-Daiya ([118], 2015) as a “new” special function extending the M-L function ($p=q=0,k=1$), the Prabhakar function (38), the M-series (76) of Sharma and Jain ([119], 2009) with $\gamma =1$, $k=1$, etc., as follows:

$$S\left[z\right]:=S_{(p,q)}^{\alpha ,\beta ,\gamma ,\tau ,k}(a_1,\dots ,a_p;b_1,\dots ,b_q;z)={\displaystyle \sum _{n=0}^{\infty }}\frac{{\left(a_1\right)}_n\dots {\left(a_p\right)}_n·{\left(\gamma \right)}_{n\tau ,k}}{{\left(b_1\right)}_n\dots {\left(b_q\right)}_n·{\Gamma }_k(n\alpha +\beta )}\frac{z^n}{n!},$$

(74)

$$\text{with}k\in ℝ;\alpha ,\beta ,\gamma ,\tau \in ࠶;\text{Re}\left(\alpha \right)>0;\text{Re}\left(\alpha \right)>k\text{Re}\left(\tau \right),p<q+1.$$

For $p=q=0$, it reduces to the generalized k-Mittag–Leffler function $E_{k,\alpha ,\beta }^{\gamma ,\tau }\left(z\right)$, a variant of (69). However, it can be easily seen to be special case of the generalized hypergeometric function of Wright of the form ${}_{p+1}{\Psi }_{q+1}$. Unfortunately, this fact has not been observed, neither by the authors introducing (74) nor by their numerous followers. Namely, one can write (74) as follows (see details in [104]):

$$S\left[z\right]=k^{1-\frac{\beta }{k}}\frac{\Gamma \left(b_1\right)\dots \Gamma \left(b_q\right)}{\Gamma \left(a_1\right)\dots \Gamma \left(a_p\right)·\Gamma \left(\frac{\gamma }{k}\right)}{}_{p+1}{\Psi }_{q+1}\left[\begin{array}{c}(a_1,1),...,(a_p,1),(\frac{\gamma }{k},\tau )\\ (b_1,1),...,(b_q,1),(\frac{\beta }{k},\frac{\alpha }{k})\end{array};z{k}^{\tau -\frac{\alpha }{k}}\right].$$

That is, the “new” special function $S\left[z\right]$ is nothing but a case of the Wright function ${}_{p+1}{\Psi }_{q+1}\left(z{k}^{\tau -\frac{\alpha }{k}}\right)$. Then, all results for images of FC operators, as R-L, E-K, Saigo, M-S-M and the Euler-transform, follow from the statements in Section 7.

• Special cases of the S-function in 6.4. are the generalized K-series and the M-series. Recently, (K.) Sharma ([120], 2012) introduced an extension of both g.h.f. ${}_p{F}_q\left(z\right)$ and Prabhakar function $E_{\alpha ,\beta }^{\gamma }\left(z\right)$:

$${}_p{\stackrel{\alpha ,\beta ;\gamma }{K}}_q\left(a_1,\dots ,a_p;b_1,\dots ,b_q;z\right)={\displaystyle \sum _{n=0}^{\infty }}\frac{{\left(a_1\right)}_n\dots {\left(a_p\right)}_n}{{\left(b_1\right)}_n\dots {\left(b_q\right)}_n}\frac{{\left(\gamma \right)}_n{z}^n}{\Gamma (\alpha n+\beta )},z,\alpha ,\beta \in ࠶,\text{Re}\alpha >0,$$

(75)

with integers $p\le q$ (and additional requirement $\left|z\right|<R={\alpha }^{\alpha }$ if $p=q+1$). For $\gamma =1$, this gives the M-series (76) of (M.) Sharma and Jain ([119], 2009):

$${}_p{\stackrel{\alpha ,\beta }{M}}_q\left(a_1,\dots ,a_p;b_1,\dots ,b_q;z\right):={}_p{\stackrel{\alpha ,\beta }{M}}_q\left(z\right)={\displaystyle \sum _{n=0}^{\infty }}\frac{{\left(a_1\right)}_n\dots {\left(a_p\right)}_n}{{\left(b_1\right)}_n\dots {\left(b_q\right)}_n}\frac{z^n}{\Gamma (\alpha n+\beta )}$$

$$=\kappa {}_{p+1}{\Psi }_{q+1}\left[\left.\begin{array}{c}(a_1,1),\dots (a_p,1),(1,1)\\ (b_1,1),\dots ,(b_q,1),(\beta ,\alpha )\end{array}\right|z\right],\text{where}\kappa ={\displaystyle \prod _{j=1}^q}\Gamma \left(b_j\right)/{\displaystyle \prod _{i=1}^p}\Gamma \left(a_i\right).$$

(76)

We can mention its particular cases, for example: (1) for $\beta =1$, the (simpler) M-series, introduced by M. Sharma (2008); (2) for $p=q=0$ (no upper and lower parameters), M-L function $E_{\alpha ,\beta }\left(z\right)$; (3) for $p=0,q=1,b_1=1$, the Wright function $\varphi (\alpha ,\beta ;z)$, or the generalized Bessel-Maitland function (57); (4) for $p=q=1$, $a_1=\gamma ,b_1=1$ in (75), the Prabhakar type function (38); and (5) for $\alpha =\beta =1$, the g.h.f. ${}_p{F}_q(a_1,\dots ,a_p;b_1,\dots ,b_q;z)$.

In the recent arXiv preprint [121], Lavault represented (75) as a Wright g.h.f.:

$${}_p{\stackrel{\alpha ,\beta ;\gamma }{K}}_q\left(a_1,\dots ,a_p;b_1,\dots ,b_q;z\right):={}_p{\stackrel{\alpha ,\beta ;\gamma }{K}}_q\left(z\right)$$

$$=\frac{{\displaystyle \prod _{j=1}^q}\Gamma \left(b_j\right)}{\Gamma \left(\gamma \right){\displaystyle \prod _{i=1}^p}\Gamma \left(a_i\right)}{}_{p+2}{\Psi }_{q+2}\left[\left.\begin{array}{c}(a_1,1),...,(a_p,1),(\gamma ,1),(1,1)\\ (b_1,1),...,(b_q,1),(1,1),(\beta ,\alpha )\end{array}\right|z\right],$$

(77)

although this can also be reduced to: $={\displaystyle \frac{{\displaystyle \prod _{j=1}^q}\Gamma \left(b_j\right)}{\Gamma \left(\gamma \right){\displaystyle \prod _{i=1}^p}\Gamma \left(a_i\right)}}{}_{p+1}{\Psi }_{q+1}\left[\left.\begin{array}{c}(a_1,1),\dots (a_p,1),(\gamma ,1))\\ (b_1,1),\dots ,(b_q,1),(\beta ,\alpha )\end{array}\right|z\right]$, since the two pairs $(1,1)$ of parameters in the upper and low rows eliminate each other. In [121] some FC operators of this K-series are calculated, as the R-L, Saigo and M-S-M operators. Naturally, a R-L integral is transforming a ${}_p{\stackrel{\alpha ,\beta ;\gamma }{K}}_q$-function into a ${}_{p+1}{\stackrel{\alpha ,\beta ;\gamma }{K}}_{q+1}$-function (Theorem 4.1, there), similarly to our Example 11 in [112] for the M-series. Next, in Theorem 4.2 of [121] for the M-series and Corollary 4.3 for the K-series, the Saigo operator (78) (with Gauss hypergeometric function (35), GFC with $m=2$) is derived,

$$I^{\alpha ,\beta ,\eta }f(z)=\frac{z^{-\alpha -\beta }}{\Gamma (\alpha )}{\displaystyle \underset{0}{\overset{z}{\int }}}{(z-\xi )}^{\alpha -1}{}_2{F}_1(\alpha +\beta ,-\eta ;\alpha ;1-\frac{\xi }{z})f\left(\xi \right)d\xi$$

$$=\frac{z^{-\beta }}{\Gamma \left(\alpha \right)}{\displaystyle \underset{0}{\overset{1}{\int }}}{(1-\sigma )}^{\alpha -1}{}_2{F}_1(\alpha +\beta ,-\eta ;\alpha ;1-\sigma )f\left(z\sigma \right)d\sigma .$$

(78)

Since the K-series (75) is a ${}_{p+1}{\Psi }_{q+1}$-function, from our results (and Corollary 3 [112]; see also Corollary 2 in the next section), it is expected that the result should be given as a ${}_{p+3}{\Psi }_{q+3}$-function (the indices are to be increased by 2), which is the result (4.10) in [121]:

$$I^{\alpha ,\beta ,\gamma }\left\{t^{\sigma -1}{}_p{\stackrel{\xi ,\eta ;\nu }{K}}_q\left(c{z}^{\mu }\right)\right\}$$

$$=\frac{{\prod }_1^q\Gamma \left(b_j\right)}{{\prod }_1^p\Gamma \left(a_i\right)}\frac{z^{\sigma -\beta -1}}{\Gamma \left(\nu \right)}{}_{p+3}{\Psi }_{q+3}\left[\left.\begin{array}{c}{(a_i,1)}_1^p,(\sigma ,\mu ),(-\beta +\gamma +\sigma ,\mu ),(\nu ,1)\\ {\left(b_j\right)}_1^q,(\beta +\sigma ,\mu ),(\alpha +\gamma +\sigma ,\mu ),(\eta ,\xi )\end{array}\right|c{z}^{\mu }\right].$$

Similarly, the M-S-M-images (68) follow as ${}_{p+4}{\Psi }_{q+4}$-functions, according to Corollary 3 in next section.

6.5.k-Wright generalized hypergeometric function ${}_p{\Psi }_q^k$. Purohit and Badguzer ([122], 2018) introduced the generalized k-Wright function, as a k-extension ($k>0$) of the Wright g.h.f. (4), by

$${}_p{\Psi }_q^k\left(z\right)={}_p{\Psi }_q^k\left[\left.\begin{array}{c}(a_1,A_1),\dots ,(a_p,A_p)\\ (b_1,B_1),\dots ,(b_q,B_q)\end{array}\right|z\right]={\displaystyle \sum _{n=0}^{\infty }}\frac{{\Gamma }_k(a_1+n{A}_1)\dots {\Gamma }_k(a_p+n{A}_p)}{{\Gamma }_k(b_1+n{B}_1)\dots {\Gamma }_k(b_q+n{B}_q)}\frac{z^n}{n!}.$$

(79)

Replacing the k-Gamma function by the classical Gamma function according to (64), it is seen that the “new” function is again a Wright generalized hypergeometric function, of the form

$$\text{const}{}_{p+1}{\Psi }_{q+1}\left[\left.\begin{array}{c}{(a_i/k,A_i/k)}_{i=1}^p\\ {(b_j/k,B_j/k)}_{j=1}^q\end{array}\right|k^{(A_1+\dots +A_p-B_1-\dots -B_q)/k}·z\right].$$

(80)

7. Results for the FC and GFC Images of SF of FC

Recently, there have appeared too many papers that deal with evaluation of FC and GFC operators of various special functions. They use the same standard techniques—replace the particular function by its power series, then interchange the orders of integration (fractional order integrals) and summation, etc. Usually only the special functions are changed and also the FC operators—with more and more general ones (but all these happen to be cases of our GFC operators). The great number of combinations “special function + particular operator” explains the dramatically increasing production of such works.

Based on our older results on GFC for SF, since the work in [9], in the papers [64,104,112,117,123], and in a recent survey paper [105] in this same journal, we propose an unified approach how this job can be done at once, for all SF of FC (we mean the H- and G-functions and in particular the Wright g.h.f., multi-index M-L functions and all their particular cases) and for all operators of GFC (we mean the generalized fractional integrals and derivatives of the form (25) and (28), thus including the R-L, E-K, Saigo, Marichev–Saigo–Maeda operators, etc.). For the initiating idea, we need to pay tribute to the initial classical results of 20th century in the Bateman Project on Integral Transforms [7] and in works by Askey [2], Lavoie–Osler–Tremblay [124], etc. for the R-L images of many elementary functions and of the simplest ${}_p{F}_q$-functions, as: ${}_0{F}_1$, ${}_1{F}_1$ and ${}_1{F}_2$. We combined these with the composition/decomposition rule (26) presenting the GFC operators as compositions of weighted R-L/E-K operators. As a recent survey on FC images of elementary functions, we mention also the work of Garrappa–Kaslik–Popolizio [125].

Below, we remind only the statements of the main results from the mentioned author’s papers, as surveyed in [105], in this same journal.

Theorem 4.

The $I_{\left({\beta }_k\right),m}^{\left({\gamma }_k\right),\left({\delta }_k\right)}$-image (23) of a H-function is also a H-function whose last three components of the order are increased by m (the multiplicity in GFC operators), and with additional parameters depending on those of the generalized fractional integration. Namely,

$$I_{\left({\beta }_k\right),m}^{\left({\gamma }_k\right),\left({\delta }_k\right)}\left\{H_{u,v}^{s,t}\left[\lambda z\left|\begin{array}{c}{(c_i,C_i)}_1^u\\ {(d_j,D_j)}_1^v\end{array}\right.\right]\right\}=H_{u+m,v+m}^{s,t+m}\left[\lambda z\left|\begin{array}{c}{(c_i,C_i)}_1^t,{(-{\gamma }_k)}_1^m,{(c_i,C_i)}_{t+1}^u\\ {(d_j,D_j)}_1^s,{(-{\gamma }_k-{\delta }_k)}_1^m,{(d_j,D_j)}_{s+1}^v\end{array}\right.\right].$$

(81)

Then, GFC images of almost all SF of FC can be evaluated from (81). This result is based on a formula for the integral of product of two arbitrary H-functions, namely for the Mellin transform of such a product ([9] (App., (E.21${}^{\prime }$), [12] (§5.1, (5.1.1)), [14] (§2.25, (1))). A similar formula presents the GFC operators (with $G_{m,m}^{m,0}$-kernel) of arbitrary G-function, in terms of another G-function with increased orders and additional parameters ((Lemma 1.2.2 in [9] and Corollary 1 in [105]).

Since most of the considered SF of FC are Wright g.h.f., the main and most useful result is as follows.

Theorem 5.

The image of a Wright g.h.f. ${}_p{\Psi }_q\left(z\right)$ by a generalized fractional integral (23) (multiple, m-tuple Erdélyi-Kober integral), provided ${\delta }_k\ge 0,{\gamma }_k>-1,k=1,\dots ,m$, $c>-1$, $\mu >0,\lambda \ne 0$, is another Wright g.h.f. with indices p and q increased by the multiplicity m and additional parameters related to these of the GFC integral:

$$I_{{\left({\beta }_k\right)}_1^m,m}^{{\left({\gamma }_k\right)}_1^m,{\left({\delta }_k\right)}_1^m}\left\{z^c{}_p{\Psi }_q\left[\left.\begin{array}{c}(a_1,A_1),\dots ,(a_p,A_p)\\ (b_1,B_1),\dots ,(b_q,B_q)\end{array}\right|\lambda z^{\mu }\right]\right\}$$

$$=z^c{}_{p+m}{\Psi }_{q+m}\left[\left.\begin{array}{c}{(a_i,A_i)}_1^p,{({\gamma }_k+1+\frac{c}{{\beta }_k},\frac{\mu }{{\beta }_k})}_1^m\\ {(b_j,B_j)}_1^q,{({\gamma }_k+{\delta }_k+1+\frac{c}{{\beta }_k},\frac{\mu }{{\beta }_k})}_1^m\end{array}\right|\lambda z^{\mu }\right].$$

(82)

Specially, for $c=0$, $\mu =1$, this result is simplified to ${}_p{\Psi }_q\left(\lambda z\right)⟼{}_{p+m}{\Psi }_{q+m}\left(\lambda z\right)$, as above.

Similarly (Theorem 4.2 in [104]; Theorem 4 in [105]),

$$D_{{\left({\beta }_k\right)}_1^m,m}^{{\left({\gamma }_k\right)}_1^m,{\left({\delta }_k\right)}_1^m}\left\{z^c{}_p{\Psi }_q\left[\left.\begin{array}{c}(a_1,A_1),\dots ,(a_p,A_p)\\ (b_1,B_1),\dots ,(b_q,B_q)\end{array}\right|\lambda z^{\mu }\right]\right\}$$

$$=z^c{}_{p+m}{\Psi }_{q+m}\left[\left.\begin{array}{c}{(a_i,A_i)}_1^p,{({\gamma }_k+{\delta }_k+1+\frac{c}{{\beta }_k},\frac{\mu }{{\beta }_k})}_1^m\\ {(b_j,B_j)}_1^q,{({\gamma }_k+1+\frac{c}{{\beta }_k},\frac{\mu }{{\beta }_k})}_1^m\end{array}\right|\lambda z^{\mu }\right].$$

(83)

The simpler results for the ${}_p{F}_q$-functions read by analogy (Corollarys 4.1 and 4.2 in [104]), for example with $\beta =1$, as:

$$I_{1,m}^{{\left({\gamma }_k\right)}_1^m,{\left({\delta }_k\right)}_1^m}\left\{z^c{}_p{F}_q\left(a_1,\dots ,a_p;b_1,\dots ,b_q;\lambda z\right)\right\}$$

$$=\left[{\displaystyle \prod _{k=1}^m}\frac{\Gamma ({\gamma }_k+c+1)}{\Gamma ({\gamma }_k+{\delta }_k+c+1)}\right]z^c{}_{p+m}{F}_{q+m}\left(a_i,\dots ,a_p,{({\gamma }_k+c+1)}_1^m;b_1,\dots ,b_q,{({\gamma }_k+{\delta }_k+c+1)}_1^m;\lambda z\right).$$

(84)

We also describe the corollaries of the results (82) and (83) for the particular cases of most often FC operators on which the other authors have exercised their evaluations, say for: $m=1$ (R-L and E-K), $m=2$ (Saigo operators) and $m=3$ (M-S-M operators). These results for arbitrary Wright g.h.f. are mentioned below.

Corollary 1.

For the Riemann–Liouville (R-L) integrals and derivatives, the simplest results are parts of Lemmas 1 and 2 in Kiryakova [105]:

$$R^{\delta }\left\{z^c{}_p{\Psi }_q\left[\left.\begin{array}{c}(a_1,A_1),\dots ,(a_p,A_p)\\ (b_1,B_1),\dots ,(b_q,B_q)\end{array}\right|\lambda z^{\mu }\right]\right\}=z^{c+\delta }{}_{p+1}{\Psi }_{q+1}\left[\left.\begin{array}{c}{(a_i,A_i)}_1^p,(c+1,\mu )\\ {(b_j,B_j)}_1^q,(c+\delta +1,\mu )\end{array}\right|\lambda z^{\mu }\right],$$

(85)

$$D^{\delta }\left\{z^c{}_p{\Psi }_q\left[\left.\begin{array}{c}(a_1,A_1),\dots ,(a_p,A_p)\\ (b_1,B_1),\dots ,(b_q,B_q)\end{array}\right|\lambda z^{\mu }\right]\right\}=z^{c-\delta }{}_{p+1}{\Psi }_{q+1}\left[\left.\begin{array}{c}{(a_i,A_i)}_1^p,(c+1,\mu )\\ {(b_j,B_j)}_1^q,(c+1-\delta ,\mu )\end{array}\right|\lambda z^{\mu }\right].$$

(86)

The results for the E-K operators have same expressions as in (82) and (83) with $m=1$.

Corollary 2.

The images of the Wright g.h.f. ${}_p{\Psi }_q$ and, in particular, of the g.h.f. ${}_p{F}_q$ under the Saigo operators (78) are given by the formulas:

$$I^{\alpha ,\beta ,\eta }\left\{z^c{}_p{\Psi }_q\left[\left.\begin{array}{c}{(a_i,A_i)}_1^p\\ {(b_j,B_j)}_1^q\end{array}\right|\lambda z^{\mu }\right]\right\}=z^{c-\beta }{}_{p+2}{\Psi }_{q+2}\left[\left.\begin{array}{c}{(a_i,A_i)}_1^p,(\eta -\beta +1+c,\mu ),(1+c,\mu )\\ {(b_j,B_j)}_1^q,(-\beta +1+c,\mu ),(\alpha +\eta +1+c,\mu )\end{array}\right|\lambda z^{\mu }\right],$$

(87)

(for $c=0$, $\mu =1$, this is Corollary 3 in [112]) and

$$I^{\alpha ,\beta ,\eta }\left\{{}_p{F}_q\left(a_1,...,a_p;b_1,...,b_q;\lambda z\right)\right\}=z^{-\beta }{}_{p+2}{F}_{q+2}\left(a_1,\dots ,a_p,\eta -\beta +1,1;b_1,\dots ,b_q,-\beta +1,\alpha +\eta +1;\lambda z\right).$$

(88)

Corollary 3.

The Marichev–Saigo–Maeda (M-S-M) operators (68) transform a Wright g.h.f. function into same kind of special function but with indices increased by 3:

$$I^{a,a^{\prime },b,b^{\prime },c}\left\{{}_p{\Psi }_q\left[\left.\begin{array}{c}{(a_i,A_i)}_1^p\\ {(b_j,B_j)}_1^q\end{array}\right|\lambda z^{\mu }\right]\right\}$$

$$=z^{c-a-a^{\prime }}{}_{p+3}{\Psi }_{q+3}\left[\left.\begin{array}{c}{(a_i,A_i)}_1^p,(a-a^{\prime }+1,1),(b-a^{\prime }+1,1),(c-2a^{\prime }-b^{\prime }+1,1)\\ {(b_j,B_j)}_1^q,(a-a^{\prime }+b+1,1),(c-2a^{\prime }+1,1),(c-a^{\prime }-b^{\prime }+1,1)\end{array}\right|\lambda z^{\mu }\right].$$

(89)

We state here also the more general result for images of arbitrary Wright generalized hypergeometric function in the case of multiple Wright–Erdélyi–Kober operators (33).

Theorem 6.

(Kiryakova, [60], Theorem 9) The image of a Wright generalized function ${}_p{\Psi }_q\left(z\right)$ by a multiple W-E-K operator (33) has the form

$$I_{\left({\beta }_k\right),\left({\lambda }_k\right),m}^{\left({\gamma }_k\right),\left({\delta }_k\right)}\left\{{}_p{\Psi }_q\left[\left.\begin{array}{c}(a_1,A_1),\dots ,(a_p,A_p)\\ (b_1,B_1),\dots ,(b_q,B_q)\end{array}\right|z\right]\right\}={}_{p+m}{\Psi }_{q+m}\left[\left.\begin{array}{c}{(a_j,A_j)}_1^p;{({\gamma }_k+1,1/{\lambda }_k)}_1^m\\ {(b_k,B_k)}_1^q;{({\gamma }_k+{\delta }_k+1,1/{\beta }_k)}_1^m\end{array}\right|z\right].$$

(90)

Conversely, the alternatively stated result reads as: each ${}_{p+m}{\Psi }_{q+m}$-function can be represented by means of a multiple (m-tuple) operator $\tilde{I}$ of GFC, of a ${}_p{\Psi }_q$-function, the orders of which are reduced by m:

$${}_{p+m}{\Psi }_{q+m}\left[\left.\begin{array}{c}{(a_j,A_j)}_{j=1}^p;{(a_{p+i},A_{p+i})}_{i=1}^m\\ {(b_k,B_k)}_{k=1}^q;{(b_{q+i},B_{q+i})}_{i=1}^m\end{array}\right|z\right]=\tilde{I}\left\{{}_p{\Psi }_q\left[\left.\begin{array}{c}{(a_j,A_j)}_{j=1}^p\\ {(b_k,B_k)}_{k=1}^q\end{array}\right|z\right]\right\},$$

(91)

with

$$\tilde{I}f\left(z\right)=I_{{(1/B_{q+1})}_{i=1}^m,{(1/A_{p+1})}_{i=1}^m,m}^{{(a_{p+i}-1)}_{i=1}^m,{(b_{q+i}-a_{p+i})}_{i=1}^m}f\left(z\right)oftheform\left(33\right).$$

A long list of examples how these general results work at once for any of the SF of FC mentioned in previous sections is provided in author’s works [104,105,112,117,123], including some of the particular cases of W.g.h.f. and of multi-index M-L f., mentioned in Section 5 and Section 6. There we also provided the details on the references items for the authors cited here only with years.

8. Theory of SF of FC in View of GFC Operators

Usually, the special functions of mathematical physics are defined by means of power series representations. However, some alternative representations can be used as their definitions. Let us mention the well-known Poisson integral (52) for the Bessel function and the analytical continuation of the Gauss hypergeometric function via the Euler integral formula. The Rodrigues differential formulas, involving repeated or fractional differentiation are also used as definitions of the classical orthogonal polynomials and their generalizations. As to the other special functions (most of them being ${}_p{F}_q$- and ${}_p{\Psi }_q$-functions), such representations have been less popular and even unknown in the general case. There exist various integral and differential formulas, but, unfortunately, quite peculiar for each corresponding special function and scattered in the literature without any common idea to relate them.

In our works since 1985 (e.g., [9] (Ch.4), [58,60]), we showed that all the classical SF and the SF of FC (in the sense of generalized hypergeometric functions ${}_p{F}_q$ and ${}_p{\Psi }_q$) can be presented by means of generalized fractional integrals or derivatives of three basic elementary functions. On this basis, these special functions have been classified into three specific classes, and several new integral and differential representations have been proposed under a unified idea. Besides, for these three classes of SF, we provide analogs of the mentioned Poisson and Euler integral formulas and of the Rodrigues differential formulas, which can also be used for alternative definitions of these special functions, their analytical extensions or for numerical algorithms.

The idea is briefly explained as follows: (i) most of the classical SF (SF of mathematical physics) and SF of FC are nothing but modifications of the g.h.f. ${}_p{F}_q$ or ${}_p{\Psi }_q$; (ii) each ${}_p{F}_q$-function or ${}_p{\Psi }_q$-function can be represented as an E-K fractional differintegral (i.e., integral or derivative) of a ${}_{p-1}{F}_{q-1}$-function or ${}_{p-1}{\Psi }_{q-1}$, respectively; (iii) a finite number of steps (ii) leads to one of the basic g.h.f. (${}_0{F}_{q-p}$ (for $q-p=1$: Bessel function); ${}_1{F}_1$ (confluent h.f.) and ${}_0{F}_0$ (exponent); and ${}_2{F}_1$ (Gauss h.f.) and ${}_1{F}_0$ (beta-distribution) to the simplest functions ${}_0{\Psi }_{q-p}$, ${}_1{\Psi }_1$, ${}_1{\Psi }_0$, respectively); (iv) the above three basic g.h.f. can be considered themselves as fractional differintegrals of the three elementary functions, depending on whether $p<q,p=q$, or $p=q+1$; and (v) the compositions of E-K operators arising in Step (iii) give generalized (q-tuple) fractional integrals or derivatives.

Thus, for the simpler case of ${}_p{F}_q$-functions, we have the following general proposition.

Theorem 7.

(Kiryakova [58]) All the generalized hypergeometric functions ${}_p{F}_q$ can be considered as generalized (q-tuple) fractional differintegrals (24), (30) (with $G_{m,m}^{m,0}$-kernels) of one of the elementary functions

$${cos}_{q-p+1}\left(z\right)(\mathrm{if}p<q),z^{\alpha }\exp z(\mathrm{if}p=q),z^{\alpha }{(1-z)}^{\beta }(\mathrm{if}p=q+1),$$

(92)

depending on whether $p<qp=qp=q+1$.

It is based on the known auxiliary result coming yet from the Bateman Project on integral transforms [7], Askey [2], Lavoie–Osler–Tremblay [124] for the R-L derivatives that we have paraphrased in terms of E-K operators (e.g., Equation (4.2.2${}^{\prime }$) in [9] and Lemma 3.2 in [58]) as follows:

$$\frac{\Gamma \left(a_p\right)}{\Gamma \left(b_q\right)}{}_p{F}_q(a_1,\dots ,a_p;b_1,\dots ,b_q;z)$$

$$=\left\{\begin{array}{c}{I}_{1,1}^{a_p-1,b_q-a_p}\left\{{}_{p-1}{F}_{q-1}(a_1,...,a_{p-1};b_1,...,b_{q-1};z)\right\}if{b}_q>a_p,\hfill \\ D_{1,1}^{b_q-1,a_p-b_q}\left\{{}_{p-1}{F}_{q-1}(a_1,...,a_{p-1};b_1,...,b_{q-1};z)\right\}if{b}_q<a_p,\hfill \end{array}\right.$$

(93)

for all complex z, and if $p=q+1$ we require additionally $\left|z\right|<1$. Then, this basic fact is to be used repeatedly, and combined with the composition/decomposition property (26) for the operators of GFC. In each of the three separate cases, we reach to one of the basic functions (92) with smallest possible first index p, namely: ${}_0{F}_{q-p}\left(z\right)={cos}_{q-p+1}\left(z\right)$; ${}_1{F}_1\left(z\right)$ and then ${}_0{F}_0\left(z\right)=expz$; and ${}_2{F}_1\left(z\right)$ and then ${}_1{F}_0(\beta ;-;z)={(1-z)}^{-\beta }$.

For the Wright generalized hypergeometric functions (4), this proposition reads almost the same, only the third basic function (for $p=q+1$) is more general, namely ${}_1{\Psi }_0=H_{1,1}^{1,1}$, and the GFC operators have as kernel the $H_{m,m}^{m,0}$-function with different parameters βs and λs in the upper and low rows.

Theorem 8.

(Kiryakova [60] (Theorem 14)) All the Wright generalized hypergeometric functions ${}_p{\Psi }_q$ can be represented as multiple (q-tuple) W-E-K fractional integrals (33), or their corresponding fractional derivatives, of one of the following three basic functions:

$${cos}_{q-p+1}\left(z\right)(ifp<q),expz(ifp=q),{}_1{\Psi }_0\left[(a,A)|z\right](ifp=q+1).$$

(94)

In this case, the basic used result is Theorem 6, following similar Steps (i)–(v) as described above.

The three cases, for both Theorems 7 and 8, are considered in detail, in separate statements.

(1) $p<q$. The Poisson integral representation (52) is extended in [9] (Ch.4) and [90] for the hyper-Bessel functions (58), $m\ge 2$, that is for the ${}_0{F}_{m-1}$-functions, via generalized fractional integrals (24) of the function ${cos}_m$, (54) as follows:

$$J_{{\nu }_1,..,{\nu }_{m-1}}^{(m-1)}\left(z\right)=\sqrt{\frac{m}{{\left(2\pi \right)}^{m-1}}}{\left(\frac{z}{m}\right)}^{{\nu }_1+\dots +{\nu }_{m-1}}{I}_{\frac{1}{m},m-1}^{(\frac{k}{m}-1),({\nu }_k-\frac{k}{m}+1)}\left\{{cos}_m\left(z\right)\right\}.$$

(95)

By analogy with the hyper-Bessel functions (58), we consider what we call the Wright hyper-Bessel functions:

$${}_0{\Psi }_m\left[\left.\begin{array}{c}-\\ (b_1,B_1),...,(b_m,B_m)\end{array}\right|-z\right]=H_{0,m+1}^{1,0}\left[z\left|\begin{array}{c}-\\ (0,1),(1-b_1,B_1),\dots ,(1-b_m,B_m)\end{array}\right.\right]$$

$$={\displaystyle \sum _{k=0}^{\infty }}\frac{z^k}{\Gamma (b_1+k{B}_1)\dots \Gamma (b_m+k{B}_m).k!}:=J_{b_1-1,\dots ,b_m-1}^{B_1,\dots ,B_m}\left(z\right).$$

(96)

The latter denotation is to remind of the analogy with the hyper-Bessel functions (58), when $\forall B_k=1$. It is easy to observe that (96) appears as special case of the multi-index Mittag–Leffler functions (39), namely: $J_{b_1-1,\dots ,b_m-1}^{B_1,\dots ,B_m}\left(z\right)=E_{(1,B_1,\dots ,B_m),(1,b_1,\dots ,b_m)}^{(m+1)}(-z)$.

We have then a result, analogous to (95), and more general than (53) for the multi-M-L functions, that: each Wright hyper-Bessel function ${}_0{\Psi }_{q-p}$, $p<q$, can be represented by means of a Poisson type integral of the ${cos}_{p-q+1}$-function, written in the form

$${}_0{\Psi }_{q-p}\left[\left.\begin{array}{c}-\\ (b_1,B_1),...,\left(b_{q-p}\right)\end{array}\right|-z\right]=J_{b_1-1,\dots ,b_m-1}^{B_1,\dots ,B_m}\left(z\right)$$

$$=I_{\left(\frac{1}{B_k}\right),\left(1\right),q-p}^{(\frac{k}{q-p+1}-1),(b_k-\frac{k}{q-p+1})}\left\{{cos}_{q-p+1}\left((q-p+1)z^{\frac{1}{q-p+1}}\right)\right\}.$$

(97)

Let us now apply to the function ${}_0{\Psi }_{q-p}$ above, p-times the results (90), (91) (Theorem 6) with $m=1$, combined with the composition rule for the W-E-K integrals (33). Then, we obtain the following:

Theorem 9.

(Kiryakova [60] (Theorem 15)) Each ${}_p{\Psi }_q$-function with $p<q$ is a generalized q-tuple W-E-K fractional (differ-)integral operator of ${cos}_{q-p+1}\left(z\right)$,

$${}_p{\Psi }_q\left[\left.\begin{array}{c}(a_1,A_1),\dots ,(a_p,A_p)\\ (b_1,B_1),\dots ,(b_q,B_q)\end{array}\right|-z\right]=I_{\left(\frac{1}{B_k}\right),\left({\lambda }_k\right),q}^{\left({\gamma }_k\right),\left({\delta }_k\right)}\left\{{cos}_{q-p+1}\left((q-p+1)z^{\frac{1}{q-p+1}}\right)\right\},$$

(98)

with the following parameters:

$${\gamma }_k=\left\{\begin{array}{c}\frac{k}{q-p+1}-1,\hfill \\ a_{k-q+p}-1,\hfill \end{array};\right.{\delta }_k=\left\{\begin{array}{c}{b}_k-\frac{k}{q-p+1},\hfill \\ b_k-a_{k-q+p},\hfill \end{array};\right.{\lambda }_k=\left\{\begin{array}{cc}1,\hfill & k=1,\dots ,q-p\hfill \\ \frac{1}{A_{k-q+p}},\hfill & k=q-p+1,\dots ,q.\hfill \end{array}\right.$$

If the following conditions are satisfied:

$$b_k>\frac{k}{q-p+1},k=1,\dots ,q-p;b_k>a_{k-q+p}>0,k=q-p+1,\dots ,q,$$

$$B_k\ge 1,k=1,\dots ,q-p;B_k\ge A_{k-q+p},k=q-p+1,\dots ,q,$$

then relation (98) gives a Poisson type integral representation; otherwise, the operator in the R.H.S. should be interpreted as a multiple W-E-K derivative (see, e.g., Definition 7 in [60]), and then (98) turns into a new Rodrigues type differential formula, or a mixed differ-integral representation.

The particular case of Poisson type representation (53) for the multi-index M-L function has been already stated as Theorem 3 in Section 4.

In the other two cases, $p=q$ and $p=q+1$, the starting results for ${}_p{\Psi }_q$ were formulated as Lemmas 11 and 12 in Kiryakova [60]:

$${}_1{\Psi }_1\left[\left.\begin{array}{c}(a_1,A_1)\\ (b_1,B_1)\end{array}\right|z\right]=W_{1/B_1,1/A_1}^{a_1-1,b_1-a_1}\left\{expz\right\},if{A}_1\ge B_1,b_1\ge a_1,for\left|z\right|<\infty ;$$

(99)

$${}_2{\Psi }_1\left[\left.\begin{array}{c}(a_1,A_1),(a_2,A_2)\\ (b_1,B_1)\end{array}\right|z\right]=W_{1/B_1,1/A_1}^{a_1-1,b_1-a_1}\left\{{}_1{\Psi }_0\left[\left.\begin{array}{c}(a_2,A_2)\\ -\end{array}\right|z\right]\right\}$$

(100)

$$=W_{1/B_1,1/A_1}^{a_1-1,b_1-a_1}\left\{H_{1,1}^{1,1}\left[-z\left|\begin{array}{c}(1-a_2,A_2)\\ (0,1)\end{array}\right.\right]\right\},$$

$$\mathrm{if}{A}_1\ge B_1,b_1\ge a_1;\mathrm{and} \mathrm{if}{A}_2<1,\mathrm{for}\left|z\right|<\infty ;\mathrm{or} \mathrm{if}{A}_2=1,\mathrm{for}\left|z\right|<1.$$

After additional $(p-1)$ steps, from ${}_p{\Psi }_q$ passing via ${}_1{\Psi }_1$ to ${}_0{\Psi }_0$, respectively, to ${}_1{\Psi }_0$, the following explicit results for the statement in Theorem 8 are provided in [60].

(2) $p=q$.

Theorem 10.

If $p=q$, each g.h.f. ${}_p{\Psi }_p\left(z\right)$ is an p-tuple W-E-K fractional integral of the exponential function, namely, if $B_k\ge A_k>0,b_k>a_k>0$, $k=1,\dots ,p$:

$${}_p{\Psi }_q\left[\left.\begin{array}{c}(a_1,A_1),...,(a_p,A_p)\\ (b_1,B_1),...,(b_p,B_p)\end{array}\right|z\right]=I_{\left(\frac{1}{B_k}\right),\left(\frac{1}{A_k}\right),p}^{(a_k-1),(b_k-a_k)}\{expz\},for\left|z\right|<\infty .$$

(101)

If for some indices k, the above inequalities for parameters are not satisfied, representation (101) turns into differ-integral one, or in special cases to purely differential one.

Theorem 10 suggests us to separate the g.h.f-s ${}_p{\Psi }_p$ with $p=q$ in a class of so-called Wright g.h.f. of confluent type, involving the confluent hypergeometric function ${}_1{F}_1(a;b;z)=\Phi (a;b;z)$ and $expz$ as the simplest cases.

(3) $p=q+1$. Analogously, we call the ${}_{q+1}{\Psi }_q$-functions with $p=q+1$ as Wright g.h.f. of Gauss type, since the simplest case of such special function is the Gauss function. We have following specific result.

Theorem 11.

Each Wright g.h.f. of Gauss type ${}_p{\Psi }_q$, that is with $p=q+1$, is a q-tuple Wright–Erdélyi–Kober fractional integral (or differ-integral) of the ${}_1{\Psi }_0$-function. Namely, for $0<A_0\le 1$ and $b_k>a_k>0,k=1,\dots ,p$:

$${}_{q+1}{\Psi }_q\left[\left.\begin{array}{c}(a_0,A_0),(a_1,A_1),\dots ,(a_q,A_q)\\ (b_1,B_1),\dots ,(b_q,B_q)\end{array}\right|z\right]=I_{\left(\frac{1}{B_k}\right),\left(\frac{1}{A_k}\right),q}^{(a_k-1),(b_k-a_k)}\left\{{}_1{\Psi }_0\left[\left.\begin{array}{c}(a_0,A_0)\\ -\end{array}\right|z\right]\right\}$$

(102)

$$=I_{\left(\frac{1}{B_k}\right),\left(\frac{1}{A_k}\right),q}\left\{H_{1,1}^{1,1}\left[-z\left|\begin{array}{c}(1-a_0,A_0)\\ (0,1)\end{array}\right.\right]\right\},if{A}_0<1for\left|z\right|<\infty ;or if{A}_0=1,for\left|z\right|<1.$$

For other arrangements between $b_k$ and $a_k$, the operator in (102) is a generalized fractional derivative.

For particular choices of parameters $b_k$ and $a_k$ not satisfying the conditions $b_k>a_k>0$, some integer order differentiations appear in place of the fractional integrals or derivatives and lead to Rodrigues type differential formulas, analogous to these for the classical orthogonal polynomials.

Note that the integral representation (102) generalizes the Euler integral formula for the Gauss hypergeometric functions that serves for its analytical extension outside $\left|z\right|<1$ to the domain $\left|\arg \right(1-z\left)\right|<\pi$:

$${}_2{F}_1(a_1,a_2;b_1;z)=\frac{\Gamma \left(b_1\right)}{\Gamma \left(a_2\right)\Gamma (b_1-a_2)}{\displaystyle \underset{0}{\overset{1}{\int }}}\frac{{(1-\sigma )}^{b_1-a_2-1}{\sigma }^{a_2-1}}{{(1-z\sigma )}^{a_1}}d\sigma ,b_1>a_2>0.$$

(103)

This gave us the reason to name ${}_p{\Psi }_q$ with $p=q+1$ as a Gauss type g.h.f.

In particular, for $A_0=1$, the basic function in the case $p=q+1$ reduces to the geometric series:

$${}_1{\Psi }_0\left[\left.\begin{array}{c}(a_0,1)\end{array}\right|z\right]=H_{1,1}^{1,1}\left[-z\left|\begin{array}{c}(1-a_0,1)\\ (0,1)\end{array}\right.\right]=G_{1,1}^{1,1}\left[-z\left|\begin{array}{c}1-a_0\\ 0\end{array}\right.\right]={}_1{F}_0(a_0;-;z)={(1-z)}^{-a_0}.$$

Therefore, based on the statements in Theorems 7–11, we suggest a classification of the classical SF and of the SF of FC into three classes, namely “Bessel”, “confluent” and “Gauss” types, depending on whether $p<q$, $p=q$ or $p=q+1$. This approach can facilitate applied scientists and engineers, escaping a deep knowledge on SF, to think of them in a very general view as similar to a cosine- (Bessel) function, exponent or geometric series, because the fractional integrations keep in some sense the asymptotic and general behavior.

The results from Theorems 7–11 for ${}_p{\Psi }_q$, and their specifications for the ${}_p{F}_q$-functions, yield also several new integral and differential formulas for them, with possible hints for computational procedures.

Below, we mention some few of them, say in the simpler cases of ${}_p{F}_q$-functions.

The case $p=q$: For the g.h.f. ${}_p{F}_p$, the integral representation can be written not only by means of $G_{p,p}^{p,0}$-functions in the kernel, but also avoiding SF due to decomposition property (26). Thus, we have an integral formula, as follows:

$${}_p{F}_p(a_1,\dots ,a_p;b_1,\dots ,b_p;z)=B{z}^{1-a_1}{I}_{1,p}^{(a_k-a_1),(b_k-a_k)}\left\{z^{a_1-1}expz\right\}$$

$$=B{\displaystyle \underset{0}{\overset{1}{\int }}}\underset{\left(p\right)}{\dots }{\displaystyle \underset{0}{\overset{1}{\int }}}{\displaystyle \prod _{k=1}^p}\left[\frac{{(1-{\sigma }_k)}^{b_k-a_k-1}{\sigma }_k^{a_k-1}}{\Gamma (b_k-a_k)}\right]exp(z{\sigma }_1\dots {\sigma }_p)d{\sigma }_1\dots d{\sigma }_p,B:={\displaystyle \prod _{j=1}^p}\frac{\Gamma \left(b_j\right)}{\Gamma \left(a_j\right)},$$

(104)

under conditions $b_k>a_k>0,k=1,\dots ,p$. If the parameters do not satisfy them, the GFC operator above is interpreted as generalized fractional derivative of the form (30).

Specially, let all the differences $a_k-b_k={\eta }_k,k=1,\dots ,p$ be non-negative integers. In this case, we call the ${}_p{F}_p$-functions as “spherical” g.h.f. of confluent type, using the analogy with the spherical Bessel, hyper-Bessel functions and spherical multi-M-L functions $E_{\left({\alpha }_i\right),\left({\beta }_i\right)}\left(z\right)$ with “semi-integer” indices, mentioned in Remark 1. Then, the operator in (104) turns into a differential operator $D_{\eta }$ of integer order $\eta ={\eta }_1+\dots +{\eta }_p\ge 0$ of the form (27), and we obtain a differential formula of the form

$$\frac{\Gamma \left(a_p\right)}{\Gamma \left(b_q\right)}{}_p{F}_q(a_1,\dots ,a_p;b_1,\dots ,b_q;z){}_p{F}_p(b_1+{\eta }_1,\dots ,b_p+{\eta }_p;b_1,\dots ,b_p;z)$$

$$={}_p{F}_p(b_1+{\eta }_1,\dots ,b_p+{\eta }_p;b_1,\dots ,b_p;z)$$

$$=\left[{\displaystyle \prod _{j=1}^p}\frac{\Gamma \left(b_j\right)}{\Gamma (b_j+{\eta }_j)}\right]\left[{\displaystyle \prod _{k=1}^p}{\displaystyle \prod _{j=1}^{{\eta }_k}}(z\frac{d}{dz}+b_k+j-1)\right]\left\{expz\right\}=Q_p\left(z\right)\left\{expz\right\}.$$

(105)

The representation (105) gives an example how differential formulas for the “spherical” g.h.f. introduced by Kiryakova [9] can be used for their explicit calculation, especially in the case $p=q$ in the form $Q_p\left(z\right)\left\{expz\right\}$, where $Q_p\left(z\right)$ is a p-degree polynomial. A special case of (105) with $b_k={\eta }_k=1,k=1,\dots ,p$ and $Q_p\left(z\right)={\displaystyle \frac{d}{dz}}{\left(z{\displaystyle \frac{d}{dz}}\right)}^p$ was presented by Prudnikov–Brychkov–Marichev [14] (p. 593).

The case $p=q+1$: For the Gauss type g.h.f. ${}_{q+1}{F}_q$, we have in the unit disk $\left|z\right|<1$ an integral representation (if written by repeated integral with no use of the kernel $G_{q,q}^{q,0}$-function), for $b_k>a_{k+1},k=1,..,q$:

$${}_{q+1}{F}_q(a_1,\dots ,a_{q+1};b_1,\dots ,b_q;z)$$

$$=\left[{\displaystyle \prod _{j=1}^q}\frac{\Gamma \left(b_j\right)}{\Gamma \left(a_{j+1}\right)\Gamma (b_j-a_{j+1})}\right]z^{1-a_2}{I}_{1,q}^{{(a_{k+1}-1)}_1^q,{(b_k-a_{k+1})}_1^q}\left\{z^{a_2-1}{(1-z)}^{-a_1}\right\}$$

$$=\left[{\displaystyle \prod _{j=1}^q}\frac{\Gamma \left(b_j\right)}{\Gamma \left(a_{j+1}\right)\Gamma (b_j-a_{j+1})}\right]{\displaystyle \underset{0}{\overset{1}{\int }}}\underset{\left(p\right)}{\dots }{\displaystyle \underset{0}{\overset{1}{\int }}}{\displaystyle \prod _{j=1}^q}\left[{(1-{\sigma }_k)}^{b_k-a_{k+1}-1}{\sigma }_k^{a_{k+1}-1}\right]{(1-z{\sigma }_1\dots {\sigma }_q)}^{-a_1}d{\sigma }_1\dots d{\sigma }_q.$$

(106)

In this form, (106) can also be found in [14] (p. 438). In the case $q=1$, this is exactly the Euler integral formula (103) for the Gauss hypergeometric function. Similarly, (106) proposes a way for an analytical continuation of the functions ${}_{q+1}{F}_q\left(z\right)$ outside the unit disk to the domain $|arg(1-z\left)\right|<\pi$.

In the case when the $a_k$’s and $b_k$’s do not satisfy the above conditions, the operator in (106) turns into a generalized fractional derivative, and this also provides useful corollaries. By analogy with the previous two cases ($p<q$ and $p=q$), we introduce the notion of spherical g.h.f. of Gauss type when all the differences $a_{k+1}-b_k={\eta }_k,k=1,\dots ,q$ are non-negative integers. Then, ${}_{q+1}{F}_q\left(z\right)$ is representable by a purely differential operator of a function ${(1-z)}^{-a_1}$, and a special case of such differential formula is presented in [14] (p. 572).

Another interesting case concerns the so-called hypergeometric polynomials

$${}_{q+1}{F}_q(-n,a_1,\dots ,a_p;b_1,\dots ,b_q;z)={\displaystyle \sum _{k=0}^n}\frac{{(-n)}_k{\left(a_1\right)}_k\dots {\left(a_p\right)}_k}{{\left(b_1\right)}_k\dots {\left(b_q\right)}_k}\frac{z^k}{k!}.$$

(107)

Taking $a_{q+1}=-n$ with integer $n\ge 0$ and $a_k>b_k>0,k=1,\dots ,q$, the fractional derivative form of the operator in (106) provides the Rodrigues type formula ([9] (Ch.4)):

$$\begin{array}{l}\left[{\displaystyle \prod _{j=1}^q}\frac{\Gamma \left(a_j\right)}{\Gamma \left(b_j\right)}\right]{}_{p+1}{F}_q(-n,a_1,\dots ,a_q;b_1,\dots ,b_q;z)=D_{1,q}^{(b_k-1),(a_k-b_k)}\left\{{(1-z)}^n\right\}\\ =z^{1-a_q}{D}_{1,q}^{(b_k-a_q),(a_k-b_k)}\left\{z^{a_1-1}{(1-z)}^n\right\}=z^{1-b_q}{D}^{a_q-b_q}{z}^{a_p-b_{q-1}}{D}^{a_{p-1}-b_{q-1}}\\ \times \dots z^{a_3-b_2}{D}^{a_2-b_2}{z}^{a_2-b_1}{D}^{a_1-b_1}\left\{z^{a_1-1}{(1-z)}^n\right\}.\end{array}$$

(108)

Special cases of (108) yield some classical Rodrigues formulas. For example, $p=q=1$ with $a_1=n+1,b_1=1$ and $z\to \frac{1-z}{2}$ gives the Rodrigues formula for the Legendre polynomials:

$$\begin{array}{ll}{P}_n\left(z\right)={(-1)}^n{}_2{F}_1(-n,n+1;1;\frac{1-z}{2})& =\frac{{(-1)}^n}{n!}\frac{d^n}{d{z}^n}\left[{\frac{1-z}{2}}^n·{\frac{1+z}{2}}^n\right]\\ & =\frac{1}{2^{n}n!}\frac{d^n}{d{z}^n}\left\{{(z^2-1)}^n\right\},\end{array}$$

and $p=q=2$ with $a_1=n+1,b_1=1,a_2=\zeta ,b_2=p(\zeta >p>0)$ gives the Rodrigues formula for the Rice polynomials, namely

$$R_n\left(z\right)={}_3{F}_2(-n,n+1,\zeta ;1,p;z)=\frac{\Gamma \left(p\right)}{n!\Gamma \left(\zeta \right)}\left[\frac{d^n}{d{z}^n}{z}^{1-p}{\left(\frac{d}{dz}\right)}^{\zeta -p}\right]\left\{z^n{(1-z)}^n\right\}.$$

9. Numerical Aspects of SF of FC

In the days before the electronic computers, the necessary complement to a special function was the computation, by hand, of extended tables of its values, intended to make the function available for users, similarly to the familiar logarithm tables. After mechanical calculators appeared and were more widespread, several huge special-function-table projects were started. Let us mention as examples the handbooks Gradshteyn–Ryzhik [126] and Magnus–Oberhettinger [127], both initiated in 1943.

R. Askey (at the Conf. “SF 2000”, ASU): “...The advent of fast computing machines was thought to have made special functions a subject of the past. The reality has been different. Continued development of older functions and the introduction of new special functions has been the reality … and still remains to be discovered … The classical handbooks as mentioned, although useful as references, maybe no longer the primary means of accessing the special functions of mathematical physics. A number of high level programs appeared that are better suited for this challenging purpose, to mention as Mathematica, Maple, Matlab, Mathcad, …”

We like to add a citation from Stephen Wolfram [128] (Wolfram Mathematica), “… and special functions became a big business. Table making had become a major activity for the governments, and was thought strategically important. Particularly for things like navigation, nuclear physics, military reasons, H-bomb, etc. And there were lots of tables …The aspects of the theory then mattered might be as two: – for numerical analysis, discovery of infinite series or other analytical expression allowing rapid calculation; and – reduction of as many functions as possible to the given (better known) function …” (Author’s comment: compare with the approach applied in works of Kiryakova as [9] (Ch.4), [58,60], discussed in Section 8). (S. Wolfram:): “…There gradually started to appear systematic reference works on the properties of special functions. Each one based on lot of work …”, “…I guess integrals are timeless. They don’t really bear the marks of the human creators. So we have the tables, but we really don’t quite know where they came from …”. (Author’s comment: However, it seems Marichev knew, and we refer to his book [11]).

In the 1960s and 1970s, a lot of efforts started for developing numerical algorithms for computers. Evaluation of special functions became a favorite area. S. Wolfram: “Well, a few years passed. And in 1986, I started designing Mathematica. I wanted to be sure to do a definitive job, and to have good numerics for all functions, for all values of parameters, anywhere in the complex plane, to any precision …And I remember very distinctly a phone call I had with someone at a government lab. And there was a silence. And then he said: “Look, you have to understand that by the end of the 1990s we hope to have the integer-order Bessel functions done to quad precision.” … (S.W., cont’d): “You know, it’s actually quite a difficult thing to put a special function into Mathematica. You don’t just have to do the numerics… So what makes a special function good? Well, we can start thinking about that question empirically. Asking what’s true about the special functions we normally use. And of course, from what we have in Mathematica and in our Wolfram Functions Site [129], we should be in the best position to answer this.”

Let us note that the standard SF—the hypergeometric functions (Gauss, ${}_p{F}_q$-), the Meijer G-function, etc.—are well presented there ([129]), but (it seems) none of the M-L type, Wright and H-functions, that is cases of SF of FC, are available yet. Meanwhile, the fractional nature of the world needs better reflection by fractional order (FO) models in whose solutions the so-called SF of FC appear. Thus, it is yet a challenging trend to be developed.

Here, we try to provide only a short information on some “recent” numerical jobs done with respect to the M-L-function, classical Wright function and only few of their extensions.

For numerical algorithms and results in the case of the Mittag–Leffler functions (one- and two-parameter and matrix analog), we start with reference to Caputo–Mainardi [130] (1971). We note that this is one of the first works to propose a plot of the M-L function. On those times, without possibilities to take advantage of software packages such as Mathematica, Maple and Matlab, this task was difficult, as it was managing series expansions convergent only in the mathematical sense but not in the numerical sense. Further, some other authors worked on similar problems either simultaneously but independently, or in years afterward: Gorenflo–Loutchko–Luchko [131], 2002; Diethelm–Ford–Freed–Luchko [132], 2005; Podlubny [133], 2005–2009–2012, (v 1.2.0.0) 2021; Hilfer–Seybold and Seybold–Hilfer [134,135], 2006–2008; Garrappa [136], 2015 and Garrappa–Popolizio [137], 2018; etc.

Numerical algorithms and results on the (classical) Wright function (56) and its special cases, including the Mainardi function, can be found in works by Luchko [138], 2008; Luchko–Trujillo–Velasco [139], 2010; Consiglio [140], 2019; Mainardi–Consiglio [97], 2020; etc.

Concerning the Prabhakar (three-parameter) M-L type function (38), see Garrappa [136], 2015; etc.

For the generalized exponential integrals as (63) and related generalized trigonometric functions involving M-L functions (in the sense of 6.2.), and shown to be Wright g.h.f. ${}_2{\Psi }_2$, one can find some tables and plots for physically interesting parameters and related models, proposed by Mainardi and Masina [102] (2018) and Paris [103] (2020).

This list can surely be extended with more information.

We would like to attract readers’ attention to the challenging Open Problem: What about possibilities for numerical and graphical interpretations, plots and tables and implementing software packages for some more general Special Functions of Fractional Calculus, such as the multi-index Mittag-Leffler functions? At least, to treat illustrative examples for few typical sets if multi-indices?

10. Conclusions

In this survey, under the notion of Special Functions of Fractional Calculus (SF of FC), we have in mind the Fox H-function and the Wright generalized hypergeometric functions ${}_p{\Psi }_q$, including the Mittag–Leffler function, its multi-index extensions and all their particular cases. The standard (classical) special functions (SF) naturally come as part of this scheme, as cases of the Meijer G-function and of the ${}_p{F}_q$-functions, including so many named SF and orthogonal polynomials. Here, we try to review some of the basic results on the theory of the SF of FC, obtained in author’s works over more than 30 years, and support the wide spreading and important role of these functions by several examples.

The short outline of the contents is as follows:

In Section 1, we start with a historical introduction to pay tribute to the older projects that gave life to the contemporary development of the topic. Some short definitions and facts on the considered basic special functions are given in Section 2. In Section 3, we pay attention to the use of the H- and G-functions, especially of orders $(m,0;0,m)$ and $(m,0;m,m)$, as kernel-functions of generalized integral transforms of Laplace type and of the operators of the so-called generalized fractional calculus (GFC). In Section 4, we introduce the Mittag–Leffler functions and the multi-index Mittag–Leffler functions, with short information on their properties derived in author’s works. Section 5 and Section 6 contain long lists of examples of SF that appear as cases of the multi-index Mittag–Leffler functions and in more general setting, of the Wright generalized hypergeometric functions ${}_p{\Psi }_q$. These include also citations to many other authors who introduced and applied such functions in their works. The author’s unified approach to evaluate images of classical SF and of SF of FC under operators of FC and GFC is shortly described in Section 7, because the details are presented in another survey paper in the same journal [105]. In Section 8, we collect some of our basic propositions on the representations of the SF and of SF of FC as operators of GFC of three basic and simplest elementary functions and propose a classification of the SF based on the cases $p<q$, $p=q$ and $p=q+1$. Thus, a new sight on the theory of SF is proposed. Since the computational aspects related to the considered SF are of important interest for their applications, in Section 9, we provide some short information on the state of affairs and some recent works on this direction by other authors. A provoking challenge in this respect is mentioned.