The previous post computed Apéry’s constant $\zeta(3)$ via the Fermi integral

That evaluation was classical analysis. This post asks a different question: can a continuous-time analog machine — a polynomial initial-value problem (PIVP), the mathematical heart of a Chemical Reaction Network — compute $\zeta(3)$? And if so, with what resources?

Two independent encodings present themselves:

1. The integral route. Turn the Fermi integral above into a bounded ODE system and integrate.
2. The series route. Start from Apéry’s series $\zeta(3) = \tfrac{5}{2}\sum_{n \geq 1}(-1)^{n-1}/(n^3\binom{2n}{n})$ — which has a holonomic generating function — and encode the generating function’s ODE.

The integral route works, and this post gives the explicit 5-variable bounded polynomial system. The series route also works numerically, but reveals a genuine obstruction rooted in the analytic structure of the generating function at the origin. That obstruction is not a bookkeeping nuisance; it is the same object that blocks the sharpest irrationality and transcendence questions about $\zeta(3)$.

Both constructions are written out in detail below, with the aim of laying groundwork for future revisits.

1. What “Encoding” Means Here

A PIVP in dimension $d$ is an autonomous system

with $p$ a vector of polynomials in the state variables and rational coefficients. A real number $\alpha \in \mathbb{R}$ is PIVP-computable if some component $y_i(t)$ of some such system satisfies $y_i(t) \to \alpha$ as $t \to \infty$ along a trajectory that remains bounded on $[0, \infty)$.

The model was introduced by Shannon (1941) as the GPAC — General-Purpose Analog Computer — and later refined by Graça, Pouly, and others. On the CRN (Chemical Reaction Network) side, Fages–Gay–Soliman and subsequent work showed that polynomial ODEs in nonnegative variables correspond precisely to CRNs via dual-rail encoding, so “PIVP-computable” and “CRN-computable” are essentially the same class.

The interesting question is not merely which real numbers are PIVP-computable — that class is enormous — but rather at what cost. The cost is measured by a time modulus $\mu(r)$: the time required to be within $2^{-r}$ of the target. When $\mu(r) = \Theta(r)$ — linear in the bit count of precision — the number is called real-time computable. This is the first floor of a larger complexity hierarchy, and it is where constants like $e$, $\pi$, and algebraic numbers live.

The open question at the center of this post: does $\zeta(3)$ live on the first floor?

2. The Integral Route — A 5-Variable Bounded PIVP

Start from the classical identity

rewritten after absorbing the rational prefactor $\tfrac{2}{3}$ into the integrand:

We will engineer a PIVP whose output variable is exactly this integral with the $\tfrac{2}{3}$ baked in, so $S(\infty) = \zeta(3)$ directly — no trailing rational-scaling step. The integrand involves $x^2$, which is unbounded, and $1/(1+e^x)$, which decays exponentially. To turn this into a bounded polynomial system one must simultaneously:

• replace $x^2$ by something that stays bounded while encoding the same growth,
• replace $e^x$ by a polynomial object that stays bounded,
• express the integrand as a polynomial in the state variables.

All three goals are met by introducing five state variables and expressing the integrand via their products.

State Variables

Define

Each of these has a closed ceiling on $[0, \infty)$:

The bounds on $b$ and $c$ come from $\max_t t\,e^{-t} = 1/e$ at $t=1$ and $\max_t t^2\,e^{-t} = 4/e^2$ at $t=2$ respectively.

Polynomial Dynamics

Direct differentiation gives:

For the output variable, notice the algebraic identity

so $\dot{S} = \tfrac{2}{3}\,c\,q$, where the rational factor $\tfrac{2}{3}$ is absorbed into the Ṡ-dynamics so the output variable hits $\zeta(3)$ exactly rather than $\tfrac{3}{2}\zeta(3)$. No post-hoc rational scaling is needed. (We have renamed the integration variable to match the state.)

Collecting everything:

This is a 5-variable polynomial system of degree at most 2 in the right-hand side. All initial conditions are rational numbers, all five state variables are bounded for all $t \geq 0$, and — as a direct consequence of the Fermi-integral evaluation —

The rational coefficient $\tfrac{2}{3}$ in the $\dot{S}$-equation is an inert number, mathematically indistinguishable from an integer coefficient — it simply rescales an accumulator. By folding it into the dynamics we remove the last cosmetic artifact between “the PIVP output” and “the target constant”. They now coincide.

Convergence Rate

The tail estimate

(two integrations by parts) is sharp up to a factor $1/(1+e^{-t}) \in (1/2, 1)$. Numerical integration with DOP853 at $t=50$ confirms $S(50) = 1.20205690316\ldots$ agreeing with $\zeta(3)$ to 15 significant digits. The modulus of convergence is thus

This construction therefore places $\zeta(3)$ directly on the first floor of the analog-time hierarchy. To our knowledge this is the cleanest known 5-variable bounded polynomial witness for $\zeta(3)$ at the real-time level.

Remarks

The construction depends on three non-trivial pieces of structure, each of which deserves notice:

1. The pair $(b, c)$ encodes the powers $t, t^2$ multiplied by $e^{-t}$. Neither $t$ nor $t^2$ is individually bounded, but damping them by $e^{-t}$ makes them bounded while preserving the information needed to reconstruct the integrand. This is the integration-by-parts bookkeeping that converts unbounded polynomial growth into bounded exponentially-damped variables.

2. The variable $q = 1/(1+e^{-t})$ satisfies the autonomous Riccati-like polynomial ODE $\dot{q} = a\,q^2$. It does not require a dedicated logistic variable; the logistic is already encoded by the combination of $a$ and $q$ themselves.

3. The integrand factors polynomially in the state: $x^2/(1+e^x) = c\cdot q$. No division, no special function, no auxiliary variable. This factorization is what makes the whole construction work.

All three features are specific to $s=2$ in the Fermi integral $\int_0^\infty x^s/(1+e^x)\,dx = (1-2^{-s})\Gamma(s+1)\zeta(s+1)$. For general $s$, the integrand $x^s/(1+e^x)$ decomposes as $(s!\cdot c_s)\cdot q$ where $c_s$ is the exponentially-damped polynomial of degree $s$, requiring $s+1$ damped-polynomial variables plus $q, a, S$. So the Fermi route gives a first-floor encoding of $\zeta(s+1)$ for every integer $s \geq 1$ using $s+3$ bounded state variables.

3. The Series Route — Apéry’s Generating Function

Apéry’s original irrationality proof for $\zeta(3)$ (1978) rested on the accelerated series

This series converges very fast — each term shrinks by roughly a factor of $4$ — and it was the crux of van der Poorten’s “A Proof that Euler Missed”. It is natural to ask whether the series itself can be encoded as a PIVP, entirely bypassing the integral representation.

Define the generating function

Then the Apéry identity rewrites as $-\tfrac{5}{2}\,f(-1) = \zeta(3)$, or equivalently $\tfrac{5}{2}\,F(1) = \zeta(3)$ where $F(x) := f(-x) = \sum (-1)^{n-1} x^n / (n^3\binom{2n}{n})$. Evaluating at $x=1$ is the heart of the matter.

Following the same principle we used for the Fermi route, we can absorb the rational prefactor $\tfrac{5}{2}$ directly into the ODE: let $Y(x) := \tfrac{5}{2}F(x)$, so $Y$ satisfies the same left-hand side as (E’) but with right-hand side $\tfrac{5}{2}$ instead of $1$. Then $Y(1) = \zeta(3)$ directly. The discussion below treats $F$ for notational clarity; the scaling to $Y$ is cosmetic and does not affect the analytic obstructions that follow.

3.1 The Recurrence

The coefficients $a_n := 1/(n^3\binom{2n}{n})$ satisfy a one-term ratio:

Rearranging,

This is a pure polynomial (in $n$) recurrence of order 1, with $\mathbb{Q}$-coefficients. One can lift such a recurrence to a linear ODE with polynomial coefficients — this is the translation the D-finite (holonomic) formalism automates.

3.2 The Holonomic ODE

Use the Euler operator $\theta := x\partial_x$, which acts diagonally on monomials: $\theta x^n = n\,x^n$. This converts polynomial-in-$n$ operations on coefficients into polynomial-in-$\theta$ operators on the series.

From (R) we get

(The boundary term $-x$ comes from $a_1 = 1/2$ being absorbed into the lower index when shifting.) Applying (R) term by term:

Expanding the Euler operator using $\theta f = xf'$, $\theta^2 f = x f' + x^2 f''$, $\theta^3 f = xf' + 3x^2 f'' + x^3 f'''$ and collecting:

This is a third-order linear ODE with integer polynomial coefficients and inhomogeneous right-hand side $1$. For $F(x) = f(-x)$ we get the sign-flipped version

The initial data at $x=0$ is rational: differentiating (E) and using $f = \sum a_n x^n$,

Literature pointer: the same ODE (E) appears in Koecher’s 1980 Math. Intelligencer letter, Borwein–Bradley’s Thirty-two Goldbach Variations, and Zagier’s survey of the “15 sporadic” Apéry-like Calabi–Yau pullbacks (Zagier 2009; Almkvist–Zudilin arXiv:math/0402386).

3.3 Singular Structure

The leading coefficient of (E) is $x^2(4-x)$. Singular points are where this vanishes: $\{0, 4\}$, plus $\infty$. Everywhere else the ODE is ordinary — the leading coefficient is nonzero, standard existence/uniqueness applies, and solutions extend analytically.

This point cannot be emphasized enough. The evaluation point $x = -1$ for $f$ (equivalently $x = +1$ for $F$) is an ordinary point of (E), at which $x^2(4-x) = -5 \neq 0$. No Frobenius machinery, no branch cuts, no exotic analytic continuation are needed at the evaluation point. The series $\sum a_n x^n$ extends analytically across $x=\pm 1$ with no drama.

The only substantive singularity — the only genuine obstruction to a naive PIVP encoding — is at $x = 0$, the same point that serves as our initial condition.

3.4 Indicial Analysis at $x = 0$

At a regular singular point, Frobenius theory predicts solutions of the form $x^r \cdot (\text{analytic})$ or $x^r \log(x) \cdot (\text{analytic})$, where the exponents $r$ are roots of the indicial polynomial. For a third-order ODE the indicial polynomial is cubic; substituting the ansatz $y(x) = x^r$ into the leading (lowest-power) part of the equation gives

The roots are

Three pieces of information are encoded here:

1. The root $r=0$ with multiplicity two means the space of analytic solutions at $x=0$ is two-dimensional (as it should be for a third-order equation with one Puiseux solution), but a double root at a regular singular point generally introduces a $\log x$ branch in one of the two solutions. The series $f$ we want is the one without the log branch — pinned by the condition $f(0) = 0$ and the absence of a $\log$ term at first order. This is fine.

2. The root $r = \tfrac{1}{2}$ produces a genuine Puiseux-series solution $y_3(x) = \sqrt{x}\cdot(b_0 + b_1 x + \cdots)$. This is not an analytic function at the origin — it has a $\sqrt{x}$ branch point. In the context of a PIVP, the solution space at $x = 0$ contains elements that cannot be approximated by polynomials in $x$ alone; no finite-dimensional polynomial ODE can track all three independent solutions through the origin.

3. The particular solution (the inhomogeneous one we want, $f(x) = \sum a_n x^n$) is analytic and pinned by $f(0) = 0$. It does not itself have a branch. But the ODE’s solution space is three-dimensional, and any PIVP lift of the equation inherits the full space’s behavior near $x=0$.

3.5 The PIVP Form of (E)

To write (E) as a first-order autonomous polynomial system, introduce state variables $(x, y_1, y_2, y_3) = (x, f, f', f'')$. Then $y_1' = y_2$ and $y_2' = y_3$ by construction, and from (E)

This is not a polynomial RHS — the denominator $x^2(4-x)$ prevents it. To polynomialize, rescale time: let $\tau$ be the new time variable with $d\tau = dt\,/\,[x^2(4-x)]$. Denote $\dot{y} := dy/d\tau$. Then:

This is a 4-variable degree-$\leq 4$ autonomous polynomial system. It is a genuine PIVP. Starting at $\tau = 0$ with $(x, y_1, y_2, y_3) = (x_0, f(x_0), f'(x_0), f''(x_0))$ for any $x_0 \in (0, 4)$, the system integrates until $x$ reaches any target in $(x_0, 4)$. To evaluate $f(1)$ (or $F(1)$, with the sign-flipped (P)), integrate until $x$ reaches 1.

3.6 Numerical Verification

A scipy DOP853 integration of (P) (sign-flipped to F-form) with Frobenius-series initial conditions at $x_0$ gives, for the arrival time $\tau$ satisfying $x(\tau) = 1$:

The system does compute $F(1) = \tfrac{2}{5}\zeta(3)$ to machine precision, for any non-zero starting point. Script: projects/Ripple/experiments/apery_gf_pivp.py.

3.7 Why This Is Not (Yet) a First-Floor Encoding

Three interlinked obstructions.

(a) $x=0$ is a fixed point of $\dot{x} = x^2(4-x)$. From $x_0 = 0$, the $x$-variable stays at zero forever; the integration cannot “start”. We must start from some $x_0 > 0$.

(b) The arrival time diverges as $x_0 \to 0^+$. Computing explicitly,

For target precision $2^{-r}$ on $F(1)$, one needs $x_0 = O(2^{-r})$ (because the Frobenius-series remainder at $x_0$ decays like $x_0^N/4^N$, and the error propagation through the flow is essentially Lipschitz). Thus

This is the same exponential blowup that occurs in any direct encoding of a holonomic series with a regular singular point at the initial condition: the time cost to “escape” the origin grows exponentially in the precision demanded.

(c) Rational initial conditions require approximation. Even if we decide to live with a positive $x_0$ — say $x_0 = 1/4$ — the initial condition $F(1/4) = \sum_{n \geq 1}(-1)^{n-1}/(4^n\,n^3\binom{2n}{n})$ is a transcendental number, not a rational. A PIVP with rational initial conditions can only approximate $F(x_0)$ via a truncated series, introducing a second source of error.

3.8 Attempted Rescue: $x = \tau^2$ Substitution

The Puiseux root $r = \tfrac{1}{2}$ at $x=0$ suggests the uniformization $x = \tau^2$. This doubles the branch cover, potentially converting the $\sqrt{x}$ branch into an integer-power analytic function in $\tau$.

Carrying out the substitution: let $g(\tau) := f(\tau^2)$. Then $g'(\tau) = 2\tau\,f'(\tau^2)$, $g''(\tau) = 2\,f'(\tau^2) + 4\tau^2\,f''(\tau^2)$, and $g'''(\tau) = 12\tau\,f''(\tau^2) + 8\tau^3\,f'''(\tau^2)$. Substituting into (E) and simplifying (a small symbolic computation — done explicitly in experiments/apery_gf_tau2_substitution.py):

Indicial analysis at $\tau = 0$ on the leading part: plugging $g = \tau^r$ and extracting the coefficient of $\tau^{r}$ gives

The indicial roots are

The $\sqrt{x}$ Puiseux branch has been absorbed (it mapped to $r = 1$ in the $\tau$ variable, an integer power). This is real progress: the three solutions of the $\tau$-ODE are all honest Laurent/analytic in $\tau$.

But the double root at $\tau = 0$ is still there, and it still implies a logarithmic branch in one of the solutions. The fundamental solution space now consists of two analytic branches and one $\log\tau$ branch — better than before (no more fractional powers), but still not a pure polynomial-in-$\tau$ system amenable to rational initial conditions at the origin.

Moreover, the coefficient $(4\tau - \tau^3)$ still vanishes at $\tau = 0$, so the PIVP lift continues to have a fixed point at the origin. The $\tau^2$ substitution does not circumvent obstruction (a).

This is the subtle point: the singularity at the origin is not merely a Puiseux issue fixable by uniformization. It is a deeper feature of the operator — the vanishing of the leading coefficient at the same point where the series starts — and no simple algebraic substitution appears to remove it.

4. Comparison and What It Means

The two routes sit in very different positions.

The integral route wins on every operational metric. But the series route remains interesting for several reasons:

1. Generality. The Apéry-generating-function ODE is a special case of a large family — Zagier’s “15 sporadic” Apéry-like sequences, the Calabi–Yau DE’s of Almkvist–Zudilin, and more broadly every holonomic (D-finite) generating function. If we can find a clean way to encode a holonomic GF as a PIVP with rational IC, the method immediately covers hundreds of constants of interest in enumerative combinatorics, mirror symmetry, and the theory of $L$-values.

2. Structural significance. The obstruction at $x=0$ — a regular singular point coinciding with the initial condition — is precisely the situation in which every holonomic generating function with nontrivial recurrence origin data lives. Resolving it (or proving it is unresolvable) would settle a broad question about the analog-computability of holonomic series in bulk, not just $\zeta(3)$.

3. Link to transcendence theory. The monodromy of the Apéry ODE at $\{0, 4, \infty\}$ is the same object studied by Beukers and others in connection with irrationality measures of $\zeta(3)$. The analog-computability question is therefore genuinely connected to classical open problems.

5. Open Questions (for Future Revisits)

Q1 (Zero-IC lift of (E)). Does there exist a polynomial change of coordinates $(x, y_1, y_2, y_3) \mapsto (X_1, \ldots, X_k)$ and a polynomial vector field $P$ on $\mathbb{R}^k$ such that (i) the image of the curve $\tau \mapsto (x(\tau), y_1, y_2, y_3)$ solves $\dot{X} = P(X)$, and (ii) the image of the initial condition $(0, 0, 1/2, 1/24)$ is a rational point of $\mathbb{R}^k$? This is the precise statement of “encode (E) as a rational-IC PIVP”. If yes, $\zeta(3)$ is first-floor via the series route, not only the integral route.

Q2 (First-floor series encoding of general holonomic GFs). Which holonomic generating functions admit rational-IC PIVP encodings of their evaluation points? Is there a structural characterization, e.g., “those whose ODE has no indicial Puiseux roots and no Fuchsian log term at $x=0$”?

Q3 (Second-floor lower bound). Is the $\Theta(2^r)$ modulus for the series route (as computed in §3.7) optimal? That is, does every PIVP encoding of $\zeta(3)$ via its generating function necessarily require exponential time to reach $2^{-r}$ precision, or is there a clever reparameterization that improves the modulus?

Q4 (Relationship to classical obstructions). The Apéry recurrence and its associated ODE are also central objects in the irrationality measure theory of $\zeta(3)$. Is there a dictionary between the analog-computability obstructions uncovered here and the Diophantine obstructions (irrationality exponent, transcendence degree) studied in classical number theory?

6. What We Actually Learned

The Fermi-integral route shows that $\zeta(3)$ is real-time analog-computable — it fits on the first floor of the analog complexity hierarchy via the explicit 5-variable bounded PIVP written in §2.

The series route shows that encoding a general holonomic generating function as a PIVP is subtle, and the subtlety sits at the origin where the Frobenius theory introduces branch cuts and log terms that do not live in the category of polynomial vector fields. This is not a failure; it is a genuine mathematical obstruction that future work should engage with head-on.

Two encodings. One works cleanly, one meets a wall. Together they map out a small but real piece of the boundary between “real-time analog-computable” and everything else.

格物致知，誠意正心。

To investigate things is to extend knowledge; to make one’s intentions sincere is to rectify one’s heart.

— 《禮記·大學》 · Book of Rites, The Great Learning