a

r
X
i
v
:
1
1
1
1
.
3
6
0
5
v
1


[
g
r
-
q
c
]


1
5

N
o
v

2
0
1
1
Resonances in Extreme Mass-Ratio Inspirals:
Asymptotic and Hyperasymptotic Analysis
Jonathan Gair,
1
Nicol as Yunes,
2, 3
and Carl M. Bender
4
1
Institute of Astronomy, Madingley Road, Cambridge, CB30HA, United Kingdom
2
Department of Physics, Montana State University, Bozeman, Montana 59717, USA
3
MIT and Kavli Institute, Cambridge, Massachusetts 02139, USA
4
Department of Physics, Washington University, St. Louis, Missouri 63130, USA
(Dated: November 16, 2011)
An expected source of gravitational waves for future detectors in space are the inspirals of small
compact objects into much more massive black holes. These sources have the potential to provide a
wealth of information about astronomy and fundamental physics. On short timescales the orbit of
the small object is approximately geodesic. Generic geodesics for a Kerr black hole spacetime have
a complete set of integrals and can be characterized by three frequencies of the motion. Over the
course of an inspiral, a typical system will pass through resonances where two of these frequencies
become commensurate. The eect of the resonance will be to alter signicantly the rate of inspiral for
the duration of the resonance. Understanding the impact of these resonances on gravitational wave
phasing is important to detect and exploit these signals for astrophysics and fundamental physics.
Two dierential equations that might describe the passage of an inspiral through such a resonance
are investigated. These dier depending on whether it is the phase or the frequency components of a
Fourier expansion of the motion that are taken to be continuous through the resonance. Asymptotic
and hyperasymptotic analysis are used to nd the late-time analytic behavior of the solution for a
system that has passed through a resonance. Linearly growing (weak resonances) or linearly decaying
(strong resonances) solutions are found depending on the strength of the resonance. In the weak-
resonance case, frequency resonances leave an imprint (a resonant memory) on the gravitational
frequency evolution. The transition between weak and strong resonances is characterized by a
square-root singularity, and as one approaches this transition from above, the solutions to the
frequency resonance equation bunch up into families exponentially fast.
PACS numbers: 04.30.-w,04.50.Kd,04.25.-g,04.25.Nx
I. INTRODUCTION
Divergent series are the invention of the devil, and
it is shameful to base on them any demonstration what-
soever. Niels Hendrik Abels 1828 statement [1] sug-
gests that asymptotic analysis, which commonly leads
to divergent series, should not be applied to problems
of physical interest. Asymptotics, however, has become
an invaluable tool for physicists seeking approximate an-
alytic solutions. Multiple-scale analysis, which includes
boundary-layer theory and WKB theory, allows us to un-
derstand diverse problems, such as semiclassical quantum
theory, airplane wing design and turbines [2]. In the con-
text of general relativity, asymptotic (post-Newtonian)
series [3, 4] constitute the basis of the lters used in cur-
rent gravitational wave detectors to extract signals from
the noise.
Resonances are a common occurrence in physical phe-
nomena. In a traditional oscillatory system, a resonance
is a point in frequency space where the system stores and
transfers energy between kinetic and potential modes, al-
lowing a small driving force to generate large amplitude
oscillations. Resonant phenomena can occur in many vi-
brational or wave-like systems and include electromag-
netic resonances, nuclear magnetic resonance, electron
spin resonance, and so on. In the realm of general rela-
tivity, black holes can sometimes be treated as resonators,
as they relax after being perturbed (see, for example, [5]).
In the context of general-relativistic orbital mechan-
ics, resonance has recently been adopted to represent a
slightly dierent phenomenon, namely, the enhancement
of gravitational-wave energy dissipation due to the lack
of cancellation of oscillatory modes that for generic in-
spirals average out [6, 7]. This is particularly relevant for
extreme mass-ratio inspirals (EMRIs), in which a small
compact object, such as a stellar-mass black hole or a
neutron star, orbits around a supermassive black hole [8].
In such a two-body system the smaller object slowly spi-
rals inwards due to gravitational-wave energy-momentum
losses, on a radiation-reaction timescale much longer
than the orbital one. This inspiral is usually modeled
by computing an orbit-average of the gravitational-wave
energy ux. This procedure discards terms proportional
to odd-powers of sines or cosines of the sum of the orbital
phases [912]. For orbital congurations or points in fre-
quency space at which the sum of the orbital frequencies
vanishes, the orbit-averaged energy ux is not equal to
the limit of the orbit-averaged uxes for nearby, nonres-
onant orbits. This is because harmonics of the frequency
that vanish on resonance contribute to the secular com-
ponent of the change in the orbital elements on resonance
[6, 7], but average to zero for o-resonance orbits. Un-
like traditional oscillators, however, there is no external
driving force in the EMRI case; the emitted gravitational
waves drive the inspiral themselves and the resonance is
2
caused by the orbital frequencies becoming commensu-
rate, that is, some linear combination of the three fre-
quencies with integer coecients vanishes at resonance.
A secularly growing radiation-reaction force can leave
strong imprints on the orbital motion, even if this secu-
lar growth is active for a very short time. These imprints
can then propagate into the gravitational waves emitted
and could have important consequences for gravitational-
wave detection. Unlike conventional telescopes operating
in the electromagnetic spectrum, current gravitational-
wave detectors will not observe signals above the average
noise. Instead, signals are expected to be buried deep
in the noise, and will be extracted using lters based on
the expected signals. Although EMRIs are not expected
to be detected with current ground-based gravitational
wave detectors, they are a key target for future space-
based detectors for which accurate EMRI lters will be
needed. It is therefore important to understand how
EMRI resonances can aect the emitted gravitational
waves.
Resonances in Extreme Mass-Ratio Inspirals
The extreme-mass-ratio (typically 10
6
10
5
) en-
sures that over short timescales the orbit of the smaller
object in an EMRI system is approximately geodesic. It
is therefore appropriate to use an osculating-element
formulation in which the EMRI is identied by a sequence
of geodesics [7, 12]. Geodesics in a Kerr background are
uniquely characterized by three constants of the motion,
energy E, z-component of angular momentum L
z
, and
Carter constant Q, and four initial phases that specify
the coordinates of the object at a particular time. The
position and velocity of an object uniquely identies a
geodesic. Since the evolution of the orbit is governed by
a second-order dierential equation, the values of these
seven geodesic parameters at each point on the inspiral
provide an alternative parametrization of the inspiral.
The time-evolution of the geodesic constants of the mo-
tion is
dJ

dt
= F
SF

(q, J) +O(
2
) , (1)
where J

= (E, L
z
, Q) is a vector of these constants,
while is the mass ratio, q is an angle phase variable,
and F
SF

is the self-force. The rate of change of J

can then be used to construct the rate of change of the

orbital frequencies in a similar form.
For any given geodesic, the self-force can be expanded
in a Fourier series in terms of the fundamental orbital
frequencies. These frequencies can be mapped to the
geodesic constants of the motion. A geodesic resonance
occurs when the ratio of the frequencies of the radial
and vertical motion is a rational number. The third fre-
quency, that of meridional motion, is not relevant for
resonances due to the axisymmetry for the background
Kerr spacetime. Henceforth, we only consider the depen-
dence of the self-force on the two frequencies that can
lead to a resonance.
Let us then expand the self-force in a two-frequency
Fourier series, where one of the frequencies () ap-
proaches zero while the other (f) remains nite:
d
dt
=

,n
G
n
cos [( +nf) t] +H
n
sin [( +nf) t] ,
(2)
where G
n
and H
n
are time-independent Fourier coef-
cients that depend on the orbital parameters. Clearly,
when +nf = 0, the cosine and sine terms average out
for suciently long integration times. However, at reso-
nance, where = nf (which in this case we take to
be n = = 0), the cosine function goes to unity, leaving
a sum of secular (zero-frequency) Fourier coecients.
Let us now further assume that the Fourier coecients
(G
n
, H
n
) vary smoothly as the resonance is approached,
such that they can be expanded as their on-resonance
values plus a correction of O(
1/2
). Such corrections can
always be made small by choosing a suciently small ,
independent of the magnitude of . Similarly, corrections
from other terms of O() on the right side of (1) can be
ignored.
We are now left with a number of rapidly oscillating
terms (those with + nf = 0) and also terms that
slowly oscillate away from resonance and then vanish at
resonance. The rapidly oscillating terms are less impor-
tant because they average to zero on a short timescale.
Changing variables to y /

and x

t, we nd
that Eq. (2) becomes
dy
dx
=

n
G
n
cos
_
xy +
nfx

_
+H
n
sin
_
xy +
nfx

_
.
(3)
Expanding the sum for the rst few (, n) modes and
rescaling by y = y/

G
00
, x =

G
00
x, we then have
y

= 1 +k cos (xy) , (4)

where we have dropped the tildes, prime denotes dier-
entiation with respect to x, and, in principle, the param-
eter k G
01
/G
00
is known. In most scenarios k 1,
but there could be orbits for which k = O(1).
In deriving (4), we have made several approximations:
(a) we have ignored the rapidly oscillating
1/2
terms,
and thus considered only the n = 0 modes; (b) we have
ignored a phase constant induced by the H
10
term; (c) we
have considered only the = 1 mode because these are
the dominant ones and are on resonance for the longest
time. Assumption (a) is justied, given that the rapidly
oscillating components tend to average out over a suf-
ciently long integration time. The relaxation of as-
sumptions (b) and (c) will be addressed more carefully
in Sec. V.
In addition to the approximations described above, (4)
also makes the critical assumption that the self-force can
be expanded as a Fourier series in the frequency with ar-
3
gument (+nf)t and that the coecients of this expan-
sion are continuous at resonance. An alternative way to
write the same equation o resonance would be as an ex-
pansion in the phase with argument +n, where (, )
are angle variables. For a geodesic, the time derivatives
of (, ) are the frequencies (

) (, f), but if one

regards these phase angles as fundamental and assumes
that the coecients of that expansion are continuous at
resonance, one ends up with a slightly dierent dieren-
tial equation:

= 1 +k cos , (5)
which admits the rst integral
1
2
(

)
2
= +k sin +

(0) , (6)
where

(0) is an integration constant.

Equations (4) and (5) give two alternative descriptions
of an EMRI resonance, but they are not equivalent. To
make this clear, we rewrite (5) in terms of y:
y

= 1 +k cos
__
ydx
_
. (7)
This equation is equivalent to (4) only in the limit
xy

1. In this paper we seek solutions to these two

dierential equations in the limit x +. The equa-
tions are deceptively simple (they are just ordinary dif-
ferential equations) but due to the nonlinearity, nding
exact solutions is impossible.
This paper describes the solution to both the
frequency-resonance dierential equation (4) and the
phase-resonance dierential equation (5) and is organ-
ised as follows: The leading-order behavior of the solu-
tion at late times (x ) for both frequency and phase
resonances are calculated in Sec. II. The higher-order
behavior in k of these solutions is given in Sec. III. Sec-
tion IV describes the qualitative change in behavior as
k transitions from k > 1 to k < 1. Section V discusses
generalisations of the resonance equations and explains
how the solutions are modied. Section VI gives some
conclusions and describes possible future work.
II. LEADING-ORDER ASYMPTOTIC
BEHAVIOR FOR LARGE x
Equations (4) and (7) describe simple models of nonlin-
ear resonant behavior and similar versions of these equa-
tions have been studied before. In fact, equations of the
form y

= f(x + y) or y

= f(y/x) have solutions in

quadrature. However, equations of the form y

= f(xy)
cannot be solved exactly. Instead, one relies on asymp-
totic techniques to understand their behavior.
The prototypical equation to study with these tools is
[2]
y

= cos xy , (8)
whose asymptotic expansion in the limit x + is
y(x)
a
x
(x +), (9)
where a = (n + 1/2) and n is an integer. For a slowly
varying solution, as x +, y

1, which implies that

cos xy 1 and thus xy = a (n+1/2). In fact, one can
show that corrections to this asymptotic solution scale
with powers of (1/x)
m
for m > 1 [2]. Similar techniques
can be used to show that the solution to equation y

=
tan 2xy also behaves as y (2n + 1)/x as x +
[13, 14].
A. Frequency Resonances
Let us rst consider the case k > 1 of (4). For slowly
varying solutions y

1, 1 +k cos xy 1 and then

y
arccos(1/k)
x
as x , k > 1. (10)
When k < 1, the above solution does not exist and more
subtle asymptotic techniques must be used.
Next, we consider the k < 1 case. The form of (4)
suggests that y y
c
ax as x +. Let us then try
this ansatz, which when inserted in (4) gives
a = 1 + cos(ax
2
) . (11)
Clearly, the y
c
ansatz is not a proper solution. We can
understand this by averaging the cosine term over all x
a = 1 +

cos(ax
2
)
_
, (12)
where the angle brackets stand for averaging. The Fresnel
cosine function is dened via the integral
C(x)
_
2a

_
1/2
_
x
0
cos(ax
2
)dx

. (13)
and as x +, C(x) 1/2+sin(x
2
/2)/(x). At very
large x then,
_
cos(ax
2
)dx
1
4
_
2
a
_
1/2
+
sin(ax
2
)
2ax
. (14)
Note, however, that the second term still depends on x,
so our ansatz y
c
is still not a valid solution.
These considerations motivate the improved ansatz
y y
1
as x +, where
y
1
(1 +c)x +b +
a
1
x
sin[(1 +c)x
2
+bx] (15)
with constants a
1
, b, c to be determined. Inserting this
ansatz into (4), we obtain
c + 2(1 +c)a
1
cos[(1 +c)x
2
+bx] +O
_
x
1
_
(16)
k cos[(1 +c)x
2
+bx] cos{a
1
sin[(1 +c)x
2
+bx)]}
k sin[(1 +c)x
2
+bx] sin{a
1
sin[(1 +c)x
2
+bx]} ,
4
where we have expanded the cosine function with stan-
dard trigonometric identities. Comparing terms of the
left and right sides of this equation, we see that a
1
=
O(k).
Let us now require that k 1. Since the sine and
cosine functions are bounded by unity and since a
1
=
O(k) 1, we know that a
1
sin[(1 + c)x
2
+ bx] 1 for
all x, and we can therefore expand the cosine and sine
functions on the righthand side of (16). Performing the
expansion, we nd that (16) becomes
c + 2(1 +c)a
1
cos[(1 +c)x
2
+bx] k cos[(1 +c)x
2
+bx]
a
1
k sin
2
[(1 +c)x
2
+bx] (17)
to O(x
1
, k
3
) from which we infer that a
1
k/[2(1 +c)].
We are then left with
c
a
1
k
2
_
1 cos[2(1 +c)x
2
+bx]
_
, (18)
which implies that c ka
1
/2, and thus, c k
2
/4.
Substituting this back into a
1
, we nd that a
1
k/2.
The second term of (18) is not included in c because it
must be canceled by terms of O(k
2
) in y, which we have
neglected here. Our solution to (4) then becomes
y
1
= b+
_
1
k
2
4
_
x+
k
2x
sin
__
1
k
2
4
_
x
2
+bx
_
, (19)
with remainders of O(x
1
, k
3
), and where b remains un-
determined and depends on the initial conditions. We
have solved (4) numerically in the range k (0, 0.5)
and x (0, 10
3
) and veried that indeed (19) is a good
approximation to the numerical solution, as we show in
Sec. III.
The frequency evolution described by (19) is particu-
larly interesting. At late times, the behavior of the fre-
quency is dominated by the term linearly proportional
to x, with all others becoming subdominant. The slope
of the frequency, however, is dependent on k. That is,
as the physical system goes through a resonance, it ac-
quires a slope correction that depends on the properties
of the resonance (that is, on k), a resonant memory of
sorts. If present in EMRIs, this resonant memory could
have a large impact on the gravitational wave phase as
the system traverses a resonance.
1. Matched Asymptotic Expansion
The constant b is xed by the initial condition imposed
at x = 0, which requires that a solution be valid in the
x 1 limit. Recall that the solution found in (19) is
valid in the x 1 limit and that it diverges as x 0.
Let us now look for a solution valid for kx 1, with
k 1, by using the ansatz
y(x) =

n=0
k
n
y
n
(x) , as kx 1, k 1 , (20)
where y
n
(x) are undetermined functions independent of
k. We could have expanded y in k(kx)
n
instead of in k
n
,
but this would lead to more complicated dierential equa-
tions, although the solutions would be the same. Here,
we choose initial conditions y
n
(0) = 0 for all n, but the
extension to more general initial conditions is trivial.
The zeroth-order solution (n = 0) satises
y

0
(x) = 1 y
0
(x) = x. (21)
The rst-order solution (n = 1) satises
y

1
(x) = cos(x
2
) y
1
(x) =
_

2
C
_
_
2

x
_
.
(22)
To next order in k, we substitute the solutions found thus
far into the dierential equation
cos x
2
+k y

2
= cos
_
x
2
+kx y
1
(x) +k
2
x y
2
(x)
_
. (23)
We can expand the cosine using the assumption kx 1
to nd that
y

2
(x) = xsin(x
2
) y
1
, (24)
which then leads to the solution
y
2
(x) =
1
2
_

2
cos(x
2
)C
_
x

8
C
_
2

x
_

x
4
.
(25)
Similarly, the equation satised by the third-order term
in the expansion is
y

3
(x) =
x
2
2
cos(x
2
) y
2
1
(x) xsin(x
2
) y
2
(x) , (26)
but this cannot be explicitly integrated. Putting together
all the pieces found so far, we get
y = x +k
_

2
C
_
_
2

x
_
+k
2
_
1
2
_

2
cos(x
2
)C(
x

8
C
_
2x

x
4
_
.
(27)
Let us now asymptotically match this solution to the
one found in (19). For such a procedure to be valid, a
buer zone must exist where both solutions are simulta-
neously valid. Since (19) was found by assuming x 1,
while (27) assumes that kx 1, this implies that a buer
zone does exist with extension 1 x k
1
. Asymp-
totic matching requires that we asymptotically expand
(19) in kx 1 and (27) in x 1, and then set these
two expansions equal order by order. To leading order,
we nd that
x +

2
4
k x +b , (28)
with remainders of relative O(1/x, kx). This immedi-
ately leads to
b

2
4
k. (29)
5
B. Phase Resonances
Let us rst consider the solution to the phase resonance
equation, (5), in the k > 1 case. A constant solution as
x 1 exists provided

= 0, which is satised if
arccos(1/k) ,

(0) arccos(1/k)
_
k
2
1 . (30)
The rst condition is the same one that xy had to satisfy
in the frequency resonance case, but the second condition
now imposes a constraint on the initial conditions. If
one chooses the initial conditions (0) =

(0) = 0, the
second constraint leads to
_
k
2
1 = 2n arccos(1/k) for n Z
+
(31)
where arccos(x) is the principal value of the inverse co-
sine, taking values in the range [0, ]. Expanding for
k 1 we nd the approximate solution
k
th
=

4
(1 n)+
1
4
_

2
8n
2
+ 16n
2

2
8

. (32)
Evaluating this expression for the rst few values of n =
1, 2, . . ., we nd that k
th
4.60378, 10.9499, . . ., which
are to be compared with the exact numerical solutions to
(31), which are k
th
= 4.60334, 10.9499, . . .. We see that
the error in the above asymptotic expansion goes roughly
as 1/k
5
.
Let us now consider the solution to (5) in the k < 1
case. As before, we concentrate on perturbative solutions
in k 1. The zeroth-order solution is found by setting
k = 0 in (5):
0
x
2
/2. The rst-order solution in
k can be found by postulating that

0
(x) +
1
(x) . (33)
Inserting this into (5) we have

1
= k cos
_
1
2
x
2
+
1
_
. (34)
We assume that
1
is subdominant relative to
0
, and so
we approximate the argument of the cosine as x
2
/2. We
can then solve exactly for
1
to nd that

1
= k
_
x

C
_
x

_
sin
_
x
2
2
__
, (35)
where C(x) is the Fresnel cosine function dened in (13).
We can now compare this solution to the one obtained
for frequency resonances. Dierentiating and taking the
x 1 limit, we nd that

x +

k
2
+O(x
1
) . (36)
Notice that as x this agrees with (19) in func-
tional form, but not in slope. That is, as the system goes
through a phase resonance, the slope is not corrected by
k. We will nd in Sec. III that this remains true as one
calculates the solution to higher order in k. Therefore,
although frequency resonances seem to induce a memory,
phase resonances do not.
III. HIGHER-ORDER ASYMPTOTIC
BEHAVIOR FOR x 1 AND k 1
In this section we consider higher-order solutions in the
x 1 limit for both phase and frequency resonances, and
then compare these to numerical solutions.
A. Frequency Resonances
In order to obtain a higher-order solution, we must
construct an ansatz that eliminates the x-dependent part
of the right side of (18). We thus pose the ansatz y
y
asy
= y
1
+y
2
, where y
1
is given in (19), while y
2
is
y
2

a
2
x
sin{2[(1 +c)x
2
+bx]} . (37)
Inserting this ansatz into (4), we nd that
c + 2(1 +c)a
1
C
1
+ 4(1 +c)a
2
C
2
+O(x
1
)
k cos[(1 +c)x
2
+bx] cos(a
1
S
1
+a
2
S
2
)
k sin[(1 +c)x
2
+bx] sin(a
1
S
1
+a
2
S
2
) , (38)
where we use the notation
S
n
sin{n[(1 +c)x
2
+bx]} ,
C
n
cos{n[(1 + c)x
2
+bx]} . (39)
As before, we note that S
n
and C
n
are bounded by unity,
and since a
1
= O(k) we expect that a
2
= O(k
2
) or
smaller. This suggests that we can expand the cosine
on the right side of the above equation as in Sec. II A to
obtain
c+(1+c) (2a
1
C
1
+ 4a
2
C
2
)
a
1
k
2
(1 +C
2
)+kC
1
, (40)
to O(x
1
, k
3
). We see then that our previous solution
still holds: c = ka
2
1
/2 and 2(1 + c)a
1
= k, implying
that c k
2
/4. We also see that 4(1 + c)a
2
= a
1
k/2,
which implies that a
2
= k
2
/[16(1 + c)
2
] or simply that
a
2
k
2
/16 when expanding in k 1. The second-order
solution therefore becomes
y
2

k
2
16x
sin
_
2
_
1
k
2
4
_
x
2
+ 2bx
_
. (41)
We can obtain the next-order solution by constructing
the ansatz y y
asy
= y
1
+y
2
+y
3
, where y
3

a3
x
S
3
and
where we assume that a
3
= O(k
3
). Inserting this into
(4), we get
c + 2(1 +c)a
1
C
1
+ 4a
2
(1 +c)C
2
+ 6a
3
(1 +c)C
3
kC
1
cos(a
1
S
1
+ a
2
S
2
+a
3
S
3
)
kS
1
sin(a
1
S
1
+a
2
S
2
+a
3
S
3
) . (42)
6
Expanding the above equations in a
1
1 and a
3
1,
we nd that
c + 2(1 +c)a
1
C
1
+ 4a
2
(1 +c)C
2
+ 6a
3
(1 +c)C
3

ka
1
2
+
_
k
ka
2
2

ka
2
1
8
_
C
1
+
ka
1
2
C
2
+
_
ka
2
1
8
+
ka
2
2
_
C
3
. (43)
Matching cosine coecients, this leads to the following
system of equations:
c
ka
1
2
,
2a
1
(1 +c) k
_
1
a
2
2

a
2
1
8
_
,
4a
2
(1 +c) k
a
1
2
,
6a
3
(1 +c) k
_
a
2
1
8
+
a
2
2
_
, (44)
which we can solve as an expansion in k to nd that
a
1

k
2
, a
2

k
2
16
,
a
3

k
3
96
, c
k
2
4

3
64
k
4
. (45)
Therefore, our solution to third order becomes
y
asy
= b +
_
1
k
2
4

3
64
k
4
_
x
+
k
2
1
x
sin
__
1
k
2
4

3
64
k
4
_
x
2
+bx
_
+
k
2
16
1
x
sin
_
2
_
1
k
2
4

3
64
k
4
_
x
2
+ 2bx
_
+
k
3
96
1
x
sin
_
3
_
1
k
2
4

3
64
k
4
_
x
2
+ 3bx
_
. (46)
Notice that this higher-order solution retains the reso-
nant memory computed in the previous section (that is,
the k
2
correction to the linear-in-x term, which is domi-
nant at late times).
This procedure can be generalized to arbitrary high
order in k 1. To leading order in 1/x, we make the
ansatz y y
(1)
, where
y
(1)
(1 +c)x +b +
k
x

n=1
a
(1)
n
S
n
+b
(1)
n
C
n
+O
_
1
x
2
_
,
(47)
where (a
(1)
n
, b
(1)
n
, b, c) are constant coecients that de-
pend on k. We can expand (a
(1)
n
, b
(1)
n
) as an expansion in
k, that is, a
(1)
n
=

m
a
(1)
nm
k
m
, and solve for these coe-
cients by equating the coecients of the S
n
s and C
n
s
at dierent orders in k. Because the derivative of (47) is
y

(1)
= 1+c+

n=1
2n(1 +c)k
x
_
a
(1)
n
C
n
b
(1)
n
S
n
_
+O
_
1
x
2
_
.
(48)
we can evaluate (4) to obtain
y

(1)
= 1 +k C
1
cos
_
k

n=1
a
(1)
n
S
n
+b
(1)
n
C
n
_
k S
1
sin
_
k

n=1
a
(1)
n
S
n
+b
(1)
n
C
n
_
. (49)
The generic sine and cosine Taylor expansion formula
allows us to rewrite the above equation as
y

(1)
= 1 +k C
1

=0
(1)

k
2
(2)!
_

n=1
a
(1)
n
S
n
+b
(1)
n
C
n
_
2
k S
1

=0
(1)
+1
k
2+1
(2 + 1)!
_

n=1
a
(1)
n
S
n
+b
(1)
n
C
n
_
2+1
. (50)
At this point, no further progress can be achieved because
one needs to evaluate the (2 +1)st and the (2)th power
of an innite series, which is not easy to do in closed form.
This is why it is more convenient to expand the rst few
terms in the series, as done earlier in this section.
We can generalize the previous procedure to higher
order in x. Postulate the ansatz y y
(1)
+y
(2)
, where
y
(2)

k
x
2

n=0
a
(2)
n
S
n
+b
(2)
n
C
n
. (51)
The derivative of y is then simply
y

= 1 +c +

n=1
2n(1 +c)k
x
_
a
(1)
n
C
n
b
(1)
n
S
n
_
(52)
k

n=1
_
2n(1 +c)b
(1)
n
+a
(1)
n
x
2
_
S
n
+ k

n=1
_
2n(1 +c)a
(1)
n
b
(1)
n
x
2
_
C
n
+O
_
1
x
3
_
,
while the right side of (4) implies that
y

= 1 +k C
1
cos
_
k

n=1
_
a
(1)
n
S
n
+b
(1)
n
C
n
_
_
k S
1
sin
_
k

n=1
_
a
(1)
n
S
n
+b
(1)
n
C
n
_
_
(53)

k
2
x
C
1
sin
_
k

n=1
_
a
(1)
n
S
n
+b
(1)
n
C
n
_
_

n=1
_
a
(2)
n
S
n
+b
(2)
n
C
n
_
+
k
2
x
S
1
cos
_
k

n=1
_
a
(1)
n
S
n
+b
(1)
n
C
n
_
_

n=1
_
a
(2)
n
S
n
+b
(2)
n
C
n
_
.
7
One could now expand these equations in k 1 and
equate coecients to get equations for the a
(2)
n
s and
b
(2)
n
s. Following this scheme, one can nd the subdom-
inant terms in the asymptotic expansion of the solution
as series in 1/x.
B. Phase Resonances
Let us now concentrate on higher-order solutions to
the phase-resonance equation in the k < 1 case. We thus
postulate the ansatz
asy
, where

asy
=
0
(x) +
1
(x) +
2
(x) , (54)
where we recall that
0
= x
2
/2, and we rewrite
1
as

1
= k

2
x +k
_

x
_
C
_
x

1
2
_
sin
_
x
2
2
__
.
(55)
We have factored out the unbounded-in-x part of
1
and
the second term is now bounded for all x and tends to 0
as x . This ensures that the term in curly brackets
is small for suciently small k as x . We can then
see that
2
must satisfy the dierential equation

2
k cos (
0
+
1
) . (56)
where we seek solutions accurate to O(k
2
) and, as in
Sec. II B, we neglect the
2
term in the source. Inserting

0
and
1
and using the fact that the bracketed term in
(55) is everywhere small to expand the cosine, we nd
that

2
= k cos
_
x
2
2
+k

2
x
_
k
2
sin
_
x
2
2
+k

2
x
_

x
_
C
_
x

1
2
_
sin
_
x
2
2
__
. (57)
Integrating this equation twice, imposing the condition
that
2
(0) =

2
(0) = 0, and ignoring terms explicitly
proportional to k
3
or higher, we obtain

asy
=
x
2
2
+k

x
_
1
2
C
_
k
2
_
+k
(

2 1)
2

2
_
+

8
k
2
+k
_

_
x +

k
2
__
C
_
x

+
k
2
_

1
2
_
sin
_
1
2
_
x +k

2
_
2
__
+k
2
_
x
_

2
_
C
_

2x

1
2
_
+
1
2
sin(x
2
)
+

2
_
C
_
x

__
2

2
C
_
x

_
+

8
_
. (58)
The slope of this solution for large x gives the k-
correction to the gradient

asy
x +

2
k
_
1
k

2
+O(k
2
)
_
. (59)
Notice that although the constant is k-corrected, the
linear-in-x term is not, showing again that phase reso-
nances do not acquire a resonant memory imprint.
We proceed to higher order, and by analogy we write
down the solution for
2
(x) as the sum of the part on the
rst line of (58) that grows linearly with x and a part that
is bounded for all x and has a convergent integral on the
range [0, ]. Schematically, this takes the form

1
+
2
=

8
k
2
+
_

2
k

2
k
2
_
x+k
1
(x)+k
2

2
(x),
(60)
where k
1
(x) is the term on the second and third lines
of (58), while k
2

2
(x) is the term on the fourth and fth
lines. Notice that k
1
(x) does contain terms of O(k
2
).
We now seek the next order solution,
3
, that solves the
equation

3
= k cos
_
1
2
x
2
+
1
+
2
_
, (61)
where again we have neglected
3
in the source term.
Inserting the solution known so far and expanding the
cosine, keeping terms up to O(k
3
), we nd that

3
= k cos
_
1
2
x
2
+
_

2
k

2
k
2
_
x +

8
k
2
_

k
3
2
cos
_
1
2
x
2
_

2
1
k
2
sin
_
1
2
x
2
+

2
kx
_

1
k
3
sin
_
1
2
x
2
_

2
. (62)
Note that in each term we have eliminated terms in the
arguments of the cosine and sine that are lower order
than k
3n
, where n is the order of the k prefactor to the
term. By integrating these terms over the range [0, ],
we can derive the O(k
3
) correction to the asymptotic gra-
dient. We nd the contributions to the O(k
3
) term in the
gradient from each line of (62) are

/2

2, 0.06202,
0, and

(1 2

2 +

3)/8 0.07844 respectively. The

nal form for the asymptotic gradient is then

asy

2
k

2
k
2
+ 0.46484k
3
. (63)
C. Comparison to Numerical Results
The left panel of Fig. 1 shows the dierence between
the numerical solution for y and the asymptotic expan-
sion in (46) to all computed orders in k (solid line) and
8
FIG. 1. Left: Dierence between the numerical solution, y, to the frequency resonance equation and the asymptotic expansion
in (46) to all computed orders in k (solid line) and to linear order in k (dotted lines) as a function of x. Right: Same dierence
as left-panel but for the solution, , to the phase resonance equation. Observe that in both cases the full solution does much
better than the O(k) truncation.
to linear order in k (dotted lines). In all plots we have
chosen b to be that found in Sec. II A1. As expected, the
solid lines are much closer to zero than the dotted ones.
Moreover, notice that the asymptotic solutions found in
the limit x 1 are already quite good (better than 1%
relative to the numerical solution) below x < 10.
The right-panel of Fig. 1 shows the same type of dif-
ference as the left-panel, but for the phase solution. As
before, observe that the full solution agrees with the nu-
merical solution much better than its truncated version.
Unlike the y
asy
solution, the
asy
is globally valid, as we
did not restrict attention to the x 1 limit, and this
can be clearly seen in Fig. 1. It is also apparent that, as
expected, in both cases the smaller k is, the better the
asymptotic solution.
Although the previous gures establish that the
asymptotic solutions are indeed accurate representations
of the numerical ones, they do not compare the y and
solutions to each other. Figure 2 shows the dierence be-
tween y and d/dx computed numerically in both cases.
As predicted by the asymptotic solutions, the dierence
is approximately k
2
x/2 (the dotted lines) in the x 1
limit.
IV. ASYMPTOTIC TRANSITION AT k = 1
As described in the previous sections, both the y and
solutions experience a transition as k 1. In this
section we discuss this transition in more detail using
hyper-asymptotic tools.
A. Frequency Resonances
Let us try to understand the fundamental change in
the asymptotic behavior of the y solution as k transitions
FIG. 2. Dierence between the numerical solution for y and
the x-derivative of the numerical solution for . For compar-
ison, we also plot the asymptotic slopes k
2
x/2.
from k > 1 to k < 1 from above. To do so, we must nd
the most important term as k approaches 1 to all orders
in powers of 1/x. Then, we must sum the series and
identify the singularity at k = 1.
Assuming that k > 1, we know that the behavior of
y(x) in (4) as x is described by a series in inverse
odd powers of x:
y
a
x
+
b
x
3
+
c
x
5
+
d
x
7
+ x , (64)
where
cos(a) =
1
k
. (65)
As discussed earlier, as k 1
+
, a (2n+1), and thus
sin a approaches 0. To higher order in 1/x, one easily
9
nds that
b =
a
k sin(a)
, (66)
c =
6a sin(a)k +a
2
2k
3
[sin(a)]
3
. (67)
Thus, when sin(a) is small, the most singular part of c is
c
a
2
2k
3
[sin(a)]
3
. (68)
Similarly, the most singular part of d, e, f, . . . is
d
a
3
2k
5
[sin(a)]
5
, e
5a
4
8k
7
[sin(a)]
7
,
f
7a
5
8k
9
[sin(a)]
9
, g
21a
6
16k
11
[sin(a)]
11
,
h
33a
7
16k
13
[sin(a)]
13
. (69)
The numerical coecients, 1,
1
2
,
1
2
,
5
8
, . . ., are given by
a very simple formula:
F(n) =
(2n)!
n!(n + 1)!2
n
(n = 0, 1, 2, 3, ...). (70)
Therefore, if we sum the most singular terms as k 1
+
to all orders in powers of 1/x, we nd that
y(x)
a
x
+
a
kx
3
sin(a)

n=0
F(n)
_
a
x
2
k
2
[sin(a)]
2
_
n
. (71)
From the Taylor expansion
1
z

1 2z
z
= 1 +
1
2
z +
1
2
z
2
+
5
8
z
3
+. . . =

n=0
F(n)z
n
(72)
one nds that
y(x)
1
x
_
a +k sin(a) k sin(a)

1
2a
x
2
(k
2
1)
_
,
(73)
where we have used the identity k
2
[sin(a)]
2
= k
2
1. This
shows that there is a square-root-branch-cut singularity
where the asymptotic behavior goes complex as k 1
+
.
In addition to this branch-cut singularity, one can also
show that the higher-order terms in 1/x bunch up into
families as x , with each pair of families separated
by an unstable separatrix curve. As shown earlier, as
k 1
+
, the leading-order slope of the solution a
(2n + 1) and sin a approaches 0. There are, however,
many solutions for a in (65) as k 1
+
: the rst lies
just below (but above /2); the second lies just above
; the third and fourth lie just below and just above 3,
and so on.
Consider two dierent solutions, y
1
(x) and y
2
(x) corre-
sponding to one of the innite number of possible values
of a and dene
Y (x) y
1
(x) y
2
(x). (74)
Observe that Y (x) satises the dierential equation
Y

(x) = k cos [xy

1
(x)] k cos [xy
2
(x)] . (75)
Using the identity
cos cos = 2 sin
_
+
2
_
sin
_

2
_
, (76)
we can rewrite (75) as
Y

(x) 2k(sina) sin

_
xY (x)
2
_
, (77)
for large x.
Let us now make the assumption that y
1
(x) and y
2
(x)
approach one another as x gets large, so that Y (x) is
small when x >> 1. Then, this dierential equation
becomes
Y

(x) x(tan a) Y (x), (78)

whose solution is
Y (x) C e
1
2
x
2
tan a
. (79)
Note that this solution is growing exponentially if tan a
is positive, and thus the assumption that Y (x) is small
as x is not valid. This is the unstable (separa-
trix) case. However, if tana is negative, then we have
the stable case, and we have shown that the family of
solutions corresponding to this case all bunch together
exponentially fast . Note that there is an alternation be-
tween stable and unstable behavior: Stable behavior oc-
curs only for the values of a that are just below (2n+1)
for n Z, that is, , 3, 5, and so on, while unstable
behavior occurs for the values of a just above (2n +1).
B. Phase Resonances
The phase solutions also show a fundamental quali-
tative change of behavior as k transitions from k > 1
to k < 1. Because the phase dierential equation, (5),
admits a rst integral of the motion, the analysis is sim-
pler than in the frequency case and it does not require
a hyper-asymptotic analysis. This is best understood if
we consider the rst-order form of (5), which we rewrite
here as
1
2
(

)
2
= E V (; k) , (80)
where V (; k) = k sin and E =

(0). The po-

tential is shown for representative values of k in the left
panel of Figure 3. Motion can only exist in regions where
E > V (). The initial conditions (0) =

(0) = 0 cor-
respond to E = 0. For k < 1, the eective potential
has no turning points and for any initial conditions the
motion will be unbounded with +. When k > 1,
the potential does have turning points. The right panel
10
-20
-15
-10
-5
0
5
10
0 2 4 6 8 10 12
V
(

;
k
)

k = 0.1
k = 1.1
k = 4.60334
k = 10.9499
-3.18
-3.16
-3.14
-3.12
-3.1
-3.08
-3.06
-3.04
-3.02
-3
2 2.5 3 3.5 4
V
(

;
k
)

E = -3.1
E = -3.11303
E=-3.15
E = -3.17015
FIG. 3. Eective potential for the phase resonance equation, as dened by Eq. (80). The left panel shows the potential for
four dierent values of k, while the right panel shows a close up of the region 2 4 when k = 1.1. Each horizontal line
corresponds to a particular choice of the constant E. Motion can only exist where E > V (; k), as indicated by the solid parts
of the lines shown.
of Figure 3 shows a close-up of the potential for k = 1.1
in the vicinity of the turning points. For E > 3.11303,
the motion is unbounded as in the k < 1 case. However,
for 3.11303 > E > 3.17015, the motion intersects the
potential twice, and it can be oscillatory for suitable ini-
tial values of , or unbounded if (0) is suciently large.
For E < 3.17015, the motion is again only unbounded
for suciently large (0).
These specications for E and (0) place restrictions
on the initial conditions, which are inconsistent with the
conditions we want to impose. For the conditions (0) =

(0) = 0, the motion is unbounded for k 4.60334, the

critical value computed in (32). When k 4.60334, the
potential intersects the E = 0 axis a second time. The
motion will asymptotically approach the limiting value

4.49341. For k 4.60334, the motion is oscillatory.

At the next limiting solution for k, k 10.9499, the
eective potential has another intersection with the E =
0 axis. However, this region is inaccessible to motion with
these initial conditions and the motion is still oscillatory.
Starting the motion with

(0) = 0 and 5.73224 (0)

10.9041 would generate a solution that asymptotically
approaches the next limiting value

10.9041.
At the threshold value k 4.60334, two solutions

1
(x),
2
(x) with 0 > E
1
= E
2
< 0 will oscillate with
dierent frequencies, and we therefore expect the dier-
ence
1
(x)
2
(x) to be oscillatory. If E
1
< 0 < E
2
,
one solution will be unbounded and therefore the dier-
ence will grow like x
2
/2. For 0 < E
1
= E
2
> 0, both
solutions are unbounded, and we expect the dierence to
grow linearly with x.
V. GENERALIZED RESONANCES
Section I derived certain equations [Eqs. (4) and (5)]
that are representative of phase and frequency resonances
in EMRIs, but in doing so we made two important sim-
plifying assumptions: (b) we ignored the sine term in the
sum given in Eq. (3) and (c) we ignored higher- terms
in this same sum. In this section we relax these two as-
sumptions and discuss how the solutions are modied.
Let us rst relax assumption (b). If the sine term is
included, (4) and (5) can be written as
dy
dx
= 1 +k cos(xy +),
d
2

dx
2
= 1 +k cos( +) ,
(81)
where is a constant. We can repeat the analysis of
the frequency evolution equation with the modication
introduced above by making the ansatz
y
1
(1 +c)x +b +
a
1
x
sin[(1 +c)x
2
+bx +] (82)
in place of (15). The analysis proceeds exactly as before,
but with the arguments of the various cosine and sine
terms modied via (1+c)x
2
+bx (1+c)x
2
+bx+. The
asymptotic slope, (1+c), is unchanged as a function of k.
The solution for b will be modied, however, because the
expansion for kx 1 described in Sec. II A1 is modied.
In particular
y
0
(x) = x,
y
1
(x) = cos()
_

2
C
_
_
2

x
_
sin()
_

2
S
_
_
2

x
_
y
2
(x) =
1
2
y
1

x
4

8
_
cos()C
_
2

x
_
sin()S
_
2

x
__
, (83)
in which S() denotes the Fresnel sine function, dened
by a similar equation to Eq. (13), but with the cosine
replace by a sine. After asymptotic matching, we then
nd that b

(k/2) cos(/4 +).

11
In the phase-resonance case, the addition of the
to the equation of motion is equivalent to solving the
original problem with a modied initial condition: (0) =
and

(0) = 0. However, the solution to the modied

equation with the standard initial condition (0) = 0
can also be found straightforwardly using the method
described in this paper. In that case, the zeroth-order-in-
k solution is unchanged,
0
= x
2
/2, but the rst-order-
in-k correction,
1
, is modied to

1
= k

x
_
cos()C
_
x

_
sin()S
_
x

__
k sin
_
x
2
2
+
_
+k sin() (84)
from which we see that the asymptotic correction to the
gradient is modied to k
_
/2 cos( + /4). Continu-
ing to the next order in k, we obtain for the asymptotic
gradient
d
dx
1 +k
_

2
cos
_
+

4
_

4
k
2
_
1 + (2

2) cos
_
2 +

4
__
(85)
for x 1. We note that the leading-order-in-k correction
to the asymptotic gradient in the phase equation and the
leading-order correction to b in the frequency equation
can be made to vanish for = /4. However, even with
this choice, the next-to-leading-order correction does not
vanish, and so although the eect of the resonance can
be suppressed for certain values of , it cannot be elim-
inated.
Let us now relax assumption (c); that is, let us in-
clude higher modes in the resonant dierential equa-
tions. Consider rst an evolution equation of the form
dy
dx
= 1 +k cos( xy) , (86)
where n is some integer. By writing Y =

y, X =

x,
we can rewrite this equation as
dY
dX
= 1 +k cos(XY ), (87)
reducing it to the same form that we considered before.
We note in particular that since the solution to this equa-
tion behaves as Y (1+c)X for x 1, this implies that
y (1 + c)x and hence the value of c is unchanged as a
function of k. The value of b would be modied, however.
In the phase case, if we modify the equation to
d
2

dx
2
= 1 +k cos( ) (88)
with n again an integer, the substitution = , X =

x gives
d
2

dX
2
= 1 +k cos(). (89)
Using the solution to this equation that we found earlier,
we nd that X
2
/2 + (k)X + and we deduce
that

1
2
x
2
+
(k)

X + . (90)
Clearly then, higher- modes in the sum of (3) are sup-
pressed by a factor of 1/

.
Let us now consider the case where we have more than
one term on the righthand side of the evolution equa-
tion. In the frequency resonance case, we would have an
equation like
dy
dx
= 1 +k
1
cos(xy) +k

cos( xy). (91)

To solve this equation, we can proceed as before, making
an ansatz of the form
y
1
(1 +c)x +b +
a
(1)
1
x
sin[(1 +c)x
2
+bx]
+
a
()
1
x
sin[((1 +c)x
2
+ bx)]. (92)
Inserting this into the dierential equation, we nd the
same source terms involving sin
2
() that we found earlier,
one from each mode. We also nd various cross source
terms, but these do not contribute to the zero-frequency
part of the solution because they are products of oscil-
latory functions with unequal frequencies. We deduce
then that the solution for c at leading order is the sum
of the solutions treating each of the modes individually.
At the next order, cross terms will come in that may be
important, but these will be subdominant and further
exploration of these is beyond the scope of this paper.
In the case of the phase equation we would have
d
2

dx
2
= 1 +k
1
cos() +k

cos(). (93)
The linear-in-k term in the solution of this equation can
readily be seen to be the sum of the linear-in-k solutions
to the equation with only one of the cosine terms on
the righthand side. Cross terms again come in at higher
order, but by the preceding argument the size of the cor-
rections from the k

term are suppressed relative to the

dominant mode by factors of 1/

.
VI. CONCLUSIONS
We have studied the behavior of the solution to two dif-
ferential equations that describe the gravitational-wave
phase and frequency evolution during an EMRI that ex-
periences a resonant transition. We have found two gen-
eral dierential equations that might describe this behav-
ior: one based on the assumption that the coecients
of a Fourier expansion of the self-force in frequency are
continuous at a resonance, and the other based on the
12
assumption that it is the Fourier phase coecients that
are continuous at resonance. We have solved both dier-
ential equations at late times using asymptotic methods.
Depending on the strength of the resonance (controlled
by the parameter k), both resonant equations lead to so-
lutions with qualitative dierent behavior: Weak reso-
nances (k < 1) lead to linear temporal growth of the fre-
quency and quadratic growth of the orbital phase; strong
resonances (k > 1) lead to linear temporal decay of the
frequency function, leading to a constant-phase oset at
late times.
Even though the dierential equations for phase and
frequency resonances might both describe the evolution
of an EMRI through a resonance, we have found that
the evolution depends on which dierential equation one
assumes. That is, at late times and in the weak reso-
nance case, the evolution of the frequency in the phase
and frequency resonant cases possess dierent asymptotic
slopes. The dierence in slope depends on the value of
k, with the frequency evolution acquiring a k-dependent
memory in the frequency resonant case that is absent for
phase resonances. Further work is required to explore
which of the two equations is in fact most applicable to
the EMRI resonance problem.
We also studied the transition between weak and
strong resonances. We found that at the transition point
k = 1, there is a square-root branch cut in the solution
to the frequency resonance equation. Close to this point,
we proved that frequency solutions bunch up into fami-
lies that decay as 1/x exponentially fast. In fact, there
is an alternation between stable (bunching up of solu-
tions) and unstable behavior, depending on the branch
of solutions considered.
Future work should concentrate on exploring which of
these equations is applicable to EMRI evolutions in prac-
tice and what the implications are for the construction
of waveform template models of EMRI signals. The ex-
istence of a memory eect in the frequency resonances
is particularly interesting and would have a profound
impact on our ability to detect EMRI signals. An ap-
proximate post-Newtonian prescription for the self-force
on resonance has been suggested [6] and would provide
a suitable framework in which to explore these ques-
tions further. Whichever of the two equations applies
to the problem that motivated this work, the results de-
scribed in this paper provide important insights into the
behaviour of the solutions to these dierential equations
and predictions for the change in the frequency and phase
of the evolution as the orbit passes through a resonance.
These results will be invaluable for constructing approx-
imate models to describe the evolution of EMRI orbits.
ACKNOWLEDGMENTS
CMB is supported by grants from the Leverhulme
Foundation and the U.S. Department of Energy. JGs
work is supported by the Royal Society. NY acknowl-
edges support from NASA grant NNX11AI49G, under
sub-award 00001944 and NASA through the Einstein
Postdoctoral Fellowship Award Number PF0-110080 is-
sued by the Chandra X-ray Observatory Center, which is
operated by the Smithsonian Astrophysical Observatory
for and on behalf of NASA under contract NAS8-03060.
JG thanks the MIT Kavli Institute for Astrophysics and
NY thanks the Yukawa Institute for Theoretical Physics
for their hospitality while this paper was being nished.
We also thank Scott Hughes for useful discussions.
[1] J. P. Boyd, Acta Applicandae 56, 1 (1999).
[2] C. M. Bender and S. A. Orszag, Advanced mathematical
methods for scientists and engineers 1, Asymptotic meth-
ods and perturbation theory (Springer, New York, 1999).
[3] N. Yunes and E. Berti, Phys. Rev. D77, 124006 (2008),
arXiv:0803.1853 [gr-qc].
[4] Z. Zhang, N. Yunes, and E. Berti,
Phys.Rev. D84, 024029 (2011).
[5] E. Berti, V. Cardoso, and A. O.
Starinets, Class.Quant.Grav. 26, 163001 (2009),
arXiv:0905.2975 [gr-qc].
[6] E. E. Flanagan and T. Hinderer, (2010),
arXiv:1009.4923 [gr-qc].
[7] J. R. Gair, E. E. Flanagan, S. Drasco, T. Hinderer, and
S. Babak, Phys.Rev. D83, 044037 (2011).
[8] P. Amaro-Seoane et al., Class. Quantum Grav. 24, 113
(2007).
[9] E. Poisson, Phys. Rev. D47, 1497 (1993).
[10] S. A. Hughes, Phys. Rev. D61, 084004 (2000).
[11] S. A. Hughes, Phys. Rev. D64, 064004 (2001).
[12] A. Pound and E. Poisson, Phys. Rev. D77, 044013
(2008).
[13] C. M. Bender, D. W. Hook, P. N. Meisinger, and Q.-h.
Wang, Phys.Rev.Lett. 104, 061601 (2010).
[14] C. M. Bender, D. W. Hook, P. N. Meisinger, and Q.-h.
Wang, Annals Phys. 325, 2332 (2010).