Third post-Newtonian effective-one-body Hamiltonian in scalar-tensor and Einstein-scalar-Gauss-Bonnet gravity

We build an effective-one-body (EOB) Hamiltonian at third post-Newtonian (3PN) order in scalar-tensor (ST) and Einstein-scalar-Gauss-Bonnet (ESGB) theories of gravity. The latter is an extension of general relativity that predicts scalar hair for black holes. We start from the known two-body Lagrangian at 3PN order, and use order-reduction methods to construct its ordinary Hamiltonian counterpart. We then reduce the conservative two-body dynamics to the (nongeodesic) motion of a test particle in an effective metric by means of canonical transformations. The resulting EOB Hamiltonian is a modification of the general relativistic Hamiltonian, and already at 3PN order, it must account for nonlocal-in-time tail contributions. We include the latter beyond circular orbits and up to sixth order in the binary's orbital eccentricity. We finally calculate the orbital frequency at the innermost stable circular orbit (ISCO) of binary black holes in the shift-symmetric ESGB model. Our work extends F.L. Juli\'e and N. Deruelle [Phys. Rev. D 95, 124054 (2017)], and it is an essential step toward the accurate modeling of gravitational waveforms beyond general relativity.


I. INTRODUCTION
The observations of gravitational waves (GWs) from coalescing binary systems composed of black holes (BHs) and neutron stars (NSs) [1][2][3][4][5] with the LIGO and Virgo detectors [6,7] offer the unique opportunity to unveil the nature of these compact objects and to test Einstein's theory of general relativity (GR) in the highly dynamical strong-field regime [8][9][10][11][12]. The GW signals are at first "chirps" produced during the long inspiral phase, where the two bodies steadily and adiabatically come closer to each other, losing energy because of GW emission. The inspiral is followed by a short plunge and merger stage, where nonlinearities prevail, and then by the so-called "ringdown" phase for binary BHs [13,14], or by more complex pre-and postmerger signals (depending on the equation of state of the NS and on the properties of the BH) for binaries comprising at least one NS [15,16].
Among the simplest modifications of GR, scalar-tensor (ST) theories add one massless scalar degree of freedom, which couples universally to matter. They were introduced by Jordan, Fierz, Thiry, Brans and Dicke [77] and put in a modern perspective in Refs. [34,78,79]. The corresponding two-body dynamics has been computed within the post-Newtonian (PN) formalism [34,[37][38][39][40][41]80]. Interestingly, compact objects in ST theories can undergo a phase transition associated with the spontaneous symmetry breaking of the scalar field near the compact object in the presence of large curvature or relativistic matter [35]. For NSs, this phase transition leads to a rapid growth of the scalar charge ("spontaneous scalarization"). An analogous nonperturbative phenomenon ("dynamical scalarization") was found in binary NS and NS-BH simulations in NR [53,55,81]. Various methods to describe these nonperturbative effects in waveform models have been proposed [25,82,83].
However in ST theories, vacuum BH solutions are the same as in GR. By contrast, Einstein-scalar-Gauss-Bonnet (ESGB) theories have attracted particular attention because they have the interesting property that (i) for certain functional forms of the coupling constant, BH solutions in ESGB gravity are different from the solutions of GR, but admit the ordinary Kerr solutions as a special limit; and (ii) there is the possibility of "spontaneous scalarization" [84][85][86] (i.e., BHs can "grow hair"). These observations opened up a much richer phenomenology for binary BHs [87]. Recent progress in gravitational waveform modeling within ESGB gravity includes the calculation of inspiral waveforms using PN theory [48,50,51] and the first calculation of quasinormal mode frequencies of rotating ESGB BHs at quadratic order in a smallspin expansion [88][89][90]. The numerical calculation of merger-ringdown waveforms in ESGB gravity has also made remarkable progress, at first using a small-coupling approximation to numerically solve the field equations in an "effective field theory" approach [57,[60][61][62]68], and then by showing that numerical evolutions are possible in the full theory, although hyperbolicity can break down in some regions of the parameter space [59,63,64,66,69].
An important step to build semianalytic IMR waveforms is to construct an accurate analytic description of the two-body conservative inspiral dynamics. We achieve this here within the effective-one-body (EOB) formalism [91][92][93]. The EOB approach builds IMR waveforms by combining analytical predictions for the inspiral, notably PN results, with perturbative calculations for the ringdown, and physically motivated ansatzes for the plunge-merger stage. The EOB waveforms are then informed and made highly accurate by calibration to NR simulations (see, e.g., Refs. [70,72]). One key ingredient of the EOB formalism is the conservative EOB Hamiltonian. The latter, for nonspinning compact objects and in GR, is built by mapping the two-body dynamics into that of an effective body moving in a deformed Schwarzschild spacetime, whose deformation parameter is the symmetric mass ratio ν = µ/M , where µ = m A m B /M is the binary's reduced mass, m A and m B are the component masses, and M = m A + m B is the total mass [91,92]. Previous work extended the EOB Hamiltonian to ST and Einstein-Maxwell-scalar theories at 2PN and 1PN, respectively [43-45, 47, 49]. In this paper we build upon Ref. [43], and take advantage of recent progress in PN calculations in ST and ESGB theories [40,41,48], to construct an EOB Hamiltonian at 3PN order for NSs and BHs in ST and ESGB theories. This paper is organized as follows. In Sec. II, starting from the two-body 3PN Lagrangian in ST and ESGB theories, we derive, using order-reduction methods, the two-body Hamiltonian at 3PN order in the Einstein frame. In Sec. III, we construct a canonical transformation that maps the two-body Hamiltonian into the EOB Hamiltonian, including nonlocal-in-time terms due to tail effects, which are already present at 3PN order in ST and ESGB theories. More specifically, we compute such tails for generic orbits in an expansion in the orbital eccentricity parameter. In Sec. IV we specify our EOB Hamiltonian to BH binaries in the shift-symmetric ESGB model, and we calculate the orbital frequency at the ISCO. In Sec. V we summarize our main conclusions and future research directions. Various technical details are relegated to the appendixes. In Appendix A we develop a dictionary to relate quantities in the Einstein and Jordan frames. In Appendix B we list the expression of the 3PN Lagrangian. In Appendix C we discuss contact transformations of the two-body Lagrangian. In Appendix D we give the twobody Hamiltonian. Finally, in Appendix E we list the coefficients of the canonical transformations. Throughout this paper we use geometrical units (G = c = 1).

A. ST and ESGB gravity
We consider the theory described by the Einstein-frame action [34,48] where R is the Ricci scalar, g = det g µν is the metric determinant, and G = R µνρσ R µνρσ − 4R µν R µν + R 2 is the Gauss-Bonnet scalar, with R µ νρσ and R µν the Riemann and Ricci tensors, respectively. The integral of the Gauss-Bonnet scalar over a four-dimensional spacetime d 4 x √ −g G is a boundary term [94]. Matter fields Ψ are minimally coupled to the Jordan metricg µν = A 2 (φ)g µν . The dimensionless functions A and f and the constant quantity ℓ (with dimensions of length) specify the theory. We recover ST theories when either ℓ = 0 or f is a constant, and GR when moreover A (and φ) are constant.
When dealing with compact bodies, we adopt the phenomenological treatment initiated in Refs. [34,95] in ST theories, and describe them as point particles: ] is the worldline of particle A. The constant GR mass is replaced by a function m A (φ) that depends on the internal structure of body A and on the value of the scalar field at x µ A (s A ). For an explicit calculation of the mass m A (φ) of an ESGB BH, see Refs. [48,51]; see also Refs. [34,35,96] for NSs in ST theories.
From now on, we will refer to the theory with action (II.1) as "ESGB gravity," but we note that the action includes ST gravity as a special case.

B. The two-body Lagrangian at 3PN
In this paper, we focus on the conservative dynamics of compact binaries on bound orbits. When the relative orbital velocity is small and the gravitational field is weak, the motion can be studied in the PN framework. 1 To do so, the field equations of the theory (II.1) with (II.2) are solved iteratively around a flat metric g µν = η µν + δg µν and a constant scalar background φ = φ 0 + δφ, where φ 0 is imposed by the binary's cosmological environment. In particular, the functions m A (φ) and m B (φ), describing bodies A and B, can be expanded at 3PN by introducing and their counterparts for body B, where from now on the superscript 0 denotes a quantity evaluated at φ = φ 0 . The ST two-body Lagrangian was derived at 1PN by Damour and Esposito-Farèse [34], at 2PN by Mirshekari and Will [37], and at 3PN by Bernard [40,41]. It was then generalized by Julié and Berti, who derived its ESGB corrections in Ref. [48]. However, the results in Refs. [37,40,41] are presented using a different, "Jordan-frame" formulation of ST theories based on a set of Brans-Dickeinspired parameters. To recover the conventions of the present paper, we must proceed as follows: 1. Translate the parameters in Refs. [37,40,41] in terms of the quantities (II.3). The conversion is detailed in Appendix A.
3. Since alsot = A 0 t, the two-body LagrangiansL given in Refs. [37,40,41] must be rescaled as L = A 0L . 1 We denote by nPN the relative O(v 2n ) ∼ O(M/r) n corrections to Newtonian gravity, with v the system's relative orbital velocity, r the orbital separation, and M the total mass.
We denote by x A the spatial position of body A, and introduce the notations The ESGB two-body Lagrangian is then, in harmonic coordinates such that ∂ µ ( √ −gg µν ) = 0: where the contributions up to 2PN were presented in Refs. [43,44] and are recalled in Appendix B. We decompose the new 3PN contribution as where the lengthy expressions of the terms L (i) 3PN are also given in Appendix B. They depend on the logarithms ln(r/r A ) and ln(r/r B ), where r A and r B are regularization lengths that we shall eliminate later.
However, one of us noticed, while performing the conversion from the Jordan to the Einstein frame, that some terms in L (i) 3PN , originated from the results of Ref. [41], must be revised. The Einstein-frame two-body dynamics is indeed described by the action (II.1) with matter explicitly accounted for by Eq. (II.2). The associated PN Lagrangian should thus not depend on A and its derivatives at infinity. Yet, the prefactors of the first lines in Eqs. (B.2d) and (B.2e) are inversely proportional tõ α = (1 + α 0 A α 0 B )/(1 + α 2 0 ), where α 0 = (d ln A/dφ) φ0 , cf. Appendix A. This issue will be addressed in an upcoming publication [97]. For now, we note that adjustingα will not affect the structure of our results.
The ESGB correction beyond ST reads [48] ∆L ESGB 3PN = with f ′ (φ 0 ) = (df /dφ) φ0 and M = m 0 A + m 0 B . It is numerically of the same order of magnitude as a 3PN term whenever ℓ 2 f ′ (φ 0 ) ≲ M 2 . It turns out that this condition is satisfied by the nonperturbative ESGB BH solutions studied in Ref. [51].
Finally, L 3PN depends on a nonlocal-in-time "tail" contribution which we converted from the Jordan-frame expression of Ref. [41], which is driven by the acceleration of the scalar dipole This tail term is absent in GR. By the same arguments we made earlier, the tail term should be independent of A 0 = A(φ 0 ). This issue will also be addressed in Ref. [97], and for now, we note that replacing A(φ 0 ) by a different constant will not change the structure of our final results. Here, "PF" denotes the Hadamard partie finie, and we follow the conventions of Refs. [40,41,98]: given a regular function f (t) vanishing sufficiently fast at infinity and a constant s, we have (II.9) The two-body Lagrangian (II.5) depends on the theorydependent combination ℓ 2 f ′ (φ 0 ) entering Eq. (II.7), and on ten body-dependent parameters: the masses of each body and their logarithmic derivatives (II.3) at infinity. It is also useful to introduce the following quantities, ordered by the PN level at which they appear, from 0PN to 3PN: and their (A ↔ B) counterparts. The quantities (II.10d) are new to this paper, and we named the first three of them according to their (field theory) diagrammatic interpretation, as was initiated at 2PN by Damour and Esposito-Farèse in Ref. [99]. We recover the ST Lagrangian in the limit ℓ 2 f ′ (φ 0 ) = 0, such that (II.7) vanishes. We recover GR when moreover m A (φ) and m B (φ) are constants: then, Eqs. (II.3) and their B counterparts are zero, so that G AB = 1 and (II.10b)-(II.10d) all vanish.

C. The order-reduced Lagrangian
The Lagrangian (II.5) is written in harmonic coordinates, and it depends on the accelerations a A and a B of the bodies, both linearly via L 2PN , L Consider a degree of freedom q(t) described by the action I = dtL[q,q,q], where L[q,q,q] = L 0 (q,q) + ϵL 1 (q,q) + ϵ 2 [L 2 (q,q) + ℓ 2 (q,q)q] + ϵ 3 L 3 (q,q) + ℓ 3 (q,q)q +q PF The Euler-Lagrange variation of where the second equality follows from the chain rule. Now introduce the notations where H F (q,q) is the Hessian of F (q,q). Then Eq. (II.13) can be rewritten as the identitÿ which reduces toq =q F (q,q) + O(ϵ 2 ) when the Euler-Lagrange equations of F (q,q) are satisfied, δF/δq = 0. We can then insert (II.15) into (II.11) to find: is an ordinary Lagrangian depending only on q andq, obtained by replacing the accelerationq in L[q,q,q] by its on-shell expression deduced from F (q,q). 2 The third line of Eq. (II.16) is doubly zero: its O(ϵ 3 ) contributions to the equations of motion are at least linear in δF/δq, which vanishes on shell. Thus it can be discarded. As for the second line of Eq. (II. 16), it can also be eliminated via a variable change q → q + δq[q,q] with δq = O(ϵ 2 ). Indeed, the Lagrangian then transforms, by definition, as modulo an irrelevant total time derivative, and we can choose This variable change belongs to the class of contact transformations introduced by Schäfer and Damour in Refs. [100,101], which we generalized to include nonlocalin-time terms for our purpose. Now return to the two-body Lagrangian (II.5). We can replace the accelerations by their on-shell expressions deduced from F = L 0PN + L 1PN : where (a F ) i A is a function of the positions and velocities given in Appendix C. Note that it is sufficient to replace the accelerations entering the Lagrangian at 3PN level by their 0PN expressions. This procedure amounts to making an implicit 4D coordinate change via a contact transformation x A → x A + δx A resembling (II. 19), from which we deduce an identity that is useful to prove Eq. (II.16): which we also give explicitly in Appendix C. We verified that applying the contact transformation (C.4) to the two-body Lagrangian (II.5) yields a result that matches, modulo total time derivatives and doubly zero terms, the order-reduced Lagrangian (II. 20). From now on we work with the order-reduced Lagrangian L red , and thus, in a coordinate system other than harmonic.

D. The two-body Hamiltonian at 3PN
From the order-reduced Lagrangian (II.20), we can infer an ordinary Hamiltonian via the Legendre transformation: A technically useful remark is that, when deriving a Hamiltonian from a nPN Lagrangian, it is sufficient to calculate (II.23) at (n − 1)PN order, since after inversion, . We thus need the momenta at 2PN order only.
In the center-of-mass frame such that The motion being planar, we use polar coordinates (r, ϕ) with conjugate momenta p r = n · p and p ϕ = r(n × p) z . We introduce the reduced mass and mass ratios together with the following dimensionless quantities: We denote by a subscript + (respectively, −) the (anti) symmetrization of the quantities (II.10), as in, e.g., [note the factor of 2 compared to Refs. [39][40][41]]. The two-body Hamiltonian is then: where the contributions up to 2PN were first derived in Refs. [43,44] and are recalled in Appendix D. The new 3PN contributions read where the lengthy expressions ofĤ 3PN are also given in Appendix D. They depend on ln ± = ln(r A ) ± ln(r B ), wherer A = r A /M andr B = r B /M are the dimensionless regularization lengths mentioned below Eq. (II.6).
For the reason given above, the tail and ESGB contributions are equal and opposite to their Lagrangian counterparts: where the cosine of ∆ϕ = ϕ(t +τ ) − ϕ(t) follows from the order-reduced, center-of-mass frame accelerationD i = −G AB ν(α 0 A − α 0 B )n i /r 2 taken att andt +τ , and (II.29) The ESGB correction beyond ST reads Finally, to prepare for the calculations of Sec. III below when the eccentricity is nonzero, and thusr is not constant, we use the first identity in footnote 2 to get . (II.34b)

A. The EOB Hamiltonian at 3PN order
In Sec. II, we derived a two-body Hamiltonian at 3PN order and in the center-of-mass frame. We shall now use canonical transformations (r, ϕ, p r , p ϕ ) → (R, Φ, P R , P Φ ) to identify it with an EOB Hamiltonian [91,102,103] whereĤ eff is an effective Hamiltonian to be constructed.
In this paper, we write the effective Hamiltonian in the same gauge as that used in GR at 2PN, 3PN and 4PN orders [91,93,102], that iŝ which depends on three potentials. Through 3PN order they can be expanded as When restricted to 2PN order, the EOB Hamiltonian above depends on five coefficients (a 1 , a 2 , a 3 ) and (d 1 , d 2 ), which were derived in Ref. [43] in ST theories. We will recall their expressions in Sec. III D for completeness. In this paper, we introduce the remaining eight coefficients to include 3PN contributions in ST-ESGB gravity.
The effective Hamiltonian (III.2) describes the motion of a test particle with mass µ in an effective static, spherically symmetric metric [in Schwarzschild-Droste coordinates with θ = π/2] but it is now deformed by a nongeodesic 3PN potential Q which vanishes for circular orbits such thatP R = 0. The reason is twofold: 1. At 3PN order and already in GR, the two-body dynamics cannot be reduced to geodesic motion [93]. Following Damour, Jaranowski and Schäfer, we thus include a postgeodesic correction controlled by the coefficient q 1 .
2. In GR, the two-body Hamiltonian depends at 4PN order on a quadrupole-driven tail [98,104], which can yet be accounted for in a local-in-time EOB Hamiltonian by extending Q to an infinite postgeodesic series [102] together with lnR-dependent EOB potentials. Our ansatz (III.3) adapts this strategy to include the dipole-driven tail entering already at 3PN order.
Note that in practice, one must truncate Q. We choose to do so at O(P 6 R ), as otherwise it would diverge at infinity [cf. theR-dependence in Eq. (III.3c)]. This will amount to including beyond-GR tails up to sixth order in the binary's orbital eccentricity in Sec. III C. For circular orbits such thatP R = 0, we have Q = 0; in this simpler subcase than what is done here,Ĥ eff depends only on A, and the tails affect a 4 and a 4,ln only.
Finally, following Ref. [102], we find it useful to split the potentials as (III.9c) The EOB Hamiltonian can then be decomposed in a similar fashion as its two-body counterpart, whereĤ I EOB is obtained by formally setting A II , D II and Q II to zero in Eq. (III.1) while, at first order, Let us first focus on the local-in-time partĤ I of the two-body Hamiltonian, cf. Eq. (II.34a). We perform a canonical transformation such that H I is a scalar and its action only changes by a boundary term: and thus For practical reasons, we will rather use G(r, ϕ, P R , which generates a canonical transformation introduced in Refs. [43,44], We choose the ansatz: which yields coordinate changes between 1PN and 3PN levels when the positive integers i, j and k satisfy 1 ⩽ i+j +k ⩽ 3. Our ansatz does not depend on ϕ to preserve isotropy, and thus p ϕ = P Φ . Moreover, for circular orbits such that p r = P R = 0, we have Φ = ϕ. From Eqs. (III.15) we can express bothĤ I andĤ I EOB in the same mixed coordinate system (r, ϕ, P R , P Φ ). We then solve, order-by-order, the equation to fix the coefficients of the potentials (III.8) and of the generating function (III. 16). The solution is unique, and the new 3PN coefficients are: The coefficients of the canonical transformation (III.16) are given in Appendix E for completeness. We explicitly checked that in the GR limit,Ĥ I can also be identified to the 3PN (ADM) Hamiltonian of Ref. [98] via a canonical transformation whose coefficients are given in Appendix E.

C. Nonlocal-in-time contributions
Let us now turn to the nonlocal-in-time 3PN Hamilto-nianĤ II , which (we recall) readŝ , (III. 21) with ∆ϕ = ϕ(t +τ ) − ϕ(t). We wish to identify it, modulo canonical transformations, to a local-in-time, ordinary EOB counterpartĤ II EOB depending on positions and momenta only. To do so, we can Taylor expandr(t +τ ) and ϕ(t +τ ) aroundτ = 0, and treatĤ II as a local-in-time function ofr(t) and ϕ(t), and their arbitrarily high-order time derivatives. The Newtonian equations of motion can then be used to order reduce these derivatives, as we now prove.
Consider a pair of phase-space variables q(t) and p(t) described by the action where ϵ ≪ 1, and where ∆H depends on arbitrarily highorder time derivatives of q and p. The Euler-Lagrange variations of with respect to p and q yield, respectively, where we have introduced the notatioṅ As usual, the system (III.24) reduces to the Hamilton We can then insert (III.24) into (III.22) and expand at first order only in δL 0 /δp and δL 0 /δq (and their time derivatives), since higher-order contributions are doubly zero. We then find, modulo total time derivatives: where we have defined the Euler-Lagrange variations of ∆H with respect toq andṗ: Now we denote the order-reduced (over) accelerations, obtained recursively from the Hamilton equations of H 0 (q, p), byq with n ⩾ 1 and, similarly, We then have, using Eqs. (III.24) again, , which we can plug into the second term in the right-hand side of Eq. (III.26). Expanding the result at first order in δL 0 /δp and δL 0 /δq (and their time derivatives) and integrating by parts finally yields: modulo doubly zero terms and total time derivatives. The subscript "red" indicates an order-reduced quantity, as in p,ṗ 0 (q, p),p 0 (q, p), · · · .
It is now elementary to eliminate the second to fifth lines of Eq. (III.31): under a phase-space contact transformation (q, p) → (q + δq, p + δp) with δq = O(ϵ 3 ) and δp = O(ϵ 3 ), the Hamiltonian transforms as modulo an irrelevant total time derivative, and we can choose to identify δp(q, p) and δq(q, p) with the long coefficients of δL 0 /δp and δL 0 /δq in Eq. (III.31), respectively. This toy model is an adaptation of Refs. [100,101,105], which we have extended to Hamiltonians depending on arbitrarily high-order time derivatives of q and p for our purpose. Indeed, return now to the nonlocal-in-time 3PN Hamil-tonianĤ II . We shall see that there exists a set of phase-space variables other than polar in which the steps above are elegantly carried out. This is the route of Refs. [102,106], which we adapt here to the ST-ESGB case.
In polar coordinates (r, ϕ, p r , p ϕ ), the Keplerian trajectory can be parametrized by the semimajor axisâ = a/M and the eccentricity e as [107] r =â(1 − e cos η) , (III.34a) (III.34b) We sett = 0 at the periastron without loss of generality, and define the eccentric anomaly η aŝ is the mean orbital frequency. Now, we observe thatâ and e can be treated as functions of the Delaunay action-angles (L, G, l, g). 3 Indeed, it is a textbook exercise to show that, on Keplerian orbits [108], which can be inverted aŝ while the conjugate angles are the mean anomaly and argument of the periastron, respectively: (III.39b) In these canonical variables, the 0PN equations of motion are particularly simple:Ĥ 0 is the Delaunay Hamiltonian, and thus while L and G are constants. The order-reduced Hamiltonian then has the structure: (III.48) Another purpose of the Delaunay action-angles is that H II red is a perturbation of the 0PN Hamiltonian (III.40), which depends only on L in these variables. That way, the l dependence inĤ D =Ĥ 0 (L) +Ĥ II red (L, G, l) can be eliminated through a canonical transformation resembling Eq. (III.15), − e 2 (−9 lnâ + 3 ln G AB + 6γ E + 6 lnŝ + 14 ln 2) where we recall thatâ and e are the functions (III.38) of L and G.
The same steps can be applied toĤ II EOB of Eq. (III.11). It is already local-in-time and ordinary, but we rewrite it in terms ofâ and e using Eq. (III.34a) withr →R, andP R = dR/dt. 4 We then use Eq. (III.42) and invoke canonical transformations to discard l-dependent terms with the same form as in Eq. (III.47). We find: In Eqs. (III.7) and below, we split the EOB potentials into their parts I and II, which we determined in Secs. III B and III C, respectively. We recall thatP R = P R /µ, and we introduce the notation Adding the results yields: and, at 3PN order, where Eqs. (III.60a)-(III.60c) were already introduced in Ref. [43], while the remaining quantities are new to this paper. The results above are available online [109].
A few comments are in order. The potentials (III.57) are a beyond-GR extension of the 3PN results of Damour, Jaranowski and Schäfer [93], which we recover in the GR limit detailed below Eqs. (II.10). Indeed, in that limit ⟨β⟩, γ AB and the coefficients (III.58) and (III.59) all vanish. We observe that Newton's constant is now substituted by the effective gravitational coupling G AB entering u at all orders. Thus this effect can be absorbed in a redefinition of the total mass M . When truncated to 2PN order, our potentials depend on five coefficients and reproduce those of Ref. [43] in ST theories. This can be checked using Eqs. (III.58) and B = (AD) −1 .
The 3PN coefficients (III.59) are the central new results of this paper. They show that among the eight coefficients in our ansatz at 3PN level [cf. Eq. (III.3)], only five are nonzero. The contribution from the nonlocal-in-time tail is driven by the constant k tail , which enters all coefficients in Eqs. (III.59). This is necessary to include the non-GR tail beyond circular orbits, and at sixth order in the orbital eccentricity. Note that the tail is fully responsible for the unique logarithmic correctionā 4,ln and the postpost-geodesic coefficientq 2 . As for the ESGB corrections beyond ST, they are driven by k ESGB , and their inclusion is particularly simple: they only enter in δā 4 .
Contrary to the two-body HamiltonianĤ, the coefficients (III.59) do not depend on ln ± = ln(r A ) ± ln(r B ), wherer A andr B are the regularization lengths mentioned in Sec. II D. As expected, they have been reabsorbed in a canonical transformation at 3PN level (see Appendix E).
Finally, at the end of Sec. II D we split the two-body Hamiltonian into its local-in-time and nonlocal-in-time partsĤ I andĤ II by introducing an arbitrary constantŝ which propagated in both A I and A II , cf. Eqs. (III.19a) and (III.54a). As expected,ŝ cancels out from our final EOB potentials.

A. Hairy BH binaries
The coefficients of the potentials (III.57), we recall, are built out of the theory-dependent product ℓ 2 f ′ (φ 0 ) and of ten body-dependent parameters: the values of the masses m A (φ) and m B (φ) and their logarithmic derivatives (II.3) evaluated at infinity (i.e., at φ = φ 0 ). Now, these quantities can be calculated once the theory and the bodies are specified. In ST theories, they were derived numerically for NSs and their scalarized counterparts (e.g., see Refs. [35,110] and references therein). They were also calculated for BHs in ESGB models, both analytically in the small-ℓ limit, and numerically for nonperturbative solutions such as scalarized BHs (cf. Refs. [48,51]).
Let us complete this paper with an explicit illustration. Consider a BH in the shift-symmetric theory f (φ) = 2φ and A(φ) = 1. For simplicity, here we consider only terms at leading order in ℓ, such that [48] α 0 , and can thus be neglected in what follows. The BH is fully described by the values of ℓ and m 0 A . In the GR limit ℓ/m 0 A = 0, the quantities above all vanish, because then the BH reduces to Schwarzschild and its mass m A (φ) is a constant (cf. Ref. [48]). We find it useful to introduce the dimensionless ratiol . This means that body A reduces to a Schwarzschild spacetime with constant scalar field, cf. Eq. (IV.1) and below. We recover the conservative sector of the extreme mass-ratio analysis of Ref. [111].

B. Orbital frequency at the ISCO
We can now evaluate the beyond-GR modifications to the dynamics, focusing on circular orbits for simplicity. Consider the motion described by the effective Hamilto-nianĤ eff given in Eq. (III.2). It does not depend ont nor on Φ, and thusĤ are constants of motion. WhenP R = 0, we have from the system above that while the circularity of the orbit also requires dP R /dt = 0, that is The ISCO is characterized by a third (inflection point) condition, Hence E and j = J/G AB relate to u [cf. Eq. (III.56)] as where the primes denote derivatives with respect to u, while u ISCO is the outermost root of Let us turn to the EOB HamiltonianĤ EOB given in Eq. (III.1). The associated Hamilton equations define a resummed two-body dynamics. In this paper, we focus on the dimensionless orbital frequencyΩ = M Ω = dΦ/dt, which readŝ where E(u) and j(u) are given by Eqs. (IV.11) on circular orbits. The orbital frequency, which we shall evaluate at the ISCO, is thus fully fixed by the effective potential A [91,112]. We follow Refs. [93,113] and resum our 3PN result (III.57a) by means of the (1, 3)-Padé approximant (IV.14) This ensures that A P (u) has a simple zero (by construction) and the presence of an ISCO, by continuity with the Schwarzschild metric recovered in the GR, test-mass limit. The Padé resummation was adopted in several studies that calibrated GR-EOB waveforms to NR [113,114]. For further discussions on the effects of the Padé resummation in the ST case at 2PN, see Ref. [43]. Figure 1 shows the ISCO location u ISCO and dimensionless frequency G ABΩ of a binary BH system in the shift-symmetric ESGB model discussed in Sec. IV A. The beyond-GR coefficients of the potential A are thus the functions of ν andl = ℓ/µ given in Eqs. (IV.3)-(IV.6), which we truncated at the leading order inl given there. We recover GR whenl = 0, and we consider four symmetric mass ratio values, ν = {0, 0.1, 0.2, 1/4}.
When ν = 0, we find that u ISCO = 1/6 and G ABΩ = 6 −3/2 reduce to their Schwarzschild values for alll, consistently with the extreme mass-ratio limit described at the end of Sec. IV A. However, when ν ̸ = 0, both u ISCO and G ABΩ increase withl. In particular, the slope (or "sensitivity") of the ISCO frequency is maximal when ν = 1/4: For such equal-mass binaries (with µ = M/4), the relative modification to the GR ISCO frequency, (G ABΩ ) ISCO /(G ABΩ )l =0 ISCO − 1, then reaches the percent level when ℓ/M = 0.132. For comparison, Ref. [115] obtained one of the most stringent constraints to this day in shift-symmetric ESGB gravity, ℓ/M < 0.344, from the BH-NS system GW200115 with total mass M = 7.1M Sun . (Note that we translate between our conventions and those of Ref. [115] by setting φ = √ 4πϕ and ℓ 2 = 2 √ 4πα GB .) The ISCO analysis above motivates the obtention of full EOB waveforms, including the dissipative sector [34,38,50,80], to be confronted to GW signals. This issue will be addressed in future work.

V. CONCLUSIONS
In this paper we have extended the work of Refs. [43,44] and built an EOB Hamiltonian in ST and ESGB gravity at 3PN order. Our main new results are:  (Fig. 1).
It is important that the EOB framework can be extended beyond GR. Here, we have reduced the 3PN dynamics to the (nongeodesic) motion in a modification of the GR EOB metric, and accounted for the beyond-GR tail effects by adapting the 4PN methods of Ref. [102].
The EOB framework is also suitable to include other modified theories of gravity, such as Einstein-Maxwellscalar models at 1PN [45,47,49]. Our work can thus be regarded as another step toward the development of a parametrized EOB framework, by providing a "dictionary" between modified gravity theories and the values of the coefficients of the effective potentials (III.57). In the future, the tools and methods we developed in this paper could be applied to other models, such as disformal ST, massive gravity, or Horndeski theories (which also predict hairy BHs [116]).
We have focused here on the conservative part of the dynamics. The corresponding EOB radiation-reaction force, to be inferred from already available energy fluxes [34,38,50,80], and gravitational waveforms, will be the topic of future work. Since NSs and BHs are, in general, spinning, it will also be important to extend the present work to include spin effects. For the PN analysis of spin-orbit effects in ST gravity some work has been done in Ref. [117]. As a first step, the beyond-GR EOB Hamiltonian derived here could be included in the stateof-the-art spinning EOB Hamiltonians in GR (see, e.g., Refs. [114,118] and references therein), and then used to generate beyond-GR inspiral waveforms.
The EOB approach uses a resummation of the twobody dynamics that can be extended through the plunge of the two BHs, after which the waveform is matched to the merger-ringdown signal. The latter should make use of the quasinormal mode spectrum of ESGB BHs, which has been computed up to second order in a slow-rotation expansion [89,90]. The Padé-resummed spectrum of Kerr BHs computed at the same order in the slow-rotation approximation is typically accurate at the percent level when evaluated at the dimensionless spins ∼ 0.7 of interest for LIGO-Virgo-KAGRA observations [90]. Therefore it is reasonable to assume that deviations induced by beyond-GR terms could be testable at the same (percent) level of accuracy. The quasinormal mode spectrum could be included in the EOB model using the parametrized spin expansion coefficient (ParSpec) framework [26], which has been used to perform theory-specific tests of GR with ringdown signals using the pyRing code in Ref. [119], and with EOB waveforms in Refs. [120].
For the case of binary BHs, once the EOB waveforms are completed with physically motivated ansatzes for the merger-ringdown in ESGB gravity, they could be compared with and informed by NR simulations (see, e.g., Ref. [66]). Developing precise and complete EOB-NR waveform models is crucial to obtain new experimental bounds on ESGB models, and more generally, on wider classes of modified gravity theories in the future.
Note added: While this project was nearing completion, we became aware of an independent effort that recently appeared on the arXiv [121]. Their work focuses on the computation of the 3PN EOB Hamiltonian in ST theories, and restricts the inclusion of tail effects to circular orbits. This limit amounts to setting formally k ESGB = 0, and k tail = 0 in Eqs. (III.59b)-(III.59d).
After both works appeared on the arXiv, we compared the results in the limit given above. We found that: 5 1. Equation (5.16) of Ref. [121] still differs from Eq. (III.59c) by an overall minus sign; 2. The term proportional to ⟨δ⟩/α AB in Eq. (5.14) of Ref. [121] differs from that of Eq. (III.59a) by a factorγ AB . Note that α AB =α in our conventions.
We explicitly checked that reexpanding our H EOB at 3PN order yields Eqs. (5.4)-(5.5) of Ref. [41] on circular orbits. By contrast, we find that it does not if we replace Eq. (III.59a) by Eq. (5.14) of Ref. [121]. After our work appeared on the arXiv, Ref. [121] was extended in Ref. [122] to include the tail effects up to O(e 4 ). This limit amounts to setting formally k ESGB = 0, and k tail = 0 in Eq. (III.59d), since we recall that our work includes the tail effects up to O(e 6 ). The tail contributions (4.12)-(4.17) in Ref. [122] all differ from ours by an overall factor A 2 0 . The latter indeed enters Eq. (II.29), and it originates from the translation of the Jordan-frame Lagrangian (A3) of Ref. [

ACKNOWLEDGMENTS
We thank Laura Bernard for sharing a Mathematica notebook containing the 3PN Lagrangian of Ref. [41]. We are also grateful to Thibault Damour  In Refs. [37,40,41], the ST two-body Lagrangian was computed up to 3PN order by adopting the Jordan-frame formulation of the theory (we use tildes for clarity): where (∂ϕ) 2 =g µν ∂ µ ϕ∂ ν ϕ, and where ω(ϕ) is a function defining the theory. As for compact bodies, they were described by performing the substitution In the present paper, we describe ST theories by the Einstein-frame action, which (setting to zero the GB coupling) reads: where (∂φ) 2 = g µν ∂ µ φ∂ ν φ, and where we account for compact bodies by the substitution where φ(ϕ) is obtained by inverting A(φ) = 1/ √ ϕ. Let us also introduce the notation: where the subscript 0 denotes a quantity evaluated at infinity, φ(ϕ 0 ) = φ 0 . The quantities above can be obtained by inserting Eq. (A.5c) into Eqs. (II.3) and taking the limitm A (ϕ) = const. In this limit, body A is said to have negligible self-gravity, and its motion reduces to geodesics ofg µν , cf. Eq. (A.2).