Multiple orthogonal analogs of classical orthogonal polynomials

Symmetry, Integrability and Geometry: Methods and Applications SIGMA 14 (2018), 063, 27 pages d-Orthogonal Analogs of Classical Orthogonal Polynomials arXiv:1609.06157v3 [math.CA] 26 Jun 2018 Emil HOROZOV †‡ † Department of Mathematics and Informatics, Sofia University, 5 J. Bourchier Blvd., Sofia 1126, Bulgaria E-mail: horozov@fmi.uni-sofia.bg ‡ Institute of Mathematics and Informatics, Bulg. Acad. of Sci., Acad. G. Bonchev Str., Block 8, 1113 Sofia, Bulgaria Received October 01, 2017, in final form June 13, 2018; Published online June 26, 2018 https://doi.org/10.3842/SIGMA.2018.063 Abstract. Classical orthogonal polynomial systems of Jacobi, Hermite and Laguerre have the property that the polynomials of each system are eigenfunctions of a second order ordinary differential operator. According to a famous theorem by Bochner they are the only systems on the real line with this property. Similar results hold for the discrete orthogonal polynomials. In a recent paper we introduced a natural class of polynomial systems whose members are the eigenfunctions of a differential operator of higher order and which are orthogonal with respect to d measures, rather than one. These polynomial systems, enjoy a number of properties which make them a natural analog of the classical orthogonal polynomials. In the present paper we continue their study. The most important new properties are their hypergeometric representations which allow us to derive their generating functions and in some cases also Mehler–Heine type formulas. Key words: d-orthogonal polynomials; finite recurrence relations; bispectral problem; generalized hypergeometric functions; generating functions 2010 Mathematics Subject Classification: 34L20; 30C15; 33E05 1 Introduction This paper is a natural continuation of the study initiated in [38] on the basis of several classes of examples in [36, 37]. There we constructed large families of polynomial systems that were called d-orthogonal polynomials with the Bochner property. The terminology “Bochner’s property” derives from the Bochner’s theorem [19] mentioned in the abstract and means that the polynomials are eigenfunction of a differential operator. We recall that, by definition, general d-orthogonal polynomials are polynomial systems Pn (x), n = 0, 1, 2, . . ., deg(Pn ) = n iff there exist d linear functionals Lj , j = 0, . . . , d − 1 on the space of all polynomials C[x] such that ( Lj (Pn Pm ) = 0, Lj (Pn Pnd+j ) 6= 0, m > nd + j, n ≥ 0, n ≥ 0, for each j ∈ Nd := {0, . . . , d − 1}. When d = 1 this is the ordinary notion of orthogonal polynomials. The orthogonality is connected with d functionals rather than with only one. According to [47, 61] the above property is equivalent to the existence of a linear recurrence This paper is a contribution to the Special Issue on Orthogonal Polynomials, Special Functions and Applications (OPSFA14). The full collection is available at https://www.emis.de/journals/SIGMA/OPSFA2017.html 2 E.I. Horozov relation of the form xPn (x) = Pn+1 + d X γj (n)Pn−j (x). j=0 Here and later we use mostly monic polynomials, i.e., whose coefficient at the highest degree is 1. The d-orthogonal polynomials and the more general class of the so-called multiple orthogonal polynomials have been intensively studied in the last 30 years due to their intriguing properties and applications, cf., e.g., [5, 8, 18] and [40, Chapter 23] and further references in the cited literature. In particular they have applications to random matrices [44, 45], simultaneous Padé approximations [25], number theory [12, 17, 53] (which in fact go back to Hermite), etc. Notice that the classical orthogonal polynomials have a number of properties that are missing in general for the rest of the polynomial systems. Here we list some of them: • they are eigenfunctions of an ordinary differential operator, • they have explicit differential ladder operators (operators raising or lowering the index), • they can be presented in terms of hypergeometric functions, • they can be presented via Rodrigues formulas, • there are Pearson’s equations for the weights of their measures, • they possess the Hahn’s property, i.e., the polynomial system of their derivatives are again orthogonal polynomials. The class of polynomial systems that we introduced in [38] also shares these properties. Some of them were established in the cited paper, e.g., explicit differential ladder operators1 and Rodrigues-like formulas, apart of the differential equation. All these properties were direct consequences of our main construction. Other properties will be found here – their hypergeometric representations, generating functions and in some cases – Mehler–Heine formulas. In another project we intend to obtain also the weights, defining the functionals Lj , together with the Pearson equations for them and show their connection to biorthogonal ensembles, cf. [20, 45]. All these properties make them close analogs of the classical orthogonal polynomials. In fact there are other polynomial systems that are analogs of classical orthogonal polynomials. In [7, 60] the authors take another direction to generalize classical orthogonal polynomials. Namely they use the weights of the latter to produce a vector of weights. Their polynomial systems were later used in the study of different random matrix models – non-intersecting Brownian motion, matrix models with an external source, two-matrix models, cf. [6, 18, 44]. The main tools we used in [38] are the automorphisms of non-commutative algebras. We explain the construction for the case of the first Weyl algebra W1 . It can be realized as the algebra of differential operators with polynomial coefficients in one variable. W1 acts on the space of polynomials C[x]. We consider the simplest polynomial system {ψn (x) = xn , n = 0, 1, . . .}, and the differential operator H = x∂ which satisfies Hψn (x) = nψn (x), ∂ψn (x) = nψn−1 (x) and xψn (x) = ψn+1 (x). d . Below for simplicity we take q(∂) = −∂ l /l, Take any polynomial q(∂) in ∂, where ∂ := dx adq(∂) l ∈ N. It defines an automorphism σ = e of W1 , acting on A ∈ W1 as j σ(A) = e − ad∂ l /l A= ∞ ad (A) X −∂ l /l j=0 1 j! , As one of the referees kindly pointed to me in fact each polynomial system has ladder operators, cf. [13]. However in [38] we have given explicit differential ones. d-Orthogonal Analogs of Classical Orthogonal Polynomials 3 where adA (B) = [A, B]. It is easy to see that the sum is finite. The images of the above operators are σ(H) = H − ∂ l , σ(x) = x − ∂ l−1 . σ(∂) = ∂, If we define the polynomials Pn (x) = e −∂ l /l ψn (x) = ∞ X (−∂ l /l)j xn j! j=0 , and put L = σ(H) we easily see that LPn (x) = nPn (x), xPn (x) = Pn+1 (x) + n(n − 1) · · · (n − l + 2)Pn−l (x). The main point is that we constructed simultaneously the polynomial system {Pn (x)}, the differential operator L and the finite-term recurrence. Notice that for l = 2 these are the Hermite polynomials. For arbitrary l these are the Gould–Hopper polynomials [33]. See also the examples. The same procedure can be repeated with other algebras. We can take instead of ∂ an operator of the form G = R(H)∂, where R(H) is any polynomial in H. Then instead of W1 we can take the algebra spanned by G, H, x. The case of discrete orthogonal polynomials can be treated exactly in the same manner, realizing W1 by difference operators, and starting with ψ(x, n) = x(x − 1) · · · (x − n + 1), cf. the next section or [38]. In the present paper we study further the d-orthogonal polynomials with the Bochner property adding to the tools some new arguments which are not present in [38]. Let us first list the new properties which we discuss below. The most fundamental one, of which the rest are consequences, is the hypergeometric representation of the class of d-orthogonal polynomials with the Bochner property corresponding to the special case q(G) = ρGl , l ∈ N, ρ ∈ C, which we obtain both in the continuous and the discrete cases. Moreover, in both cases we provide two different representations in terms of hypergeometric functions. One of these representations has the advantage that the corresponding formula has the same values of hypergeometric parameters for all positive integers n. However, the other formula, in which the corresponding parameters depend on n (mod l) seems to be more useful in applications. The hypergeometric representations show that the class of d-orthogonal polynomials under consideration is very similar to the Gould–Hopper polynomials [33], which correspond to G = ∂x and can also be included as a special case of our construction. For this reason, we call the dorthogonal polynomials introduced above the generalized Gould–Hopper (GGH) polynomials. One immediate consequence of the above hypergeometric representations is that these special polynomial systems naturally split in l families, each originating from the initial differential (difference) operator G, exactly as in the case of the Hermite polynomials, which also is a particular case of our general construction. Recall that the latter are naturally subdivided into two families: the even-indexed H2n (x), n = 0, 1, 2 . . . and the odd-indexed ones H2n+1 (x), n = 0, 1, 2 . . .. (1/2) 2
(−1/2) 2
x (up to x and H2n (x) = xLn The corresponding representations are H2n (x) = Ln (α) multiplicative constants), where Ln are the generalized Laguerre polynomials. It is worth to point out that the above mentioned hypergeometric representations also use our main construction of the d-orthogonal polynomials in [38], but they require some additional transformations. The point is that in these representations of the polynomials Pn (x) each summand naturally corresponds to a summand in a generalized hypergeometric series. In the 4 E.I. Horozov differential case when one uses q(G) = G, the hypergeometric representations were known before, see [15]. Also the case of the Gould–Hopper polynomials is known [46]. However, for q(G) = ρGl , l > 1, our formulas are new. In the discrete case, the formulas for the families corresponding to G = a∆, a ∈ C and q(G) = Gl , l > 1 can be found in [16]. All other formulas, except for the Charlier and Meixner polynomials are new. Our next goal is to find the generating functions for the GGH polynomials. We present two such formulas, both based on the second hypergeometric representation. For the first family of generating functions, we use the well-known method to obtain the generating functions for the classical orthogonal polynomials. For the second family, we apply a formula due to Srivastava [55]. Some of the d-orthogonal polynomials are known to have a generating function. These are exactly the ones mentioned above and which can be found in [15]. Finally we find Mehler–Heine type asymptotic formulas [35, 48], which are again based on the second hypergeometric representation. In case of the multiple orthogonal polynomials, such formulas can be found in [57, 59]. In particular [59] contains our first result. We do not treat the discrete case, where as kindly pointed to me one of the referees the corresponding notion is local separate convergence, although some authors use again Mehler–Heine type asymptotics, cf. [28]. It is essential to mention that the constructions of all systems of d-orthogonal polynomials in [38] are based on the so-called bispectral problem, see [10, 11, 29]. Namely, these polynomials are the eigenfunctions for two linear operators – the first one is differential (difference) L(x) in the variable x, and the second one is the operator Λ(n) in the variable n, corresponding to the finite recurence relation: L(x)Pn (x) = nPn (x), xPn (x) = Λ(n)Pn (x). It turned out that both the generating functions and the Mehler–Heine asymptotics in the continuous case are expressed in terms of hypergeometric functions of the form   − ; xt . 0 Fq α 1 , α2 , . . . , α q The latter functions in the form of Meijer’s G-functions appear in an entirely different bispectral problem – for which both variables x and n are continuous, see [11], where they are called generalized Bessel functions. It seems that this is not a mere coincidence but could be exploited further. In particular, the already known results for the Bessel bispectral functions and their Darboux transformations could be used as a model in the study of our d-orthogonal polynomials with the Bochner property. Apart from Darboux transformations, some other possible directions of continuations of the present studies might include the so-called linearization problem (Clebsch–Gordan coefficients) and the problem of “connection coefficients” for the GGH polynomials. The asymptotic formulas from the present paper seem promising in the study of the zeros of the d-orthogonal polynomials with the Bochner property. Some other classical issues could be studied as well, such as finding the measures and the corresponding Pearson equations for them. In another project we intend to pursue the connections of the polynomial systems and the corresponding measures with integrable systems such as bigraded Toda hierarchy and KP-hierarchy. Some of the well known matrix models originate from the special solutions of these hierarchies, corresponding to these polynomial systems or their measures, e.g. the generalized Kontsevich–Penner model, Brezin–Gross–Witten model, etc., see [3, 4, 49] as well as the Kontsevich [43] model itself. We finally point out that some of these polynomial systems have been studied and used for many years. Apart from the Gould–Hopper polynomials we can mention the Konhauser–Toscano polynomials, whose special cases have been d-Orthogonal Analogs of Classical Orthogonal Polynomials 5 studied and used as early as in 1951 in connection with the penetration and diffusion of the Xrays, see [54]. These polynomials found more recent applications in the random matrix theory where they have been used in the so-called Muttalib–Borodin biorthogonal ensembles. The latter appeared in the studies of disordered conductors in the metallic regime cf. [20, 31, 50, 64]. Also some of the continuous families describe products of Ginibre random matrices, cf. [2, 45] (again by making use of the hypergeometric and Meijer’s G-functions representations). Based on the rich mathematical properties of the d-orthogonal polynomials from this paper and in [38] we hope that their study will be useful also in other problems. However here we do not pursue direct application but rather to give a unified approach to all these important special cases, scattered in the literature, see, e.g., [14, 15, 16, 32, 33, 63]. The paper is organized as follows. In order to make it independent of [38], in Section 2 we recall all definitions and statements needed in the main part of the paper. We also give a proof of Hahn’s property. Section 3 is the central one; here we derive different formulas for the GGH polynomials in terms of the hypergeometric functions, including some well-known. However our approach is novel; it is based entirely on the algebraic construction from [38] and it emphasizes the common origin of all formulas, old and new. Next, in Section 4 using the hypergeometric representations, we derive the generating functions for the GGH polynomials. In Section 5 we prove some results of the Mehler–Heine type asymptotics for GGH polynomials. We finish the paper with a number of examples, see Section 6. Besides pure illustrations some of them provide new interpretations, see, e.g., certain cases of “matching polynomials of graphs”. 2 Preliminaries To make the present paper self-contained below we briefly recall some of the required notions and results obtained in [38]. Given a field F of characteristic zero, consider the Weyl algebra W1 with coefficients in F spanned by two generators Y , Z subject to the relation [Z, Y ] = 1, where [Z, Y ] := ZY − Y Z is the standard commutator. Let us introduce some subalgebras of W1 . Set H = Y Z. Fixing a nonzero polynomial R(H), put G = R(H)Z. The first subalgebra B1 ⊂ W1 is, by definition, generated by the elements H, G, Y . These elements satisfy the following relations: [H, Y ] = Y, [H, G] = −G, [G, Y ] = R(H)(H + 1) − HR(H − 1). For any polynomial q(G) without a constant term2 , define the automorphism of B1 given by the operator σq = eadq(G) . The images of the generators of B1 under the automorphism σq are given below. Lemma 2.1. In the above notation, σq (G) = G, σq (H) = H + q ′ (G)G, σq Y = Y + ld+l−1 X γj (H)Gj . j=0 for some polynomials γj (H). To move further, we need to introduce an auxiliary algebra R2 over F defined by the generators T , T −1 , n̂ subject to the relations:  −1  T · T −1 = T −1 · T = 1, [T, n̂] = T, T , n̂ = −T −1 . 2 The constant term contributes only to multiplication of the polynomials, defined below in (2.1), by a constant. On the other hand, if it vanishes, the formulas and the arguments become simpler. 6 E.I. Horozov   One can easily check that T, n̂T −1 = 1 which implies that the operators T , n̂T −1 determine a realization of the Weyl algebra W1 . We can now introduce another non-commutative algebra B2 as follows. First we define an anti-homomorphism b, i.e., a map b : B1 → R2 satisfying b(m1 · m2 ) = b(m2 ) · b(m1 ), for each m1 , m2 ∈ B1 given by b(Y ) = T, b(G) = n̂T −1 R(n̂). b(H) = n̂, By definition, the algebra B2 is the image b(B1 ). With this definition, b : B1 → B2 is an anti-isomorphism and, in particular, b−1 : B2 → B1 is well defined. Observe that we can represent the algebra B1 on the space C[x] of polynomials in one variable by realizing Y as the operator x of multiplication by x, and Z as the operator of differentiation ∂x . Consider the polynomial system ψ(x, n) = xn . Then the action of the operators H, G on C[x] is given by H → x∂x , G → R(x∂x )∂x . In the same way, we can represent in C[x] the algebra B2 by realizing T and T −1 as the shift operators acting on a function f (n) as T ± f (n) = f (n ± 1). Finally, n̂ acts on C[x] by multiplication by the number n. Using this notation, we get Lemma 2.2. Gψ(x, n) = n̂R(n̂ − 1)T −1 ψ(x, n), Hψ(x, n) = n̂ψ(x, n), xψ(x, n) = T ψ(x, n). Furthermore, using the operator q(G), we can define another polynomial system {Pnq (x)} as Pnq (x) = eq(G) ψ(x, n) = ∞ X q(G)j ψ(x, n) j! j=0 . (2.1) Notice that, in fact, the above series is always finite since the operator q(G) reduces the degree of any polynomial it is applied to. Denote by L the operator σq (H) and put d = deg R, l = deg q. Theorem 2.3. In the above notation, the polynomials Pnq (x) have the following properties: (i) They are the eigenfunctions of the differential operator L := q ′ (G)G + x∂ with the eigenvalues λ(n) = n. (ii) They satisfy a recurrence relation of the form q (x) + xPnq (x) = Pn+1 ld+l−1 X q γj (n)Pn−j (x). j=0 (iii) They possess the Hahn’s property, i.e., their derivatives are of the same class with a new Ĝ = R(H + 1)∂x . Proof . Proof of (i) and (ii) can be found in [38]. The proof of (iii) is simple and goes as follows. We have ∂x Pn (x) = ∞ X ∂x q(G)j xn j=0 j! . Notice that ∂x G = ∂x R(H)∂x = R(H + 1)∂x2 . Hence ∂x q(G)j xn = nq(Ĝ)xn−1 . This shows ′ that the system Qn−1 (x) = Pn (x)/n, n = 1, 2, . . . is the system corresponding to the operator q(Ĝ).  d-Orthogonal Analogs of Classical Orthogonal Polynomials 7 Notice that, similarly to the above, we have earlier realized the abstract construction of B1 but in terms of difference operators instead of differential ones, acting again on the space of polynomials C[x] [38]. Then we define the operators τ± acting on f ∈ C[x] by the shift of the argument τ±1 f (x) = f (x ± 1). The operator x acts on f (x) as multiplication by x. We also need the notation ∆ = τ+1 − 1, ∇ = τ−1 − 1, H = −x∇. Finally, put g = x − H. For this realization we use the following polynomials system ψ(x, n) = (−1)n (−x)n = x(x − 1) · · · (x − n + 1). Here we use the Pochhammer symbol (a)j = a(a + 1) · · · (a + j − 1), (a)0 = 1. (2.2) Lemma 2.4. The following identities hold: Hψ(x, n) = nψ(x, n), gψ(x, n) = T ψ(x, n), ∆ψ(x, n) = nT −1 ψ(x, n). Exactly as in the continuous case we define the operator G = R(H)∆. Let q(G) be a polynomial without a constant term. Then we can define the automorphism σq = eadq(G) . Also introduce the new polynomial system {Png (x)} given by Pnq (x) = eq(G) ψ(x, n).  In [38] among other things, we proved the following properties of Pnq (x) . Theorem 2.5. The polynomials Pnq (x) have the following properties: (i) They are eigenfunctions of the difference operator L := q ′ (G)G − x∇. (ii) They satisfy a recurrence relation of the form q xPnq (x) = Pn+1 + m X q , γj (n)Pn−j j=0 where m = l for d = 0 and m = ld + l − 1 for d > 0. (iii) They possess the Hahn’s property, i.e., their polynomial system def Qn (x) = ∆Pn+1 , n+1 n = 0, 1, . . . , corresponds to q(Ĝ) with a new Ĝ = R(H + 1)∆. Proof . For (i) and (ii) see [38]. The proof of (iii) is similar to the one of Theorem 2.3, using the relation [∆, H] = ∆.  3 Hypergeometric representations In this section which is central for the present paper, we derive formulas for some of the families of d-orthogonal polynomials introduced in [38] and Section 2 in terms of the generalized hypergeometric functions. These families are defined via the operators q(G) = ρGl , l ≥ 1. This property is very important and, in particular, it makes our polynomial systems similar to  (α) the classical orthogonal polynomials. Notice that the system of Laguerre polynomials Ln (x) 8 E.I. Horozov corresponds to G = (x∂x + α + 1)∂x and l = 1, ρ = 1. The system of Hermite polynomials corresponds to ∂x and q(G) = −G2 /2. Let us recall the definition of generalized hypergeometric series, cf. [51]. Given a pair of nonnegative integers (p, q), let α1 , . . . , αp and β1 , . . . , βq be complex constants and let x be a complex variable. (The parameters βi are assumed to be different from non-positive integers.) The series   ∞ j def X (α1 )j · · · (αp )j x α1 , . . . , αp ; x := , (3.1) p Fq β1 , . . . , βq (β1 )j · · · (βq )j j! j=0 where (α)n is the Pochhammer symbol (2.2) (also the symbol of the rising factorial), is called a generalized hypergeometric series. When it is convergent in some open set its analytic continuation is called generalized hypergeometric function. This function satisfies the differential equation     q p  Y  Y α , . . . , αp ∂x ; x = 0. (x∂x + βj − 1) − (x∂x + αk ) p Fq 1 β1 , . . . , β q   j=1 k=1 We have to point out that the above function does not always exist. However, when one of the parameters αi is a non-positive integer, the series terminates and the function becomes a polynomial. In this paper we mainly deal with generalized hypergeometric polynomials. In the remaining cases the series will be obviously convergent, due to the fact that p < q in the generating functions and in the Mehler–Heine asymptotic formulas. (For more details consult [51].) From now on we will omit the word “generalized” since there will be no danger of confusion. Notice however that the expression “the hypergeometric function” usually refers to 2 F1 [ a,b c ; x]. 3.1 Notation In what follows we will introduce some notation that will be used to formulate and prove the results. First we will use the shorthand notation   (αp ) ;x p Fq (βq ) for the hypergeometric series (3.1). Let us introduce αd+1 = 0 in an effort to make the notation less clumsy. By (αd+1 ) we denote the vector of parameters (α1 , . . . , αd , αd+1 ). Recall that the notation ∆(l; λ) abbreviates the vector   λ+l−1 λ λ+1 , l ∈ N, , ,..., l l l of l parameters. We will also combine the latter symbol with (αd+1 ) to write (∆(l; αd+1 )) for α1 α1 + l − 1 αd+1 αd+1 + l − 1 ,..., ,..., ,..., , l l l l in the expressions for the hypergeometric functions. We assume that ρ 6= 0. Let us put η = ld+1 ρ. Let us fix i ∈ {0, 1, . . . , l − 1} and consider the polynomials with n = ml + i, m = 0, 1, . . .. Then the parameters ∆(l; −n − αd+1 ) can be presented as follows  i + α1 i − 1 + α1 i + α1 − l + 1 ∆(l; −n − αd+1 ) = −m − , −m − , . . . , −m − ,..., l l l d-Orthogonal Analogs of Classical Orthogonal Polynomials 9 i + αd+1 i − 1 + αd+1 i + αd+1 − l + 1 −m− , −m − , . . . , −m − l l l Let us introduce the sets   αk + i − r + 1, r = 0, . . . , l − 1 , Sk (i) = l  . k = 1, . . . , d + 1 and S(i) = ∪d+1 k=1 Sk (i). By Ŝ(i) we denote the set S(i) \ {1}, 1 is the element, corresponding to k = d + 1, r = i. The elements of S(i) will be denoted by Sβ (i), β = 1, . . . , (d + 1)l. By I we denote the set of indexes β, which correspond to Ŝ(i). Notice that I depends on i. We also introduce a notation for the product of Pochhammer symbols: def [αp ]n = p Y (αk )n . k=1 Most of the notation is taken from [40, 41, 46, 56]. 3.2 Continuous GGH polynomials Some of the polynomial systems of the previous section have well known representations in terms of generalized hypergeometric functions [15], where the authors take the formulas as their definition. Below we will derive the representation from [15] on the basis of the constructions from [38]. Our proof will help us to find hypergeometric representations for the other systems studied here, that are not treated elsewhere. The discrete d-orthogonal polynomials will be treated in the same manner. We start with the case when the automorphism σ is defined by G = R(H)∂x where R(H) is m Q (H + αj + 1), a polynomial. Before that let us factor the polynomial R(H) into R(H) = ρ j=1 ρ ∈ C, ρ 6= 0. Here some of the numbers αj can be equal. In what follows we will drop the dependence of the d-orthogonal polynomials on R. Theorem 3.1. The polynomials Pn (x) have the following hypergeometric representations:   −n, (−n − αd ) d+1 −1 n , (3.2) ; (−1) ρx Pn (x) = x d+1 F0 −   x −n n (3.3) ;− . Pn (x) = ρ [αk + 1]n 1 Fd (αd + 1) ρ For the second formula the coefficients αj have to be different from negative integers or zero. Proof . Notice that Gxn = nR(H)xn−1 = nR(n − 1)xn−1 . By induction we find that for any j ∈ N j n G x = j−1 Y (n − s)R(n − 1 − s)xn−j . s=0 In terms of the Pochhammer symbol we can write the last formula as Gj xn = (−1)j(d+1) (−n)j ρj [−n − αd ]j xn−j . 10 E.I. Horozov The last formula allows to express the polynomials Pn (x) as Pn (x) = xn n X (−1)j(d+1) (−n)j [−n − αd ]j ρj x−j j! j=0 =x n d+1 F0   −n, (−n − αd ) d+1 −1 . ; (−1) ρx − (3.4) To obtain the second expression for the polynomials Pn (x), we transform the coefficients in the first line of (3.4) using the formula (−n − µ)j = (−1)j (µ + 1)n . (µ + 1)n−j (3.5) Substitute this expression for the parameter values µ = 0, α1 , . . . , αd into the sum (3.4), defining Pn (x). As a result, we obtain n Pn (x) = ρ [αd + 1]n n X n! j=0 1 ρ−n+j xn−j . j! [αd + 1]n−j (n − j)! Write the expression n!/j! as (−1)n−j (−n)n−j . After changing the summation index j → s = n − j, we find n Pn (x) = ρ [αd + 1]n 1 Fd   x −n ;− . (αd + 1) ρ Of course we need to impose the condition αk 6= −1, −2, . . ..  Remark 3.2. (1) This expression coincides with the corresponding one in [15]. The form of the roots of R was chosen for this purpose as well as to obtain the hypergeometric formula of the (α) generalized Laguerre polynomials Ln (x). (2) Formula (3.4) is valid for all values of the coefficients αj . However if some of the coefficients αj is are negative integers, the polynomial system becomes xn for n ≥ −αj and therefore in such a case it forms a finite system of d-orthogonal polynomials. Let us consider the polynomials obtained by the automorphisms σq , where q(G) is some polynomial without constant term. We recall that they are given by Pnq (x) := e q(G) n n x =x + ∞ k X q (G)xn k=1 k! . I don’t know if all these d-orthogonal polynomials have representations in terms of generalized hypergeometric functions or other special functions. However in the case when q(G) = ρGl they have. We are going to present the formula. Let us fix q(G) = ρGl . Finally we will drop the dependence on q(G) as there is no danger of confusion. Theorem 3.3. For q(G) = ρGl the polynomials Pn (x) have the following representation in hypergeometric functions: "  # d+1 η l (−1) (∆(l; −n − α )) d Pn (x) = xn dl+l F0 ; . (3.6) − x d-Orthogonal Analogs of Classical Orthogonal Polynomials 11 Proof . Again we use the formula for Pn (x): Pn (x) = ∞ X Glj xn j! j=0 . We know that Glj xn = lj−1 Y (n − s)R(n − 1 − s)xn−lj . s=0 Using that αd+1 = 0 we present the coefficient at xn−lj of Glj xn in the form  l−1  Y −n − αd+1 + r lj lj(d+1) . η (−1) l j (3.7) r=0 Then the polynomials Pn (x) are given by n Pn (x) = x + n X (Gl )j xn j=1 = xn n Y l−1  X j=0 r=0 j! =x n n X lj η (−1) lj(d+1) r=0 j=1 −n − αd+1 + r l  l−1  Y −n − αd+1 + r x−lj l j j!   j (−1)d+1 η x lj 1 . j! This proves (3.6).  Remark 3.4. The formula for the d-orthogonal polynomials, corresponding to q(G) = ρGl , in the above theorem resembles the representation of Hermite polynomials in hypergeometric   functions, see, e.g., [51]. Notice that we don’t need to sum up to n but only up to nl as the next terms are 0. As in the case l = 1 we are going to find a second representation. Theorem 3.5. The polynomials Pn (x), where n = ml + i, have the following hypergeometric representation "  l # Y
x −m ml i
; − Ŝβ (i) m 1 Fld+l Pn (x) = η x . (3.8) Ŝ(i) η β∈I Proof . Notice that (3.6) can be written in the form " #  d+1 η l (−1) (−S(i) + 1 − m) . Pn (x) = xml+i dl+l F0 ; − x
l (i) (i) Let us put y = xη and consider the polynomials Pm (y) defined by Pn (x) = xi η ml Pm (y). We use that one of the elements of the set S(i) is 1, which shows that they have the representation (3.3) from Theorem 3.1   −m, (−Ŝ(i) + 1 − m) (−1)ld+l m (i) . ; Pm (y) = y dl+l F0 − y They are exactly in the form of (3.2). This means that they can be written as   Y −m (i) ml
; −y . Pn (y) = η (Sβ (i))m 1 Fld+l−1 Ŝ(i) β∈I Returning to the polynomials Pn (x) we obtain the formula (3.8).  12 E.I. Horozov Corollary 3.6. The polynomials Pn (x) have also the representation Pn (x) = η ml Y (Sβ (i))m 2 Fld+l−1 β∈I "  l # x −m, 1 ;− . (S(i)) η Proof . Use that S(i) = Ŝ(i) ∪ {1}. Then the general formula p Fq     (αp ), 1 (αp ) ;x ; x = p+1 Fq+1 (βq ), 1 (βq ) applied to (3.8) gives the result.  We see that the GGH polynomial system Pn (x) can be split into l families of d-orthogonal polynomials with Bochner’s property exactly as Hermite polynomials split into 2 families in terms of Laguerre polynomials. More precisely (i) Corollary 3.7. The polynomial system Pn (x) consists of l families Pn (x) = Ci xi Pm (y), i = 0, . . . , l − 1, where  l x , y= η Ci = η ml Y (Sβ (i))m . β∈I (i) The family Pm (x) Q is the system of d-orthogonal polynomials with Bochner’s property corre(H + Sβ )∂. sponding to G = β∈I While the above statement does not need a proof it is worth to point out that it is a precise formulation of Remark 7.4. in [38]. In [38, Example 7.3] we demonstrated that such a splitting of the family exists for the simplest case of q(G) = G2 , where G = (x∂x + α + 1)∂x without using hypergeometric functions. 3.3 Discrete GGH Polynomials The discrete polynomial systems also have representation in terms of generalized hypergeometric functions. Below we will derive them using the approach that was exploited for continuous dorthogonal polynomials. These representations are new except for the cases of Charlier and Meixner polynomials. m Q We again factor the polynomial R(H) into R(H) = ρ (H + αj + 1). j=1 Theorem 3.8. The discrete d-orthogonal polynomials Pn (x) have the following hypergeometric representations: n Pn (x) = (−1) (−x)n d+1 F1   −n, (−n − αd ) d+1 ; (−1) ρ , x−n+1 (3.9) and Pn (x) = [αd + 1]n 2 Fd   −n, −x ;ρ . (αd + 1) (3.10) d-Orthogonal Analogs of Classical Orthogonal Polynomials 13 Proof . First we transform the formula for Pn (x) to get rid of the difference operators. We are going to exploit again the formula Pn (x) = ∞ X Gj ψ(x, n) j=0 ψ(x, n) = (−1)n (−x)n . with j! Using that Gψ(x, n) = nR(n − 1)ψ(x, n − 1), we obtain Gj ψ(x, n) = ρj (−1)n−j (−x)n−j (−1)j (−n)j (−1)j [−αd − n]j . Let us use the formula (−1)n−j (−x)n−j = (−x)n (−1)n . (x − n + 1)j (3.11) We obtain that Gj ψ(x, n) = (−1)n (−x)n ρj (−1)dj+j (−n)j [−αd − n]j . (x − n + 1)j For the polynomials Pn (x) this gives the expression n Pn (x) = (−1) (−x)n n X ρj (−1)(d+1)j (−n)j [−αd − n]j (x − n + 1)j j=0 (x − n + 1)j j!, which is (3.9). The second hypergeometric expression can be obtained as in the continuous case. Using again (3.5) we find Pn (x) = m Y (αk + 1)n n X j=0 k=1 n! (−1)n−j ρj (−x)n−j (x − n + 1)j . (n − j)! j![−αd − n]n−j After standard manipulations that we used for the continuous d-orthogonal polynomials we  come to the second formula (3.10). We can obtain hypergeometric representations for the class of discrete d-orthogonal polynomials, constructed via q(G), where G = R(H)∆. Again we will treat the case when q(G) = ρGl . Let us put η1 = (−1)(d+1)l η . l Theorem 3.9. In the case when q(G) = ρGl the polynomials Pn (x) have the following representation in hypergeometric functions  ∆(l; −n), (∆(l; −n − α)) l ; η1 . Pn (x) = (−1) (−x)n dl+l Fl ∆(l; x − n + 1) n  (3.12) 14 E.I. Horozov Proof . The proof needs a few changes in comparison with the continuous case but otherwise is straightforward. We are going to use again the series defining the polynomials Pn (x) = ∞ X Gjl ψ(x, n) j! j=0 . A formula, similar to the one in the continuous case holds Glj (−1)n (−x)n = lj−1 Y (n − s)R(n − 1 − s)(−x)n−lj (−1)n−lj . s=0 For the term (−x)n−lj we use (3.11) with lj instead of j (−1)n−lj (−x)n−lj = (−1)n (−x)n . (x − n + 1)lj The factor (x − n + 1)lj can be presented as indicated in Section 3.2 (x − n + 1)lj = l lj  l−1  Y x−n+r+1 l r=0 . j This gives (−1)n−lj (−x)n−lj = llj (−1)n (−x)n . l−1 Q x−n+r+1
r=0 Using (3.7) for lj−1 Q l (3.13) j (n − s)R(n − 1 − s) and expressing the last factor by (3.13) in the sum for s=0 the polynomials we obtain the desired formula (3.12).  Again we can find a second formula for the d-orthogonal polynomials, corresponding to q(G) = ρGl . We use the notations (3.2) from Section 3.1. With this notation we have Theorem 3.10. The polynomials Pn (x), where n = ml + i have the following hypergeometric representation Pml+i (x) = (−1) mdl+1 C(i)(−x)i ld+l−1 Y k=1 where C(i) was defined in Corollary 3.7. l Fld+l−1   −m, ∆(l; −x + i − 1) −l
; −η1 , Ŝ(i) Proof . We start with the formula Pn (x) = ∞ X j=0 lj−1 Q s=0 (n − s)R(n − 1 − s)(−1)n−lj (−x)n−lj We transform the expression lj−1 Y s=0 (n − s)R(n − 1 − s) j! . (3.14) d-Orthogonal Analogs of Classical Orthogonal Polynomials 15 using (3.7) to obtain (−1)lj η lj (−m)j Y (−m − Sβ (i))j . β∈I We further transform the last expression using (3.5) into Q (Sβ (i))m j β∈I l(d+1)j lj (−1) m! Q . (−1) η (m − j)! (Sβ (i))m−j (3.15) β∈I Next we transform the factor (−x)n−lj as follows. First (−1)n−jl (−x)n−lj = (−1)i (−x)i (−1)ml−jl (−x + i)lm−lj . Then we present the last factor in the right-hand side of the last formula as  l−1  Y −x + i − r (m−j)l (−x + i)lm−lj = l . l m−j r=0 At the end plugging the last formula and (3.15) into the sum for Pn (x) and changing the summation index j → m − j we obtain Pml+i = C(i)(−x)i m X j=0 which is (3.14). 4 (−m)j l−1 Q r=0 Q −x+i−r l j (Sβ )j β∈I
 j −η1−l , j!  Generating functions In this section we will find generating functions for all d-orthogonal polynomials for which we have earlier obtained hypergeometric representations. Generating functions for the continuous d-orthogonal polynomials, corresponding to q(G) = G, can be found in [15]. 4.1 Continuous GGH Polynomials In what follows we use a different normalization of the d-orthogonal polynomials. Namely, we introduce the polynomials " l #  x −m
; . (4.1) Qn (x) := xi 1 Fld+l−1 Ŝ(i) (−1)d+1 ηρ They differ from the polynomials Pn (x), given by (3.8) by a multiplicative constant. We also assume that ((−1)d+1 η)l = −1 which can be achieved by rescaling of x together with multiplication by a suitable constant. Our first result is as follows. Theorem 4.1. For a given positive integer i, the function Φi (x, t) defined as   − l i tl
; (xt) , Φi (x, t) := (tx) e 0 Fld+l−1 Ŝ(i) generates the polynomial system {Qlm+i (x)}∞ m=0 by means of the formula ∞ ml+i X t Qml+i (x). Φi (x, t) = m! m=0 (4.2) 16 E.I. Horozov Proof . Consider the series   ∞ ml+i ∞ ml X X t t −m l i i
; −x . Qml+i (x) = t x 1 Fld+l−1 Ŝ(i) m! m! m=0 m=0 One can present the sum in the right-hand side of the latter equation in the form of a double series
j m ∞ ml+i ∞ ml X X X −xl (−m)j t t Q Qml+i (x) = . m! m! (Sβ (i))j j! m=0 m=0 j=0 β∈I Changing the order of summation in the double series, we obtain
j ∞ ∞ X X −xl (−m)j tml Q . (Sβ (i))j j! m! j=0 m=j β∈I It is easy to see that (−m)j (−1)j = . m! (m − j)! Hence after introducing a new index s = m − j of summation, the double series becomes ∞ X j=0 ∞ Q β∈I X tsl (xt)lj . (Sβ (i))j j! s! s=0 This gives for the double sum   l − l
; (xt) et , 0 Fld+l−1 Ŝ(i) which implies (4.2) for the function Φi (x, t).  Remark 4.2. We notice that the hypergeometric functions without upper parameters appear in a quite different bispectral problem. In [11] we defined the functions Ψβ (x, z), which are solutions of an equation of the form x−N (θ − β1 ) · · · (θ − βN )Ψβ (x, z) = z N Ψβ (x, z), where θ = x∂x and βj ∈ C. We called them generalized Bessel functions. As one of the referees kindly informed me these functions have been studied by P. Delerue in [26], and are called hyperBessel functions. They are expressed in terms of the hypergeometric functions without upper parameters (in [11] we used Meijer’s G-functions). Through these functions we were able to find non-trivial bispectral operators of any rank. It is interesting to understand if this is a mere coincidence or the reasons are deeper. I hope that such a connection exists and in this case it would be useful in the studies of Darboux transformation of the generalized Gould–Hopper polynomials. Even in the case of l = 1 (for which the same formula was found by different methods in [15], see also below Corollary 4.5) the connection deserves attention. From the above formulas we can write a generating function for the entire family Qn . Let us define the function   l−1 l−1 X X − l tl i
Φ(x, t) = ; (xt) . Φi (x, t) = e (xt) · 0 Fld+l−1 Ŝ(i) i=0 i=0 d-Orthogonal Analogs of Classical Orthogonal Polynomials 17 Corollary 4.3. The function Φ(x, t) is a generating function for the polynomials Qn (x)/⌈n/l⌉!: ∞ X tn Qn (x) = Φ(x, t). ⌈n/l⌉! n=0 The proof is obvious and we omit it. Remark 4.4. Polynomial systems that have a generating function of the form ∞ X Pn (x)tn = A(t)B(xt). n=0 are called Brenke polynomials, see, e.g., [14]. The last corollary shows that the GGH polynomials are Brenke polynomials. It deserves to write separately the formula for the case l = 1. Corollary 4.5 ([15]). When l = 1 we have ∞ n X t n=0 n! Qn (x) = 0 Fd   − ; xt et . (αd + 1) We are going to obtain a second formula based on a theorem from [55]. Let us formulate the corresponding result explicitly in a slightly less general form that suffices for our purposes. Proposition 4.6. Let a ∈ C, −a ∈ / N and let α1 , . . . , αp , β1 , . . . , βq be complex numbers such that the hypergeometric function   −n, (αp ) ;x p+1 Fq+1 a + 1, (βq ) be well defined. Then the following formula holds  ∞  X a+n n n=0 p+1 Fq+1     1 (αp ) xt −n, (αp ) n ;x t = ; . p Fq a + 1, (βq ) (βq ) 1 − t (1 − t)a+1 (4.3) See [55] for a simple proof. We again use the d-orthogonal polynomials Qn (x) from (4.1) as well as the convention η = 1. Theorem 4.7. The function Gi given by   (xt)l (tx)i 1
; , Gi (x, t) = 1 Fld+l−1 Ŝ(i) 1 − tl 1 − tl generates the polynomials Qlm+i (x), m = 0, 1, . . . in the following way Gi (x, t) = ∞ X tml+i Qml+i (x). m=0 Proof . Let us multiply the polynomials Qlm+i by tml+i and sum. Thus we obtain ∞ X m=0 t ml+i i i Qml+i (x) = t x ∞ X m=0 t ml 2 Fld+l   −m, 1
; xl . Ŝ(i) , 1 18 E.I. Horozov Here we have used that     (αp ), µ (αp ) ;x , ; x = p+1 Fq+1 p Fq (βq ), µ (βq ) µ 6= 0. We apply the above cited formula (4.3) from [55] with a = 0 to obtain ∞ X t m=0 ml+i (tx)i Qml+i (x) = 1 Fld+l−1 1 − tl   (xt)l 1
; . Ŝ(i) 1 − tl  From this theorem we can write a generating function for all polynomials Qn (x). Corollary 4.8. A generating function for Qn (x) is given by ∞ X tn Qn (x) = n=0   l−1 X (tx)i (xt)l 1
; . F 1 ld+l−1 Ŝ(i) 1 − tl 1 − tl i=0 Proof . Just sum up the generating functions Φi (x, t) for i = 0, . . . , l − 1 and replace n by ml + i.  Notice that the coefficients Ŝ(i) depend on i, which makes it difficult to obtain a better formula in general. However when l = 1 the above expression gives Corollary 4.9. The d-orthogonal polynomials defined in terms of q(G) = ρG have the following generating function   ∞ X 1 xt 1 n . t Pn (x) = ; 1 Fd (αd + 1) 1 − t 1−t n=0 4.2 Discrete GGH polynomials In the discrete case there is nothing special. We follow the arguments for the continuous case. However we will keep the coefficient ρ. We again put n = ml + i and fix i. We use the following modification of the polynomials (3.14)   −m, ∆(l; −x + i − 1) l i
; η1 , Qn (x) = (−1) (−x)i 1+l Fld+l−1 Ŝ(i) which differ from Pn (x) only by a multiplicative constant. Theorem 4.10. The function Φi given by   ∆(l; −x + i) l j i tl
; −(η1 t) . Φi (x, t) = (−1) (−x)i t e l Fld+l−1 Ŝ(i) is a generating function for the polynomials Qlm+i (y) in the sense that ∞ ml+i X t Qml+i (x) = Φi (x, t). m! m=0 Proof . We write the defining series as (−1)i (−x)i ∞ ml+i X m X (−m)j (∆(l; −x + i))j η1jl t Q . m! (Sβ (i))j j! m=0 j=0 β∈I d-Orthogonal Analogs of Classical Orthogonal Polynomials 19 In the right-hand side we make the following transformations. We first change the order of the summation and then introduce a new summation variable m → s = m − j. We obtain  m ∞ X (∆(l; −x + i))j η1l ]j X (−s − j)j t(s+j)l i i Q . (−1) (−x)i t (Sβ (i))j j! (s + j)! j=0 s=0 β∈I Notice that (−s − j)j (−1)j = . (s + j)! s! This gives Φi (x, t) = (−1)i (−x)i ti m ∞ X η lj (−1)j tlj X tsl Q1 , (Sβ (i))j j! s! j=0 s=0 β∈I which is exactly the desired formula.  As in the continuous case we are going to derive a second formula. Theorem 4.11. The function Gi given by " l #  (−1)d ti (−1)i (−x)i ∆(l; −x + i), 1
;− Gi (x, t) = 1+l Fld+l−1 Ŝ(i) 1−t ρ generates the polynomials Qlm+i (y) in the sense that Gi (x, t) = ∞ X tm Qml+i (x). m=0 Proof . We are going to use again Proposition 4.6. It is obvious that we need to put a = 0. We have   ∞ ∞ X X −m, ∆(l; −x + i − 1), 1 l m i m
t Qml+i (x) = (−1) (−x)i ; η1 . t 2+l Fld+l Ŝ(i) , 1 m=0 m=0 Then (4.3) gives ∞ X  −m, ∆(l; −x + i − 1), 1 l
; η1 t (−1) (−x)i 2+l Fld+l Ŝ(i) , 1 m=0   (−1)i (−x)i ti ∆(l; −x + i − 1), 1 η1l
= . ; 1+l Fld+l−1 Ŝ(i) 1 − tl 1 − tl i ml+i   As a trivial corollary again we can write a generating function for all polynomials Qn (x). Corollary 4.12. A generating function for Qn (x) is given by   ∞ l−1 X X (−1)i (−x)i ti ∆(l; −x + i − 1), 1 η1l n
t Qn (x) = . ; 1+l Fld+l−1 Ŝ(i) 1 − tl 1 − tl n=0 i=0 Finally for l = 1 we get a better formula Corollary 4.13. The polynomials defined in terms of q(G) = ρG have the following generating function   ∞ X tρ 1 x, 1 n ; . t Qn (x) = 2 Fd (αd + 1) 1 − t 1−t n=0 20 5 E.I. Horozov Mehler–Heine type formulas In this section, we are going to consider only the continuous d-orthogonal polynomials and, without loss of generality, we will assume that ρld+1 = 1. (This assumption will make the formulas simpler.) There are many ways to write down the Mehler–Heine type formulas depending on the normalization of the polynomials of which we choose only the simplest one. Let us start with the case of continuous d-orthogonal polynomials corresponding to q(G) = G and use the polynomials Qn (x) given by (4.1). Theorem 5.1. For the d-orthogonal polynomials obtained from q(G) = G, the following Mehler– Heine type formula holds   − ;x . (5.1) lim Qn (x/n) = 0 Fd n→∞ (αd + 1) Proof . Notice that in the case l = 1, equation (4.3) can be written as   −n ; −x . Qn (x) = 1 Fd (αd + 1) Then we use the formula     (αp ), aλ x (αp ) ; lim p+1 Fq = p Fq ; ax , (βq ) λ (βq ) λ→∞ with λ = n, a = −1, see, e.g., [41, p. 5]. This immediately gives (5.1).  Remark 5.2. This theorem is proved in [59]. We present it here as an illustration of the results of the present paper. Also the Mehler–Heine formula for general GGP follows from it. As one of the referees kindly pointed to me, the polynomials Qn (x) considered here all are ∞ P n Jensen polynomials. This means that they are associated to an entire function ϕ(x) = γn xn! n=0 in the following way Qn (x) = n   X n j=0 j γ j xj . Also the function et ϕ(xt) is their generating function: et ϕ(xt) = ∞ X n=0 Qn (x) tn , n! see, e.g., [23] for a contemporary reference to properties of Jensen polynomials that we refer to here. In our case   − ;x . ϕ(x) = 0 Fd (αd + 1) By the properties of the Jensen polynomials we have   − ;x lim Qn (x/n) = ϕ(x) = 0 Fd n→∞ (αd + 1) locally uniformly as was found by Jensen himself in 1913. d-Orthogonal Analogs of Classical Orthogonal Polynomials 21 Observe that the hypergeometric function   − ;x ϕ(x) = 0 Fd (αd + 1) appears once again naturally as in the formulas for the generating functions for the polynomial system {Qn (x)}. In both cases this is connected to the fact that they are Jensen polynomials. Now consider the polynomials obtained by the automorphisms σq , where q(G) is some polynomial. We recall that they are given by Pn (x) := e q(G) n n x =x + ∞ k X q (G)xn k! k=1 . We will restrict ourselves to the case q(G) = Gl for which the corresponding hypergeometric representation was obtained in Section 3. Presenting n as n = ml + i and using (4.1) for the polynomials Qn (x), we get i Qn (x) = x 1 Fld+l−1   l −m  d+1
; (−1) x . Ŝ(i) Theorem 5.3. For the d-orthogonal polynomials obtained via the automorphisms σq with q(G) = Gl , the following Mehler–Heine type formula holds i/l lim m Qn x/m m→∞ 1/l
i = x 0 Fld+l−1    − 
; (−1)d+1 x l . Ŝ(i) Proof . We
will use the hypergeometric representation (3.8). Qn x/m1/l which gives 1/l Qml+i x/m
xi = i/l 1 Fld+l−1 m " (5.2) Consider the expression  l # −m (−1)d+1 x
; . Ŝ(i) m The latter formula can be rewritten in the form "  l #
−m (−1)d+1 x 1/l i i/l
; . = x 1 Fld+l−1 m Qml+i x/m Ŝ(i) m Then formula (5.2) follows from (5.1).  It is worth noticing that the asymptotics depends on the remainder n (mod l) which indicates that there is probably no general asymptotic formula, but it might exist for the subsequence with the same value of the remainder n (mod l). This phenomenon is well-known in the case of Hermite polynomials, where the even-indexed and the odd-indexed polynomials have different asymptotics, see, e.g., [1]. Notice that as we explained in Remark 4.2 the function 0 Fld+l−1   − 
; (−1)d+1 x l Ŝ(i)  is also a generalized Bessel function in the sense of [11]. 22 6 E.I. Horozov Examples Example 6.1 (Gould–Hopper polynomials). Consider the simplest case R(H) ∈ C, i.e., G = ∂. Then for q = τ Gl , using equation (3.6) we obtain the polynomial system " # −n −n + l − 1  l −l , . . . , Pn (x) = xn l F0 l ;τ . l x − These polynomials coincide with the well-known Gould–Hopper polynomials gnl (x, τ ), cf. [33, 46]. According to our scheme they are the eigenfunctions of the differential operator L = lτ ∂ l + x∂ and satisfy the recurrence relation l l xgnl (x, τ ) = gn+1 (x, τ ) − τ ln(n − 1) · · · (n − l + 2)gn−l+1 (x, τ ). Using the second form of hypergeometric representation (3.8) we can also express them as   x l  −m l i
; −τ gml+i (x) = x 1 Fl−1 . Ŝ(i) l Observe that for l = 2, these polynomials coincide (up to rescaling) with the classical Hermite polynomials. Eventually the Gould–Hopper polynomials turned out to be quite useful in quantum mechanics, integrable systems (Novikov–Vesselov equation), combinatorics, etc., see, e.g., [21, 24, 62]. The cases with G = R(H)∂ with arbitrary R can be considered as generalizations of the Gould–Hopper polynomials. In this situation we use q(G) = Gl , where G = R(H)∂, R being a polynomial of an arbitrary degree. The corresponding expression for these polynomials provided by (3.6) is as follows "  l # ∆(l; −n), (∆(l; −n − α)) (−l)d+1 ρ n Pn (x) = x dl+l F0 ; . − x We can explicitly write discrete analogs of the (generalized) Gould–Hopper polynomials. Namely,  
l ∆(l; −n), (∆(l; −n − α)) ; (−1)(d+1) ld ρ . Pn (x) = (x)n dl+l Fl ∆(l; x − n) The most straightforward analog, which one might call the discrete Gould–Hopper polynomials, corresponds to G = τ ∆l (d = 0). These polynomials have the following hypergeometric representation   ∆(l; −n) l ; (−lρ) . Pn (x) = (x)n l Fl ∆(l; x − n) Having in mind the existing applications of the Gould–Hopper polynomials it is worth checking whether their generalized versions have similar or other applications. Example 6.2 (Konhauser–Toscano polynomials). In [42] Konhauser has defined two families of polynomials denoted by Ynα (x; l) and Znα (x; l), n = 0, . . ., where α ∈ R, l ∈ N. The polynomials Znα (x; l) are in fact polynomials in xl . The polynomials Ynα (x; l) are polynomials in the original variable x. d-Orthogonal Analogs of Classical Orthogonal Polynomials 23 These two families are biorthogonal with respect to the weight function corresponding to the Laguerre polynomials, i.e., Z ∞ α (x; l)dx = hm δn,m xα e−x Ynα (x; l)Zm 0 with hm 6= 0. The polynomials Znα (x; l) were introduced earlier by Toscano, see [58]. Their hypergeometric representation     α + ln (ln)! −n l α Zn (x; l) = ; (x/l) , 1 Fl ∆(l; α + 1) n! ln found in [46] shows that they can be constructed using the methods of the present paper. Let l
Q . Then H + α+s us consider G = R(H)∂, withl R(H) = l s=1 Znα (x; l) = Pn (x/l) .
l In the case l = 2 the polynomials were discovered much earlier by L.V. Spencer and U. Fano [54] in their studies of the X-rays diffusion. Remark 6.3. In [42] Konhauser has proven that the polynomials Ynα (x; l) are the eigenfunctions of a differential operator of order l + 1 and that they satisfy a (l + 2)-recurrence relation of the form x l Ynα (x; l) = l X α (x; l). γj (n)Yn+j j=−1 The polynomials Znα (x; l) satisfy a relation of the form x l Znα (x; l) = l X α βj (n)Zn−j (x; l). j=−1
(This relation follows from the definition Zn (x; l) := PnR (x/l)l and the properties of
PnR (x/l)l .) As mentioned earlier the Konhauser–Toscano polynomials have applications to random matrix theory. The so-called Borodin–Muttalib ensembles [20, 50] are based on their biorthogonality. In fact the polynomials Ynα (x; l), n = 0, 1, . . . , l − 1 are closely related with the functionals, which define the d-orthogonal polynomials Pn (x). It is quite interesting to define discrete analogs for the above Konhauser polynomials. Example 6.4 (matching polynomials of graphs). Polynomial systems considered in this example are taken from [9] and they are relevant for the so-called chemical graph theory. We need some notions from the graph theory, see, e.g., [9, 27, 30] and the references therein. Let K be a connected graph with n vertices. Following [52] we define the higher Hosoya number pr (K, j) as the number of ways one can select j non-incident paths of length r in K. Using the Hosoya numbers, the higher-order matching polynomial Mr (K) of K is defined by the relation n ⌈ r+1 ⌉ Mr (K) := X j=0 (−1)j pr (K, j)xn−(r+1)j , 24 E.I. Horozov see [9, 30]. When K = Kn is the complete graph on n vertices (i.e., each pair of vertices is connected by an edge) the corresponding polynomials were explicitly computed in [9, 30]. Using combinatorial arguments it was shown that these polynomials are given by Mr (Kn ) = n X (−1)j j=0 n!xn−(r+1)j . (n − (r + 1)j)!j!2j If we compare this expression with Pn = e−G r+1 /2 xn , where G = ∂x , we see that they coincide. Hence we obtain a hypergeometric representation   ∆(r + 1; −n) (−1)r (r + 1)r+1 n . ; Mr (Kn ) = x r+1 F0 − 2xr+1 (This representation was earlier found in [9].) We see that the matching polynomials of complete graphs are the eigenfunctions of a linear differential operator and that they satisfy an (r+2)-term recurrence relation which can be useful in the studies of these polynomials. The case r = 1 deserves a special attention since it corresponds to the Hermite polynomials as was observed long ago, e.g., in [9, 34]. The corresponding coefficients p1 (K, j) are the original Hosoya numbers [39]. Let us consider another example of matching polynomials, this time of complete bipartite graphs Kn,m , n ≥ m with n + m vertices. (Recall that Kn,m is a graph whose vertices are split into two nonintersecting sets Vn and Vm with n and m elements resp. and every vertex in Vn is connected to every vertex in Vm . We consider the case when r is odd. Then [9, 30] contain the formula n    X n m (−1)j ((jr)!)2 n−2rj Mr (Kn,m ) = x , j! jr jr j=0 which can easily be transformed into   ∆(r + 1; −n), ∆(r + 1; −m) −(r + 1)r+1 n+m ; Mr (Kn,m ) = x . r+1 F0 − (2x)r+1 Set m = n − M , M ≥ 0. Using this notation, we see that   ∆(r + 1; −n), ∆(r + 1; −n + M ) −(r + 1)r+1 2n−M . ; Mr (Kn,m ) = x r+1 F0 − (2x)r+1 In other words, Mr (Kn,m ) coincide with xn−M Pn (x), where the polynomials Pn (x) are constructed via q(G) = −Gr+1 /2 and G = (x∂ − M )∂. While the hypergeometric representation is known, the properties of these polynomials listed in Theorem 2.3 are new. In particular, the differential equation and the recurrence relations they satisfy seem to be new. Notice that for r = 1 (i.e., when all the paths are edges), these polynomials coincide with (M ) 2
(α) Ln x , where Ln (y) are the generalized Laguerre polynomials. The above result suggests a conjecture about the matching polynomials of complete k-partite graphs, i.e., graphs whose vertices can be colored into k distinct colors, so that the two endpoints of every edge have different colors. By a complete k-partite graphs we mean that any two vertices with different colors are connected by an edge, see more details in [22]. Namely, consider all graphs with N = n1 + · · · + nk vertices. Denote the corresponding k-partite graph by K(n) . d-Orthogonal Analogs of Classical Orthogonal Polynomials 25 Conjecture 6.5. For odd r, the matching polynomials Mr (K(n) ) are given by Mr (K(n) ) = x N r+1 F0   ∆(r + 1; −n1 ) . . . ∆(r + 1; −nk ) −(r + 1)r+1 . ; − (2x)r+1 More examples can be found in the cited papers. Acknowledgements The author is sincerely grateful to Boris Shapiro for sharing and discussing some polynomial systems studied here. Without this the current project would probably have never seen the light of the day. Also his advises for improvement of the text are acknowledged. The author wants to thank the Mathematics Department of Stockholm University for the hospitality in April 2015 and April 2017. Last but not least the author acknowledges extremely helpful suggestions and corrections made by the referees which helped to improve considerably the text. This research has been partially supported by the Grant No DN 02-5 of the Bulgarian Fund “Scientific research”. References [1] Abramowitz M., Stegun I.A., Pocketbook of mathematical functions, Verlag Harri Deutsch, 1984. [2] 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. [3] Alexandrov A., Open intersection numbers, Kontsevich–Penner model and cut-and-join operators, J. High Energy Phys. 2015 (2015), no. 8, 028, 25 pages, arXiv:1412.3772. [4] Alexandrov A., Cut-and-join description of generalized Brezin–Gross–Witten model, arXiv:1608.01627. [5] Aptekarev A.I., Multiple orthogonal polynomials, J. Comput. Appl. Math. 99 (1998), 423–447. [6] Aptekarev A.I., Bleher P.M., Kuijlaars A.B.J., Large n limit of Gaussian random matrices with external source. II, Comm. Math. Phys. 259 (2005), 367–389, math-ph/0408041. [7] Aptekarev A.I., Branquinho A., Van Assche W., Multiple orthogonal polynomials for classical weights, Trans. Amer. Math. Soc. 355 (2003), 3887–3914. [8] Aptekarev A.I., Kuijlaars A.B.J., Hermite–Padé approximations and ensembles of multiple orthogonal polynomials, Russian Math. Surveys 66 (2011), 1133–1199. [9] Araujo O., Estrada M., Morales D.A., Rada J., The higher-order matching polynomial of a graph, Int. J. Math. Math. Sci. (2005), 1565–1576. [10] Bakalov B., Horozov E., Yakimov M., General methods for constructing bispectral operators, Phys. Lett. A 222 (1996), 59–66, q-alg/9605011. [11] Bakalov B., Horozov E., Yakimov M., Bispectral algebras of commuting ordinary differential operators, Comm. Math. Phys. 190 (1997), 331–373, q-alg/9602011. [12] Ball K., Rivoal T., Irrationalité d’une infinité de valeurs de la fonction zêta aux entiers impairs, Invent. Math. 146 (2001), 193–207. [13] Ben Cheikh Y., On obtaining dual sequences via quasi-monomiality, Georgian Math. J. 9 (2002), 413–422. [14] Ben Cheikh Y., Ben Romdhane N., On d-symmetric classical d-orthogonal polynomials, J. Comput. Appl. Math. 236 (2011), 85–93. [15] Ben Cheikh Y., Douak K., A generalized hypergeometric d-orthogonal polynomial set, C. R. Acad. Sci. Paris Sér. I Math. 331 (2000), 349–354. [16] Ben Cheikh Y., Zaghouani A., Some discrete d-orthogonal polynomial sets, J. Comput. Appl. Math. 156 (2003), 253–263. [17] Beukers F., Padé-approximations in number theory, in Padé Approximation and its Applications, Amsterdam 1980 (Amsterdam, 1980), Lecture Notes in Math., Vol. 888, Springer, Berlin – New York, 1981, 90–99. 26 E.I. Horozov [18] Bleher P., Delvaux S., Kuijlaars A.B.J., Random matrix model with external source and a constrained vector equilibrium problem, Comm. Pure Appl. Math. 64 (2011), 116–160, arXiv:1001.1238. [19] Bochner S., Über Sturm–Liouvillesche Polynomsysteme, Math. Z. 29 (1929), 730–736. [20] Borodin A., Biorthogonal ensembles, Nuclear Phys. B 536 (1999), 704–732, math.CA/9804027. [21] Chang J.-H., The Gould–Hopper polynomials in the Novikov–Veselov equation, J. Math. Phys. 52 (2011), 092703, 15 pages, arXiv:1011.1614. [22] Chartrand G., Zhang P., Chromatic graph theory, Discrete Mathematics and its Applications (Boca Raton), CRC Press, Boca Raton, FL, 2009. [23] Craven T., Csordas G., Jensen polynomials and the Turán and Laguerre inequalities, Pacific J. Math. 136 (1989), 241–260. [24] Dattoli G., Lorenzutta S., Maino G., Torre A., Cesarano C., Generalized Hermite polynomials and superGaussian forms, J. Math. Anal. Appl. 203 (1996), 597–609. [25] de Bruin M.G., Simultaneous Padé approximation and orthogonality, in Orthogonal Polynomials and Applications (Bar-le-Duc, 1984), Lecture Notes in Math., Vol. 1171, Springer, Berlin, 1985, 74–83. [26] Delerue P., Sur le calcul symbolique à n variables et les fonctions hyperbesséliennes. II. Fonctions hyperbesséliennes, Ann. Soc. Sci. Bruxelles. Sér. I. 67 (1953), 229–274. [27] Diestel R., Graph theory, 4th ed., Graduate Texts in Mathematics, Vol. 173,, Springer, Heidelberg, 2010. [28] Dominici D., Mehler–Heine type formulas for Charlier and Meixner polynomials, Ramanujan J. 39 (2016), 271–289, arXiv:1406.6193. [29] Duistermaat J.J., Grünbaum F.A., Differential equations in the spectral parameter, Comm. Math. Phys. 103 (1986), 177–240. [30] Farrell E.J., An introduction to matching polynomials, J. Combin. Theory Ser. B 27 (1979), 75–86. [31] Forrester P.J., Wang D., Muttalib–Borodin ensembles in random matrix theory – realisations and correlation functions, Electron. J. Probab. 22 (2017), 54, 43 pages, arXiv:1502.07147. [32] Genest V.X., Vinet L., Zhedanov A., d-orthogonal polynomials and su(2), J. Math. Anal. Appl. 390 (2012), 472–487, arXiv:1110.1902. [33] Gould H.W., Hopper A.T., Operational formulas connected with two generalizations of Hermite polynomials, Duke Math. J. 29 (1962), 51–63. [34] Gutman I., Polynomials in graph theory, in Complexity in chemistry. Introduction and fundamentals, Mathematical Chemistry, Vol. 7, Editors D. Bonchev, D.H. Rouvray, Taylor & Francis, London, 2003, 133–176. [35] Heine E., Handbuch der Kugelfunctionen. Zweite umgearbeitete und vermehrte Auflage, Thesaurus Mathematicae, Vol. 1, Georg Reimer, Berlin, 1861. [36] Horozov E., Automorphisms of algebras and Bochner’s property for vector orthogonal polynomials, SIGMA 12 (2016), 050, 14 pages, arXiv:1512.03898. [37] Horozov E., Automorphisms of algebras and Bochner’s property for discrete vector orthogonal polynomials, Ann. l’Unuversité de Sofia “St. Kliment Okhridski” 104 (2017), 23–38, arXiv:1602.04343. [38] Horozov E., Vector orthogonal polynomials with Bochner’s property, Constr. Approx., to appear, arXiv:1609.06151. [39] Hosoya H., Topological index. A newly proposed quantity characterizing the topological nature of structural isomers of saturated hydrocarbons, Bull. Chem. Soc. Jpn. 44 (1971), 2332–2339. [40] Ismail M.E.H., Classical and quantum orthogonal polynomials in one variable, Encyclopedia of Mathematics and its Applications, Vol. 98, Cambridge University Press, Cambridge, 2009. [41] Koekoek R., Lesky P.A., Swarttouw R.F., Hypergeometric orthogonal polynomials and their q-analogues, Springer Monographs in Mathematics, Springer-Verlag, Berlin, 2010. [42] Konhauser J.D.E., Biorthogonal polynomials suggested by the Laguerre polynomials, Pacific J. Math. 21 (1967), 303–314. [43] Kontsevich M., Intersection theory on the moduli space of curves and the matrix Airy function, Comm. Math. Phys. 147 (1992), 1–23. [44] Kuijlaars A.B.J., Multiple orthogonal polynomials in random matrix theory, in Proceedings of the International Congress of Mathematicians, Vol. III, Hindustan Book Agency, New Delhi, 2010, 1417–1432, arXiv:1004.0846. d-Orthogonal Analogs of Classical Orthogonal Polynomials 27 [45] Kuijlaars A.B.J., Zhang L., Singular values of products of Ginibre random matrices, multiple orthogonal polynomials and hard edge scaling limits, Comm. Math. Phys. 332 (2014), 759–781, arXiv:1308.1003. [46] Lin S.-D., Chao Y.-S., Srivastava H.M., Some families of hypergeometric polynomials and associated integral representations, J. Math. Anal. Appl. 294 (2004), 399–411. [47] Maroni P., Two-dimensional orthogonal polynomials, their associated sets and the co-recursive sets, Numer. Algorithms 3 (1992), 299–311. [48] Mehler F.G., Ueber die Vertheilung der statischen Elektricität in einem von zwei Kugelkalotten begrenzten Körper, J. Reine Angew. Math. 68 (1868), 134–150. [49] Mironov A., Morozov A., Semenoff G.W., Unitary matrix integrals in the framework of the generalized Kontsevich model, Internat. J. Modern Phys. A 11 (1996), 5031–5080, hep-th/9404005. [50] Muttalib K.A., Random matrix models with additional interactions, J. Phys. A: Math. Gen. 28 (1995), L159–L164. [51] Olver F.W.J., Lozier D.W., Boisvert R.F., Clark C.W. (Editors), NIST handbook of mathematical functions, U.S. Department of Commerce, National Institute of Standards and Technology, Washington, DC, Cambridge University Press, Cambridge, 2010, available at http://dlmf.nist.gov. [52] Randić M., Morales D.A., Araujo O., Higher-order Fibonacci numbers, J. Math. Chem. 20 (1996), 79–94. [53] Sorokin V.N., Hermite–Padé approximants for Nikishin systems and the irrationality of ζ(3), Russian Math. Surveys 49 (1994), 176–177. [54] Spencer L.V., Fano U., Penetration and diffusion of X-rays. Calculation of spatial distributions by polynomial expansion, J. Research Nat. Bureau Stand. 44 (1951), 446–461. [55] Srivastava H.M., Generating functions for Jacobi and Laguerre polynomials, Proc. Amer. Math. Soc. 23 (1969), 590–595. [56] Srivastava H.M., Manocha H.L., A treatise on generating functions, Ellis Horwood Series: Mathematics and its Applications, Ellis Horwood Ltd., Chichester, Halsted Press, New York, 1984. [57] Takata T., Asymptotic formulae of Mehler-Heine-type for certain classical polyorthogonal polynomials, J. Approx. Theory 135 (2005), 160–175. [58] Toscano L., Una generalizzazione dei polinomi di Laguerre, Giorn. Mat. Battaglini (5) 4 (1956), 123–138. [59] Van Assche W., Mehler–Heine asymptotics for multiple orthogonal polynomials, Proc. Amer. Math. Soc. 145 (2017), 303–314, arXiv:1408.6140. [60] Van Assche W., Coussement E., Some classical multiple orthogonal polynomials, J. Comput. Appl. Math. 127 (2001), 317–347, math.CA/0103131. [61] Van Iseghem J., Vector orthogonal relations. Vector QD-algorithm, J. Comput. Appl. Math. 19 (1987), 141–150. [62] Vignat C., Lévêque O., Proof of a conjecture by Gazeau et al. using the Gould–Hopper polynomials, J. Math. Phys. 54 (2013), 073513, 8 pages, arXiv:1203.5418. [63] Vinet L., Zhedanov A., Automorphisms of the Heisenberg–Weyl algebra and d-orthogonal polynomials, J. Math. Phys. 50 (2009), 033511, 19 pages. [64] Zhang L., Local universality in biorthogonal Laguerre ensembles, J. Stat. Phys. 161 (2015), 688–711, arXiv:1502.03160.