Symmetry, Integrability and Geometry: Methods and Applications A Discretization of the Nonholonomic Chaplygin Sphere Problem ⋆

The celebrated problem of a non-homogeneous sphere rolling over a horizontal plane was proved to be integrable and was reduced to quadratures by Chaplygin. Applying the formalism of variational integrators (discrete Lagrangian systems) with nonholonomic constraints and introducing suitable discrete constraints, we construct a discretization of the n-dimensional generalization of the Chaplygin sphere problem, which preserves the same first integrals as the continuous model, except the energy. We then study the discretization of the classical 3-dimensional problem for a class of special initial conditions, when an analog of the energy integral does exist and the corresponding map is given by an addition law on elliptic curves. The existence of the invariant measure in this case is also discussed.


Introduction
The Chaplygin problem on a non-homogeneous sphere rolling over a horizontal plane without slipping is probably one of the best known integrable systems of the classical nonholonomic mechanics. Although being non-Hamiltonian in the whole phase space (see [2]), the equations of motion possess an invariant measure, which is a rather strong property putting them rather close to Hamiltonian systems.
It is even more interesting that the Chaplygin sphere appears to be an algebraic integrable system in the sense that generic level varieties of the first integrals are open subsets of 2-dimensional Abelian varieties and, after an appropriate time reparameterization, the phase flow becomes a straight line uniform on them, [6,8].
Note that a Lax pair with a spectral parameter for the Chaplygin sphere, which could provide all the constants of motion, is still unknown and, probably does not exist. Hence, one cannot use the powerful method of Baker-Akhieser functions to find theta-function solution of the problem or to construct its integrable discretization by applying a Bäcklund transformation, as described in [12,7], or by factorization of matrix polynomials.
Contents of the paper. We briefly recall the equations of motion of the Chaplygin sphere, as well as their n-dimensional generalization. Then we construct a discretization of the problem by applying the formalism of variational integrators with nonholonomic constraints recently developed in [4,14]. Apart from the nonholonomic distribution on the tangent bundle T Q of the configuration manifold Q, the formalism requires introducing discrete constraints on Q × Q which, in a certain sense, must be consistent with the continuous ones. As one can see in the literature [4,10,15], the choice of discrete constraints is crucial for the behavior of the discretized nonholonomic system: the latter may inherit the main properties of its continuous counterpart, or may not. As an example, we note that although continuous systems with stationary nonholonomic constraints possess the energy integral, almost all their known discretizations do not enjoy this property (see e.g., [5]). Nevertheless, for a class of discretizations considered in [10] there exists a natural method to specify discrete constraints which ensures the exact preservation of energy.
Our choice of discrete nonholonomic constraints on E(n) × E(n) that mimics the condition of the sphere rolling without slipping over a horizontal plane allows us to construct a map discretizating the Chaplygin sphere problem, which has the form of a momentum conservation law and therefore preserves all the momenta integrals of the original system.
We then consider the discretization of the classical 3-dimensional problem for the case when the discrete angular momentum is vertical. In this special case the structure of the map is reminiscent to that of the Veselov-Moser discretization of the Euler top on SO(3) [19,16] and an analog of the quadratic energy integral does exist.
This implies that the invariant manifolds of the discretization map are elliptic curves E and the map is described as an addition law. However, in contrast to most known integrable discrete systems, in the discrete Chaplygin sphere the translation on E depends not only on the constants of motion but also on the point on the curve. Thus, in order to find the explicit solution, we arrive at a rather difficult problem of reparameterization of E or its real part, which would make the translation constant. We notice that this problem is equivalent to the problem of the existence of an invariant measure of the map.
2 The Chaplygin sphere and its multidimensional generalization Following Chaplygin [3], consider a dynamically non-symmetric sphere with inertia tensor J, radius ρ, and mass m rolling without slipping over a horizontal plane. Assume that the mass center and the geometric center C of the sphere coincide. Let γ = (γ 1 , γ 2 , γ 3 ) T be the vertical unit vector and ω = (ω 1 , ω 2 , ω 3 ) T , v ∈ R 3 be respectively the angular velocity of the sphere and the velocity of its center in the moving frame. The condition of non-slipping of the point P of contact of the sphere with the horizontal plane is v + ρω × γ = 0. (2.1) Here and below, × denotes the vector product in R 3 . On the configuration space of the problem, the Lie group E(3), these equations define nonholonomic constraints: for any two positions of the sphere on the plane there exists a linking trajectory that satisfies (2.1). Under these constraints the equations of motion can be reduced to the following closed system for ω, γ [3,9] I being the identity 3 × 3 matrix. Let K = Λω − D ω, γ γ be the vector of the angular momentum of the sphere with respect to the contact point P . Then the system (2.2) also admits representation in forṁ Hence, like γ, the momentum K is fixed in space. As a result, the system possesses four independent first integrals Here the last integral is the kinetic energy of the sphere and, since it can be rewritten in terms of K, γ: In addition, the system (2.2) possesses an invariant measure hence, by the Euler-Jacobi theorem (which is also often refereed to as the Jacobi last multiplier theorem, see e.g., [20]), it is integrable by quadratures and its generic invariant varieties are 2-dimensional tori. There are two special types of the initial conditions, when the equations of motion are simplified. In the first case the momentum K is horizontal, K, γ = 0. As shown by Chaplygin, after a time reparameterization and introducing spheroconic coordinates on S 2 = γ, γ the variables separate and the system reduces to hyperelliptic quadratures. Theta function solution for the unreduced Chaplygin system in this, as well as in the generic case, was obtained in [8].
The second special case is described below.
The case of the vertical momentum K. As noticed in [3] and as follows from the energy integral in (2.4), under the special initial conditions K = hγ, h = const one has K, ω = h γ, ω = l, i.e., As a result, the first vector equation in (2.2) takes the form of the Euler top equations and integrals (2.4) reduce to ω, Λω = l, Λω, Λω = k 2 , l, k = const.
Hence, for almost all initial conditions K = hγ, the variables ω i , γ i are elliptic functions of the original time t and the solution of the reduced system is periodic. Multidimensional generalization. Now, following [9], consider the generalized Chaplygin problem on an n-dimensional ball rolling without slipping on an (n − 1)-dimensional 'horizontal' hyperplane H in R n . The configuration space for the ball is the Lie group E(n) = (R, r), where R ∈ SO(n), r ∈ R n are the rotation matrix of the sphere and the position vector of its center C. For a trajectory R(t), r(t), define the Lie algebra element (ω, v C ), where are respectively the angular velocity and the velocity of C in the frame attached to the sphere. Let now γ ∈ R n be the unit vector orthogonal to the hyperplane H and directed 'upwards', i.e., from H to the center C, and, as above, ρ be radius of the ball. Then the condition for the sphere rolling without slipping that generalizes (2.1) reads v C + ρ ω γ = 0. (2.9) One can show that this vector constraint determines a non-integrable distribution on the tangent bundle T E(n), which is neither left-no right-invariant with respect to the action of E(n).
Next, introduce the angular momentum of the sphere in the body frame with respect to the center C, J being a constant diagonal n × n matrix. Then, as shown in [9], the motion of the sphere is described by the following Euler-Lagrange equationṡ where F ∈ R n is the reaction force acting at the point P of contact of the sphere with the hyperplane H. Here all the vectors are considered in the frame attached to the ball. Now, differentiating the constraints (2.9) and using the second equation in (2.11), we find F = mρωγ. Then (2.11) giveṡ where D = mρ 2 , Γ = γ ⊗ γ ≡ γγ T . Then this system can be represented in the following compact commutative form that generalizes (2.3) can be regarded as the angular momentum of the ball relative to the contact point P .
It is natural to introduce the nondegenerate inertia operator A : so(n) → so(n) * such that K = Aω. Since A is nondegenerate, equations (2.12) give a closed system for the variables ω, γ.
As follows from the form of (2.12), Ad R K is a constant tensor in the space. Hence this system has the following set of integrals (2.14) naturally, the system also possesses the energy integral and, as shown in [9], the invariant measure µ dω ∧ dγ with density µ = √ det A, which depends on the components of γ only.
In the classical case n = 3, under the isomorphism between so(3) and R 3 , the integrals (2.14), (2.15) and the measure transform to (2.4) and (2.6) respectively. Various properties of the Chaplygin sphere and its n-dimensional generalization were studied in [6,11,8,18], although the integrability or nonintegrability of this generalization was not proved.

Discrete mechanical systems with nonholonomic constraints
Before describing the discrete setting, let us recall that a continuous nonholonomic Lagrangian system is a triple (Q, L, D), where Q is a smooth n-dimensional configuration space, L : T Q → R is a smooth function called the Lagrangian, and D ⊂ T Q is a k-dimensional constraint distribution. Let q = (q 1 , . . . , q n ) be local coordinates on Q. In the induced coordinates (q,q) on the tangent bundle T Q we write L(q,q). It is assumed that the map The equations of motion are given by the Lagrange-d'Alembert principle: where δq(t) ∈ D q(t) for t ∈ (a, b) and δq(a) = δq(b) = 0. This principle is supplemented by the condition that the curve itself satisfies the constraints. Note that we take the variation before imposing the constraints.
Assuming that the constraint distribution is specified by a set of differential forms A j (q), j = 1, . . . , s < n, Coupled with (3.2), they give a complete description of the dynamics of the system. Note that equations (3.3) conserve the energy The discrete Lagrange-d'Alembert principle and equations. Let, as above, Q be a smooth manifold. According to [4], a discrete nonholonomic mechanical system on Q is defined by three ingredients (2) an (n − s)-dimensional distribution D on T Q (given by equations (3.2)); (3) a discrete constraint manifold D d ⊂ Q × Q, which has the same dimension as D and satisfies the condition (q, q) ∈ D d for all q ∈ Q.
The dynamics is given by the following discrete Lagrange-d'Alembert principle (see [4]), Here D 1 L and D 2 L denote the partial derivatives of the discrete Lagrangian with respect to the first and the second inputs, respectively. The discrete constraint manifold is specified by the discrete constraint functions which impose the above restriction on the solution sequence {(q k , q k+1 )}.
Remark 1. If the discrete Lagrangian L is obtained from a continuous one, L(q,q), via a discretization mapping Ψ : Q × Q → T Q defined in a neighborhood of the diagonal of Q × Q, i.e., L = L • Ψ, then the variety D d must be consistent with the continuous distribution D. Namely, D d is locally defined by the equations A j • Ψ = 0, j = 1, . . . , s.
We emphasize that, in general, the discretization mapping is not unique and hence there are many ways to define the discrete Lagrangian L and the discrete constraint manifold D d for a given nonholonomic system (Q, L, D). 2 The discrete principle implies that the solutions of a discrete nonholonomic system are represented by sequences {(q k , q k+1 )} that satisfy the discrete Lagrange-d'Alembert equations with multipliers The multipliers λ k j are determined from the discrete constraints (3.4), but, in general, not uniquely. Hence, the map Q × Q → Q × Q defined by (3.4), (3.5) is multi-valued.
Remark 2. According to [4], equations (3.5) introduce a well-defined mapping (q is invertible for each (q k , q k+1 ) in a neighborhood of the diagonal of Q × Q.

A discretization of the Chaplygin sphere
Now we apply the above approach to discretize the generalized Chaplygin sphere problem on E(n). The trajectory of such system is a sequence (R k , r k ), k ∈ Z. We choose the discrete Lagrangian in form Here the rotational part 1 2 Tr(R k JR T k+1 ) coincides with the Lagrangian of the Euler top on SO(n) introduced in [19].
Note that, in view of (2.8), the continuous constraints (2.9) can be rewritten aṡ where γ is the unit normal vector in the fixed frame (without loss of generality we can set it to be (0, . . . , 0, 1) T .) Then, in view of (4.1), (4.2), in our case the discrete Lagrange-d'Alembert equations (3.5) take the form where Λ k is the symmetric matrix Lagrange multiplier and f k = (f k 1 , . . . , f k n ) T , is the vector multiplier corresponding to the constraints (4.2).
In view of the property R k R T k = I, equations (4.3) yield Following [19], we introduce the discrete angular momentum of the sphere with respect to its center C in space Then (4.5), (4.4) give rise to Next, like in [19], introduce the discrete momentum of the sphere in the body frame where Ω k is the finite rotation in the body frame. In the continuous limit Ω k → I+εω, ω ∈ so(n), it transforms to the momentum (2.10).
Then the system (4.6) gives where represent the vectors f k , γ in the body frame. Equations (4.7) can be regarded as a discrete version of the equations of motion of the Chaplygin sphere (2.11). Discrete constraints. In order to determine the vector multiplier f k we must specify discrete constraints on E(n) × E(n). A naive choice that imitates the form of the continuous constraints (2.9) is ∆r k + ρΩ kγ = 0, whereΩ k = R k+1 R T k ∈ SO(n) is the finite rotation of the sphere in the space frame. However, this choice allows that the displacement ∆r k of the center C may be not "horizontal" (i.e., orthogonal to the vertical vectorγ), which is incompatible with the mechanical setting of the problem.
For this reason, our choice of the corresponding discrete constraints (3.4) will be which does ensure that ∆r k is orthogonal toγ. In the continuous limitΩ k → I + ω it transforms to the constraint (4.2). Then, in view of (4.9) and the second equation in (4.7), and, therefore, where, as in Section 2, Γ k = γ k ⊗ γ k ≡ γ k γ T k . As a result, we eliminate the multipliers F k and from (4.7), (4.8) get the equations which define an implicit map C : (Ω k−1 , γ k−1 ) → (Ω k , γ k ).
Proposition 1. The map C admits the following compact representation

12)
where K k is the discrete analog of the momentum with respect to the contact point P of the sphere, Indeed, in the continuous limit Ω k → I + εω, ε 1 the matrix K/ε tends precisely to the angular momentum (2.13) of the n-dimensional sphere with respect to its contact point, whereas the relations (4.12) transform to the continuous equations (2.12).
The latter equation, in general, has several solutions, hence the map C is multi-valued.
An an immediate corollary of Proposition 1, we obtain Proposition 2. Regardless to a branch of the map C, it preserves the set of momenta integrals (2.14) of the continuous Chaplygin sphere problem with K replaced by K.
This property gives a solid justification of our choice of the discrete constraint (4.9). Note however that a discrete analog of the energy integral (2.15) in the generic case is unknown (see also Proposition 3 below).
Proof of Proposition 1. The second equation in (4.12) follows directly from (4.11). Now replacing Γ k in the right hand side of (4.10) by Ω T k−1 Γ k−1 Ω k−1 and using the identity Ω k−1 Ω T k−1 = I, we get which, in view of (4.13), is equivalent to the first equation in (4.12).
Remark 3. The form of equations (4.12), (4.13) is reminiscent to that of the discrete Euler top on SO(n) first described in [19], however, the discrete momentum K in (4.13) depends not only on Ω, but also on γ. Note that, in view of Proposition 1, the subsequent momenta K k−1 , K k admit curious intertwining relations which are reminiscent to the following relation between the subsequent momenta of the discrete Euler top (4.14) (see [16])

Discrete Chaplygin sphere on E(3) and its particular solutions
In the case n = 3 we can make use of the following parameterization of the body finite rotations where q 0 , q i are the Euler parameters subject to relations q 2 0 + q 2 1 + q 2 2 + q 2 3 = 1 and q 0 > 0. The operator (5.1) describes a finite rotation in R 3 about the vector q = (q 1 , q 2 , q 3 ) T by the angle θ such that q 0 = cos(θ/2) and |q| = sin(θ/2) (see, e.g., [20]).
Let now M k , K k ∈ R 3 be the vector representations of M k , K k ∈ so * (3), Then, in view of (4.13), we get simple expressions In the continuous limit, when the angle θ is small, θ = εω, ε 1, we have q = ε 2 ω + O(ε 3 ), q 0 = 1 − O(ε 2 ), and, up to terms of order ε, which coincides with the expression for the angular momentum vector K in the classical Chaplygin problem. As a result, in the vector form the map (4.12) reads where Ω k is recovered from K k , γ k by solving (5.2) with respect to q and substituting the solution into (5.1). One can show that for real K, γ, the equations (5.2) have at most 4 and at least 2 real solutions q. Like in most of other discrete systems, in order to choose one of the branches of the map (4.12) one should restrict to the case of sufficiently small finite rotations Ω k . In this case only one of the above real solutions q will be small and should be taken as the appropriate branch. Like the continuous system (2.3), apart from the geometric condition γ, γ = 1, the discrete system preserves two independent integrals However, as simple numerical tests show, the energy integral (2.5) is not preserved. The special case K γ. Like the continuous Chaplygin system, the map (5.3) has the special case when K k = hγ k , h =const. In view of (5.2), this implies the following relation between γ and q 2     Then one obtains the reduced map G h : S 2 → S 2 , G h (γ k−1 ) = γ k evaluated as follows: 1) given γ k−1 , one recovers q as a solution of (5.4), 2) one calculates Ω k by (5.1) and γ k by the second equation in (5.3).
As a result, in our special case the complex invariant variety of the map G h is the spatial elliptic curve E, the intersection of the unit sphere γ, γ = 1 with the quadric given by (5.5).
The curve E is 4-fold unramified covering of the planar elliptic curve 3 )(z − l) and the points of E are parameterized by the Jacobi elliptic functions associated to E 0 . Assume for concreteness that Λ −1 1 < Λ −1 2 < l < Λ −1 3 . Then the parameterization reads (see e.g., [20,13]) γ 1 = C 1 cn(u|k), γ 2 = C 2 sn(u|k), γ 3 = C 3 dn(u|k), (5.7) where u is a complex phase parameter and k is the modulus of E 0 given by Therefore, for a fixed l, the map G h is reduced to one-dimensional map with the increment ∆u k is a function of u k , l to be determined below.
Remark 4. As seen from (5.6), the structure of the map G h is similar to that of the Veselov-Moser discretization (4.14) of the Euler top on so (3). As shown in [16,1], the latter discretization preserves the momentum and the energy integrals, whereas the corresponding increment ∆u k does not depend on the argument u k in the elliptic parameterization. The same holds for another exact discretization of the top found in [7].
In comparison with the Veselov-Moser discretization, our map G h involves the extra factor σ k , which is not constant on the orbits of G h . This gives an indication that under our map G h the increment ∆u k depends essentially on u k . Moreover, as G h is multi-valued, the function ∆u k (u k , l) is expected to be multi-valued as well. This stays in contrast with the solutions of the continuous equations (2.7) obtained from the Chaplygin sphere system under the condition K γ: the components of γ are elliptic functions whose argument changes uniformly with time t.
In this connection the following natural problem of the complete integrability of the map G h arises: is there a (multi-valued) reparameterization u = f (s, l), s ∈ C such that u k+1 = f (s k + ∆s) and the increment ∆s does not depend on s k ? Note that this property is equivalent to that the one-dimensional map (5.9) preserves the invariant measure µ(u)du, where the density µ(u) must satisfy the discrete Liouville equation and is related to f (s) by the formula df /ds = 1/µ. Due to the existence of the integral (5.5), this is also equivalent to that the map G h itself preserves an invariant measure on S 2 . Since u is the argument of an elliptic function, the restriction of the map (5.9) onto the real axis describes a diffeomorphism of a circle with the angular function ∆u k (u k , l). Hence, the answer to the above integrability question requires studying the properties of the diffeomorphism and applying known theorems in this field (see, e.g., [17] and references therein).
Determination of ∆u k (u k , l). For convenience, in the sequel we omit the discrete time index k and denote G h (γ) =γ, u + ∆u =ũ. To determine the function ∆u(u, l), we shall use the following property. Proposition 4. The map G h admits the following implicit form that involves the sum and the difference of γ,γ only: . (5.11) Proof . As follows from (5.6), and, therefore, Using (5.12), it is easy to check that and that the factor σ defined in (5.6) is the solution of the equation Then, in view of relation γ, Λ −1 (γ + γ) = 1 2 γ +γ, Λ −1 (γ +γ) , one obtains (5.10), (5.11).
Below we assume h = 1 in (5.11), which can always be made by appropriate rescaling of Λ and D.

Conclusive remarks
In given paper we constructed a discretization of the Chaplygin sphere problem on E(n) which preserves all the momenta integrals and, in the particular case of the motion, the energy. This case is reduced to the multi-valued map G h on so * (3) which is reminiscent to the Veselov-Moser discretization of the Euler top and is given by an addition law on elliptic curves. However, in contrast to known integrable maps of this type, the increment depends not only on the energy constant, but also on the point on the curve. According to [1], this implies that G h does not preserve the standard Lie-Poisson structure on so * (3), which is quite expected due to the nonholonomic origin of the map. The integrability of the map G h (which is interpreted as the possibility of writing its explicit discrete time solution) was shown to be equivalent to preservation of an invariant measure on the Lie algebra. Note that the latter property is often related to integrability in quadratures in the continuous case (recall the Jacobi last multiplier theorem).
In this connection it would be interesting to develop integrability criteria for discrete systems preserving an invariant measure, in particular to formulate a discrete version of the Jacobi last multiplier theorem.
Although the energy (2.5) of the Chaplygin sphere expressed in terms of the momentum K is not preserved in the discrete setting in general, this does not exclude the existence of another integral that transforms to the energy in the continuous limit.
Lastly, as shown in [18], the classical Chaplygin sphere problem is a natural example of a nonholonomic LR system on the direct product SO(3)×R 3 (the kinetic energy is left-invariant, while the constraint distribution is right-invariant with respect to the group action). One can show that the discrete Lagrangian (4.1) and the discrete constraint subvariety given by (4.9) possess the same properties, i.e., our discretization of the Chaplygin sphere is a discrete LR system on SO(n) × R n .
We then believe that it is worth studying properties of generic discrete LR systems and of their reductions, such as preservation of momenta, energy, and an invariant measure.