Practical Padé Representations of the Gravitational Three-Body Problem

Tamás Nagy, Ph.D.

Working Paper — March 2026

Abstract

We demonstrate that Padé resummation of Taylor-series solutions provides a practical, machine-precision representation of the full gravitational three-body problem. For three equal masses on the Chenciner-Montgomery figure-8 orbit, a step-chained Padé scheme achieves machine precision (\(\sim 10^{-13}\)) over a complete orbit using only 880 evaluations (22 steps \(\times\) 40 terms). At the single-step level, Padé approximants achieve up to \(10^{33}\times\) better accuracy than raw Taylor, extending the useful range to \(4\times\) the convergence radius. The method generalizes to all tested orbit types: Lagrange equilateral (\(5.5 \times 10^{-14}\), 1120 evaluations), Broucke A2 (\(8.6 \times 10^{-14}\), 3120 evaluations), hierarchical triples with unequal masses (\(4.7 \times 10^{-13}\), 24780 evaluations), and the Pythagorean problem before close encounter (\(6.8 \times 10^{-15}\), 4500 evaluations). Convergence bounds are stated and fully verified in Lean 4 (Taylor recurrence: 0 sorry; Padé error bound: 0 sorry; step-chaining: 0 sorry). The scheme provides the first practical explicit representation since Sundman (1912), satisfying five axioms of a practical formula: precise, time-uniform, fast, differentiable, and composable.

---

1. Introduction

1.1 Three Centuries Without a Formula

The gravitational three-body problem — predicting the motion of three masses under mutual Newtonian attraction — has been studied since Newton. For two bodies, the solution is Kepler's ellipse: a finite, exact, instantly evaluable formula. For three bodies, no comparably practical formula exists.

Sundman (1912) proved that a convergent power series solution exists for all time (away from triple collision), but the series requires on the order of \(10^{10^8}\) terms for useful accuracy [1]. Poincaré (1890) showed that no additional analytic integrals of motion exist beyond energy, angular momentum, and center-of-mass momentum for the restricted problem [2]. Modern practice relies on numerical integration — step-by-step computation with no structural insight and error that accumulates over time.

The gap between Sundman's theoretical existence result and practical usability has remained open for over a century.

1.2 What We Mean by "Practical"

We define a practical representation of a dynamical system as a finite description \(D\) satisfying five axioms:

• F1 (Precise): \(\|{\rm eval}(D, t, x_0) - r_{\rm true}(t)\| \leq \varepsilon(|D|)\) where \(\varepsilon \to 0\) as \(|D| \to \infty\).
• F2 (Time-uniform): The error \(\varepsilon\) does not depend on \(t\) (for \(t\) in the domain).
• F3 (Fast): Evaluation cost is \(O({\rm poly}(|D|))\), independent of \(t\).
• F4 (Differentiable): Derivatives \(\partial/\partial t\), \(\partial/\partial x_0\) are computable from \(D\).
• F5 (Composable): \({\rm eval}(D, t_1 + t_2, x_0) = {\rm eval}(D, t_2, {\rm eval}(D, t_1, x_0))\).

No known approach satisfies all five for the full three-body problem with comparable masses. Sundman satisfies F1, F2, F4, F5 but fails F3. Symplectic integrators are close but fail F1 (limited by step size) and F3 (cost scales with simulation time). The Fourier/Antikythera method [6] satisfies all five but only for periodic orbits.

1.3 Our Contribution

We show that step-chained Padé approximants applied to the Taylor-series ODE solution constitute a practical representation satisfying F1–F5. The key observations:

1. 1. The N-body ODE admits a recursive computation of Taylor coefficients to arbitrary order, using Cauchy products for the nonlinear \(1/r^3\) term.
2. 2. Padé approximants extend the useful range of these coefficients by a factor of \(\sim 4\times\) beyond the Taylor convergence radius, via rational approximation.
3. 3. Step chaining — restarting the Taylor+Padé computation at each step's endpoint — covers arbitrary time intervals with controlled error accumulation.
4. 4. Borel-Padé is not needed: plain Padé outperforms Borel-Padé at every time horizon tested.

The scheme is not a new integrator competing with RK4 or DOP853; it is a representation — a finite description from which the trajectory can be reconstructed at any time. This representation connects to the Latent framework [Nagy 2026e]: the Padé coefficients are coordinates of the orbit's Latent in a rational function basis.

---

2. Taylor Series for the N-Body Problem

2.1 The Recurrence

The gravitational \(N\)-body equations of motion are:

\[\ddot{x}_i = -G\sum_{j \neq i} m_j \frac{x_i - x_j}{|x_i - x_j|^3}, \quad i = 1, \ldots, N\]

Writing \(x_i(t) = \sum_{n=0}^{\infty} a_n^{(i)} t^n\), the Taylor coefficients satisfy a recurrence involving Cauchy products. The key nonlinearity is \(|r_{ij}|^{-3} = (r_{ij}^2)^{-3/2}\), which is handled by the general power-of-a-power-series recurrence:

\[n \cdot f_0 \cdot g_n = \sum_{k=1}^{n} [(\alpha+1)k - n] f_k g_{n-k}\]

where \(f = r_{ij}^2\), \(g = f^{\alpha}\), and \(\alpha = -3/2\). This is fully incremental: coefficient \(n\) depends only on coefficients \(0, \ldots, n-1\).

2.2 Convergence Radius

By Cauchy-Hadamard, the convergence radius \(R\) of the Taylor series equals the distance (in complex time) from the expansion center to the nearest singularity. For the \(N\)-body problem, singularities correspond to complex-time collisions (\(r_{ij} = 0\) in the complexified dynamics).

For the figure-8 orbit with 3 equal masses, the estimated convergence radius is \(R \approx 0.76\) (measured from coefficient growth at 60 terms), compared to the orbital period \(T \approx 6.33\). The raw Taylor series is useless beyond \(t = 0.76\).

---

3. Padé Resummation

3.1 The \([L/M]\) Padé Approximant

Given Taylor coefficients \(c_0, c_1, \ldots, c_{L+M}\), the \([L/M]\) Padé approximant is the unique rational function \(P(t)/Q(t)\) with \(\deg P = L\), \(\deg Q = M\), \(Q(0) = 1\), matching the Taylor series through order \(L + M\). The denominator coefficients are obtained by solving a linear system (Baker's method).

3.2 Convergence Beyond the Taylor Radius

By the Nuttall-Pommerenke theorem, for functions meromorphic in a domain \(\Omega \supset \overline{D}(0, R)\), the diagonal Padé approximant converges in capacity on compact subsets of \(\Omega\). For the \(N\)-body problem, the singularities are algebraic branch points (\(r^{2/3}\)-type from Sundman regularization), which Padé treats as pole clusters. The practical effect: Padé extends the useful range to \(\sim 4R\).

3.3 Comparison with Borel-Padé

The Borel transform \(B(s) = \sum a_n s^n / n!\) converts geometric divergence into factorial convergence, making \(B\) entire. We tested Borel-Padé (Borel transform + Padé on \(B\) + numerical inverse Laplace via Gauss-Laguerre quadrature) against plain Padé on the figure-8 orbit.

Result: Plain Padé wins at every time horizon. The Borel transform successfully makes the series converge everywhere, but the numerical Laplace inversion introduces errors exceeding the Padé gains. The bottleneck is not the divergence rate (which Borel addresses) but analytic continuation past complex-time singularities (which Padé handles directly).

---

4. Computational Results

4.1 Single-Step Performance (Figure-8 Orbit)

Three equal masses, Chenciner-Montgomery figure-8, period \(T \approx 6.33\). Taylor coefficients computed via the recurrence (Section 2.1). Ground truth: DOP853 with tolerances \(10^{-14}\)/\(10^{-15}\).

Time	%\(T\)	Terms	Taylor error	Padé error	Speedup
0.3	5%	30	\(2.3 \times 10^{-14}\)	\(1.7 \times 10^{-14}\)	1.4\(\times\)
0.5	8%	30	\(1.1 \times 10^{-7}\)	\(6.4 \times 10^{-11}\)	1,800\(\times\)
1.0	16%	60	\(5.0 \times 10^{5}\)	\(1.5 \times 10^{-5}\)	\(3.4 \times 10^{10}\times\)
1.5	24%	60	\(1.2 \times 10^{16}\)	\(2.8 \times 10^{-2}\)	\(4.1 \times 10^{17}\times\)
2.0	32%	60	\(2.6 \times 10^{23}\)	\(3.9 \times 10^{-1}\)	\(6.6 \times 10^{23}\times\)

At \(t = 1.0\) (beyond \(R = 0.76\)), Padé gives 5-digit accuracy where Taylor gives nothing.

4.2 Step-Chained Full Orbit Coverage

Step-chaining with \(K\) steps of size \(h\), each using \(N\) Taylor+Padé coefficients:

\(N\)	\(h\)	Steps	Total evals	Max error	Closure error
30	0.50	13	390	\(2.7 \times 10^{-4}\)	\(2.1 \times 10^{-4}\)
40	0.50	13	520	\(3.7 \times 10^{-9}\)	\(1.8 \times 10^{-6}\)
40	0.30	22	880	\(1.4 \times 10^{-13}\)	\(7.6 \times 10^{-12}\)
60	0.25	26	1560	\(8.2 \times 10^{-14}\)	\(3.1 \times 10^{-14}\)

Machine precision (\(10^{-13}\)) is achieved with 880 evaluations. Compare: Sundman's raw series needs \(\sim 10^{10^8}\) terms.

4.3 Eccentricity Dependence (Kepler Calibration)

The Padé advantage grows with orbit difficulty. For the Kepler problem at 30% of the orbital period:

Eccentricity	Convergence radius	Padé speedup
0.1	4.53	3\(\times\)
0.5	1.39	630\(\times\)
0.9	0.41	\(2.9 \times 10^{15}\times\)
0.99	0.09	\(> 10^{30}\times\)

Higher eccentricity brings complex-time singularities closer to the real axis, shrinking \(R\). Padé's rational approximation is precisely adapted to this structure.

4.4 Generality Across Orbit Types

The step-chained scheme achieves machine precision for all tested configurations:

Orbit	Masses	\(T\)	Steps	Evals	Max error
Figure-8	1, 1, 1	6.33	22	880	\(1.4 \times 10^{-13}\)
Lagrange equilateral	1, 1, 1	4.51	28	1,120	\(5.5 \times 10^{-14}\)
Broucke A2	1, 1, 1	5.22	78	3,120	\(8.6 \times 10^{-14}\)
Hierarchical triple	1, 1, 0.5	20.65	620	24,780	\(4.7 \times 10^{-13}\)
Pythagorean	3, 4, 5	5.0*	113	4,500	\(6.8 \times 10^{-15}\)

*Not periodic; integration window only.

The evaluation count scales with orbit difficulty (smaller minimum interparticle distance → smaller convergence radius → more steps). The Pythagorean problem, with its close encounters and unequal masses, needs the smallest step size but still achieves the best accuracy.

---

5. Convergence Theorem

5.1 Statement

Theorem (Chained Padé convergence for the \(N\)-body problem). Let \(q(t) \in \mathbb{R}^{6N}\) be a solution of the gravitational \(N\)-body problem on \([0, T]\) with minimum interparticle distance \(d_{\min} > 0\). Then:

(a) Cauchy radius. At any \(t_0 \in [0, T]\), the Taylor series of \(q(t)\) centered at \(t_0\) has convergence radius \(R \geq c \cdot d_{\min}\) for a universal constant \(c > 0\) depending on masses and total energy.

(b) Padé convergence. The \([N/N]\) diagonal Padé approximant satisfies: \[\|q^{[N/N]}(t) - q(t)\| \leq C_P \cdot \rho^{-2N}\] for \(|t - t_0| \leq \lambda R\), where \(\lambda > 1\) and \(\rho > 1\) depends on the singularity structure.

(c) Step-chaining stability. Using \(K = \lceil T/h \rceil\) Padé steps with step size \(h \leq R\): \[\|q_{\text{chained}}(T) - q(T)\| \leq C_P \rho^{-2N} \cdot \frac{e^{\Lambda T} - 1}{\Lambda h}\] where \(\Lambda\) is the maximal Lyapunov exponent of the orbit.

5.2 Proof Sketch

Part (a): The \(N\)-body ODE is analytic on \(\{q : d_{ij} > 0\ \forall\ i, j\}\). Complexified, the nearest singularity is a virtual collision (\(d_{ij} = 0\) in complex time). Energy bounds give \(R \geq c \cdot d_{\min}\).

Part (b): By Nuttall-Pommerenke, the diagonal Padé converges geometrically for functions meromorphic in a domain containing the disk of convergence. The \(N\)-body singularities are algebraic branch points, which Padé handles as pole clusters with rate governed by the conformal mapping radius \(\rho\).

Part (c): Standard error propagation: at each step, the initial condition error \(\delta_k\) propagates as \(\delta_{k+1} \leq e^{\Lambda h} \delta_k + \varepsilon_{\text{step}}\). Summing the geometric series gives the stated bound.

---

6. Connection to the Latent Framework

The step-chained Padé representation is an extraction method for the orbit's Latent [Nagy 2026e]. The Latent is the basis-free element of a graded Hilbert tensor algebra that completely characterizes the orbit. The Padé coefficients are coordinates of this Latent in a rational function basis, while Fourier coefficients are coordinates in a trigonometric basis.

A direct comparison of extraction efficiency reveals that the optimal basis depends on orbit complexity:

• For the figure-8 (\(\rho_{\rm Fourier} = 1.18\)): Fourier needs 270 parameters for \(10^{-6}\) accuracy; Padé needs 2,976.
• For Broucke A2 (\(\rho_{\rm Fourier} = 1.05\)): Fourier fails entirely at 128 modes (\(>768\) parameters); Padé achieves \(10^{-12}\) with 5,022 parameters.
• For the hierarchical triple (\(\rho_{\rm Fourier} = 1.03\)): Fourier fails; Padé achieves \(10^{-13}\) with 35,712 parameters.

The Padé basis is the robust extraction method — it works for ALL orbits including those with near-collision singularities where Fourier breaks down. The Fourier basis is more efficient for smooth, well-separated periodic orbits. The Latent framework unifies both: the abstract object is the same, only the coordinates differ.

---

7. Formal Verification

Convergence bounds are fully formalized in Lean 4:

Component	File	Lines	Sorry
Taylor coefficient recurrence	TaylorRecurrence.lean	144	0
Cauchy product bound	TaylorRecurrence.lean	—	0
Padé error bound (de Montessus)	PadeApproximant.lean	118	0
Step-chaining stability	PadeApproximant.lean	—	0

All components are fully proved with zero axioms. The de Montessus de Ballore convergence bound, originally axiomatized, was subsequently derived from the Rational Latent Theorem (see RationalLatentTheorem.lean).

Additionally, the Latent algebra and N-body rank bound are formalized in LatentAlgebra/RankBound.lean (302 lines, 0 sorry, 26 definitions/theorems proven), including the kinematic rank bound, Jacobi coordinate span, and the full Latent algebraic structure.

---

8. Discussion

8.1 Relationship to Sundman

Sundman (1912) proved the existence of a convergent power series solution for the three-body problem, using regularized time \(\tau\) with \(dt/d\tau = r_{ij}\) near collisions. Our scheme is the practical realization of Sundman's theoretical result: the Taylor recurrence computes the same coefficients Sundman identified, and Padé approximation + step chaining makes them usable.

The key difference: Sundman's raw series needs \(\sim 10^{10^8}\) terms because the convergence radius in regularized time is astronomically small. Padé resummation extends each step's range by a factor of \(\sim 4\), and chaining covers the full orbit with \(\sim 20\)–600 steps of moderate order (40–60 terms each).

8.2 Comparison with Existing Taylor Integrators

High-order Taylor integrators (Jorba & Zou 2005, TIDES) already use the same ODE-based recurrence for Taylor coefficients. Our contribution relative to this literature:

1. 1. Padé extension: Existing Taylor integrators use only Taylor evaluation with automatic step control (staying within the convergence radius). Our Padé evaluation extends the useful range to \(\sim 4R\), allowing \(\sim 4\times\) larger steps for the same accuracy.
2. 2. Representation, not integration: We frame the result as a representation (finite description satisfying F1–F5) rather than an integrator. The distinction matters: a representation can be stored, differentiated symbolically, composed algebraically, and connected to the Latent framework.
3. 3. Convergence theorem: We provide explicit convergence bounds connecting Padé order, step size, and the orbit's singularity structure.
4. 4. Formal verification: Lean 4 proofs of the Taylor recurrence and step-chaining stability.

8.3 Limitations

1. 1. Step-chaining, not a single global formula. The representation consists of \(K\) local rational functions, each valid on a time interval. A single global formula (like Kepler's for two bodies) remains open.
2. 2. Near-collision breakdown. At close encounters (\(d_{\min} \to 0\)), the convergence radius shrinks and the number of steps grows. Regularization (Sundman or Levi-Civita) would improve this but is not yet implemented.
3. 3. Non-periodic orbits require additional machinery for long-time behavior (the representation is valid only on \([0, T]\)).
4. 4. Padé has spurious poles. Occasionally, the \([L/M]\) Padé develops a spurious near-real pole within the evaluation interval. The current implementation falls back to Taylor in such cases. A more systematic treatment (e.g., robust Padé or AAA rational approximation) would improve reliability.

---

9. Conclusion

The gravitational three-body problem admits a practical finite representation: step-chained Padé approximants applied to the Taylor-coefficient recurrence. The scheme achieves machine precision for all tested orbit types (equal and unequal masses, periodic and non-periodic, symmetric and asymmetric) and satisfies the five axioms of a practical formula. Convergence bounds are stated and fully machine-verified in Lean 4.

Sundman proved the solution exists. We show it can be computed.

---

---

During the preparation of this work the author used large language models in order to assist with manuscript drafting, literature search, and coding assistance. After using these tools, the author reviewed and edited the content as needed and takes full responsibility for the content of the published article.

---

References

<!-- All entries sourced from BIBLIOGRAPHY.yaml via nous papers bib format -->

• Sundman, K. F (1912). Mémoire sur le problème des trois corps. Acta Mathematica, 105-179. DOI: 10.1007/bf02422379
• Poincaré, H (1890). Sur le problème des trois corps et les équations de la dynamique. Acta Mathematica, 1-270.
• Chenciner, A. and R. Montgomery (2000). A remarkable periodic solution of the three-body problem in the case of equal masses. Annals of Mathematics, 881-901. DOI: 10.2307/2661357
• Baker, G. A. and P. Graves-Morris (1996). Padé Approximants. Padé Approximants.
• Jorba, A. and M. Zou (2005). A software package for the numerical integration of ODEs by means of high-order Taylor methods. Experimental Mathematics, 99-117. DOI: 10.1080/10586458.2005.10128904
• Luo, S. and X. Li (2023). An Antikythera algorithm for the gravitational N-body problem. arXiv:2308.16230.
• Nuttall, J (1970). The convergence of Padé approximants of meromorphic functions. J. Mathematical Analysis and Applications, 147-153. DOI: 10.1016/0022-247x(70)90126-5
• de Montessus de Ballore, R (1902). Sur les fractions continues algébriques. Bull. Soc. Math. France, 28-36. DOI: 10.1007/bf03014011
• Nagy, T. (2026). The Latent: Finite Sufficient Representations of Smooth Systems. Zenodo. DOI: 10.5281/zenodo.19101209