Series Solutions of the Non-Stationary Heun Equation

We consider the non-stationary Heun equation, also known as quantum Painlev\'e 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.


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 equation 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].

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 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 1 To be more precise: the special case of the KZB equations corresponding to g = sl2 and n = 1.
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) andH(y) is an explicitly known function K(x, y) satisfying the functional identity (H(x) −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 =H) or different (H = 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 socalled balancing condition κ = 0 with κ 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, τ ) satisfying with τ the half-period ratio of the elliptic functions appearing in the CMS Hamiltonians H andH [21,23,25] (see, e.g., Lemma 4.1 where the τ -dependence of H,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 i π κ ∂ ∂τ +H(y, τ ) to an eigenfunction of i π κ ∂ ∂τ + H(x, τ ) [25].

Summary of results
The non-stationary Heun equation can be written as 2 with ℘(x) the Weierstrass elliptic function with periods (2π, 2πτ ) and (for the convenience of the reader we collect the definitions of ℘ and other well-known special functions we need in Appendix A). To simplify notation we set ω 1 = π here and in most parts of this paper; see Appendix B for how to transform our results to other values of ω 1 . The parameters g 0 , g 1 , g 2 , g 3 and κ can be arbitrary complex numbers for our general results. 3 Our aim is to construct functions ψ(x) ≡ ψ(x, τ ; {g ν } 3 ν=0 , κ) of two complex variables x and τ , (τ ) > 0, that satisfy this differential equation for some E ≡ E(τ ; {g ν } 3 ν=0 , κ); a more precise characterization of our solutions is given in (1.5)-(1.6) below. It is important to note that, for κ = 0, E can be transformed to 0, or any other convenient value, by changing the normalization of ψ(x) (this follows from the obvious invariance of (1.2) under the transformation for arbitrary analytic functions C of τ ). However, we sometimes find it convenient to impose a normalization condition on ψ(x) so that E remains significant even for κ = 0. Important special cases of (1.2) include the Heun equation (κ = 0), the non-stationary Lamé equation (g ν = g independent of ν, or g ν (g ν − 1) = 0 for three of the ν'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 (τ ) → +∞, (1.2) reduces to the stationary Schrödinger equation with the Pöschl-Teller potential ∝ 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 ψ n (x) = (2q 1/4 ) −g 0 −g 1 3 ν=0 θ ν+1 1 2 x gν P n (cos(x)) , 2 Here and in the following we often suppress the dependence of functions on the variable τ . 3 We 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.
with θ ν+1 (z) the Jacobi theta functions, q = exp(iπτ ) the nomé, η 1 /π in (A.5) and with and P g 0 − 1 2 ,g 1 − 1 2 n (z) the Jacobi polynomials (see (A.11)). Our main result provides efficient recursive procedures to compute the functions P ( ) n (z) which, as we show, are polynomials of degree n+ in z; see Propositions 5.1 and 5.2 for two complementary variants of this result. Thus one can regard 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 → 0. We conjecture that the series in (1.6) is absolutely convergent and converges to a L 2 -function on [0, π] for |q| ≤ q 0 and some q 0 > 0 depending on parameters (this is known to be true in the Heun case κ = 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 ψ n (x) provide a complete orthonormal basis in the Hilbert space of L 2 -functions on [0, π] 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 P ( ) n (z) and E ( ) n of the power series in (1.6). We note that, in the Heun case κ = 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 ψ 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 κ = 0, whereas the first applies to the case κ = 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 κ → 0.
One important feature of our method is that it provides particular q-dependent basis functions {f m (z)} m∈Z to expand the functions 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 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 κ = 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 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].

Plan
Section 2 contains preliminary material: a summary of notation (Section 2.1), a review of a well-known solution of (1.2) for (τ ) → ∞ 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 ω 1 = π to other values of ω 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).

Preliminaries
We collect definitions and preliminary results that we use.

Notation
We use the special functions (ξ and q are complex variables; |q| < 1) which all are closely related to the Jacobi theta functions (see (A.7) and (A.8)). We denote as N 0 and Z the sets of non-negative and non-zero integers, respectively. The symbol δ(m, n) for integers m, n denotes the Kronecker delta.

Trigonometric limit
The non-stationary Heun equation simplifies in the trigonometric case q = 0 to which is known to have solutions with the Jacobi polynomials P (α,β) n (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.
n (x) be a solution of (2.3) of the form for some constant E (0) , with P (z) a polynomial of degree n such that for some n ∈ N 0 . Then P (z) = P 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.

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, with the special functions Θ ν+1 (ξ), Θ(z, ξ) 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. 4 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).
have the following power series expansion In particular, with the binomial coefficient −λ m as usual (see (A.10)).
We will also use the following integral representation of these functions: with C the contour once around the circle with radius |ξ| = R, |q| < R < 1, taken counterclockwise (this is equivalent to (2.6) by Cauchy's theorem).
, only in this case the binomial coefficient on the r.h.s. in (2.9) is always non-zero). If −λ = k ∈ N 0 , one can see from (2.10) that f (0) m (z) is a polynomials of degree ≤ k for all m ∈ N 0 . Thus, to get complete results, we sometimes impose the condition −λ / ∈ N 0 .

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 (κ) = 0 and g 0 + g 1 ∈ R, or κ = 0 and g 0 + g 1 / ∈ 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, n in (2.4) the energy eigenvalues of the unperturbed problem q = 0 and = 0, 1, 2, . . .; we refer to a case (m, ) 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 κ = 0 and for fixed n ∈ N 0 , the results for the expansion coefficients P ( ) n (z) and E ( ) n in Proposition 5.1 and Theorem 6.1 are valid even in the limit (The proof is given at the end of this section.) Thus, in the Heun case κ = 0 and for all n ∈ N 0 , our results can be extended to the case g 0 + g 1 ∈ 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 κ. 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 κ = 0 is not known (to our knowledge).
Proof of Lemma 2.4. For κ = 0 and fixed n ∈ N 0 , the energy denominator in (2.11) appearing in standard perturbation theory [39] are only for m ∈ 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 = n. However, it is clear that our results for the coefficients P ( ) n (z) and E ( ) n 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 , κ = 0. The interested reader is referred to [24] for a more extensive discussion of resonances.

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 differentialdifference 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).

Dif ferential-dif ference equation
We construct solutions ψ 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 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), with the raising Pochhammer symbol (x) n in (A.9) and λ in (2.7); note that N n is finite and non-zero for all integers n if −λ / ∈ N 0 for n > 0 and −(g 0 + g 1 ) / ∈ N 0 for n < 0 (see the discussion after (A.9)).
Our key result is equations determining α n (m) and E n and which, as we will show, can be solved efficiently. To state this result we introduce the convenient shorthand notation k+2r for all integers r, s). Proposition 3.1. Let n ∈ Z, −λ / ∈ N 0 for n > 0 and −(g 0 + g 1 ) / ∈ N 0 for n < 0, and assume that E n and α n (m) for m ∈ Z satisfy the following system of equations, and the condition α n (m)| q=0 = 0, m > n, 1, m = n. (The proof is given in Section 4.) It is important to note that the conditions above do not determine P n (z) and E n uniquely: 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 α n (m) and E n of the equations in Proposition 3.1 are essentially independent of n in the following sense: they are of the form with functions a(k) andẼ depending on n only in the combination (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.

Explicit solutions by integrals
The solutions α n (m), E n of (3.4)-(3.5) are complicated in general. However, there exist nontrivial cases where Proposition 3.1 gives simple explicit solutions of the non-stationary Heun equation: , and C the integration contour defined after (2.10). Then the non-stationary Heun equation in (1.2) has solutions ψ n (x), E n as in (1.5) with and such that the conditions in Proof . The assumptions imply γ µ k = 0 for all k, µ, and the equations in (3.4)-(3.5) in this case have the solution α n (m) = δ(m, n), E n = E (0) n . Proposition 3.1 implies the result.
Note that (3.8) gives several different one-parameter families ({g ν } 3 ν=0 , κ), depending on λ, where this result provides simple integral representations of elliptic Jacobi polynomials 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, for n ∈ N 0 (the latter identities are obtained from Corollary 3.3 for the cases where (g 0 ,g 1 , λ) 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].  [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 κ 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  We emphasize that the result in Corollary 3.3 is more general in that it includes some nonstationary 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,Ñ ) = (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.

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).

Kernel function method
We introduce the notation with ω ν in (1.3). This allows us to write the non-stationary Heun equation in (1.2) as We also recall the definitions ofg ν and λ 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: obeys the identity Proof . This is the special case N = M = 1 of Corollary 3.2 in [25] (note that the symbols β and A 1,1 in [25] correspond to −2πiτ and 2κ here, respectively).
We use this to compute the action of i π κ ∂ ∂τ + H(x; {g ν } 3 ν=0 ) on the functions with f m (z) defined in (2.6). For this we note that Θ ν+1 (e iy )g ν Θ(cos x, e iy ) λ e iy(g 0 +g 1 )/2 with G defined in (A.6) (we used (A.7) and (A.8)). This and (2.6) show that K(x, y) is a generating function for the functions in (4.6): To evaluate H(y; {g ν } 3 ν=0 )K(x, y) we use the following expansions (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 (4.12) (The proof can be found at the end of this section.) To proceed we make the ansatz with N a q-independent normalization constant to be determined. Inserting this in (1.2) and using Lemma 4.2 we obtain It follows that the function in (4.13) is a solution of (1.2) if the coefficients α(m) and E satisfy (4.14) Inserting (4.9) and (4.11) and changing variables from τ to q = e iπτ allow us to write this as equivalent to (3.4). This proves that ψ(x) in (4.13) and E = E + C 0 solve (1.2) provided (4.15) is fulfilled.
We are left to show that a solution α(m) = α n (m), E = E n satisfying the condition in (3.5) is such that (1.7) holds true. This is done in Section 4.2.

Trigonometric limit
For q = 0 the equations in (3.4)-(3.5) simplify to and α (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: n (n) = 1, and provided that g 0 + g 1 is such that The latter condition is implied by our assumption that −(g 0 + g 1 ) / ∈ N.
We recall that f (0) m (z) ≡ f m (z)| q=0 is zero for m < 0 and a polynomial satisfying (2.9) for m ≥ 0 (see Lemma 2.2), and thus is a polynomial of degree n. Moreover, the special case q = 0 of results proved above implies that ψ 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 N n as in (3.2) and completes the proof of Proposition 3.1.

Recursive algorithms
We now present algorithms to compute the expansion coefficients P ( ) n (z) and E ( ) n 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 κ, including κ = 0. We then give a second algorithm, which is a variant of the first, which is simpler but requires κ = 0; see Section 5.2.

Algorithm I
We compute the functions P 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 P ( ) n (z) which make manifest that they are polynomials of degree n + in z for n + ≥ 0 and zero otherwise.
To formulate this result we find it useful to introduce the shorthand notation n in (2.4)) and recall the definition of λ in (2.7). Note that the denominators in the fractions in (2.11) are equal to b ( ) n (m − 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 α This computation can be straightforwardly implemented in a symbolic programming language like MAPLE or MATHEMATICA and, in this way, extended to higher values of . 6 As proved by Ruijsenaars in [35], the Heun equation in (1.2) for κ = 0 has solutions ψ n (x) and corresponding eigenvalues E n such that the latter are invariant under all permutations of the following affine combinations of coupling parameters, We used MAPLE to check that the coefficients E

Algorithm II
For κ = 0, the condition α  6 We computed E ( ) n up to = 6 using ordinary laptops and with computation times of the order of minutes.
Proposition 5.2. Let n ∈ Z, −λ / ∈ N 0 for n > 0 and −(g 0 +g 1 ) / ∈ N 0 for n < 0, f ( ) m (z) m∈Z the functions defined and characterized in Lemma 2.2, and assume that (κ) = 0 and g 0 + g 1 ∈ R. 7 Then the non-stationary Heun equation in (1.2) has a unique solution ψ n (x), E n as in (1.5)-(1.6) given by (3.2), (5.3) and with coefficients α ( ) n (m) determined by the following recursion relations, for all m ∈ Z if = 0 and all m ∈ Z if ≥ 1, together with the following conditions The functions P ( ) n (z) thus obtained are polynomials of degree n + in z for n + ≥ 0 and zero otherwise.
It follows from (5.1) that the first and last conditions in (5.6) are equivalent to α I n (n) = 1, which give rise to E I n = E , κq ∂ ∂q α II n (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 α II n (n) up to some order O(q N +1 ), one obtains a Padé approximation of the generalized eigenvalues E I n as follows, As an example we give n and E (2) n 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 κ = 0 in the numerator and denominator cancel. From our results above it is clear that this cancellation of singularities at κ = 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 κ → 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.

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 ∈ Z, (1.2) has solutions ψ n (x), E n as in (1.5) and (3.1)-(3.2) provided that the coefficients α n (m) and the (shifted) eigenvalue E n are (suitable) solutions of the differential-difference equations in (3.4). We observe that, by introducing the notation 9 m in (2.4) and S µ in (4.11), one can write (3.4) in a simple way as follows, (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 κ = 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 κ = 0; Theorem 6.3 generalizes Proposition 5.2.

Heun equation
The Heun equation corresponds to the special case κ = 0 of (1.2). This is an eigenvalue equation for a 1D Schrödinger operator, and ψ(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 κ = 0 to all orders [24]. To formulate this result we introduce the following two functions of a complex variable z (n ∈ N 0 and m ∈ Z are fixed), with the Kronecker delta δ(n, m), S µ in (4.11), and the convenient shorthand notations It is important to note that these formulas should be interpreted in a perturbative sense as follows: for fixed N ∈ N, there are only finitely many terms in (6.3) and (6.4) that are O(q ) with < 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 ≤ s ≤ N and −N ≤ µ j ≤ N for j = 1, 2, . . . , 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 ∈ Z, −λ / ∈ N 0 for n > 0 and −(g 0 + g 1 ) / ∈ N 0 for n < 0, f ( ) m (z) m∈Z the functions defined and characterized in Lemma 2.2. Then the Heun equation in (1.2) for κ = 0 and g 0 + g 1 / ∈ Z has a unique solution as in (1.5) and (3.1)-(3.2) provided that E n = E (0) n +Ẽ n , α n (m) = G n Ẽ n ; m withẼ n the unique solution of the equatioñ (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, µs∈Z r 1 ,...,r s−1 ∈N 0 S µ 1 · · · S µs × δ(0, µ 1 + · · · + µ s )δ(r 1 + · · · r s−1 , r) and similarly for α 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 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 g 0 − 1 2 ,g 1 − 1 2 n (z) provided by Theorem 6.1 can be formally written as It is possible to identifyP n (ξ) with a singular eigenfunction of the Inozemtsev Hamiltonian H(y; {g ν } 3 ν=0 ) 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 {α(m)} m∈Z and regard A and B in (6.1) as linear operators V → V . Define projections P on V as follows, 10 (Pα)(m) ≡ δ(n, m)α(m) ∀ α ∈ V (6.8) and P ⊥ ≡ I − P. Then A in (6.1) commutes with P, and the equation Pα 0 = α 0 is solved by (α 0 )(m) ≡ δ(n, m). Thus Lemma 4.1.1 in [24] implies that is a solution of (6.2). Expanding the resolvent in a Neumann series we obtain [24] n (m − n) −Ẽ n α(m) (6.11) using the shorthand notation in (6.5). Using (6.11) to compute (6.10) we obtain α(m) = G n (Ẽ; m) andẼ = Φ n (Ẽ) with the functions defined in (6.3)-(6.4) [24]. The latter should be interpreted as non-linear equation whose solutionẼ =Ẽ n vanishing like O(q) determines the solution of the Heun equation we are interested in; see [24] for details.

Time dependent Heun equation
We now present generalizations of the results in the previous section to the non-stationary case κ = 0. The argument to solve (6.1)-(6.2) in the proof of Theorem 6.1 can be generalized to κ = 0 if the following projection is used, where α ( ) (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 E = E using the shorthand notations with γ µ k in (3.3), and (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 E n of the time dependent Heun equation: with α n (m) in (6.13) for m = n; the limit κ → 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 κ = 0: the formula in (6.6), and the limit κ → 0 of the formula in (6.16). The former has the advantage that it is manifestly finite for κ = 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 α = {α ( ) (m)} ∈N 0 ,m∈Z identified with α = {α(m)} m∈Z as in (5.1). Then A and B in (6.1) can be written as linear operators V → V as follows, for µ ∈ Z and ∈ 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 P in (6.12) commutes with A in (6.17), and that the equation Thus Lemma 4.1.1 in [24] implies that the solution α and E of (3.4) satisfies (6.9) with that α 0 (P ⊥ = I − P). The definitions imply We showed already in Section 4.2 that the solution α of (3.4) is such that α (0) (m > n) = 0, and thus PBα = 0. With that we find that the second equation in (6.9) is solved by E = E (0) n . 12 Using this the first equation in (6.10) simplifies to Inserting definitions we find using the notation in (6.15). With that one can prove by induction that for all s = 1, 2, . . .. Inserting this in (6.20) and using (6.19) give (to simplify notation we extended the 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 ( ) µ in (6.18) suggests to change summation variables from to k such that = |µ|k (to reduce the number of terms in the formula which give zero contributions). Writing S µ (k) short for S (|µ|k) µ we obtain the result in (6.13)-(6.14).

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 ≥ 2 [12], which defines a natural many-variable generalization of the non-stationary Heun equation: note that generalized kernel function identities for this problem and a discrete set of κ-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 κ = 0, this symmetry can be exploited to construct simpler representations of the generalized eigenvalues E which are useful even in the case κ = 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: 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.

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

B Scaling
In the literature different conventions for ω 1 are used. We therefore state how our results for ω 1 = π can be transformed into results for other values of ω 1 . For arbitrary ω 1 = 0, the non-stationary Heun equation can be written as be a solution of this equation for generic values of ω 1 . Then this solution can be obtained from a corresponding one for ω 1 = π as follows, (to see this, use ℘(x|ω 1 , ω 3 ) = s 2 ℘(sx|sω 1 , sω 3 ) to show that (B.1) is invariant under the scaling transformation and set s = π/ω 1 ; we scale ψ(x) such that its L 2 -norm is invariant).
C Explicit solution to low orders for Algorithm I To shorten notation we present our results using the following notation 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 = 0, 1, 2. Note that in the computations the conditions in (5.6) are used without this being mentioned.

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 E n by Theorem 6.3 is simpler than the one provided by Theorem 6.1.
On the other hand, if one uses (6.16) to compute E ( ) n , one only needs α ( ) n (n) for = 1, 2, . . . , , and (6.13) for m = n allows to compute this by a combinatorial formula which is very similar to the one for Φ (0, ) n . It thus is clear that the formula in (6.16) gives a significantly simpler representation of E ( ) n 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 Φ (r, ) n of the solution. This suggests to us that all elliptic CMS type problems have a common underlying combinatorial structure which deserves further study.