# Characterization of classical orthogonal polynomials in two variables

M. KENFACK NANGHOa, K. JORDAANb and B. JIEJIP NKWAMOUOc

a,cDepartment of Mathematics and Computer Science,
University of Dschang, PO Box 67, Dschang, Cameroon
b Department of Decision Sciences,
University of South Africa, PO Box 392, Pretoria, 0003, South Africa

###### Abstract

For a family of polynomials in two variables, orthogonal with respect to a
weight function, we prove under some conditions, equivalence between: the
Matrix Pearson equation of the weight, the second order linear partial
differential equation, the orthogonality of the gradients, the Matrix
Rodrigues formula involving tensor product of Matrices, and the so-called
first structure relation. We then propose a definition of classical orthogonal
polynomials in two variables.

  * Keywords:Orthogonal polynomials in two variables, Characterization, Partial Differential Equations.

  * Mathematics Subject Classification (2020): Primary 42C05; Secondary 33C50

##  1 Introduction and Notations

###  1.1 Introduction

Let $\left\\{p_{n}\right\\}_{n\in\mathbb{N}}$ be a family of polynomials,
orthogonal with respect to a weight function, $\rho$, on an interval (a,b) of
$\mathbb{R}$ (set of real numbers). $(p_{n})_{n}$ is classical if and only if
its derivatives are also orthogonal. This property leads to the existence of a
sequence of real numbers $\mu_{n,m-1}$, $n\in\mathbb{N}$, $0\leq m\leq n$ and
a polynomial $\phi$ of degree at most 2 such that [16, Eq. (1.3.5)]

$\left[\rho(x)\phi^{m}p_{n}^{(m)}\right]^{(1)}(x)=-\mu_{n,m-1}\rho(x)\phi^{m-1}(x)p_{n}^{(m-1)}(x),$

where $p_{n}^{(m)}=\frac{d^{m}}{dx^{m}}p_{n}$. Iterating and taking $m=n$, one
has

$p_{n}(x)=\frac{(-1)^{n}p_{n}^{(n)}(x)}{\rho(x)\prod_{j=0}^{n-1}\mu_{n,j}}\left[\rho\,\phi^{n}\right]^{(n)}(x)$
(1)

which is the Rodrigues formula. Characterization of classical orthogonal
polynomials in one variable is obtained by proving equivalence between (1),
the Pearson equation of the weight, orthogonality of derivatives, a Sturm-
Liouville type equation, the so-called first structure relation and a non-
linear relation (see [1]). These properties are very important when using
classical orthogonal polynomials in probability, partial differential
equations and mathematical physics (see [17, 7] and references therein). Since
the work of Krall and Sheffer [11] introducing classical polynomials in two
variables as polynomial solutions to the linear partial differential equation

$\displaystyle(\alpha x^{2}+d_{1}x+e_{1}y+f_{1})\partial_{x}^{2}Y+(2\alpha
xy+d_{2}x+e_{2}y+f_{2})\partial_{x}\partial_{y}Y$ $\displaystyle+(\alpha
y^{2}+d_{3}x+e_{3}y+f_{3})\partial_{y}^{2}Y+(\delta\,x+h_{1})\partial_{x}Y+(\delta\,y+h_{2})\partial_{x}Y+\lambda_{n}Y=0,$
(2)

$\lambda_{n}\neq 0$ for $n\neq 0$ where $n$ is the degree of the polynomial
solution, theory and application of orthogonal polynomials in two variables
has been attracting researchers from various domains of mathematics [8, 6, 3].
In 1975, Koornwinder [8] studied examples of two-variable analogues of Jacobi
polynomials, and he introduced seven classes of orthogonal polynomials which
he considered to be bivariate analogues of Jacobi polynomials. Some of these
polynomials are classical according to the Krall and Sheffer definition and
others are not. Therefore the Krall and Sheffer classification seems to be
incomplete. Lyskova [13] studied conditions on the polynomial coefficients in
(1.1) in such a way that the partial derivatives of orthogonal polynomial
solutions satisfy a partial differential equation of the same type. He
obtained that $e_{1}=d_{3}=0$, i.e. the polynomial coefficients of the partial
derivatives in $x$ (respectively in $y$) depend only on $x$ (respectively on
$y$). An analogue of Rodrigues’ formula for Krall and Sheffer classical
orthogonal polynomials in two variables has been obtained by Suetin [18]. In
fact for $n$ a positive integer he defines

$p_{n-i,i}=\frac{1}{\rho}\partial_{x}^{n-i}\partial_{y}^{i}\left(p^{n-i}q^{i}\rho\right),$

where $\rho(x,\,y)$ is a weight function over a simply connected domain and
$p(x,\,y)$ $q(x,\,y)$ are polynomials related with the polynomial coefficients
of (1.1). However, in several cases, Rodrigues formula (1) provides
polynomials of total degree greater than $n$ (see [12] ) and therefore it can
not be used to construct a basis of orthogonal polynomials in two variables.
In [5], by using matrix notation of polynomials in several variables
introduced in [9, 10] and modified in [19], the authors extended the
definition by Krall and Sheffer of classical orthogonal polynomials in two
variables as matrix polynomial families satisfying a second order linear
partial differential equation analogue to (1.1), whose coefficients are some
polynomials of degree at most two without any restriction on their shape.
Using this notation, Álvarez de Morales et. al. [2] obtained the matrix
Rodrigues formula

$\mathbb{P}^{t}_{n}=\frac{1}{\rho}div^{\\{n\\}}\left(\mathbb{\varPhi}^{\\{n\\}}\rho\right),$

where $\rho$ is the nontrivial function, $C^{2}$ in some open set, solution to
the Matrix Pearson-type equation

$div(\varPhi\,\rho)=\varPsi^{t}\,\rho,\varPhi=\left(\begin{matrix}\phi_{1,1}&\phi_{1,2}\\\
\phi_{2,1}&\phi_{2,2}\end{matrix}\right),\varPsi=\left(\begin{matrix}\psi_{1}\\\
\psi_{2}\end{matrix}\right)$ (3)

$\phi_{i,j},i,j=1,2$ and $\psi_{l},\,l=1,2$ are polynomials in two variables
of degree at most 2 and 1 respectively; $\varPhi^{\\{n\\}}$ denotes the second
kind Kronecker power [4] of the matrix $\varPhi$; $div^{\\{n\\}}$ is a
$n^{th}$ order differential operator. For $n=1$, $div^{\\{1\\}}$ is the usual
divergence operator. In 2018, Marcellán et. al. [14] proved that the moment
functional is classical if and only if its moments satisfy two three-term
relations. In this work using the vector notation of polynomials in two
variables, we establish the following extension of (1) to classical orthogonal
polynomials in two variables

$\displaystyle\mathbb{P}^{t}_{n}=\frac{(-1)^{n}}{\rho}div^{(n)}\left[\rho\,\varPhi^{\otimes\,n}\right]\nabla^{(n)}\mathbb{P}_{n}^{t}\prod_{j=0}^{n-1}\varLambda_{n,j}^{-1}$
(4)

where $\varLambda_{n,j}\in\mathcal{M}_{n+1}(\mathbb{R})$, space of
$(n+1,\,n+1)$ matrices with coefficients in $\mathbb{R}$, $div^{(n)}$ is the
usual $n^{th}$ divergence and $\nabla^{(n)}$ is $n^{th}$ gradient; $\varPhi$
is the matrix in (3) and $\rho$ is solution to the same equation.
$\varPhi^{\otimes\,n}$ is the $n^{th}$ Kronecker product (i.e., the tensor
product). Under the boundary condition

$\int_{\Omega}div\left(\rho\,\varPhi\nabla P\right)dxdy=0,$

where $P$ is a polynomial in two variables, and the differential conditions on
the symmetric matrix polynomials $\Phi$

$\displaystyle\phi_{1,1}\partial_{x}\varPhi+\phi_{2,1}\partial_{y}\varPhi=\varPhi\nabla\left(\phi_{1,1},\phi_{2,1}\right),$
$\displaystyle\phi_{1,2}\partial_{x}\varPhi+\phi_{2,2}\partial_{y}\varPhi=\varPhi\nabla\left(\phi_{1,2},\phi_{2,2}\right),$

without any restriction on the shape of its coefficients, we prove equivalence
between: the Matrix Pearson equation; the second order linear partial
differential equation; the orthogonality of the gradients; the Matrix
Rodrigues formula (4), and the first structure relation.

###  1.2 Notations

Let us recall some fundamental notations and results.
$\mathcal{P}=\mathbb{R}[x,y]$ is the space of polynomials in two variables
with real coefficients. $\mathcal{P}_{n}$ is the subspace of polynomials of
degree (also called total degree ) less or equal to $n$. $\Pi_{n}$ is the
subspace of polynomials of degree $n$. Let $\mathbb{N}$ be the set of non-
negative integers. For $\alpha=(\alpha_{1},\;\alpha_{2})\in\mathbb{N}^{2}$,
$|\alpha|=\alpha_{1}+\alpha_{2}$. Let $p\in\mathcal{P}$ be a polynomial of
degree $n$,

$\displaystyle p(x,y)$ $\displaystyle=$ $\displaystyle\sum_{|\alpha|\leq
n}c_{\alpha}x^{\alpha_{1}}y^{\alpha_{2}},\;\;c_{\alpha}\in\mathbb{R},$ (6)
$\displaystyle=$ $\displaystyle\sum_{k=0}^{n}C_{\alpha}X_{k},$

where $X_{k}$ is the column vector $X_{0}=1$,
$X_{k}=(x^{k},x^{k-1}y,...,xy^{k-1},\,y^{k})^{t}$, $k\geq 1$, of size
$(k+1,1)$ and $C_{k}\in\mathcal{M}_{1,\,k+1}(\mathbb{R})$.
$\mathcal{M}_{h,\,k}(\mathbb{R})$ and $\mathcal{M}_{h,\,k}(\mathcal{P})$
denote the linear spaces of $h$ line(s) and $k$ column(s) real and polynomial
matrices respectively. When $h=k$, the second index is omitted. The degree of
matrix polynomial $\varPhi\in\mathcal{M}_{h,\,k}(\mathcal{P})$ is

${\rm degree}\,{\varPhi}=max\left\\{{\rm degree}\,\phi_{ij},1\leq i\leq
h,\;1\leq j\leq k\right\\},$

where $\phi_{ij}$ denotes (i,j)-entry of $\varPhi$. A polynomial system is a
sequence of vectors $\left\\{\mathbb{P}_{n}\right\\}_{n\geq 0}$ of increasing
size such that [10]

$\mathbb{P}_{n}=\left(P_{n,0},\;P_{n-1,1},...,P_{1,n-1},\;P_{0,n}\right)^{t}\in\mathcal{M}_{n+1,1}\left(\mathcal{P}_{n}\right),$

where $P_{i,j}$, $0\leq i\leq n$, $0\leq j\leq n$, are polynomials of degree
$n$ independent modulo $\mathcal{P}_{n-1}$. Expanding $P_{i,j}$ in the basis
$X_{k}$ and using (6), $\mathbb{P}_{n}(x,y)$ can be written as follows

$\mathbb{P}_{n}(x,y)=\sum_{k=0}^{n}G_{n,k}X_{k},$

where $G_{n,k}\in\mathcal{M}_{n+1,\,k+1}\left(\mathbb{R}\right)$ and $G_{n,n}$
is the leading coefficient. $\mathbb{P}_{n}$ is monic if $G_{n,n}$ is equal to
the matrix identity, $I_{n+1}$, of $\mathcal{M}_{n+1}\left(\mathbb{R}\right)$.
Let $\Omega$ be a domain of $\mathbb{R}^{2}$ and $\rho$ a weight function on
$\Omega$ (i.e. a non-negative and integrable function on $\Omega$). The system
of polynomials $\left\\{\mathbb{P}_{n}\right\\}_{n}$ is orthogonal with
respect to $\rho(x,y)$ if for $n\geq 0$

$\begin{cases}\int_{\Omega}X_{m}\mathbb{P}^{t}_{n}(x,y)\rho(x,y)dxdy&=0,\;\;m<n,\\\
\int_{\Omega}X_{n}\mathbb{P}^{t}_{n}(x,y)\rho(x,y)dxdy&=S_{n},\end{cases}$

where $S_{n}$ is an invertible matrix of size $(n+1,\,n+1)$. For $n\geq 0$,
the multiplication of $X_{n}$ by $x$ and $y$ are respectively given by

$xX_{n}=L_{n,1}X_{n+1}\;\;{\rm and}\;\;yX_{n}=L_{n,2}X_{n+1},$ (7)

where $L_{n,1}$ and $L_{n,2}$ are the $(n+1,\,n+2)$ matrices

$L_{n,1}=\left(\begin{matrix}1&0&\cdots&0&0\\\ 0&1&&\vdots&0\\\
\vdots&&\ddots&0&\vdots\\\ 0&\cdots&0&1&0\end{matrix}\right)\hskip
28.45274pt{\rm and}\hskip
28.45274ptL_{n,2}=\left(\begin{matrix}0&1&0&\cdots&0\\\ 0&0&1&&\vdots\\\
\vdots&\vdots&&\ddots&0\\\ 0&0&\cdots&0&1\end{matrix}\right).$ (8)

On the other hand, the partial derivative of $X_{n}$ with respect to $x$ and
$y$ are defined by [14]

$\partial_{x}X_{n}=N_{n,1}^{t}X_{n-1},\hskip
56.9055pt\partial_{y}X_{n}=N_{n,2}^{t}X_{n-1},$ (9)

where $N_{n,k}$, $k=1,2$ are $(n,\,n+1)$ matrices defined as follow
$N_{0,1}=N_{0,2}=0$, for $n\geq 1$

$N_{n,1}=\left(\begin{matrix}n&0&\cdots&0&0\\\ 0&n-1&&\vdots&0\\\
\vdots&&\ddots&0&\vdots\\\ 0&\cdots&0&1&0\end{matrix}\right)\hskip
28.45274pt{\rm and}\hskip
28.45274ptN_{n,2}=\left(\begin{matrix}0&1&0&\cdots&0\\\ 0&0&2&&\vdots\\\
\vdots&\vdots&&\ddots&0\\\ 0&0&\cdots&0&n\end{matrix}\right).$ (10)

$A=\left(a_{ij}\right)_{i,j=1}^{m}\in\mathcal{M}_{m}\left(P\right)$ and
$B=\left(b_{ij}\right)_{i,j=1}^{n}\in\mathcal{M}_{n}\left(P\right)$. The
Kronecker product of $A$ and $B$ denoted $A\otimes B$ is the block matrix [4,
p.228]

$A\otimes B=\left(\begin{matrix}a_{11}B&\cdots&a_{1n}B\\\ \vdots&&\vdots\\\
a_{n1}B&\cdots&{a_{nn}B}\end{matrix}\right).$ (11)

For a positive integer $k\geq 1$, the $k^{th}$ Kronecker product of a matrix
$A\in\mathcal{M}_{m}\left(\mathcal{P}\right)$ is the matrix of size $m^{k}$
written as

$A^{\otimes\,k}=A\otimes A^{\otimes\,k-1},\;\;A^{\otimes\,0}=1.$

Moreover, if $A$ and $C$ are matrices of the same size, $B$ and $D$ are
matrices of the same size, then [4, p.228]

$\left(AC\right)\otimes\left(BD\right)=(A\otimes\,B)(C\otimes\,D).$ (12)

Let $A,\,B\in\mathcal{M}_{m,n}(\mathcal{P})$. The gradient of $A$ and the
divergence of $\left(\begin{matrix}A\\\ B\end{matrix}\right)$ are given
respectively by $\nabla A=\left(\begin{matrix}\partial_{x}A\\\
\partial_{y}A\end{matrix}\right)$, and $div\left(\begin{matrix}A\\\
B\end{matrix}\right)=\partial_{x}A+\partial_{y}B.$ The scalar product of
$\left(\begin{matrix}A\\\ B\end{matrix}\right)$ by a vector
$\overrightarrow{n}\left(\begin{matrix}n_{x}\\\ n_{y}\end{matrix}\right)$ is

$\left(\begin{matrix}A\\\
B\end{matrix}\right)\cdot\overrightarrow{n}=An_{x}+Bn_{y}.$

The norm of a matrix $A=(a_{ij})_{i=1,j=1}^{m,n}$ is given by
$\|A\|_{max}=max_{i,j}|a_{ij}|$. Let be $\Omega$ a subset of $\mathbb{R}^{2}$.
$1_{\Omega}$ is the indicator function of $\Omega$. If the boundary of
$\Omega$ exists, we denote it by $\partial\Omega$.
$B(O,\varepsilon),\,\varepsilon>0$, is the open ball of $\mathbb{R}^{2}$
centered at $O$ with radius $\varepsilon$.

##  2 Fundamental results

This subsection contains results that will be used for the proof of the main
theorem of this work.

######  Proposition 1.

Let $m\geq 0$ be a fixed integer and $\left\\{I_{2^{m}}\otimes
X^{t}_{n},\;n\in\mathbb{N}\right\\}$ a family of matrix vector polynomials.

  1. 1.

For all $n\geq\,0$, the following algebraic properties are satisfied

$\displaystyle x\left(I_{2^{m}}\otimes X^{t}_{n}\right)=\left(I_{2^{m}}\otimes
X^{t}_{n+1}\right)\left(I_{2^{m}}\otimes
L_{n,1}^{t}\right),\;\;y\left(I_{2^{m}}\otimes
X^{t}_{n}\right)=\left(I_{2^{m}}\otimes
X^{t}_{n+1}\right)\left(I_{2^{m}}\otimes L_{n,2}^{t}\right),$ (13a)
$\displaystyle x^{2}\left(I_{2^{m}}\otimes
X^{t}_{n}\right)=\left(I_{2^{m}}\otimes
X^{t}_{n+2}\right)\left(I_{2^{m}}\otimes L_{n+1,1}^{t}L_{n,1}^{t}\right),$
(13b) $\displaystyle xy\left(I_{2^{m}}\otimes
X^{t}_{n}\right)=\left(I_{2^{m}}\otimes
X^{t}_{n+2}\right)\left(I_{2^{m}}\otimes L_{n+1,1}^{t}L_{n,2}^{t}\right),$
(13c) $\displaystyle y^{2}\left(I_{2^{m}}\otimes
X^{t}_{n}\right)=\left(I_{2^{m}}\otimes
X^{t}_{n+2}\right)\left(I_{2^{m}}\otimes L_{n+1,2}^{t}L_{n,2}^{t}\right),$
(13d) $\displaystyle\left(I_{2^{m}}\otimes X^{t}\right)A\left(I_{2^{m}}\otimes
X^{t}_{n}\right)=\left(I_{2^{m}}\otimes
X^{t}_{n+1}\right)\left(I_{2^{m}}\otimes
L_{n}^{t}\right)\left(A\otimes\,I_{n+1}\right),\;A\in\mathcal{M}_{2^{m+1},2^{m}}\left(\mathbb{R}\right),$
(13e)

  2. 2.

For all $n\geq\,1$, the following partial differentiation properties are
satisfied

$\displaystyle\partial_{x}\left(I_{2^{m}}\otimes
X^{t}_{n}\right)=\left(I_{2^{m}}\otimes
X^{t}_{n-1}\right)\left(I_{2^{m}}\otimes
N_{n,1}\right),\;\partial_{y}\left(I_{2^{m}}\otimes
X^{t}_{n}\right)=\left(I_{2^{m}}\otimes
X^{t}_{n-1}\right)\left(I_{2^{m}}\otimes N_{n,2}\right),$ (14a)
$\displaystyle\partial_{x}^{2}\left(I_{2^{m}}\otimes
X^{t}_{n}\right)=\left(I_{2^{m}}\otimes
X^{t}_{n-2}\right)\left(I_{2^{m}}\otimes N_{n-1,1}N_{n,1}\right),$ (14b)
$\displaystyle\partial_{x}\partial_{y}\left(I_{2^{m}}\otimes
X^{t}_{n}\right)=\left(I_{2^{m}}\otimes
X^{t}_{n-2}\right)\left(I_{2^{m}}\otimes N_{n-1,1}N_{n,2}\right),$ (14c)
$\displaystyle\partial_{y}^{2}\left(I_{2^{m}}\otimes
X^{t}_{n}\right)=\left(I_{2^{m}}\otimes
X^{t}_{n-2}\right)\left(I_{2^{m}}\otimes N_{n-1,2}N_{n,2}\right),$ (14d)

with

$L_{n}=\left(\begin{matrix}L_{n,1}\\\ L_{n,2}\end{matrix}\right),\;{\rm
and}\;\;N_{n}=\left(\begin{matrix}N_{n,1}\\\ N_{n,2}\end{matrix}\right).$

###### Proof.

Since

$\left(I_{2^{m}}\otimes\,X^{t}_{n}\right)=\left(\begin{matrix}X_{n}^{t}&0&\cdots&0\\\
0&X_{n}^{t}&&\vdots\\\ \vdots&&\ddots&0\\\
0&\cdots&0&X_{n}^{t}\end{matrix}\right)$ (15)

using (7), we have

$\displaystyle x\left(I_{2^{m}}\otimes\,X^{t}_{n}\right)$ $\displaystyle=$
$\displaystyle\left(\begin{matrix}X_{n+1}^{t}L_{n,1}&0&\cdots&0\\\
0&X_{n+1}^{t}L_{n,1}&&\vdots\\\ \vdots&&\ddots&0\\\
0&\cdots&0&X_{n+1}^{t}L_{n,1}\end{matrix}\right)$ $\displaystyle=$
$\displaystyle\left(\begin{matrix}X_{n+1}^{t}&0&\cdots&0\\\
0&X_{n+1}^{t}&&\vdots\\\ \vdots&&\ddots&0\\\
0&\cdots&0&X_{n+1}^{t}\end{matrix}\right)\left(\begin{matrix}L_{n,1}&0&\cdots&0\\\
0&N_{n,1}&&\vdots\\\ \vdots&&\ddots&0\\\
0&\cdots&0&L_{n,1}\end{matrix}\right)$ $\displaystyle=$
$\displaystyle\left(I_{2^{m}}\otimes\,X_{n+1}^{t}\right)\left(I_{2^{m}}\otimes\,L_{n,1}\right).$

In a similar way, we obtain $y\left(I_{2^{m}}\otimes\,X^{t}_{n}\right)$ as
well as (13b)-(13d). Let us prove (13e). Let
$A\in\mathcal{M}_{2^{m+1},2^{m}}\left(\mathbb{R}\right)$.
$A=\left(a_{ij}\right)_{i,j=1}^{2^{m}}$, with
$a_{ij}=\left(a_{ij}^{(1)},\;a_{ij}^{(2)}\right)^{t}$

$\displaystyle\left(I_{2^{m}}\otimes X^{t}\right)A\left(I_{2^{m}}\otimes
X^{t}_{n}\right)=\left(X^{t}a_{ij}X_{n}^{t}\right)_{i,j=1}^{m}=\left(a_{ij}^{(1)}xX_{n}^{t}+a_{ij}^{(2)}yX_{n}^{t}\right)_{i,j=1}^{2^{m}}$

Taking into account (7), we obtain

$\displaystyle\left(I_{2^{m}}\otimes X^{t}\right)A\left(I_{2^{m}}\otimes
X^{t}_{n}\right)$ $\displaystyle=$
$\displaystyle\left(X_{n+1}^{t}\left(a_{ij}^{(1)}L_{n,1}^{t}+a_{ij}^{(2)}L_{n,2}^{t}\right)\right)_{i,j=1}^{2^{m}}$
$\displaystyle=$
$\displaystyle\left(I_{2^{m}}\otimes\,X_{n+1}^{t}\right)\left(a_{ij}^{(1)}L_{n,1}^{t}+a_{ij}^{(2)}L_{n,2}^{t}\right)_{i,j=1}^{2^{m}}$

Observing from the transposed version of (7) that
$a_{ij}^{(1)}L_{n,1}^{t}+a_{ij}^{(2)}L_{n,2}^{t}=L_{n}^{t}\left(\begin{matrix}a_{ij}^{(1)}I_{n+1}\\\
a_{ij}^{(2)}I_{n+1}\end{matrix}\right)$, we obtain

$\displaystyle\left(I_{2^{m}}\otimes X^{t}\right)A\left(I_{2^{m}}\otimes
X^{t}_{n}\right)$ $\displaystyle=$
$\displaystyle\left(I_{2^{m}}\otimes\,X_{n+1}^{t}\right)\left(I_{2^{m}}\otimes\,L_{n}^{t}\right)\left(\begin{matrix}a_{ij}^{(1)}I_{n+1}\\\
a_{ij}^{(2)}I_{n+1}\end{matrix}\right)_{i,j=1}^{2^{m}}$ $\displaystyle=$
$\displaystyle\left(I_{2^{m}}\otimes\,X_{n+1}^{t}\right)\left(I_{2^{m}}\otimes\,L_{n}^{t}\right)\left(A\otimes\,I_{n+1}\right).$

To prove (14a), apply $\partial_{x}$ from the left of (15) and use the
transposed version of (9) to obtain the result. (14b)-(14d) are obtained in a
similar way. ∎

######  Proposition 2.

  1. 1.

Let A and B be two matrices, the following product rule is satisfied

$\partial_{x}(A\otimes B)=\partial_{x}A\otimes B+A\otimes\partial_{x}B.$ (16)

  2. 2.

For all non-negative integers $n$ and $m$, $0\leq m\leq n$, the polynomial
$Q_{n,m}=\nabla^{(m)}\mathbb{P}^{t}_{n+m}(x,\,y)$ of total degree $n$ can be
expanded in the system as $\left\\{\left(I_{2^{m}}\otimes
X^{t}_{k}\right)\right\\}_{k}$

$Q_{n,m}(x,\,y)=\sum_{k=0}^{n}\left(I_{2^{m}}\otimes
X^{t}_{k}\right)G_{n,m,k},$

where the leading coefficient of $\mathbb{Q}_{n,m}$ is the
$(2^{m}(n+1),n+m+1)$-matrix given by the recurrence relation

$G_{n,\,m,\,n}=\left(\begin{matrix}I_{2^{m-1}}\otimes\,N_{n+1,1}&0\\\
0&I_{2^{m-1}}\otimes\,N_{n+1,2}\end{matrix}\right)\left(\begin{matrix}G_{n+1,\,m-1,\,n+1}\\\
G_{n+1,\,m-1,\,n+1}\end{matrix}\right),m\geq\,1,\;\;G_{n,0,n}=I_{n+1}.$ (17)

###### Proof.

Write $A\otimes\,B$ in the form (11) and use the usual product rule to obtain
(16). Let us prove the second item by induction on $m$, $m\geq\,0$. The item
is obviously satisfied for $m=0$. Let us assume it is satisfied up to $m-1$.
Therefore

$\mathbb{Q}_{n,m}(x,y)=\nabla\left(\nabla^{(m-1)}\mathbb{P}_{n+1+m-1}^{t}\right)=\nabla\left(\sum_{k=0}^{n+1}\left(I_{2^{m-1}}\otimes\,X_{k}^{t}\right)G_{n+1,m-1,k}\right).$

Next, we use (14a) and rewrite the obtained sums from 0 to get

$\displaystyle\mathbb{Q}_{n,m}(x,y)$ $\displaystyle=$
$\displaystyle\left(\begin{matrix}\sum_{k=0}^{n}\left(I_{2^{m-1}}\otimes\,X_{k}^{t}\right)\left(I_{2^{m-1}}\otimes\,N_{k+1,1}\right)G_{n+1,m-1,k+1}\\\
\\\
\sum_{k=0}^{n}\left(I_{2^{m-1}}\otimes\,X_{k}^{t}\right)\left(I_{2^{m-1}}\otimes\,N_{k+1,2}\right)G_{n+1,m-1,k+1}\end{matrix}\right)$
$\displaystyle=$
$\displaystyle\sum_{k=0}^{n}\left(\begin{matrix}I_{2^{m-1}}\otimes\,X_{k}^{t}&0\\\
0&I_{2^{m-1}}\otimes\,X_{k}^{t}\end{matrix}\right)\left(\begin{matrix}I_{2^{m-1}}\otimes\,N_{k+1,1}&0\\\
0&I_{2^{m-1}}\otimes\,N_{k+1,2}\end{matrix}\right)\left(\begin{matrix}G_{n+1,m-1,k+1}\\\
G_{n+1,m-1,k+1}\end{matrix}\right)$ $\displaystyle=$
$\displaystyle\sum_{k=0}^{n}\left(I_{2^{m}}\otimes\,X_{k}^{t}\right)G_{n,m,k}.$

∎

######  Proposition 3.

Let $\Omega$ be a simply connected open subset of $\mathbb{R}^{2}$ and $\rho$
a weight function on $\Omega$. Let be
$A\in\mathcal{M}_{2}\left(\mathcal{P}\right)$, $\Omega_{j}=\Omega\cap
B(O,j),\,j=1,2...,$ and $\overrightarrow{n_{j}}$ an outer vector normal of
$\partial\Omega_{j}$. If the Neumann type boundary condition

$\lim\limits_{j\rightarrow\infty}1_{\partial\Omega_{j}}\left(\rho
A\,\nabla\,u\right)\cdot\overrightarrow{n_{j}}=0,\,\,\forall u\in\mathcal{P}$
(18)

is satisfied, then

$\int_{\Omega}div\left(\rho\left(I_{2}\otimes\,M\right)A^{\otimes\,m}\nabla\,N\right)dxdy=0,$
(19)

for all $M\in\mathcal{M}_{p,2^{m-1}}\left(\mathcal{P}\right)$ and
$N\in\mathcal{M}_{2^{m-1},r}\left(\mathcal{P}\right)$, where $m$, $p$ and $r$
are positive integers.

###### Proof.

First prove that
$\lim\limits_{j\rightarrow\infty}1_{\partial\Omega_{j}}\left(\rho\left(I_{2}\otimes\,M\right)A^{\otimes\,m}\nabla
N\right)\cdot\overrightarrow{n_{j}}=0$ then use the divergence theorem as well
as the Lebesgue’s dominated convergence theorem to obtain (19). For this
purpose, observe that

$\rho\left(\left(I_{2}\otimes\,M\right)A^{\otimes\,m}\nabla
N\right)\cdot\overrightarrow{n}=\rho\left(a_{11}n_{x}+a_{21}n_{y}\right)MA^{\otimes\,m-1}\partial_{x}\,N+\rho\left(a_{12}n_{x}+a_{22}n_{y}\right)MA^{\otimes\,m-1}\partial_{y}\,N$

and take $MA^{\otimes\,m-1}\partial_{x}\,N=(p_{ik})_{i=1,k=1}^{p,r}$,
$MA^{\otimes\,m-1}\partial_{y}\,N=(q_{ls})_{l=1,s=1}^{p,r}$ , where $p_{ik}$,
$q_{ls}\in\mathcal{P}$, to obtain

$\displaystyle\left(\left(I_{2}\otimes\,M\right)A^{\otimes\,m}\nabla
N\right)\cdot\overrightarrow{n}$ $\displaystyle=$
$\displaystyle\left(\rho\left(a_{11}n_{x}+a_{21}n_{y}\right)p_{ik}+\rho\left(a_{12}n_{x}+a_{22}n_{y}\right)q_{ik}\right)_{i=1,\,k=1}^{p,\,r}$
$\displaystyle=$
$\displaystyle\left(\left(\left(\rho\,A\right)\cdot\overrightarrow{n}\right)\left(\begin{matrix}p_{ik}\\\
q_{ik}\end{matrix}\right)\right)_{i=1,\,k=1}^{p,\,r}.$

Since $\left(\nabla\mathbb{P}^{t}_{n}\right)_{n\geq 1}$ is a basis of
$\mathcal{P}\times\mathcal{P}$, for fixed $i$ and $k$,

$\left(\begin{matrix}p_{ik}\\\
q_{ik}\end{matrix}\right)=\sum_{l=0}^{z}\nabla\mathbb{P}_{l}^{t}C_{n,l}^{i,\,k},\;\;z={\rm
degree}\left(\begin{matrix}p_{ik}\\\ q_{ik}\end{matrix}\right).$

Therefore,

$\displaystyle\left(\rho\left(I_{2}\otimes\,M\right)A^{\otimes\,m}\nabla
N\right)\cdot\overrightarrow{n}$ $\displaystyle=$
$\displaystyle\left(\sum_{l=0}^{z}\left(\left(\rho\,A\right)\cdot\overrightarrow{n}\right)\nabla\mathbb{P}_{l}^{t}C_{n,l}^{i,\,k}\right)_{i=1,\,k=1}^{p,\,r}.$

Replace $\overrightarrow{n}$ by $\overrightarrow{n}_{j}$, multiply from the
left by the indicator function of $\partial\Omega_{j}$,
$1_{\partial\Omega_{j}}$, and use the fact that
$\left(\left(\rho\,A\right)\cdot\overrightarrow{n_{j}}\right)\nabla\mathbb{P}_{l}^{t}=\left(\rho\,A\nabla\mathbb{P}_{l}^{t}\right)\cdot\overrightarrow{n_{j}}$
as well as the hypothesis (18) to get

$\lim\limits_{j\rightarrow\infty}1_{\partial\Omega_{j}}\left(\rho\left(I_{2}\otimes\,M\right)A^{\otimes\,m}\nabla
N\right)\cdot\overrightarrow{n_{j}}=0.$ (20)

To end this proof, observe that
$\int_{\Omega}div\left(\rho\left(I_{2}\otimes\,M\right)A^{\otimes\,m}\nabla\,N\right)dxdy=\lim\limits_{j\rightarrow\infty}\int_{{\Omega}_{j}}div\left(\rho\left(I_{2}\otimes\,M\right)A^{\otimes\,m}\nabla\,N\right)dxdy.$
Since $\Omega_{j}$ is a connected open bounded subset of $\mathbb{R}^{2}$ it
has a piecewise smooth boundary. So, from the divergence theorem,
$\int_{{\Omega}_{j}}div\left(\rho\left(I_{2}\otimes\,M\right)A^{\otimes\,m}\nabla\,N\right)dxdy=\int_{\partial{\Omega}_{j}}\left(\rho\left(I_{2}\otimes\,M\right)A^{\otimes\,m}\nabla\,N\right)\cdot\overrightarrow{n_{j}}dxdy.$
Therefore

$\int_{\Omega}div\left(\rho\left(I_{2}\otimes\,M\right)A^{\otimes\,m}\nabla\,N\right)dxdy=\lim\limits_{j\rightarrow\infty}\int
1_{\partial{\Omega}_{j}}\left(\rho\left(I_{2}\otimes\,M\right)A^{\otimes\,m}\nabla\,N\right)\cdot\overrightarrow{n_{j}}dxdy$

and
$\left\|\int_{\Omega}div\left(\rho\left(I_{2}\otimes\,M\right)A^{\otimes\,m}\nabla\,N\right)dxdy\right\|_{max}\leq\lim\limits_{j\rightarrow\infty}\int\left\|1_{\partial{\Omega}_{j}}\left(\rho\left(I_{2}\otimes\,M\right)A^{\otimes\,m}\nabla\,N\right)\cdot\overrightarrow{n_{j}}\right\|_{max}dxdy$.
From (20),
$\lim\limits_{j\rightarrow\infty}\left\|1_{\partial{\Omega}_{j}}\left(\rho\left(I_{2}\otimes\,M\right)A^{\otimes\,m}\nabla\,N\right)\cdot\overrightarrow{n_{j}}\right\|_{max}=0$.
Therefore there exists $j_{0}\geq\,1$ such that

$\left\|1_{\partial{\Omega}_{j}}\left(\rho\left(I_{2}\otimes\,M\right)A^{\otimes\,m}\nabla\,N\right)\cdot\overrightarrow{n_{j}}\right\|_{max}\leq\left\|1_{\partial{\Omega}_{j_{0}}}\left(\rho\left(I_{2}\otimes\,M\right)A^{\otimes\,m}\nabla\,N\right)\cdot\overrightarrow{n_{j_{0}}}\right\|_{max}.$

Use Lebesgue’s dominated convergence theorem to have

$\lim\limits_{j\rightarrow\infty}\int\left\|1_{\partial{\Omega}_{j}}\left(\rho\left(I_{2}\otimes\,M\right)A^{\otimes\,m}\nabla\,N\right)\cdot\overrightarrow{n_{j}}\right\|_{max}dxdy=\int\lim\limits_{j\rightarrow\infty}\left\|1_{\partial{\Omega}_{j}}\left(\rho\left(I_{2}\otimes\,M\right)A^{\otimes\,m}\nabla\,N\right)\cdot\overrightarrow{n_{j}}\right\|_{max}dxdy=0.$

Thus
$\left\|\int_{\Omega}div\left(\rho\left(I_{2}\otimes\,M\right)A^{\otimes\,m}\nabla\,N\right)dxdy\right\|_{max}=0.$
∎

######  Proposition 4.

Let $A$ be a $(2,2)$-matrix functions.

  1. 1.

Consider $M\in\mathcal{M}_{p,2^{m-1}}\left(\mathcal{P}\right)$ and
$N\in\mathcal{M}_{2^{(m-1)},q}\left(\mathcal{P}\right)$. The following
relation holds

$div\left[\left(I_{2}\otimes\,M\right)A^{\otimes\,m}\nabla\,N\right]=Mdiv(A^{\otimes\,m}\nabla
N)+\left(\partial_{x}M,\partial_{y}M\right)\left(A^{\otimes\,m}\nabla
N\right).$ (21)

  2. 2.

Let $m$ be a positive integer and
$N\in\mathcal{M}_{2^{m},q}\left(\mathcal{P}\right)$. The following relation
holds

$div\left[\left(B\otimes\,A^{\otimes
m}\right)\nabla\,N\right]=A^{\otimes\,m}\left[\left(\left(B\nabla\right)\cdot\nabla\right)N\right]+div\left(B\otimes\,A^{\otimes\,m}\right)\nabla\,N.$
(22)

###### Proof.

To prove the first item, write

$\left(I_{2}\otimes\,M\right)A^{\otimes\,m}\nabla
N=\left(\begin{matrix}a_{11}MA^{\otimes\,m-1}\partial_{x}N+a_{12}MA^{\otimes\,m-1}\partial_{y}N\\\
a_{21}MA^{\otimes\,m-1}\partial_{x}N+a_{22}MA^{\otimes\,m-1}\partial_{y}N\end{matrix}\right),\,A=\left(\begin{matrix}a_{11}&a_{12}\\\
a_{21}&a_{22}\end{matrix}\right),$

and apply $div$ on both sides to obtain the result. As for the last item of
this proposition, apply $div$ on both sides of the relation

$\left(B\otimes\,A^{\otimes\,m}\right)=\left(\begin{matrix}A^{\otimes\,m}\left(b_{11}\partial_{x}\,N+b_{12}\partial_{y}\,N\right)\\\
A^{\otimes\,m}\left(b_{21}\partial_{x}\,N+b_{22}\partial_{y}\,N\right)\end{matrix}\right),\;\,B=\left(\begin{matrix}b_{11}&b_{12}\\\
b_{21}&b_{22}\end{matrix}\right),$

and use the usual product rule to obtain the result. ∎

######  Proposition 5.

Let $m$ and $n$ be two non-negative integers, $0\leq m\leq n$. Let $\Phi$ be a
$(2,\,2)$-matrix polynomial of total degree at most two. Let
$\psi_{i}^{(m)}(x,y)=\left(I_{2^{m}}\otimes
X^{t}\right)D_{i}^{(m)}+E_{i}^{(m)}$ with
$D_{i}^{(m)}\in\mathcal{M}_{2^{m+1},\,2^{m}}\left(\mathbb{R}\right)$ and
$E_{i}^{(m)}\in\mathcal{M}_{2^{m}}\left(\mathbb{R}\right)$, $i=1,2$ be two
$2^{m}$-matrix polynomials of total degree one. Then

$\left(\left(\Phi\nabla\right)\cdot\nabla\right)\mathbb{Q}_{n,m}+\psi_{1}^{(m)}\partial_{x}\mathbb{Q}_{n,m}+\psi_{2}^{(m)}\partial_{y}\mathbb{Q}_{n,m}=\left(I_{2^{m}}\otimes\,X_{n}^{t}\right)T_{n}^{(m)}G_{n,m,n}+\cdots,$
(23)

where

$T_{n}^{(m)}=\left[\left(L_{n-1}^{(m)*}\right)^{t}\left(\left(A_{1},2A_{2},A_{3}\right)\otimes\,I_{2^{m}(n-1)}\right)N_{n}^{(m)*}+\left(L_{n-1}^{(m)}\right)^{t}\left(\left(D_{1}^{(m)},D_{2}^{(m)}\right)\otimes\,I_{n}\right)N_{n}^{(m)}\right]$

with $A_{i}=\left(a_{i1},a_{i2},a_{i3}\right)^{t}$, $i=1,2,3$,
$D_{j}=\left(d_{j1},d_{j2}\right)^{t}$, $j=1,2$,

$\displaystyle L_{n}^{(m)*}$ $\displaystyle=$
$\displaystyle\left(\begin{matrix}I_{2^{m}}\otimes\,L_{n-1,1}L_{n,1}\\\
I_{2^{m}}\otimes\,L_{n-1,2}L_{n,1}\\\
I_{2^{m}}\otimes\,L_{n-1,2}L_{n,2}\end{matrix}\right),\;\;N_{n}^{(m)*}=\left(\begin{matrix}I_{2^{m}}\otimes\,N_{n-1,1}N_{n,1}\\\
I_{2^{m}}\otimes\,N_{n-1,2}N_{n,1}\\\
I_{2^{m}}\otimes\,N_{n-1,2}N_{n,2}\end{matrix}\right),$ $\displaystyle
L_{n}^{(m)}$ $\displaystyle=$
$\displaystyle\left(\begin{matrix}I_{2^{m}}\otimes\,L_{n,1}\\\
I_{2^{m}}\otimes\,L_{n,2}\end{matrix}\right)\;{\rm
and}\;N_{n}^{(m)}=\left(\begin{matrix}I_{2^{m}}\otimes\,N_{n,1}\\\
I_{2^{m}}\otimes\,N_{n,2}\end{matrix}\right),$

where $L_{n,j}$ and $N_{n,j}$ are matrices in (8) and (10) respectively.
$L_{0}^{(0)*}\equiv\,0$, $N_{0}^{(0)*}\equiv\,0$, $L_{n}^{(0)}\equiv\,L_{n}$
and $N_{n}^{(0)}\equiv\,N_{n}$. $G_{n,m,n}$ is the
$\left(2^{m}(n+1),n+m+1\right)$-matrix given by (17) and $T_{n}^{(m)}$ a
$({2^{m}(n+1)})$-matrix with $T_{n}^{(0)}\equiv\,T_{n}$.

###### Proof.

Observe that

$\left(\left(\Phi\nabla\right)\cdot\nabla\right)\,\mathbb{Q}_{n,m}+\psi\nabla\mathbb{Q}_{n,m}=\phi_{11}\partial_{x}^{2}\mathbb{Q}^{t}_{n,m}+2\phi_{12}\partial_{xy}^{2}\mathbb{Q}_{n,m}+\phi_{22}\partial_{y}^{2}\mathbb{Q}_{n,m}+\psi_{1}^{(m)}\partial_{x}\mathbb{Q}_{n,m}+\psi_{2}^{(m)}\partial_{y}\mathbb{Q}_{n,m}$

write
$\mathbb{Q}_{n,m}=\left(I_{2^{m}}\otimes\,X^{t}_{n}\right)G_{n,m,n}+\dots$ and
use the partial derivatives of $\left(I_{2^{m}}\otimes\,X^{t}_{n}\right)$
given by (14) to obtain

$\displaystyle\left(\left(\Phi\nabla\right)\cdot\nabla\right)\,\mathbb{Q}_{n,m}+\psi_{1}^{(m)}\partial_{x}\mathbb{Q}_{n,m}+\psi_{2}^{(m)}\partial_{y}\mathbb{Q}_{n,m}$
(24)
$\displaystyle=\left[\phi_{11}\left(I_{2^{m}}\otimes\,X^{t}_{n-2}\right)\left(I_{2^{m}}\otimes\,N_{n-1,1}N_{n,1}\right)+2\phi_{12}\left(I_{2^{m}}\otimes\,X^{t}_{n-2}\right)\left(I_{2^{m}}\otimes\,N_{n-1,1}N_{n,2}\right)\right]G_{n,m,n}$
$\displaystyle+\left[\phi_{22}\left(I_{2^{m}}\otimes\,X^{t}_{n-2}\right)\left(I_{2^{m}}\otimes\,N_{n-1,2}N_{n,2}\right)+\psi_{1}^{(m)}\left(I_{2^{m}}\otimes\,X^{t}_{n-1}\right)\left(I_{2^{m}}\otimes\,N_{n,1}\right)\right]G_{n,m,n}$
$\displaystyle+\left[\psi_{2}^{(m)}\left(I_{2^{m}}\otimes\,X^{t}_{n-1}\right)\left(I_{2^{m}}\otimes\,N_{n,2}\right)\right]G_{n,m,n}+\dots.$

Since we assumed there is no restriction on the shape of the polynomials
$\phi_{ij},\,i,j=1,2$ and $\psi_{k},\,k=1,2$, expand $\phi_{11}$ as

$\phi_{11}(x,y)=a_{11}x^{2}+a_{12}xy+a_{13}y^{2}+b_{11}x+b_{12}y+c_{1}=X_{2}^{t}A_{1}+X^{t}B_{1}+c_{1},$

where $A_{1}=\left(a_{11},\,a_{12},\,a_{13}\right)$,
$B_{1}=\left(b_{11},\,b_{12}\right)$ and take into account (13a)-(13d) to
obtain

$\displaystyle\phi_{11}(x,y)\left(I_{2^{m}}\otimes\,X^{t}_{n-2}\right)$
$\displaystyle=\left(I_{2^{m}}\otimes\,X^{t}_{n}\right)\left[a_{11}\left(I_{2^{m}}\otimes\,L_{n-1,1}^{t}L_{n-2,1}^{t}\right)+a_{12}\left(I_{2^{m}}\otimes\,L_{n-1,1}^{t}L_{n-2,2}^{t}\right)+a_{13}\left(I_{2^{m}}\otimes\,L_{n-1,1}^{t}L_{n-2,2}^{t}\right)\right]$
$\displaystyle+\left(I_{2^{m}}\otimes\,X_{n-1}^{t}\right)\left[b_{11}\left(I_{2^{m}}\otimes\,L_{n-2,1}^{t}\right)+b_{12}\left(I_{2^{m}}\otimes\,L_{n-2,2}^{t}\right)\right]+c_{1}\left(I_{2^{m}}\otimes\,X_{n-1}^{t}\right)$
$\displaystyle=\left(I_{2^{m}}\otimes\,X^{t}_{n}\right)\left(L_{n-1}^{(m)*}\right)^{t}\left(A_{1}\otimes\,I_{2^{m}(n-1)}\right)+\left(I_{2^{m}}\otimes\,X^{t}_{n-1}\right)\left(L_{n-2}^{(m)}\right)^{t}\left(B_{1}\otimes\,I_{2^{m}(n-1)}\right)+c_{1}\left(I_{2^{m}}\otimes\,X^{t}_{n-2}\right).$

Replace $\left(A_{1},\,B_{1},c_{1}\right)$ by
$\left(A_{i},\,B_{i},c_{i}\right)$, $i=2,\,3$, to get
$\phi_{12}(x,y)\left(I_{2^{m}}\otimes\,X^{t}_{n-2}\right)$ and
$\phi_{22}(x,y)\left(I_{2^{m}}\otimes\,X^{t}_{n-2}\right)$ in terms of
$\left(I_{2^{m}}\otimes\,X^{t}_{n}\right)$,
$\left(I_{2^{m}}\otimes\,X^{t}_{n-1}\right)$ and
$\left(I_{2^{m}}\otimes\,X^{t}_{n-2}\right)$. Use (13e) to have

$\displaystyle\psi_{1}^{(m)}\left(I_{2^{m}}\otimes\,X^{t}_{n-1}\right)$
$\displaystyle=\left(I_{2^{m}}\otimes\,X^{t}\right)D_{1}^{(m)}\left(I_{2^{m}}\otimes\,X^{t}_{n-1}\right)+E_{1}^{(m)}\left(I_{2^{m}}\otimes\,X^{t}_{n-1}\right)$
$\displaystyle=\left(I_{2^{m}}\otimes\,X^{t}_{n}\right)\left(I_{2^{m}}\otimes\,L_{n-1}\right)\left(D_{1}^{(m)}\otimes\,I_{n}\right)+E_{1}^{(m)}\left(I_{2^{m}}\otimes\,X^{t}_{n-1}\right)$

Replace $\left(D_{1}^{(m)},\,E_{1}^{(m)}\right)$ by
$\left(D_{2}^{(m)},\,E_{2}^{(m)}\right)$ to expand
$\psi_{1}^{(m)}(x,y)\left(I_{2^{m}}\otimes\,X^{t}_{n-1}\right)$ in terms of
$\left(I_{2^{m}}\otimes\,X^{t}_{n}\right)$ and
$\left(I_{2^{m}}\otimes\,X^{t}_{n-1}\right)$. Finally, substitute into (24) to
obtain

$\displaystyle\left(\Phi\nabla\right)\cdot\nabla\,\nabla\mathbb{Q}_{n,m}+\psi_{1}^{(m)}\partial_{x}\mathbb{Q}_{n,m}+\psi_{2}^{(m)}\partial_{y}\mathbb{Q}_{n,m}$
$\displaystyle=\left(I_{2^{m}}\otimes\,X^{t}_{n}\right)\left(L_{n-1}^{(m)*}\right)^{t}\left[\left(A_{1}\otimes\,I_{2^{m}(n-1)}\right)\left(I_{2^{m}}\otimes\,N_{n-1,1}N_{n,1}\right)+2\left(A_{2}\otimes\,I_{2^{m}(n-1)}\right)\left(I_{2^{m}}\otimes\,N_{n-1,1}N_{n,2}\right)\right]$
$\displaystyle\times\,G_{n,m,n}+\left(I_{2^{m}}\otimes\,X^{t}_{n}\right)\left(L_{n-1}^{(m)*}\right)^{t}\left[\left(A_{3}\otimes\,I_{2^{m}(n-1)}\right)\left(I_{2^{m}}\otimes\,N_{n-1,2}N_{n,2}\right)\right]G_{n,m,n}$
$\displaystyle+\left(I_{2^{m}}\otimes\,L_{n-1}^{t}\right)\left[\left(D_{1}^{(m)}\otimes\,I_{n}\right)\left(I_{2^{m}}\otimes\,N_{n,1}\right)+\left(D_{2}^{(m)}\otimes\,I_{n}\right)\left(I_{2^{m}}\otimes\,N_{n,2}\right)\right]G_{n,m,n}+\dots$
$\displaystyle=\left(I_{2^{m}}\otimes\,X^{t}_{n}\right)T_{n}^{(m)}G_{n,m,n}+\dots\,.$

with
$T_{n}^{(m)}=\left[\left(L_{n-1}^{(m)*}\right)^{t}\left(\left(A_{1},\,2A_{2},\,A_{3}\right)\otimes\,I_{2^{m}(n-1)}\right)N_{n}^{(m)*}+\left(I_{2^{m}}\otimes\,L_{n-1}^{t}\right)\left(\left(D_{1}^{(m)},D_{2}^{(m)}\right)\otimes\,I_{n}\right)N_{n}^{(m)}\right]$.
∎

######  Proposition 6.

Let $\left\\{\mathbb{P}_{n}\right\\}_{n}$ be a family of vector polynomials,
orthogonal with respect to the weight function $\rho$, satisfying the Pearson
equation $div\left(\rho\Phi\right)=\rho\left(\psi_{1},\psi_{2}\right)$ with
the Neumann boundary condition (18) where
$\Phi\in\mathcal{M}_{2}\left(\mathcal{P}\right)$ and
$\psi_{i}(X)=X^{t}D_{i}+E_{i}$ is a polynomial of total degree one, $i=1,\,2$,
with $det(D_{1},\,D_{2})\neq 0$ then

${\rm
degree}\left(\left(\left(\Phi\nabla\right)\cdot\nabla\right)\,\mathbb{P}^{t}_{n}+\Psi\nabla\mathbb{P}_{n}^{t}\right)=n,\;\;\forall
n,\;n\geq 1.$

###### Proof.

For the proof, we are going to consider two cases. First case, if degree
$\Phi$ =2.
Take $M=\mathbb{P}_{n}$, $m=1$, $A=\rho\Phi$ and $N=\mathbb{P}_{j}^{t}$ into
(21) to obtain

$div\left[\left(I_{2}\otimes\mathbb{P}_{n}\right)\rho\Phi\nabla\mathbb{P}_{j}^{t}\right]=\mathbb{P}_{n}div\left(\rho\Phi\nabla\mathbb{P}_{j}^{t}\right)+\left(\nabla\mathbb{P}_{n}^{t}\right)^{t}\rho\phi\nabla\mathbb{P}_{j}^{t}$

Integrate both sides on $\Omega$ and use the boundary condition (18), (22)
with $m=0$, $B=\rho\Phi$ and $N=\mathbb{P}_{j}^{t}$ as well as the hypothesis
$div\left(\rho\Phi\right)=\rho\left(\psi_{1},\psi_{2}\right)$ to obtain

$\displaystyle-\int_{\Omega}\left(\nabla\mathbb{P}_{n}^{t}\right)^{t}\rho\Phi\nabla\mathbb{P}_{j}^{t}dxdy=\int_{\Omega}\mathbb{P}_{n}\left[\left(\rho\phi\nabla\right)\cdot\nabla\mathbb{P}_{j}^{t}+\rho\left(\psi_{1},\psi_{2}\right)\nabla\mathbb{P}_{j}^{t}\right]dxdy.$

Take into account (23) (with $m=0$ and $n=j$), write
$X_{j}^{t}=\mathbb{P}_{j}^{t}+\dots$ and use the orthogonality of the system
$\\{\mathbb{P}_{n}^{t}\\}_{n}$ with respect to $\rho$ to obtain

$\displaystyle-\int_{\Omega}\left(\nabla\mathbb{P}_{n}^{t}\right)^{t}\rho\phi\nabla\mathbb{P}_{j}^{t}dxdy=\int_{\Omega}\mathbb{P}_{n}\mathbb{P}_{n}^{t}T_{n}\delta_{n,j}dxdy=S_{n}T_{n}\delta_{n,j},$

where $\delta_{n,j}$ is the Kronecker symbol. Therefore, multiplying from the
left hand side the expansion

$\left(\begin{matrix}\mathbb{P}_{n+2}^{t}\\\
\mathbb{P}_{n+2}^{t}\end{matrix}\right)=\sum_{k=0}^{n+2}\nabla\mathbb{P}_{k+1}^{t}A_{k+1}$

by $\left(\nabla\mathbb{P}_{n+1}^{t}\right)^{t}\rho\Phi$ and integrating both
sides on $\Omega$, we obtain

$\int_{\Omega}\left(\nabla\mathbb{P}_{n+1}^{t}\right)^{t}\rho\Phi\left(\begin{matrix}\mathbb{P}_{n+2}^{t}\\\
\mathbb{P}_{n+2}^{t}\end{matrix}\right)=-S_{n+1}T_{n+1}A_{n+1}.$ (25)

Since

$\displaystyle\left(\nabla\mathbb{P}_{n+1}^{t}\right)^{t}\rho\Phi\left(\begin{matrix}\mathbb{P}_{n+2}^{t}\\\
\mathbb{P}_{n+2}^{t}\end{matrix}\right)$ $\displaystyle=$
$\displaystyle\left(\phi_{11}\partial_{x}\mathbb{P}_{n+1}+\phi_{12}\partial_{y}\mathbb{P}_{n+1}+\phi_{21}\partial_{x}\mathbb{P}_{n+1}+\phi_{22}\partial_{y}\mathbb{P}_{n+1}\right)\mathbb{P}_{n+2}^{t}$
$\displaystyle=$ $\displaystyle C_{n}X_{n+2}\mathbb{P}^{t}_{n+2}+\dots,$

and $\Phi$ is of total degree $2$, $C_{n}$ is a $(n+2,n+3)$-matrix different
from $0$. Therefore the left hand side of (25) is different from $0$. Thus
$T_{n+1}\neq\,0$.
If degree $\Phi<2$, then

$T_{n}=\left(L_{n-1}\right)^{t}\left(\left(D_{1},D_{2}\right)\otimes\,I_{n}\right)N_{n}$

Observing that
$\left(D_{1},D_{2}\right))\otimes\,I_{n}\in\mathcal{M}_{2n}(\mathbb{R})$,
$N_{n}\in\mathcal{M}_{2n,n+1}(\mathbb{R})$ and using the Sylvester inequality,
we obtain

$\displaystyle
Rank[\left(D_{1},D_{2}\right))\otimes\,I_{n}]+Rank(N_{n})-2n\leq$
$\displaystyle
Rank[\left(\left(D_{1},D_{2}\right)\otimes\,I_{n}\right)\,N_{n}]$
$\displaystyle\leq$ $\displaystyle
min\left(Rank[\left(D_{1},D_{2}\right))\otimes\,I_{n}],\,Rank(N_{n})\right).$

Moreover,
$Rank\left(\left(D_{1},D_{2}\right)\otimes\,I_{n}\right)=Rank[\left(D_{1},D_{2}\right))]Rank(I_{n})$
and $Rank[\left(D_{1},D_{2}\right)]=2$ (for $det\left(D_{1},D_{2}\right)\neq
0$). Therefore $Rank\left(\left(D_{1},D_{2}\right)\otimes\,I_{n}\right)=2n$
and $Rank[\left(\left(D_{1},D_{2}\right)\otimes\,I_{n}\right)\,N_{n}]=n+1$.
Taking $A=\left(L_{n-1}\right)^{t}$ and
$B=\left(\left(D_{1},D_{2}\right)\otimes\,I_{n}\right)N_{n}$ and using again
the Sylvester inequality, we have

$\displaystyle n+1=Rank(A^{t})+Rank(B^{t})-(n+1)\leq
Rank(B^{t}A^{t})\leq\,min(Rank(B^{t}),Rank(A^{t}))=n+1$

So, $Rank(AB)=n+1$. i.e. $Rank(T_{n})=n+1$. Therefore $T_{n}$ is invertible,
for $T_{n}$ is a $(n+1)$-matrix.

∎

We now consider a Hilbert space that will be useful for the proof of main
theorem of this work. For that purpose we consider the weighted Lebesque space

$L^{2}(\Omega,\rho)=\left\\{u:\Omega\to\mathbb{R}\,:\text{$u$ measurable and
}\;\int_{\Omega}u^{2}\rho\,dxdy<+\infty\right\\}$

equipped with the scalar product $\langle
v,u\rangle_{\rho}=\int_{\Omega}uv\rho\,dxdy.$ Clearly the corresponding norm
is $\|u\|_{\rho}=\left(\int_{\Omega}u^{2}\rho\,dxdy\right)^{1\over 2}$.

######  Theorem 7.

[15] Let $\Omega$ be a simply connected subset of $\mathbb{R}^{2}$ and $\rho$
be a weight function on $\Omega$. Let $\Phi$ be a symetric matrix of
$\mathcal{M}_{2}\left(\mathcal{P}_{2}\right)$, positive defined on $\Omega$
such that the Neumann type boundary condition (18) is satisfied. Then
$\mathcal{P}$ is dense in $L^{2}(\Omega,\rho)$.

##  3 Characterization theorem of classical orthogonal polynomials in two
variables

######  Theorem 8.

Let $\Omega$ be a simply connected open subset of $\mathbb{R}^{2}$, $\rho$ a
weight function on $\Omega$ and
$\varPhi\in\mathcal{M}_{2}\left(\mathcal{P}_{2}\right)$ a symmetric matrix
$\varPhi=\left(\begin{matrix}\phi_{1,1}&\phi_{1,2}\\\
\phi_{2,1}&\phi_{2,2}\end{matrix}\right)$ such that the Neumann type boundary
condition

$\lim\limits_{j\rightarrow\infty}1_{\partial\Omega_{j}}\left(\rho\Phi\,\nabla\,u\right)\cdot\overrightarrow{n_{j}}=0,\;\Omega_{j}=\Omega\cap
B(O,j),\;{\rm for~{}all}\;u\in\mathcal{M}_{1,n}\left(\mathcal{P}\right)$ (26)

is satisfied and

$\displaystyle\phi_{1,1}\partial_{x}\varPhi+\phi_{2,1}\partial_{y}\varPhi=\varPhi\nabla\left(\phi_{1,1},\phi_{2,1}\right),$
(27a)
$\displaystyle\phi_{1,2}\partial_{x}\varPhi+\phi_{2,2}\partial_{y}\varPhi=\varPhi\nabla\left(\phi_{1,2},\phi_{2,2}\right).$
(27b)

Let $\left(\mathbb{P}_{n}\right)_{n}$,
$\mathbb{P}_{n}\in\mathcal{M}_{n+1,1}\left(\mathcal{P}\right)$ be the family
of vector polynomials, orthogonal with respect to the weight $\rho$ and
$\left(\mathbb{Q}_{n,m}\right)_{n}$,
$\mathbb{Q}_{n,m}\in\mathcal{M}_{2^{m},n+1}\left(\mathcal{P}\right)$, the
family of monic polynomials of degree $n$

$\mathbb{Q}_{n,m}=\nabla\left(\mathbb{Q}_{n+1,m-1}\right),\;\;\;\mathbb{Q}_{n,0}=\mathbb{P}_{n}^{t}.$

The following properties are equivalent:

  1. a.

There exist two polynomials $\psi_{i}(x,y)=X^{t}D_{i}+E_{i},\,i=1,2$, such
that $det\left(D_{1},\,D_{2}\right)\neq\,0$ and

$div(\rho\varPhi)=\rho\left(\psi_{1},\,\psi_{2}\right).$

  2. b.

There exist two polynomials $\psi_{i}(x,y)=X^{t}D_{i}+E_{i},\,i=1,2$, with
$det\left(D_{1},\,D_{2}\right)\neq\,0$ such that for all positive integers
$m$, $\left\\{\mathbb{Q}_{n,m}\right\\}_{n}$ is orthogonal with respect to
$\rho_{m}=\rho\phi^{\otimes\,m}$ and

$div(\rho_{m}\otimes\varPhi)=\rho_{m}\left(\psi_{1}^{(m)},\psi^{(m)}_{2}\right),$
(28)

where $\psi_{i}^{(m)},\,i=1,2$, are $2^{m}$-matrix polynomials defined by

$\displaystyle\psi^{(m)}_{1}=I_{2}\otimes\psi^{(m-1)}_{1}+\nabla\left(\phi_{1,1},\phi_{2,1}\right)\otimes
I_{2^{m-1}},$ (29a)
$\displaystyle\psi^{(m)}_{2}=I_{2}\otimes\psi^{(m-1)}_{2}+\nabla\left(\phi_{1,2},\phi_{2,2}\right)\otimes
I_{2^{m-1}}$ (29b)

with $\psi_{1}^{(0)}=\psi_{1}$ and $\psi_{2}^{(0)}=\psi_{2}$.

  3. c.

For fixed integers $n\geq 0$ and $m\geq 0$, $\mathbb{Q}_{n,m}$ satisfies the
second order partial differential equation

$\left(\left(\varPhi\nabla\right)\cdot\nabla\right)\mathbb{Q}_{n,m}+\psi_{1}^{(m)}\partial_{x}\mathbb{Q}_{n,m}+\psi_{2}^{(m)}\partial_{y}\mathbb{Q}_{n,m}+\mathbb{Q}_{n,m}\Lambda_{n+m,m}=0,$
(30)

where the matrix polynomials $\psi_{1}^{(m)}$, $\psi_{2}^{(m)}$ are defined in
(29a) and (29b), $\Lambda_{n+m,m}$ is the $n+m+1$-matrix defined by

$\displaystyle G_{n,m,n}\Lambda_{n+m,m}$
$\displaystyle=-\left[\left(L_{n-1}^{(m)*}\right)^{t}\left(\left(A_{1},\,2A_{2},\,A_{3}\right)\otimes\,I_{2^{m}(n-1)}\right)N_{n}^{(m)*}+\left(L_{n-1}^{(m)}\right)^{t}\left(\left(D_{1}^{(m)},\,D_{2}^{(m)}\right)\otimes\,I_{n}\right)N_{n}^{(m)}\right]G_{n,m,n},$

where $G_{n,m,n}$ is the leading coefficient of $\mathbb{Q}_{n,m}$ given in
(17).

  4. d.

For an integer $n\geq 1$ the following Rodrigues formula holds

$\displaystyle\mathbb{P}^{t}_{n}=\frac{(-1)^{n}}{\rho}div^{(n)}\left[\rho\varPhi^{\otimes\,n}\right]R_{n},$

where
$R_{n}=\nabla^{(n)}\mathbb{P}_{n}^{t}\prod_{j=0}^{n-1}\left(\Lambda_{n,j}\right)^{-1}$
and $\Lambda_{n,j},\,j=0..n-1$ are the matrices in (30).

  5. e.

There exist $A^{n,m}_{n+1}$, $A^{n,m}_{n}$ and $A^{n,m}_{n-1}$, where
$A^{n,m}_{p}=\left(\begin{matrix}A^{n,m}_{p,1}&\\\
A^{n,m}_{p,2}\end{matrix}\right)$, with
$A^{n,m}_{p,i}\in\mathcal{M}_{p+1,\,n+1}\left(\mathbb{R}\right)$, $i=1$ ,$2$
and $p=n-1\,,n\,,n+1$ such that

$\left(\varPhi\otimes\,I_{2^{m}}\right)\mathbb{Q}_{n-1,m+1}=(I_{2}\otimes\mathbb{Q}_{n+1,m})A^{n,m}_{n+1}+(I_{2}\otimes\mathbb{Q}_{n,m})A^{n,m}_{n}+(I_{2}\otimes\mathbb{Q}_{n-1,m})A_{n-1}^{n,m}$
(31)

and $A^{n,m}_{n-1}$ invertible.

The following lemmas will be useful in the proof of the Theorem 8.

######  Lemma 9.

Let $D_{1}=\left(d_{11},\,d_{12}\right)^{t}$ and
$D_{2}=\left(d_{21},\,d_{22}\right)^{t}$ be two (2,1)-vectors. The determinant
of the $2m$-matrix $\left(I_{m}\otimes\,D_{1},\;I_{m}\otimes\,D_{2}\right)$
,$m\geq\,1$, is given by

$det\left(I_{m}\otimes\,D_{1},\;I_{m}\otimes\,D_{2}\right)=\left(-1\right)^{\left[m\over
2\right]}\left[det\left(D_{1},\;D_{2}\right)\right]^{m}.$

###### Proof.

The proof is done by induction on $m$. ∎

######  Lemma 10.

Polynomials defined by (29a) and (29b) can be written as

$\psi_{i}^{(m)}(x,y)=\left(I_{2^{m}}\otimes\,X^{t}\right)D_{i}^{(m)}+E_{i}^{(m)},\;i=1,2,$
(32)

where

$\displaystyle D_{i}^{(m+1)}$
$\displaystyle=H_{i}^{(m)}+I_{2}\otimes\,D_{i}^{(m)},\;\;D_{i}^{(0)}\equiv\,D_{i},$
$\displaystyle E_{i}^{(m+1)}$
$\displaystyle=K_{i}^{(m)}+I_{2}\otimes\,E_{i}^{(m)},\;\;E_{i}^{(0)}\equiv\,E_{i}$

with

$\displaystyle
H_{i}^{(m+1)}=\left(\begin{matrix}I_{2^{m}}\otimes\,\left(N_{2,1}A_{i}\right)&I_{2^{m}}\otimes\,\left(N_{2,1}A_{i+1}\right)\\\
I_{2^{m}}\otimes\,\left(N_{2,2}A_{i}\right)&I_{2^{m}}\otimes\,\left(N_{2,2}A_{i+1}\right)\end{matrix}\right),\,K_{i}^{(m+1)}=\left(\begin{matrix}\left(N_{1,1}B_{i}\right)I_{2^{m}}&\left(N_{1,1}B_{i+1}\right)I_{2^{m}}\\\
\left(N_{1,2}B_{i}\right)I_{2^{m}}&\left(N_{1,2}B_{i+1}\right)I_{2^{m}}\end{matrix}\right),$

$m\geq\,1.$

###### Proof.

We show by induction on $m$ that (32) is satisfied for $m\geq 0$. It is
obvious that (32) is satisfied for $m=0$. Assume that they are satisfied up to
a fixed integer $m>0$. We use the relation (29a) to obtain

$\displaystyle\psi_{1}^{(m+1)}$ $\displaystyle=$
$\displaystyle\left(\begin{matrix}\partial_{x}\phi_{1,1}I_{2^{m}}+\psi_{1}^{m}&\partial_{x}\phi_{2,1}I_{2^{m}}\\\
\partial_{y}\phi_{1,1}I_{2^{m}}&\partial_{y}\phi_{2,1}I_{2^{m}}+\psi_{1}^{(m)}\end{matrix}\right).$

Since $\phi_{11}$ can be written as
$\phi_{11}(x,y)=X_{2}^{t}A_{1}+X^{t}B_{1}+c_{1}$, using (9), the action of
$\partial_{x}$ on $\phi_{11}$ is
$\partial_{x}\,\phi_{11}(x,y)=X^{t}\left(N_{21}A_{1}\right)+N_{11}B_{1}$.
Observing that
$X^{t}\left(N_{21}A_{1}\right)I_{2^{m}}=\left(I_{2^{m}}\otimes\,X^{t}\right)\left(I_{2^{m}}\otimes\,\left(N_{21}A_{1}\right)\right)$
and using the fact that (32) is satisfied at order $m$ we obtain

$\partial_{x}\phi_{11}(x,y)I_{2^{m}}+\psi_{1}^{(m)}(x,y)=\left(I_{2^{m}}\otimes\,X^{t}\right)\left[I_{2^{m}}\otimes\,\left(N_{21}A_{1}\right)+D_{1}^{(m)}\right]+\left(N_{11}B_{1}\right)I_{2^{m}}+E_{1}^{(m)}.$

In a similar way, we write $\partial_{y}\phi_{11}(x,y)I_{2^{m}}$,
$\partial_{x}\phi_{21}(x,y)I_{2^{m}}$ and
$\partial_{y}\phi_{21}(x,y)I_{2^{m}}+\psi_{1}^{(m)}$ in terms of
$I_{2^{m}}\otimes X^{t}$ and take into account into the expression of
$\psi_{1}^{(m+1)}$ given previously to get

$\psi_{1}^{(m+1)}(x,y)=\left(I_{2^{m+1}}\otimes\,X^{t}\right)D_{1}^{(m+1)}+E_{1}^{(m+1)},$

where

$\displaystyle D_{1}^{(m+1)}$
$\displaystyle=\left(\begin{matrix}I_{2^{m}}\otimes\,\left(N_{2,1}A_{1}\right)+D_{1}^{(m)}&I_{2^{m}}\otimes\,\left(N_{2,1}A_{2}\right)\\\
I_{2^{m}}\otimes\,\left(N_{2,2}A_{1}\right)&I_{2^{m}}\otimes\,\left(N_{2,2}A_{2}\right)+D_{1}^{(m)}\end{matrix}\right)=H_{1}^{(m)}+I_{2}\otimes\,D_{1}^{(m)},$
$\displaystyle E_{1}^{(m+1)}$
$\displaystyle=\left(\begin{matrix}I_{2^{m}}\otimes\,\left(N_{1,1}B_{1}\right)+E_{1}^{(m)}&I_{2^{m}}\otimes\,\left(N_{1,1}B_{2}\right)\\\
I_{2^{m}}\otimes\,\left(N_{1,2}B_{1}\right)&I_{2^{m}}\otimes\,\left(N_{1,2}B_{2}\right)+D_{1}^{(m)}\end{matrix}\right)=K_{1}^{(m)}+I_{2}\otimes\,E_{1}^{(m)}.$

In a similar way, we obtain $\psi_{2}^{(m+1)}(x,y)$. ∎

######  Lemma 11.

Let $\left\\{\mathbb{Q}_{n,m}\right\\}_{n}$ be a family of vector polynomials,
orthogonal with respect to the weight function
$\rho_{m}=\rho\Phi^{\otimes\,m}$, satisfying the Pearson equation
$div\left(\rho_{m}\otimes\,\Phi\right)=\rho_{m}\left(\psi_{1}^{(m)},\psi_{2}^{(m)}\right)$
with the Neumann boundary condition (18) where
$\Phi\in\mathcal{M}_{2}\left(\mathcal{P}\right)$ and
$\psi_{i}^{(m)}(X)=X^{t}D_{i}^{(m)}+E_{i}^{(m)}$ , $i=1,\,2$, are polynomials
given by (29a)-(29b) with $det(D_{1},\,D_{2})\neq 0$ then

${\rm
degree}\left(\left(\left(\Phi\nabla\right)\cdot\nabla\right)\,\mathbb{Q}_{n,m}+\psi_{1}^{(m)}\partial_{x}\mathbb{Q}_{n,m}+\psi_{2}^{(m)}\partial_{y}\mathbb{Q}_{n,m}\right)=n,\;\;\forall
n,\;n\geq 1.$ (33)

###### Proof.

If $degree(\Phi)=2$, follow the first part of the proof of Proposition 6 to
obtain (33). If $degree(\Phi)\leq 2$, use Lemma 10 to have
$D_{i}^{(m)}=I_{2}\otimes\,D_{i}^{(m-1)}$, $i=1,2$ and iterate to obtain
$D_{i}^{(m)}=I_{2^{m}}\otimes\,D_{i}$, $i=1,2$. Next taking into account Lemma
9 yields
$det\left(D_{1}^{(m)},D_{2}^{(m)}\right)=\left(-1\right)^{\left[2^{m-1}\right]}\left[det\left(D_{1},\;D_{2}\right)\right]^{2^{m}}\neq\,0$,
for $det\left(D_{1},\,D_{2}\right)\neq\,0$. Finally follow the method
described in the second part of the proof of the Proposition 6, to get
$Rank\left(T_{n}^{(m)}\right)=2^{m}(n+1)$. Thus (33) follows. ∎

###### Proof of Theorem 8.

We organize the proof in the following scheme:

  1. Step 1

$(a)\Rightarrow\,(b)\Rightarrow\,(c)\Rightarrow(a)$ which is equivalent to
$(a)\Leftrightarrow\,(b)\Leftrightarrow\,(c)$.

  2. Step 2

$(b)$ and $(c)\Rightarrow\,(d)\Rightarrow\,(a)$ which taking into account Step
1 is equivalent to $(d)\Leftrightarrow(a)$.

  3. Step 3

$(b)\Rightarrow\,(e)\Rightarrow\,(a)$ which taking into account Step 1 is
equivalent to $(e)\Leftrightarrow\,(a)$.

  1. Step 1

$(a)\Rightarrow\,(b)\Rightarrow\,(c)\Rightarrow(a)$, which is equivalent to
$(a)\Leftrightarrow\,(b)\Leftrightarrow\,(c)$.

[Step 1.1] $(a)\Rightarrow\,(b)$
We assume that property (a) is satisfied, and we show by induction on $m$ that
(b) is satisfied for $m\geq 0$.
Let us show by induction on $m$ that (28) is satisfied for $m\geq 0$. Assume
that (28) is satisfied for a fixed integer $m\geq\,1$. Let us first observe
from the associativity of $\otimes$ that

$\rho_{m}\otimes\varPhi=\rho\varPhi\otimes\varPhi^{\otimes\,m}=\left(\begin{matrix}\rho\phi_{1,1}\varPhi^{\otimes\,m}&\rho\phi_{1,2}\varPhi^{\otimes\,m}\\\
\rho\phi_{2,1}\varPhi^{\otimes\,m}&\rho\phi_{2,2}\varPhi^{\otimes\,m}\end{matrix}\right).$
(34)

Therefore, the induction hypothesis
$div(\rho_{m}\otimes\varPhi)=\rho_{m}(\psi_{1}^{(m)},\,\psi_{2}^{(m)})$ is
equivalent to

$\displaystyle\partial_{x}\left(\rho\phi_{1,1}\varPhi^{\otimes\,m}\right)+\partial_{y}\left(\rho\phi_{2,1}\varPhi^{\otimes\,m}\right)=\rho_{m}\psi_{1}^{(m)},$
(35a)
$\displaystyle\partial_{x}\left(\rho\phi_{1,2}\varPhi^{\otimes\,m}\right)+\partial_{y}\left(\rho\phi_{2,2}\varPhi^{\otimes\,m}\right)=\rho_{m}\psi_{2}^{(m)}.$
(35b)

Replacing $m$ by $m+1$ in (34) and applying the divergence operator on both
sides we obtain

$div(\rho_{m+1}\otimes\Phi)=\left(\partial_{x}\left(\rho\phi_{1,1}\Phi^{\otimes\,m+1}\right)+\partial_{y}\left(\rho\phi_{2,1}\Phi^{\otimes\,m+1}\right),\partial_{x}\left(\rho\phi_{1,2}\Phi^{\otimes\,m+1}\right)+\partial_{y}\left(\rho\phi_{2,2}\varPhi^{\otimes\,m+1}\right)\right).$

Observing that for a polynomial $p\in\mathcal{P}$,
$p\Phi^{\otimes\,m+1}=\Phi\otimes p\varPhi^{\otimes\,m}$, taking
$p=\phi_{1,1}$ and using the product rule (16) with $A=\varPhi$ and
$B=\rho\phi_{1,1}\Phi^{\otimes\,m}$, we have
$\partial_{x}\left(\rho\phi_{1,1}\varPhi^{\otimes\,m+1}\right)=\varPhi\otimes\partial_{x}\left(\rho\phi_{1,1}\Phi^{\otimes\,m}\right)+\rho\phi_{1,1}\partial_{x}\varPhi\otimes\phi^{\otimes\,m}.$
In a similar way, we have
$\partial_{y}\left(\rho\phi_{2,1}\Phi^{\otimes\,m+1}\right)=\Phi\otimes\partial_{x}\left(\rho\phi_{2,1}\Phi^{\otimes\,m}\right)+\rho\phi_{2,1}\partial_{x}\varPhi\otimes\Phi^{\otimes
m}.$ Adding both and using (27a) and ( 35a), we obtain

$\displaystyle\partial_{x}\left(\rho\phi_{1,1}\varPhi^{\otimes\,m+1}\right)+\partial_{y}\left(\rho\phi_{2,1}\Phi^{\otimes\,m+1}\right)=\Phi\otimes\left(\rho\phi^{\otimes
m}\psi_{1}^{(m)}\right)+\left(\Phi\nabla(\phi_{1,1},\phi_{2,1})\right)\otimes\rho\phi^{\otimes
m}$

which can be rewritten as

$\displaystyle\partial_{x}\left(\rho\phi_{1,1}\Phi^{\otimes\,m+1}\right)+\partial_{y}\left(\rho\phi_{1,1}\varPhi^{\otimes\,m+1}\right)=\left(\Phi\,I_{2}\right)\otimes\left(\rho\phi^{\otimes
m}\psi_{1}^{(m)}\right)+\left(\Phi\nabla(\phi_{1,1},\phi_{2,1})\right)\otimes\left(\rho\phi^{\otimes
m}I_{2^{m}}\right)$

for $\Phi$ and $\Phi^{\otimes m}$ are $2$ and $2^{m}$ matrices. Taking
consecutively $\left(A,B,C,D\right)=\left(\Phi,\,\rho\phi^{\otimes
m},\,I_{2},\psi_{1}^{(m)}\right)$ and
$\left(A,B,C,D\right)=\left(\Phi,\,\rho\Phi^{\otimes
m},\,\nabla\left(\phi_{1,1},\,\phi_{2,1}\right),\,I_{2^{m}}\right)$ in (12),
we have

$\displaystyle\left(\Phi\,I_{2}\right)\otimes\left(\rho\phi^{\otimes
m}\psi_{1}^{(m)}\right)$ $\displaystyle=$
$\displaystyle\left(\Phi\otimes\rho\phi^{\otimes
m}\right)\left(I_{2}\otimes\psi_{1}^{(m)}\right)$
$\displaystyle\left(\Phi\nabla(\phi_{1,1},\phi_{2,1})\right)\otimes\left(\rho\phi^{\otimes
m}I_{2^{m}}\right)$ $\displaystyle=$
$\displaystyle\left(\varPhi\otimes\rho\phi^{\otimes
m}\right)\left(\nabla(\phi_{1,1},\phi_{2,1})\otimes I_{2^{m}}\right).$

Therefore

$\displaystyle\partial_{x}\left(\rho\phi_{1,1}\varPhi^{\otimes\,m+1}\right)+\partial_{y}\left(\rho\phi_{2,1}\varPhi^{\otimes\,m+1}\right)$
$\displaystyle=$
$\displaystyle\rho_{m+1}\left(I_{2}\otimes\psi_{1}^{(m)}+\nabla(\phi_{1,1},\phi_{2,1})\otimes
I_{2^{m}}\right),$ $\displaystyle=$ $\displaystyle\rho_{m+1}\psi_{1}^{(m+1)}.$

In a similar way, one proves that

$\displaystyle\partial_{x}\left(\rho\phi_{1,2}\Phi^{\otimes\,m+1}\right)+\partial_{y}\left(\rho\phi_{2,2}\Phi^{\otimes\,m+1}\right)=\rho_{m+1}\psi_{2}^{(m+1)}.$

So,
$div(\rho_{m+1}\otimes\varPhi)=\rho_{m+1}\left(\psi_{1}^{(m+1)},\,\psi_{2}^{(m+1)}\right).$
Let us prove that for a fixed $m\geq 1$ the family of polynomials
$\left\\{Q_{n,m}\right\\}_{n}$ is orthogonal with respect to
$\rho_{m}=\rho\Phi^{\otimes m}$. Taking $M=\mathbb{Q}_{n+1,m-1}^{t}$,
$N=\mathbb{Q}_{j+1,m-1}$ into (21), replacing $A^{\otimes\,m}$ by
$\rho\Phi^{\otimes\,m}$ and taking into account the definition of $\rho_{m}$
we have

$\displaystyle
div\left[\left(I_{2}\otimes\,\mathbb{Q}_{n+1,m-1}^{t}\right)\rho\Phi^{\otimes\,m}\nabla\,\mathbb{Q}_{j+1,m-1}\right]$
$\displaystyle=\mathbb{Q}_{n+1,m-1}^{t}div(\rho_{m-1}\otimes\Phi\nabla\mathbb{Q}_{j+1,m-1})+\left(\partial_{x}\mathbb{Q}_{n+1,m-1}^{t},\partial_{x}\mathbb{Q}_{n+1,m-1}^{t}\right)\rho_{m}\nabla\mathbb{Q}_{j+1,m-1}$

Integrating both sides on $\Omega$ we obtain

$\displaystyle\int_{\Omega}\left(\nabla\mathbb{Q}_{n+1,m-1}^{t}\right)^{t}\rho_{m}\nabla\mathbb{Q}_{j+1,m-1}dxdy$
$\displaystyle=$
$\displaystyle\int_{\Omega}div\left[\left(I_{2}\otimes\,\mathbb{Q}_{n+1,m-1}^{t}\right)\rho\Phi^{\otimes\,m}\nabla\,\mathbb{Q}_{j,m-1}\right]dxdy$
$\displaystyle-\int_{\Omega}\mathbb{Q}_{n+1,m-1}div(\left(\rho_{m-1}\otimes\Phi\right)\nabla\mathbb{Q}_{j+1,m-1})dxdy.$

Taking into account the boundary condition (26) as well as the second item of
Proposition 3, we obtain

$\displaystyle\int_{\Omega}\left(\nabla\mathbb{Q}_{n+1,m-1}^{t}\right)^{t}\rho_{m}\mathbb{Q}_{j+1,m-1}dxdy=-\int_{\Omega}\mathbb{Q}_{n+1,m-1}^{t}div(\left(\rho_{m-1}\otimes\Phi\right)\nabla\mathbb{Q}_{j,m-1})dxdy.$

Using (22) with $m$ taking for $m-1$, $B=\Phi,\,A^{\otimes\,m-1}$ taking for
$\rho\Phi^{\otimes\,m-1}$ and $N=\mathbb{Q}_{j+1,m-1}$, and taking into
account the fact that
$div\left(\rho_{m-1}\otimes\Phi\right)=\rho_{m-1}\left(\psi_{1}^{(m-1)},\,\psi_{1}^{(m-1)}\right)$
we obtain

$\displaystyle\int_{\Omega}\left(\nabla\mathbb{Q}_{n+1,m-1}^{t}\right)^{t}\rho_{m}\nabla\mathbb{Q}_{j+1,m-1}dxdy$
$\displaystyle=\int_{\Omega}\mathbb{Q}_{n+1,m-1}^{t}\rho_{m-1}\left[\left[\left(\left(\Phi\nabla\right)\cdot\nabla\right)\right]\mathbb{Q}_{j+1,m-1}+\left(\psi_{1}^{(m-1)},\,\psi_{2}^{(m-1)}\right)\nabla\mathbb{Q}_{j+1,m-1}\,\right]dxdy$
$\displaystyle=H_{n+1}\delta_{n,j},$

where $H_{n+1}$ is an $(n+1,n+1)$-invertible matrix because the family
$\left\\{\mathbb{Q}_{n+1,m-1}\right\\}$ is orthogonal with respect to
$\rho_{m-1}$ and the matrix polynomial
$\left[\left[\left(\left(\Phi\nabla\right)\cdot\nabla\right)\right]\mathbb{Q}_{j+1,m-1}+\left(\psi_{1}^{(m-1)},\,\psi_{2}^{(m-1)}\right)\nabla\mathbb{Q}_{j+1,m-1}\,\right]$
is of total degree $j+1$, thanks to Lemma 11.

[Step 1.2] $(b)\Rightarrow\,(c)$
We assume $(b)$ and fix two non-negative integers $n$ and $m$. Multiplying the
expansion

$\left(\left(\Phi\nabla\right)\cdot\nabla\right)\mathbb{Q}_{n,m}+\left(\psi_{1}^{(m)},\,\psi_{2}^{(m)}\right)\nabla\mathbb{Q}_{n,m}=\sum_{j=0}^{n}\mathbb{Q}_{j,m}A_{j+m,m}$

from the left side by $\rho_{m}$, and using (28) and then (22) with
$A^{\otimes\,m}=\rho\Phi^{\otimes\,m}$, $B=\Phi$ and $N=\mathbb{Q}_{n,m}$, we
obtain

$div\left[\left(\rho_{m}\otimes\Phi\right)\nabla\mathbb{Q}_{n,m}\right]=\sum_{j=0}^{n}\rho_{m}\mathbb{Q}_{j,m}A_{j+m,m}.$

Multiplying from the left side by $\mathbb{Q}_{n,k}^{t}$, $0\leq\,k\leq\,n$
and taking into account (21) with $M=\mathbb{Q}_{n,k}^{t}$,
$A^{\otimes\,m+1}=\rho\Phi^{\otimes\,m+1}=\rho_{m+1}$ and
$N=\mathbb{Q}_{n,m}$, we obtain

$div\left[\left(I_{2}\otimes\mathbb{Q}_{k,m}^{t}\right)\rho_{m+1}\nabla\mathbb{Q}_{n,m}\right]-\left(\nabla\mathbb{Q}_{k,m}\right)^{t}\rho_{m+1}\nabla\mathbb{Q}_{n,m}=\sum_{j=0}^{n}\mathbb{Q}_{k,m}^{t}\rho_{m}\mathbb{Q}_{j,m}A_{j+m,m}.$

Integrating both sides on $\Omega$ and taking into account the boundary
condition (26) as well as Proposition3, we obtain

$\int_{\Omega}\mathbb{Q}_{k,m}^{t}\rho_{m}\mathbb{Q}_{k,m}dxdyA_{k+m,m}=-\int_{\Omega}\left(\nabla\mathbb{Q}_{k,m}\right)^{t}\rho_{m+1}\nabla\mathbb{Q}_{n,m}dxdy.$

Since the family $\left\\{\nabla\mathbb{Q}_{n,m}\right\\}_{n\geq\,1}$ is
orthogonal with respect to $\rho_{m+1}$, $A_{k+m,m}=0$, $k<n$. Therefore,

$\left[\left(\Phi\nabla\right)\cdot\nabla\right]\mathbb{Q}_{n,m}+\left(\psi_{1}^{(m)},\,\psi_{2}^{(m)}\right)\nabla\mathbb{Q}_{n,m}=\mathbb{Q}_{n,m}A_{n+m,m}.$

Write
$\mathbb{Q}_{n,m}=\left(I_{2^{m}}\otimes\,X_{n}^{t}\right)G_{n,m,n}$+lower
terms and identify both sides using the relation (23) to obtain
$G_{n,m,n}\Lambda_{n+m,m}=-T_{n}^{(m)}G_{n,m,n}$.

[Step 1.3] $(c)\Rightarrow\,(a)$
Taking $m=0$ into (30) and multiplying by $\rho$ we have

$\left[\left(\rho\Phi\nabla\right)\cdot\nabla\right]\mathbb{P}_{n}^{t}+\rho\left(\psi_{1},\psi_{2}\right)\nabla\mathbb{P}_{n}^{t}+\rho\mathbb{P}_{n}^{t}\Lambda_{n,0}=0$

Taking into account the relation (22) with $m=0$, $B=\rho\Phi$ and
$N=\mathbb{P}_{n}^{t}$ we obtain

$div\left(\rho\Phi\nabla\mathbb{P}_{n}^{t}\right)+\left[\rho\left(\psi_{1},\,\psi_{2}\right)-div\left(\rho\Phi\right)\right]\nabla\mathbb{P}_{n}^{t}+\rho\mathbb{P}_{n}^{t}\Lambda_{n,0}=0.$

Integrating both sides on the domain $\Omega$ and using the boundary condition
(26), Proposition 3 as well as the orthogonality of the family
$\left\\{\mathbb{P}_{n}^{t}\right\\}_{n}$ with respect to $\rho$ ,we obtain

$\int_{\Omega}\left[\rho\left(\psi_{1},\,\psi_{2}\right)-div\left(\rho\Phi\right)\right]\nabla\mathbb{P}_{n}^{t}dxdy=0.$

Since $\\{\nabla\mathbb{P}_{n}^{t}\\}_{n\geq\,1}$ is a basis of
$\mathcal{P}\times\,\mathcal{P}$, we get

$\int_{\Omega}\left[\rho\left(\psi_{1},\,\psi_{2}\right)-div\left(\rho\Phi\right)\right]\mathbb{P}dxdy=0,\;{\rm
forall},\;\mathbb{P}=(p,q)^{t}\in\mathcal{P}\times\mathcal{P}.$

Taking $\mathbb{P}=(p,0)$ and using the fact that
$div(\rho\Phi)=\left(\partial_{x}\left(\rho\phi_{11}\right)+\partial_{y}\left(\rho\phi_{21}\right),\partial_{x}\left(\rho\phi_{12}\right)+\partial_{y}\left(\rho\phi_{22}\right)\right)$
we obtain

$\int_{\Omega}{\left(\rho\psi_{1}-\partial_{x}\left(\rho\phi_{11}\right)-\partial_{y}\left(\rho\phi_{21}\right)\right)\rho^{-1}}p\rho
dxdy=0,p\in\mathcal{P}.$

Without lost of generality, we assume
$\left(\partial_{x}\left(\rho\phi_{11}\right)+\partial_{y}\left(\rho\phi_{21}\right)\right)\rho^{-1}\in\,L^{2}(\Omega,\rho)$.
Therefore

$\left(\partial_{x}\left(\rho\phi_{11}\right)+\partial_{y}\left(\rho\phi_{21}\right)-\rho\psi_{1}\right)\rho^{-1}\in
L^{2}(\Omega,\rho).$

Thus,
$\partial_{x}\left(\rho\phi_{11}\right)+\partial_{y}\left(\rho\phi_{21}\right)=\rho\psi_{1}$,
by Theorem 7. In a similar way, taking $\mathbb{P}=(0,q),\,q\in\mathcal{P}$,
we obtain
$\partial_{x}\left(\rho\phi_{12}\right)+\partial_{y}\left(\rho\phi_{22}\right)=\rho\psi_{2}$.
Thus $div\left(\rho\phi\right)=\rho\left(\psi_{1},\,\psi_{2}\right)$.

  2. Step 2

$(b)\text{ and }(c)\Rightarrow\,(d)\Rightarrow\,(a)$ which, taking into
account Step 1, is equivalent to $(d)\Leftrightarrow(a)$
[Step 2.1] $(b)$ and $(c)\Rightarrow\,(d)$ From (30),

$\left[\left(\Phi\nabla\right)\cdot\nabla\right]\nabla^{(m)}\mathbb{P}_{n+m}^{t}+\left(\psi_{1}^{(m)},\psi_{2}^{(m)}\right)\nabla\nabla^{(m)}\mathbb{P}_{n+m}^{t}+\nabla^{(m)}\mathbb{P}_{n+m}^{t}\Lambda_{n+m,m}=0.$

Replace $n$ by $n-m$, multiply the obtained equation by $\rho_{m}$ from the
left side and use (22), with $A^{\otimes m}=\rho_{m}$, $B=\Phi$ and
$N=\mathbb{P}_{n}^{t}$, to obtain

$div\left(\left(\rho_{m}\otimes\Phi\right)\nabla^{(m+1)}\mathbb{P}_{n}^{t}\right)+\rho_{m}\nabla^{(m)}\mathbb{P}_{n}^{t}\Lambda_{n,m}=0,$

Replace $m$ by $m-1$ and iterate the obtained relation to obtain

$div^{(m)}\left(\rho_{m}\nabla^{(m)}\mathbb{P}_{n}^{t}\right)=(-1)^{n}\rho\mathbb{P}_{n}^{t}\Lambda_{n,m-1}\Lambda_{n,m-2}\dots\,\Lambda_{n,0}.$

Take $m=n$ to get the result.
[Step 2.2] $(d)\Rightarrow\,(a)$ Take $n=1$ into (d) to obtain

$div(\rho\Phi)\nabla\,\mathbb{P}_{1}^{t}=-\rho\mathbb{P}_{1}^{t}\Lambda_{1,0}$

Since $\mathbb{P}_{1}$ is monic,
$P_{1}^{t}(x,y)=X^{t}+G_{1,0,0}=(x+g_{10},\,y+g_{01})$ and
$\nabla\,P_{1}^{t}=I_{2}$. From (30),

$G_{1,0,1}\Lambda_{1,0}=-\left(L_{0}^{(0)}\right)^{t}\left(D_{1},\,D_{2}\right)N_{1}^{(0)}G_{1,0,1}.$

That is $\Lambda_{1,0}=-(D_{1},D_{2})$ for
$G_{1,0,1}=L_{0}^{(0)}=N_{1}^{(0)}=I_{2}$ see Proposition 5 and Eq. (17).
Therefore writing $D_{i}=(d_{i1},\,d_{i2}),\,i=1,2$, we obtain
$div(\rho\Phi)=\rho\mathbb{P}_{1}^{t}(D_{1},D_{2})=\rho\left(\psi_{1},\psi_{2}\right)$
with $\psi_{1}(x,y)=xd_{11}+yd_{21}+g_{10}d_{11}+g_{01}d_{21}$ and
$\psi_{2}(x,y)=xd_{12}+yd_{22}+g_{10}d_{12}+g_{01}d_{22}$ which are
polynomials of total degree one, for $det(D_{1},D_{2})\neq\,0$.

  3. Step 3

$(b)\Rightarrow\,(e)\Rightarrow\,(a)$ which taking into account Step 1 is
equivalent to $(e)\Leftrightarrow\,(a)$.
[Step 3.1] $(b)\Rightarrow\,(e)$ Use the relation (11) to rewrite the right
side of the expansion

$\left(\Phi\otimes\,I_{2^{m}}\right)\nabla\,\mathbb{Q}_{n,m}=\sum_{k=0}^{n+1}\left(I_{2}\otimes\mathbb{Q}_{k,m}\right)A_{k}^{n,m}$

and multiply, for a fixed $0\leq\,j\leq n$, both sides, from the left, by

$\left(\begin{matrix}\mathbb{Q}_{j,m}^{t}&0\\\
0&\mathbb{Q}_{j,m}^{t}\end{matrix}\right)\left(\begin{matrix}\rho_{m}&0\\\
0&\rho_{m}\end{matrix}\right).$

Observe that $\rho_{m+1}=\rho_{m}\left(\Phi\otimes\,I_{2^{m}}\right)$ to have

$\left(\begin{matrix}\mathbb{Q}_{j,m}^{t}&0\\\
0&\mathbb{Q}_{j,m}^{t}\end{matrix}\right)\rho_{m+1}\mathbb{Q}_{n-1,m+1}=\sum_{k=0}^{n+1}\left(\begin{matrix}\mathbb{Q}_{j,m}^{t}\rho_{m}\mathbb{Q}_{k,m}&0\\\
0&\mathbb{Q}_{j,m}^{t}\rho_{m}\mathbb{Q}_{k,m}\end{matrix}\right)A_{k}^{n,m}.$

Integrate both sides on the domain $\Omega$ and use the fact
$\left\\{\mathbb{Q}_{n,m}\right\\}_{n}$ (resp.
$\left\\{\mathbb{Q}_{n-1,m+1}\right\\}_{n}$) is orthogonal with respect to
$\rho_{m}$ (resp. $\rho_{m+1}$) to obtain

$\displaystyle\left(\begin{matrix}H_{j,m}&0\\\
0&H_{j,m}\end{matrix}\right)A_{j}^{n,m}$ $\displaystyle=$
$\displaystyle\int_{\Omega}\left(\begin{matrix}\mathbb{Q}_{j,m}^{t}&0\\\
0&\mathbb{Q}_{j,m}^{t}\end{matrix}\right)\rho_{m+1}\mathbb{Q}_{n-1,m+1}dxdy$
$\displaystyle=$ $\displaystyle\begin{cases}&0\;{\rm if}\;0\leq\,j<n-1\\\
&M_{n-1,m}\;{\rm if}\;j=n-1,\end{cases}$

where $H_{j,m}=\int_{\Omega}\mathbb{Q}^{t}_{j,m}\rho_{m}\mathbb{Q}_{j,m}dxdy$
and $M_{n-1,m}$ are invertible matrices. Therefore, $A_{j}^{n,m}=0$, for
$j=0,..,n-2$ and $A_{n-1}^{n,m}$ is an invertible matrix.
[Step 3.2] $(e)\Rightarrow\,(a)$ Let

$\frac{div(\rho\Phi)}{\rho}=\sum_{n=0}^{\infty}\mathbb{P}^{t}_{n}A_{n}$

be the formal Fourier expansion of the function $\frac{div(\rho\Phi)}{\rho}$
in the system $\\{\mathbb{P}_{n}\\}_{n\in\mathbb{N}}$. Since the system is
orthogonal with respect to $\rho$, the coefficients $A_{n}$ are given by

$\int_{\Omega}\mathbb{P}_{n}\mathbb{P}_{n}^{t}dxdy\,A_{n}=\int_{\Omega}\mathbb{P}_{n}div(\rho\Phi)dxdy.$

Take into account the relation

$div\left(\rho\Phi\otimes\mathbb{P}_{n}\right)=\mathbb{P}_{n}div\left(\rho\Phi\right)+\left(\partial_{x}\mathbb{P}_{n},\partial_{y}\mathbb{P}_{n}\right)\rho\Phi,$

obtained by direct computation, as well as the boundary condition (26) to
obtain

$\int_{\Omega}\mathbb{P}_{n}\mathbb{P}_{n}^{t}\rho\,dxdy\,A_{n}=-\int_{\Omega}\left(\partial_{x}\mathbb{P}_{n},\partial_{y}\mathbb{P}_{n}\right)\rho\Phi\,dxdy.$

Transpose and take into account the fact that the matrix $\Phi$ is symmetric
to obtain

$A_{n}^{t}\int_{\Omega}\rho\mathbb{P}_{n}\mathbb{P}_{n}^{t}dxdy=-\int_{\Omega}\rho\Phi\nabla\mathbb{P}_{n}^{t}\,dxdy.$

Therefore $A_{0}=0$. Use the assumption to get $A_{n}=0$, $n>1$ and

$\displaystyle
A_{1}^{t}\int_{\Omega}\rho\mathbb{P}_{1}\mathbb{P}_{1}^{t}dxdy=$
$\displaystyle-\int_{\Omega}\left(\begin{matrix}\rho\mathbb{P}_{2}^{t}&0\\\
0&\rho\mathbb{P}_{2}^{t}\end{matrix}\right)A_{2}^{1,0}dxdy-\int_{\Omega}\left(\begin{matrix}\rho\mathbb{P}_{1}^{t}&0\\\
0&\rho\mathbb{P}_{1}^{t}\end{matrix}\right)A_{1}^{1,0}dxdy$
$\displaystyle-\int_{\Omega}\left(\begin{matrix}\rho\mathbb{P}_{0}^{t}&0\\\
0&\rho\mathbb{P}_{0}^{t}\end{matrix}\right)A_{0}^{1,0}dxdy.$

Use the fact that the system $\\{\mathbb{P}_{n}\\}_{n}$ is orthogonal with
respect to $\rho$ to obtain

$A_{1}^{t}\int_{\Omega}\rho\mathbb{P}_{1}\mathbb{P}_{1}^{t}dxdy=-\left(\begin{matrix}a_{0}&0\\\
0&a_{0}\end{matrix}\right)A_{0}^{n,0},\;\;a_{0}=\int_{\Omega}\rho\mathbb{P}_{0}^{t}dxdy.$

Therefore $A_{1}=\left(\begin{matrix}b_{0}&0\\\ 0&b_{1}\end{matrix}\right)$,
with $b_{0}b_{1}\neq\,0$ and $b_{1}\neq 0$ for
$\int_{\Omega}\rho\mathbb{P}_{1}\mathbb{P}_{1}^{t}dxdy$ and $A_{0}^{1,0}$ are
invertible. Thus,

$\displaystyle{div(\rho\Phi)}$ $\displaystyle=$
$\displaystyle-\mathbb{P}_{1}^{t}\left(\begin{matrix}b_{0}&0\\\
0&b_{1}\end{matrix}\right)$ $\displaystyle=$
$\displaystyle\left(\psi_{1},\,\psi_{2}\right),$

where $\psi_{1}$ and $\psi_{2}$ are polynomials of degree 1 of the variable
$x$ and $y$ respectively.

∎

We now propose the following definition of classical orthogonal polynomials in
two variables.

######  Definition 12.

Let $\Omega$ be a domain of $\mathbb{R}^{2}$. A weight function $\rho$ on
$\Omega$ is classical if there exists a symmetric $2$-matrix vector
polynomials $\Phi$ of total degree at most $2$ satisfying the Neumann boundary
condition (26) and the differential system (29a)-(29b), and there exist two
polynomials $\psi_{i}(x,y)=X^{t}D_{i}+E_{i},\,i=1,2$ of polynomials of total
degree $1$ such that $det(D_{1},D_{2})\neq\,0$ and

$div\left(\rho\Phi\right)=\rho\left(\psi_{1},\,\psi_{2}\right).$

##  4 Concluding remarks

In this paper we characterize polynomials in two variables, orthogonal with
respect to a weight function satisfying (3), by means of five equivalent
properties and derive a definition of classical orthogonal polynomials (COP)
in two variables. The properties appearing in the main theorem of this work
are an analogue of those in [1] for COP in one variable. Our definition of
classical orthogonal in two variables includes the one in [5, Definition 3].
As examples of orthogonal polynomials in two variables that are classical in
our sense, we have the product of Hermite polynomials, the product of Laguerre
polynomials, the product of Hermite and Laguerre polynomials, orthogonal
polynomials on the triangle (see [6, p. 38]) and the product of Jacobi
polynomials. This last example is the Class V Koornwinder polynomials [8] that
are not classical in the Krall and Sheffer sense [11]. A natural question that
arises is whether a classification of all the known families of orthogonal
polynomials in two variables, that satisfy our definition, is possible. Work
on an analysis along these lines is in progress.

##  5 Acknowledgment

The first author would like to express his appreciation to the organizers of
OPSFA 2019 during which the idea of this paper was conceived. The authors are
grateful to the University of South Africa for an International General
Research grant which supported a research visit by the first author to the
University of South Africa for the period 1 Oct 2021–11 Dec 2021.

## References

  * [1] Al-Salam WA. Characterization theorems for orthogonal polynomials. In: Nevai P, editor. Orthogonal polynomials. Dordrecht: Springer; 1990. p. 1–24.
  * [2] Álvarez de Morales M, Fernández L, Pérez TE, Piñar MA. A matrix Rodrigues formula for classical orthogonal polynomials in two variables. J. Approx. Theory. 2009;157(1):32–52.
  * [3] Bakry D, Orevkov S and Zani M. Orthogonal polynomials and diffusion operators. arXiv preprint arXiv:1309.5632. 2013 Sep 22.
  * [4] Bellman R. Introduction to matrix analysis. Society for Industrial and Applied Mathematics; 1997. (Classics in Applied Mathematics).
  * [5] Fernández L, Pérez TE, Piñar MA. Classical orthogonal polynomials in two variables: a matrix approach. Numer. Algorithms. 2005;39(1):131–142.
  * [6] Dunkl CF, Xu Y. Orthogonal polynomials of several variables. Cambridge University Press; 2014.
  * [7] Fornberg B. A practical guide to pseudospectral methods. Cambridge university press; 1998.
  * [8] Koornwinder T. Two-variable analogues of the classical orthogonal polynomials. In: Askey RA, editor. Theory and application of special functions. Academic Press; 1975 p. 435-495.
  * [9] Kowalski MA. Orthogonality and recursion formulas for polynomials in $n$ variables. SIAM J. Math. Anal.. 1982;13(2):316–323.
  * [10] Kowalski MA. The recursion formulas for orthogonal polynomials in $n$ variables. SIAM J. Math. Anal.. 1982;13(2):309–315.
  * [11] Krall HL, Sheffer IM, Orthogonal polynomials in two variables. Ann. Mat. Pura Appl.. 1967;76(1):325–376.
  * [12] Littlejohn LL. Orthogonal polynomial solutions to ordinary and partial differential equations. In: Alfaro M, Dehesa JS, Marcellán FJ, Rubio de Frances JL, Vinuesa J, editors. Orthogonal Polynomials and their Applications. Proceedings of an International Symposium held in Segovia, Spain, Sept.22–27, 1986. Berlin: Springer; 1988. p.98–124. (Lecture Notes in Mathematics, volume 1329).
  * [13] Lyskova AS. Orthogonal polynomials of several variables. In: Doklady Akademii Nauk, vol. 316, no. 6, pp. 1301-1306. Russian Academy of Sciences, 1991.
  * [14] Marcellán F, Marriaga M, Pérez TE and Piñar MA. On bivariate classical orthogonal polynomials. Appl. Math. Comput.. 2018;325:340–357.
  * [15] Kenfack Nangho M. and Jiejip Nkwamouo B. Embedding theorems and completeness of orthogonal polynomials systems associated with classical Weights in two variables, In progress.
  * [16] Nikiforov AF, Suslov SK and Uvarov VB. Classical orthogonal polynomials of a discrete variable on nonuniform grids. Dokl. Akad. Nauk SSSR, 1986;291(5):1056–1059.
  * [17] Williams P. Jacobi pseudospectral method for solving optimal control problems. J. Guid. Control Dynam. 2004;27(2): 293–297.
  * [18] Suetin PK. Orthogonal Polynomials in two variables. Analytical Method and Special functions. Amsterdam: Gordon and Breach Science Publishers; 1984.
  * [19] Xu Y. On multivariate orthogonal polynomials. SIAM J. Math. Anal.. 1993;24:783–794.