On a Five-Parameter Mittag-Lefﬂer Function and the Corresponding Bivariate Fractional Operators

: Several extensions of the classical Mittag-Lefﬂer function, including multi-parameter and multivariate versions, have been used to deﬁne fractional integral and derivative operators. In this paper, we consider a function of one variable with ﬁve parameters, a special case of the Fox– Wright function. It turns out that the most natural way to deﬁne a fractional integral based on this function requires considering it as a function of two variables. This gives rise to a model of bivariate fractional calculus, which is useful in understanding fractional differential equations involving mixed partial derivatives.


Introduction
The original Mittag-Leffler function E α (z), applied to one variable according to one parameter, was first defined and studied by Gösta Mittag-Leffler in the 1900s. In the hundred years since then, many variants and extensions of this function have been defined, including functions of more than one variable and functions with arbitrarily many parameters [1,2].
One of the major motivations for studying Mittag-Leffler functions is their relationship with fractional calculus [2][3][4][5], a field of mathematics which has become very popular due to its many applications in various areas of science [6][7][8]. Mittag-Leffler functions emerge naturally as the solution to some elementary fractional differential equations, and their eigenfunction properties have led them sometimes to be called "fractional exponential functions" [9].
Among the many extensions of the Mittag-Leffler function to more variables and parameters, we mention just a few which are of particular importance or interest in motivating the present work.
• A Mittag-Leffler function of one variable with three parameters was defined by Prabhakar [10] to solve a certain singular integral equation. Its use as an integral kernel gave rise to a model of fractional calculus which has a semigroup property and which is already broad enough to include many other named fractional-calculus operators [11][12][13], although it is itself a special case of the general Fox-Wright function. • A Mittag-Leffler function of n variables with n + 1 parameters was used by Luchko et al. [14,15] to solve multi-term linear fractional differential equations involving multiple independent fractional orders. Note that this is independent from the Mittag-Leffler function of n variables with 2n + 1 parameters which was defined by Saxena et al. [16] and which gives rise to a model of fractional calculus with a semigroup property [17]: neither of these general functions is a special case of the other.
• A Mittag-Leffler function of two variables with four parameters and a Mittag-Leffler function of three variables with five parameters were recently defined [18,19] and used to solve multi-order systems of fractional differential equations [20,21]. These are closely related to the general function considered by Luchko et al., but they are not special cases of it. Used as kernels, they also give rise to new models of fractional calculus with a semigroup property [18,21]. For other types of bivariate Mittag-Leffler function that have been defined in the literature, we refer to the papers [22][23][24].
All of these Mittag-Leffler type functions have been connected in some way to fractional calculus and fractional differential equations. Indeed, the strong connection between special functions and fractional calculus has been known for decades and continues to be written about today [25][26][27], while multi-parameter and multi-variable generalisations are being studied both for special functions and for the operators of fractional calculus [17].
Fractional partial differential equations have been an interesting and challenging topic of study [3], with various methods able to be extended (with modifications) from classical partial differential equations in order to solve them, such as the unified transform method [28], distribution theory [29,30], and weak solutions [31,32]. Much attention has been paid to partial differential equations involving the fractional Laplacian operator [33,34], but less attention has gone to differential equations involving mixed partial fractional derivatives. They are mentioned in the classical textbooks [3] ( §6.1.1) and [6] ( §24.2), and in a few papers (e.g., [35][36][37]), but in general they have attracted little notice in the research literature on fractional partial differential equations.
In the work below, we study a specific type of Mittag-Leffler function, initially a function of one variable with five parameters, and define a double fractional integral operator by converting this function to a bivariate version. In this way, it is possible to involve double integrals and derivatives, with respect to two variables, while preserving the relatively simple structure of a single power series defining the Mittag-Leffler function. These operators are therefore useful in the study of fractional partial differential equations, including those involving mixed partial fractional derivatives.
As an initial motivation, let us consider the following bivariate Abel equation of the second kind: One can reformulate this equation in terms of the standard Riemann-Liouville fractional integral operators as: 1 + λI α 1 t I α 2 s u(t, s) = f (t, s).
Formally, by means of symbolic operational calculus without regard for rigour, this can be solved as follows: This last formal result can be written in terms of the following Mittag-Leffler type function: namely, in the following way, making use of the double (two-dimensional) Laplace convolution operator denoted here as " * ": Therefore, it makes sense to define an integral operator E α 1 ,α 2 acting on functions of two variables as follows: This operator emerges naturally from consideration of bivariate Abel equations, but it has the drawback of lacking a semigroup property, or index law, in any of the parameters: if we take a composition of this operator with itself, we do not find another operator in the same form for different values of the parameters.
It should be noted at this point that some useful and important operators, which have been criticised for lacking a semigroup property, can be embedded into a larger class of fractional-calculus operators which has an index law and which therefore contains both the original operators and their compositions. For example, Prabhakar fractional calculus forms a class of operators which contains various useful operators that lack semigroup properties in themselves, their compositions being different elements of the Prabhakar class [11]. We can apply the same way of thinking here, extending the basic two-parameter Mittag-Leffler function (2) and the associated integral operator to more general versions where a semigroup property can be found.
In this case, it is sufficient to add three extra parameters in the definition in order to obtain a generalised operator which has a semigroup property. We consider, throughout this paper, the following five-parameter Mittag-Leffler function of one variable: where α 1 , α 2 , β 1 , β 2 , γ are complex constants (with some constraints to be determined later) and (γ) k is the Pochhammer symbol defined by This paper is devoted to a detailed study of the five-parameter Mittag-Leffler function (4), the associated bivariate fractional integral operators, related concepts such as the corresponding fractional derivative operators, and special cases of particular interest.
Specifically, the structure of the paper is as follows. We discuss the five-parameter Mittag-Leffler function and its properties in Section 2, then pass to bivariate fractional calculus in Section 3, firstly integral operators and then derivative operators. In Section 4, we investigate the case where the operators are non-singular and expressible as finite sums of Riemann-Liouville integrals. Section 5 is devoted to discussion of applications and potential future work.

The Five-Parameter Mittag-Leffler Function
Definition 1. Let α 1 , α 2 , β 1 , β 2 , γ ∈ C be five parameters satisfying Re(α 1 + α 2 ) > 0. The five-parameter Mittag-Leffler function applied to a single variable z is defined by the following power series: This can be seen as a special case of the Fox-Wright function, which is defined [3] ( §1.11) by where, in [3], the parameters are given as a i , b j ∈ C, α i , β j ∈ R satisfying α i = β j for i = 1, · · · , p; j = 1, · · · , q, and it is proved that this power series is absolutely convergent for all z ∈ C in the case ∆ : The five-parameter Mittag-Leffler function can be written in terms of the Fox-Wright function as follows: As an immediate consequence of the convergence result shown in [3] (Theorem 1.5) and [38], we conclude that, if α 1 , α 2 ∈ R and α 1 + α 2 > 0, then the power series defining E γ α 1 ,α 2 ;β 1 ,β 2 (z) is locally uniformly convergent and therefore it is an entire function.
In fact, there is no need to assume that any of the parameters are real. This assumption is imposed in [3,38] to simplify the calculations, but the same convergence result can be proved for complex values of the parameters by using Stirling's formula and the ratio test, similarly to the second author's work in [39]. The condition then to be imposed on the parameters is Re(α 1 + α 2 ) > 0, as stated in Definition 1.
It is worth noting that a four-parameter special case of the general Fox-Wright function was recently given particular consideration by Luchko [40], namely the following function: which is the case γ = 1 of our five-parameter Mittag-Leffler function. It is also a special case of the vector-index Mittag-Leffler function studied by Al-Bassam and Luchko [41], which is defined as follows: for any n ∈ N. The case n = 1 gives the usual two-parameter Mittag-Leffler function, while the case n = 2 gives the four-parameter Wright function (5) considered by Luchko. The function (6) was further extended [42] to include a numerator (γ) κn in terms of two further parameters.
Our function is also a special case of the general Fox-Wright function, but it is different from the special cases (5) and (6) considered previously, because of the Pochhammer symbol appearing in the numerator. This extra parameter and Pochhammer symbol is important because, as we show below, it gives rise to a semigroup property for the resulting fractional integral operators -a property which is lacking for functions such as (5) and (6) which were considered before, and even for functions such as that in [42] which contain a generalised Pochhammer symbol. The semigroup property will arise from the idea of turning this univariate function into a bivariate integral operator, a notion which we justify below from consideration of the gamma functions on the denominator.
Using the Mellin-Barnes integral representation of the Fox-Wright function p Ψ q (z), proved in [3], we have for α 1 > 0, α 2 > 0, and γ, β 1 , where the integration path L is a Bromwich contour starting at a point C − i∞ and terminating at a point C + i∞ and separating the poles −m (m = 0, 1, 2, · · · ) of the gamma function Γ(s) to the left and the poles γ + l (l = 0, 1, 2, · · · ) to the right with the assumption −m = γ + l for l, m = 0, 1, 2, · · · . Another complex integral representation is given by the following Theorem.
be the contour which is the union of the following three parts oriented according to non-decreasing arg τ: 1.
the ray arg τ = φ, |τ| ≥ ε, Then, the five-parameter Mittag-Leffler function possesses the following complex integral representations: Proof. Using the known representation [43] 1 with z replaced by α 2 k + β 2 , and using the fact that the series represents an entire function, which allows the interchange of the series and the integral, gives which directly yields (7), the first of the desired formulae, in terms of the three-parameter Mittag-Leffler function of Prabhakar.
On the other hand, using (9) for both of the terms 1 Γ(α 1 k+β 1 ) and 1 Γ(α 2 k+β 2 ) , we have which directly yields (8), the second of the desired formulae. This derivation assumes that |zξ −α 1 τ −α 2 | < 1, which is uniformly true on the given contours provided that Re(α 1 ) > 0, Re(α 2 ) > 0, and |z| is sufficiently small. The bound on |z| can be removed, by analytic continuation in z, to yield the same result valid for all z.
Thus far, we have considered the function (4) as a five-parameter function of a single variable z. As a function of z, this function has several properties, such as complex integral representations, which we prove above. However, when we start to construct connections with the topic of fractional calculus, it is more natural to consider a related bivariate function instead.
To see why, consider that fractional differintegral relations for special functions, such as many found in [44,45], usually make use of gamma functions appearing in power series, in the following way: with appropriate choices of ν in order to use such identities for every term of a power series expansion. For many univariate functions defined by power series, the above relations give rise to interesting identities between special functions. In our case, however, the power series (4) has two different gamma functions in the denominator, each of them involving k times a different parameter, α 1 or α 2 . Therefore, to take full advantage of the function's symmetry in α 1 and α 2 , we should use two different variables: one raised to the power of α 1 , the other raised to the power of α 2 , and then both of them further raised to the power of k. Motivated by this discussion, we now begin to study, instead of the univariate function E γ α 1 ,α 2 ;β 1 ,β 2 (z), the closely related bivariate function E γ α 1 ,α 2 ;β 1 ,β 2 (x α 1 y α 2 ). A similar idea, substituting products of fractional powers of new variables instead of old variables, was used in a 2017 paper of the first author [23], but in that case it was used to convert a bivariate function to another bivariate function of different variables. Here, we use it to convert a univariate function to a bivariate function. This seems at first sight as an unnecessary complication, but we see that it makes many things more natural and smooth in the studies related to the five-parameter Mittag-Leffler function.

Definition 2 ([3]
). The double fractional integrals of a bivariate function f (x, t) are defined in the natural way by combining fractional integrals with respect to x and t, namely as follows for Similarly, the partial fractional derivatives of a bivariate function f (x, t) are defined as follows, for λ, µ ∈ C with Re(λ) > 0, Re(µ) > 0 and for n = Re(µ) Proof. We know from above that the infinite series defining the five-parameter Mittag-Leffler function is locally uniformly convergent, so we have the right to interchange the order of this series with fractional integral operators. Since all the exponents of (x − a) and (y − c) are greater than −1, we have: Note how the two gamma functions in the denominator of the series mesh together precisely with the gamma-function quotients arising from the two fractional integrals, in order to achieve the desired result.
Proof. This can be deduced from Lemma 1 by analytic continuation in the variables λ and µ, using the analyticity properties of fractional differintegrals. Alternatively, it can be proved from the series expansion of the function, following similar lines as the proof of Lemma 1.

Lemma 3.
Setting all parameters to 1 in the bivariate form of the five-parameter Mittag-Leffler function, we can recover a case of the modified Bessel function with parameter 0: Proof. This follows immediately from the series definition of the function.

Bivariate Operators with Five-Parameter Mittag-Leffler Kernels
In this section, we move on from functions to operators. Having established the five-parameter Mittag-Leffler function and some of its properties, we now wish to define a fractional integral operator using this function as a kernel, following in the footsteps of other papers [10,18,46] which defined new models of fractional calculus by using various types of Mittag-Leffler functions as kernels.
The five-parameter Mittag-Leffler function (4) is defined by a single power series in terms of a single variable z. When we use it to define a fractional integral operator, however, we must transform it to a bivariate function, again using a single summation, but this time in terms of two independent variables x, y. Correspondingly, we use a double integral and create an operator to be applied to bivariate functions. This is the only natural way to define a model of fractional calculus using the five-parameter Mittag-Leffler function, because the two gamma functions in its denominator will give rise to a double fractional integral in each summand when we wish to write a series representation for the new operator: we need two separate powers in the integrand, one corresponding to each of the gamma functions from the denominator. Definition 3. The fractional integral operator based on the five-parameter Mittag-Leffler function as a kernel is given by Note that these restrictions on the parameters are necessary in order to have a properly convergent singular integral for all reasonably well-behaved functions f (more details later on the function space for f ). It is clear that the two-parameter bivariate integral operator defined in (3) above is a special case of this new fractional integral operator: .
Since the five-parameter Mittag-Leffler function is an entire function when Re( Since a, b, c, d, β 1 , β 2 are fixed, the result is proved.
. Then, the operator a,b I γ;λ α 1 ,α 2 ;β 1 ,β 2 can be represented as an infinite series of double fractional integrals of Riemann-Liouville type: where the right-hand side is locally uniformly convergent for Proof. Since E γ α 1 ,α 2 ;β 1 ,β 2 (u) is an entire function defined by a locally uniformly convergent power series, and , we have the right to interchange the order of summation and integration, to yield Whence the result.
Before the next results, we need to recall the bivariate Laplace transform, or double Laplace transform, which is applied to bivariate functions f (x, y) in the following way: for p, q ∈ C such that these Laplace transforms exist.
Proof. This follows from the series representation of the operator, using the fact that the series is locally uniformly convergent to interchange the Laplace integration and summation: where we have assumed |λp −α 1 q −α 2 | < 1 for convergence of the binomial series, although this condition can be removed by analytic continuation in the variables p and q.
As an application of the above result, we can use Laplace transforms to quickly learn the result of applying the new fractional operator to the five-parameter Mittag-Leffler function itself, as follows. Example 1. Let α 1 , α 2 , β 1 , β 2 , ε 1 , ε 2 , γ, λ, σ ∈ C with Re(α i ), Re(β i ), Re(ε i ) > 0 for i = 1, 2. Then, we have L 2 a,b I γ;λ α 1 ,α 2 ;β 1 ,β 2 Taking inverse Laplace transforms on both sides of this equation, we get the following formula showing how the bivariate fractional integral operator with five-parameter Mittag-Leffler function kernel can be applied to this particular type of function: A very important property of the bivariate fractional calculus defined in this paper is that it has a semigroup property in the variables β 1 , β 2 , γ, as expressed by the following theorem. There are several different ways to prove this result, as in [18], and we mention here two of them.
Proof. In the case that a = b = 0 and f is a function whose Laplace transform exists, we know from Theorem 4 that applying the operator a,b I γ;λ α 1 ,α 2 ;β 1 ,β 2 corresponds, in the Laplace domain, to multiplication by p −β 1 q −β 2 (1 − λp −α 1 q −α 2 ) −γ . The latter operation clearly has a semigroup property in the parameters β 1 , β 2 , γ, since these appear only as exponents. Thus, the desired result is clear in this case. For the general case, we must proceed by manipulation of infinite series and gamma functions, using Theorem 3: a,b I γ;λ α 1 ,α 2 ;β 1 ,β 2 a,b I σ;λ α 1 ,α 2 ;ε 1 ,ε 2 where for the part in square brackets we use a finite-sum identity on gamma functions (see [12], Theorem 2.9)).

Remark 1.
It should be remarked that the left inverse operator is independent of the parameters ε 1 and ε 2 . The easiest way to see this is by using the series formula: which is independent of ε 1 and ε 2 . Here, we make use of the semigroup property for Riemann-Liouville fractional calculus, in the case where the inner operator is a fractional integral, and we also use the convention (valid by analytic continuation in the order of integration) that a Riemann-Liouville fractional integral to negative order means a Riemann-Liouville fractional derivative, Since we can choose any values of ε 1 and ε 2 with positive real part and get the same left inverse operator, we opt for the values which give ordinary (non-fractional) derivatives of order β 1 + ε 1 and β 2 + ε 2 . This means choosing, just like in Riemann-Liouville fractional calculus, ε 1 = N 1 − β 1 and ε 2 = N 2 − β 2 , to obtain the following definition.
Then, the fractional derivatives of Riemann-Liouville type and Caputo type, with five-parameter Mittag-Leffler kernel, are related to each other as follows: Proof. Starting from the series formulae (13) and (14), we have: where we use both the formula (12) for the double nth integral of the double nth derivative and also the well-known formulae for Riemann-Liouville fractional differintegrals of power functions.

The Non-Singular Cases of the Operators
Recently, the second author [11] made a detailed study of Prabhakar fractional calculus, separating this class of operators into several subclasses according to their properties. Of particular importance is the consideration of whether an operator is singular or nonsingular, and whether the series formula expressing it in terms of Riemann-Liouville integrals is a finite or infinite sum. The same considerations can be applied to other types of fractional calculus, such as the one we are studying here.
In the case of Prabhakar, it was proved [11] (Theorem 4.5) that the most special case is when the operators are non-singular and the sum is finite, and this subclass of Prabhakar fractional calculus consists precisely of the integer powers of the following (inverse) operators: These operators are obtained [47] by setting β = 0 (for non-singular operators) and ρ = ±1 (for finite sums) in the operator P a I α,β,ρ,ω t of Prabhakar integration. We may try to do something similar for the bivariate operator of integration with five-parameter Mittag-Leffler kernel.
Let us firstly compare the Prabhakar fractional integral with the bivariate fractional integral considered in this paper. We have It is clear that our parameters α 1 , α 2 correspond to the α of Prabhakar, our β 1 , β 2 correspond to the β of Prabhakar, our γ corresponds to the ρ of Prabhakar, and our λ corresponds to the ω of Prabhakar, using all notation as above.
The inverse operator is given by setting γ = 1 instead of γ = −1, but this one cannot be written in such an elementary way: a,b I 1;λ α 1 ,α 2 ;0,0 ( f )(x, y) = a,b D −1;λ α 1 ,α 2 ;0,0 ( f )(x, y) = ∂ 2 ∂x∂y a,b I 1;λ α 1 ,α 2 ;1,1 ( f )(x, y) or equivalently, using the series formula (13), The pair of integro-differential operators given by (16) and (17) is of course reminiscent of the so-called Atangana-Baleanu (AB) integral and derivative [48,49]. The latter are defined similarly: the integral by a linear combination of a function with its Riemann-Liouville fractional integral and the derivative by a derivative of an integral transform involving a Mittag-Leffler function kernel. The new development here is that the operators are now bivariate: we can think of the operators (16) and (17) as forming a two-dimensional Atangana-Baleanu calculus.
We also note that there is now a direct connection with the bivariate Abel equation (1) which we considered at first to motivate this paper. The aforementioned Abel equation was rewritten in a form involving the operator (16), and its solution was constructed in a form involving the inverse operator (17).
It is also important to realise that both operators (16) and (17) involve mixed partial integro-differential operators: the AB-type integral operator (16) is defined using fractional integrals with respect to x and y together, and the AB-type derivative operator (17) is defined using mixed partial derivatives with respect to x and y. Thus, these operators may be useful in the study and understanding of fractional PDEs involving mixed partial derivatives, which we note above are under-appreciated in fractional calculus.
To fully recover a bivariate analogue of the AB operators, we need to choose a value for the parameter λ so that appropriate boundary conditions are realised when the parameters α 1 , α 2 go to 0 or 1. This is the reason for the choice of multipliers in the definition of the AB integral and derivative to order α: so that the ordinary first-order integral and derivative are recovered when α → 1 and the original function itself when α → 0. Thus, the following definition is motivated.
The mixed bivariate AB derivatives, of Riemann-Liouville and Caputo type, respectively, are defined to be ABR a,b D α 1 ,α 2 x,y f (x, y) x,y f (x, y) where B(α 1 , α 2 ) is the same normalisation function as above. These definitions are valid for any α 1 , α 2 ∈ C with positive real parts, although the main case of interest is when α 1 , α 2 ∈ (0, 1).

Remark 3.
The multiplier function B(α 1 , α 2 ) is introduced by analogy with the original definition of the AB derivative and AB integral [48], and the restrictions on this function are imposed in order to ensure the following limiting behaviour as the fractional orders α 1 , α 2 approach 0 or 1: Note that if the pair (α 1 , α 2 ) takes the values (0, 1) or (1, 0), then the mixed bivariate AB integral becomes zero. This is why we call it specifically a mixed bivariate operator: in the square set [0, 1] × [0, 1] for values of the pair (α 1 , α 2 ), the values of the mixed bivariate AB integral are weighted towards the diagonal (where both the fractional orders α 1 and α 2 have more similar values) rather than towards the asymmetric corners.
The following result is the bivariate analogue of [47] (Theorem 3.1) or [11] (Proposition 2.4). It follows directly from our consideration above of the non-singular finite-sum special case of the five-parameter Mittag-Leffler kernel operators which culminated in Equations (16) and (17). Proposition 3. The mixed bivariate AB integral and derivatives can be written in terms of the five-parameter Mittag-Leffler kernel operators as follows: However, due to the non-singularity properties of the AB-type operators, they can be defined on a larger function space than the general fractional derivatives with fiveparameter Mittag-Leffler kernel: no differentiability assumptions are required. The following result is the bivariate analogue in [49] (Lemma 2.1) which was the first result to establish appropriate function spaces for the AB derivatives. x,y f (x, y) can also be defined for any function f ∈ L 1 ([a, c] × [b, d]), while its ABC counterpart ABC a,b D α 1 ,α 2 x,y f (x, y) is defined for any twice differentiable function f on [a, c] Proof. For the mixed bivariate AB integral, the result is clear, since this is just a linear combination of f (x, y) with its double Riemann-Liouville integral.
For the mixed bivariate ABC derivative, the result follows from Theorem 2, since this operator is defined by applying a special case of the five-parameter Mittag-Leffler kernel operator to the function ∂ 2 ∂x∂y f (x, y). It remains to consider the mixed bivariate ABR derivative, which we can simplify in the following way since the kernel is non-singular, as a bivariate version of the arguments used in [11,49]: x,y f (x, y) where the new kernel function is and therefore integrable, as x − t → 0 and y − s → 0. Thus, the ABR derivative of f (x, y) equals the function f (x, y) plus an integral term which behaves as a double Riemann-Liouville integral of f (x, y) in the singular limit x − t → 0, y − s → 0. This means the operator is well-defined for any , as required.
The following result is the bivariate analogue of the original series formulae for AB derivatives, [49] (Theorems 2.1 and 2.2). It follows directly from the work we did above to obtain Equation (18). x,y f (x, y) = B(α 1 , α 2 ) x,y f (x, y) = B(α 1 , where this series is locally uniformly convergent for any x,y f (x, y) = ABR a,b D Therefore, the domain of this operator can be extended from the one mentioned in Theorem 8, and the mixed bivariate ABC derivative can also be defined on the whole space L 1 ([a, c] × [b, d]), without any differentiability conditions. x,y f (x, Proof. This is immediate from Proposition 4 and Corollary 1. Note that the series for the ABR derivative consists only of Riemann-Liouville fractional integrals, which have a semigroup property unlike their fractional derivative counterparts.
The following result is the bivariate analogue of the Laplace transform for AB derivatives [48,49].

Proposition 5.
The double Laplace transforms of the mixed bivariate AB integral and derivatives, in the case a = b = 0, can be expressed as follows: Proof. For the mixed bivariate AB integral, the result follows directly from the definition and the known facts on Laplace transforms of Riemann-Liouville fractional integrals. For the mixed bivariate AB derivative of Riemann-Liouville type, the result follows from the series formula of Proposition 4 and a simple application of the binomial theorem. An important fact to notice here is that, via the series formula, this Riemann-Liouville type operator can be written solely in terms of Riemann-Liouville integrals, so there is no need to involve initial conditions here.
For the mixed bivariate AB derivative of Caputo type, derivatives and therefore initial conditions become involved. We use the correct form for the double Laplace transform of a mixed partial derivative, given in [46,50], to achieve the result: which trivially rearranges to the stated result.

Discussion and Conclusions
In this paper, we construct a type of bivariate fractional calculus based originally on a function of a single variable defined by a single power series. We prove some of the fundamental properties of this type of fractional calculus: boundedness of operators, series formulae, Laplace transforms, semigroup and inversion properties, etc. As a special case, we also consider the non-singular versions of our bivariate operators, in order to construct a two-dimensional version of the well-known Atangana-Baleanu calculus.
Some of the ideas used in the above work are similar to those for previously existing models of fractional calculus with various Mittag-Leffler type kernels; however, this is the first time that a univariate function has been used in this way to construct bivariate fractional-calculus operators. Previous contributions in this direction have included using univariate single series to construct univariate operators [10], using bivariate double series to construct univariate operators [18], and using bivariate double series to construct bivariate operators [23,24].
At first glance, it may seem that our construction is unnatural: Given a function of one variable defined by a single series, why would one make a simple thing more complicated by replacing z with x α 1 y α 2 and introducing a bivariate integral operator? The answer is that the underlying mathematical structure is still simple, because everything is defined by a single series whose convergence is easy to describe, but the range of problems that can be tackled is now richer and more diverse, because more variables and parameters allow for more flexibility in adapting the operators to particular scenarios. Our "trick" of turning a univariate function into a bivariate operator gives a shortcut, a possibility of modelling complicated problems using simpler tools.
The bivariate operators that we define lend themselves more to partial differential equations than to ordinary differential equations. Even our initial motivation for proposing them, in the first section of this paper, arose from a bivariate Abel-type equation for a function of two variables. In the literature thus far, Mittag-Leffler kernel operators have mostly been applied to modelling problems with ordinary differential equations. While the existing operators could of course be combined to trivially create bivariate derivatives and integrals with Mittag-Leffler behaviour, we believe that a richer structure will emerge from considering operators such as those here, where the behaviours with respect to x and y are intertwined so that the operator cannot simply be broken down into one with respect to x and another with respect to y. We believe that our work here may have useful ramifications in the study of fractional partial differential equations in 2 or 2+1 dimensions.
For higher dimensions, it will be possible to construct similarly a model of trivariate or multivariate fractional calculus based on a univariate function defined by a single series with 3 or n different gamma functions on the denominator. We give here, without justification, the trivariate version.
The detailed analysis and investigation of these functions and operators, including verification of their essential properties and discussion of special cases and applications, is left for a future research project.
Other related directions of research may include a deeper study of the function spaces on which the integral and derivative operators considered here can be defined. For example, Theorem 2 can be extended easily to show that the bivariate fractional integral operator is bounded on an L p space. Then, a measure-theoretical approach may yield extensions to Morrey spaces [51], or a distributional approach may yield extensions to still larger spaces related to generalised integral operators [52]. Describing more deeply the various relevant function spaces may be useful in the qualitative theory of partial differential equations related to these operators, for example regularity theory [29].
Author Contributions: All authors contributed equally to this article. All authors have read and agreed to the published version of the manuscript.