Higher-order differential operators having bivariate orthogonal polynomials as eigenfunctions

We introduce a systematic method for constructing higher-order partial differential equations for which bivariate orthogonal polynomials are eigenfunctions. Using the framework of moment functionals, the approach is independent of the orthogonality domain’s geometry, enabling broad applicability across different polynomial families. Applications to classical weight functions on the unit disk and triangle modified by measures defined on lower-dimensional manifolds are presented.