On complex oscillation theory, quasi-exact solvability and Fredholm Integral Equations
Yik-Man Chiang
Biconfluent Heun equation (BHE) is a confluent case of the general Heun equation which has one more regular singular points than the Gauss hypergeometric equation on the Riemann sphere $\hat{\mathbb{C}}$. Motivated by a Nevanlinna theory (complex oscillation theory) approach, we have established a theory of periodic BHE (PBHE) in parallel with the Lam\'e equation verses the Heun equation, and the Mathieu equation verses the confluent Heun equation. We have established condition that lead to explicit construction of eigen-solutions of PBHE, and their single and double orthogonality, and a related first-order Fredholm-type integral equation for which the corresponding eigen-solutions must satisfy. We have also established a Bessel polynomials analogue at the BHE level which is based on the observation that both the Bessel equation and the BHE have a regular singular point at the origin and an irregular singular point at infinity on the Riemann sphere $\hat{\mathbb{C}}$, and that the former equation has orthogonal polynomial solutions with respect to a complex weight. Finally, we relate our results to an equation considered by Turbiner, Bender and Dunne, etc concerning a quasi-exact solvable Schr\"odinger equation generated by first order operators such that the second order operators possess a finite-dimensional invariant subspace in a Lie algebra of $SL_2(\mathbb{C})$.
1310.5507.tex
%------------------------------------------------------------------------------
% Beginning of journal.tex
%------------------------------------------------------------------------------
%
% AMS-LaTeX 1.2 sample file for journals, based on amsart.cls.
%
% Replace amsart by the documentclass for the target journal, e.g. tran-l.
%
% 5.1 by Edmund 3rd April 2013
% 5.3 by Edmund 14th April 2013
% 5.5 by Edmund 26th July 2013
% 5.15 by Edmund 21st October 2013; Submission version to Transaction to AMS
\documentclass{amsart}
\usepackage{amscd,amssymb,array,amsmath,epsfig,booktabs}
\usepackage{graphicx,color}
\usepackage{ulem}


\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}

\theoremstyle{definition}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{example}[theorem]{Example}
\newtheorem{xca}[theorem]{Exercise}

\theoremstyle{remark}
\newtheorem{remark}[theorem]{Remark}

\numberwithin{equation}{section}

%    Absolute value notation
\newcommand{\abs}[1]{\lvert#1\rvert}

%    Blank box placeholder for figures (to avoid requiring any
%    particular graphics capabilities for printing this document).
\newcommand{\blankbox}[2]{%
  \parbox{\columnwidth}{\centering
%    Set fboxsep to 0 so that the actual size of the box will match the
%    given measurements more closely.
    \setlength{\fboxsep}{0pt}%
    \fbox{\raisebox{0pt}[#2]{\hspace{#1}}}%
  }%
}

\begin{document}
\title[Oscillation theory, Quasi-Exact Solvability, and Fredholm Integral Equations]{On complex oscillation theory, quasi-exact solvability and Fredholm Integral Equations}

%\title[Eigen-solutions of Biconfluent Heun Equation]{Eigen-solutions of Biconfluent Heun Equation related to Bessel polynomials and Quasi-exactly solvable system}

%    Information for first author
\author{Yik-Man Chiang}
%    Address of record for the research reported here
\address{Department of Mathematics, Hong Kong University of Science and
Technology, Clear Water Bay, Kowloon, Hong Kong,\ P. R. China}
%    Current address
%\address{Department of Mathematics, Hong Kong University of Science and Technology, Clear Water Bay, Kowloon, Hong Kong,\ P. R. China}
\email{machiang@ust.hk}
%    \thanks will become a 1st page footnote.
\thanks{This research was supported in part by the Research Grants Council of the Hong Kong Special Administrative Region, China
(GRF600806). The second author was also partially supported by the  National Natural Science Foundation of China under Grant no. 11371251}

%    Information for second author
\author{Guo-fu Yu}
\address{Department of Mathematics,
Shanghai Jiao Tong University, Shanghai, 200240, P. R.\ China}
\email{gfyu@sjtu.edu.cn}
%\thanks{Support information for the second author.}

%    General info
\subjclass{Primary 30D35, 34M05, 33E10; Secondary 33C47, 45B05}

\date{21 October 2013}

%\dedicatory{This paper is dedicated to our authors.}

\keywords{Biconfluent Heun equations, complex weight, Bessel polynomials, Integral equations, quasi-exact solvable systems.}

\begin{abstract}
Biconfluent Heun equation (BHE) is a confluent case of the general Heun equation which has one more regular singular points than the Gauss hypergeometric equation on the Riemann sphere $\hat{\mathbb{C}}$. Motivated by a Nevanlinna theory (complex oscillation theory) approach, we have established a theory of \textit{periodic} BHE (PBHE) in parallel with the Lam\'e equation verses the Heun equation, and the Mathieu equation verses the confluent Heun equation. We have established condition that lead to explicit construction of eigen-solutions of PBHE, and their single and double orthogonality, and a related first-order Fredholm-type integral equation for which the corresponding eigen-solutions must satisfy. We have also established a Bessel polynomials analogue at the BHE level which is based on the observation that both the Bessel equation and the BHE have a regular singular point at the origin and an irregular singular point at infinity on the Riemann sphere $\hat{\mathbb{C}}$, and that the former equation has orthogonal polynomial solutions with respect to a complex weight. Finally, we relate our results to an equation considered by Turbiner, Bender and Dunne, etc concerning a quasi-exact solvable Schr\"odinger equation generated by first order operators such that the second order operators possess a finite-dimensional invariant subspace in a Lie algebra of $SL_2(\mathbb{C})$.
\end{abstract}

\maketitle

%\section*{This is an unnumbered first-level section head}
%This is an example of an unnumbered first-level heading.

%% The correct journal style for \specialsection is all uppercase; a known bug
%% in amsart.cls prevents this, so input must be uppercase until it is fixed.
%\specialsection*{This is a Special Section Head}
%\specialsection*{THIS IS A SPECIAL SECTION HEAD}
%This is an example of a special section head%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\footnote{Here is an example of a footnote. Notice that this footnote text is running on so that it can stand as an example of how a footnote with separate paragraphs should be written.
%\par
%nd here is the beginning of the second paragraph.}%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%\section{This is a numbered first-level section head}
%This is an example of a numbered first-level heading.

%\subsection{This is a numbered second-level section head}
%This is an example of a numbered second-level heading.

%\subsection*{This is an unnumbered second-level section head}
%This is an example of an unnumbered second-level heading.

%\subsubsection{This is a numbered third-level section head}
%This is an example of a numbered third-level heading.

%\subsubsection*{This is an unnumbered third-level section head}
%This is an example of an unnumbered third-level heading.


\section{Introduction}

The \textit{Biconfluent Heun equation} (BHE)
    \begin{equation}\label{E:BHE}
        z\frac{d^2u}{dz^2}+\big(1+\alpha-\beta z-2z^2\big)\frac{du}{dz}+
        \big((\gamma-\alpha-2)z-\frac12[\delta+(1+\alpha)\beta]\big)u=0,
    \end{equation}
has a regular singular point at $z=0$ and an irregular singular point at $z=\infty$. The irregular singular point is a result of coalesces of two finite regular singular points to that at $z=\infty$ on the Riemann sphere from {the four regular singular points of the} general Heun equation (GHE) \cite{Heun1889} which includes the classical Lam\'e equation as a special case \cite{WW27}, \cite{Ince1926}. Thus, the $z=\infty$ has an higher irregularity rank than the confluent hypergeometic equation and hence its special case, the Bessel equation. Therefore, the BHE plays analogues roles as the confluent hypergeometric equation and the Bessel equation in physical sciences but at an higher level. Indeed, the BHE were studied in several \textit{disconnected cycles} over  different periods between 1930s-1980s (Dunham \cite{Dunham1932}, Masson \cite{Masson1982}, Turbiner \cite{Turb1988}, Bender and Dunne \cite{BenDunn1996}, Gonz\'alez-L\'opez, A., N. Kamran and P. J. Olver \cite{GKO1993}, \cite{GKO1994}, etc), and in the words of Masson \cite{Masson1982} that the BHE has a ``checked history". However,
the understanding of the equation as well as the others that belong to the Heun-class is still far from satisfactory.
\medskip


Hautot \cite{Hautot1969} shows that the BHE possesses a set of polynomial eigen-solutions $\{P_{m,\, \mu}(\alpha,\, \beta;\, x)\}$ that can be written as a very special linear combination of Hermite polynomials. Urwin \cite{Urwin1975} showed that these polynomials are orthogonal with respect to the real weight $x^\alpha\, e^{-\beta x-x^2}$ supported on the positive real axis  in the $L^2(0,\,+\infty)$, such that for each integer $m\ge 0$
    \begin{equation}
        \label{E:Real-orth}
            \int_0^{+\infty}t^\alpha e^{-\beta t-t^2}P_{m,\, \mu}(\alpha,\, \beta;\, t)\,P_{m,\, \nu}(\alpha,\, \beta;\, t)\, dt=\delta_{\mu,\, \nu}\, h_{m,\,\mu,\, \nu},
    \end{equation}
where $0\le \mu,\, \nu\le m$ holds and the value of $h_m$ is generally not known (see also \cite[p. 206]{Ronv1995} and \cite{DMR1978_II}). It is clear that the orthogonality of these polynomials is different from the conventional orthogonality usually found in the literature.
\medskip

The present paper takes a \textit{value distribution (Nevanlinna) theory} approach to the periodic differential equation
    \begin{equation}\label{E:BHE-periodic}
        f^{\prime\prime}(z)+\left(K_{4}e^{4z}+K_{3} e^{3z}+K_{2}e^{2z}+K_{1} e^{z}+K_0\right)f(z)=0,
    \end{equation}
$K_j\in \mathbb{C},\ j=0,\, 1,\, 2,\, 3,\, 4$. We realize that the (\ref{E:BHE-periodic}) is a periodic analogue of the BHE when we were trying to resolve a \textit{complex oscillation problem} (see below). The coefficients $\{K_j\}$ in (\ref{E:BHE-periodic}) are related to the $\alpha,\, \beta,\,\gamma,\,
\delta$ of (\ref{E:BHE}). We abbreviate the periodic BHE as PBHE. The key observation here is that eigen-solutions of the BHE must correspond to \textit{complex non-oscillatory} solutions (see \S\ref{S:complex-osc}) to the PBHE. Besides, the PBHE should be compared with the classical \textit{Lam\'e equation}
    \begin{equation}
        \label{E:lame-eqn}
            f^{\prime\prime}(z)+[h-n(n+1)\,k^2\textrm{sn}^2z]\,f(z)=0,
    \end{equation}
and the \textit{Mathieu equation}
    \begin{equation}
        \label{E:mathieu-eqn}
        f^{\prime\prime}(z)+(a+k^2\cos^2 z)\, f(z)=0
    \end{equation}
 as \textit{periodic forms} of a special case of the \textit{Heun equation} and also a special case of the \textit{confluent Heun equation} respectively. Historically, the Lam\'e and Mathieu equations which have transcendental meromorphic/entire coefficients were derived and studied separate from their algebraic forms differential equations since the 19th century. However, we emphasis that the singularity structure is fundamentally different from the BHE. For example, the BHE has a regular singular point at the origin and an irregular singular point at $\infty$. However, all points in  $\mathbb{C}$ are \textit{ordinary points} for the PBHE. Analogues situations apply to Lam\'e equation in relation to the Heun equation and Mathieu equation in relation to confluent Heun equation.
 \medskip

The advantage of having established a periodic form of the BHE enables us to develop a theory of \textit{Bessel polynomials at the} BHE \textit{level} in parallel with Krall and Frink's work \cite{KF1949}. Indeed, we have found that when the Hautot orthogonal polynomials are written in the \textit{reversed form}, then these polynomials are orthogonal with respect to a \textit{countably infinite family} of complex measures $\{\rho_n(z)\}$ that are supported on the unit circle. Our approach indicates that the BHE is an Heun level analogue of the Bessel equation, and the reversed polynomials that we have identified thus correspond to the classical orthogonal \textit{Bessel polynomials} which have a complex measure and they are supported on the unit circle $|z|=1$.
\medskip

We then show that the complex non-oscillatory solutions of the PBHE (\ref{E:BHE-periodic}) are written in terms of the reversed Hautot polynomials composed with an exponential function, which can be considered as the classical \textit{Lam\'e polynomials} analogues to the Lam\'e equation (\ref{E:lame-eqn}). Whittaker \cite{Whitt1914} appears to be the first who showed that the Lam\'e polynomials satisfy the following \textit{Fredholm integral equation of the second kind}:
        \begin{equation}
            \label{E:lame-integral}
                y(z)=\lambda\,\int_{-2K}^{2K} P_n(k\, \textrm{sn\,}z\,
                \textrm{sn\,} t)\, y(t)\, dt,
        \end{equation}
where the $P_n(t)$ is the $n-$th Legendre polynomial and the $\lambda$ the corresponding eigen-value, and $4K$ the real period of the Jacobian elliptic function $\textrm{sn\,}z$. Moreover, these Lam\'e polynomials possess both single and double orthogonalities with respect to appropriate weights. We refer the reader to Arscott's classic \cite[Chap. 9]{Arscott1964} and \cite{Ronv1995} for more recent development of the topics. In this connection we derive a Fredholm equation of the second kind for the PBHE which exhibits new phenomenon of having a sequence of \textit{eigen-value pairs} in contrast to the single-sequence eigen-values usually encountered, for example, for both the Lam\'e and Mathieu equations. We further show that our complex non-oscillatory solutions to (\ref{E:BHE-periodic}) written as a \textit{composition} of the Hautot polynomial and exponential function, like the Lam\'e polynomials, which is the composition of Legendre polynomials and a Jacobi elliptic function, are orthogonal on $[-2K,\, 2K]$ and they are also \textit{double orthogonal} over $[-2K,\, 2K]\times [K-2iK^\prime,\, K+2iK^\prime]$ (with respect to a suitable weight), where $K$ is the real period of $\textrm{sn\,}z$,  also have both the \textit{single} and \textit{double orthogonality} for the PBHE over $[0,\, 2\pi i]$ and $[0,\, 2\pi i]\times [\pi,\, \pi+2\pi i]$ respectively.
 \medskip

Finally we consider the following Schr\"odinger differential equation coming from mathematical physics:
\begin{equation}\label{E:Bender-Dunne-eqn-1}
        H\psi=E\psi,
    \end{equation}
where
    \begin{equation}\label{E:Bender-Dunne-eqn-2}
        H=-\frac{d^2}{dx^2}+\frac{(4s-1)(4s-3)}{4x^2}-(4s+4J-2)x^2+x^6.
    \end{equation}
is amongst the equations listed in \cite{Turb1988} in which Turbiner gave a classification of second order differential operators generated by first order operators such that the second order operators possess a \textit{finite-dimensional invariant subspace} in a Lie algebra $g$ of $SL_2(\mathbb{C})$. Such second order differential equations are called \textit{quasi-exactly solvable}. Physicists encountered many quasi-exactly solvable Schr\"odinger equations from a wide range of physical models \cite{Turb1996}.
\medskip

In another study, Bender and Dunne \cite{BenDunn1996} assumes when a power series solution is written in the form
    \begin{equation}
        \label{E:Bender_Dunne_series}
            \psi(x)=e^{-\frac{x^4}{4}}\, x^{2s-1/2}\, \sum_{k=0}^\infty
            \Big(-\frac14\Big)^k\frac{P_k(E)}{k!\,\Gamma(n+2s)}\, x^{2k},
    \end{equation}
then the expansion coefficients satisfy the three-term recursion
    \begin{equation}
        \label{E:BHE-recursion-3}
            P_k(E)=EP_{k-1}(E)+16 (k-1)(k-J-1)(k+2s-2)\, P_{k-2}(E) ,\qquad k\ge 2,
    \end{equation}
where
    \begin{equation}
        \label{E:initial}
            P_0(E)=1,\qquad P_1(E)=E,
    \end{equation}
which possesses a remarkable \textit{factorization property} that when the subscript exceeds the \textit{critical value} $J$, then
    \begin{equation}
        \label{E:factor}
            P_{n+J}(E)=P_J(E)\, Q_n(E),\qquad n\ge 0
    \end{equation}
hold. In addition, Bender and Dunne observe \cite{BenDunn1996} from (\ref{E:BHE-recursion-3}) (see \cite{AAR1999}) that the $\{P_n(E)\}$ forms a family of orthogonal polynomials. These polynomials are now called \textit{Bender-Dunne} polynomials.
\medskip

We identify that the above Turbiner equation (\ref{E:Bender-Dunne-eqn-1})
is essentially a special case of the BHE where the $K_3=0$ in the PBHE. In fact, the Bender-Dunne equation (\ref{E:Bender-Dunne-eqn-1}) is transformed to the equation
    \begin{equation}
        \label{E:Bender-Dunne-BHE}
            z\phi^{\prime\prime}(z)+(-2z^2+2s)\,\phi^\prime(z)+\Big((2J-2)z+\frac{\sqrt{2}}{4}E\Big)\,\phi(z)=0,
    \end{equation}
which is a special case of the BHE (\ref{E:BHE}) via the transformations
    \begin{equation}
        \label{E:BD-trans-1}
        z=\frac{1}{\sqrt{2}}\,x^2,\qquad \psi(x)=\exp\Big(-\frac{1}{4}x^2\Big)x^{2s-1/2}\phi(z).
    \end{equation}
\bigskip

Moreover, we have identified that when the termination of (\ref{E:Bender_Dunne_series}) corresponds to a special case of the Hautot polynomials. So our main focus here is properties of those terminating solutions from (\ref{E:Bender_Dunne_series}). As a results we have verified that the general BHE has quasi-exactly solvable phenomenon, and we have extended many of the results obtained by Bender and Dunne  \cite{BenDunn1996}, such as the expansion (\ref{E:Bender_Dunne_series}) and the three-term recursion (\ref{E:BHE-recursion-3}) including the factorization property (\ref{E:factor}) to those of the BHE. We have also transformed the Turbiner equation (\ref{E:Bender-Dunne-eqn-1}) to its periodic form based on our PBHE. As a result, we have found a Fredholm equation of the second kind, single and double orthogonality for the periodic Turbiner equation as corollaries of what we have established for the PBHE, which maybe of independent interest.
\medskip

 This is the first of a series of papers that studies the fundamentals of the BHE/PBHE from the viewpoint of the complex oscillation theory. The present paper shows how the application of the complex oscillation theory can lead us to (i) identify the reversed Hautot polynomials are a Heun analogue of the classical Bessel polynomials, the orthogonality is supported on the unit circle $|z|=1$ with respect to a complex measure, (ii) the eigen-solutions to the PBHE with respect to having finite \textit{exponent of convergence of zeros} $\lambda(f)<+\infty$ (see (\ref{E:lambda}) in \S2)
can be given in terms of the exponential type Hautot polynomials (it composes with exponential function), thus generating novel \textit{eigen-values pairs} that characterise the corresponding eigen-solutions, (iii) show that these exponential-type Hautot polynomials are eigen-solutions to a new Fredholm type homogeneous integral equation of the second kind, (iv) a single orthogonality of the exponential-type Hautot polynomials, and (v) a double orthogonality for products of the exponential-type Hautot polynomials.
\medskip

This paper is organized as follow. A brief outline of the complex oscillation theory and the connection to (\ref{E:BHE} will be given in \S\ref{S:complex-osc}. \S\ref{E:bessel-poly}  will be devoted to describing basic properties of the classical generalized Bessel polynomials and to show how one can construct Bessel polynomials at the Heun level. We show the reversed Hautot polynomials are orthogonal on the $|z|=1$ with respect to a complex weight and its construction in \S\ref{S:const-I}). The convergence of the complex weight is proved in \S\ref{S:convergence}. The next three sections of the paper will be devoted to study the complex oscillation theory for PBHE against the exponential type Hautot polynomials. So \S\ref{E:lambda} will give results for complex non-oscillatory solutions to the PBHE, followed by deriving a Fredholm integral equation of the second kind in \S\ref{S:integral}. The single and double orthogonality of the exponential type Hautot polynomials and products of any two such polynomials respectively will be given in next section \S\ref{S:orth}. We return to the Turbiner equation (\ref{E:Bender-Dunne-eqn-1}--\ref{E:Bender-Dunne-eqn-2}) studied by Bender and Dunne in \S\ref{S:quasi} where we derive a three-term recursion relation of the expansion coefficients similar to (\ref{E:Bender_Dunne_series}) for the BHE from which we exhibit expansion factorisation of coefficients similar to (\ref{E:factor}), thus showing that the BHE also exhibits the quasi-exact solvable property elaborated by Bender and Dunne. In fact, the BHE is amongst the equations listed in \cite[(3--5)]{Turb1988} to be quasi-exact solvable from Lie algebraic viewpoint. The \S\ref{S:P-Turbiner} rewrites previous results in \S\ref{S:lambda} and \S\ref{S:orth} for PBHE in the special case for Turbiner's equations in the periodic forms.




\bigskip
\section{Complex Oscillation Theory}\label{S:complex-osc}

The theory utilizes classical Nevanlinna's (Picard's) value distribution theory over $\mathbb{C}$ \cite{Hayman1975} to study the value distribution of solutions of the Schr\"odinger-type linear ordinary differential equations of the form
\begin{equation} \label{E:BL-general}
        f^{\prime\prime}(z)+A(z)f(z)=0,
    \end{equation}
where $A(z)$ is a transcendental entire function on $\mathbb{C}$. One can deduce that $z=0$ could be an exceptional value in the sense of Picard and hence Nevanlinna. It turns out that the study of zero distribution of $f$ is of particular importance. More specifically, let $n(r,\, f)$ counts the number of zeros of the solution $f$ of (\ref{E:BL-general}) in a disc of radius $r$. It is fundamental that
        \begin{equation}
            \label{E:lambda}
            \lambda(f):=\varlimsup_{r\to +\infty}\frac{\log n(r,\, f)}{\log r}
\le\varlimsup_{r\to +\infty}\frac{\log T(r,\,f)}{\log r}= \rho(f)=+\infty
        \end{equation}
hold whenever $A(z)$ is transcendental and for any non-trivial solution. Here the $\lambda(f)$ on the left side of (\ref{E:lambda}) is called the \textit{exponent of convergence of the zeros} of the solution $f(z)$, while the $\lambda(f)$ on the right side is called the \textit{Nevanlinna order} of $f$ \cite{Hayman1975,La93}.
Bank and Laine then \cite{BL1983(1)} developed an oscillation theory for the (\ref{E:BL-general}) with periodic coefficient $A(z)=B(e^z)$ where
    \begin{equation}\label{E:periodic-coeff}
        B(\zeta)=\frac{K_s}{\zeta^s}+\cdots+K_0+\cdots+K_\ell \zeta^\ell,
    \end{equation}
and $K_j\, (k=s,\cdots,\, \ell)$ are complex constants. They have shown by using Nevanlinna theory \cite{Hayman1975}, that $\lambda(f-a)=+\infty$ for all finite non-zero $a\in\mathbb{C}$ and for all entire solutions of (\ref{E:BL-general}) for (\ref{E:periodic-coeff}). The Picard theorem asserts that $a=0$ may be the only exceptional value. That is, this solution could have $\lambda(f)<\sigma(f)=+\infty$.  Besides, such an exceptional solution assumes the semi-explicit representation  \cite{BL1983(1)}
    \begin{equation}\label{E:BL-representation}
        f(z)=\psi(e^{z/h})\exp \big(dz+P(e^z)\big)
    \end{equation}
where $h=2$ if $\ell$ is odd and $h=1$ otherwise, $\psi(\zeta)$ is a polynomial, $P(\zeta)$ is a Laurent polynomial and $d$ is a constant.
\medskip

 Previous studies by the one of the authors with others for certain  periodic special cases of (\ref{E:BL-general}) and its non-homogeneous analogues   \cite {CI2006,CYu2011,CYu2012} only involved the use of Bessel functions, Coulomb-wave functions, Lommel functions and related functions. As a results, new value distribution (global) properties for these classical special functions were found. In the present investigation, a new class of orthogonal polynomials at the BHE level are found and this also leads to a new Fredholm-type integral equation for the reversed BHE. The main results obtained here will be used to study the corresponding complex oscillation problems elsewhere.
\medskip


In order to better appreciate the role played by the complex oscillation theory in the study of BHE, let us  review previous results in \cite{CI2006} in which Ismail and one of the authors \cite{CI2006} where able to develop suitable special function tools to identify that the polynomial $\psi(\zeta)$ corresponds to the equation (\ref{E:BL-general}) for (\ref{E:periodic-coeff}) with $(s,\, \ell)=(0,\, 1)$ considered earlier by Bank, Laine and Langley in \cite{BLL1986}  belongs to \textit{reversed Bessel polynomial}. While the $\psi(\zeta)$ in  (\ref{E:BL-representation}) of (Morse-type equation) equation (\ref{E:BL-general}) for (\ref{E:periodic-coeff}) with $(s,\, \ell)=(-2,\, 0)$, that is,
    \begin{equation}
\label{E:CI_2}
f^{\prime\prime}+\left(K_{-2}e^{-2z}+K_{-1} e^{-z}+K_0\right)\,f=0,
    \end{equation}
was considered in \cite{BL1983(1)}  is a \textit{generalized Bessel polynomial}. Indeed, the criterion in terms of the coefficients in (\ref{E:CI_2}) and the special solution representations found in \cite{CI2006} rendered the complex oscillation problem for the case $(s,\, \ell)=(-2,\, 0)$ completely solved \cite{BLL1989}.  However, the complex oscillation problems for the majority of other cases $(s,\, \ell)$ remain unsolved. The first author observed that the equation \cite[(3.3)]{CI2006} contains, as a special case, the periodic form (\ref{E:BHE-periodic}) of BHE \cite[p. 199]{Ronv1995}. More specifically, the BHE (\ref{E:BHE}) corresponds to the  case of $(s,\, \ell)=(0,\, 4)$ in (\ref{E:periodic-coeff}).
According to the theory of Bank and Laine, a solution $f$ of equation (\ref{E:BHE-periodic}) with $\lambda(f)<+\infty$ admits a solution of form (\ref{E:BL-representation}) with $h=1$. This suggests the polynomial component $\psi(\zeta)$ of those solutions $f$ for the equation (\ref{E:BHE-periodic})  with $\lambda(f)<+\infty$ may come from certain special functions with important properties, such as orthogonality, etc. We show this is indeed the case below with a ``special property" replaced by \textit{quasi-exact solvability}  \cite{BenDunn1996}, but with respect to a \textit{complex weight} supported on the unit circle $|z|=1$. In fact, these eigen-solutions of the BHE when \textit{reversed} is a new class of orthogonal polynomials, a phenomenon that is analogues in the classical theory for the (generalized) Bessel polynomials that are related to Coulomb Wave equation and Bessel equation.
%This paper is organized as follow. in order to better appreciate the orthogonality of these new orthogonal polynomials, the \S 2 will be devoted to describe basic properties of the classical generalized Bessel polynomials. \S3 will be ...


%We will explain in the next section how an observation from the complex oscillation theory helps us to find The equation (\ref{E:periodic}) can be transformed, via the Lommel transformation (see \cite{CI2006}), to the following equation



\section{Bessel polynomials at Biconfluent Heun equation Level}\label{E:bessel-poly}
Let us recall that although the \textit{Bessel equation of order $\nu$}
    \begin{equation}\label{E:bessel-eqn}
        x^2f^{\prime\prime}(x)+xf^\prime(x)+(x^2-\nu^2)f(x)=0,
    \end{equation}
is not known to admit any polynomial-type eigen-solution, but the \textit{MacDonald function} $K_\nu(z)$ \cite[p. 80]{Watson:1944} obtained from the equation  (\ref{E:bessel-eqn}) via a rotation by $\frac{\pi}{2}$ induces a polynomial $\theta_n(z)$ of degree $n$ when $\nu=n+\frac12$:
    \begin{equation}
        \label{E:induced_I}      K_{n+\frac12}(z)=\sqrt{\frac{\pi}{2}}z^{-n-\frac12}e^{-z}\theta_n(z).
        \end{equation}
The degree $n$ polynomial $\theta_n(z)$, called the \textit{reverse Bessel polynomial of degree $n$} \cite{Gr1978}, does satisfy a differential equation
    \begin{equation}
        \label{E:rev-bessel-poly-eqn}
            z\, \theta^{\prime\prime}(z)-2(z+n)\,\theta^{\prime}(z)+2n\,\theta(z)=0.
    \end{equation}
This equation has a regular singular point at $z=0$ and an irregular singular point at $z=\infty$ as is the case for the Bessel equation (\ref{E:bessel-eqn}) itself. Reversing this equation by the transformation $y_n(z)=z^n\theta_n(\frac1z)$ leads to the equation
    \begin{equation}\label{E:bessel-poly-eqn}
        z^2\frac{d^2y}{dz^2}+(2z+2)\frac{dy}{dz}-n(n+1)y=0,
    \end{equation}
where the polynomial solutions $y_n(z)$, called \textit{Bessel polynomial of degree} $n$. This equation is called \textit{Bessel polynomial equation} which has an irregular singular point at the origin and a regular singular point at $\infty$. The so-called \textit{generalized Bessel polynomial equation}
    \begin{equation}\label{E:gen-bessel-poly-eqn}
        z^2\frac{d^2y}{dz^2}+(az+b)\frac{dy}{dz}-n(n+a-1)y=0,
    \end{equation}
where $a\not=-1,\, -2,\, -3,\,\cdots$ and $b\not=0$, also shares this property. Krall and Frink \cite{KF1949} called $\{y_n(a,\,b;\, z)\}$ the \textit{generalized Bessel polynomials}, and showed that they are orthogonal on the unit circle $|z|=1$ with respect to a complex weight
    \begin{equation}\label{E:bessel-weight}
        \rho(z)=\frac{1}{2\pi i}\sum_{k=0}^\infty\frac{\Gamma(a)}{\Gamma(a+n-1)}
        \Big(-\frac{b}{z}\Big)^k,
    \end{equation}
which is given in an infinite Laurent series.
Although the Bessel polynomials were formally named by Krall and Frink in 1949 \cite{KF1949}, they appeared in an earlier work of Bochner \cite{Boch1929}  that they along with the Jacobi, Laguerre and Hermite polynomials are the only orthogonal polynomials that belong to the so-called \textit{Sturm-Liouville class polynomials}. They also appear in the work of Burchnall and Chaundy about commuting differential operators \cite{BC31} in 1931 and in \textit{complex oscillation problems} by Bank, \textit{et al} \cite{BL1983(1),BLL1986} as identified by Chiang and Ismail \cite{CI2006}, see also Andrews, Askey and Roy \cite{AAR1999}, Grosswald \cite{Gr1978} and Szeg\H{o} \cite{Sz1975} for their recent development.

It follows from the complex oscillation theory that the equation (\ref{E:BL-general}) with (\ref{E:periodic-coeff}) such that $s=0$ {in (\ref{E:periodic-coeff})} can be transformed, via $x=e^z$, $\Psi(x)=f(z)$, to the equation \cite[(3.2)]{CI2006}:
    \begin{equation}
        \label{E:Generalized-Bessel}
            x^2\Psi^{\prime\prime}(x)+x\Psi^{\prime}(x)+\Big(\sum_{j=0}^\ell K_j x^{j}\Big)\Psi(x)=0,
    \end{equation}
which we call a \textit{generalized Bessel equation}. When $\ell=4$, the function $\psi(x)$ defined by
    \begin{equation}
        \label{E:induced_II}
            \Psi(x):= x^{\frac{\alpha}{2}}\,e^{-\frac{\beta}{2}\,x-\frac{1}{2}x^2}\,\psi(x)
    \end{equation}

satisfies the equation (\ref{E:BHE}) (see \cite[p. 62]{DMR1978_II}). Hence the function $\psi(x)$ plays the role of $\theta(z)$ in (\ref{E:induced_I}). We now show that the reverse polynomial of $\psi(x)$ of degree $n$, that is
    \begin{equation}
        \label{E:rev-BHE-trans}
            Y(x):=x^n\psi(1/x),
    \end{equation}
satisfies the differential equation
    \begin{equation}
        \begin{split}
        \label{E:rev-BHE-poly-eqn}
            x^3\, Y^{\prime\prime}(x) +[(1-2n-\alpha)x^2 &+\beta x+2]\,Y^{\prime}(x)\\
            & +[(\alpha+n)nx-\beta n-\theta)]Y(x)=0,
        \end{split}
    \end{equation}
where
    \begin{equation}
             \theta=\frac12\big[\delta+\beta(1+\alpha)\big]=\frac12\big[\delta+\beta(1+\gamma-2-2n)\big].
    \end{equation}
That is, the $Y_n(x)=Y(x)$ is an BHE analogue of the generalized Bessel polynomial $y_n(a,\,b;\, z)$. Therefore, the equation (\ref{E:rev-BHE-poly-eqn}) is the BHE analogue of the generalized Bessel polynomial equation (\ref{E:gen-bessel-poly-eqn}).

If the BHE (\ref{E:BHE}) \cite[pp. 203--206]{Ronv1995} admits a polynomial solution $P_{n,\, \mu}(\alpha,\, \beta;\, x)$, then it is necessary that $\gamma-\alpha-2=2n$ holds. In fact, a power series (analytic) solution is given by
    \begin{equation}
        \label{E:infinite-sum}
        N(\alpha,\,\beta,\,\gamma,\,\delta;\, x)=\sum_{k=0}^\infty \frac{A_k}{(1+\alpha)_k\, k!}\, x^{k},
    \end{equation}
where $(\alpha)_k=\Gamma(\alpha+k)/\Gamma(\alpha)$, and the coefficients $A_k,\, k\ge0$ satisfy the three-term recursion formula:
    \begin{align}
        \label{E:BHE-recursion}
             A_{k+2}&=\Big\{(k+1)\beta+\frac12\big[\delta+\beta(1+\alpha)\big]\Big\}A_{k+1}\notag\\
                    &\qquad
                    -(k+1)(k+1+\alpha)(\gamma-\alpha-2-2k)A_k,\quad k\ge 0
    \end{align}
where
    \begin{equation}
        \label{E:recursion-initial}
            A_0=1,\qquad A_1=\frac12\Big[\delta+\beta(1+\alpha)\Big].
    \end{equation}
Clearly the three-term recursion (\ref{E:BHE-recursion}) terminates if and only if
    \begin{equation}
        \label{E:terminates}
            \gamma-\alpha-2=2m,\qquad A_{m+1}=0,
    \end{equation}
simultaneously, where $m$ is some non-negative integer. It is also easy to see from induction that $A_{m+1}(\delta)$ is a polynomial of $\delta$ of degree $m+1$, hence possessing at most $m+1$ roots $\delta_\nu^m,\ \nu=0,\, 1,\, \cdots m$ at degree $m$. When the series solution terminates, then we write
    \begin{equation}
        \label{E:eigen-soln}
        P_{m, \, \nu}(\alpha,\,\beta;\,x)=N\big(\alpha,\,\beta,\,\alpha+2(m+1),\,\delta^m_{\nu};\, x\big)\quad 0\le \nu\le m,
    \end{equation}
for $m=0,\, 1,\, 2,\, \cdots$. We note that when $\alpha+1>0$ and $\beta\in \mathbb{R}$, then the $m+1$ roots are real. When all these roots are simple, then the polynomial solutions described are precisely those orthogonal polynomials described in (\ref{E:Real-orth}) (see \cite[\S 3.3]{Ronv1995}).

\section{Construction of a complex weight I}\label{S:const-I}
     Since $P_{n,\, \mu}(\alpha,\, \beta;\, x)$ exhibits quasi-exact solvability,  so is its reversed polynomial
    \begin{equation}
            Y_{n,\, \mu}(x)=~x^nP_{n,\, \mu}(\alpha,\, \beta;\, 1/x).
    \end{equation}
We now show that these polynomials $\{Y_{n,\, \mu}(x)\}$  are orthogonal with respect to a complex weight $w(z)$ supported on the unit circle $|z|=1$, in much the same way as that the generalized Bessel polynomials do.

\begin{theorem}\label{T:orthogonality} Let $\rho(z)=\rho_n(z)=\sum_{k=0}^{\infty}a_kz^{-k}$ be a formal series such that the coefficients $a_k$ satisfies
    \begin{align}
            &(1+2n+\alpha)\,a_1-\beta a_0=0, \label{E:recu1}\\
            &(2-k+2n+\alpha)\,a_k-\beta a_{k-1}-2a_{k-2}=0,\qquad k=2,\,3,\, 4,\, \cdots   \label{E:recu2}
    \end{align}
Then for each integer $n\ge 0$, $\{Y_{n,\, \mu}(x)\}$ satisfies,  the orthogonality
    \begin{equation}
        \label{E:complex-orth}
            \int_{\mathcal{C}} Y_{n,\, \mu}(t)\, Y_{n,\, \nu}(t) \, \rho_n(t) \, dt
            =0,
    \end{equation}
whenever, ${\mu\not= \nu}$ for $0\le \mu,\, \nu\le n$, where $\mathcal{C}$ is any simple closed curve that surrounds the origin.  In particular, if $\beta=0$ in \textrm{(\ref{E:recu1}-\ref{E:recu2})} then we have the weight function $\rho_n(x)$ in closed form
    \begin{equation}
        \label{E:BHE-complex-weight}
             \rho_n(z)=\sum_{k=0}^{\infty}\frac{\Gamma(1-n-\frac{\alpha}{2})}{\Gamma(k-n-\frac{\alpha}{2})}\left(
-\frac{1}{z^2}\right)^k.
    \end{equation}
\end{theorem}

\begin{proof} We first derive the weight function $\rho_n$ (\ref{E:BHE-complex-weight}) and leave its convergence to the next section where more detailed analysis of the series is needed. We adapt a well-known procedure,  see for example the well-known \cite{KF1949}, to find a complex weight for the reverse BHE (\ref{E:rev-BHE-poly-eqn}). The first step is to construct a factor $\sigma(z)$ that makes the (\ref{E:rev-BHE-poly-eqn}) self-adjoint. That is, $\sigma(z)$ satisfies the first order differential equation
    \begin{equation}
        \label{E:self-adjoint}
            \big(z^3\sigma(z)\big)^\prime=[(1-2n-\alpha)z^2+\beta z+2]\,\sigma(z).
    \end{equation}

Hence
    \begin{equation}
        \label{E:self-adjoint-soln}
            \sigma(z)=z^{-(2+2n+\alpha)}e^{-\frac{\beta}{z}-\frac{1}{z^2}}.
    \end{equation}
However, this weight is branched unless $2+2n+\alpha$ is an integer. This suggests that we aim for a weight, depending on $n$, in the formal series
    \begin{equation}
        \label{E:weight}
            \rho(z)=\rho_n(z)=\sum_{k=0}^{\infty}a_kz^{-k},
    \end{equation}
that satisfies the differential equation
    \begin{equation}
        \label{E:self-adjoint-II}
        \big(z^3\rho(z)\big)^\prime-[(1-2n-\alpha)z^2+\beta
z+2]\, \rho(z)=(2+2n+\alpha)\big(-n-\frac{\alpha}{2}\big)\,z^2
    \end{equation}
instead. Without loss of generality, we may assume $a_0=-n-\frac{\alpha}{2}$.
Then one can easily verify that
    \begin{equation}
        \begin{split}
        (z^3\rho(z))^\prime &-[(1-2n-\alpha)z^2+\beta z+2]\rho(z)-(2+2n+\alpha)\big(-n-\frac{\alpha}{2}\big)z^2\\
        &=[(2+2n+\alpha)z^2-\beta z-2]\,\rho(z)+z^3\rho^\prime(z)-(2+2n+\alpha)\big(-n-\frac{\alpha}{2}\big)z^2\\
        &=\sum_{k=0}^{\infty}(2-k+2n+\alpha)\,a_kz^{-k+2}-\beta
         \sum_{k=0}^{\infty}a_kz^{-k+1}-2\sum_{k=0}^{\infty}\,a_k z^{-k}\\
         &\qquad  -(2+2n+\alpha)\big(-n-\frac{\alpha}{2}\big)z^2\\
        &=
        %(2+2m+\alpha)a_0z^2+
        [(1+2n+\alpha)a_1-\beta a_0]z\nonumber\\
        &\quad+\sum_{k=2}^{\infty}[(2-k+2n+\alpha)\,a_k-\beta a_{k-1}-2\,a_{k-2}]z^{-k+2},
        \end{split}
    \end{equation}
the coefficients $a_k$ satisfy the following three-term recurrence relation:
\begin{align*}
            &(1+2n+\alpha)\,a_1-\beta a_0=0,\\
            &(2-k+2n+\alpha)\,a_k-\beta a_{k-1}-2a_{k-2}=0,\qquad k=2,\,3,\, 4,\, \cdots
    \end{align*}
This gives the weight function as asserted in the Theorem when  $\beta \neq 0$. Unfortunately, as in the cases of similar consideration related to Heun equation, at least three-term recurrence will be obtained, thus making explicit expression difficult to find. The convergence of the (\ref{E:weight}) will be established in \S\ref{S:convergence} where a detailed account on series defined by three-term recurrence will be presented. In the case when $\beta=0$, we easily see from the relations (\ref{E:recu1}--\ref{E:recu2}) that the odd terms $a_{2k+1}=0$. We set $a_{2k}=(-1)^kb_k$,
then from the relation (\ref{E:recu2}) we have
    \begin{equation}
        b_k=\frac{1}{k-1-(n+\alpha/2)}\, b_{k-1},\quad k=1,\,2,\, 3,\,\cdots.
    \end{equation}
So we have
    \begin{equation}
        b_k=\frac{1}{(-n-\alpha/2)_{k}}\, b_0=\frac{\Gamma(1-n-\frac{\alpha}{2})}{\Gamma(k-n-\frac{\alpha}{2})}.
    \end{equation}
Then the complex weight (\ref{E:weight}) for the (\ref{E:rev-BHE-poly-eqn}) could explicitly be written as
    \begin{equation}
        \label{E:special-weight}
        \rho_n(z)=\sum_{k=0}^{\infty}\frac{\Gamma(1-n-\frac{\alpha}{2})}{    \Gamma(k-n-\frac{\alpha}{2})}\left(
-\frac{1}{z^2}\right)^k,
    \end{equation}
which converges trivially at any $z\not=0$.
\medskip

We next show the orthogonality. We suppose $Y_{n,\nu}$ and $Y_{n,\mu}$ are polynomial solutions of (\ref{E:rev-BHE-poly-eqn}) with $\theta=\mu$ and
$\theta=\nu$ respectively. That is,

    \begin{align}
        \label{E:poly-nu}
            x^3Y_{n,\nu}^{\prime\prime}(x) +[(1-2n &-\alpha)x^2 +\beta x+2]\, Y_{n,\nu}^\prime(x)\nonumber\\
            & +[(\alpha+n)nz-\beta n-\mu]\, Y_{n,\nu}(z)=0
    \end{align}
Multiplying (\ref{E:poly-nu}) throughout by the weight $\rho_n(x)$ and applying the self-adjoint differential equation (\ref{E:self-adjoint-II}) to the resulting equation yield
    \begin{align}\label{mp}
        \big(x^3\rho_n(x)\, Y^\prime_{n,\,\nu}\big)^\prime(x) +(2+2n &+\alpha) \left(n+\frac{\alpha}{2}\right)x^2 Y_{n,\nu}^\prime(x)\nonumber\\
        &+[(\alpha+n)nx -\beta n-\nu]\,\rho_n(x)\, Y_{n,\,\nu}(x)=0,
    \end{align}
Multiplying the (\ref{mp}) throughout by $Y_{n,\mu}(x)$ and integrating the resulting equation along the unit circle $U:=\{x:\ |x|=1\}$ yield
    \begin{align}
        \label{E:first-integration}
            \int_U \big(t^3\rho_n &Y^\prime_{n,\nu}(t)\big)^\prime\,Y_{n,\mu}(t)\, dt
            +(2+2n+\alpha)\big(n+\frac{\alpha}{2}\big)\int_U t^2\,Y^\prime_{n,\, \nu}(t)\,              Y_{m,\,\nu}(t)\,dt\nonumber\\
            &\qquad+ \int_U [(\alpha+n)n t-\beta n-\nu]\, \rho_n(t)\, Y_{n,\,\nu}(t)\,Y_{n,\,\mu}(t)\,dt\nonumber\\
            &=0-\int_U t^3\, Y_{n,\,\nu}^\prime(t)\, Y^\prime_{n,\,\mu}(t)\,\rho_n(t)\, dt+0 \nonumber\\
            &\qquad + \int_U [(\alpha+n)n t-\beta n-\nu]\,\rho_n(t)\, Y_{n,\,\nu}(t)\, Y_{n,\,\mu}(t)\, dt\nonumber\\
            &=\int_U t^3 Y^\prime_{n,\,\nu}(t)\,Y^\prime_{n,\, \mu}(t)\,\rho_n(t)\, dt+\int_U [(\alpha+n)n t-\beta n-\nu]\, \rho_n(t)\, Y_{n,\nu}(t)\,Y_{n,\mu}(t)\,dt.
    \end{align}
Interchanging the roles of $\nu$ and $\mu$ yields a similar equation as (\ref{E:first-integration}). Then subtracting this equation with (\ref{E:first-integration}) implies
    \[
        \int_U Y_{n,\nu}(t)\,Y_{n,\mu}(t)\, \rho_n(t)\,dt=0
    \]
whenever $\nu\not= \mu$. This proves the orthogonality relation (\ref{E:complex-orth}).
\end{proof}
\bigskip

\section{Construction of a complex weight II: \\ completion of the proof of Theorem \ref{T:orthogonality}}\label{S:convergence}

We now deal with the convergence of the formal series for the weight function $\rho_n(x)$ defined in (\ref{E:recu1}-\ref{E:recu2}):
    \begin{align}
            &(1+2n+\alpha)\,a_1-\beta a_0=0, \label{E:recu3}\\
            &(2-k+2n+\alpha)\,a_k-\beta a_{k-1}-2a_{k-2}=0,\qquad k=2,\,3,\, 4,\, \cdots   \label{E:recu4}
    \end{align}
\medskip

\begin{theorem} For each integer $n$, the weight function $\rho_n(x)$ in Theorem \ref{T:orthogonality}, defined by the formal power series (\ref{T:orthogonality}) where its coefficients satisfy the three-term recursion (\ref{E:recu3}-\ref{E:recu4}), converges everywhere in $\mathbb{C}$ except when $z=0$.
\end{theorem}
\bigskip

%we consider the limit of the ratio $\frac{a_{n+1}}{a_n}$, in order to do so we set $b_n=\frac{a_{n+1}}{a_n}$, then we have
%\begin{eqnarray}
%b_n=\frac{\beta}{1-n+2m+\alpha}+\frac{2}{1-n+2m+\alpha}\frac{1}{b_{n-1}},
%\end{eqnarray}
%We know that if
%\begin{eqnarray}
%\lim  b_n=0,
%\end{eqnarray}
%then the power series (\ref{ps}) would  converge in the whole plane except zero. Through numerical computation we find that the limit $b_n$ does not exist in many cases when we fix initial value $b_0$ and parameters $m, \alpha,\beta$. But we real find that when the follow relation holds,
%\begin{eqnarray}
%&& b_0=\frac{1}{2.5},\qquad
%b_{n}=\frac{1}{2.5-n}+\frac{2}{2.5-n}\frac{1}{b_{n-1}},\quad
%k=1,2,\cdots,
%\end{eqnarray}
%Then the limit exists and $b_n \rightarrow 0\quad (n\rightarrow +\infty).$

In order to establish this claim, we utilize classical asymptotic theory results of linear second order difference equations established by Poincar\'e, later refined by Perron (1911) and Kreuser (1914). A detailed survey can be found in \cite[\S 2]{Gautschi1967} or in Wimp \cite{Wimp1984}.  For the sake of completeness, we state their results as stated in \cite{Gautschi1967} but tailored to our application below.

\begin{theorem}[see \cite{Gautschi1967}] Consider the second order difference equation
    \begin{equation}
        \label{E:Perron}
            y_{n+1}+A_n\, y_n+B_n\, y_{n-1}=0,\qquad n=1,\,2,\, 3,\, \cdots,
    \end{equation}
where $B_n\not= 0$ for all $n$. Suppose
    \begin{equation}
        \label{E:asymptotic}
            A_n\sim A\,n^\alpha,\qquad B_n\sim B\,n^\beta,\qquad AB\not=0,\quad
            \alpha,\, \beta \ \textrm{real};\qquad n\to \infty.
    \end{equation}
We then construct the polygonal path $P_0P_1P_2$ in a Newton-Puiseux diagram where the points $P_0,\, P_1,\, P_2$ have coordinates $(0,\, 0),\, (1,\, \alpha),\, (2,\, \beta)$ respectively. We denote the gradients of the straight line segments of $P_0P_1,\, P_1P_2$ by $\sigma=\alpha,\, \tau=\beta-\alpha$ respectively. Then the followings hold
    \begin{enumerate}
        \item If the point $P_1$ lies above the line segment $P_0P_2$ (i.e., $\sigma>\tau$), then the difference equation (\ref{E:Perron}) has two linearly independent solutions $y_{n,\, 1}$ and $y_{n,\, 2}$, such that
        \begin{equation*}
            \frac{y_{n+1,\, 1}}{y_{n,\, 1}}\sim -An^\sigma,\qquad
            \frac{y_{n+1,\, 2}}{y_{n,\, 2}}\sim -\frac{B}{A}n^\tau,\qquad n\to\infty.
        \end{equation*}
        \item If the three points $P_j,\, j=0,\, 1,\, 2$ are collinear (i.e., $\sigma=\tau=\alpha$), let $t_1,\, t_2$ be the two roots of $t^2+At+B=0$, and $|t_1|\ge |t_2|$. Then (\ref{E:Perron}) has two linearly independent solutions $y_{n,\, 1}$ and $y_{n,\, 2}$, such that
        \begin{equation*}
            \frac{y_{n+1,\, 1}}{y_{n,\, 1}}\sim t_1n^\alpha,\qquad
            \frac{y_{n+1,\, 2}}{y_{n,\, 2}}\sim t_2n^\alpha,\qquad n\to\infty,
        \end{equation*}
        if $|t_1|>|t_2|$. While if $|t_1|=|t_2|$, then
            \[
                 \limsup_{n\to\infty}\Big(\frac{|y_n|}{(n!)^\alpha}\Big)^{\frac1n}=|t_1|
            \]
        for all non-trivial solutions of (\ref{E:Perron}).
        \item If  $P_1$ lies below the line segment $P_0P_2$, then
            \begin{equation}
                \label{E:Perron-limit}
                 \limsup_{n\to\infty}\Big(\frac{|y_n|}{(n!)^{\frac{\beta}{2}}}\Big)^{\frac1n}=\sqrt{|B|}
            \end{equation}
        for all non-trivial solutions of (\ref{E:Perron}).
        \end{enumerate}
\end{theorem}
\bigskip

\begin{proof}
Let us apply part (3) of the above theorem to show that the formal series (\ref{E:weight}) does converge. To do so let us first rewrite the three-term recurrence relation (\ref{E:recu1}-\ref{E:recu2}) into the form
    \begin{equation}
        \label{E:recur5}
            y_{k+1}+\left(\frac{\kappa}{k-1-2n-\eta}\right)
            y_{k}+\left(\frac{2}{k-1-2n-\eta}\right)y_{k-1}=0,\qquad
k=2,\,3,\, 4,\, \cdots
    \end{equation}
where $y_k=a_k$ defined above.

    It is easily seen from (\ref{E:recur5}) that we have $\alpha=-1,\, \beta=-1$ and $\sigma=-1<0=\tau$ in (\ref{E:asymptotic}). Hence the coordinates $P_1=(1,\, -1)$ ,\, $P_2=(2,\, -1)$ implying that the point $P_1$ lies below the line segment $P_0P_2$. Thus the part (3) of the above theorem applies. Thus, we deduce from (\ref{E:Perron-limit}) that
    \[
        \limsup_{k\to\infty}\left(|y_k|(k!)^{\frac{1}{2}}\right))^{\frac1k}
        =\sqrt{|B|}=\sqrt{2}.
    \]


Thus given $\varepsilon>0$ we can find an integer $k_0$ depending only on $\varepsilon$ such that
    \begin{equation}
        \label{E:convergence-1}
            |y_k|<
%               \frac{(1+\varepsilon)^k}{\sqrt{(k!)}}<
                \frac{(\sqrt{2}+\varepsilon)^k}{\sqrt{(k!)}}
%           < \left(\frac{\sqrt{2}+\varepsilon}{\sqrt{(k^2)}}\right)^k
%           =\left(\frac{\sqrt{2}+\varepsilon}{k}\right)^k
,\qquad k\ge k_0.
    \end{equation}
We deduce from Stirling's formula that when  $k>k_0$
    \[
        \frac12\,\sqrt{2\pi k}\, \left(\frac{k}{e}\right)^k<\ k!\ <2\sqrt{2\pi k}\, \left(\frac{k}{e}\right)^k.
    \]
We deduce for all finite $X$
    \begin{align*}
        \sum_{k=k_0}^\infty \frac{(2\sqrt{2})^k}{\sqrt{(k!)}}|X|^k
        &<\sum_{k=k_0}^\infty
        \frac{(2\sqrt{2})^k}{\left[\frac12\,\sqrt{2\pi k}\, \left(\frac{k}{e}\right)^k\right]^\frac12} |X|^k\\
        &=\sum_{k=k_0}^\infty \frac{\sqrt{2}}{(2\pi)^\frac14\, k^\frac14}
        \left(\sqrt{\frac{8\, e}{k}}\right)^k |X|^k.
    \end{align*}
It follows easily from the ratio test that the series
    \[
            \sum_{k=k_0}^\infty \left(\sqrt{\frac{8\, e}{k}}\right)^k |X|^k
    \]
converges uniformly for all finite $X$:
    \begin{align*}
        \frac{\left[8e/(k+1)\right]^{(k+1)/2}}{(8e/k)^{k/2}}|X|
        &=\left(\frac{8e\, k^k}{(k+1)^{k+1}}\right)^\frac12 |X|\\
         &=\left(\frac{8e}{k\left(1+\frac{1}{k}\right)^{k+1}}\right)^\frac12|X|\to 0
    \end{align*}
as $k\to \infty$ for all $X\not=0$. Hence the original weight function must converge for all $x\not=0$. This completes the proof.

%so it can be deduced that
%\begin{eqnarray}
%\lim\sup_{n\rightarrow \infty}\left(|y_n| \right)^{1/n}\leq \sqrt{2}
%\end{eqnarray}
%Then we have the conclusion that the series (\ref{ps}) would be
%convergent if $|z|> \sqrt{2}$.
\end{proof}

\section{Bessel polynomial-type solutions for periodic BHE (PBHE)}\label{S:lambda}

Without loss of generality, we may assume that $K_4=-1$ for convenience of later calculation. That is, we consider the following equation
\medskip

\begin{theorem}
    \label{T:determinant_soln}
        Let $K_3,\, K_2,\, K_1$ and $K_0$ be complex numbers. Then the equation
\begin{equation}\label{E:BHE-periodic-normal}
        f^{\prime\prime}+\left(-e^{4z}+K_{3}\, e^{3z}+K_{2}\,e^{2z}+K_{1}\, e^{z}+K_0\right)f=0,
    \end{equation}
admits an entire solution with $\lambda(f)<+\infty$ if there exists an non-negative integer $n$ such that the following equation
    \begin{equation}
        \label{E:condition_I}
            \frac{K_3^2}{4}+K_2\pm 2\sqrt{-K_0}=2(n+1),
    \end{equation}
holds amongst the $K_3,\, K_2$ and $K_0$. Moreover, when (\ref{E:condition_I}) holds, then there are
        \begin{enumerate}
            \item $n+1$\footnote{The $n+1$ roots come from choosing the ``$+$" sign from (\ref{E:condition_I}).}
, possibly repeated, choices, of $K_1$
            \item precisely $n+1$ distinct real roots if, in addition that when $K_3,\, K_2,\, K_0<0$ are real and
                \begin{equation}
                    \label{E:distinct}
                        1\pm 2\sqrt{-K_0} >0\footnote{Same as the last footnote.}
                \end{equation}
                hold
        \end{enumerate}
which consist of the roots of the determinant $D_{n+1}(K_1)=0$ where $D_{n+1}(K_1)$ equals to
%   \begin{align}
    \begin{equation}
        \label{E:Det}
    \left|
    \begin{array}{lllll}
        k_1 & -1 & & &  \\
        -(1+2\sqrt{-K_0})k_2 & k_1-K_3 & -1& &  \\
        &
            \begin{split}
                -2(2+2\sqrt{-K_0})\\
                \quad \times (k_2-2)
            \end{split}
        &k_1-2K_3 & -1 & \\
        &\ddots & \ddots & \ddots & \\
        &\ddots & \ddots & \ddots &\\
        & &
            \begin{split}
                    &-(n-1)(n-1+2\sqrt {-K_0}) \\
                    &\quad \times (k_2-2(n-2))
            \end{split}
        &  k_1-(n-1)K_3 &  -1 \\
        & & &
            \begin{split}
                &-n(n+2\sqrt{-K_0})\\
                &\quad \times (k_2-2(n-1))
            \end{split}
        & k_1-nK_3
    \end{array}
    \right|,
    \end{equation}
%   \end{align}
and the $k_1$, $k_2$ are given by
\begin{align}
& k_1=\frac{1}{2}\left(-2K_1-K_3(1+2\sqrt{-K_0})\right),\\
& k_2=K_2+\frac{K_3^2}{4}-2-2\sqrt{-K_0}.
\end{align}
Moreover, we have
	\begin{equation}
        \label{E:BHE-periodic-soln}
        \begin{split}
            BH_{n,\,\nu}(z)=f(z)&= P_{n,\, \nu}\big(2\sqrt{-K_0},\, -K_3,\, 2(n-1+\sqrt{-K_0}),\, -2K_1;\, e^{z}\big)\\
                &\qquad\qquad\qquad \cdot\exp\big[-\frac{K_3}{2}\,e^z-\frac{1}{2}e^{2z}+\sqrt{-K_0}\,z\big]\\
            &= Y_{n,\,\nu}(e^{-z})\, \exp\big[-\frac{K_3}{2}\,e^z-\frac{1}{2}e^{2z}+(n+\sqrt{-K_0})\,z\big]
        \end{split}
    \end{equation}
where $P_{n,\, \nu}$ are the orthogonal polynomials mentioned in (\ref{E:Real-orth}) and their reversed forms, namely $Y_{n\,\nu}(x)=x^n\, P_{n,\, \nu}(\alpha,\, \beta;\, 1/x)$ are {the Bessel orthogonal polynomials of BHE class,} with respect to the complex weight defined in Theorem \ref{T:orthogonality}.
\end{theorem}
\bigskip


\begin{proof}  Suppose that the equation (\ref{E:BHE-periodic-normal}) admits a solution $f$ with $\lambda(f)<+\infty$. Then according to \cite[Thm. 1]{BL1983(1)} (see also \cite[Prop. 1]{Chiang2000}) that the $f$ must assume the form (\ref{E:BL-representation}) where $h=1$ and $\ell=4$ in (\ref{E:periodic-coeff}) is even. Moreover, $s=0$ (\cite[Prop. 1]{Chiang2000}). The transformations $x=e^z$ and $\Psi(x)=f(z)$ transform the (\ref{E:BHE-periodic-normal}) into a generalised Bessel equation (\ref{E:Generalized-Bessel})
    \begin{equation}
        \label{E:Generalized-Bessel-4}
            x^2\Psi^{\prime\prime}(x)+x\Psi^{\prime}(x)+\Big(K_4 x^{4}+K_3x^3+K_2x^2+K_1x+K_0\Big)\Psi(x)=0,
    \end{equation}
when $\ell=4$ in (\ref{E:Generalized-Bessel}).  On the other hand, the transformation (\ref{E:induced_II}) transforms the (\ref{E:Generalized-Bessel-4}) into BHE (\ref{E:BHE}) via the identifications
    \begin{equation}
        \label{E:two-set-parameters-II}
        \begin{array}{lll}
            {K_4= -1}; & K_3=-\beta, & K_2=\gamma-\displaystyle\frac{\beta^2}{4},\\
            K_1=-\displaystyle\frac{\delta}{2}, & K_0=-\displaystyle\frac{\alpha^2}{4}.&\\
        \end{array}
    \end{equation}
between the two sets of parameters of the two equations respectively.
\medskip

We next show that the function defined by (\ref{E:BHE-periodic-soln}) satisfies $\lambda(f)<+\infty$ indeed. We recall from the infinite series solution (\ref{E:infinite-sum}) of the BHE with the coefficients $A_k$ defined by the three-term recursion (\ref{E:BHE-recursion}) terminates if and only if the condition (\ref{E:terminates}) is met. That is, if and only if
    \[
            \frac{K_3^2}{4}+K_2\pm 2\sqrt{-K_0}=2(n+1),
    \]
 and $A_{n+1}=0$ hold in view of (\ref{E:two-set-parameters-II}). However, the condition that $A_{n+1}=0$ can be written as the determinant $D_{n+1}(K_1)=0$ of order $n+1$. In fact, it is a polynomial equation in the variable $K_1$ of degree $n+1$, as asserted in (\ref{E:Det}). It remains to derive the solution (\ref{E:BHE-periodic-soln}). But this follows from the transformation (\ref{E:induced_II}) that the equation (\ref{E:Generalized-Bessel-4}) also has a polynomial solution. In fact, one can compute the coefficient $d$ and the coefficients of the polynomial $P(z)$ from (\ref{E:BL-representation}) by substituting the (\ref{E:BL-representation}) into (\ref{E:Generalized-Bessel-4}). 
Finally, the transformation (\ref{E:rev-BHE-trans}) and $Y_{n\,\nu}(x)=x^n\, P_{n,\, \nu}(\alpha,\, \beta;\, 1/x)$ confirm that the (\ref{E:BHE-periodic-soln}) has $\lambda(f)<+\infty$.
\medskip

We note that if the $K_3,\, K_2,\, K_0$ are arbitrary complex constants that satisfy the (\ref{E:condition_I}), then the $K_1$ could be a repeated root for the $D_{n+1}(K_1)=0$ in case (1). If, however, we have assumed that the $K_3,\, K_2,\, K_0<0$ are real and (\ref{E:distinct}) holds, then (\ref{E:two-set-parameters-II}) implies that $ 1+\alpha=1\pm 2\sqrt{K_0}>0$. We deduce from Rovder's result \cite{Rovder1974} on BHE that the corresponding polynomial solution has distinct real roots. That is, the determinant (\ref{E:Det}) has $n+1$ distinct real roots for $K_1$.
\end{proof}
\bigskip

\begin{remark} We observe from part (2) of Theorem \ref{T:determinant_soln} that for each integer $n$  where the relation (\ref{E:condition_I}) holds, there are $n+1$ choices of $K_1$ where the equation (\ref{E:BHE-periodic-normal}) admits a solution of form (\ref{E:BHE-periodic-soln}), so when we describe such \textit{generalized type eigen-solutions}, that we may write
    \begin{equation}\label{E:BHE-eigenstate}
        f^{\prime\prime}+\left({-e^{4z}}+K_{3}\, e^{3z}+K_{2}\,e^{2z}+\big(K_{1}\big)_{n,\, \nu}\, e^{z}+K_0\right)f=0,\qquad 0\le\nu\le n,
    \end{equation}
for each integer $n$.
\end{remark}
\bigskip

\section{Fredholm integral equations with symmetric periodic kernels}
\label{S:integral}
It has been conjectured that Heun's and its confluent equations, including BHE of course, have no general integral representations of its solution in terms of \textit{simpler} integrands {in terms of ${}_2\,F_1$ or its confluent forms} (see also \cite{Ronv1995}). Instead, it is observed that there are new structures of having \textit{homogeneous Fredholm-type integral equations} (or \textit{integral equation of second kind}) of the form
    \begin{equation}
        \label{E:real-integral}
            u(x)=\lambda\int_a^b K(x,\, t)\, u(t)\, dt,
    \end{equation}
where $\lambda$ is the eigenvalue of the solution $u(x)$, and {the kernel $K(x,\, t)$ is symmetric in $x$ and $t$}, if any. Such integral equations are of fundamental importance for Heun equations (see, for example, \cite[Part A]{Ronv1995}), {which play the role of integral representations for hypergeometric equation.}
\medskip

It has long been known that \textit{Laplace's} and analogous transforms method can solve {some} linear differential equations by definite integrals. The method involves first to write the differential equation in its adjoint form and use it to construct an auxiliary partial differential equation which is then solved by the method of \textit{separation of variables}. See for examples, Forsyth \cite[p. 252]{Forsyth1902} and Ince \cite[XVIII]{Ince1926}. The methods differ from Green's function consideration where the kernel of the integrals involved usually have discontinuities along the ``diagonal", while the Laplace-type integral transform method produces ``continuous" kernels. Our principal concern here is to obtain such an explicit Fredholm-type integral equation and solutions for the {periodic} BHE (\ref{E:BHE-periodic}), using the Laplace's {(integration-by-parts)} method. Of course, one should also investigate the existence and uniqueness problems of the (\ref{E:real-integral}), see for example \cite{Smithies1958}, %\cite{PortStir1990},
but the existing theories do not seem to be applicable for the BHE or the periodic BHE. On the other hand, there are non-linear integral equations and relations exist for the Heun equation, as described in \cite{Ronv1995}, that are not treated here for the periodic BHE.
\medskip

Whittaker applied the above methodology to produce Fredholm integral equations of the second kind. We have already mentioned that he found the integral equation (\ref{E:lame-integral}) for the Lam\'e equation (\ref{E:lame-eqn}) in the \S1. It turns out that another such integral equation that he derived for the Mathieu equation is closer in spirit to the Fredholm integral equation of the second kind for PBHE that we have derived below. Let us review Whittaker's result of   1912 \cite{Whitt1912}. Whittaker showed that one can write eigen-function solution to the Mathieu equation
    \begin{equation}
        \label{E:Mathieu}
        f^{\prime\prime}(z)+(a+k^2\cos^2 z) f(z)=0
    \end{equation}
to satisfy the homogeneous form of \textit{Fredholm integral equation of the second kind} \cite{Smithies1958}
    \begin{equation}
        \label{E:mathieu-integral}
        f(z)=\lambda\int_0^{2\pi}e^{k\cos z\cos s}\, f(s)\, ds,
    \end{equation}
where there is an increasing sequence of eigenvalues $\lambda_1<\lambda_2<\lambda_3<\cdots$ and their corresponding eigen-functions $y_1,\, y_2,\, y_3,\cdots$. An example, is given by Whittaker \cite{Whitt1912} {for the zeroth even Mathieu functions:}
    \begin{equation}
        \label{E:ce_0_integral}
            ce_0(z)=\lambda\int_0^{2\pi}e^{k\cos z\cos s}\, ce_0(s)\, ds,
    \end{equation}
where $\lambda=\displaystyle\frac{ce_0(\pi/2)}{2\pi}$, and
    \[
        ce_0(z)=1+\frac{k^2}{8}\cos 2z+\frac{k^4}{2^9}\cos 4z+\cdots.
    \]

The corresponding Fredholm integral equation we have found for the PBHE is given in the following result where we have assumed that all the coefficients in (\ref{E:BHE-periodic-normal}) are \textit{real}.
\bigskip


\begin{theorem}\label{T:periodic_integral_Eqn} {Let $K_3,\, K_2$ and $K_0<0$ be real and that
        \begin{equation}
                \label{E:distinct-2}
                    1\pm 2\sqrt{-K_0} >0.
        \end{equation}
}
For each non-negative integer $n$, there are $n+1$ distinct pairs of \textit{generalized eigenvalues} $\big((K_1)_{n,\, \nu},\,\lambda_{n,\, \nu}\big),\ 0\le\nu\le n$ such that the Fredholm integral equation of the second kind
    \begin{equation}
        \label{E:periodic-integral-eqn}
            f(z)=\lambda\int_0^{2\pi}
            e^{\frac{i}{2}(e^{2iz}+e^{2is})-K_1(e^{iz}+e^{is})}\, f(s)\, ds
    \end{equation}
admitting corresponding eigen-solutions
    \begin{equation}
        \label{E:BHE-periodic-soln_II}
         BH_{n,\, \nu}(iz)=Y_{n,\,\nu}(e^{-iz})\, \exp\Big(-{(K_1)_{n,\, \nu}}\,e^{iz}-\frac{1}{2}\, e^{2iz}+(n+\sqrt{-K_0})\,iz\Big),\quad
    \begin{cases}
        n=0,\,1,\, 2,\, 3,\cdots & \\
        \nu=0,\,1,\cdots,  n &
    \end{cases}
    \end{equation}
where $(n+\sqrt{-K_0}+1)^2=1$ and the $Y_{n,\,\nu}(\zeta)$ is defined in (\ref{E:BHE-periodic-soln}) and  $(K_1)_{n,\, \nu}$ satisfies the determinant
        \[
            \det\big( (K_1)_{n,\, \nu}\big)=0
        \]
given in {(\ref{E:Det})}.
\end{theorem}

\medskip

\begin{remark}
        \begin{enumerate}
            \item We first note that we will actually prove below that the eigen-solutions  (\ref{E:BHE-periodic-soln_II}) for (\ref{E:BHE-periodic-soln_II}) are actually eigen-solutions (\ref{E:BHE-periodic-soln_II}) to (\ref{E:periodic-integral-eqn}).
            \item There are two choices of the values of $n+\sqrt{-K_0}$, namely $0$ or $-2$ in (\ref{E:BHE-periodic-soln_II}) above.
			\item We further note that the kernel function under the integral sign of our result is closer in spirit to that of the Mathieu equation in (\ref{E:mathieu-integral}) then that of the Lam\'e integral equation.
We also note that the result differs from the eigen-values problem from the Mathieu equation described above where  a single increasing sequence of eigenvalues were found.
            \item Finally, we remark although the eigen-value pairs\\ $\big((K_1)_{n,\, \nu},\,\lambda_{n,\, \nu}\big),\ 0\le\nu\le n$  are distinct, we only know the $\big((K_1)_{n,\, \nu}\big),\ 0\le\nu\le n$ are distinct but we have no result for the corresponding $\lambda$'s at this moment.
        \end{enumerate}
\end{remark}

\bigskip

\begin{proof} Without loss of generality, we may consider for each integer $n\in\mathbb{N}$, and $\nu=0,\, 1,\, \cdots, n$, the coefficients $K_3,\, K_2,\, K_0$ satisfy the relation (\ref{E:condition_I}), then there are $n+1$ of $K_1=(K_1)_{n,\, \nu}$ that are roots to the determinant $D_{n+1}(K_1)=0$ from (\ref{E:Det}),
the equation
\begin{equation}
    \label{E:BHE-periodic_rotated}
        \tilde{f}^{\prime\prime}(z)+\left(\tilde{K}_{4}\,e^{4iz}+\tilde{K}_{3}\, e^{3iz}+\tilde{K}_{2}\,e^{2iz}+\tilde{K}_{1}\, e^{iz}+\tilde{K}_0\right)\tilde{f}(z)=0,
    \end{equation}
$\tilde{f}(z)=f(iz)$, where
        \begin{equation}
            \label{E:first-set}
            \tilde{K}_j=-K_j,\qquad  j=0\le j\le 4,
        \end{equation}
$K_1=(K_1)_{n,\, \nu}$ and in particular $\tilde{K}_4=1$.
\bigskip

Suppose (\ref{E:BHE-periodic_rotated}) admits an ``eigen-solution" $u(z)$. Then we define a sequence of second order partial differential operators
    \begin{equation}
        L_z:=\frac{\partial^2}{\partial z^2}+\ell(z),
    \end{equation}
where {$L_z=(L_{n,\, \nu})_z$ and}
    \begin{equation}
        \ell(z):={(\ell_{n,\, \nu})(z)=}\tilde{K}_{4}\,e^{4iz}+\tilde{K}_{3}\, e^{3iz}+\tilde{K}_{2}\,e^{2iz}+\tilde{K}_{1}\, e^{iz}.
    \end{equation}
Let $K(z,\, s)$ be a function with two complex variables. Then we construct a partial differential equation for which $K(z,\, s)$ solves:
    \begin{equation}
        \label{E:PDE}
            L_z(K)-L_s(K)=\frac{\partial^2K}{\partial z^2}-\frac{\partial^2 K}{\partial s^2}+[\ell(z)-\ell(s)]\, K.
    \end{equation}
Now let
    \begin{equation}
        \label{E:kernel}
        K(z,\, s)=\exp\big[a\,(e^{2iz}+e^{2is})+c\,(e^{iz}+e^{is})\big],
    \end{equation}
with the constants {$a,\, c$} {remain to be chosen}.
\medskip


Substituting (\ref{E:kernel}) into the  (\ref{E:PDE}) yields
    \[
        L_z(K)-L_s(K)=K_{zz}-K_{ss}+\big[\ell(z)-\ell(s)\big]\cdot K
    \]
and which vanishes identically if we set

    \begin{equation}
        \label{E:two-set-parameters-I}
        {4a^2=\tilde{K}_4=1},\quad 4ac=\tilde{K}_3,\quad 4a+c^2=\tilde{K}_2,
        \quad c=\tilde{K}_1.
    \end{equation}
    It follows that
        \begin{equation}
        \label{E:K}
            K_3=4a \, K_1.
        \end{equation}
\medskip

On the other hand, we know from \cite[p. 62]{DMR1978_II} that the transformation
    \begin{equation}
        \label{E:BHE-canonical-BHE}
             \tilde{u}(z)=z^{-\frac{\alpha}{2}}\,e^{\frac{\beta}{2}z+\frac12z^2}\,\tilde{\Psi}(z)
    \end{equation}
transforms the BHE (\ref{E:BHE}) into the equation
    \begin{equation}
        \label{E:Generalized-Bessel-5}
            x^2\tilde{\Psi}^{\prime\prime}(x)+x\tilde{\Psi}^{\prime}(x)+\Big(\tilde{K_4} x^{4}+\tilde{K_3}x^3+\tilde{K_2}x^2+\tilde{K_1}x+\tilde{K_0}\Big)\,\tilde{\Psi}(x)=0,
    \end{equation}
which is a special case of (\ref{E:Generalized-Bessel}) when $\ell=4$, and with the identifications

    \begin{equation}
        \label{E:two-set-parameters-III}
        \begin{array}{lll}
            \tilde{K}_4=1; & \tilde{K}_3=\beta, & \tilde{K}_2=-\Big(\gamma-\displaystyle\frac{\beta^2}{4}\Big),\\
            \tilde{K}_1=\displaystyle\frac{\delta}{2}, & \tilde{K}_0=\displaystyle\frac{\alpha^2}{4}.&\\
        \end{array}
    \end{equation}
between the two sets of parameters of the two equations (please also refer to (\ref{E:first-set}) and (\ref{E:two-set-parameters-II})).
\medskip


Let us first consider the {following unconventional} periodic boundary value problem of
    \begin{equation}
        \label{E:BVP_I}
            \begin{array}{lll}
                ((L_{n,\, \nu})_x+(\tilde{K_0})_n)\, \tilde{u}_{n,\, \nu}=0, & \textrm{on} & \ 0\le x\le 2\pi,\\
                &&\\
                \tilde{u}_{n,\, \nu}(0)=\tilde{u}_{n,\, \nu}(2\pi), & \textrm{and} & \tilde{u}_{n,\, \nu}^\prime(0)=\tilde{u}_{n,\, \nu}^\prime(2\pi)\\
            \end{array}
    \end{equation}
{where $0\le\nu\le n$ and $n\in\mathbb{N}$}.
%admits  an ordered sequence of eigenvalues $\{(\tilde{K_0})_n\}$ \textcolor{blue}{\cite[pp. 736--739]{MorseFesh}, whose absolute values has a lower bound but without upper bound, corresponding to eigen-functions \textcolor{blue}{$\tilde{u}_{n,\, \nu}$} that satisfies the BVP (\ref{E:BVP_I})}. Suppose \textcolor{blue}{$(\tilde{K_0})_n$ and $\tilde{u}_{n,\, \nu}$} is such a pair to the (\ref{E:BVP_I}). So we have}
%   \textcolor{blue}{
%   \begin{equation}
%       \label{E:BVP_II}
%           L_s(\tilde{u}_n)+(\tilde{K_0})_n \, \tilde{u}_{n,\, \nu}=0.
%   \end{equation}
%   }
\medskip

We now apply the so-called classical \textit{Lagrange adjoint equation} method which is
  quite effective to find integral transforms or integral equations formulae as formulated in Forsyth \cite[pp. 251--253]{Forsyth1902}. The current variation is inspired from the papers by Whittaker {\cite{Whitt1912} and \cite{Whitt1915}} and Ince \cite{Ince1922} as well as the \textit{approximate square-root} method in \cite{Bank1993} and \cite{Chiang2000}.

We note that it is clear that the kernel $K$ constructed above {depends on $n$ and $\nu$ which we de-emphasis and its partial derivative in $x$ (or its partial derivative in $s$) are equal at the periodic boundary points:}
    \begin{equation}
        \label{E:periodic-boundary-K}
            K(x,\, x_1)=K(x,\, x_2)\qquad \textrm{and}\qquad
            K_x(x,\, x_1)=K_x(x,\, x_2).
    \end{equation}
So let {$\tilde{u}_{n,\, \nu}$ be the $(n,\,\nu)$-th eigen-solution to (\ref{E:BHE-periodic_rotated}) with $(\tilde{K}_0)_{n,\, \nu}$ for an arbitrary $n=1,\, 2,\, 3,\, \cdots;\ 0\le \nu \le n$. We also let}
    \begin{equation}
        \label{E:integral}
            I(x):=\int_{0}^{2\pi} K(x,\,s)\, {\tilde{u}_{n,\, \nu}(s)}\, ds.
    \end{equation}
Applying the operator $L_x$ onto $I(x)$ yields
    \begin{align}
        \label{E:concomitant-comp}
        L_x[I(x)] &= \int_{0}^{2\pi} L_x[K(x,\,s)]\,{\tilde{u}_{n,\, \nu}(s)}\, ds\notag\\
        &=\int_{0}^{2\pi} L_s[K(x,\,s)]\, {\tilde{u}_{n,\, \nu}(s)}\, ds\notag\\
        &=\big[C(x,\, s)\big]_{0}^{2\pi}+\int_{0}^{2\pi} K(x,\, s)\,L_s[{\tilde{u}_{n,\, \nu}(s)}]\, ds\notag\\
        &={\big[C(x,\, s)\big]_{0}^{2\pi}}-{(\tilde{K}_0)_{n,\, \nu}}\, I(x),
    \end{align}
where $C$ is the \textit{bilinear concomitant} in the form
    \begin{equation}
        \label{E:concomitant}
            C(x,\,s):={\tilde{u}_{n,\, \nu}(s)}K_s(x,\, s)-{\tilde{u}_{n,\, \nu}^\prime(s)}K(x,\, s).
    \end{equation}
{In order to exhibit that the (\ref{E:integral}) is an eigen-solution to (\ref{E:BVP_I}), we need to show that the concomitant (\ref{E:concomitant-comp}) vanishes for {$x$ at both the} end points $0,\, 2\pi$. This follows because of the fact that the $K$ satisfies the periodic boundary values assumption
(\ref{E:periodic-boundary-K}) stated in (\ref{E:BVP_I}).} It remains to verify that this assumption does indeed {hold for $\tilde{u}_{n,\, \nu}(s)$ for each $n=1,\, 2,\, 3,\, \cdots;\ 0\le \nu \le n$}. . To do so, let us first recall from the Theorem \ref{T:determinant_soln} that the periodic BHE (\ref{E:BHE-periodic-normal}) admits an eigen-solution (\ref{E:BHE-periodic-soln}) when (\ref{E:condition_I}) and equivalently (\ref{E:terminates}) hold. In particular, we deduce from (\ref{E:two-set-parameters-I}) and (\ref{E:two-set-parameters-II}) that
    \begin{equation}
        \label{E:2m}
        \begin{split}
         2n &={\gamma}-{\alpha}-2 = -\tilde{K}_2+\frac{\tilde{K}_3^2}{4}-2\sqrt{\tilde{K}_0}-2\notag\\
            &=-(4a+c^2)+\frac{16\,a^2c^2}{4}-2\sqrt{\tilde{K}_0}-2\\
            &=-(4a+c^2)+4a^2c^2-2\sqrt{\tilde{K}_0}-2\\
            &=-4a-2\sqrt{\tilde{K}_0}-2
        \end{split}
    \end{equation}
and $4a^2=-\tilde{K}_4=1$ hold, implying that
		\begin{align*}
            \sqrt{\tilde{K}_0} &= -2a-1-n\\
            &=\pm 1-1-n\\
        \end{align*}
is an integer. Thus in order for the periodic BHE (\ref{E:BHE-periodic_rotated}) to admit the corresponding eigen-solution $\tilde{u}_{n,\,\nu}(z):=\tilde{f}_n(z)=f_n(iz)$  that meets the periodicity assumption (\ref{E:concomitant}) that $C(0,\, s)=C(2\pi,\, s)$, it is sufficient to show that both $\tilde{u}_{n,\,\nu}$ and $\tilde{u}^\prime_{n,\,\nu}$ are periodic of period $2\pi$ so that the (\ref{E:concomitant}) is periodic with the same period. But this now follows from the $(n+\sqrt{-K_0}+1)^2=(n+\sqrt{\tilde{K}_0}+1)^2= 1$. That is, we deduce that $n+\sqrt{-K_0}$ equals either $0$ or $-2$, that is , it equals an integer and its appearance in (\ref{E:BHE-periodic-soln}) implies that the (\ref{E:BHE-periodic-soln}) is indeed periodic. The above proves that $I(x)$ is indeed a solution to the BVP (\ref{E:BVP_I}). We also note that the $\tilde{K}_0=(n+1\pm1)^2$ is given by a \textit{monotone} sequence of eigenvalues. Finally it follows from the equation (\ref{E:K}) that the eigen-solution (\ref{E:BHE-periodic-soln}) assume the form (\ref{E:BHE-periodic-soln_II}) where  $\det( (K_1)_{n,\, \nu})=0$ holds.
\end{proof}
\medskip

{
\begin{remark}  The unusual appearance of eigen-values pairs $\big((K_1)_{k,\, \nu},\,\lambda_{k,\, \nu}\big)$, \\$0\le\nu\le k$ given in the Theorem \ref{T:periodic_integral_Eqn} appears to be a result of the BHE/PBHE possess not just one parameter but several related parameters  (i.e., $\alpha,\, \beta,\, \gamma,\, \delta$) in the differential equations that plays the role of single spectral parameter in their respective eigenvalues problems in the classical equations, such as the Hermite and Laguerre equations. Acscott \cite{Arscott1964}, Sleeman and his co-workers have developed a \textit{multiparameter spectral theory} and applied it to Lam\'e equation and \textit{ellipsoidal wave equation}. In particular, they have related the spectral problems for differential equations to their respective integral equations. See, e.g. \cite{Sleeman1978} and \cite{BrowneSleeman1980}.
 However, it is not clear if their theories could apply to the BHE/PBHE as we have formulated in the current paper.
\end{remark}
}
%\vfill\eject


\section{Single and double orthogonality properties}\label{S:orth}
     It is known that the classical Lam\'e polynomials \cite[\S9.2--9.3]{Arscott1964}, which are regarded as eigen-solutions to the Lam\'e equation
    \begin{equation}
        \label{E:lame-2}
            w^{\prime\prime}(z)+[h-n(n+1)\,k^2\textrm{sn}^2z]\,w(z)=0,
    \end{equation}
where $n$ is a chosen non-negative integer, then there are precisely $n+1$ Lam\'e polynomials that satisfies orthogonality relation \cite[\S9.4]{Arscott1964}
    \begin{equation}
        \label{E:lame-orth}
            \int_{-2K}^{2K} E_n^{m_1}(z) \, E_n^{m_2}(z)\, dz=0
    \end{equation}
whenever $m_1\not= m_2$($0\le m_1,\,m_2\le n$), where the \textit{real} and \textit{imaginary periods} of $\textrm{sn\,}z$ are denoted by $4K$ and $2iK^\prime$ respectively.
\medskip

We show that the PBHE (\ref{E:BHE-periodic}) with solutions that have polynomial component (\ref{E:BHE-periodic-soln_II}) also exhibit an orthogonality relation that is different from the one given in (\ref{E:complex-orth}).

\begin{theorem}\label{T:simple-orth} Let $n$ be a non-negative integer. Suppose the coefficients $K_3,\, K_2$ and $K_0$ of the PBHE
    \begin{equation}\label{E:BHE-eigenstate-II}
        f^{\prime\prime}+\left(-e^{4z}+K_{3}\, e^{3z}+K_{2}\,e^{2z}+\big(K_{1}\big)_{n,\, \nu}\, e^{z}+K_0\right)f=0,\qquad 0\le\nu\le n,
    \end{equation}
satisfy the relation (\ref{E:condition_I}) for each $n$. Then
        {
        \begin{enumerate}
            \item  there are $n+1$, possibly repeated, choices of $\big(K_{1}\big)_{n,\, \nu},\ (0\le \nu\le n)$
            \item and if in addition that $K_3,\, K_2$ and $K_0<0$ are real and $1\pm 2\sqrt{-K_0}>0$
, then there are {exactly $n+1$ distinct real} values of $\big(K_{1}\big)_{n,\, \nu},\ (0\le \nu\le n)$
        \end{enumerate}
        }
such that the corresponding solutions $BH_{n,\, \nu},\ (0\le\nu\le n)$ to (\ref{E:BHE-eigenstate-II}) given in (\ref{E:BHE-periodic-soln}) satisfy the orthogonality relation
    \begin{equation}
		\label{E:simple-orth-2}
            \int_0^{2\pi i} BH_{n,\, \mu}(z)\cdot BH_{n,\, \nu}(z) \,e^z\, dz=0
    \end{equation}
whenever $\mu\not=\nu$.
\end{theorem}
\medskip

\begin{proof} Apart from some minor variation specifically for the periodic BHE, the argument follows from a conventional approach of proving orthogonality. Given a non-negative integer $n$, then once the coefficients $K_3,\, K_2$ and $K_0$ are determined as stated, then we immediately determine the $n+1$ {distinct real} values of $\big(K_{1}\big)_{n,\, \nu}$ from the determinant (\ref{E:Det}). Let us consider
        \begin{equation}
            \begin{split}
                    f_{n,\, \mu}^{\prime\prime}(z) &+\left(-e^{4z}+K_{3}\, e^{3z}+K_{2}\,e^{2z}+\big(K_{1}\big)_{n,\, \mu}\, e^{z}+K_0\right)\,f_{n,\, \mu}(z)=0\\
                    f_{n,\, \nu}^{\prime\prime}(z) &+\left(-e^{4z}+K_{3}\, e^{3z}+K_{2}\,e^{2z}+\big(K_{1}\big)_{n,\, \nu}\, e^{z}+K_0\right)\,f_{n,\, \nu}(z)=0
            \end{split}
        \end{equation}
\medskip

Subtracting the two equations resulting from multiplying the first one by $f_{n,\, \nu}$ and the second one by $f_{n,\, \mu}$ respectively and integrating the resulting equation from $0$ to $2\pi i$ yields
        \begin{equation}
            \begin{split}
            &\int_0^{2\pi i} \big[f_{n,\, \nu}(z)f_{n,\, \mu}^{\prime\prime}(z)-f_{n,\, \nu}(z)f_{n,\, \mu}^{\prime\prime}(z)\big]\, dz\\
            &\qquad +\big[\big(K_{1}\big)_{n,\, \mu}-\big(K_{1}\big)_{n,\, \nu}\big]\,\int_0^{2\pi i}  f_{n,\, \nu}\,(z)f_{n,\, \mu}(z) \,e^z\, dz=0.
            \end{split}
        \end{equation}
Since the $f_{n,\, \nu}(z)$, $f_{n,\, \mu}(z)$ and their first derivatives are periodic of period $2\pi i$, so integration-by-parts yields
        \begin{equation*}
            \begin{split}
             &\int_0^{2\pi i}  f_{n,\, \nu}(z)f_{n,\, \mu}(z) \,e^z\, dz \\
            &={\big[\big(K_{1}\big)_{n,\, \mu}-\big(K_{1}\big)_{n,\, \nu}\big]^{-1}} \int_0^{2\pi i}  \big[f_{n,\, \nu}(z)f_{n,\, \mu}^{\prime\prime}(z)-f_{n,\, \nu}(z)f_{n,\, \mu}^{\prime\prime}(z)\big]\, dz\\
             &={\big[\big(K_{1}\big)_{n,\, \mu}-\big(K_{1}\big)_{n,\, \nu}\big]^{-1}} \big[f_{n,\, \nu}(z)f_{n,\, \mu}^{\prime}(z)-f_{n,\, \nu}(z)\,f_{n,\, \mu}^{\prime}(z)\big]_0^{2\pi i}\\
            &=0
            \end{split}
        \end{equation*}

\end{proof}
\bigskip

Although the above theorem proves orthogonality for the exponential-type Hautot polynomials without restricting the coefficients $\{K_3,\, K_2,\, K_1,\, K_0\}$ are real, these polynomials will only be really orthogonal when these coefficients are real. For if all the $\{K_3,\, K_2,\, K_0<0\}$  are real and $1+\pm\sqrt{-K_0}>0$, then the Theorem  \ref{T:determinant_soln} asserts that there are $n+1$ of $(K_1)_{n,\, \nu}$ so obtained will also be real and distinct. That is, there are precisely $n+1$ exponential-type Hautot polynomials solutions to the PBHE. Besides, a careful examination on the equations in (\ref{E:two-set-parameters-II}) indicates that all the corresponding coefficients $\{\alpha,\, \beta,\, \gamma,\, \delta\}$ are also real. As a result, we deduce from the three-term recursion (\ref{E:BHE-recursion}) that all the coefficients of the (exponential-type) Hautot polynomials are real. Thus, when $\nu=\mu$, we have
	\begin{equation}
		\int_0^{2\pi i} \big(BH_{n,\, \nu}(z)\big)^2\, e^z\, dz \not= 0.
	\end{equation}
Thus, we have established that under the above assumption on the reality of the coefficients on $K_j$, the exponential-type Hautot polynomials of degree $n$, $\{BH_{n,\, \nu}\}_{\nu=0}^n$ are orthogonal with respect to the the \textit{complex weight} $e^z$ over $[0,\, 2\pi i]$.
\bigskip

We further observe that one can define the inner product by integrating the Lam\'e polynomials over a \textit{complex period} $[K-2iK^\prime,\, K+2iK^\prime]$ instead of the real period $[-2K,\, 2K]$ \cite[\S9.4]{Arscott1964}
considered above. That is,
	 \begin{equation}
        \label{E:lame-orth-2}
            \int_{K-2iK^\prime}^{K+2iK^\prime} E_n^{m_1}(z) \, E_n^{m_2}(z)\, dz=0
    \end{equation}
whenever $m_1\not= m_2$. 
It is not difficult to see that analogue situation also holds for the exponential-type Hautot polynomials:	
	\begin{equation}
        \label{E:simple-orth-3}
            \int_{\pi}^{\pi+2\pi i} BH_{n,\, \mu}(z)\, BH_{n,\, \nu}(z) \,e^z\, dz=0
    \end{equation}
whenever $\mu\not=\nu$.
\bigskip


We observe from (\ref{E:lame-orth}) and (\ref{E:simple-orth-2}) that although the corresponding Lam\'e polynomial solutions and periodic BHE polynomial solutions of same degree respectively are orthogonal with respect to the parameters $\mu,\, \nu$, it is not clear if polynomials of different degrees are orthogonal to each other. It turns out that an orthogonality exists for \textit{polynomials} of different degrees when they are formed from product of two Lam\'e polynomials of same degree but in different variables. The following \textit{double orthogonality} can be found in \cite[\S9.4]{Arscott1964}:
    \begin{equation}
        \label{E:lame-double-orth}
            \int_{-2K}^{2K}\int_{K-2iK^\prime }^{K+2iK^\prime} E_m^\mu(z)\, E_m^\mu(s)\,E_n^\nu(z)\, E_n^\nu(s)\,
            (\textrm{sn\,}^2z-\textrm{sn\,}^2s)\,dz\, ds=0
    \end{equation}
\medskip
whenever $n\not= m$, while if $m=n$, then the orthogonality still hold when $\mu\not=\nu$, for
$0\le\mu,\, \nu\le n$. We show below that a corresponding double orthogonality can be extended to PBHE:
\medskip

\begin{theorem}
    \label{T:double-orth}
     Let $K_3$ and $K_2$ be given {real numbers}. Let $m$ and $n$ be non-negative integers such that {$(K_0)_n<0$ and $(K_0)_m<0$} satisfy
     \begin{equation}
        \label{E:condition_II}
            \frac{K_3^2}{4}+K_2\pm 2\sqrt{-(K_0)_n}=2(n+1),\quad 1\pm 2\sqrt{-(K_0)_n}>0
    \end{equation}
 and
     \begin{equation}
        \label{E:condition_III}
            \frac{K_3^2}{4}+K_2\pm 2\sqrt{-(K_0)_m}=2(m+1), \quad 1\pm 2\sqrt{-(K_0)_m}>0
    \end{equation}
respectively. Suppose $BH_{n,\, \mu}\ (0\le \mu\le n)$ are solutions of the differential equation (\ref{E:BHE-eigenstate-II}) as defined in the Theorem \ref{T:simple-orth}, and  $BH_{m,\,\nu}\ (0\le \nu\le m)$ are the corresponding solutions to the equation (\ref{E:BHE-eigenstate-II}) with $n$ replaced by $m$. Then we have
        \begin{equation}
            \label{E:double-orth}
             \int_0^{2\pi i}\int_\pi^{\pi+2\pi i} BH_{n,\, \mu}(z)BH_{n,\, \mu}(s)\, BH_{m,\, \nu}(z)BH_{m,\, \nu}(s)\,
              (e^z-e^s)\, dz\,ds=0
        \end{equation}
whenever $(n,\, \mu)\not=(m,\, \nu)$.
\end{theorem}
\medskip

\begin{proof} The idea is along the classical treatment by constructing a suitable pair of partial differential equations and applying integration-by-parts. Let
        \[
            F_{n,\, \mu}(z,\, s):=f_{n,\, \mu}(z) f_{n,\, \mu}(s),\qquad
            F_{m,\, \nu}(z,\, s):=f_{m,\, \nu}(z) f_{m,\, \nu}(s)
        \]
be functions of two variables $(z,\, s)$. Clearly they satisfy the following partial differential equations
        \begin{align}
                \frac{\partial^2 F_{n,\, \mu}}{\partial z^2}-\frac{\partial^2 F_{n,\, \mu}}{\partial s^2}
                    &=\Big[-(e^{4z}-e^{4s})+ K_3\,(e^{3z}-e^{3s})\notag\\
                    &\qquad+K_{2}\,(e^{2z}-e^{2s})+\big(K_{1}\big)_{n,\, \mu}\, (e^{z}-e^s)\Big] \, F_{n,\, \mu}=0\label{E:PDE_I}\\
                    \frac{\partial^2 F_{m,\, \nu}}{\partial z^2}-\frac{\partial^2 F_{m,\, \nu}}{\partial s^2}
                    &=\Big[-(e^{4z}-e^{4s})+ K_3\,(e^{3z}-e^{3s})\notag\\
                    &\qquad +K_{2}\,(e^{2z}-e^{2s})+\big(K_{1}\big)_{m,\, \nu}\, (e^{z}-e^s)\Big] \, F_{m,\, \nu}=0\label{E:PDE_II}
            \end{align}
\medskip

Subtracting the equations (\ref{E:PDE_I}) after multiplying  by $F_{m,\, \nu}$ and the equation (\ref{E:PDE_II}) after multiplying by $F_{n,\, \mu}$ and integrating the two variables from $0$ to $2\pi i$ and $\pi$ to $\pi+2\pi i$ respectively results in
	   \begin{equation*}
                \begin{split}
                     &\int_0^{2\pi i}\int_\pi^{\pi+ 2\pi i}\Big\{
                        F_{m,\, \nu}\big[(F_{n,\, \mu}\big)_{zz}-(F_{n,\, \mu})_{ss}\big]
                        -
                        F_{n,\, \mu}\big[(F_{m,\, \nu}\big)_{zz}-(F_{m,\, \nu})_{ss}\big]\Big\}\, dz\, ds\\
                        &+
                        \big[(K_1)_{n,\,\mu}-(K_1)_{m,\,\nu}\big]\int_0^{2\pi i}\int_\pi^{\pi+2\pi i}
                        F_{n,\, \mu}(z,\, s)\,F_{m,\, \nu}(z,\, s) \, (e^{z}-e^s)\, dz\, ds=0.
                \end{split}
            \end{equation*}
Integration of the above equation by parts and using the fact that both the $F_{n,\, \mu}$, $F_{m,\, \nu}$ and their partial derivatives are periodic with respective to both the variables $(z,\, s)$ yields
		\begin{equation}
            \begin{split}
                0&=\int_0^{2\pi i}\int_\pi^{\pi+2\pi i}\big[F_{m,\, \nu}\,(F_{n,\, \mu}\big)_{zz}-
                        F_{n,\, \mu}\,(F_{m,\, \nu})_{zz}\big]\\
                 &\quad+\big[F_{n,\, \mu}\, (F_{m,\, \nu})_{ss}-F_{m,\, \nu}\,(F_{n,\, \mu}\big)_{ss}\big]\,
                        dz\,ds\\
                &\quad+\big[(K_1)_{n,\,\nu}-(K_1)_{m,\,\mu}\big]\int_0^{2\pi i}\int_\pi^{\pi+2\pi i}
                        F_{n,\, \mu}(z,\, s)\,F_{m,\, \nu}(z,\, s) \, (e^{z}-e^s)\, dz\, ds\\
                &=  \int_0^{2\pi i}\big[F_{m,\, \nu}\,(F_{n,\, \mu}\big)_{z}-
                        F_{n,\, \mu}\,(F_{m,\, \nu})_{z}\big]\big|_{\pi}^{\pi+2 \pi i}\,ds\\
                &\quad -\int_0^{2\pi i}\int_\pi^{\pi+2\pi i} (F_{m,\, \nu})_z\,(F_{n,\, \mu})_{z}-
                        (F_{n,\, \mu})_z\,(F_{m,\, \nu})_{z}\, dz\, ds\\
                &\quad +
                  \int_\pi^{\pi+2\pi i}\big[F_{n,\, \mu}\, (F_{m,\, \nu})_s-F_{m,\, \nu}\,(F_{n,\, \mu})_s\big]\big|_0^{2\pi i}\, dz\\
                & \quad -\int_0^{2\pi i}\int_\pi^{\pi+2\pi i} (F_{m,\, \nu})_s\,(F_{n,\, \mu})_s-
                        (F_{n,\, \mu})_s\,(F_{m,\, \nu})_s\, dz\, ds\\
                &\quad +\big[(K_1)_{n,\,\nu}-(K_1)_{m,\,\mu}\big]\int_0^{2\pi i}\int_\pi^{\pi+2\pi i}
                        F_{n,\, \mu}(z,\, s)\,F_{m,\, \nu}(z,\, s) \, (e^{z}-e^s)\, dz\, ds\\
                &=\int_{0}^{2\pi i}\big[F_{m,\, \nu}\,(F_{n,\, \mu}\big)_{z}-
                        F_{n,\, \mu}\,(F_{m,\, \nu})_{z}\big]\big|_{\pi }^{\pi+2\pi i}\,ds-0\\
                & \quad +
                 \int_\pi^{\pi+2\pi i} \big[F_{n,\, \mu}\, (F_{m,\, \nu})_{s}-F_{m,\, \nu}\,(F_{n,\, \mu})_{s}\big]\big|_0^{2\pi i}\, dz-0\\
                 &\quad+\big[(K_1)_{n,\,\nu}-(K_1)_{m,\,\mu}\big]\int_0^{2\pi i}\int_\pi^{\pi+2\pi i}
                        F_{n,\, \mu}(z,\, s)\,F_{m,\, \nu}(z,\, s) \, (e^{z}-e^s)\, dz\, ds
            \end{split}
        \end{equation}
where the second and the fourth double integrals are obviously zero. In fact, both of the two remaining two single integrals in the last equal sign are also zero since the  $F_{n,\, \mu}$, $F_{m,\, \nu}$ and their partial derivatives are periodic. That is, we have shown that
	\[
		\big[(K_1)_{n,\,\nu}-(K_1)_{m,\,\mu}\big]\int_0^{2\pi i}\int_\pi^{\pi+2\pi i}
                        F_{n,\, \mu}(z,\, s)\,F_{m,\, \nu}(z,\, s) \, (e^{z}-e^s)\, dz\, ds=0.
    \]
 \medskip                       

 This proves the (\ref{E:double-orth}) holds whenever $n\not= m$ and irrespective to the choices of $\nu$ and $\mu$. If, however, $m=n$ in (\ref{E:double-orth}), then we may apply the known single orthogonality relation (\ref{E:simple-orth-2}) in the above argument to show that (\ref{E:double-orth}) again holds whenever $\nu\not=\mu$. We omit the straightforward argument.
\end{proof}
\bigskip

We see that if $n=m$ and $\nu=\mu$, then we have
	\begin{equation*}
		\begin{split}
			&\int_0^{2\pi i}\int_\pi^{\pi+2\pi i}
                        F_{n,\, \mu}(z,\, s)\,F_{m,\, \nu}(z,\, s) \, (e^{z}-e^s)\, dz\, ds\\
             &=\int_0^{2\pi i}\int_0^{2\pi i}
                        [F_{n,\, \mu}(z,\, \pi+s)]^2 \, (e^{z}-e^\pi\cdot e^s)\, dz\, ds\not =0           
		\end{split}
	\end{equation*}
so that the products $\{BH_{n,\, \mu}(z)\, BH_{n,\,\mu}(s)\}$ forms an orthogonal system.
\bigskip

\section{Quasi-exact solvability}\label{S:quasi}

Lie group methods have helped in computing the spectrums of the Hermite and Laguerre equations which are, respectively, the governing equations for harmonic oscillator and hydrogen atom which are two most important quantum mechanical models at the beginning of quantum mechanics. Turbiner \cite{Turb1988} and his co-workers \cite{GKO1993}, \cite{FGR1996}, etc. have apply Lie algebraic and Lie group methods to study the so-called \textit{quasi-exactly solvable} quantum mechanics problems. We shall show below that a prime example of quasi-exact solvable equation is a BHE and hence the earlier results that were deduced from complex oscillation theory become applicable.
\medskip

Turbiner applied finite group representation of the group $SL(\mathbb{C})$ to generate the so-called \textit{quasi-exactly solvable} quantum mechanics problems. As a result a few Schr\"odinger equations for which a finite number of eigenvalues and eigenfunctions can be solved explicitly are produced in \cite{Turb1988} and \cite{Turb1996}. See also \cite{GKO1993}, \cite{FGR1996} and \cite{GKO1994}.
\medskip



    We propose to consider the following Schr\"odinger equation
    \begin{equation}\label{E:Bender-Dunne-eqn-3}
        {H_c}\psi=E\psi,
    \end{equation}
where
        \begin{equation}\label{E:Bender-Dunne-eqn-4}
            {H_c}=-\frac{d^2}{dx^2}+\frac{(4s-1)(4s-3)}{4x^2}-(4s+4J-2)x^2-cx^4+x^6.
        \end{equation}
which includes the Bender-Dunne equation as a special case when $c=0$, and which can be transformed to the general BHE (\ref{E:BHE}) via two transformations. We first obtain
        \begin{equation}
            \label{E:Bender-Dunne-eqn-5}
                \begin{split}
                x\,g^{\prime\prime}(x)+(-2x^4-cx^2 &+4s-1)\, g^\prime(x)\\
                    &\qquad + \Big((-3+\frac{c^2}{4}+4J)\,x^3+(E-2cs)\, x\Big)\, g(x)=0,
                \end{split}
        \end{equation}
from (\ref{E:Bender-Dunne-eqn-4}) via the transformation
        \begin{equation}
            \label{E:BD-trans-2}
                \psi (x)=e^{-\frac14x^4-\frac{c}{4}x^2}\, x^{2s-1}\, g(x).
        \end{equation}
We then obtain from the equation (\ref{E:Bender-Dunne-eqn-5}) the general BHE with the following coefficients
        \begin{equation}
            \label{E:Bender-Dunne-eqn-6}
                \begin{split}
                z\, u^{\prime\prime}(z) &+\Big(2s-\frac{\sqrt{2}}{2}c\, z-2\, z^2\Big)\, u^\prime(z) \\
                &\qquad+\Big(\frac{\sqrt{2}}{4}\,E-\frac{\sqrt{2}}{2}\,c+\frac{c^2}{8}+(2J-2)z\Big)\, u(z) =0,
                \end{split}
        \end{equation}
via the transformation
        \begin{equation}
            \label{E:BD-trans-2}
                z=\frac{x^2}{\sqrt{2}},\qquad u(z)=g(x).
        \end{equation}
That is, the (\ref{E:Bender-Dunne-eqn-6}) is (\ref{E:BHE}) with
    \begin{equation}
    \begin{array}{ll}
        \label{E:identification}
            \alpha=2s-1, & \beta=\displaystyle\frac{\sqrt{2}}{2}\,c,\\
            \gamma=2J+2s-1, & \delta=-\displaystyle\frac{\sqrt{2}}{2}\,E-\displaystyle\frac{c^2}{4}.
    \end{array}
    \end{equation}


The coefficients of a series solution (\ref{E:infinite-sum}) of the BHE (\ref{E:BHE}) given by Maroni in \cite[p. 163]{DMR1978_I} or \cite[\S (3.1.1)]{Ronv1995} satisfies the three-term recursion (\ref{E:BHE-recursion}) which now assumes the form
    \begin{equation}
        \label{E:infinite-sum-2}
            \begin{split}
        N\Big(2s-1& ,\,\frac{\sqrt{2}c}{2},\,2J  +2s-1,\-\displaystyle\frac{\sqrt{2}}{2}\,E-\displaystyle\frac{c^2}{4};\, \frac{x^2}{\sqrt{2}}\Big)\\
            &=\sum_{k=0}^\infty \frac{A_k}{(1+\alpha)_k\, k!}\, x^{k},
            \end{split}
    \end{equation}
where  the coefficients $A_k,\, k\ge0$ satisfy the three-term recursion formula (\ref{E:BHE-recursion}) with the four parameters chosen from (\ref{E:identification}) :
    \begin{align}
        \label{E:BHE-recursion-2}
             A_{k}&=\Big[(k-1)\frac{\sqrt{2}}{2}\,c+\frac{\sqrt{2}}{2}\, s\, c-\frac{c^2}{8}-\frac{\sqrt{2}}{4}\,E\Big]A_{k-1}\notag\\
                    &\qquad\qquad
                    -(k+1)(k-2+2s)(2J+2-2k)\, A_{k-2},\quad k\ge 2
    \end{align}
where
    \begin{equation}
        \label{E:recursion-initial}
            A_0=1,\qquad A_1=2J-2.
    \end{equation}
Substituting
        \begin{equation}
            \label{E:substit-1}
                A_k=\Big(-\frac{\sqrt{2}}{4}\Big)^k\, P_k^c(E)
        \end{equation}
into (\ref{E:BHE-recursion-2}) yields
        \begin{equation}
            \label{E:BHE-recursion-4}
                \begin{split}
                P_k^c(E) &=\Big[E-\Big(2c(k-1)+2cs-\displaystyle\frac{c^2}{2\sqrt{2}}\Big)\Big]\, P_{k-1}^c(E)\\
                &\qquad\qquad -16 (k-1)(k-2+2s)(J+1-k)\, P_{k-2}^c(E)
                \end{split}
        \end{equation}
where the initial condition is
        \begin{equation}
            \label{E:initial-c}
            P_0^c (E)=1,\qquad P_1^c(E)=E-\Big(2cs-\displaystyle\frac{c^2}{2\sqrt{2}}\Big)
        \end{equation}
and which degenerates into (\ref{E:BHE-recursion-3}) when $c=0$. That is, we have
            \begin{equation}
        \label{E:Bender_Dunne_series_c}
            \psi(x)=e^{-\frac{x^4}{4}-\frac{c}{4}x^2}\, x^{2s-1/2}\, \sum_{k=0}^\infty
            \Big(-\frac14\Big)^k\frac{P_k^c(E)}{k!\,\Gamma(n+2s)}\,
            x^{2k}.
    \end{equation}


It also follows from general theory of orthogonal polynomials that the $\{P_k^c(E)\}$ is a \textit{finite family} of orthogonal polynomials whenever $J+1\le k$. This follows from Favard's theorem \cite[p. 21]{Chihara1978}. We see immediately that the above Bender-Dunne polynomials are special cases of polynomials solutions to the general BHE due to Hautot \cite{Hautot1969} by choosing $c=0$.
\medskip

Moreover, when $J$ is chosen as a positive integer and the case when $c= 0$ above in (\ref{E:factor}), the  $\{P_k^c(E)\}$ exhibits a similar \textit{factorization property} as indicated by Bender and Dunne (\ref{E:factor}) when $c=0$:
        \begin{equation}
        \label{E:factor-c}
            P_{n+J}^c(E)=P_J^c(E)\, Q_n^c(E),\qquad n\ge 0.
    \end{equation}
\medskip

The first few $P_J^c$ are given below when $J=3$,

\begin{align}
P_0^c(E)&=1,\\
P_1^c(E)&=E-2cs+\frac{\sqrt{2}}{4}c^2,\\
P_2^c(E)&=E^2+\Big(\frac{\sqrt{2}}{2}c^2-2c-4cs\Big)\,E \nonumber\\
            &\quad + 4c^2s^2+(-\sqrt{2}c^3+4c^2-64)s+\frac{1}{8}c^4-\frac{\sqrt{2}}{2}c^3,\\
P_3^c(E)&=
E^3+\Big(\frac{3\sqrt{2}}{4}c^2-6c-6cs\Big)E^2\\
    &\quad+\Big(8c^2+\frac{3}{8}c^4-128s-3\sqrt{2}c^3s
    -32+24c^2s+12c^2s^2-3\sqrt{2}c^3\Big)\,E\nonumber\\
    &\quad -8c^3s^3+(256c-24c^3+3\sqrt{2}c^4)s^2+\Big(320c
-32\sqrt{2}c^2-16c^3+6\sqrt{2}c^4-\frac{3}{4}c^5\Big)s
\nonumber\\
    &\quad -8\sqrt{2}c^2
+2\sqrt{2}c^4-\frac{3}{4}c^5+\frac{\sqrt{2}}{32}c^6,\nonumber\\
 P_4^c(E)&=\Big[E-6c-2cs+\frac{\sqrt{2}}{4}c^2\Big]P_3(E),\\
 P_5^c(E)&=\Big[E^2+\Big(-4cs+\frac{\sqrt{2}}{2}c^2-14c\Big)\,E\nonumber\\
    &   \quad +4c^2s^2
+(28c^2+128-\sqrt{2}c^3)s+\frac{1}{8}c^4-\frac{7\sqrt{2}}{2}c^3+48c^2+192\Big]\, P_3(E),
\end{align}
from which we see that
    \[
        P_4^c(E) =P_3^c(E)\, Q_1^c(E),\qquad P_5^c(E)=P_3^c(E)\, Q_2^c(E)
    \]
where
        \[
            Q_1^c(E):=E-6c-2cs+\frac{\sqrt{2}}{4}c^2,
        \]
and
        \[
            \begin{split}
            Q_2^c(E) &:=E^2+\Big(-4cs+\frac{\sqrt{2}}{2}c^2-14c\Big)\,E\\
                &   \quad +4c^2s^2+(28c^2+128-\sqrt{2}c^3)s+\frac{1}{8}c^4-\frac{7\sqrt{2}}{2}c^3+48c^2+192.
            \end{split}
        \]
We deduce immediately that the above (\ref{E:factor-c}) and expressions for $P_4^c(E)$ and $P_5^c(E)$
reduce to Bender and Dunne's (\ref{E:factor}) and (\ref{E:initial}).
\medskip

We note that Maroni has also discussed generating functions for the  $\{P_k^c(E)\}$ in \cite[Part D, \S5]{Ronv1995}.



\section{Periodic Turbiner equation}\label{S:P-Turbiner}

Since the Turbiner equation is a special case of the BHE which we have written it in a periodic form, so we shall rephrase some previous results for the equation in periodic form below. The identification (\ref{E:two-set-parameters-II}) and (\ref{E:identification}) allow to write the equation (\ref{E:Bender-Dunne-eqn-6}) into the form
        \begin{equation}
            \label{E:Periodic-BD-general}
%               \begin{split}
                f^{\prime\prime}+\Big(-e^{4z}-\frac{\sqrt{2}}{2}c\,e^{3z}
                +(2s+2J+1)\,e^{2z}+\frac{\sqrt{2}\,E}{4}e^z-\frac{(2s-1)^2}{4}\Big)
                f(z)=0
        \end{equation}
and it becomes
    \begin{equation}
            \label{E:Periodic-BD}
%               \begin{split}
                f^{\prime\prime}+\Big(-e^{4z}
                +(2s+2J+1)\,e^{2z}+\frac{\sqrt{2}\,E}{4}e^z-\frac{(2s-1)^2}{4}\Big)
                f(z)=0
    \end{equation}
when $c=0$, that is, it is the periodic form of the equation (\ref{E:Bender-Dunne-eqn-1}). We assume all the parameters $c,\, s,\, J$ and $E$ above are real values.
We note that the requirement that $\alpha+1>0$ in order apply Rovder's result to conclude that there are $J+1$ choices of $E$ in the theorems below. But this is automatically satisfied since
    \[
        1+\pm 2\sqrt{-K_0}=1+\pm 2\, \sqrt{\frac{(2s-1)^2}{4}}=2s>0
    \]
provided that $s>0$. We shall omit their proofs.
\medskip

We deduce from the Theorem \ref{E:Bender-Dunne-eqn-1} that
\medskip

\begin{theorem}\label{T:determinant_soln_2}
    Let $c$, $s>0$, and non-negative integer $J$ be arbitrarily chosen. Then  there are $J+1$\footnote{The $J+1$ solutions come from choosing the "$+$" sign in (\ref{E:condition_II})} choices of $E$ that satisfy the relation
        \begin{equation}
            \label{E:condition_II}
                4c^2+2(2s+2J+1)\pm (2s-1)=4(J+1)
        \end{equation}
such that when (\ref{E:condition_II}) holds, then there are precisely $J+1$\footnote{We have a further $J+1$ solutions if choose the ``$-$" sign in (\ref{E:condition_II})}
distinct choices of $E$ which consist of the roots of the determinant $D_{J+1}(E)=0$ where $D_{J+1}(E)$ is given by
%   \begin{align}
    \begin{equation}
        \label{E:Det-2}
        \left|
    \begin{array}{lllll}
        k_1(E) & -1 & & &   \\
        -2s\,k_2 & k_1(E)+\displaystyle\frac{\sqrt{2}c}{2} & -1 & & \\
        & -2(1+2s)(k_2-2) &k_1(E)+2\cdot\displaystyle\frac{\sqrt{2}c}{2} & -1 & \\
        &\ddots & \ddots & \ddots  &\ddots\\
        &\ddots & \ddots & \ddots &\\
        & &
            \begin{split}
                    &-(J-1)(J+2s-2) \\
                    &\quad\times (k_2-2(J-2))
            \end{split}
        & k_1(E)+(J-1)\displaystyle\frac{\sqrt{2}c}{2} & -1\\
        & & &
            \begin{split}
                &-J(J+1+2s)\\
                &\quad\times (k_2-2(J-1))
            \end{split}
        & k_1(E)+J\displaystyle\frac{\sqrt{2}c}{2}
    \end{array}
    \right|,
    \end{equation}
%   \end{align}
and the $k_1$, $k_2$ are given by
    \begin{align}
        \label{E:BD-general-para}
& k_1(E)=\frac{\sqrt{2}}{4}\big(-E+2s),\\
& k_2=2J+\frac{1}{8}\, c^2.
    \end{align}
Moreover, we have
	\begin{equation}
        \label{E:BHE-periodic-soln-4}
        \begin{split}
            BH_{J,\, \nu}(z)&= P_{J,\, \nu}\Big(2s-1,\,\frac{\sqrt{2}}{2}\,c,\, 2(J+s)-3,\, -\frac{\sqrt{2}}{2}\,E;\, e^{z}\Big)\\
                &\qquad\qquad\qquad \cdot\exp\Big[\frac{\sqrt{2}c}{4}\,e^z-\frac{1}{2}e^{2z}+\big(\frac{2s-1}{2}\big)z\Big]\\
            &= Y_{J,\,\nu}(e^{-z})\, \exp\Big[\frac{\sqrt{2}c}{4}\,e^z-\frac{1}{2}e^{2z}+\big(J+\frac{2s-1}{2}\big)z\Big]
        \end{split}
    \end{equation}
where $P_{n,\, \nu}$ are the orthogonal polynomials mentioned in (\ref{E:Real-orth}) and their reversed forms, namely $Y_{J\,\nu}(x)=x^J\, P_{J,\, \nu}(\alpha,\, \beta;\, 1/x)$ are the Bessel orthogonal polynomials of BHE class, with respect to the complex weight defined in Theorem \ref{T:orthogonality}.
\end{theorem}
\bigskip


Clearly a corresponding Fredholm integral equation of second type follows from the Theorem \ref{T:periodic_integral_Eqn} immediately:

\begin{theorem}\label{T:periodic_integral_Eqn_2} Let $c$, $s>0$, and non-negative integer $J$ be arbitrary chosen. For each non-negative integer $J$, there are $J+1$\footnote{The $J+1$ solutions come from choosing the "$+$" sign in (\ref{E:condition_II}).} pairs of \textit{generalized eigenvalues} $\big(\sqrt{2}(E)_{J,\, \nu}/4,\,\lambda_{J,\, \nu}\big),\ 0\le\nu\le J$ such that the Fredholm integral equation of the second kind
    \begin{equation}
        \label{E:periodic-integral-eqn}
            f(z)=\lambda\int_0^{2\pi}
            e^{\frac{i}{2}(e^{2iz}+e^{2is})-\frac{\sqrt{2}E}{4} (e^{iz}+e^{is})}\, f(s)\, ds
    \end{equation}
admitting corresponding eigen-solutions
	\begin{equation}
        \label{E:BHE-periodic-soln_III}
         BH_{J,\, \nu}(iz)=Y_{J,\,\nu}(e^{-iz})\, \exp\Big(\frac{\sqrt{2}(E)_{J,\, \nu}}{4}\,e^{iz}-\frac{1}{2}\, e^{2iz}+\big(J+\frac{2s-1}{2}\big)iz\Big),\quad
    \begin{cases}
        J=0,\,1,\, 2,\, 3,\cdots & \\
        \nu=0,\,1,\cdots  J, &
    \end{cases}
    \end{equation}
where the $Y_{n,\,\nu}(\zeta)$ is defined in (\ref{E:BHE-periodic-soln}) and  ${\sqrt{2}(E)_{J,\, \nu}}/{4}$ satisfies the determinant $\det( (E)_{J,\, \nu})=0$ given in {(\ref{E:Det-2})}.
\end{theorem}
\bigskip

Likewise, the orthogonality relationships for the PBHE derived in the \S8 have their counterparts as stated below. We start with single orthogonality:
\bigskip

\begin{theorem}\label{T:simple-orth-2}
    Let $c$, $s>0$, and non-negative integer $J$ be arbitrarily chosen. Then there are $J+1$ , distinct choices of $E$ that satisfy the relation (\ref{E:condition_II}). The $E$ are the roots of the determinant $D_{J+1}(E)=0$ (\ref{E:Det-2}), such that the corresponding eigen-solutions to $f_{J,\, \nu}$  of the equation (\ref{E:Periodic-BD-general}) satisfy the orthogonality relation
	   \begin{equation}
            \label{E:BD-orth}
                \int_0^{2\pi i} BH_{J,\, \mu}(z)\,BH_{J,\, \nu}(z)\,e^z\, dz=0
        \end{equation}
whenever $\mu\not= \nu$, for $0\le \nu,\, \mu\le J,\ J=1,\, 2,\, 3,\,\cdots$.
\end{theorem}
\bigskip

Just as in the discussion in the \S8, the orthogonality between two eigen-solutions $ BH_{J,\, \mu}(z)\,BH_{J,\, \nu}(z)$ are indistinguishable with respect to a single measure and single integral, we can distinguish orthogonality between the products $BH_{J,\, \mu}(z)BH_{J,\, \mu}(t)$ and $BH_{I,\, \nu}(z)BH_{I,\, \nu}(t)$ and with respect to a mixed weigh and double integral.
\medskip

\begin{theorem}
    \label{T:double-orth-2}
    Let $c$, $s>0$, and non-negative integers $J,\, I$ such that
        \begin{equation*}
%           \label{E:condition_II}
                4c^2+2(2s+2J+1)\pm (2s-1)=4(J+1)
        \end{equation*}
and
    \begin{equation}
            \label{E:condition_III}
                4c^2+2(2s+2I+1)\pm (2s-1)=4(I+1)
        \end{equation}
hold. Then there are precisely $J+1$ distinct choices of $E_1$ and $I+1$ choices of $E_2$ that satisfy, respectively, the determinants $D_{J}(E_1)=0$ given in (\ref{E:Det-2}) and the corresponding $D_{I}(E_2)=0$ (with the suitably defined $k_1$ and $k_2$ in (\ref{E:BD-general-para})). Suppose  that $BH_{J,\, \mu}\ (0\le \mu\le J)$ and $BH_{I,\, \nu}\ (0\le \nu\le I)$ are solutions to the differential equation (\ref{E:Periodic-BD-general}) and the same equation with $J$ replaced by $I$ respectively. Then we have
	\begin{equation}
            \label{E:double-orth}
             \int_0^{2\pi i}\int_\pi^{\pi+2\pi i} BH_{J,\, \mu}(z)BH_{J,\, \mu}(t)\, BH_{I,\, \nu}(z)BH_{I,\, \nu}(t)\,
              (e^z-e^t)\, dz\,dt=0
        \end{equation}
whenever $(I,\, \mu)\not=(J,\, \nu)$.
\end{theorem}
\bigskip

\section{Periodic forms of Heun's Equations}
   Periodic forms of the Heun equation and its confluent forms have been studied long before their algebraic forms. Heun appears to be the first one to study the equation name after him in 1889 \cite{Heun1889}
    \begin{equation}\label{E:GHE}
         \frac{d^2y}{dz^2}+\left(\frac{c}{z}+\frac{d}{z-1}+\frac{e}{z-t}\right)\frac{dy}{dz}+\frac{ab(z-t)-\sigma}{z(z-1)(z-t)}\,y=0,
    \end{equation}
now known as the \textit{general Heun equation} (GHE). However, the well-known  Lam\'e equation \cite[chap. XXIII]{WW27}, \cite[chap. XV]{Bateman} which is a special case of the GHE was written down by Lam\'e in 1839, see \cite[\S23.1]{WW27}. Moreover, Darboux \cite{Darboux1882} has wrote down the equation name after him
\[
                    \begin{split}
                    \frac{d^2y}{dz^2}&=\Big(\frac{\nu(\nu+1)}{\text{sn\,}^2z}
                    +\frac{\nu^\prime(\nu^\prime+1)}{\text{cn\,}^2z}\,\text{dn\,}^2z+\\
                    &\qquad +\frac{\nu^{\prime\prime}(\nu^{\prime\prime}+1)}{\text{dn\,}^2z}
                    k^2\,\text{cn\,}^2z+n(n+1)k^2\text{sn\,}^2z+h
\Big)y.
                    \end{split}
                \]
in 1882 which is a periodic form of the (\ref{E:GHE}) before Heun.
\medskip

Similarly, Mathieu wrote down the Mathieu equation  (\ref{E:Mathieu}) in 1868, see \cite[\S19.1]{WW27} which is a special case of the Prolate spheroidal equation \cite{Slep1978}, \cite[chap. XVI]{Bateman} or called \textit{confluent Heun equation} (CHE)
    \begin{equation}
        \label{E:CHE}
        \frac{d^2y}{dz^2}-\left(\frac{c}{z(z-1)}+\frac{d}{z-1}\right)
\frac{dy}{dz}+\frac{\sigma-(t+\sigma)z}{z(z-1)^2}\,y=0,
    \end{equation}
which has singularities $\{0,\, 1,\, \infty\}$ at which the $\infty$ is an irregular singular point. It results from the coalesces of $t$ and $\infty$ in the Heun equation (\ref{E:GHE}). One can also find a periodic form equation for the CHE from \cite[chap. XVI]{Bateman}.
\medskip

A further coalesce of the singularity $1$ and $\infty$ from the CHE results in the BHE (\ref{E:BHE}) which  is the main focus of this article.


\subsection*{Acknowledgments}
The authors would like to acknowledge very useful conversations with Alex Zhedanov and Peter Olver on the weight function in Theorem \ref{T:orthogonality} and the relation between quasi-exact solvability and Heun equations \cite{OlverP1997} respectively. The authors would also like to acknowledge Chiu-Yin Tsang for pointing out some inaccuracies in an earlier manuscript.



\bibliographystyle{amsplain}
\begin{thebibliography}{10}

\bibitem{AbrSte1972} Abramowitz, M. \and I. Stegun, \textrm{Handbook of Mathematical Functions}, 9th edition, Dover 1972.

\bibitem{AAR1999}  Andrews, G. E., R. Askey \and R. Roy, \textit{Special Functions}, Encyclopedia of Mathematics and its Applications Vol. \textbf{71}, Camb. Univ. Press, Cambridge 1999.

\bibitem{Arscott1964} Arscott, F. M. \textrm{Periodic Differential Equations: An Introduction to Mathieu, Lam\'e and Allied Functions}, Pergamon Press, 1964.

\bibitem{Bank1993} Bank, S. B., \textit{On the explicit determination of certain solutions of periodic differential equations}, Complex Variables Theory Appl. \textbf{23} (1993), no. 1--2, 101--121.

\bibitem{BL1982} Bank, S. B. \and I. Laine,
\textit{On the oscillation theory of {$f\sp{\prime\prime}+{A}f=0$}
where {${A}$} is entire}, Trans. Amer. Math. Soc., \textbf{273}
(1982), 351--363.

\bibitem{BL1983(1)} Bank, S. B. \and I. Laine,
\textit {Representations of solutions of periodic second order linear differential equations},  J. Reine Angew. Math., \textbf{344} (1983), 1--21.

\bibitem{BLL1986}  Bank, S. B., I. Laine \and J. K. Langley,
\textit {On the frequency of zeros of solutions of second order linear differential equations}, Resultate Math., \textbf{10} (1986), 8--24.

\bibitem{BLL1989} Bank, S. B., I. Laine and J. K. Langley,
\textit{Oscillation results for solutions of linear differential equations in the complex domain},
Resultate Math., \textbf{16} (1989), 3--15.

\bibitem{Batola1982} Batola, F., \textit{Une g\'en\'eralisation d\'une formule Meixer-Tricomi}, Canad. J. Math. \textbf{34} (1982) 411-422.
.
\bibitem{BenDunn1996} Bender, C. M. \and G. V. Dunne, \textit{Quasi-exactly solvable systems and orthogonal polynomials}, J. Math. Phys. \textbf{37} (1) (1996) 6--11.

\bibitem{Boch1929} Bochner, S., \textit{\"Uber Sturm-Liouvillesche Polynomsysteme}, Math. Z. \textbf{29} (1929) 730--736.

\bibitem{BrowneSleeman1980} {Browne, P. J.} \and {B. D. Sleeman}, \textit{Applications of multiparameter spectral theory to special functions}, Math. Proc. Camb. Phil. Soc. (1980), \textbf{87}, 275--283.


\bibitem{BC31}  Burchnall, J. L. \and T. W. Chaundy, \textit{The commutative ordinary differential operators II -- The identity} $P^n=Q^m$ Proc. Royal Soc. London Ser. A \textbf{134} (1931), 471--485.

\bibitem{Chiang2000}  {Chiang, Y. M. }, \textit{On the {zero-free}
solutions of linear periodic differential equations in the complex
plane}, Results Math., \textbf{38} (2000), no. 3--4, 213--225.

\bibitem{CI2006}
{Chiang, Y. M. } \and {M. E. H. Ismail}, \textit{On value
distribution theory of second order periodic ODEs, special functions
and orthogonal polynomials}, {Canad. J. Math.} {\bf 58} (2006),
726--767 (Corrigendum {Canad. J. Math.} {\bf 62} (2010), 261).

\bibitem{CYu2011} {Chiang, Y. M.} \and {K. W. Yu} , \textit{Subnormal solutions of non-homogeneous periodic ODEs, special functions and special polynomials} J. Reine Angew. Math., \textbf{651} (2011), 127--164.

\bibitem{CYu2012} {Chiang, Y. M.} \and {K. W. Yu} , \textit{On complex oscillation, quantization of periodic non-homogeneous ODEs and special functions} Proc. Roy. Soc. Edinb., Sect. A Math. \textbf{142A} (2012) 449--477

\bibitem{Chihara1978} Chihara, T. S., \textrm{An Introduction to Orthogonal Polynomials}, Gordon and Breach, N. Y., Lond. and Pairs 1978

\bibitem{Darboux1882} Darboux, G., \textit{Sur une \'equation lin\'eaire}, C.R. t. XCIV, no. 25, pp. 1645--1648 (1882).

\bibitem{DMR1978_I} Decarreau, A.; Maroni, P.; Robert, A., \textit{Sur les \'equations confluentes de l'\'equation de Heun}, (French)  Ann. Soc. Sci. Bruxelles, S\'er. I  \textbf{92}  (1978), no. 3, 151--189.

\bibitem{DMR1978_II} Decarreau, A.; Dumont-Lepage, M.-Cl.,  Maroni, P.; Robert, A., Ronveaux, A.,\textit{Forms canoniques des \'equations confluentes de l\'equation de Heun}, (French)  Ann. Soc. Sci. Bruxelles, S\'er. I-II  \textbf{92}  (1978), no. 3, 53--78.

\bibitem{Dunham1932} Dunham, J. L., \textit{The energy levels of a rotating vibrator}, Phys. Rev. \textbf{41} (1932), 721--731.

\bibitem{BatemanI} Erd\'elyi, A, { W. Magnus}, {F. Oberhettinger} \and {\it F. G. Tricomi}, {Higher Transcendental Functions, Vol. I}, McGraw-Hill, New York, London and Toronto 1953.

\bibitem{Bateman} Erd\'elyi, A, { W. Magnus}, {F. Oberhettinger} \and {\it F. G. Tricomi}, {Higher Transcendental Functions, Vol. III}, McGraw-Hill, New York, London and Toronto 1953.

%\bibitem{GKO1993} Gonz\'alez-L\'opez, A., N. Kamran \and P. J. Olver, \textit{Normalizability of one-dimensional quasi-exactly solvable Schr\"odinger operators}, Commun. Math. Phys. \textbf{153} (1993) 117-146.

\bibitem{Erdelyi1942} Erd\'elyi, A. \textit{Integral equations for Heun functions}, Quat. J. Math., Oxford \textbf{13} (1942) 107--112.

\bibitem{FGR1996} Finkel, F., A. Gonz\'alez-L\'opez \and M. A. Rodr\'iguez, \textit{Quasi-exactly solvable potentials on the line and orthogonal polynomials}, J. Math. Phys., \textbf{37} (8) (1996) 3954--3972.

\bibitem{Forsyth1902} Forsyth, A. R., Theory of Differential Equations, Part III, Vol. IV, Cambridge Univ. Press 1902.

\bibitem{Gautschi1967} Gautschi, W., \textit{Computational aspects of three-term recurrence relations}, SIAM Rev. \textbf{9} (1967) 24--82.

\bibitem{GKO1993} Gonz\'alez-L\'opez, A., N. Kamran \and P. J. Olver, \textit{Normalizability of one-dimensional quasi-exactly solvable Schr?dinger operators}, Comm. Math. Phys. \textbf{153} (1993), no. 1, 117--146.

\bibitem{GKO1994} Gonz\'alez-L\'opez, A., N. Kamran \and P. J. Olver, \textit{Quasi-exact solvability}, in Lie  algebras, cohomology, and new applications to quantum mechanics (Springfield, MO, 1992), 113�C140, Contemp. Math. \textbf{160} Amer. Math. Soc., Providence, RI, 1994.

\bibitem{Gr1978} Grosswald, E., \textit{Bessel Polynomials}, Lecture Notes in Math., \textbf{698} Berlin, Springer-Verlag, 1978.

\bibitem {Hayman1975} Hayman, W. K., {Meromorphic Functions}, Clarendon Press, Oxford 1964 (reprinted in 1975 with appendix).

\bibitem{Hautot1969} Hautot, A., \textit{Sur les solutions polynomiales de \'le\'quation diff\'erentielle $zP^{\prime\prime}_n+(az^2+bz+c)P^\prime_n+(d+ez+fz^2)P_n=0$}, Bull. Soc. Sci. Li\'ege, \textbf{38} (1969), 660--663
\bibitem{Heun1889} Heun, K., \textit{Zur Theorie der Riemann'schen Functionen Zweiter Ordnung mit vier Verzweigungspunkten}
Math. Ann. \textbf{33} (1889), 161--179.

\bibitem{Ince1922} Ince, E. L., \textit{On the connection between linear differential systems and integral equations}, Proc. Roy. Soc. Edinb. \textbf{42} (1922), 43-53.

\bibitem{Ince1926} Ince, E. L., Ordinary Differential Equations, Longmans, Green and Co. 1926 (reprinted by Dover Publ. Inc., New York. 1956).

%\bibitem{KazSlav1996} Kazakov, A. Ya. \and Slavyanov, S. Yu., \textit{Integral equations for special functions of Heun class}, Meth. Appl. Anal., \textbf{4} (3) (1996) 447-456.

\bibitem{KF1949} Krall, H. L., \and O. Frink, \textit{A new class of orthogonal polynomials: the Bessel polynomials}, Trans. Amer. Math. Soc. \textbf{65} no. 1, (1949) 100--115.

%\bibitem{Kreuser1914} Kreuser, P. \textit{\"Uber das Verhalten der Integrate homogener linearer Differenzengleichungen im Unendlichen}, Dissertation, University of T\"ubingen, Leipzig, 1914.

\bibitem{La93} Laine, I., Nevanlinna Theory and Complex Differential
Equations, Walter de Gruyter, Berlin, New York 1993.

%\bibitem{LambeWard1934} Lambe, C. G., \and Ward, D. R., \textit{Some differential equations and associated integral equations}, Quart. J. Math. Oxford, \textbf{5} (1934) 81-97.

%\bibitem{LKL2002} Li, Kang \and Leong, Spheroidal Wave Functions in Electromagnetic Theory, ??? 2002.

%\bibitem{Lovitt1924} Lovitt, W. V., Linear Integral Equations, McGraw-Hill Book Comp., 1924 (republished by Dover in 1950)

\bibitem{Maroni1979} Maroni, P., \textit{Sur quelques relations int\'egrales entre les solutions de l'equation biconfluente de Heun}, Ann. Inst. Henri Poincar\'e, \textbf{30} (1979) 315--332.

\bibitem{Masson1982} Masson, D., \textit{The rotating harmonic oscillator eigenvalue problem. I. Continued fractions and analytic continuation}, J. Math. Phys. \textbf{24} (1983), no. 8, 2074--2088.

%\bibitem{MorseFesh} Morse, P. M. \and Feshbach, H. ``Methods of Theoretical Physics" Part I, New York, Toronto and London {McGraw-Hill} 1953.

\bibitem{Olver1997} Olver, F. W. J., \textrm{Asymptotics and Special Functions}, A. K. Peter, Natick, Massachusetts 1997.

\bibitem{OlverP1997} Olver, P.J., \textit{A quasi-exactly solvable travel guide}, in: GROUP21: Physical Applications and Mathematical Aspects of Geometry, Groups, and Algebras, vol. 1, H.-D. Doebner, W. Scherer, and P. Nattermann, eds., World Scientific, Singapore, 1997, pp. 285--295.

\bibitem{Perron1911} Perron, O. \textit{\"Uber einen Satz des Herrn Poincar\'e} Acta. Math., \textbf{34} (1911) 109--137.

%\bibitem{PortStir1990} Porter, D. \and Stirling, D. S. G., \textrm{Integral equations. A practical treatment, from spectral theory to applications}, Camb. Texts in Applied Mathematics. Camb. Univ. Press, Cambridge, 1990.

\bibitem{Rovder1974} Rovder, J., \textit{Zeros of the polynomial solutions of the differential equation $xy^{\prime\prime}+(\beta_0+\beta_1x+\beta_2x^2)y^\prime+(\gamma-n\beta_2x)y=0$}, Mat. \v cas., \textbf{24}, (1974), 15--20.

\bibitem{Ronv1995} Ronveaux, A., Heun's Differential Equations, Oxford, N.Y. and Tokyo, Oxford University Press, 1995.


\bibitem{SlavLay2000} Slavyanov, S. Yu., \and W. Lay, Special Functions: A unified Theory Based on Singularities, Oxford Univ. Press, 2000.

\bibitem{Slep1978} Slepian, D. \textit{Prolate spheroidal wave functions, Fourier analysis, and uncertainty. V - The discrete case} The Bell System Technical J. (1978).

\bibitem{Sleeman1978} Sleeman, B. D., \textrm{Multiparameter spectral theory in Hilbert space},  London, Pitman Press, 1978

\bibitem{Smithies1958} Smithies, F., ``Integral Equations", Cambridge Tract Math., \textbf{49}, 1958, Camb. Univ. Press

\bibitem{Sz1975} Szeg\H{o}, G., \textit{Orthogonal Polynomials} (4th ed.), Amer. Math. Soc. Coll. Pulb. Vol. 23, 1975 (first edition in 1939).


\bibitem{Turb1988} Turbiner, A. V., \textit{Quasi-exactly-solvable problems and $sl(2)$ algebra}, Commun. Math. Phys. \textbf{118} (1988) 467--474.

\bibitem{Turb1996} Turbiner, A. V., \textit{Quasi-exactly-solvable differential equations}, CRC Handbook of Lie Group Analysis of Differential Equations, Vol. 3, N. H. Ibragimov, ed., CRC Press, Boca Raton, Fl., 1996, pp. 329--364.

\bibitem{Urwin1975} Urwin, K. M., \textit{Note on certain polynomials solutions of Heun's equation and two confluent forms}, Bull. Soc. Roy. Sci. Li\'ege, \textbf{44} (1975), no. 3--4, 169--173.

\bibitem{Watson:1944} Watson, G. N.,  {A Treatise on the Theory of Bessel Functions}, Camb. Univ. Press, 1944.

\bibitem{Whitt1912} Whittaker, E. T., \textit{On the functions associated with the elliptic cylinder in harmonic analysis}, Proc. Inter. Cong., Camb., Modern Analysis, 1912, 366--371.

\bibitem{Whitt1914} Whittaker, E. T., \textit{On a class of differential equations whose solutions satisfy integral equations}, Proc. Edinb. Math. Soc., \textbf{33} (1914), 14--23.

\bibitem{Whitt1915} Whittaker, E. T., \textit{On Lam\'e's differential equation and ellipsoidal harmonics}, Proc. London Math. Soc.(2) \textbf{14} (1915), 260--268

\bibitem{Wimp1984} Wimp, J., Computation with recurrence relations, Appl. Math. Ser., Pitman, Boston, MA, 1984

\bibitem{WW27} Whittaker, E. T. \and Watson, G. N. , {A Course of Modern Analysis} (4th ed.), Camb. Univ. Press 1927, reprinted in 1992.

%\bibitem {A} T. Aoki, \textit{Calcul exponentiel des op\'erateurs microdifferentiels d'ordre infini.} I, Ann. Inst. Fourier (Grenoble) \textbf{33} (1983), 227--250.

%\bibitem {B} R. Brown, \textit{On a conjecture of Dirichlet}, Amer. Math. Soc., Providence, RI, 1993.

%\bibitem {D} R. A. DeVore, \textit{Approximation of functions}, Proc. Sympos. Appl. Math., vol. 36, Amer. Math. Soc., Providence, RI, 1986, pp. 34--56.

\end{thebibliography}

\end{document}

%------------------------------------------------------------------------------
% End of journal.tex
%------------------------------------------------------------------------------

1310.5507.extracted_bib
\begin{thebibliography}{10}

\bibitem{AbrSte1972} Abramowitz, M. \and I. Stegun, \textrm{Handbook of Mathematical Functions}, 9th edition, Dover 1972.

\bibitem{AAR1999}  Andrews, G. E., R. Askey \and R. Roy, \textit{Special Functions}, Encyclopedia of Mathematics and its Applications Vol. \textbf{71}, Camb. Univ. Press, Cambridge 1999.

\bibitem{Arscott1964} Arscott, F. M. \textrm{Periodic Differential Equations: An Introduction to Mathieu, Lam\'e and Allied Functions}, Pergamon Press, 1964.

\bibitem{Bank1993} Bank, S. B., \textit{On the explicit determination of certain solutions of periodic differential equations}, Complex Variables Theory Appl. \textbf{23} (1993), no. 1--2, 101--121.

\bibitem{BL1982} Bank, S. B. \and I. Laine,
\textit{On the oscillation theory of {$f\sp{\prime\prime}+{A}f=0$}
where {${A}$} is entire}, Trans. Amer. Math. Soc., \textbf{273}
(1982), 351--363.

\bibitem{BL1983(1)} Bank, S. B. \and I. Laine,
\textit {Representations of solutions of periodic second order linear differential equations},  J. Reine Angew. Math., \textbf{344} (1983), 1--21.

\bibitem{BLL1986}  Bank, S. B., I. Laine \and J. K. Langley,
\textit {On the frequency of zeros of solutions of second order linear differential equations}, Resultate Math., \textbf{10} (1986), 8--24.

\bibitem{BLL1989} Bank, S. B., I. Laine and J. K. Langley,
\textit{Oscillation results for solutions of linear differential equations in the complex domain},
Resultate Math., \textbf{16} (1989), 3--15.

\bibitem{Batola1982} Batola, F., \textit{Une g\'en\'eralisation d\'une formule Meixer-Tricomi}, Canad. J. Math. \textbf{34} (1982) 411-422.
.
\bibitem{BenDunn1996} Bender, C. M. \and G. V. Dunne, \textit{Quasi-exactly solvable systems and orthogonal polynomials}, J. Math. Phys. \textbf{37} (1) (1996) 6--11.

\bibitem{Boch1929} Bochner, S., \textit{\"Uber Sturm-Liouvillesche Polynomsysteme}, Math. Z. \textbf{29} (1929) 730--736.

\bibitem{BrowneSleeman1980} {Browne, P. J.} \and {B. D. Sleeman}, \textit{Applications of multiparameter spectral theory to special functions}, Math. Proc. Camb. Phil. Soc. (1980), \textbf{87}, 275--283.


\bibitem{BC31}  Burchnall, J. L. \and T. W. Chaundy, \textit{The commutative ordinary differential operators II -- The identity} $P^n=Q^m$ Proc. Royal Soc. London Ser. A \textbf{134} (1931), 471--485.

\bibitem{Chiang2000}  {Chiang, Y. M. }, \textit{On the {zero-free}
solutions of linear periodic differential equations in the complex
plane}, Results Math., \textbf{38} (2000), no. 3--4, 213--225.

\bibitem{CI2006}
{Chiang, Y. M. } \and {M. E. H. Ismail}, \textit{On value
distribution theory of second order periodic ODEs, special functions
and orthogonal polynomials}, {Canad. J. Math.} {\bf 58} (2006),
726--767 (Corrigendum {Canad. J. Math.} {\bf 62} (2010), 261).

\bibitem{CYu2011} {Chiang, Y. M.} \and {K. W. Yu} , \textit{Subnormal solutions of non-homogeneous periodic ODEs, special functions and special polynomials} J. Reine Angew. Math., \textbf{651} (2011), 127--164.

\bibitem{CYu2012} {Chiang, Y. M.} \and {K. W. Yu} , \textit{On complex oscillation, quantization of periodic non-homogeneous ODEs and special functions} Proc. Roy. Soc. Edinb., Sect. A Math. \textbf{142A} (2012) 449--477

\bibitem{Chihara1978} Chihara, T. S., \textrm{An Introduction to Orthogonal Polynomials}, Gordon and Breach, N. Y., Lond. and Pairs 1978

\bibitem{Darboux1882} Darboux, G., \textit{Sur une \'equation lin\'eaire}, C.R. t. XCIV, no. 25, pp. 1645--1648 (1882).

\bibitem{DMR1978_I} Decarreau, A.; Maroni, P.; Robert, A., \textit{Sur les \'equations confluentes de l'\'equation de Heun}, (French)  Ann. Soc. Sci. Bruxelles, S\'er. I  \textbf{92}  (1978), no. 3, 151--189.

\bibitem{DMR1978_II} Decarreau, A.; Dumont-Lepage, M.-Cl.,  Maroni, P.; Robert, A., Ronveaux, A.,\textit{Forms canoniques des \'equations confluentes de l\'equation de Heun}, (French)  Ann. Soc. Sci. Bruxelles, S\'er. I-II  \textbf{92}  (1978), no. 3, 53--78.

\bibitem{Dunham1932} Dunham, J. L., \textit{The energy levels of a rotating vibrator}, Phys. Rev. \textbf{41} (1932), 721--731.

\bibitem{BatemanI} Erd\'elyi, A, { W. Magnus}, {F. Oberhettinger} \and {\it F. G. Tricomi}, {Higher Transcendental Functions, Vol. I}, McGraw-Hill, New York, London and Toronto 1953.

\bibitem{Bateman} Erd\'elyi, A, { W. Magnus}, {F. Oberhettinger} \and {\it F. G. Tricomi}, {Higher Transcendental Functions, Vol. III}, McGraw-Hill, New York, London and Toronto 1953.

%\bibitem{GKO1993} Gonz\'alez-L\'opez, A., N. Kamran \and P. J. Olver, \textit{Normalizability of one-dimensional quasi-exactly solvable Schr\"odinger operators}, Commun. Math. Phys. \textbf{153} (1993) 117-146.

\bibitem{Erdelyi1942} Erd\'elyi, A. \textit{Integral equations for Heun functions}, Quat. J. Math., Oxford \textbf{13} (1942) 107--112.

\bibitem{FGR1996} Finkel, F., A. Gonz\'alez-L\'opez \and M. A. Rodr\'iguez, \textit{Quasi-exactly solvable potentials on the line and orthogonal polynomials}, J. Math. Phys., \textbf{37} (8) (1996) 3954--3972.

\bibitem{Forsyth1902} Forsyth, A. R., Theory of Differential Equations, Part III, Vol. IV, Cambridge Univ. Press 1902.

\bibitem{Gautschi1967} Gautschi, W., \textit{Computational aspects of three-term recurrence relations}, SIAM Rev. \textbf{9} (1967) 24--82.

\bibitem{GKO1993} Gonz\'alez-L\'opez, A., N. Kamran \and P. J. Olver, \textit{Normalizability of one-dimensional quasi-exactly solvable Schr?dinger operators}, Comm. Math. Phys. \textbf{153} (1993), no. 1, 117--146.

\bibitem{GKO1994} Gonz\'alez-L\'opez, A., N. Kamran \and P. J. Olver, \textit{Quasi-exact solvability}, in Lie  algebras, cohomology, and new applications to quantum mechanics (Springfield, MO, 1992), 113�C140, Contemp. Math. \textbf{160} Amer. Math. Soc., Providence, RI, 1994.

\bibitem{Gr1978} Grosswald, E., \textit{Bessel Polynomials}, Lecture Notes in Math., \textbf{698} Berlin, Springer-Verlag, 1978.

\bibitem {Hayman1975} Hayman, W. K., {Meromorphic Functions}, Clarendon Press, Oxford 1964 (reprinted in 1975 with appendix).

\bibitem{Hautot1969} Hautot, A., \textit{Sur les solutions polynomiales de \'le\'quation diff\'erentielle $zP^{\prime\prime}_n+(az^2+bz+c)P^\prime_n+(d+ez+fz^2)P_n=0$}, Bull. Soc. Sci. Li\'ege, \textbf{38} (1969), 660--663
\bibitem{Heun1889} Heun, K., \textit{Zur Theorie der Riemann'schen Functionen Zweiter Ordnung mit vier Verzweigungspunkten}
Math. Ann. \textbf{33} (1889), 161--179.

\bibitem{Ince1922} Ince, E. L., \textit{On the connection between linear differential systems and integral equations}, Proc. Roy. Soc. Edinb. \textbf{42} (1922), 43-53.

\bibitem{Ince1926} Ince, E. L., Ordinary Differential Equations, Longmans, Green and Co. 1926 (reprinted by Dover Publ. Inc., New York. 1956).

%\bibitem{KazSlav1996} Kazakov, A. Ya. \and Slavyanov, S. Yu., \textit{Integral equations for special functions of Heun class}, Meth. Appl. Anal., \textbf{4} (3) (1996) 447-456.

\bibitem{KF1949} Krall, H. L., \and O. Frink, \textit{A new class of orthogonal polynomials: the Bessel polynomials}, Trans. Amer. Math. Soc. \textbf{65} no. 1, (1949) 100--115.

%\bibitem{Kreuser1914} Kreuser, P. \textit{\"Uber das Verhalten der Integrate homogener linearer Differenzengleichungen im Unendlichen}, Dissertation, University of T\"ubingen, Leipzig, 1914.

\bibitem{La93} Laine, I., Nevanlinna Theory and Complex Differential
Equations, Walter de Gruyter, Berlin, New York 1993.

%\bibitem{LambeWard1934} Lambe, C. G., \and Ward, D. R., \textit{Some differential equations and associated integral equations}, Quart. J. Math. Oxford, \textbf{5} (1934) 81-97.

%\bibitem{LKL2002} Li, Kang \and Leong, Spheroidal Wave Functions in Electromagnetic Theory, ??? 2002.

%\bibitem{Lovitt1924} Lovitt, W. V., Linear Integral Equations, McGraw-Hill Book Comp., 1924 (republished by Dover in 1950)

\bibitem{Maroni1979} Maroni, P., \textit{Sur quelques relations int\'egrales entre les solutions de l'equation biconfluente de Heun}, Ann. Inst. Henri Poincar\'e, \textbf{30} (1979) 315--332.

\bibitem{Masson1982} Masson, D., \textit{The rotating harmonic oscillator eigenvalue problem. I. Continued fractions and analytic continuation}, J. Math. Phys. \textbf{24} (1983), no. 8, 2074--2088.

%\bibitem{MorseFesh} Morse, P. M. \and Feshbach, H. ``Methods of Theoretical Physics" Part I, New York, Toronto and London {McGraw-Hill} 1953.

\bibitem{Olver1997} Olver, F. W. J., \textrm{Asymptotics and Special Functions}, A. K. Peter, Natick, Massachusetts 1997.

\bibitem{OlverP1997} Olver, P.J., \textit{A quasi-exactly solvable travel guide}, in: GROUP21: Physical Applications and Mathematical Aspects of Geometry, Groups, and Algebras, vol. 1, H.-D. Doebner, W. Scherer, and P. Nattermann, eds., World Scientific, Singapore, 1997, pp. 285--295.

\bibitem{Perron1911} Perron, O. \textit{\"Uber einen Satz des Herrn Poincar\'e} Acta. Math., \textbf{34} (1911) 109--137.

%\bibitem{PortStir1990} Porter, D. \and Stirling, D. S. G., \textrm{Integral equations. A practical treatment, from spectral theory to applications}, Camb. Texts in Applied Mathematics. Camb. Univ. Press, Cambridge, 1990.

\bibitem{Rovder1974} Rovder, J., \textit{Zeros of the polynomial solutions of the differential equation $xy^{\prime\prime}+(\beta_0+\beta_1x+\beta_2x^2)y^\prime+(\gamma-n\beta_2x)y=0$}, Mat. \v cas., \textbf{24}, (1974), 15--20.

\bibitem{Ronv1995} Ronveaux, A., Heun's Differential Equations, Oxford, N.Y. and Tokyo, Oxford University Press, 1995.


\bibitem{SlavLay2000} Slavyanov, S. Yu., \and W. Lay, Special Functions: A unified Theory Based on Singularities, Oxford Univ. Press, 2000.

\bibitem{Slep1978} Slepian, D. \textit{Prolate spheroidal wave functions, Fourier analysis, and uncertainty. V - The discrete case} The Bell System Technical J. (1978).

\bibitem{Sleeman1978} Sleeman, B. D., \textrm{Multiparameter spectral theory in Hilbert space},  London, Pitman Press, 1978

\bibitem{Smithies1958} Smithies, F., ``Integral Equations", Cambridge Tract Math., \textbf{49}, 1958, Camb. Univ. Press

\bibitem{Sz1975} Szeg\H{o}, G., \textit{Orthogonal Polynomials} (4th ed.), Amer. Math. Soc. Coll. Pulb. Vol. 23, 1975 (first edition in 1939).


\bibitem{Turb1988} Turbiner, A. V., \textit{Quasi-exactly-solvable problems and $sl(2)$ algebra}, Commun. Math. Phys. \textbf{118} (1988) 467--474.

\bibitem{Turb1996} Turbiner, A. V., \textit{Quasi-exactly-solvable differential equations}, CRC Handbook of Lie Group Analysis of Differential Equations, Vol. 3, N. H. Ibragimov, ed., CRC Press, Boca Raton, Fl., 1996, pp. 329--364.

\bibitem{Urwin1975} Urwin, K. M., \textit{Note on certain polynomials solutions of Heun's equation and two confluent forms}, Bull. Soc. Roy. Sci. Li\'ege, \textbf{44} (1975), no. 3--4, 169--173.

\bibitem{Watson:1944} Watson, G. N.,  {A Treatise on the Theory of Bessel Functions}, Camb. Univ. Press, 1944.

\bibitem{Whitt1912} Whittaker, E. T., \textit{On the functions associated with the elliptic cylinder in harmonic analysis}, Proc. Inter. Cong., Camb., Modern Analysis, 1912, 366--371.

\bibitem{Whitt1914} Whittaker, E. T., \textit{On a class of differential equations whose solutions satisfy integral equations}, Proc. Edinb. Math. Soc., \textbf{33} (1914), 14--23.

\bibitem{Whitt1915} Whittaker, E. T., \textit{On Lam\'e's differential equation and ellipsoidal harmonics}, Proc. London Math. Soc.(2) \textbf{14} (1915), 260--268

\bibitem{Wimp1984} Wimp, J., Computation with recurrence relations, Appl. Math. Ser., Pitman, Boston, MA, 1984

\bibitem{WW27} Whittaker, E. T. \and Watson, G. N. , {A Course of Modern Analysis} (4th ed.), Camb. Univ. Press 1927, reprinted in 1992.

%\bibitem {A} T. Aoki, \textit{Calcul exponentiel des op\'erateurs microdifferentiels d'ordre infini.} I, Ann. Inst. Fourier (Grenoble) \textbf{33} (1983), 227--250.

%\bibitem {B} R. Brown, \textit{On a conjecture of Dirichlet}, Amer. Math. Soc., Providence, RI, 1993.

%\bibitem {D} R. A. DeVore, \textit{Approximation of functions}, Proc. Sympos. Appl. Math., vol. 36, Amer. Math. Soc., Providence, RI, 1986, pp. 34--56.

\end{thebibliography}