Séminaire Lotharingien de Combinatoire 54 (2006), Article B54g

THE FERMAT CUBIC, ELLIPTIC FUNCTIONS, CONTINUED FRACTIONS, AND A COMBINATORIAL EXCURSION

ERIC VAN FOSSEN CONRAD AND PHILIPPE FLAJOLET

Kindly dedicated to Gérard ⋯Xavier Viennot on the occasion of his sixtieth birthday.

Abstract. Elliptic functions considered by Dixon in the nineteenth century and related to Fermat’s cubic, x^{3} +y^{3} = 1, lead to a new set of continued fraction expansions with sextic numerators and cubic denominators. The functions and the fractions are pregnant with interesting combinatorics, including a special Pólya urn, a continuoustime branching process of the Yule type, as well as permutations satisfying various constraints that involve either parity of levels of elements or a repetitive pattern of order three. The combinatorial models are related to but different from models of elliptic functions earlier introduced by Viennot, Flajolet, Dumont, and Françon.

In 1978, Apéry announced an amazing discovery: “ζ(3) ≡ ∑ 1/n^{3} is irrational”. This represents a great piece of Eulerian mathematics of which van der Poorten has written a particularly vivid account in [59]. At the time of Apéry’s result, nothing was known about the arithmetic nature of the zeta values at odd integers, and not unnaturally his theorem triggered interest in a whole range of problems that are now recognized to relate to much “deep” mathematics [38, 51]. Apéry’s original irrationality proof crucially depends on a continued fraction representation of ζ(3). To wit:

(1)  \begin{matrix}ζ(3) = \frac{6}{ϖ(0) − \frac{1^{6}}{ϖ(1) − \frac{2^{6}}{ϖ(2) − \frac{3^{6}}{...}}}}, \\ where ϖ(n) = (2n + 1)(17n(n + 1) + 5). \end{matrix}

(This is not of a form usually considered by number theorists.) What is of special interest to us is that the nth stage of the fraction involves the sextic numerator n^{6}, while the corresponding numerator is a cubic polynomial in n. Mention must also be made at this stage of a fraction due to Stieltjes (to be later rediscovered and extended

Date: July 9, 2005; accepted March 2, 2006; revised March 25, 2006. Key words and phrases. Combinatorial analysis, elliptic function, continued fraction, permutation,

urn process.

In Séminaire Lotharingien de Combinatoire (Journal), volume 54 (2006) [special issue, Viennot Festschrift] —available electronically. Presented at the 54th Séminaire Lotharingien de Combinatoire (Conference), Lucelle, April 2005. The work of E. Conrad was supported in part by National Security Agency grants MDA 904-93-H-3032, MDA 904-97-1-0019 and MDA 904-99-1-0003, and by National Science Foundation grant DMS-0100288.

by Ramanujan [4, Ch. 12]), namely,

(2)  \begin{matrix}∑_{k=1}^{∞} \frac{1}{(k + z)^{3}} = \frac{1}{σ(0) − \frac{1^{6}}{σ(1) − \frac{2^{6}}{σ(2) − \frac{3^{6}}{...}}}}, \\ where σ(n) = (2n + 1)(n(n + 1) + 2z(z + 1) + 1). \end{matrix}

Unfortunately that one seems to have no useful arithmetic content.

Explicit continued fraction expansions of special functions are really rare [48, 62]. Amongst the very few known, most involve numerators and denominators which, at depth n in the continued fraction, depend rationally on n in a manner that is at most quartic. In this context, the Stieltjes-Ramanujan-Apéry fractions stand out. It then came as a surprise that certain functions related to Dixon’s nineteenth century parametrization of Fermat’s cubic X^{3} + Y^{3} = 1 lead to continued fractions that precisely share with (1) and (2) the cubic–sextic dependency of the their coefficients on the depth n, for instance,

(3) \begin{matrix}∫_{0}^{∞} sm(u)e^{−u/x} du = \frac{x^{2}}{1 + b_{0}x^{3} − \frac{1 · 2^{2} · 3^{2} · 4 x^{6}}{1 + b_{1}x^{3} − \frac{4 · 5^{2} · 6^{2} · 7 x^{6}}{1 + b_{2}x^{3} − \frac{7 · 8^{2} · 9^{2} · 10 x^{6}}{...}}}} , \\ where b_{n} = 2(3n + 1)((3n + 1)^{2} + 1). \end{matrix}

There, the function sm is in essence the inverse of a _{2}F_{1}–hypergeometric of type F [\frac{1}{3}, \frac{2}{3}; \frac{4}{3}]; see below for proper definitions. This discovery, accompanied by several related continued fractions, was first reported in Conrad’s PhD thesis [12] defended in 2002. It startled the second author with a long standing interest in continued fractions [19, 20, 21], when he discovered from reading Conrad’s thesis in early 2005, that certain elliptic functions could precisely lead to a cubic–sextic fraction. This paper describes our ensuing exchanges. We propose to show that there is an interesting orbit of ideas and results surrounding the Fermat cubic, Dixon’s elliptic functions, Conrad’s fractions, as well as certain urn models of probability theory, Flajolet’s combinatorial theory of continued fractions, and finally, the elementary combinatorics of permutations.

Plan of the paper. The Dixonian elliptic functions, “sm” and “cm, are introduced in Section 1. Their basic properties are derived from a fundamental nonlinear differential system that they satisfy. It this way, one can prove simply that they parametrize the Fermat cubic and at the same time admit of representations as inverses of hypergeometric functions. Section 2 presents the complete proof of six continued fractions of the Jacobi type and three fractions of the Stieltjes type that are associated (via a Laplace transform) to Dixonian functions—this is the first appearance in print of results from Conrad’s PhD thesis [12]. Next, in Section 3, we prove that the combinatorics of the nonlinear differential system defining sm, cm is isomorphic to the stochastic evolution of a special process, which is an urn of the Pólya type. As a consequence, this urn, together with its continuous-time analogue (a two-particle version of a classical binary branching process) can be solved analytically in terms of Dixonian functions. Conversely, this isomorphism constitutes a first interpretation of the coefficients of sm, cm phrased in terms of combinatorial objects that are urn histories. Our second combinatorial interpretation is in terms of permutations and it involves peaks, valleys, double rises, and double falls as well as the parity of levels of nodes in an associated tree representation; see Section 4. (Several combinatorial models of elliptic functions due to Dumont, Flajolet, Françon, and Viennot [14, 15, 19, 21, 60] are otherwise known to involve parity-constrained permutations.) Next, the continued fraction expansions relative to Dixonian functions can be read combinatorially through the glasses of a theorem of Flajolet [19] and a bijection due to Françon and Viennot [27], relative to systems of weighted lattice paths and continued fractions. This is done in Section 5 which presents our third combinatorial model of Dixonian functions: the coefficients of sm, cm are shown to enumerate certain types of permutations involving a repetitive pattern of order three. Finally, Section 6 briefly summarizes a few other works where Dixonian functions make an appearance.

Warning. This paper corresponds to an invited lecture at the Viennotfest and, as such, its style is often informal. Given the scarcity of the literature relative to Dixonian functions, we have attempted to provide pointers to all of the relevant works available to us. Thus, our article attempts to kill three birds with one stone, namely be a tribute to Viennot, survey the area, and present original results.

Dedication. This paper is kindly dedicated to Gérard ⋯Xavier Viennot on the occasion of his sixtieth birthday, celebrated at Lucelle in April 2005. His works in lattice path enumeration, bijective combinatorics, and the combinatorics of Jacobian elliptic functions have inspired us throughout the present work.

1. The Fermat curves, the circle, and the cubic

The Fermat curve F_{m} is the complex algebraic curve defined by the equation

X^{m} + Y^{m} = 1.

(Fermat-Wiles asserts that this curve has no nontrivial rational point as soon as m ≥ 3.)

Let’s start with the innocuous looking F_{2}, that is, the circle. Of interest for this discussion is the fact that the circle can be parametrized by trigonometric functions. Consider the two functions from ℂ to ℂ defined by the linear differential system,

(4)  s′ = c, c′ = −s,  with initial conditions s(0) = 0, c(0) = 1.

Then the transcendental functions s, c do parametrize the circle, since

s(z)^{2} + c(z)^{2} = 1,

as is verified immediately from the differential system, which implies (s^{2} + c^{2})′ = 0. At this point one can switch to conventional notations and set

s(z) ≡ sin z, c(z) ≡ cos(z).

Figure 1. The two Fermat curves F_{2} and F_{3}.

It is of interest to note that these functions are also obtained by inversion from an Abelian integral^{1} on the F_{2} curve:

z = ∫_{0}^{sin z} \frac{dt}{\sqrt{1 − t^{2}}} , cos(z) = \sqrt{1 − sin(z)^{2}}

For combinatorialists, it is of special interest to note that coefficients in the series expansions of the related functions

tan z = \frac{sin z}{cos z} , sec z = \frac{1}{cos z} ,

enumerate a special class of permutations, the alternating (also known as “up-and-down” or “zig-zag”) ones. This last fact is a classic result of combinatorial analysis discovered by Désiré André around 1880.

1.1. The Fermat cubic and its Dixonian parametrization. Next to the circle, in order of complexity, comes the Fermat cubic F_{3}. Things should be less elementary since the Fermat curve has (topological) genus 1, but this very fact points to strong connections with elliptic functions.

The starting point is a clever generalization of (4). Consider now the nonlinear differential system

(5)  s′ = c^{2},  c′ = −s^{2},

with initial conditions s(0) = 0, c(0) = 1. These functions are analytic about the origin, a fact resulting from the standard existence theorem for solutions of ordinary differential equations. Then, a one line calculation similar to the trig function case shows that

s(z)^{3} + c(z)^{3} = 1,

since

(s^{3} + c^{3})′ = 3s^{2}c^{2} − 3c^{2}s^{2} = 0.

Consequently, the pair ⟨s(z), c(z)⟩ parametrizes the Fermat curve F_{3}, at least locally near the point (0, 1). The basic properties of these functions have been elicited by

1Given an algebraic curve P (z, y) = 0 (with P a polynomial), an Abelian integral is any integral ∫ R(z, y) dz, where R is a rational function.

Figure 2. Plots of sm(z) [left] and cm(z) [right] for z ∈ ℝ.

Alfred Cardew Dixon (1865–1936) in a long paper [13]. Dixon established that the functions are meromorphic in the whole of the complex plane and doubly periodic (that is, elliptic), hence they provide a global parametrization of the Fermat cubic.

From now on, we shall give the s, c functions the name introduced by Dixon (see Figure 2 for a rendering) and set

sm(z) ≡ s(z), cm(z) ≡ c(z).

Their Taylor expansions at 0 (not currently found in Sloane’s Encyclopedia of Integer Sequences (EIS) [54]) start as follows:

(6)  \{ \begin{matrix}sm(z) & = & z − 4\frac{z^{4}}{4!} + 160\frac{z^{7}}{7!} − 20800\frac{z^{10}}{10!} + 6476800\frac{z^{13}}{13!} − ⋯ \\ cm(z) & = & 1 − 2\frac{z^{3}}{3!} + 40\frac{z^{6}}{6!} − 3680\frac{z^{9}}{9!} + 8880000\frac{z^{12}}{12!} − · · · . \end{matrix}

In summary, we shall operate with the new notations sm, cm and with the defining system:

(7)  sm′(z) = cm^{2}(z), cm′(z) = − sm^{2}(z);  sm(0) = 0, cm(0) = 1.

1.2. A hypergeometric connection. At this point, it is worth noting that one can easily make the s ≡ sm and c ≡ cm functions somehow “explicit.” Start from the defining system (System (7) abbreviated here as (Σ)) and apply differentiation (∂):

(8)  s′ = c^{2} ⟹^{∂} s′′ = 2cc′ ⟹^{Σ} s′′ = −2cs^{2} ⟹^{Σ} s′′ = −2s^{2}\sqrt{s′}.

Then “cleverly” multiply by \sqrt{s′} to get, via integration (∫),

(9)  s′′\sqrt{s′} = −2s^{2}s′ \underrightarrow{∫} \frac{2}{3}(s′)^{3/2} = −\frac{2}{3}s^{3} + K,

for some integration constant K which must be equal to \frac{2}{3}, given the initial conditions. This proves in two lines(!) that sm is the inverse of an integral,

(10)  ∫_{0}^{sm(z)} \frac{dt}{(1 − t^{3})^{2/3}} = z,

this integral being at the same time an incomplete Beta integral and an Abelian integral over the Fermat curve (∫ \frac{dy}{y^{2}}). By the same devices, it is seen that the function cm(z) satisfies

(11)  z = ∫_{cm z}^{1} \frac{dt}{(1 − t^{3})^{2/3}} ,

and is thus also the inverse of an Abelian integral.

Via expansion and term-wise integration, the latest finding (10) can then be rewritten in terms of the classical hypergeometric function,

2F_{1}[α, β, γ; z] := 1 + \frac{α · β}{γ} \frac{z}{1!} + \frac{α(α + 1) · β(β + 1)}{γ(γ + 1)} \frac{z^{2}}{2!} + · · · .

with special rational values of the parameters. Letting Inv(f) denote the inverse of f with respect to composition (i.e., Inv(f) = g if f ∘ g = g ∘ f = Id) we state:

Proposition 1 (Dixon [13]). The function sm is the inverse of an Abelian integral over F_{3} and equivalently the inverse of a _{2}F_{1}:

sm(z) = Inv ∫_{0}^{z} \frac{dt}{(1 − t^{3})^{2/3}} = Inv z ·_{2}F_{1} [\frac{1}{3}, \frac{2}{3}, \frac{4}{3}; z^{3}] .

The function cm is then defined near 0 by cm(z) = \sqrt[3]{1 − sm^{3}(z)}, or alternatively by (11).

The analogy with the sine function is striking. Of course, Proposition 1 is not new and all this is related to extremely classical material. Dixon [13] discusses the implicit integral representations (10), (11) and writes concerning the prehistory of his sm, cm:

The only direct references that I have come across elsewhere are certain passages in the lectures of Professor Cayley and Mr. Forsyth where the integral ∫ \frac{dx}{(1−x^{3})^{2/3}} was used to illustrate Abel’s Theorem, in the treatises of Legendre, and Briot and Bouquet, and again in Professor Cayley’s lectures and elsewhere where it is shewn how to turn the integral into elliptic form, and lastly at the end of Bobek’s Einleitung in die Theorie der elliptischen Functionen where expressions are found for the coordinates of any point on the above cubic.

It is fascinating to be able to develop a fair amount of the theory from the differential equation (7), using only first principles of analysis. (See Dixon’s article as well as Section 6 of the present paper for more.)

Note 1. Lundberg’s hypergoniometric functions. The question of higher degree generalizations of the differential system satisfied by sm, cm is a natural one. Indeed, the system s′ = c^{p−1}, c′ = −s^{p−1} parametrizes locally near (0, 1) the Fermat curve X^{p} + Y^{p} = 1. (This parametrization ceases to be a global one, however, since F_{p} has genus g = (p−1)(p−2)/2, so that g ≥ 3 as soon as p ≥ 4; see for instance [39]. Thus, only the case p = 3 leads directly to elliptic functions.) The corresponding functions are still locally inverses of Abelian integrals over the Fermat curve, which is verified by calculations similar to Equations (8) and (9). In fact, a rather unknown high school teacher from Sweden, Erik Lundberg, developed in 1879 an elementary theory (see [44]) of what he called “hypergoniometric functions”, including a sinualis and a cosinualis that are indexed by rational numbers. Interest in these questions was recently rekindled by an insightful article of Lindqvist and Peetre published in the American Mathematical Monthly [42]. The authors discuss various connections to elliptic functions, in particular, the reduction of Fermat’s cubic to its Weierstraß normal form by elementary manipulations. See also a problem [41] in the same issue of the Monthly posed by these authors together with a remark due to Jon Borwein to the effect that

sm (\frac{π_{3}}{3} − x) = cm(x).

There π_{3} is a fundamental constant of Dixonian functions:

(12) π_{3} = 3 ∫_{0}^{1} \frac{dt}{(1 − t^{3})^{2/3}} = B (\frac{1}{3}, \frac{1}{3}) ≡ \frac{Γ (\frac{1}{3})^{2}}{Γ (\frac{2}{3})} = \frac{\sqrt{3}}{2π}Γ (\frac{1}{3})^{3} ≐ 5.29991 62508.

Many more properties of this type are to be found in Dixon’s paper. In particular, one has,

sm(π_{3} + u) = sm(u), cm(π_{3} + u) = cm(u),

that is, π_{3} is a real period. The complete lattice of periods of sm, cm is

ℤπ_{3} ⊕ ℤπ_{3}ω, ω := e^{2iπ/3},

which is consistent with the rotational invariance: sm(ωu) = ω sm(u) and cm(ωu) = cm(u). It thus corresponds to the hexagonal lattice displayed in Figure 8 and to Case C of the urn models evoked in Subsection 3.5. □

2. Some startling fractions

In 1907, L. J. Rogers [52] devised two methods to obtain continued fraction expansions of Laplace transforms of the Jacobian elliptic functions sn, cn (see for instance [63] for the basic theory of these functions). In his first method, he resorted to integration by parts. His second method involves a general “addition theorem” [48, 62]; it was to some extent a rediscovery of a technique introduced by T. J. Stieltjes in [57] that relies on diagonalization of certain infinite quadratic forms. Stieltjes and Rogers found altogether three families of continued fractions relative to sn, cn. (A fourth family, implicit in the work of Rogers, was discovered almost a century later by Ismail and Masson [33].) Such expansions turn out to be useful: S. Milne in [46] obtained additional expansions implying explicit Hankel determinant evaluations, which enabled him to prove some deep results about sums of squares and sums of triangular numbers. Milne’s results in particular include exact explicit infinite families of identities expressing the number of ways to write an integer as the sum of 4N^{2} or 4N(N + 1) squares of integers, where N is an arbitrary positive integer.

In his PhD thesis defended in 2002, Conrad [12, Ch. 3] applied the integration by parts method to develop completely new continued fractions arising from Dixon’s elliptic functions. These fractions fall into six families which are naturally grouped as two sets of three. The underlying symmetries suggest that this classification into families is fundamental and complete.

In what follows, we make use of the Laplace transform classically defined by

L(f, s) = ∫_{0}^{∞} f(u)e^{−su} du.

For our purposes, it turns out to be convenient to set s = x^{−1}. In that case, one has

x^{−1}L (u^{ν}, x^{−1}) = ν!x^{n},

which means that the Laplace transform formally maps exponential generating functions (EGFs) to ordinary generating functions (OGFs):

(13)  x^{−1}L (∑_{ν≥0} a_{ν} \frac{u^{ν}}{ν!} , x^{−1}) = ∑_{ν≥0} a_{ν}x^{ν}.

There are a number of conditions ensuring the analytic or asymptotic validity of this last equation. In what follows, we shall make use of the integral notation for the Laplace transformation, but only use it as a convenient way to represent the formal transformation from EGFs to OGFs in (13).

The continued fractions that we derive are of two types [48, 62]. The first type, called a J-fraction (J stands for Jacobi), associates to a formal power series f(z) a fraction whose denominators are linear functions of the variable z and whose numerators are quadratic monomials:

f(z) = \frac{1}{1 − c_{0}z − \frac{b_{1}z^{2}}{1 − c_{1}z − \frac{b_{2}z^{2}}{...}}} .

(Such fractions are the ones naturally associated to orthogonal polynomials [11].) The second type^{2}, called in this paper an S-fraction (S stands for Stieltjes) has denominators that are the constant 1 and numerators that are monomials of the first degree:

f(z) = \frac{1}{1 − \frac{d_{1}z}{1 − \frac{d_{2}z}{...}}} .

An S-fraction can always be contracted into a J-fraction, with the corresponding formulæ being explicit, but not conversely. Accordingly, from the point of view of the theory of special functions, an explicit S-fraction expansion should be regarded as a stronger form than its J-fraction counterparts.

2.1. J-fractions for the Dixon functions. We introduce three classes of formal Laplace transforms:

(14)

S_{n} := ∫_{0}^{∞} sm(u)^{n} e^{−u/x}du
C_{n} := ∫_{0}^{∞} sm(u)^{n} cm(u) e^{−u/x}du
D_{n} := ∫_{0}^{∞} sm(u)^{n} cm(u)^{2} e^{−u/x}du

(Since sm^{3} + cm^{3} = 1, we can reduce any polynomial in sm and cm to one which is at most of degree 2 in cm. We let powers of sm grow since sm u vanishes at u = 0.) After integrating by parts and reducing to canonical form using the Fermat cubic parametrization, we obtain the following recurrences:

2An alternative name for S-fraction is “regular C-fraction”. What we called J-fraction is also known as an “associated continued fraction”.

(15)  \begin{matrix}S_{0} = x, & S_{n} = nxD_{n−1} & (n > 0) \\ C_{0} = x − xS_{2}, & C_{n} = nxS_{n−1} − (n + 1)xS_{n+2} & (n > 0) \\ D_{0} = x − 2xC_{2}, & D_{n} = nxC_{n−1} − (n + 2)xC_{n+2} & (n > 0) \end{matrix}

These recurrences will provide six J-fractions and three S-fractions.

For Laplace transforms of powers of sm(u), the recurrences (15) yield the fundamental relation,

(16)  \frac{S_{n}}{S_{n−3}} = \frac{n(n − 1)(n − 2)x^{3}}{1 + 2n(n^{2} + 1)x^{3} − n(n + 1)(n + 2)x^{3} \frac{S_{n+3}}{S_{n}}} ,

subject to the initial conditions

S_{1} = \frac{x^{2}}{1 + 4x^{3} − 6x^{3} \frac{S_{4}}{S_{1}}} , S_{2} = \frac{2x^{3}}{1 + 20x^{3} − 24x^{3} \frac{S_{5}}{S_{2}}} , S_{3} = \frac{6x^{4}}{1 + 60x^{3} − 60x^{3} \frac{S_{6}}{S_{3}}} .

Repeated use of the relation (16) starting from the initial conditions then leads to three continued fraction expansions relative to each of S_{1}, S_{2}, S_{3}.

An entirely similar process applies to C_{n}. We find

\frac{C_{n}}{C_{n−3}} = \frac{n(n − 1)(n − 2)x^{3}}{1 + ((n − 1)n^{2} + (n + 1)^{2}(n + 2)) x^{3} − (n + 1)(n + 2)(n + 3)x^{3} \frac{C_{n+3}}{C_{n}}} ,

subject to

C_{0} = \frac{x}{1 + 2x^{3} − 6x^{3} \frac{C_{3}}{C_{0}}} , C_{1} = \frac{x^{2}}{1 + 12x^{3} − 24x^{3} \frac{C_{4}}{C_{1}}} , C_{2} = \frac{2x^{3}}{1 + 40x^{3} − 60x^{3} \frac{C_{5}}{C_{2}}} .

For D_{n}, we could handle things in a similar manner, but it turns out that there is a simpler procedure available to us. Differentiating powers of sm(u) with respect to u, we have

\frac{d}{du} sm(u)^{n+1} = (n + 1) sm(u)^{n} cm^{2}(u).

Given the effect of Laplace transforms on derivatives, it follows that

L(sm(u)^{n} cm(u)^{2}, x^{−1}) = \frac{1}{(n + 1)x}L(sm^{n+1}(u), x^{−1}).

We thus obtain no substantially new J-fractions from D_{n} by this process.

Theorem 1 (Conrad [12]). The formal Laplace transform of the functions sm^{n} (n = 1, 2, 3) and sm^{n} · cm (n = 0, 1, 2) have explicit J-fraction expansions with cubic denominators and sextic numerators as given in Figure 3.

Consider the list of positive integers in ascending order with each integer listed twice:

1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, …

There are essentially six ways to break this up into ascending 6-tuples, allowing for the possibility of missing leading entries in the first 6-tuple. The a_{n} terms in the Jfraction for the Laplace transform of sm(u), sm^{2}(u), sm^{3}(u) are seen to correspond to

10 ERIC VAN FOSSEN CONRAD AND PHILIPPE FLAJOLET

∫_{0}^{∞} sm u e^{−u/x} du = \frac{x^{2}}{1 + 4x^{3} − \frac{a_{1}x^{6}}{1 + b_{1}x^{3} − \frac{a_{2}x^{6}}{1 + b_{2}x^{3} − ...}}}

where a_{n} = (3n − 2)(3n − 1)^{2}(3n)^{2}(3n + 1), b_{n} = 2(3n + 1) ((3n + 1)^{2} + 1).

∫_{0}^{∞} sm^{2}(u) e^{−u/x} du = \frac{2x^{3}}{1 + 20x^{3} − \frac{a_{1}x^{6}}{1 + b_{1}x^{3} − \frac{a_{2}x^{6}}{1 + b_{2}x^{3} − ...}}}

where a_{n} = (3n − 1)(3n)^{2}(3n + 1)^{2}(3n + 2), b_{n} = 2(3n + 2) ((3n + 2)^{2} + 1).

∫_{0}^{∞} sm^{3}(u) e^{−u/x} du = \frac{6x^{4}}{1 + 60x^{3} + \frac{a_{1}x^{6}}{1 + b_{1}x^{3} − \frac{a_{2}x^{6}}{1 + b_{2}x^{3} − ...}}}

where a_{n} = (3n)(3n + 1)^{2}(3n + 2)^{2}(3n + 3), b_{n} = 2(3n + 3) ((3n + 3)^{2} + 1).

∫_{0}^{∞} cm(u) e^{−u/x} du = \frac{x}{1 + 2x^{3} + \frac{a_{1}x^{6}}{1 + b_{1}x^{3} − \frac{a_{2}x^{6}}{1 + b_{2}x^{3} − ...}}}

where a_{n} = (3n − 2)^{2}(3n − 1)^{2}(3n)^{2}, b_{n} = (3n − 1)(3n)^{2} + (3n + 1)^{2}(3n + 2).

∫_{0}^{∞} sm(u) cm(u) e^{−u/x} du = \frac{x^{2}}{1 + 12x^{3} + \frac{a_{1}x^{6}}{1 + b_{1}x^{3} − \frac{a_{2}x^{6}}{1 + b_{2}x^{3} − ...}}}

where a_{n} = (3n − 1)^{2}(3n)^{2}(3n + 1)^{2}, b_{n} = (3n)(3n + 1)^{2} + (3n + 2)^{2}(3n + 3).

∫_{0}^{∞} sm^{2}(u) cm(u) e^{−u/x} du = \frac{2x^{3}}{1 + 40x^{3} + \frac{a_{1}x^{6}}{1 + b_{1}x^{3} − \frac{a_{2}x^{6}}{1 + b_{2}x^{3} − ...}}}

where a_{n} = (3n)^{2}(3n + 1)^{2}(3n + 2)^{2}, b_{n} = (3n + 1)(3n + 2)^{2} + (3n + 3)^{2}(3n + 4).

Figure 3. The six basic J-fractions relative to Dixonian functions. the partitions:

\begin{matrix}sm(u) : & (1, 2, 2, 3, 3, 4), (4, 5, 5, 6, 6, 7), (7, 8, 8, 9, 9, 10), …, \\ sm^{2}(u) : & (2, 3, 3, 4, 4, 5), (5, 6, 6, 7, 7, 8), (8, 9, 9, 10, 10, 11), …, \\ sm^{3}(u) : & (3, 4, 4, 5, 5, 6), (6, 7, 7, 8, 8, 9), (9, 10, 10, 11, 11, 12), …, \end{matrix}

These are the three possible “odd” partitions, odd in the sense that the first integer in each 6-tuple appears exactly once. The three remaining continued fractions, those of cm, cm · sm, cm · sm^{2}, are associated in the same way with the three possible even partitions.

2.2. S-fractions for Dixon functions. Starting with the original recurrences of the previous section(Equation (15), we can use the relation between S_{n} and D_{n−1} to eliminate either the letter S or the letter D. These recurrences reduce to S-fraction recurrences, which we tabulate here:

\begin{aligned}\frac{S_{n}}{C_{n−2}} & = & \frac{n(n − 1)x^{2}}{1 + n(n + 1)x^{2} \frac{C_{n+1}}{S_{n}}} , & \frac{C_{n}}{S_{n−1}} & = & \frac{nx}{1 + (n + 1)x\frac{S_{n+2}}{C_{n}}} \\ \frac{C_{n}}{D_{n−2}} & = & \frac{n(n − 1)x^{2}}{1 + (n + 1)(n + 2)x^{2} \frac{D_{n+1}}{C_{n}}} , & \frac{D_{n}}{C_{n−1}} & = & \frac{nx}{1 + (n + 2)x\frac{C_{n+2}}{D_{n}}} . \end{aligned}

We need initial conditions to generate S-fractions, and the six candidate starting points give just three S-fraction initial conditions:

\begin{aligned}S_{1} & = & x^{2} − 2x^{2}C_{2} & = & \frac{x^{2}}{1 + 2x\frac{C_{2}}{S_{1}}} \\ C_{0} & = & x − xS_{2} & = & \frac{x}{1 + x\frac{S_{2}}{C_{0}}} \\ C_{1} & = & x^{2} − 2xS_{3} & = & \frac{x}{1 + 2x\frac{S_{3}}{C_{1}}} . \end{aligned}

The three that fail to give good initial conditions for an S-fraction are as follows:

\begin{aligned}S_{2} & = 2x^{3} − 2x^{3}S_{2} − 6x^{2}C_{3} & = \frac{2x^{3}}{1 + 2x^{3} + 6x^{2} \frac{C_{3}}{S_{2}}} \\ S_{3} & = 6x^{4} − 12x^{3}S_{3} − 12x^{2}C_{4} & = \frac{6x^{4}}{1 + 12x^{3} + 12x^{2} \frac{C_{4}}{S_{3}}} \\ C_{2} & = 2x^{3} − 4x^{3}C_{2} − 3x^{2}S_{4} & = \frac{2x^{3}}{1 + 4x^{3} + 3x^{2} \frac{S_{4}}{C_{2}}} . \end{aligned}

On iterating the first three relations, we obtain three S-fraction expansions:

Theorem 2 (Conrad [12]). The formal Laplace transform of the functions sm, cm and sm · cm have explicit S-fraction expansions with cubic numerators,

∫_{0}^{∞} sm u e^{−u/x} du = \frac{x^{2}}{1 + \frac{a_{1}x^{3}}{1 + \frac{a_{2}x^{3}}{1 + ...}}}

where for r ≥ 1: a_{2r−1} = (3r − 2)(3r − 1)^{2}, a_{2r} = (3r)^{2}(3r + 1);

∫_{0}^{∞} cm u e^{−u/x} du = \frac{x}{1 + \frac{a_{1}x^{3}}{1 + \frac{a_{2}x^{3}}{1 + ...}}}

where for r ≥ 1: a_{2r−1} = (3r − 2)^{2}(3r − 1), a_{2r} = (3r − 1)(3r)^{2};

∫_{0}^{∞} sm(u) cm(u) e^{−u/x} du = \frac{x^{2}}{1 + \frac{a_{1}x^{3}}{1 + \frac{a_{2}x^{3}}{1 + ...}}}

where for r ≥ 1: a_{2r−1} = (3r − 1)^{2}(3r), a_{2r} = (3r)(3r + 1)^{2}.

.

Consider now the set of positive integers taken in increasing order with each integer written twice:

1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, …

There are three ways to divide this list into ascending triples if we permit ourselves to discard leading ones, corresponding to the numerators coefficients a_{n} of the S-fractions associated to cm, sm, and sm · cm:

cm(u) :  (1, 1, 2), (2, 3, 3), (4, 4, 5), (5, 6, 6), …,
sm(u) :  (1, 2, 2), (3, 3, 4), (4, 5, 5), (6, 6, 7), …,
sm(u) · cm(u) : (2, 2, 3), (3, 4, 4), (5, 5, 6), (6, 7, 7), ….

This correspondence suggests that these three continued fractions form a complete set.

Note 2. The cubic X^{3} + Y^{3} − 3αXY = 1. Conrad in his dissertation [12], follows Dixon and examines the larger class of functions corresponding to the cubic

X^{3} + Y^{3} − 3αXY = 1,

for arbitrary α. These do give rise to continued fractions, but ones that are non-standard: they are not of the S or J type as they involve denominators that are linear (in x) and numerators that are cubic. □

An alternative derivation of the J-fraction expansions, based on a direct use of the differential system, will be given in Subsection 5.1, when we discuss a method of André (Note 8).

3. First combinatorial model: Balls games

Dixonian elliptic functions, as we shall soon see, serve to describe the evolution of a simple urn model with balls of two colours. In this perspective, the Taylor coefficients of sm, cm count certain combinatorial objects that are urn “histories”, or equivalently, weighted knight’s walks in the discrete plane—this provides our first combinatorial interpretation. Going from the discrete to the continuous, we furthermore find that Dixonian functions quantify the extreme behaviour of a classic continuous-time branching process, the Yule process. Finally, we show that the composition of the system at any instant, whether in the discrete or the continuous, case, can be fully worked out and is once more expressible in terms of Dixonian functions.

3.1. Urn models. Balls games have been of interest to probabilists since the dawn of time. For instance in his Théorie analytique des probabilités (first published in 1812), Laplace writes “Une urne A renfermant un très grand nombre n de boules blanches et noires; à chaque tirage, on en extrait une que l’on remplace par une boule noire; on demande la probabilité qu’après r tirages, le nombre des boules blanches sera x.” Such games were later systematically studied by Pólya. What is directly relevant to us here is the following version known in the standard probability literature as the Pólya urn model (also Pólya-Eggenberger):

Pólya urn model. An urn is given that contains black and white balls. At each epoch (tick of the clock), a ball in the urn is chosen at random (but not removed). If it is black, then α black and β white balls are placed into the urn; if it is white, then γ black and δ white balls are placed into the urn.

The model is fully described by the “placement matrix”,

M = ( \begin{matrix}α & β \\ γ & δ \end{matrix} ) .

The most frequently encountered models are balanced, meaning that α + β = γ + δ, and negative entries in a matrix M are interpreted as subtraction (rather than addition) of balls. For instance, Laplace’s original problem corresponds to ( \begin{matrix}0 & 0 \\ 1 & −1 \end{matrix} ), which is none other than the coupon collector’s problem in modern terminology. What is sought in various areas of science is some characterization, exact or asymptotic, of the composition of the urn at epoch n, given fixed initial conditions. The elementary introduction by Johnson and Kotz [36] mentions applications to sampling statistics, learning processes, decision theory, and genetics. Recently, Pólya urn models have been found to be of interest in the analysis of several algorithms and data structures of computer science; see especially Mahmoud’s survey [45].

The main character of this section is the special urn defined by the matrix

M_{12} = ( \begin{matrix}−1 & 2 \\ 2 & −1 \end{matrix} ) .

This can be visualized as a game with balls of either black (‘x’) or white (‘y’) colour. If a ball is chosen, it is removed [the −1 entry in the matrix] from the urn and replaced by two balls of the opposite colour according to the rules

(17)  x ⟶ yy,  y ⟶ xx.

A history^{3} of length n is, loosely speaking, any complete description of a legal sequence of n moves of the Pólya urn. Let conventionally the urn be initialized at time 0 with one black ball (x). A history (of length n) is obtained by starting with the oneletter word x at time 0 and successively applying (n times) the rules of (17). For instance,

\underline{x} ⟶ y\underline{y} ⟶ y\underline{x}x ⟶ \underline{y}yyx ⟶ x\underline{x}yyx ⟶ xyyyyx,

is the complete description of a history of length 5. (The replaced letters have been underlined for readability.) We let H_{n,k} be the number of histories that start with an x and, after n actions, result in a word having k occurrences of y (hence n + 1 − k occurrences of x). Clearly the total number of histories of length n satisfies

H_{n} := ∑_{k} H_{n,k} = n!,

since the number of choices is 1, 2, 3, … at times 1, 2, 3 . . .. Here is a small table of all histories of length ≤ 3:

n = 0 : x
n = 1 : x ⟶ yy
n = 2 : x ⟶ yy ⟶ xxy

x ⟶ yy ⟶ yxx

n = 3 :  x ⟶ yy ⟶ xxy ⟶ yyxy

x ⟶ yy ⟶ xxy ⟶ xyyy x ⟶ yy ⟶ xxy ⟶ xxxx x ⟶ yy ⟶ yxx ⟶ xxxx x ⟶ yy ⟶ yxx ⟶ yyyx x ⟶ yy ⟶ yxx ⟶ yxyy

The sequence (H_{n,0}) starts as 1, 0, 0, 2 for n = 0, 1, 2, 3 and it is of interest to characterize these combinatorial numbers.

3.2. Urns and Dixonian functions. Let us come back to Dixonian functions. Consider for notational convenience the (autonomous, nonlinear) ordinary differential system

(18)  Σ :  \frac{dx}{dt} = y^{2}, \frac{dy}{dt} = x^{2}, with x(0) = x_{0}, y(0) = y_{0},

which is the signless version of (5). In this subsection, we only need the specialization x(0) = 0, y(0) = 1, but we will make use of the general case (18) in Subsection 3.4 below.

The pair ⟨x(t), y(t)⟩ with initial conditions x(0) = 0, y(0) = 1 parametrizes the “Fermat hyperbola”,

y^{3} − x^{3} = 1,

3A history is the combinatorial analogue of a trajectory or sample path for combinatorial processes.

See Françon’s work [26] on “histoires” for this terminology and similar ideas.

which is plainly obtained from Fermat’s “circle” F_{3} by a vertical reflection. It is natural to denote this pair of functions by ⟨smh(t), cmh(t)⟩, these functions smh, cmh being trivial variants of sm, cm, where the alternation of signs has been suppressed:

(19)  smh(z) = − sm(−z),  cmh(z) = cm(−z).

A first combinatorial approach to Dixonian functions can be developed in a straightforward manner by simply looking at the basic algebraic relations induced by the system Σ of (18). For this purpose, we define a linear transformation δ acting on the vector space ℂ[x, y] of polynomials in two formal variables x, y that is specified by the rules,

(20)  δ[x] = y^{2},  δ[y] = x^{2},  δ[u · v] = δ[u] · v + u · δ[v],

u, v being arbitrary elements of ℂ[x, y]. A purely mechanical way to visualize the operation of δ comes from regarding δ as a rewriting system

(21)  x ⟶^{δ} yy,  y ⟶^{δ} xx,

according to the first two rules of (20). The third rule means that δ is a derivation, and it can be read algorithmically as follows: when δ is to be applied to a monomial w of total degree d in x, y, first arrange d copies of w, where in each copy one instance of a variable is marked (underlined), then apply the rewrite rule (21) once in each case to the marked variable, and finally collect the results. For instance,

xyy ↦ \underline{x}yy, x\underline{y}y, xy\underline{y} ⟶^{δ} (yy)yy, x(xx)y, xy(xx), so that δ[xy^{2}] = y^{4} + 2x^{3}y.

Two facts should now be clear from the description of δ and the definition of histories.

(i) Combinatorially, the nth iterate δ^{n}[x^{a}y^{b}] describes the collection of all the possible histories at time n of the Pólya urn with matrix M_{12}, when the initial configuration of the urn has a balls of the first type (x) and b balls of the second type (y). This is exactly the meaning of the replacement rule (21). In particular, the coefficient^{4},

H_{n,k}^{(a,b)} = [x^{k}y^{ℓ}] (δ^{n}[x^{a}y^{b}]), k + ℓ = n + a + b

is the number of histories of a Pólya urn that lead from the initial state x^{a}y^{b} to the final state x^{k}y^{ℓ}.

(ii) Algebraically, the operator δ does nothing but describe the “logical consequences” of the differential system Σ. In effect, the first two rules of (21) mimic the effect of a derivative applied to terms containing x = x(t) and y = y(t), “knowing” that x′ = y^{2}, y′ = x^{2}. Accordingly, the quantity δ^{n}[x^{a}y^{b}] represents an nth derivative,

δ^{n}[x^{a}y^{b}] = \frac{d^{n}}{dt^{n}} x(t)^{a}y(t)^{b} expressed in x(t), y(t),

where x(t), y(t) solve the differential system x′ = y^{2}, y′ = x^{2}.

4If f = ∑_{m,n} f_{m,n}x^{m}y^{n}, then the notation [x^{m}y^{n}]f is used to represent coefficient extraction:

[x^{m}y^{n}]f ≡ f_{m,n}.

In summary, we have a principle^{5} whose informal statement is as follows.

Equivalence principle. The algebra of a nonlinear autonomous system that is monomial and homogeneous, like Σ in (18), is isomorphic to the combinatorics of an associated Pólya urn.

This principle can be used in either direction. For us it makes it possible to analyse the Pólya urn M_{12} in terms of the functions sm, cm that have already been made explicit in Proposition 1. In so doing, we rederive and expand upon an analysis given in a recent paper on analytic and probabilistic aspects of urn models [22]. Here is one of the many consequences of this equivalence.

Theorem 3 (First combinatorial interpretation). The exponential generating function of histories of the urn with matrix M_{12} that start with one ball and terminate with balls that are all of the other colour is

∑_{n≥0} H_{n,0} \frac{z^{n}}{n!} = smh(z) = \frac{sm(z)}{cm(z)} = − sm(−z).

The exponential generating function of histories that start with one ball and terminate with balls that are all of the original colour is

∑_{n≥0} H_{n,n+1} \frac{z^{n}}{n!} = cmh(z) = \frac{1}{cm(z)} = cm(−z).

Proof. This is nothing but Taylor’s formula. Indeed, for the first equation, we have by the combinatorial interpretation of the differential system:

H_{n,0} = δ^{n}[x]|_{x=0,y=1} .

But from the algebraic interpretation,

δ^{n}[x]|_{x=0,y=1} = \frac{d^{n}}{dt^{n}} smh(t)|_{t=0} .

The result then follows. □

Theorem 3 thus provides a first combinatorial model of Dixonian functions in terms of urn histories.

Note 3. On Dumont. Though ideas come from different sources, there is a striking parallel between what we have just presented and some of Dumont’s researches in the 1980’s and 1990’s. Dumont gives an elegant presentation of Chen grammars in [18]. There, he considers chains of general substitution rules on words: such chains are partial differential operators in disguise. For instance, our δ operator is nothing but

δ = x^{2} \frac{∂}{∂y} + y^{2} \frac{∂}{∂x}.

Dumont has shown in [18] that operators of this type can be used to approach a variety of questions like rises in permutations, Bell polynomials, increasing trees, parking functions, and

5In this paper, we have chosen to develop a calculus geared for Dixonian functions. Our approach is in fact more general and it can be applied to any balanced urn model, as considered in [22] corresponding to a monomial system, of the form {x′ = x^{α}y^{β}, y′ = x^{γ}y^{δ}}, which is homogeneous in the sense that α + β = γ + δ.

[figure: P = (p, q); P = (p, q); multiplicity p; multiplicity q]

Figure 4. The Bousquet-Mélou-Petkovšek moves (left) and their version with multiplicities (right) for the M_{12} urn model.

Pólya grammars. The interest of such investigations is also reinforced by consideration of the trivariate operator

yz \frac{∂}{∂x} + zx\frac{∂}{∂y} + xy \frac{∂}{∂z} ,

itself related to the differential system

x′ = yz, y′ = zx, z′ = xy,

and to elliptic functions of the Jacobian type, which had been formerly researched by Dumont and Schett; see [14, 15, 53] and Subsection 3.5 below for further comments on the operator point of view. □

Note 4. Knight’s walks of Bousquet-Mélou and Petkovšek. Theorem 3 is somehow related to other interesting combinatorial objects (Figure 4). Say we only look at the possible evolutions of the urn, disregarding which particular ball is chosen. An urn with p black balls and q white balls may be represented in the Cartesian plane by the point P with coordinates (p, q). A sequence of moves then defines a polygonal line, P_{0}P_{1} ⋯ P_{n} with P_{0} = (1, 1), each move \overrightarrow{P_{j}P}_{j+1} being of type either β = (−1, +2) [i.e., one step West, two steps North] when a black ball is picked up or ω = (+2, −1) [i.e., two steps East, one step South] otherwise. This defines a random walk in the first quadrant that makes use of two types of knight moves on the chessboard. The enumeration of these walks is a nontrivial combinatorial problem that has been solved recently by Bousquet-Mélou and Petkovšek in [7]. They show for instance that the ordinary generating function of walks that start at (1, 1) and end on the horizontal axis is

G(x) = ∑_{i≥0}(−1)^{i} (ξ^{⟨i⟩}(x)ξ^{⟨i+1⟩}(x))^{2} ,

where ξ is a branch of the (genus 0) cubic xξ − x^{3} − ξ^{3} = 0:

ξ(x) = x^{2} ∑_{m≥0} \binom{3m}{m}\frac{x^{3m}}{2m + 1}.

There, ξ^{⟨i⟩} = ξ ∘ ⋯ ∘ ξ is the ith iterate of ξ. These walks have a merit [7]: they provide an extremely simple example of a (multivariate!) linear recurrence with constant coefficients whose generating function is highly transcendental (in fact, not even holonomic).

In order to obtain a complete history from a knight’s walk, one has to add some supplementary information, namely, which ball is chosen at each stage. This corresponds to introducing multiplicative weights. The rule is then as follows. For a walk with a β-step that starts at point (p, q) (a black ball is chosen) the weight is p; dually, for an ω-step the weight is q. □

In summary, at this stage, we have available three variants of the interpretation of the Taylor coefficients of the Dixonian functions sm, cm provided by Theorem 3: (i) the enumeration of urn histories relative to the M_{12} urn; (ii) the iterates of the special operator δ = x^{2}∂_{y} + y^{2}∂_{x} in the style of Dumont; (iii) multiplicatively weighted knight walks of the type introduced by Bousquet-Mélou and Petkovšek.

Note 5. A probabilistic consequence: extreme large deviations. From a probabilistic standpoint, the Pólya urn model is a discrete time Markov chain with a countable set of states embedded in ℤ × ℤ. The number of black balls at time n then becomes a random variable, X_{n}. Theorem 3 quantifies the probability (ℙ) of extreme large deviations of X_{n} as

ℙ(X_{n} = 0) = \frac{H_{n,0}}{n!} = [z^{n}] smh(z) = [z^{n}] − sm(−z).

Then by an easy analysis of singularities, one finds that this quantity decreases exponentially fast,

ℙ(X_{3ν+1} = 0) ∼ c (\frac{π_{3}}{3} )^{−3ν−1} , where \frac{π_{3}}{3} = \frac{\sqrt{3}}{6π}Γ (\frac{1}{3})^{3} ,

a quantity already encountered in Equation (12). For instance, we have to 10 decimal places (10D):

|\frac{[z^{28}] sm(z)}{[z^{31}] sm(z)}|^{1/3} ≐ 1.76663 87502 ⋯; \frac{π_{3}}{3} ≐ 1.76663 87490 · · · .

One can refer to the calculation of [22], but is is just as easy to note that Proposition 1 provides directly the dominant positive singularity ρ of smh(z) as a special value of the fundamental Abelian integral. □

3.3. A continuous-time branching process. The Dixonian functions also make it possible to answer questions concerning chain reactions in a certain form of particle physics. You have two types of particles, say, foatons and viennons. Any particle lives an amount of time T that is an exponentially distributed random variable (i.e., ℙ(T ≥ t) = e^{−t}), this independently of the other particles; then it disintegrates into two particles of the other type. Thus a foaton gives rise to two viennons and a viennon gives rise to two foatons (see Figure 5). What is the composition of the system at some time t ≥ 0, assuming one starts with one foaton?

Figure 5. A truncated view of a tree of particles provided by the Yule process (the ordinate of each particle represents the time at which it splits into two particles of the other type).

The resemblance with the Pólya urn should be clear (with a foaton being a black x and a viennon a white y). First, considering only the total population of particles, one gets what is perhaps the simplest continuous-time branching process, classically known as the Yule process. (For useful information regarding this process, see for instance the recent articles by Chauvin, Rouault, and collaborators [9, 10].) Let S_{k}(t) be the probability that the total population at time t is of size k. Introduce the bivariate generating function,

Ξ(t; w) = ∑_{k=1}^{∞} S_{k}(t)w^{k}.

Examination of what happens between times 0 and dt leads to the recurrence

(22)  S_{k}(t + dt) = (1 − dt)S_{k}(t) + dt ∑_{i+j=k} S_{i}(t)S_{j}(t),

which is none other than the usual “backwards equation” of Markov processes. Then, one has

S′_{k}(t) + S_{k}(t) = ∑_{i+j=k} S_{i}(t)S_{j}(t).

This induces a nonlinear equation satisfied by Ξ, namely,

Ξ′(t; w) + Ξ(t; w) = Ξ(t; w)^{2}, Ξ(0, w) = w,

where derivatives are implicitly taken with respect to the time parameter t. The solution of this ordinary differential equation is easily found by separation of variables,

Ξ(t; w) = \frac{we^{−t}}{1 − w(1 − e^{−t})},

which yields

(23)  S_{k}(t) = e^{−t} (1 − e^{−t})^{k−1} ,  k ≥ 1.

In summary, the size of the population at time t obeys a geometric law of parameter (1 − e^{−t}), with expectation e^{t}. (The previous calculations are of course extremely classical: they are disposed of in just six lines in Athreya and Ney’s treatise [3, p. 109].)

The result of (23) shows that any calculation under the discrete urn model can be automatically transferred to the continuous branching process. Precisely, let H be the set of all histories of the Pólya urn with matrix M_{12}. Let K ⊆ H be a subset of H and let K(z) be the exponential generating function of K. Then, by virtue of (23), the probability that, at time t, the Yule process has evolved according to a history that lies in K is given by

(24)  e^{−t}K (1 − e^{−t}) .

This has an immediate consequence.

Proposition 2. Consider the Yule process with two types of particles. The probabilities that particles are all of the second type at time t are

X(t) = e^{−t} smh(1 − e^{−t}), Y (t) = e^{−t} cmh(1 − e^{−t}),

depending on whether the system at time 0 is initialized with one particle of the first type (X) or of the second type (Y).

Like its discrete-time counterpart, this proposition can be used to quantify extreme deviations: the probability that all particles be of the second type at time t is asymptotic to

(25)  e^{−t} smh(1),  smh(1) ≐ 1.20541 51514.

At time t, the Yule system has expected size e^{t}. On the other hand, for the discrete process, the probability at a large discrete time that balls or particles be all of one colour is exponentially small. These two facts might lead to expect that, in continuous time, the probability of all particles to be of the same colour is doubly (and not singly, as in Equation (25)) exponentially small. The apparent paradox is resolved by observing that the extreme large deviation regime is driven by the very few cases where the system consists of O(1) particles only, an event whose probability is exponentially small (but not a double exponential).

Note 6. An alternative direct derivation. The memoryless nature of the process implies, like in (22), the differential system:

X(t) + X′(t) = Y (t)^{2}, Y (t) + Y′(t) = X(t)^{2},

with initial conditions X(0) = 0, Y (0) = 1. This system closely resembles some of our earlier equations. By simple algebra [first multiply by e^{t}, then apply the change of variables t ↦ e^{−t}], the solution in terms of Dixonian functions results. □

3.4. Dynamics of the Pólya and Yule processes. The previous sections have quantified the extreme behaviour of the processes—what is the probability that balls be all of one colour? In fact, a slight modification of previous arguments gives complete access to the composition of the system at any instant.

Consider once more the differential system

(26)  \frac{dx}{dt} = y^{2}, \frac{dy}{dt} = x^{2}, with x(0) = x_{0}, y(0) = y_{0}

(this repeats Equation (18)). The initial conditions are now treated as free parameters or formal variables. This system can be solved exactly by means of previously exposed techniques, and, by virtue of the equivalence principle, its solution describes combinatorially the composition of the urn, not just the extremal configurations (that correspond to x_{0} = 0, y_{0} = 1).

First, the solution to the system. Following the chain of (8) and (9), one finds

x′(t)^{3/2} = x(t)^{3} + Δ^{3}, Δ^{3} := y_{0}^{3} − x_{0}^{3},

where the initial conditions have been taken into account. This last equation can be put under the form

\frac{x′}{(x^{3} + Δ^{3})^{2/3}} = 1,

which, after integration, gives

(27)  Δt = ∫_{x_{0}/Δ}^{x/Δ} \frac{dw}{(1 + w^{3})^{2/3}} .

Figure 6. A plot of the coefficients of \frac{1}{n!}δ^{n}[x] for n = 50 and n = 150 illustrating the asymptotically gaussian character of coefficients.

The function smh(t) is also defined by an inversion of an integral,

t = I(smh(t)), where I(y) := ∫_{0}^{y} \frac{dw}{(1 + w^{3})^{2/3}} ,

as is easily verified by the technique of Section 1.2. Then Equation (27) provides

(28)  x(t) = Δ smh (Δt + I (\frac{x_{0}}{Δ} )) ,

which constitutes our main equation.

Second, the interpretation of δ in Section 3.2 implies that δ^{n}[x] expresses the composition of the urn at the nth stage of its operation when it has been initialized with one ball of the first type (x). Thus the quantity

(29)  F (t, x_{0}) = ∑_{n≥0} \frac{t^{n}}{n!} (δ^{n}[x])_{x↦x_{0}, y↦1}

is none other than the x(t) solution to the system (26) initialized with x(0) = x_{0} and y(0) = 1. On the other hand, δ^{n}[x] represents the nth derivative of x ≡ x(t) expressed as a function of x and y.

There results from Equations (28), (29) and the accompanying remarks an expression for the bivariate generating function of histories with t marking length and x_{0} marking the number of balls of the first kind. Switching to more orthodox notations (x_{0} ↦ x), we state:

Proposition 3. The composition of the M_{12} urn at all discrete instants is described by

\begin{aligned}∑_{n,k≥0} H_{n,k}x^{k} \frac{z^{n}}{n!} & = & (1 − x^{3})^{1/3} smh ((1 − x^{3})^{1/3}z + ∫_{0}^{x(1−x^{3})^{−1/3}} \frac{ds}{(1 + s^{3})^{2/3}} ) , \\ & = & (x + y^{2} \frac{z}{1!} + 2x^{2}y \frac{z^{2}}{2!} + (4xy^{3} + 2x^{4})\frac{z^{3}}{3!} + ⋯)_{y=1} , \end{aligned}

where H_{n,k} is the number of histories of length n of the urn initialized with x and terminating with x^{k}y^{n+1−k}.

(M)  (Σ)  (δ)  Type
A \binom{−2 3}{4 −3} x′ = x^{−1}y^{3}, y′ = x^{4}y^{−2} x^{−1}y^{3} \frac{∂}{∂x} + x^{4}y^{−2} \frac{∂}{∂y} ℘ [22, 47];

also Dixonian (§6.1)

B \binom{−1 2}{3 −2} x′ = y^{2}, y′ = x^{3}y^{−1}  y^{2} \frac{∂}{∂x} + x^{3}y^{−1} \frac{∂}{∂y}
C \binom{−1 2}{2 −1} x′ = y^{2}, y′ = x^{2}  y^{2} \frac{∂}{∂x} + x^{2} \frac{∂}{∂y}  Dixonian
D \binom{−1 3}{3 −1} x′ = y^{3}, y′ = x^{3}  y^{3} \frac{∂}{∂x} + x^{3} \frac{∂}{∂y}  lemniscatic [22]
E \binom{−1 3}{5 −3} x′ = y^{3}, y′ = x^{5}y^{−2}  y^{3} \frac{∂}{∂x} + x^{5}y^{−2} \frac{∂}{∂y}
F \binom{−1 4}{5 −2} x′ = y^{4}, y′ = x^{5}y^{−1}  y^{4} \frac{∂}{∂x} + x^{5}y^{−1} \frac{∂}{∂y}

Figure 7. The six elliptic cases of a 2 × 2 Pólya urn: matrix (M), system (Σ), operator (δ), and type, following [22].

This statement is the incarnation of Theorem 1 of [22] in the case of the M_{12} urn model. From it, one can for instance deduce by an application of the Quasi-powers Theorem of analytic combinatorics [22, 24, 31]: The distribution of the number of balls of the first kind, equivalently the sequence of coefficients of the homogeneous polynomial δ^{n}[x], is asymptotically normal. See Figure 6.

The transfer from discrete to continuous time afforded by (24) now permits us to deduce the composition of the Yule process from the Pólya urn via the transformation K(z) ↦ e^{−t}K(1 − e^{−t}).

Proposition 4. In the Yule process with two types of particles, the probability generating function of the number of particles of the first type (x) at time t is

e^{−t}(1 − x^{3})^{1/3} smh ((1 − x^{3})^{1/3}(1 − e^{−t}) + ∫_{0}^{x(1−x^{3}))^{−1/3}} \frac{ds}{(1 + s^{3})^{2/3}} ) .

By means of the continuity theorem for characteristic functions, it can be verified that the distribution of the number of particles of the first type at large times is asymptotically exponential with mean ∼ e^{t}/2. (In this case, the distribution is essentially driven by the size of the system.)

3.5. Complements regarding elliptic urn models. As already noted, the operator δ of the M_{12} urn is a partial differential operator,

δ[f] = y^{2} \frac{∂}{∂x}f + x^{2} \frac{∂}{∂y} f,

which is linear and first order. In a recent study, Flajolet, Gabarró and Pekari [22], investigate the general class of 2 × 2 balanced urn schemes corresponding to matrices of the form

(30)  M = ( \begin{matrix}−a & s + a \\ s + b & −b \end{matrix} ) .

Any such urn is modelled by a particular partial differential operator,

(31)  δ = x^{1−a}y^{s+a} \frac{∂}{∂x} + x^{s+b}y^{1−b} \frac{∂}{∂y} ,  (a, b, s > 0).

Figure 8. The six lattices \binom{A B C}{D E F} corresponding to urn models that are exactly solvable in terms of elliptic functions.

The bivariate generating function of urn histories (a prototype is provided by Proposition 3 under the M_{12} scenario) is in each case expressed in terms of a function ψ that is determined as the inverse of an Abelian integral over the Fermat curve F_{p} of degree p = s + a + b (so that the genus is (p − 1)(p − 2)/2). Reference [22] starts by building the bivariate generating function of histories as H := e^{zδ}[x]; then the method of characteristics is applied in order to solve a partial differential equation satisfied by H. Elementary conformal mapping arguments eventually lead to a complete characterization of all urn models and operators (31) such that H is expressible in terms of elliptic functions. (Caveat: In [22], a model is said to be solvable by elliptic functions if its “base” function ψ is a power, possibly fractional, of an elliptic function. Under this terminology, some of the models eventually turn out to be elliptic, though they are a priori associated to Fermat curves of genus > 1.). It is found in [22] that there are altogether only six possibilities, which correspond to regular tessellations of the Euclidean plane; see Figures 7 and 8.

This classification nicely complements some of Dumont’s researches in the 1980s; see [16]. In particular Dumont developed a wealth of combinatorial connections between the Jacobian elliptic function, Schett’s operator (which involves three variables rather than two), and permutations. The Dumont–Schett result should in particular provide a complete analytic model for the urn with balls of three colours,

M = ( \begin{matrix}−1 & 1 & 1 \\ 1 & −1 & 1 \\ 1 & 1 & −1 \end{matrix} ) ,

corresponding to the trivariate operator introduced by Schett [53]

(32)  δ = yz \frac{∂}{∂x} + zx\frac{∂}{∂y} + xy \frac{∂}{∂z} .

Observe that very few models are known to be explicitly solvable in the case of urns with three types of balls or more^{6}.

The proofs obtained in the present paper rely on the basic combinatorial properties of an associated nonlinear differential system. They are purely “conceptual”, thereby bypassing several computational steps of [22], in particular the method of characteristics. The process based on ordinary differential systems that has been developed here for the M_{12} urn is in fact applicable to all urns of type (30) that are balanced (α + β = γ + δ). This observation suggests, more generally, that interesting combinatorics is likely to be found amongst several nonlinear autonomous systems. Perhaps something along the lines of Leroux and Viennot’s combinatorial-differential calculus [40] is doable here.

4. Second model: permutations and parity of levels

Our second combinatorial model of Dixonian functions is in terms of permutations. It needs the notion of level of an element (a value) in a permutation, itself related to a basic tree representation, as well as a basic classification of elements into four local order types (peaks, valleys, double rises, and double falls).

A permutation can always be represented as a tree, which is binary, rooted, and increasing (see Stanley’s book [55, p. 23]). Precisely, let w = w_{1}w_{2} ⋯ w_{n} be a word on the alphabet ℤ_{>0} without repeated letters. A tree denoted by Tree(w) is associated to a word w by the following rules.

— If w is the empty word, then Tree(w) = ε is the empty tree.

— Else, let ξ = min(w) be the least element of w. Factor w as w′ξw′′. The tree Tree(w) is then inductively defined by

Tree(w) = ⟨Tree(w′), ξ, Tree(w′′)⟩,

that is, the root is ξ, and Tree(w′), Tree(w′′) are respectively the left and right subtrees of the root.

The tree so obtained^{7} is such that the smallest letter of the word appears at the root, and the labels of any branch stemming from the root go in increasing order. This

construction applies in particular to any permutation, σ = ( \begin{matrix}1 & 2 & ⋯ & n \\ σ_{1} & σ_{2} & ⋯ & σ_{n} \end{matrix} ), once

it is written as the equivalent word σ_{1}σ_{2} ⋯ σ_{n}. An infix order traversal of the tree, a projection really, gives back the permutation from the increasing tree, so that the correspondence is bijective; see Figure 9. Observe also that we are dealing with binary trees in the usual sense of computer science [37, §2.3], a nonempty node being of one of four types: binary, nullary (leaf), left-branching, and right-branching.

6Puyhaubert’s thesis [50] provides a discussion of 3 × 3 triangular cases. On the other hand, probabilistic techniques yielding asymptotic information are available for urns with an arbitrary number of colours: see for instance recent studies by Janson and Pouyanne [35, 49]. In particular, Pouyanne (following a suggestion of Janson) obtains results relative to the Yule process with k types of balls.

The rules are of the form x_{1} ⟶ x_{2}x_{2}, x_{2} ⟶ x_{3}x_{3}, …, x_{k−1} ⟶ x_{k}x_{k}, x_{k} ⟶ x_{1}x_{1}, and a curious “phase transition” is found to occur at k = 9.

7In [55, p. 41], Stanley writes that such models of permutations have been “extensively developed primarily by the French” and refers to Foata-Schützenberger [25].

Figure 9. A permutation of length 10 and its corresponding increasing binary tree.

Definition 1. Given a rooted tree, the level of a node ν is defined as the distance, measured in the number of edges, between ν and the root. The level of an element (letter) in a word comprised of distinct letters is the level of the corresponding node in the associated tree.

The statement below also makes use of a fundamental classification of local order types in permutations; see in particular Françon and Viennot’s paper [27].

Definition 2. Let a permutation be written as a word, σ = σ_{1} ⋯ σ_{n}. Any value σ_{j} can be classified into four local order types called peaks, valleys, double rises, and double falls, depending on its order relative to its neighbours:

Peaks  Valleys  Double rises  Double falls
σ_{j−1} < σ_{j} > σ_{j+1} σ_{j−1} > σ_{j} < σ_{j+1} σ_{j−1} < σ_{j} < σ_{j+1} σ_{j−1} > σ_{j} > σ_{j+1}

For definiteness, border conditions must be adopted, and we shall normally opt for one of the two choices,

σ_{0} = −∞, σ_{n+1} ∈ {−∞, +∞}.

(For instance, an alternating permutation is characterized by the fact that it has only peaks and valleys.) The corresponding quadrivariate exponential generating function of all permutations was first determined by Carlitz.

Theorem 4 (Second combinatorial model). Consider the class X of permutations bordered by (−∞, −∞) such that elements at any odd level are valleys only. Then the exponential generating function is

X(z) = smh(z) = − sm(−z).

For the class Y of permutations also bordered by (−∞, −∞) such that elements at any even level are valleys only, the exponential generating function is

Y (z) = cmh(z) = − cm(−z).

Proof. This statement^{8} easily results from first principles of combinatorial analysis. The interpretation is obtained by examining the differential system

X′(z) = Y (z)^{2}, Y′(z) = X(z)^{2}, X(0) = 0, Y (0) = 1,

or, under an equivalent integral form,

X(z) = ∫_{0}^{z} Y (w)^{2} dw, Y (z) = 1 + ∫_{0}^{z} X(w)^{2} dw.

It suffices to note that, if A and B are two combinatorial classes with exponential generating functions A and B, then the product A · B enumerates the labelled product A ⋆ B. Also, the integral ∫ A · B enumerates all well-labelled triples ⟨α, ξ, β⟩, where ξ is the smallest of all labels and α, β are of respective types A, B.

By the preceding remarks, we see that X, Y are the generating function of increasing trees satisfying the relations

X = ⟨Y, min, Y⟩, Y = ε + ⟨X , min, X ⟩.

It is then apparent that trees of type X have only double nodes at odd levels, while those of type Y have only double nodes at even levels. The statement of the proposition finally results from the correspondence between tree node and permutation value types, to wit,

double node ↔ valley , left-branching node ↔ double fall,
leaf ↔ peak, right-branching node ↔ double rise,

which is classical (and obvious via projection). □

Consequently, for n = 3ν, the number of permutations of type Y is Y_{n} = n![z^{n}] cmh(z). We have

Y_{0} = 1 Y_{2} = {ϵ} (the empty permutation)
Y_{3} = 2 Y_{2} = {213, 312},

which agree with (6). In order to form the shapes of trees of type Y, one can use the grammar represented graphically as (□ represents the empty tree)

Y = □ + [figure: Y; Y; Y; Y]

This automatically generates trees whose sizes are multiples of 3. (The number of such tree shapes of size n = 3ν is the same as the number of quaternary trees of size ν, namely, \frac{1}{3ν+1} \binom{4ν}{ν}). The increasing labellings of these trees in all possible ways followed by projections provide a way of listing all permutations of type Y. (See Figure 10 for an illustration of this construction.) In this way, it is found that there are 4 legal tree shapes of type Y for n = 3ν = 6, each of which admitting 10 increasing labellings, which globally corresponds to Y_{6} = 40, in agreement with (6).

8As we learnt in May 2005, this result was also obtained independently by Dumont at Ouagadougou in 1988; see his unpublished note [17].

12 \underline{6} 10 3 \underline{1} 2 5 11 \underline{8} 9 \underline{4} 7

Figure 10. The construction of a Y-permutation enumerated by cm(z). The elements at even level (underlined) are valleys only.

Similarly, trees corresponding to X are determined by the grammar

X = [figure: Y; Y]

which yields, again in agreement with (6):

\begin{matrix}X_{1} = 1 : & X_{1} = \{1\}, \\ X_{4} = 4 : & X_{4} = \{3241, 4231, 1324, 1423\}. \end{matrix}

Note 7. Parity-based permutation models for elliptic functions. There already exist several combinatorial models of the Jacobian elliptic functions sn, cn due to Viennot, Flajolet and Dumont^{9} and discovered around 1980. These all involve permutations restricted by a parity condition of sorts.

— Viennot [60] developed a model based on the differential system satisfied by the Jacobian elliptic functions, namely

x′ = yz,  y′ = −zx,  z′ = −k^{2}xy,

where x, y, z represent the classical Jacobian functions sn, cn, dn. This leads to an interpretation of the coefficients in terms of a class of permutations called by Viennot “Jacobi permutations”. Such permutations satisfy parity restrictions (mimicking the differential system) and are enumerated by Euler numbers.

— Flajolet [19] observed, from continued fraction theory, that the quantity

(−1)^{n}(2n)! · [z^{2n}α^{2k}] cn(z, α)

counts alternating permutations of length 2n that have k valleys of even value.

— Dumont [14, 15] provided an elegant interpretation of the coefficients of sn, cn, dn in terms of the parity of peaks of cycles in the cycle decomposition of permutations. Dumont’s results are based on consideration of Schett’s partial differential operator, already encountered in (32).

Theorem 4 adds another parity-based permutation model to the list. □

9These authors were answering a question of Marco Schützenberger in the 1970’s. Schützenberger first conjectured that Jacobian elliptic functions should have combinatorial content since their coefficients involve both factorial and Euler numbers.

5. Third model: permutations and repeated patterns

Our third combinatorial model of Dixonian functions is again in terms of permutations. It relies on repeated permutations, much in the style of Flajolet and Françon’s earlier interpretation of Jacobian (sn, cn) elliptic functions [21], but different. It is notable that the Flajolet-Françon permutations are based on a binary pattern, whereas those needed here involve a symmetry of order three. The type of an element (a value) in a permutation is, as in the previous section, any of peak, valley, double rise, double fall.

Definition 3. An r–repeated permutation of size n is a permutation such that for each j with j ≥ 0, the elements of value in {jr + 1, jr + 2, …, jr + r − 1} ∩ {1, …, n} are all of the same local order type, namely all peaks, valleys, double rises, or double falls.

For instance, the permutation

12 8 7 2 5 10 3 6 13 9 1 4 11_{↗↘ ↘↘ ↘↘ ↘↗ ↗↗ ↗↘ ↘↗ ↗↗ ↗↘ ↘↘ ↘↗ ↗↗ ↗↘}

bordered with (−∞, −∞) is a 3-repeated permutation of size 13, as is immediately verified by the listing of order types on the second line.

Dixonian functions will be proved to enumerate a variety of 3–repeated permutations. The proof is indirect and it first needs Flajolet’s combinatorial theory of continued fractions [19, 30] as well as a bijection between a system of weighted lattice paths and permutations, of which a first instance was discovered by Françon and Viennot in [27].

5.1. Combinatorial aspects of continued fractions. Define a lattice path, also known as a Motzkin path, of length n as a sequence of numbers s = (s_{0}, s_{1}, …, s_{n}), satisfying the conditions

s_{0} = s_{n} = 0, s_{j} ∈ ℤ_{≥0}, |s_{j+1} − s_{j}| ∈ {−1, 0, +1}.

This can be represented as a polygonal line in the Cartesian plane ℤ × ℤ. A step is an edge (s_{j}, s_{j+1}), and it is said to be (respectively) an ascent, a level, or a descent according to the value (respectively) +1, 0, −1, of s_{j+1} − s_{j}; the quantity s_{j} is called the (starting) altitude of the step. A path without level steps is a Dyck path. Motzkin paths are enumerated by Motzkin numbers; Dyck paths belong to the Catalan realm [56, pp. 219–229].

Let P (a, b, c) be the generating function of lattice paths in infinitely many indeterminates a = (a_{k}), b = (b_{k}), c = (c_{k}), with a_{k} marking an ascent from altitude k and similarly for descents marked by b_{k} and for levels marked by c_{k}. In other words, associate to each path ϖ a monomial 𝔪(ϖ), as in this example:

[figure: ϖ =; \begin{aligned}𝔪(ϖ) & = & a_{0}a_{1}c_{2}b_{2}a_{1}a_{2}b_{3}b_{2}b_{1}a_{0}b_{1} \\ & = & a_{0}^{2}a_{1}^{2}a_{2}b_{1}^{2}b_{2}^{2}b_{3}c_{1}. \end{aligned}]

One has

P (a, b, c) = ∑_{ϖ} 𝔪(ϖ),

where the sum is over all Motzkin paths ϖ.

Here is what Foata once called “the shallow Flajolet Theorem”^{10} taken from [19]:

Theorem 5 (Flajolet [19]). The generating function in infinitely many variables enumerating all Motzkin paths according to ascents, descents, levels, and corresponding altitudes is

(33)  P (a, b, c) = \frac{1}{1 − c_{0} − \frac{a_{0}b_{1}}{1 − c_{1} − \frac{a_{1}b_{2}}{1 − c_{2} − \frac{a_{2}b_{3}}{...}}}}.

Proof. All there is to it is the correspondence between a rational generating function given by a quasi-inverse and certain Motzkin paths of height ≤ 1. Figuratively, start from

\frac{1}{1 − c_{0} − a_{0}Y_{1}b_{1}} ≡ \{ [figure: c_{0}; +; a_{0}; Y_{1}; b_{1}] \}⋆ .

(These are Motzkin paths of height at most 1, with level steps at altitude 1 marked by Y_{1}.) It then suffices to apply repeatedly substitutions of the form

Y_{1} ↦ (1 − c_{1} − a_{1}Y_{2}b_{2})^{−1}, Y_{2} ↦ (1 − c_{2} − a_{2}Y_{3}b_{3})^{−1}, ⋯,

in order to generate the continued fraction of the statement. □

A bona fide generating function is obtained by making the generating function homogeneous, setting

a_{j} ↦ α_{j}z, b_{j} ↦ β_{j}z, c_{j} ↦ γ_{j}z

where α_{j}, β_{j}, γ_{j} are new formal variables. Once these variables are assigned numerical values (for instance α_{j} = β_{j} = γ_{j} = j + 1), the generating function of (33) enumerates weighted lattice paths by length (marked by z), with weights {α_{j}} ∪ {β_{j}} ∪ {γ_{j}} taken multiplicatively. For ease of reference, we encapsulate this notion into a formal definition.

Definition 4. A system of path diagrams is the class of multiplicatively weighted Motzkin paths determined by a possibility function Π, which assigns numerical values to the formal variables α_{j}, β_{j}, γ_{j}.

The values are normally integers, in which case a particular path diagram is equivalent to a lattice path augmented by a sequence of “choices”, the number of possible choices being α_{j} for an ascent from altitude j, and so on. Various systems of path diagrams are known to be bijectively associated to permutations, involutions, set partitions, and preferential arrangements, to name a few [19, 20, 27, 30].

Note 8. André’s method and an alternative derivation of the J-fractions relative to sm, cm. Désiré André published in 1877 a remarkable study [1], in which he was able to represent the coefficients of the Jacobian functions sn, cn as arising from weighted lattice paths, that is, from a system of path diagrams. We apply here André’s methodology to Dixonian functions.

10This designation stands to reason as the proof is extremely easy, so that the theorem borders on being an “observation”. However, the paper [19] is really a “framework” where orthogonal polynomial systems, lattice paths, continued fractions, Hankel determinants, etc, all find a combinatorial niche.

A direct use of the continued fraction Theorem 5, then provides a direct derivation of the J-fractions relative to sm, cm, which relies in a simple way on basic algebraic properties of the fundamental differential system (7).

Let us consider the case of smh, which is the s-component of the usual system s′ = c^{2}, c′ = s^{2}. From this system and elementary algebra induced by c^{3} − s^{3} = 1, we find that s satisfies a third-order non linear differential equation, s′′′ = 6s^{4} + 4s, and more generally (∂ represents differentiation with respect to the independent variable):

(34)  ∂^{3}s^{m} = (m(m + 1)(m + 2)s^{3} + 2m(m^{2} + 1) + \frac{m(m − 1)(m − 2)}{s^{3}}) s^{m}.

This shows that there exists a family of polynomials (P_{k}) with deg(P_{k}) = 3k + 1 such that

∂^{3k}s = P_{k}(s),

where P_{k} is of the form P_{k}(w) = wP̂_{k}(w^{3}), with P̂_{k} itself a polynomial. Then, Taylor’s formula provides

(35)  smh(z) = ∑_{k=1}^{∞} P′_{k}(0)\frac{z^{3k+1}}{(3k + 1)!}

upon taking into account the fact that only coefficients of index 1, 4, 7, … in smh are nonzero. Thus, the Taylor coefficients of smh are accessible from the coefficients of the lowest degree monomials in the P_{k} polynomials.

Introduce now a notation for the coefficients of the P_{k} polynomials:

P_{k,m} = [w^{m}]P_{k}(w).

The basic equation (34) implies the recurrence (m ≡ 1 (mod 3)):

P_{k+1,m} = (m + 1)(m + 2)(m + 3)P_{k,m+3} + 2m(m^{2} + 1)P_{k,m} + (m − 1)(m − 2)(m − 3)P_{k,m−3}.

This relation expresses precisely the fact that the coefficient [w^{m}]P_{k}(w) is the number of weighted paths starting at the point (0, 1) and ending at (k, m) in the lattice ℤ ⊕ (1 + 3ℤ), where the elementary steps are of the form a⃗ = (1, +3), b⃗ = (1, −3), c⃗ = (1, 0), and the weights of steps starting at an ordinate m ≡ 1 (mod 3) are respectively

(36)  α̂_{m} = m(m + 1)(m + 2), β̂_{m} = m(m − 1)(m − 2), γ̂_{m} = 2m(m^{2} + 1).

Up to a vertical translation of −1 and a vertical rescaling of the lattice by a factor of \frac{1}{3}, the coefficient P_{ν,1} is seen to enumerate standard weighted lattice paths (having steps (1, 1), (1, −1), and (1, 0)) of length ν with the new weights for steps starting at altitude ℓ being:

α_{ℓ} = (3ℓ + 1)(3ℓ + 2)(3ℓ + 3), β_{ℓ} = (3ℓ − 1)(3ℓ)(3ℓ + 1), γ_{ℓ} = 2(3ℓ + 1)((3ℓ + 1)^{2} + 1).

Theorem 5 then yields the continued fraction expansion of the ordinary generating function associated to smh (equivalently, to sm under a simple change of signs). Similar calculations apply to the other J-fractions of Section 2, including the one relative to cm.

It is especially interesting to observe the way in which the original differential system satisfied by smh, cmh churns out weighted lattice paths by means of higher order differential relations. This is the spirit of André’s work who played the original game on the Jacobi normal form of elliptic functions taken under the form (André’s notations):

[\frac{dφ(x)}{dx}]^{2} = D + Vφ^{2}(x) + Gφ^{4}.

Obviously, André intends D to mean “right” (“droite” in french), G to mean “left” (“gauche” in french) and V to mean “vertical”, as he had in mind paths rotated by 90^{∘}. (In one more step, André could have discovered the continued fractions associated to sn, cn had he known Theorem 5 or the earlier technique of “Stieltjes matrices” [57], whose determination is equivalent to a continued fraction expansion.) □

Note 9. On addition theorems. A standard way to derive explicit continued fraction expansions is by means of Rogers’ addition theorem (itself logically equivalent to a diagonalization technique of Stieltjes). We say that an analytic function defined near 0 (or a formal power series) f(z) satisfies an addition formula if it can be decomposed in terms of a sequence ϕ_{k} of functions as

(37)  f(z) = ∑_{k≥0} ϕ_{k}(x) · ϕ_{k}(y),  where ϕ_{k}(x) =_{x→0} O(x^{k}).

Then knowledge of the coefficients [x^{k}]ϕ_{k}(x) and [x^{k+1}]ϕ_{k}(x) implies knowledge of the coefficients in the J-fraction representation of the formal Laplace transform F of f, F (s) = L(f, s): the formulæ are simple, see Theorem 53.1 in Wall’s treatise [62, p. 203]. A good illustration is the addition formula for sec(z), namely

sec(x + y) = \frac{1}{cos x cos y − sin x sin y} = ∑_{k≥0} sec x tan^{k} x · sec y tan^{k} y,

which provides analytically the continued fraction expansion of the ordinary generating function of Euler numbers derived combinatorially below as (41).

Dixonian functions, being elliptic functions, are known to admit addition formulæ, albeit of a form different from (37). For instance, one has

cm(u + v) = \frac{c_{1}c_{2} − s_{1}s_{2}(s_{1}c_{2}^{2} + s_{2}c_{1}^{2})}{1 − s_{1}^{3}s_{2}^{3}}
with s_{1} = sm(u), c_{1} = cm(u), s_{2} = sm(v), c_{2} = cm(v),

found in Dixon’s paper (§33, Equation (39), p. 183 of [13]). It is in fact possible to combine this elliptic addition theorem with developments from the previous note in order to come up with an addition theorem of Stieltjes-Rogers type (37) relative to Dixonian functions: the form of the ϕ_{k} is essentially given by P_{k}(sm(z)) and the generating function for the P_{k} polynomials can be explicitly determined (details omitted). □

5.2. Correspondences between lattices paths and permutations. We next need bijections due to Françon-Viennot [27] between permutations and two systems of path diagrams.

It is convenient to start the discussion by introducing what V.I. Arnold [2] calls snakes. Consider piecewise monotonic smooth functions from ℝ to ℝ, such that all their critical values (i.e., the values of their maxima and minima) are different, and take the equivalence classes up to orientation preserving maps of ℝ × ℝ. It is sufficient to restrict attention to the two types,

(38)  _{(7 1 4 2 6 3 5)  (5 1 6 3 4 2)}

Figure 11. The sweepline algorithm: a snake and its associated Dyck path.

respectively called the (−∞, −∞) and (−∞, +∞) types. Clearly an equivalence class is an alternating permutation (see the second line of (38)). The exponential generating functions corresponding to the two types of (38) are by André’s classic theorem respectively,

tan(z) = \frac{sin(z)}{cos(z)}, sec(z) = \frac{1}{cos(z)},

size being measured by the number of critical points and equivalently the length of the corresponding permutation.

A simple sweepline algorithm associates to a snake a Dyck path as follows. Consider first the (−∞, −∞) case. Imagine moving a line from −∞ in the vertical direction towards +∞. At ordinates that are a midpoint between two successive critical values, associate a number which is the number of intersection points of the line with the curve. For a snake with 2n + 1 critical points, this gives us a sequence of numbers x = (x_{0}, x_{1}, …, x_{2n+2}) such that x_{j+1} − x_{j} = ±2, x_{0} = 2, and x_{2n+1} = 0. The rescaled sequence

ξ = (ξ_{0}, ξ_{1}, …, ξ_{2n+1}), ξ_{j} := (x_{j} − 2)/2,

becomes a standard Dyck path (Figure 11), which obviously depends only on the underlying alternating permutation.

Alternating permutations are far more numerous than Dyck paths, so that they must be supplemented by additional information in order to obtain a proper encoding. To this effect, introduce the system of path diagrams given by the possibility rule

(39)  Π^{odd} : Π^{(−∞,−∞)} :  α_{j} = (j + 1), β_{j} = (j + 1), γ_{j} = 0,

which, we claim, is now bijectively associated to odd-length alternating permutations,

−∞ < σ_{1} > σ_{2} < ⋯ < σ_{2m+1} > −∞ .

This is easily understood as follows. When executing the sweepline algorithm, there are several places (possibilities) at which a local maximum (or peak) and a local minimum (or valley) can be attached by “capping” or “cupping”:

These possibilities, as a function of altitude (the ξ_{j}) are seen to correspond exactly to the possibility set (39). (For instance, on the example, the rescaled Dyck path is at altitude 2, there are three possibilities for cupping and three possibilities for capping.)

For the (−∞, +∞) case, that is, even-length alternating permutations, a similar reasoning shows that the possibilities are

(40)  Π^{(−∞,+∞)} :  α_{j} = (j + 1), β_{j} = j, γ_{j} = 0,

See:

By the continued fraction theorem, these two encodings yield two continued fraction expansions originally discovered by Stieltjes (the Laplace transforms are taken as “formal” here):

(41)  \{ \begin{matrix}∫_{0}^{∞} tan(zt)e^{−t} dt & = & \frac{z}{1 − \frac{1 · 2 z^{2}}{1 − \frac{2 · 3 z^{2}}{...}}}, \\ ∫_{0}^{∞} sec(zt)e^{−t} dt & = & \frac{1}{1 − \frac{1^{2} z^{2}}{1 − \frac{2^{2} z^{2}}{...}}}. \end{matrix}

The previous two bijections can be modified so as to take into account all permutations, not just alternating ones. What this corresponds to is spotted snakes, which are snakes augmented with an arbitrary finite number of non-critical points that are distinguished. (As usual, one operates up to topological equivalence and the spotted points must have different altitudes.) It then suffices to encode nodes on upward and downward slopes by level steps of the system of path diagrams in order to get Motzkin paths [19, 27, 30]. The path diagrams so obtained have closely resembling possibility Figure 12. A figure suggesting that the Françon–Viennot correspondence and the sweepline algorithms are one and the same thing.

rules:

\{ \begin{matrix}\overline{Π}^{(−∞,−∞)} : & α_{j} = (j + 1), & β_{j} = (j + 1), & γ_{j} = 2j + 2 \\ \overline{Π}^{(−∞,+∞)} : & α_{j} = (j + 1), & β_{j} = j, & γ_{j} = 2j + 1. \end{matrix}

(In the (−∞, +∞) case, one can never insert a cap or a descending node on the extreme right, so that one possibility is suppressed for descents and level steps of the path diagrams.) The resulting continued fractions are then

(42)  \{ \begin{matrix}∑_{n=0}^{∞}(n + 1)!z^{n+1} & = & \frac{z}{1 − 2z − \frac{1 · 2 z^{2}}{1 − 4z − \frac{2 · 3 z^{2}}{...}}}, \\ ∑_{n=0}^{∞} n!z^{n} & = & \frac{1}{1 − z − \frac{1^{2} z^{2}}{1 − 3z − \frac{2^{2} z^{2}}{...}}}. \end{matrix}

These last two fractions are originally due to Euler. (The last four continued fraction expansions were first established combinatorially by Flajolet in [19].)

Note 10. Snakes of bounded width. The usual enumeration of snakes up to deformation is usually presented by the Russian School as resulting from the “Seidel-Entringer-Arnold triangle”. The continued fraction connection exposed here gives access to new parameters, and, in particular, the ones associated with convergents of a basic continued fraction. In this context, it provides the ordinary generating function of odd snakes of bounded width, where width is defined as the maximal cardinality of the image of any value:

width(s) := max_{y∈ℝ}card {x | s(x) = y} .

(Width is a trivial variant of the clustering index introduced in [19, p. 159].) Take for definiteness snakes with boundary condition (−∞, +∞). Transfer matrix methods first imply that the generating function W^{[h]}(z) of snakes having width at most 2h − 1 is a priori a rational function,

W^{[h]}(z) = \frac{\overline{P}_{h}(z)}{\overline{Q}_{h}(z)},

the first few values being, for width at most 1, 3, 5, 7,

\frac{1}{1}, \frac{1}{1 − z^{2}} , \frac{1 − 4z^{2}}{1 − 5z^{2}} , \frac{1 − 13z^{2}}{1 − 14z^{2} + 9z^{4}} .

(For instance, width ≤ 3 corresponds to the type (−∞, 2, 1, 4, 3, 6, 5 …, +∞).) Continued fraction theory [19] further implies that W^{[h]}(z) is a convergent of the continued fraction relative to sec(z) in (41), with the denominators \overline{Q}_{h}(z) being reciprocal polynomials of an orthogonal polynomial system. Here, one finds

\overline{Q}_{h}(z) = z^{h}Q_{h}(1/z), where Q_{h}(z) = [t^{h}](1 + t^{2})^{−1/2} exp(z arctan t),

and the Q_{h} are the Meixner polynomials [11]. (This can be verified via the generating function of the Q_{h}, which satisfies a differential equation of order 1; the method also adapts to the determination of the associated P_{h} polynomials.) □

Note 11. On Françon-Viennot. The bijective encoding of permutations by path diagrams as presented here is exactly the same as the one obtained from the original Françon-Viennot correspondence—only our more geometric presentation differs. Indeed, one may think of the Françon–Viennot^{11} correspondence given in [27] as the gradual construction of an increasing binary tree, upon successively appending nodes at dangling links [19, 30]; see also Figure 12, where both a snake and the underlying tree are represented. Use will made below of the tree view of the Françon-Viennot correspondence. (Biane [6] discovered a related sweepline algorithm, based on the decomposition of permutations into cycles.) □

5.3. The models of r–repeated permutations. The Françon-Viennot bijection and the Continued Fraction Theorem provide:

Proposition 5. Let R_{rn+1} be the number of r–repeated permutations of length rn+1 bordered by (−∞, −∞), and let R_{rn}^{⋆} be the number of r–repeated permutations of length rn bordered by (−∞, +∞) The corresponding ordinary generating functions admit a continued fraction expansions of the Jacobi type,

\begin{matrix}∑_{ν≥0} R_{rν+1}z^{rν+1} = \\ \frac{z}{1 − 2 · 1^{r} z^{r} − \frac{1 · 2^{2} ⋯ r^{2} · (r + 1) · z^{2r}}{1 − 2 · (r + 1)^{r} z^{r} − \frac{(r + 1) · (r + 2)^{2} ⋯ (2r)^{2} · (2r + 1) · z^{2r}}{...}}}, \end{matrix}

11Viennot (private communication, February 2006) comments that, when the Françon-Viennot paper was submitted, the formulation of the bijection between trees and path diagrams was very much like our description, but an anonymous referee insisted that all references to binary trees be deleted and be replaced by descriptions in terms of words. (Twenty years later, trees have become fully acceptable objects in the combinatorial literature!) For a tree version of the correspondence, Viennot also refers to his lecture notes: “Une théorie combinatoire des polynômes orthogonaux”, 221 pages, Publication du LACIM (1983), UQAM, Montréal, of which an accessible summary is available [61].

\begin{matrix}∑_{ν≥0} R_{rν}^{⋆}z^{rν} = \\ \frac{1}{1 − (0^{r} + 1^{r})z^{r} − \frac{1^{2} · 2^{2} ⋯ r^{2} · z^{2r}}{1 − (r^{r} + (r + 1)^{r})z^{r} − \frac{(r + 1)^{2} · (r + 2)^{2} ⋯ (2r)^{2} · z^{2r}}{...}}}, \end{matrix}

where the numerators are of degree 2r and the denominators are of degree r in the depth index.

Proof. In accordance with Note 11 and Figure 12, it is best to regard the Françon– Viennot correspondence as the inductive construction of an increasing binary tree. If at some stage there are ℓ dangling links, then there are several cases to be considered for adding a node depending on its type. Here is a table giving: the types of elements in a permutation; the types of nodes in the tree; the number of possibilities.

Perm.: peak  valley  double rise  double fall
Tree: leaf double node left branching right branching
Poss.: ℓ  ℓ  ℓ  ℓ
ℓ − 1  ℓ  ℓ − 1  ℓ

The third line corresponds to allowing the largest (and last inserted) value rn+1 to fall anywhere, the permutation being bordered by (−∞, −∞). The fourth line corresponds to a permutation to be such that rn + 1 occurs at the end, that is, a permutation bordered by (−∞, +∞). The history of the tree construction corresponds to Motzkin paths that start at altitude 1, have steps that are grouped in batches of r, all of the same type within a batch. This is converted into a standard Motzkin path by a succession of two operations: (i) shift the path down by 1; (ii) then divide the altitude by r, so that sequences of steps of a single type inherit weights multiplicatively. (This is essentially the argument of [21].) The statement results. □

The close resemblance between the case r = 3 of Proposition 5 and Theorem 2 is striking. Notice however that it is required to adjust the possibility function for level steps, i.e., correct the denominators. This is achieved by means of “polarization”:

Definition 5. A 3-repeated permutation of length 3ν + 1 is said to be polarized if some (possibly none, possibly all) of the consecutive factors in the word representation of the permutation that are of the form 3j +3, 3j +2, 3j +1 or of the form 3j +1, 3j +2, 3j +3 are marked.

We shall use a minus sign as a mark and write ¯3 ≡ −3. For instance, two polarized 3-repeated permutation of size 19 are

14 2 17 8 19 3 5 6 16 12 11 10 7 13 1 4 18 9 15 14 2 17 8 19 3 5 6 16 \overline{12} \overline{11} \overline{10} 7 13 1 4 18 9 15

(Only one factor is amenable to marking here.) We state:

Theorem 6. The exponential generating function of 3–repeated polarized permutations bordered by (−∞, −∞) is

smh(z).

Proof. By Conrad’s fractions, Françon–Viennot, and the reasoning that underlies Proposition 5. Polarization has been introduced on purpose in order to add to the number of possibilities: it contributes 2(3ℓ + 1) further possibilities to the Motzkin path when we are at altitude ℓ, corresponding to a (polarized) sequence of three consecutive one-way branching nodes of two possible types (left- or right-branching) attached to one of the 3ℓ + 1 dangling links. □

The notion of polarization can be similarly introduced to interpret the coefficients of cmh(z) (details omitted). This theorem nicely completes the picture of the relation between r-repeated permutations and special functions. r = 1: We are dealing with unconstrained permutations. In this case, Proposition 5

reduces to the two continued fraction expansions relative to ∑ n! z^{n}. Note that marking, with an additional variable t, the number of rises (i.e., double rises and valleys), leads to a continued fraction expansion of the bivariate ordinary generating function of Eulerian numbers [19], which was found analytically by Stieltjes [58].

r = 2: We are dealing with the doubled permutations of Flajolet-Françon [21]. The permutations are enumerated by Euler numbers, and when rises are taken into account, one obtains a bivariate ordinary generating function that is the Laplace transform of the Jacobian elliptic functions sn(z, t), cn(z, t).

r = 3: The generating functions are related to sm, cm, the correspondence being exact

when polarization is introduced. The Dixonian functions involve once more a third-order symmetry that is curiously

evocative of the fact that they parametrize the Fermat cubic. It would be of obvious interest (but probably difficult) to identify which special functions are associated with higher order symmetries corresponding to r ≥ 4.

6. Further connections

Our goal in this article has been to demonstrate that the Dixonian parametrization of the Fermat cubic has interesting ramifications in several different fields. The way these functions have largely independently surfaced in various domains is striking: occurrences now known include the theory of continued fractions and orthogonal polynomials, special functions (e.g, Lundberg’s hypergoniometric functions), combinatorial analysis (the elementary combinatorics of permutations), and a diversity of stochastic processes (special urn models and branching processes, but also birth and death processes). We now briefly discuss other works, some very recent, which confirm that Dixonian functions should indeed be considered as part of the arsenal of special functions.

6.1. Jacobi and Weierstraß forms. As it is well known, there are three major ways of introducing elliptic functions [63], namely, by way of the Jacobian functions sn, cn, dn (defined from inverses of an Abelian integral over a curve y^{2} = P_{4}(x)), by their Weierstraß form ℘ (associated to a curve y^{2} = P_{3}(x)), and by theta functions. The reductions of sm, cm to normal form are no surprise since they are granted by general theorems. What stands out, however, is the simplicity of the connections, which are accompanied in one case by further combinatorial connections.

6.1.1. The Jacobian connection. This is the one observed by Cayley in his two page note [8] and already alluded to in Section 1. Cayley finds the parametrization of Fermat’s cubic x^{3} + y^{3} = 1 (locally near (−1, 2^{1/3})) in the form (x, y) = (ξ(u), η(u)), where

ξ(u) = \frac{−1 + θscd}{1 + θscd}, η(u) = \frac{2^{1/3}(1 + θ^{2}s^{2})^{2}}{1 + θscd} , θ := 3^{\frac{1}{4}} e^{5iπ/12}

and here s ≡ sn(u), c ≡ cn(u), d ≡ dn(u) are Jacobian elliptic functions for the modulus k := e^{5iπ/6}. (This corrects what seem to be minor errors in Cayley’s calculations.) Cayley’s calculations imply an expression of sm, cm in terms of sn, cn, dn as simple variants of ξ, η (with π_{3} as in (12)):

sm(z) = ξ (\frac{z + π_{3}/6}{2^{1/3}θ}) , cm(z) = η (\frac{z + π_{3}/6}{2^{1/3}θ}) .

For an alternative approach based on Lundberg’s hypergoniometric functions, see [42, Sec. 5].

6.1.2. The Weierstraß connection. The relations appear to be more transparent than with Jacobian functions. As a simple illustration, consider the function

P (z) := smh(z) · cmh(z).

From the basic differential equation framework, this function satisfies the ordinary differential equations,

P′′ = 6P^{2}, P^{′2} = 4P^{3} + 1.

Thus, up to a shift of the argument, P (z) is a Weierstraß ℘ with parameters g_{2} = 0, g_{3} = −1 (corresponding to the usual hexagonal lattice). While ℘(0) = ∞, the initial conditions are here P (z) = z + O(z^{2}). In this way we find (see (12))

P (z) = P(z − ζ_{0}; 0; −1), ζ_{0} = \frac{2}{3}π_{3} = \frac{1}{3π} Γ (\frac{1}{3})^{3} .

From these calculations, there results that P = smh · cmh is implicitly defined as the solution of the equation

(43)  ∫_{0}^{Y} \frac{dw}{\sqrt{1 + 4w^{3}}} ≡ Y ·_{2}F_{1} [\frac{1}{3}, \frac{1}{2}; \frac{4}{3}; −4Y^{3}] = z.

Thus, the continued fraction expansion of sm · sm found in Section 2 can be re-expressed as follows:

Proposition 6. The Laplace transform of the compositional inverse of the function Y ·_{2}F_{1} [\frac{1}{3}, \frac{1}{2}; \frac{4}{3}; −4Y^{3}], equivalently of the special ℘(z; 0, −1) expanded near its real zero ζ_{0} := \frac{2}{3}π_{3}, admits a continued fraction expansion with sextic numerators and cubic denominators.

(The hypergeometric parameters differ from those of Proposition 1. This function otherwise constitutes a good case of application of André’s method of Note 8.)

In passing, we remark that P (z) serves to express the fundamental function (the ψ-function) of the urn model defined by the matrix

T_{23} = ( \begin{matrix}−2 & 3 \\ 4 & −3 \end{matrix} ) .

This model (corresponding to case A in Figures 7 and 8) is of interest as it describes the fringe behaviour of 2–3 trees and other locally balanced trees [22, 47]. The relation is (see [22, p. 1223] for the definition of ψ)

ψ(z) = 2^{2/3}P (2^{1/3}z) = 2^{2/3} smh(2^{1/3}z) cmh(2^{1/3}z),

as can be verified either from differential relations or from the available ℘ forms. Thus yet another combinatorial model of Dixonian functions is available in terms of histories of the T_{23} urn, corresponding to the rewrite rules

xx ⟶ yyy, yyy ⟶ xxxx.

In other words:

Proposition 7. The T_{23} urn model can be described in terms of Dixonian functions via the product P (z) = sm(z) · cm(z). In particular, the number of histories of the T_{23} urn initialized with xx such that at time 3ν + 1 all balls are of type y is

H_{3ν+1,3ν+3} = (3ν + 1)!2^{ν+1} [z^{3ν+1}] smh(z) · cmh(z).

In an unpublished manuscript Dumont [17] also observes (by way of singularities) that the function

Q(z) := \frac{sm(z)}{3(1 − cm(z))},

which now has a double pole at 0 is none other than ℘(z; 0, \frac{1}{27}). From this and similar considerations, Dumont then obtains in a simple way the identities

cm(z) = \frac{3℘′(z) + 1}{3℘′(z) − 1},  sm(z) = \frac{6℘(z)}{1 − 3℘′(z)},  ℘(z) := ℘(z, 0, \frac{1}{27}).

6.2. Laplace transforms of elliptic functions. Transforms of elliptic functions are only considered sporadically in the literature, usually in the context of continued fraction theory. Some of them are explicitly known in the classical case of Jacobian elliptic functions, which are somehow close to sm, cm given the developments of Section 6.1.1.

Stieltjes, followed by Rogers, was the first to determine continued fraction expansions for the Laplace transforms of the three fundamental Jacobian functions, sn, cn, dn. Some of these expansions were rediscovered by Ramanujan who made the Laplace transforms explicit. For instance, following Perron’s account (pp. 134–135 and 219 of [48]) and Berndt’s edition of Ramanujan’s Notebooks [5, p. 163], one has

(44)  \begin{aligned}∫_{0}^{∞} cn(tx; k)e^{−t} dt & = & \frac{2π}{kK} ∑_{ν=0}^{∞} \frac{q^{ν+ \frac{1}{2}}}{1 + q^{2ν+1}} \frac{1}{1 + x^{2} (\frac{(2ν+1)π}{2K})^{2}} \\ & = & \frac{1}{1 + \frac{1^{2}x^{2}}{1 + \frac{2^{2}k^{2}x^{2}}{1 + \frac{3^{2}x^{2}}{1 + ...}}}}. \end{aligned}

(The numerators are of the alternating form 1^{2}, 2^{2}k^{2}, 3^{2}, 4^{2}k^{2}, 5^{2}x^{2}, …; usual notations from elliptic function theory are employed [63].) The calculation of the Laplace transform is effected via the Fourier expansion of cn, which goes back to Jacobi and is detailed in [63, §22.6]. It would be of obvious interest to carry out calculations and determine in which class of special functions the Laplace transforms of the Dixonian functions live.

Note otherwise that the Hankel determinant evaluations stemming from (44) lie at the basis of Milne’s results [46] regarding sums of squares, while the continued fraction expansion is related to several permutation models of Jacobian elliptic functions that were briefly mentioned in Section 5.

6.3. Orthogonal polynomial systems. As it is well known each continued fraction of type J is associated to an orthogonal polynomial system (OPS), which provides in particular the denominators of the fraction’s convergents [11, 48, 62]. The OPS arising from the Laplace transforms of Jacobian functions have a long tradition, starting with Stieltjes who first found an explicit representation of the orthogonality measure. See the study by Ismail, Valent, and Yoon [34] for a recent perspective and pointers to the older literature as well as the book by Lomont and Brillhart [43] entirely dedicated to “elliptic polynomials”.

It is a striking fact that the orthogonal polynomial systems related to the continued fraction expansions of Section 2 have very recently surfaced in independent works of Gilewicz, Leopold, Ruffing, and Valent [28, 29]. These authors are motivated by the classification of (continuous-time) birth and death processes^{12}, especially the ones whose rates are polynomials in the size of the population. Here the birth and death rates are cubic.

From the rather remarkable calculations of [28], one gets in particular [do c → 0 and μ_{0} = 0] the generating function of a family of polynomials,

(45)  \begin{aligned}∑_{n≥0} Q_{n}(z)\frac{t^{3n}}{(3n)!} & = & (1 − t^{3})^{−1/3}E_{3,0} (zθ(t)^{3}) \\ E_{3,0}(y) & := & ∑_{n≥0} \frac{y^{n}}{(3n)!},  θ(t) := ∫_{0}^{t} \frac{dw}{(1 − w^{3})^{2/3}} . \end{aligned}

12Thanks to works of Karlin and McGregor in the 1950’s, many stochastic characteristics of a process can be described in terms of a continued fraction, its OPS, and its orthogonality measure(s); see for instance [23] for a recent review.

(This OPS resembles a hybrid of the Meixner and Brenke classifications of orthogonal polynomials [11, Ch. V].) The first few values, Q_{0}, Q_{1}, Q_{2}, Q_{3}, Q_{4}, are

\begin{matrix}1, 2 + z, 160 + 100 z + z^{2}, 62720 + 42960 z + 672 z^{2} + z^{3}, \\ 68992000 + 49755200 z + 963600 z^{2} + 2420 z^{3} + z^{4}. \end{matrix}

These polynomials are the monic versions of the denominators of the convergents in the continued fraction

\frac{1}{z + 2 − \frac{36}{z + 98 − \frac{14400}{z + 572 − ...}}},

which, up to normalization, represents the Laplace transform of cm(z). The occurrence in (45) of the fundamental Abelian integral from which sm, cm are defined is especially striking. As shown in [28, 29], the moment problem is indeterminate; see also Ismail’s recent book [32] for a perspective.

Conclusion. We have studied a class of continued fraction with sextic numerators and cubic denominators and shown that, thanks to several converging works, these are endowed with a rich set of properties. This suggests, regarding the Stieltjes-Apéry fractions described in the introduction, that it would be of great interest to determine whether these are similarly endowed with a rich structure as regards combinatorics, special functions, and orthogonality relations. The question of investigating continued fractions whose coefficients are polynomials of higher degrees is a tantalizing one, but it is likely to be very difficult.

Acknowledgements. Eric Conrad is deeply indebted to Dominique Dumont (Université Louis Pasteur, Strasbourg) for referring him to the paper of A. C. Dixon that forms a centerpiece of his PhD dissertation [12]. Both authors also express their gratitude to Dumont for many ensuing discussions relative to this work. We are thankful to Hjalmar Rosengren for first pointing us to the works of Lindqvist and Peetre—Peter Lindqvist then very kindly offered a wealth of mathematical and historical material relative to Lundberg’s hypergoniometric functions. We are grateful to Galliano Valent for sharing with us an unpublished version of his recent work. Finally, Xavier Viennot and Christian Krattenthaler both provided insightful remarks on the preprint—thanks to both!

References

[1] Désiré André, Développements en séries des fonctions elliptiques et de leurs puissances, Annales scientifiques de l’École Normale Supérieure VI (1877), 265–328, Available freely from the site Numdam.

[2] Vladimir I. Arnold, Nombres d’Euler, de Bernoulli et de Springer pour les groupes de Coxeter et les espaces de morsification: le calcul des serpents, Leçons de mathématiques d’aujourd’hui (Éric Carpentier and Nicolas Nikolski, eds.), Cassini, Paris, 2000, pp. 61–98.

[3] Krishna B. Athreya and Peter E. Ney, Branching processes, Springer-Verlag, New York, 1972, Die Grundlehren der mathematischen Wissenschaften, Band 196.

[4] Bruce C. Berndt, Ramanujan’s notebooks, part II, Springer Verlag, 1989.

[5] , Ramanujan’s notebooks, part III, Springer Verlag, 1991.

[6] Philippe Biane, Permutations suivant le type d’excedance et le nombre d’inversions et interprétation combinatoire d’une fraction continue de Heine, European Journal of Combinatorics 14 (1993), 277–284.

[7] Mireille Bousquet-Mélou and Marko Petkovšek, Walks confined in a quadrant are not always D-finite, Theoretical Computer Science 307 (2003), no. 2, 257–276.

[8] Arthur Cayley, Reduction of ∫ \frac{dx}{(1−x^{3})^{2/3}} to elliptic integrals, Messenger of Mathematics XI (1882), 142–143, Reprinted in Collected Mathematical Papers.

[9] Brigitte Chauvin, Thierry Klein, Jean-Francois Marckert, and Alain Rouault, Martingales and profile of binary search trees, Preprint, University of Versailles, 2004, Available from

Arχiv:math.PR/0410211.

[10] Brigitte Chauvin and Alain Rouault, Connecting Yule process, bisection and binary search tree via martingales, Journal of the Iranian Statistical Society 3 (2004), no. 2, 89–116, Available from

Arχiv:math.PR/0410211.

[11] T. S. Chihara, An introduction to orthogonal polynomials, Gordon and Breach, New York, 1978.

[12] Eric van Fossen Conrad, Some continued fraction expansions of Laplace transforms of elliptic functions, PhD Thesis, The Ohio State University, 2002.

[13] Alfred Cardew Dixon, On the doubly periodic functions arising out of the curve x^{3}+y^{3}−3αxy = 1, The Quarterly Journal of Pure and Applied Mathematics 24 (1890), 167–233.

[14] Dominique Dumont, A combinatorial interpretation for the Schett recurrence on the Jacobian elliptic functions, Mathematics of Computation 33 (1979), 1293–1297.

[15] , Une approche combinatoire des fonctions elliptiques de Jacobi, Advances in Mathematics 1 (1981), 1–39.

[16] , Pics de cycles et dérivées partielles, Séminaire Lotharingien de Combinatoire 37 (1986), Art. B13a, 19 pp. (electronic).

[17] , Le paramétrage de la courbe d’équation x^{3} + y^{3} = 1, Manuscript, 1988, 18 pages.

[18] , Grammaires de William Chen et dérivations dans les arbres et arborescences, Séminaire Lotharingien de Combinatoire 37 (1996), Art. B37a, 21 pp. (electronic).

[19] Philippe Flajolet, Combinatorial aspects of continued fractions, Discrete Mathematics 32 (1980), 125–161.

[20] , On congruences and continued fractions for some classical combinatorial quantities, Discrete Mathematics 41 (1982), 145–153.

[21] Philippe Flajolet and Jean Françon, Elliptic functions, continued fractions and doubled permutations, European Journal of Combinatorics 10 (1989), 235–241.

[22] Philippe Flajolet, Joaquim Gabarró, and Helmut Pekari, Analytic urns, Annals of Probability 33 (2005), no. 3, 1200–1233, Available from Arχiv:math.PR/0407098.

[23] Philippe Flajolet and Fabrice Guillemin, The formal theory of birth-and-death processes, lattice path combinatorics, and continued fractions, Advances in Applied Probability 32 (2000), 750–778.

[24] Philippe Flajolet and Robert Sedgewick, Analytic combinatorics, Cambridge University Press, November 2004, Chapters I–IX of a book to be published, 609p.+x, available electronically from P. Flajolet’s home page.

[25] Dominique Foata and Marcel-P. Schützenberger, Théorie géométrique des polynômes Euleriens, Lecture Notes in Mathematics, vol. 138, Springer Verlag, 1970.

[26] Jean Françon, Histoires de fichiers, RAIRO Informat. Théor. 12 (1978), no. 1, 49–62.

[27] Jean Françon and Gérard Viennot, Permutations selon leurs pics, creux, doubles montées et doubles descentes, nombres d’Euler et de Genocchi, Discrete Mathematics 28 (1979), 21–35.

[28] Jacek Gilewicz, Elie Leopold, Andreas Ruffing, and Galliano Valent, Some cubic birth and death processes and their related orthogonal polynomials, Preprint, 2005, 17 pages.

[29] Jacek Gilewicz, Elie Leopold, and Galliano Valent, New Nevanlinna matrices for orthogonal polynomials related to cubic birth and death processes, Journal of Computational and Applied Mathematics 178 (2005), no. 1-2, 235–245.

[30] Ian P. Goulden and David M. Jackson, Combinatorial enumeration, John Wiley, New York, 1983.

[31] Hsien-Kuei Hwang, On convergence rates in the central limit theorems for combinatorial structures, European Journal of Combinatorics 19 (1998), no. 3, 329–343.

[32] Mourad E. H. Ismail, Classical and quantum orthogonal polynomials in one variable, Encyclopedia of Mathematics and its Applications, no. 98, Cambridge University Press, 2005.

[33] Mourad E. H. Ismail and D. Masson, Some continued fractions related to elliptic functions, Continued Fractions: from Analytic Number Theory to Constructive Approximation (Providence, RI) (B. Berndt and F. Gesztesy, eds.), Contemporary Mathematics, vol. 236, AMS, 1999, pp. 149–166.

[34] Mourad E. H. Ismail, Galliano Valent, and Gang J. Yoon, Some orthogonal polynomials related to elliptic functions, Journal of Approximation Theory 112 (2001), no. 2, 251–278.

[35] Svante Janson, Functional limit theorems for multitype branching processes and generalized Pólya urns, Stochastic Processes and Applications 110 (2004), no. 2, 177–245.

[36] Norman L. Johnson and Samuel Kotz, Urn models and their application, John Wiley, 1977.

[37] Donald E. Knuth, The art of computer programming, 3rd ed., vol. 1: Fundamental Algorithms, Addison-Wesley, 1997.

[38] Maxim Kontsevich and Don Zagier, Periods, Mathematics unlimited—2001 and beyond, Springer, Berlin, 2001, pp. 771–808.

[39] Serge Lang, Introduction to algebraic and abelian functions, second ed., Graduate Texts in Mathematics, vol. 89, Springer-Verlag, New York, 1982.

[40] Pierre Leroux and Xavier G. Viennot, Combinatorial resolution of systems of differential equations. IV: Separation of variables, Discrete Mathematics 72 (1988), 237–250.

[41] Peter Lindqvist and Jaak Peetre, Generalized trigonometric functions (solution to problem 10744), The American Mathematical Monthly 108 (2001), no. 5, 473–474.

[42] , Two remarkable identities, called twos, for inverses to some Abelian integrals, The American Mathematical Monthly 108 (2001), no. 5, 403–410.

[43] J. S. Lomont and John Brillhart, Elliptic polynomials, Chapman & Hall/CRC, Boca Raton, FL, 2001.

[44] Erik Lundberg, Om hypergoniometriskafunktioner af komplexa variabla, Manuscript, 1879, English translation, “On hypergoniometric functions of complex variables” available from Jaak Peetre’s web page.

[45] Hosam Mahmoud, Urn models and connections to random trees: a review, Journal of the Iranian Mathematical Society 2 (2003), 53–114.

[46] Stephen C. Milne, Infinite families of exact sums of squares formulas, Jacobi elliptic functions, continued fractions, and Schur functions, Developments in Mathematics, vol. 5, Kluwer Academic Publishers, Boston, MA, 2002, Reprinted from Ramanujan J. 6 (2002), no. 1, 7–149.

[47] Alois Panholzer and Helmut Prodinger, An analytic approach for the analysis of rotations in fringe-balanced binary search trees, Annals of Combinatorics 2 (1998), 173–184.

[48] Oskar Perron, Die Lehre von der Kettenbrüchen, vol. 2, Teubner, 1954.

[49] Nicolas Pouyanne, An algebraic approach of large Pólya processes, Submitted, 2005, 40 pages. (Extended abstract in Discrete Mathematics and Theoretical Computer Science, 2005, in press).

[50] Vincent Puyhaubert, Modèles d’urnes et phénomènes de seuil en combinatoire analytique, PhD thesis, École Polytechnique, Palaiseau, France, 2005.

[51] Tanguy Rivoal, Irrationalité d’au moins un des neuf nombres ζ(5), ζ(7), …, ζ(21), Acta Arithmetica 103 (2002), no. 2, 157–167.

[52] L. J. Rogers, On the representation of certain asymptotic series as convergent continued fractions, Proceedings of the London Mathematical Society (Series 2) 4 (1907), 72–89.

[53] Alois Schett, Properties of the Taylor series expansion coefficients of the Jacobian elliptic functions, Mathematics of Computation 30 (1976), no. 133, 143–147.

[54] N. J. A. Sloane, The On-Line Encyclopedia of Integer Sequences, 2000, Published electronically

    at http://www.research.att.com/~njas/sequences/.

[55] Richard P. Stanley, Enumerative combinatorics, vol. I, Wadsworth & Brooks/Cole, 1986.

[56] , Enumerative combinatorics, vol. II, Cambridge University Press, 1998.

[57] T. Stieltjes, Sur la réduction en fraction continue d’une série procédant selon les puissances descendantes d’ une variable, Annales de la Faculté des Sciences de Toulouse 4 (1889), 1–17.

[58] Thomas Jan Stieltjes, Œuvres complètes, Springer Verlag, 1993, Edited by Gerrit van Dijk.

[59] Alfred van der Poorten, A proof that Euler missed … Apéry’s proof of the irrationality of ζ(3), Mathematical Intelligencer 1 (1979), 195–203.

[60] Gérard Viennot, Une interpretation combinatoire des coefficients des développements en série entière des fonctions elliptiques de Jacobi, Journal of Combinatorial Theory, Series A 29 (1980), 121–133.

[61] , A combinatorial theory for general orthogonal polynomials with extensions and applications, Orthogonal polynomials and applications (Bar-le-Duc, 1984), Lecture Notes in Math., vol. 1171, Springer, Berlin, 1985, pp. 139–157.

[62] H. S. Wall, Analytic theory of continued fractions, Chelsea Publishing Company, 1948.

[63] E. T. Whittaker and G. N. Watson, A course of modern analysis, fourth ed., Cambridge University Press, 1927, Reprinted 1973.

E.C.: Department of Mathematics, The Ohio State University, 231 W 18th Avenue, Columbus, OH 43210 USA

  E-mail address: econrad@math.ohio-state.edu

P.F.: Algorithms Project, INRIA Rocquencourt, F-78153 (France)

  E-mail address: Philippe.Flajolet@inria.fr