#  Inspiral-merger-ringdown waveforms in Einstein-scalar-Gauss-Bonnet
gravity within the effective-one-body formalism

Félix-Louis Julié<EMAIL_ADDRESS>Max Planck Institute for
Gravitational Physics (Albert Einstein Institute), Am Mühlenberg 1, Potsdam
14476, Germany  Lorenzo Pompili<EMAIL_ADDRESS>Max Planck
Institute for Gravitational Physics (Albert Einstein Institute), Am Mühlenberg
1, Potsdam 14476, Germany  Alessandra Buonanno<EMAIL_ADDRESS>Max Planck Institute for Gravitational Physics (Albert Einstein Institute), Am
Mühlenberg 1, Potsdam 14476, Germany  Department of Physics, University of
Maryland, College Park, MD 20742, USA

###### Abstract

Gravitational waves (GWs) provide a unique opportunity to test General
Relativity (GR) in the highly dynamical, strong-field regime. So far, the
majority of the tests of GR with GW signals have been carried out following
parametrized, theory-independent approaches. An alternative avenue consists in
developing inspiral-merger-ringdown (IMR) waveform models in specific beyond-
GR theories of gravity, by combining analytical and numerical-relativity
results. In this work, we provide the first example of a full IMR waveform
model in a beyond-GR theory, focusing on Einstein-scalar-Gauss-Bonnet (ESGB)
gravity. This theory has attracted particular attention due to its rich
phenomenology for binary black-hole (BH) mergers, thanks to the presence of
non-trivial scalar fields. Starting from the state-of-the-art, effective-one-
body (EOB) multipolar waveform model for spin-precessing binary BHs
SEOBNRv5PHM, we include theory-specific corrections to the EOB Hamiltonian,
the metric and scalar energy fluxes, the GW modes, the quasi-normal-mode (QNM)
spectrum and the mass and spin of the remnant BH. We also propose a way to
marginalize over the uncertainty in the merger morphology with additional
nuisance parameters. Interestingly, we observe that changes in the frequency
of the ringdown waveform due to the final mass and spin corrections are
significantly larger than those due to ESGB corrections to the QNM spectrum.
By performing Bayesian parameter estimation for the GW events GW190412,
GW190814 and GW230529$\\_$181500, we place constraints on the fundamental
coupling of the theory ($\sqrt{\alpha_{\mathrm{GB}}}\lesssim
0.31\leavevmode\nobreak\ \mathrm{km}$ at $90\%$ confidence). The bound could
be improved by one order of magnitude by observing a single “golden” binary
system with next-generation ground-based GW detectors. Our model can be used
to improve constraints on modifications of GR with upcoming GW observations,
and to provide forecasts for future GW detectors on the ground and in space.

##  I Introduction

The observations of gravitational waves (GWs) from coalescing binary systems
[1, 2, 3, 4, 5, 6] with the LIGO and Virgo detectors [7, 8] have for the first
time allowed us to test Einstein’s theory of General Relativity (GR) in the
highly nonlinear and dynamical regime [9, 10, 11, 12, 13]. Near-future
upgrades to the LIGO, Virgo and KAGRA [14] (LVK) interferometers, and upcoming
detectors on the ground, Einstein Telescope (ET) and Cosmic Explorer (CE) [15,
16], and in space (LISA) [17], will significantly increase the number of
detected sources, some of which will be observed with signal-to-noise ratios
(SNR) reaching thousands, allowing for tests of GR with exquisite precision
(see, e.g., Ref. [18]).

One can adopt either a theory-independent or theory-specific approach to
testing GR. Theory-independent tests modify GR waveform models by introducing
some generic, parameterized deviations, which are constrained by comparing the
waveform templates with an observed GW signal (a non-exhaustive list includes
Refs. [19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31]), or look for
consistency between the signal, or portions of the signal, and the data,
without invoking any parametrization of the deviations [32, 33, 34, 35].
Theory-specific tests (see, e.g., Refs. [36, 37, 38, 39, 40]) use instead
waveforms predicted in a particular alternative theory of gravity and aim at
estimating its underlying physical parameters. So far, the majority of tests
of GR with GWs have been carried out following a theory-independent approach.
While theory-independent tests can, in principle, constrain a wide range of
alternative theories, parameterized deviations are not unique, and
degeneracies among different deviation parameters complicate constraining
multiple of them at the same time unless specific combinations are selected
[41]. Moreover, there is no guarantee that parameterized deviations can
represent the potentially infinite landscape of beyond-GR theories (see Ref.
[42] for an attempt to address this issue). Furthermore, although such tests
are valuable in identifying possible deviations from GR, they do not
necessarily inform us about the new physics at play. Thus, it is relevant to
develop, both analytically (see, e.g., Refs. [43, 44, 45, 46, 47, 48, 49, 50,
51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69,
70, 71, 72]) and numerically (see, e.g., Refs. [73, 74, 75, 76, 77, 78, 79,
80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,
99, 100, 101, 102, 103, 104, 105]), waveform models in specific beyond-GR
theories of gravity.

One avenue to conduct theory-specific tests is to develop accurate inspiral-
merger-ringdown (IMR) waveform models in alternative theories of gravity,
combining analytical and numerical-relativity (NR) results, as done in GR
[106, 107, 108, 109, 110, 111, 112]. In this work, we provide the first
example of a full IMR waveform model featuring analytical beyond-GR
corrections. We focus on Einstein-scalar-Gauss-Bonnet (ESGB) theories, which
supplement GR with a massless scalar field coupled to the Gauss-Bonnet scalar
$\mathcal{G}=R^{\mu\nu\rho\sigma}R_{\mu\nu\rho\sigma}-4R^{\mu\nu}R_{\mu\nu}+R^{2}$,
with a coupling of the form $\ell_{\rm GB}^{2}f(\varphi)\mathcal{G}$ in the
action. These theories have attracted particular attention, thanks to their
rich phenomenology for black holes (BHs) with scalar “hair” [113, 114, 115,
116, 117, 118, 45, 119, 120, 121, 122, 123, 124, 61, 65]. In particular, for
certain functions $f(\varphi)$, BHs can undergo spontaneous [125, 126, 127,
128, 129, 130, 131] or dynamical [86, 93, 94, 132, 72] (de)scalarization
phenomena analogous to those predicted for neutron stars (NSs) in scalar-
tensor (ST) theories of gravity [44, 74, 76, 77]. Different methods have been
proposed to include dynamical scalarization in waveform models [133, 134, 135,
136]. Recent progress in ESGB gravity was made to obtain inspiral waveforms
within the post-Newtonian (PN) and effective-one-body (EOB) formalisms [61,
68, 62, 65, 137], to calculate quasi-normal-mode (QNM) frequencies of rotating
BHs at second order in a small-spin expansion [50, 63, 66] and very recently
for large spins [138, 139], and to numerically simulate binary black holes
(BBHs) [78, 80, 81, 84, 85, 83, 87, 88, 89, 91, 92, 101, 100, 99, 102, 103],
binary neutron stars (BNSs) [140, 141] and neutron star-black hole (NSBH)
binaries [104].

The most commonly used methods to build complete IMR waveform models for
compact binaries in GR are the NR surrogate [142, 143, 144], phenomenological
[145, 146, 147, 148, 149, 150] and EOB [108, 151, 152, 153, 154, 155]
frameworks. EOB waveform models combine and resum information from several
analytical approximation methods (PN, post-Minkowskian, gravitational self-
force) with a physically motivated ansatz for the merger, and BH perturbation
theory for the ringdown, and are made highly accurate through calibration with
NR simulations. There are currently two state-of-the-art families of EOB
waveform models: SEOBNR (see, e.g., Refs. [156, 157, 158, 159, 160]) and
TEOBResumS (see, e.g., Refs. [161, 162, 163, 164, 165]). Here, we focus on the
former.

In this work, we start from the state-of-the-art EOB waveform model
SEOBNRv5PHM [157], which includes effects of higher harmonics and spin
precession, and we include ESGB corrections to the EOB Hamiltonian, the GW and
scalar-energy fluxes, the spherical-harmonic GW modes, the QNM spectrum, and
the mass and spin of the remnant BH. To account for the uncertainty in merger
features that may only be determined when further NR simulations are available
in ESGB gravity, we use parametrized deviations inspired by the pSEOBNR
pipeline [25, 30, 166], interpreted as nuisance parameters. We implement our
waveform model within the pySEOBNR Python package [167], a code recently
developed to construct, calibrate and validate new EOB waveform models with
efficiency and flexibility. We refer to our new model with the name
SEOBNRv5PHM$\\_$ESGB.

In a second step, we use SEOBNRv5PHM$\\_$ESGB to provide constraints on the
coupling length $\ell_{\rm GB}$ of the dilatonic theory
$f(\varphi)=\exp(2\varphi)/4$ by analyzing the GW events GW190412, GW190814
and GW230529$\\_$181500 (to which we will refer as GW230529 for brevity)
through Bayesian inference, and compare our results to the literature. We
improve previous analyses [168, 169, 170] by including ESGB corrections to the
waveform for the entire coalescence, and as a consequence, by analyzing the
complete signal rather than restricting to its inspiral phase. Our results for
GW190814 are in agreement with Refs. [168, 170], with a $90\%$ bound on their
coupling length $\sqrt{\alpha_{\mathrm{GB}}}\equiv\ell_{\rm
GB}/(4\pi^{1/4})\lesssim 0.4\leavevmode\nobreak\ \mathrm{km}$. For GW230529,
which currently places the best constraint on this theory, we obtain
$\sqrt{\alpha_{\mathrm{GB}}}\lesssim 0.31\rm{km}$. This is consistent with
Refs. [171, 172], which obtain a $90\%$ bound on $\sqrt{\alpha_{\mathrm{GB}}}$
ranging from $0.26\leavevmode\nobreak\ \mathrm{km}$ to
$0.35\leavevmode\nobreak\ \mathrm{km}$ depending on the waveform approximant
employed. On the other hand, we obtain from GW190412 the bound
$\sqrt{\alpha_{\mathrm{GB}}}\lesssim 1.81\leavevmode\nobreak\ \mathrm{km}$,
roughly a factor $\sim 2$ smaller than previous analyses [168]. This is
consistent with the larger mass of the source of GW190412, resulting in a
higher relative SNR in the late-inspiral stage, for which we are including
additional corrections. Stacking the posteriors of the three events, we obtain
a combined bound $\sqrt{\alpha_{\mathrm{GB}}}\lesssim
0.30-0.32\leavevmode\nobreak\ \mathrm{km}$, depending on whether we include
the merger-ringdown portion of the signal or not. Our waveform templates also
allow for Bayesian–model-selection studies, to address whether the GR
hypothesis is favored over a specific modified gravity theory in describing GW
data [173]. We compute Bayes factors between the GR and ESGB hypotheses,
finding the three analyzed events to be consistent with GR.

The paper is organized as follows. In Sec. II we introduce ESGB theories and
their PN description. In Sec. III, we review the structure of the EOB
Hamiltonian in GR, and include its ESGB corrections at 3PN. In Sec. IV we
present the ESGB metric and scalar energy fluxes on circular orbits and at
relative 2PN and 2.5PN orders. In Section V, we recall the structure of the
IMR spherical-harmonic waveform modes in GR, and obtain their ESGB
corrections. In Sec. VI we present our results on waveform morphology and
Bayesian parameter estimation. In Section VII, we summarize our main
conclusions and discuss future work. Various technical details are relegated
to the appendices. In Appendix A we present a dictionary relating quantities
in the Einstein and Jordan frames. In Appendix B we recall the ESGB
corrections to the EOB metric. In Appendix C we gather factors entering the
metric modes. In Appendix D we list the coefficients of the scalar flux.
Finally, in Appendix E we provide the spherical-harmonic scalar modes at
relative 1.5PN order. Throughout this paper we use geometrical units
($c=G=1$).

##  II Einstein-scalar-Gauss-Bonnet reminder

###  II.1 ESGB and ST gravity

We study the class of theories described by the Einstein-frame action [61, 43]

$\displaystyle I=\frac{1}{16\pi}\int\\!d^{4}x$
$\displaystyle\sqrt{-g}\big{(}R-2g^{\mu\nu}\partial_{\mu}\varphi\partial_{\nu}\varphi+\ell_{\rm
GB}^{2}f(\varphi)\mathcal{G}\big{)}$ $\displaystyle+I_{\rm
m}[\Psi,\mathcal{A}^{2}(\varphi)g_{\mu\nu}]\,,$ (1)

where $R$ is the Ricci scalar, $g={\rm det}\,g_{\mu\nu}$ is the metric
determinant, and
$\mathcal{G}=R^{\mu\nu\rho\sigma}R_{\mu\nu\rho\sigma}-4R^{\mu\nu}R_{\mu\nu}+R^{2}$
is the Gauss-Bonnet scalar, where $R^{\mu}{}_{\nu\rho\sigma}$ and $R_{\mu\nu}$
are the Riemann and Ricci tensors, respectively. The integral $\int
d^{4}x\sqrt{-g}\,\mathcal{G}$ over a four-dimensional spacetime is a boundary
term [174]. As for matter fields $\Psi$, they couple minimally to the Jordan
metric $\tilde{g}_{\mu\nu}=\mathcal{A}^{2}(\varphi)g_{\mu\nu}$. The length
$\ell_{\rm GB}$ and the dimensionless functions $\mathcal{A}$ and $f$ specify
the theory. We recover ST theories in the limit $\ell_{\rm GB}=0$ or if $f$ is
constant, and GR if moreover $\mathcal{A}$ (and $\varphi$) is constant.

To model compact bodies, we use the phenomenological treatment initiated in
Refs. [175, 43] in ST theories, and describe them as point particles:

$\displaystyle I_{\rm m}\to I_{\rm m}^{\rm
pp}[g_{\mu\nu},\varphi,\\{x_{A}^{\mu}\\}]=-\sum_{A}\int
m_{A}(\varphi)\,ds_{A}\,,$ (2)

where $ds_{A}=\sqrt{-g_{\mu\nu}dx_{A}^{\mu}dx_{A}^{\nu}}$ and
$x_{A}^{\mu}[s_{A}]$ is the worldline of particle $A$. The constant GR mass
becomes a function $m_{A}(\varphi)$ that depends on the internal structure of
body $A$ and on the value of the scalar field at $x_{A}^{\mu}(s_{A})$.

In this paper, we will refer to the theory with action (1) as ESGB gravity,
but note that ST theories are included as a particular subcase.

0PN | $G_{AB}=1+\alpha_{A}^{0}\alpha_{B}^{0}$
---|---
1PN |  $\bar{\gamma}_{AB}=\frac{-2\alpha_{A}^{0}\alpha_{B}^{0}}{1+\alpha_{A}^{0}\alpha_{B}^{0}}$ $\bar{\beta}_{A}=\frac{\beta_{A}^{0}(\alpha_{B}^{0})^{2}}{2(1+\alpha_{A}^{0}\alpha_{B}^{0})^{2}}$
2PN |  $\delta_{A}=\frac{(\alpha_{A}^{0})^{2}}{(1+\alpha_{A}^{0}\alpha_{B}^{0})^{2}}$ $\epsilon_{A}=\frac{{\beta^{\prime}}_{A}^{0}(\alpha_{B}^{0})^{3}}{(1+\alpha_{A}^{0}\alpha_{B}^{0})^{3}}$ $\zeta_{AB}=\frac{{\beta}_{A}^{0}{\beta}_{B}^{0}\alpha_{A}^{0}\alpha_{B}^{0}}{(1+\alpha_{A}^{0}\alpha_{B}^{0})^{3}}$
3PN |  $\kappa_{A}=\frac{(\alpha_{B}^{0})^{4}{\beta^{\prime\prime}_{A}}^{0}}{8(1+\alpha_{A}^{0}\alpha_{B}^{0})^{4}}$ $\psi_{A}=\frac{\alpha_{A}^{0}\alpha_{B}^{0}\beta_{A}^{0}}{(1+\alpha_{A}^{0}\alpha_{B}^{0})^{3}}$ $\xi_{A}=\frac{(\alpha_{A}^{0})^{2}\alpha_{B}^{0}\beta_{A}^{0}{\beta^{\prime}}_{B}^{0}}{(1+\alpha_{A}^{0}\alpha_{B}^{0})^{4}}$ $\omega_{A}=\frac{(\alpha_{A}^{0})^{2}\beta_{A}^{0}(\beta_{B}^{0})^{2}}{(1+\alpha_{A}^{0}\alpha_{B}^{0})^{4}}$
|  $k_{\rm
tail}=-\frac{2}{3}\frac{\alpha_{-}^{2}\,\nu}{(1+\alpha_{A}^{0}\alpha_{B}^{0})^{2}}$
$k_{\rm ESGB}=-\frac{\ell_{\rm
GB}^{2}f^{\prime}(\varphi_{0})}{2M^{2}}\frac{3\alpha_{+}+m_{-}\alpha_{-}}{(1+\alpha_{A}^{0}\alpha_{B}^{0})^{4}}$
(Anti) symmetry |  $\alpha_{\pm}=\alpha^{0}_{A}\pm\alpha^{0}_{B}$ $X_{\pm}=X_{A}\pm X_{B}$ with $X_{A}=\\{\bar{\beta}_{A},\delta_{A},\epsilon_{A},\kappa_{A},\psi_{A},\xi_{A},\omega_{A}\\}$
Mean value | $\langle\bar{\beta}\rangle=\frac{m_{A}^{0}\bar{\beta}_{B}+m_{B}^{0}\bar{\beta}_{A}}{M}$
|
$\langle\delta\rangle=\frac{m_{A}^{0}\bar{\delta}_{A}+m_{B}^{0}\delta_{B}}{M}\qquad\langle\epsilon\rangle=\frac{m_{A}^{0}\epsilon_{B}+m_{B}^{0}\epsilon_{A}}{M}$
|
$\langle\kappa\rangle=\frac{m_{A}^{0}\kappa_{B}+m_{B}^{0}\kappa_{A}}{M}\qquad\langle\psi\rangle=\frac{m_{A}^{0}\psi_{B}+m_{B}^{0}\psi_{A}}{M}\qquad\langle\xi\rangle=\frac{m_{A}^{0}\xi_{A}+m_{B}^{0}\xi_{B}}{M}\qquad\langle\omega\rangle=\frac{m_{A}^{0}\omega_{A}+m_{B}^{0}\omega_{B}}{M}$
Table 1: Notations of this paper, following Refs. [43, 176, 52, 53, 61, 68].
The quantities in the first five lines are ordered by the PN level at which
they enter the conservative dynamics. Note the factor of 2 compared to Refs.
[51, 57, 58] in our (anti) symmetrization.

###  II.2 Post-Newtonian framework

We focus on the dynamics of compact binary systems on bound orbits. When the
relative orbital velocity is small and the gravitational field is weak, the
motion can be studied in the PN framework.111We denote by $n$PN the
$\mathcal{O}(v^{2n})\sim\mathcal{O}(M/r)^{n}$ corrections to the Newtonian
dynamics, with $v$ the system’s relative orbital velocity, $M$ the total mass,
and $r$ the orbital radius. To this aim, the field equations of the theory (1)
with the substitution (2) are solved iteratively around a flat metric
$g_{\mu\nu}=\eta_{\mu\nu}+\delta g_{\mu\nu}$, and a constant scalar field
background $\varphi=\varphi_{0}+\delta\varphi$, where $\varphi_{0}$ is imposed
by the binary’s cosmological environment.

From now on, a superscript $0$ denotes a quantity evaluated at
$\varphi=\varphi_{0}$. We consider a binary system with $m^{0}_{A}\geq
m^{0}_{B}$, and define its total mass, reduced mass, and the dimensionless
ratios

$\begin{gathered}M\equiv
m^{0}_{A}+m^{0}_{B},\qquad\mu\equiv\frac{m^{0}_{A}m^{0}_{B}}{M},\qquad\nu\equiv\frac{\mu}{M},\\\
m_{-}\equiv\frac{m^{0}_{A}-m^{0}_{B}}{M},\qquad
q\equiv\frac{m^{0}_{A}}{m^{0}_{B}},\qquad\hat{\ell}_{\rm
GB}\equiv\frac{\ell_{\rm GB}}{\mu}.\end{gathered}$ (3)

The functions $m_{A}(\varphi)$ and $m_{B}(\varphi)$, describing bodies $A$ and
$B$, can then be further expanded at 3PN by introducing the following
quantities:

$\displaystyle\alpha_{A}^{0}$ $\displaystyle=\frac{d\ln
m_{A}}{d\varphi}(\varphi_{0}),$ (4a) $\displaystyle\beta_{A}^{0}$
$\displaystyle=\frac{d\alpha_{A}}{d\varphi}(\varphi_{0}),$ (4b)
$\displaystyle{\beta^{\prime}}_{A}^{0}$
$\displaystyle=\frac{d\beta_{A}}{d\varphi}(\varphi_{0}),$ (4c)
$\displaystyle{\beta^{\prime\prime}}_{A}^{0}$
$\displaystyle=\frac{d{\beta^{\prime}}_{A}}{d\varphi}(\varphi_{0}),$ (4d)

and their counterparts for body $B$. We also introduce the notation
$f^{\prime}(\varphi_{0})=(df/d\varphi)_{\varphi_{0}}$.

The ESGB corrections included in this work will depend only on the theory-
dependent product $\ell_{\rm GB}^{2}f^{\prime}(\varphi_{0})$ and on ten body-
dependent parameters: the masses of bodies $A$ and $B$ and their logarithmic
derivatives (4) at infinity (i.e., at $\varphi=\varphi_{0}$), via the
combinations gathered in Table 1. We recover ST theories in the limit
$\ell_{\rm GB}^{2}f^{\prime}(\varphi_{0})=0$, and GR when moreover
$m_{A}(\varphi)$ and $m_{B}(\varphi)$ are constants: then, Eqs. (4) and their
$B$-counterparts are zero, so that only $G_{AB}=1$ is nonzero in Table 1.

Note that $k_{\rm tail}$ differs from its original definition in Ref. [68] by
a factor $\mathcal{A}_{0}^{2}$. This factor originates from a typographical
error in Ref. [58], which the author communicated to us, after it was pointed
out in Ref. [68] [cf. Sec. II.B there] that the terms in Ref. [58] depending
on $\mathcal{A}_{0}$ or $\alpha_{0}=(d\ln\mathcal{A}/\varphi)_{\varphi_{0}}$
(in the Einstein frame) should be revised. In practice, the factor
$1/\phi_{0}$ in Eq. (A3) of Ref. [58] must be replaced by $1$. The terms
depending on $\alpha_{0}$ in Ref. [58] still require revision, and they enter
the present paper via
$\tilde{\alpha}=(1+\alpha_{A}^{0}\alpha_{B}^{0})/(1+\alpha_{0}^{2})$ in Eq.
(85b).

###  II.3 The example of dilatonic ESGB gravity

The quantities (4) can be calculated once the theory and the bodies are
specified. In ST theories, they were derived numerically for NSs in Refs. [43,
177, 44, 44, 178]. They were also calculated in ESGB models for BHs, including
spontaneously and dynamically scalarized ones, in Refs. [61, 65, 72].

For the sake of illustration, we will apply our results explicitly to the
dilatonic ESGB theory222We translate between our conventions and those of
Refs. [168, 170] by setting $\varphi=\sqrt{4\pi}\phi$ and $\ell_{\rm
GB}=\sqrt{\alpha_{\rm GB}}\,4\pi^{1/4}$.

$\displaystyle
f(\varphi)=\frac{1}{4}\exp(2\varphi),\quad\mathcal{A}(\varphi)=1,$ (5)

for which the action (1) becomes invariant under the simultaneous
redefinitions $\varphi\to\varphi+\Delta\varphi$ and $\ell_{\rm GB}\to\ell_{\rm
GB}\exp(-\Delta\varphi)$, where $\Delta\varphi$ is a constant. In this theory,
the quantities (4) were derived analytically and numerically for nonspinning
BHs in Refs. [61, 65]. Given a BH $A$, they are functions of two parameters
only: $\varphi_{0}$, which we set here to zero without loss of generality,
using the symmetry of the theory; and the dimensionless ratio $\ell_{\rm
GB}/m_{A}^{0}\lesssim 0.831$, where the bound ensures the regularity of the
horizon [65].

In this paper, we will describe BBHs and NSBH binaries as follows. For BBHs,
we use the analytical, nonperturbative $(5,5)$-Padé approximant of
$\alpha_{A}^{0}$ (and $\alpha_{B}^{0}$) discussed in Sec. III.B of Ref. [65],
which exhibits remarkable agreement with its numerical counterpart. We then
infer $\beta_{A}^{0}$, ${\beta^{\prime}}_{A}^{0}$ and
${\beta^{\prime\prime}}_{A}^{0}$ (and their counterparts for BH $B$) using
Eqs. (4b)–(4d). Since $m^{0}_{A}\geq m^{0}_{B}$, it is sufficient to impose
the horizon regularity bound on BH $B$.

By contrast, Eqs. (4) were so far not calculated for NSs in the dilatonic ESGB
theory. However, we can describe NSBH binaries at leading order in the
small-$\ell_{\rm GB}$ approximation. We indeed have, for a BH $A$ [61, 65],

$\displaystyle\left.\alpha_{A}^{0}\right|_{\rm
BH}=-\frac{1}{4}\left(\frac{\ell_{\rm
GB}}{m_{A}^{0}}\right)^{2}+\mathcal{O}\left(\frac{\ell_{\rm
GB}}{m_{A}^{0}}\right)^{4}.$ (6)

At this order, the theory reduces to the shift-symmetric ESGB model [61], in
which NSs are known to carry no scalar monopole [46, 179]. Since Eq. (4a) is
proportional to the latter [43, 54], we conclude that for a NS $B$ [170],

$\displaystyle\left.\alpha_{B}^{0}\right|_{\rm
NS}=\mathcal{O}\left(\frac{\ell_{\rm GB}}{m_{B}^{0}}\right)^{4}.$ (7)

In the small-$\ell_{\rm GB}$ approximation, the contributions from
$\beta_{A}^{0}$, ${\beta^{\prime}}_{A}^{0}$ and
${\beta^{\prime\prime}}_{A}^{0}$ (and their counterparts for NS $B$) are
subleading and can be discarded, because they only enter the dynamics via the
combinations gathered in Table 1, which are already quadratic in
$\alpha_{A}^{0}$ or $\alpha_{B}^{0}$.

The ST two-body Lagrangian was derived at 1PN by Damour and Esposito-Farèse
[43], at 2PN by Mirshekari and Will [47] and at 3PN by Bernard [58]. It was
then generalized by Julié and Berti, who derived its ESGB corrections in Ref.
[61]. The resulting conservative dynamics was included within an EOB
Hamiltonian at 2PN by Julié and Deruelle [52, 53, 68], and at 3PN by Julié et
al. [68]; see also Jain et al. [67, 69] in the ST limit. We include the main
results as corrections to the SEOBNRv5PHM Hamiltonian in Sec. III.1, and
extend those works to the radiative sector in Secs. III.2, IV and V.

##  III Effective-one-body Hamiltonian and equations of motion

###  III.1 The EOB Hamiltonian

In the EOB formalism [108, 151, 152, 153, 154], the two-body Hamiltonian in
the center-of-mass frame is mapped, through canonical transformations, to the
effective Hamiltonian $H_{\rm eff}$ of a test mass in a deformed BH spacetime,
the deformation parameter being the symmetric mass-ratio $\nu$, as [108]

$\hat{H}_{\mathrm{EOB}}=\frac{H_{\mathrm{EOB}}}{\mu}=\frac{1}{\nu}\sqrt{1+2\nu\left(\hat{H}_{\mathrm{eff}}-1\right)}\,,$
(8)

where $\hat{H}_{\rm eff}=H_{\rm eff}/\mu$. For nonspinning (noS) binaries, the
motion is planar and we can use polar coordinates $(r,\phi,p_{r},p_{\phi})$.
We introduce the following reduced variables:

$\displaystyle\hat{r}=\frac{r}{M},\quad\hat{p}_{r}=\frac{p_{r}}{\mu},\quad\hat{p}_{\phi}=\frac{p_{\phi}}{M\mu},\quad\hat{t}=\frac{t}{M}.$
(9)

The effective Hamiltonian then takes the form

$\hat{H}_{\mathrm{eff}}^{\mathrm{noS}}=\frac{{H}_{\mathrm{eff}}^{\mathrm{noS}}}{\mu}=\sqrt{\hat{p}_{r_{*}}^{2}+A_{\mathrm{noS}}(\hat{r})\left[1+\frac{\hat{p}_{\phi}^{2}}{\hat{r}^{2}}+Q_{\mathrm{noS}}\left(\hat{r},\hat{p}_{r_{*}}\right)\right]}\,,$
(10)

where we use the tortoise-coordinate $\hat{p}_{r_{*}}$ instead of
$\hat{p}_{r}$ to facilitate the numerical integration of the equations of
motion close to merger. Indeed, $\hat{p}_{r}$ diverges when approaching the
zero of $A_{\mathrm{noS}}(\hat{r})$, while $\hat{p}_{r_{*}}$ is finite [180,
181]. For nonspinning binaries, the tortoise-coordinate $\hat{r}_{*}$ is
defined by

$\frac{d\hat{r}_{*}}{d\hat{r}}=\frac{1}{\xi_{\mathrm{noS}}(\hat{r})}\,,\quad\xi_{\mathrm{noS}}(\hat{r})\equiv
A_{\mathrm{noS}}(\hat{r})\sqrt{\bar{D}_{\mathrm{noS}}(\hat{r})}\,,$ (11)

and the conjugate momentum $\hat{p}_{r_{*}}$ is

$\hat{p}_{r_{*}}=\hat{p}_{r}\,\xi_{\mathrm{noS}}(\hat{r})\,.$ (12)

In the general relativistic SEOBNRv5PHM Hamiltonian [159], the potentials
$A^{\rm GR}_{\text{noS}}$ and $\bar{D}^{\rm GR}_{\text{noS}}$ are based on the
5PN results of Refs. [182, 183], which are complete modulo two unknown
constants at 5PN order, denoted $a_{6}^{\nu^{2}}$ and $\bar{d}_{5}^{\nu^{2}}$
there. In SEOBNRv5PHM, one sets $\bar{d}_{5}^{\nu^{2}}=0$ and replaces the 5PN
coefficient in $A^{\rm GR}_{\text{noS}}$, with the exception of the part
proportional to $\ln(\hat{r})$, by a function $\nu\,a_{6}(\nu)$ calibrated to
NR simulations [156]. The potential $Q^{\rm GR}_{\text{noS}}$ uses the full
5.5PN expression of Refs. [183, 184], once $\hat{p}_{r}$ is traded for
$\hat{p}_{r_{*}}$ by inverting Eq. (12) at the same order. This potential is
known to $\operatorname{\mathcal{O}}(\hat{p}_{r}^{8})$, thus accounting for
nonlocal-in-time “tail” effects (entering at 4PN in GR) to eighth order in the
binary’s orbital eccentricity. The complete potentials are given explicitly in
Ref. [159]

To include the ESGB effects up to 3PN order in the SEOBNRv5PHM Hamiltonian, we
follow the prescription of Ref. [68]:

  1. 1.

Formally replace $\hat{r}^{-1}$ by $u$, defined as

$\displaystyle u\equiv\frac{G_{AB}M}{r},$ (13)

at all PN orders in the potentials $A^{\rm GR}_{\text{noS}}$, $\bar{D}^{\rm
GR}_{\text{noS}}$ and $Q_{\text{noS}}$, where $G_{AB}$ is defined in Table 1.

  2. 2.

Supplement the resulting potentials, which are for now Taylor-expanded (Tay)
with respect to $u$, with the corrections found in Refs. [52, 68]:

$\displaystyle A^{\text{Tay}}_{\text{noS}}(u)$ $\displaystyle=A^{\rm
GR}_{\text{noS}}(u)+2(\langle\bar{\beta}\rangle-\bar{\gamma}_{AB})u^{2}+\delta\bar{a}_{3}u^{3}$
$\displaystyle+[\delta\bar{a}_{4}+\bar{a}_{4,\ln}\ln u]u^{4},$ (14a)
$\displaystyle\bar{D}^{\text{Tay}}_{\text{noS}}(u)$
$\displaystyle=\bar{D}^{\text{\rm
GR}}_{\text{noS}}(u)-2\bar{\gamma}_{AB}u+\delta\bar{d}_{2}u^{2}+\delta\bar{d}_{3}u^{3},$
(14b) $\displaystyle Q^{\text{Tay}}_{\text{noS}}(u,\hat{p}_{r_{*}})$
$\displaystyle=Q^{\rm
GR}_{\text{noS}}(u,\hat{p}_{r_{*}})+\delta\bar{q}_{1}\hat{p}_{r_{*}}^{4}u^{2}+\bar{q}_{2}\hat{p}_{r_{*}}^{6}u,$
(14c)

where $\langle\bar{\beta}\rangle$ and $\bar{\gamma}_{AB}$ are defined in Table
1, and where the explicit expressions of
$(\delta\bar{a}_{3},\,\delta\bar{a}_{4},\,\bar{a}_{4,\ln})$,
$(\delta\bar{d}_{2},\,\delta\bar{d}_{3})$ and
$(\delta\bar{q}_{1},\,\bar{q}_{2})$ are given in Appendix B.

We note that the potential $Q$ in Ref. [68] is a function of $\hat{p}_{r}$,
and it is linked to our potentials through $Q=A_{\text{noS}}Q_{\text{noS}}$.
Since it was calculated at leading PN order, our
$(\delta\bar{q}_{1},\,\bar{q}_{2})$ are the same as those presented in Ref.
[68]. The ESGB corrections in Eq. (14c) are known to
$\mathcal{O}(\hat{p}_{r}^{6})$, which amounts to including the ESGB “tail”
effects to sixth order in the eccentricity. These effects already appear at
3PN, due to (dipolar) scalar radiation that is absent in GR.

To improve the accuracy of the model against NR waveforms in the GR limit, we
resum $A_{\text{noS}}^{\text{Tay}}(u)$ using a (1,5)-Padé approximant with
respect to $u$ [152, 185, 186, 161, 156]

$A_{\text{noS}}(u)=P^{1}_{5}[A_{\text{noS}}^{\text{Tay}}(u)],$ (15)

which features a simple zero (by construction), consistently with the
Schwarzschild metric recovered in the GR, test-mass limit. We also resum
$\bar{D}_{\text{noS}}(u)$ as a (2,3)-Padé approximant [156, 187]

$\bar{D}_{\text{noS}}(u)=P^{2}_{3}[\bar{D}_{\text{noS}}^{\text{Tay}}(u)],$
(16)

which remains numerically close to the Schwarzschild value 1 in the range
$u\lesssim 1$. Note that the ESGB corrections in the $A_{\text{noS}}$ and
$\bar{D}_{\text{noS}}$ potentials are added before performing the Padé
resummations (which are the same as in GR) so that the potentials remain
qualitatively consistent with their GR counterparts as $\ell_{\rm GB}$ is
increased.

Figure 1:  Padé-resummed EOB potentials $A_{\rm noS}(u)u^{2}$ and
$\bar{D}_{\rm noS}(u)$ for a nonspinning BBH with $q=4$. We include ESGB
corrections for the dilatonic theory $f(\varphi)=\exp(2\varphi)/4$ and
$\mathcal{A}(\varphi)=1$, and vary the Gauss-Bonnet coupling
$\hat{\ell}_{\mathrm{GB}}$ between $0$ and $1$. We add the ESGB corrections
before performing the Padé resummations, so that the potentials remain
qualitatively consistent with their GR counterparts as
$\hat{\ell}_{\mathrm{GB}}$ increases. The peak of the $A_{\rm noS}(u)u^{2}$
function corresponds to the value of $u=G_{AB}M/r$ at the effective light-
ring, located at $u=1/3$ in the Schwarzschild metric (red dot).

For generic nonzero spins, the SEOBNRv5PHM Hamiltonian takes inspiration from
that of a test mass in a deformed Kerr background [153, 188, 189, 159]. It
includes spin-orbit (SO) information up to next-to-next-to-leading order
(NNLO), spin-spin (SS) information to NNLO, as well as cubic- and quartic-in-
spin terms at leading order (LO), corresponding to the 4PN order for generic,
precessing spins. For efficiency, some spin-precession effects are
analytically averaged out over circular orbits [159]. The tortoise coordinate
is also extended as

$\xi(\hat{r})=\frac{\sqrt{\bar{D}_{\mathrm{noS}}}\left(A_{\mathrm{noS}}+\chi_{+}^{2}/\hat{r}^{2}\right)}{1+\chi_{+}^{2}/\hat{r}^{2}},$
(17)

with $\chi_{+}=(m_{A}^{0}/M)\chi_{A}+(m_{B}^{0}/M)\chi_{B}$, where $\chi_{A}$
and $\chi_{B}$ are the dimensionless spins of the bodies projected along the
orbital angular momentum at Newtonian order. More details and full expressions
can be found in Refs. [159, 156, 157]. We do not add any ESGB correction in
the spin contributions to the Hamiltonian, because they are currently unknown.

The ESGB corrections to the EOB potentials in Eqs. (15) and (16) are valid for
arbitrary choices of the functions $f(\varphi)$ and $\mathcal{A}(\varphi)$
entering the action (1). They can be specified to a particular theory and
compact binary once the corresponding values of the quantities (4) are given.
As an explicit illustration, we show in Fig. 1 the potentials $A_{\text{noS}}$
and $\bar{D}_{\text{noS}}$ for a nonspinning BBH with $q=4$ in the dilatonic
ESGB case (5), varying the rescaled Gauss-Bonnet coupling [cf. Eq. (3)] in the
range $0<\hat{\ell}_{\mathrm{GB}}<1$. The value $\hat{\ell}_{\mathrm{GB}}=1$
is just below that at which the BH $B$ becomes a naked singularity. Adding
back the units, this corresponds to a coupling $\sqrt{\alpha_{\rm{GB}}}\simeq
0.044\,\mathrm{km}\leavevmode\nobreak\ (M/M_{\odot})$, in the notation of Ref.
[170]. We note that the peak of $A_{\rm noS}(u)\,u^{2}$ corresponds to the
value of $u$ at the effective light-ring. While the differences in $A_{\rm
noS}(u)$ may seem small, especially below $u\simeq 0.4$, we recall that
$\Delta A_{\rm noS}\gtrsim 10^{-4}$ can still lead to significantly different
conservative dynamics [187], as is the case here.

###  III.2 The EOB equations of motion and radiation-reaction force

For nonspinning binaries or in the spin-aligned case, the SEOBNRv5PHM
equations of motion are given by Eqs. (15) of Ref. [156]. They read:

$\displaystyle\begin{aligned} \dot{\hat{r}}&=\xi\frac{\partial\hat{H}_{\rm
EOB}}{\partial\hat{p}_{r_{*}}},\quad&\dot{\hat{p}}_{r_{*}}&=-\xi\frac{\partial\hat{H}_{\rm
EOB}}{\partial\hat{r}}+\hat{\mathcal{F}}_{r},\\\
\dot{\phi}&=\frac{\partial\hat{H}_{\rm
EOB}}{\partial\hat{p}_{\phi}},\quad&\dot{\hat{p}}_{\phi}&=\hat{\mathcal{F}}_{\phi},\end{aligned}$
(18)

where the dot indicates a derivative with respect to $\hat{t}$, while
$\hat{\mathcal{F}}_{r}=\mathcal{F}_{r}/\mu$ and
$\hat{\mathcal{F}}_{\phi}=\mathcal{F}_{\phi}/\mu$. The radiation reaction (RR)
force is computed as [154]

$\mathcal{F}_{\phi}=-\frac{\mathcal{F}}{\Omega},\quad\mathcal{F}_{r}=\mathcal{F}_{\phi}\frac{\hat{p}_{r_{*}}}{\hat{p}_{\phi}},$
(19)

where $M\Omega=\dot{\phi}$ is the reduced orbital frequency, and $\mathcal{F}$
is the energy flux radiated by the binary on quasi-circular orbits. The metric
flux is obtained by summing over the contributions of the factorized metric
modes $h_{\ell m}^{\mathrm{F}}$ [180, 186, 155, 190], which we will define in
Sec. V.1, that is

$\mathcal{F}^{\mathrm{metric}}\equiv\frac{\Omega^{2}}{8\pi}\sum_{\ell=2}^{8}\sum_{m=1}^{\ell}m^{2}\left|Rh_{\ell
m}^{\mathrm{F}}\right|^{2},$ (20)

where $R$ is the distance to the source. In ESGB theories, we have to account
for an additional contribution to the flux from the scalar field,

$\mathcal{F}=\mathcal{F}^{\mathrm{metric}}+\mathcal{F}^{\mathrm{scalar}},$
(21)

which we will discuss in Sec. IV.2. In the following, we express $\mathcal{F}$
in terms of $\Omega$, and insert the result back into the system (18).

In the spin-precessing case, the SEOBNRv5PHM model uses the same equations of
motion as in Eqs. (18) in the co-precessing frame, in which the GW radiation
resembles that of an aligned-spin binary [191, 192, 193, 194, 195], together
with separate PN-expanded evolution equations for the spins and angular
momentum inferred from $\hat{H}_{\rm EOB}$ [159, 157]. For this work, we
employ the same description of the precession dynamics as in the SEOBNRv5PHM
model and refer the reader to Ref. [159] and Sec. II of Ref. [157] for more
details. We point out that in the future, the spin and angular momentum
evolution equations should be generalized to account for their ESGB
corrections.

##  IV Gravitational radiation in ESGB gravity

In this section, we adopt the conventions of Refs. [196, 197], and define the
metric perturbation:

$\displaystyle h^{\mu\nu}=\sqrt{-g}g^{\mu\nu}-\eta^{\mu\nu},$ (22)

whose indices are manipulated with the Minkowski metric. Let us introduce a
Cartesian orthonormal basis $(\mathbf{e}_{X},\mathbf{e}_{Y},\mathbf{e}_{Z})$.
In the center-of-mass frame of a spin-aligned binary, the relative motion is
planar and we set the source position as
$r^{i}=r\left(\cos\phi,\,\sin\phi,\,0\right)$. We locate the observation point
with spherical coordinates, such that
$R^{i}=R\left(\sin\Theta\,\cos\Phi,\,\sin\Theta\,\sin\Phi,\,\cos\Theta\right)$.
When $R\gg r$, the metric waveform is the transverse, traceless part of the
metric perturbation,

$\displaystyle h_{ij}^{\rm
TT}=\left(P^{k}_{i}P^{l}_{j}-\frac{1}{2}P^{kl}P_{ij}\right)\,h_{kl}+\mathcal{O}\left(\frac{1}{R^{2}}\right)\,,$
(23)

where $P^{ij}=\delta^{ij}-N^{i}N^{j}$ and $N^{i}=R^{i}/R$. Finally, we
introduce the polarization vectors $P^{i}=-\partial_{\Phi}N^{i}/\sin\Theta$
and $Q^{i}=\partial_{\Theta}N^{i}$ and define:

$\displaystyle h_{+}$
$\displaystyle=\frac{1}{2}\left(P^{i}P^{j}-Q^{i}Q^{j}\right)h_{ij}^{\rm TT},$
(24a) $\displaystyle h_{\times}$
$\displaystyle=\frac{1}{2}\left(P^{i}Q^{j}+Q^{i}P^{j}\right)h_{ij}^{\rm TT},$
(24b)

which can in turn be decomposed in spin-weighted spherical harmonics of weight
$-2$,

$\displaystyle
h_{+}-ih_{\times}=\sum_{\ell=2}^{\infty}\sum_{m=-\ell}^{\ell}{}_{-2}Y_{\ell
m}(\Theta,\Phi)\,h_{\ell m}\,.$ (25)

###  IV.1 Spherical-harmonic metric modes at 2PN

In ST theories, the metric waveform (23) was derived for nonspinning binaries
at 0PN by Damour and Esposito-Farèse [43] and at 2PN by Lang [48, 49]
relatively to the quadrupolar order. The corresponding metric modes $h_{\ell
m}$ were then inferred by Sennett et al. in Ref. [51] for circular orbits. The
corrections beyond ST from the Gauss-Bonnet coupling $\ell_{\rm GB}$ were
calculated at leading PN order by Shiralilou et al. [198, 62], but they enter
at 3PN. This means that the metric modes are complete only at 2PN in ST and
ESGB gravity, an order at which their expressions are the same. Rather, they
differ by the specific values taken by the quantities (4).

Yet, the results in Ref. [51] were presented using another, Jordan-frame
formulation of ST theories based on a set of Brans-Dicke-inspired parameters.
To recover the conventions of the present paper, we must proceed as follows:

  1. 1.

Translate the parameters in Ref. [51] in terms of the quantities (4). The
conversion is detailed in Appendix A.

  2. 2.

Notice that Ref. [51] uses a coordinate system $\\{{x^{\prime}}^{\mu}\\}$ such
that the Jordan metric $\tilde{g}_{\mu\nu}=\mathcal{A}^{2}(\varphi)g_{\mu\nu}$
is Minkowski at infinity, $\tilde{g}^{\prime}_{\mu\nu}\to\eta_{\mu\nu}$. By
contrast, we use here coordinates $\\{x^{\mu}\\}$ such that
$g_{\mu\nu}\to\eta_{\mu\nu}$. The coordinate systems are thus related by

$\displaystyle{x^{\prime}}^{\mu}=\mathcal{A}_{0}x^{\mu},$ (26)

with $\mathcal{A}_{0}=\mathcal{A}(\varphi_{0})$. This means that $R^{\prime}$,
$t^{\prime}$, and the orbital velocity $\omega^{\prime}$ entering Ref. [51]
must be rescaled as $R^{\prime}=\mathcal{A}_{0}R$,
$t^{\prime}=\mathcal{A}_{0}t$, and $\omega^{\prime}=\Omega/\mathcal{A}_{0}$.

  3. 3.

Observe that Ref. [51] studies the perturbation of
$(\phi/\phi_{0})\,\tilde{g}_{\mu\nu}^{\prime}$ where
$\phi=\mathcal{A}(\varphi)^{-2}$ is the scalar field in the Jordan frame. This
conformal metric identifies to our $g_{\mu\nu}$ given the coordinate change
above. The perturbation (22) also differs from that of Ref. [51] by an overall
minus sign. However, their polarization vectors are $\tilde{P}^{i}=Q^{i}$ and
$\tilde{Q}^{i}=-P^{i}$, so that our $h_{+}$ and $h_{\times}$ are the same.

In harmonic coordinates $(t,R^{i})$ such that $\partial_{\mu}h^{\mu\nu}=0$,
the metric waveform exhibits a logarithmic correction to its $\sim 1/R$ decay,
due to tail effects at 1.5PN [48, 49]. To eliminate it, we can introduce the
radiative time

$\displaystyle T=t-2M\ln\left(\frac{R}{R_{0}}\right),$ (27)

where $R_{0}$ is an arbitrary length. The ESGB metric modes then read, at 2PN
and in the radiative coordinate system $(T,R^{i})$:

$\displaystyle h_{\ell m}(T,R)=\frac{2M\nu
x}{R}\sqrt{\frac{16\pi}{5}}\hat{H}_{\ell m}e^{-im\psi},$ (28)

where

$\displaystyle x=(G_{AB}M\Omega)^{2/3},$ (29)

with, recall, $M\Omega=\dot{\phi}$ [cf. below Eq. (18)], and where

$\displaystyle\psi$
$\displaystyle=\phi-\Phi+\frac{\pi}{2}-\frac{2x^{3/2}}{1+\alpha_{A}^{0}\alpha_{B}^{0}}\left(\gamma_{\rm
E}-\frac{11}{12}+\ln\left(4\Omega R_{0}\right)\right),$ (30)

where $\gamma_{E}$ is Euler’s constant. The effective phase (30) was first
introduced in GR to factor out $\ln(\Omega)$-dependent tail corrections to the
modes [199], and it is now dressed-up by the scalar field effects. As for the
dimensionless quantities $\hat{H}_{\ell m}$, we gather their explicit
expressions in Appendix C. Note that the orbit being planar, we have by
symmetry [196]

$\displaystyle h_{\ell,-m}=(-1)^{\ell}(h_{\ell m})^{*},$ (31)

where $*$ denotes complex conjugation. Thus, it is sufficient to present the
modes with $m\geq 0$.

In Sec. V.1, we shall resum the modes (28) by rewriting them in a factorized
form $h_{\ell m}^{\mathrm{F}}$, allowing in turn to compute the metric flux
[43, 200] (cf. also Refs. [60, 56, 62]),

$\displaystyle\mathcal{F}^{\rm
metric}=\frac{R^{2}}{32\pi}\int_{R\to\infty}\\!\\!\\!\\!\\!\\!d^{2}\Omega\,\left(\frac{\partial
h_{ij}^{\rm TT}}{\partial U}\right)^{2}\,,$ (32)

where $U=T-R$. This flux formula returns Eq. (20) on circular orbits, given
the orthogonality of the spherical harmonics $\int
d^{2}\Omega\,({}_{-2}Y_{\ell
m})({}_{-2}Y_{\ell^{\prime}m^{\prime}})^{*}=\delta_{\ell\ell^{\prime}}\delta_{mm^{\prime}}$,
the contribution of each mode $m<0$ being the same as that of its counterpart
$-m$.

###  IV.2 Scalar flux at relative 2.5PN order

In ST theories, the scalar waveform

$\displaystyle\delta\varphi=\varphi-\varphi_{0}+\mathcal{O}\left(\frac{1}{R^{2}}\right)$
(33)

was derived for nonspinning binaries at 2PN by Lang [49] and at 2.5PN by
Bernard et al. [64] relatively to the (scalar) dipolar order. The
decomposition of $\delta\varphi$ on spherical harmonics was performed in Ref.
[64], but the associated scalar modes are presented there through 1.5PN only.
In this paper, we thus start from the scalar flux itself, which was derived at
2.5PN in Ref. [64]. As for the corrections beyond ST from the Gauss-Bonnet
coupling $\ell_{\rm GB}$ to the scalar flux, they were calculated at leading
PN order by Shiralilou et al. [198, 62], but they enter at 3PN, an order at
which even the ST subcase is unknown.

Once again, the results in Refs. [49, 64] were presented using the Jordan-
frame formulation of the theory. We recover our conventions as follows:

  1. 1.

Repeat steps 1 and 2 from Sec. IV.1.

  2. 2.

Notice that Refs. [49, 64] study the scalar perturbation
$\Psi=\phi/\phi_{0}-1$, where $\phi=\mathcal{A}(\varphi)^{-2}$ is the scalar
field in the Jordan frame. Our $\delta\varphi$ is thus related to $\Psi$
through $\delta\varphi=-\Psi/(2\alpha_{0})$, where $\alpha_{0}$ is defined in
Appendix A.

  3. 3.

Since $t^{\prime}=\mathcal{A}_{0}t$ and the total flux is identified to the
mechanical energy loss $\mathcal{F}=-dH/dt$ of the system, we must rescale the
scalar fluxes of Refs. [49, 64] as $\mathcal{F}^{\rm
scalar}=\mathcal{A}_{0}^{2}\,\tilde{\mathcal{F}}^{\rm scalar}$.

As a result of this procedure, the scalar flux formula in Refs. [49, 64]
becomes

$\displaystyle\mathcal{F}^{\rm scalar}=\frac{R^{2}}{4\pi}\int
d^{2}\Omega\left(\frac{\partial\delta\varphi}{\partial U}\right)^{2},$ (34)

which coincides with our Einstein-frame definition, cf. also Refs. [43, 56,
60]. For circular orbits, the ESGB scalar flux then reads:

$\displaystyle\mathcal{F}^{\rm
scalar}=\frac{\nu^{2}x^{5}}{\left(1+\alpha_{A}^{0}\alpha_{B}^{0}\right)^{2}}$
$\displaystyle\Big{[}x^{-1}f_{-1\rm{PN}}+f_{0\rm{PN}}+x^{1/2}f_{0.5\rm{PN}}$
$\displaystyle+xf_{1\rm{PN}}+x^{3/2}f_{1.5\rm{PN}}+\cdots\Big{]},$ (35)

where $x$ was defined in Eq. (29), and where the lengthy expressions of the
coefficients $f_{i\rm{PN}}$ are relegated to Appendix D. To make the
comparison with the metric flux transparent, the labels $i$PN are here
relative to the quadrupolar order.

We point out that $f_{-1\rm{PN}}$ is proportional to
$(\alpha_{A}-\alpha_{B})^{2}$. Thus, when restricted to the dipolar order, the
scalar flux vanishes for symmetric binaries. In Fig. 2, we show the behavior
of the scalar flux (IV.2) truncated at successive PN orders, as a function of
$x$. As an illustration, we consider a BBH with $q=4$ and
$\hat{\ell}_{\mathrm{GB}}=1$ in the dilatonic ESGB theory (5). The curves
terminate at merger, taken as the peak of the $(2,2)$ mode amplitude when all
PN orders of the scalar flux contribute to the binary’s evolution. We also
show the corresponding number of GW cycles. The ratio between the scalar flux
truncated at different PN orders and the leading order -1PN term highlights
the poor convergence of the low PN order contributions in the late inspiral,
which seems, however, improved when including higher-order PN corrections. We
leave the comparison of potential resummation methods for the scalar flux to
future work, when more NR waveforms become available in ESGB gravity.

Figure 2: Top panel: The PN-expanded scalar flux given in Eq. (IV.2),
truncated at successive PN orders. We consider a BBH with $q=4$ and
$\hat{\ell}_{\mathrm{GB}}=1$ in the dilatonic theory
$f(\varphi)=\exp(2\varphi)/4$ and $\mathcal{A}(\varphi)=1$. The curves
terminate at merger, taken as the peak of the $(2,2)$ mode when all PN orders
of the scalar flux contribute to the binary’s evolution, and the corresponding
number of GW cycles. Bottom panel: Ratio between the scalar flux truncated at
different PN orders and the leading -1PN term.

For completeness, we present in Appendix E the decomposition of the scalar
waveform in spherical harmonics, once translated in our Einstein-frame
conventions. As explained below Eq. (33), we do not infer the scalar flux from
these modes, because the result would be one PN order lower than in Eq.
(IV.2).

##  V EOB inspiral-merger-ringdown metric waveform in ESGB gravity

In the EOB framework, the metric modes defined in Eq. (25) are decomposed into
inspiral-plunge and merger-ringdown modes. The SEOBNRv5PHM model includes the
modes $(\ell,|m|)=$ (2,2), (2,1), (3,3), (3,2), (4,4), (4,3) and (5,5) in the
co-precessing frame. We have:

$h_{\ell m}(\hat{t})=\begin{cases}h_{\ell m}^{\text{insp-plunge
}}(\hat{t}),&\hat{t}<\hat{t}_{\text{match }}^{\ell m}\\\ h_{\ell
m}^{\text{merger-RD }}(\hat{t}),&\hat{t}>\hat{t}_{\text{match }}^{\ell
m}\end{cases},$ (36)

where $\hat{t}_{\text{match }}^{\ell m}$ is the peak time of the $(2,2)$-mode
amplitude ($\hat{t}_{\text{peak }}^{22}$) for all the modes except $(5,5)$,
for which it is taken as the peak time of the $(2,2)$ mode minus $10M$ [156].
From now on, we shall refer to $\hat{t}_{\text{peak }}^{22}$ as the merger
time. Note that the other modes with $\ell\leq 8$ are known in their
factorized form (cf. below), but they are only used to compute the fluxes
until merger.

In GR, $\hat{t}_{\text{peak}}^{22}$ is suitably chosen to agree with the
corresponding time in NR simulations. Specifically, in the SEOBNRv5PHM model,
one defines

$\hat{t}_{\text{peak
}}^{22}=\hat{t}_{\mathrm{ISCO}}+\Delta\hat{t}_{\mathrm{NR}}$ (37)

where $\hat{t}_{\rm{ISCO}}$ is the time at which $u=u_{\rm{ISCO}}$ reaches the
inverse radius of the Kerr innermost stable circular orbit (ISCO) [201] with
the mass and spin of the remnant BH, as predicted by NR fitting formulas [202,
203]. The calibration parameter $\Delta\hat{t}_{\rm NR}$ is determined via
comparisons to NR simulations [156]. In this work, we also use Eq. (37),
formally replacing $u$ by its ESGB expression (13), to determine the merger-
ringdown attachment time. We note that $\hat{t}_{\mathrm{ISCO}}$ is affected
by ESGB corrections, as they give rise to differences in the orbital dynamics,
but we don’t include additional corrections to $\Delta\hat{t}_{\rm NR}$, which
should be determined from NR simulations in ESGB gravity.

###  V.1 Inspiral-plunge $h_{\ell m}$ modes

The inspiral-plunge EOB modes can be written as

$h_{\ell m}^{\text{insp-plunge }}=h_{\ell m}^{\mathrm{F}}N_{\ell m}\,,$ (38)

where $h_{\ell m}^{\mathrm{F}}$ is a factorized, resummed form of the PN-
expanded metric modes for aligned-spin binaries on circular orbits [180, 155,
190], while $N_{\ell m}$ is the non-quasi-circular (NQC) correction [204]
which does not enter the flux (20) in the SEOBNRv5PHM model. More precisely,
the factorized metric modes are defined as

$h_{\ell m}^{\mathrm{F}}=h_{\ell m}^{\text{N}}\hat{S}_{\mathrm{eff}}T_{\ell
m}f_{\ell m}e^{i\delta_{tm}}.$ (39)

In GR, the first factor $h_{\ell m}^{\text{N,GR}}$ is the mode at leading PN
order (called Newtonian). Its explicit expression can be found in Refs. [155,
190] for all $(\ell,m)$, with an overall minus sign compared to our
conventions. The second factor $\hat{S}_{\text{eff}}^{\rm GR}$ takes
inspiration from the Regge-Wheeler-Zerilli equation sourced by a test mass
$\mu$, on circular orbits of the Schwarzschild spacetime of mass $M$ [155].
Depending on the parity of $\ell+m$, it is either identified to the (reduced)
effective energy of the EOB framework, or to the orbital angular momentum,
normalized by its Newtonian value on circular orbits:

$\hat{S}_{\text{eff}}^{\rm
GR}=\left\\{\begin{array}[]{ll}\hat{H}_{\text{eff}}^{\rm
GR},&\quad\ell+m\text{ even},\\\ \\\ \sqrt{x_{\rm
GR}}{\hat{p}_{\phi}},&\quad\ell+m\text{ odd}.\end{array}\right.$ (40)

Here, $x_{\rm GR}=(M\Omega)^{2/3}$ and we recall that $M\Omega=\dot{\phi}$ was
introduced below Eq. (18). The third factor $T_{\ell m}^{\rm GR}$ resums
leading logarithms entering the tail contributions [205]:

$\displaystyle T_{\ell m}^{\rm GR}=\frac{\Gamma\left(\ell+1-2i\hat{k}_{\rm
GR}\right)}{\Gamma(\ell+1)}e^{\pi\hat{k}_{\rm GR}}e^{2i\,\hat{k}_{\rm
GR}\ln(2m\Omega R_{0})},$ (41)

where $\Gamma(.)$ is the Euler gamma function, while $\hat{k}_{\rm
GR}=m\,H_{\rm EOB}^{\rm GR}\Omega$ is proportional to the binary’s resummed
EOB energy. The constant $R_{0}$ is the same as that entering Eq. (30). In GR,
it is usually set to $R_{0}=2M/\\!\sqrt{e}$, where $e\simeq 2.718$ is Euler’s
number, to reproduce the correct test-particle limit waveforms [190]. The
remaining factor depends on an amplitude $f_{\ell m}^{\rm GR}$ and a phase
$\delta_{\ell m}^{\rm GR}$, whose PN expansions in $x_{\rm GR}$ can be
determined by identifying order-by-order the mode $h^{\rm GR}_{\ell m}$ to its
resummation $h_{\ell m}^{\text{F,GR}}$ on circular orbits.

The variable $x_{\rm GR}$ can be traded for $y_{\rm GR}=\big{(}H_{\rm
EOB}^{\rm GR}\Omega\big{)}^{2/3}$ in $\delta_{\ell m}^{\rm GR}$ to gather
relativistic corrections sourced by the binary’s Arnowitt-Deser-Misner (ADM)
mass [155]. For $m$ even, or for nonspinning binaries, $f_{\ell m}^{\rm GR}$
is further resummed as $f_{\ell m}^{\rm GR}=\big{(}\rho_{\ell m}^{\rm
GR}\big{)}^{\ell}$ to reduce the magnitude of the 1PN coefficient, which grows
linearly with $\ell$.

By contrast, for $m$ odd, the spins deserve a particular treatment. Indeed,
the procedure above yields spin contributions to $\delta_{\ell m}^{\rm GR}$
inversely proportional to $m_{-}$ [cf. Eq. (3)], thus making $h_{\ell
m}^{\text{F,GR}}$ ill-defined in the equal-mass case. This issue traces back
to $h_{\ell m}^{\rm GR}$ being zero when $m_{-}=0$ at Newtonian order, while
some of its spin contributions, entering at subleading PN orders, are not. We
thus adopt the treatment of Refs. [190, 206, 207], and redefine the amplitude
and phase as follows:

$\displaystyle\bar{f}_{\ell m}^{\rm GR}=\left\\{\begin{array}[]{ll}(\rho_{\ell
m}^{\rm GR})^{\ell},&\quad m\text{ even},\\\ \\\ (\rho_{\ell
m}^{\text{noS,GR}})^{\ell}+\bar{f}_{\ell m}^{\text{S,GR}},&\quad m\text{
odd},\end{array}\right.$ (45) $\displaystyle\bar{\delta}_{\ell m}^{\rm
GR}=\left\\{\begin{array}[]{ll}\delta_{\ell m}^{\rm GR},&\quad m\text{
even},\\\ \\\ \delta_{\ell m}^{\rm noS,GR},&\quad m\text{
odd},\end{array}\right.$ (49)

where $\rho_{\ell m}^{\text{noS,GR}}$ and $\delta_{\ell m}^{\rm noS,GR}$ are
inferred from the nonspinning limit of $f_{\ell m}^{\rm GR}$ and $\delta_{\ell
m}^{\rm GR}$, while the additive constant $\bar{f}_{\ell m}^{\rm S,GR}$
reabsorbs their spin contributions, by expanding the (complex) exponential in
Eq. (39). The amplitude and phase being replaced by their “barred” versions,
Eq. (39) has a non-ambiguous equal-mass limit, because the terms inversely
proportional to $m_{-}$ only enter $\bar{f}_{\ell m}^{\rm GR}$ (and not
$\bar{\delta}_{\ell m}^{\rm GR}$), which is multiplied by $h_{\ell
m}^{\text{N,GR}}\propto m_{-}$.

The SEOBNRv5PHM model also incorporates information from second-order self-
force (2GSF) [208, 209] in the GW mode amplitudes and radiation-reaction
force, by introducing new calibration parameters in the $\rho_{\ell m}^{\rm
GR}$ coefficients and matching them to the 2GSF energy-flux multipolar data
for quasicircular nonspinning BBHs [158]. The explicit expressions of
$\bar{f}_{\ell m}^{\rm GR}$ and $\bar{\delta}_{\ell m}^{\rm GR}$ used in the
SEOBNRv5PHM model are provided in Appendix B of Ref. [156].

The steps described above can now be generalized to ESGB theories, starting
from the 2PN metric modes discussed in Sec. IV.1 [cf. also Appendix C]. To do
so, we must Taylor-expand $\hat{H}_{\rm eff}$ and $\hat{H}_{\rm EOB}$ [defined
in Sec. III.1] at 1PN, as well as $\hat{p}_{\phi}$ at 2PN on circular orbits.
Now, observe that when $m$ is odd, the modes $h_{\ell m}$ defined in Eq. (28)
and Appendix C are imaginary and proportional to $m_{-}$ at Newtonian order.
However, some real ESGB corrections, entering at subleading PN orders, are
not. The situation is reminiscent of that with spins, and thus, we shall treat
ESGB effects in a similar fashion.

More explicitly, we include the ESGB corrections to the metric modes up to 2PN
within SEOBNRv5PHM as follows:

  1. 1.

In the first three factors $h_{\ell m}^{\text{N,GR}}$,
$\hat{S}_{\text{eff}}^{\rm GR}$ and $T_{\ell m}^{\rm GR}$ of Eq. (39),
formally replace the effective and EOB energies by their ESGB expressions
$\hat{H}_{\text{eff}}$ and $\hat{H}_{\text{EOB}}$ [cf. Sec. III.1], and
$x_{\rm GR}$ by $x$ [cf. Eq. (29)].

  2. 2.

For $\ell+m$ odd, normalize $\hat{S}_{\text{eff}}^{\rm GR}(x,\hat{p}_{\phi})$
by its Newtonian value on circular orbits in ESGB gravity, that is

$\displaystyle\hat{S}_{\rm
eff}=\frac{\sqrt{x}\hat{p}_{\phi}}{1+\alpha_{A}^{0}\alpha_{B}^{0}},\qquad\ell+m\text{
odd}.$ (50)

  3. 3.

In $\rho_{\ell m}^{\rm noS,GR}$ and $\delta_{\ell m}^{\rm noS,GR}$ and in the
nonspinning parts of $\rho_{\ell m}^{\rm GR}$ and $\delta_{\ell m}^{\rm GR}$,
formally replace $x_{\rm GR}$ by $x$ [cf. Eq. (29)], and $y_{\rm GR}$ by

$\displaystyle y=\big{(}H_{\text{EOB}}\Omega\big{)}^{2/3}.$ (51)

  4. 4.

Deform the resulting amplitude and phase by ESGB corrections:

$\displaystyle\bar{f}_{\ell m}=\left\\{\begin{array}[]{ll}(\rho_{\ell m}^{\rm
GR}+\delta\rho_{\ell m})^{\ell},&\quad m\text{ even},\\\ \\\ (\rho_{\ell
m}^{\text{noS,GR}})^{\ell}+\bar{f}_{\ell m}^{\text{S,GR}}+\delta\bar{f}_{\ell
m},&\quad m\text{ odd},\end{array}\right.$ (55)
$\displaystyle\bar{\delta}_{\ell m}=\left\\{\begin{array}[]{ll}\delta_{\ell
m}^{\rm GR}+\delta\delta_{\ell m},&\quad m\text{ even},\\\ \\\ \delta_{\ell
m}^{\rm noS,GR},&\quad m\text{ odd}.\end{array}\right.$ (59)

At 2PN, the modes with $\ell\leq 6$ contribute, and we find that the nonzero
ESGB corrections read:

$\displaystyle\delta\rho_{22}$
$\displaystyle=-\frac{1}{3}x\left(\bar{\gamma}_{AB}+2\langle\bar{\beta}\rangle\right)+x^{2}\left[\frac{1}{504}\bigg{(}-8\langle\bar{\beta}\rangle\left(14\bar{\gamma}_{AB}+42\bar{\beta}_{+}-3\right)\right.$
$\displaystyle+77\bar{\gamma}_{AB}^{2}-352\bar{\gamma}_{AB}+84\left(\bar{\beta}_{+}^{2}-\bar{\beta}_{-}^{2}+\langle\delta\rangle+\langle\epsilon\rangle\right)-112\langle\bar{\beta}\rangle^{2}\bigg{)}$
$\displaystyle+\frac{\nu}{252}\bigg{(}84\left(\zeta_{AB}\bar{\gamma}_{AB}-\delta_{+}\right)-42\left(\bar{\gamma}_{AB}^{2}+\epsilon_{+}\right)-389\bar{\gamma}_{AB}+168\bar{\beta}_{+}^{2}$
$\displaystyle+378\bar{\beta}_{+}+252\zeta_{AB}-232\langle\bar{\beta}\rangle\bigg{)}\Bigg{]}+\mathcal{O}\left(x^{5/2}\right),$
(60a) $\displaystyle\delta\rho_{32}$
$\displaystyle=-\frac{4}{9}x\left(\bar{\gamma}_{AB}+\langle\bar{\beta}\rangle\right)+\mathcal{O}\left(x^{3/2}\right),$
(60b) $\displaystyle\delta\rho_{44}$
$\displaystyle=-\frac{1}{3}x\left(\bar{\gamma}_{AB}+2\langle\bar{\beta}\rangle\right)+\mathcal{O}\left(x^{3/2}\right),$
(60c) $\displaystyle\delta\rho_{42}$
$\displaystyle=-\frac{1}{3}x\left(\bar{\gamma}_{AB}+2\langle\bar{\beta}\rangle\right)+\mathcal{O}\left(x^{3/2}\right),$
(60d)
$\displaystyle\delta\delta_{22}=\frac{1}{24}y^{3/2}\left[2\left(\alpha_{+}^{2}-\alpha_{-}^{2}\right)-9\nu\alpha_{-}^{2}\right]+\mathcal{O}\left(y^{5/2}\right),$
(61) $\displaystyle\delta\bar{f}_{21}$
$\displaystyle=-\frac{1}{6}x\left(5\bar{\gamma}_{AB}+4\langle\bar{\beta}\rangle\right)-\frac{i\nu}{3}y^{3/2}\left(\alpha_{-}^{2}+\frac{\alpha_{-}\alpha_{+}}{m_{-}}\right)+\mathcal{O}\left(x^{2}\right),$
(62a) $\displaystyle\delta\bar{f}_{33}$
$\displaystyle=-x\left(\bar{\gamma}_{AB}+2\langle\bar{\beta}\rangle\right)$
$\displaystyle+iy^{3/2}\left[\frac{3}{40}\left(\alpha_{+}^{2}-\alpha_{-}^{2}\right)-\frac{2\nu}{9}\left(\alpha_{-}^{2}-\frac{\alpha_{-}\alpha_{+}}{m_{-}}\right)\right]+\mathcal{O}\left(x^{2}\right),$
(62b) $\displaystyle\delta\bar{f}_{31}$
$\displaystyle=-x\left(\bar{\gamma}_{AB}+2\langle\bar{\beta}\rangle\right)$
$\displaystyle+iy^{3/2}\left[\frac{1}{40}\left(\alpha_{+}^{2}-\alpha_{-}^{2}\right)-\frac{2\nu}{3}\left(5\alpha_{-}^{2}-\frac{\alpha_{-}\alpha_{+}}{m_{-}}\right)\right]+\mathcal{O}\left(x^{2}\right).$
(62c)

We do not add any ESGB correction in the spin contributions to the phase and
amplitude, because they are currently unknown.

The remaining $N_{\ell m}$ factor in the inspiral-plunge modes (38) is the NQC
correction and reads

$\displaystyle N_{\ell m}$
$\displaystyle=\left[1+\frac{{\hat{p}}_{r_{*}}^{2}}{(\hat{r}M\Omega)^{2}}\left(a_{1}^{h_{\ell
m}}+\frac{a_{2}^{h_{\ell m}}}{{\hat{r}}}+\frac{a_{3}^{h_{\ell
m}}}{{\hat{r}}^{3/2}}\right)\right]$ (63)
$\displaystyle\quad\times\exp\left[i\left(b_{1}^{h_{\ell
m}}\frac{{\hat{p}}_{r_{*}}}{\hat{r}M\Omega}+b_{2}^{h_{\ell
m}}\frac{{\hat{p}}_{r_{*}}^{3}}{\hat{r}M\Omega}\right)\right].$

This factor guarantees that the modes’ amplitude and frequency agree with NR
input values, given by NR fits, at the matching point $\hat{t}^{\ell
m}_{\text{match}}$. In practice, one first solves for the binary’s trajectory
to infer $h_{\ell m}^{\rm F}(\hat{t})$. Second, one fixes the 5 constants
$(a_{1}^{h_{\ell m}}$, $a_{2}^{h_{\ell m}}$, $a_{3}^{h_{\ell m}}$,
$b_{1}^{h_{\ell m}}$, $b_{2}^{h_{\ell m}})$ by requiring that the amplitude of
the EOB modes $\left|h_{\ell m}^{\text{insp-plunge
}}\left(\hat{t}_{\text{match }}^{\ell m}\right)\right|$ and its first two
derivatives, and the frequency of the EOB modes $\omega_{\ell m}^{\text{insp-
plunge }}\left(\hat{t}_{\text{match }}^{\ell m}\right)$ and its first
derivative, are the same as that of the NR modes at the matching point
$\hat{t}_{\text{match }}^{\ell m}$ [207, 210, 156]. In particular, the NQC
correction enforces that the (2,2)-mode’s derivative vanishes (i.e., it peaks)
at the merger time defined in Eq. (37). Since the factorized modes include
ESGB corrections, the five constants above are changed so that the inspiral-
plunge waveform (38) still matches the NR predictions in GR at $\hat{t}_{\ell
m}^{\rm match}$. In Sec. V.2, we shall discuss how to relax this assumption.

Figure 3: Top panel: Fluxes in GR and ESGB gravity normalized by the leading
order GR contribution, considering only the ESGB corrections to the metric, or
also adding the scalar flux. We consider a BBH with $q=4$ and
$\hat{\ell}_{\mathrm{GB}}=1$ in the dilatonic theory
$f(\varphi)=\exp(2\varphi)/4$ and $\mathcal{A}(\varphi)=1$. Bottom panel:
Fractional difference between the ESGB fluxes and the GR flux. We note that
the difference in the metric flux at low frequencies is constant with $x$
(e.g., at leading PN order), and is due to the redefinition of $x$ as in Eq.
(29), while the difference when including the scalar flux increases at low
frequencies as $x^{-1}$, due to the dipolar -1PN term.

In Fig. 3, we show the GR and ESGB fluxes normalized by the leading order GR
contribution

$\bar{\mathcal{F}}=\frac{\mathcal{F}}{\frac{32}{5}\nu^{2}x_{\rm GR}^{5}},$
(64)

considering only the ESGB corrections to the metric flux, or also adding the
scalar flux. As in Fig. 2, we consider a binary with $q=4$ and
$\hat{\ell}_{\mathrm{GB}}=1$ in the dilatonic theory (5). Looking at the
fractional difference between the ESGB fluxes and the GR flux, we note that
the difference in the metric flux at low frequencies is constant with $x$
(e.g., at leading PN order), and is due to the replacement of
$x_{\mathrm{GR}}$ with $x$ as in Eq. (29), while the difference when including
the scalar flux increases at low frequencies as $x^{-1}$, due to the dipolar
-1PN term.

###  V.2 Merger-ringdown $h_{\ell m}$ modes

For the merger-ringdown waveform in ESGB, we assume the same phenomenological
ansätze as in GR [211, 210, 212, 156], but introduce an appropriate
parameterization for the modes’ amplitude and frequency [25, 30, 166]. We also
use the QNMs derived recently in ESGB gravity [138].

More precisely, for all harmonics except $(\ell,|m|)=(3,2)$ and $(4,3)$, which
exhibit post-merger oscillations due to mode-mixing [213, 214], the merger-
ringdown GR waveform employs the following ansatz [215, 216, 210, 156]:

$h_{\ell m}^{\textrm{merger-RD}}(\hat{t})=\nu\ \tilde{\mathcal{A}}_{\ell
m}(\hat{t})\ e^{i\tilde{\phi}_{\ell m}(\hat{t})}\ e^{-iM\sigma_{\ell
m0}(\hat{t}-\hat{t}_{\textrm{match}}^{\ell m})},$ (65)

where $\sigma_{\ell m0}$ is the complex frequency of the least-damped QNM
(with overtone number zero) of the remnant BH. The corresponding oscillation
frequency $f_{\ell m0}$ and damping time $\tau_{\ell m0}$ read

$\displaystyle f_{\ell m0}$
$\displaystyle=\frac{1}{2\pi}\operatorname{Re}\left(\sigma_{\ell m0}\right),$
(66a) $\displaystyle\tau_{\ell m0}$
$\displaystyle=-\frac{1}{\operatorname{Im}\left(\sigma_{\ell m0}\right)}.$
(66b)

The functions $\tilde{\mathcal{A}}_{\ell m}(\hat{t})$ and $\tilde{\phi}_{\ell
m}(\hat{t})$ are given by [156]:

$\displaystyle\tilde{\mathcal{A}}_{\ell m}(\hat{t})$
$\displaystyle=c_{1,c}^{\ell m}\tanh[c_{1,f}^{\ell m}\
(\hat{t}-\hat{t}_{\textrm{match}}^{\ell m})\ +\ c_{2,f}^{\ell m}]\ +\
c_{2,c}^{\ell m},$ (67a) $\displaystyle\tilde{\phi}_{\ell m}(\hat{t})$
$\displaystyle=\phi_{\textrm{match}}^{\ell m}-d_{1,c}^{\ell
m}\log\left[\frac{1+d_{2,f}^{\ell m}e^{-d_{1,f}^{\ell
m}(\hat{t}-\hat{t}_{\textrm{match}}^{\ell m})}}{1+d_{2,f}^{\ell m}}\right],$
(67b)

where $\phi_{\textrm{match}}^{\ell m}$ is the phase of the inspiral-plunge
mode $(\ell,m)$ at $\hat{t}=\hat{t}_{\textrm{match}}^{\ell m}$ [cf. below Eq.
(36)]. The coefficients $d_{1,c}^{\ell m}$ and $c_{i,c}^{\ell m}$ ($i=1,2$)
are constrained by the requirement that the amplitude and phase of $h_{\ell
m}(\hat{t})$ are continuously differentiable at
$\hat{t}=\hat{t}_{\text{match}}^{\ell m}$, while the coefficients
$c_{i,f}^{\ell m}$ and $d_{i,f}^{\ell m}$ are obtained through fits to a large
set of $\sim 440$ aligned-spin NR waveforms, spanning mass ratios up to 20 and
spins up to 0.998, and BH perturbation-theory merger-ringdown waveforms with
mass ratio 1000. For the $(3,2)$ and $(4,3)$ modes, the mode-mixing behavior
is modeled by applying the previous construction to the spheroidal harmonics
[217] $(3,2,0)$ and $(4,3,0)$, which feature a monotonic amplitude and
frequency evolution [218]. The $(3,2)$ and $(4,3)$ modes are then obtained by
combining the $(3,2,0)$ and $(4,3,0)$ modes with the $(2,2)$ and $(3,3)$ modes
using the appropriate mixing coefficients [219].

In the SEOBNRv5PHM model [156, 157], the complex QNM frequencies in GR are
obtained for each $(\ell,m)$ harmonic as a function of the BH’s
(dimensionless) final mass and spin $(\hat{M}_{f}=M_{f}/M,\,\chi_{f})$ using
the qnm Python package [220]. Note that the qnm package provides
$M_{f}\sigma_{\ell m0}$, as a function of the final spin $\chi_{f}$, so in Eq.
(65) we use

$M\sigma_{\ell m0}=\frac{1}{\hat{M}_{f}}M_{f}\sigma_{\ell m0}(\chi_{f}).$ (68)

The BH’s mass and spin are in turn computed using the fitting formulas of
Refs. [202] and [203], respectively.

We now describe how we account for ESGB effects in the merger-ringdown
$h_{\ell m}$ modes. For the merger, as reminded on Sec. V.1, the NQC
corrections impose that the modes’ amplitude and frequency agree with fits to
NR waveforms in GR. In recent years, there has been a significant effort to
produce NR simulations in beyond-GR theories [79, 82, 90, 95, 97, 98], with a
few examples of NR evolutions of BBHs in ESGB gravity [80, 85, 83, 94, 91,
99]. For example, Ref. [91] found a negligible effect on the GW amplitude at
merger with varying ESGB coupling, in contrast with the large effect found in
the order-reduced simulation of Ref. [83], possibly due to the presence of
secularly growing errors present in the latter perturbative approach (see
[103] for a detailed comparison of different numerical approaches). In any
case, making parameter-space fits (calibration) of such merger properties, as
done in GR, requires a large number of NR waveforms, which are not yet
available.

To account for possible deviations in the waveform at merger, we take
inspiration from the parametrized SEOBNR approach (pSEOBNR) [25, 30, 166], and
introduce fractional deviations to the NR predictions in GR for the amplitude
and frequency at $\hat{t}_{\text{match }}^{\ell m}$ for each mode. These
deviations parameterize our ignorance of the merger morphology with extra
nuisance parameters, to be marginalized over in parameter estimation. Given
this idea, we want the corrections to vanish in the GR limit, so we rescale
them as

$\displaystyle|h_{\ell m}^{\rm NR,GR}|$ $\displaystyle\rightarrow|h_{\ell
m}^{\rm NR,GR}|\,\left(1+{\hat{\ell}_{\mathrm{GB}}}^{4}\delta A_{\ell
m}\right),$ (69a) $\displaystyle\omega_{\ell m}^{\rm NR,GR}$
$\displaystyle\rightarrow\omega_{\ell m}^{\rm
NR,GR}\,\left(1+{\hat{\ell}_{\mathrm{GB}}}^{4}\delta\omega_{\ell m}\right),$
(69b)

where the ${\hat{\ell}_{\mathrm{GB}}}^{4}$ dependence is inspired by the
scaling of the PN corrections in the inspiral modes in the
small-${\hat{\ell}_{\mathrm{GB}}}$ regime. In turn, this rescaling changes the
values taken by the constants $(a_{i},b_{i})$ and $(c^{\ell m}_{i,c},d^{\ell
m}_{i,c})$ entering Eqs. (63) and (67), respectively. The $\hat{\ell}_{\rm
GB}$-dependence in Eqs. (69), and how the merger morphology is modified for
different modes, could be further studied using NR waveforms in ESGB gravity.
Still, current NR simulations [91, 99] suggest that corrections at merger are
small, for values of the coupling allowed by current observations.

In Sec. VI.3, we quantify the impact of marginalizing over these corrections,
applied for simplicity only to the dominant $(2,2)$ mode, on the analysis of
GW signals, and find it to be small.

###  V.3 Corrections to the quasi-normal-mode frequencies

To incorporate ESGB corrections in the merger-ringdown waveform using the
ansatz of Sec. V.2, one should account for the following effects:

  1. 1.

The complex frequency of the least-damped QNM for each $(\ell,m)$ harmonic
$\sigma_{\ell m0}$ should vary, given the BH’s final mass and spin.

  2. 2.

In several modified gravity theories, the BH metric perturbations with even
and odd parity have different QNM frequencies, breaking the isospectrality of
GR [221].

  3. 3.

The final mass and spin of the remnant BH should deviate from the GR
predictions due to differences in both the conservative and dissipative binary
dynamics. For example, the final mass and spin are smaller in ESGB gravity in
the NR simulations of Ref. [99].

  4. 4.

The $c_{i,f}^{\ell m}$ and $d_{i,f}^{\ell m}$ coefficients, modeling the early
ringdown, are obtained through fits of NR waveforms in GR, and should also be
corrected. Since these could only be obtained through fits of beyond-GR NR
waveforms, we neglect this last correction.

QNM frequencies for nonspinning ESGB BHs were first computed in Refs. [45,
50]. They were then extended up to the second order in a spin expansion in
Refs. [63, 66], focusing on the polar sector. Binary coalescences typically
lead to BHs with final spin $\chi_{f}\sim 0.7$, making the accuracy of an
expansion to second order in the spin limited. While Ref. [66] suggested that
its domain of validity could be extended using Padé approximants, Refs. [222,
223] argued that higher orders are necessary to achieve accuracy better than
the statistical uncertainties of current and future GW detectors. Various
groups developed modified Teukolsky formalisms [224, 225, 226, 227, 228] and a
novel approach based on perturbative spectral expansions [229, 230] to
calculate QNM frequencies in a wide class of modified gravity theories for BHs
with generic spins. Ref. [138] applied the latter formalism to ESGB gravity at
leading order in $\ell_{\rm GB}/M_{f}$, to compute corrections to the
$(\ell,m,n)=(2,2,0)$ QNM accurately for up to $\chi_{f}\sim 0.8$, which we
implement in this work.333We translate between the conventions of Ref. [138]
and ours by setting $\zeta^{1/4}=\ell_{\rm GB}/(4\pi^{1/4}M_{f})$.

Ref. [138] computed corrections for both the polar and axial sectors of the
QNMs. However, understanding how axial and polar QNMs appear in the GW signal
is a complex task outside the scope of this paper, to be answered on a theory-
by-theory basis. Following [39], we only apply corrections coming from the
polar sector, which for typical spin magnitudes $\chi_{f}\sim 0.7$ is the
least damped parity, as well as the one giving rise to larger frequency
corrections.

We note that Ref. [138] focuses on the shift-symmetric theory, while we
consider here the dilatonic case (5). However, when analyzing real events,
which also include the inspiral phase, we will obtain constraints that are
within the small coupling limit. In this case, the dilatonic theory reduces to
the shift-symmetric one at leading order in $\hat{\ell}_{\rm GB}$ [cf. Sec.
II.3], which justifies a-posteriori to use the results of Ref. [138].

###  V.4 Corrections to the final mass and spin

In this section, we discuss the ESGB corrections to the mass and spin of the
BH remnant as a function of the coupling, for a BBH in the dilatonic ESGB case
(5). While a precise estimate of these corrections relies on NR simulations,
in EOB models they can be approximated from the binary’s dynamics. For
simplicity, we consider spin-aligned binaries and only study the final mass
and spin magnitude. We also compare our results to a self-consistent estimate
following Refs. [33, 34].

Coefficient | $\Delta\hat{M}_{f}^{\mathrm{ESGB}}$ | $\Delta\hat{J}_{f}^{\mathrm{ESGB}}$
---|---|---
$a_{0}$ | $-4.1011\times 10^{-5}$ | $-5.1696\times 10^{-4}$
$a_{1}$ | $-21.4091$ | $-21.7572$
$a_{2}$ | $-1.5879$ | $-1.2942$
$a_{3}$ | $155.6233$ | $179.7129$
$a_{4}$ | $0.0340$ | $-0.0853$
$a_{5}$ | $27.9436$ | $21.2925$
$a_{6}$ | $7.9436$ | $5.2301$
$a_{7}$ | $-48.0261$ | $-30.6130$
$a_{8}$ | $-230.0593$ | $-537.8520$
Table 2: Coefficients for the fits of the final mass and final spin
corrections (71), for a BBH in the dilatonic theory
$f(\varphi)=\exp(2\varphi)/4$ and $\mathcal{A}(\varphi)=1$.

A simple phenomenological estimate for the final mass and spin of a BBH
coalescence can be made under the assumption that the system radiates most of
its energy and angular momentum in the inspiral stage until the ISCO, while
during the merger and ringdown stages the mass and angular momentum of the
system are roughly conserved [231]. A similar procedure was applied to ESGB
gravity in Refs. [232, 233]. Here, we estimate the (dimensionless) final mass
and orbital angular momentum from the values of
$\nu\hat{H}_{\mathrm{EOB}}=H_{\mathrm{EOB}}/M$ and
$\nu\hat{p}_{\phi}=p_{\phi}/M^{2}$ at the merger time, associated in the model
with the peak of the (2,2)-mode waveform (37). We compute the quantities above
in both ESGB and GR, and add the difference to the GR estimates given by the
NR fitting formulas of Refs. [202] and [203], used in SEOBNRv5PHM. This
amounts to neglecting the effects of ESGB corrections to the post-merger
radiation. Proceeding in an iterative manner, we evaluate the ESGB corrections
to the remnant mass and spin by first inserting their GR values in Eq. (37).
In a second step, in Sec. VI we build waveforms by evaluating $\hat{t}_{\rm
ISCO}$ in Eq. (37) using the corrected remnant mass and spin. Another
limitation of our estimates is that the merger time uses the
$\Delta\hat{t}_{\mathrm{NR}}$ parameter (37), which is calibrated from NR
simulations in GR, and should thus also be corrected in the future.

To avoid generating an additional GR waveform each time we wish to include
this correction, we provide fits as a function of the intrinsic parameters and
of the coupling $\hat{\ell}_{\mathrm{GB}}$. For simplicity, we explore
configurations with equal-spin magnitude $\chi\in[-0.9,0.9]$ with increment
$\Delta\chi=0.18$, for several mass ratios $q\in\\{1,2,4,10\\}$, and coupling
$\hat{\ell}_{\rm GB}^{4}\in[0,0.07]$ with increment $\Delta\hat{\ell}_{\rm
GB}^{4}=0.007$. We use a polynomial ansatz in $(\nu,\chi)$ (with $\chi$ being
promoted to $\chi_{+}=(m_{A}^{0}/M)\chi_{A}+(m_{B}^{0}/M)\chi_{B}$ for
unequal-spin cases below), while the dependence on $\hat{\ell}_{\mathrm{GB}}$
is well approximated by a simple $\hat{\ell}_{\mathrm{GB}}^{4}$ scaling. We
define:

$\displaystyle\hat{M}_{f}$
$\displaystyle=\hat{M}_{f}^{\mathrm{GR}}+\Delta\hat{M}_{f}^{\mathrm{ESGB}},$
(70a) $\displaystyle\hat{J}_{f}$
$\displaystyle=\hat{J}_{f}^{\mathrm{GR}}+\Delta\hat{J}_{f}^{\mathrm{ESGB}},$
(70b) $\displaystyle\chi_{f}$
$\displaystyle=\frac{\hat{J}_{f}}{\hat{M}_{f}^{2}}=\chi_{f}^{\mathrm{GR}}+\Delta\chi_{f}^{\mathrm{ESGB}},$
(70c)

where, recall, $\hat{M}_{f}=M_{f}/M$ and $\hat{J}_{f}=J_{f}/M^{2}$ are the
dimensionless final mass and total angular momentum.
$\hat{M}_{f}^{\mathrm{GR}}$, $\chi_{f}^{\mathrm{GR}}$ are given by the NR
fitting formulas of Refs. [202, 203]. We obtain the following fits for the
ESGB corrections:

$\displaystyle\Delta X^{\mathrm{ESGB}}=$
$\displaystyle\,a_{0}\hat{\ell}_{\mathrm{GB}}^{4}\left(1+a_{1}\nu+a_{2}\chi+a_{3}\nu^{2}+a_{4}\chi^{2}\right.$
(71)
$\displaystyle\left.+a_{5}\nu\chi+a_{6}\nu\chi^{2}+a_{7}\nu^{2}\chi+a_{8}\nu^{3}\right),$

where $X=(\hat{M}_{f},\hat{J}_{f})$. The respective coefficients are listed in
Table 2. By combining Eqs. (70), we can compute the linear correction to the
final spin as

$\Delta\chi_{f}^{\mathrm{ESGB}}=\frac{\Delta\hat{J}_{f}^{\mathrm{ESGB}}}{(\hat{M}_{f}^{\mathrm{GR}})^{2}}-2\chi_{f}^{\mathrm{GR}}\frac{\Delta\hat{M}_{f}^{\mathrm{ESGB}}}{\hat{M}_{f}^{\mathrm{GR}}}.$
(72)

Figure 4: ESGB corrections to the (dimensionless) final mass (top panel) and
angular momentum (bottom panel), as given by the fits (71). We consider BBHs
with equal-spin magnitude $\chi$ and mass ratio $q=4$, in the dilatonic theory
$f(\varphi)=\exp(2\varphi)/4$ and $\mathcal{A}(\varphi)=1$.

Figure 4 shows the values of the ESGB corrections to the (dimensionless) final
mass and angular momentum for a set of BBHs with equal-spin magnitude $\chi$
and mass ratio $q=4$, as a function of the ESGB coupling
$\hat{\ell}_{\mathrm{GB}}^{4}$ in the dilatonic theory (5). These corrections
are quite modest compared to those affecting the inspiral dynamics, being
$\mathcal{O}(10^{-4})$ for $\hat{\ell}_{\mathrm{GB}}^{4}\sim 1$. We find that
the corrections to the final mass are always negative, while the final angular
momentum corrections can take both signs depending on the mass ratio and
spins. Note that the value of $\hat{p}_{\phi}$ at the ISCO can be larger when
adding ESGB corrections to the EOB Hamiltonian [cf. Fig. 6 below].

One can also estimate the final mass and angular momentum by subtracting the
energy and angular momentum radiated during the entire coalescence (including
the merger-ringdown stage) from their initial values. These can in turn be
obtained by integrating the energy and angular momentum fluxes, in both the
metric and the scalar channels. The energy emitted through GWs depends on the
final mass and spin $(\hat{M}_{f},\chi_{f})$, as they enter the computation of
the merger-ringdown waveform modes used to compute the GW flux. One can then
numerically solve for self-consistent values of the final mass and spin [33,
34], which minimize the difference between the values set for the remnant BH
and those obtained by energy and angular-momentum balance. Using this method,
we find that the corrections show a less clear dependence on intrinsic
parameters, due to errors in the numerical integration and root-finding,
making a parameter-space fit harder to perform. However, the median values of
the corrections (across 10,000 configurations uniform in $\nu$ with
$q\in[1,20]$ and uniform in $\chi_{i}\in[-0.99,0.99]$) still show a clear
scaling with $\hat{\ell}_{\mathrm{GB}}^{4}$, with the average correction being
consistent with Eqs. (71).

Figure 5: Fractional change between ESGB and GR in the real and imaginary part
of the $(2,2)$-mode’s ringdown frequency, as a function of
$\hat{\ell}_{\mathrm{GB}}$, for a $q=4$ nonspinning BBH ($\chi_{f}\simeq
0.47$) in the dilatonic theory $f(\varphi)=\exp(2\varphi)/4$ and
$\mathcal{A}(\varphi)=1$. We show the change only including corrections to the
QNM spectrum, but fixing the final mass and spin to the GR values, only
changing the final mass and spin using Eq. (71), as well as considering both
corrections.

Figure 6: The orbital separation $\hat{r}$ (left panel) and angular momentum
$\hat{p}_{\phi}$ (right panel) as a function of time for a nonspinning BBH
with $q=4$, starting at the same dimensionless orbital frequency
$M\Omega_{i}=0.015$, for different values of the ESGB coupling
$\hat{\ell}_{\mathrm{GB}}\in\\{0,0.7,1.0\\}$ in the dilatonic theory
$f(\varphi)=\exp(2\varphi)/4$ and $\mathcal{A}(\varphi)=1$. The insets
highlight the differences in the dynamics at the start of the evolution, and
close to the ISCO, indicated by the dots. Each curve ends at merger, taken as
the peak of the (2,2) mode. We note that: (i) The binary merges at earlier
times for larger values of the ESGB coupling, since more energy is radiated
away through the additional scalar field, (ii) ESGB binaries at the same
initial dimensionless frequency start at a larger radius and angular momentum.
This is due to the leading-order correction $1/\hat{r}\rightarrow
G_{AB}/\hat{r}$ in the EOB potentials, with $G_{AB}\geq 1$ (iii) The ISCO is
at a smaller separation and slightly larger angular momentum compared to GR.
This is a non-trivial effect due to higher-order PN corrections in the EOB
Hamiltonian.

In Fig. 5, we show the fractional change between ESGB and GR in the real and
imaginary part of the $(2,2)$ mode ringdown frequency, as a function of
$\hat{\ell}_{\mathrm{GB}}$, for a $q=4$ nonspinning BBH ($\chi_{f}\simeq
0.47$). We either include corrections to the QNM spectrum, change the final
mass and spin using fits from Eq. (71), or consider both corrections. As noted
in previous works [66, 40], typical corrections to the QNMs are very small.
While corrections to the remnant mass and spin are also modest, they can cause
changes in the frequency of the ringdown waveform from one to two orders of
magnitude larger than the ESGB corrections to the QNM spectrum. The precise
estimation of such corrections would be useful for high-precision tests of GR
with next-generation detectors, since prospects for BH spectroscopy in ESGB
could be much more optimistic than originally thought. Reference [40] found
that, even with next-generation GW detectors, constraints on ESGB gravity from
ringdown observations alone are weaker than current bounds from LVK inspiral
observations, when only accounting for corrections to the QNM spectrum. Our
results suggest that much better constraints might be achieved by exploiting
complementary information from the full IMR signal.

##  VI Results

Given the coupling of matter fields to the metric in ESGB gravity (1), the
mirrors of a GW interferometer must follow the geodesics of the perturbed
Jordan metric $\tilde{g}_{\mu\nu}=\mathcal{A}^{2}(\varphi)g_{\mu\nu}$ which
reads, in the Solar System:

$\displaystyle\tilde{g}_{\mu\nu}=\mathcal{A}^{2}_{\odot}\left[\eta_{\mu\nu}(1+2\alpha_{\odot}\delta\varphi)-h_{\mu\nu}^{\rm
TT}\right]+\mathcal{O}\left(\frac{1}{R^{2}}\right)\,,$ (73)

where $R$ is the distance to the source, while
$\mathcal{A}_{\odot}=\mathcal{A}(\varphi_{\odot})$ and
$\alpha_{\odot}=(d\ln\mathcal{A}/d\varphi)(\varphi_{\odot})$ are the couplings
of matter to the scalar field in the Solar System. Note that for the Jordan
metric to reduce to Minkowski in the absence of metric and scalar waves, we
must perform the local coordinate change $d\tilde{t}=\mathcal{A}_{\odot}dt$
and $d\tilde{x}^{i}=\mathcal{A}_{\odot}dx^{i}$.

In what follows, we focus on the metric waveform $h_{\mu\nu}^{\rm TT}$, which
was discussed in the previous sections. Indeed, the contribution of the scalar
waveform $\delta\varphi=\varphi-\varphi_{\odot}$ is cut off by the numerical
factor $\alpha_{\odot}^{2}\lesssim 10^{-5}$ in the Solar System [234].
Moreover, we illustrate our results in the dilatonic ESGB case (5) for which
$\mathcal{A}(\varphi)=1$, and thus $\mathcal{A}_{\odot}=1$ and
$\alpha_{\odot}=0$.

###  VI.1 Waveform morphology

In this section, we describe the impact of the ESGB corrections on the
dynamics and waveforms of BBH systems, modeled with SEOBNRv5PHM$\\_$ESGB. As a
concrete illustration, we consider a nonspinning binary with mass ratio of
$q=4$, and $\hat{\ell}_{\mathrm{GB}}\in\\{0,0.7,1.0\\}$. Our conclusions
remain qualitatively unchanged for different mass ratios and non-zero spins.
However, in the equal-mass limit, the deviations from the GR flux are smaller,
because the -1PN leading-order term in the scalar flux (IV.2) then vanishes.

First, we examine in Fig. 6 the orbital separation $\hat{r}$ and angular
momentum $\hat{p}_{\phi}$ as a function of time. The insets highlight the
differences in the dynamics at the start of the evolution, which for all cases
is taken at a dimensionless orbital frequency $M\Omega_{i}=0.015$, and close
to the ISCO, indicated by the dots. The values of $\hat{r}$ and
$\hat{p}_{\phi}$ at the ISCO are determined by setting both the first and
second derivatives of the Hamiltonian with respect to $\hat{r}$ to zero,
taking $\hat{p}_{r_{*}}=0$ (circular orbits) (i.e., $\partial
H_{\mathrm{EOB}}(\hat{p}_{r_{*}}=0)/\partial\hat{r}=0=\partial^{2}H_{\mathrm{EOB}}(\hat{p}_{r_{*}}=0)/\partial\hat{r}^{2}$).
Each curve ends at merger, taken as the peak of the (2,2) mode waveform (37).
Observe that, on the examples considered:

  1. 1.

The binary merges at earlier times for larger values of the ESGB coupling,
since more energy is radiated away due to the additional scalar field.

  2. 2.

Given a fixed initial dimensionless frequency, binaries start at a larger
radius and angular momentum for larger values of the ESGB coupling. This is
due to the replacement $1/\hat{r}\to G_{AB}/\hat{r}$ in the EOB potentials,
with $G_{AB}\geq 1$. Note that the model uses the quasi-circular, adiabatic,
initial condition formulas derived in Ref. [154], also used in the SEOBNRv5PHM
model.

  3. 3.

The ISCO is at a smaller separation [cf. also Ref. [68]] and slightly larger
angular momentum compared to GR. This is a non-trivial effect due to higher-
order PN corrections in the EOB Hamiltonian.

In Fig. 7 we show, for the same set of binaries, the amplitude of the (2,2)
mode waveform, to highlight how a non-zero value of the coupling changes the
time-to-merger when including different beyond-GR corrections. The dash-dotted
vertical lines indicate the merger, taken as the peak of the (2,2) mode, and
the dotted vertical lines correspond to the time at which
$\hat{r}=\hat{r}_{\rm{ISCO}}$. We set $\hat{t}=0$ at the GR merger time. In
the top panel, we incorporate ESGB corrections in the Hamiltonian, while using
the GR flux, and mostly observe that modifications to the conservative
dynamics tend to delay the merger. In the middle panel, we include instead
ESGB corrections in the dissipative sector and use the GR Hamiltonian. As
expected, these corrections accelerate the inspiral, mostly due to the
additional energy dissipation via the scalar field. In the bottom panel, we
combine both contributions. Changes to the dissipative dynamics are
predominant, and the binary merges earlier than GR for a non-zero value of the
ESGB coupling. The trend seen in this example is consistent with other cases
we considered, and agrees qualitatively with NR simulations [91, 99]. Finally,
Fig. 8 shows the real part of the (2,2) mode for the same set of binaries. The
waveforms are aligned over a low-frequency interval, by shifting their time of
coalescence $t_{c}$ and phase $\Phi$ [cf. Sec. IV] as done in Ref. [206], and
display a difference in the time-to-merger and a noticeable dephasing as they
approach the merger.

Note that the ESGB corrections to the remnant mass and spin [entering the
merger time and ringdown frequency, cf. Sec. V.4] and to the QNM spectrum [cf.
Sec. V.3] are included in our waveforms. However, their impact is subdominant
compared to the effects above.

Figure 7: Amplitude of the (2,2) mode of a nonspinning BBH with $q=4$,
starting at a dimensionless frequency $M\Omega_{i}=0.015$, for
$\hat{\ell}_{\mathrm{GB}}\in\\{0,0.7,1.0\\}$ in the dilatonic theory
$f(\varphi)=\exp(2\varphi)/4$ and $\mathcal{A}(\varphi)=1$. The dash-dotted
vertical lines indicate the merger, taken as the peak of the (2,2) mode, and
the dotted vertical lines correspond to the time at which
$\hat{r}=\hat{r}_{\rm{ISCO}}$. We set $\hat{t}=0$ at the GR merger time. Top
panel: we include ESGB corrections in the Hamiltonian, and use the GR flux.
Middle panel: we include ESGB corrections in the flux, and use the GR
Hamiltonian. Bottom panel: we include ESGB corrections in both the Hamiltonian
and the flux. We see that corrections to the conservative dynamics tend to
delay the merger, while corrections to the dissipative sector accelerate the
inspiral. The trend seen here is consistent with other cases we considered:
when combining all contributions, the changes to the dissipative dynamics are
predominant, and the binary merges earlier than in GR for a non-zero value of
the coupling.

###  VI.2 Quantifying measurability of deviations from GR

We now quantify which deviation magnitudes are detectable by GW detectors,
and, if so, estimate the SNR required for such measurements. For this purpose,
we employ a common metric in GW data analysis - the mismatch. Given two
waveforms $h_{1}(t)$ and $h_{2}(t)$, we introduce the overlap, defined as the
noise-weighted inner product [235, 236]

$\left(h_{1}\mid h_{2}\right)\equiv
4\operatorname{Re}\int_{f_{l}}^{f_{h}}\frac{\tilde{h}_{1}(f)\tilde{h}_{2}^{*}(f)}{S_{n}(f)}\mathrm{d}f,$
(74)

where $\tilde{h}(f)$ is the Fourier transform of the time-domain signal, and
$S_{n}(f)$ is the one-sided power spectral density (PSD) of the detector
noise. The faithfulness is then defined as the overlap between the normalized
waveforms, maximized over the relative time of coalescence and phase, that is

$\left\langle h_{1}\mid
h_{2}\right\rangle=\max_{\Phi,t_{c}}\frac{\left(h_{1}\left(\Phi,t_{c}\right)\mid
h_{2}\right)}{\sqrt{\left(h_{1}\mid h_{1}\right)\left(h_{2}\mid
h_{2}\right)}}.$ (75)

We finally define the mismatch, or unfaithfulness, as

$\mathcal{M}=1-\left\langle h_{1}\mid h_{2}\right\rangle.$ (76)

In Eq. (74), we fix $f_{h}=2048\leavevmode\nobreak\ \text{Hz}$ and choose
different values for $f_{l}$ depending on the detector. We consider the
advanced LIGO (aLIGO) detector, for which we assume the design zero-detuned
high-power noise PSD [237] and $f_{l}=20\mathrm{Hz}$, and the ET detector, for
which we assume the EinsteinTelescopeP1600143 PSD implemented in pyCBC [238]
and $f_{l}=5\mathrm{Hz}$. For simplicity, we only consider the dominant (2,2)
mode of the waveform when computing the mismatch.

Figure 8: (2,2) mode of the SEOBNRv5PHM$\\_$ESGB waveform of a nonspinning BBH
with $q=4$, starting at a dimensionless orbital frequency $\Omega_{i}=0.015$,
for different values of the ESGB coupling
$\hat{\ell}_{\mathrm{GB}}\in\\{0,0.7,1.0\\}$ in the dilatonic theory
$f(\varphi)=\exp(2\varphi)/4$ and $\mathcal{A}(\varphi)=1$. We set $\hat{t}=0$
at the peak of the (2,2) mode of the GR waveform. In this example, the merger
parameters $\delta A_{\ell m}$ and $\delta\omega_{\ell m}$ entering Eqs. (69)
are set to zero.

The mismatch can be used to estimate at which SNR two signals could be
distinguishable, given the detector sensitivity. This is usually formulated in
terms of an indistinguishability criterion [239, 240, 241, 242], which states
that if two waveforms fulfill the condition

$\mathcal{M}<\frac{D}{2\leavevmode\nobreak\ \mathrm{SNR}^{2}},$ (77)

for a given PSD and SNR, then these waveforms are considered
indistinguishable, and the differences in the recovered parameters are
expected to be smaller than statistical errors. The prefactor $D$ can be
estimated as the number of intrinsic parameters whose measure is affected by
model differences [241], or it can be tuned by considering synthetic
injections at increasing SNR [242]. In Fig. 9 we show the mismatch of
SEOBNRv5PHM$\\_$ESGB against its GR limit as a function of
$\hat{\ell}_{\mathrm{GB}}$, for the aLIGO and ET detectors. We consider a
nonspinning binary with total mass $20M_{\odot}$ and mass ratios
$q\in\\{2,4,8\\}$ as a signal. The horizontal lines correspond to the
indistinguishability thresholds obtained by taking $D=9$ (our templates depend
on two masses, $6$ spin components and $\hat{\ell}_{\mathrm{GB}}$) for some
realistic SNRs ($10,30$ for aLIGO, $100,300$ for ET). As expected, we find a
larger mismatch for more asymmetric binaries, due to the scalar dipole
radiation. This results in a stronger bound on $\hat{\ell}_{\mathrm{GB}}$ if
the data is consistent with GR. The minimum value of
$\hat{\ell}_{\mathrm{GB}}$ that could be distinguished from the GR case goes
from around $0.8$ ($q=2$, $\rm{SNR}=10$) to $0.3$ ($q=8$, $\rm{SNR}=30$) for
aLIGO and from $0.3$ ($q=2$, $\rm{SNR}=100$) to $0.1$ ($q=8$, $\rm{SNR}=300$)
for ET. The main source of improvement for ET is its larger low-frequency
bandwidth, where the scalar flux dominates (see Fig. 3). Mismatch values
$\lesssim 10^{-5}$ are not to be attributed to a non-zero coupling value, as
differences of this order can be simply due to errors in the numerical
integration of Hamilton’s equations with finite tolerance. Comparing our
results to current constraints on ESGB gravity with LVK observations, the
bound on GW190814 from Refs. [168, 170] (assuming the event is from a BBH), is
given as a function of the coupling
$\sqrt{\alpha_{\rm{GB}}}=\ell_{\rm{GB}}/(4\pi^{1/4})$ with $\ell_{\rm
GB}/M=\nu\hat{\ell}_{\rm{GB}}$ by $\sqrt{\alpha_{\rm{GB}}}\lesssim
0.4\mathrm{km}$. Taking a total mass of $\sim 26M_{\odot}$ and a mass ratio
$q=10$, this corresponds to $\hat{\ell}_{\rm{GB}}\lesssim 0.67$. Our
prospective bound $\hat{\ell}_{\rm{GB}}\lesssim 0.1$ for an ET source with
$q=8$, $M=20M_{\odot}$ translates instead to $\sqrt{\alpha_{\rm{GB}}}\lesssim
0.055\mathrm{km}$.

Figure 9: Mismatch of the (2,2) mode waveform of SEOBNRv5PHM$\\_$ESGB against
its GR limit, as a function of the ESGB coupling $\hat{\ell}_{\mathrm{GB}}$
for the aLIGO (left panel) and ET (right panel) detectors. We consider a
nonspinning BBH with total mass $20M_{\odot}$ and mass ratios
$q\in\\{2,4,8\\}$ in the dilatonic theory $f(\varphi)=\exp(2\varphi)/4$ and
$\mathcal{A}(\varphi)=1$. The horizontal lines show the indistinguishability
thresholds given by Eq. (77) for different SNRs.

The potential constraints on $\hat{\ell}_{\mathrm{GB}}$ should only be
regarded as an order of magnitude estimate. Since the indistinguishability
criterion offers a necessary, but not sufficient measure, if the mismatch
exceeds the threshold, the differences are not necessarily measurable [243,
244]. Furthermore, only varying $\hat{\ell}_{\mathrm{GB}}$ does not account
for possible correlations between the ESGB coupling and GR parameters. For
example, the main effect of a non-zero coupling being to shorten the waveform,
one could expect it to be anti-correlated with the chirp mass
$\mathcal{M}_{c}=\mu^{3/5}M^{2/5}$ of the binary [245],which affects the GW
phase at leading PN order. Spin effects also affect the waveform in a
qualitatively similar fashion, although they enter the GW phase starting from
1.5 PN order [107], and are thus expected to be less degenerate with ESGB
corrections (starting from -1 PN) than the chirp mass. Therefore, this
provides an optimistic estimate of the constraint derived from a single event.
On the other hand, these bounds could be significantly improved by combining
results from multiple observations. Since the coupling is common to all
analyzed events, one can stack constraints from multiple observations by
simply multiplying the marginalized likelihoods [173, 169], instead of
combining the results using a hierarchical approach, assuming that beyond-GR
parameters in each event are drawn from a common underlying distribution
[246], as typically done for parametrized deviations.

More sophisticated predictions of the measurability of the ESGB coupling could
be either based on the Fisher matrix formalism or employ full Bayesian
inference. Forecasts of constraints with future detectors on ESGB gravity, as
well as on other beyond-GR theories, have been obtained in Ref. [18]. This
work performed a Fisher analysis including corrections to the inspiral
waveform based on the parameterized post-Einsteinian (ppE) framework [21],
mapping then the ppE results to theory-specific constraints, also including
the effect of stacking multiple observations. A similar study, focused instead
on the ringdown stage, was recently presented in Ref. [40], using the
parametrized spin-expansion coefficients (ParSpec) framework [28].

Similar analyses could be refined by repeating them with SEOBNRv5PHM$\\_$ESGB
. Ref. [18] employed as baseline GR waveform IMRPhenomPv2 [147], which is a
less accurate model than SEOBNRv5PHM. It does not include higher modes and
employs a simplified modeling of spin precession. It also focused on
estimating modifications at a single PN order at a time, while in our model,
modifications at each PN order are specified as a function of the coupling
parameter of the theory. Furthermore, we include beyond-GR corrections in both
the inspiral and the ringdown stages. One could employ Bayesian parameter
estimation to provide forecasts for the next-generation GW detectors ET and
CE, examining correlations among different parameters. We leave such analyses
for future work. However, to show a practical application of
SEOBNRv5PHM$\\_$ESGB to GW data analysis, we perform full Bayesian parameter
estimation on some of the GW events reported by the LVK collaboration.

###  VI.3 Constraining the ESGB coupling

As an application of SEOBNRv5PHM$\\_$ESGB to GW data analysis, we perform full
Bayesian parameter estimation for some GW events reported by the LVK
collaboration, also including the ESGB coupling constant $\ell_{\mathrm{GB}}$
[or equivalently $\sqrt{\alpha_{\mathrm{GB}}}$, cf. Footnote 2] in the
analysis. We consider three signals, GW190814 [247], GW190412 [248] and
GW230529 [249], which were recently analyzed in Refs. [169, 168, 170, 171,
172] to provide constraints on the coupling of the dilatonic ESGB theory (5),
using waveform models featuring ESGB corrections in the inspiral phase.

We expect larger corrections to the waveform for systems with smaller masses,
as the leading ESGB corrections to the phase are proportional to
$\hat{\ell}_{\rm GB}^{4}$ [61, 68, 170, 171], as well as for large mass
ratios, due to the increased scalar dipole radiation [cf. Sec IV.2]. As the
sources of GW190412 and GW190814 are relatively low-mass binaries with
asymmetric masses, these events are expected to offer better constraints on
this theory. The source of GW190814, having a secondary with a mass around
$2.6M_{\odot}$, is consistent with being either a BBH or a NSBH. Under the
hypothesis of GW190814 being a BBH, it imposes stringent bounds on ESGB
gravity around $\sqrt{\alpha_{\mathrm{GB}}}\lesssim 0.4\leavevmode\nobreak\
\mathrm{km}$ [168, 170]. However, this bound is not necessarily robust as it
changes significantly whether one assumes the event to be an NSBH or a BBH.
Indeed, the scalar monopole vanishes for a NS in shift-symmetric ESGB gravity,
so the correction can be much larger for a BBH than for a NSBH, whose hair is
carried by the BH (which is more massive) [169, 170].

The most probable interpretation of the source of GW230529 is that of a NSBH
with the primary being a BH with mass $2.5\text{--}4.5\leavevmode\nobreak\
M_{\odot}$. Under this hypothesis, it currently places the best constraint on
ESGB gravity [171, 172], with a $90\%$ bound on $\sqrt{\alpha_{\mathrm{GB}}}$
ranging from $0.26\leavevmode\nobreak\ \mathrm{km}$ to
$0.35\leavevmode\nobreak\ \mathrm{km}$ depending on the waveform approximant
employed, when including ESGB corrections in the inspiral phase. We note that
for GW230529 one cannot determine from GW data alone if either component of
the source is a NS or a BH, and the primary analysis from the LVK
collaboration [249] employed BBH waveforms, including the SEOBNRv5PHM model
used here as a baseline GR approximant. Furthermore, Ref. [171] finds that,
when considering parameterized deviations to the inspiral phase, NSBH and BNS
results with the tides constrained to realistic values are very similar to the
BBH results. Therefore, in our analysis, we also neglect tidal effects, but we
take the sensitivities for the secondary to be 0, as is the case for NSs in
shift-symmetric ESGB theories and in the dilatonic ESGB theory at leading
order in the coupling [cf. Sec II.3]. In this work we concentrate on GW190814
and GW190412 under the BBH assumption, and on GW230529 under the NSBH
assumption, for simplicity.

Figure 10: 1-dimensional marginalized posterior of $\ell_{\rm GB}$ ($\sqrt{\alpha_{\mathrm{GB}}}$) obtained from the analysis of GW230529, GW190814 and GW190412, using the SEOBNRv5PHM$\\_$ESGB model in the dilatonic case $f(\varphi)=\exp(2\varphi)/4$ and $\mathcal{A}(\varphi)=1$. The vertical dashed lines indicate the $90\%$ confidence bound. We show the posteriors for the full IMR analyses, without adding the merger corrections (69) (light blue) or marginalizing over them (yellow), and for inspiral-only analyses (dark blue). We perform a single analysis for GW230529, since it has negligible SNR in the merger-ringdown stage. Event |  GW230529 |  GW190814 |  GW190412 |  Combined
---|---|---|---|---
IMR |  $1.66\leavevmode\nobreak\ \text{km}$ ($0.31\leavevmode\nobreak\ \text{km}$) |  $2.66\leavevmode\nobreak\ \text{km}$ ($0.50\leavevmode\nobreak\ \text{km}$) |  $9.65\leavevmode\nobreak\ \text{km}$ ($1.81\leavevmode\nobreak\ \text{km}$) |  $1.72\leavevmode\nobreak\ \text{km}$ ($0.32\leavevmode\nobreak\ \text{km}$)
IMR with merger deviations | – |  $2.55\leavevmode\nobreak\ \text{km}$ ($0.48\leavevmode\nobreak\ \text{km}$) |  $9.42\leavevmode\nobreak\ \text{km}$ ($1.77\leavevmode\nobreak\ \text{km}$) |  $1.65\leavevmode\nobreak\ \text{km}$ ($0.31\leavevmode\nobreak\ \text{km}$)
Inspiral only | – |  $2.19\leavevmode\nobreak\ \text{km}$ ($0.41\leavevmode\nobreak\ \text{km}$) |  $11.42\leavevmode\nobreak\ \text{km}$ ($2.15\leavevmode\nobreak\ \text{km}$) |  $1.60\leavevmode\nobreak\ \text{km}$ ($0.30\leavevmode\nobreak\ \text{km}$)
Table 3: Summary of the $90\%$ confidence bound on $\ell_{\rm GB}$
($\sqrt{\alpha_{\mathrm{GB}}}$) from the analysis of GW230529, GW190814 and
GW190412, using the SEOBNRv5PHM$\\_$ESGB model in the dilatonic case
$f(\varphi)=\exp(2\varphi)/4$ and $\mathcal{A}(\varphi)=1$. We list the bounds
from the full IMR analyses, without adding deviations to the merger and
marginalized over parameterized merger deviations, and for inspiral-only
analyses. We perform a single analysis for GW230529, since it has negligible
SNR in the merger-ringdown stage.

The goal of Bayesian parameter estimation is to infer the posterior
distribution $p(\theta|d)$ for the parameters $\theta$ given the observed data
$d$, using Bayes theorem

$p(\theta\mid d)=\frac{\mathcal{L}(d\mid\theta)\pi(\theta)}{Z}$ (78)

where $\mathcal{L}(d|\theta)$ is the likelihood of the data $d$ given the
parameters $\theta$, $\pi(\theta)$ is the prior on $\theta$, and
$Z\equiv\int\mathrm{d}\theta\mathcal{L}(d\mid\theta)\pi(\theta)$ is the
evidence. To determine whether a model $A$ is preferred over a model $B$, one
can look at the Bayes factor, defined as the ratio of the evidence for the two
different models $\mathrm{BF}^{A}_{B}=Z_{A}/Z_{B}$ (see, e.g., Ref. [250] for
a review).

We perform the analysis using parallel Bilby [251], a highly parallelized
version of the Bilby [252, 253] parameter-estimation code, and the nested
sampler dynesty [254] using the acceptance-walk stepping method. We use the
recommended setting for the number of accepted MCMC-chains
$\mathrm{naccept}=60$, number of live points $\mathrm{nlive}=1000$, and set
the remaining sampling parameters to their default values. For GW190814, we
find that these settings were not sufficient to fully resolve the
$\ell_{\mathrm{GB}}$ posterior distribution, and use $\mathrm{naccept}=100$,
and $\mathrm{nlive}=2000$. We employ strain data from the Gravitational Wave
Open Source Catalog (GWOSC) [255] and the released PSD and calibration
envelopes. We use the same priors for the GR parameters as those employed in
Refs. [157, 249], in which these events have been analyzed using the
SEOBNRv5PHM waveform model, and uniform priors on $\ell_{\rm GB}$, with the
constraint $\ell_{\mathrm{GB}}/m^{0}_{B}\lesssim 0.831$ in the BBH case [cf.
Sec. II.3]. In the NSBH case, we rather impose
$\ell_{\mathrm{GB}}/m^{0}_{A}\lesssim 0.831$. To facilitate a more direct
comparison with other results in the literature, we also perform inspiral-only
analyses, taking a maximum frequency for the likelihood of
$137\leavevmode\nobreak\ \mathrm{Hz}$ for GW190814 and $83\leavevmode\nobreak\
\mathrm{Hz}$ for GW190412, following Ref. [168]. These frequencies coincide
with the ISCO frequency in GR, and thus, doing so amounts to cutting off the
late stages of our inspiral-plunge waveforms above. We don’t perform an
inspiral-only analysis for GW230529 since the cutoff frequency for the
inspiral stage would be around $600\leavevmode\nobreak\ \mathrm{Hz}$, and
while the default analysis extends up to $1792\leavevmode\nobreak\
\mathrm{Hz}$, we don’t expect significant differences due to the reduced
sensitivity of the LIGO detectors in this high-frequency range.

As detailed in Sec. V.2, our model allows for parametrized deviations in the
merger waveform amplitude and frequency (69), which can be varied during
sampling. We perform runs assuming that they take GR values, as well as
allowing them to vary for the $(2,2)$ mode. We use uniform priors in the range
$[-1,1]$ for both $\delta A_{22}$ and $\delta\omega_{22}$. We introduce
additional merger parameters only in the IMR analysis, as the merger
deviations do not significantly impact the results when only analyzing the
inspiral. For the same reason, we don’t include these corrections in the
analysis of GW230529, which has negligible SNR in the merger-ringdown stage.

Figure 11: Posterior distribution of $\ell_{\mathrm{GB}}$, inverse mass ratio
$1/q$ and chirp mass $\mathcal{M}_{c}$ obtained from the analysis of GW190814
using the SEOBNRv5PHM$\\_$ESGB model in the dilatonic case
$f(\varphi)=\exp(2\varphi)/4$ and $\mathcal{A}(\varphi)=1$. For the masses, we
also show the values recovered assuming GR using the SEOBNRv5PHM model. We
notice a slight anti-correlation between $\ell_{\mathrm{GB}}$ and the chirp
mass, as well as a wider chirp mass posterior shifted to larger values for the
inspiral-only analysis. This might contribute to a marginally better
constraint on $\ell_{\mathrm{GB}}$ for the inspiral-only analysis, despite the
lower recovered SNR.

Figure 10 presents the 1-dimensional posteriors of $\ell_{\mathrm{GB}}$ (or,
equivalently, $\sqrt{\alpha_{\mathrm{GB}}}$), marginalized over all other
parameters, resulting from the analysis of GW230529, GW190814 and GW190412.
The vertical dashed lines indicate the $90\%$ confidence bounds, also
summarized in Table 3. We show the posteriors for the full IMR analyses,
without adding corrections to the merger amplitude and frequency (light blue)
and marginalized over parameterized merger deviations (yellow), and for
inspiral-only analyses (dark blue). In all cases, the posteriors are
compatible with $\ell_{\rm GB}=0$.

For the NSBH event GW230529, we find a $90\%$ bound on
$\sqrt{\alpha_{\mathrm{GB}}}$ of $0.31\leavevmode\nobreak\ \mathrm{km}$. Our
bound is slightly less stringent than those of Refs. [171, 172] when relying
on a BBH baseline model in GR, as they found
$\sqrt{\alpha_{\mathrm{GB}}}\lesssim 0.28\leavevmode\nobreak\ \mathrm{km}$ and
$\sqrt{\alpha_{\mathrm{GB}}}\lesssim 0.26\leavevmode\nobreak\ \mathrm{km}$,
respectively. However, our analysis differs in three ways: (i) we include ESGB
corrections at higher PN orders. The references above indeed include
corrections to the conservative sector at 2PN, and to the dissipative sector
at 1.5PN and 1PN relatively to the quadrupolar order, respectively; (ii) our
model resums these ESGB corrections within the EOB framework; (iii) the
waveform models of these works are not the same as ours in the GR limit, i.e.
SEOBNRv4HM$\\_$ROM [256] and IMRPhenomXPHM [148] respectively. This last point
alone can induce comparable bound discrepancies between these works.

For GW190814, our results are also in agreement with those of Refs. [168,
170], with a $90\%$ bound on $\sqrt{\alpha_{\mathrm{GB}}}$ ranging from
$0.41\leavevmode\nobreak\ \mathrm{km}$ to $0.48\leavevmode\nobreak\
\mathrm{km}$. We notice that, for this event, the inspiral-only analysis
provides a slightly better constraint compared to the IMR analyses. Figure 11
shows the posterior distribution of $\ell_{\mathrm{GB}}$, inverse mass ratio
$1/q$ and chirp mass $\mathcal{M}_{c}$ inferred from GW190814, along with the
mass values recovered assuming GR. We observe: (i) a slight anti-correlation
between $\ell_{\mathrm{GB}}$ and the chirp mass, whether the plunge-merger-
ringdown is included or not, as expected for inspiral-dominated signals [cf.
Sec. VI.2]; and (ii) that neglecting the plunge-merger-ringdown signal results
in a wider chirp mass posterior shifted towards larger values. These two
points might contribute to the slightly better constraint from the inspiral-
only analysis. Another factor could be that the merger-ringdown waveform is
practically indistinguishable from GR for this event. Therefore, its inclusion
reduces the mismatch of the waveform with its GR counterpart at fixed
$\ell_{\rm GB}$.

For the BBH event GW190412, we obtain a constraint
$\sqrt{\alpha_{\mathrm{GB}}}=2.15\leavevmode\nobreak\ \mathrm{km}$ from the
inspiral-only analysis and a slightly tighter bound
$\sqrt{\alpha_{\mathrm{GB}}}=1.77\leavevmode\nobreak\ \mathrm{km}$ when
including the plunge-merger-ringdown portions of the signal, thanks to the
additional SNR. This is roughly a factor $\sim 2$ smaller compared to previous
analyses [168]. The larger chirp mass of the source ($\mathcal{M}_{c}\sim
15\leavevmode\nobreak\ M_{\odot}$, versus $\mathcal{M}_{c}\sim
6\leavevmode\nobreak\ M_{\odot}$ for GW190814 and $\mathcal{M}_{c}\sim
2\leavevmode\nobreak\ M_{\odot}$ for GW230529) results in a higher relative
SNR in the late-inspiral stage, increasing the impact of higher order PN
corrections and the EOB resummation. For GW190412 and GW190814, the bounds on
$\sqrt{\alpha_{\mathrm{GB}}}$ are robust with respect to uncertainties in the
merger morphology, with a slight decreas when marginalized over parameterized
merger deviations.

Combining the constraints from individual events by multiplying the
marginalized posteriors, we obtain a combined bound on
$\sqrt{\alpha_{\mathrm{GB}}}$ ranging from $0.30\leavevmode\nobreak\
\mathrm{km}$ (for the inspiral-only analyses) to
$0.31-0.32\leavevmode\nobreak\ \mathrm{km}$ (for the IMR analyses, with and
without merger deviations respectively). The combined result is very similar
to the one obtained by analyzing GW230529 alone, and if anything slightly
worsened in the IMR analyses by the fact that the
$\sqrt{\alpha_{\mathrm{GB}}}$ posterior for GW190814 peaks slightly away from
0.

We find, a posteriori, that our bounds are within the small-coupling limit.
This justifies using the QNM expressions from Ref. [138], as well as employing
the leading-order expressions (6) and (7) for the quantities in Eqs. (4) when
analyzing GW230529 [cf. Sec. II.3]. For GW230529, we find
$\ell_{\mathrm{GB}}/m_{A}^{0}\lesssim 0.28$ at $90\%$ confidence. For GW190412
we obtain $\ell_{\mathrm{GB}}/M_{f}\lesssim 0.15$, while for GW190814
$\ell_{\mathrm{GB}}/M_{f}\lesssim 0.065$.

Table 4 summarizes the natural log-Bayes factors between the ESGB and GR
hypotheses, inferred from the full IMR analyses of the signals. Without adding
deviations to the merger amplitude and frequency, we find a log Bayes factor
of $-0.94$ for GW230529, of $-0.02$ for GW190814 and of $-0.52$ GW190412,
indicating no substantial preference for either hypothesis [257]. When adding
parametrized merger deviations, we obtain lower Bayes factors, of $-3.0$ for
GW190814 and of $-4.1$ GW190412. This means that they do not improve how well
our model fits the data, which is consistent with GW190412 and GW190814 being
inspiral-dominated signals with low SNR in the merger-ringdown. The decrease
in the Bayes factor is then the manifestation of the Occam penalty for
introducing additional parameters. Accurate modeling of the merger phase with
NR simulations in ESGB gravity would be valuable to improve this analysis.

We note that the computational cost of these runs is similar to the
corresponding GR analyses with SEOBNRv5PHM, slightly higher due to additional
parameters and a marginally longer waveform evaluation time. They took around
36 hours for GW190412 (256 cores), 3 days for GW190814 (384 cores) and 5 days
for GW230529 (640 cores), with each analysis requiring $\mathcal{O}(10^{8})$
waveform evaluations.

Our combined result could be further improved by also including the NSBH
events GW200105, GW200115, which also place competitive bounds on this theory
[170], as well as selected BBH events. Since the NSBHs are low-mass binaries,
leading to inspiral-dominated signals, the constraints on ESGB gravity would
likely be similar to those obtained by Ref. [170], which found the result to
be robust with respect to the inclusion of higher-order PN corrections. On the
other hand, slightly improved bounds could be obtained for heavier binaries
when using SEOBNRv5PHM$\\_$ESGB . However, the result obtained by combining
several events would likely improve only marginally, as GW230529 provides a
much more stringent bound compared to other events. This is a computationally
expensive but still feasible analysis that we leave for future work. Another
way to improve the combined constraint we present here would be to incorporate
information about the underlying astrophysical population following Ref.
[258], by simultaneously inferring the astrophysical distribution to correctly
account for correlations between intrinsic and beyond-GR parameters.

Event |  GW230529 |  GW190814 |  GW190412
---|---|---|---
IMR | $-0.94$ | $-0.02$ | $-0.52$
IMR with merger deviations | — | $-3.0$ | $-4.1$
Table 4: Natural log Bayes factors between the dilatonic ESGB and GR
hypotheses from the analysis of GW230529, GW190814 and GW190412, using the
SEOBNRv5PHM$\\_$ESGB model. We consider the IMR analyses in ESGB, both with
and without parameterized deviations to the waveform amplitude and frequency
at merger.

##  VII Conclusions

Using the pySEOBNR code [167], we built the first IMR waveform model in a
beyond-GR theory within the EOB approach, focusing on ESGB gravity for BBHs.
The SEOBNRv5PHM$\\_$ESGB waveforms reduce to the SEOBNRv5PHM model [157] in
the GR limit. The EOB Hamiltonian and fluxes in the pySEOBNR code can be used
to build IMR waveforms for BBHs and inspiral waveforms for BNS and NSBHs in
any ESGB (or ST) theory. The EOB Hamiltonian [cf. Sec. III] includes ESGB (and
ST) corrections through 3PN order [68], while the metric modes in factorized
form and the scalar flux entering the RR force contain ESGB (and ST)
corrections through 2PN and 2.5PN (relative to the scalar dipolar) orders,
respectively [cf. Sec. IV]. We recall that the factorized ESGB (and ST) metric
modes [cf. Sec. V.1] are a new result of this paper. While the ESGB inspiral
waveform is theory specific, the merger-ringdown waveform follows
phenomenological ansätze used in GR [211, 210, 212, 156], augmented with an
agnostic parameterization for the modes’ amplitude and frequency at merger,
inspired by Refs. [25, 30, 166]. It also uses corrections to the QNMs
frequencies and decay times derived in ESGB theory [138] and corrections to
the remnant mass and spin, which we find can cause changes in the frequency of
the ringdown waveform from one to two orders of magnitude larger than the ESGB
corrections to the QNM spectrum. Although our model can be applied to
different ESGB coupling functions, we focused our investigation to dilatonic
ESGB gravity.

In Sec. VI.1, we analyzed the morphology of the dilatonic ESGB waveforms, and
contrasted them with the GR ones, finding that the ESGB corrections in the
conservative and dissipative dynamics cause a BBH to merge in a shorter time
than in GR, when starting at the same orbital frequency (see Fig. 8). We
provided measurability estimates of the ESGB coupling [cf. Sec. VI.2] based on
the mismatch criterion for upcoming runs of the LVK collaboration and future
ground-based detectors, such as the Einstein Telescope [cf. Fig. 9]. The bound
from a single asymmetric ($q\sim 8$) BBH observed with ET could be one order
of magnitude better than current constraints using LVK data. Low-mass NSBHs
could also be “golden” binary sources, as exemplified by GW230529. This
indicates that stellar-mass binaries observed with next-generation ground-
based detectors are ideal candidates to test this theory, also thanks to their
high expected rates [259], and would likely provide better constraints than
extreme-mass-ratio-inspirals seen with LISA [260].

We also used SEOBNRv5PHM$\\_$ESGB to place constraints on the coupling of the
dilatonic ESGB theory using GW observations and performed Bayesian model
selection between ESGB and GR [cf. Sec. VI.3]. We found that the GW events
GW190412, GW190814 and GW230529 are consistent with GR, and we constrained the
fundamental coupling of the theory to be $\ell_{\mathrm{GB}}\lesssim
1.65\leavevmode\nobreak\ \mathrm{km}$ ($\sqrt{\alpha_{\mathrm{GB}}}\lesssim
0.31\leavevmode\nobreak\ \mathrm{km}$) at 90 $\%$ confidence. Among the events
we analyzed, GW230529 provides the most stringent constraint, consistently
with other recent analyses [171, 172].

It will be very beneficial to validate SEOBNRv5PHM$\\_$ESGB against NR
simulations in ESGB gravity and to calibrate it using NR data, particularly by
including NR information for the merger phase, and by guiding resummations of
the PN scalar flux. Despite remarkable progress in NR evolutions of BBHs in
ESGB gravity [80, 84, 85, 83, 87, 88, 91, 100, 99] and NSBH binaries [104],
the number and length of current simulations, as well as their accuracy, have
so far prevented us from carrying out these studies.

We limited our investigation to BHs that have already scalarized. However,
ESGB gravity also allows for dynamical scalarization [135, 72]. Interesting
work has been done to model dynamical scalarization [133, 134, 135, 136], and
we plan to extend SEOBNRv5PHM$\\_$ESGB to include these nonperturbative
effects, which can largely affect the waveforms toward merger, producing
distinctive signatures that are absent in GR.

Furthermore, it should be straightforward to include in our waveform model
ESGB corrections in the spinning sector, once they are available. Another
extension of our model would involve describing complete waveforms for BNSs
and NSBHs in both ST and ESGB theories, once these models have been finalized
in GR within the SEOBNRv5 family.

On the data-analysis side, the simple forecast we presented for upcoming runs
of the LVK collaboration and next-generation ground-based GW detectors can be
extended to a broader region of the binary’s parameter space by using linear-
signal approximation methods, such as the Fisher matrix formalism.
Furthermore, they can be made more realistic by stacking many events and
taking into account astrophysical populations [258]. Results for selected
events can be validated by performing synthetic injection and analyzing them
with Bayesian inference. Furthermore, our model can be used to validate
results obtained with parametrized tests of GR [23, 31] and to test their
effectiveness in different regions of the parameter space. It would also be
interesting to analyze more LVK events, including the NSBHs GW200105, GW200115
and selected low-mass BBHs, to improve the combined constraints on $\ell_{\rm
GB}$. However, the results would likely only improve marginally, as GW230529
provides a much more stringent bound compared to other events.

## Acknowledgments

It is a pleasure to thank Maxence Corman, Quentin Henry, Elisa Maggio, Elise
Sänger, Hector Okada da Silva, and David Trestini for useful discussions. We
also thank Raffi Enficiaud for his assistance with scientific computing.

The computational work for this manuscript was carried out on the Hypatia
computer cluster at the Max Planck Institute for Gravitational Physics in
Potsdam. This work made use of the python package pySEOBNR [167]. Stable
versions of pySEOBNR are published through the Python Package Index (PyPI),
and can be installed via pip install pyseobnr.

This research used data or software obtained from the Gravitational Wave Open
Science Center (gwosc.org), a service of the LIGO Scientific Collaboration,
the Virgo Collaboration, and KAGRA. This material is based upon work supported
by NSF’s LIGO Laboratory, a major facility fully funded by the National
Science Foundation, as well as the Science and Technology Facilities Council
(STFC) of the United Kingdom, the Max-Planck-Society (MPS), and the State of
Niedersachsen/Germany for support of the construction of Advanced LIGO and
construction and operation of the GEO600 detector. Additional support for
Advanced LIGO was provided by the Australian Research Council. Virgo is
funded, through the European Gravitational Observatory (EGO), by the French
Centre National de Recherche Scientifique (CNRS), the Italian Istituto
Nazionale di Fisica Nucleare (INFN) and the Dutch Nikhef, with contributions
by institutions from Belgium, Germany, Greece, Hungary, Ireland, Japan,
Monaco, Poland, Portugal, Spain. KAGRA is supported by Ministry of Education,
Culture, Sports, Science and Technology (MEXT), Japan Society for the
Promotion of Science (JSPS) in Japan; National Research Foundation (NRF) and
Ministry of Science and ICT (MSIT) in Korea; Academia Sinica (AS) and National
Science and Technology Council (NSTC) in Taiwan.

##  Appendix A Einstein and Jordan frames

In Refs. [47, 57, 58], the ST two-body Lagrangian was derived up to 3PN order
by adopting the Jordan-frame formulation of the theory:

$\displaystyle I_{\rm
ST}=\int\frac{d^{4}x\sqrt{-\tilde{g}}}{16\pi}\left(\phi\tilde{R}-\frac{\omega(\phi)}{\phi}(\partial\phi)^{2}\right)+I_{\rm
m}[\Psi,\tilde{g}_{\mu\nu}],$ (79)

using tildes for clarity. Here,
$(\partial\phi)^{2}=\tilde{g}^{\mu\nu}\partial_{\mu}\phi\partial_{\nu}\phi$,
and $\omega(\phi)$ is a function defining the theory. As for compact bodies,
they were described by the point particle action

$\displaystyle I_{\rm m}\to I_{\rm m}^{\rm
pp}=-\sum_{A}\int\tilde{m}_{A}(\phi)d\tilde{s}_{A},$ (80)

with $d\tilde{s}_{A}=\sqrt{-\tilde{g}_{\mu\nu}dx_{A}^{\mu}dx_{A}^{\nu}}$. By
contrast, we adopt here the Einstein-frame formulation of ST theories:

$\displaystyle I_{\rm
ST}=\int\frac{d^{4}x\sqrt{-g}}{16\pi}\left(R-2(\partial\varphi)^{2}\right)+I_{\rm
m}[\Psi,\mathcal{A}^{2}g_{\mu\nu}],$ (81)

where
$(\partial\varphi)^{2}=g^{\mu\nu}\partial_{\mu}\varphi\partial_{\nu}\varphi$,
and we account for compact bodies by the substitution

$\displaystyle I_{\rm m}\to I_{\rm m}^{\rm pp}=-\sum_{A}\int
m_{A}(\varphi)ds_{A},$ (82)

with $ds_{A}=\sqrt{-g_{\mu\nu}dx_{A}^{\mu}dx_{A}^{\nu}}$. The actions
(79)-(80) and (81)-(82) are the same, modulo boundary terms and the following
redefinitions:

$\displaystyle\tilde{g}_{\mu\nu}$ $\displaystyle=\mathcal{A}^{2}g_{\mu\nu},$
(83a) $\displaystyle 3+2\omega(\phi)$
$\displaystyle=\left(\frac{d\ln\mathcal{A}}{d\varphi}\right)^{-2},$ (83b)
$\displaystyle m_{A}(\varphi)$
$\displaystyle=\mathcal{A}(\varphi)\tilde{m}_{A}(\varphi),$ (83c)

where $\varphi(\phi)$ is obtained by inverting
$\mathcal{A}(\varphi)=1/\sqrt{\phi}$. We also introduce the notation:

$\displaystyle\alpha_{0}$
$\displaystyle=\frac{d\ln\mathcal{A}}{d\varphi}(\varphi_{0}),$ (84a)
$\displaystyle\beta_{0}$
$\displaystyle=\frac{d\alpha}{d\varphi}(\varphi_{0}),$ (84b)
$\displaystyle{\beta^{\prime}}_{0}$
$\displaystyle=\frac{d\beta}{d\varphi}(\varphi_{0}),$ (84c)
$\displaystyle{\beta^{\prime\prime}}_{0}$
$\displaystyle=\frac{d{\beta^{\prime}}}{d\varphi}(\varphi_{0}),$ (84d)

where the subscript $0$ denotes a quantity evaluated at infinity,
$\varphi(\phi_{0})=\varphi_{0}$. The quantities above can be obtained by
inserting Eq. (83c) into Eqs. (4) and taking the limit
$\tilde{m}_{A}(\phi)=\rm const.$, in which case body $A$ is said to have
negligible self-gravity.

We can now translate the parameters of Refs. [47, 57, 58] into our conventions
using Eqs. (83) and below. The results are gathered in Table 5.

Refs. [47, 51, 57, 58, 64] | Refs. [43, 176, 52, 53, 61, 68] and this paper
---|---
Theory dependent |
$\tilde{G}$ | $\mathcal{A}_{0}^{2}(1+\alpha_{0}^{2})$
$\zeta$ | $\frac{\alpha_{0}^{2}}{1+\alpha_{0}^{2}}$
$\lambda_{1}$ | $\frac{\beta_{0}}{2(1+\alpha_{0}^{2})}$
$\lambda_{2}$ | $\frac{1}{4(1+\alpha_{0}^{2})^{2}}(-\beta^{\prime}_{0}\alpha_{0}-2\beta_{0}\alpha_{0}^{2}+4\beta_{0}^{2})$
$\lambda_{3}$ | $\frac{1}{8(1+\alpha_{0}^{2})^{3}}(6\beta^{\prime}_{0}\alpha_{0}^{3}-24\beta_{0}^{2}\alpha_{0}^{2}+8\beta_{0}\alpha_{0}^{4}-13\beta^{\prime}_{0}\beta_{0}\alpha_{0}+24\beta_{0}^{3}+{\beta^{\prime\prime}}_{0}\alpha_{0}^{2})$
Body dependent |
$m_{A}$ | $m_{A}^{0}/\mathcal{A}_{0}$
$s_{A}$ | $\frac{1}{2}(1-\alpha_{A}^{0}/\alpha_{0})$
$s_{A}^{\prime}$ | $\frac{1}{4\alpha_{0}}(\beta_{A}^{0}/\alpha_{0}-\alpha_{A}^{0}\beta_{0}/\alpha_{0}^{2})$
$s_{A}^{\prime\prime}$ | $\frac{1}{8\alpha_{0}^{5}}(3\beta_{A}^{0}\alpha_{0}\beta_{0}-3\alpha_{A}^{0}\beta_{0}^{2}-{\beta^{\prime}}_{A}^{0}\alpha_{0}^{2}+\alpha_{A}^{0}\beta^{\prime}_{0}\alpha_{0})$
$s_{A}^{\prime\prime\prime}$ |  $\frac{1}{16\alpha_{0}^{7}}(15\beta_{A}^{0}\alpha_{0}\beta_{0}^{2}-15\alpha_{A}^{0}\beta_{0}^{3}-5{\beta^{\prime}}_{A}^{0}\beta_{0}\alpha_{0}^{2}+10\alpha_{A}^{0}\beta_{0}^{\prime}\beta_{0}\alpha_{0}-{\beta^{\prime}}_{A}^{0}\alpha_{0}^{2}\beta_{0}$ $-4\beta_{A}^{0}\alpha_{0}^{2}\beta_{0}^{\prime}+{\beta^{\prime\prime}}_{A}^{0}\alpha_{0}^{3}-\alpha_{A}^{0}{\beta^{\prime\prime}}_{A}^{0}\alpha_{0}^{2})$
0PN |
$\tilde{G}\tilde{\alpha}$ | $\mathcal{A}_{0}^{2}G_{AB}=\mathcal{A}_{0}^{2}(1+\alpha_{A}^{0}\alpha_{B}^{0})$
1PN |
$\bar{\gamma}$ | $\bar{\gamma}_{AB}=\frac{-2\alpha_{A}^{0}\alpha_{B}^{0}}{1+\alpha_{A}^{0}\alpha_{B}^{0}}$
$\bar{\beta}_{A}$ | $\bar{\beta}_{A}=\frac{\beta_{A}^{0}(\alpha_{B}^{0})^{2}}{2(1+\alpha_{A}^{0}\alpha_{B}^{0})^{2}}$
2PN |
$\bar{\delta}_{A}$ | $\delta_{A}=\frac{(\alpha_{A}^{0})^{2}}{(1+\alpha_{A}^{0}\alpha_{B}^{0})^{2}}$
$-4\,\bar{\chi}_{A}$ | $\epsilon_{A}=\frac{{\beta^{\prime}}_{A}^{0}(\alpha_{B}^{0})^{3}}{(1+\alpha_{A}^{0}\alpha_{B}^{0})^{3}}$
$-8\,\bar{\beta}_{A}\bar{\beta}_{B}/\bar{\gamma}$ | $\zeta_{AB}=\frac{{\beta}_{A}^{0}{\beta}_{B}^{0}\alpha_{A}^{0}\alpha_{B}^{0}}{(1+\alpha_{A}^{0}\alpha_{B}^{0})^{3}}$
3PN |
$\bar{\kappa}_{A}$ | $\kappa_{A}=\frac{(\alpha_{B}^{0})^{4}{\beta^{\prime\prime}_{A}}^{0}}{8(1+\alpha_{A}^{0}\alpha_{B}^{0})^{4}}$
$-4\,\bar{\beta}_{A}\bar{\delta}_{A}/\bar{\gamma}$ | $\psi_{A}=\frac{\alpha_{A}^{0}\alpha_{B}^{0}\beta_{A}^{0}}{(1+\alpha_{A}^{0}\alpha_{B}^{0})^{3}}$
$16\,\bar{\beta}_{A}\bar{\chi}_{B}/\bar{\gamma}$ | $\xi_{A}=\frac{(\alpha_{A}^{0})^{2}\alpha_{B}^{0}\beta_{A}^{0}{\beta^{\prime}}_{B}^{0}}{(1+\alpha_{A}^{0}\alpha_{B}^{0})^{4}}$
$32\,\bar{\beta}_{A}(\bar{\beta}_{B})^{2}/\bar{\gamma}^{2}$ | $\omega_{A}=\frac{(\alpha_{A}^{0})^{2}\beta_{A}^{0}(\beta_{B}^{0})^{2}}{(1+\alpha_{A}^{0}\alpha_{B}^{0})^{4}}$
Other |
$\sqrt{\zeta}\,\mathcal{S}_{\pm}$ | $\frac{\alpha_{\pm}}{2\sqrt{1+\alpha_{A}^{0}\alpha_{B}^{0}}}$
Table 5: Translation of the parameters from Refs. [47, 51, 57, 58, 64]. We
rename their $\alpha$ as $\tilde{\alpha}$ here to avoid confusing it with
ours.

##  Appendix B ESGB corrections to the EOB potentials

The coefficients $(\delta\bar{a}_{3},\,\delta\bar{a}_{4},\,\bar{a}_{4,\ln})$,
$(\delta\bar{d}_{2},\,\delta\bar{d}_{3})$ and
$(\delta\bar{q}_{1},\,\bar{q}_{2})$ entering Eqs. (14) were calculated in
Refs. [52, 68]. They read:

$\displaystyle\delta\bar{a}_{3}$
$\displaystyle=\frac{1}{12}\Big{[}\nu(-24\zeta_{AB}-36\bar{\beta}_{+}+4\bar{\gamma}_{AB}^{2}+40\bar{\gamma}_{AB}+8\delta_{+}+4\epsilon_{+})-24\langle\bar{\beta}\rangle\left(1-2\bar{\gamma}_{AB}\right)-35\bar{\gamma}_{AB}^{2}-20\bar{\gamma}_{AB}+4\langle\delta\rangle-4\langle\epsilon\rangle\Big{]},$
(85a) $\displaystyle\delta\bar{a}_{4}$ $\displaystyle=2k_{\rm ESGB}-4k_{\rm
tail}(\gamma_{E}+\ln
2)+\frac{\nu\bar{\gamma}_{AB}\big{(}11\left(\bar{\gamma}_{AB}+2\right)^{2}-4\langle\delta\rangle\big{)}}{4\tilde{\alpha}\left(\bar{\gamma}_{AB}+2\right)}$
$\displaystyle+\frac{1}{12}\Big{[}2\langle\bar{\beta}\rangle\left(47\bar{\gamma}_{AB}^{2}-28\bar{\gamma}_{AB}+6\bar{\beta}_{+}+4\delta_{+}-12\langle\delta\rangle-28\right)-60\bar{\gamma}_{AB}^{3}-78\bar{\gamma}_{AB}^{2}-24\bar{\gamma}_{AB}+\langle\delta\rangle\left(16\bar{\gamma}_{AB}+8\right)-12\langle\epsilon\rangle\bar{\gamma}_{AB}$
$\displaystyle-4\delta_{-}\bar{\beta}_{-}-4\delta_{+}\bar{\beta}_{+}+3\bar{\beta}_{-}^{2}-3\bar{\beta}_{+}^{2}+8\langle\delta\rangle\bar{\beta}_{+}+60\langle\bar{\beta}\rangle^{2}+8\langle\kappa\rangle\Big{]}$
$\displaystyle+\frac{\nu}{1152}\Big{[}-288\bar{\beta}_{+}\left(40\bar{\gamma}_{AB}-3\right)+192\langle\bar{\beta}\rangle\left(-4\bar{\gamma}_{AB}^{2}-64\bar{\gamma}_{AB}+30\bar{\beta}_{+}-8\delta_{+}-19\right)$
$\displaystyle+\delta_{+}\left(126\pi^{2}\bar{\gamma}_{AB}+768\bar{\gamma}_{AB}+252\pi^{2}-6080\right)-6912\zeta_{AB}\bar{\gamma}_{AB}+1152\epsilon_{+}\bar{\gamma}_{AB}+63\pi^{2}\bar{\gamma}_{AB}^{3}-432\bar{\gamma}_{AB}^{3}-180\pi^{2}\bar{\gamma}_{AB}^{2}$
$\displaystyle+15296\bar{\gamma}_{AB}^{2}-1350\pi^{2}\bar{\gamma}_{AB}+37184\bar{\gamma}_{AB}+192\langle\delta\rangle\left(3\bar{\gamma}_{AB}+2\right)+1536\delta_{-}\bar{\beta}_{-}-864\bar{\beta}_{-}^{2}-288\bar{\beta}_{+}^{2}+3456\zeta_{AB}$
$\displaystyle+1152\langle\bar{\beta}\rangle^{2}+1152\langle\epsilon\rangle-768\kappa_{+}-768\langle\kappa\rangle-1536\langle\psi\rangle+1152\langle
w\rangle-768\psi_{+}+1152\langle\xi\rangle\Big{]},$ (85b)
$\displaystyle\bar{a}_{4,\ln}$ $\displaystyle=-2k_{\rm tail},$ (85c)

$\displaystyle\delta\bar{d}_{2}$
$\displaystyle=\frac{1}{4}\Big{[}-3\bar{\gamma}_{AB}^{2}-12\bar{\gamma}_{AB}+4\langle\delta\rangle-24\langle\bar{\beta}\rangle+8\nu(2\bar{\gamma}_{AB}-\langle\bar{\beta}\rangle)\Big{]},$
$\displaystyle\delta\bar{d}_{3}$ $\displaystyle=k_{\rm tail}(21-32\ln
2)+\frac{1}{12}\Big{[}4(\langle\delta\rangle-6\langle\bar{\beta}\rangle)\left(3\bar{\gamma}_{AB}+8\right)-9\bar{\gamma}_{AB}^{3}-52\bar{\gamma}_{AB}^{2}-64\bar{\gamma}_{AB}+16\langle\epsilon\rangle\Big{]}$
$\displaystyle+\frac{\nu}{12}\Big{[}4\left(-\delta_{+}\left(6\bar{\gamma}_{AB}+11\right)+9\bar{\gamma}_{AB}^{3}+77\bar{\gamma}_{AB}^{2}+173\bar{\gamma}_{AB}+27\zeta_{AB}-6\langle\delta\rangle+3\langle\epsilon\rangle-4\epsilon_{+}\right)$
$\displaystyle+3\bar{\beta}_{+}\left(4\bar{\gamma}_{AB}+69\right)-18\langle\bar{\beta}\rangle\left(8\bar{\gamma}_{AB}+11\right)\Big{]}+\nu^{2}\Big{[}-\bar{\gamma}_{AB}^{2}-10\bar{\gamma}_{AB}+9\bar{\beta}_{+}+6\zeta_{AB}-2\delta_{+}-\epsilon_{+}\Big{]},$
(86a)

and, finally,

$\displaystyle\delta\bar{q}_{1}$ $\displaystyle=\frac{1}{6}k_{\rm
tail}(93+1753\ln 2-729\ln