Scalar self-force on eccentric geodesics in Schwarzschild spacetime: A time-domain
computation
Roland Haas
Department of Physics, University of Guelph, Guelph, Ontario, Canada N1G 2W1
(Dated: April 3, 2007)
We calculate the self-force acting on a particle with scalar charge moving on a generic geodesic
around a Schwarzschild black hole. This calculation requires an accurate computation of the retarded
scalar field produced by the moving charge; this is done numerically with the help of a fourth-order
convergent finite-difference scheme formulated in the time domain. The calculation also requires
a regularization procedure, because the retarded field is singular on the particle’s world line; this
is handled mode-by-mode via the mode-sum regularization scheme first introduced by Barack and
Ori. This paper presents the numerical method, various numerical tests, and a sample of results for
mildly eccentric orbits as well as “zoom-whirl” orbits.
PACS numbers: 04.25.-g, 04.40.-b, 41.60.-m, 45.50.-j, 02.60.Cb, 02.70.Bf
I. INTRODUCTION
The inspiral and capture of solar-mass compact objects
by supermassive black holes is one of the most promis-
ing and interesting sources of gravitational radiation to
be detected by the future space-based gravitational-wave
antenna LISA [1]. For these extreme mass-ratio inspirals,
one can treat the compact object as a point mass and de-
scribe its influence on the spacetime perturbatively. Go-
ing beyond the test mass limit, its motion is no longer
along a geodesic of the unperturbed spacetime of the cen-
tral black hole; it is a geodesic of the perturbed space-
time created by the presence of the moving body. When
viewed from the unperturbed spacetime, the small body
is said to move under the influence of its gravitational
self-force. The self-force induces radiative losses of energy
and angular momentum, which will eventually drive the
object into the black hole. To describe the motion of the
body, including its inspiral toward the black hole, we seek
to evaluate the self-force and calculate its effect on the
motion. One way of doing this uses the mode-sum reg-
ularization procedure introduced by Barack and Ori [2].
(For a comprehensive introduction of the problem, see
the special issue of Classical and Quantum Gravity [3].)
In this paper, in an effort to build expertise to calculate
the gravitational self-force, we retreat to the technically
simpler problem of a point particle of mass m endowed
with a scalar charge q orbiting a Schwarzschild black hole
of mass M . Following up on a previous paper [4], we
implement the numerical part of the regularization pro-
cedure for generic orbits with a time-domain integration
of the scalar-wave equation.
A. The problem
Our goal is to calculate the regularized self-force acting
on a scalar point charge in orbit around a Schwarzschild
black hole. In analogy with the gravitational case, where
in a first-order (in m/M) perturbative calculation the
particle moves on a geodesic of the background space-
time, we take the orbit of the particle to be a geodesic and
calculate the self-force as a vector field on this geodesic.
We start by writing the Schwarzschild metric using the
tortoise coordinate r∗ = r + 2M ln
ds2 = f
−dt2 + dr∗2
+ r2dΩ2, (1.1)
where f =
1− 2M
, dΩ2 =
dθ2 + sin2 θdφ2
is the
metric on a two-sphere, and t, r, θ and φ are the usual
Schwarzschild coordinates. Our task is to solve the scalar
wave equation
gαβ∇α∇βΦ(x) = −4πµ(x), (1.2)
µ(x) = q
δ4(x, z(τ))dτ , (1.3)
where ∇α is the covariant derivative compatible with the
metric gαβ , Φ(x) is the scalar field created by a scalar
charge q which moves along a world line γ : τ 7→ z(τ)
parametrized by proper time τ . The source term µ(x)
appearing on the right-hand side is written in terms of a
scalarized four-dimensional Dirac δ-function δ4(x, x
′) :=
δ(x0 − x′0)δ(x1 − x′1)δ(x2 − x′2)δ(x3 − x′3)/
− det(gαβ).
Because of the singularity in the source term, the re-
tarded solution to Eq. (1.2) is singular on the world line,
and the näıve expression for the self-force,
Fα(τ) = q∇αΦ(z(τ)), (1.4)
must be regularized. Following DeWitt and Brehme [5],
Mino, Sasaki, Tanaka [6], Quinn and Wald [7], Quinn [8]
carried out this regularization for the electromagnetic,
scalar and gravitational radiation reaction. In later work,
Detweiler and Whiting [9] introduced a very useful de-
composition of the retarded solution of Eq. (1.2) in terms
of a singular part ΦS and a regular remainder ΦR:
Φ = ΦS +ΦR. (1.5)
ΦR is regular and differentiable at the position of the par-
ticle, satisfies the homogeneous wave equation associated
http://arxiv.org/abs/0704.0797v2
with Eq. (1.2), and is solely responsible for the self-force
acting on the particle. ΦS , on the other hand, satisfies
Eq. (1.2), is just as singular at the particle’s position as
the retarded solution, and produces no force on the par-
ticle. Rearranging Eq. (1.5) and differentiating once, we
can write the regularized self-force as
Fα := q∇αΦR = q
∇αΦ−∇αΦS
. (1.6)
In a previous paper [4], we described our implemen-
tation of the regularization procedure to find a mode-
sum representation of ∇αΦS along a generic geodesic of
the Schwarzschild spacetime. Schematically, we intro-
duce a tetrad eα
and decompose the tetrad components
Φ(µ) := e
∇αΦ of the field gradient in terms of ordinary
scalar spherical harmonics Yℓm:
Φ(µ)(t, r, θ, φ) =
Φℓm(µ)(t, r)Yℓm(θ, φ). (1.7)
Each mode Φℓm
(t, r) is finite at the position of the par-
ticle, but their sum diverges on the world line. In [4], we
derive analytic expressions for the mode-sum decompo-
sition of ΦS(µ),
ΦS(µ) =q
ΦS(µ),ℓ (1.8)
ΦS(µ),ℓ = A(µ)
+B(µ) +
(ℓ− 1
)(ℓ + 3
+ · · · , (1.9)
where the coefficients A(µ), B(µ), C(µ), and D(µ) are in-
dependent of ℓ; they are listed in Appendix B for conve-
nience.
As each mode of Φ is finite, it is straightforward to
compute the modes of the retarded solution using nu-
merical methods, and we will describe how this was done
in Sec. IV. We use the numerical solutions in Eq. (1.6)
to calculate the regularized self-force, regularizing mode-
by-mode:
ΦR(µ) =
Φ(µ),ℓ − ΦS(µ),ℓ
, (1.10)
where Φ(µ),ℓ :=
Yℓm (no summation over ℓ im-
plied).
For numerical purposes it is convenient to define ψℓm
Φ(x) =
ℓm, (1.11)
where Yℓm are the usual scalar spherical harmonics. Af-
ter substituting in Eq. (1.2), this yields a reduced wave
equation for the multipole moments ψℓm:
−∂2tψℓm + ∂2r∗ψℓm − Vℓψℓm =
− 4πq
Ȳℓm(π/2, φ0)δ(r
∗ − r∗0), (1.12)
where
Vℓ = f
ℓ (ℓ+ 1)
. (1.13)
An overbar denotes complex conjugation, E = −ut is the
particle’s conserved energy per unit mass, and uα = dz
is its four velocity. Quantities bearing a subscript “0” are
evaluated at the particle’s position; they are functions of
τ that are obtained by solving the geodesic equation
uβ∇βuα = 0 (1.14)
in the background spacetime. Without loss of general-
ity, we have confined the motion of the particle to the
equatorial plane θ = π
Once we have numerically solved Eq. (1.12), we ex-
tract numerical estimates for ψℓm, ∂tψℓm and ∂r∗ψℓm,
which can then be used to find Φℓm, ∂tΦℓm and ∂rΦℓm.
These—together with the translation table displayed in
Eqs. (1.23)–(1.26) of [4], reproduced in Appendix A—
allow us to find the tetrad components Φ(µ)ℓm with re-
spect to the tetrad defined by Eqs. (1.18)–(1.21) of [4].
Eventually we regularize the multipole coefficients
Φ(µ)ℓ =
Φ(µ)ℓm(t0, r0)Y
ℓm(π/2, φ0) (1.15)
using Eq. (1.10); this involves the regularization param-
eters listed in Eqs. (1.30)–(1.45) of [4], which are repro-
duced in Appendix B.
B. Organization of this paper
In Sec. II we introduce the main ideas behind the
discretization scheme used in the numerical simulation.
Sec. III describes the choices we make in order to handle
the problems of specifying initial data and proper bound-
ary conditions. The next section—Sec. IV—provides de-
tails on the concrete implementation of the ideas put
forth in Secs. II and III. In Sec. V we describe the tests
we performed in order to validate our implementation of
the numerical method. Sec. VI finally presents sample
results for a small number of representative simulations.
C. Future work
This work, which deals with a scalar charge moving in
the Schwarzschild spacetime, is not intended to produce
physically or astrophysically interesting results. Instead,
its goal is to help us evaluate the merits of several strate-
gies that could be used to tackle the more interesting (and
difficult) problems of electromagnetism and gravity.
One future project we are currently exploring is to ap-
ply the formalism developed so far to the electromagnetic
self-force acting on an electric charge. Beyond the tech-
nical complication of having to deal with a vector field
instead of a single scalar quantity, we are also faced with
the reality of having to impose a gauge (in our case: the
Lorenz gauge) and to eliminate (or at least control) gauge
violations in the numerical simulation. The first step,
namely, the calculation of the regularization parameters
A(µ), B(µ), C(µ), and D(µ) for the self-force, is currently
underway. Also underway is the calculation of the regu-
larization parameters for he gravitational self-force.
Another project is the implementation of a scheme to
use the calculated self-force to update the orbital pa-
rameters of a particle on its inspiral toward the black
hole. The standard proposed approach to this problem
in the past has been to calculate the self-force on a set
of geodesics which are momentarily tangent to the par-
ticle’s trajectory. The self-force calculated in this way
is then used to update the orbital elements. This “after
the fact” calculation of the motion requires one to build
(in advance) a large database of self-force values for the
anticipated set of orbital parameters that the particle’s
trajectory will assume during its inspiral. Alternatively,
and conceptually more simply, the self-force could be cal-
culated self-consistently along the real, accelerated tra-
jectory. Such an approach requires changes in the expres-
sions of the regularization parameters, which so far have
been derived only for geodesic orbits. We are currently
investigating the merits of such an approach.
II. NUMERICAL METHOD
In this section we describe the algorithm used to inte-
grate the reduced wave equation [Eq. (1.12)] numerically.
For the most part we use the fourth-order algorithm in-
troduced by Lousto [10], with some modifications to suit
our needs. We choose to implement a fourth-order con-
vergent code because second-order convergence for the
potential Φ, while much easier to achieve, would guaran-
tee only first-order convergence for ∇αΦ, the quantity in
which we are ultimately interested. With a fourth-order
convergent code we can expect to achieve third-order con-
vergence for ∇αΦ, which is required for an accurate es-
timation of the self-force. Numerical experiments, how-
ever, show that in practice we do achieve fourth-order
convergence for the derivatives of Φ, a fortunate outcome
that we exploit but cannot explain.
From now on, we will suppress the subscripts ℓ and m
on Vℓ and ψℓm for convenience of notation. The wave
equation consists of three parts: the wave-operator term
(∂2r∗ − ∂2t )ψ and the potential term V ψ on the left-hand
side, and the source term on the right-hand side of the
equation. Of these, the wave operator turns out to be
easiest to handle, and the source term does not create a
substantial difficulty. The term involving the potential
V turns out to be the most difficult one to handle.
Following Lousto we introduce a staggered grid with
step sizes ∆t = 1
∆r∗ ≡ h, which follows the characteris-
tic lines of the wave operator in Schwarzschild spacetime;
see Fig. 1 for a sketch of a typical grid cell. The basic
idea behind the method is to integrate the wave equation
over a unit cell of the grid, which nicely deals with the
Dirac-δ source term on the right-hand side. To this end,
we introduce the Eddington-Finkelstein null coordinates
v = t + r∗ and u = t − r∗ and use them as integration
variables.
A. Differential operator
Rewriting the wave operator in terms of u and v, we
find −∂2t + ∂2r∗ = −4∂u∂v, which allows us to evaluate
the integral involving the wave operator exactly. We find
−4∂u∂vψ du dv =− 4[ψ(t+ h, r∗) + ψ(t− h, r∗)
− ψ(t, r∗ − h)− ψ(t, r∗ + h)].
(2.1)
B. Source term
If we integrate over a cell traversed by the particle, then
the source term on the right-hand side of the equation
will have a non-zero contribution. Writing the source
term as G(t, r∗)δ(r∗ − r∗0(t)) with
G(t, r∗) = −4πq f
Ȳℓm(π/2, φ0), (2.2)
we find
Gδ(r∗ − r∗0(t)) du dv =−
f0(t)
r0(t)
× Ȳℓm(π/2, φ0(t)) dt,
(2.3)
where t1 and t2 are the times at which the particle enters
and leaves the cell, respectively. While we do not have
an analytic expression for the trajectory of the particle
(except when the particle follows a circular orbit), we can
numerically integrate the first-order ordinary differential
equations that govern the particle’s motion to a precision
that is much higher than that of the partial differential
equation governing ψ. In this sense we treat the integral
over the source term as exact. To evaluate the integral
we adopt a four-point Gauss-Legendre scheme, which has
an error of order h8.
C. Potential term
The most problematic term—from the point of view of
implementing an approximation of sufficiently high or-
der in h—turns out to be the term V ψ in Eq. (1.12).
Since this term does not contain a δ-function, we have to
approximate the double integral
V ψ du dv (2.4)
t0 − h
t0 + h
t0 − 2h
0 − 2h r
0 + 2h
0 − 3h r
0 + 3hr
0 + hr
0 − h
13 34
87 9 10
FIG. 1: Points used to calculate the integral over the potential
term for vacuum cells. Grid points are indicated by blue cir-
cles while red cross-hairs indicate points in between two grid
points. We calculate field values at points that do not lie on
the grid by employing the second-order algorithm described
in [10].
up to terms of order h6 for a generic cell in order to
achieve an overall O(h4) convergence of the scheme.
Here we have to treat cells traversed by the particle
(“sourced” cells) differently from the generic (“vacuum”)
cells. While much of the algorithm can be transferred
from the vacuum cells to the sourced cells, some modifica-
tions are required. We will describe each case separately
in the following subsections.
1. Vacuum case
To implement Lousto’s algorithm to evolve the field
across the vacuum cells, we use a double Simpson rule to
compute the integral Eq. (2.4). We introduce the nota-
g(t, r∗) = V (r∗)ψ(t, r∗) (2.5)
and label our points in the same manner (see Fig. 1) as
in [10]:
g du dv =
[g1 + g2 + g3 + g4 + 4(g12+
g24 + g34 + g13) + 16g0] +O(h
6). (2.6)
Here, for example, g1 is the value of g at the grid point
labeled 1, and g12 is the value of g at the off-grid point
labeled 12, etc. Deviating from Lousto’s algorithm, we
choose to calculate g0 using an expression different from
that derived in [10]. Unlike Lousto’s approach, our ex-
pression exclusively involves points that are within the
past light cone of the current cell. We find
8V4 ψ4 + 8V1 ψ1 + 8V2 ψ2 − 4V6 ψ6 − 4V5 ψ5
+ V10 ψ10 + V7 ψ7 − V9 ψ9 − V8 ψ8
+O(h4). (2.7)
In order to evaluate the term in parentheses in
Eq. (2.6), we again use a variant of the equations given
in [10]. Lousto’s equations (33) and (34),
g13 + g12 =V (r
0 − h/2) (ψ1 + ψ0)
V (r∗0 − h/2)
+O(h4),
(2.8)
g24 + g34 =V (r
0 + h/2) (ψ0 + ψ4)
V (r∗0 + h/2)
+O(h4)
(2.9)
contain isolated occurrences of ψ0, the value of the field
at the central point. Since Eq. (2.7) only allows us to
find g0 = V0ψ0, finding ψ0 would involve a division by
V0, which will be numerically unstable very close to the
event horizon where V0 ≈ 0. Instead we choose to express
the potential term appearing in the square brackets as a
Taylor series around r∗0 . This allows us to eliminate the
isolated occurrences of ψ0, and we find
g13+g12 + g24 + g34 = 2V (r
V (r∗0)
+ V (r∗0 − h/2)ψ1
V (r∗0 − h/2)
+ V (r∗0 + h/2)ψ4
V (r∗0 + h/2)
V (r∗0 − h/2)− 2V (r∗0)
+ V (r∗0 + h/2)
(ψ1 + ψ4) +O(h
4). (2.10)
Because of the
factor in Eq. (2.6), this allows us
to reach the required O(h6) convergence for a generic
vacuum cell. This—given that there is a number of order
N = 1/h2 of such cells—yields the desired overall O(h4)
convergence of the full algorithm, at the end of the N
steps required to finish the simulation.
2. Sourced cells
For vacuum cells, the algorithm described above is
the complete algorithm used to evolve the field forward
in time. For cells traversed by the particle, however,
we have to reconsider the assumptions used in deriving
Eqs. (2.7) and (2.10). When deriving Eq. (2.10) we have
employed the second-order evolution algorithm (see [10]),
in which the single step equation
ψ3 =− ψ2 +
(ψ1 + ψ4) (2.11)
O(h3)O(h4)O(h5) O(h5)O(h4)
traje
tory
t0 + h
t0 − h
t0 − 2h
FIG. 2: Cells affected by the passage of the particle, showing
the reduced order of the single step equation
is accurate only to O(h3) for cells traversed by the parti-
cle. For these cells, therefore, the error term in Eq. (2.10)
is O(h3) instead of O(h4). As there is a number of or-
der N ′ = 1/h of cells that are traversed by the parti-
cle in a simulation run, the overall error—after including
factor in Eq. (2.6)—is of order h4. We can
therefore afford this reduction of the convergence order
in Eq. (2.10)
Equation (2.7), however, is accurate only to O(h) for
cells traversed by the particle. Again taking the
factor into account, this renders the overall algorithm
O(h2). Figure 2 shows the cells affected by the particle’s
traversal and the reduced order of the single step equa-
tion for each cell. Cells whose convergence order is O(h5)
or higher do not need modifications, since there is only a
number N ′ = 1/h of such cells in the simulation. We are
therefore concerned about cells neighboring the particle’s
trajectory and those traversed by the particle.
a. Cells neighboring the particle These cells are not
traversed by the particle, but the particle might have
traversed cells in their past light-cone, which are used in
the calculation of g0 in Eq. (2.7). For these cells, we use
a one-dimensional Taylor expansion of g(t, r∗) within the
current time-slice t = t0,
5V (r∗0 − h)ψ(t0, r∗0 − h)
+ 15V (r∗0 − 3h)ψ(t0, r∗0 − 3h)
− 5V (r∗0 − 5h)ψ(t0, r∗0 − 5h)
+ V (r∗0 − 7h)ψ(t0, r∗0 − 7h)
+O(h4) (2.12)
for the cell on the left-hand side, and
5V (r∗0 + h)ψ(t0, r
0 + h)
+ 15V (r∗0 + 3h)ψ(t0, r
0 + 3h)
− 5V (r∗0 + 5h)ψ(t0, r∗0 + 5h)
+ V (r∗0 + 7h)ψ(t0, r
0 + 7h)
+O(h4) (2.13)
for the cell on the right-hand side, where (t0, r
0) is the
center of the cell traversed by the particle. Both of these
are more accurate than is strictly necessary; we would
t0 − h
t0 + h
0 + 2hr
0 − 2h r
0 − h r
0 + h
(3a) (3b)
(2a) (2b)
FIG. 3: Typical cell traversal of the particle. We split the
domain into sub-parts indicated by the dotted line based on
the time the particle enters (at t1) and leaves (at t2) the cell.
The integral over each sub-part is evaluated using an iterated
two-by-two point Gauss-Legendre rule.
need error terms of order h3 to achieve the desired over-
all O(h4) convergence of the algorithm. Keeping the ex-
tra terms, however, improves the numerical convergence
slightly.
b. Cell traversed by the particle We choose not to
implement a fully explicit algorithm to handle cells tra-
versed by the particle, because this would increase the
complexity of the algorithm by a significant factor. In-
stead we use an iterative approach to evolve the field
using the integrated wave equation
−4(ψ3+ψ2 − ψ1 − ψ4)−
V ψ du dv =
− 8πq
f0(t)
r0(t)
Ȳℓm(π/2, φ0(t)) dt. (2.14)
In this equation the integral involving the source term
can be evaluated to any desired accuracy at the begin-
ning of the iteration, because the motion of the particle
is determined by a simple system of ordinary differential
equations, which are easily integrated with reliable nu-
merical methods. It remains to evaluate the integral over
the potential term, which we do iteratively. Schemati-
cally the method works as follows:
• Make an initial guess for ψ3 using the second-order
algorithm. This guess is correct up to terms of
O(h3).
• Match a second-order piecewise interpolation poly-
nomial to the six points that make up the past light-
cone of the future grid point, including the future
point itself.
• Use this approximation for ψ to numerically calcu-
V ψ du dv,
using two-by-two point Gauss-Legendre rules for
the six sub-parts indicated in Fig. 3.
• Update the future value of the field and repeat the
process until the iteration has converged to a re-
quired degree of accuracy.
trajectory
FIG. 4: Numerical domain evolved during the simulation. We
impose an inner boundary condition close to the black whole
where we can implement it easily to the accuracy of the un-
derlying floating point format. Far away from the black hole,
we evolve the full domain of dependence of the initial data
domain without imposing boundary conditions.
III. INITIAL VALUES AND BOUNDARY
CONDITIONS
As is typical for numerical simulations, we have to pay
careful attention to specifying initial data and appropri-
ate boundary conditions. These aspects of the numerical
method are highly non-trivial problems in full numerical
relativity, but they can be solved or circumvented with
moderate effort in the present work.
A. Initial data
In this work we use a characteristic grid consisting of
points lying on characteristic lines of the wave operator
to evolve ψ forward in time. As such, we need to specify
characteristic initial data on the lines u = u0 and v = v0
shown in Fig. 4. We choose not to worry about specifying
“correct” initial data, but instead arbitrarily choose ψ to
vanish on u = u0 and v = v0:
ψ(u = u0, v) = ψ(u, v = v0) = 0. (3.1)
This is equivalent to adding spurious initial waves in the
form of a homogeneous solution of Eq. (1.12) to the cor-
rect solution. This produces an initial wave burst that
moves away from the particle with the speed of light,
and quickly leaves the numerical domain. Any remain-
ing tails of the spurious initial data decay as t−(2ℓ+2) as
shown in [11] and become negligible after a short time.
We conclude that the influence of the initial-wave con-
tent on the self-force becomes negligible after a time of
the order of the light-crossing time of the particle’s orbit.
B. Boundary conditions
On the analytical side we would like to impose ingoing
boundary conditions at the event horizon r∗ → −∞ and
outgoing boundary conditions at spatial infinity r∗ → ∞,
r∗→−∞
∂uψ =0, lim
∂vψ =0. (3.2)
Because of the finite resources available to a computer
we can only simulate a finite region of the spacetime,
and are faced with the reality of implementing boundary
conditions at finite values of r∗. Two solutions to this
problem present themselves:
1. choose the numerical domain to be the domain of
dependence of the initial data surface. Since the
effect of the boundary condition can only propagate
forward in time with at most the speed of light,
this effectively hides any influence of the boundary.
This is what we choose to do in order to deal with
the outer boundary condition.
2. implement boundary conditions sufficiently “far
out” so that numerically there is no difference be-
tween imposing the boundary condition there or at
infinity. Since the boundary conditions depend on
the vanishing of the potential V (r) appearing in the
wave equation, this will happen once 1−2M/r ≈ 0.
Near the horizon r ≈ 2M(1 + exp(r∗/2M)), so
this will happen—to numerical accuracy—for mod-
estly large (negative) values of r∗ ≈ −73M . We
choose to implement the ingoing waves condition
∂uψℓm = 0 there.
IV. IMPLEMENTATION
Making more precise the ideas developed in the pre-
ceding sections, we implement the following numerical
scheme.
A. Particle motion
Following Darwin [12] we introduce the dimensionless
semi-latus rectum p and the eccentricity e such that for
a bound orbit around a Schwarzschild black hole of mass
1 + e
, r2 =
(4.1)
are the radial positions of the periastron and apastron,
respectively. Energy per unit mass and angular momen-
tum per unit mass are then given by
(p− 2− 2e)(p− 2 + 2e)
p (p− 3− e2)
, L2 =
p− 3− e2
(4.2)
Together with these definitions it is useful to introduce
an orbital parameter χ such that along the trajectory of
the particle,
r(χ) =
1 + e cosχ
, (4.3)
where χ is single-valued along the orbit. We can then
write down first-order differential equations for χ(t) and
the azimuthal angle φ(t) of the particle,
(p− 2− 2e cosχ)(1 + e cosχ)(1 + e cosχ)
(Mp2)
p− 6− 2e cosχ
(p− 2− 2e)(p− 2 + 2e)
, (4.4)
(p− 2− 2e cosχ)(1 + e cosχ)2
p3/2M
(p− 2− 2e)(p− 2 + 2e)
. (4.5)
We use the embedded Runge-Kutta-Fehlberg (4, 5) algo-
rithm provided by the GNU Scientific Library routine
gsl odeiv step rkf45 and an adaptive step-size control
to evolve the position of the particle forward in time.
Intermediate values of the particle’s position are found
using a Hermite interpolation of the nearest available cal-
culated positions.
B. Initial data
We do not specify initial data. The field is set to zero
on the initial characteristic slices, u = u0 and v = v0.
C. Boundary conditions
We adjust the outer boundary of the numerical do-
main at each time-step so that we cover the domain of
dependence of the initial characteristic surfaces and the
particle’s world line. The resulting numerical domain was
already shown in Fig. 4.
Near the event horizon, at r∗ ≈ −73M , we implement
an ingoing-wave boundary condition by imposing
ψ(t+ h, r∗) = ψ(t, r∗ − h). (4.6)
This allows us to drastically reduce the number of cells
in the numerical domain, and consequently the running
time of the simulation.
D. Evolution in vacuum
Cells not traversed by the particle are evolved using
Eqs. (2.1), (2.6) – (2.10). Explicitly written out, we use
ψ3 = −ψ2
(V0 + V1) +
V0 (V0 + V1)
(V0 + V4) +
V0 (V0 + V4)
(g12 + g24 + g34 + g13 + 4g0),
(4.7)
where g0 is given by Eq. (2.7) and the sum g12 + g24 +
g34 + g13 is given by Eq. (2.10).
E. Cells next to the particle
Vacuum cells close to the current position of the parti-
cle require a different approach to calculate g0, since the
cells in their past light cone could have been traversed
by the particle. We use Eqs. (2.12) and (2.13) to find
g0 in this case. Other than this modification, the same
algorithm as for generic vacuum cells is used.
F. Cells traversed by the particle
We evolve cells traversed by the particle using the it-
erative algorithm described in Sec. II C 2. Here
ψ3 =− ψ1 + ψ2
+ ψ4 −
V ψ du dv
f0(t)
r0(t)
Ȳℓm(π/2, π0(t)) dt, (4.8)
where the initial guess for the iterative evolution of∫∫
V ψ du dv is obtained using the second order algo-
rithm of Lousto and Price [13],
ψ3 =− ψ1 +
× [ψ2 + ψ4]
f0(t)
r0(t)
Ȳℓm(π/2, π0(t)) dt. (4.9)
Successive iterations use a four-point Gauss-Legendre
rule to evaluate the integral of V ψ; this requires a second-
order polynomial interpolation of the current field values
as described in Appendix C.
G. Extraction of the field data at the particle
In order to extract the value of the field and its first
derivatives at the position of the particle, we again use
a polynomial interpolation at the points surrounding the
particle’s position. Using a fourth-order polynomial, as
described in Appendix C, we can estimate ψ, ∂tψt, and
∂r∗ψ at the position of the particle up to errors of order
h4. As was briefly mentioned in Sec. II, we would expect
an error term of order h3 for ∂tψt and ∂r∗ψ. The O(h
accuracy we actually achieve by using a fourth-order (in-
stead of a third-order) piecewise polynomial shows up
clearly in a regression plot such as Fig. 7.
H. Regularization of the mode sum
We use the calculated multipole moments ψℓm to con-
struct the multipole moments Φℓm, and first derivatives
∂tΦℓm and ∂rΦℓm, of the scalar field. These, in turn, are
used to calculate the tetrad components Φ(0)ℓm, Φ(+)ℓm,
Φ(−)ℓm, and Φ(3)ℓm of the field gradient according to
Eqs. (1.23)–(1.26) of [4], which are reproduced in Ap-
pendix A. These multipoles then give rise to the multi-
pole coefficients of the retarded field,
Φ(µ)ℓ(t, r, θ, φ) =
Φ(µ)ℓm(t, r)Yℓm(θ, φ), (4.10)
which are subjected to the regularization procedure de-
scribed by Eq. (1.29) of [4],
ΦR(µ)(t, r0, π/2, φ0) = lim
Φ(µ)ℓ(t, r0 +∆, π/2, φ0)
(ℓ+ 1/2)A(µ) +B(µ)
(ℓ + 1/2)
(ℓ− 1/2)(ℓ+ 3/2)
+ · · ·
, (4.11)
using the regularization parameters A(µ), B(µ), C(µ), and
D(µ) tabulated in Appendix B.
Finally we reconstruct the vector components of the
field gradient using Eqs. (1.47)–(1.48) of [4],
ΦRt =
(0), (4.12)
ΦRr =
ΦR(+)e
−iφ0 +ΦR(−)e
, (4.13)
ΦRθ = −r0ΦR(3), (4.14)
ΦRφ = −
ΦR(+)e
−iφ0 − ΦR(−)e
, (4.15)
and calculate the self-force
Fα = qΦ
α . (4.16)
We recall the discussion in Sec. I A concerning the def-
inition of ΦR, its connection to the self-force acting on
the particle, and its regularity at the particle’s position.
V. NUMERICAL TESTS
In this section we present the tests we have performed
to validate our numerical evolution code. First, in order
to check the fourth-order convergence rate of the code,
we perform regression runs with increasing resolution for
both a vacuum test case, where we seeded the evolution
with a Gaussian wave packet, and a case where a particle
is present. As a second test, we compute the regularized
self-force for several different combinations of orbital el-
ements p and e and check that the multipole coefficients
decay with ℓ as expected. This provides a very sensi-
tive check on the overall implementation of the numerical
scheme, as well as the analytical calculations that lead to
the regularization parameters. Finally, we calculate the
self-force for a particle on a circular orbit and show that
it agrees with the results presented in [4, 14].
A. Convergence tests: Vacuum
As a first test of the validity of our numerical code we
estimate the convergence order by removing the particle
and performing regression runs for several resolutions.
We use a Gaussian wave packet as initial data,
ψ(u = u0, v) = exp(−[v − vp]2/[2σ2]), (5.1)
ψ(u, v = v0) = 0, (5.2)
where vp = 75M and σ = 10M , v0 = −u0 = 6M +
2M ln 2, and we extract the field values at r∗ = 20M .
Several such runs were performed, with varying resolu-
tion of 2, 4, 8, 16, and 32 grid points per M . Figure 5
shows ψ(2h)−ψ(h) rescaled by appropriate powers of 2,
so that in the case of fourth-order convergence the curves
would lie on top of each other. As can be seen from the
plots, they do, and the vacuum portion of the code is
indeed fourth-order convergent.
B. Convergence tests: Particle
While the convergence test described in section VA
clearly shows that the desired convergence is achieved
for vacuum evolution, it does not test the parts of the
code that are used in the integration of the inhomoge-
neous wave equation. To test these we perform a second
set of regression runs, this time using a non-zero charge
q. We extract the field at the position of the particle,
thus also testing the implementation of the extraction
algorithm described in section IVG. For this test we
choose the ℓ = 6, m = 4 mode of the field generated
by a particle on a mildly eccentric geodesic orbit with
p = 7, e = 0.3. As shown in Fig. 6 the convergence
is still of fourth order, but the two curves no longer lie
precisely on top of each other at all times. The region
before t ≈ 100M is dominated by the initial wave burst
and therefore does not scale as expected, yielding two
very different curves. In the region 300M . t . 400M
the two curves lie on top of each other, as expected for a
fourth-order convergent algorithm. In the region between
t ≈ 200M and t ≈ 300M , however, the dashed curves
-8.0e-07
-6.0e-07
-4.0e-07
-2.0e-07
0.0e+00
2.0e-07
4.0e-07
6.0e-07
8.0e-07
 0  20  40  60  80  100  120
δ16-8
δ32-16
FIG. 5: Convergence test of the numerical algorithm in the
vacuum case. We show differences between simulations using
different step sizes h = 0.5M (ψ2), h = 0.25M (ψ4), h =
0.125M (ψ8), h = 0.0625M (ψ16), and h = 0.03125M (ψ32).
Displayed are the rescaled differences δ4−2 = ψ4 −ψ2, δ8−4 =
24(ψ8 − ψ4), δ16−8 = 4
4(ψ8 − ψ4), and δ32−16 = 8
4(ψ8 − ψ4)
for the real part of the ℓ = 2, m = 2 mode at r∗ ≈ 20M .
The maximum value of the field itself is of the order of 0.1,
so that the errors in the field values are roughly five orders
of magnitude smaller than the field values themselves. We
can see that the convergence is in fact of fourth-order, as the
curves lie nearly on top of each other, with only the lowest
resolution curve δ4−2 deviating slightly.
have slightly smaller amplitudes than the solid one, indi-
cating an order of convergence different from (but close
to) four.
To explain this behavior we have to examine the terms
that contribute significantly to the error in the simula-
tion. The numerical error is almost completely domi-
nated by that of the approximation of the potential term∫∫
V ψ du dv in the integrated wave equation. For vac-
uum cells the error in this approximation scales as h6,
where h is the step size. For cells traversed by the parti-
cle, on the other hand, the approximation error depends
also on the difference t2−t1 of the times at which the par-
ticle enters and leaves the cell. This difference is bounded
by h but does not necessarily scale as h. For example, if
a particle enters a cell at its very left, then scaling h by 1
would not change t2 − t1 at all, thus leading to a scaling
behavior that differs from expectation.
To investigate this further we conducted test runs
of the simulation for a particle on a circular orbit at
r = 6M . In order to observe the expected scaling behav-
ior, we have to make sure that the particle passes through
the tips of the cell it traverses. When this is the case, then
t2 − t1 ≡ h and a plot similar to the one shown in Fig. 6
shows the proper scaling behavior. As a further test we
artificially reduced the convergence order of the vacuum
algorithm to two by implementing the second-order algo-
rithm described in [10]. By keeping the algorithm that
deals with sourced cells unchanged, we reduced the rela-
-4e-05
-2e-05
 2e-05
 4e-05
 100  200  300  400  500
δ16-8
FIG. 6: Convergence test of the numerical algorithm in the
sourced case. We show differences between simulations using
different step sizes of 4 (ψ4), 8 (ψ8), 16 (ψ16), and 32 (ψ32)
cells per M . Displayed are the rescaled differences δ8−4 =
ψ8 −ψ4, etc. (see caption of Fig. 5 for definitions) of the field
values at the position of the particle for a simulation with
ℓ = 6, m = 4 and p = 7, e = 0.3. We see that the convergence
is approximately fourth-order.
-4e-05
-2e-05
 2e-05
 4e-05
 100  200  300  400  500
δ32-16
δ16-8
FIG. 7: Convergence test of the numerical algorithm in the
sourced case. We show differences between ∂rΦ for simula-
tions using different step sizes of 4 (Φr,4), 8 (Φr,8), 16 (Φr,16),
and 32 (Φr,32) cells per M . Displayed are the rescaled differ-
ences δ8−4 = Φr,8 − Φr,4 etc. of the values at the position of
the particle for a simulation with ℓ = 6, m = 4 and p = 7,
e = 0.3. Although there is much noise caused by the piece-
wise polynomials used to extract the data, we can see that
the convergence is approximately fourth-order.
tive impact on the numerical error. This, too, allows us
to recover the expected (second-order) convergence. Fig-
ures 8 and 9 illustrate the effects of the measures taken
to control the convergence behavior.
-3.000e-05
-2.000e-05
-1.000e-05
0.000e+00
1.000e-05
2.000e-05
3.000e-05
 200  210  220  230  240  250  260  270  280
δ64-32
δ32-16
δ16-8
-3.000e-05
-2.000e-05
-1.000e-05
0.000e+00
1.000e-05
2.000e-05
3.000e-05
 200  210  220  230  240  250  260  270  280
FIG. 8: Behavior of convergence tests for a particle in circular
orbit at r = 6M . We show differences between simulations of
the ℓ = 2, m = 2 multipole moment using different step sizes
of 2 (ψ2), 4 (ψ4), 8 (ψ8), 16 (ψ16), 32 (ψ32) and 64 (ψ64) cells
per M . Displayed are the real part of the rescaled differences
δ4−2 = (ψ4 −ψ2) etc. of the field values at the position of the
particle, defined as in Fig. 5. The values have been rescaled
so that—for fourth order convergence—the curves should all
coincide. The upper panel corresponds to a set of simulations
where the particle traverses the cells away from their tips.
The curves do not coincide perfectly with each other, seem-
ingly indicating a failure of the convergence. The lower panel
was obtained in a simulation where the particle was carefully
positioned so as to pass through the tips of each cell it tra-
verses. This set of simulations passes the convergence test
more convincingly.
C. High-ℓ behavior of the multipole coefficients
Inspection of Eq. (4.11) reveals that a plot of Φ(µ)ℓ as
a function of ℓ (for a selected value of t) should display
a linear growth in ℓ for large ℓ. Removing the A(µ) term
should produce a constant curve, removing the B(µ) term
(given that C(µ) = 0) should produce a curve that decays
as ℓ−2, and finally, removing the D(µ) term should pro-
duce a curve that decays as ℓ−4. It is a powerful test of
the numerical methods to check whether these expecta-
tions are borne out by the numerical data. Fig. 10 plots
the remainders as obtained from our numerical simula-
tion, demonstrating the expected behavior. It displays,
on a logarithmic scale, the absolute value of ReΦR(+)ℓ, the
real part of the (+) component of the self-force. The orbit
is eccentric (p = 7.2, e = 0.5), and all components of the
self-force require regularization. The first curve (in trian-
gles) shows the unregularized multipole coefficients that
increase linearly in ℓ, as confirmed by fitting a straight
line to the data. The second curve (in squares) shows par-
tially regularized coefficients, obtained after the removal
of (ℓ + 1/2)A(µ); this clearly approaches a constant for
large values of ℓ. The curve made up of diamonds shows
the behavior after removal of B(µ); because C(µ) = 0, it
decays as ℓ−2, a behavior that is confirmed by a fit to
-8.000e-05
-4.000e-05
0.000e+00
4.000e-05
8.000e-05
 200  210  220  230  240  250  260  270  280
-8.000e-05
-4.000e-05
0.000e+00
4.000e-05
8.000e-05
 200  210  220  230  240  250  260  270  280
-8.000e-05
-4.000e-05
0.000e+00
4.000e-05
8.000e-05
 200  210  220  230  240  250  260  270  280
δ64-32
δ32-16
δ16-8
FIG. 9: Behavior of convergence tests for a particle in circular
orbit at r = 6M . We show differences between simulations
of the ℓ = 2, m = 2 multipole moment using different step
sizes of 8 (ψ8), 16 (ψ16), 32 (ψ32), and 64 (ψ64) cells per M .
Displayed are the real part of the rescaled differences δ16−8 =
ψ16−ψ8 etc. of the field values at the position of the particle,
defined as in Fig. 5. The values have been rescaled so that—
for second order convergence—the curves should all coincide.
The upper two panels correspond to simulations where the
second order algorithm was used throughout. For the topmost
one, care was taken to ensure that the particle passes through
the tip of each cell it traverses, while in the middle one no
such precaution was taken. Clearly the curves in the middle
panel do not coincide with each other, indicating a failure
of the second-order convergence of the code. The lower panel
was obtained in a simulation using the mixed-order algorithm
described in the text. While the curves still do not coincide
precisely, the observed behavior is much closer to the expected
one than for the purely second order algorithm.
the ℓ ≥ 5 part of the curve. Finally, after removal of
D(µ)/[(ℓ − 12 ) (ℓ +
)] the terms of the sum decrease in
magnitude as ℓ−4 for large values of ℓ, as derived in [15].
Each one of the last two curves would result in a con-
verging sum, but the convergence is much faster after
subtracting the D(µ) terms. We thereby gain more than
2 orders of magnitude in the accuracy of the estimated
Figure 10 provides a sensitive test of the implemen-
tation of both the numerical and analytical parts of the
calculation. Small mistakes in either one will cause the
difference in Eq. (4.11) to have a vastly different behav-
 0  2  4  6  8  10  12  14
ReΦ(+)
ReΦ(+)-A
ReΦ(+)-A-B
ReΦ(+)-A-B-D
FIG. 10: Multipole coefficients of the dimensionless self-force
ReΦR(+) for a particle on an eccentric orbit (p = 7.2, e =
0.5). The coefficients are extracted at t = 500M along the
trajectory shown in Fig. 12. The plots show several stages of
the regularization procedure, with a closer description of the
curves to be found in the text.
 0  2  4  6  8  10  12  14  16  18  20
FIG. 11: Multipole coefficients of ΦR(0) for a particle on a circu-
lar orbit. Note that ΦR(0)ℓ is linked to Φ
t via Φ
The multipole coefficients decay exponentially with ℓ until
ℓ ≈ 16, at which point numerical errors start to dominate.
D. Self-force on a circular orbit
For the case of a circular orbit, the regularization pa-
rameters A(0), B(0), and D(0) all vanish identically, so
that the (0) (or alternatively the t) component of the
self-force does not require regularization. Figure 11 thus
shows only one curve, with the magnitude of the multi-
pole coefficients decaying exponentially with increasing
As a final test, in Table I we compare our result for the
self-force on a particle in a circular orbit at r = 6M to
those obtained in [4, 14] using a frequency-domain code.
For a circular orbit, a calculation in the frequency domain
TABLE I: Results for the self-force on a scalar particle with
scalar charge q on a circular orbit at r0 = 6M . The
first column lists the results as calculated in this work us-
ing time-domain numerical methods, while the second and
third columns list the results as calculated in [4, 14] using
frequency-domain methods. For the t and φ components the
number of digits is limited by numerical roundoff error. For
the r component the number of digits is limited by the trun-
cation error of the sum of multipole coefficients.
This work: Previous work: Diaz-Rivera
time-domain frequency-domain [4] et. al. [14]
ΦRt 3.60339 × 10
−4 3.60907254 × 10−4
1.6767 × 10−4 1.67730 × 10−4 1.6772834 × 10−4
ΦRφ −5.30424 × 10
−3 −5.30423170 × 10−3
is more efficient, and we expect the results of [4, 14] to
be much more accurate than our own results. This fact
is reflected in the number of regularization coefficients
we can reliably extract from the numerical data, before
being limited by the accuracy of the numerical method:
the frequency-domain calculation found usable multipole
coefficients up to ℓ = 20, whereas our data for ΦR(0)ℓ is
dominated by noise by the time ℓ reaches 16. Figure 11
shows this behavior.
E. Accuracy of the numerical method
Several figures of merit can be used to estimate the
accuracy of numerical values for the self-force.
An estimate for the truncation error arising from cut-
ting short the summation in Eq. (4.11) at some ℓmax can
be calculated by considering the behavior of the remain-
ing terms for large ℓ. Detweiler et. al. [15] showed that
the remaining terms scale as ℓ−4 for large ℓ. They find
the functional form of the terms to be
EP3/2
(2ℓ− 3)(2ℓ− 1)(2ℓ+ 3)(2ℓ+ 5)
, (5.3)
where P3/2 = 36
2. We fit a function of this form to the
tail end of a plot of the multipole coefficients to find the
coefficient E in Eq. (5.3). Extrapolating to ℓ → ∞ we
find that the truncation error is
ℓ=ℓmax
[Eq. (5.3)] (5.4)
2Eℓmax
(2ℓmax + 3)(2ℓmax + 1)(2ℓmax − 1)(2ℓmax − 3)
(5.5)
where ℓmax is the value at which we cut the summation
short. For all but the special case of the (0) component
for a circular orbit, for which all regularization parame-
ters vanish identically, we use this approach to calculate
an estimate for the truncation error.
A second source of error lies in the numerical calcula-
tion of the retarded solution to the wave equation. This
error depends on the step size h used to evolve the field
forward in time. For a numerical scheme of a given con-
vergence order, we can estimate this discretization error
by extrapolating the differences of simulations using dif-
ferent step sizes down to h = 0. This is what was done
in the graphs shown in Sec. VB.
We display results for mildly eccentric orbits. A high
eccentricity causes ∂rΦ (displayed in Fig. 7) to be plagued
by high frequency noise produced by effects similar to
those described in Sec. VB. This makes it impossible to
reliably estimate the discretization error for these orbits.
We do not expect this to be very different from the errors
for mildly eccentric orbits.
Finally we compare our final results for the self-force
Fα to “reference values”. For circular orbits, frequency-
domain calculations are much more accurate than our
time-domain computations. We thus compare our results
to the results obtained in [4]. Table II lists typical values
for the various errors listed above.
error estimation mildly eccentric orbit
truncation error (M
Φ(+)) ≈ 2× 10
discretization error (M
∂rΦℓm) ≈ 10
comparison with reference values circular orbit
Ft 0.2%
Fr 0.04%
Fφ 2× 10
TABLE II: Estimated values for the various errors in the com-
ponents of the self-force as described in the text. We show
the truncation and discretization errors for a mildly eccentric
orbit and the total error for a circular orbit. The truncation
error is calculated using a plot similar to the one shown in
Fig. 16. The discretization error is estimated using a plot
similar to that in Fig. 7 for the ℓ = 2, m = 2 mode, and the
total error is estimated as the difference between our values
and those of [4]. We use p = 7.2 , e = 0.5 for the mildly
eccentric orbit. Note that we use the tetrad component Φ(+)
for the truncation error and the vector component ∂rΦ for
the discretization error. Both are related by the translation
table Eqs. (A6) – (A9), we expect corresponding errors to be
comparable for Φ(+) and ∂rΦ.
VI. SAMPLE RESULTS
In this section we describe some results of our numer-
ical calculation.
A. Mildly eccentric orbit
We choose a particle on an eccentric orbit with p = 7.2,
e = 0.5 which starts at r = pM/(1−e2), halfway between
 15  10  5  0  5  10
trajectory for p=7.2, e=0.5
FIG. 12: Trajectory of a particle with p = 7.2, e = 0.5. The
cross-hair indicates the point where the data for Fig. 10 was
extracted.
-0.014
-0.012
-0.01
-0.008
-0.006
-0.004
-0.002
 0.002
 100  200  300  400  500  600  700  800  900  1000
time/M
FIG. 13: Regularized dimensionless self-force M
and M
Fφ on a particle on an eccentric orbit with p = 7.2,
e = 0.5.
periastron and apastron. The field is evolved for 1000M
with a resolution of 16 grid points per M , both in the t
and r∗ directions, for ℓ = 0. Higher values of ℓ (and thus
m) require a corresponding increase in the number of
grid points used to achieve the same fractional accuracy.
Multipole coefficients for 0 ≤ ℓ ≤ 15 are calculated and
used to reconstruct the regularized self-force Fα along
the geodesic. Figure 13 shows the result of the calcula-
tion. For the choice of parameters used to calculate the
force shown in Fig. 13, the error bars corresponding to
the truncation error (which are already much larger than
than the discretization error) would be of the order of
the line thickness and have not been drawn.
Already for this small eccentricity, we see that the self-
force is most important when the particle is closest to the
black hole (ie. for 200M . t . 400M and 600M . t .
 60  50  40  30  20  10  0  10  20  30  40  50
trajectory for p=7.8001, e=0.9
FIG. 14: Trajectory of a particle on a zoom-whirl orbit with
p = 7.8001, e = 0.9. The cross-hairs indicate the positions
where the data shown in Fig. 16 and 17 was extracted.
800M); the self-force acting on the particle is very small
once the particle has moved away to r ≈ 15M .
B. Zoom-whirl orbit
Highly eccentric orbits are of most interest as sources
of gravitational radiation. For nearly parabolic orbits
with e . 1 and p & 6+2e, a particle revolves around the
black hole a number of times, moving on a nearly circu-
lar trajectory close to the event horizon (“whirl phase”),
before moving away from the black hole (“zoom phase”).
During the whirl phase the particle is in the strong field
region of the black hole, emitting copious amounts of
radiation. Figures 14 and 15 show the trajectory of a
particle and the force on such an orbit with p = 7.8001,
e = 0.9. Even more so than for the mildly eccentric
orbit discussed in Sec. VIA, the self-force (and thus the
amount of radiation produced) is much larger while the
particle is close to the black hole than when it zooms out.
Defining energy E per unit mass and angular momen-
tum L per unit mass in the usual way,
E = −
uα, L =
uα, (6.1)
and following eg. the treatment of Wald [16], Ap-
pendix C, it is easy to see that the rates of change Ė
and L̇ (per unit proper time) are directly related to com-
ponents of the acceleration aα (and therefore force) ex-
perienced by the particle via
Ė = −at, L̇ = aφ. (6.2)
The self-force shown in Fig. 15 therefore confirms our
näıve expectation that the self-force should decrease both
the energy and angular momentum of the particle as ra-
diation is emitted.
-0.025
-0.02
-0.015
-0.01
-0.005
 0.005
 500  1000  1500  2000  2500  3000
time/M
-0.002
-0.001
 0.001
 0.002
 0.003
 0.004
 2000 2050 2100 2150 2200
FIG. 15: Self-force acting on a particle. Shown is the dimen-
sionless self-force M
Fr and
Fφ on a zoom-whirl
orbit with p = 7.8001, e = 0.9. The inset shows a magni-
fied view of the self-force when the particle is about to enter
the whirl phase. No error bars showing an estimate error are
shown, as the errors shown eg. in Table II are to small to
show up on the graph. Notice that the self-force is essentially
zero during the zoom phase 500M . t . 2000M and reaches
a constant value very quickly after the particle enters into the
whirl phase.
It is instructive to have a closer look at the force acting
on the particle when it is within the zoom phase, and also
when it is moving around the black hole on the nearly cir-
cular orbit of the whirl phase. In Fig. 16 and Fig. 17 we
show plots of Φ(0)ℓ vs. ℓ after the removal of the A(µ),
B(µ), and D(µ) terms. While the particle is still zooming
in toward the black hole, Φ(0)ℓ behaves exactly as for the
mildly eccentric orbit described in Sec. VIA over the full
range of ℓ plotted; ie. the magnitude of each term scales
as ℓ0, ℓ−2 and ℓ−4, after removal of the A(µ), B(µ), and
D(µ) terms respectively. Close to the black hole, on the
other hand, the particle moves along a nearly circular tra-
jectory. If the orbit were perfectly circular for all times,
ie. ṙ ≡ 0, then the (0) component would not require reg-
ularization at all, and the multipole coefficients would
decay exponentially, resulting in a straight line on the
semi-logarithmic plot shown in Fig. 17. As the real orbit
is not precisely circular, curves eventually deviate from a
straight line. Removal of the A(µ) term is required almost
immediately (beginning with ℓ ≈ 3), while the D(µ) term
starts to become important only after ℓ ≈ 11. This shows
that there is a smooth transition from the self-force on a
circular orbit, which does not require regularization for
the t and φ components, to that of a generic orbit, for
which all components of the self-force require regulariza-
tion.
 0  2  4  6  8  10  12  14
Φ(0)-A
Φ(0)-A-B
Φ(0)-A-B-D
FIG. 16: Multipole coefficients of M
ReΦR(0) for a particle on
a zoom-whirl orbit (p = 7.8001, e = 0.9). The coefficients are
extracted at t = 2000M as the particle is about to enter the
whirl phase. As ṙ is non-zero, all components of the self-force
require regularization and we see that the dependence of the
multipole coefficients on ℓ is as predicted by Eq. 1.9. After the
removal of the regularization parameters A(µ), B(µ), and D(µ)
the remainder is proportional to ℓ0, ℓ−2 and ℓ−4 respectively.
 0  2  4  6  8  10  12  14
Φ(0)-A
Φ(0)-A-B
Φ(0)-A-B-D
FIG. 17: Multipole coefficients of ReΦR(0) for a particle on
a zoom-whirl orbit (p = 7.8001, e = 0.9). The coefficients
are extracted at t = 2150M while the particle is in the whirl
phase. The orbit is nearly circular at this time, causing the
dependence on ℓ after removal of the regularization parame-
ters to approximate that of a true circular orbit.
Acknowledgments
We thank Eric Poisson and Eran Rosenthal for useful
discussions and suggestions. This work was supported by
the Natural Sciences and Engineering Council of Canada.
APPENDIX A: TRANSLATION TABLES
We quote the results of [4] for the translation table be-
tween the modes Φℓm and the tetrad components Φ(µ)ℓm
with respect to the pseudo-Cartesian basis
eα(0) =
, 0, 0, 0
, (A1)
eα(1) =
f sin θ cosφ,
cos θ cosφ,− sinφ
r sin θ
, (A2)
eα(2) =
f sin θ sinφ,
cos θ sinφ,
r sin θ
, (A3)
eα(3) =
f cos θ,−1
sin θ, 0
, (A4)
and the complex combinations eα
:= eα
± ieα
eα(±) =
f sin θe±iφ,
cos θe±iφ,
±ie±iφ
r sin θ
. (A5)
With these, the spherical-harmonic modes Φ(µ)ℓm(t, r)
are given in terms of Φℓm(t, r) by
Φ(0)ℓm =
Φℓm, (A6)
Φ(+)ℓm =−
(ℓ+m− 1)(ℓ +m)
(2ℓ− 1)(2ℓ+ 1)
− ℓ− 1
Φℓ−1,m−1
(ℓ−m+ 1)(ℓ −m+ 2)
(2ℓ+ 1)(2ℓ+ 3)
Φℓ+1,m−1, (A7)
Φ(−)ℓm =
(ℓ −m− 1)(ℓ−m)
(2ℓ− 1)(2ℓ+ 1)
− ℓ− 1
Φℓ−1,m+1
(ℓ+m+ 1)(ℓ +m+ 2)
(2ℓ+ 1)(2ℓ+ 3)
Φℓ+1,m+1, (A8)
Φ(3)ℓm =
(ℓ −m)(ℓ+m)
(2ℓ− 1)(2ℓ+ 1)
Φℓ−1,m
(ℓ−m+ 1)(ℓ +m+ 1)
(2ℓ+ 1)(2ℓ+ 3)
Φℓ+1,m. (A9)
APPENDIX B: REGULARIZATION
PARAMETERS
For completeness we list the regularization parameters
as calculated in [4]. Quantities bearing a subscript “0”
are evaluated at the particle’s position.
A(0) =
0 + L
sign(∆), (B1)
A(+) = −eiφ0
0 + L
sign(∆), (B2)
A(3) = 0, (B3)
where f0 := 1 − 2M/r0 and sign(∆) is equal to +1 if
∆ > 0 and to −1 if ∆ < 0. We have, in addition, A(−) =
Ā(+), A(1) = Re[A(+)], and A(2) = Im[A(+)].
We also use
B(0) = −
Er0ṙ0√
0 + L
2)3/2
E + Er0ṙ0
0 + L
2)3/2
B(+) = e
Bc(+) − iB
, (B5)
Bc(+) =
0 + L
2)3/2
r20 + L
0 + L
2)3/2
f0 − 1
r20 + L
K, (B6)
Bs(+) = −
f0)ṙ0
r20 + L
E + (2−
f0)ṙ0
r20 + L
B(3) = 0. (B8)
In addition, B(−) = B̄(+), B(1) = Re[B(+)] =
cosφ0 + B
sinφ0, and B(2) = Im[B(+)] =
Bc(+) sinφ0 −B
(+) cosφ0.
Here, the rescaled elliptic integrals E and K are defined
E := 2
∫ π/2
(1− k sin2 ψ)1/2 dψ = F
; 1; k
K := 2
∫ π/2
(1− k sin2 ψ)−1/2 dψ = F
; 1; k
(B10)
in which k := L2/(r20 + L
We also use
C(µ) = 0 (B11)
D(0) = −
Er30(r
0 − L2)ṙ30
0 + L
2)7/2
E(r70 + 30Mr
0 − 7L2r50 + 114ML2r40 + 104ML4r20 + 36ML6)ṙ0
16r40
0 + L
2)5/2
Er30(5r
0 − 3L2)ṙ30
0 + L
2)7/2
E(r50 + 16Mr
0 − 3L2r30 + 42ML2r20 + 18ML4)ṙ0
16r20
0 + L
2)5/2
K, (B12)
D(+) = e
Dc(+) − iD
, (B13)
Dc(+) =
r30(r
0 − L2)ṙ40
0 + L
2)7/2
− r0ṙ
4(r20 + L
2)3/2
(3r70 + 6Mr
0 − L2r50 + 31ML2r40 + 26ML4r20 + 9ML6)ṙ20
0 + L
2)5/2
(3r70 + 8Mr
0 + L
2r50 + 26ML
2r40 + 22ML
4r20 + 8ML
16r60(r
0 + L
2)3/2
0 + 2Mr
0 + 4ML
r20 + L
0 − 3L2)ṙ40
0 + L
2)7/2
8(r20 + L
2)3/2
− (7r
0 + 12Mr
0 − L2r30 + 46ML2r20 + 18ML4)ṙ20
16r20
0 + L
2)5/2
− (7r
0 + 6Mr
0 + 6L
2r30 + 12ML
2r20 + 4ML
16r40(r
0 + L
2)3/2
r20 + L
K, (B14)
Ds(+) =
r20(r
0 − 7L2)(
f0 − 2)ṙ30
0 + L
2)5/2
− (2r
0 +Mr
0 + 5L
2r50 + 10ML
2r40 + 29ML
4r20 + 14ML
6)ṙ0
8r50L(r
0 + L
2)3/2
(r50 −Mr40 + 4L2r30 − 5ML2r20 + 2ML4)ṙ0
4r30L
0 + L
2)3/2
r20(r
0 − 3L2)(
f0 − 2)ṙ30
0 + L
2)5/2
(4r50 + 2Mr
0 + 7L
2r30 + 10ML
2r20 + 14ML
4)ṙ0
16r30L(r
0 + L
2)3/2
− (2r
0 − 2Mr20 + 5L2r0 − 8ML2)ṙ0
0 + L
2)3/2
K, (B15)
D(3) = 0. (B16)
And finally, D(−) = D̄(+), D(1) = Re[D(+)] =
cosφ0 + D
sinφ0, and D(2) = Im[D(+)] =
Dc(+) sinφ0 −D
(+) cosφ0.
APPENDIX C: PIECEWISE POLYNOMIALS
In two places in the numerical simulation we introduce
piecewise polynomials to approximate the scalar field ψℓm
across the world line, where it is continuous but not dif-
ferentiable. By a piecewise polynomial we mean a poly-
nomial of the form
p(t, r∗) =


n,m=0
unvm if r∗(u, v) > r∗0
n,m=0
unvm if r∗(u, v) < r∗0
, (C1)
where u = t−r∗, v = t+r∗ are characteristic coordinates,
r∗0 is the position of the particle at the time t(u, v), and
N is the order of the polynomial, which for our purposes
is N = 4 or less. The two sets of coefficients cnm and
c′nm are not independent of each other, but are linked
via jump conditions that can be derived from the wave
equation [Eq. (1.12)]. To do so, we rewrite the wave
equation in the characteristic coordinates u and v and
reintroduce the integral over the world line on the right-
hand side,
−4∂u∂vψ − V ψ =
Ŝ(τ)δ(u − up)δ(v − vp) dτ , (C2)
where Ŝ(τ) = −8πq Ȳℓm(π/2,φp(τ))
rp(τ)
is the source term and
quantities bearing a subscript p are evaluated on the
world line at proper time τ .
Here and in the following we use the notation
[∂nu∂
v ψ] = lim
[∂nu∂
v ψ(t0, r
0 + ǫ)− ∂nu∂mv ψ(t0, r∗0 − ǫ)]
to denote the jump in ∂nu∂
v ψ across the world line. First,
we notice that the source term does not contain any
derivatives of the Dirac δ-function, causing the solution
ψ to be continuous. This means that the zeroth-order
jump vanishes: [ψ] = 0. Our task is then to find the re-
maining jump conditions at a point (t0, r
0) for n,m ≤ 4.
Alternatively, instead of crossing the world line along a
line t = t0 = const we can also choose to cross along
lines of u = u0 = const or v = v0 = const, noting that
for a line of constant v the coordinate u runs from u0+ ǫ
to u0 − ǫ to cross from the left to the right of the world
line. Figure 18 provides a clearer description of the paths
taken.
(u0 − ǫ, v0)
(u0, v0 + ǫ)(u0 + ǫ, v0)
(u0, v0 − ǫ)
(t0, r
0) = (u0, v0)
FIG. 18: Paths taken in the calculation of the jump condi-
tions. (u0, v0) denotes an arbitrary but fixed point along the
world line γ. The wave equation is integrated along the lines
of constant u or v indicated in the sketch. Note that in order
to move from the domain on the left to the domain on the
right, u has to run from u0 + ǫ to u0 − ǫ. Where appropriate
we label quantities connected to the domain on the left by a
subscript “−” and quantities connected to the domain on the
right by “+”.
In order to find the jump [∂uψ] we integrate the wave
equation along the line u = u0 from v0 − ǫ to v0 + ǫ
∫ v0+ǫ
∂u∂vψdv −
∫ v0+ǫ
V ψdv =
Ŝ(τ)δ(u0 − up)
∫ v0+ǫ
δ(v − vp)dv dτ ,
which, after involving
∫ v0+ǫ
δ(v − vp)dv = θ(vp − v0 +
ǫ)θ(v0 − vp + ǫ) and δ(g(x)) = δ(x− x0)/ |g′(x0)|, yields
[∂uψ] = −
E − ṙ0
Ŝ(τ0), (C5)
where the overdot denotes differentiation with respect to
proper time τ .
Similarly, after first taking a derivative of the wave
equation with respect to v and integrating from u0+ ǫ to
u0 − ǫ, we obtain
∫ u0−ǫ
vψdu−
∫ u0−ǫ
V ψdu =
Ŝ(τ)
∫ u0−ǫ
δ(u− up)du δ′(v0 − vp)dτ .
We find
E + ṙ0
E + ṙp
Ŝ(τ)
|τ=τ0
. (C7)
Systematically repeating this procedure we find expres-
sions for the jumps in all the derivatives that are purely
in the u or v direction. Table III lists these results. Jump
[ψ] =0
[∂uψ] =−
bS(τ0), [∂vψ] =
bS(τ0)
bS(τ )
|τ=τ0
bS(τ )
|τ=τ0
V ξ0ξ̄
0 [∂uψ]−
bS(τ )
|τ=τ0
V ξ̄0ξ
0 [∂vψ] +
bS(τ )
|τ=τ0
|τ=τ0
+ 3ξ0ξ̄
0 ∂uV + ξ
0 ∂vV
[∂uψ] +
bS(τ )
|τ=τ0
|τ=τ0
+ 3ξ̄0ξ
0 ∂vV + ξ̄
0 ∂uV
[∂vψ]−
bS(τ )
|τ=τ0
TABLE III: Jump conditions for the derivatives purely in the
u or v directions. ṙ and r̈ are the particle’s radial velocity and
acceleration, respectively. They are obtained from the equa-
tion of motion for the particle. ξ̄ := E−ṙ
and ξ := E+ṙ
introduced for notational convenience. Quantities bearing a
subscript p are evaluated on the particle’s world line, while
quantities bearing a subscript 0 are evaluated at the parti-
cle’s current position. Derivatives of V with respect to either
u or v are evaluated as ∂uV = −
f∂rV and ∂vV =
f∂rV ,
respectively.
conditions for derivatives involving both u and v are ob-
tained directly from the wave equation [Eq. (C2)]. We
see that
[∂u∂vψ] = 0, (C8)
and taking an additional derivative with respect to u on
both sides reveals that
∂2u∂vψ
V [∂uψ] . (C9)
Systematically repeating this procedure we can find jump
conditions for each of the mixed derivatives by evaluating
∂n+1u ∂
[∂nu∂
v (V ψ)] , (C10)
where n,m ≥ 0 and derivatives of V with respect to
either u or v are evaluated as ∂uV = − 12f∂rV and ∂vV =
f∂rV , respectively.
The results of Table III and Eq. (C10) allow us to
express the coefficients of the left-hand polynomial in
Eq. (C1) in terms of the jump conditions and the co-
efficients of the right-hand side:
c′nm = cnm − [∂nu∂mv ψ] . (C11)
For N = 4 this leaves us with 25 unknown coefficients
cnm which can be uniquely determined by demanding
that the polynomial match the value of the field on the
25 grid points surrounding the particle. When we are
interested in integrating the polynomial, as in the case of
the potential term in the fourth-order algorithm, we do
not need all these terms. Instead, in order to calculate
e.g. the integral
V ψ du dv up to terms of order h5,
as is needed to achieve overall O(h4) convergence, it is
sufficient to include only terms such that n+m ≤ 2, thus
reducing the number of unknown coefficients to 6. In this
case Eq. (C1) becomes
p(t, r∗) =
m+n≤2
unvm if r∗(u, v) > r∗0
m+n≤2
unvm if r∗(u, v) < r∗0
. (C12)
The six coefficients can then be determined by matching
the polynomial to the field values at the six grid points
which lie within the past light cone of the grid point
whose field value we want to calculate.
[1] The LISA web site is located at
http://lisa.jpl.nasa.gov/.
[2] L. Barack and A. Ori, Phys. Rev. D 61, 061502 (2000),
gr-qc/9912010.
[3] C. Lousto, Class. Quantum Grav. 22, S543 (2005).
[4] R. Haas and E. Poisson, Phys. Rev. D 74,
044009 (pages 29) (2006), gr-qc/0605077, URL
http://link.aps.org/abstract/PRD/v74/e044009 .
[5] B. S. DeWitt and R. W. Brehme, Annals of Physics 9,
220 (1960).
[6] Y. Mino, M. Sasaki, and T. Tanaka, Phys. Rev. D 55,
3457 (1997), gr-qc/9606018.
[7] T. C. Quinn and R. M. Wald, Phys. Rev. D 56, 3381
(1997), gr-qc/9610053.
[8] T. C. Quinn, Phys. Rev. D 62, 064029 (2000), gr-
qc/0005030.
[9] S. Detweiler and B. F. Whiting, Phys. Rev. D 67, 024025
(2003), gr-qc/0202086.
[10] C. O. Lousto, Class. Quant. Grav. 22, S543 (2005), gr-
qc/0503001.
[11] R. H. Price, Phys. Rev. D 5, 2419 (1972).
[12] C. G. Darwin, Proc. R. Soc. A 249, 180 (1959).
http://lisa.jpl.nasa.gov/
http://link.aps.org/abstract/PRD/v74/e044009
[13] C. O. Lousto and R. H. Price, Phys. Rev. D 56, 6439
(1997), gr-qc/9705071.
[14] L. M. Diaz-Rivera, E. Messaritaki, B. F. Whit-
ing, and S. Detweiler, Physical Review D (Par-
ticles, Fields, Gravitation, and Cosmology) 70,
124018 (pages 14) (2004), gr-qc/0410011, URL
http://link.aps.org/abstract/PRD/v70/e124018.
[15] S. Detweiler, E. Messaritaki, and B. F. Whiting, Phys.
Rev. D 67, 104016 (2003), gr-qc/0205079.
[16] R. M. Wald, General relativity (University of Chicago
Press, Chicago, 1984), ISBN 0226870324.
http://link.aps.org/abstract/PRD/v70/e124018

  We calculate the self-force acting on a particle with scalar charge moving on
a generic geodesic around a Schwarzschild black hole. This calculation requires
an accurate computation of the retarded scalar field produced by the moving
charge; this is done numerically with the help of a fourth-order convergent
finite-difference scheme formulated in the time domain. The calculation also
requires a regularization procedure, because the retarded field is singular on
the particle's world line; this is handled mode-by-mode via the mode-sum
regularization scheme first introduced by Barack and Ori. This paper presents
the numerical method, various numerical tests, and a sample of results for
mildly eccentric orbits as well as ``zoom-whirl'' orbits.

Scalar self-force on eccentric geodesics in Schwarzschild spacetime: A time-domain
computation
Roland Haas
Department of Physics, University of Guelph, Guelph, Ontario, Canada N1G 2W1
(Dated: April 3, 2007)
We calculate the self-force acting on a particle with scalar charge moving on a generic geodesic
around a Schwarzschild black hole. This calculation requires an accurate computation of the retarded
scalar field produced by the moving charge; this is done numerically with the help of a fourth-order
convergent finite-difference scheme formulated in the time domain. The calculation also requires
a regularization procedure, because the retarded field is singular on the particle’s world line; this
is handled mode-by-mode via the mode-sum regularization scheme first introduced by Barack and
Ori. This paper presents the numerical method, various numerical tests, and a sample of results for
mildly eccentric orbits as well as “zoom-whirl” orbits.
PACS numbers: 04.25.-g, 04.40.-b, 41.60.-m, 45.50.-j, 02.60.Cb, 02.70.Bf
I. INTRODUCTION
The inspiral and capture of solar-mass compact objects
by supermassive black holes is one of the most promis-
ing and interesting sources of gravitational radiation to
be detected by the future space-based gravitational-wave
antenna LISA [1]. For these extreme mass-ratio inspirals,
one can treat the compact object as a point mass and de-
scribe its influence on the spacetime perturbatively. Go-
ing beyond the test mass limit, its motion is no longer
along a geodesic of the unperturbed spacetime of the cen-
tral black hole; it is a geodesic of the perturbed space-
time created by the presence of the moving body. When
viewed from the unperturbed spacetime, the small body
is said to move under the influence of its gravitational
self-force. The self-force induces radiative losses of energy
and angular momentum, which will eventually drive the
object into the black hole. To describe the motion of the
body, including its inspiral toward the black hole, we seek
to evaluate the self-force and calculate its effect on the
motion. One way of doing this uses the mode-sum reg-
ularization procedure introduced by Barack and Ori [2].
(For a comprehensive introduction of the problem, see
the special issue of Classical and Quantum Gravity [3].)
In this paper, in an effort to build expertise to calculate
the gravitational self-force, we retreat to the technically
simpler problem of a point particle of mass m endowed
with a scalar charge q orbiting a Schwarzschild black hole
of mass M . Following up on a previous paper [4], we
implement the numerical part of the regularization pro-
cedure for generic orbits with a time-domain integration
of the scalar-wave equation.
A. The problem
Our goal is to calculate the regularized self-force acting
on a scalar point charge in orbit around a Schwarzschild
black hole. In analogy with the gravitational case, where
in a first-order (in m/M) perturbative calculation the
particle moves on a geodesic of the background space-
time, we take the orbit of the particle to be a geodesic and
calculate the self-force as a vector field on this geodesic.
We start by writing the Schwarzschild metric using the
tortoise coordinate r∗ = r + 2M ln
ds2 = f
−dt2 + dr∗2
+ r2dΩ2, (1.1)
where f =
1− 2M
, dΩ2 =
dθ2 + sin2 θdφ2
is the
metric on a two-sphere, and t, r, θ and φ are the usual
Schwarzschild coordinates. Our task is to solve the scalar
wave equation
gαβ∇α∇βΦ(x) = −4πµ(x), (1.2)
µ(x) = q
δ4(x, z(τ))dτ , (1.3)
where ∇α is the covariant derivative compatible with the
metric gαβ , Φ(x) is the scalar field created by a scalar
charge q which moves along a world line γ : τ 7→ z(τ)
parametrized by proper time τ . The source term µ(x)
appearing on the right-hand side is written in terms of a
scalarized four-dimensional Dirac δ-function δ4(x, x
′) :=
δ(x0 − x′0)δ(x1 − x′1)δ(x2 − x′2)δ(x3 − x′3)/
− det(gαβ).
Because of the singularity in the source term, the re-
tarded solution to Eq. (1.2) is singular on the world line,
and the näıve expression for the self-force,
Fα(τ) = q∇αΦ(z(τ)), (1.4)
must be regularized. Following DeWitt and Brehme [5],
Mino, Sasaki, Tanaka [6], Quinn and Wald [7], Quinn [8]
carried out this regularization for the electromagnetic,
scalar and gravitational radiation reaction. In later work,
Detweiler and Whiting [9] introduced a very useful de-
composition of the retarded solution of Eq. (1.2) in terms
of a singular part ΦS and a regular remainder ΦR:
Φ = ΦS +ΦR. (1.5)
ΦR is regular and differentiable at the position of the par-
ticle, satisfies the homogeneous wave equation associated
http://arxiv.org/abs/0704.0797v2
with Eq. (1.2), and is solely responsible for the self-force
acting on the particle. ΦS , on the other hand, satisfies
Eq. (1.2), is just as singular at the particle’s position as
the retarded solution, and produces no force on the par-
ticle. Rearranging Eq. (1.5) and differentiating once, we
can write the regularized self-force as
Fα := q∇αΦR = q
∇αΦ−∇αΦS
. (1.6)
In a previous paper [4], we described our implemen-
tation of the regularization procedure to find a mode-
sum representation of ∇αΦS along a generic geodesic of
the Schwarzschild spacetime. Schematically, we intro-
duce a tetrad eα
and decompose the tetrad components
Φ(µ) := e
∇αΦ of the field gradient in terms of ordinary
scalar spherical harmonics Yℓm:
Φ(µ)(t, r, θ, φ) =
Φℓm(µ)(t, r)Yℓm(θ, φ). (1.7)
Each mode Φℓm
(t, r) is finite at the position of the par-
ticle, but their sum diverges on the world line. In [4], we
derive analytic expressions for the mode-sum decompo-
sition of ΦS(µ),
ΦS(µ) =q
ΦS(µ),ℓ (1.8)
ΦS(µ),ℓ = A(µ)
+B(µ) +
(ℓ− 1
)(ℓ + 3
+ · · · , (1.9)
where the coefficients A(µ), B(µ), C(µ), and D(µ) are in-
dependent of ℓ; they are listed in Appendix B for conve-
nience.
As each mode of Φ is finite, it is straightforward to
compute the modes of the retarded solution using nu-
merical methods, and we will describe how this was done
in Sec. IV. We use the numerical solutions in Eq. (1.6)
to calculate the regularized self-force, regularizing mode-
by-mode:
ΦR(µ) =
Φ(µ),ℓ − ΦS(µ),ℓ
, (1.10)
where Φ(µ),ℓ :=
Yℓm (no summation over ℓ im-
plied).
For numerical purposes it is convenient to define ψℓm
Φ(x) =
ℓm, (1.11)
where Yℓm are the usual scalar spherical harmonics. Af-
ter substituting in Eq. (1.2), this yields a reduced wave
equation for the multipole moments ψℓm:
−∂2tψℓm + ∂2r∗ψℓm − Vℓψℓm =
− 4πq
Ȳℓm(π/2, φ0)δ(r
∗ − r∗0), (1.12)
where
Vℓ = f
ℓ (ℓ+ 1)
. (1.13)
An overbar denotes complex conjugation, E = −ut is the
particle’s conserved energy per unit mass, and uα = dz
is its four velocity. Quantities bearing a subscript “0” are
evaluated at the particle’s position; they are functions of
τ that are obtained by solving the geodesic equation
uβ∇βuα = 0 (1.14)
in the background spacetime. Without loss of general-
ity, we have confined the motion of the particle to the
equatorial plane θ = π
Once we have numerically solved Eq. (1.12), we ex-
tract numerical estimates for ψℓm, ∂tψℓm and ∂r∗ψℓm,
which can then be used to find Φℓm, ∂tΦℓm and ∂rΦℓm.
These—together with the translation table displayed in
Eqs. (1.23)–(1.26) of [4], reproduced in Appendix A—
allow us to find the tetrad components Φ(µ)ℓm with re-
spect to the tetrad defined by Eqs. (1.18)–(1.21) of [4].
Eventually we regularize the multipole coefficients
Φ(µ)ℓ =
Φ(µ)ℓm(t0, r0)Y
ℓm(π/2, φ0) (1.15)
using Eq. (1.10); this involves the regularization param-
eters listed in Eqs. (1.30)–(1.45) of [4], which are repro-
duced in Appendix B.
B. Organization of this paper
In Sec. II we introduce the main ideas behind the
discretization scheme used in the numerical simulation.
Sec. III describes the choices we make in order to handle
the problems of specifying initial data and proper bound-
ary conditions. The next section—Sec. IV—provides de-
tails on the concrete implementation of the ideas put
forth in Secs. II and III. In Sec. V we describe the tests
we performed in order to validate our implementation of
the numerical method. Sec. VI finally presents sample
results for a small number of representative simulations.
C. Future work
This work, which deals with a scalar charge moving in
the Schwarzschild spacetime, is not intended to produce
physically or astrophysically interesting results. Instead,
its goal is to help us evaluate the merits of several strate-
gies that could be used to tackle the more interesting (and
difficult) problems of electromagnetism and gravity.
One future project we are currently exploring is to ap-
ply the formalism developed so far to the electromagnetic
self-force acting on an electric charge. Beyond the tech-
nical complication of having to deal with a vector field
instead of a single scalar quantity, we are also faced with
the reality of having to impose a gauge (in our case: the
Lorenz gauge) and to eliminate (or at least control) gauge
violations in the numerical simulation. The first step,
namely, the calculation of the regularization parameters
A(µ), B(µ), C(µ), and D(µ) for the self-force, is currently
underway. Also underway is the calculation of the regu-
larization parameters for he gravitational self-force.
Another project is the implementation of a scheme to
use the calculated self-force to update the orbital pa-
rameters of a particle on its inspiral toward the black
hole. The standard proposed approach to this problem
in the past has been to calculate the self-force on a set
of geodesics which are momentarily tangent to the par-
ticle’s trajectory. The self-force calculated in this way
is then used to update the orbital elements. This “after
the fact” calculation of the motion requires one to build
(in advance) a large database of self-force values for the
anticipated set of orbital parameters that the particle’s
trajectory will assume during its inspiral. Alternatively,
and conceptually more simply, the self-force could be cal-
culated self-consistently along the real, accelerated tra-
jectory. Such an approach requires changes in the expres-
sions of the regularization parameters, which so far have
been derived only for geodesic orbits. We are currently
investigating the merits of such an approach.
II. NUMERICAL METHOD
In this section we describe the algorithm used to inte-
grate the reduced wave equation [Eq. (1.12)] numerically.
For the most part we use the fourth-order algorithm in-
troduced by Lousto [10], with some modifications to suit
our needs. We choose to implement a fourth-order con-
vergent code because second-order convergence for the
potential Φ, while much easier to achieve, would guaran-
tee only first-order convergence for ∇αΦ, the quantity in
which we are ultimately interested. With a fourth-order
convergent code we can expect to achieve third-order con-
vergence for ∇αΦ, which is required for an accurate es-
timation of the self-force. Numerical experiments, how-
ever, show that in practice we do achieve fourth-order
convergence for the derivatives of Φ, a fortunate outcome
that we exploit but cannot explain.
From now on, we will suppress the subscripts ℓ and m
on Vℓ and ψℓm for convenience of notation. The wave
equation consists of three parts: the wave-operator term
(∂2r∗ − ∂2t )ψ and the potential term V ψ on the left-hand
side, and the source term on the right-hand side of the
equation. Of these, the wave operator turns out to be
easiest to handle, and the source term does not create a
substantial difficulty. The term involving the potential
V turns out to be the most difficult one to handle.
Following Lousto we introduce a staggered grid with
step sizes ∆t = 1
∆r∗ ≡ h, which follows the characteris-
tic lines of the wave operator in Schwarzschild spacetime;
see Fig. 1 for a sketch of a typical grid cell. The basic
idea behind the method is to integrate the wave equation
over a unit cell of the grid, which nicely deals with the
Dirac-δ source term on the right-hand side. To this end,
we introduce the Eddington-Finkelstein null coordinates
v = t + r∗ and u = t − r∗ and use them as integration
variables.
A. Differential operator
Rewriting the wave operator in terms of u and v, we
find −∂2t + ∂2r∗ = −4∂u∂v, which allows us to evaluate
the integral involving the wave operator exactly. We find
−4∂u∂vψ du dv =− 4[ψ(t+ h, r∗) + ψ(t− h, r∗)
− ψ(t, r∗ − h)− ψ(t, r∗ + h)].
(2.1)
B. Source term
If we integrate over a cell traversed by the particle, then
the source term on the right-hand side of the equation
will have a non-zero contribution. Writing the source
term as G(t, r∗)δ(r∗ − r∗0(t)) with
G(t, r∗) = −4πq f
Ȳℓm(π/2, φ0), (2.2)
we find
Gδ(r∗ − r∗0(t)) du dv =−
f0(t)
r0(t)
× Ȳℓm(π/2, φ0(t)) dt,
(2.3)
where t1 and t2 are the times at which the particle enters
and leaves the cell, respectively. While we do not have
an analytic expression for the trajectory of the particle
(except when the particle follows a circular orbit), we can
numerically integrate the first-order ordinary differential
equations that govern the particle’s motion to a precision
that is much higher than that of the partial differential
equation governing ψ. In this sense we treat the integral
over the source term as exact. To evaluate the integral
we adopt a four-point Gauss-Legendre scheme, which has
an error of order h8.
C. Potential term
The most problematic term—from the point of view of
implementing an approximation of sufficiently high or-
der in h—turns out to be the term V ψ in Eq. (1.12).
Since this term does not contain a δ-function, we have to
approximate the double integral
V ψ du dv (2.4)
t0 − h
t0 + h
t0 − 2h
0 − 2h r
0 + 2h
0 − 3h r
0 + 3hr
0 + hr
0 − h
13 34
87 9 10
FIG. 1: Points used to calculate the integral over the potential
term for vacuum cells. Grid points are indicated by blue cir-
cles while red cross-hairs indicate points in between two grid
points. We calculate field values at points that do not lie on
the grid by employing the second-order algorithm described
in [10].
up to terms of order h6 for a generic cell in order to
achieve an overall O(h4) convergence of the scheme.
Here we have to treat cells traversed by the particle
(“sourced” cells) differently from the generic (“vacuum”)
cells. While much of the algorithm can be transferred
from the vacuum cells to the sourced cells, some modifica-
tions are required. We will describe each case separately
in the following subsections.
1. Vacuum case
To implement Lousto’s algorithm to evolve the field
across the vacuum cells, we use a double Simpson rule to
compute the integral Eq. (2.4). We introduce the nota-
g(t, r∗) = V (r∗)ψ(t, r∗) (2.5)
and label our points in the same manner (see Fig. 1) as
in [10]:
g du dv =
[g1 + g2 + g3 + g4 + 4(g12+
g24 + g34 + g13) + 16g0] +O(h
6). (2.6)
Here, for example, g1 is the value of g at the grid point
labeled 1, and g12 is the value of g at the off-grid point
labeled 12, etc. Deviating from Lousto’s algorithm, we
choose to calculate g0 using an expression different from
that derived in [10]. Unlike Lousto’s approach, our ex-
pression exclusively involves points that are within the
past light cone of the current cell. We find
8V4 ψ4 + 8V1 ψ1 + 8V2 ψ2 − 4V6 ψ6 − 4V5 ψ5
+ V10 ψ10 + V7 ψ7 − V9 ψ9 − V8 ψ8
+O(h4). (2.7)
In order to evaluate the term in parentheses in
Eq. (2.6), we again use a variant of the equations given
in [10]. Lousto’s equations (33) and (34),
g13 + g12 =V (r
0 − h/2) (ψ1 + ψ0)
V (r∗0 − h/2)
+O(h4),
(2.8)
g24 + g34 =V (r
0 + h/2) (ψ0 + ψ4)
V (r∗0 + h/2)
+O(h4)
(2.9)
contain isolated occurrences of ψ0, the value of the field
at the central point. Since Eq. (2.7) only allows us to
find g0 = V0ψ0, finding ψ0 would involve a division by
V0, which will be numerically unstable very close to the
event horizon where V0 ≈ 0. Instead we choose to express
the potential term appearing in the square brackets as a
Taylor series around r∗0 . This allows us to eliminate the
isolated occurrences of ψ0, and we find
g13+g12 + g24 + g34 = 2V (r
V (r∗0)
+ V (r∗0 − h/2)ψ1
V (r∗0 − h/2)
+ V (r∗0 + h/2)ψ4
V (r∗0 + h/2)
V (r∗0 − h/2)− 2V (r∗0)
+ V (r∗0 + h/2)
(ψ1 + ψ4) +O(h
4). (2.10)
Because of the
factor in Eq. (2.6), this allows us
to reach the required O(h6) convergence for a generic
vacuum cell. This—given that there is a number of order
N = 1/h2 of such cells—yields the desired overall O(h4)
convergence of the full algorithm, at the end of the N
steps required to finish the simulation.
2. Sourced cells
For vacuum cells, the algorithm described above is
the complete algorithm used to evolve the field forward
in time. For cells traversed by the particle, however,
we have to reconsider the assumptions used in deriving
Eqs. (2.7) and (2.10). When deriving Eq. (2.10) we have
employed the second-order evolution algorithm (see [10]),
in which the single step equation
ψ3 =− ψ2 +
(ψ1 + ψ4) (2.11)
O(h3)O(h4)O(h5) O(h5)O(h4)
traje
tory
t0 + h
t0 − h
t0 − 2h
FIG. 2: Cells affected by the passage of the particle, showing
the reduced order of the single step equation
is accurate only to O(h3) for cells traversed by the parti-
cle. For these cells, therefore, the error term in Eq. (2.10)
is O(h3) instead of O(h4). As there is a number of or-
der N ′ = 1/h of cells that are traversed by the parti-
cle in a simulation run, the overall error—after including
factor in Eq. (2.6)—is of order h4. We can
therefore afford this reduction of the convergence order
in Eq. (2.10)
Equation (2.7), however, is accurate only to O(h) for
cells traversed by the particle. Again taking the
factor into account, this renders the overall algorithm
O(h2). Figure 2 shows the cells affected by the particle’s
traversal and the reduced order of the single step equa-
tion for each cell. Cells whose convergence order is O(h5)
or higher do not need modifications, since there is only a
number N ′ = 1/h of such cells in the simulation. We are
therefore concerned about cells neighboring the particle’s
trajectory and those traversed by the particle.
a. Cells neighboring the particle These cells are not
traversed by the particle, but the particle might have
traversed cells in their past light-cone, which are used in
the calculation of g0 in Eq. (2.7). For these cells, we use
a one-dimensional Taylor expansion of g(t, r∗) within the
current time-slice t = t0,
5V (r∗0 − h)ψ(t0, r∗0 − h)
+ 15V (r∗0 − 3h)ψ(t0, r∗0 − 3h)
− 5V (r∗0 − 5h)ψ(t0, r∗0 − 5h)
+ V (r∗0 − 7h)ψ(t0, r∗0 − 7h)
+O(h4) (2.12)
for the cell on the left-hand side, and
5V (r∗0 + h)ψ(t0, r
0 + h)
+ 15V (r∗0 + 3h)ψ(t0, r
0 + 3h)
− 5V (r∗0 + 5h)ψ(t0, r∗0 + 5h)
+ V (r∗0 + 7h)ψ(t0, r
0 + 7h)
+O(h4) (2.13)
for the cell on the right-hand side, where (t0, r
0) is the
center of the cell traversed by the particle. Both of these
are more accurate than is strictly necessary; we would
t0 − h
t0 + h
0 + 2hr
0 − 2h r
0 − h r
0 + h
(3a) (3b)
(2a) (2b)
FIG. 3: Typical cell traversal of the particle. We split the
domain into sub-parts indicated by the dotted line based on
the time the particle enters (at t1) and leaves (at t2) the cell.
The integral over each sub-part is evaluated using an iterated
two-by-two point Gauss-Legendre rule.
need error terms of order h3 to achieve the desired over-
all O(h4) convergence of the algorithm. Keeping the ex-
tra terms, however, improves the numerical convergence
slightly.
b. Cell traversed by the particle We choose not to
implement a fully explicit algorithm to handle cells tra-
versed by the particle, because this would increase the
complexity of the algorithm by a significant factor. In-
stead we use an iterative approach to evolve the field
using the integrated wave equation
−4(ψ3+ψ2 − ψ1 − ψ4)−
V ψ du dv =
− 8πq
f0(t)
r0(t)
Ȳℓm(π/2, φ0(t)) dt. (2.14)
In this equation the integral involving the source term
can be evaluated to any desired accuracy at the begin-
ning of the iteration, because the motion of the particle
is determined by a simple system of ordinary differential
equations, which are easily integrated with reliable nu-
merical methods. It remains to evaluate the integral over
the potential term, which we do iteratively. Schemati-
cally the method works as follows:
• Make an initial guess for ψ3 using the second-order
algorithm. This guess is correct up to terms of
O(h3).
• Match a second-order piecewise interpolation poly-
nomial to the six points that make up the past light-
cone of the future grid point, including the future
point itself.
• Use this approximation for ψ to numerically calcu-
V ψ du dv,
using two-by-two point Gauss-Legendre rules for
the six sub-parts indicated in Fig. 3.
• Update the future value of the field and repeat the
process until the iteration has converged to a re-
quired degree of accuracy.
trajectory
FIG. 4: Numerical domain evolved during the simulation. We
impose an inner boundary condition close to the black whole
where we can implement it easily to the accuracy of the un-
derlying floating point format. Far away from the black hole,
we evolve the full domain of dependence of the initial data
domain without imposing boundary conditions.
III. INITIAL VALUES AND BOUNDARY
CONDITIONS
As is typical for numerical simulations, we have to pay
careful attention to specifying initial data and appropri-
ate boundary conditions. These aspects of the numerical
method are highly non-trivial problems in full numerical
relativity, but they can be solved or circumvented with
moderate effort in the present work.
A. Initial data
In this work we use a characteristic grid consisting of
points lying on characteristic lines of the wave operator
to evolve ψ forward in time. As such, we need to specify
characteristic initial data on the lines u = u0 and v = v0
shown in Fig. 4. We choose not to worry about specifying
“correct” initial data, but instead arbitrarily choose ψ to
vanish on u = u0 and v = v0:
ψ(u = u0, v) = ψ(u, v = v0) = 0. (3.1)
This is equivalent to adding spurious initial waves in the
form of a homogeneous solution of Eq. (1.12) to the cor-
rect solution. This produces an initial wave burst that
moves away from the particle with the speed of light,
and quickly leaves the numerical domain. Any remain-
ing tails of the spurious initial data decay as t−(2ℓ+2) as
shown in [11] and become negligible after a short time.
We conclude that the influence of the initial-wave con-
tent on the self-force becomes negligible after a time of
the order of the light-crossing time of the particle’s orbit.
B. Boundary conditions
On the analytical side we would like to impose ingoing
boundary conditions at the event horizon r∗ → −∞ and
outgoing boundary conditions at spatial infinity r∗ → ∞,
r∗→−∞
∂uψ =0, lim
∂vψ =0. (3.2)
Because of the finite resources available to a computer
we can only simulate a finite region of the spacetime,
and are faced with the reality of implementing boundary
conditions at finite values of r∗. Two solutions to this
problem present themselves:
1. choose the numerical domain to be the domain of
dependence of the initial data surface. Since the
effect of the boundary condition can only propagate
forward in time with at most the speed of light,
this effectively hides any influence of the boundary.
This is what we choose to do in order to deal with
the outer boundary condition.
2. implement boundary conditions sufficiently “far
out” so that numerically there is no difference be-
tween imposing the boundary condition there or at
infinity. Since the boundary conditions depend on
the vanishing of the potential V (r) appearing in the
wave equation, this will happen once 1−2M/r ≈ 0.
Near the horizon r ≈ 2M(1 + exp(r∗/2M)), so
this will happen—to numerical accuracy—for mod-
estly large (negative) values of r∗ ≈ −73M . We
choose to implement the ingoing waves condition
∂uψℓm = 0 there.
IV. IMPLEMENTATION
Making more precise the ideas developed in the pre-
ceding sections, we implement the following numerical
scheme.
A. Particle motion
Following Darwin [12] we introduce the dimensionless
semi-latus rectum p and the eccentricity e such that for
a bound orbit around a Schwarzschild black hole of mass
1 + e
, r2 =
(4.1)
are the radial positions of the periastron and apastron,
respectively. Energy per unit mass and angular momen-
tum per unit mass are then given by
(p− 2− 2e)(p− 2 + 2e)
p (p− 3− e2)
, L2 =
p− 3− e2
(4.2)
Together with these definitions it is useful to introduce
an orbital parameter χ such that along the trajectory of
the particle,
r(χ) =
1 + e cosχ
, (4.3)
where χ is single-valued along the orbit. We can then
write down first-order differential equations for χ(t) and
the azimuthal angle φ(t) of the particle,
(p− 2− 2e cosχ)(1 + e cosχ)(1 + e cosχ)
(Mp2)
p− 6− 2e cosχ
(p− 2− 2e)(p− 2 + 2e)
, (4.4)
(p− 2− 2e cosχ)(1 + e cosχ)2
p3/2M
(p− 2− 2e)(p− 2 + 2e)
. (4.5)
We use the embedded Runge-Kutta-Fehlberg (4, 5) algo-
rithm provided by the GNU Scientific Library routine
gsl odeiv step rkf45 and an adaptive step-size control
to evolve the position of the particle forward in time.
Intermediate values of the particle’s position are found
using a Hermite interpolation of the nearest available cal-
culated positions.
B. Initial data
We do not specify initial data. The field is set to zero
on the initial characteristic slices, u = u0 and v = v0.
C. Boundary conditions
We adjust the outer boundary of the numerical do-
main at each time-step so that we cover the domain of
dependence of the initial characteristic surfaces and the
particle’s world line. The resulting numerical domain was
already shown in Fig. 4.
Near the event horizon, at r∗ ≈ −73M , we implement
an ingoing-wave boundary condition by imposing
ψ(t+ h, r∗) = ψ(t, r∗ − h). (4.6)
This allows us to drastically reduce the number of cells
in the numerical domain, and consequently the running
time of the simulation.
D. Evolution in vacuum
Cells not traversed by the particle are evolved using
Eqs. (2.1), (2.6) – (2.10). Explicitly written out, we use
ψ3 = −ψ2
(V0 + V1) +
V0 (V0 + V1)
(V0 + V4) +
V0 (V0 + V4)
(g12 + g24 + g34 + g13 + 4g0),
(4.7)
where g0 is given by Eq. (2.7) and the sum g12 + g24 +
g34 + g13 is given by Eq. (2.10).
E. Cells next to the particle
Vacuum cells close to the current position of the parti-
cle require a different approach to calculate g0, since the
cells in their past light cone could have been traversed
by the particle. We use Eqs. (2.12) and (2.13) to find
g0 in this case. Other than this modification, the same
algorithm as for generic vacuum cells is used.
F. Cells traversed by the particle
We evolve cells traversed by the particle using the it-
erative algorithm described in Sec. II C 2. Here
ψ3 =− ψ1 + ψ2
+ ψ4 −
V ψ du dv
f0(t)
r0(t)
Ȳℓm(π/2, π0(t)) dt, (4.8)
where the initial guess for the iterative evolution of∫∫
V ψ du dv is obtained using the second order algo-
rithm of Lousto and Price [13],
ψ3 =− ψ1 +
× [ψ2 + ψ4]
f0(t)
r0(t)
Ȳℓm(π/2, π0(t)) dt. (4.9)
Successive iterations use a four-point Gauss-Legendre
rule to evaluate the integral of V ψ; this requires a second-
order polynomial interpolation of the current field values
as described in Appendix C.
G. Extraction of the field data at the particle
In order to extract the value of the field and its first
derivatives at the position of the particle, we again use
a polynomial interpolation at the points surrounding the
particle’s position. Using a fourth-order polynomial, as
described in Appendix C, we can estimate ψ, ∂tψt, and
∂r∗ψ at the position of the particle up to errors of order
h4. As was briefly mentioned in Sec. II, we would expect
an error term of order h3 for ∂tψt and ∂r∗ψ. The O(h
accuracy we actually achieve by using a fourth-order (in-
stead of a third-order) piecewise polynomial shows up
clearly in a regression plot such as Fig. 7.
H. Regularization of the mode sum
We use the calculated multipole moments ψℓm to con-
struct the multipole moments Φℓm, and first derivatives
∂tΦℓm and ∂rΦℓm, of the scalar field. These, in turn, are
used to calculate the tetrad components Φ(0)ℓm, Φ(+)ℓm,
Φ(−)ℓm, and Φ(3)ℓm of the field gradient according to
Eqs. (1.23)–(1.26) of [4], which are reproduced in Ap-
pendix A. These multipoles then give rise to the multi-
pole coefficients of the retarded field,
Φ(µ)ℓ(t, r, θ, φ) =
Φ(µ)ℓm(t, r)Yℓm(θ, φ), (4.10)
which are subjected to the regularization procedure de-
scribed by Eq. (1.29) of [4],
ΦR(µ)(t, r0, π/2, φ0) = lim
Φ(µ)ℓ(t, r0 +∆, π/2, φ0)
(ℓ+ 1/2)A(µ) +B(µ)
(ℓ + 1/2)
(ℓ− 1/2)(ℓ+ 3/2)
+ · · ·
, (4.11)
using the regularization parameters A(µ), B(µ), C(µ), and
D(µ) tabulated in Appendix B.
Finally we reconstruct the vector components of the
field gradient using Eqs. (1.47)–(1.48) of [4],
ΦRt =
(0), (4.12)
ΦRr =
ΦR(+)e
−iφ0 +ΦR(−)e
, (4.13)
ΦRθ = −r0ΦR(3), (4.14)
ΦRφ = −
ΦR(+)e
−iφ0 − ΦR(−)e
, (4.15)
and calculate the self-force
Fα = qΦ
α . (4.16)
We recall the discussion in Sec. I A concerning the def-
inition of ΦR, its connection to the self-force acting on
the particle, and its regularity at the particle’s position.
V. NUMERICAL TESTS
In this section we present the tests we have performed
to validate our numerical evolution code. First, in order
to check the fourth-order convergence rate of the code,
we perform regression runs with increasing resolution for
both a vacuum test case, where we seeded the evolution
with a Gaussian wave packet, and a case where a particle
is present. As a second test, we compute the regularized
self-force for several different combinations of orbital el-
ements p and e and check that the multipole coefficients
decay with ℓ as expected. This provides a very sensi-
tive check on the overall implementation of the numerical
scheme, as well as the analytical calculations that lead to
the regularization parameters. Finally, we calculate the
self-force for a particle on a circular orbit and show that
it agrees with the results presented in [4, 14].
A. Convergence tests: Vacuum
As a first test of the validity of our numerical code we
estimate the convergence order by removing the particle
and performing regression runs for several resolutions.
We use a Gaussian wave packet as initial data,
ψ(u = u0, v) = exp(−[v − vp]2/[2σ2]), (5.1)
ψ(u, v = v0) = 0, (5.2)
where vp = 75M and σ = 10M , v0 = −u0 = 6M +
2M ln 2, and we extract the field values at r∗ = 20M .
Several such runs were performed, with varying resolu-
tion of 2, 4, 8, 16, and 32 grid points per M . Figure 5
shows ψ(2h)−ψ(h) rescaled by appropriate powers of 2,
so that in the case of fourth-order convergence the curves
would lie on top of each other. As can be seen from the
plots, they do, and the vacuum portion of the code is
indeed fourth-order convergent.
B. Convergence tests: Particle
While the convergence test described in section VA
clearly shows that the desired convergence is achieved
for vacuum evolution, it does not test the parts of the
code that are used in the integration of the inhomoge-
neous wave equation. To test these we perform a second
set of regression runs, this time using a non-zero charge
q. We extract the field at the position of the particle,
thus also testing the implementation of the extraction
algorithm described in section IVG. For this test we
choose the ℓ = 6, m = 4 mode of the field generated
by a particle on a mildly eccentric geodesic orbit with
p = 7, e = 0.3. As shown in Fig. 6 the convergence
is still of fourth order, but the two curves no longer lie
precisely on top of each other at all times. The region
before t ≈ 100M is dominated by the initial wave burst
and therefore does not scale as expected, yielding two
very different curves. In the region 300M . t . 400M
the two curves lie on top of each other, as expected for a
fourth-order convergent algorithm. In the region between
t ≈ 200M and t ≈ 300M , however, the dashed curves
-8.0e-07
-6.0e-07
-4.0e-07
-2.0e-07
0.0e+00
2.0e-07
4.0e-07
6.0e-07
8.0e-07
 0  20  40  60  80  100  120
δ16-8
δ32-16
FIG. 5: Convergence test of the numerical algorithm in the
vacuum case. We show differences between simulations using
different step sizes h = 0.5M (ψ2), h = 0.25M (ψ4), h =
0.125M (ψ8), h = 0.0625M (ψ16), and h = 0.03125M (ψ32).
Displayed are the rescaled differences δ4−2 = ψ4 −ψ2, δ8−4 =
24(ψ8 − ψ4), δ16−8 = 4
4(ψ8 − ψ4), and δ32−16 = 8
4(ψ8 − ψ4)
for the real part of the ℓ = 2, m = 2 mode at r∗ ≈ 20M .
The maximum value of the field itself is of the order of 0.1,
so that the errors in the field values are roughly five orders
of magnitude smaller than the field values themselves. We
can see that the convergence is in fact of fourth-order, as the
curves lie nearly on top of each other, with only the lowest
resolution curve δ4−2 deviating slightly.
have slightly smaller amplitudes than the solid one, indi-
cating an order of convergence different from (but close
to) four.
To explain this behavior we have to examine the terms
that contribute significantly to the error in the simula-
tion. The numerical error is almost completely domi-
nated by that of the approximation of the potential term∫∫
V ψ du dv in the integrated wave equation. For vac-
uum cells the error in this approximation scales as h6,
where h is the step size. For cells traversed by the parti-
cle, on the other hand, the approximation error depends
also on the difference t2−t1 of the times at which the par-
ticle enters and leaves the cell. This difference is bounded
by h but does not necessarily scale as h. For example, if
a particle enters a cell at its very left, then scaling h by 1
would not change t2 − t1 at all, thus leading to a scaling
behavior that differs from expectation.
To investigate this further we conducted test runs
of the simulation for a particle on a circular orbit at
r = 6M . In order to observe the expected scaling behav-
ior, we have to make sure that the particle passes through
the tips of the cell it traverses. When this is the case, then
t2 − t1 ≡ h and a plot similar to the one shown in Fig. 6
shows the proper scaling behavior. As a further test we
artificially reduced the convergence order of the vacuum
algorithm to two by implementing the second-order algo-
rithm described in [10]. By keeping the algorithm that
deals with sourced cells unchanged, we reduced the rela-
-4e-05
-2e-05
 2e-05
 4e-05
 100  200  300  400  500
δ16-8
FIG. 6: Convergence test of the numerical algorithm in the
sourced case. We show differences between simulations using
different step sizes of 4 (ψ4), 8 (ψ8), 16 (ψ16), and 32 (ψ32)
cells per M . Displayed are the rescaled differences δ8−4 =
ψ8 −ψ4, etc. (see caption of Fig. 5 for definitions) of the field
values at the position of the particle for a simulation with
ℓ = 6, m = 4 and p = 7, e = 0.3. We see that the convergence
is approximately fourth-order.
-4e-05
-2e-05
 2e-05
 4e-05
 100  200  300  400  500
δ32-16
δ16-8
FIG. 7: Convergence test of the numerical algorithm in the
sourced case. We show differences between ∂rΦ for simula-
tions using different step sizes of 4 (Φr,4), 8 (Φr,8), 16 (Φr,16),
and 32 (Φr,32) cells per M . Displayed are the rescaled differ-
ences δ8−4 = Φr,8 − Φr,4 etc. of the values at the position of
the particle for a simulation with ℓ = 6, m = 4 and p = 7,
e = 0.3. Although there is much noise caused by the piece-
wise polynomials used to extract the data, we can see that
the convergence is approximately fourth-order.
tive impact on the numerical error. This, too, allows us
to recover the expected (second-order) convergence. Fig-
ures 8 and 9 illustrate the effects of the measures taken
to control the convergence behavior.
-3.000e-05
-2.000e-05
-1.000e-05
0.000e+00
1.000e-05
2.000e-05
3.000e-05
 200  210  220  230  240  250  260  270  280
δ64-32
δ32-16
δ16-8
-3.000e-05
-2.000e-05
-1.000e-05
0.000e+00
1.000e-05
2.000e-05
3.000e-05
 200  210  220  230  240  250  260  270  280
FIG. 8: Behavior of convergence tests for a particle in circular
orbit at r = 6M . We show differences between simulations of
the ℓ = 2, m = 2 multipole moment using different step sizes
of 2 (ψ2), 4 (ψ4), 8 (ψ8), 16 (ψ16), 32 (ψ32) and 64 (ψ64) cells
per M . Displayed are the real part of the rescaled differences
δ4−2 = (ψ4 −ψ2) etc. of the field values at the position of the
particle, defined as in Fig. 5. The values have been rescaled
so that—for fourth order convergence—the curves should all
coincide. The upper panel corresponds to a set of simulations
where the particle traverses the cells away from their tips.
The curves do not coincide perfectly with each other, seem-
ingly indicating a failure of the convergence. The lower panel
was obtained in a simulation where the particle was carefully
positioned so as to pass through the tips of each cell it tra-
verses. This set of simulations passes the convergence test
more convincingly.
C. High-ℓ behavior of the multipole coefficients
Inspection of Eq. (4.11) reveals that a plot of Φ(µ)ℓ as
a function of ℓ (for a selected value of t) should display
a linear growth in ℓ for large ℓ. Removing the A(µ) term
should produce a constant curve, removing the B(µ) term
(given that C(µ) = 0) should produce a curve that decays
as ℓ−2, and finally, removing the D(µ) term should pro-
duce a curve that decays as ℓ−4. It is a powerful test of
the numerical methods to check whether these expecta-
tions are borne out by the numerical data. Fig. 10 plots
the remainders as obtained from our numerical simula-
tion, demonstrating the expected behavior. It displays,
on a logarithmic scale, the absolute value of ReΦR(+)ℓ, the
real part of the (+) component of the self-force. The orbit
is eccentric (p = 7.2, e = 0.5), and all components of the
self-force require regularization. The first curve (in trian-
gles) shows the unregularized multipole coefficients that
increase linearly in ℓ, as confirmed by fitting a straight
line to the data. The second curve (in squares) shows par-
tially regularized coefficients, obtained after the removal
of (ℓ + 1/2)A(µ); this clearly approaches a constant for
large values of ℓ. The curve made up of diamonds shows
the behavior after removal of B(µ); because C(µ) = 0, it
decays as ℓ−2, a behavior that is confirmed by a fit to
-8.000e-05
-4.000e-05
0.000e+00
4.000e-05
8.000e-05
 200  210  220  230  240  250  260  270  280
-8.000e-05
-4.000e-05
0.000e+00
4.000e-05
8.000e-05
 200  210  220  230  240  250  260  270  280
-8.000e-05
-4.000e-05
0.000e+00
4.000e-05
8.000e-05
 200  210  220  230  240  250  260  270  280
δ64-32
δ32-16
δ16-8
FIG. 9: Behavior of convergence tests for a particle in circular
orbit at r = 6M . We show differences between simulations
of the ℓ = 2, m = 2 multipole moment using different step
sizes of 8 (ψ8), 16 (ψ16), 32 (ψ32), and 64 (ψ64) cells per M .
Displayed are the real part of the rescaled differences δ16−8 =
ψ16−ψ8 etc. of the field values at the position of the particle,
defined as in Fig. 5. The values have been rescaled so that—
for second order convergence—the curves should all coincide.
The upper two panels correspond to simulations where the
second order algorithm was used throughout. For the topmost
one, care was taken to ensure that the particle passes through
the tip of each cell it traverses, while in the middle one no
such precaution was taken. Clearly the curves in the middle
panel do not coincide with each other, indicating a failure
of the second-order convergence of the code. The lower panel
was obtained in a simulation using the mixed-order algorithm
described in the text. While the curves still do not coincide
precisely, the observed behavior is much closer to the expected
one than for the purely second order algorithm.
the ℓ ≥ 5 part of the curve. Finally, after removal of
D(µ)/[(ℓ − 12 ) (ℓ +
)] the terms of the sum decrease in
magnitude as ℓ−4 for large values of ℓ, as derived in [15].
Each one of the last two curves would result in a con-
verging sum, but the convergence is much faster after
subtracting the D(µ) terms. We thereby gain more than
2 orders of magnitude in the accuracy of the estimated
Figure 10 provides a sensitive test of the implemen-
tation of both the numerical and analytical parts of the
calculation. Small mistakes in either one will cause the
difference in Eq. (4.11) to have a vastly different behav-
 0  2  4  6  8  10  12  14
ReΦ(+)
ReΦ(+)-A
ReΦ(+)-A-B
ReΦ(+)-A-B-D
FIG. 10: Multipole coefficients of the dimensionless self-force
ReΦR(+) for a particle on an eccentric orbit (p = 7.2, e =
0.5). The coefficients are extracted at t = 500M along the
trajectory shown in Fig. 12. The plots show several stages of
the regularization procedure, with a closer description of the
curves to be found in the text.
 0  2  4  6  8  10  12  14  16  18  20
FIG. 11: Multipole coefficients of ΦR(0) for a particle on a circu-
lar orbit. Note that ΦR(0)ℓ is linked to Φ
t via Φ
The multipole coefficients decay exponentially with ℓ until
ℓ ≈ 16, at which point numerical errors start to dominate.
D. Self-force on a circular orbit
For the case of a circular orbit, the regularization pa-
rameters A(0), B(0), and D(0) all vanish identically, so
that the (0) (or alternatively the t) component of the
self-force does not require regularization. Figure 11 thus
shows only one curve, with the magnitude of the multi-
pole coefficients decaying exponentially with increasing
As a final test, in Table I we compare our result for the
self-force on a particle in a circular orbit at r = 6M to
those obtained in [4, 14] using a frequency-domain code.
For a circular orbit, a calculation in the frequency domain
TABLE I: Results for the self-force on a scalar particle with
scalar charge q on a circular orbit at r0 = 6M . The
first column lists the results as calculated in this work us-
ing time-domain numerical methods, while the second and
third columns list the results as calculated in [4, 14] using
frequency-domain methods. For the t and φ components the
number of digits is limited by numerical roundoff error. For
the r component the number of digits is limited by the trun-
cation error of the sum of multipole coefficients.
This work: Previous work: Diaz-Rivera
time-domain frequency-domain [4] et. al. [14]
ΦRt 3.60339 × 10
−4 3.60907254 × 10−4
1.6767 × 10−4 1.67730 × 10−4 1.6772834 × 10−4
ΦRφ −5.30424 × 10
−3 −5.30423170 × 10−3
is more efficient, and we expect the results of [4, 14] to
be much more accurate than our own results. This fact
is reflected in the number of regularization coefficients
we can reliably extract from the numerical data, before
being limited by the accuracy of the numerical method:
the frequency-domain calculation found usable multipole
coefficients up to ℓ = 20, whereas our data for ΦR(0)ℓ is
dominated by noise by the time ℓ reaches 16. Figure 11
shows this behavior.
E. Accuracy of the numerical method
Several figures of merit can be used to estimate the
accuracy of numerical values for the self-force.
An estimate for the truncation error arising from cut-
ting short the summation in Eq. (4.11) at some ℓmax can
be calculated by considering the behavior of the remain-
ing terms for large ℓ. Detweiler et. al. [15] showed that
the remaining terms scale as ℓ−4 for large ℓ. They find
the functional form of the terms to be
EP3/2
(2ℓ− 3)(2ℓ− 1)(2ℓ+ 3)(2ℓ+ 5)
, (5.3)
where P3/2 = 36
2. We fit a function of this form to the
tail end of a plot of the multipole coefficients to find the
coefficient E in Eq. (5.3). Extrapolating to ℓ → ∞ we
find that the truncation error is
ℓ=ℓmax
[Eq. (5.3)] (5.4)
2Eℓmax
(2ℓmax + 3)(2ℓmax + 1)(2ℓmax − 1)(2ℓmax − 3)
(5.5)
where ℓmax is the value at which we cut the summation
short. For all but the special case of the (0) component
for a circular orbit, for which all regularization parame-
ters vanish identically, we use this approach to calculate
an estimate for the truncation error.
A second source of error lies in the numerical calcula-
tion of the retarded solution to the wave equation. This
error depends on the step size h used to evolve the field
forward in time. For a numerical scheme of a given con-
vergence order, we can estimate this discretization error
by extrapolating the differences of simulations using dif-
ferent step sizes down to h = 0. This is what was done
in the graphs shown in Sec. VB.
We display results for mildly eccentric orbits. A high
eccentricity causes ∂rΦ (displayed in Fig. 7) to be plagued
by high frequency noise produced by effects similar to
those described in Sec. VB. This makes it impossible to
reliably estimate the discretization error for these orbits.
We do not expect this to be very different from the errors
for mildly eccentric orbits.
Finally we compare our final results for the self-force
Fα to “reference values”. For circular orbits, frequency-
domain calculations are much more accurate than our
time-domain computations. We thus compare our results
to the results obtained in [4]. Table II lists typical values
for the various errors listed above.
error estimation mildly eccentric orbit
truncation error (M
Φ(+)) ≈ 2× 10
discretization error (M
∂rΦℓm) ≈ 10
comparison with reference values circular orbit
Ft 0.2%
Fr 0.04%
Fφ 2× 10
TABLE II: Estimated values for the various errors in the com-
ponents of the self-force as described in the text. We show
the truncation and discretization errors for a mildly eccentric
orbit and the total error for a circular orbit. The truncation
error is calculated using a plot similar to the one shown in
Fig. 16. The discretization error is estimated using a plot
similar to that in Fig. 7 for the ℓ = 2, m = 2 mode, and the
total error is estimated as the difference between our values
and those of [4]. We use p = 7.2 , e = 0.5 for the mildly
eccentric orbit. Note that we use the tetrad component Φ(+)
for the truncation error and the vector component ∂rΦ for
the discretization error. Both are related by the translation
table Eqs. (A6) – (A9), we expect corresponding errors to be
comparable for Φ(+) and ∂rΦ.
VI. SAMPLE RESULTS
In this section we describe some results of our numer-
ical calculation.
A. Mildly eccentric orbit
We choose a particle on an eccentric orbit with p = 7.2,
e = 0.5 which starts at r = pM/(1−e2), halfway between
 15  10  5  0  5  10
trajectory for p=7.2, e=0.5
FIG. 12: Trajectory of a particle with p = 7.2, e = 0.5. The
cross-hair indicates the point where the data for Fig. 10 was
extracted.
-0.014
-0.012
-0.01
-0.008
-0.006
-0.004
-0.002
 0.002
 100  200  300  400  500  600  700  800  900  1000
time/M
FIG. 13: Regularized dimensionless self-force M
and M
Fφ on a particle on an eccentric orbit with p = 7.2,
e = 0.5.
periastron and apastron. The field is evolved for 1000M
with a resolution of 16 grid points per M , both in the t
and r∗ directions, for ℓ = 0. Higher values of ℓ (and thus
m) require a corresponding increase in the number of
grid points used to achieve the same fractional accuracy.
Multipole coefficients for 0 ≤ ℓ ≤ 15 are calculated and
used to reconstruct the regularized self-force Fα along
the geodesic. Figure 13 shows the result of the calcula-
tion. For the choice of parameters used to calculate the
force shown in Fig. 13, the error bars corresponding to
the truncation error (which are already much larger than
than the discretization error) would be of the order of
the line thickness and have not been drawn.
Already for this small eccentricity, we see that the self-
force is most important when the particle is closest to the
black hole (ie. for 200M . t . 400M and 600M . t .
 60  50  40  30  20  10  0  10  20  30  40  50
trajectory for p=7.8001, e=0.9
FIG. 14: Trajectory of a particle on a zoom-whirl orbit with
p = 7.8001, e = 0.9. The cross-hairs indicate the positions
where the data shown in Fig. 16 and 17 was extracted.
800M); the self-force acting on the particle is very small
once the particle has moved away to r ≈ 15M .
B. Zoom-whirl orbit
Highly eccentric orbits are of most interest as sources
of gravitational radiation. For nearly parabolic orbits
with e . 1 and p & 6+2e, a particle revolves around the
black hole a number of times, moving on a nearly circu-
lar trajectory close to the event horizon (“whirl phase”),
before moving away from the black hole (“zoom phase”).
During the whirl phase the particle is in the strong field
region of the black hole, emitting copious amounts of
radiation. Figures 14 and 15 show the trajectory of a
particle and the force on such an orbit with p = 7.8001,
e = 0.9. Even more so than for the mildly eccentric
orbit discussed in Sec. VIA, the self-force (and thus the
amount of radiation produced) is much larger while the
particle is close to the black hole than when it zooms out.
Defining energy E per unit mass and angular momen-
tum L per unit mass in the usual way,
E = −
uα, L =
uα, (6.1)
and following eg. the treatment of Wald [16], Ap-
pendix C, it is easy to see that the rates of change Ė
and L̇ (per unit proper time) are directly related to com-
ponents of the acceleration aα (and therefore force) ex-
perienced by the particle via
Ė = −at, L̇ = aφ. (6.2)
The self-force shown in Fig. 15 therefore confirms our
näıve expectation that the self-force should decrease both
the energy and angular momentum of the particle as ra-
diation is emitted.
-0.025
-0.02
-0.015
-0.01
-0.005
 0.005
 500  1000  1500  2000  2500  3000
time/M
-0.002
-0.001
 0.001
 0.002
 0.003
 0.004
 2000 2050 2100 2150 2200
FIG. 15: Self-force acting on a particle. Shown is the dimen-
sionless self-force M
Fr and
Fφ on a zoom-whirl
orbit with p = 7.8001, e = 0.9. The inset shows a magni-
fied view of the self-force when the particle is about to enter
the whirl phase. No error bars showing an estimate error are
shown, as the errors shown eg. in Table II are to small to
show up on the graph. Notice that the self-force is essentially
zero during the zoom phase 500M . t . 2000M and reaches
a constant value very quickly after the particle enters into the
whirl phase.
It is instructive to have a closer look at the force acting
on the particle when it is within the zoom phase, and also
when it is moving around the black hole on the nearly cir-
cular orbit of the whirl phase. In Fig. 16 and Fig. 17 we
show plots of Φ(0)ℓ vs. ℓ after the removal of the A(µ),
B(µ), and D(µ) terms. While the particle is still zooming
in toward the black hole, Φ(0)ℓ behaves exactly as for the
mildly eccentric orbit described in Sec. VIA over the full
range of ℓ plotted; ie. the magnitude of each term scales
as ℓ0, ℓ−2 and ℓ−4, after removal of the A(µ), B(µ), and
D(µ) terms respectively. Close to the black hole, on the
other hand, the particle moves along a nearly circular tra-
jectory. If the orbit were perfectly circular for all times,
ie. ṙ ≡ 0, then the (0) component would not require reg-
ularization at all, and the multipole coefficients would
decay exponentially, resulting in a straight line on the
semi-logarithmic plot shown in Fig. 17. As the real orbit
is not precisely circular, curves eventually deviate from a
straight line. Removal of the A(µ) term is required almost
immediately (beginning with ℓ ≈ 3), while the D(µ) term
starts to become important only after ℓ ≈ 11. This shows
that there is a smooth transition from the self-force on a
circular orbit, which does not require regularization for
the t and φ components, to that of a generic orbit, for
which all components of the self-force require regulariza-
tion.
 0  2  4  6  8  10  12  14
Φ(0)-A
Φ(0)-A-B
Φ(0)-A-B-D
FIG. 16: Multipole coefficients of M
ReΦR(0) for a particle on
a zoom-whirl orbit (p = 7.8001, e = 0.9). The coefficients are
extracted at t = 2000M as the particle is about to enter the
whirl phase. As ṙ is non-zero, all components of the self-force
require regularization and we see that the dependence of the
multipole coefficients on ℓ is as predicted by Eq. 1.9. After the
removal of the regularization parameters A(µ), B(µ), and D(µ)
the remainder is proportional to ℓ0, ℓ−2 and ℓ−4 respectively.
 0  2  4  6  8  10  12  14
Φ(0)-A
Φ(0)-A-B
Φ(0)-A-B-D
FIG. 17: Multipole coefficients of ReΦR(0) for a particle on
a zoom-whirl orbit (p = 7.8001, e = 0.9). The coefficients
are extracted at t = 2150M while the particle is in the whirl
phase. The orbit is nearly circular at this time, causing the
dependence on ℓ after removal of the regularization parame-
ters to approximate that of a true circular orbit.
Acknowledgments
We thank Eric Poisson and Eran Rosenthal for useful
discussions and suggestions. This work was supported by
the Natural Sciences and Engineering Council of Canada.
APPENDIX A: TRANSLATION TABLES
We quote the results of [4] for the translation table be-
tween the modes Φℓm and the tetrad components Φ(µ)ℓm
with respect to the pseudo-Cartesian basis
eα(0) =
, 0, 0, 0
, (A1)
eα(1) =
f sin θ cosφ,
cos θ cosφ,− sinφ
r sin θ
, (A2)
eα(2) =
f sin θ sinφ,
cos θ sinφ,
r sin θ
, (A3)
eα(3) =
f cos θ,−1
sin θ, 0
, (A4)
and the complex combinations eα
:= eα
± ieα
eα(±) =
f sin θe±iφ,
cos θe±iφ,
±ie±iφ
r sin θ
. (A5)
With these, the spherical-harmonic modes Φ(µ)ℓm(t, r)
are given in terms of Φℓm(t, r) by
Φ(0)ℓm =
Φℓm, (A6)
Φ(+)ℓm =−
(ℓ+m− 1)(ℓ +m)
(2ℓ− 1)(2ℓ+ 1)
− ℓ− 1
Φℓ−1,m−1
(ℓ−m+ 1)(ℓ −m+ 2)
(2ℓ+ 1)(2ℓ+ 3)
Φℓ+1,m−1, (A7)
Φ(−)ℓm =
(ℓ −m− 1)(ℓ−m)
(2ℓ− 1)(2ℓ+ 1)
− ℓ− 1
Φℓ−1,m+1
(ℓ+m+ 1)(ℓ +m+ 2)
(2ℓ+ 1)(2ℓ+ 3)
Φℓ+1,m+1, (A8)
Φ(3)ℓm =
(ℓ −m)(ℓ+m)
(2ℓ− 1)(2ℓ+ 1)
Φℓ−1,m
(ℓ−m+ 1)(ℓ +m+ 1)
(2ℓ+ 1)(2ℓ+ 3)
Φℓ+1,m. (A9)
APPENDIX B: REGULARIZATION
PARAMETERS
For completeness we list the regularization parameters
as calculated in [4]. Quantities bearing a subscript “0”
are evaluated at the particle’s position.
A(0) =
0 + L
sign(∆), (B1)
A(+) = −eiφ0
0 + L
sign(∆), (B2)
A(3) = 0, (B3)
where f0 := 1 − 2M/r0 and sign(∆) is equal to +1 if
∆ > 0 and to −1 if ∆ < 0. We have, in addition, A(−) =
Ā(+), A(1) = Re[A(+)], and A(2) = Im[A(+)].
We also use
B(0) = −
Er0ṙ0√
0 + L
2)3/2
E + Er0ṙ0
0 + L
2)3/2
B(+) = e
Bc(+) − iB
, (B5)
Bc(+) =
0 + L
2)3/2
r20 + L
0 + L
2)3/2
f0 − 1
r20 + L
K, (B6)
Bs(+) = −
f0)ṙ0
r20 + L
E + (2−
f0)ṙ0
r20 + L
B(3) = 0. (B8)
In addition, B(−) = B̄(+), B(1) = Re[B(+)] =
cosφ0 + B
sinφ0, and B(2) = Im[B(+)] =
Bc(+) sinφ0 −B
(+) cosφ0.
Here, the rescaled elliptic integrals E and K are defined
E := 2
∫ π/2
(1− k sin2 ψ)1/2 dψ = F
; 1; k
K := 2
∫ π/2
(1− k sin2 ψ)−1/2 dψ = F
; 1; k
(B10)
in which k := L2/(r20 + L
We also use
C(µ) = 0 (B11)
D(0) = −
Er30(r
0 − L2)ṙ30
0 + L
2)7/2
E(r70 + 30Mr
0 − 7L2r50 + 114ML2r40 + 104ML4r20 + 36ML6)ṙ0
16r40
0 + L
2)5/2
Er30(5r
0 − 3L2)ṙ30
0 + L
2)7/2
E(r50 + 16Mr
0 − 3L2r30 + 42ML2r20 + 18ML4)ṙ0
16r20
0 + L
2)5/2
K, (B12)
D(+) = e
Dc(+) − iD
, (B13)
Dc(+) =
r30(r
0 − L2)ṙ40
0 + L
2)7/2
− r0ṙ
4(r20 + L
2)3/2
(3r70 + 6Mr
0 − L2r50 + 31ML2r40 + 26ML4r20 + 9ML6)ṙ20
0 + L
2)5/2
(3r70 + 8Mr
0 + L
2r50 + 26ML
2r40 + 22ML
4r20 + 8ML
16r60(r
0 + L
2)3/2
0 + 2Mr
0 + 4ML
r20 + L
0 − 3L2)ṙ40
0 + L
2)7/2
8(r20 + L
2)3/2
− (7r
0 + 12Mr
0 − L2r30 + 46ML2r20 + 18ML4)ṙ20
16r20
0 + L
2)5/2
− (7r
0 + 6Mr
0 + 6L
2r30 + 12ML
2r20 + 4ML
16r40(r
0 + L
2)3/2
r20 + L
K, (B14)
Ds(+) =
r20(r
0 − 7L2)(
f0 − 2)ṙ30
0 + L
2)5/2
− (2r
0 +Mr
0 + 5L
2r50 + 10ML
2r40 + 29ML
4r20 + 14ML
6)ṙ0
8r50L(r
0 + L
2)3/2
(r50 −Mr40 + 4L2r30 − 5ML2r20 + 2ML4)ṙ0
4r30L
0 + L
2)3/2
r20(r
0 − 3L2)(
f0 − 2)ṙ30
0 + L
2)5/2
(4r50 + 2Mr
0 + 7L
2r30 + 10ML
2r20 + 14ML
4)ṙ0
16r30L(r
0 + L
2)3/2
− (2r
0 − 2Mr20 + 5L2r0 − 8ML2)ṙ0
0 + L
2)3/2
K, (B15)
D(3) = 0. (B16)
And finally, D(−) = D̄(+), D(1) = Re[D(+)] =
cosφ0 + D
sinφ0, and D(2) = Im[D(+)] =
Dc(+) sinφ0 −D
(+) cosφ0.
APPENDIX C: PIECEWISE POLYNOMIALS
In two places in the numerical simulation we introduce
piecewise polynomials to approximate the scalar field ψℓm
across the world line, where it is continuous but not dif-
ferentiable. By a piecewise polynomial we mean a poly-
nomial of the form
p(t, r∗) =


n,m=0
unvm if r∗(u, v) > r∗0
n,m=0
unvm if r∗(u, v) < r∗0
, (C1)
where u = t−r∗, v = t+r∗ are characteristic coordinates,
r∗0 is the position of the particle at the time t(u, v), and
N is the order of the polynomial, which for our purposes
is N = 4 or less. The two sets of coefficients cnm and
c′nm are not independent of each other, but are linked
via jump conditions that can be derived from the wave
equation [Eq. (1.12)]. To do so, we rewrite the wave
equation in the characteristic coordinates u and v and
reintroduce the integral over the world line on the right-
hand side,
−4∂u∂vψ − V ψ =
Ŝ(τ)δ(u − up)δ(v − vp) dτ , (C2)
where Ŝ(τ) = −8πq Ȳℓm(π/2,φp(τ))
rp(τ)
is the source term and
quantities bearing a subscript p are evaluated on the
world line at proper time τ .
Here and in the following we use the notation
[∂nu∂
v ψ] = lim
[∂nu∂
v ψ(t0, r
0 + ǫ)− ∂nu∂mv ψ(t0, r∗0 − ǫ)]
to denote the jump in ∂nu∂
v ψ across the world line. First,
we notice that the source term does not contain any
derivatives of the Dirac δ-function, causing the solution
ψ to be continuous. This means that the zeroth-order
jump vanishes: [ψ] = 0. Our task is then to find the re-
maining jump conditions at a point (t0, r
0) for n,m ≤ 4.
Alternatively, instead of crossing the world line along a
line t = t0 = const we can also choose to cross along
lines of u = u0 = const or v = v0 = const, noting that
for a line of constant v the coordinate u runs from u0+ ǫ
to u0 − ǫ to cross from the left to the right of the world
line. Figure 18 provides a clearer description of the paths
taken.
(u0 − ǫ, v0)
(u0, v0 + ǫ)(u0 + ǫ, v0)
(u0, v0 − ǫ)
(t0, r
0) = (u0, v0)
FIG. 18: Paths taken in the calculation of the jump condi-
tions. (u0, v0) denotes an arbitrary but fixed point along the
world line γ. The wave equation is integrated along the lines
of constant u or v indicated in the sketch. Note that in order
to move from the domain on the left to the domain on the
right, u has to run from u0 + ǫ to u0 − ǫ. Where appropriate
we label quantities connected to the domain on the left by a
subscript “−” and quantities connected to the domain on the
right by “+”.
In order to find the jump [∂uψ] we integrate the wave
equation along the line u = u0 from v0 − ǫ to v0 + ǫ
∫ v0+ǫ
∂u∂vψdv −
∫ v0+ǫ
V ψdv =
Ŝ(τ)δ(u0 − up)
∫ v0+ǫ
δ(v − vp)dv dτ ,
which, after involving
∫ v0+ǫ
δ(v − vp)dv = θ(vp − v0 +
ǫ)θ(v0 − vp + ǫ) and δ(g(x)) = δ(x− x0)/ |g′(x0)|, yields
[∂uψ] = −
E − ṙ0
Ŝ(τ0), (C5)
where the overdot denotes differentiation with respect to
proper time τ .
Similarly, after first taking a derivative of the wave
equation with respect to v and integrating from u0+ ǫ to
u0 − ǫ, we obtain
∫ u0−ǫ
vψdu−
∫ u0−ǫ
V ψdu =
Ŝ(τ)
∫ u0−ǫ
δ(u− up)du δ′(v0 − vp)dτ .
We find
E + ṙ0
E + ṙp
Ŝ(τ)
|τ=τ0
. (C7)
Systematically repeating this procedure we find expres-
sions for the jumps in all the derivatives that are purely
in the u or v direction. Table III lists these results. Jump
[ψ] =0
[∂uψ] =−
bS(τ0), [∂vψ] =
bS(τ0)
bS(τ )
|τ=τ0
bS(τ )
|τ=τ0
V ξ0ξ̄
0 [∂uψ]−
bS(τ )
|τ=τ0
V ξ̄0ξ
0 [∂vψ] +
bS(τ )
|τ=τ0
|τ=τ0
+ 3ξ0ξ̄
0 ∂uV + ξ
0 ∂vV
[∂uψ] +
bS(τ )
|τ=τ0
|τ=τ0
+ 3ξ̄0ξ
0 ∂vV + ξ̄
0 ∂uV
[∂vψ]−
bS(τ )
|τ=τ0
TABLE III: Jump conditions for the derivatives purely in the
u or v directions. ṙ and r̈ are the particle’s radial velocity and
acceleration, respectively. They are obtained from the equa-
tion of motion for the particle. ξ̄ := E−ṙ
and ξ := E+ṙ
introduced for notational convenience. Quantities bearing a
subscript p are evaluated on the particle’s world line, while
quantities bearing a subscript 0 are evaluated at the parti-
cle’s current position. Derivatives of V with respect to either
u or v are evaluated as ∂uV = −
f∂rV and ∂vV =
f∂rV ,
respectively.
conditions for derivatives involving both u and v are ob-
tained directly from the wave equation [Eq. (C2)]. We
see that
[∂u∂vψ] = 0, (C8)
and taking an additional derivative with respect to u on
both sides reveals that
∂2u∂vψ
V [∂uψ] . (C9)
Systematically repeating this procedure we can find jump
conditions for each of the mixed derivatives by evaluating
∂n+1u ∂
[∂nu∂
v (V ψ)] , (C10)
where n,m ≥ 0 and derivatives of V with respect to
either u or v are evaluated as ∂uV = − 12f∂rV and ∂vV =
f∂rV , respectively.
The results of Table III and Eq. (C10) allow us to
express the coefficients of the left-hand polynomial in
Eq. (C1) in terms of the jump conditions and the co-
efficients of the right-hand side:
c′nm = cnm − [∂nu∂mv ψ] . (C11)
For N = 4 this leaves us with 25 unknown coefficients
cnm which can be uniquely determined by demanding
that the polynomial match the value of the field on the
25 grid points surrounding the particle. When we are
interested in integrating the polynomial, as in the case of
the potential term in the fourth-order algorithm, we do
not need all these terms. Instead, in order to calculate
e.g. the integral
V ψ du dv up to terms of order h5,
as is needed to achieve overall O(h4) convergence, it is
sufficient to include only terms such that n+m ≤ 2, thus
reducing the number of unknown coefficients to 6. In this
case Eq. (C1) becomes
p(t, r∗) =
m+n≤2
unvm if r∗(u, v) > r∗0
m+n≤2
unvm if r∗(u, v) < r∗0
. (C12)
The six coefficients can then be determined by matching
the polynomial to the field values at the six grid points
which lie within the past light cone of the grid point
whose field value we want to calculate.
[1] The LISA web site is located at
http://lisa.jpl.nasa.gov/.
[2] L. Barack and A. Ori, Phys. Rev. D 61, 061502 (2000),
gr-qc/9912010.
[3] C. Lousto, Class. Quantum Grav. 22, S543 (2005).
[4] R. Haas and E. Poisson, Phys. Rev. D 74,
044009 (pages 29) (2006), gr-qc/0605077, URL
http://link.aps.org/abstract/PRD/v74/e044009 .
[5] B. S. DeWitt and R. W. Brehme, Annals of Physics 9,
220 (1960).
[6] Y. Mino, M. Sasaki, and T. Tanaka, Phys. Rev. D 55,
3457 (1997), gr-qc/9606018.
[7] T. C. Quinn and R. M. Wald, Phys. Rev. D 56, 3381
(1997), gr-qc/9610053.
[8] T. C. Quinn, Phys. Rev. D 62, 064029 (2000), gr-
qc/0005030.
[9] S. Detweiler and B. F. Whiting, Phys. Rev. D 67, 024025
(2003), gr-qc/0202086.
[10] C. O. Lousto, Class. Quant. Grav. 22, S543 (2005), gr-
qc/0503001.
[11] R. H. Price, Phys. Rev. D 5, 2419 (1972).
[12] C. G. Darwin, Proc. R. Soc. A 249, 180 (1959).
http://lisa.jpl.nasa.gov/
http://link.aps.org/abstract/PRD/v74/e044009
[13] C. O. Lousto and R. H. Price, Phys. Rev. D 56, 6439
(1997), gr-qc/9705071.
[14] L. M. Diaz-Rivera, E. Messaritaki, B. F. Whit-
ing, and S. Detweiler, Physical Review D (Par-
ticles, Fields, Gravitation, and Cosmology) 70,
124018 (pages 14) (2004), gr-qc/0410011, URL
http://link.aps.org/abstract/PRD/v70/e124018.
[15] S. Detweiler, E. Messaritaki, and B. F. Whiting, Phys.
Rev. D 67, 104016 (2003), gr-qc/0205079.
[16] R. M. Wald, General relativity (University of Chicago
Press, Chicago, 1984), ISBN 0226870324.
http://link.aps.org/abstract/PRD/v70/e124018
