In particular, the two milestones that have been adapted to the study of the gravitational two-body problem are Effective Field Theories[1]} (EFT)(see [2]}, [3]}, [4]} for introductions) and modern Scattering Amplitudes Multi-loop techniques.
Effective Field Theories allow us to study a physical system at a given length scale by focusing only on the minimal set of degrees of freedom, which are relevant to describe a specific problem.
They have been first applied in the case of a coalescing binary system in [5]}, and later systematized in [6]}, and in [7]} in the case of binaries with spin (see [8]}, [9]}, [10]}, [11]}, [12]}, [13]} for modern reviews). They allow us to study the dynamics of a binary coalescing system using a tower of effective theories, each one focused on a specific length scale, and to compute Post-Newtonian corrections to General Relativity with a Feynman diagrammatic approach, characterized by a clear power counting in the expansion parameter \(v/c\)  [14]}, [15]}, [16]}, [17]}, and connections between physics at separate scales are made with renormalization group techniques. Thanks to the so-called method of regions [18]}, [19]}, [20]}, [21]}, one can organize the computation of PN corrections into two different EFTs, by splitting the carriers of the gravitational force, the gravitons, into potential gravitons, which are the short distance components that mediate the gravitational binding between the binary components, and radiation gravitons, which are the long-range modes emitted by the system. In the so-called near zone the system can be described by two point-like particles exchanging potential gravitons in a flat spacetime, whereas in the far zone we can detail the problem with a single source expanded in multipole moments and emitting radiation gravitons [22]}, [23]}, [24]}. The contributions from the two EFTs have to be summed in order to get consistently PN corrections.
The second milestone is the use of Quantum Field Theory techniques for the evaluations of Feynman diagrams. PN contributions, which are computationally organized into near and far zone contributions, can be topologically mapped into massless multi-loop two points and one point functions, respectively, with a relation that has been first introduced in [25]} in the case of potential gravitons. Exploiting this relation, one can perform computations using the machinery of modern amplitudes multi-loop techniques. Instead of directly evaluating all the Feynman diagrams appearing in the theory, we can divide the computations in three different stages: first one can make a tensor decomposition of an amplitude containing different free tensorial indices, by writing it as a linear combination of Lorentz vectors and tensors, multiplied by coefficients called form factors, where all indices are contracted. Then, one can use some special identities of dimensional regularization relating different integrals, known as Integration by parts identities (IBPs) [26]}, [27]}, Lorentz invariance identities, and other symmetries, to decompose all the scalar integrals appearing in the form factors in a minimal set of integrals, known as Master Integrals (MIs), which form some sort of basis of the vector space of integrals. Eventually, one has to actually evaluate only the MIs, a procedure that can be done with direct evaluation or with the use of difference [28]}, [29]} or differential equations [30]}, [31]} (for reviews see [32]}, [33]}).
These new methods allowed rapidly to obtain higher order Post-Newtonian corrections in this new approach: at 2PN in [34]}, at 3PN in [35]}, at 4PN in [36]}, [25]}, and at 5PN in [38]}, [39]}, [40]}, and up to NNNLO including Spin-Orbit interactions [41]}, [42]}. An alternative derivation up to 6PN can be found in [43]}, [44]}, [45]}, [46]}, [47]}, [48]}. Starting at 4PN order, there are some radiative terms contributing to the conservative dynamics appearing in the far zone known as Hereditary effects[49]}, [50]}, [51]}, [52]}, [53]}, [54]}, [8]}, [56]}, [57]}, [47]}, [59]}, [60]}, [61]}, [62]}, [63]}, [64]}, [65]}, [66]}, as first noticed in [49]}. These terms are due to GWs emitted by the system and then scattered-back into the same GWs source. At the moment there are different derivations of these terms at 5PN, made by different communities and using different techniques [68]}, [57]}, [47]}, [71]}, which do not give the same predictions for physical observables, such as the scattering angle. Within the EFT approach that we have so far described, radiative terms have to be studied using the Schwinger-Keldysh (or In-In or CTP) formalism [72]}([73]}, [74]} for reviews), which was first developed by Schwinger in [75]}, to compute expectation values in quantum mechanics from a path integral approach, and then adapted to NRGR to get the radiative terms at 2.5PN in [52]}, at 3.5PN in [53]}, at 4PN in [54]}, [8]}, [56]}, at 5PN in [57]}, [47]}, [59]}, [60]}, [61]}, and partially up to 7PN order in [62]}. These terms are not fully understood at the moment, and this puzzle represents one of the main obstacle for the computation of higher order PN corrections.
