Researchers have developed a semi-analytic and semi-numerical formulation for the second post-Newtonian (2PN) equations of motion for point-mass N-body systems in harmonic gauge. This advance is crucial for describing with greater accuracy the dynamics of complex gravitational systems, where relativistic effects are significant but not so extreme as to require a purely relativistic description. The formulation employs Hadamard regularization to handle the singularities inherent in point-mass interactions.

The resulting equations of motion are separated into a closed analytic contribution and a non-closed integral contribution. The integral part, which contains the main complexity, has been analyzed for its singular structure and further regularized into a numerically evaluable representation. This hybrid approach allows combining the precision of analytical solutions with the flexibility of numerical methods, facilitating application to real astrophysical systems.

The formulation was applied to two benchmark systems: Sun-Jupiter-Saturn and Sun-Mercury-Venus. The instantaneous non-closed 2PN acceleration along Newtonian trajectories and its leading finite-time relative-distance response were evaluated. In both benchmarks, the non-closed acceleration remained a small fraction of the complete 2PN acceleration. The induced relative-distance perturbation remained oscillatory, with oscillation amplitudes that could reach larger values at later times. This work lays the groundwork for more precise calculations in the dynamics of stellar and planetary systems.