Characteristic determinant and Manakov triple for the double elliptic integrable system

Using the intertwining matrix of the IRF-Vertex correspondence we propose a determinant representation for the generating function of the commuting Hamiltonians of the double elliptic integrable system. More precisely, it is a ratio of the normally ordered determinants, which turns into a single determinant in the classical case. With its help we reproduce the recently suggested expression for the eigenvalues of the Hamiltonians for the dual to elliptic Ruijsenaars model. Next, we study the classical counterpart of our construction, which gives expression for the spectral curve and the corresponding $L$-matrix. This matrix is obtained explicitly as a weighted average of the Ruijsenaars and/or Sklyanin type Lax matrices with the weights as in the theta function series definition. By construction the $L$-matrix satisfies the Manakov triple representation instead of the Lax equation. Finally, we discuss the factorized structure of the $L$-matrix.

List of main notations: q j , j = 1, ..., N -positions of particles; q j = q j − q 0 -positions of particles in the center of mass frame, q 0 = (1/N ) k q k ; q ij = q i − q j x j = e qj or x j = e 2πıqj in trigonometric or elliptic cases respectively; p j , j = 1, ..., N -the classical momenta of particles; ω = e 2πıτ -the elliptic modular parameter, controlling the ellipticity in momenta; p = e 2πıτ -the modular parameter, controlling the ellipticity in coordinates; q = e -exponent of the Planck constant; t = e η -exponent of the coupling constant; λ -the spectral parameter (1.1) (sometimes called also u); z -the second spectral parameter (1.5); A x,p -the space of operators, generated by {x 1 , .., x N , q x1∂1 , ..., q xN ∂N }; : : -normal ordering on A x,p , moving all shift of operators in each monomial to the right (2. All products of non-commuting operators should be understood as left ordered products.
1 Introduction and summary 1

.1 Brief review
The double elliptic (or Dell) model [8] is an integrable system with an elliptic dependence on both -positions of particles and their momenta. It extends the widely known Calogero-Moser-Sutherland [9,19] and Ruijsenaars-Schneider [29] families of many-body integrable systems. Historically, the model was first derived as the elliptic self-dual system with respect to the Ruijsenaars (or equivalently, p-q or action-angle) duality interchanging positions of particles and action variables [28]. At the classical level the original group-theoretical Ruijsenaars construction was not applicable to the elliptic case. Instead, a geometrical approach was used based on the studies of spectral curves and Seiberg-Witten differentials [14]. In this way the Dell Hamiltonians where proposed in terms of higher genus theta-functions with a dynamical period matrices. For this reason a definition of the standard set of algebraic tools for integrable systems (including Lax pairs, R-matrix structures, exchange relations etc) appeared to be a complicated problem. The classical Poisson structures underlying the Dell model were studied in [7,2].
An alternative version of the Dell Hamiltonians was suggested recently in [18]. The authors exploited the explicit form of the 6d Supersymmetric Yang-Mills partition functions with surface defects compactified on torus, which are conjectured to serve as the wavefunctions for the corresponding Seiberg-Witten intergable systems [25,26,27,1]. The exact correspondence of their results with the previous studies is an interesting open problem though the matching has being already verified in a few simple cases. In this paper we deal with the Koroteev-Shakirov version of the generating function for commuting Hamiltonians. Namely, for the N -body system consider the operator of complex variables: This is a definition of the infinite set of (non-commuting) operatorsÔ k . The positions of particles q i enter through x i = e q i ; t = e η -is exponent of the coupling constant η; q = e -is exponent of the Planck constant ; and ∂ i = ∂ x i , so that ∂ q i = x i ∂ i . The constant ω is the second modular parameter (controlling the ellipticity in momenta) and λ is the (spectral) parameter of the generating function. The definition of the theta-function θ p (x) with the constant modular parameter τ (p = e 2πiτ ) (controlling the ellipticity in coordinates) is given in (A.1). The commuting Hamiltonians of the Dell system were conjectured and argued to be of the form: H n =Ô −1 0Ô n , n = 1, ..., N. (1.2) Solution to the eigenvalue problem forĤ n was suggested in [3,4] by extending the Shiraishi functions [32] -solutions to a non-stationary Macdonald-Ruijsenaars quantum problem.
Our study, on the contrary does not appeal to the explicit form of the wavefunctions and is mostly focused on the generating function itself. It is based on the usage of the intertwining matrix Ξ(z) of the IRF-Vertex correspondence (see (9) for their explicit form) and the Hasegawa's factorization formula [16,17] for the gl N elliptic Ruijsenaars-Schneider Lax operator with spectral parameter z [29] The matrix Ξ(z) = Ξ(z, x 1 , ..., x N |p) enters the normalized intertwining matrix g(z, τ ) = Ξ(z) is a diagonal matrix used for convenient normalization only, see (6.10).
A key property of these matrices, which will use, is that det Ξ is proportional to the Vandermonde determinant. These intertwining matrices are known from the IRF-Vertex correspondence at quantum and classical levels [6,20,21,35]. The IRF-Vertex correspondence provides relation between dynamical and non-dynamical quantum (or classical) R-matrices as a special twisted gauge transformation with the matrix g(z), thus relating the Lax operator (1.4) with the one of the Sklyanin type [33].

Outline of the paper and summary of results
In this paper, using the Hamiltonians (1.1), we construct a generalization of the Macdonald determinant operator for the Dell system and study its applications.
We use a slightly modified and extended version of the generating functionÔ ′ (z, λ) (1.1), which depends on additional spectral parameter z, and generates an equivalent 3 set of operatorsÔ ′ k : The paper is organized as follows.
In Section 2 we derive the expression for the generalized Macdonald determinant: where q 0 is the center of mass coordinate. The determinant is well defined as the columns of the matrix commute. For the precise form of the matrix Ξ ij = Ξ i (q j , z) see (9).
In Section 3 we express the generating function (1.5) in terms of the Lax matrix of the Ruijsenaars-Schneider model: and the normal ordering is defined in (2.24). The trigonometric and rational limits (for coordinate dependence) of (1.5)-(1.8) are described as well.
In Section 4 we study the eigenvalue problem for the operatorÔ(u) (1.1) in the (coordinate) trigonometric limit p = 0, which corresponds to the dual to elliptic Ruijsenaars model, and compare our results to the known in the literature [18,3].
The main statement here is that the eigenvalues ofÔ(u) in the limit p = 0 are labelled by Young diagrams λ = (λ 1 , ..., λ N ), and have the form: (1.9) In Section 5 we study the classical limit of the Dell system. Using the classical analogue L(z, λ) of (1.8) we show that the L-matrix satisfies the Manakov triple representation [23,12] (instead of the Lax equation): The conservation laws are generated by the function det L(z, λ) only. It reduces to expression for the spectral curve of the Ruijsenaars-Schneider model in the ω → 0 limit.
In Section 6 we describe the factorized structure for the L-matrix (1.10) L(z, λ). Up to an inessential modification it is presented in the form, which is similar to the elliptic Kronecker function 4 (A.12): thus generalizing the classical version of the factorization (1.3) to the double elliptic case. The elliptic moduliτ appears as ω = e 2πıτ . It is responsible for the ellipticity in momenta, while τ controls the ellipticity in positions of particles.
We also describe connection of the L-matrix with the Sklyanin Lax operators, and propose its quantization in terms of the elliptic quantum R-matrix in the fundamental representation of GL N .
Possible applications of the obtained results and future plans are discussed in the end of the paper. Appendices contain the elliptic functions definitions and properties, description of the intertwining matrices Ξ, computations of GL 2 examples and relations between different forms of the generating functions.

Characteristic Macdonald determinant for the Dell system
In this Section we express the generating functionÔ(λ) (1.1) as a determinant of N × N matrix. The main idea is to introduce the bilinear map (2.6) | : Where A x and A p are the spaces of operators, depending only on coordinate and momenta correspondingly. And express the generating function as an image of such map, with the first argument being the Vandermonde function. Then use the property of the IRF-Vertex correspondence intertwining matrices to have determinant proportional to this function. And finally swap the orders of the calculation of the determinant with the bilinear map, using the special properties of this map and the row linearity of the determinant.
The set of intertwining matrices which we are going to use in different cases is given in the Appendix B . The elliptic coordinate case could not be treated without spectral parameter, while it is possible in the rational and trigonometric cases. The result will be proven in all the details in the trigonometric case. The rational case could be done absolutely analogously. The main statement of this Section is the following.

5)
Proof: First, notice that the above determinant is well defined since any two elements from different columns of the corresponding matrix commute. Consider the space of difference operators, generated by {x 1 , ..., x N , q x 1 ∂ 1 , ..., q x N ∂ N }. We refer to it as A x,p . Similarly, denote the spaces, generated by {x 1 , ..., x N } and {q x 1 ∂ 1 , ..., q x N ∂ N } only as A x and A p respectively.
Introduce the bilinear pairing: by defining it on the basis elements Due to x i q x j ∂ j = q x j ∂ j x i for i = j, the pairing (2.7) satisfies an important property: Then, the generating function (1.1) is represented as follows 6 : Next, we use the determinant property (2.2) and the linearity property of (2.7). From (2.9) we concludê We knew this result from [31]. and the property (2.8) provideŝ (2.14) Finally,Ô Expression under the determinant is easily calculated: This yields (2.4). Plugging the explicit expression for Ξ into the r.h.s. of (2.16) and summing up over n we get that is (2.5). Now let us write the answer for the rational case: Then the generating function O(λ) (1.1) in the coordinate rational limit : is represented as follows: Proof: The proof is word by word repetition of the trigonometric case.
These theorems generalize the generating functions for commuting Hamiltonians for the quantum Calogero-Ruijsenaars family [30,16].

The case with the spectral parameter
Let us proceed to the case with the spectral parameter. First of all we need to introduce a convenient notation. In this Section in place of the symbol θ p (x) any of the three functions could be substituted We will also use the odd theta function: For the precise relation between them, see (A.4) and (11).
In this Section, we will also need the normal ordering on A x,p , which moves all the shift operators to the right of all coordinates. Namely, on monomials where I -multi-index.
Let us formulate the main statement. Define the new generating function: (2.25) Its relation to the previous one is also explained in (11).
Then the generating functionÔ ′ (z, λ) (2.25) is represented as follows: or, equivalently, where the following expansions for the functions Ξ ij (q j , z) are assumed: for some C-numbers α and σ ij .
Proof: In order to use the trick from the Theorem 2.1 let us define the shifted Ξ matrix, calledΞ: Now a matrix elementΞ ij depends on the coordinate q j only. Therefore, the following determinants can be calculated as ordinary determinants since the elements from different columns commute. Indeed, Applying the pairing trick to them as in the proof of the Theorem 2.1, we arrive at Substituting the explicit expression for the determinant ofΞ, one obtains which, after taking the pairing equals So, one obtains: By the same argument as above, we could restore the normal ordering: Finally, by shifting the parameter z to z + N q 0 , we obtain the desired identity.

Determinant representation in terms of the Ruijsenaars-Schneider L-matrix
In this paragraph we will derive one more useful representation for the generating function. Consider (2.4 is the (quantum) trigonometric Ruijsenaars-Schneider Lax matrix. In order to witness it in its convenient form one should also perform the gauge transformation with the diagonal matrix . See details in [17,35]. With the normal ordering the gauge transformed Lax operator has the same determinant. Hence, we arrive to the following determinant representation: is the averaged sum of the Ruijsenaars Lax matrices. The averaging is over Z with the theta-function weights. Explicit form ofL RS (t, q) to be substituted into (2.44) is given in Section 3, expression (3.10). And its generation to the rational case is given by the expressions (3.13) and (3.14).
In the case with the spectral parameter the generating function depends on z. However, the arguments above could be repeated without any complications. Thus, its determinant representation is: whereL RS (z, q, t) is the elliptic Ruijsenaars-Schneider Lax matrix given by (3.2).
We present an alternative direct proof of the statements (2.43), (2.45) without usage of intertwining matrix in the next Section.

Double elliptic GL N model
The definition (1.1) can be alternatively written in terms of the standard odd Jacobi theta-function, see (11): Therefore, the HamiltoniansĤ n = (Ô ′ 0 ) −1Ô′ n also commute. Its extension to the case with spectral parameter z is given in (2.25). We are in position to represent (2.25) in terms of the (quantum) elliptic Ruijsenaars-Schneider GL N Lax matrix with spectral parameter [29,16]: Theorem 3.1. LetL RS ij (z, q, t) be the quantum Lax matrix for the elliptic Ruijsenaars-Schneider model (z -its spectral parameter), then the generating function (2.25) forÔ ′ n operators acquires the form: and substitute it into (3.3). Let us represent it as a sum of determinants. For this purpose collect all the terms with N i q n i x i ∂ i : where the matrixL RS ij (z, q i − q j , n j η, n j ) is constructed by combining rows from different terms of the sum (3.4). Using its explicit form (3.2) let us rewrite it through the elliptic Cauchy matrix: (3.9) Plugging the Cauchy determinant (A.11) into (3.9) we get (2.25).
The result of this Theorem is valid for the trigonometric and rational cases as well. The degenerations are obtained by substitutions ϑ(u) → sinh(u) → u.

Dual to the elliptic Ruijsenaars model
In the GL N case the relations (2.43)-(2.44) hold true for the Ruijsenaars-Schneider Lax matrix The intertwining matrix (B.9) provides the Lax matrix with spectral parameter (see details in [35]): To get this Lax matrix one should substitute y j = e −2q j +2q 0 +2z/N into (B.10). Then from (B.11) we obtain the following generating function of the Hamiltonians related to the Lax matrix (3.11): Being substituted into the averaged matrix it provides the following rational analogue for (1.1): In the case with spectral parameter z we deal with intertwining matrix (B.6), which leads to the following Lax operator (see details in [35]): Then The limit z → ∞ of this answer yields the expression (3.15).

Eigenvalues for the dual to elliptic Ruijsenaars model
Let us now proceed to the first possible application of our result. We are going to derive the general formula for the eigenvalues of the dual to elliptic Ruijsenaars model. It corresponds to the trigonometric limit p = 0 of the Dell system. So, the generating function looks as follows: By introducing standard notations δ j = N − j and ∆ for the Vandermonde determinant we write (4.1) asÔ By the same arguments as for ω = 0 case, we can see that the operatorÔ(u) preserves the space Λ N of symmetric functions of the variables x 1 , ..., x N . So, we can consider the eigenvalue problem for the operatorÔ(u): where Ψ is an element of Λ N . The generating function of the eigenvalues takes the form:

.4)
Proof: Let m λ be the monomial symmetric function: They form a basis in the space of symmetric functions. Following the Macdonald's book [22] let us calculate its image under the action of the operatorÔ(u): Next, make a change in the summation of variables by introducing π with the property σ ′ = σπ. Then by rearranging the factors in the product we get or in terms of the Schur functions: Recall that the Schur functions have these properties: 1) s π(λ) is either zero or equal to s µ for some µ ≤ λ; 2) their relation to the monomial symmetric functions is given by for some numbers u λµ . Therefore, the operatorÔ(u) is upper triangular in the basis {m λ }, and its eigenvalues have the form: This finishes the proof.

Eigenvalues in the GL 2 case and comparison to the known answer
Let us write down the explicit expression for the eigenvalue E 1,λ of the first Hamiltonian (1.2)Ĥ 1 = O −1 0Ô 1 . By expanding the result for the eigenvalues ofÔ(u) (4.11) in powers of u (and to the first order in ω) we obtain the eigenvalue ofĤ 1 : Up to the factor t 1 2 the results (4.11)-(4.13) coincide with those obtained in [18] and [3]. The factor t 1 2 comes from a slightly different definition of the Hamiltonians.

Classical mechanics: Manakov representation
In this section we will describe the classical limit of our construction and derive the Manakov L-A-B triple representation from it. The first step is to express the generating function of the Hamiltonians as the ratio of the two determinants. In the classical limit then, these two determinants could be combined into one, thus giving the expression for the classical spectral curve and the corresponding L-matrix.

One more generating function for the Dell Hamiltonians
So to fulfil the first step, described above, let us introduce an alternative version for the generating function of the commuting Dell Hamiltonians: put the operatorÔ (1) is also a generating function of the commuting Hamiltonians, so that commutativity ofĤ n follows from commutativity ofĤ n .
Proof: First, let us notice that the operatorsĤ k n =Ô −1 kÔ n =Ĥ −1 kĤ n also commute with each other due to commutativity ofĤ k . Therefore,Ĥ mkĤnk =Ĥ nkĤmk , or acting on this equality byÔ −1 k from the rightÔ On the one handÔ(1), compared toÔ 0 is hard to invert as its Taylor series expansion in ω starts not with 1. On the other hand the advantage ofÔ(1) is that its determinant representation, while there is no natural way to find a determinant representation forÔ 0 .

Spectral L-matrix
The operatorÔ(1) −1 in (5.1) really acts onÔ(λ) as a quantum operator, so that we can not unify the normal orderings in (5.6). At the same time in the classical limit (5.6) reduces to L RS (z, q n , t n ) (5.9) arises, which determinant H(z, λ) is the generating function of the classical Hamiltonians. They commute with respect to canonical Poisson structure Expression H(z, λ) can be considered as an analogue of the expression det(λ − l(z)) for spectral curve of an integrable system with the Lax matrix l(z). This is easy to see in the limit ω = 0. Due to (A.5) we have L(z, λ)| ω=0 = 1 N − λL RS (z, q, t) , (5.11) where 1 N is the identity N × N matrix. Plugging (5.11) into (5.8) we get Therefore, equation H(z, λ)| ω=0 = 0 is indeed the spectral curve of the elliptic Ruijsenaars-Schneider model (written in some complicated way). In the general case L(z, λ) is not a Lax matrix. Its eigenvalues do not commute with respect to (5.10).
Let us remark that the existence of the Manakov's L-A-B triple does not contradict the possible simultaneous existence of some Lax pair. If we had a true Lax matrix for the Dell model then det L(z, λ) should represent its spectral curve. So, if the Lax representation exists we need to find a matrixL of a size M × M (as was mention in [8] it is natural to expect M = ∞) and a change of variables u = u(z, λ), ζ = ζ(z, λ) satisfying Another comment is about the geometrical meaning of L(z, λ) (5.9) and the L-matrix (5.8). In the Krichever-Hitchin approach to integrable systems these matrices are sections of (the Higgs) bundles over a base spectral curve with a coordinate z. The classical analogue of the Ruijsenaars-Schneider Lax matrix (1.4) takes the form c -the "speed of light" constant of the classical Ruijsenaars model. Its quasi-periodic behaviour is as follows: L RS (z + 1) = L RS (z) and The first factor in (5.15) means that all terms in the sum (5.9) have different quasi-periodic behaviour on the lattice of periods 1, τ . Therefore, L(z, λ) is not a section of a bundle over the elliptic curve. This can be easily corrected by the substitution Then the first factor in (5.15) get canceled, and we come to the quasi-periodic matrix L(z, λ). The L-matrix (5.8) is quasi-periodic as well. The price for this change of variables (5.16) is as follows.
Initially we had the matrix L(z, λ) to be not quasi-periodic, but with a single simple pole at z = 0. After the change of variables we come to a quasi-periodic matrix function, but having higher order poles. The terms with positive n in the sum (5.9) acquire the n-th order pole at z = η, and the terms with negative n acquire the −n + 1-th order pole at z = 0. In what follows we do not use the substitution (5.16) keeping in mind that it can be done.

L-A-B triple
Consider the L-matrix (5.8) It is easy to see that this matrix satisfies identically the following equations known as the Manakov representation [23]: where trB k (z, λ) = 0 . (5.18) and the "time" derivatives will be specified later. Indeed, by differentiating (5.8) we get and The property (5.18) follows from The l.h.s. of (5.21) equals zero since det L(z, λ) = H(z, λ), while the r.h.s. equals to the trace of B k (z, λ) (5.20). Alternatively, introduce Then it follows from the conservation of det L(z, λ) that trM k (z, λ) = trM k (z, 1) .

Factorization of the L-matrices
In this section we will show how our result naturally embeds the Dell model into the standard factorized Lax matrix approach, to the description of which we proceed in the next paragraph.

Classification of the factorized L-matrices
In [8] integrable many-body systems of the Calogero-Ruijsenaars family were naturally classified by the types of dependence on the coordinates and/or momenta. Both types are numerated by three possibilities. Each can be either rational, trigonometric or elliptic. For example, the choice (rational coordinate, trigonometric momenta) corresponds to the rational Ruijsenaars-Schneider model, while the choice (rational momenta, trigonometric coordinate) imply the trigonometric Calogero-Sutherland system. In the coordinate part this classification follows from solutions of the underlying functional equations [9,29] -the Fay identity (A.14) and its degenerations. By interchanging the types of coordinate and momenta dependence (as in the example pair above) one gets a pair of systems related by the Ruijsenaars duality transformation [28]. When both types coincide the corresponding model is self-dual. These are the rational Calogero-Moser system, the trigonometric Ruijsenaars-Schneider model, and finally the double elliptic model, which existence was predicted by these arguments.
Here we supply the upper classification with precise substitutions corresponding to the factorized Lax (or Manakov) L-matrices. As was discussed in [35] the factorized Lax matrices (with and without spectral parameter) for the systems of Calogero-Ruijsenaars type can be specified by a choice of two ingredients: the function f , and the intertwining matrix Ξ(z): where the matrix G(z) is defined in terms of Ξ(z): with some diagonal matrix D (see (6.10) below) and c -the light speed parameter 7 in the Ruijsenaars-Schneider model. The function f (w) is either: 2) exponent: f (w) = e w (6.4) The first choice of the function f being substituted into (6.1) provides the Lax matrix of the Calogero-Moser-Sutherland systems [9,19]. The second choice of f gives rise to the Lax matrices of the Ruijsenaars-Schneider models [29]. The choices of Ξ(z) in (6.1) are given by (B.12) in the elliptic case, in the trigonometric case it is (B.9), and in the rational case it is (B.6). In trigonometric and rational cases one can also use the Vandermonde matrices (B.3) and (B.1) as Ξ ij (z) = (z − q j ) i−1 and Ξ ij (z) = (e z x j ) N −i respectively. The spectral parameter is cancelled out in these cases, and we get the Lax pairs of the Calogero-Ruijsenaars models without spectral parameter. See the review [35] for details.
Based on (1.1) and the Manakov L-matrix structure (5.8)-(5.9) we come to elliptic version for the function f : The latter result is explained in the next paragraph in detail. A more universal classification picture arises if we slightly change the definition of f λ (w) as in transition from θ p (e w ) to ϑ(w) (A.4) together with additional normalization factor ϑ ′ (0)/ϑ(log(λ)). Then the function f λ (w) turns into the Kronecker elliptic function [36] depending on the moduliτ (defined through ω = e 2πıτ ): where u = log(λ). In trigonometric and rational limits (when Im(τ ) → 0) it is as follows: This function was used by I. Krichever [19] to construct the Lax representation with spectral parameter for elliptic Calogero-Moser system. It is widely used in elliptic integrable systems due to an addition Theorem known also as the genus one Fay identity (A.14). Considered as functional equation its solutions (including degenerated versions) were extensively studied [9].
In this way we come to the classification for the function f u (w) (responsible for the momenta type dependence), which is parallel to the well known classification of the coordinates dependence without spectral parameter [9] and with spectral parameter [19].

Factorized structure of the Dell L-matrix
Recall that the Lax matrix of the Ruijsenaars-Schneider model (3.2) is factorized as follows (1.3): where g(z) = g(z, τ ) = Ξ(z)D −1 (6.9) with the intertwining matrix Ξ ij (z) (B.12) and the diagonal matrix A conjugation with the latter diagonal matrix D is performed in order to have a convenient form for L RS (z). Consider the matrix G(z) (6.2). The Ruijsenaars-Schneider Lax matrix (6.8) takes the form up to gauge transformation with the diagonal matrix exp ( z N cη P ): Let us proceed to the double elliptic case. Plugging (6.8) into the matrix L(z, λ) (5.9) we get Ncη P L(z, λ)e z Ncη P = G −1 (z)θ ω λAd e −Nη∂z G(z) . (6.14) By introducing also we come to the following expression for the Manakov L-matrix (5.8): It is also gauge equivalent to bothL(z, λ) and L(z, λ). In terms of the Kronecker function (6.6) we may write the Manakov L-matrix aš where u = log(λ). From the point of view of the classification of the factorized L-matrices discussed above, the expressions (6.15) and (6.17) can be considered as matrix analogues for theta-function and the Kronecker elliptic function respectively. The Kronecker function in (6.17) is constructed by means of the theta-operator Θ(z, λ) understood in a "plethystic sense" (6.15).

Relation to Sklyanin Lax operators
Due to the IRF-Vertex relation one can equivalently use the Sklyanin type Lax operators [33] instead of the Ruijsenaars-Schneider one in (1.8). Consider the gauge transformed Ruijsenaars Lax matrix (6.8): It is the classical analogue for the representation of the quantum Sklyanin Lax operator [17]: Consider first the classical case (6.18). Its generalization to the double elliptic model is performed similarly to the previous paragraph. Define So, it could be expressed as a sum of the Sklyanin type Lax matrices with different coupling constants. For each of them, we know from [21] that it can be alternatively represented in terms of the underlying elliptic R-matrix:  Consider the matrix operator :L Skl (z, λ) : (6.19). In the special case η = − /c the Sklyanin Lax operator is represented as the quantum Baxter-Belavin R-matrix in the fundamental representation [5,17]. The R-matrix coefficients are of the form: where using the definition (A.6) we introduced (6.28) In this way we come to the matrix representation forL Dell (z, λ) (6.26): L Dell (z, λ) = R 12 (z, λ) = R 12 (z, 1) −1 R 12 (z, λ) ∈ Mat(N, C) .

Discussion
• The Manakov L-matrix (6.24) can be used as a building block for the (partial) monodromy as it happens in integrable chains [12]. Namely, construct monodromy for a chain of length L: with the L Skl (6.24), and the L-matrix at each site depends on its own set of canonical variables. Then, det T (z, λ) provides a product of generating functions for each site. In [18] this was called the spin generalization of the double-elliptic model. Besides this construction one can also study the averaging of the spin Ruijsenaars-Schneider Lax operators.
• The quantization of (7.1) can be performed in a usual way with R(z, λ) defined in (6.29)-(6.31). Presumably, a properly defined quantum determinant of T(z, λ) could be a generating function of commuting operators. Notice that, we did not prove the Yang-Baxter equation for R(z, λ) (it is rather not fulfilled since the traces of T (z, λ) do not commute), so that R(z, λ) is not an R-matrix. Finding equation for R(z, λ) is another interesting problem.
• Orthogonality of the eigenvectors. To proceed further in our construction along with the analogous treatment for the usual Macdonald's symmetric functions we need to know the analog of the Macdonald's measure, with respect to which the operatorsĤ n are self-conjugate. Probably, expressing the eigenvectors ofĤ n as vector-valued characters of some elliptic algebras might work, in the analogy with the paper [11].
• The integrable many-body systems can be also described via commuting differential or difference KZ-type connections or Dunkl operators using the Matsuo-Cherednik (or Heckman) projections respectively. The Ruijsenaars duality in many-body problems then turns (or rather embeds) into the spectral duality interchanging canonical coordinates and momenta being written in separated variables of the corresponding Gaudin models and/or spin chains. Much progress has been achieved in studies of these relations including its elliptic version [24]. An interesting problem is to find the Dunkl-Cherednik like description for the double-elliptic models (and define a doubleelliptic version of the qKZ equations). We discuss these topics in our forthcoming paper [15].
• Let us mention a recent paper [34], where an elliptic integrable many-body system of Calogero type with the Manakov representation was obtained instead of the Lax representation. It was derived through a reduction from an integrable hierarchy. Presumably, it is a special limiting case, which can be deduced from the double-elliptic Manakov L-matrix.
• Further study of algebraic structures underlying the double-elliptic model is another set of important problems. This includes r-matrix structures at classical and quantum levels (RTT relations) and extensions of the Sklyanin quadratic algebras by means of R-matrix type operators (6.29)-(6.30).
• The determinant formula (1.7)-(1.8) can be naturally extended to the cases of many-body systems associated to the root systems of BC N type by substituting the Lax operators for the Ruijsenaars-Schneider-van Diejen type. This can be performed using results of [13]. We will discuss it in our future publication.
• Another open problem is to describe the classical L-matrix (5.8)-(5.9) in the group-theoretical (or Krichever-Hitchin) approach. As we mentioned in (5.16), one way is to consider the matrix valued function with higher order poles at a pair of marked points. Another possibility is to consider (block-matrix) direct sum k∈Z L RS (z, q k , t k ) embedded into GL(∞). Each block is well defined, and the weighted sum of all blocks can be viewed as a matrix valued character.
We will discuss these questions in future publications.

Appendix A: Elliptic functions
We use several different theta-functions. The first is the one used in [18]: where the moduli of the elliptic curve τ ∈ C, Im τ > 0 enters through Another theta-function is the standard odd Jacobi one: They are easily related: In the trigonometric limit p → 0 The Riemann theta-functions with characteristics are defined as follows: In particular, where η D (τ ) is the Dedekind eta-function: The determinant of the elliptic Cauchy matrix is given by Define the elliptic Kronecker function and the first Eisenstein function [36] They satisfy the (genus one) Fay trisecant identity and its degeneration Also, (A. 16) 9 Appendix B: Intertwining matrices Following [35] we consider the intertwining matrices in two cases: with a spectral parameter and without a spectral parameter. These matrices lead to the Ruijsenaars-Schneider Lax matrices via (1.3) with and without spectral parameter respectively.
The cases without spectral parameter. Here we deal with the Vandermonde matrices in the rational and trigonometric coordinates.