  Accurate and efficient modeling of the dynamics of binary black holes (BBHs)
is crucial to their detection and parameter estimation through gravitational
waves, both with LIGO/Virgo and LISA. General BBH configurations will have
misaligned spins and eccentric orbits, eccentricity being particularly relevant
at early times. Modeling these systems is both analytically and numerically
challenging. Even though the 1.5 post-Newtonian (PN) order is Liouville
integrable, numerical work has demonstrated chaos at 2PN order, which impedes
the existence of an analytic solution. In this article we revisit integrability
at both 1.5PN and 2PN orders. At 1.5PN, we construct four (out of five) action
integrals. At 2PN, we show that the system is indeed integrable - but in a
perturbative sense - by explicitly constructing five mutually-commuting
constants of motion. Because of the KAM theorem, this is consistent with the
past numerical demonstration of chaos. Our method extends to higher PN orders,
opening the door for a fully analytical solution to the generic eccentric,
spinning BBH problem.
