* 1 GRAPPA, Anton Pannekoek Institute for Astronomy and Institute of High-Energy Physics, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands

  * 2 Nikhef, Science Park 105, 1098 XG Amsterdam, The Netherlands

  * 3 Institute for Theoretical Physics, Utrecht University, Princetonplein 5, 3584 CC Utrecht, The Netherlands

  * 4 Instituto de Ciencias Nucleares (ICN), Universidad Nacional Autónoma de México (UNAM), Circuito Exterior C.U., A.P. 70-543, México D.F. 04510, México.

  * 5 Illinois Center for Advanced Studies of the Universe & Department of Physics, University of Illinois at Urbana-Champaign, Urbana, Illinois 61801, USA

#  Post-Newtonian Gravitational and Scalar Waves in Scalar-Gauss-Bonnet
Gravity

Banafsheh Shiralilou1,2, Tanja Hinderer1,3, Samaya M. Nissanke1,2, Néstor
Ortiz4, Helvi Witek5

###### Abstract

Gravitational waves emitted by black hole binary inspiral and mergers enable
unprecedented strong-field tests of gravity, requiring accurate theoretical
modelling 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 first post-Newtonian
order beyond the quadrupole approximation and calculate new contributions from
nonlinear curvature terms. We quantify the effect of these terms and provide
ready-to-implement gravitational wave and scalar waveforms as well as the
Fourier domain phase for quasi-circular binaries. We also perform a parameter
space study, which indicates that the values of black hole scalar charges play
a crucial role in the detectability of deviation from General Relativity. We
also compare the scalar waveforms to numerical relativity simulations to
assess the impact of the relativistic corrections to the scalar radiation. Our
results provide important foundations for future precision tests of gravity.

###### Contents

  1. 1 Introduction
  2. 2 Einstein-Scalar-Gauss-Bonnet gravity
    1. 2.1 Gravitational action and the field equations
    2. 2.2 Different viewpoints of the action and their consequences
    3. 2.3 Matter action
  3. 3 Relaxed-field equations and wave generation
  4. 4 Calculating the tensor and scalar waveforms
    1. 4.1 General structure of near-zone calculation
    2. 4.2 Formal structure of near-zone source
    3. 4.3 Evaluation of two-body Epstein-Wagoner moments
      1. 4.3.1 Two-index moment $I^{ij}_{EW}$
      2. 4.3.2 Three-index moment $I^{ijk}_{EW}$
      3. 4.3.3 Four-index moment $I^{ijkl}_{EW}$
    4. 4.4 Scalar multipole moments and waveform
  5. 5 Two-body waveforms in the center of mass frame
    1. 5.0.1 Center of mass and relative acceleration to 1 Post-Newtonian order
    2. 5.0.2 Final expressions for the waveform
  6. 6 Energy loss rate
    1. 6.1 Tensor mode
    2. 6.2 Scalar mode
  7. 7 Orbit equations and waveforms in the time domain
    1. 7.1 Dynamics of quasi-circular inspirals
    2. 7.2 Gravitational wave phase evolution
      1. 7.2.1 Dipolar driven regime
      2. 7.2.2 Quadrupolar driven regime
    3. 7.3 Gravitational waveforms in the time domain
    4. 7.4 Scalar waveforms in time domain
  8. 8 Phase evolution in Fourier domain
    1. 8.1 Dipolar-driven regime
    2. 8.2 Quadrupolar-driven regime
    3. 8.3 Interpretation of GW constraints for fundamental parameters
  9. 9 Conclusions
  10. 10 Acknowledgment
  11. A Curvature Tensors in terms of gothic metric
    1. A.1 Transformation rules
    2. A.2 Christoffel symbols and curvature tensors
    3. A.3 Weak-field limit
  12. B Field integrals and calculation techniques
  13. C Far-zone contribution to waveforms
  14. D Calculation of coupling-dependant energy flux terms
  15. E  Final polarization waveforms and scalar waveform

##  1 Introduction

The breakthrough discovery of gravitational waves (GWs) emitted by merging
black holes (BHs) Abbott  _et al._ (2016a) has opened a unique new avenue for
probing unexplored regimes of the universe. GWs have already provided
unprecedented insights into the origin and population of BHs, including
observations closing the mass-gaps Abbott  _et al._ (2019a, 2020a, 2020b).
They also hold the key to test gravity in its strong-field, nonlinear regime
Yunes and Siemens (2013); Berti  _et al._ (2015); Yunes  _et al._ (2016); Yagi
and Stein (2016); Abbott  _et al._ (2019b, 2020c); Carson and Yagi (2020);
Abbott  _et al._ (2016b).

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 Yunes  _et al._ (2016); Abbott  _et
al._ (2019b, 2020c, 2016b); Isi  _et al._ (2019). These tests are performed
only on single coefficients associated to a scaling with frequency Yunes and
Pretorius (2009); Will (2014); Abbott  _et al._ (2019b, 2020c) and thus, the
interpretation of theoretical constraints remains limited (Yunes  _et al._ ,
2016). 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
$\mathcal{R}_{\mathrm{G}B}=R_{\mu\nu\sigma\rho}R^{\mu\nu\sigma\rho}-4R_{\mu\nu}R^{\mu\nu}+R^{2}$,
with a coupling parameter $\alpha$. 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 Boulware and Deser
(1985); Gross and Sloan (1987); Kanti  _et al._ (1996); (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 Charmousis  _et al._ (2012a, b); (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 (Stelle, 1977); 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 Kovács (2019); Kovács and Reall (2020a, b)).

The coupling, $f(\phi)$, between the dynamical scalar field and the GB
invariant determines the “flavour” of sGB gravity Antoniou  _et al._ (2018a,
b). For the sake of the summary we label them as “type I” and “type II.” In
our notation, type I sGB corresponds to the subclass of theories for which
$f^{\prime}(\phi)$ never vanishes. Representative examples include shift-
symmetric ($f=2\phi$) or dilatonic ($f\sim e^{2\phi}$) couplings. In this case
the space-time curvature always sources the scalar field and, thus, inevitably
yields hairy black holes Mignemi and Stewart (1993); Kanti  _et al._ (1996);
Sotiriou and Zhou (2014); Benkel  _et al._ (2017, 2016); Ripley and Pretorius
(2019); Yunes and Stein (2011); Pani and Cardoso (2009); Barausse and Yagi
(2015); Kleihaus  _et al._ (2011, 2016) that can form dynamically Benkel  _et
al._ (2016, 2017); Ripley and Pretorius (2020). Type II sGB corresponds to the
subclass of theories for which $f^{\prime}(\phi)$ can vanish for some values
of the scalar field. Examples include a quadratic or Gaussian coupling
function. Then, the Schwarzschild and Kerr solutions of GR still exist, and
they can spontaneously scalarize Silva  _et al._ (2018); Doneva and Yazadjiev
(2018); Berti  _et al._ (2020); Silva  _et al._ (2020); Dima  _et al._ (2020);
Herdeiro  _et al._ (2021); Berti  _et al._ (2020); Collodel  _et al._ (2020);
Doneva and Yazadjiev (2021); Doneva  _et al._ (2020); Silva  _et al._ (2020);
East and Ripley (2021a); Annulli (2021). Likewise, neutron stars can
spontaneously scalarize in Type II sGB gravity Silva  _et al._ (2018). Such
scalarized black holes or neutron stars can form dynamically Ripley and
Pretorius (2020); Kuan  _et al._ (2021).

Type I sGB gravity has been tested extensively against observations of low-
mass x-ray binaries Yagi (2012) and a Bayesian parameter estimation for GW
detections Nair  _et al._ (2019). The most recent observational bounds of
$\sqrt{\alpha}\lesssim 1.7$ km were obtained by analysing signals of
LIGO/VIRGO’s first two GW catalogs Perkins  _et al._ (2021); Wang  _et al._
(2021). Type II sGB gravity has, so far, remained unconstrained from
observations.

To-date, modelling 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 Witek  _et al._ (2019); Okounkova (2020)
or the decoupling limit Silva  _et al._ (2020). The numerical study of fully
non-linear field equations for general couplings has recently been advanced in
the direction of initial data construction Kovacs (2021), and numerical
simulations East and Ripley (2021b) through a modified generalized harmonic
formulation of the evolution equations Kovács and Reall (2020a). This has also
enabled the recent study of spontaneous BH scalarization East and Ripley
(2021a). The ringdown spectrum has been studied extensively in Ref. Blázquez-
Salcedo  _et al._ (2016); Pierini and Gualtieri (2021). However, the early
inspiral waveform, typically first constructed with Post-Newtonian (PN)
methods, has only been modelled at leading (Newtonian) order Yagi  _et al._
(2012). 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 Julié and
Berti (2019).

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 (Epstein and Wagoner, 1975) and extended by Will, Wiseman
and Pati Will and Wiseman (1996); Pati and Will (2000, 2002). The PN
assumption of a gravitationally bound source in the weak-field, low-velocity
regime $Gm/rc^{2}\approx v^{2}/c^{2}\ll 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 Landau and Lifshitz (1975), 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 Fig. 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 to
provide non-divergent results in agreement with other approaches such as the
multipolar post-Minkowskian formalism (Blanchet, 2014). We note that this
method has also been applied to massless scalar-tensor (ST) gravity
theories—which admits scalarized neuron star solutions Damour and Esposito-
Farese (1993)—providing the equation of motion for compact binaries to 3PN
order (Bernard, 2018; Mirshekari and Will, 2013), GWs to 2PN (Lang, 2014;
Sennett  _et al._ , 2016), and scalar waves to 1.5PN (Lang, 2015) 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 Shiralilou  _et al._ (2020) 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 Kanti  _et al._ (1996); Sotiriou and Zhou (2014), 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 Ref.
Julié and Berti (2019). 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 Tahura and Yagi (2018) 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 a 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 Carson  _et al._ (2020). 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 Sec. 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 Sec. 3,
and outline the procedure for finding the solution at large distances from the
source. In Sec. 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 Sec. 5 we derive the equations of motion using
the two-body Lagrangian from Ref. Julié and Berti (2020), and transform the
waveforms to the center-of-mass (CM) frame. In Sec. 6, we solve for the
gravitational and scalar energy loss rate, from which we find the time domain
and frequency domain waveforms in Sec. 7 and Sec. 8, respectively. In Sec. 7,
we also compare our results for the scalar waveform to the numerical
relativity waveforms of Ref. Witek  _et al._ (2019) computed in the small
coupling approximation. Finally, in Sec. 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 multi-index notation for products of vector
components: $x^{ijk}\equiv x^{i}x^{j}x^{k}$, and a capital letter superscript
denotes a product of that dimensionality: $x^{L}\equiv
x^{k_{1}}x^{k_{2}}\cdots x^{k_{l}}$. A subscript preceded by a comma such as
$A,_{\mu}$ stands for the partial derivative $\partial_{\mu}A$, while
$\nabla_{\mu}A$ indicates the covariant derivative associated to the metric.
We label the two bodies in a binary system by $A$, $B$.

##  2 Einstein-Scalar-Gauss-Bonnet gravity

###  2.1 Gravitational action and the field equations

The gravitational action of sGB theory, in the presence of matter, is written
in the Einstein frame as

$S=\frac{c^{3}}{16\pi
G}\int_{M}d^{4}x\sqrt{-g}\left[R-2(\nabla\phi)^{2}+\alpha
f(\phi)\mathcal{R}^{2}_{GB}\right]+S_{m}[\Psi_{m},\mathcal{A}^{2}(\phi)g_{\mu\nu}]\,,$
(1)

where $R$ is the Ricci scalar on the 4-dimensional manifold $M$ with the
metric $g_{\mu\nu}$ in the Einstein frame. $S_{m}$ is the matter action with
$\Psi_{m}$ denoting the matter fields with a generic non-minimal coupling to
the metric. The fundamental coupling constant $\alpha$ has dimensions of
length squared, and $f(\phi)$ is the dimensionless coupling function.

In this paper, the choice $f(\phi)=e^{2\phi}/4$ corresponds to Einstein
dilaton Gauss Bonnet (EdGB) gravity Kanti  _et al._ (1996), and
$f(\phi)=2\phi$ to shift symmetric sGB (ssGB) gravity Sotiriou and Zhou
(2014). It is useful to rewrite the GB scalar as

$\mathcal{R}^{2}_{GB}=^{*}R_{abcd}^{*}R^{abdc}\,,$ (2)

with ${}^{*}R_{c\mu\nu
d}^{*}=\frac{1}{4}\epsilon^{abef}R_{efgh}\epsilon^{ghcd}$ the dual of Riemann
tensor and $\epsilon_{\alpha\beta\mu\nu}$ the anti-symmetric Levi Civita
tensor.

The field equations resulting from action in Eq. (1) are

$\displaystyle\Box\phi=-\frac{\alpha}{4}f^{\prime}(\phi)\mathcal{R}^{2}_{GB}+\frac{4\pi
G}{c^{4}\sqrt{-g}}\frac{\partial S_{m}}{\partial\phi}\,,$ (3a) $\displaystyle
G^{\mu\nu}=\frac{8\pi
G}{c^{4}}T^{\mu\nu}+2\nabla^{\mu}\phi\nabla^{\nu}\phi-g^{\mu\nu}(\nabla\phi)^{2}-4\alpha\left({}^{*}R^{*c\mu\nu
d}\nabla_{cd}f(\phi)\right)\,,$ (3b)

with $\square\equiv g^{\alpha\beta}\nabla_{\alpha}\nabla_{\beta}$ being the
d’Alembertian operator, $G^{\mu\nu}$ the Einstein tensor and
$T^{\mu\nu}=(-2/\sqrt{-g})(\partial S_{m}/\partial g_{\mu\nu})$ the
distributional matter energy-momentum tensor.

###  2.2 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. Refs. Yunes
and Stein (2011); Yagi  _et al._ (2012)) and corresponds to setting the
conformal factor $\mathcal{A}(\phi)=1$ in Eq. (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, Eq. (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(\phi)=e^{-2\phi}$ for which the conformal factor is
$\mathcal{A}(\phi)=e^{\phi}$ (see e.g. Ref. Metsaev and Tseytlin (1987)). Both
of these interpretations of the action become reconciled in the weak-field
regime, where $\mathcal{A}(\phi)=e^{\phi}=1+\mathcal{O}(\phi)$. 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.

###  2.3 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 Eardley (1975a) for ST theories (Eardley, 1975b; Damour
and Esposito-Farese, 1992a), and later generalized to Einstein-Maxwell-dilaton
theories Julié (2018), 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

$S_{m}[g_{\mu\nu},\phi,x_{A}^{\mu}]=-c\int
M_{A}(\phi)\sqrt{-g_{\mu\nu}dx^{\mu}_{A}dx^{\nu}_{A}}+(A\leftrightarrow B)\,.$
(4)

Here $dx_{A}^{\mu}$ 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}(\phi)$.

Note that the point particle description of Eq. (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}(\phi)$ can be
parametrized by its expansion about the background scalar field $\phi_{0}$:

$M_{A}=m_{A}\Big{[}1+\alpha_{A}^{0}\delta\phi+\frac{1}{2}\left((\alpha^{0}_{A})^{2}+\beta_{A}^{0}\right)\delta\phi^{2}\Big{]}+\mathcal{O}(\delta\phi^{3})\,,$
(5)

with $\delta\phi=\phi-\phi_{0}$, and the strong-field parameters

$\left.\alpha_{A}^{0}=\frac{d\,ln\,M_{A}(\phi)}{d\phi}\right|_{\phi=\phi_{0}},\hskip
25.0pt\hskip
25.0pt\left.\beta_{A}^{0}=\frac{d\,\alpha_{A}(\phi)}{d\phi}\right|_{\phi=\phi_{0}}.$
(6)

Also, $m_{A}$ is the asymptotic value of the particle’s mass in the Einstein
frame. The parameter $\alpha_{A}^{0}$ 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 $\alpha_{A}^{0}$
corresponding to static, spherically symmetric BHs has been derived to fourth
order in the coupling parameter in Ref. (Julié and Berti, 2019). For example,
to first order, the scalar charge is

$\alpha_{A}^{0}\equiv-\alpha f^{\prime}(\phi_{0})/2m_{A}^{2}\,.$ (7)

##  3 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 $\mathfrak{g}^{ab}=\sqrt{-g}g^{ab}$ called gothic
metric. By defining

$H^{\alpha\mu\beta\nu}\equiv\mathfrak{g}^{\mu\nu}\mathfrak{g}^{\alpha\beta}-\mathfrak{g}^{\alpha\nu}\mathfrak{g}^{\beta\mu}\,,$
(8)

it can be shown that the following is an identity, valid for any space-time:

$\partial_{\mu\nu}H^{\mu\alpha\nu\beta}=(-g)\left(2G^{\alpha\beta}+\frac{16\pi
G}{c^{4}}\,t^{\alpha\beta}_{LL}\right)\,,$ (9)

where $t^{\alpha\beta}_{LL}$ is the Landau-Lifshitz energy-momentum pseudo-
tensor as given in Ref. Landau and Lifshitz (1975).

To incorporate sGB gravity into this framework, we rewrite the expression for
$G^{\mu\nu}$ given by Eq. (3b) in terms of the gothic metric. In order to
analyze the behaviour 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
$\phi_{0}$. The perturbation variables are thus defined through

$\mathfrak{g}^{\mu\nu}=\eta^{\mu\nu}+h^{\mu\nu}\,,\quad\quad\mathrm{and}\quad\quad\phi=\phi_{0}+\delta\phi\,.$
(10)

By choosing the Harmonic gauge defined as
$\partial_{\nu}\mathfrak{g}^{\mu\nu}=0$, Eq. (9) can be written as a wave
equation for the perturbations, with sources

$\begin{split}\Box h^{\alpha\beta}&=\frac{16\pi
G}{c^{4}}\mu^{\alpha\beta}\,,\\\
\mu^{\alpha\beta}&=(-g)T^{\alpha\beta}_{m}+\frac{c^{4}}{16\pi
G}(\Lambda_{GB}^{\alpha\beta}+\Lambda^{\alpha\beta}_{GR})\,,\\\
\Lambda^{\alpha\beta}_{GR}&=\frac{16\pi
G}{c^{4}}(-g)t_{LL}^{\alpha\beta}+h^{\alpha\nu}_{,\mu}h^{\beta\mu}_{,\nu}-h^{\mu\nu}h^{\alpha\beta}_{\,\,,\mu\nu}\,,\\\
\Lambda_{GB}^{\alpha\beta}&=-8\alpha(-g)\left({}^{*}\hat{R}^{*c\alpha\beta
d}f(\phi)_{,cd}\right)+4\phi_{,c}\phi_{,d}\big{(}\mathfrak{g}^{\alpha
c}\mathfrak{g}^{\beta
d}-\frac{1}{2}\mathfrak{g}^{\alpha\beta}\mathfrak{g}^{cd}\big{)}\,,\end{split}$
(11)

where $\Lambda_{GB}^{\alpha\beta}$ captures the scalar field and the nonlinear
GB contributions, with ${}^{*}\hat{R}^{*c\alpha\beta d}$ the gauge fixed, dual
Riemann tensor written in terms of gothic variables. The explicit expression
for $\Lambda_{GB}$ in terms of the metric and scalar perturbations is given in
Appendix (A.3).
The scalar field equation, Eq. (3a), is itself a wave equation. In terms of
the gothic metric we have

$\begin{split}\Box\delta\phi&=\frac{4\pi G}{c^{4}}\mu_{s}\,,\\\
\mu_{s}&=\frac{S_{m,\phi}}{\sqrt{-g}}-\frac{c^{4}}{16\pi G}\alpha
f^{\prime}(\phi)\hat{\mathcal{R}}^{2}_{GB}\,,\end{split}$ (12)

with $\hat{\mathcal{R}}^{2}_{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

$\displaystyle h^{\mu\nu}(x)$
$\displaystyle=\frac{4G}{c^{4}}\int\frac{\mu^{\mu\nu}\left(t^{\prime},\mathbf{x}^{\prime}\right)\delta\left(t^{\prime}-t+\left|\mathbf{x}-\mathbf{x}^{\prime}\right|/c\right)}{\left|\mathbf{x}-\mathbf{x}^{\prime}\right|}d^{4}x^{\prime}\,,$
(13a) $\displaystyle\delta\phi(x)$
$\displaystyle=\frac{G}{c^{4}}\int\frac{\mu_{s}\left(t^{\prime},\mathbf{x}^{\prime}\right)\delta\left(t^{\prime}-t+\left|\mathbf{x}-\mathbf{x}^{\prime}\right|/c\right)}{\left|\mathbf{x}-\mathbf{x}^{\prime}\right|}d^{4}x^{\prime}\,.$
(13b)

The integrals here are over the past null cone of the field point
$x=(ct,\mathbf{x})$.
In principle, the explicit solution of Eq. (13) depends on the position of the
field point $\mathbf{x}$ relative to the source $\mathbf{x^{\prime}}$, i.e.
the size of the quantity $\|\mathbf{x}-\mathbf{x}^{\prime}\|=R$ relative to
other scales in the system. Defining the characteristic size of the source to
be $\mathcal{S}$, the near-zone is defined as the region with $R<\mathcal{R}$,
where $\mathcal{R}\sim\mathcal{S}/v$ is the characteristic wavelength of GWs
from the system. The region outside the near-zone is called far-zone (see Fig.
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.

Given the position of the field point x, one can evaluate each integral in two
separate pieces: an integral over the source point $\mathbf{x}^{\prime}$ 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 Refs. (Will and
Wiseman, 1996; Pati and Will, 2000).

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 $\mu^{\mu\nu}$, and its scalar
analogue $\mu_{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 Will and Wiseman (1996).

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 Sec. 5.

##  4 Calculating the tensor and scalar waveforms

###  4.1 General structure of near-zone calculation

Having the field point in the far-away zone, with $R\gg\mathcal{R}$, it can be
shown Will and Wiseman (1996) that the near-zone contribution to $h^{ij}$ ,
named $h_{\mathcal{N}}^{ij}$, is given as

$h^{ij}_{\mathcal{N}}(x)=\frac{2G}{Rc^{4}}\frac{d^{2}}{dt^{2}}\sum_{l=0}^{\infty}\hat{n}_{k_{1}...k_{l}}I_{EW}^{ijk_{1}...k_{m}}(\tau)+\mathcal{O}(R^{-2})\,,\hskip
25.0pt$ (14)

where $\hat{n}_{k_{i}}$ is the unit normal vector pointing from the source to
the detector. The tensor multipole moments $I_{\mathrm{EW}}^{ijk_{1}...k_{m}}$
are known as Epstein-Wagoner (EW) moments and are functions of retarded time
$\tau=t-R/c$. They are given by

$\begin{split}&I_{\mathrm{EW}}^{ij}\equiv\frac{1}{c^{2}}\left[\int_{\mathcal{M}}\mu^{00}x^{ij}d^{3}x+I_{\mathrm{EW}(\mathrm{surf})}^{ij}\right]\,,\\\
&I_{\mathrm{EW}}^{ijk}\equiv\frac{1}{c^{3}}\left[\int_{\mathcal{M}}\left(2\mu^{0(i}x^{j)k}-\mu^{0k}x^{ij}\right)d^{3}x+I_{\mathrm{EW}(\mathrm{surf})}^{ijk}\right]\,,\\\
&I_{\mathrm{EW}}^{ijk_{1}\cdots
k_{m}}\equiv\frac{2}{m!c^{2}}\frac{d^{m-2}}{d(ct)^{m-2}}\int_{\mathcal{M}}\mu^{ij}x^{k_{1}\cdots
k_{m}}d^{3}x\quad(m\geq 2)\,.\end{split}$ (15)

with the so-called surface moments being

$\begin{split}&\frac{d^{2}}{dt^{2}}I_{\mathrm{EW}(\mathrm{surf})}^{ij}=\oint_{\partial\mathcal{M}}\left(4\mu^{l(i}x^{j)}-\left(\mu^{kl}x^{ij}\right)_{,k}\right)\mathcal{R}^{2}\hat{n}^{l}d^{2}\Omega\,,\\\
&\frac{d}{dt}I_{\mathrm{EW}(\mathrm{surf})}^{ijk}=\oint_{\partial\mathcal{M}}\left(2\mu^{l(i}x^{j)k}-\mu^{kl}x^{ij}\right)\mathcal{R}^{2}\hat{n}^{l}d^{2}\Omega\,.\,\,\end{split}$
(16)

These two- and three- index surface moments are the results of using the
conservation law $\partial_{\alpha}\mu^{\alpha\beta}=0$ in the derivation of
Eq. (14) (see Ref. Will and Wiseman (1996)). As before, $\mathcal{R}$ is the
radius of the boundary region $\partial\mathcal{M}$ of the near-zone, and we
have used the fact that the surface element is
$dS^{k}=\mathcal{R}^{2}d\Omega^{2}\hat{n}^{k}$ 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

$h^{ij}_{\mathcal{N},_{TT}}(x)=\frac{2G}{Rc^{4}}\frac{d^{2}}{dt^{2}}\left\\{I^{ij}+\hat{n}_{k}I^{ijk}+\hat{n}_{kl}I^{ijkl}\right\\}_{{}_{TT}}.$
(17)

In the above equation, the TT projection operator acting on a tensor $A^{kl}$
is such that

$h_{\mathrm{TT}}^{ij}=\left(P^{ik}P^{jl}-\frac{1}{2}P^{ij}P^{kl}\right)A^{kl}\,,\hskip
25.0ptP^{ij}=\delta^{ij}-\hat{n}^{i}\hat{n}^{j}\,,$ (18)

with the properties $P^{ii}=2,P^{ij}P^{ij}=2,$ and $P^{ij}P^{ik}=P^{jk}$.

###  4.2 Formal structure of near-zone source

To compute EW moments requires expanding the components of the source
$\mu^{\alpha\beta}$ to 1PN order. The 1PN expansions will depend on the
solution of Eq. (13) found in the near-zone. To 1PN order, these potentials
are computed and can be found in Ref. Julié and Berti (2019).

For the expansion of the source terms, we begin by defining a simplifying
notation for the different components of the fields, following Ref. (Pati and
Will, 2000),

$\begin{split}&N\equiv h^{00}\sim\mathcal{O}(c^{-2})\,,\quad\,\,\,K^{i}\equiv
h^{0i}\sim\mathcal{O}(c^{-3})\,,\\\ &B^{ij}\equiv
h^{ij}\sim\mathcal{O}(c^{-4})\,,\quad B\equiv
h^{ii}\sim\mathcal{O}(c^{-4})\,,\\\
&\Phi\equiv\delta\phi\sim\mathcal{O}(c^{-2})\,.\end{split}$ (19)

Using the same strategies used in GR calculations, we find that the weak-field
expansion of the source $\mu^{\alpha\beta}$ to 1PN order (see Appendix. (A.3)
for the expansion of the GB terms), is

$\begin{split}\mu^{00}&=m_{A}c^{2}\left[1+\left(\frac{v_{A}^{2}}{2c^{2}}+\frac{3}{4}N+\alpha^{0}_{A}\Phi\right)\right]\delta^{3}(\mathbf{x}-\mathbf{x}_{A}(t))+(A\leftrightarrow
B)\\\ &-\frac{c^{4}}{16\pi G}\left[\frac{7}{8}(\nabla
N)^{2}-2(\nabla\Phi)^{2}\right]-2\frac{\alpha f^{\prime}(\phi_{0})c^{4}}{16\pi
G}\left(\partial_{ij}\Phi\partial_{ij}N\right)+\mathcal{O}(c^{-2})\,,\hskip
25.0pt\hskip 25.0pt\hskip 25.0pt\end{split}$ (20a)
$\mu^{0j}=m_{A}c\,v_{A}^{j}\delta^{3}(\mathbf{x}-\mathbf{x}_{A}(t))+(A\leftrightarrow
B)+\mathcal{O}(c^{-1})\,,\hskip 25.0pt\hskip 25.0pt\hskip 25.0pt\hskip
25.0pt\hskip 25.0pt\hskip 25.0pt\hskip 25.0pt$ (20b)
$\begin{split}\mu^{ij}&=m_{A}v_{A}^{i}v_{A}^{j}\delta^{3}(\mathbf{x}-\mathbf{x}_{A}(t))+(A\leftrightarrow
B)\\\ &+\frac{c^{4}}{16\pi
G}\Bigg{[}\frac{1}{4}\bigg{(}\partial^{i}N\partial^{j}N-\frac{\delta^{ij}(\nabla
N)^{2}}{2}\bigg{)}+4\Big{(}\partial^{i}\Phi\partial^{j}\Phi-\frac{\delta^{ij}}{2}(\nabla\Phi)^{2}\Big{)}\Bigg{]}\\\
&-2\frac{\alpha f^{\prime}(\phi_{0})c^{4}}{16\pi
G}\Big{[}\Big{(}\partial_{ki}\Phi\partial_{kj}N-\partial_{kj}\Phi\partial_{ki}N\Big{)}+\Big{(}\partial_{ij}\Phi\nabla^{2}N+\partial_{ij}N\nabla^{2}\Phi\Big{)}\Big{]}+\mathcal{O}(c^{-2})\,.\quad\,\,\,\,\end{split}$
(20c)

As pointed out, for a binary system, the expressions for the fields in the
near-zone are given in Ref. Julié and Berti (2019). In the above expressions
we only need to substitute for $N$ and $\Phi$ to the lowest PN order, which is
given by

$\begin{split}&N=4\frac{G(m_{A}+m_{B})}{rc^{2}}+\mathcal{O}(c^{-4})\,,\hskip
25.0pt\quad\Phi=-\frac{G(m_{A}\alpha^{0}_{A}+m_{B}\alpha^{0}_{B})}{rc^{2}}+\mathcal{O}(c^{-4})\,,\end{split}$
(21)

where $r$ is the relative separation between the masses such that
$\mathbf{r}^{i}=\mathbf{x}_{A}^{i}-\mathbf{x}_{B}^{i}$.

###  4.3 Evaluation of two-body Epstein-Wagoner moments

As can be seen from Eqs. (15) and (16), the EW moments are integrals over a
sphere of radius $\mathcal{R}$, about the CM of the system, with the variables
entering the integrands to be evaluated at retarded time $\tau$. We repeatedly
integrate by parts using the identity (see Ref. (Will and Wiseman, 1996), Eq.
(4.1) )

$\int_{\mathcal{M}}\partial_{k}F^{ij\ldots
m}d^{3}x=\left.\oint_{\partial\mathcal{M}}F^{ij\cdots
m}\right|_{\mathcal{R}}\hat{n}^{k}\mathcal{R}^{2}d^{2}\Omega\,.$ (22)

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 Eq. (18)

$\left(\delta^{ij}\right)_{TT}=0\,,\hskip
25.0pt\quad\left(\hat{n}^{i}F^{j}\right)_{TT}=0\,,$ (23)

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 Ref. (Will and Wiseman, 1996) and can be summarized to
the following steps: (i) use integration by parts to leave one potential
undifferentiated using Eq. (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\Omega$, (v) use the integration formula (Eq. (112)) for $f(\frac{1}{r})$ to
integrate over $dr$.
In the following calculations, we only report the calculations that are new to
the sGB theory, and we refer the reader to Ref. Will and Wiseman (1996) 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 Eq. 21
is such that $\Phi_{A}=-\alpha_{A}U_{A}$, with $U_{A}=Gm_{A}/rc^{2}$ the
Newtonian-order potential of body A. In the following calculations, we will
use this relation to substitute for $\Phi$ so that the source terms simplify
and that the functional form of the integrals mimic the usual GR terms.

####  4.3.1 Two-index moment $I^{ij}_{EW}$

We split the calculation into GB coupling dependent and independent terms,
such that $I^{ij}_{EW}=I^{ij}+I^{ij}_{GB}$.
The calculation of $I^{ij}$ 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 ($\alpha_{A}^{0}\Phi$) as well as the
$(\nabla\Phi)^{2}$ contribution, which follow the same calculation as the GR
ones, resulting in

$I^{ij}=m_{A}x_{A}^{ij}\left[1+\frac{1}{c^{2}}\left(\frac{v_{A}^{2}}{2}-\frac{Gm_{B}}{2r}-\frac{Gm_{B}\alpha^{0}_{A}\alpha^{0}_{B}}{2r}\right)\right]+(A\leftrightarrow
B)+{\mathcal{}O}(c^{-3}).$ (24)

For $I^{ij}_{GB}$ to 1PN order, substituting GB terms of Eq. (20a) in Eq. (15)
gives

$I^{ij}_{GB}\propto-\frac{c^{2}}{16\pi
G}\int_{\mathcal{M}}\partial_{ij}\Phi\partial_{ij}Nx^{ij}d^{3}x\,,$ (25)

where we have omitted the coupling factor for simplicity. We use integration-
by-parts several times on this expression and use Eq. (21) to simplify the
potentials. Ignoring constants and scalar-charge factors in front, we find

$\begin{split}I^{ij}_{GB}\propto\Big{\\{}&\oint_{\partial\mathcal{M}}\partial_{m}U\partial_{lm}Ux^{ij}dS^{l}-\oint_{\partial\mathcal{M}}\partial_{ll}U\partial_{m}Ux^{ij}dS^{m}+\int_{\mathcal{M}}\partial_{ll}U\partial_{m}U\partial_{m}(x^{ij})d^{3}x\\\
&-\int_{\mathcal{M}}\partial_{m}U\partial_{lm}U\partial_{l}(x^{ij})d^{3}x+\int_{\mathcal{M}}\partial_{ll}U\partial_{mm}Ux^{ij}d^{3}x\Big{\\}}\,,\end{split}$
(26)

where the Newtonian potential $U=\sum_{A}U_{A}=\sum_{A}Gm_{A}/rc^{2}$ is used
to simplify the expressions. It can be shown that the resulting surface
integrals vanish as they either do not have $\mathcal{R}$ independent
integrands or involve integration of a delta function
$\delta^{3}(\mathbf{x_{A}}-\mathbf{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
$\partial_{ll}U=\nabla^{2}U=-4\pi\sum_{A}m_{A}\delta^{3}(\mathbf{x}-\mathbf{x_{A}})$.
This results in a term proportional to
$m_{A}m_{B}\delta^{3}(\mathbf{x_{A}}-\mathbf{x_{B}}).$ As the calculations are
being done in the inspiral regime, where $\mathbf{x_{A}}\neq\mathbf{x_{B}}$ by
definition, this term is zero and thus, we are left with

$I^{ij}_{GB}\propto\frac{1}{Gc^{2}}\int_{\mathcal{M}}\partial_{ll}U\partial_{m}U\partial_{m}(x^{ij})d^{3}x-\frac{1}{Gc^{2}}\int_{\mathcal{M}}\partial_{m}U\partial_{lm}U\partial_{l}(x^{ij})d^{3}x\,.$
(27)

The first terms involves $\nabla^{2}U$ and can be readily evaluated. For the
second term, we further apply integration by parts and get

$\begin{split}-\int_{\mathcal{M}}\partial_{m}U\partial_{lm}U\partial_{l}(x^{ij})d^{3}x=&-\oint\partial_{l}U\partial_{m}U\partial_{l}(x^{ij})dS^{m}+\int\partial_{mm}U\partial_{l}U\partial_{l}(x^{ij})d^{3}x\\\
&+\int\partial_{m}U\partial_{l}U\partial_{ml}(x^{ij})d^{3}x\end{split}$ (28)

It can be shown again that the surface term vanishes and that the second term
is similar to the first term of Eq. (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

$\begin{split}\int\partial_{m}U\partial_{l}U\partial_{ml}(x^{ij})=\oint
U\partial_{l}U\partial_{ml}(x^{ij})dS^{m}-\int
U\partial_{lm}U\partial_{ml}(x^{ij})d^{3}x\,.\end{split}$ (29)

The surface integral does not have a $\mathcal{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

$I_{GB}^{ij}=-\frac{\alpha
f^{\prime}(\phi_{0})}{r^{2}}\frac{Gm_{A}m_{B}\alpha^{0}_{B}}{rc^{2}}\left[2r^{(i}x^{j)}_{A}+r^{ij}-2\delta^{ij}r^{2}\right]+(A\leftrightarrow
B)+{\mathcal{}O}(c^{-4})\,.$ (30)

Regarding the two-index surface moment $I^{ij}_{EW\operatorname{(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
$\mathcal{R}$-independent integrands.

Dropping the last term of Eq. (30) that does not have a TT contribution and
adding Eq. (24), the final expression for $I_{EW}^{ij}$ becomes

$\begin{split}I_{EW}^{ij}&=m_{A}x_{A}^{ij}\left[1+\frac{1}{c^{2}}\left(\frac{v_{A}^{2}}{2}-\frac{Gm_{B}(1+\alpha_{A}^{0}\alpha_{B}^{0})}{2r}\right)\right]-\frac{\alpha
f^{\prime}(\phi_{0})}{r^{2}}\frac{Gm_{A}m_{B}\alpha^{0}_{B}}{rc^{2}}\left(2r^{(i}x_{A}^{j)}+r^{ij}\right)\\\
&+(A\leftrightarrow B)+\mathcal{O}(c^{-4})\,.\end{split}$ (31)

An important point to note is the scaling of the GB terms. Since the terms
enter at ${\mathcal{}O}(c^{-2})$ in the PN expansion, they are a 1PN effect,
yet, due to the dimensionless ratio $\alpha/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
${\mathcal{}O}(c^{0})$ have the same scaling with separation $r$ as 5PN point-
mass terms Flanagan and Hinderer (2008). Calculations of tidal effects in
compact binaries involve a double expansion in PN and finite-size corrections
Flanagan (1998). Likewise, the GB term in Eq. (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.

####  4.3.2 Three-index moment $I^{ijk}_{EW}$

In order to obtain the waveform to relative 1PN order, it is sufficient to
calculate $I_{EW}^{ijk}$ to 0.5PN order. As there are no GB dependent terms at
this order (see Eq. (20b)), the three-index moment is the same as in GR

$I^{ijk}_{EW}=\frac{m_{A}}{c}\left(2v_{A}^{(i}x_{A}^{j)}x_{A}^{k}-v_{A}^{k}x_{A}^{ij}\right)+(A\leftrightarrow
B)+\mathcal{O}(c^{-3})\,.$ (32)

Also, the three-index surface moment $I_{EW\mathrm{(surf)}}^{ijk}$ is shown to
be zero in GR to 1PN order.

####  4.3.3 Four-index moment $I^{ijkl}_{EW}$

We again split the calculation into coupling dependent and independent terms,
namely $I_{EW}^{ijkl}=I^{ijkl}+I^{ijkl}_{GB}$.
The moment $I^{ijkl}$ has a GR part and a contribution from the scalar field
that is similar to the GR one (Will and Wiseman, 1996) up to a factor
difference. Including the factors, these terms together give

$I^{ijkl}=\frac{Gm_{A}m_{B}(1+\alpha_{A}^{0}\alpha_{B}^{0})r^{ij}}{12\,r\,c^{2}}\left(\frac{r^{kl}}{r^{2}}-\delta^{kl}-6\frac{x_{A}^{kl}}{r^{2}}\right)+\frac{m_{A}}{c^{2}}v_{A}^{ij}x_{A}^{kl}+(A\leftrightarrow
B)+{\mathcal{}O}(c^{-4})\,,$ (33)

where the terms proportional to $\delta^{ij}$ are dropped as they do not
produce TT components. For $I^{ijkl}_{GB}$ we have

$I^{ijmn}_{GB}\propto-\frac{c^{2}}{16\pi
G}\int_{\mathcal{M}}\Big{\\{}\Big{(}\partial_{ki}\Phi\partial_{kj}N-\partial_{kj}\Phi\partial_{ki}N\Big{)}+\Big{(}\partial_{ij}\Phi\nabla^{2}N+\partial_{ij}N\nabla^{2}\Phi\Big{)}\Big{\\}}x^{mn}d^{3}x\,.$
(34)

Using Eq. (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 Eq.
(21) to remove the integrals. Ignoring the $\delta^{ij}$ dependent terms that
do not give TT components we find

$\begin{split}I^{ijkl}_{EW}=&\frac{1}{c^{2}}m_{A}v_{A}^{ij}x_{A}^{kl}+\frac{Gm_{A}m_{B}(1+\alpha_{A}^{0}\alpha_{B}^{0})r^{ij}}{12\,r\,c^{2}}\left(\frac{r^{kl}}{r^{2}}-\delta^{kl}-\frac{6x_{A}^{kl}}{r^{2}}\right)\\\
&-\frac{3G\alpha
f^{\prime}(\phi_{0})m_{A}m_{B}\alpha^{0}_{B}r^{ij}}{c^{2}\,r^{5}}\,x^{kl}_{A}+(A\leftrightarrow
B)+{\mathcal{}O}(c^{-4}).\end{split}$ (35)

###  4.4 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 Eq. (12), can be written as

$\begin{split}&\Phi_{\mathcal{N}}(x)=\sum_{l=0}^{\infty}\Phi_{l}(t,\textbf{x})\,,\hskip
25.0pt\Phi_{l}(x)=\frac{G}{R\,c^{4}}\frac{\hat{n}_{L}}{l!}\left(\frac{\partial}{c\,\partial
t}\right)^{l}I^{L}_{s}+\mathcal{O}(R^{-2})\,,\end{split}$ (36)

with the scalar multipole moments defined as

$I_{s}^{L}=\int_{\mathcal{M}}\mu_{s}(\tau,x^{\prime})x^{\prime
L}d^{3}x^{\prime}\,.$ (37)

It is important at this stage to discuss the counting of PN orders. By
comparing Eq. (36) to Eq. (14), we see that the lowest order term in the
tensor wave, containing the quadrupole moment, is of order
$\mathcal{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 $\mathcal{R}^{2}_{GB}$ (see Appendix A.3, Eq. (104)) to
leading order in the weak-field limit, we are limited to having $\mu_{s}$ to
$\mathcal{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 Eq. (19), it can be easily shown that the PN
expansion of $\mathcal{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

$\begin{split}\mu_{s}=&+m_{A}c^{2}\,\alpha_{A}^{0}\delta^{3}(\mathbf{x}-\mathbf{x}_{A}(t))\left[1-\frac{v^{2}_{A}}{2c^{2}}-\frac{3N}{4}+\left(\alpha_{A}^{0}+\frac{\beta_{A}^{0}}{\alpha_{A}^{0}}\right)\Phi\right]\\\
&-\frac{\alpha f^{\prime}(\phi_{0})c^{4}}{32\pi
G}\Big{[}(\partial_{kl}N)^{2}-(\partial_{kk}N)^{2}\Big{]}+(A\leftrightarrow
B)+\mathcal{O}(c^{-2})\,.\end{split}$ (38)

With this result, we can now calculate the scalar multipole moments of Eq. 37.
We define $I_{s}^{L}=I_{s,c}^{L}+I_{s,GB}^{L}$ 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.5PN GB
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

$I_{s,GB}^{i}=-\frac{\alpha f^{\prime}(\phi_{0})c^{4}}{32\pi
G}\int_{\mathcal{M}}\Big{[}(\nabla^{2}N)^{2}-\partial_{i}N\nabla^{2}N-\partial_{l}N\partial_{il}N\Big{]}d^{3}x\,.$
(39)

The non-vanishing part of the above equation involves the integral
$\int(\nabla N)^{2}d^{3}x$. As already discussed (see discussion below Eq.
(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_{s,c}^{L}$ are terms depending on
the matter sources. It is straightforward to see that:

$\begin{split}I_{s,c}&=m_{A}c^{2}\alpha_{A}^{0}\Big{\\{}1-\frac{v^{2}_{A}}{2c^{2}}-\frac{Gm_{B}\alpha_{AB}}{r\,c^{2}}+\mathcal{O}(c^{-4})\Big{\\}}+(A\leftrightarrow
B),\\\
I_{s,c}^{i}&=x^{i}_{A}m_{A}c^{2}\alpha_{A}^{0}\Big{\\{}1-\frac{v^{2}_{A}}{2c^{2}}-\frac{Gm_{B}\alpha_{AB}}{r\,c^{2}}+\mathcal{O}(c^{-4})\Big{\\}}\
+(A\leftrightarrow B)\,,\\\
I_{s,c}^{ij}&=x^{ij}_{A}m_{A}c^{2}\alpha_{A}^{0}\\{1+\mathcal{O}(c^{-2})\\}+(A\leftrightarrow
B)\,,\\\
I_{s,c}^{ijk}&=x^{ijk}_{A}m_{A}c^{2}\alpha_{A}^{0}\\{1+\mathcal{O}(c^{-2})\\}+(A\leftrightarrow
B)\,,\end{split}$ (40)

where we have defined
$\alpha_{AB}=(1+\alpha_{A}\alpha_{B}+\beta_{A}\alpha_{B}\alpha_{A}^{-1})$ .

##  5 Two-body waveforms in the center of mass frame

####  5.0.1 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 Ref. (Julié and Berti, 2019), where the authors
constructed the Fokker Lagrangian using the near-zone solution to the field
equations [Eq. (3)], finding

$\begin{split}&L_{AB}=-m_{A}c^{2}+\frac{1}{2}m_{A}\mathbf{v}_{A}^{2}+\frac{G\bar{\alpha}\,m_{A}m_{B}}{2r}+\frac{1}{8c^{2}}m_{A}\mathbf{v}_{A}^{4}\\\
&+\frac{G\bar{\alpha}\,m_{A}m_{B}}{rc^{2}}\left[-\frac{G\bar{\alpha}m_{A}}{2r}\left(1+2\bar{\beta}_{B}\right)+\frac{3}{2}\left(\mathbf{v}_{A}^{2}\right)-\right.\\\
&\left.\frac{7}{4}\left(\mathbf{v}_{A}\cdot\mathbf{v}_{B}\right)-\frac{1}{4}\left(\hat{\mathbf{n}}_{AB}\cdot\mathbf{v}_{A}\right)\left(\hat{\mathbf{n}}_{AB}\cdot\mathbf{v}_{B}\right)+\frac{\bar{\gamma}}{2}\left(\mathbf{v}_{A}-\mathbf{v}_{B}\right)^{2}\right]\\\
&+\frac{\alpha
f^{\prime}\left(\varphi_{0}\right)}{r^{2}}\frac{G^{2}m_{A}m_{B}}{r^{2}c^{2}}\left[m_{A}\left(\alpha_{B}^{0}+2\alpha_{A}^{0}\right)\right]+(A\leftrightarrow
B)\,.\end{split}$ (41)

As before, $r=\|\mathbf{r}\|=\|\mathbf{x}_{A}-\mathbf{x}_{B}\|$ so that
$\hat{\mathbf{n}}_{AB}=\mathbf{r}/r$ and in anologous to ST theories Damour
and Esposito-Farese (1992a), we define the binary parameters

$\begin{split}\bar{\alpha}\equiv\left(1+\alpha_{A}^{0}\alpha_{B}^{0}\right)\,,\quad\bar{\gamma}\equiv-2\frac{\alpha_{A}^{0}\alpha_{B}^{0}}{\bar{\alpha}}\,,\quad\bar{\beta}_{A}\equiv\frac{1}{2}\frac{\beta_{A}^{0}(\alpha_{B}^{0})^{2}}{\bar{\alpha}^{2}}.\end{split}$
(42)

The parameters $\bar{\gamma}$ and $\bar{\beta}$ are generalizations of 1PN
Eddington parameters, and $\bar{\alpha}$ appears as a re-scaling factor for
the gravitational constant $G$. It is useful to also introduce the
combinations

$\mathcal{S}_{\pm}=\frac{\alpha_{A}^{0}\pm\alpha_{B}^{0}}{2\sqrt{\bar{\alpha}}},\hskip
25.0pt\hskip 25.0pt\beta_{\pm}=\frac{\bar{\beta}_{A}\pm\bar{\beta}_{B}}{2}\,.$
(43)

We see that the Lagrangian $L_{AB}$ decomposes into a piece that is
structurally equivalent to the two-body Lagrangian in massless ST theories
(Damour and Esposito-Farese, 1992b) plus a distinct GB-coupling dependent
contribution. The effects of the scalar field are thus entirely contained in
the charge-dependent binary parameters of Eq. (42). In massless ST theories,
the gravitational constant $G$ always appears in combination with ST analogous
parameter of $\bar{\alpha}$, which implies that the effective gravitational
coupling in these theories is $G\bar{\alpha}$. However in sGB gravity, the GB
higher curvature terms break this rescaling, as made explicit in the last line
of Eq. (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

$\begin{split}&X^{i}_{CM}=\frac{1}{m}\int\mu^{00}r^{i}d^{3}\mathbf{r}\,.\\\
\end{split}$ (44)

By computing the above integral with the same methods as in Sec. 4.3, we find
that, to 1PN order, the coordinates of each body in the CM frame are

$\begin{split}&\mathbf{x}_{A}=\left[\frac{m_{B}}{m}+\frac{\mu\Delta
m}{2m^{2}c^{2}}\left(v^{2}-\frac{Gm\bar{\alpha}}{r}\right)\right]\mathbf{r}+\delta+{\mathcal{}O}(c^{-3})\,,\\\
&\mathbf{x}_{B}=\left[-\frac{m_{A}}{m}+\frac{\mu\Delta
m}{2m^{2}c^{2}}\left(v^{2}-\frac{Gm\bar{\alpha}}{r}\right)\right]\mathbf{r}+\delta+{\mathcal{}O}(c^{-3}),\end{split}$
(45)

where $m=m_{A}+m_{B}$ is total mass of the system,
$\eta=\mu/m=m_{A}m_{B}/m^{2}$ is the symmetric mass ratio,
$\mathbf{v}=\mathbf{v}_{A}-\mathbf{v}_{B}$ is the relative velocity, $\Delta
m=m_{A}-m_{B}$, and $\dot{r}=(\mathbf{v}\cdot\mathbf{r})/r$.
The GB coupling dependent term is given by

$\delta=-2\eta\left(\frac{Gm\bar{\alpha}}{r\,c^{2}}\frac{\alpha
f^{\prime}(\phi_{0})}{\sqrt{\bar{\alpha}}r^{2}}\mathcal{S}_{+}\right)\mathbf{r}\,.$
(46)

Taking a time derivative of Eq. (45), we obtain the velocities

$\begin{split}\mathbf{v}_{A}&=\frac{m_{B}}{m}\mathbf{v}+\frac{\mu\Delta
m}{2m^{2}c^{2}}\left[\left(v^{2}-\frac{Gm\bar{\alpha}}{r}\right)\mathbf{v}-\frac{Gm\bar{\alpha}}{r^{2}}\dot{r}\mathbf{r}\right]+\dot{\delta}+{\mathcal{}O}(c^{-3})\,,\\\
\mathbf{v}_{B}&=-\frac{m_{A}}{m}\mathbf{v}+\frac{\mu\Delta
m}{2m^{2}\,c^{2}}\left[\left(v^{2}-\frac{Gm\bar{\alpha}}{r}\right)\mathbf{v}-\frac{Gm\bar{\alpha}}{r^{2}}\dot{r}\mathbf{r}\right]+\dot{\delta}+{\mathcal{}O}(c^{-3})\,,\end{split}$
(47)

with the GB coupling dependent contribution being

$\dot{\delta}=2\eta\frac{Gm\bar{\alpha}}{rc^{2}}\frac{\alpha
f^{\prime}(\phi_{0})}{\sqrt{\bar{\alpha}}r^{2}}\mathcal{S}_{+}(3\dot{r}\mathbf{r}-\mathbf{v})\,.$
(48)

We substitute these relations in the Euler-Lagrange equations corresponding to
Eq. (41), and find the 1PN relative matter equation of motion

$\begin{split}\frac{d^{2}\mathbf{r}}{dt^{2}}=&\bm{a}=-\frac{G\bar{\alpha}m}{r^{2}}\hat{\mathbf{n}}_{AB}+\frac{G\bar{\alpha}m}{r^{2}\,c^{2}}\Big{\\{}2\dot{r}\mathbf{v}\Big{[}2-\eta+\bar{\gamma}\Big{]}+\Big{[}-(1+3\eta+\bar{\gamma})\mathbf{v}^{2}+\frac{3}{2}\eta\dot{r}^{2}\\\
&+\frac{2G\bar{\alpha}m}{r}\Big{(}2+\eta+\bar{\gamma}+\beta_{+}-\frac{\Delta
m}{m}\beta_{-}+\frac{2\alpha
f^{\prime}(\phi_{0})}{\bar{\alpha}^{3/2}r^{2}}(3\mathcal{S}_{+}+\frac{\Delta
m}{m}\mathcal{S}_{-})\Big{)}\Big{]}\hat{\mathbf{n}}_{AB}\Big{\\}}+\mathcal{O}(c^{-3})\,,\end{split}$
(49)

From the Lagrangian, we can also read the conserved binding energy of the
system to 1PN order,

$\begin{split}E&=m_{A}\left(\frac{v_{A}^{2}}{2}+\frac{3v_{A}^{4}}{8c^{2}}\right)+m_{B}\left(\frac{v_{B}^{2}}{2}+\frac{3v_{B}^{4}}{8c^{2}}\right)-\frac{G\mu
m\bar{\alpha}}{r}\\\ &+\frac{G\mu
m\bar{\alpha}}{rc^{2}}\left[\frac{3+2\bar{\gamma}}{2}(v_{A}^{2}+v_{B}^{2})-\frac{1}{2}(v_{A}\cdot
n)(v_{B}\cdot n)-\frac{7+4\bar{\gamma}}{2}(v_{A}\cdot v_{B})\right]\\\
&+\frac{\mu}{c^{2}}\left(\frac{Gm\bar{\alpha}}{r}\right)^{2}\left[\frac{1}{2}+\beta_{+}-\frac{\Delta
m}{m}\beta_{-}\right]-\frac{\alpha
f^{\prime}(\phi_{0})\mu}{\bar{\alpha}^{3/2}r^{2}c^{2}}\left(\frac{Gm\bar{\alpha}}{r}\right)^{2}\left(3\mathcal{S}_{+}+\frac{\Delta
m}{m}\mathcal{S}_{-}\right)\,.\end{split}$ (50)

In the CM frame, we have

$\begin{split}E=&\mu\Bigg{\\{}\frac{1}{2}v^{2}-\frac{G\bar{\alpha}m}{r}+\frac{3}{8c^{2}}(1-3\eta)v^{4}+\frac{G\bar{\alpha}m}{2rc^{2}}\left[(3+2\bar{\gamma}+\eta)v^{2}+\eta\dot{r}^{2}\right]\\\
&+\frac{1}{c^{2}}\left(\frac{G\bar{\alpha}m}{r}\right)^{2}\left[\frac{1}{2}+\beta_{+}-\frac{\Delta
m}{m}\beta_{-}-\frac{\alpha
f^{\prime}(\phi_{0})}{\bar{\alpha}^{3/2}r^{2}}\left(3\mathcal{S}_{+}+\frac{\Delta
m}{m}\mathcal{S}_{-}\right)\right]\Bigg{\\}}\,.\end{split}$ (51)

####  5.0.2 Final expressions for the waveform

By substituting the expressions for EW moments (see Eqs. (31), (32), and (35))
in Eq. (17), and simplifying further using Eqs. (45) and (47), we find to
$\mathcal{O}(c^{-3})$ that

$\begin{split}&h^{ij}_{TT}=\frac{2G\mu}{Rc^{4}}\frac{d^{2}}{dt^{2}}\Bigg{\\{}\left[1+\frac{1}{c^{2}}\left(1-\frac{3\mu}{m}\right)\tilde{E}+\frac{Gm\bar{\alpha}}{3r\,c^{2}}\left(1-\frac{6\mu}{m}\right)\right]r^{ij}\\\
&-\frac{\Delta
m}{m\,c^{2}}\left(2v^{(i}r^{j)}(\hat{\mathbf{n}}\cdot\mathbf{r})-(\hat{\mathbf{n}}\cdot\mathbf{v})r^{ij}\right)+\frac{1}{c^{2}}\left(1-\frac{3\mu}{m}\right)(\mathbf{r}\cdot\hat{\mathbf{n}})^{2}\left(v^{ij}-\frac{Gm\bar{\alpha}r^{ij}}{3r^{3}}\right)\\\
&-\frac{\alpha
f^{\prime}(\phi_{0})}{\sqrt{\bar{\alpha}}r^{2}}\frac{G\bar{\alpha}mr^{ij}}{rc^{2}}\left[4\left(\mathcal{S}_{+}+\frac{\Delta
m}{m}\mathcal{S}_{-}\right)+3\big{[}\mathcal{S}_{+}(1-2\eta)+\mathcal{S}_{-}(1+2\eta)\big{]}\frac{(\mathbf{r}\cdot\hat{\mathbf{n}})^{2}}{r^{2}}\right]\Bigg{\\}}_{TT},\end{split}$
(52)

with $\tilde{E}=v^{2}/2-Gm\bar{\alpha}/r+\mathcal{O}(c^{-2})$ being the
orbital binding energy per reduced mass at Newtonian order.
Taking time derivatives and perturbatively using the relative acceleration
when needed (see Eq.(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

$\begin{split}&h^{ij}_{{}_{TT}}=\frac{2G\mu}{Rc^{4}}\\{Q^{ij}\\}_{{}_{TT}}=\frac{2G\mu}{Rc^{4}}\left\\{\tilde{Q}^{ij}+\frac{1}{c}P^{1/2}\tilde{Q}^{ij}+\frac{1}{c^{2}}\left(P\tilde{Q}^{ij}+P\tilde{Q}^{ij}_{GB}\right)+{\mathcal{}O}(c^{-3})\right\\}_{{}_{TT}}\,,\\\
\\\
&\tilde{Q}^{ij}=2\left[v^{ij}-\frac{Gm\bar{\alpha}r^{ij}}{r^{3}}\right]\,,\\\
\\\ &P^{1/2}\tilde{Q}^{ij}=\frac{\Delta
m}{m}\left[3\frac{Gm\bar{\alpha}}{r^{3}}(\hat{\mathbf{n}}\cdot\mathbf{r})\left(2v^{(i}r^{j)}-\frac{\dot{r}r^{ij}}{r}\right)-(\hat{\mathbf{n}}\cdot\mathbf{v})\left(2v^{ij}-\frac{Gm\bar{\alpha}r^{ij}}{r^{3}}\right)\right]\,,\\\
\\\
&P\tilde{Q}^{ij}=\frac{1-3\eta}{3}\Bigg{\\{}(\mathbf{r}\cdot\hat{\mathbf{n}})^{2}\,\frac{G\bar{\alpha}m}{r^{3}}\left[\left(6\bar{E}-15\dot{r}^{2}+13\frac{G\bar{\alpha}m}{r}\right)\frac{r^{ij}}{r^{2}}+30\dot{r}\frac{r^{(i}v^{j)}}{r}-14v^{ij}\right]\\\
&\quad+(\hat{\mathbf{n}}\cdot\mathbf{v})^{2}\left[6v^{ij}-2\frac{G\bar{\alpha}m}{r^{3}}r^{ij}\right]+(\mathbf{r}\cdot\hat{\mathbf{n}})(\hat{\mathbf{n}}\cdot\mathbf{v})\frac{G\bar{\alpha}m}{r^{2}}\left[12\frac{\dot{r}r^{ij}}{r^{2}}-32\frac{r^{(i}v^{j)}}{r}\right]\Bigg{\\}}\\\
&\quad+\frac{1}{3}\left\\{\left[3(1-3\eta)v^{2}-2(2-3\eta)\frac{G\bar{\alpha}m}{r}\right]v^{ij}+4\frac{G\bar{\alpha}m}{r}(5+3\eta+3\bar{\gamma})\frac{\dot{r}}{r}r^{(i}v^{j)}\right.\\\
&\quad\left.+\frac{G\bar{\alpha}m}{r^{3}}r^{ij}\left[3(1-3\eta)\dot{r}^{2}-(10+3\eta+6\bar{\gamma})v^{2}+\left(29+12\bar{\gamma}+12\beta_{+}-12\frac{\Delta
m}{m}\beta_{-}\right)\frac{G\bar{\alpha}m}{r}\right]\right\\}\,,\\\ \\\
&P\tilde{Q}_{GB}^{ij}=4\frac{G\bar{\alpha}m}{r}\frac{\alpha
f^{\prime}(\phi_{0})}{\sqrt{\bar{\alpha}}r^{2}}\left[\left(\mathcal{S}_{+}+\frac{\Delta
m}{m}\mathcal{S}_{-}\right)\left(\frac{6\dot{r}r^{(i}v^{j)}}{r}+\frac{r^{ij}}{r^{2}}\left(6\tilde{E}+\frac{5G\bar{\alpha}m}{r}-15\dot{r}^{2}\right)\right.\right.\\\
&\quad\left.\left.-2v^{ij}\right)+2\frac{G\bar{\alpha}m}{r}\frac{r^{ij}}{r^{2}}\left(3\mathcal{S}_{+}+\frac{\Delta
m}{m}\mathcal{S}_{-}\right)\right]+3\frac{G\bar{\alpha}m}{r}\frac{\alpha
f^{\prime}(\phi_{0})}{\sqrt{\bar{\alpha}}r^{4}}\bigg{(}\mathcal{S}_{+}(1-2\eta)+\mathcal{S}_{-}(1+2\eta)\bigg{)}\\\
&\quad\Bigg{\\{}(\mathbf{r}\cdot\hat{\mathbf{n}})^{2}\left[\frac{r^{ij}}{r^{2}}\left(10\tilde{E}-35\dot{r}^{2}+\frac{9G\bar{\alpha}m}{r}\right)+\frac{10\dot{r}r^{(i}v^{j)}}{r}-2v^{ij}\right]-2(\hat{\mathbf{n}}\cdot\mathbf{v})^{2}r^{ij}\\\
&\hskip
25.0pt+(\mathbf{r}\cdot\hat{\mathbf{n}})(\hat{\mathbf{n}}\cdot\mathbf{v})\left[20\frac{\dot{r}r^{ij}}{r}-4r^{(i}v^{j)}\right]\Bigg{\\}}\,,\end{split}$
(53)

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
$\Phi_{l}$, including up to 0.5PN terms. By rewriting scalar multipoles (Eq.
(40)) in the CM frame using Eqs. (45) - (47), and taking time derivatives, we
find

$\begin{split}&\Phi_{0}=\frac{2Gm\sqrt{\bar{\alpha}}}{Rc^{2}}\left(\mathcal{S}_{+}+\frac{\Delta
m}{m}\mathcal{S}_{-}\right)\\\
&\quad+\frac{2G\mu\sqrt{\bar{\alpha}}}{Rc^{4}}\left\\{-\frac{v^{2}}{2}\left(\mathcal{S}_{+}-\frac{\Delta
m}{m}\mathcal{S}_{-}\right)+\left(\frac{8}{\bar{\gamma}}\left(\mathcal{S}_{+}\beta_{+}+\mathcal{S}_{-}\beta_{-}\right)-2\mathcal{S}_{+}\right)\frac{G\bar{\alpha}m}{r}\right\\}\,,\\\
&\Phi_{1}=\frac{2G\mu\sqrt{\bar{\alpha}}}{Rc^{3}}\Bigg{\\{}(\hat{\mathbf{n}}\cdot\mathbf{v})\Bigg{[}2\mathcal{S}_{-}+\frac{v^{2}}{c^{2}}\left(\frac{\Delta
m}{m}\mathcal{S}_{+}-\eta\mathcal{S}_{-}\right)\\\
&\quad+\frac{G\bar{\alpha}m}{rc^{2}}\Bigg{(}\frac{\Delta
m}{2m}\mathcal{S}_{+}+\left(2\eta-\frac{3}{2}\right)\mathcal{S}_{-}-\frac{4}{\bar{\gamma}}\frac{\Delta
m}{m}\left(\mathcal{S}_{+}\beta_{+}+\mathcal{S}_{-}\beta_{-}\right)+\frac{4}{\bar{\gamma}}\left(\mathcal{S}_{-}\beta_{+}+\mathcal{S}_{+}\beta_{-}\right)\\\
&\quad-2\eta\mathcal{S}_{+}\frac{\Delta
m}{m}\left(\mathcal{S}_{+}+\frac{\Delta
m}{m}\mathcal{S}_{-}\right)\frac{\alpha
f^{\prime}(\phi_{0})}{\sqrt{\bar{\alpha}}r^{2}}\Bigg{)}\Bigg{]}+\frac{G\bar{\alpha}m}{c^{3}r^{2}}\dot{r}(\hat{\mathbf{n}}\cdot\mathbf{r})\Bigg{[}-\frac{5}{2}\frac{\Delta
m}{m}\mathcal{S}_{+}+\frac{3}{2}\mathcal{S}_{-}\\\
&\quad+\frac{4}{\bar{\gamma}}\frac{\Delta
m}{m}\left(\mathcal{S}_{+}\beta_{+}+\mathcal{S}_{-}\beta_{-}\right)-\frac{4}{\bar{\gamma}}\left(\mathcal{S}_{-}\beta_{+}+\mathcal{S}_{+}\beta_{-}\right)+6\eta\mathcal{S}_{+}\frac{\Delta
m}{m}\left(\mathcal{S}_{+}+\frac{\Delta
m}{m}\mathcal{S}_{-}\right)\frac{\alpha
f^{\prime}(\phi_{0})}{\sqrt{\bar{\alpha}}r^{2}}\Bigg{]}\Bigg{\\}}\,,\\\
&\Phi_{2}=\frac{2G\mu\sqrt{\bar{\alpha}}}{Rc^{4}}\left(\mathcal{S}_{+}-\frac{\Delta
m}{m}\mathcal{S}_{-}\right)\left\\{(\hat{\mathbf{n}}\cdot\mathbf{v})^{2}-\frac{G\bar{\alpha}m}{r}\left(\frac{\hat{\mathbf{n}}\cdot\mathbf{r}}{r}\right)^{2}\right\\}\,,\\\
&\Phi_{3}=\frac{2G\mu\sqrt{\bar{\alpha}}}{Rc^{5}}\left((1-2\eta)\mathcal{S}_{-}-\frac{\Delta
m}{m}\mathcal{S}_{+}\right)\left\\{\frac{3}{2}\frac{G\bar{\alpha}m}{r^{4}}\dot{r}(\hat{\mathbf{n}}\cdot\mathbf{r})^{3}-\frac{7}{2}\frac{G\bar{\alpha}m}{r^{3}}(\hat{\mathbf{n}}\cdot\mathbf{v})(\hat{\mathbf{n}}\cdot\mathbf{r})^{2}+(\hat{\mathbf{n}}\cdot\mathbf{v})^{3}\right\\}\,.\end{split}$
(54)

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 equal-mass binaries,
for which the dipole radiation is suppressed. We will come back to this point
in Sec. 7.4, where we compare the PN scalar waveform with scalar waves from
numerical relativity 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

$\begin{split}&\Phi=\frac{2G\mu\sqrt{\bar{\alpha}}}{Rc^{3}}\Big{\\{}P^{-1/2}\tilde{\Phi}+\frac{1}{c}\tilde{\Phi}+\frac{1}{c^{2}}P^{1/2}\tilde{\Phi}+{\mathcal{}O}(c^{-3})\Big{\\}}\,,\\\
&P^{-1/2}\tilde{\Phi}=2\mathcal{S}_{-}(\hat{\mathbf{n}}\cdot\mathbf{v})\,,\\\
&\tilde{\Phi}=\left(\mathcal{S}_{+}-\frac{\Delta
m}{m}\mathcal{S}_{-}\right)\left[-\frac{G\bar{\alpha}m}{r}\left(\frac{\hat{\mathbf{n}}\cdot\mathbf{r}}{r}\right)^{2}+(\hat{\mathbf{n}}\cdot\mathbf{v})^{2}-\frac{1}{2}v^{2}\right]\\\
&\quad+\frac{G\bar{\alpha}m}{r}\left[-2\mathcal{S}_{+}+\frac{8}{\bar{\gamma}}\left(\mathcal{S}_{+}\beta_{+}+\mathcal{S}_{-}\beta_{-}\right)\right]\,,\\\
&P^{1/2}\tilde{\Phi}=\left(-\frac{\Delta
m}{m}\mathcal{S}_{+}+(1-2\eta)\mathcal{S}_{-}\right)\left[\frac{3}{2}\frac{G\bar{\alpha}m}{r^{4}}\dot{r}(\hat{\mathbf{n}}\cdot\mathbf{r})^{3}-\frac{7}{2}\frac{G\bar{\alpha}m}{r^{3}}(\hat{\mathbf{n}}\cdot\mathbf{v})(\hat{\mathbf{n}}\cdot\mathbf{r})^{2}+(\hat{\mathbf{n}}\cdot\mathbf{v})^{3}\right]\\\
&\quad+(\hat{\mathbf{n}}\cdot\mathbf{v})\left\\{\left(\frac{\Delta
m}{m}\mathcal{S}_{+}-\eta\mathcal{S}_{-}\right)v^{2}+\frac{G\bar{\alpha}m}{r}\left[\frac{1}{2}\frac{\Delta
m}{m}\mathcal{S}_{+}+\left(2\eta-\frac{3}{2}\right)\mathcal{S}_{-}\right.\right.\\\
&\quad\left.\left.-\frac{4}{\bar{\gamma}}\frac{\Delta
m}{m}\left(\mathcal{S}_{+}\beta_{+}+\mathcal{S}_{-}\beta_{-}\right)+\frac{4}{\bar{\gamma}}\left(\mathcal{S}_{-}\beta_{+}+\mathcal{S}_{+}\beta_{-}\right)\right]\right\\}\\\
&\quad+\frac{G\bar{\alpha}m}{r^{2}}\dot{r}(\hat{\mathbf{n}}\cdot\mathbf{r})\left[\frac{3}{2}\mathcal{S}_{-}-\frac{5}{2}\frac{\Delta
m}{m}\mathcal{S}_{+}+\frac{4}{\bar{\gamma}}\frac{\Delta
m}{m}\left(\mathcal{S}_{+}\beta_{+}+\mathcal{S}_{-}\beta_{-}\right)-\frac{4}{\bar{\gamma}}\left(\mathcal{S}_{-}\beta_{+}+\mathcal{S}_{+}\beta_{-}\right)\right]\\\
&\quad+2\eta\frac{\Delta m}{m}\frac{G\bar{\alpha}m}{r}\frac{\alpha
f^{\prime}(\phi_{0})}{\sqrt{\bar{\alpha}}r^{2}}\mathcal{S}_{+}\left(\mathcal{S}_{+}+\frac{\Delta
m}{m}\mathcal{S}_{-}\right)\left[3\frac{\dot{r}}{r}(\hat{\mathbf{n}}\cdot\mathbf{r})-(\hat{\mathbf{n}}\cdot\mathbf{v})\right]\,,\end{split}$
(55)

For both the scalar and tensor waveforms, we obtain similar results in
structure to those of ST gravity Lang (2014, 2015) (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 $\bar{\alpha}$ multiplying the
gravitational constant. At 1PN order, we see dependency on the new set of
parameters and explicit GB coupling-dependent terms.

In Sec. 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 Eq. (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 Eq. (49)) at 1PN
order. Also note that, in the GB coupling-dependent terms of the gravitational
waveform, we see dependency on the parameters $\mathcal{S}_{+}$ and
$\mathcal{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 radiation-reaction terms in the
matter equation of motion.

##  6 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 Sec. 7, the energy flux
rate is used to determine the phase evolution of GWs.

###  6.1 Tensor mode

The energy loss from tensor waves is given by

$\dot{E}_{T}=\frac{c^{3}\,R^{2}}{32\pi
G}\oint\dot{h}_{\mathrm{TT}}^{ij}\dot{h}_{\mathrm{TT}}^{ij}d^{2}\Omega\,.$
(56)

We simplify this expression starting from the effects of TT projection
operators. Based on the definition given in Eq. (18), it is straightforward to
show that

$\begin{split}\left(P^{ik}P^{jl}-\frac{1}{2}P^{ij}P^{kl}\right)&\left(P^{im}P^{jn}-\frac{1}{2}P^{ij}P^{mn}\right)=P^{km}P^{ln}-\frac{1}{2}P^{kl}P^{mn}\,.\end{split}$
(57)

Using this identity, the calculation simplifies to

$\begin{split}\dot{E}_{T}=\frac{\mu^{2}c^{3}}{32\pi
G}\oint&\big{(}4\dot{Q}^{ij}\dot{Q}^{ij}-8n^{ln}\dot{Q}^{kl}\dot{Q}^{kn}+2n^{klmn}\dot{Q}^{kl}\dot{Q}^{mn}\big{)}d^{2}\Omega\,,\end{split}$
(58)

where we use the same notation as introduced in the introduction. Taking a
time derivative and evaluating angular integrals using Eq.(107), we find that

$\begin{split}\dot{E}_{T}&=\frac{8}{15}\frac{\eta^{2}}{G\bar{\alpha}^{2}c^{5}\,}\left(\frac{G\bar{\alpha}m}{r}\right)^{4}\Bigg{\\{}(12v^{2}-11\dot{r}^{2})\\\
&+\frac{1}{28c^{2}}\bigg{[}-16\left(170-10\eta+63\bar{\gamma}+84\beta_{+}-84\frac{\Delta
m}{m}\beta_{-}\right)v^{2}\frac{G\bar{\alpha}m}{r}\\\
&+(785-852\eta+336\bar{\gamma})v^{4}-2(1487-1392\eta+616\bar{\gamma})v^{2}\dot{r}^{2}+3(687-620\eta+280\bar{\gamma})\dot{r}^{4}\\\
&+8\left(367-15\eta+140\bar{\gamma}+168\beta_{+}-168\frac{\Delta
m}{m}\beta_{-}\right)\dot{r}^{2}\frac{G\bar{\alpha}m}{r}+16(1-4\eta)\left(\frac{G\bar{\alpha}m}{r}\right)^{2}\bigg{]}\\\
&+\frac{3f^{\prime}(\phi_{0})\alpha}{4\sqrt{\bar{\alpha}}r^{2}c^{2}}\bigg{[}4\left(\mathcal{S}_{+}+\frac{\Delta
m}{m}\mathcal{S}_{-}\right)\Big{(}-4v^{2}\left(18\tilde{E}+13\frac{G\bar{\alpha}m}{r}-45\dot{r}^{2}\right)\\\
&-\dot{r}^{2}\left(108\tilde{E}+85\frac{G\bar{\alpha}m}{r}-150\dot{r}^{2}\right)+54\dot{r}^{4}\Big{)}\frac{\bar{G\alpha}m}{r}\left(3\mathcal{S}_{+}+\frac{\Delta
m}{m}\mathcal{S}_{-}\right)\left(32v^{2}+56\dot{r}^{2}\right)\bigg{]}\\\
&+\frac{45f^{\prime}(\phi_{0})\alpha}{8\sqrt{\bar{\alpha}}r^{2}\,c^{2}}\left[\mathcal{S}_{+}(1-2\eta)+\mathcal{S}_{-}(1+2\eta)\right]\bigg{[}-\frac{4}{7}v^{2}\left(22\tilde{E}+77\frac{G\bar{\alpha}m}{r}-\frac{199}{3}\dot{r}^{2}+\frac{18}{5}v^{2}\right)\\\
&+\frac{24}{5}\dot{r}^{4}-\frac{16}{5}v^{2}\dot{r}^{2}-\frac{\dot{r}^{2}}{7}\left(992\tilde{E}+737\frac{G\bar{\alpha}m}{r}-\frac{2404}{3}\dot{r}^{2}+\frac{8}{5}v^{2}\right)\bigg{]}+\mathcal{O}(c^{-3})\Bigg{\\}}\,.\end{split}$
(59)

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 $0\mathrm{PN}$ piece of $h_{TT}^{ij}$ by itself. At
$0.5\mathrm{PN}$ order, there is no tensor flux as the product of
$0.5\mathrm{PN}$ and $0\mathrm{PN}$ terms has an odd number of unit vectors
and thus vanishes upon angular integration. The remaining terms are at
$1\mathrm{PN}$ order, comprising $0\mathrm{PN}-1\mathrm{PN}$, and
$0.5\mathrm{PN}-0.5\mathrm{PN}$ products.

###  6.2 Scalar mode

For the scalar field, the energy loss is evaluated from

$\dot{E}_{s}=\frac{c^{3}R^{2}}{32\pi G}\oint\dot{\Phi}^{2}d^{2}\Omega\,.$ (60)

Since we have defined the lowest-order tensor flux term as a $0\mathrm{PN}$
term, the lowest-order piece of the scalar flux is a -1PN term, resulting from
multiplying the -0.5 PN piece of $\Phi$ by itself. We find the scalar energy
flux to 0PN order to be

$\begin{split}&\dot{E_{S}}=\frac{\eta^{2}}{G\bar{\alpha}c^{3}}\left(\frac{G\bar{\alpha}m}{r}\right)^{4}\Bigg{\\{}\frac{4}{3}\mathcal{S}_{-}^{2}+\frac{8}{15\,c^{2}}\left(\frac{G\bar{\alpha}m}{r}\left[\left(-23+\eta-10\bar{\gamma}-10\beta_{+}+10\frac{\Delta
m}{m}\beta_{-}\right)\mathcal{S}_{-}^{2}\right.\right.\\\
&\left.\left.-2\frac{\Delta
m}{m}\mathcal{S}_{+}\mathcal{S}_{-}\right]\right.+v^{2}\left[2\mathcal{S}_{+}^{2}+2\frac{\Delta
m}{m}\mathcal{S}_{+}\mathcal{S}_{-}+(6-\eta+5\bar{\gamma})\mathcal{S}_{-}^{2}-\frac{10}{\bar{\gamma}}\frac{\Delta
m}{m}\mathcal{S}_{-}\left(\mathcal{S}_{+}\beta_{+}+\mathcal{S}_{-}\beta_{-}\right)\right.\\\
&\left.+\frac{10}{\bar{\gamma}}\mathcal{S}_{-}\left(\mathcal{S}_{-}\beta_{+}+\mathcal{S}_{+}\beta_{-}\right)\right]+\dot{r}^{2}\left[\frac{23}{2}\mathcal{S}_{+}^{2}-8\frac{\Delta
m}{m}\mathcal{S}_{+}\mathcal{S}_{-}+\left(9\eta-\frac{37}{2}-10\bar{\gamma}\right)\mathcal{S}_{-}^{2}\right.-\frac{80}{\bar{\gamma}}\mathcal{S}_{+}\left(\mathcal{S}_{+}\beta_{+}\right.\\\
&\left.+\mathcal{S}_{-}\beta_{-}\right)+\frac{30}{\bar{\gamma}}\frac{\Delta
m}{m}\mathcal{S}_{-}\left(\mathcal{S}_{+}\beta_{+}+\mathcal{S}_{-}\beta_{-}\right)\left.\left.-\frac{10}{\bar{\gamma}}\mathcal{S}_{-}\left(\mathcal{S}_{-}\beta_{+}+\mathcal{S}_{+}\beta_{-}\right)+\frac{120}{\bar{\gamma}^{2}}\left(\mathcal{S}_{+}\beta_{+}+\mathcal{S}_{-}\beta_{-}\right)^{2}\right]\right)\\\
&-\frac{\Delta m}{m}\frac{\eta}{6\,c^{2}}\left(\frac{\alpha
f^{\prime}(\phi_{0})\mathcal{S}_{-}\mathcal{S}_{+}}{\sqrt{\bar{\alpha}}r^{2}}\right)\left(\mathcal{S}_{+}+\frac{\Delta
m}{m}\mathcal{S}_{-}\right)\left[-9\dot{r}^{2}+3v^{2}-\frac{2G\bar{\alpha}m}{r}\right]+\mathcal{O}(c^{-3})\Bigg{\\}}\,.\end{split}$
(61)

In the above expression, there is no contribution to the flux at
$-0.5\mathrm{PN}$ order because the product of the $-0.5\mathrm{PN}$ and
$0\mathrm{PN}$ pieces of $\Phi$ has an odd number of $\hat{n}^{i}$. At 0 PN
order, we have $(-0.5\mathrm{PN})-(+0.5\mathrm{PN})$, and
$0\mathrm{PN}-0\mathrm{PN}$ contributions.

##  7 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 Sec. 7.4, we assume the small coupling limit and
compare our results for scalar waveforms against the NR results of Ref. (Witek
_et al._ , 2019).

###  7.1 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 Sec. 5, does not explicitly appear in the equations of motion
to 1PN order.

Using the relative acceleration (see Eq. (49)) to solve for the circular
condition $\ddot{r}=\dot{r}=0$, we derive the angular velocity $\omega$ in
terms of the relative separation to be

$\omega^{2}=\frac{G\bar{\alpha}m}{r^{3}}\Bigg{\\{}1-\frac{G\bar{\alpha}m}{rc^{2}}\Bigg{[}3-\eta+\bar{\gamma}+2\beta_{+}-2\frac{\Delta
m}{m}\beta_{-}+\frac{4\alpha
f^{\prime}(\phi_{0})}{\bar{\alpha}^{3/2}r^{2}}\left(3\mathcal{S}_{+}+\frac{\Delta
m}{m}\mathcal{S}_{-}\right)\Bigg{]}+\mathcal{O}(c^{-4})\Bigg{\\}}\,.$ (62)

It is useful at this stage to distinguish between the two commonly used PN
parameters, namely

$\begin{split}\gamma_{PN}=\frac{G\bar{\alpha}m}{c^{2}r}\,,\hskip 25.0pt\quad
x=\left(\frac{G\bar{\alpha}m\omega}{c^{3}}\right)^{2/3}\,,\end{split}$ (63)

which differ from their GR definition by an additional factor of
$\bar{\alpha}$. At leading order, one has
$r^{3}\omega^{2}=Gm\bar{\alpha}+\mathcal{O}(c^{-2})$. From Eq. (62), we find
the relation between PN parameters to next-to-leading order,

$x=\gamma_{PN}\Bigg{\\{}1-\frac{\gamma_{PN}}{3}\Bigg{[}3-\eta+\bar{\gamma}+2\beta_{+}-2\frac{\Delta
m}{m}\beta_{-}+\gamma_{PN}^{2}\frac{4c^{4}}{G^{2}}\frac{\alpha
f^{\prime}(\phi_{0})}{\bar{\alpha}^{7/2}m^{2}}\left(3\mathcal{S}_{+}+\frac{\Delta
m}{m}\mathcal{S}_{-}\right)\Bigg{]}+\mathcal{O}(c^{-4})\Bigg{\\}}\,,\hskip
25.0pt$ (64)

where we have substituted
$r=Gm\bar{\alpha}/(c^{2}\gamma_{PN})+\mathcal{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 Eq. (64) we can express the total binding energy (see Eq. (51)) of
circular orbits in terms of the orbital frequency,

$\displaystyle E$ $\displaystyle=-\frac{\mu
c^{2}x}{2}\Bigg{\\{}1+x\left[-\frac{3}{4}-\frac{\eta}{12}-\frac{2\bar{\gamma}}{3}+\frac{2\beta_{+}}{3}-\frac{\Delta
m}{m}\frac{2\beta_{-}}{3}\right.$ (65) $\displaystyle\left.\hskip
25.0pt+\frac{22c^{4}}{3G^{2}}\frac{\alpha
f^{\prime}(\phi_{0})}{\bar{\alpha}^{7/2}m^{2}}x^{2}\left(3\mathcal{S}_{+}+\frac{\Delta
m}{m}\mathcal{S}_{-}\right)\right]+\mathcal{O}(c^{-4})\Bigg{\\}}.\quad$

###  7.2 Gravitational wave phase evolution

In order to compute the full time-dependent waveforms, we need the evolution
of the orbital phase angle $\varphi$, obtained from the angular velocity. In
the adiabatic approximation, i.e. $\dot{\omega}/\omega^{2}\ll 1$, the
luminosity of GWs is equal to the change in orbital energy averaged over a
period, leading to the energy balance equation

$\frac{dE(x)}{dt}=-\mathcal{F}(x)\,,$ (66)

with $\mathcal{F}(x)$ being the total energy flux rate. Using
$\dot{\varphi}=\omega$, the balance equation can be reformulated to

$\frac{d\varphi}{dt}-\frac{c^{3}x^{3}}{G\bar{\alpha}m}=0\,,\hskip
25.0pt\frac{dx}{dt}+\frac{\mathcal{F}(x)}{E^{\prime}(x)}=0\,,$ (67)

where $E^{\prime}(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
$\mathcal{F}/E^{\prime}$. 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 Ref. Buonanno
_et al._ (2009).

From Eq. (65), we derive $E^{\prime}(x)$ to relative 1PN to be

$\begin{split}&E^{\prime}(x)=-\frac{\mu
c^{2}}{2}\left[1+E^{\prime}_{2}x+\mathcal{O}(c^{-4})\right]\,,\\\
&E^{\prime}_{2}=-\frac{3}{2}-\frac{\eta}{6}-\frac{4\bar{\gamma}}{3}+\frac{4}{3}\left(\beta_{+}-\frac{\Delta
m}{m}\beta_{-}\right)+\frac{88c^{4}}{3G^{2}}\frac{\alpha
f^{\prime}(\phi_{0})}{m^{2}\bar{\alpha}^{7/2}}\,x^{2}\left(3\mathcal{S}_{+}+\frac{\Delta
m}{m}\mathcal{S}_{-}\right)\,.\end{split}$ (68)

The total energy flux, constituting Eqns. (59) and (61), has the overall
structure

$\mathcal{F}=\mathcal{F}_{-1,S}+\mathcal{F}_{0,S}+\mathcal{F}_{0,T}+\mathcal{F}_{1,T}\,,$
(69)

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 $\mathcal{F}(x)/E^{\prime}(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 $\mathcal{F}_{-1,S}$ as the
dominant term, and the tensor quadrupolar driven (QD) regime where, instead,
$\mathcal{F}_{0,T}$ dominates. The DD regime is relevant for low frequencies

$x_{\mathrm{D}D}\ll\frac{5}{24}{\mathcal{}S}_{-}^{2}\bar{\alpha},\quad{\mathrm{o}r}\quad
f_{\mathrm{G}W}^{\mathrm{D}D}\ll\left(\frac{5}{24}\right)^{3/2}\frac{c^{3}\,{\mathcal{}S}_{-}^{2}\,\sqrt{\bar{\alpha}}}{\pi
Gm}.$ (70)

At much higher frequencies than this condition, the system is in the QD
regime. Below, we find the 1PN expansion of $\mathcal{F}(x)/E^{\prime}(x)$ in
the two regimes.

####  7.2.1 Dipolar driven regime

For systems with large scalar dipole or very large separation, factoring out
the dipolar scalar flux in Eq. (69) gives,

$\begin{split}\mathcal{F}^{DD}(x)=\frac{4\eta^{2}\mathcal{S}^{2}_{-}c^{5}x^{4}}{3\bar{\alpha}G}&\left\\{1+f^{DD}_{2}x+f^{DD}_{3}x^{3/2}+f^{DD}_{4}x^{2}+\mathcal{O}(c^{-5})\right\\}\,.\end{split}$
(71)

For 1PN corrections to the phase at this regime, it is sufficient to keep the
factor $f_{2}^{DD}=(\mathcal{F}_{0,S}+\mathcal{F}_{0,T})/\mathcal{F}_{-1,S}$,
being

$\begin{split}&f_{2}^{DD}=\frac{24}{5\bar{\alpha}\mathcal{S}_{-}^{2}}+\frac{4}{5}\left(\frac{\mathcal{S}_{+}}{\mathcal{S}_{-}}\right)^{2}-\frac{20}{3}\left(\beta_{+}-\frac{\Delta
m}{m}\beta_{-}\right)-\frac{54}{5}+\frac{4}{3}\eta-\frac{10}{3}\bar{\gamma}+\frac{4}{\bar{\gamma}}\frac{\mathcal{S}_{+}}{\mathcal{S}_{-}}\left(\beta_{-}-\frac{\Delta
m}{m}\beta_{+}\right)\\\ &+\frac{4}{\bar{\gamma}}\left(\beta_{+}-\frac{\Delta
m}{m}\beta_{-}\right)-x^{2}\frac{c^{4}}{G^{2}}\frac{2\alpha
f^{\prime}(\phi_{0})}{m^{2}\bar{\alpha}^{5/2}}\left[\frac{13}{3\bar{\alpha}}\left(3\mathcal{S}_{+}+\frac{\Delta
m}{m}\mathcal{S}_{-}\right)+\eta\frac{\Delta
m}{16m}\frac{\mathcal{S}_{+}}{\mathcal{S}_{-}}\left(\mathcal{S}_{+}+\frac{\Delta
m}{m}\mathcal{S}_{-}\right)\right]\,.\end{split}$ (72)

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 $\mathcal{F}^{DD}(x)/E^{\prime}(x)$ to 1PN order we find

$\frac{\mathcal{F}^{DD}(x)}{E^{\prime}(x)}=-\frac{8\eta
c^{3}\mathcal{S}_{-}^{2}x^{4}}{3G\bar{\alpha}m}\left[1+(f_{2}^{DD}-E^{\prime}_{2})x+\mathcal{O}(c^{-4})\right]\,.$
(73)

Note that this expression is complete at the relative 1PN order.

####  7.2.2 Quadrupolar driven regime

For quadrupolar driven systems, the flux is expanded about the Newtonian-order
term $\mathcal{F}_{0,T}$ such that

$\mathcal{F}^{QD}(x)=\frac{32\eta^{2}c^{5}x^{5}}{5G\bar{\alpha}^{2}}\left[\xi+f_{2,T}^{nd}x+\mathcal{O}(c^{-3})\right]+x^{4}\frac{4\eta^{2}\mathcal{S}_{-}^{2}}{3G\bar{\alpha}}\,.$
(74)

The factor $\xi=5\bar{\alpha}f_{2}^{DD}\mathcal{S}_{-}^{2}/24$, captures the
0PN flux terms. For the coefficient
$f_{2}^{nd}=\mathcal{F}_{1,T}/\mathcal{F}_{0,T}\equiv
f_{2,T}^{nd}+f_{2,GB}^{nd}$ we have

$\displaystyle f_{2,T}^{nd}=$
$\displaystyle\frac{-1}{336}\left[1247+980\eta+448\bar{\gamma}+896\left(\beta_{+}-\frac{\Delta
m}{m}\beta_{-}\right)\right]$ $\displaystyle f_{2,GB}^{nd}=$
$\displaystyle-\frac{c^{4}}{G^{2}}\frac{f^{\prime}(\phi_{0})\alpha\mathcal{S}_{+}}{m^{2}\bar{\alpha}^{5/2}}x^{2}\Bigg{[}12+\frac{495}{28}(1-2\eta)+\frac{32}{3}\frac{\mathcal{S}_{-}}{\mathcal{S}_{+}}\frac{\Delta
m}{m}+\frac{495}{28}\frac{\mathcal{S}_{-}}{\mathcal{S}_{+}}(1+2\eta)\Bigg{]}\,,$
(75)

As mentioned above, these expansions miss a contribution from
$\mathcal{F}_{1,S}$ which we are unable to compute within our approximations,
having $\dot{E}_{S}$ to 0PN order. The inclusion of the missing term would, in
principle, increase the overall energy flux, and thus the difference between
sGB and GR waveforms. In the following, we set these contributions to zero,
similar to what was done in a similar situation for tidal effects in GR Damour
_et al._ (2012).
Overall, the expansion of $\mathcal{F}^{QD}(x)/E^{\prime}(x)$ to 1PN order
becomes

$\displaystyle\frac{\mathcal{F}^{QD}(x)}{E^{\prime}(x)}=$
$\displaystyle-\frac{64\eta
c^{3}x^{5}}{5Gm\bar{\alpha}^{2}}\left[\left(\xi-E^{\prime}_{2}\frac{5\bar{\alpha}\mathcal{S}_{-}^{2}}{24}\right)+(f_{2,T}^{nd}-\xi
E^{\prime}_{2})x+\frac{5\bar{\alpha}\mathcal{S}_{-}^{2}}{24}x^{-1}+\mathcal{O}(c^{-4})\right]\,.$
(76)

Having Eqns. (73) and (76), together with an adequate choice of parameters, we
numerically solve Eq. (67) to determine the time evolution of $\varphi$ and
$\omega$ in different regimes, and therefore, time-domain waveforms.

###  7.3 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
$\\{\hat{\mathbf{n}},\hat{\mathbf{p}},\hat{\mathbf{q}}\\}$, with
$\hat{\mathbf{n}}$ being the radial direction to the observer,
$\hat{\mathbf{p}}$ lying along the intersection of the orbital plane with the
plane of the sky, and
$\hat{\mathbf{q}}=\hat{\mathbf{n}}\times\hat{\mathbf{p}}$. A schematic view of
this choice can be found in Fig. 2.

Figure 2: Orientation of the orbital plane and the sky plane based on the
orthonormal triad $\\{\hat{\mathbf{n}},\hat{\mathbf{p}},\hat{\mathbf{q}}\\}$.
$\mathbf{\hat{p}}$ lie along the lines of nodes and defines an origin for the
orbital phase angle $\varphi$. $\hat{\lambda}$ is the unit vector of the
orbital angular momentum.

The normal to the orbit is inclined at an angle $i$ relative to
$\hat{\mathbf{n}}$. The orbital phase angle of body $A$ is measured from the
line of nodes in a positive sense. We have that

$\begin{split}&\hat{\mathbf{n}}_{AB}=\hat{\mathbf{p}}\cos\varphi+(\hat{\mathbf{q}}\cos
i+\hat{\mathbf{n}}\sin i)\sin\varphi\,,\\\
&\hat{\bm{\lambda}}=-\hat{\mathbf{p}}\sin\varphi+\left(\hat{\mathbf{q}}\cos
i+\hat{\mathbf{n}}\sin i\right)\cos\varphi\,,\end{split}$ (77)

where $\bm{v}=r\omega\hat{\bm{\lambda}}$ for circular orbits.

Gravitational wave detectors measure a linear combination of polarization
waveforms $h_{+}(t)$ and $h_{\times}(t)$ that are defined by the projections

$\begin{split}&h_{+}=\frac{1}{2}\left(\hat{p}_{i}\hat{p}_{j}-\hat{q}_{i}\hat{q}_{j}\right)h^{ij}\,,\hskip
25.0pth_{\times}=\frac{1}{2}\left(\hat{p}_{i}\hat{q}_{j}+\hat{q}_{i}\hat{p}_{j}\right)h^{ij}\,.\end{split}$
(78)

Applying Eqns. (77) and (78) on Eq. (53) and simplifying combinations such as
$\hat{n}^{i}\hat{n}^{j},\hat{\lambda}^{i}\hat{\lambda}^{j},$ and
$\hat{n}^{(i}\hat{\lambda}^{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 Fig. 3 by plotting analytical GWs of BH binary
systems in EdGB theory and comparing them with 1PN GR waveforms. Using Eq.
(7), we write BH scalar charges, and the related parameters, in terms of
$\alpha$. The GB coupling parameter $\alpha$ is thus the only parameter that
we have to specify. Exemplarily, we choose $\alpha=0.01m^{2}$. Also note that,
as shown In Sec. 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_{\odot}$ 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 Fig.2 of Ref. Shiralilou  _et al._ (2020).

Figure 3 shows waveforms in the aforementioned two cases. The observer is
viewing the orbit edge on, so that $i=\pi/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 Fig. 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 Eq. (31)).
Despite the amplitudes, the dephasing between EdGB and GR waveforms starts
early in the evolution and could thus be a detectable feature.

Figure 3: Time evolution of gravitational waveform and orbital frequency for a
$m=15M_{\odot}$ binary with $q=1/2$ (left) and $q=1/4$ (right), and
$\alpha=0.01m^{2}$. Blue dashed curves correspond to EdGB waveforms and black
curves correspond to 1PN GR waveforms. Orbits are viewed edge-on ($i=\pi/2$),
and $t=0$ corresponds to $f\approx 50Hz$. The green shaded regions (top) show
$1s$ interval of the waveforms at intermediate times, whilst red shaded
regions (middle) show the $0.3s$ intervals where EdGB binaries are close to
the merger.

###  7.4 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 Eq. (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 Fig. 4 we display the scalar waveform emitted by quasi-circular binaries,
corresponding to the gravitational waveforms shown in Fig. 3. In these
examples, for the mass ratio $q=1/2$, the scalar wave amplitude is suppressed
with respect to the gravitational wave 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

$\Phi_{lm}(t,R)=\int
d\Omega\,\Phi(t,R,\hat{\Theta},\hat{\Phi})Y^{*}_{lm}(\hat{\Theta},\hat{\Phi})\,,$
(79)

where, in order to express $\Phi$ in terms of spherical coordinate variables,
we choose to work with $\hat{\bm{p}}=-\mathbf{e}_{\hat{\Phi}}$,
$\hat{\bm{q}}=\mathbf{e}_{\hat{\Theta}}$, such that $\hat{\Theta}=i$ and
$\hat{\Phi}=\pi/2-\omega t$. It is easy to verify that the leading-order
contribution to the radiation comes from $\Phi_{l=m}$ modes, $i.e.$,

$\Phi_{ll}(t,R)=\int d\Omega\,\Phi_{l}Y^{*}_{ll}(\hat{\Theta},\hat{\Phi})\,.$
(80)

We use Eq. (7) to re-write the scalar-charge dependent parameters in terms of
the GB coupling. For ssGB gravity we have

$\displaystyle\alpha_{A}^{0}$
$\displaystyle=-\frac{\alpha}{m_{A}^{2}}+\mathcal{O}(\alpha^{2})\,,\quad\bar{\alpha}=1+\frac{\alpha^{2}}{2m_{A}^{2}m_{B}^{2}}+\mathcal{O}(\alpha^{3})\,,$
$\displaystyle\mathcal{S}_{+}$
$\displaystyle=-\frac{\alpha}{2\sqrt{\bar{\alpha}}}\left(\frac{1}{m_{A}^{2}}+\frac{1}{m_{B}^{2}}\right)+\mathcal{O}(\alpha^{2})\,,$
$\displaystyle\mathcal{S}_{-}$
$\displaystyle=-\frac{\alpha}{2\sqrt{\bar{\alpha}}}\left(\frac{1}{m_{A}^{2}}-\frac{1}{m_{B}^{2}}\right)+\mathcal{O}(\alpha^{2})\,,$
(81)

with the rest of the parameters being identically zero.
Specifying Eq. (54) to quasi-circular orbits with the above-mentioned
parameters, we find the various $\Phi_{l}$ modes, required for Eq. (80), to be

$\begin{split}\Phi_{1}&=\frac{G\alpha}{c^{2}R}\left(\frac{\Delta
m}{m_{A}m_{B}}x_{0}^{1/2}-\frac{m\Delta
m}{2m_{A}m_{B}}x_{0}^{3/2}\right)\,,\\\
\Phi_{2}&=-\frac{G\alpha}{c^{2}R}\frac{(m_{B}^{2}-m_{A}m_{B}+m_{A}^{2})}{mm_{A}m_{B}}x_{0}\,,\\\
\Phi_{3}&=-\frac{G\alpha}{c^{2}R}\frac{\Delta
m(m_{A}^{2}+m_{B}^{2})}{8m^{2}m_{A}m_{B}}x_{0}^{3/2}\,,\end{split}$ (82)

where the parameter $x_{0}=(Gm\omega/c^{3})^{2/3}+\mathcal{O}(\alpha)$ is the
leading order term of the PN parameter $x$.

Figure 4: The time evolution of scalar field waveform for $m=15M_{\odot}$
binary with (left) $q=1/2$ and (right) $q=1/4$, and $\alpha=0.01m^{2}$,
corresponding to Fig. 3 binaries and time intervals. Orbits are viewed edge-on
($i=\pi/2$) and thus the waves deviate from a full sinusoidal.

We compare our PN scalar waveforms against previous PN calculations of Ref.
Yagi  _et al._ (2012) and also the results from a numerical relativity
simulation reported in Ref. Witek  _et al._ (2019), 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 $\alpha$
in ssGB theory, and it covers about ten orbits before the binary BHs merge. As
expected, our results agree with those of Ref. (Witek  _et al._ , 2019). 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
$\Phi_{11}$ at order $\mathcal{O}(\alpha)$. In general, we also derive 0.5PN
corrections to $\Phi_{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 $\mathcal{O}(\alpha)$ corrections.

In order to have a meaningful comparison between the new PN dipolar radiation
terms and the results of Ref. (Witek  _et al._ , 2019), we obtain the orbital
frequency evolution from the derivative of the GW phase of numerical data,
instead of using the analytical results of Sec. 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 Fig. 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 Ref. (Yagi  _et al._ , 2012) (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_{\text{merger}}:\\{-1410,-900\\}$. In this way, we avoid spliting the
waves by introducing arbitrary fudge factors.

Figure 5: Time evolution of the scalar mode $\Phi_{11}$ (left) and its
absolute value (right), re-scaled by the extraction radius $R=100m$. Black
curves indicate results from NR simulations. The blue dashed curve indicates
the analytical PN inspiral results of Ref. Yagi  _et al._ (2012) to -0.5 PN
order. The red dashed curve shows the 0.5PN order results of this paper.

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 Ref. Yagi  _et al._ (2012). 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.

##  8 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., Refs. Abbott  _et
al._ (2019b, 2016b), 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 Sec. 7.2. Using the stationary phase approximation (SPA) and
focusing on waveforms with GW Fourier phase $\psi$ at frequency
$f=\omega/\pi$, the waveform in the frequency domain can be written as

$\begin{split}\tilde{h}^{\mathrm{spa}}(f)&=a\left(t_{f}\right)\sqrt{\frac{2\pi}{m\dot{\omega}\left(t_{f}\right)}}e^{-i\left[\psi_{f}\left(t_{f}\right)+\pi/4\right]}\,,\hskip
25.0pt\psi_{f}(t)\equiv 2\phi(t)-2\pi ft\,,\end{split}$ (83)

where $t_{f}$ is the time when the GW frequency becomes equal to the Fourier
variable $f$, by solving $d\psi_{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

$\begin{split}t(\bar{v})&=t(\bar{v}_{\mathrm{ref}})+\int_{\bar{v}}^{\bar{v}_{\mathrm{ref}}}\frac{E^{\prime}(\bar{v})}{\mathcal{F}(\bar{v})}d\bar{v}\,,\\\
\phi(\bar{v})&=\phi(\bar{v}_{\mathrm{ref}})+\frac{c^{3}}{G\bar{\alpha}m}\int_{\bar{v}}^{\bar{v}_{\mathrm{ref}}}\bar{v}^{3}\frac{E^{\prime}(\bar{v})}{\mathcal{F}(\bar{v})}d\bar{v}\,,\end{split}$
(84)

where $\bar{v}=x^{1/2}=(G\bar{\alpha}m\omega/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

$\psi_{f}(t_{f})=2\left.\left(\varphi(v)-\frac{c^{3}}{G\bar{\alpha}m}\bar{v}^{3}t(\bar{v})\right)\right|_{\bar{v}=\bar{v}_{f}}\,,$
(85)

where $\bar{v}_{f}\equiv(\pi Gm\bar{\alpha}f/c^{3})^{1/3}$. As before, we
split the calculation into DD regime and QD regime, and evaluate Eq. (84)
term-by-term by re-expanding the expressions in the PN variable $\bar{v}$ and
truncating them at 1PN order.

###  8.1 Dipolar-driven regime

In the DD regime, expanding the ratio
$E^{\prime}(\bar{v})/\mathcal{F}^{DD}(\bar{v})$ to 1PN order gives

$\frac{E^{\prime}(\bar{v})}{\mathcal{F}^{DD}(\bar{v})}=-\frac{3G\bar{\alpha}m}{4c^{3}\eta\mathcal{S}_{-}^{2}\bar{v}^{7}}\left[1+\left(E^{\prime}_{0}-f_{2}^{DD}\right)\bar{v}^{2}+\mathcal{O}\left(c^{-4}\right)\right]\,.$
(86)

By substituting Eq. (86) into Eq. (84), and integrating term-by-term we find
the phasing to be

$\begin{split}\psi(t_{f})&=-\frac{1}{4\eta\mathcal{S}_{-}^{2}\bar{v}_{f}^{3}}\Bigg{[}1+\frac{9}{2}\rho^{DD}\bar{v}_{f}^{2}+6\rho_{GB}^{DD}\log(\bar{v}_{f})\bar{v}_{f}^{6}+\left(\frac{\bar{v}_{f}}{\bar{v}_{ref}}\right)^{6}\Big{(}1+\frac{3}{2}\bar{v}_{ref}^{2}-6\rho_{GB}^{DD}\bar{v}_{ref}^{6}\log(\bar{v}_{ref})\Big{)}\\\
&-2\left(\frac{\bar{v}_{f}}{\bar{v}_{ref}}\right)^{3}\Big{(}1+3\rho^{DD}\bar{v}_{ref}^{2}-\rho^{DD}_{GB}\bar{v}_{ref}^{6}\Big{)}\Bigg{]}+\phi(\bar{v}_{ref})-2\pi
ft(\bar{v}_{ref})\,,\end{split}$ (87)

where one can easily verify that

$\displaystyle\rho^{DD}$
$\displaystyle=\frac{93}{10}-\frac{3}{2}\eta+2\bar{\gamma}+\frac{24}{3}\left(\beta_{+}-\frac{\Delta
m}{m}\beta_{-}\right)-\frac{4}{5}\left(\frac{\mathcal{S}_{+}}{\mathcal{S}_{-}}\right)^{2}$
$\displaystyle-\frac{24}{5\bar{\alpha}\mathcal{S}_{-}^{2}}-\frac{4}{\bar{\gamma}}\frac{\mathcal{S}_{+}}{\mathcal{S}_{-}}\left(\beta_{-}-\frac{\Delta
m}{m}\beta_{+}\right)-\frac{4}{\bar{\gamma}}\left(\beta_{+}-\frac{\Delta
m}{m}\beta_{-}\right)\,,$ $\displaystyle\rho^{DD}_{GB}$
$\displaystyle=\frac{c^{4}}{G^{2}}\frac{\alpha
f^{\prime}(\phi_{0})}{m^{2}\bar{\alpha}^{5/2}}\left[\frac{104}{3\bar{\alpha}}\left(3\mathcal{S}_{+}+\frac{\Delta
m}{m}\mathcal{S}_{-}\right)+\eta\frac{\Delta
m}{16m}\frac{\mathcal{S}_{+}}{\mathcal{S}_{-}}\left(\mathcal{S}_{+}+\frac{\Delta
m}{m}\mathcal{S}_{-}\right)\right]\,.$ (88)

###  8.2 Quadrupolar-driven regime

For the QD regime, by definition (see Eq. (70)), we expect the dipole
radiation to be negligible. This means that we can use the parameter
$\mathcal{S}_{-}$ as an identifier of small terms, and thus, split the flux
into two pieces as follows,

$\begin{split}&\mathcal{F}^{\mathrm{QD}}=\mathcal{F}_{\mathrm{non}-\mathrm{dip}}+\mathcal{F}_{\mathrm{dip}}\,,\\\
\\\ &\mathcal{F}_{\text{non-dip }}\equiv\lim_{\mathcal{S}_{-}\rightarrow
0}\mathcal{F}\,,\hskip 25.0pt\mathcal{F}_{\text{dip
}}\equiv\mathcal{F}-\mathcal{F}_{\text{non-dip}}\,.\end{split}$ (89)

where $\mathcal{F}_{\mathrm{d}ip}$ indicates the tensor and scalar flux terms
that depend on $\mathcal{S}_{-}$ and $\mathcal{F}_{\mathrm{non-dip}}$ denotes
the terms that do not depend on $\mathcal{S}_{-}$.

With this division, to first order in the small quantity
$\mathcal{F}_{\text{dip }}/\mathcal{F}_{\text{non-dip}}$, the ratio
$E^{\prime}(\bar{v})/\mathcal{F}(\bar{v})$ can be approximated as

$\frac{E^{\prime}(\bar{v})}{\mathcal{F}(\bar{v})}\simeq\frac{E^{\prime}(\bar{v})}{\mathcal{F}_{\text{non-
dip }}(\bar{v})}\left(1-\frac{\mathcal{F}_{\text{dip
}}(\bar{v})}{\mathcal{F}_{\text{non-dip }}(\bar{v})}\right)\,.$ (90)

We obtain the dipolar and non-dipolar parts to be

$\displaystyle\mathcal{F}_{\text{non-dip
}}(\bar{v})=\frac{32\eta^{2}\bar{\xi}c^{5}}{5G\bar{\alpha}^{2}}\bar{v}^{10}\left[1+f_{2}^{nd}\bar{v}^{2}+\mathcal{O}\left(c^{-3}\right)\right],$
$\displaystyle\mathcal{F}_{\text{dip
}}(\bar{v})=\frac{4\mathcal{S}_{-}^{2}\eta^{2}c^{5}}{3G\bar{\alpha}}\bar{v}^{8}\left[1+f_{2}^{d}\bar{v}^{2}+\mathcal{O}\left(c^{-3}\right)\right]\,,$
(91)

where $\bar{\xi}=(1+\mathcal{S}_{+}^{2}\bar{\alpha}/6)$ comes from the non-
dipolar flux terms at Newtonian order. Note that this factor differs from the
previously defined factor $\xi$ by the fact that it only includes the flux
terms that do not depend on $\mathcal{S}_{-}$ explicitly. Similarly, the
factor $f^{nd}_{2}$ equals those terms in Eq. (7.2.2) that do not depend on
$\mathcal{S}_{-}$, replacing $\xi$ by $\bar{\xi}$. 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/\mathcal{S}_{-}^{2}$. Overall, for Eq. (90) we find

$\begin{split}&\frac{E^{\prime}(\bar{v})}{\mathcal{F}(\bar{v})}\simeq\frac{5Gm\bar{\alpha}^{2}}{32c^{3}\eta\bar{\xi}\bar{v}^{9}}\left[1+\left(E^{\prime}_{0}-f^{nd}_{2}\right)\bar{v}^{2}+\mathcal{O}\left(c^{-4}\right)\right]\\\
&\hskip
25.0pt-\frac{25Gm\bar{\alpha}^{3}\mathcal{S}_{-}^{2}}{768c^{3}\bar{\xi}^{2}\eta\bar{v}^{11}}\left[1+\left(E^{\prime}_{0}-2f^{nd}_{2}+f^{d}\right)\bar{v}^{2}+\mathcal{O}\left(c^{-4}\right)\right]\,.\end{split}$
(92)

Integrating Eq. (84) by using the above expression we obtain

$\displaystyle\Psi(t_{f})$ $\displaystyle=\psi_{\text{non-
dip}}+\psi_{\text{dip}}+\phi(\bar{v}_{ref})-2\pi ft(\bar{v}_{ref})\,\hskip
25.0pt$ (93) $\displaystyle\psi_{\text{non-
dip}}(t_{f})=-\frac{6\bar{\alpha}}{256\eta\bar{\xi}\bar{v}_{f}^{5}}\Bigg{[}1+\frac{20}{9}\rho^{nd}\bar{v}_{f}^{2}-20\rho^{nd}_{GB}\bar{v}_{f}^{6}$
$\displaystyle+\frac{5}{3}\left(\frac{\bar{v}_{f}}{\bar{v}_{ref}}\right)^{8}\left(1+\frac{4}{3}\rho^{nd}\bar{v}_{ref}^{2}+4\rho^{nd}_{GB}\bar{v}_{ref}^{6}\right)-\frac{8}{3}\left(\frac{\bar{v}_{f}}{\bar{v}_{ref}}\right)^{5}\left(1+\frac{5}{3}\rho^{nd}\bar{v}_{ref}^{2}-5\rho^{nd}_{GB}\bar{v}_{ref}^{6}\right)\Bigg{]}\,,\hskip
25.0pt$ (94)
$\displaystyle\psi_{\text{dip}}(t_{f})=\frac{10\mathcal{S}_{-}^{2}\bar{\alpha}^{2}}{3584\eta\bar{\xi}^{2}\bar{v}_{f}^{7}}\left[1+\frac{7}{4}\rho^{d}\bar{v}_{f}^{2}+\frac{70}{4}\rho^{d}_{GB}\bar{v}_{f}^{6}\right.$
$\displaystyle+\left.\frac{7}{3}\left(\frac{\bar{v}_{f}}{\bar{v}_{ref}}\right)^{10}\left(1+\frac{5}{4}\rho^{d}\bar{v}_{f}^{2}+\frac{5}{2}\rho_{GB}^{d}\bar{v}_{f}^{6}\right)\frac{10}{3}\left(\frac{\bar{v}_{f}}{\bar{v}_{ref}}\right)^{7}\left(1+\frac{7}{3}\rho^{d}\bar{v}_{f}^{2}+7\rho_{GB}^{d}\bar{v}_{f}^{6}\right)\right]\,,\hskip
25.0pt\quad$ (95)

with the coefficients being

$\begin{split}&\rho^{nd}=\frac{1247}{336\bar{\xi}}-\frac{3}{2}+\left(\frac{980}{336\bar{\xi}}-\frac{1}{6}\right)\eta+\left(\frac{448}{336\bar{\xi}}-\frac{4}{3}\right)\bar{\gamma}+\left(\frac{896}{336\bar{\xi}}+\frac{4}{3}\right)\left(\beta_{+}-\frac{\Delta
m}{m}\beta_{-}\right)\,,\\\
&\rho^{nd}_{GB}=\frac{c^{4}}{G^{2}}\frac{f^{\prime}(\phi_{0})\alpha}{m^{2}\bar{\alpha}^{5/2}}\left[\frac{12\mathcal{S}_{+}}{\bar{\xi}}+\frac{495\mathcal{S}_{+}}{28\bar{\xi}}(1-2\eta)+\frac{88}{3\bar{\alpha}}\left(3\mathcal{S}_{+}+\frac{\Delta
m}{m}\mathcal{S}_{-}\right)\right]\,,\\\
&\rho^{d}=\frac{1247}{168\bar{\xi}}-\frac{123}{10}+\left(\frac{980}{168\bar{\xi}}+\frac{7}{6}\right)\eta+\left(\frac{448}{168\bar{\xi}}-\frac{14}{3}\right)\bar{\gamma}+\left(\frac{896}{168\bar{\xi}}-\frac{32}{3}\right)\left(\beta_{+}-\frac{\Delta
m}{m}\beta_{-}\right)\\\ &\hskip
25.0pt-\frac{4}{\bar{\gamma}}\left(\frac{\mathcal{S}_{+}}{\mathcal{S}_{-}}\right)\left(\beta_{-}-\frac{\Delta
m}{m}\beta_{+}\right)-\frac{4}{\bar{\gamma}}\left(\beta_{+}-\frac{\Delta
m}{m}\beta_{-}\right)\,,\\\
&\rho^{d}_{GB}=\frac{c^{4}}{G^{2}}\frac{f^{\prime}(\phi_{0})\alpha}{m^{2}\bar{\alpha}^{5/2}}\left[\frac{24\mathcal{S}_{+}}{\bar{\xi}}+\frac{495\mathcal{S}_{+}}{14\bar{\xi}}(1-2\eta)\right.\\\
&\hskip
25.0pt\left.+\frac{104}{3\bar{\alpha}}\left(3\mathcal{S}_{+}+\frac{\Delta
m}{m}\mathcal{S}_{-}\right)+\frac{\Delta
m}{16m}\frac{\mathcal{S}_{+}}{\mathcal{S}_{-}}\eta\left(\mathcal{S}_{+}+\frac{\Delta
m}{m}\mathcal{S}_{-}\right)\right]\,.\end{split}$ (96)

Note, that the current mapping between ppE templates and waveforms in sGB
theory is known for the waveform at Newtonian order (corresponding to
$b_{ppE}=-7$) Yunes  _et al._ (2016); Yagi  _et al._ (2012). The results
presented 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 Nair  _et al._ (2019); Perkins  _et al._
(2021).

In Fig. 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 $\alpha=0.01m^{2}$ and $m=15M_{\odot}$, as it
corresponds to the systems considered in Fig. 3. The upper bound on frequency
is chosen as $f_{\text{max}}=2(6^{3/2}\pi m)^{-1}\approx 586\,\mathrm{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\neq 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 $\mathcal{O}(1)$ cycles and thus is more
difficult to distinguish from the deviations from GR.

Figure 6: Left: The GW phase difference between EdGB $\psi_{\text{EdGB}}$ and
1PN GR $\psi_{\text{GR}}$ shown for the systems with $\alpha=0.01m^{2}$ and
$m=15M_{\odot}$, as considered in Fig. 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.

###  8.3 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_{\mu\nu}$, 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 Eq. (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
$\mathcal{A}(\phi)=1+\mathcal{O}(\phi)$, we urge the reader (and tester of
strong-field gravity) to caution outside this approximation.

##  9 Conclusions

We have extensively studied the generation of GWs and scalar radiation in
inspiralling black-hole 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 post-
Newtonian expansion, we have computed the GWs and energy flux to relative 1PN
order beyond the quadrupole emission, and the scalar analogue to relative 0.5
PN order. We advanced the knowledge of GW and scalar waveforms in sGB theories
beyond the leading results of Ref. Yagi  _et al._ (2012) 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 post-Newtonian scheme of Ref.
Yagi  _et al._ (2012) and obtained with numerical relativity simulations Witek
_et al._ (2019).

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 Ref. (Julié and Berti, 2019). 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 Eq. (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, Eqs. (126) to (128), or the
GW phasing in the Fourier domain, Eqs. (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 Shiralilou  _et al._ (2020), 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 Aasi  _et al._ (2015); Acernese  _et al._ (2015); Akutsu  _et al._
(2020), 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, ground- and space-based GW
detector networks combining high-frequency observations by, e.g., the Einstein
Telescope Punturo  _et al._ (2010) or Cosmic Explorer Abbott  _et al._ (2017),
with low-frequency observations by, e.g., the space-based LISA mission Amaro-
Seoane  _et al._ (2017), 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 Silva  _et al._ (2020) 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.

##  10 Acknowledgment

We thank L. Gualtieri and H. O. Silva for useful discussions, and A. Samajdar
for comments on the manuscript. BS, TH and SMN are grateful for financial
support from the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO)
through the Projectruimte and VIDI grants (Nissanke). TH and SMN also
acknowledges financial support from the NWO sectorplan. NO acknowledges
financial support by the CONACyT grants “Ciencia de frontera” 140630 and
376127, and by the UNAM-PAPIIT grant IA100721. HW acknowledges financial
support provided by NSF Grant No. OAC-2004879 and the Royal Society Research
Grant No. RGF\R1\180073. TH also acknowledges Cost Actions CA18108 - Quantum
gravity phenomenology in the multi-messenger approach and CA16104 -
Gravitational waves, black holes and fundamental physics. Some calculations
used the computer algebra system Mathematica, in combination with the
xAct/xTensor suite Brizuela  _et al._ (2009).

##  Appendix A Curvature Tensors in terms of gothic metric

###  A.1 Transformation rules

In order to re-write curvature tensors in terms of the gothic metric, defined
by $\mathfrak{g}^{ab}=\sqrt{-g}g^{ab}$, we first define its inverse
$\mathfrak{g}_{ab}=g_{ab}/\sqrt{-g}$, which is a tensor density of weight -1.
To find Christoffel symbols, we first need to convert the derivative of the
metric to the derivative of the gothic one. This is done by differentiating
the relation $\mathfrak{g}^{ab}g_{ac}=\sqrt{-g}\delta_{c}^{b}$. In the case of
$b=c$, we have

$\begin{split}0&=\partial_{c}\left(\mathfrak{g}^{ab}g_{ab}\right)-4\partial_{c}\sqrt{-g}=g_{ab}\partial_{c}\mathfrak{g}^{ab}+\mathfrak{g}^{ab}g_{ab,c}-2\sqrt{-g}g^{ab}g_{ab,c}=g_{ab}\partial_{c}\mathfrak{g}^{ab}-\mathfrak{g}^{ab}\partial_{c}g_{ab}\,,\end{split}$
(97)

and in the case of $b\neq c$,

$\begin{split}0&=g_{bj}\left[\partial_{c}\left(\mathfrak{g}^{ab}g_{af}\right)-\delta_{f}^{b}\partial_{c}\sqrt{-g}\right]=\sqrt{-g}g_{jf,c}+(-g)\mathfrak{g}_{af}\mathfrak{g}_{bj}\partial_{c}\mathfrak{g}^{ab}-\frac{(-g)}{2}\mathfrak{g}_{fj}\mathfrak{g}_{ab}\partial_{c}\mathfrak{g}^{ab}\\\
&=g_{jf,c}+\sqrt{-g}\left(\mathfrak{g}_{af}\mathfrak{g}_{bj}-\frac{1}{2}\mathfrak{g}_{ab}\mathfrak{g}_{fj}\right)\partial_{c}\mathfrak{g}^{ab}\,.\end{split}$
(98)

Equation (98) can be re-written to relate the partial derivative of the gothic
metric to that of its inverse,

$\partial_{c}\mathfrak{g}_{jf}=-\mathfrak{g}_{af}\mathfrak{g}_{bj}\partial_{c}\mathfrak{g}^{ab}.$
(99)

###  A.2 Christoffel symbols and curvature tensors

By using Eqns. (97)-(99), we can express quantities in terms of gothic
variables. In particular, the Christoffel symbols become

$\Gamma^{a}_{ij}=\frac{1}{4}\left(2\mathfrak{g}^{ah}\mathfrak{g}_{dj}\mathfrak{g}_{gi}\partial_{h}\mathfrak{g}^{dg}-\mathfrak{g}^{ah}\mathfrak{g}_{dg}\mathfrak{g}_{ji}\partial_{h}\mathfrak{g}^{dg}-2\mathfrak{g}_{ej}\partial_{i}\mathfrak{g}^{ae}+\delta_{j}^{a}\mathfrak{g}_{be}\partial_{i}\mathfrak{g}^{be}-2\mathfrak{g}_{fi}\partial_{j}\mathfrak{g}^{af}+\delta_{i}^{a}\mathfrak{g}_{cf}\partial_{j}\mathfrak{g}^{cf}\right).$
(100)

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

$\begin{split}R=-\frac{1}{\sqrt{-\mathfrak{g}}}\Big{(}&-\frac{1}{4}\mathfrak{g}^{{ab}}\mathfrak{g}_{{cd}}\mathfrak{g}_{{ef}}\partial_{{a}}\mathfrak{g}^{ce}\partial_{{b}}\mathfrak{g}^{{df}}-\frac{1}{8}\mathfrak{g}^{{ab}}\mathfrak{g}_{{cd}}\mathfrak{g}_{{ef}}\partial_{{a}}\mathfrak{g}^{{cd}}\partial_{{b}}\mathfrak{g}^{{ef}}+\partial_{{b}}\partial_{{a}}\mathfrak{g}^{{ab}}\\\
&+\frac{1}{2}\mathfrak{g}^{{ab}}\mathfrak{g}_{{cd}}\partial_{b}\partial_{{a}}\mathfrak{g}^{{cd}}-\frac{1}{2}\mathfrak{g}_{{ab}}\partial_{{c}}\mathfrak{g}^{{bd}}\partial_{{d}}\mathfrak{g}^{{ac}}+\frac{1}{2}\mathfrak{g}_{{ab}}\partial_{c}\mathfrak{g}^{{ab}}\partial_{{d}}\mathfrak{g}^{{cd}}\Big{)}\,.\end{split}$
(101)

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*}\propto\frac{1}{\sqrt{-\mathfrak{{g}}^{3}}}\Big{(}\mathfrak{g}^{..}\mathfrak{g}^{..}\partial_{.}\partial_{.}\mathfrak{g}^{..},\,\mathfrak{g}^{..}\partial_{.}\mathfrak{g}^{..}\partial_{.}\mathfrak{g}^{..},\,\mathfrak{g}^{..}\mathfrak{g}^{..}\mathfrak{g}_{..}\partial_{.}\mathfrak{g}^{..}\partial_{.}\mathfrak{g}^{..},\,\mathfrak{g}^{..}\mathfrak{g}^{..}\mathfrak{g}^{..}\mathfrak{g}_{..}\partial_{.}\partial_{.}\mathfrak{g}^{..},\,\mathfrak{g}^{..}\mathfrak{g}^{..}\mathfrak{g}^{..}\mathfrak{g}_{..}\mathfrak{g}_{..}\partial_{.}\mathfrak{g}^{..}\partial_{.}\mathfrak{g}^{..}\Big{)}\,.$
(102)

In the main text, we called this quantity ${}^{*}\hat{R}^{abcd*}$. The
expression for the $\mathcal{R}_{GB}$ is derived from the above mentioned
curvature quantities through Eq. (2).

###  A.3 Weak-field limit

The gauge-fixed GB coupling-dependent quantities $(^{*}R^{*\alpha a\beta
d})\nabla_{ad}f(\phi)$ and $\mathcal{R}^{2}_{GB}$ enter as part of the source
of the wave equations (11) and (12), respectively. Using Eq. (10), We expand
these quantities to second order in the weak-field limit as higher order terms
are not relevant for 1PN calculations. It is easy to verify that the expansion
of these quantities to first order in the weak-field limit is identically
zero. The second order expansions are

$\begin{split}{}^{*}&\hat{R}^{*\alpha a\beta
d}\nabla_{ad}f(\phi)=f^{\prime}(\phi_{0})\nabla_{ad}\delta\phi\Big{(}-4\eta^{\alpha
e}\eta^{df}\partial_{fe}h^{a\beta}+4\eta^{\alpha
d}\eta^{ef}\partial_{fe}h^{a\beta}+4\eta^{\alpha e}\eta^{\beta
f}\partial_{fe}h^{ad}\\\
&-4\eta^{\alpha\beta}\eta^{ef}\partial_{f}\partial_{e}h^{ad}+4\eta^{ae}\eta^{df}\partial_{fe}h^{\alpha\beta}-4\eta^{ad}\eta^{ef}\partial_{fe}h^{\alpha\beta}-4\eta^{ae}\eta^{\beta
f}\partial_{fe}h^{\alpha d}+4\eta^{a\beta}\eta^{ef}\partial_{fe}h^{\alpha
d}\\\ &+2\eta^{ae}\eta^{\alpha d}\eta^{\beta
f}\eta_{gh}\partial_{fe}h^{gh}-2\eta^{ad}\eta^{\alpha e}\eta^{\beta
f}\eta_{gh}\partial_{fe}h^{gh}-2\eta^{ae}\eta^{\alpha\beta}\eta^{df}\eta_{gh}\partial_{fe}h^{gh}+2\eta^{a\beta}\eta^{\alpha
e}\eta^{df}\eta_{gh}\partial_{fe}h^{gh}\\\
&+2\eta^{ad}\eta^{\alpha\beta}\eta^{ef}\eta_{gh}\partial_{fe}h^{gh}-2\eta^{a\beta}\eta^{\alpha
d}\eta^{ef}\eta_{gh}\partial_{fe}h^{gh}\Big{)}+\mathcal{O}(h^{2}\delta\phi^{2})\,,\end{split}$
(103)

and

$\begin{split}&\hat{\mathcal{R}}^{2}_{GB}=2\partial_{bc}h\,\partial^{bc}h-8\partial_{bd}h_{c}^{d}\,\partial^{bc}h+2\partial^{bc}h\,\partial_{cb}h-4\partial^{bc}h\,\partial_{cd}h_{b}^{d}+\partial_{ca}h^{ac}\,\partial_{db}h^{bd}+8\partial_{c}\partial^{c}h\,\partial_{db}h^{bd}\\\
&+8\partial_{a}\partial^{b}h^{ac}\,\partial_{db}h_{c}^{d}+8\partial^{b}\partial_{a}h^{ac}\,\partial_{db}h_{c}^{d}+4\partial^{bc}h\,\partial_{db}h_{c}^{d}+8\partial_{a}\partial^{b}h^{ac}\,\partial_{dc}h_{b}^{d}+8\partial^{b}\partial_{a}h^{ac}\,\partial_{dc}h_{b}^{d}\\\
&+8\partial_{b}\partial^{b}h^{ac}\,\partial_{d}\partial^{d}h_{ac}+4\partial_{c}\partial^{c}h_{a}^{a}\,\partial_{d}\partial^{d}h_{b}^{b}-24\partial_{a}\partial^{b}h^{ac}\,\partial_{d}\partial^{d}h_{cb}-8\partial^{b}\partial_{a}h^{ac}\,\partial_{d}\partial^{d}h_{cb}\\\
&+8\partial^{bc}h_{a}^{a}\,\partial_{d}\partial^{d}h_{cb}+8\partial_{ca}h_{bd}\,\partial^{d}\partial^{b}h^{ac}+8\partial_{cd}h_{ab}\,\partial^{db}h^{ac}-8\partial_{db}h_{ac}\,\partial^{db}h^{ac}+8\partial_{dc}h_{ab}\,\partial^{db}h^{ac}+\mathcal{O}(h^{3}).\end{split}$
(104)

The expressions, once expanded and simplified using the PN parameter
$1/c^{2}$, give rise to 1PN source terms in the near-zone, which we have
reported in Eq. (20) and Eq. (38), together with the 1PN expansion of matter
sources.

##  Appendix B Field integrals and calculation techniques

In this Appendix we present the calculation of two of the key field integrals
for EW moments. The rest of the EW integrals follow the same methodology and
reasoning. As many of these integrals involve surface integration over product
on unit vectors $\hat{n}^{ij...m}$, it is useful to convert these products to
STF products for which the following simplifying relations holds

$\hat{n}^{\langle L\rangle}\equiv\sum_{p=0}^{\lfloor
l/2\rfloor}(-1)^{p}\frac{(2l-2p-1)!!}{(2l-1)!!}\left[\hat{n}^{L-2P}\delta^{P}+\operatorname{sym}(q)\right]\,.$
(105)

Angle braces on indices define a STF tensor. We use $\lfloor l/2\rfloor$ to
denote the largest integer less than or equal to $l/2$. The expression
$\operatorname{sym}(q)$ stands for all the distinct terms which result from
permuting the indices on $\hat{n}^{L-2P}\delta^{P}$. As an example, this
relation gives

$\begin{array}[]{l}\hat{n}^{\langle
ij\rangle}=\hat{n}^{ij}-\frac{1}{3}\delta^{ij},\\\ \hat{n}^{\langle
ijk\rangle}=\hat{n}^{ijk}-\frac{1}{5}\left(\hat{n}^{i}\delta^{jk}+\hat{n}^{j}\delta^{ik}+\hat{n}^{k}\delta^{ij}\right)\end{array}$
(106)

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

$\begin{split}&\int\hat{n}_{k_{1}...k_{m}}d^{2}\Omega=0\hskip 25.0pt(m\text{
odd })\,,\\\
&\int\hat{n}_{k_{1}...k_{m}}d^{2}\Omega=[4\pi/(m+1)!!]\,\times\left[\delta_{k_{1}k_{2}}\ldots\delta_{k_{m-1}k_{m}}+\text{distinct
permutations}\right]\,\,\,(m\text{ even })\,.\end{split}$ (107)

The simplest integral that appears in the two-index EW moment calculation is
$\int(\nabla U)^{2}d^{3}x$. Applying integration-by-parts on this term gives

$\int_{\mathcal{M}}(\nabla
U)^{2}x^{ij}d^{3}x=\oint_{\partial\mathcal{M}}UU^{,k}x^{ij}d^{2}S^{k}-\int_{\mathcal{M}}UU^{,kk}x^{ij}d^{3}x-2\int_{\mathcal{M}}UU^{(i}x^{j)}d^{3}x.$
(108)

The surface integral is evaluated on the boundary of the near zone
$\partial\mathcal{M}$ as a sphere with radius $\mathcal{R}$, such that we can
write $x^{i}=\mathcal{R}\hat{n}^{i}$ and
$d^{2}S^{k}=\mathcal{R}^{2}\hat{n}^{k}$. To evaluate this term we shall expand
$U$ and $U^{,k}$ in inverse powers of $\mathcal{R}$ and ignore the integrands
that depend on $\mathcal{R}$. The $\mathcal{R}$-independent piece of this
specific surface integral contains the $\mathcal{O}(\mathcal{R}^{-4})$
expansion of $UU^{,k}$ which gives an integrand with odd number of unit
vectors and thus vanishes [see Eq. (105)]. For the first volume integral,
substituting the Laplacian of $U$ easily shows that
$-\int_{\mathcal{M}}UU^{,kk}x^{ij}d^{3}x=4\pi\sum_{A,B\neq
A}\frac{m_{A}m_{B}}{r}x_{A}^{ij}$. The contribution from the second volume
integrals is zero. This can be shown by integrating-by-parts again,
$-2\int_{\mathcal{M}}UU^{,(i}x^{j)}d^{3}x=-\oint_{\partial\mathcal{M}}U^{2}x^{(i}d^{2}S^{j)}+\int_{\mathcal{M}}U^{2}\delta^{ij}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 $\int
U^{,m}U^{,lm}d^{3}x$) is the $\int UU^{,lm}d^{3}x$ term, which can be written
as Will and Wiseman (1996):

$\begin{split}\int UU^{,lm}d^{3}x&=\sum_{A,B\neq
A}\int\frac{m_{A}m_{B}}{\|\mathbf{x}-\mathbf{x}_{B}\|}\left(\frac{3(\mathbf{x}-\mathbf{x}_{A})^{lm}}{\|\mathbf{x}-\mathbf{x_{A}}\|^{5}}-\frac{\delta^{lm}}{\|\mathbf{x}-\mathbf{x}_{A}\|^{3}}-\frac{4\pi}{3}\delta^{lm}\delta^{3}(\mathbf{x}-\mathbf{x}_{B})\right)d^{3}x\\\
&=\sum_{A,B\neq
A}\int\frac{m_{A}m_{B}}{\|\mathbf{x}-\mathbf{r}_{AB}\|}\left(\frac{3y^{\langle
lm\rangle}+\delta^{lm}}{y^{3}}-\frac{\delta^{lm}}{y^{3}}\right)d^{3}y-\sum_{A,B\neq
A}\frac{4\pi}{3}\frac{\delta^{lm}m_{A}m_{B}}{r}\\\
&-\oint\frac{1}{\|\mathbf{y}-\mathbf{r}_{AB}\|}\frac{y^{lm}}{y^{3}}\left(\mathbf{y}\cdot\mathbf{x}_{A}\right)\mathcal{R}^{2}d^{2}\Omega_{y}+...\end{split}$
(109)

where in the second line we have changed the integration variables from
$\mathbf{x}$ to $\mathbf{y}=\mathbf{x}-\mathbf{x}_{A}$ and have changed
$y^{lm}$ to its STF form. There are two cases to consider: $A=B$ and $A\neq
B$. For both, the infinite series of surface integrals vanishes as the terms
either depend on $\mathcal{R}$ or average to zero because of an odd number of
unit vectors. For the $A=B$ case, $r_{AB}=0$, so the volume integral can be
evaluated easily. When $A\neq B$, we may use of the following expansions:

$\begin{split}&\frac{1}{\|\mathbf{y}+\mathbf{r}_{AB}\|}=\sum_{l,m}\frac{4\pi}{2l+1}\frac{\left(-r_{<}\right)^{l}}{r_{>}^{ll+1}}Y_{lm}^{*}\left(\hat{\mathbf{n}}_{AB}\right)Y_{lm}(\mathbf{y})\,,\\\
&\sum_{m}\int
Y_{lm}^{*}\left(\hat{\mathbf{n}}_{AB}\right)Y_{lm}(\hat{\mathbf{y}})\hat{y}^{\left(L^{\prime}\right\rangle}d^{2}\Omega_{y}=\hat{n}_{AB}^{\langle
L\rangle}\delta^{ll^{\prime}}\,,\end{split}$ (110)

where $Y_{lm}$ are the spherical harmonics, and $r_{<(>)}$ denotes the lesser
(greater) of $r_{AB}$ and $y$. We substitute the first expansion into the
volume integral and use the second identity to evaluate the terms

$\sum_{A,B\neq A}\frac{3m_{A}m_{B}}{\|x-r_{AB}\|}\frac{y^{\langle
lm\rangle}}{y^{3}}d^{3}y=\sum_{A,B\neq A}4\pi
m_{A}m_{B}\int\sum_{m^{\prime}}\frac{3}{5}\frac{r_{<}^{2}}{r_{>}^{3}}\frac{y^{\langle
lm\rangle}}{y^{3}}Y^{*}_{2m^{\prime}}(n_{AB})Y_{2m}(y)d^{3}y\,.$ (111)

The radial integral is then evaluated by

$\int_{0}^{\mathcal{R}}\frac{r_{<}^{l}}{r_{>}^{l+1}}y^{q}dy=\frac{2l+1}{(l+q+1)(l-q)}r_{AB}^{q}$
(112)

Overall we find

$\sum_{A,B\neq A}3\frac{m_{A}m_{B}}{\|x-r_{AB}\|}\frac{y^{\langle
lm\rangle}}{y^{3}}d^{3}y=\sum_{A,B\neq A}\frac{2\pi m_{A}m_{B}}{r}n_{AB}^{lm}$
(113)

##  Appendix C Far-zone contribution to waveforms

In this Appendix, we focus on the solution of Eq. (13) with far-zone field
point and show that the contributions are beyond 1PN order. In the far-zone
region, the source terms $\mu_{s}$ and $\mu^{ij}$ 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., (Will and Wiseman, 1996)) the fields at any
intermediate distance R are given by:

$h_{\mathcal{NW}}^{\alpha\beta}(x)=\frac{4G}{c^{4}}\sum_{q=0}^{\infty}\frac{(-1)^{q}}{q!}\left(\frac{M^{\alpha\beta
k_{1}\cdots k_{q}}}{R}\right)_{,k_{1}\cdots k_{q}},\hskip
25.0pt\Phi_{\mathcal{NW}}(x)=\frac{2G}{c^{4}}\sum_{q=0}^{\infty}\frac{(-1)^{q}}{q!}\left(\frac{M_{s}^{k_{1}\cdots
k_{q}}}{R}\right)_{,k_{1}\cdots k_{q}},$ (114)

with a new set of multipole moments defined as:

$M^{\alpha\beta k_{1}\cdots
k_{q}}(\tau)\equiv\int_{\mathcal{M}}\mu^{\alpha\beta}\left(\tau,\mathbf{x}^{\prime}\right)x^{\prime
k_{1}}\cdots x^{k_{q}}d^{3}x^{\prime}\,,\hskip 25.0ptM^{k_{1}\cdots
k_{q}}(\tau)\equiv\int_{\mathcal{M}}\mu_{s}\left(\tau,\mathbf{x}^{\prime}\right)x^{\prime
k_{1}}\cdots x^{k_{q}}d^{3}x^{\prime}\,.$ (115)

Note that Eq. (114) reduces to Eqns. (14) and (36) in the limit where
$R\gg\mathcal{R}$.

As can be seen from the structure of terms in Eq. (20) and Eq. (38), the
sources in the far-zone are composed only of the field components $N$ and
$\Phi$. Finding these components by using Eq. (114) requires computing the
following multipoles:

$\begin{split}&M^{00}=m_{A}c^{2}\left[1+\frac{v^{2}_{A}}{2c^{2}}-\frac{Gm_{B}(1+\alpha_{A}^{0}\alpha_{B}^{0})}{2rc^{2}}+{\mathcal{}O}(c^{-4})\right]+(A\leftrightarrow
B)\,,\\\
&M^{00i}=m_{A}c^{2}x^{i}_{A}\left[1+\frac{v^{2}_{A}}{2c^{2}}-\frac{Gm_{B}(1+\alpha_{A}^{0}\alpha_{B}^{0})}{2rc^{2}}+{\mathcal{}O}(c^{-4})\right]+(A\leftrightarrow
B)\,,\\\
&M^{00ij}=m_{A}c^{2}x^{ij}_{A}\left[1+{\mathcal{}O}(c^{-2})\right]+(A\leftrightarrow
B)\,,\end{split}$ (116)

and

$\begin{split}&M_{s}=-m_{A}c^{2}\alpha^{0}_{A}\left[1-\frac{v^{2}_{A}}{c^{2}}-\frac{Gm_{B}\alpha_{AB}}{r\,c^{2}}+{\mathcal{}O}(c^{-4})\right]+(A\leftrightarrow
B)\,,\\\
&M_{s}^{i}=-m_{A}c^{2}x^{i}_{A}\alpha^{0}_{A}\left[1-\frac{v^{2}_{A}}{c^{2}}-\frac{Gm_{B}\alpha_{AB}}{r\,c^{2}}+{\mathcal{}O}(c^{-4})\right]+(A\leftrightarrow
B),\\\
&M^{ij}_{s}=-m_{A}c^{2}\alpha^{0}_{A}x^{ij}_{A}\left[1+{\mathcal{}O}(c^{-2})\right]+(A\leftrightarrow
B)\,.\end{split}$ (117)

These moments are plugged into Eq. (114) to give far-zone fields, which
further generate the wave equations sources. To 1PN order, it is easy to
confirm that the expressions for these fields are identical to Eq. (21) by
changing $r$ to $R$, and thus give

$\begin{split}\mu^{ij}&=\frac{G}{4\pi}\left[\frac{m^{2}}{R^{4}}+\frac{m\left(m_{A}\alpha_{A}^{0}+m_{B}\alpha_{B}^{0}\right)}{R^{4}}\right]\left(\hat{n}^{\langle
ij\rangle}-\frac{1}{6}\delta^{ij}\right)\\\ &-\frac{G\alpha
f^{\prime}(\phi_{0})}{\pi}\frac{m\left(m_{A}\alpha_{A}^{0}+m_{B}\alpha_{B}^{0}\right)}{R^{6}}\left(\hat{n}^{\langle
ij\rangle}+\frac{1}{3}\delta^{ij}\right).\end{split}$ (118)

A similar calculation for $\mu_{s}$ shows that it equals zero, and thus there