Fundamental Gravity and Gravitational Waves

While being as old as general relativity itself, the gravitational two-body problem has never been under so intense investigation as it is today, spurred by both phenomenological and theoretical motivations. The observations of gravitational waves emitted by compact binary coalescences bear the imprint of the source dynamics, and as the sensitivity of detectors improve over years, more accurate modeling is being required. The analytic modeling of classical gravitational dynamics has been enriched in this century by powerful methods borrowed from field theory. Despite being originally developed in the context of fundamental particle quantum scatterings, their applications to classical, bound system problems have shown that many features usually associated with quantum field theory, such as, e.g., divergences and counterterms, renormalization group, loop expansion, and Feynman diagrams, have only to do with field theory, be it quantum or classical. The aim of this work is to present an overview of this approach, which models massive astrophysical objects as nonrelativistic particles and their gravitational interactions via classical field theory, being well aware that while the introductory material in the present article is meant to represent a solid background for newcomers in the field, the results reviewed here will soon become obsolete, as this field is undergoing rapid development.


Introduction
The recent detections of gravitational waves (GWs) [1][2][3] by the large interferometric detectors LIGO [4] and Virgo [5] have revived the scientific interest on the gravitational two-body problem. All events detected so far have been originated by coalescences of compact binary systems, which are best searched for by correlating noisy observational data with precomputed waveform models, or templates, via standard matched-filtering [6] technique. As a result, the output of the filtering process is particularly sensitive to the phase of GW signal; high accuracy modeling is required to enable faithful source parameter reconstruction and to maximize the physics output of the detections in general.
The waveform templates used in the LIGO/Virgo data analysis pipeline [7] have been coded over decades building on several different approaches, relying both on the analytic and numerical knowledge of the source gravitational dynamics and on phenomenological methods to describe the spacetime deformation around them.
For source modeling, one of the most successful perturbative analytic method has been the post-Newtonian (PN) approximation to GR, see [8] for a review, together with the self-force (SF) approach, see [9] for a recent review, whose small expansion parameter is the ratio of masses of the two objects composing the binary system. Exact numerical methods solving for the entire spacetime had tremendous success (see, e.g., [10] for one of the most recent and complete catalogs of numerical gravitational waveforms from binary systems). Building on results from all these approaches, two main families of accurate templates have been built: one using the effective one body framework [11], which mapped the two-body dynamics into that of a test body immersed in an effective Schwarzschildlike metric (see, e.g., [12] for the latest version), and the other leading to phenomenological waveforms (see, e.g., [13] for latest waveform developments). See also [14] for a review of waveforms used in data analysis.
The present work is intended as overview of the PN approximation framework treated with effective field theory methods, as introduced in [15] and known as nonrelativistic general relativity (NRGR), complementing existing reviews (see, e.g., [16][17][18][19] for the latest development in this rapidly evolving field of research).
In the PN approximation of the two-body problem, the dynamics is expanded around the Newtonian result, with expansion parameter being the relative velocity v, with v 2 ∼ G N M/r according to the third Kepler's law (G N is the Newton's constant, M is the total mass of the binary system, r is the binary constituents' mutual distance, using natural units for the speed of light c = 1), and n-PN corrections correspond to terms of the order G n−j+1 N v 2j , with 0 ≤ j ≤ n + 1. Note that from a point of view of particle physics scattering, a natural expansion parameter would be the coupling constant G N , or rather, its dimensionless version G N √ sq, with s being the invariant mass squared and q being the modulus of the exchanged momentum. Such expansion is known as post-Minkowskian (PM) [20], which has received a great amount of attention recently after the successful derivation of the 3PM dynamics [21,22] and partial results at 4PM order [23,24]. An important feature that PM results have highlighted [25,26] is that computation of the deviation angle in two-body unbound scattering processes seem to be more fundamental that effective action obtained adding the contributions of gravitational mode exchanges, as it is done, e.g., for deriving the Newtonian potential from general relativity, as pointed out in [27].
Another issue that is currently under vigorous investigation is the contribution of radiative modes to conservative dynamics, both in the PN and PM approaches, whose quantitative role is not completely understood yet. In particular, exchange of radiative modes is known to give rise to causal, nonlocal effects, which have been first identified at O(G N ) beyond leading order in emission processes, which have been name hereditary as they depend on the history of the source rather than on its instantaneous state at retarded time. Historically, these have been divided into memory and tail effects [28], the former arising from scattering of radiation onto radiation [29] and the latter from scattering of radiation onto the static background curvature sourced by the total mass of the system [30]. The denominations are related to the nature of the phenomenological effects they have on the waveform: the tail part of the waveform arrives later than the "wavefront", being delayed by the scattering, and then smoothly fades off with time; the memory part is a persistent zero-frequency effect, which is still present well after the wavefront has passed. While hereditary in the waveform, tail and memory contributions are vanishing in the emitted flux [8], and instantaneous (i.e., nonhereditary) in the conservative energy [31] for circular orbits.
Overall, interest in the gravitational two-body problem goes beyond its applications to GW physics, due to its richness in intriguing theoretical aspects, representing a highly nontrivial test bed for classical field theory, including expressing GR as a classical double copy of a non-Abelian gauge theory [32][33][34][35], applying seminal ideas by [36]. The clear separation of scales of the underlying physical problem allows a paradigmatic application of effective field theory techniques and a systematic expansions of interactions in terms of Feynman diagrams. Divergences are encountered in purely classical computations, which are regularized and renormalized, and applications of the renormalization group flow equations can be used to derive nonperturbative results.
The focus of this review is on the conservative, spinless sector for pointlike sources. As far as dissipative effects are concerned, like computation of the gravitational luminosity sourced by a coalescing binary, within the NRGR approach, it has been completed to the 2PN level [37], and the tail-squared quadrupole contribution to the flux has been obtained at 3PN level in [38]. Spin degrees of freedom can also straightforwardly be included in the NRGR approach, as per the pioneer work of [39]. Recent results in the conservative sector reached next 3 -to-leading order level for linear-in spin [40], spin-squared subsec-tors [41], next-to-leading order in the spin-cube [42], and spin-to-the-fourth subsectors [43], whereas those from the flux the investigations reached next-to-leading order in linear-inspin terms [44,45] and spin-squared [46] sectors. (The NRGR method has also been applied to non-GR theories; see, e.g., [47][48][49].) Another topic not treated here is represented by finite size/absorption effects, which have been investigated within NRGR both in the nonspinning [50] and in the spinning sector [51]. Finite size effects have received attention on the phenomenological side as they could provide a handle to distinguish between GR black holes, which, at leading order, absorb energy but do not exhibit tidal deformation, and other exotic types of compact objects [52].
This review is structured as follows: NRGR is introduced in Section 2, with emphasis being given to classical path integral and how to split the integration region of Green's function momenta into a near and a far zone. Then, a summary of the methods permitting the derivation of the conservative dynamics withing NRGR is presented in Section 3, followed by Section 4, which is dedicated to the radiative sector of the problem. The final Section 5 is dedicated to a brief outlook on the currently open problems and imminent future line of investigations in the fields.

Notation
We adopt the mostly plus metric diag(−, +, +, +) and natural units for the speed of light c = 1 but keep the Newton's constant G N explicit. Vectors are indicated in boldface character, e.g., x, and for integration over momenta, the notation k ≡ d d k (2π) d is adopted, with d denoting the number of purely spatial dimensions (d = 3 in our world). Given a quantity A(t, x) in direct space, its Fourier transform and inverse Fourier transform are defined asÃ For spatial indices, we will not make distinction between covariant and contravariant indices, e.g., k 2 = k i k i = k i k i .

The Setup
The starting point to describe the gravitational dynamics of any system is the bulk Lagrangian of gravity, which is given by the gauge-fixed Einstein-Hilbert action where the gauge fixing term involves Γ µ ≡ Γ µ αβ g αβ , with g αβ being the metric and Γ µ αβ being the Christoffel symbol. Λ is a constant with dimension of (length d−2 /mass) 1/2 defined by , with L 0 an arbitrary length needed to ensure the correct dimension of G d .
Action (1) needs to be complemented with the coupling of gravity to world-line sources, whose most general form in terms of mass m, spin S µν and multipole moments of electric type I ijL and of magnetic type J a|Lb can be written as The spin antisymmetric tensor S µν provides a redundant description of the physical spin S i = 1 2 ijk S jk (see [39]); Ω µν are the components of the relativistic generalization of the angular velocity of the object, electric multipoles of 2 (2+l) -order couple to the lth gradient of the electric part of the Riemann, and symbol L denotes a set of l symmetric, traceless indices. Note that for multipoles of the magnetic type, we have adopted the notation introduced in [53], which avoids the use of the Levi-Civita tensor that is not straightforwardly generalizable to the generic d-dimension. Explicit expression for the coefficients c E,Bl can be found in [54,55]. The multipole tensors can be intrinsic to the object or induced by the external field. It is useful to separate the latter from the former case by adding them explicitly to Equation (2), e.g., for the quadrupole where c El,i,Bl,i are constant quantities scaling as mr 2+l+i s , r s is the size of the source, 2G N M for black holes, and ellipses stand for terms with higher number of time derivatives. Actually, all terms with an odd number of time derivatives in Equation (3) are total derivatives, hence they do not contribute to equation of motions derived with time-symmetric boundary conditions. When dealing with an action principle to derive equation of motions, usually the in-out formalism, with time-symmetric boundary conditions, is understood, which, however, cannot account for dissipation, which is exactly what terms with odd number of time derivatives are accounting for (see, e.g., [50,56] for the values of c E2,1 , c B2,1 ). The correct treatment of dissipative terms consists in adopting the in-in formalism [57,58]; however, its description goes beyond the scope of this review. The method for extracting information about time averaged dissipation from the imaginary part of the Fourier-transformed Lagrangian in the in-out formalism will be shown later in this work.
Series (2) and (3) are expansions in terms of r s k and r s ω; k = |k| and ω are the typical wave-number and frequency of the Riemann tensor, and thus the multipolar action expansion parameter is the size of the source divided by the background curvature length, and in the case in which Riemann curvature is determined by radiation, one has k = ω.
The following step is to map actions in Equations (2) and (3) into the description of a definite physical system. When applied to individual black holes, which can be described solely by mass and spin (neglecting charge which is astrophysically insignificant), one can drop all terms in the last integral of Equation (2) but the electric quadrupole one, since Kerr black holes have a spin-induced quadrupole [59] I (spin) ij The induced quadrupoles c E20,B20 have been shown to vanish for a Schwarzschild black hole [60][61][62], and their value for Kerr black holes is subject to controversy at the moment, as there are indications that they are vanishing [63,64], and that they are not [65]. For matter sources such as neutron stars, independently on their spin, individual objects can be endowed with permanent multipoles, as well as induced multipoles of any kind, coupling to the companion's gravitational field. The spin-induced quadrupole multipole can be larger than the black hole one with the same mass by a factor ∼4-8 depending on the equation of state [66], and none of the c E2i , c B2i coefficients are expected to vanish in general. When applied to the binary system as a whole, one can discard Equation (3) and consider the action in Equation (2) as the coupling of the binary system multipoles to emitted radiation. In this case, the multipole expanded action (2) is an expansion over the internal velocity of the source, as k = ω ∼ v/r s , so truncation of the expansion to a finite order correctly represent source with nonrelativistic or mildly relativistic internal velocities.

Computation of the Effective Potential
One can use the coupling of the world-line action for either point particles with no permanent multipoles, i.e., the first two terms in Equation (2), or extended objects, to derive the effective potential ruling the dynamics of a binary system via Gaussian integration. Formally, the effective action S e f f is obtained by integrating the gravitational field mediating the interaction out: where S wl can be taken from Equation (2), possibly adding Equation (3) to take finite size effects into account, and ellipses stand for all possible world-line degrees of freedom (such as spin, acceleration, eventually higher derivatives of the trajectory, multipoles, tidal coefficients, etc.). While our approach is completely classical, it is necessary to normalize the action in the exponent of Equation (5) by a constant with dimensions of an angular momentum, that is,h. However, we will work here at the classical level, i.e., by considering only terms homogeneous in 1/h; thus, all physical results will be independent on the normalization of the action.
Neglecting gravitational interactions, i.e., keeping in S EH only the quadratic, gaugefixed, kinetic term, the Gaussian integral for the continuously infinite variablesh µν (ω, k) can be performed in analogy to the fundamental one-dimensional case According to standard path integral procedure, one can substitute the simple variable x with the gravitational field, , integrating over space-time, and finally taking the logarithm, one obtains the effective action for matter particles, which schematically is In the case of gravity, the source charge J a actually stands for the energy-momentum tensor of particle a, G is the gravitational field Green's function, the t a s are understood to be functions of proper times τ a , and Lorentz indices have been suppressed here for simplicity. (As usual, the energy-momentum tensor is defined by T µν = − 2 √ −g δS wl δg µν .) The effective potential derived here is manifestly symmetric under 1 ↔ 2, as this construction is tailored for determining the conservative dynamics. The natural choice for the Green's function boundary conditions is the Feynman one where "a" stands for an arbitrary small positive number, necessary to ensure convergence of the Gaussian integral with the imaginary unit at the exponent as in Equation (6). Adopting (a)causal boundary conditions, one can define (advanced) retarded Green's functions where the last equality holds for d = 3 only. Introducing the Wightman functions one finds that [67] (To derive (11) the representation of the Heaviside theta function from which one can derive One then finds that Feynman Green's function is complex (both in real and Fourier space), its real part being given by one half of the time symmetric combination G A + G R , and its imaginary part is given by the symmetric combination of Wightman functions known as the Hadamard function whose support in Fourier space is only on the light cone, as it can be seen from Equation (10), whereasG R,A,F , being given by convolutions ofθ(ω) and∆ ± (ω), have support also inside the light cone. Plugging the expression in Equation (12) for G F into Equation (7), one sees that the real part of G F gives the conservative, time-symmetric potential, and the imaginary part of G F determines the imaginary part of the effective action, which is related to the "decay width" Γ of the process, i.e., it represents the rate of particle emission, which is a genuinely quantum concept. However, by writing S e f f has been conveniently expressed in terms of an integral over Fourier variable ω conjugate to time, enabling us to relate the intrinsically quantum objectΓ(ω), which is the ω-spectrum of number of emitted particles, to the energy spectrum of emitted particles, that is,hωΓ(ω), which has a clear classical counterpart in the case of macroscopically large number of emitted gravitons. Equivalently, one can determine the emitted flux by first obtaining the classical gravitational perturbation, which can be computed by solving the linearized Einstein equation for metric perturbations around Minkowski: g µν = η µν + h µν , using the transverse-traceless (TT) gauge for convenience, where Λ ij,kl is the standard TT projector, and then compute the energy flux F via the standard formula See Section 3 of [18] for an explicit demonstration. At the hearth of the equivalence between the two methods lies the optical theorem, which equates the imaginary part of the process of emission and absorption of the same radiation mode with the modulus square of the emission process, i.e., the product of amplitude for emission by the amplitude of absorption. In summary, the imaginary part of the effective action in Equation (7) depends on the imaginary part of the Green's function, which, for the Feynman choice, is Im(G F (ω)) ∝ δ(ω 2 − k 2 ).
To conclude this subsection, we find it useful to recall the explicit form of Equation (15) expanded for small internal velocities of the source: where |x − x | has been expanded in the argument of the energy momentum tensor as r −n · x and collapsed to r in the denominator. The expansion parameter of Equation (17) is v, the same as in the multipolar world-line action Equation (2).

Method of Regions
The main goal of the two-body problem modeling within the PN approximation is determining the dynamics in terms of the two body trajectories (and their derivatives) by integrating the gravitational field mediating the interactions out, i.e., by substituting the gravitational fields with their values determined by the field equations of motion involving the sources. This is efficiently performed via path-integral Gaussian integration as sketched in Section 2.1.
As it is well known [68], gravitational modes separate into longitudinal modes, which are completely specified by the sources, and radiative ones, which are the ones actually propagating with the speed of light and taking energy and angular momentum out of the source.
At leading order, only longitudinal modes contribute to the conservative dynamics, a fact used in standard approach to PN approximation to decompose the contributions of the gravitational interactions into a near and a far zone, corresponding to the scale separation between the size of the binary system r and the gravitational wavelength λ ∼ T ∼ r/v r of emitted radiation, with T being the orbital period.
In the near zone, i.e., at distances shorter than or of the order of r, the source appears as composed of particles, eventually with finite size, whereas in the far zone, i.e., at distances larger than r/v, the binary system is described as a composite object of negligible size and endowed with mass, spin, and all sort of multipoles.
Assuming that the world-line action in Equation (2) describes the source dynamics, the effective potential between two particles with trajectories x a (t a ) can be expressed in terms of the particle energy-momentum tensors (in mixed time/Fourier-space representation) and considering for simplicity only the contribution of the 00 gravitational polarization to the effective action, the only one surviving in the static case, one has The O(G N ) potential in the small velocity limit can be reproduced by expading the Green's function for ω k, i.e., by replacing the expression inside the k integral in Equation (19) with After this substitution, the ω integration becomes straightforward, and it leads to δ(t 1 − t 2 ) terms and derivatives. Note that the expansion parameter (ω/k) 2 of Equation (20) then becomes equivalent to (k · v) 2 /k 2 ∼ v 2 . Expansion in Equation (20) matches the qualitative expectation that the effective potential is mediated by the exchange of off-shell longitudinal modes with ω k, not by on-shell radiative modes with k = ω. The result is an instantaneous potential with v corrections arising from the repeated use of i.e., relativstic retardation effects are reconstructed by the derivatives of the trajectories. Approximation in Equation (20) represents a dramatic simplification in the computations but leaves out the contribution of on-shell modes with k ∼ ω, for which a complementary expansion can be taken in the far zone, characterized by kr < 1: Note that both integrals in Equations (20) and (22) are performed over all values of k, hence none of the two is expected to reproduce quantitatively the exact result of Equation (19) for any finite order truncation of the infinite sums. As it will be shown in a simplified example (see [69] for a formal proof based on [70,71]), the exact result is given by the sum of Equations (20) and (22), provided dimensional regularization [72] is adopted in computing the integrals, which consists of extending the three dimensional integrals to arbitrary d = 3 + dimensions and taking the limit → 0 only at the end of the computation.
The sketch of the proof goes as follows. Let us define and rewrite the original k integral of Equation (19) as the trivial identity withk being a conveniently chosen wavenumber such that ω <k < 1/r, hence requiring that the near and far zone (lenth) scale r and ω −1 must be well separated. Since for k >k (near zone), one has k > ω, and for k <k (far zone), kr < 1, the full integral I in each region can be written as Hence, to demonstrate that the left hand side of Equation (24) is given by the first two terms only on the right hand side, i.e., by the sum of the integral of the near and far zone expansions over the full momentum space, one just needs to show that which holds because both integrands in the left hand side of Equation (26) admit the same parametrization which vanishes in dimensional regularization. We have seen that Equation (20) leads to an expansion in terms of v 2 ; in the far zone, the expansion parameter is ∼ ωr s (with r s being the size of the source, whether it is a single fat object or a binary system), which is the expansion parameter of the multipole part of the action in Equation (2), and for binary systems, one has ωr s ∼ v/r × r = v.

The Conservative Sector
Having presented an overview of the calculation framework, we give in this section some more details about the tools needed to compute the conservative dynamics. We refer to [73][74][75][76][77] for the actual near zone computations at respectively 2, 3, 4, and 5PN order.

Expansion in Terms of Feynman Diagrams
When gravitational self-interactions or nonlinear coupling of the source to gravity are taken into account, the Gaussian integral in Equation (5) involves bulk terms of the type h a , with a ≥ 3, or world-line terms of the type Jh b , with b ≥ 2. It is then convenient to to treat all interactions, including linear ones, perturbatively.
In case of the one-dimensional Gaussian variable in the example of Equation (6), one Taylor expands from the exponential everything but the quadratic term Ax 2 to be left with integration of the type Performing the x integration is equivalent to substituting pairs of x variables with one Green's function for each pair, summing over all the possible ways to form pairs. When dealing with infinite number of fields, such as forh µν (ω, k), the path integral integration has the effect of replacing each pair of fields in the integrand outside the exponential with the inverse of what multiplies the quadratic term |h µν | 2 in the exponent, i.e., the Green's function. This substitution is called Wick contraction.
Examples of representations in terms of Feynman graphs of processes involving the direct linear interaction of two sources, interaction mediated by a bulk, cubic self coupling, and a nonlinear source interaction are shown in Figure 1. Another fundamental concept for the derivation of the effective action in terms of perturbative path integral is the definition of a connected diagram: if by following Green's function lines, all the vertices (both gravitational self-interaction and gravity source ones) can be connected, the diagram is said to be connected; otherwise, it is disconnected. Only connected diagrams contribute to the effective action. This statement will not be demonstrated here, but the following is a sketch of the proof.
Let us collectively indicate with J all the source degrees of freedom in Equation (6) In the expression for S e f f , all terms with n > 1 describe disconnected diagrams, and some disconnected diagrams are also present, along with connected ones, in the n = 1 term. However, the n = 1 disconnected contribution is precisely canceled by the n = 2 terms, and an exact cancellation of all disconnected diagrams happens at all n.
Note that among the connected diagrams, one has to consider both nonfactorizable and factorizable ones, such as the ones in Figure 1, whereas disconnected diagrams such as the one in Figure 2 do not contribute to the effective action. Factorizable diagrams have the advantage that they can be computed via simple algebraic manipulations from lower order ones, a well-known fact that has been exploited, e.g., in [78,79] to show that near zone static contributions to the effective potential at odd PN orders are completely ascribable to factorizable diagrams only, whereas at generic order, factorizable diagrams still make a non-negligible share of the complete set of diagrams, their share growing with PN order [80]. Factorizable diagrams are then intrinsically less fundamental, implying that the effective potential is not the most fundamental object to compute, as part of it can be derived by lower order computations. Remarkably, it was observed in [81,82] that the functional dependence of the scattering angle χ on the total energy and the total angular momentum in a two-body (unbound) scattering process is a more fundamental object to deal with. Indeed, the 3PM dynamics have been rederived in terms of scattering angles [22,25] and a partial result at 4PM order (confined to the pontential mode exchange [24]) confirmed earlier results for the effective Hamiltonian [23].
A crucial observation in the scattering angle approach [27] is that the χ dependence on the symmetric mass ratio ν ≡ m 1 m 2 /M 2 is surprisingly (at least for the author) simple. The Hamiltonian at n-PN order involves terms up to ν n+1 , whereas ν dependence drastically simplifies for scattering angle, as the exchanged momentum ∆p can be schematically written as with b being the impact parameter of the scattering, and q (nPM) i are generic functions not depending on m 1,2 but on the relative velocity only (in the spinless case). Using symmetry under 1 ↔ 2 in the center of mass frame, and substituting the identities m 1,2 /M = (1 ± √ 1 − 4ν)/2, one obtains for appropriate functions X (nPM) i that do not depend of ν nor m 1,2 , thus showing that terms involving ν n+1 , or equivalently nSF results, allow to determine up (2n + 2)PM order scatter-ing angle dynamics, being the 0SF order the limit of a test mass in a Schwarzschild background.
In particular, this simple argument predicts that one can determine up to 4PM dynamics using only the 1SF results, which are well available [83]. This observation has profound consequences for constraining PM and PM results, as it will be shown in Section 3.4, yet it is unclear to the author if the argument of [27] connecting nSF results to (2n + 2)PM ones can be extended to the effect of radiative modes at arbitrary PN order.

Classical vs. Quantum
All the examples in Figure 1 correspond to connected diagrams, two of them being loop diagrams, i.e., involving an integration over internal momenta. However they all describe classical processes.
For a generic graph the number of loops, L is given by L = I − V + 1, with V being the number of vertices and I the number of internal lines (including the solid line for nonpropagating massive particles). Let us denote with I h the number of internal gravitational lines and by I m the number of massive source internal lines, with I = I h + I m . Each Green's function, being the inverse of the quadratic part of (S EG + S GF )/h, brings a factorh; each vertex, as it comes from the expansion of S EH /h, has a factorh −1 . Note that in NRGR, the massive bodies are nondynamical source/sinks of gravitational mode; they are not associated to any Green's function. Thus, for any graph, we have the quantum scalingh whereh −1 is the scaling for classical graphs, and it can be straightforwardly verified that all three diagrams in Figure 1 are classical. The two-body potential V(x) is a classical concept, and one can derive it from the computation of relativistic scattering amplitudes A i→ f for a process taking from some two-body initial state i to a two-body final state f via (see, e.g., chapter 6 of [84]) where q is a wave number (i.e., momentum divided byh) and the nonrelativistic normalization for the quantum states is understood. In the computation of a matrix element between states |p and |p mediated by the energy momentum tensor T µν for a particle of mass m, the relativistically and nonrelativistically normalized amplitudes are related by Then, for the tree-level graph in Figure 1, assuming for simplicity a matter scalar field φ with action S φ given by and using that p|T 00 |p (R) ∼ p 0 /h 2 ∼ (m/h) 2 for small velocities, one finds the scaling of the first diagram in Figure 1 to be which correctly reproduces the classical scaling withh (and masses) of the diagram. The diagram (b) in Figure 1 is slightly more involved, as it includes a loop integral k dω 2π 1 where (0, q) is the overall wave number exchanged between the two objects, and the massive body with two gravitational insertions has initial wave number (E, p). In the nonrelativistic limit E ∼ m/h,hp/m → 0, expression (35) gives [85] h 2m k dω 2π 1 which can be integrated by closing the integration contour in the upper half-plane to pick the contribution from the residues at the poles ω = −|k| − ia and ω = −|k + q| − ia to get This result is proportional to the loop integral in the potential region of momenta (see the Green's function representation in Equation (20)), and it is understood that the region of integration is |k| m/h. This exercise shows that the loop integral in expression (35) has both classical and quantum pieces, the former being isolated by the limithq/m → 0; taking such a limit before performing the loop integration over k kills the quantum part of the diagram.
Going back toh power counting, one gets hence it scales as a classical diagram with appropriate mass powers. As a result, we have shown how classical contributions from loop diagram arise by taking the large mass limit.

Ultraviolet Divergences
When dealing with loop integrals, it is not unlikely to encounter divergent quantities, which need regularization to save predictivity of the computations, and the regularization mechanism naturally embedded in effective field theory methods is dimensional regularization.
As an introductory example, let us analyze the divergences of the fundamental 1-loop diagrmam, which is the standard textbook integral [86] and remember that the Gamma function has simple poles for all nonpositive integers. The poles for a = d/2 and b = d/2 correspond to infrared (IR) divergences of the type ∼ k k 2 −d/2 for k → 0, whereas an ultraviolet (UV) divergence appears for a + b = d/2. The integral vanishes identically when either a = 0 or b = 0, but care is needed for the scaleless integral k k 2 −a , which is always vanishing, but it can be considered as the sum of mutually compensating UV and IR divergences for the specific case a = (3 − )/2 [87]: where the arbitrary length scale L 0 has been added on the right-hand side to adjust dimensions. Any dependence on L 0 vanishes when → 0, i.e., d → 3, and the notation UV,IR has been introduced to identify from which integration region the divergences arise, as it is important to keep track of the nature of divergences, which have different origins and different physical meanings, as we now discuss. (We could have also introduced aL 0IR and L 0UV to clarify which limit the arbitrary L 0 scale originates from.) Note that the divergences of the scaless integral in Equation (40) does not depend on r, so one may be tempted to simply discard it as it does not affect the potential; however it could still enter as a subdiagram of a larger, factorizable diagram overall dependent on r, or involve derivatives of the particle velocity, such as |v| ∼ v 2 /r, thus having an implicit r dependence.
In general, in gauge theories, UV divergences signal the failure to describe the shortdistance features of the physical problem being modeled: they point to some short distance incompleteness of the model. Divergences in dimensional regularization appear as a pole for d → 3 (poles for → 0 can be reabsorbed in dimensional parameter redefinitions); however, since the gravitational coupling has a dependence on the arbitrarily length scale L 0 of the type G d ∝ L d−3 0 , one has so that the effective potential acquires a spurious dependence on the arbitrary quantity L 0 in presence of divergences.
Since the physical result cannot depend on L 0 , there are two main possibilities: 1. the divergences are unphysical, as one can get rid of them by a field redefinition, i.e., a mere change of (field) coordinates, which does not affect observables; 2.
a physical parameter (such as the mass or a multipole moment) acquires a logarithmic dependence on the typical source-observer distance scale, say r, at which the physical parameter is measured, i.e., ∝ log(r/L 0 ), to compensate for the log L 0 dependence from Equation (41).
It turns out the the UV divergences of the near zone description up to 5PN are of type 1. Considering a generic field redefinition such as that changes the action by the amount A necessary condition for UV divergences to be absorbed by such redefinition is that they vanish on the equation of motions, and this is indeed what happens in the conservative sector of the two-body dynamics as explicitly checked up to 4PN order in [88].
Equivalently, using the language of field theory, it is possible to add to the action counterterms to cancel the UV divergences, and since they vanish on the equation of motions, they do not affect any physical observables.
The explicit counterterms respecting the symmetry of the problem are made of curvature invariants and particle trajectory derivatives, and at lowest order, they are with R and R µν being the Ricci scalar and tensor, respectively, and where a term a iµ a µ i has not been included since it is proportional to the square of the equations of motion (as long as geodesic motion is considered), hence it can be considered exactly zero.
The Ricci scalar and tensor appearing in Equation (44) on the ith particle world-line are the ones sourced by the other particle hence they vanish on the equation of motions (i.e., Einstein equation), and a µ i = 0 because of the geodesic equation. (Note that for a spinning particle, one can indeed build a term linear in the Riemann of the type I ij R 0i0j with a permanent quadrupole given by Equation (4).) In [88], the values of the counterterm coefficients have been computed, resulting in the coefficientsc i corresponding to arbitrary finite term addition of the counter terms, which, in general, must be fixed by matching to an observable whose value is affected by their values. However, as stated earlier, these terms are proportional to equation of motions, hence the values of thec i are completely irrelevant.
At second order in the curvature, a possible invariant term is which is exactly of the type of the static finite size effect for spinless object introduced in Equation (3). However, as it is well known [60][61][62], static tidal effects that would be generated by this term for nonspinning black holes vanish in 3 + 1 dimensions. The lowest order diagram in which the Riemann-squared term can enter has the same topology as the (c) one in Figure 1, and considering that R µ νρσ ∼ ∂ 2 h, power counting shows that this diagram with a Riemann-squared insertion would contribute to S e f f /h as Lv 10 /h, i.e., a 5PN term, in the static case, in agreement with the effacement principle [89]. Indeed, the near zone contributions to the two-body conservative dynamics have been derived up to 5PN [77], and no finite size effect has appeared.

Infrared Divergences and Radiative Contributions to the Effective Potential
Infrared (IR) divergences are of a complete different nature, showing that the modeling is inconsistent at large distances. The source of the inconsistency can be readily identified in the near zone Taylor expansion Equation (20): the leading k → 0 behavior of the nth order element of the expansion is ∼ ∂ 2n t e ik·x /k 2n ∼ k · (d 2n−1 x/dt 2n−1 )/k 2n , hence the approximation used in the near zone is bound to create infrared divergences for large enough n. Equivalently, one can see that expansion in Equation (20) fails to capture the contribution of modes with k ∼ ω, which are precisely the radiative ones.
The problem admits a straightforward fix, named zero-bin subtraction [70], which consists of deleting all IR divergences because they do not belong to the near zone region (k ∼ 1/r > v/r), by combining them with the UV divergences of the far zone. Zero-bin subtraction can be adopted as a simple recipe justified by the general description of adding near and far zone integrals outlined in Section 2.2.
The far zone diagrams involve a source endowed with multipoles, which describes the composite object made by the binary system.
The diagram involving the lowest number of G N factors is given in Figure 3a, where a gravitational mode is emitted and reabsorbed, its contribution to the effective action being [31] The diagram in Figure 3a (that contributes to iS e f f ) gives a purely imaginary contribution to the effective action, which, as discussed in Section 2.1, can be related to the emitted energy since we used the Green's function with Feynman boundary conditions in Equation (47). Had we used the causal (retarded or advanced) Green's function, the result of Figure 3a would have been identically zero, as in the limit for a → 0 one has having assumed the branch cut on the negative real semiaxis of the ω complex plane.
Note that had we inserted in the same diagrams mass M or angular momentum L at the source-gravity interaction vertex instead of the electric quadrupole I ij , the amplitude would have vanished in both the real and imaginary part because bothM(ω) andL(ω) have support only for ω = 0, being conserved quantities. SubstitutingĨ ij (ω) with higher order multipoles would give their contribution to the emitted flux, where L denotes collectively l (symmetric and traceless on any pair) space indices, with the numerical coefficients f E l ,B l given by (including those for magnetic multipoles) [54] f E l ≡ (n + 1)(n + 2) n(n − 1)n!(2n + 1)!! , At the next order in G N , one hits on the tail (Figure 3b; see [30]) for the first ever computation of a tail process, which describes the emission, scattering, and absorption of radiation [90,91] This diagram contains both a real and an imaginary part, contributing to the real part of S e f f a UV divergent quantity with its related log(L 0 ω) and a finite piece. The result in Equation (51) has to be combined with the near zone result as per the discussion of the previous section to obtain the complete binary dynamics and give a finite answer by canceling the IR divergence there.
To combine near and far zone results, one needs to express the quadrupole in Equation (51) in terms of individual binary components parameters, according to a standard procedure that goes under the name of matching, that will be deferred to next section. However, the leading order quadrupole expression in terms of binary constituents can be immediately written as which, once plugged into Equation (51), allows to recover the correct finite result for the conservative dynamics, with the UV divergence in Equation (51) exactly compensating the IR divergence from the near zone and the L 0 dependence of the logarithm being canceled, as explicitly demonstrated in [88]. I ij I mn I kl (c) Figure 3. Diagrams describing far zone processes. Subdiagram (a) is the leading order emission/absorption process, which is purely imaginary. Subdiagram (b) represents emission/scattering/absorption of radiation, known as the tail process, where scattering is off the static background curvature generated by total mass M. Subdiagram (c) involves a three-radiation-mode interaction, and it is the leading order memory process [31].
Note that the coefficients of the divergence, of the related logarithm, and of the relative imaginary term are in a simple relationship with one another. This is not an accident of the quadrupole case, but the relationship can be generalized to arbitrarily multipole (see [31,55] for details). This implies that the flux coefficients in Equation (50) fix the tail correction to the flux and also the logarithmic contribution to the conservative dynamics; the logarithmic term coming from high PN order from the near zone dynamics is fixed because its L 0 dependence must match the L 0 dependence in Equation (51).
This shows that the near zone IR and far zone UV divergences are two sides of the same coin; the arbitrary replacement of the full theory integral with the sum of two partial representations introduces spurious divergences that cancel out in the full theory.
In this case, the near zone act as a UV completion of the far zone theory. There is an additional lesson that we can take from this: what if we had at our disposal only the far zone, i.e., the low-energy theory, without knowing its UV completion? If that were the case, we could have rather imposed L 0 -independence of the far zone result and obtain renormalization group equations for physical parameters such as mass, angular momentum, and multipoles. To describe how this works, we need to compute processes involving emission of external gravitational modes; that will be the topic of the next section.
What about the imaginary term in Equation (51)? As it should be clear by now, it appears because of the choice of Feynman boundary conditions for the Green's functions of the two radiative modes. Had one used causal boundary conditions, one would have obtained 1/2(G 2 R (ω) +G 2 A (ω)) instead of G 2 F (ω), which would have led to the same result for the real part, i.e., the contribution to the conservative action. Indeed, using definitions in Equation (11), one has with∆ + (ω)∆ − (ω) = 0 as∆ ± (ω) have no common support, and no imaginary part would have been obtained using causal boundary conditions. For the case of the diagram in Figure 3c, involving the interaction of three radiative modes, things are more complicated and still under investigation at the moment of writing the present work, with apparently mutually inconsistent results in [31,92,93].
In particular, using the argument outlined around Equation (30), at O(G 4 N ), there should be no ν 2 term in the scattering angle formula. However, using the near zone result at 5PN of [77] and the far zone one of [31] or [92], such a condition is not satisfied. This disagreement represents the last obstacle preventing a satisfactorily completion of the 5PN conservative dynamics. (Note that, formally, the tail result of Equation (51) is O(G 2 N ); however, it involves the square of the quadrupole derived three times, hence once expressed in terms of position and velocities of binary constituents, it becomes O(G 4 N ).)

The Radiative Sector
All processeses under consideration so far had the same initial and final states: either two fundamental particles exchanging gravitational interactions or a composite object inter-acting with itself. This section is dedicated to processes involving one external gravitational mode, such as the ones represented in Figure 4.
(c) Figure 4. Leading order (a,b) and next-to-leading order (c) coupling of an external gravitational mode to a composite source indicated by a double line. Black dashed lines indicate longitudinal gravitational modes; green wavy lines represent radiative modes. The (a) coupling is non-radiative, the (b) one is radiative, whose leading tail correction is given by (c).
To express the multipoles in terms of binary constituent parameters, it is necessary to map diagrams of the type of Figure 5 that represent coupling to radiation as seen in the near zone into those of Figure 4. . Diagrams describing leading order radiative process from a spinless binary system. The first one must be supplemented by its mirror image under particle exchange. Green wavy lines represent radiation; the black dashed line stands for longitudinal modes. Figure 5 give the following coupling between the radiative field h ij and the binary system:

Diagrams in
Using the Newtonian equation of motion a 1 = −G N m 2 r/r 3 , one can recast Equation (54) into the standard 1 2 I ij R 0i0j of Equation (2), where I ij at lowest order is given as expected by Equation (52). The derivation of the equivalence between Equation (54) and 1/2I ij R 0i0j , with I ij expressed by Equation (52), is a standard textbook exercise that can be obtained in complete generality, i.e., not necessarily for the binary system case, by repeated use of the energy-momentum conservation law T µν ,ν = 0. (See [37] for matching calculation of multipoles up to electric exadecapole and the magnetic octupole.) The diagram in Figure 4c is the radiative analog of the tail process Figure 3b, and it gives O(G N Mω) correction to the leading order radiative coupling of Equation (54). Relative to the leading order process represented in Figure 4b, it has a finite real part and an IR divergent imaginary part [44] The imaginary part can actually be resummed to a phase and hence drops out of the modulus squared of the (corrected) quadrupole to determine the emitted flux. However, such a phase can in principle leave an observable imprint on the waveform, as its shift is not universal in its finite rational part for all multipoles, while the universal unphysical divergent (and the irrational) part(s) are unobservable.
To compute radiation processes at next order, one has to evaluate the diagrams in Figure 6, which are (G N Mω) 2 corrections to leading order emission process. Such diagrams also present a UV divergence, which makes manifest the failure of the far zone description to treat the binary as a single object up arbitrary small scales. The modulus squared of the sum of diagrams in Figures 4c and 6 gives the (G N Mω) 2 ∼ v 6 , i.e., 3PN gravitational self-interaction (or tail-squared) corrected quadrupole moment entering the flux formula Equation (16), which is [38] |I where the ellipses stand for a finite part and β I l s for any l have been computed in [94] for the electric case and in [95] for the magnetic ones. Tail corrections represent the leading order gravitational corrections to the multipoles scaling with an odd power of v: G N Mω ∼ v 3 . Additional corrections to the multipoles come from near zone nPN corrections to the "matching" diagrams in Figure 5, adding v 2n corrections in the expansion of Equation (17).
Differently from the conservative sector, the near zone computation to resolve this far zone UV divergence has not been performed; however, we have the possibility to use renormalization group equation to know about logarithmic terms in multipoles. After subtracting the infinite part, one can write the renormalization group flow equation for the renormalized quadrupole (which will still be indicated by I ij ) by imposing µ ≡ L −1 0independence of the flux, which is a physical observable, to obtain [38] µ d dµĨ ij (ω, µ) = −β I 2 (G N Mω) 2Ĩij (ω, µ) , which can be solved immediately for constant M to give [38] I ij (ω, µ) = µ µ 0

2I
(1) ij I (5) ij − 2I (2) ij I (4) ij + I By observing that µ dependence in M starts only at G 2 N order, Equation (58) turns out to be a consistent solution of Equation (57). Short-circuiting Equations (58) and (59), one finds [98]: where the brackets denote the time average,μ ≡ µ/µ 0 and M and I ij without arguments are understood to be evaluated at the scale µ 0 . Such a result can be cast into a more interesting form by finding the leading logarithms at all order in v of the energy of circular orbit, which is a relationship between two physical observables: the energy and the frequency of circular orbits. To do so, we must express the energy and the quadrupoles in Equation (60) not as functions of the velocity v and radius r of the orbit, which are coordinate-dependent and have no invariant meaning, but as a function of the invariant orbital frequency ω. To obtain a relationship between r, v, and ω, this last one usually traded for the PN parameter x ≡ (G N Mω) 2/3 , we need an extra piece of information, which is usually provided by the binary constituents equation of motion. Here, however, it is more convenient to use the information that on circular orbits [99,100], the differential of the binding energy E and the differential of the angular momentum modulus J are related by also known as thermodynamic relation, derived from the first law of binary black hole mechanics for constant individual masses. Equation (61), together with the two Equations (60), permit to determine r(x) and v(x), leading to invariant functions M(x) and J(x) containing all the leading logarithms of the invariant energy and angular momentum functions. Remarkably, these functions admit a resummation, and we report here only the one relative to the energy [98] E(x) = − 8mν 2 x 2 15β I 2 1 + 24β I x 3 log x x 4β I x 3 − 1 .
This result for the leading logarithmic terms in the energy of circular orbits involves only up to the next-to-leading order in ν (1SF), hence it can be checked with self-force calculations, from which information about up to the 21PN order has been extracted [101], showing perfect agreement.

Conclusions
The effective field theory modeling of compact binary dynamics is a program currently under active development whose main phenomenological motivation lies in improving the templates used for gravitational wave detections. The program has been developed successfully up to the fourth perturbative post-Newtonian order, where a complete understanding of the interplay of near and far zone dynamics has been obtained, and it has been extended to the radiative sector up to third post-Newtonian order and to the spinning sector.
At fifth post-Newtonian order, the conservative dynamics has been obtained for both contributions of the potential and radiative modes, but discrepancies recently emerged by comparing post-Newtonian, scattering angle, and self-force computations show that a complete understanding of either the interplay between potential and radiative modes, or of the relationship between scattering angle and bound orbit dynamics, is still lacking.
Anyway, it is clear that cross-pollination among different investigation methods is the key to a deeper understanding of the problem that is advancing rapidly, spurred by the observations that are accumulating at ever-increasing pace.
Funding: The work of the author is partially supported by by CNPq under grant n.312320/2018-3.