Integrable Mappings from a Unified Perspective

TL;DR

This paper explores two integrable discrete dynamical systems that reveal geometric origins of their integrability, providing insights into combinatorial probability problems and connecting diverse mathematical fields.

Contribution

It offers a unified geometric framework explaining the integrability and closed-form solutions of models in combinatorial probability and dynamical systems.

Findings

Identification of geometric sources of integrability

Explicit closed-form expressions for probability distributions

Connections to asymptotic behavior and diverse mathematical fields

Abstract

Two discrete dynamical systems are discussed and analyzed whose trajectories encode significant explicit information about a number of problems in combinatorial probability, including graphical enumeration on Riemann surfaces and random walks in random environments. The two models are integrable and our analysis uncovers the geometric sources of this integrability and uses this to conceptually explain the rigorous existence and structure of elegant closed form expressions for the associated probability distributions. Connections to asymptotic results are also described. The work here brings together ideas from a variety of fields including dynamical systems theory, probability theory, classical analogues of quantum spin systems, addtion laws on elliptic curves, and links between randomness and symmetry.

Equations

$$x_{n+1}$$

$$y_{n+1}$$

$$x_n=1+g{x}_n\left(x_{n+1}+c{x}_n+x_{n-1}\right),$$

$$x^{\ast }=1+(2+c)g{x}^{\ast 2}.$$

$$x_{\pm }^{\ast }=\frac{1\pm \sqrt{1-4(2+c)g}}{2(2+c)g}.$$

$$J(x^{\ast },x^{\ast })=\left(\begin{array}{cc}\frac{1}{g{x}^{\ast 2}}-c& -1\\ 1& 0\end{array}\right)=\left(\begin{array}{cc}\frac{1-2(c+1)g{x}^{\ast }}{g{x}^{\ast }}& -1\\ 1& 0\end{array}\right).$$

$${\lambda }_{\pm }=\frac{\tau \pm \sqrt{{\tau }^2-4\Delta }}{2}.$$

$$I_{c=1}(x,y)=xy\left(x+y-\frac{1}{g}\right)+\frac{1}{g}(x+y)-\frac{1}{g^2}$$

$$I_{c=0}(x,y)=xy(1-gx)(1-gy)+gxy-x-y+\frac{1}{g}.$$

$$x_{n+1}+x_{n-1}$$

$$R(g)=\sum _{k\ge 0}{r}_k{g}^k$$

$$r_k=3\sum _{i=0}^{k-1}{r}_i{r}_{k-1-i}.$$

$$R=1+3g{R}^2$$

$$R(g)=\frac{1-\sqrt{1-12g}}{6g}.$$

$$R_n(g):=\sum _{k=0}^{\infty }{r}_{n,k}{g}^k$$

$$R_n=1+g{R}_n\left(R_{n-1}+R_n+R_{n+1}\right).$$

$$R_{-1}$$

$$\underset{n\to \infty }{\lim }{R}_n$$

$$R_n=R\frac{(1-x^{n+1})(1-x^{n+4})}{(1-x^{n+2})(1-x^{n+3})}.$$

$$x+\frac{1}{x}=\frac{1-4gR}{gR}.$$

$$T(g)=\sum _{k=0}^{\infty }{t}_k{g}^k\text{mbox}and{T}_j(g)=\sum _{k=0}^{\infty }{t}_{j,k}{g}^k$$

$$t_k=2\sum _{j=0}^{k-1}{t}_j{t}_{k-1-j},$$

$$T=1+2g{T}^2.$$

$$T_j=1+g{T}_j\left(T_{j-1}+T_{j+1}\right)$$

$$T_{-1}=0$$

$$\underset{j\to \infty }{\lim }{T}_j=T.$$

$$T_j=T\frac{(1-Z^{j+1})(1-Z^{j+5})}{(1-Z^{j+2})(1-Z^{j+4})},$$

$$x_{\ast }=\frac{(1+Z{)}^2}{1+Z^2},$$

$$f(x,y)=xy(1-gx)(1-gy)+gxy-x-y$$

$$x_{n+1}$$

$$y_{n+1}$$

$$I(x,y)=xy\left(x+y-\frac{1}{g}\right)+\frac{1}{g}(x+y)-\frac{1}{g^2}$$

Full text

Supported by NSF grant DMS-1615921

The authors thank Joceline Lega for orignally encouraging us to study the problems considered here from a dynamical systems perspective, both theoretically and numerically, and for many helpful conversations in this regard along the way.

1. 1 Introduction

2. 1.1 Outline

3. 2 Background

4. 2.1 The QRT Mapping and the Sakai Classification

5. 2.2 Numerical study

6. 2.3 Features of the $c=1$ system, dPI

7. 2.4 Features of the $c=0$ system

8. 3 Combinatorial and probabilistic applications

9. 3.1 Two combinatorial problems

10. 3.1.1 Quadrangulations

11. 3.1.2 Triangulations

12. 3.2 Combinatorics context and applications

13. 4 Elliptic Parametrization, c=1

14. 4.1 The Generic Orbit

15. 4.2 Separatrix Degeneration to the Combinatorial Orbit

16. 5 Elliptic Parametrization, c=0

17. 5.1 The Generic Orbit

18. 5.2 Separatrix Degeneration to the Combinatorial Orbit

19. 6 Concluding Remarks

20. 6.1 Combinatorial Problems in Non-Autonomous Extensions

21. 6.2 Elliptic Combinatorics

1. Introduction

The modern subject of dynamical systems can often be described in terms of a number of fundamental dichotomies such as discrete versus continuous or integrable versus chaotic. The diversity that these dichotomies represent often invite comparisons across the divides that actually provide new insights into the respective extremes that can, in turn, lead to new mathematical developments. For example in the case of the former dichotomy, methods stemming from combinatorial analysis or arithmetic can be brought to bear on discrete systems that have no immediate analogue in the continuous setting but may still suggest mathematical themes that transcend the dichotomy. Such was the case, in the middle of the last century, when dynamic mappings on the interval and associated period doubling cascades led to notions of universality and phase transitions that were reminiscent of themes from statistical mechanics. Similarly, studies of chaotic systems or systems of more general ergodic type were spurred on by studying perturbations of nonlinear integrable dynamical systems, providing many concrete examples of so-called KAM phenomena with suggestions of connections to stochastic analysis.

More recently new bridges across these divides have emerged stimulated by developments in random matrix theory, analytical combinatorics and related areas of probability theory and mathematical physics. In broad terms the realization that underlies these developments is that random settings in the presence of symmetry leads to new universality classes in the probabilistic or statistical mechanical sense; moreover, the symmetries present here frequently lead to the asymptotic distributions of these universality classes being describable in terms of integrable systems theory. A perhaps by now classical example of this is the asymptotic analysis of the longest increasing sequence in a random permutation [1], a long-standing problem in the area of asymptotic combinatorial probability, first posed by Erdös. It was solved in terms of powerful analytical methods of Riemann-Hilbert analysis arising in the theory of integrable PDE and whose universal distribution, due to Tracy and Widom, is expressible in terms of a Painlevé transcendent. This celebrated result has had ramifications in the study of random walks on Lie algebras[6], quantum integrable systems [23] and non-equilibrium statistical mechanics [24, 20].

What we will discuss in this paper is certainly less extensive in scope than the results just mentioned; however, it is illustrative of these trends and it does bring together in a novel and unified way themes from combinatorial probability theory, dynamical systems theory, symmetry groups, algebraic geometry and asymptotic analysis.

1.1. Outline

We will begin in Section 2 with a careful description of the two dynamical systems we want to study, together with an indication of the source of their integrabilty and a numerical study of their trajectories, particularly the ones of interest to us. Then, in Section 3, we will reveal the combinatorial and probabilistic significance of the special orbits we have identified, as well as the hierarchy of combinatorial generating functions whose elegant closed form expressions will be the object of our analysis in the remainder of the paper. That analysis, in terms of elliptic function theory, its solitonic degenerations, and associated algebro-geometric features of the dynamical phase space, will be detailed in Sections 4 and 5. Finally in Section 6 we will amplify upon the combinatorial and probabilistic relevance of our work and indicate some directions for related future research.

2. Background

The explicit focus of this paper is on two two-dimensional autonomous (that is, whose coefficients are $n$-independent) discrete dynamical systems of the form

[TABLE]

equivalently given by the second-order recurrences

[TABLE]

for $c$ equal to either [math] or $1$, and with system parameter $g$. Both of these systems are integrable in a sense that will be described shortly; but despite that and the fact that superficially they look quite similar, the structure of this integrability looks quite different. However, in the remainder of this subsection we will focus on the similarities.

An orbit, or trajectory, in either of these systems refers to the infinite sequence of points $\left\{(x_{n},y_{n})\right\}$ generated from an initial condition $(x_{0},y_{0})$ according to equations (1). Everything is considered as a function of the system parameter $g$.

Both systems have the same fixed point structure: two fixed points, one hyperbolic and the other a center. A fixed point is of the form $(x^{*},x^{*})$ where $x^{*}$ satisfies

[TABLE]

Therefore, the two fixed points $(x^{*}_{+},x^{*}_{+})$ and $(x^{*}_{-},x^{*}_{-})$ are located at the solutions to this quadratic equation:

[TABLE]

Classification of the fixed points is determined from the Jacobian evaluated at these points:

[TABLE]

With $\tau=\mbox{Tr}(J)$ and $\Delta=\mbox{det}(J)$, the eigenvalues of the Jacobian satisfy $|\lambda_{+}||\lambda_{-}|=\Delta=1$ and are given by

[TABLE]

At $(x^{*}_{+},x^{*}_{+})$ the determinant $\tau^{2}-4\Delta<0$ for all appropriate values of $g$ (that is, $0<g<\frac{1}{4(2+c)}$, in order that the fixed points are real), so both eigenvalues are complex. Therefore they are complex conjugates of each other, forcing that $|\lambda_{+}|=|\lambda_{-}|=1$, therefore the fixed point $(x^{*}_{+},x^{*}_{+})$ is a center of the system. At $(x^{*}_{-},x^{*}_{-})$, the determinant is positive for all $g$, so the eigenvalues are real and distinct. Their magnitudes are thus not constrained to be equal but are inversely proportional, and so there are attracting and repelling directions and the fixed point $(x^{*}_{-},x^{*}_{-})$ is hyperbolic, a saddle.

Besides their fixed point structure, both systems also share the key property of being integrable, in the sense that each one has an integral of motion, or invariant, $I(x,y)$ satisfying $I(x_{n},y_{n})=I(x_{n+1},y_{n+1})$ for all $n$. The invariants have different degrees:

[TABLE]

and

[TABLE]

Level sets $I(x,y)=e$ for constants $e$ we call “energies” foliate the plane with curves; an orbit of the system is restricted to lie on one such curve. The level sets are in general the real loci of elliptic curves (degree three and four curves respectively) with no singularities in the affine (finite) plane; there are also special energy values for which the curves become singular in the affine plane. Among these latter we will be primarily interested in those special values for which the singular degeneration of the curve corresponds to an irreducible curve containing a single separatrix of the dynamical system. In Table 1 these correspond to the energy value $e_{1}$.

One should also note that there are various compactifications of the phase plane which yield an extension of the level sets to infinity where the associated curve can become singular. If one completes to the projective plane $\mathbb{P}^{2}$ then, in the $c=1$ case, the curve extends to be smooth at infinity. However, in the $c=0$ case there are two singularities at infinity for all values of the energy. However, if one extends the affine plane to $\mathbb{P}^{1}\times\mathbb{P}^{1}$ then, for $c=0$, there are no singularities on this completion of the level curve to infinity. This will play a key role in our analysis of the $c=0$ case with more details provided in Section 5.

2.1. The QRT Mapping and the Sakai Classification

There is a systematic construction that one may look to in trying to understand the integrable structure which underlies the types of dynamical systems we consider here. This construction grew out of earlier studies of classical analogues of quantum spin systems of Heisenberg type which, in its current form, is usually attributed to Quispel, Roberts and Thompson from which it derives its name, the QRT mapping. The mapping provides a construction of integral invariants of the type (5) for a large class of systems of discrete Painlevé type. Our case of $c=1$ is in fact the discrete Painlevé I equation (dPI). For full details about the QRT mapping we refer the reader to [25], but for our purposes here we just state the general QRT form for systems of dPI type:

[TABLE]

The case of $c=0$ does not fit directly within this scheme, but it does correspond to a degeneration in which the numerator and denominator on the RHS of (7) have a common factor [26]. It is not immediately clear how to find an invariant when the QRT mapping degenerates, but in our case it was possible to do so.

QRT turns out to be the first step to an even deeper insight into the integrability of the discrete Painlevé systems due to Sakai [11] who established a correspondence between the classification of rational and elliptically fibered algebro-geometric surfaces on the one hand and Painlevé equations, both discrete and continuous, on the other. This correspondence encodes relations between the dynamics, Bäcklund transformations, and generalized Lie algebras. The $c=0$ case of our study does not fit immediately into this corresponence either, but just recently a more elaborate treatement of Sakai’s point configuration method was put forward [21] which can be applied to this case. It will be of future interest to relate this extension to the results we present in this paper.

2.2. Numerical study

Numerical investigations of the systems (1) were done in Python, with simulations taking as input a value for the system parameter $g$ lying in the region $0<g<\frac{1}{4(2+c)}$, an initial condition $(x_{0},y_{0})$ lying in the finite plane, and the maximum value of $n$ to which the calculation should be run. A sequence of points $\left\{(x_{n},y_{n})\right\}_{n=0}^{n_{max}}$ is calculated using (1), and the numbers are output as scatter plots of the orbit. Invariant level sets are viewed via contour plots.

We note that orbits near the hyperbolic fixed point experience a strong pull in the unstable directions, rapidly accumulating numerical error and producing artifacts in the image. In order to combat this source of error, a multiprecision library in Python [19] was used so that computations could all be done to any level of precision desired (100 decimal places of accuracy was sufficient to balance survival of an accurate simulation with computation time, in contrast to the standard approximately 13 digits of accuracy).

In all simulation images, the fixed points are plotted as black dots. The initial condition for every orbit is plotted in yellow, and all other points are plotted as blue dots. Only one orbit is shown per image. Level sets of the integral of motion are plotted as red curves.

Due to the integrable nature of the systems, the shape of each orbit is constrained to a curve given by the level set $I(x,y)=I(x_{0},y_{0})$. Initial conditions which are close to one another yield orbit shapes which are close to one another, except where a separatrix lies between them due to a singularity occurring on an intermediate curve.

2.3. Features of the $c=1$ system, dPI

Four generic orbits are displayed in Figure 1. We observe the presence of closed orbits encircling the center, and of three-branched (cubic elliptic curve) orbits in this system.

The three singular level sets of this system are shown in Figure 2, separately and superimposed on one set of axes.

For sufficiently small positive values of the system parameter $g$, the qualitative shape of the entire system is preserved by continuous deformation. As $g$ approaches a critical value $g_{c}=\frac{1}{12}=0.08\bar{3}$, the nearby elliptic and hyperbolic fixed points approach each other. Figure 3 illustrates this movement and the changing shape of orbits via plotting of the $e_{1}$ separatrix.

2.4. Features of the $c=0$ system

In comparing the following figures with those of the previous section, we point out that, in terms of what immediately meets the eye, the respective trajectories show many similarities with the exception of the four-branched nature of this system, in contrast to the three-branched nature of the $c=1$ system.

3. Combinatorial and probabilistic applications

In this section, we introduce a primary motivation for studying the dynamical systems described in the previous sections. The two systems are directly related to two combinatorial problems coming from the world of map enumeration. A map is a graph which is embedded into a surface so that (i) the vertices are distinct points in the surface, (ii) the edges are curves on the surface intersecting only at vertices, and (iii) if the surface is cut along the graph, what remains is a disjoint union of connected components, called faces, each homeomorphic to an open disk. By graph we mean a collection of vertices, edges, and incidence relations, with loops and multi-edges allowed, taken to be connected unless stated otherwise. By surface we mean a compact connected complex-analytic manifold of complex dimension one, up to orientation-preserving homeomorphism.

Map enumeration has a history dating back to the early 1960s, which is briefly summarized in section 3.2. We begin in section 3.1 by introducing two problems and solutions appearing in the literature which motivated the present study. These problems are basic in the sense that they treat triangulations and quadrangulations, that is, maps of the lowest face degrees. They are interesting because both families of maps have a notion of a distance, and the enumeration knows about this distance, opening up opportunities for questions about random metrics, statistical mechanics on random lattices, random walks in random environments and random combinatorial structures in general. These problems are special in that they have very nice closed form combinatorial solutions which exhibit hints of integrability, something that was not expected during their initial study but is turning out to characterize a whole collection of problems related to map enumeration (see section 3.2).

3.1. Two combinatorial problems

3.1.1. Quadrangulations

The material in this section is from the seminal paper [2], where this problem was first introduced and a solution was proposed. It concerns the enumeration of 4-valent 2-legged (that is, there are also two 1-valent vertices) planar maps, sorted according to a notion of geodesic distance that will be defined presently. It is equivalent to counting rooted planar quadrangulations (that is, the face degree, rather than the vertex degree, is 4) with an origin vertex, again counting according to an appropriate distance [3].

Let $R=R(g)$ denote the generating function for the family of 2-legged 4-valent planar maps:

[TABLE]

where $r_{k}$ denotes the number of 2-legged 4-valent planar maps on $k$ vertices. By Schaeffer’s bijection [27], this is also the generating function for 4-valent blossom trees. Blossom trees are rooted trees whose external vertices (called leaves) are either black or white, and in which a certain number of leaves at every internal vertex are black, depending on the valence of the trees (1 in a 4-valent blossom tree). The subtree structure, shown in Figure 7, is that each internal vertex of a 4-valent $k$-vertex blossom tree has exactly one descendent black leaf and two descendent subtrees on $k-1$ vertices, where a single white leaf is a subtree with zero vertices. The number of ways to choose two such subtrees and arrange them in the ordered tree is

[TABLE]

Therefore, the generating function $R(g)$ satisfies

[TABLE]

with analytic solution

[TABLE]

The geodesic distance of a 2-legged map is the minimal number of edges crossed by a path on the surface connecting the two legs. Schaeffer’s bijection was extended in [2] to keep track of the geodesic distance, preserved through the bijection and called root depth on the tree side, of even-valent blossom trees. The subtree structure decomposition in Figure 7 was again applied to recursively build a blossom tree from its subtrees, and the root depth of the pieces and the whole kept track of using contour walks; the resulting generating function is for maps with an upper bound on the geodesic distance. Furthermore, defining

[TABLE]

for r_{n,k}=\#\{\mbox{2-leg 4-valent k$-vertex planar maps with goedesic distance$\leq n}\} the derivation yields a second-order three-term recurrence relation for the generating functions instead of the quadratic functional relation that determined $R(g)$:

[TABLE]

The recurrence relation has two boundary conditions:

[TABLE]

the initial condition stating that there are no maps with negative geodesic distance, and the limit condition stating that the large $n$ limit amounts to removing the upper bound on geodesic distance.

The recurrence relation (10) with boundary conditions (11) was solved in [2]. The closed-form solution to the recursion is given by

[TABLE]

where $x$ is the solution with modulus less than 1 to the characteristic equation

[TABLE]

The original paper did not contain a proof of this result; it can be directly checked that the closed-form solves the recursion, and the authors provided the formal means by which they discovered this amazing formula. The first proof of this closed-form appeared in [5], by continued fraction expansions and exploiting a remarkable connection to the problem of enumerating planar maps with a boundary. This leads to an expression for the generating functions in terms of symplectic Schur functions. Our approach provides an alternative proof that emphasizes the connection to integrable dynamical systems and natural extensions to Painlevé equations and Riemann-Hilbert analysis (see section 3.2); however, it will be very interesting to study relations between these two approaches.

3.1.2. Triangulations

The second problem comes from the paper [7] and is an enumeration problem of labeled plane trees which are bijectively equivalent to certain triangulations of the sphere (see section 6.2). A plane tree is a planar map with only one face, and the enumeration is of the family of labeled plane trees, rooted with the root labeled 0 and with the labels of two adjacent vertices differing by $\pm 1$. The valence of these trees is unrestricted. Figure 8 shows a tree in this family, with its root labeled in boldface. Such trees are also called embedded, because trees labeled in this way can be naturally embedded into a one-dimensional integer lattice, with the labeling of each vertex dictating its location in the lattice. The maximal label, or equivalently the largest point in the lattice having positive mass, is called the label height. The label height is the appropriate notion of distance for this enumeration problem, corresponding to the geodesic distance of the previous problem.

Define the generating functions

[TABLE]

where t_{k}=\#\{\mbox{trees in the family with k edges}\} and t_{j,k}=\#\{\mbox{trees in the family}\linebreak\mbox{having k$edges and label height$\leq j}\}. A rooted tree either has zero edges (and is simply the root), or it has at least one edge; in the latter case, it can be decomposed (see Figure 9) into two subtrees connected by an edge incident to the root, and the sum of the edges in the two subtrees will be one less than the number of edges in the overall tree. The subtree which doesn’t contain the root is still in the family, by shifting all its labels by $+1$ or $-1$.

The recurrence relation for $t_{k}$,

[TABLE]

is summed to see that the generating function for this family of trees is also given by a quadratic equation:

[TABLE]

Keeping track of label height through this enumeration, the recursion relation obeyed by $T_{j}$ is

[TABLE]

with boundary conditions

[TABLE]

Then [7] contains the following closed-form solution (formally derived, and verifiable by insertion into the recursion) for this family of generating functions:

[TABLE]

where $Z\equiv Z(g)$ is implicitly given by

[TABLE]

where $x_{*}=x^{*}_{-}$ as defined by (4) which also specifies its $g$-dependence.

A discrete integral of motion

[TABLE]

satisfying $f(T_{n-1},T_{n})=f(T_{n},T_{n+1})$ with convergence of $T_{n}$ to its limit $T$ assured by $f(T_{n},T_{n+1})=f(T,T)$ was given in [4].

3.2. Combinatorics context and applications

Systematic enumeration of families of maps began with Tutte in [28, 29, 30, 31]. He studied rooted maps (maps having a distinguished directed edge) and, by decomposing maps at their root, derived functional equations for generating functions which are often multi-variate due to the combinatorial complexity of the maps. Tutte’s Quadratic Method along with Lagrange Inversion gives the asymptotic behavior of many families of maps (see [8] for a summary). This foundation has made the study of map enumeration a signature component of the modern theory of random combinatorial structures. Further research has followed two methodological streams, one algebraic and the other analytic. In the later category one has the subject of analytic combinatorics amply illustrated in [17] on the one hand and methods of Riemann-Hilbert Analysis coming from integrable systems theory that solved the Erdös problem on longest increasing sequences in permutations mentioned earlier. However we believe the most promising developments will come from a combination of both methods. The problems discussed in this paper are a perfect example of that. On the algebraic side was the breakthrough result of Schaeffer [27], already mentioned, which built upon pioneering work by Cori and Vacquelin [10] to establish a bijection between the families of planar quadrangulations/triangulations and appropriate families of trees. On the analytic side was the asymptotic analysis of families of orthogonal polynomials which motivated the postulation in [2] of formulas like (12). The rigorous asymptotic analysis of the Riemann-Hilbert problem for orthogonal polynomials with exponential weights that has been carried out in [12, 13, 14, 15] provides a rigorous background for the analysis and extension of (12). Indeed our original motivation was to understand non-autonomous extensions of (1) that arose from the analysis in the above cited papers. More will be said about that in section 6.1. Here we just close by noting that besides providing a natural derivation of the formulas (12) and (17), the results in this paper provide for a rigorous global analysis of these generating functions as well as the determination of their maximal domains of holomorphy.

4. Elliptic Parametrization, c=1

4.1. The Generic Orbit

In this section we show how the combinatorial problem, corresponding to $c=1$, embeds into a more general integrable dynamical system and provide an elliptic parametrization of the general orbits of this system. This is a discrete analogue of action-angle variables of continuous integrable systems. The approach used here is primarily that of classical function theory using canonical coordinates based on the Weierstrass $\wp$ function. In this perspective the separation of variables provided by (21) may be viewed as a discrete analogue of the Hamilton-Jacobi equation. In the next theorem we will make use of such an elliptic parametrization for generic orbits of the dPI system that has already been developed by [9]. Our main results concern the degeneration of this parametrization and are presented in the next subsection.

Theorem 1**.**

The discrete dynamical system

[TABLE]

with invariant

[TABLE]

satisfying $I(x_{n+1},x_{n})=I(x_{n},x_{n-1})$ has an elliptic parametrization through the Weierstrass $\wp$-function as

[TABLE]

where $v$ is such that $\wp(v)=\frac{4g+1}{12g^{2}}$. The discrete dynamics are recovered in this continuous parametrization via addition of $v$.

Remark: $x_{n}=x(\xi+nv)$.

Proof.

On a level set $I(x,y)=e$ of the invariant, we have the factorization

[TABLE]

into which can be introduced a new parameter, $t$, so that

[TABLE]

Then

[TABLE]

so $x$ and $y$ can be recovered as the roots of the quadratic with coefficients the elementary symmetric polynomials $xy$ and $x+y$:

[TABLE]

Thus

[TABLE]

Set

[TABLE]

This curve is put into Weierstrass normal form,

[TABLE]

by the translation $t=z+\alpha$ for $\alpha=\frac{1}{12g^{2}}(4g+1)$, giving the following for Weierstrass invariants:

[TABLE]

The Weierstrass normal form of an elliptic curve is equivalent to the differential equation

[TABLE]

defining the Weierstrass $\wp$-function and its derivative, under the identification $z=-\wp(\xi)$ and $w=\wp^{\prime}(\xi)$. Make the change of variables

[TABLE]

where (23c) is a definition of $v$ and (23d) follows from the calculation

[TABLE]

Under this change of variables, points $(x,y)=\left(x(\xi),y(\xi)\right)$ on the curve $I(x,y)=e$ are given by the parametrization

[TABLE]

equivalent to equation (20). We separate out the proof of discrete dynamics recovery into Lemma 1. ∎

The following corollary is for computational use later on, and is simply a reformulation of the elliptic parametrization in terms of a different function in the Weierstrass function family.

Corollary 1**.**

The function $x(\xi)$ of Theorem 1 can be expressed entirely in terms of the Weierstrass $\sigma$-function:

[TABLE]

Proof.

We make use of the following two identities ([22], equation III.5.52 and exercise III.5.16) which relate the Weierstrass $\sigma$-function to the $\wp$-function:

[TABLE]

if $Z$ is not a period of $\wp(z)$ and

[TABLE]

Then

[TABLE]

the last line since $\sigma$ is an odd function. ∎

Lemma 1**.**

Under the parametrization of Theorem 1, the discrete dynamics of the system are recovered via addition of $v$. That is,

[TABLE]

where $v$ is such that $\wp(v)=\frac{4g+1}{12g^{2}}$.

Proof.

The lemma rests solely on the use of certain properties of the Weierstrass functions, including the Weierstrass $\zeta$-function given by $\zeta^{\prime}(u)=-\wp(u)$: the fact that $\zeta$ is an odd function, the identity ([22] equation (5.59))

[TABLE]

and the addition law for $\wp$ ([22] Theorem 5.15)

[TABLE]

We first calculate that

[TABLE]

starting from (20), expanding, and using (29):

[TABLE]

Demonstration of the result (28b) is now direct via evaluation at $\xi+v$:

[TABLE]

Demonstration of the other part of the result proceeds in a few more steps. First compute the product $x(\xi)y(\xi)$ using the forms in (20), expanding along the way $\wp^{\prime}(v)^{2}-\wp^{\prime}(\xi)^{2}=-4(2\wp(v)+\wp(\xi))(\wp(v)-\wp(\xi))^{2}-\frac{2\wp^{\prime}(v)}{g}(\wp(v)-\wp(\xi))$ using the differential equation; thus

[TABLE]

Starting with the left-hand side of the following equality, expanding $\wp(\xi+v)$ using the addition formula (30), and then making a substitution based on (20a) give that

[TABLE]

which is inserted into

[TABLE]

the result of summing $x$ and $y$ from (20) evaluated at $\xi+v$, along with (28b). Inserting (32) into (33) yields $\wp(v)-\wp(\xi+v)=-x(\xi)\left(x(\xi)+y(\xi)-\frac{1}{g}\right)$, and recall that

[TABLE]

Putting all of these together into (34) the remaining result (28a) is concluded:

[TABLE]

∎

4.2. Separatrix Degeneration to the Combinatorial Orbit

In the previous section, elliptic parametrizations were given for the $c=1$ system, relying on the invariant $I$. Orbits of the system lie on level sets $I(x,y)=e$, and we have so far treated the generic case in which $I(x,y)=e$ is a nonsingular curve. It is in this generic setting that results of the previous section hold. In this section, we study certain cases in which $I(x,y)=e$ is singular, and we extend the parametrizations to certain such curves, in particular, the ones of combinatorial interest. Recall that Table 1 contains all the energy values resulting in singular level sets for the two invariants. We will show at the appropriate time that the energy $e_{1}$ in the $c=1$ case is the one corresponding to geodesic distance for quadrangulations. However, this can be observed in Figure 2 as it is the $e_{1}$ separatrix which contains the hyperbolic fixed point, with coordinates given by $(x^{*}_{-},x^{*}_{-})$ and in the combinatorial problems the limit $x^{*}=\lim_{n\rightarrow\infty}x_{n}$ of the generating functions is the analytic function given by the root $x^{*}_{-}$.

Theorem 2**.**

The discrete integrable recurrence

[TABLE]

with invariant

[TABLE]

and limit $\lim_{n\rightarrow\infty}x_{n}=x_{-}^{*}$ satisfying $x_{-}^{*}=1+3gx_{-}^{*2}$ possesses a solution

[TABLE]

where $\nu$ is such that $\sinh(\nu\sqrt{-3r_{2}})=\sqrt{\frac{-3r_{2}}{\alpha+r_{2}}}$ and $\chi=\exp\{-2\nu\sqrt{-3r_{2}}\}$.

In particular, the combinatorial generating functions $R_{n}(g)=x_{n}$ and $R(g)=x_{-}^{*}$ specialized by the additional boundary condition $R_{-1}=0$ have the closed form

[TABLE]

for $\chi$ satisfying $\chi+\frac{1}{\chi}=\frac{1-4gR}{gR}$, by specializing $\xi=-2\nu$.

An orbit $\{(x_{n},y_{n})\}$ satisfying $\lim_{n\rightarrow\infty}x_{n}=x^{*}_{-}$ limits to the hyperbolic fixed point $(x^{*}_{-},x^{*}_{-})$, and thus the energy of the invariant level set on which this orbit lies can be calculated by $e=I(x^{*}_{-},x^{*}_{-})=I(R,R)$.

[TABLE]

In the following Proposition, the elliptic parametrization is extended to orbits on the level set $I(x,y)=e_{1}$, and the proof of the Theorem is deferred until after this degeneration of the parametrization is established.

Proposition 1**.**

The elliptic parametrization of Theorem 1 and Corollary 1 extends to a parametrization of the singular curve

[TABLE]

as

[TABLE]

where the function $\sigma_{d}$ is a degeneration of the Weierstrass $\sigma$-function, given explicitly by the formula

[TABLE]

where

[TABLE]

The related functions in this degeneration of the family of Weierstrass functions are given explicitly by

[TABLE]

Proof.

With the energy $e$ specialized to $e_{1}(g)$, the two variables $g$ and $e$ in the system are tied together, and we pursue further calculations after expressing everything in sight in terms of a new fundamental variable, $\gamma$, which suppresses the square root appearing in $e_{1}$:

[TABLE]

On the singular curve $I(x,y)=e_{1}$, at least two roots of $-4z^{3}+g_{2}z-g_{3}$ coalesce. We therefore seek a factorization

[TABLE]

where $r_{1}$ and $r_{2}$ are not assumed to be distinct. Equating coefficients in the above equation gives that $r_{1}+2r_{2}=0$, $g_{2}=-4(2r_{1}r_{2}+r_{2}^{2})$, and $g_{3}=-4r_{1}r_{2}^{2}$. Therefore, $r_{1}=-2r_{2}$, $g_{2}=12r_{2}^{2}$, and $g_{3}=8r_{2}^{3}$. Thus,

[TABLE]

The identification $z=-\wp_{d}(\xi)$ and $w=\wp_{d}^{\prime}(\xi)$ is made just as before, to a function $\wp_{d}(\xi)$ and its derivative satisfying the same differential equation as the Weierstrass $\wp$-function,

[TABLE]

The $\wp$-function is given as the inverse of an elliptic integral ([22], Eq. (5.28))

[TABLE]

and the related, degenerate, function $\wp_{d}(\xi)$ is calculated directly by the same integral:

[TABLE]

The result of the integration is

[TABLE]

and inverting the equation gives an explicit formula for this degeneration of the Weierstrass elliptic function:

[TABLE]

The derivative (with respect to $\xi$) is given by

[TABLE]

The Weierstrass $\zeta$-function, an anti-derivative of $\wp$, is defined by ([22] 5.40)

[TABLE]

Integrating $\wp_{d}(\xi)$ and pinning down the integration constant by the limit condition using the first few terms of a Taylor expansion for $\coth(y)$ about $y=0$, we have the degeneration

[TABLE]

The Weierstrass $\sigma$-function is a logarithmic antiderivative of $\zeta$ ([22] 5.44) given by

[TABLE]

which can be calculated by the integral ([22] 5.45)

[TABLE]

Thus by integrating the degenerate $\zeta_{d}$,

[TABLE]

with integration constant $C=\log\left(\frac{1}{\sqrt{-3r_{2}}}\right)$, and the degenerate $\sigma_{d}(\xi)$ is given by

[TABLE]

Theorem 1 now goes through with all of the Weierstrass functions replaced by their appropriate degenerations. The proofs of Corollary 1 and Lemma 1 rely on three identities of the Weierstrass functions: equations (26), (29), (30), and the oddness of $\zeta$, which can all be verified to hold in the degenerate case as well. The proofs are just algebra along with hyperbolic trigonometric identities; one could, for example, use the Pythagorean Theorem to write all hyperbolic trigonometric functions in terms of the hyperbolic cotangent, and then use simply the addition laws for the hyperbolic cotangent. Once these identities are established, all of the general parametrization and dynamics machinery extend directly to the $I(x,y)=e_{1}$ separatrix and Equation (40) is immediate. ∎

Proof of Theorem 2.

The results of Proposition 1 give the degeneration of the elliptic parametrization appropriate for the Theorem, as established by the calculation of the energy $I(x^{*}_{-},x^{*}_{-})=e_{1}$ in equation (39). Inserting the explicit formula for $\sigma_{d}$ given in equation (45) into the parametrization of equation (40), observe that all of the exponentials cancel out and an overall factor of $\sqrt{-3r_{2}}$ rests in the numerator of the large fraction. Then, scaling the variables temporarily to suppress factors of $\sqrt{-3r_{2}}$ appearing in every argument, via

[TABLE]

and simplifying using the double-angle formula for the hyperbolic sine, we have

[TABLE]

By definition of each hyperbolic trigonometric function in terms of exponentials, and defining

[TABLE]

then

[TABLE]

We Taylor expand for large $n$ (note that $\chi^{n}$ is small). First, the various pieces:

[TABLE]

The product of the last two simplifies as

[TABLE]

since for $k\geq 2$ the following sum holds by telescoping the series:

[TABLE]

Now putting all these pieces together, we have

[TABLE]

The constants out front are simplified by first observing that

[TABLE]

by writing the fraction in $\chi$ as a product of the hyperbolic sine and cosine and squaring both sides of the equation, so that everything can be reduced to its expression in terms of $\gamma$. Next, we claim that

[TABLE]

which is easily proved by expressing $\sinh\left(v\sqrt{-3r_{2}}\right)$ in terms of $\chi$ and then reducing everything in sight to its expression as a function of $\gamma$. Using this second identity, it is also easily proved that

[TABLE]

from which the expression for $x(\xi+nv)$ is finally re-summable:

[TABLE]

∎

5. Elliptic Parametrization, c=0

5.1. The Generic Orbit

The discrete system (1) with $c=0$ has the invariant, (6),

[TABLE]

which is quartic, or in fact bi-quadratic. This does not lend itself naturally to a canonical Weierstrass parametrization as in the case of $c=1$ where the invariant was cubic. Therefore we will be taking a more geometric approach to the analysis of the $c=0$ case.

For notational convenience we will sometimes pass to the translated energy variable

[TABLE]

To understand the elliptic parametrization of a generic level set of (6) we pass to a projective completion of the $x-y$ plane as $\mathbb{P}^{1}\times\mathbb{P}^{1}$ in which the dynamical system becomes

[TABLE]

This stems from its more usual form in affine coordinates:

[TABLE]

This mapping is reversible with inverse given by

[TABLE]

The invariant (46) for this system can also be expressed in homogeneous coordinates:

[TABLE]

One may check directly that the level sets of this invariant are generically smooth and, in particular, they are all smooth along the lines at infinity: $x_{0}=0$ and $y_{0}=0$.

Now we may consider the Segre embedding of $\mathbb{P}^{1}\times\mathbb{P}^{1}$ into $\mathbb{P}^{3}$ given by

[TABLE]

It is immediate from this representation that the Segre map embeds $\mathbb{P}^{1}\times\mathbb{P}^{1}$ as a quadratic surface in $\mathbb{P}^{3}$ whose equation is given by

[TABLE]

It is also straightforward to check that, for each value of $E$, the invariant curve (5.1) corresponds to the intersection of the surface $\{F=0\}$ with another quadric surface in $\mathbb{P}^{3}$ explicitly given by

[TABLE]

It is well known that the intersection of two smooth quadrics in $\mathbb{P}^{3}$ is, generically, a smooth elliptic curve: a smooth space curve of degree 4, which we will sometimes refer to as a space quartic.

Returning to energy parameter $e$, the invariant is seen to factor as

[TABLE]

We work with two additional models of the level set associated to this fundamnetal model of the elliptic space curve. These are based on two projections:

[TABLE]

[TABLE]

To understand the geometric meaning of these projections, set $t=\left(gZ_{3}+Z_{0}\right)$, and also note that from (51 - 52) one has

[TABLE]

(53 - 54) show that the image of the elliptic space curve under $\pi_{1}$ is nothing but the invariant level set in the original coordinates, (46), completed in the projective plane $\mathbb{P}^{2}$. This image is, for generic values of $E$, a smooth quartic curve in the affine plane but with two double points on the line at infinity at $[0:1:0]$ and $[0:0:1]$ respectively. Since the space quartic is, generically, a smooth space curve, it follows that the pull-back $\pi^{-1}_{1}$ effectively desingularizes the plane quartic; in fact, it amounts to a blow-up of $\mathbb{P}^{2}$ at $[0:1:0]$ and $[0:0:1]$.

On the other hand, (55) shows that $\pi_{2}$ is projection onto the $t$-line where $t=1+gxy$ which presents the space curve as a double cover of $\mathbb{P}^{1}$. To see this we note from (55) that fixing $t$ determines $\frac{Z_{3}}{Z_{0}}=xy$; then, (56) determines $\frac{Z_{1}+Z_{2}}{Z_{0}}=x+y$. So then $Z_{2}/Z_{0}=x$ and $Z_{1}/Z_{0}=y$ are determined up to two choices by the sign of $\Delta$ in (57). This shows that the elliptic involution of the underlying abstract elliptic curve is realized on the space quartic by the rational involution which interchanges $Z_{1}$ and $Z_{2}$ in $\mathbb{P}^{3}$. (Consistent with this, we note that this involution preserves both $F$ and $G_{E}$ and hence must induce an involution of the space quartic.) We can make use of this involution to explicitly define the double cover of the $t$-line.

Recall that the Weierstrass points of an elliptic curve are the fixed points of its elliptic involution. This is realized in different ways in the various models of the curve that we have been describing. In the Segre (space) model we have just seen that the involution is realized by exchanging $Z_{1}$ and $Z_{2}$. Hence the fixed points in this model are the four points of intersection of the space quartic with the plane $Z_{1}-Z_{2}=0$. Making the consequent substitutions $Z_{2}=Z_{1},Z_{0}=1,Z_{3}=Z_{1}^{2}$ in the factored form $G_{e}$, one reduces to

[TABLE]

Setting the first factor equal to $t$ and the second factor to $e/t$ and solving for $t$ in terms of $e$ and $g$ one derives the equation for the Weierstrass points in terms of the zeroes of

[TABLE]

More commonly these are referred to as the branch points of the projection $\pi_{2}$ from the space quartic onto the $t$-line. We have denoted the function in (58) by $\Delta(t)$ because it is straightforward to check that this quartic is proportional to $\Delta$ defined in (57), consistent with the fact that the elliptic involution is given by interchanging $Z_{1}$ and $Z_{2}$. Indeed, the elliptic irrationality corresponds to $\sqrt{\Delta}$.

The elliptic involution in the plane curve model is realized by exchanging $x$ and $y$: $(x,y)\to(y,x)$ on the curve. The four Weierstrass points of the invariant curve in the plane are the fixed points determined by setting $x=y$ in (46). So these fixed points are of the form $(x,x)$ where $x$ solves

[TABLE]

We can relate these points to the zeroes of $\Delta(t)$ which we denote by $t_{i},i=0,\dots,3$. Using (56) these respective coordinates of the Weierstrass points can be related by setting $x=y$:

[TABLE]

which then must also be the solutions of (59).

We may further make use of these observations about the elliptic involution to directly parametrize $x$ and $y$ in terms of elliptic functions. The function $Z_{i}$ restricted to the quartic curve which is the intersection $\{F=0\}\cap\{G_{E}=0\}$ is necessarily a fourth order theta function on the curve [16]. It has four zeroes on the curve which we can determine by algebra. We illustrate this in the case of $Z_{0}$. If $Z_{0}=0$, it follows from the form of $F$ that one must also have that either $Z_{1}=0$ or $Z_{2}=0$. It then also follows from the form of $G_{e}$ that either $Z_{3}=0$ or $gZ_{3}-Z_{1}-Z_{2}=0$. Hence the points of intersection in $\{Z_{0}=0\}\cap\{F=0\}\cap\{G_{e}=0\}$ are $[0:0:1:0]+[0:1:0:0]+[0:0:g:1]+[0:g:0:1]$ representing them as a formal linear combination of the intersection points. We call this the divisor of $Z_{0}$ and denote it by $(Z_{0})$. One can similarly work out the divisors of the other $Z_{i}$. Here are the results:

[TABLE]

From (53 - 54) we can then also determine the divisors of $x$ and $y$ for which there is significant cancellation:

[TABLE]

Now let us pass to the projection $\pi_{2}$ and consider the uniformizing variable $\xi$ given by the abelian integral of the first kind

[TABLE]

This sets up a mapping

[TABLE]

from a point on the space quartic to the parameter $\xi$ on the universal cover of the elliptic curve. Modulo the periods of $\xi$, this map is a globally analytic homeomorphism.

By standard results on the Weierstrass sigma function, $\sigma$ [22], the uniformization (62) and the divisor information (60 - 61) it follows that the functions $x$ and $y$ can be paramerized as

[TABLE]

where we have the correspondence

[TABLE]

with $\xi$-coordinates on the left and points on the space quartic on the right. The $c_{i}$ are some constants to be determined. We mentioned earlier that the $\pi_{1}^{-1}$ desingularizes the plane quartic by pulling back to the Segre embedding. We can now make this more precise by observing from the above that it is the two points $a_{1}$ and $b_{3}$ that lie over $[0:1:0]$, and $b_{1}$ and $b_{2}$ that lie over $[0:0:1]$. (We will observe in the next paragraph that $b_{1}=-a_{1}$ and $b_{3}=-b_{2}$ so that these pairs are elliptic involutes of one another just as their base points are.)

Using the fact that the elliptic involution on the space quartic is given by the interchange of $Z_{1}$ and $Z_{2}$, it follows that the following pairs are in involution on the $\xi$-plane

[TABLE]

But in $\xi$ the elliptic involution is given by sign change, meaning, that

[TABLE]

so that the sigma function representations reduce to the following fundamental parametrization

[TABLE]

We also record here, for later use, an application of the fundamental divisor relation for elliptic functions [22]:

[TABLE]

where $2\omega_{1},2\omega_{2}$ are fundamental periods for the elliptic space curve.

To help pin down the coefficients $c_{i}$, let us observe from the $\xi-(x,y)$ correspondence above that

[TABLE]

Using the fact that $\sigma$ is an odd function in (63 - 64) one has

[TABLE]

Then setting $\xi=a_{2}$ in this relation and using the previous identifications yields

[TABLE]

Hence, $c_{1}=c_{2}\doteq c$. It follows that $x(\xi)$ and $y(\xi)$ are in fact elliptic involutes of each other. As a corollary we see that elliptic involution of the level set in the $(x,y)$-plane is induced by the coordinate exchange $(x,y)\to(y,x)$ and that the second component (49) of the dynamic map is itself induced by elliptic involution. This calculation shows that

[TABLE]

Let us now relate the above observations to the parametrization of the map (48 - 49). The fact is that this dynamical system expresses an elegant geometric construction.

Consider an initial condition $(x_{0},y_{0})$ that corresponds to a point $\xi_{0}$ on the abstract elliptic curve and its first iterate $(x_{1},y_{1})$ corresponding to another point $\xi_{1}$ on the abstract elliptic curve. In other words

[TABLE]

where $x(\xi)$ is defined by (63).

Starting with the initial point one fixes the vertical line $x=x_{0}$. This line meets the plane quartic (46) in four points, two of which are the (multiplicity 2) double point $[0:1:0]$ at infinity. A third point is, of course, the initial point $(x_{0},y_{0})$, which we know lies on the plane quartic. To find the remaining point is a matter of algebra; its $y$-coordinate is the other root of

[TABLE]

After cancellations we find that this root is simply

[TABLE]

So the residual point of intersection of the line $x=x_{0}$ with the invariant curve is

[TABLE]

Finally, applying the elliptic involution to this point yields our map

[TABLE]

So, in summary, our dynamic map can be described in completely geometric terms: Starting with an initial point $(x_{0},y_{0})$ on the planar invariant curve one finds the residual point of intersection of the line $x=x_{0}$ with the curve and then flips it across the diagonal $x=y$ to produce $(x_{1},y_{1})$. Clearly this process can be arbitrarily iterated. It can also be reversed to yield the inverse map by starting with the vertical line $y=y_{0}$ and proceeding analogously.

Hence our map determines an algebro-geometric addition law on each of the affinely smooth level sets of the invariant. (We shall see in the next section that this extends to singular affine level sets as well.) Indeed, this type of addition law has been studied before [18] and shown to linearize on the abstract elliptic curve. Hence, there is a phase $\nu$ such that $\xi_{1}=\xi_{0}+\nu$ and the higher iterates $x_{n}=x(\xi_{n})$ are determined by $\xi_{n}=\xi_{0}+n\nu$. To pin down this phase we apply Abel’s theorem, as described in [18], which in our setting implies that if $A,B,C,D$ are four points of intersection of a line in the plane with the planar quartic (46) then

[TABLE]

where $K$ is a constant independent of the planar line one considers. So for our case, in choosing the line $x=x_{0}$ one has

[TABLE]

(The factor of 2 here comes from the fact that $[0:1:0]$ is a double point on the curve.) But in fact, as we observed, this linear relation holds for any iterate and so one has in general

[TABLE]

To pin down the constant $K$ one can consider the line $x=0$ in $\mathbb{P}^{2}$. Starting with the standard projectivization in coordinates $[x:y;z]$ of (46),

[TABLE]

Considering this in the vicinity of the point $[0:1:0]$ and setting $y=1$, the invariant has the form

[TABLE]

To leading order near $[0:1:0]$ this has the form $0=gx(gx-z)$ + higher order terms, and so this point is a double point of the curve (as we had already observed) and the two branches of the curve there are $x=0$ and $z=gx$. Hence the line $x=0$ is tangent to the first branch and the line $x=0$ must meet this double point with multiplicity 3. Indeed, this is confirmed by setting $x=0$ in the LHS of (69) which evaluates to $-z^{3}$. There is therefore just one further point of intersection of the line with the affine part of the curve (68). By setting $z=1$ and $x=0$ to get

[TABLE]

one sees that the third point of intersection is $[0:-E:1]$ which corresponds to $a_{2}$ in the Segre model. The points at infinity in the intersection of $x=0$ with the curve correspond to $2a_{1}$, because of the tangency, and $b_{3}$. It then follows from Abel’s theorem that

[TABLE]

where the second line follows from the involution pairings and the third line follows from (65). Therefore the RHS of (67) vanishes modulo periods.

Putting all this together one has

[TABLE]

which is clearly independent of $n$. The evaluations in the second and third equalities are made using the definitions of the $a_{j},b_{j}$ and the relations among them are determined by involution and (65).

Finally, we are in a position to present an elliptic parametrization of the map (48 - 49):

[TABLE]

Making use of (63) in (71) we see that the poles on the LHS of (71) are located at $\xi=-\nu-a_{1}=a_{2}$ and $\xi=-\nu+b_{2}=3a_{1}+2a_{2}$ while those on the RHS must occur among the values $\xi=a_{1},\xi=a_{2}$ and $\xi=-b_{2}$. (One of these latter values cannot occur as a pole in order to maintain the required balance between the number of poles on the two sides of equation (71); i.e. upon passing to a common denominator on the LHS, one of the zeroes of the numerator must cancel one these possible polar values.) So upon comparison we see that $a_{2}$ is a pole. On the other hand $a_{1}$ must be the potential pole on the RHS that gets cancelled, since otherwise we would deduce that $a_{2}=-a_{1}$ which we know from our earlier analysis is not the case. Hence we may deduce that

[TABLE]

5.2. Separatrix Degeneration to the Combinatorial Orbit

We now specialize our considerations to the case of the invariant level curve which passes through the unique hyperbolic fixed point of the map (47) which has energy level $e_{1}$. This energy is given as a function of $g$ in Table 1:

[TABLE]

which we will continue to denote by $e$ in this section. Setting $g=\dfrac{1-\gamma^{2}}{8}$ we may express this and other parameters relevant to the degenerate level set purely in terms of $\gamma$:

[TABLE]

where $x_{*}$ denotes the fixed point, (4), of the separatrix. Given these definitions we also have a concise expression for the eigenvalues of the linearized combinatorial orbit at $x_{*}$

[TABLE]

In addition, the $t$-equation degenerates to

[TABLE]

whose roots are

[TABLE]

By completing the square in the quadratic factor, the equation for $t$ may be rewritten as

[TABLE]

Here $\sigma_{d}$ denotes the generalized sigma function that $\sigma$ degenerates to on the singular curve that (46) limits to. We need to determine the form of $\sigma_{d}$. We observe that the equation (59), determining the Weierstrass points, degenerates in this limit to

[TABLE]

The sigma function may be explicitly expressed in terms of the odd Jacobi theta function $\vartheta_{1}$ [22] as

[TABLE]

where $v=\frac{\xi}{2\omega_{1}}$ and $\vartheta_{1}$ depends on the modular parameter $\tau=\omega_{2}/\omega_{1}$ where $\omega_{1}$ is real. $\eta_{1}$ is an integration constant of the Weierstrass $\zeta$-function. The separatrix limit here in which $t_{2}$ and $t_{3}$ coalesce (or equivalently $x(\omega_{2})$ and $x(\omega_{3})$ coalesce) corresponds to the infinite period limit in which $\tau\to\infty$. In this limit, $\vartheta_{1}$ must converge to a linear combination of exponentials of the form $e^{s\xi}$ and $e^{-s\xi}$ which is odd since $\vartheta_{1}$ remains odd throughout the limit. Hence in fact it must be a multiple of $\sinh(s\xi)$. Then, setting $s=\lim\frac{1}{2\omega_{1}}$ and $\rho=\lim\frac{\eta_{1}}{2\omega_{1}}$, the limiting form is

[TABLE]

One has the implicit integral inversion:

[TABLE]

Setting $z=r_{2}\cosh(x)$ this transforms to

[TABLE]

where $u=e^{x}$. Unravelling the substitutions and setting $z=r_{2}w$ yields

[TABLE]

This last expression presents $\xi$ as a closed form branched 2:1 cover of the $w$ plane or, by scaling and translation, the $t$-plane. We can also pin down the values at branch points:

[TABLE]

We can apply the above inversion analysis to the near-separatix linearization. Setting $z=-r_{1}+q$ in (87) gives

[TABLE]

Comparing to the inversion of (87) detailed above, one may infer that

[TABLE]

and relate this to the eigenvalues at the linearization at the separatrix: the linearization realized in (5.1) presents the dynamics as a phase increment on the covering space of the planar quartic. This fact extends to the degenerate limit of the separatrix. This flow is the tangent space flow near the fixed point and so must coincide with the linearized flow in the directions of the stable or the unstable manifolds. But these flows, we know, are just the power flows of the respective eigenvalues, which are the exponentials of these phase translations. It follows that

[TABLE]

One can now say more about the dynamical system (48 - 49) on the separatrix level set. The parametrization on this level set becomes

[TABLE]

We note that $2a_{1}-b_{2}+a_{2}$ is the limit of divisors of an elliptic function as noted in (65) and therefore must equal the limit of a linear combination of periods. By re-centering the origin of $\xi$ we may assume this sum is zero. Independently, this is also verified from the dynamical systems perspective since, otherwise, $x$ would blow up or vanish as $n$ goes to infinity which our phase plane analysis has shown is not the case. (Indeed, this argument shows that the real part of $2a_{1}-b_{2}+a_{2}$ must already be zero, before any centering.) Hence, the expression for the first component reduces (using $b_{2}=2a_{1}+a_{2}$) to

[TABLE]

In the separatrix limit, as $\xi\to\infty$, we must have

[TABLE]

Hence, one finally has

[TABLE]

To establish (17) it remains to impose the initial condition (16). In our current framework this means that we need to find $\xi_{0}$ such that $x(\xi_{0}-\nu)=0$. This is achieved in terms of (91) by setting

[TABLE]

With this in place, the combinatorial orbit is given by

[TABLE]

where the last line is justified by (73) which persists under degeneration because Abel’s theorem does. In particular one sees that $\frac{a_{2}-a_{1}}{a_{1}+a_{2}}=4$ is equivalent to $5a_{1}+3a_{2}=0$.

Finally we need to relate this to Bousquet-Melou’s $Z$ in (17). Comparing to (89) one has $\chi=\chi^{+}=e^{-2s\nu}$ and so one may rewrite the combinatorial flow as

[TABLE]

It then suffices to show that the eigenvalue $\chi$ satisfies the equation (18) for $Z$. Starting with (79), and using the identity $2gx^{2}_{*}=x_{*}-1$ for the fixed point one has

[TABLE]

which is easily inverted by appropriate squaring to yield

[TABLE]

which finally establishes (17).

We close with a summary of the key results in this section:

Theorem 3**.**

The discrete integrable recurrence

[TABLE]

with invariant

[TABLE]

and limit $\lim_{n\rightarrow\infty}x_{n}=x_{-}^{*}$ satisfying $x_{-}^{*}=1+2gx_{-}^{*2}$ possesses a solution

[TABLE]

where $\nu=2\xi([0:1:0])$ under the composite Abel map determined by (87), $\pi_{1}$ and $\pi_{2}$. Here $\chi=(\frac{1+\sqrt{\gamma}}{1-\sqrt{\gamma}})=e^{-2s\nu}$ where $\gamma=\sqrt{1-8g}$ and $s=\frac{1}{2}\sqrt{r_{1}^{2}-r_{2}^{2}}$. $\xi_{0}$ is determined by the initial condition $x(\xi_{0}-\nu)=0$.

Finally, making the combinatorial identification $T_{n}(g)=x_{n}$ and $T(g)=x^{*}_{-}$ one has

[TABLE]

for $\chi$ satisfying $T=\frac{(1+\chi)^{2}}{(1+\chi^{2})}$.

6. Concluding Remarks

In this last section we mention a few directions that are currently under investigation for building on what was derived in this paper.

6.1. Combinatorial Problems in Non-Autonomous Extensions

The continuous Painlevé equations that give rise to Painlevé’s famous six transcendents are non-autonomous differential equations. In the autonomous limits of these equations the transcendent limits to a classical function such as the elliptic functions discussed in this paper.

The discrete Painlevé equations also have natural non-autonomous extensions. In the case of dPI, we consider such an extension of (1), with $c=1$ having the general form

[TABLE]

$g,N$ and $\zeta$ are parameters. Setting $\zeta=-1=\frac{n}{N}$ one recovers our autonomous dPI equation. One may now ask, are there orbits of this non-autonomous system having combinatorial significance related to what we saw in the autonomous case? The answer is yes and it comes from the analysis of the asymptotics of orthogonal polynomials and their relation to the enumeration of quadrangulations of surfaces of general genus. (The combinatorial problem studied in the autonomous case was for just planar maps; i.e., maps of genus 0.) This connection is mediated by the analysis of the Riemann-Hilbert problem for orthogonal polynomials as it relates to hermitian random matrix models, carried out in [12, 13, 15]. We briefly outline the essentials of this as it relates to (97).

Define the orthogonal polynomials $p_{n}(\lambda)=\gamma_{n}\lambda^{n}+\cdots$ for positive $\gamma_{n}$ with respect to the exponential weight

[TABLE]

for $g>0$, so that

[TABLE]

for $n,m\geq 0.$ These polynomials satisfy the three-term recurrence relation

[TABLE]

for $n\geq 0$ and $p_{-1}=0$. Since the weight is even, the recurrence coefficients $a_{n,N}$ are all zero, and the polynomials are entirely defined by the recurrence coefficients $b_{n,N}$. These coefficients satisfy a nonlinear recurrence of their own:

[TABLE]

In approximation theory (98) is referred to as Freud’s equation. It is straightforward to see that this coincides with (97) if one sets $x_{n}=b^{2}_{n,N}$. The initial conditions that then relate these orthogonal polynomials recurrence coefficients to an orbit within our non-autonomous dPI system, (97), are $x_{0}=0$ and $x_{1}=\frac{c_{2}}{c_{0}}$ where

[TABLE]

is the $j^{th}$ moment of the measure. The latter is derived from the first equation in the recurrence relations,

[TABLE]

by the normalization requiremnt that $p_{1}(\lambda)=\frac{\lambda\gamma_{0}}{b_{1}}$ has norm 1:

[TABLE]

so that

[TABLE]

Similarly, it is required that $p_{0}(\lambda)\gamma_{0}$ have norm 1:

[TABLE]

So finally one has the second initial condition:

[TABLE]

The connection to the combinatorial problem of enumerating quadrangulations now stems from a result in [13] where it is shown that for $N$ large and with $X\doteq\frac{n}{N}\sim 1$, one has a full asymptotic expansion of $x_{n}=b^{2}_{n,N}$:

[TABLE]

where $s=-X\frac{g}{4}$. This expansion is uniformly valid on compact subsets of complex $s$ with $\Re s<0$. The coefficients $z_{g}(s)$ are the generating functions for 4-valent maps (whose dual maps are the quadrangulations in question) of genus $g$.

It will be of interest to determine how the combinatorial features just described are related to the dynamic properties of the orbit, corresponding to the orthogonal polynomial recurrence relations, of non-autonomous dPI. It will also be interesting to see if this can be related to the analysis of the autonomous case discussed in this paper. There is already one indication that such relations do hold: the generating function, $z_{0}(s)$, which enumerates planar 4-valent maps solves the same functional equation, (8), as does the limit $R(g)$ of the distance generating funcitons $R_{n}(g)$. Thus the non-autonomous $x_{n}$ and the autonomous $R_{n}$ agree at leading order in large $n$. These questions are currently under further investigation and will be reported on elsewhere.

As indicated at the start of section 3.1.2 there is an analogous relation between the combinatorics of (Eulerian) triangulations and the $T_{j}$ generating functions [3]. The connection to recurrence relations for (generalized) orthogonal polynomials is more complicated on several levels in this case. For one thing, the coefficients $a_{n}$ will now no longer automatically be zero. However, the analysis in [14] has shown how to handle this additional complication. Based on this we are currently exploring non-autonomous extensions of the $c=0$ system.

6.2. Elliptic Combinatorics

The combinatorial focus of this paper has been on the separatrix orbit of (1). Is there a related combinatorial significance for the other, generic, elliptically parametrized orbits? One affirmative answer to this question has been provided in [3]. There the study of random embedded trees, corresponding to the sytem (1) with $c=0$ but with the added structure of walls is discussed. A wall introduces a conditioning on the random system that strictly bounds the size of the labels/ positions that the random process can attain. A “one-wall” case which is that the half-line bounded below at 0 or -1 reproduces the model we have been studying in this paper (introduced in section 3.1.2) that corresponds to the separatrix orbit. By contrast, a“two-wall” condiitoning would require labels to take their values in a finite subinterval of $\mathbb{Z}$; i.e., it would replace the asymptotic boundary condition (16 b) by another finite boundary condition. From the dynamical point of view this corresponds to a two-point finite boundary value problem. The one-wall case we studied in this paper corresponded to a two-point semi-infinite boundary value problem. One may derive recursive formulas for the generating functions of such trees in the general two point boundary value problem of a two-wall conditioning. The effect of this is to select one of the generic orbits we described in section 5 whose closed form solutions for the generating functions will now be in terms of elliptic functions. These generating functions may be used to study the continuum scaling limit in terms of the periods of the associated elliptic curve, which may be applied to the probabilistic analysis of population spreading. We are studying how the geometric analysis developed in section 5 might be used in this application.