# Orthogonal polynomials on a class of planar algebraic curves

Marco Fasondini  School of Computing and Mathematical Sciences
University of Leicester
United Kingdom<EMAIL_ADDRESS>,  Sheehan Olver  Department of
Mathematics
Imperial College
London
United Kingdom<EMAIL_ADDRESS>and  Yuan Xu  Department of Mathematics
University of Oregon
Eugene, Oregon 97403-1222<EMAIL_ADDRESS>
###### Abstract.

We construct bivariate orthogonal polynomials (OPs) on algebraic curves of the
form $y^{m}=\phi(x)$ in ${\mathbb{R}}^{2}$ where $m=1,2$ and $\phi$ is a
polynomial of arbitrary degree $d$, in terms of univariate semiclassical OPs.
We compute connection coefficients that relate the bivariate OPs to a
polynomial basis that is itself orthogonal and whose span contains the OPs as
a subspace. The connection matrix is shown to be banded and the connection
coefficients and Jacobi matrices for OPs of degree $0,\ldots,N$ are computed
via the Lanczos algorithm in $\mathcal{O}(Nd^{4})$ operations.

##  1\. Introduction

In previous works, we studied orthogonal polynomials (OPs) on planar curves
and on quadratic surfaces of revolution. The OPs were constructed explicitly
on the wedge and the square [11], planar quadratic curves [14], cubic curves
[6], as well as on and inside quadratic surfaces of revolution [12]. These OPs
were used in applications, such as the statistics of determinantal point
processes, the computation of Stieltjes transforms, multivariate function
approximation and orthogonal series [21, 22, 23, 24, 25, 26], the
approximation of singular or nearly singular functions of one variable
(including those that arise as solutions to differential equations), the
explicit solution to wave equations [15], spectral methods for partial
differential equations on trapeziums, disk slices and spherical caps [17, 18]
and a spectral method for fractional integral equations [16].

In this paper we generalise both the types of curves on which OPs are
constructed as well as their method of construction, in particular to OPs
defined on curves of the form $y^{m}=\phi(x)$ in ${\mathbb{R}}^{2}$ where
$m=1,2$ and $\phi$ is a polynomial of arbitrary degree $d$. Unlike the
previously mentioned studies, our constructions of OPs are computational, not
explicit. However, on quadratic and cubic curves, our computational method
reproduces the OP bases given in [14, 6] and in section 4 we give explicit
bases for a few new special cases. Our approach is first to explicitly
construct a polynomial basis whose span contains the bivariate OPs as a
subspace, then we cast the computation of the OPs as a numerical linear
algebra problem, namely, simultaneously block-tridiagonalize two symmetric,
commuting matrices, which we solve with the Lanczos algorithm. This is a very
general approach that allows us to construct OPs (i) on certain algebraic
curves (the topic of this paper) and also (ii) inside those curves (two-
dimensional domains), (iii) their surfaces of revolution (three-dimensional
regions) and (iv) associated star domains (2D and 3D regions). The
construction of OPs on (ii)–(iv) will be considered elsewhere. Sparse spectral
methods for PDEs is an application of the OPs we construct that is of special
interest and will be pursued in future research.

Acknowledgment. The first and second authors were supported by the Leverhulme
Trust Research Project Grant RPG-2019-144 “Constructive approximation theory
on and inside algebraic curves and surfaces”. The second author was also
supported by the Engineering and Physical Sciences Research Council (EPSRC)
Standard Grant EP/T022132/1. The third author was partially supported by
Simons Foundation Grant #849676.

##  2\. Univariate OPs

Throughout, we shall denote a univariate, orthonormal, degree-$n$ OP that is
orthogonal with respect to the inner product

(2.1) ${\langle}f,g{\rangle}_{\omega}=\int f(x)g(x)\omega(x)\mathrm{d}x,$

where $\omega$ is a non-negative weight function, by $p_{n}(\omega;x)$, or
$p_{n}(\omega)$. We shall also make use of quasimatrix [19] notation. For
example, the quasimatrix

(2.3) $\displaystyle\mathbf{P}(\phi w)=\left(\begin{array}[]{c c c
c}p_{0}(\phi w)&p_{1}(\phi w)&p_{2}(\phi w)&\cdots\end{array}\right),$

where $\phi=\phi(x)$ is polynomial that is strictly positive on the support of
the weight $w$, is a ‘matrix’ with a countable infinity of columns consisting
of OPs with respect to the weight $\phi w$, ordered by degrees.

We shall let $w$ be a classical weight, i.e., a Jacobi, Laguerre or Hermite
weight [4, Ch. 18] with support $\mathrm{supp}(w)=(-1,1),(0,\infty)$ and
$(-\infty,\infty)$, respectively. If all the roots of the polynomial $\phi$
(with $\phi(x)>0$ for $x\in\mathrm{supp}(w)$) are at the endpoints of
$\mathrm{supp}(w)$, then $\phi w$ is also a classical weight, otherwise at
least one of the roots of $\phi$ are off $\mathrm{supp}(w)$ and $\phi w$ is a
semiclassical weight [10] and we refer to the $p_{n}(\phi w)$ as (orthonormal)
semiclassical OPs.

The fundamental building blocks of our constructions of bivariate OPs will be
Jacobi matrices and a ‘raising matrix’ of univariate OPs, which we shall now
define. Recall that if a family of OPs, e.g., $p_{n}(\phi w)$, $n\geq 0$, is
orthonormal with respect to an inner product
${\langle}\cdot,\cdot{\rangle}_{\phi w}$, then they satisfy a three-term
recurrence of the form

(2.4) $\displaystyle xp_{n}(\phi w)=\beta_{n-1}p_{n-1}(\phi
w)+\alpha_{n}p_{n}(\phi w)+\beta_{n}p_{n+1}(\phi w),\qquad n\geq 0,$

where $\beta_{-1}=0=p_{-1}(\phi w)$ and the recurrence coefficients are

(2.5) $\alpha_{n}=\alpha_{n}(\phi w)={\langle}xp_{n}(\phi w),p_{n}(\phi
w){\rangle}_{\phi w},\qquad\beta_{n}=\beta_{n}(\phi w)={\langle}xp_{n}(\phi
w),p_{n+1}(\phi w){\rangle}_{\phi w}.$

Since $w$ is a classical weight and $\phi(x)>0$ for $x\in\mathrm{supp}(w)$,
$\alpha_{n}(\phi w),\beta_{n}(\phi w)\in\mathbb{R}$ for $n\geq 0$. In
quasimatrix notation, (2.4) becomes

(2.6) $x\mathbf{P}(\phi w)=\mathbf{P}(\phi w)J(\phi w),$

where $J(\phi w)$ is an infinite, symmetric, tridiagonal matrix known as the
Jacobi matrix of the OPs,

(2.7) $J(\phi w)=\left(\begin{array}[]{c c c c}\alpha_{0}(\phi
w)&\beta_{0}(\phi w)&&\\\ \beta_{0}(\phi w)&\alpha_{1}(\phi w)&\beta_{1}(\phi
w)&\\\ &\beta_{1}(\phi w)&\alpha_{2}(\phi w)&\ddots\\\
&&\ddots&\ddots\end{array}\right).$

We say that a (finite or infinite) matrix $A$ is banded and has bandwidths
$(\lambda,\mu)$, $\lambda,\mu\geq 0$ if $A_{i,j}=0$ for $i-j>\lambda$ and
$j-i>\mu$ and the same applies for block-matrices (in which case $A_{i,j}$
refers to block $(i,j)$ of $A$). Hence, $J(\phi w)$ has bandwidths $(1,1)$.

######  Proposition 2.1.

If $\deg\phi=d$, then the OPs $\\{p_{n}(w)\\}$ and $\\{p_{n}(\phi w)\\}$ are
related via a banded raising matrix $R$ with bandwidths $(0,d)$ according to

(2.8) $\mathbf{P}(w)=\mathbf{P}(\phi w)R,$

where $\mathbf{P}(\phi w)$ is given in (2.3) and $\mathbf{P}(w)$ is defined
analogously and

(2.9) $R=\left(\begin{array}[]{c c c c c c
c}r_{0,0}&\cdots&\cdots&r_{0,d}&&&\\\ &r_{1,1}&\cdots&\cdots&r_{1,d+1}&&\\\
&&r_{2,2}&\cdots&\cdots&r_{2,d+2}&\\\
&&&\ddots&\ddots&\ddots&\ddots\end{array}\right),$

where

(2.10) $r_{k,n}=\langle p_{n}(w),p_{k}(\phi w)\rangle_{\phi w}=\langle
p_{n}(w),\phi p_{k}(\phi w)\rangle_{w}.$

###### Proof.

Since any polynomial of degree $n$ is a linear combination of the $p_{k}(\phi
w)$ with $0\leq k\leq n$, we have $p_{n}(w)=\sum_{k=0}^{n}r_{k,n}p_{k}(\phi
w)$, where the expression (2.10) and the fact that $r_{k,n}=0$ for $0\leq
k<n-d$ follow since the $p_{k}(\phi w)$ and $p_{n}(w)$ are (orthonormal) OPs.
∎

In quasimatrix notation, we write

(2.11) $R=\langle\mathbf{P}(\phi w),\mathbf{P}(w)\rangle_{\phi
w}=\int\mathbf{P}^{\top}(\phi w)\mathbf{P}(w)\phi w\,\mathrm{d}x,$

where $\mathbf{P}^{\top}(\phi w)\mathbf{P}(w)$ is interpreted as an outer
product of two (infinite) row vectors whose entries are functions and
integration is performed entry-wise. Therefore, $\langle\mathbf{P}(\phi
w),\mathbf{P}(w)\rangle_{\phi w}$ is an infinite matrix whose entries are
defined as follows: Let $\bm{e}_{k}$ be the standard basis vector

(2.12) $\bm{e}_{k}=\Big{(}\underbrace{0\>\cdots\>0}_{k\text{
zeros}}\>1\>0\>0\>\cdots\Big{)}^{\top}$

then the $(k,n)$ entry of $\langle\mathbf{P}(\phi
w),\mathbf{P}(w)\rangle_{\phi w}$ for $k,n\geq 0$ is (cf. (2.10))

(2.13) $\left(\langle\mathbf{P}(\phi w),\mathbf{P}(w)\rangle_{\phi
w}\right)_{k,n}=\int\mathbf{P}(\phi w)\bm{e}_{k}\mathbf{P}(w)\bm{e}_{n}\phi
w\,dx=\langle p_{n}(w),p_{k}(\phi w)\rangle_{\phi w}=r_{k,n}.$

Similarly, we can express the Jacobi matrix as $J(\phi
w)=\langle\mathbf{P}(\phi w),x\mathbf{P}(\phi w)\rangle_{\phi w}$ and since
the OPs $\\{p_{n}(w)\\}$ and $\\{p_{n}(\phi w)\\}$ are orthonormal,
$\langle\mathbf{P}(\phi w),\mathbf{P}(\phi w)\rangle_{\phi
w}=I=\langle\mathbf{P}(w),\mathbf{P}(w)\rangle_{w}$, where $I$ is an infinite
identity matrix.

The coefficients $r_{k,n}$ of the raising matrix $R$ can also be viewed as the
connection coefficients relating the classical OP basis $\\{p_{n}(w)\\}$ to
the semiclassical OP basis $\\{p_{n}(\phi w)\\}$. Algorithms for computing
connection coefficients between (mainly univariate) classical OP families with
different weights are discussed in [13]. The more general problem of relating
the OP families $\\{p_{n}(rw)\\}$ and $\\{p_{n}(w)\\}$, where $r=r(x)$ is a
rational function, and determining the recurrence coefficients of the former
OP family is discussed in [7, sect. 2.4] and associated algorithms are
considered in [9, Ch. 5]. For our purposes (i.e., the case $r(x)=\phi(x)$), we
note that one method for computing $R$ and $J(\phi w)$ is via the Cholesky
factor $L$ of the symmetric, positive definite matrix

$\Phi=\langle\mathbf{P}(w),\mathbf{P}(w)\rangle_{\phi
w}=\int\mathbf{P}(w)^{\top}\mathbf{P}(w)\phi w\,\mathrm{d}x=LL^{\top}$

because from (2.8) we have that

$I=\langle\mathbf{P}(\phi w),\mathbf{P}(\phi w)\rangle_{w}=\int\mathbf{P}(\phi
w)^{\top}\mathbf{P}(\phi w)\phi w\,\mathrm{d}x=R^{-\top}\Phi
R^{-1}=R^{-\top}LL^{\top}R^{-1}$

from which we conclude that

(2.14) $R=L^{\top},\qquad J(\phi w)=R^{\top}J(w)R,$

where the second identity follows from (2.6) and (2.8). Another method for
computing $R$ and $J(\phi w)$ given $J(w)$ is to use Christoffel’s theorem [7,
sect. 2.4].

The topic of this paper is also the computation of connection coefficients
(and, as a byproduct, recurrence coefficients will also be computed), however
we shall relate bivariate OPs (restricted to curves) to a polynomial basis
that is orthogonal with respect to the _same_ weight as the bivariate OPs. As
we shall see, not only is our connection problem complicated by the fact that
it is bivariate but also because in general the degree of a bivariate
polynomial on $\mathbb{R}^{2}$ is different from its degree when restricted to
an algebraic curve in $\mathbb{R}^{2}$.

##  3\. Bivariate OPs on $y=\phi$ and $y^{2}=\phi$

Curves defined by $p(x,y)=0$, where $p$ is an irreducible bivariate quadratic
or cubic polynomial, can be simplified to the form $y^{2}=\phi(x)$, with
$\deg\phi\leq 2$ for quadratics and $\deg\phi=3$ for cubics, see [14] and [6],
respectively. Through affine transformations, non-degenerate quadratics can be
reduced to either a parabola ($\phi=x$), a circle ($\phi=1-x^{2}$), or a
hyperbola ($\phi=x^{2}-1$). A cubic curve in the standard form
$y^{2}=\phi(x)$, $\deg\phi=3$ is an elliptic curve if it has no cusps,
isolated points or self-intersections; the OPs we constructed in [6] are valid
on both elliptic and non-elliptic curves. In this paper we consider bivariate
OPs on the following more general curves

(3.1)
$\displaystyle{\gamma}_{m,d}=\left\\{(x,y)\in{\mathbb{R}}^{2}:\begin{cases}y=\phi(x)&\text{if
}m=1\\\ y^{2}=\phi(x)>0&\text{if
}m=2\end{cases},x\in\mathrm{supp}(w),\deg\phi=d>0\right\\}$

where $\mathrm{supp}(w)$ denotes the support of a non-negative weight function
$w(x)$; either $\mathrm{supp}(w)=(a,b)$, $(a,\infty)$ or $(-\infty,\infty)$,
where $-\infty<a<b<\infty$ and $\phi$ may vanish at $a$ or $b$. The
restriction $\phi>0$ does not apply for $m=1$ but for $m=2$, the positivity of
$\phi$ is necessary for the curve to be defined in $\mathbb{R}^{2}$.

For the special case $m=1$, $d=0$ and a classical weight $w$, the classical
OPs are OPs on the curves ${\gamma}_{1,0}$. The case $m=2$ and $d=0$ is
treated in [14]. OPs on curves of the form $y^{m}=\phi$ with $m>2$, which we
leave for future work, will be discussed briefly in section 6.

We consider OPs in $(x,y)$ that are restricted to the curve ${\gamma}_{m,d}$
and which are orthogonal with respect to the following inner product

${\langle}f,g{\rangle}_{{\gamma}_{m,d},w}=\int_{{\gamma}_{m,d}}f(x,y)g(x,y)\omega(x,y)\mathrm{d}\sigma(x,y),$

where $\mathrm{d}\sigma$ is the arc length measure on ${\gamma}_{m,d}$. We
incorporate the arc length into the non-negative weight function by setting
$\omega(x,y)=w(x)/\sqrt{1+(y^{\prime}(x))^{2}}$ so that the inner product
simplifies to

(3.2)
$\displaystyle{\langle}f,g{\rangle}_{{\gamma}_{m,d},w}=\begin{cases}&\displaystyle{\int}f\left(x,\phi(x)\right)g\left(x,\phi(x)\right)w(x)\mathrm{d}x,\hskip
61.17325ptm=1,\\\
&\displaystyle{\int}\left[f\left(x,\sqrt{\phi(x)}\right)g\left(x,\sqrt{\phi(x)}\right)\right.\\\
&\left.\>\>\>+f\left(x,-\sqrt{\phi(x)}\right)g\left(x,-\sqrt{\phi(x)}\right)\right]w(x)\mathrm{d}x,\>\>\>m=2,\end{cases}$

where the integrals are over $\mathrm{supp}(w)$. The bilinear form
${\langle}\cdot,\cdot{\rangle}_{{\gamma}_{m,d},w}$ defines an inner product on
the space of polynomials of two variables restricted to the curve
${\gamma}_{m,d}$,

$\Pi({\gamma}_{m,d}):={\mathbb{R}}[x,y]/{\langle}y^{m}-\phi(x){\rangle},$

which denotes the space of bivariate polynomials modulo the ideal generated by
$y^{m}-\phi(x)=0$.

######  Proposition 3.1.

Let ${\mathcal{V}}_{n}({\gamma}_{m,d},w)$ (or ${\mathcal{V}}_{n}$) be the
space of bivariate OPs of degree $n$ with respect to the inner product (3.2)
defined on $\Pi({\gamma}_{m,d})$, then

(3.3) $\dim{\mathcal{V}}_{n}({\gamma}_{m,d},w)=\min\\{n+1,\max\\{d,m\\}\\}.$

###### Proof.

We first consider the dimension of the space of monomials of degree $n$ on the
curve ${\gamma}_{m,d}$, which we denote by
${\mathcal{B}}_{n}({\gamma}_{m,d})$. Three cases arise: (i)
$n\leq\max\\{d,m\\}-1$, (ii) $m>d,n\geq m$ and (iii) $m\leq d$, $n\geq d$.
Since we only let $m=1,2$ and $d>0$, case (ii) only arises if $m=2$ and $d=1$.

For case(i), we have that

(3.4)
${\mathcal{B}}_{n}({\gamma}_{m,d})=\displaystyle{\left\\{x^{n-k}y^{k}\right\\}_{k=0}^{n}},\qquad
n\leq\max\\{d,m\\}-1,$

because the monomials $x^{k}y^{n-k}$ are linearly independent for $0\leq k\leq
n$ and thus $\dim{\mathcal{B}}_{n}({\gamma}_{m,d})=n+1$.

For case (ii),

(3.5)
${\mathcal{B}}_{n}({\gamma}_{m,d})=\left\\{x^{n-k}y^{k}\right\\}_{k=0}^{m-1},\qquad
m>d,\qquad n\geq m,$

because we can set $y^{m}=\phi(x)$ and express the monomials $x^{n-k}y^{k}$
with $k\geq m$ in terms of the monomials in
${\mathcal{B}}_{k}({\gamma}_{m,d})$, with $0\leq k\leq n$, hence only the
monomials $x^{n-k}y^{k}$ with $0\leq k\leq m-1$ are linearly independent and
thus $\dim{\mathcal{B}}_{n}({\gamma}_{m,d})=m$.

For the final case (iii),

(3.6)
${\mathcal{B}}_{n}({\gamma}_{m,d})=\left\\{x^{d-1}y^{n+1-d},x^{d-2}y^{n+2-d},\ldots,y^{n}\right\\},\qquad
d\geq m,\qquad n\geq d,$

because we can use the fact that $y^{m}=c_{d}x^{d}+\cdots+c_{1}x+c_{0}$, with
$c_{d}\neq 0$ to set $x^{d}=\psi(x,y)$, where $\deg\psi<d$, and express
$y^{n-k}x^{k}$ with $k\geq d$ in terms of lower degree monomials, therefore
$\dim{\mathcal{B}}_{n}({\gamma}_{m,d})=d$

Applying the Gram–Schmidt process to the basis
${\mathcal{B}}_{n}({\gamma}_{m,d})$ in cases (i)–(iii), it follows that
$\dim{\mathcal{B}}_{n}({\gamma}_{m,d})=\dim{\mathcal{V}}_{n}({\gamma}_{m,d},w)$,
which completes the proof. ∎

######  Proposition 3.2.

A bivariate polynomial $p(x,y)\in\Pi({\gamma}_{m,d})$ can be represented in
terms of polynomials of one variable as follows,

(3.7) $p(x,y)=\begin{cases}q(x)&\text{if }m=1,\\\ r(x)+ys(x)&\text{if
}m=2,\end{cases}$

where $q(x),r(x),s(x)$ are polynomials in $x$.

###### Proof.

Any $p(x,y)\in\Pi({\gamma}_{m,d})$ can be expressed as a linear combination of
the elements of the monomial bases in (3.4)–(3.6). The result follows
immediately upon setting $y=\phi(x)$ (if $m=1$) or $y^{2}=\phi(x)$ (if $m=2$)
in the monomial expansion of $p(x,y)$. ∎

In view of this result, we introduce the “$P$-polynomials” or “$P$-basis”

(3.8) $\displaystyle P^{(m,d)}_{n,k}(x,y)=\begin{cases}p_{n}(w;x)&\text{for
}m=1,k=0,n\geq 0\\\ p_{n}(w;x)/\sqrt{2}&\text{for }m=2,k=0,n\geq 0\\\
yp_{n-1}(\phi w;x)/\sqrt{2}&\text{for }m=2,k=1,n\geq 1\end{cases},$

which clearly form a basis for $\Pi({\gamma}_{m,d})$ and, moreover:

######  Proposition 3.3.

The $P$-polynomials are orthonormal with respect to
$\langle\cdot,\cdot\rangle_{{\gamma}_{m,d},w}$:

$\langle
P^{(m,d)}_{n,k},P^{(m,d)}_{\ell,i}\rangle_{{\gamma}_{m,d},w}=\delta_{n,\ell}\delta_{k,i},\qquad
m=1,2.$

###### Proof.

Since the $p_{n}(w)$ are orthonormal with respect to the univariate inner
product $\langle\cdot,\cdot\rangle_{w}$ we have, for $m=1$,

$\langle P^{(1,d)}_{n,0},P^{(1,d)}_{\ell,0}\rangle_{{\gamma}_{1,d},w}=\langle
p_{n}(w),p_{\ell}(w)\rangle_{w}=\delta_{n,\ell}.$

For $m=2$, three orthogonality relations need to be demonstrated. First,
orthogonality between $P$-polynomials with $k=0$ in (3.8):

$\langle P^{(2,d)}_{n,0},P^{(2,d)}_{\ell,0}\rangle_{{\gamma}_{2,d},w}=\langle
p_{n}(w),p_{\ell}(w)\rangle_{w}=\delta_{n,\ell}.$

Next, between polynomials with $k=1$,

$\langle
P^{(2,d)}_{n,1},P^{(2,d)}_{\ell,1}\rangle_{{\gamma}_{2,d},w}=\langle\phi
p_{n-1}(\phi w),p_{\ell-1}(\phi w)\rangle_{w}=\langle p_{n-1}(\phi
w),p_{\ell-1}(\phi w)\rangle_{\phi w}=\delta_{n,\ell}.$

The final orthogonality relation is due to the reflection symmetry of the
curve for $m=2$:

$\langle
P^{(2,d)}_{n,1},P^{(2,d)}_{\ell,0}\rangle_{{\gamma}_{2,d},w}=\int\left[yp_{n-1}(\phi
w)p_{\ell}(w)-yp_{n-1}(\phi w)p_{\ell}(w)\right]w\mathrm{d}x=0.$

∎

###  3.1. Connection coefficients

We denote the degree-$n$, orthonormal OPs on ${\gamma}_{m,d}$ by
$Y^{(m,d)}_{n,k}(x,y)$, where $k=1,\ldots,\dim{\mathcal{V}}_{n}$, with
$\dim{\mathcal{V}}_{n}$ defined in (3.3). Hence, the OPs satisfy

$\langle
Y^{(m,d)}_{n,k},Y^{(m,d)}_{\ell,i}\rangle_{{\gamma}_{m,d},w}=\delta_{n,\ell}\delta_{k,i},\qquad
m=1,2,$

for $n,\ell\geq 0$, $1\leq k\leq\dim{\mathcal{V}}_{k}$, $1\leq
i\leq\dim{\mathcal{V}}_{i}$ and $Y^{(m,d)}_{n,k}$ has degree $n$ on
${\gamma}_{m,d}$.

It follows from Proposition 3.2 that the OPs can be expressed as linear
combinations of the $P$-polynomials. That is, there exist _connection
coefficients_ $C^{(m,d)}_{j,\ell}$, where $\ell=\ell(n,k)\in\mathbb{N}_{0}$
will be defined in section 5.5, such that

(3.9) $\displaystyle Y_{n,k}^{(1,d)}(x,y)$ $\displaystyle=\sum_{j\geq
1}C^{(1,d)}_{j,\ell}P_{j,0}^{(1,d)},$ (3.10) $\displaystyle
Y_{n,k}^{(2,d)}(x,y)$ $\displaystyle=\sum_{j\geq
1}C^{(2,d)}_{2j-1,\ell}P_{j,0}^{(2,d)}+\sum_{j\geq
1}C^{(2,d)}_{2j,\ell}P_{j,1}^{(2,d)},\quad$

for $n\geq 1$, $k=1,\ldots,\dim{\mathcal{V}}_{n}$ and with
$Y_{0,1}^{(m,d)}=C^{(m,d)}_{0,0}P_{0,0}^{(m,d)}$ for $m=1,2$. The upper and
lower bounds of the summation indices in (3.9)–(3.10) will be considered in
section 5.6.

###  3.2. Minimal and maximal degrees of polynomials on ${\gamma}_{m,d}$

If $p$ is a polynomial in $\mathbb{R}[x,y]$, then we define the degree of the
restriction of $p$ to ${\gamma}_{m,d}$ as the lowest degree of $q$ among all
$q$ satisfying $p=q$ modulo $\langle y^{m}-\phi(x)\rangle$. By construction,
the monomial bases in (3.4)–(3.6) are minimal degree representations of
polynomials in $\Pi({\gamma}_{m,d})$. Therefore, the degree of $p(x,y)$ on the
curve ${\gamma}_{m,d}$ is the degree of its expansion in the monomial basis
given in (3.4)–(3.6).

If $d>m$ and $p(x,y)\in\Pi({\gamma}_{m,n})$, then a maximal degree
representation of $p$ is obtained by setting $y=\phi$ (if $m=1$) or
$y^{2}=\phi$ (if $m=2$), in which case $p$ takes the form given in (3.7). For
example, if $m=1<d$, then $y^{n}=\phi^{n}$ is a polynomial of degree $n$ in
$\Pi({\gamma}_{m,n})$ whose maximal degree on ${\gamma}_{1,d}$ is $nd$.
Observe from (3.8) that the total degree of a $P$-polynomial $P^{(m,d)}_{n,k}$
is $n$, which is also its maximal degree since it is in the form (3.7).
Therefore, the maximal degree of $p(x,y)\in\Pi({\gamma}_{m,n})$ is the total
degree of its expansion in the $P$-basis. We denote the maximal degree of
$p(x,y)$ on ${\gamma}_{m,d}$ by $\deg p$ and its minimal degree by
$\deg_{{\gamma}_{m,d}}(p)$. For example, if $m=1<d$, then $\deg y^{n}=nd$ and
$\deg_{{\gamma}_{1,d}}(y^{n})=n$ and for $n>d$, $\deg
P^{(m,d)}_{n,k}=n>\deg_{{\gamma}_{m,d}}(P^{(m,d)}_{n,k})$ while, by
definition, $\deg_{{\gamma}_{m,d}}(Y^{(m,d)}_{n,k})=n$.

##  4\. Curves with explicit OP bases

Before recursively constructing OPs on the general curves (3.1), we consider
cases for which the OP basis can be given explicitly in terms of the
$P$-polynomials.

######  Proposition 4.1.

The following are orthonormal OP bases on curves ${\gamma}_{m,d}$ defined in
(3.1): if $m=1$ and $\deg\phi=1$, then

(4.1) $\displaystyle Y_{n,1}^{(1,1)}(x,y)=P_{n,0}^{(1,1)}(x),\qquad n\geq 0;$

if $m=1$ and $\deg\phi=2$, then

(4.6) $\displaystyle
Y_{0,1}^{(1,2)}(x,y)=P_{0,0}^{(1,2)}(x),\qquad\left(\begin{array}[]{c}Y^{(1,2)}_{n,1}(x,y)\\\
Y^{(1,2)}_{n,2}(x,y)\end{array}\right)=\left(\begin{array}[]{l}P_{2n-1,0}^{(1,2)}(x)\\\
P_{2n,0}^{(1,2)}(x)\end{array}\right).$

If $m=2$ and $\deg\phi=1$ or $\deg\phi=2$, then

(4.11) $\displaystyle
Y_{0,1}^{(2,d)}(x,y)=P_{0,0}^{(2,d)}(x),\qquad\left(\begin{array}[]{c}Y^{(2,d)}_{n,1}(x,y)\\\
Y^{(2,d)}_{n,2}(x,y)\end{array}\right)=\left(\begin{array}[]{l}P_{n,0}^{(2,d)}(x)\\\
P_{n,1}^{(2,d)}(x,y)\end{array}\right),\qquad d=1,2.$

If $m=2$ and $\deg\phi=3$, then $Y_{0,1}^{(2,3)}=P_{0,0}^{(2,3)}$,
$Y_{1,1}^{(2,3)}=P_{1,0}^{(2,3)}$, $Y_{1,2}^{(2,3)}=P_{1,1}^{(2,3)}$ and

(4.24) $\displaystyle\left(\begin{array}[]{c}Y_{2n,1}^{(2,3)}\\\
Y_{2n,2}^{(2,3)}\\\
Y_{2n,3}^{(2.3)}\end{array}\right)=\left(\begin{array}[]{l}P_{3n-1,0}^{(2,3)}\\\
P_{3n-1,1}^{(2,3)}\\\
P_{3n,0}^{(2,3)}\end{array}\right),\qquad\left(\begin{array}[]{c}Y_{2n+1,1}^{(2,3)}\\\
Y_{2n+1,2}^{(2,3)}\\\
Y_{2n+1,3}^{(2,3)}\end{array}\right)=\left(\begin{array}[]{l}P_{3n,1}^{(2,3)}\\\
P_{3n+1,0}^{(2,3)}\\\ P_{3n+1,1}^{(2,3)}\end{array}\right),\quad n\geq 1.$

If $m=2$, $\deg\phi=4$, $\phi$ is even, $\mathrm{supp}(w)=(-a,a)$, $a>0$ and
the weight $w$ is even, so that $P^{(2,4)}_{n,0}$ has the same parity as $n$,
then $Y_{n,i}^{(2,4)}=P_{n,i-1}^{(2,4)}$ for $(n,i)=(0,1)$, $(1,1)$, $(1,2)$,
$(2,1)$, $(2,2)$, $Y_{2,3}^{(2,4)}=P_{4,0}^{(2,4)}$ and for $n\geq 1$

(4.41) $\displaystyle\left(\begin{array}[]{c}Y_{2n+1,1}^{(2,4)}\\\
Y_{2n+1,2}^{(2,4)}\\\ Y_{2n+1,3}^{(2,4)}\\\
Y_{2n+1,4}^{(2,4)}\end{array}\right)=\left(\begin{array}[]{l}P_{4n-1,0}^{(2,4)}\\\
P_{4n-1,1}^{(2,4)}\\\ P_{4n+1,0}^{(2,4)}\\\
P_{4n+1,1}^{(2,4)}\end{array}\right),\qquad\left(\begin{array}[]{c}Y_{2n+2,1}^{(2,4)}\\\
Y_{2n+2,2}^{(2,4)}\\\ Y_{2n+2,3}^{(2,4)}\\\
Y_{2n+2,4}^{(2,4)}\end{array}\right)=\left(\begin{array}[]{l}P_{4n,1}^{(2,4)}\\\
P_{4n+2,0}^{(2,4)}\\\ P_{4n+2,1}^{(2,4)}\\\
P_{4n+4,0}^{(2,4)}\end{array}\right).$

###### Proof.

The sets of $Y^{(m,d)}_{n,k}$ given above are orthonormal (see Proposition
3.3) and the dimensions of ${\mathcal{V}}_{n}$ are as in Proposition 3.1. What
remains to be shown is that the $Y^{(m,d)}_{n,k}$ have degree $n$ on
${\gamma}_{m,d}$. We only prove the latter case (4.41) since the proofs for
the remaining cases are similar.

Since $\phi$ is assumed to be even and of degree 4,
$y^{2}=c_{4}x^{4}+c_{2}x^{2}+c_{0}$, with $c_{4}\neq 0$. To show that
$Y_{2n+1,1}^{(2,4)}=P^{(2,4)}_{4n-1,0}$ has degree $2n+1$ on ${\gamma}_{2,4}$,
it suffices to consider the highest degree monomial of $P^{(2,4)}_{4n-1,0}$,
viz. $x^{4n-1}$ and note that

$x^{4n-1}=x^{3}\left(x^{4}\right)^{n-1}=x^{3}\left(\frac{y^{2}-c_{2}x^{2}-c_{0}}{c_{4}}\right)^{n-1},$

which is an expression of degree $2n+1$ whose degree cannot be further reduced
on the curve ${\gamma}_{2,4}$. We conclude that
$\deg_{{\gamma}_{2,4}}\left(P^{(2,4)}_{4n-1,0}\right)=2n+1$. Similarly,
considering $Y_{2n+1,2}^{(2,4)}=P^{(2,4)}_{4n-1,1}$, the highest degree
monomial is
$yx^{4n-2}=yx^{2}\left(x^{4}\right)^{n-1}=yx^{2}[(y^{2}-c_{2}x^{2}-c_{0})/c_{4}]^{n-1}$
and it follows that
$\deg_{{\gamma}_{2,4}}\left(P^{(2,4)}_{4n-1,1}\right)=2n+1$. The remaining
cases can be similarly confirmed.

∎

The OP bases on quadratic curves (4.11) and cubic curves (4.24) are the cases
studied in, respectively, [14] and [6].

##  5\. OPs on $y=\phi$, $\deg\phi\geq 3$ and $y^{2}=\phi$, $\deg\phi\geq 4$

We now consider the construction of OPs on $y=\phi$ and $y^{2}=\phi$ in cases
for which we are unable to give an explicit OP basis, which are the curves
$y=\phi$, $\deg\phi\geq 3$ and $y^{2}=\phi$, $\deg\phi\geq 4$.

By the same reasoning used in the proof of Proposition 4.1, we deduce that

(5.1) $Y_{n,k}^{(1,d)}=P_{n,k-1}^{(1,d)},\qquad(n,k)=(0,1),(1,1),$

for $d\geq 3$ and

(5.2)
$Y_{n,k}^{(2,d)}=P_{n,k-1}^{(2,d)},\qquad(n,k)=(0,1),(1,1),(1,2),(2,1),(2,2),$

for $d\geq 4$. To generate higher degree OPs, we shall apply the Gram–Schmidt
process and show that it is equivalent to the Lanczos algorithm applied to two
symmetric, commuting matrices.

We introduce the following infinite quasimatrices for the $P$-bases,

(5.4) $\displaystyle\mathbf{P}^{(1,d)}=\mathbf{P}^{(1,d)}(x)$ $\displaystyle=\left(\begin{array}[]{c c c c}P^{(1,d)}_{0,0}(x)&P^{(1,d)}_{1,0}(x)&P^{(1,d)}_{2,0}(x)&\cdots\end{array}\right),$ (5.6) $\displaystyle\mathbf{P}^{(2,d)}=\mathbf{P}^{(2,d)}(x,y)$ $\displaystyle=\left(\begin{array}[]{c | c | c | c}{\mathbb{P}}^{(2,d)}_{0}(x)&{\mathbb{P}}^{(1,d)}_{1}(x,y)&{\mathbb{P}}^{(1,d)}_{2}(x,y)&\cdots\end{array}\right),$

where the ${\mathbb{P}}^{(2,d)}_{n}$ denote the block-quasimatrices

(5.7)
${\mathbb{P}}^{(2,d)}_{0}(x)=\left(\begin{array}[]{c}P_{0,0}^{(2,d)}(x)\end{array}\right),\quad{\mathbb{P}}^{(2,d)}_{n}(x,y)=\left(\begin{array}[]{c
c}P_{n,0}^{(2,d)}(x)&P_{n,1}^{(2,d)}(x,y)\end{array}\right),\quad n\geq 1.$

Similarly, for the OPs,

(5.8) $\mathbf{Y}^{(m,d)}=\mathbf{Y}^{(m,d)}(x,y)=\left(\begin{array}[]{c | c | c | c}{\mathbb{Y}}^{(m,d)}_{0}&{\mathbb{Y}}^{(m,d)}_{1}&{\mathbb{Y}}^{(m,d)}_{2}&\cdots\end{array}\right),$

with

(5.9) ${\mathbb{Y}}^{(m,d)}_{n}(x,y)=\left(\begin{array}[]{c c
c}Y_{n,1}^{(m,d)}(x,y)&\cdots&Y_{n,\dim{\mathcal{V}}_{n}}^{(m,d)}(x,y)\end{array}\right),\qquad
n\geq 0,$

where, since $m<d=\deg\phi$ for the cases we consider in this section,

$\dim{\mathcal{V}}_{n}=\min\\{n+1,d\\},$

see Proposition 3.1.

We can express (3.9)–(3.10) as

(5.10) $\mathbf{Y}^{(m,d)}=\mathbf{P}^{(m,d)}C^{(m,d)},$

where $C^{(m,d)}$ is referred to as the connection matrix and its entries
$C^{(m,d)}_{i,j}$, $i,j\geq 0$ are the connection coefficients. For all the
values of $(m,d)$ considered in Proposition 4.1, $C^{(m,d)}=I$, except for the
final case (4.41), where there exists a permutation matrix $\mathcal{P}$ such
that $C^{(2,4)}\mathcal{P}=I$. Figures 1 and 2 depict examples of connection
matrices for $m=1,\deg\phi\geq 3$ and $m=2,\deg\phi\geq 4$, the construction
of which we shall describe in the rest of this section.

Figure 1. Plots of $\log_{10}$ of the entries of the connection matrices
$C^{(1,d)}$ for OPs on $y=\phi(x)$. The curves are defined by
$\phi(x)=x^{3}-x^{2}+x+1$ (top-left), $\phi(x)=x^{4}+x^{3}-x^{2}+x+1$ (top-
right), $\phi(x)=(x^{6}+x^{5}-x^{4}+x^{3}-x^{2}+x+1)/10$ (bottom-left) and
$\phi(x)=x^{6}-x^{4}-x^{2}+1$ (bottom right). The uniform weight $w(x)=1$ on
$[-1,1]$ was used in all cases.

Figure 2. Plots of $\log_{10}$ of the entries of the connection matrices
$C^{(2,d)}$ for OPs on $y^{2}=\phi(x)$. The curves in these examples are
$\phi=(2+x)(1-x^{2})(5/4-x)$ (top left),
$\phi=\phi_{7}:=(1.15+x)(1.05+x)(1-x^{2})(1.1-x)(1.2-x)(1.3-x)$ (top right),
$\phi=(1.25+x)\phi_{7}$ (bottom left) and
$\phi=(1-x^{2})(1.05^{2}-x^{2})(1.1^{2}-x^{2})(1.15^{2}-x^{2})$ (bottom
right). The uniform weight $w=1$ on $[-1,1]$ was used in all cases.

######  Proposition 5.1.

The connection matrix $C^{(m,d)}$ for $m=1,2$ is orthogonal.

###### Proof.

Since the $\mathbf{Y}^{(m,d)}$ are orthonormal OPs and the
$\mathbf{P}^{(m,d)}$ are orthonormal for $m=1,2$, it follows that
${\langle}\mathbf{Y}^{(m,d)},\mathbf{Y}^{(m,d)}{\rangle}_{{\gamma}_{m,d},w}=I$
and
${\langle}\mathbf{P}^{(m,d)},\mathbf{P}^{(m,d)}{\rangle}_{{\gamma}_{m,d},w}=I$,
where $\langle\cdot,\cdot\rangle_{{\gamma}_{m,d},w}$ for quasimatrices is
defined as in (2.11)–(2.13), mutatis mutandis. Thus, the connection matrix is
orthogonal because

$\displaystyle I$
$\displaystyle={\langle}\mathbf{Y}^{(m,d)},\mathbf{Y}^{(m,d)}{\rangle}_{{\gamma}_{m,d},w}=\left(C^{(m,d)}\right)^{\top}{\langle}\mathbf{P}^{(m,d)},\mathbf{P}^{(m,d)}{\rangle}_{{\gamma}_{m,d},w}C^{(m,d)}$
$\displaystyle=\left(C^{(m,d)}\right)^{\top}C^{(m,d)},\qquad m=1,2.$

∎

###  5.1. Multiplication matrices for the ${P}$-bases

To construct the connection matrices we shall need matrices
$\mathcal{X}^{(m,d)}$ and $\mathcal{Y}^{(m,d)}$ representing multiplication of
$\mathbf{P}^{(m,d)}$ by $x$ and $y$:

(5.11) $x\mathbf{P}^{(m,d)}=\mathbf{P}^{(m,d)}\mathcal{X}^{(m,d)},\qquad
y\mathbf{P}^{(m,d)}=\mathbf{P}^{(m,d)}\mathcal{Y}^{(m,d)}.$

######  Proposition 5.2.

The multiplication matrices of the $\mathbf{P}^{(1,d)}$ basis are

(5.12) $\mathcal{X}^{(1,d)}=J(w),\qquad\mathcal{Y}^{(1,d)}=\phi(J(w)),$

where $J(w)$ is the Jacobi matrix of the orthonormal OPs $\\{p_{n}(w)\\}$,
which is defined as in (2.4)–(2.7) but with $\phi w$ replaced by $w$.
Furthermore, $\mathcal{X}^{(1,d)}$ and $\mathcal{Y}^{(1,d)}$ commute,
$\mathcal{X}^{(1,d)}$ is symmetric and tridiagonal and $\mathcal{Y}^{(1,d)}$
is symmetric with bandwidths $(d,d)$.

For $m=2$, $\mathcal{X}^{(2,d)}$ and $\mathcal{Y}^{(2,d)}$ commute,
$\mathcal{X}^{(2,d)}$ is symmetric, block-tridiagonal with diagonal blocks:

$\displaystyle\mathcal{X}^{(2,d)}=\left(\begin{array}[]{c c c c
c}\mathcal{X}_{0,0}^{(2,d)}&\left(\mathcal{X}_{1,0}^{(2,d)}\right)^{\top}&&&\\\
\mathcal{X}_{1,0}^{(2,d)}&\mathcal{X}_{1,1}^{(2,d)}&\left(\mathcal{X}_{2,1}^{(2,d)}\right)^{\top}&&\\\
&\mathcal{X}_{2,1}^{(2,d)}&\mathcal{X}_{2,2}^{(2,d)}&\left(\mathcal{X}_{3,2}^{(2,d)}\right)^{\top}&\\\
&&\ddots&\ddots&\ddots\end{array}\right),$

where
$\mathcal{X}_{0,0}^{(2,d)}=\left(\begin{array}[]{c}{\alpha}_{0}(w)\end{array}\right)\in{\mathbb{R}}^{1\times
1}$, $\mathcal{X}_{1,0}^{(2,d)}=\left(\begin{array}[]{c
c}{\beta}_{0}(w)&0\end{array}\right)^{\top}\in{\mathbb{R}}^{2\times 1}$, and
for $n\geq 1$

$\displaystyle\mathcal{X}_{n,n}^{(2,d)}=\left(\begin{array}[]{c
c}{\alpha}_{n}(w)&0\\\ 0&{\alpha}_{n-1}(\phi
w)\end{array}\right),\qquad\mathcal{X}_{n+1,n}^{(2,d)}=\left(\begin{array}[]{c
c}{\beta}_{n}(w)&0\\\ 0&{\beta}_{n-1}(\phi w)\end{array}\right).$

The matrix $\mathcal{Y}^{(2,d)}$ is symmetric with block-bandwiths
$(d-1,d-1)$:

(5.19) $\displaystyle\mathcal{Y}^{(2,d)}=\left(\begin{array}[]{ c c c c c c c
c c}0&\left(\mathcal{Y}_{1,0}^{(2,d)}\right)^{\top}&&&&&&&\\\
\mathcal{Y}_{1,0}^{(2,d)}&\mathcal{Y}_{1,1}^{(2,d)}&\left(\mathcal{Y}_{2,1}^{(2,d)}\right)^{\top}&\cdots&\left(\mathcal{Y}_{d,1}^{(2,d)}\right)^{\top}&&&&\\\
&\mathcal{Y}_{2,1}^{(2,d)}&\mathcal{Y}_{2,2}^{(2,d)}&\cdots&\cdots&\left(\mathcal{Y}_{d+1,2}^{(2,d)}\right)^{\top}&&&\\\
&\vdots&\ddots&\ddots&\ddots&\ddots&\ddots&&\\\
&\mathcal{Y}_{d,1}^{(2,d)}&\cdots&\cdots&\mathcal{Y}_{d,d}^{(2,d)}&\left(\mathcal{Y}_{d+1,d}^{(2,d)}\right)^{\top}&\cdots&\left(\mathcal{Y}_{2d-1,d}^{(2,d)}\right)^{\top}&\\\
&&\ddots&\ddots&\ddots&\ddots&\ddots&\ddots&\ddots\end{array}\right),$

where $\mathcal{Y}_{1,0}^{(2,d)}=\left(\begin{array}[]{c
c}0&r_{0,0}\end{array}\right)^{\top}$,

$\displaystyle\mathcal{Y}^{(2,d)}_{n,n}=\left(\begin{array}[]{c
c}0&r_{n-1,n}\\\
r_{n-1,n}&0\end{array}\right),\quad\mathcal{Y}^{(2,d)}_{n+1,n}=\left(\begin{array}[]{c
c}0&r_{n-1,n+1}\\\ r_{n,n}&0\end{array}\right),\quad n\geq 1,$
$\displaystyle\mathcal{Y}^{(2,d)}_{n+k,n}=\left(\begin{array}[]{c
c}0&r_{n-1,n+k}\\\ 0&0\end{array}\right),\qquad 2\leq k\leq d-1,\qquad n\geq
1,$

and the $r_{i,j}$ are defined in (2.10).

###### Proof.

Since $x^{k}\mathbf{P}(w)=\mathbf{P}(w)\left(J(w)\right)^{k}$, $k\geq 0$ and
$\left(J(w)\right)^{k}$ is symmetric with bandwidths $(k,k)$, it follows that
$y\mathbf{P}(w)=\phi(x)\mathbf{P}(w)=\mathbf{P}(w)\phi(J(w))$, where
$\phi(J(w))$ is symmetric with bandwidths $(d,d)$ because $\phi$ is a
polynomial of degree $d$. The result (5.12) then follows since
$\mathbf{P}^{(1,d)}=\mathbf{P}(w)$.

The commutativity of $\mathcal{X}^{(m,d)}$ and $\mathcal{Y}^{(m,d)}$ for
$m=1,2$ follows immediately from their definitions in (5.11) and the
commutativity of multiplication by $x$ and $y$,
$xy\mathbf{P}^{(m,d)}=yx\mathbf{P}^{(m,d)}$.

The entries of $\mathcal{X}^{(2,d)}$ follow straightforwardly from (3.8) and
the three-term recurrences of the OP families $\\{p_{n}(w)\\}$ and
$\\{p_{n}(\phi w)\\}$.

To construct $\mathcal{Y}^{(2,d)}$, we first note from (3.8) and (2.8)–(2.9)
that

(5.20) $\displaystyle yP_{n,0}^{(2,d)}$
$\displaystyle=r_{n-d,n}P_{n-d+1,1}^{(2,d)}+\cdots+r_{n,n}P_{n+1,1}^{(2,d)},$

for $n\geq 0$, where $r_{i,j}=0$ if $i<0$. Since
$yP_{n+1,1}^{(2,d)}=y^{2}p_{n}(\phi w)/\sqrt{2}=\phi p_{n}(\phi w)/\sqrt{2}$,
we can express $yP_{n+1,1}^{(2,d)}$ in the $\mathbf{P}^{(2,d)}$ basis using
the entries of the lowering matrix $L$ defined by

$\phi(x)\mathbf{P}(\phi w)=\mathbf{P}(w)L,$

where (cf. (2.11))

$\displaystyle L=\langle\mathbf{P}(w),\phi\mathbf{P}(\phi
w)\rangle_{w}=\int\mathbf{P}^{\top}(w)\mathbf{P}(\phi w)\phi
w\mathrm{d}x=R^{\top},$

and therefore

(5.21) $\displaystyle yP_{n+1,1}^{(2,d)}$
$\displaystyle=r_{n,n}P_{n,0}^{(2,d)}+\cdots+r_{n,n+d}P_{n+d,0}^{(2,d)}.$

The entries and structure of $\mathcal{Y}^{(2,d)}$ then follow from (5.20) and
(5.21). ∎

###  5.2. The Gram–Schmidt process viewed as a block Lanczos algorithm

Since the $Y_{n,k}^{(m,d)}$ are orthonormal OPs, the
${\mathbb{Y}}_{n}^{(m,d)}$ satisfy three-term recurrences of the form [5]

(5.22) $\displaystyle
x{\mathbb{Y}}_{n}^{(m,d)}={\mathbb{Y}}_{n-1}^{(m,d)}\left(B_{n-1,x}^{(m,d)}\right)^{\top}+{\mathbb{Y}}_{n}^{(m,d)}A_{n,x}^{(m,d)}+{\mathbb{Y}}_{n+1}^{(m,d)}B_{n,x}^{(m,d)},$

and

(5.23) $\displaystyle
y{\mathbb{Y}}_{n}^{(m,d)}={\mathbb{Y}}_{n-1}^{(m,d)}\left(B_{n-1,y}^{(m,d)}\right)^{\top}+{\mathbb{Y}}_{n}^{(m,d)}A_{n,y}^{(m,d)}+{\mathbb{Y}}_{n+1}^{(m,d)}B_{n,y}^{(m,d)},$

for $n\geq 0$ with ${\mathbb{Y}}_{-1}:=0$, where

(5.24) $\displaystyle
A_{n,x}^{(m,d)}={\langle}{\mathbb{Y}}_{n}^{(m,d)},x{\mathbb{Y}}_{n}^{(m,d)}{\rangle}_{{\gamma}_{m,d},w},\quad
B_{n,x}^{(m,d)}={\langle}{\mathbb{Y}}_{n+1}^{(m,d)},x{\mathbb{Y}}_{n}^{(m,d)}{\rangle}_{{\gamma}_{m,d},w},$

and $A_{n,y}^{(m,d)}$, $B_{n,y}^{(m,d)}$ are defined similarly. Hence,
$A_{n,x}^{(m,d)},A_{n,y}^{(m,d)}\in{\mathbb{R}}^{\dim{\mathcal{V}}_{n}\times\dim{\mathcal{V}}_{n}}$,
$B_{n,x}^{(m,d)},B_{n,y}^{(m,d)}\in{\mathbb{R}}^{\dim{\mathcal{V}}_{n+1}\times\dim{\mathcal{V}}_{n}}$
and $A_{n,x}^{(m,d)}$ and $A_{n,y}^{(m,d)}$ are symmetric. In quasimatrix
notation, the recurrences (5.22) and (5.23) can be expressed as

$x\mathbf{Y}^{(m,d)}=\mathbf{Y}^{(m,d)}J_{x}^{(m,d)},\qquad
y\mathbf{Y}^{(m,d)}=\mathbf{Y}^{(m,d)}J_{y}^{(m,d)},$

where $J_{x}^{(m,d)}$ and $J_{y}^{(m,d)}$ are the symmetric, commuting, block-
tridiagonal Jacobi matrices of the OPs $\mathbf{Y}^{(m,d)}$.

To express the recurrences (5.22) and (5.23) in the $\mathbf{P}^{(m,d)}$
bases, we write the connection matrix in block form

(5.25) $C^{(m,d)}=\left(\begin{array}[]{c | c | c | c}{\mathbb{C}}_{0}^{(m,d)}&{\mathbb{C}}_{1}^{(m,d)}&{\mathbb{C}}_{2}^{(m,d)}&\cdots\end{array}\right),$

where
${\mathbb{C}}_{n}^{(m,d)}\in\mathbb{R}^{\infty\times\dim{\mathcal{V}}_{n}}$,
then

(5.26) ${\mathbb{Y}}_{n}^{(m,d)}=\mathbf{P}^{(m,d)}{\mathbb{C}}_{n}^{(m,d)}.$

######  Proposition 5.3.

The Jacobi matrices of the OPs $\mathbf{Y}^{(m,d)}$ and the multiplication
matrices of the $\mathbf{P}^{(m,d)}$ bases are related as follows

(5.27)
$J_{x}^{(m,d)}=\left(C^{(m,d)}\right)^{\top}\mathcal{X}^{(m,d)}C^{(m,d)},\qquad
J_{y}^{(m,d)}=\left(C^{(m,d)}\right)^{\top}\mathcal{Y}^{(m,d)}C^{(m,d)},\quad
m=1,2.$

Furthermore, the ${\mathbb{C}}_{n}^{(m,d)}$ satisfy the recurrence relations

(5.28)
$\displaystyle\mathcal{X}^{(m,d)}{\mathbb{C}}_{n}^{(m,d)}={\mathbb{C}}_{n-1}^{(m,d)}\left(B_{n-1,x}^{(m,d)}\right)^{\top}+{\mathbb{C}}_{n}^{(m,d)}A_{n,x}^{(m,d)}+{\mathbb{C}}_{n+1}^{(m,d)}B_{n,x}^{(m,d)},$

and

(5.29)
$\displaystyle\mathcal{Y}^{(m,d)}{\mathbb{C}}_{n}^{(m,d)}={\mathbb{C}}_{n-1}^{(m,d)}\left(B_{n-1,y}^{(m,d)}\right)^{\top}+{\mathbb{C}}_{n}^{(m,d)}A_{n,y}^{(m,d)}+{\mathbb{C}}_{n+1}^{(m,d)}B_{n,y}^{(m,d)}.$

###### Proof.

Since
$x\mathbf{Y}^{(m,d)}=x\mathbf{P}^{(m,d)}C^{(m,d)}=\mathbf{P}^{(m,d)}\mathcal{X}^{(m,d)}C^{(m,d)}$,
we have that

$\displaystyle
J_{x}^{(m,d)}={\langle}\mathbf{Y}^{(m,d)},x\mathbf{Y}^{(m,d)}{\rangle}_{{\gamma}_{m,d},w}=$
$\displaystyle\left(C^{(m,d)}\right)^{\top}{\langle}\mathbf{P}^{(m,d)},\mathbf{P}^{(m,d)}{\rangle}_{{\gamma}_{m,d},w}\mathcal{X}^{(m,d)}C^{(m,d)}=\left(C^{(m,d)}\right)^{\top}\mathcal{X}^{(m,d)}C^{(m,d)},$

for $m=1,2$.

Using (5.26) and setting
$x{\mathbb{Y}}_{n}^{(m,d)}=\mathbf{P}^{(m,d)}\mathcal{X}^{(m,d)}{\mathbb{C}}_{n}^{(m,d)}$
in (5.22), we arrive at (5.28) and (5.29) can be derived similarly. ∎

The standard Stieltjes method for computing OPs and their recurrence
coefficients uses a discretised quadrature rule to compute the required inner
products, see [2, 7]. Hence, the (discretised) OPs need to be evaluated at a
discrete set of points. By contrast, we compute OPs in “coefficient space” (as
opposed to “function-value space”, as is done in the standard approach) via
(5.28) and (5.29) to generate subsequent blocks of the connection matrices.
Computation of the inner products (5.24) is accomplished via

(5.30) $\displaystyle
B_{n-1,x}^{(m,d)}=\left({\mathbb{C}}_{n}^{(m,d)}\right)^{\top}\mathcal{X}^{(m,d)}{\mathbb{C}}_{n-1}^{(m,d)},\qquad
A_{n,x}^{(m,d)}=\left({\mathbb{C}}_{n}^{(m,d)}\right)^{\top}\mathcal{X}^{(m,d)}{\mathbb{C}}_{n}^{(m,d)},$

which follows from (5.24) and (5.26). The inner products required for
$B_{n-1,y}^{(m,d)}$ and $A_{n,y}^{(m,d)}$ are computed similarly. Hence, once
we have the multiplication matrices $\mathcal{X}^{(m,d)}$ and
$\mathcal{Y}^{(m,d)}$ (which are built from the Jacobi matrices $J(w)$,
$J(\phi w)$ and the raising matrix $R$), we compute inner products via only
finite matrix multiplications, with no need for quadrature of the evaluation
of OPs.

Since $J_{x}^{(m,d)}$ and $J_{y}^{(m,d)}$ are symmetric, block-tridiagonal
matrices, we note from (5.27) that the connection matrices simultaneously
block-tridiagonalise the symmetric, commuting matrices $\mathcal{X}^{(m,d)}$
and $\mathcal{Y}^{(m,d)}$ via orthogonal similarity transformations. Hence,
our approach, which is essentially the Gram–Schmidt process applied in
coefficient space to OPs in the $\mathbf{P}^{(m,d)}$ bases via the connection
coefficients, can be viewed as the Lanczos procedure applied to
$\mathcal{X}^{(m,d)}$ and $\mathcal{Y}^{(m,d)}$.

###  5.3. Construction of the initial columns of the connection matrices

Let $\bm{e}_{k}$ be defined as in (2.12) and let $\mathbf{C}^{(m,d)}_{j}$
denote column $j$ of the connection matrix $C^{(m,d)}$, with $j\geq 0$.

######  Proposition 5.4.

The first few columns of the connection matrices are given by

(5.31)
${\mathbb{C}}_{0}^{(1,d)}=\mathbf{C}^{(1,d)}_{0}=\bm{e}_{0},\qquad{\mathbb{C}}_{1}^{(1,d)}=\left(\begin{array}[]{c
c}\bm{e}_{1}&\mathbf{C}^{(1,d)}_{2}\end{array}\right),$

and for $m=2$,

(5.32)
${\mathbb{C}}_{0}^{(2,d)}=\bm{e}_{0},\qquad{\mathbb{C}}_{1}^{(2,d)}=\left(\begin{array}[]{c
c}\bm{e}_{1}&\bm{e}_{2}\end{array}\right),\qquad{\mathbb{C}}_{2}^{(2,d)}=\left(\begin{array}[]{c
c c}\bm{e}_{3}&\bm{e}_{4}&\mathbf{C}^{(2,d)}_{5}\end{array}\right),$

where

(5.33)
$\mathbf{C}^{(1,d)}_{2}=\frac{1}{b_{2,1}^{0,y}}\left(\mathcal{Y}_{2,0}^{(1,d)}\bm{e}_{2}+\cdots+\mathcal{Y}_{d,0}^{(1,d)}\bm{e}_{d}\right),\quad\left(b_{2,1}^{0,y}\right)^{2}=\left(\mathcal{Y}_{2,0}^{(1,d)}\right)^{2}+\cdots+\left(\mathcal{Y}_{d,0}^{(1,d)}\right)^{2},$
(5.34)
$\mathbf{C}^{(2,d)}_{5}=\frac{1}{b_{3,2}^{1,y}}\left(r_{0,3}\bm{e}_{5}+\cdots+r_{0,d}\bm{e}_{2d-1}\right),\quad\left(b_{3,2}^{1,y}\right)^{2}=r_{0,3}^{2}+\cdots+r_{0,d}^{2},$

and $\mathcal{Y}_{i,j}^{(1,d)}$, $i,j\geq 0$ denotes the $(i,j)$ entry of
$\mathcal{Y}^{(1,d)}$ (defined in (5.12)). Furthermore, the first few block-
entries of the Jacobi matrices $J_{x}^{(m,d)}$ and $J_{y}^{(m,d)}$ are given
by (see (5.22)–(5.23))

(5.35)
$A_{0,x}^{(1,d)}=\left(\begin{array}[]{c}\mathcal{X}_{0,0}^{(1,d)}\end{array}\right),\qquad
B_{0,x}^{(1,d)}=\left(\begin{array}[]{c
c}\mathcal{X}_{1,0}^{(1,d)}&0\end{array}\right)^{\top},$ (5.36)
$A_{0,y}^{(1,d)}=\left(\begin{array}[]{c}\mathcal{Y}_{0,0}^{(1,d)}\end{array}\right),\qquad
B_{0,y}^{(1,d)}=\left(\begin{array}[]{c
c}\mathcal{Y}_{1,0}^{(1,d)}&b_{2,1}^{0,y}\end{array}\right)^{\top},$

and for $m=2$,

(5.37) $A_{0,x}^{(2,d)}=\mathcal{X}^{(2,d)}_{0,0},\quad
B_{0,x}^{(2,d)}=\mathcal{X}^{(2,d)}_{1,0},\quad
A_{1,x}^{(2,d)}=\mathcal{X}^{(2,d)}_{1,1},\quad
B_{1,x}^{(2,d)}=\left(\mathcal{X}^{(2,d)}_{2,1}\begin{array}[]{c}0\\\
0\end{array}\right)^{\top},$

and

(5.38) $A_{0,y}^{(2,d)}=\left(\begin{array}[]{c}0\end{array}\right),\quad
B_{0,y}^{(2,d)}=\mathcal{Y}^{(2,d)}_{1,0},\quad
A_{1,y}^{(2,d)}=\mathcal{Y}^{(2,d)}_{1,1},\quad
B_{1,y}^{(2,d)}=\left(\mathcal{Y}^{(2,d)}_{2,1}\begin{array}[]{c}0\\\
b_{3,2}^{1,y}\end{array}\right)^{\top}.$

###### Proof.

The first few columns of $C^{(m,d)}$ as well as the matrices (5.35)–(5.38)
follow from (5.1), (5.2) and (5.24). We shall obtain $\mathbf{C}^{(1,d)}_{2}$,
$\mathbf{C}^{(2,d)}_{5}$ and the entries $b_{2,1}^{0,y}$ and $b_{3,2}^{1,y}$
by constructing, respectively, $Y_{1,2}^{(1,d)}$ and $Y_{2,3}^{(2,d)}$.

The recurrence (5.23) and (5.36) imply

$yY_{0,1}^{(1,d)}=\mathcal{Y}_{0,0}^{(1,d)}Y_{0,1}^{(1,d)}+\mathcal{Y}_{1,0}^{(1,d)}Y_{1,1}^{(1,d)}+b_{2,1}^{0,y}Y_{1,2}^{(1,d)},$

on the other hand, using the multiplication matrix $\mathcal{Y}^{(1,d)}$ in
(5.12),

$\displaystyle yY_{0,1}^{(1,d)}$
$\displaystyle=\phi(x)P_{0,1}^{(1,d)}=\mathcal{Y}^{(1,d)}_{0,0}P_{0,0}^{(1,d)}+\cdots+\mathcal{Y}^{(1,d)}_{d,0}P_{d,0}^{(1,d)}$
$\displaystyle=\mathcal{Y}_{0,0}^{(1,d)}Y_{0,1}^{(1,d)}+\mathcal{Y}_{1,0}^{(1,d)}Y_{1,1}^{(1,d)}+\mathcal{Y}_{2,0}^{(1,d)}P_{2,0}^{(1,d)}+\cdots+\mathcal{Y}_{d,0}^{(1,d)}P_{d,0}^{(1,d)},$

comparing these two expressions, we have that

$\displaystyle
Y_{1,2}^{(1,d)}=\frac{1}{b_{2,1}^{0,y}}\left(\mathcal{Y}_{2,0}^{(1,d)}P_{2,0}^{(1,d)}+\cdots+\mathcal{Y}_{d,0}^{(1,d)}P_{d,0}^{(1,d)}\right)=\mathbf{P}^{(1,d)}\mathbf{C}^{(1,d)}_{2},$

from which (5.33) follows and the expression for $b_{2,1}^{0,y}$ follows from
the requirement that
$\|Y_{1,2}^{(1,d)}\|^{2}_{{\gamma}_{1,d},w}={\langle}Y_{1,2}^{(1,d)},Y_{1,2}^{(1,d)}{\rangle}_{{\gamma}_{1,d},w}=1$,
which implies that

${\langle}Y_{1,2}^{(1,d)},Y_{1,2}^{(1,d)}{\rangle}_{{\gamma}_{1,d},w}={\langle}\mathbf{P}^{(1,d)}\mathbf{C}^{(1,d)}_{2},\mathbf{P}^{(1,d)}\mathbf{C}^{(1,d)}_{2}{\rangle}_{{\gamma}_{1,d},w}=\left(\mathbf{C}^{(1,d)}_{2}\right)^{\top}\mathbf{C}^{(1,d)}_{2}=1;$

$\mathbf{C}^{(1,d)}_{5}$ and $b_{3,2}^{1,y}$ can be derived similarly by using
(5.38) and (5.21). ∎

Note from (5.33)–(5.34) that the entries $b^{0,y}_{2,1}$ and $b^{1,y}_{3,2}$,
whose values are determined by the normalisation of the OPs, are unique up to
a choice of sign. Throughout, we shall let entries that arise from
normalisation be positive.

###  5.4. Sequences of orthogonalisation for the subsequent columns of the
connection matrices

We now consider the sequence of multiplications by $x$ and $y$ to be used for
constructing the OPs. We first note that different sequences of multiplication
can result in different Jacobi matrices $J_{x}^{(m,d)}$, $J_{y}^{(m,d)}$ and
connection matrices $C^{(m,d)}$. For example, suppose $m=1$, we set
$Y^{(1,d)}_{0,1}=P_{0,0}^{(1,d)}$, first multiply by $y$ and orthogonalise to
obtain $Y^{(1,d)}_{1,1}$, and then multiply by $x$ and orthogonalise to obtain
$Y^{(1,d)}_{1,2}$, then, cf. (5.35) and (5.36),

$A_{0,y}^{(1,d)}=\left(\begin{array}[]{c}\mathcal{Y}_{0,0}^{(1,d)}\end{array}\right),\qquad
B_{0,y}^{(1,d)}=\left(\begin{array}[]{c
c}b_{1,1}^{0,y}&0\end{array}\right)^{\top},$
$A_{0,x}^{(1,d)}=\left(\begin{array}[]{c}\mathcal{X}_{0,0}^{(1,d)}\end{array}\right),\qquad
B_{0,x}^{(1,d)}=\left(\begin{array}[]{c
c}b^{0,x}_{1,1}&b^{0,x}_{2,1}\end{array}\right)^{\top},$

and

$\displaystyle
Y_{1,1}^{(1,d)}=\frac{1}{b_{1,1}^{0,y}}\left(\mathcal{Y}_{1,0}^{(1,d)}P_{1,0}^{(1,d)}+\cdots+\mathcal{Y}_{d,0}^{(1,d)}P_{d,0}^{(1,d)}\right)=\mathbf{P}^{(1,d)}\mathbf{C}_{1}^{(1,d)},$

where

$\displaystyle\left(b_{1,1}^{0,y}\right)^{2}=\left(\mathcal{Y}_{1,0}^{(1,d)}\right)^{2}+\cdots+\left(\mathcal{Y}_{d,0}^{(1,d)}\right)^{2},\>b_{1,1}^{0,x}={\langle}xY_{0,1}^{(1,d)},Y_{1,1}^{(1,d)}{\rangle}_{{\gamma}_{1,d},w}=\frac{\mathcal{X}^{(1,d)}_{1,0}\mathcal{Y}^{(1,d)}_{1,0}}{b_{1,1}^{0,y}},$

and there are (generally nonzero) constants
$C^{(1,d)}_{1,2},\ldots,C^{(1,d)}_{d,2}$ such that

$\displaystyle Y_{1,2}^{(1,d)}$
$\displaystyle=\frac{1}{b_{2,1}^{0,x}}\left(xY_{0,1}^{(1,d)}-\mathcal{X}_{0,0}^{(1,d)}Y_{0,1}^{(1,d)}-b_{1,1}^{0,x}Y_{1,1}^{(1,d)}\right)$
$\displaystyle=C^{(1,d)}_{1,2}P_{1,0}^{(1,d)}+\cdots+C^{(1,d)}_{d,2}P_{d,0}^{(1,d)}=\mathbf{P}^{(1,d)}\mathbf{C}_{2}^{(1,d)}.$

Hence (cf. (5.31) and (5.33)), ${\mathbb{C}}^{(1,d)}_{0}=\bm{e}_{0}$ and
${\mathbb{C}}^{(1,d)}_{1}=\left(\begin{array}[]{c
c}\mathbf{C}_{1}^{(1,d)}&\mathbf{C}_{2}^{(1,d)}\end{array}\right)$, where

$\mathbf{C}_{1}^{(1,d)}=\frac{1}{b_{1,1}^{0,y}}\left(\mathcal{Y}_{1,0}^{(1,d)}\bm{e}_{1}+\cdots+\mathcal{Y}_{d,0}^{(1,d)}\bm{e}_{d}\right),\qquad\mathbf{C}_{2}^{(1,d)}=C^{(1,d)}_{1,2}\bm{e}_{1}+\cdots+C^{(1,d)}_{d,2}\bm{e}_{d}.$

This illustrates that the block ${\mathbb{C}}^{(1,d)}_{1}$ in (5.31), obtained
by first multiplying by $x$ and orthogonalising, is sparser (has fewer nonzero
elements) than the one obtained by first multiplying by $y$. The same
conclusion holds true for the case $m=2$.

For ease of reference, we restate the monomial basis of degree $n$ (given in
(3.4) and (3.6)) on the curves considered in this section:

(5.39)
$\displaystyle{\mathcal{B}}_{n}({\gamma}_{m,d})=\left\\{x^{n-k}y^{k}\right\\}_{k=0}^{n},\quad
0\leq n\leq
d-1,\quad{\mathcal{B}}_{n}({\gamma}_{m,d})=\left\\{x^{d-k}y^{n+k-d}\right\\}_{k=1}^{d},\quad
n\geq d.$

To ensure the OPs span the same space as the monomial basis and in light of
the example above, to construct ${\mathbb{Y}}_{n+1}^{(m,d)}$ for $n\leq d-2$,
we first multiply by $x$ (and orthogonalise) and then multiply by $y$. That
is, $Y_{n+1,k}^{(m,d)}$, $k=1,\ldots,n+1$ are generated by orthogonalising
$xY_{n,k}^{(m,d)}$ and $Y^{(m,d)}_{n+1,n+2}$ is found by orthogonalising
$yY^{(m,d)}_{n,n+1}$. We refer to this sequence as the _primary sequence of
orthogonalisation_. This implies that for $0\leq n\leq d-1$, the expansions of
the OPs in the monomial basis (5.39) take the form

(5.40)
$Y^{(m,d)}_{n,k}={\alpha}_{k,k}^{(n)}x^{n-k+1}y^{k-1}+\cdots+{\alpha}_{1,k}^{(n)}x^{n-1}y+{\alpha}_{0,k}^{(n)}x^{n}+\mathcal{O}(x^{\ell_{1}}y^{\ell_{2}}),\qquad
0\leq\ell_{1}+\ell_{2}<n,$

for $k=1,\ldots,n+1$, with ${\alpha}_{k,k}^{(n)}\neq 0$, and where
$\mathcal{O}(x^{\ell_{1}}y^{\ell_{2}})$ represents the remainder of the
monomial expansion, which consists of terms of degree less than $n$ (hence
$\ell_{1}+\ell_{2}<n$). Hence, for the $Y_{n,k}^{(m,d)}$ generated via the
primary sequence of orthogonalization, it follows from (5.39) that
$\deg_{{\gamma}_{m,d}}(Y_{n,k}^{(m,d)})=n$.

To generate ${\mathbb{Y}}_{d}^{(m,d)}$, we note that $Y_{d,1}^{(m,d)}$ cannot
be constructed by orthogonalising $xY_{d-1,1}^{(m,d)}$ since
$xY_{d-1,1}^{(m,d)}$ has degree less than $d$ on ${\gamma}_{m,d}$. This is
because its monomial expansion is
$xY_{d-1,1}^{(m,d)}={\alpha}_{0,1}^{(d-1)}x^{d}+\mathcal{O}(x^{\ell_{1}}y^{\ell_{2}})$
with $\ell_{1}+\ell_{2}<d$ and since $y=\phi$ (if $m=1$) or $y^{2}=\phi$ (if
$m=2$), with $\phi=c_{d}x^{d}+\cdots+c_{1}x+c_{0}$, $c_{d}\neq 0$,
$\deg_{{\gamma}_{m,d}}(x^{d})=\deg_{{\gamma}_{m,d}}(xY_{d-1,1}^{(m,d)})<d$.
However, ${\mathbb{Y}}_{d}^{(m,d)}$ (or, more generally,
${\mathbb{Y}}_{n}^{(m,d)}$) can be constructed by orthogonalising
$xY_{n-1,k}^{(m,d)}$, $k=2,\ldots,d$, and then $yY_{n-1,d}^{(m,d)}$ for $n\geq
d$ because these polynomials are linearly independent and it follows from
(5.39) that
$\deg_{{\gamma}_{m,d}}(yY_{n-1,d}^{(m,d)})=\deg_{{\gamma}_{m,d}}(xY_{n-1,k}^{(m,d)})=n$,
$k=2,\ldots,d$. Alternatively, and for the same reasons, we can orthogonalise
$yY_{n-1,k}^{(m,d)}$, $k=1,\ldots,d$ to construct ${\mathbb{Y}}_{n}^{(m,d)}$,
$n\geq d$. It can be verified (using the commutativity of multiplication by
$x$ and $y$ and applying the Gram–Schmidt procedure to the monomial basis)
that both these sequences yield the same OP basis. Hence, we refer to either
of these sequences as the _secondary sequence of orthogonalisation_. This
implies that for $n\geq d$, the monomial expansions of the OPs are

(5.41)
$Y^{(m,d)}_{n,k}={\alpha}_{k,k}^{(n)}x^{d-k}y^{n-d+k}+\cdots+{\alpha}_{2,k}^{(n)}x^{d-2}y^{n-d+2}+{\alpha}_{1,k}^{(n)}x^{d-1}y^{n-d+1}+\mathcal{O}(x^{\ell_{1}}y^{\ell_{2}}),$

for $k=1,\ldots,d$, with ${\alpha}_{k,k}^{(n)}\neq 0$, and where
$\ell_{1}+\ell_{2}<n$.

To summarize: the $Y_{n,k}^{(m,d)}$ generated via the primary and secondary
sequences of orthogonalization span the same space as the monomial basis
(5.39) on the curve ${\gamma}_{m,d}$ and
$\deg_{{\gamma}_{m,d}}(Y_{n,k}^{(m,d)})=n$. Since the $Y_{n,k}^{(m,d)}$ thus
generated are orthogonal by construction, they are (orthonormal) OPs on
${\gamma}_{m,d}$ with respect to
${\langle}\cdot,\cdot{\rangle}_{{\gamma}_{m,d},w}$. The
${\mathbb{Y}}_{n}^{(m,d)}$ constructed in this manner are in
${\mathcal{V}}_{n}$ (the space of bivariate OPs of degree $n$ with respect to
$\langle\cdot,\cdot\rangle_{{\gamma}_{m,d},w}$) but note that if
$\mathcal{C}\in{\mathbb{R}}^{d\times d}$ is an orthogonal matrix, then
${\mathbb{Y}}_{n}^{(m,d)}\mathcal{C}$ (for $n\geq d-1$) are also in
${\mathcal{V}}_{n}$ since
${\langle}{\mathbb{Y}}_{n}^{(m,d)}\mathcal{C},{\mathbb{Y}}_{n}^{(m,d)}\mathcal{C}{\rangle}_{{\gamma}_{m,d},w}=I\in{\mathbb{R}}^{d\times
d}$.

We now describe the computation of the connection coefficients and the
structure of the connection matrices (examples of which are shown in Figures 1
and 2) in more detail.

###  5.5. Recursive formulae for the connection coefficients

Let $C^{(m,d)}_{r,c}$, $r,c\geq 0$ denote the entries of the connection
matrix. We introduce the map $\ell=\ell(n,k)$, from the indices $(n,k)$ of the
OP $Y^{(m,d)}_{n,k}$ for $m=1,2$ to the column index of the quasimatrix

$\mathbf{Y}^{(m,d)}=\left(\begin{array}[]{c c c
c}Y^{(m,d)}_{0,1}&Y^{(m,d)}_{1,1}&Y^{(m,d)}_{1,2}&\cdots\end{array}\right).$

The OP $Y^{(m,d)}_{0,1}$ is in column $0$ of $\mathbf{Y}^{(m,d)}$,
$Y^{(m,d)}_{1,1}$ is in column $1$ of $\mathbf{Y}^{(m,d)}$, etc. Hence,

(5.42)
$\displaystyle\ell(n,k)=\sum_{\ell=0}^{n-1}\dim{\mathcal{V}}_{\ell}+k-1=\begin{cases}\vspace{0.25cm}\displaystyle{\frac{n(n+1)}{2}}+k-1,&0\leq
n\leq d-1,\qquad 1\leq k\leq n+1\\\
\displaystyle{dn-\frac{d(d-1)}{2}}+k-1,&n\geq d,\qquad 1\leq k\leq
d\end{cases},$

where $\dim{\mathcal{V}}_{\ell}=\min\\{\ell+1,d\\}$. Let
$\mathbf{C}^{(m,d)}_{\ell(n,k)}\in{\mathbb{R}}^{\infty\times 1}$ denote column
number $\ell(n,k)\geq 0$ of the connection matrix $C^{(m,d)}$, then since
$\mathbf{Y}^{(m,d)}=\mathbf{P}^{(m,d)}C^{(m,d)}$, we have that

(5.43) $Y_{n,k}^{(m,d)}=\mathbf{P}^{(m,d)}\mathbf{C}^{(m,d)}_{\ell(n,k)},$

which expresses (3.9)–(3.10) in quasimatrix notation. For $m=1$, it follows
from (3.9) and the orthogonality of the $P^{(1,d)}_{k,0}$ polynomials with
respect to ${\langle}\cdot,\cdot{\rangle}_{{\gamma}_{1,d},w}$ that

(5.44)
$C^{(1,d)}_{r,\ell(n,k)}={\langle}Y^{(1,d)}_{n,k},P^{(1,d)}_{r,0}{\rangle}_{{\gamma}_{1,d},w},\qquad
r\geq 0,$

and similarly, for $m=2$, (3.10) and the orthogonality of the
$P^{(2,d)}_{k,0}$, $k\geq 0$ and $P^{(2,d)}_{k,1}$, $k\geq 1$ polynomials
imply that
$C^{(2,d)}_{0,0}={\langle}Y^{(2,d)}_{0,1},P^{(2,d)}_{0,0}{\rangle}_{{\gamma}_{2,d},w}$
and

(5.45)
$C^{(2,d)}_{2r-1,\ell(n,k)}={\langle}Y^{(2,d)}_{n,k},P^{(2,d)}_{r,0}{\rangle}_{{\gamma}_{2,d},w},\qquad
C^{(2,d)}_{2r,\ell(n,k)}={\langle}Y^{(2,d)}_{n,k},P^{(2,d)}_{r,1}{\rangle}_{{\gamma}_{2,d},w},\qquad
r\geq 1.$

Let $a_{i,j}^{n,x}$ for $i,j\geq 1$ denote an entry of $A_{n,x}^{(1,d)}$ or
$A_{n,x}^{(2,d)}$, where it should be clear from the context whether $m=1$ or
$m=2$, and similarly, $b_{i,j}^{n,x}$ denotes an entry of $B_{n,x}^{(1,d)}$ or
$B_{n,x}^{(2,d)}$. Likewise, $a_{i,j}^{n,y}$ and $b_{i,j}^{n,y}$ denote
entries of, respectively, $A_{n,y}^{(1,d)}$ or $A_{n,y}^{(2,d)}$ and
$B_{n,y}^{(1,d)}$ or $B_{n,y}^{(2,d)}$.

Using the known initial columns of the connection matrices given in
Proposition 5.4, we can recursively compute entries in the subsequent columns
using the results in Propositions 5.5, A.1, 5.6 and B.1

######  Proposition 5.5.

For $m=1$, the connection coefficients satisfy the following equations:
suppose first that $0\leq n\leq d-2$, then for $k=1,\ldots,n+1$

(5.46) $\displaystyle b^{n,x}_{k,k}C_{j,\ell(n+1,k)}^{(1,d)}=$
$\displaystyle\sum_{r=j-1}^{j+1}\mathcal{X}^{(1,d)}_{r,j}C_{r,\ell(n,k)}^{(1,d)}-\sum_{r=k}^{n}b^{n-1,x}_{k,r}C_{j,\ell(n-1,r)}^{(1,d)}$
$\displaystyle-\sum_{r=1}^{n+1}a^{n,x}_{r,k}C_{j,\ell(n,r)}^{(1,d)}-\sum_{r=1}^{k-1}b^{n,x}_{r,k}C_{j,\ell(n+1,r)}^{(1,d)},$

and we obtain the column $\mathbf{C}^{(1,d)}_{\ell(n+1,n+2)}$ from

(5.47) $\displaystyle b^{n,y}_{n+2,n+1}C_{j,\ell(n+1,n+2)}^{(1,d)}=$
$\displaystyle\sum_{r=j-d}^{j+d}\mathcal{Y}^{(1,d)}_{r,j}C_{r,\ell(n,n+1)}^{(1,d)}-b^{n-1,y}_{n+1,n}C_{j,\ell(n-1,n)}^{(1,d)}$
$\displaystyle-\sum_{r=1}^{n+1}a^{n,y}_{r,n+1}C_{j,\ell(n,r)}^{(1,d)}-\sum_{r=1}^{n+1}b^{n,y}_{r,n+1}C_{j,\ell(n+1,r)}^{(1,d)};$

the entries of $\mathbf{C}^{(1,d)}_{\ell(d,k)}$ for $k=1,\ldots,d$ follow from

(5.48) $\displaystyle b^{d-1,y}_{k,k}C_{j,\ell(d,k)}^{(1,d)}=$
$\displaystyle\sum_{r=j-d}^{j+d}\mathcal{Y}^{(1,d)}_{r,j}C_{r,\ell(d-1,k)}^{(1,d)}-\sum_{r=\max\\{k-1,1\\}}^{d-1}b^{d-2,y}_{k,r}C_{j,\ell(d-2,r)}^{(1,d)}$
$\displaystyle-\sum_{r=1}^{d}a^{d-1,y}_{r,k}C_{j,\ell(d-1,r)}^{(1,d)}-\sum_{r=1}^{k-1}b^{d-1,y}_{r,k}C_{j,\ell(d,r)}^{(1,d)}$

and for $n\geq d$, we compute $\mathbf{C}^{(1,d)}_{\ell(n+1,k)}$ for
$k=1,\ldots,d$ using

(5.49) $\displaystyle b^{n,y}_{k,k}C_{j,\ell(n+1,k)}^{(1,d)}=$
$\displaystyle\sum_{r=j-d}^{j+d}\mathcal{Y}^{(1,d)}_{r,j}C_{r,\ell(n,k)}^{(1,d)}-\sum_{r=k}^{d}b^{n-1,y}_{k,r}C_{j,\ell(n-1,r)}^{(1,d)}$
(5.50)
$\displaystyle-\sum_{r=1}^{d}a^{n,y}_{r,k}C_{j,\ell(n,r)}^{(1,d)}-\sum_{r=1}^{k-1}b^{n,y}_{r,k}C_{j,\ell(n+1,r)}^{(1,d)}.$

The coefficients $b^{n,x}_{k,k}$, $b^{n,y}_{n+2,n+1}$, $b^{d-1,y}_{k,k}$ and
$b^{n,y}_{k,k}$ on the left-hand sides of, respectively, (5.46)–(5.49) are
determined by the requirement that the OPs have unit norm, e.g., in (5.46) we
require
$\left(\mathbf{C}^{(1,d)}_{\ell(n+1,k)}\right)^{\top}\mathbf{C}^{(1,d)}_{\ell(n+1,k)}=1={\langle}Y_{n+1,k}^{(1,d)},Y_{n+1,k}^{(1,d)}{\rangle}_{{\gamma}_{1,d},w}$.

###### Proof.

It follows from the definition (5.24) that $A_{n,x}^{(1,d)}$ and
$A_{n,y}^{(1,d)}$ are dense symmetric matrices, e.g.,

(5.56) $\displaystyle A_{n,x}^{(1,d)}=\left(\begin{array}[]{c c c c
c}a_{1,1}^{n,x}&a_{2,1}^{n,x}&a_{3,1}^{n,x}&\cdots&a_{\dim{\mathcal{V}}_{n},1}^{n,x}\\\
a_{2,1}^{n,x}&a_{2,2}^{n,x}&a_{3,2}^{n,x}&\cdots&a_{\dim{\mathcal{V}}_{n},2}^{n,x}\\\
a_{3,1}^{n,x}&a_{3,2}^{n,x}&a_{3,3}^{n,x}&\cdots&a_{\dim{\mathcal{V}}_{n},3}^{n,x}\\\
\vdots&\ddots&\ddots&\ddots&\vdots\\\
a_{\dim{\mathcal{V}}_{n},1}^{n,x}&a_{\dim{\mathcal{V}}_{n},2}^{n,x}&\cdots&\cdots&a_{\dim{\mathcal{V}}_{n},\dim{\mathcal{V}}_{n}}^{n,x}\end{array}\right)\in{\mathbb{R}}^{\dim{\mathcal{V}}_{n}\times\dim{\mathcal{V}}_{n}},$

where $\dim{\mathcal{V}}_{n}=\dim{\mathcal{V}}_{n}({\gamma}_{1,d},w)=n+1$ for
$0\leq n\leq d-1$ and $\dim{\mathcal{V}}_{n}=d$ for $n\geq d$. The three-term
recurrences (5.22)–(5.23) and the monomial expansion (5.40) imply that for
$0\leq n\leq d-2$:

$xY_{n,k}^{(1,d)}\in\mathrm{span}\\{{\mathbb{Y}}_{n-1}^{(1,d)},{\mathbb{Y}}_{n}^{(1,d)},Y_{n+1,1}^{(1,d)},\ldots,Y_{n+1,k}^{(1,d)}\\},\qquad
1\leq k\leq n+1,$

and

$yY_{n,k}^{(1,d)}\in\mathrm{span}\\{{\mathbb{Y}}_{n-1}^{(1,d)},{\mathbb{Y}}_{n}^{(1,d)},Y_{n+1,1}^{(1,d)},\ldots,Y_{n+1,k+1}^{(1,d)}\\},\qquad
1\leq k\leq n+1.$

Therefore, the matrices $B_{n,x}^{(1,d)}$ and $B_{n,y}^{(1,d)}$ have the
following structures for $0\leq n\leq d-2$:

(5.62) $\displaystyle B_{n,x}^{(1,d)}=\left(\begin{array}[]{c c c
c}b_{1,1}^{n,x}&b_{1,2}^{n,x}&\cdots&b_{1,n+1}^{n,x}\\\
0&b_{2,2}^{n,x}&\cdots&b_{2,n+1}^{n,x}\\\ \vdots&\ddots&\ddots&\vdots\\\
0&\cdots&0&b_{n+1,n+1}^{n,x}\\\
0&\cdots&0&0\end{array}\right)\in{\mathbb{R}}^{n+2\times n+1},$

and

(5.68) $\displaystyle B_{n,y}^{(1,d)}=\left(\begin{array}[]{c c c
c}b_{1,1}^{n,y}&b_{1,2}^{n,y}&\cdots&b_{1,n+1}^{n,y}\\\
b_{2,1}^{n,y}&b_{2,2}^{n,y}&\cdots&b_{2,n+1}^{n,y}\\\
0&b_{3,2}^{n,y}&\cdots&b_{3,n+1}^{n,y}\\\ \vdots&\ddots&\ddots&\vdots\\\
0&\cdots&0&b_{n+2,n+1}^{n,y}\end{array}\right)\in{\mathbb{R}}^{n+2\times
n+1}.$

The monomial expansions (5.41) imply that for $n\geq d-1$,

$xY_{n,k}^{(1,d)}\in\mathrm{span}\\{{\mathbb{Y}}_{n-1}^{(1,d)},{\mathbb{Y}}_{n}^{(1,d)},Y_{n+1,1}^{(1,d)},\ldots,Y_{n+1,k-1}^{(1,d)}\\},\qquad
1\leq k\leq d,$

and

$yY_{n,k}^{(1,d)}\in\mathrm{span}\\{{\mathbb{Y}}_{n-1}^{(1,d)},{\mathbb{Y}}_{n}^{(1,d)},Y_{n+1,1}^{(1,d)},\ldots,Y_{n+1,k}^{(1,d)}\\},\qquad
1\leq k\leq d.$

Hence, for $n\geq d-1$,

(5.74) $\displaystyle B_{n,x}^{(1,d)}=\left(\begin{array}[]{c c c c
c}0&b_{1,2}^{n,x}&b_{1,3}^{n,x}&\cdots&b_{1,d}^{n,x}\\\
0&0&b_{2,3}^{n,x}&\cdots&b_{2,d}^{n,x}\\\
\vdots&\ddots&\ddots&\ddots&\vdots\\\ 0&\cdots&0&0&b_{d-1,d}^{n,x}\\\
0&\cdots&\cdots&0&0\end{array}\right)\in{\mathbb{R}}^{d\times d}$

and

(5.79) $\displaystyle B_{n,y}^{(1,d)}=\left(\begin{array}[]{c c c c
}b_{1,1}^{n,y}&b_{1,2}^{n,y}&\cdots&b_{1,d}^{n,y}\\\
0&b_{2,2}^{n,y}&\cdots&b_{2,d}^{n,y}\\\ \vdots&\ddots&\ddots&\vdots\\\
0&\cdots&0&b_{d,d}^{n,y}\end{array}\right)\in{\mathbb{R}}^{d\times d}.$

The three-term recurrences (5.22) and (5.23) satisfied by the OPs and the
structures of the matrices given above imply that for $0\leq n\leq d-2$ and
$k=1,\ldots,n+1$

(5.80) $\displaystyle xY_{n,k}^{(1,d)}$
$\displaystyle=\sum_{r=k}^{n}b^{n-1,x}_{k,r}Y_{n-1,r}^{(1,d)}+\sum_{r=1}^{n+1}a^{n,x}_{r,k}Y_{n,r}^{(1,d)}+\sum_{r=1}^{k}b^{n,x}_{r,k}Y_{n+1,r}^{(1,d)}$
(5.81) $\displaystyle yY_{n,k}^{(1,d)}$
$\displaystyle=\sum_{r=\max\\{k-1,1\\}}^{n}b^{n-1,y}_{k,r}Y_{n-1,r}^{(1,d)}+\sum_{r=1}^{n+1}a^{n,y}_{r,k}Y_{n,r}^{(1,d)}+\sum_{r=1}^{k+1}b^{n,y}_{r,k}Y_{n+1,r}^{(1,d)}.$

The recurrences satisfied by the OPs for the cases $n=d-1$, $k=1,\ldots,d$ and
$n\geq d$, $k=1,\ldots,d$ can be found similarly using (5.56)–(5.79).

Let $\mathcal{X}^{(1,d)}_{r,c}$ and $\mathcal{Y}^{(1,d)}_{r,c}$, $r,c\geq 0$
denote the entries of $\mathcal{X}^{(1,d)}$ and $\mathcal{Y}^{(1,d)}$,
respectively, which are defined in (5.12). Since $\mathcal{X}^{(1,d)}$ is
tridiagonal and $\mathcal{Y}^{(1,d)}$ has bandwidths $(d,d)$, we have for
$j\geq 1$ that

(5.82) $\displaystyle
xP_{j,0}^{(1,d)}=\sum_{r=j-1}^{j+1}\mathcal{X}^{(1,d)}_{r,j}P_{r,0}^{(1,d)},\qquad
yP_{j,0}^{(1,d)}=\sum_{r=j-d}^{j+d}\mathcal{Y}^{(1,d)}_{r,j}P_{r,0}^{(1,d)}.$

It now follows from (5.44), (5.80) and (5.82) that for $0\leq n\leq d-2$ and
$k=1,\ldots,n+1$

$\displaystyle{\langle}xY^{(1,d)}_{n,k},P^{(1,d)}_{j,0}{\rangle}_{{\gamma}_{1,d},w}=\sum_{r=k}^{n}b^{n-1,x}_{k,r}C_{j,\ell(n-1,r)}^{(1,d)}+\sum_{r=1}^{n+1}a^{n,x}_{r,k}C_{j,\ell(n,r)}^{(1,d)}+\sum_{r=1}^{k}b^{n,x}_{r,k}C_{j,\ell(n+1,r)}^{(1,d)}$
$\displaystyle={\langle}Y^{(1,d)}_{n,k},xP^{(1,d)}_{j,0}{\rangle}_{{\gamma}_{1,d},w}=\sum_{r=j-1}^{j+1}\mathcal{X}^{(1,d)}_{r,j}C_{r,\ell(n,k)}^{(1,d)},$

from which we obtain (5.46). The recurrence (5.47) can similarly be derived
using (5.44), (5.81) and (5.82). The remaining equations (5.48)–(5.49) follow
from the recurrences satisfied by the OPs for the cases $n=d-1$,
$k=1,\ldots,d$ and $n\geq d$, $k=1,\ldots,d$. ∎

The formulae from which the connection coefficients can be computed for $m=2$
are given in Appendix A.

###  5.6. Structure of the connection matrices

We now aim to deduce the entries of the connection matrices that are zero,
which will specify the ranges over which the row indices should vary in the
formulae for the connection coefficients for $m=1$ (given in (5.46)–(5.49))
and for $m=2$ (given in Appendix A).

######  Proposition 5.6.

For $m=1$, if $0\leq n\leq d-2$, then $C^{(1,d)}_{j,\ell(n+1,k)}=0$ for

(5.83) $\displaystyle j<n+1\qquad$ $\displaystyle\text{and}\qquad 1\leq k\leq
n+2$ (5.84) $\displaystyle j>nd\qquad$ $\displaystyle\text{and}\qquad 1\leq
k\leq n$ (5.85) $\displaystyle j>nd+1\qquad$ $\displaystyle\text{and}\qquad
k=n+1$ (5.86) $\displaystyle j>(n+1)d\qquad$ $\displaystyle and\qquad k=n+2.$

If $n\geq d-1$, then $C^{(1,d)}_{j,\ell(n+1,k)}=0$ for

(5.87) $\displaystyle j<d-1+d(n+2-d)\qquad$ $\displaystyle and\qquad 1\leq
k\leq d$ (5.88) $\displaystyle j>nd\qquad$ $\displaystyle and\qquad 1\leq
k\leq d-2$ (5.89) $\displaystyle j>nd+1\qquad$ $\displaystyle and\qquad k=d-1$
(5.90) $\displaystyle j>(n+1)d\qquad$ $\displaystyle and\qquad k=d$

###### Proof.

Recall that the (maximal) degree of an OP $Y^{(1,d)}_{n,k}$ on
${\gamma}_{1,d}$ is obtained by setting $y=\phi(x)$ and finding the degree of
the resulting polynomial, which is also the degree of $Y^{(1,d)}_{n,k}$ in the
$P^{(1,d)}_{n,0}$ basis. We denote this maximal degree by $\deg
Y^{(1,d)}_{n,k}$, which is in general different from the (minimal) degree of
$Y^{(1,d)}_{n,k}$ on ${\gamma}_{1,d}$, viz. $n$, which is denoted by
$\deg_{{\gamma}_{1,d}}(Y^{(1,d)}_{n,k})$.

To determine the degrees of the OPs $Y^{(1,d)}_{n+1,k}$ in the
$P^{(1,d)}_{j,0}$ basis, we let $n\mapsto n+1$ in (5.40) and (5.41) and note
that among the degree $n+1$ monomial terms, the leading term has the largest
maximal degree, which is, for $0\leq n\leq d-2$,

(5.91) $\deg x^{n+2-k}y^{k-1}=n+2-k+d(k-1),\qquad k=1,\ldots,n+2,$

and for $n\geq d-1$,

(5.92) $\deg x^{d-k}y^{n+1+k-d}=d-k+d(n+1+k-d),\qquad k=1,\ldots,d.$

All the OPs of degree $n$ (${\mathbb{Y}}_{n}^{(1,d)}$) appear in the
recurrences satisfied by the OPs of degree $n+1$ (see (5.22)–(5.23)). Among
all monomials of degree $n$ on ${\gamma}_{1,d}$, $y^{n}$ has the largest
maximal degree, with $\deg y^{n}=nd$. We conclude from this and (5.91) that
for $0\leq n\leq d-2$, $\deg Y^{(1,d)}_{n+1,k}=\max\\{nd,n+2-k+d(k-1)\\}$ and
from (5.92) it follows that for $n\geq d-1$, $\deg
Y^{(1,d)}_{n+1,k}=\max\\{nd,d-k+d(n+1+k-d)\\}$. Hence,

$\displaystyle\deg Y^{(1,d)}_{n+1,k}=\begin{cases}nd&\text{for }0\leq n\leq
d-2,1\leq k\leq n\text{ and }n\geq d-1,1\leq k\leq d-2\\\ nd+1&\text{for
}0\leq n\leq d-2,k=n+1\text{ and }n\geq d-1,k=d-1\\\ (n+1)d&\text{for }0\leq
n\leq d-2,k=n+2\text{ and }n\geq d-1,k=d\end{cases}.$

Since these maximal degrees are also the degrees of the $Y^{(1,d)}_{n+1,k}$ in
the $P^{(1,d)}_{j,0}$ basis, they are the upper bounds on $j$ in the
expansions $Y^{(1,d)}_{n+1,k}=\sum_{j\geq
0}C^{(1,d)}_{j,\ell(n+1,k)}P^{(1,d)}_{j,0}$. Thus, we have shown that
$C^{(1,d)}_{j,\ell(n+1,k)}=0$ for the cases (5.84)–(5.86) and (5.88)–(5.90).

Note from (5.91)–(5.92) that the smallest maximal degree of a monomial of
degree $n+1$ on ${\gamma}_{1,d}$ is attained for $k=1$, which is $n+1$ for
$0\leq n\leq d-2$ and $d-1+d(n+2-d)$ for $n\geq d-1$. This implies that

(5.93) $\deg_{{\gamma}_{1,d}}(P^{(1,d)}_{j,0})<n+1\text{ for
}\begin{cases}j<n+1&\text{if }0\leq n\leq d-2\\\ j<d-1+d(n+2-d)&\text{if
}n\geq d-1\end{cases}.$

Since the $Y^{(1,d)}_{n+1,k}$ are OPs of degree $n+1$ on ${\gamma}_{1,d}$, it
follows that

(5.94)
$C^{(1,d)}_{j,\ell(n+1,k)}={\langle}Y^{(1,d)}_{n+1,k},P^{(1,d)}_{j,0}{\rangle}_{{\gamma}_{1,d},w}=0\quad\text{if}\quad\deg_{{\gamma}_{1,d}}(P^{(1,d)}_{j,0})<n+1.$

Hence, combining (5.93) and (5.94), we have shown that
$C^{(1,d)}_{j,\ell(n+1,k)}=0$ for the cases (5.83) and (5.87). ∎

For $m=2$, the indices for which $C^{(2,d)}_{j,\ell(n+1,k)}=0$ are given in
Appendix B.

###  5.7. Numerical stability and computational complexity

The recursive formulae for the connection coefficients given in section 5.5
are numerically unstable. This is to be expected since the formulae arise from
the Gram–Schmidt procedure, which is well-known to be unstable [1] and the
instability of the Lanczos algorithm is also well-established (see [7, section
2.2.3] and [3, Supplement D]). We shall consider the computational cost of our
Lanczos (or, equivalently, Gram–Schmidt) algorithm when it is stabilized via
reorthogonalization [1, 3] (the Lanczos algorithm can also be stabilized if it
is performed via Householder transformations [20] or Givens rotations [8]).

Due to the instability of the Lanczos procedure, the OPs are not exactly
orthogonal in floating point arithmetic. After a new OP has been generated via
a new column of the connection matrix (using the formulae in section 5.5 and
Appendix A), we orthogonalize this OP against lower degree OPs to ensure
orthogonality up to the level of roundoff errors. That is, we replace the OP
$Y_{n+1,k}^{(m,d)}$ with its re-orthogonalized version:

(5.95) $Y_{n+1,k}^{(m,d)}\leftarrow
Y_{n+1,k}^{(m,d)}-\sum_{\begin{subarray}{c}0\leq s\leq n+1\\\ 1\leq
i\leq\dim{\mathcal{V}}_{s}\end{subarray}}\langle
Y_{n+1,k}^{(m,d)},Y^{(m,d)}_{s,i}\rangle_{{\gamma}_{m,d},w}\,Y^{(m,d)}_{s,i},$

which we express in terms of the columns of the connection matrix by setting
$Y_{n,k}^{(m,d)}=\mathbf{P}^{(m,d)}\mathbf{C}^{(m,d)}_{\ell(n,k)}$:

$\mathbf{C}^{(m,d)}_{\ell(n+1,k)}\leftarrow\mathbf{C}^{(m,d)}_{\ell(n+1,k)}-\sum_{\begin{subarray}{c}0\leq
s\leq n+1\\\ 1\leq
i\leq\dim{\mathcal{V}}_{s}\end{subarray}}\left[\left(\mathbf{C}^{(m,d)}_{\ell(n+1,k)}\right)^{\top}\mathbf{C}^{(m,d)}_{\ell(s,i)}\right]\mathbf{C}^{(m,d)}_{\ell(s,i)}.$

We can use the results in section 5.6 and Appendix B (i.e., the fact that the
connection matrices are banded) to show that we only need to re-orthogonalize
against a constant number of lower degree OPs (as opposed to _all_ lower
degree OPs). For example, for $m=1$, we have

$Y^{(1,d)}_{n+1,k}=\sum_{j=j_{\min}}^{j_{\max}}C^{(1,d)}_{j,\ell(n+1,k)}P^{(1,d)}_{j,0},$

where $j_{\min}=j_{\min}(n+1,k)$ and $j_{\max}=j_{\max}(n+1,k)$ are given in
Proposition 5.6. For instance, for $n\geq d-1$ and $1\leq k\leq d$,
$j_{\min}(n+1,k)=d-1+d(n+2-d)$ and for $k=d$, $j_{\max}(n+1,d)=(n+1)d$.
Likewise, for $m=2$, in (A.1), $j$ ranges from $j_{\min}$ to $j_{\max}$, which
are given in Appendix B. Hence, if in (5.95) we have
$j_{\min}(n+1,k)>j_{\max}(s,i)$, then $\langle
Y_{n+1,k}^{(m,d)},Y^{(m,d)}_{s,i}\rangle_{{\gamma}_{m,d},w}$ is identically
zero, also in floating point arithmetic, and therefore we only need to re-
orthogonalize against OPs $Y^{(m,d)}_{s,i}$ with $j_{\min}(n+1,k)\leq
j_{\max}(s,i)$. For $m=1$, we have to re-orthogonalise $Y^{(m,d)}_{n+1,k}$
against $Y^{(m,d)}_{s,i}$ with $n-d+3\leq s\leq n+1$ because

(5.96)
$\displaystyle\begin{split}&\min\left\\{s\in{\mathbb{N}}_{0}\>:j_{\min}(n+1,k)\leq
j_{\max}(s,i)\right\\}\\\
&=\min\left\\{s\in{\mathbb{N}}_{0}\>:\>d-1+d(n+2-d)\leq
sd\right\\}=n-d+3\end{split}$

and for $m=2$, $n-d+4\leq s\leq n+1$ since

(5.97)
$\displaystyle\begin{split}&\min\left\\{s\in{\mathbb{N}}_{0}\>:j_{\min}(n+1,k)\leq
j_{\max}(s,i)\right\\}\\\
&=\min\left\\{s\in{\mathbb{N}}_{0}\>:\>2d+d(n+1-d)\leq
sd-1\right\\}=n-d+4.\end{split}$

######  Proposition 5.7.

For $m=1,2$, the connection coefficients and the entries of the Jacobi
matrices for OPs of degrees $0,\ldots,N$ (i.e., ${\mathbb{C}}^{(m,d)}_{n}$,
$A_{n,x}^{(m,d)}$, $A_{n,y}^{(m,d)}$ for $0\leq n\leq N$ and
$B_{n,x}^{(m,d)}$, $B_{n,y}^{(m,d)}$ for $0\leq n\leq N-1$) can be computed in
$\mathcal{O}(Nd^{4})$ operations without re-orthogonalisation and in
$\mathcal{O}(Nd^{5})$ operations with re-orthogonalisation.

###### Proof.

Starting with the first few columns of $C^{(m,d)}$ given in Proposition 5.4,
the subsequent columns of the connection matrix are computed using (for $m=1$)
(5.46)–(5.49), with $j$ ranging between the bounds given in (5.83)–(5.90) and
(for $m=2$) using the formulae in Appendix A, with the row indices ranging
between the bounds given in Appendix B.

To compute the entries of $A_{n,x}^{(m,d)}$, $B_{n,x}^{(m,d)}$,
$A_{n,y}^{(m,d)}$ and $B_{n,y}^{(m,d)}$, we use the fact that

(5.98)
$\displaystyle{\langle}xY_{n_{1},i_{1}}^{(m,d)},Y_{n_{2},i_{2}}^{(m,d)}{\rangle}_{{\gamma}_{m,d},w}=\left(\mathbf{C}^{(m,d)}_{\ell(n_{1},i_{1})}\right)^{\top}\mathcal{X}^{(m,d)}\mathbf{C}^{(m,d)}_{\ell(n_{2},i_{2})}=\left(\mathbf{C}^{(m,d)}_{\ell(n_{2},i_{2})}\right)^{\top}\mathcal{X}^{(m,d)}\mathbf{C}^{(m,d)}_{\ell(n_{1},i_{1})},$

which follows from (5.43), the orthonormality of the $P$-polynomials and the
symmetry of $\mathcal{X}^{(m,d)}$. A similar identity holds for inner products
of the form
${\langle}yY_{n_{1},i_{1}}^{(m,d)},Y_{n_{2},i_{2}}^{(m,d)}{\rangle}_{{\gamma}_{m,d},w}$.

The most computationally expensive part of the procedure is the computation of
$\mathcal{Y}^{(m,d)}\mathbf{C}^{(m,d)}_{\ell(n,k)}$. For $m=1$, since the
indices of the nonzero rows of $\mathbf{C}^{(1,d)}_{\ell(n,k)}$ range from
$d-1+d(n-d+1)$ to $nd$, see (5.87) and (5.90) with $n+1\mapsto n$, the column
$\mathbf{C}^{(1,d)}_{\ell(n,k)}$ has $\mathcal{O}(d^{2})$ nonzero entries.
Since $\mathcal{Y}^{(1,d)})$ has bandwidths $(d,d)$, the computation of
$\mathcal{Y}^{(1,d)}\mathbf{C}^{(1,d)}_{\ell(n,k)}$ requires
$\mathcal{O}(d^{3})$ operations, which has to be computed $\mathcal{O}(Nd)$
times, hence the computational complexity is $\mathcal{O}(Nd^{4})$ without re-
orthogonalisation. By the same reasoning, for $m=2$ the computational
complexity is also $\mathcal{O}(Nd^{4})$, however the computational cost is
reduced by approximately a factor of $4$ since $\mathcal{Y}^{(2,d)}$ and
$\mathbf{C}^{(2,d)}_{\ell(n,k)}$ have roughly half the number of nonzero
entries compared to $\mathcal{Y}^{(1,d)}$ and
$\mathbf{C}^{(1,d)}_{\ell(n,k)}$.

The re-orthogonalisation step (see (5.96)–(5.97)) costs $\mathcal{O}(d^{4})$
operations per OP for $m=1,2$ and hence the overall computational cost of our
Lanczos procedure with re-orthogonalisation is $\mathcal{O}(Nd^{5})$. ∎

The linear computational complexity (in $N$) of the Lanczos method applied to
an orthogonal basis compares favourably to the cubic $\mathcal{O}(N^{3}d^{3})$
complexity (without re-orthogonalisation) of orthogonalising the bivariate
monomial bases111One orthogonalisation method is to compute the Cholesky
factor $L$ of the Gram matrix $\mathcal{M}$, where
$\mathcal{M}=\langle\mathbf{M},\mathbf{M}\rangle_{{\gamma}_{m,d},w}=LL^{\top}$,
and where $\mathbf{M}$ is a quasimatrix whose columns are the monomial basis
functions in (5.39). Then the OPs are given by
$\mathbf{Y}^{(m,d)}=\mathbf{M}L^{-\top}$. For OPs of degrees $0,\ldots,N$ this
has $\mathcal{O}(N^{3}d^{3})$ complexity since it requires the Cholesky
decomposition of a $\nu\times\nu$ finite section of the dense matrix
$\mathcal{M}$ with $\nu=\mathcal{O}(Nd)$.  in (5.39).

We remark that if the results of Propositions 5.5, A.1, 5.6 and B.1 are used
to compute the connection coefficients for the cases in Proposition 4.1, then
the same OP bases as in (4.1)–(4.41) are obtained.

###  5.8. Numerical experiments

For $m=1$, if we compare the upper and lower bounds on $j$ in Proposition 5.6
to the connection matrices shown in Figure 1, we find that the bounds for the
cases (5.84)–(5.86) and (5.88)–(5.90) are sharp, while the bounds for the
cases (5.83) and (5.87) are sharp only for $k=1$. That is, the entries of the
first three connection matrices in Figure 1 are nonzero for the following rows
and columns: for $0\leq n\leq d-2$, $C^{(1,d)}_{j,\ell(n+1,k)}\neq 0$ for

$\displaystyle n+1\leq j\leq nd\qquad$ $\displaystyle\text{and}\qquad k=1$
$\displaystyle n+2\leq j\leq nd\qquad$ $\displaystyle\text{and}\qquad 2\leq
k\leq n$ $\displaystyle n+2\leq j\leq nd+1\qquad$
$\displaystyle\text{and}\qquad k=n+1$ $\displaystyle n+2\leq
j\leq(n+1)d\qquad$ $\displaystyle\text{and}\qquad k=n+2$

and for $n\geq d-1$, $C^{(1,d)}_{j,\ell(n+1,k)}\neq 0$ for

(5.99) $\displaystyle d-1+d(n+2-d)\leq j\leq nd\qquad$
$\displaystyle\text{and}\qquad k=1$ $\displaystyle 2d-2+d(n+2-d)\leq j\leq
nd\qquad$ $\displaystyle\text{and}\qquad 2\leq k\leq d-2$ $\displaystyle
2d-1+d(n+2-d)\leq j\leq nd+1\qquad$ $\displaystyle\text{and}\qquad k=d-1$
(5.100) $\displaystyle 2d-1+d(n+2-d)\leq j\leq(n+1)d\qquad$
$\displaystyle\text{and}\qquad k=d$

These bounds on the indices of the nonzero entries of $C^{(1,d)}$ imply that
the maximum bandwidths of $C^{(1,d)}$ are attained at the entries
$C^{(1,d)}_{d-1+d(n+2-d),\ell(n+1,1)}$ (see (5.99)) and
$C^{(1,d)}_{(n+1)d,\ell(n+1,d)}$ (see (5.100)). Therefore the bandwidths of
$C^{(1,d)}$ are $(\lambda,\mu)$ with $\lambda=\mu=\lambda_{1}$, where

$\lambda_{1}=\frac{(d-2)(d-1)}{2}$

because

$\lambda=\ell(n+1,1)-\left[d-1+d(n+2-d)\right]=\frac{(d-2)(d-1)}{2},$

where $\ell$ is defined in (5.42), and

$\mu=(n+1)d-\ell(n+1,d)=\frac{(d-2)(d-1)}{2}.$

Similar to the case $m=1$, for $m=2$ we find that in the connection matrices
in Figure 2, the lower bounds on the row indices $j$ with
$C^{(2,d)}_{j,\ell(n+1,k)}=0$ given in Appendix B are sharp, while the upper
bounds (B.1), (B.7) and (B.8) are sharp only in the case $k=1$. The maximum
bandwidths are therefore attained at the connection coefficients corresponding
to $P^{(2,d)}_{j,0}$ with $j=(n+1)d/2$ (which is the largest among the maximal
degrees in (B.10)–(B.11), which is attained for $k=d$ and $n$ odd) and
$P^{(2,d)}_{j,1}$ with $j=d+d(n+1-d)/2$ (which is the smallest among the
maximal degrees in (B.10)–(B.11), which is attained for $k=1$, with
$\mathrm{mod}(n+1-d,2)=0$). That is, the maximum bandwidths are attained at
$C_{(n+1)d-1,\ell(n+1,d)}^{(2,d)}$ and $C_{2d+d(n+1-d),\ell(n+1,1)}^{(2,d)}$.
Hence the bandwidths of the connection matrices in Figure 2 are
$(\lambda,\mu)=(\lambda_{2},\lambda_{2})$, where

$\lambda_{2}=\frac{d(d-3)}{2}$

This is because the number of rows from the main diagonal to the former entry
is

$(n+1)d-1-\ell(n+1,d)=\frac{d(d-3)}{2},$

and for the latter entry, the number is

$\ell(n+1,1)-\left[2d+d(n+1-d)\right]=\frac{d(d-3)}{2}.$

The bottom-right frame of Figures 1 and 2 show that if $w$ and $\phi$ are even
functions, then the connection matrix has roughly half the number of nonzero
entries compared to the case where either $w$ or $\phi$ is not even. This is
because the multiplication matrices $\mathcal{X}^{(m,d)}$ and
$\mathcal{Y}^{(m,d)}$ have roughly half the number of nonzero entries for even
$w$ and $\phi$.

##  6\. Conclusion

For the cases $m=1,2$, we have constructed bivariate OPs on the algebraic
curves $y^{m}=\phi$. For $m>2$, if we generalise the $P$-polynomials

$P^{(m,d)}_{n,k}=y^{k}p_{n-k}(\phi^{2k/m}w;x),\qquad k=0,\ldots,m-1,$

and the inner product

$\displaystyle{\langle}f,g{\rangle}_{{\gamma}_{m,d},w}=\begin{cases}\displaystyle{\int}f\left(x,\sqrt[m]{\phi(x)}\right)g\left(x,\sqrt[m]{\phi(x)}\right)w(x)\mathrm{d}x,&m\text{
odd},\\\
\displaystyle{\int}\left[f\left(x,\sqrt[m]{\phi(x)}\right)g\left(x,\sqrt[m]{\phi(x)}\right)\right.&\\\
\left.\>\>\>+f\left(x,-\sqrt[m]{\phi(x)}\right)g\left(x,-\sqrt[m]{\phi(x)}\right)\right]w(x)\mathrm{d}x,&m\text{
even},\end{cases}$

then the $P$-polynomials are not orthogonal with respect to
${\langle}\cdot,\cdot{\rangle}_{{\gamma}_{m,d},w}$ for $m>2$. We can
nevertheless still construct the OPs on $y^{m}=\phi$ as linear combinations of
the $P$-polynomials. However, the computation of the connection matrix would
require the Cholesky factorisation of a dense Gram matrix, which has cubic
complexity, see footnote 1. For $m=1,2$, we showed that the connection matrix
is banded but for $m>2$, the connection matrix will in general be upper
triangular.

Topics for future research include the development of the theory of
interpolation, quadrature and orthogonal series on the curves $y^{m}=\phi$, as
was done for quadratic and cubic curves in [14, 6]. As an application, OPs on
$y^{m}=\phi(x)$ could be used to approximate functions with algebraic
singularities of the form $\sqrt[m]{\phi(x)}$ at the zeros of $\phi$, again as
was done in [14, 6] for square-root and cubic singularities. Another
application of the OPs on $y^{m}=\phi$ is the construction of OPs inside
$y^{m}=\phi$ and their associated surfaces of revolution for the purpose of
spectral methods for partial differential equations.

## References

  * [1] A. Björck, Numerical Methods in Matrix Computations, Springer, 2015.
  * [2] P.D. Brubeck, Y. Nakatsukasa and L.N. Trefethen, Vandermonde with Arnoldi, SIAM Review, 63(2) (2021), 405–415.
  * [3] T. Chen, T. Trogdon and S. Ubaru, Analysis of stochastic Lanczos quadrature for spectrum approximation, ICML, 2021, 1728–1739.
  * [4] Library of Mathematical Functions, http://dlmf.nist.gov/, F. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller, B. V. Saunders, H. S. Cohl, and M. A. McClain, eds.
  * [5] C. F. Dunkl and Y. Xu, Orthogonal Polynomials of Several Variables, Second edition, Encyclopedia of Mathematics and its Applications 155, Cambridge University Press, Cambridge, 2014.
  * [6] M. Fasondini, S. Olver and Y. Xu, Orthogonal polynomials on planar cubic curves. Found. Comp. Math., to appear.
  * [7] W. Gautschi, Orthogonal Polynomials: Computation and Approximation, Oxford University Press, 2004.
  * [8] W.B. Gragg and W.J. Harrod, The numerically stable reconstruction of Jacobi matrices from spectral data, Numer. Math., 44(3) (1984), 317–335.
  * [9] G.H. Golub and G. Meurant, Matrices, moments and quadrature with applications, Princeton University Press, 2009.
  * [10] E. Hendriksen, H. van Rossum, Semi-classical orthogonal polynomials, pp.354–361 in Polynômes Orthogonaux et Applications, Proceedings, Bar-le-Duc 1984, (C.Brezinski et al., editors), Lecture Notes Math. 1171, Springer, Berlin 1985\.
  * [11] S. Olver and Y. Xu, Orthogonal structure on a wedge and on the boundary of a square, Found. Comp. Math., 19 (2019), 561–589.
  * [12] S. Olver and Y. Xu, Orthogonal polynomials in and on a quadratic surface of revolution, Maths Comp., 89 (2020), 2847–2865.
  * [13] S. Olver, R.M. Slevinsky and A. Townsend, Fast algorithms using orthogonal polynomials, Acta Numer., 29 (2020), 573–699.
  * [14] S. Olver and Y. Xu, Orthogonal structure on a quadratic curve, IMA J. Numer. Anal., 41 (2021), 206–246.
  * [15] S. Olver and Y. Xu, Non-homogeneous wave equation on a cone, Int. Trans. Spec. Funcs., 32 (2021), 604–619.
  * [16] T. Pu and M. Fasondini, The numerical solution of fractional integral equations via orthogonal polynomials in fractional powers, arXiv:2206.14280, 2022\.
  * [17] B. Snowball and S. Olver, Sparse spectral and p-finite element methods for partial differential equations on disk slices and trapeziums, Stud. Appl. Math., 145(1) (2020), 3–35.
  * [18] B. Snowball and S. Olver, Sparse spectral methods for partial differential equations on spherical caps, to appear in Trans. Maths Appl.
  * [19] L.N. Trefethen, Householder triangularization of a quasimatrix, IMA J. Numer. Anal., 30(4), (2010), 887–897.
  * [20] H.F. Walker, Implementation of the GMRES method using Householder transformations, SIAM J. Sci. Comput., 9 (1) (1988), 152–163.
  * [21] Y. Xu, Orthogonal polynomials and Fourier orthogonal series on a cone, J. Fourier Anal. Appl., 26(3), 2020, 1–42.
  * [22] Y. Xu, Orthogonal structure and orthogonal series in and on a double cone or a hyperboloid, Trans. Amer. Math. Soc., 374 (2021), 3603–3657.
  * [23] Y. Xu, Laguerre expansions on conic domains, J. Fourier Anal. Appl., 27(64), 2021, 36 pp.
  * [24] Y. Xu, Approximation and localized polynomial frame on conic domains,  J. Functional Anal., 281(12), 2021, paper No. 109257, 94 pp.
  * [25] Y. Xu, Fourier orthogonal series on a paraboloid, Submitted.
  * [26] Y. Xu, Approximation and localized polynomial frame on double hyperbolic and conic domains, Submitted.

##  Appendix A Recurrence relations for the connection coefficients for $m=2$

######  Proposition A.1.

For $m=2$, the connection coefficients satisfy the following equations:
suppose that $0\leq n\leq d-2$, then for $k=1,3,\ldots,n+1$

$\displaystyle
b^{n,x}_{k,k}C^{(2,d)}_{2j-1,\ell(n+1,k)}={\beta}_{j-1}(w)C^{(2,d)}_{2j-3,\ell(n,k)}+{\alpha}_{j}(w)C^{(2,d)}_{2j-1,\ell(n,k)}+{\beta}_{j}(w)C^{(2,d)}_{2j+1,\ell(n,k)}$
$\displaystyle-\sum_{\begin{subarray}{c}r=k\\\
\mathrm{mod}(r,2)=\mathrm{mod}(k,2)\end{subarray}}^{n}b^{n-1,x}_{k,r}C^{(2,d)}_{2j-1,\ell(n-1,r)}-\sum_{\begin{subarray}{c}r>0\\\
\mathrm{mod}(r,2)=\mathrm{mod}(k,2)\end{subarray}}^{n+1}a^{n,x}_{r,k}C^{(2,d)}_{2j-1,\ell(n,r)}$
$\displaystyle-\sum_{\begin{subarray}{c}r>0\\\
\mathrm{mod}(r,2)=\mathrm{mod}(k,2)\end{subarray}}^{k-2}b^{n,x}_{r,k}C^{(2,d)}_{2j-1,\ell(n+1,r)};$

and for $k=2,4,\ldots,n+1$,

$\displaystyle b^{n,x}_{k,k}C^{(2,d)}_{2j,\ell(n+1,k)}={\beta}_{j-2}(\phi
w)C^{(2,d)}_{2j-2,\ell(n,k)}+{\alpha}_{j-1}(\phi
w)C^{(2,d)}_{2j,\ell(n,k)}+{\beta}_{j-1}(\phi w)C^{(2,d)}_{2j+2,\ell(n,k)}$
$\displaystyle-\sum_{\begin{subarray}{c}r=k\\\
\mathrm{mod}(r,2)=\mathrm{mod}(k,2)\end{subarray}}^{n}b^{n-1,x}_{k,r}C^{(2,d)}_{2j,\ell(n-1,r)}-\sum_{\begin{subarray}{c}r>0\\\
\mathrm{mod}(r,2)=\mathrm{mod}(k,2)\end{subarray}}^{n+1}a^{n,x}_{r,k}C^{(2,d)}_{2j,\ell(n,r)}$
$\displaystyle-\sum_{\begin{subarray}{c}r>0\\\
\mathrm{mod}(r,2)=\mathrm{mod}(k,2)\end{subarray}}^{k-2}b^{n,x}_{r,k}C^{(2,d)}_{2j,\ell(n+1,r)}.$

Suppose that $0\leq n\leq d-2$ and $k=n+1$, if $n$ is odd

$\displaystyle
b^{n,y}_{k+1,k}C_{2j-1,\ell(n+1,k+1)}^{(2,d)}=\sum_{r=j-d+1}^{j+1}r_{r-1,j}C_{2r,\ell(n,k)}^{(2,d)}-\sum_{\begin{subarray}{c}r=\max\\{k-1,1\\}\\\
\mathrm{mod}(r,2)=\mathrm{mod}(k-1,2)\end{subarray}}^{n}b^{n-1,y}_{k,r}C_{2j-1,\ell(n-1,r)}^{(2,d)}$
$\displaystyle-\sum_{\begin{subarray}{c}r>0\\\
\mathrm{mod}(r,2)=\mathrm{mod}(k-1,2)\end{subarray}}^{n+1}a^{n,y}_{r,k}C_{2j-1,\ell(n,r)}^{(2,d)}-\sum_{\begin{subarray}{c}r>0\\\
\mathrm{mod}(r,2)=\mathrm{mod}(k-1,2)\end{subarray}}^{k-1}b^{n,y}_{r,k}C_{2j-1,\ell(n+1,r)}^{(2,d)}$

and if $n$ is even

$\displaystyle
b^{n,y}_{k+1,k}C_{2j,\ell(n+1,k+1)}^{(2,d)}=\sum_{r=j-1}^{j-1+d}r_{j-1,r}C_{2r-1,\ell(n,k)}^{(2,d)}-\sum_{\begin{subarray}{c}r=\max\\{k-1,1\\}\\\
\mathrm{mod}(r,2)=\mathrm{mod}(k-1,2)\end{subarray}}^{n}b^{n-1,y}_{k,r}C_{2j,\ell(n-1,r)}^{(2,d)}$
$\displaystyle-\sum_{\begin{subarray}{c}r>0\\\
\mathrm{mod}(r,2)=\mathrm{mod}(k-1,2)\end{subarray}}^{n+1}a^{n,y}_{r,k}C_{2j,\ell(n,r)}^{(2,d)}-\sum_{\begin{subarray}{c}r>0\\\
\mathrm{mod}(r,2)=\mathrm{mod}(k-1,2)\end{subarray}}^{k-1}b^{n,y}_{r,k}C_{2j,\ell(n+1,r)}^{(2,d)}.$

Suppose that $n=d-1$, then for $k=2,4,\ldots,d$

$\displaystyle
b^{n,y}_{k,k}C_{2j-1,\ell(n+1,k)}^{(2,d)}=\sum_{r=j-d+1}^{j+1}r_{r-1,j}C_{2r,\ell(n,k)}^{(2,d)}-\sum_{\begin{subarray}{c}r=\max\\{k-1,1\\}\\\
\mathrm{mod}(r,2)=\mathrm{mod}(k-1,2)\end{subarray}}^{d-1}b^{n-1,y}_{k,r}C_{2j-1,\ell(n-1,r)}^{(2,d)}$
$\displaystyle-\sum_{\begin{subarray}{c}r>0\\\
\mathrm{mod}(r,2)=\mathrm{mod}(k-1,2)\end{subarray}}^{d}a^{n,y}_{r,k}C_{2j-1,\ell(n,r)}^{(2,d)}-\sum_{\begin{subarray}{c}r>0\\\
\mathrm{mod}(r,2)=\mathrm{mod}(k,2)\end{subarray}}^{k-2}b^{n,y}_{r,k}C_{2j-1,\ell(n+1,r)}^{(2,d)}$

and for $k=1,3,\ldots,d$

$\displaystyle
b^{n,y}_{k,k}C_{2j,\ell(n+1,k)}^{(2,d)}=\sum_{r=j-1}^{j-1+d}r_{j-1,r}C_{2r-1,\ell(n,k)}^{(2,d)}-\sum_{\begin{subarray}{c}r=\max\\{k-1,1\\}\\\
\mathrm{mod}(r,2)=\mathrm{mod}(k-1,2)\end{subarray}}^{d-1}b^{n-1,y}_{k,r}C_{2j,\ell(n-1,r)}^{(2,d)}$
$\displaystyle-\sum_{\begin{subarray}{c}r>0\\\
\mathrm{mod}(r,2)=\mathrm{mod}(k-1,2)\end{subarray}}^{d}a^{n,y}_{r,k}C_{2j,\ell(n,r)}^{(2,d)}-\sum_{\begin{subarray}{c}r>0\\\
\mathrm{mod}(r,2)=\mathrm{mod}(k,2)\end{subarray}}^{k-2}b^{n,y}_{r,k}C_{2j,\ell(n+1,r)}^{(2,d)}.$

Suppose that $n\geq d$. If $\mathrm{mod}(n+1-d,2)=1$, then for
$k=1,3,\ldots,d$

$\displaystyle
b^{n,y}_{k,k}C_{2j-1,\ell(n+1,k)}^{(2,d)}=\sum_{r=j-d+1}^{j+1}r_{r-1,j}C_{2r,\ell(n,k)}^{(2,d)}-\sum_{\begin{subarray}{c}r=k\\\
\mathrm{mod}(r,2)=\mathrm{mod}(k,2)\end{subarray}}^{d}b^{n-1,y}_{k,r}C_{2j-1,\ell(n-1,r)}^{(2,d)}$
$\displaystyle-\sum_{\begin{subarray}{c}r>0\\\
\mathrm{mod}(r,2)=\mathrm{mod}(k-1,2)\end{subarray}}^{d}a^{n,y}_{r,k}C_{2j-1,\ell(n,r)}^{(2,d)}-\sum_{\begin{subarray}{c}r>0\\\
\mathrm{mod}(r,2)=\mathrm{mod}(k,2)\end{subarray}}^{k-2}b^{n,y}_{r,k}C_{2j-1,\ell(n+1,r)}^{(2,d)}$

and this equation also holds if $\mathrm{mod}(n+1-d,2)=0$ and for
$k=2,4,\ldots,d$. If $\mathrm{mod}(n+1-d,2)=0$, then for $k=1,3,\ldots,d$

$\displaystyle
b^{n,y}_{k,k}C_{2j,\ell(n+1,k)}^{(2,d)}=\sum_{r=j-1}^{j-1+d}r_{j-1,r}C_{2r-1,\ell(n,k)}^{(2,d)}-\sum_{\begin{subarray}{c}r=k\\\
\mathrm{mod}(r,2)=\mathrm{mod}(k,2)\end{subarray}}^{d}b^{n-1,y}_{k,r}C_{2j,\ell(n-1,r)}^{(2,d)}$
$\displaystyle-\sum_{\begin{subarray}{c}r>0\\\
\mathrm{mod}(r,2)=\mathrm{mod}(k-1,2)\end{subarray}}^{d}a^{n,y}_{r,k}C_{2j,\ell(n,r)}^{(2,d)}-\sum_{\begin{subarray}{c}r>0\\\
\mathrm{mod}(r,2)=\mathrm{mod}(k,2)\end{subarray}}^{k-2}b^{n,y}_{r,k}C_{2j,\ell(n+1,r)}^{(2,d)}$

and this equation also holds if $\mathrm{mod}(n+1-d,2)=1$ and for
$k=2,4,\ldots,d$.

###### Proof.

Note from (5.2), (5.32) and (5.34) that for $m=2$, each of the OPs of degrees
$0$, $1$ and $2$ can be expressed as a linear combination of either (i) the
$P^{(2,d)}_{n,0}$ polynomials or (ii) the $P^{(2,d)}_{n,1}$ polynomials. In
the former case, it is expressed as a polynomial in $x$ and in the latter, as
a $y$ times a polynomial in $x$. Since OPs of degree $n+1$ are generated by
orthogonalizing either $xY_{n,k}^{(2,d)}$ or $yY_{n,k}^{(2,d)}$ for $n\geq 2$
and $y^{2}=\phi(x)$, it follows that every $Y_{n,k}^{(2,d)}$ can be expanded
in terms of either the $P^{(2,d)}_{n,0}$ or the $P^{(2,d)}_{n,1}$ polynomials.
In the former case, $C^{(2,d)}_{2j,\ell(n,k)}=0$, $j\geq 1$ and in the latter,
$C^{(2,d)}_{2j-1,\ell(n,k)}=0$, $j\geq 1$ (see (3.10) and (5.45)). We also
note that in the former case, the monomial expansion of $Y_{n,k}^{(2,d)}$
contains only even powers of $y$ and in the latter case, only odd powers of
$y$.

From the monomial expansions (5.40), for $0\leq n\leq d-1$, we conclude (since
${\alpha}_{k,k}^{(n)}\neq 0$) that if $k$ is odd, then $Y_{n,k}^{(2,d)}$ has
an expansion in the $P^{(2,d)}_{n,0}$ polynomials, and hence
$C^{(2,d)}_{2j,\ell(n,k)}=0$, $j\geq 1$. Using similar reasoning and also
(5.41), we list all the cases:

(A.1) $\begin{split}&1\leq n\leq d-1,\qquad\mathrm{mod}(k,2)=1,\qquad
Y_{n,k}^{(2,d)}(x,y)=\sum C^{(2,d)}_{2j-1,\ell}P_{j,0}^{(2,d)}\\\ &1\leq n\leq
d-1,\qquad\mathrm{mod}(k,2)=0,\qquad Y_{n,k}^{(2,d)}(x,y)=\sum
C^{(2,d)}_{2j,\ell}P_{j,1}^{(2,d)}\\\ &n\geq
d,\qquad\mathrm{mod}(n-d,2)=0,\mathrm{mod}(k,2)=1,\qquad
Y_{n,k}^{(2,d)}(x,y)=\sum C^{(2,d)}_{2j,\ell}P_{j,1}^{(2,d)}\\\ &n\geq
d,\qquad\mathrm{mod}(n-d,2)=0,\mathrm{mod}(k,2)=0,\qquad
Y_{n,k}^{(2,d)}(x,y)=\sum C^{(2,d)}_{2j-1,\ell}P_{j,0}^{(2,d)}\\\ &n\geq
d,\qquad\mathrm{mod}(n-d,2)=1,\mathrm{mod}(k,2)=1,\qquad
Y_{n,k}^{(2,d)}(x,y)=\sum C^{(2,d)}_{2j-1,\ell}P_{j,0}^{(2,d)}\\\ &n\geq
d,\qquad\mathrm{mod}(n-d,2)=1,\mathrm{mod}(k,2)=0,\qquad
Y_{n,k}^{(2,d)}(x,y)=\sum C^{(2,d)}_{2j,\ell}P_{j,1}^{(2,d)}\end{split}$

where $j\geq 1$, $\ell=\ell(n,k)$; in the first two cases, $1\leq k\leq n+1$
and in the final four cases, $1\leq k\leq d$.

Observe that if $Y_{n,k}^{(2,d)}$ has an expansion in the $P_{j,0}^{(2,d)}$ or
$P_{j,1}^{(2,d)}$ polynomials, then so does $xY_{n,k}^{(2,d)}$ but then
$yY_{n,k}^{(2,d)}$ has an expansion in, respectively, the $P_{j,1}^{(2,d)}$ or
$P_{j,0}^{(2,d)}$ polynomials. We conclude from this and the cases listed
above that
$A_{n,x}^{(2,d)}\in{\mathbb{R}}^{\dim{\mathcal{V}}_{n}\times\dim{\mathcal{V}}_{n}}$
is a symmetric matrix whose odd diagonals are zero since
${\langle}xY^{(2,d)}_{n,k_{1}},Y^{(2,d)}_{n,k_{2}}{\rangle}_{{\gamma}_{2,d}}=0$
if $k_{1}$ and $k_{2}$ have opposite parity. This is because, for example, if
$Y^{(2,d)}_{n,k_{1}}$ has an expansion in the $P_{j,1}^{(2,d)}$ polynomials,
then $Y^{(2,d)}_{n,k_{2}}$ has an expansion in the $P_{j,0}^{(2,d)}$
polynomials, which are mutually orthogonal. Similarly,
$A_{n,y}^{(2,d)}\in{\mathbb{R}}^{\dim{\mathcal{V}}_{n}\times\dim{\mathcal{V}}_{n}}$
is a symmetric matrix whose even diagonals are zero since
${\langle}yY^{(2,d)}_{n,k_{1}},Y^{(2,d)}_{n,k_{2}}{\rangle}_{{\gamma}_{2,d}}=0$
if $k_{1}$ and $k_{2}$ have the same parity. It also follows that for $0\leq
n\leq d-1$, the structures of $B^{(2,d)}_{n,x}$ and $B^{(2,d)}_{n,y}$ are the
same as those in, respectively, (5.62) and (5.68) but, additionally, the odd
super-diagonals of $B^{(2,d)}_{n,x}$ are zero and the even super-diagonals of
$B^{(2,d)}_{n,y}$ are zero. For $n\geq d$, the structures of $B^{(2,d)}_{n,x}$
and $B^{(2,d)}_{n,y}$ are the same as those in, respectively, (5.74) and
(5.79) but the even super-diagonals of $B^{(2,d)}_{n,x}$ are zero and the odd
super-diagonals of $B^{(2,d)}_{n,y}$ are zero. To summarise, the structures of
the matrices $A^{(2,d)}_{n,x}$, $A^{(2,d)}_{n,y}$, $B^{(2,d)}_{n,x}$ and
$B^{(2,d)}_{n,y}$ are the same to those in (5.56)–(5.79), except that for
$m=2$ they have an additional checkerboard pattern.

The recurrence relations for the connection coefficients can be derived as in
the case for $m=1$ in Proposition 5.5, but using instead (5.45),

$\displaystyle xP_{j,0}^{(2,d)}$
$\displaystyle={\beta}_{j-1}(w)P_{j-1,0}^{(2,d)}+{\alpha}_{j}(w)P_{j,0}^{(2,d)}+{\beta}_{j}(w)P_{j+1,0}^{(2,d)}$
$\displaystyle xP_{j,1}^{(2,d)}$ $\displaystyle={\beta}_{j-2}(\phi
w)P_{j-1,1}^{(2,d)}+{\alpha}_{j-1}(\phi w)P_{j,1}^{(2,d)}+{\beta}_{j-1}(\phi
w)P_{j+1,1}^{(2,d)}$

(see section 5.1) and (5.20)–(5.21), which we combine with the fact that the
OPs with $m=2$ satisfy the same recurrence relations as the OPs for $m=1$
(e.g., (5.80)–(5.81)), but with the additional checkerboard pattern of
$A_{n,x}^{(2,d)}$, $A_{n,y}^{(2,d)}$, $B_{n,x}^{(2,d)}$ and $B_{n,y}^{(2,d)}$.

∎

##  Appendix B The cases for which $C^{(2,d)}_{j,\ell(n+1,k)}=0$

######  Proposition B.1.

For $m=2$, $C^{(2,d)}_{j,\ell(n+1,k)}=0$ for the following indices: suppose
$0\leq n\leq d-2$, then $C^{(2,d)}_{j,\ell(n+1,k)}=0$ for

(B.1) $j<2n+1\qquad\text{and }\qquad 1\leq k\leq n+2;$

if, additionally, $n$ is even, then $C^{(2,d)}_{j,\ell(n+1,k)}=0$ for

(B.2) $\displaystyle\mathrm{mod}(j,2)=0\quad\text{or}$ $\displaystyle\quad
j>nd-1\quad$ $\displaystyle\text{and}\quad k=1,3,\ldots,n-1$ (B.3)
$\displaystyle\mathrm{mod}(j,2)=1\quad\text{or}$ $\displaystyle\quad
j>(n-2)d+4\quad$ $\displaystyle\text{and }\quad k=2,4,\ldots,n-2$ (B.4)
$\displaystyle\mathrm{mod}(j,2)=1\quad\text{or}$ $\displaystyle\quad
j>(n-2)d+6\quad$ $\displaystyle\text{and }\quad k=n$ (B.5)
$\displaystyle\mathrm{mod}(j,2)=0\quad\text{or}$ $\displaystyle\quad
j>nd+1\quad$ $\displaystyle\text{and }\quad k=n+1$ (B.6)
$\displaystyle\mathrm{mod}(j,2)=1\quad\text{or}$ $\displaystyle\quad
j>nd+2\quad$ $\displaystyle\text{and }\quad k=n+2$

if $n$ is odd (with $0\leq n\leq d-2$), then $C^{(2,d)}_{j,\ell(n+1,k)}=0$ for

$\displaystyle\mathrm{mod}(j,2)=0\qquad\text{or}$ $\displaystyle\qquad
j>1+(n-1)d\qquad$ $\displaystyle\text{and }\qquad k=1,3,\ldots,n-2$
$\displaystyle\mathrm{mod}(j,2)=1\qquad\text{or}$ $\displaystyle\qquad
j>2+(n-1)d\qquad$ $\displaystyle\text{and }\qquad k=2,4,\ldots,n-1$
$\displaystyle\mathrm{mod}(j,2)=0\qquad\text{or}$ $\displaystyle\qquad
j>3+(n-1)d\qquad$ $\displaystyle\text{and }\qquad k=n$
$\displaystyle\mathrm{mod}(j,2)=1\qquad\text{or}$ $\displaystyle\qquad
j>4+(n-1)d\qquad$ $\displaystyle\text{and }\qquad k=n+1$
$\displaystyle\mathrm{mod}(j,2)=0\qquad\text{or}$ $\displaystyle\qquad
j>(n+1)d-1\qquad$ $\displaystyle\text{and }\qquad k=n+2$

Suppose now that $n\geq d-1$. If $\mathrm{mod}(n+1-d,2)=0$, then
$C^{(2,d)}_{j,\ell(n+1,k)}=0$ for

(B.7) $j<2d+d(n+1-d)\qquad\text{and}\qquad 1\leq k\leq d$

and if $\mathrm{mod}(n+1-d,2)=1$, then $C^{(2,d)}_{j,\ell(n+1,k)}=0$ for

(B.8) $j<2d-3+d(n+2-d)\qquad\text{and}\qquad 1\leq k\leq d.$

Suppose $d$ is even and $n\geq d-1$. If $\mathrm{mod}(n+1-d,2)=0$, then $n$ is
odd and $C^{(2,d)}_{j,\ell(n+1,k)}=0$ for

$\displaystyle\mathrm{mod}(j,2)=1\qquad\text{or}$ $\displaystyle\qquad
j>2+(n-1)d\qquad$ $\displaystyle\text{and }\qquad k=1,3,\ldots,d-3$
$\displaystyle\mathrm{mod}(j,2)=0\qquad\text{or}$ $\displaystyle\qquad
j>1+(n-1)d\qquad$ $\displaystyle\text{and }\qquad k=2,4,\ldots,d-4$
$\displaystyle\mathrm{mod}(j,2)=0\qquad\text{or}$ $\displaystyle\qquad
j>3+(n-1)d\qquad$ $\displaystyle\text{and }\qquad k=d-2$
$\displaystyle\mathrm{mod}(j,2)=1\qquad\text{or}$ $\displaystyle\qquad
j>4+(n-1)d\qquad$ $\displaystyle\text{and }\qquad k=d-1$
$\displaystyle\mathrm{mod}(j,2)=0\qquad\text{or}$ $\displaystyle\qquad
j>(n+1)d-1\qquad$ $\displaystyle\text{and }\qquad k=d$

and if $\mathrm{mod}(n+1-d,2)=1$, then $n$ is even and
$C^{(2,d)}_{j,\ell(n+1,k)}=0$ for

$\displaystyle\mathrm{mod}(j,2)=0\qquad\text{or}$ $\displaystyle\qquad
j>nd-1\qquad$ $\displaystyle\text{and }\qquad k=1,3,\ldots,d-3$
$\displaystyle\mathrm{mod}(j,2)=1\qquad\text{or}$ $\displaystyle\qquad
j>(n-2)d+4\qquad$ $\displaystyle\text{and }\qquad k=2,4,\ldots,d-4$
$\displaystyle\mathrm{mod}(j,2)=1\qquad\text{or}$ $\displaystyle\qquad
j>(n-2)d+6\qquad$ $\displaystyle\text{and }\qquad k=d-2$
$\displaystyle\mathrm{mod}(j,2)=0\qquad\text{or}$ $\displaystyle\qquad
j>nd+1\qquad$ $\displaystyle\text{and }\qquad k=d-1$
$\displaystyle\mathrm{mod}(j,2)=1\qquad\text{or}$ $\displaystyle\qquad
j>nd+2\qquad$ $\displaystyle\text{and }\qquad k=d$

Suppose $d$ is odd and $n\geq d-1$. If $\mathrm{mod}(n+1-d,2)=0$, then $n$ is
even and $C^{(2,d)}_{j,\ell(n+1,k)}=0$ for

$\displaystyle\mathrm{mod}(j,2)=1\qquad\text{or}$ $\displaystyle\qquad
j>(n-2)d+4\qquad$ $\displaystyle\text{and }\qquad k=1,3,\ldots,d-4$
$\displaystyle\mathrm{mod}(j,2)=0\qquad\text{or}$ $\displaystyle\qquad
j>nd-1\qquad$ $\displaystyle\text{and }\qquad k=2,4,\ldots,d-3$
$\displaystyle\mathrm{mod}(j,2)=1\qquad\text{or}$ $\displaystyle\qquad
j>(n-2)d+6\qquad$ $\displaystyle\text{and }\qquad k=d-2$
$\displaystyle\mathrm{mod}(j,2)=0\qquad\text{or}$ $\displaystyle\qquad
j>nd+1\qquad$ $\displaystyle\text{and }\qquad k=d-1$
$\displaystyle\mathrm{mod}(j,2)=1\qquad\text{or}$ $\displaystyle\qquad
j>nd+2\qquad$ $\displaystyle\text{and }\qquad k=d$

If $\mathrm{mod}(n+1-d,2)=1$, then $n$ is odd and
$C^{(2,d)}_{j,\ell(n+1,k)}=0$ for

$\displaystyle\mathrm{mod}(j,2)=0\qquad\text{or}$ $\displaystyle
j>1+(n-1)d\qquad$ $\displaystyle\text{and }\qquad k=1,3,\ldots,d-4$
$\displaystyle\mathrm{mod}(j,2)=1\qquad\text{or}$ $\displaystyle
j>2+(n-1)d\qquad$ $\displaystyle\text{and }\qquad k=2,4,\ldots,d-3$
$\displaystyle\mathrm{mod}(j,2)=0\qquad\text{or}$ $\displaystyle
j>3+(n-1)d\qquad$ $\displaystyle\text{and }\qquad k=d-2$
$\displaystyle\mathrm{mod}(j,2)=1\qquad\text{or}$ $\displaystyle
j>4+(n-1)d\qquad$ $\displaystyle\text{and }\qquad k=d-1$
$\displaystyle\mathrm{mod}(j,2)=0\qquad\text{or}$ $\displaystyle
j>(n+1)d-1\qquad$ $\displaystyle\text{and }\qquad k=d$

###### Proof.

It follows from the cases listed in (A.1) whether, for a given column
$\mathbf{C}^{(2,d)}_{\ell(n,k)}$, the even-indexed rows or the odd-indexed
rows are zero. For example, for $0\leq n\leq d-1$, if $k$ is odd, then
$C^{(2,d)}_{2j,\ell(n,k)}=0$ for $j\geq 1$.

To determine the upper and lower bounds on the summation indices in (A.1), we
use the same arguments as in Proposition 5.6, which we present here with more
brevity.

For $0\leq n\leq d-2$,

(B.9) $\deg
x^{n+2-k}y^{k-1}=\begin{cases}n+2-k+\displaystyle{d\left(\frac{k-1}{2}\right)},\qquad
k=1,3,\ldots,\begin{cases}n+1&n\text{ even}\\\ n+2&n\text{ odd}\end{cases}\\\
n+3-k+\displaystyle{d\left(\frac{k-2}{2}\right)},\qquad
k=2,4,\ldots,\begin{cases}n+2&n\text{ even}\\\ n+1&n\text{
odd}\end{cases}.\end{cases}$

For $n\geq d-1$ with $\mathrm{mod}(n+1-d,2)=0$,

(B.10) $\deg
x^{d-k}y^{n+1+k-d}=\begin{cases}d-k+1+\displaystyle{d\left(\frac{n+k-d}{2}\right)},\qquad
k=1,3,\ldots,\begin{cases}d&n\text{ even}\\\ d-1&n\text{ odd}\end{cases}\\\
d-k+\displaystyle{d\left(\frac{n+1+k-d}{2}\right)},\qquad
k=2,4,\ldots,\begin{cases}d-1&n\text{ even}\\\ d&n\text{
odd}\end{cases}\end{cases}$

and if $n\geq d-1$ and $\mathrm{mod}(n+1-d,2)=1$, then

(B.11) $\deg
x^{d-k}y^{n+1+k-d}=\begin{cases}d-k+\displaystyle{d\left(\frac{n+1+k-d}{2}\right)},\qquad
k=1,3,\ldots,\begin{cases}d-1&n\text{ even}\\\ d&n\text{ odd}\end{cases}\\\
d-k+1+\displaystyle{d\left(\frac{n+k-d}{2}\right)},\qquad
k=2,4,\ldots,\begin{cases}d&n\text{ even}\\\ d-1&n\text{
odd}\end{cases}.\end{cases}$

By the same reasoning used in Proposition 5.6, we conclude from (B.9) that if
$n$ is even, with $0\leq n\leq d-2$, then

(B.12) $\deg Y^{(2,d)}_{n+1,k}=\begin{cases}dn/2&k=1,3,\ldots,n-1\\\
2+d(n-2)/2&k=2,4,\ldots,n-2\\\ 3+d(n-2)/2&k=n\\\ 1+dn/2&k=n+1,n+2\end{cases}.$

For $0\leq n\leq d-2$, it follows from (B.9) that

(B.13)
$\displaystyle\deg_{{\gamma}_{2,d}}(P^{(2,d)}_{j,0})<n+1\quad\text{and}\quad\deg_{{\gamma}_{2,d}}(P^{(2,d)}_{j,1})<n+1\quad\text{for}\quad
j<n+1.$

Combining the first and second cases in (A.1) with (B.12) and (B.13), it
follows that if $n$ is even, with $0\leq n\leq d-2$, then

(B.14)
$Y^{(2,d)}_{n+1,k}=\sum_{j=n+1}^{dn/2}C^{(2,d)}_{2j-1,\ell(n+1,k)}P_{j,0}^{(2,d)},\qquad
k=1,3,\ldots,n-1\\\ $ (B.15)
$Y^{(2,d)}_{n+1,k}=\sum_{j=n+1}^{2+d(n-2)/2}C^{(2,d)}_{2j,\ell(n+1,k)}P_{j,1}^{(2,d)},\qquad
k=2,4,\ldots,n-2\\\ $ (B.16)
$Y^{(2,d)}_{n+1,k}=\sum_{j=n+1}^{3+d(n-2)/2}C^{(2,d)}_{2j,\ell(n+1,k)}P_{j,1}^{(2,d)},\qquad
k=n\\\ $ (B.17)
$Y^{(2,d)}_{n+1,k}=\sum_{j=n+1}^{1+dn/2}C^{(2,d)}_{2j-1,\ell(n+1,k)}P_{j,0}^{(2,d)},\qquad
k=n+1\\\ $ (B.18)
$Y^{(2,d)}_{n+1,k}=\sum_{j=n+1}^{1+dn/2}C^{(2,d)}_{2j,\ell(n+1,k)}P_{j,1}^{(2,d)},\qquad
k=n+2.$

Note that in (B.14), the row indices of $C^{(2,d)}_{j,\ell(n+1,k)}$ range from
$2n+1$ to $dn-1$, and hence the statements in (B.1) and (B.2) that
$C^{(2,d)}_{j,\ell(n+1,k)}=0$ for $j<2n+1$ and $j>dn-1$. Similarly,
(B.3)–(B.6) follow from, respectively, (B.15)–(B.18).

The proof that $C^{(2,d)}_{j,\ell(n+1,k)}=0$ for the remaining cases is
similar. ∎