Sobolev orthogonality of polynomial solutions of second-order partial differential equations

Given a second-order partial differential operator L\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathscr {L}}$$\end{document} with nonzero polynomial coefficients of degree at most 2, and a Sobolev bilinear form (P,Q)S=∑i=0N∑j=0iu(i,j),∂xi-j∂yjP∂xi-j∂yjQ,N⩾0,\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} (P,Q)_S\,=\,\sum _{i=0}^N\sum _{j=0}^i\left\langle {\textbf{u}}^{(i,j)}, \partial _x^{i-j}\partial _y^jP\,\,\partial _x^{i-j}\partial _y^{j}Q\right\rangle , \quad N\geqslant 0, \end{aligned}$$\end{document}where u(i,j)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\textbf{u}}^{(i,j)}$$\end{document}, 0⩽j⩽i⩽N\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0\leqslant j \leqslant i \leqslant N$$\end{document}, are linear functionals defined on the space of bivariate polynomials, we study the orthogonality of the polynomial solutions of the partial differential equation L[p]=λn,mp\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathscr {L}}[p]=\lambda _{n,m}\,p$$\end{document} with respect to (·,·)S\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\cdot ,\cdot )_S$$\end{document}, where λn,m\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _{n,m}$$\end{document} are eigenvalue parameters depending on the total and partial degree of the solutions. We show that the linear functionals in the bilinear form must satisfy Pearson equations related to the coefficients of L\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathscr {L}}$$\end{document}. Therefore, we also study solutions of the Pearson equations that can be obtained from univariate moment functionals. In fact, the involved univariate functionals must satisfy Pearson equations in one variable. Moreover, we study polynomial solutions of L[p]=λn,mp\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathscr {L}}[p]=\lambda _{n,m}\,p$$\end{document} obtained from univariate sequences of polynomials satisfying second-order ordinary differential equations.