\documentclass[reqno]{amsart} \usepackage{hyperref} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2016 (2016), No. 22, pp. 1--21.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu} \thanks{\copyright 2016 Texas State University.} \vspace{8mm}} \begin{document} \title[\hfilneg EJDE-2016/22\hfil Higher order self-adjoint operators] {Higher order self-adjoint operators with polynomial coefficients} \author[H. Azad, A. Laradji, M. T. Mustafa \hfil EJDE-2016/22\hfilneg] {Hassan Azad, Abdallah Laradji, Muhammad Tahir Mustafa} \address{Hassan Azad \newline Department of Mathematics and Statistics, King Fahd University of Petroleum and Minerals, Dhahran 31261, Saudi Arabia} \email{hassanaz@kfupm.edu.sa} \address{Abdallah Laradji \newline Department of Mathematics and Statistics, King Fahd University of Petroleum and Minerals, Dhahran 31261, Saudi Arabia} \email{alaradji@kfupm.edu.sa} \address{Muhammad Tahir Mustafa \newline Department of Mathematics, Statistics and Physics, Qatar University, Doha 2713, Qatar} \email{tahir.mustafa@qu.edu.qa} \thanks{Submitted October 21, 2016. Published January 18, 2017.} \subjclass[2010]{33C45, 34A05, 34A30, 34B24, 42C05} \keywords{Self-adjoint operators; polynomial coefficients} \begin{abstract} We study algebraic and analytic aspects of self-adjoint operators of order four or higher with polynomial coefficients. As a consequence, a systematic way of constructing such operators is given. The procedure is applied to obtain many examples up to order 8; similar examples can be constructed for all even order operators. In particular, a complete classification of all order 4 operators is given. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{lemma}[theorem]{Lemma} \newtheorem{proposition}[theorem]{Proposition} \newtheorem{corollary}[theorem]{Corollary} \newtheorem{example}[theorem]{Example} \allowdisplaybreaks \section{Introduction} The classification of self-adjoint second order operators with polynomial coefficients is a classical subject going back to Brenke \cite{brenke}. This paper is a contribution to certain algebraic and analytic aspects of higher order self-adjoint operators with polynomial coefficients. Its main aim is to construct such operators. This involves determining the explicit differential equations for the polynomial coefficients of the operators and the boundary conditions which ensure self-adjointness. These operators are not in general iterates of second order classical operators - as stated in \cite{littlejohnpams1994}; (cf. \cite{littlejohn-krall89}). As the weight function which makes these operators self-adjoint depends only on the first two leading terms of the operator, therefore, if one can find a second operator with the same weight function, the eigenpolynomials for both operators would be the same; see Section 4. We should point out that some of the most important recent contributions to this subject are due to Kwon, Littlejohn and Yoon \cite{kwon2001}, Bavnick \cite{bavinck95, bavinck98, bavinck2000}, Koekoek \cite{koekoek93} and Koekoek-Koekoek \cite{kockock2000}; see also the references therein. A classical reference for higher order Sturm-Liouville theory is the book of Ince \cite[Chap.IX]{ince}. This theory was revived by Everitt in \cite{everitt57}; see also Everitt et al.\ \cite{everitt2001}. Classical references that deal with various aspects of polynomial solutions of differential equations are the references \cite{ Boc, brenke,lesky62, Kra40, Rou}. More recent papers that deal with the same subject are \cite{azad1, azad2, koornwinder84, krall81, littlejohn84, Turbiner}. A reference for the related topic of orthogonal polynomials is \cite{Sze}. A more recent reference for this topic which also has an extensive bibliography is the book \cite{Ism}. We should point out that the classification problem we address is related to the Bochner-Krall problem of finding all families of orthogonal polynomials that are eigenfunctions of a differential operator (see \cite{Boc, brenke, everitt2001, krall81, kwon2006}). Here is a more detailed description of the results of this paper. We consider linear differential operators with polynomial coefficients that map the space of polynomials of degree at most $k$ into itself - for all $k$. Proposition \ref{prop2.1} gives necessary and sufficient conditions for such an $n^{\mathrm{th}}$ order linear operator to be self-adjoint. The integrability and asymptotic properties of the weight function and its derivatives near the zeroes of the leading term given in Propositions \ref{prop3.2} and \ref{prop3.3} carry enough information to determine the form of the first two terms of the operator in specific cases. This is used further to determine the full operator, using the boundary conditions and the determining equations of Proposition \ref{prop2.1}. Although such determining equations are known (see for example \cite{littlejohn-race1990}), boundary conditions involving all polynomial coefficients of the linear operator do not seem to have been considered earlier in the general case and these are equally crucial to the construction of operators given here. We should point out that the $( n+1) $ determining equations in Proposition \ref{prop2.1} are equivalent to a system of $\frac{n}{2}$ equations, as shown in \cite{littlejohn-race1990}. The classical second order operators are completely determined by the integrability of the associated weight function. In general, if the operator is of the form $L(y)= a y^{(n)}+b y^{(n-1)}+ \dots$ then, as shown in Proposition \ref{prop3.4}, for $n>2$, the multiplicity of each root of $a$ is at least $2$ and its multiplicity in $b$ is less than its multiplicity in $a$ and it is of multiplicity at least $1$ in $b$. In particular, if the operator is of fourth order and the leading term has distinct roots, then every root occurs with multiplicity 2 and therefore the leading term must have exactly two distinct roots and their multiplicity in the next term is $1$ - as shown in Proposition \ref{prop3.4}. In \cite{littlejohn-krall89}, fourth order Sturm-Liouville systems were given for ordinary weights and for weights involving distributions. In particular, for ordinary weights, the authors found operators that are iterates of second order ones. In this paper, a systematic way of producing Sturm-Liouville systems with ordinary weights for all even orders is given. Although it is believed that all fourth order operators with classical weights are iterates of second order ones (see \cite{littlejohnpams1994, littlejohn-krall89}), the results we obtain show this is in fact incorrect. Indeed, the classification of fourth order operators we provide here (see Section 4.2) includes examples of operators that cannot be iterates of second order ones, as an eigenvalue argument shows. The examples II.6, II.7 of fourth order operators in \cite{littlejohn-krall89} with weights involving delta and Heaviside functions and the sixth order operator in Example 4.1 in \cite{kwon2001} with weight involving delta functions are solutions of the determining equations given in Section 4. For these 4th order operators of \cite{littlejohn-krall89}, the 3rd boundary condition in the sense of \eqref{boundary-n-4} fails at one or both the boundary points. For the sixth order operator of Example 4.1 in \cite{kwon2001} the last three boundary conditions in \eqref{boundary-n-6} do not hold. Many such examples can be constructed but the solutions of the determining equations are too many to be listed efficiently (cf. Example 4.1, Section 4.3). We also obtain examples of sixth and eighth order operators where all the boundary conditions hold. Due to space constraint, these examples are provided in the expanded online version of the paper at http://arxiv.org/abs/1409.2523. In fact similar examples can be constructed for any even order; such constructions involve increased computational complexity. All the solutions presented in the examples have been first verified using Mathematica and then directly generated (from the same file) by Mathematica as LaTeX output for the paper. \section{Algebraic aspects of higher order Sturm-Liouville theory} Consider, on the space $C^{\infty}$, the $n$th-order linear operator $L=\sum_{k=1}^na_{k}(x)D^{k}$, where $D$ is the usual differential operator and each $a_{k}:=a_{k}(x)$ is a polynomial of degree at most $k$. In this way, for each natural number $N$, the vector space $\mathbb{P}_{N}$ of all polynomials of degree at most $N$ is $L$-invariant. Our first objective is to obtain conditions on the polynomials $a_{k}$ for the existence of an inner product $\langle u,v\rangle =\int_Ipuvdx$ on $C^{\infty}$ for which $L$ is self-adjoint, and where the weight $p$ is sufficiently differentiable in a real interval $I$ where $a_n$ does not vanish. This smoothness assumption is reasonable since it is satisfied by all weights of classical orthogonal polynomials. For a function $f$ and an interval $J$ with (possibly infinite) endpoints $\alpha<\beta$, $\partial J$ denotes the boundary $\{\alpha,\beta\}$ of $J$, and $[ f]_{J}$ means $\lim_{x\to\beta^{-}}f(x)-\lim_{x\to\alpha^{+}}f(x)$, where both limits are finite. The statement that a function vanishes at an endpoint of the interval is to be understood in the sense of limits. For notational convenience, let $b_{j}:=pa_{j}$ ($1\leq j\leq n)$, and $b_0$ be the zero function. Our main result in this section is the following proposition. As mentioned in the introduction, the determining equations in (i) are known (see for example \cite{littlejohn-race1990}). \begin{proposition} \label{prop2.1} With the above notation, for $L$ to be self-adjoint with respect to the inner product $\langle u,v\rangle =\int_I p(x)u(x)v(x)dx$, it is necessary that \[ p(x)=\frac{\exp\big(\frac{2}{n}\int\frac{a_{n-1}(x)}{a_n(x)}dx\big)}{| a_n(x)| } \] on $I$ and that $n$ be even. Conversely if $p(x)=\frac{\exp\big(\frac{2}{n}\int\frac{a_{n-1}(x)}{a_n(x)}dx\big)}{| a_n(x)| }$ and $b_{j}=pa_{j}$, it is necessary and sufficient for $L$ to be self adjoint that the following conditions hold. \begin{itemize} \item[(i)] $(-1)^{j}\binom{j}{j}b_{j}+\dots+( -1) ^n \binom{n}{j}b_n^{(n-j)}=b_{j}$ on $I$ for $0\leq j\leq n$; \item[(ii)] \[ \begin{bmatrix} \binom{n-1}{j-1} & -\binom{n-2}{j-1} & \dots & (-1)^{n-j}\\ \binom{n-2}{j-2} & -\binom{n-3}{j-2} & \dots & (-1)^{n-j}\\ \dots & \dots & \dots & \dots\\ \binom{n-j}{0} & -\binom{n-j-1}{0} & \dots & (-1)^{n-j} \end{bmatrix} \begin{bmatrix} b_n^{(n-j)}\\ b_{n-1}^{(n-j-1)}\\ \dots\\ b_{j} \end{bmatrix} =0 \] on $\partial I$ for $1\leq j\leq n$. In particular, $n$ must be even. \end{itemize} \end{proposition} To prove the above proposition we need the following lemmas. The first is a formula for repeated integration by parts, while the last one may be of independent interest. \begin{lemma}[\cite{naimark}] \label{lem2.2} Let $f$ and $y$ be functions $k$ times differentiable on some interval $I$. Then \[ \int_Ify^{(k)}dx=(-1)^{k}\int_If^{(k)}y\,dx+[ \sum_{j=0}^{k-1}(-1)^{j}f^{(j)}y^{(k-1-j)}]_I. \] \end{lemma} \begin{lemma} \label{lem2.3} Let $f_{j}$ $(1\leq j\leq r)$ be functions continuous on some interval $(a,b)$, where $b$ may be infinite. If there exists a non-singular square matrix $A=[a_{ij}]_{1\leq i,j\leq r}$ such that $\lim_{x\to b^{-}}\sum_{j=1}^{r}a_{ij}f_{j}(x)=0$ for $1\leq i\leq r$, then $\lim_{x\to b^{-}}f_{j}(x)=0$ for $1\leq j\leq r$. \end{lemma} \begin{proof} Put $g_{i}(x)=\sum_{j=1}^{r}a_{ij}f_{j}(x)$ $(1\leq i\leq r)$, so that for all $x\in(a,b)$, \[ A(f_1(x),f_2(x),\dots,f_{r}(x))^{T}=(g_1(x),g_2(x),\dots,g_{r}(x))^{T}. \] The conclusion follows from the fact that $\lim_{x\to b^{-}}g_{i}(x)=0$ and that \[ (f_1(x),f_2(x),\dots,f_{r}(x))^{T}=A^{-1}(g_1(x),g_2(x),\dots,g_{r} (x))^{T}. \] \end{proof} \begin{lemma} \label{lem2.4} Let $v_{i}$ $(0\leq i\leq r)$ be functions continuous on a real interval $I$ such that for all polynomials $u$, $[ \sum_{i=0}^{r} v_{i}u^{(i)}]_I=0$. Then $v_{i}=0$ $(0\leq i\leq r)$ at each endpoint of $I$. \end{lemma} \begin{proof} Suppose first that the endpoints $\alpha,\beta$ ($\alpha<\beta)$ of $I$ are finite. We need only show that $v_{r}=0$ at each endpoint of $I$, the statement for the remaining $v_{i}$ would then follow by straightforward (reverse) induction. From $u=(x-\alpha)^{r}(x-\beta)^{r}z$ $(z$ a polynomial$)$ we get $[ r!v_{r}z]_I=0$ i.e. $[ v_{r}z]_I=0$. Then, from $z=1$ and $z=x$ respectively, we get $v_{r}(\beta)-v_{r}(\alpha)=0$ and $\beta v_{r}(\beta)-\alpha v_{r} (\alpha)=0$. These equations imply $v_{r}(\beta)=v_{r}(\alpha)=0$, as required. Suppose now that $\beta=\infty$ (with $\alpha$ possibly infinite). Put $u_{j}=x^{r+j}$ $(1\leq j\leq r+1)$. Then, since $\lim_{x\to\infty}x^{j}=\infty$ and $\lim_{x\to \infty}\sum_{i=0}^{r}\frac{(r+j)!}{(r+j-i)!}x^{r+j-i}v_{i}(x)$ is finite for each $j$ ($1\leq j\leq r+1)$, we obtain $\sum_{i=0}^{r}\frac{(r+j)!x^{r-i}v_{i}(x)}{(r+j-i)!}\to0$ as $x\to\infty$. If we now put $f_{i}(x)=x^{r-i}v_{i}(x)$, $0\leq i\leq r$, we get $\lim_{x\to\infty}\sum_{i=0}^{r}a_{ij}f_{i}(x)=0$ for $1\leq j\leq r+1$, where $a_{ij}=\frac{(r+j)!}{(r+j-i)!}$. In view of Lemma \ref{lem2.3}, we need only show that the matrix $A=[a_{ij}]_{0\leq i\leq r,1\leq j\leq r+1}$ is non-singular. We have $a_{ij}=i!\binom{r+j}{i}$ hence $A$ is non-singular if and only if the matrix $B=[ \binom{r+j}{i}]_{0\leq i\leq r,1\leq j\leq r+1}$ is non-singular. If, more generally, we let $D(m,n):=\det[\binom{m+j-1}{i-1}]_{1\leq i,j\leq n+1}$ for $m\geq n\geq0$, then it is easy to show that $D(m,n)=D(m,n-1)$. This gives $D(m,n)=D(m,1)=1$, so that $\det B=D(r+1,r)\neq0$, i.e. $A$ is non-singular. The case when $\alpha$ is infinite is dealt with in a similar way. \end{proof} \begin{lemma} \label{lem2.5} Let $c$ and $c_{j}$ $(0\leq j\leq n)$ be functions continuous on an interval $I$. If $[ \sum_{j=0}^nc_{j}y^{(j)}]_I=\int_Icy\,dx$ for all infinitely differentiable functions $y$, then $c=0$ on $I$ and each $c_{j}$ vanishes at each endpoint of $I$. \end{lemma} \begin{proof} We first prove that $c=0$ on $I$. Suppose on the contrary that $c(\gamma)\neq0$ for some $\gamma$ in $I$. We can assume that $c>0$ on some subinterval $[\delta,\varepsilon] $ of $I$ containing $\gamma$. Let \[ \phi( x) =\begin{cases} \exp\big(\frac{1}{( x-\delta) ( x-\varepsilon)}) & \text{if }\delta\leq x\leq\varepsilon\\ 0 & \text{otherwise} \end{cases} \] Then, putting $y=\phi(x)$, we get $[ \sum_{j=0}^nc_{j}y^{(j)}]_I=0$ and so $\int_Icy\,dx=\int_{\delta}^{\varepsilon}c\phi(x)dx=0$. This is impossible since the integrand in this last integral is positive. This shows that $c=0$ on $I$. By Lemma \ref{lem2.4}, each $c_{j}$ equals zero at each endpoint of $I$, and the proof is complete. \end{proof} \begin{proof}[Proof of Proposition \ref{prop2.1}] By Lemma \ref{lem2.2}, for any functions $y$ and $u$ in $C^{\infty}$ we have \begin{equation} \label{e2.1} \begin{aligned} \langle Ly,u\rangle &=\int_Ip(Ly)udx=\sum_{k=1}^n \int_Ipa_{k}y^{(k)}udx =\sum_{k=1}^n\int_I( pa_{k}u) y^{(k)}dx \\ & =\sum_{k=1}^n\Big( \Big[ \sum_{j=0}^{k-1} (-1)^{j}(pa_{k}u)^{(j)}y^{(k-1-j)}\Big]_I+(-1)^{k}\int_I (pa_{k}u)^{(k)}y\,dx\Big). \end{aligned} \end{equation} Suppose first that $L$ is self-adjoint. Then $\langle Ly,u\rangle =\langle y,Lu\rangle $, i.e. \begin{equation} \label{e2.2} \begin{aligned} &\sum_{k=1}^n\Big( \Big[ \sum_{j=0}^{k-1}(-1)^{j} (pa_{k}u)^{(j)}y^{(k-1-j)}\Big]_I+(-1)^{k}\int_I(pa_{k}u)^{(k)} y\,dx\Big) \\ & =\sum_{k=1}^n\int_Ipa_{k}u^{(k)}y\,dx. \end{aligned} \end{equation} Fix $u$ and put \[ c=\sum_{k=1}^n(-1)^{k}(b_{k}u)^{(k)} -\sum_{k=1}^nb_{k}u^{(k)}, \quad c_{j}=\sum_{k=j+1}^n(-1)^{k-1-j}( b_{k}u) ^{(k-1-j)} \] $(0\leq j\leq n-1)$. Then \eqref{e2.2} gives $[\sum_{j=0}^{n-1}c_{j}y^{(j)}]_I=\int_Icy\,dx$, which by Lemma \ref{lem2.5} implies $c=0$ on $I$ and $c_{j}=0$ $(0\leq j\leq n-1)$ at each endpoint of $I$. Applying Leibniz rule to the terms $(b_{k}u) ^{(k)}$ of $c$ we obtain \begin{equation} \sum_{k=1}^nb_{k}u^{(k)}=\sum_{k=1}^n( -1) ^{k}\sum_{j=0}^{k}\binom{k}{j}b_{k}^{(k-j)}u^{(j)}. \label{e2.3} \end{equation} Since this is true for all $u$ in $C^{\infty}$, we can equate coefficients of $u^{(k)}$ and get from $k=n$ the classical fact that $n$ is even and that for $0\leq k\leq n-1$, \begin{equation} \begin{gathered} 2b_{k}=\binom{k+1}{k}b_{k+1}'-\binom{k+2}{k}b_{k+2}^{'' }+\dots+\binom{n}{k}b_n^{( n-k) }\quad \text{if $k$ is odd}\\ 0=-\binom{k+1}{k}b_{k+1}'+\binom{k+2}{k}b_{k+2}''-\dots +\binom{n}{k}b_n^{( n-k) }\quad \text{ if $k$ is even} \end{gathered} \label{e2.4} \end{equation} From $k=n-1$, we obtain the equation $2pa_{n-1}=n(pa_n)'$, which gives the well-known form of the weight $p=\exp(\frac{2}{n}\int\frac{a_{n-1}}{a_n}dx)/| a_n|$. Applying Leibniz rule again to $c_{j}=\sum_{k=j+1} ^n(-1)^{k-1-j}( b_{k}u) ^{(k-1-j)}$ $(0\leq j\leq n-1)$, we obtain in a similar manner, but this time on $\partial I$ (i.e. at the endpoints of $I$) the following equations for $1\leq j\leq n$, \begin{gather*} \binom{n-1}{j-1}b_n^{(n-j)}-\binom{n-2}{j-1}b_{n-1}^{(n-j-1)}+\dots +(-1)^{n-j}b_{j}=0\\ \binom{n-2}{j-2}b_n^{(n-j)}-\binom{n-3}{j-2}b_{n-1}^{(n-j-1)}+\dots +(-1)^{n-j}b_{j}=0\\ \dots\\ \binom{n-j}{0}b_n^{(n-j)}-\binom{n-j-1}{0}b_{n-1}^{(n-j-1)}+\dots +(-1)^{n-j}b_{j}=0 \end{gather*} This can be put in matrix form $A_{j} \begin{bmatrix} b_n^{(n-j)} & b_{n-1}^{(n-j-1)} & \dots & b_{j} \end{bmatrix} ^{T}=0$ $(1\leq j\leq n)$, where $A_{j}$ is the $j\times(n-j+1)$ matrix \[ A_{j}= \begin{bmatrix} \binom{n-1}{j-1} & -\binom{n-2}{j-1} & \dots & (-1)^{n-j}\\ \binom{n-2}{j-2} & -\binom{n-3}{j-2} & \dots & (-1)^{n-j}\\ \dots & \dots & \dots & \dots\\ \binom{n-j}{0} & -\binom{n-j-1}{0} & \dots & (-1)^{n-j} \end{bmatrix}. \] We thus have $(-1)^{j}\binom{j}{j}b_{j}+\dots+( -1) ^n\binom{n}{j}b_n^{(n-j)}=b_{j}$ on $I$ for $0\leq j\leq n$ and \[ A_{j} \begin{bmatrix} b_n^{(n-j)}\\ b_{n-1}^{(n-j-1)}\\ \dots\\ b_{j} \end{bmatrix} =0 \] on $\partial I$ for $1\leq j\leq n$. Conversely, it is clear that if $(-1)^{j}\binom{j}{j}b_{j} +\dots+( -1) ^n\binom{n}{j}b_n^{(n-j)}=b_{j}$ on $I$ for $0\leq j\leq n$ and \[ A_{j} \begin{bmatrix} b_n^{(n-j)}\\ b_{n-1}^{(n-j-1)}\\ \dots\\ b_{j} \end{bmatrix} =0 \] on $\partial I$, then equation \eqref{e2.2} holds for all functions $y$, $u$ in $C^{\infty}$, and therefore $L$ is self-adjoint. \end{proof} The following observation is particularly useful. When, in the above proof, $j\geq n-j+1$ i.e. $j\geq1+n/2$, we get more equations than ``unknowns'' $b_n^{(n-j)}$, $b_{n-1}^{(n-j-1)},\dots$, $b_{j}$. Thus, deleting the first $2j-n-1$ rows of $A_{j}$ and putting $k=n-j$ $(0\leq k\leq n/2-1)$, we obtain the equations $B_{k} \begin{bmatrix} b_n^{(k)} & b_{n-1}^{(k-1)} & \dots & b_{n-k} \end{bmatrix} ^{T}=0$ where $B_{k}$ is the $(k+1)\times(k+1)$ matrix \[ B_{k}= \begin{bmatrix} \binom{2k}{k} & -\binom{2k-1}{k-1} & \dots & (-1)^{k}\\ \binom{2k-1}{k} & -\binom{2k-2}{k-1} & \dots & (-1)^{k}\\ \dots & \dots & \dots & \dots\\ \binom{k}{k} & -\binom{k-1}{k-1} & \dots & (-1)^{k} \end{bmatrix}. \] Clearly, $\det B_{k}=\pm\det E_{k}$ where \[ E_{k}=\begin{bmatrix} \binom{0}{0} & \binom{1}{0} & \dots & \binom{k}{0}\\ \binom{1}{1} & \binom{2}{1} & \dots & \binom{k+1}{1}\\ \dots & \dots & \dots & \dots\\ \binom{k}{k} & \binom{k+1}{k} & \dots & \binom{2k}{k} \end{bmatrix}. \] As in the proof of Lemma \ref{lem2.4}, if we let \[ E(m,k)= \begin{vmatrix} \binom{m-k}{0} & \binom{m-k+1}{0} & \dots & \binom{m}{0}\\ \binom{m-k+1}{1} & \binom{m-k+2}{1} & \dots & \binom{m+1}{1}\\ \dots & \dots & \dots & \dots\\ \binom{m}{k} & \binom{m+1}{k} & \dots & \binom{m+k}{k} \end{vmatrix}, \] then elementary column operations give $E(m,k)=E(m,k-1)=\dots=E(m,1)=1$. Therefore $\det E_{k}=E(k,k)\neq0$, i.e. $B_{k}$ is non-singular, and we obtain \begin{equation} \label{e2.23} b_n^{(k)}=b_{n-1}^{(k-1)}=\dots=b_{n-k}=0\quad \text{on $\partial I$ for } 0\leq k\leq n/2-1. \end{equation} An interesting consequence of this is that if the weight $p=1$, then $a_n^{(k)}=a_{n-1}^{(k-1)}=0$ for $0\leq k\leq n/2-1$ on $\partial I$, and so, if $a_n$ is not constant, $I$ must be a finite interval $[\alpha,\beta]$ with $a_n=A(x-\alpha)^{n/2}(x-\beta)^{n/2}$ for some non-zero constant $A$ and $a_{n-1}=\frac{An^2}{2}(x-(\alpha+\beta)/2)(x-\alpha)^{n/2-1} (x-\beta)^{n/2-1}$ (recall that the degree of $a_{k}$ is at most $k $ and that $2pa_{n-1}=n(pa_n)'$). We thus obtain the form of the two leading polynomial coefficients of what may be considered as the $n$-th order Legendre differential equation. Proposition \ref{prop2.1} gives a necessary and sufficient set of conditions under which the operator $L$ is self-adjoint with respect to an inner product of the form $\langle u,v\rangle =\int_Ipuvdx$. This was achieved under the assumption that $p$ is an admissible weight, that is $\int_Ipudx$ is finite for all $C^{\infty}$ functions $u$, and that $p$ satisfies certain differentiability conditions. It is therefore highly desirable that we obtain conditions under which this assumption holds. The next section is devoted to such an analysis. \section{Analytic aspects of Sturm-Liouville theory} Keeping the same notation as before, let $ p=\exp\big(\frac{2}{n}\int\frac{b}{a}dx\big)/|a|$ be the weight function, where, for brevity, $a=a_n$ is a polynomial of degree at most $n$ and $b=a_{n-1}$ is a polynomial of degree at most $n-1$. Without loss of generality, we may assume $a$ to be monic. The weight function $p$ is a priori defined on an interval where $a$ does not vanish. However, on an interval $I$ where $a$ may have roots, it is clear that we can take $p$ to be defined piecewise, with $p=\exp(\frac{2}{n}\int\frac{b}{a}dx)/|a|$ except at the roots of $a$. Thus it is important to discuss the integrability of the weight function near the roots of $a$ and the differentiability at these roots. The integrability is basically a consequence of the following lemma. \begin{lemma} \label{lem3.1} For $\epsilon>0$ and $d, \alpha$ integers with $\alpha>0$, $\int_0^{\epsilon} \frac{e^{kx^{d}}}{x^{\alpha}}dx$ exists if and only if $k<0$ and $d<0$. \end{lemma} \begin{proof} If $d=0$ and $\alpha=1$, then clearly the integral is $+\infty$. So assume that $d>0$. Then for sufficiently small positive $x$, $1/21$, \[ 0<\eta<\epsilon \int_{\eta}^{\epsilon}\,x^{-\alpha}\,dx =\frac{\epsilon^{-\alpha+1}}{-\alpha +1}+\frac{1}{\alpha-1}\,\eta^{-(\alpha-1)}, \] so the integral $\int _0^{\epsilon}\frac{e^{kx^{d}}}{x^{\alpha}}dx$ is infinite. Hence a necessary condition for the integral to be finite is that $d<0$. Let us write $d=-\delta$ with $\delta>0$. Therefore, \[ \int_0^{\epsilon} \frac{e^{kx^{d}}}{x^{\alpha}}dx=\int_0^{\epsilon} \frac{e^{k( \frac {1}{x}) ^{\delta}}}{x^{\alpha}}dx. \] If $k=0$, this is integrable if and only if $-\alpha+1>0$, which is not the case. If $k>0$, then \[ \int_0^{\epsilon}\frac{e^{k^{(\frac{1}{x})^{\delta}}}}{x^{\alpha}}\,dx =\int_{\infty}^{1/\epsilon} u^{\alpha} e^{ku^{\delta}}(-1) u^{-2}\,du =\int_{1/\epsilon}^{\infty} u^{\alpha-2}\,e^{k u^{\delta}}\,du >\int_{\frac{1}{\epsilon}+1}^{\infty}e^{ku} u^{\alpha-2}\,du \] (as $\delta$ is a positive integer). If $\alpha-2\geq0$, this is clearly infinite. For the remaining case $\alpha=1$, substituting $ku=v$, the integral $\int_{\frac{1}{\epsilon}+1}^{\infty}e^{ku} u^{\alpha-2}\,du$ becomes $\int_{k( \frac {1}{\epsilon}+1) }^{\infty}\frac{e^{v}}{v}dv$. As $\frac{e^{v}}{v}>\frac{v}{2}$ (for $v>0$), this integral is clearly infinite. Thus for $\epsilon>0$, a necessary condition for $\int_0 ^{\epsilon}\frac{e^{kx^{d}}}{x^{\alpha}}dx$ to be finite is that $k<0$, $d<0$. It remains to show that in this case, the above integral is finite. Let $d=-\delta$, $k=-l$, $\delta>0$, $l>0$. So $\frac{e^{kx^{d}}}{x^{\alpha} }=u^{\alpha}e^{-lu^{\delta}}$, where $u=\frac{1}{x}$. Now $\lim_{u\to\infty}\frac{e^{lu^{\delta}}}{u^{\alpha} }=\infty$, so $\lim_{u\to\infty}u^{\alpha}e^{-lu^{\delta }}=0$. Thus $\int_0^{\epsilon}\frac{e^{kx^{d}}}{x^{\alpha}}dx$ is bounded near $0$. Hence the function $f(\eta)=\int_{\eta}^{\epsilon} \frac{e^{kx^{d}}}{x^{\alpha}}dx$, with $0<\eta<\epsilon$ is monotonic bounded, hence its limit as $\eta\to0^{+}$ exists. This completes the proof of the lemma. \end{proof} \begin{corollary} \label{coro1} For $\epsilon>0$, $d, \alpha$ integers with $\alpha>0$, $\int_{-\epsilon}^{0} \frac{e^{kx^{d}}}{|x^{\alpha}|}dx$ exists if and only if $k(-1)^{d}<0$ and $d<0$. \end{corollary} The proof of the above corollary follows by substituting $u=-x$ and applying Lemma \ref{lem3.1}. We say that a function $f(x)$ defined and continuous on an open interval $I$ containing $0$ as a left end point is left integrable at $0$ if for any $\eta\in I$, $ \lim_{\epsilon\to0^{+}} \int_{\epsilon }^{\eta}f(x)dx$ exists. Similarly if $f(x)$ is a function defined and continuous on an an open interval $I$ containing $0$ as a right end point is right integrable at $0$ if for any $\eta\in I$, $ \lim_{\epsilon\to0^{-}} \int_{\eta}^{\epsilon}f(x)dx$ exists. Clearly this is equivalent to saying that the function $g(x)=f(-x)$ is left integrable at $0$. By suitable translations, one can define the concept of left and right integrability at the end points of an interval $I$ on which the given function is defined and continuous. Let $r$ be a zero of $a(x)=a_n(x)$ and let $m_{a}(r)=\alpha$ and $m_{b}(r)=\beta$ be the multiplicities of $r$ as a root of $a(x)$ and $b(x)=b_{n-1}(x)$. Thus $\frac{b(x)}{a(x)}=(x-r)^{\beta-\alpha}\phi(x)$ where $\phi(x)$ is a rational function with $\phi(r)\ne0$. Hence $\frac{b(x)}{a(x)}=\phi(r) (x-r)^{\beta-\alpha}\psi(x)$ where $\psi(r)=1$. \smallskip \noindent\textbf{Definition.} We say that a root $r$ of $a(x)$ is an ordinary root if $m_{b}(r)-m_{a}(r)+1\ne0$, and it is a logarithmic root if $m_{b}(r)-m_{a}(r)+1 = 0$. Using Lemma \ref{lem3.1} we have the following result which is one of the main tools for explicit determination of self-adjoint operators. \begin{proposition} \label{prop3.2} Let $\beta=m_{b}(r),\alpha =m_{a}(r)$, so that $\frac{b(x)}{a(x)}=(x-r)^{\beta-\alpha}\phi(x)$, where $\phi(x)$ is a rational function with $ \phi(r)=\lim_{x\to r}(x-r)^{\alpha-\beta}\frac{b(x)}{a(x)} \neq 0$. Then \begin{itemize} \item[(i)] For an ordinary root $r$ of $a$, the weight function $p(x)$ is integrable from the right at $r$ if and only if $\alpha-\beta\geq2$ and $\phi(r)>0$. It is integrable from the left at $r$ if and only if $\alpha -\beta\geq2$ and $(-1)^{\alpha-\beta}\phi(r)<0$. In this case, the weight function $p(x)$ is respectively right/left $C^{\infty}$ differentiable at $r$ and $p$ and all its (one sided) derivatives vanish at $r$. \item[(ii)] For a logarithmic root $r$ of $a$, the weight function $p(x)$ is right/left integrable near $r$ if and only if $\frac{|x-r|^{\frac{2}{n}\phi(r)}}{|x-r|^{\alpha}}$ is integrable near $r$ if and only if $\frac{2}{n}\phi(r)-\alpha+1>0$. \end{itemize} \end{proposition} \begin{proof} Using the notation in this proposition, we can write \begin{gather*} a(x)=(x-r)^{\alpha}(a_0 + a_1 (x-r) + \dots), \\ b(x)=(x-r)^{\beta}(b_0 + b_1 (x-r) + \dots), \end{gather*} where $a_0$, $b_0$ are not zero. For convenience of notation, we may suppose that $r=0$. The weight function, near the root $r=0$ can then be written as \[ p(x)=\frac{\exp(\int\frac{2}{n}\frac {b(x)}{a(x)}\,dx)}{|a_0||x^{\alpha}|}\,\psi(x), \] where $\psi(0)=1$ and $\psi$ is infinitely differentiable near $0$. Hence the order of left or right differentiability at $0$ of $p(x)$ is the same as that of $\frac {\exp(\int\frac{2}{n}\frac{b(x)}{a(x)}\,dx)}{|x^{\alpha}|}$. We will prove differentiability at $0$ for all orders after the proof of Proposition \ref{prop3.3}. For the integrability, we may suppose that near $0$, $1/2<\psi(x)<3/2$ and therefore near $0$, we have the estimate \[ \frac{1}{2}\frac{e^{\int\frac{2}{n}\frac{b(x)}{a(x)}\,dx}} {|a_0||x^{\alpha}|}0$ and $\alpha- \beta\geq2$. It is integrable from the left if and only if $\frac{d_0 (-1)^{\beta-\alpha+1} }{\beta-\alpha+1}<0$ and $\beta-\alpha+1 <0$. Thus, the requirement for left integrability becomes $d_0 (-1)^{\beta-\alpha+1}>0$ and $\alpha- \beta\geq 2$. Hence the weight is both left and right integrable if and only if $d_0>0$, $\beta-\alpha+1<0$, $(-1)^{\beta-\alpha+1}=1$. Thus $\beta -\alpha+1$ is an even negative integer and $d_0 >0$ are the requirements for both integrability from the left and right. Recall that $\frac{b(x)}{a(x)}= d_0 x^{\beta- \alpha} \psi(x)$, with $\psi(0)=1$. Thus $x^{\alpha-\beta} \frac{b(x)}{a(x)}=d_0 \psi(x)$. Hence $ \lim_{x \to0} x^{\alpha-\beta} \frac{b(x)}{a(x)}=d_0$. This gives part (i) of the Proposition - as far as integrability is concerned. Now assume, with the same notation as above and with $r=0$, that $\beta-\alpha+1=0$. In this case, we need to investigate the integrability of $\exp(\int_{p_0}^{x}\frac{2}{n}t^{-1}\psi(t)\,d)/|x^{\alpha}|$. As $\psi(0)=1$, given any $\epsilon>0$, for sufficiently small $x$, $1-\epsilon<\psi(x)<1+\epsilon$. Therefore \[ \frac{2}{n}\,d_0(1-\epsilon)|x|^{-1}<\frac{2}{n}\,d_0 \,|x|^{-1}<\frac{2}{n}\,d_0\,(1+\epsilon)|x|^{-1}. \] Integrating - from the right near $0$, we get \[ K_1\frac{e^{\frac{2}{n}d_0(1-\epsilon)\,ln\,|x|}} {|x^{\alpha}|}<\frac{e^{\int_{p_0}^{x}\frac{2}{n}\,d_0\,t^{-1}\psi(t)dt} }{|x^{\alpha}|}< K_2\frac{e^{\frac{2}{n}d_0(1+\epsilon )\,ln\,|x|}}{|x^{\alpha}|}, \] where $K_1, K_2$ are positive constants. This gives \[ K_1\,|x|^{\frac{2}{n}d_0\,(1-\epsilon)-\alpha}<\frac {e^{\int_{p_0}^{x}( \frac{2}{n}\frac{b(t)}{a(t)}) \,dt} }{|x^{\alpha}|}<2K_2\,|x|^{\frac{2}{n}d_0(1+\epsilon)-\alpha.} \] If the weight function is right integrable near $0$, then necessarily $\frac{2}{n}d_0(1-\epsilon)-\alpha+1>0$. Hence $\frac{2}{n}d_0-\alpha +1>0$. If this holds then the displayed inequalities above establish the integrability of the weight function. Similar arguments give the same condition for integrability from the left- namely $\frac{2}{n}d_0 -\alpha+1 >0$. This completes the proof of the proposition, except for the differentiability of the weight, which is discussed after Proposition \ref{prop3.3} \end{proof} Using lower and upper bounds on the asymptotic form of the weight function (as $x\to\infty$), or partial fraction decomposition of $b/a$, we have the following result. \begin{proposition} \label{prop3.3} (i) If $a$ has no real roots and $p(x)=\exp(\frac{2}{n}\int_{p_0}^{x} \frac{b(t)}{a(t)}dt)/|a(x)|$ then the weight function $p(x)$ gives finite norm for all polynomials if and only if $\deg b-\deg a$ is an odd positive integer and the leading term of $b$ is negative. (ii) If $a$ has only one root, say $0$, then $p(x)=\exp(\frac{2}{n}\int_{p_0}^{x} \frac{b(t)}{a(t)}dt)/|a(x)|$ gives a finite norm for all polynomials on $(0,\infty)$ if and only if \begin{itemize} \item[(a)] $\deg b-\deg a\geq0$ and the leading term of $b$ is negative. \item[(b)] If $a=x^{\alpha}( A_0+A_1x+\dots)$ and $b=x^{\beta}( B_0+B_1x+\dots)$, where $A_0$ and $B_0$ are nonzero constants, then $\alpha-\beta\geq1$, and $\frac{B_0}{A_0}>0$ for $\alpha-\beta\geq2$ whereas $\frac{2B_0}{nA_0}-\alpha+1>0$ for $\alpha-\beta=1$. \end{itemize} \end{proposition} \begin{proof} Assume that $a(x)$ has no real roots. We may assume that the leading term of $a(x)$ is $1$. Thus \[ a(x)=x^n+a_{n-1}x^{n-1}+\dots,\quad b(x)=kx^{m}+\dots, \] where $n, m$ are the degrees of $a,\ b$. Thus $\frac{b(x)}{a(x)}=kx^{m-n}\psi(x)$, with $\lim_{x\to\infty} \psi(x)=1$. Also $a(x)=x^n\phi(x)$ and $\lim_{x\to\infty}\phi(x)=1$. Hence for sufficiently large positive $x$, there are positive constants $d_1, d_2, c_1, c_2$ with \begin{equation} \frac{1}{|x|^nd_1}\,e^{\int_{p_0}^{x}\frac{2}{n} kc_1t^{m-n}dt}0$, $e^{\int_{p_0}^{x}\frac{2}{n}kct^{m-n} dt}=A|x|^{\frac{2}{n}kc}$. Therefore $\int_{M}^{\infty}Ax^{N}\frac{x^{\frac{2}{n}kc}}{x^n}dx$ cannot be finite if $N$ is large enough. Thus $m-n\neq-1$. Therefore from \eqref{east} we obtain the estimate \[ \frac{A_1}{|x|^n\,d_1}\,e^{\frac{2}{n}kc_1\frac {x^{m-n+1}}{m-n+1}}0$, where $k$ is the leading term of $b(x)$: recall that we have taken the leading term of $a(x)$ to be $1$. Similarly if we require finiteness of the integrals $\int_{-\infty}^{-M}x^{N} \,e^{\frac{2}{n}kc\frac{x^{m-n+1}}{m-n+1}}\,dx$, then the requirements are $\frac{k(-1)^{m-n+1}}{m-n+1}<0$, $-(m-n+1)<0$ Therefore the conditions for $p(x)=\exp(\int\frac{2}{n}\frac{b(x)}{a(x)}dx)/|x^{\alpha}|$ to be a weight on $(-\infty,\infty)$ are that $k<0$, and $(m-n+1)$ should be an even positive number. This completes the proof of (i). Part (ii) follows from the integrability of the weight on $(0,\infty)$ and Proposition \ref{prop3.2}. \end{proof} \subsection*{Differentiability properties of the weight function} The differentiability properties of the weight function at zeroes of the leading term $a(x)$ follow from the following observation: if $d$ is a positive integer and $k$ is a positive number, then $ \lim_{x \to0^{+}} \frac{e^{-kx^{-d}}}{P(x)}=0$ for all polynomials P(x). Let $N$ be the order of $0$ in $P(x)$. Then $P(x)$ can be written as $P(x)=ax^{N} Q(x)$ where $Q$ is a polynomial with $Q(0)=1$. Thus it suffices to show that $ \lim_{x \to0^{+}} \frac{e^{-kx^{-d}}}{x^{N}}=0$. This is the same as $ \lim_{u \to\infty} \frac{u^{N}}{e^{ku^{d}}}=0$, which is obviously true. Therefore $ \lim_{x \to0^{+}} \frac {e^{-kx^{-d}}}{P(x)} Q(x)=0$ for polynomials $P, Q$ with $P\ne0$. As in the proof of proposition \ref{prop3.2}, the weight function, near the root $r=0$ can then be written as $p(x)=\psi(x) \exp(\int\frac{2}{n} \frac{b(x)}{a(x)}dx)/(|a_0| |x^{\alpha}|)$ with $\psi(0)=1$. Thus, if $\alpha, \beta$ are the multiplicities of the root $0$ in $a, b$ respectively and $\beta- \alpha+1 \ne0$, we have the estimate \[ \frac{A_1e^{(k_1\frac{2}{n}d_0)( \frac{x^{\beta-\alpha+1}} {\beta-\alpha+1}) }}{|x^{\alpha}|}\leq\frac{e^{\int^{x}_{p_0}( \frac{2}{n}\,\frac{b(x)}{a(x)}) \,dt}}{|x^{\alpha}|} \leq\frac {A_2e^{(k_2\frac{2}{n}d_0)( \frac{x^{\beta-\alpha+1}}{\beta -\alpha+1}) }}{|x^{\alpha}|} \] where $A_1, A_2, k_1, k_2$ are positive constants. We discuss the right hand limit, as the left hand limit is treated similarly. As the weight is integrable, we must have $\beta-\alpha+1<0$, $d_0>0$. Thus, by the above observation, $ \lim_{x \to0^{+}}p(x)=0$, where $p(x)$ is the weight function. Moreover, by the same observation, $ \lim_{x \to0^{+}}p(x) R(x)=0$ for any rational function. Let $\phi(x)=\int^{x}_{p_0}( \frac{2}{n} \frac {b(x)}{a(x)} ) dt $. Therefore $p(x)=e^{\phi(x)}/|a(x)|$. Then all derivatives of $p(x)$ are of the form $e^{\phi}(x) R(x)$, where $R(x)$ is a rational function. Therefore the right-hand limits at $0$ of all the derivatives of the weight function are $0$. In case $\beta-\alpha+1 = 0$, the weight can be written near a zero of the leading term - by change of notation- as \[ p(x) = \frac{c}{|a|} |x|^{( \frac{2}{n}d_0 -\alpha) } \frac{e^{\phi(x)}}{1+ \psi(x)}, \] where $\phi,\ \psi$ are analytic functions near $0$, and \begin{gather*} a(x)=(x)^{\alpha}(a_0 + a_1 (x) + \dots), \\ b(x)=(x)^{\beta}(b_0 + b_1 (x) + \dots), \end{gather*} where $a_0, b_0$ are not $0$ and $d_0=\frac{b_0}{a_0}$. Here $\frac{2}{n} d_0-\alpha+1>0$ - by the assumption of integrability of the weight. This completes our discussion of the differentiability properties of the weight function near a zero of the leading term. The differentiability properties of ordinary roots have already been discussed. We now assume that the multiplicity of a root $r$ of $a(x)=a_n (x)$ is $\alpha$ and its multiplicity in $b(x)=a_{n-1}(x)$ is $\beta$. For convenience of notation we assume that $r$ is zero. As above, we have $a_n(x)=x^{\alpha}(A_0 + A_1 x+ \dots)$, $a_{n-1}(x)=x^{\beta}(B_0 + B_1 x+ \dots)$ with $(\beta-\alpha)=-1$. Thus near $x=0$, the weight is of the form $p(x)=\frac{1}{|A_0|}|x|^{(\frac{2}{n}\frac{B_0}{A_0}-\alpha) } \frac{e^{\phi(x)}}{1+\psi(x)}$ where $\phi$, $\psi$ are analytic near zero and $\psi(0)=0$. When there is no danger of confusion we will write $p(x)\sim|x|^{( \frac{2}{n}\frac{B_0}{A_0}-\alpha) }$. Now $p'=p( \frac{2}{n}\frac{b}{a}-\frac{a'}{a} ) $. Therefore, all higher derivatives of $p$ are of the form $p \rho$ where $\rho$ is a rational function and all higher derivatives of $p \rho$ are also multiples of $p$ by rational functions. For later use we record the asymptotic behavior of $p'$ near a zero of $a_n$. \[ p'(x)=\frac{1}{|A_0|}|x|^{( \frac{2}{n}\frac{B_0}{A_0 }-\alpha) } \frac{e^{\phi(x)}}{1+\psi(x)} \Big( \frac{2}{n}\frac {b}{a}-\frac{a'}{a} \Big) \] The weight $p$ is integrable near zero if and only if $( \frac{2} {n}\frac{B_0}{A_0}-\alpha+1) >0$. Moreover $ \lim _{x\to0}p(x) a_n(x) =0$ if and only if $\frac{2}{n}\frac{B_0}{A_0} >0$. By the integrability of the weight $\frac{2}{n}\frac{B_0}{A_0} >\alpha-1 \geq0$. Thus the boundary condition $ \lim_{x\to0}p(x) a_n(x)=0$ is a consequence of the integrability of the weight near zero. Similarly $p(x)a_{n-1}(x)=\frac{1}{|A_0|}|x|^{( \frac{2}{n}\frac{B_0}{A_0 }-1) } \frac{e^{\phi(x)}}{1+\psi(x)} (B_0 + B_1 x + \dots)$, keeping in mind that $\alpha-\beta= 1$. Hence $ \lim_{x\to0}p(x) a_{n-1}(x)=0$ if and only if $( \frac{2}{n}\frac{B_0}{A_0}-1)>0$. \subsection{Higher order operators} The principal aim of this section is to prove the following result. \begin{proposition} \label{prop3.4} Let $L=a_n(x) y^{(n)}+a_{n-1}(x)y^{(n)} +\dots+ a_2(x) y''+a_1(x) y'$ be a self-adjoint operator of order $n$ with $n>2$. If $a_n$ has a real root then the multiplicity of the root is at least $2$ and the multiplicity of the same root in $a_{n-1}$ is positive and less than its multiplicity in $a_n$. \end{proposition} \begin{proof} Let $r$ be a real root of $a_n$ and assume that it is a simple root. It is then a logarithmic root. Therefore, near $r$, we have \[ p(x)\sim|x-r|^{( \frac{2}{n}\frac{B_0}{A_0}-1) }, \] where $a_n(x)=(x-r) (A_0 + A_1 (x-r)+ \dots)$, $a_{n-1}(x)= (B_0 +B_1 (x-r)+ \dots)$ and $A_0$, $B_0$ are not zero. Now $p'=p\big( \frac{2}{n}\frac{a_{n-1}}{a_n}-\frac{a_n'}{a_n} \big) $. Therefore \[ p'(x)\sim|x-r|^{( \frac{2}{n}\frac{B_{0}}{A_{0}}-1)} \Big( \frac{2}{n}\frac{a_{n-1}(x)}{a_{n}(x)}-\frac{a_{n}'(x)} {a_{n}(x)}\Big). \] Similarly \[ p(x) a_{n-1}(x)\sim|x-r|^{( \frac{2}{n}\frac{B_0}{A_0}-1) } (B_0 + B_1 (x-r)+ \dots) \] The boundary conditions and the determining equations in Proposition \ref{prop2.1} imply that $(a_n p)$, $(a_{n-1} p)$ and $(a_{n-1}p)'$ vanish on the boundary. Now $ \lim_{x\to r}a_n(x)p(x)=0$ is a consequence of the integrability of the weight near $r$. Similarly $ \lim_{x\to r}a_{n-1}(x)p(x)=0$ if and only if $( \frac{2}{n}\frac{B_0}{A_0}-1) >0$. Let $ l_{r} =\lim_{x\to r} (x-r)^{\alpha-\beta} \frac{a_{n-1}(x)}{a_n(x)}$. Clearly $l_{r}=\frac{B_0}{A_0}=\frac{a_{n-1}(r)}{a_n'(r)}$, as $\alpha-\beta=1$. Since $ \lim_{x\to r}pa_{n-1}=0$ and $a_{n-1}(r)\ne0$ we see that $p$ must vanish at $r$ in the sense that its limit at $r$ is zero. The boundary condition $ \lim_{x\to r}(a_{n-1}p)'=0$ now implies that $ \lim_{x\to r}p'(x)=0$. Now $p'=p( \frac{2}{n}\frac{a_{n-1}}{a_n}-\frac{a'_n}{a_n}) $. Thus near the root $r$, \[ p'\sim|x-r|^{(\frac{2}{n}l_{r}-2\alpha)} \Big( \frac{2}{n}a_{n-1}-a_{n}'\Big). \] If $a_{n-1}-\frac{n}{2}a'_n\equiv0$ then in particular $( \frac{2}{n}\frac{a_{n-1}(r)}{a'_n(r)}-1)=0$. This means that $ \lim_{x\to r}(x-r)\frac{2}{n} \frac{a_{n-1}(x)}{a_n(x)}-1=0$ i.e. $( \frac{2}{n}\frac{B_0}{A_0}-1) =0$. This contradicts the boundary condition $ \lim_{x\to r}a_{n-1}(x)p(x)=0$. Let $(a_{n-1} - \frac{n}{2}a'_n)=(x-r)^{\lambda} H(x)$ where $\lambda\geq0$ and $H(r)\ne0$. If $\lambda> 0$ then $ \lim_{x\to r}(x-r)\frac{2}{n}\frac{a_{n-1}(x)}{a_n(x)}-1=0$ which again contradicts the boundary condition $ \lim_{x\to r}a_{n-1}(x)p(x)=0$. Hence $p'\sim|x-r|^{(\frac{2}{n}l_{r}-2\alpha)}H(x)$ so $p'\to0$ at $r$ if and only if $( \frac{2}{n}l_{r}-2\alpha) >0$. By Proposition \ref{prop2.1}, the operator must satisfy - beside other equations - the determining equations \begin{gather} \label{det1-gen} n(a_np)'=2(a_{n-1}p), \\ \label{det2-gen}\frac{( n-1) (n-2)}{6}(a_{n-1}p)^{\prime\prime }-( n-2) (a_{n-2}p)'+2(a_{n-3}p)=0 \end{gather} Equation \eqref{det2-gen} is equivalent to \begin{equation} \label{rational-gen} \begin{aligned} &C_1( \frac{a_{n-1}}{a_n}) ^3 +C_2( \frac{a_{n-1}}{a_n}) ( \frac{a_{n-1}}{a_n}) ' + C_3( \frac{a_{n-1}}{a_n}) ''\\ &+C_4 ( \frac{a_{n-2}}{a_n}) '+ C_5 \frac{a_{n-1}}{a_n}\frac{a_{n-2}}{a_n}+ C_{6}\frac{a_{n-3}}{a_n}=0 \end{aligned} \end{equation} where \begin{gather*} C_1=\frac{2(n-1)(n-2)}{3n^2},\quad C_2=\frac{(n-1)(n-2)}{n}, \quad C_3=\frac{(n-1)(n-2)}{6},\\ C_4= -(n-2),\quad C_5=-\frac{2(n-2)}{n}, \quad C_{6}=2. \end{gather*} This implies the identity \begin{equation}\label{divide-general} a_{n-1}\big( a_{n-1}-na_n'\big) \big(a_{n-1}-\frac{n}{2}a'_n\big) \equiv0 \quad \mod{a_n} \end{equation} Using this identity, as $(x-r)$ divides $a_n$ but it does not divide $a_{n-1}$ nor $( a_{n-1}-\frac{n}{2}a'_n) $, it must divide $(a_{n-1}-na_n')$. But then $ \lim_{x\to r}(a_{n-1}-na_n') =0$. This means that $ \lim_{x\to r}\frac{2}{n}\frac{a_{n-1} }{a'_n}-2=0$ i.e. $\frac{2}{n}\frac{B_0}{A_0}-2=0$. As seen above $p'\to0$ at $r$ if and only if $( \frac{2}{n}l_{r} -2\alpha) >0$. Since $\alpha=1$ we have a contradiction. Therefore $a_n$ cannot have a simple root and its multiplicity $\alpha$ in $a_n$ is at least $2$. Suppose that the multiplicity $\beta$ of $r$ in $a_{n-1}$ is zero. By considering the order of poles of $a_n$ in \eqref{rational-gen} we see that $\beta$ cannot be zero. This completes the proof. \end{proof} This result has an important consequence for fourth order self-adjoint operators. \begin{corollary} \label{coro3.5} Let $L$ be a self-adjoint operator of order $4$ and $a_4$ be its leading term. If $a_4$ has more than one real root then it has exactly two real roots with multiplicity $2$. Moreover the multiplicity of each real root of $a_4$ in $a_3$ is $1$. \end{corollary} \begin{proposition} \label{prop3.6} Let $n>2$ and suppose that $a=a_n$ has at most one real root. Then $2\deg b-\deg a\leq n-2$ or $3\deg b-2\deg a\leq n-3$, where $b=a_{n-1}$. \begin{itemize} \item[(i)] If $a$ has no real root then $\deg a<\deg b\leq n-3$; \item[(ii)] Suppose $a$ has only one real root $r$ with multiplicity $\alpha$, let $\beta$ be the multiplicity of $r$ as a root of $b$, and let $a=(x-r) ^{\alpha}u$, $b=( x-r) ^{\beta}v$. Then \[ 2\leq\deg a\leq\deg b\leq n-2\quad \text{and}\quad 1+\deg u\leq\deg v\leq n-3 \] \end{itemize} \end{proposition} \begin{proof} First, in all cases, $\deg b\geq1$. This is because if $a$ has no real root then $\deg b$ is odd and if $a$ has (at least) one real root then this will also be a root for $b$ (by Proposition \ref{prop3.4}). If we multiply by $a^3$ both sides of \eqref{rational-gen}, then the six terms on the left-hand side will be polynomials with respective degrees \begin{gather*} 3\deg b,\quad 2\deg b+\deg a-1,\quad 2\deg a+\deg b-2,\\ 2\deg a+\deg a_{n-2}-1,\quad \deg a+\deg b+\deg a_{n-2},\quad 2\deg a+\deg a_{n-3} \end{gather*} A comparison of degrees shows that $a_{n-2}$ and $a_{n-3}$ cannot be both zero and that \[ 2\deg b \leq\deg a+\deg a_{n-2},\quad \text{or}\quad 3\deg b \leq2\deg a+\deg a_{n-3} \] Using the fact that $\deg a_{j}\leq j$, we obtain \[ 2\deg b-\deg a\leq n-2\quad\text{or}\quad 3\deg b-2\deg a\leq n-3 \] If $a$ has no real roots then, by Proposition \ref{prop3.3} (i), $\deg b-\deg a\geq1$ and hence \[ \deg a<\deg b\leq n-3 \] If $a$ has only one real root $r$ with multiplicity $\alpha$, then $\alpha\geq2$ and $b$ has $r$ as a root with multiplicity $\beta$, where $1\leq\beta<\alpha$ (by Proposition \ref{prop3.4}). Since $\deg b\geq\deg a$, we obtain that \[ 2\leq\deg a\leq\deg b\leq n-2 \] Let $a=( x-r) ^{\alpha}u$, $b=( x-r)^{\beta}v$. Then $\deg a=\alpha+\deg u$ and $\deg b=\beta+\deg v$, and we obtain $\deg v-\deg u\geq\alpha-\beta\geq1$. Now $\beta+\deg v\leq n-2$, so $1+\deg u\leq\deg v\leq n-3$ and thus $\deg u\leq n-4$. \end{proof} \section{Examples of higher order operators, their eigenvalues and orthogonal eigenfunctions} Let $L$ be an operator of the form $L(y)=\sum_{k=0}^na_{k}(x)y^{(k)}$, where $\deg a_{k}\leq k$; then the eigenvalues of $L$ are the coefficients of $x^n$ in $L(x^n)$, $n=0,1,2,\dots$. \begin{proposition} \label{prop4.1} Let $L$ be a linear operator that maps the space $\mathbb{P}_n$ of all polynomials of degree at most $n$ into itself for all $n\leq N$. If the eigenvalues of $L$ are distinct or if $L$ is a self-adjoint operator then there is an eigenpolynomial of $L$ in every degree $\leq N$. \end{proposition} This means that if two operators leave the space of polynomials of degree at most $n$ invariant for all $n$ and the weight function which makes the two operators self-adjoint is the same, then they have the same eigenfunctions. The eigenvalues in general are not simple. The proof of the above proposition is left to the reader. Let $\lambda$ be an eigenvalue of $L$ and $\mathbb{P}_n(\lambda)$ the corresponding eigenspace in the space $\mathbb{P}_n$ of all polynomials of degree less than or equal to $n$. If $n_0$ is the minimal degree in $\mathbb{P}_n(\lambda)$, then there is, up to a scalar only one polynomial in $\mathbb{P}_n(\lambda)$ of degree $n_0$. Choose a monic polynomial $Q_1$ in $\mathbb{P}_{n_0}(\lambda)$. Let $n_1$ be the smallest degree, if any, greater than $n_0$ in $\mathbb{P}_n(\lambda)$. The codimension of $\mathbb{P}_{n_0}(\lambda)$ in $\mathbb{P}_n(\lambda)$ is $1$. Therefore, in the orthogonal complement of $\mathbb{P}_{n_0}(\lambda)$ in $\mathbb{P}_{n_1}(\lambda)$, choosing a monic polynomial $Q_2$, which will necessarily be of degree $n_1$, the polynomials $Q_1,\ Q_2$ give an orthogonal basis of $\mathbb{P}_{n_1 }(\lambda)$. Continuing this process, we eventually get an orthogonal basis of $\mathbb{P}_n(\lambda)$ consisting of monic polynomials. We illustrate this by an example of a fourth order self-adjoint operator that has repeated eigenvalues but which is not an iterate of a second order operator. Consider the operator \begin{equation}\label{op-special-eg} L=(1-x^2)^2 y^{(4)}-8x(1-x^2) y'''+ 8 y''-24x y'. \end{equation} Its eigenvalues are $\lambda_n = n[ (n-1)(n-2)(n+5)-24]$. The eigenvalue $\lambda= -24$ is repeated in degrees $n=1$ and $n=3$. The weight function for which this operator is self-adjoint is $p(x)=1$. The eigenpolynomials of degree at most $3$ are \[ y_0(x)=1,\quad y_1(x)=x,\quad y_2(x)=x^2-\frac{1}{3}, \quad y_3(x)=x^3. \] This gives the set of orthogonal polynomials $\{1,x,x^2-\frac{1}{3} ,x^3-\frac{3}{5}x\}$. Since the weight function is the same as that for the classical Legendre polynomials, this family up to scalars is the same as the corresponding classical Legendre polynomials. We now return to examples of higher order operators. The restrictions on the parameters appearing in all the examples come from integrability of the weight and boundary conditions. Before giving examples of higher order operators, it is instructive to consider the classical case of second order operators in the frame work of section 3. \subsection{Self-adjoint operators of order $2$} Assume $n=2$ and that $a_2 (x)$ has distinct roots, which we may assume to be $-1$ and $1$. If $\alpha$ is the multiplicity of a root $r$ of $a_2 (x)$ and $\beta$ is its multiplicity in $a_1 (x)$ then the integrability of the associated weight gives the equation $\alpha=\beta+ 1 + \delta$, with $\delta\geq0$. As $\alpha=1$, we must have $\beta=0$, $\delta=0$. Thus, only the logarithmic case can occur. Let $a_1(x)= cx+d$. The integrability condition at a root $r$ reads $ \lim_{x\to r}(x-r) a_1(x)/ a_2(x)>0$. As we are assuming that $a_2 (x)=x^2-1$, the integrability conditions at both the roots gives $c+d>0$, $-c+d<0$, so $-c0$. \subsection{Self-adjoint operators of order $4$} In this section we determine all self-adjoint operators \begin{equation}\label{main-op-n-4} L=a_4(x) y^{(4)}+a_3(x) y'''+ a_2(x) y''+a_1(x) y', \end{equation} with an admissible weight $p(x)=\frac{e^{\frac{1}{2}\int\frac{a_3(x)} {a_4(x)}dx}}{| a_4(x)| }$, satisfying the differential equation \begin{equation}\label{main-eq-n-4} L(y)=\lambda y. \end{equation} By Proposition \ref{prop2.1}, the operator $L$ must satisfy the determining equations \begin{gather} (a_4p)'= \frac{1}{2} (a_3 p)\label{det1-n-4}\\ (a_3p)''- 2(a_2p)'+ 2(a_1p)=0 \label{det2-n-4} \end{gather} on $I$, subject to the vanishing of \begin{equation}\label{boundary-n-4} (a_4 p),\quad (a_3 p) \quad\text{and}\quad ((a_2p) - \frac{(a_3p)'}{2}) \end{equation} on the boundary $\partial I$. \subsection*{Case: $a_4(x)$ has no real roots} If $a_4(x)$ is a monic polynomial having no real roots then $I = (-\infty,\infty)$ and, by Proposition \ref{prop3.6}, we have $a_4(x)=1$, $a_3(x)$ linear. Considering $a_4(x)=1$, $a_3(x)$ linear and solving determining equations \eqref{det1-n-4}, \eqref{det2-n-4} subject to the constraints \eqref{boundary-n-4} determines the fourth order self-adjoint operators \eqref{main-op-n-4} and the differential equations \eqref{main-eq-n-4} as \begin{gather*} a_4(x)=1, \quad a_3(x) = 2 (m_1-2 m_2^2 x),\\ a_2(x) = 4 m_2^{4} x^2-4 m_1 m_2^2 x+{A}, \quad a_1(x) = (-m_1^2+2 m_2^2+{A}) (m_1-2 {m_2}^2 x) \end{gather*} with the weight function \[ p(x)= e^{-m_2^2 x^2 +{m_1}x + {m_0}} \quad (m_2 \ne0) \] and the eigenvalues \[ \lambda_n = 2 m_2^2 ( m_1^2-{A}+2 m_2^2 (n-2)) n . \] Depending on the choice of $A, m_1, m_2$ one can get repeated eigenvalues. The specific case $m_1=0, m_2=1, A=-4$ gives the standard fourth order Hermite operator \cite{littlejohn-krall89} \[ a_4(x)=1, \quad a_3(x) =-4x, \quad a_2(x)=4 (x^2-1), \quad a_1(x)=4x \] with the weight $p(x)= e^{-x^2}$ and non-repeated eigenvalues $\lambda_n = 4n^2$. It is worth noticing that this specific case is the only operator in the class of fourth order operators with $a_4(x)=1$ that is an iterate of the second order Hermite operator, and in general, this class is not obtainable as iteration of the second order case. \subsection*{Case: $a_4(x)$ has only one real root} In this case we have, by Proposition \ref{prop3.6}, $a_4(x)=x^2$ and $a_3 (x)=x(a+bx)$ with $b\ne0$. The weight is determined as $p(x) = |x|^{\frac {a}{2}-2} e^{\frac{b}{2}x}$ with $b<0$ and $I = (0,\infty)$. Solving the determining equation \eqref{det2-n-4} subject to the constraints \eqref{boundary-n-4} determines the fourth order self-adjoint operators \eqref{main-op-n-4} and the differential equations \eqref{main-eq-n-4} as \begin{gather*} p(x) = |x|^{\frac{a}{2}-2} e^{\frac{b}{2}x} \quad \text{with } a>2,\; b<0\\ a_4(x) = x^2, \quad a_3(x) = x(a+bx),\\ a_2(x) = \frac{1}{4}[-2a+a^2+x(4A+ b^2 x)], \quad a_1(x) = \frac {1}{4}( 2A-ab)(-2+a+bx) \end{gather*} with the eigenvalues \[ \lambda_n= \frac{b}{4} n (2 {A}-a b-b+b n), \] which in general are not iterates of second order operator. The special case $a=4,b=-2, A=-5$ gives the classical fourth order Laguerre operator \cite{littlejohn-krall89} with the weight $p(x)=e^{-x}$ as \[ a_4(x) = x^2, \quad a_3(x) = -2(-2+x)x, \quad a_2(x) = x^2-5 x+2, \quad a_1(x) = -1+x \] and the eigenvalues $\lambda_n= n^2$. This coincides with the second iteration of the classical second order Laguerre operator % \subsection*{Case: $a_4(x)$ has more that one real root} By Corollary \ref{coro3.5}, $a_4(x)$ must have exactly two real roots with multiplicity $2$ and the multiplicity of each real root of $a_4$ in $a_3$ is $1$. By a linear change of variables and scaling, we may assume that the roots are $-1$ and $1$, and take $a_4(x)=(1-x^2)^2$. Now $a_3(x) = K(1-x^2)$ is ruled out by Proposition \ref{prop3.2} (ii). So $a_3(x)$ must have the form $(k_1+k_2 x)(1-x^2)$ with $k_2\ne0$. Without loss of generality, we consider the form $a_3(x) = -2(b+(-2+a)x)(-1+x^2)$ which determines the weight function. \[ p(x) = \frac{(1+x)^{\frac{b-a-2}{2}}}{(1-x)^{\frac{b+a+2}{2}}} \quad\text{with $b-a > 0$ and $b+a < 0$}. \] (1) If $b=0$, then $a_3(x)=2 (-2+a) x(1-x^2)$ with $a < 0$. Solving the determining equation \eqref{det2-n-4} subject to the constraints \eqref{boundary-n-4} determines the fourth order self-adjoint operators \eqref{main-op-n-4} and the differential equations \eqref{main-eq-n-4} as \begin{gather*} p(x) = (1-x^2)^{-1-\frac{a}{2}} \qquad(a < 0)\\ a_4(x) = (1-x^2)^2, \quad a_3(x) = -2 (-2+a) x(-1+x^2)\\ a_2(x) = -2a+a^2+A(-1+x^2), \quad a_1(x) = a(2-3a+a^2-A) x \end{gather*} with the eigenvalues \[ \lambda_n = n (-a+n-1) (-a^2-n a+4 a+n^2+{A}-n-2). \] The particular cases that are iterates of corresponding second order operators are given below: \begin{itemize} \item The values $a=-2$, $A=14$ lead to the Legendre operator \cite{littlejohn-krall89} with the weight $p(x)=1$ as \[ a_4(x) = (1-x^2)^2,\quad a_3(x) = -8x (1-x^2),\quad a_2(x) = 14 x^2-6 ,\quad a_1(x) =4x \] with the eigenvalues $\lambda_n = n^2 (n+1)^2$. \item The special case $a=-1$, $A=7$ is the Chebychev operator of first kind \cite{littlejohn-krall89} with the weight $p(x)=\frac{1}{\sqrt{1-x^2}}$ as \[ a_4(x) = (1-x^2)^2,\quad a_3(x) = -6x (1-x^2),\quad a_2(x) = 7 x^2-4,\quad a_1(x) =x \] with the eigenvalues $\lambda_n = n^{4}$. \item The special case $a=-3$, $A=23$ is the Chebychev operator of second kind \cite{littlejohn-krall89} with the weight $p(x)=\sqrt{1-x^2}$ as \[ a_4(x) = (1-x^2)^2,\quad a_3(x) = -10x (1-x^2),\quad a_2(x) = 23 x^2-8 ,\quad a_1(x) =9x \] with the eigenvalues $\lambda_n = n^2 (n+2)^2. $ \end{itemize} (2) If $b\ne0$, then $a_3(x) = -2(b+(-2+a)x)(-1+x^2)$ with $b-a >0$ and $b+a < 0$. So $a 0, \ b+a<0 \text{ and } b\ne0)\\ a_4(x) = (1-x^2)^2\\ a_3(x) = -2(b+(-2+a)x)(-1+x^2)\\ a_2(x) = \frac{b^3+B+2(-1+a)b^2x-Bx^2+b(-2+a+2x^2-3ax^2 +a^2x^2)}{b}\\ a_1(x)= B + \frac{aB x}{b} \end{gather*} with the eigenvalues \[ \lambda_n = \frac{(a-n+1) n ( -b n^2+a b n+b n-a b+{B}) } {b}. \] \subsection{Self-adjoint operators of order $6$} Consider the self-adjoint operators \begin{equation} \label{main-op-n-6} L=a_{6}(x) y^{(6)}+a_5(x) y^{(5)}+a_4(x) y^{(4)} +a_3(x) y'''+ a_2(x) y''+a_1(x) y' \end{equation} with an admissible weight $p(x)=\exp(\frac{1}{3}\int\frac{a_5(x)} {a_{6}(x)}dx)/| a_{6}(x)|$, satisfying the differential equation \begin{equation} \label{main-eq-n-6} L(y)=\lambda y. \end{equation} By Proposition \ref{prop2.1}, the determining equations in this case are \begin{gather} 3(a_{6}p)'= (a_5 p)\label{det1-n-6}\\ 5(a_5p)''- 6(a_4p)'+ 3(a_3p)=0 \label{det2-n-6n}\\ (a_4p)'''-3(a_3p)''+5(a_2p)'- 5(a_1p)=0 \label{det3-n-6n} \end{gather} on $I$, subject to the vanishing of \begin{equation}\label{boundary-n-6} (a_{6} p),\quad(a_5 p), \quad (a_5p)', \quad(a_4p), \quad(a_4p)'- 3(a_3p), \quad (a_4p)''-3(a_3p)'+5(a_2p) \end{equation} on the boundary $\partial I$. Equations \eqref{det2-n-6n} and \eqref{det3-n-6n} are equivalent to \begin{equation} \frac{10}{27}\big(\frac{a_5}{a_{6}}\big) ^3 +\frac{10}{3}\big(\frac{a_5}{a_{6}}\big) \big( \frac{a_5}{a_{6}}\big) '+ \frac{10}{3}\big( \frac{a_5}{a_{6}}\big) '' -4 \big(\frac{a_4}{a_{6}}\big) '-\frac{4}{3} \frac{a_5}{a_{6}} \frac{a_4}{a_{6}}+ 2\frac{a_3}{a_{6}}=0 \end{equation} and \begin{align*} & -81\big( \frac{{a_3} }{{a_{6}} }\big) '' +\frac{{a_4}}{{a_{6}} } \big( \frac{{a_5} }{{a_{6}} }\big) ^3 -9 \frac{{a_3}}{{a_{6}} } \big( \frac{{a_5} }{{a_{6}} }\big) ^2 +45 \frac{{a_2} }{{a_{6}} } \big( \frac{{a_5} }{{a_{6}} }\big) -135 \frac{{a_1}}{{a_{6}} }+135 \big( \frac{{a_2} }{{a_{6}} }\big) '\\ &+27 \big( \frac{{a_4} }{{a_{6}} }\big) ''' +9 \frac{{a_4} }{{a_{6}} } \big( \frac{{a_5} }{{a_{6}} }\big)'' +27 \frac{{a_5} }{{a_{6}} } \big( \frac{{a_4} }{{a_{6}}}\big) '' -54 \big( \frac{{a_5} }{{a_{6}} }\big) \big(\frac{{a_3} }{{a_{6}} }\big) '\\ & +9 \big( \frac{{a_5} }{{a_{6}}}\big) ^2 \big( \frac{{a_4} }{{a_{6}} }\big) ' +9\frac{{a_4} }{{a_{6}} } \frac{{a_5} }{{a_{6}} } \big( \frac{{a_5} }{{a_{6}} }\big) ' -27 \big( \frac{{a_3} }{{a_{6}} }\big) \big( \frac{{a_5} }{{a_{6}} }\big) '+27 \big( \frac{{a_4} }{{a_{6}} }\big)' \big( \frac{{a_5} }{{a_{6}} }\big) '=0. \end{align*} These identities give the congruences \begin{gather} \label{divide-n-6-1} 5 {a_5} ( {a_5} -6 {a_{6}}') ( {a_5} -3 {a_{6}}') \equiv0 \quad\mod{a_{6}}\\ \label{divide-n-6-2}{a_4} ( {a_5} -9 {a_{6}}') ( {a_5} -6 {a_{6}}') ( {a_5} -3 {a_{6}}') \equiv0 \quad \mod{a_{6}}. \end{gather} Before presenting examples, we note that, in case the leading term has a real root, there are many operators that satisfy the determining equations but for which one of the boundary condition fails at one or both end points of interval $I$; in case the leading term has no real roots, the boundary conditions will be satisfied - because of the form of the weight - but one of the determining equations will not be satisfied. Here are some typical examples. \begin{example} \label{examp4.1}\rm The operator in \eqref{main-op-n-6} with \begin{gather*} p(x) = (x-1)^2 (x+1),\\ a_{6} (x) = (x-1)^2 (x+1)^{4}, \quad a_5 (x) = 3 (x-1) (x+1)^3 (9 x-1),\\ a_4(x) = 60 x (x+1)^2 (5 x-3), \quad a_3(x) = 240 ( 7 x^3+6 x^2-2 x-1) ,\\ a_2(x) = 720 x (5 x+3), \quad a_1(x) = 360 (x+1) \end{gather*} satisfies the determining equations \eqref{det1-n-6}, \eqref{det2-n-6n}, \eqref{det3-n-6n} and all the boundary conditions in \eqref{boundary-n-6} except that the last boundary condition fails at the end point $1$ of $I=[-1,1]$. Therefore, these coefficients and weights would not constitute a self-adjoint operator. \end{example} \begin{example} \label{examp4.2} \rm The operator in \eqref{main-op-n-6} with \begin{gather*} p(x) = e^{-m^2 x} x^2 \quad (m\ne0)\\ a_{6} (x) = x^2, \quad a_5 (x) = -3 x (m^2 x-4),\\ a_4(x) = -30 m^2 x+{A} x^2+30, \quad a_3(x) = 5 x^2 m^{6}-2 ({A} x^2+30 ) m^2+8 {A} x,\\ a_2(x) = x ({C}-3 m^{8} x)+{A} (m^{4} x^2 +12),\\ a_1(x) = 18 m^{8} x-60 m^{6}-8 {A} m^{4} x+24 {A} m^2+{C} (3-m^2 x) \end{gather*} satisfies the determining equations \eqref{det1-n-6}, \eqref{det2-n-6n}, \eqref{det3-n-6n} and all the boundary conditions in \eqref{boundary-n-6} except that the last boundary condition fails at the end point $0$ of $I = [0,\infty)$. Again these coefficients and weights would not constitute a self-adjoint operator. \end{example} \begin{example} \label{examp4.3} \rm The operator in \eqref{main-op-n-6} with \begin{gather*} p(x) = \frac{e^{-m^2 x^2}}{x^2+1} \quad(m\ne0)\\ a_{6} (x) = x^2+1, \quad a_5 (x) = -6 m^2 x (x^2+1),\\ a_4(x) = 10 x^{4} m^{4}-10 m^{4}+{A} x^2+{A}, \\ a_3(x) = -4 m^2 (-10 m^{4}+5 m^2+{A}) x (x^2+1),\\ a_2(x) = {C_2} x^2+{C_1} x+{C_0}, \quad a_1(x) = {D_0 }+{D_1} x \end{gather*} satisfies the determining equations \eqref{det1-n-6}, \eqref{det2-n-6n} and all the boundary conditions in \eqref{boundary-n-6} for $I = (-\infty,\infty)$ but fails to satisfy the remaining determining equation \eqref{det3-n-6n}. \end{example} Examples of sixth order self-adjoint operators with the weights of the form $p(x)=e^{-x^2}$, $p(x)=|x|^n e^{m x}$ and $p(x)=\frac{(1+x)^{m} }{(1-x)^n}$ can be found, similar to fourth order case, by solving the determining equations and boundary conditions using Mathematica. However, because of space constraint, the long expressions for operators and eigenvalues are not reproduced here, and are provided in expanded online version of the paper at http://arxiv.org/abs/1409.2523. \subsection{Self-adjoint operators of order $8$} Consider the self-adjoint operator \begin{equation} \label{main-op-n-8} \begin{aligned} L&=a_{8}(x) y^{(8)}+a_{7}(x) y^{(7)}+a_{6}(x) y^{(6)} +a_5(x) y^{(5)}+a_4(x) y^{(4)}\\ &\quad +a_3(x) y'''+ a_2(x) y''+a_1(x) y' \end{aligned} \end{equation} with an admissible weight $p(x)=\exp(\frac{1}{4}\int\frac{a_{7}(x)} {a_{8}(x)}dx)/| a_{8}(x)|$, satisfying the differential equation \begin{equation} \label{main-eq-n-8}L(y)=\lambda y. \end{equation} By Proposition \ref{prop2.1}, the determining equations in this case are \begin{gather} 4(a_{8}p)'= (a_{7} p)\label{det1-n-8}\\ 7(a_{7}p)''- 6(a_{6}p)'+ 2(a_5p)=0 \label{det2-n-8n}\\ (a_{6}p)'''-2(a_5p)''+2(a_4p)'-(a_3p)=0\label{det3-n-8n}\\ (a_5p)^{(4)} - 4 (a_4p)'''+ 9(a_3p)''-14(a_2p)'+ 14 (a_1p)=0 \label{det4-n-8n} \end{gather} on $I$, subject to the vanishing of \begin{equation} \label{boundary-n-8} \begin{gathered} (a_{8} p),\quad (a_{7} p), \quad (a_{7}p)', \quad (a_{6}p), \quad (a_{6}p)', \quad (a_5 p), \\ 5(a_{6}p)''-11(a_5p)' +14(a_4p),\quad 9(a_{6}p)''-17(a_5p)'+14(a_4p), \\ (a_5p)''-4(a_4p)'+9(a_3p), \quad (a_5p)'''- 4(a_4p)''+9(a_3p)'-14(a_2p) \end{gathered} \end{equation} on the boundary $\partial I$. The examples of eighth order self-adjoint operators with the weights of the form $p(x)=e^{-x^2}$, $p(x)=|x|^n e^{m x}$ and $p(x)=\frac{(1+x)^{m}}{(1-x)^n}$ can be found, similar to fourth order case, by solving the determining equations and boundary conditions using Mathematica. Because of space constraint, the long expressions for operators and eigenvalues are not reproduced here, and instead are provided in expanded online version of the paper at http://arxiv.org/abs/1409.2523. \subsection*{Acknowledgements} We thank Professor Mourad Ismail for his valuable comments and suggestions during the preparation of this paper. \begin{thebibliography}{99} \bibitem{azad1} H. Azad, A. Laradji, M. T. Mustafa; Polynomial solutions of differential equations. \emph{Advances in Difference Equations} (2011) 2011:58 (doi:10.1186/1687-1847-2011-58) \bibitem{azad2} H. Azad, A. Laradji, M. T. Mustafa; Polynomial solutions of certain differential equations arising in physics, \emph{Mathematical Methods in The Applied Sciences} (2012) DOI: 10.1002/mma.2711. \bibitem{bavinck95} H. Bavinck; A direct approach to Koekoek's differential equation for generalized Laguerre polynomials, \emph{Acta Math. Hungar.} 66 (1995), 247--253. \bibitem{bavinck98} H. Bavinck; Differential operators having Sobolev type Laguerre polynomials as eigenfunctions, \emph{Proc. Amer. Math. Soc.} 125 (1997), 3561--3567. \bibitem{bavinck2000} H. Bavinck; A note on the Koekoeks' differential equation for generalized Jacobi polynomials, \emph{J. Comput. Appl. Math.} 115 (2000), 87--92. \bibitem{Boc} S. Bochner; \"{U}ber Sturm--Liouvillesche polynomsysteme, \emph{Math. Zeit.} 29 (1929), 730--736. \bibitem{brenke} W. C. Brenke; On polynomial solutions of a class of linear differential equations of the second order, \emph{Bulletin A. M. S.} 36 (1930), 77-84. \bibitem{everitt57} W. N. Everitt; The Sturm-Liouville problem for fourth-order differential equations, \emph{Quart. J. Math. Oxford Ser. (2)} 8 (1957), 146--160. \bibitem{everitt2001} W. N. Everitt, K. H. Kwon, L. L. Littlejohn, R. Wellman; Orthogonal polynomial solutions of linear ordinary differential equations, \emph{J. Comput. Appl. Math.} 133 (2001), 85--109. \bibitem{ince} Edward L. Ince; \emph{Ordinary Differential Equations}, Dover Publications, New York (1956). \bibitem{Ism}M. E. H. Ismail; \emph{Classical and Quantum Orthogonal Polynomials in one Variable}, paperback edition, Cambridge University Press, Cambridge, 2009. \bibitem{koekoek93} R. Koekoek; The search for differential equations for certain sets of orthogonal polynomials, \emph{J. Comput. Appl. Math.} 49 (1993), 111--119. \bibitem{kockock2000} J. Koekoek, R. Koekoek; Differential equations for generalized Jacobi polynomials, \emph{J. Comput. Appl. Math.} 126 (2000), 1--31. \bibitem{koornwinder84} T. H. Koornwinder; Orthogonal polynomials with weight function $(1-x)^{\alpha}(1 +x)^{\beta}+M \delta(x+ 1) +N\delta(x-1)$, \emph{Canad. Math. Bull.} 27 (2) (1984), 205--214. \bibitem{krall81} A. M. Krall; Orthogonal polynomials satisfying fourth order differential equations, \emph{Proc. Roy. Soc. Edinburgh, Sect. A} 87 (1981), 271--288. \bibitem{Kra40} H. L. Krall; Orthogonal polynomials satisfying a certain fourth order differential equation, \emph{Penn State College Bull.} No 6 (1940), 24 pages. \bibitem{littlejohnpams1994} K. H. Kwon, L. L. Littlejohn, J. K. Lee, B. H. Yoo; Characterization of classical type orthogonal polynomials, Proc. Amer. Math. Soc. 120 (1994), 485--493. \bibitem{kwon2001} K. H. Kwon, L. L. Littlejohn, G. J. Yoon; Orthogonal Polynomial Solutions of Spectral Type Differential Equations: Magnus' Conjecture, \emph{J. Approx. Theory} 112 (2001), 189--215. \bibitem{kwon2006} K. H. Kwon, L. L. Littlejohn, G.J. Yoon; Construction of differential operators having Bochner-Krall orthogonal polynomials as eigenfunctions, J. Math. Anal. Appl. 324 (2006), 285--303. \bibitem{lesky62} P. Lesky; Die Charakterisierung der klassischen orthogonalen Polynome durch Sturm--Liouvillesche Differentialgleichungen, \emph{Arch. Rational Mech. Anal.} 10 (1962), 341--352. \bibitem{littlejohn84} L. L. Littlejohn; On the classiffication of differential equations having orthogonal polynomial solutions, \emph{Ann. Mat. Pura Appl.} 138 (1984), 35--53. \bibitem{littlejohn-krall89} L. L. Littlejohn, A. M. Krall; Orthogonal polynomials and higher order singular Sturm-Liouville systems, \emph{Acta Applicandae Mathematicae,} 17 (1989), 99--170. \bibitem{littlejohn-race1990} L.L. Littlejohn, D. Race; Symmetric and symmetrisable differential expressions, \emph{Proc. Lond. Math. Soc}. 60 (1990), 344-364. \bibitem{naimark} M. A Naimark, \emph{Linear differential operators}, F. Ungar Pub. Co (1967). \bibitem{Rou} E. Routh; On some properties of certain solutions of a differential equation of the second order, \emph{Proc. London Math. Soc.} 16 (1884), 245--261. \bibitem{Sze} G. Szeg\"o; \emph{Orthogonal Polynomials}, fourth edition, Amer. Math. Soc., Providence, 1975. \bibitem{Turbiner} A. Turbiner; On polynomial solutions of differential equations, \emph{J. Math. Phys.} 33 (1992), 3989. \end{thebibliography} \end{document}