text:
# Generalized Lamé operators

## 1. Introduction

In this paper we consider higher-dimensional analogues of the classical Lamé operator

(1.1)
$$L=\frac{d^2}{dz^2}+m(m+1)\mathrm{}(z),m_+.$$

Here $`\mathrm{}(z)=\mathrm{}(z|1,\tau )`$ is the Weierstrass $`\mathrm{}`$-function with periods $`1,\tau `$. More generally, we are interested in multivariable analogues of the so-called elliptic algebro-geometric operators $`L=d^2/dz^2+u(z)`$, which appeared in the finite-gap theory initiated in 70’s by Novikov . This theory provides a beautiful interplay between the spectral theory and algebraic geometry, and the Lamé operator is the simplest and best known member of this family of operators (see  for a survey). Since then there have been several attempts to generalize some parts of that theory to higher dimensions, most notably , see also . We should stress, however, that in general this leads to differential operators with matrix coefficients. On the other hand, in  it was suggested to consider the quantum elliptic Calogero–Moser problem and its versions related to the root systems  as natural multidimensional analogues of the Lamé operator. More specifically, a conjecture from  says that for integer values of the coupling parameters the corresponding Schrödinger operators are algebraically integrable (this is a proper generalization of the properties of the algebro-geometric operators to higher dimensions, see  and Section 3 below). For the rational and trigonometric versions of the Calogero–Moser problem this was proved in . Elliptic version, however, turned out to be more difficult: until now it was known for $`A_n`$ case only, due to .

One of the results of the present paper is a proof of that conjecture of  for all root systems. In fact, our approach applies to a wider class of Schrödinger operators, which is an elliptic version of the class introduced in , see also . Their singularities are the second order poles along a set of hyperplanes satisfying some special conditions, which encode the triviality of the local monodromy around each of the poles. We call the corresponding Schrödinger operators the generalized Lamé operators. The Calogero–Moser operators with integer coupling parameters give particular examples of such operators. Our main result says that for a generalized Lamé operator its algebraic integrability follows from the usual, Liouville integrability (in a slightly stronger sense). The proof uses a criterion from , based on differential Galois theory. In dimension one we recover in this way the main result of . For the elliptic Calogero–Moser problem the complete integrability was proved (for all root systems) by Cherednik , and this allows us to prove the conjecture of .

A complete description of all generalized Lamé operators is an open problem. All known (irreducible) examples in dimension $`>1`$ are related to the root systems and their deformations which appeared in ; we list them all in Section 4. We conjecture that they are all algebraically integrable. Using our main result, we check this for all two-dimensional examples, since the complete integrability is relatively easy to work out in that case.

One important property of the algebraically integrable operators is that their eigenfunctions can be calculated explicitly (at least, in principle). In section 5 we explain how to find the Bloch eigenfunctions for a given integrable generalized Lamé operator, in particular, for the elliptic Calogero–Moser problem. As a result, we will see that the Bloch solutions are parametrized by the points of an algebraic variety, which is a covering of a product of elliptic curves (in a perfect agreement with the situation in dimension one, due to Krichever ). Let us mention that for the elliptic Calogero–Moser problem (in the $`A_n`$ case) the Bloch eigenfunctions were calculated by Felder and Varchenko . Our procedure is different and more general (at cost of being less effective). We also explain how these Bloch solutions can be used to construct the discrete spectrum eigenstates for the Calogero–Moser problem.

In the last three sections of the paper we consider two particular examples of the generalized Lamé operators in dimension $`2`$, for which we make the formulas for the Bloch solutions very explicit. The first example is

(1.2)
$$L=\mathrm{\Delta }+2\mathrm{}(x)+2\mathrm{}(y)+4\mathrm{}(xy)+4\mathrm{}(x+y).$$

This is a special case of the elliptic Calogero–Moser problem of the $`B_2`$-type. Our second example is

(1.3)
$$L=\mathrm{\Delta }+2\mathrm{}(x)+2(a^2+b^2)\mathrm{}(ax+by)+2(\overline{a}^2+\overline{b}^2)\mathrm{}(\overline{a}x+\overline{b}y),$$

$$a=\frac{1+i\sqrt{3}\mathrm{sin}\alpha }{2},\overline{a}=\frac{1i\sqrt{3}\mathrm{sin}\alpha }{2},b=\overline{b}=\frac{i\sqrt{3}\mathrm{cos}\alpha }{2},$$

where $`\alpha `$ is a complex parameter. Such operator was considered by Hietarinta  who showed that it admits a commuting operator of order $`3`$. Its rational version $`\mathrm{}(z)=z^2`$ corresponds to a specific choice of the parameters in the family of two-dimensional Schrödinger operators introduced by Berest and Lutsenko  in connection with Huygens’ Principle (see section 4 of  for details). Note that the potential in (1.3) is real-valued (for real $`x,y`$) if $`\alpha `$ is real and the period $`\tau `$ is pure imaginary.

It is more convenient to work with the following $`3`$-dimensional version of (1.3):

(1.4)
$$\begin{array}{c}L=_1^2_2^2_3^2+2(a_1^2+a_2^2)\mathrm{}(a_1x_1a_2x_2)\hfill \\ \hfill +2(a_2^2+a_3^3)\mathrm{}(a_2x_2a_3x_3)+2(a_3^2+a_1^2)\mathrm{}(a_3x_3a_1x_1),\end{array}$$

where $`a_1^2+a_2^2+a_3^2=0`$, $`_i=/x_i`$. Then it is easy to see that $`L`$ commutes with the operator

$$L_0=a_1^1_1+a_2^1_2+a_3^1_3$$

and after restriction to the plane

$$a_1^1x_1+a_2^1x_2+a_3^1x_3=0$$

it reduces to the operator (1.3) with proper $`a,b`$.

We calculate explicitly the Bloch eigenfunctions of the operators (1.2), (1.4). Notice a certain similarity between our approach and the one used by Inozemtsev for $`A_2`$ case . Let us also mention that in dimension two there is a nice theory of Schrödinger operators which are finite-gap at a fixed energy level, see . It would be interesting to analyze our results from that point of view.

In the last two sections we also calculate the Bloch solutions for the discrete versions of (1.2), (1.4), which are given by certain difference operators of Macdonald–Ruijsenaars type. This raises a natural question about generalizing our results to the difference setting. We hope to return to this problem in future.

Acknowledgments. We are grateful to A.P.Veselov for stimulating discussions and to Yu.Berest, S.Ruijsenaars and K.Takemura for useful comments. P.E. thanks C. De Concini for a discussion which was useful for the proof of Theorem 3.8. The work of O.C. was supported by EPSRC. The work of A.O. was supported by Russian Foundation for Basic Research (grant RFBR-01-01-00803). The work of P.E. was partially supported by the NSF grant DMS-9988796, and partly done for the Clay Mathematics Institute.

## 2. Generalized Lamé operators

Let $`V=^n`$ be a complex Euclidean space with the scalar product denoted by $`(,)`$, and $`𝒜=\{\alpha \}`$ be a given finite set of affine-linear functions on $`V`$. Let us consider a Schrödinger operator

(2.1)
$$L=\mathrm{\Delta }+u(x),\mathrm{\Delta }=^2/x_1^2+\mathrm{}+^2/x_n^2,$$

with the elliptic potential $`u`$ of the following form:

(2.2)
$$u(x)=\underset{\alpha 𝒜}{}c_\alpha \mathrm{}(\alpha (x)|\tau ).$$

Here $`\mathrm{}(z|\tau )`$ is the Weierstrass $`\mathrm{}`$-function with the periods $`1,\tau `$, $`\mathrm{Im}(\tau )>0`$, and $`c_\alpha `$ are the parameters which will be specified later (below we will mostly suppress $`\tau `$, denoting $`\mathrm{}(z|\tau )`$ simply by $`\mathrm{}(z)`$). Each of the functions $`\alpha (x)`$ in standard coordinates on $`V`$ looks as

$$\alpha (x)=a_0+a_1x_1+\mathrm{}+a_nx_n,$$

so $`\alpha `$ is, effectively, a pair $`(\alpha _0,a_0)`$ where $`a_0=\alpha (0)`$ and $`\alpha _0=\mathrm{grad}\alpha =(a_1,\mathrm{},a_n)`$. We assume that each $`\alpha _0`$ is non-isotropic (co)vector, i.e. $`(\alpha _0,\alpha _0)=a_1^2+\mathrm{}+a_n^20`$. Let us project $`𝒜`$ onto $`V^{}`$, $`\alpha \alpha _0`$, denoting by $`𝒜_0`$ the resulting set.

We want $`u`$ to be periodic, more precisely, we assume that the lattice $`V^{}`$, generated over $``$ by the set $`𝒜_0`$, has rank $`n`$. For simplicity, let us assume that $`\mathrm{rk}=n`$ (the general case can be reduced to that by passing to the factor space $`V/\mathrm{Ann}`$). In that case the lattice $`+\tau `$ is the period lattice for $`u`$, with $`:=\mathrm{Hom}(,)V`$. Thus, $`u`$ may be considered as a meromorphic function on a (compact) torus $`T=V/+\tau `$, which is isomorphic to the product of $`n`$ copies of the elliptic curve $`=/+\tau `$. Singularities of $`u(x)`$ are the second order poles along the following set of the hyperplanes:

$$Sing=\underset{\alpha 𝒜}{}\underset{m,n}{}\pi _\alpha ^{m,n},\pi _\alpha ^{m,n}:=\{x:\alpha (x)=m+n\tau \}.$$

We will assume that all hyperplanes $`\pi _\alpha ^{m,n}`$ are pairwise different; this can always be achieved by rearranging the terms in (2.2). This will imply that

$$u(x)c_\alpha (\alpha (x)mn\tau )^2+O(1)\text{near}\pi _\alpha ^{m,n}.$$

Our next important assumption is that each $`c_\alpha `$ in (2.2) has a form

$$c_\alpha =m_\alpha (m_\alpha +1)(\alpha _0,\alpha _0)\text{with some}m_\alpha _{>0}.$$

Now we are going to put some more restrictions on $`u`$ demanding its quasi-invariance in the following sense.

###### Definition.

Let us say that the potential (2.2) as above is quasi-invariant if for any hyperplane $`\pi =\pi _\alpha ^{m,n}Sing`$ the (meromorphic) function $`u(x)u(s_\pi x)`$ is divisible by $`(\alpha (x)mn\tau )^{2m_\alpha +1}`$, where $`s_\pi `$ denotes the orthogonal reflection with respect to $`\pi `$.

Here are two important examples of such potentials (more examples will appear later).

###### Example 2.1.

Consider the following potential $`u(x)`$ in $`^n`$:

(2.3)
$$u=\underset{i<j}{}2m(m+1)\mathrm{}(x_ix_j),m_{>0}.$$

It is invariant under any permutation of the coordinates $`x^1,\mathrm{},x^n`$, and also under translations $`xx+l`$ for $`l^n+\tau ^n`$ ($`^n`$ is the standard integer lattice in $`n`$ dimensions). As a result, $`u`$ will be symmetric with respect to any of the hyperplanes $`x_ix_j=m+n\tau `$. Thus, its Laurent expansion in the normal direction will have no odd terms at all, hence $`u`$ is quasi-invariant. The corresponding Schrödinger operator (2.1) is the Hamiltonian of the quantum elliptic Calogero–Moser problem. More generally, the elliptic Calogero–Moser problems related to other root systems  also lead, in the same way, to quasi-invariant potentials. In particular, for the rank-one system $`A_1`$ we have the classical Lamé operator (1.1)

###### Example 2.2.

This example is in dimension one. The potential $`u`$ has $`N`$ poles $`x_1,\mathrm{},x_N`$ and looks as

(2.4)
$$u=\underset{i=1}{\overset{N}{}}2\mathrm{}(xx_i).$$

To ensure its quasi-invariance, one has to impose the condition that $`u`$ has zero derivative at each of its poles, more explicitly:

(2.5)
$$\underset{ji}{}\mathrm{}^{}(x_ix_j)=0\text{for all}i=1,\mathrm{},N.$$

This system of equations describes the so-called ’elliptic locus’ from , which has an intimate connection with the classical elliptic Calogero–Moser system and the KdV hierarchy.

Let us call a Schrödinger operator $`L`$ with quasi-invariant elliptic potential $`u(x)`$ a generalized Lamé operator. In trigonometric and rational versions ($`\mathrm{}(x)=\mathrm{sin}^2x`$ or $`x^2`$) such operators were considered in , where their eigenfunctions were effectively constructed. From the results of  the so-called algebraic integrability of $`L`$ follows (see the paper  for the rational case). We recall the definition of the algebraic integrability in the next section, following ; let us just remark that in dimension one this coincides with the class of algebro-geometric operators which appear in the finite-gap theory, see . This motivates the following

Conjecture.The generalized Lamé operators are all algebraically integrable.

As a particular case, this contains a conjecture of  about the algebraic integrability of the elliptic Calogero–Moser problems. As we already mentioned in the introduction, for the $`A_n`$-case (2.3) this has been proved by Braverman, Etingof and Gaitsgory in . It is also known to be true in dimension one, due to Gesztesy and Weikard . In the next section we prove this conjecture under additional assumption of the usual, Liouville integrability of $`L`$ (in a slightly stronger sense). As a corollary, we will obtain the algebraic integrability of the quantum Calogero–Moser problems for integer coupling parameters.

## 3. Monodromy and algebraic integrability

Let $`L`$ be a generalized Lamé operator as defined previously. Recall that $`L`$ is completely integrable if it is a member of a commutative family of differential operators $`L_1=L,L_2,\mathrm{},L_n`$ which are algebraically independent. We assume that the $`L_i`$’s have meromorphic coefficients and are periodic with respect to the same lattice, which makes them (singular) differential operators on the torus $`T=^n/+\tau `$. The following proposition shows that possible singularities of $`L_i`$ are contained in the singular locus $`Sing`$ of the Schrödinger operator $`L`$.

###### Proposition 3.1.

Let $`L`$ be a Schrödinger operator regular in an open set $`U`$. Then any differential operator $`M`$ on $`U`$ with meromorphic coefficients commuting with $`L`$ is regular in $`U`$.

To prove the proposition, we will need the following lemma.

###### Lemma 3.2.

Let $`S(x,p)`$ be a meromorphic function on $`T^{}U=U\times V^{}`$ which is polynomial in the momentum $`p`$, and $`\{p^2,S(x,p)\}`$ is a regular function. Then $`S`$ is a regular function.

###### Proof.

Assume the contrary. The function $`S`$ is a finite sum $`S_𝐤(x)p^𝐤`$, where $`p^𝐤`$ are monomials. Let $`DU`$ be the divisor of poles of $`S`$. Take a generic point $`z_0`$ of this divisor. Near this point $`D`$ is given by an equation $`f=0`$, where $`f`$ is analytic at $`z_0`$ and $`df(z_0)0`$. Let $`S=\frac{Q}{f^k}(1+O(f))`$ near $`z_0`$ ($`Q`$ is regular at $`z_0`$, with $`Q(z_0)0`$). Then

$$\{p^2,S\}=2kp_i\frac{f}{x_i}f^{k1}Q+O(f^k).$$

But $`p_i\frac{f}{x_i}(z_0,p)0`$ for generic $`p`$. Thus, $`\{p^2,S\}`$ is singular, which is a contradiction. ∎

Now we prove the proposition. Suppose $`[M,L]=0`$. Assume $`M`$ is not regular. Then we can write $`M`$ as $`M^{}+M^{\prime \prime }`$, where $`M^{\prime \prime }`$ is regular, and $`M^{}`$ has a singular highest symbol. Then $`[L,M^{}]=[L,M^{\prime \prime }]`$ is regular. Let $`S(x,p)`$ be the symbol of $`M^{}`$. We have $`\{p^2,S\}`$ is regular (since if it is nonzero, it is the symbol of $`[L,M^{}]`$). Then the lemma implies that $`S`$ is regular. This contradiction proves the proposition.

Notice also that if $`S(x,p)`$ is the highest symbol of $`L_i`$, then from $`[L,L_i]=0`$ it follows that $`\{p^2,S\}=0`$. Thus, by Lemma 3.2, $`S`$ must be regular everywhere. Hence, each $`L_i`$ must have constant highest symbol, i.e. $`L_i=s_i()+\mathrm{}`$ for some polynomial $`s_i`$.

###### Definition.

Let us say that a Schrödinger operator $`L=\mathrm{\Delta }+u`$ in $`V=^n`$ is strongly integrable if the commuting operators $`L_1=L,\mathrm{},L_n`$ have algebraically independent homogeneous constant highest symbols $`s_1,\mathrm{},s_n`$ and if $`[V]`$ is finitely generated as a module over the ring generated by $`s_1,\mathrm{},s_n`$ or, equivalently, if the system $`s_1(\xi )=0,\mathrm{},s_n(\xi )=0`$ has the unique solution $`\xi =0`$.

Now let $`L`$ be a strongly integrable generalised Lamé operator, so we have the operators $`L_1=L,\mathrm{},L_n`$ with meromorphic coefficients on the torus $`T=^n/+\tau `$, and $`L_i=s_i()+\mathrm{}`$. First of all, $`[V]`$ is locally free as a module over $`[s_1,\mathrm{},s_n]`$. This follows from the fact, due to Serre , that if $`f:XY`$ is a finite map of smooth affine varieties of the same dimension, then $`𝒪(X)`$ is a locally free $`𝒪(Y)`$-module. Further, since $`s_i`$ are homogeneous, this module is graded, hence, it must be free. Denote by $`N`$ the rank of this free module.

Consider now the eigenvalue problem

(3.1)
$$L_i\psi =\lambda _i\psi ,i=1,\mathrm{},n,$$

with $`\lambda =(\lambda _1,\mathrm{},\lambda _n)^n`$. Then the space of solutions of this system in any simply connected domain in $`TSing`$ is N-dimensional. So we have a holonomic system of rank $`N`$ on $`T`$ with singularities along a finite union of hypertori $`SingT`$.

###### Theorem 3.3.

The holonomic system (3.1) has regular singularities.

###### Proof.

Regularity of singularities is a codimension one condition. Thus, it is sufficient to restrict our attention to a neighborhood $`U`$ of a point which lies on exactly one subtorus from the pole divisor. Let us assume that this point is $`0`$, and the subtorus is locally defined by the equation $`x_1=0`$.

###### Lemma 3.4.

For any $`r0`$, the singular part of the differential operator $`(x_1)^rL_i`$ has the order at most $`d_ir1`$, where $`d_i`$ is the degree of $`s_i`$.

###### Proof.

Let us introduce new variables $`y_j=tx_j`$, and write our operators with respect to them. Then $`L=t^2L^{(0)}+O(1)`$, $`t0`$, where $`L^{(0)}=\mathrm{\Delta }+\frac{m_\alpha (m_\alpha +1)}{x_1^2}`$ with $`m_\alpha _+`$ corresponding to the chosen hyperplane $`x_1=0`$. Let $`L_i^{(0)}`$ be the coefficient at the lowest power of $`t`$ in the expression for $`L_i`$. Then $`L_i^{(0)}`$ is homogeneous, of degree $`d_i`$ in $`t`$ (since the symbol of $`L_i`$ has this degree). On the other hand, the order of this operator is at most $`d_i`$. Thus, if this degree of $`L_i^{(0)}`$ is less than $`d_i`$, then its symbol would have to be singular. But $`L_i^{(0)}`$ commutes with $`L^{(0)}`$, so this contradicts Lemma 3.2. Thus, the degree is exactly $`d_i`$, which implies the statement. ∎

###### Lemma 3.5.

Consider a system of differential equations $`df/dz=A(z)f(z)`$ with holomorphic coefficients on a punctured disk $`0<|z|<ϵ`$ ($`f`$ is a vector function, $`A`$ is a matrix function). Assume that $`A(z)`$ is meromorphic at $`z=0`$ and that there exists an integer-valued function $`D(i)`$ such that the order of $`a_{ij}(z)`$ at $`z=0`$ is at least $`D(j)D(i)1`$. Then the system has a regular singularity at $`z=0`$.

Indeed, let us make a change of variable $`g_i=z^{D(i)}f_i`$. This will change the matrix $`A`$ into a new matrix $`\stackrel{~}{A}`$ which obviously has at most simple pole at $`z=0`$.

Now let us prove the theorem. First, let us reduce (3.1) to a first order holonomic system in a standard way. Choose a collection of $`N`$ homogeneous polynomials $`q_1,\mathrm{},q_N`$ in $`_i`$, which form a basis in $`[_1,\mathrm{},_n]`$ as a free module over $`[s_1(),\mathrm{},s_n()]`$. Let $`\psi `$ be a solution the eigenvalue problem (3.1), and consider the functions

$$q_1()\psi ,\mathrm{},q_N()\psi .$$

Then, for any polynomial $`q`$ of $`_i`$, the function $`q\psi `$ can be expressed via $`q_i\psi `$, with the coefficients depending on $`x`$ and $`\lambda `$. Indeed, assume that $`q`$ is homogeneous of degree $`d`$. We know that $`q`$ can be uniquely represented in the form $`g_jq_j`$, where $`g_j[s_1,\mathrm{},s_n]`$. Thus, $`q\psi q_jg_j(L_1,..,L_n)\psi `$ is expressed through differential polynomials of $`\psi `$ of degree smaller than $`d`$. But $`g_j(L_1,\mathrm{},L_n)\psi =g_j(\lambda _1,\mathrm{},\lambda _n)\psi `$, so, eventually, by induction we get a desired representation of $`q\psi `$.

Now observe that by Lemma 3.4, the coefficient of $`q_j\psi `$ in the expansion of $`q\psi `$ has a pole in $`x_1`$ of order at most $`\mathrm{deg}(q)\mathrm{deg}(q_j)`$. Thus, we have a holonomic system of matrix partial differential equations

$$_kq_i\psi =a_{ij}^{(k)}(x)q_j\psi k=1,\mathrm{},n,$$

which is equivalent to the system (3.1). Each of the matrices $`a_{ij}^{(k)}`$ ($`k=1,\mathrm{},n`$) satisfies the conditions of Lemma 3.5 with respect to $`z=x_1`$, for $`D(i)=\mathrm{deg}(q_i)`$. In particular, using the lemma, we conclude from the first equation $`_1q_i\psi =a_{ij}^{(1)}\psi `$ that all the solutions have at most power growth in $`x_1`$ when approaching $`x_1=0`$, hence the system (3.1) has a regular singularity at $`x_1=0`$. ∎

Now take any point $`x_0TSing`$ and consider the monodromy of the system (3.1) at point $`x_0`$, this gives an $`N`$-dimensional representation of the fundamental group $`\pi _1(TSing)`$.

###### Proposition 3.6.

The monodromy group of the system (3.1) is commutative for any $`\lambda `$.

###### Proof.

Take a hyperplane $`\pi =\pi _\alpha ^{m,n}Sing`$ and let $`P`$ be a generic point of $`\pi `$. Changing the coordinates if necessary, we may assume that $`P=(0,\mathrm{},0)`$ is the origin and $`\pi `$ is given by equation $`x_1=0`$. From the regularity of singularities it follows that there exist $`\gamma _1,\mathrm{},\gamma _r`$ such that any solution $`\psi `$ of the system (3.1) near $`P\pi `$ has a convergent series expansion in the subspace

$$\underset{j=1}{\overset{r}{}}x_1^{\gamma _j}[[x_1,\mathrm{},x_n]][\mathrm{log}(x_1)].$$

Substituting such a series into the first equation $`L\psi =\lambda _1\psi `$, one arrives at certain recurrence relations from which follows that (1) we have two ’leading exponents’ $`\gamma _1=m_\alpha `$, $`\gamma _2=m_\alpha +1`$, both integer; (2) there will be no $`\mathrm{log}(x_1)`$ involved. The latter fact is due to the quasi-invariance of $`u`$, see  for the one-dimensional case and , section 2 for a discussion in the multivariable setting.

As a result, we see that all the solutions are single-valued near $`\pi `$, so the local monodromy corresponding to a loop around $`\pi `$ is trivial. Consequently, the global monodromy group will be a homomorphic image of the commutative group $`\pi _1(T)`$. ∎

###### Corollary 3.7.

The differential Galois group of the system (3.1) is commutative for any $`\lambda `$.

###### Proof.

It is known (see ) that for a regular holonomic system on a smooth projective variety the differential Galois group coincides with the Zariski closure of the monodromy group. In our situation the proof is simple: first, we know that all solutions of the system (3.1) are meromorphic in $`^n`$. Indeed, this is true outside the set of codimension $`2`$, by the proposition above, hence it holds everywhere by Hartogs’ theorem. Now let $`F`$ be the field of ’elliptic functions’, i.e. meromorphic functions on the torus $`T`$, and $`L`$ is the solution field of (3.1) (on some simply connected domain). The monodromy gives the homomorphism

$$\pi _1(T)Dgal(L:F)$$

of the fundamental group of the torus $`T`$ into the differential Galois group. Let $`G`$ be the Zariski closure of the image of this homomorphism inside $`Dgal(L:F)`$. By the main theorem of the differential Galois theory , to prove that $`G=Dgal(L:F)`$ it suffices to show that in the solution field any $`G`$-invariant function is actulaly $`Dgal`$-inavriant, i.e. belongs to $`F`$. But any meromorphic $`G`$-invariant function is $`\pi _1(T)`$-invariant, hence elliptic. Thus, it is $`Dgal`$-invariant by definition.

This proves that $`Dgal(L:F)=G`$, and the latter is obviously commutative. ∎

Now let us remark (see ) that a quantum completely integrable system (QCIS) on a smooth algebraic variety $`X`$ of dimension $`n`$ naturally defines an embedding

$$\theta :𝒪(\mathrm{\Lambda })D(X),$$

where $`\mathrm{\Lambda }𝔸^n`$ is an affine space and $`D(X)`$ denotes the ring of differential operators on $`X`$. More generally, $`\mathrm{\Lambda }`$ can be any affine variety with $`dim(\mathrm{\Lambda })=n`$. Then we have an analogous eigenvalue problem $`\theta (g)\psi =g(\lambda )\psi `$, $`g𝒪(\mathrm{\Lambda })`$. The dimension of the local solution space of this system at generic point of $`X`$ is called the rank of a QCIS. Recall further, that a QCIS $`S=(\mathrm{\Lambda },\theta )`$ is algebraically integrable if it is dominated by another QCIS $`S^{}=(\mathrm{\Lambda }^{},\theta ^{})`$ of rank one ($`S`$ is dominated by $`S^{}`$ if there is a map of algebras $`h:𝒪(\mathrm{\Lambda })𝒪(\mathrm{\Lambda }^{})`$ such that $`\theta =\theta ^{}h`$). In our situation, this implies that apart from the operators $`L_1,\mathrm{},L_n`$, we have additional commuting operators which are not algebraic combinations of $`L_i`$ (though, of course, they are algebraically dependent with $`L_i`$). In dimension $`n=1`$ this is equivalent to saying that $`L`$ is a member of a maximal commutative ring in $`D(X)`$ of rank one; this is known to coincide with the class of algebro-geometric operators.

###### Theorem 3.8.

Any generalised Lamé operator $`L`$ which is strongly integrable is algebraically integrable.

This follows immediately from the result above and the criterion from . Moreover, according to , for generic $`\lambda `$ the solution space is generated by the quasiperiodic solutions:

###### Corollary 3.9.

There exist meromorphic $`1`$-forms $`\omega _j`$ on the torus $`T`$, with first order poles and depending analytically on $`\lambda `$, such that the functions $`\psi _j=e^{{\scriptscriptstyle \omega _j}}`$ give a basis of the solution space of (3.1) for generic $`\lambda `$.

Each of these functions will be double-Bloch, in terminology of :

(3.2)
$$\psi _j(x+l)=e^{2\pi ia_j,l}\psi _j(x),\psi _j(x+\tau l)=e^{2\pi ib_j,l}\psi _j(x),$$

for appropriate $`a_j,b_jV^{}`$ and for all $`l`$. Namely,

$$a_j,l=\frac{1}{2\pi i}_z^{z+l}\omega _j\text{and}b_j,l=\frac{1}{2\pi i}_z^{z+\tau l}\omega _i$$

(since $`\psi _j`$ has no branching along $`Sing`$, these are well defined modulo $``$). In Section 5 below we explain how one can calculate these double-Bloch solutions for the generalized Lamé operators.

###### Remark 3.10.

Our argument applies to a more general situation when one has commuting operators $`L_1,\mathrm{},L_n`$ on an abelian variety, such that the system $`L_i\psi =\lambda _i\psi `$, $`i=1,\mathrm{},n`$, is regular holonomic. Then the triviality of its local monodromy around singularities implies that this system is algebraically integrable.

## 4. Examples

### 4.1. One-dimensional case

Let $`L`$ be a one-dimensional Schrödinger operator $`L=\frac{d^2}{dx^2}+u(x)`$, $`x`$. Consider the eigenvalue problem $`L\psi =\lambda \psi `$. Now, assuming that $`u`$ is meromorphic and has a pole at $`x=0`$, let us consider the local monodromy of the solutions of this second order differential equation. Suppose that this monodromy is trivial for all $`\lambda `$, in other words, all the solutions are meromorphic at $`x=0`$. Then, using the classical Frobenius analysis, one can see that $`u`$ must have a pole of the second order, with no residue: $`u=m(m+1)x^2+O(1)`$, with integer $`m`$. Furthermore, in the series for $`u`$ at $`x=0`$ there must be no terms of order $`2j1`$ for all $`j=1,\mathrm{},m`$ (see ). This is exactly the quasi-invariance of $`u`$ with respect to the symmetry $`xx`$.

Now let $`u`$ be elliptic, with periods $`1,\tau `$. Let us demand all the solutions of the equation $`L\psi =\lambda \psi `$ to be meromorphic in $``$ (such $`u`$ are called Picard potentials in ). By the discussion above, this is equivalent to the quasi-invariance of $`u`$ at each pole. More explicitly,

$$u(x)=\underset{i=1}{\overset{M}{}}m_i(m_i+1)\mathrm{}(xx_i),\text{with}m_i_+,$$

and the poles $`x_i`$ must satisfy the following system of equations, which generalizes (2.5):

(4.1)
$$\underset{ji}{}m_j(m_j+1)\mathrm{}^{(2s1)}(x_ix_j)=0$$

for all $`i=1,\mathrm{},M`$ and $`s=1,\mathrm{},m_i`$.

Such $`L`$ automatically defines a $`QCIS`$ (of rank 2), and the regularity of singularities is obvious. Applying the results of the previous section, we conclude that $`L`$ is algebraically integrable. Algebraic integrability of $`L`$ implies the existence of a differential operator $`P`$, commuting with $`L`$. Since $`P`$ is not a polynomial of $`L`$, we may assume that it is of odd order. Thus, $`L`$ must be algebro-geometric, according to the Burchnall–Chaundy–Krichever theory, see e.g.. Thus, our result in this case is equivalent to the main result of . Notice that our approach easily extends to the case of operators of any order (cf. the remark at the end of ).

### 4.2. Quantum Calogero–Moser system

Let $`V`$ be a complex Euclidean space, $`dimV=n`$. Let $`R=\{\alpha \}`$ be a reduced irreducible root system in $`V^{}`$, $`W`$ be the corresponding Weyl group, and the parameters $`m_\alpha `$ be chosen in a $`W`$-invariant way. The corresponding Calogero–Moser operator  looks as

(4.2)
$$L=\mathrm{\Delta }+\underset{\alpha R_+}{}m_\alpha (m_\alpha +1)(\alpha ,\alpha )\mathrm{}(\alpha ,x).$$

Let $`[V]^W`$ be the ring of $`W`$-invariant polynomials. By the Chevalley theorem, it is freely generated by $`n`$ elements $`p_1,\mathrm{},p_n`$, and $`[V]`$ is a free module over $`[V]^W`$, of rank $`|W|`$. The following result has been proved in  for any $`W`$-invariant $`m_\alpha `$.

###### Theorem 4.1 (Cherednik).

For each homogenious $`p[V]^W`$ there exists a differential operator $`L_p`$ with the highest symbol $`p`$, commuting with $`L`$: $`[L,L_p]=0`$. The family $`\{L_p\}`$ is commutative.

For integer $`m_\alpha `$ the potential $`u`$ in (4.2) is quasi-invariant (cf. Example 2.1). Altogether this proves the conjecture of :

###### Corollary 4.2.

The Calogero–Moser operator (4.2) is algebraically integrable for integer $`m_\alpha `$.

Let us mention that for the classical root systems $`R=A\mathrm{}D`$ the complete set of commuting operators $`L_1,\mathrm{},L_n`$ for (4.2) was explicitly found by Ochima and Sekiguchi . Their results cover the $`BC_n`$ case, too. In that case we have $`5`$ parameters $`m,g_0,g_1,g_2,g_3`$ and the Inozemtsev operator

(4.3)
$$L=\mathrm{\Delta }+2m(m+1)\underset{i<j}{}(\mathrm{}(x_ix_j)+\mathrm{}(x_i+x_j))+\underset{i=1}{\overset{n}{}}\underset{s=0}{\overset{3}{}}g_s(g_s+1)\mathrm{}(x_i+\omega _s),$$

with $`\omega _s(s=0\mathrm{}3)`$ denoting the half-periods $`0,1/2,\tau /2,(1+\tau )/2`$. The Weyl group $`W`$ for this case is generated by the permutations of $`x_i`$ and sign flips, and Theorem 4.1 still holds true, due to . Applying Theorem 3.8 we obtain

###### Corollary 4.3.

The Inozemtsev operator (4.3) is algebraically integrable for any integer $`m,g_0,g_1,g_2,g_3`$.

### 4.3. Deformed root systems

Other known examples of the generalized Lamé operators in dimension $`>1`$ are related to deformed root systems, which appeared in . Below we describe the set of linear functionals $`𝒜=\{\alpha \}`$ and the corresponding multiplicities $`m_\alpha `$.

(1) $`A_{n,1}(m)`$ system .

It consists of the following covectors in $`^{n+1}`$:

$$\{\begin{array}{ccc}x_ix_j,\hfill & 1i<jn,\hfill & \text{with multiplicity }m,\hfill \\ x_i\sqrt{m}x_{n+1},\hfill & i=1,\mathrm{},n\hfill & \text{with multiplicity }1.\hfill \end{array}$$

Here $`m`$ is an integer parameter, and $`m`$ stands for $`m=\mathrm{max}\{m,1m\}`$.

(2) $`C_{n,1}(m,l)`$ system .

It consists of the following covectors in $`^{n+1}`$:

$$\{\begin{array}{ccc}x_i\pm x_j,\hfill & 1i<jn,\hfill & \text{with multiplicity }k,\hfill \\ 2x_i,\hfill & i=1,\mathrm{},n\hfill & \text{with multiplicity }m,\hfill \\ x_i\pm \sqrt{k}x_{n+1},\hfill & i=1,\mathrm{},n\hfill & \text{with multiplicity }1,\hfill \\ 2\sqrt{k}x_{n+1}\hfill & & \text{with multiplicity }l.\hfill \end{array}$$

Here $`k,l,m`$ are integer parameters related as $`k=\frac{2m+1}{2l+1}`$, and $`k,l,m`$ have the same meaning as in $`A_{n,1}(m)`$ case. In two-dimensional case $`n=1`$ the first group of roots is absent and there is no restriction for $`k`$ to be integer.

(3) Here is a $`BC_n`$-type generalization of the previous example. Let $`\omega _s(s=0\mathrm{}3)`$ denote the half-periods, as in $`BC_n`$ case (4.3). The set of linear functionals $`\alpha 𝒜`$ and the corresponding multiplicities look as follows:

$$\{\begin{array}{ccc}x_i\pm x_j,\hfill & 1i<jn,\hfill & \text{with multiplicity }k,\hfill \\ x_i+\omega _s,\hfill & i=1,\mathrm{},n\hfill & \text{with multiplicity }m_s,(s=0\mathrm{}3),\hfill \\ x_i\pm \sqrt{k}x_{n+1},\hfill & i=1,\mathrm{},n\hfill & \text{with multiplicity }1,\hfill \\ \sqrt{k}x_{n+1}+\omega _s\hfill & & \text{with multiplicity }l_s,(s=0\mathrm{}3).\hfill \end{array}$$

Here $`k,l_s,m_s`$ are nine integer parameters related through $`k=\frac{2m_s+1}{2l_s+1}`$ for all $`s=0\mathrm{}3`$. The previous case corresponds to $`m_sm`$ and $`l_sl`$. Again, in case $`n=1`$ the first group of roots is absent and $`k`$ may not be integer.

(4) Hietarinta operator .

In this case we have three covectors in $`^3`$,

$$\alpha =a_1x_1a_2x_2,\beta =a_2x_2a_3x_3,\gamma =a_1x_1a_3x_3,$$

with $`m_\alpha =m_\beta =m_\gamma =1`$. Here $`a_i`$ are arbitrary complex parameters such that $`a_1^2+a_2^2+a_3^2=0`$. Notice that the system is essentially two-dimensional since $`\alpha +\beta +\gamma =0`$.

(5) $`A_{n1,2}(m)`$ system .

It consists of the following covectors in $`^{n+2}`$:

$$\{\begin{array}{ccc}x_ix_j,\hfill & 1i<jn,\hfill & \text{with multiplicity }m,\hfill \\ x_i\sqrt{m}x_{n+1},\hfill & i=1,\mathrm{},n\hfill & \text{with multiplicity }1,\hfill \\ x_i\sqrt{1m}x_{n+2},\hfill & i=1,\mathrm{},n\hfill & \text{with multiplicity }1,\hfill \\ \sqrt{m}x_{n+1}\sqrt{1m}x_{n+2}\hfill & & \text{with multiplicity }1.\hfill \end{array}$$

Notice that for $`m=1`$ this system coincides with the system $`A_{n,1}(2)`$ above.

In all these cases a direct check shows that the corresponding potential $`u`$ is quasi-invariant. Notice that in all cases $`u`$ is symmetric with respect to $`s_\alpha `$ as soon as $`m_\alpha >1`$. Thus, one have to check the quasi-invariance only for those $`\alpha `$ where $`m_\alpha =1`$.

We believe that all these operators are algebraically integrable. In the next section we check this for all two-dimensional examples.

### 4.4. Two-dimensional case

Let us check that all known two-dimensional generalized Lamé operators are algebraically integrable. Apart from the root systems $`A_2,BC_2`$ and $`G_2`$ considered previously, we have three deformed cases, namely, the $`A_{1,1}(m)`$ case, the Hietarinta operator and the deformed $`BC_2`$ case.

First, let us consider the $`A_{1,1}(m)`$ case:

(4.4)
$$\begin{array}{c}L=_1^2_2^2_3^2+2m(m+1)\mathrm{}(x_1x_2)+\hfill \\ \hfill 2(m+1)\mathrm{}(x_1\sqrt{m}x_3)+2(m+1)\mathrm{}(x_2\sqrt{m}x_3).\end{array}$$

In this case we can use the result of  where the complete integrability of (4.4) was established. First, $`L`$ obviously commutes with $`L_0=_1+_2+\frac{1}{\sqrt{m}}_3`$.

###### Proposition 4.4 ().

For any $`m`$ there exists a third order operator $`L_2`$ commuting with $`L_0,L`$, and its highest symbol is $`_1^3+_2^3+\sqrt{m}_3^3`$.

It is easy to check that as soon as $`m1`$, the highest symbols of $`L_0,L,L_2`$ will satisfy the requirements of the strong integrabilty (notice that for $`m=1`$ the operator $`L`$ is trivial). As a result, for integer $`m`$ we obtain the algebraic integrability of the operator (4.4). Note that for the special case $`m=2`$ the algebraic integrability of (4.4) was demonstrated in  by presenting an explicit extra operator commuting with $`L_0,L,L_2`$.

Now let us consider the Hietarinta operator

(4.5)
$$\begin{array}{c}L=_1^2_2^2_3^2+2(a_1^2+a_2^2)\mathrm{}(a_1x_1a_2x_2)+\hfill \\ \hfill 2(a_2^2+a_3^3)\mathrm{}(a_2x_2a_3x_3)+2(a_3^2+a_1^2)\mathrm{}(a_3x_3a_1x_1),\end{array}$$

where $`a_1^2+a_2^2+a_3^2=0`$. Such $`L`$ commutes with $`L_0=a_1^1_1+a_2^1_2+a_3^1_3`$, and we need one more operator for the complete integrability (we assume that all $`a_i`$ are nonzero, otherwise $`L`$ is reducible). Such operator was found in .

###### Proposition 4.5 ().

There exists a third order operator $`L_2`$ commuting with $`L_0,L`$ above, and its highest symbol is $`a_1_1^3+a_2_2^3+a_3_3^3`$.

If we denote the highest symbols of these three operators as $`s_1,s_2,s_3`$, then it is easy to check that the system $`s_1(\xi )=s_2(\xi )=s_3(\xi )=0`$ has the only solution $`\xi =0`$; the only exception is the case when $`a_1^3=a_2^3=a_3^3`$. As a result, we conclude that for all values of the parameters (apart from the case $`a_1^3=a_2^3=a_3^3`$) the Hietarinta operator is strongly integrable, thus, it is algebraically integrable. Note that this also follows from  where one more operator commuting with $`L_0,L,L_2`$ was found.

###### Remark 4.6.

In the case $`a_1^3=a_2^3=a_3^3`$ the operator (4.5) is still algebraically integrable, although it is no longer strongly integrable.

Finally, let us consider the deformed $`BC_2`$ case. The Schrödinger operator $`L`$ has the following form:

(4.6)
$$L=_x^2_y^2+U(x,y),$$

where $`U=2(k+1)(\mathrm{}(x+\sqrt{k}y)+\mathrm{}(x\sqrt{k}y))+v(x)+w(y)`$ and $`v,w`$ are given by the expressions

(4.7)
$$v=\underset{s}{}m_s(m_s+1)\mathrm{}(x+\omega _s),w=k\underset{s}{}l_s(l_s+1)\mathrm{}(\sqrt{k}y+\omega _s).$$

Here $`\omega _0,\omega _1,\omega _2,\omega _3`$ are the half-periods and $`m_s,l_s`$ and $`k`$ are nine parameters such that $`k=(2m_s+1)/(2l_s+1)`$ for all $`s=0,1,2,3`$ (thus, effectively, $`L`$ contains five independent parameters).

###### Proposition 4.7.

For any values of the parameters $`m_s,l_s,k`$ such that $`k=\frac{2m_s+1}{2l_s+1}`$ the following operator commutes with $`L`$:

$`M=`$ $`_x^4k_y^4+2(Uw)_x^24\sqrt{k}(k+1)(\mathrm{}(x+\sqrt{k}y)\mathrm{}(x\sqrt{k}y))_x_y+`$
$`2k(Uv)_y^2+(2v^{}+2(k+1)(2k)(\mathrm{}^{}(x+\sqrt{k}y)+\mathrm{}^{}(x\sqrt{k}y)))_x+`$
$`(2k^2w^{}+2\sqrt{k}(k+1)(2k1)(\mathrm{}^{}(x\sqrt{k}y)\mathrm{}^{}(x+\sqrt{k}y)))_y`$
$`(k+1)^3\mathrm{}(x+\sqrt{k}y)\mathrm{}(x\sqrt{k}y)+8(k^3+1)(\mathrm{}^2(x+\sqrt{k}y)+\mathrm{}^2(x\sqrt{k}y))`$
$`4(k+1)(v+kw)(\mathrm{}(x+\sqrt{k}y)+\mathrm{}(x\sqrt{k}y))+v^{\prime \prime }v^2+k(w^{\prime \prime }w^2).`$
Here $`v^{},w^{}`$ and so on are the derivatives with respect to the corresponding variable.

Now since the system $`\xi _1^2+\xi _2^2=\xi _1^4+k\xi _2^4=0`$ does not have nontrivial solutions as soon as $`k1`$, we conclude that $`L`$ is strongly integrable (for $`k=1`$ this is also true because $`L`$ is reducible in that case). Hence, for any integer values of the parameters $`l_s,m_s`$ the operator (4.6)–(4.7) is algebraically integrable.

## 5. Bloch solutions

Let $`L`$ be a generalized Lamé operator which is strongly integrable, thus algebraically integrable. We know already that for generic $`\lambda `$ the solution space of (3.1) is spanned by the meromorphic double-Bloch solutions. Now we are going to explain how one can, in principle, calculate them.

Let $`𝒲`$ denote the following linear subspace in the space of meromorphic functions on $`^n`$. First, its elements are holomorphic everywhere apart from the singular locus $`Sing=\pi _\alpha ^{m,n}`$ of the operator $`L`$, where they may have poles, of order $`m_\alpha `$ along $`\pi _\alpha ^{m,n}`$. Next, take any hyperplane $`\pi =\pi _\alpha ^{m,n}`$ with $`s_\pi `$ denoting the orthogonal reflection with respect to $`\pi `$. Then any function $`\phi 𝒲`$ must have the following property:

(5.1)
$$\phi (x)(1)^{m_\alpha }\phi (s_\pi x)\text{ is divisible by}(\alpha (x)mn\tau )^{m_\alpha +1}.$$

###### Proposition 5.1.

The subspace $`𝒲`$ is stable under the action of $`L`$: $`L(𝒲)𝒲`$. Furthermore, any meromorphic eigenfunction $`\phi `$ of $`L`$ must belong to $`𝒲`$. The same is true for any of the commuting operators $`L_2,\mathrm{},L_n`$: $`L_i(𝒲)𝒲`$.

###### Proof.

Consider any hyperplane $`\pi Sing`$ and adjust the coordinates in such a way that $`\pi `$ is given by equation $`x_n=0`$. Take a generic point in $`\pi `$ and expand $`\phi `$ in Laurent series in normal direction to $`\pi `$, i.e. $`\phi =_ja_j(x_n)^j`$, $`a_j=a_j(x_1,\mathrm{},x_{n1})`$. Now put

$$M:=\{m_\alpha +2_0\}\{m_\alpha +1+2_0\}$$

and split the series into two parts, with $`jM`$ and with $`j\backslash M`$: $`\phi =\phi _1+\phi _2`$. First claim is that an application of $`L`$ to $`\phi _1`$ will produce a series of a similar kind. This follows directly from the quasi-invariance of $`u`$, proving the first part of the proposition. On the other hand, if $`a_j(x_n)^j`$ is the first nonzero term in $`\phi _2`$, then applying $`L`$ to $`\psi _2`$ will give a series starting from $`(x_n)^{j2}`$, which would contradict the equation $`L\phi =\lambda \phi `$, thus proving the second claim.

In a similar way, if $`L^{}L=LL^{}`$ for some other operator $`L^{}`$, then $`L^{}(L)^r=(L)^rL^{}`$ for any $`r1`$. Thus, if $`𝒲^{}:=L^{}(𝒲)`$ then $`L^r(𝒲^{})𝒲^{}`$ for any $`r`$. Now suppose we could find a function $`\phi `$ in $`𝒲^{}`$ which is not in $`𝒲`$. Then, consider a series expansion of $`\phi `$ in the direction, normal to $`\pi `$, and split it into $`\phi =\phi _1+\phi _2`$ as above. If $`\phi _20`$, then $`L^r\phi `$ would have a pole of an arbitrarily high order along $`\pi `$ (as $`r`$ increases), which is impossible for an element in $`𝒲^{}`$ (since $`𝒲^{}=L^{}𝒲`$ and $`L^{}`$ has meromorphic coefficients). This would contradict the inclusion $`L^r(𝒲^{})𝒲^{}`$. Thus, $`𝒲^{}=𝒲`$. ∎

Now let $`\psi `$ be a double-Bloch solution of (3.1), so for appropriate $`a,bV^{}`$ and for any $`l`$ we have:

(5.2)
$$\psi (x+l)=\psi (x)e^{2\pi ia,l},\psi (x+\tau l)=\psi (x)e^{2\pi ib,l}.$$

We know that $`\psi `$ is meromorphic in $`^n`$ with possible poles along the hyperplanes $`\pi _\alpha ^{m,n}`$, of order $`m_\alpha `$. In our discussion below we restrict ourselves to the case when all the linear functions $`\alpha 𝒜`$ have zero constant term, i.e. $`\alpha =\alpha _0V^{}`$ for all $`\alpha `$, so $`\alpha (x)=\alpha ,x`$. Everything extends to the general case with obvious modifications.

As a result, we see that $`\psi `$ can be presented in the form

(5.3)
$$\psi =\mathrm{\Phi }/\delta ,\delta (x)=\underset{\alpha 𝒜}{}\theta (\alpha ,x)^{m_\alpha },$$

for some holomorphic $`\mathrm{\Phi }(x)`$. Here $`\theta =\theta _1`$ is the classical (odd) Jacobi theta function,

(5.4)
$$\theta (z)=\underset{n}{}\mathrm{exp}(\pi i(n+1/2)^2\tau +2\pi i(n+1/2)(z+1/2)).$$

Recall that $`\theta (z)`$ has the following translation properties in $`z`$:

$$\theta (z+1)=\theta (z),\theta (z+\tau )=e^{2\pi iz\pi i\tau }\theta (z).$$

This determines the translation properties of $`\delta `$. To write them down, it is convenient to introduce the following linear map $`\mathrm{\Omega }:VV^{}`$ which is defined as

(5.5)
$$\mathrm{\Omega }:x\underset{\alpha 𝒜}{}m_\alpha \alpha ,x\alpha .$$

Note that $`\mathrm{\Omega }`$ maps the lattice $``$ to a sublattice in $`=\mathrm{Hom}(,)`$. We also need a covector $`\varrho =\frac{1}{2}_{\alpha 𝒜}m_\alpha \alpha `$.

Under these notations, we have the following translation formulas for any $`l`$:

(5.6)
$$\delta (x+l)=e^{2\pi i\varrho ,l}\delta (x),$$

(5.7)
$$\delta (x+l\tau )=e^{2\pi i\varrho ,l2\pi i\mathrm{\Omega }l,x\pi i\mathrm{\Omega }l,l\tau }\delta (x).$$

As a corollary of (5.2) and (5.6)-(5.7), we conclude that the numerator $`\mathrm{\Phi }`$ in (5.3) must have the translation properties as follows:

$$\mathrm{\Phi }(x+l)=e^{2\pi ia+\varrho ,l}\mathrm{\Phi }(x),$$

$$\mathrm{\Phi }(x+l\tau )=e^{2\pi ib+\varrho ,l2\pi i\mathrm{\Omega }l,x\pi i\mathrm{\Omega }l,l\tau }\mathrm{\Phi }(x).$$

The vector space of entire functions with such properties is finite-dimensional. Indeed, let $`\omega (x,y)`$ denote the bilinear form on $`V`$ associated with the operator $`\mathrm{\Omega }`$:

(5.8)
$$\omega (x,y)=\underset{\alpha 𝒜}{}m_\alpha \alpha ,x\alpha ,y.$$

It is symmetric positive and integer-valued on the lattice $``$. Let us define $`\mathrm{\Theta }\left[\genfrac{}{}{0.0pt}{}{p}{q}\right]`$ by the following series:

(5.9)
$$\mathrm{\Theta }\left[\genfrac{}{}{0.0pt}{}{p}{q}\right](x)=\underset{l}{}\mathrm{exp}(2\pi i\omega (l+p,x+q)+\pi i\tau \omega (l+p,l+p)).$$

It has the following translation properties:

(5.10)
$$\mathrm{\Theta }\left[\genfrac{}{}{0.0pt}{}{p}{q}\right](x+l)=e^{2\pi i\omega (l,p)}\mathrm{\Theta }\left[\genfrac{}{}{0.0pt}{}{p}{q}\right](x),$$

(5.11)
$$\mathrm{\Theta }\left[\genfrac{}{}{0.0pt}{}{p}{q}\right](x+l\tau )=e^{2\pi i\omega (l,x+q)\pi i\tau \omega (l,l)}\mathrm{\Theta }\left[\genfrac{}{}{0.0pt}{}{p}{q}\right](x).$$

It is easy to show (see e.g.) that the space of holomorphic functions with such translation properties has dimension equal to $`[:\mathrm{\Omega }]`$; this is equal to $`det(\mathrm{\Omega }_{ij})`$ where $`\mathrm{\Omega }_{ij}=\mathrm{\Omega }e_i,e_j`$ for some basis $`e_1,\mathrm{},e_n`$ of $``$. A natural basis in this space is given by the functions $`\mathrm{\Theta }\left[\genfrac{}{}{0.0pt}{}{p}{q+r}\right]`$ with $`r\mathrm{\Omega }^1`$ running over the set of representatives in $`\mathrm{\Omega }^1()/`$.

For the later purposes, let us use slightly different basis, namely, the functions

(5.12)
$$\mathrm{\Phi }_r=e^{k,x}\mathrm{\Theta }(x+\gamma +r),\mathrm{\Theta }:=\mathrm{\Theta }\left[\genfrac{}{}{0.0pt}{}{0}{0}\right],r\mathrm{\Omega }^1()/.$$

It is easy to relate the parameters $`kV^{}`$ and $`\gamma V`$ to $`p,q`$:

$$k=2\pi i\mathrm{\Omega }p,\gamma =q+\tau p.$$

Let us denote the linear space generated by the functions (5.12) as $`𝒰_{k,\gamma }`$.

We conclude that $`\psi `$ must belong to linear space $`\delta ^1𝒰_{k,\gamma }`$ with $`k,\gamma `$ related in a simple way to the ’quasimomenta’ $`a,b`$ in (5.2). Now recall Proposition 5.1. It implies that $`\psi `$ must belong to the (finite-dimensional) subspace $`𝒲_{k,\gamma }=\left(\delta ^1𝒰_{k,\gamma }\right)𝒲`$. It also implies that $`L`$ (as well as any of $`L_i`$’s) preserves this finite-dimensional space, so we can eventually find $`\psi `$ by diagonalizing the action of $`L`$ on $`𝒲_{k,\gamma }`$. Note that since the double-Bloch solutions $`\psi `$ form an $`n`$-parametric family (with $`\lambda =(\lambda _1,\mathrm{},\lambda _n)`$ being the parameters), this space $`𝒲_{k,\gamma }`$ will be nonzero only for $`(k,\gamma )`$ belonging to a certain $`n`$-dimensional subvariety. In most cases $`dim𝒲_{k,\gamma }1`$, so the (unique) function $`\psi _{k,\gamma }`$ generating $`W_{k,\gamma }0`$ will be an eigenfunction for $`L`$ automatically.

All this simplifies a little for the Calogero–Moser models, so let us consider this case in more detail. For a given reduced, irreducible root system $`R`$ in $`V^{}`$ and a fixed $`W`$-invariant $`m(\alpha )`$, we consider the Calogero–Moser operator

(5.13)
$$L=\mathrm{\Delta }+\underset{\alpha R_+}{}m_\alpha (m_\alpha +1)(\alpha ,\alpha )\mathrm{}(\alpha ,x|\tau ).$$

The Bloch solutions must be of the form $`\psi =\delta ^1\mathrm{\Phi }`$ as in (5.3). Note that $`\delta `$ in this case has the following symmetry:

(5.14)
$$\delta (wx)=\epsilon _m(w)\delta (x)\text{for any}wW,$$

where $`\epsilon _m`$ is the one-dimensional character of $`W`$ such that

(5.15)
$$\epsilon _m(s_\alpha )=(1)^{m_\alpha }.$$

The bilinear symmetric form (5.8) is obviously $`W`$-invariant, and since $`R`$ is irreducible, $`\omega `$ must be proportional to the $`W`$-invariant scalar product $`(x,y)`$. So, $`\omega (x,y)=\kappa (x,y)`$ for some $`\kappa `$ which depends on $`m_\alpha `$. (For instance, if $`R`$ consists of one $`W`$-orbit only and $`m_\alpha m`$, one has $`\omega (x,y)=mh(x|y)`$ where $`h=h(R)`$ is the Coxeter number and the form $`(x|y)`$ on $`V`$ is normalized in such a way that $`(\alpha ^{}|\alpha ^{})=2`$ for all $`\alpha ^{}R^{}`$.)

The lattice $``$ in this case is the coweight lattice $`P^{}`$ of $`R`$, while $``$ is the root lattice $`Q`$. Still, the numerator $`\mathrm{\Phi }`$ must belong to the finite-dimensional space $`𝒰_{k,\gamma }`$ spanned by the functions (5.12).

We already know that a Bloch solution $`\psi =\delta ^1\mathrm{\Phi }`$ appears only for those $`k,\gamma `$ when $`\delta ^1𝒰_{k,\gamma }𝒲0`$. Such $`k,\gamma `$ can be effectively determined. Indeed, due to (5.14), the conditions (5.1) for $`\psi `$ reduce to the quasi-invariance of $`\mathrm{\Phi }`$:

(5.16)
$$\mathrm{\Phi }(x)\mathrm{\Phi }(s_\alpha x)\text{is divisible by}\alpha ,x^{2m_\alpha +1}\text{for all}\alpha R.$$

(These are local conditions near $`\pi _\alpha =\{x:\alpha ,x=0\}`$, similar conditions for other hyperplanes $`\pi _\alpha ^{m,n}Sing`$ will follow because $`\psi `$ is quasiperiodic.)

These conditions can be rewritten as

(5.17)
$$\alpha ^{},^{2j1}\mathrm{\Phi }0\text{for}\alpha ,x=0\text{and}j=1,\mathrm{},m_\alpha ,$$

with $`\alpha ^{},`$ denoting the derivative in $`\alpha ^{}`$-direction.

Now recall that we have the period lattice $`=P^{}`$, and $`\mathrm{\Phi }`$ belongs to the linear space $`𝒰_{k,\gamma }`$ of the functions with the translation properties (5.10)–(5.11). Consider the following sublattice $`^\alpha `$:

(5.18)
$$^\alpha :=^{}^{\prime \prime },^{}=\mathrm{ker}\alpha ,^{\prime \prime }=\alpha ^{}.$$

Let $`𝒰_{k,\gamma }^\alpha `$ denote the space of theta functions with the same translation properties, but for the translations $`l`$ from $`^\alpha `$ only. Obviously, we have a natural inclusion map $`𝒰_{k,\gamma }𝒰_{k,\gamma }^\alpha `$. It is possible to describe this linear map explicitly, using the standard bases in both spaces (look at the formula (6.4) below which is a particular example of such relation). According to (5.18), the lattice $`^\alpha `$ is the direct orthogonal sum of two sublattices. Thus, the corresponding theta functions from $`𝒰_{k,\gamma }^\alpha `$ will be the products of the $`(n1)`$-dimensional theta functions related to $`^{}`$ and the one-dimensional theta functions related to the lattice $`^{\prime \prime }`$. This corresponds to the decomposition of $`𝒰_{k,\gamma }^\alpha `$ into a tensor product: $`𝒰_{k,\gamma }^\alpha =(𝒰_{k,\gamma }^\alpha )^{}(𝒰_{k,\gamma }^\alpha )^{\prime \prime }`$. Applying derivative in $`\alpha ^{}`$ direction will affect the one-dimensional theta functions only. As a result, for each $`j`$ we have an explicit linear map $`\mathrm{\Gamma }^{\alpha ,j}`$ from $`𝒰_{k,\gamma }`$ to $`(𝒰_{k,\gamma }^\alpha )^{}`$ given by

$$\mathrm{\Gamma }^{\alpha ,j}\phi =\alpha ^{},^{2j1}\phi |_{\alpha ,x=0}.$$

This map is given by a matrix whose entries are certain combinations of one-dimensional theta-functions and their derivatives. Now we can organize all these maps for $`\alpha R_+`$ into one big linear map

(5.19)
$$\mathrm{\Gamma }:𝒰_{k,\gamma }\underset{\genfrac{}{}{0.0pt}{}{\alpha R_+}{j=1\mathrm{}m_\alpha }}{}(𝒰_{k,\gamma }^\alpha )^{},\mathrm{\Gamma }=(\mathrm{\Gamma }^{\alpha ,j})_{\alpha R_+,j=1,\mathrm{},m_\alpha },$$

with $`\mathrm{\Gamma }^{\alpha ,j}`$ defined above. We can think of $`\mathrm{\Gamma }`$ as an $`M\times N`$ matrix, where $`N,M`$ are the dimensions of the source and the target spaces, respectively.

The outcome is the following: a double-Bloch solution $`\psi `$ appears exactly for those $`k,\gamma `$ where this linear map has nontrivial kernel. This gives equations on such $`k,\gamma `$ (by equating to zero all $`N\times N`$ minors of the matrix $`\mathrm{\Gamma }`$). In its turn, the kernel will determine a corresponding Bloch eigenfunction. (If the kernel has dimension $`>1`$, it still defines an invariant subspace for the action of $`L`$, so we have at least one double-Bloch solution.)

So, let $`\stackrel{~}{𝒞}`$ denote an analytic subvariety in $`V^{}\times V`$ given by

(5.20)
$$\stackrel{~}{𝒞}=\{(k,\gamma )|\mathrm{ker}\mathrm{\Gamma }0\}.$$

Formulas (5.10)-(5.12) make clear that $`\stackrel{~}{𝒞}`$ is invariant under the following transformations of $`k,\gamma `$:

(5.21)
$$(k,\gamma )(k,\gamma +l),l\mathrm{\Omega }^1,$$

(5.22)
$$(k,\gamma )(k+2\pi i\mathrm{\Omega }l,\gamma +\tau l),l\mathrm{\Omega }^1.$$

Now, for a function of one variable $`f(z)=e^{kz}\theta (z+c)`$ its derivatives at $`z=0`$ are obviously polynomial in $`k`$. Therefore $`\stackrel{~}{𝒞}`$, after being factored by the translations above, can be considered as an algebraic covering of an abelian variety (a product of elliptic curves) $`^n/^{}+\tau ^{}`$ where $`^{}=\mathrm{\Omega }^1`$.

###### Proposition 5.2.

The double-Bloch eigenfunctions of the Calogero–Moser operator (5.13) are parametrized by the points of an algebraic variety which is a covering of an abelian variety $`^n/^{}+\tau ^{}`$ where $`^{}=\mathrm{\Omega }^1Q`$, $`Q`$ is the root lattice and the map $`\mathrm{\Omega }`$ is defined by (5.5).

A similar analysis applies to any (integrable) generalized Lamé operator, so the double-Bloch eigenfunctions are also parametrized by the points of an algebraic variety covering a product of elliptic curves.

Now let $`𝒞`$ denote the result of factoring the variety (5.20) by the translations (5.21)–(5.22). It is an algebraic variety parametrizing the double-Bloch eigenfunctions of $`L`$. Below, following , we will refer to it as the Hermite–Bloch variety for $`L`$. It differs from the complex Bloch variety, traditionally defined as the set of $`(\mu ,E)(^\times )^n\times `$ such that there exists $`\psi `$ with $`L\psi =E\psi `$ and $`\psi (x+l_i)=\mu _i\psi (x)`$ where $`l_1,\mathrm{}l_n`$ is a basis in $``$. Note that the latter is a transcendental complex analytic variety.

###### Remark 5.3.

In case of the Calogero–Moser operator related to a root system $`R`$, there is a natural action of the Weyl group $`W`$ on the Hermite–Bloch variety $`𝒞`$. Also, there is a natural projection of $`𝒞`$ onto $`^n`$ sending $`\psi =\psi _{k,\gamma }`$ to the set of eigenvalues $`\lambda _i`$, $`L_i\psi =\lambda _i\psi `$. This is a $`|W|`$-sheeted covering and the Weyl group acts on $`𝒞`$ by permuting the points in the fiber.

###### Remark 5.4.

Note that our results do not contradict the theorem of Feldman–Knörrer–Trubowitz  in dimension two, since their result only applies to a real-valued smooth potential $`u`$ in $`^2`$.

The Hermite–Bloch variety $`𝒞`$ is a subvariety in the total space of a certain bundle over the product of elliptic curves defined by (5.21)–(5.22). This bundle naturally compactifies to a bundle with the fibers isomorphic to the projective space $`^n`$. As a result, $`𝒞`$ compactifies to a projective variety, covering the product of elliptic curves.

The variety $`𝒞`$ is closely related to the so-called spectral variety which is defined as follows. Suppose $`L`$ is a strongly integrable generalized Lamé operator, so we have $`n`$ commuting operators $`L_1=L,\mathrm{},L_n`$, which generate a commutative subalgebra in the ring of PDO with meromorphic elliptic coefficients. Then by , theorem 2.2, the centralizer of this subalgebra will be a maximal commutative ring which we will denote by $`𝒵(L)`$ (using Proposition 5.1 one can show that the operators in this ring will share a common family of the double-Bloch eigenfunctions). Each operator in $`𝒵(L)`$ must have constant highest symbols, by Lemma 3.2. Then from the strong integrability we immediately derive that $`𝒵(L)`$ is finitely generated. Thus, $`\mathrm{Spec}𝒵(L)`$ defines an affine algebraic variety, which we call the spectral variety. It is not quite clear whether the spectral variety is isomorphic to the Hermite–Bloch variety (for instance, the latter may not be affine), but at least they must be birationally equivalent.

Finally, let us remark on some algebraic geometry behind the double-Bloch solutions and the Hermite–Bloch variety for the Calogero–Moser system (5.13). We consider the torus $`T=^n/Q^{}+\tau Q^{}`$ where $`Q^{}=R^{}`$ is the coroot lattice. Let us define the following subsheaf $`𝒬𝒪(T)`$ of the structure sheaf of $`T`$ by requiring its local sections to have zero normal derivatives of order $`1,3,\mathrm{},2m_\alpha 1`$ along each of the hyperplanes $`\alpha ,x=0`$ (considered as hypertori in $`T`$). The sheaf $`𝒬`$ can be considered as the structure sheaf $`𝒪(X)`$ of a singular variety $`X`$, with $`T`$ being its injective normalization (cf. ). Such $`X`$ is projective; it is a $`|W|`$-sheeted covering of the weighted projective space $`T/W`$ considered by Looijenga , see also . Notice that from the results of  it follows that $`X`$ is Cohen-Macaulay and Gorenstein. Let us consider now the group $`\mathrm{Pic}(X)`$ of invertible sheaves on $`X`$. Then each of the double-Bloch solutions $`\psi _{k,\gamma }`$ represents a meromorphic section of a degree zero line bundle on $`X`$ (to define degree, we use the pull-back to the torus $`T`$). In this way the Hermite–Bloch variety for the Calogero–Moser system becomes an $`n`$-dimensional subvariety in $`\mathrm{Pic}^0(X)`$. It would be interesting to study this relation in more detail.

###### Remark 5.5.

An interesting thing is to analyse how the spectral variety changes when $`\tau `$ goes to $`+i\mathrm{}`$ (trigonometric limit). In this limit the spectral variety becomes rational and is relatively well understood. Thus, one could think of the whole family depending on $`\tau `$ as a deformation of this rational variety. This point of view was used in  to construct the spectral surface in the simplest $`A_2`$ case.

### 5.1. Discrete spectrum eigenstates

Let us explain how the Bloch solutions can be used to construct the discrete spectrum eigenstates of $`L`$. Our discussion is strictly confined to the Calogero–Moser operator (5.13). We take a purely imaginary $`\tau `$, this ensures that the potential in (5.13) is real-valued for $`x^n`$. The Calogero–Moser operator $`L`$ is defined on a dense subset of $`L^2(^n)`$ and it is self-adjoint only formally, and its Bloch solutions are singular. It has square-integrable eigenstates, though. Namely, let $`\psi =\psi _{k,\gamma }`$ be one of the double-Bloch solutions constructed in the previous section. Given such a $`\psi `$, let us symmetrize it as follows:

(5.23)
$$\mathrm{\Psi }(x)=\underset{wW}{}\epsilon _m(w)detw\psi (wx),$$

where $`W`$ is the Weyl group of the root system $`R`$ and $`\epsilon _m`$ is the character (5.15). The Calogero–Moser operator $`L`$ is $`W`$-invariant, thus the constructed $`\mathrm{\Psi }`$ will be again its eigenfunction (by the same reason, it will be an eigenfunction for all commuting operators $`L_i`$). A priori, $`\mathrm{\Psi }`$ might have poles in $`^n`$ along the hyperplanes $`\alpha ,x=c`$, $`c`$. However, it is easy to see that $`\mathrm{\Psi }`$ has no poles along the hyperplanes $`\alpha ,x=0`$. This follows immediately from the properties (5.1) of $`\psi `$. To avoid the appearance of singularities on other hyperplanes $`\alpha ,x=c`$, one has to impose the condition that all the terms $`\psi (wx)`$ in the sum (5.23) have the same Bloch–Floquet multipliers with respect to a shift $`xx+l`$ with $`l=P^{}`$. This means that $`\mathrm{exp}k,l=\mathrm{exp}wk,l`$ for all $`wW`$ and $`lP^{}`$, which in its turn implies that $`k`$ belongs to the lattice $`2\pi iP`$. So, we have the following result.

###### Proposition 5.6.

Let $`L`$ be the Calogero–Moser operator (5.13). Then for any point $`(k,\gamma )`$ of its Bloch–Hermite variety which satisfies an additional condition $`k2\pi iP`$ (with $`P`$ being the weight lattice for $`R`$), the corresponding function (5.23) (if nonzero) will be a nonsingular in $`^n`$ eigenfunction of the Calogero–Moser operator (5.13) and of the higher operators $`L_2,\mathrm{},L_n`$.

By construction, $`\mathrm{\Psi }`$ vanishes along the hyperplanes $`\alpha ,x`$, and it gets a factor of $`(1)^{m_\alpha +1}`$ under the orthogonal reflection with respect to such a hyperplane. Since these are the reflection hyperplanes of the affine Weyl group of $`R`$, they cut $`^n`$ into its fundamental domains (alcoves), so the restriction of $`\mathrm{\Psi }`$ to each alcove will be, essentially, the same. We can restrict $`\mathrm{\Psi }`$ to any alcove, extending it by zero outside, and this gives us a finitely supported smooth eigenfunction of $`L`$ (notice that in the complex domain it still has poles). We see from this that the discrete spectrum of $`L`$ in $`L^2(^n)`$ is infinitely degenerate (one says that the spectral problem for $`L`$ splits into identical spectral problems on each of the alcoves). Morally, this is the reason why one should expect the same spectrum considering $`L`$ not on $`L^2(^n)`$ but on the space $`L^2(T)^W`$ of $`W`$-invariant functions on the torus $`T=^n/Q^{}`$, as in . The latter case is simpler from the technical point of view, since the operator $`L`$ is essentially self-adjoint on $`L^2(T)^W`$, see  for the details.

In  Komori and Takemura considered the elliptic Calogero–Moser problems as a perturbation (in $`\tau `$) of the trigonometric case $`\tau =+i\mathrm{}`$, and Theorem 3.7 of  claims that for sufficiently small $`p=e^{2\pi i\tau }`$ the family of eigenfunctions (Jack polynomials) which corresponds to $`p=0`$, admits analytic continuation in $`p`$, and the resulting functions will give rise to a complete orthogonal family of eigenfunctions of $`L`$ in $`L^2(T)^W`$. One can show that our family in the limit $`\tau +i\mathrm{}`$ specializes to the Jack polynomials. Comparing this with the previous discussion, we conclude that our family must coincide with the one considered in .

## 6. Calogero–Moser model of $`B_2`$ type

In this section we consider the following 2-dimensional Schrödinger operator

(6.1)
$$L=\mathrm{\Delta }+2\mathrm{}(x_1)+2\mathrm{}(x_2)+4\mathrm{}(x_1x_2)+4\mathrm{}(x_1+x_2),$$

where $`\mathrm{}(z)=\mathrm{}(z|\tau )`$ is the Weierstrass $`\mathrm{}`$-function with the periods $`1,\tau `$ ($`\mathrm{Im}\tau >0`$). Our goal is to calculate its double-Bloch eigenfunctions, i.e. such $`\psi `$ that

(6.2) $`\psi (x+e_j)`$ $`=\lambda _j\psi (x)`$
(6.3) $`\psi (x+\tau e_j)`$ $`=\mu _j\psi (x)(j=1,2),`$
where $`(e_1,e_2)`$ is the standard basis in $`^2`$ and $`(\lambda _1,\lambda _2,\mu _1,\mu _2)`$ are fixed Bloch–Floquet multipliers.

First we recall some standard definitions and formulas from the theory of theta-functions, see . Let $`\theta \left[\genfrac{}{}{0.0pt}{}{\alpha }{\beta }\right]`$ be the one-dimensional theta-function (with characteristics), defined by the following series:

$$\theta \left[\genfrac{}{}{0.0pt}{}{\alpha }{\beta }\right](z|\tau )=\underset{n}{}\mathrm{exp}(\pi i(n+\alpha )^2\tau +2\pi i(n+\alpha )(z+\beta )).$$

Notice that $`\alpha `$ and $`\beta `$ are defined modulo 1:

$$\theta \left[\genfrac{}{}{0.0pt}{}{\alpha +1}{\beta }\right]=\theta \left[\genfrac{}{}{0.0pt}{}{\alpha }{\beta }\right],\theta \left[\genfrac{}{}{0.0pt}{}{\alpha }{\beta +1}\right]=e^{2\pi i\alpha }\theta \left[\genfrac{}{}{0.0pt}{}{\alpha }{\beta }\right].$$

Later we will need the following formula which can be easily derived from the definitions:

(6.4)
$$\begin{array}{c}\theta \left[\genfrac{}{}{0.0pt}{}{\alpha _1}{0}\right](x_1|\tau )\theta \left[\genfrac{}{}{0.0pt}{}{\alpha _2}{0}\right](x_2|\tau )=\theta \left[\genfrac{}{}{0.0pt}{}{\alpha _+}{0}\right](x_+|2\tau )\theta \left[\genfrac{}{}{0.0pt}{}{\alpha _{}}{0}\right](x_{}|2\tau )\hfill \\ \hfill +\theta \left[\genfrac{}{}{0.0pt}{}{\alpha _++\frac{1}{2}}{0}\right](x_+|2\tau )\theta \left[\genfrac{}{}{0.0pt}{}{\alpha _{}+\frac{1}{2}}{0}\right](x_{}|2\tau ),\end{array}$$

where $`\alpha _\pm =\frac{1}{2}(\alpha _1\pm \alpha _2),x_\pm =x_1\pm x_2`$.

We will mostly use $`\theta \left[\genfrac{}{}{0.0pt}{}{1/2}{1/2}\right]`$ which we will denote simply by $`\theta (z)`$, which will always stand for the odd Jacobi theta function (5.4).

According to the previous section, $`\psi `$ must have the form

(6.5)
$$\psi =\frac{\mathrm{\Phi }(x_1,x_2)}{\theta (x_1|\tau )\theta (x_2|\tau )\theta (x_1x_2|\tau )\theta (x_1+x_2|\tau )}.$$

Here $`\mathrm{\Phi }`$ is nonsingular in $`^2`$. The translation properties for $`\psi `$ easily translate into the properties of $`\mathrm{\Phi }`$:

$$\mathrm{\Phi }(x+e_j)=\lambda _j\mathrm{\Phi }(x),$$

$$\mathrm{\Phi }(x+\tau e_j)=\mu _je^{3\pi i\tau 6\pi ix_j}\mathrm{\Phi }(x).$$

Standard considerations from the theory of theta-functions show that the linear space of functions with these properties has dimension $`9`$ and $`\mathrm{\Phi }`$ must be of the form

(6.6)
$$\mathrm{\Phi }=\mathrm{exp}(K_1x_1+K_2x_2)\underset{0i,j2}{}c_{ij}\theta \left[\genfrac{}{}{0.0pt}{}{i/3}{0}\right](3x_1+\gamma _1|3\tau )\theta \left[\genfrac{}{}{0.0pt}{}{j/3}{0}\right](3x_2+\gamma _2|3\tau ),$$

where $`c_{ij}`$ are arbitrary constants and parameters $`\gamma _j,K_j`$ relate to $`\lambda _j,\mu _j`$ as follows:

(6.7)
$$\lambda _j=e^{K_j},\mu _j=e^{2\pi i\gamma _j+K_j\tau }.$$

###### Remark 6.1.

The shifting

(6.8)
$$\gamma _j\gamma _j+1$$

does not change the space (6.6), the same applies to the shifts

(6.9)
$$(K_j,\gamma _j)(K_j+2\pi i,\gamma _j+\tau ).$$

Conversely, for any given $`\lambda _j,\mu _j`$ the corresponding $`(\gamma _j,K_j)`$ are determined uniquely modulo shifts (6.8)–(6.9).

Now, in accordance with Proposition 5.1, we impose the following ’vanishing’ conditions on $`\mathrm{\Phi }`$:

(6.10) $`_1\mathrm{\Phi }0`$ $`\text{for}x_1=0,`$
(6.11) $`_2\mathrm{\Phi }0`$ $`\text{for}x_2=0,`$
(6.12) $`(_1+_2)\mathrm{\Phi }0`$ $`\text{for}x_1+x_2=0,`$
(6.13) $`(_1_2)\mathrm{\Phi }0`$ $`\text{for}x_1x_2=0.`$
As we will see below, for a certain $`2`$-dimensional surface in $`4`$-dimensional space of parameters $`(k_j,a_j)`$, the conditions (6.10)–(6.13) cut a one-dimensional subspace in $`9`$-dimensional space (6.6). Thus, the corresponding $`\psi `$ will be an eigenfunction for $`L`$ automatically.

To determine the corresponding $`(K_j,\gamma _j)`$, let us rewrite $`\mathrm{\Phi }`$ using (6.4) and making identification $`\theta \left[\genfrac{}{}{0.0pt}{}{\alpha +1}{0}\right]=\theta \left[\genfrac{}{}{0.0pt}{}{\alpha }{0}\right]`$:

$$\begin{array}{c}\mathrm{\Phi }=e^{K_+x_++K_{}x_{}}\underset{0l,m2}{}\stackrel{~}{c}_{lm}\{\theta \left[\genfrac{}{}{0.0pt}{}{l/3}{0}\right](3x_++\gamma _+|6\tau )\theta \left[\genfrac{}{}{0.0pt}{}{m/3}{0}\right](3x_{}+\gamma _{}|6\tau )\hfill \\ \hfill +\theta \left[\genfrac{}{}{0.0pt}{}{l/3+1/2}{0}\right](3x_++\gamma _+|6\tau )\theta \left[\genfrac{}{}{0.0pt}{}{m/3+1/2}{0}\right](3x_{}+\gamma _{}|6\tau )\},\end{array}$$

where $`K_\pm =\frac{1}{2}(K_1\pm K_2)`$$`\gamma _\pm =\gamma _1\pm \gamma _2`$ and

(6.14)
$$\stackrel{~}{c}_{lm}=c_{ij}\text{with}il+m(mod3),jlm(mod3).$$

It is easy to see now that (6.12) leads to six linear equations on $`\stackrel{~}{c}_{lm}`$ of the form

$$\underset{l=0}{\overset{2}{}}A_l\stackrel{~}{c}_{lm}=0,\underset{l=0}{\overset{2}{}}B_l\stackrel{~}{c}_{lm}=0(m=0,1,2)$$

for certain explicitly given $`A=(A_l)`$, $`B=(B_l)`$. For generic parameters $`\gamma _1,\gamma _2`$ the vectors $`A,B`$ will be linearly independent. Therefore, these $`6`$ equations determine the $`2`$-dimensional kernel of the matrix $`\stackrel{~}{C}=(\stackrel{~}{c}_{lm})`$. Similarly, (6.13) gives six more equations for $`\stackrel{~}{c}_{lm}`$, which determine the cokernel of $`\stackrel{~}{C}`$. Thus, for any $`K_1,K_2`$ and generic $`\gamma _1,\gamma _2`$ the vanishing conditions (6.12)–(6.13) determine $`\stackrel{~}{c}_{lm}`$ (and, hence, $`\mathrm{\Phi }`$) uniquely up to a common factor. In principle, it is straightforward to write down explicit expressions for the coefficients $`c_{ij}`$ but they are cumbersome and not very useful. However, there is a better way of getting an expression for $`\mathrm{\Phi }`$, by taking a limit $`\omega 0`$ in the formula for the difference case, see (7.7) below. It turns out that $`\mathrm{\Phi }`$ has the form

(6.15)
$$\mathrm{\Phi }=\frac{e^{(k,x)}}{\theta (a_1|\tau )\theta (a_2|\tau )}\underset{0i,j2}{}b_{ij}(k_1+k_2)^i(k_1k_2)^j,$$

where the parameters $`k_1,k_2`$ and $`a_1,a_2`$ relate to $`K_j,\gamma _j`$ as

(6.16)
$$k_j=K_j+\pi i,a_j=\gamma _j(1+\tau )/2.$$

The coefficients $`b_{ij}`$ depend on $`x`$ and $`a_1,a_2`$ and are given by the following recipe. Let us introduce formal commutative variables $`A,B,C,D`$. We also need the following scalar $`\lambda `$ depending on $`\tau `$:

$$\lambda :=\theta ^{\prime \prime \prime }(0|\tau )/\theta ^{}(0|\tau ).$$

Now introduce $`U,V`$ as $`U:=A+BC`$, $`V:=A+BD`$ and put

$`b_{22}=`$ $`1,b_{21}=2U,b_{12}=2V,`$
(6.17) $`b_{20}=`$ $`\lambda +U^2,b_{02}=\lambda +V^2,b_{11}=4UV,`$
$`b_{10}=`$ $`2V(\lambda +U^2),b_{01}=2U(\lambda +V^2),`$
$`b_{00}=`$ $`(\lambda +U^2)(\lambda +V^2).`$
After that one opens the brackets, so each of $`b_{ij}`$ becomes a sum of monomials in $`A,B,C,D`$ with scalar coefficients, and then replaces each monomial using the following rule:

(6.18)
$$cA^pB^qC^rD^sc\theta ^{(p)}(x_1+a_1)\theta ^{(q)}(x_2+a_2)\theta ^{(r)}(x_1x_2)\theta ^{(s)}(x_1+x_2).$$

Here $`\theta =\theta (z|\tau )`$, as before, denotes the odd Jacobi theta function, and the upper index in brackets refers to taking derivatives in $`z`$. We treat a scalar as a multiple of $`A^0B^0C^0D^0`$ assuming, as usual, that $`f^{(0)}=f`$. To illustrate this, we present below some first of the coefficients:

$`b_{22}=`$ $`\theta (x_1+a_1)\theta (x_2+a_2)\theta (x_1x_2)\theta (x_1+x_2),`$
$`b_{21}=`$ $`2\theta ^{}(x_1+a_1)\theta (x_2+a_2)\theta (x_1x_2)\theta (x_1+x_2)`$
(6.19) $`2\theta (x_1+a_1)\theta ^{}(x_2+a_2)\theta (x_1x_2)\theta (x_1+x_2)`$
$`2\theta (x_1+a_1)\theta (x_2+a_2)\theta _1^{}(x_1x_2)\theta (x_1+x_2),`$
$`b_{20}=`$ $`\lambda \theta (x_1+a_1)\theta (x_2+a_2)\theta (x_1x_2)\theta (x_1+x_2)+\mathrm{}`$
###### Proposition 6.2.

For any $`K_1,K_2`$ and generic $`\gamma _1,\gamma _2`$ there is a unique (up to a constant factor) function $`\mathrm{\Phi }(x)`$ of the form (6.6) with the properties (6.10)–(6.13). It is described by the formulas (6.15)–(6.18) above.

To prove the formula (6.15), one goes to the limit $`\omega 0`$ in the formula (7.7) from the next section, picking up the first nonzero term (of order $`6`$ in $`\omega `$).

Formulas (6.15)–(6.18) fix the dependence of $`\mathrm{\Phi }(x)`$ on $`4`$ parameters $`a_1,a_2,k_1,k_2`$. Let us consider now the translation properties of $`\mathrm{\Phi }`$ regarded as a function of these parameters. Recall that for generic $`a_1,a_2`$ the function $`\mathrm{\Phi }`$ was determined uniquely up to a factor by (6.6) and (6.12)–(6.13). Thus, it follows immediately from Remark 6.1 that under the shifts (6.8)–(6.9) $`\mathrm{\Phi }`$ must remain the same, up to a factor independent on $`x`$. To find this factor, it is sufficient to look at the formula (6.19) for the leading coefficient $`b_{22}`$. As a result, we conclude that the function $`\mathrm{\Phi }`$, given by formulas (6.6), (6.15)–(6.18) is invariant with respect to the shifts (6.8)–(6.9):

(6.20)
$$\mathrm{\Phi }(a_1+1,a_2,k_1,k_2)=\mathrm{\Phi }(a_1,a_2+1,k_1,k_2)=\mathrm{\Phi }(a_1,a_2,k_1,k_2),$$

(6.21)
$$\mathrm{\Phi }(a_1+\tau ,a_2,k_1+2\pi i,k_2)=\mathrm{\Phi }(a_1,a_2+\tau ,k_1,k_2+2\pi i)=\mathrm{\Phi }(a_1,a_2,k_1,k_2).$$

Now let us find the interrelations between the parameters $`a_j,k_j`$ which will guarantee two remaining vanishing conditions (6.10)–(6.11). Let $`G_1`$ and $`G_2`$ denote the derivatives $`_1_2\mathrm{\Phi }`$ and $`_1_2^3\mathrm{\Phi }`$ evaluated at $`x_1=x_2=0`$:

(6.22)
$$G_1=\frac{^2\mathrm{\Phi }}{x_1x_2}(0,0),G_2=\frac{^4\mathrm{\Phi }}{x_1x_2^3}(0,0).$$

Thus defined $`G_1,G_2`$ will be regarded as functions of the parameters $`a_1,a_2,k_1,k_2`$. Consider now the following two equations on these $`4`$ parameters:

(6.23)
$$G_1=0,G_2=0.$$

It is clear that conditions (6.10) imply both of the equations (6.23). Indeed, they guarantee that the function

(6.24)
$$f(t)=_1\mathrm{\Phi }(0,t)$$

is identically zero, in particular, $`f^{}(0)=G_1=0`$ and $`f^{\prime \prime \prime }(0)=G_2=0`$.

More interestingly, (6.23) are ’almost’ equivalent to (6.10). To see this, note that the function (6.24) has the following translation properties in $`t`$:

$$f(t+1)=e^{k_2}f(t),$$

$$f(t+\tau )=e^{k_2\tau 2\pi i(3t+a_2)3\pi i\tau }f(t).$$

From (6.12)–(6.13) we know that

$$(_1\pm _2)\mathrm{\Phi }(0,0)=0.$$

This gives that $`f(0)=_1\mathrm{\Phi }(0,0)=0`$, while the first equation in (6.23) gives that $`f^{}(0)=0`$. Together with the translation properties above this implies that $`f`$ is proportional to the following theta function:

$$e^{k_2t}\theta (t+a_2|\tau )[\theta (t|\tau )]^2.$$

This function has nonzero third derivative at $`t=0`$ as soon as

(6.25)
$$k_2\theta (a_2|\tau )+\theta ^{}(a_2|\tau )0.$$

Thus, the conditions $`f(0)=f^{}(0)=f^{\prime \prime \prime }(0)=0`$ imply that $`f`$ is identically zero provided that (6.25) is true. The outcome is the following: for all $`a_1,a_2,k_1,k_2`$ satisfying the condition (6.25) the equations (6.23) imply (6.10).

Our next remark is that for any $`a_1,a_2,k_1,k_2`$ we have the relation $`_1_2^3\mathrm{\Phi }(0,0)=_1^3_2\mathrm{\Phi }(0,0)`$. This is due to the identity

$$4_1^3_2\mathrm{\Phi }4_1_2^3\mathrm{\Phi }=(_1+_2)^3(_1_2)\mathrm{\Phi }(_1_2)^3(_1+_2)\mathrm{\Phi }$$

where the right-hand side vanishes at $`x_1=x_2=0`$ due to (6.12)–(6.13). Thus, repeating the same arguments as above, we conclude that the equations (6.23) imply also (6.11), as soon as

(6.26)
$$k_1\theta (a_1|\tau )+\theta ^{}(a_1|\tau )0.$$

As a result, we see that under assumptions (6.25)–(6.26) both vanishing conditions (6.10)–(6.11) are equivalent to one system (6.23). To remove the restrictions (6.25)–(6.26), let us consider equations (6.23) in more details.

First notice that the functions $`G_1,G_2`$ share the same translation properties (6.20)–(6.21) in $`(a_1,a_2,k_1,k_2)`$ with $`\mathrm{\Phi }`$ (since differentiating in $`x`$ doesn’t affect these properties). Notice also that $`G_1,G_2`$ are polynomials in $`k_1,k_2`$. Let us pass from $`k_1,k_2`$ to another variables $`p_1,p_2`$ as follows:

(6.27)
$$p_1=k_1+\zeta (a_1),p_2=k_2+\zeta (a_2),$$

where we used $`\zeta =\zeta (z|\tau )`$ to denote the logarithmic derivative of $`\theta (z)`$:

$$\zeta (z)=\frac{\theta ^{}(z|\tau )}{\theta (z|\tau )}.$$

This is slightly different from the Weierstrass $`\zeta `$-function.

Clearly, $`G_1,G_2`$ are still polynomials in $`p_1,p_2`$ with the coefficients depending on $`a_1,a_2`$. The translation properties (6.20)–(6.21) imply that the coefficients in these polynomials will be elliptic functions of $`a_1`$ and $`a_2`$. In fact, one can write down $`G_1,G_2`$ quite explicitly. The following result follows from our calculations for the difference case from the next section.

###### Proposition 6.3.

The system $`G_1=G_2=0`$ is equivalent to the following system:

(6.28)
$$\begin{array}{c}p_1(p_2^3+3\zeta _2^{}p_2+\zeta _2^{\prime \prime })=p_2(p_1^3+3\zeta _1^{}p_1+\zeta _1^{\prime \prime }),\hfill \\ p_1(p_2^5+10\zeta _2^{}p_2^3+10\zeta _2^{\prime \prime }p_2^2+(5\zeta _2^{\prime \prime \prime }+15(\zeta _2^{})^2)p_2+\zeta _2^{\prime \prime \prime \prime }+10\zeta _2^{}\zeta _2^{\prime \prime })\hfill \\ =p_2(p_1^5+10\zeta _1^{}p_1^3+10\zeta _1^{\prime \prime }p_1^2+(5\zeta _1^{\prime \prime \prime }+15(\zeta _1^{})^2)p_1+\zeta _1^{\prime \prime \prime \prime }+10\zeta _1^{}\zeta _1^{\prime \prime }),\hfill \end{array}$$

where $`\zeta _1,\zeta _2`$ stand for $`\zeta (a_1)`$ and $`\zeta (a_2)`$ while the prime denotes taking the derivative with respect to the corresponding variable. (Notice that $`\mathrm{}(z)=\zeta ^{}+\mathrm{const}`$, so $`\zeta ^{\prime \prime }=\mathrm{}^{}`$ and so on.)

From the discussion below will follow that for generic $`a_1,a_2`$ the system (6.28) has a finite number of solutions $`(p_1,p_2)`$. Thus, we can think of (6.28) as a finite covering of the product $`\times `$ of two copies of an elliptic curve $`=_\tau =/(+\tau )`$. In fact, the only $`(a_1,a_2)`$ where the fiber is infinite are those with $`a_1=\pm a_2(\mathrm{mod}1,\tau )`$. This corresponds to the following ’vertical’ components of (6.28):

(6.29)
$$\{a_1=a_2,p_1=p_2\},\{a_1=a_2,p_1=p_2\}.$$

Another ’trivial’ component is, obviously,

(6.30)
$$\{p_1=p_2=0\}.$$

If we delete these three components from (6.28), the remaining part will be, in fact, a $`13`$-fold covering of $`\dot{}\times \dot{}`$. Since we deleted the component (6.30), the conditions (6.25)–(6.26) and, hence, (6.10)–(6.11) are satisfied on the remaining part. Thus, we arrive at the following theorem.

###### Theorem 6.4.

Let $`𝒞`$ be the finite covering of the product of two (punctured) elliptic curves which is obtained from (6.28) by deleting the components (6.29) and (6.30). Then $`𝒞`$ is the Hermite–Bloch variety for the operator (6.1) and a double-Bloch solution $`\psi (x)`$ corresponding to a point $`(a_1,a_2,p_1,p_2)`$ in $`𝒞`$ is given by the formulas (6.5), (6.15)–(6.18) and (6.27).

We still have to explain why $`𝒞`$ is a $`13`$-fold covering. To this end let us consider a family of plane rational curves $`\phi :^1^2`$ of degree $`5`$, depending on a parameter $`a=/+\tau `$ and defined as follows: if $`a`$ and $`u=(u_0:u_1)^1`$ then

$`\phi (a,u)=(\phi _0:\phi _1:\phi _2),\text{where}`$
$`\phi _0=u_1u_0^4,`$
$`\phi _1=u_1^3u_0^2+3\zeta ^{}(a)u_1u_0^4+\zeta ^{\prime \prime }(a)u_0^5,`$
$`\phi _2=u_1^5+10\zeta ^{}(a)u_1^3u_0^2+(5\zeta ^{\prime \prime \prime }(a)+15(\zeta ^{}(a))^2)u_1u_0^4+(\zeta ^{\prime \prime \prime \prime }(a)+10\zeta ^{}(a)\zeta ^{\prime \prime }(a))u_0^5.`$
Then the solutions $`(p_1,p_2)`$ of (6.28) correspond to the intersection points of two curves $`C_1=\phi (a_1,)`$, $`C_2=\phi (a_2,)`$ from our family. Namely, if $`\phi (a_1,u)=\phi (a_2,v)`$ then $`p_1=u_1/u_0`$ and $`p_2=v_1/v_0`$ obviously satisfy (6.28) and vice versa, provided $`p_1,p_20`$. We should, however, exclude from consideration points with $`p_1,p_2=\mathrm{}`$. Namely, all the curves from our family pass through $`(0:0:1)=\phi (a,\mathrm{})`$. It is a standard exercise in basic algebraic geometry to show that $`\mathrm{mult}(C_1C_2)`$ at this point equals $`12`$. In doing this local analysis, one immediately observes that the condition $`\zeta (a_1)=\zeta (a_2)`$ is necessary for $`C_1`$ and $`C_2`$ to coincide near the point $`(0:0:1)`$. This proves that the covering (6.28) is finite apart from the components (6.29). After that Bezout’s theorem tells us that the number of common points, apart from $`(0:0:1)`$, equals $`5\times 512=13`$.

###### Remark 6.5.

Operator (6.1) is symmetric under the Weyl group $`W`$ of $`B_2`$, whose $`8`$ elements act by permuting coordinates and/or changing their signs. This induces the action of $`W`$ on Bloch solutions and, therefore, on the Hermite–Bloch variety (6.28). This action, in terms of $`(a_1,a_2,p_1,p_2)`$, is generated by involutions $`(a_1,a_2,p_1,p_2)(a_2,a_1,p_2,p_1)`$ and $`(a_1,a_2,p_1,p_2)(a_1,a_2,p_1,p_2)`$.

### 6.1. Algebraic integrability

The operator $`L`$ (6.1) is completely integrable. According to Theorem 4.1, it has a commuting operator of order four,

$$L_1=_1^2_2^2+\mathrm{}.$$

Using Proposition 5.1, we see that our $`\psi `$ is a common eigenfunction for $`L,L_1`$:

(6.31)
$$L\psi =E\psi ,L_1\psi =E_1\psi ,L_2\psi =E_2\psi .$$

Here $`E,E_1`$ are some functions of the parameters $`a_1,a_2,k_1,k_2`$ which, in principle, can be calculated explicitly (though we didn’t have enough energy to perform such a calculation).

In  it was shown that apart from $`L_1`$, there is another operator

$$L_2=_1^55_1^3_2^2+\mathrm{}$$

which commutes with $`L,L_1`$ (see  for the explicit expression for $`L_1,L_2`$). The existence of a fifth order quantum integral $`L_2`$ means in this case that the Schrödinger operator $`L`$ is algebraically integrable. To see this directly, let us consider one more operator $`L_3`$, obtained from $`L_2`$ by interchanging $`x_1`$ and $`x_2`$. Then one easily checks that the common eigenspace (6.31) of $`L`$ and $`L_1`$ is $`8`$-dimensional, and for generic $`E`$ and $`E_1`$ it is spanned by the double-Bloch solutions $`\psi `$, constructed previously. On the other hand, by proposition 5.1, each $`\psi `$ will be an eigenfunction of $`L_2`$ and $`L_3`$ as well. So, the only thing to check is that the eigenvalues $`E_2,E_3`$ separate all $`8`$ solutions of the system (6.31). This is enough to check in the trigonometric limit $`\tau +i\mathrm{}`$, which is not difficult.

###### Remark 6.6.

We do not give the precise relation between the Hermite–Bloch variety and the spectral surface. To find such a relation, a careful analysis of the structure of the divisor at infinity is needed. Let us remark that in  two algebraic relations between the $`4`$ operators $`L,\mathrm{},L_3`$ were calculated explicitly. Thus, they determine a $`2`$-dimensional affine algebraic variety in $`^4`$. However, it is not isomorphic to the spectral surface as an affine variety (though they are birationally equivalent). This can be seen already in the trigonometric limit $`\tau +i\mathrm{}`$, by using the information about the spectral variety from . Namely, the results of  imply that these four operators $`L,\mathrm{},L_3`$ are not enough to generate the whole commutative ring (which is isomorphic to the coordinate ring of the spectral surface).

### 6.2. Spectrum of $`L`$

Throughout this subsection we assume that the parameter $`\tau `$ is pure imaginary, so the potential $`u(x)`$ of the Schrödinger operator (6.1) is real-valued for $`x^2`$. Its singularities is the family of lines

(6.32)
$$x_1,x_2,x_1+x_2,x_1x_2.$$

These lines cut $`^2`$ into triangles and the spectral problem for $`L`$ splits into separate spectral problems for each triangle.

Let $`\psi (x)=\psi (x;a_1,a_2,p_1,p_2)`$ be a double-Bloch eigenfunction for $`L`$, which corresponds to a point $`(a_1,a_2,p_1,p_2)`$ of the surface (6.28) in accordance with the formulas (6.5), (6.15)–(6.18) and (6.27). Given such a $`\psi `$, let us symmetrize it in the following way:

(6.33)
$$\mathrm{\Psi }(x)=\underset{wW}{}\psi (wx),$$

where $`W`$ denotes the Weyl group for the system $`B_2`$. To get a square-integrable eigenfunction, according to Proposition 5.6, one takes $`k2\pi iP`$, that is $`k=(i\pi m,i\pi n)`$ with integer $`m,n`$ having the same parity. For such $`k`$ the substitution of (6.27) into (6.28) leads to a system of (transcendental) equations on $`a_1,a_2`$, and by solving it one eventually finds the corresponding eigenfunctions. Note that because of the invariance of the system (6.28) (and of $`\psi `$) under the shifts (6.8)–(6.9), it is enough to look for the solutions $`(a_1,a_2)`$ with $`a_i`$ lying inside the fundamental parallelogram with the vertices $`\pm (1+\tau )/2`$. Also, taking into account the Weyl group action, we can restrict ourselves to the dominant weights, i.e. $`0mn`$.

Now let us consider the trigonometric limit $`\tau +i\mathrm{}`$, then one can show that for any $`m,n`$ the corresponding system (6.28) will have the unique solution $`a_1,a_2`$ inside the fundamental parallelogram (at least, for sufficiently big $`\tau `$). Moreover, if one fixes $`k=(k_1,k_2)`$ and takes then the limit $`\tau +\mathrm{}`$, then $`\psi `$ will go to the Baker–Akhiezer function $`\psi (k,x)`$ considered in  (this can be seen directly from the formulas for $`\psi `$). Now, for the Baker–Akhiezer function $`\psi (k,x)`$ it is known (see Theorem 6.7 of ) that the formula (6.33) will produces all the Jack polynomials if

(6.34)
$$k=2\pi i(\lambda +\rho )\text{with}\lambda P_+\text{ and}\rho =\frac{1}{2}\underset{\alpha R_+}{}(m_\alpha +1)\alpha .$$

For others $`k2\pi iP_+`$ the symmetrized $`\mathrm{\Psi }`$ will be zero.

In our situation this means that the function $`\mathrm{\Psi }=\mathrm{\Psi }_{m,n}`$ defined by (6.33) will be non-zero as soon as

(6.35)
$$n4m2,\text{with}nm(\mathrm{mod2}).$$

For other $`k2\pi iP_+`$ it will be zero in the trigonometric limit. But the result of  cited in Section 5.1 claims that the family of the eigenfunctions of $`L`$ is analytic in $`p=e^{\pi i\tau }`$ and specializes to the Jack polynomials at $`p=0`$. Hence, if $`\mathrm{\Psi }_{m,n}=0`$ in the trigonometric limit, it must be zero for all $`\tau `$ identically.

We conclude that the eigenfunctions $`\mathrm{\Psi }_{m,n}`$ of $`L`$ are labeled by $`m,n`$ satysfying (6.35). In particular, the ground state corresponds to $`(m,n)=(2,6)`$. The constructed solutions $`\mathrm{\Psi }_{m,n}`$ will have second order zeros along the lines (6.32) and will be invariant under orthogonal reflections with respect to these lines. According to the results of , the resulting family is complete in $`L^2(T)^W`$ (see Section 5.1 above).

###### Remark 6.7.

For certain $`(a_1,a_2,k_1,k_2)`$ the corresponding Bloch solution $`\psi `$ has a nontrivial symmetry in $`x`$, which must be a subgroup of the Weyl group $`W`$. Our formulas for $`\psi `$ do not work directly for some of these cases, because of the presence of an extra component (6.29). Of particular interest are those of the points which correspond to the solutions which are double-(anti)periodic (in each of the variables $`x_1,x_2`$). These can be viewed as multidimensional analogues of the classical Lamé polynomials .

## 7. Difference $`B_2`$ case

In this section we will generalize the results above to the following difference version of the operator (6.1):

(7.1)
$$L=a_0+a_+T_1^{2\omega }+a_{}T_1^{2\omega }+b_+T_2^{2\omega }+b_{}T_2^{2\omega },$$

where $`T_i^ϵ`$ stands for a shift in $`x_i`$ by $`ϵ`$, and the coefficients $`a_\pm ,b_\pm `$ are

$`a_\pm =`$ $`{\displaystyle \frac{\theta (x_1\omega )\theta (x_1+x_22\omega )\theta (x_1x_22\omega )}{\theta (x_1\pm \omega )\theta (x_1+x_2)\theta (x_1x_2)}},`$
$`b_\pm =`$ $`{\displaystyle \frac{\theta (x_2\omega )\theta (x_1+x_22\omega )\theta (x_1x_2\pm 2\omega )}{\theta (x_2\pm \omega )\theta (x_1+x_2)\theta (x_1x_2)}},`$
while $`a_0`$ has the form $`a_0=c_++c_{}+d_++d_{}`$ with $`c_\pm ,d_\pm `$ given by

$`c_\pm =`$ $`{\displaystyle \frac{\theta (2\omega )\theta (x_1\pm 5\omega )\theta (x_1+x_22\omega )\theta (x_1x_22\omega )}{\theta (4\omega )\theta (x_1\pm \omega )\theta (x_1+x_2)\theta (x_1x_2)}},`$
$`d_\pm =`$ $`{\displaystyle \frac{\theta (2\omega )\theta (x_2\pm 5\omega )\theta (x_1+x_22\omega )\theta (x_1x_2\pm 2\omega )}{\theta (4\omega )\theta (x_2\pm \omega )\theta (x_1+x_2)\theta (x_1x_2)}}.`$
In all formulas $`\theta (z)=\theta (z|\tau )`$ is the odd Jacobi theta function (5.4).

This is a very special case of the so-called $`BC_n`$ generalization of the quantum Ruijsenaars model . In trigonometric case it has been introduced by Koornwinder . Elliptic version was first suggested by van Diejen  and later extended in , where its complete integrability has been proven. This also can be viewed as an elliptic generalization of one of the Macdonald operators  for $`B_2`$. In what follows we assume that the parameter $`\omega `$ is generic. Note that the operator (6.1) can be restored (up to a certain gauge) in the limit $`\omega 0`$.

It is worth mentioning that coefficients of $`L`$ are not periodic, so instead of (double-) Bloch solutions one should look for the eigenfunctions in a certain $`\theta `$-functional space. Another way of putting it is to observe that $`L`$ can be reduced to elliptic form using proper gauge. For instance, consider

$$\stackrel{~}{L}=\delta ^1L\delta ,\delta =\theta (x_1)\theta (x_2)\theta (x_1+x_2)\theta (x_1x_2).$$

Then $`\stackrel{~}{L}`$ will have elliptic coefficients, so we can look for its Bloch solutions $`\psi (x)`$. Correspondingly, $`\mathrm{\Phi }=\delta \psi `$ will be an eigenfunction for $`L`$ and it will have translation properties similar to those of $`\delta `$. Abusing the language, below we refer to $`\mathrm{\Phi }`$ as a Bloch solution for $`L`$.

### 7.1. Bloch solutions

We are going to construct eigenfunctions of $`L`$ similar to the differential case above. Our ansatz for $`\mathrm{\Phi }`$ remains unchanged:

(7.2)
$$\mathrm{\Phi }=\mathrm{exp}(K_1x_1+K_2x_2)\underset{0i,j2}{}c_{ij}\theta \left[\genfrac{}{}{0.0pt}{}{i/3}{0}\right](3x_1+\gamma _1|3\tau )\theta \left[\genfrac{}{}{0.0pt}{}{j/3}{0}\right](3x_2+\gamma _2|3\tau ).$$

An analogue of the vanishing conditions (6.10)–(6.13) is dictated by the singularities of $`L`$ and is the following:

(7.3) $`\mathrm{\Phi }(\omega ,x_2)`$ $`\mathrm{\Phi }(\omega ,x_2)`$ $`\text{for all}x_2,`$
(7.4) $`\mathrm{\Phi }(x_1,\omega )`$ $`\mathrm{\Phi }(x_1,\omega )`$ $`\text{for all}x_1,`$
(7.5) $`\mathrm{\Phi }(x_1+\omega ,x_2+\omega )`$ $`\mathrm{\Phi }(x_1\omega ,x_2\omega )`$ $`\text{for}x_1+x_2=0,`$
(7.6) $`\mathrm{\Phi }(x_1+\omega ,x_2\omega )`$ $`\mathrm{\Phi }(x_1\omega ,x_2+\omega )`$ $`\text{for}x_1x_2=0.`$
We are going to show that for a certain two-dimensional variety in the space of parameters $`\gamma _1,\gamma _2,K_1,K_2`$ there is only one (up to a factor) such $`\mathrm{\Phi }`$. As a consequence, $`\mathrm{\Phi }`$ will be an eigenfunction of $`L`$, due to a natural analogue of Proposition 5.1.

We start from conditions (7.5), (7.6). Using the formula (6.4) and repeating the arguments used in case $`\omega =0`$, we obtain a linear system for $`c_{ij}`$ and can see that (for generic $`\gamma _1,\gamma _2`$) it defines $`c_{ij}`$ uniquely, up to a factor. However, solving this system leads to a very cumbersome formula. Instead, let us define $`\mathrm{\Phi }`$ by the following formula:

(7.7)
$$\mathrm{\Phi }=\mathrm{exp}(k_1x_1+k_2x_2)\underset{i,j}{}b_{ij}\theta (x_1+a_1+i\omega )\theta (x_2+a_2+j\omega )e^{\omega (ik_1+jk_2)},$$

where the summation is taken over the following set of indices:

$$(i,j)=(0,4),(4,0),(0,4),(4,0),(2,2),(2,2),(2,2),(2,2),(0,0),$$

and the coefficients $`b_{ij}=b_{ij}(x)`$ look as follows:

$$b_{ij}=\beta _{ij}\theta \left(x_1+x_2\frac{i+j}{2}\omega \right)\theta \left(x_1x_2\frac{ij}{2}\omega \right)$$

with

$`\beta _{04}=\beta _{40}=\beta _{0,4}=\beta _{4,0}=(\theta (2\omega ))^2`$
$`\beta _{22}=\beta _{2,2}=\beta _{2,2}=\beta _{2,2}=\theta (2\omega )\theta (4\omega )`$
$`\beta _{00}=(\theta (4\omega ))^2.`$
###### Proposition 7.1.

For any $`K_1,K_2`$ and generic $`\gamma _1,\gamma _2`$ there exists unique (up to a factor) function $`\mathrm{\Phi }`$ of the form (7.2) satisfying conditions (7.5)–(7.6). It is given by the formula (7.7), where $`k_j=K_j+\pi i`$ and $`a_j=\gamma _j(1+\tau )/2`$ ($`j=1,2`$).

To prove the proposition one first checks that this $`\mathrm{\Phi }`$ has the needed translation properties in $`x_1,x_2`$, then an elementary check shows that conditions (7.5)–(7.6) are satisfied.

Let us turn now to conditions (7.3)–(7.4). First, let us remark that (7.5), (7.6) at $`x=(0,0)`$ imply that

(7.8)
$$\mathrm{\Phi }(\omega ,\omega )\mathrm{\Phi }(\omega ,\omega )=\mathrm{\Phi }(\omega ,\omega )\mathrm{\Phi }(\omega ,\omega )=0.$$

Introduce now

$$G_1=\mathrm{\Phi }(\omega ,\omega )\mathrm{\Phi }(\omega ,\omega )\mathrm{\Phi }(\omega ,\omega )+\mathrm{\Phi }(\omega ,\omega ),$$

$$G_2=\mathrm{\Phi }(\omega ,3\omega )\mathrm{\Phi }(\omega ,3\omega )\mathrm{\Phi }(\omega ,3\omega )+\mathrm{\Phi }(\omega ,3\omega ),$$

and consider the system

(7.9)
$$G_1=0,G_2=0.$$

Obviously, the vanishing condition (7.3) implies (7.9). Conversely, from the system (7.9) we deduce immediately that the function

$$f(t)=\mathrm{\Phi }(\omega ,t)\mathrm{\Phi }(\omega ,t)$$

satisfies the conditions $`f(\omega )=f(\omega )`$ and $`f(3\omega )=f(3\omega )`$. Together with (7.8) this gives that $`f(\omega )=f(\omega )=0`$. Since $`f`$ is a theta function of order $`3`$, it must have the form $`f=\theta (t\omega )\theta (t+\omega )g(t)`$ with $`g(3\omega )=g(3\omega )`$. Now, since $`g`$ is a theta function of the first order with known characteristics (expressed in terms of $`a_2,k_2`$), the condition $`g(3\omega )=g(3\omega )`$ implies $`g0`$ as soon as

(7.10)
$$e^{6\omega k_2}\frac{\theta (a_23\omega |\tau )}{\theta (a_2+3\omega |\tau )}.$$

The latter condition is a difference version of (6.25). Now let us collect some corollaries of (7.5)–(7.6):

$$\mathrm{\Phi }(\omega ,3\omega )=\mathrm{\Phi }(3\omega ,\omega ),\mathrm{\Phi }(\omega ,3\omega )=\mathrm{\Phi }(3\omega ,\omega ),$$

$$\mathrm{\Phi }(\omega ,3\omega )=\mathrm{\Phi }(3\omega ,\omega ),\mathrm{\Phi }(\omega ,3\omega )=\mathrm{\Phi }(3\omega ,\omega ).$$

Thus, $`G_2`$ can be rewritten as

$$G_2=\mathrm{\Phi }(3\omega ,\omega )\mathrm{\Phi }(3\omega ,\omega )\mathrm{\Phi }(3\omega ,\omega )+\mathrm{\Phi }(3\omega ,\omega ).$$

Hence, we can repeat the same arguments with respect to $`x_2`$-variable and conclude that the system (7.9) implies (7.4) as soon as

(7.11)
$$e^{6\omega k_1}\frac{\theta (a_13\omega |\tau )}{\theta (a_1+3\omega |\tau )}.$$

Summing up, we see that the system (7.9) implies both of the conditions (7.3)–(7.4) under assumptions (7.10)–(7.11).

To get rid of restrictions (7.10)–(7.11) let us calculate $`G_1,G_2`$. First, introduce the notation $`\xi _1,\xi _2`$ for

$$\xi _1=e^{\omega k_1},\xi _2=e^{\omega k_2}.$$

A direct substitution gives that

$`\mathrm{\Phi }(\omega ,\omega )=`$ $`\beta _{0,4}\theta (4\omega )\theta (2\omega )\theta (a_1+\omega )\theta (a_23\omega )\xi _1\xi _2^3+`$
$`\beta _{4,0}\theta (4\omega )\theta (2\omega )\theta (a_13\omega )\theta (a_2+\omega )\xi _1^3\xi _2+`$
$`\beta _{2,2}\theta (2\omega )\theta (2\omega )\theta (a_1\omega )\theta (a_2+3\omega )\xi _1^1\xi _2^3+`$
$`\beta _{2,2}\theta (2\omega )\theta (2\omega )\theta (a_1+3\omega )\theta (a_2\omega )\xi _1^3\xi _2^1.`$
In a similar way we calculate $`\mathrm{\Phi }(\omega ,\omega )`$ and $`G_1=2\mathrm{\Phi }(\omega ,\omega )2\mathrm{\Phi }(\omega ,\omega )`$. As a result, the equation $`G_1=0`$ takes the following form (up to a nonessential factor):

(7.12)
$$\begin{array}{c}(\theta (a_1+3\omega )\xi _1^3\theta (a_13\omega )\xi _1^3)(\theta (a_2+\omega )\xi _2\theta (a_2\omega )\xi _2^1)\hfill \\ \hfill =(\theta (a_1+\omega )\xi _1\theta (a_1\omega )\xi _1^1)(\theta (a_2+3\omega )\xi _2^3\theta (a_23\omega )\xi _2^3).\end{array}$$

In a similar way one can compute $`G_2`$. It turns out that it is a linear combination of $`8`$ terms of the form $`\xi _1^{\pm 3}\xi _2^{\pm 5}`$, $`\xi _1^{\pm 5}\xi _2^{\pm 3}`$ and $`8`$ terms of the form $`\xi _1^{\pm 1}\xi _2^{\pm 3}`$, $`\xi _1^{\pm 3}\xi _2^{\pm 1}`$. Moreover, the combination of the last $`8`$ terms is proportional to $`G_1`$, so we can get rid of them by subtracting $`G_1`$, this does not affect the system (7.9). After these transformations equation $`G_2=0`$ takes the following nice form:

(7.13)
$$\begin{array}{c}(\theta (a_1+5\omega )\xi _1^5\theta (a_15\omega )\xi _1^5)(\theta (a_2+3\omega )\xi _2^3\theta (a_23\omega )\xi _2^3)\hfill \\ \hfill =(\theta (a_1+3\omega )\xi _1^3\theta (a_13\omega )\xi _1^3)(\theta (a_2+5\omega )\xi _2^5\theta (a_25\omega )\xi _2^5).\end{array}$$

Summarizing, we see that the system (7.9) is equivalent to the equations (7.12)–(7.13). These equations are obviously invariant under the transformations (6.8)–(6.9), in this way they define a covering over the product $`\times `$ of two elliptic curves $`=/+\tau `$. It has two ’vertical’ components

(7.14)
$$\{a_1=a_2,\xi _1=\xi _2\},\{a_1=a_2,\xi _1=(\xi _2)^1\}.$$

Another ’trivial’ component is given by

(7.15)
$$\theta (a_1+3\omega )\xi _1^3\theta (a_13\omega )\xi _1^3=\theta (a_2+3\omega )\xi _2^3\theta (a_23\omega )\xi _2^3=0.$$

After deleting these three components, one gets a finite (in fact, $`17`$-fold) covering of $`\times `$, let us denote it by $`𝒞`$. We can conclude now that for any point in $`𝒞`$ the corresponding function $`\mathrm{\Phi }`$ will satisfy the vanishing conditions (7.3)–(7.6) and, hence, it will be an eigenfunction of the difference operator $`L`$.

###### Theorem 7.2.

Formulas (7.7) and equations (7.12)–(7.13) describe the (double-) Bloch eigenfunctions of the difference operator (7.1). The Bloch–Hermite variety $`𝒞`$, which is obtained by deleting components (7.14)–(7.15) from the variety (7.12)–(7.13), is a $`17`$-fold covering of the product $`\times `$ of two elliptic curves.

As a corollary, considering the limit $`\omega 0`$ we can calculate explicitly the Hermite–Bloch variety for the operator (6.1). Namely, one picks up the terms of order $`4`$ and $`6`$ in $`\omega `$ in equation (7.12).

###### Corollary 7.3.

In the limit $`\omega 0`$ the system (7.12)–(7.13) goes to (6.28).

The only thing we still have to explain is why the degree of the covering is $`17`$. To this end, let us define a family of plane rational curves $`\phi :^1^2`$ of degree $`5`$, depending on parameter $`a`$. Namely, for $`a`$ and $`u=(u_0:u_1)^1`$ put

$`\phi (a,u)=(\phi _0:\phi _1:\phi _2),\text{where}`$
$`\phi _0={\displaystyle \frac{u_0^2u_1^3\theta (a+\omega )\theta (a\omega /2)}{\theta (a+\omega /2)}}{\displaystyle \frac{u_0^3u_1^2\theta (a\omega )\theta (a+\omega /2)}{\theta (a\omega /2)}},`$
$`\phi _1={\displaystyle \frac{u_0u_1^4\theta (a+3\omega )\theta ^3(a\omega /2)}{\theta ^3(a+\omega /2)}}{\displaystyle \frac{u_0^4u_1\theta (a3\omega )\theta ^3(a+\omega /2)}{\theta ^3(a\omega /2)}},`$
$`\phi _2={\displaystyle \frac{u_1^5\theta (a+5\omega )\theta ^5(a\omega /2)}{\theta ^5(a+\omega /2)}}{\displaystyle \frac{u_0^5\theta (a5\omega )\theta ^5(a+\omega /2)}{\theta ^5(a\omega /2)}}.`$
Then the solutions $`(\xi _1,\xi _2)`$ of (6.28) correspond to the intersection points of two curves $`C_1=\phi (a_1,)`$, $`C_2=\phi (a_2,)`$ from our family. Namely, if $`\phi (a_1,u)=\phi (a_2,v)`$ then $`\xi _1,\xi _2`$ with

$$(\xi _1)^2=\frac{u_1\theta ^2(a_1\omega /2)}{u_0\theta ^2(a_1+\omega /2)}\text{and}(\xi _2)^2=\frac{v_1\theta ^2(a_2\omega /2)}{v_0\theta ^2(a_2+\omega /2)}$$

clearly satisfy (7.12)–(7.13) and vice versa, provided (7.10)–(7.11). We should, however, exclude from consideration points with $`\xi _1,\xi _2=0,\mathrm{}`$ since $`\xi _i=e^{\omega k_i}`$. Namely, all the curves from our family pass through $`(0:0:1)=\phi (a,0)=\phi (a,\mathrm{})`$. A little difference with $`\omega =0`$ case is that now we have $`\mathrm{mult}(C_1C_2)=8`$ at this point. By Bezout’s theorem, the number of intersection points of $`C_1,C_2`$ , apart from $`(0:0:1)`$, equals $`5\times 58=17`$.

###### Remark 7.4.

Note that for given $`\xi _1^2,\xi _2^2`$ the corresponding quasimomenta $`k_1,k_2`$ seem to be non-unique, with the ambiguity of adding some multiple of $`i\pi /\omega `$. However, this would result in multiplying $`\mathrm{\Phi }`$ by $`e^{i\pi x_1/\omega },e^{i\pi x_2/\omega }`$ which are quasi-constant on the lattice $`2\omega ^2`$. Thus, this leads to the same eigenfunction, so the constructed Bloch solutions are in one-to-one correspondence with the points of the surface $`𝒞`$.

###### Remark 7.5.

It is clear that the Weyl group action on $`𝒞`$ is generated by two involutions

$$(a_1,\xi _1)(a_1,\xi _1^1)\text{and }(a_1,\xi _1)(a_2,\xi _2).$$

### 7.2. Structure of the solution space

As we mentioned above, the difference operator (7.1) is completely integrable, i.e. there exists another difference operator $`L_1`$ which commutes with $`L`$. It is given by the following expression :

(7.16)
$$L_1=c_{++}T_1^\omega T_2^\omega +c_+T_1^\omega T_2^\omega +c_+T_1^\omega T_2^\omega +c_{}T_1^\omega T_2^\omega ,$$

where the coefficients $`c_{ϵ_1,ϵ_2}`$ look as follows (we treat $`\pm `$ as $`\pm 1`$):

$$c_{ϵ_1,ϵ_2}=\frac{\theta (x_1\omega ϵ_1)\theta (x_2\omega ϵ_2)\theta (x_1+x_22\omega (ϵ_1+ϵ_2))\theta (x_1x_22\omega (ϵ_1ϵ_2))}{\theta (x_1)\theta (x_2)\theta (x_1+x_2)\theta (x_1x_2)}.$$

Now let us consider the system of two partial difference equations:

(7.17)
$$Lf=Ef,L_1f=E_1f,$$

defined on the lattice $`2\omega `$ where

$$=\{(m,n)m\pm n\}.$$

More precisely, we fix generic $`x^0^2`$ as a base point and regard a function $`f`$ in (7.17) as being defined on $`x^0+2\omega ^2`$. The base point $`x^0`$ must be outside the singular locus of $`L,L_1`$, i.e. such that $`L,L_1`$ are nonsingular on $`x^0+2\omega `$.

The Bloch solutions $`\mathrm{\Phi }`$, constructed above, are common eigenfunctions of $`L`$ and $`L_1`$. (The proof for $`L_1`$ is the same: the only thing to check is an analogue of Proposition 5.1.) Since $`L`$ and $`L_1`$ are $`W`$-symmetric, each of the $`8`$ functions $`\mathrm{\Phi }(wx)`$ ($`wW`$) will solve the system (7.17). We know that for generic point of the spectral surface $`X`$ (thus, for generic $`E,E_1`$) all $`8`$ functions $`\mathrm{\Phi }(wx)`$ are linearly independent (as functions on $`^2`$), because they have different translation properties with respect to the shifts (6.8)–(6.9). Hence, their restriction to $`x^0+2\omega `$ also gives $`8`$ linearly independent solutions of (7.17)(at least, for generic base point $`x^0`$). On the other hand, it is not difficult to see that any solution $`f`$ is uniquely determined by its values at eight points $`x^0+\nu `$ with the following $`\nu `$:

$$\nu =(0,0),(\pm \omega ,\omega ),(\omega ,\omega ),(\pm 2\omega ,0),(0,2\omega ),(\omega ,3\omega ).$$

This implies that the dimension of the solution space of (7.17) is at most $`8`$. Thus, we conclude that any solution of (7.17) (for generic $`E,E_1`$) is a linear combination of $`8`$ Bloch solutions $`\{\mathrm{\Phi }(wx)\}_{wW}`$.

###### Proposition 7.6.

The space of solutions of the system (7.17) has dimension $`8`$ and for generic $`E,E_1`$ is generated by the Bloch solutions $`\mathrm{\Phi }(wx)`$ ($`wW`$).

###### Remark 7.7.

Above we associated a double-Bloch eigenfunction $`\mathrm{\Phi }`$ to a solution $`(a_1,a_2,\xi _1,\xi _2)`$ of the equations (7.12)– (7.13). Note that these equations are invariant under $`\xi _j\xi _j`$, but this does not lead to another eigenfunction, since they will differ by a quasiconstant factor. Situation is different for the system (7.17), since it is defined on a different lattice. It is easy to see that $`(\xi _1,\xi _2)`$ and $`(\xi _1,\xi _2)`$ still lead to the same solution modulo quasiconstants, the same is true for $`(\xi _1,\xi _2)`$ and $`(\xi _1,\xi _2)`$. The resulting two functions have the same eigenvalue $`E`$ in (7.17), but opposite values of $`E_1`$. Thus, the Bloch variety for the system (7.17) is a double covering of the surface $`𝒞`$ introduced above.

We will not go into discussing the spectral properties of the difference operator $`L`$. See papers  devoted to this rather delicate matter. Let us just remark on some special solutions analogous to the ’discrete spectrum’ considered in section 6.2. Namely, let us consider the following anti-invariant solution of the system (7.17):

$$\mathrm{\Phi }_{skew}(x)=\underset{wW}{}(detw)\mathrm{\Phi }(wx).$$

The vanishing conditions (7.3)–(7.6) imply that $`\mathrm{\Phi }_{skew}`$ vanishes along lines $`x_1\pm x_2=\pm 2\omega `$ and $`x_j=\pm \omega `$. It also vanishes if $`x_1=\pm x_2`$ and $`x_j=0`$ due to anti-invariance. Let us require now for all $`8`$ functions $`\mathrm{\Phi }(wx)`$ to have the same Floquet–Bloch multipliers with respect to the shifts by $`e_1`$ and $`e_2`$, which is equivalent to the conditions

$$\mathrm{exp}(k_1)=\mathrm{exp}(k_2)=\pm 1.$$

Then $`\mathrm{\Phi }_{skew}`$ will vanish also along the shifted lines

$$x_1\pm x_2=m,m\pm 2\omega ,x_j=n,n\pm \omega (m,n).$$

In the limit $`\omega 0`$ these solutions go to those $`\mathrm{\Psi }`$ constructed in section 6.2, more precisely,

$$\omega ^6\mathrm{\Phi }_{skew}64(\theta ^{}(0))^2\theta (a_1)\theta (a_2)\theta (x_1)\theta (x_2)\theta (x_1+x_2)\theta (x_1x_2)\mathrm{\Psi }\text{as}\omega 0.$$

## 8. Hietarinta operator and its discretization

### 8.1. Continuous case

We consider now the Schrödinger operator (1.4) but first let us rescale the coordinates $`x_ia_ix_i`$, so instead of (1.4) we will consider

(8.1)
$$\begin{array}{c}L=a_1^2_1^2a_2^2_2^2a_3^2_3^2\hfill \\ \hfill +2(a_1^2+a_2^2)\mathrm{}(x_1x_2)+2(a_2^2+a_3^3)\mathrm{}(x_2x_3)+2(a_3^2+a_1^2)\mathrm{}(x_3x_1),\end{array}$$

where, as before, $`\mathrm{}(z)=\mathrm{}(z|\tau )`$ is the Weierstrass $`\mathrm{}`$-function and $`a_1^2+a_2^2+a_3^2=0`$. We are going to calculate the double-Bloch eigenfunctions of $`L`$. More specifically, we are looking for the solutions $`\psi `$ of the equation $`L\psi =E\psi `$ with the following properties:

(i) $`\psi `$ is of the form

(8.2)
$$\psi (x)=\frac{\mathrm{\Phi }(x)}{\theta (x_{12})\theta (x_{23})\theta (x_{31})}\mathrm{exp}(k_1x_1+k_2x_2+k_3x_3),$$

where $`x_{ij}:=x_ix_j`$$`\theta =\theta \left[\genfrac{}{}{0.0pt}{}{1/2}{1/2}\right]`$ and $`\mathrm{\Phi }`$ is holomorphic in $`^3`$ and depends on the differences $`x_{ij}`$ only, in other words, $`(_1+_2+_3)\mathrm{\Phi }=0`$ ;

(ii) $`\psi `$ has the following translation properties:

(8.3)
$$\psi (x+e_j)=e^{k_j}\psi (x),\psi (x+\tau e_j)=e^{\mu _j}\psi (x)(j=1,2,3),$$

where $`(e_1,e_2,e_3)`$ is the standard basis in $`^3`$.

It is not difficult to conclude that for fixed $`k_j,\mu _j`$ the conditions above determine a three-dimensional functional space, and the corresponding $`\mathrm{\Phi }(x)`$ in (8.2) must be of the form

(8.4)
$$\mathrm{\Phi }=\underset{l=0}{\overset{2}{}}c_l\theta (x_{12}+b_{12}+l\tau /3)\theta (x_{23}+b_{23}+l\tau /3)\theta (x_{31}+b_{31}+l\tau /3),$$

where $`c_0,c_1,c_2`$ are arbitrary constants and the parameters $`b_{12},b_{23},b_{31}`$ are related to $`\mu _j`$ above in the following way:

$$e^{\mu _1}=e^{k_1\tau +2\pi ib_{31}2\pi ib_{12}},$$

(8.5)
$$e^{\mu _2}=e^{k_2\tau +2\pi ib_{12}2\pi ib_{23}},$$

$$e^{\mu _3}=e^{k_3\tau +2\pi ib_{23}2\pi ib_{31}}.$$

This shows that the three-dimensional space (8.4) depends, essentially, on the pairwise differences of the parameters $`b_{lm}`$ only. Thus, without loss of generality we may assume that

(8.6)
$$b_{12}+b_{23}+b_{31}=0.$$

In formulas below we will also use $`b_{21},b_{32},b_{13}`$ under the convention that $`b_{ij}=b_{ji}`$.

Now, in accordance with Proposition 5.1, we impose certain vanishing conditions on $`\psi `$ which are motivated by the structure of the singularities of the operator (8.1). Namely, for any $`i=1,2,3`$ consider the function $`f(t)=\psi (x+ta_{i1}^2e_{i1}ta_{i+1}^2e_{i+1})`$ for $`x`$ such that $`x_{i1}=x_{i+1}`$ (we treat indices modulo $`3`$, so $`x_0=x_3`$). Our assumptions about $`\psi `$ imply that for such $`x`$ the function $`f`$ will have a pole at $`t=0`$, so its Laurent expansion will look as $`f=a_1t^1+a_0+a_1t+\mathrm{}`$. The coefficients in this expansion depend on $`x`$. Let us require that $`a_0=0`$ for all $`x`$ such that $`x_{i1}=x_{i+1}`$. Using (8.2) one rewrites this condition as follows:

(8.7) $`a_{i1}^2k_{i1}a_{i+1}^2k_{i+1}+F_i=0,`$
(8.8) $`F_i:={\displaystyle \frac{a_{i1}^2_{i1}\mathrm{\Phi }a_{i+1}^2_{i+1}\mathrm{\Phi }}{\mathrm{\Phi }}}a_{i1}^2{\displaystyle \frac{\theta ^{}(x_{i1,i})}{\theta (x_{i1,i})}}+a_{i+1}^2{\displaystyle \frac{\theta ^{}(x_{i+1,i})}{\theta (x_{i+1,i})}},`$
with (8.7) to be valid for all $`x`$ such that $`x_{i1}=x_{i+1}`$.

The following lemma follows from Proposition 5.1.

###### Lemma 8.1.

If $`\psi `$ has the form (8.2), (8.4) and satisfies the vanishing conditions (8.7), then the same will be true for its image $`\stackrel{~}{\psi }=L\psi `$ under the action of the operator (8.1).

As we will see below, for a certain three-dimensional subvariety in the space of the parameters $`k_j,b_{lm}`$ the vanishing conditions cut a one-dimensional subspace in the space (8.4). Thus, the lemma ensures that the corresponding $`\psi (x)`$ will be an eigenfunction of $`L`$.

We may regard the restriction of the expression $`F_i`$ on the plane $`x_{i1}=x_{i+1}`$ as a function of $`z=x_{i1,i}=x_{i+1,i}`$. It is easy to check then that $`F_i`$ will be an elliptic function of $`z`$ with periods $`1,\tau `$. So, first of all we have to choose the parameters $`b_{lm},c_j`$ in such a way that $`F_i(z)`$ would be non-singular. Let us assume that the parameters $`a_1,a_2,a_3`$ are generic enough, i.e. that $`a_i^2a_j^2`$. Then for $`F_i`$ to be non-singular at $`z=0`$ we need $`\mathrm{\Phi }(0)=0`$. This gives the following condition:

(8.9)
$$\stackrel{~}{c}_0+\stackrel{~}{c}_1+\stackrel{~}{c}_2=0,\stackrel{~}{c}_l=c_l\theta (b_{12}+\frac{l\tau }{3})\theta (b_{23}+\frac{l\tau }{3})\theta (b_{31}+\frac{l\tau }{3}).$$

Now we note that

$$\mathrm{\Phi }(z):=\mathrm{\Phi }_{x_{i1}=x_{i+1}}$$

is a one- dimensional $`\theta `$-function of order $`2`$, hence it has two zeros (modulo $`1,\tau `$). First zero is $`z=0`$ (due to condition (8.9)). An easy check shows that the second zero is $`z=b_{i,i1}+b_{i,i+1}`$. So, up to a constant factor,

$$\mathrm{\Phi }(z)=\theta (z)\theta (zb_{i,i1}b_{i,i+1}).$$

To get rid of a possible pole at $`z=b_{i,i1}+b_{i,i+1}`$ in (8.8) we must require that $`a_{i1}^2_{i1}\mathrm{\Phi }a_{i+1}^2_{i+1}\mathrm{\Phi }=0`$ for $`x_{i1,i}=x_{i+1,i}=b_{i,i1}+b_{i,i+1}`$. This leads to the following relation:

(8.10) $`{\displaystyle \underset{l=0}{\overset{2}{}}}\stackrel{~}{c}_l(a_1^2\zeta (b_{23}+{\displaystyle \frac{l\tau }{3}})+a_2^2\zeta (b_{31}+{\displaystyle \frac{l\tau }{3}})+a_3^2\zeta (b_{12}+{\displaystyle \frac{l\tau }{3}}))=0.`$
Here and below $`\zeta (z):=\frac{d}{dz}log\theta (z)`$. Notice that the relation (8.10) is symmetric with respect to indices $`1,2,3`$ (so we have just one condition instead of possible three!).

We use the relations (8.9), (8.10) to express (up to a common factor) the parameters $`c_l`$ in terms of the parameters $`b_{ij}`$. These relations imply that each of $`F_i`$ (8.8) is nonsingular in $`z=x_{i1,i}=x_{i+1,i}`$, therefore they are some constants depending on $`c_l,b_{ij}`$. Thus, (8.7) leads to the expressions for the differences $`a_{i1}^2k_{i1}a_{i+1}^2k_{i+1}`$ in terms of $`c_l`$ and $`b_{lm}`$. One can check that the resulting system is always compatible (i.e. that $`F_1+F_2+F_3=0`$). This follows, for instance, from the compatibility of the system (8.18) below by going to the limit $`\omega 0`$. As a corollary, the formulas (8.2),(8.4), together with (8.9), (8.10) and (8.7) deliver the expression for the double-Bloch eigenfunctions of the operator $`L`$.

Finally, let us discuss the structure of the Hermite–Bloch variety of the operator (8.1). The double-Bloch solutions are parametrized by $`b_{12},b_{23},b_{31}`$ with $`b_{12}+b_{23}+b_{31}=0`$, and the corresponding $`k_1,k_2,k_3`$ are determined from (8.7). Denote by $`b`$ and $`k`$ the three-component vectors $`b=(b_{12},b_{23},b_{31})`$ and $`k=(k_1,k_2,k_3)`$. Then two different points in the parameter space $`(b,k)`$ lead to the same solution iff the corresponding Floquet multipliers in (8.3) are the same. Taking into account relations (8.5), we conclude that the following transformations do not lead to another Bloch solution:

$`bb+\epsilon _1,kk(\epsilon _1=(2/3,1/3,1/3))`$
(8.11) $`bb+\epsilon _2,kk(\epsilon _2=(1/3,2/3,1/3))`$
$`bb+\tau \epsilon _1,kk+2\pi i(1,1,0)`$
$`bb+\tau \epsilon _2,kk+2\pi i(0,1,1).`$
Thus, $`b`$ is effectively represented by a point of a factor $`^2/+\tau `$ where $`^2=\{z_1+z_2+z_3=0\}^3`$ and the lattice $``$ is generated by $`\epsilon _1,\epsilon _2`$. This factor is isomorphic to the product of two elliptic curves with parameter $`\tau `$. Above each point $`b`$ we have a complex line of double-Bloch solutions, because equations (8.7) determine $`k`$ up to adding any multiple of $`(a_1^2,a_2^2,a_3^2)`$.

### 8.2. Discrete case

We keep the notation $`x_{ij}`$ for $`x_ix_j`$. The discrete version of the operator (8.1) looks as follows:

(8.12)
$$D=\underset{i=1}{\overset{3}{}}\frac{\theta (\omega )\theta (x_{i1,i}+\omega a_i^2)\theta (x_{i,i+1}\omega a_i^2)}{\theta (\omega a_i^2)\theta (x_{i1,i})\theta (x_{i,i+1})}T_i^{\omega a_i^2},$$

where $`a_1^2+a_2^2+a_3^2=0`$ and $`T_i^ϵ`$ stands for a shift by $`ϵ`$ in $`x_i`$. Its rational version $`\theta (z)=z`$ was communicated to us by M.Feigin who found it to be dual (in bispectral sense) to the trigonometric version $`\mathrm{}=\mathrm{sin}^2`$ of the Hietarinta operator.

The difference operator (8.12) relates to (8.1) in the following way:

$$D=a_1^2+a_2^2+a_3^2+\omega (_1+_2+_3)+\frac{\omega ^2}{2}(\mathrm{const}\stackrel{~}{L})+o(\omega ^2)\text{as}\omega 0,$$

where $`\stackrel{~}{L}`$ is gauge-equivalent to $`L`$,

$$\stackrel{~}{L}=\delta L\delta ^1,\delta =\theta (x_1x_2)\theta (x_2x_3)\theta (x_3x_1).$$

Unlike $`L`$, the operator $`D`$ is not periodic. As a result, instead of the double Bloch eigenfunctions, we will look for eigenfunctions with the translation properties similar to those of $`\delta `$. Apart from that, our ansatz for the eigenfunctions $`\phi `$ of the operator $`D`$ remains the same:

(8.13) $`\phi =\mathrm{exp}(k_1x_1+k_2x_2+k_3x_3)\mathrm{\Phi },`$
(8.14) $`\mathrm{\Phi }={\displaystyle \underset{l=0}{\overset{2}{}}}c_l\theta (x_{12}+b_{12}+l\tau /3)\theta (x_{23}+b_{23}+l\tau /3)\theta (x_{31}+b_{31}+l\tau /3),`$
(8.15) $`b_{12}+b_{23}+b_{31}=0.`$
The vanishing conditions now look as follows: for each $`i=1,2,3`$

(8.16)
$$F_i:=\theta (x_{i,i1}+\omega a_{i1}^2)T_{i1}^{\omega a_{i1}^2}(\phi )\theta (x_{i,i+1}+\omega a_{i+1}^2)T_{i+1}^{\omega a_{i+1}^2}(\phi )=0$$

identically for all $`x`$ with $`x_{i+1}=x_{i1}`$.

We have then a straightforward analog of Lemma 8.1, so the same approach as above will give us the eigenfunctions for $`D`$.

Let us first formulate the result. Namely, we consider the following two conditions on the function $`\mathrm{\Phi }(x_1,x_2,x_3)`$ given by (8.14):

(8.17)
$$\mathrm{\Phi }(\omega a_1^2,0,\omega a_2^2)=0,\mathrm{\Phi }(\omega a_2^2,0,\omega a_3^2)=0.$$

This gives us two linear equations on $`c_l`$ and we use them to express $`c_l`$ (up to a factor) through $`b_{ij}`$.

Secondly, we impose the following three relations:

$`e^{\omega a_1^2k_1\omega a_3^2k_3}=`$ $`{\displaystyle \frac{\theta (\omega a_3^2)}{\theta (\omega a_1^2)}}{\displaystyle \frac{\mathrm{\Phi }(0,0,\omega a_3^2)}{\mathrm{\Phi }(\omega a_1^2,0,0)}}`$
(8.18) $`e^{\omega a_2^2k_2\omega a_1^2k_1}=`$ $`{\displaystyle \frac{\theta (\omega a_1^2)}{\theta (\omega a_2^2)}}{\displaystyle \frac{\mathrm{\Phi }(\omega a_1^2,0,0)}{\mathrm{\Phi }(0,\omega a_2^2,0)}}`$
$`e^{\omega a_3^2k_3\omega a_2^2k_2}=`$ $`{\displaystyle \frac{\theta (\omega a_2^2)}{\theta (\omega a_3^2)}}{\displaystyle \frac{\mathrm{\Phi }(0,\omega a_2^2,0)}{\mathrm{\Phi }(0,0,\omega a_3^2)}}.`$
We use these formulas to express $`k_1,k_2,k_3`$ through $`\mathrm{\Phi }`$. The solution is not unique, and in fact we have a one-parameter family of $`k_i`$. Altogether, formulas (8.17)–(8.18) fix the dependence of $`c_l`$ and $`k_i`$ and hence of $`\phi `$ on three parameters $`b_{ij}`$ (related by (8.15)). The resulting family of functions $`\phi (x)`$ depends on three parameters: two of $`b_{ij}`$ and one more due to the freedom in resolving (8.18), see more comments below.

###### Theorem 8.2.

The formulas (8.13)–(8.15) and the relations (8.17)–(8.18) give a three-parameter family of eigenfunctions for the difference operator $`D`$.

To prove the theorem, let us first notice that each of $`F_i`$ in (8.16), being regarded as a function of $`z=x_{i,i1}=x_{i,i+1}`$, is a one-dimensional theta-function of order $`3`$, so if it doesn’t vanish, it must have three zeros (modulo $`1,\tau `$). Moreover, a simple count shows that the sum of these zeros will be equal to $`b_{i,i1}+b_{i,i+1}`$. On the other hand, a direct substitution into (8.16) shows that the relations (8.17) imply that $`F_2`$ vanishes for $`z=\omega a_1^2`$ and $`z=\omega a_3^2`$. Further, the first relation in (8.18) simply encodes the fact that $`F_2`$ vanishes at $`z=0`$. Since the sum of these three zeros is $`\omega a_2^2`$ which, generically, is not $`b_{21}+b_{23}`$, we conclude that $`F_2(z)`$ is zero identically.

At first glance it seems that we need to add four more conditions to ensure all the vanishing properties (8.16). Namely, one needs also

(8.19)
$$\begin{array}{c}\mathrm{\Phi }(0,\omega a_2^2,\omega a_1^2)=\mathrm{\Phi }(0,\omega a_1^2,\omega a_3^2)\hfill \\ \hfill =\mathrm{\Phi }(\omega a_3^2,\omega a_2^2,0)=\mathrm{\Phi }(\omega a_1^2,\omega a_3^2,0)=0.\end{array}$$

However, since $`\mathrm{\Phi }`$ depends on the pairwise differences of $`x_i`$ only, we will have that

$$\mathrm{\Phi }(0,\omega a_2^2,\omega a_1^2)=\mathrm{\Phi }(\omega a_2^2,\omega a_2^2+\omega a_2^2,\omega a_2^2\omega a_1^2)=\mathrm{\Phi }(\omega a_2^2,0,\omega a_3^2)=0.$$

In the same way other relations in (8.19) follow from (8.17).

This demonstrates that the three-parameter family constructed in the theorem satisfies the vanishing conditions (8.16), thus proving the theorem.

Finally, let us comment on the structure of the Hermite–Bloch variety. Similarly to the case $`\omega =0`$ above, Bloch solutions are parametrized by $`b=(b_{12},b_{23},b_{31})`$ with $`b_{12}+b_{23}+b_{31}=0`$. This determines the corresponding $`c_l`$ by (8.17). After that $`k=(k_1,k_2,k_3)`$ are determined from (8.18). At this point we have certain freedom: if $`k=(k_1,k_2,k_3)`$ is a solution of (8.18), then any

$$k^{}=k+\frac{t}{\omega }(a_1^2,a_2^2,a_3^2)+\frac{2\pi i}{\omega }(n_1a_1^2,n_2a_2^2,n_3a_3^2)$$

with any $`t`$ and integer $`n_1,n_2,n_3`$ will be a solution, too. However, the last term is not essential since it results in multiplying $`\mathrm{\Phi }`$ by a quasiconstant. For the same reason, the factor $`t`$ in the second term is essential modulo $`2\pi i`$ only. Besides, we still have the translation invariance of $`\mathrm{\Phi }`$ with respect to the transformations (8.11). Thus, the Hermite–Bloch variety is fibered over the product of two elliptic curves with the fibers isomorphic to $`/2\pi i`$.