# Measuring scalar charge with compact binaries: High accuracy modelling with
self-force

Andrew Spiers  Nottingham Centre of Gravity & School of Mathematical Sciences,
University of Nottingham, University Park, Nottingham, NG7 2RD, UK  Andrea
Maselli  Gran Sasso Science Institute (GSSI), I-67100 L’Aquila, Italy  INFN,
Laboratori Nazionali del Gran Sasso, I-67100 Assergi, Italy  Thomas P.
Sotiriou  Nottingham Centre of Gravity & School of Mathematical Sciences,
University of Nottingham, University Park, Nottingham, NG7 2RD, UK  School of
Physics and Astronomy, University of Nottingham, University Park, Nottingham,
NG7 2RD, UK

###### Abstract

Using the self-force approach, we present the premier first-post-adiabatic
accuracy formalism for modelling compact binaries in theories with a massless
scalar field non-minimally coupled to gravity. We limit the binary secondary
to being a non-spinning compact body with no scalar dipole (we will address
the spinning and scalar dipole cases in an upcoming paper). By producing an
ansatz for the scalar charged point particle action, we derive first- and
second-order perturbative field equations and equations of motion for the
secondary compact object. Under our assumptions, implementing this formalism
will produce sufficiently accurate waveform templates for precision
measurements of the scalar charge of the secondary with LISA data on extreme-
mass-ratio inspirals. Our formalism is consistent with almost general scalar-
tensor theories of gravity. Implementing our formalism builds on self-force
models in General Relativity; we show the incorporation into the two-timescale
formalism is straightforward. Excitingly, implementation poses no
significantly more challenging barriers than computing first-post adiabatic
waveforms in General Relativity.

##  I  Introduction

The detection of a new fundamental field via its imprint on compact objects
and the gravitational wave signals they produce would be a ground-breaking
discovery. Indeed, such searches are among the key goals of current and future
gravitational wave detectors [1, 2, 3, 4, 5]. Asymmetric binaries are
particularly promising sources in this context [6, 7, 8, 9, 10, 11, 12, 13,
14, 15, 16, 17]. They consist of a larger body, the primary, of mass $M$, and
a much smaller body, the secondary, of mass $\mu$111In this paper $\mu$ is the
zeroth-order in $\varepsilon$ mass, we allow for higher order corrections to
the mass which come from the presence of the scalar field.. Their mass ratio

$\displaystyle\varepsilon=\frac{\mu}{M}\ ,$ (1)

is then very small, with the more extreme cases, known as extreme mass ratio
inspirals (EMRIs), reaching $\varepsilon\sim 10^{-6}$. EMRIs could execute
some $10^{5}$ orbits while in LISA’s band and be continuously observed for
very long periods — several months or even years. Hence, they are expected to
be excellent probes of the properties of their source [18].

Until recently attention was strongly focused on EMRIs potential to probe the
spacetime around, and hence the properties of, the primary black hole [19, 20,
21]. However, it was recently pointed out in [22] that the scalar charge of
the secondary could leave a strong imprint on the EMRI waveforms — strong
enough to make the charge measurable by LISA [23]. It has also been shown in
[22], using an effective field theory framework, that, to leading order in
$\varepsilon$, the charge of the primary is negligible. Hence, the primary is
adequately approximated by a Kerr black hole, and the charge for the secondary
fully controls the deviations from GR. This drastically simpler framework has
already been extended to eccentric orbits [24] and massive scalar fields [25].
We will discuss the assumptions underpinning this framework in detail below,
but its most important practical limitation is that it does not include any
post-adiabatic corrections that appear at second order in $\varepsilon$.

Indeed, a burgeoning method for modelling compact binaries is the self-force
approach in black hole perturbation theory. This method provides an accurate
approximation in the extreme-mass-ratio limit, $\varepsilon\lesssim 10^{-4}$,
where it is exorbitantly computationally expensive to produce full inspiral
waveforms using Numerical Relativity due to the binary’s disparate length
scales. Tackling this problem using perturbation theory is advantageous
because $\varepsilon$ is a naturally small expansion parameter. World-leading
perturbative self-force models are to first-post-adiabatic accuracy. Such
accuracy provides waveforms with $\mathcal{O}(\varepsilon)$ phase error over
the course of an inspiral of $\mathcal{O}(\frac{1}{\varepsilon})$ orbits.
Currently, waveforms to first-post adiabatic accuracy have been computed for
quasi-circular inspirals of Schwarzschild black holes in GR. These waveforms
show impressive agreement with Numerical Relativity waveforms for mass-ratios
smaller than $\mathcal{O}(\frac{1}{10})$. Significant efforts are ongoing to
extend these results to generic orbits around a Kerr primary black hole [26,
27, 28, 29, 30, 31, 32, 33], and include effects such as resonances [34] and
the spin of the secondary [35, 36, 37, 38, 39].

Extending the perturbative self-force approach to theoretical scenarios that
include new fundamental fields is a fledgling field of research of clear
importance. Constraining the existence of new fundamental fields with LISA
measurements on EMRIs requires high-accuracy waveform templates, including the
effects of these new fundamental fields. As in GR, one has to reach first-post
adiabatic accuracy. As a first, critical step in this direction, in this
paper, we push the approach of [22, 23] to post-adiabatic order. By
generalising the perturbative self-force approach and using a new ansatz for
the point-particle action, we derive the field equations and equations of
motion for the secondary. That is, we provide the required equations for
modelling inspirals to first-post adiabatic accuracy in scenarios with a
massless scalar field nonminimally coupled to gravity. This is the first
scheme for modelling binaries perturbatively to first-post adiabatic accuracy
beyond GR and the Standard Model and provides a roadmap to full calculations.

The paper is organised as follows. In Sec. II, we discuss the theoretical
ground in which we model asymmetric binaries beyond GR. In Sec. III, we derive
a formalism for calculating the first- and second-order self-force in a large
class of scalar-tensor theories of gravity. We exploit our formalism to build
a two-timescale approximation for efficient first-post adiabatic accurate
modelling in Sec. V. In appendix A, we briefly review the self-force approach
within black hole perturbation theory in GR, derive the metric perturbation
field equations and perturbative equations of motion from our action approach,
and demonstrate why the first-post adiabatic models provide high-accuracy
waveforms ready for LISA observations.

##  II Asymmetric binaries and scalar fields

###  II.1 Action and field equations

Following Ref. [22], our starting point will be a general action which
describes a real scalar field $\varphi$ non-minimally coupled to gravity:

$\displaystyle
S[\mathbf{g}_{ab},\varphi,\Psi]=S_{0}[\mathbf{g}_{ab},\varphi]+\alpha
S_{\mathrm{c}}[\mathbf{g}_{ab},\varphi]+S_{\mathrm{m}}[\mathbf{g}_{ab},\varphi,\Psi]\
,$ (2)

where

$\displaystyle
S_{0}[\mathbf{g}_{ab},\varphi]=\int\frac{\sqrt{-\mathbf{g}}}{16\pi}\Big{(}R-\frac{1}{2}\partial_{\mu}\varphi\partial^{\mu}\varphi\Big{)}d^{4}x\
,$ (3)

$\mathbf{g}$ is the metric determinant, and we work in geometric units
$G=c=\hbar=1$. $S_{\mathrm{c}}$ contains any additional (self)-interactions
for the scalar, whereas $S_{\mathrm{m}}$ describes the matter fields $\Psi$.
For concreteness, we consider here a massless scalar, omit the bare mass term
in $S_{0}$, and consistently assume that $S_{\mathrm{c}}$ respects shift
symmetry. However, this approach can be straightforwardly generalized to
massive scalars [25], and we will return to this point later.

We take the coupling constant $\alpha$ to have dimensions $[{\rm length}]^{n}$
where $n\geq 2$. This corresponds to negative mass dimensions in particle
physics units; i.e., we assume that these interactions are suppressed by some
characteristic energy scale. Note that there are no terms that would
correspond to $n\leq 1$ that are consistent with local Lorentz symmetry [40].
We apply no other conditions to $S_{\mathrm{c}}$, and hence, our approach
covers a very broad range of theoretical scenarios. For the time being, we
leave $S_{m}$ generic.

To derive field equations for the metric, one varies the action (2) with
respect to the metric, that yields

$\displaystyle G_{ab}[\mathbf{g}_{cd}]=T_{ab}^{\mathrm{scal}}[\varphi]+\alpha
T_{ab}^{\mathrm{c}}[\mathbf{g}_{cd},\varphi]+T_{ab}^{\mathrm{m}}[\mathbf{g}_{cd},\varphi,\Psi]\
,$ (4)

where $G_{ab}[\mathbf{g}_{cd}]$ is the Einstein operator,

$\displaystyle T^{(i)}_{ab}=-\frac{8\pi}{\sqrt{-\mathbf{g}}}\frac{\delta
S_{(i)}}{\delta g^{ab}}\ ,$ (5)

where $i\in\\{\mathrm{c},\mathrm{m}\\}$, and

$\displaystyle T_{ab}^{\mathrm{scal}}[\varphi]$
$\displaystyle=\frac{1}{2}\partial_{a}\varphi\partial_{b}\varphi-\frac{1}{4}\mathbf{g}_{ab}\partial_{c}\varphi\partial^{c}\varphi\
.$ (6)

To derive the scalar field equation one varies Eq. (2) with respect to
$\varphi$, yielding

$\displaystyle\Box_{\mathbf{g}}\varphi=\alpha\Sigma_{\mathrm{c}}+\Sigma\,,$
(7)

where
$\Box_{\mathbf{g}}=\mathbf{g}^{ab}\nabla^{\mathbf{g}}_{a}\nabla^{\mathbf{g}}_{b}$,
$\nabla^{\mathbf{g}}_{a}$ is the covariant derivative associated with the
metric $\mathbf{g}_{ab}$, and

$\displaystyle\Sigma=-\frac{16\pi}{\sqrt{-\mathbf{g}}}\frac{\delta
S_{\mathrm{m}}}{\delta\varphi},\qquad\Sigma_{\mathrm{c}}=-\frac{16\pi}{\sqrt{-\mathbf{g}}}\frac{\delta
S_{\mathrm{c}}}{\delta\varphi}.$ (8)

###  II.2 The primary

We assume that the primary is a black hole. The action $S_{0}$ is covered by
no-hair theorems [41, 42], so any scalar hair would have to be introduced by
terms in $S_{\mathrm{c}}$. For shift-symmetric theories, however,
$S_{\mathrm{c}}$ is also covered by the no-hair theorem, for static and
spherically symmetric [43] and for slowly rotating [44] asymptotically flat
black holes. The only interaction that evades this theorem is a linear
coupling between $\varphi$ and the Gauss-Bonnet invariant ${\cal G}\equiv
R_{abcd}R^{abcd}-4R_{ab}R^{ab}+R^{2}$ [45]. Adding the term $\alpha_{\rm
GB}\,\varphi\,{\cal G}$ to $S_{0}$ introduces a non-constant scalar field to
black holes, described by the black hole having a scalar charge; however, the
scalar charge (which can be thought of as the scalar monopole) is not an
independent parameter, instead, the scalar charge is fixed by a regularity
condition on the horizon and is determined by the mass and spin of the black
hole and $\alpha_{\rm GB}$ [45, 44, 46] (see also [47] for earlier work
without shift symmetry). The charge per unit mass $d$ scales as $\alpha_{\rm
GB}/M^{2}$ in geometric units [44].

Adding additional shift-symmetric interactions will change the scalar
configuration [48], but the regularity conditions that fix the charge
persists. The charge per unit mass is then given by an integral over the
horizon ${\cal H}$, [40]

$\displaystyle d=\frac{\alpha_{\rm GB}}{4\pi M}\int_{\cal H}n^{a}{\cal G}_{a}\
d\Omega$ (9)

where $n^{a}$ is the horizon generator and ${\cal G}=\nabla_{a}{\cal G}^{a}$.
This implies that any terms in $S_{\mathrm{c}}$ other than the linear coupling
with ${\cal G}$, controlled by a coupling $\alpha_{i}$, contribute to $d$ with
a factor of $\alpha_{\rm GB}\alpha_{i}$.

###  II.3 The secondary

To initially define the matter action for (2), which describes the secondary
body in the EMRI, we use the conventional skeletonized approach222The
skeletonized formalism was first developed for electromagnetism and gravity
[49] and has previously been extended to scalar-tensor theories with multiple
fields [50]. [50]. The skeletonized description of a compact object replaces
the matter action, $S_{\mathrm{m}}$ (assuming no other matter fields are
present), with a point particle action $S_{\mathrm{p}}$. Ref. [50] presented a
point-particle action for a massive, scalar-charged, compact object:

$\displaystyle
S_{\mathrm{p}}=-\int_{\gamma}m[\varphi]ds=-\int_{\gamma}m[\varphi]\sqrt{\mathbf{g}_{ab}\mathbf{u}^{a}\mathbf{u}^{b}}d\tau,$
(10)

where $\gamma$ is the worldline of the compact object,
$\mathbf{u}^{\alpha}=\frac{dz^{\alpha}}{d\tau}$ is the four-velocity of the
compact object in $\mathbf{g}_{ab}$ and $\tau$ is the proper time in
$\mathbf{g}_{ab}$. Eq. (10) introduces a mass function $m[\varphi]$, which
depends on the scalar field, generating the scalar charge. We will show that
Eq. (10) is sufficient for deriving the linear field equations but encounters
issues beyond linear order.

With our point-particle action in hand, we can now derive the stress-energy
tensor and scalar charge density that will appear in the field equations.
Varying Eq. (10) with respect to the metric, $\mathbf{g}_{ab}$, yields

$\displaystyle
T_{\mathrm{m}}^{ab}=8\pi\int_{\gamma}m[\varphi]\frac{\delta^{4}[x^{\mu}-z_{\mathrm{p}}^{\mu}[\tau]]}{\sqrt{-\mathbf{g}}}\mathbf{u}^{a}\mathbf{u}^{b}d\tau,$
(11)

the scalar charged point particle stress-energy tensor. Varying Eq. (10) with
respect to the scalar field, $\varphi$, yields

$\displaystyle\Sigma=16\pi\int_{\gamma}m^{\prime}[\varphi]\frac{\delta^{4}[x^{\mu}-z_{\mathrm{p}}^{\mu}[\tau]]}{\sqrt{-\mathbf{g}}}d\tau,$
(12)

the point particle scalar-charge density.

##  III Perturbative expansion

The contribution $\alpha$ makes to $\mathbf{g}_{ab}$ and $\varphi$ must be
dimensionless as $\mathbf{g}_{ab}$ and $\varphi$ are dimensionless. As $M$ is
the only length scale associated with the background spacetime, $\alpha$ must
be accompanied by a factor of $\frac{1}{M^{n}}$ for the leading-order
contribution. We introduce the dimensionless coupling

$\displaystyle\zeta=\frac{\alpha}{M^{n}}.$ (13)

This can be expressed in terms of the mass ratio $\varepsilon$ as

$\displaystyle\zeta=\frac{\alpha}{M^{n}}=\varepsilon^{n}\frac{\alpha}{\mu^{n}}\
.$ (14)

$\zeta$ represents the non-minimal coupling perturbation parameter for
$\mathbf{g}_{ab}$ and $\varphi$.

If we assume that solutions of the field equations are continuously connected
to GR solutions as $\alpha\to 0$, our earlier assumptions for
$S_{\mathrm{c}}$, that $[\alpha]=[{\rm length}]^{n}$, with $n\geq 2$, and the
expression for $d$ in Eq. (9), imply that deviations from GR are controlled by
$\zeta$. In particular, so long as the length-(energy)-scale associated with a
particular coupling $\alpha_{i}$ of a term in $S_{\mathrm{c}}$ ($\alpha$
denotes them collectively) is not significantly larger (smaller) than the
scale associated with $\alpha_{\rm GB}$, $d\propto\alpha_{\rm GB}/M^{2}$ to
leading order in $M^{-1}$. Existing constraints inferred from astrophysical
observations imply that $\frac{\alpha}{\mu^{n}}=\mathcal{O}(1)$ or smaller
[51], as $\mu$ corresponds to solar mass bodies. Therefore, the mass-ratio,
$\varepsilon$, the natural bookkeeping parameter for the self-force approach,
can be used as the sole perturbative parameter for the problem at hand. That
is, conservatively, $\zeta\sim\varepsilon^{n}$.

To build our formalism, we will consistently expand the field’s equations
(4)-(7), as well as the metric and the scalar field, up to the second order in
$\varepsilon$. An expansion of a generic tensor, ${\bf A}$, takes the form,

${\bf A}={\bf A}^{(0)}+\varepsilon{\bf A}^{(1)}+\varepsilon^{2}{\bf
A}^{(2)}+\mathcal{O}(\varepsilon^{3})\ .$ (15)

The metric expansion is written as

$\displaystyle\mathbf{g}_{ab}=g_{ab}+\varepsilon
h^{(1)}_{ab}+\varepsilon^{2}h^{(2)}_{ab}+\mathcal{O}(\varepsilon^{3}),$ (16)

where $g_{ab}$ is the background metric (which we take to be the Kerr metric)
and $h^{(n)}_{ab}$ are the metric perturbations. The background metric is used
to raise and lower the indices of all tensors and tensor perturbations, such
as $h^{(n)}_{ab}$. We label all the tensor perturbations with a subscript or
superscript number in brackets, which denotes the order in $\varepsilon$ of
the perturbation. In practice, the $\varepsilon$ dependence is implicit in the
labelled tensors, and the explicit factors of $\varepsilon$ in Eq. (16) are
used as a counting parameter. That is, $\varepsilon$ is set to $1$ before
computing calculations.

The scalar field expansion takes the form

$\displaystyle\varphi=\varphi^{(0)}+\varepsilon\varphi^{(1)}+\varepsilon^{2}\varphi^{(2)}+\mathcal{O}(\varepsilon^{3})\
.$ (17)

Note that $\varphi^{(0)}$ corresponds to the contribution from the action
$S_{0}$ for an isolated black hole, which is covered by no-hair theorems.
Hence, $\varphi^{(0)}$ is constant and can be set to zero by a constant shift
without loss of generality.

Our procedure requires a specific treatment for the mass function $m[\varphi]$
that appears in the secondary’s stress-energy tensor (23). To this aim we
expand $m[\varphi]$ as:

$\displaystyle
m[\varphi]=m_{[0]}+m_{[1]}\varphi+m_{[2]}\varphi^{2}+\mathcal{O}[\varphi^{3}],$
(18)

where $m_{[0]}$, $m_{[1]}$ and $m_{[2]}$ are constant coefficients. In our
setup, $m_{[0]}=\mu$; $\mu$ can be considered as the secondary mass in GR
(that is, for $\varphi=0$). Note $\mu$ is not the total mass of the secondary
in scalar-tensor theories of gravity as $m_{[1]}$ and $m_{[2]}$ can contribute
to the mass when $\varphi\neq 0$, as we will show. We assume $m_{[1]}$ and
$m_{[2]}$ have the same (stellar mass) scale of $m_{[0]}$; that is
$m_{[0,1,2]}/M=\mathcal{O}(\varepsilon)$.

We can derive the explicit expression of $m[\varphi]$ up to the second order
in $\varepsilon$ by replacing (17) (with $\varphi^{(0)}=0$) into Eq. (18),
obtaining:

$\displaystyle m[\varphi]$ $\displaystyle=\mu+\varepsilon
m_{[1]}\varphi_{(1)}+\varepsilon^{2}(m_{[2]}\varphi_{(1)}^{2}+m_{[1]}\varphi_{(2)})\
.$ (19)

###  III.1 Field regularization

We now return to the matter action and show how the conventional skeletonized
point particle action is problematic in the self-force context. Equation (10)
poses a typical problem within self-force: near the worldline ($\gamma$), the
singular nature of the metric and scalar field makes $\mathbf{g}_{ab}$ and
$m[\varphi]$ ill-defined, whereas the four-velocity is only defined on the
worldline [52]. Hence, Eq. (10) is ill-defined.

To solve this problem and extend Eq. (10) beyond linear order, we assume the
existence of a singular (${\cal S}$) and regular (${\cal R}$) split of the
metric and scalar field perturbations [53, 54]:
$h^{(n)}_{ab}=h^{(n)\mathcal{S}}_{ab}+h^{(n)\mathcal{R}}_{ab}$,
$\varphi^{(n)}=\varphi^{(n)\mathcal{S}}+\varphi^{(n)\mathcal{R}}$. A
motivation for this assumption is such a split exists in the decoupled case
[53]. We also define

$\displaystyle
h^{\mathcal{R}}_{ab}=\sum_{n=1}\varepsilon^{n}h^{\mathcal{R}(n)}_{ab}\quad,\quad\varphi^{\mathcal{R}}=\sum_{n=1}\varepsilon^{n}\varphi^{(n)}_{\mathcal{R}}.$
(20)

$h^{\mathcal{R}}_{ab}$ and $\varphi^{\mathcal{R}}$ identify the regular part
of the field perturbations (the part that generates the self-force) [52]. Ref.
[53] provides an extensive analysis of the singular-regular decomposition to
linear order, which has been generalised to the non-linear regime in GR [55,
56, 57].

The decomposition into singular and regular pieces of the metric perturbations
has been studied extensively in GR in [53]. Conventionally,
$G_{ab}^{(1)}[h^{(1)\mathcal{S}}_{cd}]=8\pi T^{\rm m(1)}_{ab}$ and
$G_{ab}^{(1)}[h^{(1)\mathcal{R}}_{cd}]=0$ are satisfied. This definition does
not fully fix the field; additionally, $h^{\mathcal{S}}_{ab}$ is chosen to
depend only on the instantaneous state and position of the particle, and
$h^{\mathcal{R}}_{ab}$ on the compact objects’ causal past [53]. We expect
similar definitions to hold for the metric perturbation and scalar field in
our formalism.

With Eqns. (20) in hand we define an effective metric and scalar field:

$\displaystyle\tilde{g}_{ab}=g_{ab}+h^{\mathcal{R}}_{ab}\quad,\quad\tilde{\varphi}=\varphi^{\mathcal{R}}\
.$ (21)

First and second-order self-force calculations in GR have found that compact
objects move as a test body of the effective metric [58, 59, 60, 61].

Replacing $\mathbf{g}_{ab}$ and $\phi$ with $\tilde{g}_{ab}$ and
$\tilde{\varphi}$ in Eq. (10), we obtain our effective point-particle action
of a scalar-charged compact object:

$\displaystyle
S_{\mathrm{p}}=-\int_{\gamma}m[\tilde{\varphi}]d\tilde{s}=-\int_{\gamma}m[\tilde{\varphi}]\sqrt{\tilde{g}_{ab}\tilde{u}^{a}\tilde{u}^{b}}d\tilde{\tau},$
(22)

where $\tilde{\tau}$ is the proper time in the effective spacetime and
$\tilde{u}^{\alpha}=dz^{\alpha}/d\tilde{\tau}$. Equation (22) represents our
ansatz for the point particle action. For $\tilde{\varphi}=0$ it reduces to
the point particle action in GR, Eq. (90), which we show to be consistent with
self-force calculations up to second-order [31, 60, 61] in App. A.

By varying the action (22) with respect to the effective metric,
$\tilde{g}_{ab}$, yields

$\displaystyle
T_{\mathrm{m}}^{ab}=8\pi\int_{\gamma}m[\tilde{\varphi}]\frac{\delta^{4}[x^{\mu}-z_{\mathrm{p}}^{\mu}[\tilde{\tau}]]}{\sqrt{-\tilde{g}}}\tilde{u}^{a}\tilde{u}^{b}d\tilde{\tau},$
(23)

varying with respect to the effective scalar field, $\tilde{\varphi}$, yields

$\displaystyle\Sigma=16\pi\int_{\gamma}m^{\prime}[\tilde{\varphi}]\frac{\delta^{4}[x^{\mu}-z_{\mathrm{p}}^{\mu}[\tilde{\tau}]]}{\sqrt{-\tilde{g}}}d\tilde{\tau},$
(24)

which replace Eqs. (23) and (24) respectively. Eqs. (23) and (24) inform us
that $m[\tilde{\varphi}]$ and $m^{\prime}[\tilde{\varphi}]$ are the Bondi mass
and scalar charge of the secondary compact object respectively.

We remark that our approach is motivated by Ref. [62], which applied the
effective metric approach to the point-particle stress-energy of a massive
compact object in GR. Their main result was a conjecture for the second-order
stress-energy tensor of a compact object in self-force in GR. This conjecture
was later proven to hold in regular gauges by Ref. [31].

###  III.2 The equation of motion for the secondary

Using our point particle action, Eq. (22), we can also derive the equation of
motion of the compact object. We follow the approach of [54], varying the
whole action (2) with respect to the body’s path, $x^{\mu}\rightarrow
x^{\mu}+\delta x^{\mu}$. This yields:

$\displaystyle\delta_{x^{\mu}}S_{\mathrm{p}}$
$\displaystyle=\int_{\gamma}\bigg{[}-(\tilde{g}^{\alpha}_{\
\mu}+\tilde{u}^{\alpha}\tilde{u}_{\mu})\frac{\partial
m[\varphi]}{\partial\varphi}\frac{\partial\varphi}{\partial x^{\alpha}}\delta
x^{\mu}$ $\displaystyle+\delta
x^{\mu}m[\varphi]\Big{(}\tilde{\Gamma}_{\mu\alpha\nu}\tilde{u}^{\alpha}\tilde{u}^{\nu}-\tilde{g}_{\mu\nu}\frac{d^{2}x^{\nu}}{d\tilde{\tau}^{2}}\Big{)}\bigg{]}d\tilde{\tau}\
+\mathcal{O}(\varepsilon^{3})\ .$ (25)

Note that the contribution coming from the non-minimal action is at least
order $\mathcal{O}(\varepsilon^{3})$ \- based on our earlier assumptions that
$\zeta=\mathcal{O}(\varepsilon^{2})$ and $S_{\mathrm{c}}$ contains at least
one copy of $\varphi$, which is itself order $\varepsilon$.

Requiring stationarity under first-order variations yields

$\displaystyle
m[\tilde{\varphi}]\tilde{a}^{a}=m^{\prime}[\tilde{\varphi}](\tilde{g}^{ab}+\tilde{u}^{a}\tilde{u}^{b})\partial_{b}\tilde{\varphi}+\mathcal{O}(\varepsilon^{3})\
,$ (26)

where $\tilde{a}^{a}=\tilde{u}^{b}\tilde{\nabla}_{b}\tilde{u}^{a}$ and
$\tilde{\nabla}_{b}$ is the covariant derivative of the effective metric. Eq.
(26) is equivalent to the standard self-force equation of motion for a point
scalar charge [63, 64, 65] but extended to at least second-order. The charge
moves as a point-particle being pushed away from geodesic motion in the
effective spacetime by a self-force generated by the effective scalar field.

We can also derive evolution equations for the mass and scalar charge of the
secondary compact object [54]:

$\displaystyle\frac{dm}{d\tau}=\frac{\partial
m}{\partial\tilde{\varphi}}\frac{\partial\tilde{\varphi}}{\partial\tau}=m^{\prime}[\tilde{\varphi}]u^{a}\nabla_{a}\tilde{\varphi},$
(27) $\displaystyle\frac{dm^{\prime}}{d\tau}=\frac{\partial
m^{\prime}}{\partial\tilde{\varphi}}\frac{\partial\tilde{\varphi}}{\partial\tau}=m^{\prime\prime}[\tilde{\varphi}]u^{a}\nabla_{a}\tilde{\varphi}.$
(28)

Note, Eqns. (23)-(24), and (26) reduce to the GR limit, Eq. (91), when
$\tilde{\varphi}\rightarrow 0$. The validity of the equations of motion for
black holes and self-gravitating extended compact objects in an effective
spacetime has been assessed in GR up to the second order in the self-force
expansion [55, 58, 59, 60, 61]. Our results here are analogous but in scalar-
tensor theories of gravity, up to our assumptions.

##  IV Perturbed equations

To isolate first- and second-order contributions of the field’s equations (4)
in orders of $\varepsilon$ we need to expand the non-linear behaviour of the
Einstein tensor, $G_{ab}[\mathbf{g}_{ab}]$. For a generic tensor
$g^{(0)}_{ab}+x_{ab}$, the expansion can be expressed as

$\displaystyle G_{ab}[x_{ab}]=\delta G_{ab}[x_{ab}]$
$\displaystyle+\delta^{2}G_{ab}[x_{ab},x_{ab}]$
$\displaystyle+\delta^{3}G_{ab}[x_{ab},x_{ab},x_{ab}]+\ldots$ (29)

where

$\delta^{n}G=\frac{1}{n!}\frac{d^{n}}{d\lambda^{n}}G\left[g^{(0)}_{ab}+\lambda
x_{ab}\right]\big{|}_{\lambda=0}\ .$ (30)

Replacing the expansion (15) for the metric into Eq. (IV) gives:

$\displaystyle G_{ab}[\mathbf{g}_{ab}]$ $\displaystyle=\varepsilon\delta
G_{ab}[h^{(1)}_{ab}]$ $\displaystyle\ \ \ +\varepsilon^{2}\Big{(}\delta
G_{ab}[h^{(2)}_{ab}]+\delta^{2}G_{ab}[h^{(1)}_{ab},h^{(1)}_{ab}]\Big{)}+\mathcal{O}(\varepsilon^{3}),$
$\displaystyle=\varepsilon
G^{(1)}_{ab}[\mathbf{g}_{ab}]+\varepsilon^{2}G^{(2)}_{ab}[\mathbf{g}_{ab}]+\mathcal{O}(\varepsilon^{3})$
(31)

where we have used $G_{ab}[g^{(0)}_{ab}]=0$.

We shall now derive the explicit form for the first- and second-order field
equations which determine the evolution of the metric and scalar
perturbations. We also derive the first- and second-order equations of motion
of the compact body, which depend on said perturbations. Deriving these
equations requires expanding all the quantities defined in the metric
($\mathbf{g}_{ab}$) and effective metric ($\tilde{g}_{ab}$) as perturbative
expansions around the background metric $g^{(0)}_{ab}$. The required
expansions are Eqs. (17), (16), (21), and (26), plus some non-trivial
expansions given in App. B.

###  IV.1 First order perturbations

Expanding Eq.(4) to linear order in the mass-ratio, the field equation for
$h_{ab}^{(1)}$ are given by

$G^{(1)}_{ab}=\delta
G_{ab}[h^{(1)}_{ab}]=8\pi\mu\int_{\gamma}\frac{\delta^{4}[x^{\mu}-z_{\mathrm{p}}^{\mu}[\tau]]}{\sqrt{-g}}u_{a}u_{b}d\tau\
,$ (32)

which are the same as in GR333In order for Eq (32) to satisfy the Bianchi
identity $\delta G_{ab}$ must be gauge fixed (e.g., to the Lorenz gauge) in
the self-consistent approximation [66]. This would similarly impose
restrictions on our action (2). Alternatively, one can adopt the Gralla–Wald
approximation [67], expanding the worldline in powers of $\varepsilon$, but
this approach does not admit a multi-scale expansion.. Hence, the standard
methods of metric reconstruction [68, 69, 27, 28, 70, 29] hold. Note that Eq.
(32) can be derived using either (11) or (23) for $T^{\mathrm{m}}_{ab}$
because this is linear order; that is, the result is equivalent to the metric
field equation in Refs. [22, 23, 71].

The wave equation for the scalar field perturbation $\varphi^{(1)}$ can be
obtained by expanding Eq. (7) to linear order,

$\displaystyle\Box\varphi^{(1)}$
$\displaystyle=16\pi\int_{\gamma}m_{[1]}\frac{\delta^{(4)}[x^{\mu}-z_{\mathrm{p}}^{\mu}[\tau]]}{\sqrt{-g}}d\tau\
,$ (33)

where $\Box=g^{ab}_{(0)}\nabla_{a}\nabla_{b}$ and $\nabla_{a}$ is the
covariant derivative associated with the background metric $g^{(0)}_{ab}$.
Similarly, using either (12) or (24) for $\Sigma$ is viable at this order.

Eqs. (32) and (33) are equivalent to the decoupled first-order field
equations. Hence, we can use the regular and singular field split for
$h^{(1)}_{ab}$ and $\varphi^{(1)}$ developed in Ref. [53].

As shown in [22, 23, 24], at first-order in $\varepsilon$ the mass parameter
$m_{[1]}$ can be interpreted in terms of the scalar charge of the secondary.
Expanding the scalar field in a buffer region inside the world-tube containing
the compact object, such that $\mu\ll\hat{r}\ll M$, with $\hat{r}$ being the
distance from the worldline, gives

$\displaystyle\varphi^{(1)}=\frac{\mu
d}{\hat{r}}+\mathcal{O}\Big{(}\frac{\mu^{2}}{\hat{r}^{2}}\Big{)},$ (34)

where $d$ is the dimensionless scalar charge. Replacing Eq. (33) into the
field’s equation (34), and matching the solution in the buffer zone leads to

$\displaystyle m_{[1]}=-\frac{\mu d}{4}\ .$ (35)

Equations (32)-(33) have been numerically integrated in [22, 24] for binaries
on circular and eccentric orbits, to compute the leading order dissipative
correction to the energy and the angular momentum gravitational wave fluxes.
These results are equivalent the dissipative piece of the first order self
force and are sufficient for evolving the binary dynamics to adiabatic
accuracy.

The solution for $h^{(1)}_{ab}$ and $\varphi^{(1)}$ can be used to compute the
full (dissipative plus conservative) first-order self-force using the first-
order equation of motion. Expanding Eq. (26) gives the first-order equations
of motion444Note, Eq. (36) cannot be derived directly from Eq. (10) because of
the lack of regularisation; that is, the effective action, Eq. (22), is
required. The expansion also requires Eq. (102)

$\displaystyle a^{a}_{(1)}=a^{a}_{(1)\mathrm{grav}}+a^{a}_{(1)\mathrm{scal}}\
,$ (36)

where the gravitational and the scalar components are given by:

$a^{a}_{(1)\mathrm{grav}}=-\frac{1}{2}(g^{ab}+u^{a}u^{b})(2h^{(1)\mathcal{R}}_{bd;e}-h^{(1)\mathcal{R}}_{de;b})u^{d}u^{e}\
,$ (37)

which is identical to Eq. (92), and

$a_{(1)\mathrm{scal}}^{a}=m_{[1]}(g^{ab}+u^{a}u^{b})\nabla_{b}\varphi^{(1)}_{\mathcal{R}}\
.$ (38)

Calculations for both $a^{(1)}_{\mathrm{grav}}$ and $a^{a}_{(1)\mathrm{scal}}$
have been carried out for generic orbits in the Kerr background in the
literature [72, 73, 74, 75]. These results can be exploited to derive the full
first-order self-force in scalar-tensor theories, which includes conservative
corrections to an EMRI’s evolution.

###  IV.2 Second order perturbations

Using our expansion of Eq.(4) (including the expansions in App. B), we can
express the field’s equations for the second-order metric perturbation

$\displaystyle\delta
G_{ab}[h^{(2)}_{cd}]=-\delta^{2}G_{ab}[h^{(1)}_{cd},h^{(1)}_{cd}]+\frac{1}{2}\partial_{a}\varphi^{(1)}\partial_{b}\varphi^{(1)}$
$\displaystyle-\frac{g_{ab}}{4}\partial_{c}\varphi^{(1)}\partial^{c}\varphi^{(1)}+4\pi\int_{\gamma}\frac{\delta^{4}[x^{\mu}-z_{\mathrm{p}}^{\mu}[\tau]]}{\sqrt{-g}}\bigg{[}2m_{[1]}\varphi^{(1)}_{\mathcal{R}}u_{a}u_{b}$
$\displaystyle+\mu\big{(}4h^{\mathcal{R}(1)}_{ac}u^{c}u_{b}-u_{a}u_{b}(g^{cd}_{(0)}-u^{c}u^{d})h^{\mathcal{R}(1)}_{cd}\big{)}\bigg{]}d\tau\
.$ (39)

Note, $G^{(2)}_{ab}=\delta
G_{ab}[h^{(2)}_{cd}]+\delta^{2}G_{ab}[h^{(1)}_{cd},h^{(1)}_{cd}]$. For
$m_{[1]}=0$ and $\varphi^{(1)}_{\mathcal{R}}=0$ the right hand side of Eq.
(39) reduces to the GR form, Eq. (96).

Expanding the scalar field equation, Eq. (7) (including the expansions in App.
B), to second-order, we obtain:

$\displaystyle\Box\varphi^{(2)}$
$\displaystyle=-\frac{8\pi\alpha^{(2)}}{\sqrt{-g}}\mathcal{G}^{(0)}-h^{ab}_{(1)}\nabla_{a}\nabla_{b}\varphi^{(1)}-(\nabla^{a}h^{(1)}_{ab})\nabla^{b}\varphi^{(1)}$
$\displaystyle\ \ \
+\frac{1}{2}(\nabla^{b}h_{(1)})\nabla_{b}\varphi^{(1)}+16\pi\int_{\gamma}\Big{[}m_{[2]}\varphi^{(1)}_{R}$
$\displaystyle\ \ \
-\frac{1}{2}m_{[1]}(g^{ab}+u^{a}u^{b})h_{ab}^{R(1)}\Big{]}\frac{\delta^{(4)}[x^{\mu}-z_{\mathrm{p}}^{\mu}[\tau]]}{\sqrt{-g}}d\tau.$
(40)

where $h_{(1)}=g^{ab}_{(0)}h^{(1)}_{ab}$555Eq. (39) and (40) are both
problematic in the sense they are not well defined as distributions. This is
due to both equations containing products of quantities that are singular
$h^{(1)}_{ab}$ and $\varphi^{(1)}$ on the worldline. However, it is possible
to construct well-defined equations by using a puncture scheme or
regularisation prescription (such as Hadamard finite part [76]) following Ref.
[31]. In calculating the second-order puncture [77], one will also define the
existence of an appropriate singular and regular splits of $h^{(2)}_{ab}$ and
$\varphi^{(2)}$, which we leave for future work.. As discussed earlier, $\zeta
S_{\mathrm{c}}$ is at least $\mathcal{O}(\varepsilon^{3})$, as $\zeta$ is at
least order $\varepsilon^{2}$ and $S_{\mathrm{c}}$ contains at least one copy
of $\varphi$. Due to this, there are no contributions from $S_{\mathrm{c}}$ to
Eq. (39) and the only contribution to Eq. (40) comes from a linear coupling
between $\varphi$ and ${\cal G}$ [45]. That is

$\displaystyle T_{ab}^{\mathrm{c}}$
$\displaystyle=\mathcal{O}(\varepsilon^{3}),$ (41)
$\displaystyle\Sigma_{\mathrm{c}}$
$\displaystyle=\varepsilon^{2}\frac{8\pi\alpha^{(2)}}{\sqrt{-g}}\mathcal{G}^{(0)}+\mathcal{O}(\varepsilon^{3}).$
(42)

For the leading order contribution in Eq. (42), $\alpha$ has dimension $[{\rm
length}]^{2}$, so we label $\alpha\rightarrow\alpha^{(2)}$ as
$\alpha\mathcal{G}^{(0)}=\mathcal{O}(\varepsilon^{2})$. In a Kerr background

$\displaystyle\mathcal{G}^{(0)}=24M^{2}\big{(}(r-ia\cos[\theta])^{-6}+(r+ia\cos[\theta])^{-6}\big{)}.$
(43)

The expansion of the equation of motion, Eq. (26) (using Eq. (102)), to
second-order gives

$\displaystyle a^{a}_{(2)}=a^{a}_{(2)\mathrm{grav}}+a^{a}_{(2)\mathrm{scal}}\
,$ (44)

with

$\displaystyle
a^{a}_{(2)\mathrm{grav}}=-\frac{1}{2}\Big{[}(g^{ab}+u^{a}u^{b})(2h^{(2)\mathcal{R}}_{bd;e}-h^{(2)\mathcal{R}}_{de;b})$
$\displaystyle\ \ \ \ \ -(g^{ab}+u^{a}u^{b})h_{b(1)\mathcal{R}}^{\
c}(2h^{(1)\mathcal{R}}_{cd;e}-h^{(1)\mathcal{R}}_{de;c})\Big{]}u^{d}u^{e}\ ,$
(45)

which is identical to Eq. (A.1), and

$\displaystyle\mu a_{(2)\mathrm{scal}}^{a}$
$\displaystyle=(g_{(0)}^{ab}+u^{a}u^{b})\Big{(}m_{[1]}\nabla_{b}\varphi^{(2)}_{\mathcal{R}}$
$\displaystyle\ \ \
+2m_{[2]}\varphi^{(1)}_{\mathcal{R}}\nabla_{b}\varphi^{(1)}_{\mathcal{R}}-\frac{m_{[1]}^{2}}{\mu}\varphi^{(1)}_{\mathcal{R}}\nabla_{b}\varphi^{(1)}_{\mathcal{R}}\Big{)}$
$\displaystyle\ \ \
+m_{[1]}\big{(}h_{cd}^{\mathcal{R}}u^{c}u^{d}u^{a}u^{b}-h^{ab}_{\mathcal{R}}\big{)}\nabla_{b}\varphi^{(1)}_{\mathcal{R}}.$
(46)

Eq. (IV.2) is equivalent to the second-order self-force in GR, Eq. (A.1)666The
second-order self-force equations, Eqs. (IV.2) and (46) (and similarly
calculating $h^{(2)\mathcal{R}}_{ab}$ and $\varphi^{(2)}_{\mathcal{R}}$) may
be unnecessary if flux balance laws can be derived to extract the dissipative
piece of the second-order self-force directly from $h^{(2)}_{ab}$ and
$\varphi^{(2)}$..

The interpretation of $m_{[2]}$ is similar to $m_{[1]}$ except its
contribution to the scalar charge is suppressed by an order in $\varepsilon$
as there is further coupling to the scalar field
$\varphi^{(1)}_{\mathcal{R}}$. Examining the $m_{[2]}$ piece in Eq. (40) we
see it will also generates a term in $\varphi^{(2)}$ equivalent to Eq. (34).
Note, the other delta function term (which is coupled to
$h^{\mathcal{R}(1)}_{ab}$) in Eq. (40) will also contribute a term equivalent
to Eq. (34). Therefore,

$\displaystyle\varepsilon\mu d$ $\displaystyle=-4\bigg{(}\varepsilon m_{[1]}$
$\displaystyle\ \ \
+\varepsilon^{2}\Big{(}m_{[2]}\varphi^{(1)}_{R}-\frac{1}{2}m_{[1]}(g^{ab}+u^{a}u^{b})h_{ab}^{\mathcal{R}(1)}\Big{)}\bigg{)}.$
(47)

That is, $m_{[2]}$, which we call hereafter "charge coupling", and the
coupling to $h_{ab}^{R(1)}$, provide an $\mathcal{O}(\varepsilon)$ correction
to the scalar charge $d$. We can expand $d$,

$\displaystyle d=d^{(0)}+\varepsilon d^{(1)}+\mathcal{O}(\varepsilon^{2}).$
(48)

We re-define Eq.(35) as

$\displaystyle m_{[1]}=-\frac{\mu d^{(0)}}{4},$ (49)

and define

$\displaystyle d^{(1)}$
$\displaystyle=-\frac{4}{\mu}\Big{(}m_{[2]}\varphi^{(1)}_{R}-\frac{1}{2}m_{[1]}(g^{ab}+u^{a}u^{b})h_{ab}^{\mathcal{R}(1)}\Big{)}.$
(50)

Interestingly, the only piece in the second-order equations that does not
derive from a minimally coupled scalar field to the metric is the Gauss-Bonnet
invariant term in Eq. (40), which is stationary in Kerr spacetime. Its effect
on $\varphi^{(2)}$ in Eq. (40) is then also stationary and, hence, it will not
affect the dissipative piece of the second-order self-force. It can,
therefore, be neglected for first-post adiabatic accurate modelling. Hence,
our formalism is independent of the choice of scalar-tensor theories which
obey our assumptions.

In general (including in the GR limit), the most challenging term to solve in
Eqs. (39) is computing the mode decomposition of
$\delta^{2}G_{ab}[h^{(1)}_{ab},h^{(1)}_{ab}]$ near the worldline. This problem
derives from $\delta^{2}G_{ab}[h^{(1)}_{ab},h^{(1)}_{ab}]$ being quadratic in
the first-order metric perturbation, which is singular on $\gamma$. Ref. [78]
describes and addresses this problem in detail by splitting
$\delta^{2}G_{ab}[h^{(1)}_{ab},h^{(1)}_{ab}]$ into three pieces: (i) the
regular term
$\delta^{2}G_{ab}[h^{(1)\mathcal{R}}_{ab},h^{(1)\mathcal{R}}_{ab}]$, (ii) the
mildly singular term,
$\delta^{2}G_{ab}[h^{(1)\mathcal{S}}_{ab},h^{(1)\mathcal{R}}_{ab}]$ computed
by casting $h^{(1)\mathcal{S}}_{ab}$ and $h^{(1)\mathcal{R}}_{ab}$ as a sum of
modes, and (iii) the very singular term,
$\delta^{2}G_{ab}[h^{(1)\mathcal{S}}_{ab},h^{(1)\mathcal{S}}_{ab}]$, which can
be calculated using a 4-dimensional expression for $h^{(1)\mathcal{S}}_{ab}$
associated with the mode expansion of $h^{(1)\mathcal{S}}_{ab}$. While this
approach has been recently used to calculate the first-post adiabatic quasi-
circular inspiral of a Schwarzschild binary system [79, 80, 81], the method
remains computationally expensive and highly technical. Overcoming this issue
remains a major obstacle for self-force waveform modelling in GR.

One may have expected this problem to be exacerbated in scalar-tensor theories
of gravity as introducing additional degrees of freedom may have resulted in
$\delta^{2}G_{ab}[h^{(1)}_{ab},h^{(1)}_{ab}]$ being different, for different
theories. However, the decoupling of scales on which our approach builds
provides key simplifications. As the first-order metric perturbation,
$h^{(1)}_{ab}$, is the same in GR and in all the theories of gravity specified
by the action (2), the expression for
$\delta^{2}G_{ab}[h^{(1)}_{ab},h^{(1)}_{ab}]$ remain the same. Therefore, we
can use the same $\delta^{2}G_{ab}[h^{(1)}_{ab},h^{(1)}_{ab}]$ constructed in
GR, significantly reducing the computational burden of our formalism.

Nonetheless, additional singular terms appear at the second order, as
quadratic functions of the scalar field contribute in Eq. (39). Additionally,
mixed products between the metric and the scalar linear perturbations appear
in Eq. (40). These terms are similarly challenging to calculate as they
contain products of singular functions on the worldline. However, their mode
decomposition can be computed adopting the same method used in [78].
Remarkably, since $\varphi^{(1)}$ is independent of the specific theory of
gravity, up to an amplitude rescaling given by the scalar charge, all these
contributions are scalar-tensor theory invariant.

##  V Two-timescale expansion

The two-timescale expansion is an example of a multi-scale expansion [82]. In
the EMRI context, Ref. [83] used the two-timescale expansion to argue that
first-post adiabatic models were necessary to model inspirals accurately.
Since then, there has been growing interest in applying the two-timescale
approximation to the EMRI problem [84, 85, 86, 36]. A two-timescale expansion
was implemented to produce the first-post-adiabatic waveform models for
Schwarzschild black holes in a quasi-circular inspiral in GR [79, 80, 81].
Here, we apply the two-timescale approximation to our formalism similarly. We
show that the resemblance to the two-timescale framework in GR allows
calculations in scalar-tensor theories of gravity to require only
supplementary terms.

As discussed in Appendix A.1, EMRI dynamics allow us to identify two distinct
timescales [83, 84]: (i) the fast-timescale over which the orbital phases
evolve and (ii) the slow-timescale that dictates the change of the orbital
frequencies and of the physical parameters of the system. Due to the near
periodicity of each EMRI orbit, the time evolution on the fast-timescale is
effectively periodic. In contrast, the slow-time evolution is non-trivial but
contributes beyond the leading order (as
$\tilde{t}=\mathcal{O}(\frac{1}{\varepsilon})$) [84].

The position of the compact object at any given time can be represented by the
three orbital phases,

$\displaystyle\phi_{i}:=\\{\phi_{r},\phi_{\theta},\phi_{\phi}\\}.$ (51)

The evolution of the phases, $\phi_{i}$, can be expressed in terms of their
respective frequencies
$\Omega_{i}:=\\{\Omega_{r},\Omega_{\theta},\Omega_{\phi}\\}$,

$\displaystyle\phi_{i}=\int\Omega_{i}[\tilde{t}]dt,$ (52)

where $\Omega_{i}$ evolves on the slow-time $\tilde{t}$. We assume
$\Omega_{i}$ is approximately constant on the fast-timescale. This assumption
is valid because the background spacetime is Kerr and geodesics in Kerr are
triperiodic. An appropriate choice of fast-times are the phases, $\phi_{i}$,
as they evolve on the orbital timescale.

The frequencies, $\Omega_{i}$, can be expressed in terms of the three
constants of motion of Kerr geodesics: energy, angular momentum and the Carter
constant, $J_{i}:=\\{E,J_{z},Q\\}$ [87]. That is,

$\displaystyle\Omega_{i}:=\Omega_{i}[J_{i}].$ (53)

As the compact object does not remain on a geodesic over the course of the
inspiral, the three constants of motion evolve over the slow-timescale. Their
evolution can be computed through the self-force. In turn, the evolution of
the frequencies of motion, $\Omega_{i}$ can be constructed from the self-
force,

$\displaystyle\frac{d\Omega_{i}}{dt}$ $\displaystyle=\varepsilon
F_{\Omega_{i}}^{\\{0\\}}[\Omega_{i}]+\varepsilon^{2}F_{\Omega_{i}}^{\\{1\\}}[\Omega_{i},\delta
M_{\mu}]+{\cal O}(\varepsilon^{3}),$ (54)

where $F_{i}^{\\{n\\}}$ are $n$th-post adiabatic-order self-force
coefficients. Note the adiabatic self-force coefficients only depend on the
frequencies. In contrast, the first-post adiabatic coefficients depend on the
frequencies and the change in the system’s physical parameters $\delta
M_{\mu}$, which we explain next.

Additional physical parameters of the binary evolve on the slow-timescale of
an EMRI; one must account for their evolution to achieve first-post adiabatic
accurate models. In the GR case, the mass and spin of the primary black hole
evolve as gravitational waves pass into the primary black hole horizon. The
change in mass and spin are labelled as

$\displaystyle\delta M_{A}[\tilde{t}]:=\\{\delta M,\delta J\\}.$ (55)

We now assess whether additional physical parameters of an EMRI need to be
evolved in scalar-tensor theories of gravity. In scalar-tensor theories of
gravity, the primary black hole also absorbs scalar radiation. This scalar
radiation carries energy and angular momentum, which must be accounted for in
the evolution of the $\delta M_{A}[\tilde{t}]$. One may question if the scalar
field can carry scalar charge into the supermassive black hole. However, in
general, scalar charge does not tend to be a free parameter in scalar-tensor
theories of gravity. Instead, when present, its value is fixed in terms of the
mass and spin of the black hole by regularity conditions on the horizon (see
e.g. [47, 45, 44, 46]). Hence, we expect the scalar charge of the supermassive
black hole to evolve consistently with the evolution of the mass and angular
momentum of the black hole. The scaling arguments in Secs. II.2 and III, that
multiple orders of $\varepsilon$ suppress the scalar charge of the
supermassive black hole, still hold. As the evolution of the mass is an
$\mathcal{O}(\varepsilon)$ effect, we can neglect the evolution of the scalar
charge of the supermassive black hole as it is a higher-order effect (at least
an $\mathcal{O}(\varepsilon^{3})$ effect for $n\geq 2$).

A potential caveat of this argument may be that the expectation that the
scalar charge is fixed with respect to the mass and spin is based on the
properties of stationary black holes. It is known that relaxing the assumption
of stationarity can allow for hair formation in principle [88], but it is
reasonable to expect that the charge per unit mass of the primary introduced
by the absorption of scalar radiation will be negligible in the present
scenario. Further investigation of whether regular (near the horizon)
perturbative modes could excite further independent scalar charge degrees of
freedom within this multi-scale formalism would test this hypothesis. We can
also test this hypothesis against time-domain self-force evolutions, for which
the formalism in the main body of this paper is also consistent. Fully
nonlinear numerical investigations of a supermassive black hole absorbing
scalar radiation arising from an asymmetric binary would also be of interest
in this respect.

We next turn our attention to whether the physical characteristics of the
secondary evolve on the orbital timescale. For EMRIs in GR, the evolution of
the mass of the secondary object is a high-order effect [36]. This is due to
the length scales of the orbital dynamics and the radiation wavelength being
much larger than the scale of the secondary. This argument extends to scalar-
tensor theories of gravity, so we expect the secondary object’s parameters
$\mu$ and $m_{[1]}$ (related to $\mu$ by Eq. (49)) to remain constant
throughout the inspiral for first-post adiabatic modelling. This is not to say
that the mass ($m[\tilde{\varphi}]$) and scalar charge
($\mathrm{m}^{\prime}[\tilde{\varphi}]$) of the secondary remain constant,
they evolve via Eqs. (27) and (28) respectively. That is, their evolution is
determined solely by the evolution of the scalar field, which is accounted for
in the two-timescale formalism. Hence, no additional parameters evolve on the
slow-timescale in our formalism are required; that is, Eq. (55) holds.

The evolution of the EMRI parameters, $\delta M_{A}$, can also be constructed
from the self-force,

$\displaystyle\frac{d\delta M_{A}}{dt}$ $\displaystyle=\varepsilon
F_{A}^{\\{1\\}}[\Omega_{j}]+{\cal O}(\varepsilon^{2}).$ (56)

More precisely, the evolution of $\delta M$ and $\delta J$ can be calculated
from the first-order scalar and gravitational fluxes that pass through the
supermassive black hole horizon.

To implement the two-timescale approximation we need to re-express the first-
and second-order field equations (Eqs. (32), (33), (39), and (40)) and the
equations of motion (Eqs. (36) and (44)). In the two-timescale approximation,
the field variables are expressed as [84]

$\displaystyle h^{(1)}_{ab}$
$\displaystyle=\sum_{p,q,m}h^{(1),\omega_{p,q,m}}_{ab}[\Omega_{i},x^{i}]e^{-ik^{i}\phi_{i}},$
(57) $\displaystyle\varphi^{(1)}$
$\displaystyle=\sum_{p,q,m}\varphi^{(1),\omega_{p,q,m}}[\Omega_{i},x^{i}]e^{-ik^{i}\phi_{i}},$
(58)

and

$\displaystyle h^{(2)}_{ab}$
$\displaystyle=\sum_{p,q,m}h^{(2),\omega_{p,q,m}}_{ab}[\Omega_{i},\delta
M_{A},x^{i}]e^{-ik^{i}\phi_{i}},$ (59) $\displaystyle\varphi^{(2)}$
$\displaystyle=\sum_{p,q,m}\varphi^{(2),\omega_{p,q,m}}[\Omega_{i},\delta
M_{A},x^{i}]e^{-ik^{i}\phi_{i}},$ (60)

where $k^{i}:=\\{p,q,m\\}$, $\omega_{p,q,m}=k^{i}\Omega_{i}$, and $p$, $q$,
and $m$ are integers to sum over. Eqs. (57), (57), (59) and (60) are a
discrete Fourier series in terms of the phases of motion, whose coefficients
evolve on the slow-timescale.

The advantage of the two-timescale approximation becomes apparent when you act
with a time derivative on the field variables; for example, a time derivative
of $\varphi^{(n)}$ gives,

$\displaystyle\frac{d\varphi^{(n)}}{dt}=\sum_{p,q,m}\Bigg{(}\varphi^{(n),\omega_{p,q,m}}[\tilde{t},x^{i}]\frac{\partial\phi_{j}}{\partial
t}(-ik^{j})e^{-ik^{i}\phi_{i}}$
$\displaystyle+\Big{(}\frac{\partial\Omega_{j}}{\partial
t}\frac{\partial\varphi^{(n),\omega_{p,q,m}}}{\partial\Omega_{j}}+\frac{\partial\delta
M_{\mu}}{\partial t}\frac{\partial\varphi^{(n),\omega_{p,q,m}}}{\partial\delta
M_{\mu}}\Big{)}e^{-ik^{i}\phi_{i}}\Bigg{)}.$ (61)

Evaluating the partial derivatives in Eq. (61):
$\frac{\partial\phi_{j}}{\partial t}=\Omega_{j}$ using Eq. (52);
$\frac{\partial\Omega_{j}}{\partial t}=\varepsilon
F^{\\{0\\}}_{\Omega_{i}}[\Omega_{i}]+\mathcal{O}(\varepsilon^{2})$ from Eq.
(54); and $\frac{\partial\delta M_{\mu}}{\partial t}=\varepsilon
F_{\mu}^{\\{1\\}}[\Omega]+{\cal O}(\varepsilon^{2})$ from Eq. (56). As the
background spacetime and background scalar field are stationary, and the
perturbations are of the form in Eqs. (57)-(60), we can split all time
derivatives into an algebraic $\mathcal{O}(\varepsilon^{0})$ piece and a
differential $\mathcal{O}(\varepsilon)$ piece,

$\displaystyle\frac{\partial}{\partial t}\rightarrow-
ik^{j}\Omega_{j}+\varepsilon\frac{\partial}{\partial\tilde{t}}\ ,$ (62)

where we have defined the slow-time derivatives,
$\frac{\partial}{\partial\tilde{t}}$, such that when it acts on a field
variable

$\displaystyle\frac{\partial}{\partial\tilde{t}}:=F^{\\{0\\}}_{\Omega_{j}}[\Omega_{i}]\frac{\partial\varphi^{(n),\omega_{p,q,m}}}{\partial\Omega_{j}}+F_{A}^{\\{1\\}}[\Omega_{j}]\frac{\partial\varphi^{(n),\omega_{p,q,m}}}{\partial\delta
M_{A}}.$ (63)

Note, the slow-time derivatives in Eq. (62) contribute at one order in
$\varepsilon$ higher than the fast-time derivatives (the algebraic part of Eq.
(62)). We now use a labelling convention [89] to denote the number of slow-
time derivatives in a differential operator. Take, for a simple example, the
operator

$\displaystyle A[\varphi^{(n)}]:=\frac{\partial^{2}\varphi^{(n)}}{\partial
t^{2}}.$ (64)

As $\varphi^{(n)}$ can be expressed using Eqs. (57) and (60), we replace the
time derivatives with fast and slow-time derivatives using Eq. (62), giving

$\displaystyle A[\varphi^{(n)}]$
$\displaystyle=\Big{(}-ik^{j}\Omega_{j}+\varepsilon\frac{\partial}{\partial\tilde{t}}\Big{)}\Big{(}-ik^{j}\Omega_{j}+\varepsilon\frac{\partial}{\partial\tilde{t}}\Big{)}\varphi^{(n)}$
(65) $\displaystyle=-(k^{j}\Omega_{j})^{2}\varphi^{(n)}-2i\varepsilon
k^{j}\Omega_{j}\frac{\partial\varphi^{(n)}}{\partial\tilde{t}}+\varepsilon^{2}\frac{\partial^{2}\varphi^{(n)}}{\partial\tilde{t}^{2}}.$
(66)

$A[\varphi^{(n)}]$ can be expressed in orders of slow-time derivatives:

$\displaystyle A[\varphi^{(n)}]=A^{\langle 0\rangle}$
$\displaystyle[\varphi^{(n)}]+A^{\langle 1\rangle}[\varphi^{(n)}]+A^{\langle
2\rangle}[\varphi^{(n)}],$ (67)

where

$\displaystyle A^{\langle 0\rangle}[\varphi^{(n)}]:$
$\displaystyle=-(k^{j}\Omega_{j})^{2}\varphi^{(n)},$ (68) $\displaystyle
A^{\langle 1\rangle}[\varphi^{(n)}]:$ $\displaystyle=-2i\varepsilon
k^{j}\Omega_{j}\frac{\partial\varphi^{(n)}}{\partial\tilde{t}},$ (69)
$\displaystyle A^{\langle 2\rangle}[\varphi^{(n)}]:$
$\displaystyle=\varepsilon^{2}\frac{\partial^{2}\varphi^{(n)}}{\partial\tilde{t}^{2}},$
(70)

where the number in angular brackets denotes the number of slow-time
derivatives in the differential operator.

We also need to expand additional quantities that appear in our equations:

$\displaystyle z^{\alpha}$
$\displaystyle=z_{(0)}^{\alpha}+z_{(1)}^{\alpha}+...$ (71) $\displaystyle
u^{\alpha}$ $\displaystyle=u_{(0)}^{\alpha}+u_{(1)}^{\alpha}+...$ (72)
$\displaystyle\Big{(}\frac{d\tau}{dt}\Big{)}$
$\displaystyle=\Big{(}\frac{d\tau}{dt}\Big{)}_{(0)}+\Big{(}\frac{d\tau}{dt}\Big{)}_{(1)}+...$
(73)

with the usual multi-scale expansion definitions [85].

The operators in our field equations, Eqs. (32), (39), (33), and (40),
separate into pieces with various numbers of slow-time derivatives. Take Eq.
(32), we can re-express it in the two-timescale approximation as

$\displaystyle\delta G_{ab}^{\langle 0\rangle}[h^{(1)}_{ab}]+\delta
G_{ab}^{\langle 1\rangle}[h^{(1)}_{ab}]$
$\displaystyle=8\pi\varepsilon\int_{\gamma}\mu\Bigg{(}\Big{(}u^{a}_{(0)}u^{b}_{(0)}+\varepsilon
2u^{(a}_{(1)\perp}u^{b)}_{(0)}\Big{)}\frac{\delta^{4}[x^{\mu}-z_{\mathrm{p}}^{\mu}[\tau]]}{\sqrt{-g}}$
$\displaystyle-
u_{a}^{(0)}u_{b}^{(0)}z^{\gamma}_{(1)\perp}\nabla_{\gamma}\frac{\delta^{(4)}[x^{\mu}-z_{\mathrm{p}}^{\mu}[\tau]]}{\sqrt{-g}}\Bigg{)}d\tau_{(0)}+\mathcal{O}(\varepsilon^{3}).$
(74)

where $z^{\gamma}_{(1)\perp}=(g^{\gamma}_{\
\alpha}+u^{\gamma}_{(0)}u_{\alpha}^{(0)})z^{\alpha}_{(1)}$ and
$u^{\gamma}_{(1)\perp}=(g^{\gamma}_{\
\alpha}+u^{\gamma}_{(0)}u_{\alpha}^{(0)})u^{\alpha}_{(1)}$. Note, Eq. (32) is
a first order in $\varepsilon$ equation, but implementing the two-timescale
approximation has introduced $\varepsilon^{2}$ pieces in Eq. (V). The second
order in $\varepsilon$ pieces can be promoted to the second-order field
equation, Eq. (39), giving

$\displaystyle\delta G_{ab}^{\langle
0\rangle}[h^{(1)}_{ab}]=8\pi\int_{\gamma}\mu\frac{\delta^{4}[x^{\mu}-z_{\mathrm{p}}^{\mu}[\tau]]}{\sqrt{-g}}u^{a}_{(0)}u^{b}_{(0)}d\tau_{(0)},$
(75) $\displaystyle\delta G_{ab}^{\langle
0\rangle}[h^{(2)}_{ab}]=-\delta^{2}G_{ab}^{\langle
0\rangle}[h^{(1)}_{ab},h^{(1)}_{ab}]-\delta G_{ab}^{\langle
1\rangle}[h^{(1)}_{ab}]$ $\displaystyle+\frac{1}{2}\partial^{\langle
0\rangle}_{a}\varphi^{(1)}\partial^{\langle
0\rangle}_{b}\varphi^{(1)}-\frac{1}{4}g_{ab}(\partial^{\langle
0\rangle}\varphi^{(1)})^{2}$
$\displaystyle+4\pi\int_{\gamma}\frac{\delta^{4}[x^{\mu}-z_{\mathrm{p}}^{\mu}[\tau]]}{\sqrt{-g}}\Big{(}2m_{[1]}\varphi^{(1)}_{R}u^{(0)}_{a}u^{(0)}_{b}+$
$\displaystyle\mu\big{(}4h^{R(1)}_{ac}{(0)}^{c}{(0)}_{b}-{(0)}_{a}{(0)}_{b}(g^{cd}_{(0)}-u_{(0)}^{c}u_{(0)}^{d})h^{R(1)}_{cd}\big{)}\Big{)}d\tau_{(0)}$
$\displaystyle+16\pi\int_{\gamma}\mu\bigg{(}2u^{(0)\perp}_{(a}u^{(0)}_{b)}\frac{\delta^{(4)}[x^{\mu}-z_{\mathrm{p}}^{\mu}[\tau]]}{\sqrt{-g}}$
$\displaystyle-
u_{a}^{(0)}u_{b}^{(0)}z^{\gamma}_{(1)\perp}\nabla_{\gamma}\frac{\delta^{(4)}[x^{\mu}-z_{\mathrm{p}}^{\mu}[\tau]]}{\sqrt{-g}}\bigg{)}d\tau_{(0)},$
(76)

where we have also expanded $\delta^{2}G_{ab}[h^{(1)}_{ab}]$ and
$\partial_{\mu}$ in orders of slow-time derivatives.

As we are interested in first-post adiabatic modelling here, we have neglected
any terms that are $\mathcal{O}(\varepsilon^{3})$. For example,
$\delta^{2}G^{\langle
1\rangle}_{ab}[h^{(1)}_{ab}]=\mathcal{O}(\varepsilon^{3})$; generally
$\delta^{n}G^{\langle
m\rangle}_{ab}[h^{(i)}_{ab}]=\mathcal{O}(\varepsilon^{(ni+m)})$. Deriving
higher-order equations with this algorithm is trivial.

Applying the slow-time derivative expansion algorithm to the scalar
perturbation field equations, Eqs. (33) and (40), give

$\displaystyle\Box$ ${}^{\langle
0\rangle}\varphi^{(1)}=16\pi\int_{\gamma}m_{[1]}\frac{\delta^{(4)}[x^{\mu}-z_{\mathrm{p}}^{\mu}[\tau]]}{\sqrt{-g}}d\tau_{(0)},$
(77) $\displaystyle\Box$ ${}^{\langle
0\rangle}\varphi^{(2)}=-\frac{8\pi\alpha^{(n)}}{\sqrt{-g}}\mathcal{G}^{(0)}-\Box^{\langle
1\rangle}\varphi^{(1)}$
$\displaystyle-h^{ab}_{(1)}\nabla_{a}\nabla_{b}\varphi^{(1)}-(\nabla^{a}h^{(1)}_{ab})\nabla^{b}\varphi^{(1)}+\frac{1}{2}(\nabla^{b}h_{(1)})\nabla_{b}\varphi^{(1)}$
$\displaystyle+16\pi\Bigg{(}\int_{\gamma}m_{[2]}\varphi^{(1)}_{R}\frac{\delta^{(4)}[x^{\mu}-z_{\mathrm{p}}^{\mu}[\tau]]}{\sqrt{-g}}d\tau_{(0)}$
$\displaystyle-\frac{1}{2}\int_{\gamma}m_{[1]}\frac{\delta^{(4)}[x^{\mu}-z_{\mathrm{p}}^{\mu}[\tau]]}{\sqrt{-g}}(g^{ab}+u_{(0)}^{a}u_{(0)}^{b})h_{ab}^{R(1)}d\tau_{(0)}$
$\displaystyle-\int_{\gamma}m_{[1]}z^{\gamma}_{(1)\perp}\nabla_{\gamma}\frac{\delta^{(4)}[x^{\mu}-z_{\mathrm{p}}^{\mu}[\tau]]}{\sqrt{-g}}d\tau_{(0)}\Bigg{)}\
.$ (78)

From the field variable coefficients, $\varphi^{(1),\omega_{p,q,m}}$,
$\varphi^{(2),\omega_{p,q,m}}$, $h^{(1),\omega_{p,q,m}}_{ab}$, and
$h^{(2),\omega_{p,q,m}}_{ab}$, one can calculate the orbit averaged self-
force:

$\displaystyle F^{a}_{SF}$ $\displaystyle=\varepsilon
F_{(1)}^{a}[h^{(1),\omega_{p,q,m}}_{ab},\varphi^{(1),\omega_{p,q,m}}]$
$\displaystyle\ \ \
+\varepsilon^{2}F_{(2)}^{a}[h^{(2),\omega_{p,q,m}}_{ab},\varphi^{(2),\omega_{p,q,m}}]+\mathcal{O}(\varepsilon^{3}).$
(79)

Where we have extended Eq. (86) to include the self-force from the scalar
field. We use the arguments in Ref. [83] to split the self-force into
dissipative and conservative pieces in a post-adiabatic expansion,

$\displaystyle F^{a}_{SF}$
$\displaystyle=F_{\\{0\\}(1)\mathrm{diss}}^{a}[h^{(1),\omega_{p,q,m}}_{ab},\varphi^{(1),\omega_{p,q,m}}]$
$\displaystyle\ \ \
+F_{\\{1\\}(1)\mathrm{cons}}^{a}[h^{(1),\omega_{p,q,m}}_{ab},\varphi^{(1),\omega_{p,q,m}}]$
$\displaystyle\ \ \
+F_{\\{1\\}(2)\mathrm{diss}}^{a}[h^{(2),\omega_{p,q,m}}_{ab},\varphi^{(2),\omega_{p,q,m}}]+...$
(80)

where the numbers in the curly brackets denotes the post-adiabatic order of
the self-force contribution.

We can input Eq. (V) into Eq. (99) to calculate the phase evolution. In
practice, the phases evolve via Eq. (52); that is, the evolution of the
orbital frequencies is determined by the self-force coefficients
$F_{\Omega_{i}}^{\\{0\\}}$ and $F_{\Omega_{i}}^{\\{1\\}}$ in Eq. (54). The
self-force coefficients $F_{\Omega_{i}}^{\\{n\\}}$ and $F_{A}^{\\{1\\}}$ can
be determined from the field variable coefficients (similarly to Eq. (V)),

$\displaystyle F_{\Omega_{i}}^{\\{0\\}}[\Omega_{j}]$
$\displaystyle=F_{\Omega_{i}}^{\\{0\\}}[h^{(1),\omega_{p,q,m}}_{ab},\varphi^{(1),\omega_{p,q,m}}],$
(81) $\displaystyle F_{\Omega_{i}}^{\\{1\\}}[\Omega_{j},\delta M_{A}]$
$\displaystyle=F_{\Omega_{i}}^{\\{1\\}}[h^{(1),\omega_{p,q,m}}_{ab},\varphi^{(1),\omega_{p,q,m}},$
$\displaystyle\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
h^{(2),\omega_{p,q,m}}_{ab},\varphi^{(2),\omega_{p,q,m}}],$ (82)
$\displaystyle F_{A}^{\\{1\\}}[\Omega_{j}]$
$\displaystyle=F_{A}^{\\{1\\}}[h^{(1),\omega_{p,q,m}}_{ab},\varphi^{(1),\omega_{p,q,m}}].$
(83)

The force coefficients in Eqs. (81) and (82) can be derived from the
perturbative equations of motion, Eq. (36) and (44). However, one must account
for slow-time derivative contributions to Eq. (44). To find the additional
contribution one can examine Eq. (36) using Eqs. (62), (71), and (72), giving

$\displaystyle
a_{(2)\mathrm{slow}}^{a}=m_{[1]}^{(1)}(g^{ab}+u_{(0)}^{a}u_{(0)}^{b}+2u^{(a}_{(1)}u^{b)}_{(0)})\partial^{\langle
1\rangle}_{b}\varphi^{(1)}_{\mathcal{R}}$
$\displaystyle-\frac{1}{2}(g^{ab}+u_{(0)}^{a}u_{(0)}^{b})\big{(}2\partial^{\langle
1\rangle}_{e}h^{(1)\mathcal{R}}_{bd}-\partial^{\langle
1\rangle}_{b}h^{(1)\mathcal{R}}_{be}\big{)}u_{(0)}^{d}u_{(0)}^{e}$
$\displaystyle-\frac{1}{2}(g^{ab}+u_{(0)}^{a}u_{(0)}^{b})(2h^{(1)\mathcal{R}}_{bd;e}-h^{(1)\mathcal{R}}_{de;b})2u^{(d}_{(1)}u^{e)}_{(1)}.$
(84)

The second-order equation of motion then becomes

$\displaystyle
a^{a}_{(2)}=a^{a}_{(2)\mathrm{grav}}+a^{a}_{(2)\mathrm{scal}}+a_{(2)\mathrm{slow}}^{a}.$
(85)

From Eqs. (36) and (85) one can derive precise relations for Eqs. (81) and
(82). The force coefficients in Eq. (83) can be determined from the orbit
averaged fluxes entering the primary black hole horizon in terms of
$h^{(1),\omega_{p,q,m}}_{ab}$, and $\varphi^{(1),\omega_{p,q,m}}$.

This completes the two-timescale expansion of our first-post adiabatic
accuracy modelling scheme for non-spinning binaries in scalar-tensor theories
of gravity.

##  VI Summary and Conclusions

In this paper, we have derived a framework for computing the premiere first-
post adiabatic binary modelling scheme in theories with a massless scalar
field non-minimally coupled to gravity. As discussed in detail in Sec. II, our
main assumptions are that any interaction terms for the scalar are suppressed
by a coupling with mass-dimension-2 or higher, the scalar field is shift
symmetric, and the solutions of the theory continuously connect to those of
general relativity with a minimally coupled scalar field. We have used
perturbative scaling arguments, using the suppression of black hole scalar
charge by its mass, to isolate the adiabatic and first-post adiabatic
contributions. This allows us to ignore most of the general coupling between
the metric and scalar field. This approach is similar to Refs. [22, 23, 24]
who developed and implemented this method to adiabatic order.

We have produced an ansatz for the matter action of a compact object with a
scalar charge, which is the foundation of our formalism. The long-established
scalar charged point particle action, Eq. (10) [50], becomes problematic in
the self-force approach because of divergences in the field perturbations on
the worldline. We have combined the scalar charged point particle action with
the effective metric approach [62] to produce our ansatz for the particle
action, Eq. (22). The effective metric and effective scalar field are regular
on the worldline, making our action well defined. In the limit
$\varphi\rightarrow 0$, Eq. (22) is equivalent to a point-mass effective
action in GR, which is consistent with the first- and second-order self-force
method [58, 59, 60, 61, 31, 62]. Our ansatz and formalism could be checked by
calculating a matched asymptotic expansion [90, 91, 92, 93, 60] for a scalar
charged compact binary in scalar-tensor theories of gravity.

In Sec. IV, we derived field equations for the first- and second-order metric
and scalar perturbations, Eqs. (32), (39), (33), and (40). We also derive the
first- and second-order equations of motion for the secondary object, Eq. (36)
and (44). Our formalism also trivially extends to higher orders. Because Eq.
(32) is identical in GR, the first-order Teukolsky equation holds in our
formalism [94, 95]. It is straightforward to convert the second-order metric
perturbation field equations (Eqs. (75) or (76)) and convert them to second-
order Teukolsky equations by applying the methods used in Refs. [26, 27, 96].
Ref. [97] presents an alternative method for deriving a Teukolsky equation in
alternative theories of gravity. Also, the metric reconstruction methods
developed in GR hold in our formalism (see CCK [68, 69], GHZ [27, 28], AAB
[70], and Lorenz gauge [29] metric reconstruction).

In Sec. V, we integrate our formalism into the two-timescale approximation,
allowing for efficient first-post adiabatic calculations. Our formalism is
also consistent in the time-domain [84] and the self-consistent approximation
[98]. An additional modelling problem that needs addressing is resonances
[99]. We expect that the methods used to address resonances in GR, such as
those implemented in Ref. [34], will also be applicable to our formalism.

Remarkably, our formalism has added only one additional parameter,
$m_{[2]}^{(1)}$, to the adiabatic order formalism of Refs. [22]. Hence, it
appears that the scalar charge $d$ and $m_{[2]}^{(1)}$ capture the effects of
the scalar field to first post-adiabatic order in a very large class of
theories. Our understanding of $m_{[2]}^{(1)}$ is currently limited, and it
would be interesting to investigate how different theories of gravity generate
a non-zero $m_{[2]}^{(1)}$ and how significant this contribution is to an EMRI
model. It is conceivable that in a subset of theories, the only significant
extra parameter to the first post-adiabatic order is indeed $d$. It would also
be interesting to extend our formalism to include a mass of the scalar field
and interactions that do not respect shift symmetry. The effect of the former
has already been studied at adiabatic order in [24]. We also intend to publish
a follow-up paper which will extend our formalism to include the first-post
adiabatic effects of the spin and scalar dipole of the secondary compact
object.

The main motivation for this work is to model EMRI waveforms to first-post
adiabatic accuracy for LISA data analysis. Accuracy requirement arguments
suggest that calculating the second-order self-force to $\sim 1\%$ accuracy
will likely be sufficient for LISA data analysis [100]. If the second-order
scalar self-force (and the effect of the scalar field on the gravitational
second-order self-force) are suppressed by two orders of magnitude compared to
the gravitational self-force, then their effect may be neglected. Ref. [22]
found the adiabatic scalar self-force is $\mathcal{O}(1\%)$ of the
gravitational self-force for $d=0.3$; for smaller $d$ the scalar self-force is
further suppressed. If a similar relationship is found at first-post adiabatic
order, then the conservative piece of the first-order scalar self-force and
the dissipative piece of the second-order scalar self-force may be negligible.
Nevertheless, the adiabatic contribution (the dissipative piece of the first-
order self-force) will still be significant. That is, waveforms in scalar-
tensor theories of gravity will be significantly different from those in GR,
but it may be even easier to model binaries in scalar-tensor theories of
gravity than our formalism suggests.

Implementing our formalism will be a similar difficulty as computing first-
post adiabatic models in GR. The most challenging part of such calculations is
calculating $\delta^{2}G[h^{(1)}_{ab},h^{(1)}_{ab}]$ near the worldline. We
have shown that this calculation is identical in GR and our formalism. A
method for calculating $\delta^{2}G[h^{(1)}_{ab},h^{(1)}_{ab}]$ near the
worldline is given in Ref. [78], but it is inefficient and highly technical.
There are additional, similarly challenging to calculate, pieces in our
formalism (again containing products of divergences on the worldline):
$\frac{1}{2}\partial_{a}\varphi^{(1)}\partial_{b}\varphi^{(1)}-\frac{1}{4}g_{ab}(\partial\varphi^{(1)})^{2}$
in Eq. (39) and
$h^{ab}_{(1)}\nabla_{a}\nabla_{b}\varphi^{(1)}-(\nabla^{a}h^{(1)}_{ab})\nabla^{b}\varphi^{(1)}+\frac{1}{2}(\nabla^{b}h_{(1)})\nabla_{b}\varphi^{(1)}$
in Eq. (40). We propose again using the method in Ref. [78], or any new
methods developed to tackle the problem in GR (this is currently an active
area of research in the self-force community).

Our formalism could also be used to produce intermediate-mass-ratio inspiral
waveforms. The results in Refs. [79, 80, 81] show encouraging agreement with
Numerical Relativity and first-post adiabatic self-force waveforms in GR for
quasi-circular inspirals of Schwarzschild black holes in the mass-ratio regime
of $1:10$. These results will soon be used to help future LVK data analysis
(as LVK begins to probe deeper into the disparate mass ratio regime).
Implementing our formalism, even for the simpler case of non-spinning (or
linear in spin) black hole primary, would be an important step towards
detecting or constraining the existence of a new fundamental scalar field.

##  Appendix A Self-force formalism in GR

In this appendix, we summarise the self-force approach for binary black holes
in GR. We focus on the relevance of first-post adiabatic waveform models for
gravitational wave observations. Additionally, we show how the first- and
second-order field equations and equations of motion can be derived from an
effective action.

In black hole perturbation theory, the metric ($\mathbf{g}_{ab}$) is expressed
as an expansion in orders of a small parameter, as in Eq. (16). In the
perturbative self-force approach, the small parameter is the mass ratio,
$\varepsilon=\mu/M$. Working in GR, $g^{(0)}_{ab}$ is a solution to the
Einstein field’s equations. Here, we take $g^{(0)}_{ab}$ to be the Kerr metric
[101].

The presence of the secondary compact object produces the metric
perturbations. These perturbations impart a force on the compact object,
causing it to move off geodesic motion in the background spacetime. The so-
called self-force per unit mass ($F^{a}_{SF}$) can also be expressed as an
expansion in orders of the mass ratio:

$\displaystyle F^{a}_{\textnormal{SF}}=\varepsilon
F_{(1)}^{a}\big{[}h^{(1)}_{ab}\big{]}+\varepsilon^{2}F_{(2)}^{a}\big{[}h^{(2)}_{ab}\big{]}+\mathcal{O}(\varepsilon^{3}).$
(86)

The effect of the self-force on the motion of the inspiraling object is
described by the equation of motion,

$\displaystyle u^{b}\nabla_{b}u^{a}=a^{a}=(\varepsilon
F_{(1)}^{a}+\varepsilon^{2}F_{(2)}^{a})+\mathcal{O}(\varepsilon^{3}),$ (87)

where $u^{a}$ is the four-velocity of the compact object in the background
spacetime, $a^{a}$ is the four-acceleration, and $\nabla_{a}$ is the covariant
derivative of the background metric. Eq. (87) is given by the MiSaTaQuWa
equation [58, 59] to first order, and by the second-order equation of motion
[60, 61] to second order.

###  A.1 The self-force action approach

Here, we show that it is possible to derive the GR self-force equations of
motion and field equations, up to second-order, directly from an action. We
begin with the action,

$\displaystyle
S[\mathbf{g}_{ab},\Psi]=S_{EH}[\mathbf{g}_{ab}]+S_{\mathrm{m}}[\mathbf{g}_{ab},\Psi]\
,$ (88)

where $S_{\mathrm{m}}$ is the matter action, $\Psi$ are the matter fields, and

$\displaystyle
S_{EH}[\mathbf{g}_{ab},\varphi]=\int\frac{\sqrt{-\mathbf{g}}}{16\pi}R\ d^{4}x\
,$ (89)

is the Einstein–Hilbert action.

In the compact binary problem, spacetime is a vacuum everywhere except on the
position of the compact object (the worldline). The matter action in Eq. (88)
can be replaced by the effective point particle action

$\displaystyle S_{\mathrm{p}}=-\int_{\gamma}\mu
d\tilde{s}=-\int_{\gamma}\mu\sqrt{\tilde{g}_{ab}\tilde{u}^{a}\tilde{u}^{b}}d\tilde{\tau}.$
(90)

where $\tilde{g}_{ab}$ is the effective metric defined in Eq. (21),
$\tilde{\tau}$ is the proper time in the effective spacetime, and
$\tilde{u}^{\alpha}=dz^{\alpha}/d\tilde{\tau}$.

Varying the action (88) with respect to the body’s path, $x^{\mu}\rightarrow
x^{\mu}+\delta x^{\mu}$, results in the equation of motion

$\displaystyle\mu\tilde{a}^{a}=0,$ (91)

where $\tilde{a}^{a}=\tilde{u}^{b}\tilde{\nabla}_{b}\tilde{u}^{a}$.

Expanding Eq. (91) using the expansions Eq. (21), and those in App. B, one
recovers the MiSaTaQuWa equation,

$a^{a}_{(1)}=-\frac{1}{2}(g^{a}_{(0)c}+u^{a}u_{c})(2h^{(1)\mathcal{R}}_{d(b;e)}-h^{(1)\mathcal{R}}_{be;d})u^{b}u^{e}\
,$ (92)

and the second-order equation of motion,

$\displaystyle
a^{a}_{(2)}=-\frac{1}{2}\Big{[}(g^{a}_{(0)c}+u^{a}u_{c})(2h^{(2)\mathcal{R}}_{d(b;e)}-h^{(2)\mathcal{R}}_{be;d})$
$\displaystyle\ \ \ \ \
+(g^{a}_{(0)c}+u^{a}u_{c})h_{(1)\mathcal{R}}^{cd}(2h^{(1)\mathcal{R}}_{d(b;e)}-h^{(1)\mathcal{R}}_{be;d})\Big{]}u^{b}u^{e}\
.$ (93)

Additionally, the field equations can be derived by varying the action (88).
This process is more delicate because one must vary the
$S_{EH}[\mathbf{g}_{ab}]$ with respect to the metric ($\mathbf{g}_{ab}$) and
$S_{\mathrm{m}}[\mathbf{g}_{ab}]=S_{p}[\tilde{g}_{ab}]$ with respect to
$\tilde{g}_{ab}$. The resulting field equation is [62, 31]

$\displaystyle
G[\mathbf{g}_{ab}]=8\pi\int_{\gamma}\mu\frac{\delta^{4}[x^{\mu}-z_{\mathrm{p}}^{\mu}[\tilde{\tau}]]}{\sqrt{-\tilde{g}}}\tilde{u}^{a}\tilde{u}^{b}d\tilde{\tau}.$
(94)

Geroch’s theorem [102] tells us Eq. (94) is not formally well defined as a
partial differential equation, but it will be useful for deriving the correct
perturbative equations (to at least second-order). A perturbative expansion of
Eq. (94) using the expansions in Eq. (16), (21), and the those in App. B,
gives the first- and second order field equations:

$\displaystyle\delta G_{ab}[h^{(1)}_{ab}]$
$\displaystyle=8\pi\int_{\gamma}\mu\frac{\delta^{4}[x^{\mu}-z_{\mathrm{p}}^{\mu}[\tau]]}{\sqrt{-g}}u^{a}u^{b}d\tau,$
(95) $\displaystyle\delta G_{ab}[h^{(2)}_{ab}]$
$\displaystyle=-\delta^{2}G_{ab}[h^{(1)}_{ab},h^{(1)}_{ab}]$
$\displaystyle+4\pi\int_{\gamma}\frac{\delta^{4}[x^{\mu}-z_{\mathrm{p}}^{\mu}[\tau]]}{\sqrt{-g}}$
$\displaystyle\mu\bigg{(}4h^{\mathcal{R}(1)}_{ac}u^{c}u_{b}-u_{a}u_{b}(g^{cd}_{(0)}-u^{c}u^{d})h^{R(1)}_{cd}\bigg{)}d\tau.$
(96)

###  A.2 Why second-order self-force and first-post adiabatic accuracy

Next, we summarise which parts of the self-force expansion are required for an
accurate model. Summarising Ref. [83], we characterise an accurate model as
having a small error on the final position of the compact object over the
course of the inspiral. An inspiral is generally considered to evolve on a so-
called slow-timescale ($\tilde{t}$) related to the radiation reaction
timescale, $\tilde{t}\sim
t_{rr}\sim\frac{M}{\varepsilon}=\mathcal{O}(\frac{1}{\varepsilon})$ [84]. This
relation can be derived by considering the rate at which orbital energy is
dissipated from the system through gravitational waves. The rate of energy
dissipation scales as
$\dot{E}=\frac{dE}{dt}\sim(h^{(1)}_{ab})^{2}\sim\varepsilon^{2}$ [103]. The
orbital energy of the compact object is $E\sim\mu$. As
$t\sim\frac{E}{\dot{E}}$ and
$\frac{E}{\dot{E}}\sim\frac{\mu}{\varepsilon^{2}}\sim\frac{M}{\varepsilon}$;
therefore, $t\sim\frac{M}{\varepsilon}$.

The error in the final position ($\delta z^{\mu}$) relates to the slow-
timescale and the error in the acceleration ($\delta a^{\mu}$) [56]:

$\displaystyle\delta z^{\mu}\sim t^{2}\delta
a^{\mu}\sim\tfrac{1}{\varepsilon^{2}}\delta a^{\mu}.$ (97)

For $\delta z^{\mu}$ to be small, $\delta
a^{\mu}=\mathcal{O}(\varepsilon^{3})$ is necessary. Hence, Eq. (87) must be
calculated through second order to achieve an accurate model over an entire
inspiral. That is, the first- and second-order self-force is required.

In practice, the position of the compact object is expressed using orbital
phases, $\phi_{i}$. For generic orbits in Kerr, there are three orbital
phases, that is, $i=\\{1,2,3\\}$, or $i=\\{r,\theta,\phi\\}$. The phases obey
the expansion [83]

$\displaystyle\phi_{i}[t,\varepsilon]=\frac{1}{\varepsilon}\phi^{\\{0\\}}_{i}[t,\varepsilon]+\phi^{\\{1\\}}_{i}[t,\varepsilon]+\mathcal{O}(\varepsilon).$
(98)

For the error in Eq. (98) to be small as $\varepsilon\rightarrow 0$ the
coefficients $\phi^{\\{1\\}}_{i}[t,\varepsilon]$ and
$\phi^{\\{0\\}}_{i}[t,\varepsilon]$ must be known. The appearance of
$\frac{1}{\varepsilon}$ terms in the expansion of $\phi_{i}[t,\varepsilon]$ is
acquired from the inspiral evolving over a slow-timescale.

Ref. [83] showed that $\phi^{\\{0\\}}_{i}[t,\varepsilon]$, the so-called
adiabatic contribution, depends on the dissipative piece of the first-order
self-force. They additionally showed that $\phi^{\\{1\\}}_{i}[t,\varepsilon]$
depends on the conservative piece of the first-order self-force and the
dissipative piece of the second-order self-force:

$\displaystyle\phi^{\\{0\\}}_{i}[t,\varepsilon]$
$\displaystyle=\phi^{\\{0\\}}_{i}\big{[}F_{(1)\mathrm{diss}}^{a}[h^{(1)}_{ab}]\big{]},$
(99) $\displaystyle\phi^{\\{1\\}}_{i}[t,\varepsilon]$
$\displaystyle=\phi^{\\{1\\}}_{i}\big{[}F_{(1)\mathrm{cons}}^{a}[h^{(1)}_{ab}],F_{(2)\mathrm{diss}}^{a}[h^{(2)}_{ab}]\big{]}.$
(100)

The reasoning for the conservative self-force being suppressed by one order in
$\varepsilon$ is the conservative self-force averages out over a generic Kerr
geodesic [83]. Eqs. (98), (99), and (100) show first-post adiabatic accurate
models require the full first-order self-force and the dissipative piece of
the second-order self-force.

Ref. [83] also provides the framework for implementing a two-timescale
approximation to produce first-post adiabatic self-force binary waveforms
models in GR. Recently, first-post adiabatic waveforms have shown incredible
agreement with Numerical Relativity waveforms for quasi-circular inspirals of
Schwarzschild black holes, even in the $1:10$ mass ratio regime [79, 80, 81].
These ground-breaking results suggest that first-post adiabatic models will
play a key role in future gravitational wave science, across a mass-ratio
range much wider than expected.

##  Appendix B Necessary expansions

In Eq. (23), $\tilde{g}$ appears explicitly in $\sqrt{-\tilde{g}}$ and
implicitly in $\tilde{u}^{a}$ and $d\tilde{\tau}$. We expand each $\tilde{g}$
dependence perturbatively around $g_{ab}^{(0)}$ as follows [31],

$\displaystyle\frac{1}{\sqrt{-\tilde{g}}}$
$\displaystyle=\frac{1}{\sqrt{-g}}(1-\frac{\varepsilon}{2}g^{ab}h^{R(1)}_{ab})+\mathcal{O}(\varepsilon^{2}),$
(101) $\displaystyle\frac{d\tau}{d\tilde{\tau}}$
$\displaystyle=\frac{1}{\sqrt{1-h^{R}_{ab}u^{a}u^{b}}}$
$\displaystyle=1+\frac{\varepsilon}{2}h^{R(1)}_{ab}u^{a}u^{b}+\frac{3\varepsilon^{2}}{8}[h^{R(1)}_{ab}u^{a}u^{b}]^{2}+\mathcal{O}(\varepsilon^{3}),$
(102) $\displaystyle\frac{d\tilde{\tau}}{d\tau}$
$\displaystyle=\sqrt{1-h^{R}_{ab}u^{a}u^{b}}=1-\frac{\varepsilon}{2}h^{R(1)}_{ab}u^{a}u^{b}+\mathcal{O}(\varepsilon^{2}),$
(103)

noting777We require the expansion in Eq. (102) to an order higher than Eqs.
(101) and (103) for the expansion of Eq. (26).,
$\tilde{u}^{a}=\frac{d\tau}{d\tilde{\tau}}u^{a}$ and
$\tilde{\tau}=\frac{d\tilde{\tau}}{d\tau}d\tau$.

$T_{ab}^{\mathrm{m}}$ appears in Eq. (4) with indices down whereas in Eq. (23)
it is expressed with indices up. Ref. [31] showed that the indices of the
stress-energy tensor are raised and lowered by the effective metric (not the
background metric). That is,

$\displaystyle T^{\mathrm{m}}_{ab}$
$\displaystyle=\tilde{g}_{ac}\tilde{g}_{bd}T_{\mathrm{m}}^{ab}=\varepsilon
g^{(0)}_{ac}g^{(0)}_{bd}T_{(1)}^{ab}$
$\displaystyle+\varepsilon^{2}\left[g^{(0)}_{ac}g^{(0)}_{bd}T_{\mathrm{m}(2)}^{ab}+2h^{R(1)}_{a(c}g_{d)b}^{(0)}T^{ab}_{\mathrm{m}(1)}\right]+\mathcal{O}(\varepsilon^{3}).$
(104)

###### Acknowledgements.

AS and TS acknowledge the partial support from the STFC Consolidated Grant no.
ST/V005596/1. AS would like to thank Adam Pound for helpful discussions and
comments.

## References

  * Barausse  _et al._ [2020] E. Barausse _et al._ , Gen. Rel. Grav. 52, 81 (2020), arXiv:2001.09793 [gr-qc] .
  * Kalogera  _et al._ [2021] V. Kalogera _et al._ , (2021), arXiv:2111.06990 [gr-qc] .
  * Arun  _et al._ [2022] K. G. Arun _et al._ (LISA), Living Rev. Rel. 25, 4 (2022), arXiv:2205.01597 [gr-qc] .
  * Branchesi  _et al._ [2023] M. Branchesi _et al._ , JCAP 07, 068 (2023), arXiv:2303.15923 [gr-qc] .
  * Abbott  _et al._ [2021] R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), (2021), arXiv:2112.06861 [gr-qc] .
  * Yunes  _et al._ [2012] N. Yunes, P. Pani, and V. Cardoso, Phys. Rev. D 85, 102003 (2012), arXiv:1112.3351 [gr-qc] .
  * Barausse  _et al._ [2016] E. Barausse, N. Yunes, and K. Chamberlain, Phys. Rev. Lett. 116, 241104 (2016), arXiv:1603.04075 [gr-qc] .
  * Datta  _et al._ [2020] S. Datta, R. Brito, S. Bose, P. Pani, and S. A. Hughes, Phys. Rev. D 101, 044004 (2020), arXiv:1910.07841 [gr-qc] .
  * Annulli  _et al._ [2020] L. Annulli, V. Cardoso, and R. Vicente, Phys. Lett. B 811, 135944 (2020), arXiv:2007.03700 [astro-ph.HE] .
  * Destounis  _et al._ [2020] K. Destounis, A. G. Suvorov, and K. D. Kokkotas, Phys. Rev. D 102, 064041 (2020), arXiv:2009.00028 [gr-qc] .
  * Maggio  _et al._ [2021] E. Maggio, M. van de Meent, and P. Pani, Phys. Rev. D 104, 104026 (2021), arXiv:2106.07195 [gr-qc] .
  * Sago and Tanaka [2021] N. Sago and T. Tanaka, Phys. Rev. D 104, 064009 (2021), arXiv:2106.07123 [gr-qc] .
  * Collodel  _et al._ [2022] L. G. Collodel, D. D. Doneva, and S. S. Yazadjiev, Phys. Rev. D 105, 044036 (2022), arXiv:2108.11658 [gr-qc] .
  * Liang  _et al._ [2023] D. Liang, R. Xu, Z.-F. Mai, and L. Shao, Phys. Rev. D 107, 044053 (2023), arXiv:2212.09346 [gr-qc] .
  * Zhang  _et al._ [2023] C. Zhang, H. Guo, Y. Gong, and B. Wang, JCAP 06, 020 (2023), arXiv:2301.05915 [gr-qc] .
  * Fell  _et al._ [2023] S. D. B. Fell, L. Heisenberg, and D. Veske, (2023), arXiv:2304.14129 [gr-qc] .
  * Brito and Shah [2023] R. Brito and S. Shah, (2023), arXiv:2307.16093 [gr-qc] .
  * Amaro-Seoane  _et al._ [2015] P. Amaro-Seoane, J. R. Gair, A. Pound, S. A. Hughes, and C. F. Sopuerta, in  _Journal of Physics: Conference Series_ , Vol. 610 (IOP Publishing, 2015) p. 012002.
  * Ryan [1995] F. D. Ryan, Phys. Rev. D 52, 5707 (1995).
  * Barack and Cutler [2007] L. Barack and C. Cutler, Phys. Rev. D 75, 042003 (2007), arXiv:gr-qc/0612029 .
  * Glampedakis and Babak [2006] K. Glampedakis and S. Babak, Class. Quant. Grav. 23, 4167 (2006), arXiv:gr-qc/0510057 .
  * Maselli  _et al._ [2020] A. Maselli, N. Franchini, L. Gualtieri, and T. P. Sotiriou, Physical Review Letters 125, 141101 (2020).
  * Maselli  _et al._ [2022] A. Maselli, N. Franchini, L. Gualtieri, T. P. Sotiriou, S. Barsanti, and P. Pani, Nature Astron. 6, 464 (2022), arXiv:2106.11325 [gr-qc] .
  * Barsanti  _et al._ [2022a] S. Barsanti, N. Franchini, L. Gualtieri, A. Maselli, and T. P. Sotiriou, Phys. Rev. D 106, 044029 (2022a), arXiv:2203.05003 [gr-qc] .
  * Barsanti  _et al._ [2023] S. Barsanti, A. Maselli, T. P. Sotiriou, and L. Gualtieri, Phys. Rev. Lett. 131, 051401 (2023), arXiv:2212.03888 [gr-qc] .
  * Spiers  _et al._ [2023] A. Spiers, A. Pound, and J. Moxon, arXiv preprint arXiv:2305.19332 (2023).
  * Green  _et al._ [2020] S. R. Green, S. Hollands, and P. Zimmerman, Classical and Quantum Gravity 37, 075001 (2020).
  * Toomani  _et al._ [2021] V. Toomani, P. Zimmerman, A. Spiers, S. Hollands, A. Pound, and S. R. Green, Classical and Quantum Gravity 39, 015019 (2021).
  * Dolan  _et al._ [2022] S. R. Dolan, C. Kavanagh, and B. Wardell, Physical Review Letters 128, 151101 (2022).
  * [30] B. Leather and N. Warburton, “ Applying the effective-source approach to frequency-domain self-force calculations for eccentric orbits,” in preparation.
  * Upton and Pound [2021] S. D. Upton and A. Pound, Physical Review D 103, 124016 (2021).
  * Upton [2023] S. D. Upton,  (2023), Arxiv:2309.03778v1 .
  * Osburn and Nishimura [2022] T. Osburn and N. Nishimura, arXiv preprint arXiv:2206.07031 (2022).
  * Nasipak and Evans [2021] Z. Nasipak and C. R. Evans, Physical Review D 104, 084011 (2021).
  * Piovano  _et al._ [2020] G. A. Piovano, A. Maselli, and P. Pani, Phys. Rev. D 102, 024041 (2020), arXiv:2004.02654 [gr-qc] .
  * Mathews  _et al._ [2022] J. Mathews, A. Pound, and B. Wardell, Physical Review D 105, 084031 (2022).
  * Drummond and Hughes [2022a] L. V. Drummond and S. A. Hughes, Physical Review D 105, 124040 (2022a).
  * Drummond and Hughes [2022b] L. V. Drummond and S. A. Hughes, Physical Review D 105, 124041 (2022b).
  * Drummond  _et al._ [2023] L. V. Drummond, A. G. Hanselman, D. R. Becker, and S. A. Hughes, arXiv preprint arXiv:2305.08919 (2023).
  * Saravani and Sotiriou [2019] M. Saravani and T. P. Sotiriou, Phys. Rev. D 99, 124004 (2019), arXiv:1903.02055 [gr-qc] .
  * Hawking [1972] S. W. Hawking, Commun. Math. Phys. 25, 167 (1972).
  * Sotiriou and Faraoni [2012] T. P. Sotiriou and V. Faraoni, Phys. Rev. Lett. 108, 081103 (2012), arXiv:1109.6324 [gr-qc] .
  * Hui and Nicolis [2013] L. Hui and A. Nicolis, Phys. Rev. Lett. 110, 241104 (2013), arXiv:1202.1296 [hep-th] .
  * Sotiriou and Zhou [2014a] T. P. Sotiriou and S.-Y. Zhou, Phys. Rev. D 90, 124063 (2014a), arXiv:1408.1698 [gr-qc] .
  * Sotiriou and Zhou [2014b] T. P. Sotiriou and S.-Y. Zhou, Phys. Rev. Lett. 112, 251102 (2014b), arXiv:1312.3622 [gr-qc] .
  * Thaalba  _et al._ [2023a] F. Thaalba, G. Antoniou, and T. P. Sotiriou, Class. Quant. Grav. 40, 155002 (2023a), arXiv:2211.05099 [gr-qc] .
  * Kanti  _et al._ [1996] P. Kanti, N. E. Mavromatos, J. Rizos, K. Tamvakis, and E. Winstanley, Phys. Rev. D 54, 5049 (1996), arXiv:hep-th/9511071 .
  * Thaalba  _et al._ [2023b] F. Thaalba, M. Bezares, N. Franchini, and T. P. Sotiriou, (2023b), arXiv:2306.01695 [gr-qc] .
  * Mathisson [2010] M. Mathisson, General relativity and gravitation 42, 1011 (2010).
  * Damour and Esposito-Farese [1992] T. Damour and G. Esposito-Farese, Classical and Quantum Gravity 9, 2093 (1992).
  * Nair  _et al._ [2019] R. Nair, S. Perkins, H. O. Silva, and N. Yunes, Physical review letters 123, 191101 (2019).
  * Dirac [1938] P. A. M. Dirac, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 167, 148 (1938).
  * Detweiler and Whiting [2003] S. Detweiler and B. F. Whiting, Physical Review D 67, 024025 (2003).
  * Zimmerman [2015] P. Zimmerman, Physical Review D 92, 064051 (2015).
  * Harte [2012] A. I. Harte, Classical and quantum gravity 29, 055012 (2012).
  * Pound [2012a] A. Pound, Physical Review D 86, 084019 (2012a).
  * Harte [2015] A. I. Harte, Equations of Motion in Relativistic Gravity , 327 (2015).
  * Y Mino and Tanaka [1997] M. S. Y Mino and T. Tanaka, Physical Review D 55, 3457 (1997).
  * Quinn and Wald [1997] T. C. Quinn and R. M. Wald, Physical Review D 56, 3381 (1997).
  * Pound [2012b] A. Pound, Physical Review Letters 109, 051101 (2012b).
  * Gralla [2012] S. E. Gralla, Physical Review D 85, 124011 (2012).
  * Detweiler [2012] S. Detweiler, Physical Review D 85, 044048 (2012).
  * Quinn [2000] T. C. Quinn, Physical Review D 62, 064029 (2000).
  * Gralla [2010] S. E. Gralla, Physical Review D 81, 084060 (2010).
  * Poisson  _et al._ [2011] E. Poisson, A. Pound, and I. Vega, Relativity 14 (2011).
  * Pound [2010a] A. Pound, Physical Review D 81, 024023 (2010a).
  * Gralla and Wald [2008] S. E. Gralla and R. M. Wald, Classical and Quantum Gravity 25, 205009 (2008).
  * Chrzanowski [1975] P. L. Chrzanowski, Physical Review D 11, 2042 (1975).
  * Cohen and Kegeles [1975] J. M. Cohen and L. S. Kegeles, Physics Letters A 54, 5 (1975).
  * Aksteiner  _et al._ [2019] S. Aksteiner, L. Andersson, and T. Bäckdahl, Physical Review D 99, 044043 (2019).
  * Barsanti  _et al._ [2022b] S. Barsanti, N. Franchini, L. Gualtieri, A. Maselli, and T. P. Sotiriou, arXiv preprint arXiv:2203.05003 (2022b).
  * Meent and Shah [2015] M. V. D. Meent and A. G. Shah, Physical Review D 92, 064025 (2015).
  * Meent [2016] M. V. D. Meent, Physical Review D 94, 044034 (2016).
  * Meent [2018] M. V. D. Meent, Physical Review D 97, 104033 (2018).
  * Heffernan [2022] A. Heffernan, Physical Review D 106, 064031 (2022).
  * Hadamard [1932] J. Hadamard, 11, 243 (1932).
  * Pound and Miller [2014] A. Pound and J. Miller, Physical Review D 89, 104020 (2014).
  * Miller  _et al._ [2016] J. Miller, B. Wardell, and A. Pound, Physical Review D 94, 104018 (2016).
  * Pound  _et al._ [2020] A. Pound, B. Wardell, N. Warburton, and J. Miller, Physical Review Letters 124, 021101 (2020).
  * Wardell  _et al._ [2021] B. Wardell, A. Pound, N. Warburton, J. Miller, L. Durkan, and A. L. Tiec, arXiv preprint arXiv:2112.12265 (2021).
  * Warburton  _et al._ [2021] N. Warburton, A. Pound, B. Wardell, J. Miller, and L. Durkan, Physical review letters 127, 151102 (2021).
  * Kevorkian and Cole [1996] J. Kevorkian and J. Cole, in  _Multiple Scale and Singular Perturbation Methods_ (Springer, 1996) pp. 267–409.
  * Flanagan and Hinderer [2008] E. E. Flanagan and T. Hinderer, Physical Review D 78, 064028 (2008).
  * Miller and Pound [2020] J. Miller and A. Pound, arXiv preprint arXiv:2006.11263 (2020).
  * Pound and Wardell [2021] A. Pound and B. Wardell, arXiv preprint arXiv:2101.04592 (2021).
  * Durkan and Warburton [2022] L. Durkan and N. Warburton, arXiv preprint arXiv:2206.08179 (2022).
  * Barack and Pound [2018] L. Barack and A. Pound, Reports on Progress in Physics 82, 016904 (2018).
  * Jacobson [1999] T. Jacobson, Physical Review Letters 83, 2699 (1999).
  * Spiers [2022] A. R. C. Spiers,  _Second-order gravitational self-force in Kerr spacetime_ , Ph.D. thesis, University of Southampton (2022).
  * Kevorkian and Cole [2012] J. K. Kevorkian and J. D. Cole,  _Multiple scale and singular perturbation methods_ , Vol. 114 (Springer Science & Business Media, 2012).
  * Eckhaus [2011] W. Eckhaus,  _Asymptotic analysis of singular perturbations_ (Elsevier, 2011).
  * Kates [1981] R. E. Kates, Annals of Physics 132, 1 (1981).
  * Pound [2010b] A. Pound, Physical Review D 81, 124009 (2010b).
  * Teukolsky [1972] S. A. Teukolsky, Physical Review Letters 29, 1114 (1972).
  * Teukolsky [1973] S. A. Teukolsky, Astrophys. J. 185, 635 (1973).
  * Hussain and Zimmerman [2022] A. Hussain and A. Zimmerman, Physical Review D 106, 104018 (2022).
  * Li  _et al._ [2023] D. Li, P. Wagle, Y. Chen, and N. Yunes, Physical Review X 13, 021029 (2023).
  * Pound [2015] A. Pound, in  _Equations of Motion in Relativistic Gravity_ (Springer, 2015) pp. 399–486.
  * van de Meent [2014] M. van de Meent, Physical Review D 89, 084033 (2014).
  * Miller and Pound [2021] J. Miller and A. Pound, Phys. Rev. D 103, 064048 (2021), arXiv:2006.11263 [gr-qc] .
  * Kerr [1963] R. P. Kerr, Physical review letters 11, 237 (1963).
  * Geroch [1970] R. Geroch, Journal of Mathematical Physics 11, 437 (1970).
  * Andersson [2019] N. Andersson,  _Gravitational-Wave Astronomy: Exploring the Dark Side of the Universe_ (Oxford Graduate Texts, 2019).