UvA-DARE (Digital Academic Repository) Post-Newtonian gravitational and scalar waves in scalar-Gauss-Bonnet gravity

Gravitational waves emitted by black hole binary inspiral and mergers enable unprecedented strong-field tests of gravity, requiring accurate theoretical modeling of the expected signals in extensions of general relativity. In this paper we model the gravitational wave emission of inspiralling binaries in scalar Gauss–Bonnet gravity theories. Going beyond the weak-coupling approximation, we derive the gravitational waveform to relative first post-Newtonianorder beyond the quadrupole approximation and calculate new contributions from nonlinear curvature terms. We also compute the scalar waveform to relative 0.5PN order beyond the leading − 0.5PN order terms. We quantify the effect of these terms and provide ready-to-implement gravitational wave and scalar waveforms as well as the Fourier domain

The GW measurements of BH binaries rely on accurate template waveforms that cover the inspiral, merger and ringdown of compact binaries.Significant recent progress has been made on constructing accurate templates in general relativity (GR).However, testing the underlying theory of gravity and searching for signatures of quantum gravity requires theory-specific waveforms beyond GR.Therefore, GW-based tests of gravity have focused mostly on parameterized or null tests against GR [7,9,10,12,13].These tests are performed only on single coefficients associated to a scaling with frequency [9,10,14,15] and thus, the interpretation of theoretical constraints remains limited [7].Given the plethora of proposed gravity theories beyond GR, it is not feasible to compute accurate template waveforms for all of them.This makes it important to consider well-motivated classes of theories that capture broadly applicable features for the GW modeling.
One of the most compelling beyond-GR class of theories is scalar Gauss-Bonnet (sGB) gravity, which extends GR by a dynamical scalar field non-minimally coupled to the Gauss-Bonnet (GB) invariant R 2 GB = R μνσρ R μνσρ − 4R μν R μν + R 2 , with a coupling parameter α. sGB is particularly well motivated for several reasons (i) it is inspired by the low energy limit of quantum gravity paradigms such as string theory after compactification [16][17][18]; (ii) it can be derived from Lovelock gravity, the most general theory of gravity in D spacetime dimensions that yields second order field equations, via a dimensional reduction [19,20]; (iii) it corresponds to the lowest order in a series expansion in curvature terms and thus, it is representative of a more general class of higher derivative theories [21]; and (iv) its field equations contain at most second derivatives of the metric, so sGB can potentially be made mathematically well-posed (as shown for weak couplings to the GB invariant in [22][23][24]).
Type I sGB gravity has been tested extensively against observations of low-mass x-ray binaries [50] and a Bayesian parameter estimation for GW detections [51].The most recent observational bounds of √ α 1.7 km were obtained by analysing signals of LIGO/VIRGO's first two GW catalogs [52,53].Type II sGB gravity has, so far, remained unconstrained from observations.
To-date, modeling the dynamics of compact binaries in sGB gravity has mainly been restricted to the small coupling regime.This approximation has enabled the first numerical relativity (NR) simulations of BH binaries with an effective-field theoretical treatment [54,55] or the decoupling limit [41].The numerical study of fully non-linear field equations for general couplings has recently been advanced in the direction of initial data construction [56], and numerical simulations [57] through a modified generalized harmonic formulation of the evolution equations [23].This has also enabled the recent study of spontaneous BH scalarization [47].The ringdown spectrum has been studied extensively in references [58,59].However, the early inspiral waveform, typically first constructed with post-Newtonian (PN) methods, has only been modeled at leading (Newtonian) order [60].At this order, the dynamics and the radiation are affected only by the scalar field, and consequently the waveforms are missing the effect of the curvature nonlinearities.Recently, the Lagrangian for the matter dynamics has been computed to 1PN order, where the nonlinearities first appear [61].
In this work we compute, for the first time, analytical waveforms for the inspiral stage of the binary evolution with the effect of nonlinear higher curvature corrections.We use the PN expansion with the direct integration of the relaxed Einstein equations (DIRE) approach, originally developed by Epstein and Wagoner [62] and extended by Will, Wiseman and Pati [63][64][65].The PN assumption of a gravitationally bound source in the weak-field, low-velocity regime Gm/rc 2 ≈ v 2 /c 2  1, where m, r, and v are the characteristic mass, size, and velocity of the source, enables asymptotic expansions of the solutions to the field equations with 1/c 2 treated as a formal small expansion parameter.Each factor of c −2 then corresponds to one PN order.In the DIRE approach to PN calculations, the field equations are re-expressed in a relaxed form [66], namely as an inhomogeneous, sourced wave equation for field perturbations, accompanied by a harmonic gauge condition that is imposed after solving the field equations.The equations are solved in two space-time regions, (i) the near-zone, where the separation between source and field point is less than the characteristic wavelength of the GWs; and (ii) the far zone at larger distances, see figure 1.In both cases, the solutions are computed as retarded integrals over the past null cone of the field point.This integration method has been shown, for example, to provide non-divergent results for gravitational waveforms to 2PN order in agreement with results obtained using the multipolar post-Minkowskian formalism [67].We note that this method has also been applied to massless scalar-tensor (ST) gravity theories-which admits scalarized neuron star solutions [68]-providing the equation of motion for compact binaries to 3PN order [69,70], GWs to 2PN [71,72], and scalar waves to 1.5PN [73] beyond the leading order.Specifically, in this paper, we calculate scalar and tensor waves of sGB gravity to half and one relative PN-order, respectively.This is the methods paper to our shorter companion paper [74] in which we solely concentrated on the application of our results and assessed the detectability of GB deviations in the Fourier-domain GW phase.In this paper, we present the explicit computations to produce the GW waveforms and phasing as well as the scalar waveform.Our methods and results are applicable to arbitrary sGB coupling strengths that lie within the theoretical and observational bounds [18,28], as well as general couplings that have remained unconstrained.We also calculate the equations of motion of the sources to 1PN order, based on the results of reference [61].Having calculated the waveforms, we then focus on quasi-circular inspirals and calculate the energy flux carried off by radiation, as well as the GW phasing to which both the measurements and parametrized test of gravity such as [75] are very sensitive.Our results include higher order strong-field effects than previously computed, which may mimic biases in fundamental source parameters of BH binaries when analysing with GR-only GW waveforms.We also show that the effect of the GB scalar is distinct from the scalar field of an ST theory due to the presence of explicit 1PN GB coupling dependent terms.This has consequences for interpreting GW signals from BH-neutron star binaries when analyzed by beyond-GR theories [76].Further, we discuss the importance of distinguishing different perspectives on the sGB action (i.e. as a phenomenological extension of GR or as resulting from string theory) when interpreting GW constraints as they correspond to different physical frames.
This paper is organized as follows.In section 2, we present the field equations for tensor and scalar perturbations and define the matter action for self-gravitating binaries.We cast the field equations into relaxed form in section 3, and outline the procedure for finding the solution at large distances from the source.In section 4, we calculate the near-zone contribution to gravitational and scalar waves for binary systems.We defer the calculation of far-zone contribution to appendix C because, as we show, it does not contribute to the final waveform at 1PN order.In section 5 we derive the equations of motion using the two-body Lagrangian from reference [77], and transform the waveforms to the center-of-mass (CM) frame.In section 6, we solve for the gravitational and scalar energy loss rate, from which we find the time domain and frequency domain waveforms in sections 7 and 8, respectively.In section 7, we also compare our results for the scalar waveform to the NR waveforms of reference [54] computed in the small coupling approximation.Finally, in section 9, we present the discussion and conclusions.
In this paper we use the standard notation for symmetrized and antisymmetrized indices, namely x (i y j) denotes symmetrization and x [i y j] anti-symmetrization.We use a multiindex notation for products of vector components: x i jk ≡ x i x j x k , and a capital letter superscript denotes a product of that dimensionality: x L ≡ x k 1 x k 2 . . .x k l .A subscript preceded by a comma such as A, μ stands for the partial derivative ∂ μ A, while ∇ μ A indicates the covariant derivative associated to the metric.We label the two bodies in a binary system by A, B.

Gravitational action and the field equations
The gravitational action of sGB theory, in the presence of matter, is written in the Einstein frame as where R is the Ricci scalar on the four-dimensional manifold M with the metric g μν in the Einstein frame.S m is the matter action with Ψ m denoting the matter fields with a generic nonminimal coupling to the metric.The fundamental coupling constant α has dimensions of length squared, and f (φ) is the dimensionless coupling function.

Different viewpoints of the action and their consequences
We note that action (1) may be interpreted in two distinct ways: (i) one may postulate that the action describes quadratic gravity corrections to GR in the physical (i.e.Jordan) frame, where matter fields are minimally coupled to the metric and results directly correspond to observables.This viewpoint is typically adopted in the literature on testing gravity (see e.g.references [32,60]) and corresponds to setting the conformal factor A(φ) = 1 in equation ( 1).(ii) If the action, instead, is regarded as the low-energy effective action of string theories, the conformal factor is nontrivial.From this perspective, equation (1) denotes the action in the Einstein frame, and one must transform any results to the Jordan (or 'string') frame in which observables are measured.For example, the bosonic sector of heterotic string theory is characterized by f (φ) = e −2φ for which the conformal factor is A(φ) = e φ (see e.g.reference [78]).Both of these interpretations of the action become reconciled in the weak-field regime, where A(φ) = e φ = 1 + O(φ).Nevertheless, it is important to remain aware of these assumptions when interpreting observational bounds on higher-curvature effects and their possible implications (or bounds) on underlying quantum gravity theories.In particular, a non-minimal coupling factor in Einstein frame would lead to additional polarization modes when converted to the Jordan frame.Ignoring such couplings would lead to the conclusion that a breathing scalar mode is absent in the sGB theory, in contrast to, for example, ST theories.

Matter action
To describe the action of compact self-gravitating binaries, one has to take the internal gravity of each object into consideration.Based on the approach pioneered by Eardley [79] for ST theories [79,80], and later generalized to Einstein-Maxwell-dilaton theories [81], this can be done by describing the compact objects as point particles but with scalar-dependent total masses.Such a skeletonized matter action for binary systems is thus implicitly dependent on the scalar field and is given by the ansatz Here dx μ A is the world line of particle A and the internal self-gravity of the compact objects is incorporated through the scalar-dependent mass M A (φ).
Note that the point particle description of equation ( 4) does not depend on any gradients of scalar field, which results in neglecting finite size effects (e.g.tidal effects).Inclusion of such terms and calculation of the corresponding scalar-induced tidal effects is subject of future work.
Within the PN approximation, the expression for M A (φ) can be parametrized by its expansion about the background scalar field φ 0 : with δφ = φ − φ 0 , and the strong-field parameters Also, m A is the asymptotic value of the particle's mass in the Einstein frame.The parameter α 0 A is called the scalar charge and measures the coupling of the physically measurable mass M A to the background scalar field.Within the small-coupling approximation, the explicit form of α 0 A corresponding to static, spherically symmetric BHs has been derived to fourth order in the coupling parameter in reference [61].For example, to first order, the scalar charge is

Relaxed-field equations and wave generation
To solve the field equations in the weak-field limit, we use the DIRE approach adapted to modified theories of gravity.In order to do so, it is common to introduce the tensor-density g ab = √ −gg ab called gothic metric.By defining it can be shown that the following is an identity, valid for any space-time: where t αβ LL is the Landau-Lifshitz energy-momentum pseudo-tensor as given in reference [66].
To incorporate sGB gravity into this framework, we rewrite the expression for G μν given by equation (3b) in terms of the gothic metric.In order to analyze the behavior of fields outside the sources and to look for the generated tensor and scalar radiation, we expand the gothic metric around Minkowski space-time and the scalar field around the background value φ 0 .The perturbation variables are thus defined through g μν = η μν + h μν , and φ = φ 0 + δφ.(10) By choosing the Harmonic gauge defined as ∂ ν g μν = 0, equation ( 9) can be written as a wave equation for the perturbations, with sources where Λ αβ GB captures the scalar field and the nonlinear GB contributions, with * R * cαβd the gauge fixed, dual Riemann tensor written in terms of gothic variables.The explicit expression for Λ GB in terms of the metric and scalar perturbations is given in appendix A.3.The scalar field equation, equation (3a), is itself a wave equation.In terms of the gothic metric we have with R2 GB being the gauge-fixed, GB scalar written in terms of gothic variables.For its expression in terms of metric perturbation see appendix A.3.
Formal solutions to the wave equations ( 11) and ( 12) can be computed in all regions of space-time by using the appropriate retarded Green's function, giving The integrals here are over the past null cone of the field point x = (ct, x).
In principle, the explicit solution of equation ( 13) depends on the position of the field point x relative to the source x , i.e. the size of the quantity x − x = R relative to other scales in the system.Defining the characteristic size of the source to be S, the near-zone is defined as the region with R < R, where R ∼ S/v is the characteristic wavelength of GWs from the system.The region outside the near-zone is called far-zone (see figure 1).Given the position of the field point x, one can evaluate each integral in two separate pieces: an integral over the source point x in the near-zone and another over the source point in the far-zone.This leads to four different integrals in total, as described in detail in references [63,64].
We are particularly interested in the two calculations based on a far-zone field point.Assuming weak-field and low-velocity approximations, we perturbatively expand the nonlinear terms in μ μν , and its scalar analogue μ s , using the formal PN expansion parameter 1/c 2 .We keep terms up to the relative first PN order.Our matching conditions between the near and far zone are chosen such that terms depending on the boundary radius that separates the two zones cancel, i.e. the final answer is independent of this parameter.This property is generally expected and is explicitly shown to be correct within GR [63].
In the next section, we first find the near-zone contribution to the waveforms, by defining the tensor and scalar multipole moments and explicitly calculating them for binary systems.We outline the general formalism for evaluating the far-zone contribution to waveforms in appendix C, and show that the results are of higher PN order than considered in this paper.In both cases, we would only be concerned with the results in a subset of the far-zone named the far-away zone, where the GW detectors are located.We hold off on presenting the final expressions of waveforms until section 5.

General structure of near-zone calculation
Having the field point in the far-away zone, with R R, it can be shown [63] that the near-zone contribution to h i j , named h i j N , is given as where nk i is the unit normal vector pointing from the source to the detector.The tensor multipole moments I i jk 1 ...k m EW are known as Epstein-Wagoner (EW) moments and are functions of retarded time τ = t − R/c.They are given by with the so-called surface moments being These two-and three-index surface moments are the results of using the conservation law ∂ α μ αβ = 0 in the derivation of equation ( 14) (see reference [63]).As before, R is the radius of the boundary region ∂M of the near-zone, and we have used the fact that the surface element is dS k = R 2 dΩ 2 nk at this boundary region.
For computing GWs to relative 1PN order higher than the quadrupole, it is sufficient to consider up to the four-index EW moments and to only look for transverse-traceless (TT) parts of the spatial tensor, such that In the above equation, the TT projection operator acting on a tensor A kl is such that with the properties P ii = 2, P i j P i j = 2, and P i j P ik = P jk .

Formal structure of near-zone source
To compute EW moments requires expanding the components of the source μ αβ to 1PN order.The 1PN expansions will depend on the solution of equation ( 13) found in the near-zone.To 1PN order, these potentials are computed and can be found in reference [61].
For the expansion of the source terms, we begin by defining a simplifying notation for the different components of the fields, following reference [64], Using the same strategies used in GR calculations, we find that the weak-field expansion of the source μ αβ to 1PN order (see appendix A.3 for the expansion of the GB terms), is As pointed out, for a binary system, the expressions for the fields in the near-zone are given in reference [61].In the above expressions we only need to substitute for N and Φ to the lowest PN order, which is given by where

Evaluation of two-body Epstein-Wagoner moments
As can be seen from equations ( 15) and ( 16), the EW moments are integrals over a sphere of radius R, about the CM of the system, with the variables entering the integrands to be evaluated at retarded time τ .We repeatedly integrate by parts using the identity (see reference [63], equation (4.1)) Moreover, as mentioned earlier, we are only interested in the physically measurable, TT components of the far-zone tensor perturbation.We therefore make frequent use of the identities below, which follow from the definition of the projection operator P i j in equation ( 18) where the indices i and j apply to the final components of the waveform (and not the integrands), and F denotes a general term.The general method for evaluating the field integrals of EW moments is explained in detail in reference [63] and can be summarized to the following steps: (i) use integration by parts to leave one potential undifferentiated using equation (22) and keep track of the resulting surface terms, (ii) change integration variables in such a way to put the center of the differentiated potentials at the origin, (iii) expand the undifferentiated potential in spherical harmonics, (iv) express all unit vector products in terms of symmetric-trace-free (STF) unit vectors and use the relations between spherical harmonics and STF unit vectors (see appendix B) to integrate over dΩ, (v) use the integration formula (equation ( 112)) for f ( 1 r ) to integrate over dr with r = r A − r B .
In the following calculations, we only report the calculations that are new to the sGB theory, and we refer the reader to reference [63] for a detailed description of integration methods, as well as the calculation of GR contributions.In appendix B, we calculate two of the key integrals and explain their main steps that are then used throughout the rest of the calculations.
We note that the specific form of the near-zone potentials as given in equation ( 21) is such that Newtonian-order potential of body A. In the following calculations, we will use this relation to substitute for Φ so that the source terms simplify and that the functional form of the integrals mimic the usual GR terms.

Two-index moment I ij
EW .We split the calculation into GB coupling dependent and independent terms, such that I i j EW = I i j + I i j GB .The calculation of I i j is similar to that of GR at 1PN order.The only difference is the additional scalar charge dependent term in the energy-momentum tensor of matter (α 0 A Φ) as well as the (∇Φ) 2 contribution, which follow the same calculation as the GR ones, resulting in where r = r A − r B is the relative separation between the two masses.
For I i j GB to 1PN order, substituting GB terms of equation (20a) in equation (15) gives where we have omitted the coupling factor for simplicity.We use integration-by-parts several times on this expression and use equation ( 21) to simplify the potentials.Ignoring constants and scalar-charge factors in front, we find where the Newtonian potential U = A U A = A Gm A /r A c 2 is used to simplify the expressions.It can be shown that the resulting surface integrals vanish as they either do not have R independent integrands or involve integration of a delta function δ 3 (x A − x B ) that is zero on the boundary of the near-zone.For the evaluation of the last volume integral, we use the fact that As the calculations are being done in the inspiral regime, where x A = x B by definition, this term is zero and thus, we are left with The first terms involves ∇ 2 U and can be readily evaluated.For the second term, we further apply integration by parts and get It can be shown again that the surface term vanishes and that the second term is similar to the first term of equation (27).The new contribution is the volume integral in the last line that involves integration of a U ,m U ,l term for which we can again use integration-by-parts to leave one potential undifferentiated.This gives The surface integral does not have a R-independent term and thus does not contribute.The remaining volume integral can be evaluated using methods within GR and can be found in appendix B.
Putting together the non-vanishing contributions, the final expression for the GB coupling dependent two-index moment is Regarding the two-index surface moment I i j EW(surf ) , we can show that the coupling independent terms do not have any TT contribution at 1PN order.The contributions from terms depending on the GB coupling also vanish by taking surface integration, as they do not give R-independent integrands.
Dropping the last term of equation ( 30) that does not have a TT contribution and adding equation ( 24), the final expression for I i j EW becomes An important point to note is the scaling of the GB terms.Since the terms enter at O(c −2 ) in the PN expansion, they are a 1PN effect, yet, due to the dimensionless ratio α/r 2 , the GB contribution is suppressed at large separation compared to the other 1PN terms.This is not surprising as it encapsulates a different physical effect.This is similar to the case of tidal interactions in GR, where Newtonian tides at O(c 0 ) have the same scaling with separation r as 5PN point-mass terms [82].Calculations of tidal effects in compact binaries involve a double expansion in PN and finite-size corrections [83].Likewise, the GB term in equation ( 30) is the leading-order gravitational effect of the higher curvature corrections.However, since the GB effect first appears only at 1PN order, we can here keep the full dependence on the GB coupling without requiring an explicit double expansion, nor any assumptions on the coupling strength.

Three-index moment I ijk
EW .In order to obtain the waveform to relative 1PN order, it is sufficient to calculate I i jk EW to 0.5PN order.As there are no GB dependent terms at this order (see equation (20b)), the three-index moment is the same as in GR Also, the three-index surface moment I i jk EW(surf) is shown to be zero in GR to 1PN order.

Four-index moment I ijkl
EW .We again split the calculation into coupling dependent and independent terms, namely I i jkl EW = I i jkl + I i jkl GB .The moment I i jkl has a GR part and a contribution from the scalar field that is similar to the GR one [63] up to a factor difference.Including the factors, these terms together give where the terms proportional to δ i j are dropped as they do not produce TT components.For I i jkl GB we have Using equation (21) we see that the terms in the first parentheses vanish.For the terms in the second parentheses, we can directly insert the Laplacian of equation ( 21) to remove the integrals.Ignoring the δ i j dependent terms that do not give TT components we find

Scalar multipole moments and waveform
The scalar field perturbation, solving its field equation (13b) in the near zone and with the source term given in equation ( 12), can be written as with the scalar multipole moments defined as It is important at this stage to discuss the counting of PN orders.By comparing equation (36) to equation ( 14), we see that the lowest order term in the tensor wave, containing the quadrupole moment, is of order O(c −2 ) higher than that of the scalar wave due to the double time derivative.Thus, by assuming the leading (quadrupole) part of the tensor radiation to be relative 0PN order, the scalar monopole and dipole moments will be relatively −1PN and −0.5PN, respectively.This implies that finding the scalar radiation to 1PN requires the scalar moments to be calculated to relative 2PN order.
As we have expanded R 2 GB (see appendix A.3, equation ( 104)) to leading order in the weakfield limit, we are limited to having μ s to O(c −3 ) at most, and therefore, our final expression for the scalar waveform will be at most 0.5PN order.
Using the scalings of terms in equation ( 19), it can be easily shown that the PN expansion of R 2 GB does not contain 0.5PN and 1.5PN terms.The 1PN term, put together with the expansion of the matter source, gives the scalar source With this result, we can now calculate the scalar multipole moments of equation (37).We define I L s = I L s,c + I L s,GB in order to separate the calculation of GB coupling dependent and independent terms again.As the 0.5PN expansion of GB source is zero, the scalar quadrupole moment at this order does not have a GB contribution.Also, integration-by-parts shows that the 1PN GB contribution to the scalar monopole moment vanishes.The 1.5PN contribution to the scalar monopole is also zero as we mentioned that there is no 1.5PNGB source term.This leaves only the 1PN scalar dipole calculation.After integrating-by-parts and neglecting the surface terms that do not contribute, we obtain The non-vanishing part of the above equation involves the integral (∇N) 2 d 3 x.As already discussed (see discussion below equation ( 30)), this term is also zero and thus leaves us with no explicit GB corrections to the multipole moments.
Thus, the only remaining contributions to I L s,c are terms depending on the matter sources.It is straightforward to see that: where we have defined

Center of mass and relative acceleration to 1 post-Newtonian order
In the previous subsections, we explicitly derived the EW moments and the scalar multipoles for a binary system.In this section, we calculate their corresponding waveforms in the CM frame.
The conservative dynamics of a two-body system to 1PN order is governed by the Lagrangian derived in reference [61], where the authors constructed the Fokker Lagrangian using the near-zone solution to the field equations (equation ( 3)), finding As before, r = r = x A − x B so that nAB = r/r and in analogous to ST theories [80], we define the binary parameters The parameters γ and β are generalizations of 1PN Eddington parameters, and ᾱ appears as a re-scaling factor for the gravitational constant G.It is useful to also introduce the combinations We see that the Lagrangian L AB decomposes into a piece that is structurally equivalent to the two-body Lagrangian in massless ST theories [80] plus a distinct GB-coupling dependent contribution.The effects of the scalar field are thus entirely contained in the charge-dependent binary parameters of equation (42).Currently, the expansion of the Lagrangian in massless ST theory is known up to 3PN order [69].Up to 2PN order terms, the gravitational constant G always appears in combination with ST parameter that is analogous to ᾱ, which implies that the effective gravitational coupling in these theories is Gᾱ to this order.However in sGB gravity, the GB higher curvature terms break this rescaling at 1PN order, as made explicit in the last line of equation (41).Thus, starting at the 1PN relative order, the effect of the scalar on the two-body dynamics in the two theories can in fact be distinguished.
Our final aim in this section is to express the waveforms in terms of relative variables, by transforming to the CM frame with X i CM = 0 defined as By computing the above integral with the same methods as in section 4.3, we find that, to 1PN order, the coordinates of each body in the CM frame are where The GB coupling dependent term is given by Taking a time derivative of equation ( 45), we obtain the velocities with the GB coupling dependent contribution being We substitute these relations in the Euler-Lagrange equations corresponding to equation (41), and find the 1PN relative matter equation of motion From the Lagrangian, we can also read the conserved binding energy of the system to 1PN order, In the CM frame, we have

Final expressions for the waveform
By substituting the expressions for EW moments (see equations (31), (32), and ( 35)) in equation (17), and simplifying further using equations ( 45) and (47), we find to O(c −3 ) that with ) being the orbital binding energy per reduced mass at Newtonian order.
Taking time derivatives and perturbatively using the relative acceleration when needed (see equation ( 49)), we end up with the final expression of the near-zone contribution to the gravitational waveform.As shown in appendix C, the far-zone contribution is beyond 1PN order, and thus the overall expression for the gravitational waveform becomes where P is a book keeping parameter with its superscripts denoting the PN order of each expression.
For the scalar waveform (36), we first find the various contributions Φ l , including up to 0.5PN terms.By rewriting scalar multipoles (equation ( 40)) in the CM frame using equations ( 45)- (47), and taking time derivatives, we find As expected, the scalar dipole moment does not vanish by choosing CM coordinates.This is a direct consequence of sGB theories violating the strong equivalence principle.Also, up to the order we are considering, the terms that explicitly depend on GB coupling only lead to scalar dipole radiation.Higher scalar multipoles, such as the quadrupole radiation, appear at higher PN orders.These modes are the dominant contribution for equalmass binaries, for which the dipole radiation is suppressed.We will come back to this point in section 7.4, where we compare the PN scalar waveform with scalar waves from NR simulations.
Next, we consider the scalar waves, arranging terms together based on their PN order.The final expression for the scalar waveform is given by For both the scalar and tensor waveforms, we obtain similar results in structure to those of ST gravity [71,73] (i.e. through redefined parameters in ( 42)), apart from new terms that explicitly depend on the GB coupling parameter.This feature allows to potentially distinguish the two theories in BH-neutron star binaries, where in ST theories only neutron stars can develop scalar hair.For GWs up to 0.5PN order, the only difference with respect to GR is the presence of the factor ᾱ multiplying the gravitational constant.At 1PN order, we see dependency on the new set of parameters and explicit GB coupling-dependent terms.
In section 4.4, we saw that the scalar multipoles consist only of compactly supported terms at 0.5PN order, and that the GB contribution to multipoles vanishes (see equation (39) and its subsequent paragraph).Thus, the deviation of the GB scalar waveform from that of ST gravity arises only from the novel GB coupling dependent terms of matter equation of motion (see equation ( 49)) at 1PN order.Also note that, in the GB coupling-dependent terms of the gravitational waveform, we see dependency on the parameters S + and S − , which are absent in the gravitational waveform of ST theory at 1PN order and start to appear at 1.5PN order.
Overall, most of the terms have the same form as in GR with modified coefficients.We expect more complicated structures to arise at 1.5PN order and beyond, with contributions from far-zone integrals (i.e.hereditary terms), surface EW moments, as well as dipole radiationreaction terms in the matter equation of motion.

Energy loss rate
In this section, we compute the energy dissipation due to tensor and scalar radiation from BH binary systems.As we will see in section 7, the energy flux rate is used to determine the phase evolution of GWs.

Tensor mode
The energy loss from tensor waves is given by We simplify this expression starting from the effects of TT projection operators.Based on the definition given in equation (18), it is straightforward to show that P ik P jl − 1 2 P i j P kl P im P jn − 1 2 P i j P mn = P km P ln − 1 2 Using this identity, the calculation simplifies to where we use the same notation as introduced in the introduction.Taking a time derivative and evaluating angular integrals using equation (107), we find that The explicit calculations of GB contributions to the above result can be found in appendix D. Following the PN convention of the waveforms, we call the lowest order piece of the tensor flux to be 0PN, which results from multiplying the 0PN piece of h i j TT by itself.At 0.5PN order, there is no tensor flux as the product of 0.5PN and 0PN terms has an odd number of unit vectors and thus vanishes upon angular integration.The remaining terms are at 1PN order, comprising 0PN-1PN, and 0.5PN-0.5PNproducts.

Scalar mode
For the scalar field, the energy loss is evaluated from Since we have defined the lowest-order tensor flux term as a 0PN term, the lowest-order piece of the scalar flux is a −1PN term, resulting from multiplying the −0.5PN piece of Φ by itself.We find the scalar energy flux to 0PN order to be In the above expression, there is no contribution to the flux at −0.5PN order because the product of the −0.5PN and 0PN pieces of Φ has an odd number of ni .At 0PN order, we have (−0.5PN)− (+0.5PN), and 0PN-0PN contributions.

Orbit equations and waveforms in the time domain
Having the energy flux and the conserved energy at hand, we determine the evolution of the orbital frequency and phase, which is needed to generate inspiral-waveform templates.Such templates are important for parameter inference studies to search for deviations from GR and quantify possible biases in source parameters that may mimic beyond-GR effects.
In this section, we first present the time-domain evolution of tensor and scalar waveforms for arbitrary GB coupling parameters but focusing on quasi-circular binary systems.In section 7.4, we assume the small coupling limit and compare our results for scalar waveforms against the NR results of reference [54].

Dynamics of quasi-circular inspirals
Here, we focus on the dynamics of orbits that are quasi-circular when they enter the sensitivity band of GW detectors.For such orbits, the only departure from circular motion is induced by radiation reaction, which, as we saw earlier in section 5, does not explicitly appear in the equations of motion to 1PN order.
Using the relative acceleration (see equation ( 49)) to solve for the circular condition r = ṙ = 0, we derive the angular velocity ω in terms of the relative separation to be It is useful at this stage to distinguish between the two commonly used PN parameters, namely which differ from their GR definition by an additional factor of ᾱ.At leading order, one has r 3 ω 2 = Gmᾱ + O(c −2 ).From equation ( 62), we find the relation between PN parameters to next-to-leading order, where we have substituted r = Gmᾱ/(c 2 γ PN ) + O(c −4 ) in the GB coupling dependent term.As mentioned earlier, despite the fact that the GB term here renders a term having a similar scaling with the orbital parameters as 3PN terms, it is a 1PN correction and indicates a different physical effect.Using equation (64) we can express the total binding energy (see equation ( 51)) of circular orbits in terms of the orbital frequency,

Gravitational wave phase evolution
In order to compute the full time-dependent waveforms, we need the evolution of the orbital phase angle ϕ, obtained from the angular velocity.In the adiabatic approximation, i.e. ω/ω 2 1, the luminosity of GWs is equal to the change in orbital energy averaged over a period, leading to the energy balance equation dE(x) dt = −F(x), (66) with F (x) being the total energy flux rate.Using φ = ω, the balance equation can be reformulated to where E (x) is the binding energy derivative with respect to the PN parameter x.In a PN approximation, this pair of differential equations can be solved in different ways, depending on how one chooses to expand the ratio F /E .Here, we choose the so-called Taylor T4 approximant, which is obtained from expanding the aforementioned ratio to the consistent PN order as a whole.For a review on different approximants, see reference [84].
From equation ( 65), we derive E (x) to relative 1PN to be The total energy flux, constituting equations ( 59) and ( 61), has the overall structure where the lower-index indicates the PN order of each term.Note that we calculated 0.5PN corrections to the leading order scalar waveform, which corresponds to a scalar energy flux that is complete at 0PN order, as the leading term is of −1PN order.Since we calculated the tensor flux to 1PN order, the 0.5 and 1PN scalar contributions to the total flux remain undetermined.The inclusion of these terms would increase the overall energy flux and thus the phase differences between sGB and GR.
To proceed with the expansion of F (x)/E (x), we shall distinguish between the regime where the scalar dipole flux dominates over the leading-order tensor flux, and vice versa.These regimes are commonly referred to as the scalar dipolar driven (DD) regime with F −1,S as the dominant term, and the tensor quadrupolar driven (QD) regime where, instead, F 0,T dominates.The DD regime is relevant when F −1,S (x) F 0,T (x), thus for frequencies as low as At much higher frequencies than this condition, the system is in the QD regime.Below, we find the 1PN expansion of F (x)/E (x) in the two regimes.

Dipolar driven regime.
For systems with large scalar dipole or very large separation, factoring out the dipolar scalar flux in equation ( 69) gives, For 1PN corrections to the phase at this regime, it is sufficient to keep the factor f DD 2 = (F 0,S + F 0,T )/F −1,S , being The last term in the above equation comes from the leading order term of the tensor flux.The rest of the terms are the 0PN scalar flux terms.Expanding the ratio F DD (x)/E (x) to 1PN order we find Note that this expression is complete at the relative 1PN order.

Quadrupolar driven regime.
For quadrupolar driven systems, the flux is expanded about the Newtonian-order term F 0,T such that The factor ξ = 5 ᾱ f DD 2 S 2 − /24, captures the 0PN flux terms.For the coefficient f nd 2 = F 1,T /F 0,T ≡ f nd 2,T + f nd 2,GB we have As mentioned above, these expansions miss a contribution from F 1,S which we are unable to compute within our approximations, having ĖS to 0PN order.Note that, as a result of these missing 1PN terms, the overall energy loss rate and hence the time evolution of waveforms in QD regime is only partially 1PN order, despite the waveform expressions being systematically computed to 1PN relative order in section 5.In the following, we set these contributions to zero, similar to what was done in a similar situation for tidal effects in GR [85].
Overall, the expansion of F QD (x)/E (x) to 1PN order becomes Figure 2. Orientation of the orbital plane and the sky plane based on the orthonormal triad {n, p, q}.p lie along the lines of nodes and defines an origin for the orbital phase angle ϕ. λ is the unit vector of the orbital angular momentum.
Having equations ( 73) and ( 76), together with an adequate choice of parameters, we numerically solve equation ( 67) to determine the time evolution of ϕ and ω in different regimes, and therefore, time-domain waveforms.

Gravitational waveforms in the time domain
To derive time-domain tensor and scalar waveforms, we parametrize the orbital motion by choosing the standard convention for the direction and orientation of the orbits, namely the orthonormal triad {n, p, q}, with n being the radial direction to the observer, p lying along the intersection of the orbital plane with the plane of the sky, and q = n × p.A schematic view of this choice can be found in figure 2.
The normal to the orbit is inclined at an angle i relative to n.The orbital phase angle of body A is measured from the line of nodes in a positive sense.We have that nAB = p cos ϕ + (q cos i + n sin i) sin ϕ , λ = −p sin ϕ + q cos i + n sin i cos ϕ, (77) where v = rω λ for circular orbits.GW detectors measure a linear combination of polarization waveforms h + (t) and h × (t) that are defined by the projections Applying equations ( 77) and ( 78) on equation ( 53) and simplifying combinations such as ni n j , λi λ j , and n(i λ j) , we find the 1PN polarization waveforms as functions of angular phase and orientation.The corresponding ready-to-implement waveforms can be found in appendix E.
We first display our results in figure 3 by plotting analytical GWs of BH binary systems in EdGB theory and comparing them with 1PN GR waveforms.Using equation (7), we write BH scalar charges, and the related parameters, in terms of α.The GB coupling parameter α is thus the only parameter that we have to specify.Exemplarily, we choose α = 0.01m 2 .Also note that, as shown in section (13b), the difference between the ssGB and EdGB phase evolution is relatively small compared to their overall deviations from the GR phasing.As a result, we only present the waveforms for the EdGB case.
By its very nature, the deviation from GR is most prominent in the high curvature regime, i.e. for small-mass BHs.Therefore, we consider quasi-circular BH binaries with a total mass of m = 15M and mass ratios q = 1/2 and q = 1/4.For q = 1 binaries, the deviations from GR waveforms are expected to be smaller due to the suppressed dipole radiation, as can also be seen in figure 2 of reference [74].
Figure 3 shows waveforms in the aforementioned two cases.The observer is viewing the orbit edge on, so that i = π/2, and thus only the h + (t) polarization is present.Due to the dissipation of energy in scalar field radiation, the inspiral of BH binaries in sGB is accelerated as compared to their GR counterpart and results in a GW phase shift.The scalar charge, and hence the scalar radiation, increases as the mass ratio decreases for fixed total mass of the binary.As can be seen in figure 3, the mass ratio has a high impact on the GW signal, which exhibits a larger phase shift for smaller mass ratios.Throughout the early inspiral evolution of the binaries, in both cases, the amplitudes of EdGB waves do not deviate significantly from the GR waveforms.This is because the GB coupling dependent terms are suppressed at large distances compared to other 1PN terms (see the discussion below equation ( 31)).Despite the amplitudes, the dephasing between EdGB and GR waveforms starts early in the evolution and could thus be a detectable feature.

Scalar waveforms in time domain
For scalar waveforms, similar steps to the ones described in the previous subsection should be taken in order to rewrite equation (55) in terms of the orbital phase and inclination angle.The final ready-to-implement 0.5PN order scalar waveform is reported in appendix E.
In figure 4 we display the scalar waveform emitted by quasi-circular binaries, corresponding to the gravitational waveforms shown in figure 3.In these examples, for the mass ratio q = 1/2, the scalar wave amplitude is suppressed with respect to the GW amplitude by an order of magnitude.As we decrease the mass ratio to q = 1/4, the increase of the scalar field amplitude during the late inspiral becomes stronger and, in fact, the amplitude becomes comparable to that of the gravitational radiation.In general, one expects the scalar radiation to increase with decreasing the mass ratio while keeping the total mass fixed.However, we should note that the PN approach is not reliable for extreme-mass-ratio systems.
We next compare PN scalar waveforms against those resulting from NR simulations of BH binary systems during the inspiral.For this purpose, we use spherical harmonics Y lm to decompose the scalar radiation into its radiative modes where, in order to express Φ in terms of spherical coordinate variables, we choose to work with p = −e Φ, q = e Θ, such that Θ = i and Φ = π/2 − ωt.It is easy to verify that the leading-order contribution to the radiation comes from Φ l=m modes, i.e.
We use equation ( 7) to re-write the scalar-charge dependent parameters in terms of the GB coupling.For ssGB gravity we have with the rest of the parameters being identically zero.Specifying equation ( 54) to quasi-circular orbits with the above-mentioned parameters, we find the various Φ l modes, required for equation (80), to be where the parameter x 0 = (Gmω/c 3 ) 2/3 + O(α) is the leading order term of the PN parameter x.We compare our PN scalar waveforms against previous PN calculations of reference [60] and also the results from a NR simulation reported in reference [54], with data kindly provided by the authors.This NR simulation implements an effective-field theoretical approach, which is valid to first order in the GB coupling parameter α in ssGB theory, and it covers about ten orbits before the binary BHs merge.Note that, for the purpose of this comparison, we have further incorporated a small coupling expansion on our general scalar waveform results at this stage, keeping only the terms up to O(α).As expected, our results agree with those of reference [54].For the leading PN terms of each mode.Also, our calculations expand these results to next-to-leading order, deriving a 0.5PN correction to the dipole amplitude Φ 11 at order O(α).In general, we also derive 0.5PN corrections to Φ 11 that are of higher order in the GB coupling constant.We neglect these corrections here as the comparison versus NR results only requires the O(α) corrections.
In order to have a meaningful comparison between the new PN dipolar radiation terms and the results of reference [54], we obtain the orbital frequency evolution from the derivative of the GW phase of numerical data, instead of using the analytical results of section 7.2.This means that the comparison between PN and NR results concerns only the amplitudes of the waveforms, while their phases agree by definition.
In figure 5 we show the NR scalar waveform for the l = m = 1 mode during the inspiral phase, and its comparison to the analytical results of this paper (red dashed curve) and that of reference [60] (blue dashed curve).The binary has a mass ratio of q = 1/4 and extraction radius R = 100 m.The waveforms are shifted in time such that t = 0 marks the merger time.We initially align the waves by minimizing the squared difference between the amplitude of the waveforms over an extended window in time, t − t merger : {−1410, −900}.In this way, we avoid splitting the waves by introducing arbitrary fudge factors.Overall, we see an approximate factor of 1.5 improvement in the matching of analytical results to numerical data due to the new 0.5PN contributions, as compared to the results of reference [60].Comparing the q = 1/2 waveforms shows a similar overall factor of improvement, and thus we do not show the plot here.The noise in the analytical amplitude is due to the precision of numerical data from which the phase evolution is extracted and does not affect the overall results.Based on these results, we claim that our analytical waveforms can thus be used as benchmarks for NR simulations.

Phase evolution in Fourier domain
For the purpose of data analysis, and since GW measurements are mainly sensitive to the phase evolution of waveforms, it is useful to provide waveforms in the Fourier domain.For instance, in theory-agnostic tests of gravity, the parametrized models used by the LVC in e.g.references [9,12], are also based on the analysis of GWs in the Fourier domain.
Here, we derive the 1PN Fourier phase analytically by following roughly the same steps as in section 7.2.Using the stationary phase approximation (SPA) and focusing on waveforms with GW Fourier phase ψ at frequency f = ω/π, the waveform in the frequency domain can be written as where t f is the time when the GW frequency becomes equal to the Fourier variable f , by solving dψ f (t)/dt = 0. a(t f ) represents the amplitude of the mode with frequency f .Using the so-called Taylor F2 approximation, we can solve for the time and orbital phase through where v = x 1/2 = (Gᾱmω/c 3 ) 1/3 and the subscript ref refers to the choice of reference point in the evolution.The gravitational phase is thus given by where v f ≡ (πGmᾱ f /c 3 ) 1/3 .As before, we split the calculation into DD regime and QD regime, and evaluate equation ( 84) term-by-term by re-expanding the expressions in the PN variable v and truncating them at 1PN order.

Dipolar-driven regime
In the DD regime, expanding the ratio E (v)/F DD (v) to 1PN order gives By substituting equation ( 86) into equation ( 84), and integrating term-by-term we find the phasing to be where one can easily verify that

Quadrupolar-driven regime
For the QD regime, by definition (see equation ( 70)), we expect the dipole radiation to be negligible.This means that we can use the parameter S − as an identifier of small terms, and thus, split the flux into two pieces as follows, where F dip indicates the tensor and scalar flux terms that depend on S − and F non-dip denotes the terms that do not depend on S − .
With this division, to first order in the small quantity F dip /F non-dip , the ratio E (v)/F (v) can be approximated as We obtain the dipolar and non-dipolar parts to be where ξ = (1 + S 2 + ᾱ/6) comes from the non-dipolar flux terms at Newtonian order.Note that this factor differs from the previously defined factor ξ by the fact that it only includes the flux terms that do not depend on S − explicitly.Similarly, the factor f nd 2 equals those terms in equation ( 75) that do not depend on S − , replacing ξ by ξ.Also, f d 2 can be found from the factor f DD 2 , and it equals the terms of f DD 2 that do not depend on 1/S 2 − .Overall, for equation (90) we find Integrating equation ( 84) by using the above expression we obtain Figure 6.Left: the GW phase difference between EdGB ψ EdGB and 1PN GR ψ GR shown for the systems with α = 0.01m 2 and m = 15M , as considered in figure 3, as well as for an equal-mass case.The blue dashed-dotted and red dashed lines correspond to q = 1/2 and q = 1/4, respectively.The black solid lines correspond to q = 1.Right: GW phase difference between EdGB and ssGB gravity for the aforementioned binaries.
with the coefficients being A useful tool for placing theory-agnostic constraints on deviations from GR is the parameterized post-Einsteinian (ppE) framework, which adds arbitrary parameters b ppE to the coefficients of powers of v in the expansion of the GW phasing.The mapping between ppE templates and waveforms in sGB theory is currently only known for the waveform at Newtonian order [7,60] (the first term in equation (95) corresponding to b ppE = 7).The higher-order terms computed here can be directly mapped to the ppE template with b ppE = −5, and thus potentially further improve the bounds on the coupling parameter e.g. by extending the work of [51,52].
In figure 6, we show an example of the Fourier-domain GW phase evolution in order to study the effect of dilatonic and shift-symmetric coupling functions on the phasing.We choose α = 0.01m 2 and m = 15M , as it corresponds to the systems considered in figure 3. The upper bound on frequency is chosen as f max = 2(6 3/2 πm) −1 ≈ 586 Hz and to simplify the comparison, all phases are aligned with the 1PN equal-mass phase in GR at the minimum frequency limit.
As the left panel shows, the difference between EdGB gravity and 1PN GR phase evolution is significant for the q = 1 binaries.The phase difference between EdGB and ssGB gravity is shown on the right panel, where we see that the overall difference is small compared to the difference with GR.For this specific example, the difference is O(1) cycles and thus is more difficult to distinguish from the deviations from GR.

Interpretation of GW constraints for fundamental parameters
There is an important point to note regarding the interpretation of results in quadratic gravity.If the action ( 1) is taken to be in the Jordan frame, i.e. if matter is minimally coupled to the metric g μν , our results directly provide the waveforms that would be measured by GW detectors and they can be employed to derive observational bounds on quadratic gravity.However, if equation ( 1) represents the low-energy effective action, of a string theory then it is, strictly speaking, in the Einstein frame.That is, the waveforms would have to be transformed to the physical (i.e.Jordan or string) frame to connect between GW detections and observational bounds on string theory.While there is no distinction between the two frames to leading order in the weak-field approximation, where the conformal factor transforming between them is A(φ) = 1 + O(φ), we urge the reader (and tester of strong-field gravity) to caution outside this approximation.

Conclusions
We have extensively studied the generation of GWs and scalar radiation in inspiralling blackhole binaries in sGB theories, which are characterized by the coupling of a scalar field to higher curvature terms, namely the GB invariant.Using the direct integration of Einstein equations in the PN expansion, we have computed the GWs and energy flux to relative 1PN order beyond the quadrupole emission, and the scalar analogue to relative 0.5PN order.We advanced the knowledge of GW and scalar waveforms in sGB theories beyond the leading results of reference [60] where the only effects came from the scalar field, and capture new effects due to curvature nonlinearities.
Furthermore, we compared our results in the small coupling regime against the scalar waveforms obtained with the lower-order PN scheme of reference [60] and obtained with NR simulations [54].
In order to compute the waveforms, we have skeletonized the compact bodies to effective point particles with scalar-field dependent masses to consistently include the effect of BH scalar-charges.We derived the two-body equations of motion to 1PN order based on the solutions of near-zone fields and the 1PN Lagrangian given in reference [61].In both the Lagrangian and the waveforms, we see that results differ from those of 1PN ST theory by (i) specific combinations of scalar-charge dependent parameters (see equation ( 42)), and (ii) additional GB coupling-dependent terms which have no analogue in ST theory and encode the impact of higher-curvature contributions.
This has consequences for a future effective-one-body (EOB) description of sGB theories, and suggests that the EOB results for ST theories can not be trivially extended to sGB theories.In addition, this difference would allow us to distinguish the two classes of theories through the analysis of mixed BH-neutron star binary inspirals which can exhibit only one scalarized body in each theory.
Focusing on compact binary systems in quasi-circular orbits, we have presented the tensor and scalar waveforms in a ready-to-implement form.That is, the GW polarizations and scalar wave in the time domain, equation (126) to (128), or the GW phasing in the Fourier domain, equations ( 87) and (93), can be employed directly to construct (phenomenological or EOB) waveform templates.Our results also provide a critical first step for constructing inspiral-merger-ringdown GW templates at high curvature regimes.
We have employed the SPA approximation to derive analytical expressions for the phasing in Fourier domain for systems whose inspiral is driven by the emission of scalar dipolar radiation, as well as those driven by tensor quadrupolar flux.In a companion paper [74], we have quantified the detectability of ssGB and EdGB phase deviations from GR, varying the binary parameters and the GB coupling.We have shown that the GB phase deviations are potentially detectable by A+LIGO/Virgo/KAGRA sensitivity bands [86][87][88], and thus the analytical results presented here can be used to put tight constraints on sGB theories, either through matching to ppE waveforms or through direct comparison against GW data.Future, groundand space-based GW detector networks combining high-frequency observations by, e.g. the Einstein Telescope [89] or Cosmic Explorer [90], with low-frequency observations by, e.g. the space-based LISA mission [91], open new avenues for multi-band tests of gravity.Our results, in principle, are extendable to cover the phase evolution across the frequency bands and enable novel, multi-band tests of gravity.A detailed analysis is left for future work.
We emphasize that our results and the methods applied here are not restricted to specific choices of the coupling function, nor to the weak coupling limit, and thus can potentially be extended to explore dynamical scalarization or de-scalarization effects of BH binaries [41] during the early inspiral phase.In particular, in the current treatment we have neglected the finite-size effects which may account for both the dynamical (de-)scalarization as well as scalar-induced tidal interaction.We leave the investigation of such effects for future work.
we have and in the case of b = c, Equation ( 98) can be re-written to relate the partial derivative of the gothic metric to that of its inverse, A.2. Christoffel symbols and curvature tensors By using equations ( 97)-(99), we can express quantities in terms of gothic variables.In particular, the Christoffel symbols become To derive the curvature tensors, we compute the partial derivative of the Christoffel symbols.In terms of these, the Ricci scalar in gothic formulation is The rest of the expressions for curvature tensors can be found in the corresponding Mathematica notebook.
The two important curvature combinations that we need for the field equations are the dual Riemann tensor and the GB scalar.When written in terms of the gothic metric, these quantities become very long.As we work in the weak-field limit of these quantities, we do not report their full expressions here.We consider that it is sufficient to describe their overall structure.In particular, for the dual Riemann tensor, the metric and its derivatives appear in the following combinations: * R abcd * ∝ 1 −g 3 (g In the main text, we called this quantity * Rabcd * .The expression for the R GB is derived from the above mentioned curvature quantities through equation (2).
from permuting the indices on nL−2P δ P .As an example, this relation gives The STF product of unit vectors are such that their integration over solid angle is zero.
Converting back to non-STF form would lead to the following identities The simplest integral that appears in the two-index EW moment calculation is (∇U) 2 d 3 x.Applying integration-by-parts on this term gives The surface integral is evaluated on the boundary of the near zone ∂M as a sphere with radius R, such that we can write x i = Rn i and d 2 S k = R 2 nk .To evaluate this term we shall expand U and U ,k in inverse powers of R and ignore the integrands that depend on R. The R-independent piece of this specific surface integral contains the O(R −4 ) expansion of UU ,k which gives an integrand with odd number of unit vectors and thus vanishes (see equation ( 105)).For the first volume integral, substituting the Laplacian of U easily shows that − M UU ,kk x i j d 3 x = 4π A,B =A m A m B r x i j A .The contribution from the second volume integrals is zero.This can be shown by integrating-by-parts again, −2 M UU ,(i x j) d 3 x = − ∂M U 2 x (i d 2 S j) + M U 2 δ i j d 3 x, which shows that the surface integral vanishes due to an odd number of unit vectors.The volume integral is also ignored as it vanishes through TT projection.
Other important integrals that we encounter in the calculations (such as U ,m U ,lm d 3 x) is the UU ,lm d 3 x term, which can be written as [63]: ,where in the second line we have changed the integration variables from x to y = x − x A and have changed ŷlm to its STF form.Also we note that r = x A − x B and nAB = r/r.There are two cases to consider: A = B and A = B.For both, the infinite series of surface integrals vanishes as the terms either depend on R or average to zero because of an odd number of unit vectors.
For the A = B case, r = 0, so the volume integral can be evaluated easily.When A = B, we may use of the following expansions: > y q dy = 2l + 1 (l + q + 1)(l − q) r q AB .(112) Overall we find

Appendix C. Far-zone contribution to waveforms
In this appendix, we focus on the solution of equation ( 13) with far-zone field point and show that the contributions are beyond 1PN order.In the far-zone region, the source terms μ s and μ i j are composed purely of field terms as, by definition, there is no matter source present in this region.It can be shown that (see e.g.[63]) the fields at any intermediate distance R are given by: with a new set of multipole moments defined as: Note that equation (114) reduces to equations ( 14) and (36) in the limit where R R. As can be seen from the structure of terms in equations ( 20) and (38), the sources in the far-zone are composed only of the field components N and Φ. Finding these components by Applying the same procedure on equation ( 55) gives the 0.5PN scalar waveform: x −1/2 Φ−1/2 + Φ0 + x

Figure 1 .
Figure 1.The past null cone of the field point x and its intersection with the near-zone, for a field point located in the far-zone.

Figure 3 .
Figure 3.Time evolution of gravitational waveform and orbital frequency for an m = 15M binary with q = 1/2 (left) and q = 1/4 (right), and α = 0.01m 2 .Blue dashed curves correspond to EdGB waveforms and black curves correspond to 1PN GR waveforms.Orbits are viewed edge-on (i = π/2), and t = 0 corresponds to f ≈ 50 Hz.The green shaded regions (top) show 1 s interval of the waveforms at intermediate times, while red shaded regions (middle) show the 0.3 s intervals where EdGB binaries are close to the merger.

Figure 4 .
Figure 4.The time evolution of scalar field waveform for m = 15M binary with (left) q = 1/2 and (right) q = 1/4, and α = 0.01m 2 , corresponding to figure 3 binaries and time intervals.Orbits are viewed edge-on (i = π/2) and thus the waves deviate from a full sinusoidal.

Figure 5 .
Figure 5.Time evolution of the scalar mode Φ 11 (left) and its absolute value (right), re-scaled by the extraction radius R = 100 m.Black curves indicate results from NR simulations.The blue dashed curve indicates the analytical PN inspiral results of reference [60] to −0.5PN order.The red dashed curve shows the 0.5PN order results of this paper.