Discrete Versions of Some Classical Integrable Systems and Factorization of Matrix Polynomials

. Discrete versions of several classical integrable systems are investigated, such as a discrete analogue of the higher dimensional force-free spinning top (Euler-Arnold equations), the Heisenberg chain with classical spins and a new discrete system on the Stiefel manifold. The integrability is shown with the help of a Lax-pair representation which is found via a factorization of certain matrix polynomials. The complete description of the dynamics is given in terms of Abelian functions; the flow becomes linear on a Prym variety corresponding to a spectral curve. The approach is also applied to the billiard problem in the interior of an ΛΓ-dimensional ellipsoid.


Introduction
In this paper we study a class of discrete integrable systems which are closely related to problems occurring in mathematical physics such as the Heisenberg model for classical spins or the billiard problem in the interior of an ellipsoid. A discrete system can be viewed as the iterates of a symplectic mapping, the time t EΈ being the number of iterations. Such a system will be called integrable if it possesses sufficiently many integrals which are in involution with respect to a symplectic structure.
To describe such a discrete system we take as starting point a variational principle for a functional S = S(X) defined on the space of sequences X = (X k ), keΈ by a formal sum s= Σ &{x k ,x k+1 ). keZ Here X k are points on a manifold Jί n and if is a function onQ 2n = Jί n x J( n . The Euler-Lagrange equation of such a functional are second order difference equations (see Sect. 1) and keZ plays the role of the discrete time.
This description is to be viewed as the discrete analogue of the Hamilton principle δS = 0 for S=i&(q 9 q)dt, and the related symplectic flow can, as usual, be defined via the Legendre transform provided det(c£? 44 )φ0. Similarly, one can define a symplectic structure on Q 2n under an appropriate nondegeneracy assumption (see Sect. 1, Sect. 2, (8) and [1]). We will call such a system "integrable" if there are sufficiently many integrals which are in involution with respect to this symplectic structure.
As an example of this setup we mention the Heisenberg chain with classical spins, where and where J is a symmetric matrix, which we may take to be diagonal.
The corresponding chain of quantum spins \ (so-called XYZ Heisenberg model) was investigated by Faddeev and Takhtajan [2] in the framework of the quantum inverse scattering method, using the fundamental results by Baxter [3]. As it was shown by Pokrovsky and Khokhlachev [4], the problem of finding some special eigenfunctions in the quantum XYZ-model leads to the stationary equation δS = 0 for the Heisenberg chain with classical spins. For this discrete system Granovsky and Zhedanov [5,6] found two algebraic integrals and special solutions. The integrability of these systems, even for arbitrary dimension n of the sphere: M = S n , was shown by one of the authors ( [7], see also [1]), where the general solution was described in terms of ^-functions, generalizing the connection between the spectral theory of one-dimensional Schrόdinger operators and the classical Neumann systems derived in [8,9]).
The main problem to be discussed in the first part of this paper is a chain of JV(iV-l) orthogonal matrices: We take Ji n = O{N\ n= , and where J 1 is a positive symmetric matrix. This problem was introduced in [1] where it was shown that in the continuous limit this problem leads to the Euler problem for force-free motion of a rigid body as it was generalized by Arnold to arbitrary dimensions. Alternately, one could use the positive Lagrangian ^tr((X-Y)J(X -Y) τ ). Indeed, for X, YeΘ(N) this expression agrees with trJtΐ(XJY τ ), i.e. differs from -JSf (X, Y) only by a constant.
For the cases JV = 3, JV = 4 this system was shown to be integrable by explicit construction of commuting integrals [1]. In Sect. 3 we will establish the integrability of this system for all JV by using the discrete version of the isospectral technique, which leads to the complete description of the dynamics of this system in terms of Abelian functions. The flow is quasi-periodic in the discrete time parameter k and is linear on a Prym variety.
We indicate the approach underlying the solution of this problem. It is based on the construction of an isospectral mapping on a class of matrices (L) into itself, analogous to the Lax approach in the continuous case in which the differential equation is cast in the form L=\L,A~] for some class of linear operators (L).
Finding the class of matrices (L) and the isospectral deformation is a hit-ormiss game and depends on good guesses. In the discrete case it turns out to be connected with a factorization of matrices. We recall the beautiful observation by Symes [11] that the QR-algorithm of Jacobi matrices is closely linked to the Toda flow: The QR-algorithm, an important device in numerical analysis for diagonalization of matrices, consists in factoring a real nondegenerate matrix L into a product L=QR of an orthogonal matrix Q and an upper triangle matrix R with positive diagonal elements. Now the mapping gives rise to an isospectral map since L ί = Q~1LQ. It was Symes' observation [11] that the application of this process to L=expK, where K is a symmetric tridiagonal (or Jacobi-) matrix leads to an integrable mapping which is interpolated by the Toda flow. This idea has been extended to more general classes of matrices by Deift et al., see [12].
Although this mapping has no variational description, this idea turns out to be fruitful also for the problems at hand. The main new feature is that we start with a class of certain quadratic matrix polynomials and a suitable factorization 1 We denote the "moment of inertia" by J and not, as is customary, by /, to avoid confusion with the identity matrix If such a factorization exists and can be defined in a unique way then it gives rise to an isospectral mapping Liλj^L^λ) by exchanging the factors: This can be viewed as a discrete analogue of the Lax representation.
The main difficulty is, of course, to find the class of matrix polynomials together with a factorization in such a way that it corresponds to the dynamics of the given problem. Moreover, even if one has such a factorization it is generally not unique and the above procedure gives rise only to a correspondence, i.e. a multiple valued mapping; this is in good agreement with the fact that frequently the difference equation δS = 0 gives rise to such correspondences.
In the problem of orthogonal chains we are able to describe such a class of matrix polynomials and a corresponding factorization which can be made unique by specifying a suitable splitting of the spectrum. From this Lax representation we will find the integrals as well as the algebraic curve on whose Jacobian variety the flow becomes linear in k. A general theory of factorization of matrix polynomials can be found in [13]. The special factorization derived here (see Sect. 1) uses standard ideas from [13]; it involves orthogonal matrices and may be of interest in itself. Using the ideas of "finite-gap" integration [14] in the matrix case developed by Dubrovin [15] we exhibit explicit formulas for the dynamics in the case n = 3 in terms of classical elliptic functions.
In Sect. 2 we discuss some generalizations of this system, including the above mentioned Heisenberg model with classical spins. We take M as the manifold of rectangular nxN matrices X, nrg JV for which XX τ = I n and define the Lagrange function by where J is a symmetric NxN matrix. In other words M is the Stiefel manifold V n N of orthonormal n-frames in RΛ For n = l, N = 3 this represents the Heisenberg model and for n = N we obtain the chain of orthogonal matrices described above. For n = 1 we show how the factorization procedure leads to the hyperelliptic curves and formulas involving ^-functions as in [1]. In the general case ί<n<N we exhibit the corresponding factorization problem without full treatment.
The last section (Sect. 3) is devoted to the billiard problem in the interior of an ellipsoid in KΛ This problem also fits into the above framework. The relevant class of matrices for the Lax representation agree with those introduced in [10] for the study of the geodesic flow on an ellipsoid -probably the oldest nontrivial integrable system in arbitrary dimensions. This class of matrices fits into the procedure of Sect. 1. Finally, we will establish a connection between this billiard problem and a discrete version of the Neumann problem, where a certain symmetry of this system, found in [1], will play a crucial role. In the continuous case such a connection between the geodesic flow on the ellipsoid and the Neumann system was discovered by Knόrrer with the help of the Gauss mapping of the ellipsoid [16]. However, the connection described here is of a different nature.
In the above discussion we referred frequently to a discrete version of a continuous system, as in the case of the billiard problem inside an ellipsoid and the geodesic flow on the same ellipsoid -whose orbits are obtained as limits of tangential billiard shots. Another example is the above mentioned chain of orthogonal matrices and the corresponding continuous system of the force free top. Without trying to be precise we require in all these cases a) that the discrete system tends to the continuous system under a limit process and b) that both systems are integrable and are given by "natural" variational problems.
As a rule it is easy to go from a discrete system to a continuous one without destroying integrability. However, the converse is much more difficult, as is often the case if one wants to preserve some structure under discretization of a continuous system. Of course, one could take the "time ε" mapping of a flow but that is usually not described by a simple variational problem. The difficulty is to preserve the integrability under discretization. In this sense our method may be of interest since it provides an approach to the construction of an integrable discretization for the continuous system with known Lax representation, polynomially depending on an additional "spectral" parameter λ. The importance of representations of this type became clear after the paper of Novikov [17]. They exist for most of the known integrable hamiltonian systems and are related to the theory of Lie algebras (see [18,19]  In a forthcoming paper [29] it is shown that the factorization of certain linear (!) matrix polynomials, introduced in [10], leads to the billiard problems in domains on the sphere and the Lobachevsky space, bounded by conic sections.
The present paper was completed in February 1989 and has been circulated as a preprint of the Forschungsinstitut fur Mathematik Zurich. For various reasons its publication has been delayed. We added some relevant new references at the end of this revised paper. In particular, we draw attention to the interesting work by Deift, Li, and Tomei [32] in which the systems considered in the present paper are related to loop groups. Moreover, it is shown how the discrete mappings considered here can be interpolated by integrable Hamiltonian flows.

The Equations of "Motion"
We consider the functional S(X\ determined by a formal sum on the sequences X = (X k ) with X k e0(N% i.e. orthogonal N by N matrices. The stationary points of S are described by the equation δS = 0 or (2) where Λ k = AI is a matrix Lagrange multiplier, which is determined in such a way that X k X k =L Λ k is uniquely determined by X k -ί9 X k , X k +i> but as we will see later, not uniquely by X k -U X k \ it is a complicated function of X k -l9 X k .
Therefore we use the Euler description of the dynamics. This can be done in the following way. Rewrite (2) as Xl^. (3) Introducing m k = X k JX\γ -X k _ γ JX\ we see that the last Eq. (3) means that m fc+1 -m k-The conservation of m h which is the discrete analogue of the angular momentum in space [20] is the consequence of the left-invariance of £?(X, Y) (see [1]). In the variables fixed relative to the body we have the "angular velocity" ω k = X k Γ X kί =X k~ί X kί €θ(N) and "angular momentum with respect to the body" M k = X k } 1 m k X k^1 =ω k J-Jω k eo(N)* and thus Eq. (3) can be rewritten as a "discrete Euler-Arnold equation" [1] M k+ί =ω k M k ω k~19 Jω k , ω k e0(N). l ; In the continuous limit when X k = X(t k ) 9 t k = t 0 + kε, ω k = X k i X k _ 1 πI -εΩ(t k \ ω k = X k 1 X kί and M k = ωlJ~Jω k^ε (JΩ + ΩJ) = εM(tk k % M = JΩ + ΩJ, this Eq. (4) becomes the usual Euler-Arnold equations for the motion of the JV-dimensional rigid body Ωeo(N). (5) The main new feature of the discrete system (4) is the connection between M and ω: which we need to solve to find ω. In fact such ω is not unique (see below) and ω k+ι is not uniquely determined by (4), which therefore leads to a correspondence and not to a mapping. We discuss the symplectic properties of this correspondence (see also [1]). The Eq. (2) is a particular case of the Lagrangian equations <5*S = 0 for the functional S= Σ JS?(X*,X k+1 ), keZ (see the Introduction), which in an appropriate coordinate system (x, y) on Q 2n = Jl n x Jί n can be written as (X k ,X k+i n(X^uX k ) = 0.
The submanifold Γ 2n in Q 2n x β 2n , defined by determines generally some correspondence between subsets of Q 2n . On Q 2n one can define a closed 2-form σ by

d2j ? A A σ = -z-Γ-ax A ay
The submanifold Γ 2n is isotropic for the form σ'-σon Q 2n x Q 2n . Indeed We see that if is the generating function of the mapping, determined locally by (8) in the domain of nondegeneracy of σ: det dxdy *0, and therefore this mapping is symplectic with respect to the symplectic structure σ.
In this connection it is useful to introduce the discrete version of the Legendre transformation τ of Q 2n into T*Jί. It is defined by where p is the fiber coordinate and a = pdx the standard 1-form on T*Jί.
Thus the pullback of this form is τ*a = β = £? x dx, and the standard symplectic form da on T*Jί pulls back to τ*da = dβ = σ, which is nondegenerate whenever τ is noncritical. Generally, τ has only a local inverse.
We where S==S T is so chosen that X T P is skew symmetric, i.e.

P = ±{YJ-XJY T X).
The standard 1-form a = tv(P τ dX) is taken into β = tτ(dXJY τ \ since X τ dX is skew symmetric. The standard sympletic form da on T*O(N) is mapped into We record that this 2-form σ is the pullback of the standard symplectic form on T*O(N) under τ; it is nondegenerate at all noncritical points of τ. If we define locally a mapping φ: (X, Y)->(X f , Y') by selecting a branch of the correspondence where X, Y 9 X', Y'eO(N), then this mapping preserves σ. Equivalently, the mapping τφτ~ί, locally defined near regular values of τ, preserves the standard symplectic structure as well as the corresponding Poisson structure on T*O(N).
Here it is crucial to solve the matrix equation M = ω τ J-Jω for ωeO(N\ a question which will be discussed completely in Subsect. 1.2. Now it is well known that the reduction of T*0(N) to o*(iV) takes the standard Poisson structure of T*0(N) (up to a constant nonzero factor) into the Lie-Poisson structure on o*(JV) which we identify again with o(N); here f M denotes the skew-symmetric matrix of partial derivatives df/dM^. This proves that the mapping ψ: M-*M' of (4) preserves this Poisson structure, which agrees with the Poisson structure preserved by the usual continuous rigid body motion given by (5). This reduction is the discrete version of the well known reduction procedure [20] for Hamiltonian systems. For the derivation of (11) see also [30,31].
Our next goal is to show that this mapping is integrable, i.e. preserves sufficiently many functions F i9 which are in involution with respect to the above Poisson structure. As a matter of fact it turns out that these "integrals" have the same form as in the continuous case which are known to be in involution.

The Solution of the Matrix Eq. (6): ω τ J -Jω -M
We have to solve two crucial problems: a) to define the mapping φ in a unique way by selecting a branch of the correspondence and b) to verify that this mapping is integrable. Both these problems can be reduced to an appropriate factorization problem for a matrix polynomial, as we will show now.
The first problem, to construct a well defined map φ:(X 9 Y)->(X', Y') whose graph belongs to the correspondence (9) reduces to finding a well defined solution ωeO(N) of the matrix equation: for a given skew-symmetric matrix M. Indeed, setting ω=Y τ X, M = ω τ J-Jω, This problem is, in fact, equivalent to finding a definite inverse for the Legendre transformation τ:(X, Y)-+(X,P) 9 since Hence a solution ω of this equation gives rise to Y=Xω τ , thus defining τ~ι.
The crucial observation is contained in the following The proof is an obvious verification, which shows also that the solution ω is necessarily an orthogonal matrix. It turns out that the choice of the solution ω is fixed by the corresponding factorization of the determinant We will prove below Theorem 1. Assume that for the real skew symmetric matrix M the polynomial P(λ) = P( -λ) admits a splitting with a real polynomial p(λ) satisfying for all λe<E then there exists a unique matrix ωeO(N) satisfying (12) and We postpone the proof of this theorem to Subsect. 1.4. We discuss the splitting of the determined P(λ). Since M + M Γ = 0 one has P{λ) = P(-λ) and the set Σ of all roots of P(λ) satisfies Σ= -Σ. The factorization (13) corresponds to a splitting Σ = Σ + uΣ_ into disjoint sets Σ+, Σ_ satisfying where Σ + is the zero set of the real polynomial p(λ). Here we denote for any set A C C by A the set of all ά,aeA and by -A the set of (-a\ a e A. Any such splitting gives rise to such a factorization (12) and thus to a solution of (6). Obviously the possibility of such a factorization requires that P(λ) has no roots on the imaginary axis. In this case one factorization is obtained by taking for Σ + the roots of P(λ) in the right half plane, and £"_ = -Σ+.
We give an outline of the proof of Theorem 1 under the assumption that the roots of P(λ) are distinct, leaving the complete proof of the general case for later.
Then there exist eigenvectors ψ k :(l -λ k M -Λ|J 2 )ψ fc = 0. Because of the nondegeneracy of ω τ + λiJ we have from the factorization (12) (ω-λ k J)ψ k = 0, or, equivalently where ψ is the N by N matrix with columns ψ k and Λ = diag(Λ, 1? λ 2 , ...,Λ, N ). If ψ is invertible then defines the desired solution. It can be shown that ψ is indeed nondegenerate and that (14) actually defines the solution of (6) which, moreover, is real and orthogonal. But in Subsect. 1.4 we present another approach to the solution of (6) which is a bit more general. This proof will also provide the nondegeneracy of ψ and complete the above considerations. In this connection the concepts of symplectic geometry will turn out to be useful. We note that the solutions ω e O(N) so obtained have the property that the polynomials p{λ)= ±det(ω -λJ) and p( -λ) have no common roots. In other words, any two eigenvalues λ, λ' oϊωJ~ι satisfy λ + X φ0. We will denote the set of these matrices by E, i.e.
This is clearly an open subset of O(N) containing a neighborhood of the identity.
Since p(λ) does not vanish on the imaginary axis E decomposes into several components, depending on how many roots of p(λ) lie in the left half plane. With the aid of Theorem 1 it is easy to define a mapping φ in a unique way. We do this in the reduced form and rewrite the above factorization (12) in the form Then the image point M' = ψM is given by where the two factors were exchanged. This equation is readily verified from (10). Thus the determinants P(λ), P\λ) of (16), (17) respectively are identical. By Theorem 1 any splitting P(λ)=p(λ)p( -λ) gives rise to a unique factorization. Hence for (17) there exists a unique ω' and A'(λ) = ω f -λJ with This gives rise to a well defined mapping ω-+ω r taking E into itself. The uniqueness is achieved by the requirement (18)

Isospectral Deformations
From the above considerations we obtain the desired integrals. For this purpose we write the mapping in terms of an isospectral deformation. The Eq. (4) is already ΓΛΠ in this form but yields only k= -"trivial" integrals tr(M 2v ), v = l,2, ...,/c, (in fact, these are coadjoint invariants of O(N)) which is not sufficient for complete integrability. As was pointed out by Novikov [17] it is crucial to have such a representation for a matrix depending polynomially on a parameter λ.
Consequently the polynomials f k (M, λ) = tr(M + λJ 2 ) k are integrals of φ. In other words, the characteristic polynomial det(L(/l)μl) is preserved by ψ, or in homogeneous form the coefficients Q aβγ (M) for α^ l, 2oc-\-β + y = N provide k 2 integrals if N = 2k, or k(k + l) integrals if N = 2k +1. These integrals f k (M,λ) or Q aβγ (M) are precisely the same as for the Euler-Arnold Eq. (5). Indeed, for these equations the Lax representation was found by Manakov [22] in the form showing that also f k (M, λ), or Q aβγ (M) are integrals of the motion.
It is well known that these functions are in involution with respect to the Poisson structure (11) and independent, making the system (5) completely integrable. Since our discrete map ψ:M-^M f preserves the same Poisson structure (11) as well as the functions f k {M,λ) we conclude that ψ is also integrable. We summarize these results in Theorem 2. The discrete Euler Eq. (4) is equivalent to the isospectral deformation where L k = M k + λJ 2 , A k = ω k -λJ, M k+ί =xp{M k ). This mapping ψ preserves the Poisson structure (11) and is completely integrable. It preserves the same Poisson structure and integrals F t as the continuous system (5) for the motion of the rigid body.
The "integration" of this system is now rather straightforward, since the integration of the continuous case is known: The nonsingular compact level sets T c = f]{F~ c t ) consist of a finite union of tori, according to well known arguments i [20]. Since our mapping ψ preserves the Poisson structure (11) as well as the functions F t it commutes with all commuting Hamiltonian flows generated by the F i9 defined by M=[M, VFJ. On each such torus our mapping ψ must be a translation with respect to the affine structure, determined by these flowί In this case this mapping can be represented as a shift along the trajectory of a certain integral H (see [1] and Subsect. 1.5).
We will show that in our case T c is the real part of a complex Abelian variety of the curve det(M + λJ 2 -μl) = 0, and Eq. (4) determines a translation on it. In fact, it turns out to be the same Prym variety occurring in the integration of the Euler-Arnold equation (see, for example, [21] and [33]).

The Symplectic Geometry of Eq. (6)
First of all we write (6) as ω~ίJ -Jω = M, ωω τ = I, and introducing W=ω~ιJ we obtain the quadratic matrix equation with the additional condition W T W=J 2 . If v is an eigenvalue of W then Comparing with (13) we see that Q(v) = v 2N P{v " x ). Since 0 φ Σ the splitting Σ = Σ + KJΣ_ defines the splitting S = S + uS_ of the set S of roots (21): We require that this splitting satisfies the following conditions: Such splitting exists if (21) has no purely imaginary roots. Notice that now we do allow multiple roots but we do suppose that no root belongs to both components S+ and S_. In particular, purely imaginary roots are excluded. We formulate Theorem 1 in the equivalent form:

Theorem Y. For any splitting S = S+vS-with the properties (22) there exists a unique solution W of (20) (and therefore the solution of (6) ω = JW~ι) with specW=S + .
For the proof the solution of (20) will be played back to a problem of symplectic geometry, namely the determination of invariant subspaces of a linear Hamiltonian vector field.
The real 2N x 2N matrix in question is We look for an iV-dimensional invariant subspace of A: i.e. if the invariant subspace can be viewed as the graph y = YX 1 x,z=l j then we obtain for W= YX~\ Eq. (20).
To prove that W T W=J 2 , we note that A is antisymmetric with respect to the symplectic bilinear form (24) and Az can be viewed as the Hamiltonian vector field with Hamiltonian Jf = !(H 2 ,z) = i(-|Jx| 2 + |#), #=(~0 J2 j), since BA = H and Since ( vI we have and the spectrum of A is S + uS_. Denote the iV-dimensional eigenspaces of A with respect to S+, S_ by F+, F_, respectively. Because of S + =S+, S_ =S_ they are real and since μ t -f μ/φO for μ /? μ ; eS + they are Lagrangian, isotropic spaces with respect to the symplectic form [ , ] and the symmetric form J f respectively, as follows from the following lemma. Note that V+, V-are real subspaces since S + =S + , <S_=5_, while E μ are generally complex. Now we return to the proof of Theorem 2. Let z l5 ...,z n be any basis in V + ; combining these as column vectors of an JV x 2N matrix of rank N we have from AV + CV + for some real NxN matrix C+. To prove (23) we use the relations (26), yielding for any w, ι;eR N , Assuming v is so chosen that X+v = 0, we set u = C+υ, so that X + u = X + Cv= Y + v and we find from the above identity i.e. Y + v = 0, hence Z+v = 0 9 i.e. £> = 0. Therefore detX+ ΦO, proving (23). Thus F + is given by y = W + x. Since V + lies on the zero energy surface it follows \Jx\ 2 -\W + x\ 2 = 0 for all xeWL N proving W+W + =J 2 , hence ω = JW+ 1 is orthogonal. Moreover spec VF+ =iS+ proving Theorem Γ.

y.5. T/zβ Integration of the Discrete Euler Equation
Now we apply our results to finding the solution of (4), following the procedure which was described in the continuous case by Dubrovin [9,10]. For the initial data X θ9 X x e O(N) we define ω x =X\X Q = Xϊ ι X 0 e O(N) and M x = ω[ J -Jω v As follows from the previous considerations Eqs. (4) define only a correspondence, but if we fix the splitting S = S + uS_ of the roots of the polynomial β(v) = (v 2 / -μM k -J 2 ), which in fact does not depend on k, then we have a well defined mapping f s+tS -: ( ω k> M k )-^(ω k+x , M k+ J. In order to "integrate" this dynamics consider the spectral curve Γ: We will assume that J 2 has distinct eigenvalues different from zero: Jf φ J 2 for i Φ j, is meromorphic on Γ whose poles define a divisor @ = @ 1 + ...+@ g+N _ ί (see [10,16]). In the points at infinity P t eΓ, where μ&λJf, Λ->oo (i = 1,...,N) we have V>\Pj) = #j. (30) This means that ψ\λ, μ) is the basis of the linear space of meromorphic functions with the pole divisor ^ 2, determined by the conditions (30). In our case M is skew symmetric, therefore Γ has a symmetry σ: Γ-+Γ, σ 2 = id: The divisor Q) also is not arbitrary because of the following proposition. Denote ψ τ (λ, -μ) by tp*(/ί, μ) and fix λ e C such that the eigenvalues μ l5 ..., μ N of M + ΛJ 2 determined by (27) are distinct.

where B is the set of branch points of μ as a function of λ, and « means linear equivalence of divisors.
This equivalence is given by the function F(λ, μ) = ψ*{λ, μ)ψ(λ, μ) as follows from the proposition. Thus 2 belongs to the shifted Prym variety PC J{Γ). We restrict ourselves to these considerations because the detailed discussion of the algebraicgeometric aspects of this spectral problem can be found in the literature (see [21] and references therein). The corresponding problems of real algebraic geometry are considered in [23]. Now we use the representation (19) for describing the analytic properties on Γ of ψ k for arbitrary L Fix some splitting Σ = Σ + KJΣ_. As follows from (19) ψ(λ, μ) = (ω k -λJ)ψ k (λ, μ) is an eigenvalue of M k+ί + λJ 2 \ (M h +! + λ J 2 ) ψ(λ, μ) = μφ(λ 9 μ).
This means that we can define ψ k+ί as (36) Notice that ψ k+1 does not satisfy Dubrovin's normalization (29) which is required only for ψ v One can see from (36) that ψ k+x has N new poles at the "infinities" P l9 ..., P N . In order to find the new zeros consider the hyperbola Jf, determined by the equation and the intersection JtfnΓ. This intersection is described by the equations μ = λ~x and 1 2 2 which coincides with (11). So we have 2N points of intersection which we denote QX,..., ξ)χ, Qΐ,..., Qΰ in agreement with the splitting Σ = Z + uΓ~.As follows from the construction of ω k (see Sect. 3) n This means that QX, . ,6JV are the new zeros of ψ k+ί . Thus we have proven

Lemma. For a given splitting Σ = Σ + ul_
the vector eigenfunction ψ k+1 (36) of the matrix M k+ 1 +λJ 2 has on the spectral curve Γ the following analytical properties, which determine ψ k+1 uniquely: 1. ψ k+1 has a simple pole in Sf depending on the initial data M 1 and the poles at the "infinities" P u ...,P N with asymptotics in 2. ψ k+ί has a zero of order k in Q^, ...,6^.
For ψ k+ί one can write the explicit formulas in terms of Prym's ^-functions as it was done, for example, by Bobenko in [24]. But here we restrict ourselves to the example of 0(3) (see below). We summarize the result of this section in the following theorem.  (4) and (2) can be expressed as some abelian function on P in the points z k = ί. 6

. Explicit Formulas for the Discrete Dynamics of the 3-Dimensional Rigid Body
We consider here Eqs. gives the equivalence of Γ and C/Zτ 1 +Zτ 2 , where τ x = 2 J X3 dx τ 2 -2\ τ x is real, τ 2 is purely imaginary.

χ2 γR(x)
The equation μλ = 1 in the variables x, w has the form it determines on Γ the set of six points, which we denote as Qϊ, Q 2 , Q 3 , Qϊ, Q 2 , Qΐ according to the splitting Σ = Σ + vΣ_ (see Fig. 1) This figure depicts the situation, corresponding to J\<H/M 2 <J\<3\ and sufficiently small M 2 and H, when all roots of (43) and P(λ) are real numbers. The condition that P(λ) has no purely imaginary roots, which is sufficient for the   This means that Qf, Q}, Qi tend to P l9 P 2 , P 3 when e->0 and the shift U = ΣQ i 0. One should observe that the continuous limit corresponds to the special splitting Σ 9 where: Σ + contains all roots in the right half plane of P(λ) of (13) Sect. 1.2, which becomes here

(λ 2 J 2 -l)(λ 2 J 2 -l)(λ 2 J 2 -l) + Hλ*-
As was shown before this continuous limit coincides with the classical problem about the force free motion of a rigid body, the explicit solution for which in terms of elliptic functions were found by Jacobi [25]. The comparison of this formula and the continuous limit of (46) may be complicated.

The Discrete Dynamics on Stiefel Manifolds and the Heisenberg Chain with Classical Spins
In this section we generalize the previous results and consider the functional S=£tr(X fc JX fc+1 ) on the sequences X = (X k ) of n + N matrices X k , l^n^N k satisfying the condition For n = ί our results agree with the solution of this problem proposed in [1]. Notice that for n = 1, N = 3 we have the problem about the stationary states of a Heisenberg chain with classical spins (see [1,[5][6][7] and Introduction).

The Equation of the Dynamics and Isospectral Deformations
The variation of S under the condition (1) leads to the equation where Λ k = Λ k is an (n x n) matrix multiplier. We assume J to be symmetric and nondegenerate.
so that The proof follows by a direct calculation. The determinant of L k has the form where we used the Weinstein-Aronszajn formula, see [10]. Therefore P(λ) is an even polynomial in λ of degree 2n. The factorization (4) If such a splitting is fixed and all zeros of P(λ) are distinct then for given X k _ ί ,X k the matrix X k+1 can be defined as follows.
Leaving the further discussion of the general case 1 ^ n < N we now turn to a detailed treatment of the case n = 1 in the next section. In this connection we want to mention the work of Adams, Harnad, and Previato [26], in which the matrices introduced in [10] are generalized.

Discrete Version of the Neumann System and the Heisenberg Chain With Classical Spins
For n = 1 we have the functional feeZ where x = (x k ), x k belong to the unit sphere S^" 1 in IR^: |x fc | = 1 and J = diag(J 1? ..., J n ). For N = 3 we have unit vectors in R 3 which can be interpreted as classical spins. In this case the functional S (10) defines the energy of a spin chain in the Heisenberg model (see [1,[6][7][8]). The equations of the stationary configuration have the form The multipliers λ k (1 x 1 matrix) can be found from therefore A k μ k |J-1 x k | 2 -2(J-1 x to x fc . 1 )) = 0.
In the literature (see, for example, [5,6]) usually the second possibility is considered because only in this case it is possible to have the continuous limit. This can be seen as follows. Let J = I + ε 2 J c and x k = x(t 0 + kε) for small ε, then from (11) we have for x(t) the equations or x" + 2J c xπμx (14) with μ = (λ -2)ε~2, i.e. the Neumann system describing the motion on the unit sphere under the influence of the potential U(x) = (J c x, x) [8,9,27]. For this reason we call the system (11), (13) a discrete version of the Neumann system. Now we explain how this system can be integrated with the help of Theorem 4. The matrix L k (3) has the form where x = x k -ί9 y = x k , XΛy = x®y-y®x. This is a special case of the matrix introduced and investigated in connection with the classical integrable systems in [10]. The determinant L has the form and the zeros λ= ±(x,J~ιy)~1. So we have two possible splittings of this set Σ = Σ + KJΣ_ :Σ + =(x,J~ίy)~i or Σ + = -(x,J~ίy)~1. It is easy to see that they correspond to two solutions of (12). This means that the procedure proposed in the previous section actually determines the dynamics.

N-ί
where v = A fj (μ-ju f ). We see that the mapping (11) is the shift on the Jacobi i=ί variety of the hyperelliptic curve (23). To find this shift we follow the same procedure as in Sect. 1.5; namely the eigenfunction ψ k+ί of the matrix L k+1 can be obtained as ψ k+ί = A k (λ)ψ k or

Ψk + i={J-λXk®Xk-i)Ψk-
We see that ψ k + 1 has a new pole at "infinity" P^ corresponding to μ= oo for (23) and a new zero in one of the points P+ : μ = 0, λ = + (x, J ~ * y) ~1, determined by the splitting Γ = Σ + uI'_. So the shift is l/ = P 00 -P + . Following the procedure proposed in the previous paragraph one can write the explicit formulas for ψ k and x k in terms of hyperelliptic β-functions. Such formulas were found first in [7] (see also [1]) by another approach, requiring some unmotivated steps. These results are in a good agreement with those of [1,7], as one can verify.

The Billiard Inside an Ellipsoid
Now we apply the above procedure to describe the billiard motion inside a domain ΩclR 1 * bounded by the ellipsoid Q = dΩ given by the equation where A is a positive symmetric N by N matrix. This application has some surprising features. In particular, it leads in a natural way to a mapping φ which is not directly related to the system under consideration while φ o φ = φ 2 is. In other words this mapping can be viewed as a "square root" of the billiard mapping. Moreover, this mapping φ commutes with the billiard mapping, hence takes a billiard orbit into another billiard orbit. Such a symmetry of the billiard problem was found in [1], We begin with the construction of this mapping φ.
The dynamics of the billiard in the domain Ω can be described as the stationary points of the functional S: The equation of motion can be written in the form (see Fig. 2): where y k = (x k -x k -ι)/\x k -x k -J, \y k \ = 1 be the momentum, the multipliers μ k , v k are determined from the conditions \y k \ = 1, (Ax h x k ) = l: with C = A ~1 /2 (compare with [1]). One can easily check that if (x k5 y k ) is a solution of (3), so is (x fc , yj. Moreover, x^r = Cy k + 1 = -x k + u y' k ' = -C " ^ =y k +1? showing that the dynamics of φ "contains" the billiard dynamics. In the continuous limit this symmetry means only that if x(s) is a geodesic on the ellipsoid (1) parametrized by the length then the trajectory of the vector x(s) -Cx(s) is also a geodesic line on the same ellipsoid (but in general s is not the length of x(s)). This geometrical fact was observed in [28]. Notice that the operator C transforms the unit sphere into the ellipsoid Q, whose equation can be rewritten as
The explicit formulas for the ellipsoidal billiard were found in [1] by using the connection with the spectral theory of the difference operator. Now we are also able to incorporate these formulas into the framework of Sect. 2.

Connection Between the Ellipsoidal Billiard and the Discrete Neumann System
Replacing x' k , y' k by x k+1 , y k+1 we rewrite the system φ of (4) as or as Xfc+i + Xfc-i^fcC" 1 **. (12) Denoting C~1x k = q k , \q k \ = l we have ft + i+ft-i^CΓ 1^ \q k \ = \ (13) which coincides with (11)  This theorem, reminiscent of the connection between the geodesic flow and the Neumann system of [9], can, however, not be considered as a "discretization" of Knόrrer's result [16]. First of all instead of the Gaussian mapping we have a linear mapping C. The second new feature is the appearance of the mapping φ. Moreover, this connection is unexpected because of the different character of the transition to the continuous limit for this discrete system: for the billiard it corresponds simply to the dynamics near the diagonal in Q x Q, but for the system (11), Sect. 2.2, the diagonal is not invariant, so we have to use another limit process. There is another possibility in considering the dynamics near the zero level of the integral F = (J~1x k ,x k - 1 ) which corresponds to λ k &0, x k+ί -3 rx k^i ι x0 [see formulas (5), (6), Sect. 2]. This possibility leads to the ellipsoidal billiard, as follows from Theorem 6.

Remark.
A natural question arises whether or not the solutions for our discrete systems have the form x k = x(kA), where x(i) is a solution of the corresponding classical problem. The answer is: in general not. This can be seen in the example of the ellipsoidal billiard. Suppose that every orbit lies on one geodesic. Fix q 0 on a symmetry plane Oxy: q Ό = (x θ9 y 0 ,0), x 0 : y 0 φ 0, and choose q _ t on the intersection of the ellipsoid with the plane, orthogonal to Oxy and containing the normal to the ellipsoid in q 0 . Obviously q λ = σ(q. ι \ σ(x,y,z) = (x,y, -z). For q 0 and q^x being close there exists a unique geodesic line y, passing through q 0 and g_ v From our assumption it follows that it passes also through q v From the uniqueness σ(y) = y and therefore the tangent to y in q 0 is orthogonal to Oxy. Thus y passes through q 0 in a definite direction and therefore γ does not depend on g _ v This means that γ coincides with our planar curve and therefore cannot be a geodesic.