Continued fractions and Hankel determinants from hyperelliptic curves

Following van der Poorten, we consider a family of nonlinear maps which are generated from the continued fraction expansion of a function on a hyperelliptic curve of genus $\mathrm{g}$. Using the connection with the classical theory of J-fractions and orthogonal polynomials, we show that in the simplest case $\mathrm{g}=1$ this provides a straightforward derivation of Hankel determinant formulae for the terms of a general Somos-4 sequence, which were found in a particular form by Chang, Hu and Xin, We extend these formulae to the higher genus case, and prove that generic Hankel determinants in genus two satisfy a Somos-8 relation. Moreover, for all $\mathrm{g}$ we show that the iteration for the continued fraction expansion is equivalent to a discrete Lax pair with a natural Poisson structure, and the associated nonlinear map is a discrete integrable system. This paper is dedicated to the memory of Jon Nimmo.


Introduction
The Somos-4 recurrence is given by The surprising observation of Somos was that when α = β = 1 and the four initial values τ 0 , τ 1 , τ 2 , τ 3 are all 1, the recurrence (1.1) generates a sequence of integers [49], beginning with A proof of this fact was eventually published [36], but a better understanding of the mechanism by which such rational recurrences can yield integer sequences came from the observation that (1.1) exhibits the Laurent property [21]: the iterates are Laurent polynomials in the initial values with integer coefficients, that is to say τ n ∈ Z[τ ±1 0 , τ ±1 1 , τ ±1 2 , τ ±1 3 , α, β] for all n, which makes it obvious why (1.2) consists entirely of integers. The Laurent phenomenon [19] eventually appeared as a key property of the distinguished generators (cluster variables) in Fomin and Zelevinsky's cluster algebras [18], which are constructed by a recursive process called mutation, and Fordy and Marsh showed how the Somos-4 recurrence and various higher order analogues arise from cluster mutations starting from quivers of a particular type [20].
Cluster algebras fit within a broader setting of Laurent phenomenon algebras [33], leading to a wide variety of nonlinear recurrences that exhibit the Laurent property [2]; and it is possible to reverse engineer a rational recurrence to generate an integer sequence [15]. However, Somos-4 sequences have some very special features which are a consequence of the fact that each such sequence is associated with a sequence of points P 0 + nP on an elliptic curve E, and this leads to an analytic formula for the terms of the sequence. The following result was proved in [23].
The formula (1.3) also makes sense in the degenerate case when the discriminant g 3 2 −27g 2 3 = 0. Although the above result is formulated over the complex numbers, its algebraic content -associating a solution of (1.1) with a sequence of points on an elliptic curve -is valid in any field over which the initial values and coefficients α, β are defined (up to appropriate adjustments in characteristic 2 or 3); this was described independently by Swart [51], and also, in terms of a quartic model for E, by van der Poorten [43]. This underlying algebraic structure has many consequences, including the existence of higher order relations between the terms [24,34,44], and more refined versions of the Laurent property which produce large families of integer sequences [25]. From this point of view, Somos-4 sequences are natural extensions of Ward's elliptic divisibility sequences [55], which correspond to the special case P 0 = ∞ (the identity element in the group law of E), and generalize the arithmetical properties of Fibonacci or Lucas sequences to a nonlinear setting [16]. Aside from their intrinsic interest for certain problems of an arithmetical [48] or Diophantine nature [9], Somos sequences and their higher order analogues appear in discrete integrable systems, underlying many integrable maps [28], especially via reductions of the discrete Hirota equation (bilinear discrete KP, also known as the octahedron recurrence) [27] or Miwa's equation (bilinear discrete BKP, or the cube recurrence) [17]. They also arise in solvable models of statistical mechanics and quantum field theory, such as the hard hexagon model, as mentioned in [47], or dimer models and quiver gauge theory [14].
It was conjectured by Somos, and later proved by Xin [56], that the terms of the sequence (1.2) have another explicit expression that is rather different from (1.3), being given by the Hankel determinants D n = det(s i+j−2 ) i,j=1,...,n , (1.4) where the entriess j are obtained from the function η = η(x) satisfying the algebraic equation To be precise, solving (1.5) for η with a fixed choice of square root, one should take the functionG = η/x − 1, and expand it as η x − 1 = j≥1s j−1 x j = x + x 2 + 3x 3 + 8x 4 + 23x 5 + · · · , (1. and so on, where the matrix entries are generated by the recursioñ s j =αs j−1 +βs j−2 +γ j−2 i=0s isj−2−i , j ≥ 2, (1.8) withα = 2,β = 0,γ = 1 ands 0 =s 1 = 1. It was further conjectured by Barry [4] that the Hankel determinants D n formed from a particular family of sequences (s j ), defined by the recursion (1.8) with s 0 = 1,s 1 =α, satisfy the Somos-4 recurrence (1.1) with coefficients α =α 2γ2 , β =γ 2 (β +γ) 2 −α 2γ3 , and this was proved by Chang and Hu using identities for block Hankel determinants [10]. The latter result does not overlap with that of Xin, since the conditions on the coefficients and initial conditions do not include the original sequence (1.2). However, it was subsequently shown by Chang, Hu and Xin that, for any Somos-4 sequence with two adjacent initial values equal to 1, the terms with positive index n are given by a Hankel determinant of the form (1.4), where the entriess j satisfy a recursion of the form (1.8), for a suitable choice ofα,β,γ,s 0 ,s 1 [11].
In this paper we start from van der Poorten's construction in [43] for Somos-4, based on the continued fraction expansion of a function on a quartic curve of genus one, and the results of [45,46]. In the latter work, the continued fraction approach was extended to hyperelliptic curves of higher genus g, defined by a polynomial of even degree 2g + 2, but with partial success: a Somos-6 relation was obtained in genus two, but only in a special case. The continued fraction expansion of a hyperelliptic function of a certain type is described in section 2 (for other related results, and the connection with the geometry of the Jacobian of the curve, see [1,5,6,22,42]). Next, in section 3, the recursion for the continued fraction is reformulated as a discrete dynamical system defined by a matrix linear problem (a discrete Lax pair), and we state the first main result, Theorem 3.1, which says that this nonlinear dynamical system is integrable in the sense that it satisfies a discrete analogue of the Liouville-Arnold theorem from classical mechanics [3], having an invariant symplectic structure and a sufficient number of first integrals in involution with respect to the corresponding Poisson bracket (see [8,35,53] for the precise definition of a discrete integrable system). We present explicit details of the symplectic map and first integrals for the cases g = 1 and g = 2, while the complete proof for any g is deferred until section 6. Section 4 is concerned with the derivation of Hankel determinant formulae for the solutions of the nonlinear system, based on the classical theory of J-fractions and orthogonal polynomials, which are presented in a uniform fashion for any genus g. (Note that Hankel determinants and continued fractions have appeared in the solution of many other integrable systems, particularly those of Toda type, and Painlevé equations [12,31,50], and there are more recent results in the broader context of Padé approximants associated with isomonodromic deformations [37].) Subsequently, in section 5, we show how the Hankel formulae obtained generalize the results of Chang, Hu and Xin on Somos-4 sequences: even in the elliptic case g = 1, the results are more general than [11], since the Hankel determinants depend on an additional free parameter, and the formulae from the previous section extend to negative indices n. In genus two we prove that, for generic parameter values, the corresponding Hankel determinants satisfy a Somos-8 relation (Theorem 5.4), and indicate how van der Poorten's Somos-6 recurrence arises as a special case. We also present a precise conjecture that provides an analytic formula analogous to (1.3) for the Hankel determinants in genus two, and briefly explain how an appropriate higher genus analogue of this conjecture implies the existence of Somos recurrences for all values of g.
In section 6 we employ a space of 2 × 2 Lax matrices, related to those in section 3 by a gauge transformation, which admits a natural Poisson structure, and construct a completely integrable system on this phase space, given by a set of commuting flows defined by suitable Hamiltonian functions. We then show how the nonlinear map coming from the continued fraction arises from a Poisson map on this phase space, which preserves the same Hamiltonians and Casimirs as the continuous system. The map we obtain is somewhat reminiscent of the Bäcklund transformation (BT) for the even Mumford systems, introduced in [32], except that the entries of the Lax pair have a different degree structure, and (in contast with the BT, which is a multivalued correspondence) it is an explicit birational map. Some determinantal identities that directly yield the formulae for the coefficients of J-fractions are also presented in an appendix.

Continued fractions for hyperelliptic functions
Following van der Poorten [43,45,46], we consider a hyperelliptic curve defined by C : for a pair of polynomials In addition to the affine points (X, Y ) ∈ C 2 satisfying (2.1), one can adjoin two points at infinity, ∞ 1 and ∞ 2 , such that, in terms of a local parameter t at each of these points, X = 1/t and Y ∼ ±1/t g+1 respectively. Thus, in the generic case that all the roots of F (X) are distinct, one obtains a compact Riemann surface of genus g, also denoted C. In the associated function field for certain polynomials P 0 , Q 0 in X of degrees g + 1, g respectively, taking the form and we impose the additional requirement that To compute the continued fraction of any element Y 0 ∈ F, we take its expansion in the neighbourhood of the point ∞ 1 ∈ C, given by a power series in X −1 ; this can be viewed as an element of C((X −1 )). Then the continued fraction is where, for any element of C((X −1 )), the floor symbol denotes the polynomial part, and the remainder is a series in positive powers of X −1 . Thus, by iterating the standard recursion Y n = a n + 1 Y n+1 , a n = Y n (2.7) for n = 0, 1, 2, . . ., one obtains the successive partial quotients a n (X) in the continued fraction (2.6) above, which are polynomials in X. Now let us describe in detail the form of the continued fraction expansion for the particular type of function Y 0 specified by (2.3). Because the neighbourhood of ∞ 1 is being considered, with X → ∞, it follows that Y ∼ A(X) ∼ X g+1 , so Y + P 0 ∼ 2A(X) and hence Y 0 ∼ 2X/u 0 . Thus a 0 is linear in X. Moreover, by (2.5), there is some polynomial Q −1 (X) of degree g such that If a polynomial P 1 (X) of degree g + 1 is defined by then it follows from (2.5) that Q 0 |Y 2 − P 2 1 , so there is some polynomial Q 1 (X) of degree g such that and also Q 1 |Y 2 − P 2 1 . Then by induction, at each stage of the recursion (2.7) we find linear partial quotients a n (X) = 2(X + v n )/u n , (2.9) for a sequence of polynomials P n , Q n of degrees g + 1 and g respectively, with P n (X) = A(X) + 2d n X g−1 + · · · , Q n (X) = u n X g + · · · . (2.11) Note from (2.9) that u n = 0 is required for the recursion to make sense at each stage. Moreover, at each stage, Y n has positive degree in X, andȲ n = (−Y + P n )/Q n (its image under the hyperelliptic involution) has negative degree; in the terminology of van der Poorten, Y n is reduced [46]. Observe that in the equation (2.1) for C there is always the freedom to shift X → X+ const, which replaces F (X) by another monic polynomial of the same degree. Henceforth we will exploit this freedom in order to remove the coefficient at order X g in F , which means that for some constants k (j) . This choice is convenient because in the continued fraction expansion it means that In other words, modulo factors of 2 and u n , the next to leading order term in Q n completely fixes the constant term in the partial quotient (2.9), and we will always assume that A has the form (2.12) so that this is the case. (We will make another comment about this later, when we discuss Poisson brackets.) As we shall see in the next section, the recursion (2.7) and the relations (2.10) together yield a set of coupled nonlinear recurrences for the coefficients appearing in P n and Q n . For the time being we derive just one such relation, by considering the second equality in (2.10), which is equivalent to the identity Y 2 = P 2 n + Q n Q n−1 . (2.14) If we make use of (2.1) together with (2.11), and cancel A 2 from both sides, then we have so the leading order term, at order X 2g , gives the formula 4d n + u n u n−1 = 0. (2.16) The above identity can be used to eliminate all prefactors involving u n wherever they appear, so that the interesting dynamical relations that remain will only involve coefficients of P n and Q n /u n . Note also that, from (2.9), the recursion breaks down if d n vanishes at some stage. There is a natural geometrical interpretation of the iteration that produces the continued fraction expansion (2.6) for (2.3), which goes back to work of Adams and Razar on the elliptic case [1], and was further generalized by Bombieri and Cohen in the setting of Padé approximation of functions on algebraic curves of general type [6]. Let us consider the function 1 the roots of Q 0 (X) and Q 1 (X), respectively. If we also set y (j) n = P n (x (j) n ), j = 1, . . . , g for n = 0, 1, then under the Abel map each of the degree g divisors corresponds to a point on the Jacobian variety Jac(C) ∼ = Sym g (C), identified with the g-fold symmetric product of the curve [38]. From its expression as Y 0 −a 0 , the poles of G lie at the points (x 1 ) and ∞ 1 , where Y 1 has poles. Therefore D 0 and D 1 are related by the linear equivalence and by the same argument the shift n → n + 1 in each line of the continued fraction is equivalent to a translation in the Jacobian by the divisor class of ∞ 2 − ∞ 1 .

Lax pair and nonlinear system
The recursion (2.7) and the relations (2.10) can be reformulated in terms of a linear system, which makes their structure much more transparent. To do this, it suffices to introduce projective coordinates, setting Y n = ψ n /φ n , Ψ n = (ψ n , φ n ) T , and substituting into (2.10), which leads to the eigenvalue problem for the Lax matrix Upon substituting the ratio of the projective coordinates into (2.7), and fixing an arbitrary multiplier, the fractional linear relation between Y n and Y n+1 separates into two linear equations, which can be written as where, taking the standard formula from the classical theory of continued fractions for real numbers, we may set M n = a n 1 1 0 .
(3.4) (The choice of multiplier means that the matrix M n is only defined up to overall scaling M n → λ n M n for some arbitrary n-dependent quantity λ n .) The compatibility of the linear system consisting of (3.1) and (3.3) is the discrete Lax equation which produces two non-trivial conditions, namely where we have substituted the expression (2.9) for the partial quotient a n . Note also that, because it is equivalent to conjugation by the nonsingular matrix M n , (3.5) is an isospectral evolution, preserving the spectral curve det(Y 1 − L n (X)) = 0, which reproduces the formula (2.14) for all n. Now, from the terms in the continued fraction, we can expand with two particular coefficients being specified as in terms of the notation used previously. Upon substituting (3.7) into (3.6), the terms involving u n can be replaced using (2.16) in the first relation, and cancelled from the second relation, to yield a set of recurrences for the coefficients π Let us introduce the g-tuples of affine coordinates. Due to (2.12) and (3.8), the coefficients at order X g in (3.9) and at orders X g+1 and X g in (3.10) provide only tautologies, so that altogether there are 2g non-trivial relations between the components of π n , π n+1 , and ρ n , ρ n±1 . To be precise, these relations mean that the 2g quantities π n+1 and ρ n+1 can be calculated as rational functions of the components of π n , ρ n−1 , and ρ n , which (together with the relation (2.16) for the prefactors) shows how the entries of L n+1 are determined from those of L n . Similarly, in the reverse direction n + 1 → n, these relations mean that the entries of L n can be obtained as rational functions of the entries of L n+1 . In the above form, the map corresponding to the shift n → n + 1 from one line of the continued fraction to the next can be interpreted as a discrete dynamical system, where (ignoring the prefactors u n ) this can be viewed as a birational map (π n , ρ n−1 , ρ n ) → (π n+1 , ρ n , ρ n+1 ) in dimension 3g. However, at the expense of introducing more parameters, one can use the equation for the spectral curve (2.14) to eliminate g coordinates and rewrite this in terms of a birational map in dimension 2g. In particular, in the explicit formula (2.15) the leading order (X 2g ) term gives (2.16), while the coefficients at each order from X 2g−1 down to X g can be used to rewrite ρ n−1 in terms of the components of π n , ρ n , as well as the coefficients in appearing in A(X), and also u, the leading coefficient of R(X) on the left-hand side. The remaining g coefficients of R(X), which appear at orders X j , j = 0, . . . , g − 1, can then be written as rational functions of the components of π n , ρ n and the other parameters, and these g quantities are independent of n. In this way, we arrive at a birational map ϕ : (π n , ρ n ) → (π n+1 , ρ n+1 ) (3.11) which has g conserved quantities. In fact there is more that one can say about this map: it turns out to be symplectic, and integrable in the sense of a suitable discrete analogue of Liouville's theorem [8,35,53]. Observe that the expression (2.14) is symmetrical in Q n−1 and Q n , so one can just as well use it eliminate the components of ρ n from (3.9) and (3.10), to obtain a birational map ϕ : (π n , ρ n−1 ) → (π n+1 , ρ n ). (3.12) Clearly the latter map is conjugate to ϕ, in the sense that there is a birational transformation χ such that ϕ = χ −1 • ϕ • χ, and the above theorem applies equally well toφ. The general proof this theorem is given in section 6, where we make use of a Poisson structure for the Lax matrices (3.2). For now, we just give explicit details for g = 1 and 2.
Example 3.1. The case g = 1: In the genus one case, following [43], we write for arbitrary parameters f, u, v defining the quartic curve in the (X, Y ) plane. There are only two non-trivial relations from (3.9) and (3.10), given by which define a birational map in 3 dimensions, that is However, using the equation for the curve, and removing an overall factor of 4, the formula (2.15) becomes The first non-trivial relation, at order X, gives which allows v n+1 to be rewritten as a function of d n , v n and the parameters f, u. Hence, making use of (3.14), this yields a map in the plane, that is The above map preserves the symplectic form so replacing v n−1 as before and setting H = −uv we see that is a conserved quantity for ϕ, so it is an integrable symplectic map in two dimensions.
Example 3.2. The case g = 2: In the genus two case, adopting the notation in [45], we write for arbitrary parameters f, g, u, v, w defining the sextic curve C : From (3.9) and (3.10) there are four relations that define a birational map in 6 dimensions, namely is used to obtain d n+1 , and then (3.23) produces an expression for e n+1 , which allows v n+1 and w n+1 to be calculated from (3.20) and (3.21), respectively. In order to obtain a map in 4 dimensions, one can use (2.15) to eliminate v n−1 and w n−1 , giving (3.11), or instead eliminate v n and w n , to obtain (3.12); here we take the latter option. To be precise, compared with the quantities used in (3.12) we have made an invertible change of coordinates, that is (π , but by a slight abuse of of notation we will use the same symbolφ to denote the map that describes the shift n → n + 1 in terms of (d n , e n , v n−1 , w n−1 ). There are four non-trivial relations coming from (2.15), given by Using (3.24) and (3.25), together with (3.22) and (3.23), we find the map where the shifted variables are given explicitly in terms of the previous ones by This four-dimensional map preserves a nondegenerate Poisson bracket, given by hence it is symplectic. (As we shall see in section 6, in a different set of coordinates, with e n replaced by π  also rewrite this as a map (d n−1 , d n , v n−1 , v n ) → (d n , d n+1 , v n , v n+1 ), so that it takes a simpler form as a pair of coupled recurrence relations of second order, that is (3.29) and in these coordinates (up to an overall choice of scaling) the symplectic form is

Orthogonal polynomials and Hankel determinants
After removing the first term a 0 , the continued fraction expansion (2.6) becomes (4.1) A continued fraction of this form, where each partial quotient a j = a j (X) is a linear function of X, is called a J-fraction [54]. If we multiply the main numerator and denominator by u 1 /2 and apply (2.16), then this can be rewritten in its more classical form, that is from which we see that, up to an overall prefactor of −2/u 0 , it is completely determined by the quantities d n and v n . In this section we apply standard results on J-fractions and associated orthogonal polynomials, which lead directly to formulae for the quantities d n and v n in terms of ratios of Casorati determinants, and Hankel determinants in particular. This generalizes certain results obtained for g = 1 by Chang, Hu and Xin [11] to hyperelliptic curves of any genus g.
In the neighbourhood of the point ∞ 1 ∈ C, the function G defined by (2.17) has the series expansion which can be used to define a linear functional < · >: F → C according to for any function Φ ∈ F, where the integral is taken along any sufficiently small closed contour around ∞ 1 (anticlockwise in the X plane) that does not encircle the poles of G at (x and ∞ 2 . In other words, G can be regarded as a moment generating function with moments < X j >= s j , j = 0, 1, 2, . . . , (4.5) although in general, for complex s j , this may not be associated with some positive measure. The linear functional (4.4) also defines a scalar product between any pair of functions Φ, Ψ, that is The convergents of (4.1) are the sequence of rational functions of X given by where by convention one can take p 0 = 0, q 0 = 1, (4.7) followed by and for all j ≥ 1 p j (X) and q j (X) are polynomials of degrees j − 1 and j in X, respectively, where without loss of generality each q j is taken to be monic. The recursion for the convergents is essentially controlled by the same matrix (3.4) that appears in the classical theory of continued fraction expansions of numbers in R. However, the entries of this matrix must be scaled in order to ensure that the q j are monic, making use of (2.16) to remove dependence on the prefactors u j , to yield which is equivalent to the three-term linear recurrence relation for the sequence of polynomials (q n ), and the same for the sequence (p n ), with the initial conditions (4.7) and (4.8). From the linear recurrence it is clear that p n /s 0 is monic for each n.
Due to the fact that each a j = Y j is linear in X, it is straightforward to see by induction that the nth approximant p n /q n satisfies G − p n q n = 1 and since q n has degree n this implies that Given the requirement on the degrees of p n and q n , by considering the terms at orders X n−1 , X n−2 , . . . , X −n , the equation (4.10) provides n linear equations that determine the non-trivial coefficients of q n in terms of the coefficients s j in (4.3), and a further n linear equations for the non-trivial coefficients of p n in terms of those of q n and the s j . This leads to a standard formula for q n , given explicitly in determinantal form as s n−2 · · · · · · s 2n−4 s 2n−2 s n−1 · · · · · · s 2n−3 s 2n−1 . (4.13) We refer to the latter as a shifted Hankel determinant. By convention we set ∆ 0 = 1 and ∆ * 0 = 0. We can now state our main result about Hankel determinants and orthogonal polynomials. 14) where the entries are determined recursively by s 0 = u 1 /2 and, for j ≥ 1, Furthermore, the polynomials q j that appear as the denominators of the convergents of the J-fraction are orthogonal with respect to the scalar product (4.6) associated with the series (4.3), that is Proof: The formulae in (4.14) are classical expressions for the coefficients appearing in the linear recurrence (4.9) for orthogonal polynomials. A direct proof is obtained by substituting (4.11) into the three-term recurrence and expanding in powers of X: the determinantal expression for v n appears immediately at order X n , while at order X n−1 one finds a formula for d n as a combination of four terms that can be condensed into a single ratio by applying various identities for determinants (which are collected in an appendix for completeness). To find the entries of the Hankel determinants recursively, note from (2.8) that using (2.14), so that G satisfies the quadratic equation Upon substituting the series (4.3) for G and expanding the polynomial coefficients in (4.17) with the notation in (2.12) and (3.7), and making use of (2.16) to replace u 0 , the recursion (4.15) results. Finally, the orthogonality of the sequence of polynomials (q n ) follows by a standard inductive argument using the three-term recurrence, also making use of the moments (4.5) with (4.11) to expand ∆ 2 n < q 2 n > as a sum of n+1 products of determinants, only one of which is non-vanishing, which yields < q n , q n >= ∆ n+1 /∆ n = h n .
Example 4.1. The recursion for moments in the elliptic case: For g = 1 we use the same notation as in Example 3.1, and for the recursion (4.15) we have for all j ≥ 2. To illustrate this with a particular numerical example, let us pick the curve so that f = −3, u = −1, v = −2, and set u 0 = −2, d 0 = 1, v 0 = −1. This fixes the function This produces the sequence of moments which should remind the reader of (1.2).
Example 4.2. The recursion for moments in genus 2: With the notation of Example 3.2, the recursion (4.15) for g = 2 has initial values and subsequent coefficients in the series expansion (4.3) are determined for j ≥ 3 by As a particular example, consider the curve  Fig. 1, which gives and the recursion for the coefficients (moments) in the above expansion is with initial values s 0 = 1, s 1 = 0, s 2 = 2. In this case, the sequence of moments (s j ) begins with The map (3.11), or equivalently (3.12), corresponds to the recursion for the continued fraction expansion, and since this map is birational, it is also possible to reverse the direction of iteration and extend to all negative indices n (again, this is always possible subject to the condition that d n does not vanish for some n). This immediately leads to a J-fraction expression for Y 0 , that is (4.24) which corresponds to a power series expansion around ∞ 2 , This means that the quantities d n , v n can also be written in terms of ratios of determinants when n is negative, but involving the Hankel determinant as well as the associated shifted Hankel determinant ∆ † * n , which is just the analogue of (4.13) built from the coefficients in the series (4.25).
Proof: Since G † = Y 0 = (Y + P 0 )/Q 0 , the generating function for the moments s † j satisfies the quadratic and subsequent coefficients in the series expansion (4.3) are determined for j ≥ 3 by In particular, taking the specific curve Y 2 = (X 3 − 5X − 1) 2 − 4(X 2 + 2X + 3) that was used for illustration in Example 4.2, with the same function Y 0 , as before we have f = −5, g = u = −1, u 0 = −4, d 0 = 5/4, e 0 = 3/5, v 0 = −1/2, w 0 = −3/2, and also v −1 = −1/10, w −1 = −3/2, which gives with X → ∞ and Y ∼ −X 3 , and the recursion for the coefficients (moments) in the above expansion is with initial values s † 0 = −5/8, s † 1 = −1/16, s † 2 = 45/32. In this case, the sequence of moments (s † j ) begins  for all n ∈ Z, for some set of quantities τ n , τ * n . However, in general one cannot just take τ n = ∆ n for non-negative n and τ n = ∆ † −n−1 for negative n (and similarly for τ * n ) , because there will be a mismatch at the values of d 0 , d 1 and v 0 which are left unspecified by Theorems 4.1 and 4.2. Nevertheless, one can make use of the fact that the expressions for d n and v n in (4.30) are left invariant by the three-parameter group of gauge transformations τ n → ab n τ n , τ * n → ab n (τ * n + cτ n ).

The Somos connection
In this section, we explain how Somos sequences naturally arise from the continued fraction expansion, as quadratic relations for the Hankel determinants ∆ n . This is most straightforward to describe in the genus one case, as it follows from the fact that, for a fixed value of the first integral H = −uv given by (3.17), each orbit of (3.16) coincides with an orbit of a symmetric QRT map, and, as was already noted in [47] in an example related to the hard hexagon model, the bilinear form of the latter is precisely (1.1). (For a detailed discussion of normal forms of QRT maps restricted to fixed invariant curves, see [29,30].) Proof: Putting X = v n into (3.15) and then applying the first equation in (3.14) yields Then putting X = v into (3.15) and using the above result gives so that (5.1), which is an example of a symmetric QRT map [47], follows immediately.
The connection with Somos-4 is almost immediate, since if d n = τ n τ n−2 /(τ n−1 ) 2 , then τ n satisfies (1.1) whenever d n satisfies (5.1). So in particular, by Theorem 4.1, the Hankel determinants for g = 1 satisfy a Somos-4 relation, and since the latter is invariant under the first of the gauge transformations (4.31), any Somos-4 sequence can be expressed in terms of Hankel determinants. More precisely, starting from any Somos-4 sequence, one can always make a gauge transformation to a sequence with τ 0 = 1, and then use the coefficients α, β and the other initial conditions τ 1 , τ 2 , τ 3 to specify the values of s 0 = ∆ 1 = τ 1 and v 0 , v 1 , d 1 , f , so that the values of s 1 and the coefficients in (4.18) are fixed. (In fact, since the gauge transformation involves two parameters a, b, there is also the freedom to fix τ 1 = 1, which corresponds to taking s 0 = 1.) Thus we arrive at the following result.
Theorem 5.2. In the case g = 1, the Hankel determinants (4.12) with moments defined recursively by (4.18) satisfy the Somos-4 recurrence Moreover, every solution of the Somos-4 recurrence (1.1) can be written in the form where ∆ † n is constructed from moments that satisfy (4.27) with g = 1, for suitable constantsâ,b, a † , b † . There is an apparent mismatch between the Hankel determinants in (1.7), which were shown by Xin to yield the terms of the Somos-4 sequence (1.2), and those in (4.21) above. We now explain the relation between these two sets of Hankel determinant formulae, and see how the results of [11] are a consequence of the continued fraction expansion for g = 1.
Theorem 5.3. For n ≥ 2, the quantity d n that satisfies (3.16) is given by in terms of the Hankel determinant (1.4) defined in terms of momentss j that satisfy the recursion (1.8) for j ≥ 2, withs Moreover, the sequence (D n ) n≥0 is identical to the sequence of Hankel determinants (∆ n ) n≥0 with moments satisfying (4.18), hence satisfies the Somos-4 recurrence (5.2) with coefficients as in (5.3).
Proof: By replacing X with the shifted variableX = X − v 0 , and lettingG(X) = G(X), we obtain the J-fractionG where, from (4.17) in the case g = 1, the generating functionG satisfies In the first line of the continued fraction (5.5), we are at liberty to chooseũ 0 = u 0 , which implies that d 1 = d 1 , and then in each subsequent line we haveṽ n = v n + v 0 for n ≥ 1, and alsod n = d n for n ≥ 2. Hence, by the same argument as in the proof of Theorem 4.1, the formula (5.4) holds for n ≥ 2. The moments s j are obtained from the series expansion ofG in powers ofX, with the leading order term (orderX) in (5.6) givings 0 = s 0 = u 1 /2, and the next to leading order term (orderX 0 ) givings 1 + 2v 0s0 = u1 , we find the recursion relation (1.8) with the stated values ofα,β,γ. By convention we have D 0 = ∆ 0 = 1, and also D 1 =s 0 = s 0 = ∆ 1 , and then it follows by induction from (5.4) that D n = ∆ n for all n ≥ 0.
Remark 5.1. In order to see how Xin's result [56] follows from the above, it is sufficient to note that the quartic curve (4.19) is isomorphic to (1.5), via the birational equivalence Y = x −2 (1 − 2η), X =X − 1 = x −1 − 1, so that the expansion (1.6) in powers of x is equivalent to an expansion in powers ofX −1 . Also, by setting η = (1 − y)/2, the Weierstrass cubic y 2 = 4x 3 − 4x + 1 derived from analytic formulae in [23] is seen to be isomorphic to the curve (1.5); over Q, this is known as 37a1, the elliptic curve of minimal conductor with positive rank. (See www.lmfdb.org/EllipticCurve/Q/37/a/1 in the online database of L-functions, modular forms, and related objects.) The analogue of Theorem 5.2 in genus two is more difficult to state explicitly, due to the size of the expressions for the coefficients, and at present we are only able to prove it with the use of computer algebra.
Theorem 5.4. In the case g = 2, the Hankel determinants (4.12) with moments defined recursively by (4.22) satisfy a Somos-8 recurrence of general type, that is where the coefficients α 1 , . . . , α 5 are certain first integrals of the map (3.29).
Proof: The recurrence (5.7) is equivalent to a relation for the iterates d n , that is along an orbit of the 4D map (3.29). Equivalently, writing a solution of this map in the form (4.30), it means that τ n should satisfy the Somos-8 recurrence (5.7) with some coefficients α j that are constant on each orbit. This requires the vanishing of a 5 × 5 determinant, namely τ n+4 τ n−4 τ n+3 τ n−3 τ n+2 τ n−2 τ n+1 τ n−1 τ 2 n τ n+5 τ n−3 · · · · · · · · · τ 2 n+1 . . . . . . τ n+8 τ n · · · · · · · · · τ 2 n+4 = 0 (5.8) for all n, and also that the ratios of certain 4 × 4 minors should be independent of n. In particular, denoting byD j,n the minor formed from the first four rows in (5.8) with the jth column removed, so that , and so on, the existence of the relation (5.7) for τ n is equivalent to the requirement that the ratios α j α 1 =D j,n /D 1,n , j = 2, 3, 4, 5 are independent of n (noting the possibility of a vanishing denominator, in the case that α 1 = 0). So in order to verify the statement, it is sufficient to check that D j,nD1,n+1 −D j,n+1D1,n = 0, (5.9) holds for each j = 2, 3, 4, 5, and for all n. In fact, for each j, it is enough to check that this holds for a single shift n → n + 1 with arbitrary initial data in the map (3.29), which guarantees that each of the ratios α j /α 1 is a first integral, and then it automatically holds for all n. Even with the help of computer algebra, this is not a completely straightforward task, and in order to do it as efficiently as possible it is convenient to make a gauge transformation (4.31) to fix s 0 = 1, and then note that there is a one-to-one correspondence between two sets of seven parameters: the four initial values d 0 , d 1 , v 0 , v 1 and three parameters f, g, u needed to iterate the map (3.29), and the two initial values s 1 , s 2 and five coefficientsα,β,γ,δ,˜ that specify the genus two recursion (4.22) in the form Now, using the above recursion, one can calculate the sequence (s j ) j≥0 , and then compute the Hankel determinants τ n = ∆ n for n = 0, 1, 2, . . ., which are polynomials in Z[s 1 , s 2 ,α,β,γ,δ,˜ ], but this rapidly becomes very computationally intensive as n increases. More efficient is to rewrite the map (3.29) as an equivalent pair of coupled recurrence relations of degree 6 for τ n , τ * n , which are of overall order 7. To iterate the latter, one needs 7 initial values (four adjacent τ n and three adjacent τ * n ), and it is convenient to take τ 0 = ∆ 1 = 1, τ 1 = ∆ 0 = s 0 = 1 but also τ −1 = d 1 =γ and τ −2 =γ 2 d 0 =αγ +γ 3 −δ 2 +γ˜ , together with The verification of (5.9) requires 13 adjacent values of τ n , but due to the size of the expressions involved it is best to compute only up to τ 6 using five forward steps of the coupled recurrence for τ n , τ * n , and then apply this recurrence in reverse, making four backward steps to go back as far as τ −6 , so that the adjacent values τ −6 , τ −5 , . . . , τ 5 , τ 6 are obtained as explicit polynomials in Z[s 1 , s 2 ,α,β,γ,δ,˜ ]. This means that the minorsD j,−2 andD j,−1 can be computed explicitly, which allows (5.9) to be checked directly when n = −2. The formulae for the first integrals α j /α 1 , j = 2, 3, 4, 5 as rational functions of s 1 , s 2 ,α,β,γ,δ,˜ are so large that they are difficult to display even on a computer screen, but if α 1 is regarded as the first integral that is the lowest common denominator of these four quantities, then we arrive at the Somos-8 relation (5.7).
Remark 5.2. As already noted, there is the possibility of a vanishing denominator α 1 = 0 in the ratios α j /α 1 , j = 2, 3, 4, 5. Given that the map (3.29) only has two independent first integrals, which can be specified by H 1 = −uv, H 2 = uw as in (3.26) and (3.27), it follows that these four ratios are rational functions of f, g, u, H 1 , H 2 , with coefficients in Q, so that α 1 can be fixed as the polynomial in Z[f, g, u, H 1 , H 2 ] that is the lowest common denominator of these four rational functions. Thus it can happen that α 1 = 0 for certain combinations of f, g, u, H 1 , H 2 , in which case ∆ n (or τ n ) satisfies a Somos-6 relation, rather than a Somos-8. Numerical experiments suggest that u|α 1 (f, g, u, H 1 , H 2 ), consistent with a result of van der Poorten, who showed that there is a Somos-6 relation in the special case u = 0 [45,46].
Before concluding this section, we state a conjecture which is the genus two analogue of Theorem 1.1.
Conjecture 5.5. When g = 2, the Hankel determinants (4.12) with moments defined by (4.22) are given by 10) in terms of the genus two Kleinian sigma function σ(z) = σ(z;c 0 ,c 1 ,c 2 ,c 3 ) associated with a quintic curveC that is isomorphic to the sextic C in (3.19), given byC : where ∞ is the unique point at infinity onC,P 2 ∈C is the point corresponding to ∞ 2 ∈ C, andâ,b are certain non-zero constants.
A proof of the above result would follow from an analytic solution for the iterates of the map (3.29), which we propose to consider elsewhere. However, to see why this result is plausible we note that since the class of the divisorP 2 −P 1 ∼ lin 2(P 2 − ∞) in the Jacobian ofC corresponds to that of the divisor ∞ 2 − ∞ 1 on C, each shift n → n + 1 increases the argument of the numerator in (5.10) by z, which is consistent with (2.18). Furthermore, if the formula (5.10) is correct then, by essentially the same analytical calculations as those in [7,26], it follows immediately that ∆ n satisfies the Somos-8 relation (5.7), or a Somos-6 relation when a certain constraint on z holds.
In the higher genus case, we further conjecture that there should be an analytical formula analogous to (5.10). Equivalently, there should be an expression in terms of the Riemann theta function associated with the Jacobian of the curve (2.1), of the form ∆ n =âb nĉn 2 Θ(v 0 + nv) for someâ,b,ĉ ∈ C * and v 0 , v ∈ C g . If this expression holds, then by counting the dimension of the vector space of quasiperiodic functions of weight 2 with respect to the period lattice (see [39]), it follows that ∆ n satisfies a Somos-k relation for some k ≤ 2 g+1 . We have verified this in numerical examples for g = 3, 4.

Poisson structure and integrability
In this section we slightly change our notation, and consider a modified family of Lax matrices, given by where P and R are monic polynomials of degrees g + 1 and g in ζ, respectively, and Q is a polynomial of degree g with non-constant leading coefficient, which we write as The set of all such matrices forms an affine space of dimension 3g + 2, with coordinates given by the nontrivial coefficients of P, Q and R. We will endow this space with a particular Poisson structure of rank 2g, and show how this leads to the construction of an associated set of Hamiltonian vector fields that are completely integrable in the Liouville sense. Then we will present a compatible discrete integrable system on the same phase space, and show that this is equivalent to the iteration of the recursion for the continued fraction expansion of the hyperelliptic function considered previously.
The Poisson brackets between the entries of L are specified by In terms of the coefficients appearing in L, this is a linear bracket, since the right-hand sides are linear in P, Q, R. In order for this to define a Poisson bracket, it must satisfy the Jacobi identity, and although this can be verified directly this requires many tedious calculations; we set this question aside for now, and a simpler argument will be presented in due course.
To begin with, we consider the function Using the above bracket relations, we find that It is straightforward to check that the right-hand sides of the above expressions are polynomials in η of degrees g − 2, g − 1, g − 1 respectively. Thus, if we expand then these expressions imply that { · , c j } = 0 for j = g, g + 1, . . . , 2g + 1, or in other words the leading g + 2 non-trivial coefficients of F are Casimirs. In fact c g , . . . , c 2g+1 provide the full set of Casimirs, and the symplectic leaves have dimension 2g. To see this, factorize R as and set y i = P(x i ), i = 1, . . . , g. (6.11) Then P(ζ) = ζ g+1 + 1 2 where, from (6.11), the g coefficients π j can be found explicitly as functions of the x i , the y i and the Casimir c 2g+1 by solving a linear system. The two relations in (6.2) then imply that for all i, j, while from the second bracket in (6.3) it follows that so that (up to scaling) the pairs (x i , y i ) i=1,...,g provide a set of canonical coordinates on a symplectic manifold of dimension 2g, and by (6.5), for fixed values of the coefficients c j in (6.9), each pair (ζ, µ) = (x i , y i ) is a point on the spectral curve µ 2 = F (ζ). (6.12) Then by using (6.10), (6.11) and the leading terms in (6.5) up to and including order ζ g , all of the coefficients in L can be expressed as functions of the x i , y i and the Casimirs, so the whole Poisson algebra is expressed in terms of these coordinates. The latter argument begs the question of whether the Poisson brackets for the entries of L satisfy the Jacobi identity in the first place, but this is easily seen by reversing the direction of the preceding argument. Indeed, starting with the canonical bracket between the x i and y j , one extends it with the Casimirs c g , . . . , c 2g+1 to obtain a Poisson algebra of dimension 3g + 2, where the entries of L are defined in terms of these canonical coordinates as above, and by construction they satisfy the linear bracket relations given before. Hence the Jacobi identity is trivially satisfied.
We now consider a family of vector fields on the space of Lax matrices, defined by the flow From the bracket relations (6.6), (6.7), (6.8), it can be verified directly that this can be written in the form of a Lax equation, that is d dt with the matrix Moreover, all the flows in this family commute with one another, because the same bracket relations imply that Thus it follows that all of the coefficients c j in (6.9) are in involution, and in order to get non-trivial flows we can take the non-Casimir functions (Hamiltonians) Then since { H j , H k } = 0 for all j, k, this gives g commuting flows which can be written in Lax form as where, for j = 0, . . . , g − 1, the matrix P (j) is defined from (6.14) by P (j) (ζ) = res P(ζ, η) η j+1 η=0 .
Hence we have integrability in the sense of Liouville [1].
Theorem 6.1. The Hamiltonians (6.16) define a completely integrable system on the space of Lax matrices (6.1).
In fact, there is more that one can say: for fixed values of c j in (6.9), the set of triples of polynomials P(ζ), Q(ζ), R(ζ) of the specified form that satisfy (6.5) for a fixed set of coefficients c j in (6.9) is an affine algebraic variety of dimension g, that is canonically isomorphic to the affine Jacobian of the corresponding hyperelliptic spectral curve (6.12), by associating each such triple with the degree g divisor D = (x 1 , y 1 ) + · · · + (x g , y g ) defined by (6.10) and (6.11). As is explained in [40], the analogous construction of the Jacobian variety in the case of odd hyperelliptic curves goes back to Jacobi, and arises in the context of finite gap solutions of the Korteweg-deVries equation. The Poisson brackets and first integrals H j are all algebraic -in actual fact, they are given by polynomial functions of the coefficients of the polynomials P(ζ), Q(ζ), R(ζ) -and over C the generic common level set of these first integrals is an affine part of a complex algebraic torus, with the Hamiltonian flows being linear on the torus, so this is what is known as an algebraic completely integrable system (see [52] and references for details).
We now proceed to describe how the Lax matrices (6.1) are related to the discrete Lax pair and nonlinear map introduced in section 3.
First of all, observe that completing the square in (6.9) means that F = A 2 +c g ζ g +c g−1 ζ g−1 · · · +c 0 , A(ζ) = ζ g+1 +c 2g+1 ζ g + · · · +c g+1 , wherec 2g+1 = 1 2 c 2g+1 and all the coefficientsc 2g+2−j for j = 1, . . . , g + 2 are polynomial functions of the g + 2 Casimirsc 2g+2−j for the same range of j, and there is a bijection between these two sets of Casimir functions. In particular, all of the coefficients of A are Casimirs, and from (6.5) we may write and take the g non-trivial coefficients of R(ζ) and the quantities π (j) 0 for j = 0, . . . , g − 1 as coordinates on each symplectic leaf. We shall see shortly that the latter is consistent with the notation in (3.7), but before getting to this we must restrict to a particular set of symplectic leaves, by fixing the value of the top Casimir to be zero, i.e. c 2g+1 = 0, so thatc 2g+1 = 0, and the coefficients of F , A and P at next to leading order are zero; this slightly simplifies the formulae for the nonlinear map, and agrees with our previous conventions for the continued fraction expansion (but if necessary the case of non-zero c 2g+1 can always be obtained by making a shift in the spectral parameter ζ).
Next, in order to obtain (6.1), we wish to remove the multipliers u n that appear in the off-diagonal terms of (3.2), since although they provide an arbitrary choice of scale in the continued fraction, they do not behave well from the point of view of the Poisson structure. Thus we consider diagonal gauge transformations applied to the eigenvector in (3.1). These have the effect of changing the prefactors in the off-diagonal entries of L n while leaving the diagonal terms the same. Hence, for a suitable choice of λ (1) n , λ n , upon setting n = 0 we obtain L(ζ) = G 0 L 0 (ζ)G −1 0 with L(ζ) being given by (6.1), and from (3.3) we find Comparing the effect of the gauge transformation with the notation used in section 3, it is apparent that P(ζ) = P 0 (ζ), Q(ζ) = u −1 Q 0 (ζ) and R(ζ) = Q −1 (ζ)/u −1 . from which it follows that { v 0 ,R(ζ) } = 0, (6.23) and, making use of (6.21), we also have { P(ζ),R(η) } = −4d 0 { P(ζ), Q(η) } = 2 R (ζ) −R(η) ζ − η , (6.24) which then shows that { v 0 , P(ζ) } = res 2(R(ζ) −R(η)) η g (ζ − η) ζ − η , (6.26) corresponding to the shifted version of the second bracket relation in (6.3). It is not necessary to verify directly that the remaining three bracket relations are preserved by the map, since they follow from observing that, given the g pairs of coordinates (x i ,ỹ i ) defined bỹ R(ζ) = g j=1 (ζ −x j ),ỹ i =P(x i ), for i = 1, . . . , g, the already verified relations (6.20) and (6.26) imply that Hence the map defined by (6.19) restricts to a canonical transformation on each symplectic leaf, and it preserves all the Casimirs, so it is a Poisson map. This also proves Theorem 3.1, since π (j) 0 together with ρ (j) −1 (the coefficients of R) for j = 1, . . . , g also provide coordinates on each symplectic leaf, so that the map (3.12) which is written in these coordinates is symplectic and has g first integrals in involution, as in (6.16). Thusφ is integrable, as is the conjugate map ϕ given by (3.11).

Appendix: Identities for determinants of Hankel type
There are various ways to derive the classical formulae (4.14), and the expression for d n in particular; see the proof of Theorem A in [54], for instance. However, here we present determinantal formulae that yield the latter expression directly from the three-term relation (4.9).