Hypergeometric Differential Equation and New Identities for the Coefficients

Hypergeometric Differential Equation
and New Identities for the Coefficients
of Nørlund and Bühring\star\star$\star$This paper is a contribution to the
Special Issue on Orthogonal Polynomials, Special Functions and Applications.
The full collection is available at
http://www.emis.de/journals/SIGMA/OPSFA2015.html

Dmitrii KARP †‡ and Elena PRILEPKINA †‡ D. Karp and E. Prilepkina † Far
Eastern Federal University, 8 Sukhanova Str., Vladivostok, 690950, Russia ‡
Institute of Applied Mathematics, Far Eastern Branch of the Russian Academy of
Sciences,
‡ 7 Radio Str., Vladivostok, 690041, Russia<EMAIL_ADDRESS>pril-
<EMAIL_ADDRESS>http://dmkrp.wordpress.com

Received February 25, 2016, in final form May 15, 2016; Published online May
21, 2016

The fundamental set of solutions of the generalized hypergeometric
differential equation in the neighborhood of unity has been built by Nørlund
in 1955. The behavior of the generalized hypergeometric function in the
neighborhood of unity has been described in the beginning of 1990s by Bühring,
Srivastava and Saigo. In the first part of this paper we review their results
rewriting them in terms of Meijer’s $G$-function and explaining the
interconnections between them. In the second part we present new formulas and
identities for the coefficients that appear in the expansions of Meijer’s
$G$-function and generalized hypergeometric function around unity. Particular
cases of these identities include known and new relations for Thomae’s
hypergeometric function and forgotten Hermite’s identity for the sine
function.

generalized hypergeometric function; hypergeometric differential equation;
Meijer’s $G$-function; Bernoulli polynomials; Nørlund’s coefficients;
Bühring’s coefficients

33C20; 33C60; 34M35

##  1 Introduction

We will use standard notation $\mathbb{Z}$, $\mathbb{N}$ and $\mathbb{C}$ to
denote integer, natural and complex numbers, respectively;
$\mathbb{N}_{0}=\mathbb{N}\cup\\{0\\}$. The hypergeometric differential
equation ($D=z\frac{d}{dz}$)

$\displaystyle\left\\{(D-a_{1})(D-a_{2})\cdots(D-a_{p})-z(D+1-b_{1})(D+1-b_{2})\cdots(D+1-b_{p})\right\\}y=0$
(1.1)

for $p=2$ was first considered by Euler and later studied by Gauss, Kummer,
Riemann, Papperitz and Schwarz, among others, see [2, Section 2.3]. For
general $p>2$ it was probably first investigated by Thomae [42]. Note that our
choice of parameters differs slightly from that of [3, (16.8.3)]. As will
become apparent in the sequel, this choice is more convenient if the solution
is to be built in terms of Meijer’s $G$-functions and not generalized
hypergeometric functions. Equation (1.1) is of Fuchsian type and has three
regular singularities located at the points $0$, $1$, $\infty$. The local
exponents read [5, (2.6)–(2.8)]

$\displaystyle a_{1},a_{2},\ldots,a_{p}\quad$ $\displaystyle\text{at}~{}z=0,$
$\displaystyle 0,1,2,\ldots,p-2,\sum(b_{i}-a_{i})-1\quad$
$\displaystyle\text{at}~{}z=1,$ $\displaystyle
1-b_{1},1-b_{2},\ldots,1-b_{p}\quad$ $\displaystyle\text{at}~{}z=\infty.$

The fundamental sets of solutions around the points $z=0$ and $z=\infty$ were
found by Thomae in [42] and are expressed in terms of the generalized
hypergeometric series

$\displaystyle{{}_{p}F_{q}}\left(\begin{matrix}\mathbf{a}\\\
\mathbf{b}\end{matrix}\,\vline\,z\right)={{}_{p}F_{q}}\left(\mathbf{a};\mathbf{b};z\right)=\sum\limits_{n=0}^{\infty}\frac{(a_{1})_{n}(a_{2})_{n}\cdots(a_{p})_{n}}{(b_{1})_{n}(b_{2})_{n}\cdots(b_{q})_{n}n!}z^{n},$

where $(a)_{n}=\Gamma(a+n)/\Gamma(a)$ denotes the rising factorial and
$\mathbf{a}=(a_{1},\ldots,a_{p})$, $\mathbf{b}=(b_{1},\ldots,b_{q})$ are
complex vectors such that $-b_{i}\notin\mathbb{N}_{0}$ for $i=1,\ldots,q$. If
the components of $\mathbf{a}$ are distinct modulo $\mathbb{Z}$, then the
basis of solutions near $z=0$ is given by [33, (1.13)]

$\displaystyle
z^{a_{k}}{{}_{p}F_{p-1}}\left(\begin{matrix}1-\mathbf{b}+a_{k}\\\
1-\mathbf{a}_{[k]}+a_{k}\end{matrix}\,\vline\,z\right),\qquad k=1,\ldots,p,$

where $\mathbf{a}_{[k]}$ signifies the vector $\mathbf{a}$ with the element
$a_{k}$ omitted and $\mathbf{a}+\alpha$ is understood as $(a_{1}+\alpha$,
$\ldots,a_{p}+\alpha)$. If the components of $\mathbf{b}$ are distinct modulo
$\mathbb{Z}$ then the fundamental system of solutions near $z=\infty$ is given
by [33, (1.14)]

$\displaystyle
z^{b_{k}-1}{{}_{p}F_{p-1}}\left(\begin{matrix}1+\mathbf{a}-b_{k}\\\
1+\mathbf{b}_{[k]}-b_{k}\end{matrix}\,\vline\,\frac{1}{z}\right),\qquad
k=1,\ldots,p.$

Finally, the fundamental set of solutions around the point $z=1$ was found by
the remarkable Danish mathematician Niels Erik Nørlund in his milestone work
[33]. These facts are well-known and have been frequently cited in the
literature. It seems to be less known that the fundamental solutions
constructed by Nørlund can be expressed in terms of $G$-function introduced
some 15 years earlier by Meijer, see [30]. In fact, Meijer himself studied the
hypergeometric differential equation more general than (1.1) and built the
basis of solutions in the neighborhood of zero and infinity in terms of
$G$-functions. The connection between Nørlund’s solutions and Meijer’s
$G$-function was observed by Marichev and Kalla [27] and Marichev [26] but
these two papers remained largely unnoticed. The fact that each solution
around $z=1$ must equal to a linear combination of fundamental solutions
around $z=0$ is reflected for $p=2$ in the following identity due to Gauss [2,
Theorem 2.3.2]

$\displaystyle{{}_{2}F_{1}}\left(\begin{matrix}\alpha_{1},\alpha_{2}\\\
\beta\end{matrix};1-z\right)=\frac{\Gamma(\beta)\Gamma(\beta-\alpha_{1}-\alpha_{2})}{\Gamma(\beta-\alpha_{1})\Gamma(\beta-\alpha_{2})}{{}_{2}F_{1}}\left(\begin{matrix}\alpha_{1},\alpha_{2}\\\
\alpha_{1}+\alpha_{2}-\beta+1\end{matrix};z\right)$
$\displaystyle\hphantom{{{}_{2}F_{1}}\left(\begin{matrix}\alpha_{1},\alpha_{2}\\\
\beta\end{matrix};1-z\right)=}{}+\frac{\Gamma(\beta)\Gamma(\alpha_{1}+\alpha_{2}-\beta)}{\Gamma(\alpha_{1})\Gamma(\alpha_{2})}z^{\beta-\alpha_{1}-\alpha_{2}}{{}_{2}F_{1}}\left(\begin{matrix}\beta-\alpha_{1},\beta-\alpha_{2}\\\
\beta-\alpha_{1}-\alpha_{2}+1\end{matrix};z\right).$ (1.2)

The fact that each solution around $z=0$ must equal to a linear combination of
fundamental solutions around $z=1$ leads to essentially the same identity with
$z$ replaced by $1-z$. For $p>2$, however the above connection formula has
_two different generalizations_. One of them is the expansion (2.4) of
$G^{p,0}_{p,p}$ (a solution around $z=1$) into a sum of ${}_{p}F_{p-1}$
(fundamental solutions around $z=0$) which can also be obtained by applying
the residue theorem to the definition of $G$-function and is found in standard
references on $G$-function, see, for instance, [3, (16.17.2)]. The other one
is the expansion of ${}_{p}F_{p-1}$ into a sum of $G^{p,0}_{p,p}$ with $p-1$
instances of $G^{2,p}_{p,p}$ as given by formula (2.3) below. This expansion,
discovered by Nørlund without any mentioning of $G$-function has been
reproduced in [27] in a slightly different form but seems to be forgotten
afterwards. A related expansion, again with no reference to $G$-function, has
been then found in [7, 8] by Bühring whose goal was to describe the behavior
of ${}_{p}F_{p-1}(z)$ near $z=1$. The logarithmic cases have been studied by
Saigo and Srivastava in [37]. Let us also mention that the monodromy group of
the generalized hypergeometric differential equation has been constructed by
Beukers and Heckman in [5].

This paper is organized as follows. Section 2 is of survey nature: we review
the results of Nørlund and Bühring and rewrite them in terms of $G$-function.
We also reveal connections between Nørlund’s and Bühring’s expansions and
relate them to various results obtained in statistics literature. In Section 3
we present some new formulas for Nørlund’s coefficients and utilize the above
mentioned connections to derive various new identities for Nørlund’s and
Bühring’s coefficients which reduce to hypergeometric and trigonometric
identities for small $p$. The most striking of these identities,

$\displaystyle\sum\limits_{k=1}^{p}\frac{\prod\limits_{i=1}^{p}\sin(\beta_{i}-\alpha_{k})}{\prod\limits_{\begin{subarray}{c}i=1\\\
i\neq{k}\end{subarray}}^{p}\sin(\alpha_{i}-\alpha_{k})}=\sin\left(\sum\limits_{k=1}^{p}(\beta_{k}-\alpha_{k})\right),$
(1.3)

can be viewed as a generalization of Ptolemy’s theorem: for a quadrilateral
inscribed in a circle the product of the lengths of its diagonals is equal to
the sum of the products of the lengths of the pairs of opposite sides. This
theorem can be written in trigonometric form which yields precisely the above
identity for $p=2$. For $p=3$ formula (1.3) was first presented without proof
by Glaisher at a 1880 conference and then proved in full generality by Hermite
in his 1885 paper [20]. The same paper contains a number of other amazing
identities completely forgotten until the recent article [21] by Johnson,
where many interesting mathematical and historical details can be found.
Furthermore, if we substitute $\sin(z)$ by $z$ in every occurrence of $\sin$
in (1.3), we obtain an identity discovered by Gosper, Ismail and Zhang in [19]
and known as non-local derangement identity. It has been recently used by
Feng, Kuznetsov and Yang to find new formulas for sums of products of
generalized hypergeometric functions [18].

When this paper was nearly finished E. Scheidegger published a preprint [38]
that also deals with the hypergeometric differential equation (1.1) and builds
largely on the works of Nørlund [33] and Bühring [8]. Scheidegger suggests a
new basis of solutions around $1$ and presents its series expansion. Main
emphasis in [38] is on the logarithmic cases (parameter differences are
integers). Scheidegger’s work is motivated by applications in certain one-
parameter families of Calabi–Yau manifolds, known as the mirror quartic and
the mirror quintic.

##  2 Results of Nørlund and Bühring revisited

###  2.1 Fundamental solutions around unity

The simple observation that nevertheless seems to be largely overlooked in the
literature (except for [26, 27]) lies in the fact that Nørlund’s solutions to
(1.1) can be expressed in terms of Meijer’s $G$-function. Let us remind its
definition first. Suppose $0\leq{m}\leq{q}$, $0\leq{n}\leq{p}$ are integers
and $\mathbf{a}$, $\mathbf{b}$ are arbitrary complex vectors, such that
$a_{i}-b_{j}\notin\mathbb{N}$ for all $i=1,\ldots,n$ and $j=1,\ldots,m$.
Meijer’s $G$-function is defined by the Mellin–Barnes integral of the form
(see [17, Section 5.3], [25, Chapter 1], [36, Section 8.2] or [3, Section
16.17]),

$\displaystyle G^{m,n}_{p,q}\left(z\,\vline\,\begin{matrix}\mathbf{a}\\\
\mathbf{b}\end{matrix}\right):=\frac{1}{2\pi{i}}\int_{\mathcal{L}}\frac{\Gamma(b_{1}+s)\cdots\Gamma(b_{m}+s)\Gamma(1-a_{1}-s)\cdots\Gamma(1-a_{n}-s)z^{-s}}{\Gamma(a_{n+1}+s)\cdots\Gamma(a_{p}+s)\Gamma(1-b_{m+1}-s)\cdots\Gamma(1-b_{q}-s)}ds,$
(2.1)

where the contour $\mathcal{L}$ is a simple loop that separates the poles of
the integrand of the form $b_{jl}=-b_{j}-l$, $l\in\mathbb{N}_{0}$, leaving
them on the left from the poles of the form $a_{ik}=1-a_{i}+k$,
$k\in\mathbb{N}_{0}$, leaving them on the right [25, Section 1.1]. It may have
one of the three forms $\mathcal{L}_{-}$, $\mathcal{L}_{+}$ or
$\mathcal{L}_{i\gamma}$ described below. Choose any

$\displaystyle\begin{split}&\varphi_{1}<\min\\{-\Im{b_{1}},\ldots,-\Im{b_{m}},\Im(1-a_{1}),\ldots,\Im(1-a_{n})\\},\\\
&\varphi_{2}>\max\\{-\Im{b_{1}},\ldots,-\Im{b_{m}},\Im(1-a_{1}),\ldots,\Im(1-a_{n})\\}\end{split}$

and arbitrary real $\gamma$. The contour $\mathcal{L}_{-}$ is a left loop
lying in the horizontal strip $\varphi_{1}\leq\Im{s}\leq\varphi_{2}$. It
starts at the point $-\infty+i\varphi_{1}$, terminates at the point
$-\infty+i\varphi_{2}$ and coincides with the boundary of the strip for
sufficiently large $|s|$. Similarly, the contour $\mathcal{L}_{+}$ is a right
loop lying in the same strip, starting at the point $+\infty+i\varphi_{1}$ and
terminating at the point $+\infty+i\varphi_{2}$. Finally, the contour
$\mathcal{L}_{i\gamma}$ starts at $\gamma-i\infty$, terminates at
$\gamma+i\infty$ and coincides with the line $\Re{s}=\gamma$ for all
sufficiently large $|s|$. The power function $z^{-s}$ is defined on the
Riemann surface of the logarithm, so that

$\displaystyle z^{-s}=\exp(-s\\{\log|z|+i\arg(z)\\})$

and $\arg(z)$ is allowed to take any real value. Hence, $G^{m,n}_{p,q}(z)$ is
also defined on the Riemann surface of the logarithm. In order that the above
definition be consistent one needs to prove that the value of the integral
remains intact if convergence takes place for several different contours.
Alternatively, one may split the parameter space into nonintersecting subsets
and stipulate which contour should be used in each subset. Another key issue
that must be addressed with regard to the above definition is whether the
integral in (2.1) equals the sum of residues of the integrand and, if yes, on
which side of $\mathcal{L}$ the residues are to be counted. This is important
for both theoretical considerations (expressing $G$-function in terms of
hypergeometric functions) and especially for actually computing the value of
$G$-function (although numerical contour integration can also be employed). In
this paper we will only need the $G$-function of the form $G^{m,n}_{p,p}$. For
this particular type of $G$-function the solutions to the above problems seem
to be rather complete. We placed further details regarding the definition and
answers to the above questions for $G^{m,n}_{p,p}$ in Appendix A.

A simple property of Meijer’s $G$-function implied by its definition (2.1)
which will be frequently used without further mentioning is given by [36,
(8.2.2.15)]

$\displaystyle
z^{\alpha}G^{m,n}_{p,q}\left(z\,\vline\,\begin{matrix}\mathbf{a}\\\
\mathbf{b}\end{matrix}\right)=G^{m,n}_{p,q}\left(z\,\vline\,\begin{matrix}\mathbf{a}+\alpha\\\
\mathbf{b}+\alpha\end{matrix}\right),$ (2.2)

where $\alpha\in\mathbb{C}$ and $\mathbf{a}+\alpha$ is understood as
$(a_{1}+\alpha,\ldots,a_{p}+\alpha)$.

The functions $G^{p,0}_{p,p}$ and $G^{2,p}_{p,p}$ will play a particularly
important role in this paper, so that we found it useful to cite the known
explicit expressions for small $p$. If $p=1$ then [36, formula (8.4.2.3)]

$\displaystyle G^{1,0}_{1,1}\left(z\,\vline\,\begin{matrix}b\\\
a\end{matrix}\right)=\frac{z^{a}(1-z)_{+}^{b-a-1}}{\Gamma(b-a)},$

where $(x)_{+}=x$ for $x\geq 0$ and $0$ otherwise. If $p=2$ we have [36,
formula (8.4.49.22)]

$\displaystyle G^{2,0}_{2,2}\left(z\,\vline\,\begin{matrix}b_{1},b_{2}\\\
a_{1},a_{2}\end{matrix}\right)=\frac{z^{a_{2}}(1-z)_{+}^{b_{1}+b_{2}-a_{1}-a_{2}-1}}{\Gamma(b_{1}+b_{2}-a_{1}-a_{2})}\,{}_{2}F_{1}\left(\begin{matrix}b_{1}-a_{1},b_{2}-a_{1}\\\
b_{1}+b_{2}-a_{1}-a_{2}\end{matrix};1-z\right),$

and [36, formula (8.4.49.20)]

$\displaystyle G^{2,2}_{2,2}\left(z\,\vline\,\begin{matrix}b_{1},b_{2}\\\
a_{1},a_{2}\end{matrix}\right)=\frac{z^{a_{1}}\Gamma(1+a_{1}-b_{1})\Gamma(1+a_{1}-b_{2})\Gamma(1+a_{2}-b_{2})\Gamma(1+a_{2}-b_{1})}{\Gamma(2+a_{1}+a_{2}-b_{1}-b_{2})}$
$\displaystyle\hphantom{G^{2,2}_{2,2}\left(z\,\vline\,\begin{matrix}b_{1},b_{2}\\\
a_{1},a_{2}\end{matrix}\right)=}{}\times{{}_{2}F_{1}}\left(\begin{matrix}1+a_{1}-b_{1},1+a_{1}-b_{2}\\\
2+a_{1}+a_{2}-b_{1}-b_{2}\end{matrix};1-z\right).$

If $p=3$ we have [36, formula (8.4.51.2)]

$\displaystyle\begin{split}&G^{3,0}_{3,3}\left(z\,\vline\,\begin{matrix}b_{1},b_{2},b_{3}\\\
a_{1},a_{2},a_{3}\end{matrix}\right)=\frac{z^{a_{1}+a_{2}-b_{1}-1}(1-z)_{+}^{b_{1}+b_{2}+b_{3}-a_{1}-a_{2}-a_{3}-1}}{\Gamma(b_{1}+b_{2}+b_{3}-a_{1}-a_{2}-a_{3})}\\\
&\qquad{}\times{F_{3}}(b_{1}-a_{2},b_{3}-a_{3};b_{1}-a_{1},b_{2}-a_{3};b_{1}+b_{2}+b_{3}-a_{1}-a_{2}-a_{3};1-1/z,1-z),\end{split}$

where $F_{3}$ is Appell’s hypergeometric function of two variables [36, index
of functions].

Having made these preparations we can formulate Nørlund’s result regarding a
fundamental solution of (1.1) in the neighborhood of $1$. First, introduce the
following notation:

$\displaystyle\psi_{m}=\sum_{i=1}^{m}(b_{i}-a_{i}),\qquad 1\leq m\leq p,$
$\displaystyle\mathbf{a}=(a_{1},\ldots,a_{p}),\qquad\mathbf{b}=(b_{1},\ldots,b_{p}),\qquad\mathbf{a}_{[k]}=(a_{1},\ldots,a_{k-1},a_{k+1},\ldots,a_{p}),$
$\displaystyle\sin(\mathbf{a})=\sin{a_{1}}\sin{a_{2}}\cdots\sin{a_{p}},\qquad\Gamma(\mathbf{a})=\Gamma(a_{1})\Gamma(a_{2})\cdots\Gamma(a_{p})$

and let $\mathbf{a}_{[k,s]}$ denote the vector $\mathbf{a}$ with the elements
$a_{k}$ and $a_{s}$ removed. We write $\Re(\mathbf{a})>0$ for $\Re(a_{i})>0$,
$i=1,\dots,p$. The next theorem is implicit in [33].

######  Theorem 2.1 (Nørlund).

Suppose $k,s\in\\{1,\ldots,p\\}$ and

$\displaystyle
u_{s}(z)=G^{p,0}_{p,p}\left(z\,\vline\,\begin{matrix}\mathbf{b}\\\
\mathbf{a}\end{matrix}\right),\qquad
u_{k}(z)=G^{2,p}_{p,p}\left(z\,\vline\,\begin{matrix}\mathbf{b}\\\
a_{k},a_{s},\mathbf{a}_{[k,s]}\end{matrix}\right)\qquad\text{for}\quad
k\neq{s}.$

Then the set $\\{u_{k}(z)\\}_{k=1}^{p}$ forms the fundamental system of
solutions of (1.1) in the neighborhood of $z=1$ and

$\displaystyle\sin(\pi\psi_{p})\frac{z^{a_{s}}\Gamma(1-\mathbf{b}+a_{s})}{\Gamma(1-\mathbf{a}_{[s]}+a_{s})}\,{}_{p}F_{p-1}\left(\begin{matrix}1-\mathbf{b}+a_{s}\\\
1-\mathbf{a}_{[s]}+a_{s}\end{matrix};z\right)$ $\displaystyle\qquad{}=\pi
G^{p,0}_{p,p}\left(z\,\vline\,\begin{matrix}\mathbf{b}\\\
\mathbf{a}\end{matrix}\right)-\frac{1}{\pi}\sum_{\begin{subarray}{c}k=1\\\
k\neq{s}\end{subarray}}^{p}\frac{\sin(\pi(\mathbf{b}-a_{k}))}{\sin(\pi(\mathbf{a}_{[k,s]}-a_{k}))}G^{2,p}_{p,p}\left(z\,\vline\,\begin{matrix}\mathbf{b}\\\
a_{k},a_{s},\mathbf{a}_{[k,s]}\end{matrix}\right).$ (2.3)

###### Proof.

Formula (2.3) is a rewriting of [33, (5.40)]. Indeed, first compare [33,
(2.44)] with the definition of $G$-function (2.1) to see how Nørlund’s
$\xi_{n}$ is expressed by $G^{p,0}_{p,p}$. Alternatively, this connection
follows on comparing the Mellin transforms (A.3) and [33, (2.18)]. Further,
use [33, (5.45)] to expand Nørlund’s $\varphi_{s}$ and [33, (5.7)] to express
$y_{k,s}$ in terms of $G^{2,p}_{p,p}$ (see also (2.10) below). ∎

######  Remark 2.2.

Formula (2.3) can now be viewed as the reflection of the fact that any $p+1$
solutions must be linearly dependent so that any solution in the neighborhood
of $0$ can be expressed in terms of the fundamental set of solutions around
$1$. This formula extends the connection formula (1.2) for the Gauss
hypergeometric function to which it reduces when $p=2$. A closely related
formula has been also discovered by Bühring in two papers [7] (for $p=3$) and
[8] (for general $p$). See formula (2.19) below.

######  Remark 2.3.

Marichev and Kalla reproduced formula [33, (5.40)] in [27, (11)] and gave
expressions for its components in terms of $G$-function in [27, (27), (33),
(36)].

Expansion of the solution $G^{p,0}_{p,p}$ in terms of the fundamental
solutions around $z=0$ that compliments (2.3) coincides with the well-known
expansion obtained from the definition of $G$-function by the residue theorem.
Namely, if the elements of the vector $\mathbf{a}$ are different modulo
integers, formula (A.1) applied to the function $G^{p,0}_{p,p}$ takes the form
(see also [26, formula (34)] and [36, 8.2.2.3]):

$\displaystyle G^{p,0}_{p,p}\left(z\,\vline\,\begin{matrix}\mathbf{b}\\\
\mathbf{a}\end{matrix}\right)=\sum\limits_{k=1}^{p}z^{a_{k}}\frac{\Gamma(\mathbf{a}_{[k]}-a_{k})}{\Gamma(\mathbf{b}-a_{k})}\,{}_{p}F_{p-1}\left(\begin{matrix}1-\mathbf{b}+a_{k}\\\
1-\mathbf{a}_{[k]}+a_{k}\end{matrix}\,\vline\,z\right).$ (2.4)

According to Theorem A.1(a) this formula holds for $|z|<1$.

###  2.2 Expansion of $\boldsymbol{G^{p,0}_{p,p}}$ in the neighborhood of
unity

Among many other results contained in [33] Nørlund showed that the series

$\displaystyle G^{p,0}_{p,p}\left(z\,\vline\,\begin{matrix}\mathbf{b}\\\
\mathbf{a}\end{matrix}\right)=\frac{z^{a_{k}}(1-z)^{\psi_{p}-1}}{\Gamma(\psi_{p})}\sum_{n=0}^{\infty}\frac{g_{p}^{k}(n)}{(\psi_{p})_{n}}(1-z)^{n}$
(2.5)

represents a solution in the neighborhood of $z=1$ corresponding to the local
exponent $\psi_{p}-1$ (see (2.1) for the definition of $\psi_{p}$) if this
number is not a negative integer. Formula (2.5) holds in the disk $|1-z|<1$
for all $-\psi_{p}\notin\mathbb{N}_{0}$ and each $k=1,2,\ldots,p$. For
$-\psi_{p}=l\in\mathbb{N}_{0}$ we have by taking limit in (2.5) (see [33,
(1.34)]):

$\displaystyle G^{p,0}_{p,p}\left(z\,\vline\,\begin{matrix}\mathbf{b}\\\ \
a\end{matrix}\right)=z^{a_{k}}\sum\limits_{n=0}^{\infty}\frac{g^{k}_{p}(n+l+1)}{n!}(1-z)^{n},\qquad
k=1,2,\ldots,p.$ (2.6)

The value of $k$ in (2.5) and (2.6) can be chosen arbitrarily from the set
$\\{1,2,\ldots,p\\}$. This choice affects the first factor $z^{a_{k}}$ and the
coefficients $g_{p}^{k}(n)$ while the left-hand side is of course independent
of $k$.

Let us present Nørlund’s formulas for the coefficients $g^{k}_{p}(n)$. First,
by applying the Frobenius ansatz to the differential equation (1.1) Nørlund
demonstrated that $g^{k}_{p}(n)$ satisfy the $p$-th order difference equation
with polynomial coefficients given by

$\displaystyle\sum_{i=0}^{p-1}P_{p-i}(n\,|\,\mathbf{a},\mathbf{b})g^{k}_{p}(n+i)=0,\qquad
n=1,2,\ldots,$ (2.7)

where $P_{m}(z\,|\,\mathbf{a},\mathbf{b})$ is a polynomial in $z$ of degree
$m$ with coefficients dependent on the parameters $\mathbf{a}$ and
$\mathbf{b}$. It is expressed in terms of the polynomials

$\displaystyle Q(z)=\prod_{i=1}^{p}(z-a_{i}),\qquad
R(z)=\prod_{i=1}^{p}(z+1-b_{i})$

as follows [33, (1.28)]:

$\displaystyle P_{1}(z\,|\,\mathbf{a},\mathbf{b})=p-1+z,$ $\displaystyle
P_{j}(z\,|\,\mathbf{a},\mathbf{b})=\frac{(-1)^{j}}{(p-j-1)!}\Delta^{p-j-1}Q(\psi_{p}+a_{k}+z)-\frac{(-1)^{j}}{(p-j)!}\Delta^{p-j}R(\psi_{p}-1+a_{k}+z)$

for $j=2,\ldots,p-1$, and

$\displaystyle
P_{p}(z\,|\,\mathbf{a},\mathbf{b})=-(-1)^{p}R(\psi_{p}-1+a_{k}+z).$

Here the difference is understood as the forward difference
$\Delta{Q(z)}=Q(z+1)-Q(z)$, $\Delta^{m}{Q(z)}=\Delta(\Delta^{m-1}{Q(z)})$. The
initial values $g^{k}_{p}(0)$, $g^{k}_{p}(1)$, $\ldots$, $g^{k}_{p}(p-1)$ for
the recurrence (2.7) are found by solving the next triangular system:

$\displaystyle\begin{split}&g^{k}_{p}(0)=1,\\\
&g^{k}_{p}(1)+P_{2}(2-p\,|\,\mathbf{a},\mathbf{b})g^{k}_{p}(0)=0,\\\
&2g^{k}_{p}(2)+P_{2}(3-p\,|\,\mathbf{a},\mathbf{b})g^{k}_{p}(1)+P_{3}(3-p\,|\,\mathbf{a},\mathbf{b})g^{k}_{p}(0)=0,\\\
&3g^{k}_{p}(3)+P_{2}(4-p\,|\,\mathbf{a},\mathbf{b})g^{k}_{p}(2)+P_{3}(4-p\,|\,\mathbf{a},\mathbf{b})g^{k}_{p}(1)+P_{4}(4-p\,|\,\mathbf{a},\mathbf{b})g^{k}_{p}(0)=0,\\\
&\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\\\
&(p-1)g^{k}_{p}(p-1)+P_{2}(0\,|\,\mathbf{a},\mathbf{b})g^{k}_{p}(p-2)+P_{3}(0\,|\,\mathbf{a},\mathbf{b})g^{k}_{p}(p-3)+\\!\cdots\\!+P_{p}(0\,|\,\mathbf{a},\mathbf{b})g^{k}_{p}(0)=0.\end{split}$

Another method to compute the coefficients $g^{k}_{p}(n)$ discovered by
Nørlund is the following recurrence in $p$ [33, (2.7)]:

$\displaystyle g^{1}_{1}(n)=\begin{cases}1,&n=0,\\\ 0,&n\geq
1,\end{cases}\qquad
g^{p}_{p}(n)=\sum_{j=0}^{n}\frac{(b_{p}-a_{k})_{n-j}}{(n-j)!}(\psi_{p-1}+j)_{n-j}g^{k}_{p-1}(n).$
(2.8)

The last formula can be applied for any $k=1,2,\ldots,p-1$ without affecting
the left-hand side. The value of $g^{k}_{p}(n)$ for $k\neq{p}$ can then be
obtained by exchanging the roles of $a_{p}$ and $a_{k}$ in the resulting
expression for $g^{p}_{p}(n)$. Alternatively, Nørlund gives the following
connection formula [33, (1.35)]:

$\displaystyle
g^{k}_{p}(n)=\sum\limits_{j=0}^{n}\frac{(a_{k}-a_{l})_{n-j}}{(n-j)!}(\psi_{p}+j)_{n-j}g^{l}_{p}(j),\qquad
k,l=1,\ldots,p,\qquad k\neq{l}.$

Furthermore, he solved the above recurrence to obtain [33, (2.11)]:

$\displaystyle
g^{p}_{p}(n)=\sum\limits_{0\leq{j_{1}}\leq{j_{2}}\leq\cdots\leq{j_{p-2}}\leq{n}}\prod\limits_{m=1}^{p-1}\frac{(\psi_{m}+j_{m-1})_{j_{m}-j_{m-1}}}{(j_{m}-j_{m-1})!}(b_{m+1}-a_{m})_{j_{m}-j_{m-1}},$
(2.9)

where, as before, $\psi_{m}=\sum_{i=1}^{m}(b_{i}-a_{i})$ and $j_{0}=0$,
$j_{p-1}=n$. This formula shows explicitly that $g^{p}_{p}(n)$ does not depend
on $a_{p}$. Observe that

  * •

the coefficient $g^{k}_{p}(n)$ is a symmetric polynomial in the components of
$\mathbf{b}$ and $\mathbf{a}_{[k]}$ (separately);

  * •

the summation in (2.9) is over all Young diagrams that fit $n\times{(p-2)}$
box.

Symmetry follows from expansion (2.5) and the invariance of $G$-function with
respect to permutation the elements of $\mathbf{a}$ and $\mathbf{b}$. In
Section 3 below we find another way to compute the coefficients $g_{p}^{k}(n)$
in terms of generalized Bernoulli polynomials and we further give explicit
formulas for $g_{p}^{k}(n)$, $n=1,2,3$. Connection of these coefficients to
combinatorics deserves further investigation.

######  Remark 2.4.

Formulas (2.3) and (2.5) have been reproduced by Marichev [26, (4), (17)] and
Marichev and Kalla in [27, (12), (27)]. Surprisingly, these important formulas
seem to have been largely overlooked in the special function literature and
have never been included in any textbook on hypergeometric functions, except
for the reference book [36] by Prudnikov, Brychkov and Marichev. However, even
in this book the description of the behavior of $G^{p,0}_{p,p}(z)$ in the
neighborhood of $z=1$ contains an incorrect assertion in the case of non-
positive integer $\psi_{p}$, see [36, Section 8.2.2.59]. None of the
identities presented in Section 3 below are contained in [26, 27, 36].

######  Remark 2.5.

Another very prolific line of research that involves $G$-function forms an
important part of the statistics literature and began with 1932 paper of Wilks
[44], where he observed that the moments of many likelihood ratio criteria in
multivariate hypothesis testing are expressed in terms of product ratios of
gamma functions. Wilks introduced two types of integral equations, of which
“type B” is essentially equation (A.3) of the Appendix. He also noticed that
the solution of “type B integral equation” represents the probability density
of the product of independent beta distributed random variables. Wilks’ ideas
were elaborated in dozens of papers that followed, mainly concerned with
calculating and approximating the solution of “type B integral equation”. We
will just mention a few key contributions, where an interested reader may find
further references. In his 1939 paper [32] Nair derived the differential
equation (1.1) satisfied by the solution of “type B integral equation” thus
demonstrating that $G^{p,0}_{p,p}$ satisfies (1.1). This happened long before
Nørlund wrote his paper [33] and was also independent of Meijer’s work,
although Meijer already introduced $G$-function as a linear combination of
hypergeometric functions in 1936. Nair also was the first (among researchers
in statistics) to apply the inverse Mellin transform to the right-hand side of
(A.3) – the approach further developed by Consul in a series of papers between
1964 and 1969, where the connection to Meijer’s $G$-function was first
observed. See [13] and references therein. Mellin transform technique was then
utilized in a number of papers by Mathai who later rediscovered Nørlund’s
coefficients in a form similar to (2.9) in [29]. Mathai’s and other
contributions until 1973 are described in his survey paper [28]. In the same
period Springer and Thompson independently expressed the densities of products
and ratios of gamma and beta distributed random variables in terms of Meijer’s
$G$-function, see [39]. Davis [14] presented the matrix form of the
differential equation for $G^{p,0}_{p,p}$ and suggested the series solution
similar to (2.5) with coefficients found by certain recursive procedure.
Another noticeable contribution is due to Gupta and Tang who found two series
expansions for $G^{p,0}_{p,p}$ (again using “type B integral equation”
terminology), one of them equivalent to (2.5), and rediscovered recurrence
relation (2.8). See [40, 41] and references there. This line of research
continues until today, as evidenced, for example, by a series of papers by
Carlos Coelho with several co-authors, see [12] and references there.
Furthermore, Charles Dunkl rediscovered Nørlund’s recurrence (2.7) in his
preprint [16] dated 2013, where he again considers the probability density of
a product of beta distributed random variables. Independently, probability
distribution with $G$-function density has been found to be the stationary
distribution of certain Markov chains considered, for example, in actuarial
science and is known as Dufresne law in this context. See [11, 15] and
references therein. Let us also mention that $G$-function popped up recently
in the random matrix theory as correlation kernel of a determinantal point
process that governs singular values of products of $M$ rectangular random
matrices with independent complex Gaussian entries [1]. Curiously enough, none
of these authors cited Nørlund’s work. Let us also mention that in our recent
paper [23] conditions are given under which $G_{p,p}^{p,0}(e^{-x})$ is
infinitely divisible distribution on $[0,\infty)$.

###  2.3 Expansion of $\boldsymbol{G^{2,p}_{p,p}}$ in the neighborhood of
unity

Formula [33, (5.7)] shows that Nørlund’s function $y_{1,2}(x)$ defined by [33,
formula (5.2)] is expressed in terms of Meijer’s $G$-function as follows

$\displaystyle
y_{1,2}(x)=G^{2,p}_{p,p}\left(x\,\vline\,\begin{matrix}1-\alpha_{1},\ldots,1-\alpha_{p}\\\
\gamma_{1},\gamma_{2},\gamma_{3},\ldots,\gamma_{p}\end{matrix}\right).$ (2.10)

To define $y_{k,s}(x)$ the roles of $\gamma_{1}$, $\gamma_{2}$ are exchanged
with those of $\gamma_{k}$, $\gamma_{s}$. Nørlund found several expansions of
the functions $y_{k,s}(x)$ in hypergeometric polynomials which we cite below.
Changing Nørlund’s notation to ours according to the rule $x\mapsto{z}$,
$\gamma_{i}\mapsto{a_{i}}$, $\alpha_{i}\mapsto{1-b_{i}}$, formula [33, (5.3)]
takes the form

$\displaystyle\sin\pi(a_{s}-a_{i})G^{2,p}_{p,p}\left(z\,\vline\,\begin{matrix}\mathbf{b}\\\
a_{s},a_{i},\mathbf{a}_{[s,i]}\end{matrix}\right)+\sin\pi(a_{i}-a_{k})G^{2,p}_{p,p}\left(z\,\vline\,\begin{matrix}\mathbf{b}\\\
a_{i},a_{k},\mathbf{a}_{[i,k]}\end{matrix}\right)$
$\displaystyle\hphantom{\sin\pi(a_{s}-a_{i})G^{2,p}_{p,p}\left(z\,\vline\,\begin{matrix}\mathbf{b}\\\
a_{s},a_{i},\mathbf{a}_{[s,i]}\end{matrix}\right)}{}+\sin\pi(a_{k}-a_{s})G^{2,p}_{p,p}\left(z\,\vline\,\begin{matrix}\mathbf{b}\\\
a_{k},a_{s},\mathbf{a}_{[k,s]}\end{matrix}\right)=0$ (2.11)

for any distinct values of $s,i,k\in\\{1,\ldots,p\\}$, where
$\mathbf{a}_{[k,s]}$ denotes the vector $\mathbf{a}$ with elements $a_{k}$ and
$a_{s}$ removed. Expansions [33, (5.20), (5.22), (5.23), (5.31)] written in
terms of $G$-function take the form

$\displaystyle G^{2,p}_{p,p}\left(z\,\vline\,\begin{matrix}\mathbf{b}\\\
a_{k},a_{s},\mathbf{a}_{[k,s]}\end{matrix}\right)=\frac{z^{a_{k}}\Gamma(1-\mathbf{b}+a_{k})\Gamma(1-b_{1}+a_{s})\Gamma(1-b_{2}+a_{s})}{\Gamma(1-\mathbf{a}_{[k,s]}+a_{k})\Gamma(2+a_{k}+a_{s}-b_{1}-b_{2})}$
$\displaystyle\quad{}\times\sum_{n=0}^{\infty}\frac{(1-b_{1}+a_{k})_{n}(1-b_{2}+a_{k})_{n}}{(2+a_{k}+a_{s}-b_{1}-b_{2})_{n}}\,{}_{p-1}F_{p-2}\left(\begin{matrix}-n,1-\mathbf{b}_{[1,2]}+a_{k}\\\
1-\mathbf{a}_{[k,s]}+a_{k}\end{matrix}\,\vline\,z\right)$ (2.12)
$\displaystyle{}=\frac{z^{a_{k}}\Gamma(1-a_{k}+a_{s})\Gamma(1-\mathbf{b}+a_{k})}{\Gamma(1-\mathbf{a}_{[k,s]}+a_{k})}\sum_{n=0}^{\infty}\frac{(1-b_{1}+a_{k})_{n}}{(1-b_{1}+a_{s})_{n+1}}\,{}_{p}F_{p-1}\left(\begin{matrix}-n,1-\mathbf{b}_{[1]}+a_{k}\\\
1,1-\mathbf{a}_{[k,s]}+a_{k}\end{matrix}\,\vline\,z\right)\\!\\!\\!\\!\\!$
(2.13)
$\displaystyle{}=\frac{z^{a_{k}}\Gamma(1-b_{1}+a_{s})\Gamma(1-\mathbf{b}_{[1]}+a_{k})}{\Gamma(1-\mathbf{a}_{[k,s]}+a_{k})}\sum_{n=0}^{\infty}\frac{(1-a_{s}+a_{k})_{n}}{n!(1-b_{1}+a_{k}+n)}$
$\displaystyle\quad{}\times{}_{p}F_{p-1}\left(\begin{matrix}-n,1-\mathbf{b}_{[1]}+a_{k}\\\
1-\mathbf{a}_{[k]}+a_{k}\end{matrix}\,\vline\,z\right)$ (2.14)
$\displaystyle{}=\frac{z^{a_{k}}\pi(a_{k}-a_{s})\Gamma(1-\mathbf{b}+a_{k})}{\sin(\pi(a_{k}-a_{s}))\Gamma(1-\mathbf{a}_{[k]}+a_{k})}\sum_{n=0}^{\infty}\frac{(1-a_{s}+a_{k})_{n}}{(n+1)!}\,{}_{p+1}F_{p}\left(\begin{matrix}-n,1-\mathbf{b}+a_{k}\\\
1,1-\mathbf{a}_{[k]}+a_{k}\end{matrix}\,\vline\,z\right).$ (2.15)

The series in (2.12) converges in the disk $|z-1|<1$ if
$\Re(1-\mathbf{b}_{[1,2]}+a_{s})>0$, the series in (2.13), (2.14) converge in
the same disk if $\Re(1-\mathbf{b}_{[1]}+a_{s})>0$ and, finally, (2.15)
converges in $|z-1|<1$ if $\Re(1-\mathbf{b}+a_{s})>0$. In all cases, it is
also required that none of the gamma functions in the numerator had poles.
Further, Nørlund found two expansions of his function $y_{k,s}(x)$ in powers
of $1-z$, which we will need below. Written in our notation, formulas [33,
(5.35), (5.36)] read

$\displaystyle G^{2,p}_{p,p}\left(z\,\vline\,\begin{matrix}\mathbf{b}\\\
a_{k},a_{s},\mathbf{a}_{[k,s]}\end{matrix}\right)=z^{a_{s}}\sum_{n=0}^{\infty}D_{n}^{[k,s]}(1-z)^{n},$
(2.16)

where the coefficients $D_{n}^{[k,s]}$ are given by

$\displaystyle
D_{n}^{[k,s]}=\frac{\Gamma(1-\mathbf{b}+a_{k})\Gamma(1-b_{1}+a_{s}+n)\Gamma(1-b_{2}+a_{s}+n)}{\Gamma(1-\mathbf{a}_{[k,s]}+a_{k})\Gamma(2+a_{k}+a_{s}-b_{1}-b_{2}+n)n!}$
$\displaystyle\hphantom{D_{n}^{[k,s]}=}{}\times\sum_{j=0}^{\infty}\frac{(1-b_{1}+a_{k})_{j}(1-b_{2}+a_{k})_{j}}{j!(2+a_{k}+a_{s}-b_{1}-b_{2}+n)_{j}}\,{}_{p-1}F_{p-2}\left(\begin{matrix}-j,1-\mathbf{b}_{[1,2]}+a_{k}\\\
1-\mathbf{a}_{[k,s]}+a_{k}\end{matrix}\right)$ (2.17)
$\displaystyle\hphantom{D_{n}^{[k,s]}}{}=\frac{\Gamma(1-\mathbf{b}+a_{k})\Gamma(1-a_{k}+a_{s}+n)}{\Gamma(1-\mathbf{a}_{[k,s]}+a_{k})n!}\sum_{j=0}^{\infty}\frac{(1-b_{1}+a_{k})_{j}}{(1-b_{1}+a_{s}+n)_{j+1}}$
$\displaystyle\hphantom{D_{n}^{[k,s]}=}{}\times{}_{p}F_{p-1}\left(\begin{matrix}-j,1-\mathbf{b}_{[1]}+a_{k}\\\
1,1-\mathbf{a}_{[k,s]}+a_{k}\end{matrix}\right).$ (2.18)

Here and below ${}_{p}F_{p-1}$ without an argument is understood as
${}_{p}F_{p-1}(1)$. The series (2.16) converges in $|z-1|<1$ if
$\Re(1-\mathbf{b}_{[1,2]}+a_{s})>0$. The same condition suffices for
convergence of (2.17), while (2.18) converges if
$\Re(1-\mathbf{b}_{[1]}+a_{s})>0$. Uniqueness of power series coefficients
implies equality of the coefficients in (2.17) and (2.18). For $p=2$ this
equality is nothing but Gauss formula for ${}_{2}F_{1}(1)$. For $p=3$ we get
after some simplifications and renaming variables

$\displaystyle{}_{3}F_{2}\left(\begin{matrix}\alpha_{1},\alpha_{2},\alpha_{3}\\\
\beta_{1},\beta_{2}\end{matrix}\right)=\frac{\Gamma(\beta_{2}-\alpha_{1}-\alpha_{2}+1)\Gamma(\beta_{2})}{\Gamma(\beta_{2}-\alpha_{1})\Gamma(\beta_{2}-\alpha_{2})}\sum_{j=0}^{\infty}\frac{(\alpha_{1})_{j}}{(\beta_{2}-\alpha_{2})_{j+1}}\,{}_{3}F_{2}\left(\begin{matrix}-j,\alpha_{2},\beta_{1}-\alpha_{3}\\\
\beta_{1},1\end{matrix}\right).$

###  2.4 Connection to Bühring expansion of $\boldsymbol{{}_{p}F_{p-1}}$

In his 1992 paper Bühring found the representation [8, Theorem 2]

$\displaystyle\frac{\Gamma(\boldsymbol{\alpha})}{\Gamma(\boldsymbol{\beta})}{}_{p}F_{p-1}\left(\begin{matrix}\boldsymbol{\alpha}\\\
\boldsymbol{\beta}\end{matrix};z\right)=(1-z)^{\nu}\sum_{n=0}^{\infty}f_{p}(n|\boldsymbol{\alpha},\boldsymbol{\beta})(1-z)^{n}+\sum_{n=0}^{\infty}h_{p}(n|\boldsymbol{\alpha},\boldsymbol{\beta})(1-z)^{n},$
(2.19)

where $\nu=\sum\limits_{k=1}^{p-1}\beta_{k}-\sum\limits_{k=1}^{p}\alpha_{k}$.
He also derived explicit formulas for the coefficients
$f_{p}(n|\boldsymbol{\alpha},\boldsymbol{\beta})$ and
$h_{p}(n|\boldsymbol{\alpha},\boldsymbol{\beta})$ to be given below. On
setting $\boldsymbol{\alpha}=1-\mathbf{b}+a_{s}$,
$\boldsymbol{\beta}=1-\mathbf{a}_{[s]}+a_{s}$, $\nu=\psi_{p}-1$ and denoting

$\displaystyle f_{p}^{s}(n)\equiv
f_{p}(n|1-\mathbf{b}+a_{s},1-\mathbf{a}_{[s]}+a_{s}),$ $\displaystyle
h_{p}^{s}(n)\equiv h_{p}(n|1-\mathbf{b}+a_{s},1-\mathbf{a}_{[s]}+a_{s}),$

Bühring’s formula takes the form

$\displaystyle\frac{\Gamma(1-\mathbf{b}+a_{s})}{\Gamma(1-\mathbf{a}_{[s]}+a_{s})}\,{}_{p}F_{p-1}\left(\begin{matrix}1-\mathbf{b}+a_{s}\\\
1-\mathbf{a}_{[s]}+a_{s}\end{matrix};z\right)$
$\displaystyle\qquad{}=(1-z)^{\psi_{p}-1}\sum_{n=0}^{\infty}f_{p}^{s}(n)(1-z)^{n}+\sum_{n=0}^{\infty}h_{p}^{s}(n)(1-z)^{n}.$
(2.20)

Expansion of the form (2.19) is unique as long as $\nu$ is not an integer.
Hence, in view of (2.4) and (2.17), identity (2.20) is equivalent to formula
(2.3). From (2.3) we conclude that

$\displaystyle(1-z)^{\psi_{p}-1}\sum_{n=0}^{\infty}f_{p}^{s}(n)(1-z)^{n}=\frac{\pi
z^{-a_{s}}}{\sin(\pi\psi_{p})}G^{p,0}_{p,p}\left(z\,\vline\,\begin{matrix}\mathbf{b}\\\
\mathbf{a}\end{matrix}\right)\qquad\text{and}$ (2.21)
$\displaystyle\sum_{n=0}^{\infty}h_{p}^{s}(n)(1-z)^{n}=\frac{-z^{-a_{s}}}{\pi\sin(\pi\psi_{p})}\sum_{\begin{subarray}{c}k=1\\\
k\neq{s}\end{subarray}}^{p}\frac{\sin(\pi(\mathbf{b}-a_{k}))}{\sin(\pi(\mathbf{a}_{[k,s]}-a_{k}))}G^{2,p}_{p,p}\left(z\,\vline\,\begin{matrix}\mathbf{b}-a_{s}\\\
a_{k},a_{s},\mathbf{a}_{[k,s]}\end{matrix}\right)$ (2.22)

for each $s\in\\{1,\ldots,p\\}$. Comparing (2.21) with (2.5) we arrive at

$\displaystyle
f^{s}_{p}(n)=\frac{\Gamma(1-\psi_{p})}{(\psi_{p})_{n}}g^{s}_{p}(n),$ (2.23)

where Euler’s reflection formula has been used. Thus, Bühring coefficients
$f_{p}(n|\boldsymbol{\alpha},\boldsymbol{\beta})$ (denoted by $g_{n}(s)$ in
[8]) are essentially the same as Nørlund coefficients $g^{s}_{p}(n)$ (denoted
by $c^{(s)}_{n,p}$ in [33]). Explicit representation for the coefficients
$f_{p}(n|\boldsymbol{\alpha},\boldsymbol{\beta})$ found in [8, (2.9), (2.16)]
after some renaming of variables and setting $s=p$ can be put into the form

$\displaystyle
f_{p}^{p}(n)=\frac{\Gamma(1-\psi_{p})}{(\psi_{p})_{n}}\sum^{n}_{j_{p-2}=0}\frac{(-1)^{j_{p-2}}(\psi_{p}+a_{p}-b_{p-1}+j_{p-2})_{n-j_{p-2}}}{(n-j_{p-2})!}$
$\displaystyle\hphantom{f_{p}^{p}(n)=}{}\times\sum^{j_{p-2}}_{j_{p-3}=0}\cdots\sum^{j_{3}}_{j_{2}=0}\sum^{j_{2}}_{j_{1}=0}\prod\limits_{m=1}^{p-1}\frac{(\psi_{m}+j_{m-1})_{j_{m}-j_{m-1}}}{(j_{m}-j_{m-1})!}\prod\limits_{m=1}^{p-2}(b_{m}-a_{m+1})_{j_{m}-j_{m-1}}$
$\displaystyle\hphantom{f_{p}^{p}(n)}{}=\frac{\Gamma(1-\psi_{p})}{(\psi_{p})_{n}}\sum^{\infty}_{j_{1},j_{2},\ldots,j_{p-2}=0}\frac{(\psi_{p}+a_{p}-b_{p-1})_{n}}{(\psi_{p}+a_{p}-b_{p-1})_{j_{p-2}}}\prod\limits_{m=1}^{p-1}\frac{(\psi_{m})_{j_{m}}(-j_{m})_{j_{m-1}}}{(\psi_{m})_{j_{m-1}}j_{m}!}$
$\displaystyle\hphantom{f_{p}^{p}(n)=}{}\times\prod\limits_{m=1}^{p-2}\frac{(b_{m}-a_{m+1})_{j_{m}}}{(1-b_{m}+a_{m+1}-j_{m})_{j_{m-1}}},$

where $\psi_{m}=\sum\limits_{i=1}^{m}(b_{i}-a_{i})$ and $j_{0}=0$,
$j_{p-1}=n$. In view of (2.9), formula (2.23) then leads to the following
(presumably new) transformation for multiple hypergeometric series

$\displaystyle\sum^{n}_{j_{p-2}=0}\frac{(\psi_{p}+a_{p}-b_{p-1}+j_{p-2})_{n-j_{p-2}}}{(-1)^{j_{p-2}}(n-j_{p-2})!}$
$\displaystyle\qquad\quad{}\times\sum^{j_{p-2}}_{j_{p-3}=0}\cdots\sum^{j_{3}}_{j_{2}=0}\sum^{j_{2}}_{j_{1}=0}\prod\limits_{m=1}^{p-1}\frac{(\psi_{m}+j_{m-1})_{j_{m}-j_{m-1}}}{(j_{m}-j_{m-1})!}\prod\limits_{m=1}^{p-2}(b_{m}-a_{m+1})_{j_{m}-j_{m-1}}$
$\displaystyle\qquad{}=\sum^{n}_{j_{p-2}=0}\cdots\sum^{j_{3}}_{j_{2}=0}\sum^{j_{2}}_{j_{1}=0}\prod\limits_{m=1}^{p-1}\frac{(\psi_{m}+j_{m-1})_{j_{m}-j_{m-1}}}{(j_{m}-j_{m-1})!}(b_{m+1}-a_{m})_{j_{m}-j_{m-1}},$

where $j_{0}=0$ and $j_{p-1}=n$. For $p=3$ this formula reduces to the
identity

$\displaystyle{}_{3}F_{2}\left(\begin{matrix}-n,\alpha_{1},\alpha_{2}\\\
\beta_{1},\beta_{2}\end{matrix}\right)=\frac{(\beta_{2}-\alpha_{2})_{n}}{(\beta_{2})_{n}}\,{}_{3}F_{2}\left(\begin{matrix}-n,\beta_{1}-\alpha_{1},\alpha_{2}\\\
\beta_{1},1-\beta_{2}+\alpha_{2}-n\end{matrix}\right),$

which is a guise of Sheppard’s transformation [2, Corollary 3.3.4], also
rediscovered by Bühring [9, (4.1)]. For $p=4$ the corresponding result is
given by [9, (4.2)]

$\displaystyle\sum\limits_{k=0}^{n}\frac{(-k)_{n}(\alpha_{1})_{k}(\alpha_{2})_{k}}{(\beta_{1})_{k}(\beta_{2})_{k}k!}\,{}_{3}F_{2}\left(\begin{matrix}-k,\gamma_{1},\gamma_{2}\\\
\alpha_{1},\alpha_{2}\end{matrix}\right)$
$\displaystyle\qquad{}=\frac{(\beta_{2}-\alpha_{2})_{n}}{(\beta_{2})_{n}}\sum\limits_{k=0}^{n}\frac{(-k)_{n}(\alpha_{2})_{k}(\beta_{1}-\alpha_{1})_{k}}{(\beta_{1})_{k}(1+\alpha_{2}-\beta_{2}-n)_{k}k!}\,{}_{3}F_{2}\left(\begin{matrix}-k,\gamma_{1},\gamma_{2}\\\
\alpha_{2},1+\alpha_{1}-\beta_{1}-k\end{matrix}\right).$

In his paper [8] Bühring mentioned that the structure of (2.19) “was already
given in the classic paper by Nørlund [33], but the coefficients were not all
known”. In his subsequent joint work [10] with Srivastava is it said that the
coefficients $f_{p}^{s}(n)$ can be computed using Nørlund’s recurrence, “but
it is desirable to get an explicit representation” which is indeed derived in
[10] in terms of some limit relations. Explicit formulas for the coefficients
$h_{p}^{s}(n)$ have been first found in [8]. Adopted to our notation Bühring’s
formulas [8, (2.9), (2.16)], [10, (2.13), (2.15)] read

$\displaystyle
h_{p}^{p}(n)=-\frac{\Gamma(\psi_{p})\Gamma(1-b_{p}+a_{p}+n)\Gamma(1-b_{p-1}+a_{p}+n)}{(1-\psi_{p})_{n+1}\Gamma(\psi_{p-1})\Gamma(\psi_{p}-b_{p-1}+a_{p})n!}$
$\displaystyle\hphantom{h_{p}^{p}(n)=}{}\times\sum\limits_{j_{p-2}=0}^{\infty}\frac{(\psi_{p}-n-1)_{j_{p-2}}}{(\psi_{p-1})_{j_{p-2}}(\psi_{p}-b_{p-1}+a_{p})_{j_{p-2}}}$
$\displaystyle\hphantom{h_{p}^{p}(n)=}{}\times\sum\limits_{j_{p-3}=0}^{j_{p-2}}\cdots\sum\limits_{j_{1}=0}^{j_{2}}\prod\limits_{m=1}^{p-2}\frac{(\psi_{m}+j_{m-1})_{j_{m}-j_{m-1}}}{(j_{m}-j_{m-1})!}(b_{m}-a_{m+1})_{j_{m}-j_{m-1}},$
(2.24)

where, as before, $\psi_{m}=\sum\limits_{i=1}^{m}(b_{i}-a_{i})$ and $j_{0}=0$.
Conditions for convergence for the outer series above have also been found in
[8] and [10] and are given by

$\displaystyle\Re(1-b_{i}+a_{p}+n)>0\qquad\text{for}\quad i=1,2,\ldots,p-2.$

Bühring and Srivastava found expressions for the coefficients
$f_{p}(n|\boldsymbol{\alpha},\boldsymbol{\beta})$,
$h_{p}(n|\boldsymbol{\alpha},\boldsymbol{\beta})$ as limits of hypergeometric
polynomials [10, (3.13)]. They computed this limit for $p=3,4,5$ yielding [10,
(3.15), (3.17)]

$\displaystyle h^{s}_{3}(n)\equiv
h_{3}(n|1-\mathbf{b}+a_{s},1-\mathbf{a}_{[s]}+a_{s})$
$\displaystyle\hphantom{h^{s}_{3}(n)}{}=-\frac{\Gamma(\psi_{3})\Gamma(1-b_{2}+a_{s}+n)\Gamma(1-b_{3}+a_{s}+n)}{(1-\psi_{3})_{n+1}n!\Gamma(\psi_{3}-b_{2}+a_{s})\Gamma(\psi_{3}-b_{3}+a_{s})}\,{}_{3}F_{2}\left(\begin{matrix}\psi_{3}-1-n,b_{1}-\mathbf{a}_{[s]}\\\
\psi_{3}-\mathbf{b}_{[1]}+a_{s}\end{matrix}\right),\\!\\!\\!\\!$ (2.25)
$\displaystyle
h^{s}_{4}(n)=-\frac{\Gamma(\psi_{4})\Gamma(1-b_{3}+a_{s}+n)\Gamma(1-b_{4}+a_{s}+n)}{(1-\psi_{4})_{n+1}n!\Gamma(\psi_{4}-b_{3}+a_{s})\Gamma(\psi_{4}-b_{4}+a_{s})}$
$\displaystyle\hphantom{h^{s}_{4}(n)=}{}\times\sum\limits_{k=0}^{\infty}(\psi_{4}-1-n)_{k}\frac{(b_{1}+b_{2}-a_{i_{1}}-a_{i_{2}})_{k}(b_{1}+b_{2}-a_{i_{1}}-a_{i_{3}})_{k}}{(\psi_{4}-b_{3}+a_{s})_{k}(\psi_{4}-b_{4}+a_{s})_{k}}$
$\displaystyle\hphantom{h^{s}_{4}(n)=}{}\times{}_{3}F_{2}\left(\begin{matrix}-k,b_{1}-a_{i_{1}},b_{2}-a_{i_{1}}\\\
b_{1}+b_{2}-a_{i_{1}}-a_{i_{2}},b_{1}+b_{2}-a_{i_{1}}-a_{i_{3}}\end{matrix}\right),$

where $\\{a_{i_{1}},a_{i_{2}},a_{i_{3}}\\}=\mathbf{a}_{[s]}$ (i.e., if $s=1$
then $\\{a_{i_{1}},a_{i_{2}},a_{i_{3}}\\}=\\{a_{2},a_{3},a_{4}\\}$, if $s=2$
then $\\{a_{i_{1}},a_{i_{2}},a_{i_{3}}\\}=\\{a_{1},a_{3},a_{4}\\}$, etc.).
Expression for $p=5$ is quite cumbersome and can be found in [10, (3.19)]. In
spite of their non-symmetric appearance, both formulas are invariant with
respect to permutation of the elements of $\mathbf{b}$.

Let us now cite the formulas for the coefficients $g^{s}_{p}(n)$ for
$p=2,3,4$. For $p=2$ the corresponding formulas can be read off formula (1.2).
For both $p=2$ and $p=3$ expressions for $g^{s}_{p}(n)$ have been found by
Nørlund, see [33, (2.10)]. They are

$\displaystyle
g^{s}_{2}(n)=\frac{(b_{1}-a_{3-s})_{n}(b_{2}-a_{3-s})_{n}}{n!},\qquad s=1,2,$

and

$\displaystyle
g^{s}_{3}(n)=\frac{(\psi_{3}-b_{2}+a_{s})_{n}(\psi_{3}-b_{3}+a_{s})_{n}}{n!}\,{}_{3}F_{2}\left(\begin{matrix}-n,b_{1}-a_{i_{1}},b_{1}-a_{i_{2}}\\\
\psi_{3}-b_{2}+a_{s},\psi_{3}-b_{3}+a_{s}\end{matrix}\right),\qquad s=1,2,3,$

where $\\{a_{i_{1}},a_{i_{2}}\\}=\mathbf{a}_{[s]}$. Notwithstanding the non-
symmetric appearance, the last formula is symmetric with respect to the
elements of $\mathbf{b}$. For $p=4$ we convert the expression for
$f_{p}^{s}(n)$ calculated in [10, (3.17)] into expressions for $g^{s}_{p}(n)$
using (2.23) and the necessary renaming of parameters. This yields

$\displaystyle
g_{4}^{s}(n)=\frac{(\psi_{4}-b_{3}+a_{s})_{n}(\psi_{4}-b_{4}+a_{s})_{n}}{n!}$
$\displaystyle\hphantom{g_{4}^{s}(n)=}{}\times\sum\limits_{k=0}^{n}\frac{(-n)_{k}(b_{1}+b_{2}-a_{i_{1}}-a_{i_{2}})_{k}(b_{1}+b_{2}-a_{i_{1}}-a_{i_{3}})_{k}}{(\psi_{4}-b_{3}+a_{s})_{k}(\psi_{4}-b_{4}+a_{s})_{k}}$
$\displaystyle\hphantom{g_{4}^{s}(n)=}{}\times{}_{3}F_{2}\left(\begin{matrix}-k,b_{1}-a_{i_{1}},b_{2}-a_{i_{1}}\\\
b_{1}+b_{2}-a_{i_{1}}-a_{i_{2}},b_{1}+b_{2}-a_{i_{1}}-a_{i_{3}}\end{matrix}\right),$

where $\\{a_{i_{1}},a_{i_{2}},a_{i_{3}}\\}=\mathbf{a}_{[s]}$. This formula is
also invariant with respect to permutation of the elements of $\mathbf{b}$.

##  3 Main results

Having made these preparations, we are ready to formulate our main results.
First, we give explicit formulas for the Nørlund’s coefficients
$g^{s}_{p}(1)$, $g^{s}_{p}(2)$, $g^{s}_{p}(3)$ with arbitrary $p$. Note that
formula (2.9) contains $p-2$ summations even for small $n$. Combined with
(2.5) and $g^{p}_{p}(0)=1$ our formulas essentially provide first four terms
in asymptotic expansion of $G^{p,0}_{p,p}$ as $z\to 1$. Then we derive a new
way to compute the coefficients $g^{s}_{p}(n)$ for all $n$ in terms of
generalized Bernoulli polynomials in Theorem 3.3. Theorems 3.5 and 3.10
contain new identities for the coefficients $g^{s}_{p}(n)$, $h_{p}^{s}(n)$ and
$D_{n}^{[k,s]}$ defined by expansions (2.5), (2.16) and (2.20), respectively.
For small values of $n$ these identities lead to new and known relations
involving hypergeometric and trigonometric functions, in particular, to
identity (1.3).

######  Theorem 3.1.

The first four coefficients of the expansion (2.5) are given by

$\displaystyle g^{p}_{p}(0)=1,\qquad
g^{p}_{p}(1)=\sum_{m=1}^{p-1}(b_{m+1}-a_{m})\psi_{m},$ (3.1) $\displaystyle
g^{p}_{p}(2)=\frac{1}{2}\sum_{m=1}^{p-1}(b_{m+1}-a_{m})_{2}(\psi_{m})_{2}+\sum_{k=2}^{p-1}(b_{k+1}-a_{k})(\psi_{k}+1)\sum_{m=1}^{k-1}(b_{m+1}-a_{m})\psi_{m}$
(3.2)

and

$\displaystyle
g^{p}_{p}(3)=\frac{1}{6}\sum_{m=1}^{p-1}(b_{m+1}-a_{m})_{3}(\psi_{m})_{3}+\frac{1}{2}\sum_{k=2}^{p-1}(b_{k+1}-a_{k})(\psi_{k}+2)\sum_{m=1}^{k-1}(b_{m+1}-a_{m})_{2}(\psi_{m})_{2}$
$\displaystyle\hphantom{g^{p}_{p}(3)=}{}+\frac{1}{2}\sum_{k=2}^{p-1}(\psi_{k}+1)_{2}(b_{k+1}-a_{k})_{2}\sum_{m=1}^{k-1}(b_{m+1}-a_{m})\psi_{m}$
$\displaystyle\hphantom{g^{p}_{p}(3)=}{}+\sum_{n=3}^{p-1}(b_{n+1}-a_{n})(\psi_{n}+2)\sum_{k=2}^{n-1}(b_{k+1}-a_{k})(\psi_{k}+1)\sum_{m=1}^{k-1}(b_{m+1}-a_{m})\psi_{m},$
(3.3)

where, as before, $\psi_{m}=\sum\limits_{j=1}^{m}(b_{j}-a_{j})$.

###### Proof.

We note that for $n=1$ summation in (2.9) is over the index sets of the form

$\displaystyle\\{j_{0},j_{1},\ldots_{,}j_{p-1}\\}=\\{0,\ldots,0,1,\ldots,1\\},$

where the number of ones changes from $1$ to $p-1$. In view of this
observation rearrangement of the formula (2.9) yields (3.1). Next, we prove
(3.3). Summation in (2.9) is over all Young diagrams that fit $3\times(p-2)$
box. We break all possible diagrams in four disjoint groups as follows (by
definition $j_{0}=0$, $j_{p-1}=3$):

$\displaystyle(1)\quad$ $\displaystyle j_{0}=j_{1}=\cdots={j_{m-1}}=0,\qquad
j_{m}=\cdots=j_{p-1}=3,\qquad m\in\\{1,2,\ldots,p-1\\},$
$\displaystyle(2)\quad$ $\displaystyle j_{0}=j_{1}=\cdots={j_{m-1}}=0,\qquad
j_{m}=\cdots=j_{k-1}=2,\qquad j_{k}=\cdots=j_{p-1}=3,$ $\displaystyle
k\in\\{2,\ldots,p-1\\},\qquad m\in\\{1,\ldots,k-1\\},\qquad m<k,$
$\displaystyle(3)\quad$ $\displaystyle j_{0}=j_{1}=\cdots={j_{m-1}}=0,\qquad
j_{m}=\cdots=j_{k-1}=1,\qquad j_{k}=\cdots=j_{p-1}=3,$ $\displaystyle
k\in\\{2,\ldots,p-1\\},\qquad m\in\\{1,\ldots,k-1\\},\qquad m<k,$
$\displaystyle(4)\quad$ $\displaystyle j_{0}=\cdots={j_{m-1}}=0,\qquad
j_{m}=\cdots=j_{k-1}=1,\qquad j_{k}=\cdots=j_{n-1}=2,$ $\displaystyle
j_{n}=\cdots=j_{p-1}=3,\qquad n\in\\{3,\ldots,p-1\\},\qquad
k\in\\{2,\ldots,n-1\\},$ $\displaystyle m\in\\{1,\ldots,k-1\\},\qquad m<k<n.$

Summation over the first type of diagrams leads to the first term (3.3).
Similarly, summation over the $i$-th type of diagrams leads to the $i$-th term
(3.3) for $i=2,3,4$. Analogous considerations lead to (3.2). ∎

######  Remark 3.2.

Exchanging the roles of $a_{p}$ and $a_{s}$ formulas (3.1), (3.2) and (3.3)
lead to expressions for $g^{s}_{p}(1)$, $g^{s}_{p}(2)$, $g^{s}_{p}(3)$. Since
each coefficient $g^{s}_{p}(n)$ is a symmetric polynomial in the elements of
$\mathbf{a}_{[s]}$ and $\mathbf{b}$, such expressions can be expanded in terms
of some basis of symmetric polynomials. In particular, denoting by
$e_{k}(x_{1},\ldots,x_{m})$ the $k$-th elementary symmetric polynomial of
$x_{1},\ldots,x_{m}$, we get

$\displaystyle
g^{s}_{p}(1)=e_{2}(\mathbf{b})-e_{2}(\mathbf{a}_{[s]})+e_{1}(\mathbf{a}_{[s]})(e_{1}(\mathbf{a}_{[s]})-e_{1}(\mathbf{b}))$

and (with the help of Mathematica)

$\displaystyle
g^{s}_{p}(2)=\frac{1}{2}e_{1}(\mathbf{a}_{[s]})^{4}-e_{1}(\mathbf{a}_{[s]})^{3}+\frac{1}{2}e_{1}(\mathbf{a}_{[s]})^{2}+\frac{3}{2}e_{1}(\mathbf{a}_{[s]})e_{2}(\mathbf{a}_{[s]})+\frac{1}{2}e_{2}(\mathbf{a}_{[s]})^{2}-\frac{1}{2}e_{2}(\mathbf{a}_{[s]})$
$\displaystyle\hphantom{g^{s}_{p}(2)=}{}-e_{1}(\mathbf{a}_{[s]})^{2}e_{2}(\mathbf{a}_{[s]})-\frac{1}{2}e_{3}(\mathbf{a}_{[s]})+\frac{1}{2}e_{2}(\mathbf{b})^{2}+\frac{1}{2}e_{2}(\mathbf{b})+\frac{1}{2}e_{3}(\mathbf{b})-\frac{1}{2}e_{1}(\mathbf{a}_{[s]})e_{1}(\mathbf{b})$
$\displaystyle\hphantom{g^{s}_{p}(2)=}{}+\frac{3}{2}e_{1}(\mathbf{a}_{[s]})^{2}e_{1}(\mathbf{b})-e_{1}(\mathbf{a}_{[s]})^{3}e_{1}(\mathbf{b})-\frac{1}{2}e_{1}(\mathbf{a}_{[s]})e_{1}(\mathbf{b})^{2}+\frac{1}{2}e_{1}(\mathbf{a}_{[s]})^{2}e_{1}(\mathbf{b})^{2}$
$\displaystyle\hphantom{g^{s}_{p}(2)=}{}+e_{1}(\mathbf{a}_{[s]})^{2}e_{2}(\mathbf{b})+\frac{1}{2}e_{1}(\mathbf{b})e_{2}(\mathbf{b})-e_{1}(\mathbf{a}_{[s]})e_{1}(\mathbf{b})e_{2}(\mathbf{b})-e_{2}(\mathbf{a}_{[s]})e_{2}(\mathbf{b})-e_{1}(\mathbf{b})e_{2}(\mathbf{a}_{[s]})$
$\displaystyle\hphantom{g^{s}_{p}(2)=}{}+e_{1}(\mathbf{a}_{[s]})e_{2}(\mathbf{a}_{[s]})e_{1}(\mathbf{b})-e_{1}(\mathbf{a}_{[s]})e_{2}(\mathbf{b}).$

Next theorem gives a presumably new expression for Nørlund’s coefficients
$g_{p}^{k}(n)$ in terms of the generalized Bernoulli polynomials. Let us start
by recalling that the Bernoulli–Nørlund (or the generalized Bernoulli)
polynomial $\mathcal{B}^{(\sigma)}_{k}(x)$ is defined by the generating
function [34, (1)]:

$\displaystyle\frac{t^{\sigma}e^{xt}}{(e^{t}-1)^{\sigma}}=\sum\limits_{k=0}^{\infty}\mathcal{B}^{(\sigma)}_{k}(x)\frac{t^{k}}{k!}.$

In particular, $\mathcal{B}^{(1)}_{k}(x)=\mathcal{B}_{k}(x)$ is the classical
Bernoulli polynomial.

######  Theorem 3.3.

Coefficients $g_{p}^{k}(n)$, defined in (2.5), are given by any of the
following formulas

$\displaystyle
g_{p}^{k}(n)=\sum\limits_{r=0}^{n}\frac{(-1)^{n-r}(r+1)_{n-r}}{(n-r)!}\tilde{l}_{r}\mathcal{B}^{(n+1)}_{n-r}(2-a_{k}-\psi_{p})$
$\displaystyle\hphantom{g_{p}^{k}(n)}{}=\sum\limits_{r=0}^{n}\frac{(-1)^{n-r}(\psi_{p}+r)_{n-r}}{(n-r)!}l_{r}\mathcal{B}^{(n+\psi_{p})}_{n-r}(1-a_{k}).$
(3.4)

Here $\tilde{l}_{0}=1$ and $\tilde{l}_{r}$, $r\geq 1$, are found from the
recurrence

$\displaystyle\tilde{l}_{r}=\frac{1}{r}\sum\limits_{m=1}^{r}\widetilde{q}_{m}\tilde{l}_{r-m},$
(3.5)

where

$\displaystyle\widetilde{q}_{m}=\frac{(-1)^{m+1}}{m+1}\left[\mathcal{B}_{m+1}(a_{k}+\psi_{p}-1)-{\mathcal{B}}_{m+1}(a_{k})+\sum\limits_{j=1}^{p}{(\mathcal{B}_{m+1}(a_{j})}-{\mathcal{B}_{m+1}(b_{j}))}\right].$

Similarly, the coefficients $l_{r}$ satisfy the recurrence relation

$\displaystyle
l_{r}=\frac{1}{r}\sum\limits_{m=1}^{r}q_{m}l_{r-m}\qquad\text{with}\quad
l_{0}=1$ (3.6)

and

$\displaystyle
q_{m}=\frac{(-1)^{m+1}}{m+1}\sum\limits_{j=1}^{p}\big{(}\mathcal{B}_{m+1}(a_{j})-\mathcal{B}_{m+1}(b_{j})\big{)}.$

######  Remark 3.4.

It is known [22, Lemma 1] that the recurrences (3.5) and (3.6) can be solved
to give the following explicit expressions for $l_{r}$:

$\displaystyle
l_{r}=\sum\limits_{k_{1}+2k_{2}+\cdots+rk_{r}=r}\frac{q_{1}^{k_{1}}(q_{2}/2)^{k_{2}}\cdots(q_{r}/r)^{k_{r}}}{k_{1}!k_{2}!\cdots
k_{r}!}=\sum\limits_{n=1}^{r}\frac{1}{n!}\sum\limits_{k_{1}+k_{2}+\cdots+k_{n}=r}\prod\limits_{i=1}^{n}\frac{q_{k_{i}}}{k_{i}}.$

Similar formula is of course true for $\tilde{l}_{r}$ once we write
$\widetilde{q}_{m}$ instead of $q_{m}$. Moreover, Nair [32, Section 8] found a
determinantal expression for such solution which in our notation takes the
form

$\displaystyle
l_{r}=\frac{\det(\Omega_{r})}{r!},\qquad\Omega_{r}=[\omega_{i,j}]_{i,j=1}^{r},\qquad\omega_{i,j}=\begin{cases}q_{i-j+1}(i-1)!/(j-1)!,&i\geq{j},\\\
-1,&i=j-1,\\\ 0,&i<j-1.\end{cases}$

###### Proof.

The theorem is a corollary of an expansion of the $H$-function of Fox found in
our recent paper [24]. Since Meijer’s $G$-function is a particular case of
Fox’s $H$-function, formula (3.4) is a particular case of [24, Theorem 1] once
we set $p=q$, $\mathbf{A}=\mathbf{B}=(1,\ldots,1)$, $\nu=1$, $\mu=\psi_{p}$,
$\theta=a_{k}-1$ in that theorem. ∎

As before, we use the shorthand notation

$\displaystyle\sin(\pi(\mathbf{b}-a_{k}))=\prod_{j=1}^{p}\sin(\pi(b_{j}-a_{k})),\qquad\sin(\pi(\mathbf{a}_{[k]}-a_{k}))=\prod_{\begin{subarray}{c}j=1\\\
j\neq{k}\end{subarray}}^{p}\sin(\pi(a_{j}-a_{k})),$

for the products and $[a]_{j}=a(a-1)\cdots(a-j+1)$ for the falling factorial.

######  Theorem 3.5.

For each nonnegative integer $m$ the following identities holds

$\displaystyle\sum\limits_{j=0}^{m}\frac{(-1)^{j}}{j!}\sum\limits_{k=1}^{p}[a_{k}]_{j}h^{k}_{p}(m-j)\frac{\sin(\pi(\mathbf{b}-a_{k}))}{\sin(\pi(\mathbf{a}_{[k]}-a_{k}))}=0,$
(3.7)

where the numbers $h^{p}_{p}(n)$ are defined by (2.20) and given explicitly by
(2.24), $h^{k}_{p}(n)$ is obtained from $h^{p}_{p}(n)$ by exchanging the roles
of $a_{p}$ and $a_{k}$; furthermore,

$\displaystyle\sum\limits_{j=0}^{m}\frac{(-1)^{j}}{(\psi_{p})_{m-j}j!}\left\\{[a_{s}]_{j}g^{s}_{p}(m-j)\sin(\pi\psi_{p})-\sum\limits_{k=1}^{p}[a_{k}]_{j}g^{k}_{p}(m-j)\frac{\sin(\pi(\mathbf{b}-a_{k}))}{\sin(\pi(\mathbf{a}_{[k]}-a_{k}))}\right\\}=0,$
(3.8)

where $s\in\\{1,\ldots,p\\}$ is chosen arbitrarily and the numbers
$g^{k}_{p}(n)$ are defined by expansion (2.5) and solve the recurrence (2.7)
in $n$ and the recurrence (2.8) in $p$. They are given explicitly by (2.9).

###### Proof.

Assume that the components of the vector $\mathbf{a}$ are distinct modulo $1$.
Substituting expansion (2.20) into formula (2.4) and taking account of (2.23)
we obtain

$\displaystyle G^{p,0}_{p,p}\left(z\,\vline\,\begin{matrix}\mathbf{b}\\\
\mathbf{a}\end{matrix}\right)=(1-z)^{\psi_{p}-1}\sum\limits_{k=1}^{p}z^{a_{k}}\frac{\Gamma(\mathbf{a}_{[k]}-a_{k})\Gamma(1-\mathbf{a}_{[k]}+a_{k})}{\Gamma(\mathbf{b}-a_{k})\Gamma(1-\mathbf{b}+a_{k})}\sum_{n=0}^{\infty}\frac{\Gamma(1-\psi_{p})}{(\psi_{p})_{n}}g^{k}_{p}(n)(1-z)^{n}$
$\displaystyle\hphantom{G^{p,0}_{p,p}\left(z\,\vline\,\begin{matrix}\mathbf{b}\\\
\mathbf{a}\end{matrix}\right)=}{}+\sum\limits_{k=1}^{p}z^{a_{k}}\frac{\Gamma(\mathbf{a}_{[k]}-a_{k})\Gamma(1-\mathbf{a}_{[k]}+a_{k})}{\Gamma(\mathbf{b}-a_{k})\Gamma(1-\mathbf{b}+a_{k})}\sum_{n=0}^{\infty}h_{p}^{k}(n)(1-z)^{n}$
$\displaystyle\hphantom{G^{p,0}_{p,p}\left(z\,\vline\,\begin{matrix}\mathbf{b}\\\
\mathbf{a}\end{matrix}\right)}{}=\frac{1}{\pi}(1-z)^{\psi_{p}-1}\Gamma(1-\psi_{p})\sum_{n=0}^{\infty}\frac{(1-z)^{n}}{(\psi_{p})_{n}}\sum\limits_{k=1}^{p}z^{a_{k}}\frac{\sin(\pi(\mathbf{b}-a_{k}))}{\sin(\pi(\mathbf{a}_{[k]}-a_{k}))}g^{k}_{p}(n)$
$\displaystyle\hphantom{G^{p,0}_{p,p}\left(z\,\vline\,\begin{matrix}\mathbf{b}\\\
\mathbf{a}\end{matrix}\right)=}{}+\frac{1}{\pi}\sum_{n=0}^{\infty}(1-z)^{n}\sum\limits_{k=1}^{p}z^{a_{k}}\frac{\sin(\pi(\mathbf{b}-a_{k}))}{\sin(\pi(\mathbf{a}_{[k]}-a_{k}))}h_{p}^{k}(n),$

where we applied Euler’s reflection formula
$\Gamma(z)\Gamma(1-z)=\pi/\sin(\pi{z})$. Further, substitute Nørlund’s
expansion (2.5) in place of $G$-function on the left-hand side and rearrange
terms to get

$\displaystyle(1-z)^{\psi_{p}-1}\sum\limits_{n=0}^{\infty}(1-z)^{n}\left\\{\frac{z^{a_{s}}g^{s}_{p}(n)}{\Gamma(\psi_{p}+n)}-\frac{\Gamma(1-\psi_{p})}{(\psi_{p})_{n}}\sum\limits_{k=1}^{p}\frac{z^{a_{k}}\sin(\pi(\mathbf{b}-a_{k}))}{\pi\sin(\pi(\mathbf{a}_{[k]}-a_{k}))}g^{k}_{p}(n)\right\\}$
$\displaystyle\qquad{}=\sum_{n=0}^{\infty}(1-z)^{n}\sum\limits_{k=1}^{p}\frac{z^{a_{k}}\sin(\pi(\mathbf{b}-a_{k}))}{\pi\sin(\pi(\mathbf{a}_{[k]}-a_{k}))}h_{p}^{k}(n),$

where $s\in\\{1,\ldots,p\\}$ can be chosen arbitrarily. Denote for brevity

$\displaystyle\gamma_{n,s}(z)=\frac{z^{a_{s}}g^{s}_{p}(n)}{\Gamma(\psi_{p}+n)}-\frac{\Gamma(1-\psi_{p})}{(\psi_{p})_{n}}\sum\limits_{k=1}^{p}\frac{z^{a_{k}}\sin(\pi(\mathbf{b}-a_{k}))}{\pi\sin(\pi(\mathbf{a}_{[k]}-a_{k}))}g^{k}_{p}(n),$
$\displaystyle\chi_{n}(z)=\sum\limits_{k=1}^{p}\frac{z^{a_{k}}\sin(\pi(\mathbf{b}-a_{k}))}{\pi\sin(\pi(\mathbf{a}_{[k]}-a_{k}))}h_{p}^{k}(n),$

so that the above equality reduces to

$\displaystyle(1-z)^{\psi_{p}-1}\sum\limits_{n=0}^{\infty}(1-z)^{n}\gamma_{n,s}(z)=\sum_{n=0}^{\infty}(1-z)^{n}\chi_{n}(z).$
(3.9)

Assume for a moment that $\psi_{p}$ is not an integer. We know that both
series in (3.9) converge in $|z-1|<1$. Furthermore, all functions
$\gamma_{n,s}(z)$ and $\chi_{n}(z)$ are analytic the same disk. This implies
that (3.9) is only possible for all $z$ in a disk centered at $1$ if

$\displaystyle\sum\limits_{n=0}^{\infty}(1-z)^{n}\gamma_{n,s}(z)\equiv
0\qquad\text{and}\qquad\sum_{n=0}^{\infty}(1-z)^{n}\chi_{n}(z)\equiv 0.$
(3.10)

In terms of the functions $\chi_{n}(z)$ the claimed identity (3.7) takes the
form

$\displaystyle\sum\limits_{j=0}^{m}\frac{(-1)^{j}}{j!}\chi_{m-j}^{(j)}(1)=0,$
(3.11)

which we prove by induction in $m$.

Letting $z\to 1$ in the second identity in (3.10) we get $\chi_{0}(1)=0$ which
establishes our claim for $m=0$. Next, suppose (3.11) holds for
$m=0,1,2,\ldots,r-1$. Divide the second identity in (3.10) by $(1-z)^{r}$ and
expand each $\chi_{n}(z)$, $n=0,1,\ldots,r$, in Taylor series around $z=1$:

$\displaystyle
0=\sum_{n=0}^{r}(1-z)^{n-r}\bigl{\\{}\chi_{n}(1)+\chi_{n}^{\prime}(1)(z-1)+\cdots$
$\displaystyle\hphantom{0=}{}+\chi_{n}^{(r-n)}(z-1)^{r-n}/(r-n)!+O\big{(}(z-1)^{r-n+1}\big{)}\bigr{\\}}+\sum_{n=r+1}^{\infty}(1-z)^{n-r}\chi_{n}(z)$
$\displaystyle\hphantom{0}{}=\sum\limits_{k=0}^{r}(1-z)^{-k}\sum\limits_{i=0}^{r-k}\frac{(-1)^{i}}{i!}\chi_{r-k-i}^{(i)}(1)+O(z-1)=\sum\limits_{i=0}^{r}\frac{(-1)^{i}}{i!}\chi_{r-k-i}^{(i)}(1)+O(z-1),$

where the last equality is by induction hypothesis. We now obtain (3.11) for
$m=r$ on letting $z\to 1$ in this formula. It is immediate to check that the
claimed identity (3.8) takes the form

$\displaystyle\sum\limits_{j=0}^{m}\frac{(-1)^{j}}{j!}\gamma_{m-j,s}^{(j)}(1)=0,$

which can be demonstrated in a similar fashion starting with the first
identity in (3.10).

Finally, we remove the assumption that $\psi_{p}$ is not an integer. The left-
hand sides of (3.7) and (3.8) are analytic functions of, say, parameter
$b_{1}$ except for possible poles. Identities (3.7) and (3.8) for non-integer
$\psi_{p}$ then clearly imply by analytic continuation that these poles are
removable and both identities hold for all $\psi_{p}$. ∎

######  Corollary 3.6.

Suppose for some $i\in\\{1,2,3\\}$ inequality $\Re(b_{i})<\Re(a_{k}+1)$ holds
for $k\in\\{1,2,3\\}$. Then

$\displaystyle\sum\limits_{k=1}^{3}\frac{\sin(\pi(\mathbf{b}-a_{k}))\Gamma(1-\mathbf{b}_{[i]}+a_{k})}{\sin(\pi(\mathbf{a}_{[k]}-a_{k}))\Gamma(\psi-\mathbf{b}_{[i]}+a_{k})}\,{}_{3}F_{2}\left(\begin{matrix}\psi-1,b_{i}-\mathbf{a}_{[k]}\\\
\psi-\mathbf{b}_{[i]}+a_{k}\end{matrix}\right)=0$ (3.12)

and

$\displaystyle\sum\limits_{k=1}^{3}\frac{\sin(\pi(\mathbf{b}-a_{k}))\Gamma(1-\mathbf{b}_{[i]}+a_{k})}{\sin(\pi(\mathbf{a}_{[k]}-a_{k}))\Gamma(\psi-\mathbf{b}_{[i]}+a_{k})}\left\\{(1-\mathbf{b}_{[i]}+a_{k})\,{}_{3}F_{2}\left(\begin{matrix}\psi-2,b_{i}-\mathbf{a}_{[k]}\\\
\psi-\mathbf{b}_{[i]}+a_{k}\end{matrix}\right)\right.$
$\displaystyle\left.\qquad{}-a_{k}(2-\psi)\,{}_{3}F_{2}\left(\begin{matrix}\psi-1,b_{i}-\mathbf{a}_{[k]}\\\
\psi-\mathbf{b}_{[i]}+a_{k}\end{matrix}\right)\right\\}=0,$ (3.13)

where $\psi=\sum\limits_{k=1}^{3}(b_{k}-a_{k})$ and ${}_{p}F_{p-1}$ without
argument is understood as ${}_{p}F_{p-1}(1)$.

###### Proof.

Put $p=3$. Identities (3.12) and (3.13) are now reformulations of (3.7) for
$m=0$ and $m=1$, respectively. The coefficients $h_{3}^{k}(0)$ and
$h_{3}^{k}(1)$ have been computed by formula (2.25). ∎

The next corollary is a rewriting of (3.8) for $m=0$ in view of
$g^{k}_{p}(0)=1$.

######  Corollary 3.7.

For any complex vectors $\mathbf{a}$, $\mathbf{b}$ the following identity
holds

$\displaystyle\sum\limits_{k=1}^{p}\frac{\sin(\pi(\mathbf{b}-a_{k}))}{\sin(\pi(\mathbf{a}_{[k]}-a_{k}))}=\sin(\pi\psi_{p}).$
(3.14)

The right-hand side gives a continuous extension of the left-hand side if
$\mathbf{a}_{[k]}-a_{k}$ contains integers.

######  Remark 3.8.

For $p=2$ this is equivalent to Ptolemy’s theorem: if a quadrilateral is
inscribed in a circle then the product of the lengths of its diagonals is
equal to the sum of the products of the lengths of the pairs of opposite
sides, which can be written as

$\displaystyle\sin(\theta_{3}-\theta_{1})\sin(\theta_{4}-\theta_{2})=\sin(\theta_{2}-\theta_{1})\sin(\theta_{4}-\theta_{3})+\sin(\theta_{4}-\theta_{1})\sin(\theta_{3}-\theta_{2}).$

Further details regarding the history behind the identity (3.14) can be found
in the introduction and [21].

######  Corollary 3.9.

For each $m\in\mathbb{N}_{0}$ and each $p\in\mathbb{N}$ the function

$\displaystyle
F_{p,m}(\mathbf{a},\mathbf{b})=\sum\limits_{j=0}^{m}\frac{(-1)^{j}}{j!}[a_{k}]_{j}[\psi_{p}+m-1]_{j}g^{k}_{p}(m-j)$
$\displaystyle\hphantom{F_{p,m}(\mathbf{a},\mathbf{b})}{}=\sum\limits_{j=0}^{m}\sum\limits_{r=0}^{m-j}\frac{(-1)^{m-j}(\psi_{p}+r)_{m-r}}{j!(m-j-r)!}l_{r}[a_{k}]_{j}\mathcal{B}^{(m-j+\psi_{p})}_{m-j-r}(1-a_{k})$
(3.15)

is independent of $k$ and represents a symmetric polynomial in the components
of $\mathbf{a}$ and $\mathbf{b}$ $($separately$)$. Here $l_{r}$ is defined by
the recurrence relation (3.6).

###### Proof.

Indeed, the first formula in (3.15) follows from (3.8) once we open the braces
and apply the obvious relation
$1/(\psi_{p})_{m-j}=[\psi_{p}+m-1]_{j}/(\psi_{p})_{m}$. Substitution of (3.4)
for $g^{k}_{p}$ leads to the second formula. ∎

The following theorem can be viewed as a new method for computing the
coefficients $h_{p}(n|\boldsymbol{\alpha},\boldsymbol{\beta})$ in expansion
(2.19) given by the multiple sum (2.24) by relating them to the numbers
$D_{n}^{[k,s]}$ given by the single sums (2.17) and (2.18).

######  Theorem 3.10.

For each nonnegative integer $n$ and arbitrary $s\in\\{1,\ldots,p\\}$ the
following identity holds

$\displaystyle
h_{p}^{s}(n)=-\frac{1}{\pi\sin(\pi\psi_{p})}\sum_{\begin{subarray}{c}k=1\\\
k\neq{s}\end{subarray}}^{p}\frac{\sin(\pi(\mathbf{b}-a_{k}))}{\sin(\pi(\mathbf{a}_{[k,s]}-a_{k}))}D_{n}^{[k,s]},$
(3.16)

where $h_{p}^{s}(n)$ is defined by expansion (2.20) and given explicitly by
(2.24), while $D_{n}^{[k,s]}$ are given by (2.17) or (2.18). Moreover, for
arbitrary distinct integers $s$, $i$, $k$ from the set $\\{1,2,\ldots,p\\}$
the following identity holds

$\displaystyle\sum\limits_{j=0}^{n}\frac{(-1)^{j}}{j!}\left([a_{i}]_{j}\sin(\pi(a_{s}-a_{i}))D_{n-j}^{[s,i]}+[a_{k}]_{j}\sin(\pi(a_{i}-a_{k}))D_{n-j}^{[i,k]}\right.$
$\displaystyle\left.\hphantom{\sum\limits_{j=0}^{n}\frac{(-1)^{j}}{j!}}{}+[a_{s}]_{j}\sin(\pi(a_{k}-a_{s}))D_{n-j}^{[k,s]}\right)=0.$
(3.17)

###### Proof.

To prove (3.16) it suffices to substitute expansion (2.16) into formula (2.22)
and equate coefficients. Identity (3.17) is a direct consequence of (2.11). ∎

######  Corollary 3.11.

For each $n\in\mathbb{N}_{0}$ the following identity holds true

$\displaystyle\frac{\Gamma(a_{3}-a_{1})}{\Gamma(\mathbf{b}-a_{1})\Gamma(2+a_{1}+a_{2}-b_{1}-b_{2}+n)}\,{}_{3}F_{2}\left(\begin{matrix}1+a_{1}-b_{1},1+a_{1}-b_{2},b_{3}-a_{3}\\\
1+a_{1}-a_{3},2+a_{1}+a_{2}-b_{1}-b_{2}+n\end{matrix}\right)$
$\displaystyle\quad{}+\frac{\Gamma(a_{1}-a_{3})}{\Gamma(\mathbf{b}-a_{3})\Gamma(2+a_{2}+a_{3}-b_{1}-b_{2}+n)}\,{}_{3}F_{2}\left(\begin{matrix}1+a_{3}-b_{1},1+a_{3}-b_{2},b_{3}-a_{1}\\\
1+a_{3}-a_{1},2+a_{2}+a_{3}-b_{1}-b_{2}+n\end{matrix}\right)$
$\displaystyle{}=\frac{1}{\Gamma(2-\psi+n)\Gamma(\psi+a_{2}-b_{1})\Gamma(\psi+a_{2}-b_{2})}\,{}_{3}F_{2}\left(\begin{matrix}b_{3}-a_{1},b_{3}-a_{3},\psi-1-n\\\
\psi+a_{2}-b_{1},\psi+a_{2}-b_{2}\end{matrix}\right).$

###### Proof.

For $p=3$, $s=2$ formula (3.16) takes the form

$\displaystyle
h_{3}^{2}(n)=-\frac{1}{\pi\sin(\pi\psi)}\left(\frac{\sin(\pi(\mathbf{b}-a_{1}))}{\sin(\pi(a_{3}-a_{1})}D_{n}^{[1,2]}+\frac{\sin(\pi(\mathbf{b}-a_{3}))}{\sin(\pi(a_{1}-a_{3})}D_{n}^{[3,2]}\right).$
(3.18)

Now, write $h_{3}^{2}(n)$ according to (2.25) and exchange the roles of
$b_{1}$ and $b_{3}$. Next, apply Chu–Vandermonde identity on the right-hand
side of (2.17) to get

$\displaystyle
D_{n}^{[k,s]}=\frac{\Gamma(1-\mathbf{b}+a_{k})\Gamma(1-b_{1}+a_{s}+n)\Gamma(1-b_{2}+a_{s}+n)}{\Gamma(1-\mathbf{a}_{[k,s]}+a_{k})\Gamma(2+a_{k}+a_{s}-b_{1}-b_{2}+n)n!}$
$\displaystyle\hphantom{D_{n}^{[k,s]}=}{}\times{}_{3}F_{2}\left(\begin{matrix}1-b_{1}+a_{k},1-b_{2}+a_{k},b_{3}-\mathbf{a}_{[k,s]}\\\
2+a_{k}+a_{s}-b_{1}-b_{2}+n,1-\mathbf{a}_{[k,s]}+a_{k}\end{matrix}\right).$

Substituting this into (3.18), applying Euler’s reflection formula for the
gamma function and rearranging we get the claimed identity. ∎

##  Appendix A Definition of Meijer’s $\boldsymbol{G}$-function revisited

Meijer’s $G$-function has been defined in the introduction, where we mentioned
various aspects that need to be clarified in order that this definition be
consistent. Most accurate information with proofs regarding $G$-function’s
definition is contained, in our opinion, in the series of papers of Meijer
himself [30], the paper [6] by Braaksma and in the first chapters of the books
[35] and [25]. Further facts are scattered in the literature with most
comprehensive collection being [36, Chapter 8], [3] and especially [31]. An
accessible introduction to $G$-function can be found in a nice recent survey
by Beals and Szmigielski [4]. In this paper we only deal with the function
$G^{m,n}_{p,p}$. For convenience, we have gathered all the necessary
information regarding its definition in the following theorem.

######  Theorem A.1.

Denote $a^{*}=m+n-p$, $\psi=\sum\limits_{k=1}^{p}(b_{k}-a_{k})$ and

$\displaystyle\mathcal{G}(s)=\frac{\Gamma(b_{1}+s)\cdots\Gamma(b_{m}+s)\Gamma(1-a_{1}-s)\cdots\Gamma(1-a_{n}-s)}{\Gamma(a_{n+1}+s)\cdots\Gamma(a_{p}+s)\Gamma(1-b_{m+1}-s)\cdots\Gamma(1-b_{p}-s)}.$

  1. $(a)$

If $|z|<1$ then the integral in (2.1) converges for
$\mathcal{L}=\mathcal{L}_{-}$ and

$\displaystyle G^{m,n}_{p,p}\left(z\,\vline\,\begin{matrix}\mathbf{a}\\\
\mathbf{b}\end{matrix}\right)=\sum\limits_{j=1}^{m}\sum\limits_{l=0}^{\infty}\operatorname*{\mathrm{res}}_{s=b_{jl}}\mathcal{G}(s)z^{-s},\qquad
b_{jl}=-b_{j}-l.$ (A.1)

If, in addition, $a^{*}>0$ and $|\arg(z)|<a^{*}\pi$ or $a^{*}=0$,
$\Re(\psi)<0$ and $0<z<1$ then the integral in (2.1) also converges for
$\mathcal{L}=\mathcal{L}_{i\gamma}$ and has the same value.

  2. $(b)$

If $|z|>1$ then the integral in (2.1) converges for
$\mathcal{L}=\mathcal{L}_{+}$ and

$\displaystyle G^{m,n}_{p,p}\left(z\,\vline\,\begin{matrix}\mathbf{a}\\\
\mathbf{b}\end{matrix}\right)=-\sum\limits_{i=1}^{n}\sum\limits_{k=0}^{\infty}\operatorname*{\mathrm{res}}_{s=a_{ik}}\mathcal{G}(s)z^{-s},\qquad
a_{ik}=1-a_{i}+k.$

If, in addition, $a^{*}>0$ and $|\arg(z)|<a^{*}\pi$ or $a^{*}=0$,
$\Re(\psi)<0$ and $z>1$ then the integral in (2.1) also converges for
$\mathcal{L}=\mathcal{L}_{i\gamma}$ and has the same value.

######  Remark A.2.

If $|z|=1$ and $\Re(\psi)<-1$ then according to [25, Theorem 1.1] both
integrals over $\mathcal{L}_{-}$ and over $\mathcal{L}_{+}$ exist. If, in
addition, $a^{*}>0$ and $|\arg(z)|<a^{*}\pi$, then the integral over
$\mathcal{L}_{i\gamma}$ also exists. We found no proof in the literature that
these integrals are equal.

###### Proof.

The claims regarding the contours $\mathcal{L}_{+}$ and $\mathcal{L}_{-}$ have
been demonstrated in [25, Theorems 1.1, 1.2]. They follow in a relatively
straightforward manner from Stirling’s asymptotic formula for the gamma
function and have been observed by Meijer himself in [30, Section 1]. If
$a^{*}=0$, $\Re\psi<0$ the result was proved by Kilbas and Saigo [25, Theorem
3.3]. For $a^{*}>0$ and $|z|<1$ [30, Theorem D] states that the integrals over
$\mathcal{L}_{-}$ and $\mathcal{L}_{i\gamma}$ coincide for any real $\gamma$
which covers case (a) of the theorem. For the proof Meijer refers to [43,
Section 14.51], where it is essentially demonstrated that

$\displaystyle\lim_{R\to\infty}\int_{l^{-}_{\gamma}(R)}\frac{\prod\limits_{j=1}^{m}\Gamma(s+b_{j})\prod\limits_{i=1}^{n}\Gamma(1-s-a_{i})}{\prod\limits_{i=n+1}^{p}\Gamma(s+a_{i})\prod\limits_{j=m+1}^{p}\Gamma(1-s-b_{j})}z^{-s}ds=0,$
(A.2)

where $l^{-}_{\gamma}(R)$ is the shortest arc of the circle $|s|=R$ which
connects the contours $\mathcal{L}_{-}$ and $\mathcal{L}_{i\gamma}$ (in fact,
Whittaker and Watson proved a particular case, but the proof for the general
case goes along exactly the same lines). Since $\mathcal{G}(s)$ has no poles
in the domain bounded by $\mathcal{L}_{-}$ and $\mathcal{L}_{i\gamma}$ this
leads to equality of the integrals over $\mathcal{L}_{-}$ and
$\mathcal{L}_{i\gamma}$ stated by Meijer. It remains to consider the case
$|z|>1$, $\mathcal{L}=\mathcal{L}_{i\gamma}$. Denote by $l^{+}_{\gamma}(R)$
the reflection of $l^{-}_{\gamma}(R)$ with respect to $\mathcal{L}_{i\gamma}$.
For $|z|>1$ we get by changing $s$ to $-s$

$\displaystyle\int_{l^{+}_{\gamma}(R)}\frac{\prod\limits_{j=1}^{m}\Gamma(s+b_{j})\prod\limits_{i=1}^{n}\Gamma(1-s-a_{i})z^{-s}}{\prod\limits_{i=n+1}^{p}\Gamma(s+a_{i})\prod\limits_{j=m+1}^{p}\Gamma(1-s-b_{j})}ds$
$\displaystyle\qquad{}=-\int_{l^{-}_{\gamma}(R)}\frac{\prod\limits_{j=1}^{m}\Gamma(-s+b_{j})\prod\limits_{i=1}^{n}\Gamma(1+s-a_{i})z^{s}}{\prod\limits_{i=n+1}^{p}\Gamma(-s+a_{i})\prod\limits_{j=m+1}^{p}\Gamma(1+s-b_{j})}ds$
$\displaystyle\qquad{}=-\int_{l^{-}_{\gamma}(R)}\frac{\prod\limits_{i=1}^{n}\Gamma(s+b_{i}^{\prime})\prod\limits_{j=1}^{m}\Gamma(1-s-a_{j}^{\prime})}{\prod\limits_{j=m+1}^{p}\Gamma(s+a_{j}^{\prime})\prod\limits_{i=n+1}^{p}\Gamma(1-s-b_{i}^{\prime})}(1/z)^{-s}ds.$

Here $b_{i}^{\prime}=1-a_{i}$, $a_{j}^{\prime}=1-b_{j}$. In view of (A.2), we
immediately conclude that for $|z|>1$

$\displaystyle\lim_{R\to\infty}\int_{l^{+}_{\gamma}(R)}\frac{\prod\limits_{j=1}^{m}\Gamma(s+b_{j})\prod\limits_{i=1}^{n}\Gamma(1-s-a_{i})}{\prod\limits_{i=n+1}^{p}\Gamma(s+a_{i})\prod\limits_{j=m+1}^{p}\Gamma(1-s-b_{j})}z^{-s}ds=0,$

which implies that the integrals over $\mathcal{L}_{+}$ and
$\mathcal{L}_{i\gamma}$ coincide. ∎

######  Remark A.3.

If $p>q$ ($q>p$) the integral in (2.1) exists for
$\mathcal{L}=\mathcal{L}_{+}$ ($\mathcal{L}=\mathcal{L}_{-}$) and all complex
$z\neq{0}$ and is equal to the corresponding sum of residues [25, Theorems 1.1
and 1.2]. At the same time if $a^{*}>0$ and $|\arg(z)|<a^{*}\pi$ or $a^{*}=0$
and $z>0$, $z\neq{1}$, the integral in (2.1) also exists for
$\mathcal{L}=\mathcal{L}_{i\gamma}$. Most authors assume in this situation
that $\mathcal{L}_{i\gamma}$ can be deformed into
$\mathcal{L}=\mathcal{L}_{+}$ if $p>q$ or $\mathcal{L}=\mathcal{L}_{-}$ if
$q>p$ without altering the value of the integral. However, we were unable to
find any proof of this claim in the literature.

######  Remark A.4.

It follows from the above theorem that $G^{m,n}_{p,p}(z)$ is analytic in the
sector $|\arg(z)|$ $<a^{*}\pi$ if $a^{*}>0$ (since the integral converges
uniformly in $z$ for $\mathcal{L}=\mathcal{L}_{i\gamma}$), while for
$a^{*}\leq 0$ we get two different analytic functions – one defined inside and
the other outside of the unit circle, see [36, (8.2.2.7)].

In the proof of the above theorem we essentially used the next well-known
reflection property of $G$-function

$\displaystyle
G^{m,n}_{p,q}\left(\frac{1}{z}\,\vline\,\begin{matrix}\mathbf{a}\\\
\mathbf{b}\end{matrix}\right)=G^{n,m}_{q,p}\left(z\,\vline\,\begin{matrix}1-\mathbf{b}\\\
1-\mathbf{a}\end{matrix}\right).$

It is important to note that by Theorem A.1(b)

$\displaystyle G^{p,0}_{p,p}\left(z\,\vline\,\begin{matrix}\mathbf{b}\\\
\mathbf{a}\end{matrix}\right)=0\qquad\text{for}\quad|z|>1,$

which is, of course, different from the analytic continuation of the right-
hand side of (2.4). The Mellin transform of $G^{p,0}_{p,p}$ exists if either
$\Re(\psi)>0$ or $\psi=-m\in\mathbb{N}_{0}$. In the former case [25, Theorem
2.2]

$\displaystyle\int_{0}^{\infty}x^{s-1}G^{p,0}_{p,p}\left(x\,\vline\,\begin{matrix}\mathbf{b}\\\
\mathbf{a}\end{matrix}\right)dx=\int_{0}^{1}x^{s-1}G^{p,0}_{p,p}\left(x\,\vline\,\begin{matrix}\mathbf{b}\\\
\mathbf{a}\end{matrix}\right)dx=\frac{\Gamma(\mathbf{a}+s)}{\Gamma(\mathbf{b}+s)}$
(A.3)

for $\Re(s)>-\Re(\mathbf{a})$. If $\psi=-m\in\mathbb{N}_{0}$ then [33, (2.28)]

$\displaystyle\int_{0}^{\infty}x^{s-1}G^{p,0}_{p,p}\left(x\,\vline\,\begin{matrix}\mathbf{b}\\\
\mathbf{a}\end{matrix}\right)dx=\int_{0}^{1}x^{s-1}G^{p,0}_{p,p}\left(x\,\vline\,\begin{matrix}\mathbf{b}\\\
\mathbf{a}\end{matrix}\right)dx=\frac{\Gamma(\mathbf{a}+s)}{\Gamma(\mathbf{b}+s)}-q(s)$

for $\Re(s)>-\Re(\mathbf{a})$. Here $q(s)$ is a polynomial of degree $m$ given
by

$\displaystyle q(s)=\sum\limits_{j=0}^{m}g^{k}_{p}(m-j)(s+a_{k}-j)_{j},\qquad
k=1,2,\ldots,p,$

where the coefficients $g^{k}_{p}(n)$ are defined in expansion (2.5) and are
given explicitly by (2.9). Note that $g^{k}_{p}(n)$ depends on $k$ while the
polynomial $q(s)$ is the same for each $k\in\\{1,\ldots,p\\}$.

### Acknowledgements

This work has been supported by the Russian Science Foundation under project
14-11-00022. We are also indebted to anonymous referees for a number of useful
remarks that helped to improve the exposition.

## References

  * [1] Akemann G., Ipsen J.R., Kieburg M., Products of rectangular random matrices: singular values and progressive scattering, Phys. Rev. E 88 (2013), 052118, 13 pages, arXiv:1307.7560.
  * [2] Andrews G.E., Askey R., Roy R., Special functions, Encyclopedia of Mathematics and its Applications, Vol. 71, Cambridge University Press, Cambridge, 1999.
  * [3] Askey R.A., Olde Daalhuis A.B., Generalized hypergeometric functions and Meijer $G$-function, in NIST Handbook of Mathematical Functions, Editors F.W.J. Olver, D.W. Lozier, R.F. Boisvert, C.W. Clark, U.S. Department of Commerce, National Institute of Standards and Technology, Washington, DC, Cambridge University Press, Cambridge, 2010, Chapter 16, 403–418, available at http://dlmf.nist.gov/.
  * [4] Beals R., Szmigielski J., Meijer $G$-functions: a gentle introduction, Notices Amer. Math. Soc. 60 (2013), 866–872.
  * [5] Beukers F., Heckman G., Monodromy for the hypergeometric function ${}_{n}F_{n-1}$, Invent. Math. 95 (1989), 325–354.
  * [6] Braaksma B.L.J., Asymptotic expansions and analytic continuations for a class of Barnes-integrals, Compositio Math. 15 (1964), 239–341.
  * [7] Bühring W., The behavior at unit argument of the hypergeometric function ${}_{3}F_{2}$, SIAM J. Math. Anal. 18 (1987), 1227–1234.
  * [8] Bühring W., Generalized hypergeometric functions at unit argument, Proc. Amer. Math. Soc. 114 (1992), 145–153.
  * [9] Bühring W., Transformation formulas for terminating Saalschützian hypergeometric series of unit argument, J. Appl. Math. Stochastic Anal. 8 (1995), 189–194.
  * [10] Bühring W., Srivastava H.M., Analytic continuation of the generalized hypergeometric series near unit argument with emphasis on the zero-balanced series, in Approximation Theory and Applications, Hadronic Press, Palm Harbor, FL, 1998, 17–35, math.CA/0102032.
  * [11] Chamayou J.F., Letac G., Additive properties of the Dufresne laws and their multivariate extension, J. Theoret. Probab. 12 (1999), 1045–1066.
  * [12] Coelho C.A., Arnold B.C., Instances of the product of independent beta random variables and of the Meijer $G$ and Fox $H$ functions with finite representations, AIP Conf. Proc. 1479 (2012), 1133–1137.
  * [13] Consul P.C., The exact distributions of likelihood criteria for different hypotheses, in Multivariate Analysis, II (Proc. Second Internat. Sympos., Dayton, Ohio, 1968), Editor P.R. Krishnaian, Academic Press, New York, 1969, 171–181.
  * [14] Davis A.W., On the differential equation for Meijer’s $G^{p,0}_{p,p}$ function, and further tables of Wilks’s likelihood ratio criterion, Biometrika 66 (1979), 519–531.
  * [15] Dufresne D., $G$ distributions and the beta-gamma algebra, Electron. J. Probab. 15 (2010), no. 71, 2163–2199.
  * [16] Dunkl C.F., Products of beta distributed random variables, arXiv:1304.6671.
  * [17] Erdélyi A., Magnus W., Oberhettinger F., Tricomi F.G., Higher transcendental functions, Vol. I, McGraw-Hill Book Company, Inc., New York – Toronto – London, 1953.
  * [18] Feng R., Kuznetsov A., Yang F., New identies for finite sums of products of generalized hypergeometric functions, arXiv:1512.01121.
  * [19] Gosper R.W., Ismail M.E.H., Zhang R., On some strange summation formulas, Illinois J. Math. 37 (1993), 240–277.
  * [20] Hermite C., Sur une identité trigonométrique, Nouv. Ann. Mat. 4 (1885), 57–59.
  * [21] Johnson W.P., Trigonometric identities à la Hermite, Amer. Math. Monthly 117 (2010), 311–327.
  * [22] Kalinin V.M., Special functions and limit properties of probability distributions. I, Zap. Naučn. Sem. Leningrad. Otdel. Mat. Inst. Steklov. (LOMI) 13 (1969), 5–137.
  * [23] Karp D.B., Prilepkina E.G., Completely monotonic gamma ratio and infinitely divisible $H$-function of Fox, Comput. Methods Funct. Theory 16 (2016), 135–153, arXiv:1501.05388.
  * [24] Karp D.B., Prilepkina E.G., Some new facts around the delta neutral $H$ function of Fox, arXiv:1511.06612.
  * [25] Kilbas A.A., Saigo M., $H$-transforms. Theory and applications, Analytical Methods and Special Functions, Vol. 9, Chapman & Hall/CRC, Boca Raton, FL, 2004.
  * [26] Marichev O.I., On the representation of Meijer’s $G$-function in the vicinity of singular unity, in Complex Analysis and Applications ’81 (Varna, 1981), Publ. House Bulgar. Acad. Sci., Sofia, 1984, 383–398.
  * [27] Marichev O.I., Kalla S.L., Behaviour of hypergeometric function ${}_{p}F_{p-1}(z)$ in the vicinity of unity, Rev. Técn. Fac. Ingr. Univ. Zulia 7 (1984), no. 2, 1–8.
  * [28] Mathai A.M., A review of the different techniques used for deriving the exact distributions of multivariate test criteria, Sankhyā Ser. A 35 (1973), 39–60.
  * [29] Mathai A.M., Extensions of Wilks’ integral equations and distributions of test statistics, Ann. Inst. Statist. Math. 36 (1984), 271–288.
  * [30] Meijer C.S., On the $G$-function, I–VIII, Nederl. Akad. Wetensch. Proc. Ser. A. 49 (1946), 227–237, 344–356, 457–469, 632–641, 765–772, 936–943, 1063–1072, 1165–1175.
  * [31] Meijer $G$-functions, Wolfram Functions Site, http://functions.wolfram.com/HypergeometricFunctions/MeijerG/.
  * [32] Nair U.S., The application of the moment function in the study of distribution laws in statistics, Biometrika 30 (1939), 274–294.
  * [33] Nørlund N.E., Hypergeometric functions, Acta Math. 94 (1955), 289–349.
  * [34] Nörlund N.E., Sur les valeurs asymptotiques des nombres et des polynômes de Bernoulli, Rend. Circ. Mat. Palermo 10 (1961), 27–44.
  * [35] Paris R.B., Kaminski D., Asymptotics and Mellin–Barnes integrals, Encyclopedia of Mathematics and its Applications, Vol. 85, Cambridge University Press, Cambridge, 2001.
  * [36] Prudnikov A.P., Brychkov Y.A., Marichev O.I., Integrals and series. Vol. 3. More special functions, Gordon and Breach Science Publishers, New York, 1990.
  * [37] Saigo M., Srivastava H.M., The behavior of the zero-balanced hypergeometric series ${}_{p}F_{p-1}$ near the boundary of its convergence region, Proc. Amer. Math. Soc. 110 (1990), 71–76.
  * [38] Scheidegger F., Analytic continuation of hypergeometric functions in the resonant case, arXiv:1602.01384.
  * [39] Springer M.D., Thompson W.E., The distribution of products of beta, gamma and Gaussian random variables, SIAM J. Appl. Math. 18 (1970), 721–737.
  * [40] Tang J., Gupta A.K., On the distribution of the product of independent beta random variables, Statist. Probab. Lett. 2 (1984), 165–168.
  * [41] Tang J., Gupta A.K., Exact distribution of certain general test statistics in multivariate analysis, Austral. J. Statist. 28 (1986), 107–114.
  * [42] Thomae J., Ueber die höheren hypergeometrischen Reihen, insbesondere über die Reihe: $1+\frac{{a_{0}a_{1}a_{2}}}{{1.b_{1}b_{2}}}x+\frac{{a_{0}(a_{0}+1)a_{1}(a_{1}+1)a_{2}(a_{2}+1)}}{{1.2.b_{1}(b_{1}+1)b_{2}(b_{2}+1)}}x^{2}+\cdots$, Math. Ann. 2 (1870), 427–444.
  * [43] Whittaker E.T., Watson G.N., A course of modern analysis. An introduction to the general theory of infinite processes and of analytic functions; with an account of the principal transcendental functions, Cambridge Mathematical Library, Cambridge University Press, Cambridge, 1996.
  * [44] Wilks S.S., Certain generalizations in the analysis of variance, Biometrika 24 (1932), 471–494.