$d$-Orthogonal Analogs of Classical Orthogonal Polynomials

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.


Introduction
This paper is a continuation of [22]. See also [23,24]. However we have made efforts to be read independently by including a short description of the needed results from [22]. There we constructed large families of polynomial systems that were called vector orthogonal polynomials (VOP) with Bochner's property. We recall that the VOP are polynomial systems P n (x), n = 0, 1, 2, . . . , deg(P n ) = n, which satisfy a recurrence relation of order d + 2 with d ≥ 1, cf. [39]. The VOP and the more general class of multiple orthogonal polynomials are an area of intensive research in the last 30 years due to their interesting properties and their applications, cf. [3,4,10]. See also the the cited references in [22]. The Bochner's property means that they are eigenfunctions of a linear differential (or difference) operator. Their name is motivated by the classical theorem of Bochner [11], which gives a complete classification of the case d = 1, corresponding to the orthogonal polynomials, and to order of the differential operator equal to 2. Bochner proved that all polynomials with the listed properties are the classical orthogonal polynomials of Hermite, Laguerre, Jacobi. If we treat the notion of orthogonality broader, i.e. with respect to an arbitrary linear functional on the space of all polynomials in one variable C[x], we have to add also the Bessel polynomials (see [26]. It is important to emphasize that together with the construction of the VOP we found a number of their properties analogous to the properties of the classical orthogonal polynomials -ladder operators, Rodrigues-like formulas. All these came automatically with their construction. Our main tools are based on automorphisms of algebras [6], which allow to construct the VOP from the simplest families -{ψ n (x) = x n , n = 0, 1, . . .} in the case of differential operator, and {ψ n (x) = (x) n , n = 0, 1, . . .}, where (x) n = x(x − 1) . . . (x − n + 1) in the case of difference operator. We start with a polynomial R and construct an operator G = R(x∂ x ) in the differential case and with G = R(x∆)∇, where the difference operators are ∆ = D − 1, ∇ = D −1 − 1 and Then for a polynomial q(G) we define the polynomials P n (x) by "an operational formula" of the form: P n (x) = e q(G) x n = ∞ j=0 (q(G)) j ψ n (x) j! .
We point out that the series is in fact a finite sum which justifies the definition of the polynomials P n (x). For a brief review, see the next section and for a detailed exposition see [22]. In the present paper we find more properties, which however require some extra arguments to prove in comparison with [22]. Let us list them. The most basic one, of which the rest are consequences, is the hypergeometric representations of the class of VOP corresponding to q(G) = ρG l , l ∈ N, ρ ∈ C, both in the continuous and the discrete case. For each case we obtain representations in two different hypergeometric functions. One of the representaions has the advantage that the corresponding formula is with the same parameters for all n. However the formula, which depends on n ≡ (mod l) seems to be more useful for the rest of the applications.
The hypergeometric representations show that the class is very similar to the Gould-Hopper polynomials [19], corresponding to G = ∂ x , the latter being also part of our construction. For these reasons we call these VOP generalized Gould-Hopper polynomials (GGHP), which has motivated the title of the paper. One almost obvious consequence of the hypergeometric representations is that the polynomial systems in this case can be split into l families of VOP, each originating from different operators G, exactly as in the case of Hermite polynomials. Recall that the latter can be presented as H 2n (x), n = 0, 1, 2 . . . and H 2n+1 (x), n = 0, 1, 2 . . .. The corresponding representations are H 2n (x) = c n L (−1/2) n (x 2 ) and H 2n (x) = c ′ n xL (1/2) n (x 2 ), where c n and c ′ n are constants.
It is worth mentioning that these representations also follow quite directly from our construction of the VOP in [22] but still require some transformations. The point is that in the above representations of the polynomials P n (x) each summand corresponds naturally to a summand in a generalized hypergeometric function. For the differential case when q(G) = G the hypergeometric representations are known, see [8] as well as for the Gould-Hopper polynomials. The rest of the formulas are new. For the discrete case the formulas for the families corresponding to G = ∇ and any l > 1 are found in [9]. All of the other formulas, except for the Charlier and Meixner polynomials are new.
Our next goal is to find generating functions for the GGHP. We present two formulas, both based on the second hypergeometric representation. For the first family of generating functions we use the well known method by which one obtains the generating functions for the classical orthogonal polynomials. For the second family we apply a formula by Srivastava [33]. Some of the VOP are known to have generating functions. These are exactly the ones mentioned above from [8].
Finally we find Mehler-Heine type asymptotic formulas, which are again based on the second hypergeometric representations. Such formulas in cases of multiple orthogonal polynomials have been found in [35,38]. However the discrete case apparently requires more efforts. Nevertheless, see [15], where the cases of Charlier and Meixner polynomials are successfully studied, which maybe is an indication that the discrete GGHP can also be treated.
An important remark follows. We recall that the constructions of all systems of VOP in [22] are based on the so called bispectral problem, see [17,6,7]. Namely, they are eigenfunctions for two operators -the first one is L(x) in the variable x, the other Λ(n) in the variable n: L(x)P n (x) = nP n (x), xP n (x) = Λ(n)P n (x). These type of VOP were called VOP with Bochner's property. For more details, see [22]. 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 The latter functions appear in entirely different bispectral problem -for which both variables x and n are continuous, see [7]. They, after certain rescaling, coincide with the generalized Bessel functions from [7]. It seems that this is not a mere coincidence but could be exploited further. In particular the results already known for the Bessel bispectral functions and their Darboux transformations could be used (probably rather as a model, than directly) in the study of our VOP with Bochner's property.
Other possible continuations of the present studies, apart from Darboux transformations, could be the so called linearization problem (Clebsch-Gordan coefficients) and the problem of "connection coefficients" for the GGHP. The asymptotics formulas from the paper seem promising for the study of the zeros of the VOP with Bochner's property.
The paper is organized as follows. In order to make the paper independent from [22] at least for those readers that are ready to accept it without proofs, we state the needed definitions and statements in Sect. 2. The next Sect. 3 is the central for the paper -there we derive different formulas for the GGHP in hypergeometric functions, including some well known. However our approach even for the latter is different. It is based entirely on the algebraic construction from [22] and emphasizes the common origin of all formulas, new and old. Next, in Sect. 4 from the hypergeometric representations we derive generating functions for the GGHP. In Sect. 5 we prove some results of the Mehler-Heine type asymptotics for GGHP. We finish the paper with a number of examples in Sect. 6. They include illustrative examples but also some new interpretations, e.g. for certain cases of "matching polynomials of graphs".
Acknowledgements. The author is deeply grateful to Boris Shapiro for showing and discussing some examples of systems of polynomials studied here and in particular the examples from [36]. Without this probably the project would have never be even started. The author is grateful to the Mathematics Department of Stockholm University for the hospitality in April 2015. My thanks also go to the organizers and the participants of the International Conference Groups and Rings -Theory and Applications (GRiTA2016) and in particular to V. Drensky for the invitation to the conference, where some of the results were presented. Also I appreciate the comments made by the participants of the Seminar "Dinamical systems and number theory" at Sofia University, where parts of this paper were presented.

Preliminaries
In this section we briefly recall some of the notions and the results from [22]. Let F be a field of characteristics zero. We start with the first Weyl algebra W 1 , defined over For any polynomial q(G) without a constant term 1 we defined the automorphism σ q = e q(G) of B 1 .
Then we found the images of the generators of B 1 under the automorphism σ q : with some polynomials γ j (H).
Let us introduce an auxiliary algebra R 2 defined over F with generators T, T −1 ,n subject to the relations One easily sees that [A,nB] = 1. Hence the operators A,nB define a realization of the Weyl algebra. We can introduce another algebra B 2 as follows.
First we define an anti-homomorphism b, i.e. a map b : B 1 → R 2 satisfying b(m 1 .m 2 ) = b(m 2 ).b(m 1 )), for each m 1 , m 2 ∈ B 1 by: The algebra B 2 will be the image b(B 1 ). Then b : B 1 → B 2 is an anti-isomorphism and in particular b −1 : B 2 → B 1 is well defined.
We can represent the algebra B 1 on the space of the polynomials in one variable C[x] by realizing Y as the operatorx of multiplication by x, and Z as the operator of differentiation ∂ x . Let us define the polynomial system ψ(x, n) = x n . Then the corresponding operators H, G,x act on In the same way we represent the algebra B 2 in C[x] by realizing T and T −1 as the shift operators acting on functions f (n) by T ± f (n) = f (n ± 1). Finallyn the multiplication by n. With this notation we have Using q(G) we define another polynomial system P q n (x) by Notice that the above series is in fact finite as the operator q(G) reduces the degree of any polynomial. Let us denote the operator σ q (H) by L. With these notions and notations we proved among other statements the following Theorem 2.3. The polynomials P q n (x) have the following properties: (i) They are eigenfunctions of the differential operator with eigenvalues λ(n) = n.
(ii) They satisfy a recurrence relation of the form In a similar way we realized the abstract construction of B 1 but in terms of difference operators instead of differential ones, acting again on the space of polynomials C[x]. Then we define the operators D ± acting on f ∈ C[x] by shift of the argument D ± f (x) = f (x − 1). The operatorx acts on f (x) as multiplication by x. We also need ∆ = D − 1, ∇ = D −1 − 1 and H = −x∇. Finally we put g =x − H. For this realization we use the following polynomials system ψ( is the the symbol of the falling factorial.
Lemma 2.4. The following identities hold: 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 = e ad q(G) . Also introduce the new polynomial system P q n (x) = e q(G) ψ(x, n). With this notation we proved in [22], among the other, the following properties of the polynomial system P q n (x): Theorem 2.5. The polynomials P q n (x) have the following properties: (i) They are eigenfunctions of the difference operator (ii) They satisfy a recurrence relation of the form where m = l for d = 0 and m = ld + l − 1 for d > 0.

Hypergeometric representations
This section is central for the paper. Here we derive formulas in generalized hypergeometric functions for some of the families of VOP, introduced in [22] and Sect. 2. These are the families defined via q(G) = ρG l , l ≥ 1. This property is one of those, that make our systems similar to the classical orthogonal polynomials. Notice that Laguerre polynomials L (α) n (x) correspond to G = (x∂ x + α + 1)∂ x and l = 1, ρ = 1. Hermite polynomials correspond to ∂ x and q(G) = −G 2 /2. Our treatment includes these polynomial systems, too.
Let us recall the definition of generalized hypergeometric function. Fix two 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 negative integers or zero. The function is the symbol of the rising factorial, is called generalized hypergeometric function. We have to point out that the above function does not always exist. However, when one of the parameters α i is nonpositive integer it is well defined as the series terminates and the function is a polynomial. In this paper we deal mostly with such generalized hypergeometric functions. In the rest of the 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. From now on we will skip the word "generalized", as there is no danger of confusion, although the expression "the hypergeometric function" usually means the 3.1. Continuous VOP. Some of the polynomial systems of the previous section have well known representations in terms of generalized hypergeometric functions [8], which the authors of this paper use as a definition. Below we will derive the representation from [8] on the basis of the constructions from [22]. Our proof will help us to find hypergeometric representations for the other systems studied here, that are not treated elsewhere. The discrete VOP will be treated in the same manner.
We start with the case when the automorphism σ is defined by Here some of the numbers α j can be equal. In what follows we will drop the dependence of the VOP on R.
Theorem 3.1. The polynomials P n (x) have the following hypergeometric representations: Proof. Notice that The formula for P n (x) can easily be transformed into one without differentiations using the symbol of the falling factorial.
This expression is almost in the form of a generalized hypergeometric function. For the purposes of hypergeometric functions we need to use also the notation of the rising factorial (3.4). The connection between the falling and the rising factorials is Using the symbol of rising factorial we can express the polynomials P n (x) as follows To obtain the second expression we transform the coefficients in (3.7) making use of the formula .
Let us plug this expressions for values of α = 0, α 1 , . . . , α d into the sum, defining P n (x). We obtain Write the expression n!/j! as (−1) n−j (−n) (n−j) . After changing the summation index j → s = n − j we find Of course we need to impose the condition α j = −1, −2, . . ..

Remark 3.2.
(1) This expression coincides with the corresponding one in [8]. The form of the roots of R was chosen for this purpose as well as for the hypergeometric formula of the generalized Laguerre polynomials L (α) n (x). (2) Our formula (3.8) is valid without any restrictions for the coefficients α j . However if some of the coefficients α j is a negative integer, the polynomial system becomes x n for n ≥ −α j . Hence in such a case it forms a finite VOP.
Let us consider the polynomials obtained by the automorhisms σ q , where q(G) is some polynomial without constant term. We recall that they are given by I don't know if these VOP have representations in terms of generalized hypergeometric functions or other special functions. However in the case when q(G) = G l they have. We are going to present the formula. Let us fix q(G) = G l . Finally we will drop the dependence on q(G) as there is no danger of misunderstanding.
In what follows we use the notation ∆(l; λ), which abbreviates the array of l parameters: We will also combine the latter symbol with (α) to write (∆(l; α)) for l in the expressions for the hypergeometric functions. We assume that ρ = 0. Theorem 3.3. For q(G) = G l the polynomials P n (x) have the following representation in hypergeometric functions: In particular for k = lj we come to the expression We are going to split the product into several products of falling factorials where the indexes run up to j. This can be achieved in the following way (cf. e.g. [34]). Let us start with lj−1 s=0 (n − s). We have Repeating the above argument for R(n − 1) = ρ d k=1 (n + α k ) we present it as .
For later use we write the explicit form of the coefficient at x n−lj of G lj x n : .
Then the polynomials P n (x) are given by This proves (3.10).
Remark 3.4. The formula for the VOP, corresponding to q(G) = G l , in the above theorem resembles the representation of Hermite polynomials in hypergeometric functions, see e.g [30] . Notice that we don't need to sum up to n but only up to n l as the next terms are 0.
As in the case l = 1 we are going to find a second representation. For this we need some preparations that will be useful also in the discrete case.
Let us fix i ∈ {0, 1, . . . , l − 1} and consider the polynomials with n = ml + i, m = 0, 1, . . .. Then the above parameters ∆(l; −n) can be presented as follows: Introduce the constant C ′ by and put C mi = C ′ /m!. Theorem 3.5. The polynomials P n (x), where n = ml + i have the following hypergeometric representation (3.14) Proof. We start with the expression (3.13) for P n (x) written in the form: We see that the formula for the polynomials contains several rising factorials of the type (−m − µ) (j) . They can be presented as We are going to use this presentation for values of µ = −r+i l and µ = −r+α k +i l , i = 0, . . . , l − 1. With the above notation we can transform the formula for the polynomials P n (x) into We are going to use the formula m! = (−1) m−j (−m) (m−j) j! Recall that C ′ = m!C mi and apply the last formula. It will give where we have inserted (m − j)! in the numerator and the denominator. For the numerator we will use that (m − j)! = (1) (m−j) . Also we multiply the sum by l (d+1)ml ρ ml and divide by the same expression in the sum. Put the factor x ml into the sum. After a few manipulations and changing the summation index j → s = m − j the expression takes the form which gives (3.14) Remark 3.6. Notice that b rr = 0. Hence we can cancel the equal entries in the upper and lower parameters in (3.14). This will be used later. We keep the above formula due to its symmetry.
Denote by b * the vector b(i) without the entry b ii = 0.
Corollary 3.7. The polynomials P n (x) have the hypergeometric representation While the above statement does not need a proof it is worth to point out that the representation shows that the VOP originating from G l are modifications of the VOP constructed from another algebra with R 1 (H) = d k=1 l−1 r=0 (H + c irk + 1) r =i (H + b ir + 1). Let us introduce the polynomials Q ml+i (x l ) = P ml+i (x)x −i , m = 0, 1, . . .. We see that Q ml+i (y), y = x l form VOP for each fixed i. This again shows the resemblance to the Hermite polynomials, which can be represented as Laguerre polynomials L (α) n (x) for certain values of the parameter α.
In [22] (see Example 7.3) we showed that such a splitting of the family exists for the simple case of q(G) = G 2 , where G = (x∂ x + α + 1)∂ x without using hypergeometric functions. Also in Remark 7.4. of the same paper we pointed that the fact is true for any G l with arbitrary G. In Corollary 3.7 we have given a proof of this result based on hypergeometric functions.
Remark 3.8. Notice that the statements of Theorem 3.3 and Theorem 3.5 differ considerably in form, when l > 1. While Theorem 3.3 presents a formula valid for each n, the formulas in Theorem 3.5 depend on the residue of n mod l. The same is true for the discrete case.
3.2. Discrete VOP. 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 VOP. These representations are new except for the cases of Charlier and Meixner polynomials.
We again factor the polynomial R(H) into R(H) = ρ m j=1 (H + α j + 1). Theorem 3.9. The discrete VOP P n (x) have the following hypergeometric representations: Proof. First we transform the formula for P B n (x) to get rid of the difference operators. We are going to exploit again the formula Using that G(x) n = nR(n − 1)(x) n−1 , we obtain Hence making use of the falling factorial symbol (a) j we obtain With the help of the formula The second hypergeometric expression can be obtained as in the continuous case. Using again (3.9) we find After standard manipulations that we used for the continuous VOP we come to the second formula (3.19).
We can obtain hypergeometric representations for the class of discrete VOP, constructed via q(G), where G = R(H)∆. Again we will treat the case when q(G) = G l . Theorem 3.10. In the case when q(G) = G l the polynomials P n (x) have the following representation in hypergeometric functions: 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: A formula, similar to the one in the continuous case holds: Only the term (x) n−lj needs some attention. We have (x) n = (x) n−lj (x − n + 1) (lj) .
The factor (x − n + 1) (lj) can be presented as indicated in section 3.1: . This gives Now we can express the polynomials as Using the last expression and (3.12) for we obtain the desired formula Again we can find a second formula for the VOP, corresponding to q(G) = G l . We use the notations (3.1) from subsection 3.1. With this notation we have Theorem 3.11. The polynomials P n (x), where n = ml + i have the following hypergeometric representation Proof. We start with the formula For the expression lj−1 s=0 (n − s)R(n − 1 − s) j! we use (3.12) (see also the proof of Theorem 3.5) to transform it to (3.22) ρ where v = m − j.
Next we transform the factor (x) n−lj as follows. First Then we present the second factor in the right-hand side of the last formula as .
In this way we obtain . .
At the end multiplying (3.23) and (3.22) and summing up we obtain which is (3.21).

Generating functions
In this section we will find generating functions for all VOP for which we obtained hypergeometric representations.
Generating functions for the continous VOP, corresponding to q(G) = G, were found in [8].

Continuous GGH Polynomials.
In what follows we use different normalization for the VOP. Namely we introduce the polynomials They differ from P n (x) by a multiplicative constant. We also put ((−l) d+1 ρ) l = −1, which can be achieved by rescaling x together with multiplying by a suitable constant. Let us fix i. Our first result is Proof. Consider the series ; −x l .
Let us present the sum in the right-hand side of the last equation in the form of a double series.
Changing the order of summation of the double series we obtain It is easy to see that Hence by introducing a new index of summation s = m − j we obtain that the double series becomes This gives for the double sum which implies (4.25) for the function Φ i (x, t).
Notice that the numbers b ii = 0. This shows that one of the numbers in (b(i) + 1) is equal to 1. Let us denote the array (b(i) + 1) without the element equal to 1 by (b * (i) + 1).
Proof. follows immediately from the obvious equality where µ = 0 is otherwise an arbitrary number.

Remark 4.3.
We notice that the hypergeometric functions without upper parameters appear in a quite different bispectral problem. In [7] we defined the functions Ψ β (x, z), which are solutions of an equation of the form , z), where θ = x∂ x and β j ∈ N. We called them generalized Bessel functions. They are expressed in terms of the hypergeometric functions without upper parameters (or Meijer's G-functiones). 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 the case of l = 1 (for which the same formula was found by different methods in [8], see also below Corollary 4.5) the connection deserves attention.
From the above formulas we can write a generating function for the entire family Q n . Let us define the function Proof. is obvious.
It deserves to write separately the formula for the case l = 1. ; xt e t , We are going to obtain a second formula based on a theorem from [33]. 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 , . . . , β p be complex numbers such that the hypergeometric function be well defined. Then the following formula holds We again use the VOP Q n (x) from (4.24) as well as the convention l d+1 ρ = 1 Theorem 4.7. The function G i given by generates the polynomials Q lm+i (x), m= 0, 1, . . . in the following way Proof. Let us multiply the polynomials Q lm+i by t n and sum. Notice that in our case a = 0, which implies m m = 1.
Thus we obtain Next we use that to represent the polynomials Q n as We apply the above cited formula (4.26) from [33] with a = 0 to obtain From this theorem we can write a generating function for all polynomials Q n (x) Corollary 4.8. A generating function for Q n (x) is given by Proof. is obvious. Just sum up the generating functions Φ i (x, t) for i = 0, . . . , l − 1 and replace n by ml + i.
Notice that the coefficients c and b 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 VOP defined in terms of q(G) = ρG have the following generating function

4.2.
Discrete generalized GH 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.21) which differ only by a multiplicative constant from P n (x).
Theorem 4.10. The function Φ i given by is a generating function for the polynomials Q lm+i (y) in the sense that 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 Notice that This gives ; l d+1 ρt l .
As in the continuous case we are going to derive a second formula.
Theorem 4.11. The function G i given by generates the polynomials Q lm+i (y) in the sense that Proof. We are going to use again Proposition 4.6 . It is obvious that we need to put a = 0. We have Then (4.26) gives As a trivial corollary again we can write a generating function for all polynomials Q n (x). Corollary 4.12. A generating function for Q n (x) is given by 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

Mehler-Heine type formulas
We are going to consider only the continuous VOP. In what follows we will assume without loss of generality that ρl d+1 = 1. This will make the formulas simpler. There are many ways to write down the Mehler-Heine type formulas depending on the normalization of the polynomials. We choose some of the simplest. Let us start with the case of continuous VOP corresponding to q(G) = G. Also we use the polynomials Q n (x) given by (4.24).
Proof. Notice that in the case l = 1 (4.26) can be written as ; −x .
Remark 5.2. As in the formulas for the generating functions for GGHP the hypergeometric function ; x d appears once again. It is interesting to understand if this connection is deeper or just a coincidence.
Let us consider the polynomials obtained by the automorphisms σ q , where q(G) is some polynomial. We recall that they are given by (5.28) ψ 1 (x, n) := P n (x) := e q(G) x n = x n + ∞ k=1 q k (G)x n k! .
We will be concerned with the case when q(G) = G l for which we found a hypergeometric representation in section 3. We will modify slightly this formula to make it convenient for our purposes. We will introduce some notation in order to simplify the formulas. Now we are ready to give the main result. Let us again present n as n = ml + i. We use the polynomials Q(x) introduced in the previous section: Theorem 5.3. For the VOP obtained via the automorphisms σ q with q(G) = G l the following Mehler-Heine type formula holds Proof. We are going to use the the hypergeometric representation (3.14).
Consider the expression Q n (x/m 1/l ). Explicitly it is Rewrite the last formula in the form We again use the formula from [26], p. 5, (1.4.3). Hence we obtain It is worth noticing that the asymptotics depends on the residue of n(mod l). This points to the fact that there is no general asymptotic formula but only in the above form. i.e. only for subsequence with the same residue i. This is exactly as in the case of Hermite polynomials, where the even-indexed and the odd-indexed ones have different asymptotics, see, e.g. [1].
Notice that the function is also a generalized Bessel function in the sense of [7] as we explained in Corollary 4.2.
According to our scheme they are eigenfunctions of the differential operator L = lτ ∂ l + x∂ and satisfy the recurrence relations xg l n (x, τ ) = g l n+1 (x, τ ) − τ ln(n − 1) . . . (n − l + 2)g l n−l+1 (x, τ ). Using our second form of hypergeometric representation (3.14) we can write them also as where c * is the vector c(i) without the entry c ii . Notice that for l = 2 these are the classical Hermite polynomials (up to rescaling). During the years after their discovery the Gould-Hopper polynomials turned out to be quite useful in studies of quantum mechanics, integrable systems (Novikov-Vesselov equation), combinatorics, etc. see, e.g. [12,13,40].
The cases with G = R(H)∂ with arbitrary R can be considered generalizations of Gould-Hopper polynomials. We only have to take q(G) = G l , where G = R(H)∂ with R a polynomial of any degree. Then the corresponding expression for these polynomials take the form from (3.10) We can write explicit discrete analogs of the (generalized) Gould-Hopper polynomials using the map M * . We have The most direct analog, which could be called discrete Gould-Hopper polynomials, is given by G = τ ∆ l (d = 0). They have the explicit representation Having in mind the applications of the Gould-Hopper polynomials it would be interesting to find if the generalized versions have similar or other applications. The two families are bi-orthogonal with respect to the measure of Laguerre polynomials: with h m = 0. Earlier the polynomials Z α n (x; d) were found by Toscano [37]. Their hypergeometric representation [28] Z α n (x; d) = α + dn dn ; (x/d) d shows that they can be constructed via the methods of the present paper. Let us consider G = R(H)∂, where the polynomial R(H) = d s=1 (H + α+s d ). Then we see that Z α n (x; d) = P n ((x/d) d )). In the case d = 2 the polynomials were found much earlier by L. V. Spencer and U. Fano [32] in their studies of X-rays diffusion. Remark 6.3. Konhauser has proved in [27] that the polynomials of the family Y α n (x; d) are eigenfunctions of a differential operator of order d + 1 and that they satisfy d + 2 recurrence relation of the form The polynomials Z α n (x; d) satisfy a relation of the form which follows from the definition of Z n (x; d), i.e Z n (x; d) = P R n ((x/d) d )) and the properties of P R n ((x/d) d )). It looks interesting to define both continuous and discrete analogs for the Konhauser polynomials from the VOP, studied here. Example 6.4. Matching polynomials of graphs.
The polynomial systems considered in this example are relevant for a branch of research, known under the name of chemical graph theory. They are taken from [5].
In this example we use some notions from graph theory. One can consult [14,18,5] and the references therein. Let K be a connected graph with n vertices. Following [31] we define the higher Hosoya numbers p r (K, j) as the number of ways that one can select j nonincident paths of length r in K. A higher-order matching polynomial M r (K n ) of the graph K n via Hosoya numbers is defined in [18] , see also [5] by In the case when K = K n is a complete graph (i.e. each pair of vertices is connected with exactly one edge) the corresponding polynomials have been computed explicitly in [18,5]. It was found by combinatorial arguments that the polynomials are M r (K n ) = n j=0 (−1) j n!x n−(r+1)j (n − (r + 1)j)!j!2 j .
This representation was found in [5]. We see that the matching polynomials of complete graphs are eigenfunctions of a differential operator and they satisfy an (r + 2)-term recurrence relation. This facts could be useful in the studies of the polynomials. Of course, the case r = 1 deserves special attention; it corresponds to Hermite polynomials as observed long ago, see e.g. [5,20]. The corresponding coefficients p 1 (K, j) are the original Hosoya numbers.
Let us consider another example of matching polynomials, this time of complete bipartite graphs K n,m , n ≥ m with n+m vertices. We recall that this means a complete graph whose vertices are a union of two nonintersecting sets V n and V m with n and m elements and any edge connects a vertix in V n to a vertex in V m . We consider the case when r is odd. Then in [18,5] one can find the formula Let us put m = n − M, M ≥ 0. In this notation we see that the polynomials M r (K n,m ) become M r (K n,m ) = x 2n−M r+1 F 0 ∆(r + 1; −n) ∆(r + 1; −n + M ) − ; −(r + 1) r+1 (2x) r+1 .
They coincide with x n−M P n (x), where the polynomials P n (x) are constructed via q(G) = −G r+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 are new.
Notice that for r = 1 (when all the paths are edges) these are the polynomials L