Hypergeometric Differential Equation and New Identities for the Coefficients of N{\o}rlund and B\"uhring

The fundamental set of solutions of the generalized hypergeometric differential equation in the neighborhood of unity has been built by N{\o}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\"uhring, 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.


Introduction
We will use standard notation Z, N and C to denote integer, natural and complex numbers, respectively; N 0 = N ∪ {0}. The hypergeometric differential equation (D = z d dz ) 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, ∞. The local exponents read [5, (2.6)-(2.8)] 1 z , k = 1, . . . , 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] 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, 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

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 ≤ m ≤ q, 0 ≤ n ≤ p are integers and a, b are arbitrary complex vectors, such that a i − b j / ∈ N for all i = 1, . . . , n and j = 1, . . . , 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]), where the contour L is a simple loop that separates the poles of the integrand of the form b jl = −b j − l, l ∈ N 0 , leaving them on the left from the poles of the form a ik = 1 − a i + k, k ∈ N 0 , leaving them on the right [25, Section 1.1]. It may have one of the three forms L − , L + or L iγ described below. Choose any and arbitrary real γ. The contour L − is a left loop lying in the horizontal strip ϕ 1 ≤ s ≤ ϕ 2 . It starts at the point −∞ + iϕ 1 , terminates at the point −∞ + iϕ 2 and coincides with the boundary of the strip for sufficiently large |s|. Similarly, the contour L + is a right loop lying in the same strip, starting at the point +∞ + iϕ 1 and terminating at the point +∞ + iϕ 2 . Finally, the contour L iγ starts at γ − i∞, terminates at γ + i∞ and coincides with the line s = γ for all sufficiently large |s|. The power function z −s is defined on the Riemann surface of the logarithm, so that 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 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 where α ∈ C and a + α is understood as (a 1 + α, . . . , a p + α).
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.
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: sin(a) = sin a 1 sin a 2 · · · sin a p , Γ(a) = Γ(a 1 )Γ(a 2 ) · · · Γ(a p ) and let a [k,s] denote the vector a with the elements a k and a s removed. We write (a) > 0 for (a i ) > 0, i = 1, . . . , p. The next theorem is implicit in [33].
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] [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 Gfunction by the residue theorem. Namely, if the elements of the vector a are different modulo integers, formula (A.1) applied to the function G p,0 p,p takes the form (see also [26, According to Theorem A.1(a) this formula holds for |z| < 1.
2.2 Expansion of G p,0 p,p in the neighborhood of unity Among many other results contained in [33] Nørlund showed that the series represents a solution in the neighborhood of z = 1 corresponding to the local exponent ψ p − 1 (see (2.1) for the definition of ψ p ) if this number is not a negative integer. Formula (2.5) holds in the disk |1 − z| < 1 for all −ψ p / ∈ N 0 and each k = 1, 2, . . . , p. For −ψ p = l ∈ N 0 we have by taking limit in (2.5) (see [33, (1.34)]): The value of k in (2.5) and (2.6) can be chosen arbitrarily from the set {1, 2, . . . , p}. This choice affects the first factor z a k and the coefficients g k p (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 where P m (z | a, b) is a polynomial in z of degree m with coefficients dependent on the parameters a and b. It is expressed in terms of the polynomials as follows [33, (1.28)]: for j = 2, . . . , p − 1, and Here the difference is understood as the forward difference The initial values g k p (0), g k p (1), . . ., g k p (p − 1) for the recurrence (2.7) are found by solving the next triangular system: )g k p (0) = 0, · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · Another method to compute the coefficients g k p (n) discovered by Nørlund is the following recurrence in p [33, (2.7)]: The last formula can be applied for any k = 1, 2, . . . , p − 1 without affecting the left-hand side. The value of g k p (n) for k = 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)]: Furthermore, he solved the above recurrence to obtain [33, (2.11)]: where, as before, 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 b and a [k] (separately); • the summation in (2.9) is over all Young diagrams that fit n × (p − 2) box.
Symmetry follows from expansion (2.5) and the invariance of G-function with respect to permutation the elements of a and b. In Section 3 below we find another way to compute the coefficients g k p (n) in terms of generalized Bernoulli polynomials and we further give explicit formulas for g k p (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 ψ 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,0 p,p (e −x ) is infinitely divisible distribution on [0, ∞).

Expansion of G 2,p
p,p in the neighborhood of unity To define y k,s (x) the roles of γ 1 , γ 2 are exchanged with those of γ k , γ 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 The series in ( 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 Here and below p F p−1 without an argument is understood as p F p−1 (1). The series (2.16) con-

Connection to Bühring expansion of p F p−1
In his 1992 paper Bühring found the representation [8, Theorem 2] α k . He also derived explicit formulas for the coefficients f p (n|α, β) and h p (n|α, β) to be given below. On setting α = 1 − b + a s , β = 1 − a [s] + a s , ν = ψ p − 1 and denoting Bühring's formula takes the form  where Euler's reflection formula has been used. Thus, Bühring coefficients f p (n|α, β) (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|α, β) found in [8, (2.9), (2.16)] after some renaming of variables and setting s = p can be put into the form (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 where j 0 = 0 and j p−1 = n. For p = 3 this formula reduces to the identity which is a guise of Sheppard's transformation [ 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 s p (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 s p (n) have been first found in [8]. Adopted to our notation Bühring's formulas [ where, as before, ψ m = m i=1 (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  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 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 . Notwithstanding the non-symmetric appearance, the last formula is symmetric with respect to the elements of b. For p = 4 we convert the expression for f s p (n) calculated in [10, (3.17)] into expressions for g s p (n) using (2.23) and the necessary renaming of parameters. This yields . This formula is also invariant with respect to permutation of the elements of b.

4)
Herel 0 = 1 andl r , r ≥ 1, are found from the recurrencẽ Similarly, the coefficients l r satisfy the recurrence relation and 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 : Similar formula is of course true forl r once we write 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 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, A = B = (1, . . . , 1), ν = 1, µ = ψ p , θ = a k − 1 in that theorem.
where we applied Euler's reflection formula Γ(z)Γ(1 − z) = π/ sin(πz). Further, substitute so that the above equality reduces to Assume for a moment that ψ p is not an integer. We know that both series in (3.9) converge in |z − 1| < 1. Furthermore, all functions γ n,s (z) and χ 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 which we prove by induction in m.
Letting z → 1 in the second identity in (3.10) we get χ 0 (1) = 0 which establishes our claim for m = 0. Next, suppose (3.11) holds for m = 0, 1, 2, . . . , r − 1. Divide the second identity in (3.10) by (1 − z) r and expand each χ n (z), n = 0, 1, . . . , r, in Taylor series around z = 1: where the last equality is by induction hypothesis. We now obtain (3.11) for m = r on letting z → 1 in this formula. It is immediate to check that the claimed identity Finally, we remove the assumption that ψ 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 ψ p then clearly imply by analytic continuation that these poles are removable and both identities hold for all ψ p . Corollary 3.6. Suppose for some i ∈ {1, 2, 3} inequality (b i ) < (a k +1) holds for k ∈ {1, 2, 3}. Then 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 k 3 (0) and h k 3 (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. sin(π(b − a k )) sin(π(a [k] − a k )) = sin(πψ p ). (3.14) The right-hand side gives a continuous extension of the left-hand side if 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 Further details regarding the history behind the identity (3.14) can be found in the introduction and [21].
is independent of k and represents a symmetric polynomial in the components of a and 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/(ψ p ) m−j = [ψ p + m − 1] j /(ψ 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|α, β) in expansion (2.19) given by the multiple sum (2.24) by relating them to the numbers D 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 ∈ N 0 the following identity holds true Proof . For p = 3, s = 2 formula (3.16) takes the form Now, write h 2 3 (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 Substituting this into (3.18), applying Euler's reflection formula for the gamma function and rearranging we get the claimed identity.

A Def inition of Meijer's 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. .
Remark A.2. If |z| = 1 and (ψ) < −1 then according to [25,Theorem 1.1] both integrals over L − and over L + exist. If, in addition, a * > 0 and | arg(z)| < a * π, then the integral over L iγ also exists. We found no proof in the literature that these integrals are equal. where l − γ (R) is the shortest arc of the circle |s| = R which connects the contours L − and L iγ (in fact, Whittaker and Watson proved a particular case, but the proof for the general case goes along exactly the same lines). Since G(s) has no poles in the domain bounded by L − and L iγ this leads to equality of the integrals over L − and L iγ stated by Meijer. It remains to consider the case |z| > 1, L = L iγ . Denote by l + γ (R) the reflection of l − γ (R) with respect to L iγ . For |z| > 1 we get by changing s to −s Here b i = 1 − a i , a j = 1 − b j . In view of (A.2), we immediately conclude that for |z| > 1 which implies that the integrals over L + and L iγ coincide.
Remark A.3. If p > q (q > p) the integral in (2.1) exists for L = L + (L = L − ) and all complex z = 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 * π or a * = 0 and z > 0, z = 1, the integral in (2.1) also exists for L = L iγ . Most authors assume in this situation that L iγ can be deformed into L = L + if p > q or L = 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 * π if a * > 0 (since the integral converges uniformly in z for L = L iγ ), while for a * ≤ 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 G m,n p,q 1 z It is important to note that by Theorem A.1(b) G p,0 p,p z b a = 0 for |z| > 1, which is, of course, different from the analytic continuation of the right-hand side of (2.4  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 ∈ {1, . . . , p}.