Three-dimensional discrete systems of Hirota-Kimura type and deformed Lie-Poisson algebras

Recently Hirota and Kimura presented a new discretization of the Euler top with several remarkable properties. In particular this discretization shares with the original continuous system the feature that it is an algebraically completely integrable bi-Hamiltonian system in three dimensions. The Hirota-Kimura discretization scheme turns out to be equivalent to an approach to numerical integration of quadratic vector fields that was introduced by Kahan, who applied it to the two-dimensional Lotka-Volterra system. The Euler top is naturally written in terms of the $\mathfrak{so}(3)$ Lie-Poisson algebra. Here we consider algebraically integrable systems that are associated with pairs of Lie-Poisson algebras in three dimensions, as presented by G\"umral and Nutku, and construct birational maps that discretize them according to the scheme of Kahan and Hirota-Kimura. We show that the maps thus obtained are also bi-Hamiltonian, with pairs of compatible Poisson brackets that are one-parameter deformations of the original Lie-Poisson algebras, and hence they are completely integrable. For comparison, we also present analogous discretizations for three bi-Hamiltonian systems that have a transcendental invariant, and finally we analyze all of the maps obtained from the viewpoint of Halburd's Diophantine integrability criterion.


Introduction
The problem of numerical integration, namely that of approximating the flow of a smooth vector field by an iterative scheme given in terms of a difference equation or a map, is one of the central problems of numerical analysis. If the underlying differential equation is Hamiltonian, or volume-preserving, or has some other important geometrical feature (such as being invariant under the action of a Lie group of symmetries), then as far as possible one would like to select a discretization scheme which preserves this feature, and this has led to the development of geometrical integration methods [4]. For the special case of completely integrable systems, ideally one would like to obtain discretizations which are themselves completely integrable. The area of integrable discretization has been developed quite extensively, especially from the Hamiltonian viewpoint, and a comprehensive review of the field can be found in the monograph [39]. In this paper we are concerned with a novel approach to discretization, which was used by Hirota and Kimura to obtain new integrable discrete analogues of the Euler and Lagrange tops [15,22].
The discretization method studied in this paper seems to be introduced in the geometric integration literature by W. Kahan in the unpublished notes [20]. It is applicable to any system of ordinary differential equations for x : R → R n with a quadratic vector fielḋ where each component of Q : R n → R n is a quadratic form, while B ∈ Mat n×n (R) and c ∈ R n . Kahan's discretization reads as x where is the symmetric bilinear form corresponding to the quadratic form Q. Here and below we use the following notational convention which will allow us to omit a lot of indices: for a sequence x : Z → R we write x for x k and x for x k+1 . Eq. (1) is linear with respect to x and therefore defines a rational map x = f (x, ǫ). Clearly, this map approximates the time-(2ǫ)-shift along the solutions of the original differential system, so that x k ≈ x(2kǫ). (We have chosen a slightly unusual notation 2ǫ for the time step, in order to avoid appearance of powers of 2 in numerous formulae; a more standard choice would lead to changing ǫ → ǫ/2 everywhere.) Since Eq. (1) remains invariant under the interchange x ↔ x with the simultaneous sign inversion ǫ → −ǫ, one has the reversibility property f −1 (x, ǫ) = f (x, −ǫ). In particular, the map f is birational. Kahan applied this discretization scheme to the famous Lotka-Volterra system and showed that in this case it possesses a very remarkable non-spiralling property. Some further applications of this discretization have been explored in [21,36].
The next, even more intriguing, appearance of this discretization was in the two papers by R. Hirota and K. Kimura who (being apparently unaware of the work by Kahan) applied it to two famous integrable systems of classical mechanics, the Euler top and the Lagrange top [15,22]. Surprisingly, the Kahan-Hirota-Kimura discretization scheme produced integrable maps in both the Euler and the Lagrange cases of rigid body motion. Even more surprisingly, the mechanism which assures integrability in these two cases seems to be rather different from the majority of examples known in the area of integrable discretizations, and, more generally, integrable maps, cf. [39]. We shall use the term "Hirota-Kimura type discretization" for Kahan's discretization in the context of integrable systems.
In the recent paper [32] the Hirota-Kimura integrability mechanism has been further investigated and its application to the integrable (six-dimensional) Clebsch system has been considered. The integrability of the Hirota-Kimura type discretization of the Clebsch system has been established, in the sense of: i) existence, for every initial point, of a fourdimensional pencil of quadrics containing the orbit of this point; ii) existence of four functionally independent integrals of motion. Note that for the purposes of paper [32], integrability of a dynamical system is synonymous with the existence of a sufficient number of functionally independent conserved quantities, or integrals of motion, that is, functions constant along the orbits. Other aspects of the notion of integrability, such as Hamiltonian properties or explicit solutions, still require further investigation. However, it is known that algebraically completely integrable cases of geodesic flow on SO(4) are related to the intersection of four quadrics in P 6 [2]. The Hirota-Kimura method of discretization has been recently applied to the classical three-dimensional nonholonomic Suslov problem in [7].
The above examples of Hirota-Kimura type discretizations suggested the following Conjecture [32]. For any algebraically completely integrable system with a quadratic vector field, its Hirota-Kimura type discretization is algebraically completely integrable. Since algebraically completely integrable systems generically correspond to linear flows on abelian varieties [40], this statement should be related to addition theorems for multidimensional theta-functions.
The aim of this paper is both to study how this novel method of discretization applies to a set of algebraically integrable systems in three dimensions, and to see how these results compare with the analogous discretizations of some quadratic vector fields with transcendental invariants. The former set of systems considered are algebraically integrable in the sense that they have a sufficient number of algebraic integrals in involution; however, there are various other (more stringent) notions of algebraic complete integrability. Through this study we are able both to verify the above conjecture for a new set of examples, and to gain further understanding of how the integrability of the discretization depends on the algebraic nature (or otherwise) of the integrals of motion in the original continuous system. Kahan's discrete Lotka-Volterra system illustrates the subtlety of this dependence, as we now describe.
Kahan used his approach to discretize the Lotka-Volterra systeṁ which preserves the Poisson bracket {x, y} = xy, or equivalently, the symplectic form 1/(xy) dx ∧ dy. This is an integrable system with one degree of freedom,ẋ = {x, H}, y = {y, H} with Hamiltonian H = log xy − (x + y).
Kahan's discretization for this system reads [20] x This discretization preserves the same symplectic structure as the original system of ordinary differential equations, for which it provides a numerically stable integration scheme which appears to retain the qualitative features of the continuous orbits (which are closed curves H = constant in the positive quadrant x > 0, y > 0) [37]. In fact, as noted in [31], Kahan's discrete Lotka-Volterra system is algebraically integrable for ǫ = ±1. To be precise, when ǫ = 1, it reduces to the second order recurrence x n+1 x n−1 = x n (2 − x n ) for the x coordinate, which belongs to the class of antisymmetric QRT maps studied in [43]. This recurrence is linearizable (the iterates satisfy a linear recurrence of sixth order [19]), and the map has the integral which (for fixed H) defines a quartic curve of genus zero; this is also an integral for ǫ = −1 (which can be seen immediately from the reversibility property of Kahan's discretization scheme). However, for other non-zero values of ǫ, Kahan's discrete Lotka-Volterra system should not be algebraically integrable; we present some numerical evidence for this in section 7 below. Indeed, since the integral H for the original system is transcendental, from continuity arguments one would expect that (at least for small enough ǫ) any integral of the the discretization should be transcendental as well. Further numerical studies, as mentioned in [32], indicate that this discrete system may well be non-integrable, with characteristics of chaos only evident by zooming in deeply on regions of the phase plane. The outline of the paper is as follows. In the next section, we briefly review the Euler top together with the discrete Euler top found by Hirota and Kimura. In section 3 we describe six quadratic bi-Hamiltonian flows in three dimensions, which were presented in [11] (extending results in [3]), and are associated with pairs of real three-dimensional Lie algebras. Moreover, each of these systems, which we denote by E i for i = 1, . . . , 6, is algebraically integrable; the system E 6 is equivalent to a special case of the Euler top. The fourth section is devoted to applying the Hirota-Kimura discretization scheme to these six systems, to obtain discrete systems (or maps) in three dimensions which we denote by dE i , and in section 5 we present the explicit solutions of these maps for i = 1, . . . , 5 (the case i = 6 being already included in the work of Hirota and Kimura [15]). Section 6 is concerned with applying the same discretization method to three other bi-Hamiltonian systems from [11] which have transcendental integrals. In section 7 we present the results of applying Halburd's Diophantine integrability test to each of the maps obtained, and prove that all but one of them are Diophantine integrable in the sense of [12]. The final section is devoted to some conclusions.

Euler top and its Hirota-Kimura type discretization
The so(3) Euler top is a well-known three-dimensional bi-Hamiltonian system belonging to the realm of classical mechanics [34]. The differential equations of motion of the Euler top readẋ with α i being real parameters of the system. We recall that this system can be explicitly integrated in terms of elliptic functions, and admits two functionally independent integrals of motion. Indeed, a quadratic function H(x) = γ 1 x 2 + γ 2 y 2 + γ 3 z 2 is an integral for Eqs.
(3), if γ 1 α 1 + γ 2 α 2 + γ 2 α 2 = 0. In particular, the following three functions are integrals of motion: Clearly, only two of them are functionally independent because of α 1 H 1 +α 2 H 2 +α 3 H 3 = 0. x − x = ǫα 1 ( yz + y z), Thus, the map f : x → x obtained by solving (4) for x, is given by: Apart from the Lax representation which is still unknown, the discretization (5) exhibits all the usual features of an integrable map: an invariant volume form, a bi-Hamiltonian structure (that is, two compatible invariant Poisson structures), two functionally independent conserved quantities in involution, and solutions in terms of elliptic functions. For further details about the properties of this discretization we refer to [15] and [30].

Some bi-Hamiltonian flows related to real three-dimensional Lie algebras
Hamiltonian systems in three dimensions provide the simplest non-trivial examples of degenerate Poisson structures, where the rank of the Poisson tensor is less than the dimension of the phase space. In three dimensions, a non-trivial Poisson tensor P has rank two at generic points of the phase space, which means that (at least locally) there exists a Casimir function K and another function φ such that in local coordinates x 1 , x 2 , x 3 ; cf. Theorem 5 in [9]. This can be expressed in invariant form, by using the standard volume three-form Ω = dx 1 ∧ dx 2 ∧ dx 3 in R 3 to associate P with the one-form J = P Ω = φ dK. An important thing to observe from the form of the Poisson bracket (6) is that in three dimensions the Poisson tensor can be multiplied by an arbitrary function while preserving the Jacobi identity. Given an Hamiltonian systemẋ = {x, H} defined in terms of the bracket (6) with an Hamiltonian function H (functionally independent of K), it is clear that the equations of motion have two independent integrals, namely H and K. Moreover, by fixing the value of the Casimir function (which may not be defined everywhere), we can regard this locally as a system with one degree of freedom which is integrable on each of the two-dimensional symplectic leaves K = constant. However, for complete integrability the global existence of H and K is required. Gümral and Nutku made a detailed study of the geometry of three-dimensional Poisson structures, and considered the conditions for the existence of globally integrable bi-Hamiltonian structures [11]. For a given three-dimensional system to be bi-Hamiltonian it is necessary and sufficient that the Jacobian at an arbitrary point be a Poisson tensor and that there exist two globally defined and (almost everywhere) functionally independent integrals of motion. Associated with two independent integrals K and H, there are two compatible Poisson tensors, such that K is the Casimir for one Poisson structure while H is the Casimir for the other. In other words, if the dimension is three then two compatible Poisson tensors are completely determined by the constants of motion, and according to a relevant Theorem by Magri [24], provided certain technical conditions are satisfied, bi-Hamiltonian systems are completely integrable in the sense of Liouville-Arnold. Furthermore, in this setting there is an invariant volume form Ω (not necessarily canonical) which is preserved by the bi-Hamiltonian flow. An important example of such flows corresponds to Nambu mechanics [28], given bẏ which in these coordinates gives a divergenceless vector field (divẋ = 0); this means that the canonical measure is preserved by the flow. In particular, the Euler top is an example of Nambu mechanics in three dimensions; for other examples of Nambu mechanics in optics and elsewhere, see [17].
In [11] the authors present a list of all non-trivial bi-Hamiltonian flows that are associated with pairs of real three-dimensional Lie algebras and their Casimir invariants (as described in [29]); this list extends results in [3]. To be more precise, they consider pairs of real Lie-Poisson algebras defined by pairs of linear Poisson structures P, Q, and write down vector fieldsẋ = V satisfying where K is the Casimir for Q, while H is the Casimir for P (and minus signs are included in order to be consistent with the conventions of Gümral and Nutku). In [11] twelve such systems are presented, extending a list in [3], and each flow preserves a corresponding measure given in coordinates (x 1 , x 2 , x 3 ) = (x, y, z) by related to the standard volume form by the conformal factor (or multiplier) c. 1 To begin with, we shall be concerned with only six out of the twelve systems in Gümral and Nutku's list, namely the ones which have non-transcendental integrals of motion. 1 In fact, on page 5704 of [11] the authors state that the given systems are all "with multiplier unity", and denoting the multiplier by M they say "these equations have M = 1 [...] they are Nambu mechanics representatives", but as should be clear from Table 1 this is not the case: two of the systems given there have a non-constant multiplier. For systems with non-constant multiplier, it is remarked in [11] that they may be only locally (but not globally) equivalent to Nambu mechanics, by a suitable change of coordinates.
They read  Table  1. For instance, the flow E 1 , given in Eq. (7), admits the bi-Hamiltonian structure given by the compatible pair (P (1) 31 = {z, x} = 2y, with conformal factor c 1 = 1/x 2 . The quantities H 1 = y/x and K 1 = zx + y 2 , preserved by the flow, are respectively the Casimir functions of P (1) and Q (1) . This is equivalent to say that the following Lenard-Magri chain [24] is satisfied: where dH 1 and dK 1 denote the differentials of the functions H 1 and K 1 respectively. The same scheme holds for the flows E i with 2 ≤ i ≤ 6.

Hirota-Kimura type discretization of the flows E i
The goal of this section is to show that Hirota-Kimura type discretizations of the bi-Hamiltonian flows E i , 1 ≤ i ≤ 6, provide completely integrable discrete-time systems. The following result holds.

Theorem 1. The Hirota-Kimura type discretizations of the bilinear flows
given in Eqs. (7)(8)(9)(10)(11)(12), read where the matrices A i (x; ǫ) are given in Table 2. The quantities H i (ǫ), K i (ǫ), given in Table 2, are integrals of motion for the maps (13). Moreover the maps (13) preserve the volume form Note that for small ǫ the birational maps (13) approximate the time shift along the trajectories of the corresponding continuous equations of motion (7)(8)(9)(10)(11)(12). The same invariant volume form (14), which is independent of ǫ, is preserved by both the continuous and the discrete systems. Proof: We shall prove Theorem 1 for just one case, namely i = 5. The remaining cases can be proved by similar straightforward computations.
The Hirota-Kimura discretization of the flow E 5 , given by Eq. (11), reads explicitly as that is The fact that the quantities Table 2. Matrices A i (x; ǫ) and discrete integrals of motion are integrals of motion of the map (15) is proved by the following computation. Equation that is, using Eq. (15), which is an algebraic identity. A similar computation shows that equation K 5 (ǫ) = K 5 (ǫ) is identically satisfied. We now prove that the map (15) preserves the volume form which is equivalent to saying that First of all we note that differentiating Eq. (15) with respect to x, y, z one obtains the columns of the matrix equation Computing determinants lead to Now, by using the map (15), a straightforward computation shows that the relation holds identically.
In the construction of an invariant Poisson structure for the maps (13) we shall make use of results from [5] (Proposition 15 and Corollary 16 there), which we restate here. Suppose that f : M → M is a smooth mapping of an n-dimensional manifold M, with an invariant volume form Ω (that is, f * Ω = Ω). Define ω to be the dual n-vector field to Ω such that ω Ω = 1, where as usual the symbol denotes the contraction between multivector fields and forms. It follows that if I 1 , . . . , I n−2 are integrals of f with dI 1 ∧· · ·∧I n−2 = 0, then the bivector field σ = ω dI 1 · · · dI n−2 is an invariant Poisson structure for f . If J 1 , . . . , J n−2 is another set of independent integrals and τ = ω dJ 1 · · · dJ n−2 is the corresponding Poisson structure, then σ and τ are compatible, i.e. for any constants a, b, the bivector field aσ + bτ is again a Poisson structure.
In particular for n = 3, if a three-form Ω, given by Eq. (14) in our case, is invariant under a map f defined by (13), we can define the dual trivector field so that for any integral I of f the bivector field is an invariant Poisson structure for f , as well as any linear combination of such bivector fields. Explicitly, the Poisson brackets of coordinate functions are given by Note that the inverse volume density φ(x, y, z) can be multiplied by an arbitrary integral of f without violating the Poisson property.
For the maps (13), the invariant Poisson structures P (i) , c −1 i Q (i) can be computed according to the following formulae: with 1 ≤ i ≤ 6, the summation convention is assumed for the index l, and above we have used (x 1 , x 2 , x 3 ) to denote (x, y, z). (The reader should note that these indices 1, 2, 3 for the coordinates in R 3 should not be confused with the index n used to denote iterates of maps in subsequent sections.) This corresponds to taking φ = H i K i /c i above, and then rescaling by the inverse of the product of the integrals, 1/H i (ǫ)K i (ǫ), in each case. Thus the following statement holds. Tables 3 and 4. The conformal factors c i are the same as in the continuous case, given in Table 1.

Theorem 2. The maps (13) admit the compatible pair of invariant Poisson structures
Note that Eqs. (16)(17) provide one-parameter deformations of the Lie-Poisson tensors P (i) , Q (i) given in Table 1. This is equivalent to saying that Tables 3 and 4 provide deformations of the real three-dimensional Lie algebras A 3,3 , sl(2, R), so(3), e(1, 1), e(2). Finally we note that the integrable discrete-time system dE 6 is just a particular case of the Hirota-Kimura discretization of the so(3) Euler top [15], whose bi-Hamiltonian structure has been presented recently in [30].
As shown in [15,22], and recently in [30], the integrable discrete-time systems obtained through the Hirota-Kimura type discretization seem to admit a straightforward construction of their explicit solutions, at least for the case of three-dimensional maps. Here we provide the explicit solutions for the discrete-time integrable systems dE i with 1 ≤ i ≤ 5. The cases i = 4, 6 are each special cases of the so(3) Euler top, whose solutions, both continuous and discrete, are investigated in [15,30], so here we present the solution only for i = 4, since i = 6 is similar. For comparison, in Table 5 we give the explicit solutions for the continuous-time flows E i with 1 ≤ i ≤ 5. The parameters α, β, γ, θ, λ, µ, k appearing in the table can easily be expressed in terms of the initial conditions and/or the integrals of motion.
We now construct the explicit solutions to the discrete-time systems dE i with 1 ≤ i ≤ 5, thus providing the discrete counterpart of Table 5. Let us recall that we consider each of x, y, z as functions on ǫZ. To simplify the notation we set x = x n , y = y n , z = z n , so that x = x n+1 , y = y n+1 , z = z n+1 . For the sake of brevity, henceforth the discrete integrals of motion H i (ǫ) and K i (ǫ) in Table 2 will be denoted respectively by H i and K i , 1 ≤ i ≤ 6.
The following statement holds.
The system dE 3 reads: It has two integrals of motion, , in involution with respect to the pair (P (3) (ǫ), c −1 3 Q (3) (ǫ)), as given in Tables 3 and 4. The solution to the continuous-time flow E 3 , as in Table 5, suggests the following ansatz for the solution of the map (34-36): x n = ν cos θ cosh T n , y n = ν sin θ cosh T n , z n = λ tanh T n , with constant parameters λ, ν, θ. By substituting the ansatz (37) into the formulae for the integrals, we see that H 3 = tan θ, while is constant (for all T n ) if and only if ν 2 = λ 2 /(1 − ǫ 2 λ 2 ). Upon setting λ = ǫ −1 tanh δ, in terms of another parameter δ (with δ/ǫ = O(1) in the continuum limit ǫ → 0) this gives ν 2 = ǫ −2 sinh 2 δ and K 3 = sinh 2 δ/(2ǫ 2 ). Substituting the ansatz into the third part of the map, namely (36), and using the addition formulae for hyperbolic functions, one can see that this equation implies that hence T n+1 − T n = 2δ. This implies that T n = 2δn + κ for some constant κ, and then it is straightforward to verify that Eqs. (34) and (35) are also satisfied identically.

Discretization of three-dimensional bi-Hamiltonian flows with one transcendental invariant
There have been several studies of integrable Hamiltonian systems which have transcendental invariants [8,13]. Among the six bi-Hamiltonian flows with transcendental invariants listed in [11] we select the following ones: In [29] the real parameter ξ is restricted to the range |ξ| ∈ (0, 1), but here we need not impose this requirement. Observe that the equations of motion (38) for i = 7, 8, 9 respectively, with Q (7) = Q (8) = P (6) and Q (9) = Q (6) (related to sl(2, R) and to so(3) respectively, see Table 1), : P 12 = {x, y} = 0, P 23 = {y, z} = ξy, P and P (9) = P (7) , with Thus the transcendental invariants in each case are given by H 7 , H 8 and H 9 respectively. (Strictly speaking, H 7 = H 9 is only transcendental when ξ ∈ Q, otherwise it is algebraic.) Moreover, note that the Lie algebra related to P (7) is actually a one-parameter family of Lie algebras, parametrized by ξ; see [29] for more details. It can also be regarded as a four-dimensional Lie algebra, by taking y = log y as a new coordinate and regarding ξ as a central element.
We now construct the Hirota-Kimura type discretizations of the flows E 7 , E 8 and E 9 ; these are denoted using the notation introduced in section 5.
The decoupled equation for x n can be rewritten as a total difference, 1 from which it follows by summation that this is the discrete version of the first equation in (41), to which it tends in the continuum limit By substituting x n given by Eq. (45) into Eq. (43) we get a difference equation for the variable y n , whose solution reads where Γ(z) is the complete gamma function. We can now solve Eq. (46) for the constant β (up to scale) to write it as a function of x n and y n , which gives an explicit transcendental integral: .
Now inserting x n and y n , given respectively by Eqs. (45-46) into Eq. (44) we find a difference equation for z n . Its solution is where In principle, Eq. (47) can implicitly be solved for γ (after first replacing j + τ by (2ǫx j ) −1 everywhere to remove explicit dependence on the parameter τ ), to give another transcendental invariant K 7 , in that case the bi-Hamiltonian structure can be reconstructed by the same formulae as above in cases 1-6; this means that the system dE 7 is completely integrable. Using the formula for H 7 above we can reconstruct one invariant Poisson bracket for this map explicitly, as where Ψ is the digamma function. This bracket has H 7 as a Casimir, and for ξ = ±1 (up to scaling) it reduces to the brackets P (1) (ǫ) and P (2) (ǫ) respectively. It is straightforward to verify that the explicit form of the solution for x n , y n , z n given by Eqs. (45-47) can be used to recover the previous formulae for the discrete systems dE 1 and dE 2 given in Eqs. (18)(19)(20) and (21)(22)(23) by setting ξ = ±1 in the respective cases. 6.2. Explicit solutions to dE 8 . The explicit solution to the equations of motion (39) is given by with H 8 = e −β and K 8 = γ. The discrete-time version of the flow E 8 reads: y n+1 − y n = −ǫ(y n+1 x n + y n x n+1 ) − 2ǫx n x n+1 , (49) z n+1 − z n = ǫ(x n+1 z n + x n z n+1 ) + 2ǫ(y n+1 x n + y n x n+1 ) + 4ǫy n y n+1 . (50) The first equation for x n is identical to that in the previous case, and has the solution x n = (n + τ ) −1 /(2ǫ) as before. By substituting x n into Eq. (49) we get a difference equation for the variable y n , whose solution reads where U n = Ψ (n + τ + 1/2) − Ψ (τ + 1/2) , with Ψ(z) denoting the digamma function as before. This leads to the transcendental invariant H 8 = y n x n + Ψ 1/2 + (2ǫx n ) −1 .
Upon inserting x n as in (45) and y n given by Eq. (51) into Eq. (50) we find a difference equation for z n , whose solution is given by where Similarly to the situation for dE 7 , the system dE 8 has another transcendental integral K 8 which is given implicitly by solving Eq. (52) for γ. This implies that dE 8 is also bi-Hamiltonian and hence completely integrable.
6.3. The system dE 9 . For all values of the parameter ξ, the equations of motion (40) can be reduced to a quadrature, namely Given x(t) determined by this quadrature, y and z are then given by The constants H and K are respectively the values of H 9 and K 9 along an orbit. For certain values of ξ the quadrature can be performed explicitly; for instance, when ξ = 1 it becomes an elementary integral, and the problem reduces to the solution of E 3 , while for when ξ = −1 it becomes an elliptic integral, corresponding to the solution of E 4 , as given in Table 5. The case ξ = 1/2 is also an elementary one, while ξ = −1/2 and ξ = 2 also give elliptic integrals (of the first and third kind, respectively). More generally, for all rational values of ξ this quadrature is an hyperelliptic integral. However, it is straightforward to check that the cases ξ = ±1, which were solved already, are the only ones for which the system has the Painlevé property (i.e. all solutions are meromorphic functions of t in these cases only). In general the solutions have movable algebraic branch points in the complex t plane when ξ ∈ Q, and movable logarithmic branch points when ξ ∈ Q.
The qualitative nature of the solutions is fairly insensitive to the parameter ξ. In fact, for ξ > 0 the trajectories interpolate between the two fixed points (x, y, z) = (0, 0, ± √ 2H), at the north/south poles of the sphere x 2 + y 2 + z 2 = 2H, while for ξ < 0 there are closed periodic orbits. These two types of behaviour are exemplified by each of the explicitly solvable cases ξ = ±1.
We have not attempted to solve this discrete system in the case ξ = ±1. In fact, numerical results for the latter case (as described in the next section) provide evidence for the nonintegrability of the system for generic values of ξ.

Diophantine integrability test
Over the past fifteen years or so there has been a gradual development of methods for testing integrability of maps or difference equations, using such concepts as singularity confinement [10], algebraic entropy [14], Nevanlinna theory [1] and orbit counting over finite fields [35]. In certain limited cases it has been proved that these tests provide necessary conditions for integrability of a map, in a suitable sense, most usually in the setting of algebraic integrability (see [23], for instance), but in general it is an open problem to determine when these tests are effective.
Most recently Halburd proposed an extremely simple criterion for integrability which applies to rational maps defined over Q (or more generally over a number field), which he named the Diophantine integrability test [12]. For a map whose n-th iterate has components x n ∈ Q, written as a fraction x n = p n /q n in lowest terms, the height of x n is defined to be H(x n ) = max(|p n |, |q n |); this is the archimidean height of x n , and the logarithmic height is h(x n ) = log H(x n ). For a map in dimension N, with N components, the height H n of the n-th point on an orbit is defined to be the maximum of the heights of all the components, with h n = log H n being the logarithmic height. Halburd defined a map to be Diophantine integrable if the logarithmic height h n of the iterates of all orbits has at most polynomial growth in n. If we define the Diophantine entropy along an orbit O to be then a Diophantine integrable map is one for which E(O) = 0 for all orbits. Diophantine entropy is somewhat similar to algebraic entropy [14], which measures the height growth of rational functions generated by rational maps. In the latter setting the height of each iterate is just the maximum of the degrees of the polynomials in the numerator and denominator, considered as a rational function in the initial data. However, a huge disadvantage of using algebraic entropy is that one must usually try to guess a recursive relation to generate the degrees of these polynomials. The great advantage of Halburd's test is that it is extremely quick and straightforward to implemement numerically with a computer, and if the map is Diophantine integrable then a plot of log h n against log n should look asymptotically like a straight line (see Figure 1), otherwise it will have an exponential shape (see Figure 7). The main drawback of using the test is that at present it has the status of a distinct definition of integrability, and it is not clear how it is related to other such definitions, like complete integrability in the Liouville-Arnold sense. Despite these drawbacks, it is worth remarking that, at least for maps in two or three dimensions, Diophantine integrability is a necessary condition for algebraic integrability. For example, a two-dimensional map which is algebraically integrable has a conserved quantity whose level sets are algebraic curves. Assuming that each of these curves is irreducible, and that not all orbits of the map are periodic, it was observed by Veselov [41] that they must all have genus zero or one; this follows from a theorem of Hurwitz which says that curves of genus two or more have automorphism groups of finite order [27]. (This argument also extends to the case when the level curves are reducible.) If the curve is rational (genus zero), then the map can be linearized, in which case the logarithmic heights grow linearly, h n ∼ Cn for some constant C, while a curve of genus one is birationally equivalent to an elliptic curve, for which the heights grow as h n ∼ Cn 2 . (See chapter 17 in [6] for an introduction to archimidean heights on elliptic curves, or chapter VIII in [38] for a more general discussion of heights.) Similar considerations apply to algebraically integrable maps in three dimensions, where the algebraic curves are the level sets of two independent integrals, or to systems with N − 1 algebraic integrals in N dimensions (as considered in [23] from the viewpoint of singularity confinement). However, in general these level sets can have two or more irreducible components; see [16] for several examples with two components in three dimensions.
Here we prove that all of the discrete systems constructed here, except for dE 9 , pass the Diophantine integrability test, before presenting numerical results which show more detailed behaviour of the growth of heights for some of these systems. For the theoretical and numerical analysis here it is convenient to set ǫ = 1/2; since the right hand sides of the difference equations are homogeneous (of degree two), this can always be achieved by scaling x n , y n , z n by the same factor. Proof: Without loss of generality we set ǫ = 1/2, as mentioned above, and consider each of the maps with rational initial data x 0 , y 0 , z 0 (and parameter ξ ∈ Q for the case of dE 7 ). This implies that all of the iterates (x n , y n , z n ) of these birational maps are also rational numbers for all n (except on a set of initial data where these maps become singular).
For the maps dE 1 and dE 2 it is clear from the explicit solutions, as given in Eqs. (18)(19)(20) and Eqs. (21)(22)(23) respectively, that in each case the iterates are given in terms of parameters α, β, γ ∈ Q, and these rational iterates have numerators and denominators which grow linearly in n. Hence the logarithmic height satisfies h n = log n + O(1) (subpolynomial growth) for these two maps.
The maps dE 3 and dE 5 are naturally considered together, because their explicit solutions given in Theorem 3 are in terms of hyperbolic functions (or equivalently, exponential functions of n) in each case, which means that the intersections of the level sets of their two integrals are curves of genus zero. This implies that the heights of iterates should grow like h n ∼ Cn. To prove this directly for dE 3 , note that one can eliminate y n from Eq. (36) by setting y n = H 3 x n , and then further eliminate z n between that equation and Eq. (34) to get an expression of the form z n = F (x n , x n+1 , H 3 ) with F being a rational function. This leads to a single recurrence of second order for w n = 1/x n , namely .
The latter recurrence has the conserved quantity , and furthermore admits the linearization which linearizes the system dE 3 ; in terms of the original integrals and solution parameters we find L = (2 + K 3 )/(2 − K 3 ) = cosh 2δ. From the second order linear recurrence (53) it follows directly that the height H(w n ) grows exponentially with n, and hence h(w n ) = h(x n ) ∼ Cn (cf. Figure 1) for some C > 0. Since y n = H 3 x n , and z n can be written as a rational function of x n and x n+1 , it follows that h(y n ) and h(z n ) also have linear growth in n. Analogous arguments apply to dE 5 . Similarly, it is natural to consider the maps dE 4 and dE 6 together, because the intersections of the level sets of their two integrals are curves of genus one; the details for dE 4 are given in the Appendix. For dE 4 each of the coordinates x n , y n , z n of a point on an orbit can be written in terms of Jacobi functions, which are related by a Möbius transformation to the Weierstrass ℘ function. For instance, the solution for z n in (29) is linear in the Jacobi sine, which is an elliptic function of order two with two simple poles in each period parallelogram; this implies that a relation of the form z n = (aX n + b)/(cX n + d) holds, for some constants a, b, c, d, where X n is the nth term in a sequence of X coordinates of points P 0 + nP ∈ E, for an elliptic curve E given in Weierstrass form as Y 2 = X 3 + AX + B (for some A, B). It is known that, as long as P is not a torsion point (which would correspond to a periodic orbit), the height grows like h(X n ) ∼ Cn 2 as n → ∞, where the constant C > 0 only depends on the height of the point P [38]. Since z n is related to X n by a rational map of degree one, it follows that h(z n ) has the same quadratic growth in n, and similarly for h(x n ) and h(y n ). The same arguments apply to dE 6 , this being a special case of the Hirota-Kimura discrete Euler top, whose solutions are most naturally written in terms of Jacobi functions. Finally, for the systems dE 7 and dE 8 we make use of direct estimates of the growth of heights, based on the original maps. For both these systems, note that from the explicit solution we have h(x n ) = h(n + τ ) = log n + O(1). It is convenient to define Y n = y n /x n and Z n = z n /x n in each case, and then note that h(y n ) = h(Y n ) + O(log n), and similarly for h(z n ). From the second part of the map dE 7 we have which implies where the second implication follows by summing over n. Thus h(Y n ) has weaker than quadratic growth in n. Similarly for Z n we have (n + τ − 1/2)Z n+1 = (n + τ + 3/2)Z n + 2ξY n Y n+1 , which implies that and hence h(Z n ) ≤ 1 2 n 2 log n + O(n 2 ), which is weaker than cubic in n. For dE 8 , analogous estimates show that h(Y n ) ≤ n log n + O(n) and h(Z n ) ≤ 2n 2 log n + O(n 2 ), so this system is Diophantine integrable as well.
Having proved that the systems are all Diophantine integrable, we can compare the theoretical results with some numerical experiments. For the system dE 3 we see that the log-log plot gives what we expect: genus zero means linear growth of logarithmic height, so log h(x n ) = log n + O(1); this is evident from the plot of points in Figure 1, which lie asymptotically on a straight line of slope 1. Similarly for the genus one case, we expect log h(x n ) = 2 log n + O(1), and Figure 1 shows points which asymptote to a line with slope 2. In this case the offset, corresponding to the correction at O(1), is function of the height of a point on an associated elliptic curve, and both the point and the curve vary with the initial data of the map.  The theoretical results on the growth of heights for the algebraically integrable systems dE i for 1 ≤ i ≤ 6, as detailed in the above proof, are confirmed by the numerical calculations, and for those cases we have an exact expression for the leading order asymptotic behaviour. Moreover, one can also look at how the asymptote is approached. Taking the system dE 3 for example, h(x n )/n approaches a constant as n → ∞, and from the numerical plot in Figure 4 one can see that this limit is reached in a very uniform manner, in keeping with a correction of O(1/n) to this constant. Similarly, in the case of dE 4 , corresponding to motion on an elliptic curve, we see from Figure 5 that once again the convergence of h(x n )/n 2 to a constant appears to be almost monotone.  Figure 5. Plot of log(h(x n )/n 2 ) versus n for dE 4 with the same data as Figure 2.
The non-algebraically integrable cases, dE 7 and dE 8 , have some extremely interesting features compared with the others. First of all, the method of proof used in Theorem 4 above has not necessarily provided the leading order asymptotics of the logarithmic heights, but has merely given upper bounds on the growth of the form Cn j log n with j = 1 for h(y n ) and j = 2 for h(z n ) in each case. Let us focus on the case of dE 7 . Upon looking more closely at Eq. (54), it would appear that the upper bound for h(Y n ) might be sharp, so that h(Y n ) ∼ n log n (and h(y n ) would have the same leading order asymptotics). However, studies of particular sequences of rational iterates show that cancellations occur between the numerator and denominator of Y n and the prefactor (n + α + 1 − ξ/2)/(n + α + ξ/2), which means that the height of Y n+1 is therefore smaller than the crudest estimate for the upper bound. This weaker growth has a knock-on effect, meaning that the growth of heights of z n also seems to be much weaker than expected. Indeed, Figure 3 suggests that the correct asymptotics should be linear growth in n for the logarithmic heights of both y n and z n , i.e. h(y n ) ∼ C 1 n and h(z n ) ∼ C 2 n for positive constants C 1 , C 2 . For the particular sequence of heights plotted in that figure, a numerical fit shows that log C 1 ≈ 1.64 and log C 2 ≈ 2.64 ≈ log C 1 + 1. We have also plotted log n + log log n for comparison, to show how the upper bound for h(y n ) fails to be sharp. Another surprising feature of this system is that for different choices of initial data we find (to within numerical accuracy) the same values of C 1 and C 2 ; this is in contrast to the algebraically integrable setting described above, where the coefficient in front of the leading order term is dependent on the initial data. Thus we might conjecture that for this map C 1 , C 2 are independent of initial data, and also that C 2 = e C 1 holds identically, in which case there should be some deeper arithmetical explanation for this asymptotic behaviour. Similarly to the case of dE 7 , numerical results for the system dE 8 also show linear growth of logarithmic heights.
Supposing that the numerical observation of linear growth of h n for these discrete systems with transcendental invariants is indeed the correct asymptotic behaviour, it is then interesting to look at how h n /n approaches a constant value. The results we find are in  Figure 6. Plot of 400 iterates of log(h(x n )/n) versus n for dE 7 with the same initial data and parameters as in Figure 3.
stark contrast to the algebraic setting: rather than the almost monotone convergence seen in the previous examples, for dE 7 we find that h n /n shows rapid fluctuations which persist for increasing values of n. These fluctuations in the asymptotics are somewhat reminiscent of the "random"-looking error terms that appear in some famous arithmetical functions, such as the difference between the prime-counting function π(n) and the logarithmic integral [25]. It would be interesting to know whether these fluctuations might provide a means of characterizing the difference between discrete systems which are algebraically integrable and those with transcendental invariants. For comparison with the Diophantine integrable examples above, in Figure 7 above we have plotted the growth of h n for a particular case of the discrete Lotka-Volterra system due to Kahan, which is the degree two birational map given in Eq. (2). This figure shows that the logarithmic height seems to grow exponentially, indicating non-integrability of this system. Indeed, the heights of iterates grow so fast that even on a fairly new computer it took 1 hour to calculate the heights of 17 rational iterates with Maple; the value of h(x 17 ) is of the order of 10 325009 in this case. Upon examining the data used in Figure 7 more carefully, it is apparent that h(x n+1 ) ≈ 2 h(x n ) to very good accuracy, so we expect h n ∼ C 2 n . This would mean that the logarithmic height essentially doubles with each step, giving a Diophantine entropy of log 2 for generic (aperiodic) orbits. This appears to be the same as the algebraic entropy of the map when ǫ 2 = 1, which was calculated by A. Ramani [33]. In any case, this is the entropy value that one would expect for a generic (non-integrable) birational map of degree two.
Finally we should mention the results of numerical calculations of the growth of heights for the map dE 9 for various cases with ξ = ±1 (that is, excluding the two special cases where the map is already known to be algebraically integrable). For generic rational values of the parameter ξ we find that the map dE 9 is not Diophantine integrable, but rather the Diophantine entropy is log 3 for generic orbits, this being the typical value to be expected for a non-integrable birational map of degree three. (See Figure 8 for an illustration of the exponential growth of logarithmc heights when ξ = −1/2.) These numerical results suggest that while that the original continuous system E 9 is algebraically integrable (in the sense of having two independent algebraic integrals), the corresponding discrete system is not. We shall return to this point in our conclusions.

Concluding remarks
In this paper we studied three-dimensional birational maps which provide integrable time-discretizations of quadratic bi-Hamiltonian flows associated with pairs of real threedimensional Lie algebras, as presented in [3,11]. We have shown that for the six cases of continuous flows which are algebraically integrable, the Hirota-Kimura type discretization provides maps admitting two independent rational integrals of motion, in involution with respect to a pair of compatible Poisson tensors. We have also provided explicit solutions of the resulting discrete systems dE i for i = 1, . . . , 6, which are given in terms of either rational, hyperbolic or elliptic functions in each case. These results confirm the conjecture that the property of algebraic integrability is preserved by this discretization scheme.
We have also applied the same procedure to two cases of integrable continuous flows in three dimensions having one rational and one transcendental integral of motion, for which the resulting maps, dE 7 and dE 8 , admit explicit solutions in terms of rational functions and either gamma or digamma functions. In each of these cases, we have found an explicit formula for one transcendental integral of the map, but the second integral is only defined implicitly by the solution. Nevertheless, this is sufficient to assert that the latter two cases are also completely integrable in the Liouville-Arnold sense. Therefore the Kahan-Hirota-Kimura discretization scheme preserves integrability even in these transcendental cases. However, for another example of a continuous integrable system with one transcendental integral, it appears likely that the corresponding discretization dE 9 is not integrable for generic values of the parameter in the map.
In an attempt to gain a better understanding of the difference between the algebraic and transcendental cases, we have also analyzed all of these discrete systems from a different viewpoint, within the arithmetical setting of Diophantine integrability. So far, the Diophantine integrability test has been applied to various algebraically integrable systems and discrete Painlevé equations [12], as well as to certain birational maps that are not algebraically integrable and fail the test [18,19]. Our theoretical results show that dE 7 and dE 8 provide examples of discrete integrable systems with transcendental invariants that are also Diophantine integrable, in the sense defined by Halburd. Moreover, more detailed numerical results suggest that these discrete integrable systems might be distinguished from algebraically integrable maps by the manner in which the logarithmic heights converge to their leading order asymptotics. This asymptotic behaviour deserves to be studied more carefully in the future.
A further interesting point is the comparison with Kahan's discretization of the Lotka-Volterra system, which is an integrable flow in the plane with a transcendental integral.
The discrete system provides a non-standard symplectic integrator of this flow, and seems to preserve the qualitative features of the continuous counterpart. However, the numerical results indicate that the discrete system is not Diophantine integrable, which adds further evidence to the conjecture that (for generic values of ǫ) it should not be Liouville integrable either.
Similar considerations apply to the discrete system dE 9 : numerically it appears that for generic values of the parameter ξ, it fails the Diophantine integrability test. Since the continuous system E 9 has algebraic integrals for ξ ∈ Q, this means that in general algebraic integrability (in the weakest sense of the term) is not preserved by the Kahan-Hirota-Kimura discretization. Actually one could already observe this for the system E 7 , since it has the integral H 7 which is algebraic for all ξ ∈ Q, but the integral H 7 for dE 7 is transcendental unless ξ ∈ Z. This suggests that one should impose a much stronger notion of algebraic complete integrability (a.c.i.) if this is to be preserved by the discretization scheme. For instance, one can require that the generic level sets of the integrals are smooth abelian varieties, possibly extended by (C * ) m for some m. An excellent discussion of various different definitions of a.c.i. can be found in chapter V of [40].
As for the case of the discrete three-dimensional Euler top (see [15,30]), there is one other standard attribute of integrable systems that remains to be found for the maps dE i for i = 1, . . . , 8, namely their Lax representation. This is an open problem which deserves further investigation.
There are three other bi-Hamiltonian flows in the list of Gümral and Nutku, all of which have one transcendental invariant. Preliminary results suggest that their Kahan-Hirota-Kimura type discretizations are qualitatively similar to the system dE 9 , and we expect that these maps are not Liouville integrable. We reserve the study of these systems for future work.
Finally we remark that another non-standard symplectic integrator for the Lotka-Volterra model, with similar numerical properties, has been given by Mickens in [26]. It would be very interesting to see if the approach to discretization proposed by Mickens shares some of the remarkable properties of Kahan's.

Acknowledgments
AH and MP are grateful to the organisers for supporting their attendance at the Miniworkshop on Integrable Systems at Università di Milano-Bicocca in September 2007, where this collaboration began. Both authors are grateful to Yuri Suris for helpful correspondence. AH would also like to thank Kim Towler for a careful reading of the text, and Gavin Brown for useful discussions.

Appendix: solution of the system dE 4 in elliptic functions
Here we derive the formulae (27)(28)(29) corresponding to the solution of dE 4 , as in Theorem 3. The first observation to make is that the curve V in affine three-dimensional space defined by the equations V : xy = H 1 + ǫ 2 (x 2 + y 2 ) , x 2 + y 2 + z 2 = 2K 1 + ǫ 2 (x 2 + y 2 ) , corresponding to the intersection of the level sets H 4 = H, K 4 = K, has genus one (at least for generic values of H and K). To see this, note that V is a double cover of the curve C : xy = H 1 + ǫ 2 (x 2 + y 2 ) in two dimensions, via the covering map π : V → C (x, y, z) → (x, y) which is ramified over the four points (x, y) ∈ C obtained from the simultaneous solutions of xy = H(1 + ǫ 2 (x 2 + y 2 )), (1 − 2Kǫ 2 )(x 2 + y 2 ) = 2K (when z = 0). Since the curve C is a conic (genus zero), it follows from the Riemann-Hurwitz formula [27] that V is (the affine part of) a curve of genus one. The first order recurrence relations for x n , y n , z n , namely x n+1 − x n = −ǫ(x n+1 z n + x n z n+1 ), (55) y n+1 − y n = ǫ(y n+1 z n + y n z n+1 ), (56) z n+1 − z n = 2ǫ(x n+1 x n − y n y n+1 ), correspond to a birational map from this curve to itself, inducing an automorphism of an isomorphic elliptic curve (i.e. a plane curve defined by a Weierstrass cubic), and it follows that x n = X(u + nv) for a suitable elliptic function X, and similarly for y n , z n . One can see some points on a real connected component of such a curve in Figure 9.  Figure 9. Plot of the first 1000 points on the orbit of the integrable map dE 4 with initial conditions x 0 = 7/3, y 0 = 11/13, z 0 = 23/47.
From the equations for V it is easy to see that the functions corresponding to x n , y n , z n each have simple poles at the same places, and they are elliptic functions of order two. These facts suggest that it may be most convenient to write the formulae in terms of Jacobian (rather than Weierstrassian) elliptic functions. Indeed, if we set then the equations for V become (1 − 2Hǫ 2 )s 2 − (1 + 2Hǫ 2 )d 2 = 4H, (1 − 2Kǫ 2 )(s 2 + d 2 ) + 4z 2 = 4K, which are reminiscent of (linear combinations of) the quadratic relations sn 2 (u) + cn 2 (u) = 1, k 2 sn 2 (u) + dn 2 (u) = 1, for Jacobi functions. In order to obtain the formulae (27)(28)(29), it is instructive to take a detour through Jacobi theta functions, by deriving bilinear equations from Eqs. (55-57). Upon setting x n = A n + B n 2ǫD n , y n = A n − B n 2ǫD n , z n = C n ǫD n , the system (55) is equivalent to the following three bilinear equations: One should hesitate to call (58) the Hirota bilinear form of Eqs. (55-57), because there are four unknowns (tau-functions) A n , B n , C n , D n but only three equations, so the system is underdetermined. Despite this apparent problem, we can solve this bilinear system in terms of Jacobi theta functions ϑ j , j = 1, . . . , 4, by comparing these equations with the identities in exercise number 3 on page 488 of [42]; the first of these is the relation ϑ 1 (u ± v)ϑ 2 (u ∓ v)ϑ 3 ϑ 4 = ϑ 1 (u)ϑ 2 (u)ϑ 3 (v)ϑ 4 (v) ± ϑ 3 (u)ϑ 4 (u)ϑ 1 (v)ϑ 2 (v) the last is and there are four other relations of this kind, for different permutations of the four indices.
Here ϑ j without argument denotes a theta constant (i.e. ϑ j = ϑ j (0), which depends on the modulus k). By taking the sum of the two equations given in (59) with opposite choices of ± signs, and similarly taking the difference of the two equations specified by (60), one sees that both ϑ 1 (u+v)ϑ 2 (u−v)+ϑ 1 (u−v)ϑ 2 (u+v) and ϑ 3 (u+v)ϑ 4 (u−v)−ϑ 3 (u−v)ϑ 4 (u+v) are proportional to ϑ 1 (u)ϑ 2 (u), modulo v-dependent factors. Thus if v is regarded as a fixed constant, and the shift u → u+2v is identified with n → n+1, then the first equation in (58)  , , it follows that the solution of the difference equations (55) has the form x n = λ dn(κ + 2nδ) + µ cn(κ + 2nδ), y n = λ dn(κ + 2nδ) − µ cn(κ + 2nδ), z n = ν sn(κ + 2nδ), for constants δ, κ and suitable prefactors λ, µ, ν which are given in terms of δ and the theta constants. The expressions (27)(28)(29) can also be verified directly from the addition formula for sn, namely sn(u + v) = sn(u)cn(v)dn(v) + sn(v)cn(u)dn(u) 1 − k 2 sn 2 (u)sn 2 (v) , as well as the analogous formulae for cn and dn. Using (61) to calculate sn(u+v)−sn(u−v), and then setting u → κ+ (2n+ 1)δ, v → δ, gives an expression for the left hand side of the third equation in (55), and performing analogous computations for the right hand side and for the other two difference equations allows the prefactors λ, µ, ν to be determined directly in terms of Jacobi functions with argument δ, in agreement with (27)(28)(29). Finally, note that the solution depends on the required number of arbitrary constants, namely the three parameters δ, κ, k. The parameter δ and the modulus k are determined by the values of the integrals H 4 and K 4 , by solving the relations (30) as a system for k and snδ and then performing the elliptic integral δ =