Closed-form solutions of spinning, eccentric binary black holes at 1.5 post-Newtonian order

The closed-form solution of the 1.5 post-Newtonian (PN) accurate binary black hole (BBH) Hamiltonian system has proven to be evasive for a long time since the introduction of the system in 1966. Solutions of the PN BBH systems with arbitrary parameters (masses, spins, eccentricity) are required for modeling the gravitational waves emitted by them. Accurate models of gravitational waves are crucial for their detection by LIGO/Virgo and LISA. Only recently, two solution methods for solving the BBH dynamics were proposed in Ref. [G. Cho and H. M. Lee, Phys. Rev. D 100, 044046 (2019)] (without using action-angle variables), and (action-angle based). This paper combines the ideas laid out in the above articles, fills the missing gaps and compiles the two solutions which are fully 1.5PN accurate. We also present a public Mathematica package BBHPNTOOLKIT which implements these two solutions and compares them with the result of numerical integration of the evolution equations. The level of agreement between these solutions provides a numerical verification for all the five action variables constructed in Refs. [S. Tanay et al., Phys. Rev. D 103, 064066 (2021), S. Tanay et al., Phys. Rev. D 107, 103040 (2023)]. This paper hence serves as a stepping stone for pushing the action-angle-based solution to 2PN order via canonical perturbation theory.