# The two-dimensional fractional orthogonal derivative.

Enno Diekema 111email adress<EMAIL_ADDRESS>
###### Abstract

This paper is an edited and shortened version of Chapter 6 from the thesis of
the author [4]. First the one dimensional orthogonal derivative will be
extended to the two-dimensional case. In the two-dimensional case we have to
define the region of integration. In this paper we treat the integration over
the square region and over the triangle region where in the last case we use
biorthogonal polynomials expressed in terms of Appell functions. Next the two-
dimensional orthogonal derivative will be extended to the two-dimensional
fractional orthogonal derivative. The results are highly dependent on the
$F_{1}$, $F_{2}$ and $F_{3}$ Appell functions.

##  1 Introduction

In Chapter 2 of the thesis of the author [4, Ch. 2] the orthogonal derivative
of a function of one variable was defined. Just as this definition used
orthogonal polynomials in one variable, we can try to approximate a partial
derivative of a function in several variables by using orthogonal polynomials
in several variables. In this paper we will consider the two-variable case.
Similarly, we will extend the definition in [4, Ch. 3] of the fractional
orthogonal derivative to the two-variable case.
The contents of this paper are as follows.
In section 2 we apply the method of [4, Ch. 2] to the two-variable case. We
want to approximate $\frac{\partial^{n}f(x,y)}{\partial x^{n-k}\partial
y^{k}}$. Only in the trivial case of a direct product of two orthogonality
measures in one variable we can use fully orthogonal polynomials in two
variables (again direct products) for this purpose. For other orthogonality
measures on $\mathbb{R}^{2}$ we have to work with biorthogonal polynomials.
The most obvious regions for the supports of these measures are the triangular
region and the disk. We will only treat the triangular region with orthogonal
measure already considered by Appell and yielding biorthogonal polynomials
expressed in terms of Appell hypergeometric functions.

In section 3 we apply the theory of [4, Ch. 3] to the two-variable case, by
which we obtain a two-dimensional fractional orthogonal derivative.

In section 4 we briefly treat the trivial case of the fractional Jacobi
derivative with weight function on the square region.

In section 5 we treat the fractional biorthogonal derivative with weight
function on the triangular region.

In the appendix we give an overview of the functions we need for the
derivation of the formulas for the fractional orthogonal derivatives. We treat
the $F_{1}$, $F_{2}$, $F_{3}$ Appell functions, the $H_{2}$ Horn function and
the $F_{P}$, $F_{Q}$ and the $F_{PR}$ functions of Olsson.

##  2 The two-dimensional orthogonal derivatives

In [4, Ch. 2] we derived a formula for the orthogonal derivative

$f^{(n)}(x)=\dfrac{k_{n}n!}{h_{n}}\lim\limits_{\delta\downarrow
0}\dfrac{1}{\delta^{n}}\int\nolimits_{\mathbb{R}}f(x+\delta u)p_{n}(u)d\mu(u)$

where the $p_{n}$ are orthogonal polynomials. This formula can be considered,
for instance, for Jacobi, Laguerre or Hermite polynomials. We can also use
discrete orthogonal polynomials, for example the Hahn polynomials.

In the two-variable case we will work with Taylor series, as we did in the
one-dimensional case, but now there is more freedom. Write the Taylor
expansion of a function $f$ around $(x,y)$ as

$f(x+\delta u,y+\delta
v)=\delta^{n}\sum\limits_{m=0}^{n}\sum\limits_{l=0}^{m}\dfrac{\partial^{m}}{\partial
x^{m-l}\partial y^{l}}f(x,y)\dfrac{u^{m-l}v^{l}}{(m-l)!\,l!}+o(\delta^{n})$
(2.1)

where $\delta$ is small. We want to approximate the partial derivative
$\frac{\partial^{n}}{\partial x^{n-k}\partial y^{k}}f(x,y)$. For this purpose
we look for a polynomial $p_{n,k}(u,v)$ of degree $n$ such that, when both
sides of (2.1) are multiplied by $p_{n,k}(u,v)$ and next both sides are
integrated over $\mathbb{R}^{2}$ with respect to a suitable orthogonality
measure $d\mu(u,v)$, all terms in the double summation except for the term
with $(m,l)=(n,k)$ become zero. So first we want

$\int\nolimits_{\mathbb{R}}\int\nolimits_{\mathbb{R}}p_{n,k}(u,v)u^{m-l}v^{l}d\mu(u,v)=0\qquad\text{if}\
\ m<n$

This forces $p_{n,k}$ to be in the ($n+1$)-dimensional space ${\cal V}_{n}$ of
orthogonal polynomials of degree $n$ with respect to the measure $\mu$. Next
we require that

$\int\nolimits_{\mathbb{R}}\int\nolimits_{\mathbb{R}}p_{n,k}(u,v)u^{n-l}v^{l}d\mu(u,v)=0\qquad\text{if}\
\ k\neq l$ (2.2)

This can be equivalently written as

$\int_{\mathbb{R}^{2}}p_{n,k}(u,v)q_{m,l}(u,v)d\mu(u,v)=0\qquad\text{if}\ \
(m,l)\neq(n,k)$

where the polynomials $q_{n,l}$ are elements of ${\cal V}_{n}$ such that
$q_{n,l}(u,v)=u^{n-l}v^{l}$ \+ polynomial of degree less than $n$. The
polynomials $q_{n,l}$ form the so-called monomial or monic basis of ${\cal
V}_{n}$. The above condition on $p_{n,k}\in{\cal V}_{n}$ now defines the
polynomials $p_{n,k}$ uniquely, up to a constant factor, as the basis of
${\cal V}_{n}$ which is biorthogonal to the monomial basis. Then we get from
(2.1) and the properties of $p_{n,k}$ that

$\dfrac{\partial^{n}}{\partial x^{n-k}\partial
y^{k}}f(x,y)=(n-k)!\,k!\lim\limits_{\delta\downarrow
0}\dfrac{1}{\delta^{n}}\dfrac{\int\nolimits_{\mathbb{R}}\int\nolimits_{\mathbb{R}}f(x+\delta
u,y+\delta
v)p_{n,k}(u,v)d\mu(u,v)}{\int\nolimits_{\mathbb{R}}\int\nolimits_{\mathbb{R}}p_{n,k}(u,v)u^{n-k}v^{k}d\mu(u,v)}$
(2.3)

under some conditions on $f$, for instance that $f$ is a $C^{n}$ function on
some neighbourhood of $(x,y)$ and, if $\mu$ has unbounded support, that
$|f(u,v)|$ tends sufficiently fast to $0$ as $(u,v)$ tends to $\infty$. More
generally, we call the right-hand side of (2.3) a two-dimensional partial
orthogonal derivative of $f$ at $(x,y)$ if the limit exists. But we will in
particular be interested in the numerator on the right-hand side (after the
limit sign) for fixed $\delta$. We will only treat the square and the
triangular region. For the disk see [6, Section 5.2.2 with $d=2$], [16,
Chapter 3].

###  2.1 The two-dimensional Jacobi derivative with weight function on the
square region

An easy example of (2.3) can be obtained by letting
$d\mu(x,y)=d\mu_{1}(x)d\mu_{2}(y)$ and $p_{n,k}(x,y)=r_{n-k}(x)s_{k}(y)$ with
the $r_{n}$ and $s_{n}$ one-variable orthogonal polynomials for orthogonality
measures $\mu_{1}$ and $\mu_{2}$, respectively. Then the polynomials $p_{n,k}$
are two-variable orthogonal polynomials for the orthogonality measure $\mu$
while also $p_{n,k}(x,y)=$C$\,x^{n-k}y^{k}$ \+ polynomial of degree less than
$n$. So the $p_{n,k}$ form the monomial basis, but this basis is also
orthogonal, This is also biorthogonal to the monomial basis, by which the
$p_{n,k}$ are the polynomials desired in (2.3). Now the denominator after the
limit sign on the right-hand side of (2.3) factorizes as the product

$\int_{\mathbb{R}}r_{n-k}(u)u^{n-k}d\mu_{1}(u)\times\int_{\mathbb{R}}s_{k}(v)v^{k}d\mu_{2}(v)=\frac{h_{n-k}^{\prime}}{k_{n-k}^{\prime}}\frac{h_{k}^{\prime\prime}}{k_{k}^{\prime\prime}}$

by [4, (2.2.4)], where the accented and double accented $h_{n}$’s and
$k_{n}$’s are in the obvious way, following [4, (2.2.3)], related to the
orthogonal polynomials $r_{n}$ and $s_{n}$, respectively. Specialization of
(2.3) gives:

$\dfrac{\partial^{m+n}}{\partial^{m}x\partial^{n}y}f(x,y)=\dfrac{k_{m}\,m!}{h_{m}}\dfrac{k_{n}\,n!}{h_{n}}\lim\limits_{\delta\rightarrow
0}\dfrac{1}{\delta^{m+n}}\int\nolimits_{-1}^{1}\int\nolimits_{-1}^{1}f(x+\delta
v,y+\delta u)p_{m}(u)q_{n}(v)d\mu(u)d\mu(v)$ (2.4)

This is a very trivial case. If the measures $\mu_{1}$ and $\mu_{2}$ are
supported on the segment $[-1,1]$ then the polynomials $p_{n,k}$ are
orthogonal on the square $[-1,1]\times[-1,1]$. For example, take the product
of the Jacobi measures $d\mu_{1}(x)=(1-x)^{\alpha}(1+x)^{\beta}dx$ and
$d\mu_{2}(y)=(1-y)^{\gamma}(1+y)^{\delta}dy$. Substitution in (2.4) gives at
last

$\dfrac{\partial^{n+k}}{\partial x^{n}\partial
y^{k}}f(x,y)=\lim\limits_{\delta\rightarrow
0}\dfrac{1}{\delta^{n+k}}\dfrac{k_{n,\alpha,\beta}\
k_{k,\gamma,\delta}}{h_{n,\alpha,\beta}\ h_{k,\gamma,\delta}}\\\
\times\int\nolimits_{-1}^{1}\int\nolimits_{-1}^{1}f(x+\delta u,y+\delta
v)P_{n}^{(\alpha,\beta)}(u)P_{k}^{(\gamma,\delta)}(v)(1-u)^{\alpha}(1+u)^{\beta}(1-v)^{\gamma}(1+v)^{\delta}\,du\,dv$

###  2.2 The two-dimensional biorthogonal derivative with weight function on
the triangular region

The triangular region of ${\mathbb{R}}^{2}$ is defined as
$T^{2}:=\left\\{(x,y):0\leq x,y\leq 1\wedge x+y\leq 1\right\\}$, on which we
consider the Jacobi weight function defined as:

$W_{\alpha,\beta,\gamma}(x,y):=\dfrac{1}{B(\alpha+1,\beta+1,\gamma+1)}x^{\alpha}y^{\beta}(1-x-y)^{\gamma}\qquad\alpha,\beta,\gamma>-1$

which is normalized in such a way that its integral over $T^{2}$ is $1$. The
function $B(x,y,z)$ is defined as an extension of the Beta function

$B(x,y,z):=\dfrac{\Gamma(x)\Gamma(y)\Gamma(z)}{\Gamma(x+y+z)}$

It is easy to prove

$B(x,y,z)=B(x,y+z)B(y,z)=B(y,x+z)B(x,z)=B(z,x+y)B(x,y)$

Let

$\left\langle
f,g\right\rangle_{\alpha,\beta,\gamma}=\int\nolimits_{T^{2}}f(x,y)g(x,y)W_{\alpha,\beta,\gamma}(x,y)\;dx\;dy$

Classical orthogonal polynomials of one variable all satisfy a Rodrigues’
formula. For obtaining something similar in connection with the triangular
region Appell considered the following analogue of Rodrigues’ formula [6,
Section 2.4]:

$U_{k,n}^{\alpha,\beta,\gamma}(x,y)=\big{[}x^{\alpha}y^{\beta}(1-x-y)^{\gamma}\big{]}^{-1}\dfrac{\partial^{n}}{\partial
x^{k}\partial
y^{n-k}}\left[x^{k+\alpha}y^{n-k+\beta}(1-x-y)^{n+\gamma}\right]\qquad 0\leq
k\leq n$ (2.5)

The set $\\{U_{k,n}^{\alpha,\beta,\gamma}(x,y):0\leq k\leq n\\}$ is a basis of
$\mathcal{V}_{n}^{2}(W_{\alpha,\beta,\gamma})$ where $\mathcal{V}_{n}^{2}$ is
the space of orthogonal polynomials of degree $n$. See [6, Prop. 2.4.3].

To write $U_{k,n}^{\alpha,\beta,\gamma}(x,y)$ as an $F_{2}$ Appell function we
can use [7, 5.13(1)]. There follows:

$U_{k,n}^{\alpha,\beta,\gamma}(x,y)=(\alpha+1)_{k}(\beta+1)_{n-k}(1-x-y)^{n}\\\
F_{2}\left(-\gamma-n;-k,-n+k;\alpha+1,\beta+1;\dfrac{x}{x+y-1},\dfrac{y}{x+y-1}\right)$

From this expression we see that $U_{k,n}^{\alpha,\beta,\gamma}(x,y)$ are
polynomials of degree less or equal $n$. Using [7, 5.11(8)] we get the much
simpler form

$U_{k,n}^{\alpha,\beta,\gamma}(x,y)=(\alpha+1)_{k}(\beta+1)_{n-k}(1-x-y)^{-\gamma}\\\
F_{2}(-\gamma-n;\alpha+1+k,\beta+1+n-k;\alpha+1,\beta+1;x,y)$ (2.6)

This basis is not orthogonal. It can be shown that $\left\langle
U_{k,n}^{\alpha,\beta,\gamma},U_{j,n}^{\alpha,\beta,\gamma}\right\rangle_{\alpha,\beta,\gamma}\neq
0$ for $j\neq k$. If $n\neq m$ then $\left\langle
U_{k,n}^{\alpha,\beta,\gamma},U_{j,m}^{\alpha,\beta,\gamma}\right\rangle_{\alpha,\beta,\gamma}=0$.

Let $V_{m,n}^{\alpha,\beta,\gamma}$ be defined with $0\leq m\leq n$ by

$V_{m,n}^{\alpha,\beta,\gamma}(x,y):=\\\
=\sum\limits_{i=0}^{m}\sum\limits_{j=0}^{n}(-1)^{m+n+i+j}\dbinom{m}{i}\dbinom{n}{j}\dfrac{(\alpha+1)_{m}(\beta+1)_{n}(\alpha+\beta+\gamma+2)_{m+n+i+j}}{(\alpha+1)_{i}(\beta+1)_{j}(\alpha+\beta+\gamma+2)_{2m+2n}}x^{i}y^{j}$
(2.7)

This function can be written as an $F_{2}$ Appell function.

$V_{m,n}^{\alpha,\beta,\gamma}(x,y)=(-1)^{n+m}\dfrac{(\alpha+1)_{m}(\beta+1)_{n}(\alpha+\beta+\gamma+2)_{m+n}}{(\alpha+\beta+\gamma+2)_{2n+2m}}\\\
F_{2}(\alpha+\beta+\gamma+2+m+n;-m,-n;\alpha+1,\beta+1;x,y)$ (2.8)

Then $V_{m,n}$ is the orthogonal projection of $x^{m}y^{n}$ in
$L^{2}\left\langle T^{2},W_{\mu}\right\rangle$ onto the space
$\mathcal{V}_{m+n}$ of orthogonal polynomials of degree $m+n$ (so they are a
monomial basis) and the two families of polynomials are biorthogonal with

$\left\langle
U_{k,n}^{\alpha,\beta,\gamma},V_{j,n-j}^{\alpha,\beta,\gamma}\right\rangle_{\alpha,\beta,\gamma}=(-1)^{n}\dfrac{(\alpha+1)_{k}(\beta+1)_{n-k}(\gamma+1)_{n}k!\,(n-k)!}{(\alpha+\beta+\gamma+3)_{2n}}\delta_{j,k}\qquad
0\leq j,k\leq n$ (2.9)

We can write this last inner product more generally by replacing the parameter
$n-j$ by $m-j$. There follows

$\left\langle
U_{k,n}^{\alpha,\beta,\gamma},V_{j,m-j}^{\alpha,\beta,\gamma}\right\rangle_{\alpha,\beta,\gamma}=0\qquad
m\neq n$

The basis of polynomials $V_{m,n}^{\alpha,\beta,\gamma}$ as well as the above
biorthogonality is already given by Appell for $\gamma=0$ and by Fackerel &
Littler [9] for the general case. If we compare with the notation for the
general case at the beginning of Section 2 then we see that
$U_{k,n}^{\alpha,\beta,\gamma}=p_{n,n-k}$ and
$V_{l,m-l}^{\alpha,\beta,\gamma}=q_{m,m-l}$. Hence, for the partial orthogonal
derivative (2.3) we have to use the $U$-basis in the present case.

Using the triangular region $T^{2}=\left\\{(x,y):0\leq x,y,x+y<1\right\\}$ we
get for the denominator of (2.3)

$\displaystyle I(m,l)$ $\displaystyle=\iint\nolimits_{0\leq
u,v,u+v<1}U_{k,n}^{\alpha,\beta,\gamma}(u,v)u^{m-l}v^{l}W_{\alpha,\beta,\gamma}(u,v)dudv=$
$\displaystyle=\big{[}B(\alpha+1,\beta+1,\gamma+1)\big{]}^{-1}$
$\displaystyle\qquad\qquad\times\int\nolimits_{0}^{1}v^{l}\int\nolimits_{0}^{1-v}u^{m-l}\dfrac{\partial^{k}}{\partial
u^{k}}\dfrac{\partial^{n-k}}{\partial
v^{n-k}}\left[u^{k+\alpha}v^{n-k+\beta}(1-u-v)^{n+\gamma}\right]\,du\,dv$
(2.10)

By integration by parts with respect to $u$ we see that $I\left(m,l\right)=0$
unless $l+k\leq m$. By integration by parts with respect to $v$ we see that
$I(m,l)=0$ unless $n\leq l+k$. Together with $m\leq n$ there follows: $m=n$
and $l=n-k$. Substitution in (2.10) gives

$\displaystyle I(n,n-k)$ $\displaystyle B(\alpha+1,\beta+1,\gamma+1)=$
$\displaystyle=\int\nolimits_{0}^{1}v^{n-k}\int\nolimits_{0}^{1-v}u^{k}\dfrac{\partial^{k}}{\partial
u^{k}}\dfrac{\partial^{n-k}}{\partial
v^{n-k}}\left[u^{k+\alpha}v^{n-k+\beta}(1-u-v)^{n+\gamma}\right]\,du\,dv$
(2.11)

Let

$I_{1}=\int\nolimits_{0}^{1-v}u^{k}\dfrac{\partial^{k}}{\partial
u^{k}}\left[\dfrac{\partial^{n-k}}{\partial
v^{n-k}}\left[u^{k+\alpha}v^{n-k+\beta}(1-u-v)^{n+\gamma}\right]\right]\,du$

After $k$-fold integration by parts there remains

$I_{1}=(-1)^{k}\Gamma(k+1)\int\nolimits_{0}^{1-v}u^{k+\alpha}\dfrac{\partial^{n-k}}{\partial
v^{n-k}}\left[v^{n-k+\beta}(1-u-v)^{n+\gamma}\right]du$

Substitution in (2.11) gives

$\displaystyle I(n,n-k)$
$\displaystyle\big{[}(-1)^{k}\Gamma(k+1)\big{]}^{-1}B(\alpha+1,\beta+1,\gamma+1)=$
$\displaystyle=\int\nolimits_{0}^{1}v^{n-k}\int\nolimits_{0}^{1-v}u^{k+\alpha}\dfrac{\partial^{n-k}}{\partial
v^{n-k}}\left[v^{n-k+\beta}(1-u-v)^{n+\gamma}\right]\,du\,dv$
$\displaystyle=\int\nolimits_{0}^{1}u^{k+\alpha}\int\nolimits_{0}^{1-u}v^{n-k}\dfrac{\partial^{n-k}}{\partial
v^{n-k}}\left[v^{n-k+\beta}(1-u-v)^{n+\gamma}\right]\,dv\,du$ (2.12)

Let

$I_{2}=\int\nolimits_{0}^{1-u}v^{n-k}\dfrac{\partial^{n-k}}{\partial
v^{n-k}}\left[v^{n-k+\beta}(1-u-v)^{n+\gamma}\right]\,dv$

After $(n-k)$-fold integration by parts we get

$I_{2}=(-1)^{n-k}\Gamma(n-k+1)\int\nolimits_{0}^{1-u}v^{n-k+\beta}(1-u-v)^{n+\gamma}dv$

The integral is a Beta function.

$\int\nolimits_{0}^{1-u}v^{n-k+\beta}(1-u-v)^{n+\gamma}dv=(1-u)^{k+2n+\beta+\gamma+1}B(n-k+\beta+1,n+\gamma+1)$

with convergence conditions: $\beta,\gamma>-1$. So for $I_{2}$ we get

$I_{2}=(-1)^{n-k}\Gamma(n-k+1)(1-u)^{k+2n+\beta+\gamma+1}B(n-k+\beta+1,n+\gamma+1)$

Substitution in (2.12) gives

$I(n,n-k)\big{[}(-1)^{n}\Gamma(k+1)\Gamma(n-k+1)\big{]}^{-1}B(\alpha+1,\beta+1,\gamma+1)=\\\
=B(n-k+\beta+1,n+\gamma+1)\int\nolimits_{0}^{1}u^{k+\alpha}(1-u)^{k+2n+\beta+\gamma+1}du$

The integral is again a Beta function. So we get

$\displaystyle I(n,n-k)\big{[}(-1)^{n}$
$\displaystyle\Gamma(k+1)\Gamma(n-k+1)\big{]}^{-1}B(\alpha+1,\beta+1,\gamma+1)=$
$\displaystyle=B(n-k+\beta+1,n+\gamma+1)B(k+\alpha+1,2n-k+\beta+\gamma+2)$
$\displaystyle=B(k+\alpha+1,n-k+\beta+1,n+\gamma+1)$

with convergence conditions:
$\operatorname{Re}{\alpha},\operatorname{Re}{\beta},\operatorname{Re}{\gamma}>-1$.
After some manipulations we get

$I(n,n-k)=(-1)^{n}\dfrac{(\alpha+1)_{k}(\beta+1)_{n-k}(\gamma+1)_{n}k!\,(n-k)!}{(\alpha+\beta+\gamma+3)_{2n}}$

and this is equal to the constant in (2.9). This also follows from the fact
that

$\displaystyle I(n,n-k)$ $\displaystyle=<U_{k,n},\,u^{k}\,v^{n-k}>$
$\displaystyle=<U_{k,n},\,V_{k,n-k}+\text{polynomial of degree less than
$n$}>$ $\displaystyle=<U_{k,n},\,V_{k,n-k}>+<U_{k,n},\text{polynomial of
degree less than $n$}>$ $\displaystyle=<U_{k,n},\,V_{k,n-k}>$

Substitution in (2.3) gives

$\dfrac{\partial^{n}}{\partial x^{k}\partial
y^{n-k}}f(x,y)=\dfrac{(-1)^{n}}{B(\alpha+1+k,\beta+1+n-k,\gamma+1+n)}\\\
\lim\limits_{\delta\downarrow 0}\dfrac{1}{\delta^{n}}\iint\nolimits_{0\leq
u,v,u+v<1}f(x+\delta u,y+\delta
v)U_{k,n}^{\alpha,\beta,\gamma}(u,v)u^{\alpha}v^{\beta}(1-u-v)^{\gamma}\,dv\,du$
(2.13)

with convergence conditions: $\alpha,\beta,\gamma>-1$. Using (2.6) gives

$\dfrac{\partial^{n}}{\partial x^{k}\partial
y^{n-k}}f(x,y)=\dfrac{(-1)^{n}\Gamma(\alpha+\beta+\gamma+3+2n)}{\Gamma(\alpha+1)\Gamma(\beta+1)\Gamma(\gamma+1+n)}\\\
\lim\limits_{\delta\downarrow 0}\dfrac{1}{\delta^{n}}\iint\nolimits_{0\leq
u,v,u+v<1}f(x+\delta u,y+\delta
v)F_{2}\left(\begin{array}[]{c}-\gamma-n;\alpha+1+k,\beta+1+n-k\\\
\alpha+1,\beta+1\end{array};u,v\right)u^{\alpha}v^{\beta}\,dv\,du$

##  3 The two-dimensional fractional orthogonal derivative

We have seen that for the orthogonal derivative in two dimensions there are
several choices of region and orthogonality measures for which explicit
results can be obtained. For the fractional derivative in two dimensions we
have the same possibilities. So we treat first the trivial case of the
fractional orthogonal derivative with weight function on the square. Then we
treat the fractional biorthogonal derivative with weight function on the
triangular region.

Analogous to the Weyl fractional integral in one dimension we define the Weyl
fractional integral in two dimensions as:

$W^{-\mu,-\nu}[f](x,y):=\dfrac{1}{\Gamma(\mu)\Gamma(\nu)}\int\nolimits_{x}^{\infty}\int\nolimits_{y}^{\infty}f(U,V)(U-x)^{\mu-1}(V-y)^{\nu-1}dVdU$

The fractional partial derivative will be given by

$W^{\mu,\nu}[f](x,y)=(-1)^{m}\dfrac{\partial^{m}}{\partial x^{m-l}\partial
y^{l}}W^{\mu-m+l,}{}^{\nu-l}[f](x,y)\\\
=\dfrac{(-1)^{m}}{\Gamma(m-l-\mu)\Gamma(l-\nu)}\dfrac{\partial^{m}}{\partial
x^{m-l}\partial
y^{l}}\left[\int\nolimits_{x}^{\infty}\int\nolimits_{y}^{\infty}f(U,V)(U-x)^{m-l-\mu-1}(V-y)^{l-\nu-1}dVdU\right]$
(3.1)

Analogous to [4, (3.3.1)] and in view of (2.3) we can write

$W^{\mu,\nu}[f](x,y)=\lim\limits_{\delta\downarrow
0}W^{\mu,\nu,m,l}_{\delta}[f](x,y)$

where

$\displaystyle W^{\mu,\nu,m,l}_{\delta}[f](x,y):$
$\displaystyle=(-1)^{m}D^{m,l}_{\delta}\big{[}W^{\mu-m+l,\nu-l}[f]\big{]}(x,y)$
$\displaystyle=(-1)^{m}W^{\mu-m+l,\nu-l}\big{[}D^{m,l}_{\delta}[f]\big{]}(x,y)$
(3.2)

with

$D^{m,l}_{\delta}[g](x,y):=\dfrac{(m-l)!\,l!}{\delta^{m}}\dfrac{}{}\dfrac{\int\nolimits_{\mathbb{R}}\int\nolimits_{\mathbb{R}}g(x+\delta
u,y+\delta
v)p_{m,l}(u,v)d\mu(u,v)}{\int\nolimits_{\mathbb{R}}\int\nolimits_{\mathbb{R}}p_{m,l}(u,v)u^{m-l}v^{l}d\mu(u,v)}$
(3.3)

After application of (3.2) and (3.3) to (3.1) we obtain

$W^{\mu,\nu,m,l}_{\delta}[f](x,y)=\dfrac{(-1)^{m}}{\Gamma(m-l-\mu)\Gamma(l-\nu)}\dfrac{(m-l)!\,l!}{\int\nolimits_{\mathbb{R}}\int\nolimits_{\mathbb{R}}p_{m,l}(u,v)u^{m-l}v^{l}d\mu(u,v)}\dfrac{1}{\delta^{m}}\\\
\times\int\nolimits_{\mathbb{R}}\int\nolimits_{\mathbb{R}}\left[\int\nolimits_{x+\delta
u}^{\infty}\int\nolimits_{y+\delta v}^{\infty}f(U,V)(U-x-\delta
u)^{m-l-\mu-1}(V-y-\delta v)^{l-\nu-1}dVdU\right]\\\ \times
p_{m,l}(u,v)d\mu(u,v)$

With $U=x+\delta s$ and $V=y+\delta t$ there follows

$W^{\mu,\nu,m,l}_{\delta}[f](x,y)=\dfrac{(-1)^{m}}{\Gamma(m-l-\mu)\Gamma(l-\nu)}\dfrac{(m-l)!\,l!}{\int\nolimits_{\mathbb{R}}\int\nolimits_{\mathbb{R}}p_{m,l}(u,v)u^{m-l}v^{l}d\mu(u,v)}\dfrac{1}{\delta^{\mu+\nu}}\\\
\times\int\nolimits_{\mathbb{R}}\int\nolimits_{\mathbb{R}}\left[\int\nolimits_{x}^{\infty}\int\nolimits_{y}^{\infty}f(x+\delta
s,y+\delta t)(s-u)^{m-l-\mu-1}(t-v)^{l-\nu-1}dtds\right]p_{m,l}(u,v)d\mu(u,v)$

Interchange of the inner double integral with the outer double integral gives

$W^{\mu,\nu,m,l}_{\delta}[f](x,y)=\dfrac{(-1)^{m}}{\Gamma(m-l-\mu)\Gamma(l-\nu)}\dfrac{(m-l)!\,l!}{\int\nolimits_{\mathbb{R}}\int\nolimits_{\mathbb{R}}p_{m,l}(u,v)u^{m-l}v^{l}d\mu(u,v)}\dfrac{1}{\delta^{\mu+\nu}}\\\
\times\int\nolimits_{\mathbb{R}}\int\nolimits_{\mathbb{R}}f(x+\delta
s,y+\delta
t)\left[\int\nolimits_{u=-\infty}^{s}\int\nolimits_{v=-\infty}^{t}(s-u)^{m-l-\mu-1}(t-v)^{l-\nu-1}p_{m,l}(u,v)d\mu(u,v)\right]dsdt$
(3.4)

In subsequent subsections we will write (3.4) in more explicit form for the
cases of the square and the triangular regions.

##  4 The two-dimensional fractional Jacobi derivative with weight function on
the square.

Let $d\mu(x,y)=d\mu_{1}(x)d\mu_{2}(y)$ and $p_{n,k}(x,y)=r_{n-k}(x)s_{k}(y)$
as in the beginning of section 2.1. Then by (2.4) formula (3.4) takes the form

$W^{\mu,\nu,m,l}_{\delta}[f](x,y)=\dfrac{(-1)^{m}(m-l)!\ l!\ k^{\prime}_{m-l}\
k^{\prime\prime}_{l}}{\Gamma(m-l-\mu)\Gamma(l-\nu)h^{\prime}_{m-l}\
h^{\prime\prime}_{l}}\dfrac{1}{\delta^{\mu+\nu}}\int\nolimits_{\mathbb{R}}\int\nolimits_{\mathbb{R}}f(x+\delta
s,y+\delta t)\\\
\times\left[\int\nolimits_{u=-\infty}^{s}(s-u)^{m-l-\mu-1}r_{m-l}(u)d\mu_{1}(u)\right]\left[\int\nolimits_{v=-\infty}^{t}(t-v)^{l-\nu-1}s_{l}(v)d\mu_{2}(v)\right]dsdt$

Now specialize to $d\mu_{1}(x)=(1-x)^{\alpha}(1+x)^{\beta}dx$ on $[-1,1]$ for
$r_{n}(x)=P^{(\alpha,\beta)}_{n}(x)$ and
$d\mu_{2}(y)=(1-y)^{\gamma}(1+y)^{\delta}\ dy$ for
$s_{n}(y)=P^{(\gamma,\delta)}_{n}(y)$. Then we get for the fractional Jacobi
derivative

$W^{\mu,\nu,m,l}_{\delta}[f](x,y)=\dfrac{(-1)^{m}(m-l)!\ l!\ k^{\prime}_{m-l}\
k^{\prime\prime}_{l}}{\delta^{\mu+\nu}\ h^{\prime}_{m-l}\
h^{\prime\prime}_{l}}\\\
\times\int\nolimits_{-1}^{1}\int\nolimits_{-1}^{1}f(x+\delta s,y+\delta
t)\left[\int\nolimits_{-\infty}^{s}(s-u)^{m-l-\mu-1}P^{(\alpha,\beta)}_{m-l}(u)(1-u)^{\alpha}(1+u)^{\beta}du\right]\\\
\times\left[\int\nolimits_{-\infty}^{t}(t-v)^{l-\nu-1}P^{(\gamma,\delta)}_{l}(v)(1-v)^{\gamma}(1+v)^{\delta}\
dv\right]dsdt\qquad$ (4.1)

We now look at the two-dimensional analogue of [4, (3.3.6)]. There we had
transformed the 2-fold integral into a sum of two single integrals, where the
integrands are a product of the function $f(x+\delta s)$ with a factor
depending on a convolution integral of an orthogonal polynomial. In the two-
dimensional case we can transform the $4$-fold integral in such a manner that
there results a sum of four $2$-fold integrals with a product of the function
$f(x+\delta s,y+\delta t)$ and a factor depending on the $2$-fold convolution
integral of the orthogonal polynomials. In the special case of a product of
Jacobi polynomials we can use the same computation as in the one-dimensional
case. We give here only the result.

$\displaystyle W^{\mu,\nu,m+n,n}[f](x,y)$
$\displaystyle\left[(-1)^{m+n}\dfrac{k_{m}m!}{h_{m}}\dfrac{k_{n}n!}{h_{n}}\dfrac{1}{\delta^{\mu+\nu}}\dfrac{1}{\Gamma(m-\mu)\Gamma(n-\nu)}\right]^{-1}=$
$\displaystyle=\int\nolimits_{-1}^{1}\int\nolimits_{-1}^{1}f(x+\delta
s,y+\delta t)J_{1}(t,n,\alpha,\beta,\nu)J_{1}(s,m,\gamma,\delta,\mu)\,dt\,ds+$
$\displaystyle+\int\nolimits_{1}^{\infty}\int\nolimits_{-1}^{1}f(x+\delta
s,y+\delta t)J_{1}(t,n,\alpha,\beta,\nu)J_{2}(s,m,\gamma,\delta,\mu)\,dt\,ds+$
$\displaystyle+\int\nolimits_{-1}^{1}\int\nolimits_{1}^{\infty}f(x+\delta
s,y+\delta t)J_{2}(t,n,\alpha,\beta,\nu)J_{1}(s,m,\gamma,\delta,\mu)\,dt\,ds+$
$\displaystyle+\int\nolimits_{1}^{\infty}\int\nolimits_{1}^{\infty}f(x+\delta
s,y+\delta t)J_{2}(t,n,\alpha,\beta,\nu)J_{2}(s,m,\gamma,\delta,\mu)\,dt\,ds$

The integrals $J_{1}$ and $J_{2}$ are known [4, (3.5.9)] and [4, (3.5.10)]. We
obtain

$J_{1}(\xi,\lambda,\alpha,\beta,\nu)=(-1)^{\lambda}\dfrac{\Gamma(\lambda+\beta+1)\Gamma(\lambda-\nu)}{2^{\lambda-\nu}\lambda!\,\Gamma(\lambda-\nu+\beta+1)}(1-\xi)^{\lambda+\alpha-\nu}(1+\xi)^{\lambda+\beta-\nu}\\\
\,\mbox{}_{2}F_{1}\\!\left(\genfrac{}{}{0.0pt}{}{-\nu,2\lambda-\nu+\alpha+\beta+1}{\lambda-\nu+\beta+1};\dfrac{1+\xi}{2}\right)$
$J_{2}(\xi,\lambda,\alpha,\beta,\nu)=(-1)^{\lambda}\dfrac{2^{\lambda+\alpha+\beta+1}}{(\xi+1)^{\nu+1}}\dfrac{\Gamma(\lambda-\nu)}{\Gamma(-\nu)\lambda!}\dfrac{\Gamma(\lambda+\alpha+1)\Gamma(\lambda+\beta+1)}{\Gamma(2\lambda+\alpha+\beta+2)}\\\
\,\mbox{}_{2}F_{1}\\!\left(\genfrac{}{}{0.0pt}{}{\nu+1,\lambda+\beta+1}{2\lambda+\alpha+\beta+2};\dfrac{2}{\xi+1}\right)$

##  5 The fractional biorthogonal derivative with weight function on the
triangular region

To develop a formula for the fractional biorthogonal derivative on the
triangular region, we use (3.4) for the triangular region $T^{2}$.

$W^{\mu,\nu,n,n-k}_{\delta}[f](x,y)\left[\dfrac{h_{1}(\alpha,\beta,\gamma,k,n)}{\Gamma(k-\mu)\Gamma(n-k-\nu)}\dfrac{1}{\delta^{n}}\right]^{-1}=\int\nolimits_{\mathbb{R}}\int\nolimits_{\mathbb{R}}f(x+\delta
s,y+\delta t)\\\
\times\left[\iint\nolimits_{T^{2}\cap\big{(}(-\infty,s]\times(-\infty,t]\big{)}}(s-u)^{k-\mu-1}(t-v)^{n-k-\nu-1}U^{\alpha,\beta,\gamma}_{k,n}(u,v)u^{\alpha}v^{\beta}(1-u-v)^{\gamma}dudv)\right]dsdt$
(5.1)

with

$h_{1}(\alpha,\beta,\gamma,k,n)=\dfrac{1}{B(\alpha+1+k,\beta+1+n-k,\gamma+1+n)}$
(5.2)

We see that the integration domain $T^{2}\cap((-\infty,s]\times(-\infty,t])$
of the inner double integral in (5.1) depends on the position of the point
$(s,t)$ with respect to the triangle $\left\\{\left(u,v\right)\text{ }|\text{
}0<v<u<1\right\\}$. Clearly $T^{2}\cap((-\infty,s]\times(-\infty,t])$ is empty
if $(s,t)$ is outside the first quadrant, while the first quadrant in the
$(s,t)$ plane can be split up as the union of five regions with disjoint
interiors as in Figure 1 such that the corresponding regions
$T^{2}\cap((-\infty,s]\times(-\infty,t])$ inside $T^{2}$ in the $(u,v)$ plane
are given by Figure 2.

Figure 1: The five possible regions of the position of the point $(s,t)$.
Figure 2: The five possibilities of the position of the point $(s,t)$ with
respect to the $u,v$ axes.

I.  The point $(s,t)$ lies inside the triangle with $0<s<1,0<t<1-s$. The
appropriate part of the integral for $W$ becomes

$I=\int_{s=0}^{1}\int_{t=0}^{1-s}\left(\int_{v=0}^{t}\int_{u=0}^{s}(.)\,du\,dv\right)\,dt\,ds$

II.  The point $(s,t)$ lies outside the triangle with $1<s<\infty,1<t<\infty$.
The appropriate part of the integral for $W$ becomes

$I=\int_{s=1}^{\infty}\int_{t=1}^{\infty}\left(\int_{v=0}^{1}\int_{u=0}^{1-v}(.)\,du\,dv\right)\,dt\,ds$

III.  The point $(s,t)$ lies outside the triangle with $0<s<1<t<\infty$. The
appropriate part of the integral for $W$ becomes

$I=\int_{t=1}^{\infty}\int_{s=0}^{1}\left(\int_{u=0}^{s}\int_{v=0}^{1-u}(.)\,dv\,du\right)\,ds\,dt$

IV.  The point $(s,t)$ lies outside the triangle with $0<t<1<s<\infty$. The
appropriate part of the integral for $W$ becomes

$I=\int_{s=1}^{\infty}\int_{t=0}^{1}\left(\int_{v=0}^{t}\int_{u=0}^{1-v}(.)\,du\,dv\right)\,dt\,ds$

V.  The point $(s,t)$ lies outside the triangle with $0<s<1,\ 1-s<t<1$. We
divide the region of integration in two parts. The appropriate part of the
integral for $W$ becomes

$I=\int_{s=0}^{1}\int_{t=1-s}^{1}\left(\int_{v=0}^{1-s}\int_{u=0}^{s}(.)\,du\,dv\right)\,dt\,ds+\int_{s=0}^{1}\int_{t=1-s}^{1}\left(\int_{v=1-s}^{t}\int_{u=0}^{1-v}(.)\,du\,dv\right)\,dt\,ds$

For the fractional derivative we get (the roman numbers refer to the cases in
Figure 2)

$\displaystyle W^{\mu,\nu}$
$\displaystyle[f](x,y)\left[\dfrac{h_{1}(\alpha,\beta,\gamma,k,n)}{\Gamma(k-\mu)\Gamma(n-k-\nu)}\dfrac{1}{\delta^{\nu+\mu}}\right]^{-1}=$
$\displaystyle=\int_{s=0}^{1}\int_{t=0}^{1-s}f(x+\delta s,y+\delta
t)\left(\int_{v=0}^{t}\int_{u=0}^{s}\text{ \
}(.)\,du\,dv\right)\,dt\,ds+\qquad\qquad\quad\ \text{(I)}$
$\displaystyle+\int_{s=1}^{\infty}\int_{t=1}^{\infty}f(x+\delta s,y+\delta
t)\left(\int_{v=0}^{1}\int_{u=0}^{1-v}(.)\,du\,dv\right)\,dt\,ds+\qquad\qquad\qquad\
\text{(II)}$ $\displaystyle+\int_{t=1}^{\infty}\int_{s=0}^{1}\text{
}f(x+\delta s,y+\delta
t)\left(\int_{u=0}^{s}\int_{v=0}^{1-u}(.)\,dv\,du\right)\,ds\,dt+\qquad\qquad\quad\
\ \text{(III)}$ $\displaystyle+\int_{s=1}^{\infty}\int_{t=0}^{1}f(x+\delta
s,y+\delta
t)\left(\int_{v=0}^{t}\int_{u=0}^{1-v}(.)\,du\,dv\right)\,dt\,ds+\qquad\qquad\qquad\text{
(IV)}$ $\displaystyle+\int_{s=0}^{1}\int_{t=1-s}^{1}f(x+\delta s,y+\delta
t)\left(\int_{v=0}^{1-s}\int_{u=0}^{s}(.)\,du\,dv\right)\,dt\,ds+\qquad\qquad\quad\
\text{ (V)}$
$\displaystyle\quad\quad\quad+\int_{s=0}^{1}\int_{t=1-s}^{1}f(x+\delta
s,y+\delta t)\left(\int_{v=1-s}^{t}\int_{u=0}^{1-v}(.)\,du\,dv\right)\,dt\,ds$

The abbreviation $(.)$ stands for

$(.)=(s-u)^{k-\mu-1}(t-v)^{n-k-\nu-1}U_{k,n}^{\alpha,\beta,\gamma}(u,v)u^{\alpha}v^{\beta}(1-u-v)^{\gamma}$
(5.4)

For $U_{k,n}^{\alpha,\beta,\gamma}(u,v)$ we repeat formulas (2.5) and (2.6):

$U_{k,n}^{\alpha,\beta,\gamma}(u,v)=\dfrac{1}{u^{\alpha}v^{\beta}(1-u-v)^{\gamma}}\\\
\dfrac{\partial^{n}}{\partial u^{k}\partial
v^{n-k}}\left[u^{k+\alpha}v^{n-k+\beta}(1-u-v)^{n+\gamma}\right]$ (5.5)
$\qquad\qquad\
U_{k,n}^{\alpha,\beta,\gamma}(x,y)=(\alpha+1)_{k}(\beta+1)_{n-k}(1-x-y)^{-\gamma}\\\
F_{2}(-\gamma-n;\alpha+1+k,\beta+1+n-k;\alpha+1,\beta+1;x,y)$

for $0\leq k\leq n$. We can use both formulas to compute the fractional
derivatives. The results will be of course the same, but with (5.5) there
appear some nice results and so we continue with (5.5).

For the integrals inside the brackets in (LABEL:4.9.11) we can use (5.4) and
(5.5) and write in general

$\iint_{V}(.)\,du\,dv=\\\ =\iint_{V}\dfrac{\partial^{k}}{\partial
u^{k}}\dfrac{\partial^{n-k}}{\partial
v^{n-k}}\left[u^{k+\alpha}v^{n-k+\beta}(1-u-v)^{n+\gamma}\right](s-u)^{k-\mu-1}(t-v)^{n-k-\nu-1}\,du\,dv$

After $k$-fold integration by parts with respect to the variable $u$ we get

$\iint_{V}(.)\,du\,dv\left[(-1)^{k}\dfrac{\Gamma(k-\mu)}{\Gamma(-\mu)}\right]^{-1}=\\\
=\iint\nolimits_{V}(s-u)^{-\mu-1}\dfrac{\partial^{n-k}}{\partial
v^{n-k}}\left[u^{k+\alpha}v^{n-k+\beta}(1-u-v)^{n+\gamma}\right](t-v)^{n-k-\nu-1}\,du\,dv$

After $(n-k)$-fold integration by parts with respect to the variable $v$ we
get

$\iint_{V}(.)\,du\,dv\left[(-1)^{n}\dfrac{\Gamma(k-\mu)}{\Gamma(-\mu)}\dfrac{\Gamma(n-k-\nu)}{\Gamma(-\nu)}\right]^{-1}=\\\
=\iint\nolimits_{V}u^{k+\alpha}(s-u)^{-\mu-1}(1-u-v)^{n+\gamma}v^{n-k+\beta}(t-v)^{-\nu-1}\,du\,dv$

Combination of (5.2) and (LABEL:4.9.11) gives

$\displaystyle W^{\mu,\nu}[f](x,y)$
$\displaystyle\left[\dfrac{1}{\Gamma(-\mu)\Gamma(-\nu)}\dfrac{1}{B(k+\alpha+1,n+\gamma+1,n-k+\beta+1)}\dfrac{1}{\delta^{\nu+\mu}}\right]^{-1}=$
$\displaystyle=\int_{0}^{1}\int_{0}^{1-s}f(x+\delta s,y+\delta
t)I_{1}(s,t)\,dt\,ds+\qquad\qquad\text{ \ \ \ \ (I)}$
$\displaystyle+\int_{1}^{\infty}\int_{1}^{\infty}f(x+\delta s,y+\delta
t)I_{3}(s,t)\,dt\,ds+\qquad\qquad\qquad\text{(II)}$
$\displaystyle+\int_{1}^{\infty}\int_{0}^{1}\text{ }f(x+\delta s,y+\delta
t)I_{2}(s,t)\,ds\,dt+\qquad\qquad\qquad\text{(III)}$
$\displaystyle+\int_{1}^{\infty}\int_{0}^{1}\text{ }f(x+\delta s,y+\delta
t)I_{4}(s,t)\,dt\,ds+\qquad\qquad\qquad\text{(IV)}$
$\displaystyle+\int_{0}^{1}\int_{1-s}^{1}\text{ }f(x+\delta s,y+\delta
t)I_{5}(s,t)\,dt\,ds+\qquad\qquad\qquad\text{(V)}$
$\displaystyle\qquad\qquad\qquad+\int_{0}^{1}\int_{1-s}^{1}f(x+\delta
s,y+\delta t)I_{6}(s,t)\,dt\,ds$

Putting $a=k+\alpha+1$, $b=n-k+\beta+1$, $c=-\mu$, $d=-\nu$ and $e=-n-\gamma$
we get for the functions $I_{1..6}(s,t)$

$I_{1}(s,t)=\int_{0}^{t}\left(\int_{0}^{s}u^{a-1}(s-u)^{c-1}(1-u-v)^{-e}du\right)v^{b-1}(t-v)^{d-1}dv$
(5.6)

with region for $s$ and $t$:$\qquad 0<s<1,0<t<1-s$

$I_{2}(s,t)=\int_{0}^{1}\left(\int_{0}^{1-v}u^{a-1}(s-u)^{c-1}(1-u-v)^{-e}du\right)v^{b-1}(t-v)^{d-1}dv$
(5.7)

with region for $s$ and $t$:$\qquad 1<s<\infty,1<t<\infty$

$I_{3}(s,t)=\int_{0}^{s}\left(\int_{0}^{1-u}v^{b-1}(t-v)^{d-1}(1-u-v)^{-e}dv\right)u^{a-1}(s-u)^{c-1}du$
(5.8)

with region for $s$ and $t$:$\qquad 0<s<1<t<\infty$

$I_{4}(s,t)=\int_{0}^{t}\left(\int_{0}^{1-v}u^{a-1}(s-u)^{c-1}(1-u-v)^{-e}du\right)v^{b-1}(t-v)^{d-1}dv$
(5.9)

with region for $s$ and $t$:$\qquad 0<t<1<s<\infty$

$I_{5}(s,t)=\int_{0}^{1-s}\left(\int_{0}^{s}u^{a-1}(s-u)^{c-1}(1-u-v)^{-e}du\right)v^{b-1}(t-v)^{d-1}dv$
(5.10)
$I_{6}(s,t)=\int_{1-s}^{t}\left(\int_{0}^{1-v}u^{a-1}(s-u)^{c-1}(1-u-v)^{-e}du\right)v^{b-1}(t-v)^{d-1}dv$
(5.11)

with region for $s$ and $t$:$\qquad 0<s<1,\ 1-s<t<1$

We will treat the evaluation of the integrals $I_{1,\dots,6}$ in the next
subsections.

###  5.1 The integral $I_{1}(s,t)$

For the computation of the integral $I_{1}(s,t)$ we write the integral (5.6)
as

$I_{1}(s,t)=s^{a+c-1}t^{b+d-1}\int_{0}^{1}\int_{0}^{1}u^{a-1}v^{b-1}(1-u)^{c-1}(1-v)^{d-1}(1-su-
tv)^{-e}\,du\,dv$ (5.12)

The integrals are convergent if
$0<\operatorname{Re}(a)<\operatorname{Re}(a+c)$ and
$0<\operatorname{Re}(b)<\operatorname{Re}(b+d)$. Formula (5.12) is the
integral representation for the $F_{2}$ Appell function [7, 5.8.1.(2)]. We get

$I_{1}(s,t)=B(a,c)B(b,d)s^{a+c-1}t^{b+d-1}F_{2}(e;a,b;a+c,b+d;s,t)$

###  5.2 The integral $I_{2}(s,t)$

For the computation of the integral $I_{2}(s,t)$ we write the integral (5.7)
as

$I_{2}(s,t)=s^{c-1}t^{d-1}\int_{0}^{1}\int_{0}^{1-v}u^{a-1}v^{b-1}(1-u-v)^{-e}\left(1-\dfrac{1}{s}u\right)^{c-1}\left(1-\dfrac{1}{t}v\right)^{d-1}\,du\,dv$
(5.13)

The integrals are convergent if
$\operatorname{Re}(a),\operatorname{Re}(b),\operatorname{Re}(1-e)>0$. Formula
(5.13) is the integral representation for the $F_{3}$ Appell function [7,
5.8.1.(3)]. We get

$I_{2}(s,t)=B(a,b,1-e)s^{c-1}t^{d-1}F_{3}\left(a,b;1-c,1-d;a+b+1-e;\dfrac{1}{s},\dfrac{1}{t}\right)$
(5.14)

###  5.3 The integral $I_{3}(s,t)$

For the computation of the integral $I_{3}(s,t)$ we write the integral (5.8)
as

$I_{3}(s,t)=s^{a+c-1}t^{d-1}\\\
\int_{0}^{1}\int_{0}^{1}u^{a-1}v^{b-1}(1-u)^{c-1}(1-v)^{-e}(1-su)^{b-e}\left(1-\dfrac{1}{t}v+\dfrac{s}{t}uv\right)^{d-1}\,du\,dv$
(5.15)

Formula (5.15) is an integral representation of the $H_{2}$ function [5,
(3.4)]. We get

$I_{3}(s,t)=B(a,c)B(b,1-e)s^{a+c-1}t^{d-1}H_{2}\left(e-b,a,1-d,b;a+c;s,-\dfrac{1}{t}\right)$
(5.16)

###  5.4 The integral $I_{4}(s,t)$

For $I_{4}(s,t)$ we look at $I_{2}(s,t)$. Interchanging in $I_{2}(s,t)$ the
parameters $a$ and $b$, $c$ and $d$ and the variables $s$ and $t$ we get
$I_{4}(s,t)$. So from (5.16) we get directly

$I_{4}(s,t)=B(a,1-e)B(b,d)t^{b+d-1}s^{c-1}H_{2}\left(e-a,b,1-c,a;b+d;t,-\dfrac{1}{s}\right)$

###  5.5 The integral $I_{5}(s,t)$

For $I_{5}(s,t)$ we write the inner integral of (5.10) as

$\int_{0}^{s}u^{a-1}(s-u)^{c-1}(1-v-u)^{-e}du=B(a,c)s^{a+c-1}(1-v)^{-e}\,\mbox{}_{2}F_{1}\\!\left(\genfrac{}{}{0.0pt}{}{a,e}{a+c};\dfrac{s}{1-v}\right)$

with convergence conditions
$\operatorname{Re}(a),\operatorname{Re}(c),\operatorname{Re}(c-e)>0$ and
$\dfrac{s}{1-v}\leq 1$. Substitution in (5.10) and replacing $v$ by $(1-s)v$
gives

$I_{5}(s,t)=B(a,c)s^{a+c-1}t^{d-1}(1-s)^{b}\\\
\times\int_{0}^{1}v^{b-1}\big{(}1-(1-s)v\big{)}^{-e}\left(1-\dfrac{1-s}{t}v\right)^{d-1}\,\mbox{}_{2}F_{1}\\!\left(\genfrac{}{}{0.0pt}{}{a,e}{a+c};\dfrac{s}{1-(1-s)v}\right)dv$
(5.17)

Because the computation of the integral is very long we only give the result
from [4, section 6.5.5.]

$\displaystyle I_{5}(s,t)$
$\displaystyle=B(c,e-c)B(b,c-e+1)s^{a-1}t^{d-1}(1-s)^{b+c-e}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\times
F_{3}\left(\begin{array}[]{c}1-a,b;c,1-d\\\
b+c-e+1\end{array};\dfrac{s-1}{s},\dfrac{1-s}{t}\right)+$ (5.20)
$\displaystyle+B(a,c-e)B(b,d)s^{a+c-1}t^{b+d-1}F_{P}(e,b,a,b+d,a+c;t,s)-$
$\displaystyle-\dfrac{B(a,c-e)}{d}s^{a+c-e-1}t^{b-1}(s+t-1)^{d}$
$\displaystyle\qquad\times\widetilde{F_{3}}\left(\begin{array}[]{c}1,1-b;e-a-c+1,d\\\
d+1\end{array};\begin{array}[]{c}e\\\
e-c+1\end{array};\dfrac{s+t-1}{s},\dfrac{s+t-1}{t}\right)$ (5.25)

with convergence conditions for the $F_{P}$ function: $\left|t\right|<1$,
$\left|s-1\right|<1$. The convergence conditions of the $F_{3}$ function are:
$\left|\dfrac{s-1}{s}\right|<1$ and $\left|\dfrac{1-s}{t}\right|<1$. This
gives a region of convergence $\dfrac{1}{2}<s<1$. But after analytical
continuation (see (B.1)) this function is well-defined for $\dfrac{s-1}{s}<1$
and $\dfrac{1-s}{t}<1$. For the extended $F_{3}$ function we have with
analytical continuation convergence for the same region as for the $F_{3}$
function (see appendix B.2).This includes the $(s,t)$ domain on which the
integral $I_{5}(s,t)$ is defined. For the parameters we have the convergence
conditions:
$0<\operatorname{Re}(a),\operatorname{Re}(b),\operatorname{Re}(c),\operatorname{Re}(d),\operatorname{Re}(e)<\operatorname{Re}(c)+1$.

###  5.6 The integral $I_{6}(s,t)$

For $I_{6}(s,t)$ we write the inner integral of (5.11) as

$\int_{0}^{1-v}u^{a-1}(s-u)^{c-1}(1-v-u)^{-e}du=B(a,1-e)s^{c-1}(1-v)^{a-e}\,\mbox{}_{2}F_{1}\\!\left(\genfrac{}{}{0.0pt}{}{a,1-c}{a-e+1};\dfrac{1-v}{s}\right)$

with convergence conditions $\operatorname{Re}(a),\operatorname{Re}(1-e)>0$
and $\dfrac{1-v}{s}<1$.

Because the computation of the integral is also very long we only give the
result from [4, section 6.5.6.]

$\displaystyle I_{6}(s,t)$
$\displaystyle=\dfrac{B(a,c-e)}{d}s^{a+c-e-1}t^{b-1}(s+t-1)^{d}$
$\displaystyle\qquad\qquad\times\widetilde{F_{3}}\left(\begin{array}[]{c}1,1-b;e-a-c+1,d\\\
d+1\end{array};\begin{array}[]{c}e\\\
e-c+1\end{array};\dfrac{s+t-1}{s},\dfrac{s+t-1}{t}\right)+$ (5.30)
$\displaystyle+B(1-e,e-c)B(d,1-e+c)s^{a-1}t^{b-1}(s+t-1)^{c+d-e}$
$\displaystyle\qquad\qquad\times F_{3}\left(\begin{array}[]{c}1-a,1-b;c,d\\\
c+d-e+1\end{array};\dfrac{s+t-1}{s},\dfrac{s+t-1}{t}\right)$ (5.33)

with convergence conditions:
$\operatorname{Re}(a),\operatorname{Re}(d),\operatorname{Re}(1-e),\operatorname{Re}(c-e+1)>0$.

###  5.7 The summation of the integrals $I_{5}(s,t)$ and $I_{6}(s,t)$

When adding the functions $I_{5}(s,t)$ (5.25) and $I_{6}(s,t)$ (5.33) the
extended $F_{3}$ functions will be eliminated. We get

$\displaystyle I_{5}(s,t)+I_{6}(s,t)$
$\displaystyle=B(a,c-e)B(b,d)s^{a+c-1}t^{b+d-1}F_{P}(e,b,a,b+d,a+c;t,s)+$
$\displaystyle+B(c,e-c)B(b,c-e+1)s^{a-1}(1-s)^{b+c-e}t^{d-1}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad
F_{3}\left(\begin{array}[]{c}1-a,b;c,1-d\\\
b+c-e+1\end{array};\dfrac{s-1}{s},\dfrac{1-s}{t}\right)+$ (5.36)
$\displaystyle+B(1-e,e-c)B(d,1-e+c)s^{a-1}t^{b-1}(s+t-1)^{c+d-e}$
$\displaystyle\qquad\qquad\qquad\qquad
F_{3}\left(\begin{array}[]{c}1-a,1-b;c,d\\\
c+d-e+1\end{array};\dfrac{s+t-1}{s},\dfrac{s+t-1}{t}\right)$ (5.39)

The power series of the functions $F_{P}$ (B.7) and $F_{3}$ (B.1) are
convergent in the region $0<s<1,1-s<t<1$. For the first $F_{3}$ function we
use the analytical continuation over the whole $(s,t)$ region. For the
parameters we get the convergence conditions:
$0<~{}\operatorname{Re}(a),\operatorname{Re}(b),\operatorname{Re}(c),\operatorname{Re}(d)$
and $\operatorname{Re}(e)<1$. From [13, p 441] we get

$\displaystyle F_{3}\left(\begin{array}[]{c}a_{0},b_{1};b_{2},c_{1}\\\
c_{2}\end{array};x_{1},x_{2}\right)$
$\displaystyle=\dfrac{\Gamma(c_{2})\Gamma(c_{2}-b_{1}-c_{1})}{\Gamma(c_{2}-b_{1})\Gamma(c_{2}-c_{1})}(x_{1})^{-a_{0}}(x_{2})^{-b_{1}}$
(5.42) $\displaystyle
F_{Q}\left(b_{2}+c_{1}-c_{2}+1,b_{2},c_{1};1-a_{0}+b_{2},c_{1}-b_{1}+1;\dfrac{1}{x_{1}},\dfrac{1}{x_{2}}\right)+$
$\displaystyle+\dfrac{\Gamma(c_{2})\Gamma(b_{1}+c_{1}-c_{2})}{\Gamma(b_{1})\Gamma(c_{1})}(x_{2})^{1-c_{2}}(1-x_{2})^{c_{2}-b_{1}-c_{1}}$
$\displaystyle\qquad\quad
F_{3}\left(\begin{array}[]{c}a_{0},1-b_{1};b_{2},1-c_{1}\\\
1-b_{1}-c_{1}+c_{2}\end{array};\dfrac{x_{1}(x_{2}-1)}{x_{2}},1-x_{2}\right)$
(5.45)

After application of (5.45) to the second $F_{3}$ function of (5.39) we get

$\displaystyle I_{5}(s,t)+I_{6}(s,t)$
$\displaystyle=B(a,c-e)B(b,d)s^{a+c-1}t^{b+d-1}F_{P}(e,b,a,b+d,a+c;t,s)+$
$\displaystyle+B(1-e,e-c)B(d,c-e+b)(s+t-1)^{a+b+c+d-e-2}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad
F_{Q}\left(e,c,d;a+c,b+d;\dfrac{s}{s+t-1},\dfrac{t}{s+t-1}\right)+$
$\displaystyle+\Gamma(e-c)\Gamma(c-e+1)\left[\dfrac{\Gamma(b)\Gamma(c)}{\Gamma(e)\Gamma(c+b-e+1)}+\dfrac{\Gamma(1-e)\Gamma(e-b-c)}{\Gamma(1-c)\Gamma(1-b)}\right]$
$\displaystyle\qquad\qquad\qquad\qquad
s^{a-1}(1-s)^{b+c-e}t^{d-1}F_{3}\left(\begin{array}[]{c}1-a,b;c,1-d\\\
b+c-e+1\end{array};\dfrac{s-1}{s},\dfrac{1-s}{t}\right)$

The function $F_{Q}$ is defined in (B.8). After some simplification of the
terms within the square brackets we get

$\displaystyle I_{5}(s,t)+I_{6}(s,t)$
$\displaystyle=B(a,c-e)B(b,d)s^{a+c-1}t^{b+d-1}F_{P}(e,b,a,b+d,a+c;t,s)+$
$\displaystyle+B(1-e,e-c)B(d,c-e+b)s^{a+c-1}t^{b+d-1}(s+t-1)^{-e}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad
F_{Q}\left(e,c,d;a+c,b+d;\dfrac{s}{s+t-1},\dfrac{t}{s+t-1}\right)+$
$\displaystyle+B(e-b-c,c)B(1-e,b)s^{a-1}\left(1-s\right)^{b+c-e}t^{d-1}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad
F_{3}\left(\begin{array}[]{c}1-a,b;c,1-d\\\
b+c-e+1\end{array};\dfrac{s-1}{s},\dfrac{1-s}{t}\right)$

The power series of the $F_{P}$ and the $F_{Q}$ functions are convergent in
the region $\\{0<s<1\wedge 1-s<t<1\\}$. For the $F_{3}$ function we use an
analytical continuation so we get the convergence in the right region. For the
parameters we get the convergence conditions:
$0<\operatorname{Re}(a,b,c,d,e)<1$.

###  5.8 Summary of the integrals $I_{1..6}(s,t)$

We give a summary of the functions $I_{1..9}(s,t)$

$I_{1}(s,t)=B(a,c)B(b,d)s^{a+c-1}t^{b+d-1}F_{2}(e;a,b;a+c,b+d;s,t)$ (5.46)

with $0<s<1,0<t<1-s$ and $0<\operatorname{Re}(a)<\operatorname{Re}(a+c)$,
$0<\operatorname{Re}(b)<\operatorname{Re}(b+d)$.

$I_{2}(s,t)=B(a,b,1-e)s^{c-1}t^{d-1}F_{3}\left(a,b;1-c,1-d;a+b+1-e;\dfrac{1}{s},\dfrac{1}{t}\right)$

with $1<s<\infty,1<t<\infty$ and
$\operatorname{Re}(a),\operatorname{Re}(b),\operatorname{Re}(1-e)>0$.

$I_{3}(s,t)=B(a,c)B(b,1-e)s^{a+c-1}t^{d-1}H_{2}\left(e-b,a,1-d,b;a+c;s,-\dfrac{1}{t}\right)$
(5.47)

with $0<s<1<t<\infty$ and
$\operatorname{Re}(a),\operatorname{Re}(b)>0,\operatorname{Re}(c)>0$.

$I_{4}(s,t)=B(a,1-e)B(b,d)t^{b+d-1}s^{c-1}H_{2}\left(e-a,b,1-c,a;b+d;t,-\dfrac{1}{s}\right)$
(5.48)

with $0<t<1<s<\infty$ and
$\operatorname{Re}(a),\operatorname{Re}(b),\operatorname{Re}(d)>0$.

$\displaystyle I_{5}(s,t)+I_{6}(s,t)$
$\displaystyle=B(a,c-e)B(b,d)s^{a+c-1}t^{b+d-1}F_{P}(e,b,a,b+d,a+c;t,s)+$
$\displaystyle+B(c,e-c)B(b,c-e+1)s^{a-1}(1-s)^{b+c-e}t^{d-1}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad
F_{3}\left(\begin{array}[]{c}1-a,b;c,1-d\\\
c+b-e+1\end{array};\dfrac{s-1}{s},\dfrac{1-s}{t}\right)+$ (5.51)
$\displaystyle+B(1-e,e-c)B(d,1-e+c)t^{b-1}s^{a-1}(s+t-1)^{c+d-e}$
$\displaystyle\qquad\qquad\qquad\qquad
F_{3}\left(\begin{array}[]{c}1-a,1-b;c,d\\\
c+d-e+1\end{array};\dfrac{s+t-1}{s},\dfrac{s+t-1}{t}\right)$ (5.54)

with $0<s<1,\ 1-s<t<1$ and
$0<\operatorname{Re}(a),\operatorname{Re}(c),\operatorname{Re}(d),\operatorname{Re}(e)<1$.
The functions $F_{2}$, $F_{3}$, $H_{2}$ and $F_{P}$ are all five-parametric
double-hypergeometric functions which are solutions of a system of two partial
differential equations of second order [14]. This system is associated with
the Appell function $F_{2}(a,b_{1},b_{2},c_{1},c_{2};x,y)$ [7, 5.9(10)]. Note
that the orthogonal polynomials (2.6) and (2.8) used for the fractional
orthogonal derivative are also expressed in terms of $F_{2}$ functions.

###  5.9 Remarks

######  Remark 5.1.

It is well-known [7, 5.9(10)] that the
$F_{2}(a_{0};b_{1},b_{2};c_{1},c_{2};x,y)$ function is a solution of the
following system of partial differential equations

$x(1-x)r-xys+\Big{[}c_{1}-(a_{0}+b_{1}+1)x\Big{]}p-b_{1}yq-
a_{0}b_{1}F_{2}(x,y)=0$ (5.55)
$y(1-y)t-xys+\Big{[}c_{2}-(a_{0}+b_{2}+1)y\Big{]}q-b_{2}xp-
a_{0}b_{2}F_{2}(x,y)=0$ (5.56)

with

$p=\dfrac{\partial F_{2}(x,y)}{\partial x}\ \ \ q=\dfrac{\partial
F_{2}(x,y)}{\partial y}\ \ \ r=\dfrac{\partial^{2}F_{2}(x,y)}{\partial x^{2}}\
\ \ s=\dfrac{\partial^{2}F_{2}(x,y)}{\partial x\partial y}\ \ \
t=\dfrac{\partial^{2}F_{2}(x,y)}{\partial y^{2}}$

Olsson treats all 36 explicit solutions of (5.55) and (5.56) [14]. Using the
substitutions $a_{0}=e-a-b-c-d+2,\ b_{1}=1-c,\ b_{2}=1-d,\ c_{1}=2-a-c,\
c_{2}=2-b-d$ and $x\rightarrow s,\ y\rightarrow t$ he gives among others the
following solutions of (5.55) and (5.56)

$s^{a+c-1}t^{b+d-1}F_{2}\left(\begin{array}[]{c}e;a,b\\\
a+c,b+d\end{array};s,t\right)$ (5.57)
$s^{a+c-1}F_{2}\left(\begin{array}[]{c}e-b-d+1;a,1-d\\\
a+c,2-b-d\end{array};s,t\right)$ (5.58)
$s^{c-1}t^{d-1}F_{3}\left(\begin{array}[]{c}a,b;1-c,1-d\\\
a+b+1-e\end{array};\dfrac{1}{s},\dfrac{1}{t}\right)$ (5.59)
$t^{b+d-1}F_{P}(e-c-a+1,b,1-c,b+d,2-a-c;t,s)$ (5.60)
$s^{a-1}t^{d-1}(1-s)^{b+c-e}F_{3}\left(\begin{array}[]{c}1-a,b;c,1-d\\\
b+c-e+1\end{array};\dfrac{s-1}{s},\dfrac{1-s}{t}\right)$ (5.61)
$s^{a-1}t^{b-1}(1-s-t)^{c+d-e}F_{3}\left(\begin{array}[]{c}1-a,1-b;c,d\\\
c+d-e+1\end{array};\dfrac{s+t-1}{s},\dfrac{s+t-1}{t}\right)$ (5.62)
$s^{a+c-1}t^{d-1}H_{2}\left(e-b,a,1-d,b;a+c;s,-\dfrac{1}{t}\right)$ (5.63)

For the integral $I_{1}(s,t)$ we see that the function
$s^{a+c-1}t^{b+d-1}F_{2}\left(\begin{array}[]{c}e;a,b\\\
a+c,b+d\end{array};s,t\right)$ equals the solution (5.57).

For the integral $I_{2}(s,t)$ we see that the function
$s^{c-1}t^{d-1}F_{3}\left(a,b;1-c,1-d;a+b+1-e;\dfrac{1}{s},\dfrac{1}{t}\right)$
equals the solution (5.59).

For the integral $I_{3}(s,t)$ we see that the function
$s^{a+c-1}t^{d-1}H_{2}\left(e-b,a,1-d,b;a+c;s,-\dfrac{1}{t}\right)$ equals the
solution (5.63).

To compute $I_{4}(s,t)$ we can interchange in $I_{3}(s,t)$ the parameters $a$
and $b$, $c$ and $d$ and the variables $s$ and $t$. Then we see that (5.48) is
a solution of (5.55) and (5.56).

For the integral $I_{5}(s,t)+I_{6}(s,t)$ we have the functions

$L_{1}(s,t)=s^{a+c-1}t^{b+d-1}F_{P}(e,b,a,b+d,a+c;t,s)$ (5.64)

$L_{2}(s,t)=s^{a-1}(1-s)^{b+c-e}t^{d-1}F_{3}\left(\begin{array}[]{c}1-a,b;c,1-d\\\
b+c-e+1\end{array};\dfrac{s-1}{s},\dfrac{1-s}{t}\right)$

$L_{3}\left(s,t\right)=s^{a-1}t^{b-1}(s+t-1)^{c+d-e}F_{3}\left(\begin{array}[]{c}1-a,1-b;c,d\\\
c+d-e+1\end{array};\dfrac{s+t-1}{s},\dfrac{s+t-1}{t}\right)$

For the function $L_{1}(s,t)$ we use the solution [14, (26)] and interchange
$b_{1}\longleftrightarrow b_{2},c_{1}\longleftrightarrow c_{2}$ and
$x_{1}\longleftrightarrow x_{2}$. Then the solution becomes

$\left(x_{2}\right)^{1-c_{2}}F_{P}(a_{0}-c_{2}+1,b_{1},b_{2}c_{2}+1,c_{1},2-c_{2};x_{2},x_{1})$

Application to (5.64) gives as result

$t^{b+d-1}F_{P}(e-c-a+1,b,1-c,b+d,2-a-c;t,s)$

This function is equal (5.60). So $L_{1}(s,t)$ is a solution of (5.57) and
(5.58).

The function $L_{2}(s,t)$ equals the solution (5.61).

The function $L_{3}(s,t)$ equals the solution (5.62).

######  Remark 5.2.

It is clear that the inner double integrals in (LABEL:4.9.11) are in fact
double integrals transformations which transform the solutions of the $F_{2}$
equations (5.55) and (5.56) into a solution of the $F_{2}$ equations with
other parameters.

######  Remark 5.3.

When looking at Figure 2 case Va we can consider the integral when the
arguments are on the boundaries. Then we distinguish four possibilities.

I. $s=1$.

For $s=1$ it can be shown that $I_{5}(1,t)+I_{6}(1,t)=I_{4}(1,t)$.
II. $t=1$.

For $t=1$ it can be shown that $I_{5}(s,1)+I_{6}(s,1)=I_{3}(s,1)$.
III. $s=t=1$.

From (5.17) it follows $I_{5}(1,t)=0$ and $I_{6}(1,t)=I_{4}(1,t)$. For $s=t=1$
it can be shown that $I_{6}(1,1)=I_{2}(1,1)=I_{3}(1,1)$.
IV. $s+t=1$.

For $s+t=1$ it can be shown that $I_{5}(s,1-s)+I_{6}(s,1-s)=I_{1}(s,1-s)$.

######  Remark 5.4.

When looking at case V in Figure 2 we see directly that when interchanging
$a\leftrightarrow b$, $c\leftrightarrow d$, $s\leftrightarrow t$ and the
integration variables $u\leftrightarrow v$ we get

$I_{5}(a,b,c,d,e;s,t)+I_{6}(a,b,c,d,e;s,t)=I_{5}(b,a,d,c;e;t,s)+I_{6}(b,a,d,c;e;t,s)$
(5.65)

From the symmetry there follows that $I_{5}(s,t)+I_{6}(s,t)$ can be written as
half the sum of the two sides of (5.65).

Looking at equation (5.54) we see that only the second $F_{3}$ function in the
right-hand side is symmetrical in the parameters and the arguments. Now write
(5.65) as

$\displaystyle I_{5}(a,b,c,d,e;s,t)+I_{6}(a,b,c,d,e;s,t)$
$\displaystyle=\lambda\big{(}I_{5}(a,b,c,d,e;s,t)+I_{6}(a,b,c,d,e;s,t)\big{)}+$
$\displaystyle+(1-\lambda)\big{(}I_{5}(b,a,d,c;e;t,s)+I_{6}(b,a,d,c;e;t,s)\big{)}$

Now we can choose $\lambda$ so that the second $F_{3}$ functions will be
eliminated. There results

$\displaystyle I_{5}(s,t)+I_{6}(s,t)$
$\displaystyle=A_{1}\;s^{a+c-1}t^{b+d-1}F_{P}(e,b,a,b+d,a+c;t,s)+$
$\displaystyle+A_{2}\;s^{a+c-1}t^{b+d-1}F_{P}(e,a,b,a+c,b+d;s,t)+$
$\displaystyle+A_{3}\;t^{b-1}(1-t)^{a+d-e}s^{c-1}F_{3}\left(\begin{array}[]{c}d,a;1-b,1-c\\\
a+d-e+1\end{array};\dfrac{t-1}{t},\dfrac{1-t}{s}\right)+$
$\displaystyle+A_{4}\;s^{a-1}(1-s)^{b+c-e}t^{d-1}F_{3}\left(\begin{array}[]{c}c,b;1-a,1-d\\\
b+c-e+1\end{array};\dfrac{s-1}{s},\dfrac{1-s}{t}\right)$

with

$A_{1}=\dfrac{\Gamma(a)\Gamma(b)\Gamma(e)\Gamma(c-d)\Gamma(d-c+1)\Gamma(1-e)}{\Gamma(a+c-e)\Gamma(b+d)\Gamma(1-d)\Gamma(e-c+1)}$
$A_{2}=\dfrac{\Gamma(a)\Gamma(b)\Gamma(e)\Gamma(d-c)\Gamma(c-d+1)\Gamma(1-e)}{\Gamma(b+d-e)\Gamma(a+c)\Gamma(1-c)\Gamma(e-d+1)}$
$A_{3}=\dfrac{\Gamma(a)\Gamma(d)\Gamma(c-d)\Gamma(d-c+1)\Gamma(1-e)}{\Gamma(c)\Gamma(1-c)\Gamma(a+d-e+1)}\qquad$
$A_{4}=\dfrac{\Gamma(b)\Gamma(c)\Gamma(d-c)\Gamma(c-d+1)\Gamma(1-e)}{\Gamma(d)\Gamma(1-d)\Gamma(c+b-e+1)}\qquad$

Application of Olsson [14, (26)] gives

$\displaystyle I_{5}(s,t)+I_{6}(s,t)$
$\displaystyle=A_{1}\;s^{a+c-1}t^{b+d-1}(1-t)^{-e}F_{P}\left(e,d,a,b+d,a+c;\dfrac{t}{t-1},\dfrac{s}{1-t}\right)+$
$\displaystyle+A_{3}\;s^{c-1}t^{b-1}(1-t)^{a+d-e}F_{3}\left(\begin{array}[]{c}d,a;1-b,1-c\\\
a+d-e+1\end{array};\dfrac{t-1}{t},\dfrac{1-t}{s}\right)+$
$\displaystyle+A_{2}\;s^{a+c-1}t^{b+d-1}(1-s)^{-e}F_{P}\left(e,c,b,a+c,b+d;\dfrac{s}{s-1},\dfrac{t}{1-s}\right)+$
$\displaystyle+A_{4}\;s^{a-1}t^{d-1}(1-s)^{b+c-e}F_{3}\left(\begin{array}[]{c}c,b;1-a,1-d\\\
b+c-e+1\end{array};\dfrac{s-1}{s},\dfrac{1-s}{t}\right)$

Note that the arguments of the $F_{P}$ functions are the inverses of the
arguments of the $F_{3}$ functions.

Appendices: Hypergeometric functions of two variables

When deriving formulas for the two-dimensional fractional orthogonal
derivatives we encounter various hypergeometric functions of two variables,
some of which are beyond the familiar Appell and Horn hypergeometric
functions. Here we give an overview of the functions we need in this thesis.

##  Appendix A The Pochhammer symbol

We have already frequently used the Pochhammer symbol, which we defined for
$a\in\mathbb{C}$ by

$(a)_{i}:=a(a+1)\ldots(a+i-1)=\frac{\Gamma(a+i)}{\Gamma(a)}\quad(i=0,1,\ldots),$

and more generally by

$(a)_{b}:=\frac{\Gamma(a+b)}{\Gamma(a)}\quad(b\in\mathbb{C}).$

We will use this second definition in particular for $b\in\mathbb{Z}$. Then

$(a)_{-i}=\frac{(-1)^{i}}{(1-a)_{i}}\quad(i\in\mathbb{Z}).$

Further properties for $i,j\in\mathbb{Z}$ are:

$\Gamma(a+i)=\Gamma(a)(a)_{i}$
$\Gamma(a-i)=(-1)^{i}\dfrac{\Gamma(a)}{(1-a)_{i}}$
$(a+i)_{j}=\dfrac{(a)_{i+j}}{(a)_{i}}=\dfrac{(a+j)_{i}(a)_{j}}{(a)_{i}}$
$(a-i)_{j}=(-1)^{i}(a)_{j-i}(1-a)_{i}=(-1)^{j}\dfrac{(1-a)_{i}}{(1-a)_{i-j}}=(a)_{j}\dfrac{(1-a)_{i}}{(1-a-j)_{i}}$

There are a lot more. See for example [15, Appendix I].

##  Appendix B Definitions

For the basics of hypergeometric functions of two variables see for instance
[3, Ch. 9], [7, sections 5.7-5.12], [8, Ch. 1], [15, Ch. 8].

Appell [1], [2] defines four functions $F_{1},F_{2},F_{3},F_{4}$ in two
variables which are generalizations of the ${}_{2}F_{1}$ hypergeometric
function. We use only the first three functions [7, 5.7].

$F_{1}(a;b_{1},b_{2};c;x,y)=\sum\limits_{i=0}^{\infty}\sum\limits_{j=0}^{\infty}\dfrac{(a)_{i+j}(b_{1})_{i}(b_{2})_{j}}{(c)_{i+j}}\dfrac{1}{i!\,j!}x^{i}y^{j}$

with convergence region:
$\left\\{\left|x\right|<1\wedge\left|y\right|<1\right\\}$.

$F_{2}(a;b_{1},b_{2};c_{1},c_{2};x,y)=\sum\limits_{i=0}^{\infty}\sum\limits_{j=0}^{\infty}\dfrac{(a)_{i+j}(b_{1})_{i}(b_{2})_{j}}{(c_{1})_{i}(c_{2})_{j}}\dfrac{1}{i!\,j!}x^{i}y^{j}$

This power series has as convergence region: $\\{|x|+|y|<1\\}$. From the known
transformations of the $F_{2}$ function [7, 5.11. (6),(7),(8)] and from the
union of the four partially overlapping regions of convergence (or from the
integral representation of the $F_{2}$ function [7, 5.8(2)]) we get a region
of unique analytic convergence $\left\\{x<1\wedge y<1\wedge x+y<1\right\\}$.

$F_{3}(a_{1},a_{2};b_{1},b_{2};c;x,y)=\sum\limits_{i=0}^{\infty}\sum\limits_{j=0}^{\infty}\dfrac{(a_{1})_{i}(a_{2})_{j}(b_{1})_{i}(b_{2})_{j}}{(c)_{i+j}}\dfrac{1}{i!\,j!}x^{i}y^{j}$
(B.1)

with convergence region: $\left\\{|x|<1\wedge|y|<1\right\\}$. With the
integral representation of the $F_{3}$ function [7, 5.8(3)]222Note that the
power of the factor $(1-u-v)$ should be $\gamma-\beta-\beta^{\prime}.$ we get
a region of unique analytical continuation of $\\{x<1,\,y<1\\}$.

Another integral representation is

$F_{3}(a_{1},a_{2};b_{1},b_{2};c;x,y)=\dfrac{\Gamma(c)}{\Gamma(a_{1})\Gamma(a_{2})\Gamma(c-a_{1}-a_{2})}\\\
\times\int_{0}^{1}\int_{0}^{1}u^{a_{1}-1}v^{a_{2}-1}(1-u)^{c-a_{2}-a_{1}-1}(1-v)^{c-a_{2}-1}(1-y\,v)^{-b_{2}}(1-x\,u+x\,u\,v)^{-b_{1}}dudv$
(B.2)

with $\operatorname{Re}(a_{1},a_{2},c-a_{2}-a_{1},c-a_{2})>0$.

We define the extended $F_{3}$ function as

$\widetilde{F_{3}}\left(\begin{array}[]{c}a_{1},a_{2},b_{1},b_{2}\\\
c\end{array};\begin{array}[]{c}d_{1}\\\
d_{2}\end{array};x,y\right)=\sum\limits_{i=0}^{\infty}\sum\limits_{j=0}^{\infty}\dfrac{(a_{1})_{i}(a_{2})_{j}(b_{1})_{i}(b_{2})_{j}}{(c)_{i+j}}\dfrac{(d_{1})_{i}}{(d_{2})_{i}}\dfrac{1}{i!\,j!}x^{i}y^{j}$
(B.3)

This extended $F_{3}$ function is of order three (see below) and can be
written as a triple integral

$\widetilde{F_{3}}\left(\begin{array}[]{c}a_{1},a_{2},b_{1},b_{2}\\\
c\end{array};\begin{array}[]{c}d_{1}\\\
d_{2}\end{array};x,y\right)=\dfrac{\Gamma(c)}{\Gamma(a_{1})\Gamma(a_{2})\Gamma(c-a_{1}-a_{2})}\dfrac{\Gamma(d_{2})}{\Gamma(d_{1})\Gamma(d_{2}-d_{1})}\\\
\qquad\qquad\times\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}u^{a_{1}-1}v^{a_{2}-1}w^{d_{1}-1}(1-u)^{c-a_{1}-1}(1-v)^{c-a_{1}-a_{2}-1}(1-w)^{d_{2}-d_{1}-1}\\\
\times(1-yv+yuv)^{-b_{2}}(1-x\,u\,w)^{-b_{1}}dudvdw$

Then we have the convergence conditions for the parameters:
$\operatorname{Re}(c)>\operatorname{Re}(a_{1})>0,\operatorname{Re}(c-~{}a_{1})>\operatorname{Re}(a_{2})>0,\operatorname{Re}(d_{2})>\operatorname{Re}(d_{1})$.
For the arguments we have the analytical continuation of convergence:
$\\{x<1\wedge y<1\\}$.

Horn [12] more generally calls a double power series

$\sum\limits_{i,j=0}^{\infty}A\left(i,j\right)x^{i}y^{j}$

hypergeometric if the two quotients $\dfrac{A(i+1,j)}{A(i,j)}$ and
$\dfrac{A(i,j+1)}{A(i,j)}$ are rational functions of $i$ and $j$. Write the
two rational functions as quotients of polynomials in $i$ and $j$ without
common factors.

$\dfrac{A(i+1,j)}{A(i,j)}=\dfrac{F(i,j)}{F^{\prime}(i,j)}\qquad$and
$\dfrac{A(i,j+1)}{A(i,j)}=\dfrac{G(i,j)}{G^{\prime}(i,j)}$

In addition it is assumed that $F^{\prime}(i,j)$ contains the factor $i+1$ and
$G^{\prime}(i,j)$ the factor $j+1$. Then the highest degree in $i,j$ of the
four polynomials $F,F^{\prime},G,G^{\prime}$ is called the $order$ of the
hypergeometric series. Horn also gives a rule (see [7, Section 5.7.2]) how to
determine the region of convergence of the power series from the two above
quotients. Then Horn classifies the hypergeometric functions of order two. He
gives a list of 34 such functions. The first four items in the list are the
Appell functions $F_{1},\ F_{2},\ F_{3}$ and $F_{4}$.

Number nine in the list of Horn is the $H_{2}$ function.

$H_{2}(a,b_{1},b_{2},c_{1},c_{2};x,y)=\sum\limits_{i=0}^{\infty}\sum\limits_{j=0}^{\infty}\dfrac{(a)_{i-j}(b_{1})_{i}(b_{2})_{j}(c_{1})_{j}}{(c_{2})_{i}}\dfrac{1}{i!\,j!}x^{i}y^{j}$
(B.4)

with convergence condition (determined by Horn’s rule):

$\left\\{|x|<1,|y|<(|x|+1)^{-1}\right\\}$ (B.5)

In [5, Section 3] Diekema and Koornwinder showed that

$\\{(x<0\wedge(x-1)y<1)\vee(0\leq x<1\wedge y>-1)\\}$

is a region in $\mathbb{R}^{2}$ where $H_{2}$ has unique analytical
continuation.

Figure 3: Convergence region of the $H_{2}$ function

Figure 4: Region of analytical continuation of $H_{2}$

In [7, section 5.7] for each function in Horn’s list a system of two pde’s is
given which has this function as a solution. This follows a method already
indicated by Horn in [11], [12].

A comprehensive list of the solutions of the system of pde’s (5.55),(5.56) for
the $F_{2}$ function was given by Olsson, [14, p.1289, table I]. Apart from
$F_{2},\ F_{3}$ and $H_{2}$ there are solutions expressed in terms of
functions $F_{P},\ F_{Q}$ and $F_{R}$ which do not occur in Horn’s list [13].
Here $F_{P}$ and $F_{Q}$ have order three, while $F_{R}$ is not even of
hypergeometric type. In his paper he also mentioned an $F_{PR}$ function [14,
(57)].

Olsson defines his $F_{P}$ function as:

$F_{P}(a,b_{1},b_{2},c_{1},c_{2};x,y)=y^{-a}\sum\limits_{i=0}^{\infty}\sum\limits_{j=0}^{\infty}\dfrac{(a)_{i+j}(a-c_{2}+1)_{i+j}(b_{1})_{i}}{(a+b_{2}-c_{2}+1)_{i+j}(c_{1})_{i}}\dfrac{1}{i!\,j!}\left(\dfrac{x}{y}\right)^{i}\left(\dfrac{y-1}{y}\right)^{j}$
(B.6)

with convergence region: $\\{|xy^{-1}|+|1-y^{-1}|<1\\}$. In [13] (in a
different notation) he gives

$F_{P}(a,b_{1},b_{2},c_{1},c_{2};x,y)=y^{-a}\sum\limits_{i=0}^{\infty}\sum\limits_{j=0}^{\infty}\dfrac{(a)_{i+j}(a-c_{2}+1)_{i}(b_{2})_{j}}{(a+b_{2}-c_{2}+1)_{i+j}}\dfrac{(b_{1})_{i}}{(c_{1})_{i}}\dfrac{1}{i!\,j!}x^{i}(1-y)^{j}$
(B.7)

with convergence region $\\{|x|<1,|y-1|<1\\}$. In [5, Theorem 4.1] we showed
that $\\{x<1\wedge y>0\\}$ is a region in $\mathbb{R}^{2}$ where $F_{P}$ has
unique analytical continuation.

Olsson [13], [14] defined his $F_{Q}$ function as

$F_{Q}(a,b_{1},b_{2},c_{1},c_{2};x,y)=x^{-b_{1}}y^{b_{1}-a}\\\
\times\sum\limits_{i=0}^{\infty}\sum\limits_{j=0}^{\infty}\dfrac{(b_{1}+c_{2}-b_{2}-a)_{i-j}(b_{1})_{i}(b_{1}-c_{1}+1)_{i}}{(b_{1}-a+1)_{i-j}(b_{1}+c_{2}-a)_{i-j}}\dfrac{1}{i!\,j!}\left(\dfrac{y}{x}\right)^{i}\left(\dfrac{1-y}{y}\right)^{j}$
(B.8)

where the convergence conditions (by Horn’s rule) are the same as for the
function

$H_{2}\left(1-\dfrac{1}{y},-\dfrac{y}{x}\right)$ [14, (3)]. So the convergence
conditions following (B.5) are:

$\left\\{(x,y)\in\mathbb{R}^{2}\text{ }|\text{
}|y-1|<|y|,|y-1|+|y|<|x|\right\\}$

Figure 5: Convergence region of the $F_{Q}$ function

Olsson gives formula [14, (56)]

$F_{P}(a,b_{1},b_{2};c_{1},c_{2};x,y)=F_{PR}(a,b_{1},b_{2};c_{1},c_{2};x,y)+\\\
+\dfrac{\Gamma(a+b_{2}-c_{2}+1)\Gamma(a+b_{1}-c_{1}-b_{2})\Gamma(c_{1})}{\Gamma(a)\Gamma(b_{1})\Gamma(a-c_{2}+1)}x^{b_{1}-c_{1}}y^{-b_{2}}(1-x)^{c_{1}-b_{1}+b_{2}-a}\\\
\times F_{3}\left(\begin{array}[]{c}1-b_{1},b_{2};c_{1}-b_{1},b_{2}-c_{2}+1\\\
c_{1}-b_{1}+b_{2}-a+1\end{array};\dfrac{x-1}{x},\dfrac{1-x}{y}\right)$ (B.9)

where $F_{PR}$ is the part of the function $F_{P}$ which is regular at
$(1,1)$. For the function $F_{PR}$ he gave the formula [14, (57)]

$F_{PR}(a,b_{1},b_{2};c_{1},c_{2};x,y)=\dfrac{\Gamma(a+b_{2}-c_{2}+1)\Gamma(c_{1}-b_{1}+b_{2}-a)\Gamma(c_{1})}{\Gamma(a)\Gamma(c_{1}-b_{1}+b_{2}-c_{2}+1)\Gamma(c_{1}+b_{2}-a)}\\\
\qquad\qquad\qquad\qquad\times\sum\limits_{i=0}^{\infty}\sum\limits_{j=0}^{\infty}\dfrac{(a-c_{2}+1)_{i}(b_{1})_{i}(b_{2})_{j}}{(a+b_{1}-c_{1}-b_{2}+1)_{i}}\dfrac{1}{i!\,j!}(1-x)^{i}(1-y)^{j}\\\
\times\,\mbox{}_{3}F_{2}\\!\left(\genfrac{}{}{0.0pt}{}{b_{2}-c_{2}+1,c_{1}-b_{1}+b_{2}-a-i,c_{1}-a-j}{c_{1}-b_{1}+b_{2}-c_{2}+1,c_{1}+b_{2}-a};1\right)$
(B.10)

where the ${}_{3}F_{2}$ function is absolute convergent if
$\operatorname{Re}(a)>0$. He then referred to the paper [10, (12)] in which
the convergence region of the double summation is given as

$|x-1|+|y-1|<1$

The region of convergence of the $F_{3}$ function is:

$\left\\{\left|\dfrac{x-1}{x}\right|<1\wedge\left|\dfrac{1-x}{y}\right|<1\right\\}$

The point $(1,1)$ is inside both regions.

Figure 6: Convergence region of the $F_{PR}$ function

Figure 7: Convergence region of the $F_{3}$ function

## References

  * [1] Appell, P. Sur les séries hypergéometriques de deux variables et sur de équations différentielles linéaires aux dérivées partielles. C.R. Acad. Sci. Paris 90 (1880), 296–298.
  * [2] Appell, P., Kampé de Fériet, M.J. Fonctions hypergéometriques et hypersphériques polynomes d’Hermite. Gauthier-Villars, 1926.
  * [3] Bailey, W.N. Generalized hypergeometric series. Cambridge University Press, 1935.
  * [4] Diekema, E. The fractional orthogonal derivative for functions of one and two variables. Thesis, University of Amsterdam (2018). https://pure.uva.nl/ws/files/29401750/Thesis.pdf
  * [5] Diekema, E., Koornwinder, T.H. Integral representations for Horn’s $H_{2}$ function and Olsson’s $F_{P}$ function. Kyushu J. Math. 73(1) (2019) 1–24.
  * [6] Dunkl, C.F., Xu, Y. Orthogonal polynomials in several variables. Cambridge University Press, Second edition, 2014.
  * [7] Erdélyi, A. Higher transcendental functions, Vol. I. McGraw-Hill, 1953.
  * [8] Exton, H. Multiple hypergeometric functions and applications. Ellis Horwood, 1976.
  * [9] Fackerel, E.D., Littler, R.A. Polynomials biorthogonal to Appell’s polynomials. Bull. Austr. Math. Soc. 11 (1974) 181–195.
  * [10] Hahne, G.E. Analytic continuation of Appell’s hypergeometric series $F_{2}$ to the vincinity of the singular point $x=1,y=1$. J. Math. Phys. (10) 3 (1969) 524–531.
  * [11] Horn, J. Ueber die Convergenz der Hypergeometrische Reihen Zweier und Dreier Veränderlichen. Math. Ann. 34 (1889), 544–600.
  * [12] Horn, J. Hypergeometrische Funktionen Zweier Veränderlichen. Math. Ann. 105 (1931), 381–407.
  * [13] Olsson, P.O.M. A higher order solution of the partial differential equations associated with Appell’s hypergeometric function $F_{2}(a,b_{1},b_{2},c_{1},c_{2};x_{1},x_{2})$. Ark. Fys. B33 (1967), 433–442.
  * [14] Olsson, P.O.M. On the integration of the differential equations of five-parametric double hypergeometric functions of second order. J. Math. Phys. 18 (1977), 1285–1294.
  * [15] Slater, L.J. Generalized hypergeometric functions. Cambridge University Press, 1966.
  * [16] Suetin, P.K. Orthogonal polynomials in two variables. English translation of Russian original of 1988. Gordon and Breach, 1999.