High precision predictions are needed in the upcoming LVK runs, in space-based detectors such as LISA [1]}, and in future ground-based detectors such as LIGO-India [2]}, Cosmic Explorer[3]} and Einstein Telescope[4]}. Precise studies of GWs waveforms are very cumbersome due to the multiple physical scales non-linearly coupled through General Relativity, and in order to obtain precision predictions for GWs observables one often finds approximate solutions of Einstein's field equations, expanded around the Newtonian results. The coalescence of a compact binary system can be divided in three distinctive phases: inspiral, merger, and ring-down, and GWs models are a mix of analytical techniques (see [5]}, [6]}, [7]}, [8]} for reviews) used for the inspiral phase, of numerical-relativity methods used for the merger and the ring-down phases, which can be resummed within the effective-one-body approach [9]}.
During the inspiral phase, the two binary constituents are very far away one from each other, they are slowly moving and weakly interacting, and they can be studied within the so-called Post-Newtonian (PN) perturbative scheme.
The PN approximation is a low-velocity weak-field expansion, that allow us to study the dynamics of general-relativistic bound systems by adding perturbative correction in the dimensionless parameter \(v/c\) , where \(v\)  is the typical velocity of the system. Since for self-gravitating binary systems the virial theorem holds,
the potential energy of the binary is expected to be of the same order of magnitude as the kinetic energy, that is:
\(v^2\approx \frac{G_N m }{r} \ ,\) 
