For reference, the following erratum corrects the published version of the paper. These errors have been fixed in
this arxiv-version (the article starting on page 2 has the corrected expressions).
Erratum: Evolution of the Carter constant for inspirals into a black hole: Effect of the
black hole quadrupole
[Phys. Rev. D 75, 124007 (2007)]
Éanna É. Flanagan, Tanja Hinderer
In Eqs. (3.16), (3.17), (3.18), (3.24), (3.25) and (3.26) of this paper, the variable r should be replaced everywhere
by the variable r̃, and the variable θ should be replaced everywhere by the variable θ̃. The definitions of r̃ and θ̃ are
given in Eq. (2.11). These replacements do not affect the any of the subsequent results in the paper.
Also, the right hand side of Eq. (B3) is missing a term −4SLzr̃ and Eq. (2.24) is missing a factor of dϕ/dt̃ in front
of Q.
Some terms are missing in Eqs. (3.18), (3.26) and (3.30) - (3.33). The additional terms in Eqs. (3.18) and (3.26)
15r̃7
−75K2 + 2Kr̃(51r̃E + 50) + 8r̃2(r̃E + 1)(3r̃E + 5)
15p2r̃7
25p3(3p− 4r̃) + p2r̃2
11− 51e2
+ 32pr̃3
1− e2
+ 6r̃4
1− e2
respectively. These result in additional fractional corrections to Eq. (3.30) given by
and the full expression replacing the O(Q) terms in Eq. (3.30) is then
〈K̇〉 = −
(1− e2)3/2
cos(2ι)
+O(S), O(S2)− terms.
Equations (3.31), (3.32) and (3.33) contain typos in the O(S) and O(Q) terms, the corrected expressions are given
below. We thank P. Komorowski for pointing this out. Equation (3.31) should be replaced by
〈ṗ〉 = −64
(1− e2)3/2
− S cos(ι)
96p3/2
1064 + 1516e2 + 475e4
149e2
469e2
227e4
cos(2ι)
+ e2 +
[13− cos(2ι)]
, (0.1)
Equation (3.32) should be replaced by
〈ė〉 = −304
e(1− e2)3/2
121e2
Se(1− e2)3/2 cos(ι)
5p11/2
1172 + 932e2 +
1313e4
Q(1− e2)3/2
785e2
− 219e
+ 13e6 +
2195e2
+ 251e4 +
218e6
cos(2ι)
2e(1− e2)3/2
2 + 3e2 +
[13− cos(2ι)] , (0.2)
and the corrected Eq. (3.33) is
〈ι̇〉 = S sin(ι)(1 − e
2)3/2
p11/2
1− e2
S2 sin(2ι)
240p6
8 + 3e2
8 + e2
Q cot(ι)(1 − e2)3/2
312 + 736e2 − 83e4 −
408 + 1268e2 + 599e4
cos(2ι)
. (0.3)
http://arxiv.org/abs/0704.0389v8
Evolution of the Carter constant for inspirals into a black hole: effect of the black hole
quadrupole
Éanna É. Flanagan1,2 and Tanja Hinderer1
Center for Radiophysics and Space Research, Cornell University, Ithaca, NY 14853, USA
Laboratory for Elementary Particle Physics, Cornell University, Ithaca, NY 14853, USA
(Dated: November 4, 2018)
We analyze the effect of gravitational radiation reaction on generic orbits around a body with an
axisymmetric mass quadrupole moment Q to linear order in Q, to the leading post-Newtonian order,
and to linear order in the mass ratio. This system admits three constants of the motion in absence
of radiation reaction: energy, angular momentum along the symmetry axis, and a third constant
analogous to the Carter constant. We compute instantaneous and time-averaged rates of change
of these three constants. For a point particle orbiting a black hole, Ryan [15] has computed the
leading order evolution of the orbit’s Carter constant, which is linear in the spin. Our result, when
combined with an interaction quadratic in the spin (the coupling of the black hole’s spin to its own
radiation reaction field), gives the next to leading order evolution. The effect of the quadrupole, like
that of the linear spin term, is to circularize eccentric orbits and to drive the orbital plane towards
antialignment with the symmetry axis.
In addition we consider a system of two point masses where one body has a single mass multipole
or current multipole of order l. To linear order in the mass ratio, to linear order in the multipole,
and to the leading post-Newtonian order, we show that there does not exist an analog of the Carter
constant for such a system (except for the cases of an l = 1 current moment and an l = 2 mass
moment). Thus, the existence of the Carter constant in Kerr depends on interaction effects between
the different multipoles. With mild additional assumptions, this result falsifies the conjecture that
all vacuum, axisymmetric spacetimes posess a third constant of the motion for geodesic motion.
PACS numbers: 04.25.Nx, 04.30.Db
I. INTRODUCTION AND SUMMARY
The inspiral of stellar mass compact objects with
masses µ in the range µ ∼ 1 − 100M⊙ into massive
black holes with masses M ∼ 105 − 107M⊙ is one of
the most important sources for the future space-based
gravitational wave detector LISA. Observing such events
will provide a variety of information: (i) the masses and
spins of black holes can be measured to high accuracy
(∼ 10−4); which can constrain the black hole’s growth
history [1]; (ii) the observations will give a precise test
of general relativity in the strong field regime and unam-
biguously identify whether the central object is a black
hole [2]; and (iii) the measured event rate will give in-
sight into the complex stellar dynamics in galactic nu-
clei [1]. Analogous inspirals may also be interesting for
the advanced stages of ground-based detectors: it has
been estimated that advanced LIGO could detect up to
∼ 10 − 30 inspirals per year of stellar mass compact
objects into intermediate mass black holes with masses
M ∼ 102 − 104M⊙ in globular clusters [3]. Detect-
ing these inspirals and extracting information from the
datastream will require accurate models of the gravita-
tional waveform as templates for matched filtering. For
computing templates, we therefore need a detailed un-
derstanding of the how radiation reaction influences the
evolution of bound orbits around Kerr black holes [4–7].
There are three dimensionless parameters characteriz-
ing inspirals of bodies into black holes:
• the dimensionless spin parameter a = |S|/M2 of
the black hole, where S is the spin.
• the strength of the interaction potential ǫ2 =
GM/rc2, i.e. the expansion parameter used in post-
Newtonian (PN) theory.
• the mass ratio µ/M .
For LISA data analysis we will need waveforms that are
accurate to all orders in a and ǫ2, and to leading order
in µ/M . However, it is useful to have analytic results in
the regimes a ≪ 1 and/or ǫ2 ≪ 1. Such approximate
results can be useful as a check of numerical schemes
that compute more accurate waveforms, for scoping out
LISA’s data analysis requirements [1, 6], and for assessing
the accuracy of the leading order in µ/M or adiabatic
approximation [8–10]. There is substantial literature on
such approximate analytic results, and in this paper we
will extend some of these results to higher order.
A long standing difficulty in computing the evolution
of generic orbits has been the evolution of the orbit’s
”Carter constant”, a constant of motion which governs
the orbital shape and inclination. A theoretical prescrip-
tion now exists for computing Carter constant evolution
to all orders in ǫ and a in the adiabatic limit µ ≪ M
[9, 11–13], but it has not yet been implemented numer-
ically. In this paper we focus on computing analytically
the evolution of the Carter constant in the regime a≪ 1,
ǫ ≪ 1, µ/M ≪ 1, extending earlier results by Ryan
[14, 15].
We next review existing analytical work on the effects
of multipole moments on inspiral waveforms. For non-
spinning point masses, the phase of the l = 2 piece of
the waveform is known to O(ǫ7) beyond leading order
[16], while spin corrections are not known to such high
order. To study the leading order effects of the central
body’s multipole moments on the inspiral waveform, in
the test mass limit µ ≪ M , one has to correct both the
conservative and dissipative pieces of the forces on the
bodies. For the conservative pieces, it suffices to use the
Newtonian action for a binary with an additional multi-
pole interaction potential. For the dissipative pieces, the
multipole corrections to the fluxes at infinity of the con-
served quantities can simply be added to the known PN
point mass results. The lowest order spin-orbit coupling
effects on the gravitational radiation were first derived by
Kidder [17], then extended by Ryan [14, 15], Gergely [18],
and Will [19]. Recently, the corrections of O(ǫ2) beyond
the leading order to the spin-orbit effects on the fluxes
were derived [20, 21]. Corrections to the waveform due
to the quadrupole - mass monopole interaction were first
considered by Poisson [22], who derived the effect on the
time averaged energy flux for circular equatorial orbits.
Gergely [23] extended this work to generic orbits and
computed the radiative instantaneous and time averaged
rates of change of energy E, magnitude of angular mo-
mentum |L|, and the angle κ = cos−1(S ·L) between the
spin S and orbital angular momentum L. Instead of the
Carter constant, Gergely identified the angular average
of the magnitude of the orbital angular momentum, L̄, as
a constant of motion. The fact that to post-2-Newtonian
(2PN) order there is no time averaged secular evolution of
the spin allowed Gergely to obtain expressions for L̇ and
κ̇ from the quadrupole formula for the evolution of the
total angular momentum J = L+S. In a different paper,
Gergely [18] showed that in addition to the quadrupole,
self-interaction spin effects also contribute at 2PN order,
which was seen previously in the black hole perturbation
calculations of Shibata et al. [24]. Gergely calculated
the effect of this interaction on the instantaneous and
time-averaged fluxes of E and |L| but did not derive the
evolution of the third constant of motion.
In this paper, we will re-examine the effects of the
quadrupole moment of the black hole and of the leading
order spin self interaction. For a black hole, our analysis
will thus contain all effects that are quadratic in spin to
the leading order in ǫ2 and in µ/M . Our work will extend
earlier work by
• Considering generic orbits.
• Using a natural generalization of the Carter-type
constant that can be defined for two point particles
when one of them has a quadrupole. This facilitates
applying our analysis to Kerr inspirals.
• Computing instantaneous as well as time-averaged
fluxes for all three constants of motion: energy
E, z-component of angular momentum Lz, and
Carter-type constant K. For most purposes, only
time-averaged fluxes are needed as only they are
gauge invariant and physically relevant. However,
there is one effect for which the time-averaged
fluxes are insufficient, namely transient resonances
that occur during an inspiral in Kerr in the vicin-
ity of geodesics for which the radial and azimuthal
frequencies are commensurate [10, 25]. The instan-
taneous fluxes derived in this paper will be used in
[10] for studying the effect of these resonances on
the gravitational wave phasing.
We will analyze the effect of gravitational radiation re-
action on orbits around a body with an axisymmetric
mass quadrupole moment Q to leading order in Q, to the
leading post-Newtonian order, and to leading order in the
mass ratio. With these approximations the adiabatic ap-
proximation holds: gravitational radiation reaction takes
place over a timescale much longer than the orbital pe-
riod, so the orbit looks geodesic on short timescales. We
follow Ryan’s method of computation [14]: First, we cal-
culate the orbital motion in the absence of radiation re-
action and the associated constants of motion. Next, we
use the leading order radiation reaction accelerations that
act on the particle (given by the Burke-Thorne formula
[26] augmented by the relevant spin corrections [14]) to
compute the evolution of the constants of motion. In the
adiabatic limit, the time-averaged rates of change of the
constants of motion can be used to infer the secular or-
bital evolution. Our results show that a mass quadrupole
has the same qualitative effect on the evolution as spin: it
tends to circularize eccentric orbits and drive the orbital
plane towards antialignment with the symmetry axis of
the quadrupole.
The relevance of our result to point particles inspi-
ralling into black holes is as follows. The vacuum space-
time geometry around any stationary body is completely
characterized by the body’s mass multipole moments
IL = Ia1,a2...al and current multipole moments SL =
Sa1,a2...al [27]. These moments are defined as coefficients
in a power series expansion of the metric in the body’s
local asymptotic rest frame [28]. For nearly Newtonian
sources, they are given by integrals over the source as
IL ≡ Ia1,...al =
ρx<a1 . . . xal>d
3x, (1.1)
SL ≡ Sa1,...al =
ρxpvqǫpq<a1xa2 . . . xal>d
3x.(1.2)
Here ρ is the mass density and vq is the velocity, and ”<
· · · >” means ”symmetrize and remove all traces”. For
axisymmetric situations, the tensor multipole moments
IL (SL) contain only a single independent component,
conventionally denoted by Il (Sl) [27]. For a Kerr black
hole of mass M and spin S, these moments are given by
Il + iSl =M
l+1(ia)l, (1.3)
where a is the dimensionless spin parameter defined by
a = |S|/M2. Note that Sl = 0 for even l and Il = 0 for
odd l.
Consider now inspirals into an axisymmetric body
which has some arbitrary mass and current multipoles
Il and Sl. Then we can consider effects that are linear in
Il and Sl for each l, effects that are quadratic in the mul-
tipoles proportional to IlIl′ , IlSl′ , SlSl′ , effects that are
cubic, etc. For a general body, all these effects can be sep-
arated using their scalings, but for a black hole, Il ∝ al
for even l and Sl ∝ al for odd l [see Eq.(1.3)], so the ef-
fects cannot be separated. For example, a physical effect
that scales as O(a2) could be an effect that is quadratic
in the spin or linear in the quadrupole; an analysis in
Kerr cannot distinguish these two possibilities. For this
reason, it is useful to analyze spacetimes that are more
general than Kerr, characterized by arbitrary Il and Sl,
as we do in this paper. For recent work on computing
exact metrics characterized by sets of moments Il and Sl,
see Refs. [29, 30] and references therein.
The leading order effect of the black hole’s multipoles
on the inspiral is the O(a) effect computed by Ryan [15].
This O(a) effect depends linearly on the spin S1 and is
independent of the higher multipoles Sl and Il since these
all scale as O(a2) or smaller. In this paper we compute
the O(a2) effect on the inspiral, which includes the lead-
ing order linear effect of the black hole’s quadrupole (lin-
ear in I2 ≡ Q) and the leading order spin self-interaction
(quadratic in S1).
We next discuss how these O(a2) effects scale with the
post-Newtonian expansion parameter ǫ. Consider first
the conservative orbital dynamics. Here it is easy to see
that fractional corrections that are linear in I2 scale as
O(a2ǫ4), while those quadratic in S1 scale as O(a
2ǫ6).
Thus, the two types of terms cleanly separate. We com-
pute only the leading order, O(a2ǫ4), term. For the dissi-
pative contributions to the orbital motion, however, the
scalings are different. There are corrections to the radi-
ation reaction acceleration whose fractional magnitudes
are O(a2ǫ4) from both types of effects linear in I2 and
quadratic in S1. The effects quadratic in S1 are due to
the backscattering of the radiation off the piece of space-
time curvature due to the black hole’s spin. This effect
was first pointed out by Shibata et al. [24], who com-
puted the time-averaged energy flux for circular orbits
and small inclination angles based on a PN expansion of
black hole perturbations. Later, Gergely [18] analyzed
this effect on the instantaneous and time-averaged fluxes
of energy and magnitude of orbital angular momentum
within the PN framework.
The organization of this paper is as follows. In Sec.
II, we study the conservative orbital dynamics of two
point particles when one particle is endowed with an ax-
isymmetric quadrupole, in the weak field regime, and to
leading order in the mass ratio. In Sec. III, we com-
pute the radiation reaction accelerations and the instan-
taneous and time-averaged fluxes. In order to have all
the contributions at O(a2ǫ4) for a black hole, we include
in our computations of radiation reaction acceleration
the interaction that is quadratic in the spin S1. The ap-
plication to black holes in Sec. IV briefly discusses the
qualitative predictions of our results and also compares
with previous results.
The methods used in this paper can be applied only
to the black hole spin (as analyzed by Ryan [14]) and
the black hole quadrupole (as analyzed here). We show
in Sec. V that for the higher order mass and current
multipole moments taken individually, an analog of the
Carter constant cannot be defined to the order of our
approximations. We then show that under mild assump-
tions, this non-existence result can be extended to exact
spacetimes, thus falsifying the conjecture that all vac-
uum axisymmetric spacetimes possess a third constant
of geodesic motion.
II. EFFECT OF AN AXISYMMETRIC MASS
QUADRUPOLE ON THE CONSERVATIVE
ORBITAL DYNAMICS
Consider two point particles m1 and m2 interacting in
Newtonian gravity, where m2 ≪ m1 and where the mass
m1 has a quadrupole moment Qij which is axisymmetric:
Qij =
d3xρ(r)
xixj −
r2δij
(2.1)
ninj −
. (2.2)
For a Kerr black hole of mass M and dimensionless spin
parameter a with spin axis along n, the quadrupole scalar
is Q = −M3a2.
The action describing this system, to leading order in
m2/m1, is
µv2 − µΦ(r)
, (2.3)
where v = ṙ is the velocity, the potential is
Φ(r) = −M
xixjQij , (2.4)
µ is the reduced mass and M the total mass of the bi-
nary, and we are using units with G = c = 1. We work to
linear order in Q, to linear order in m2/m1, and to lead-
ing order in M/r. In this regime, the action (2.3) also
describes the conservative effect of the black hole’s mass
quadrupole on bound test particles in Kerr, as discussed
in the introduction. We shall assume that the quadrupole
Qij is constant in time. In reality, the quadrupole will
evolve due to torques that act to change the orientation
of the central body. An estimate based on treating m1 as
a rigid body in the Newtonian field ofm2 gives the scaling
of the timescale for the quadrupole to evolve compared to
the radiation reaction time as (see Appendix I for details)
Tevol
(2.5)
Here, we have denoted the dimensionless spin and
quadrupole of the body by S̄ and Q̄ respectively, and
the last relation applies for a Kerr black hole. Since
µ/M ≪ 1, the first factor in Eq. (2.5) will be large, and
since 1/a ≥ 1 and for the relativistic regime M/r ∼ 1,
the evolution time is long compared to the radiation re-
action time. Therefore we can neglect the evolution of
the quadrupole at leading order.
This system admits three conserved quantities, the en-
µv2 + µΦ(r), (2.6)
the z-component of angular momentum
Lz = ez · (µr× v), (2.7)
and the Carter-type constant
K = µ2(r× v)2 − 2Qµ
(n · r)2
(n · v)2 −
. (2.8)
(See below for a derivation of this expression for K).
A. Conservative orbital dynamics in a
Boyer-Lindquist-like coordinate system
We next specialize to units where M = 1. We also
define the rescaled conserved quantities by Ẽ = E/µ,
L̃z = Lz/µ, K̃ = K/µ
2, and drop the tildes. These spe-
cializations and definitions have the effect of eliminating
all factors of µ and M from the analysis. In spherical
polar coordinates (r, θ, ϕ) the constants of motion E and
Lz become
(ṙ2 + r2θ̇2 + r2 sin2 θϕ̇2)−
(1− 3 cos2 θ), (2.9)
Lz = r
2 sin2 θϕ̇. (2.10)
In these coordinates, the Hamilton-Jacobi equation is not
separable, so a separation constant K cannot readily be
derived. For this reason we switch to a different coordi-
nate system (r̃, θ̃, ϕ) defined by
r cos θ = r̃ cos θ̃
r sin θ = r̃ sin θ̃
. (2.11)
We also define a new time variable t̃ by
cos(2θ̃)
dt̃. (2.12)
The action (2.3) in terms of the new variables to linear
order in Q is
r̃2 sin2 θ̃
sin2 θ̃
. (2.13)
However, a difficulty is that the action (2.13) does not
give the same dynamics as the original action (2.3). The
reason is that for solutions of the equations of motion for
the action (2.3), the variation of the action vanishes for
paths with fixed endpoints for which the time interval ∆t
is fixed. Similarly, for solutions of the equations of motion
for the action (2.13), the variation of the action vanishes
for paths with fixed endpoints for which the time interval
∆t̃ is fixed. The two sets of varied paths are not the
same, since ∆t 6= ∆t̃ in general. Therefore, solutions of
the Euler-Lagrange equations for the action (2.3) do not
correspond to solutions of the Euler-Lagrange equations
for the action (2.13). However, in the special case of zero-
energy motions, the extra terms in the variation of the
action vanish. Thus, a way around this difficulty is to
modify the original action to be
µv2 − µΦ(r) + E
. (2.14)
This action has the same extrema as the action (2.3),
and for motion with physical energy E, the energy com-
puted with this action is zero. Transforming to the new
variables yields, to linear order in Q:
r̃2 sin2 θ̃
sin2 θ̃
+ E − QE
cos(2θ̃)
. (2.15)
The zero-energy motions for this action coincide with the
zero energy motions for the action (2.14). We use this
action (2.15) as the foundation for the remainder of our
analysis in this section.
The z-component of angular momentum in terms of
the new variables (r̃, θ̃, ϕ, t̃) is
Lz = r̃
2 sin2 θ̃
sin2 θ̃
. (2.16)
We now transform to the Hamiltonian:
p2r̃ −
− E − Q
sin2 θ̃
+QE cos(2θ̃)
(2.17)
and solve the Hamiltonian Jacobi equation. Denoting
the separation constant by K we obtain the following
two equations for the r̃ and θ̃ motions:
= 2E +
, (2.18)
= K − L
sin2 θ̃
−QE cos(2θ̃). (2.19)
Note that the equations of motion (2.18) and (2.19) have
the same structure as the equations of motion for Kerr
geodesic motion. Using Eqs. (2.18), (2.19) and (2.16)
together with the inverse of the transformation (2.11)
to linear order in Q, we obtain the expression for K in
spherical polar coordinates:
K = r4(θ̇2 + sin2 θϕ̇2) +Q(ṙ cos θ − rθ̇ sin θ)2 + Q
(ṙ2 + r2 θ̇2 + r2 sin2 θϕ̇2)− 2Q
cos2 θ. (2.20)
This is equivalent to the formula (2.8) quoted earlier.
B. Effects linear in spin on the conservative orbital
dynamics
To include the linear in spin effects, we repeat Ryan’s
analysis [14, 15] (he only gives the final, time averaged
fluxes; we will also give the instantaneous fluxes). We
can simply add these linear in spin terms to our results
because any terms of order O(SQ) will be higher than
the order a2 to which we are working. The correction to
the action (2.3) due to spin-orbit coupling is
Sspin−orbit =
−2µSn
iǫijkxj ẋk
. (2.21)
We will restrict our analysis to the case when the unit
vectors ni corresponding to the axisymmetric quadrupole
Qij and to the spin Si coincide, as they do in Kerr.
Including the spin-orbit term in the action (2.3) results
in the following modified expressions for Lz and K:
Lz = n · (µr× v)−
[r2 − (n · r)2], (2.22)
K = (r× v)2 − 4S
n · (r× v)− 2Q
(n · r)2
(n · v)2 − 1
. (2.23)
In terms of the Boyer-Lindquist like coordinates, the con-
served quantities with the linear in spin terms included
Lz = r̃
2 sin2 θ̃
sin2 θ̃ −Q sin4 θ̃
(2.24)
K = r4(θ̇2 + sin2 θϕ̇2)− 4Sr sin2 θϕ̇
cos2 θ +Q(ṙ cos θ − rθ̇ sin θ)2 + QM
(ṙ2 + r2θ̇2 + r2 sin2 θϕ̇2). (2.25)
The equations of motion are
= 2E+
− 4SLz
, (2.26)
= K −
sin2 θ̃
−QE cos(2θ̃). (2.27)
III. EFFECTS LINEAR IN QUADRUPOLE AND
QUADRATIC IN SPIN ON THE EVOLUTION OF
THE CONSTANTS OF MOTION
A. Evaluation of the radiation reaction force
The relative acceleration of the two bodies can be writ-
ten as
a = −∇Φ(r) + arr, (3.1)
where arr is the radiation-reaction acceleration. Combin-
ing this with Eqs. (2.6), (2.22) and (2.23) for E, Lz and
K gives the following formulae for the time derivatives of
the conserved quantities:
Ė = v · arr, (3.2)
L̇z = n · (r× arr), (3.3)
K̇ = 2(r× v) · (r× arr)−
n · (r× arr)
+2Q(n · v) (n · arr)−Qv · arr. (3.4)
The standard expression for the leading order radiation
reaction acceleration acting on one of the bodies is [31]:
ajrr = −
jk xk +
ǫjpqS
pk xkxq +
ǫjpqS
pk xkvq
ǫpq[jS
xqvk. (3.5)
Here the superscripts in parentheses indicate the number
of time derivatives and square brackets on the indices
denote antisymmetrization.
The multipole moments Ijk(t) and Sjk(t) in Eq. (3.5)
are the total multipole moments of the spacetime, i.e.
approximately those of the black hole plus those due to
the orbital motion. The expression (3.5) is formulated
in asymptotically Cartesian mass centered (ACMC) co-
ordinates of the system, which are displaced from the
coordinates used in Sec. II by an amount [28]
δr(t) = − µ
r(t). (3.6)
This displacement contributes to the radiation reaction
acceleration in the following ways:
1. The black hole multipole moments Il and Sl, which
are time-independent in the coordinates used in
Sec. II, will be displaced by δr and thus will con-
tribute to the (l + 1)th ACMC radiative multipole
[28].
2. The constants of motion are defined in terms of the
black hole centered coordinates used in Sec. II, so
the acceleration arr we need in Eqs. (3.2) – (3.4)
is the relative acceleration. This requires calculat-
ing the acceleration of both the black hole and the
point mass in the ACMC coordinates using (3.5),
and then subtracting to find arr = a
rr − aMrr [14].
To leading order in µ, the only effect of the acceler-
ation of the black hole is via a backreaction of the
radiation field: the lth black hole moments couple
to the (l+1)th radiative moments, thus producing
an additional contribution to the acceleration.
For our calculations at O(S1ǫ
3), O(I2ǫ
4), O(S21ǫ
4), we
can make the following simplifications:
• quadrupole corrections: The fractional corrections
linear in I2 = Q that scale as O(a
2ǫ4) require only
the effect of I2 on the conservative orbital dynamics
as computed in Sec. IIA and the Burke-Thorne for-
mula for the radiation reaction acceleration [given
by the first term in Eq. (3.5)].
• spin-spin corrections: As discussed in the intro-
duction, the fractional corrections quadratic in S1
to the conservative dynamics scale as O(a2ǫ6) and
are subleading order effects which we neglect. At
O(a2ǫ4), the only effect quadratic in S1 is the
backscattering of the radiation off the spacetime
curvature due to the spin. As discussed in item 1.
above, the black hole’s current dipole Si = S1δi3
(taking the z-axis to be the symmetry axis) will
contribute to the radiative current quadrupole an
amount
ij = −
S1xiδj3. (3.7)
The black hole’s current dipole Si will couple to
the gravitomagnetic radiation field due to Sij as
discussed in item 2. above, and contribute to the
relative acceleration as [14]:
aj spinrr =
S1δi3S
ij . (3.8)
For our purposes of computing terms quadratic in
the spin, we substitute S
ij for Sij in Eq. (3.8).
Evaluating these quadratic in spin terms requires
only the Newtonian conservative dynamics, i.e. the
results of Sec. II and Eqs. (3.2) – (3.4) with the
quadrupole set to zero.
• linear in spin corrections: Contributions to these
effects are from Eq. (3.5) with the current
quadrupole replaced by just the spin contribution
(3.7), and from Eq. (3.8) evaluated using only the
orbital current quadrupole.
With these simplifications, we replace the expression
(3.5) for the radiation reaction acceleration with
ajrr = −
jk xk +
ǫjpqS
(6) spin
pk xkxq
ǫjpqS
(5) spin
pk xkvq +
ǫpq[jS
(5) spin
S1δi3
(5) orbit
ij + S
(5) spin
. (3.9)
To justify these approximations, consider the scaling of
the contribution of black hole’s acceleration to the orbital
dynamics. The mass and current multipoles of the black
hole contribute terms to the Hamiltonian that scale with
∆H ∼ Slǫ2l+3 & Ilǫ2l+2. (3.10)
Since the Newtonian energy scales as ǫ2, the fractional
correction to the orbital dynamics scale as
∆H/E ∼ Slǫ2l+1 & Ilǫ2l. (3.11)
To O(ǫ4), the only radiative multipole moments that con-
tribute to the acceleration (3.5) are the mass quadrupole
I2, the mass octupole I3, and the current quadrupole S2
(cf. [17]). Since we are focusing only on the leading or-
der terms quadratic in spin (these can simply be added
to the known 2PN point particle and 1.5PN linear in spin
results), the only terms in Eq. (3.5) relevant for our pur-
poses are those given in Eq. (3.9). The results from a
computation of the fully relativistic metric perturbation
for black hole inspirals [24] show that quadratic in spin
corrections to the l = 2 piece compared to the flat space
Burke-Thorne formula first appear at O(a2ǫ4), which is
consistent with the above arguments.
B. Instantaneous fluxes
We evaluate the radiation reaction force as follows.
The total mass and current quadrupole moment of the
system are
QTij = Qij + µxixj , (3.12)
STij = S
ij + xiǫjkmxkẋm, (3.13)
where from Eq. (2.11)
r̃ sin θ̃
cosϕ, r̃ sin θ̃
sinϕ,
r̃ cos θ̃
. (3.14)
Only the second term in Eq. (3.12) contributes to the
time derivative of the quadrupole. We differentiate five
times by using
cos(2θ̃)
, (3.15)
to the order we are working as discussed above. Af-
ter each differentiation, we eliminate any occurrences of
dϕ/dt̃ using Eq. (2.24), and we eliminate any occurrences
of the second order time derivatives d2r̃/dt̃2 and d2θ̃/dt̃2
in favor of first order time derivatives using (the time
derivatives of) Eqs. (2.26) and (2.27). For computing the
terms linear and quadratic in S1, we set the quadrupole
Q to zero in all the formulae. We insert the resulting ex-
pression into the formula (3.9) for the self-acceleration,
and then into Eqs. (3.2) – (3.4). We eliminate (dr̃/dt̃)2,
(dθ̃/dt̃)2, and (dϕ/dt̃) in favor of E, Lz, and K using
Eqs. (2.24) – (2.27). In the final expressions for the in-
stantaneous fluxes, we keep only terms that are of O(S),
O(Q) and O(S2) and obtain the following results:
15r̃4
− 40K
272KE
196K2 +
r̃2 − 3668
Kr̃ − 352KEr̃2 + 1024
Er̃3 +
E2r̃4
−49K2 − 169KL2z + r̃
+ 2r̃2
+ 47KE +
− 152
r̃3E − 16r̃4E2
−562K2 +
Kr̃ −
r̃2 +
KEr̃2 −
r̃3E − 160r̃4E2
cos(2θ̃)
sin(2θ̃)
439K − 926
r̃ − 1528
θ ˙̃r
−K2 + 22
Kr̃ − 28
r̃2 +
KEr̃2 − 236
r̃3E − 32
r̃4E2
cos(2θ̃)− r̃3 sin(2θ̃)
θ ˙̃r
−49K2 + 6KL2z + 2r̃
63K − 16
L2z −
+ r̃2
112KE − 48
− 1652
r̃3E − 224
r̃4E2
, (3.16)
L̇z =
144LzE
− 24KLz
−50K2 + 240KL2z +
Kr̃ − 7376
L2z r̃ +
r̃2 + 56KEr̃2 − 1824
EL2z r̃
Er̃3 +
E2r̃4
50K2 − 62
Kr̃ − 316
r̃2 − 56KEr̃2 − 624
Er̃3 − 128
E2r̃4
cos(2θ̃)
−104K + 64r̃ + 64Er̃2
sin(2θ̃) ˙̃r
660Er̃2 + 753r̃− 360L2z − 435K +
1601r̃+ 1512r̃2E − 1185K
cos 2θ̃
174QLz
sin(2θ̃) ˙̃r
2S2Lz
Er̃2 + 16r̃ − 9K
, (3.17)
20r̃ + 18r̃2E − 15K
280K2 − 14008
Kr̃ +
r̃2 +
Er̃3 − 2528
KEr̃2 +
E2r̃4
−45K2 + r̃L2z(83 + 80r̃E)− 115KL2z + 14Kr̃(6 + 5r̃E)
15r̃7
cos(2θ̃)
−2175K2 + 2975Kr̃+ 80r̃2 + 3012KEr̃2 − 112Er̃3 − 168E2r̃4
15r̃4
3075K − 20r̃ − 192Er̃2
sin(2θ̃)
θ ˙̃r
7K − 2L2z
−3K + 16
+K cos(2θ̃)
3K − 16
r̃ − 24
sin(2θ̃)
−4K +
θ ˙̃r. (3.18)
C. Alternative set of constants of the motion
A body in a generic bound orbit in Kerr traces an
open ellipse precessing about the hole’s spin axis. For
stable orbits the motion is confined to a toroidal region
whose shape is determined by E, Lz, K. The motion
can equivalently be characterized by the set of constants
inclination angle ι, eccentricity e, and semi-latus rectum
p defined by Hughes [32]. The constants ι, p and e are
defined by cos ι = Lz/
K, and by r̃± = p/(1± e), where
r̃± are the turning points of the radial motion, and r̃
is the Boyer-Lindquist radial coordinate. This param-
eterization has a simple physical interpretation: in the
Newtonian limit of large p, the orbit of the particle is an
ellipse of eccentricity e and semilatus rectum p on a plane
whose inclination angle to the hole’s equatorial plane is
ι. In the relativistic regime p ∼M , this interpretation of
the constants e, p, and ι is no longer valid because the
orbit is not an ellipse and ι is not the angle at which the
object crosses the equatorial plane (see Ryan [14] for a
discussion).
We adopt here analogous definitions of constants of
motion ι, e and p, namely
cos(ι) = Lz/
K, (3.19)
= r̃±. (3.20)
Here K is the conserved quantity (2.23) or (2.25), and r̃±
are the turning points of the radial motion using the r̃
coordinate defined by Eq. (2.11), given by the vanishing
of the right-hand side of Eq. (2.26).
We now rewrite our results in terms of the new con-
stants of the motion e, p and ι. We can use Eq. (2.26)
together with the equations (3.19) and (3.20) to write E,
Lz and K as functions of p, e and ι. To leading order in
Q and S we obtain
K = p
1− 2S cos ι
3 + e2
1 + e2
) 2Q cos2 ι
3 + e2
, (3.21)
E = − (1− e
2S cos ι
1− e2
1− e2
cos2 ι− 1
, (3.22)
p cos ι
1− S cos ι
(3 + e2)−
1 + e2
) Q cos2 ι
3 + e2
. (3.23)
As discussed in the introduction, the effects quadratic in
S on the conservative dynamics scale as O(a2ǫ6) and thus
are not included in this analysis to O(a2ǫ4).
Inserting these relations into the expressions (3.16)–
(3.18) gives, dropping terms of O(QS), O(Q2) and
O(QS2):
Ė = −
15p2r̃7
75p4 − 100p3r̃ + p2r̃2
11− 51e2
+ 32pr̃3
1− e2
)− 6r̃4
1− e2
4S cos ι
15p7/2r̃9
735p6 − 2751p5r̃ + 10p4r̃2(365− 6e2)− 128pr̃5(1− e2)2 − 48r̃6(e2 − 1)3
64S cos ι
15p3/2r̃6
5p(−23 + 3e2)− 3r̃(−9 + e2 + 8e4)
15p4r̃9
4005p6 − 6499p5r̃ + 2p4r̃2
1577− 1977e2
− 24r̃6
1− e2
)3 − 32p3r̃3
8− 33e2
+ 64pr̃5
1− 2e2 + e4
15p4r̃9
24p2r̃4
5− 27e2 + 22e4
− pr̃3 sin(2θ̃)
6585p2 − 4630pr̃+ 2292r̃2(1 − e2)
θ ˙̃r
15p4r̃9
2p2 cos(2θ̃)
4215p4 − 7495p3r̃ + 4p2r̃2(1151− 951e2)− 1012pr̃3(1− e2) + 300r̃4(1− 2e2 + e4)
15p4r̃9
cos(2ι)
2535p6 − 3307p5r̃ + 12p4r̃2(37− 237e2)− 48r̃6(1− e2)3 + 800p3r̃3(1 + e2) + 128pr̃5(1− 2e2 + e4)
15p2r̃5
cos(2ι)
1 + 2e2 − 3e4
15p2r̃9
84r̃4(1− e2)2(1 + e2)2 + 345p4 − 905p3r̃ − 413pr̃3(1− e2) + 2p2r̃2(446− 201e2)
15p2r̃9
cos(2θ̃)
15p4 − 110p3r̃ + 4p2r̃2(47− 12e2)− 118pr̃3(1− e2) + 24r̃4(1− e2)2(1 + e2)2
15r̃9
cos(2ι)
45p2 − 80pr̃ + 36r̃2(1− e2)
15pr̃6
sin(2θ̃) ˙̃r
15p2 + 10pr̃ − 12r̃2(1− e2)
, (3.24)
L̇z = −
8 cos ι
15p2 − 20pr̃ + 9r̃2(1 − e2)
15p2r̃7
525p4 − 1751p3r̃ + 34p2r̃2(61− 6e2) + 12pr̃3(−69 + 29e2) + 6r̃4(17 + 2e2 − 19e4)
15p2r̃7
375p4 − 93p3r̃ + 468pr̃3(1− e2)− 10p2r̃2(58 + 21e2)− 48r̃4(1 − 2e2 + e4)
cos(2θ̃)
15p2r̃7
450p4 − 922p3r̃ − 60pr̃3(3 + e2)− 9p2r̃2(−83 + 23e2) + 27r̃4(1 + 2e2 − 3e4)
cos(2ι)
13p2 − 8pr̃ + 4r̃2(1− e2)
sin(2θ̃) ˙̃r
− Q cos ι
5p5/2r̃7
615p4 − 753p3r̃ + 15p2r̃2
19− 31e2
+ 20pr̃3
1 + 3e2
+ 9r̃4
1− 6e2 + 5e4
− Q cos ι
5p1/2r̃7
cos(2θ̃)
1185p2 − 1601pr̃+ 756r̃2(1− e2)
− 2Q cos ι
5p5/2r̃7
2 cos(2ι)
45p4 − 18r̃4e2(1− e2)− 45p2r̃2(1 + e2) + 20pr̃3(1 + e2)
− 435p3r̃3 sin(2θ̃) ˙̃θ ˙̃r
2 cos ι
p1/2r̃7
9p2 − 16pr̃ + 36
r̃2(1− e2)
, (3.25)
20pr̃ − 15p2 − 9r̃2(1− e2)
8S cos ι
15p3/2r̃7
525p4 − 1751p3r̃ + 2p2r̃2(1172− 57e2) + 12pr̃3(−99 + 19e2)− 24r̃4(−11 + 4e2 + 7e4)
5p2r̃7
−615p4 + 753p3r̃ + 30p2r̃2(17e2 − 9) + 72r̃4e2(1− e2)− 40pr̃3(1 + 3e2)
5p2r̃7
cos(2ι)
−345p4 + 249p3r̃ − 160pr̃3(1 + e2) + 120p2r̃2(1 + 3e2) + 36r̃4(1 + 2e2 − 3e4)
15p2r̃7
2 cos(2θ̃)
2175p4 − 2975p3r̃ − 56pr̃3(1 − e2) + 2p2r̃2(713− 753e2) + 42r̃4(1 − 2e2 + e4)
15pr̃4
sin(2θ̃)
3075p2 − 20pr̃ + 96r̃2(1 − e2)
−9p2 + 16pr̃ − 36
r̃2(1− e2)
cos(2θ̃) + cos(2ι)
3p2 − 16
pr̃ +
r̃2(1 − e2)
sin(2θ̃) ˙̃r
−2p2 + 7
pr̃ − 4
r̃2(1− e2)
. (3.26)
D. Time averaged fluxes
In this section we will compute the infinite time-
averages 〈Ė〉, 〈L̇z〉 and 〈K̇〉 of the fluxes. These averages
are defined by
〈Ė〉 ≡ lim
∫ T/2
Ė(t)dt. (3.27)
These time-averaged fluxes are sufficient to evolve or-
bits in the adiabatic regime (except for the effect of res-
onances) [12, 25]. In Appendix II, we present two dif-
ferent ways of computing the time averages. The first
approach is based on decoupling the r̃ and θ̃ motion us-
ing the analog of the Mino time parameter for geodesic
motion in Kerr [12]. The second approach uses the ex-
plicit Newtonian parameterization of the orbital motion.
Both averaging methods give the following results:
〈Ė〉 = −32
(1− e2)3/2
e4 − S
cos(ι)
cos(2ι)
cos(2ι)
,(3.28)
〈L̇z〉 = −
(1 − e2)3/2
cos ι
e2 − S
2p3/2 cos ι
+ 7e2 +
cos(2ι)
45 + 148e2 +
cos(2ι)
1 + 3e2 +
, (3.29)
〈K̇〉 = −64
(1 − e2)3/2
e2 − S
2p3/2
+ 37e2 +
cos(ι)
cos(2ι)
cos(2ι)
. (3.30)
Using Eqs. (3.21) and (3.23), we obtain from (3.28) – (3.30) the following time averaged rates of change of the
orbital elements e, p, ι:
〈ṗ〉 = −64
(1− e2)3/2
− S cos(ι)
96p3/2
1064 + 1516e2 + 475e4
149e2
469e2
227e4
cos(2ι)
+ e2 +
[13− cos(2ι)]
, (3.31)
〈ė〉 = −
e(1− e2)3/2
121e2
Se(1− e2)3/2 cos(ι)
5p11/2
1172 + 932e2 +
1313e4
Q(1− e2)3/2
785e2
219e4
+ 13e6 +
2195e2
+ 251e4 +
218e6
cos(2ι)
S2e(1− e2)3/2
2 + 3e2 +
[13− cos(2ι)] , (3.32)
〈ι̇〉 = S sin(ι)(1 − e
2)3/2
p11/2
1− e2
S2 sin(2ι)
240p6
8 + 3e2
8 + e2
Q cot(ι)(1 − e2)3/2
312 + 736e2 − 83e4 −
408 + 1268e2 + 599e4
cos(2ι)
. (3.33)
IV. APPLICATION TO BLACK HOLES
A. Qualitative discussion of results
The above results for the fluxes, Eqs. (3.31), (3.32)
and (3.33) show that the correction terms at O(a2ǫ4)
due to the quadrupole have the same type of effect on the
evolution as the linear spin correction computed by Ryan:
they tend to circularize eccentric orbits and change the
angle ι such as to become antialigned with the symmetry
axis of the quadrupole.
The effects of the terms quadratic in spin are quali-
tatively different. In the expression (3.28) for 〈Ė〉, the
coefficient of cos(2ι) due to the spin self-interaction has
the same sign as the quadrupole term, while the terms
not involving ι have the opposite sign. The terms in-
volving cos(2ι) in Eq. (3.30) for 〈K̇〉 of O(Q) and O(S2)
terms have the same sign, while the terms not involving
ι have the opposite sign. The fractional spin-spin cor-
rection to 〈L̇z〉, Eq. (3.29), has no ι-dependence, and in
expression (3.33) for 〈ι̇〉, the dependence on ι of the two
effects O(Q) and O(S2) is different, too. This is not sur-
prising as the O(Q) effects included here are corrections
to the conservative orbital dynamics, while the effects of
O(S2) that we included are due to radiation reaction.
B. Comparison with previous results
The terms linear in the spin in our results for the time
averaged fluxes, Eqs. (3.28) – (3.33), agree with those
computed by Ryan, Eqs. (14a) – (15c) of [15], and with
those given in Eqs. (2.5) – (2.7) of Ref. [33], when we use
the transformations to the variables used by Ryan given
in Eqs. (2.3) – (2.4) in [33].
Equation (3.28) for the time averaged energy flux
agrees with Eq. (3.10) of Gergely [23] and Eq. (4.15)
of [18] when we use the following transformations:
K = L̄2
Ā2 sin2 κ cos δ − (1− Ā2) cos2 κ
= L̄2
E cos2 κ
(1 + 2L̄2) sin2 κ cos δ
, (4.1)
cos ι = cosκ
E cos2 κ
(1 + 2L̄2) sin2 κ cos δ
, (4.2)
(δ + κ), (4.3)
ξ0 = (ψ0 − ψi) +
, (4.4)
where Ā, L̄, κ, δ, ψ0 and ψi are the quantities used by
Gergely. The first relation here is obtained from the turn-
ing points of the radial motion as follows. We compute
r̃± in terms of E and K and map these expressions back
to r using Eqs. (2.11). The result can then be com-
pared with the turning points in Gergely’s variables, Eq.
(2.19) of [23], using the fact that E is the same in both
cases. Instead of the evolution of the constants of motion
K and Lz, Gergely computes the rates of change of the
magnitude L of the orbital angular momentum and of the
angle κ defined by cosκ = (L · S)/L. Using the trans-
formations (4.1) – (4.4) and the definition of κ we verify
that our Eq. (3.29) agrees with the 〈L̇z〉 computed using
Gergely’s Eqs. (3.23) and (3.35) in [23] and Eq. (4.30)
of [18].
In the limit of the circular equatorial orbits analyzed
by Poisson [22], our Eq. (3.28) agrees with Poisson’s Eq.
(22) when we use the transformations and specializations:
, (4.5)
ι = 0, (4.6)
e2 = 0, (4.7)
cosαA = 1, (4.8)
where v and αA are the variables used by Poisson and the
relation (4.5) is obtained by comparing the expressions
for the constants of motion in the two sets of variables.
The main improvement of our analysis over Gergely’s
is that we express the results in terms of the Carter-type
constant K, which facilitates comparing our results with
other analyses of black hole inspirals. Our computations
also include the spin curvature scattering effects for all
three constants of motion; Gergely [18] only considers
these effects for two of them: the energy and magnitude
of angular momentum, not for the third conserved quan-
tity.
When we expand Eq. (3.28) for small inclination an-
gles and specialize to circular orbits, then after converting
p to the parameter v using Eq. (4.5), we obtain
〈Ė〉 = − 32
11Q− S
= − 32
33− 527
. (4.9)
This result agrees with the terms atO(a2v4) of Eq. (3.13)
of Shibata et al. [24], whose calculations were based on
the fully relativistic expressions. This agreement is a
check that we have taken into account all the contribu-
tions at O(a2ǫ4). The analysis in Ref. [24] could not dis-
tinguish between effects due to the quadrupole and those
due curvature scattering, but we can see from Eq. (4.9)
that those two interactions have the opposite dependence
on ι. Comparing (4.9) with Eq. (3.7) of [24] (which gives
the fluxes into the different modes (l = 2,m, n), where m
and n are the multiples of the ϕ and θ frequencies), we see
that the terms in the (2,±2, 0) and the (2,±1,±1) modes
are entirely due to the quadrupole, while the spin-spin in-
teraction effects are fully contained in the (2,±1, 0) and
(2, 0,±1) modes.
V. NON-EXISTENCE OF A CARTER-TYPE
CONSTANT FOR HIGHER MULTIPOLES
In this section, we show that for a single axisymmetric
multipole interaction, it is not possible to find an ana-
log of the Carter constant (a conserved quantity which
does not correspond to a symmetry of the Lagrangian),
except for the cases of spin (treated by Ryan [15]) and
mass quadrupole moment (treated in this paper). Our
proof is valid only in the approximations in which we
work – expanding to linear order in the mass ratio, to
the leading post-Newtonian order, and to linear order in
the multipole. However we will show below that with
very mild additional smoothness assumptions, our non-
existence result extends to exact geodesic motion in exact
vacuum spacetimes.
We start in Sec. VA by showing that there is no co-
ordinate system in which the Hamilton-Jacobi equation
is separable. Now separability of the Hamilton-Jacobi
equation is a sufficient but not a necessary condition for
the existence of a additional conserved quantity. Hence,
this result does not yield information about the existence
or non-existence of an additional constant. Nevertheless
we find it to be a suggestive result. Our actual derivation
of the non-existence is based on Poisson bracket compu-
tations, and is given in Sec. VB.
A. Separability analysis
Consider a binary of two point masses m1 and m2,
where the mass m1 is endowed with a single axisymmet-
ric current multipole moment Sl or axisymmetric mass
multipole moment Il. In this section, we show that the
Hamilton-Jacobi equation for this motion, to linear order
in the multipoles, to linear order in the mass ratio and to
the leading post-Newtonian order, is separable only for
the cases S1 and I2.
We choose the symmetry axis to be the z-axis and write
the action for a general multipole as
ṙ2 + r2θ̇2 + r2 sin2 θϕ̇2
+ f(r, θ) + g(r, θ)ϕ̇+ E] . (5.1)
For mass moments, g(r, θ) = 0, while for current mo-
ments f(r, θ) = 0. For an axisymmetric multipole of
order l, the functions f and g will be of the form
f(r, θ) =
clIlPl(cos θ)
, g(r, θ) =
dlSl sin θ∂θPl(cos θ)
(5.2)
where Pl(cos θ) are the Legendre polynomials and cl and
dl are constants. We will work to linear order in f and g.
In Eq. (5.1), we have added the energy term needed when
doing a change of time variables, cf. the discussion before
Eq. (2.14) in section III. Since ϕ is a cyclic coordinate,
pϕ = Lz is a constant of motion and the system has
effectively only two degrees of freedom. Note that in the
case of a current moment, there will be correction term
in Lz:
Lz = r
2 sin2 θϕ̇+ g(r, θ). (5.3)
Next, we switch to a different coordinate system
(r̃, θ̃, ϕ) defined by
r = r̃ + α(r̃, θ̃, Lz), (5.4)
θ = θ̃ + β(r̃, θ̃, Lz), (5.5)
where the functions α and β are yet undetermined. We
also define a new time variable t̃ by
1 + γ(r̃, θ̃, Lz)
dt̃. (5.6)
Since we work to linear order in f and g, we can work
to linear order in α, β, and γ. We then compute the
action in the new coordinates and drop the tildes. The
Hamiltonian is given by
p2r(1 + γ − 2α,r) +
(1− 2α
− 2β,θ + γ)
(−α,θ − r2β,r)− E(1 + γ)
2r2 sin2 θ
(1 + γ − 2α
− 2β cot θ)
(1− α
+ γ)− f − gLz
r2 sin2 θ
(5.7)
and the corresponding Hamilton-Jacobi equation is
Ĉ1 +
+ 2V̂ , (5.8)
where we have denoted
Ĉ1 = J(r, θ) [1 + γ − 2α,r] = 1 + γ − 2α,r + j, (5.9)
Ĉ2 = J(r, θ)
1− 2α
− 2β,θ + γ
= 1− 2α
− 2β,θ + γ + j, (5.10)
Ĉ3 = J(r, θ)
−α,θ − r2β,r
= −α,θ − r2β,r, (5.11)
V̂ = J(r, θ)
2r2 sin2 θ
(1 + γ − 2α
− 2β cot θ)
+ γ)− E(1 + γ)
− f − gLz
r2 sin2 θ
2r2 sin2 θ
(1 + γ − 2α
− 2β cot θ + j)
−E(1 + γ + j)− 1
(1− α
+ γ + j)
−f − gLz
r2 sin2 θ
. (5.12)
The unperturbed problem is separable, so make the
perturbed problem separable, we have multiplied the
Hamilton-Jacobi equation by an arbitrary function
J(r, θ), which can be expanded as J(r, θ) = 1 + j(r, θ),
where j(r, θ) is a small perturbation.
To find a solution of the form W =Wr(r)+Wθ(θ), we
first specialize to the case where Ĉ3 = 0:
− Ĉ3 = β,rr2 + α,θ = 0. (5.13)
We differentiate Eq. (5.8) with respect to θ, using Eq.
(5.8) to write (dWr/dr)
2 in terms of (dWθ/dθ)
2 and then
differentiate the result with respect to r to obtain
∂θĈ2
∂θĈ1
2V̂ ∂θĈ1
Ĉ1Ĉ2
. (5.14)
Expanding Eq. (5.14) to linear order in the small quan-
tities then yields the two conditions for the kinetic and
the potential part of the Hamiltonian to be separable:
0 = ∂r∂θ
2α,r −
− 2β,θ
, (5.15)
sin2 θ
2β,r cot
2 θ − 3β,rθ cot θ + β,r csc2 θ
sin2 θ
+ α,rθ
−∂r∂θ
Pl(cos θ) +
dlSlLz
rl sin θ
∂θPl(cos θ)
2α,rθ −
+ 2Er2α,rθ
, (5.16)
where we have used Eq. (5.2) for f and g. Therefore, the
following conditions must be satisfied:
M4(θ)−N(r) =
+ β,θ − 2α,r, (5.17)
M1(θ) = 2β cot
2 θ + β csc2 θ + β,θθ
−3β,θ cot θ, (5.18)
M2(θ) = r
2∂r(r
2β,r), (5.19)
M3(θ) = 2rα,rθ − α,θ +
∂θPl(cos θ)
−SlLz
∂θ(csc θ ∂θPl(cos θ)). (5.20)
Here, the functions M and N are arbitrary integration
constants.
Solving the condition for the kinetic term to be sep-
arable, Eq. (5.17), together with Eq. (5.13) gives the
general solution that goes to zero at large r as
cos(nθ + ν), (5.21)
β = − A
sin(nθ + ν), (5.22)
where A and ν are arbitrary and n is an integer. These
functions must satisfy the conditions (5.18) – (5.20) in
order for the potential term to be separable as well. To
see when this will be the case, we start by considering Eq.
(5.20). Substituting the general ansatz α = a1(r)a2(θ)
shows that a′2 = P
l or a
2 = (cscθ P
′ depending on
whether a mass or a current multipole is present. The
function a1(r) is then determined from
0 = 2ra′1 − a1 +
clIl/r
(l−1)
dlSlLz/r
(5.23)
Hence,
[clIl/(2l)] r
(1−l)
[dlSlLz/(2l+ 1)] r
(5.24)
so that we obtain for mass moments
Pl(cos θ)
, β =
P ′l (cos θ)
(5.25)
and for current moments
dlSlLz
2l+ 1
csc θP ′l (cos θ)
, (5.26)
dlSlLz
(2l+ 1)(l + 1)
(csc θ P ′l (cos θ))
, (5.27)
where we have used the condition (5.13) to solve for β.
Substituting this in Eq. (5.19) determines that l = 2
for mass moments and l+1 = 2 for current moments. For
an l = 2 mass moment, conditions (5.17) and (5.18) are
satisfied as well, with n = 2 and ν = 0. For the case of an
l = 1 current moment, the extra term inH is independent
of θ anyway. But for any other multipole interaction,
the Hamilton-Jacobi equation will not be separable. For
example, for the current octupole Sijk, the last term in
Eq. (5.7) is proportional to S3Lz(5 cos
2 θ − 1)/r5 and
is therefore not separable. From Eq. (5.2) one can see
that, for a general multipole, the functions f or g contain
different powers of cos θ appearing with the same power
of r since the Legendre polynomials can be expanded as
[34]:
Pl(cos θ) =
(−1)n(2l− 2n)!
2ln!(l − n)!(l − 2n)!
(cos θ)l−2n, (5.28)
where N = l/2 for even l and N = (l + 1)/2 for odd l.
It will not be possible to cancel all of these terms with
(5.21) – (5.22) for l > 2.
The case when Ĉ3 is non-vanishing will only be sepa-
rable if all the coefficients are functions of r or of θ only,
and if in addition, the potential also depends only on r or
on θ. Achieving this for our problem will not be possible
because the potential cannot be transformed to the form
required for separability.
B. Derivation of non-existence of additional
constants of the motion
In this subsection, we show using Poisson brackets that
for a single axisymmetric multipole interaction, to linear
order in the multipole and the mass ratio, a first integral
analogous to the Carter constant does not exist, except
for the cases of mass quadrupole and spin.
Suppose that such a constant does exist. We write the
Hamiltonian corresponding to the action (5.1) as H =
H0 + δH and the Carter-type constant as K = K0 +
δK(pr, pθ, Lz, r, θ), where
2r2 sin2 θ
, (5.29)
δH = − clIl
Pl(cos θ)−
dlSlLz
rl+2 sin θ
∂θPl(cos θ),(5.30)
K0 = p
sin2 θ
. (5.31)
Computing the Poisson bracket gives, to linear order in
the perturbations
0 = {H0, δK}+ {δH,K0} (5.32a)
δK + {δH,K0}, (5.32b)
where we have used that {H0,K0} = 0 and the fact that
{H0, δK} = d(δK)/dt. Here, d/dt denotes the total time
derivative along an orbit (r(t), θ(t), pr(t), pθ(t)) of H0 in
phase space. The partial differential equation (5.32a) for
δK thus reduces to a set of ordinary differential equa-
tions that can be integrated along the individual orbits
in phase space.
The unperturbed motion for a bound orbit is in a
plane, so we can switch from spherical to plane polar co-
ordinates (r, ψ). In terms of these coordinates, we have
H0 = p
r/2+p
ψ/2, K0 = p
ψ, and cos θ = sin ι sin(ψ+ψ0),
with cos ι = Lz/
K and the constant ψ0 denoting the
angle between the direction of the periastron and the
intersection between the orbital and equatorial plane.
Then Eq. (5.32) becomes
δK = η(t), (5.33)
η(t) = − 2pψ dlSlLz
sin ι rl+2(t)
∂ψPl(sin ι sin(ψ(t) + ψ0))
cos(ψ(t) + ψ0)
2pψ clIl
rl+1(t)
∂ψPl(sin ι sin(ψ(t) + ψ0)). (5.34)
For unbound orbits, one can always integrate Eq.
(5.33) to determine δK. However, for bound periodic
orbits there is a possible obstruction: the solution for
the conserved quantity K0 + δK will be single valued if
and only if the integral of the source over the closed orbit
vanishes,
∮ Torb
η(t)dt = 0. (5.35)
Here, Torb is the orbital period. In other words, the par-
tial differential equation (5.32) has a solution δK if and
only if the condition (5.35) is satisfied. This is the same
condition as obtained by the Poincare-Mel’nikov-Arnold
method, a technique for showing the non-integrability
and existence of chaos in certain classes of perturbed dy-
namical systems [35].
Thus, it suffices to show that the condition (5.35) is
violated for all multipoles other than the spin and mass
quadrupole. To perform the integral in Eq. (5.35), we use
the parameterization for the unperturbed motion, r =
K/(1+ e cosψ) and dt/dψ = K3/2/(1+ e cosψ)2, so that
the condition for the existence of a conserved quantity
K0 + δK becomes
clIl(1 + e cosψ)
l−1∂ψPl(sin ι sin(ψ + ψ0))−
dlSlLz
K sin ι
(1 + e cosψ)l∂ψ
∂ψPl(sin ι sin(ψ + ψ0))
cos(ψ + ψ0)
(5.36)
In terms of the variable χ = ψ + ψ0 − π/2, Eq. (5.36) can be written as
dχclIl [1 + e(sinψ0 cosχ− cosψ0 sinχ)]l−1
Pl(sin ι cosχ)
dlSlLz
sin ι
[1 + e(sinψ0 cosχ− cosψ0 sinχ)]l
Pl(sin ι cosχ)
. (5.37)
Inserting the expansion (5.28) for Pl(cosχ), taking the derivatives, and using the binomial expansion for the first term
in Eq. (5.37), we get
0 = clIl
Alnjk e
j(sin ι)l−2n(sinψ0)
k(cosψ0)
dχ (sinχ)j−k+1(cosχ)k+l−2n−1
dlSlLz
Blnjk e
j(sin ι)l−2n−1(sinψ0)
k(cosψ0)
dχ (sinχ)j−k+1(cosχ)k+l−2n−2. (5.38)
The coefficients Alnkj and Blnkj are
Alnkj =
(−1)n+k+1(l − 1)!(2l − 2n)!
2ln!(l − 1− j)!k!(j − k)!(l − n)!(l − 2n− 1)!
, Blnkj =
(−1)n+kl!(2l− 2n)!
2ln!(l − j)!k!(j − k)!(l − n)!(l − 2n− 2)!
. (5.39)
The only non-vanishing contribution to the integrals in Eq. (5.38) will come from terms with even powers of both
cosχ and sinχ. These can be evaluated as multiples of the beta function:
0 = clIl
Clnjk e
j(sin ι)l−2n(sinψ0)
k(cosψ0)
j−k δ(j−k+1),even δ(l+k−1),even
dlSlLz
Dlnjk e
j(sin ι)l−2n−1(sinψ0)
k(cosψ0)
j−k δ(j−k+1),even δ(l+k),even. (5.40)
Here, the coefficients are
Clnjk =
2Γ( j
+ 1)Γ(k
− n+ 1)
Alnkj , Dlnjk =
2Γ( j
+ 1)Γ(k
− n− 1
− n+ 3
Blnkj (5.41)
Eq. (5.40) shows that for even l, terms with j =even
(odd) and k =odd (even) give a non-vanishing contribu-
tion for the case of a mass (current) multipole, and hence
K0+δK is not a conserved quantity for the perturbed mo-
tion. Note that terms with j =even and k =odd for even
l occur only for l > 3, so for l = 2 the mass quadrupole
term in Eq. (5.40) vanishes and therefore there exists an
analog of the Carter constant, which is consistent with
our results of Sec. II and our separability analysis. For
odd l, terms with j =odd (even) and k =even (odd) are
finite for Il (Sl). Note that for the case l = 1 of the spin,
the derivatives with respect to χ in Eq. (5.37) evaluate to
zero, so in this case there also exists a Carter-type con-
stant. These results show that for a general multipole
other than I2 and S1, there will not be a Carter-type
constant for such a system.
1. Exact vacuum spacetimes
Our result on the non-existence of a Carter-type con-
stant can be extended, with mild smoothness assump-
tions, to falsify the conjecture that all exact, axisymmet-
ric vacuum spacetimes posess a third constant of the mo-
tion for geodesic motion. Specifically, we fix a multipole
order l, and we assume:
• There exists a one parameter family
(M, gab(λ))
of spacetimes, which is smooth in the parameter λ,
such that λ = 0 is Schwarzschild, and each space-
time gab(λ) is stationary and axisymmetric with
commuting Killing fields ∂/∂t and ∂/∂φ, and such
that all the mass and current multipole moments of
the spacetime vanish except for the one of order l.
On physical grounds, one expects a one parameter
family of metrics with these properties to exist.
• We denote by H(λ) the Hamiltonian on the tan-
gent bundle overM for geodesic motion in the met-
ric gab(λ). By hypothesis, there exists for each λ
a conserved quantity M(λ) which is functionally
independent of the conserved energy and angular
momentum. Our second assumption is that M(λ)
is differentiable in λ at λ = 0. One would expect
this to be true on physical grounds.
• We assume that the conserved quantity M(λ) is
invariant under the symmetries of the system:
L~ξM(λ) = L~ηM(λ) = 0,
where ~ξ and ~η are the natural extensions to the 8
dimensional phase space of the Killing vectors ∂/∂t
and ∂/∂φ. This is a very natural assumption.
These assumptions, when combined with our result of
the previous section, lead to a contradiction, showing
that the conjecture is false under our assumptions.
To prove this, we start by noting that M(0) is a con-
served quantity for geodesic motion in Schwarzschild, so
it must be possible to express it as some function f of
the three independent conserved quantities:
M(0) = f(E,Lz,K0). (5.42)
Here E is the energy, Lz is the angular momentum, and
K0 is the Carter constant. Differentiating the exact re-
lation {H(λ),M(λ)} = 0 and evaluating at λ = 0 gives
{H0,M1} =
{E,H1}+
{Lz, H1}+
{K0, H1},
(5.43)
where H0 = H(0), H1 = H
′(0), and M1 = M
′(0). As
before, we can regard this is a partial differential equa-
tion that determines M1, and a necessary condition for
solutions to exist and be single valued is that the integral
of the right hand side over any closed orbit must vanish:
{E,H1}+
{Lz, H1}+
{K0, H1}
(5.44)
Now strictly speaking, there are no closed orbits in
the eight dimensional phase space. However, the ar-
gument of the previous section applies to orbits which
are closed in the four dimensional space with coordinates
(r, θ, pr, pθ), since by the third assumption above every-
thing is independent of t and φ, and pt and pφ are con-
served. Here (t, r, θ, φ) are Schwarzschild coordinates and
(pt, pr, pθ, pφ) are the corresponding conjugate momenta.
Next, we can pull the partial derivatives ∂f/∂E etc.
outside of the integral. It is then easy to see that the first
two terms vanish, since there do exist a conserved energy
and a conserved z-component of angular momentum for
the perturbed system. Thus, Eq. (5.44) reduces to
{K0, H1} = 0. (5.45)
Since M(0) is functionally independent of E and Lz, the
prefactor ∂f/∂K0 must be nonzero, so we obtain
{K0, H1} = 0. (5.46)
The result (5.46) applies to fully relativistic orbits in
Schwarzschild. We need to take the Newtonian limit of
this result in order to use the result we derived in the
previous section. However, the Newtonian limit is a lit-
tle subtle since Newtonian orbits are closed and generic
relativistic orbits are not closed. We now discuss how the
limit is taken.
The integral (5.46) is taken over any closed orbit in
the four dimensional phase space (r, θ, pr, pθ) which cor-
responds to a geodesic in Schwarzschild. Such orbits are
non generic; they are the orbits for which the ratio be-
tween the radial and angular frequencies ωr and ωθ is a
rational number. We denote by qr and qθ the angle vari-
ables corresponding to the r and θ motions [36]. These
variables evolve with proper time τ according to
qr = qr,0 + ωrτ, (5.47a)
qθ = qθ,0 + ωθτ, (5.47b)
where qr,0 and qθ,0 are the initial values. We denote the
integrand in Eq. (5.46) by
I(qr, qθ, a, ε, ι),
where I is some function, and a, ε and ι are the parame-
ters of the geodesic defined by Hughes [32] (functions of
E, Lz and K0). The result (5.46) can be written as
∫ T/2
dτ I[qr(τ), qθ(τ), a, ε, ι] = 0, (5.48)
where T = T (a, ε, ι) is the period of the r, θ motion.
Since the variables qr and qθ are periodic with period
2π, we can express the function I as a Fourier series
I(qr, qθ, a, ε, ι) =
n,m=−∞
Inm(a, ε, ι)einqr+imqθ . (5.49)
Now combining Eqs. (5.47), (5.48) and (5.49) gives
n,m=−∞
Inm(a, ε, ι)einqr,0+imqθ,0
×Si [(nωr +mωθ)T/2] , (5.50)
where Si(x) = sin(x)/x. Since the initial conditions qr,0
and qθ,0 are arbitrary, it follows that
Inm(a, ε, ι)Si [(nωr +mωθ)T/2] = 0 (5.51)
for all n, m.
Next, for closed orbits the ratio of the frequencies must
be a rational number, so
, (5.52)
where p and q are integers with no factor in common.
These integers depend on a, ε and ι. The period T is
given by 2π/T = qωr = pωθ. The second factor in Eq.
(5.51) now simplifies to
(np+mq)π
, (5.53)
which vanishes if and only if
n = n̄q, m = m̄p, n̄+ m̄ 6= 0, (5.54)
for integers n̄, m̄. It follows that
Inm(a, ε, ι) = 0 (5.55)
for all n, m except for values of n, m which satisfy the
condition (5.54)
Consider now the Newtonian limit, which is the limit
a → ∞ while keeping fixed ε and ι and the mass of the
black hole. We denote by IN(qr , qθ, a, ε, ι) the Newtonian
limit of the function I(qr , qθ, a, ε, ι). The integral (5.48)
in the Newtonian limit is given by the above computation
with p = q = 1, since ωr = ωθ in this limit. This gives
dτIN =
INn,−n(a, ε, ι) ein(qr,0−qθ,0), (5.56)
where INnm are the Fourier components of IN. In the
previous subsection, we showed that this function is non-
zero, which implies that there exists a value k of n for
which IN k,−k 6= 0.
Now as a → ∞, we have ωr/ωθ → 1, and hence from
Eq. (5.52) there exists a critical value ac of a such that
the values of p and q exceed k for all closed orbits with
a > ac. (We are keeping fixed the values of ε and ι). It
follows from Eqs. (5.54) and (5.55) that
Ik,−k(a, ε, ι)
IN k,−k(a, ε, ι)
= 0 (5.57)
for all such values of a. However this contradicts the fact
Ik,−k(a, ε, ι)
IN k,−k(a, ε, ι)
→ 1 (5.58)
as a→ ∞. This completes the proof.
Hence, if the three assumptions listed at the start of
this subsection are satisfied, then the conjecture that all
vacuum, axisymmetric spacetimes possess a third con-
stant of the motion is false.
Finally, it is sometimes claimed in the classical dynam-
ics literature that perturbation theory is not a sufficiently
powerful tool to assess whether the integrability of a sys-
tem is preserved under deformations. An example that
is often quoted is the Toda lattice Hamiltonian [38, 39].
This system is integrable and admits a full set of con-
stants of motion in involution. However, if one approx-
imates the Hamiltonian by Taylor expanding the poten-
tial about the origin to third order, one obtains a sys-
tem which is not integrable. This would seem to indicate
that perturbation theory can indicate a non-integrability,
while the exact system is still integrable.
In fact, the Toda lattice example does not invalidate
the method of proof we use here. If we write the Toda
lattice Hamiltonian as H(q,p), then the situation is that
H(λq,p) is integrable for λ = 1, but it is not integrable
for 0 < λ < 1. Expanding H(λq,p) to third order in λ
gives a non-integrable Hamiltonian. Thus, the perturba-
tive result is not in disagreement with the exact result
for 0 < λ < 1, it only disagrees with the exact result for
λ = 1. In other words, the example shows that pertur-
bation theory can fail to yield the correct result for finite
values of λ, but there is no indication that it fails in ar-
bitrarily small neighborhoods of λ = 0. Our application
is qualitatively different from the Toda lattice example
since we have a one parameter family of Hamiltonians
H(λ) which by assumption are integrable for all values
of λ.
VI. CONCLUSION
We have examined the effect of an axisymmetric
quadrupole moment Q of a central body on test parti-
cle inspirals, to linear order in Q, to the leading post-
Newtonian order, and to linear order in the mass ratio.
Our analysis shows that a natural generalization of the
Carter constant can be defined for the quadrupole inter-
action. We have also analyzed the leading order spin self-
interaction effect due to the scattering of the radiation off
the spacetime curvature due to the spin. Combining the
effects of the quadrupole and the leading order effects
linear and quadratic in the spin, we have obtained ex-
pressions for the instantaneous as well as time-averaged
evolution of the constants of motion for generic orbits un-
der gravitational radiation reaction, complete at O(a2ǫ4).
We have also shown that for a single multipole interaction
other than Q or spin, in our approximations, a Carter-
type constant does not exist. With mild additional as-
sumptions, this result can be extended to exact space-
times and falsifies the conjecture that all axisymmetric
vacuum spacetimes possess a third constant of motion for
geodesic motion.
VII. ACKNOWLEDGMENTS
This research was partially supported by NSF grant
PHY-0457200. We thank Jeandrew Brink for useful cor-
respondence.
Appendix A: Time variation of quadrupole: order of
magnitude estimates
In this appendix, we give an estimate of the timescale
Tevol for the quadrupole to change. The analysis in the
body of this paper is valid only when Tevol ≫ Trr, where
Trr is the radiation reaction time, since we have neglected
the time evolution of the quadrupole. We distinguish be-
tween two cases: (i) when the central body is exactly non-
spinning but has a quadrupole, and (ii) when the central
body has finite spin in addition to the quadrupole.
1. Estimate of the scaling for the nonspinning case
For the purpose of a crude estimate, the relevant in-
teraction is the tidal interaction with energy
QijEij ∼ −
Q̄I cos2 θ, (A1)
where Eij is the tidal field, θ is the angle between the
symmetry axis and the normal to the orbital plane of
m2, and we have written the quadrupole as Q ∼ Q̄I,
where Q̄ is dimensionless and I is the moment of inertia.
For small deviations from equilibrium, the relevant piece
of the Lagrangian is schematically
L ∼ Iψ̇2 + Q̄I m2
ψ2. (A2)
We define the evolution timescale Tevol to be the time
it takes for the angle to change by an amount of order
unity, and since the amplitude of the oscillation scales
roughly as ∼ m2/m1, the evolution time scales as
T−2evol ∼
ω2orbit, (A3)
where ω2orbit = M/r
3. Thus, the ratio of the evolution
timescale compared to the radiation reaction timescale
scales as
Tevol/Trr ∼
. (A4)
2. Estimate of the scaling for the spinning case
When the body is spinning the effect of the tidal cou-
pling is to cause a precession. For the purpose of this
estimate, we calculate the torque on m1 due to the com-
panion’s Newtonian field. The torque N scales as
Ni ∼ ǫimjQmkEjk. (A5)
We assume that the precession is slow, i.e.
ωprec ≪ S̄/m1
, (A6)
where ωprec is the precession frequency and S̄ = S/m
is the dimensionless spin. This gives the approximate
scaling of the precession timescale as (cf. [37])
Tprec/Trr ∼
. (A7)
and the evolution timescale is thus
Tevol/Trr ∼
. (A8)
Because of our assumption (A6) that the precession is
slow, equation (A8) is valid only when
) S̄2
. (A9)
When S̄ is sufficiently small that the condition (A9) is
violated, the relevant timescale is instead given by Eq.
(A3).
3. Application to Kerr inspirals
For Kerr inspirals,
S̄ ∼ a, Q̄ ∼ a2, µ/M ≪ 1 and r ∼M. (A10)
Therefore, the condition (A9) is satisfied, and the pre-
cession time is longer than the radiation reaction time
Tprec/Trr ∼
. (A11)
Note that for Kerr inspirals, since r ∼ M both formulas
(A3) and (A7) give the same scaling.
Moreover, for Kerr inspirals, the amplitude of the pre-
cession will be small, of order the mass ratio µ/M . This is
because of angular momentum conservation: in the rela-
tivistic regime, the orbital angular momentum is a factor
of µ/M smaller than the angular momentum of the black
hole and can therefore not cause a large precession ampli-
tude. Even if the orbital angular momentum at infinity
is large, most of it will be radiated away as outgoing
gravitational waves during the earlier phase of the inspi-
ral. This factor of µ/M is taken into account when we
consider the evolution timescale, which for Kerr inspirals
reduces to
Tevol/Trr ∼
. (A12)
Since 1/a ≥ 1, M/r ∼ 1 and M/µ ≪ 1, the evolution
time is long compared to the radiation reaction time and
we can neglect the time variation of the quadrupole at
leading order.
Appendix B: Computation of time averaged fluxes
1. Averaging method that parallels fully
relativistic averaging
We start by noting that the differential equations
(2.26) and (2.27) governing the r̃ and θ̃ motions decouple
if we define a new time parameter t̂ by
dt̂ =
dt̃. (B1)
This is the analog of the Mino time parameter for
geodesic motion in Kerr [12]. The equations of motion
(2.26)–(2.24) then become
= V̂r̃(r̃), (B2)
V̂r̃(r̃) = 2Er̃
4 + 2r̃3 −Kr̃2 − 4SLzr̃
r̃ − 2L2z
, (B3)
= V̂θ̃(θ̃), (B4)
V̂θ̃(θ̃) = K −
sin2 θ̃
−QE cos 2θ̃, (B5)
= V̂ϕr̃(r̃) + V̂ϕθ̃(θ̃), (B6)
V̂ϕr̃(r̃) =
, V̂ϕθ̃(θ̃) =
sin2 θ̃
. (B7)
The parameters t and t̂ are related by:
= V̂tr̃(r̃) + V̂tθ̃(θ̃) (B8)
V̂tr̃(r̃) = r̃
2, V̂tθ̃(θ̃) =
cos 2θ̃. (B9)
It follows from Eqs. (B2) and (B4) that the functions
r̃(t̂) and θ̃(t̂) are periodic; and we denote their periods
by Λr̃ and Λθ̃. We define the fiducial motion associated
with the constants of motion E, Lz and K to be the
motion with the initial conditions r̃(0) = r̃min and θ̃(0) =
θ̃min, where r̃min and θ̃min are given by the vanishing of
the right-hand sides of Eqs. (B2) and (B4) respectively.
The functions r̂(t̂) and θ̂(t̂) associated with this fiducial
motion are given by
∫ r̂(t̂)
r̃min
V̂r̃(r̃)
= t̂, (B10)
∫ θ̂(t̂)
θ̃min
V̂θ̃(θ̃)
= t̂. (B11)
From Eq. (B8) it follows that
t(t̂) = t0 +
V̂tr̃[r̃(t
′)] + V̂tθ̃[θ̃(t
, (B12)
where t0 = t(0). Next, we define the constant Γ to be
the following average value:
∫ Λr̃
dt′V̂tr̃[r̂(t
′)] +
dt′V̂tθ̃[θ̂(t
′)]. (B13)
Then we can write t(t̂) as a sum of a linear term and
terms that are periodic:
t(t) = t0 + Γt̂+ δt(t̂), (B14)
where δt(t̂) denotes the oscillatory terms in Eq. (B12).
To average a function over the time parameter t̂, it is
convenient to parameterize r̃ and θ̃ in terms of angular
variables as follows. For the average over θ̃ we introduce
the parameter χ by
cos2 θ̂(t̂) = z− cos
2 χ, (B15)
where z− = cos
2 θ̃− with z− being the smaller root of Eq.
(B4):
K + 3QE ±
(K −QE)2 + 4QEL2z
(B16)
and where β = 2QE. Then from the definition (B11)
of θ̂ together with Eq. (B4) and the requirement that χ
increases monotonically with t̂ we obtain
β (z+ − z− cos2 χ). (B17)
Then we can write the average over t̂ of a function Fθ̃(t̂)
which is periodic with period Λθ̃ in terms of χ as
〈Fθ̃〉t̂ =
dt̂Fθ̃(t̂)
Fθ̃[t̂(χ)]
β (z+ − z− cos2 χ)
, (B18)
where
Λθ̃ =
β (z+ − z− cos2 χ)
. (B19)
Similarly, to average a function Fr̃(t̂) that is periodic with
period Λr̃, we introduce a parameter ξ via
1 + e cos ξ
, (B20)
where the parameter ξ varies from 0 to 2π as r̃ goes
through a complete cycle. Then,
= P (ξ), (B21)
P (ξ) ≡
V̂r̃[r̃(ξ)]
pe | sin ξ |
(1 + e cos ξ)
(B22)
The average over t̂ of Fr̃(t̂) can then be computed from
〈Fr̃〉t̂ =
dξ Fr̃/P (ξ)
dξ/P (ξ)
. (B23)
Now, a generic function Fr̃,θ̃[r̃(t̂), θ̃(t̂)] will be biperiodic
in t̂: Fr̃,θ̃[r̃(t̂+Λr̃), θ̃(t̂+Λθ̃)] = Fr̃,θ̃[r̃(t̂), θ̃(t̂)]. Combin-
ing the results (B18) and (B23) we can write its average
as a double integral over χ and ξ as
〈Fr̃,θ̃〉t̂ =
Λθ̃Λr̃
Fr̃,θ̃[r̃(ξ), θ̃(χ)]
β (z+ − z− cos2 χ)P (ξ)
(B24)
To compute the time average of Ė, L̇z, and K̇, we need
to convert the average of a function over t̂ calculated from
(B24) to the average over t. As explained in detail in
[9], in the adiabatic limit we can choose a time interval
∆t which is long compared to the orbital timescale but
short compared to the radiation reaction time. From
Eq. (B12) we have ∆t = Γt̂+ osc.terms. The oscillatory
terms will be bounded and will therefore be negligible in
the adiabatic limit, so we have to a good approximation
〈Ė〉t =
〈Ė V̂t〉t̂, (B25)
where V̂t ≡ V̂tr̃ + V̂tθ̃, cf. Eq. (B8), and similarly for L̇z
and K̇.
The explicit results we obtain using this method are
given in section III, Eqs. (3.28), (3.29), and (3.30).
2. Averaging method using the explicit
parameterization of Newtonian orbits
To perform the time-averaging using this method, we
define a parameter ξ via
1 + e cos ξ
, (B26)
where the parameter ξ varies from 0 to 2π as r̃ goes
through a complete cycle. Note that θ appears in Eqs.
(3.16) – (3.18) only in terms that are linear in Q, so we
can write θ in terms of ξ using the Newtonian relation
x3 = r cos θ = r sin ι sin(ξ + ξ0). (B27)
Here, ξ0 is the angle between the direction of the peri-
helion and the intersection of the orbital and equatorial
plane. Similarly, for the ṙθ̇ terms in Eqs. (3.17) and
(3.26) we can use the Newtonian relations ṙ = e/
p sin ξ
and ξ̇ =
p/r2. From Eqs. (2.27) and (B20) it follows
(1 + e cos ξ)2
−3 + e2 − 2e cos ξ + 2 cos2 ι(8 − e2 + 8e cos ξ + e2 cos 2ξ)
, (B28)
and from Eq. (2.12)
(1 + e cos ξ)
2 sin2 ι sin2(ξ + ξ0)− 1
. (B29)
Using these expressions, we compute the time-averaged fluxes from
〈Ė〉 =
dξ Ė (dt/dt̃) (dt̃/dξ)
dξ (dt/dt̃) (dt̃/dξ)
(B30)
and obtain:
〈Ė〉 = −32
(1− e2)3/2
e4 − S
cos(ι)
cos(2ι)
cos(2ι)
cos(2ξ0) sin
cos(2ξ0) sin
, (B31)
〈L̇z〉 = −
(1− e2)3/2
cos ι
e2 − S
2p3/2 cos ι
+ 7e2 +
cos(2ι)
−3− 45
45 + 148e2 +
cos(2ι)
1 + 3e2 +
e2 cos(2ξ0) sin
, (B32)
〈K̇〉 = −64
(1 − e2)3/2
e2 − S
2p3/2
+ 37e2 +
cos(ι)
cos(2ι)
cos(2ι)
e2 cos(2ξ0) sin
. (B33)
In the adiabatic limit, the terms involving cos(2ξ0) can
be omitted because they average to zero. As explained
by Ryan [15], the radiation reaction timescale for terms
involving ξ0 is much longer than the precession timescale
for most orbits, so the terms involving ξ0 will average
away. This is consistent with our results for the adia-
batic infinite time-averaged fluxes using the Mino time
parameter. The Mino-time averaging method was based
on the assumption that the fundamental frequencies are
incommensurate and the motion fills up the whole torus,
which is equivalent to averaging over ξ0.
[1] L. Barack and C. Cutler, Phys. Rev. D 69, 082005 (2004) [2] K. Glampedakis and S. Babak, Class. Quantum Grav.
23, 4167 (2006)
[3] D. A. Brown, et al. gr-qc/0612060
[4] E. Poisson, Living Rev. Relat. 7, 6 (2004),
http://relativity.livingreviews.org/Articles/lrr-2004-6/index.html
[5] Special Issue: Gravit. Rad. from Binary Black Holes: Ad-
vances in the perturbative approach, Class. Quant. Grav.
22 (2005)
[6] K. Glampedakis, Class. Quantum Grav. 22, S605 (2005)
[7] S. Drasco, Class. Quantum Grav. 23, S769 (2006)
[8] M. Favata and É. É. Flanagan Accuracy of adiabatic
waveforms for eccentric orbits, (in preparation)
[9] S. Drasco, É. É. Flanagan, and S. A. Hughes, Class.
Quantum Grav. 22, S801 (2005)
[10] É. É. Flanagan and T. Hinderer, Two timescale analysis
of extreme mass ratio inspirals in Kerr. II. Numerical
integration through resonances, (in preparation)
[11] N. Sago, et al., Progr. Theor. Phys. 114, 509 (2005)
[12] Y. Mino, Phys. Rev. D 67, 084027 (2003)
[13] S. A. Hughes, S. Drasco, É. É. Flanagan, and J. Franklin,
Phys. Rev. Lett. 94, 221101 (2005)
[14] F. D. Ryan, Phys. Rev. D 52, R3159 (1995)
[15] F. D. Ryan, Phys. Rev. D 53, 3064 (1996)
[16] L. Blanchet, T. Damour, G. Farese, and B. Iyer, Phys.
Rev. D 71, 124004 (2005)
[17] L. Kidder, C. Will, and A. Wiseman, Phys. Rev. D 47,
R4183 (1993)
[18] L. A. Gergely, Phys. Rev. D 61, 024035 (1999)
[19] C. Will, Phys. Rev. D 71, 084027 (2005)
[20] G. Faye, L. Blanchet, and A. Buonanno, Phys. Rev. D
74, 104033 (2006)
[21] L. Blanchet, A. Buonanno, and G. Faye, Phys. Rev. D
74, 104034 (2006)
[22] E. Poisson, Phys. Rev. D 57, 5287 (1998)
[23] L. A. Gergely and Z. Keresztes, Phys. Rev. D 67, 024020
(2003)
[24] M. Shibata, M. Sasaki, H. Tagoshi, and T. Tanaka, Phys.
Rev. D 51, 1646 (1995)
[25] É. É. Flanagan and T. Hinderer, Two timescale analy-
sis of extreme mass ratio inspirals in Kerr. I. General
formalism, (in preparation)
[26] C. Misner, K. Thorne, and J. Wheeler, Gravitation (W.H.
Freeman and Co., San Francisco, 1973)
[27] R. O. Hansen, J. Math. Phys. 15, 46 (1974)
[28] K. S. Thorne, Rev. Mod. Phys 52, 300 (1980)
[29] T. Bäckdahl, gr-qc/0612043 (2006); Class. Quantum
Grav. 22, 3585 (2005)
[30] C. Li and G. Lovelace, gr-qc/0702146 (2007)
[31] L. Blanchet and T. Damour, Phys. Lett. 104A, 82 (1984)
[32] S. A. Hughes, Phys. Rev. D 61, 084004 (2000)
[33] K. Glampedakis, S. Hughes, and D. Kennefick, Phys.
Rev. D 66, 064005 (2002)
[34] G. Arfken, Mathematical Methods for Physicists (Aca-
demic Press, CA, 1985)
[35] V. K. Melnikov, Trans. Moscow Math. Soc. 12, 1-56
(1956)
[36] W. Schmidt, Celestial mechanics in Kerr spacetime
Class. Quantum Grav. 19 (2002) 2743-2764
[37] H. Goldstein, C. Poole, and J. Safko, Classical Mechanics
(Addison-Wesley, 2002)
[38] Ý. Birol and A. Hacinliyan, Phys. Rev. E 52, 4750 (1995)
[39] H. Yoshida, Comm. Math. Phys. 116, 529 (1988)
http://arxiv.org/abs/gr-qc/0612060
http://relativity.livingreviews.org/Articles/lrr-2004-6/index.html
http://arxiv.org/abs/gr-qc/0612043
http://arxiv.org/abs/gr-qc/0702146

  We analyze the effect of gravitational radiation reaction on generic orbits
around a body with an axisymmetric mass quadrupole moment Q to linear order in
Q, to the leading post-Newtonian order, and to linear order in the mass ratio.
This system admits three constants of the motion in absence of radiation
reaction: energy, angular momentum, and a third constant analogous to the
Carter constant. We compute instantaneous and time-averaged rates of change of
these three constants. For a point particle orbiting a black hole, Ryan has
computed the leading order evolution of the orbit's Carter constant, which is
linear in the spin. Our result, when combined with an interaction quadratic in
the spin (the coupling of the black hole's spin to its own radiation reaction
field), gives the next to leading order evolution. The effect of the
quadrupole, like that of the linear spin term, is to circularize eccentric
orbits and to drive the orbital plane towards antialignment with the symmetry
axis. In addition we consider a system of two point masses where one body has a
single mass multipole or current multipole. To linear order in the mass ratio,
to linear order in the multipole, and to the leading post-Newtonian order, we
show that there does not exist an analog of the Carter constant for such a
system (except for the cases of spin and mass quadrupole). With mild additional
assumptions, this result falsifies the conjecture that all vacuum, axisymmetric
spacetimes posess a third constant of geodesic motion.

For reference, the following erratum corrects the published version of the paper. These errors have been fixed in
this arxiv-version (the article starting on page 2 has the corrected expressions).
Erratum: Evolution of the Carter constant for inspirals into a black hole: Effect of the
black hole quadrupole
[Phys. Rev. D 75, 124007 (2007)]
Éanna É. Flanagan, Tanja Hinderer
In Eqs. (3.16), (3.17), (3.18), (3.24), (3.25) and (3.26) of this paper, the variable r should be replaced everywhere
by the variable r̃, and the variable θ should be replaced everywhere by the variable θ̃. The definitions of r̃ and θ̃ are
given in Eq. (2.11). These replacements do not affect the any of the subsequent results in the paper.
Also, the right hand side of Eq. (B3) is missing a term −4SLzr̃ and Eq. (2.24) is missing a factor of dϕ/dt̃ in front
of Q.
Some terms are missing in Eqs. (3.18), (3.26) and (3.30) - (3.33). The additional terms in Eqs. (3.18) and (3.26)
15r̃7
−75K2 + 2Kr̃(51r̃E + 50) + 8r̃2(r̃E + 1)(3r̃E + 5)
15p2r̃7
25p3(3p− 4r̃) + p2r̃2
11− 51e2
+ 32pr̃3
1− e2
+ 6r̃4
1− e2
respectively. These result in additional fractional corrections to Eq. (3.30) given by
and the full expression replacing the O(Q) terms in Eq. (3.30) is then
〈K̇〉 = −
(1− e2)3/2
cos(2ι)
+O(S), O(S2)− terms.
Equations (3.31), (3.32) and (3.33) contain typos in the O(S) and O(Q) terms, the corrected expressions are given
below. We thank P. Komorowski for pointing this out. Equation (3.31) should be replaced by
〈ṗ〉 = −64
(1− e2)3/2
− S cos(ι)
96p3/2
1064 + 1516e2 + 475e4
149e2
469e2
227e4
cos(2ι)
+ e2 +
[13− cos(2ι)]
, (0.1)
Equation (3.32) should be replaced by
〈ė〉 = −304
e(1− e2)3/2
121e2
Se(1− e2)3/2 cos(ι)
5p11/2
1172 + 932e2 +
1313e4
Q(1− e2)3/2
785e2
− 219e
+ 13e6 +
2195e2
+ 251e4 +
218e6
cos(2ι)
2e(1− e2)3/2
2 + 3e2 +
[13− cos(2ι)] , (0.2)
and the corrected Eq. (3.33) is
〈ι̇〉 = S sin(ι)(1 − e
2)3/2
p11/2
1− e2
S2 sin(2ι)
240p6
8 + 3e2
8 + e2
Q cot(ι)(1 − e2)3/2
312 + 736e2 − 83e4 −
408 + 1268e2 + 599e4
cos(2ι)
. (0.3)
http://arxiv.org/abs/0704.0389v8
Evolution of the Carter constant for inspirals into a black hole: effect of the black hole
quadrupole
Éanna É. Flanagan1,2 and Tanja Hinderer1
Center for Radiophysics and Space Research, Cornell University, Ithaca, NY 14853, USA
Laboratory for Elementary Particle Physics, Cornell University, Ithaca, NY 14853, USA
(Dated: November 4, 2018)
We analyze the effect of gravitational radiation reaction on generic orbits around a body with an
axisymmetric mass quadrupole moment Q to linear order in Q, to the leading post-Newtonian order,
and to linear order in the mass ratio. This system admits three constants of the motion in absence
of radiation reaction: energy, angular momentum along the symmetry axis, and a third constant
analogous to the Carter constant. We compute instantaneous and time-averaged rates of change
of these three constants. For a point particle orbiting a black hole, Ryan [15] has computed the
leading order evolution of the orbit’s Carter constant, which is linear in the spin. Our result, when
combined with an interaction quadratic in the spin (the coupling of the black hole’s spin to its own
radiation reaction field), gives the next to leading order evolution. The effect of the quadrupole, like
that of the linear spin term, is to circularize eccentric orbits and to drive the orbital plane towards
antialignment with the symmetry axis.
In addition we consider a system of two point masses where one body has a single mass multipole
or current multipole of order l. To linear order in the mass ratio, to linear order in the multipole,
and to the leading post-Newtonian order, we show that there does not exist an analog of the Carter
constant for such a system (except for the cases of an l = 1 current moment and an l = 2 mass
moment). Thus, the existence of the Carter constant in Kerr depends on interaction effects between
the different multipoles. With mild additional assumptions, this result falsifies the conjecture that
all vacuum, axisymmetric spacetimes posess a third constant of the motion for geodesic motion.
PACS numbers: 04.25.Nx, 04.30.Db
I. INTRODUCTION AND SUMMARY
The inspiral of stellar mass compact objects with
masses µ in the range µ ∼ 1 − 100M⊙ into massive
black holes with masses M ∼ 105 − 107M⊙ is one of
the most important sources for the future space-based
gravitational wave detector LISA. Observing such events
will provide a variety of information: (i) the masses and
spins of black holes can be measured to high accuracy
(∼ 10−4); which can constrain the black hole’s growth
history [1]; (ii) the observations will give a precise test
of general relativity in the strong field regime and unam-
biguously identify whether the central object is a black
hole [2]; and (iii) the measured event rate will give in-
sight into the complex stellar dynamics in galactic nu-
clei [1]. Analogous inspirals may also be interesting for
the advanced stages of ground-based detectors: it has
been estimated that advanced LIGO could detect up to
∼ 10 − 30 inspirals per year of stellar mass compact
objects into intermediate mass black holes with masses
M ∼ 102 − 104M⊙ in globular clusters [3]. Detect-
ing these inspirals and extracting information from the
datastream will require accurate models of the gravita-
tional waveform as templates for matched filtering. For
computing templates, we therefore need a detailed un-
derstanding of the how radiation reaction influences the
evolution of bound orbits around Kerr black holes [4–7].
There are three dimensionless parameters characteriz-
ing inspirals of bodies into black holes:
• the dimensionless spin parameter a = |S|/M2 of
the black hole, where S is the spin.
• the strength of the interaction potential ǫ2 =
GM/rc2, i.e. the expansion parameter used in post-
Newtonian (PN) theory.
• the mass ratio µ/M .
For LISA data analysis we will need waveforms that are
accurate to all orders in a and ǫ2, and to leading order
in µ/M . However, it is useful to have analytic results in
the regimes a ≪ 1 and/or ǫ2 ≪ 1. Such approximate
results can be useful as a check of numerical schemes
that compute more accurate waveforms, for scoping out
LISA’s data analysis requirements [1, 6], and for assessing
the accuracy of the leading order in µ/M or adiabatic
approximation [8–10]. There is substantial literature on
such approximate analytic results, and in this paper we
will extend some of these results to higher order.
A long standing difficulty in computing the evolution
of generic orbits has been the evolution of the orbit’s
”Carter constant”, a constant of motion which governs
the orbital shape and inclination. A theoretical prescrip-
tion now exists for computing Carter constant evolution
to all orders in ǫ and a in the adiabatic limit µ ≪ M
[9, 11–13], but it has not yet been implemented numer-
ically. In this paper we focus on computing analytically
the evolution of the Carter constant in the regime a≪ 1,
ǫ ≪ 1, µ/M ≪ 1, extending earlier results by Ryan
[14, 15].
We next review existing analytical work on the effects
of multipole moments on inspiral waveforms. For non-
spinning point masses, the phase of the l = 2 piece of
the waveform is known to O(ǫ7) beyond leading order
[16], while spin corrections are not known to such high
order. To study the leading order effects of the central
body’s multipole moments on the inspiral waveform, in
the test mass limit µ ≪ M , one has to correct both the
conservative and dissipative pieces of the forces on the
bodies. For the conservative pieces, it suffices to use the
Newtonian action for a binary with an additional multi-
pole interaction potential. For the dissipative pieces, the
multipole corrections to the fluxes at infinity of the con-
served quantities can simply be added to the known PN
point mass results. The lowest order spin-orbit coupling
effects on the gravitational radiation were first derived by
Kidder [17], then extended by Ryan [14, 15], Gergely [18],
and Will [19]. Recently, the corrections of O(ǫ2) beyond
the leading order to the spin-orbit effects on the fluxes
were derived [20, 21]. Corrections to the waveform due
to the quadrupole - mass monopole interaction were first
considered by Poisson [22], who derived the effect on the
time averaged energy flux for circular equatorial orbits.
Gergely [23] extended this work to generic orbits and
computed the radiative instantaneous and time averaged
rates of change of energy E, magnitude of angular mo-
mentum |L|, and the angle κ = cos−1(S ·L) between the
spin S and orbital angular momentum L. Instead of the
Carter constant, Gergely identified the angular average
of the magnitude of the orbital angular momentum, L̄, as
a constant of motion. The fact that to post-2-Newtonian
(2PN) order there is no time averaged secular evolution of
the spin allowed Gergely to obtain expressions for L̇ and
κ̇ from the quadrupole formula for the evolution of the
total angular momentum J = L+S. In a different paper,
Gergely [18] showed that in addition to the quadrupole,
self-interaction spin effects also contribute at 2PN order,
which was seen previously in the black hole perturbation
calculations of Shibata et al. [24]. Gergely calculated
the effect of this interaction on the instantaneous and
time-averaged fluxes of E and |L| but did not derive the
evolution of the third constant of motion.
In this paper, we will re-examine the effects of the
quadrupole moment of the black hole and of the leading
order spin self interaction. For a black hole, our analysis
will thus contain all effects that are quadratic in spin to
the leading order in ǫ2 and in µ/M . Our work will extend
earlier work by
• Considering generic orbits.
• Using a natural generalization of the Carter-type
constant that can be defined for two point particles
when one of them has a quadrupole. This facilitates
applying our analysis to Kerr inspirals.
• Computing instantaneous as well as time-averaged
fluxes for all three constants of motion: energy
E, z-component of angular momentum Lz, and
Carter-type constant K. For most purposes, only
time-averaged fluxes are needed as only they are
gauge invariant and physically relevant. However,
there is one effect for which the time-averaged
fluxes are insufficient, namely transient resonances
that occur during an inspiral in Kerr in the vicin-
ity of geodesics for which the radial and azimuthal
frequencies are commensurate [10, 25]. The instan-
taneous fluxes derived in this paper will be used in
[10] for studying the effect of these resonances on
the gravitational wave phasing.
We will analyze the effect of gravitational radiation re-
action on orbits around a body with an axisymmetric
mass quadrupole moment Q to leading order in Q, to the
leading post-Newtonian order, and to leading order in the
mass ratio. With these approximations the adiabatic ap-
proximation holds: gravitational radiation reaction takes
place over a timescale much longer than the orbital pe-
riod, so the orbit looks geodesic on short timescales. We
follow Ryan’s method of computation [14]: First, we cal-
culate the orbital motion in the absence of radiation re-
action and the associated constants of motion. Next, we
use the leading order radiation reaction accelerations that
act on the particle (given by the Burke-Thorne formula
[26] augmented by the relevant spin corrections [14]) to
compute the evolution of the constants of motion. In the
adiabatic limit, the time-averaged rates of change of the
constants of motion can be used to infer the secular or-
bital evolution. Our results show that a mass quadrupole
has the same qualitative effect on the evolution as spin: it
tends to circularize eccentric orbits and drive the orbital
plane towards antialignment with the symmetry axis of
the quadrupole.
The relevance of our result to point particles inspi-
ralling into black holes is as follows. The vacuum space-
time geometry around any stationary body is completely
characterized by the body’s mass multipole moments
IL = Ia1,a2...al and current multipole moments SL =
Sa1,a2...al [27]. These moments are defined as coefficients
in a power series expansion of the metric in the body’s
local asymptotic rest frame [28]. For nearly Newtonian
sources, they are given by integrals over the source as
IL ≡ Ia1,...al =
ρx<a1 . . . xal>d
3x, (1.1)
SL ≡ Sa1,...al =
ρxpvqǫpq<a1xa2 . . . xal>d
3x.(1.2)
Here ρ is the mass density and vq is the velocity, and ”<
· · · >” means ”symmetrize and remove all traces”. For
axisymmetric situations, the tensor multipole moments
IL (SL) contain only a single independent component,
conventionally denoted by Il (Sl) [27]. For a Kerr black
hole of mass M and spin S, these moments are given by
Il + iSl =M
l+1(ia)l, (1.3)
where a is the dimensionless spin parameter defined by
a = |S|/M2. Note that Sl = 0 for even l and Il = 0 for
odd l.
Consider now inspirals into an axisymmetric body
which has some arbitrary mass and current multipoles
Il and Sl. Then we can consider effects that are linear in
Il and Sl for each l, effects that are quadratic in the mul-
tipoles proportional to IlIl′ , IlSl′ , SlSl′ , effects that are
cubic, etc. For a general body, all these effects can be sep-
arated using their scalings, but for a black hole, Il ∝ al
for even l and Sl ∝ al for odd l [see Eq.(1.3)], so the ef-
fects cannot be separated. For example, a physical effect
that scales as O(a2) could be an effect that is quadratic
in the spin or linear in the quadrupole; an analysis in
Kerr cannot distinguish these two possibilities. For this
reason, it is useful to analyze spacetimes that are more
general than Kerr, characterized by arbitrary Il and Sl,
as we do in this paper. For recent work on computing
exact metrics characterized by sets of moments Il and Sl,
see Refs. [29, 30] and references therein.
The leading order effect of the black hole’s multipoles
on the inspiral is the O(a) effect computed by Ryan [15].
This O(a) effect depends linearly on the spin S1 and is
independent of the higher multipoles Sl and Il since these
all scale as O(a2) or smaller. In this paper we compute
the O(a2) effect on the inspiral, which includes the lead-
ing order linear effect of the black hole’s quadrupole (lin-
ear in I2 ≡ Q) and the leading order spin self-interaction
(quadratic in S1).
We next discuss how these O(a2) effects scale with the
post-Newtonian expansion parameter ǫ. Consider first
the conservative orbital dynamics. Here it is easy to see
that fractional corrections that are linear in I2 scale as
O(a2ǫ4), while those quadratic in S1 scale as O(a
2ǫ6).
Thus, the two types of terms cleanly separate. We com-
pute only the leading order, O(a2ǫ4), term. For the dissi-
pative contributions to the orbital motion, however, the
scalings are different. There are corrections to the radi-
ation reaction acceleration whose fractional magnitudes
are O(a2ǫ4) from both types of effects linear in I2 and
quadratic in S1. The effects quadratic in S1 are due to
the backscattering of the radiation off the piece of space-
time curvature due to the black hole’s spin. This effect
was first pointed out by Shibata et al. [24], who com-
puted the time-averaged energy flux for circular orbits
and small inclination angles based on a PN expansion of
black hole perturbations. Later, Gergely [18] analyzed
this effect on the instantaneous and time-averaged fluxes
of energy and magnitude of orbital angular momentum
within the PN framework.
The organization of this paper is as follows. In Sec.
II, we study the conservative orbital dynamics of two
point particles when one particle is endowed with an ax-
isymmetric quadrupole, in the weak field regime, and to
leading order in the mass ratio. In Sec. III, we com-
pute the radiation reaction accelerations and the instan-
taneous and time-averaged fluxes. In order to have all
the contributions at O(a2ǫ4) for a black hole, we include
in our computations of radiation reaction acceleration
the interaction that is quadratic in the spin S1. The ap-
plication to black holes in Sec. IV briefly discusses the
qualitative predictions of our results and also compares
with previous results.
The methods used in this paper can be applied only
to the black hole spin (as analyzed by Ryan [14]) and
the black hole quadrupole (as analyzed here). We show
in Sec. V that for the higher order mass and current
multipole moments taken individually, an analog of the
Carter constant cannot be defined to the order of our
approximations. We then show that under mild assump-
tions, this non-existence result can be extended to exact
spacetimes, thus falsifying the conjecture that all vac-
uum axisymmetric spacetimes possess a third constant
of geodesic motion.
II. EFFECT OF AN AXISYMMETRIC MASS
QUADRUPOLE ON THE CONSERVATIVE
ORBITAL DYNAMICS
Consider two point particles m1 and m2 interacting in
Newtonian gravity, where m2 ≪ m1 and where the mass
m1 has a quadrupole moment Qij which is axisymmetric:
Qij =
d3xρ(r)
xixj −
r2δij
(2.1)
ninj −
. (2.2)
For a Kerr black hole of mass M and dimensionless spin
parameter a with spin axis along n, the quadrupole scalar
is Q = −M3a2.
The action describing this system, to leading order in
m2/m1, is
µv2 − µΦ(r)
, (2.3)
where v = ṙ is the velocity, the potential is
Φ(r) = −M
xixjQij , (2.4)
µ is the reduced mass and M the total mass of the bi-
nary, and we are using units with G = c = 1. We work to
linear order in Q, to linear order in m2/m1, and to lead-
ing order in M/r. In this regime, the action (2.3) also
describes the conservative effect of the black hole’s mass
quadrupole on bound test particles in Kerr, as discussed
in the introduction. We shall assume that the quadrupole
Qij is constant in time. In reality, the quadrupole will
evolve due to torques that act to change the orientation
of the central body. An estimate based on treating m1 as
a rigid body in the Newtonian field ofm2 gives the scaling
of the timescale for the quadrupole to evolve compared to
the radiation reaction time as (see Appendix I for details)
Tevol
(2.5)
Here, we have denoted the dimensionless spin and
quadrupole of the body by S̄ and Q̄ respectively, and
the last relation applies for a Kerr black hole. Since
µ/M ≪ 1, the first factor in Eq. (2.5) will be large, and
since 1/a ≥ 1 and for the relativistic regime M/r ∼ 1,
the evolution time is long compared to the radiation re-
action time. Therefore we can neglect the evolution of
the quadrupole at leading order.
This system admits three conserved quantities, the en-
µv2 + µΦ(r), (2.6)
the z-component of angular momentum
Lz = ez · (µr× v), (2.7)
and the Carter-type constant
K = µ2(r× v)2 − 2Qµ
(n · r)2
(n · v)2 −
. (2.8)
(See below for a derivation of this expression for K).
A. Conservative orbital dynamics in a
Boyer-Lindquist-like coordinate system
We next specialize to units where M = 1. We also
define the rescaled conserved quantities by Ẽ = E/µ,
L̃z = Lz/µ, K̃ = K/µ
2, and drop the tildes. These spe-
cializations and definitions have the effect of eliminating
all factors of µ and M from the analysis. In spherical
polar coordinates (r, θ, ϕ) the constants of motion E and
Lz become
(ṙ2 + r2θ̇2 + r2 sin2 θϕ̇2)−
(1− 3 cos2 θ), (2.9)
Lz = r
2 sin2 θϕ̇. (2.10)
In these coordinates, the Hamilton-Jacobi equation is not
separable, so a separation constant K cannot readily be
derived. For this reason we switch to a different coordi-
nate system (r̃, θ̃, ϕ) defined by
r cos θ = r̃ cos θ̃
r sin θ = r̃ sin θ̃
. (2.11)
We also define a new time variable t̃ by
cos(2θ̃)
dt̃. (2.12)
The action (2.3) in terms of the new variables to linear
order in Q is
r̃2 sin2 θ̃
sin2 θ̃
. (2.13)
However, a difficulty is that the action (2.13) does not
give the same dynamics as the original action (2.3). The
reason is that for solutions of the equations of motion for
the action (2.3), the variation of the action vanishes for
paths with fixed endpoints for which the time interval ∆t
is fixed. Similarly, for solutions of the equations of motion
for the action (2.13), the variation of the action vanishes
for paths with fixed endpoints for which the time interval
∆t̃ is fixed. The two sets of varied paths are not the
same, since ∆t 6= ∆t̃ in general. Therefore, solutions of
the Euler-Lagrange equations for the action (2.3) do not
correspond to solutions of the Euler-Lagrange equations
for the action (2.13). However, in the special case of zero-
energy motions, the extra terms in the variation of the
action vanish. Thus, a way around this difficulty is to
modify the original action to be
µv2 − µΦ(r) + E
. (2.14)
This action has the same extrema as the action (2.3),
and for motion with physical energy E, the energy com-
puted with this action is zero. Transforming to the new
variables yields, to linear order in Q:
r̃2 sin2 θ̃
sin2 θ̃
+ E − QE
cos(2θ̃)
. (2.15)
The zero-energy motions for this action coincide with the
zero energy motions for the action (2.14). We use this
action (2.15) as the foundation for the remainder of our
analysis in this section.
The z-component of angular momentum in terms of
the new variables (r̃, θ̃, ϕ, t̃) is
Lz = r̃
2 sin2 θ̃
sin2 θ̃
. (2.16)
We now transform to the Hamiltonian:
p2r̃ −
− E − Q
sin2 θ̃
+QE cos(2θ̃)
(2.17)
and solve the Hamiltonian Jacobi equation. Denoting
the separation constant by K we obtain the following
two equations for the r̃ and θ̃ motions:
= 2E +
, (2.18)
= K − L
sin2 θ̃
−QE cos(2θ̃). (2.19)
Note that the equations of motion (2.18) and (2.19) have
the same structure as the equations of motion for Kerr
geodesic motion. Using Eqs. (2.18), (2.19) and (2.16)
together with the inverse of the transformation (2.11)
to linear order in Q, we obtain the expression for K in
spherical polar coordinates:
K = r4(θ̇2 + sin2 θϕ̇2) +Q(ṙ cos θ − rθ̇ sin θ)2 + Q
(ṙ2 + r2 θ̇2 + r2 sin2 θϕ̇2)− 2Q
cos2 θ. (2.20)
This is equivalent to the formula (2.8) quoted earlier.
B. Effects linear in spin on the conservative orbital
dynamics
To include the linear in spin effects, we repeat Ryan’s
analysis [14, 15] (he only gives the final, time averaged
fluxes; we will also give the instantaneous fluxes). We
can simply add these linear in spin terms to our results
because any terms of order O(SQ) will be higher than
the order a2 to which we are working. The correction to
the action (2.3) due to spin-orbit coupling is
Sspin−orbit =
−2µSn
iǫijkxj ẋk
. (2.21)
We will restrict our analysis to the case when the unit
vectors ni corresponding to the axisymmetric quadrupole
Qij and to the spin Si coincide, as they do in Kerr.
Including the spin-orbit term in the action (2.3) results
in the following modified expressions for Lz and K:
Lz = n · (µr× v)−
[r2 − (n · r)2], (2.22)
K = (r× v)2 − 4S
n · (r× v)− 2Q
(n · r)2
(n · v)2 − 1
. (2.23)
In terms of the Boyer-Lindquist like coordinates, the con-
served quantities with the linear in spin terms included
Lz = r̃
2 sin2 θ̃
sin2 θ̃ −Q sin4 θ̃
(2.24)
K = r4(θ̇2 + sin2 θϕ̇2)− 4Sr sin2 θϕ̇
cos2 θ +Q(ṙ cos θ − rθ̇ sin θ)2 + QM
(ṙ2 + r2θ̇2 + r2 sin2 θϕ̇2). (2.25)
The equations of motion are
= 2E+
− 4SLz
, (2.26)
= K −
sin2 θ̃
−QE cos(2θ̃). (2.27)
III. EFFECTS LINEAR IN QUADRUPOLE AND
QUADRATIC IN SPIN ON THE EVOLUTION OF
THE CONSTANTS OF MOTION
A. Evaluation of the radiation reaction force
The relative acceleration of the two bodies can be writ-
ten as
a = −∇Φ(r) + arr, (3.1)
where arr is the radiation-reaction acceleration. Combin-
ing this with Eqs. (2.6), (2.22) and (2.23) for E, Lz and
K gives the following formulae for the time derivatives of
the conserved quantities:
Ė = v · arr, (3.2)
L̇z = n · (r× arr), (3.3)
K̇ = 2(r× v) · (r× arr)−
n · (r× arr)
+2Q(n · v) (n · arr)−Qv · arr. (3.4)
The standard expression for the leading order radiation
reaction acceleration acting on one of the bodies is [31]:
ajrr = −
jk xk +
ǫjpqS
pk xkxq +
ǫjpqS
pk xkvq
ǫpq[jS
xqvk. (3.5)
Here the superscripts in parentheses indicate the number
of time derivatives and square brackets on the indices
denote antisymmetrization.
The multipole moments Ijk(t) and Sjk(t) in Eq. (3.5)
are the total multipole moments of the spacetime, i.e.
approximately those of the black hole plus those due to
the orbital motion. The expression (3.5) is formulated
in asymptotically Cartesian mass centered (ACMC) co-
ordinates of the system, which are displaced from the
coordinates used in Sec. II by an amount [28]
δr(t) = − µ
r(t). (3.6)
This displacement contributes to the radiation reaction
acceleration in the following ways:
1. The black hole multipole moments Il and Sl, which
are time-independent in the coordinates used in
Sec. II, will be displaced by δr and thus will con-
tribute to the (l + 1)th ACMC radiative multipole
[28].
2. The constants of motion are defined in terms of the
black hole centered coordinates used in Sec. II, so
the acceleration arr we need in Eqs. (3.2) – (3.4)
is the relative acceleration. This requires calculat-
ing the acceleration of both the black hole and the
point mass in the ACMC coordinates using (3.5),
and then subtracting to find arr = a
rr − aMrr [14].
To leading order in µ, the only effect of the acceler-
ation of the black hole is via a backreaction of the
radiation field: the lth black hole moments couple
to the (l+1)th radiative moments, thus producing
an additional contribution to the acceleration.
For our calculations at O(S1ǫ
3), O(I2ǫ
4), O(S21ǫ
4), we
can make the following simplifications:
• quadrupole corrections: The fractional corrections
linear in I2 = Q that scale as O(a
2ǫ4) require only
the effect of I2 on the conservative orbital dynamics
as computed in Sec. IIA and the Burke-Thorne for-
mula for the radiation reaction acceleration [given
by the first term in Eq. (3.5)].
• spin-spin corrections: As discussed in the intro-
duction, the fractional corrections quadratic in S1
to the conservative dynamics scale as O(a2ǫ6) and
are subleading order effects which we neglect. At
O(a2ǫ4), the only effect quadratic in S1 is the
backscattering of the radiation off the spacetime
curvature due to the spin. As discussed in item 1.
above, the black hole’s current dipole Si = S1δi3
(taking the z-axis to be the symmetry axis) will
contribute to the radiative current quadrupole an
amount
ij = −
S1xiδj3. (3.7)
The black hole’s current dipole Si will couple to
the gravitomagnetic radiation field due to Sij as
discussed in item 2. above, and contribute to the
relative acceleration as [14]:
aj spinrr =
S1δi3S
ij . (3.8)
For our purposes of computing terms quadratic in
the spin, we substitute S
ij for Sij in Eq. (3.8).
Evaluating these quadratic in spin terms requires
only the Newtonian conservative dynamics, i.e. the
results of Sec. II and Eqs. (3.2) – (3.4) with the
quadrupole set to zero.
• linear in spin corrections: Contributions to these
effects are from Eq. (3.5) with the current
quadrupole replaced by just the spin contribution
(3.7), and from Eq. (3.8) evaluated using only the
orbital current quadrupole.
With these simplifications, we replace the expression
(3.5) for the radiation reaction acceleration with
ajrr = −
jk xk +
ǫjpqS
(6) spin
pk xkxq
ǫjpqS
(5) spin
pk xkvq +
ǫpq[jS
(5) spin
S1δi3
(5) orbit
ij + S
(5) spin
. (3.9)
To justify these approximations, consider the scaling of
the contribution of black hole’s acceleration to the orbital
dynamics. The mass and current multipoles of the black
hole contribute terms to the Hamiltonian that scale with
∆H ∼ Slǫ2l+3 & Ilǫ2l+2. (3.10)
Since the Newtonian energy scales as ǫ2, the fractional
correction to the orbital dynamics scale as
∆H/E ∼ Slǫ2l+1 & Ilǫ2l. (3.11)
To O(ǫ4), the only radiative multipole moments that con-
tribute to the acceleration (3.5) are the mass quadrupole
I2, the mass octupole I3, and the current quadrupole S2
(cf. [17]). Since we are focusing only on the leading or-
der terms quadratic in spin (these can simply be added
to the known 2PN point particle and 1.5PN linear in spin
results), the only terms in Eq. (3.5) relevant for our pur-
poses are those given in Eq. (3.9). The results from a
computation of the fully relativistic metric perturbation
for black hole inspirals [24] show that quadratic in spin
corrections to the l = 2 piece compared to the flat space
Burke-Thorne formula first appear at O(a2ǫ4), which is
consistent with the above arguments.
B. Instantaneous fluxes
We evaluate the radiation reaction force as follows.
The total mass and current quadrupole moment of the
system are
QTij = Qij + µxixj , (3.12)
STij = S
ij + xiǫjkmxkẋm, (3.13)
where from Eq. (2.11)
r̃ sin θ̃
cosϕ, r̃ sin θ̃
sinϕ,
r̃ cos θ̃
. (3.14)
Only the second term in Eq. (3.12) contributes to the
time derivative of the quadrupole. We differentiate five
times by using
cos(2θ̃)
, (3.15)
to the order we are working as discussed above. Af-
ter each differentiation, we eliminate any occurrences of
dϕ/dt̃ using Eq. (2.24), and we eliminate any occurrences
of the second order time derivatives d2r̃/dt̃2 and d2θ̃/dt̃2
in favor of first order time derivatives using (the time
derivatives of) Eqs. (2.26) and (2.27). For computing the
terms linear and quadratic in S1, we set the quadrupole
Q to zero in all the formulae. We insert the resulting ex-
pression into the formula (3.9) for the self-acceleration,
and then into Eqs. (3.2) – (3.4). We eliminate (dr̃/dt̃)2,
(dθ̃/dt̃)2, and (dϕ/dt̃) in favor of E, Lz, and K using
Eqs. (2.24) – (2.27). In the final expressions for the in-
stantaneous fluxes, we keep only terms that are of O(S),
O(Q) and O(S2) and obtain the following results:
15r̃4
− 40K
272KE
196K2 +
r̃2 − 3668
Kr̃ − 352KEr̃2 + 1024
Er̃3 +
E2r̃4
−49K2 − 169KL2z + r̃
+ 2r̃2
+ 47KE +
− 152
r̃3E − 16r̃4E2
−562K2 +
Kr̃ −
r̃2 +
KEr̃2 −
r̃3E − 160r̃4E2
cos(2θ̃)
sin(2θ̃)
439K − 926
r̃ − 1528
θ ˙̃r
−K2 + 22
Kr̃ − 28
r̃2 +
KEr̃2 − 236
r̃3E − 32
r̃4E2
cos(2θ̃)− r̃3 sin(2θ̃)
θ ˙̃r
−49K2 + 6KL2z + 2r̃
63K − 16
L2z −
+ r̃2
112KE − 48
− 1652
r̃3E − 224
r̃4E2
, (3.16)
L̇z =
144LzE
− 24KLz
−50K2 + 240KL2z +
Kr̃ − 7376
L2z r̃ +
r̃2 + 56KEr̃2 − 1824
EL2z r̃
Er̃3 +
E2r̃4
50K2 − 62
Kr̃ − 316
r̃2 − 56KEr̃2 − 624
Er̃3 − 128
E2r̃4
cos(2θ̃)
−104K + 64r̃ + 64Er̃2
sin(2θ̃) ˙̃r
660Er̃2 + 753r̃− 360L2z − 435K +
1601r̃+ 1512r̃2E − 1185K
cos 2θ̃
174QLz
sin(2θ̃) ˙̃r
2S2Lz
Er̃2 + 16r̃ − 9K
, (3.17)
20r̃ + 18r̃2E − 15K
280K2 − 14008
Kr̃ +
r̃2 +
Er̃3 − 2528
KEr̃2 +
E2r̃4
−45K2 + r̃L2z(83 + 80r̃E)− 115KL2z + 14Kr̃(6 + 5r̃E)
15r̃7
cos(2θ̃)
−2175K2 + 2975Kr̃+ 80r̃2 + 3012KEr̃2 − 112Er̃3 − 168E2r̃4
15r̃4
3075K − 20r̃ − 192Er̃2
sin(2θ̃)
θ ˙̃r
7K − 2L2z
−3K + 16
+K cos(2θ̃)
3K − 16
r̃ − 24
sin(2θ̃)
−4K +
θ ˙̃r. (3.18)
C. Alternative set of constants of the motion
A body in a generic bound orbit in Kerr traces an
open ellipse precessing about the hole’s spin axis. For
stable orbits the motion is confined to a toroidal region
whose shape is determined by E, Lz, K. The motion
can equivalently be characterized by the set of constants
inclination angle ι, eccentricity e, and semi-latus rectum
p defined by Hughes [32]. The constants ι, p and e are
defined by cos ι = Lz/
K, and by r̃± = p/(1± e), where
r̃± are the turning points of the radial motion, and r̃
is the Boyer-Lindquist radial coordinate. This param-
eterization has a simple physical interpretation: in the
Newtonian limit of large p, the orbit of the particle is an
ellipse of eccentricity e and semilatus rectum p on a plane
whose inclination angle to the hole’s equatorial plane is
ι. In the relativistic regime p ∼M , this interpretation of
the constants e, p, and ι is no longer valid because the
orbit is not an ellipse and ι is not the angle at which the
object crosses the equatorial plane (see Ryan [14] for a
discussion).
We adopt here analogous definitions of constants of
motion ι, e and p, namely
cos(ι) = Lz/
K, (3.19)
= r̃±. (3.20)
Here K is the conserved quantity (2.23) or (2.25), and r̃±
are the turning points of the radial motion using the r̃
coordinate defined by Eq. (2.11), given by the vanishing
of the right-hand side of Eq. (2.26).
We now rewrite our results in terms of the new con-
stants of the motion e, p and ι. We can use Eq. (2.26)
together with the equations (3.19) and (3.20) to write E,
Lz and K as functions of p, e and ι. To leading order in
Q and S we obtain
K = p
1− 2S cos ι
3 + e2
1 + e2
) 2Q cos2 ι
3 + e2
, (3.21)
E = − (1− e
2S cos ι
1− e2
1− e2
cos2 ι− 1
, (3.22)
p cos ι
1− S cos ι
(3 + e2)−
1 + e2
) Q cos2 ι
3 + e2
. (3.23)
As discussed in the introduction, the effects quadratic in
S on the conservative dynamics scale as O(a2ǫ6) and thus
are not included in this analysis to O(a2ǫ4).
Inserting these relations into the expressions (3.16)–
(3.18) gives, dropping terms of O(QS), O(Q2) and
O(QS2):
Ė = −
15p2r̃7
75p4 − 100p3r̃ + p2r̃2
11− 51e2
+ 32pr̃3
1− e2
)− 6r̃4
1− e2
4S cos ι
15p7/2r̃9
735p6 − 2751p5r̃ + 10p4r̃2(365− 6e2)− 128pr̃5(1− e2)2 − 48r̃6(e2 − 1)3
64S cos ι
15p3/2r̃6
5p(−23 + 3e2)− 3r̃(−9 + e2 + 8e4)
15p4r̃9
4005p6 − 6499p5r̃ + 2p4r̃2
1577− 1977e2
− 24r̃6
1− e2
)3 − 32p3r̃3
8− 33e2
+ 64pr̃5
1− 2e2 + e4
15p4r̃9
24p2r̃4
5− 27e2 + 22e4
− pr̃3 sin(2θ̃)
6585p2 − 4630pr̃+ 2292r̃2(1 − e2)
θ ˙̃r
15p4r̃9
2p2 cos(2θ̃)
4215p4 − 7495p3r̃ + 4p2r̃2(1151− 951e2)− 1012pr̃3(1− e2) + 300r̃4(1− 2e2 + e4)
15p4r̃9
cos(2ι)
2535p6 − 3307p5r̃ + 12p4r̃2(37− 237e2)− 48r̃6(1− e2)3 + 800p3r̃3(1 + e2) + 128pr̃5(1− 2e2 + e4)
15p2r̃5
cos(2ι)
1 + 2e2 − 3e4
15p2r̃9
84r̃4(1− e2)2(1 + e2)2 + 345p4 − 905p3r̃ − 413pr̃3(1− e2) + 2p2r̃2(446− 201e2)
15p2r̃9
cos(2θ̃)
15p4 − 110p3r̃ + 4p2r̃2(47− 12e2)− 118pr̃3(1− e2) + 24r̃4(1− e2)2(1 + e2)2
15r̃9
cos(2ι)
45p2 − 80pr̃ + 36r̃2(1− e2)
15pr̃6
sin(2θ̃) ˙̃r
15p2 + 10pr̃ − 12r̃2(1− e2)
, (3.24)
L̇z = −
8 cos ι
15p2 − 20pr̃ + 9r̃2(1 − e2)
15p2r̃7
525p4 − 1751p3r̃ + 34p2r̃2(61− 6e2) + 12pr̃3(−69 + 29e2) + 6r̃4(17 + 2e2 − 19e4)
15p2r̃7
375p4 − 93p3r̃ + 468pr̃3(1− e2)− 10p2r̃2(58 + 21e2)− 48r̃4(1 − 2e2 + e4)
cos(2θ̃)
15p2r̃7
450p4 − 922p3r̃ − 60pr̃3(3 + e2)− 9p2r̃2(−83 + 23e2) + 27r̃4(1 + 2e2 − 3e4)
cos(2ι)
13p2 − 8pr̃ + 4r̃2(1− e2)
sin(2θ̃) ˙̃r
− Q cos ι
5p5/2r̃7
615p4 − 753p3r̃ + 15p2r̃2
19− 31e2
+ 20pr̃3
1 + 3e2
+ 9r̃4
1− 6e2 + 5e4
− Q cos ι
5p1/2r̃7
cos(2θ̃)
1185p2 − 1601pr̃+ 756r̃2(1− e2)
− 2Q cos ι
5p5/2r̃7
2 cos(2ι)
45p4 − 18r̃4e2(1− e2)− 45p2r̃2(1 + e2) + 20pr̃3(1 + e2)
− 435p3r̃3 sin(2θ̃) ˙̃θ ˙̃r
2 cos ι
p1/2r̃7
9p2 − 16pr̃ + 36
r̃2(1− e2)
, (3.25)
20pr̃ − 15p2 − 9r̃2(1− e2)
8S cos ι
15p3/2r̃7
525p4 − 1751p3r̃ + 2p2r̃2(1172− 57e2) + 12pr̃3(−99 + 19e2)− 24r̃4(−11 + 4e2 + 7e4)
5p2r̃7
−615p4 + 753p3r̃ + 30p2r̃2(17e2 − 9) + 72r̃4e2(1− e2)− 40pr̃3(1 + 3e2)
5p2r̃7
cos(2ι)
−345p4 + 249p3r̃ − 160pr̃3(1 + e2) + 120p2r̃2(1 + 3e2) + 36r̃4(1 + 2e2 − 3e4)
15p2r̃7
2 cos(2θ̃)
2175p4 − 2975p3r̃ − 56pr̃3(1 − e2) + 2p2r̃2(713− 753e2) + 42r̃4(1 − 2e2 + e4)
15pr̃4
sin(2θ̃)
3075p2 − 20pr̃ + 96r̃2(1 − e2)
−9p2 + 16pr̃ − 36
r̃2(1− e2)
cos(2θ̃) + cos(2ι)
3p2 − 16
pr̃ +
r̃2(1 − e2)
sin(2θ̃) ˙̃r
−2p2 + 7
pr̃ − 4
r̃2(1− e2)
. (3.26)
D. Time averaged fluxes
In this section we will compute the infinite time-
averages 〈Ė〉, 〈L̇z〉 and 〈K̇〉 of the fluxes. These averages
are defined by
〈Ė〉 ≡ lim
∫ T/2
Ė(t)dt. (3.27)
These time-averaged fluxes are sufficient to evolve or-
bits in the adiabatic regime (except for the effect of res-
onances) [12, 25]. In Appendix II, we present two dif-
ferent ways of computing the time averages. The first
approach is based on decoupling the r̃ and θ̃ motion us-
ing the analog of the Mino time parameter for geodesic
motion in Kerr [12]. The second approach uses the ex-
plicit Newtonian parameterization of the orbital motion.
Both averaging methods give the following results:
〈Ė〉 = −32
(1− e2)3/2
e4 − S
cos(ι)
cos(2ι)
cos(2ι)
,(3.28)
〈L̇z〉 = −
(1 − e2)3/2
cos ι
e2 − S
2p3/2 cos ι
+ 7e2 +
cos(2ι)
45 + 148e2 +
cos(2ι)
1 + 3e2 +
, (3.29)
〈K̇〉 = −64
(1 − e2)3/2
e2 − S
2p3/2
+ 37e2 +
cos(ι)
cos(2ι)
cos(2ι)
. (3.30)
Using Eqs. (3.21) and (3.23), we obtain from (3.28) – (3.30) the following time averaged rates of change of the
orbital elements e, p, ι:
〈ṗ〉 = −64
(1− e2)3/2
− S cos(ι)
96p3/2
1064 + 1516e2 + 475e4
149e2
469e2
227e4
cos(2ι)
+ e2 +
[13− cos(2ι)]
, (3.31)
〈ė〉 = −
e(1− e2)3/2
121e2
Se(1− e2)3/2 cos(ι)
5p11/2
1172 + 932e2 +
1313e4
Q(1− e2)3/2
785e2
219e4
+ 13e6 +
2195e2
+ 251e4 +
218e6
cos(2ι)
S2e(1− e2)3/2
2 + 3e2 +
[13− cos(2ι)] , (3.32)
〈ι̇〉 = S sin(ι)(1 − e
2)3/2
p11/2
1− e2
S2 sin(2ι)
240p6
8 + 3e2
8 + e2
Q cot(ι)(1 − e2)3/2
312 + 736e2 − 83e4 −
408 + 1268e2 + 599e4
cos(2ι)
. (3.33)
IV. APPLICATION TO BLACK HOLES
A. Qualitative discussion of results
The above results for the fluxes, Eqs. (3.31), (3.32)
and (3.33) show that the correction terms at O(a2ǫ4)
due to the quadrupole have the same type of effect on the
evolution as the linear spin correction computed by Ryan:
they tend to circularize eccentric orbits and change the
angle ι such as to become antialigned with the symmetry
axis of the quadrupole.
The effects of the terms quadratic in spin are quali-
tatively different. In the expression (3.28) for 〈Ė〉, the
coefficient of cos(2ι) due to the spin self-interaction has
the same sign as the quadrupole term, while the terms
not involving ι have the opposite sign. The terms in-
volving cos(2ι) in Eq. (3.30) for 〈K̇〉 of O(Q) and O(S2)
terms have the same sign, while the terms not involving
ι have the opposite sign. The fractional spin-spin cor-
rection to 〈L̇z〉, Eq. (3.29), has no ι-dependence, and in
expression (3.33) for 〈ι̇〉, the dependence on ι of the two
effects O(Q) and O(S2) is different, too. This is not sur-
prising as the O(Q) effects included here are corrections
to the conservative orbital dynamics, while the effects of
O(S2) that we included are due to radiation reaction.
B. Comparison with previous results
The terms linear in the spin in our results for the time
averaged fluxes, Eqs. (3.28) – (3.33), agree with those
computed by Ryan, Eqs. (14a) – (15c) of [15], and with
those given in Eqs. (2.5) – (2.7) of Ref. [33], when we use
the transformations to the variables used by Ryan given
in Eqs. (2.3) – (2.4) in [33].
Equation (3.28) for the time averaged energy flux
agrees with Eq. (3.10) of Gergely [23] and Eq. (4.15)
of [18] when we use the following transformations:
K = L̄2
Ā2 sin2 κ cos δ − (1− Ā2) cos2 κ
= L̄2
E cos2 κ
(1 + 2L̄2) sin2 κ cos δ
, (4.1)
cos ι = cosκ
E cos2 κ
(1 + 2L̄2) sin2 κ cos δ
, (4.2)
(δ + κ), (4.3)
ξ0 = (ψ0 − ψi) +
, (4.4)
where Ā, L̄, κ, δ, ψ0 and ψi are the quantities used by
Gergely. The first relation here is obtained from the turn-
ing points of the radial motion as follows. We compute
r̃± in terms of E and K and map these expressions back
to r using Eqs. (2.11). The result can then be com-
pared with the turning points in Gergely’s variables, Eq.
(2.19) of [23], using the fact that E is the same in both
cases. Instead of the evolution of the constants of motion
K and Lz, Gergely computes the rates of change of the
magnitude L of the orbital angular momentum and of the
angle κ defined by cosκ = (L · S)/L. Using the trans-
formations (4.1) – (4.4) and the definition of κ we verify
that our Eq. (3.29) agrees with the 〈L̇z〉 computed using
Gergely’s Eqs. (3.23) and (3.35) in [23] and Eq. (4.30)
of [18].
In the limit of the circular equatorial orbits analyzed
by Poisson [22], our Eq. (3.28) agrees with Poisson’s Eq.
(22) when we use the transformations and specializations:
, (4.5)
ι = 0, (4.6)
e2 = 0, (4.7)
cosαA = 1, (4.8)
where v and αA are the variables used by Poisson and the
relation (4.5) is obtained by comparing the expressions
for the constants of motion in the two sets of variables.
The main improvement of our analysis over Gergely’s
is that we express the results in terms of the Carter-type
constant K, which facilitates comparing our results with
other analyses of black hole inspirals. Our computations
also include the spin curvature scattering effects for all
three constants of motion; Gergely [18] only considers
these effects for two of them: the energy and magnitude
of angular momentum, not for the third conserved quan-
tity.
When we expand Eq. (3.28) for small inclination an-
gles and specialize to circular orbits, then after converting
p to the parameter v using Eq. (4.5), we obtain
〈Ė〉 = − 32
11Q− S
= − 32
33− 527
. (4.9)
This result agrees with the terms atO(a2v4) of Eq. (3.13)
of Shibata et al. [24], whose calculations were based on
the fully relativistic expressions. This agreement is a
check that we have taken into account all the contribu-
tions at O(a2ǫ4). The analysis in Ref. [24] could not dis-
tinguish between effects due to the quadrupole and those
due curvature scattering, but we can see from Eq. (4.9)
that those two interactions have the opposite dependence
on ι. Comparing (4.9) with Eq. (3.7) of [24] (which gives
the fluxes into the different modes (l = 2,m, n), where m
and n are the multiples of the ϕ and θ frequencies), we see
that the terms in the (2,±2, 0) and the (2,±1,±1) modes
are entirely due to the quadrupole, while the spin-spin in-
teraction effects are fully contained in the (2,±1, 0) and
(2, 0,±1) modes.
V. NON-EXISTENCE OF A CARTER-TYPE
CONSTANT FOR HIGHER MULTIPOLES
In this section, we show that for a single axisymmetric
multipole interaction, it is not possible to find an ana-
log of the Carter constant (a conserved quantity which
does not correspond to a symmetry of the Lagrangian),
except for the cases of spin (treated by Ryan [15]) and
mass quadrupole moment (treated in this paper). Our
proof is valid only in the approximations in which we
work – expanding to linear order in the mass ratio, to
the leading post-Newtonian order, and to linear order in
the multipole. However we will show below that with
very mild additional smoothness assumptions, our non-
existence result extends to exact geodesic motion in exact
vacuum spacetimes.
We start in Sec. VA by showing that there is no co-
ordinate system in which the Hamilton-Jacobi equation
is separable. Now separability of the Hamilton-Jacobi
equation is a sufficient but not a necessary condition for
the existence of a additional conserved quantity. Hence,
this result does not yield information about the existence
or non-existence of an additional constant. Nevertheless
we find it to be a suggestive result. Our actual derivation
of the non-existence is based on Poisson bracket compu-
tations, and is given in Sec. VB.
A. Separability analysis
Consider a binary of two point masses m1 and m2,
where the mass m1 is endowed with a single axisymmet-
ric current multipole moment Sl or axisymmetric mass
multipole moment Il. In this section, we show that the
Hamilton-Jacobi equation for this motion, to linear order
in the multipoles, to linear order in the mass ratio and to
the leading post-Newtonian order, is separable only for
the cases S1 and I2.
We choose the symmetry axis to be the z-axis and write
the action for a general multipole as
ṙ2 + r2θ̇2 + r2 sin2 θϕ̇2
+ f(r, θ) + g(r, θ)ϕ̇+ E] . (5.1)
For mass moments, g(r, θ) = 0, while for current mo-
ments f(r, θ) = 0. For an axisymmetric multipole of
order l, the functions f and g will be of the form
f(r, θ) =
clIlPl(cos θ)
, g(r, θ) =
dlSl sin θ∂θPl(cos θ)
(5.2)
where Pl(cos θ) are the Legendre polynomials and cl and
dl are constants. We will work to linear order in f and g.
In Eq. (5.1), we have added the energy term needed when
doing a change of time variables, cf. the discussion before
Eq. (2.14) in section III. Since ϕ is a cyclic coordinate,
pϕ = Lz is a constant of motion and the system has
effectively only two degrees of freedom. Note that in the
case of a current moment, there will be correction term
in Lz:
Lz = r
2 sin2 θϕ̇+ g(r, θ). (5.3)
Next, we switch to a different coordinate system
(r̃, θ̃, ϕ) defined by
r = r̃ + α(r̃, θ̃, Lz), (5.4)
θ = θ̃ + β(r̃, θ̃, Lz), (5.5)
where the functions α and β are yet undetermined. We
also define a new time variable t̃ by
1 + γ(r̃, θ̃, Lz)
dt̃. (5.6)
Since we work to linear order in f and g, we can work
to linear order in α, β, and γ. We then compute the
action in the new coordinates and drop the tildes. The
Hamiltonian is given by
p2r(1 + γ − 2α,r) +
(1− 2α
− 2β,θ + γ)
(−α,θ − r2β,r)− E(1 + γ)
2r2 sin2 θ
(1 + γ − 2α
− 2β cot θ)
(1− α
+ γ)− f − gLz
r2 sin2 θ
(5.7)
and the corresponding Hamilton-Jacobi equation is
Ĉ1 +
+ 2V̂ , (5.8)
where we have denoted
Ĉ1 = J(r, θ) [1 + γ − 2α,r] = 1 + γ − 2α,r + j, (5.9)
Ĉ2 = J(r, θ)
1− 2α
− 2β,θ + γ
= 1− 2α
− 2β,θ + γ + j, (5.10)
Ĉ3 = J(r, θ)
−α,θ − r2β,r
= −α,θ − r2β,r, (5.11)
V̂ = J(r, θ)
2r2 sin2 θ
(1 + γ − 2α
− 2β cot θ)
+ γ)− E(1 + γ)
− f − gLz
r2 sin2 θ
2r2 sin2 θ
(1 + γ − 2α
− 2β cot θ + j)
−E(1 + γ + j)− 1
(1− α
+ γ + j)
−f − gLz
r2 sin2 θ
. (5.12)
The unperturbed problem is separable, so make the
perturbed problem separable, we have multiplied the
Hamilton-Jacobi equation by an arbitrary function
J(r, θ), which can be expanded as J(r, θ) = 1 + j(r, θ),
where j(r, θ) is a small perturbation.
To find a solution of the form W =Wr(r)+Wθ(θ), we
first specialize to the case where Ĉ3 = 0:
− Ĉ3 = β,rr2 + α,θ = 0. (5.13)
We differentiate Eq. (5.8) with respect to θ, using Eq.
(5.8) to write (dWr/dr)
2 in terms of (dWθ/dθ)
2 and then
differentiate the result with respect to r to obtain
∂θĈ2
∂θĈ1
2V̂ ∂θĈ1
Ĉ1Ĉ2
. (5.14)
Expanding Eq. (5.14) to linear order in the small quan-
tities then yields the two conditions for the kinetic and
the potential part of the Hamiltonian to be separable:
0 = ∂r∂θ
2α,r −
− 2β,θ
, (5.15)
sin2 θ
2β,r cot
2 θ − 3β,rθ cot θ + β,r csc2 θ
sin2 θ
+ α,rθ
−∂r∂θ
Pl(cos θ) +
dlSlLz
rl sin θ
∂θPl(cos θ)
2α,rθ −
+ 2Er2α,rθ
, (5.16)
where we have used Eq. (5.2) for f and g. Therefore, the
following conditions must be satisfied:
M4(θ)−N(r) =
+ β,θ − 2α,r, (5.17)
M1(θ) = 2β cot
2 θ + β csc2 θ + β,θθ
−3β,θ cot θ, (5.18)
M2(θ) = r
2∂r(r
2β,r), (5.19)
M3(θ) = 2rα,rθ − α,θ +
∂θPl(cos θ)
−SlLz
∂θ(csc θ ∂θPl(cos θ)). (5.20)
Here, the functions M and N are arbitrary integration
constants.
Solving the condition for the kinetic term to be sep-
arable, Eq. (5.17), together with Eq. (5.13) gives the
general solution that goes to zero at large r as
cos(nθ + ν), (5.21)
β = − A
sin(nθ + ν), (5.22)
where A and ν are arbitrary and n is an integer. These
functions must satisfy the conditions (5.18) – (5.20) in
order for the potential term to be separable as well. To
see when this will be the case, we start by considering Eq.
(5.20). Substituting the general ansatz α = a1(r)a2(θ)
shows that a′2 = P
l or a
2 = (cscθ P
′ depending on
whether a mass or a current multipole is present. The
function a1(r) is then determined from
0 = 2ra′1 − a1 +
clIl/r
(l−1)
dlSlLz/r
(5.23)
Hence,
[clIl/(2l)] r
(1−l)
[dlSlLz/(2l+ 1)] r
(5.24)
so that we obtain for mass moments
Pl(cos θ)
, β =
P ′l (cos θ)
(5.25)
and for current moments
dlSlLz
2l+ 1
csc θP ′l (cos θ)
, (5.26)
dlSlLz
(2l+ 1)(l + 1)
(csc θ P ′l (cos θ))
, (5.27)
where we have used the condition (5.13) to solve for β.
Substituting this in Eq. (5.19) determines that l = 2
for mass moments and l+1 = 2 for current moments. For
an l = 2 mass moment, conditions (5.17) and (5.18) are
satisfied as well, with n = 2 and ν = 0. For the case of an
l = 1 current moment, the extra term inH is independent
of θ anyway. But for any other multipole interaction,
the Hamilton-Jacobi equation will not be separable. For
example, for the current octupole Sijk, the last term in
Eq. (5.7) is proportional to S3Lz(5 cos
2 θ − 1)/r5 and
is therefore not separable. From Eq. (5.2) one can see
that, for a general multipole, the functions f or g contain
different powers of cos θ appearing with the same power
of r since the Legendre polynomials can be expanded as
[34]:
Pl(cos θ) =
(−1)n(2l− 2n)!
2ln!(l − n)!(l − 2n)!
(cos θ)l−2n, (5.28)
where N = l/2 for even l and N = (l + 1)/2 for odd l.
It will not be possible to cancel all of these terms with
(5.21) – (5.22) for l > 2.
The case when Ĉ3 is non-vanishing will only be sepa-
rable if all the coefficients are functions of r or of θ only,
and if in addition, the potential also depends only on r or
on θ. Achieving this for our problem will not be possible
because the potential cannot be transformed to the form
required for separability.
B. Derivation of non-existence of additional
constants of the motion
In this subsection, we show using Poisson brackets that
for a single axisymmetric multipole interaction, to linear
order in the multipole and the mass ratio, a first integral
analogous to the Carter constant does not exist, except
for the cases of mass quadrupole and spin.
Suppose that such a constant does exist. We write the
Hamiltonian corresponding to the action (5.1) as H =
H0 + δH and the Carter-type constant as K = K0 +
δK(pr, pθ, Lz, r, θ), where
2r2 sin2 θ
, (5.29)
δH = − clIl
Pl(cos θ)−
dlSlLz
rl+2 sin θ
∂θPl(cos θ),(5.30)
K0 = p
sin2 θ
. (5.31)
Computing the Poisson bracket gives, to linear order in
the perturbations
0 = {H0, δK}+ {δH,K0} (5.32a)
δK + {δH,K0}, (5.32b)
where we have used that {H0,K0} = 0 and the fact that
{H0, δK} = d(δK)/dt. Here, d/dt denotes the total time
derivative along an orbit (r(t), θ(t), pr(t), pθ(t)) of H0 in
phase space. The partial differential equation (5.32a) for
δK thus reduces to a set of ordinary differential equa-
tions that can be integrated along the individual orbits
in phase space.
The unperturbed motion for a bound orbit is in a
plane, so we can switch from spherical to plane polar co-
ordinates (r, ψ). In terms of these coordinates, we have
H0 = p
r/2+p
ψ/2, K0 = p
ψ, and cos θ = sin ι sin(ψ+ψ0),
with cos ι = Lz/
K and the constant ψ0 denoting the
angle between the direction of the periastron and the
intersection between the orbital and equatorial plane.
Then Eq. (5.32) becomes
δK = η(t), (5.33)
η(t) = − 2pψ dlSlLz
sin ι rl+2(t)
∂ψPl(sin ι sin(ψ(t) + ψ0))
cos(ψ(t) + ψ0)
2pψ clIl
rl+1(t)
∂ψPl(sin ι sin(ψ(t) + ψ0)). (5.34)
For unbound orbits, one can always integrate Eq.
(5.33) to determine δK. However, for bound periodic
orbits there is a possible obstruction: the solution for
the conserved quantity K0 + δK will be single valued if
and only if the integral of the source over the closed orbit
vanishes,
∮ Torb
η(t)dt = 0. (5.35)
Here, Torb is the orbital period. In other words, the par-
tial differential equation (5.32) has a solution δK if and
only if the condition (5.35) is satisfied. This is the same
condition as obtained by the Poincare-Mel’nikov-Arnold
method, a technique for showing the non-integrability
and existence of chaos in certain classes of perturbed dy-
namical systems [35].
Thus, it suffices to show that the condition (5.35) is
violated for all multipoles other than the spin and mass
quadrupole. To perform the integral in Eq. (5.35), we use
the parameterization for the unperturbed motion, r =
K/(1+ e cosψ) and dt/dψ = K3/2/(1+ e cosψ)2, so that
the condition for the existence of a conserved quantity
K0 + δK becomes
clIl(1 + e cosψ)
l−1∂ψPl(sin ι sin(ψ + ψ0))−
dlSlLz
K sin ι
(1 + e cosψ)l∂ψ
∂ψPl(sin ι sin(ψ + ψ0))
cos(ψ + ψ0)
(5.36)
In terms of the variable χ = ψ + ψ0 − π/2, Eq. (5.36) can be written as
dχclIl [1 + e(sinψ0 cosχ− cosψ0 sinχ)]l−1
Pl(sin ι cosχ)
dlSlLz
sin ι
[1 + e(sinψ0 cosχ− cosψ0 sinχ)]l
Pl(sin ι cosχ)
. (5.37)
Inserting the expansion (5.28) for Pl(cosχ), taking the derivatives, and using the binomial expansion for the first term
in Eq. (5.37), we get
0 = clIl
Alnjk e
j(sin ι)l−2n(sinψ0)
k(cosψ0)
dχ (sinχ)j−k+1(cosχ)k+l−2n−1
dlSlLz
Blnjk e
j(sin ι)l−2n−1(sinψ0)
k(cosψ0)
dχ (sinχ)j−k+1(cosχ)k+l−2n−2. (5.38)
The coefficients Alnkj and Blnkj are
Alnkj =
(−1)n+k+1(l − 1)!(2l − 2n)!
2ln!(l − 1− j)!k!(j − k)!(l − n)!(l − 2n− 1)!
, Blnkj =
(−1)n+kl!(2l− 2n)!
2ln!(l − j)!k!(j − k)!(l − n)!(l − 2n− 2)!
. (5.39)
The only non-vanishing contribution to the integrals in Eq. (5.38) will come from terms with even powers of both
cosχ and sinχ. These can be evaluated as multiples of the beta function:
0 = clIl
Clnjk e
j(sin ι)l−2n(sinψ0)
k(cosψ0)
j−k δ(j−k+1),even δ(l+k−1),even
dlSlLz
Dlnjk e
j(sin ι)l−2n−1(sinψ0)
k(cosψ0)
j−k δ(j−k+1),even δ(l+k),even. (5.40)
Here, the coefficients are
Clnjk =
2Γ( j
+ 1)Γ(k
− n+ 1)
Alnkj , Dlnjk =
2Γ( j
+ 1)Γ(k
− n− 1
− n+ 3
Blnkj (5.41)
Eq. (5.40) shows that for even l, terms with j =even
(odd) and k =odd (even) give a non-vanishing contribu-
tion for the case of a mass (current) multipole, and hence
K0+δK is not a conserved quantity for the perturbed mo-
tion. Note that terms with j =even and k =odd for even
l occur only for l > 3, so for l = 2 the mass quadrupole
term in Eq. (5.40) vanishes and therefore there exists an
analog of the Carter constant, which is consistent with
our results of Sec. II and our separability analysis. For
odd l, terms with j =odd (even) and k =even (odd) are
finite for Il (Sl). Note that for the case l = 1 of the spin,
the derivatives with respect to χ in Eq. (5.37) evaluate to
zero, so in this case there also exists a Carter-type con-
stant. These results show that for a general multipole
other than I2 and S1, there will not be a Carter-type
constant for such a system.
1. Exact vacuum spacetimes
Our result on the non-existence of a Carter-type con-
stant can be extended, with mild smoothness assump-
tions, to falsify the conjecture that all exact, axisymmet-
ric vacuum spacetimes posess a third constant of the mo-
tion for geodesic motion. Specifically, we fix a multipole
order l, and we assume:
• There exists a one parameter family
(M, gab(λ))
of spacetimes, which is smooth in the parameter λ,
such that λ = 0 is Schwarzschild, and each space-
time gab(λ) is stationary and axisymmetric with
commuting Killing fields ∂/∂t and ∂/∂φ, and such
that all the mass and current multipole moments of
the spacetime vanish except for the one of order l.
On physical grounds, one expects a one parameter
family of metrics with these properties to exist.
• We denote by H(λ) the Hamiltonian on the tan-
gent bundle overM for geodesic motion in the met-
ric gab(λ). By hypothesis, there exists for each λ
a conserved quantity M(λ) which is functionally
independent of the conserved energy and angular
momentum. Our second assumption is that M(λ)
is differentiable in λ at λ = 0. One would expect
this to be true on physical grounds.
• We assume that the conserved quantity M(λ) is
invariant under the symmetries of the system:
L~ξM(λ) = L~ηM(λ) = 0,
where ~ξ and ~η are the natural extensions to the 8
dimensional phase space of the Killing vectors ∂/∂t
and ∂/∂φ. This is a very natural assumption.
These assumptions, when combined with our result of
the previous section, lead to a contradiction, showing
that the conjecture is false under our assumptions.
To prove this, we start by noting that M(0) is a con-
served quantity for geodesic motion in Schwarzschild, so
it must be possible to express it as some function f of
the three independent conserved quantities:
M(0) = f(E,Lz,K0). (5.42)
Here E is the energy, Lz is the angular momentum, and
K0 is the Carter constant. Differentiating the exact re-
lation {H(λ),M(λ)} = 0 and evaluating at λ = 0 gives
{H0,M1} =
{E,H1}+
{Lz, H1}+
{K0, H1},
(5.43)
where H0 = H(0), H1 = H
′(0), and M1 = M
′(0). As
before, we can regard this is a partial differential equa-
tion that determines M1, and a necessary condition for
solutions to exist and be single valued is that the integral
of the right hand side over any closed orbit must vanish:
{E,H1}+
{Lz, H1}+
{K0, H1}
(5.44)
Now strictly speaking, there are no closed orbits in
the eight dimensional phase space. However, the ar-
gument of the previous section applies to orbits which
are closed in the four dimensional space with coordinates
(r, θ, pr, pθ), since by the third assumption above every-
thing is independent of t and φ, and pt and pφ are con-
served. Here (t, r, θ, φ) are Schwarzschild coordinates and
(pt, pr, pθ, pφ) are the corresponding conjugate momenta.
Next, we can pull the partial derivatives ∂f/∂E etc.
outside of the integral. It is then easy to see that the first
two terms vanish, since there do exist a conserved energy
and a conserved z-component of angular momentum for
the perturbed system. Thus, Eq. (5.44) reduces to
{K0, H1} = 0. (5.45)
Since M(0) is functionally independent of E and Lz, the
prefactor ∂f/∂K0 must be nonzero, so we obtain
{K0, H1} = 0. (5.46)
The result (5.46) applies to fully relativistic orbits in
Schwarzschild. We need to take the Newtonian limit of
this result in order to use the result we derived in the
previous section. However, the Newtonian limit is a lit-
tle subtle since Newtonian orbits are closed and generic
relativistic orbits are not closed. We now discuss how the
limit is taken.
The integral (5.46) is taken over any closed orbit in
the four dimensional phase space (r, θ, pr, pθ) which cor-
responds to a geodesic in Schwarzschild. Such orbits are
non generic; they are the orbits for which the ratio be-
tween the radial and angular frequencies ωr and ωθ is a
rational number. We denote by qr and qθ the angle vari-
ables corresponding to the r and θ motions [36]. These
variables evolve with proper time τ according to
qr = qr,0 + ωrτ, (5.47a)
qθ = qθ,0 + ωθτ, (5.47b)
where qr,0 and qθ,0 are the initial values. We denote the
integrand in Eq. (5.46) by
I(qr, qθ, a, ε, ι),
where I is some function, and a, ε and ι are the parame-
ters of the geodesic defined by Hughes [32] (functions of
E, Lz and K0). The result (5.46) can be written as
∫ T/2
dτ I[qr(τ), qθ(τ), a, ε, ι] = 0, (5.48)
where T = T (a, ε, ι) is the period of the r, θ motion.
Since the variables qr and qθ are periodic with period
2π, we can express the function I as a Fourier series
I(qr, qθ, a, ε, ι) =
n,m=−∞
Inm(a, ε, ι)einqr+imqθ . (5.49)
Now combining Eqs. (5.47), (5.48) and (5.49) gives
n,m=−∞
Inm(a, ε, ι)einqr,0+imqθ,0
×Si [(nωr +mωθ)T/2] , (5.50)
where Si(x) = sin(x)/x. Since the initial conditions qr,0
and qθ,0 are arbitrary, it follows that
Inm(a, ε, ι)Si [(nωr +mωθ)T/2] = 0 (5.51)
for all n, m.
Next, for closed orbits the ratio of the frequencies must
be a rational number, so
, (5.52)
where p and q are integers with no factor in common.
These integers depend on a, ε and ι. The period T is
given by 2π/T = qωr = pωθ. The second factor in Eq.
(5.51) now simplifies to
(np+mq)π
, (5.53)
which vanishes if and only if
n = n̄q, m = m̄p, n̄+ m̄ 6= 0, (5.54)
for integers n̄, m̄. It follows that
Inm(a, ε, ι) = 0 (5.55)
for all n, m except for values of n, m which satisfy the
condition (5.54)
Consider now the Newtonian limit, which is the limit
a → ∞ while keeping fixed ε and ι and the mass of the
black hole. We denote by IN(qr , qθ, a, ε, ι) the Newtonian
limit of the function I(qr , qθ, a, ε, ι). The integral (5.48)
in the Newtonian limit is given by the above computation
with p = q = 1, since ωr = ωθ in this limit. This gives
dτIN =
INn,−n(a, ε, ι) ein(qr,0−qθ,0), (5.56)
where INnm are the Fourier components of IN. In the
previous subsection, we showed that this function is non-
zero, which implies that there exists a value k of n for
which IN k,−k 6= 0.
Now as a → ∞, we have ωr/ωθ → 1, and hence from
Eq. (5.52) there exists a critical value ac of a such that
the values of p and q exceed k for all closed orbits with
a > ac. (We are keeping fixed the values of ε and ι). It
follows from Eqs. (5.54) and (5.55) that
Ik,−k(a, ε, ι)
IN k,−k(a, ε, ι)
= 0 (5.57)
for all such values of a. However this contradicts the fact
Ik,−k(a, ε, ι)
IN k,−k(a, ε, ι)
→ 1 (5.58)
as a→ ∞. This completes the proof.
Hence, if the three assumptions listed at the start of
this subsection are satisfied, then the conjecture that all
vacuum, axisymmetric spacetimes possess a third con-
stant of the motion is false.
Finally, it is sometimes claimed in the classical dynam-
ics literature that perturbation theory is not a sufficiently
powerful tool to assess whether the integrability of a sys-
tem is preserved under deformations. An example that
is often quoted is the Toda lattice Hamiltonian [38, 39].
This system is integrable and admits a full set of con-
stants of motion in involution. However, if one approx-
imates the Hamiltonian by Taylor expanding the poten-
tial about the origin to third order, one obtains a sys-
tem which is not integrable. This would seem to indicate
that perturbation theory can indicate a non-integrability,
while the exact system is still integrable.
In fact, the Toda lattice example does not invalidate
the method of proof we use here. If we write the Toda
lattice Hamiltonian as H(q,p), then the situation is that
H(λq,p) is integrable for λ = 1, but it is not integrable
for 0 < λ < 1. Expanding H(λq,p) to third order in λ
gives a non-integrable Hamiltonian. Thus, the perturba-
tive result is not in disagreement with the exact result
for 0 < λ < 1, it only disagrees with the exact result for
λ = 1. In other words, the example shows that pertur-
bation theory can fail to yield the correct result for finite
values of λ, but there is no indication that it fails in ar-
bitrarily small neighborhoods of λ = 0. Our application
is qualitatively different from the Toda lattice example
since we have a one parameter family of Hamiltonians
H(λ) which by assumption are integrable for all values
of λ.
VI. CONCLUSION
We have examined the effect of an axisymmetric
quadrupole moment Q of a central body on test parti-
cle inspirals, to linear order in Q, to the leading post-
Newtonian order, and to linear order in the mass ratio.
Our analysis shows that a natural generalization of the
Carter constant can be defined for the quadrupole inter-
action. We have also analyzed the leading order spin self-
interaction effect due to the scattering of the radiation off
the spacetime curvature due to the spin. Combining the
effects of the quadrupole and the leading order effects
linear and quadratic in the spin, we have obtained ex-
pressions for the instantaneous as well as time-averaged
evolution of the constants of motion for generic orbits un-
der gravitational radiation reaction, complete at O(a2ǫ4).
We have also shown that for a single multipole interaction
other than Q or spin, in our approximations, a Carter-
type constant does not exist. With mild additional as-
sumptions, this result can be extended to exact space-
times and falsifies the conjecture that all axisymmetric
vacuum spacetimes possess a third constant of motion for
geodesic motion.
VII. ACKNOWLEDGMENTS
This research was partially supported by NSF grant
PHY-0457200. We thank Jeandrew Brink for useful cor-
respondence.
Appendix A: Time variation of quadrupole: order of
magnitude estimates
In this appendix, we give an estimate of the timescale
Tevol for the quadrupole to change. The analysis in the
body of this paper is valid only when Tevol ≫ Trr, where
Trr is the radiation reaction time, since we have neglected
the time evolution of the quadrupole. We distinguish be-
tween two cases: (i) when the central body is exactly non-
spinning but has a quadrupole, and (ii) when the central
body has finite spin in addition to the quadrupole.
1. Estimate of the scaling for the nonspinning case
For the purpose of a crude estimate, the relevant in-
teraction is the tidal interaction with energy
QijEij ∼ −
Q̄I cos2 θ, (A1)
where Eij is the tidal field, θ is the angle between the
symmetry axis and the normal to the orbital plane of
m2, and we have written the quadrupole as Q ∼ Q̄I,
where Q̄ is dimensionless and I is the moment of inertia.
For small deviations from equilibrium, the relevant piece
of the Lagrangian is schematically
L ∼ Iψ̇2 + Q̄I m2
ψ2. (A2)
We define the evolution timescale Tevol to be the time
it takes for the angle to change by an amount of order
unity, and since the amplitude of the oscillation scales
roughly as ∼ m2/m1, the evolution time scales as
T−2evol ∼
ω2orbit, (A3)
where ω2orbit = M/r
3. Thus, the ratio of the evolution
timescale compared to the radiation reaction timescale
scales as
Tevol/Trr ∼
. (A4)
2. Estimate of the scaling for the spinning case
When the body is spinning the effect of the tidal cou-
pling is to cause a precession. For the purpose of this
estimate, we calculate the torque on m1 due to the com-
panion’s Newtonian field. The torque N scales as
Ni ∼ ǫimjQmkEjk. (A5)
We assume that the precession is slow, i.e.
ωprec ≪ S̄/m1
, (A6)
where ωprec is the precession frequency and S̄ = S/m
is the dimensionless spin. This gives the approximate
scaling of the precession timescale as (cf. [37])
Tprec/Trr ∼
. (A7)
and the evolution timescale is thus
Tevol/Trr ∼
. (A8)
Because of our assumption (A6) that the precession is
slow, equation (A8) is valid only when
) S̄2
. (A9)
When S̄ is sufficiently small that the condition (A9) is
violated, the relevant timescale is instead given by Eq.
(A3).
3. Application to Kerr inspirals
For Kerr inspirals,
S̄ ∼ a, Q̄ ∼ a2, µ/M ≪ 1 and r ∼M. (A10)
Therefore, the condition (A9) is satisfied, and the pre-
cession time is longer than the radiation reaction time
Tprec/Trr ∼
. (A11)
Note that for Kerr inspirals, since r ∼ M both formulas
(A3) and (A7) give the same scaling.
Moreover, for Kerr inspirals, the amplitude of the pre-
cession will be small, of order the mass ratio µ/M . This is
because of angular momentum conservation: in the rela-
tivistic regime, the orbital angular momentum is a factor
of µ/M smaller than the angular momentum of the black
hole and can therefore not cause a large precession ampli-
tude. Even if the orbital angular momentum at infinity
is large, most of it will be radiated away as outgoing
gravitational waves during the earlier phase of the inspi-
ral. This factor of µ/M is taken into account when we
consider the evolution timescale, which for Kerr inspirals
reduces to
Tevol/Trr ∼
. (A12)
Since 1/a ≥ 1, M/r ∼ 1 and M/µ ≪ 1, the evolution
time is long compared to the radiation reaction time and
we can neglect the time variation of the quadrupole at
leading order.
Appendix B: Computation of time averaged fluxes
1. Averaging method that parallels fully
relativistic averaging
We start by noting that the differential equations
(2.26) and (2.27) governing the r̃ and θ̃ motions decouple
if we define a new time parameter t̂ by
dt̂ =
dt̃. (B1)
This is the analog of the Mino time parameter for
geodesic motion in Kerr [12]. The equations of motion
(2.26)–(2.24) then become
= V̂r̃(r̃), (B2)
V̂r̃(r̃) = 2Er̃
4 + 2r̃3 −Kr̃2 − 4SLzr̃
r̃ − 2L2z
, (B3)
= V̂θ̃(θ̃), (B4)
V̂θ̃(θ̃) = K −
sin2 θ̃
−QE cos 2θ̃, (B5)
= V̂ϕr̃(r̃) + V̂ϕθ̃(θ̃), (B6)
V̂ϕr̃(r̃) =
, V̂ϕθ̃(θ̃) =
sin2 θ̃
. (B7)
The parameters t and t̂ are related by:
= V̂tr̃(r̃) + V̂tθ̃(θ̃) (B8)
V̂tr̃(r̃) = r̃
2, V̂tθ̃(θ̃) =
cos 2θ̃. (B9)
It follows from Eqs. (B2) and (B4) that the functions
r̃(t̂) and θ̃(t̂) are periodic; and we denote their periods
by Λr̃ and Λθ̃. We define the fiducial motion associated
with the constants of motion E, Lz and K to be the
motion with the initial conditions r̃(0) = r̃min and θ̃(0) =
θ̃min, where r̃min and θ̃min are given by the vanishing of
the right-hand sides of Eqs. (B2) and (B4) respectively.
The functions r̂(t̂) and θ̂(t̂) associated with this fiducial
motion are given by
∫ r̂(t̂)
r̃min
V̂r̃(r̃)
= t̂, (B10)
∫ θ̂(t̂)
θ̃min
V̂θ̃(θ̃)
= t̂. (B11)
From Eq. (B8) it follows that
t(t̂) = t0 +
V̂tr̃[r̃(t
′)] + V̂tθ̃[θ̃(t
, (B12)
where t0 = t(0). Next, we define the constant Γ to be
the following average value:
∫ Λr̃
dt′V̂tr̃[r̂(t
′)] +
dt′V̂tθ̃[θ̂(t
′)]. (B13)
Then we can write t(t̂) as a sum of a linear term and
terms that are periodic:
t(t) = t0 + Γt̂+ δt(t̂), (B14)
where δt(t̂) denotes the oscillatory terms in Eq. (B12).
To average a function over the time parameter t̂, it is
convenient to parameterize r̃ and θ̃ in terms of angular
variables as follows. For the average over θ̃ we introduce
the parameter χ by
cos2 θ̂(t̂) = z− cos
2 χ, (B15)
where z− = cos
2 θ̃− with z− being the smaller root of Eq.
(B4):
K + 3QE ±
(K −QE)2 + 4QEL2z
(B16)
and where β = 2QE. Then from the definition (B11)
of θ̂ together with Eq. (B4) and the requirement that χ
increases monotonically with t̂ we obtain
β (z+ − z− cos2 χ). (B17)
Then we can write the average over t̂ of a function Fθ̃(t̂)
which is periodic with period Λθ̃ in terms of χ as
〈Fθ̃〉t̂ =
dt̂Fθ̃(t̂)
Fθ̃[t̂(χ)]
β (z+ − z− cos2 χ)
, (B18)
where
Λθ̃ =
β (z+ − z− cos2 χ)
. (B19)
Similarly, to average a function Fr̃(t̂) that is periodic with
period Λr̃, we introduce a parameter ξ via
1 + e cos ξ
, (B20)
where the parameter ξ varies from 0 to 2π as r̃ goes
through a complete cycle. Then,
= P (ξ), (B21)
P (ξ) ≡
V̂r̃[r̃(ξ)]
pe | sin ξ |
(1 + e cos ξ)
(B22)
The average over t̂ of Fr̃(t̂) can then be computed from
〈Fr̃〉t̂ =
dξ Fr̃/P (ξ)
dξ/P (ξ)
. (B23)
Now, a generic function Fr̃,θ̃[r̃(t̂), θ̃(t̂)] will be biperiodic
in t̂: Fr̃,θ̃[r̃(t̂+Λr̃), θ̃(t̂+Λθ̃)] = Fr̃,θ̃[r̃(t̂), θ̃(t̂)]. Combin-
ing the results (B18) and (B23) we can write its average
as a double integral over χ and ξ as
〈Fr̃,θ̃〉t̂ =
Λθ̃Λr̃
Fr̃,θ̃[r̃(ξ), θ̃(χ)]
β (z+ − z− cos2 χ)P (ξ)
(B24)
To compute the time average of Ė, L̇z, and K̇, we need
to convert the average of a function over t̂ calculated from
(B24) to the average over t. As explained in detail in
[9], in the adiabatic limit we can choose a time interval
∆t which is long compared to the orbital timescale but
short compared to the radiation reaction time. From
Eq. (B12) we have ∆t = Γt̂+ osc.terms. The oscillatory
terms will be bounded and will therefore be negligible in
the adiabatic limit, so we have to a good approximation
〈Ė〉t =
〈Ė V̂t〉t̂, (B25)
where V̂t ≡ V̂tr̃ + V̂tθ̃, cf. Eq. (B8), and similarly for L̇z
and K̇.
The explicit results we obtain using this method are
given in section III, Eqs. (3.28), (3.29), and (3.30).
2. Averaging method using the explicit
parameterization of Newtonian orbits
To perform the time-averaging using this method, we
define a parameter ξ via
1 + e cos ξ
, (B26)
where the parameter ξ varies from 0 to 2π as r̃ goes
through a complete cycle. Note that θ appears in Eqs.
(3.16) – (3.18) only in terms that are linear in Q, so we
can write θ in terms of ξ using the Newtonian relation
x3 = r cos θ = r sin ι sin(ξ + ξ0). (B27)
Here, ξ0 is the angle between the direction of the peri-
helion and the intersection of the orbital and equatorial
plane. Similarly, for the ṙθ̇ terms in Eqs. (3.17) and
(3.26) we can use the Newtonian relations ṙ = e/
p sin ξ
and ξ̇ =
p/r2. From Eqs. (2.27) and (B20) it follows
(1 + e cos ξ)2
−3 + e2 − 2e cos ξ + 2 cos2 ι(8 − e2 + 8e cos ξ + e2 cos 2ξ)
, (B28)
and from Eq. (2.12)
(1 + e cos ξ)
2 sin2 ι sin2(ξ + ξ0)− 1
. (B29)
Using these expressions, we compute the time-averaged fluxes from
〈Ė〉 =
dξ Ė (dt/dt̃) (dt̃/dξ)
dξ (dt/dt̃) (dt̃/dξ)
(B30)
and obtain:
〈Ė〉 = −32
(1− e2)3/2
e4 − S
cos(ι)
cos(2ι)
cos(2ι)
cos(2ξ0) sin
cos(2ξ0) sin
, (B31)
〈L̇z〉 = −
(1− e2)3/2
cos ι
e2 − S
2p3/2 cos ι
+ 7e2 +
cos(2ι)
−3− 45
45 + 148e2 +
cos(2ι)
1 + 3e2 +
e2 cos(2ξ0) sin
, (B32)
〈K̇〉 = −64
(1 − e2)3/2
e2 − S
2p3/2
+ 37e2 +
cos(ι)
cos(2ι)
cos(2ι)
e2 cos(2ξ0) sin
. (B33)
In the adiabatic limit, the terms involving cos(2ξ0) can
be omitted because they average to zero. As explained
by Ryan [15], the radiation reaction timescale for terms
involving ξ0 is much longer than the precession timescale
for most orbits, so the terms involving ξ0 will average
away. This is consistent with our results for the adia-
batic infinite time-averaged fluxes using the Mino time
parameter. The Mino-time averaging method was based
on the assumption that the fundamental frequencies are
incommensurate and the motion fills up the whole torus,
which is equivalent to averaging over ξ0.
[1] L. Barack and C. Cutler, Phys. Rev. D 69, 082005 (2004) [2] K. Glampedakis and S. Babak, Class. Quantum Grav.
23, 4167 (2006)
[3] D. A. Brown, et al. gr-qc/0612060
[4] E. Poisson, Living Rev. Relat. 7, 6 (2004),
http://relativity.livingreviews.org/Articles/lrr-2004-6/index.html
[5] Special Issue: Gravit. Rad. from Binary Black Holes: Ad-
vances in the perturbative approach, Class. Quant. Grav.
22 (2005)
[6] K. Glampedakis, Class. Quantum Grav. 22, S605 (2005)
[7] S. Drasco, Class. Quantum Grav. 23, S769 (2006)
[8] M. Favata and É. É. Flanagan Accuracy of adiabatic
waveforms for eccentric orbits, (in preparation)
[9] S. Drasco, É. É. Flanagan, and S. A. Hughes, Class.
Quantum Grav. 22, S801 (2005)
[10] É. É. Flanagan and T. Hinderer, Two timescale analysis
of extreme mass ratio inspirals in Kerr. II. Numerical
integration through resonances, (in preparation)
[11] N. Sago, et al., Progr. Theor. Phys. 114, 509 (2005)
[12] Y. Mino, Phys. Rev. D 67, 084027 (2003)
[13] S. A. Hughes, S. Drasco, É. É. Flanagan, and J. Franklin,
Phys. Rev. Lett. 94, 221101 (2005)
[14] F. D. Ryan, Phys. Rev. D 52, R3159 (1995)
[15] F. D. Ryan, Phys. Rev. D 53, 3064 (1996)
[16] L. Blanchet, T. Damour, G. Farese, and B. Iyer, Phys.
Rev. D 71, 124004 (2005)
[17] L. Kidder, C. Will, and A. Wiseman, Phys. Rev. D 47,
R4183 (1993)
[18] L. A. Gergely, Phys. Rev. D 61, 024035 (1999)
[19] C. Will, Phys. Rev. D 71, 084027 (2005)
[20] G. Faye, L. Blanchet, and A. Buonanno, Phys. Rev. D
74, 104033 (2006)
[21] L. Blanchet, A. Buonanno, and G. Faye, Phys. Rev. D
74, 104034 (2006)
[22] E. Poisson, Phys. Rev. D 57, 5287 (1998)
[23] L. A. Gergely and Z. Keresztes, Phys. Rev. D 67, 024020
(2003)
[24] M. Shibata, M. Sasaki, H. Tagoshi, and T. Tanaka, Phys.
Rev. D 51, 1646 (1995)
[25] É. É. Flanagan and T. Hinderer, Two timescale analy-
sis of extreme mass ratio inspirals in Kerr. I. General
formalism, (in preparation)
[26] C. Misner, K. Thorne, and J. Wheeler, Gravitation (W.H.
Freeman and Co., San Francisco, 1973)
[27] R. O. Hansen, J. Math. Phys. 15, 46 (1974)
[28] K. S. Thorne, Rev. Mod. Phys 52, 300 (1980)
[29] T. Bäckdahl, gr-qc/0612043 (2006); Class. Quantum
Grav. 22, 3585 (2005)
[30] C. Li and G. Lovelace, gr-qc/0702146 (2007)
[31] L. Blanchet and T. Damour, Phys. Lett. 104A, 82 (1984)
[32] S. A. Hughes, Phys. Rev. D 61, 084004 (2000)
[33] K. Glampedakis, S. Hughes, and D. Kennefick, Phys.
Rev. D 66, 064005 (2002)
[34] G. Arfken, Mathematical Methods for Physicists (Aca-
demic Press, CA, 1985)
[35] V. K. Melnikov, Trans. Moscow Math. Soc. 12, 1-56
(1956)
[36] W. Schmidt, Celestial mechanics in Kerr spacetime
Class. Quantum Grav. 19 (2002) 2743-2764
[37] H. Goldstein, C. Poole, and J. Safko, Classical Mechanics
(Addison-Wesley, 2002)
[38] Ý. Birol and A. Hacinliyan, Phys. Rev. E 52, 4750 (1995)
[39] H. Yoshida, Comm. Math. Phys. 116, 529 (1988)
http://arxiv.org/abs/gr-qc/0612060
http://relativity.livingreviews.org/Articles/lrr-2004-6/index.html
http://arxiv.org/abs/gr-qc/0612043
http://arxiv.org/abs/gr-qc/0702146
