\documentclass[12pt,reqno]{article} \usepackage[usenames]{color} \usepackage{amssymb} \usepackage{graphicx} \usepackage{amscd} \usepackage[colorlinks=true, linkcolor=webgreen, filecolor=webbrown, citecolor=webgreen]{hyperref} \definecolor{webgreen}{rgb}{0,.5,0} \definecolor{webbrown}{rgb}{.6,0,0} \usepackage{color} \usepackage{fullpage} \usepackage{float} \usepackage{psfig} \usepackage{graphics,amsmath,amssymb} \usepackage{amsthm} \usepackage{amsfonts} \usepackage{latexsym} \usepackage{epsf} \setlength{\textwidth}{6.5in} \setlength{\oddsidemargin}{.1in} \setlength{\evensidemargin}{.1in} \setlength{\topmargin}{-.1in} \setlength{\textheight}{8.4in} \newcommand{\seqnum}[1]{\href{http://oeis.org/#1}{\underline{#1}}} \begin{document} \begin{center} \epsfxsize=4in \leavevmode\epsffile{logo129.eps} \end{center} \theoremstyle{plain} \newtheorem{theorem}{Theorem} \newtheorem{corollary}[theorem]{Corollary} \newtheorem{lemma}[theorem]{Lemma} \newtheorem{proposition}[theorem]{Proposition} \theoremstyle{definition} \newtheorem{definition}[theorem]{Definition} \newtheorem{example}[theorem]{Example} \newtheorem{conjecture}[theorem]{Conjecture} \theoremstyle{remark} \newtheorem{remark}[theorem]{Remark} \begin{center} \vskip 1cm{\LARGE\bf Classical and Semi-Classical Orthogonal \\ \vskip 0.01in Polynomials Defined by Riordan Arrays, \\ \vskip 0.06in and Their Moment Sequences} \vskip 1cm \large Paul Barry and Arnauld Mesinga Mwafise\\ School of Science\\ Waterford Institute of Technology\\ Ireland\\ \href{mailto:pbarry@wit.ie}{\tt pbarry@wit.ie} \end{center} \vskip .2 in \begin{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. \end{abstract} \section{Introduction} Riordan arrays \cite{Book, SGWW, Survey, Shapiro_bij, Spru} 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. \begin{enumerate} \item This Introduction \item Preliminaries on ordinary Riordan arrays and orthogonal polynomials \item Classical and semi-classical orthogonal polynomials \item The classical orthogonal polynomials defined by ordinary Riordan arrays \item The non-classical case \item Complementary orthogonal polynomials defined by Riordan arrays \item A note on the INVERT transform \item A special product sequence \item Exponential Riordan arrays and classical orthogonal polynomials \item Bessel and related polynomials \item Conclusion \item Declaration \item Appendix --- The Stieltjes transform of a measure. \end{enumerate} We recall the following well-known results (the first is known as ``Favard's Theorem''), which we essentially reproduce from \cite{Kratt}, 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. \begin{theorem} \label{ThreeT}\cite{Kratt} (Cf.~\cite[Th\'eor\`eme $9$, p.~I-4]{Viennot} or \cite[Theorem 50.1]{Wall}). Let $(p_n(x))_{n\ge 0}$ be a sequence of monic polynomials, the polynomial $p_n(x)$ having degree $n=0,1,\ldots$ Then the sequence $(p_n(x))$ is (formally) orthogonal if and only if there exist sequences $(\alpha_n)_{n\ge 0}$ and $(\beta_n)_{n\ge 1}$ with $\beta_n \neq 0$ for all $n\ge 1$, such that the three-term recurrence $$p_{n+1}=(x-\alpha_n)p_n(x)-\beta_n p_{n-1}(x), \quad \text{for}\quad n\ge 1, $$ holds, with initial conditions $p_0(x)=1$ and $p_1(x)=x-\alpha_0$. \end{theorem} \begin{theorem} \label{CF} \cite{Kratt} (Cf.~\cite[Prop.~1 (7), p.~V-5]{Viennot} or \cite[Theorem~51.1]{Wall}). Let $(p_n(x))_{n\ge 0}$ be a sequence of monic polynomials, which is orthogonal with respect to some functional $\mathcal{L}$. Let $$p_{n+1}=(x-\alpha_n)p_n(x)-\beta_n p_{n-1}(x), \quad \text{for}\quad n\ge 1, $$ be the corresponding three-term recurrence which is guaranteed by Favard's theorem. Then the generating function $$g(x)=\sum_{k=0}^{\infty} \mu_k x^k $$ for the moments $\mu_k=\mathcal{L}(x^k)$ satisfies $$g(x)=\cfrac{\mu_0}{1-\alpha_0 x- \cfrac{\beta_1 x^2}{1-\alpha_1 x - \cfrac{\beta_2 x^2}{1-\alpha_2 x - \cfrac{\beta_3 x^2}{1-\alpha_3 x -\cdots}}}}.$$ \end{theorem} The {\it 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,\dots \}$ where $h_{n}=|a_{i+j}|_{i,j=0}^{n}$, i.e \begin{center} \begin{equation} \label{gen1} A=\{a_n\}_{n\in\mathbb N_0}\quad \rightarrow \quad h=\{h_n\}_{n\in\mathbb N_0}:\quad h_n=\left| \begin{array}{ccccc} a_0\ & a_1\ & \cdots & a_n & \\ a_1\ & a_2\ & & a_{n+1} \\ \vdots & & \ddots & \\ a_n\ & a_{n+1}\ & & a_{2n} \end{array} \right|. \end{equation} \end{center} 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)=\cfrac{a_0}{1-\alpha_0 x- \cfrac{\beta_1 x^2}{1-\alpha_1 x- \cfrac{\beta_2 x^2}{1-\alpha_2 x- \cfrac{\beta_3 x^2}{1-\alpha_3 x-\cdots}}}}$$ then we have \cite{Kratt} \begin{equation}\label{Kratt} h_n = a_0^{n+1} \beta_1^n\beta_2^{n-1}\cdots \beta_{n-1}^2\beta_n=a_0^{n+1}\prod_{k=1}^n \beta_k^{n+1-k}.\end{equation} Note that this is independent of $\alpha_n$. \section{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 + \cdots$$ and $$f(x)=x+f_2 x^2+ f_3 x^3 + \cdots,$$ 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 \cite{MC}. 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 \emph{ordinary} form $g(x)=\sum_{n=0}^{\infty} g_n x^n$ and $f(x)=\sum_{n=0}^{\infty} 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}{\bar{f}(x)},$$ and $$Z(x)=\frac{1}{\bar{f}(x)}\left(1-\frac{1}{g(\bar{f}(x))}\right).$$ Note that the notation $\bar{f}(x)$ denotes the compositional inverse of the power series $f$. Thus we have $$\bar{f}(f(x))=x, \quad\text{and}\quad f(\bar{f}(x))=x.$$ We shall also use the notation $\bar{f}(x)=\text{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$ \cite{ProdMat_0, ProdMat, P_W}. It is equal to the matrix $$M^{-1}\bar{M}$$ where $\bar{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}=\left(1-x \frac{Z(x)}{A(x)}, \frac{x}{A(x)}\right).$$ If $Z(x)=\gamma + \delta x$, and $A(x)=1+\alpha x + \beta x^2$ (where we assume that $\delta \ne 0$ and $\beta \ne 0$), then $P$ will be tri-diagonal. The matrix $P$ begins $$\left( \begin{array}{cccccccc} \gamma & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ \delta & \alpha & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & \beta & \alpha & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & \beta & \alpha & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & \beta & \alpha & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & \beta & \alpha & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & \beta & \alpha & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & \beta & \alpha \\ \end{array} \right).$$ In this case, we have $$M^{-1}=\left(1-x \frac{\gamma+\delta x}{1+\alpha x + \beta x^2}, \frac{x}{1+\alpha x + \beta x^2}\right)= \left(\frac{1+(\alpha-\gamma)x+(\beta-\delta)x^2}{1+\alpha x+\beta x^2},\frac{x}{1+\alpha x + \beta x^2}\right),$$ and $R=M^{-1}$ will be the coefficient array of a family of orthogonal polynomials $P_n(x)$ \cite{Meixner} where $$P_n(x)=\sum_{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 \emph{moment matrix} of the family of orthogonal polynomials $P_n(x)$ \cite{Barry_moments}. 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)) \cdot (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)) \cdot 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} = \left(\frac{1}{g(\bar{f}(x))}, \bar{f}(x)\right).$$ The identity element is $(1,x)$, which is represented by the usual identity matrix. \section{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)^{\alpha}(1+x)^{\beta}$ on $[-1,1]$, $w_L(x)=x^{\alpha}e^{-x}$ on $[0,\infty)$, and $w_H(x)=e^{-x^2}$ on $(-\infty, \infty)$, respectively. In particular, these orthogonal polynomials \cite{Chihara, Gautschi, Szego} are associated with measures that are absolutely continuous. We have $$\frac{w'_J(x)}{w_J(x)}=\frac{x(\alpha+\beta)+\alpha-\beta}{x^2-1},$$ $$ \frac{w'_L(x)}{w_L(x)}=\frac{\alpha-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 \emph{classical} \cite{Suetin} 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 $\text{deg}(V)>2$ and/or $\text{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 \emph{monic} (the coefficient of $x^n$ in $P_n(x)$ is $1$). \section{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=\left(\frac{1+ cx + dx^2}{1+ax+ bx^2}, \frac{x}{1+ax +b x^2}\right).$$ 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{equation*}\begin{split} \biggl(-\frac{(b-d)\sqrt{1-2ax+x^2(a^2-4b)}+x(a(b + d)-2bc)-b-d}{2(x^2(a^2d - 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}\biggr). \end{split}\end{equation*} The first element of this array is the generating function $\mu(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^2c + a(3b + 3c^2 - 3d) - c(2b + c^2 - 2d),\ldots$$ and their generating function is given by $$\mu(x)=\cfrac{1}{1-(a-c)x- \cfrac{(b-d)x^2}{1-ax- \cfrac{bx^2}{1-ax- \cfrac{bx^2}{1-\cdots}}}}.$$ From this we can see that the Hankel transform \cite{Kratt, Layman} 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 \cite{Barry_moments, Henrici} (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 \pi} \frac{(b-d)\sqrt{4b-(x-a)^2}}{ dx^2+x(c(b+d)-2ad)+a^2d-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^2d-b^2-b(c^2+6d)-d^2)-a^3d+a(b^2+b(c^2+6d)+d^2)-4bc(b+d)} {((x-a)^2-4b)(dx^2+x(c(b+d)-2ad)+a^2d-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 \emph{semi-classical}. Inspection of the above ratio allows us to announce the following results. \begin{proposition} The ordinary Riordan array $$\left(\frac{1+cx +dx^2}{1+ax+bx^2}, \frac{x}{1+ax + bx^2}\right)$$ defines a family of classical orthogonal polynomials in the case that either $c=d=0$ or $c=0, d=-b$.\end{proposition} \begin{corollary} When $c=d=0$, we have $$w(x)=\frac{1}{2\pi} \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 $\mu_n$ have integral representation $$\mu_n = \frac{1}{2\pi} \int_{a-2 \sqrt{b}}^{a+2 \sqrt{b}} x^n \frac{\sqrt{4b-(x-a)^2}}{b}\,dx.$$ The moments have generating function $$\mu(x)=\frac{1-ax-\sqrt{(1-ax)^2-4bx^2}}{2bx^2}$$ given by $$\mu(x)=\cfrac{1}{1-ax- \cfrac{bx^2}{1-ax- \cfrac{bx^2}{1-ax- \cfrac{bx^2}{1-\cdots}}}}.$$ By an application of Lagrange inversion \cite{Merlini_When}, we obtain \begin{eqnarray*} \mu_n &=& \frac{1}{n+1}[x^n](1+ax+bx^2)^{n+1}\\ &=& \frac{1}{n+1} \sum_{k=0}^n \binom{n+1}{j}\binom{j}{n-j}a^{2j-n}b^{n-j}\\ &=& \frac{1}{n+1} \sum_{k=0}^n \binom{n+1}{n-k}\binom{n-k}{k} a^{n-2k}b^k.\end{eqnarray*} 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),\quad 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.$$ \end{corollary} The form of the differential equation satisfied by these polynomials follows from the form of the weight function \cite{Suetin}. For a proof of the other statements, the reader is referred to \cite{Meixner}. \begin{corollary} When $c=0$ and $d=-b$, we have $$w(x)=\frac{1}{\pi} \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 $\mu_n$ have integral representation $$\mu_n=\frac{1}{\pi} \int_{a-2 \sqrt{b}}^{a+2 \sqrt{b}}x^n \frac{1}{\sqrt{4b-(x-a)^2}}\,dx.$$ The moments have generating function $$\mu(x)=\frac{1}{\sqrt{(1-ax)^2-4bx^2}}$$ given by $$\mu(x)=\cfrac{1}{1-ax- \cfrac{2bx^2}{1-ax- \cfrac{bx^2}{1-ax- \cfrac{bx^2}{1-\cdots}}}}.$$ We have the closed form expression for the moments \begin{eqnarray*}\mu_n&=& \sum_{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} \sum_{k=0}^n \binom{2n-2k}{n-k}\binom{2k}{k}(a+2\sqrt{b})^k (a-2 \sqrt{b})^{n-k}.\end{eqnarray*} 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),\quad 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.$$ \end{corollary} 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}=\sum_{n=0}^{\infty} \mathcal{P}_n(x)t^n $$ of the Legendre polynomials $\mathcal{P}_n(x)$. Then we get \cite{Noe} $$\mu_n=(a^2-4b)^{n/2} \mathcal{P}_n\left(\frac{a}{\sqrt{a^2-4b}}\right).$$ For the next result, we note that $$C_n=\frac{1}{n+1} \binom{2n}{n}$$ is the $n$-th Catalan number \seqnum{A000108}. The generating function of the Catalan numbers is given by $$c(x)=\frac{1-\sqrt{1-4x}}{2x}.$$ \begin{proposition} The ordinary Riordan array $\left(\frac{1}{1+ax}, \frac{x}{(1+ax)^2}\right)$ ($a \ne 0$) is the coefficient array of a family of classical orthogonal polynomials. We have $$w(x)=\frac{1}{2 \pi} \frac{\sqrt{x(4 a-x)}}{2 a x}$$ on the interval $$[0, 4a],$$ with $$\frac{w'(x)}{w(x)}=\frac{2a}{x(x-4a)}.$$ The moments have integral representation $$\mu_n = \frac{1}{2 \pi} \int_0^{4a} x^n \frac{\sqrt{x(4 a-x)}}{2 a x}\,dx=a^n C_n.$$ The moments have generating function $$\mu(x)=\frac{1-\sqrt{1-4ax}}{2a x},$$ with $$\mu(x)=\cfrac{1}{1-\cfrac{ax}{1-\cfrac{ax}{1-\cfrac{ax}{1-\cdots}}}},$$ or equivalently, $$\mu(x)= \cfrac{1}{1-a x - \cfrac{a^2 x^2}{1-2 a x- \cfrac{a^2 x^2}{1-2 ax-\cdots}}}.$$ The moments $\mu_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} \quad 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.$$ \end{proposition} \begin{example} The Riordan array $\left(\frac{1}{1+x^2}, \frac{x}{1+x^2}\right)$ 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,\ldots$$ given by $$\mu_n=\frac{1}{2 \pi} \int_{-2}^{2} \sqrt{4-x^2}\,dx.$$ This is closely related to Wigner's semicircle distribution \cite{Wigner_1, Wigner_2}. Note that if we define $$P_{n-1}^{(1)}=\frac{1}{2\pi}\int_{-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)=\sum_{k=0}^{\lfloor \frac{n}{2} \rfloor} \binom{n-k}{k} (-1)^k (2x)^{n-2k}$$ are the Chebyshev polynomials of the second kind \cite{Chebyshev}. \end{example} \begin{example} The Riordan array $\left(\frac{1-x^2}{1+x^2}, \frac{x}{1+x^2}\right)$ 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 (\seqnum{A000984}) $$1, 0, 2, 0, 6, 0, 20, 0, 70, 0, 252,\ldots$$ given by $$\frac{1}{\pi} \int_{-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}{\pi}\int_{-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).$$ \end{example} \section{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{\tilde{w}'}{\tilde{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'' \cite{CB}. These are ordinary Riordan arrays of the form $$\left(\frac{1+rx^2}{1+x^2}, \frac{x}{1+x^2}\right).$$ 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)=\sum_{k=0}^{\lfloor \frac{n}{2} \rfloor} \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 \cite{Agida, NMR, Boubaker, BPES, Boubaker07, Dada, Labadiah1, Labadiah2, Milanovic, Slama, Zhao}. The moments $\mu_n(r)$ of this family of orthogonal polynomials have generating function $$\mu(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 \cite{Wall} $$\mu(x;r)= \cfrac{1}{1- \cfrac{(1-r)x^2}{1- \cfrac{x^2}{1- \cfrac{x^2}{1-\cdots}}}}.$$ We note that the Hankel transform of $\mu_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,\ldots$$ are themselves moments for the family of orthogonal polynomials that have coefficient matrix given by $$\left(\frac{(1+x)(1+rx)}{(1+x)^2}, \frac{x}{(1+x)^2}\right).$$ We have the following integral representation of the moment sequence $\mu_n(r)$. $$\mu_n(r)=\frac{-1}{\pi} \int_{2}^2 x^n \frac{\sqrt{4-x^2}(r-1)}{2(rx^2+(r-1)^2)}\,dx +\frac{r+1}{2r}\left(-\frac{r-1}{\sqrt{-r}}\right)^n + \frac{r+1}{2r}\left(\frac{r-1}{\sqrt{-r}}\right)^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{\tilde{w}'(x)}{\tilde{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. \begin{example} We consider the Riordan array $$\left(\frac{1-x-x^2}{1-3x-4x^2}, \frac{x}{1-3x-4x^2}\right)$$ which begins $$\left( \begin{array}{ccccccc} 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{array} \right),$$ with inverse $$\left(\frac{5+7x+\sqrt{1+6x+25x^2}}{2(1+x-11x^2}, \frac{\sqrt{1+6x+25x^2}-3x-1}{8x}\right)$$ that begins $$\left( \begin{array}{ccccccc} 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{array} \right).$$ The moment sequence $\mu_n$ thus begins $$1,-2,1,13,-62,97,457,\ldots.$$ The inverse matrix has a production matrix that begins $$\left( \begin{array}{ccccccc} -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{array} \right),$$ 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 $$\cfrac{1}{1+2x+ \cfrac{3x^2}{1+3x+ \cfrac{4x^2}{1+3x+ \cfrac{4x^2}{1+\cdots}}}}.$$ This implies that the moments $\mu_n$ have Hankel transform given by $$h_n = (-3)^n (-4)^{\binom{n}{2}}.$$ We find that the moments have integral representation $$\mu_n = \frac{1}{\pi} \int_{-3-4i}^{-3+4i} x^n \frac{3i\sqrt{x^2+6x+25}}{2(x^2+x-11)}\,dx+\left(\frac{3\sqrt{5}}{10}+\frac{5}{2}\right)\left(\frac{3\sqrt{5}}{2}-\frac{1}{2}\right)^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{\tilde{w}'(x)}{\tilde{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)+4 P_{n-2},$$ with $P_0(x)=1$, $P_1(x)=x+2$. \end{example} \section{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 $$\mu(x)=\frac{1}{\sqrt{1-2ax+x^2(a^2-4b)}}.$$ We now claim that the generating function $$\tilde{\mu}(x)=\frac{1}{x} \text{Rev} (x\mu(x))$$ is the generating function of the moments for another family of orthogonal polynomials, which we will call the \emph{complementary orthogonal polynomials} to the family of polynomials defined by the Riordan array $$\left(\frac{1-bx^2}{1+ax+bx^2}, \frac{x}{1+ax+bx^2}\right).$$ By solving the equation $$ \frac{u}{\sqrt{1-2au+u^2(a^2-4b)}}=x,$$ and taking the appropriate branch, we find that $$\tilde{\mu}(x)=\frac{1}{x} \text{Rev} (x\mu(x))=\frac{\sqrt{1+4bx^2}-ax}{1-x^2(a^2-4b)}=\frac{1}{ax+\sqrt{1+4bx^2}}.$$ But now $$\left(\frac{1+ax+bx^2}{1-bx^2}, \frac{x}{1-bx^2}\right)^{-1}=(\tilde{\mu}(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 $\tilde{\mu}(x)$ is the generating function of the moments of the family of orthogonal polynomials whose coefficient array is given by the Riordan array $$R=\left(\frac{1+ax+bx^2}{1-bx^2}, \frac{x}{1-bx^2}\right).$$ We regard the Riordan arrays $\tilde{R}=\left(\frac{1-bx^2}{1+ax+bx^2}, \frac{x}{1+ax+bx^2}\right)$ and $\left(\frac{1+ax+bx^2}{1-bx^2}, \frac{x}{1-bx^2}\right)$ as being complementary to each other. The production matrix of the inverse matrix $(\tilde{\mu}(x),xc(-bx^2))$ begins $$\left( \begin{array}{ccccccc} -a & 1 & 0 & 0 & 0 & 0 & 0 \\ -2 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 & 1 \\ 0 & 0 & 0 & 0 & 0 & -b & 0 \\ \end{array} \right)$$ 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)=\mathcal{P}_n(x)+a\mathcal{P}_{n-1}(x)+b\mathcal{P}_{n-2}(x)$$ where $$\mathcal{P}_n(x)=(-b)^{n/2} U\left(n, \frac{x}{2 \sqrt{-b}}\right).$$ This follows from the the factorisation of Riordan arrays given by $$\left(\frac{1+ax+bx^2}{1-bx^2}, \frac{x}{1-bx^2}\right)=(1+ax+bx^2,x)\left(\frac{1}{1-bx^2}, \frac{x}{1-bx^2}\right),$$ 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}{\pi} \frac{-1}{\sqrt{-x^2-4b}}.$$ For instance, the polynomials with coefficient array $$\left(\frac{1+2x^2}{1-2x^2}, \frac{x}{1-2x^2}\right)$$ have moments that begin $$1, 0, -4, 0, 24, 0, -160, 0, 1120, 0, -8064,0\ldots$$ given by $$\mu_n = \frac{1}{\pi}\int_{-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 \seqnum{A059304}. In the semi-classical case of $$\left(\frac{1+2x+2x^2}{1-2x^2}, \frac{x}{1-2x^2}\right),$$ we find that the moments $\mu_n$, which begin $$1, -2, 0, 8, -8, -32, 64, 128, -416, -512, 2560,\ldots,$$ have integral expression $$\mu_n = \frac{1}{\pi}\int_{-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 $\mu_n$ is given by $$h_n = 2^n (-2)^{\binom{n+1}{2}}.$$ \section{A note on the INVERT transform} For a given sequence with generating function $g(x)$, the sequence with generating function $$\frac{g(x)}{1+ a x g(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 $$\tilde{\mu}(x)=\frac{1}{ax+\sqrt{1+4bx^2}}.$$ We note that $$\tilde{\mu}(x)=\frac{\frac{1}{\sqrt{1+4bx^2}}}{1+ax\frac{1}{\sqrt{1+4bx^2}}}.$$ In other words, $\tilde{\mu}(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}} \Longrightarrow \frac{w'(x)}{w(x)}=\frac{-x}{x^2+4b},$$ so this is a classical case. Corresponding to $\tilde{\mu}(x)$ we get a weight with absolutely continuous part equal to $$\tilde{w}(x)=\frac{(x^2+a^2+4b)\sqrt{-x^2-4b^2}}{(x^2-a^2+4b)^2}.$$ We find that $\frac{\tilde{w}'(x)}{\tilde{w}(x)}$ is of classical type only if $a=0$. It is instructive to consider some more basic examples. \begin{example} We consider $w(x)=\frac{1}{\pi} \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^2x^2+4x-1}.$$ Then the absolutely continuous part of the corresponding weight function is given by $$w_a(x)=-\frac{1}{\pi} \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 (\seqnum{A081696}) $$1, 1,3,9, 29,97,333 \ldots$$ given by $$\frac{1}{\pi} \int_0^4 x^n\left(\frac{-\sqrt{x(4-x)}}{x^2-4x-1}\right)\,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,\ldots$$ given by $$\frac{1}{\pi} \int_{0}^4 x^n\left(\frac{-\sqrt{x(4-x)}}{x^2-4x-1}\right)\,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]$. \end{example} \begin{example} 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 \pi} \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=\left(\frac{1}{1+x}, \frac{x}{(1+x)^2}\right)=(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+a xc(x)}=\frac{1-\sqrt{1-4x}}{x(2+a-a\sqrt{1-4x})}=\frac{1+2ax-\sqrt{1-4x}}{2x(1+a+a^2x)}.$$ 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$, $-\infty$, 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, \ldots.$$ Thus $$C_{n+1} = \frac{1}{2\pi}\int_0^4 x^n \sqrt{x(4-x)}\,dx.$$ The coefficient array of the corresponding orthogonal polynomials is the Riordan array $$\left(\frac{1}{(1+x)^2}, \frac{x}{(1+x)^2}\right).$$ The case $a=-2$ is the sequence $\binom{2n+1}{n+1}$ with generating function $$\frac{c(x)}{1-2 xc(x)}=\frac{1-\sqrt{1-4x}}{2x\sqrt{1-4x}}.$$ We then have $$\binom{2n+1}{n+1}=\frac{1}{\pi} \int_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 $$\left(\frac{1-x}{(1+x)^2}, \frac{x}{(1+x)^2}\right).$$ For other values of $a$, we have an absolutely continuous part and an atomic measure part. $$\mu_n(a)=\frac{1}{\pi} \int_0^n x^n \frac{\sqrt{x(4-x)}}{2((1+a)x+a^2)}+\frac{(a+1)^2-1}{(a+1)^2} \left(-\frac{a^2}{a+1}\right)^n.$$ For instance, when $a=3$, the moment sequence $\mu_n(3)$ has generating function $$\mu_3(x)=\frac{1-\sqrt{1-4x}}{x(3\sqrt{1-4x}-1)},$$ and we have $$\mu_n(3)=\frac{1}{\pi} \int_0^4 x^n \frac{x(4-x)}{2(4x+9)}\,dx+\frac{15}{16}(-\frac{9}{4})^n.$$ This sequence \seqnum{A049027} begins $$ 1, 4, 17, 74, 326, 1446, 6441, 28770, 128750, 576944, \ldots.$$ The corresponding Riordan array is $$\left(\frac{1-2x}{(1+x)^2}, \frac{x}{(1+x)^2}\right).$$ In general, we have a one-parameter family of orthogonal polynomials whose coefficient arrays are given by $$\left(\frac{1-(r-1)x}{(1+x)^2}, \frac{x}{(1+x)^2}\right).$$ The moment sequence $\mu_n(r)$ then has generating function $$\mu_r(x)=\cfrac{1}{1-(r+1)x- \cfrac{x^2}{1-2x- \cfrac{x^2}{1-2x-\cdots}}}.$$ \end{example} More generally, in the context of Favard's theorem, we can assume that $$\mu(x)=\cfrac{1}{1-\alpha_0 x- \cfrac{\beta_1 x^2}{1-\alpha_1 x-\cfrac{\beta_2 x^2}{1-\cdots}}}.$$ In this case, the INVERT$(a)$ transform has generating function $$\mu_a(x)=\cfrac{1}{1-(a+\alpha_0) x- \cfrac{\beta_1 x^2}{1-\alpha_1 x-\cfrac{\beta_2 x^2}{1-\cdots}}}.$$ It is then clear that both sequences will have the same Hankel transforms. \section{A special product sequence} In this section, we investigate the sequence generated by $$G(x)=\frac{1}{x}\text{Rev}\left(\frac{x}{\tilde{\mu}(x)}\right),$$ where $$\tilde{\mu}(x)=\frac{1}{ax+\sqrt{1+4bx^2}}.$$ Thus we have $$G(x)=\frac{1}{x}\text{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}]\text{Rev}(x (ax+\sqrt{1+4bx^2}))=\frac{1}{n+1}[x^n]\left(\frac{1}{ax+\sqrt{1+4bx^2}}\right)^{n+1}.$$ Then we have \begin{eqnarray*} [x^n]G(x)&=&\frac{1}{n+1}\sum_{i=0}^{\lfloor \frac{n}{2} \rfloor}\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}.\end{eqnarray*} 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 $$\pm \sqrt{\frac{1+2ax \pm \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$ (\seqnum{A000108}) and the non-negative Fibonacci numbers $F_{n+1}$ (\seqnum{A000045}) is given by ($a=b=-1$) $$ \frac{1}{x} \sqrt{\frac{1-2x-\sqrt{1-4x-16x^2}}{10}}.$$ This is the sequence \seqnum{A098614} in the On-Line Encyclopedia of Integer Sequences \cite{SL1, SL2}. This sequence, contributed by Paul D. Hanna, begins $$1, 1, 4, 15, 70, 336, 1716, 9009, 48620,\ldots.$$ 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}$ (\seqnum{A001045}). This sequence (\seqnum{A200375}) begins $$1, 1, 6, 25, 154, 882, 5676, 36465, 244530, 1657942,\ldots$$ 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. \section{Exponential Riordan arrays and classical orthogonal polynomials} An \emph{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!} + \cdots$$ and $$f(x)=\frac{x}{1!}+f_2 \frac{x^2}{2!}+ f_3 \frac{x^3}{3!} + \cdots,$$ 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 \cite{MC}. 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 \emph{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'(\bar{v}(x)), \quad Z(x)=\frac{u'(\bar{v}(x))}{u(\bar{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} (\alpha+ \beta x + y(1+ \gamma x + \delta x^2)),$$ where we have $$Z(x)=\alpha + \beta x,\quad A(x)=1+ \gamma x + \delta x^2.$$ This leads to a production matrix that begins $$\left( \begin{array}{ccccccc} \alpha & 1 & 0 & 0 & 0 & 0 & 0 \\ \beta & \alpha +\gamma & 1 & 0 & 0 & 0 & 0 \\ 0 & 2 (\beta +\delta ) & \alpha +2 \gamma & 1 & 0 & 0 & 0 \\ 0 & 0 & 3 \beta +6 \delta & \alpha +3 \gamma & 1 & 0 & 0 \\ 0 & 0 & 0 & 4 \beta +12 \delta & \alpha +4 \gamma & 1 & 0 \\ 0 & 0 & 0 & 0 & 5 \beta +20 \delta & \alpha +5 \gamma & 1 \\ 0 & 0 & 0 & 0 & 0 & 6 \beta +30 \delta & \alpha +6 \gamma \\ \end{array} \right).$$ 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 \begin{eqnarray*}R&=&\left[e^{\int_0^x \frac{Z(t)}{A(t)}\,dt}, \int_0^x \frac{dt}{A(t)}\right]\\ &=&\left[e^{\int_0^x \frac{\alpha+\beta t}{1+\gamma t+ \delta t^2}\,dt}, \int_0^x \frac{dt}{1+\gamma t + \delta t^2}\right]\\ &=&\left[\frac{e^{\frac{\beta \gamma/\delta-2 \alpha}{4 \delta - \gamma^2} \left(2 \tan^{-1}\left(\frac{\gamma+ 2 \delta x}{\sqrt{4 \delta - \gamma^2}}\right)-2 \sin^{-1}\left(\frac{\gamma}{2 \sqrt{\delta}}\right)\right)}}{(1+ \gamma x + \delta x^2)^{\beta/(2\delta)}},\frac{1}{4 \delta - \gamma^2} \left(2 \tan^{-1}\left(\frac{\gamma+ 2 \delta x}{\sqrt{4 \delta - \gamma^2}}\right)-2 \sin^{-1}\left(\frac{\gamma}{2 \sqrt{\delta}}\right)\right)\right].\end{eqnarray*} 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-(\alpha+(n-1)\gamma))P_{n-1}(x)-(n-1)(\beta+(n-2)\delta)P_{n-2}(x),$$ with $P_0(x)=1$, $P_1(x)=x-\alpha$. \begin{example} We let $A(x)=1+x+x^2/2$, $Z(x)=1+x$. We find that $R$ is given by $$R=M^{-1}=\left[\frac{2}{2+2x+x^2}, 2 \tan^{-1}(1+x)-\frac{\pi}{2}\right].$$ The moment matrix $M$ is then given by $$M=\left[(1+\sin(x))\sec(x)^2, \tan(x)+\sin(x)-1\right].$$ The moments in this case are the shifted Euler or up/down numbers $$1, 1, 2, 5, 16, 61, 272, 1385, 7936, 50521,\ldots.$$ 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 $$\left( \begin{array}{cccccc} 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{array} \right).$$ This indicates that the Hankel transform of the moment sequence is given by $$h_n = \prod_{k=0}^n \binom{k+2}{2}^{n-k}.$$ \end{example} \begin{example} We let $A(x)=1+2x+x^2$, $Z(x)=1+x$. We find that $R$ is given by $$R=M^{-1}=\left[\frac{1}{1+x}, \frac{x}{1+x}\right]$$ with moment matrix $$M=\left[\frac{1}{1-x}, \frac{x}{1-x}\right].$$ The moments are thus $n!$ and the polynomials are the scaled Laguerre polynomials $$n! \sum_{k=0}^n \binom{n}{k}\frac{(-1)^{n-k}}{k!}x^k.$$ Classically, we have $$n! = \int_0^{\infty} 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)^2P_{n-1}(x).$$ The production matrix in this case begins $$\left( \begin{array}{cccccc} 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{array} \right).$$ This indicates that the Hankel transform of the moment sequence is given by $$h_n = \prod_{k=1}^n k^{2(n-k+1)}= \prod_{k=0}^n k!^2.$$ \end{example} \begin{example} 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=\left[ \frac{3}{2\left(\cos\left(\sqrt{3}x+\frac{\pi}{3}\right)+1\right)}, \frac{\sqrt{3}}{2} \tan\left(\frac{\sqrt{3}x}{2}+\frac{\pi}{6}\right)-\frac{1}{2}\right].$$ 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 $$\left( \begin{array}{cccccc} 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{array} \right).$$ The inverse array $R=M^{-1}$ is given by $$\left[\frac{1}{1+x+x^2}, \frac{2}{\sqrt{3}}\tan^{-1}\left(\frac{1+2x}{\sqrt{3}}-\frac{\pi}{3\sqrt{3}}\right)\right].$$ 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$. \end{example} In general, we have $$h_n = \prod_{k=1}^n (k(\beta+(k-1)\delta))^{n-k+1}.$$ \begin{proposition} The exponential Riordan array $\left[\frac{1}{(1+x)^r}, \frac{x}{1+x}\right]$ 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, \infty).$$ The moments have integral representation $$\mu_n = \int_0^{\infty} x^n e^{-x} \frac{x^{r-1}}{(r-1)!} \,dx=n!\binom{n+r-1}{n}=\prod_{k=0}^{n-1} r+k.$$ The moments have generating function given by $$\mu(x)=\cfrac{1}{1-rx-\cfrac{rx^2}{1-(r+2)x-\cfrac{2(r+1)x^2}{1-(r+4)x-\cfrac{3(r+2)x^2}{1-(r+6)x-\cdots}}}}.$$ The Hankel transform of the moments is given by $$h_n=r^n \prod_{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),\quad 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'+n y=0.$$ \end{proposition} \begin{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=\left[\frac{1}{(1+x)^r}, \frac{x}{1+x}\right]$, we have $$M^{-1}=\left[\frac{1}{(1-x)^r}, \frac{x}{1-x}\right].$$ We find that $$A(x)=(1+x)^2, \quad 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 $$\left( \begin{array}{ccccccc} r & 1 & 0 & 0 & 0 & 0 & 0 \\ r & r+2 & 1 & 0 & 0 & 0 & 0 \\ 0 & 2 r+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{array} \right).$$ The continued fraction, the Hankel transform and three-term recurrence now follow. \end{proof} \begin{proposition} The exponential Riordan array $\left[\frac{1}{(1+x)^{r+1}}, \frac{x}{1+x}\right]$ 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, \infty).$$ The moments have integral representation $$\mu_n = \int_0^{\infty} x^n e^{-x} \frac{x^r}{r!} \,dx=n!\binom{n+r}{r}=\prod_{k=1}^n r+k.$$ The moments have generating function given by $$\mu(x)=\cfrac{1}{1-(r+1)x-\cfrac{(r+1)x^2}{1-(r+3)x-\cfrac{2(r+2)x^2}{1-(r+5)x-\cfrac{3(r+3)x^2}{1-(r+7)x-\cdots}}}}.$$ The Hankel transform of the moments is given by $$h_n=\prod_{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),\quad 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'+n y=0.$$ \end{proposition} \begin{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=\left[\frac{1}{(1+x)^{r+1}}, \frac{x}{1+x}\right]$ we have $M^{-1}=\left[\frac{1}{(1-x)^{r+1}}, \frac{x}{1-x}\right]$. 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 $$\left( \begin{array}{ccccccc} r+1 & 1 & 0 & 0 & 0 & 0 & 0 \\ r+1 & r+3 & 1 & 0 & 0 & 0 & 0 \\ 0 & 2 r+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{array} \right).$$ The continued fraction, the Hankel transform and three-term recurrence now follow. \end{proof} \begin{proposition} The exponential Riordan array $\left[e^{-\frac{rx^2}{2}},x\right]$ is the coefficient array of a family of classical orthogonal polynomials. We have $$w(x)=e^{-\frac{x^2}{2r}}$$ on the interval $$(-\infty, \infty).$$ The moments have integral representation $$\mu_n = \int_{-\infty}^{\infty} x^n e^{-\frac{x^2}{2r}} \,dx.$$ These begin $$1,0,r,0,3r^2,0,15r^3,0,105 r^4,0,945r^5,0,\ldots.$$ The moments have generating function given by $$\mu(x)=\cfrac{1}{1-\cfrac{rx^2}{1-\cfrac{2rx^2}{1-\cfrac{3rx^2}{1-\cdots}}}}.$$ The Hankel transform of the moments is given by $$h_n=r^{\binom{n+1}{2}}\prod_{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),\quad n>1,$$ with $P_0(x)=1$, $P_1(x)=x$. If $y=P_n(x)$, then $y$ satisfies the differential equation $$ry''-xy'+n y=0.$$ \end{proposition} \begin{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]=\left[e^{\frac{rx^2}{2}},x\right]$, and hence the production matrix is generated by $e^{xy}(rx+y)$. This matrix begins $$\left( \begin{array}{ccccccccc} 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ r & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 2 r & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 3 r & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 4 r & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 5 r & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 6 r & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 7 r & 0 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 8 r & 0 \\ \end{array} \right).$$ The continued fraction, the Hankel transform and three-term recurrence now follow. \end{proof} The case $r=1$ is related to the Hermite polynomials. The Riordan array $\left[e^{-\frac{t^2}{2}}, t\right]$ is the coefficient array of the (probabilist) Hermite polynomials $He(n,x)$ given by $$He_n(x)=\sum_{k=0}^n \frac{n!}{(-2)^{\frac{n-k}{2}} k! \left(\frac{n-k}{2}\right)!} \frac{1+(-1)^{n-k}}{2} x^k.$$ The Physicists' Hermite polynomials are given by $$H_n(x)=n! \sum_{k=0}^{\lfloor \frac{n}{2} \rfloor}\frac{(-1)^k (2x)^{n-2k}}{k!(n-2k)!}.$$ The coefficient array of the Physicists' Hermite polynomials is the exponential Riordan array $$\left[e^{-t^2}, 2t\right].$$ We have $$He_n(x)=2^{-\frac{n}{2}}H_n(\sqrt{2}x).$$ We have $$He_n(x)=x He_{n-1}(x)-(n-1)He_{n-2}(x),$$ and $$H_n(x) = 2x H_{n-1}(x)-2(n-1) H_{n-1}(x).$$ \section{Bessel and related polynomials} In this section we briefly look at two families of polynomials related to the Bessel polynomials \cite{Grosswald}. We start with the exponential Riordan array $$R=\left[\frac{1}{\sqrt{1+2x}}, \frac{x}{1+2x}\right].$$ This array begins $$\left( \begin{array}{ccccccc} 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{array} \right).$$ The production matrix of the inverse $M=R^{-1}=\left[\frac{1}{\sqrt{1-2x}}, \frac{x}{1-2x}\right]$ of this array is tri-diagonal, beginning $$\left( \begin{array}{ccccccc} 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{array} \right).$$ It follows that the exponential Riordan array $\left[\frac{1}{\sqrt{1+2x}}, \frac{x}{1+2x}\right]$ 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 $\mu_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 (\seqnum{A001147}) $$1, 1, 3, 15, 105, 945, 10395, 135135, 2027025, 34459425, \ldots.$$ We have $$(2n-1)!! = \frac{n!}{2 \pi} \int_0^2 \frac{x^n}{\sqrt{x(2-x)}}\,dx=\frac{1}{\sqrt{2 \pi}}\int_0^{\infty} x^n \frac{e^{-\frac{x}{2}}}{\sqrt{x}}\,dx.$$ The Hankel transform of the moments $\mu_n$ is then given by $$h_n=\prod_{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 \pi}} \frac{e^{-\frac{x}{2}}}{\sqrt{x}}.$$ The inverse matrix $M=\left[\frac{1}{\sqrt{1-2x}}, \frac{x}{1-2x}\right]$ 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 $$\tilde{M}=\left[\frac{d}{dx}\frac{1}{\sqrt{1-2x}}, \frac{x}{1-2x}\right]=\left[\frac{1}{(1-2x)^{3/2}}, \frac{x}{1-2x}\right]$$ has $$Z(x)=3+6x, \quad A(x)=(1+2x)^2.$$ Its production matrix thus begins $$\left( \begin{array}{ccccccc} 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{array} \right).$$ The array $\tilde{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 $$\tilde{\mu}_n = \frac{1}{\sqrt{2 \pi}}\int_0^{\infty} x^n x\frac{e^{-\frac{x}{2}}}{\sqrt{x}}\,dx.$$ Thus $$\tilde{w}(x)=\frac{1}{\sqrt{2 \pi}}\sqrt{x} e^{-\frac{x}{2}},$$ and $$\frac{\tilde{w}'(x)}{\tilde{w}(x)}=\frac{1-x}{2x}.$$ The Hankel transform of the moments $\tilde{\mu}_n$, which begin $$1, 3, 15, 105, 945, 10395, 135135, 2027025, 34459425, \ldots,$$ is given by $$h_n=\prod_{k=0}^n (2(k+1)(2k+3))^{n-k}.$$ The exponential Riordan array $\left[e^{\frac{x^2}{2}},x\right]$ 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,\ldots$$ with exponential generating function $e^{x^2/2}$. The Hankel transform of these moments is given by $$h_n = \prod_{k=0} (k+1)^{n-k}=\prod_{k=0}^n k!$$ This follows since for the moment matrix $\left[e^{x^2/2}, x\right]$ we have $A(x)=1$ and $Z(x)=x$. We now consider a family of Bessel polynomials. The exponential Riordan array $B=\left[\frac{1}{\sqrt{1-2x}},1-\sqrt{1-2x}\right]$ is the coefficient array of the \emph{reverse Bessel polynomials} $\theta_n(x)$ \cite{Grosswald} with general term $$\theta_n(x)=\sum_{k=0}^n \frac{(2n-k)!}{2^{n-k}k!(n-k)!} x^k.$$ These polynomials are orthogonal on the circle, satisfying the recurrence $$\theta_n(x)=(2n-1) \theta_{n-1}(x)+x^2 \theta_{n-2}(x).$$ We calculate $RB$. We get \begin{eqnarray*} RB&=&\left[\frac{1}{\sqrt{1+2x}}, \frac{x}{1+2x}\right]\cdot \left[\frac{1}{\sqrt{1-2x}},1-\sqrt{1-2x}\right]\\ &=&\left [\frac{1}{\sqrt{1+2x}} \frac{1}{\sqrt{1-2\frac{x}{1+2x}}},1-\sqrt{1-2\frac{x}{1+2x}}\right]\\ &=&\left[1, 1-\frac{1}{\sqrt{1+2x}}\right].\end{eqnarray*} This last array has general term $$T_{n,k}=\frac{(n-1)!}{(k-1)!} \sum_{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} \cdot \left[1, 1-\frac{1}{\sqrt{1+2x}}\right].$$ 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 $\theta_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 $\left[e^{-\frac{x^2}{2}}, x\right]$, then $T_{2n-k,k}$ is the general term of the coefficient array $\left[\frac{1}{\sqrt{1+2x}}, \sqrt{1+2x}-1\right]$. More generally, if $T_{n,k}$ is the general term of the coefficient array $\left[e^{-\frac{rx^2}{2}}, x\right]$, then $T_{2n-k,k}$ is the general term of the coefficient array $\left[\frac{1}{\sqrt{1+2xr}}, \sqrt{1+2xr}-1\right]$. The matrix $$R=\left[\frac{1}{\sqrt{1+2x}}, \frac{x}{1+2x}\right]$$ is a member of a one-parameter family of coefficient arrays of orthogonal polynomials. The general element of this family is $$\left[\frac{1}{(1+rx)^{1/r}}, \frac{x}{1+rx}\right]=\left[\frac{1}{(1-rx)^{1/r}}, \frac{x}{1-rx}\right]^{-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),\ldots,$$ or $$\mu_n(r)=\prod_{k=0}^{n-1} kr+1.$$ The Hankel transform of $\mu_n$ is then $$h_n=r^{\binom{n+1}{2}} \prod_{k=0}^n k!(kr+1)^{n-k}.$$ The corresponding family of generalized Bessel polynomials will then have coefficient array given by $$ \left[\frac{1}{(1-rx)^{1/r}}, 1-(1-rx)^{1/r}\right] = \left[1-x, \frac{1-(1-x)^r}{r}\right]^{-1}.$$ \section{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 \cite{Toda}. 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. \section{Appendix --- The Stieltjes transform of a measure} The \emph{Stieltjes transform} of a measure $\mu$ on $\mathbb{R}$ is a function $G_{\mu}$ defined on $\mathbb{C}\setminus \mathbb{R}$ by $$G_{\mu}(z)=\int_{\mathbb{R}}\frac{1}{z-t}\mu(t).$$ If $f$ is a bounded continuous function on $\mathbb{R}$, we have $$\int_{\mathbb{R}}f(x)\mu(x)=-\lim_{y \to 0^{+}}\int_{\mathbb{R}}f(x)\Im G_{\mu}(x+iy)dx.$$ If $\mu$ has compact support, then $G_{\mu}$ is holomorphic at infinity and for large $z$, $$G_{\mu}(z)=\sum_{n=0}^{\infty} \frac{a_n}{z^{n+1}},$$ where $a_n=\int_{\mathbb{R}}t^n \mu(t)$ are the moments of the measure. If $\mu(t)=d\psi(t)=\psi'(t)dt$ then (Stieltjes-Perron) $$\psi(t)-\psi(t_0)=-\frac{1}{\pi}\lim_{y \to 0^{+}}\int_{t_0}^t \Im G_{\mu}(x+iy)dx.$$ If now $g(x)$ is the generating function of a sequence $a_n$, with $g(x)=\sum_{n=0}^{\infty}a_n x^n$, then we can define $$G(z)=\frac{1}{z}g\left(\frac{1}{z}\right)=\sum_{n=0}^{\infty}\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. \begin{example} 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\left(\frac{1}{z}\right)=\frac{1}{2}\left(1-\sqrt{\frac{x-4}{x}}\right).$$ Then $$\Im G_{\mu}(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 \begin{eqnarray*}\psi'(x)&=&-\frac{1}{\pi}\lim_{y \to 0^{+}}\left\{-\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}}\right\}\\ &=& \frac{1}{2\pi}\frac{\sqrt{x(4-x)}}{x}.\end{eqnarray*} \end{example} \begin{thebibliography}{99} \bibitem{Agida} M. Agida and A.~S. Kumar, A Boubaker polynomials expansion scheme solution to random Love's equation in the case of a rational kernel, \emph{Elec. J. Theoretical Phys.}, \textbf{7} (2010), 319--326. \bibitem{NMR} 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, \emph{Current Applied Physics}, \textbf{9} (2009), 278--283. \bibitem{CB} P. Barry, On the restricted Chebyshev-Boubaker polynomials, \emph{Integral Transforms Spec. Funct.}, \textbf{28} (2017), 1--16. \bibitem{Book} P. Barry, \emph{Riordan Arrays: a Primer}, Logic Press, 2017. \bibitem{Barry_moments} P. Barry, Riordan arrays, orthogonal polynomials as moments, and Hankel transforms, \emph{J. Integer Sequences}, \textbf{14} (2011), \href{https://cs.uwaterloo.ca/journals/JIS/VOL14/Barry1/barry97r2.pdf} {Article 11.2.2}. \bibitem{Meixner} P. Barry and A. Hennessy, Meixner-type results for Riordan arrays and associated integer sequences, \emph{J. Integer Sequences}, \textbf{13} (2010), \href{https://cs.uwaterloo.ca/journals/JIS/VOL13/Barry5/barry96s.pdf} {Article 10.9.4}. \bibitem{Toda} P. Barry, The restricted Toda chain, exponential Riordan arrays, and Hankel transforms, \emph{J. Integer Sequences}, \textbf{13} (2010), \href{https://cs.uwaterloo.ca/journals/JIS/VOL13/Barry3/barry100r.pdf} {Article 10.8.4}. \bibitem{Boubaker} K. Boubaker and L. Zhang, Fermat-linked relations for the Boubaker polynomial sequences via Riordan matrices analysis, \emph{J. Assoc. Arab Univ. Basic Appl. Sci.}, \textbf{12} (2012), 74--78. \bibitem{BPES} K. Boubaker, Analytical initial-guess-free solution to Kepler's transcendental equation using Boubaker polynomials expansion scheme BPES \emph{Apeiron: Studies in Infinite Nature}, \textbf{17} (2010), 1--12. \bibitem{Boubaker07} K. Boubaker, A. Chaouachi, M. Amlouk, and H. Bouzouita, Enhancement of pyrolysis spray dispersal performance using thermal time-response to precursor uniform deposition, \emph{Eur. Phys. J. AP}, \textbf{37} (2007), 105--109. \bibitem{Chihara} T. S. Chihara, {\it An Introduction to Orthogonal Polynomials}, Dover Publicatons, 2011. \bibitem{Dada} O.~M. Dada, O.~B. Awojoyogbe, M. Agida, and K. Boubaker, Variable separation and Boubaker polynomial expansion scheme for solving the neutron transport equation, \emph{Physics International}, \textbf{2} (2011), 25--30. \bibitem{ProdMat_0} E. Deutsch, L. Ferrari, and S. Rinaldi, Production matrices, \emph{Adv. in Appl. Math.}, \textbf{34} (2005), 101--122. \bibitem{ProdMat} E. Deutsch, L. Ferrari, and S. Rinaldi, Production matrices and Riordan arrays, \emph{Ann. Comb.}, \textbf{13} (2009), 65--85. \bibitem{Gautschi} W. Gautschi, {\it Orthogonal Polynomials: Computation and Approximation}, Clarendon Press, 2004. \bibitem{Grosswald} E. Grosswald, \emph{Bessel Polynomials}, Springer, 1978. \bibitem{Henrici} P. Henrici, \emph{Applied and Computational Complex Analysis, Vol.~3}, John Wiley \& Sons, 1993. \bibitem{Kratt} C. Krattenthaler, Advanced determinant calculus, in D. Foata, and G.-N. Han, eds., \emph{The Andrews Festschrift: Seventeen Papers on Classical Number Theory and Combinatorics}, Springer, 2001, pp.~349--426. \bibitem{Labadiah1} 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, \emph{J. of Diff. Eq. and C. Proc.}, (2007) 117-133. \bibitem{Labadiah2} 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, \emph{J. of Diff. Eq. and C. Proc.}, (2008), 51--66. \bibitem{Layman} J. W. Layman, The Hankel transform and some of its properties, \emph{J. Integer Sequences}, \textbf{4} (2001), \href{https://www.cs.uwaterloo.ca/journals/JIS/VOL4/LAYMAN/hankel.html} {Article 01.1.5}. \bibitem{Chebyshev} J. C. Mason and D. C. Handscomb, \emph{Chebyshev Polynomials}, Chapman and Hall/CRC, 2002. \bibitem{Merlini_When} D. Merlini, R. Sprugnoli, and M.~C.~Verri, Lagrange inversion: when and how, \emph{Acta Appl. Math.}, \textbf{94} (2006), 233--249. \bibitem{MC} D. Merlini, R. Sprugnoli, and M.~C. Verri, The method of coefficients, \emph{Amer. Math. Monthly}, \textbf{114} (2007), 40--57. \bibitem{Milanovic} G. V. Milanovic and D. Joksimovic, Properties of Boubaker polynomials and an application to Love's integral equation, \emph{Appl. Math. Comput.}, \textbf{224} (2013), 74--87. \bibitem{Noe} T. D Noe, On the divisibility of generalized central trinomial coefficients, \emph{J. Integer Sequences}, \textbf{9} (2006), \href{https://cs.uwaterloo.ca/journals/JIS/VOL9/Noe/noe35.html} {Article 06.2.7}. \bibitem{P_W} P. Peart and W.-J. Woan, Generating functions via Hankel and Stieltjes matrices, \emph{J. Integer Sequences}, \textbf{3} (2000), \href{http://www.cs.uwaterloo.ca/journals/JIS/VOL3/PEART/peart1.html} {Article 00.2.1}. \bibitem{SGWW} L.~W.~Shapiro, S. Getu, W.-J. Woan, and L.~C. Woodson, The Riordan group, \emph{Discr. Appl. Math.} \textbf{34} (1991), 229--239. \bibitem{Survey} L. Shapiro, A survey of the Riordan group, available electronically at \href{https://tinyurl.com/ycc8vy34}, Center for Combinatorics, Nankai University, 2005. \bibitem{Shapiro_bij} L.~W.~Shapiro, Bijections and the Riordan group, \emph{Theoret. Comput. Sci.} \textbf{307} (2003), 403--413. \bibitem{Slama} 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, \emph{Eur. Phys. J. AP}, \textbf{44} (2008), 317--322. \bibitem{SL1} N. J. A.~Sloane, \emph{The On-Line Encyclopedia of Integer Sequences}. Published electronically at \url{http://oeis.org}, 2016. \bibitem{SL2} N. J. A.~Sloane, The on-line encyclopedia of integer sequences, \emph{Notices Amer. Math. Soc.}, \textbf{50} (2003), 912--915. \bibitem{Spru} R. Sprugnoli, Riordan arrays and combinatorial sums, \emph{Discrete Math.} \textbf{132} (1994), 267--290. \bibitem{Suetin} P. K. Suetin, Classical orthogonal polynomials, entry in the on-line \emph{Encyclopedia of Mathematics}, 2016. Available at \\ \url{https://www.encyclopediaofmath.org/index.php/Classical_orthogonal_polynomials}. \bibitem{Szego} G. Szeg\"o, \emph{Orthogonal Polynomials}, 4th edition, Amer. Math. Soc., 1975. \bibitem{Viennot} G. Viennot, Une th\'eorie combinatoire des polyn\^omes orthogonaux g\'en\'eraux, UQAM, Montreal, Quebec, 1983. \bibitem{Wall} H.~S. Wall, \emph{Analytic Theory of Continued Fractions}, AMS Chelsea Publishing, 2000. \bibitem{Wigner_1} E. P. Wigner, Characteristic vectors of bordered matrices with infinite dimensions, \emph{Ann. of Math.}, \textbf{62} (1955), 548--564. \bibitem{Wigner_2} E. P. Wigner, On the distribution of the roots of certain symmetric matrices, \emph{Ann. of Math.}, \textbf{67} (1958), 325--327. \bibitem{Zhao} T. G. Zhao, L. Naing, and W. X. Yue, Some new features of the Boubaker polynomials expansion scheme BPES, \emph{J. Assoc. Arab Univ. Basic Appl. Sci.}, \textbf{12} (2012) 74--78. \end{thebibliography} \bigskip \hrule \bigskip \noindent 2010 {\it Mathematics Subject Classification}: Primary 15B36; Secondary 33C45, 11B83, 11C20, 05A15. \noindent \emph{Keywords:} Riordan array, classical orthogonal polynomial, Chebyshev polynomial, Bessel polynomial, moment sequence. \bigskip \hrule \bigskip \noindent (Concerned with sequences \seqnum{A000045}, \seqnum{A000108}, \seqnum{A000984}, \seqnum{A001045}, \seqnum{A001147}, \seqnum{A049027}, \seqnum{A059304}, \seqnum{A081696}, \seqnum{A098614}, and \seqnum{A200375}.) \bigskip \hrule \bigskip \vspace*{+.1in} \noindent Received January 18 2017; revised versions received December 29 2017; December 30 2017. Published in {\it Journal of Integer Sequences}, January 21 2018. \bigskip \hrule \bigskip \noindent Return to \htmladdnormallink{Journal of Integer Sequences home page}{http://www.cs.uwaterloo.ca/journals/JIS/}. \vskip .1in \end{document} .