Series Solutions of the Non-Stationary Heun Equation

Series Solutions of the Non-Stationary Heun Equation††This paper is a
contribution to the Special Issue on Elliptic Hypergeometric Functions and
Their Applications. The full collection is available at
https://www.emis.de/journals/SIGMA/EHF2017.html

Farrokh ATAI †‡ and Edwin LANGMANN †

F. Atai and E. Langmann

† Department of Physics, KTH Royal Institute of Technology, SE-10691
Stockholm, Sweden<EMAIL_ADDRESS>‡ Department of Mathematics, Kobe
University, Rokko, Kobe 657-8501, Japan

<EMAIL_ADDRESS>
Received October 10, 2017, in final form February 08, 2018; Published online
February 16, 2018

We consider the non-stationary Heun equation, also known as quantum Painlevé
VI, which has appeared in different works on quantum integrable models and
conformal field theory. We use a generalized kernel function identity to
transform the problem to solve this equation into a differential-difference
equation which, as we show, can be solved by efficient recursive algorithms.
We thus obtain series representations of solutions which provide elliptic
generalizations of the Jacobi polynomials. These series reproduce, in a
limiting case, a perturbative solution of the Heun equation due to Takemura,
but our method is different in that we expand in non-conventional basis
functions that allow us to obtain explicit formulas to all orders; in
particular, for special parameter values, our series reduce to a single term.

Heun equation; Lamé equation; Kernel functions; quantum Painlevé VI;
perturbation theory

33E20; 81Q05; 16R60

##  1 Introduction

This paper is devoted to the study of a certain non-stationary Schrödinger
equation with elliptic potential (see (1.2) below) which has appeared in
different works on quantum integrable models and conformal field theory [7,
14, 34]. This so-called non-stationary Heun equation [25], also known as
quantum Painlevé VI [27], is a generalization of a second-order Fuchsian
differential equation with four regular singular points known as Heun equation
[32, 37]. Another important special case is the non-stationary Lamé equation
[1, 5], which is also known as KZB heat equation111To be more precise: the
special case of the KZB equations corresponding to
$\mathfrak{g}=\mathfrak{sl}_{2}$ and $n=1$. [9]. Our main result is a method
to construct a solution of the non-stationary Heun equation which is complete
(in a sense explained below) and which is complementary to results obtained by
other methods in that it can be used for generic parameter values. We also
discuss various interesting special cases of our results. As will be
explained, this is part of a research program using kernel functions to solve
quantum models of Calogero–Moser–Sutherland (CMS) type; see [10, 19, 20, 22,
24].

###  1.1 Background

There is a long tradition in mathematical physics that regards integrable
models and the mathematical theory of special functions as two sides of the
same coin. While most textbook examples of such models can be solved using
results about special functions known since a long time (they can be found in
[41], for example), modern developments in conformal field theory and quantum
statistical physics have led to special functions which belong to function
classes which are the subject of ongoing research. As an outstanding example
we mention the representation theoretic approach to the solution of certain
hyperbolic differential equations due to Etingof, Frenkel and Kirillov [4, 5],
which is a beautiful generalization of classical works on the Gauss
hypergeometric equation pioneered by Gelfand and motivated by conformal field
theory. A complementary approach to this inspired by string theory was
developed in works by Felder, Varchenko, and others [6, 8, 9]. The non-
stationary Lamé equation is the simplest non-trivial example to which these
results apply. More recently this equation appeared in the construction of
exact 4-point correlation functions of the quantum Liouville model [7], and in
the exact solution of the eight-vertex model [1]. It was conjectured in [7]
that results about the non-stationary Lamé equation should have natural
generalizations to the non-stationary Heun equation. This was recently
confirmed in important examples by Rosengren [33, 34], who proved and
generalized conjectures in [1], and Kolb [14], who generalized the
representation theoretic approach to the non-stationary Lamé equation in [4,
5] to the Heun case. We also mention another important approach due to
Nekrasov and Shatiashvili allowing to construct functions in the Heun class
and which is based on supersymmetric gauge theories [29]; see also [17] for
recent related work. Another related subject is the AGT conjecture which has
led to explicit combinatorial expressions for conformal blocks related to the
non-stationary Heun equation; see, e.g., [28] and references therein.

The Heun equation has received considerable interest as a natural equation
defining a function class generalizing the Gauss hypergeometric functions;
see, e.g., [26, 32, 37, 38, 40] and references therein. The examples mentioned
in the previous paragraph motivate us to extend this scope: our work is
intended as a contribution towards a general theory of special functions
defined by the non-stationary Heun equation. The approach we use is based on a
kernel function method developed in a series of papers in order to construct
exact eigenfunctions of quantum models of CMS type; see [10, 19, 20, 22, 24].
The relation of this to other approaches to the Heun equation based on kernel
functions [3, 13, 18, 31, 35] is discussed in Section 7.

To explain our method we recall that a kernel function for a pair of
Schrödinger operators $H(x)$ and $\tilde{H}(y)$ is an explicitly known
function $K(x,y)$ satisfying the functional identity
$(H(x)-\tilde{H}(y)-c)K(x,y)=0$ for some constant $c$; in the examples of
interest to us the Schrödinger operators are Hamiltonians defining a CMS-type
model which can have equal ($H=\tilde{H}$) or different ($\tilde{H}\neq H$)
parameters, and we write $H(x)$ to indicate that this differential operator
acts on functions depending on the variable $x$. A basic observation
underlying the kernel function method we use is that CMS-type Hamiltonians
have eigenfunctions which are easy to construct but uninteresting for
applications, and that kernel functions provide a tool to transform such
uninteresting eigenfunctions to interesting ones [22]. It was shown in [10]
that all classical CMS-type models possess kernel functions which allow the
construction of particular series representations of standard eigenfunctions
in this way (by classical CMS-type models we mean those whose eigenfunctions
provide natural many-variable generalizations of classical orthogonal
polynomials, including the corresponding deformed CMS models in the sense of
Chalykh, Feigin, Veselov and Sergeev [2, 36]). We also mention recent related
work of Hallnäs and Ruijsenaars [11] constructing eigenfunctions of the CMS
model with $1/\sinh^{2}$-interactions using such a kernel function and which
goes beyond the paradigm of polynomial eigenfunctions. It was found that the
elliptic generalizations of the CMS-models, which can be regarded as many-
variable generalizations of the Lamé and Heun equations, possess kernel
functions [21, 23, 25] generalizing those in [10], and in [19] one such kernel
function was used to construct series representations of eigenfunctions of the
elliptic quantum Calogero–Sutherland (eCS) model [24]. We mention in passing a
series solution of the eCS model by Komori and Takemura [16] which, while
using the same expansion parameter, is different in important details from the
result in [19, 24] (this difference is analogous to the difference between
perturbative results of Takemura on the Heun equation in [39] and results in
the present paper, as discussed in the last paragraph of Section 1.2 below).

It is known that the kernel functions allowing for elliptic generalizations
are restricted by a so-called balancing condition $\kappa=0$ with $\kappa$ a
constant depending on the model parameters [15, 21]. If this condition is not
fulfilled one can often find a generalized kernel function $K(x,y,\tau)$
satisfying

$\displaystyle\left(\frac{{\rm
i}}{\pi}\kappa\frac{\partial}{\partial\tau}+H(x,\tau)-\tilde{H}(y,\tau)-c(\tau)\right)K(x,y,\tau)=0$
(1.1)

with $\tau$ the half-period ratio of the elliptic functions appearing in the
CMS Hamiltonians $H$ and $\tilde{H}$ [21, 23, 25] (see, e.g., Lemma 4.1 where
the $\tau$-dependence of $H$, $\tilde{H}$, $c$ and $K$ is suppressed). In the
present paper we use a generalized kernel function found in [25] to solve the
non-stationary Heun equation, following the approach in [24]. The basic
observation underlying our approach is that one can use the generalized kernel
function in (1.1) to transform an eigenfunction of the differential operator
$\frac{{\rm i}}{\pi}\kappa\frac{\partial}{\partial\tau}+\tilde{H}(y,\tau)$ to
an eigenfunction of $\frac{{\rm
i}}{\pi}\kappa\frac{\partial}{\partial\tau}+H(x,\tau)$ [25].

###  1.2 Summary of results

The non-stationary Heun equation can be written as222Here and in the following
we often suppress the dependence of functions on the variable $\tau$.

$\displaystyle\left(\frac{{\rm
i}}{\pi}\kappa\frac{\partial}{\partial\tau}-\frac{\partial^{2}}{\partial
x^{2}}+\sum_{\nu=0}^{3}g_{\nu}(g_{\nu}-1)\wp(x+\omega_{\nu})\right)\psi(x)=E\psi(x)$
(1.2)

with $\wp(x)$ the Weierstrass elliptic function with periods $(2\pi,2\pi\tau$)
and

$\displaystyle\omega_{0}=0,\qquad\omega_{1}=\pi,\qquad\omega_{2}=-\pi-\pi\tau,\qquad\omega_{3}=\pi\tau$
(1.3)

(for the convenience of the reader we collect the definitions of $\wp$ and
other well-known special functions we need in Appendix A). To simplify
notation we set $\omega_{1}=\pi$ here and in most parts of this paper; see
Appendix B for how to transform our results to other values of $\omega_{1}$.
The parameters $g_{0}$, $g_{1}$, $g_{2}$, $g_{3}$ and $\kappa$ can be
arbitrary complex numbers for our general results.333We have to make some
restrictions on parameters due to a technical problem referred to as
resonances but, as discussed in Section 2.4, many of these restrictions are
irrelevant in practice. Our aim is to construct functions
$\psi(x)\equiv\psi(x,\tau;\\{g_{\nu}\\}_{\nu=0}^{3},\kappa)$ of two complex
variables $x$ and $\tau$, $\Im(\tau)>0$, that satisfy this differential
equation for some $E\equiv E(\tau;\\{g_{\nu}\\}_{\nu=0}^{3},\kappa)$; a more
precise characterization of our solutions is given in (1.5)–(1.6) below. It is
important to note that, for $\kappa\neq 0$, $E$ can be transformed to 0, or
any other convenient value, by changing the normalization of $\psi(x)$ (this
follows from the obvious invariance of (1.2) under the transformation

$\displaystyle\psi(x)\to C\psi(x),\qquad E\to E+\frac{{\rm
i}}{\pi}\kappa\frac{1}{C}\frac{\partial C}{\partial\tau}$ (1.4)

for arbitrary analytic functions $C$ of $\tau$). However, we sometimes find it
convenient to impose a normalization condition on $\psi(x)$ so that $E$
remains significant even for $\kappa\neq 0$. Important special cases of (1.2)
include the Heun equation ($\kappa=0$), the non-stationary Lamé equation
($g_{\nu}=g$ independent of $\nu$, or $g_{\nu}(g_{\nu}-1)=0$ for three of the
$\nu$’s), and the Lamé equation (both specializations). Many of our results
are new even for these special cases (to our knowledge).

To explain the nature of our solutions we recall that, in the trigonometric
limit $\Im(\tau)\to+\infty$, (1.2) reduces to the stationary Schrödinger
equation with the Pöschl–Teller potential $\propto
g_{0}(g_{0}-1)/\sin^{2}(x/2)+g_{1}(g_{1}-1)/\cos^{2}(x/2)$ which has
explicitly known solutions equal to Jacobi polynomials up to a common factor
$\sin(x/2)^{g_{0}}\cos(x/2)^{g_{1}}$ (see (2.3)–(2.4) for details). The
solutions of (1.2) that we construct are a generalization of this: they are of
the form

$\displaystyle\psi_{n}(x)=(2q^{1/4})^{-g_{0}-g_{1}}\left(\prod_{\nu=0}^{3}\theta_{\nu+1}\big{(}\tfrac{1}{2}x\big{)}^{g_{\nu}}\right)\mathcal{P}_{n}\left(\cos(x)\right),$
$\displaystyle
E_{n}=\kappa^{2}\left(\frac{1}{12}-\frac{\eta_{1}}{\pi}\right)-\sum_{\nu=0}^{3}g_{\nu}(g_{\nu}-1)\frac{\eta_{1}}{\pi}+\mathcal{E}_{n}$
(1.5)

with $\theta_{\nu+1}(z)$ the Jacobi theta functions, $q=\exp({\rm i}\pi\tau)$
the nomé, $\eta_{1}/\pi$ in (A.5) and

$\displaystyle\mathcal{P}_{n}(z)=\sum_{\ell=0}^{\infty}\mathcal{P}^{(\ell)}_{n}(z)q^{\ell},\qquad\mathcal{E}_{n}=\sum_{\ell=0}^{\infty}\mathcal{E}^{(\ell)}_{n}q^{\ell},$
(1.6)

with

$\displaystyle\mathcal{P}^{(0)}_{n}(z)=P_{n}^{\big{(}g_{0}-\frac{1}{2},g_{1}-\frac{1}{2}\big{)}}(z),\qquad\mathcal{E}^{(0)}_{n}=\left(n+\frac{g_{0}+g_{1}}{2}\right)^{2},\qquad
n\in\mathbb{N}_{0},$ (1.7)

and $P^{\big{(}g_{0}-\frac{1}{2},g_{1}-\frac{1}{2}\big{)}}_{n}(z)$ the Jacobi
polynomials (see (A.11)). Our main result provides efficient recursive
procedures to compute the functions $\mathcal{P}^{(\ell)}_{n}(z)$ which, as we
show, are polynomials of degree $n+\ell$ in $z$; see Propositions 5.1 and 5.2
for two complementary variants of this result. Thus one can regard
$\mathcal{P}_{n}(z)$ in (1.6) as an elliptic generalization of Jacobi
polynomials. It is interesting to note that these elliptic generalizations
exist even for negative integers $n$ if $q$ is non-zero, but they vanish like
$O(q^{-n})$ for $n<0$ as $q\to 0$. We conjecture that the series in (1.6) is
absolutely convergent and converges to a $L^{2}$-function on $[0,\pi]$ for
$|q|\leq q_{0}$ and some $q_{0}>0$ depending on parameters (this is known to
be true in the Heun case $\kappa=0$ from work by Takemura [39]; see also [24]
for a convergence proof for the Lamé case which, as we believe, can be
generalized). However, this question is left for future work: for simplicity
we treat series like in (1.6) as formal power series.

Our solution (1.5)–(1.6) of (1.2) is complete in the sense that the
$\psi_{n}(x)$ provide a complete orthonormal basis in the Hilbert space of
$L^{2}$-functions on $[0,\pi]$ for $g_{0}+g_{1}>0$ in the trigonometric case
$q=0$ (we believe that this is true even for $q>0$). Moreover, we give an
efficient recursive procedure to compute the coefficients
$\mathcal{P}^{(\ell)}_{n}(z)$ and $\mathcal{E}^{(\ell)}_{n}$ of the power
series in (1.6). We note that, in the Heun case $\kappa=0$, the $E_{n}$
correspond to the eigenvalues of the $BC_{1}$ elliptic CMS Hamiltonian
discovered by Inozemtsev [12]. We thus refer to the $E_{n}$ as generalized
eigenvalues in the following.

It is important to note that, by exploiting the invariance of (1.2) under the
transformation in (1.4), we obtain two complementary variants of results: in
the first variant we impose a normalization conditions on $\psi_{n}(x)$ such
that the generalized eigenvalues $E_{n}$ are significant (see Proposition 5.1
and Theorem 6.1), and in the second we fix $E_{n}$ by a convenient condition
(see Proposition 5.2 and Theorem 6.3). These two variants of results are
complementary to each other in that the second is somewhat simpler but
restricted to $\kappa\neq 0$, whereas the first applies to the case $\kappa=0$
as well. However, as will be discussed after Proposition 5.2, the second
variant of results implies an interesting representation of the eigenvalues
$E_{n}$ in the limit $\kappa\to 0$.

One important feature of our method is that it provides particular
$q$-dependent basis functions $\\{f_{m}(z)\\}_{m\in\mathbb{Z}}$ to expand the
functions $\mathcal{P}_{n}(z)$ in; see (2.6) and (3.1). This is useful since
these functions $f_{m}(z)$ take into account much of the complexity of the
problem; for example, in special cases the expansion coefficients are trivial,
and in these cases our method gives explicit integral representations of the
solutions $\mathcal{P}_{n}(z)$ (see Section 3.2). For general parameter
values, we obtain a system of equations for these expansion coefficients
which, in the trigonometric case $q=0$, can be solved by diagonalizing a
triangular matrix and which, for non-zero $q$, can be solved by efficient
perturbative algorithms; see Propositions 5.1 and 5.2. A perturbative solution
of these algorithms to all orders in $q$ is obtained in Section 6; see
Theorems 6.1 and 6.3. We note that results for the elliptic
Calogero–Sutherland model corresponding to the ones in Sections 4.2, 5 and 6
were obtained by one of us in [19, 20], and [24], respectively. We also
mention that, in the special case $\kappa=0$, we obtain results for the Heun
equation in Section 6.1 which differ from the ones by Takemura who used Jacobi
polynomials as basis to expand the functions $\mathcal{P}_{n}(z)$ [39]. As
already mentioned, this difference is analogous to the difference between the
perturbative results for the eCS model in [19] and the one by Komori and
Takemura in [16]: in the latter work the eigenfunctions are expanded in Jack
polynomials, whereas in the former an unconventional basis is used which
allows for an explicit solution to all orders [24].

###  1.3 Plan

Section 2 contains preliminary material: a summary of notation (Section 2.1),
a review of a well-known solution of (1.2) for $\Im(\tau)\to\infty$ in terms
of Jacobi polynomials (Section 2.2), the definition and properties of our
basis functions $f_{m}(z)$ (Section 2.3), and a discussion of a technicality
referred to as resonances (Section 2.4). In Section 3 we present our key
result, which is a transformation of the problem to solve (1.2) into a
differential-difference equation (Proposition 3.1), together with a discussion
of special cases (Section 3.2); the proof of this key result is given in
Section 4. In Section 5 we present two complementary recursive algorithms to
solve this differential-difference equation, and Section 6 contains the
corresponding explicit solutions to all orders. We conclude with final remarks
in Section 7. Five appendices contain definitions and properties of special
functions we use (Appendix A), details on how to translate our results for
$\omega_{1}=\pi$ to other values of $\omega_{1}$ (Appendix B), explicit
results from one of our recursive algorithms at low order (Appendix C),
derivations of results needed in proofs (Appendix D), and a short discussion
of the combinatorial structure of our solution (Appendix E).

##  2 Preliminaries

We collect definitions and preliminary results that we use.

###  2.1 Notation

We use the special functions ($\xi$ and $q$ are complex variables; $|q|<1$)

$\displaystyle\Theta_{1}(\xi)\equiv(1-\xi)\prod_{n=1}^{\infty}\big{(}1-q^{2n}\xi\big{)}\big{(}1-q^{2n}\xi^{-1}\big{)},\qquad\Theta_{2}(\xi)\equiv\Theta_{1}(-\xi),$
$\displaystyle\Theta_{3}(\xi)\equiv\prod_{n=1}^{\infty}\big{(}1+q^{2n-1}\xi\big{)}\big{(}1+q^{2n-1}\xi^{-1}\big{)},\qquad\Theta_{4}(\xi)\equiv\Theta_{3}(-\xi)$
(2.1)

and

$\displaystyle\Theta(z,\xi)\equiv\big{(}1-2z\xi+\xi^{2}\big{)}\prod_{n=1}^{\infty}\big{(}1-2q^{2n}\xi
z+q^{4n}\xi^{2}\big{)}\big{(}1-2q^{2n}\xi^{-1}z+q^{4n}\xi^{-2}\big{)},$ (2.2)

which all are closely related to the Jacobi theta functions (see (A.7) and
(A.8)).

We denote as $\mathbb{N}_{0}$ and $\mathbb{Z}^{\prime}$ the sets of non-
negative and non-zero integers, respectively. The symbol $\delta(m,n)$ for
integers $m$, $n$ denotes the Kronecker delta.

###  2.2 Trigonometric limit

The non-stationary Heun equation simplifies in the trigonometric case $q=0$ to

$\displaystyle\left(-\frac{\partial^{2}}{\partial
x^{2}}+\frac{g_{0}(g_{0}-1)}{4\sin^{2}\frac{1}{2}x}+\frac{g_{1}(g_{1}-1)}{4\cos^{2}\frac{1}{2}x}\right)\psi^{(0)}(x)=E^{(0)}\psi^{(0)}(x),$
(2.3)

which is known to have solutions

$\displaystyle\psi_{n}^{(0)}(x)=\big{(}\sin\tfrac{1}{2}x\big{)}^{g_{0}}\big{(}\cos\tfrac{1}{2}x\big{)}^{g_{1}}P^{\big{(}g_{0}-\frac{1}{2},g_{1}-\frac{1}{2}\big{)}}_{n}(\cos
x),\qquad E^{(0)}_{n}=\left(n+\frac{g_{0}+g_{1}}{2}\right)^{2}$ (2.4)

with the Jacobi polynomials $P_{n}^{(\alpha,\beta)}(z)$ in (A.11) (see, e.g.,
[30, Table 18.8.1, 2nd line]). We will use a well-known uniqueness result
about this solution in the following form.

######  Lemma 2.1.

Let $\psi_{n}^{(0)}(x)$ be a solution of (2.3) of the form

$\displaystyle\psi_{n}^{(0)}(x)=\big{(}\sin\tfrac{1}{2}x\big{)}^{g_{0}}\big{(}\cos\tfrac{1}{2}x\big{)}^{g_{1}}P(\cos
x)$

for some constant $E^{(0)}$, with $P(z)$ a polynomial of degree $n$ such that

$\displaystyle
P(z)=\frac{(n+g_{0}+g_{1})_{n}}{2^{n}n!}z^{n}+O\big{(}z^{n-1}\big{)}$ (2.5)

for some $n\in\mathbb{N}_{0}$. Then
$P(z)=P^{\big{(}g_{0}-\frac{1}{2},g_{1}-\frac{1}{2}\big{)}}_{n}(z)$ and
$E^{(0)}=E^{(0)}_{n}$.

The proof is standard and therefore omitted. (Note that the normalization in
(2.5) follows from (A.12).)

We will use this result to fix the normalization of our solutions so as to get
Jacobi polynomials in the trigonometric case $q=0$.

###  2.3 Basis functions

As mentioned in the introduction, one important feature of our method is that
we use non-trivial basis functions $f_{m}(z)$. These functions are defined by
the following generating function,

$\displaystyle\frac{\prod\limits_{\nu=0}^{3}\Theta_{\nu+1}(\xi)^{\tilde{g}_{\nu}}}{\Theta(z,\xi)^{\lambda}}\equiv\sum_{m\in\mathbb{Z}}f_{m}(z)\xi^{m},\qquad|q|<|\xi|<1,$
(2.6) $\displaystyle\tilde{g}_{\nu}\equiv\lambda-
g_{\nu},\qquad\lambda\equiv\tfrac{1}{2}(g_{0}+g_{1}+g_{2}+g_{3}-\kappa)$ (2.7)

with the special functions $\Theta_{\nu+1}(\xi)$, $\Theta(z,\xi)$ defined in
(2.1)–(2.2). It is easy to check that the series on the r.h.s. in (2.6) is
absolutely convergent in the region indicated, and thus the functions
$f_{m}(z)$ are well-defined.444One can show that, for fixed complex parameters
$\\{g_{\nu}\\}_{\nu=0}^{3}$ and $\kappa$, and all $m\in\mathbb{Z}$, $f_{m}(z)$
is analytic for $|q|<1$ and $z$ in some $q$-dependent open domain which
includes the interval $[-1,1]$. Since we restrict ourselves to results in the
sense of formal power series in $q$, we only need the following
characterization of these functions (the proof of this is technical and thus
deferred to an appendix).

######  Lemma 2.2.

The functions $f_{m}(z)$ defined in (2.6)–(2.7) have the following power
series expansion

$\displaystyle f_{m}(z)=\sum_{\ell=0}^{\infty}f_{m}^{(\ell)}(z)q^{\ell}$ (2.8)

with $f_{m}^{(\ell)}(z)=0$ for $m+\ell<0$ and $f_{m}^{(\ell)}(z)$ a polynomial
of degree $m+\ell$ in $z$ for $m+\ell\geq 0$. In particular,

$\displaystyle
f_{m}^{(0)}(z)=\binom{-\lambda}{m}(-2z)^{m}+O\big{(}z^{m-1}\big{)}$ (2.9)

with the binomial coefficient $\binom{-\lambda}{m}$ as usual $($see (A.10)$)$.

###### Proof.

See Appendix D.1. ∎

We will also use the following integral representation of these functions:

$\displaystyle f_{m}(z)=\oint_{\mathcal{C}}\frac{d\xi}{2\pi{\rm
i}\xi}\xi^{-m}\frac{\prod\limits_{\nu=0}^{3}\Theta_{\nu+1}(\xi)^{\tilde{g}_{\nu}}}{\Theta(z,\xi)^{\lambda}}$
(2.10)

with $\mathcal{C}$ the contour once around the circle with radius $|\xi|=R$,
$|q|<R<1$, taken counterclockwise (this is equivalent to (2.6) by Cauchy’s
theorem).

######  Remark 2.3.

It is clear from (2.9) that the functions $f^{(0)}_{m}(z)$ are non-trivial
polynomials in $z$ of degree $m$ for all $m\in\mathbb{N}_{0}$ only if
$-\lambda\notin\mathbb{N}_{0}$ (i.e., only in this case the binomial
coefficient on the r.h.s. in (2.9) is always non-zero). If
$-\lambda=k\in\mathbb{N}_{0}$, one can see from (2.10) that $f^{(0)}_{m}(z)$
is a polynomials of degree $\leq k$ for all $m\in\mathbb{N}_{0}$. Thus, to get
complete results, we sometimes impose the condition
$-\lambda\notin\mathbb{N}_{0}$.

###  2.4 Resonances

We discuss a technical issue encountered in Sections 5 and 6: to prove
Propositions 5.1 and 5.2 and Theorems 6.1 and 6.3 we find it convenient to
impose the following no-resonance conditions: either $\Im(\kappa)\neq 0$ and
$g_{0}+g_{1}\in\mathbb{R}$, or $\kappa=0$ and $g_{0}+g_{1}\notin\mathbb{Z}$.
At first sight this seems to exclude many cases of interests in applications
but, at closer inspection, one finds that this is not the case: as explained
in this section, our results can be used even in cases where these conditions
fail.

We first explain the reason for these conditions: our solutions are obtained
by an unconventional variant of perturbation theory leading to series
solutions which are linear combinations of products of the following
generalized energy difference fractions,

$\displaystyle\frac{1}{E^{(0)}_{m}-E^{(0)}_{n}-\kappa\ell}=\frac{1}{(m-n)(m+n+g_{0}+g_{1})-\kappa\ell}$
(2.11)

with $E^{(0)}_{n}$ in (2.4) the energy eigenvalues of the unperturbed problem
$q=0$ and $\ell=0,1,2,\ldots$; we refer to a case $(m,\ell)$ where, for fixed
$n$, the denominator in (2.11) is zero as resonance. The reason for the
conditions on parameters above is that they are a simple means to rule out
resonances, and this guarantees that our series are well-defined. However,
while these conditions are sufficient, they are not necessary: one peculiar
feature of our perturbation theory is that the series we obtain have
singularities coming from energy difference denominators as in (2.11), but
many of these singularities are removable. Thus our results can be extended to
cases where our no-resonance conditions fail. We now give a precise
formulation for one important such case.

######  Lemma 2.4.

In the Heun case $\kappa=0$ and for fixed $n\in\mathbb{N}_{0}$, the results
for the expansion coefficients $\mathcal{P}_{n}^{(\ell)}(z)$ and
$\mathcal{E}_{n}^{(\ell)}$ in Proposition 5.1 and Theorem 6.1 are valid even
in the limit $g_{0}+g_{1}\to k\in\mathbb{Z}$ with $k>-(2n+1)$.

(The proof is given at the end of this section.)

Thus, in the Heun case $\kappa=0$ and for all $n\in\mathbb{N}_{0}$, our
results can be extended to the case $g_{0}+g_{1}\in\mathbb{N}_{0}$ of interest
in applications. We believe that, in a similar manner, our results can be
extended to interesting cases with non-zero real $\kappa$. As will become
clear in our proof of Lemma 2.4 below, the challenge to make precise and prove
such a result is that the generalization of standard perturbation theory [39]
to $\kappa\neq 0$ is not known (to our knowledge).

###### Proof of Lemma 2.4.

For $\kappa=0$ and fixed $n\in\mathbb{N}_{0}$, the energy denominator in
(2.11) appearing in standard perturbation theory [39] are only for
$m\in\mathbb{N}_{0}$ different from $n$, and thus all pertinent fractions in
(2.11) are finite if $g_{0}+g_{1}>-(2n+1)$. In our perturbative expansions we
encounter fractions in (2.11) for arbitrary integers $m\neq n$. However, it is
clear that our results for the coefficients $\mathcal{P}_{n}^{(\ell)}(z)$ and
$\mathcal{E}_{n}^{(\ell)}$ of the perturbative solution defined in (1.6) must
be identical with the corresponding results obtained in standard perturbation
theory. Thus the singularities coming from resonance fractions in our
perturbation expansion must cancel: our perturbative results remain finite
even in the limit when $g_{0}+g_{1}$ becomes integer. ∎

######  Remark 2.5.

The resonance problem that we encounter in this paper is very similar to the
one which appeared in the treatment of the eCS model in [23, 24]. This is no
coincidence: results in the special case $N=2$ in op. cit. correspond to ours
in the Lamé case $g_{0}=g_{1}=g_{2}=g_{3}$, $\kappa=0$. The interested reader
is referred to [24] for a more extensive discussion of resonances.

##  3 Key result and special cases

In Section 3.1 we present our key result, which is a transformation of the
problem to solve the non-stationary Heun equation in (1.2) with the ansatz in
(1.5) to a problem to solve a differential-difference equation; as we show in
subsequent sections, the latter problem allows for efficient solutions. In
Section 3.2 we point out special non-trivial cases where our key result leads
to explicit integral representations of solutions of (1.2).

###  3.1 Differential-difference equation

We construct solutions $\psi_{n}(x)$, $E_{n}$ of the non-stationary Heun
equation in (1.2) of the form (1.5)–(1.7) in the sense of formal power series
in $q$. One important feature of our method is that we expand

$\displaystyle\mathcal{P}_{n}(z)={\mathcal{N}}_{n}\sum_{m\in\mathbb{Z}}\alpha_{n}(m)f_{m}(z)$
(3.1)

with non-trivial basis functions $f_{m}(z)$ given in (2.6)–(2.7) and
characterized in Lemma 2.2. As will be shown, the following constant ensures
the normalization in (1.7),

$\displaystyle{\mathcal{N}}_{n}=\frac{(n+g_{0}+g_{1})_{n}}{4^{n}(\lambda)_{n}}$
(3.2)

with the raising Pochhammer symbol $(x)_{n}$ in (A.9) and $\lambda$ in (2.7);
note that ${\mathcal{N}}_{n}$ is finite and non-zero for all integers $n$ if
$-\lambda\notin\mathbb{N}_{0}$ for $n>0$ and
$-(g_{0}+g_{1})\notin\mathbb{N}_{0}$ for $n<0$ (see the discussion after
(A.9)).

Our key result is equations determining $\alpha_{n}(m)$ and $\mathcal{E}_{n}$
and which, as we will show, can be solved efficiently. To state this result we
introduce the convenient shorthand notation

$\displaystyle\gamma_{k}^{\mu}\equiv\begin{cases}\tilde{g}_{0}(\tilde{g}_{0}-1)+(-1)^{\mu}\tilde{g}_{1}(\tilde{g}_{1}-1)&\text{if
}\frac{1}{2}k\in\mathbb{N}_{0},\\\
(-1)^{\mu}\tilde{g}_{2}(\tilde{g}_{2}-1)+\tilde{g}_{3}(\tilde{g}_{3}-1)&\text{if
}\frac{1}{2}(k-1)\in\mathbb{N}_{0}\end{cases}$ (3.3)

for $(\mu,k)\in\mathbb{Z}\times\mathbb{N}_{0}$ (recall that
$\tilde{g}_{0}=\frac{1}{2}(-g_{0}+g_{1}+g_{2}+g_{3}-\kappa)$,
$\tilde{g}_{1}=\frac{1}{2}(g_{0}-g_{1}+g_{2}+g_{3}-\kappa)$ etc.; note that
$\gamma_{k}^{\mu}=\gamma_{k+2r}^{\mu+2s}$ for all integers $r$, $s$).

######  Proposition 3.1.

Let $n\in\mathbb{Z}$, $-\lambda\notin\mathbb{N}_{0}$ for $n>0$ and
$-(g_{0}+g_{1})\notin\mathbb{N}_{0}$ for $n<0$, and assume that
$\mathcal{E}_{n}$ and $\alpha_{n}(m)$ for $m\in\mathbb{Z}$ satisfy the
following system of equations,

$\displaystyle\left[-\kappa q\frac{\partial}{\partial
q}+E^{(0)}_{m}-\mathcal{E}_{n}\right]\alpha_{n}(m)$
$\displaystyle\qquad{}=\sum_{\mu=1}^{\infty}\mu\gamma_{0}^{\mu}\alpha_{n}(m+\mu)+\sum_{\mu=1}^{\infty}\frac{\mu
q^{\mu}}{1-q^{2\mu}}\big{(}\gamma_{0}^{\mu}q^{\mu}+\gamma_{1}^{\mu}\big{)}[\alpha_{n}(m+\mu)+\alpha_{n}(m-\mu)]$
(3.4)

and the condition

$\displaystyle\alpha_{n}(m)|_{q=0}=\begin{cases}0,&m>n,\\\ 1,&m=n.\end{cases}$
(3.5)

Then $\psi_{n}(x)$, $E_{n}$ in (1.5)–(1.6) and (3.1)–(3.2) satisfy the non-
stationary Heun equation in (1.2), and the conditions in (1.7) hold true
provided $-(g_{0}+g_{1})\notin\mathbb{N}$.

(The proof is given in Section 4.)

It is important to note that the conditions above do not determine
$\mathcal{P}_{n}(z)$ and $\mathcal{E}_{n}$ uniquely: (3.4)–(3.5) are invariant
under

$\displaystyle\alpha_{n}(m)\to
C\alpha_{n}(m),\qquad\mathcal{E}_{n}\to\mathcal{E}_{n}-\kappa\frac{1}{C}q\frac{\partial}{\partial
q}C$ (3.6)

for any change of normalization $C=1+O(q)$ analytic in $q$ (this is a
consequence of the invariance of (1.2) under (1.4)). This ambiguity can be
fixed for generic parameter values by replacing (3.5) by a stronger condition;
see Propositions 5.1 and 5.2 for two different ways to do this.

This result also provides simple explicit solutions of (1.2) for special
particular parameter values, as elaborated in Section 3.2.

######  Remark 3.2.

It interesting to note that the solutions $\alpha_{n}(m)$ and
$\mathcal{E}_{n}$ of the equations in Proposition 3.1 are essentially
independent of $n$ in the following sense: they are of the form

$\displaystyle\alpha_{n}(m)=a(m-n),\qquad\mathcal{E}_{n}=(P/2)^{2}+\tilde{\mathcal{E}}$

with functions $a(k)$ and $\tilde{\mathcal{E}}$ depending on $n$ only in the
combination

$\displaystyle P\equiv 2n+g_{0}+g_{1}$ (3.7)

(this is easy to check). This and the notation introduced here are useful in
computations and in the presentation of results; we use this in (5.8)–(5.9)
and in Appendix C.

###  3.2 Explicit solutions by integrals

The solutions $\alpha_{n}(m)$, $\mathcal{E}_{n}$ of (3.4)–(3.5) are
complicated in general. However, there exist non-trivial cases where
Proposition 3.1 gives simple explicit solutions of the non-stationary Heun
equation:

######  Corollary 3.3.

Let $n\in\mathbb{Z}$, $\lambda$ a complex parameter such that
$-\lambda\notin\mathbb{N}_{0}$ for $n>0$ and
$-(g_{0}+g_{1})\notin\mathbb{N}_{0}$ for $n<0$,

$\displaystyle\tilde{g}_{\nu}\in\\{0,1\\},\qquad
g_{\nu}=\lambda-\tilde{g}_{\nu},\qquad\nu=0,1,2,3,\qquad\kappa=2\lambda-\sum_{\nu=0}^{3}\tilde{g}_{\nu},$
(3.8)

${\mathcal{N}}_{n}$ in (3.2), and $\mathcal{C}$ the integration contour
defined after (2.10). Then the non-stationary Heun equation in (1.2) has
solutions $\psi_{n}(x)$, $E_{n}$ as in (1.5) with

$\displaystyle\mathcal{P}_{n}(z)={\mathcal{N}}_{n}\oint_{\mathcal{C}}\frac{d\xi}{2\pi{\rm
i}\xi}\xi^{-n}\frac{\prod\limits_{\nu=0}^{3}\Theta_{\nu+1}(\xi)^{\tilde{g}_{\nu}}}{\Theta(z,\xi)^{\lambda}},\qquad\mathcal{E}_{n}=\left(n+\frac{g_{0}+g_{1}}{2}\right)^{2}$
(3.9)

and such that the conditions in (1.7) hold true provided
$-(g_{0}+g_{1})\notin\mathbb{N}$.

###### Proof.

The assumptions imply $\gamma_{k}^{\mu}=0$ for all $k$, $\mu$, and the
equations in (3.4)–(3.5) in this case have the solution
$\alpha_{n}(m)=\delta(m,n)$, $\mathcal{E}_{n}=E^{(0)}_{n}$. Proposition 3.1
implies the result. ∎

Note that (3.8) gives several different one-parameter families
$(\\{g_{\nu}\\}_{\nu=0}^{3},\kappa)$, depending on $\lambda$, where this
result provides simple integral representations of elliptic Jacobi polynomials
$\mathcal{P}_{n}(z)$. Moreover, these formulas are non-trivial even in the
trigonometric case $q=0$: the results above imply the following integral
representations of Jacobi polynomials,

$\displaystyle
P_{n}^{\big{(}g-\frac{1}{2},g-\frac{1}{2}\big{)}}(z)=\frac{(n+2g)_{n}}{4^{n}(g)_{n}}\oint_{\mathcal{C}}\frac{d\xi}{2\pi{\rm
i}\xi}\xi^{-n}\frac{1}{(1-2z\xi+\xi^{2})^{g}},$ $\displaystyle
P_{n}^{\big{(}g-\frac{1}{2},g+\frac{1}{2}\big{)}}(z)=\frac{(n+2g+1)_{n}}{4^{n}(g+1)_{n}}\oint_{\mathcal{C}}\frac{d\xi}{2\pi{\rm
i}\xi}\xi^{-n}\frac{(1-\xi)}{(1-2z\xi+\xi^{2})^{g+1}},$ $\displaystyle
P_{n}^{\big{(}g+\frac{1}{2},g-\frac{1}{2}\big{)}}(z)=\frac{(n+2g+1)_{n}}{4^{n}(g+1)_{n}}\oint_{\mathcal{C}}\frac{d\xi}{2\pi{\rm
i}\xi}\xi^{-n}\frac{(1+\xi)}{(1-2z\xi+\xi^{2})^{g+1}},$ $\displaystyle
P_{n}^{\big{(}g-\frac{1}{2},g-\frac{1}{2}\big{)}}(z)=\frac{(n+2g)_{n}}{4^{n}(g+1)_{n}}\oint_{\mathcal{C}}\frac{d\xi}{2\pi{\rm
i}\xi}\xi^{-n}\frac{(1-\xi^{2})}{(1-2z\xi+\xi^{2})^{g+1}}$ (3.10)

for $n\in\mathbb{N}_{0}$ (the latter identities are obtained from Corollary
3.3 for the cases where $(\tilde{g}_{0},\tilde{g}_{1},\lambda)$ is $(0,0,g)$,
$(1,0,g+1)$, $(0,1,g+1)$, and $(1,1,g+1)$, respectively). The first identity
in (3.10) is equivalent to a well-known generating function for the Gegenbauer
polynomials (see (A.13)), and also the others can be found in [30].

######  Remark 3.4.

Fateev et al. gave integral representations of solutions of the non-stationary
Lamé equation [7] which, in a special case, are similar to the one above for
$g_{0}=g_{1}=g_{2}=g_{3}$. More specifically, the solution given in equation
(3.11) of [7] can be proved by a simple variant of the argument that we used
in order to obtain our solution in (3.9) (note that our parameter $\kappa$
corresponds to $-2/b^{2}$ in [7]). We also mention similar integral
representations of solutions of the non-stationary Lamé equation appearing in
works by Etingof and Kirillov (see [5, Theorem 5.1]) and Felder and Varchenko
(see, e.g., [9, Example 1.2]).

We emphasize that the result in Corollary 3.3 is more general in that it
includes some non-stationary Heun cases that cannot be reduced to a non-
stationary Lamé case.

######  Remark 3.5.

Corollary 3.3 can be obtained as special case $(N,\tilde{N})=(1,0)$ from
Proposition 4.1 in [25]. However, this is not easy to see, and it is therefore
worthwhile to emphasize this result here.

##  4 Proof of key result

We turn to the proof of Proposition 3.1. In Section 4.1 we derive the key
identity using the kernel function method. In Section 4.2 we consider the
trigonometric case $q=0$ to prove that the conditions in (3.5) and the
normalization condition in (3.2) yield a solution satisfying (1.7).

###  4.1 Kernel function method

We introduce the notation

$\displaystyle
H\big{(}x;\\{g_{\nu}\\}_{\nu=0}^{3}\big{)}\equiv-\frac{\partial^{2}}{\partial
x^{2}}+\sum_{\nu=0}^{3}g_{\nu}(g_{\nu}-1)\wp(x+\omega_{\nu})$ (4.1)

with $\omega_{\nu}$ in (1.3). This allows us to write the non-stationary Heun
equation in (1.2) as

$\displaystyle\left(\frac{{\rm
i}}{\pi}\kappa\frac{\partial}{\partial\tau}+H\big{(}x;\\{g_{\nu}\\}_{\nu=0}^{3}\big{)}-E\right)\psi(x)=0.$
(4.2)

We also recall the definitions of $\tilde{g}_{\nu}$ and $\lambda$ in (2.7).
Note that $H$ in (4.1) is the Hamiltonian defining the $BC_{1}$ Inozemtsev
model [39].

We obtain our result from the following generalized kernel function identity:

######  Lemma 4.1.

The function

$\displaystyle\mathcal{K}(x,y)\equiv\frac{\prod\limits_{\nu=0}^{3}\theta_{\nu+1}\big{(}\frac{1}{2}x\big{)}^{g_{\nu}}\theta_{\nu+1}\big{(}\frac{1}{2}y\big{)}^{\tilde{g}_{\nu}}}{\theta_{1}\big{(}\frac{1}{2}(x+y)\big{)}^{\lambda}\theta_{1}\big{(}\frac{1}{2}(x-y)\big{)}^{\lambda}}$
(4.3)

obeys the identity

$\displaystyle\left(\frac{{\rm
i}}{\pi}\kappa\frac{\partial}{\partial\tau}+H\big{(}x;\\{g_{\nu}\\}_{\nu=0}^{3}\big{)}-H\big{(}y;\\{\tilde{g}_{\nu}\\}_{\nu=0}^{3}\big{)}-C_{1,1}\right)\mathcal{K}(x,y)=0$
(4.4)

with

$\displaystyle C_{1,1}=2\kappa(1-\lambda)\frac{\eta_{1}}{\pi}.$ (4.5)

###### Proof.

This is the special case $N=M=1$ of Corollary 3.2 in [25] (note that the
symbols $\beta$ and $A_{1,1}$ in [25] correspond to $-2\pi{\rm i}\tau$ and
$2\kappa$ here, respectively). ∎

We use this to compute the action of $\frac{{\rm
i}}{\pi}\kappa\frac{\partial}{\partial\tau}+H(x;\\{g_{\nu}\\}_{\nu=0}^{3})$ on
the functions

$\displaystyle
F_{m}(x)\equiv\big{(}2q^{1/4}\big{)}^{-(g_{0}+g_{1})}\left(\prod_{\nu=0}^{3}\theta_{\nu+1}\big{(}\tfrac{1}{2}x\big{)}^{g_{\nu}}\right)f_{m}(\cos
x),\qquad m\in\mathbb{Z}$ (4.6)

with $f_{m}(z)$ defined in (2.6). For this we note that

$\displaystyle\mathcal{K}(x,y)=2^{g_{0}+g_{1}}{\rm e}^{{\rm
i}\pi\tilde{g}_{0}/2}G^{-\kappa}\big{(}2q^{1/4}\big{)}^{-(g_{0}+g_{1})}\left(\prod_{\nu=0}^{3}\theta_{\nu+1}\big{(}\tfrac{1}{2}x\big{)}^{g_{\nu}}\right)\frac{\prod\limits_{\nu=0}^{3}\Theta_{\nu+1}({\rm
e}^{{\rm i}y})^{\tilde{g}_{\nu}}}{\Theta(\cos x,{\rm e}^{{\rm
i}y})^{\lambda}}{\rm e}^{{\rm i}y(g_{0}+g_{1})/2}$

with $G$ defined in (A.6) (we used (A.7) and (A.8)). This and (2.6) show that
$\mathcal{K}(x,y)$ is a generating function for the functions in (4.6):

$\displaystyle\mathcal{K}(x,y)=2^{g_{0}+g_{1}}{\rm e}^{{\rm
i}\pi\tilde{g}_{0}/2}G^{-\kappa}\sum_{m\in\mathbb{Z}}F_{m}(x){\rm e}^{{\rm
i}\big{(}m+\frac{1}{2}(g_{0}+g_{1})\big{)}y},\qquad 0<\Im(y)<\pi\Im(\tau).$
(4.7)

To evaluate $H(y;\\{\tilde{g}_{\nu}\\}_{\nu=0}^{3})\mathcal{K}(x,y)$ we use
the following expansions

$\displaystyle\wp(y+\omega_{\nu})=-\frac{\eta_{1}}{\pi}-\sum_{\mu\in{\mathbb{Z}^{\prime}}}(S_{\nu})_{\mu}{\rm
e}^{{\rm i}\mu y},\qquad 0<\Im(y)<\pi\Im(\tau)$ (4.8)

with

$\displaystyle(S_{0})_{\mu}=\mu\frac{1}{1-q^{2\mu}}=|\mu|\frac{q^{|\mu|-\mu}}{1-q^{2|\mu|}},$
$\displaystyle(S_{1})_{\mu}=(-1)^{\mu}\mu\frac{1}{1-q^{2\mu}}=(-1)^{\mu}|\mu|\frac{q^{|\mu|-\mu}}{1-q^{2|\mu|}},$
$\displaystyle(S_{2})_{\mu}=(-1)^{\mu}\mu\frac{q^{\mu}}{1-q^{2\mu}}=(-1)^{\mu}|\mu|\frac{q^{|\mu|}}{1-q^{2|\mu|}},$
$\displaystyle(S_{3})_{\mu}=\mu\frac{q^{\mu}}{1-q^{2\mu}}=|\mu|\frac{q^{|\mu|}}{1-q^{2|\mu|}}$
(4.9)

(see Appendix D.2 for derivations of these formulas). From this the following
result is obtained by straightforward computations.

######  Lemma 4.2.

The functions in (4.6) satisfy

$\displaystyle\left(\frac{{\rm
i}}{\pi}\kappa\frac{\partial}{\partial\tau}+H\big{(}x;\\{g_{\nu}\\}_{\nu=0}^{3}\big{)}\right)F_{n}(x)=\big{(}C_{0}+E^{(0)}_{n}\big{)}F_{n}(x)-\sum_{\mu\in{\mathbb{Z}^{\prime}}}S_{\mu}F_{n-\mu}(x)$
(4.10)

with $E^{(0)}_{n}$ in (2.4),

$\displaystyle
S_{\mu}\equiv\sum_{\nu=0}^{3}\gamma_{\nu}(S_{\nu})_{\mu},\qquad\gamma_{\nu}\equiv\tilde{g}_{\nu}(\tilde{g}_{\nu}-1),$
(4.11)

$\tilde{g}_{\nu}\equiv\lambda-g_{\nu}$, and

$\displaystyle
C_{0}=\kappa^{2}\left(\frac{1}{12}-\frac{\eta_{1}}{\pi}\right)-\sum_{\nu=0}^{3}g_{\nu}(g_{\nu}-1)\frac{\eta_{1}}{\pi}.$
(4.12)

(The proof can be found at the end of this section.)

To proceed we make the ansatz

$\displaystyle\psi(x)=\mathcal{N}\sum_{m\in\mathbb{Z}}\alpha(m)F_{m}(x)$
(4.13)

with $\mathcal{N}$ a $q$-independent normalization constant to be determined.
Inserting this in (1.2) and using Lemma 4.2 we obtain

$\displaystyle\sum_{m\in\mathbb{Z}}\left(\left(\frac{{\rm
i}}{\pi}\kappa\frac{\partial}{\partial\tau}+E^{(0)}_{m}-\mathcal{E}\right)\alpha(m)-\sum_{\mu\in{\mathbb{Z}^{\prime}}}S_{\mu}\alpha(m+\mu)\right)F_{m}(x)=0$

with $\mathcal{E}\equiv E-C_{0}$. It follows that the function in (4.13) is a
solution of (1.2) if the coefficients $\alpha(m)$ and $\mathcal{E}$ satisfy

$\displaystyle\left(\frac{{\rm
i}}{\pi}\kappa\frac{\partial}{\partial\tau}+E^{(0)}_{m}-\mathcal{E}\right)\alpha(m)=\sum_{\mu\in{\mathbb{Z}^{\prime}}}S_{\mu}\alpha(m+\mu).$
(4.14)

Inserting (4.9) and (4.11) and changing variables from $\tau$ to $q={\rm
e}^{{\rm i}\pi\tau}$ allow us to write this as

$\displaystyle\left(-\kappa q\frac{\partial}{\partial
q}+E_{m}^{(0)}-\mathcal{E}\right)\alpha(m)=\sum_{\mu=1}^{\infty}\mu(\gamma_{0}+(-1)^{\mu}\gamma_{1})\alpha(m+\mu)$
$\displaystyle\qquad{}+\sum_{\mu=1}^{\infty}\mu\left(\frac{[\gamma_{0}+(-1)^{\mu}\gamma_{1})]q^{2\mu}}{1-q^{2\mu}}+\frac{[(-1)^{\mu}\gamma_{2}+\gamma_{3}]q^{\mu}}{1-q^{2\mu}}\right)[\alpha(m+\mu)+\alpha(m-\mu)]$
(4.15)

equivalent to (3.4). This proves that $\psi(x)$ in (4.13) and
$E=\mathcal{E}+C_{0}$ solve (1.2) provided (4.15) is fulfilled.

We are left to show that a solution $\alpha(m)=\alpha_{n}(m)$,
$\mathcal{E}=\mathcal{E}_{n}$ satisfying the condition in (3.5) is such that
(1.7) holds true. This is done in Section 4.2.

###### Proof of Lemma 4.2.

It follows from (4.4) that

$\displaystyle\left(\frac{{\rm
i}}{\pi}\kappa\frac{\partial}{\partial\tau}+H\big{(}x;\\{g_{\nu}\\}_{\nu=0}^{3}\big{)}\right)G^{\kappa}\mathcal{K}(x,y)=\big{(}H\big{(}y;\\{\tilde{g}_{\nu}\\}_{\nu=0}^{3}\big{)}+C_{1,1}^{\prime}\big{)}G^{\kappa}\mathcal{K}(x,y)$
(4.16)

with

$\displaystyle C_{1,1}^{\prime}\equiv
C_{1,1}+\kappa^{2}\left(\frac{1}{12}-\frac{\eta_{1}}{\pi}\right)$ (4.17)

(we computed

$\displaystyle\frac{{\rm
i}}{\pi}\frac{1}{G}\frac{\partial}{\partial\tau}G=\sum_{n=1}^{\infty}\frac{2nq^{2n}}{1-q^{2n}}=\sum_{n=1}^{\infty}\sum_{\ell=1}^{\infty}2nq^{2n\ell}=\sum_{\ell=1}^{\infty}\frac{2q^{2\ell}}{(1-q^{2\ell})^{2}}=\frac{1}{12}-\frac{\eta_{1}}{\pi}$

using (A.6) and (A.5)). To compute the r.h.s. in (4.16) we use

$\displaystyle
H\big{(}y;\\{\tilde{g}_{\nu}\\}_{\nu=0}^{3}\big{)}=-\frac{\partial^{2}}{\partial
y^{2}}-\sum_{\mu\in{\mathbb{Z}^{\prime}}}S_{\mu}{\rm e}^{{\rm i}\mu
y}-\sum_{\nu=0}^{3}\gamma_{\nu}\frac{\eta_{1}}{\pi},\qquad
0<\Im(y)<\pi\Im(\tau)$

obtained by inserting (4.8) in (4.1) and using (4.11). Using (4.7) and
equating terms with the same factor ${\rm e}^{{\rm
i}(n+\frac{1}{2}(g_{0}+g_{1}))y}$ give (4.10) and

$\displaystyle
C_{0}=C_{1,1}^{\prime}-\sum_{\nu=0}^{3}\gamma_{\nu}\frac{\eta_{1}}{\pi}.$

Using (4.5) and (4.17), inserting (4.11), and recalling (2.7), one obtains the
formula in (4.12). ∎

###  4.2 Trigonometric limit

For $q=0$ the equations in (3.4)–(3.5) simplify to

$\displaystyle\big{[}E^{(0)}_{m}-\mathcal{E}^{(0)}_{n}\big{]}\alpha^{(0)}_{n}(m)=\sum_{\mu=1}^{n-m}\mu\gamma_{0}^{\mu}\alpha^{(0)}_{n}(m+\mu)$
(4.18)

and $\alpha^{(0)}_{n}(n)=1$; we use the superscript “$(0)$” to indicate that a
quantity is for $q=0$. This system of equations is an eigenvalue problem for a
non-degenerate triangular matrix which can be solved by recursion:
$\mathcal{E}_{n}^{(0)}=E^{(0)}_{n}$, $\alpha^{(0)}_{n}(n)=1$, and

$\displaystyle\alpha_{n}^{(0)}(m<n)=\frac{1}{b^{(0)}_{n}(m-n)}\sum_{\mu=1}^{n-m}\mu\gamma_{0}^{\mu}\alpha^{(0)}_{n}(m+\mu)$
(4.19)

with

$\displaystyle b_{n}^{(0)}(m)\equiv
E^{(0)}_{m+n}-E^{(0)}_{n}=m(m+2n+g_{0}+g_{1}),$ (4.20)

provided that $g_{0}+g_{1}$ is such that

$\displaystyle b_{n}^{(0)}(-m)\neq 0\qquad\forall\,m=1,2,\ldots,n.$ (4.21)

The latter condition is implied by our assumption that
$-(g_{0}+g_{1})\notin\mathbb{N}$.

We recall that $f^{(0)}_{m}(z)\equiv\left.f_{m}(z)\right|_{q=0}$ is zero for
$m<0$ and a polynomial satisfying (2.9) for $m\geq 0$ (see Lemma 2.2), and
thus

$\displaystyle\mathcal{P}^{(0)}_{n}(z)\equiv\left.\mathcal{P}_{n}(z)\right|_{q=0}={\mathcal{N}}_{n}\sum_{m=0}^{n}\alpha_{n}^{(0)}(m)f_{m}^{(0)}(z)={\mathcal{N}}_{n}\binom{-\lambda}{n}(-2z)^{n}+O\big{(}z^{n-1}\big{)}$
(4.22)

is a polynomial of degree $n$. Moreover, the special case $q=0$ of results
proved above implies that
$\psi_{n}^{(0)}(x)\equiv(\sin\frac{1}{2}x)^{g_{0}}(\cos\frac{1}{2}x)^{g_{1}}\mathcal{P}_{n}^{(0)}(\cos
x)$ is a solution of (2.3) with $E^{(0)}=E^{(0)}_{n}$. Thus, by Lemma 2.1,
(1.7) holds true provided that the coefficient of the leading term in (4.22)
agrees with the one in (2.5). This fixes the normalization constant
$\mathcal{N}_{n}$ as in (3.2) and completes the proof of Proposition 3.1.

##  5 Recursive algorithms

We now present algorithms to compute the expansion coefficients
$\mathcal{P}_{n}^{(\ell)}(z)$ and $\mathcal{E}_{n}^{(\ell)}$ of our solution
defined in (1.6), which are based on Proposition 3.1. The first algorithm
given in Section 5.1 is such that it can be used for all values of $\kappa$,
including $\kappa=0$. We then give a second algorithm, which is a variant of
the first, which is simpler but requires $\kappa\neq 0$; see Section 5.2.

###  5.1 Algorithm I

We compute the functions $\mathcal{P}_{n}^{(\ell)}(z)$ and
$\mathcal{E}_{n}^{(\ell)}$ defined in (1.6) by solving the equations in
(3.4)–(3.5) with the ansatz

$\displaystyle\alpha_{n}(m)=\sum_{\ell=0}^{\infty}\alpha_{n}^{(\ell)}(m)q^{\ell},\qquad\mathcal{E}_{n}=\sum_{\ell=0}^{\infty}\mathcal{E}_{n}^{(\ell)}q^{\ell},$
(5.1)

which leads to recursive relations allowing a straightforward solution.
Inserting this solution in (3.1) and using Lemma 2.2 one obtains
representations of the functions $\mathcal{P}_{n}^{(\ell)}(z)$ which make
manifest that they are polynomials of degree $n+\ell$ in $z$ for $n+\ell\geq
0$ and zero otherwise.

To formulate this result we find it useful to introduce the shorthand notation

$\displaystyle b_{n}^{(\ell)}(k)\equiv
E^{(0)}_{n+k}-E^{(0)}_{n}-\kappa\ell=k(k+2n+g_{0}+g_{1})-\kappa\ell$ (5.2)

($E^{(0)}_{n}$ in (2.4)) and recall the definition of $\lambda$ in (2.7). Note
that the denominators in the fractions in (2.11) are equal to
$b_{n}^{(\ell)}(m-n)$.

######  Proposition 5.1.

Let $n\in\mathbb{Z}$, $-\lambda\notin\mathbb{N}_{0}$ for $n>0$ and
$-(g_{0}+g_{1})\notin\mathbb{N}_{0}$ for $n<0$,
$\\{f^{(\ell)}_{m}(z)\\}_{m\in\mathbb{Z}}$ the functions defined and
characterized in Lemma 2.2, and assume that either $\Im(\kappa)\neq 0$ and
$g_{0}+g_{1}\in\mathbb{R}$ or $\kappa=0$ and
$g_{0}+g_{1}\notin\mathbb{Z}$.555Note that the latter are the no-resonance
conditions discussed in Section 2.4. Then the non-stationary Heun equation in
(1.2) has a unique solution as in (1.5)–(1.6) given by

$\displaystyle\mathcal{P}_{n}^{(\ell)}(z)=\mathcal{N}_{n}\sum_{\ell^{\prime}=0}^{\ell}\sum_{m=-\ell^{\prime}}^{n+\ell-\ell^{\prime}}\alpha_{n}^{(\ell-\ell^{\prime})}(m)f_{m}^{(\ell^{\prime})}(z)$
(5.3)

with $\mathcal{N}_{n}$ in (3.2), and $\alpha_{n}^{(\ell)}(m)$,
$\mathcal{E}^{(\ell)}_{n}$ are determined by the following recursion
relations,

$\displaystyle\alpha_{n}^{(\ell)}(m)=\frac{1}{b_{n}^{(\ell)}(m-n)}\Biggl{(}\sum_{\ell^{\prime}=1}^{\ell}\mathcal{E}_{n}^{(\ell^{\prime})}\alpha_{n}^{(\ell-\ell^{\prime})}(m)+\sum_{\mu=1}^{n-m+\ell}\mu\gamma_{0}^{\mu}\alpha_{n}^{(\ell)}(m+\mu)$
$\displaystyle\hphantom{\alpha_{n}^{(\ell)}(m)=}{}+\sum_{\ell^{\prime}=0}^{\ell-1}\sum_{\mu=1}^{\ell-\ell^{\prime}}\sum_{k=1}^{\left\lfloor{\frac{\ell-\ell^{\prime}}{\mu}}\right\rfloor}\mu\gamma_{k}^{\mu}\delta_{\ell,\ell^{\prime}+k\mu}\big{[}\alpha_{n}^{(\ell^{\prime})}(m+\mu)+\alpha_{n}^{(\ell^{\prime})}(m-\mu)\big{]}\Biggr{)}$
(5.4)

for $m\neq n$, $\ell\geq 0$, and

$\displaystyle\mathcal{E}^{(\ell)}_{n}=-\sum_{\mu=1}^{\ell}\mu\gamma_{0}^{\mu}\alpha_{n}^{(\ell)}(n+\mu)-\sum_{\ell^{\prime}=0}^{\ell-1}\sum_{\mu=1}^{\ell-\ell^{\prime}}\sum_{k=1}^{\left\lfloor{\frac{\ell-\ell^{\prime}}{\mu}}\right\rfloor}\mu\gamma_{k}^{\mu}\delta_{\ell,\ell^{\prime}+k\mu}\big{[}\alpha_{n}^{(\ell^{\prime})}(n+\mu)+\alpha_{n}^{(\ell^{\prime})}(n-\mu)\big{]}$
(5.5)

for $\ell\geq 1$, together with the following conditions

$\displaystyle\alpha_{n}^{(0)}(n)=1,\qquad\alpha_{n}^{(\ell)}(n)=0\qquad\forall\,\ell\geq
1,\qquad\alpha_{n}^{(\ell)}(m)=0\qquad\forall\,m>n+\ell,\quad\ell\geq 0.$
(5.6)

The functions $\mathcal{P}_{n}^{(\ell)}(z)$ thus obtained are polynomials of
degree $n+\ell$ in $z$ for $n+\ell\geq 0$ and zero otherwise, and (1.7) holds
true provided $-(g_{0}+g_{1})\notin\mathbb{N}$.

(The proof can be found at the end of this section.)

The equations in (5.4)–(5.6) comprise the first of our perturbative solution
algorithms. It is important to note that it has a triangular structure and is
finite, which implies that it determines each coefficient
$\alpha_{n}^{(\ell)}(m)$ and $\mathcal{E}_{n}^{(\ell)}$ as a sum of finitely
many terms. To be specific: The set of pairs $(\ell,m)$ is partially ordered
as follows,

$\displaystyle(\ell^{\prime},m^{\prime})\prec(\ell,m)\Leftrightarrow(\ell^{\prime}<\ell)\qquad\text{or}\qquad(\ell^{\prime}=\ell\mbox{
and }m^{\prime}>m),$ (5.7)

and (5.4) and (5.5) have the form

$\displaystyle\alpha_{n}^{(\ell)}(m)=\sum_{(\ell^{\prime},m^{\prime})\prec(\ell,m)}B_{(\ell,m),(\ell^{\prime},m^{\prime})}\alpha_{n}^{(\ell^{\prime})}(m^{\prime}),$
$\displaystyle\mathcal{E}^{(\ell)}_{n}=\sum_{(\ell^{\prime},m^{\prime})\prec(\ell,n)}C_{(\ell,n),(\ell^{\prime},m^{\prime})}\alpha_{n}^{(\ell^{\prime})}(m^{\prime})$

with explicitly known coefficients $B_{(\ell,m),(\ell^{\prime},m^{\prime})}$
and $C_{(\ell,n),(\ell^{\prime},m^{\prime})}$ which, for fixed first argument,
are non-zero only for finitely many values of the second argument
$(\ell^{\prime},m^{\prime})$. For the convenience of the reader we compute the
solution $\alpha_{n}^{(\ell)}(m)$, $\mathcal{E}_{n}^{(\ell)}$ of this system
of equations analytically for $\ell=0,1,2$ in Appendix C, and we obtain the
following result for the generalized eigenvalues,
$\mathcal{E}_{n}=(P/2)^{2}+\mathcal{E}_{n}^{(1)}q+\mathcal{E}_{n}^{(2)}q^{2}+O(q^{3})$
with $P=2n+g_{0}+g_{1}$ and

$\displaystyle\mathcal{E}_{n}^{(1)}=\gamma_{0}^{1}\gamma_{1}^{1}\left(\frac{1}{P-1}-\frac{1}{P+1-\kappa}\right),$
(5.8)
$\displaystyle\mathcal{E}_{n}^{(2)}=\left(\gamma_{0}^{1}\right)^{2}\left(\frac{1}{P-1}-\frac{1}{P+1-2\kappa}\right)+\left(\gamma_{1}^{1}\right)^{2}\left(\frac{1}{P-1+\kappa}-\frac{1}{P+1-\kappa}\right)$
$\displaystyle\hphantom{\mathcal{E}_{n}^{(2)}=}{}+4\gamma_{0}^{0}\gamma_{1}^{0}\left(\frac{1}{2(P-2)}-\frac{1}{2(P+2)-2\kappa}\right)-2(\gamma_{0}^{1})^{2}\gamma_{1}^{0}\left(\frac{1}{[2(P-2)][P-1]}\right.$
$\displaystyle\left.\hphantom{\mathcal{E}_{n}^{(2)}=}{}-\frac{1}{[P-1][P+1-2\kappa]}+\frac{1}{[P+1-2\kappa][2(P+2)-2\kappa]}\right)$
$\displaystyle\hphantom{\mathcal{E}_{n}^{(2)}=}{}-2\gamma_{0}^{0}\left(\gamma_{1}^{1}\right)^{2}\left(\frac{1}{[2(P-2)][P-1+\kappa]}-\frac{1}{[P-1+\kappa][P+1-\kappa]}\right.$
$\displaystyle\left.\hphantom{\mathcal{E}_{n}^{(2)}=}{}+\frac{1}{[P+1-\kappa][2(P+2)-2\kappa]}\right)-(\gamma_{0}^{1})^{2}(\gamma_{1}^{1})^{2}\left(\frac{1}{[P-1+\kappa][P-1]^{2}}\right.$
$\displaystyle\hphantom{\mathcal{E}_{n}^{(2)}=}{}-\frac{1}{[P+1-2\kappa][P+1-\kappa]^{2}}-\frac{1}{[P-1+\kappa][P-1][P+1-\kappa]}$
$\displaystyle\hphantom{\mathcal{E}_{n}^{(2)}=}{}+\frac{1}{[P-1][P+1-\kappa][P+1-2\kappa]}-\frac{1}{[2(P-2)][P-1][P-1+\kappa]}$
$\displaystyle\left.\hphantom{\mathcal{E}_{n}^{(2)}=}{}+\frac{1}{[P+1-\kappa][P+1-2\kappa][2(P+2)-2\kappa]}\right).$
(5.9)

This computation can be straightforwardly implemented in a symbolic
programming language like MAPLE or MATHEMATICA and, in this way, extended to
higher values of $\ell$.666We computed $\mathcal{E}_{n}^{(\ell)}$ up to
$\ell=6$ using ordinary laptops and with computation times of the order of
minutes.

As proved by Ruijsenaars in [35], the Heun equation in (1.2) for $\kappa=0$
has solutions $\psi_{n}(x)$ and corresponding eigenvalues $E_{n}$ such that
the latter are invariant under all permutations of the following affine
combinations of coupling parameters,

$\displaystyle c_{0}\equiv g_{0}+g_{2}-1,\qquad c_{1}\equiv
g_{1}+g_{3}-1,\qquad c_{2}\equiv g_{1}-g_{3},\qquad c_{3}\equiv g_{0}-g_{2}.$
(5.10)

We used MAPLE to check that the coefficients $\mathcal{E}_{n}^{(\ell)}$,
$\ell=1,2,3,4,5$, determined by the algorithm in Proposition 5.2 have this
permutation symmetry for $\kappa=0$ (but not for $\kappa\neq 0$).

###### Proof of Proposition 5.1.

Insert the ansatz in (5.1) and the geometric series for $q^{\mu}/(1-q^{2\mu})$
into (3.4), compare terms which have the same power in $q$, and obtain

$\displaystyle
b_{n}^{(\ell)}(m-n)\alpha_{n}^{(\ell)}(m)=\sum_{\ell^{\prime}=1}^{\ell}\mathcal{E}_{n}^{(\ell^{\prime})}\alpha_{n}^{(\ell-\ell^{\prime})}(m)+\sum_{\mu=1}^{\infty}\mu\gamma_{0}^{\mu}\alpha_{n}^{(\ell)}(m+\mu)$
$\displaystyle\hphantom{b_{n}^{(\ell)}(m-n)\alpha_{n}^{(\ell)}(m)=}{}+\sum_{\ell^{\prime}=0}^{\ell-1}\sum_{\mu=1}^{\infty}\sum_{k=1}^{\infty}\mu\gamma_{k}^{\mu}\delta_{\ell,\ell^{\prime}+k\mu}\big{[}\alpha_{n}^{(\ell^{\prime})}(m+\mu)+\alpha_{n}^{(\ell^{\prime})}(m-\mu)\big{]},$
(5.11)

where $b_{n}^{(\ell)}(k)$ is short for
$-\kappa\ell+E_{n+k}^{(0)}-\mathcal{E}_{n}^{(0)}$ for $\ell=0,1,2,\ldots$. The
first condition in (5.6) implies $\mathcal{E}_{n}^{(0)}=E^{(0)}_{n}$, and one
obtains the formula for $b_{n}^{(\ell)}(k)$ in (5.2).

One can check that the first two conditions in (5.6) are consistent with
(5.11) (details are given in Section 4.2). Moreover, they imply that it is
consistent to set

$\displaystyle\alpha^{(\ell)}_{n}(m)=0\qquad\forall\,m>n+\ell,\quad\ell\geq 0$
(5.12)

(this can be proved by induction: (5.12) is true by assumption for $\ell=0$;
for $\ell\geq 1$ it follows from (5.11) that $\alpha_{n}^{(\ell)}(m)$ for
$m>n+\ell$ is a linear combination of
$\alpha_{n}^{(\ell^{\prime})}(m^{\prime})$ with $\ell^{\prime}<\ell$,
$m^{\prime}\geq m-\mu$ and $\mu\geq 1$ constrained by the Kronecker deltas in
(5.11), i.e., $m^{\prime}\geq m-(\ell-\ell^{\prime})>n+\ell^{\prime}$, which
proves the claim). Thus one can restrict the second sum in (5.11) to
$1\leq\mu\leq n-m+\ell$. One can also check that the infinite sums in the
terms in the second line in (5.11) can be replaced by finite ones: the
Kronecker deltas there are non-zero for $\ell-\ell^{\prime}=2k\mu$ and
$\ell-\ell^{\prime}=(2k-1)\mu$, and since $\mu,k\geq 1$ this is possible only
if $1\leq\mu\leq\ell-\ell^{\prime}$ and $1\leq
k\leq(\ell-\ell^{\prime})/(2\mu)+\frac{1}{2}$. Since $b_{n}^{(\ell)}(k)$ is
non-zero for $m\neq n$ by assumption, we obtain (5.4) from (5.11). For $m=n$
and $\ell\geq 1$, the first and last conditions in (5.6) make (5.11) into an
equation determining $\mathcal{E}^{(\ell)}_{n}$ as in (5.5). ∎

###  5.2 Algorithm II

For $\kappa\neq 0$, the condition $\alpha_{n}^{(\ell\geq 1)}(n)=0$ in (5.6)
can be replaced by the condition $\mathcal{E}_{n}^{(\ell\geq 1)}=0$ and, in
this way, a somewhat simpler solution algorithm is obtained.

Recall the definitions of $\lambda$ in (2.7)$,\gamma_{\mu}^{k}$ in (3.3) and
$b_{n}^{(\ell)}(k)$ in (5.2).

######  Proposition 5.2.

Let $n\in\mathbb{Z}$, $-\lambda\notin\mathbb{N}_{0}$ for $n>0$ and
$-(g_{0}+g_{1})\notin\mathbb{N}_{0}$ for $n<0$,
$\big{\\{}f^{(\ell)}_{m}(z)\big{\\}}_{m\in\mathbb{Z}}$ the functions defined
and characterized in Lemma 2.2, and assume that $\Im(\kappa)\neq 0$ and
$g_{0}+g_{1}\in\mathbb{R}$.777The latter are no-resonance conditions discussed
in Section 2.4. Then the non-stationary Heun equation in (1.2) has a unique
solution $\psi_{n}(x)$, $E_{n}$ as in (1.5)–(1.6) given by (3.2), (5.3) and

$\displaystyle\mathcal{E}_{n}=\left(n+\frac{g_{0}+g_{1}}{2}\right)^{2},$
(5.13)

with coefficients $\alpha_{n}^{(\ell)}(m)$ determined by the following
recursion relations,

$\displaystyle\alpha_{n}^{(\ell)}(m)=\frac{1}{b_{n}^{(\ell)}(m-n)}\Biggl{(}\sum_{\mu=1}^{n-m+\ell}\mu\gamma_{0}^{\mu}\alpha_{n}^{(\ell)}(m+\mu)$
$\displaystyle\hphantom{\alpha_{n}^{(\ell)}(m)=}{}+\sum_{\ell^{\prime}=0}^{\ell-1}\sum_{\mu=1}^{\ell-\ell^{\prime}}\sum_{k=1}^{\left\lfloor{\frac{\ell-\ell^{\prime}}{\mu}}\right\rfloor}\mu\gamma_{k}^{\mu}\delta_{\ell,\ell^{\prime}+k\mu}\big{[}\alpha_{n}^{(\ell^{\prime})}(m+\mu)+\alpha_{n}^{(\ell^{\prime})}(m-\mu)\big{]}\Biggr{)}$
(5.14)

for all $m\in{\mathbb{Z}^{\prime}}$ if $\ell=0$ and all $m\in\mathbb{Z}$ if
$\ell\geq 1$, together with the following conditions

$\displaystyle\alpha_{n}^{(0)}(n)=1,\qquad\alpha_{n}^{(\ell)}(m)=0\qquad\forall\,m>n+\ell,\quad\ell\geq
0.$ (5.15)

The functions $\mathcal{P}_{n}^{(\ell)}(z)$ thus obtained are polynomials of
degree $n+\ell$ in $z$ for $n+\ell\geq 0$ and zero otherwise.

It follows from (5.1) that the first and last conditions in (5.6) are
equivalent to $\alpha_{n}^{I}(n)=1$, which give rise to
$\mathcal{E}_{n}^{I}=E^{(0)}_{n}+O(q)$, whereas the condition in (5.15) and
$\mathcal{E}_{n}^{(\ell\geq 1)}=0$ corresponding to relaxing the former
conditions to $\alpha_{n}^{II}(n)=1+O(q)$ and instead requiring
$\mathcal{E}_{n}^{II}=E^{(0)}_{n}$ (the superscripts $I$ and $II$ here are to
distinguish the results from the algorithms in Proposition 5.1 and 5.2,
respectively). From this and (3.6) we conclude that the results of the two
algorithms are related as follows,

$\displaystyle\alpha_{n}^{I}(m)=\frac{\alpha_{n}^{II}(m)}{\alpha_{n}^{II}(n)},\qquad\mathcal{E}_{n}^{I}=E^{(0)}_{n}+\frac{1}{\alpha_{n}^{II}(n)}\kappa
q\frac{\partial}{\partial q}\alpha_{n}^{II}(n).$ (5.16)

This shows that the algorithm in Proposition 5.2 is a re-summation of the one
in Proposition 5.1: it follows from (5.16) that, by computing
$\alpha^{II}_{n}(n)$ up to some order $O(q^{N+1})$, one obtains a Padé
approximation of the generalized eigenvalues $\mathcal{E}_{n}^{I}$ as follows,

$\displaystyle\mathcal{E}_{n}^{I}=E^{(0)}_{n}+\frac{\kappa\sum\limits_{\ell=1}^{N}\ell\alpha_{n}^{(II,\ell)}(n)q^{\ell}}{1+\sum\limits_{\ell=1}^{N-1}\alpha_{n}^{(II,\ell)}(n)q^{\ell}}+O\big{(}q^{N+1}\big{)}.$
(5.17)

As an example we give

$\displaystyle\alpha_{n}^{(II,1)}(n)=\frac{\mathcal{E}_{n}^{(1)}}{\kappa},\qquad\alpha_{n}^{(II,2)}(n)=\frac{\mathcal{E}_{n}^{(2)}}{2\kappa}+\frac{\big{(}\mathcal{E}_{n}^{(1)}\big{)}^{2}}{2\kappa^{2}},$

with $\mathcal{E}_{n}^{(1)}$ and $\mathcal{E}_{n}^{(2)}$ in (5.8) and (5.9),
which is obtained by straightforward computations using the algorithm in
Proposition 5.2; note that, when this is inserted in (5.17) for $N=2$, the
singularities at $\kappa=0$ in the numerator and denominator cancel. From our
results above it is clear that this cancellation of singularities at
$\kappa=0$ happens for arbitrary $N$ (this can be proved by the same argument
we used to prove Lemma 2.4 in Section 2.4). This implies that the limit
$\kappa\to 0$ is well-defined (in the sense made precise in Lemma 2.4), and we
thus obtain from (5.17) an interesting representation of the eigenvalues of
the $BC_{1}$ Inozemtsev model.

######  Remark 5.3.

It is interesting to note that (5.17) provides a simpler representation of the
Taylor coefficients $\mathcal{E}_{n}^{(\ell)}$ of the eigenvalues of the Heun
equation for larger values of $\ell$: we computed $\mathcal{E}_{n}^{(I,\ell)}$
and $\alpha_{n}^{(II,\ell)}(n)$ using MAPLE for $\ell=1,2,3,4,5$ and found
that the former have $2$, $18$, $162$, $1776$, $21002$ distinct terms,888The
numbers $2$ and $18$ in comparison with (5.8) and (5.9) explain what we mean
by “number of distinct terms”. whereas the latter have $2$, $18$, $148$,
$1298$, $11632$ distinct terms, respectively. Moreover, as discussed in
Appendix E, there exist explicit formulas for the coefficients
$\mathcal{E}_{n}^{(I,\ell)}$ and $\alpha_{n}^{(II,\ell)}(n)$, and the ones for
the latter are much simpler than the ones for the former.

##  6 Solution to all orders

In this section we use results in [24] to obtain perturbative solutions of the
equation in (1.2) to all orders, in generalization of the results in Section
4.2 for $q=0$.

As shown in Section 4, for all $n\in\mathbb{Z}$, (1.2) has solutions
$\psi_{n}(x)$, $E_{n}$ as in (1.5) and (3.1)–(3.2) provided that the
coefficients $\alpha_{n}(m)$ and the (shifted) eigenvalue $\mathcal{E}_{n}$
are (suitable) solutions of the differential-difference equations in (3.4). We
observe that, by introducing the notation999As explained further below, the
following defines operators $\mathbb{A}$ and $\mathbb{B}$ on the vector space
of sequences $\alpha=\\{\alpha(m)\\}_{m\in\mathbb{Z}}$.

$\displaystyle(\mathbb{A}\alpha)(m)\equiv\left(-\kappa
q\frac{\partial}{\partial
q}+E^{(0)}_{m}\right)\alpha(m),\qquad(\mathbb{B}\alpha)(m)\equiv\sum_{\mu\in{\mathbb{Z}^{\prime}}}S_{\mu}\alpha(m+\mu)$
(6.1)

with $E^{(0)}_{m}$ in (2.4) and $S_{\mu}$ in (4.11), one can write (3.4) in a
simple way as follows,

$\displaystyle[\mathbb{A}-\mathcal{E}]\alpha(m)=(\mathbb{B}\alpha)(m)$ (6.2)

(we find it convenient to suppress the dependence on $n$ here). We note that
the notation introduced here is such that we can use results in [24, Section
4.1] as they stand. We first consider the special case $\kappa=0$ (Heun
equation) in Section 6.1, and then the general case in Section 6.2. We note
that Theorem 6.1 is a generalization of Proposition 5.1 restricted to
$\kappa=0$; Theorem 6.3 generalizes Proposition 5.2.

###  6.1 Heun equation

The Heun equation corresponds to the special case $\kappa=0$ of (1.2). This is
an eigenvalue equation for a 1D Schrödinger operator, and $\psi(x)$ and $E$
have a quantum mechanical interpretation as energy eigenfunction and
eigenvalue, respectively. In this section we show how to use results in [24]
to obtain a perturbative solution for this case to all orders.

Using the Feshbach method and expanding a Neumann series yields implicit
solutions of (6.1)–(6.2) for $\kappa=0$ to all orders [24]. To formulate this
result we introduce the following two functions of a complex variable $z$
($n\in\mathbb{N}_{0}$ and $m\in\mathbb{Z}$ are fixed),

$\displaystyle\Phi_{n}(z)\equiv-\sum_{s=2}^{\infty}\sum_{\mu_{1}\in{\mathbb{Z}^{\prime}}}S_{\mu_{1}}\cdots\sum_{\mu_{s}\in{\mathbb{Z}^{\prime}}}S_{\mu_{s}}$
$\displaystyle\hphantom{\Phi_{n}(z)\equiv}{}\times\frac{\delta(0,\mu_{1}+\cdots+\mu_{s})}{\big{[}b^{(0)}_{n}(\mu_{1})-z\big{]}_{n}\big{[}b^{(0)}_{n}(\mu_{1}+\mu_{2})-z\big{]}_{n}\cdots\big{[}b^{(0)}_{n}(\mu_{1}+\cdots+\mu_{s-1})-z\big{]}_{n}}$
$\displaystyle\hphantom{\Phi_{n}(z)}{}=-\sum
S_{\mu_{1}}S_{\mu_{2}}\frac{\delta(0,\mu_{1}+\mu_{2})}{\big{[}b^{(0)}_{n}(\mu_{1})-z\big{]}_{n}}$
$\displaystyle\hphantom{\Phi_{n}(z)\equiv}{}-\sum
S_{\mu_{1}}S_{\mu_{2}}S_{\mu_{3}}\frac{\delta(0,\mu_{1}+\mu_{2}+\mu_{3})}{\big{[}b^{(0)}_{n}(\mu_{1})-z\big{]}_{n}\big{[}b^{(0)}_{n}(\mu_{1}+\mu_{2})-z\big{]}_{n}}+\cdots$
(6.3)

and

$\displaystyle
G_{n}(z;m)\equiv\delta(m,n)+\sum_{s=1}^{\infty}\sum_{\mu_{1}\in{\mathbb{Z}^{\prime}}}S_{\mu_{1}}\cdots\sum_{\mu_{s}\in{\mathbb{Z}^{\prime}}}S_{\mu_{s}}$
$\displaystyle\hphantom{G_{n}(z;m)\equiv}{}\times\frac{\delta(m,n+\mu_{1}+\cdots+\mu_{s})}{\big{[}b^{(0)}_{n}(m\\!-\\!n)\\!-\\!z\big{]}_{n}\big{[}b^{(0)}_{n}(m\\!-\\!n\\!+\\!\mu_{1})\\!-\\!z\big{]}_{n}\cdots\big{[}b^{(0)}_{n}(m\\!-\\!n\\!+\\!\mu_{1}\\!+\\!\cdots\\!+\\!\mu_{s-1})\\!-\\!z\big{]}_{n}}$
$\displaystyle\hphantom{G_{n}(z;m)}{}=\delta(n,m)+\sum
S_{\mu}\frac{\delta(n,m+\mu)}{\big{[}b^{(0)}_{n}(m-n)-z\big{]}_{n}}$
$\displaystyle\hphantom{G_{n}(z;m)\equiv}{}+\sum
S_{\mu_{1}}S_{\mu_{2}}\frac{\delta(n,m+\mu_{1}+\mu_{2})}{\big{[}b^{(0)}_{n}(m-n)-z\big{]}_{n}\big{[}b^{(0)}_{n}(m-n+\mu_{1})-z\big{]}_{n}}+\cdots$
(6.4)

with the Kronecker delta $\delta(n,m)$, $S_{\mu}$ in (4.11), and the
convenient shorthand notations

$\displaystyle\frac{1}{[b^{(0)}_{n}(k)-z]_{n}}\equiv\begin{cases}1/\big{[}E^{(0)}_{n+k}-E^{(0)}_{n}-z\big{]},&k\neq
0,\\\ 0,&k=0.\end{cases}$ (6.5)

It is important to note that these formulas should be interpreted in a
perturbative sense as follows: for fixed $N\in\mathbb{N}$, there are only
finitely many terms in (6.3) and (6.4) that are $O(q^{\ell})$ with $\ell<N$
and thus, to obtain results up to $O(q^{N})$-corrections, all infinite series
in our formulas can be truncated to finite series. The reason why it is
convenient to write infinite series is that it is difficult to give a simple
general recipe for finite-$N$ truncations; for example, in (6.3) one can
restrict to $2\leq s\leq N$ and $-N\leq\mu_{j}\leq N$ for $j=1,2,\ldots,s$,
but this is somewhat arbitrary since the resulting finite sum still contains
many terms which do not contribute to $O(q^{N})$.

######  Theorem 6.1.

Let $n\in\mathbb{Z}$, $-\lambda\notin\mathbb{N}_{0}$ for $n>0$ and
$-(g_{0}+g_{1})\notin\mathbb{N}_{0}$ for $n<0$,
$\big{\\{}f^{(\ell)}_{m}(z)\big{\\}}_{m\in\mathbb{Z}}$ the functions defined
and characterized in Lemma 2.2. Then the Heun equation in (1.2) for $\kappa=0$
and $g_{0}+g_{1}\notin\mathbb{Z}$ has a unique solution as in (1.5) and
(3.1)–(3.2) provided that

$\displaystyle\mathcal{E}_{n}=E^{(0)}_{n}+\tilde{\mathcal{E}}_{n},\qquad\alpha_{n}(m)=G_{n}\big{(}\tilde{\mathcal{E}}_{n};m\big{)}$

with $\tilde{\mathcal{E}}_{n}$ the unique solution of the equation

$\displaystyle\tilde{\mathcal{E}}_{n}=\Phi_{n}\big{(}\tilde{\mathcal{E}}_{n}\big{)}$

vanishing like $O(q)$ as $q\to 0$.

(The proof can be found at the end of this section.)

From this result one can obtain the following fully explicit formula for the
generalized eigenvalues of the non-stationary Heun equation,

$\displaystyle\mathcal{E}_{n}=E^{(0)}_{n}+\sum_{m=1}^{\infty}\sum_{k_{0},k_{1},\ldots,k_{m-1}\in\mathbb{N}}\\!\delta\left(\sum_{r=0}^{m-1}k_{r},m\right)\\!\delta\left(\sum_{r=1}^{m-1}rk_{r},m-1\right)\\!(m-1)!\prod_{r=1}^{m}\frac{[\Phi_{n}^{(r)}]^{k_{r}}}{k_{r}!}\\!\\!\\!$
(6.6)

with

$\displaystyle\Phi^{(r)}_{n}(z)\equiv-\sum_{s=2}^{\infty}\sum_{\mu_{1},\ldots,\mu_{s}\in{\mathbb{Z}^{\prime}}}\sum_{r_{1},\ldots,r_{s-1}\in\mathbb{N}_{0}}S_{\mu_{1}}\cdots
S_{\mu_{s}}$
$\displaystyle\hphantom{\Phi^{(r)}_{n}(z)\equiv}{}\times\frac{\delta(0,\mu_{1}+\cdots+\mu_{s})\delta(r_{1}+\cdots
r_{s-1},r)}{[b^{(0)}_{n}(\mu_{1})]^{1+r_{1}}_{n}[b^{(0)}_{n}(\mu_{1}+\mu_{2})]^{1+r_{2}}_{n}\cdots[b^{(0)}_{n}(\mu_{1}+\cdots+\mu_{s-1})]^{1+r_{s-1}}_{n}},$
(6.7)

and similarly for $\alpha_{n}(m)$; see Theorem 4.3.1 in [24] (one can check
that the proof of the later theorem applies as it stands to the present case).
As explained in Appendix E, this formula can be used to turn the computation
of the generalized eigenvalues $\mathcal{E}_{n}$ of the non-stationary Heun
equation into a combinatorial problem; see (E.1) ff.

######  Remark 6.2.

The elliptic generalizations of the Jacobi polynomials
$P^{\big{(}g_{0}-\frac{1}{2},g_{1}-\frac{1}{2}\big{)}}_{n}(z)$ provided by
Theorem 6.1 can be formally written as

$\displaystyle\mathcal{P}_{n}(z)={\mathcal{N}}_{n}\oint_{\mathcal{C}}\frac{d\xi}{2\pi{\rm
i}\xi}\frac{\prod\limits_{\nu=0}^{3}\Theta_{\nu+1}(\xi)^{\tilde{g}_{\nu}}}{\Theta(z,\xi)^{\lambda}}\tilde{\mathcal{P}}_{n}(\xi)$

with

$\displaystyle\tilde{\mathcal{P}}_{n}(\xi)=\xi^{-n}+\sum_{s=1}^{\infty}\prod_{\kappa=1}^{s}\left(\sum_{\mu_{\kappa}\in{\mathbb{Z}^{\prime}}}\frac{S_{\mu_{\kappa}}}{\Big{[}b_{n}^{(0)}(m-n+\sum\limits_{\ell=1}^{\kappa-1}\mu_{\ell})-\tilde{\mathcal{E}}_{n}\Big{]}_{n}}\right)\xi^{-n-\mu_{1}-\cdots-\mu_{s}}.$

It is possible to identify $\tilde{\mathcal{P}}_{n}(\xi)$ with a singular
eigenfunction of the Inozemtsev Hamiltonian
$H(y;\\{\tilde{g}_{\nu}\\}_{\nu=0}^{3})$ appearing in the generalized kernel
function identity in Lemma 4.1 and, in this way, extend the interpretation of
the kernel function method given in [22] to the present case.

###### Proof of Theorem 6.1.

(We are brief since interested readers can find further details in [24]. In
particular, it is explained in [24] why we can ignore questions of convergence
in the argument below.)

Let $V$ be the vector space of sequences $\\{\alpha(m)\\}_{m\in\mathbb{Z}}$
and regard $\mathbb{A}$ and $\mathbb{B}$ in (6.1) as linear operators $V\to
V$. Define projections $\mathbb{P}$ on $V$ as follows,101010We write
$\mathbb{P}$ short for $\mathbb{P}_{n}$.

$\displaystyle(\mathbb{P}\alpha)(m)\equiv\delta(n,m)\alpha(m)\qquad\forall\,\alpha\in
V$ (6.8)

and $\mathbb{P}^{\perp}\equiv I-\mathbb{P}$. Then $\mathbb{A}$ in (6.1)
commutes with $\mathbb{P}$, and the equation $\mathbb{P}\alpha_{0}=\alpha_{0}$
is solved by $(\alpha_{0})(m)\equiv\delta(n,m)$. Thus Lemma 4.1.1 in [24]
implies that

$\displaystyle\alpha=\big{[}I+\big{(}\mathbb{A}-\mathcal{E}-\mathbb{P}^{\perp}\mathbb{B}\big{)}^{-1}\mathbb{P}^{\perp}\mathbb{B}\big{]}\alpha_{0},\qquad\mathcal{E}\alpha_{0}=\mathbb{A}\alpha_{0}-\mathbb{P}\mathbb{B}\alpha$
(6.9)

is a solution of (6.2). Expanding the resolvent in a Neumann series we obtain
[24]

$\displaystyle\alpha(m)=\sum_{s=0}^{\infty}\big{(}\big{(}[\mathbb{A}-\mathcal{E}]^{-1}\mathbb{P}^{\perp}\mathbb{B}\big{)}^{s}\alpha_{0}\big{)}(m),\\!\\!\\!\qquad\mathcal{E}=E^{(0)}_{n}-\sum_{s=0}^{\infty}\big{(}\mathbb{B}\big{(}[\mathbb{A}-\mathcal{E}]^{-1}\mathbb{P}^{\perp}\mathbb{B}\big{)}^{s}\alpha_{0}\big{)}(n).\\!\\!\\!\\!\\!$
(6.10)

The ansatz $\mathcal{E}=E^{(0)}_{n}+\tilde{\mathcal{E}}$ implies

$\displaystyle\big{(}[\mathbb{A}-\mathcal{E}]^{-1}\mathbb{P}^{\perp}\alpha\big{)}(m)=\frac{1}{\big{[}b^{(0)}_{n}(m-n)-\tilde{\mathcal{E}}\big{]}_{n}}\alpha(m)$
(6.11)

using the shorthand notation in (6.5). Using (6.11) to compute (6.10) we
obtain $\alpha(m)=G_{n}(\tilde{\mathcal{E}};m)$ and
$\tilde{\mathcal{E}}=\Phi_{n}(\tilde{\mathcal{E}})$ with the functions defined
in (6.3)–(6.4) [24]. The latter should be interpreted as non-linear equation
whose solution $\tilde{\mathcal{E}}=\tilde{\mathcal{E}}_{n}$ vanishing like
$O(q)$ determines the solution of the Heun equation we are interested in; see
[24] for details. ∎

###  6.2 Time dependent Heun equation

We now present generalizations of the results in the previous section to the
non-stationary case $\kappa\neq 0$.

The argument to solve (6.1)–(6.2) in the proof of Theorem 6.1 can be
generalized to $\kappa\neq 0$ if the following projection is used,

$\displaystyle(\mathbb{P}\alpha)^{(\ell)}(m)=\delta(\ell,0)\delta(n,m)\alpha^{(\ell)}(m)$
(6.12)

where $\alpha^{(\ell)}(m)$ are defined as coefficients of the formal power
series in $q$ (see (5.1)); with that, the results in (6.9) hold true as they
stand. However, in the present case, the second of these equations is
trivially solved by $\mathcal{E}=E^{(0)}_{n}$, and this implies a stronger
result:

######  Theorem 6.3.

Let $n\in\mathbb{Z}$, $-\lambda\notin\mathbb{N}_{0}$ for $n>0$ and
$-(g_{0}+g_{1})\notin\mathbb{N}_{0}$ for $n<0$,
$\big{\\{}f^{(\ell)}_{m}(z)\big{\\}}_{m\in\mathbb{Z}}$ the functions defined
and characterized in Lemma 2.2, and assume that $\Im(\kappa)\neq 0$ and
$g_{0}+g_{1}\in\mathbb{R}$.111111The latter are no-resonance conditions
discussed in Section 2.4. Then the non-stationary Heun equation in (1.2) has a
unique solution $\psi_{n}(x)$, $E_{n}$ as in (1.5)–(1.6) given by (3.2),
(5.3), (5.13) and the coefficients

$\displaystyle\alpha^{(\ell)}_{n}(m)=\sum_{s=0}^{\infty}\sum_{k_{1},\ldots,k_{s}\in\mathbb{N}_{0}}\sum_{\mu_{1},\ldots,\mu_{s}\in{\mathbb{Z}^{\prime}}}S_{\mu_{1}}(k_{1})\cdots
S_{\mu_{s}}(k_{s})$
$\displaystyle\times\frac{\delta(\ell,|\mu_{1}|k_{1}+\cdots+|\mu_{s}|k_{s})\delta(n,m+\mu_{1}+\cdots+\mu_{s})}{\big{[}b^{(\ell)}_{n}(m\\!-\\!n)\big{]}\big{[}b^{(\ell-|\mu_{1}|k_{1})}_{n}(m\\!-\\!n\\!+\\!\mu_{1})\big{]}\cdots\big{[}b^{(\ell-|\mu_{1}|k_{1}\\!-\\!\cdots\\!-\\!|\mu_{s-1}|k_{s-1})}_{n}(m\\!-\\!n\\!+\\!\mu_{1}\\!+\\!\cdots\\!+\\!\mu_{s-1})\big{]}}$
$\displaystyle=\delta(\ell,0)\delta(n,m)+\sum
S_{\mu}(k)\frac{\delta(\ell,|\mu|k)\delta(n,m+\mu)}{[b_{n}^{(\ell)}(n-m)]}$
$\displaystyle+\sum
S_{\mu_{1}}(k_{1})S_{\mu_{2}}(k_{2})\frac{\delta(\ell,|\mu_{1}|k_{1}+|\mu_{2}|k_{2})\delta(n,m+\mu_{1}+\mu_{2})}{\big{[}b_{n}^{(\ell)}(n-m)\big{]}\big{[}b_{n}^{(\ell-|\mu_{1}|k_{1})}(n-m+\mu_{1})\big{]}}+\cdots$
(6.13)

using the shorthand notations

$\displaystyle S_{\mu}(k)\equiv\begin{cases}\mu\gamma_{0}^{\mu}&\mbox{if $k=0$
and $\mu>0$},\\\ |\mu|\gamma_{k}^{\mu}&\mbox{if $k\in\mathbb{N}$},\\\
0&\mbox{otherwise }\end{cases}$ (6.14)

with $\gamma_{k}^{\mu}$ in (3.3), and

$\displaystyle\frac{1}{[b^{(\ell)}_{n}(k)]}\equiv\begin{cases}0&\mbox{if
}\ell=0\mbox{ and }k=0,\\\ 1/b_{n}^{(\ell)}(k)&\mbox{otherwise}\end{cases}$
(6.15)

with $b_{n}^{(\ell)}(k)$ in (5.2).

(The proof can be found at the end of this section.)

It is not difficult to convince one-selves that the coefficients in (6.13) are
identical with the ones determined by the algorithm in Proposition 5.2. Thus,
by setting $m=n$ and using the second identity in (5.16), we obtain a formula
to all orders for the generalized eigenvalues $\mathcal{E}_{n}$ of the time
dependent Heun equation:

$\displaystyle\mathcal{E}_{n}=E^{(0)}_{n}+\frac{\sum\limits_{\ell=1}^{\infty}\kappa\ell\alpha_{n}^{(\ell)}(n)q^{\ell}}{1+\sum\limits_{\ell=1}^{\infty}\ell\alpha_{n}^{(\ell)}(n)q^{\ell}}$
(6.16)

with $\alpha_{n}(m)$ in (6.13) for $m=n$; the limit $\kappa\to 0$ of this
formula is non-trivial but well-defined. It would be interesting to find a re-
summation which makes this manifest.

We thus obtained two explicit formulas for the eigenvalues of the Heun
equation in (1.2) for $\kappa=0$: the formula in (6.6), and the limit
$\kappa\to 0$ of the formula in (6.16). The former has the advantage that it
is manifestly finite for $\kappa=0$, whereas the latter requires a non-trivial
limit. However, as explained in Appendix E, the latter formula is much simpler
than the former.

###### Proof of Theorem 6.3.

Let $V$ be the vector space of all sequences
$\alpha=\\{\alpha^{(\ell)}(m)\\}_{\ell\in\mathbb{N}_{0},m\in\mathbb{Z}}$
identified with $\alpha=\\{\alpha(m)\\}_{m\in\mathbb{Z}}$ as in (5.1). Then
$\mathbb{A}$ and $\mathbb{B}$ in (6.1) can be written as linear operators
$V\to V$ as follows,

$\displaystyle(\mathbb{A}\alpha)^{(\ell)}(m)=\big{(}{-}\kappa\ell+E^{(0)}_{m}\big{)}\alpha^{(\ell)}(m),\qquad(\mathbb{B}\alpha)^{(\ell)}(m)=\sum_{\ell^{\prime}=0}^{\ell}\sum_{\mu\in{\mathbb{Z}^{\prime}}}S^{(\ell^{\prime})}_{\mu}\alpha^{(\ell-\ell^{\prime})}(m+\mu)$
(6.17)

with

$\displaystyle S_{\mu}^{(\ell)}\equiv\begin{cases}\mu\gamma_{0}^{\mu}&\mbox{if
$\ell=0$, $\mu>0$},\\\ |\mu|\gamma_{\ell/|\mu|}^{\mu}&\mbox{if
$\ell/|\mu|\in\mathbb{N}$},\\\ 0&\mbox{otherwise }\end{cases}$ (6.18)

for $\mu\in{\mathbb{Z}^{\prime}}$ and $\ell\in\mathbb{N}_{0}$ (the latter
follows with (3.3), (4.9), and (4.11)). This allows to rewrite (3.4) as in
(6.2).

It is clear that the projection $\mathbb{P}$ in (6.12) commutes with
$\mathbb{A}$ in (6.17), and that the equation
$\mathbb{P}\alpha_{0}=\alpha_{0}$ is solved by

$\displaystyle(\alpha_{0})^{(\ell)}(m)=\delta(\ell,0)\delta(n,m).$ (6.19)

Thus Lemma 4.1.1 in [24] implies that the solution $\alpha$ and $\mathcal{E}$
of (3.4) satisfies (6.9) with that $\alpha_{0}$
($\mathbb{P}^{\perp}=I-\mathbb{P}$).

The definitions imply

$\displaystyle(\mathbb{P}\mathbb{B}\alpha)^{(\ell)}(m)=\delta(\ell,0)\delta(n,m)\sum_{\mu=1}^{\infty}\mu\gamma_{0}^{\mu}\alpha^{(0)}(n+\mu).$

We showed already in Section 4.2 that the solution $\alpha$ of (3.4) is such
that $\alpha^{(0)}(m>n)=0$, and thus $\mathbb{P}\mathbb{B}\alpha=0$. With that
we find that the second equation in (6.9) is solved by
$\mathcal{E}=\mathcal{E}^{(0)}_{n}$.121212This result was already obtained by
a different argument in the proof of Proposition 5.2. Using this the first
equation in (6.10) simplifies to

$\displaystyle\alpha(m)=\sum_{s=0}^{\infty}\big{(}\big{(}\big{[}\mathbb{A}-E^{(0)}_{n}\big{]}^{-1}\mathbb{P}^{\perp}\mathbb{B}\big{)}^{s}\alpha_{0}\big{)}(m)$
(6.20)

Inserting definitions we find

$\displaystyle\big{(}\big{[}\mathbb{A}-E^{(0)}_{n}\big{]}^{-1}\mathbb{P}^{\perp}\mathbb{B}\alpha\big{)}^{(\ell)}(m)=\frac{1}{\big{[}b^{(\ell)}_{n}(m-n)\big{]}}\sum_{\ell^{\prime}=0}^{\ell}\sum_{\mu\in{\mathbb{Z}^{\prime}}}S^{(\ell^{\prime})}_{\mu}\alpha^{(\ell-\ell^{\prime})}(m+\mu)$

using the notation in (6.15). With that one can prove by induction that

$\displaystyle\big{(}\big{(}\big{[}\mathbb{A}-E^{(0)}_{n}\big{]}^{-1}\mathbb{P}^{\perp}\mathbb{B}\big{)}^{s}\alpha\big{)}^{(\ell)}(m)=\sum_{\ell_{1}=0}^{\ell}\sum_{\ell_{2}=0}^{\ell-\ell_{1}}\cdots\sum_{\ell_{s}=0}^{\ell-\ell_{1}-\cdots-\ell_{s-1}}\sum_{\mu_{1},\ldots,\mu_{s}\in{\mathbb{Z}^{\prime}}}S^{(\ell_{1})}_{\mu_{1}}\cdots
S^{(\ell_{s})}_{\mu_{s}}$
$\displaystyle\qquad{}\times\frac{1}{\big{[}b^{(\ell)}_{n}(m-n)\big{]}\big{[}b^{(\ell-\ell_{1})}_{n}(m-n+\mu_{1})\big{]}\cdots\big{[}b^{(\ell-\ell_{1}-\cdots-\ell_{s-1})}_{n}(m-n+\mu_{1}+\cdots+\mu_{s-1})\big{]}}$
$\displaystyle\qquad{}\times\alpha^{(\ell-\ell_{1}-\cdots-\ell_{s})}(m+\mu_{1}+\cdots+\mu_{s})$

for all $s=1,2,\ldots$. Inserting this in (6.20) and using (6.19) give

$\displaystyle\alpha^{(\ell)}_{n}(m)=\sum_{s=0}^{\infty}\sum_{\ell_{1},\ldots,\ell_{s}\in\mathbb{N}_{0}}\sum_{\mu_{1},\ldots,\mu_{s}\in{\mathbb{Z}^{\prime}}}S^{(\ell_{1})}_{\mu_{1}}\cdots
S^{(\ell_{s})}_{\mu_{s}}$ (6.21)
$\displaystyle\hphantom{\alpha^{(\ell)}_{n}(m)=}{}\times\frac{\delta(\ell,\ell_{1}+\cdots+\ell_{s})\delta(n,m+\mu_{1}+\cdots+\mu_{s})}{b^{(\ell)}_{n}(m-n)b^{(\ell-\ell_{1})}_{n}(m-n+\mu_{1})\cdots
b^{(\ell-\ell_{1}-\cdots-\ell_{s-1})}_{n}(m-n+\mu_{1}+\cdots+\mu_{s-1})}$

(to simplify notation we extended the $\ell_{j}$-summations to all non-
negative integers, which is possible due to the first Kronecker delta in the
sum; the term for $s=0$ is by definition equal to the r.h.s. in (6.19)).

The definition of $S_{\mu}^{(\ell)}$ in (6.18) suggests to change summation
variables from $\ell$ to $k$ such that $\ell=|\mu|k$ (to reduce the number of
terms in the formula which give zero contributions). Writing $S_{\mu}(k)$
short for $S_{\mu}^{(|\mu|k)}$ we obtain the result in (6.13)–(6.14). ∎

##  7 Final remarks

We showed that a solution method based on kernel functions and developed to
solve the elliptic Calogero–Sutherland (eCS) model [19, 24] generalizes to the
non-stationary Heun equation in (1.2). This suggests to us that also other
elliptic problems can be treated by this method; as a possible candidate we
mention the non-stationary generalization of the eigenvalue problem for the
$BC_{N}$ Inozemtsev model for $N\geq 2$ [12], which defines a natural many-
variable generalization of the non-stationary Heun equation:

$\displaystyle\Biggl{(}\frac{{\rm
i}}{\pi}\kappa\frac{\partial}{\partial\tau}+\sum_{j=1}^{N}\Biggl{(}-\frac{\partial^{2}}{\partial
x_{j}^{2}}+\sum_{\nu=0}^{3}g_{\nu}(g_{\nu}-1)\wp(x_{j}+\omega_{\nu})\Biggr{)}$
$\displaystyle\qquad{}+2\lambda(\lambda-1)\sum_{1\leq j<k\leq
N}\left\\{\wp(x_{j}-x_{k})+\wp(x_{j}+x_{k})\right\\}\Biggr{)}\psi(x)=E\psi(x);$
(7.1)

note that generalized kernel function identities for this problem and a
discrete set of $\kappa$-values (depending on the other model parameters) are
known [25].

We also obtained results beyond previous results about the eCS model: the non-
stationary Heun equation in (1.2) is invariant under the transformations in
(1.4), and we found that, for $\kappa\neq 0$, this symmetry can be exploited
to construct simpler representations of the generalized eigenvalues $E$ which
are useful even in the case $\kappa=0$; see (6.16). We expect that a formula
similar to the one in (6.16) for the eigenvalues of the eCS model can be
found, and that this would be interesting since it might allow to better
understand the relations between the solutions of the eCS model in [24] and
the one by Nekrasov and Shatiashvili [29].

As discussed in Appendix E, the explicit formulas for the solutions of the
Heun equation in Theorem 6.1 can be regarded as translations of the problem to
solve this equation into a combinatorial problem. It is remarkable that
different elliptic problems lead to the same combinatorial problem: we find
that the combinatorial structure of the solution of the eCS model and the non-
stationary Heun equation are the same, and model details only affect the basic
building blocks of the solutions. We expect the same is true for other
elliptic problems like the one in (7.1). We believe that a similar remark
applies to non-stationary extensions of elliptic problems. It would be
interesting to explore these combinatorial aspects of our solution in future
work.

One important question about the non-stationary Heun equation is uniqueness:
which conditions determine a unique solution? Our results shed light on this
question: we find that the conditions in (1.5)–(1.7) do not fix the solution
uniquely, and our results suggest the following further requirements:

$\displaystyle\mathcal{E}_{n}^{(\ell)}=0,\qquad\mathcal{P}_{n}^{(\ell)}(z)=O\big{(}z^{n+\ell}\big{)}\qquad\forall\,n+\ell\geq
0,\;\ell\geq 1.$ (7.2)

It would be interesting to prove that the conditions (1.5)–(1.7) and (7.2)
imply uniqueness.

We finally note that kernel functions have been used since a long time to
transform the Heun equation into integral equations [3, 18]; see also [13,
31]. This has provided powerful tools to study analytical properties of
solutions; for example, this was used by Ruijsenaars in his work on the hidden
permutation symmetry mentioned after (5.10) [35]. Our approach is similar in
that the kernel function we use determines the analytic structure of the
solutions we construct. However, there are also important differences. For
example, our kernel functions are not given by hypergeometric functions as the
ones in [3, 18], and they are singular and not $L^{2}$ as the ones used in
[35]. Moreover, our emphasis is on constructing explicit representations of
solutions.

##  Appendix A Special functions

For the convenience of the reader we collect definitions and properties of
standard special functions that we use.

###  A.1 Elliptic functions

The Weierstrass $\wp$-function with periods $(2\omega_{1},2\omega_{3})$ is
defined as follows:

$\displaystyle\wp(x\,|\,\omega_{1},\omega_{3})\equiv\frac{1}{x^{2}}+\sum_{(m,n)\in{\mathbb{Z}}\times{\mathbb{Z}}\setminus\\{(0,0)\\}}\left\\{\frac{1}{(x-\Omega_{m,n})^{2}}-\frac{1}{\Omega_{m,n}^{2}}\right\\}$
(A.1)

with $\Omega_{m,n}\equiv 2m\omega_{1}+2n\omega_{3}$. We also recall the
definition of the theta functions,

$\displaystyle\theta_{1}(x)\equiv
2\sum_{n=0}^{\infty}(-1)^{n}q^{(n+1/2)^{2}}\sin(2n+1)x,$
$\displaystyle\theta_{2}(x)\equiv
2\sum_{n=0}^{\infty}q^{(n+1/2)^{2}}\cos(2n+1)x,$
$\displaystyle\theta_{3}(x)\equiv 1+2\sum_{n=1}^{\infty}q^{n^{2}}\cos 2nx,$
$\displaystyle\theta_{4}(x)\equiv 1+2\sum_{n=1}^{\infty}(-1)^{n}q^{n^{2}}\cos
2nx$ (A.2)

with

$\displaystyle q={\rm e}^{{\rm
i}\pi\tau},\qquad\tau=\frac{\omega_{3}}{\omega_{1}}.$ (A.3)

We need the identity

$\displaystyle\wp(z|\pi,\pi\tau)=-\frac{\eta_{1}}{\pi}+\sum_{n\in\mathbb{Z}}\frac{1}{4\sin^{2}\big{(}\frac{1}{2}[z+2n\pi\tau]\big{)}}$
(A.4)

with

$\displaystyle\frac{\eta_{1}}{\pi}=\frac{1}{12}-\sum_{n=1}^{\infty}\frac{2q^{2n}}{(1-q^{2n})^{2}}$
(A.5)

(see [30, 23.8.3 and 23.8.5]).

Let

$\displaystyle G\equiv\prod_{n=1}^{\infty}\big{(}1-q^{2n}\big{)}.$ (A.6)

From the product representations of the theta functions (see [30,
20.5.1–20.5.4]) we obtain

$\displaystyle\theta_{1}\big{(}\tfrac{1}{2}y\big{)}=q^{1/4}{\rm e}^{{\rm
i}\pi/2}{\rm e}^{-{\rm i}y/2}G\Theta_{1}\big{(}{\rm e}^{{\rm i}y}\big{)},$
$\displaystyle\theta_{2}\big{(}\tfrac{1}{2}y\big{)}=q^{1/4}{\rm e}^{-{\rm
i}y/2}G\Theta_{2}\big{(}{\rm e}^{{\rm i}y}\big{)},$
$\displaystyle\theta_{3}\big{(}\tfrac{1}{2}y\big{)}=G\Theta_{3}\big{(}{\rm
e}^{{\rm i}y}\big{)},$
$\displaystyle\theta_{4}\big{(}\tfrac{1}{2}y\big{)}=G\Theta_{4}\big{(}{\rm
e}^{{\rm i}y}\big{)}$ (A.7)

and

$\displaystyle\theta_{1}\big{(}\tfrac{1}{2}(x+y)\big{)}\theta_{1}\big{(}\tfrac{1}{2}(x-y)\big{)}=q^{1/2}G^{2}{\rm
e}^{-{\rm i}y}\Theta\big{(}\cos(x),{\rm e}^{{\rm i}y}\big{)}$ (A.8)

with the functions $\Theta_{\nu}(\xi)$, $\nu=1,2,3,4$, and $\Theta(z,\xi)$
defined in (2.1) and (2.2).

###  A.2 Jacobi polynomials

We use the following definitions of the Pochhammer symbol

$\displaystyle(x)_{n}\equiv x(x+1)\cdots(x+n-1)$ (A.9)

and binomial coefficients

$\displaystyle\binom{x}{n}\equiv\frac{(x-n+1)_{n}}{n!}$ (A.10)

for complex $x$ and $n\in\mathbb{N}_{0}$. We will also use the extensions of
these definitions to all $n\in\mathbb{Z}$ provided by the $\Gamma$ function,
i.e., $(x)_{n}=\Gamma(x+n)/\Gamma(x)$ and similarly for the binomial
coefficients.

The Jacobi polynomials are given by

$\displaystyle
P_{n}^{(\alpha,\beta)}(z)\equiv\sum_{\ell=0}^{n}\frac{(n+\alpha+\beta+1)_{\ell}(\alpha+\ell+1)_{n-\ell}}{\ell!(n-\ell)!}\left(\frac{z-1}{2}\right)^{\ell}$
(A.11)

for $\alpha,\beta\geq-1$, $n=0,1,2,\ldots$, and complex $z$ (see [30,
18.5.7]). Thus

$\displaystyle
P_{n}^{(\alpha,\beta)}(z)=\frac{(n+\alpha+\beta+1)_{n}}{2^{n}n!}z^{n}+O\big{(}z^{n-1}\big{)}.$
(A.12)

We recall the definition of the Gegenbauer (or ultraspherical) polynomials by
their well-known generating function:

$\displaystyle\big{(}1-2z\xi+\xi^{2}\big{)}^{-\lambda}=\sum_{n=0}^{\infty}C^{(\lambda)}_{n}(z)\xi^{n}=\sum_{n=0}^{\infty}\frac{(2\lambda)_{n}}{\left(\lambda+\frac{1}{2}\right)_{n}}P_{n}^{(\lambda-\frac{1}{2},\lambda-\frac{1}{2})}(z)\xi^{n}$
(A.13)

(see [30, 18.12.4]). This is equivalent to the first identity in (3.10) (to
see this, note that $(n+2g)_{n}/[4^{n}(g)_{n}]=(g+\frac{1}{2})_{n}/(2g)_{n}$).

##  Appendix B Scaling

In the literature different conventions for $\omega_{1}$ are used. We
therefore state how our results for $\omega_{1}=\pi$ can be transformed into
results for other values of $\omega_{1}$.

For arbitrary $\omega_{1}\neq 0$, the non-stationary Heun equation can be
written as

$\displaystyle\left(\frac{{\rm
i}\pi}{\omega_{1}^{2}}\kappa\frac{\partial}{\partial\tau}-\frac{\partial^{2}}{\partial
x^{2}}+\sum_{\nu=0}^{3}g_{\nu}(g_{\nu}-1)\wp(x+\omega_{\nu}|\omega_{1},\omega_{3})\right)\psi(x)=E\psi(x)$
(B.1)

with $\omega_{0}=0$, $\omega_{2}=-\omega_{1}-\omega_{1}\tau$,
$\omega_{3}=\omega_{1}\tau$ and $\Im(\tau)>0$. Let

$\displaystyle\psi(x)\equiv\psi\big{(}x,\tau;\omega_{1},\\{g_{\nu}\\}_{\nu=0}^{3},\kappa\big{)},\qquad
E\equiv E\big{(}\tau;\omega_{1},\\{g_{\nu}\\}_{\nu=0}^{3},\kappa\big{)}$

be a solution of this equation for generic values of $\omega_{1}$. Then this
solution can be obtained from a corresponding one for $\omega_{1}=\pi$ as
follows,

$\displaystyle\psi\big{(}x,\tau;\omega_{1},\\{g_{\nu}\\}_{\nu=0}^{3},\kappa\big{)}=\left(\frac{\pi}{\omega_{1}}\right)^{1/2}\psi\left(\pi
x/\omega_{1},\tau;\pi,\\{g_{\nu}\\}_{\nu=0}^{3},\kappa\right),$ $\displaystyle
E\big{(}\tau;\omega_{1},\\{g_{\nu}\\}_{\nu=0}^{3},\kappa\big{)}=\left(\frac{\pi}{\omega_{1}}\right)^{2}E\left(\tau;\pi,\\{g_{\nu}\\}_{\nu=0}^{3},\kappa\right)$

(to see this, use
$\wp(x|\omega_{1},\omega_{3})=s^{2}\wp(sx|s\omega_{1},s\omega_{3})$ to show
that (B.1) is invariant under the scaling transformation

$\displaystyle\big{(}x,\tau,\omega_{1},\kappa,\\{g_{\nu}\\}_{\nu=0}^{3},E\big{)}\to\big{(}sx,\tau,s\omega_{1},\kappa,\\{g_{\nu}\\}_{\nu=0}^{3},E/s^{2}\big{)},$

and set $s=\pi/\omega_{1}$; we scale $\psi(x)$ such that its $L^{2}$-norm is
invariant).

##  Appendix C Explicit solution to low orders for Algorithm I

To shorten notation we present our results using the following notation

$\displaystyle\mathcal{E}_{n}^{(\ell)}=\mathcal{E}^{(\ell)},\qquad
a^{(\ell)}(k)=\alpha_{n}^{(\ell)}(n+k)$

using that

$\displaystyle b_{n}^{(\ell)}(k)=k(k+P)-\kappa\ell\equiv b^{(\ell)}(k)$

only depends on $P$ in (3.7); see (5.2) and Remark 3.2 for further
explanations.

We give explicit formulas for the perturbative solution described in the main
text for $\ell=0,1,2$. Note that in the computations the conditions in (5.6)
are used without this being mentioned.

#### Zeroth order

For $\ell=0$ we get (here and in the following, $k=m-n$),

$\displaystyle
b^{(0)}(k)a^{(0)}(k)=\sum_{\mu=1}^{-k}\mu\gamma_{0}^{\mu}a^{(0)}(k+\mu)$ (C.1)

for $k=-1,-2,-3,\ldots$ and $a(0)=1$, which has the solution

$\displaystyle a^{(0)}(-1)=\frac{\gamma_{0}^{1}}{b^{(0)}(-1)},$ $\displaystyle
a^{(0)}(-2)=\frac{1}{b^{(0)}(-2)}\big{(}\gamma_{0}^{1}a^{(0)}(-1)+2\gamma_{0}^{0}\big{)}=\frac{1}{b^{(0)}(-2)}\left(2\gamma_{0}^{0}+\frac{(\gamma_{0}^{1})^{2}}{b^{(0)}(-1)}\right),$
$\displaystyle
a^{(0)}(-3)=\frac{1}{b^{(0)}(-3)}\big{(}\gamma_{0}^{1}a^{(0)}(-2)+2\gamma_{0}^{0}a^{(0)}(-1)+3\gamma_{0}^{1}\big{)}$
$\displaystyle\hphantom{a^{(0)}(-3)}{}=\frac{\gamma_{0}^{1}}{b^{(0)}(-3)}\left(3+2\gamma_{0}^{0}\left(\frac{1}{b^{(0)}(-2)}+\frac{1}{b^{(0)}(-1)}\right)+\frac{(\gamma_{0}^{1})^{2}}{b^{(0)}(-2)b^{(0)}(-1)}\right)$

etc.

#### First order

For $\ell=1$ we get

$\displaystyle
b^{(1)}(k)a^{(1)}(k)=\mathcal{E}^{(1)}a^{(0)}(k)+\sum_{\mu=1}^{1-k}\mu\gamma_{0}^{\mu}a^{(1)}(k+\mu)+\gamma_{1}^{1}\big{[}a^{(0)}(k+1)+a^{(0)}(k-1)\big{]}$
(C.2)

for $k=1,0,-1,\ldots$, which has the solution

$\displaystyle a^{(1)}(1)=\frac{\gamma_{1}^{1}}{b^{(1)}(1)},$
$\displaystyle\mathcal{E}^{(1)}=-\gamma_{0}^{1}a^{(1)}(1)-\gamma_{1}^{1}a^{(0)}(-1)=-\gamma_{0}^{1}\gamma_{1}^{1}\left(\frac{1}{b^{(0)}(-1)}+\frac{1}{b^{(1)}(1)}\right)$

implying the result in (5.8),

$\displaystyle
a^{(1)}(-1)=\frac{1}{b^{(1)}(-1)}\big{(}\mathcal{E}^{(1)}a^{(0)}(-1)+2\gamma_{0}^{0}a^{(1)}(1)+\gamma_{1}^{1}\big{(}1+a^{(0)}(-2)\big{)}\big{)}$
$\displaystyle=\frac{\gamma_{1}^{1}}{b^{(1)}(-1)}\left(1+2\gamma_{0}^{0}\left(\frac{1}{b^{(0)}(-2)}+\frac{1}{b^{(1)}(1)}\right)+\frac{(\gamma_{0}^{1})^{2}}{b^{(0)}(-1)}\left(\frac{1}{b^{(0)}(-2)}-\frac{1}{b^{(0)}(-1)}-\frac{1}{b^{(1)}(1)}\right)\right)$

etc.

#### Second order

For $\ell=2$ we get

$\displaystyle
b^{(2)}(k)a^{(2)}(k)=\mathcal{E}^{(2)}a^{(0)}(k)+\mathcal{E}^{(1)}a^{(1)}(k)+\sum_{\mu=1}^{2-k}\mu\gamma_{0}^{\mu}a^{(2)}(k+\mu)+\gamma_{0}^{1}\big{[}a^{(0)}(k+1)$
$\displaystyle\qquad{}+a^{(0)}(k-1)\big{]}+\gamma_{1}^{1}\big{[}a^{(1)}(k+1)+a^{(1)}(k-1)\big{]}+2\gamma_{1}^{0}\big{[}a^{(0)}(k+2)+a^{(0)}(k-2)\big{]}$
(C.3)

for $k=2,1,0,\ldots$, which has the following solution

$\displaystyle
a^{(2)}(2)=\frac{1}{b^{(2)}(2)}\big{(}\gamma_{1}^{1}a^{(1)}(1)+2\gamma_{1}^{0}\big{)}=\frac{1}{b^{(2)}(2)}\left(2\gamma_{1}^{0}+\frac{(\gamma_{1}^{1})^{2}}{b^{(2)}(1)}\right),$
$\displaystyle
a^{(2)}(1)=\frac{1}{b^{(2)}(1)}\big{(}\mathcal{E}^{(1)}a^{(1)}(1)+\gamma_{0}^{1}a^{(2)}(2)+\gamma_{0}^{1}+2\gamma_{1}^{0}a^{(0)}(-1)\big{)}=\frac{\gamma_{0}^{1}}{b^{(2)}(1)}\Biggl{(}1+2\gamma_{1}^{0}$
$\displaystyle\hphantom{a^{(2)}(1)=}{}\times\left(\frac{1}{b^{(0)}(-1)}+\frac{1}{b^{(2)}(2)}\right)+(\gamma_{1}^{1})^{2}\left(\frac{1}{b^{(2)}(2)b^{(1)}(1)}-\frac{1}{b^{(1)}(1)^{2}}-\frac{1}{b^{(1)}(1)b^{(0)}(-1)}\right)\Biggr{)},$

and

$\displaystyle\mathcal{E}^{(2)}=-\gamma_{0}^{1}a^{(2)}(1)-2\gamma_{0}^{0}a^{(2)}(2)-\gamma_{0}^{1}a^{(0)}(-1)-\gamma_{1}^{1}\big{[}a^{(1)}(1)+a^{(1)}(-1)\big{]}-2\gamma_{1}^{0}a^{(0)}(-2)$
$\displaystyle\hphantom{\mathcal{E}^{(2)}}{}=-\left(\gamma_{0}^{1}\right)^{2}\left(\frac{1}{b^{(0)}(-1)}+\frac{1}{b^{(2)}(1)}\right)-\left(\gamma_{1}^{1}\right)^{2}\left(\frac{1}{b^{(1)}(-1)}+\frac{1}{b^{(1)}(1)}\right)$
$\displaystyle\hphantom{\mathcal{E}^{(2)}=}{}-4\gamma_{0}^{0}\gamma_{1}^{0}\left(\frac{1}{b^{(0)}(-2)}+\frac{1}{b^{(2)}(2)}\right)$
$\displaystyle\hphantom{\mathcal{E}^{(2)}=}{}-2(\gamma_{0}^{1})^{2}\gamma_{1}^{0}\left(\frac{1}{b^{(0)}(-2)b^{(0)}(-1)}+\frac{1}{b^{(0)}(-1)b^{(2)}(1)}+\frac{1}{b^{(2)}(1)b^{(2)}(2)}\right)$
$\displaystyle\hphantom{\mathcal{E}^{(2)}=}{}-2\gamma_{0}^{0}\left(\gamma_{1}^{1}\right)^{2}\left(\frac{1}{b^{(0)}(-2)b^{(1)}(-1)}+\frac{1}{b^{(1)}(-1)b^{(1)}(1)}+\frac{1}{b^{(1)}(1)b^{(2)}(2)}\right)$
$\displaystyle\hphantom{\mathcal{E}^{(2)}=}{}+(\gamma_{0}^{1})^{2}(\gamma_{1}^{1})^{2}\left(\frac{1}{b^{(1)}(-1)b^{(0)}(-1)^{2}}+\frac{1}{b^{(2)}(1)b^{(1)}(1)^{2}}+\frac{1}{b^{(1)}(-1)b^{(0)}(-1)b^{(1)}(1)}\right.$
$\displaystyle\hphantom{\mathcal{E}^{(2)}=}{}+\left.\frac{1}{b^{(0)}(-1)b^{(1)}(1)b^{(2)}(1)}-\frac{1}{b^{(0)}(-2)b^{(0)}(-1)b^{(1)}(-1)}-\frac{1}{b^{(1)}(1)b^{(2)}(1)b^{(2)}(2)}\right)$

equal to the result in (5.9).

######  Remark C.1.

The results for the Heun equation ($\kappa=0$) are obtained from this by
replacing

$\displaystyle b^{(\ell)}(k)\to b^{(0)}(k)\qquad\forall\,k,\;\ell.$

For this case we found empirically that there is a useful symmetry as follows:
one can obtain the formula for $a^{(\ell)}(k)$, $k=1,2,\ldots,\ell$, from the
one for $a^{(\ell-k)}(-k)$ by the following transformation:

$\displaystyle b^{(0)}(k^{\prime})\leftrightarrow
b^{(0)}(-k^{\prime}),\qquad\gamma_{0}^{\mu}\leftrightarrow\gamma_{1}^{\mu}$

(i.e., swap the sign of the argument in all $b^{(0)}(k^{\prime})$’s, and
exchange all $\gamma_{0}^{\mu}$ and $\gamma_{1}^{\mu}$), and
$\mathcal{E}^{(\ell)}$ is invariant under this transformation $\forall\ell$.
This allows to reduce the computations considerably.

##  Appendix D Details and proofs

We collect computational details and technical proofs.

###  D.1 Proof of Lemma 2.2

We first consider the case $q=0$. We insert the definitions in (2.1)–(2.2) and
use binomial series to expand

$\displaystyle\sum_{n\in\mathbb{Z}}f_{n}^{(0)}(z)\xi^{n}=\frac{(1-\xi)^{\tilde{g}_{0}}(1+\xi)^{\tilde{g}_{1}}}{(1-2\xi
z+\xi^{2})^{\lambda}}=\sum_{\nu_{0},\nu_{1},\nu_{2}=0}^{\infty}\binom{\tilde{g}_{0}}{\nu_{0}}\binom{\tilde{g}_{1}}{\nu_{1}}\binom{-\lambda}{\nu_{2}}(-\xi)^{\nu_{0}}\xi^{\nu_{1}}\big{(}\xi^{2}-2\xi
z\big{)}^{\nu_{2}}$
$\displaystyle\hphantom{\sum_{n\in\mathbb{Z}}f_{n}^{(0)}(z)\xi^{n}}{}=\sum_{\nu_{0},\nu_{1},\nu_{2}=0}^{\infty}\sum_{m=0}^{\nu_{2}}\binom{\tilde{g}_{0}}{\nu_{0}}\binom{\tilde{g}_{1}}{\nu_{1}}\binom{-\lambda}{\nu_{2}}\binom{\nu_{2}}{m}(-1)^{\nu_{0}}\xi^{\nu_{0}+\nu_{1}+2\nu_{2}-m}(-2z)^{m}.$
(D.1)

Comparing both sides we obtain

$\displaystyle
f_{n}^{(0)}(z)=\sum_{\nu_{1},\nu_{2},m}\binom{\tilde{g}_{0}}{n+m-\nu_{1}-2\nu_{2}}\binom{\tilde{g}_{1}}{\nu_{1}}\binom{-\lambda}{\nu_{2}}\binom{\nu_{2}}{m}(-1)^{\nu_{0}}(-2z)^{m}$
(D.2)

where the sum is over all non-negative integer triples $(\nu_{1},\nu_{2},m)$
such that $n+m-\nu_{1}-2\nu_{2}\geq 0$ and $m\leq\nu_{2}$. From this it is
obvious that $f_{n}^{(0)}(z)$ is a polynomial of degree $n$ in $z$ satisfying
(2.9) (the highest power of $z$ is obtained for
$(\nu_{1},\nu_{2},m)=(0,n,n)$).

We now consider the case $q\neq 0$. We rewrite (2.6) using the result for
$q=0$ above:

$\displaystyle\sum_{n\in\mathbb{Z}}f_{n}(z)\xi^{n}=\sum_{m=0}^{\infty}f_{m}^{(0)}(z)\xi^{m}\prod_{k=1}^{\infty}\frac{\prod\limits_{\epsilon=\pm}\big{[}\big{(}1-\epsilon
q^{k}\xi\big{)}\big{(}1-\epsilon
q^{k}\xi^{-1}\big{)}\big{]}^{\tilde{g}_{k,\epsilon}}}{\big{[}\big{(}1-2q^{2k}\xi
z+q^{4k}\xi^{2}\big{)}\big{(}1-2q^{2k}\xi^{-1}z+q^{4k}\xi^{-2}\big{)}\big{]}^{\lambda}}$

with $(\tilde{g}_{k,+},\tilde{g}_{k,-})$ equal to
$(\tilde{g}_{0},\tilde{g}_{1})$ and $(\tilde{g}_{3},\tilde{g}_{2})$ for even-
and odd $k$, respectively. Expanding all factors in the product in binomial
series the r.h.s. of this becomes a linear combination of terms

$\displaystyle
f_{m}^{(0)}(z)\xi^{m}\prod_{k=1}^{\infty}\big{(}q^{k}\xi\big{)}^{\nu_{k,1}}\big{(}q^{k}\xi^{-1}\big{)}^{\nu_{k,2}}\big{[}q^{2k}\xi\big{(}q^{2k}\xi-z\big{)}\big{]}^{\nu_{k,3}}\big{[}q^{2k}\xi^{-1}\big{(}q^{2k}\xi^{-1}-z\big{)}\big{]}^{\nu_{k,4}}$

with non-negative integers $m$ and $\nu_{k,j}$, which can be further expanded
to a linear combination of terms

$\displaystyle
f_{m}^{(0)}(z)\xi^{m}\prod_{k=1}^{\infty}\big{(}q^{k}\xi\big{)}^{\nu_{k,1}}\big{(}q^{k}\xi^{-1}\big{)}^{\nu_{k,2}}\big{(}q^{2k}\xi\big{)}^{2\nu_{k,3}-m_{k,3}}z^{m_{k,3}}\big{(}q^{2k}\xi^{-1}\big{)}^{2\nu_{k,4}-m_{k,4}}z^{m_{k,4}}$
$\displaystyle\qquad{}\equiv f_{m}^{(0)}(z)q^{L}z^{M}\xi^{m+N}$

with non-negative integers $m_{k,3}\leq\nu_{k,3}$, $m_{k,4}\leq\nu_{k,4}$ and

$\displaystyle
L=\sum_{k=1}^{\infty}k(\nu_{k,1}+\nu_{k,2}+4\nu_{k,3}-2m_{k,3}+4\nu_{k,4}-2m_{k,4}),$
$\displaystyle M=\sum_{k=1}^{\infty}(m_{k,3}+m_{k,4}),\qquad
N=\sum_{k=1}^{\infty}(\nu_{k,1}-\nu_{k,2}+2\nu_{k,3}-m_{k,3}-2\nu_{k,4}+m_{k,4}).$

This implies that $f_{n}(z)$ is a linear combination of terms
$f_{m}^{(0)}(z)q^{L}z^{M}$ with $m+N=n$. We need to prove that, for all these
terms, $m+M\leq n+L$ or, equivalently,

$\displaystyle L-M+N\geq 0.$

We insert the constraints on these integers given above to write $L-M+N$ as

$\displaystyle\sum_{k=1}^{\infty}[(k+1)\nu_{k,1}+(k-1)\nu_{k,2}+2(k+1)\nu_{k,3}+2k(\nu_{k,3}-m_{k,3})+2(k-1)(2\nu_{k,4}-m_{k,4})],$

which is manifestly non-negative.

###  D.2 Series representations of $\boldsymbol{\wp}$

We derive the formulas in (4.8)–(4.9). Recall that $\omega_{\nu}$ equals $0$,
$\pi$, $-\pi-\pi\tau$ and $\pi\tau$ for $\nu=0,1,2$ and $3$, respectively.

We use the absolutely convergent series representation

$\displaystyle\frac{1}{4\sin^{2}\big{(}\frac{1}{2}z\big{)}}=-\sum_{\mu=0}^{\infty}\mu({\rm
e}^{\pm{\rm i}z})^{\mu}\qquad\mbox{if}\quad\Im(z)\gtrless 0$

and (A.4) to compute, assuming $2\pi\Im(\tau)>\Im(x)>0$ (so that
$\Im(x+2n\pi\tau)>0$ for $n\geq 0$ and $\Im(x-2n\pi\tau)<0$ for $n>0$),

$\displaystyle\wp(x|\pi,\pi\tau)+\frac{\eta_{1}}{\pi}=\sum_{n=0}^{\infty}\frac{1}{4\sin^{2}\big{(}\frac{1}{2}[x+2n\pi\tau]\big{)}}+\sum_{n=1}^{\infty}\frac{1}{4\sin^{2}\big{(}\frac{1}{2}[x-2n\pi\tau]\big{)}}$
$\displaystyle=-\sum_{n=0}^{\infty}\sum_{\mu=1}^{\infty}\mu\big{(}{\rm
e}^{{\rm
i}(x+2n\pi\tau)}\big{)}^{\mu}-\sum_{n=1}^{\infty}\sum_{\mu=1}^{\infty}\mu\big{(}{\rm
e}^{-{\rm i}(x-2n\pi\tau)}\big{)}^{\mu}$
$\displaystyle=-\sum_{\mu=1}^{\infty}\mu\left(\sum_{n=0}^{\infty}{\rm e}^{{\rm
i}\mu x}q^{2n\mu}+\sum_{n=1}^{\infty}{\rm e}^{-{\rm i}\mu
x}q^{2n\mu}\right)=-\sum_{\mu=1}^{\infty}\mu\left(\frac{1}{1-q^{2\mu}}{\rm
e}^{{\rm i}\mu x}+\frac{q^{2\mu}}{1-q^{2\mu}}{\rm e}^{-{\rm i}\mu x}\right)$
(D.3)

(we inserted (A.3) and summed up geometric series; the interchange of
summations in the third equality is justified by absolute convergence). This
implies (4.8)–(4.9) for $\nu=0$.

From this one gets (4.8)–(4.9) also in the other cases $\nu=1,2,3$: for
$2\pi\Im(\tau)>\Im(x)>0$,

$\displaystyle\wp(x+\pi|\pi,\pi\tau)=-\frac{\eta_{1}}{\pi}-\sum_{\mu=1}^{\infty}\mu(-1)^{\mu}\left(\frac{1}{1-q^{2\mu}}{\rm
e}^{{\rm i}\mu x}+\frac{q^{2\mu}}{1-q^{2\mu}}{\rm e}^{-{\rm i}\mu x}\right)$

implies the result for $\nu=1$; for $\pi\Im(\tau)>\Im(x)>-\pi\Im(\tau)$,

$\displaystyle\wp(x+\pi\tau|\pi,\pi\tau)=-\frac{\eta_{1}}{\pi}-\sum_{\mu=1}^{\infty}\mu\left(\frac{q^{\mu}}{1-q^{2\mu}}{\rm
e}^{{\rm i}\mu x}+\frac{q^{\mu}}{1-q^{2\mu}}{\rm e}^{-{\rm i}\mu x}\right)$

implies the result for $\nu=3$; and for $\pi\Im(\tau)>\Im(x)>-\pi\Im(\tau)$,

$\displaystyle\wp(x+\pi+\pi\tau|\pi,\pi\tau)=-\frac{\eta_{1}}{\pi}-\sum_{\mu=1}^{\infty}\mu(-1)^{\mu}\left(\frac{q^{\mu}}{1-q^{2\mu}}{\rm
e}^{{\rm i}\mu x}+\frac{q^{\mu}}{1-q^{2\mu}}{\rm e}^{-{\rm i}\mu x}\right)$

and the double periodicity of $\wp$ imply the result for $\nu=2$.

##  Appendix E Relation to combinatorics

In this appendix we give examples that our results allow to translate the
problem to solve the Heun equation into combinatorial problems. We also
explain in which sense the representation of the eigenvalues $\mathcal{E}_{n}$
by Theorem 6.3 is simpler than the one provided by Theorem 6.1.

We note that (6.6) implies the following formula for the Taylor coefficients
$\mathcal{E}_{n}^{(\ell)}$ of the generalized eigenvalues,

$\displaystyle\mathcal{E}_{n}^{(\ell\geq
1)}=\sum_{m}\sum_{k_{0},k_{1},\ldots,k_{m-1}}\sum_{\ell_{1},\cdots,\ell_{m}}\delta\left(\sum_{r=0}^{m-1}k_{r},m\right)\delta\left(\sum_{r=1}^{m-1}rk_{r},m-1\right)(m-1)!$
$\displaystyle\hphantom{\mathcal{E}_{n}^{(\ell\geq
1)}=}{}\times\delta\left(\ell_{1}+\cdots+\ell_{m},\ell\right)\frac{\big{[}\Phi_{n}^{(r_{1},\ell_{1})}\big{]}^{k_{1}}}{k_{1}!}\cdots\frac{\big{[}\Phi_{n}^{(r_{m},\ell_{m})}\big{]}^{k_{m}}}{k_{m}!}$
(E.1)

with

$\displaystyle\Phi^{(r,\ell)}_{n}\equiv-\sum_{s}\sum_{\mu_{1},\ldots,\mu_{s}}\sum_{r_{1},\ldots,r_{s-1}}\sum_{k_{1},\ldots,k_{s}}S_{\mu_{1}}(k_{1})\cdots
S_{\mu_{s}}(k_{s})$
$\displaystyle\hphantom{\Phi^{(r,\ell)}_{n}\equiv}{}\times\frac{\delta(0,\mu_{1}+\cdots+\mu_{s})\delta(\ell,|\mu_{1}|k_{1}+\cdots+|\mu_{s}|k_{s})\delta(r_{1}+\cdots
r_{s-1},r)}{\big{[}b^{(0)}_{n}(\mu_{1})\big{]}^{1+r_{1}}_{n}\big{[}b^{(0)}_{n}(\mu_{1}+\mu_{2})\big{]}^{1+r_{2}}_{n}\cdots\big{[}b^{(0)}_{n}(\mu_{1}+\cdots+\mu_{s-1})\big{]}^{1+r_{s-1}}_{n}},$
(E.2)

using notation introduced above (this follows by straightforward
computations); for fixed $\ell$, there are only finitely many integer vectors
$(k_{0},\ldots,k_{m-1})$ and $(\ell_{1},\ldots,\ell_{m})$ that give non-zero
contributions to the sums in (E.1) (due to the constraints in the Kronecker
deltas), and, similarly, (E.2) represents each of the finitely many different
$\Phi_{n}^{(r^{\prime},\ell^{\prime})}$ needed to evaluate (E.1) as a finite
number of terms. It is clear that this combinatorial problem is complicated.

On the other hand, if one uses (6.16) to compute $\mathcal{E}_{n}^{(\ell)}$,
one only needs $\alpha_{n}^{(\ell^{\prime})}(n)$ for
$\ell^{\prime}=1,2,\ldots,\ell$, and (6.13) for $m=n$ allows to compute this
by a combinatorial formula which is very similar to the one for
$\Phi_{n}^{(0,\ell^{\prime})}$. It thus is clear that the formula in (6.16)
gives a significantly simpler representation of $\mathcal{E}_{n}^{(\ell)}$
than the formula in (6.6), as noted already in Remark 5.3.

It is interesting to note that, from the combinatorial point of view, the
solutions of the eCS model in [24] and of the Heun equation in the present
paper are very similar; the only differences are the details of the building
blocks $\Phi^{(r,\ell)}_{n}$ of the solution. This suggests to us that all
elliptic CMS type problems have a common underlying combinatorial structure
which deserves further study.

### Acknowledgements

We thank M. Hallnäs, O. Chalykh, and H. Rosengren for helpful discussions and
comments, as well as an anonymous referee for carefully reading our paper. We
gratefully acknowledge partial financial support by the _Stiftelse Olle
Engkvist Byggmästare_ (contract 184-0573).

## References

  * [1] Bazhanov V.V., Mangazeev V.V., Eight-vertex model and non-stationary Lamé equation, J. Phys. A: Math. Gen. 38 (2005), L145–L153, hep-th/0411094.
  * [2] Chalykh O., Feigin M., Veselov A., New integrable generalizations of Calogero–Moser quantum problem, J. Math. Phys. 39 (1998), 695–703.
  * [3] Erdélyi A., Integral equations for Heun functions, Quart. J. Math., Oxford Ser. 13 (1942), 107–112.
  * [4] Etingof P.I., Frenkel I.B., Kirillov Jr. A.A., Spherical functions on affine Lie groups, Duke Math. J. 80 (1995), 59–90, hep-th/9407047.
  * [5] Etingof P.I., Kirillov Jr. A.A., Representations of affine Lie algebras, parabolic differential equations, and Lamé functions, Duke Math. J. 74 (1994), 585–614, hep-th/9310083.
  * [6] Falceto F., Gawȩdzki K., Unitarity of the Knizhnik–Zamolodchikov–Bernard connection and the Bethe ansatz for the elliptic Hitchin systems, Comm. Math. Phys. 183 (1997), 267–290, hep-th/9604094.
  * [7] Fateev V.A., Litvinov A.V., Neveu A., Onofri E., A differential equation for a four-point correlation function in Liouville field theory and elliptic four-point conformal blocks, J. Phys. A: Math. Theor. 42 (2009), 304011, 29 pages, arXiv:0902.1331.
  * [8] Felder G., Varchenko A., Integral representation of solutions of the elliptic Knizhnik–Zamolodchikov–Bernard equations, Int. Math. Res. Not. 1995 (1995), 221–233, hep-th/9502165.
  * [9] Felder G., Varchenko A., Special functions, conformal blocks, Bethe ansatz and ${\rm SL}(3,{\mathbb{Z}})$, Phil. Trans. R. Soc. Lond. A 359 (2001), 1365–1373, math.QA/0101136.
  * [10] Hallnäs M., Langmann E., A unified construction of generalized classical polynomials associated with operators of Calogero–Sutherland type, Constr. Approx. 31 (2010), 309–342, math-ph/0703090.
  * [11] Hallnäs M., Ruijsenaars S., A recursive construction of joint eigenfunctions for the hyperbolic nonrelativistic Calogero–Moser Hamiltonians, Int. Math. Res. Not. 2015 (2015), 10278–10313, arXiv:1305.4759.
  * [12] Inozemtsev V.I., Lax representation with spectral parameter on a torus for integrable particle systems, Lett. Math. Phys. 17 (1989), 11–17.
  * [13] Kazakov A.Y., Slavyanov S.Yu., Integral relations for Heun-class special functions, Theoret. and Math. Phys. 107 (1996), 733–739.
  * [14] Kolb S., Radial part calculations for $\widehat{\mathfrak{sl}}_{2}$ and the Heun-KZB heat equation, Int. Math. Res. Not. 2015 (2015), 12941–12990, arXiv:1310.0782.
  * [15] Komori Y., Noumi M., Shiraishi J., Kernel functions for difference operators of Ruijsenaars type and their applications, SIGMA 5 (2009), 054, 40 pages, arXiv:0812.0279.
  * [16] Komori Y., Takemura K., The perturbation of the quantum Calogero–Moser–Sutherland system and related results, Comm. Math. Phys. 227 (2002), 93–118, math.QA/0009244.
  * [17] Koroteev P., Sciarappa A., On elliptic algebras and large-$n$ supersymmetric gauge theories, J. Math. Phys. 57 (2016), 112302, 32 pages, arXiv:1601.08238.
  * [18] Lambe C.G., Ward D.R., Some differential equations and associated integral equations, Quart. J. Math. os-5 (1934), 81–97.
  * [19] Langmann E., Anyons and the elliptic Calogero-Sutherland model, Lett. Math. Phys. 54 (2000), 279–289, math-ph/0007036.
  * [20] Langmann E., Algorithms to solve the (quantum) Sutherland model, J. Math. Phys. 42 (2001), 4148–4157, math-ph/0104039.
  * [21] Langmann E., Remarkable identities related to the (quantum) elliptic Calogero–Sutherland model, J. Math. Phys. 47 (2006), 022101, 18 pages, math-ph/0406061.
  * [22] Langmann E., Singular eigenfunctions of Calogero–Sutherland type systems and how to transform them into regular ones, SIGMA 3 (2007), 031, 18 pages, math-ph/0702089.
  * [23] Langmann E., Source identity and kernel functions for elliptic Calogero–Sutherland type systems, Lett. Math. Phys. 94 (2010), 63–75, arXiv:1003.0857.
  * [24] Langmann E., Explicit solution of the (quantum) elliptic Calogero–Sutherland model, Ann. Henri Poincaré 15 (2014), 755–791, math-ph/0407050.
  * [25] Langmann E., Takemura K., Source identity and kernel functions for Inozemtsev-type systems, J. Math. Phys. 53 (2012), 082105, 19, arXiv:1202.3544.
  * [26] Maier R.S., The 192 solutions of the Heun equation, Math. Comp. 76 (2007), 811–843, math.CA/0408317.
  * [27] Nagoya H., Hypergeometric solutions to Schrödinger equations for the quantum Painlevé equations, J. Math. Phys. 52 (2011), 083509, 16 pages, arXiv:1109.1645.
  * [28] Nawata S., Givental $J$-functions, quantum integrable systems, AGT relation with surface operator, Adv. Theor. Math. Phys. 19 (2015), 1277–1338, arXiv:1408.4132.
  * [29] Nekrasov N.A., Shatashvili S.L., Quantization of integrable systems and four dimensional gauge theories, in XVIth International Congress on Mathematical Physics, World Sci. Publ., Hackensack, NJ, 2010, 265–289, arXiv:0908.4052.
  * [30] NIST digital library of mathematical functions, Release date 2010-05-07, available at http://dlmf.nist.gov/.
  * [31] Novikov D.P., Integral transformation of solutions for a Fuchsian-class equation corresponding to the Okamoto transformation of the Painlevé VI equation, Theoret. and Math. Phys. 146 (2006), 295–303.
  * [32] Ronveaux A., Arscott F. (Editors), Heun’s differential equations, Oxford Science Publications, Oxford University Press, Oxford, 1995.
  * [33] Rosengren H., Special polynomials related to the supersymmetric eight-vertex model. III. Painlevé VI equation, arXiv:1405.5318.
  * [34] Rosengren H., Special polynomials related to the supersymmetric eight-vertex model: a summary, Comm. Math. Phys. 340 (2015), 1143–1170, arXiv:1503.02833.
  * [35] Ruijsenaars S.N.M., Hilbert–Schmidt operators vs. integrable systems of elliptic Calogero–Moser type. III. The Heun case, SIGMA 5 (2009), 049, 21 pages, arXiv:0904.3250.
  * [36] Sergeev A.N., Calogero operator and Lie superalgebras, Theoret. and Math. Phys. 131 (2002), 747–764.
  * [37] Slavyanov S.Yu., Lay W., Special functions. A unified theory based on singularities, Oxford Mathematical Monographs, Oxford University Press, Oxford, 2000.
  * [38] Smirnov A.O., Elliptic solitons and Heun’s equation, in The Kowalevski Property (Leeds, 2000), CRM Proc. Lecture Notes, Vol. 32, Amer. Math. Soc., Providence, RI, 2002, 287–305, math.CA/0109149.
  * [39] Takemura K., The Heun equation and the Calogero–Moser–Sutherland system. II. Perturbation and algebraic solution, Electron. J. Differential Equations 2004 (2004), 15, 30 pages, math.CA/0112179.
  * [40] Takemura K., The Heun equation and the Calogero–Moser–Sutherland system. IV. The Hermite–Krichever ansatz, Comm. Math. Phys. 258 (2005), 367–403, math.CA/0406141.
  * [41] Whittaker E.T., Watson G.N., A course of modern analysis. 4th ed., Cambridge Mathematical Library, Cambridge University Press, Cambridge, 1996.