The N-Body 2PN Hamiltonian and Numerical Integration of the Equations of Motion

Technical Summary: The N-Body 2PN Hamiltonian and Numerical Integration of the Equations of Motion

Problem Statement
The dynamics of $N$ gravitating bodies are a fundamental problem in astrophysics, yet exact analytic solutions exist only for the Newtonian two-body case. While Post-Newtonian (PN) expansions have successfully modeled relativistic corrections for compact binaries up to high orders (5PN/6PN), extending these formalisms to general $N$-body systems ($N > 2$) remains challenging. Specifically, the complete second-order Post-Newtonian (2PN) Hamiltonian has been known in closed analytic form only for systems of up to three point masses. The primary obstacle for $N \ge 4$ is the analytical complexity of the four-point correlation function, $U^{TT}_{(4)}$, which arises from the transverse-traceless static potential term $U^{TT}$ in the 2PN Hamiltonian. Previous attempts to derive this term resulted in expressions containing integrals that lacked known closed-form analytic solutions, effectively halting the analytic description of 2PN dynamics for systems with four or more bodies.

Methodology
The authors present a hybrid analytical-numerical approach to resolve this impasse:

1. Analytic Derivation: The paper derives an explicit analytic expression for the general $N$-body 2PN Hamiltonian in the ADM gauge. The Hamiltonian is decomposed into Newtonian ($H_N$), 1PN ($H_{1PN}$), and 2PN ($H_{2PN}$) parts. The 2PN part is further broken down into two-point, three-point, and four-point correlation functions. While the two- and three-point functions have known closed forms, the four-point function $U^{TT}_{(4)}$ is expressed as a sum of analytic terms and a single remaining integral term, denoted $I^{\ln}_{ab;cd}$.
2. Numerical Evaluation: Recognizing that $I^{\ln}_{ab;cd}$ (Equation 17) has no known closed algebraic form, the authors employ the open-source Cuba library, specifically the deterministic Cuhre routine, to evaluate this 3D integral numerically. The integrand involves infinite domains and singularities, which are handled via coordinate transformation to a unit cube and adaptive subdivision. The authors utilize complex-step differentiation to compute the gradients of this integral with respect to particle positions, ensuring machine-precision accuracy for the forces without introducing significant finite-difference errors.
3. Integration Schemes: To solve the equations of motion, the authors propose efficient numerical integration strategies. They utilize a multiple time-stepping approach (Strang splitting or impulse methods) to separate the computationally expensive $U^{TT}_{(4)}$ term from the dominant parts of the Hamiltonian ($H_0$). This allows the use of adaptive high-order methods for close encounters and symplectic integrators for long-term stability, minimizing the number of expensive integral evaluations required per timestep.

Key Contributions

• General N-Body 2PN Hamiltonian: The paper provides the first analytic expression for the 2PN Hamiltonian of an arbitrary number of point masses in the ADM gauge, valid up to a single integral term.
• Numerical Solution of $I^{\ln}_{ab;cd}$: The authors demonstrate that the previously intractable integral $I^{\ln}_{ab;cd}$ can be evaluated to machine precision. They validate this by cross-checking numerical results against known analytic solutions for related integrals and verifying the convergence behavior.
• Force and Energy Computation: The work details efficient algorithms for computing the gradients of the integral term, reducing the computational cost through symmetry exploitation (reducing evaluations by a factor of 4) and parallelization. The computational complexity scales as $O(N^4)$, which is manageable for few-body systems.
• Validation of Three-Body Results: The numerical framework was used to verify the correctness of the existing closed-form expression for the three-point correlation function $U^{TT}_{(3)}$ derived in prior literature.

Results
The authors tested their methods on two specific four-body systems:

1. Close Binary-Binary Encounter: A chaotic resonance interaction involving four equal-mass bodies. The simulations showed that while the trajectories of the full 2PN system (including $U^{TT}_{(4)}$) and the system without $U^{TT}_{(4)}$ diverge significantly during strong, close multi-body encounters, the contribution of the four-point term to the total energy and forces is generally small (on the order of a few percent) except when all bodies are in close proximity. The system exhibited chaotic behavior where tiny perturbations from the four-point term led to vastly different long-term outcomes.
2. Hierarchical Binary-Binary System: A system of two binaries orbiting each other at a larger distance. In this configuration, the inclusion of $U^{TT}_{(4)}$ resulted primarily in small phase shifts within the individual binaries, with negligible effects on the center-of-mass motion or the orientation of the binaries. The relative contribution of $U^{TT}_{(4)}$ to the total Hamiltonian was found to be approximately $5 \times 10^{-3}$.

In both cases, the conservation of energy, linear momentum, and angular momentum was monitored. The authors observed that while energy errors remained bounded without secular drift, linear and angular momentum exhibited a diffusive, random-walk-like accumulation of error. This is attributed to the fact that the numerically evaluated gradients of $U^{TT}_{(4)}$ do not satisfy translation and rotation invariance to machine precision, introducing small spurious perturbations at each integration step.

Significance and Claims
The paper claims to bridge the gap between analytical treatments and numerical approaches for $N$-body dynamics at the 2PN order. By providing a method to compute the full 2PN Hamiltonian for $N \ge 4$ bodies, the authors enable self-consistent simulations of systems such as binary-binary encounters, hierarchical triples/quadruples, and dense star cluster cores.

The authors modestly state that while the four-point correlation function is essential for a complete analytic description, its practical impact on the dynamics is often small, becoming significant primarily during close multi-body encounters. Consequently, the paper suggests that for many applications, the computational cost of evaluating $U^{TT}_{(4)}$ might be managed by switching it off or adjusting error tolerances based on the average inter-particle distance. The work facilitates systematic studies of 2PN effects in $N > 3$ systems and allows for direct comparisons with numerical relativity predictions, though the authors note that future work is required to find analytic approximations for $I^{\ln}_{ab;cd}$ or alternative numerical strategies to further improve efficiency.