Effective two-body approach to the hierarchical three-body problem

The motion of three bodies can be solved perturbatively when a tightly bound inner binary is orbited by a distant perturber, giving rise for example to the well-known Kozai-Lidov oscillations. We propose to study the relativistic hierarchical three-body orbits by adapting the Effective Field Theory techniques used in the two-body problem. This allows us to conveniently treat the inner binary as an effective point-particle, thus reducing the complexity of the three-body problem to a simpler spinning two-body motion. We present in details the mapping between the inner binary osculating elements and the resulting spin of the effective point-particle. Our study builds towards a derivation of three-body analytic waveforms.


I. INTRODUCTION
Triple system in nature often come in hierarchical configurations [1]. In this kind of setting, a close (or "inner") binary m 1 -m 2 is orbited by a distant perturber m 3 (the "outer" object). Studies of systems of this kind include satellites and asteroids in the Solar system [2,3], triple stars [4][5][6][7], exoplanets [8][9][10][11] or triple black holes and neutron stars systems [12][13][14][15][16]. The first study of hierarchical triples dates back to Lidov and Kozai [17,18]. They discovered that the triple system evolves also on a characteristic timescale much longer than the period of the two orbits, which is now referred to as the Kozai-Lidov (KL) timescale: where T 3 (resp. T ) is the period of the outer (resp. inner) orbit. On this long timescale, the inner system can experience eccentricity and inclination oscillations, if it starts from a configuration with large relative inclination. Thus, hierarchical systems of black holes can feature dramatically reduced merger timescales [19]. Since these systems are quite common in dense stellar environments [20,21], this makes them especially relevant to gravitational wave astronomy. The conventional treatment of the hierarchical problem proceeds by expanding the Hamiltonian in the ratio of semimajor axes, which we denote by ε: where a an a 3 are the semimajor axes of the inner and outer orbit respectively. Then, one can average the expanded Hamiltonian on both orbital timescales to obtain a set of long-timescale evolution equations of orbital * adrien.kuntz@sns.it † francesco.serra@sns.it ‡ enrico.trincherini@sns.it quantities, a procedure known as adiabatic or secular approximation. The set of equations thus obtained is commonly referred to as the Lagrange planetary equations. The quadrupole term of this averaged expansion gives rise to the KL oscillations. The next level of approximation in a/a 3 , namely the octupole, leads to even richer possibilities like orbital flips, extreme eccentricities and chaotic evolution [22][23][24].
On the other hand, General Relativity (GR) brings corrections to the motion proportional to the relative velocity between the bodies v, which can be regarded as an expansion parameter (we will use units in which c = 1). The interplay between two-body GR effects and Kozai-Lidov oscillations has been studied by numerous authors, e.g. [25][26][27][28]. Generically, the inner binary precession tends to suppress the eccentricity oscillations if the GR timescale is much shorter than the KL one [24]. In other parts of the phase space, though, post-Newtonian corrections combined with the three-body ones can excite eccentricities [29].
However, much less is known concerning the corrections brought by genuine three-body relativistic effects. These terms are essential to derive a waveform of a hierarchical three-body system (they can give rise to the so-called "tidal resonances" [30]) or to obtain the correct time-evolution of the system in some parts of the parameter space [31]. One can approach the problem with a numerical relativistic three-body solver [32][33][34][35][36], however this method is usually time-consuming and inadequate for the derivation of an analytic inspiral waveform model to be used in matched filter analysis [37]. Another valid approach is to further assume a hierarchy of masses m 3 m 1 , m 2 so that the system can be studied with black hole perturbation theory [30,[38][39][40][41][42][43][44][45]. However, in this article we will rather focus on the similar-mass case m 1 ∼ m 2 ∼ m 3 .
Most previous studies on the fully relativistic hierarchical three-body problem use a combination of the post-Newtonian formalism [46] and the quadrupole expansion at the level of the equations of motion [31,[47][48][49][50][51]. Some other authors use instead a Hamiltonian formalism [29,52]. However, their studies do not take into account the so-called adiabatic corrections computed by Will [49,50] and Lim and Rodriguez [51]. The computations needed to get the "cross-terms" representing the interaction of the multipole expansion with relativistic corrections can be quite cumbersome, and it may be difficult to gather physical intuition from the results of the calculations.
In this paper, instead, we begin exploring three-body systems following a new approach, based on a number of powerful Effective Field Theories (EFT) techniques that have been developed for the relativistic two-body problem in recent years. The framework we build on goes under the name of Non-Relativistic General Relativity (NRGR) [53,54], and its generalization to include the spin of the constituents [55][56][57], which leads to the computations of new post-Newtonian terms related to spin [58][59][60][61][62][63][64][65][66][67][68][69]. In the EFT language, the gravitational multipole expansion is implemented at the level of the action [70], using symmetries to restrict the form of the allowed terms [71]. The multipole expansion derived in this way has then been used to compute the gravitational dissipative dynamics in the GR two-body problem.
In the following, we will apply similar ideas to the hierarchical three-body problem. Such a system is particularly suited to an EFT description for two reasons. First, its dynamics is characterized by two small dimensionless ratios of scales, v and ε. The EFT power counting rules allow to estimate easily the sizes of different contributions, thus dictating to what order in perturbation they have to be computed, at a given experimental accuracy. Ratios of different scales can be kept to different orders, depending on their numerical values. In the second place, symmetries are manifest at the level of the EFT Lagrangian and they restrict the form of the allowed terms. As we will see, this considerably simplifies the form of the cross-terms compared to the existing literature and it allows to gather some physical intuition about the effects of relativistic multipole corrections to the dynamics.
The very nature of the Effective Field Theory framework requires to first identify the hierarchy of wellseparated length scales involved in a system and remove (integrate out) each of them, one at a time, starting form the smallest. In this way a tower of EFTs is obtained, that eventually leads to the infrared (long-distance) description of the problem one is interested in.
Thus, we will first focus on the inner binary and integrate out the gravitational field in the presence of an external perturbation, which will be ultimately generated by the third body. The resulting theory will match onto a composite particle, endowed with spin and multipole moments, coupled to gravity. Such a treatment will be valid away from resonances 1 , and as long as the ratio of semimajor axes remains small at all times. This proce- 1 If the perturbation was in resonance with the modes of the inner orbital motion, then it would be much more difficult to integrate out these modes and to describe the inner binary as an effective dure means replacing a three-body problem with a simpler two-body one, where one of the two point-particles is the inner binary.
Here we describe briefly the distinct steps to guide the reader in the rest of the paper.
1. We will start from a system of three worldlines minimally coupled to gravity, where we have already integrated out the modes whose wavelengths are comparable with the size of the bodies. From this starting point, we will integrate out the off-shell modes that contribute to the gravitational potential of the inner binary, having momenta k µ ∼ (v/a, 1/a). Thus we will obtain an action describing the gravitational interaction of the two inner bodies in the presence of an external gravitational field. This will be done in Section II.
2. Then we will first expand the Lagrangian in multipoles and, after that, since we are interested in long-time scale evolution, we integrate out the point-particle orbital modes with frequencies ω > v/a. In practice this will be done by averaging over the period of the inner orbit. Doing so, we will obtain the action of a composite particle, whose spin is simply the orbital angular momentum of the inner binary, coupled to an external gravitational field. This step will be carried out in Section III. Although the final result may seem straightforward from an EFT perspective (gauge invariance fixes all the terms in the action to dipolar order without any free parameter, so that the matching might seem superfluous), the computation will allow us to find the exact relation between the center-of-mass choice and the so-called "spin supplementary condition" (SSC), which is a particular gauge choice for the spin tensor.
3. Similarly to the first step, we will then consider the two worldlines, one for the third body and the other for the composite spinning particle representing the inner binary. We introduce the "effective two-body" EFT and show explicitly its powercounting rules in both expansion parameters v and ε. Integrating out the off-shell modes with momenta k µ ∼ (V /a 3 , 1/a 3 ), we will obtain an action describing the gravitational interaction between the inner binary and the third body.
4. Finally, we will integrate out the remaining pointparticle orbital modes with frequencies ω > V /a 3 , doing an average over the period of the outer orbit. In this way we will get to a Lagrangian representing the dynamics of the 3-body system as an interaction between the composite particle representing point particle. The same obstacle is encountered in the double averaging procedure, see for instance Appendix A2 of Ref. [23].
the inner binary and the outer body. These last two steps will be carried out in Section IV.
Besides these points, in Section II we will also comment on the relativistic definition of the center of mass and introduce the osculating elements that describe the perturbed motion of the binary. We elaborate on the relation between the spin kinetic term and the Lagrange planetary equations in Appendix A, while in Appendix B we provide details about the specific spin supplementary condition used in this article.
For the moment being, we will carry out our computations up to dipole and 1PN order. Already at this stage, we will highlight a number of conceptual clarifications arising from the EFT treatment. However, several interesting new terms also arise at quadrupolar order, related to the corrections to adiabatic approximation [49,51]. While we will briefly comment on the allowed form of these terms (restricted by symmetries) in this paper, we will defer a complete study of them to further work.
We will use the mostly positive metric signature. Planck's mass is defined by M 2 P = 1/(8πG N ). Given the numerous different symbols used in this article, we provide here a dictionary of our notation: • y 1 , v 1 , y 2 , v 2 : positions and velocities of the two constituents of the inner orbit, of masses m 1 and m 2 ; • y 3 , v 3 : position and velocity of the external perturber, of mass m 3 ; • Y CM , V CM : position and velocity of the center-ofmass of the inner binary, defined in Section II B • m = m 1 + m 2 is the mass of the inner binary, M = m 1 + m 2 + m 3 is the total mass of the system, µ = m 1 m 2 /m is the reduced mass of the inner and ν = µ/m its symmetric mass ratio. Similarly, µ 3 = m 3 m/M and ν 3 = µ 3 /M are the reduced mass and symmetric mass ratio of the outer; • Ω, ω, ι: angles characterizing the orientation of the inner orbit (the "orbital elements"), defined byα = R z (Ω)R x (ι)R z (ω)û x where the R xi 's are rotation matrices along the given axis x i ; • u [η]: mean [eccentric] anomaly of the inner orbit; Figure 1: Structure of the "effective two-body" description. The inner binary is replaced by an effective point-particle whose mass and spin are respectively the relativistic binding energy and the angular momentum of the binary system.

II. A BINARY SYSTEM IN AN EXTERNAL FIELD
In this Section, we will obtain the effective Lagrangian at the first post-Newtonian order for the inner twobody system using the background field method, which amounts to integrate out the metric fluctuations in the presence of an arbitrary external field. Up to dipole order, we will then explicitly match this Lagrangian to the one of a spinning point-particle coupled to gravity. This spin coupling induces the dominant non-trivial post-Newtonian evolution of the inner binary parameters in the hierarchical three-body problem.
Before integrating out the gravitational field, let us introduce a convenient notation. We will write the Lagrangian of the binary as where µ = m 1 m 2 /(m 1 + m 2 ) is the reduced mass, r = y 1 − y 2 , v = v 1 − v 2 and L 1 is called the perturbing function. For instance, considering only the Newtonian-order perturbation due to the additional Newtonian potential Φ from the third body, the perturbing function reads where V CM is the (Newtonian) center-of-mass velocity. The aim of this Section is to compute the 1PN terms in the perturbing function.
A. The Lagrangian up to 1PN order In order to make the computations as simple as possible, we will use the Kaluza-Klein decomposition space+time of the metric presented in [72,73], since in the NR regime the time dimension can be considered as compact in comparison to the spatial dimensions. The full metric is decomposed into a scalar φ, a spatial vector A i and a spatial metric γ ij in the following way: We take the field action to be the standard Einstein-Hilbert term with a harmonic gauge-fixing term [53], where Γ µ is the harmonic gauge condition, In the non-relativistic limit and in the conservative sector of the dynamics, temporal derivatives are treated as an interaction term. Up to 1PN order we will only need the part of the action defining the φ and A i propagators, so that the action simplifies to Consequently, the propagators of the (Fourier-space) fields are given in the non-relativistic regime by and there is an additional scalar temporal vertex whose To the Einstein-Hilbert term we add two pointparticles A = 1, 2 whose action is where v µ A = (1, v A ) is the coordinate velocity of the pointparticle. We have set γ ij = δ ij since the fluctuations of γ ij contribute only starting from 2PN order [72]. We expand the point-particle action for weak-field values. At 1PN order, the only vertices contributing are: We now use the background field method by splitting the fields according to φ =φ +φ, A i =Ā i +Ã i . The tilde quantities correspond to an external arbitrary field (later on, we will relate this field to the one generated by the third point-particle), while we integrate out the barred quantities corresponding to gravitons exchanges between the two bodies. The part of the Lagrangian which does not depend onφ andÃ i is the so-called EIH Lagrangian [74]. Since it has already been computed in this framework by several references [53,72,73], we will simply give its expression without explicitly computing the relevant Feynman diagrams: where r = y 1 − y 2 , r = |r| and n = r/r. Next, including background fields, we can compute the perturbing function L 1 defined in Eq. (3), integrating out φ andĀ i . At 1PN order the result is given by: where L 0 was introduced in Eq. (3), and the last term comes from the Feynman diagram with one externalφ and one internalφ, represented in Figure 2.

B. Center-of-mass coordinates
Given the full 1PN two-body Lagrangian in Eq. (14), there remains to expand the two point-particle positions relatively to their common center-of-mass (CM). It is a well-known fact that there is no universal CM definition in General Relativity [75]. For example, the ambiguities in the choice of the CM are related to the socalled "spin supplementary condition" for spinning pointparticles, which is a gauge choice for the spin degree of freedom [56,76]. We provide a discussion about the spin of our system and its relation to the center of mass in Appendix B. In our case, we will adopt the standard post-Newtonian definition of the CM, i.e at 1PN order: Conversely, one can express the coordinates y A using the relative separation r and the CM position X CM : where we have defined and to 1PN order we have: In the absence of any external field, the CM follows a straight line in the post-Newtonian coordinates. However, in the hierarchical three-body problem the binary CM will not follow such a trajectory even at the Newtonian level.
where the field is now evaluated at the CM position Y CM . The monopole corresponds to the term involving no derivatives of the fields, the dipole to the term involving first derivatives of the fields and so on.

C. Osculating elements
Before expanding the Lagrangian (14) into multipoles and perform a matching with an effective point-particle action, we must eliminate an unwanted degree of freedom from the full theory. Indeed, we want to describe the evolution of the binary over a secular timescale, i.e. a time much longer than the period of the binary itself. In order to do so we average all quantities over the quick periodic motion of the binary, which can be approximated with an ellipse. Indeed, if the motion was purely Newtonian, the trajectory would be described by five constants of motion (six, if we count the initial time), which are nicely packaged in a set of geometrical elements. These are respectively the semimajor axis of the ellipse, the unit vector along the angular momentum and the Runge-Lenz vector: There are two angles in the unit vectorγ; furthermore e is orthogonal toγ (it points towards the perihelion) and its norm is equal to the eccentricity e. Conversely, the position and velocity vectors can be written as whereα = e/e,β =γ×α and η is the eccentric anomaly, defined by where φ is an arbitrary initial phase, and u is called the mean anomaly. Now, if the motion is slightly perturbed by post-Newtonian or quadrupolar corrections these constant elements will generically vary slowly with time (compared to the orbital frequency). Thus, in this generic case, we define the osculating elements as the (time-dependent) values of a, e, φ,α andγ such that the instantaneous position and velocity of the binary is given by the formulae (21). This physically corresponds to drawing at each point the ellipse defined by the instantaneous position and velocity of the binary. We have mapped the six components of r, v into six elements a, e, φ,α andγ.
The equations of motion for the binary system can then be translated in a set of first-order equations on the osculating elements, called the Lagrange planetary equations (LPE). For completeness, we recall them in Appendix A. For our present purposes, though, it will be sufficient to state the result of Eq.(A12), i.e. that the orbit-averaged LPE are completely equivalent to a spin kinetic term in the Lagrangian in flat spacetime: where J is the total angular momentum of the binary and Ω is an angular velocity defined by Finally, we will average all quantities in the Lagrangian over one period of the binary, using the formula valid to lowest order for any quantity A (we recall that η is the eccentric anomaly defined in Eq. (22)). Thus, we will have removed from the Lagrangian the high-energy degree of freedom contained in the mean anomaly. As a consequence of the LPE (A4), the semimajor axis a will be conserved. This can be understood as deriving from the fact that the conserved conjugate momentum associated to the mean anomaly depends only on a. As a side remark, notice that Eq. (25) is valid only if we assume that the binary exactly follows an ellipse. As explained in App. A, there will be higher-order corrections to this formula, which however are not needed for our present purposes.

A. The internal Lagrangian
To begin with, let us deal with the very first term in the Lagrangian (14), namely the EIH Lagrangian. This term does not contain any coupling to the gravitational field. As explained before, it still contains the shortdistance degree of freedom from the Kepler trajectory of the binary system. In order to remove it and keep only the long-distance degrees of freedom which can be excited by the external field (in other words, the osculating elements), we should average the Lagrangian over the inner binary timescale, splitting the variables between the center-of-mass and the relative variables.
A priori, one should be careful about the fact that in the Newtonian kinetic energy one should use the relativistic center-of-mass definition in Eq. (15). However, the meaning of the supplementary 1PN term will be better understood in terms of spin coupling, so we defer its calculation to a later Subsection. Thus, in this Subsection we stick to the Newtonian definition of the center-ofmass. Carrying out the heavy but straightforward computations, we find by using Eq (25): where we have dropped an unimportant constant term in the average (depending on the semimajor axis a only, which is constant in the adiabatic approximation). Each term in Eq. (26) lends itself to a very simple interpretation. The two first terms are just the usual relativistic expansion of the center-of-mass velocity −m 1 − V 2 CM . The third term is the average of the EIH Lagrangian of a binary system in isolation: used in the LPE equation (A8), it gives rise to the celebrated perihelion precession formula. We call such a term the "internal" Lagrangian L internal : Finally, the meaning of the last term in Eq.(26) will become clearer in the next Subsection.

B. Monopole
Starting from Eq. (14) we can collect all the terms coupling the binary system to the monopole of the external gravitational field: whereφ andÃ are evaluated at the CM position Y CM . Note that in the term multiplyingφ in the above equation, we have used the Newtonian version of the CM (i.e., we have set δ = 0 in (16)) since the terms involving δ are of higher post-Newtonian order. Averaging over the binary orbital timescale, we find Let us now gather this monopole coupling together with the average of the EIH Lagrangian (26) computed in the last Subsection. To 1PN order, we find To 1PN order, the last term can be exactly accounted for by replacing the mass m of the binary system (which is now treated as an effective point-particle) with its total energy in the worldline Lagrangian: where E is defined as Thus, the binary moves in the external field with a total mass equal to its binding energy, as could have been anticipated from an EFT perspective [57]. However, our computation highlights the fact that one should also include the internal Lagrangian in the effective action so that the binary PN precession effects are taken into account.

C. Dipole
Expanding the Lagrangian (14) to dipole order (i.e, to first derivatives in the external fields) by taking into account the relativistic CM definition (15), we find at 1PN: where X A = m A /m. As before, one should average this Lagrangian over the inner binary timescale. We find that the term proportional to the difference of masses averages out, leaving us with an averaged Lagrangian From this expression one easily recognizes the coupling of a spinning point-particle to gravity given in e.g. [56,57]. In our case, the spin tensor J µν depends on the total orbital angular momentum of the binary system J = µ G N ma(1 − e 2 )γ through the relations The second condition is called a spin supplementary condition, removing the unwanted degrees of freedom from the full spin tensor J µν . As mentioned before, this gauge condition is related to the choice of a center-of-mass of the binary system; our particular CM choice in Eq. (15) has selected the spin supplementary condition J 0i = 0, which has already been discussed e.g. in [56,77]. We further elaborate on this in Appendix B. Furthermore, note that in Eq. (35) the spin tensor has been projected to a locally flat frame through J ab = e a µ e b ν J µν , where we have introduced the worldline tetrads defined over all spacetime byg µν e µ a e ν b = η ab . As a side remark, note that on top of the spin supplementary condition, the components of the spin vector are not all independent degrees of freedom following the remark below Eq. (A13). This reflects the fact that the spin of the inner binary contains two degrees of freedom once the orbital timescale has been integrated out, instead of the three degrees of freedom contained in the Euler angles of a generic spin.
To 1PN order, the spin coupling (34) can be written in a compact form using the Ricci rotation coefficients: This formula gives back our previous equation (34) when expanded for weak-field values [57,67]. We may be tempted to add to this spin coupling the kinetic term for the spin in Eq. (23) to obtain the minimal gravitational spin coupling which has been discussed at length in the NRGR formalism [55,56]: In this equation the total angular velocity Ω µν includes both the Ricci rotation coefficients from Eq. (36) and the locally flat angular velocity from Eq. (23). It is defined through Here Ω ab flat is related to the tensor Ω ij = ijk Ω k by a relation that we discuss in Appendix B, and the rotation vector Ω =α ×α has been defined in Eq (24). However, there is a small piece that is still missing to obtain the full Eq. (37), related to the choice of the center-of-mass. As we show in Appendix B, in the spin gauge we are using (J 0i = 0), there should be a supplementary spin kinetic term related to Thomas precession, which is 1PN order higher than the kinetic term (23): where A CM is the acceleration of the center-of-mass. Such a term is related to the PN corrections to the centerof-mass position (and speed) which we ignored in Section III A. Indeed, using the full CM definition (15) in the Newtonian part of the EIH Lagrangian (13) gives a supplementary 1PN order term, where δ has been defined in Eq. (17). At first sight, we may be tempted to discard such a term when averaging out the internal binary timescale. However, one should not forget to take also the time derivative acting on V CM in δ, giving rise to which is exactly the additional spin kinetic term shown in Eq. (39).

D. Quadrupole
From the EFT point of view, at 1PN quadrupolar order the couplings to gravity are contained in two non-minimal worldline operators [71]: where I ij and J ij are the electric-type and magnetictype quadrupole moments of the source, coupled to the corresponding parts of the Weyl tensor C µναβ : Furthermore,in Eq. (42) the tensors have been projected to the locally flat frame defined below Eq. (35). We could proceed as before and carry out the integration procedure to obtain the quadrupole moment of the effective point-particle. However, at this order the procedure is somewhat more involved than one could naively expect. The first complication comes from the corrections to the time averages introduced in Eq. (25). Indeed, post-Newtonian corrections to the period of the system will matter when taking the average of the Newtonian quadrupole moment, combining to produce a quadrupolar 1PN term. In the same way, the Newtonian quadrupolar corrections to the motion of the inner binary should be taken into account in the average of the EIH Lagrangian.
The second complication comes from the corrections to the adiabatic approximation mentioned in the introduction. Indeed, in our analysis we are assuming that all the variables of the inner binary vary on long timescales (except of course the mean anomaly). This neglects shorttimescale oscillations, which can ultimately have an effect on long-wavelength modes [78,79]. It turns out that at lowest order this effect produces cross-terms of 1PN quadrupolar order [49,51] (no such corrections appear at lower multipole orders). While noting in passing that these kind of corrections have a very transparent meaning in the EFT language (they are high-energy corrections to an effective low-energy action), we will defer their complete calculation to further work.

IV. INTEGRATING OUT THE OUTER BINARY TIMESCALE
Now that we have replaced the binary system with an effective point-particle, we can integrate out the external fieldsφ,Ã in the presence of a third point-particle of mass m 3 . For simplicity, in the following we will assume this mass to be of the same order of the mass of the inner binary: m 3 ∼ m. We will first derive the Feynman rules of the effective point-particle; then, in a second step, we will integrate out the outer binary timescale and comment on the different terms obtained in the expansion of the Lagrangian. For the 1PN precision we aim to, it will be sufficient to set the total Newtonian center-of-mass of the three-body system to the origin of coordinates (it will accelerate only at 2PN order [80]). Thus, we will have the expressions where we recall that Y CM is the position of the centerof-mass of the inner binary, and we have defined R = Y CM − y 3 , R = |R|, N = R/R, M = m 1 + m 2 + m 3 , X 3 = m 3 /M and X CM = m/M . The averages over the outer binary timescale are then taken in the same way than in the preceding Section.

A. Power-counting rules
Let us recap what we have learned so far and set up power-counting rules for the vertex coupling the binary system (now treated as an effective point-particle) to gravity. Up to dipole order, the Lagrangian of the binary system can be written as Note that this Lagrangian has not yet been averaged over the period of the outer orbit T 3 , and can therefore describe the secular dynamics on timescales shorter than T 3 . With such a simple Lagrangian, one can assign the standard power-counting rules of NRGR which have been described in e.g [53,54,81], considering the motion of the effective point-particle and the third mass (the outer orbit) for which one has V 2 CM ∼ v 2 3 ∼ GM/a 3 where a 3 is the semimajor axis of the outer orbit. Thus, spatial derivatives are treated as ∂ i ∼ a −1 3 . Time intervals scale as t ∼ a 3 /V CM and the metric perturbations scale asφ ∼Ã i ∼ V 1/2 CM (M P a 3 ) −1 . As usual in NRGR, the lowest-order Lagrangian scales as the orbital angular momentum of the outer orbit J 3 ∼ M V CM a 3 which is treated non-perturbatively, higher-order corrections coming with higher powers of V CM . However, one difference with respect to the standard NRGR power-counting rules is evidently the presence of two expansions, the first one in v and the second one in ε ≡ a/a 3 . A priori, we could also have an expansion in the post-Newtonian parameter of the outer orbit V CM . However, not all these parameters are independent. We choose to write all the post-Newtonian corrections as an expansion in the velocity of the inner binary v, converting the center-of-mass velocity by means of the relation V CM ∼ vε 1/2 , which holds when m ∼ m 3 . In Table I we give the power-counting rules of the monopole and dipole vertex which we computed in Sections III B and III C. The effect of post-Newtonian corrections on the dynamics of the system is highly non-trivial, as it can lead to suppression as well as enhancement of the Kozai-Lidov oscillations depending on the part of parameter space explored [29]; we expect that our power-counting scheme will help in discriminating between the different behaviours observed.

Operator
Rule For convenience, the integral over time is not displayed, although it should be included to obtain a dimensionless rule.
Notice that the scaling of the spin is somewhat different than the one usually presented in NRGR [55,67]. Indeed, when taking compact objects as point-particles the spin is given as an order-of-magnitude by where v rot is the rotation velocity of the object and r s ∼ G N m its size. As a consequence, the ratio of the spin coupling presented in Eq. (34) to the Newtonian gravitational coupling is of v 3 (1.5PN) order. However, in our case the spin order-of-magnitude is given by J ∼ µ √ G N ma so that the ratio of (34) to the Newtonian coupling is Thus, the inner binary angular momentum coupling is formally of 1PN order, although it is suppressed by the small ratio ε 3/2 . This power-counting is different from the one of the Lense-Thirring precession caused by the intrinsic spin of the objects, which has been studied in [47,82] and enters at 1.5 PN.

B. Monopole
Let us begin by integrating out the vertex contained in the monopole operators of the effective binary system, i.e in the square root appearing in Eq. (45). The final effective action, including orders of J 3 v 2 ε, is given by: where L CM,3 EIH is the EIH Lagrangian of the system composed by the CM (of mass E, defined in Eq. (31)) and the third particle.
The Lagrangian in Eq. (48) involves a non-trivial coupling between the variables of the inner and outer binaries, given by This contribution is of order v 2 with respect to the standard Newtonian term L 0 ∼ G N M/a 3 . We average this term over one orbit of the outer binary, which gives where M = m 3 + m, X 3 = m 3 /M and ν 3 = mm 3 /M . This new monopole coupling has no effect on the dynamics. Indeed it depends only on the semimajor axes a and a 3 . Consequently, in the Lagrange planetary equations this term will only enter in the equation for the mean anomaly (A7), which is irrelevant in the adiabatic approximation. Therefore, at the level of the monopole, the resulting motion is the one of two ellipses precessing because of standard two-body GR effects. In fact, one can be quite generic about the monopole terms. Indeed, the only planetary elements upon which the monopole terms could depend are the semimajor axes a, a 3 and the eccentricities e, e 3 (they do not involve angles). In the LPE the derivatives with respect to these elements enter only in the equations for the mean anomaly (A7) and the perihelion angle (A8). Thus, the only effect that monopole terms can have is to make the ellipses precess.

C. Dipole
In order to integrate out modes contributing to the potential at dipole order, we have to compute the diagrams related to spin-orbit coupling. These are shown in Figure 3. Using the Lagrangian averaged over the inner orbit Eq. (34), we find: At this order of approximation however, we should also take into account the Thomas precession term of Eq. (41). This gives a contribution of the same size of the terms in Eq. (51). We can replace the center-of-mass acceleration in Eq. (41) using the equation of motion, since the difference between the two terms would contribute at a higher PN order (this is usually called the "double zero trick" [83][84][85]). Thus, at order J 3 v 2 ε 3/2 the full Lagrangian is given by which recovers the result already known in the NRGR approach [67]. Carrying out the average over the outer binary timescale in a way very similar to the previous Section, we find where J 3 is the angular momentum vector of the outer orbit, J 3 = µ 3 (G N M a 3 (1 − e 2 3 )) 1/2γ 3 (hereγ 3 is the unit vector along the outer orbit angular momentum, and µ 3 = mm 3 /M ). Thus, this term is indeed a coupling between the angular momentum vectors of the two orbits.
From this expression one can obtain a precession equation for the inner orbit angular momentum. Indeed, varying the kinetic term for the spin with respect to the canonical variablesα and W = J ×α, one obtains the equations of motion where the precession frequency is equal to From these two equations, and using the Jacobi identity for the cross product, one obtains the precession equation which is in complete accordance with earlier results on the hierarchical three-body problem [31,51]. Notice that conservation of the total angular momentum requires that J 3 satisfies an analogous equation, In particular, it was shown that this angular momentum precession may play an important role for stellar-mass binary mergers near a supermassive BH [31]. Quadrupolar terms would lead to further precession effects, of order J 3 v 2 ε 2 in the Lagrangian. We leave the computation and the astrophysical implications of such terms to further work.

V. CONCLUSIONS
The NRGR approach to the two-body problem was designed to deal with extended compact objects. In this article, we have extended NRGR to the setting of a hierarchical three-body problem. In the approximation that the inner orbit is much smaller in amplitude than the outer one, the inner binary system can be replaced by an effective point-particle endowed with multipole moments, which we explicitly computed up to dipole order. This is very natural from the EFT perspective and provides a new specific example of how an extended (and not so compact) system can be accounted for by means of a point-particle operator.
Our procedure consists in integrating out the short timescales associated with the period of the two hierarchical orbits. One notable result of our study is to make explicit the link between the Lagrange planetary equations, describing the long-time evolution of the inner binary Keplerian parameters, and the kinetic term for a spin in the EFT language. We have also clarified the relation between the post-Newtonian definition of the center-of-mass and the spin supplementary condition for the angular momentum of the inner binary. The computation of quadrupolar post-Newtonian terms including the corrections to the adiabatic approximation will be the subject of a future publication.
Our study moves towards a more systematic characterization of the relativistic hierarchical three-body problem. Indeed, the EFT techniques that we employed can be applied to study efficiently three-body trajectories to higher orders in both PN and multipole expansions. Two immediate fields of application will be the study of the influence of relativistic three-body interactions on the Kozai-Lidov mechanism, and the production of threebody analytic waveforms in the PN regime using the effective two-body description. Another interesting followup would be to obtain (in a matching procedure) the multipole structure of the inner binary system to higher orders from a numerical relativistic three-body solver. Finally, while we have restricted here the discussion to objects of similar mass, it could also be interesting to generalize our work to the case where m 3 m 1 , m 2 , which is particularly relevant for binary BHs orbiting a supermassive BH at the core of a galaxy. While it would be very easy to include this new large parameter in the powercounting rules in Table I, it would probably be more efficient to take advantage of the large ratio of masses in order to explore the interplay between BH perturbation theory (for the outer orbit) and post-Newtonian EFT treatment (for the inner orbit). We plan to explore these avenues in a near future.
This Appendix introduces a set of equations initially introduced by Lagrange. To begin with, note that the time-dependence of the osculating elements defined by Eq. (21) cannot be arbitrary. We must impose a gaugefixing condition such that the velocity is indeed given by Eq. (21). We denote such a condition by where the expression for the vector v was given in Eq. (21). Thus, there is a relation between time derivatives of the osculating elements. This gauge-fixing condition removes three degrees of freedom (equivalently, six variables in phase space) from the six degrees of freedom contained in the six osculating elements (equivalently, twelve variables in phase space). Now, we could write a Lagrangian for the osculating elements by implementing this constraint with a Lagrange multiplier λ, so that where L 1 has been defined in Eq (3). From there one could deduce the Lagrange planetary equations (LPE) which relate time derivatives of the osculating elements to the perturbing function L 1 . However, it is much easier to derive them in a Hamiltonian formalism, see e.g. [86] to which we refer the reader interested in the details of the derivation. The LPE are traditionally expressed using the following angles: ι is the inclination, ω the argument of periapsis, and Ω the longitude of the ascending node. In term of these, the unit vectorsα andγ are expressed aŝ whereû x ,û y ,û z are the Cartesian basis vectors. Using these angles, the LPE are given by [86] whereL 1 = L 1 /µ. It can be easily checked that the LPE can be derived from the following fist-order Lagrangian 2 The conjugate momenta are given by This Lagrangian has the nice property to be exact (no secular approximation has been done, it is completely equivalent to the original Lagrangian (3)). However, it is not manifestly invariant under a rotation of the basis vectors; such a manifest invariance can be recovered by noticing that the angular part can be rewritten as where J is the total angular momentum of the binary and Ω is an angular velocity defined by Thus, the angular kinetic term can be identified with a spin coupling in flat space (note that our sign convention for the metric is different from the one used in e.g Refs [56,57], which explains the sign difference of the kinetic term). However, note that not all the components of the spin vector are independent, since the Lagrangian shown in (A12) displays only two degrees of freedom (corresponding to four equations in phase space once a variational principle is applied). Indeed, notice that if one wants to vary the Lagrangian with respect toα andβ in order to keep a manifest rotational invariance, one should also impose that these vectors should be unitary and orthogonal in order to preserve the right number of degrees of freedom. Finally, the LPE are often averaged over the periodic motion of the binary system: this is called the adiabatic or secular approximation. This corresponds to eliminating the short-distance degree of freedom contained in the mean anomaly u; as a consequence, since the perturbing function does not depend on u any more, the semimajor axis a is constant through time from Eq. (A4). Thus, the two-body Lagrangian shown in Eq.(A10) is indeed equivalent to a spin kinetic term, since the term G N m/2a+Lu becomes an irrelevant constant in the adiabatic approximation.
After this elimination, the binary system is described by four dynamical quantities (the eccentricity e and the three Euler angles defined above) which vary over a timescale much greater than the period of the binary. Technically, we use the formula valid for any quantity of interest A: However there will be higher-order corrections to these quantities as is implied by Eqs. (A5)-(A7). These corrections contribute at the quadrupolar 1PN level, which is beyond the scope of this paper. They will be investigated in more details in a forthcoming publication, along with the corrections to the adiabatic approximation. In this appendix we provide some details of the computation of the spin kinetic term (37) as a function of the intrinsic angular momentum of the inner binary. The computations are analogous to those carried out in [56], with the difference that we specialize to the no mass dipole gauge in which the time components of the spin tensor are set to zero. This choice will make simple to connect the spin tensor to the orbital angular momentum.
First of all, it is useful to introduce a worldline tetrad e µ A (σ) defined only on the worldline y µ (σ) (σ being the affine parameter of the curve) which represents a choice of axes in the rest-frame of the body and satisfies: g µν (y(σ))e µ A (σ)e ν B (σ) = η AB . This tetrad can be used to define the angular velocity of the body: Ω µν = e ν A (De µA /Dσ), whose conjugate is the spin tensor J µν = 2∂L/∂Ω µν . Both these tensors contain gauge degrees of freedom, since only the spatial orientation of the worldline tetrad has a physical meaning. In fact we can choose arbitrarily its time-like direction, encoded in e µ [0] . This gauge choice corresponds to a redundant boost transformation of the worldline tetrad (in order to avoid ambiguities between the different set of indices, we are using square brackets to distinguish the flat indices of the worldline tetrad from the others).
The gauge fixing of e µ [0] must be supplemented with a gauge fixing of the conjugate variables in J µν , the so-called Spin Supplementary Condition (SSC). Starting from a covariant gauge choice in which e µ [0] = p µ / −p 2 and the spin tensor satisfies the covariant SSC J µν p ν = 0 , the action of a boost will change the worldline tetrad and the angular velocity tensor. The changes produced by this transformation in the Lagrangian can be interpreted by means of a redefinition of the spin tensor and a consequent change of the SSC. We will therefore use this boost degree of freedom to first pick the no mass dipole SSC for the spin tensor and then to fix the canonical gauge for the angular velocity vector. In this way we will get an expression dependent only on the intrinsic angular momentum of the binary.
As computed in [56], the transformation of the spin kinetic term of the Lagrangian under a boost of the worldline tetrad (starting from the covariant gauge and SSC) is the following: where we have used hatted symbols to label boosted variables and in particular we have defined the boosted spin tensor to be J µν =Ĵ µν − δz µ p ν + δz ν p µ , with δz µ =Ĵ µρ p ρ /(−p 2 ). We can interpret this change of the spin tensor as due to a shift of the center of the body rotation, that is the point where the worldline intersects the body. In the case of the no mass dipole gauge, in which the spin tensor is purely spatial, this shift corresponds to setting the center of the worldline on the relativistic center of mass, as shown in the main text. The second term in Eq.(B1) will instead contribute to the Thomas precession, which we can understand as due to a gravitational torque associated to the finite size of rotating objects in GR.
Before specifying the boost needed to get to the desired SSC, it is useful to disentangle the gravitational field from the spinning degrees of freedom. We can do so by introducing the gravitational tetrad field g µν (x)ẽ µ a (x)ẽ ν b (x) = η ab , which is defined on the whole space-time. This tetrad can be related to the worldline tetrad by means of a Lorentz transformation:ẽ µ a (y(σ)) = Λ A a (σ)e µ A (σ), being Λ A a (σ) a Lorentz matrix dependent on the affine parameter of the worldline. As for the worldline tetrad flat indices, when needed we will use round brackets to distinguish the flat indices of the tetrad field from the others.
In this notation, once the gauge of the tetrad field is fixed, we can fix the time-like vector of the worldline tetrad by choosing the boosted zero components of the Lorentz matrices:Λ a . Moreover, introducing the tetrad field will make possible to write all the objects in the right hand side of Eq.(B1) in terms of their counterparts with flat indices. Such quantities correspond to those computed in terms of the intrinsic angular momentum of the binary, as they are independent on the external gravitational field. In particular we have: where ω ab µ =ẽ a ν ∇ µẽ bν is the spin connection of the tetrad field, u µ = dy µ /dσ is the worldline speed and we have definedΩ ab f lat =Λ b A dΛ aA /dσ.
At this point we can fix the gauge boost of the worldline tetrad. In order to set the time components of the spin tensor to zero,Ĵ a(0) = 0 , we need to choose a boost such that p 2Λ [0]a = 2p 0 δ 0a −p a (this can be understood by inspecting the generic expression forĴ µν , as discussed in [56]). Doing so, we obtain the following: This gauge choice makes possible to unpackΩ (i)(j) f lat and expressΛ a [0] dΛ b[0] /dσ in terms of the momentum of the worldline, leaving to compute only on the spatial part of the Lorentz matrices. However, these spatial leftovers won't be SO(3) matrices, since they need to satisfy the conditionΛ a A η ABΛb B = η ab and will carry a dependence on the worldline momentum, due to the gauge condition onΛ a [0] . In order to obtain an angular velocity tensor defined in terms of rotation matrices and to remove its dependence on the worldline momentum, we can take a further boost of the worldline tetrad. This time however, we will not use a redefinition of the spin tensor to absorb the new terms appearing in the Lagrangian after the transformation. Rather, we will retain the spin tensor satisfying the no mass dipole SSC and we will keep track of the new terms explicitly.
In order to makeΛ a A an SO(3) matrix, we need to choose a gauge in whichΛ a [0] = δ a 0 . Therefore we implement a boost of the worldline tetrad that sends the time-like unit vector (2p 0 δ a 0 − p a )/ −p 2 to δ a 0 . This transformation will change only the first term in Eq.(B3) as follows: where now Ω (i)(j) SO (3) is build out of rotation matrices and we have used p a = mu a / √ −u 2 , with −u 2 = 1 at leading order in the PN expansion.
Having fixed the gauge for both angular velocity and spin tensor, we can carry out the explicit computation of the last two terms in Eq.(B3). In order to do so, we pick the tetrad field in such a way to haveẽ (i) 0 = 0. Then, at 1PN order we obtain: Plugging these results into Eq.(B1), and identifying the worldline with the trajectory of the center of mass, u (i) = V (i) CM we finally get: Then, with a mild abuse of notation, we can drop the index brackets and the hats so as to match the expressions used (for simplicity) in the main text: → Ω ij . We stress however that these are different from the (µ, ν) = (i, j) components of J µν and Ω µν , which depend on the external gravitational field.
Using this notation and the definitions J ij = ijk J k , Ω ij = ijk Ω k , we can rewrite Eq.(B6) as: which is the equation used in the main text.