
\documentstyle{amsppt} \magnification=\magstephalf \NoRunningHeads 
 \hsize=16.8 true cm \vsize=22.3 true cm  \NoBlackBoxes \TagsOnRight \nologo
\input ruled \tableinfofalse \thicksize=0pt \TightTables
\def\tstrut{\vrule height 2.7ex depth 1.2ex width 0pt}


\def\={\,=\,} \def\f{\frac} \def\bn{\binom} \def\cd{\cdot} \def\Z{\Bbb Z} \def\Q{\Bbb Q} \def\N{\Bbb N}
\def\a{\alpha} \def\b{\beta} \def\l{\lambda} \def\th{\theta} \def\v{\varepsilon} \def\g{\gamma} 
\def\z{\zeta}  \def\G{\Gamma} \def\tg{\widetilde g}

\topmatter\title Integral solutions of Ap\'ery-like recurrence equations\endtitle
\author Don Zagier\endauthor \abstract In [4], Beukers studies the differential
equation $$ \bigl((t^3+at^2+bt)F'(t)\bigr)'+(t-\l)F(t)\=0\,,\tag *$$ where $a$, $b$ and $\l$ are
rational parameters, and asks for which values of these parameters this equation has a solution
in $\Z[[t]]$, the motivating example being the Ap\'ery sequence with $a=11$, $b=-1$, $\l=-3$.
We describe a search over a suitably chosen domain of 100 million triples $(a,b,\l)$. 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 $\Bbb P^1$ with four singular fibres. \endabstract \endtopmatter


\bigskip
{\bf 1. Beukers's recurrence equation.} Beukers [4] considers the differential equation
$$ \bigl(tP(t)F'(t)\bigr)'+(t-\l)F(t)\=0\tag 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)=\sum\limits_{n=0}^\infty u_nt^n$ with $u_0=1$ and
$$ b(n+1)^2u_{n+1}+(an^2+an-\l)u_n+n^2u_{n-1}=0\qquad(n\ge0)\tag2$$
so that $u_1=\l/b$, $u_2=(\l^2-2a\l+b)/4b^2$, etc.  We are interested in finding
values of $(a,b,\l)$ for which $F(Dt)\in\Z[[t]]$ or $D^nu_n\in\Z$ ($\forall n\ge0$)
for some $D\in\N$. After a rescaling $(a,b,\l)\mapsto(a/D,b/D^2,\l/D)$, 
$u_n\mapsto D^nu_n$ we can assume $D=1$ or $u_n\in\Z$,
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\ge0$.  Beukers further assumes that $b=-1$, but this 
scaling can only be made over $\Bbb C$ 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 $\l$ leading to $u_n\in\Z$,
we observe that equation (2) for $n\in\{0,1,2\}$ gives three linear 
equations for $(a,b,\l)$ 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)\in\Z^3$, compute $(a,b,\l)$ 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\le u_1\le30$, $|u_2|\le 100$, $|u_3|\le2000$
(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)\in\Z^2$, since these values already determine $a$, $\l$, and
the further $u_n$.  A search in the domain $0<u_1\le100$, $|u_2|\le 3000$ yielded 
only a single case $(a,\l)=(11,-3)$ (the Ap\'ery 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 $\Lambda=\l/b$
were integral.  In terms of the new parameters $A$ and $B$, the recurrence relation (2) becomes
$$ (n+1)^2u_{n+1}-An(n+1)u_n+Bn^2u_{n-1}=\l u_n\qquad(n\ge0)\tag3$$
where we have changed the name of the eigenvalue from $\Lambda$ back to $\l$. (The
corresponding changes in the differential equation~(1) are to take $P(t)=Bt^2-At+1$
and replace the factor $t-\l$ by $Bt-\l$.) The observation just made then becomes
\proclaim{Conjecture} If $(3)$ has a solution with $u_0=1$ and $u_n\in\Z$ for all $n$, 
then $A$, $B$ and $\l$ are integral.\endproclaim
\noindent Notice that it is equivalent to conjecture that just $A$ is integral, since (3) 
for $n=0$ and $n=1$ then gives $\l=u_1\in\Z$, $B=(2A+\l)u_1-4u_2\in\Z$. Also, it seems 
reasonable to guess that the Conjecture is true even if we do not suppose that $u_0=1$.

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


\bigskip
{\bf 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\ge0$. 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\equiv-u_1/3\pmod3$ if $u_1$ is divisible by 3 but $u_2$ is
not, and $A\equiv 4u_1+u_1u_2(2+3u_2+2u_1^2)\!\pmod9$ if $u_1$ is not divisible~by~3)
and the integrality of $u_4$ is equivalent to the simple congruence $u_1\equiv u_2\pmod2$.
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|\le250$, $0\le u_1\le100$, $|u_2|\le1000$  (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\mid u_n\notin\Z\}$ 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,\ldots,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$. 


\bigskip
{\bf 3. Analysis of the data.} Looking at the table, we observe that the solutions
found can be divided into five classes.
\smallskip\noindent
{\bf 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,\l)=(-1,0,d^2+d),\qquad u_n=\bn dn\bn{d+n}n\,,\tag4$$
where $d$ is a non-negative integer.  The corresponding function $F(t)=\sum u_nt^n$ is given by
$$ F(t)\=F(-d,1+d;1;t)\=P_d(1-2t)\,.\tag5$$   
Here $F(a,b;c;t)$ denotes the hypergeometric function and $P_d(X)$ the $d$th Legendre polynomial.
\smallskip\noindent
{\bf 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 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\le d\le7$ of an infinite family given by
 $$(A,B,\l)=(2,1,d^2+d+1),\qquad 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\,+\cdots+\,\bigl(\sum_{k=1}^d\frac2k\bigr)\,n\,+\,1\,.$$
Looking at the first few polynomials $H_d(x)$, we find that they satisfy 
$H_d(-x-1)=(-1)^dH_d(x)$ and $H_{2r}(-\frac12)=2^{-4r}\bn{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)=\sum u_nt^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 $\bn dn$, i.e., these solutions
are given by
$$F(t)\=\sum_{n=0}^\infty H_d(n)t^n\=\frac{F(-d,-d;1;t)}{(1-t)^{d+1}}
\=\frac1{1-t}\,P_d\bigl(\frac{1+t}{1-t}\bigr)\,,\tag6$$
where again $P_d(X)$ denotes the $d$th 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.
\smallskip\noindent
{\bf 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=\l=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,
$$ \text{Sequence \#2:}\quad(A,B,\l)=(0,-16,0),\qquad 
 u_n=\cases\bn{2r}r^2 &\text{if $n=2r$,}\\0&\text{if $n$ is odd.}\endcases$$
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.:
$$ \align \text{Sequence \#11:}&\quad(A,B,\l)=(16,0,4),\qquad 
   u_n\=\bn{2n}n^2\=\frac{(2n)!^2}{n!^4}\,,\\
\text{Sequence \#14:}&\quad(A,B,\l)=(27,0,6),\qquad 
   u_n\=\bn{3n}n\bn{2n}n\=\frac{(3n)!}{n!^3}\,,\\
\text{Sequence \#20:}&\quad(A,B,\l)=(64,0,12),\!\qquad 
   u_n\=\bn{4n}{2n}\bn{2n}n\=\frac{(4n)!}{(2n)!\,n!^2}\,,\endalign $$
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\to DA$, $\l\to D\l$, $u_n\to D^nu_n$ for some positive integer $D$ in order to achieve
integrality.  (The smallest choice of $D$ is $M^2\prod_{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\ge-1/2$.)



\smallskip\noindent
{\bf 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\mapsto Dt$, $u_n\mapsto D^nu_n$ with suitable $D\in\N$)
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\in\Q$ making the coefficients
$u_n=D^nH_d(n)$ integral for all $n$ is given by $\,M^2\prod_{p|M}p^*\,$, where
$M$ is the denominator of $d$ and $p^*=4$ if $p=2$, $p^*=p$ if $p$ is odd.
\smallskip\noindent
{\bf 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 {\bf A}--{\bf F}. We also add \#2 to
this list, with label {\bf G}, since although it is hypergeometric it does not belong to the 
1-parameter family of hypergeometric solutions discussed above.
\medskip \begintable 
new label & index | $\;A\;$ & $\;B\;$ & $\;\l\;$ | $\;u_0$ & $u_1$ & $u_2$ & $u_3$ & $u_4$ & $u_5$ & $u_6$\cr 
{\bf A} & \#5 | 7 & $-8$ & 2 | 1 & 2 & 10 & 56 & 346 & 2252 & 15184 \nr
{\bf B} & \#7 | 9 & 27 & 3 | 1 & 3 & 9 & 21 & 9 & $-297$ & $-2421$ \nr
{\bf C} & \#8 | 10 & 9 & 3 | 1 & 3 & 15 & 93 & 639 & 4653 & 35169 \nr
{\bf D} & \#9 | 11 & $-1$ & 3 | 1 & 3 & 19 & 147 & 1251 & 11253 & 104959 \nr
{\bf E} & \#10 | 12 & 32 & 4 | 1 & 4 & 20 & 112 & 676 & 4304 & 28496 \nr
{\bf F} & \#13 | 17 & 72 & 6 | 1 & 6 & 42 & 312 & 2394 & 18756 & 149136 \nr
{\bf G} & \#2 | 0 & $-16$ & 0 | 1 & 0 & 4 & 0 & 36 & 0 & 400 \nr
\endtable\noindent
Since all of these are contained in the domain $0<A\le17$, $0<u_1\le6$, $0<u_2\le42$,
which comprises less than 0,005\% of the domain of our search, one could hazard a rather optimistic
\proclaim{Conjecture} Up to normalizations, the only cases where the recursion $(3)$
has an integral solution are the hypergeometric family
$$ \text{\rm HG}_d\;:\qquad (A,\,B,\,\l)\= (-1,\,0,\,d^2+d)\qquad\quad(d\in\Q,\;d\ge-1/2) $$
(which for $d\in\Z_{\ge0}$ is terminating), the Legendre
 family
$$ \text{\rm Leg}_d\;:\qquad (A,\,B,\,\l)\= (2,\,1,\,d^2+d+1)\qquad\quad(d\in\Q,\;d\ge-1/2) $$
(which for $d\in\Z_{\ge0}$ is polynomial), and the seven sporadic solutions {\bf A--G}. \endproclaim


\bigskip
{\bf 4. Binomial coefficient sums.} The sequence ``D" in the table just given
is immediately recognized to be the Ap\'ery 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
$\dsize\frac{(3n)!}{n!^3}=\sum\limits_{k=0}^{2n}(-1)^{n-k}\bn{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 ({\bf G} is trivial):
$$\align \text{Sequence {\bf A}:}\quad& u_n=\sum_{k=0}^n\bn nk^3
  =\sum_{k=0}^n\bn nk^2\bn{2k}n\,,\\
\text{Sequence {\bf B}:}\quad& u_n=\sum_{k=0}^{[n/3]}(-1)^k\,3^{n-3k}\bn n{3k}\bn{3k}k\bn{2k}k\,,\\
\text{Sequence {\bf C}:}\quad& u_n=\sum_{k=0}^n\bn nk^2\bn{2k}k\,,\\ \allowdisplaybreak
\text{Sequence {\bf D}:}\quad& u_n=\sum_{k=0}^n\bn nk^2\bn{n+k}k
  =\sum_{k=0}^n(-1)^{n-k}\bn nk\bn{n+k}k^2 \,,\\
\text{Sequence {\bf E}:}\quad& u_n=\sum_{k=0}^{[n/2]}4^{n-2k}\,\bn n{2k}\bn{2k}k^2\,.\endalign$$
The formulas for the sequences {\bf B} and {\bf E} say that the corresponding generating
functions $\sum u_nt^n$ have hypergeometric representations:
$$ F_{\bold B}(t)\=\frac1{1-3t}\,F\bigl(\frac13,\frac23;1;\bigl(\frac{-3t}{1-3t}\bigr)^3\bigr)\,,
\qquad F_{\bold E}(t)\=\frac1{1-4t}\,F\bigl(\frac12,\frac12;1;\bigl(\frac{4t}{1-4t}\bigr)^2\bigr)\,.$$

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 \S7, a simple example (courtesy of the referee) being
the formula $u_n = \sum_{0\le j\le i\le n} (-1)^i 8^{n-i} \binom ni \binom ij^3$ in case {\bf F}. 

\bigskip
{\bf 5. Modular properties.} Thanks to Beukers [3], one knows that the Ap\'ery 
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))=\sum u_nt(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 {\it 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 {\it 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
\proclaim{Conjecture} Any integral solution of the differential equation~$(1)$,
where $P(t)$ is a non-degenerate quadratic polynomial, has a modular parametrization.
\endproclaim \noindent
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\le A\le54$, $0\le u_1\le21$,
$0<u_2\le500$, 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.

\smallskip\noindent
{\bf 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~\S5 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\pi i)^{-1}\,dt/dz$ is a (meromorphic) modular 
form of weight~2 and the function $2{f'}^2-ff''$ is a modular form of weight~6, so we can write
$$\frac{t'(z)}{f(z)^2}\=\a(t(z))\,,\qquad \frac{2f'(z)^2-f(z)f''(z)}{t'(z)f(z)^4}\=\b(t(z))\tag7$$
where $\a(t)$ and $\b(t)$ are algebraic (or even, if $t(z)$ is a Hauptmodul, rational)
functions of $t$. From the calculation
$$\frac1{t'(z)}\,\frac d{dz}\biggl(\frac{\a(t(z))}{t'(z)}\,\frac{dF(t(z))}{dz}\biggr)
+\b(t(z))\,F(t(z))\=\frac1{t'}\,\bigl(\frac1{f^2}\,f'\bigr)'+\frac{2{f'}^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 $(\a F')'+\b F=0$. Comparing this with~(1), we find that in our
case the function $\a(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)=\a(t(z))f(z)^2=\a(t(z))F(t(z))^2$
to get $z=\int dt/\a(t)F(t)^2$ or, in terms of the standard expansion parameter $q=e^{2\pi iz}$
of modular functions and forms at infinity, 
$$ q \= t\,\exp\biggl(\int_0^t\biggl(\frac1{\a(s)\,F(s)^2}-\frac1s\biggr)\,ds\biggr)
\= t\,\exp\biggl(\int_0^t\biggl(\frac{F(s)^{-2}}{1-As+Bs^2}-1\biggr)\,\frac{ds}s\biggr)\,.\tag8$$
This is a power series with leading term $t$ (recall that $F(t)=u_0+u_1t+\cdots$ 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 {\it 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.

\medskip\noindent
{\bf 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.
\smallskip\noindent{\bf Sporadic cases.}  In four of the six sporadic cases, those with labels
{\bf A}, {\bf C}, {\bf E} and {\bf F}, the quadratic polynomial $x^2-Ax+B$ splits over $\Z$ as
$(x-\a)(x-\b)$. In each of these cases, the modular functions $t(z)$, $1-\a t(z)$ and $1-\b 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^32^{-3}3^{-9}6^9$ is a shorthand for the
eta-product $\eta(z)^3\eta(2z)^{-3}\eta(3z)^{-9}\eta(6z)^9$ and the functions $\th_3$ and $\th_4$
are the Eisenstein or theta series
$$ \align \th_3(z) &\= 1+6\sum_{n=1}^\infty\biggl(\sum_{d|n}\bigl(\frac{-3}d\bigr)\biggr)\,q^n 
 \=\sum_{m,n\in\Z}q^{m^2+mn+n^2}\,,\\
\th_4(z) &\= 1+4\sum_{n=1}^\infty\biggl(\sum_{d|n}\bigl(\frac{-4}d\bigr)\biggr)\,q^n 
 \=\sum_{m,n\in\Z}q^{m^2+n^2}\,,\endalign$$
i.e., the modular forms of weight 1 whose Mellin transforms are the Dedekind zeta-functions of
the two imaginary quadratic fields $\Q(\sqrt{-3})$ and $\Q(i)$ possessing non-trivial units.
\medskip \begintable 
Case | $\;\a\;$ & $\;\b\;$ | $\;t(z)\;$ & $1-\a t(z)$ & $1-\b t(z)$ & $f(z)$ | $f(z)$ \cr 
{\bf A} | $-1$ & 8 | $\frac{1^36^9}{2^33^9}$ & $\frac{2^56^1}{1^13^5}$ & $\frac{1^86^4}{2^43^8}$ &
  $\frac{2^13^6}{1^26^3}$ | $\frac13\th_3(z)+\frac23\th_3(2z)$ \nr
{\bf C} | 1 & 9 | $\frac{1^46^8}{2^83^4}$ & $\frac{1^13^5}{2^56^1}$ & $\frac{1^96^3}{2^93^3}$
   & $\frac{2^63^1}{1^36^2}$ | $\frac12\th_3(z)+\frac12\th_3(2z)$  \nr
{\bf E} | 4 & 8 | $\frac{1^44^28^4}{2^{10}}$ & $\frac{1^44^{14}}{2^{14}8^4}$ & $\frac{1^84^4}{2^{12}}$ 
  & $\frac{2^{10}}{1^44^4}$ | $\th_4(z)$ \nr
{\bf F} | 8 & 9 | $\frac{1^53^14^56^212^1}{2^{14}}$ & $\frac{1^84^86^{20}}{2^{20}3^812^8}$ 
 & $\frac{1^94^96^6}{2^{18}3^312^3}$ & $\frac{2^{15}3^212^2}{1^64^66^5}$ |
  $\th_3(z)+2\th_3(2z)-2\th_3(4z)$ \nr
\endtable\noindent
Case {\bf B} is fairly similar: here the polynomial $x^2-Ax+B=x^2-9x+27$ does not split over $\Z$,
but we still have the eta-product representations 
$$\text{Case {\bf B}:}\qquad t\=\frac{1^34^318^9}{2^99^336^3}\,,\qquad 
1-9t+27t^2\=\frac{1^94^96^{36}9^336^3}{2^{27}3^{12}12^{12}18^9}\,,\qquad f\=\frac{2^93^112^1}{1^34^36^3}$$
and the Eisenstein series representation 
$$f(z)=\frac12\,\th_3(z)-\frac32\,\th_3(3z)-\th_3(4z)+3\,\th_3(12z)\,.$$
Finally, in the last case {\bf D} (corresponding to the Ap\'ery 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$:
$$ \text{Case {\bf D}:}\qquad t\= q\,\prod_{n=1}^\infty\bigl(1-q^n\bigr)^{5(\frac n5)}\,,\qquad
 1-11t-t^2\=t\cd\frac{1^6}{5^6}\,,\qquad f^2\=t^{-1}\cd\frac{5^5}{1^1}\;.$$
The Eisenstein series representation of $f$ in this case is also more complicated:
$$ f(z) \= 1\,+\,\sum_{n=1}^\infty\,\sum_{d|n}\biggl(\frac{3-i}2\,\chi(d)
  +\frac{3+i}2\,\bar\chi(d)\biggr)\,q^n\,,$$
where $\chi$ is the Dirichlet character of conductor 5 with $\chi(2)=i$, $\chi^2=\bigl(\frac\cd5\bigr)$.
\smallskip\noindent{\bf Hypergeometric cases.}  Here two cases, the exceptional case \#2 when
$A=\l=0$ and the sequence \#11 with $(A,B,\l)=(16,0,4)$, behave like the first four sporadic cases:
\medskip \begintable 
Case | $\;\a\;$ & $\;\b\;$ | $\;t(z)\;$ & $1-\a t(z)$ & $1-\b t(z)$ & $f(z)$ | $f(z)$ \cr 
{\#2} | $-4$ & 4 | $\frac{2^48^8}{4^{12}}$ & $\frac{2^{14}8^4}{1^44^{14}}$ & $\frac{1^42^28^4}{4^{10}}$ &
  $\frac{4^{10}}{2^48^4}$ | $\th_4(2z)$ \nr
{\#11} | 0 & 16 | $\frac{1^84^{16}}{2^{24}}$ & --- & $\frac{1^{16}4^8}{2^{24}}$
   & $\frac{2^{10}}{1^44^4}$ | $\th_4(z)$  \nr
\endtable \vskip-.4cm \noindent
(Note that, amusingly enough, $f(z)$ for case \#11 is the {\it same} as for the sporadic case
{\bf 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,\l)=(27,0,6)$, is more similar
to the sporadic case {\bf 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: 
$$ \text{Case \#14:}\qquad \frac t{1-27t}\=\frac{3^{12}}{1^{12}}\,,\qquad
   f^3\cd t\=\frac{3^9}{1^3}\,,\qquad f(z)=\th_3(z)\,.$$
Finally, for the sequence \#20 with $(A,B,\l)=(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:~
$$ \text{Case \#20:}\qquad \frac t{1-64t}\=\frac{2^{24}}{1^{24}}\,,\qquad
 f^4\cd t\=\frac{2^{16}}{1^8}\,,\qquad f(z )^2 \= -E_2(z)\,+\,2\,E_2(2z)\,,$$ 
where $E_2(z)=1-24\sum\sigma_1(n)\,q^n$ is the usual not-quite-modular Eisenstein series of
weight~2 on the full modular group.
\smallskip\noindent{\bf 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:
$$ \quad\text{Case \#19:}\quad\;\,t\=\frac{1^84^{16}}{2^{24}}\,,\qquad 1-16t\=\frac{1^{16}4^8}{2^{24}}\,,
\qquad f\=\frac{2^{22}}{1^{12}4^8}\,,\qquad f(z)\=\frac{-\th_4(z)+2\th_4(2z)}{1-16t(z)}\,.$$
$$ \text{Case \#25:}\qquad t\=\frac{1^{12}4^{12}6^{36}}{2^{36}3^{12}12^{12}}\,,\qquad 
(1-27t)^2f^3\=\frac{2^{27}3^312^3}{1^94^96^9}\,,\qquad f(z)\=\frac{-\th_3(z)+2\th_3(4z)}{1-27t(z)}\,.$$


\smallskip\noindent{\bf Discussion.} We make a few remarks about the modular
parametrizations just given.  We consider only the non-degenerate case when $B\ne0$ and $A^2\ne4B$,
so that $x^2-Ax+B=(x-\a)(x-\b)$ with $\a$ and $\b$ distinct and non-zero. We further assume
that $t$ is a Hauptmodul, i.e., that $t$ and $f$ are modular forms on some group $\G\subset SL(2,\Bbb R)$
of genus~0 and that $z\mapsto t(z)$ gives an isomorphism of $\,\frak H/\G\cup\{\text{cusps}\}\,$ to
$\Bbb P^1(\Bbb C)$. From equation (7) we see that the modular forms $f(z)$, $t(z)$ and $t'(z)$ are
related by $f^2=t'/t(1-\a t)(1-\b t)$, so that the $\G$-invariant differential form $f(z)^2dz$ is
(or pushes forward to) the differential form $dt/t(1-\a t)(1-\b t)$ on $\Bbb P^1$.  This form has
three simple poles, at $t=0$, $t=1/\a$ and $t=1/\b$, and one simple zero at $t=\infty$. 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+\cdots$).  None of the functions
$t(z)$, $1-\a t(z)$ or $1-\b 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 {\bf B} and {\bf D} when $\a$
and $\b$ are not in $\Z$, but must be interpreted as the statement that the logarithmic 
derivatives of $1-\a t$ and $1-\b 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.


\bigskip
{\bf 6. The periods associated to the six sporadic solutions.} The most famous property of the
original Ap\'ery recursion $(n+1)^2u_{n+1}-(11n^2+11n+3)u_n-n^2u_{n-1}=0$ is, of course, that 
if we start with the two linearly independent solutions $\{u_n\}=\{1,\,3,\,19,\,\ldots\}$ and
$\{v_n\}=\{0,\,1,\,\f{25}4,\ldots\}$, then the ratio $v_n/u_n$ tends to the limit $L=\f15\z(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 {\bf A}--{\bf 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\'ery equation. 

We denote by $\{u_n\}$, as up to now, the solution of equation~(3) with $u_0=1$ (and hence
$u_1=\l$, $u_n\in\Z$ 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\ge1$. 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 $\a$ and $\b$ as in \S5 denote the roots of $x^2-Ax+B=0$, chosen (except in case {\bf B}
where $\a,\,\b=(9\pm3i\sqrt3)/2$) so that $|\b|>|\a|$.  Then the sequence $\{u_n\}$ has the
asymptotic behavior
$$  u_n\= C\,\f{\b^n}{n}\,\bigl(1-\f{c_1}{n}+\text O\bigl(\f1{n^2}\bigr)\bigr)
  \qquad\qquad(n\to\infty)\,,\tag 9$$
where $c_1=\dfrac{\l-\a}{\b-\a}$, while every other element of the 2-dimensional vector space of
real-valued solutions of the recursion (3) (for $n\ge n_0\ge1$) has asymptotics given either by the
same formula (with a possibly different value of $C$) or else by 
$\dsize\,(\text{const.})\f{\a^n}{n}\,\bigl(1-\f{c_2}{n}+\text O(\f1{n^2}\bigr)\bigr)$,
where  $c_2=\dfrac{\l-\b}{\a-\b}\,$, 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=(\text{const.})\,n^\v x^n(1-cn^{-1}+\text O(n^{-2}))$; then the ratio 
$\rho_n:=u_{n+1}/u_n$ equals $x\bigl(1+\v n^{-1}+\bigl(\bn\v2+c\bigr)n^{-2}+\ldots\bigr)$, and
rewriting equation~(3) in the form $(1+n^{-1})^2\rho_n-A(1+n^{-1})-\l n^{-2}+B\rho_{n-1}^{-1}=0$
we find $x\in\{\a,\b\}$, $\v=-1$ and $c=(\l-A+x)/(2x-A)$.)  It follows that there are real constants
$C'$ and $L$ such that $v_n-Lu_n=C'\f{\a^n}n(1-\f{c_2}n+\text O(\f1{n^2}\bigr)\bigr)$ and that
the quotient $v_n/u_n$ tends to $L$ exponentially quickly:
$$ \f{v_n}{u_n}\=L\,+\,C''\,\bigl(\f\a\b\bigr)^n\,\bigl(1+\f{c_1-c_2}n
   +\text O\bigl(\f1{n^2}\bigr)\bigr)  \qquad\qquad(n\to\infty)\,.\tag 10$$
The values of the constants $\b$, $\a$, $C$, $c_1$, $c_2$ and $L$ occurring in equations (9) and (10)
are given for each of the five cases {\bf A}, {\bf C}, {\bf D}, {\bf E} and {\bf F} in the table
below, in which $\phi$ in the third row denotes the golden ratio $(1+\sqrt 5)/2$ and the numbers
$\z(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,\chi)$ associated to the quadratic characters $\chi(n)=\bigl(\f{-3}n\bigr)$
and $\chi(n)=\bigl(\f{-4}n\bigr)$, each of which is a famous number: $\z(2)$ equals $\pi^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~{\bf F}) in \S4 or else from the modular form expressions 
of \S5. The evaluation of the limiting ratio $L$ will be discussed at the end of this section.


\medskip \begintable 
Case | $\;\b\;$ & $\;\a\;$ | $\;\,C\;$ | $\;c_1\;$ & $\;c_2\;$ | $L$ \cr 
{\bf A} | 8 & $-1$ | $\f2{\pi\sqrt3}$ |  $\f13$ & $\f23$ | $\;\f14\z(2)=0.4112335\cdots$ \nr
{\bf C} | 9 & 1 | $\f{3\sqrt3}{4\pi}$ |  $\f14$ & $\f34$ | $\;\f12L_{-3}(2)=0.3906512\cdots$ \nr
{\bf D} | $\;\phi^5$ & $-\phi^{-5}$ | $\;\f{\phi^{5/2}}{2\pi\root4\of5}$ | 
    $\,\f1{\phi\sqrt5}$ & $\f\phi{\sqrt5}$ | $\;\f15\z(2)=0.3289868\cdots$ \nr
{\bf E} | 8 & 4 | $\f2\pi$ | $0$ & $1$ | $\;\f12L_{-4}(2)=0.4579827\cdots$ \nr
{\bf F} | 9 & 8 | $\f{3\sqrt3}\pi$ |  $-2$ & $3$ | $\;\f58L_{-3}(2)=0.4883140\cdots$ \nr
\endtable\noindent

In each of these five cases, formula (10) gives an interesting series of rational approximations for one
of the numbers $\z(2)$, $L_{-3}(2)$ or $L_{-4}(2)$, but, as already mentioned, except in case {\bf D}
these do not converge quickly enough to yield the irrationality of the limit.  Indeed, from the
table we see that in cases {\bf E} and {\bf F} the quantity $v_n-Lu_n$ blows up exponentially like
$4^n/n$ or $8^n/n$, respectively, and even in cases {\bf A} and {\bf C}, where $|\a|=1$ and hence
$v_n-Lu_n$ tends to zero like $\,\text 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 {\bf E} (Catalan's constant) and {\bf 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 {\bf C}
and {\bf E} being 
$$ L_{-3}(2)\=\cfrac2\\P(0)-\cfrac{9\cd1^4}\\P(1)-\cfrac{9\cd2^4}\\P(2)-
   \cfrac{9\cd3^4}\\{\ddots}\endcfrac\;, \qquad L_{-4}(2)\= \cfrac{1/2}\\Q(0)-\cfrac{2\cd1^4}
   \\Q(1)-\cfrac{2\cd2^4}\\Q(2)- \cfrac{2\cd3^4}\\{\ddots}\endcfrac$$
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\'ery's irrationality proofs for $\z(2)$ and $\z(3)$,
using the modular descriptions of the differential equations associated to the recursion~(3)  which
were given in~\S5.  Consider for instance case {\bf C}.  We saw that the integer solution $\{u_n\}$
of~(3) in this case is given by the generating function expansion $f(z)=\sum u_nt(z)^n$, where 
$t(z)=\eta(z)^4\eta(6z)^8/\eta(2z)^8\eta(3z)^4$ is the ``Hauptmodul" for $\G_0(6)$ and $f(z)$ is
the theta series (or Eisenstein series) of weight 1 for $\G_0(6)$ given in~\S5.  The generating function
$h(z)=\sum v_nt(z)^n$ for the second solution of the recurrence~(3) has a more complicated modular
interpretation as $f(z)\tg(z)$, where $\tg(z)=q-\f54q^2+q^3-\f{11}{16}q^4+\cdots$ ($q=e^{2\pi iz}$)
is the ``Eichler integral" $\sum b(n)q^n/n^2$ associated to the weight~3 Eisenstein series
$g(z)=\sum b(n)q^n=q-5q^2+9q^3-11q^4+\cdots$ with $b(n)=\sum_{d|n}(-1)^{d-1}\bigl(\f{-3}{n/d}\bigr)d^2$.
(To prove this, one simply verifies that the differential equation satisfied by $f\tg$ as a function
of $t(z)$ coincides with the one corresponding to the recursion defining the $v_n$.)  The group
$\G_0(6)$ has four cusps, at $z=0$, $\f12$, $\f13$ and $\infty$, the corresponding values of the
modular function $t(z)$ being $\f19$, $\infty$, 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\to0$ (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\to0$ of the ratio
$h(z)/f(z)=\tg(z)$.  (Note that, since $\sum\bigl(v_n-Lu_n)t^n$ has a larger radius of convergence
than $\sum v_nt^n$ or $\sum u_nt^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 $\tg(z)$ as $\dsize\sum\limits_{m\ge1}\f{(\f{-3}m)}
{m^2}\f{q^m}{1+q^m}$, one sees that its limit for $z\to0$ or $q\to1$ equals $\f12L_{-3}(2)$, as claimed.
The proof in all the other cases is similar, the function $\tg(z)$ in case {\bf E}, for instance, 
being the Eichler integral $\dsize\sum\limits_{m\ge1}\f{(\f{-4}m)}{m^2}\f{q^m}{1+q^m}$ of a weight~3
Eisenstein series on $\G_0(8)$.

Finally, we say a few words about case {\bf B}, which has been omitted up to now.  Here the two roots
$\a=(9+3\sqrt{-3})/2$ and $\b=\bar\a$ of $x^2-Ax+B=0$ have the same absolute value, so all solutions
of~(3) have the same growth $\,\text O(27^{n/2}\!/n)$, and the numbers $v_n/u_n$ do not tend to a limit 
as $n\to\infty$.  Instead, the asymptotic expansions of $\{u_n\}$ and $\{v_n\}$ as $n\to\infty$ are given by
$$ \align u_n&\=\f Un\,\a^n\,\bigl(1-\f{\g_1}n+\cdots\bigr)
  \,+\,\f{\bar U}n\,\bar\a^n\,\bigl(1-\f{\bar\g_1}n+\cdots\bigr)\,,\\
  v_n&\=\f Vn\,\a^n\,\bigl(1-\f{\g_1}n+\cdots\bigr)
\,+\,\f{\bar V}n\,\bar\a^n\,\bigl(1-\f{\bar\g_1}n+\cdots\bigr)\,,\endalign$$
where $U=\dfrac{9\sqrt3}{2\pi\a}$, $\dsize\g_1=\f\a9$ and $\dsize \f VU=\f12L_{-3}(2)-\f{4i}{9\sqrt3}\z(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/\a$ and $1/\bar\a$ are the values of the Hauptmodul $t(z)$ given in \S5 at the cusps $z=\pm1/6$ and the function
$\tg(z)=h(z)/f(z)$ has the form $\tg(z)=-\sum_{n=1}^\infty\bigl(\f{-3}n\bigr)n^{-2}(-q)^n/(1-(-q)^n)$, which for
$z=1/6$, $-q=e^{-2\pi i/3}$ has the value given.)   Since $\a/\sqrt{27}=e^{\pi i/6}$, the quotients $v_n/u_n$
do tend to a limit $\f12L_{-3}(2)+\f4{9\sqrt3}\z(2)\tan\f{\pi(j-1)}6$ as $n$ tends to infinity in
any fixed residue class $n\equiv j\not\equiv4$ (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.

\bigskip
{\bf 7. The modular cases as Picard-Fuchs differential equations.}  The fact that the
non-degenerate 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 $\Bbb P^1$
having exactly four singular fibres (a stable family of elliptic curves over $\Bbb P^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(\Z)$
describing the family. We reproduce the table from~[2] here.


\medskip \begintable 
Case | Equation of the family | Group \cr 
I | $X^3+Y^3+Z^3+tXYZ=0$ | $\G(3)$ \nr 
II | $X(X^2+Z^2+2ZY)+tZ(X^2-Y^2)=0$ | $\G_1(4)\cap\G(2)$ \nr
III | $X(X-Z)(Y-Z)+tZY(X-Y)=0$ | $\G_1(5)$ \nr
IV | $(X+Y)(Y+Z)(Z+X)+tXYZ=0$ | $\G_1(6)$ \nr
V | $(X+Y)(XY-Z^2)+tXYZ=0$ | $\G_0(8)\cap\G_1(4)$ \nr
VI | $X^2Y+Y^2Z+Z^2X+tXYZ=0$ | $\G_0(9)\cap\G_1(3)$ \endtable  \noindent

Comparing the number six of cases found by Beauville with the number six of sporadic cases {\bf A--F} in \S5
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 {\bf 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 {\it 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 \S5, 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 $\infty$, 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 $\infty$, giving 12 choices (the remaining freedom of interchanging or
rescaling $t_1$ and $t_2$ has no effect on the parameters $A$, $B$, $\l$ 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:   \def\p{\phantom{X}}
\
\smallskip
\begintable Beauville Family: | IV | \p I, VI \p \hfil | \p III \p |  \p II, V \p \cr
Picard-Fuchs equation: | \p {\bf A}, {\bf C}, {\bf F}  \p | \p {\bf B}, \#14 \p | {\bf D} | \p {\bf E}, {\bf G}, \#11 \endtable

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:
\medskip \begintable 
  | Transformation to Weierstrass form | $j$-invariant \cr 
I | $X,\,Y=t(x\pm y),\;Z=2(1+3x)$ | $-t^3(t^3-216)^3/(t^3+27)^3$ \nr
II | $X=xt,\;Y=x+y,\;Z=t$ | $2^8(t^4-t^2+1)^3/t^4(t^2-1)^2$ \nr
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)$ \nr
IV | $X,\,Y=1+tx\pm y,\;Z=2x(1+tx)$ | $(t+2)^3(t^3+6t^2-12t+8)^3/t^3(t-1)^2(t+8)$ \nr
V | $X,\,Y=1+tx\pm y,\;Z=2x(1+tx)$ | $(t^4+16t^2+16)^3/t^2(t^2+16)$ \nr
VI | $X=2x,\; Y=y-xt-1,\;Z=2x^2$ | $-t^3(t^3+24)^3/(t^3+27)$ \nr
\endtable\noindent

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\"obius transformation
$z\mapsto M(z)$, $M\in\text{GL}_2(\Z)$, of the $z$-variable) up to a fractional linear transformation. 
For instance, from the above table one sees that {\bf A} and {\bf C} correspond to the same fibration, and
from the eta-product expansions given in \S5 one finds that $t_{\bold C}(z)=t_{\bold A}(z)/(1+t_{\bold 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\mapsto t^N$.  Since these
transformations change the nature and number of the singularities, one can get examples of Ap\'ery-like
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~\S4) 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 $n$th power of
$-(X^3+Y^3+Z^3)/XYZ$, 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)/XYZ$, corresponding to the family~IV, are the 
numbers 1, 2, 10, 56, \dots\ of sequence~\#5 in our list.



\bigskip 
{\bf 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 \S6 on the periods associated to Ap\'ery-like
recursions (which was written in answer to a question posed by A.~Connes during my course at the Coll\`ege de
France in 2001) and some additions to the discussion of the connection with Beauville's classification in \S7.  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 {\bf D}, {\bf A} and {\bf C} (=  \#9, \#5 and \#8) are discussed
in the paper [7] by Stienstra and Beukers, where two of them are ascribed to Ap\'ery 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 \S5 is not complete.  No further
sporadic Ap\'ery-like equations have been discovered, and the conjecture formulated in \S3 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,\l)=(-1,0,d(d+1))$  ($d\in\Q$, $d\ge-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,\l)=(2,1,d(d+1)+1)$  ($d\in\Q$, $d\ge-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\'ery-like 
differential equations similar to the one given in~\S6 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.

\vskip2cm

\Refs
\smallskip

\ref\key 1 \by G.~Almkvist and W.~Zudilin, {\it Differential equations, mirror maps and zeta values.}
In {\it Mirror symmetry.~V},  AMS/IP Stud. Adv. Math. {\bf 38}, Amer. Math. Soc. Providence, RI (2006), 481--515 \endref 
\smallskip
\ref \key 2 \by A. Beauville, {\it Les familles stables de courbes elliptiques sur $\bold P^1$ admettant quatre
fibres singuli\`eres.} C. R. Acad. Sc. Paris {\bf 294} (1982), 657--660\endref 
\smallskip
\ref \key 3 \by F. Beukers, {\it Irrationality proofs using modular forms.}  Journ\'ees Arithm\'etiques
de Besan\c con. Ast\'erisque {\bf 147--148} (1987), 271--283\endref 
\smallskip
\ref \key 4 \by F. Beukers, {\it On B.~Dwork's accessory parameter problem.} Math.~Z.~{\bf 241} (2002), 425--444. \endref 
\smallskip
\ref\key 5 \by V. Golyshev, {\it Classification problems and mirror duality.} Preprint (2005), arXiv:math/0510287 
\endref
\smallskip
\ref\key 6 \by U.~Schmickler-Hirzebruch, {\it Elliptische Fl\"achen \"uber $\Bbb P_1\Bbb C$ mit drei Ausnahmefasern und
die hypergeometrische Differentialgleichung.}  Schriftenreihe des Mathematischen Instituts der Universit\"at M\"unster {\bf 33} (1985).
\endref
\smallskip
\ref\key 7 \by J.~Stienstra and F.~Beukers, {\it On the Picard-Fuchs equation and the formal Brauer group
of certain elliptic $K3$-surfaces.}  Math.~Ann.~{\bf 271} (1985), 269--304\endref
\smallskip
\ref\key 8 \by J.~Stienstra, {\it Mahler measure variations, Eisenstein series and instanton expansions.}
In {\it Mirror symmetry.~V},  AMS/IP Stud. Adv. Math. {\bf 38}, Amer. Math. Soc. Providence, RI (2006), 139--150. \endref 
\smallskip
\ref \key 9 \by H. Verrill, {\it Picard-Fuchs equations of some families of elliptic curves.} Proceedings
  on Moonshine and Related Topics (Montreal, 1999), CRM Proc.~Lecture Notes {\bf 30}, AMS (2001), 253--268 \endref 
\smallskip
\ref \key 10 \by Y.~Yang, {\it Ap\'ery limits and special values of $L$-functions}.  Preprint (2007), arXiv:0709:1968v1
\endref\smallskip
\ref \key 11 \by D.~Zagier, {\it Elliptic modular forms and their applications.} In {\it 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
\endref
\smallskip
\ref \key 12 \by W. Zudilin, {\it Ap\'ery-like difference equation for Catalan's constant.} 
Electron.~J.~Combin.~{\bf 10} (2003), \#R14, 10~pp \endref 

\end
