The Apéry numbers are famous for proving that $\zeta(3)$ is irrational. Less well known is that they also give a way to compute $\zeta(3)$ from a polynomial ODE. The key is an adaptation law that tracks the ratio of two blowing-up solutions — and a subtle inhomogeneous term that took several failed attempts to identify. Here is the full story, including the failures.

The Challenge

We want to construct a polynomial ODE system — one whose right-hand sides are polynomials and whose initial conditions are rational numbers — such that one of its variables converges to $\zeta(3) = 1.2020569\ldots$ as $t \to \infty$.

This question arises from the theory of Chemical Reaction Networks (CRNs). The mass-action kinetics of any CRN is a polynomial ODE, and the class of real numbers computable as limits of such ODEs ($R_{\mathrm{RTCRN}}$) is a natural object of study. Can $\zeta(3)$ be computed this way?

The Apéry Sequences

In 1978, Roger Apéry stunned the mathematical world by proving that $\zeta(3) = \sum_{k=1}^{\infty} 1/k^3$ is irrational. His proof rests on a pair of integer sequences.

A companion sequence $b_n$ satisfies the same recurrence but with different initial conditions: $b_0 = 0$, $b_1 = 6$. The sequence $b_n$ involves sums weighted by partial sums of $1/k^3$: roughly, $b_n$ is the “harmonic companion” of $a_n$.

Apéry’s fundamental insight: the ratio $b_n/a_n$ converges to $\zeta(3)$ as $n \to \infty$, and the convergence is fast enough (exponential) to prove irrationality via the standard criterion. Explicitly, $|b_n/a_n - \zeta(3)| = O((\sqrt{2}-1)^{8n})$.

From Sequences to Generating Functions

Here $z$ is the generating function parameter — a real variable that will later serve as the independent variable of an ODE.

(Here $y = y(z)$, and primes denote $d/dz$.)

is called the conifold point. It is the radius of convergence of $A(z)$: the power series converges for $|z| < z_1$ and diverges beyond. Since $a_n$ grows like $(1/z_1)^n$, the function $A(z)$ blows up as $z \to z_1^-$. The same is true for $B(z)$.

As we will discover, $A(z)$ satisfies the homogeneous Apéry ODE above, but $B(z)$ does not — it satisfies an inhomogeneous variant. This distinction turned out to be the crux of the entire problem.

A Remark on the Name “Calabi-Yau”

Readers browsing the literature may encounter this ODE under the heading of Calabi-Yau differential operators. The name sounds intimidating, but the connection is a historical accident of discovery, not a prerequisite for understanding.

In the 1990s, physicists studying mirror symmetry for Calabi-Yau manifolds encountered certain differential equations — called Picard-Fuchs equations — that govern how geometric periods vary across moduli spaces. Almkvist, van Enckevort, van Straten, and Zudilin later catalogued all such fourth-order operators with integral power-series solutions, creating the AESZ database. The Apéry operator for $\zeta(3)$ appears as AESZ #1: the simplest entry.

What does this mean for us? Very little, in practice. For computing $\zeta(3)$, we need only three facts about the ODE:

1. It has a MUM point (maximally unipotent monodromy) at $z=0$, where all three indicial exponents coincide. This is what produces the logarithmic solutions and, ultimately, the transcendental number $\zeta(3)$ in the monodromy data.
2. The conifold point $z_1 = 17 - 12\sqrt{2}$ is a regular singularity where $A''$ and $B''$ diverge but their ratio converges.
3. The polynomial $p(z) = z^2(1-34z+z^2)$ that appears as the leading coefficient is all we need for the $\tau$-time reparametrization.

The Calabi-Yau geometry explains why the Apéry ODE has this beautiful structure — the integrality of $a_n$, the special monodromy, the connection to modular forms. But the computation stands entirely on its own: sequences, generating functions, ODE, PIVP.

The Idea: From Error Function to Adaptation Law

The Apéry theorem says $b_n/a_n \to \zeta(3)$ for the coefficients. Can the generating functions $A(z)$ and $B(z)$ be used to extract $\zeta(3)$ as well?

The error function $E(z)$

If $\zeta$ happens to equal $\zeta(3)$, then $E$ is the difference between $B$ and the “correct multiple” of $A$. The question is: what is special about $E$ when $\zeta = \zeta(3)$?

Local solutions at the conifold

The conifold $z_1$ is a regular singular point of the Apéry ODE. What does this mean? The leading coefficient $p(z) = z^2(1-34z+z^2)$ vanishes at $z_1$, so dividing by $p(z)$ introduces a pole. But the pole is “mild” (at most first-order after normalization), which makes $z_1$ a regular singular point, as opposed to an irregular one.

At a regular singular point, there is a standard technique — the Frobenius method — for finding local solutions. The idea: try a solution of the form $y = (z_1-z)^{\rho}\,(a_0 + a_1(z_1-z) + \cdots)$ and substitute into the ODE. The lowest-order terms yield a polynomial equation in $\rho$ alone, called the indicial equation. Its roots — the local exponents — tell us what singularities the solutions can have.

For the Apéry ODE at $z_1$, the indicial equation turns out to be $\rho(2\rho - 1)(\rho - 1) = 0$, giving three roots: $\rho = 0$, $\rho = 1/2$, and $\rho = 1$. (This can be verified by direct computation: substitute $y = (z_1-z)^{\rho}$ into the ODE, expand near $z = z_1$, and set the leading coefficient to zero.)

Each exponent corresponds to a local solution of the form $(z_1-z)^{\rho} \cdot (\text{analytic})$:

The other two solutions $\varphi_0, \varphi_2$ and all their derivatives remain bounded at $z_1$.

Why the singular parts cancel

At this value of $\zeta$, the error function $E$ is a combination of $\varphi_0$ and $\varphi_2$ only, both analytic at $z_1$. So $E(z;\,\beta_1/\alpha_1)$ is analytic at $z_1$.

The punchline: $\beta_1/\alpha_1 = \zeta(3)$. This is not a coincidence — it is the analytic content of Apéry’s irrationality proof, translated into the language of the ODE’s singular structure.

Why $B''/A'' \to \zeta(3)$

They blow up at the same rate — the dominant $(z_1-z)^{-3/2}$ factor is identical — so their ratio converges:

because $E''$ stays bounded (analytic) while $A'' \to \infty$. In short: $A''$ and $B''$ both explode because of the same $\sqrt{z_1-z}$ singularity, and the ratio of the explosions is $\zeta(3)$.

The adaptation law

This analysis suggests a dynamical mechanism: introduce a variable $\zeta(\tau)$ that tracks $B''/A''$ in real time. Define

(Here $\tau$ is a reparametrized time variable; the reason for using $\tau$ instead of $z$ directly will be explained below.)

This is an instance of the classical method of unknowns (variation of parameters): we promote the constant $\zeta$ to a function $\zeta(\tau)$ and design an ODE that drives it toward the correct value. The update is gradient-descent style: $\zeta$ moves toward the value that makes $B'' - \zeta A'' = 0$, i.e., the instantaneous equilibrium is $\zeta^* = B''(z)/A''(z)$. Crucially, the right-hand side $\varepsilon(A'' B'' - \zeta (A'')^2)$ is polynomial in the state variables.

The stability is easy to check: linearizing around $\zeta = B''/A''$, the Jacobian is $-\varepsilon (A'')^2 < 0$, so the equilibrium is attracting. As $z \to z_1$, we showed above that $B''/A'' \to \zeta(3)$, so the adaptation variable $\zeta$ converges to $\zeta(3)$.

The key remaining question: what ODE does $B(z)$ actually satisfy?

The Failed Attempts

Getting the adaptation law to work required identifying the correct ODE for $B(z)$. Here is what went wrong along the way.

Attempt 1: Assuming $B$ Satisfies the Homogeneous ODE (Failed)

The most natural assumption: since $a_n$ and $b_n$ satisfy the same recurrence, perhaps $A(z)$ and $B(z)$ satisfy the same ODE.

Under this assumption, we evolve both $A$ and $B$ using the homogeneous Apéry ODE and compute $B''(z)/A''(z)$ as $z \to z_1$.

Result: $B''/A'' \to 0.0296 \approx z_1$, nowhere near $\zeta(3) = 1.202$. The adaptation law faithfully tracks this ratio — and converges to the wrong value.

Diagnosis: The assumption is false. $B(z)$ does not satisfy the homogeneous Apéry ODE. (The same recurrence for coefficients does not imply the same ODE for generating functions, because the generating function ODE also encodes the initial conditions through its singularity structure.)

Attempt 2: Using $A$ and $B$ Directly Instead of Derivatives (Failed)

Result: Stable (the Jacobian at equilibrium is $-\gamma A^2 < 0$), but it converges to $B(z_1)/A(z_1) \approx 0.018$, not $\zeta(3)$.

Diagnosis: The ratio $B(z)/A(z)$ does not converge to $\zeta(3)$ at the conifold. The Apéry theorem $b_n/a_n \to \zeta(3)$ is a statement about coefficients (the asymptotic ratio of the $n$-th terms), not about function values (the ratio of summed series). These are fundamentally different — the generating functions weight the coefficients by $z^n$, which alters the limiting ratio.

Attempt 3: Wrong Sign (Failed)

Result: The Jacobian at equilibrium is $+\gamma A^2 > 0$, so the equilibrium is repelling instead of attracting. $\zeta$ blows up.

Lesson: The sign of the adaptation law matters. The update $d\zeta/d\tau \propto +f \cdot (g - \zeta f)$ gives a stable equilibrium at $\zeta = g/f$; a minus sign makes it unstable.

The Breakthrough: The Inhomogeneous Term

The failures above all stem from one source: using the wrong ODE for $B(z)$. What ODE does $B(z)$ actually satisfy? We can check numerically: compute $b_n$ as exact rational numbers (using the recurrence), form the truncated power series $B(z)$, substitute into the Apéry differential operator, and see what comes out.

using the truncated power series for $B(z)$ with exact rational coefficients.

(Why 6? The constant comes from $b_1 = 6$: the initial condition for the companion sequence injects an inhomogeneity. One can verify this algebraically: if $B = \sum b_n z^n$ where $b_n$ satisfies the Apéry recurrence with $b_0 = 0, b_1 = 6$, then the action of the Apéry differential operator on $B$ telescopes to a boundary term determined by $b_0$ and $b_1$, yielding the constant 6.)

This changes everything. With the corrected inhomogeneous ODE for $B$, the second-derivative ratio converges rapidly:

$B''(z)/A''(z) \to \zeta(3)$ as $z \to z_1$.

The numerical table confirms this convergence with the corrected inhomogeneous ODE.

The Complete Polynomial ODE System

Since $p(z_1) = 0$, the conifold $z_1$ is a fixed point of this ODE: $z$ approaches $z_1$ exponentially in $\tau$-time ($z_1 - z \sim e^{-\lambda\tau}$ as $\tau \to \infty$). This gives the adaptation law infinite time to converge.

The chain rule converts $z$-derivatives to $\tau$-derivatives: $\dot{y} = y' \cdot p(z)$ (where dots are $d/d\tau$ and primes are $d/dz$). For the third derivative, $\dot{y}'' = y''' p(z)$, and we eliminate $y'''$ using the ODE itself: $y''' = (-q\,y'' - r\,y' - s\,y)/p$ for $A$, and $y''' = (-q\,y'' - r\,y' - s\,y + 6)/p$ for $B$. After cancellation of $p$, the result is polynomial.

The complete 8-variable polynomial ODE:

Initial conditions at $z_0 = 1/1000$: evaluate the Taylor series with exact rational coefficients to get $A(z_0)$, $A'(z_0)$, $A''(z_0)$, $B(z_0)$, $B'(z_0)$, $B''(z_0)$ — all rational numbers. Set $\zeta_0 = 0$ (no initial knowledge of $\zeta(3)$).

All right-hand sides are polynomial in the 8 state variables. All initial conditions are rational.

Numerical Verification

Starting from $\zeta_0 = 0$, the adaptation variable tracks $B''/A''$ and converges to $\zeta(3)$:

The result is robust: all values of $\varepsilon$ from $10^{-6}$ to $10^{-2}$ give the same 14-digit answer. The solver (Radau, implicit Runge-Kutta) uses 2579 steps.

The following plot shows the dynamics of the unbounded system. On the left, all variables blow up exponentially as $z \to z_1$. On the right, the ratio $B''/A''$ converges to $\zeta(3)$ — the blow-up is “proportional,” with $\zeta(3)$ as the proportionality constant.

The Bounded System: $\zeta(3) \in R_{\mathrm{RTCRN}}$

The 8-variable system above computes $\zeta(3)$, but the variables $A', A'', B', B''$ all diverge as $z \to z_1$. For membership in $R_{\mathrm{RTCRN}}$ (the class of CRN-computable reals), all variables must remain bounded.

The fix is simple: divide everything by $A''$ (the fastest-growing variable) to form ratio variables:

All six are bounded: $\alpha, \sigma_A, \beta, \sigma_B, w \to 0$ as $\tau \to \infty$, and $\rho \to \zeta(3)$.

Applying the quotient rule to each ratio (e.g., $\dot\alpha = (\dot{A}' \cdot A'' - A' \cdot \dot{A}'')/(A'')^2$) and substituting the $\tau$-time equations, all factors of $A''$ cancel. Define the common subexpression $Q = q + r\alpha + s\sigma_A$. The resulting system is:

The adaptation law simplifies from $\dot\zeta = \varepsilon A''(B'' - \zeta A'')$ to $\dot\zeta = \varepsilon(\rho - \zeta)$. Since $\rho$ itself converges to $\zeta(3)$, the simple exponential tracking $\zeta \to \rho \to \zeta(3)$ suffices — the unbounded forcing $(A'')^2$ in the original law is no longer needed.

All right-hand sides are polynomial. All initial conditions are rational. All variables are bounded. Numerically, $\zeta$ converges to $\zeta(3)$ with error $4.4 \times 10^{-16}$ (machine precision) in 4985 steps.

The contrast between the two systems is striking. The bounded system’s variables all remain order-1 throughout, with the auxiliary ratios $\alpha, \sigma_A, \beta, \sigma_B, w$ decaying to zero while $\rho$ and $\zeta$ converge to $\zeta(3)$:

And here is the side-by-side comparison, showing how the same number $\zeta(3)$ emerges from two very different dynamical regimes — one with exponential blow-up, one with everything bounded:

This establishes:

Further Reading

• R. Apéry, Irrationalité de $\zeta(2)$ et $\zeta(3)$, Astérisque 61 (1979), 11–13.
• A. van der Poorten, A proof that Euler missed… Apéry’s proof of the irrationality of $\zeta(3)$, Math. Intelligencer 1 (1979), 195–203.
• H.-L. Chen and X. Huang, Bounded Analog Complexity (2026).

The Apéry sequences contain $\zeta(3)$ in the ratio of their asymptotic growth rates. The adaptation law is the dynamical mechanism that extracts it. The hidden constant 6 in the inhomogeneous ODE was the key that unlocked the computation — without it, the ratio converges to a meaningless number.

三川汇海知归处，百转千回终见真。
Three streams meet the sea and know their destination; a hundred turns, a thousand detours, and truth is found at last.