Classical and Semi-Classical Orthogonal Polynomials Defined by Riordan Arrays,

and Their Moment Sequences

Paul Barry and Arnauld Mesinga Mwafise School of Science Waterford Institute of Technology Ireland

                         pbarry@wit.ie

Abstract

We study the orthogonal polynomials of classical and semi-classical types that can be defined by ordinary and exponential Riordan arrays. We identify their moment sequences, giving their integral representations and Hankel transforms. For a special class of classical orthogonal polynomials defined by Riordan arrays, we identify a complementary family of orthogonal polynomials defined by reversion of moment sequences. Special product sequences arise and their generating functions are calculated.

1 Introduction

Riordan arrays [4, 28, 29, 30, 34] are simple to define (see below), providing a bridge between elements of algebra, group theory and linear algebra. This combination can shed light on other areas of mathematics. In this note, we show how Riordan arrays can yield fresh perspectives on the area of orthogonal polynomials. It is straight-forward to classify those Riordan arrays that define orthogonal polynomials - essentially, they are the ordinary Riordan arrays whose production matrices are tri-diagonal. Nevertheless, it is interesting to note that only in a limited number of cases are the associated orthogonal polynomials of classical type. This note explores this fact. We classify those ordinary Riordan arrays that define classical polynomials, and we study some integral representations of the moment sequences associated with semi-classical orthogonal polynomials defined by ordinary Riordan arrays.

We also begin the investigation of classical orthogonal polynomials defined by exponential Riordan arrays, and we look at families of such polynomials that are related to the reverse Bessel polynomials.

The paper is laid out in the following sections.

1. This Introduction

2. Preliminaries on ordinary Riordan arrays and orthogonal polynomials

3. Classical and semi-classical orthogonal polynomials

4. The classical orthogonal polynomials defined by ordinary Riordan arrays

5. The non-classical case

6. Complementary orthogonal polynomials defined by Riordan arrays

7. A note on the INVERT transform

8. A special product sequence

9. Exponential Riordan arrays and classical orthogonal polynomials

10. Bessel and related polynomials

11. Conclusion

12. Declaration

13. Appendix — The Stieltjes transform of a measure.

We recall the following well-known results (the first is known as “Favard’s Theorem”), which we essentially reproduce from [18], to specify the links between orthogonal polynomials, three term recurrences, and the recurrence coefficients and the generating function of the moment sequence of the orthogonal polynomials. Theorem 1. [18] (Cf. [37, Théorème 9, p. I-4] or [38, Theorem 50.1]). Let (p_{n}(x))_{n≥0} be a sequence of monic polynomials, the polynomial p_{n}(x) having degree n = 0, 1, … Then the sequence (p_{n}(x)) is (formally) orthogonal if and only if there exist sequences (α_{n})_{n≥0} and (β_{n})_{n≥1} with β_{n} ≠ 0 for all n ≥ 1, such that the three-term recurrence

p_{n+1} = (x − α_{n})p_{n}(x) − β_{n}p_{n−1}(x), for n ≥ 1,

holds, with initial conditions p_{0}(x) = 1 and p_{1}(x) = x − α_{0}.

Theorem 2. [18] (Cf. [37, Prop. 1 (7), p. V-5] or [38, Theorem 51.1]). Let (p_{n}(x))_{n≥0} be a sequence of monic polynomials, which is orthogonal with respect to some functional L. Let

p_{n+1} = (x − α_{n})p_{n}(x) − β_{n}p_{n−1}(x), for n ≥ 1,

be the corresponding three-term recurrence which is guaranteed by Favard’s theorem. Then the generating function

g(x) = ∑_{k=0}^{∞} μ_{k}x^{k}

for the moments μ_{k} = L(x^{k}) satisfies

g(x) = \frac{μ_{0}}{1 − α_{0}x − \frac{β_{1}x^{2}}{1 − α_{1}x − \frac{β_{2}x^{2}}{1 − α_{2}x − \frac{β_{3}x^{2}}{1 − α_{3}x − ⋯}}}}.

The Hankel transform of a given sequence A = {a_{0}, a_{1}, a_{2}, ...} is the sequence of Hankel determinants {h_{0}, h_{1}, h_{2}, … } where h_{n} = |a_{i+j}|_{i,j=0}^{n}, i.e

A = \{a_{n}\}_{n∈ℕ_{0}} → h = \{h_{n}\}_{n∈ℕ_{0}} : h_{n} = | \begin{matrix}a_{0} & a_{1} & ⋯ & a_{n} \\ a_{1} & a_{2} & & a_{n+1} \\ ⋮ & & ... \\ a_{n} & a_{n+1} & & a_{2n} \end{matrix} | .  (1)

The Hankel transform of a sequence a_{n} and its binomial transform are equal.

In the case that a_{n} has a generating function g(x) expressible in the form

g(x) = \frac{a_{0}}{1 − α_{0}x − \frac{β_{1}x^{2}}{1 − α_{1}x − \frac{β_{2}x^{2}}{1 − α_{2}x − \frac{β_{3}x^{2}}{1 − α_{3}x − ⋯}}}}

then we have [18]

h_{n} = a_{0}^{n+1}β_{1}^{n}β_{2}^{n−1} ⋯ β_{n−1}^{2}β_{n} = a_{0}^{n+1} ∏_{k=1}^{n} β_{k}^{n+1−k}.  (2)

Note that this is independent of α_{n}.

2 Preliminaries on Riordan arrays and orthogonal polynomials

An ordinary Riordan array M is an invertible lower-triangular matrix defined by two power series

g(x) = 1 + g_{1}x + g_{2}x^{2} + ⋯

and

f(x) = x + f_{2}x^{2} + f_{3}x^{3} + ⋯,

where the (n, k)-th element of the corresponding matrix is given by

m_{n,k} = [x^{n}]g(x)f(x)^{k},

where [x^{n}] is the operator that extracts the coefficient of x^{n} in the power series upon which it operates [24]. The variable “x” here is a dummy or synthetic variable, in the sense that we have m_{n,k} = [t^{n}]g(t)f(t)^{k} using another designation for this variable. Note that we have chosen g_{0} = 1 and f_{1} = 1 here, to simplify the exposition. The power series above are in ordinary form g(x) = ∑_{n=0}^{∞} g_{n}x^{n} and f(x) = ∑_{n=0}^{∞} f_{n}x^{n} (with f_{0} = 0) and as such the associated arrays are called ordinary Riordan arrays.

Along with the two power series g(x) and f(x), we can define two associated power series

A(x) = \frac{x}{f̄(x)},

and

Z(x) = \frac{1}{f̄(x)} (1 − \frac{1}{g(f̄(x))}) .

Note that the notation f̄(x) denotes the compositional inverse of the power series f. Thus we have

f̄(f(x)) = x, and f(f̄(x)) = x.

We shall also use the notation f̄(x) = Rev(f)(x).

The (infinite) matrix whose bivariate generating function is given by

Z(x) + \frac{A(x)y}{1 − xy}

is called the production matrix of M [13, 14, 27]. It is equal to the matrix

M^{−1}M̄

where M̄ is the matrix M with its first row removed. It is an infinite lower Hessenberg matrix.

If the Riordan array M has a production matrix P defined by A(x) and Z(x), then we can show that

M^{−1} = (1 − x\frac{Z(x)}{A(x)}, \frac{x}{A(x)}) .

If Z(x) = γ + δx, and A(x) = 1 + αx + βx^{2} (where we assume that δ ≠ 0 and β ≠ 0), then P will be tri-diagonal. The matrix P begins

( \begin{matrix}γ & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ δ & α & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & β & α & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & β & α & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & β & α & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & β & α & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & β & α & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & β & α \end{matrix} ) .

In this case, we have

M^{−1} = (1 − x\frac{γ + δx}{1 + αx + βx^{2}} , \frac{x}{1 + αx + βx^{2}} ) = (\frac{1 + (α − γ)x + (β − δ)x^{2}}{1 + αx + βx^{2}} , \frac{x}{1 + αx + βx^{2}} ) ,

and R = M^{−1} will be the coefficient array of a family of orthogonal polynomials P_{n}(x) [6] where

P_{n}(x) = ∑_{k=0}^{n} p_{n,k}x^{k},

where the general (n, k)-th term of R is p_{n,k}. In this case, we call M the moment matrix of the family of orthogonal polynomials P_{n}(x) [5].

An advantage of using Riordan arrays wherever possible is that the set of Riordan arrays is a group for matrix multiplication, which in terms of the power series definition of a Riordan array is translated as

(g(x), f(x)) · (u(x), v(x)) = (g(x)u(f(x)), v(f(x)).

In conformity with this definition of multiplication, we have the rule

(g(x), f(x)) · h(x) = g(x)h(f(x)),

which mirrors the operation of multiplying the vector of elements whose generating function is h(x) by the Riordan array M. The vector resulting, when regarded as a sequence of elements, will then have generating function g(x)h(f(x)). The inverse of the array M defined by (g(x), f(x)) is given by

(g(x), f(x))^{−1} = (\frac{1}{g(f̄(x))}, f̄(x)) .

The identity element is (1, x), which is represented by the usual identity matrix.

3 Classical and semi-classical orthogonal polynomials

The classical orthogonal polynomials of mathematical science are the Jacobi, Laguerre and Hermite polynomials, defined by the weights w_{J}(x) = (1 − x)^{α}(1 + x)^{β} on [−1, 1], w_{L}(x) = x^{α}e^{−x} on [0, ∞), and w_{H}(x) = e^{−x^{2}} on (−∞, ∞), respectively. In particular, these orthogonal polynomials [11, 15, 36] are associated with measures that are absolutely continuous. We have

\frac{w′_{J}(x)}{w_{J}(x)} = \frac{x(α + β) + α − β}{x^{2} − 1} ,

\frac{w′_{L}(x)}{w_{L}(x)} = \frac{α − x}{x},

and

\frac{w′_{H}(x)}{w_{H}(x)} = −2x.

In general, we shall define a family of orthogonal polynomials P_{n}(x) to be classical [35] if the associated measure is absolutely continuous with weight function w(x) satisfying

\frac{w′(x)}{w(x)} = \frac{U(x)}{V (x)} = \frac{u_{0} + u_{1}x}{v_{0} + v_{1}x + v_{2}x^{2}} .

The polynomials y = P_{n}(x) will then satisfy the differential equation

V (x)y′′ + (U(x) + V′(x))y′ − n(u_{1} + (n + 1)v_{2})y = 0.

If deg(V) > 2 and/or deg(U) > 1 then we say that the family of polynomials is semi-classical.

Note that all orthogonal polynomials that we shall consider later will be monic (the coefficient of x^{n} in P_{n}(x) is 1).

4 The classical orthogonal polynomials defined by ordinary Riordan arrays

We have seen that for an ordinary Riordan array to define a family of orthogonal polynomials, it must be of the form

R = (\frac{1 + cx + dx^{2}}{1 + ax + bx^{2}} , \frac{x}{1 + ax + bx^{2}} ) .

This matrix will then be the coefficient array of the family of polynomials. The procedure to find the weight function associated with this family is as follows. First, we form the moment matrix M = R^{−1} given by

\begin{matrix}(−\frac{(b − d)\sqrt{1 − 2ax + x^{2}(a^{2} − 4b)} + x(a(b + d) − 2bc) − b − d}{2(x^{2}(a^{2}d − ac(b + d) + b^{2} + b(c^{2} − 2d) + d^{2}) + x(c(b + d) − 2ad) + d)}, \\ \frac{1 − ax − \sqrt{1 − 2ax + x^{2}(a^{2} − 4b)}}{2bx}). \end{matrix}

The first element of this array is the generating function μ(x) of the moments of the family of orthogonal polynomials P_{n}(x). These moments begin

1, a − c, a^{2} − 2ac + b + c^{2} − d, a^{3} − 3a^{2}c + a(3b + 3c^{2} − 3d) − c(2b + c^{2} − 2d), …

and their generating function is given by

μ(x) = \frac{1}{1 − (a − c)x − \frac{(b − d)x^{2}}{1 − ax − \frac{bx^{2}}{1 − ax − \frac{bx^{2}}{1 − ⋯}}}}.

From this we can see that the Hankel transform [18, 21] of this sequence of moments is given by

h_{n} = (b − d)^{n}b^{\binom{n}{2}}.

Our next step is to use the Stieltjes-Perron theorem [5, 17] (see Appendix) to derive the associated measure. We find that the measure sought is given by w(x)dx where

w(x) = \frac{1}{2π} \frac{(b − d)\sqrt{4b − (x − a)^{2}}}{dx^{2} + x(c(b + d) − 2ad) + a^{2}d − ac(b + d) + b^{2} + b(c^{2} − 2d) + d^{2}} .

Finally, we form the ratio \frac{w′(x)}{w(x)} to obtain the expression

−\frac{dx^{3} − 3adx^{2} + x(3a^{2}d − b^{2} − b(c^{2} + 6d) − d^{2}) − a^{3}d + a(b^{2} + b(c^{2} + 6d) + d^{2}) − 4bc(b + d)}{((x − a)^{2} − 4b)(dx^{2} + x(c(b + d) − 2ad) + a^{2}d − ac(b + d) + b^{2} + b(c^{2} − 2d) + d^{2})}.

The form of this ratio now tells us that the orthogonal polynomials defined by ordinary Riordan arrays are at least semi-classical. Inspection of the above ratio allows us to announce the following results.

Proposition 3. The ordinary Riordan array

(\frac{1 + cx + dx^{2}}{1 + ax + bx^{2}} , \frac{x}{1 + ax + bx^{2}} )

defines a family of classical orthogonal polynomials in the case that either c = d = 0 or c = 0, d = −b.

Corollary 4. When c = d = 0, we have

w(x) = \frac{1}{2π} \frac{\sqrt{4b − (x − a)^{2}}}{b}

on the interval

[a − 2\sqrt{b}, a + 2\sqrt{b}],

with

\frac{w′(x)}{w(x)} = \frac{x − a}{(x − a)^{2} − 4b}.

The moments μ_{n} have integral representation

μ_{n} = \frac{1}{2π} ∫_{a−2\sqrt{b}}^{a+2\sqrt{b}} x^{n} \frac{\sqrt{4b − (x − a)^{2}}}{b} dx.

The moments have generating function

μ(x) = \frac{1 − ax − \sqrt{(1 − ax)^{2} − 4bx^{2}}}{2bx^{2}}

given by

μ(x) = \frac{1}{1 − ax − \frac{bx^{2}}{1 − ax − \frac{bx^{2}}{1 − ax − \frac{bx^{2}}{1 − ⋯}}}}.

By an application of Lagrange inversion [23], we obtain

μ_{n} = \frac{1}{n + 1}[x^{n}](1 + ax + bx^{2})^{n+1}
= \frac{1}{n + 1} ∑_{k=0}^{n} \binom{n + 1}{j}\binom{j}{n − j}a^{2j−n}b^{n−j}

= \frac{1}{n + 1} ∑_{k=0}^{n} \binom{n + 1}{n − k}\binom{n − k}{k}a^{n−2k}b^{k}.

The moments have Hankel transform

h_{n} = b^{\binom{n+1}{2}}.

The polynomials P_{n}(x) satisfy the three-term recurrence

P_{n}(x) = (x − a)P_{n−1}(x) − bP_{n−2}(x), n > 1,

with P_{0}(x) = 1, P_{1}(x) = x − a.

If y = P_{n}(x) then y satisfies the differential equation

(x^{2} − 2ax + a^{2} − 4b)y′′ + 3(x − a)y′ − n(n + 2)y = 0.

The form of the differential equation satisfied by these polynomials follows from the form of the weight function [35]. For a proof of the other statements, the reader is referred to [6].

Corollary 5. When c = 0 and d = −b, we have

w(x) = \frac{1}{π} \frac{1}{\sqrt{4b − (x − a)^{2}}}

on the interval

[a − 2\sqrt{b}, a + 2\sqrt{b}],

with

\frac{w′(x)}{w(x)} = \frac{a − x}{(x − a)^{2} − 4b}.

The moments μ_{n} have integral representation

μ_{n} = \frac{1}{π} ∫_{a−2\sqrt{b}}^{a+2\sqrt{b}} x^{n} \frac{1}{\sqrt{4b − (x − a)^{2}}} dx.

The moments have generating function

μ(x) = \frac{1}{\sqrt{(1 − ax)^{2} − 4bx^{2}}}

given by

μ(x) = \frac{1}{1 − ax − \frac{2bx^{2}}{1 − ax − \frac{bx^{2}}{1 − ax − \frac{bx^{2}}{1 − ⋯}}}}.

We have the closed form expression for the moments

μ_{n} = ∑_{i=0}^{n} \binom{n − i}{i}\binom{n − i − 1/2}{n − i}(−1)^{i}(a^{2} − 4b)^{i}(2a)^{n−2i}

= \frac{1}{4^{n}} ∑_{k=0}^{n} \binom{2n − 2k}{n − k}\binom{2k}{k}(a + 2\sqrt{b})^{k}(a − 2\sqrt{b})^{n−k}.

The moments have Hankel transform

h_{n} = 2^{n}b^{\binom{n+1}{2}}.

The polynomials P_{n}(x) satisfy the three-term recurrence

P_{n}(x) = (x − a)P_{n−1}(x) − bP_{n−2}(x), n > 2,

with P_{0}(x) = 1, P_{1}(x) = x − a, and P_{2}(x) = (x − a)^{2} − b(b + 1).

If y = P_{n}(x) then y satisfies the differential equation

(x^{2} − 2ax + a^{2} − 4b)y′′ + (x − a)y′ − n^{2}y = 0.

We note that in the case c = 0 and d = −b, the generating function

\frac{1}{\sqrt{(1 − ax)^{2} − 4bx^{2}}} = \frac{1}{\sqrt{1 − 2ax + x^{2}(a^{2} − 4b)}}

can be compared with the generating function

\frac{1}{1 − 2xt + t^{2}} = ∑_{n=0}^{∞} P_{n}(x)t^{n}

of the Legendre polynomials P_{n}(x). Then we get [26]

μ_{n} = (a^{2} − 4b)^{n/2}P_{n} (\frac{a}{\sqrt{a^{2} − 4b}}) .

For the next result, we note that

C_{n} = \frac{1}{n + 1}\binom{2n}{n}

is the n-th Catalan number \underline{A000108}. The generating function of the Catalan numbers is given by

c(x) = \frac{1 − \sqrt{1 − 4x}}{2x}.

Proposition 6. The ordinary Riordan array (\frac{1}{1+ax}, \frac{x}{(1+ax)^{2}} ) (a ≠ 0) is the coefficient array of a family of classical orthogonal polynomials. We have

w(x) = \frac{1}{2π} \frac{\sqrt{x(4a − x)}}{2ax}

on the interval

[0, 4a],

with

\frac{w′(x)}{w(x)} = \frac{2a}{x(x − 4a)}.

The moments have integral representation

μ_{n} = \frac{1}{2π} ∫_{0}^{4a} x^{n} \frac{\sqrt{x(4a − x)}}{2ax} dx = a^{n}C_{n}.

The moments have generating function

μ(x) = \frac{1 − \sqrt{1 − 4ax}}{2ax},

with

μ(x) = \frac{1}{1 − \frac{ax}{1 − \frac{ax}{1 − \frac{ax}{1 − ⋯}}}},

or equivalently,

μ(x) = \frac{1}{1 − ax − \frac{a^{2}x^{2}}{1 − 2ax − \frac{a^{2}x^{2}}{1 − 2ax − ⋯}}}.

The moments μ_{n} have Hankel transform

h_{n} = a^{n(n+1)}.

The polynomials P_{n}(x) satisfy the three-term recurrence

P_{n}(x) = (x − 2a)P_{n−1}(x) − a^{2}P_{n−2} n > 1,

with P_{0}(x) = 1, P_{1}(x) = x − a. If y = P_{n}(x) then y satisfies the differential equation

x(x − 4a)y′′ + 2(x − a)y′ − n(n + 1)y = 0.

Example 7. The Riordan array (\frac{1}{1+x^{2}} , \frac{x}{1+x^{2}} ) is the coefficient array of the scaled Chebyshev polynomials of the second kind P_{n}(x) = U_{n}(x/2), with their moments being the aerated Catalan numbers

1, 0, 1, 0, 2, 0, 5, 0, 14, 0, …

given by

μ_{n} = \frac{1}{2π} ∫_{−2}^{2} \sqrt{4 − x^{2}} dx.

This is closely related to Wigner’s semicircle distribution [39, 40]. Note that if we define

P_{n−1}^{(1)} = \frac{1}{2π} ∫_{−2}^{2} \frac{P_{n}(z) − P_{n}(x)}{z − x}\sqrt{4 − x^{2}} dx

then we find that

P_{n}^{(1)}(x) = P_{n}(x) = U_{n}(x/2).

Here,

U_{n}(x) = ∑_{k=0}^{⌊ \frac{n}{2} ⌋} \binom{n − k}{k}(−1)^{k}(2x)^{n−2k}

are the Chebyshev polynomials of the second kind [22].

Example 8. The Riordan array (\frac{1−x^{2}}{1+x^{2}} , \frac{x}{1+x^{2}} ) is closely related to the Chebyshev polynomials of the first kind. This array is the coefficient array of a family of orthogonal polynomials P_{n}(x) whose moments are the aerated central binomial numbers (\underline{A000984})

1, 0, 2, 0, 6, 0, 20, 0, 70, 0, 252, …

given by

\frac{1}{π} ∫_{−2}^{2} x^{n} \frac{1}{\sqrt{4 − x^{2}}} dx.

In this case we have

P_{n}(x) = U_{n}(x/2) − U_{n−2}(x/2).

We also find that P_{n}^{(1)}(x), where

P_{n−1}^{(1)} = \frac{1}{π} ∫_{−2}^{2} \frac{P_{n}(z) − P_{n}(x)}{z − x} \frac{1}{\sqrt{4 − x^{2}}} dx

also satisfies

P_{n}^{(1)}(x) = P_{n}(x) = U_{n}(x/2).

5 The semi-classical case

By the results of the last section, if an ordinary Riordan array is the coefficient array of a family of polynomials that is not of classical type, then the ratio \frac{w̃′}{w̃} is of semi-classical type, for the absolutely continuous part of the measure.

We begin this section by looking at a family of orthogonal polynomials said to be of “restricted Chebyshev Boubaker type” [3]. These are ordinary Riordan arrays of the form

(\frac{1 + rx^{2}}{1 + x^{2}} , \frac{x}{1 + x^{2}} ) .

We exclude the case r = −1, which is of classical type. The polynomials P_{n}(x; r) = P_{n}(x) defined by these arrays are given by

P_{n}(x; r) = ∑_{k=0}^{⌊ \frac{n}{2} ⌋} \binom{n − k}{k}\frac{n − (r + 1)k}{n − k} (−1)^{k}x^{n−2k}.

The case r = 3 corresponds to the family of Boubaker polynomials [1, 2, 8, 9, 10, 12, 19, 20, 25, 31, 41]. The moments μ_{n}(r) of this family of orthogonal polynomials have generating function

μ(x; r) = \frac{\sqrt{1 − 4x^{2}}(r − 1) + r + 1}{2(r + x^{2}(r − 1)^{2})},

which can be expressed as the continued fraction [38]

μ(x; r) = \frac{1}{1 − \frac{(1 − r)x^{2}}{1 − \frac{x^{2}}{1 − \frac{x^{2}}{1 − ⋯}}}}.

We note that the Hankel transform of μ_{n}(r), which by the above is an aerated sequence, and that of its un-aerated version, is given by

h_{n}(r) = (1 − r)^{n}.

Further, the un-aerated moments

1, 1 − r, r^{2} − 3r + 2, −r^{3} + 5r^{2} − 9r + 5, …

are themselves moments for the family of orthogonal polynomials that have coefficient matrix given by

(\frac{(1 + x)(1 + rx)}{(1 + x)^{2}}, \frac{x}{(1 + x)^{2}} ) .

We have the following integral representation of the moment sequence μ_{n}(r).

μ_{n}(r) = \frac{−1}{π} ∫_{2}^{2} x^{n} \frac{\sqrt{4 − x^{2}}(r − 1)}{2(rx^{2} + (r − 1)^{2})} dx + \frac{r + 1}{2r} (−\frac{r − 1}{\sqrt{−r}})^{n} + \frac{r + 1}{2r} (\frac{r − 1}{\sqrt{−r}})^{n} .

This shows that in this case, the measure defining the orthogonal polynomial is no longer absolutely continuous, but it takes into account the zeros of the denominator term rx^{2} + (r − 1)^{2}. Note that we have

\frac{w̃′(x)}{w̃} = \frac{x(rx^{2} − r^{2} − 6r − 1)}{(4 − x^{2})(rx^{2} + (r − 1)^{2})}

in this case.

We now move to a more general example.

Example 9. We consider the Riordan array

(\frac{1 − x − x^{2}}{1 − 3x − 4x^{2}} , \frac{x}{1 − 3x − 4x^{2}} )

which begins

( \begin{matrix}1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 2 & 1 & 0 & 0 & 0 & 0 & 0 \\ 9 & 5 & 1 & 0 & 0 & 0 & 0 \\ 35 & 28 & 8 & 1 & 0 & 0 & 0 \\ 141 & 139 & 56 & 11 & 1 & 0 & 0 \\ 563 & 670 & 339 & 93 & 14 & 1 & 0 \\ 2253 & 3129 & 1911 & 662 & 139 & 17 & 1 \end{matrix} ) ,

with inverse

(\frac{5 + 7x + \sqrt{1 + 6x + 25x^{2}}}{2(1 + x − 11x^{2}} , \frac{\sqrt{1 + 6x + 25x^{2}} − 3x − 1}{8x})

that begins

( \begin{matrix}1 & 0 & 0 & 0 & 0 & 0 & 0 \\ −2 & 1 & 0 & 0 & 0 & 0 & 0 \\ 1 & −5 & 1 & 0 & 0 & 0 & 0 \\ 13 & 12 & −8 & 1 & 0 & 0 & 0 \\ −62 & 9 & 32 & −11 & 1 & 0 & 0 \\ 97 & −217 & −43 & 61 & −14 & 1 & 0 \\ 457 & 920 & −332 & −170 & 99 & −17 & 1 \end{matrix} ) .

The moment sequence μ_{n} thus begins

1, −2, 1, 13, −62, 97, 457, ….

The inverse matrix has a production matrix that begins

( \begin{matrix}−2 & 1 & 0 & 0 & 0 & 0 & 0 \\ −3 & −3 & 1 & 0 & 0 & 0 & 0 \\ 0 & −4 & −3 & 1 & 0 & 0 & 0 \\ 0 & 0 & −4 & −3 & 1 & 0 & 0 \\ 0 & 0 & 0 & −4 & −3 & 1 & 0 \\ 0 & 0 & 0 & 0 & −4 & −3 & 1 \\ 0 & 0 & 0 & 0 & 0 & −4 & −3 \end{matrix} ) ,

and hence the generating function \frac{5+7x+\sqrt{1+6x+25x^{2}}}{2(1+x−11x^{2}} of the moment sequence has a continued fraction expression as

\frac{1}{1 + 2x + \frac{3x^{2}}{1 + 3x + \frac{4x^{2}}{1 + 3x + \frac{4x^{2}}{1 + ⋯}}}}.

This implies that the moments μ_{n} have Hankel transform given by

h_{n} = (−3)^{n}(−4)^{\binom{n}{2}}.

We find that the moments have integral representation

μ_{n} = \frac{1}{π} ∫_{−3−4i}^{−3+4i} x^{n} \frac{3i\sqrt{x^{2} + 6x + 25}}{2(x^{2} + x − 11)} dx + (\frac{3\sqrt{5}}{10} + \frac{5}{2}) (\frac{3\sqrt{5}}{2} − \frac{1}{2})^{n} .

Thus the support for the measure for the corresponding family of orthogonal polynomials has an absolutely continuous part supported by the imaginary line segment [−3−4i, −3+4i] and an atomic mass on the real axis at x = \frac{3\sqrt{5}}{2} − \frac{1}{2}. In this case we have

\frac{w̃′(x)}{w̃(x)} = −\frac{x^{3} + 9x^{2} + 64x + 58}{(x^{2} + x − 11)(x^{2} + 6x + 25)}.

The corresponding family of orthogonal polynomials P_{n}(x) satisfy the three-term recurrence

P_{n}(x) = (x + 3)P_{n−1}(x) + 4P_{n−2},

with P_{0}(x) = 1, P_{1}(x) = x + 2.

6 Complementary orthogonal polynomials

In this section, we consider the classical case where c = 0 and d = −b. In this case we have seen that the moments have generating function given by

μ(x) = \frac{1}{\sqrt{1 − 2ax + x^{2}(a^{2} − 4b)}}.

We now claim that the generating function

μ̃(x) = \frac{1}{x}Rev(xμ(x))

is the generating function of the moments for another family of orthogonal polynomials, which we will call the complementary orthogonal polynomials to the family of polynomials defined by the Riordan array

(\frac{1 − bx^{2}}{1 + ax + bx^{2}} , \frac{x}{1 + ax + bx^{2}} ) .

By solving the equation

\frac{u}{\sqrt{1 − 2au + u^{2}(a^{2} − 4b)}} = x,

and taking the appropriate branch, we find that

μ̃(x) = \frac{1}{x}Rev(xμ(x)) = \frac{\sqrt{1 + 4bx^{2}} − ax}{1 − x^{2}(a^{2} − 4b)} = \frac{1}{ax + \sqrt{1 + 4bx^{2}}} .

But now

(\frac{1 + ax + bx^{2}}{1 − bx^{2}} , \frac{x}{1 − bx^{2}} )^{−1} = (μ̃(x), xc(−bx^{2})),

where

c(x) = \frac{1 − \sqrt{1 − 4x}}{2x}

is the generating function of the Catalan numbers C_{n} = \frac{1}{n+1} \binom{2n}{n}. Thus μ̃(x) is the generating function of the moments of the family of orthogonal polynomials whose coefficient array is given by the Riordan array

R = (\frac{1 + ax + bx^{2}}{1 − bx^{2}} , \frac{x}{1 − bx^{2}} ) .

We regard the Riordan arrays R̃ = (\frac{1−bx^{2}}{1+ax+bx^{2}} , \frac{x}{1+ax+bx^{2}} ) and (\frac{1+ax+bx^{2}}{1−bx^{2}} , \frac{x}{1−bx^{2}} ) as being complementary to each other. The production matrix of the inverse matrix (μ̃(x), xc(−bx^{2})) begins

( \begin{matrix}−a & 1 & 0 & 0 & 0 & 0 & 0 \\ −2b & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & −b & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & −b & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & −b & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & −b & 0 & 1 \\ 0 & 0 & 0 & 0 & 0 & −b & 0 \end{matrix} )

and hence these complementary orthogonal polynomials Q_{n}(x) are defined by the three-term recurrence

Q_{n}(x) = xQ_{n−1} + bQ_{n−2}(x),

Q_{0}(x) = 1, Q_{1}(x) = x + a.

In fact, we have

Q_{n}(x) = P_{n}(x) + aP_{n−1}(x) + bP_{n−2}(x)

where

P_{n}(x) = (−b)^{n/2}U (n, \frac{x}{2\sqrt{−b}}) .

This follows from the the factorisation of Riordan arrays given by

(\frac{1 + ax + bx^{2}}{1 − bx^{2}} , \frac{x}{1 − bx^{2}} ) = (1 + ax + bx^{2}, x) (\frac{1}{1 − bx^{2}} , \frac{x}{1 − bx^{2}} ) ,

where the right-most matrix is closely related to the Chebyshev polynomials of the first kind (in \frac{x}{2\sqrt{−b}}).

In line with our previous results, when a = 0, we have a family of classical orthogonal polynomials with associated measure

\frac{1}{π} \frac{−1}{\sqrt{−x^{2} − 4b}}.

For instance, the polynomials with coefficient array

(\frac{1 + 2x^{2}}{1 − 2x^{2}} , \frac{x}{1 − 2x^{2}} )

have moments that begin

1, 0, −4, 0, 24, 0, −160, 0, 1120, 0, −8064, 0 …

given by

μ_{n} = \frac{1}{π} ∫_{−i2\sqrt{2}}^{i2\sqrt{2}} \frac{−x^{n}}{\sqrt{−x^{2} − 8}} dx = (−2)^{n/2}\binom{n}{\frac{n}{2}} \frac{1 + (−1)^{n}}{2} .

This is essentially \underline{A059304}. In the semi-classical case of

(\frac{1 + 2x + 2x^{2}}{1 − 2x^{2}} , \frac{x}{1 − 2x^{2}} ) ,

we find that the moments μ_{n}, which begin

1, −2, 0, 8, −8, −32, 64, 128, −416, −512, 2560, …,

have integral expression

μ_{n} = \frac{1}{π} ∫_{−2\sqrt{2}i}^{2\sqrt{2}i} \frac{−x^{n}i\sqrt{x^{2} + 8}}{x^{2} + 4} dx + i(2i)^{n}.

The Hankel transform of μ_{n} is given by

h_{n} = 2^{n}(−2)^{\binom{n+1}{2}}.

7 A note on the INVERT transform

For a given sequence with generating function g(x), the sequence with generating function

\frac{g(x)}{1 + axg(x)}

is known as the INVERT(a) transform of the first sequence, with a similar designation for the generating functions. We recall that in the last section we met the generating function

μ̃(x) = \frac{1}{ax + \sqrt{1 + 4bx^{2}}} .

We note that

μ̃(x) = \frac{\frac{1}{\sqrt{1+4bx^{2}}}}{1 + ax\frac{1}{\sqrt{1+4bx^{2}}}} .

In other words, μ̃(x) is the INVERT(a) transform of \frac{1}{\sqrt{1+4bx^{2}}} . The weight function corresponding to \frac{1}{\sqrt{1+4bx^{2}}} is given by

w(x) = \frac{1}{\sqrt{−x^{2} − 4b^{2}}} ⟹ \frac{w′(x)}{w(x)} = \frac{−x}{x^{2} + 4b},

so this is a classical case. Corresponding to μ̃(x) we get a weight with absolutely continuous part equal to

w̃(x) = \frac{(x^{2} + a^{2} + 4b)\sqrt{−x^{2} − 4b^{2}}}{(x^{2} − a^{2} + 4b)^{2}} .

We find that \frac{w̃′(x)}{w̃(x)} is of classical type only if a = 0.

It is instructive to consider some more basic examples.

Example 10. We consider w(x) = \frac{1}{π} \frac{1}{\sqrt{x(4−x)}}, which is the weight function for the family

of orthogonal polynomials whose moments are the central binomial coefficients \binom{2n}{n} with generating function \frac{1}{\sqrt{1−4x}}. We have

\frac{w′(x)}{w(x)} = \frac{2 − x}{x(x − 4)}.

The INVERT(a) transform of this generating function is given by

\frac{1}{ax + \sqrt{1 − 4x}} = \frac{ax − \sqrt{1 − 4x}}{a^{2}x^{2} + 4x − 1}.

Then the absolutely continuous part of the corresponding weight function is given by

w_{a}(x) = −\frac{1}{π} \frac{\sqrt{x(4 − x)}}{x^{2} − 4x − a^{2}} .

We have

\frac{w′_{a}(x)}{w_{a}(x)} = \frac{(2 − x)(x^{2} − 4x + a^{2})}{x(x − 4)(x^{2} − 4x − a^{2})},

which will be of classical form only if a = 0. For instance, when a = 1 we obtain moments (\underline{A081696})

1, 1, 3, 9, 29, 97, 333 …

given by

\frac{1}{π} ∫_{0}^{4} x^{n} (\frac{−\sqrt{x(4 − x)}}{x^{2} − 4x − 1}) dx + \frac{1}{\sqrt{5}}(2 − \sqrt{5})^{n}.

Similarly, for a = −1, we obtain moments that begin

1, 3, 11, 43, 173, 707, 3917, …

given by

\frac{1}{π} ∫_{0}^{4} x^{n} (\frac{−\sqrt{x(4 − x)}}{x^{2} − 4x − 1}) dx + \frac{1}{\sqrt{5}}(2 + \sqrt{5})^{n}.

Here, the zeros of the denominator x^{2} − 4x − 1 are given by 2 − \sqrt{5} and 2 + \sqrt{5}. We note that these occur outside the interval of integration [0, 4].

Example 11. In this example, we consider the generating function c(x) = \frac{1−\sqrt{1−4x}}{2x} of the Catalan numbers. The associated weight function is

w(x) = \frac{1}{2π} \frac{\sqrt{x(4 − x)}}{x},

with

\frac{w′(x)}{w(x)} = \frac{2}{x(x − 4)}.

The coefficient array of the corresponding classical orthogonal polynomials is the Riordan array

R = (\frac{1}{1 + x}, \frac{x}{(1 + x)^{2}} ) = (c(x), c(x) − 1)^{−1} = (c(x), xc(x)^{2})^{−1}.

We now consider the INVERT(a) transform of c(x). This is

\frac{c(x)}{1 + axc(x)} = \frac{1 − \sqrt{1 − 4x}}{x(2 + a − a\sqrt{1 − 4x})} = \frac{1 + 2ax − \sqrt{1 − 4x}}{2x(1 + a + a^{2}x)}.

The corresponding weight function is then given by

w_{a}(x) = \frac{\sqrt{x(4 − x)}}{2(x(a + 1) + a^{2})}.

We note that the zero of the denominator is x = \frac{−a^{2}}{a+1}. In this case, we see that \frac{w′_{a}(x)}{w_{a}(x)} is of classical form for a = −2, −1, 0, being equal to

\frac{2}{x(4 − x)}, \frac{x − 2}{x(x − 4)}, \frac{2}{x(x − 4)},

respectively. For these values of a, the roots of (1 + a)x + a^{2} = 0 are, respectively, 4, −∞, and 0. All other values lie outside the interval of integration [0, 4], and hence contribute an atomic part. The case a = 0 is just the Catalan numbers. The case a = −1 is the once shifted Catalan numbers C_{n+1},

1, 2, 5, 14, 42, 132, 429, 1430, 4862, ….

Thus

C_{n+1} = \frac{1}{2π} ∫_{0}^{4} x^{n}\sqrt{x(4 − x)} dx.

The coefficient array of the corresponding orthogonal polynomials is the Riordan array

(\frac{1}{(1 + x)^{2}} , \frac{x}{(1 + x)^{2}} ) .

The case a = −2 is the sequence \binom{2n+1}{n+1} with generating function

\frac{c(x)}{1 − 2xc(x)} = \frac{1 − \sqrt{1 − 4x}}{2x\sqrt{1 − 4x}}.

We then have

\binom{2n + 1}{n + 1} = \frac{1}{π} ∫_{0}^{4} x^{n} \frac{\sqrt{x(4 − x)}}{2(4 − x)} dx.

The coefficient array of the corresponding orthogonal polynomials is the Riordan array

(\frac{1 − x}{(1 + x)^{2}} , \frac{x}{(1 + x)^{2}} ) .

For other values of a, we have an absolutely continuous part and an atomic measure part.

μ_{n}(a) = \frac{1}{π} ∫_{0}^{n} x^{n} \frac{\sqrt{x(4 − x)}}{2((1 + a)x + a^{2})} + \frac{(a + 1)^{2} − 1}{(a + 1)^{2}} (−\frac{a^{2}}{a + 1})^{n} .

For instance, when a = 3, the moment sequence μ_{n}(3) has generating function

μ_{3}(x) = \frac{1 − \sqrt{1 − 4x}}{x(3\sqrt{1 − 4x} − 1)},

and we have

μ_{n}(3) = \frac{1}{π} ∫_{0}^{4} x^{n} \frac{x(4 − x)}{2(4x + 9)} dx + \frac{15}{16}(−\frac{9}{4})^{n}.

This sequence \underline{A049027} begins

1, 4, 17, 74, 326, 1446, 6441, 28770, 128750, 576944, ….

The corresponding Riordan array is

(\frac{1 − 2x}{(1 + x)^{2}} , \frac{x}{(1 + x)^{2}} ) .

In general, we have a one-parameter family of orthogonal polynomials whose coefficient arrays are given by

(\frac{1 − (r − 1)x}{(1 + x)^{2}}, \frac{x}{(1 + x)^{2}} ) .

The moment sequence μ_{n}(r) then has generating function

μ_{r}(x) = \frac{1}{1 − (r + 1)x − \frac{x^{2}}{1 − 2x − \frac{x^{2}}{1 − 2x − ⋯}}}.

More generally, in the context of Favard’s theorem, we can assume that

μ(x) = \frac{1}{1 − α_{0}x − \frac{β_{1}x^{2}}{1 − α_{1}x − \frac{β_{2}x^{2}}{1 − ⋯}}}.

In this case, the INVERT(a) transform has generating function

μ_{a}(x) = \frac{1}{1 − (a + α_{0})x − \frac{β_{1}x^{2}}{1 − α_{1}x − \frac{β_{2}x^{2}}{1 − ⋯}}}.

It is then clear that both sequences will have the same Hankel transforms.

8 A special product sequence

In this section, we investigate the sequence generated by

G(x) = \frac{1}{x}Rev (\frac{x}{μ̃(x)}) ,

where

μ̃(x) = \frac{1}{ax + \sqrt{1 + 4bx^{2}}} .

Thus we have

G(x) = \frac{1}{x}Rev(x(ax + \sqrt{1 + 4bx^{2}})).

We can use Lagrange Inversion to explore this quantity. Thus we have

[x^{n}]G(x) = [x^{n+1}]Rev(x(ax + \sqrt{1 + 4bx^{2}})) = \frac{1}{n + 1}[x^{n}] (\frac{1}{ax + \sqrt{1 + 4bx^{2}}} )^{n+1} .

Then we have

[x^{n}]G(x) = \frac{1}{n + 1} ∑_{i=0}^{⌊ \frac{n}{2} ⌋} \binom{2i − n − 1}{n}\binom{i − n − 1/2}{i}4^{i}b^{i}a^{n−2i}

= \frac{1}{n + 1}\binom{2n}{n}[x^{n}]\frac{1}{1 + ax + bx^{2}}
= C_{n}[x^{n}]\frac{1}{1 + ax + bx^{2}} .

In order to revert the expression x(ax + a\sqrt{1 + 4bx^{2}}), we must solve the equation

u(au + \sqrt{1 + 4bu^{2}}) = x.

Simplifying leads to the quadratic equation in u^{2}

(a^{2} − 4b)u^{4} − (1 + 2ax)u^{2} + x = 0

whose solution is given by

±\sqrt{\frac{1 + 2ax ± \sqrt{1 + 4ax + 16bx^{2}}}{2(a^{2} − 4b)}} .

Thus for instance, we find that the generating function of the product C_{n}F_{n+1} of the Catalan numbers C_{n} (\underline{A000108}) and the non-negative Fibonacci numbers F_{n+1} (\underline{A000045}) is given by (a = b = −1)

\frac{1}{x}\sqrt{\frac{1 − 2x − \sqrt{1 − 4x − 16x^{2}}}{10}} .

This is the sequence \underline{A098614} in the On-Line Encyclopedia of Integer Sequences [32, 33]. This sequence, contributed by Paul D. Hanna, begins

1, 1, 4, 15, 70, 336, 1716, 9009, 48620, ….

In similar fashion the sequence with generating function

\frac{1}{x}\sqrt{\frac{1 − 2x − \sqrt{1 − 4x − 32x^{2}}}{18}}

is given by the product of the Catalan numbers C_{n} and the Jacobsthal numbers J_{n+1} = \frac{2^{n+1}}{3} + \frac{(−1)^{n}}{3} (\underline{A001045}). This sequence (\underline{A200375}) begins

1, 1, 6, 25, 154, 882, 5676, 36465, 244530, 1657942, …

It should be noted that the On-Line Encyclopedia of Integer Sequences is a rich repository of sequences including Riordan arrays, coefficient arrays of orthogonal polynomials, and significant moment sequences.

9 Exponential Riordan arrays and classical orthogonal polynomials

An exponential Riordan array R is an invertible lower-triangular matrix defined by two power series

g(x) = 1 + g_{1} \frac{x}{1!} + g_{2} \frac{x^{2}}{2!} + ⋯

and

f(x) = \frac{x}{1!} + f_{2} \frac{x^{2}}{2!} + f_{3} \frac{x^{3}}{3!} + ⋯,

where the (n, k)-th element of the corresponding matrix is given by

r_{n,k} = \frac{n!}{k!}[x^{n}]g(x)f(x)^{k},

where [x^{n}] is the operator that extracts the coefficient of x^{n} in the power series upon which it operates [24]. Note that we have chosen g_{0} = 1 and f_{1} = 1 here, to simplify the exposition. We denote the exponential array defined by the pair g, f by [g, f]. Note that all orthogonal polynomials in this section will be monic (the coefficient of x^{n} in P_{n}(x) is 1).

In order for a Riordan array R to be the coefficient array of a family of orthogonal polynomials, we require that the production matrix P_{M} = M^{−1}\overline{M} of the inverse matrix M = R^{−1} be tri-diagonal. If M = [u, v] then this production matrix is generated by two power series, the A series and the Z series. We have

A(x) = v′(v̄(x)), Z(x) = \frac{u′(v̄(x))}{u(v̄(x))}.

The matrix P_{M} then has its bivariate generating function given by

e^{xy}(Z(x) + yA(x)).

The most general bivariate generating function of the production matrix of an exponential Riordan array M for that matrix to have a tri-diagonal production matrix is given by

e^{xy}(α + βx + y(1 + γx + δx^{2})),

where we have

Z(x) = α + βx, A(x) = 1 + γx + δx^{2}.

This leads to a production matrix that begins

( \begin{matrix}α & 1 & 0 & 0 & 0 & 0 & 0 \\ β & α + γ & 1 & 0 & 0 & 0 & 0 \\ 0 & 2(β + δ) & α + 2γ & 1 & 0 & 0 & 0 \\ 0 & 0 & 3β + 6δ & α + 3γ & 1 & 0 & 0 \\ 0 & 0 & 0 & 4β + 12δ & α + 4γ & 1 & 0 \\ 0 & 0 & 0 & 0 & 5β + 20δ & α + 5γ & 1 \\ 0 & 0 & 0 & 0 & 0 & 6β + 30δ & α + 6γ \end{matrix} ) .

In this case where the production matrix of M = R^{−1} is tri-diagonal, we call M the moment matrix of the family of orthogonal polynomials whose coefficient array is given by the Riordan array R.

We then have that R = M^{−1}, the coefficient array of the associated orthogonal polynomials, will be defined by

R = [e^{∫_{0}^{x} \frac{Z(t)}{A(t)} dt}, ∫_{0}^{x} \frac{dt}{A(t)}]
= [e^{∫_{0}^{x} \frac{α+βt}{1+γt+δt^{2}} dt}, ∫_{0}^{x} \frac{dt}{1 + γt + δt^{2}} ]

= [\frac{e^{^{\frac{βγ/δ−2α}{4δ−γ^{2}}} (2 tan^{−1}( \frac{γ+2δx}{\sqrt{4δ−γ^{2}}} )−2 sin^{−1}( \frac{γ}{2\sqrt{δ}} ))}}{(1 + γx + δx^{2})^{β/(2δ)}} , \frac{1}{4δ − γ^{2}} (2 tan^{−1} (\frac{γ + 2δx}{\sqrt{4δ − γ^{2}}} ) − 2 sin^{−1} (\frac{γ}{2\sqrt{δ}} ))] .

This is the most general form that an exponential Riordan array can have for it to be the coefficient array of a family P_{n}(x) of orthogonal polynomials. These polynomials satisfy the three-term recurrence

P_{n}(x) = (x − (α + (n − 1)γ))P_{n−1}(x) − (n − 1)(β + (n − 2)δ)P_{n−2}(x),

with P_{0}(x) = 1, P_{1}(x) = x − α.

Example 12. We let A(x) = 1 + x + x^{2}/2, Z(x) = 1 + x. We find that R is given by

R = M^{−1} = [\frac{2}{2 + 2x + x^{2}} , 2 tan^{−1}(1 + x) − \frac{π}{2} ] .

The moment matrix M is then given by

M = [(1 + sin(x)) sec(x)^{2}, tan(x) + sin(x) − 1] .

The moments in this case are the shifted Euler or up/down numbers

1, 1, 2, 5, 16, 61, 272, 1385, 7936, 50521, ….

We have

P_{n}(x) = (x − n)P_{n−1}(x) − \frac{n(n − 1)}{2}P_{n−2}(x).

The production matrix of the moment matrix begins

( \begin{matrix}1 & 1 & 0 & 0 & 0 & 0 \\ 1 & 2 & 1 & 0 & 0 & 0 \\ 0 & 3 & 3 & 1 & 0 & 0 \\ 0 & 0 & 6 & 4 & 1 & 0 \\ 0 & 0 & 0 & 10 & 5 & 1 \\ 0 & 0 & 0 & 0 & 15 & 6 \end{matrix} ) .

This indicates that the Hankel transform of the moment sequence is given by

h_{n} = ∏_{k=0}^{n} \binom{k + 2}{2}^{n−k}.

Example 13. We let A(x) = 1 + 2x + x^{2}, Z(x) = 1 + x. We find that R is given by

R = M^{−1} = [\frac{1}{1 + x}, \frac{x}{1 + x}]

with moment matrix

M = [\frac{1}{1 − x}, \frac{x}{1 − x}] .

The moments are thus n! and the polynomials are the scaled Laguerre polynomials

n! ∑_{k=0}^{n} \binom{n}{k}\frac{(−1)^{n−k}}{k!} x^{k}.

Classically, we have

n! = ∫_{0}^{∞} x^{n}e^{−x} dx.

Then w(x) = e^{−x} and \frac{w′(x)}{w(x)} = −1. The polynomials satisfy the recurrence

P_{n}(x) = (x − (2n − 1))P_{n−1}(x) − (n − 1)^{2}P_{n−1}(x).

The production matrix in this case begins

( \begin{matrix}1 & 1 & 0 & 0 & 0 & 0 \\ 1 & 3 & 1 & 0 & 0 & 0 \\ 0 & 4 & 5 & 1 & 0 & 0 \\ 0 & 0 & 9 & 7 & 1 & 0 \\ 0 & 0 & 0 & 16 & 9 & 1 \\ 0 & 0 & 0 & 0 & 25 & 11 \end{matrix} ) .

This indicates that the Hankel transform of the moment sequence is given by

h_{n} = ∏_{k=1}^{n} k^{2(n−k+1)} = ∏_{k=0}^{n} k!^{2}.

Example 14. Let Z(x) = 1 + 2x and A(x) = 1 + x + x^{2}. The exponential Riordan array whose production matrix is defined by A(x) and Z(x) is given by

M = [\frac{3}{2 (cos (\sqrt{3}x + \frac{π}{3}) + 1)}, \frac{\sqrt{3}}{2} tan (\frac{\sqrt{3}x}{2} + \frac{π}{6} ) − \frac{1}{2}] .

The (n, k)-th element of this array counts k forests of planer increasing unary-binary trees with n nodes. The production matrix of this array begins

( \begin{matrix}1 & 1 & 0 & 0 & 0 & 0 \\ 2 & 2 & 1 & 0 & 0 & 0 \\ 0 & 6 & 3 & 1 & 0 & 0 \\ 0 & 0 & 12 & 4 & 1 & 0 \\ 0 & 0 & 0 & 20 & 5 & 1 \\ 0 & 0 & 0 & 0 & 30 & 6 \end{matrix} ) .

The inverse array R = M^{−1} is given by

[\frac{1}{1 + x + x^{2}} , \frac{2}{\sqrt{3}} tan^{−1} (\frac{1 + 2x}{\sqrt{3}} − \frac{π}{3\sqrt{3}})] .

This is the coefficient array of the family of orthogonal polynomials

P_{n}(x) = (x − n)P_{n−1}(x) − n(n − 1)P_{n−2}(x),

with P_{0}(x) = 1, P_{1}(x) = x − 1.

In general, we have

h_{n} = ∏_{k=1}^{n}(k(β + (k − 1)δ))^{n−k+1}.

Proposition 15. The exponential Riordan array [\frac{1}{(1+x)^{r}} , \frac{x}{1+x}] is the coefficient array of a family of classical orthogonal polynomials. We have

w(x) = e^{−x} \frac{x^{r−1}}{(r − 1)!}

on the interval

[0, ∞).

The moments have integral representation

μ_{n} = ∫_{0}^{∞} x^{n}e^{−x} \frac{x^{r−1}}{(r − 1)!} dx = n!\binom{n + r − 1}{n} = ∏_{k=0}^{n−1} r + k.

The moments have generating function given by

μ(x) = \frac{1}{1 − rx − \frac{rx^{2}}{1 − (r + 2)x − \frac{2(r + 1)x^{2}}{1 − (r + 4)x − \frac{3(r + 2)x^{2}}{1 − (r + 6)x − ⋯}}}}.

The Hankel transform of the moments is given by

h_{n} = r^{n} ∏_{k=1}^{n} k(k(r + k))^{n−k}.

The polynomials P_{n}(x) satisfy the three-term recurrence

P_{n}(x) = (x − (r + 2(n − 1)))P_{n−1}(x) − (n − 1)(r + n − 2)P_{n−2}(x), n > 1,

with P_{0}(x) = 1, P_{1}(x) = x − r.

If y = P_{n}(x), then y satisfies the differential equation

xy′′ + (r − x)y′ + ny = 0.

Proof. The main conclusion follows from the fact that

\frac{w′(x)}{w(x)} = \frac{r − 1 − x}{x},

where

w(x) = e^{−x} \frac{x^{r−1}}{(r − 1)!}.

We note that for M = [\frac{1}{(1+x)^{r}} , \frac{x}{1+x}], we have

M^{−1} = [\frac{1}{(1 − x)^{r}} , \frac{x}{1 − x}] .

We find that

A(x) = (1 + x)^{2}, Z(x) = r(1 + x).

The corresponding production matrix P_{M^{−1}} is generated by

e^{xy}(r(1 + x) + y(1 + x)^{2}).

This matrix is therefore tri-diagonal and begins

( \begin{matrix}r & 1 & 0 & 0 & 0 & 0 & 0 \\ r & r + 2 & 1 & 0 & 0 & 0 & 0 \\ 0 & 2r + 2 & r + 4 & 1 & 0 & 0 & 0 \\ 0 & 0 & 3(r + 2) & r + 6 & 1 & 0 & 0 \\ 0 & 0 & 0 & 4(r + 3) & r + 8 & 1 & 0 \\ 0 & 0 & 0 & 0 & 5(r + 4) & r + 10 & 1 \\ 0 & 0 & 0 & 0 & 0 & 6(r + 5) & r + 12 \end{matrix} ) .

The continued fraction, the Hankel transform and three-term recurrence now follow. Proposition 16. The exponential Riordan array [\frac{1}{(1+x)^{r+1}} , \frac{x}{1+x}] is the coefficient array of a family of classical orthogonal polynomials. We have

w(x) = e^{−x} \frac{x^{r}}{r!}

on the interval

[0, ∞).

The moments have integral representation

μ_{n} = ∫_{0}^{∞} x^{n}e^{−x} \frac{x^{r}}{r!} dx = n!\binom{n + r}{r} = ∏_{k=1}^{n} r + k.

The moments have generating function given by

μ(x) = \frac{1}{1 − (r + 1)x − \frac{(r + 1)x^{2}}{1 − (r + 3)x − \frac{2(r + 2)x^{2}}{1 − (r + 5)x − \frac{3(r + 3)x^{2}}{1 − (r + 7)x − ⋯}}}}.

The Hankel transform of the moments is given by

h_{n} = ∏_{k=0}^{n}(k(r + k))^{n−k+1}.

The polynomials P_{n}(x) satisfy the three-term recurrence

P_{n}(x) = (x − (r + 2n − 1))P_{n−1}(x) − (n − 1)(r + n − 1)P_{n−2}(x), n > 1,

with P_{0}(x) = 1, P_{1}(x) = x − (r + 1).

If y = P_{n}(x), then y satisfies the differential equation

xy′′ + (r + 1 − x)y′ + ny = 0.

Proof. The main conclusion follows from the fact that

\frac{w′(x)}{w(x)} = \frac{r − x}{x},

where

w(x) = e^{−x} \frac{x^{r}}{r!} .

For M = [\frac{1}{(1+x)^{r+1}} , \frac{x}{1+x}] we have M^{−1} = [\frac{1}{(1−x)^{r+1}} , \frac{x}{1−x}]. We find that P_{M^{−1}} is generated by

e^{xy}((r + 1)(1 + x) + y(1 + x)^{2}).

This matrix is therefore tri-diagonal and begins

( \begin{matrix}r + 1 & 1 & 0 & 0 & 0 & 0 & 0 \\ r + 1 & r + 3 & 1 & 0 & 0 & 0 & 0 \\ 0 & 2r + 4 & r + 5 & 1 & 0 & 0 & 0 \\ 0 & 0 & 3(r + 3) & r + 7 & 1 & 0 & 0 \\ 0 & 0 & 0 & 4(r + 4) & r + 9 & 1 & 0 \\ 0 & 0 & 0 & 0 & 5(r + 5) & r + 11 & 1 \\ 0 & 0 & 0 & 0 & 0 & 6(r + 6) & r + 13 \end{matrix} ) .

The continued fraction, the Hankel transform and three-term recurrence now follow.

Proposition 17. The exponential Riordan array [e^{− \frac{rx^{2}}{2}} , x] is the coefficient array of a family of classical orthogonal polynomials. We have

w(x) = e^{− \frac{x^{2}}{2r}}

on the interval

(−∞, ∞).

The moments have integral representation

μ_{n} = ∫_{−∞}^{∞} x^{n}e^{− \frac{x^{2}}{2r}} dx.

These begin

1, 0, r, 0, 3r^{2}, 0, 15r^{3}, 0, 105r^{4}, 0, 945r^{5}, 0, ….

The moments have generating function given by

μ(x) = \frac{1}{1 − \frac{rx^{2}}{1 − \frac{2rx^{2}}{1 − \frac{3rx^{2}}{1 − ⋯}}}}.

The Hankel transform of the moments is given by

h_{n} = r^{\binom{n+1}{2}} ∏_{k=1}^{n} k^{n−k+1}.

The polynomials P_{n}(x) satisfy the three-term recurrence

P_{n}(x) = xP_{n−1}(x) − r(n − 1)P_{n−2}(x), n > 1,

with P_{0}(x) = 1, P_{1}(x) = x.

If y = P_{n}(x), then y satisfies the differential equation

ry′′ − xy′ + ny = 0.

Proof. The main conclusion follows from the fact that

\frac{w′(x)}{w(x)} = \frac{−x}{r},

where

w(x) = e^{− \frac{x^{2}}{r}} .

The inverse coefficient matrix [u, v] = [e^{\frac{rx^{2}}{2}} , x], and hence the production matrix is generated by e^{xy}(rx + y). This matrix begins

( \begin{matrix}0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ r & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 2r & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 3r & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 4r & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 5r & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 6r & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 7r & 0 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 8r & 0 \end{matrix} ) .

The continued fraction, the Hankel transform and three-term recurrence now follow.

The case r = 1 is related to the Hermite polynomials. The Riordan array [e^{− \frac{t^{2}}{2}} , t] is the coefficient array of the (probabilist) Hermite polynomials He(n, x) given by

He_{n}(x) = ∑_{k=0}^{n} \frac{n!}{(−2)^{\frac{n−k}{2}} k! (\frac{n−k}{2})!} \frac{1 + (−1)^{n−k}}{2} x^{k}.

The Physicists’ Hermite polynomials are given by

H_{n}(x) = n! ∑_{k=0}^{⌊ \frac{n}{2} ⌋} \frac{(−1)^{k}(2x)^{n−2k}}{k!(n − 2k)!} .

The coefficient array of the Physicists’ Hermite polynomials is the exponential Riordan array

[e^{−t^{2}}, 2t] .

We have

He_{n}(x) = 2^{− \frac{n}{2}} H_{n}(\sqrt{2}x).

We have

He_{n}(x) = xHe_{n−1}(x) − (n − 1)He_{n−2}(x),

and

H_{n}(x) = 2xH_{n−1}(x) − 2(n − 1)H_{n−1}(x).

10 Bessel and related polynomials

In this section we briefly look at two families of polynomials related to the Bessel polynomials [16]. We start with the exponential Riordan array

R = [\frac{1}{\sqrt{1 + 2x}}, \frac{x}{1 + 2x}] .

This array begins

( \begin{matrix}1 & 0 & 0 & 0 & 0 & 0 & 0 \\ −1 & 1 & 0 & 0 & 0 & 0 & 0 \\ 3 & −6 & 1 & 0 & 0 & 0 & 0 \\ −15 & 45 & −15 & 1 & 0 & 0 & 0 \\ 105 & −420 & 210 & −28 & 1 & 0 & 0 \\ −945 & 4725 & −3150 & 630 & −45 & 1 & 0 \\ 10395 & −62370 & 51975 & −13860 & 1485 & −66 & 1 \end{matrix} ) .

The production matrix of the inverse M = R^{−1} = [\frac{1}{\sqrt{1−2x}}, \frac{x}{1−2x}] of this array is tri-diagonal, beginning

( \begin{matrix}1 & 1 & 0 & 0 & 0 & 0 & 0 \\ 2 & 5 & 1 & 0 & 0 & 0 & 0 \\ 0 & 12 & 9 & 1 & 0 & 0 & 0 \\ 0 & 0 & 30 & 13 & 1 & 0 & 0 \\ 0 & 0 & 0 & 56 & 17 & 1 & 0 \\ 0 & 0 & 0 & 0 & 90 & 21 & 1 \\ 0 & 0 & 0 & 0 & 0 & 132 & 25 \end{matrix} ) .

It follows that the exponential Riordan array [\frac{1}{\sqrt{1+2x}}, \frac{x}{1+2x}] is the coefficient array of the family of orthogonal polynomials P_{n}(x) defined by the three-term recurrence

P_{n}(x) = (x − (4n − 3))P_{n−1}(x) − 2(n − 1)(2n − 3)P_{n−2}(x),

with P_{0}(x) = 1, P_{1}(x) = x − 1.

The moments μ_{n} of these orthogonal polynomials are the double factorials of the odd numbers (2n − 1)!! with exponential generating function \frac{1}{\sqrt{1−2x}} that begin (\underline{A001147})

1, 1, 3, 15, 105, 945, 10395, 135135, 2027025, 34459425, ….

We have

(2n − 1)!! = \frac{n!}{2π} ∫_{0}^{2} \frac{x^{n}}{\sqrt{x(2 − x)}} dx = \frac{1}{\sqrt{2π}} ∫_{0}^{∞} x^{n} \frac{e^{− \frac{x}{2}}}{\sqrt{x}} dx.

The Hankel transform of the moments μ_{n} is then given by

h_{n} = ∏_{k=0}^{n}(2(k + 1)(2k + 1))^{n−k}.

We have

\frac{w′(x)}{w(x)} = −\frac{1 + x}{x},

where

w(x) = \frac{1}{\sqrt{2π}} \frac{e^{− \frac{x}{2}}}{\sqrt{x}} .

The inverse matrix M = [\frac{1}{\sqrt{1−2x}}, \frac{x}{1−2x}] is the coefficient array of the related polynomial family P^{∗}(n, x) = (−1)^{n}P_{n}(−x) that satisfies the three-term recurrence

P_{n}^{∗}(x) = (x + (4n − 3))P_{n−1}^{∗}(x) − 2(n − 1)(2n − 3)P_{n−2}^{∗}(x),

with P_{0}^{∗}(x) = 1, P_{1}^{∗}(x) = x + 1.

We note further that the exponential Riordan array

M̃ = [\frac{d}{dx} \frac{1}{\sqrt{1 − 2x}}, \frac{x}{1 − 2x}] = [\frac{1}{(1 − 2x)^{3/2}} , \frac{x}{1 − 2x}]

has

Z(x) = 3 + 6x, A(x) = (1 + 2x)^{2}.

Its production matrix thus begins

( \begin{matrix}3 & 1 & 0 & 0 & 0 & 0 & 0 \\ 6 & 7 & 1 & 0 & 0 & 0 & 0 \\ 0 & 20 & 11 & 1 & 0 & 0 & 0 \\ 0 & 0 & 42 & 15 & 1 & 0 & 0 \\ 0 & 0 & 0 & 72 & 19 & 1 & 0 \\ 0 & 0 & 0 & 0 & 110 & 23 & 1 \\ 0 & 0 & 0 & 0 & 0 & 156 & 27 \end{matrix} ) .

The array M̃ is thus the inverse of the coefficient array of the polynomials that satisfy the recurrence relation

P_{n}(x) = (x − (4n − 1))P_{n−1}(x) − (n − 1)(4n − 2)P_{n−2}(x),

with P_{0}(x) = 1, P_{1}(x) = x − 3. We have

μ̃_{n} = \frac{1}{\sqrt{2π}} ∫_{0}^{∞} x^{n}x\frac{e^{− \frac{x}{2}}}{\sqrt{x}} dx.

Thus

w̃(x) = \frac{1}{\sqrt{2π}} \sqrt{x}e^{− \frac{x}{2}} ,

and

\frac{w̃′(x)}{w̃(x)} = \frac{1 − x}{2x}.

The Hankel transform of the moments μ̃_{n}, which begin

1, 3, 15, 105, 945, 10395, 135135, 2027025, 34459425, …,

is given by

h_{n} = ∏_{k=0}^{n}(2(k + 1)(2k + 3))^{n−k}.

The exponential Riordan array [e^{\frac{x^{2}}{2}} , x] is the coefficient array of the family of orthogonal polynomials P_{n}(x) that satisfy the recurrence

P_{n}(x) = xP_{n−1}(x) − (n − 1)P_{n−2}(x),

with P_{0}(x) = 1, P_{1}(x) = x. The moments of these polynomials are the aerated double factorials

1, 0, 1, 0, 3, 0, 15, 0, 105, 0, 945, …

with exponential generating function e^{x^{2}/2}. The Hankel transform of these moments is given by

h_{n} = ∏_{k=0}(k + 1)^{n−k} = ∏_{k=0}^{n} k!

This follows since for the moment matrix [e^{x^{2}/2}, x] we have A(x) = 1 and Z(x) = x.

We now consider a family of Bessel polynomials. The exponential Riordan array B = [\frac{1}{\sqrt{1−2x}}, 1 − \sqrt{1 − 2x}] is the coefficient array of the reverse Bessel polynomials θ_{n}(x) [16] with general term

θ_{n}(x) = ∑_{k=0}^{n} \frac{(2n − k)!}{2^{n−k}k!(n − k)!}x^{k}.

These polynomials are orthogonal on the circle, satisfying the recurrence

θ_{n}(x) = (2n − 1)θ_{n−1}(x) + x^{2}θ_{n−2}(x).

We calculate RB. We get

RB = [\frac{1}{\sqrt{1 + 2x}}, \frac{x}{1 + 2x}] · [\frac{1}{\sqrt{1 − 2x}}, 1 − \sqrt{1 − 2x}]

= [\frac{1}{\sqrt{1 + 2x}} \frac{1}{\sqrt{1 − 2\frac{x}{1+2x}}} , 1 − \sqrt{1 − 2\frac{x}{1 + 2x}}]

= [1, 1 − \frac{1}{\sqrt{1 + 2x}}] .

This last array has general term

T_{n,k} = \frac{(n − 1)!}{(k − 1)!} ∑_{j=0}^{2n} \binom{2n}{j}\binom{2n − k − j − 1}{n − k − j}\frac{(−1)^{j}}{2^{n−k−j}} .

We have

B = R^{−1} · [1, 1 − \frac{1}{\sqrt{1 + 2x}}] .

Now R^{−1} is the coefficient array of the polynomials (−1)^{n}P_{n}(−x). Thus the last two expressions allow us to express the Bessel polynomials θ_{n}(x) in terms of the polynomials P_{n}(x).

We note further that if T_{n,k} is the general term of the coefficient array [e^{− \frac{x^{2}}{2}} , x], then

T_{2n−k,k} is the general term of the coefficient array [\frac{1}{\sqrt{1+2x}}, \sqrt{1 + 2x} − 1]. More generally, if T_{n,k} is the general term of the coefficient array [e^{− \frac{rx^{2}}{2}} , x], then T_{2n−k,k} is the general term

of the coefficient array [\frac{1}{\sqrt{1+2xr}}, \sqrt{1 + 2xr} − 1].

The matrix

R = [\frac{1}{\sqrt{1 + 2x}}, \frac{x}{1 + 2x}]

is a member of a one-parameter family of coefficient arrays of orthogonal polynomials. The general element of this family is

[\frac{1}{(1 + rx)^{1/r}} , \frac{x}{1 + rx}] = [\frac{1}{(1 − rx)^{1/r}} , \frac{x}{1 − rx}]^{−1} .

The family of orthogonal polynomials P_{n}(x; r) defined by this matrix satisfies the three-term recurrence

P_{n}(x; r) = (x − (2(n − 1)r + 1))P_{n−1}(x; r) − (n − 1)r((n − 2)r + 1)P_{n−2}(x; r),

with P_{0}(x; r) = 1 and P_{1}(x; r) = x−1. The moments of this family of orthogonal polynomials have exponential generating function \frac{1}{(1−rx)^{1/r}}. They begin

1, 1, r + 1, (r + 1)(2r + 1), (r + 1)(2r + 1)(3r + 1), …,

or

μ_{n}(r) = ∏_{k=0}^{n−1} kr + 1.

The Hankel transform of μ_{n} is then

h_{n} = r^{\binom{n+1}{2}} ∏_{k=0}^{n} k!(kr + 1)^{n−k}.

The corresponding family of generalized Bessel polynomials will then have coefficient array given by

[\frac{1}{(1 − rx)^{1/r}} , 1 − (1 − rx)^{1/r}] = [1 − x, \frac{1 − (1 − x)^{r}}{r} ]^{−1} .

11 Conclusion

While Riordan arrays can only define a limited number of families of orthogonal polynomials, they offer fresh perspectives on these families and their study can lead to interesting results. For instance, solutions to the restricted Toda chain equations are provided by the Jacobi parameters of exponential Riordan arrays that are the moment arrays of families of paramterised orthogonal polynomials [7]. Thus links exist between orthogonal polynomials defined by appropriate Riordan arrays and integrable systems. Further research is warranted, particularly in the realm of exponential Riordan arrays and q-Riordan arrays, and their interaction with appropriate families of orthogonal polynomials.

12 Appendix — The Stieltjes transform of a measure

The Stieltjes transform of a measure μ on ℝ is a function G_{μ} defined on ℂ \ ℝ by

G_{μ}(z) = ∫_{ℝ} \frac{1}{z − t}μ(t).

If f is a bounded continuous function on ℝ, we have

∫_{ℝ} f(x)μ(x) = − lim_{y→0^{+}} ∫_{ℝ} f(x)ℑG_{μ}(x + iy)dx.

If μ has compact support, then G_{μ} is holomorphic at infinity and for large z,

G_{μ}(z) = ∑_{n=0}^{∞} \frac{a_{n}}{z^{n+1}} ,

where a_{n} = ∫_{ℝ} t^{n}μ(t) are the moments of the measure. If μ(t) = dψ(t) = ψ′(t)dt then (Stieltjes-Perron)

ψ(t) − ψ(t_{0}) = −\frac{1}{π} lim_{y→0^{+}} ∫_{t_{0}}^{t} ℑG_{μ}(x + iy)dx.

If now g(x) is the generating function of a sequence a_{n}, with g(x) = ∑_{n=0}^{∞} a_{n}x^{n}, then we can define

G(z) = \frac{1}{z}g (\frac{1}{z}) = ∑_{n=0}^{∞} \frac{a_{n}}{z^{n+1}} .

By this means, under the right circumstances we can retrieve the density function for the measure that defines the elements a_{n} as moments.

Example 18. We let g(z) = \frac{1−\sqrt{1−4z}}{2z} be the generating function of the Catalan numbers. Then

G(z) = \frac{1}{z}g (\frac{1}{z}) = \frac{1}{2} (1 − \sqrt{\frac{x − 4}{x}}) .

Then

ℑG_{μ}(x + iy) = −\frac{\sqrt{2}\sqrt{\sqrt{x^{2} + y^{2}}\sqrt{x^{2} − 8x + y^{2} + 16} − x^{2} + 4x − y^{2}}}{4\sqrt{x^{2} + y^{2}}} ,

and so we obtain

ψ′(x) = −\frac{1}{π} lim_{y→0^{+}} \{−\frac{\sqrt{2}\sqrt{\sqrt{x^{2} + y^{2}}\sqrt{x^{2} − 8x + y^{2} + 16} − x^{2} + 4x − y^{2}}}{4\sqrt{x^{2} + y^{2}}} \}

= \frac{1}{2π} \frac{\sqrt{x(4 − x)}}{x}.

References

[1] M. Agida and A. S. Kumar, A Boubaker polynomials expansion scheme solution to random Love’s equation in the case of a rational kernel, Elec. J. Theoretical Phys., 7 (2010), 319–326.

[2] O. B. Awojoyogbe and K. Boubaker, A solution to Bloch NMR flow equations for the analysis of hemodynamic functions of blood flow system using m-Boubaker polynomials, Current Applied Physics, 9 (2009), 278–283.

[3] P. Barry, On the restricted Chebyshev-Boubaker polynomials, Integral Transforms Spec. Funct., 28 (2017), 1–16.

[4] P. Barry, Riordan Arrays: a Primer, Logic Press, 2017.

[5] P. Barry, Riordan arrays, orthogonal polynomials as moments, and Hankel transforms, J. Integer Sequences, 14 (2011), Article 11.2.2.

[6] P. Barry and A. Hennessy, Meixner-type results for Riordan arrays and associated integer sequences, J. Integer Sequences, 13 (2010), Article 10.9.4.

[7] P. Barry, The restricted Toda chain, exponential Riordan arrays, and Hankel transforms, J. Integer Sequences, 13 (2010), Article 10.8.4.

[8] K. Boubaker and L. Zhang, Fermat-linked relations for the Boubaker polynomial sequences via Riordan matrices analysis, J. Assoc. Arab Univ. Basic Appl. Sci., 12 (2012), 74–78.

[9] K. Boubaker, Analytical initial-guess-free solution to Kepler’s transcendental equation using Boubaker polynomials expansion scheme BPES Apeiron: Studies in Infinite Nature, 17 (2010), 1–12.

[10] K. Boubaker, A. Chaouachi, M. Amlouk, and H. Bouzouita, Enhancement of pyrolysis spray dispersal performance using thermal time-response to precursor uniform deposition, Eur. Phys. J. AP, 37 (2007), 105–109.

[11] T. S. Chihara, An Introduction to Orthogonal Polynomials, Dover Publicatons, 2011.

[12] O. M. Dada, O. B. Awojoyogbe, M. Agida, and K. Boubaker, Variable separation and Boubaker polynomial expansion scheme for solving the neutron transport equation, Physics International, 2 (2011), 25–30.

[13] E. Deutsch, L. Ferrari, and S. Rinaldi, Production matrices, Adv. in Appl. Math., 34 (2005), 101–122.

[14] E. Deutsch, L. Ferrari, and S. Rinaldi, Production matrices and Riordan arrays, Ann. Comb., 13 (2009), 65–85.

[15] W. Gautschi, Orthogonal Polynomials: Computation and Approximation, Clarendon Press, 2004.

[16] E. Grosswald, Bessel Polynomials, Springer, 1978.

[17] P. Henrici, Applied and Computational Complex Analysis, Vol. 3, John Wiley & Sons, 1993.

[18] C. Krattenthaler, Advanced determinant calculus, in D. Foata, and G.-N. Han, eds., The Andrews Festschrift: Seventeen Papers on Classical Number Theory and Combinatorics, Springer, 2001, pp. 349–426.

[19] H. Labadiah and K. Boubaker, A Sturm-Liouville shaped characteristic differential equation as a guide to establish a quasi-polynomial expression to the Boubaker polynomials, J. of Diff. Eq. and C. Proc., (2007) 117-133.

[20] H. Labadiah, M. Dada, B. Awojoyogbe, B. Mahmoud, and A. Bannour, Establishment of an ordinary generating function and a Christoffel-Darboux type first-order differential equation for the heat equation related Boubaker-Turki polynomials, J. of Diff. Eq. and C. Proc., (2008), 51–66.

[21] J. W. Layman, The Hankel transform and some of its properties, J. Integer Sequences, 4 (2001), Article 01.1.5.

[22] J. C. Mason and D. C. Handscomb, Chebyshev Polynomials, Chapman and Hall/CRC, 2002.

[23] D. Merlini, R. Sprugnoli, and M. C. Verri, Lagrange inversion: when and how, Acta Appl. Math., 94 (2006), 233–249.

[24] D. Merlini, R. Sprugnoli, and M. C. Verri, The method of coefficients, Amer. Math. Monthly, 114 (2007), 40–57.

[25] G. V. Milanovic and D. Joksimovic, Properties of Boubaker polynomials and an application to Love’s integral equation, Appl. Math. Comput., 224 (2013), 74–87.

[26] T. D Noe, On the divisibility of generalized central trinomial coefficients, J. Integer Sequences, 9 (2006), Article 06.2.7.

[27] P. Peart and W.-J. Woan, Generating functions via Hankel and Stieltjes matrices, J. Integer Sequences, 3 (2000), Article 00.2.1.

[28] L. W. Shapiro, S. Getu, W.-J. Woan, and L. C. Woodson, The Riordan group, Discr. Appl. Math. 34 (1991), 229–239.

[29] L. Shapiro, A survey of the Riordan group, available electronically at , Center for Combinatorics, Nankai University, 2005.

[30] L. W. Shapiro, Bijections and the Riordan group, Theoret. Comput. Sci. 307 (2003), 403–413.

[31] S. Slama, J. Bessrour, K. Boubaker, and M. Bouhafs, A dynamical model for investigation of a 3 point maximal spatial evolution during resistance spot welding using Boubaker polynomials, Eur. Phys. J. AP, 44 (2008), 317–322.

[32] N. J. A. Sloane, The On-Line Encyclopedia of Integer Sequences. Published electroni-

    cally at http://oeis.org, 2016.

[33] N. J. A. Sloane, The on-line encyclopedia of integer sequences, Notices Amer. Math. Soc., 50 (2003), 912–915.

[34] R. Sprugnoli, Riordan arrays and combinatorial sums, Discrete Math. 132 (1994), 267– 290.

[35] P. K. Suetin, Classical orthogonal polynomials, entry in the on-line Encyclopedia of Mathematics, 2016. Available at

    https://www.encyclopediaofmath.org/index.php/Classical_orthogonal_polynomials.

[36] G. Szegö, Orthogonal Polynomials, 4th edition, Amer. Math. Soc., 1975.

[37] G. Viennot, Une théorie combinatoire des polynômes orthogonaux généraux, UQAM, Montreal, Quebec, 1983.

[38] H. S. Wall, Analytic Theory of Continued Fractions, AMS Chelsea Publishing, 2000.

[39] E. P. Wigner, Characteristic vectors of bordered matrices with infinite dimensions, Ann. of Math., 62 (1955), 548–564.

[40] E. P. Wigner, On the distribution of the roots of certain symmetric matrices, Ann. of Math., 67 (1958), 325–327.

[41] T. G. Zhao, L. Naing, and W. X. Yue, Some new features of the Boubaker polynomials expansion scheme BPES, J. Assoc. Arab Univ. Basic Appl. Sci., 12 (2012) 74–78.

2010 Mathematics Subject Classification: Primary 15B36; Secondary 33C45, 11B83, 11C20, 05A15. Keywords: Riordan array, classical orthogonal polynomial, Chebyshev polynomial, Bessel polynomial, moment sequence.

(Concerned with sequences \underline{A000045}, \underline{A000108}, \underline{A000984}, \underline{A001045}, \underline{A001147}, \underline{A049027}, \underline{A059304}, \underline{A081696}, \underline{A098614}, and \underline{A200375}.)

Received January 18 2017; revised versions received December 29 2017; December 30 2017. Published in Journal of Integer Sequences, January 21 2018.

Return to Journal of Integer Sequences home page.