2. Analytic formulae

To begin with, we paraphrase the main result of [Reference Hone10] and briefly describe it in a way that incorporates certain later observations made in [Reference Hone11] and [Reference Hone and Swart14].

Theorem 2.1. The general solution of the initial value problem for (1.1) over ${\mathbb{C}}$ is

(2.1) \begin{equation}x_n = A\, B^n \, \frac{{\sigma} (z_0+nz)}{{\sigma}(z)^{n^2}},\end{equation}

where ${\sigma}(z)={\sigma}(z;g_2,g_3)$ is the Weierstrass sigma function for an associated elliptic curve:

(2.2) \begin{equation}{y}^2=4{x}^3-g_2{x}-g_3,\end{equation}

where $A,B\in{\mathbb{C}}^*$ and $z,z_0,g_2,g_3\in{\mathbb{C}}$ are explicitly determined by four nonzero initial values and coefficients $x_0,x_1,x_2,x_3,{\alpha},{\beta}$ .

More precisely, the solution of the initial value problem for (1.1) is achieved by considering the sequence of ratios:

(2.3) \begin{equation}d_n = \frac{x_{n+1}x_{n-1}}{x_n^2},\end{equation}

which satisfy a recurrence of second order, namely

(2.4) \begin{equation}d_{n+1}d_{n-1} = \frac{{\alpha} d_n+{\beta}}{d_n^2}.\end{equation}

The latter recurrence is equivalent to a birational map of the plane, being an example of a symmetric QRT map [Reference Quispel, Roberts and Thompson24], and it has a rational conserved quantity:

(2.5) \begin{equation} J=d_{n}d_{n-1}+ {\alpha} \left(\frac{1}{d_n}+\frac{1}{d_{n-1}}\right) + \frac{{\beta}}{d_nd_{n-1}}\end{equation}

that defines a pencil of biquadratic plane curves, each of which is generically of genus 1 (except for a finite number of singular fibres, for certain special values of J [Reference Duistermaat3]) and hence isomorphic to a cubic curve of the form (2.2). Explicit computations with Weierstrass functions then show that the formula (2.1) does indeed provide a solution of the Somos-4 recurrence (1.1), with the coefficients being given by:

(2.6) \begin{equation}{\alpha} =\frac{{\sigma} (2z)^2}{{\sigma} (z)^8}, \qquad {\beta} = -\frac{{\sigma}(3z)}{{\sigma}(z)^9},\end{equation}

allowing $g_2,g_3$ to be obtained in terms of ${\alpha},{\beta},J$ via the expressions:

(2.7) \begin{equation}g_2 =12{\lambda}^2-2J, \quad g_3=4{\lambda}^3-g_2{\lambda}-{\alpha}, \quad {\lambda}=\frac{1}{3{\alpha}}\left(\frac{J^2}{4}-{\beta}\right),\end{equation}

and then (up to fixing an overall sign) z and $z_0$ are determined modulo periods by the elliptic integrals:

(2.8) \begin{equation}z=\int^{\lambda}_\infty \frac{{d}{x}}{{y}}, \qquad z_0=\int^{{\lambda}-d_0}_\infty \frac{{d}{x}}{{y}},\end{equation}

so that once these are given the values of A,B can be found from the initial data. (Note that when $g_2^3-27g_3^2=0$ , corresponding to a singular fibre, the formula (2.1) is still valid in terms of suitable trigonometric/rational limits of the sigma function.)

For the purposes of what follows, it is important to note that from (2.3) the expression (2.5) is equivalent to a rational conserved quantity for the Somos-4 recurrence (1.1) itself, that is,

(2.9) \begin{equation}J = \frac{x_n^2x_{n+3}^2 +{\alpha}\, \left(x_{n+1}^3x_{n+3}+x_nx_{n+2}^3\right) +{\beta}\,x_{n+1}^2x_{n+2}^2 }{x_nx_{n+1}x_{n+2}x_{n+3}}.\end{equation}

An elementary observation, which is the basis for the application of dual numbers to automatic differentiation, is that for any differentiable function $\Phi$ and $X=x+y\varepsilon\in{\mathbb{D}}$ , the value of $\Phi$ at X is

\begin{align*}\Phi(X)=\Phi(x)+\Phi^{\prime}(x)y\varepsilon.\end{align*}

We can use this result to describe solutions of (1.7) in analytic form.

Proposition 2.2. Given parameters $A=A^{(0)}+A^{(1)}\varepsilon$ , $B=B^{(0)}+B^{(1)}\varepsilon$ , $G_2=g_2+\tilde{g}_2\varepsilon$ , $G_3=g_3+\tilde{g}_3\varepsilon$ , $Z_0=z_0+\tilde{z}_0\varepsilon$ , $Z=z+\tilde{z}\varepsilon\in{\mathbb{D}}$ , the expression

(2.10) \begin{equation}X_n = A\, B^n \, \frac{{\sigma} (Z_0+nZ)}{{\sigma}(Z)^{n^2}}\end{equation}

satisfies the dual Somos-4 recurrence (1.7) with coefficients:

\begin{align*}{\alpha} ={\alpha}^{(0)}+ {\alpha}^{(1)}\varepsilon=\frac{{\sigma} (2Z)^2}{{\sigma} (Z)^8}, \quad{\beta} ={\beta}^{(0)}+ {\beta}^{(1)}\varepsilon = -\frac{{\sigma}(3Z)}{{\sigma}(Z)^9},\end{align*}

provided that $A,B\in{\mathbb{D}}^*$ and ${\sigma}(Z)={\sigma}(Z;G_2,G_3)\in{\mathbb{D}}^*$ . In terms of even/odd components, this can be written as:

(2.11) \begin{equation}X_n =x_n +x_n\left(\frac{A^{(1)}}{A^{(0)}}+ \frac{B^{(1)}}{B^{(0)}} \, n+\tilde{z}_0\zeta(z_0+nz) + \left(\tilde{z}\partial_z+\tilde{g}_2\partial_{g_2}+\tilde{g}_3\partial_{g_3}\right)\log x_n\right)\,\varepsilon,\end{equation}

where $x_n$ denotes the right-hand side of (2.1) with the replacement $A\to A^{(0)}$ , $B\to B^{(0)}$ , and the other parameters $z,z_0,g_2,g_3$ being the same, $\zeta(z)=\zeta(z;g_2,g_3)$ is the Weierstrass zeta function, and $\partial$ denotes a partial derivative.

Proof. As pointed out in [Reference Hone11], the fact that the analytic expression (2.1) satisfies a Somos-4 relation with coefficients given by (2.6) is a direct consequence of the three-term identity for the sigma function, which can be viewed as an algebraic identity between power series. The sigma function ${\sigma}(z;g_2,g_3)$ is a holomorphic function of the argument z and is also holomorphic in the two invariants $g_2,g_3$ . Therefore, the same proof carries over to the commutative algebra ${\mathbb{D}}$ , with the requirement that ${\sigma}(Z)\in{\mathbb{D}}^*$ so that one can divide by suitable powers of this quantity, and $A,B\in{\mathbb{D}}^*$ so that the even part $x_n\neq 0$ . The solution depends analytically on six dual parameters, so it can be expanded as $X_n=x_n+y_n\varepsilon$ by differentiating with respect to each of these parameters, giving

\begin{align*}y_n =(A^{(1)}\partial_{A^{(0)}}+ B^{(1)}\partial_{B^{(0)}}++\tilde{z}_0\partial_{z_0}+\tilde{z}\partial_{z}+\tilde{g}_2\partial_{g_2}+\tilde{g}_3\partial_{g_3})\,x_n.\end{align*}

It is convenient to rewrite this in terms of logarithmic derivatives of $x_n$ , so that $y_n$ is written as $x_n$ times a sum of such derivatives, and the first three derivatives are straightforward to calculate explicitly, as in (2.11).

We believe that the analytic formula (2.10) represents the general solution of the dual Somos-4 recurrence (1.7). However, a complete proof would require showing that, given ${\alpha},{\beta}\in{\mathbb{D}}$ , one can solve the initial value problem for (1.7) by determining the four dual parameters $G_2,G_3,Z_0,Z$ . The difficulty lies in constructing a suitable theory of elliptic functions and uniformization of elliptic curves over the dual numbers, which in particular would provide analogues of the elliptic integrals (2.8). An abstract theory of super Riemann surfaces is available [Reference Witten28], but to the best of our knowledge the genus 1 case has not been developed explicitly in terms of Weierstrass functions with odd components. Rather than trying to pursue this further here, we will take the formula (2.11) as the starting point for a more straightforward algebraic approach in the next section.

3. Solution of linear difference equation

Since the system consisting of (1.8) and (1.9) is triangular, and we already have the general solution of (1.8), given by Theorem 2.1, the main outstanding problem is to understand the solutions of (1.9). We beginning by considering the special case ${\alpha}^{(1)} ={\beta}^{(1)}=0$ , which is just a homogeneous linear equation for $y_n$ , namely

(3.1) \begin{equation}x_n y_{n+4} - {\alpha}^{(0)} x_{n+1}y_{n+3}-2{\beta}^{(0)}x_{n+2}y_{n+2} - {\alpha}^{(0)} x_{n+3}y_{n+1} +x_{n+4}y_n = 0. \end{equation}

The solutions of the latter equation are the shadow Somos-4 sequences, in the sense of [Reference Ovsienko21].

However, the equation (3.1) is just the linearization of the Somos-4 recurrence (1.8), and we have the explicit solution of this, which (for fixed coefficients ${\alpha}^{(0)}$ , ${\beta}^{(0)}$ ) depends on four arbitrary parameters; these can be taken to be $A^{(0)}$ , $B^{(0)}$ , $z_0$ , and $J^{(0)}$ , where the last is the value of the first integral (2.9). Thus, by adapting a technique that goes back to Lie’s work on ordinary differential equations, we can observe that the derivative of the solution $x_n$ with respect to each of these parameters provides a solution of the linearized equation, and thus we obtain four linearly independent solutions, which span the whole vector space of solutions of (3.1).

Lemma 3.1. The Somos-4 shadow equation (1.9) has four linearly independent solutions, given by $y_n^{(i)}=x_n$ , $y_n^{(ii)}=nx_n$ , $y_n^{(iii)}=x_n\zeta(z_0+nz)$ and

(3.2) \begin{equation}y_n^{(iv)}=x_n\partial_{J^{(0)}}\log x_n=x_n\left(\frac{{d} z}{{d} J^{(0)}}\partial_z+\frac{{d} g_2}{{d} J^{(0)}}\partial_{g_2}+ +\frac{{d} g_3}{{d} J^{(0)}}\partial_{g_3} \right)\log x_n .\end{equation}

Proof. The derivative of the formula (2.1) with respect to the scaling parameter A is just proportional to $x_n$ , while the derivative with respect to B is proportional to $nx_n$ . Differentiating the numerator in (2.1) with respect to $z_0$ produces $x_n$ times the Weierstrass zeta function, but the dependence on $J^{(0)}$ is more complicated, as both z and the invariants $g_2$ and $g_3$ vary with this parameter.

The parameters A,B in (2.1) correspond to the action of a two-parameter scaling group, $x_n\to AB^n\,x_n$ for $(A,B)\in({\mathbb{C}}^*)^2$ , which leaves the Somos-4 recurrence invariant. The freedom to vary $z_0$ corresponds to an arbitrary choice of initial point on the elliptic curve (2.2), viewed as a complex torus ${\mathbb{C}}/\Lambda$ where $\Lambda$ is the period lattice. The first three terms in the odd part of (2.11) are clearly proportional to $y_n^{(i)}$ , $y_n^{(ii)}$ , and $y_n^{(iii)}$ . We will shortly present an algebraic method to calculate the third independent shadow solution $y_n^{(iii)}$ .

As for the fourth solution $y_n^{(iv)}$ , it is much more complicated, because it involves the modular derivatives $\partial_{g_2} $ and $\partial_{g_3}$ . Indeed, if we fix the coefficients ${\alpha}^{(0)}$ and ${\beta}^{(0)}$ , then we can consider the following system of three equations in terms of the Weierstrass $\wp$ function and its derivatives:

(3.3) \begin{equation}{\alpha}^{(0)}=\wp^{\prime}(z)^2, \qquad \frac{{\beta}^{(0)}}{{\alpha}^{(0)}}=\wp(2z)-\wp(z),\qquad J^{(0)}=\wp^{\prime\prime}(z);\end{equation}

the first two equations are just equivalent to (2.6), corresponding to identities between elliptic functions (see [Reference Hone10, Reference Hone11]). Differentiating each of these equations with respect to $J^{(0)}$ and applying the chain rule gives a system of three linear equations for the derivatives of $z,g_2,g_3$ with respect to $J^{(0)}$ , appearing in (3.2). As for the modular derivatives of the sigma function, note that it is a weighted homogeneous function of $z,g_2,g_3$ , so it satisfies the linear partial differential equation (PDE):

\begin{align*}\left(4g_2\partial_{g_2}+6g_3\partial_{g_3} -z\partial_z+1\right)\sigma\left(z;g_2,g_3\right)=0.\end{align*}

As was shown by Weierstrass [Reference Weierstrass27], it satisfies another linear PDE which is first order in $g_2,g_3$ and second order in z: this is equivalent to the well-known heat equation for the Jacobi theta function. (It also satisfies a fourth-order nonlinear ODE in z that can be written in Hirota bilinear form [Reference Eilbeck and Enolskii4], which can be viewed as the continuous analogue of Somos-4.) In principle, one could combine all these facts to rewrite (3.2) in terms of z derivatives only. However, we do not need to calculate these analytic expressions in more detail, because in due course we will present an algebraic way to characterize a fourth independent shadow solution, which is much more straightforward.

In order to derive the third shadow solution in a more algebraic way, we consider the following coupled pair of recurrence relations with parameters u,f, which was considered in [Reference Hone12] as the simplest example of a map arising from Jacobi continued fractions in hyperelliptic function fields, based on a construction introduced by van der Poorten [Reference van der Poorten26]:

(3.4) \begin{align}v_n & = -v_{n-1}+u/d_n, \nonumber\\d_{n+1} & = -d_n-v_n^2-f.\end{align}

The above system defines an integrable map, and here we will present its analytic solution.

Proposition 3.3. The map $(v_{n-1},d_n)\mapsto (v_n,d_{n+1})$ defined by (3.4) is integrable: it preserves the symplectic form ${d} v_{n-1}\wedge {d} d_n$ , and has the first integral:

\begin{align*}H=d_n\left(v_{n-1}^2+d_n+f\right)-uv_{n-1}.\end{align*}

The general solution of the map can be written in terms of elliptic functions, up to an overall choice of sign as:

(3.5) \begin{equation}d_n=\wp (z)-\wp(z_0+nz), \qquad v_n=\pm \big({\zeta} (z_0+(n+1)z)- {\zeta} (z_0+nz)-{\zeta}(z)\big),\end{equation}

with the parameters given by:

\begin{align*}u=\pm \wp^{\prime}(z), \qquad f=-3\wp(z),\end{align*}

on the level set $H=-\frac{1}{2}\wp^{\prime\prime}(z)$ .

Proof. It is straightforward to verify from the map (3.4) that ${\omega} = {d} v_{n-1}\wedge {d} d_n={d} v_{n}\wedge {d} d_{n+1}$ , so it is a symplectic map, and a short calculation also shows that H is independent of n, so this is a discrete integrable system with one degree of freedom. By Proposition 5.1 of [Reference Hone12], on each level set $H=\,$ const, the sequence of values of $d_n$ coincides with an orbit of (2.4), if the coefficients and the value of the first integral (2.5) are identified as:

(3.6) \begin{equation}{\alpha}=u^2, \qquad {\beta}=u^2(v^2+f), \qquad J=2uv=-2H.\end{equation}

As for the analytic solution of the map, it was shown in [Reference Hone10] that when $x_n$ is given by (2.1), the corresponding solution of (2.4), associated with it via (2.3), is given in terms of the Weierstrass $\wp$ function by the first formula in (3.5), and when this expression for $d_n$ is substituted into the second component of (3.4), the explicit formula for $v_n$ is found from taking the square root in the left-hand side of the elliptic function identity:

\begin{align*}\big({\zeta}(a)+{\zeta}(b)+{\zeta}(c)\big)^2=\wp(a)+\wp(b)+\wp(c),\quad\mathrm{for}\quad a+b+c\equiv 0 \bmod{\Lambda}.\end{align*}

Then given the two formulae in (3.5), another identity between elliptic functions verifies that the first component of the system (3.4) is satisfied for all $n\in{\mathbb{Z}}$ , with the parameter $u=\pm \sqrt{{\alpha}}=\pm \wp^{\prime}(z)$ , making the same choice of sign as for $v_n$ .

Remark 3.4. Given the canonical Poisson bracket $\{\,d_n,v_{n-1}\,\}=1$ associated with ${\omega}$ , the first integral H defines the Hamiltonian vector field $\{ \,\cdot, H\,\}$ , whose flow commutes with the shift $n\to n+1$ corresponding to the map. A direct calculation shows that, for a fixed choice of scale, this flow implies that the sequence of $d_n$ and $v_n$ satisfy the set of differential equations:

(3.7) \begin{align}\dot{v}_n & = d_n-d_{n+1}, \nonumber\\\dot{d}_n & = d_n(v_{n-1}-v_n),\end{align}

for all $n\in{\mathbb{Z}}$ , which (up to squaring one of the variables) is just the infinite Toda lattice written in Flaschka coordinates. Thus, the simultaneous solutions of the map and the Hamiltonian flow provide genus 1 solutions of the Toda lattice, and for genus $g>1$ an analogous statement holds for the other maps constructed in [Reference Hone12]; further details will be presented elsewhere.

Corollary 3.5. Up to subracting off multiples of $y_n^{(i)}=x_n$ and $y_n^{(ii)}=nx_n$ and overall scale, the third independent shadow solution of Somos-4 can be obtained from the associated solution of the map (3.4), in the form:

(3.8) \begin{equation}y_n^{(iii)}=-x_n\sum_{j=0}^{n-1}v_j, \qquad n\geq 0\end{equation}

(where the case $n=0$ corresponds to an empty sum).

Proof. If we take the plus sign in the formula for $v_n$ in (3.5), then we have a telescopic sum:

\begin{align*}-x_n\sum_{j=0}^{n-1}v_j =x_n\big(n{\zeta}(z)+{\zeta}(z_0)-{\zeta}(z_0+nz)\big)\end{align*}

for $n\geq 0$ , and after subtracting off ${\zeta}(z_0) y_n^{(i)}$ and ${\zeta}(z)y_n^{(ii)}$ this is just the third independent solution $y_n^{(iii)}$ in Lemma 3.1, up to an overall sign.

We now present an observation that dramatically simplifies the solution of the dual Somos-4 equation and allows the order of (1.9) to be reduced from 4 to 3. The point is that, since ${\mathbb{D}}$ is a commutative algebra, the calculation which shows that (2.9) is a conserved quantity for (1.1) carries over directly to the dual numbers, so the same expression with the replacement $x_n\to X_n$ provides a first integral $J=J^{(0)}+J^{(1)}\varepsilon\in{\mathbb{D}}$ for (1.7). This implies that the system of even/odd equations has a pair of rational first integrals, namely $J^{(0)}$ , which is just the first integral for the even part (1.8) given by the original expression (2.9) in terms of $x_n$ alone, but with parameters ${\alpha}\to{\alpha}^{(0)}$ , ${\beta}\to{\beta}^{(0)}$ , and $J^{(1)}$ , which depends on both sets of variables $x_n$ and $y_n$ , and is linear in $y_n$ . The most efficient way to calculate $J^{(1)}$ is to clear the denominator in the formula for $J\in{\mathbb{D}}$ , multiplying both sides of the relation by $X_nX_{n+1}X_{n+2}X_{n+3}$ , to obtain a polynomial relation, and then the leading order (even) part of the resulting expression just returns the usual formula for $J^{(0)}$ , while the $O(\varepsilon)$ (odd) part gives a linear equation in $J^{(1)}$ , but also contains $J^{(0)}$ . Upon eliminating $J^{(0)}$ from the latter expression, the desired formula for $J^{(1)}$ is obtained.

Lemma 3.6. In addition to the conserved quantity $J^{(0)}$ , given by replacing the parameters ${\alpha}\to{\alpha}^{(0)}$ , ${\beta}\to{\beta}^{(0)}$ in the right-hand side of (2.9), the system consisting of (1.8) and (1.9) has a second rational first integral, namely

(3.9) \begin{equation}J^{(1)} = \frac{D_n-\sum_{j=0}^3C^{(j)}_n x_{n+j}^{-1} y_{n+j}}{x_nx_{n+1}x_{n+2}x_{n+3}},\end{equation}

where

\begin{align*} C_n^{(0)} & = {\alpha}^{(0)}x_{n+1}^3x_{n+3}+{\beta}^{(0)}x_{n+1}^2x_{n+2}^2-x_n^2x_{n+3}^2, \nonumber\\[4pt]C_n^{(1)} & = {\alpha}^{(0)}x_nx_{n+2}^3-2{\alpha}^{(0)}x_{n+1}^3x_{n+3}-{\beta}^{(0)}x_{n+1}^2x_{n+2}^2+x_n^2x_{n+3}^2, \nonumber\\[4pt]C_n^{(2)} & = -2{\alpha}^{(0)}x_nx_{n+2}^3+{\alpha}^{(0)}x_{n+1}^3x_{n+3}-{\beta}^{(0)}x_{n+1}^2x_{n+2}^2+x_n^2x_{n+3}^2, \\[4pt]C_n^{(3)} & = {\alpha}^{(0)}x_nx_{n+2}^3+{\beta}^{(0)}x_{n+1}^2x_{n+2}^2-x_n^2x_{n+3}^2, \nonumber\\[4pt]D_n & = {\alpha}^{(1)}x_nx_{n+2}^3+{\alpha}^{(1)}x_{n+1}^3x_{n+3}+{\beta}^{(1)}x_{n+1}^2x_{n+2}^2. \nonumber\end{align*}

For fixed $J^{(1)}$ , the equation (3.9) can be rearranged as a linear inhomogeneous difference equation of third order for $y_n$ , that is,

(3.10) \begin{equation}L_n (y_n) = F_n,\end{equation}

where $L_n$ is the linear difference operator:

(3.11) \begin{equation}L_n = C_n^{(3)} x_{n+3}^{-1} {\mathcal{S}}^3+ C_n^{(2)} x_{n+2}^{-1} {\mathcal{S}}^2+ C_n^{(1)} x_{n+1}^{-1} {\mathcal{S}}+ C_n^{(0)} x_{n}^{-1},\end{equation}

given in terms of the shift operator ${\mathcal{S}}$ that sends $n\to n+1$ , and

(3.12) \begin{equation}F_n = D_n -J^{(1)}x_nx_{n+1}x_{n+2}x_{n+3}.\end{equation}

Thus, we have succeeded in reducing the order of (1.9) by 1, as claimed. The corresponding third-order homogeneous equation which arises when ${\alpha}^{(1)}={\beta}^{(1)}=J^{(1)}=0$ , namely

(3.13) \begin{equation}L_n (y_n) =0,\end{equation}

is nothing other than the linearization of the equation defining $J^{(0)}$ , and three linearly independent solutions are obtained by varying the analytic solution (2.1) with respect to the three parameters that do not depend on this modular quantity: in other words, up to rescaling and taking linear combinations, they are given by the first three shadow solutions $y^{(i)}_n,y^{(ii)}_n,y^{(iii)}_n$ in Lemma 3.1. Then, a fourth independent Somos-4 shadow solution is found by taking ${\alpha}^{(1)}={\beta}^{(1)}=0$ but $J^{(1)}\neq0$ .

The form of the solutions of the homogeneous equation (3.13) is made more transparent by setting $y_n=x_nY_n$ . This gives

\begin{align*}L_n(y_n)=x_nx_{n+1}x_{n+2}x_{n+3}\tilde{L}_n (Y_n)=0,\end{align*}

where

\begin{align*} \tilde{L}_n = \left(\left(\frac{{\beta}^{(0)}}{d_{n+1}d_{n+2}}-d_{n+1}d_{n+2}\right)({\mathcal{S}}+1) +\frac{{\alpha}^{(0)}}{d_{n+2}} \, {\mathcal{S}}+\frac{{\alpha}^{(0)}}{d_{n+1}}\right)({\mathcal{S}} -1)^2.\end{align*}

So clearly, $Y_n=1$ and $Y_n=n$ lie in the kernel of the latter operator, corresponding to $y^{(i)}_n$ and $y^{(ii)}_n$ .

Finally, we can construct the general solution of (1.9) by applying the discrete analogue of the method of variation of parameters to find an arbitrary element in the affine space of solutions of (3.10). In other words, we can write the solution as the sum of a particular integral plus an arbitrary linear combination of solutions of the homogeneous equation (3.13). For variation of parameters, the initial ansatz is to write the solution of (3.10) in the form:

(3.14) \begin{equation}y_n = \sum_j f_n^{(j)} y_n^{(j)},\end{equation}

where the index j ranges over the three lower-case Roman numerals i,ii,iii, and then impose the constraints that

(3.15) \begin{equation} \sum_j \left(f_{n+1}^{(j)}-f_n^{(j)}\right) y_{n+1}^{(j)}=0= \sum_j \left(f_{n+1}^{(j)}-f_n^{(j)}\right) y_{n+2}^{(j)},\end{equation}

which together imply that $y_{n+1}= \sum_j f_n^{(j)} y_{n+1}^{(j)}$ , $y_{n+2}= \sum_j f_n^{(j)} y_{n+2}^{(j)}$ . Putting all this into (3.10) gives

\begin{align*}L_n(y_n) = C_n^{(3)}x_{n+3}^{-1} \sum_j \left(f_{n+1}^{(j)}-f_n^{(j)}\right) y_{n+3}^{(j)}+ \sum_j f_n^{(j)}L_n \left( y_{n}^{(j)}\right)=C_n^{(3)}x_{n+3}^{-1} \sum_j \left(f_{n+1}^{(j)}-f_n^{(j)}\right) y_{n+3}^{(j)} ,\end{align*}

which must equal $F_n$ . Combining the latter equality with the two constraints (3.15) gives a linear system for the three differences $f_{n+1}^{(j)}-f_n^{(j)}$ , which is solved to yield

Theorem 3.7. The general solution of (3.10) can be written in the form (3.14), where

(3.16) \begin{equation}f_n^{(j)}=f_0^{(j)}+\sum_{k=0}^{n-1}v_k^{(j)}, \qquad j=i,ii,iii, \qquad \mathrm{for}\,\,n\geq 0,\end{equation}

with

(3.17) \begin{equation}\left(\begin{array}{c} v_n^{(i)} \\[4pt] v_n^{(ii)} \\[4pt] v_n^{(iii)} \end{array}\right)=\frac{x_{n+3}F_n}{C_n^{(3)}}\left|\begin{array}{c@{\quad}cc@{\quad}c} y_{n+1}^{(i)} & y_{n+1}^{(ii)} & y_{n+1}^{(iii)} \\[4pt] y_{n+2}^{(i)} & y_{n+2}^{(ii)} & y_{n+2}^{(iii)} \\[4pt] y_{n+3}^{(i)} & y_{n+3}^{(ii)} & y_{n+3}^{(iii)}\end{array}\right|^{-1}\left(\begin{array}{c} y_{n+1}^{(ii)} y_{n+2}^{(iii)}-y_{n+1}^{(iii)} y_{n+2}^{(ii)} \\[4pt] y_{n+1}^{(iii)} y_{n+2}^{(i)}-y_{n+1}^{(i)} y_{n+2}^{(iii)}\\[4pt] y_{n+1}^{(i)} y_{n+2}^{(ii)}-y_{n+1}^{(ii)} y_{n+2}^{(i)}\end{array} \right).\end{equation}

The three parameters $f_0^{(j)}$ for $j=i,ii,iii$ in (3.16) are arbitrary, and together with the freedom to choose $J^{(1)}$ arbitrarily, this in turn provides the general solution of (1.9).

4. Hankel determinant formulae

It was conjectured by Barry and proved by various authors that certain Somos-4 sequences could be expressed as Hankel determinants [Reference Chang and Hu1, Reference Chang, Hu and Xin2, Reference Xin29]. In [Reference Hone12], we showed that these results can be unified and further generalized by applying van der Poorten’s work on Jacobi continued fractions (J-fractions) in hyperelliptic function fields [Reference van der Poorten26]. In the genus 1 case, one expands a certain function G on a quartic curve as a J-fraction, that is,

(4.1) \begin{equation}G=\cfrac{s_0}{X+v_1 -\cfrac {d_2}{ X+v_2-\cfrac{d_3}{X+v_3 -\cdots}} },\end{equation}

where $X^{-1}$ is a local parameter around one of the points at infinity, and the recursion relation for the continued fraction leads to the map (3.4) for $(v_{n-1},d_n)$ . The numerators and denominators of the convergents provide associated orthogonal polynomials, and standard results imply that the power series expansion of the generating function near $X=\infty$ , that is $G=\sum_{j\geq 1} s_{j-1}X^{-j}$ , provides a sequence of moments $(s_j)$ such that $d_n$ and $v_n$ can be written in terms of ratios of the corresponding determinants $\Delta_n$ or $\Delta_n^*$ of Hankel/bordered Hankel type, respectively. In particular, the solution of (2.4) is given by $d_n=\Delta_{n-2}\Delta_n/\Delta_{n-1}^2$ , and in [Reference Hone12] we further showed how this result extends to the solution of an integrable symplectic map associated with the J-fraction expansion of a function on a hyperelliptic curve of any genus $g>1$ .

Table 1. Four independent shadows of the original Somos-4 sequence (1.3).

Since ${\mathbb{D}}$ is a commutative algebra, identities for Hankel determinants carry over directly to suitable sequences of moments $s_j\in{\mathbb{D}}$ and allow the solution of (1.7) to be expressed in the same form as for Somos-4 sequences over ${\mathbb{C}}$ . Thus, simply by writing a dual number version of the statement of Theorem 5.2 in [Reference Hone12] (also making use of the formulae in Theorem 4.1 therein for the particular case of genus 1), and taking care that certain parameters should be units, we arrive at the following result.

Theorem 4.1. Given arbitrary $\hat{\alpha},\hat{\beta},\hat{\gamma},s_0\in{\mathbb{D}}^*$ and $s_1\in{\mathbb{D}}$ , define a sequence of dual numbers $(s_j)_{j\geq 0}$ by the recursion:

(4.2) \begin{equation}s_j = \hat{\alpha} s_{j-2} +\hat{\beta}\sum_{i=0}^{j-2}s_is_{j-2-i}+\hat{\gamma} \sum_{i=0}^{j-3}s_is_{j-3-i}, \qquad j\geq 2,\end{equation}

and form the associated sequence of Hankel determinants:

(4.3) \begin{equation}\Delta_n =\left| \begin{array}{c@{\quad}c@{\quad}c@{\quad}c}s_0 & s_1 & \cdots & s_{n-1} \\[3pt]s_1 & & & \vdots \\[4pt]\vdots & \unicode{x22F0} & & \vdots \\[4pt]s_{n-1} & \cdots & \cdots & s_{2n-2}\end{array}\right| =\det (s_{i+j-2})_{i,j=1,\ldots, n}\end{equation}

for $n\geq 1$ , with the usual convention that $\Delta_0=1$ . Then,

\begin{align*}X_n = \Delta_{n-1} \qquad \mathrm{for}\,\, n\geq 1\end{align*}

is a solution of the dual Somos-4 recurrence (1.7) with coefficients given by:

\begin{align*}{\alpha} = U^2, \qquad {\beta} ={\alpha} F + \frac{1}{4}J^2,\end{align*}

where

\begin{align*}U=-s_0\hat{\gamma}-s_1\hat{\beta}, \qquad F=-\hat{\alpha}-2s_0\hat{\beta},\end{align*}

and the value of the first integral is fixed to be

\begin{align*}J=2\left(s_0\hat{\alpha}\hat{\beta} +s_0^2\hat{\beta}^2+s_1 \hat{\gamma}\right).\end{align*}

Example 4.3. To obtain Hankel determinant formulae for the fourth shadow sequence, as in Example 3.8, we take dual numbers that are $O(\varepsilon)$ perturbations of the corresponding quantities in Example 4.2 of [Reference Hone12]. Setting

(4.4) \begin{equation}\hat{\alpha}= 1-4\varepsilon, \quad \hat{\beta}=1, \quad \hat{\gamma}=1-\frac{3}{2}\varepsilon, \quad s_0=1+\varepsilon, \quad s_1=\frac{1}{2} \varepsilon\end{equation}

gives $U=-1$ , $F=-3+2\varepsilon$ , and hence

\begin{align*}{\alpha}=1, \qquad {\beta}=1, \qquad J=4-\varepsilon,\end{align*}

and this corresponds to ${\alpha}^{(0)}={\beta}^{(0)}=1$ , $J^{(0)}=4$ and ${\alpha}^{(1)}={\beta}^{(1)}=0$ , $J^{(1)}=-1$ as required. Then upon iterating the recursion (4.2) with parameters and initial values as in (4.4), the sequence $(s_j)$ is found to be

\begin{align*}1+\varepsilon, \frac{1}{2} \varepsilon, 2-\varepsilon, 1+2\varepsilon, 6-6\varepsilon,7+2\varepsilon,24-28\varepsilon,41-23\varepsilon, 115-154\varepsilon,236-\frac{527}{2}\varepsilon, \ldots ,\end{align*}

and the corresponding sequence of Hankel determinants begins with $\Delta_0=1$ , $\Delta_1=1+\varepsilon$ ,

\begin{align*}\Delta_2=\left|\begin{array}{c@{\quad}c} 1+\varepsilon & \frac{1}{2}\varepsilon \\[4pt] \frac{1}{2}\varepsilon & 2-\varepsilon \end{array}\right|=2+\varepsilon,\quad\Delta_3= \left|\begin{array}{c@{\quad}c@{\quad}c} 1+\varepsilon & \frac{1}{2}\varepsilon & 2-\varepsilon \\[4pt] \frac{1}{2}\varepsilon & 2-\varepsilon & 1+2\varepsilon \\[4pt] 2-\varepsilon & 1+2\varepsilon & 6-6\varepsilon \end{array}\right|=3+3\varepsilon,\end{align*}

\begin{align*} \Delta_4= \left|\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} 1+\varepsilon & \frac{1}{2}\varepsilon & 2-\varepsilon & 1+2\varepsilon \\[4pt] \frac{1}{2}\varepsilon & 2-\varepsilon & 1+2\varepsilon & 6-6\varepsilon \\[4pt] 2-\varepsilon & 1+2\varepsilon & 6-6\varepsilon & 7+2\varepsilon \\[4pt] 1+2\varepsilon& 6-6\varepsilon & 7+2\varepsilon & 24-28\varepsilon \end{array}\right|=7+10\varepsilon, \ldots,\end{align*}

which (after shifting the index) gives the correct values of $X_1,X_2,X_3,X_4,X_5,\ldots$ for the combination of the original Somos-4 sequence (1.3) together with its odd part, namely the shadow sequence $y_n^{(iv)}$ as in the fourth row of Table 1.

Example 4.4. For the recurrence (1.4), the sequence with initial data given by four 1s with zero odd part was considered in [Reference Ovsienko and Tabachnikov23], which corresponds to $X_{-1}=X_0=X_1=X_2=1$ and ${\alpha}^{(0)}={\beta}^{(0)}={\beta}^{(1)}=1$ , ${\alpha}^{(1)}=0$ with our notation and indexing conventions. In this case, we take

\begin{align*}\hat{\alpha}= 1+\varepsilon, \quad \hat{\beta}=1, \quad \hat{\gamma}=1+\frac{1}{2}\varepsilon, \quad s_0=1, \quad s_1=-\frac{1}{2} \varepsilon,\end{align*}

giving $U=-1$ , $F=-3-\varepsilon$ , so that

\begin{align*}{\alpha}=1, \qquad {\beta}=1+\varepsilon, \qquad J=4+\varepsilon.\end{align*}

Then, the recursion (4.2) yields the sequence $(s_j)$ as:

\begin{align*}1, -\frac{1}{2} \varepsilon, 2+\varepsilon, 1-\varepsilon, 6+4\varepsilon,7,24+18\varepsilon,41+18\varepsilon, 115+98\varepsilon,236+\frac{345}{2}\varepsilon, \ldots ,\end{align*}

and the corresponding sequence of Hankel determinants begins with $\Delta_0=1$ , $\Delta_1=1$ ,

\begin{align*}\Delta_2=\left|\begin{array}{c@{\quad}c} 1 & -\frac{1}{2}\varepsilon \\[4pt] -\frac{1}{2}\varepsilon & 2+\varepsilon \end{array}\right|=2+\varepsilon,\quad\Delta_3= \left|\begin{array}{c@{\quad}c@{\quad}c} 1& -\frac{1}{2}\varepsilon & 2+\varepsilon \\[4pt] -\frac{1}{2}\varepsilon & 2+\varepsilon & 1-\varepsilon \\[4pt] 2+\varepsilon & 1-\varepsilon & 6+4\varepsilon \end{array}\right|=3+2\varepsilon,\end{align*}

\begin{align*} \Delta_4= \left|\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} 1 & -\frac{1}{2}\varepsilon & 2+\varepsilon & 1-\varepsilon \\[4pt] -\frac{1}{2}\varepsilon & 2+\varepsilon & 1-\varepsilon & 6+4\varepsilon \\[4pt] 2+\varepsilon & 1-\varepsilon & 6+4\varepsilon & 7 \\[4pt] 1-\varepsilon& 6+4\varepsilon & 7 & 24+18\varepsilon \end{array}\right|=7+10\varepsilon, \ldots,\end{align*}

of which the odd parts for index $n\geq 2$ give the (conjecturally) positive sequence of integers $1,2,10,48,160,1273,7346,51,394,645,078,...$ as found in [Reference Ovsienko and Tabachnikov23].

Example 4.5. The analogous sequence for the recurrence (1.5) was also presented in [Reference Ovsienko and Tabachnikov23], corresponding to $X_{-1}=X_0=X_1=X_2=1$ and ${\alpha}^{(0)}={\beta}^{(0)}={\alpha}^{(1)}=1$ , ${\beta}^{(1)}=0$ . For this example, we take

\begin{align*}\hat{\alpha}= 1+\varepsilon, \quad \hat{\beta}=1, \quad \hat{\gamma}=1+\frac{1}{2}\varepsilon, \quad s_0=1, \quad s_1=0,\end{align*}

giving $U=-1-\frac{1}{2}\varepsilon$ , $F=-3-\varepsilon$ , so that

\begin{align*}{\alpha}=1+\varepsilon, \qquad {\beta}=1, \qquad J=4+2\varepsilon.\end{align*}

Thus, from the recursion (4.2), the sequence of dual number moments $(s_j)$ is found to be

\begin{align*}1, 0, 2+\varepsilon, 1+\frac{1}{2}\varepsilon, 6+5\varepsilon,7+\frac{13}{2}\varepsilon,24+27\varepsilon,41+\frac{105}{2}\varepsilon, 115+164\varepsilon,236+378\varepsilon, \ldots ,\end{align*}

hence, the corresponding sequence of Hankel determinants begins with $\Delta_0=1$ , $\Delta_1=1$ ,

\begin{align*}\Delta_2=\left|\begin{array}{c@{\quad}c} 1 & 0 \\[5pt] 0 & 2+\varepsilon \end{array}\right|=2+\varepsilon,\quad\Delta_3= \left|\begin{array}{c@{\quad}c@{\quad}c} 1& 0 & 2+\varepsilon \\[5pt] 0 & 2+\varepsilon & 1+\dfrac{1}{2}\varepsilon \\[5pt] 2+\varepsilon & 1+\dfrac{1}{2}\varepsilon & 6+5\varepsilon \end{array}\right|=3+3\varepsilon,\end{align*}

\begin{align*} \Delta_4= \left|\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} 1 &0 & 2+\varepsilon & 1+\dfrac{1}{2}\varepsilon \\[5pt] 0 & 2+\varepsilon & 1+\dfrac{1}{2}\varepsilon & 6+5\varepsilon \\[5pt] 2+\varepsilon & 1+\dfrac{1}{2}\varepsilon & 6+5\varepsilon & 7 +\dfrac{13}{2}\varepsilon\\[5pt] 1+\dfrac{1}{2}\varepsilon& 6+5\varepsilon & 7 +\dfrac{13}{2}\varepsilon& 24+27\varepsilon \end{array}\right|=7+10\varepsilon, \ldots,\end{align*}

so the integer sequence beginning with $1,3,10,59,198,1387,9389,57,983,752,301,...$ , as in [Reference Ovsienko and Tabachnikov23], is obtained from the odd parts for index $n\geq 2$ , and this is also conjectured to be positive.

There is one more use for the $s_j$ , obtained from the moment generating function G in (4.1), that is relevant to shadow sequences, but now it is the classical case of ${\mathbb{C}}$ -valued moments that concerns us. In that setting, the quantities $v_n$ appearing in the continued fraction (4.1) can be written in the form:

(4.5) \begin{equation}v_n = \frac{\Delta_{n-1}^*}{ \Delta_{n-1}}- \frac{\Delta_{n}^*}{ \Delta_{n}} \quad \mathrm{for}\,\,n\geq 1,\end{equation}

where $\Delta_0^*=0$ and

(4.6) \begin{equation} \Delta_n^* =\left| \begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c}s_0 & s_1 & \cdots & s_{n-2} & s_{n} \\[3pt]s_1 & & \unicode{x22F0} & \vdots & \vdots \\[3pt]\vdots & \unicode{x22F0} & & \vdots & \vdots \\[3pt]s_{n-2} & \cdots & \cdots & s_{2n-4} & s_{2n-2} \\[3pt]s_{n-1} & \cdots & \cdots & s_{2n-3} & s_{2n-1}\end{array}\right|\end{equation}

is a bordered Hankel determinant for $n\geq 1$ . Thus, as pointed out in [Reference Hone12], the solution of the system (3.4) is expressed in terms of ratios of these Hankel/bordered Hankel determinants, whenever the moment sequence $(s_j)$ is generated by a recursion of the form (4.2) over ${\mathbb{C}}$ . This yields a bordered Hankel formula for the third shadow Somos-4 sequence.

Proof. Substitution of (4.5) into the formula (3.8) gives a telescopic sum, which simplifies to

\begin{align*} y_n^{(iii)}=-x_n\left(v_0 -\frac{\Delta_{n-1}^*}{\Delta_{n-1}}\right),\end{align*}

and then using $x_n=\Delta_{n-1}$ this gives the required result after removing the multiple $-v_0 x_n$ from the front.