INTEGRAL SOLUTIONS OF APÉRY-LIKE RECURRENCE EQUATIONS

Don Zagier

Abstract. In [4], Beukers studies the differential equation

((t^{3} + at^{2} + bt)F′(t))′ + (t − λ)F (t) = 0 ,  (*)

where a, b and λ are rational parameters, and asks for which values of these parameters this equation has a solution in ℤ[[t]], the motivating example being the Apéry sequence with a = 11, b = −1, λ = −3. We describe a search over a suitably chosen domain of 100 million triples (a, b, λ). In this domain there are 36 triples yielding integral solutions of (*). These can be further subdivided into members of four infinite classes, two of which are degenerate special cases of the other two, and seven sporadic solutions. Of these solutions, twelve, including all the sporadic ones, have parametrizations of Beukers type in terms of modular forms and functions. These solutions are related to elliptic curves over ℙ^{1} with four singular fibres.

1. Beukers’s recurrence equation. Beukers [4] considers the differential equation

(tP (t)F′(t))′ + (t − λ)F (t) = 0  (1)

with P (t) a quadratic polynomial, which we can take to have the form t^{2} + at + b. This equation has

a unique solution which is regular at the origin with F (0) = 1, given by F (t) = ∑_{n=0}^{∞} u_{n}t^{n} with u_{0} = 1

and

b(n + 1)^{2}u_{n+1} + (an^{2} + an − λ)u_{n} + n^{2}u_{n−1} = 0  (n ≥ 0)  (2)

so that u_{1} = λ/b, u_{2} = (λ^{2}−2aλ+b)/4b^{2}, etc. We are interested in finding values of (a, b, λ) for which F (Dt) ∈ ℤ[[t]] or D^{n}u_{n} ∈ ℤ (∀n ≥ 0) for some D ∈ ℕ. After a rescaling (a, b, λ) ↦ (a/D, b/D^{2}, λ/D), u_{n} ↦ D^{n}u_{n} we can assume D = 1 or u_{n} ∈ ℤ, and in future we shall assume after such a rescaling that the u_{n} are all integral, that there is no D > 1 such that D^{n}|u_{n} for all n, and that u_{1} ≥ 0. Beukers further assumes that b = −1, but this scaling can only be made over ℂ and destroys the desired integrality property of the u_{n}, and in fact there seem to be almost no cases of integral solutions with b = −1. (See below.)

To search for rational values of a, b and λ leading to u_{n} ∈ ℤ, we observe that equation (2) for n ∈ {0, 1, 2} gives three linear equations for (a, b, λ) in terms of (u_{1}, u_{2}, u_{3}) which can be solved uniquely if the corresponding determinant is non-zero. We can therefore search over a domain of (u_{1}, u_{2}, u_{3}) ∈ ℤ^{3}, compute (a, b, λ) for each choice of these three initial values, and then see whether the further u_{n} (up to some pre-assigned search limit) are also integral. This search was performed in the range 0 ≤ u_{1} ≤ 30, |u_{2}| ≤ 100, |u_{3}| ≤ 2000 (time on a Sun workstation: 52 hours) and yielded 19 solutions, listed below. If we restrict to the case b = −1 considered by Beukers, then we need only search over (u_{1}, u_{2}) ∈ ℤ^{2}, since these values already determine a, λ, and the further u_{n}. A search in the domain 0 < u_{1} ≤ 100, |u_{2}| ≤ 3000 yielded only a single case (a, λ) = (11, −3) (the Apéry numbers) with integral values of u_{n}.

In all 19 cases found in the range searched, the numbers A = −a/b, B = 1/b and Λ = λ/b were integral. In terms of the new parameters A and B, the recurrence relation (2) becomes

(n + 1)^{2}u_{n+1} − An(n + 1)u_{n} + Bn^{2}u_{n−1} = λu_{n}  (n ≥ 0)  (3)

where we have changed the name of the eigenvalue from Λ back to λ. (The corresponding changes in the differential equation (1) are to take P (t) = Bt^{2} − At + 1 and replace the factor t − λ by Bt − λ.) The observation just made then becomes

Conjecture. If (3) has a solution with u_{0} = 1 and u_{n} ∈ ℤ for all n, then A, B and λ are integral.

Notice that it is equivalent to conjecture that just A is integral, since (3) for n = 0 and n = 1 then gives λ = u_{1} ∈ ℤ, B = (2A + λ)u_{1} − 4u_{2} ∈ ℤ. Also, it seems reasonable to guess that the Conjecture is true even if we do not suppose that u_{0} = 1.

2. Numerical data. Assuming the above conjecture, we can do a new and much more rapid search over integral values of A, u_{1} and u_{2}. As explained above, we can suppose that u_{1} ≥ 0. The search can also be speeded up by noticing that the integrality of u_{3} is equivalent to a congruence for A modulo 1, 3 or 9 (namely: A is arbitrary if u_{1} and u_{2} are both divisible by 3, A ≡ −u_{1}/3 (mod 3) if u_{1} is divisible by 3 but u_{2} is not, and A ≡ 4u_{1}+u_{1}u_{2}(2+3u_{2}+2u_{1}^{2}) (mod 9) if u_{1} is not divisible by 3) and the integrality of u_{4} is equivalent to the simple congruence u_{1} ≡ u_{2} (mod 2). Inserting these two congruences into the triple loop results in a reduction of the search domain by a factor 7/54. We performed the search for the 100 million triples in the domain |A| ≤ 250, 0 ≤ u_{1} ≤ 100, |u_{2}| ≤ 1000 (time on the Sun workstation: 19 hours), finding altogether 36 solutions, listed in the table below. The integrality of the u_{n} was checked up to the limit n = 25 (it was subsequently verified that in each of the 36 “solutions” satisfying this criterion the u_{n} are indeed integral for all n) and this was by no means unnecessarily large, since the maximum value of n_{0} := min{n | u_{n} ∉ ℤ} in this range was in fact 23. This value occurred exactly once, for the triple (a, u_{1}, u_{2}) = (229, 4, −660). The other values of n_{0} which occurred were 5, 7, 8, 9, 11, 13, 16, 17 and 19, the last four occurring for 5926, 45, 344 and 17 triples, respectively. (Note that n_{0} is automatically at least 5 since we have done the search in such a way that u_{0}, …, u_{4} are always integral.) The fact that n_{0} is always a prime power is a consequence of a result of Beukers ([4], Prop. 3.3), according to which the smallest n (if any) for which a given prime p occurs in the denominator of u_{n} is always a power of p.

3. Analysis of the data. Looking at the table, we observe that the solutions found can be divided into five classes. Terminating solutions. The sequences numbered 1, 4, 12, 18, 23, 28, 30, and 32 in the table have u_{n} = 0 for all sufficiently large n (although in the last two cases the table does not extend to the first vanishing value). By inspection, we discover that these values are given by the formula

(A, B, λ) = (−1, 0, d^{2} + d),  u_{n} = \binom{d}{n}\binom{d + n}{n} ,  (4)

where d is a non-negative integer. The corresponding function F (t) = ∑ u_{n}t^{n} is given by

F (t) = F (−d, 1 + d; 1; t) = P_{d}(1 − 2t) .  (5)

Here F (a, b; c; t) denotes the hypergeometric function and P_{d}(X) the dth Legendre polynomial. Polynomial solutions. The solutions numbered 3, 6 and 15 in the table are also easily recognized: they are given by u_{n} = 1, u_{n} = 2n + 1, and u_{n} = 3n^{2} + 3n + 1, respectively. Looking for further

\begin{matrix}index & A & B & λ & u_{0} & u_{1} & u_{2} & u_{3} & u_{4} & u_{5} & u_{6} \\ #1 & −1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ #2 & 0 & −16 & 0 & 1 & 0 & 4 & 0 & 36 & 0 & 400 \\ #3 & 2 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\ #4 & −1 & 0 & 2 & 1 & 2 & 0 & 0 & 0 & 0 & 0 \\ #5 & 7 & −8 & 2 & 1 & 2 & 10 & 56 & 346 & 2252 & 15184 \\ #6 & 2 & 1 & 3 & 1 & 3 & 5 & 7 & 9 & 11 & 13 \\ #7 & 9 & 27 & 3 & 1 & 3 & 9 & 21 & 9 & −297 & −2421 \\ #8 & 10 & 9 & 3 & 1 & 3 & 15 & 93 & 639 & 4653 & 35169 \\ #9 & 11 & −1 & 3 & 1 & 3 & 19 & 147 & 1251 & 11253 & 104959 \\ #10 & 12 & 32 & 4 & 1 & 4 & 20 & 112 & 676 & 4304 & 28496 \\ #11 & 16 & 0 & 4 & 1 & 4 & 36 & 400 & 4900 & 63504 & 853776 \\ #12 & −1 & 0 & 6 & 1 & 6 & 6 & 0 & 0 & 0 & 0 \\ #13 & 17 & 72 & 6 & 1 & 6 & 42 & 312 & 2394 & 18756 & 149136 \\ #14 & 27 & 0 & 6 & 1 & 6 & 90 & 1680 & 34650 & 756756 & 17153136 \\ #15 & 2 & 1 & 7 & 1 & 7 & 19 & 37 & 61 & 91 & 127 \\ #16 & −27 & 0 & 12 & 1 & 12 & −126 & 2100 & −40950 & 864864 & −19171152 \\ #17 & −16 & 0 & 12 & 1 & 12 & −60 & 560 & −6300 & 77616 & −1009008 \\ #18 & −1 & 0 & 12 & 1 & 12 & 30 & 20 & 0 & 0 & 0 \\ #19 & 32 & 256 & 12 & 1 & 12 & 164 & 2352 & 34596 & 516912 & 7806224 \\ #20 & 64 & 0 & 12 & 1 & 12 & 420 & 18480 & 900900 & 46558512 & 2498640144 \\ #21 & 2 & 1 & 13 & 1 & 13 & 55 & 147 & 309 & 561 & 923 \\ #22 & −64 & 0 & 20 & 1 & 20 & −540 & 21840 & −1021020 & 51459408 & −2715913200 \\ #23 & −1 & 0 & 20 & 1 & 20 & 90 & 140 & 70 & 0 & 0 \\ #24 & 2 & 1 & 21 & 1 & 21 & 131 & 471 & 1251 & 2751 & 5321 \\ #25 & 54 & 729 & 21 & 1 & 21 & 495 & 12171 & 305919 & 7794171 & 200412801 \\ #26 & 32 & 256 & 28 & 1 & 28 & 580 & 10992 & 199524 & 3530352 & 61417616 \\ #27 & −27 & 0 & 30 & 1 & 30 & −180 & 2640 & −48510 & 989604 & −21441420 \\ #28 & −1 & 0 & 30 & 1 & 30 & 210 & 560 & 630 & 252 & 0 \\ #29 & 2 & 1 & 31 & 1 & 31 & 271 & 1281 & 4251 & 11253 & 25493 \\ #30 & −1 & 0 & 42 & 1 & 42 & 420 & 1680 & 3150 & 2772 & 924 \\ #31 & 2 & 1 & 43 & 1 & 43 & 505 & 3067 & 12559 & 39733 & 104959 \\ #32 & −1 & 0 & 56 & 1 & 56 & 756 & 4200 & 11550 & 16632 & 12012 \\ #33 & 2 & 1 & 57 & 1 & 57 & 869 & 6637 & 33111 & 124223 & 380731 \\ #34 & −16 & 0 & 60 & 1 & 60 & 420 & −1680 & 13860 & −144144 & 1681680 \\ #35 & −64 & 0 & 84 & 1 & 84 & −924 & 30800 & −1316700 & 62990928 & −3212537328 \\ #36 & −27 & 0 & 84 & 1 & 84 & 630 & −5460 & 81900 & −1493856 & 30126096 \end{matrix}

3 polynomial solutions, we find that the sequences 21, 24, 29, 31 and 33 are also of this type, and that these 8 solutions are the cases 0 ≤ d ≤ 7 of an infinite family given by

(A, B, λ) = (2, 1, d^{2} + d + 1), u_{n} = H_{d}(n) ,

where H_{d}(n) is a polynomial of degree d of the form

H_{d}(n) = \frac{(2d)!}{d!^{3}} n^{d} + ⋯ + (∑_{k=1}^{d} \frac{2}{k} ) n + 1 .

Looking at the first few polynomials H_{d}(x), we find that they satisfy H_{d}(−x − 1) = (−1)^{d}H_{d}(x) and H_{2r}(−\frac{1}{2}) = 2^{−4r}\binom{2r}{r}^{2}, but beyond this they are not easy to recognize. However, the fact that u_{n} is a polynomial of degree d in n means that the corresponding solution F_{d}(t) = ∑ u_{n}t^{n} of the differential equation (1) is a rational function with denominator (1 − t)^{d+1}, and looking at the first few values of the polynomial (1 − t)^{d+1}F_{d}(t) we immediately find that its coefficients are the squares of the binomial coefficients \binom{d}{n}, i.e., these solutions are given by

F (t) = ∑_{n=0}^{∞} H_{d}(n)t^{n} = \frac{F (−d, −d; 1; t)}{(1 − t)^{d+1}} = \frac{1}{1 − t} P_{d}(\frac{1 + t}{1 − t}) ,  (6)

where again P_{d}(X) denotes the dth Legendre polynomial. Notice that for these solutions the discriminant of the polynomial P (t) = Bt^{2} − At + 1 in equation (1) is 0, a degenerate case which Beukers excluded from consideration. Hypergeometric solutions. We next find that the table contains four solutions where u_{n} has a simple multiplicative expression as a multinomial coefficient. These are of two types. The first corresponds to A = λ = 0, so that the recursion (3) relates u_{n+1} and u_{n−1}. Up to normalization, there is only one solution of this kind, namely the sequence #2 in the table,

Sequence #2: (A, B, λ) = (0, −16, 0), u_{n} = \{ \begin{matrix}\binom{2r}{r}2 & if n = 2r, \\ 0 & if n is odd. \end{matrix}

The other type occurs when B = 0, so that the recursion (3) relates only u_{n+1} and u_{n}. There are 10 sequences of this type in our table, with the numbers 11, 14, 16, 17, 20, 22, 27, and 34–36. Some of these are easy to recognize by inspection, e.g.:

Sequence #11: (A, B, λ) = (16, 0, 4),  u_{n} = \binom{2n}{n}^{2} = \frac{(2n)!^{2}}{n!^{4}} ,
Sequence #14: (A, B, λ) = (27, 0, 6),  u_{n} = \binom{3n}{n}\binom{2n}{n} = \frac{(3n)!}{n!^{3}} ,
Sequence #20: (A, B, λ) = (64, 0, 12),  u_{n} = \binom{4n}{2n}\binom{2n}{n} = \frac{(4n)!}{(2n)! n!^{2}} ,

and there are similar formulas in the other cases, e.g. the values of u_{n} for the two sequences #16 and #17 are given by (−1)^{n−1}(3n + 1)!/(3n − 1)n!^{3} and (−1)^{n−1}(2n)!(2n + 1)!/(2n − 1)n!^{4}, respectively. The general case is easily seen to be given by the same formula as in (4), except that d is now an arbitrary rational number rather than a non-negative integer and we must make a corresponding rescaling A → DA, λ → Dλ, u_{n} → D^{n}u_{n} for some positive integer D in order to achieve integrality. (The smallest choice of D is M^{2} ∏_{p|M} p, where M is the denominator of d.) The corresponding generating function is given by the same formula as in (5), except that the variable t must be replaced by Dt and the function P_{d}(X) is now the Legendre function of order d, which is no longer a polynomial when d is not an integer. The index d for the above-listed cases # 11, 14, 16, 17, 20, 22, 27, 34, 35 and 36 are −1/2, −1/3, 1/3, 1/2, −1/4, 1/4, 2/3, 3/2, 3/4 and 4/3, respectvely. (One has a choice here between d and −1 − d as parameters, and we have chosen d ≥ −1/2.) Legendrian solutions. Next, we see that (apart from the polynomial solutions considered above) there are three triples in the table for which A^{2} = 4B, namely, those with labels 19, 25 and 26. Here there is no evident pattern for the numbers u_{n}, but after a bit of thought we realize that these are (up to a normalization t ↦ Dt, u_{n} ↦ D^{n}u_{n} with suitable D ∈ ℕ) given by the same formula as in (6), but where now d is allowed to be a rational number rather than an integer and P_{d}(X) is again the Legendre function of index d, which for non-integral d is not a polynomial in X. (The coefficient H_{d}(n) of t^{n} in F_{d}(t) is also no longer a polynomial in n for fixed d, though it is a polynomial in d for fixed n.) This time the smallest choice of D for given d ∈ ℚ making the coefficients u_{n} = D^{n}H_{d}(n) integral for all n is given by M^{2} ∏_{p|M} p^{∗}, where M is the denominator of d and p^{∗} = 4 if p = 2, p^{∗} = p if p is odd. Sporadic solutions. Finally, the table contains six solutions which do not fall into any of the four infinite families above, and these are the only really interesting ones. (They are also the only ones Beukers considered, since he assumed that the polynomial P had non-vanishing quadratic term, constant term, and discriminant.) For convenience we list these 6 solutions again, with labels A–F. We also add #2 to this list, with label G, since although it is hypergeometric it does not belong to the 1-parameter family of hypergeometric solutions discussed above.

new label index A B λ u_{0} u_{1} u_{2} u_{3} u_{4}  u_{5}  u_{6}
A  #5 7 −8 2 1 2 10 56 346 2252 15184
B  #7 9 27 3 1 3 9 21  9 −297 −2421
C  #8 10 9 3 1 3 15 93 639 4653 35169
D  #9 11 −1 3 1 3 19 147 1251 11253 104959
E  #10 12 32 4 1 4 20 112 676 4304 28496
F  #13 17 72 6 1 6 42 312 2394 18756 149136
G  #2 0 −16 0 1 0 4 0  36  0  400

Since all of these are contained in the domain 0 < A ≤ 17, 0 < u_{1} ≤ 6, 0 < u_{2} ≤ 42, which comprises less than 0,005% of the domain of our search, one could hazard a rather optimistic

Conjecture. Up to normalizations, the only cases where the recursion (3) has an integral solution are the hypergeometric family

HG_{d} :  (A, B, λ) = (−1, 0, d^{2} + d)  (d ∈ ℚ, d ≥ −1/2)

(which for d ∈ ℤ_{≥0} is terminating), the Legendre family

Leg_{d} :  (A, B, λ) = (2, 1, d^{2} + d + 1)  (d ∈ ℚ, d ≥ −1/2)

(which for d ∈ ℤ_{≥0} is polynomial), and the seven sporadic solutions A–G.

4. Binomial coefficient sums. The sequence “D” in the table just given is immediately recognized to be the Apéry numbers, which are given by a well-known expression as a sum of products of binomial coefficients, so it is reasonable to look for similar expressions for the other solutions. Apart from relatively uninteresting expressions of this form for certain hypergeometric solutions such as

\frac{(3n)!}{n!^{3}} = ∑_{k=0}^{2n}(−1)^{n−k}\binom{2n}{k}^{3}, and the evident binomial coefficient sum expression for H_{d}(n) obtained

by equating the coefficients of t^{n} on both sides of (6), we have the following binomial coefficient representations, covering all but one of the sporadic solutions (G is trivial):

Sequence A: u_{n} = ∑_{k=0}^{n} \binom{n}{k}^{3} = ∑_{k=0}^{n} \binom{n}{k}^{2}\binom{2k}{n} ,

Sequence B: u_{n} = ∑_{k=0}^{[n/3]}(−1)^{k} 3^{n−3k}\binom{n}{3k}\binom{3k}{k}\binom{2k}{k} ,

Sequence C: u_{n} = ∑_{k=0}^{n} \binom{n}{k}^{2}\binom{2k}{k} ,

Sequence D: u_{n} = ∑_{k=0}^{n} \binom{n}{k}^{2}\binom{n + k}{k} = ∑_{k=0}^{n}(−1)^{n−k}\binom{n}{k}\binom{n + k}{k}^{2} ,

Sequence E: u_{n} = ∑_{k=0}^{[n/2]} 4^{n−2k} \binom{n}{2k}\binom{2k}{k}^{2} .

The formulas for the sequences B and E say that the corresponding generating functions ∑ u_{n}t^{n} have hypergeometric representations:

F_{B}(t) = \frac{1}{1 − 3t} F (\frac{1}{3}, \frac{2}{3}; 1; (\frac{−3t}{1 − 3t})^{3}) , F_{E}(t) = \frac{1}{1 − 4t} F (\frac{1}{2}, \frac{1}{2}; 1; (\frac{4t}{1 − 4t})^{2}) .

We have listed here only expressions for the u_{n} as simple sums of products of binomial coefficients. One can always expressions as multiple sums, e.g. by the method described at the end of §7, a simple example (courtesy of the referee) being the formula u_{n} = ∑_{0≤j≤i≤n}(−1)^{i}8^{n−i}\binom{n}{i}\binom{i}{j}^{3} in case F.

5. Modular properties. Thanks to Beukers [3], one knows that the Apéry sequence has a modular interpretation, i.e., that for this sequence there is a modular function t(z) such that the function F (t(z)) = ∑ u_{n}t(z)^{n} is a modular form of weight 1. We can therefore ask whether a similar property holds for the other solutions in our list. The case when {u_{n}} terminates or is a polynomial in n is uninteresting since then F (t(z)) is a rational function of t(z) and hence is a modular function for any choice of t(z). However, all the other examples (hypergeometric, Legendrian or sporadic) are interesting, since it is not a priori evident when they can be parametrized by modular forms. It turns out that there is an algorithm to investigate this question. We explain this below. Applying it to each of the 20 non-trivial sequences in the table, we find that all of the sporadic solutions are modular, and that there are also four hypergeometric and two Legendrian sequences which give modular forms. (The specific parametrizations are given below.) This suggests the following

Conjecture. Any integral solution of the differential equation (1), where P (t) is a non-degenerate quadratic polynomial, has a modular parametrization.

Whether or not this is true, the natural question arises whether it is possible to give a complete classification of the modular solutions of the differential equation (1). This probably can be done, since the restricted nature of the singularities of the differential equation should put such strong restrictions on the modular side of the picture (e.g., the corresponding modular curve should have genus 0 and a small number of cusps and elliptic fixed points) as to make a complete listing possible. The most optimistic guess would be that the 12 modular solutions already found are a complete list. (This is supported by the case that they are contained in the domain 0 ≤ A ≤ 54, 0 ≤ u_{1} ≤ 21, 0 < u_{2} ≤ 500, which is less than 0,5% of the domain of our search.) If this and the above conjecture are true, then our list of 6 sporadic solutions is indeed complete. How to recognize modularity. A general fact about modular forms is the following: if f(z) is an arbitrary modular form of positive weight k and t(z) a modular function (i.e., a meromorphic modular form of weight 0), then the power series F (t) obtained by expressing f(z) locally as a power series in t(z) always satisfies a linear differential equation of order k + 1 with algebraic (or, if t(z) is a Hauptmodul, even polynomial) coefficients. A discussion of this phenomenon in the general case, and an algorithm to find the corresponding modular parametrization, if one exists, of a given linear differential equation, is explained in §5 of [11]. Here we are interested only in k = 1, since the differential equation (1) has second order (the case k = 0 is uninteresting, corresponding to the “terminating” and “polynomial” cases above when (1) is not the differential equation of lowest order satisfied by f), and the easiest is simply to write down the differential equation. Suppose that t(z) is a modular function and f(z) a modular form of weight 1. Then the function t′ := (2πi)^{−1} dt/dz is a (meromorphic) modular form of weight 2 and the function 2f′^{2} − ff′′ is a modular form of weight 6, so we can write

\frac{t′(z)}{f(z)^{2}} = α(t(z)) ,  \frac{2f′(z)^{2} − f(z)f′′(z)}{t′(z)f(z)^{4}} = β(t(z))  (7)

where α(t) and β(t) are algebraic (or even, if t(z) is a Hauptmodul, rational) functions of t. From the calculation

\frac{1}{t′(z)} \frac{d}{dz} (\frac{α(t(z))}{t′(z)} \frac{dF (t(z))}{dz}) + β(t(z)) F (t(z)) = \frac{1}{t′} (\frac{1}{f^{2}} f′)′ + \frac{2f′^{2} − ff′′}{t′f^{4}} f = 0

we then see that F (t) is a solution of the linear second order differential equation with algebraic coefficients (αF′)′ + βF = 0. Comparing this with (1), we find that in our case the function α(t), up to a constant factor which we can normalize to be 1, is equal to tP (t) = t(1 − At + Bt^{2}). We can then integrate the equation t′(z) = α(t(z))f(z)^{2} = α(t(z))F (t(z))^{2} to get z = ∫ dt/α(t)F (t)^{2} or, in terms of the standard expansion parameter q = e^{2πiz} of modular functions and forms at infinity,

q = t exp(∫_{0}^{t}(\frac{1}{α(s) F (s)^{2}} − \frac{1}{s}) ds) = t exp(∫_{0}^{t}(\frac{F (s)^{−2}}{1 − As + Bs^{2}} − 1) \frac{ds}{s}) .  (8)

This is a power series with leading term t (recall that F (t) = u_{0} + u_{1}t + ⋯ has leading coefficient 1 by our choice of normalization) and hence can be inverted to compute t = t(z) as a power series in q, after which the function f(z) = F (t(z)) can also be computed. We then only have to look at these power series to see whether they are in fact the q-expansions of a modular function and a weight 1 modular form, respectively. This can often be done by inspection (in many of the cases below, t(z) is a simple product of eta-functions and f(z) an Eisenstein series of weight 1, which are easily recognized from their q-expansions); when this fails, one can test numerically whether the function t(z) is related algebraically on the classical j-invariant j(z). Once one has found a relation by looking at the first few Fourier coefficients of the candidate modular forms, the validity of the asserted differential equation can be checked algorithmically, since any equality between modular forms can be verified by verifying the equality of a specific number of terms of their q-expansions. In this way, the modular parametrizations given below were found. That the other solutions (excluding the trivial ones where the sequence {u_{n}} is either polynomial or terminating and F (t) is a rational function) are not modular is harder to verify algorithmically, but in each case the numerical calculation indicates that their q-expansions have very large denominators involving infinitely many prime numbers, and this is not a possible behavior for modular forms.

Description of all known modular solutions of (1). As already mentioned, we found altogether 12 cases in our table which had modular parametrizations—the six sporadic ones, four hypergeometric ones, and two Legendrian ones. We now discuss each of these individually. Sporadic cases. In four of the six sporadic cases, those with labels A, C, E and F, the quadratic polynomial x^{2} − Ax + B splits over ℤ as (x − α)(x − β). In each of these cases, the modular functions t(z), 1−αt(z) and 1−βt(z) and the modular form f(z) have representations as products of Dedekind eta-functions, and the modular form f(z) can also be represented as an Eisenstein series of weight 1. These representations are given in the following table, in which a notation like 1^{3}2^{−3}3^{−9}6^{9} is a shorthand for the eta-product η(z)^{3}η(2z)^{−3}η(3z)^{−9}η(6z)^{9} and the functions θ_{3} and θ_{4} are the Eisenstein or theta series

θ_{3}(z) = 1 + 6 ∑_{n=1}^{∞}(∑_{d|n}(\frac{−3}{d})) q^{n} = ∑_{m,n∈ℤ} q^{m^{2}+mn+n^{2}} ,

θ_{4}(z) = 1 + 4 ∑_{n=1}^{∞}(∑_{d|n}(\frac{−4}{d})) q^{n} = ∑_{m,n∈ℤ} q^{m^{2}+n^{2}} ,

i.e., the modular forms of weight 1 whose Mellin transforms are the Dedekind zeta-functions of the two imaginary quadratic fields ℚ(\sqrt{−3}) and ℚ(i) possessing non-trivial units.

Case α β  t(z)  1 − αt(z) 1 − βt(z) f(z)  f(z)
A −1 8  \frac{1^{3}6^{9}}{2^{3}3^{9}}  \frac{2^{5}6^{1}}{1^{1}3^{5}}  \frac{1^{8}6^{4}}{2^{4}3^{8}}  \frac{2^{1}3^{6}}{1^{2}6^{3}}  \frac{1}{3}θ_{3}(z) + \frac{2}{3}θ_{3}(2z)
C  1 9  \frac{1^{4}6^{8}}{2^{8}3^{4}}  \frac{1^{1}3^{5}}{2^{5}6^{1}}  \frac{1^{9}6^{3}}{2^{9}3^{3}}  \frac{2^{6}3^{1}}{1^{3}6^{2}}  \frac{1}{2}θ_{3}(z) + \frac{1}{2}θ_{3}(2z)
E  4 8  \frac{1^{4}4^{2}8^{4}}{2^{10}}  \frac{1^{4}4^{14}}{2^{14}8^{4}}  \frac{1^{8}4^{4}}{2^{12}}  \frac{2^{10}}{1^{4}4^{4}}  θ_{4}(z)
F  8 9 \frac{1^{5}3^{1}4^{5}6^{2}12^{1}}{2^{14}}  \frac{1^{8}4^{8}6^{20}}{2^{20}3^{8}12^{8}}  \frac{1^{9}4^{9}6^{6}}{2^{18}3^{3}12^{3}}  \frac{2^{15}3^{2}12^{2}}{1^{6}4^{6}6^{5}} θ_{3}(z) + 2θ_{3}(2z) − 2θ_{3}(4z)

Case B is fairly similar: here the polynomial x^{2} − Ax + B = x^{2} − 9x + 27 does not split over ℤ, but we still have the eta-product representations

Case B:  t = \frac{1^{3}4^{3}18^{9}}{2^{9}9^{3}36^{3}} ,  1 − 9t + 27t^{2} = \frac{1^{9}4^{9}6^{36}9^{3}36^{3}}{2^{27}3^{12}12^{12}18^{9}} ,  f = \frac{2^{9}3^{1}12^{1}}{1^{3}4^{3}6^{3}}

and the Eisenstein series representation

f(z) = \frac{1}{2} θ_{3}(z) − \frac{3}{2} θ_{3}(3z) − θ_{3}(4z) + 3 θ_{3}(12z) .

Finally, in the last case D (corresponding to the Apéry numbers), the function t still has a product representation, but no longer as a product of eta-functions, and similarly for 1 − At + Bt^{2} and f:

Case D:  t = q ∏_{n=1}^{∞}(1 − q^{n})^{5( \frac{n}{5} )} ,  1 − 11t − t^{2} = t · \frac{1^{6}}{5^{6}} ,  f^{2} = t^{−1} · \frac{5^{5}}{1^{1}} .

The Eisenstein series representation of f in this case is also more complicated:

f(z) = 1 + ∑_{n=1}^{∞} ∑_{d|n}(\frac{3 − i}{2} χ(d) + \frac{3 + i}{2} χ̄(d)) q^{n} ,

where χ is the Dirichlet character of conductor 5 with χ(2) = i, χ^{2} = (\frac{·}{5}). Hypergeometric cases. Here two cases, the exceptional case #2 when A = λ = 0 and the sequence #11 with (A, B, λ) = (16, 0, 4), behave like the first four sporadic cases:

\begin{matrix}Case & α & β & t(z) & 1 − αt(z) & 1 − βt(z) & f(z) & f(z) \\ #2 & −4 & 4 & \frac{2^{4}8^{8}}{4^{12}} & \frac{2^{14}8^{4}}{1^{4}4^{14}} & \frac{1^{4}2^{2}8^{4}}{4^{10}} & \frac{4^{10}}{2^{4}8^{4}} & θ_{4}(2z) \\ #11 & 0 & 16 & \frac{1^{8}4^{16}}{2^{24}} & — & \frac{1^{16}4^{8}}{2^{24}} & \frac{2^{10}}{1^{4}4^{4}} & θ_{4}(z) \end{matrix} (Note that, amusingly enough, f(z) for case #11 is the same as for the sporadic case E = #10,

although of course the parameter t(z) and hence also the power series F (t) are different.) Another case, corresponding to sequence #14 with (A, B, λ) = (27, 0, 6), is more similar to the sporadic case D, since we do not have product expansions for t, 1 − At and f separately, but only for certain multiplicative combinations, while f has an ordinary (indeed, a particularly simple) Eisenstein series representation:

Case #14:  \frac{t}{1 − 27t} = \frac{3^{12}}{1^{12}} ,  f^{3} · t = \frac{3^{9}}{1^{3}} ,  f(z) = θ_{3}(z) .

Finally, for the sequence #20 with (A, B, λ) = (64, 0, 12), we again find eta-product representations only for certain multiplicative combinations of t, 1 − At and f, but this time f is not a modular form at all, but only the square-root of a modular form:

Case #20:  \frac{t}{1 − 64t} = \frac{2^{24}}{1^{24}} ,  f^{4} · t = \frac{2^{16}}{1^{8}} ,  f(z)^{2} = −E_{2}(z) + 2 E_{2}(2z) ,

where E_{2}(z) = 1 − 24 ∑ σ_{1}(n) q^{n} is the usual not-quite-modular Eisenstein series of weight 2 on the full modular group. Legendrian cases. The last two cases which turned out to be modular are again slightly different, since now f is no longer a holomorphic modular form, or even the square-root of one, but rather is a meromorphic modular form of weight one:

Case #19: t = \frac{1^{8}4^{16}}{2^{24}} ,  1 − 16t = \frac{1^{16}4^{8}}{2^{24}} ,  f = \frac{2^{22}}{1^{12}4^{8}} ,  f(z) = \frac{−θ_{4}(z) + 2θ_{4}(2z)}{1 − 16t(z)} .

Case #25:  t = \frac{1^{12}4^{12}6^{36}}{2^{36}3^{12}12^{12}} ,  (1 − 27t)^{2}f^{3} = \frac{2^{27}3^{3}12^{3}}{1^{9}4^{9}6^{9}} ,  f(z) = \frac{−θ_{3}(z) + 2θ_{3}(4z)}{1 − 27t(z)} .

Discussion. We make a few remarks about the modular parametrizations just given. We consider only the non-degenerate case when B ≠ 0 and A^{2} ≠ 4B, so that x^{2} − Ax + B = (x − α)(x − β) with α and β distinct and non-zero. We further assume that t is a Hauptmodul, i.e., that t and f are modular forms on some group Γ ⊂ SL(2, ℝ) of genus 0 and that z ↦ t(z) gives an isomorphism of ℌ/Γ ∪ {cusps} to ℙ^{1}(ℂ). From equation (7) we see that the modular forms f(z), t(z) and t′(z) are related by f^{2} = t′/t(1 − αt)(1 − βt), so that the Γ-invariant differential form f(z)^{2}dz is (or pushes forward to) the differential form dt/t(1 − αt)(1 − βt) on ℙ^{1}. This form has three simple poles, at t = 0, t = 1/α and t = 1/β, and one simple zero at t = ∞. In terms of the “modular” variable z, the point t = 0 corresponds to the cusp at infinity (since we have constructed t in such a way that it has an expansion q + ⋯). None of the functions t(z), 1 − αt(z) or 1 − βt(z) had a pole or zero in the upper half-plane, since then the function f^{2} would have a (simple) pole at that point, which is impossible. This explains why these three modular functions always had product expansions in the cases looked at above. (This actually remains formally true even in the two sporadic cases B and D when α and β are not in ℤ, but must be interpreted as the statement that the logarithmic derivatives of 1 − αt and 1 − βt are weight 2 Eisenstein series.) Also, the fact that the genus is zero means that there are no cusp forms of weight 1 or 2, so that f (and, for that matter, also f^{2}) must be an Eisenstein series, as indeed we found in all the examples. It seems plausible that by continuing this analysis one could get enough information about the necessary forms of f and t to classify at least all those cases when t is a Hauptmodul and f a true modular form of weight 1, but I have not tried to do this.

6. The periods associated to the six sporadic solutions. The most famous property of the original Apéry recursion (n + 1)^{2}u_{n+1} − (11n^{2} + 11n + 3)u_{n} − n^{2}u_{n−1} = 0 is, of course, that if we start with the two linearly independent solutions \{u_{n}\} = \{1, 3, 19, … \} and \{v_{n}\} = \{0, 1, \frac{25}{4}, … \}, then the ratio v_{n}/u_{n} tends to the limit L = \frac{1}{5}ζ(2) and this convergence is rapid enough to prove the irrationality of L. In this section we look at the corresponding questions for each of the six sporadic sequences A–F. We will see that the limit exists and can be computed in five of the six cases, but that the convergence, though always exponential, is only rapid enough to imply the irrationality of the limit in the case of the Apéry equation.

We denote by {u_{n}}, as up to now, the solution of equation (3) with u_{0} = 1 (and hence u_{1} = λ, u_{n} ∈ ℤ for all n by the fundamental property of our sequences), and by {v_{n}} the sequence defined by the initial values v_{0} = 0, v_{1} = 1 and the requirement that (3) holds for all n ≥ 1. We want to determine whether the sequence v_{n}/u_{n} tends to a limit L and, if so, to evaluate this limit and specify the speed of the convergence.

Let α and β as in §5 denote the roots of x^{2} − Ax + B = 0, chosen (except in case B where α, β = (9 ± 3i\sqrt{3})/2) so that |β| > |α|. Then the sequence \{u_{n}\} has the asymptotic behavior

u_{n} = C \frac{β^{n}}{n} (1 − \frac{c_{1}}{n} + O(\frac{1}{n^{2}} ))  (n → ∞) ,  (9)

where c_{1} = \frac{λ − α}{β − α}, while every other element of the 2-dimensional vector space of real-valued solutions

of the recursion (3) (for n ≥ n_{0} ≥ 1) has asymptotics given either by the same formula (with a possibly

different value of C) or else by (const.)\frac{α^{n}}{n} (1 − \frac{c_{2}}{n} + O(\frac{1}{n^{2}} )), where c_{2} = \frac{λ − β}{α − β} , with the space of

solutions of the latter type being 1-dimensional. (To obtain these asymptotic formulas, consider a general solution {u_{n}} of (3) and make the Ansatz u_{n} = (const.) n^{ε}x^{n}(1 − cn^{−1} + O(n^{−2})); then the ratio ρ_{n} := u_{n+1}/u_{n} equals x(1 + εn^{−1} + (\binom{ε}{2} + c)n^{−2} + …), and rewriting equation (3) in the form (1+n^{−1})^{2}ρ_{n}−A(1+n^{−1})−λn^{−2}+Bρ_{n−1}^{−1} = 0 we find x ∈ {α, β}, ε = −1 and c = (λ−A+x)/(2x−A).) It follows that there are real constants C′ and L such that v_{n} − Lu_{n} = C′ \frac{α^{n}}{n} (1 − \frac{c_{2}}{n} + O(\frac{1}{n^{2}} )) and that the quotient v_{n}/u_{n} tends to L exponentially quickly:

\frac{v_{n}}{u_{n}} = L + C′′ (\frac{α}{β})^{n} (1 + \frac{c_{1} − c_{2}}{n} + O(\frac{1}{n^{2}} ))  (n → ∞) .  (10)

10 The values of the constants β, α, C, c_{1}, c_{2} and L occurring in equations (9) and (10) are given for each of the five cases A, C, D, E and F in the table below, in which ϕ in the third row denotes the golden ratio (1 + \sqrt{5})/2 and the numbers ζ(2), L_{−3}(2) and L_{−4}(2) are the values at s = 2 of the Riemann zeta function and of the Dirichlet L-series L(s, χ) associated to the quadratic characters χ(n) = (\frac{−3}{n}) and χ(n) = (\frac{−4}{n}), each of which is a famous number: ζ(2) equals π^{2}/6, L_{−4}(2) is “Catalan’s constant,” which occurs in the evaluation of many classical definite integrals, and L_{−3}(2) is, up to a factor 3^{3/2}/4, the maximum volume of a tetrahedron in hyperbolic 3-space. The constant C in (9) can be evaluated using either the explicit binomial coefficient expressions given (except for case F) in §4 or else from the modular form expressions of §5. The evaluation of the limiting ratio L will be discussed at the end of this section.

Case β  α  C  c_{1} c_{2}  L
A  8 −1  \frac{2}{π\sqrt{3}}  \frac{1}{3}  \frac{2}{3}  \frac{1}{4}ζ(2) = 0.4112335 ⋯
C  9  1  \frac{3\sqrt{3}}{4π}  \frac{1}{4}  \frac{3}{4}  \frac{1}{2}L_{−3}(2) = 0.3906512 ⋯
D ϕ^{5} −ϕ^{−5}  \frac{ϕ^{5/2}}{2π \sqrt[4]{5}}  \frac{1}{ϕ\sqrt{5}} \frac{ϕ}{\sqrt{5}}  \frac{1}{5}ζ(2) = 0.3289868 ⋯
E  8  4  \frac{2}{π}  0  1 \frac{1}{2}L_{−4}(2) = 0.4579827 ⋯
F  9  8  \frac{3\sqrt{3}}{π} −2 3 \frac{5}{8}L_{−3}(2) = 0.4883140 ⋯

In each of these five cases, formula (10) gives an interesting series of rational approximations for one of the numbers ζ(2), L_{−3}(2) or L_{−4}(2), but, as already mentioned, except in case D these do not converge quickly enough to yield the irrationality of the limit. Indeed, from the table we see that in cases E and F the quantity v_{n} − Lu_{n} blows up exponentially like 4^{n}/n or 8^{n}/n, respectively, and even in cases A and C, where |α| = 1 and hence v_{n} − Lu_{n} tends to zero like O(1/n), this is not enough to give the irrationality of L because v_{n} itself has a denominator which blows up like e^{2n}. The speed of convergence is best in these two cases, with v_{n}/u_{n} − L being of the order of (1/8)^{n} and (1/9)^{n}, respectively, while the convergence in cases E (Catalan’s constant) and F is only like (1/2)^{n} and (8/9)^{n}, respectively. The approximations lead in each of the five cases to a simple infinite continued fraction expansion of the number L, the formulas in cases C and E being

L_{−3}(2) = \frac{2}{P (0) − \frac{9 · 1^{4}}{P (1) − \frac{9 · 2^{4}}{P (2) − \frac{9 · 3^{4}}{...}}}} , L_{−4}(2) = \frac{1/2}{Q(0) − \frac{2 · 1^{4}}{Q(1) − \frac{2 · 2^{4}}{Q(2) − \frac{2 · 3^{4}}{...}}}}

with P (n) = 10n^{2} + 10n + 3, Q(n) = 3n^{2} + 3n + 1. Continued fraction expansions of a similar type, though with more complicated rational functions, have been given recently by Zudilin [12].

To find the value of L listed in the above table, we follow the method devised by Beukers [3] in 1987 to obtain modular interpretations of Apéry’s irrationality proofs for ζ(2) and ζ(3), using the modular descriptions of the differential equations associated to the recursion (3) which were given in §5. Consider for instance case C. We saw that the integer solution {u_{n}} of (3) in this case is given by the generating function expansion f(z) = ∑ u_{n}t(z)^{n}, where t(z) = η(z)^{4}η(6z)^{8}/η(2z)^{8}η(3z)^{4} is the “Hauptmodul” for Γ_{0}(6) and f(z) is the theta series (or Eisenstein series) of weight 1 for Γ_{0}(6) given in §5. The generating function h(z) = ∑ v_{n}t(z)^{n} for the second solution of the recurrence (3) has a more complicated modular interpretation as f(z)g̃(z), where g̃(z) = q − \frac{5}{4}q^{2} + q^{3} − \frac{11}{16}q^{4} + ⋯ (q = e^{2πiz}) is the “Eichler integral” ∑ b(n)q^{n}/n^{2} associated to the weight 3 Eisenstein series g(z) = ∑ b(n)q^{n} = q − 5q^{2} + 9q^{3} − 11q^{4} + ⋯ with b(n) = ∑_{d|n}(−1)^{d−1}(\frac{−3}{n/d})d^{2}. (To prove this, one simply verifies that the differential equation satisfied by fg̃ as a function of t(z) coincides with the one corresponding to the recursion defining the v_{n}.) The group Γ_{0}(6) has four cusps, at z = 0, \frac{1}{2}, \frac{1}{3} and ∞, the corresponding values of the modular function t(z) being \frac{1}{9}, ∞, 1 and 0, respectively. In particular, the asymptotics of the u_{n} and the v_{n} are determined by the singularities of the functions f(z) and h(z) as z → 0 (this explains once again why both of these sequences grow roughly like 9^{n}), and hence the limiting value L of v_{n}/u_{n} is the limiting value as z → 0 of the ratio h(z)/f(z) = g̃(z). (Note that, since ∑(v_{n}−Lu_{n})t^{n} has a larger radius of convergence than ∑ v_{n}t^{n} or ∑ u_{n}t^{n} separately, the value of L is automatically equal to the limit of the ratio of the two latter series as t approaches

the smallest singularity.) Writing the formula defining g̃(z) as ∑_{m≥1} \frac{(\frac{−3}{m})}{m^{2}} \frac{q^{m}}{1 + q^{m}}, one sees that its

limit for z → 0 or q → 1 equals \frac{1}{2}L_{−3}(2), as claimed. The proof in all the other cases is similar,

the function g̃(z) in case E, for instance, being the Eichler integral ∑_{m≥1} \frac{(\frac{−4}{m})}{m^{2}} \frac{q^{m}}{1 + q^{m}} of a weight 3

Eisenstein series on Γ_{0}(8).

Finally, we say a few words about case B, which has been omitted up to now. Here the two roots α = (9 + 3\sqrt{−3})/2 and β = ᾱ of x^{2} − Ax + B = 0 have the same absolute value, so all solutions of (3) have the same growth O(27^{n/2}/n), and the numbers v_{n}/u_{n} do not tend to a limit as n → ∞. Instead, the asymptotic expansions of {u_{n}} and {v_{n}} as n → ∞ are given by

u_{n} = \frac{U}{n} α^{n} (1 − \frac{γ_{1}}{n} + ⋯) + \frac{Ū}{n} ᾱ^{n} (1 − \frac{γ̄_{1}}{n} + ⋯) ,
v_{n} = \frac{V}{n} α^{n} (1 − \frac{γ_{1}}{n} + ⋯) + \frac{V̄}{n} ᾱ^{n} (1 − \frac{γ̄_{1}}{n} + ⋯) ,

where U = \frac{9\sqrt{3}}{2πα}, γ_{1} = \frac{α}{9} and \frac{V}{U} = \frac{1}{2}L_{−3}(2)− \frac{4i}{9\sqrt{3}}ζ(2). (The ratio V/U, which is the analogue of the

limit L in the previous examples, is obtained just as before: the numbers 1/α and 1/ᾱ are the values of the Hauptmodul t(z) given in §5 at the cusps z = ±1/6 and the function g̃(z) = h(z)/f(z) has the form g̃(z) = − ∑_{n=1}^{∞}(\frac{−3}{n})n^{−2}(−q)^{n}/(1 − (−q)^{n}), which for z = 1/6, −q = e^{−2πi/3} has the value given.) Since α/\sqrt{27} = e^{πi/6}, the quotients v_{n}/u_{n} do tend to a limit \frac{1}{2}L_{−3}(2) + \frac{4}{9\sqrt{3}}ζ(2) tan \frac{π(j−1)}{6} as n tends to infinity in any fixed residue class n ≡ j ≢ 4 (mod 6), but the convergence is very slow, with error term of the order of 1/n instead of exponentially small as before, so that from the point of view of rational approximation this case is of less interest than the other five.

7. The modular cases as Picard-Fuchs differential equations. The fact that the nondegenerate solutions of our problem have parametrizations by modular forms means that they have an algebraic-geometric interpretation as the Picard-Fuchs differential equations of certain families of elliptic curves. In fact, we can be much more precise. In 1982, Beauville [2] classified completely the stable families of elliptic curves over ℙ^{1} having exactly four singular fibres (a stable family of elliptic curves over ℙ^{1} is a family in which the singular fibres have only double points; the minimal number of singular fibres is then 4) and found that there were exactly 6. For each family he gave a defining equation, a description of the singular fibres, and the subgroup of SL_{2}(ℤ) describing the family. We reproduce the table from [2] here.

Case  Equation of the family  Group
I  X^{3} + Y^{3} + Z^{3} + tXY Z = 0  Γ(3)
II X(X^{2} + Z^{2} + 2ZY) + tZ(X^{2} − Y^{2}) = 0 Γ_{1}(4) ∩ Γ(2)
III X(X − Z)(Y − Z) + tZY (X − Y) = 0  Γ_{1}(5)
IV (X + Y)(Y + Z)(Z + X) + tXY Z = 0  Γ_{1}(6)
V  (X + Y)(XY − Z^{2}) + tXY Z = 0  Γ_{0}(8) ∩ Γ_{1}(4)
VI  X^{2}Y + Y^{2}Z + Z^{2}X + tXY Z = 0  Γ_{0}(9) ∩ Γ_{1}(3)

Comparing the number six of cases found by Beauville with the number six of sporadic cases A–F in §5 makes it tempting to suppose that there is a 1:1 correspondence between these, with the Picard-Fuchs equations of Beauville’s six families of elliptic curves being precisely the equations (1) corresponding to A–F in some order. But this is not the case, for three reasons. First of all, two of Beauville’s families (V and VI) are isogenous to two of the others (II and I, respectively) and give the same Picard-Fuchs equations. This reduces the first “six” above to “four.” Secondly, our hypergeometric (and hence by definition non-sporadic) solutions #2, #11 and #14, whose modularity was noted in §5, are also Picard-Fuchs equations of families on Beauville’s list. This increases the second “six” above to “nine.” And thirdly, each of Beauville’s surfaces can give rise to several different Picard-Fuchs equations, depending on how one chooses the Hauptmodul t(z). To get a differential equation of the form (1), we must require that the singular fibres correspond to the four values 0, t_{1}, t_{2} and ∞, where t_{1} and t_{2} are the roots of P (t), since these are the singular points of the differential equation (1). So we must put one fibre at t = 0 and one at ∞, giving 12 choices (the remaining freedom of interchanging or rescaling t_{1} and t_{2} has no effect on the parameters A, B, λ of this paper). Of course, these do not necessarily all give rise to equations of the form (1) with rational coefficients, nor (since the elliptic surfaces have non-trivial automorphisms) need they all be different. The actual correspondence turns out to be as follows:

Beauville Family:  IV  I, VI  III  II, V Picard-Fuchs equation:  A, C, F  B, #14  D  E, G, #11

One way to see this is to compute the j-invariant of each of Beauville’s family as a rational function of t, invert this equation to compute t as an explicit modular function, and compare with the list of t’s which we found before. The j-invariant in turn is most conveniently found by rewriting Beauville’s equations in Weierstrass form y^{2} = f(x) (f a polynomial of degree 3) . For the reader’s convenience, we give the transformation needed in each of the 6 cases do this, and the resulting formula for the j-invariant:

Transformation to Weierstrass form  j-invariant
I  X, Y = t(x ± y), Z = 2(1 + 3x)  −t^{3}(t^{3} − 216)^{3}/(t^{3} + 27)^{3}
II  X = xt, Y = x + y, Z = t  2^{8}(t^{4} − t^{2} + 1)^{3}/t^{4}(t^{2} − 1)^{2}
III X = 1 + (y + tx)/(1 − x), Y = 2x, Z = 2 (t^{4} − 12t^{3} + 14t^{2} + 12t + 1)^{3}/t^{5}(t^{2} − 11t − 1)
IV  X, Y = 1 + tx ± y, Z = 2x(1 + tx)  (t + 2)^{3}(t^{3} + 6t^{2} − 12t + 8)^{3}/t^{3}(t − 1)^{2}(t + 8)
V  X, Y = 1 + tx ± y, Z = 2x(1 + tx)  (t^{4} + 16t^{2} + 16)^{3}/t^{2}(t^{2} + 16)
VI  X = 2x, Y = y − xt − 1, Z = 2x^{2}  −t^{3}(t^{3} + 24)^{3}/(t^{3} + 27)

The necessary calculations are explained and carried out in detail in the recent paper [9] by Helena Verrill, in which the reader can also find a lot of other information about these elliptic surfaces and their corresponding Picard-Fuchs equations. A further discussion of some of the above differential equations and their relation to Beauville’s classification and to the earlier classification by Schmickler-Hirzebruch [6] of elliptic pencils with exactly three singular fibers is given in [8], to which the reader is referred for further details. I mention briefly a few further aspects which were pointed out to me by Jan Stienstra. First of all, when different differential equations are associated to the same family of elliptic curves, then the corresponding modular functions have to be the same (possibly after a Möbius transformation z ↦ M(z), M ∈ GL_{2}(ℤ), of the z-variable) up to a fractional linear transformation. For instance, from the above table one sees that A and C correspond to the same fibration, and from the eta-product expansions given in §5 one finds that t_{C}(z) = t_{A}(z)/(1 + t_{A}(z)). Secondly, just as one has isogenies between different families of elliptic curves in Beauville’s list, in certain cases one has further transformations induced by pull-backs under t ↦ t^{N}. Since these transformations change the nature and number of the singularities, one can get examples of Apérylike differential equations coming from the pencils with three exceptional fibers classified in [6] as well as from the ones with four exceptional fibers classified in [2]. This is discussed in detail in [8]. Finally, and most interesting, in the modular cases one can always get explicit formulas for the coefficients u_{n} as sums of multinomial coefficients (although these will in general be multiple sums rather than simple sums like the ones listed in §4) from the equation of the corresponding pencil: if the equation of this pencil is written as t = F (X, Y, Z), where F (X, Y, Z) is a Laurent polynomial in X, Y and Z of degree 0 (as is the case for all of the families in Beauville’s list as reproduced above, after making a linear change of variables in cases II and III), then the coefficient u_{n} of t^{n} in the corresponding power series is equal simply to the constant term of F (X, Y, Z)^{n}. For instance, the constant term in the nth power of −(X^{3} + Y^{3} + Z^{3})/XY Z, corresponding to the “t” in Beauville’s family I, is equal (up to sign) to n!/(n/3)!^{3} if 3|n and to 0 otherwise, and we recover the hypergeometric solution #14, and similarly the constant terms of the powers of (X + Y)(Y + Z)(Z + X)/XY Z, corresponding to the family IV, are the numbers 1, 2, 10, 56, … of sequence #5 in our list.

8. Final remarks. The numerical experiments described in the first sections of this paper were carried out in 1997–98, inspired by the beautiful lecture given by Beukers on his work [4] at the conference in honor of Schinzel in Zakopane in 1997. A first version of the paper was written in 2000 and has been circulating as an informal preprint ever since, with few changes except for the addition of §6 on the periods associated to Apéry-like recursions (which was written in answer to a question posed by A. Connes during my course at the Collège de France in 2001) and some additions to the discussion of the connection with Beauville’s classification in §7. The main reason for not publishing it earlier was that I thought that more explanations of the connections to geometry should be added, but did not understand these well enough. In the meantime, several other related papers, both on the geometric and on the purely differential equations/modular forms side have appeared (in some cases quoting the informal preprint version of this paper), and part of what is presented here is perhaps obsolete. I have nevertheless kept the entire text since this paper in any case contains no real theorems but is to be seen more as an informal discussion of various experimental and theoretical aspects of the three-way connection between algebraic geometry, linear differential equations, and the theory of modular forms. For further material on this subject, the reader is referred in particular to the three papers [1], [8] and [9], all of whose authors—Gert Almkvist, Wadim Zudilin, Jan Stienstra and Helena Verrill, I would also like to thank here for many useful discussions and for their patience.

I should perhaps also mention explicitly that (of course) not all of the examples occurring in this paper are new. In particular, the equations D, A and C (= #9, #5 and #8) are discussed in the paper [7] by Stienstra and Beukers, where two of them are ascribed to Apéry and Cusick but in at least one case had already been found in the 19th century (J. Franel), as I have been informed by Zudilin. I have made no attempt at all to find the earliest reference for each differential equation. In a different direction, it should also be mentioned that the list of modular cases given in §5 is not complete. No further sporadic Apéry-like equations have been discovered, and the conjecture formulated in §3 still stands, but the infinite families of equations which we called “hypergeometric” and “Legendrian” both contain further modular cases which I overlooked because their parameters were beyond the bounds of my numerical search. Specifically, the hypergeometric case (A, B, λ) = (−1, 0, d(d + 1)) (d ∈ ℚ, d ≥ −1/2) gives a modular function for the parameter value d = −1/6 as well as for the three values d = −1/2, −1/3 and −1/4 (= #11, #14 and #20, respectively) which we discussed, the corresponding coefficients u_{n} being up to a scaling factor the multinomial coefficients (6n)!/n!(2n)!(3n)!, and the Legendre case (A, B, λ) = (2, 1, d(d + 1) + 1) (d ∈ ℚ, d ≥ −1/2) is modular not only for the two cases d = −1/2 and d = −1/3 (= #19 and #25) which we gave, but also for d = −1/4 and d = −1/6, found by van Enckevort and van Straten. All of these examples, of course, are contained in the much more extensive tables found in the paper [1], which also gives relevant references and a discussion of how equations of the type discussed here can be used to construct (fourth order) Calabi-Yau differential euquations by means of Hadamard products. Finally, I mention two recent preprints on closely related subjects. In [10], Y. Yang gives an analysis of the periods associated to Apéry-like differential equations similar to the one given in §6 here. And in [5], Vasily Golyshev describes a class of interesting higher order generalizations of the differential equations (1), which he calls equations of type DN. Here N is the order of the differential equation, so that our equation (1) is an example of his D2. For D2 equations he finds a connection with the quantum cohomology of Del Pezzo surfaces, and for D3 equations he gives both a connection with the quantum cohomology of certain Fano threefolds and a relation to modular forms similar to our Section 5.

References

[1] G. Almkvist and W. Zudilin, Differential equations, mirror maps and zeta values. In Mirror symmetry. V, AMS/IP Stud. Adv. Math. 38, Amer. Math. Soc. Providence, RI (2006), 481–515.

[2] A. Beauville, Les familles stables de courbes elliptiques sur P^{1} admettant quatre fibres singulières. C. R. Acad. Sc. Paris 294 (1982), 657–660.

[3] F. Beukers, Irrationality proofs using modular forms. Journées Arithmétiques de Besançon. Astérisque 147–148 (1987), 271–283.

[4] F. Beukers, On B. Dwork’s accessory parameter problem. Math. Z. 241 (2002), 425–444..

[5] V. Golyshev, Classification problems and mirror duality. Preprint (2005), arXiv:math/0510287.

[6] U. Schmickler-Hirzebruch, Elliptische Flächen über ℙ_{1}ℂ mit drei Ausnahmefasern und die hypergeometrische Differentialgleichung. Schriftenreihe des Mathematischen Instituts der Universität Münster 33 (1985)..

[7] J. Stienstra and F. Beukers, On the Picard-Fuchs equation and the formal Brauer group of certain elliptic K3surfaces. Math. Ann. 271 (1985), 269–304.

[8] J. Stienstra, Mahler measure variations, Eisenstein series and instanton expansions. In Mirror symmetry. V, AMS/IP Stud. Adv. Math. 38, Amer. Math. Soc. Providence, RI (2006), 139–150..

[9] H. Verrill, Picard-Fuchs equations of some families of elliptic curves. Proceedings on Moonshine and Related Topics (Montreal, 1999), CRM Proc. Lecture Notes 30, AMS (2001), 253–268.

[10] Y. Yang, Apéry limits and special values of L-functions. Preprint (2007), arXiv:0709:1968v1.

[11] D. Zagier, Elliptic modular forms and their applications. In The 1–2–3 of Modular Forms : Lectures at a Summer School in Nordfjordeid, Norway, K. Ranestad (ed.), Universitext, Springer-Verlag, Berlin (2007), 1–103.

[12] W. Zudilin, Apéry-like difference equation for Catalan’s constant. Electron. J. Combin. 10 (2003), #R14, 10 pp.