# Self-Force Calculations with a Spinning Secondary

Josh Mathews  0000-0002-5477-8470 School of Mathematics & Statistics,
University College Dublin, Belfield, Dublin 4, Ireland, D04 V1W8  Adam Pound
0000-0001-9446-0638 School of Mathematical Sciences and STAG Research Centre,
University of Southampton, Southampton, United Kingdom, SO17 1BJ  Barry
Wardell  0000-0001-6176-9006 School of Mathematics & Statistics, University
College Dublin, Belfield, Dublin 4, Ireland, D04 V1W8

###### Abstract

We compute the linear metric perturbation to a Schwarzschild black hole
generated by a spinning compact object, specialising to circular equatorial
orbits with an (anti-)aligned spin vector. We derive a two-timescale expansion
of the field equations, with an attendant waveform-generation framework, that
includes all effects through first post-adiabatic order, and we use the Regge-
Wheeler-Zerilli formalism in the frequency domain to generate waveforms that
include the complete effect of the spin on the waveform phase. We perform the
calculations using expansions at fixed orbital frequency, increasing the
computational efficiency and simplifying the procedure compared to previous
approaches. Finally, we provide the first fully relativistic, first-principles
regularisation procedure for gauge invariant self-force quantities to linear
order in spin. We use this procedure to produce the first strong-field,
conservative self-force calculation including the spin of the secondary —
computing Detweiler’s redshift invariant.

##  I Introduction

The era of gravitational wave astronomy was beckoned by the first detection of
gravitational waves by the LIGO and VIRGO collaborations [1]. The
collaborations’ ground-based detectors now frequently detect signals from the
mergers of binary systems of neutron stars and stellar-mass black holes [2].
The advent of future space-based gravitational wave detectors will bring the
ability to observe gravitational waves in lower frequency bands. The European
Space Agency is currently scheduled to launch LISA, a space-based
gravitational wave detector, in 2034 [3]. Extreme mass ratio inspirals (EMRIs)
are an important source of gravitational waves for LISA, with frequencies
optimally placed in LISA’s sensitivity band [4]. EMRIs consist of a compact
‘secondary body’ (of mass $\mu$) such as a stellar mass black hole or neutron
star, orbiting a significantly more massive ‘primary body’ (of mass $M$) such
as a supermassive black hole. The system radiates energy as gravitational
waves, leading the secondary’s orbit to gradually inspiral before it
eventually collides and merges with the primary body. EMRI configurations
evolve slowly while the secondary body effectively surveys the primary’s
spacetime. This information is imprinted on the resulting gravitational
waveforms, including data in the strong gravitational field as the separation
of the bodies diminishes. Thus the scientific potential of EMRI detection is
particularly alluring and can be used as an accurate test of General
Relativity in the strong field, as well as giving insight into compact object
populations [4]. The detection of EMRI waveforms and parameter estimation of
the binary’s physical attributes will require precise waveform models.

The extreme mass ratio of EMRI systems can be exploited by modelling their
waveforms using the self-force approach of black hole perturbation theory (see
for example, [5, 6, 7]). In this approach, the primary body’s spacetime is
perturbed by the secondary and the resulting perturbed spacetime metric is
expanded in powers of the mass ratio, $\mu/M$. The Einstein field equations
can then be solved order by order in the mass ratio for the metric
perturbations, which in turn can be used to calculate EMRI inspirals and
waveforms. Calculations involving the first order in mass ratio metric
perturbations are referred to as first-order self-force (1SF) calculations,
while calculations involving the second-order metric perturbations are
referred to as second-order self-force (2SF) calculations, and so on.

Leading order approximate EMRI waveforms, known as ‘adiabatic’ waveforms,
require only partial 1SF calculations. The phase evolution of these adiabatic
waveforms is obtained from the 1SF dissipative gravitational self-force. The
remaining conservative 1SF contributes to the waveform at next-to-leading
order or first post-adiabatic order (1PA). The other contributions to 1PA
waveforms are the dissipative part of 2SF (without spin effects) and terms
related to the spin of the secondary body. These additional spin terms
comprise of the spin-curvature coupling contributions and the linear in spin,
linear in mass ratio dissipative self-force and self-torque. Detection of
weaker EMRI signals and accurate parameter estimation necessitates 1PA
waveform modelling. Practical calculations of 1PA waveforms take advantage of
the slow evolution of EMRI configurations compared to the timescale of the
orbital motion of the secondary — the two-timescale expansion [8, 5]. In this
paper we develop a two-timescale expansion including all spin effects through
1PA order, in the special case of quasicircular orbits and (anti-)aligned
spin. This expansion provides a practical method of solving the field
equations, a rapid, first-principles waveform-generation framework, and a
simple, modular way of including spin effects into 1PA waveform calculations.
When combined with recent progress on 2SF calculations [9, 10], this will
enable production of complete 1PA waveforms, including (anti-)aligned spin for
quasicircular inspirals into non-spinning (or slowly spinning [10]) black
holes.

Although the conservative linear in spin self-force only contributes to the
waveform at post-2-adiabatic order (and is thus not required for 1PA waveform
modelling), the complete leading order self-force including the leading spin
effects is interesting in its own right. Conservative self-force dynamics can
be used to inform the effective-one-body formalism [11] or ‘EOB’ for short
(see [12, 13, 14, 15, 16, 17, 18]), especially with the discovery of the
‘First Law of Binary Black Hole Mechanics’ [19] and the development of
powerful scattering techniques. As detector sensitivity improves, it will
become increasingly critical to make use of this SF information to improve EOB
models of spin effects, eccentricity, and small mass ratios. Useful
ingredients to do so are gauge invariant quantities such as Detweiler’s
redshift invariant — first introduced for non-spinning secondary bodies in
circular equatorial orbits in Schwarzschild spacetime in [20] — and the spin
precession invariant [21]. The definition of the redshift invariant was
generalised to include spin effects in [22], and in the same work a Post-
Newtonian (PN) approximation was calculated for the case of an aligned-spin
secondary in a circular equatorial orbit around a Schwarzschild primary. This
calculation was later extended to include the body’s spin-induced quadrupole
moment [23]. There have been several calculations of the spin precession
invariant [21, 24, 25, 26], though to our knowledge these are yet to be
extended to include the spin of the secondary.

In black hole perturbation theory, a small extended ‘test’ body experiences an
acceleration driven by perturbation forces induced by the body’s own multipole
moments coupling to the background spacetime. The pole-dipole approximation
for a spinning test body assumes the body can be described by its mass-
monopole and spin-dipole moments only, neglecting the quadrupole and higher
moments111Note that the term ‘multipoles’ is also used separately to refer to
spherical harmonic modes (such as ‘the low multipoles’ in completing the
metric perturbation).. Omitting perturbations to the background spacetime, the
motion of a spinning test body is governed by the Mathisson-Papetrou-Dixon
(MPD) equations [27, 28, 29] with a supplementary spin condition (SSC). The
MPD equations assume that the length scale associated with the spin of the
test body are much smaller than the length scale of the curvature of the
background spacetime — this leads us to include only linear in spin effects. A
thorough study of motion governed by the MPD equations was produced in [30,
31]. The consistency of different physical results for solutions of the MPD
equations with several common SSCs was investigated in [32, 33]. The results
in [33] showed, for example, that varying the SSC can change the expression
for the orbital frequency of a spinning body in a circular orbit at quadratic
and higher order in spin. As the secondary’s spin is proportional to the
square of its mass, the leading linear spin terms enter the metric
perturbation at second order in the mass ratio, and linearizing in spin is
necessary to consistently satisfy the equations of motion at each order in the
mass ratio. Defining the dimensionless spin parameter $\chi\equiv S/\mu^{2}$
where $S$ is the magnitude of the spin vector, $S^{2}\equiv
S^{\beta}S_{\beta}=\frac{1}{2}S_{\mu\nu}S^{\mu\nu}$, this motivates us to work
with the dimensionless spin tensor, $\tilde{S}^{\mu\nu}\equiv
S^{\mu\nu}/\mu^{2}$. If the compact secondary is a black hole, for example,
then $|\chi|\leq 1$. The spin-linearized metric perturbation has the
additional advantage that it can be computed without fixing a value of $\chi$,
making it a very efficient way of filling the parameter space of $\chi$
values.

First order gravitational self-force calculations began by examining specific
EMRI models with a non-spinning secondary body and progressed towards generic
EMRI configurations [34]. Recently, there has been progress on incorporating
the secondary’s spin in specific EMRI configurations. The effect of the spin-
curvature coupling (MPD force) was incorporated into self-force calculations
in [35] for a Schwarzschild black hole primary, and waveforms were produced
including the coupling for secondaries in eccentric equatorial motion with the
spin vector aligned to the orbital angular momentum. That work neglected the
spin contribution to the gravitational wave fluxes. The spin-curvature
coupling in a Kerr background spacetime was studied in detail in the frequency
domain in [36] with a view to informing perturbative EMRI calculations. Fluxes
for a Kerr black hole primary and an aligned spinning secondary in a circular
equatorial orbit were produced in the frequency domain in [37] and in the time
domain in [38], without truncating quadratic and higher order in spin terms. A
flux balance law was derived in [39] that holds for generic orbits in Kerr to
linear order in spin and the law was demonstrated for aligned-spin circular
equatorial orbits in Schwarzschild spacetime, truncating at linear order in
spin. Linearized-in-spin fluxes were computed for a secondary in a spin-
aligned circular equatorial orbit around a Kerr primary in [40] and used to
calculate waveforms. That work was later extended to investigate aspects of
the detectability of the spin of the secondary in [41]. Flux calculations
including spin effects for the most general EMRI configuration that we are
aware of were performed in [42], with a Kerr primary and an aligned-spinning
secondary in an eccentric equatorial orbit, without truncating quadratic in
spin and higher terms.

Of all the self-force flux calculations including spin effects, we are not
aware of any that generalise the spin of the secondary away from the
(anti-)aligned case. The numerical calculations in this work will continue
this legacy, specialising to secondaries with (anti-)aligned spin in circular
equatorial orbits around a Schwarzschild primary black hole. We calculate both
the dissipative and conservative leading self-force including the next-to-
leading-order spin effects, and derive from first principles a fully
relativistic regularization procedure for gauge invariant self-force
quantities to linear order in spin. In doing so, we provide a covariant
expansion of the corresponding singular field. We also introduce a fixed
frequency parameterisation where the (quasi)circular motion of the spinning
body can be parameterised by a spin independent frequency. In frequency domain
calculations, this has a distinct advantage in linearizing numerical results
in spin — the numerically integrated perturbation equations depend on the
frequency and if there is spin dependence in the frequency, then one must find
a way to truncate numerical results at linear order in spin. Typically this
has been done by numerically fitting the spin dependence [39, 40]
(significantly increasing the computational cost) or more recently by
linearizing the perturbation equations and solving a coupled system with an
extended source [41]. Removing the spin dependence from the frequency (and
therefore from the numerical integration) entirely avoids the issues that come
with a non-compact source, and produces results that are readily comparable
with their post-Newtonian equivalents (calculated at fixed frequency) and that
are simple to implement in the two-timescale expansion for producing
waveforms.

In this paper we adopt the metric signature $-+++$ and work with geometrized
units such that $G=c=1$. Symmetrisation of indices is denoted using
parentheses while square brackets represent anti-symmetrisation. We use
$\varepsilon$ as a formal order-counting parameter to keep track of powers of
the mass ratio. Finally, when referring to gauge invariance, we mean
invariance within a class of gauges that do not interfere with the helical
symmetry of the perturbed spacetime.

In Section II we review the equations of motion for a spinning secondary body.
In Section III we specialise the motion to circular equatorial orbits in
Schwarzschild spacetime with the secondary body’s spin vector (anti-)aligned
to its total angular momentum. In Section IV we develop the two-timescale
expansion of the field equations and equations of motion. In Section V we
review the Regge-Wheeler-Zerilli formalism for solving the perturbation
equations to obtain the leading order metric perturbation (including sub-
leading spin effects) and the gravitational energy and angular momentum fluxes
associated with a spinning body. In Section VI we review Detweiler’s redshift
invariant for a spinning secondary. In Section VII we derive a regularisation
procedure for the treatment of the singular field associated with a spinning
point particle and produce regularisation parameters for the redshift
invariant. In Section VIII we present our numerical results and some
gravitational waveforms.

##  II Equations of motion with a spinning secondary

###  II.1 Mathisson-Papapetrou-Dixon-Harte equations

The equations of motion for a generic, spinning compact body are well
established through order $\varepsilon$ [7] but only partially known at order
$\varepsilon^{2}$. The most far-reaching results are due to Harte [43], who
showed that in fully nonlinear gravity, a self-gravitating material body obeys
the same equations of motion as a test body, but a test body in an effective
metric $\hat{g}_{\alpha\beta}$. If we ignore the effects of the body’s
quadrupole and higher moments, then these test-body equations, referred to as
the Mathisson-Papapetrou-Dixon (MPD) equations [27, 28, 29], take a simple
form:

$\displaystyle\frac{\hat{D}p^{\alpha}}{d\hat{\tau}}$
$\displaystyle=-\frac{1}{2}\hat{R}^{\alpha}{}_{\beta\gamma\delta}\hat{u}^{\beta}S^{\gamma\delta},$
(1a) $\displaystyle\frac{\hat{D}S^{\gamma\delta}}{d\hat{\tau}}$
$\displaystyle=2p^{[\gamma}\hat{u}^{\delta]},$ (1b)

where $p^{\mu}$ is the body’s linear momentum,
$\hat{u}^{\mu}=dz^{\mu}/d\hat{\tau}$ is its four-velocity, $S^{\gamma\delta}$
is its spin tensor, and
$\hat{D}/d\hat{\tau}=\hat{u}^{\alpha}\hat{\nabla}_{\alpha}$. $\hat{\tau}$ and
$\hat{\nabla}$ are the proper time and covariant derivative compatible with
$\hat{g}_{\alpha\beta}$. Several decades earlier, Thorne and Hartle [44]
showed that these equations also hold for a black hole. However, the results
of both Harte and Thorne and Hartle (particularly the latter) are limited by
an incomplete characterization of $\hat{g}_{\alpha\beta}$. One of our
ancillary goals will be to partly solidify their status.

For any given $\hat{g}_{\alpha\beta}$, Eq. (1) is an underdetermined system of
10 equations for 13 unknowns: the six components of the antisymmetric tensor
$S^{\gamma\delta}$, the 4 components of $p^{\mu}$, and the three independent
components of $\hat{u}^{\nu}$. A “spin supplementary condition” is therefore
required to uniquely determine the solution. Common spin supplementary
conditions are $\hat{g}_{\alpha\beta}\hat{u}^{\alpha}S^{\beta\gamma}=0$ and
$\hat{g}_{\alpha\beta}p^{\alpha}S^{\beta\gamma}=0$. If we work to linear order
in the spin, these conditions are equivalent, and they imply
$p^{\alpha}=\mu\hat{u}^{\alpha}$. The MPD equations in $\hat{g}_{\alpha\beta}$
then reduce to

$\displaystyle\frac{\hat{D}^{2}z^{\mu}}{d\hat{\tau}^{2}}=\frac{\hat{D}\hat{u}^{\alpha}}{d\hat{\tau}}$
$\displaystyle=-\frac{\mu}{2}\hat{R}^{\alpha}{}_{\beta\gamma\delta}\hat{u}^{\beta}\tilde{S}^{\gamma\delta},$
(2a) $\displaystyle\frac{\hat{D}\tilde{S}^{\gamma\delta}}{d\hat{\tau}}$
$\displaystyle=0,$ (2b)

where we have used $S^{\mu\nu}=\mu^{2}\tilde{S}^{\mu\nu}$.

###  II.2 Effective metric, self-force, and self-torque

In self-force theory, the physical metric of the spacetime is expanded in the
form

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

where $\varepsilon$ is an order-counting parameter that counts powers of
$\mu$. In our context, $g_{\alpha\beta}$ is the metric of the central black
hole and $h_{\alpha\beta}\equiv\sum\varepsilon^{n}h^{(n)}_{\alpha\beta}$ is
the perturbation due to the presence of the secondary. Since
$S^{\alpha\beta}\sim\mu^{2}$, the secondary’s spin contributes (linearly) to
$h^{(2)}_{\alpha\beta}$. Higher moments would contribute at order
$\varepsilon^{3}$ and above.

The effective metric is given by $\hat{g}_{\alpha\beta}=g_{\alpha\beta}+h^{\rm
R}_{\alpha\beta}$, where $h^{\rm R}_{\alpha\beta}=\sum\varepsilon^{n}h^{{\rm
R}(n)}_{\alpha\beta}$ is a certain piece of the physical perturbation
$h_{\alpha\beta}$. We can expand the MPD equations (2) in powers of $h^{\rm
R}_{\alpha\beta}$ following Sec. IIIA of Ref. [45] (for example). Defining the
difference between the connections on $\hat{g}_{\alpha\beta}$ and
$g_{\alpha\beta}$ as
$C^{\alpha}{}_{\beta\gamma}\equiv\hat{\Gamma}^{\alpha}_{\beta\gamma}-\Gamma^{\alpha}_{\beta\gamma}$,
we have the standard relations

$\displaystyle C^{\alpha}{}_{\beta\gamma}$
$\displaystyle=\frac{1}{2}\hat{g}^{\alpha\delta}\left(2h^{\rm
R}_{\delta(\beta;\gamma)}-h^{\rm R}_{\beta\gamma;\delta}\right),$ (4)
$\displaystyle\hat{R}^{\alpha}{}_{\beta\gamma\delta}$
$\displaystyle=R^{\alpha}{}_{\beta\gamma\delta}+2C^{\alpha}{}_{\beta[\delta;\gamma]}+2C^{\alpha}{}_{\rho[\gamma}C^{\rho}{}_{\delta]\beta},$
(5)

where a semicolon denotes the covariant derivative compatible with
$g_{\alpha\beta}$. The proper times in the two metrics are related by
$d\hat{\tau}/d\tau=\sqrt{1-h^{\rm R}_{\alpha\beta}u^{\alpha}u^{\beta}}$, where
$u^{\alpha}=dz^{\alpha}/d\tau$. Substituting these relations into Eq. (2) and
expanding, we obtain222This corrects the analogous equations in Ref. [39].
Since those equations were the starting point for Ref. [39]’s derivation of
the flux-balance law, we have independently re-derived the balance law,
including the key intermediate result (69).

$\displaystyle\frac{D^{2}z^{\mu}}{d\tau^{2}}$
$\displaystyle=-\frac{1}{2}P^{\mu\nu}(g_{\nu}{}^{\lambda}-h^{\mathrm{R}\,\lambda}_{\nu})\left(2h_{\lambda\rho;\sigma}^{\mathrm{R}}-h_{\rho\sigma;\lambda}^{\mathrm{R}}\right)u^{\rho}u^{\sigma}$
$\displaystyle\quad-\frac{\mu}{2}R^{\mu}{}_{\alpha\beta\gamma}\left(1-\frac{1}{2}h^{\mathrm{R}}_{\rho\sigma}u^{\rho}u^{\sigma}\right)u^{\alpha}\tilde{S}^{\beta\gamma}$
$\displaystyle\quad+\frac{\mu}{2}P^{\mu\nu}(2h^{\mathrm{R}}_{\nu(\alpha;\beta)\gamma}-h^{\mathrm{R}}_{\alpha\beta;\nu\gamma})u^{\alpha}\tilde{S}^{\beta\gamma}+\mathcal{O}(\varepsilon^{3})$
$\displaystyle\equiv F^{\mu},$ (6a)
$\displaystyle\frac{D\tilde{S}^{\mu\nu}}{d\tau}$
$\displaystyle=u^{(\rho}\tilde{S}^{\sigma)[\mu}g^{\nu]\lambda}\left(2h_{\lambda\rho;\sigma}^{\mathrm{R}}-h_{\rho\sigma;\lambda}^{\mathrm{R}}\right)+\mathcal{O}(\varepsilon^{2})$
$\displaystyle\equiv N^{\mu\nu},$ (6b)

where $P^{\mu\nu}\equiv g^{\mu\nu}+u^{\mu}u^{\nu}$. The spin-independent terms
in $F^{\mu}$ are referred to as the self-force (per unit $\mu$) and
$N^{\mu\nu}$ as the self-torque (per unit $\mu^{2}$). Note that Eq. (6b) is
expanded to one order lower than Eq. (6a) because the spin itself enters into
the metric at one order higher than the trajectory $z^{\mu}$.

We can also extract an evolution equation for the scaled spin vector
$\tilde{S}^{\mu}=-\frac{1}{2}\epsilon^{\mu}{}_{\alpha\beta\gamma}u^{\alpha}\tilde{S}^{\beta\gamma}$.
Substituting
$\tilde{S}^{\mu\nu}=-\epsilon^{\mu\nu}{}_{\alpha\beta}\tilde{S}^{\alpha}u^{\beta}$
into Eq. (6b) and contracting the equation with $P^{\alpha}_{\ \beta}$ (to
project out components tangent to $u^{\alpha}$), we find

$P^{\alpha}_{\ \beta}\frac{D\tilde{S}^{\beta}}{d\tau}=N^{\alpha},$ (7)

where
$N^{\mu}\equiv-\frac{1}{2}\epsilon^{\mu}{}_{\nu\rho\sigma}u^{\nu}N^{\rho\sigma}$.

The secondary’s spin contributes to the above equations in three ways: through
the standard MPD terms, which are independent of $h^{\rm R}_{\mu\nu}$; through
terms of the form $h^{\rm R}\cdot S$, which can be considered as spin-induced
self-forces and self-torques; and by contributing to $h^{{\rm R}(2)}_{\mu\nu}$
(via the spin’s contribution to $h^{(2)}_{\mu\nu}$). We can write these
contributions explicitly as

$\displaystyle F^{\mu}_{(\chi)}$
$\displaystyle=-\frac{\varepsilon^{2}}{2}P^{\mu\nu}\left(2h_{\nu\rho;\sigma}^{\mathrm{R}(\chi)}-h_{\rho\sigma;\nu}^{\mathrm{R}(\chi)}\right)u^{\rho}u^{\sigma}$
$\displaystyle\quad-\frac{\mu}{2}R^{\mu}{}_{\alpha\beta\gamma}\left(1-\frac{\varepsilon}{2}h^{\mathrm{R}(1)}_{\rho\sigma}u^{\rho}u^{\sigma}\right)u^{\alpha}\tilde{S}^{\beta\gamma}$
$\displaystyle\quad+\frac{\mu\varepsilon}{2}P^{\mu\nu}(2h^{\mathrm{R}(1)}_{\nu(\alpha;\beta)\gamma}-h^{\mathrm{R}(1)}_{\alpha\beta;\nu\gamma})u^{\alpha}\tilde{S}^{\beta\gamma}+\mathcal{O}(\varepsilon^{3}),$
(8a) $\displaystyle N^{\mu}_{(\chi)}$ $\displaystyle=\varepsilon
u^{(\rho}\tilde{S}^{\sigma)[\mu}g^{\nu]\lambda}\left(2h_{\lambda\rho;\sigma}^{\mathrm{R}(1)}-h_{\rho\sigma;\lambda}^{\mathrm{R}(1)}\right)+\mathcal{O}(\varepsilon^{2}),$
(8b)

where we introduce $h^{{\rm R}(\chi)}_{\mu\nu}$ as the linear-in-spin piece of
$h^{{\rm R}(2)}_{\mu\nu}$.

However, we note that the system of equations (6) is incomplete for two
reasons: it omits terms of the same order as it keeps, specifically test-body
spin-squared and quadrupole terms that first appear at
$\mathcal{O}(\varepsilon^{2})$; and we have not defined $h^{\rm R}_{\mu\nu}$.
Since we restrict our attention in this paper to linear spin effects, we
freely skip over the first omission. The definition of $h^{\rm R}_{\mu\nu}$ is
thornier. Harte defines a class of effective metrics in which the MPD
equations hold, but the specific choice he makes becomes singular on $z^{\mu}$
at order $\varepsilon^{2}$ in a perturbative expansion.333This can be deduced
from Eqs. (82)–(84) in Ref. [43], which express $\hat{g}_{\mu\nu}$ as the
solution to
$\hat{g}_{\mu\nu}=\mathbf{g}_{\mu\nu}-(\hat{g}_{\mu\rho}\hat{g}_{\nu\sigma}-\tfrac{1}{2}\hat{g}_{\mu\nu}\hat{g}_{\rho\sigma})H^{\rho\sigma}$,
where $H^{\rho\sigma}$ is a solution to a linear differential equation in the
spacetime $\hat{g}_{\mu\nu}$. For $\hat{g}_{\mu\nu}$ to be regular on the
particle’s worldline, $H^{\mu\nu}$ must cancel any singularities in
$\mathbf{g}_{\mu\nu}$. This is impossible beyond linear order because if we
substitute $\hat{g}_{\mu\nu}=g_{\mu\nu}+\varepsilon
h^{\mathrm{R}(1)}_{\mu\nu}+\varepsilon^{2}\hat{g}^{(2)}_{\mu\nu}+\mathcal{O}(\varepsilon^{3})$
and Eq. (3), we find that $\hat{g}^{(2)}_{\mu\nu}$ cannot contain terms any
more singular than $\sim h^{(1)}_{\mu\nu}h^{\mathrm{R}(1)}_{\rho\sigma}$,
while $h^{(2)}_{\mu\nu}$ contains terms as singular as
$h^{(1)}_{\mu\nu}h^{(1)}_{\rho\sigma}$. Thorne and Hartle effectively define
$\hat{g}_{\mu\nu}$ as the metric that exerts tidal fields on the body, leaving
it open how to determine the body’s own contribution to those fields. In Ref.
[46], one of us (AP) defined a $\hat{g}_{\mu\nu}$ that is a smooth vacuum
metric on $z^{\mu}$ and is well defined for both material bodies and black
holes. For a nonspinning body, the equation of motion (6a) is valid for AP’s
definition of $\hat{g}_{\mu\nu}$ [47, 48]. But in the case of a spinning body,
it has not been shown that the $O(\varepsilon^{2})$ and $O(\varepsilon)$ terms
in Eqs. (6a) and (6b) (respectively) are correct for AP’s $\hat{g}_{\mu\nu}$.

To skirt this issue, we note that Harte’s definition is well defined and
smooth for the effects we focus on: linear effects, whether linear in $\mu$ or
linear in spin. For those contributions to $h^{{\rm R}(1)}_{\mu\nu}$ and
$h^{{\rm R}(2)}_{\mu\nu}$, Harte’s definition reduces to the more familiar
Detweiler-Whiting definition [49].444The only quantity in (6) that is not
covered by this definition is the spin-independent piece of $h^{{\rm
R}(2)}_{\mu\nu}$, for which one can use AP’s definition. We can therefore
assume with some confidence that the linear-in-spin terms in (6) are correct
with this definition of $h^{{\rm R}}_{\mu\nu}$.

We return to the regular field in Sec. VII. There, as a byproduct of our
concrete calculations, we show that Harte’s definition of $h^{{\rm
R}(\chi)}_{\mu\nu}$ agrees with AP’s, giving us additional confidence in our
assumption.

###  II.3 Stress-energy of a spinning body

Since we are only interested in the gravitational field on scales much larger
than the body’s size, we can replace the body with a “gravitational skeleton”
[27]—a point singularity equipped with the body’s multipole moments.
Concretely, through second order in $\varepsilon$, a generic compact object
can be modelled as a spinning point particle in $\hat{g}_{\mu\nu}$, with a
stress-energy tensor

$T^{\alpha\beta}=\varepsilon
T_{(\mu)}^{\alpha\beta}+\varepsilon^{2}T_{(\chi)}^{\alpha\beta}+\mathcal{O}(\varepsilon^{3}),$
(9)

where $T_{(\mu)}^{\alpha\beta}$ is a mass-monopole term and
$T_{(\chi)}^{\alpha\beta}$ is a spin-dipole term. Explicitly, the two
contributions are

$\displaystyle T_{(\mu)}^{\alpha\beta}(x)$ $\displaystyle=\mu\int
d\hat{\tau}^{\prime}\,\frac{\delta^{4}\left[x^{\mu}-z^{\mu}(\tau^{\prime})\right]}{\sqrt{-\hat{g}^{\prime}}}\hat{u}^{\alpha}(\tau^{\prime})\hat{u}^{\beta}(\tau^{\prime}),$
(10a) $\displaystyle T_{(\chi)}^{\alpha\beta}(x)$
$\displaystyle=\nabla_{\rho}\left[\int
d\tau^{\prime}\,\frac{\delta^{4}\left[x^{\mu}-z^{\mu}(\tau^{\prime})\right]}{\sqrt{-g^{\prime}}}u^{(\alpha}(\tau^{\prime}){S}^{\beta)\rho}(\tau^{\prime})\right],$
(10b)

where $\delta^{4}$ is the four-dimensional Dirac delta function, and the
covariant derivative $\nabla_{\rho}$ is with respect to the arbitrary field
point $x^{\mu}$ and not the worldline point $z^{\mu}$. Note that both
$T_{(\mu)}^{\alpha\beta}$ and $T_{(\chi)}^{\alpha\beta}$ have subleading
dependence on $\chi$ (and on $h^{\rm R}_{\alpha\beta}$) via their dependence
on the worldline. However, the quantity $h^{\mathrm{R}(\chi)}_{\alpha\beta}$
defined in the previous section corresponds only to the regular field
associated with $T_{(\chi)}^{\alpha\beta}$. Also note that no hats appear on
quantities in $T_{(\chi)}^{\alpha\beta}$ since the difference would only
contribute at order $\varepsilon^{3}$.

Unlike the equations of motion (6), the stress-energy tensor (9) has been
rigorously derived from the method of matched asymptotic expansions [46, 50].
It holds for black holes, material bodies, and exotic compact objects.

##  III Circular orbits and aligned spins: test-spin effects

Before considering the full system of equations (6), we consider the orbit and
spin with $h_{\mu\nu}$ set to zero. The equations of motion are then given by
the MPD equations Eq. (2) with all hats removed; the orbit is accelerated in
$g_{\mu\nu}$ by the MPD spin force on the right-hand side of Eq. (2a), and the
spin is parallel-propagated in $g_{\mu\nu}$. The results for this case will
carry over directly to the full problem.

###  III.1 Fixed-frequency parametrisation

We specialise to the case of a non-spinning black hole primary, in which case
the background is the Schwarzschild spacetime with line element

$ds^{2}=-fdt^{2}+f^{-1}dr^{2}+r^{2}d\Omega^{2},$ (11)

where $f\equiv 1-\frac{2M}{r}$ and $d\Omega^{2}\equiv
d\theta^{2}+\sin^{2}\theta d\phi^{2}$. The spacetime admits two Killing
vectors, $\xi_{(t)}^{\mu}=\partial_{t}$ and
$\xi_{(\phi)}^{\mu}=\partial_{\phi}$, which in turn give rise to two constants
of motion preserved by the MPD equations. In terms of a generic Killing vector
$\xi^{\mu}$, the corresponding constant of motion is given by

$\Xi=u^{\alpha}\xi_{\alpha}+\frac{\mu}{2}\tilde{S}^{\alpha\beta}\nabla_{\alpha}\xi_{\beta}.$
(12)

If $\xi^{\mu}=\xi_{(t)}^{\mu}$, then $\Xi$ is (minus) the particle’s specific
energy; if $\xi^{\mu}=\xi_{(\phi)}^{\mu}$, then it is the particle’s angular
momentum.

When the secondary’s spin vector is aligned (or anti-aligned) with its orbital
angular momentum, the MPD equations admit a solution for circular equatorial
orbits. Specialising to that case, the orbit is described by

$r=r_{p},\quad\theta=\frac{\pi}{2},\quad\phi=\Omega t,$ (13)

where $\Omega=\frac{d\phi}{dt}$ is the constant orbital frequency and where
the radial and polar components of the secondary’s four-velocity vanish,
$u^{r}=0=u^{\theta}$.

For a non-spinning test body in circular, equatorial geodesic motion, the
orbital radius and frequency are related by

$r_{p}=r_{0},\quad{\Omega}=\Omega_{0}\equiv\sqrt{\frac{M}{r_{0}^{3}}}.$ (14)

The MPD spin force accelerates the orbit, altering this relationship. To
determine the change, one can fix the orbital radius and examine how the spin
alters the orbital frequency, or one can fix the frequency and examine how the
spin alters the radius. Previous numerical self-force calculations which
included the spin of the secondary parametrised the test-body motion with a
fixed-radius parametrisation [37, 38, 39, 40, 41]

$r_{p}=r_{0},\quad\Omega(r_{p},\varepsilon)={\Omega}_{0}+\varepsilon\Omega_{\chi}(r_{p})+\mathcal{O}(\varepsilon^{2}).$
(15)

However, as illustrated in Sec. V, it is better to instead adopt a fixed-
frequency parametrisation,

$r_{p}(\Omega,\varepsilon)=r_{\Omega}+\varepsilon
r_{\chi}(\Omega)+\mathcal{O}(\varepsilon^{2}),\quad\Omega={\Omega}_{0}.$ (16)

To emphasize that $r_{0}=r_{0}(\Omega)$ in Eq. (16), we have labelled it as
$r_{\Omega}$. Throughout this paper we work exclusively with the fixed-
frequency parametrisation, allowing us to refer simply to the physical
frequency $\Omega$ without a 0 adornment. (The notable expception is Appendix
C, where we discuss the alternative fixed-radius parametrisation.)

To derive the fixed-frequency parametrisation from the equations of motion
(2a) and (2b) (with hats removed), we take advantage of the fact that we are
only looking for solutions valid through order $\varepsilon$ (i.e., linear
order in the spin force). We hence seek a solution of the form

$z^{\alpha}(t,\varepsilon)=z^{\alpha}_{0}(t)+\varepsilon
z^{\alpha}_{1}+\mathcal{O}(\varepsilon^{2}),$ (17)

where $z^{\mu}_{0}=\big{(}t,r_{\Omega},\pi/2,\Omega t\big{)}$ and
$z^{\alpha}_{1}=r_{\chi}\delta^{\alpha}_{r}$.

Substituting this ansatz into (2a) and solving order by order in
$\varepsilon$, we obtain the relations

$r_{\Omega}=\frac{M}{(M\Omega)^{2/3}},\quad
r_{\chi}=\frac{-F_{(1,\chi)}^{r}}{3\left({u}_{0}^{t}\right)^{2}f_{\Omega}\Omega^{2}},$
(18)

where $f_{\Omega}\equiv f(r_{\Omega})$, $u_{0}^{\beta}$ is the zeroth-order
four-velocity, and
$F_{(1,\chi)}^{\alpha}=-\frac{\mu}{2}R^{r}{}_{\beta\gamma\delta}u_{0}^{\beta}\tilde{S}^{\gamma\delta}\delta^{\alpha}_{r}$
is the (purely radial) leading-order spin force. The first equality in (18) is
simply a restatement of the geodesic relationship (14).

The non-zero components of the four-velocity are given by the circular-orbit
condition, $u^{\phi}=u^{t}\Omega$, along with the normalisation
$u^{\alpha}u_{\alpha}=-1$. Conveniently, this implies that when parametrised
at fixed frequency the four-velocity of the spinning secondary is equal to the
corresponding non-spinning geodesic four-velocity,
$u^{\alpha}={u}_{0}^{\alpha}+\mathcal{O}(\varepsilon^{2})$.

Writing Eq. (18) more explicitly requires an explicit form for the spin
tensor. In the aligned-spin case, the spin vector is given by
$S^{\mu}=S^{\theta}\delta^{\mu}_{\theta}$, and the corresponding scaled spin
vector $\tilde{S}^{\mu}$ has only one non-zero component,
$\tilde{S}^{\theta}=-\frac{\chi}{r_{\Omega}}$. Introducing a unit vector
$\hat{z}^{\alpha}\equiv-\frac{1}{r_{\Omega}}$ along the zeroth-order
worldline, we write

$\tilde{S}^{\mu}=\chi\hat{z}^{\mu}.$ (19)

The scaled spin tensor
$\tilde{S}^{\mu\nu}=-\epsilon^{\mu\nu\alpha\beta}\tilde{S}_{\alpha}u_{\beta}$
then has two independent non-zero components [39],

$\tilde{S}^{tr}=-\frac{\chi}{r_{\Omega}}u_{\phi}=-\tilde{S}^{rt},\quad\tilde{S}^{r\phi}=-\frac{\chi}{r_{\Omega}}u_{t}=-\tilde{S}^{\phi
r}.$ (20)

and we can evaluate the spin force to get

$F_{(1,\chi)}^{\alpha}=3\mu\chi
f_{\Omega}r_{\Omega}\Omega^{3}({u}_{0}^{t})^{2}\delta^{\alpha}_{r}.$ (21)

Given this spin force, we find that in the aligned-spin case with a fixed-
frequency parametrisation, we have motion described by

$\displaystyle z^{\alpha}$ $\displaystyle=\big{(}t,r_{\Omega}-\mu\chi\Omega
r_{\Omega},\pi/2,\Omega t\big{)}+\mathcal{O}(\varepsilon^{2}),$ (22)
$\displaystyle u^{\alpha}$
$\displaystyle=\sqrt{\frac{r_{\Omega}}{r_{\Omega}-3M}}\left(1,0,0,\Omega\right)+\mathcal{O}(\varepsilon^{2}),$
(23)

where we note that the sole spin dependence is in the radial position of the
worldline, and where the contravariant components of the four-velocity are
independent of $\chi$. Note that the covariant components $u_{\alpha}$ are not
spin independent as $\chi$ enters via the metric components evaluated on the
worldline.

Evaluating Eq. (12) for the two conserved quantities associated with the
timelike and angular Killing vectors yields the specific energy and specific
angular momentum, which are given respectively by

$\displaystyle E$ $\displaystyle=E_{0}+\varepsilon
E_{\chi}=f(r_{\Omega})u^{t}_{0}-\frac{\mu\chi}{M}u^{t}_{0}\left(\frac{M}{r_{\Omega}}\right)^{5/2},$
(24) $\displaystyle J$ $\displaystyle=J_{0}+\varepsilon
J_{\chi}=\sqrt{Mr_{\Omega}}u^{t}_{0}+\frac{\mu\chi}{M}u^{t}_{0}\left(\frac{M}{r_{\Omega}}\right)(r_{\Omega}-4M),$
(25)

omitting ${\cal O}(\varepsilon^{2})$ terms.

###  III.2 Stress-energy tensor

The stress-energy tensor for the circular-orbit, aligned-spin case can be
written in the form

$\displaystyle T_{(\mu)}^{\mu\nu}$ $\displaystyle=\frac{\mu
K_{0}^{\mu\nu}}{r^{2}\sin\theta}\delta_{r}\delta_{\theta}\delta_{\phi},$ (26a)
$\displaystyle T_{(\chi)}^{\mu\nu}$
$\displaystyle=\frac{\chi\mu^{2}}{r^{2}\sin\theta}[K_{1}^{\mu\nu}\delta_{r}\delta_{\theta}\delta_{\phi}+K_{2}^{\mu\nu}\delta_{r}\delta_{\theta}\delta_{\phi}^{\prime}+K_{3}^{\mu\nu}\delta_{r}^{\prime}\delta_{\theta}\delta_{\phi}],$
(26b)

where we use the shorthand notation $\delta_{r}\equiv\delta(r-r_{p})$,
$\delta_{\theta}\equiv\delta(\theta-\frac{\pi}{2})$ and
$\delta_{\phi}\equiv\delta(\phi-\phi_{p})$. The non-zero components of the
constant tensors $K_{0}^{\mu\nu},K_{1}^{\mu\nu},K_{2}^{\mu\nu}$ and
$K_{3}^{\mu\nu}$ are the same as those defined in Ref. [39] with the
identification $r_{p}\leftrightarrow r_{\Omega}$ (along with an additional
factor of $M^{-1}$ in our definition of $K_{1}^{\mu\nu},K_{2}^{\mu\nu}$, and
$K_{3}^{\mu\nu}$); for completeness we give their full expressions in Appendix
A.

When we substitute the expansions $r_{p}=r_{\Omega}+\varepsilon
r_{\chi}+\mathcal{O}(\varepsilon^{2})$ and
$u^{\alpha}=u^{\alpha}_{0}+\mathcal{O}(\varepsilon^{2})$ into
$T^{\mu\nu}_{(\mu)}$ and $T^{\mu\nu}_{(\chi)}$, we can make the trivial
replacements $r_{p}\to r_{\Omega}$ in (26b) (and we have preemptively done so
in the expressions in Appendix A) but we must keep the subleading term in
(26a), yielding

$T_{(\mu)}^{\mu\nu}=\frac{\mu
K_{0}^{\mu\nu}}{r^{3}_{\Omega}}\delta_{\theta}\delta_{\phi}\left(r_{\Omega}\delta_{r_{\Omega}}-2\varepsilon
r_{\chi}\delta_{r_{\Omega}}-\varepsilon
r_{\chi}\delta^{\prime}_{r_{\Omega}}\right)+\mathcal{O}(\varepsilon^{2}),$
(27)

where $\delta_{r_{\Omega}}\equiv\delta(r-r_{\Omega})$ and
$\delta^{\prime}_{r_{\Omega}}\equiv\partial_{r}\delta(r-r_{\Omega})$. The
total stress-energy (9) then becomes

$T^{\mu\nu}=\varepsilon
T^{\mu\nu}_{(1)}+\varepsilon^{2}T^{\mu\nu}_{(2)}+\mathcal{O}(\varepsilon^{3}),$
(28)

where

$\displaystyle T^{\mu\nu}_{(1)}$ $\displaystyle=\frac{\mu
K_{0}^{\mu\nu}}{r^{2}_{\Omega}\sin\theta}\delta_{r_{\Omega}}\delta_{\theta}\delta_{\phi},$
(29a) $\displaystyle T^{\mu\nu}_{(2)}$ $\displaystyle=-\frac{\mu
r_{\chi}K_{0}^{\mu\nu}}{r^{3}_{\Omega}}\delta_{\theta}\delta_{\phi}\left(2\delta_{r_{\Omega}}+\delta^{\prime}_{r_{\Omega}}\right)+T_{(\chi)}^{\mu\nu}.$
(29b)

##  IV Two-timescale expansion

When the metric perturbations $h^{(n)}_{\alpha\beta}$ are accounted for, the
binary system slowly evolves due to dissipation. In this section we show how
the linear effects of the particle’s spin can easily be incorporated into the
two-timescale evolution scheme of Ref. [8]. We closely follow the particular
formulation in Appendix A of that reference.

Our method assumes the particle’s trajectory, its spin, and the spacetime
metric only depend on $t$ through the $t$ dependence of a set of mechanical
variables $(\phi_{p},{\cal J}_{I})$. The slow evolution is captured by the
parameters ${\cal J}_{I}=(\Omega,\chi,\delta M,\delta J)$, which evolve on the
radiation-reaction timescale $\sim\Omega/\dot{\Omega}\sim 1/\varepsilon$.
(Though we find below that $\chi$ is constant at 1PA order.) Here
$(\mu\,\delta M,\mu\,\delta J)$ represent corrections to the central black
hole’s mass and spin, which evolve due to the flux of energy and angular
momentum into the black hole; we pull out an overall factor of $\mu$ to make
$(\delta M,\delta J)$ order unity. During the slow evolution, the system is
assumed to retain a periodic dependence on the particle’s orbital phase
$\phi_{p}$, which evolves on the fast timescale $\sim 1/\Omega$.

Given that the method requires a choice of time coordinate, it will be
convenient to adopt the 3+1 split $x^{\mu}=(t,x^{i})$.

###  IV.1 Orbit, spin, and metric

We first consider the particle’s orbit. In place of Eq. (17), we write the
coordinate trajectory as $z^{\mu}(t,\varepsilon)=(t,z^{i}(t,\varepsilon))$ and
assume that $z^{i}(t,\varepsilon)=z^{i}(\phi_{p}(t,\varepsilon),{\cal
J}_{I}(t,\varepsilon),\varepsilon)$. Expanding in powers of $\varepsilon$ at
fixed $(\phi_{p},{\cal J}_{I})$, we write

$z^{i}(\phi_{p},{\cal
J}_{I},\varepsilon)=z_{0}^{i}(\phi_{p},\Omega)+\varepsilon z_{1}^{i}({\cal
J}_{I})+O(\varepsilon^{2}),$ (30)

where the leading-order trajectory is

$z_{0}^{i}(\phi_{p},\Omega)=(r_{0}(\Omega),\pi/2,\phi_{p}),$ (31)

and the subleading term is a purely radial correction

$z_{1}^{i}({\cal J}_{I})=(r_{1}({\cal J}_{I}),0,0).$ (32)

We continue to define the frequency as

$\frac{d\phi_{p}}{dt}\equiv\Omega.$ (33)

The above ansatz represents an orbit that remains in the equatorial plane,
with a slowly evolving radius and frequency. Accordingly, we seek a solution
in which the spin remains orthogonal to the equatorial plane,

$\tilde{S}^{\alpha}({\cal
J}_{I})=\chi\hat{z}^{\alpha}(\Omega)+O(\varepsilon),$ (34)

in analogy with (19).

Following the same pattern, we expand the metric
as555$\varepsilon^{2}\ln\varepsilon$ terms also appear, at least in the Lorenz
and similar gauges. We hide that logarithmic dependence in $h^{(2)}_{\mu\nu}$.
As discussed in Ref. [8], it is also preferable to adopt a hyperboloidal time
coordinate $s=t-k(r^{*})$. We elide that detail here, as it is not important
for the spin contributions at 1PA order.

${\sf g}_{\mu\nu}=g_{\mu\nu}(x^{i})+\varepsilon
h^{(1)}_{\mu\nu}(x^{i},\phi_{p},\Omega,\delta M,\delta J)\\\
+\varepsilon^{2}h^{(2)}_{\mu\nu}(x^{i},\phi_{p},{\cal
J}_{I})+O(\varepsilon^{3}),$ (35)

with the assumption that each term is periodic in
$\phi_{p}=\phi_{p}(t,\varepsilon)$. For simplicity, we suppress dependence on
$M$ and $\mu$, with the understanding that $h^{(n)}_{\mu\nu}$ is equal to
$\mu^{n}$ times a $\mu$-independent function of $(x^{i},\phi_{p},{\cal
J}_{I})$.

The first-order perturbation $h^{(1)}_{\mu\nu}(x^{i},\phi_{p},\Omega,\delta
M,\delta J)$ is linear in $\delta M$ and $\delta J$, and it will be convenient
to peel off that dependence, writing

$h^{(1)}_{\mu\nu}=h^{(1)}_{\mu\nu}(x^{i},\phi_{p},\Omega,0,0)+\delta
M\,h^{(\delta M)}_{\mu\nu}(x^{i})\\\ +\delta J\,h^{(\delta
J)}_{\mu\nu}(x^{i}).$ (36)

If we replaced $\phi_{p}$ with its geodesic expression $\Omega t$, then the
first term would be the standard perturbation due to a point mass on a
circular geodesic with frequency $\Omega$. The terms $\delta M\,h^{(\delta
M)}_{\mu\nu}(x^{i})$ and $\delta J\,h^{(\delta J)}_{\mu\nu}(x^{i})$ are linear
perturbations toward a Kerr black hole with mass $M+\mu\,\delta M$ and angular
momentum $\mu\,\delta J$. We show below that at 1PA order, $\delta M$ and
$\delta J$ do not couple to the spin $\chi$.

When substituting these expansions into the equations of motion and field
equations, we apply the chain rule

$\frac{\partial}{\partial x^{\alpha}}=e^{i}_{\alpha}\frac{\partial}{\partial
x^{i}}+t_{\alpha}\left(\frac{d\phi_{p}}{dt}\frac{\partial}{\partial\phi_{p}}+\frac{d{\cal
J}_{I}}{dt}\frac{\partial}{\partial{\cal J}_{I}}\right),$ (37)

where $e^{i}_{\alpha}\equiv\frac{\partial x^{i}}{\partial x^{\alpha}}$ and
$t_{\alpha}\equiv\partial_{\alpha}t$. This implies the expansion

$\nabla_{\alpha}=\nabla^{(0)}_{\alpha}+\varepsilon
t_{\alpha}\vec{\partial}_{\cal V}+{\cal O}(\varepsilon^{2}),$ (38)

where the zeroth-order covariant derivative is

$\nabla^{(0)}_{\alpha}=e^{i}_{\alpha}\frac{\partial}{\partial
x^{i}}+t_{\alpha}\Omega\frac{\partial}{\partial\phi_{p}}+\text{Christoffel
terms}.$ (39)

${\cal V}_{I}=(F^{(0)}_{\Omega},F^{(1)}_{\chi},F^{(1)}_{\delta
M},F^{(1)}_{\delta J})$ is the leading-order velocity through parameter space,
and

$\displaystyle\vec{\partial}_{\cal V}$ $\displaystyle\equiv{\cal
V}_{I}\frac{\partial}{\partial J_{I}}$ (40)
$\displaystyle=F_{\Omega}^{(0)}\frac{\partial}{\partial\Omega}+F_{\chi}^{(1)}\frac{\partial}{\partial\chi}+F_{\delta
M}^{(1)}\frac{\partial}{\partial\delta M}+F_{\delta
J}^{(1)}\frac{\partial}{\partial\delta J}.$ (41)

$\nabla^{(0)}_{\alpha}$ acts at fixed parameter values; its action on
$h^{(1)}_{\alpha\beta}(x^{i},\phi_{p},\Omega,0,0)$ is identical to the action
of $\nabla_{\alpha}$ on the linear metric perturbation from a point mass on a
circular geodesic. The directional derivative $\vec{\partial}_{\cal V}$ then
accounts for the system’s slow movement through the parameter space.

###  IV.2 Evolution equations

The rates of change of the parameters ${\cal J}_{I}$ are likewise expanded in
powers of $\varepsilon$ at fixed $(\phi_{p},{\cal J}_{I})$:

$\displaystyle\frac{d\Omega}{dt}$ $\displaystyle=\varepsilon
F_{\Omega}^{(0)}(\Omega)+\varepsilon^{2}F_{\Omega}^{(1)}({\cal J}_{I})+{\cal
O}(\varepsilon^{3}),$ (42) $\displaystyle\frac{d\chi}{dt}$
$\displaystyle=\varepsilon F_{\chi}^{(1)}({\cal J}_{I})+{\cal
O}(\varepsilon^{2}),$ (43) $\displaystyle\frac{d\delta M}{dt}$
$\displaystyle=\varepsilon F_{\delta M}^{(1)}(\Omega)+{\cal
O}(\varepsilon^{2}),$ (44) $\displaystyle\frac{d\delta J}{dt}$
$\displaystyle=\varepsilon F_{\delta J}^{(1)}(\Omega)+{\cal
O}(\varepsilon^{2}).$ (45)

From these expansions we obtain the expansion for the coordinate velocity,

$\dot{z}^{\alpha}\equiv\frac{dz^{\alpha}}{dt}=\dot{z}^{\alpha}_{0}(\Omega)+\varepsilon\dot{z}^{\alpha}_{1}(\Omega)+O(\varepsilon^{2}),$
(46)

where

$\dot{z}_{0}^{\alpha}\equiv\frac{dz_{0}^{\alpha}}{dt}=(1,0,0,\Omega),$ (47)

and

$\dot{z}_{1}^{\alpha}=(0,\dot{r}_{0},0,0),$ (48)

with $\dot{r}_{0}(\Omega)=\frac{dr_{0}}{d\Omega}F_{\Omega}^{(0)}$; the proper
four-velocity is $u^{\alpha}=u^{\alpha}_{0}(\Omega)+\varepsilon
u^{t}_{0}(\Omega)\dot{z}^{\alpha}_{1}(\Omega)+\mathcal{O}(\varepsilon^{2})$,
with $u^{\alpha}_{0}=u^{t}_{0}\dot{z}^{\alpha}_{0}$ as in (23). The driving
forces $F_{Y}^{(n)}$, which govern the evolution, are to be determined from
the equations of motion (6) and from the Einstein field equations.

In the above expansions, the numerical labels within parentheses denote the
post-adiabatic order at which the quantity enters, while the numeric labels
without parentheses correspond to the explicit powers of $\varepsilon$.666This
statement assumes that we calculate $F^{(0)}_{\Omega}$ using the local self-
force. If we instead calculate it from energy fluxes to infinity and into the
horizon, then the leading horizon fluxes $F_{\delta M}^{(1)}(\Omega)$ and
$F_{\delta J}^{(1)}(\Omega)$ enter at 0PA order, and the first subleading
horizon fluxes enter at 1PA order. However, $\delta M$ and $\delta J$
themselves only enter at 1PA order in either approach. We have foreshadowed
the structure of the solution by indicating that 0PA quantities only depend on
$(\phi_{p},\Omega)$, and that $d{\cal J}_{I}/dt$ is independent of $\phi_{p}$
(an essential requirement for the separation between slow and fast evolution).

To expand the equations of motion (6), we also expand the forces and torque as

$\displaystyle F^{\alpha}$ $\displaystyle=\varepsilon F^{\alpha}_{1}({\cal
J}_{I})+\varepsilon^{2}F^{\alpha}_{2}({\cal J}_{I})+{\cal
O}(\varepsilon^{3}),$ (49) $\displaystyle N^{\alpha}$
$\displaystyle=\varepsilon N^{\alpha}_{1}({\cal J}_{I})+{\cal
O}(\varepsilon^{2}).$ (50)

Here the numeric labels correspond to the explicit powers of $\varepsilon$,
following the usual nomenclature for “first-order” and “second-order” self-
forces. The spin-dependent contributions are given by the expansions of (8).
More precisely, the first-order spin force
$F^{\alpha}_{(1,\chi)}(\Omega,\chi)$ is given by the test-body force (21);
$F^{\alpha}_{(2,\chi)}({\cal J}_{I})$ is given by the remainder of Eq. (8a);
and $N^{\alpha}_{1}$ is given by Eq. (8b) evaluated at $r_{p}=r_{\Omega}$ and
$u^{\alpha}=u^{\alpha}_{0}$. We will not require a more explicit expression
for $F^{\alpha}_{(2,\chi)}$. The torque reduces to

$N^{\mu}_{1}=-\frac{\chi}{2}\hat{z}^{\mu}(g^{\alpha\beta}-\hat{z}^{\alpha}\hat{z}^{\beta})u_{0}^{\gamma}\nabla^{0}_{\gamma}h^{\mathrm{R}(1)}_{\alpha\beta}=0.$
(51)

To see why this vanishes, note that it can be written as
$-\frac{\chi}{2}\hat{z}^{\mu}u_{0}^{\gamma}\nabla^{0}_{\gamma}\left[(g^{\alpha\beta}-\hat{z}^{\alpha}\hat{z}^{\beta})h^{\mathrm{R}(1)}_{\alpha\beta}\right]$.
For our quasicircular, spin-aligned system, the quantity in square brackets is
constant along the zeroth-order worldline $z^{\mu}_{0}$ (at fixed $\Omega$).
The derivative therefore vanishes.

Substituting all of the above expansions into Eqs. (6a) and (6b), we can
straightforwardly solve order by order in $\varepsilon$, equating coefficients
of powers of $\varepsilon$ at fixed ${\cal J}_{I}$ rather than at fixed $t$.
We obtain

$\displaystyle r_{0}(\Omega)=r_{\Omega},\quad r_{1}({\cal
J}_{I})=-\frac{F_{1}^{r}({\cal J}_{I})}{3(u_{0}^{t})^{2}f_{\Omega}\Omega^{2}}$
(52)

from the conservative sector [the radial component of Eq. (6a)], in perfect
analogy with Eq. (18). From the dissipative sector [Eq. (6b) and the $t$ or
$\phi$ component of (6a)], we obtain

$\displaystyle F_{\Omega}^{(0)}$ $\displaystyle=-\frac{3f_{\Omega}\Omega
F^{t}_{1}(\Omega)}{y(u_{0}^{t})^{4}(1-6y)},$ (53) $\displaystyle
F_{\Omega}^{(1)}$ $\displaystyle=-\frac{3f_{\Omega}\Omega F^{t}_{2}({\cal
J}_{I})}{y(u_{0}^{t})^{4}(1-6y)}-\frac{2\vec{\partial}_{\cal V}F_{1}^{r}({\cal
J}_{I})}{\sqrt{y}(u_{0}^{t})^{4}f_{\Omega}(1-6y)}$
$\displaystyle\quad-\frac{4(1-6y+12y^{2})F_{1}^{r}({\cal
J}_{I})F_{1}^{t}(\Omega)}{y^{3/2}(u_{0}^{t})^{6}f_{\Omega}(1-6y)^{2}},$ (54)
$\displaystyle F^{(1)}_{\chi}$ $\displaystyle=N^{z}_{1}=0.$ (55)

Here

$y\equiv M/r_{\Omega}=(M\Omega)^{2/3}.$ (56)

Equation (55) shows that the spin magnitude is constant at 1PA order.

The equations of motion do not determine the driving forces $F_{\delta
M}^{(1)}(\Omega)$ and $F_{\delta J}^{(1)}(\Omega)$. However, as shown in Ref.
[8], the second-order Einstein equations [(64c) below] dictate that these are
the usual fluxes of energy and angular momentum through the horizon due to a
point mass on a circular geodesic orbit of frequency $\Omega$ (reviewed in
Sec. V.4 below).

It is easy to see that if we rewrite Eqs. (33) and (42)–(45) in terms of a
“slow time” variable $\tilde{t}\equiv\varepsilon t$, then the equations have
asymptotic solutions

$\displaystyle\phi_{p}$
$\displaystyle=\varepsilon^{-1}\phi^{(0)}_{p}(\tilde{t})+\phi^{(1)}_{p}(\tilde{t})+{\cal
O}(\varepsilon),$ (57) $\displaystyle\Omega$
$\displaystyle=\Omega^{(0)}(\tilde{t})+\varepsilon\,\Omega^{(1)}(\tilde{t})+{\cal
O}(\varepsilon^{2}),$ (58) $\displaystyle\chi$ $\displaystyle=\chi^{(1)}+{\cal
O}(\varepsilon),$ (59) $\displaystyle\delta M$
$\displaystyle=M^{(1)}(\tilde{t})+{\cal O}(\varepsilon),$ (60)
$\displaystyle\delta J$ $\displaystyle=J^{(1)}(\tilde{t})+{\cal
O}(\varepsilon),$ (61)

with constant $\chi^{(1)}$, with
$d\phi^{(n)}_{p}/d\tilde{t}=\Omega^{(n)}(\tilde{t})$, and with easily worked
out equations for $d\Omega^{(n)}/d\tilde{t}$. These expansions in powers of
$\varepsilon$ at fixed slow time make clear the structure of the solution on
the radiation-reaction timescale $t\sim M/\varepsilon$ ($\tilde{t}\sim M$).
However, they are primarily useful at the final waveform-generation stage,
where they allow one to solve for the coefficients $\phi^{(n)}_{p}$,
$\Omega^{(n)}$, $M^{(1)}$, and $J^{(1)}$ without specifying a value of the
mass ratio. Prior to that stage, we treat $(\phi_{p},\Omega,\chi,\delta
M,\delta J)$ as independent variables on the binary’s phase space.

###  IV.3 Stress-energy tensor

Substituting the expansions (30)–(48) into the stress-energy (9), we obtain

$\displaystyle\\!\\!T^{\mu\nu}=\varepsilon
T_{(1)}^{\mu\nu}(x^{i},\phi_{p},\Omega)+\varepsilon^{2}T_{(2)}^{\mu\nu}(x^{i},\phi_{p},{\cal
J}_{I})+O(\varepsilon^{3}).\\!$ (62)

The leading term is still given by the leading term in Eq. (28) except that in
$\delta_{\phi}$ we do not replace $\phi_{p}$ with $\Omega t$. The subleading
term is now

$\\!\\!\\!\\!\\!T_{(2)}^{\alpha\beta}=\frac{u_{0}^{t}}{r_{0}^{3}}\left[2\left(r_{0}\dot{z}^{(\alpha}_{0}\dot{z}_{1}^{\beta)}-\dot{z}^{\alpha}_{0}\dot{z}^{\beta}_{0}r_{1}\right)\delta_{r}-\dot{z}_{0}^{\alpha}\dot{z}_{0}^{\beta}r_{0}r_{1}\delta^{\prime}_{r}\right]\delta_{\theta}\delta_{\phi}\\\
+T_{(\chi)}^{\alpha\beta},$ (63)

which is identical to Eq. (29b) except for (i) the change $r_{\chi}\to r_{1}$,
(ii) the addition of the $\dot{z}_{1}^{\beta}$ term, and (iii) the fact that
we again do not replace $\phi_{p}$ with $\Omega t$. In this expression,
$r_{1}$ receives a contribution from the first-order self-force as well as
from the first-order spin force (21). Since the spin does not contribute to
$\dot{z}_{1}^{\beta}$, its total contribution to $T^{\alpha\beta}_{(2)}$
therefore remains precisely (29b).

###  IV.4 Field equations and Fourier expansions

Substituting the expansions (35) and (62) into the Einstein equations
$G_{\mu\nu}[{\sf g}]=T_{\mu\nu}$ and equating coefficients of powers of
$\varepsilon$ at fixed $(\phi_{p},{\cal J}_{I})$, we obtain the hierarchy

$\displaystyle G_{\mu\nu}[g]$ $\displaystyle=0,$ (64a) $\displaystyle
G^{(1,0)}_{\mu\nu}[h^{(1)}]$ $\displaystyle=8\pi T^{(1)}_{\mu\nu},$ (64b)
$\displaystyle G^{(1,0)}_{\mu\nu}[h^{(2)}]$ $\displaystyle=8\pi
T^{(2)}_{\mu\nu}-G^{(2,0)}_{\mu\nu}[h^{(1)},h^{(1)}]$
$\displaystyle\quad-G^{(1,1)}_{\mu\nu}[h^{(1)}].$ (64c)

The operators $G^{(n,j)}_{\mu\nu}$ act on functions of $(x^{i},\phi_{p},{\cal
J}_{I})$. If we first expand $G_{\mu\nu}[g+h]$ in powers of $h_{\mu\nu}$, as
$G_{\mu\nu}+G^{(1)}_{\mu\nu}[h]+G^{(2)}_{\mu\nu}[h,h]+\ldots$, then
$G^{(n,j)}_{\mu\nu}$ is derived from $G^{(n)}_{\mu\nu}$ using the expansion of
$\nabla_{\alpha}$ in Eq. (38). Using that expansion, we see that
$G^{(1,0)}_{\mu\nu}$ is the standard linearized Einstein tensor with
$\nabla_{\alpha}\to\nabla^{(0)}_{\alpha}$. $G^{(n,1)}_{\mu\nu}$ is given by
the terms in $G^{(n)}_{\mu\nu}$ that are linear in the velocity ${\cal
V}_{I}$; since these do not couple to spin terms at 1PA order, we will not
display them here, but they can be extracted from Ref. [8].

The linear-in-spin equations are

$G^{(1,0)}_{\mu\nu}[h^{(2,\chi)}]=8\pi T^{(2,\chi)}_{\mu\nu},$ (65)

where $T^{(2,\chi)}_{\mu\nu}$ is given by Eq. (29b) as described around Eq.
(63). Note that $h^{(2,\chi)}_{\mu\nu}$, and the associated regular field
$h^{\mathrm{R}(2,\chi)}_{\mu\nu}$, includes the effect of $r_{\chi}$ in
addition to the effect of $T^{\mu\nu}_{(\chi)}$. This means that
$h^{\mathrm{R}(2,\chi)}_{\mu\nu}$ differs from the field
$h^{\mathrm{R}(\chi)}_{\mu\nu}$ in Sec. II.2 by terms proportional to
$r_{\chi}$ (plus order-$\varepsilon^{3}$ differences).

Since all functions of $\phi_{p}$ are periodic, we can expand them in Fourier
series. For example

$h^{(n)}_{\alpha\beta}(x^{i},\phi_{p},{\cal
J}_{I})=\sum_{m=-\infty}^{\infty}h^{(n,m)}_{\alpha\beta}(x^{i},{\cal
J}_{I})e^{-im\phi_{p}}.$ (66)

We then have $\frac{\partial}{\partial\phi_{p}}\to-im$ when acting on
individual modes, implying

$\nabla^{(0)}_{\alpha}\to e^{i}_{\alpha}\frac{\partial}{\partial
x^{i}}-it_{\alpha}\omega_{m}+\text{Christoffel terms},$ (67)

where $\omega_{m}\equiv m\Omega$. The label $m$ here also serves as the
azimuthal mode number, such that

$h^{(n)}_{\alpha\beta}(x^{i},\phi_{p},{\cal
J}_{I})=\sum_{m=-\infty}^{\infty}h^{(n,m)}_{\alpha\beta}(r,\theta,{\cal
J}_{I})e^{im(\phi-\phi_{p})}.$ (68)

(We abuse notation by using $h^{(n,m)}_{\alpha\beta}$ for the coefficients in
both decompositions.)

The action of $\nabla^{(0)}_{\alpha}$ in Eq. (67) is identical to the action
of an ordinary covariant derivative acting on a Fourier series
$\sum_{m=-\infty}^{\infty}h^{(n,m)}_{\alpha\beta}(x^{i})e^{-im\Omega t}$, even
though $\phi_{p}$ in Eq. (66) is not equal to $\Omega t$. Analogously, when
acting on a Fourier expansion of the form (66), the leading-order linearized
Einstein tensor $G^{(1,0)}_{\mu\nu}$ is identical to the ordinary linearized
Einstein tensor acting on an ordinary Fourier series with modes $e^{-im\Omega
t}$. Equations (64b) and (64c) therefore reduce to the familiar form of
linearized Einstein equations in the frequency domain.

###  IV.5 Summary: two-timescale evolution with spin

Our analysis has shown that linear-in-spin effects are easily incorporated
into the two-timescale evolution and waveform-generation scheme of Ref. [8].
That scheme can be summarized in two conceptually simple steps:777There are
also two simplifications that are important in practice but not conceptually
essential. First, the offline computations are done on the one-dimensional
space of $\Omega$ values rather than the three-dimensional space of
$(\Omega,\delta M,\delta J)$ values because we only require effects linear in
$\delta M$ and $\delta J$, allowing us to compute the coefficients of those
effects without specifying values of $(\delta M,\delta J)$. Second, by
assuming the ansatzes (57)–(61), we can convert Eqs. (33) and (42)–(45) into
equations for the coefficients $\phi^{(0)}_{p}(\tilde{t})$,
$\phi^{(1)}_{p}(\tilde{t})$, etc., which can be solved for and stored without
specifying $\mu$. The inspiral trajectory and waveform for a given $\mu$ can
then be generated effectively instantly using the stored solutions.

  1. 1.

Offline computations. On a grid of $\Omega$ values, solve the field equations
(64b) and (64c) for the mode amplitudes $h^{(n,m)}_{\alpha\beta}$. From them,
compute the forcing functions $F^{(0)}_{\Omega}(\Omega)$, $F^{(0)}_{\delta
M}(\Omega)$, $F^{(0)}_{\delta J}(\Omega)$, and
$F^{(1)}_{\Omega}(\Omega,\chi=0,\delta M,\delta J)$. These calculations only
involve the coefficients of the phase factors $e^{-im\phi_{p}}$, never the
orbital phase $\phi_{p}(t,\varepsilon)$ itself; the orbital phase factors out
of the computations.

  2. 2.

Online simulation. Using the stored forcing functions, choose a value of $\mu$
and solve Eqs. (33) and (42)–(45) for the phase-space trajectory

$(\phi_{p}(t,\varepsilon),\Omega(t,\varepsilon),\chi=0,\delta
M(t,\varepsilon),\delta J(t,\varepsilon)).$

From that trajectory and the mode amplitudes, $h^{(n,m)}_{\alpha\beta}$,
generate the waveform
$\lim_{r\to\infty}r\sum_{n,m}\varepsilon^{n}h^{(n,m)}_{\alpha\beta}(x^{i},{\cal
J}_{I}(t,\varepsilon))e^{-im\phi_{p}(t,\varepsilon)}$.

To incorporate the spin into this framework, we simply add the following to
the first step:

  * 1.*

On the grid of $\Omega$ values, solve the field equation (65) for the linear-
in-spin contribution to the mode amplitudes $h^{(2,m)}_{\alpha\beta}$. From
them and $h^{(1,m)}_{\alpha\beta}$, compute the linear-in-spin contribution to
$F^{(1)}_{\Omega}$.

This can be done without specifying a value of $\chi$. In Step 2, we then
simply set a freely specifiable nonzero value $\chi$, which remains constant
by virtue of Eq. (55), add the linear-in-spin term to $F^{(1)}_{\Omega}$ and
$h^{(2,m)}_{\alpha\beta}$, and proceed as above to generate the waveform.

In principle, the spin contribution to $F^{(1)}_{\Omega}$ can be computed
directly from the local expression (54). However, we can also extract it from
the asymptotic fluxes of gravitational waves. As shown in [39], the quantity
$\Xi$ defined in Eq. (12) (the energy or angular momentum) satisfies a flux-
balance law: neglecting higher-order spin effects, and neglecting all other
${\cal O}(\varepsilon^{2})$ terms in the equation of motion, the rate of
change of $\Xi$ due to the local force and torque, which can be written as

$\left\langle\frac{D\Xi}{d\tau}\right\rangle=\frac{1}{2}\left\langle
u^{\alpha}u^{\beta}\mathcal{L}_{\xi}h_{\alpha\beta}^{\rm
R}+\mu\tilde{S}^{\alpha\beta}u^{\gamma}\nabla_{\alpha}\mathcal{L}_{\xi}h_{\beta\gamma}^{\rm
R}\right\rangle,$ (69)

is equal (with opposite sign and up to a factor of $u^{t}$) to the total flux
of gravitational-wave energy (or angular momentum) out to infinity and into
the black hole.888The balance law holds for a generic bound orbit in Kerr
spacetime, in which case this statement applies on average, with the average
over radial and polar oscillations denoted by angle brackets. For
quasicircular, equatorial orbits, no average is required. For example, with
the caveats of what is being neglected, the rate of change of the energy $E$
is

$\frac{dE}{dt}=-(\dot{E}^{+}+\dot{E}^{-})\equiv-{\cal F},$ (70)

where $\dot{E}^{+}$ and $\dot{E}^{-}$ are the gravitational-wave energy fluxes
(per unit $\mu$) to infinity and into the horizon, respectively.

The energy has the form $E({\cal J}_{I})=E_{0}(\Omega)+\varepsilon
E_{\chi}(\Omega,\chi)+\varepsilon E_{SF}(\Omega,\delta M,\delta J)+{\cal
O}(\varepsilon^{2})$, where the first two terms are given by the test-body
energy (24), and the final term comes from the first-order radial self-force
term in $-u^{\alpha}\xi_{\alpha}=f(r_{p})u^{t}$. Substituting this expression
for $E$ into the balance law (70), applying the chain rule, substituting
$d\Omega/dt=\varepsilon F^{(0)}_{\Omega}+{\cal O}(\varepsilon^{2})$ and
$E_{0}=f_{\Omega}u^{t}_{0}(\Omega)$, and rearranging for $F^{(0)}_{\Omega}$,
we immediately find

$F^{(0)}_{\Omega}(\Omega)=\frac{3\Omega}{y(u^{t}_{0})^{3}(1-6y)}{\cal
F}_{(1)}(\Omega).$ (71)

This is the standard adiabatic, flux-driven evolution of the frequency, and
${\cal F}_{(1)}(\Omega)$ is identical to the standard leading-order flux due
to a point mass on a geodesic circular orbit of frequency $\Omega$.

We can similarly pick off the linear-in-spin, ${\cal O}(\varepsilon^{2})$
terms in the balance law (70). We write their contribution as

$\left(\frac{dE}{dt}\right)_{\chi}=-{\cal F}_{\chi},$ (72)

noting that these denote coefficients of $\chi$ rather than including the
linear factor of $\chi$. Since we are not neglecting any linear-in-spin,
${\cal O}(\varepsilon^{2})$ terms, this formula is exact, unlike (70). Again
applying the chain rule, using $d\chi/dt={\cal O}(\varepsilon^{2})$, and
picking out the linear-in-spin terms, we find the left-hand side of Eq. (72)
evaluates to

$\left(\frac{dE_{\chi}}{dt}\right)_{\chi}=\frac{\partial
E_{0}}{\partial\Omega}F^{(1)}_{\Omega}+\left(\frac{\partial
E_{\chi}}{\partial\Omega}\right)_{\chi}F^{(0)}_{\Omega}.$ (73)

Substituting this into Eq. (72) and rearranging, we obtain our desired result:

$F^{(1,\chi)}_{\Omega}(\Omega)=\frac{3\Omega}{y(u^{t}_{0})^{3}(1-6y)}{\cal
F}_{\chi}(\Omega)\\\
-\frac{3\mu\Omega^{2}(5-12y)}{(u^{t}_{0})^{3}y(1-6y)^{2}}{\cal
F}_{(1)}(\Omega).$ (74)

This is the most essential input for a 1PA evolution including spin. Its main
ingredient, ${\cal F}_{\chi}(\Omega)$, has the distinct advantage of being
more easily computatble than the local metric perturbation, as it is derived
from the amplitudes $h^{(n,m)}_{\alpha\beta}$ at infinity and the black hole’s
horizon. It has been calculated previously by several groups using a fixed-
radius formulation [39, 40, 41]; we review its computation using a more
efficient fixed-frequency formulation in the next section.

Before moving to the next section, we note that although we have not included
quadratic-in-spin terms in our calculations, they would not enter into the 1PA
waveform. This follows from the fact that at ${\cal O}(\varepsilon^{2})$, they
only enter in the form of test-body terms in the background spacetime. Such
terms are known to be purely conservative [51, 52], and only dissipative
${\cal O}(\varepsilon^{2})$ terms can enter at 1PA (specifically, through
$F_{\Omega}^{(1)}$).

##  V Regge-Wheeler-Zerilli formalism

The Regge-Wheeler-Zerilli (RWZ) formalism [53, 54] specialises the linearized
Einstein field equations to perturbations of a Schwarzschild background
spacetime. The RWZ equations follow from a tensor spherical harmonic
decomposition of the field equations which separates the odd and even parity
perturbations. Physically speaking, the even parity sector describes fields
that are invariant under the transformation
$(\theta,\phi)\rightarrow(\pi-\theta,\pi+\phi)$ while the odd parity sector
describes fields that are changed by a factor of $-1$ under the same
transformation.

In the RWZ formalism, two master functions satisfying the RWZ equations are
defined in terms of the metric perturbation amplitudes, one for the odd parity
sector and one for the even parity sector. The master functions can then be
used to reconstruct the metric perturbation amplitudes and thus the metric
perturbation. We make a gauge choice in this final step: the odd parity
perturbations are fixed into the Regge-Wheeler gauge while the even parity
perturbations are fixed into the Zerilli gauge.

In this work we follow the gauge-invariant approach to the RWZ formalism
detailed in [55, 56]. In the even parity sector we use the Zerilli-Moncrief
[57] master function, while in the odd parity sector we use the Cunningham-
Price-Moncrief [58] master function.

###  V.1 The Regge-Wheeler-Zerilli equation

Both the Zerilli-Moncrief and the Cunningham-Price-Moncrief master functions
satisfy a wave equation (the RWZ equation), which in the frequency domain is
of the form

$\left[\dfrac{\partial^{2}}{\partial
r_{*}^{2}}-V_{\ell}(r)+\omega^{2}\right]\psi_{\ell m\omega}(r)=Z_{\ell
m\omega}(r),$ (75)

where $r_{*}=r+2M\ln(r/2M-1)$ is the Schwarzschild tortoise coordinate such
that $\dfrac{dr}{dr_{*}}=f(r)$, and where $\omega=m\Omega$ for circular
orbits. The form of the potential $V_{\ell}(r)$ and the source term $Z_{\ell
m\omega}$ depends on the master function, $\psi_{\ell m\omega}$, and therefore
is different in each parity sector, with the Zerilli potential in the even
sector and the Regge-Wheeler potential in the odd sector. In each parity
sector, the source term derives from the stress-energy tensor due to a small
spinning point particle, Eq. (28), and is of the form

$Z_{\ell m\omega}=\left(\bar{G}_{\ell
m\omega}\delta_{r_{\Omega}}+\bar{F}_{\ell
m\omega}\delta_{r_{\Omega}}^{\prime}+\bar{H}_{\ell
m\omega}\delta_{r_{\Omega}}^{\prime\prime}\right),$ (76)

where as previously we use the shorthand
$\delta_{r_{\Omega}}\equiv\delta(r-{r_{\Omega}})$. The functions
$\bar{G}_{\ell m\omega}$, $\bar{F}_{\ell m\omega}$ and $\bar{H}_{\ell
m\omega}$ depend only on ${r_{\Omega}}$, $\chi$, $\mu$ and $M$ 999In the non-
spinning limit, $\bar{H}_{\ell m\omega}\rightarrow 0$. When writing the source
at fixed frequency, in the even parity sector it also turns out that
$\bar{H}^{\rm even}_{\ell m\omega}=0.$; explicit expressions for these and for
the potentials $V_{\ell}(r)$ are given in Appendix B.

###  V.2 Retarded solutions to the RWZ equations

The RWZ equation, Eq. (75), admits two linearly independent homogeneous
solutions. There is flexibility in the particular choice of basis of
homogeneous solutions. For radiative ($\omega\neq 0$) modes we choose to work
with “in” and “up” solutions, which satisfy boundary conditions representing
radiation that is purely ingoing into the future horizon and purely outgoing
to future null infinity, respectively. Using hats to signify that they are
unit-normalised homogeneous solutions, the “in” and “up” solutions therefore
have the asymptotic behaviour

$\displaystyle\hat{R}_{\ell m\omega}^{-}(r\rightarrow 2M)$
$\displaystyle=e^{-i\omega_{m}r_{*}},$ (77a) $\displaystyle\hat{R}_{\ell
m\omega}^{+}(r\rightarrow\infty)$ $\displaystyle=e^{i\omega_{m}r_{*}}.$ (77b)

For the non-radiative ($\omega=0$) modes we use analytic solutions in terms of
hypergeometric functions,

$\displaystyle\hat{R}_{\ell m0}^{-}(r)$
$\displaystyle=\,_{2}F_{1}\left(l-s+1,l+s+1;1;f(r)\right)$ (78a)
$\displaystyle\hat{R}_{\ell m0}^{+}(r)$
$\displaystyle=\,_{2}F_{1}\left(l-s+1,l+s+1;2(l+1);2M/r\right),$ (78b)

where $s=2$ for the gravitational RWZ equations. Here, we adopt the same
notation as for the radiative modes but our solutions are now chosen on the
basis that they are regular at the horizon and at infinity, respectively.

In this work we make use of the ReggeWheeler package from the Black Hole
Perturbation Toolkit [59] to obtain numerical solutions to the homogeneous RWZ
equations. We then use the method of variation of parameters to find the
inhomogeneous solution of (75) in terms of these homogeneous solutions,

$\psi_{\ell m\omega}(r)=c^{+}_{\ell m\omega}(r)\hat{R}_{\ell
m\omega}^{+}(r)+c^{-}_{\ell m\omega}(r)\hat{R}_{\ell m\omega}^{-}(r),$ (79)

where

$\displaystyle c^{+}_{\ell m\omega}(r)\equiv\frac{1}{W_{\ell
m\omega}}\int^{r}\frac{\hat{R}_{\ell m\omega}^{-}(r^{\prime})Z_{\ell
m\omega}(r^{\prime})}{f(r^{\prime})}dr^{\prime},$ (80a) $\displaystyle
c^{-}_{\ell m\omega}(r)\equiv\frac{1}{W_{\ell
m\omega}}\int_{r}\frac{\hat{R}_{\ell m\omega}^{+}(r^{\prime})Z_{\ell
m\omega}(r^{\prime})}{f(r^{\prime})}dr^{\prime},$ (80b)

and the factor of $f(r)$ in the denominator of the integrand comes from
changing integration variable from $r_{*}$ to $r$. In the method of variation
of parameters, the Wronskian

$W_{\ell m\omega}\equiv\hat{R}_{\ell m\omega}^{-}\frac{d\hat{R}_{\ell
m\omega}^{+}}{dr_{*}}-\hat{R}_{\ell m\omega}^{+}\frac{d\hat{R}_{\ell
m\omega}^{-}}{dr_{*}},$ (81)

typically appears inside the integral. However, since there is no first-
derivative term in Eq. (75), by Abel’s identity the Wronskian is a constant so
it may be taken outside of the integral and evaluated at any convenient
radius.

Substituting the source, $Z_{\ell m\omega}(r)$, given by Eq. (76) into Eqs.
(80a) and (80b), integrating by parts and paying careful attention to the
boundary terms, the inhomogeneous solution in Eq. (79) becomes

$\psi_{\ell m\omega}(r)=R_{\ell m\omega}^{+}(r)\Theta^{+}_{r_{\Omega}}+R_{\ell
m\omega}^{-}(r)\Theta^{-}_{r_{\Omega}}+X_{\ell m\omega}\delta_{r_{\Omega}},$
(82)

where we have denoted the Heaviside step functions by
$\Theta^{+}_{r_{\Omega}}\equiv\Theta[r-r_{\Omega}]$ and
$\Theta^{-}_{r_{\Omega}}\equiv\Theta[r_{\Omega}-r]$, and $X_{\ell
m\omega}\equiv\frac{\bar{H}_{\ell m\omega}}{f^{2}_{\Omega}}$ is a constant. In
the above expressions we have introduced the shorthand $R^{\pm}_{\ell
m\omega}(r)\equiv C^{\pm}_{\ell m\omega}\hat{R}^{\pm}_{\ell m\omega}(r)$. The
constant matching coefficients are given by

$\displaystyle C_{\ell m\omega}^{\pm}$ $\displaystyle=\frac{1}{W_{\ell
m\omega}}\frac{\hat{R}_{\ell
m\omega}^{\mp}(r_{\Omega})}{f_{\Omega}}\bar{G}_{\ell m\omega}$
$\displaystyle\qquad-\frac{1}{W_{\ell
m\omega}}\frac{d}{dr}\left.\left(\frac{\hat{R}_{\ell
m\omega}^{\mp}(r)}{f(r)}\right)\right\rvert_{r_{\Omega}}\bar{F}_{\ell
m\omega}$ $\displaystyle\qquad+\frac{1}{W_{\ell
m\omega}}\frac{d^{2}}{dr^{2}}\left.\left(\frac{\hat{R}_{\ell
m\omega}^{\mp}(r)}{f(r)}\right)\right\rvert_{r_{\Omega}}\bar{H}_{\ell
m\omega}.$ (83)

In addition to the delta singularity, the resulting master function and its
derivative have jump discontinuities at $r_{\Omega}$. The expressions for the
jumps are obtained by substituting (82) and (76) into (75). Defining
$[[\psi(r_{\Omega})]]_{\ell m\omega}\equiv R_{\ell
m\omega}^{+}(r_{\Omega})-R_{\ell m\omega}^{-}(r_{\Omega})$ and
$[[\psi^{\prime}(r_{\Omega})]]_{\ell m\omega}\equiv\partial_{r}R_{\ell
m\omega}^{+}(r_{\Omega})-\partial_{r}R_{\ell m\omega}^{-}(r_{\Omega})$;

$\displaystyle f_{\Omega}^{2}[[\psi(r_{\Omega})]]_{\ell m\omega}$
$\displaystyle=\bar{F}_{\ell
m\omega}+3\frac{f^{\prime}_{\Omega}}{f_{\Omega}}\bar{H}_{\ell m\omega},$ (84)
$\displaystyle f_{\Omega}^{2}[[\psi^{\prime}(r_{\Omega})]]_{\ell m\omega}$
$\displaystyle=\left(V_{l}(r_{\Omega})-\omega^{2}+2(f^{\prime}_{\Omega})^{2}-f_{\Omega}f^{\prime\prime}_{\Omega}\right)\frac{\bar{H}_{\ell
m\omega}}{f_{\Omega}^{2}}$ $\displaystyle\qquad+\bar{G}_{\ell
m\omega}+\frac{f^{\prime}_{\Omega}}{f_{\Omega}}\bar{F}_{\ell m\omega},$ (85)

where quantities with the subscript $\Omega$ are evaluated at $r_{\Omega}$,
e.g. $f_{\Omega}\equiv f(r_{\Omega})$.

Finally, obtaining the time domain master functions from their frequency
domain counterparts is trivial for the case of circular equatorial orbits,
since there is only a single frequency per $m$ mode and thus the time domain
master functions are given by

$\Psi_{\ell m}(t,r)=\psi_{\ell m\omega}(r)e^{-i\omega_{m}t}.$ (86)

In our two-timescale expansion, this instead becomes

$\Psi_{\ell m}(t,r)=\psi_{\ell m\omega}(r)e^{-im\phi_{p}(t)}.$ (87)

###  V.3 Metric reconstruction

We next obtain the actual metric perturbation sourced by the stress-energy
(28) due to a spinning point particle. This metric perturbation has two
components: a _radiative_ piece which can be reconstructed from the RWZ master
functions; and a _completion_ piece that fully captures the mass and angular
momentum content of the perturbation.

####  V.3.1 Reconstruction from RWZ master functions

The radiative part of the metric perturbation can be reconstructed by applying
differential operators to the RWZ master functions along with source terms
involving the stress-energy tensor. The metric perturbation one obtains is in
the Regge-Wheeler/ Zerilli gauges. In terms of the complex vector
$m^{\mu}=\frac{1}{\sqrt{2}r}\\{0,0,1,i\csc\theta\\}$ (with complex conjugate
$\bar{m}^{\mu}$) the Regge-Wheeler gauge condition is equivalent to the
conditions $h_{mm}=0=h_{tm}=h_{rm}$ in the even-parity sector and to the
condition $h_{mm}=0$ in the odd-parity sector. This reconstruction procedure
is by now well-established and we omit the details here as they are given in
full in Refs. [56] and [5].

####  V.3.2 Metric completion at fixed frequency

The completion part of the metric perturbation can be obtained by solving the
harmonic decomposed linearized Einstein equations directly for the $\ell=0$
and $\ell=1$ modes. The previous RWZ gauge conditions do not fully fix the
gauge in this case since by definition $h_{mm}=0$ for $\ell=0,1$ and
$h_{tm}=0=h_{rm}$ for $\ell=0$. Instead we fix the residual gauge freedom by
working in a “RWZ-like” gauge (see Appendix G for further details).

The metric completion pieces for a secondary with aligned spin in a circular
orbit in a Schwarzschild background spacetime were first given in a RWZ-like
gauge in [22]. These results were derived at fixed radius and their
contribution to the redshift was later rewritten at fixed frequency. For
completeness in Appendix G we derive the completion pieces directly within a
fixed frequency calculation.

The key result from this derivation is that for $\ell=0$ the only non-zero
components of the metric perturbation for a spin-aligned secondary in a
circular orbit in a Zerilli-like gauge are:

$\displaystyle h^{\ell=0}_{tt}$ $\displaystyle=\frac{2\mu
E}{r}\bigg{\\{}\Theta[r-r_{\Omega}]$
$\displaystyle\quad+\frac{rf}{r_{\Omega}f_{\Omega}}\left[1-\chi\mu\Omega\left(\frac{r_{\Omega}-3M}{r_{\Omega}-2M}\right)\right]\Theta[r_{\Omega}-r]\bigg{\\}}$
$\displaystyle h^{\ell=0}_{rr}$ $\displaystyle=\frac{2\mu
E}{rf^{2}}\Theta[r-r_{\Omega}].$ (88)

Note that the expression in Ref. [22] for $h_{rr}$ featured a Dirac delta that
does not appear when expressed at fixed frequency.

The $\ell=1$ contribution is comprised of the odd parity $m=0$ mode and the
even parity $m=\pm 1$ modes. The even parity modes are pure gauge modes away
from the worldine and also do not contribute to the redshift gauge invariant
in Section VI. The odd parity mode does contribute to the redshift invariant.
In a Regge-Wheeler-like gauge, the only non-zero retarded metric perturbation
component for a spin-aligned secondary in a circular orbit is

$\displaystyle h_{t\phi}$ $\displaystyle=-\frac{2\mu
J}{r}\bigg{\\{}\Theta[r-r_{\Omega}]$
$\displaystyle\quad+\frac{r^{3}}{r_{\Omega}^{3}}\left[1+\frac{3}{2}\chi\mu\frac{\Omega}{M}\left(r_{\Omega}-3M\right)\right]\Theta[r_{\Omega}-r]\bigg{\\}}\sin^{2}(\theta).$
(89)

####  V.3.3 Evaluation on the worldline

There is a subtlety in the fixed-frequency formulation in that when we wish to
evaluate “on the worldline” (for example, to compute the local force or
Detweiler’s redshift) we must evaluate the metric perturbation at $r=r_{p}$.
When evaluating a given expression, we must therefore substitute the expansion
$r_{p}=r_{\Omega}+\varepsilon r_{\chi}$, re-expand and truncate in order to
obtain a consistent result involving the metric perturbation and its
derivatives evaluated at $r=r_{\Omega}$. An example of this final step is
given in Sec. VI, in which which we use it when computing Detweiler’s
redshift.

###  V.4 Gravitational wave fluxes

The gravitational wave energy and angular momentum fluxes at $r_{*}=\pm\infty$
(more formally, at future null infinity and at the event horizon) are in
general gauge invariant. They also involve only dissipative contributions to
the metric perturbation and thus have the advantages of not requiring
regularization and of involving a rapidly-convergent sum over spherical
harmonic modes.

Provided the homogeneous radial solutions $\hat{R}_{\ell m\omega}^{\mp}(r)$
have been normalised as unit in-going/outgoing waves at $r_{*}=\mp\infty$, the
specific (per unit $\mu$) energy and specific angular momentum fluxes (with
respect to coordinate time $t$) are [56]

$\dot{E}^{\pm}=\displaystyle\sum_{\ell=2}^{\infty}\displaystyle\sum^{\ell}_{m=-\ell}\dot{E}_{\ell
m}^{\pm},\quad\dot{J}^{\pm}=\displaystyle\sum_{\ell=2}^{\infty}\displaystyle\sum^{\ell}_{m=-\ell}\dot{J}_{\ell
m}^{\pm}.$

where

$\displaystyle\dot{E}_{\ell
m\omega}^{\pm}=\frac{1}{64\pi\mu}\frac{(\ell+2)!}{(\ell-2)!}\omega_{m}^{2}\left|C_{\ell
m\omega}^{\pm}\right|^{2},$ (90) $\displaystyle\dot{J}_{\ell
m\omega}^{\pm}=\frac{m}{64\pi\mu}\frac{(\ell+2)!}{(\ell-2)!}\omega_{m}\left|C_{\ell
m\omega}^{\pm}\right|^{2},$ (91)

are the harmonic modes of the energy and angular momentum fluxes at
$r_{*}=\pm\infty$ respectively. The $C_{\ell m\omega}^{\pm}$ are those given
by Eqn. (V.2) and also depend on $\mu$. Note that neither the zero frequency
$m=0$ modes nor the conservative $\ell=0$ and $\ell=1$ modes contribute to the
fluxes so the mode-sum starts at $\ell=2$.

##  VI Detweiler’s redshift invariant

Detweiler’s gauge invariant redshift —
$z\equiv\frac{d\hat{\tau}}{dt}=1/\hat{u}^{t}$— was generalised to the case of
an (anti-)aligned spinning secondary in a circular orbit on a Schwarzschild
background spacetime in Ref. [22]. Written in a “mixed” form in terms of the
invariant frequency variable $y=M/r_{\Omega}$ and orbital radius
$r_{p}=r_{\Omega}+\varepsilon r_{\chi}(y)$, the explicit form of the redshift
is given by101010This expression is valid for (anti-)aligned spinning
secondaries in a circular orbit parameterised at fixed frequency only.

$z(y,r_{p},\varepsilon)=z_{(0)}(y)+\varepsilon
z_{(1)}(y,r_{p})+\varepsilon^{2}z_{(2)}(y,r_{p})+\mathcal{O}(\varepsilon^{3}),$
(92)

where we have defined

$\displaystyle z_{(0)}(y)$ $\displaystyle\equiv\sqrt{1-3y},$ (93)
$\displaystyle z_{(1)}(y,r_{p})$
$\displaystyle\equiv-\frac{1}{2\sqrt{1-3y}}h^{\rm R\,(1)}_{kk}(r_{p}),$ (94)

and as usual we restrict our analysis of the $O(\varepsilon^{2})$ term to spin
effects only, which contribute

$z_{(2)}(y,r_{p})\equiv-\frac{1}{2\sqrt{1-3y}}\Big{[}h^{\rm
R\,(\chi)}_{kk}(r_{p})\\\ +\mu\chi y^{1/2}\partial_{r}h^{\rm
R\,(1)}_{kk}(r_{p})\Big{]}.$ (95)

We have also defined the helical Killing vector
$k^{\alpha}\equiv\xi_{(t)}^{\alpha}+\Omega\xi_{(\phi)}^{\alpha}$ and $h^{\rm
R}_{kk}\equiv h^{\rm R}_{\alpha\beta}k^{\alpha}k^{\beta}$. The radial
derivative term in Eq. (95) follows from the MPD equations in the perturbed
spacetime, and can be interpreted as ensuring the redshift is gauge invariant
through linear order in spin.

From Eq. (92) we can obtain a simpler expression in which $z$ is fully
expanded in powers of $\varepsilon$ at fixed $\Omega$:

$z(y,\varepsilon)=z_{(0)}(y)+\varepsilon
z_{(1)}(y)+\varepsilon^{2}z_{(2)}(y)+\mathcal{O}(\varepsilon^{3}).$ (96)

The single-argument functions here are related to the two-argument functions
in Eq. (92) by $z_{(1)}(y)=z_{(1)}(y,r_{\Omega})$ and
$z_{(2)}(y)=z_{(2)}(y,r_{\Omega})+r_{\chi}(y)\partial_{r_{\Omega}}z_{(1)}(y,r_{\Omega})$.
To understand the expansion in more detail, note that in the equations above,
$h^{\rm R\,(1)}_{kk}(r_{p})$ and $h^{\rm R\,(\chi)}_{kk}(r_{p})$ are evaluated
at the field point $r=r_{p}$, but they also denote fields that are generated
by a particle at $r_{p}$. We can make this explicit by writing the metric
perturbation’s dependence on the field point $x^{\alpha}$ and on the secondary
body’s worldline $z^{\alpha^{\prime}}$ as
$h_{kk}=h_{kk}(z^{\alpha^{\prime}},x^{\alpha})$. Substituting
$r_{p}=r_{\Omega}+\varepsilon r_{\chi}$ and expanding both arguments, we
obtain

$\displaystyle h^{\rm R\,(1)}_{kk}(z^{\alpha^{\prime}},z^{\alpha})$
$\displaystyle=h^{\rm
R\,(1)}_{kk}(z^{\alpha^{\prime}}_{0},z^{\alpha}_{0})+\varepsilon
r_{\chi}\partial_{r_{\Omega}^{\prime}}h_{kk}^{\mathrm{R}\,(1)}(z^{\alpha^{\prime}}_{0},z^{\alpha}_{0})$
$\displaystyle\quad+\varepsilon
r_{\chi}\partial_{r_{\Omega}}h_{kk}^{\mathrm{R}\,(1)}(z^{\alpha^{\prime}}_{0},z^{\alpha}_{0})+\mathcal{O}(\varepsilon^{2}),$
(97)

and $h^{\rm R\,(\chi)}_{kk}(z^{\alpha^{\prime}},z^{\alpha})=h^{\rm
R\,(\chi)}_{kk}(z_{0}^{\alpha^{\prime}},z^{\alpha}_{0})+\mathcal{O}(\varepsilon)$.
Inspecting Eqs. (92), (94), and (95) again, and noting $r_{\chi}=-\mu\chi
y^{1/2}$, we see that the final term in Eq. (VI) cancels the term involving
$\mu\chi y^{1/2}\partial_{r}h^{\rm R\,(1)}_{kk}$ in $z_{(2)}(y,r_{p})$. The
second term in Eq. (VI) combines with the $h^{\rm R\,(\chi)}_{kk}$ term in
$z_{(2)}(y,r_{p})$ to give the total linear-in-spin contribution to $h^{\rm
R(2)}_{kk}$:

$h_{kk}^{\mathrm{R}(2,\chi)}(z_{0}^{\alpha^{\prime}},z_{0}^{\alpha})=h_{kk}^{\mathrm{R}(\chi)}(z^{\alpha^{\prime}}_{0},z^{\alpha}_{0})\\\
+r_{\chi}\partial_{r_{\Omega}^{\prime}}h_{kk}^{\mathrm{R}(1)}(z^{\alpha^{\prime}}_{0},z^{\alpha}_{0}).$
(98)

The two terms in $h_{kk}^{\mathrm{R}(2,\chi)}$ correspond to the two terms in
$T^{\mu\nu}_{(2)}$ from Eq. (29b).

The explicit version of Eq. (96) is then simply

$z=\sqrt{1-3y}-\frac{1}{2\sqrt{1-3y}}\left[\varepsilon h^{\rm
R\,(1)}_{kk}(r_{\Omega})+\varepsilon^{2}h_{kk}^{\rm
R\,(2,\chi)}(r_{\Omega})\right],$ (99)

where we only keep the linear-in-spin second-order terms,

$z_{(2,\chi)}=-\frac{h_{kk}^{\rm R\,(2,\chi)}(r_{\Omega})}{2\sqrt{1-3y}}.$
(100)

This formula can also be deduced from Eq. (100) in Ref. [60], which is valid
for any radial perturbing force and includes the complete
order-$\varepsilon^{2}$ term.

##  VII Regularization

A consequence of modelling a compact object by a Dirac delta function and its
derivatives is that the retarded metric perturbation is singular on the
worldline. This leads to discontinuities in the RWZ master functions across
the worldline and the sum over spherical harmonic modes will not converge
there. In fact, equation (82) shows that for a spinning body modelled in this
way there is also a delta singularity in the master functions on the worldline
111111As $X_{\ell m\omega}\equiv\frac{\bar{H}_{\ell m\omega}}{f^{2}_{\Omega}}$
and $\bar{H}^{\rm even}_{\ell m\omega}=0$ when parameterised at fixed
frequency, the delta singularity on the worldline vanishes in the even parity
master function.. These spurious divergences are not fundamental, and can be
unambiguously avoided by a more careful treatment that uses matched asymptotic
expansions instead of distributions [5]. The net result of such an analysis is
that at leading order in perturbation theory (including subleading order spin
terms) we recover the point particle approximation with distributional
sources, along with a well-defined regularization procedure that involves
subtracting an appropriate singular field from the retarded field to produce a
so-called residual field. That is,

$h_{\mu\nu}^{\rm R}=h_{\mu\nu}^{\rm ret}-h_{\mu\nu}^{\rm S},$ (101)

where the superscripts ‘$\rm R$’ and ‘$\rm S$’ denote the regular and singular
pieces, respectively. It is this residual field that appears in local
quantities evaluated on the worldline.

The spin’s contribution to the local field $h_{\mu\nu}^{(2)\rm ret}$ near the
particle was derived in Ref. [61] as a local expansion in powers of distance
from the worldline, through order (distance)0, and in Ref. [46] through linear
order in distance. Ref. [46] also defined a split into singular and regular
pieces. However, the field was expressed in a local coordinate system; some
additional work is required to put $h_{\mu\nu}^{\rm S(2,\chi)}$ in a practical
form that can be used to calculate $h^{\rm R}_{kk}$. Moreover, as explained in
Sec. II, it was not shown that the regular field defined in Ref. [46] is the
one that enters into the equations of motion.

Here, we instead adopt the Detweiler-Whiting [49] approach to definining a
singular-regular split. While it has not been rigorously shown that the
resulting regular metric produces the correct force for a spinning body, we
again point out that Harte’s [43] choice of effective metric has been shown to
do so and is consistent with Detewiler and Whiting’s choice through linear
order in the secondary’s mass and spin — to which our calculations are
restricted. The resulting approximated singular field is also consistent with
that of Ref. [46] at least through the orders required for the calculation of
Detweiler’s redshift.

###  VII.1 Detweiler-Whiting singular field

In defining a singular-regular split our main criteria are that the singular
field has the same singular structure as the retarded field in the vicinity of
the secondary’s worldline, and that it must not contribute to the equations of
motion. Detweiler and Whiting [49] have shown that in the non-spinning case an
appropriate singular field can be defined in terms of a Green function
decomposition, which is best understood in the Lorenz gauge. Furthermore, the
singular field they identified has the property that when subtracted from the
retarded field, the residual regular field satisfies the homogeneous Lorenz-
gauge wave equation.

The trace-reversed Detweiler-Whiting singular field is defined by

$\bar{h}_{\alpha\beta}^{\rm S}(x)=4\int G^{\rm
S}_{\alpha\beta\alpha^{\prime}\beta^{\prime}}(x,x^{\prime})T^{\alpha^{\prime}\beta^{\prime}}(x^{\prime})\sqrt{-g^{\prime}}d^{4}x^{\prime},$
(102)

where $G^{\rm S}_{\alpha\beta\alpha^{\prime}\beta^{\prime}}(x,x^{\prime})$ is
the Detweiler-Whiting singular Green function. Within a normal neighbourhood
the singular Green function can be expressed in its Hadamard form [62, 49],

$G_{\alpha\beta\alpha^{\prime}\beta^{\prime}}^{\rm
S}=\frac{1}{2}\Big{[}U_{\alpha\beta\alpha^{\prime}\beta^{\prime}}\delta\left(\sigma\right)+V_{\alpha\beta\alpha^{\prime}\beta^{\prime}}\theta\left(\sigma\right)\Big{]},$
(103)

where $\sigma(x,x^{\prime})$ is the Synge world function, and where
$U(x,x^{\prime})_{aba^{\prime}b^{\prime}}$ and
$V(x,x^{\prime})_{aba^{\prime}b^{\prime}}$ are symmetric bi-tensors. The
singular metric perturbation for a spinning body in the Lorenz gauge can be
expressed in terms of these fundamental bi-tensors by substituting the stress-
energy in Eqs. (9), (10a) and (10b), and the singular Green function in Eq.
(103) into Eq. (102). Considering the mass-monopole and spin-dipole
contributions to the singular field separately,

$\bar{h}_{\alpha\beta}^{\rm S}=\bar{h}_{\alpha\beta}^{\rm
S(\mu)}+\bar{h}_{\alpha\beta}^{\rm S(\chi)},$

where

$\displaystyle\bar{h}_{\alpha\beta}^{\rm S(\mu/\chi)}(x)=4\int G^{\rm
S}_{\alpha\beta\alpha^{\prime}\beta^{\prime}}(x,x^{\prime})T^{(\mu/\chi)\alpha^{\prime}\beta^{\prime}}(x^{\prime})\sqrt{-g^{\prime}}d^{4}x^{\prime},$

the mass-monopole contribution to the singular field is well known [49, 63]
and is given by

$\displaystyle\bar{h}_{\alpha\beta}^{\rm S(\mu)}(x)$
$\displaystyle=2\mu\left.\left[\frac{U\left(x,x^{\prime}\right)_{\alpha\beta\alpha^{\prime}\beta^{\prime}}u^{\alpha^{\prime}}u^{\beta^{\prime}}}{\left|\sigma_{\gamma^{\prime}}u^{\gamma^{\prime}}\right|}\right]\right|_{x^{\prime}=x_{A/R}}$
$\displaystyle\qquad+2\mu\int_{\tau_{R}}^{\tau_{A}}V\left(x,z(\tau^{\prime})\right)_{\alpha\beta\alpha^{\prime}\beta^{\prime}}u^{\alpha^{\prime}}u^{\beta^{\prime}}d\tau^{\prime}.$
(104)

Here, the shorthand $[\cdots]|_{x^{\prime}=x_{A/R}}$ corresponds to
$[\cdots]|_{x^{\prime}=x_{A}}+[\cdots]|_{x^{\prime}=x_{R}}$ where $x_{A}$ is
the advanced point at which the future lightcone of $x$ intersects the
worldline and $x_{R}$ is the retarded point at which the past lightcone of $x$
intersects the worldline. The spin-dipole contribution (derived in Appendix D)
is given by

$\displaystyle\bar{h}_{\alpha\beta}^{\rm S(\chi)}(x)$
$\displaystyle=2\mu^{2}\Bigg{[}\bigg{(}\frac{u^{\rho^{\prime}}\nabla_{\rho^{\prime}}U_{\alpha\beta\alpha^{\prime}\beta^{\prime}}\sigma_{\rho^{\prime}}+U_{\alpha\beta\alpha^{\prime}\beta^{\prime}}u^{\kappa^{\prime}}\sigma_{\rho^{\prime}\kappa^{\prime}}}{\sigma_{\delta^{\prime}}u^{\delta^{\prime}}}$
$\displaystyle-(\nabla_{\rho^{\prime}}U_{\alpha\beta\alpha^{\prime}\beta^{\prime}}+V_{\alpha\beta\alpha^{\prime}\beta^{\prime}}\sigma_{\rho^{\prime}})$
$\displaystyle-\frac{U_{\alpha\beta\alpha^{\prime}\beta^{\prime}}\sigma_{\rho^{\prime}}\sigma_{\gamma^{\prime}\kappa^{\prime}}u^{\gamma^{\prime}}u^{\kappa^{\prime}}}{(\sigma_{\delta^{\prime}}u^{\delta^{\prime}})^{2}}\bigg{)}\frac{u^{(\alpha^{\prime}}\tilde{S}^{\beta^{\prime})\rho^{\prime}}}{|\sigma_{\gamma^{\prime}}u^{\gamma^{\prime}}|}\Bigg{]}\Bigg{|}_{x^{\prime}=x_{A/R}}$
$\displaystyle-2\mu^{2}\int_{\tau_{R}}^{\tau_{A}}\nabla_{\rho^{\prime}}V\left(x,z(\tau^{\prime})\right)_{\alpha\beta\alpha^{\prime}\beta^{\prime}}u^{(\alpha^{\prime}}\tilde{S}^{\beta^{\prime})\rho^{\prime}}d\tau^{\prime},$
(105)

where this result is not yet specialised to a given spacetime, spin alignment
or orbital configuration.

###  VII.2 Tensor harmonic regularisation parameters

The forms of equation (VII.1) and (105) are not yet suitable to subtract the
singular metric perturbation from the retarded metric perturbation on the
worldline, since subtracting infinity from infinity is not well defined.
Instead, we use the mode-sum regularisation approach originally pioneered in
the case of a scalar charged particle in Schwarzschild spacetime in [64, 65].
The idea behind the approach is that while the mode-sum producing the retarded
metric perturbation is singular on the secondary body’s worldline, the
discrete modes in the sum are themselves finite. Thus, one can subtract off
the singular metric perturbation mode-by-mode to leave a regular metric
perturbation for which the sum converges to a finite result.

In fact, in the mode sum approach only an approximation to the singular metric
perturbation is required. The mode decomposed singular field can be
represented as an infinite series of so-called regularisation parameters and
only the first few regularisation parameters need to be included to achieve a
finite mode-sum; the inclusion of successive regularisation parameters
(associated with higher-order approximations to the singular field) only
serves to speed up the rate of convergence of the sum. In the case of the
redshift invariant for a spinning body discussed in Sec. VI, the subtraction
of the first two regularisation parameters is sufficient to render the mode-
sum convergent.

In order to derive mode-sum regularisation parameters, we start from a
suitable local expansion of the singular field. High order expansions of the
singular field for a non-spinning secondary in Schwarzschild spacetime were
produced in [63] and used to derive regularisation parameters for various
self-force quantities. That work was later extended to Kerr spacetime [66] and
then to accelerated bodies in the scalar-field case [67]. Applying the
covariant expansion techniques developed in those earlier works to the
singular field in Eqs. (VII.1) and (105) and imposing the condition
$u_{\alpha}\tilde{S}^{\alpha\beta}=\mathcal{O}(\varepsilon^{2})$ leads to an
approximation for the Detweiler-Whiting singular field given by

$\displaystyle\bar{h}_{\mu\nu}^{\rm S}$
$\displaystyle=\mu^{2}\frac{4g_{\mu}^{\bar{\alpha}}g_{\nu}^{\bar{\beta}}u_{(\bar{\alpha}}\tilde{S}_{\bar{\beta})\bar{\gamma}}\sigma^{\bar{\gamma}}}{\bar{\lambda}^{2}\bar{s}^{3}}+\mu\frac{4g_{\mu}^{\bar{\alpha}}g_{\nu}^{\bar{\beta}}u_{\bar{\alpha}}u_{\bar{\beta}}}{\bar{\lambda}\bar{s}}$
$\displaystyle\qquad+\bar{h}_{\mu\nu}^{\rm
S(0)}+\bar{\lambda}\bar{h}_{\mu\nu}^{\rm
S(1)}+\mathcal{O}(\bar{\lambda}^{2},\varepsilon^{2}),$ (106)

where
$\bar{s}\equiv(g_{\bar{a}\bar{b}}+u_{\bar{a}}u_{\bar{b}})\sigma^{\bar{a}}\sigma^{\bar{b}}$.
Here, $\bar{\lambda}$ is simply used as an order counting parameter measuring
distance from the worldline. The bar over the indices represent evaluation at
$\bar{x}$, an arbitrary point on the worldline. The leading-order spin term,
of order $\bar{\lambda}^{-2}$, has been independently derived using matched
asymptotic expansions [5], and the first subleading term, of order
$\bar{\lambda}^{-1}$, derives purely from the mass-monopole and is already
well-established [63]. The subleading terms, of order $\bar{\lambda}^{0}$ and
$\bar{\lambda}^{1}$, are expressed in covariant form here (see Appendix E) for
the first time, though their equivalents were produced in Ref. [46] in Fermi-
Walker coordinates.

The next step in deriving regularization parameters is to perform a coordinate
series expansion of Eq. (VII.2) and to decompose the result into a basis of
spherical-harmonic modes. As our retarded metric perturbation is decomposed
into a basis of scalar, vector and tensor spherical harmonics, to use the
traditional mode-sum approach of using scalar-harmonic regularisation
parameters would require us to project our tensor harmonic modes onto a basis
of scalar spherical harmonics. To avoid this issue, Wardell and Warburton [68]
derived a tensor-harmonic mode-sum regularisation procedure and applied it to
the non-spinning case. We now follow their methodology to derive tensor-
harmonic regularisation parameters for the case of a spinning body. The
process is technically involved, but follows exactly the procedure described
in Ref. [68] so we only briefly summarise the key results here.

The final form of the tensor-harmonic-mode decomposed singular field is

$h_{\mu\nu}^{\rm
S,\ell}=\pm(2\ell+1)h_{\mu\nu}^{[-1]}+h_{\mu\nu}^{[0]}+\mathcal{O}(\ell^{-2}),$
(107)

where $h_{\mu\nu}^{[-1]}$ and $h_{\mu\nu}^{[0]}$ are the leading two
regularisation parameters. We then obtain the regular metric perturbation via
a mode-sum regularization procedure,

$h_{\mu\nu}^{\rm R}=\sum_{\ell=0}^{\infty}\left[h_{\mu\nu}^{{\rm
ret},\ell}\mp(2\ell+1)h_{\mu\nu}^{[-1]}-h_{\mu\nu}^{[0]}\right].$ (108)

At this point, it is important to point out that since Eqs. (VII.1) and 105
are derived in the Lorenz gauge, the regularisation parameters for the
components of the metric perturbation are, in general, only suitable for self-
force calculations in the Lorenz gauge. It is possible to transform the
regularisation parameters to the Regge-Wheeler gauge [69], but in our case
this is not necessary as we are ultimately interested in computing gauge
invariant quantities and thus may use the parameters calculated in the Lorenz
gauge. In particular, inserting the regularization parameters for the Lorenz
gauge metric perturbation into the expression for the redshift, Eq. (99), we
arrive at a mode-sum formula for the gauge-invariant redshift,

$z=\sum_{\ell=0}^{\infty}\left[z^{{\rm
ret},\ell}\mp(2\ell+1)z^{[-1]}-z^{[0]}\right].$ (109)

with regularization parameters given by

$\displaystyle z^{[-1]}$
$\displaystyle=\begin{dcases}\mu^{2}\chi\frac{M^{1/2}(r_{\Omega}-3M)}{2r_{\Omega}^{5/2}(r_{\Omega}-2M)},&\ell\geq
2,\\\ \mu^{2}\chi\frac{M^{1/2}}{2r_{\Omega}^{5/2}},&\ell=1,\\\
-\mu^{2}\chi\frac{M^{1/2}}{2r_{\Omega}^{5/2}},&\ell=0,\end{dcases}$ (110a)
$\displaystyle z^{[0]}$
$\displaystyle=-\frac{2\mu\mathcal{K}(r_{\Omega}^{2}f_{\Omega}-16Mr_{\Omega}f_{\Omega}\Lambda_{1}+16M^{2}\Lambda_{2})}{\pi
r_{\Omega}^{3/2}(r_{\Omega}-3M)(r_{\Omega}-2M)^{1/2}}$
$\displaystyle+\frac{\mu^{2}\chi}{\pi
M^{1/2}r_{\Omega}^{3}\left(r_{\Omega}-3M\right)(r_{\Omega}-2M)^{1/2}}\times$
$\displaystyle\quad\Big{\\{}M\big{[}2\mathcal{K}\left(2r_{\Omega}-3M\right)\left(4M-r_{\Omega}\right)$
$\displaystyle\qquad\quad+\mathcal{E}\left(9M-5r_{\Omega}\right)\left(2M-r_{\Omega}\right)\big{]}$
$\displaystyle\qquad+8\Lambda_{1}\big{[}\mathcal{E}\left(2M-r_{\Omega}\right)\left(13M^{2}-9Mr_{\Omega}+2r_{\Omega}^{2}\right)$
$\displaystyle\qquad\quad+2\mathcal{K}\left(-14M^{3}+16M^{2}r_{\Omega}-7Mr_{\Omega}^{2}+r_{\Omega}^{3}\right)\big{]}$
$\displaystyle\qquad+16\Lambda_{2}M\big{[}\mathcal{E}\left(11M^{2}-9Mr_{\Omega}+2r_{\Omega}^{2}\right)$
$\displaystyle\qquad\quad-2\mathcal{K}\left(5M^{2}-5Mr_{\Omega}+r_{\Omega}^{2}\right)\big{]}\Big{\\}}.$
(110b)

Here,
$\mathcal{K}\equiv\int_{0}^{\pi/2}\left(1-\frac{M}{r_{\Omega}-2M}\sin^{2}x\right)^{-1/2}dx$
and
$\mathcal{E}\equiv\int_{0}^{\pi/2}\left(1-\frac{M}{r_{\Omega}-2M}\sin^{2}x\right)^{1/2}dx$
are elliptic integrals of the first and second kind, respectively, and we have
introduced $\Lambda_{1}\equiv\frac{\ell(\ell+1)}{(2\ell-1)(2\ell+3)}$ and
$\Lambda_{2}\equiv\frac{(\ell-1)\ell(\ell+1)(\ell+2)}{(2\ell-3)(2\ell-1)(2\ell+3)(2\ell+5)}$.

##  VIII Numerical results

Our numerical results naturally divide into three sections. In Section VIII.1
we demonstrate the flux balance law, showing agreement between the asymptotic
gravitational wave energy fluxes and the local rate of change of energy. In
Section VIII.2 we give results for the Detweiler’s redshift invariant. In
Section VIII.3 we produce gravitational waveforms that are complete at
adiabatic order and include post adiabatic spin effects. Unless otherwise
stated, all numerical results are given adimensionalized in $M$, $\mu$ and
$\chi$.

###  VIII.1 Flux balance

In Table 1, we reproduce the flux balance law results of [39] at fixed
frequency, demonstrating Eqn. 70 by comparing the asymptotic energy flux at
$r_{*}=\pm\infty$ with the rate of change of local energy at the worldline. In
all cases we have summed up to $\ell_{\mathrm{max}}=30$ and have set the
tolerances in our numerical integration of the retarded field equations such
that the local rate of change of energy is accurate in all digits shown. We
find that the flux balance law is satisfied to all significant digits in our
calculation, with the asymptotic fluxes agreeing with the local rate of change
of energy to an absolute accuracy of $10^{-28}$ or better.

$\begin{array}[]{c|ccc|c|c}y&\dot{E}_{0}&\dot{E}^{-}_{\chi}&\dot{E}^{+}_{\chi}&dE/d\tau_{\chi}&\Delta^{\mathrm{rel}}_{\chi}\\\
\hline\cr 0.2&2.79273701868\times 10^{-3}&3.77193403195\times
10^{-7}&-6.10406021099\times 10^{-4}&-9.64540266941\times 10^{-4}&1.210\times
10^{-28}\\\ 0.18&1.46844806236\times 10^{-3}&7.60541476292\times
10^{-8}&-2.60585846715\times 10^{-4}&-3.84100734136\times 10^{-4}&1.029\times
10^{-28}\\\ 0.16&7.46754277822\times 10^{-4}&1.08980506901\times
10^{-8}&-1.05064301974\times 10^{-4}&-1.45682859427\times 10^{-4}&2.6\times
10^{-30}\\\ 0.14&3.58765894169\times 10^{-4}&8.06926440798\times
10^{-10}&-3.89407471318\times 10^{-5}&-5.11306464414\times
10^{-5}&1.7713\times 10^{-28}\\\ 0.12&1.58228153292\times
10^{-4}&-6.53905197867\times 10^{-11}&-1.28067951208\times
10^{-5}&-1.60085756392\times 10^{-5}&1.144\times 10^{-29}\\\
0.1&6.15163167846\times 10^{-5}&-2.66993570598\times
10^{-11}&-3.54917559346\times 10^{-6}&-4.24210812069\times 10^{-6}&4.627\times
10^{-29}\\\ 0.09&3.59063362311\times 10^{-5}&-1.01476993750\times
10^{-11}&-1.71031987628\times 10^{-6}&-2.00178988091\times 10^{-6}&4.20\times
10^{-30}\\\ 0.08&1.97579085327\times 10^{-5}&-3.10096178158\times
10^{-12}&-7.62066085171\times 10^{-7}&-8.74153307983\times
10^{-7}&4.4989\times 10^{-29}\\\ 0.07&1.00797672995\times
10^{-5}&-7.55072229184\times 10^{-13}&-3.07211805328\times
10^{-7}&-3.45641134719\times 10^{-7}&1.45815\times 10^{-28}\\\
0.06&4.65287054407\times 10^{-6}&-1.40588119661\times
10^{-13}&-1.08551794346\times 10^{-7}&-1.19875558331\times
10^{-7}&1.11333\times 10^{-28}\\\ 0.05&1.87147091142\times
10^{-6}&-1.85060798132\times 10^{-14}&-3.20089991680\times
10^{-8}&-3.47186542918\times 10^{-8}&4.890\times 10^{-30}\\\
0.04&6.15791960326\times 10^{-7}&-1.49663127140\times
10^{-15}&-7.25545365705\times 10^{-9}&-7.73434118125\times 10^{-9}&2.074\times
10^{-29}\\\ 0.03&1.47265886605\times 10^{-7}&-5.67900033298\times
10^{-17}&-1.08380957000\times 10^{-9}&-1.13614119765\times 10^{-9}&4.583\times
10^{-30}\\\ 0.02&1.96245785614\times 10^{-8}&-5.49135672040\times
10^{-19}&-7.55124235208\times 10^{-11}&-7.78851185422\times
10^{-11}&3.2001\times 10^{-29}\\\ 0.015&4.69335489271\times
10^{-9}&-2.02390121354\times 10^{-20}&-1.14903370686\times
10^{-11}&-1.17579357813\times 10^{-11}&1.0545\times 10^{-28}\\\
0.01&6.23820347340\times 10^{-10}&-1.91947958972\times
10^{-22}&-8.14067891602\times 10^{-13}&-8.26560712092\times
10^{-13}&9.494\times 10^{-30}\\\ \end{array}$

Table 1:  Total non-spin energy flux ($\dot{E}_{0}$), linear-in-spin
contributions to the asymptotic energy flux through the event horizon
($\dot{E}^{-}_{\chi}$) and future null infinity ($\dot{E}^{+}_{\chi}$), local
rate of change of energy, $dE/d\tau_{\chi}$, and relative error in the linear-
in-spin flux balance,
$\Delta^{\mathrm{rel}}_{\chi}\equiv\left|1-\frac{u^{t}\dot{E}_{\chi}}{dE/d\tau_{\chi}}\right|$,
as a function of the frequency of the circular orbit (represented by $y$).

###  VIII.2 Detweiler’s redshift invariant

In Table 2, we give results for Detweiler’s redshift, which is a conservative
gauge invariant quantity and requires regularization. While in Sec. VII.2 we
have only derived the first two regularisation parameters analytically, we can
improve the convergence of the mode-sum by fitting for higher regularisation
parameters (see Fig. 1) whose successive $\ell$ dependence is well known [68].
We can confidently subtract the fitted parameters since we know their
contribution formally sums to zero in the complete mode-sum.

$\begin{array}[]{c|cc|cc}r_{\Omega}&z_{(1)}&z_{(2,\chi)}&\Delta_{z}^{(1)}&\Delta_{z}^{(2,\chi)}\\\
\hline\cr 4&3.0467428778824171300\times 10^{-1}&-1.4586300570011681310\times
10^{-2}&-6.421372315\times 10^{-11}&1.040868132\times 10^{-10}\\\
5&1.8666094967945293038\times 10^{-1}&-1.4359997170303332565\times
10^{-3}&-8.501492777\times 10^{-13}&6.285718934\times 10^{-13}\\\
6&1.4801375464498224547\times 10^{-1}&1.2778851458222201614\times
10^{-4}&-8.494359409\times 10^{-14}&4.205017886\times 10^{-14}\\\
7&1.2619858710413698659\times 10^{-1}&3.2895525026652979198\times
10^{-4}&-1.486953002\times 10^{-14}&6.440450797\times 10^{-15}\\\
8&1.1107483972099400145\times 10^{-1}&3.0130832079882955894\times
10^{-4}&-2.789239601\times 10^{-15}&1.538046174\times 10^{-15}\\\
9&9.9573739278604293268\times 10^{-2}&2.4150088443239778085\times
10^{-4}&-7.542900971\times 10^{-17}&4.914937959\times 10^{-16}\\\
10&9.0385592074434347608\times 10^{-2}&1.8779394202853647999\times
10^{-4}&5.560967658\times 10^{-16}&1.940250420\times 10^{-16}\\\
11&8.2817927552432716594\times 10^{-2}&1.4571779302612887645\times
10^{-4}&6.457897963\times 10^{-16}&8.995489317\times 10^{-17}\\\
12&7.6451679289700237541\times 10^{-2}&1.1390103980806758238\times
10^{-4}&5.909892489\times 10^{-16}&4.722720931\times 10^{-17}\\\
13&7.1010239380097475857\times 10^{-2}&8.9977582894836949305\times
10^{-5}&5.057588184\times 10^{-16}&2.730507042\times 10^{-17}\\\
14&6.6300106250377251019\times 10^{-2}&7.1894261756329048211\times
10^{-5}&4.235268423\times 10^{-16}&1.701375794\times 10^{-17}\\\
15&6.2180255400520542626\times 10^{-2}&5.8095022692680346999\times
10^{-5}&3.528715259\times 10^{-16}&1.123780963\times 10^{-17}\\\
16&5.8544734161342785070\times 10^{-2}&4.7448661525134242815\times
10^{-5}&2.945648894\times 10^{-16}&7.770696488\times 10^{-18}\\\
17&5.5312035461636850090\times 10^{-2}&3.9141832229418841873\times
10^{-5}&2.471167798\times 10^{-16}&5.572694481\times 10^{-18}\\\
18&5.2418249523532483206\times 10^{-2}&3.2588897056415252246\times
10^{-5}&2.086117042\times 10^{-16}&4.115745383\times 10^{-18}\\\
19&4.9812462417061066657\times 10^{-2}&2.7365239308320767569\times
10^{-5}&1.772872625\times 10^{-16}&3.113963219\times 10^{-18}\\\
20&4.7453560812119756566\times 10^{-2}&2.3160083573548132458\times
10^{-5}&1.516776702\times 10^{-16}&2.403886377\times 10^{-18}\\\
25&3.8367861668465180133\times 10^{-2}&1.1094904818711662769\times
10^{-5}&7.613023266\times 10^{-17}&8.160586074\times 10^{-19}\\\
30&3.2200482214638454294\times 10^{-2}&6.0251503215045887450\times
10^{-6}&4.312790601\times 10^{-17}&3.486595670\times 10^{-19}\\\
35&2.7740018849300996376\times 10^{-2}&3.5788991032583934219\times
10^{-6}&2.666889209\times 10^{-17}&1.718004154\times 10^{-19}\\\
40&2.4364201607676976872\times 10^{-2}&2.2730939533525315130\times
10^{-6}&1.759954961\times 10^{-17}&9.349372617\times 10^{-20}\\\
45&2.1720433941951952235\times 10^{-2}&1.5205343987439192133\times
10^{-6}&1.220841693\times 10^{-17}&5.478475527\times 10^{-20}\\\
50&1.9593977138759696222\times 10^{-2}&1.0599790877771606396\times
10^{-6}&8.808541368\times 10^{-18}&3.400305245\times 10^{-20}\\\
55&1.7846580120368287542\times 10^{-2}&7.6418367844689967876\times
10^{-7}&6.560715313\times 10^{-18}&2.210127995\times 10^{-20}\\\
60&1.6385213628045196549\times 10^{-2}&5.6652220737162069625\times
10^{-7}&5.016155327\times 10^{-18}&1.492031137\times 10^{-20}\\\
65&1.5144980956649394828\times 10^{-2}&4.2998668462990765279\times
10^{-7}&3.920290881\times 10^{-18}&1.039710466\times 10^{-20}\\\
70&1.4079232946219360584\times 10^{-2}&3.3298551248965351373\times
10^{-7}&3.121481287\times 10^{-18}&7.443015862\times 10^{-21}\\\
75&1.3153576784133680096\times 10^{-2}&2.6238746624051624458\times
10^{-7}&2.525601325\times 10^{-18}&5.453231584\times 10^{-21}\\\
80&1.2342099554180064613\times 10^{-2}&2.0991759354724383914\times
10^{-7}&2.072137298\times 10^{-18}&4.076764659\times 10^{-21}\\\
85&1.1624906883132949330\times 10^{-2}&1.7019819208348517223\times
10^{-7}&1.720986920\times 10^{-18}&3.102161030\times 10^{-21}\\\
90&1.0986471753831045984\times 10^{-2}&1.3963904283654329585\times
10^{-7}&1.444862734\times 10^{-18}&2.397815109\times 10^{-21}\\\
95&1.0414498580229646206\times 10^{-2}&1.1578567408369181383\times
10^{-7}&1.224760538\times 10^{-18}&1.879447870\times 10^{-21}\\\
100&9.8991242481267771860\times 10^{-3}&9.6924289089700379549\times
10^{-8}&1.047169541\times 10^{-18}&1.491706402\times 10^{-21}\\\ \end{array}$

Table 2:  Numerical results for the non-spin ($z_{(1)}$) and the linear-in-
spin ($z_{(2,\chi)}$) contributions to the redshift invariant. The error on
$z_{(1)}$ and $z_{(2,\chi)}$ are quantified by $\Delta_{z}^{(1)}$ and
$\Delta_{z}^{(2,\chi)}$ respectively, which correspond to the error introduced
for truncating the mode-sum at $\ell_{\mathrm{max}}$. In our calculations the
truncation error was always greater than the numerical error — we summed up to
$\ell_{\mathrm{max}}=40$ and included the first seven regularisation
parameters.

Strong field values of $z_{(1)}$ were computed in Refs. [70, 69] (although
both give the values for $-(u^{t})^{2}z_{(1)}$ and truncate the mode-sum at a
larger $\ell_{\mathrm{max}}$); these agree with our non-spinning result to
within the errors given. Strong field values of the linear in spin
contribution to the redshift, $z_{(2,\chi)}$, have not been previously
computed, although a post-Newtonian expansion was produced in Ref. [22]. In
Fig. 2, we compare our numerical results for $z_{(2,\chi)}$ against the
equivalent post-Newtonian series from Ref. [22]. As expected, the absolute
error between our fully relativistic numerical results and the PN expanded
redshift is higher for large values of $y$ (in the strong field) where the PN
expansion breaks down. The error improves significantly further towards the
weak field and improves again when adding higher order terms in the PN
expansion. In the weak field, the leading contribution to the residual goes as
the next (unknown) PN term beyond the order where the series was truncated.

Figure 1: An example of the mode-sum regularisation procedure used to
calculate the linear-in-spin part of the redshift, $z_{(2,\chi)}$, for
$r_{\Omega}=100M$. The data points show the absolute value of the successively
regularized $\ell$ modes of $z_{(2,\chi)}$ when using different numbers of
regularization parameters; the joined lines are the fitted regularisation
parameters. The leading behaviour in $\ell$ of the regularized $z_{(2,\chi)}$
modes goes as the next leading regularisation parameter that was not included
in the regularisation scheme. Hence the joined lines overlay the plotted
points and qualitatively verify the fitted parameters. Figure 2: The absolute
error $\Delta
z_{(2,\chi)}\equiv\left|z_{(2,\chi)}-z_{(2,\chi)}^{\mathrm{PN}}\right|$
between our numerical results and the PN series on a logarithmic scale as a
function of $y$. Each set of plot markers show the comparison including
increasingly high order PN terms, with the highest included order in $y$
labelled in the plot legend. The joined lines give the individual $y^{5.5}$,
$y^{7.5}$, $y^{8.5}$ and $y^{9.5}$ terms in the PN series; towards the weak
field these approach the leading contributions to the residuals, $\Delta
z_{(2,\chi)}$, when the PN series is truncated to the preceding order.

###  VIII.3 Waveforms

To produce a gravitational waveform incorporating effects from a spinning
secondary, we follow the two-timescale evolution and waveform-generation
procedure summarised in Sec. IV.5. On a grid of $\Omega$ values, we solve the
RWZ equations for the master functions and calculate $C_{\ell
m\omega}^{\pm}(\Omega)$ in doing so. From the $C_{\ell m\omega}^{\pm}$ values
we obtain the asymptotic energy fluxes via Eq. (90) (summing to some
$\ell_{\rm max}$ determined by the accuracy requirements) and in turn produce
the forcing functions, $F_{\Omega}^{(0)}(\Omega)$ and
$F_{\Omega}^{(1,\chi)}(\Omega)$ via Eqs. (71) and (74). We then fix the value
of $\mu\chi$ and solve for the orbital phase $\phi_{p}(t)$ and frequency
$\Omega(t)$.

To obtain a waveform, we note that the gravitational wave strain in the RWZ
formalism is given by [5]

$r\left(h_{+}-ih_{\times}\right)=\sum_{\ell=2}^{\infty}\sum_{m=-\ell}^{\ell}\frac{D}{2}\left(\psi_{\ell
m}^{\mathrm{even}}-i\psi_{\ell m}^{\mathrm{odd}}\right)_{-2}Y_{\ell
m}(\theta,\phi),$ (111)

where the equality holds in the limit $r\rightarrow\infty$ and the constant
$D\equiv\sqrt{(\ell-1)\ell(\ell+1)(\ell+2)}$. The function ${}_{-2}Y_{\ell
m}(\theta,\phi)$ is the spin-weight $s=-2$ spherical harmonic. Conveniently,
$\psi_{\ell m}^{\mathrm{even/odd}}\rightarrow C_{\ell
m\omega}^{+\mathrm{even/odd}}e^{-im\phi_{p}(t)}$ as $r\rightarrow\infty$, so
defining $h\equiv r\left(h_{+}-ih_{\times}\right)$ and explicitly taking the
limit in Eq. (111) we can write

$h=\sum_{\ell=2}^{\infty}\sum_{m=-\ell}^{\ell}h_{\ell m}(t)\,_{-2}Y_{\ell
m}(\theta,\phi),$ (112)

where the spherical harmonic modes of the waveform are given by

$\displaystyle h_{\ell m}(t)$ $\displaystyle=\begin{dcases}\frac{D}{2}C_{\ell
m\omega}^{+}\left(\Omega(t)\right)e^{-im\phi_{p}(t)}&\ell+m\quad\text{even},\\\
-\frac{iD}{2}C_{\ell
m\omega}^{+}\left(\Omega(t)\right)e^{-im\phi_{p}(t)}&\ell+m\quad\text{odd}.\end{dcases}$
(113)

In Figures 3 and 4 we plot the $(\ell,m)=(2,2)$ mode of the waveform for mass
ratio $1:1$ and $1:10^{5}$ binaries, respectively. Figures 5 and 6 demonstrate
the impact of the spin on the waveform’s phase. The leading phase difference
between a waveform including the secondary’s spin and a non-spinning waveform
comes from the linear in spin contribution to
$\tilde{\Delta\phi}=\tilde{\phi}_{0}-\tilde{\phi}\sim-\varepsilon\chi\tilde{\phi}_{\chi}=-\chi\phi_{\chi}.$
Thus the leading phase difference is independent of the mass ratio and Fig. 6
can be rescaled for any mass ratio and $\chi$ by simply multiplying
$\tilde{\phi}$ by $\chi$ and $\tilde{t}$ by $\varepsilon^{-1}$, giving the
error to expect on the waveform’s phase accumulated over time if neglecting
the spin of the secondary. This highlights the importance of including the
secondary’s spin in waveform models, as a key requirement to test fundamental
physics with EMRI waveforms is to accurately track the phase over the full
inspiral.

Figure 3: Comparison of the $(\ell,m)=(2,2)$ mode of the waveform with
$\chi=1$ (blue) and $\chi=0$ (orange) for an equal mass binary. Both inspirals
begin at the same orbital frequency ($r_{\Omega}=20M$) and phase before
evolving. Figure 4: Comparison as in Fig. 3 but for a binary with a more
realistic EMRI mass ratio of $1:10^{5}$.  Figure 5: Comparison of the waveform
phase, $\tilde{\phi}$, as a function of orbital frequency with $\chi=1$ (blue)
and $\chi=0$ (orange) for an equal mass binary. Figure 6: Phase difference (in
radians) between the waveform with a spinning seocondary and with a non-
spinning secondary, $\Delta\tilde{\phi}=\tilde{\phi}_{0}-\tilde{\phi}$.

##  IX Conclusions

We have performed calculations of the leading order self-force including the
sub-leading linear-in-spin effects in the two-timescale expansion and
demonstrated how these spin effects can be included in 1PA waveforms. In doing
so we have produced waveforms complete at adiabatic order and including the
spin’s 1PA contribution. We have developed a fixed frequency approach for
solving the perturbation equations which significantly improves the
computational efficiency of linear-in-spin self-force quantities, including
the asymptotic gravitational wave fluxes and the metric perturbation. We
provided the first fully relativistic regularisation scheme to treat the
singular field associated with a spinning secondary body, using a first
principles approach — having derived the spin-dipole contributions to the
Detweiler-Whiting singular field and produced a covariant expansion of the
approximated field. We have performed the first fully relativistic
calculations of a conservative self-force quantity with a spinning secondary —
having computed the redshift invariant and checked agreement with the
equivalent PN expansion.

There are several obvious future extensions to this work. Firstly, we could
repeat our calculations for more general EMRI configurations, allowing for
eccentric equatorial or generic motion of the secondary body and performing
the computations in a Kerr background spacetime to include the spin of the
primary. Further generalization beyond equatorial motion requires treatment of
the secondary body’s precessing spin (instead of the spin aligned/anti-aligned
case). Now that the regularisation scheme has been developed for a spinning
body, we could extend current calculations of other gauge invariants such as
the spin precession invariant to include the spin of the secondary. We could
also study more extended body effects, beginning with including the
secondary’s quadrupole moment in our computations. Finally, as progress on 2SF
calculations continue, it would be natural to compute complete 1PA waveforms,
including all the conservative 1SF and dissipative 2SF effects as well as the
spin effects highlighted in this work. Aside from their immediate application
for LISA, complete 1PA waveforms could be compared with numerical relativity
counterparts and used to establish the domain of validity of perturbation
theory in less extreme mass ratios.

###### Acknowledgements.

The authors thank Gabriel Piovano, Vojtech Witzany and Chris Kavanagh for
helpful discussions. JM acknowledges support from the Irish Research Council
under grant GOIPG/2018/448. AP acknowledges the support of a Royal Society
University Research Fellowship. This work makes use of the Black Hole
Perturbation Toolkit [59] and xAct [71].

##  Appendix A Stress-energy tensor for a point mass with aligned spin on a
circular orbit

The coefficients appearing in Eq. (26) are given explicitly by

$\displaystyle K_{0}^{tt}=u^{t},\quad K_{0}^{t\phi}=u^{\phi},\quad
K_{0}^{\phi\phi}=\left(u^{\phi}\right)^{2}/u^{t},$ $\displaystyle
K_{1}^{tt}=\frac{-M^{3/2}}{r_{\Omega}\left(r_{\Omega}-2M\right)\sqrt{r_{\Omega}-3M}},\quad
K_{1}^{t\phi}=K_{1}^{\phi t}=-\frac{M}{r_{\Omega}^{5/2}\sqrt{r_{\Omega}-3M}},$
$\displaystyle
K_{1}^{rr}=-\frac{M^{1/2}\left(r_{\Omega}-2M\right)\sqrt{r_{\Omega}-3M}}{r_{\Omega}^{3}},\quad
K_{1}^{\phi\phi}=-\frac{M^{1/2}\left(r_{\Omega}-2M\right)}{r_{\Omega}^{4}\sqrt{r_{0}-3M}},$
$\displaystyle
K_{2}^{tr}=K_{2}^{rt}=\frac{\sqrt{r_{\Omega}-3M}}{2r_{\Omega}^{3/2}},\quad
K_{2}^{r\phi}=K_{2}^{\phi
r}=\frac{M^{1/2}\sqrt{r_{\Omega}-3M}}{2r_{\Omega}^{3}},$ $\displaystyle
K_{3}^{tt}=-\frac{M^{1/2}}{\sqrt{r_{\Omega}-3M}},\quad
K_{3}^{t\phi}=K_{3}^{\phi
t}=-\frac{(r_{\Omega}-M)}{2r_{\Omega}^{3/2}\sqrt{r_{\Omega}-3M}},\quad
K_{3}^{\phi\phi}=-\frac{M^{1/2}\left(r_{\Omega}-2M\right)}{r_{\Omega}^{3}\sqrt{r_{\Omega}-3M}}.$
(114)

##  Appendix B Sources for the RWZ equations at fixed frequency

###  B.1 Odd parity

The Cunningham-Price-Moncrief master function is defined in [55] and satisfies
the RWZ equation with the Regge-Wheeler potential;

$V_{\text{odd}}\equiv\frac{f}{r^{2}}\left[\ell(\ell+1)-\frac{6M}{r}\right].$
(115)

For a spin-aligned secondary in a circular equatorial orbit parameterised at
fixed frequency, the frequency domain source (76) in the odd parity sector has
the following coefficients of the distributional functions:

$\displaystyle H^{(2,\chi)}_{\ell m\omega}$
$\displaystyle=-\frac{4\pi\mu^{2}\chi
f_{\Omega}^{2}}{\lambda(\lambda+1)}\sqrt{1-\frac{3M}{r_{\Omega}}}\sqrt{(\ell-m)(1+\ell+m)}Y_{\ell(m+1)}\left(\frac{\pi}{2},0\right),$
$\displaystyle F^{(1)}_{\ell m\omega}$ $\displaystyle=\frac{8\pi\mu
f_{\Omega}^{2}}{\lambda(\lambda+1)}\sqrt{\frac{M}{r_{\Omega}-3M}}\sqrt{(\ell-m)(1+\ell+m)}Y_{\ell(m+1)}\left(\frac{\pi}{2},0\right),$
$\displaystyle F^{(2,\chi)}_{\ell m\omega}$
$\displaystyle=-\frac{8\pi\mu^{2}\chi
f_{\Omega}}{\lambda(\lambda+1)}\frac{\left(M^{2}+3Mr_{\Omega}-r_{\Omega}^{2}\right)}{r_{\Omega}^{3}}\sqrt{\frac{r_{\Omega}}{r_{\Omega}-3M}}\sqrt{(\ell-m)(1+\ell+m)}Y_{\ell(m+1)}\left(\frac{\pi}{2},0\right),$
$\displaystyle G^{(1)}_{\ell m\omega}$ $\displaystyle=-\frac{8\pi\mu
f_{\Omega}}{\lambda(\lambda+1)r_{\Omega}}\sqrt{\frac{M}{r_{\Omega}-3M}}\sqrt{(\ell-m)(1+\ell+m)}Y_{\ell(m+1)}\left(\frac{\pi}{2},0\right),$
$\displaystyle G^{(2,\chi)}_{\ell m\omega}$
$\displaystyle=\frac{4\pi\mu^{2}\chi
M\left(m^{2}r_{\Omega}+(2-3m^{2})M\right)}{\lambda(\lambda+1)r_{\Omega}^{4}}\sqrt{\frac{r_{\Omega}}{r_{\Omega}-3M}}\sqrt{(\ell-m)(1+\ell+m)}Y_{\ell(m+1)}\left(\frac{\pi}{2},0\right),$

where $\lambda\equiv(\ell-1)(\ell+2)/2$ and we have separated the spin
independent and linear-in-spin contributions as $\bar{F}_{\ell m\omega}\equiv
F^{(1)}_{\ell m\omega}+F^{(2,\chi)}_{\ell m\omega}$, $\bar{G}_{\ell
m\omega}\equiv G^{(1)}_{\ell m\omega}+G^{(2,\chi)}_{\ell m\omega}$ and
$\bar{H}_{\ell m\omega}\equiv H^{(2,\chi)}_{\ell m\omega}$.

###  B.2 Even parity

The Zerilli-Moncrief master function is defined in [55] and satisfies the RWZ
equation with the Zerilli potential;

$V_{\text{even}}\equiv\frac{f}{r^{2}\Lambda^{2}}\left[2\lambda^{2}\left(\lambda+1+\frac{3M}{r}\right)+\frac{18M^{2}}{r^{2}}\left(\lambda+\frac{M}{r}\right)\right].$
(116)

For a spin-aligned secondary in a circular equatorial orbit parameterised at
fixed frequency, the frequency domain source (76) in the even parity sector
has the following coefficients of the distributional functions:

$\displaystyle H^{(2,\chi)}_{\ell m\omega}=$ $\displaystyle 0,$ $\displaystyle
F^{(1)}_{\ell m\omega}=$ $\displaystyle\frac{8\pi\mu
r_{\Omega}f_{\Omega}^{3}}{(\lambda
r_{\Omega}+3M)(\lambda+1)}\sqrt{\frac{r_{\Omega}}{r_{\Omega}-3M}}Y_{\ell
m}\left(\frac{\pi}{2},0\right),$ $\displaystyle F^{(2,\chi)}_{\ell m\omega}=$
$\displaystyle-\frac{8\pi\mu^{2}\chi f_{\Omega}(r_{\Omega}-3M)}{(\lambda
r_{\Omega}+3M)(\lambda+1)r_{\Omega}^{2}}\sqrt{\frac{M}{r_{\Omega}-3M}}Y_{\ell
m}\left(\frac{\pi}{2},0\right)\left(3\left(1-m^{2}\right)\frac{M}{\lambda}+r_{\Omega}\left(\lambda-m^{2}\right)+\frac{\left(-7M^{2}-Mr_{\Omega}+r_{\Omega}^{2}\right)}{(r_{\Omega}-3M)}\right),$
$\displaystyle G^{(1)}_{\ell m\omega}=$ $\displaystyle\frac{8\pi\mu
f_{\Omega}}{(\lambda+1)r_{\Omega}^{2}(3M+\lambda
r_{\Omega})^{2}}\sqrt{\frac{r_{\Omega}}{r_{\Omega}-3M}}Y_{\ell
m}\left(\frac{\pi}{2},0\right)\left(3M^{3}\left(5+(3/\lambda)(m^{2}-1)\right)+2M^{2}r_{\Omega}\left(2\lambda+3m^{2}-9\right)\right.$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\left.+\lambda
Mr_{\Omega}^{2}\left(\lambda+m^{2}-4\right)-\lambda(\lambda+1)r_{\Omega}^{3}\right),$
$\displaystyle G^{(2,\chi)}_{\ell m\omega}=$
$\displaystyle-\frac{8\pi\mu^{2}\chi}{\lambda(\lambda+1)r_{\Omega}^{4}(3M+\lambda
r_{\Omega})^{2}}\sqrt{\frac{M}{r_{\Omega}-3M}}Y_{\ell
m}\left(\frac{\pi}{2},0\right)\left(3\lambda
m^{2}Mr_{\Omega}\left(2M^{2}+3Mr_{\Omega}-r_{\Omega}^{2}\right)\right.$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\left.+(4M-r_{\Omega})\left(9\left(m^{2}-1\right)M^{3}-3\lambda
M\left(M^{2}+Mr_{\Omega}+r_{\Omega}^{2}\right)+\lambda^{3}r_{\Omega}^{2}(M-r_{\Omega})\right)\right.$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad\left.+\lambda^{2}r_{\Omega}\left(-m^{2}r_{\Omega}\left(2M^{2}-4Mr_{\Omega}+r_{\Omega}^{2}\right)+4M^{3}-10M^{2}r_{\Omega}-2Mr_{\Omega}^{2}+r_{\Omega}^{3}\right)\right).$

##  Appendix C Fixed-radius approaches to numerical RWZ solutions

Earlier works that parameterise the secondary body at fixed radius extract out
the linear in spin terms from numerical results by repeating the calculations
for multiple values of $\chi$ and fitting a polynomial to the results [39,
40]. This method is undesirable as it is very computationally expensive. It is
also possible to linearize the perturbation equations themselves and solve a
coupled system with an extended source (as in Ref. [41]). That way, numerical
results relying on the integration of the perturbation equations may be
obtained directly to linear order in spin. In the RWZ formalism, the
linearized perturbation equations in the frequency domain are

$\displaystyle\mathcal{L}_{0}\psi^{0}_{\ell m\omega}(r)=Z^{0}_{\ell
m\omega}(r),$ (117) $\displaystyle\mathcal{L}_{0}\psi^{\chi}_{\ell
m\omega}(r)=Z^{\chi}_{\ell m\omega}(r)-2\omega_{0}\omega_{\chi}\psi^{0}_{\ell
m\omega}(r),$ (118)

where we have defined the RWZ ‘non-spin’ operator as:

$\mathcal{L}_{0}\equiv\left[\dfrac{\partial^{2}}{\partial
r_{*}^{2}}-V_{\ell}(r)+\omega_{0}^{2}\right],$

and separated the spin independent and linear-in-spin terms as $\psi_{\ell
m\omega}(r)=\psi^{0}_{\ell m\omega}(r)+\psi^{\chi}_{\ell m\omega}(r)$,
$Z_{\ell m\omega}(r)=Z^{0}_{\ell m\omega}(r)+Z^{\chi}_{\ell m\omega}(r)$ and
$\omega=\omega_{0}+\omega_{\chi}$. Eq. (117) is the usual RWZ equation for a
non-spinning secondary body. Eq. (118) for the linear in spin master function
has the same differential operator as Eq. (117) but sourced by both the
compact spin source $Z^{\chi}_{\ell m\omega}(r)$ and an extended source term
from the coupling to $\psi^{0}_{\ell m\omega}(r)$.

In general, we wish to avoid modelling with an extended source. This becomes
especially important when generalising to eccentric (or more generic) motion
of the secondary and using the method of extended homogeneous solutions (EHS)
to solve the perturbation equations [72, 56] instead of regular variation of
parameters. Motivated by this, we split the extended and compact parts of Eq.
(117) into two separate equations — by linearity we can write
$\psi^{\chi}_{\ell m\omega}(r)=\psi^{\chi,\rm ext}_{\ell
m\omega}(r)+\psi^{\chi,\rm c}_{\ell m\omega}(r)$. Acting on the extended
equation again, having identified the RWZ operator as a ‘partial annihilator’
[73], we obtain

$\displaystyle\mathcal{L}_{0}\psi^{\chi,\rm c}_{\ell
m\omega}(r)=Z^{\chi}_{\ell m\omega}(r),$ (119)
$\displaystyle\mathcal{L}_{0}^{2}\psi^{\chi,\rm ext}_{\ell
m\omega}(r)=-2\omega_{0}\omega_{\chi}Z^{0}_{\ell m\omega}(r),$ (120)

where now the entire system involves only compact sources and we have
decoupled the equations. The trade off is having to solve the fourth order Eq.
(120) which has four homogeneous solutions to find numerically. Conveniently,
only two new homogeneous solutions are required in practice as the other two
are shared with Eq. (117) — the ‘in’ and ‘up’ solutions. The two new
homogeneous solutions to Eq. (120) can be considered the extended ‘in’ and
‘up’ solutions and the boundary conditions are similarly posed by considering
their asymptotic behaviour

$\displaystyle\hat{R}_{\ell m\omega}^{\chi,\rm
ext-}(r_{*}\rightarrow-\infty)\sim r_{*}e^{-i\omega_{0}r_{*}},$ (121)
$\displaystyle\hat{R}_{\ell m\omega}^{\chi,\rm
ext+}(r_{*}\rightarrow\infty)\sim r_{*}e^{i\omega_{0}r_{*}}.$ (122)

Once all the homogeneous solutions are obtained numerically, the general
solutions for the spin-linearized master functions are easily obtained via
standard variation of parameters. The method can be extended to eccentric
orbits using EHS.

However, the issue of linearizing numerical quantities such as the master
functions is entirely avoided with the fixed frequency parametrisation where
$\omega_{\chi}=0$. Then, the homogeneous RWZ equation is independent of spin
and may be solved in the same way as for a non-spinning secondary. The
linearized inhomogeneous solution for the retarded master function is obtained
analytically from the homogeneous solutions as described in Section V.2 in
Eqs. (82, V.2) — in practice this is equivalent to solving Eqs. (117) and
(119) with $\psi^{\chi,\rm ext}_{\ell m\omega}(r)=0$.

##  Appendix D Derivation of the spin-dipole singular field

Substituting the stress-energy in Eq.(10b) and the singular Green function in
Eq. (103) into Eq. (102),

$\displaystyle\bar{h}_{\alpha\beta}^{\rm S(\chi)}(x)$ $\displaystyle=4\int
G^{\rm
S}_{\alpha\beta\alpha^{\prime}\beta^{\prime}}(x,x^{\prime})\nabla_{\rho^{\prime}}\left(\int
d\tau^{\prime}\,\frac{\delta^{4}\left[x^{\prime}-z(\tau^{\prime})\right]}{\sqrt{-g^{\prime}}}u^{(\alpha^{\prime}}\tilde{S}^{\beta^{\prime})\rho^{\prime}}\right)\sqrt{-g^{\prime}}d^{4}x^{\prime}$
$\displaystyle=-2\int\delta[\sigma]\left[\nabla_{\rho^{\prime}}U_{\alpha\beta\alpha^{\prime}\beta^{\prime}}+V_{\alpha\beta\alpha^{\prime}\beta^{\prime}}\sigma_{\rho^{\prime}}\right]u^{(\alpha^{\prime}}\tilde{S}^{\beta^{\prime})\rho^{\prime}}d\tau^{\prime}-2\int\delta^{\prime}\left[\sigma\right]U_{\alpha\beta\alpha^{\prime}\beta^{\prime}}\sigma_{\rho^{\prime}}u^{(\alpha^{\prime}}\tilde{S}^{\beta^{\prime})\rho^{\prime}}d\tau^{\prime}$
$\displaystyle\quad-2\int\theta\left[\sigma\right]\nabla_{\rho^{\prime}}V_{\alpha\beta\alpha^{\prime}\beta^{\prime}}u^{(\alpha^{\prime}}\tilde{S}^{\beta^{\prime})\rho^{\prime}}d\tau^{\prime},$
(123)

where in this context $\delta[\sigma]=\delta[\sigma(x,z(\tau^{\prime}))]$.
Considering the three different distributional integrals separately, the first
integral is

$\displaystyle-2\int\left(\frac{\delta\left[\tau^{\prime}-\tau_{A}\right]}{\left|\sigma_{c^{\prime}}u^{c^{\prime}}\right|}+\frac{\delta\left[\tau^{\prime}-\tau_{R}\right]}{\left|\sigma_{c^{\prime}}u^{c^{\prime}}\right|}\right)\left(\nabla_{\rho^{\prime}}U_{\alpha\beta
a^{\prime}b^{\prime}}+V_{\alpha\beta
a^{\prime}b^{\prime}}\sigma_{\rho^{\prime}}\right)u^{(\alpha^{\prime}}\tilde{S}^{b^{\prime})\rho^{\prime}}d\tau^{\prime}$
$\displaystyle=-2\left.\left[\frac{\left(\nabla_{\rho^{\prime}}U_{\alpha\beta
a^{\prime}b^{\prime}}+V_{\alpha\beta
a^{\prime}b^{\prime}}\sigma_{\rho^{\prime}}\right)u^{(\alpha^{\prime}}\tilde{S}^{b^{\prime})\rho^{\prime}}}{\left|\sigma_{c^{\prime}}u^{c^{\prime}}\right|}\right]\right|_{x^{\prime}=x_{A/R}},$
(124)

where $\delta\left[\sigma\right]$ was rewritten using equation (132), the fact
that the two simple roots of $\sigma\left(x,z(\tau^{\prime})\right)$ along the
worldline are $x_{A}/x_{R}$, and that
$\frac{d}{d\tau^{\prime}}\sigma(x,z(\tau^{\prime}))=\sigma_{c^{\prime}}u^{c^{\prime}}$.
The second integral is

$\displaystyle-2\int\left(\frac{\delta^{\prime}\left[\tau^{\prime}-\tau_{A}\right]+\delta^{\prime}\left[\tau^{\prime}-\tau_{R}\right]}{\left|\sigma_{c^{\prime}}u^{c^{\prime}}\right|\sigma_{d^{\prime}}u^{d^{\prime}}}+\frac{\left(\delta\left[\tau^{\prime}-\tau_{A}\right]+\delta\left[\tau^{\prime}-\tau_{A}\right]\right)\sigma_{p^{\prime}q^{\prime}}u^{p^{\prime}}u^{q^{\prime}}}{\left|\sigma_{c^{\prime}}u^{c^{\prime}}\right|(\sigma_{d^{\prime}}u^{d^{\prime}})^{2}}\right)U_{\alpha\beta
a^{\prime}b^{\prime}}\sigma_{\rho^{\prime}}u^{(\alpha^{\prime}}\tilde{S}^{b^{\prime})\rho^{\prime}}d\tau^{\prime}$
$\displaystyle=2\left.\left[\left(\frac{u^{p^{\prime}}\nabla_{p^{\prime}}U_{\alpha\beta
a^{\prime}b^{\prime}}\sigma_{\rho^{\prime}}+U_{\alpha\beta
a^{\prime}b^{\prime}}u^{q^{\prime}}\sigma_{\rho^{\prime}q^{\prime}}}{\left|\sigma_{c^{\prime}}u^{c^{\prime}}\right|\sigma_{d^{\prime}}u^{d^{\prime}}}-\frac{U_{\alpha\beta
a^{\prime}b^{\prime}}\sigma_{\rho^{\prime}}\sigma_{p^{\prime}q^{\prime}}u^{p^{\prime}}u^{q^{\prime}}}{\left|\sigma_{c^{\prime}}u^{c^{\prime}}\right|(\sigma_{d^{\prime}}u^{d^{\prime}})^{2}}\right)u^{(\alpha^{\prime}}\tilde{S}^{b^{\prime})\rho^{\prime}}\right]\right|_{x^{\prime}=x_{A/R}},$
(125)

where the $\delta^{\prime}\left[\sigma\right]$ has been rewritten using
equation (133). Also, although strictly $\frac{d^{2}}{d\tau^{\prime