# Time domain calculation of the electromagnetic self-force on eccentric
geodesics in Schwarzschild spacetime

Roland Haas  Theoretical Astrophysics, California Institute of Technology,
Pasadena, CA 91125  Center for Relativistic Astrophysics, Georgia Institute of
Technology, Atlanta, GA 30332

###### Abstract

I calculate the self-force acting on a particle with electric charge $q$
moving on a generic geodesic around a Schwarzschild black hole. Using methods
similar to those developed for the scalar field case discussed in Haas (2007),
I investigate the relative sizes of the conservative (half-advanced plus half-
retarded) and dissipative (half-advanced minus half-retarded) pieces of the
self-force. I also display the regularization parameters used in the mode-sum
regularization scheme.

###### pacs:

04.25.-g, 04.40.-b, 41.60.-m, 45.50.-j, 02.60.Cb

##  I Introduction

This is the second paper of a series of papers studying the self-force on a
point particle in generic geodesic orbit around a Schwarzschild black hole. I
extend the previous calculation of the scalar self-force Haas and Poisson
(2006) to electromagnetism, studying in particular the effects of the
conservative part of the self-force.

A test particle in orbit around a black hole will follow a geodesic. Going
beyond the test mass limit, this is no longer true and the particle’s path
will deviate from a geodesic of the background spacetime. As seen from the
background spacetime, the particle is said to experiences a self-force due to
its interaction with its own field. In order to accurately model the motion of
the body, including its inspiral toward the black hole, I seek to evaluate the
self-force and calculate its effect on the motion. Several methods to achieve
this have been proposed in the literature Barack and Ori (2000); Vega and
Detweiler (2008); Barack et al. (2007). I elect to use the mode-sum
regularization scheme introduced by Barack and Ori Barack and Ori (2000),
which been proven to be highly accurate.

In this paper, rather than dealing with the gravitational problem, I focus on
the technically simpler problem of a point particle endowed with an electric
charge $q$ orbiting a Schwarzschild black hole of mass $M$. In this context I
use a numerical simulation to check the analytically calculated regularization
parameters used in the mode-sum regularization scheme, which I calculate in a
manner analogous to Haas and Poisson (2006). This calculation also makes it
possible to investigate the behaviour of the conservative (half-advanced plus
half-retarded) part of the self-force in the strong-field limit, extending
previous work by Pound and Poisson Pound and Poisson (2008). Different from
the scalar case calculation, where the conservative self-force is suppressed,
the conservative electromagnetic self-force appears at the same post Newtonian
order as the gravitational conservative self-force. Agreement, even if only
qualitative, between the results for the electromagnetic problem, where our
physical intuition allows us to understand the mechanisms at work, and those
for a point mass recently explored by Barack and Sago (2011); Warburton et al.
(2011) can thus help provide a clearer understanding of the mechanisms at work
in the gravitational case as well.

Throughout the paper I use geometrized units in which $G=c=1$ and the sign
conventions of Misner et al. (1973).

###  I.1 The problem

Since my approach is essentially identical to that described in Haas and
Poisson (2006) and Haas (2007) (paper I and paper II from now on), I will only
briefly introduce the required notation.

The first order self-force is calculated on a geodesic of Schwarzschild
spacetime, whose metric is written in Schwarzschild coordinates as

$\mathrm{d}s^{2}=-f\mathrm{d}t^{2}+f^{-1}\mathrm{d}r^{2}+r^{2}\mathrm{d}\Omega\text{,}$
(1)

where $f=\left(1-\frac{2M}{r}\right)$,
$\mathrm{d}\Omega=\left(\mathrm{d}\theta^{2}+\sin^{2}\theta\mathrm{d}\phi^{2}\right)$
is the metric on a two-sphere, and $t$, $r$, $\theta$ and $\phi$ are the usual
Schwarzschild coordinates. I numerically solve the Maxwell equations

$\displaystyle g^{\beta\gamma}\nabla_{\gamma}F_{\alpha\beta}(x)=4\pi
j_{\alpha}(x)\text{,}$ (2)
$\displaystyle\nabla_{[\gamma}F_{\alpha\beta]}(x)=0\text{,}$ (3)
$\displaystyle
j_{\alpha}(x)=q\int_{\gamma}u_{\alpha}(\tau)\delta_{4}\bm{(}x,z(\tau)\bm{)}\mathrm{d}\tau\text{,}$
(4)

where $\nabla_{\alpha}$ is the covariant derivative compatible with the metric
$g_{\alpha\beta}$, $F_{\alpha\beta}$ is the Faraday field tensor sourced by a
charge $q$ which moves along a world line $\gamma:\tau\mapsto z(\tau)$
parametrized by proper time $\tau$. The current density $j_{\alpha}(x)$
appearing on the right-hand side is written in terms of a scalarized four-
dimensional Dirac $\delta$-function
$\delta_{4}(x,x^{\prime})\equiv\delta(x_{0}-x^{\prime}_{0})\delta(x_{1}-x^{\prime}_{1})\delta(x_{2}-x^{\prime}_{2})\delta(x_{3}-x^{\prime}_{3})/\sqrt{-\det(g_{\alpha\beta})}$.
After having obtained the Faraday tensor I regularize it using the mode-sum
regularization scheme introduced by Barack and Ori Barack and Ori (2000)

$F^{R}_{(\mu)(\nu)}=F^{\text{ret}}_{(\mu)(\nu)}-q\sum_{\ell}\left[A_{(\mu)(\nu)}\Bigl{(}\ell+\frac{1}{2}\Bigr{)}+B_{(\mu)(\nu)}\right.\\\
\mbox{}\left.+\frac{C_{(\mu)(\nu)}}{\ell+\frac{1}{2}}+\frac{D_{(\mu)(\nu)}}{(\ell-\frac{1}{2})(\ell+\frac{3}{2})}+\cdots\right]\text{,}$
(5)

where indices in parenthesis $(\mu)$ signify components with respect to an
orthonormal tetrad $e^{\alpha}_{\ (\mu)}$ and the coefficients
$A_{(\mu)(\nu)}$, $B_{(\mu)(\nu)}$, $C_{(\mu)(\nu)}$, and $D_{(\mu)(\nu)}$ are
independent of $\ell$; they are listed in Appendix B. Finally I compute the
regularized self-force

$F^{\text{self}}_{\alpha}\equiv qF^{R}_{\alpha\beta}u^{\beta}$ (6)

from the regularized Faraday tensor and the four velocity of the particle.

###  I.2 Organization of this paper

In Sec. II I introduce the ideas behind the discretization scheme used in the
numerical simulation. Sec. III describes the choices I make in order to handle
the problems of specifying initial data and proper boundary conditions. In
Sec. VII I describe the tests I performed in order to validate my
implementation of the numerical method. Sec. VIII contains sample results for
a small number of representative simulations. Finally in Sec. IX I calculate
the conservative self-force for the same set of simulations. The appendices
contain technical details and an alternative calculation using the vector
potential instead of the Faraday tensor.

##  II Numerical method

In this section I describe the algorithm used to integrate the Maxwell
equations numerically. I use the second-order algorithm introduced by Lousto
and Price Lousto and Price (1997) suitably extended to handle a coupled system
of equations.

###  II.1 Wave equations for the Faraday tensor

I introduce a vector potential $A_{\alpha}$ in terms of which the Faraday
tensor is given by

$\displaystyle F_{\alpha\beta}$
$\displaystyle=A_{\beta,\alpha}-A_{\alpha,\beta}\text{,}$ (7)

where a comma denotes an ordinary derivative. I use vector spherical harmonics
$Z^{\ell m}_{A}=\partial_{A}Y^{\ell m}$ and $X^{\ell
m}_{A}={\epsilon_{A}}^{B}\partial_{B}Y^{\ell m}$, where $\epsilon_{AB}$ is the
Levi-Civita tensor associated with the metric $\Omega_{AB}$ on the two-sphere
($\epsilon_{\theta\phi}=\sin\theta$), as introduced in Regge (1957); Martel
and Poisson (2005). I decompose the vector potential and the current density
into

$\displaystyle A_{a}(t,r,\theta,\phi)$ $\displaystyle=A^{\ell
m}_{a}(t,r)Y_{\ell m}(\theta,\phi)\text{,}$ (8a) $\displaystyle
j_{a}(t,r,\theta,\phi)$ $\displaystyle=j^{\ell m}_{a}(t,r)Y_{\ell
m}(\theta,\phi)$ for $a=t,r$, (8b) $\displaystyle A_{A}(t,r,\theta,\phi)$
$\displaystyle=v_{\ell m}(t,r)Z_{A}^{\ell m}(\theta,\phi)$
$\displaystyle\quad+\tilde{v}_{\ell m}(t,r)X_{A}^{\ell
m}(\theta,\phi)\text{,}$ (8c) $\displaystyle j_{A}(t,r,\theta,\phi)$
$\displaystyle=j^{\text{even}}_{\ell m}(t,r)Z_{A}^{\ell m}(\theta,\phi)$
$\displaystyle\quad+j^{\text{odd}}_{\ell m}(t,r)X_{A}^{\ell m}(\theta,\phi)$
for $A=\theta,\phi$, (8d)

where a summation over $\ell$ and $m$ is implied. Substituting these into Eq.
(2) I arrive at two sets of coupled equations for the even ($A^{\ell m}_{a}$,
$v_{\ell m}$) and odd ($\tilde{v}_{\ell m}$) modes

$\displaystyle-f\frac{\partial^{2}A_{t}^{\ell
m}}{\partial{r}^{2}}+f\frac{\partial^{2}A_{r}^{\ell m}}{\partial t\partial
r}-\frac{2f}{r}\frac{\partial A_{t}^{\ell m}}{\partial
r}+\frac{2f}{r}\frac{\partial A_{r}^{\ell m}}{\partial t}$
$\displaystyle\mbox{}-\frac{\ell(\ell+1)}{r^{2}}\frac{\partial v_{\ell
m}}{\partial t}+\frac{\ell(\ell+1)}{r^{2}}A_{t}^{\ell m}=4\pi j_{t}^{\ell
m}\text{,}$ (9a) $\displaystyle f^{-1}\frac{\partial^{2}A_{r}^{\ell
m}}{\partial{t}^{2}}-f^{-1}\frac{\partial^{2}A_{t}^{\ell m}}{\partial
t\partial r}-\frac{\ell(\ell+1)}{r^{2}}\frac{\partial v_{\ell m}}{\partial r}$
$\displaystyle\mbox{}+\frac{\ell(\ell+1)}{r^{2}}A_{r}^{\ell m}=4\pi
j_{r}^{\ell m}\text{,}$ (9b) $\displaystyle f^{-1}\frac{\partial^{2}v_{\ell
m}}{\partial{t}^{2}}-f\frac{\partial^{2}v_{\ell
m}}{\partial{r}^{2}}-\frac{2M}{r^{2}}\frac{\partial v_{\ell m}}{\partial
r}+f\frac{\partial A_{r}^{\ell m}}{\partial r}$
$\displaystyle\mbox{}-f^{-1}\frac{\partial A_{t}^{\ell m}}{\partial
t}+\frac{2M}{r^{2}}A_{r}^{\ell m}=4\pi j^{\text{even}}_{\ell m}\text{,}$ (9c)
$\displaystyle f^{-1}\frac{\partial^{2}\tilde{v}_{\ell
m}}{\partial{t}^{2}}-f\frac{\partial^{2}\tilde{v}_{\ell
m}}{\partial{r}^{2}}-\frac{2M}{r^{2}}\frac{\partial\tilde{v}_{\ell
m}}{\partial r}$
$\displaystyle\mbox{}+\frac{\ell(\ell+1)}{r^{2}}\tilde{v}_{\ell m}=4\pi
j^{\text{odd}}_{\ell m}\text{,}$ (9d)

where

$\displaystyle j^{\ell m}_{t}=-\frac{qf}{r_{0}^{2}}\bar{Y}^{\ell
m}(\frac{\pi}{2},\varphi_{0})\delta(r-r_{0})\text{,}$ (10a) $\displaystyle
j^{\ell m}_{r}=\frac{q\dot{r}_{0}}{Er_{0}^{2}}\bar{Y}^{\ell
m}(\frac{\pi}{2},\varphi_{0})\delta(r-r_{0})\text{,}$ (10b) $\displaystyle
j^{\text{even}}_{\ell m}=-\frac{imqfJ}{\ell(\ell+1)Er_{0}^{2}}\bar{Y}^{\ell
m}(\frac{\pi}{2},\varphi_{0})\delta(r-r_{0})\text{,}$ (10c) $\displaystyle
j^{\text{odd}}_{\ell
m}=-\frac{qfJ}{\ell(\ell+1)Er_{0}^{2}}\partial_{\theta}\bar{Y}^{\ell
m}(\frac{\pi}{2},\varphi_{0})\delta(r-r_{0})\text{.}$ (10d)

In the equation above an overbar denotes complex conjugation, an overdot
denotes differentiation with respect to $\tau$, $E=-u_{t}$ is the particle’s
conserved energy per unit mass, $J=u_{\phi}$ its conserved angular momentum
per unit mass, and $u^{\alpha}=\frac{\mathrm{d}z^{\alpha}}{\mathrm{d}\tau}$ is
its four velocity. Quantities bearing a subscript “$0$” are evaluated at the
particle’s position; they are functions of $\tau$ that are obtained by solving
the geodesic equation

$u^{\beta}\nabla_{\beta}u^{\alpha}=0$ (11)

in the background spacetime. Without loss of generality, I have confined the
motion of the particle to the equatorial plane $\theta=\frac{\pi}{2}$.

The three even mode equations Eq. (9a) – Eq. (9c) are not yet amenable to a
numerical treatment, as they are highly coupled. In order to obtain a more
convenient set of equation I define the auxiliary fields

$\displaystyle\psi^{\ell m}\equiv-r^{2}\left(\frac{\partial A_{t}^{\ell
m}}{\partial r}-\frac{\partial A_{r}^{\ell m}}{\partial t}\right)\text{,}$
(12) $\displaystyle\chi^{\ell m}\equiv f\,\left(A_{r}^{\ell m}-\frac{\partial
v^{\ell m}}{\partial r}\right)\text{,}$ (13) $\displaystyle\xi^{\ell m}\equiv
A_{t}^{\ell m}-\frac{\partial v^{\ell m}}{\partial t}\text{,}$ (14)

which, up to scaling factors, are just the even multipole moments of the $tr$,
$r\phi$ and $t\phi$ components of the Faraday tensor

$\displaystyle F_{tr}$ $\displaystyle=\sum_{\ell,m}\frac{\psi^{\ell
m}}{r^{2}}\,Y^{\ell m}\text{,}$ (15) $\displaystyle F_{tA}$
$\displaystyle=\sum_{\ell,m}(-\xi^{\ell m}\,Z^{\ell m}_{A}+\tilde{v}^{\ell
m}_{,t}\,X^{\ell m}_{A})\text{,}$ (16) $\displaystyle F_{rA}$
$\displaystyle=\sum_{\ell,m}(\frac{\chi^{\ell m}}{f}\,Z^{\ell
m}_{A}+\tilde{v}^{\ell m}_{,r}\,X^{\ell m}_{A})\text{,}$ (17) $\displaystyle
F_{\theta\varphi}$ $\displaystyle=\sum_{\ell,m}\tilde{v}_{\ell m}\,(X^{\ell
m}_{\phi,\theta}-X^{\ell m}_{\theta,\phi})$
$\displaystyle=-\sum_{\ell,m}\ell(\ell+1)\tilde{v}_{\ell
m}\,\sin(\theta)Y^{\ell m}\text{.}$ (18)

I note that the three fields $\psi^{\ell m}$, $\chi^{\ell m}$ and $\xi^{\ell
m}$ are not independent of each other, in fact knowledge of $\psi^{\ell m}$ is
sufficient to reconstruct $\chi^{\ell m}$ and $\xi^{\ell m}$. Eq. (9a) can be
rearranged to yield

$\displaystyle\xi^{\ell m}=-\frac{f}{\ell(\ell+1)}\frac{\partial\psi^{\ell
m}}{\partial r}-\frac{4\pi}{\ell(\ell+1)}j_{t}^{\ell m}\text{,}$ (19) and
similarly from Eq. (9b) $\displaystyle\chi^{\ell
m}=-\frac{1}{\ell(\ell+1)}\frac{\partial\psi^{\ell m}}{\partial t}-\frac{4\pi
f}{\ell(\ell+1)}j_{r}^{\ell m}\text{,}$ (20)

showing that knowledge of $\psi^{\ell m}$ is sufficient to reconstruct the
even multipole components of the Faraday tensor. In this work however I choose
to solve for $\chi^{\ell m}$ and $\xi^{\ell m}$ directly, rather than to
numerically differentiate $\psi^{\ell m}$ to obtain them. The gain in speed
from reducing the number of equations does not seem to offset the additional
time required to calculate $\psi^{\ell m}$ accurately enough to obtain good
approximations for its derivatives at the location of the particle. In this
approach Eqs. (19) and (20) are treated as constraints that the dynamical
variables have to satisfy.

Dropping the superscripts $\ell$, $m$ for notational convenience and following
Cunningham and Price (1979) I form linear combinations of derivatives of
Eqs.(9a) – (9c). I use
$[\partial_{r}(r^{2}\,\text{\eqref{eqn:mw-r}})-\partial_{t}(r^{2}\,\text{\eqref{eqn:mw-t}})]$
for $\psi$ and find

$\displaystyle\frac{\partial^{2}\psi}{\partial{r^{*}}^{2}}-\frac{\partial^{2}\psi}{\partial{t}^{2}}-V\psi=S_{\psi}\text{,}$
(21a) $\displaystyle S_{\psi}=4\pi f\left[\frac{\partial(r^{2}j_{t}^{\ell
m})}{\partial r}-\frac{\partial(r^{2}j_{r}^{\ell m})}{\partial
t}\right]\text{,}$ (21b) where $V=\ell(\ell+1)\frac{r-2M}{r^{3}}$ and
$r^{*}=r+2M\ln(\frac{r}{2M}-1)$ is the Regge-Wheeler tortoise coordinate.
Similarly I use
$[f\,\text{\eqref{eqn:mw-r}}-\partial_{r}(f\,\text{\eqref{eqn:mw-v}})]$ for
$\chi$ and $[\text{\eqref{eqn:mw-t}}-\partial_{t}\text{\eqref{eqn:mw-v}}]$ for
$\xi$. I find
$\displaystyle\frac{\partial^{2}\chi}{\partial{r^{*}}^{2}}-\frac{\partial^{2}\chi}{\partial{t}^{2}}-V\chi$
$\displaystyle=S_{\chi}\text{,}$ (21c) $\displaystyle S_{\chi}=4\pi
f\Biggl{[}\frac{\partial(fj^{\text{even}}_{\ell m})}{\partial r}$
$\displaystyle-fj_{r}^{\ell m}\Biggr{]}\text{,}$ (21d)
$\displaystyle\frac{\partial^{2}\xi}{\partial{r^{*}}^{2}}-\frac{\partial^{2}\xi}{\partial{t}^{2}}-V\xi-
V_{\xi}\psi$ $\displaystyle=S_{\xi}\text{,}$ (21e) $\displaystyle S_{\xi}=4\pi
f\Biggl{[}\frac{\partial(fj^{\text{even}}_{\ell m})}{\partial t}$
$\displaystyle-fj_{t}^{\ell m}\Biggr{]}\text{,}$ (21f)

where $V_{\xi}=\frac{2(r-3M)(r-2M)}{r^{5}}$. While still partially coupled
Eqs. (21b) – (21f) are much easier to deal with than the original set Eqs.
(9a) – (9c). The coupling is in the form of a staggering, which allows me to
first solve for $\psi$ and use this result in the calculation of $\xi$. On the
other hand, the source terms appearing on the right-hand side contain
derivatives of Dirac’s $\delta$-function resulting in fields that are
discontinuous at the location of the particle. Lousto’s scheme is designed to
cope with precisely this situation.

I derive explicit expressions for the source terms $S_{\alpha}$ on the right
hand sides

$\displaystyle
S_{\alpha}=G_{\alpha}(t)f_{0}\delta(r-r_{0})+F_{\alpha}(t)f\delta^{\prime}(r-r_{0})\text{,}$
(22a) $\displaystyle G_{\psi}(t)=-\frac{4\pi
q}{E^{2}}f_{0}\,\left(\ddot{r}_{0}-\frac{im\dot{r}_{0}J}{r_{0}^{2}}\right)\bar{Y}_{\ell
m}(\frac{\pi}{2},\varphi_{0})\text{,}$ (22b) $\displaystyle F_{\psi}(t)=4\pi
qf_{0}\left(\frac{\dot{r}_{0}^{2}}{E^{2}}-1\right)\bar{Y}_{\ell
m}(\frac{\pi}{2},\varphi_{0})\text{,}$ (22c) $\displaystyle
G_{\chi}(t)=-\frac{4\pi q\dot{r}_{0}}{Er_{0}^{2}}f_{0}\bar{Y}_{\ell
m}(\frac{\pi}{2},\varphi_{0})\text{,}$ (22d) $\displaystyle
F_{\chi}(t)=-\frac{4\pi qJim}{E\ell(\ell+1)r_{0}^{2}}f_{0}^{2}\bar{Y}_{\ell
m}(\frac{\pi}{2},\varphi_{0})\text{,}$ (22e)
$\displaystyle\begin{split}G_{\xi}(t)&=-4\pi
q\biggl{\\{}\frac{Jim}{E^{2}\ell(\ell+1)r_{0}^{2}}\biggl{[}\biggl{(}\frac{2M}{r_{0}^{2}}-\frac{2f_{0}}{r_{0}}\biggr{)}\dot{r}_{0}\\\
&\qquad-\frac{imJ}{r_{0}^{2}}\biggr{]}-\frac{1}{r_{0}^{2}}\biggr{\\}}f_{0}\bar{Y}_{\ell
m}(\frac{\pi}{2},\varphi_{0})\text{,}\end{split}$ $\displaystyle
F_{\xi}(t)=\frac{4\pi
qJim\dot{r}_{0}}{E^{2}\ell(\ell+1)r_{0}^{2}}f_{0}^{2}\bar{Y}_{\ell
m}(\frac{\pi}{2},\varphi_{0})\text{.}$ (22f)

My functions $G_{\alpha}$ and $F_{\alpha}$ correspond to $G/f_{0}$ and $F/f$
in Lousto and Price (1997), respectively, they are independent of $r$ (but do
contain terms in $r_{0}(t)$). I prefer this form of the source terms over the
form given in Lousto and Price (1997) since it simplifies the integral over
the source term Eq. (3.6) of Lousto and Price (1997)

$\iint\mathrm{d}AS=2\int_{t_{1}}^{t_{2}}\biggl{[}\frac{G\bm{(}r_{0}(t),t\bm{)}}{1-2M/r_{0}(t)}\\\
-\frac{\partial}{\partial
r}\biggl{(}\frac{F(r,t)}{1-2M/r}\biggr{)}\bigg{|}_{r=r_{0}(t)}\biggr{]}\mathrm{d}t\\\
\pm
2\frac{F\bm{(}r_{0}(t_{1}),t_{1}\bm{)}}{[1-2M/r_{0}(t_{1})]^{2}}[1\mp\dot{r}^{*}_{0}(t_{1})]^{-1}\\\
\pm
2\frac{F\bm{(}r_{0}(t_{2}),t_{2}\bm{)}}{[1-2M/r_{0}(t_{2})]^{2}}[1\pm\dot{r}^{*}_{0}(t_{2})]^{-1}\text{.}$
(23)

Since $G^{\text{Lousto}}=f_{0}G_{\alpha}(t)$ and
$G^{\text{Lousto}}=fF_{\alpha}(t)$, the first term in square brackets inside
the integral simplifies, while the second term vanishes completely.
$F_{\alpha}$ only appears in the boundary terms.

###  II.2 Constraint equations

The full set of Maxwell equations consists of the inhomogeneous equations Eq.
(2) as well as the homogeneous constraints Eq. (3) which have to be satisfied
by a solution to Eq. (2). In the usual approach introducing a vector potential
$A_{\alpha}$ implies that the constraints are identically satisfied since they
reduce to the Bianchi identities for the second derivatives of $A_{\alpha}$.
When solving for the components of the Faraday tensor directly there is no a
priory guarantee that a solution to Eq. (21b) – (21f), and (9d) satisfies Eq.
(3). It turns out, however, that a decomposition into spherical harmonics is
sufficient to show that all but one of the constraints are identically
satisfied. The one that is not identically true is the $tr\varphi$ (or
$tr\theta$) equation, which in terms of $\psi$, $\chi$ and $\xi$ reads

$\frac{\psi}{r^{2}}-\frac{\chi_{,t}}{f}+\xi_{,r}=0\text{.}$ (24)

If the fields satisfy the sourced Maxwell equations Eqs. (9a), (9b), then Eq.
(24) is just the evolution equation for $\psi$. Thus Eq. (24) is valid
whenever $\psi$ satisfies the consistency relations Eq. (19) and (20).

Analytically then, the situation is clear. Given a set of compatible initial
conditions for $\psi$, $\chi$ and $\xi$ which initially satisfy the constraint
equations, a solution to the system of Eq. (21b) – (21f), (9d) satisfies the
full set of Maxwell equations at all later times, too.

Numerically I monitor but do not enforce Eq. (19) and (20). I generally find
that violations of the constraints are at least three orders of magnitude
smaller than the field quantities themselves. Figures 1 and 2 compare $\chi$
obtained from its evolution equations to that obtained from Eq. (20).

Figure 1: Violations of the constraint
$Z_{\chi}=\chi+\frac{1}{\ell(\ell+1)}\frac{\partial\psi}{\partial t}=0$ in the
vacuum region away from the location of particle. I plot the $\chi$ and
$\log_{10}\left\lvert{Z_{\chi}}\right\rvert$ as obtained on a spatial slice at
time $t=600\,M$. For this slightly eccentric orbit ($p=7.0$, $e=0.3$) using a
stepsize $h=1/512M$ the errors in the $\ell=2$, $m=2$ mode are at least three
orders of magnitude smaller than the field values. The exponentially growing
signal between $300M\lesssim r^{*}\lesssim 500$ is a remnant of the initial
data pulse travelling outward.  Figure 2: Violations of the constraint
$Z_{\chi}=\chi+\frac{1}{\ell(\ell+1)}\frac{\partial\psi}{\partial t}=0$ at the
location of the particle as a function of time. I display $\chi$ and
$\log_{10}\left\lvert{Z_{\chi}}\right\rvert$ for the $\ell=5$, $m=3$ mode of a
particle on an eccentric orbit with $p=7.8001$, $e=0.9$ with stepsize
$h=1/256\,M$. During the time $400\,M\lesssim t\lesssim 800\,M$ the particle
is in the whirl phase. The exponentially decaying signal before $t\approx
250M$ is the initial data pulse.

###  II.3 Monopole mode

For the electromagnetic field, the monopole mode $\ell=0$ is non-radiative.
The vector harmonics $Z^{\ell m}_{A}$ and $X^{\ell m}_{A}$ cannot be defined
in this case and the only surviving multipole mode is $\psi$. For the monopole
case Eq. (21b) reduces to a wave equation in flat space

$\frac{\partial^{2}\psi}{\partial{r^{*}}^{2}}-\frac{\partial^{2}\psi}{\partial{t}^{2}}=4\pi
f\left[\frac{\partial(r^{2}j_{t}^{0,0})}{\partial
r}-\frac{\partial(r^{2}j_{r}^{0,0})}{\partial t}\right]\text{,}$ (25)

which is simple enough so that I can solve it analytically. A straightforward
calculation shows that

$\psi(t,r^{*})=-\sqrt{4\pi}q\theta(r^{*}-r^{*}_{0}(t))$ (26)

satisfies Eq. (25) and corresponds to no outgoing radiation
$(\partial_{t}-\partial_{r^{*}})\psi=0$ at the event horizon and no ingoing
radiation $(\partial_{t}+\partial_{r^{*}})\psi=0$ at spatial infinity.

###  II.4 Discretization—even sector

Lousto’s method is directly applicable to terms of the form
$-\frac{\partial^{2}\psi}{\partial{t}^{2}}+\frac{\partial^{2}\psi}{\partial{r^{*}}^{2}}$,
$V(r)\psi$ (ie. the wave operator and potential terms) on the left-hand side
of the equation and the source terms $S_{\alpha}(t)$ on the right hand side.
Here $\psi$ is used as a placeholder for any one of $\psi$, $\chi$ or $\xi$;
$V(r)$ is an expression depending only on $r$. I discretize these as

$\displaystyle\iint_{\text{\hbox
to0.0pt{cell\hss}}}\,\mathrm{d}u\,\mathrm{d}v\,\left(-\frac{\partial^{2}\psi}{\partial{t}^{2}}+\frac{\partial^{2}\psi}{\partial{r^{*}}^{2}}\right)=-4\left[\psi_{3}+\psi_{2}-\psi_{1}-\psi_{4}\right]\text{,}$
(27) $\displaystyle\iint_{\text{\hbox
to0.0pt{cell\hss}}}\,\mathrm{d}u\,\mathrm{d}v\,V(r)\psi=\begin{cases}h^{2}V_{0}\,\sum_{i}\psi_{i}+O(h^{4})&\text{vacuum
cells}\\\ V_{0}\,\sum_{i}A_{i}\psi_{i}+O(h^{3})&\text{sourced
cells,}\end{cases}$ (28)

and

$\displaystyle\iint_{\text{\hbox
to0.0pt{cell\hss}}}\,\mathrm{d}u\,\mathrm{d}v\,S_{\alpha}(t)$
$\displaystyle=2\int_{t_{1}}^{t_{2}}G_{\alpha}\bm{(}t,r_{0}(t)\bm{)}\,\mathrm{d}t$
$\displaystyle\qquad\pm\frac{2F_{\alpha}\bm{(}t_{1},r_{0}(t_{1})\bm{)}}{1-2M/r(t_{1})}[1\mp\dot{r}_{0}(t_{1})/E]^{-1}$
$\displaystyle\qquad\pm\frac{2F_{\alpha}\bm{(}t_{2},r_{0}(t_{2})\bm{)}}{1-2M/r(t_{2})}[1\pm\dot{r}_{0}(t_{2})/E]^{-1}\text{,}$
(29)

where $u=t-r^{*}$, $v=t+r^{*}$ are null coordinates, $\psi_{1}$,…,$\psi_{4}$
refer to values of the field at the points labelled $1$,…,$4$ in Fig. 3,
$h=\Delta_{t}=\Delta_{r^{*}}/2$ is the step size, $V_{0}$ is the value of the
potential at the centre of the cell, $A_{1}$,…,$A_{4}$ are the areas indicated
in Fig. 3 and $t_{1}$ and $t_{2}$ are the times at which the particle enters
and leaves the cell, respectively.

Figure 3: Points used to calculate the integral over the potential terms. Grid
points are indicated by blue circles.

Spelled out explicitly the evolution equations for vacuum cells are

$\displaystyle\psi_{3}$
$\displaystyle=-\psi_{2}+(1-\frac{h^{2}}{2}V_{0})(\psi_{1}+\psi_{4})\text{,}$
(30a) $\displaystyle\chi_{3}$
$\displaystyle=-\chi_{2}+(1-\frac{h^{2}}{2}V_{0})(\chi_{1}+\chi_{4})\text{,}$
(30b) $\displaystyle\xi_{3}$
$\displaystyle=-\xi_{2}+(1-\frac{h^{2}}{2}V_{0})(\xi_{1}+\xi_{4})$
$\displaystyle\quad-\frac{h^{2}}{4}V_{\xi,0}(\psi_{1}+\psi_{2}+\psi_{3}+\psi_{4})\text{,}$
(30c)

and for sourced cells

$\displaystyle\psi_{3}$
$\displaystyle=-[1+\frac{V_{0}}{4}(A_{2}-A_{3})]\psi_{2}+[1-\frac{V_{0}}{4}(A_{4}+A_{3})]\psi_{4}$
$\displaystyle\quad+[1-\frac{V_{0}}{4}(A_{1}+A_{3})]\psi_{1}$
$\displaystyle\quad-\frac{1}{4}(1-\frac{V_{0}}{4}A_{3})\iint\mathrm{d}u\,\mathrm{d}v\,S_{\psi}(t)\text{,}$
(31a) $\displaystyle\chi_{3}$
$\displaystyle=-[1+\frac{V_{0}}{4}(A_{2}-A_{3})]\chi_{2}+[1-\frac{V_{0}}{4}(A_{4}+A_{3})]\chi_{4}$
$\displaystyle\quad+[1-\frac{V_{0}}{4}(A_{1}+A_{3})]\chi_{1}$
$\displaystyle\quad-\frac{1}{4}(1-\frac{V_{0}}{4}A_{3})\iint\mathrm{d}u\,\mathrm{d}v\,S_{\chi}(t)\text{,}$
(31b) $\displaystyle\xi_{3}$
$\displaystyle=-[1+\frac{V_{0}}{4}(A_{2}-A_{3})]\xi_{2}+[1-\frac{V_{0}}{4}(A_{4}+A_{3})]\xi_{4}$
$\displaystyle\quad+[1-\frac{V_{0}}{4}(A_{1}+A_{3})]\xi_{1}$
$\displaystyle\quad-\frac{1}{4}V_{\xi,0}(A_{1}\psi_{1}+A_{2}\psi_{2}+A_{3}\psi_{3}+A_{4}\psi_{4})$
$\displaystyle\quad-\frac{1}{4}(1-\frac{V_{0}}{4}A_{3})\iint\mathrm{d}u\,\mathrm{d}v\,S_{\xi}(t)\text{.}$
(31c)

###  II.5 Discretization–odd sector

When written in terms of $r^{*}$, Eq. (9d), which governs the odd modes
$\tilde{v}^{\ell m}$, is

$\displaystyle\frac{\partial^{2}\tilde{v}_{\ell
m}}{\partial{r^{*}}^{2}}-\frac{\partial^{2}\tilde{v}_{\ell
m}}{\partial{t}^{2}}-\frac{\ell(\ell+1)(r-2M)}{r^{3}}\tilde{v}_{\ell m}=-4\pi
fj^{\text{odd}}_{\ell m}\text{,}$ (32) $\displaystyle j^{\text{odd}}_{\ell
m}=-\frac{qJ}{\ell(\ell+1)Er_{0}^{2}}\partial_{\theta}\bar{Y}^{\ell
m}(\frac{\pi}{2},\varphi_{0})\delta(r^{*}-r^{*}_{0})\text{.}$ (33)

Eq. (32) is of the form of the scalar wave equation discussed in paper II. I
re-use the fourth order numerical code described there with
$V=\frac{\ell(\ell+1)(r-2M)}{r^{3}}$,
$S=4\pi\frac{qfJ}{\ell(\ell+1)Er_{0}^{2}}\partial_{\theta}\bar{Y}^{\ell
m}(\frac{\pi}{2},\varphi_{0})$. This yields accurate results for $\tilde{v}$
and its derivatives.

##  III Initial values and boundary conditions

I follow the approach detailed in paper II for the scalar self-force and do
not specify physical initial data or an outer boundary condition. I
arbitrarily choose the fields to vanish on the characteristic slices
$u=u_{0}=t_{0}-r^{*}_{0}$ and $v=v_{0}=t_{0}+r^{*}_{0}$

$\psi(u=u_{0})=\psi(v=v_{0})=0\text{,}$ (34)

thereby adding a certain amount of spurious waves to the solution which show
up as an initial burst.

I implement ingoing wave boundary conditions near the event horizon,
sufficiently close that _numerically_ $r\approx 2M$, so that the potential
terms in Eqs. (21b) – (21f) vanish. This happens at $r^{*}\approx-73\,M$ and I
implement the ingoing waves condition $\partial_{u}\psi=0$ there. Near the
outer boundary this is not possible, since the potential decays slowly.
Instead I choose to evolve the full domain of dependence of the initial data
surface, hiding the effects of the boundary.

##  IV Particle motion

I use the same approach as described in paper II to evolve the particle’s
motion, i.e. I introduce the semi-latus rectum $p$, the eccentricity $e$ and a
fictitious angle $\chi$, not to be confused with the Faraday tensor component
$\chi$ defined in Eq. (14), such that

$r(\tau)=\frac{pM}{1+e\cos\chi(\tau)}\text{.}$ (35)

The evolution is then governed by

$\displaystyle\begin{split}\frac{\mathrm{d}\chi}{\mathrm{d}t}=\frac{(p-2-2e\cos\chi)(1+e\cos\chi)^{2}}{(Mp^{2})}\\\
\mbox{}\times\sqrt{\frac{p-6-2e\cos\chi}{(p-2-2e)(p-2+2e)}}\text{,}\end{split}$
(36)
$\displaystyle\frac{\mathrm{d}\varphi}{\mathrm{d}t}=\frac{(p-2-2e\cos\chi)(1+e\cos\chi)^{2}}{p^{3/2}M\sqrt{(p-2-2e)(p-2+2e)}}\text{.}$
(37)

I use the embedded Runge-Kutta-Fehlberg (4, 5) algorithm 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.

##  V Extraction of field data at the particle

I use a straightforward one-sided extrapolation of field values to the right
of the particle’s position to extract values for $\psi$ and
$\partial_{r^{*}}\psi$. Specifically I fit a fourth order polynomial

$p(x)=\sum_{n=0}^{4}\frac{c_{i}}{n!}x^{n}\text{,}$ (38)

where $x=r^{*}-r^{*}_{0}$ to the five points to the right of the particle’s
current position and extract $\psi$ and $\partial_{r^{*}}\psi$ as $c_{0}$ and
$c_{1}$, respectively. In order to calculate
$\frac{\partial\psi(t_{0},r^{*}_{0})}{\partial t}$ I follow Sago  and
calculate $\frac{\mathrm{d}\psi\bm{(}t,r^{*}(t)\bm{)}}{\mathrm{d}t}$ on the
world line of the particle. Since this can be calculated using either the
field values on the world line

$\frac{\mathrm{d}\psi\bm{(}t,r^{*}(t)\bm{)}}{\mathrm{d}t}=\\\
\mbox{}\frac{\psi\bm{(}t+h,r^{*}(t+h)\bm{)}-\psi\bm{(}t-h,r^{*}(t-h)\bm{)}}{2h}+O(h^{2})\text{,}$
(39)

or as

$\frac{\mathrm{d}\psi\bm{(}t,r^{*}(t)\bm{)}}{\mathrm{d}t}=\frac{\partial\psi}{\partial
t}+\frac{\partial\psi}{\partial
r^{*}}\frac{\mathrm{d}r^{*}_{0}}{\mathrm{d}t}\text{,}$ (40)

where both $\frac{\partial\psi}{\partial r^{*}}$ and
$\frac{\mathrm{d}r^{*}_{0}}{\mathrm{d}t}=\dot{r}_{0}/E$ are known, this allows
me to find

$\frac{\partial\psi}{\partial
t}=\frac{\mathrm{d}\psi\bm{(}t,r^{*}(t)\bm{)}}{\mathrm{d}t}-\frac{\partial\psi}{\partial
r^{*}}\frac{\mathrm{d}r^{*}_{0}}{\mathrm{d}t}\text{.}$ (41)

I repeat this procedure to the left of the particle. As a check for the
extraction procedure, I compare the difference between the right hand and left
hand values $\left[\psi\right]=\psi_{\text{right}}-\psi_{\text{left}}$ with
the analytically calculated jump conditions of appendix D.1. Similarly I check
whether the numerical solutions obtained for $\chi$ and $\xi$ directly are
consistent with Eqs. (20) and (19), which give them in terms of derivatives of
$\psi$.

##  VI Regularization of the mode sum

The regularization procedure operates on scalar spherical harmonic modes of
the multipole coefficients $F^{\ell m}_{(\mu)(\nu)}$ of the Faraday tensor. As
a first step I use the auxiliary fields $\psi$, $\chi$ and $\xi$ to
reconstruct

$\displaystyle
A^{\ell^{\prime}m^{\prime}}_{r,t}-A^{\ell^{\prime}m^{\prime}}_{t,r}=\frac{\psi}{r^{2}}\text{,}$
(42a)
$\displaystyle\partial_{t}v^{\ell^{\prime}m^{\prime}}-A^{\ell^{\prime}m^{\prime}}_{t}=-\xi\text{,}$
(42b) and
$\displaystyle\partial_{r}v^{\ell^{\prime}m^{\prime}}-A^{\ell^{\prime}m^{\prime}}_{r}=-\frac{\chi}{f}\text{,}$
(42c)

the combinations of the vector potential modes needed to obtain the even
sector of a tensor spherical harmonic decomposition of the Faraday tensor. The
auxiliary field $\tilde{v}$ and its derivatives provide the odd sector of the
decomposition.

Using the comples pseudo-Cartesian tetrad $e^{\alpha}_{\ (0)}$, $e^{\alpha}_{\
(\pm)}$ and $e^{\alpha}_{\ (3)}$ introduced in paper I, I define tetrad
components

$F^{\text{ret}}_{(\mu)(\nu)}\equiv F^{\text{ret}}_{\alpha\beta}e^{\alpha}_{\
(\mu)}e^{\beta}_{\ (\nu)}$ (43)

of the Faraday tensor.

I construct the spherical harmonic modes of $F^{\text{ret}}_{(\mu)(\nu)}$
using the coupling coefficients displayed in Eq. (77).

$F^{\ell
m,\text{ret}}_{(\mu)(\nu)}=\sum_{\ell^{\prime},m^{\prime}}\left[C^{ab}_{(\mu)(\nu)}(\ell^{\prime}m^{\prime}|\ell
m)\left(A^{\ell^{\prime}m^{\prime}}_{b,a}-A^{\ell^{\prime}m^{\prime}}_{a,b}\right)\right.\\\
\left.+D^{a}_{(\mu)(\nu)}(\ell^{\prime}m^{\prime}|\ell
m)\left(\partial_{a}v^{\ell^{\prime}m^{\prime}}-A^{\ell^{\prime}m^{\prime}}_{a}\right)\right.\\\
\left.+E^{a}_{(\mu)(\nu)}(\ell^{\prime}m^{\prime}|\ell
m)\partial_{a}\tilde{v}^{\ell^{\prime}m^{\prime}}\right.\\\
\left.+E_{(\mu)(\nu)}(\ell^{\prime}m^{\prime}|\ell
m)\tilde{v}^{\ell^{\prime}m^{\prime}}\right]$ (44)

I calculate the multipole coefficients of $F^{\ell,\text{ret}}_{(\mu)(\nu)}$
as

$F^{\ell,\text{ret}}_{(\mu)(\nu)}=\sum_{m}F^{\ell
m,\text{ret}}_{(\mu)(\nu)}(t,r_{0})Y_{\ell
m}(\frac{\pi}{2},\varphi_{0})\text{,}$ (45)

and regularize them as in Eq. (5).

$F^{R}_{(\mu)(\nu)}=\sum_{\ell}\biggl{\\{}F^{\ell,\text{ret}}_{(\mu)(\nu)}-q\Bigl{[}A_{(\mu)(\nu)}\Bigl{(}\ell+\frac{1}{2}\Bigr{)}+B_{(\mu)(\nu)}+\\\
\frac{C_{(\mu)(\nu)}}{\ell+\frac{1}{2}}+\frac{D_{(\mu)(\nu)}}{(\ell-\frac{1}{2})(\ell+\frac{3}{2})}\Bigr{]}\biggr{\\}}$
(46)

I calculate the regularized self-force using
$F^{R}_{(\mu)}=qF^{R}_{(\mu)(\nu)}u^{(\nu)}$. Finally I reconstruct the vector
components of the self-force by from the tetrad components

$\displaystyle F^{R}_{t}=\sqrt{f_{0}}F^{R}_{(0)}\text{,}$ (47a) $\displaystyle
F^{R}_{r}=\frac{1}{\sqrt{f_{0}}}\operatorname{Re}\left(F^{R}_{(+)}e^{-i\varphi_{0}}\right)\text{,}$
(47b) $\displaystyle
F^{R}_{\phi}=r_{0}\operatorname{Im}\left(F^{R}_{(+)}e^{-i\varphi_{0}}\right)\text{.}$
(47c)

##  VII Numerical tests

In this section I present the tests I performed to validate my numerical
evolution code. I performed the same set of tests as described in paper II.
First, in order to check the second-order convergence rate of the code, I
performed regression runs with increasing resolution. As a second test, I
computed the regularized self-force for several different combinations of
orbital elements $p$ and $e$ and checked that the multipole coefficients decay
with $\ell$ as expected. This provided a very sensitive check on the overall
implementation of the numerical scheme as well as the analytical calculations
that lead to the regularization parameters.

###  VII.1 Convergence tests

Convergence tests are a straightforward way to test the implementation of a
numerical scheme. I performed regression runs for my second-order convergent
code using a non-zero charge $q$ and an eccentric orbit. I extract the field
at the position of the particle, and thus also test the implementation of the
extraction algorithm described in section V. I choose the $\ell=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. 4 the convergence is approximately of second
order.

Figure 4: Convergence test of the numerical algorithm in the sourced case. I
show differences between simulations using different step sizes of 16, 32 and
64 cells per $M$. Displayed are the rescaled differences
$\delta_{32-16}=\xi(h=1/32M)-\xi(h=1/16M)$ etc. of the field values at the
position of the particle for a simulation with $\ell=6$, $m=4$ and $p=7$,
$e=0.3$. I see that the convergence is approximately second-order. The curves
are rescaled in such a way as to provide an estimate for the error of the
highest resolution run compared to the real ($h\equiv 0$) solution.

In the region $150\,M\lesssim t\lesssim 400\,M$ the two curves lie on top of
each other, as expected for a second-order convergent algorithm. In the region
from $400\,M$ to $450\,M$ there is some difference between the two lines,
caused by cell crossing effects similar to those discussed in paper II.

###  VII.2 Discontinuity across the world line

The singular source term on the right hand side of Eqs. (21b) – (21f) implies
that the fields $\psi$, $\chi$ and $\xi$ are discontinuous across the world
line. Since the jump conditions can be calculated analytically as done in
appendix D.1, I can check whether the numerical results faithfully reproduce
the expected behaviour. Using the methods described in section V I obtain one-
sided extrapolation for the field values and their spatial derivatives. For
the highest resolution run used in the regression analysis in section VII.1 I
find that the numerical results for $\xi$ agree with the analytical
calculation of the jump condition up to terms of the order of $10^{-8}$; two
orders of magnitude smaller than the estimated numerical error of $10^{-6}$.
For $\partial_{r^{*}}\xi$ the situation is reversed, with the numerical error
in the jump condition being about an order of magnitude larger than the
numerical error in the field derivative itself. The accuracy of the numerical
derivatives is therefore limited by the accuracy of the extraction scheme,
resulting in about three significant figures for the set of parameters
displayed in Fig. 4. However the regularization calculation is constructed in
such a way that no derivatives of the fields need be obtained in order to
calculate the self-force. I therefore feel that I can accept the reduced
accuracy provided by the simple extraction scheme.

###  VII.3 High-$\ell$ behaviour of the multipole coefficients

Inspection of Eq. (5) reveals that a plot of $F^{\ell}_{(\mu)(\nu)}$ as a
function of $\ell$ (for a fixed value of $t$) should display a linear growth
in $\ell$ for large $\ell$. Removing the $A_{(\mu)(\nu)}$ term should produce
a constant curve, removing the $B_{(\mu)(\nu)}$ term (given that
$C_{(\mu)(\nu)}=0$) should produce a curve that decays as $\ell^{-2}$, and
finally, removing the $D_{(\mu)(\nu)}$ term should produce a curve that decays
as $\ell^{-4}$. It is a powerful test of the overall implementation to check
whether the numerical data behaves as expected. Fig. 5 plots the remainders as
obtained from my numerical simulation, demonstrating the expected behaviour.

Figure 5: Multipole coefficients of the dimensionless Faraday tensor component
$\frac{M^{2}}{q}\operatorname{Im}F^{R}_{(+)(-)}$ for a particle on an
eccentric orbit ($p=7.2$, $e=0.5$). The coefficients are extracted at
$t=500\,M$ along the trajectory shown in Fig. 6. The plots show several stages
of the regularization procedure, with a closer description of the curves to be
found in the text. A uniform stepsize of $h=1/512\,M$ was used.

It displays, on a logarithmic scale, the absolute value of
$\operatorname{Im}F^{\ell,R}_{(+)(-)}$, the imaginary part of the
$F^{R}_{(+)(-)}$ tetrad component of the Faraday tensor. The orbit is
eccentric ($p=7.2$, $e=0.5$), and all components of the self-force require
regularization. The first curve (in triangles) shows the unregularized
multipole coefficients that increase linearly in $\ell$, as confirmed by
fitting a straight line to the data. The second curve (in squares) shows
partially regularized coefficients, obtained after the removal of
$(\ell+1/2)A_{(\mu)(\nu)}$; this clearly approaches a constant for large
values of $\ell$. The curve made up of diamonds shows the behaviour after
removal of $B_{(\mu)(\nu)}$; because $C_{(\mu)(\nu)}=0$, it decays as
$\ell^{-2}$, a behaviour that is confirmed by a fit to the $\ell\geq 5$ part
of the curve. Finally, after removal of
$D_{(\mu)(\nu)}/[(\ell-\frac{1}{2})\,(\ell+\frac{3}{2})]$ the terms of the sum
decrease in magnitude approximately as $\ell^{-4}$ when fitting to the data
points $\ell\geq 7$. This result depends slightly on the range of points used
for the fit. I expect this to be due to the fact that I stop at $\ell=15$,
which seems to be not large enough to show the asymptotic behaviour. Extending
the range to very high values of $\ell$ proved to be very difficult, since the
numerical code is only second order convergent, so that the numerical errors
become dominant by the time the asymptotic behaviour begins to show.

Each one of the last two curves would result in a converging sum, but the
convergence is faster after subtracting the $D_{(\mu)(\nu)}$ terms. I thereby
gain about one order of magnitude in the accuracy of the estimated sum.

Figure 5 provides a sensitive test of the implementation of both the numerical
and analytical parts of the calculation. Small mistakes in either one will
cause the difference in Eq. (5) to have a vastly different behaviour.

###  VII.4 Accuracy of the numerical method

In this work I are less demanding with the numerical accuracy then I were in
paper II, where I describe a very high accuracy numerical code. Implementing
suach a code is very tedious even for the scalar case, and much more so for
the electromagnetic case treated here. Therefore I implement a simpler method
that allows me to access the physics of the problem without being hindered by
technical problems due to a complicated numerical method.

An estimate for the truncation error arising from cutting short the summation
in Eq. (5) at some $\ell_{\text{max}}$ can be calculated by considering the
behaviour of the remaining terms for large $\ell$. Detweiler et. al. Detweiler
et al. (2003) showed that the remaining terms scale as $\ell^{-4}$ for large
$\ell$. They find the functional form of the terms to be

$\frac{E\mathcal{P}_{3/2}}{(2\ell-3)(2\ell-1)(2\ell+3)(2\ell+5)}\text{,}$ (48)

where $\mathcal{P}_{3/2}=36\sqrt{2}$. I fit a function of this form to the
tail end of a plot of the multipole coefficients to find the coefficient $E$
in Eq. (48). Extrapolating to $\ell\rightarrow\infty$ I find that the
truncation error is

$\displaystyle\varepsilon$
$\displaystyle=\sum_{\ell=\ell_{\text{max}}}^{\infty}[\text{Eq.~{}\eqref{eqn:Eterm}}]$
(49)
$\displaystyle=\frac{12\sqrt{2}E\ell_{\text{max}}}{(2\ell_{\text{max}}+3)(2\ell_{\text{max}}+1)(2\ell_{\text{max}}-1)(2\ell_{\text{max}}-3)}\text{,}$
(50)

where $\ell_{\text{max}}$ is the value at which I cut the summation short.

A second source of error lies in the numerical calculation 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
convergence order, I can estimate this discretization error by extrapolating
from simulations using different step sizes down to $h=0$. This is what was
done in the graphs shown in Sec. VII.1.

I display results for the mildly eccentric orbit shown in Fig. 6 with data
extracted at $t=500\,M$, that is at the instant shown in Fig. 5. At this
moment, the multipole coefficients of $\operatorname{Re}(F^{R}_{(+)})$ decay
as expected, but e.g. the $\operatorname{Im}(F^{R}_{(+)})$ component decays
faster with $\ell$ for the range of modes $0\leq\ell\leq 13$ modes that were
calculated. I choose an orbit of low eccentricity as high eccentricity causes
the field values to be plagued by high frequency noise, as discussed in paper
II. This makes it impossible to reliably estimate the discretization error for
these orbits.

Table 1 lists typical values for the errors discussed above.

error estimation | mildly eccentric orbit
---|---
relative truncation error in $\frac{M^{2}}{q^{2}}\operatorname{Re}(F^{R}_{(+)})$ | $2\times 10^{-4}$
relative discretization error in $\frac{M^{2}}{q}\psi$ | $\approx 10^{-7}$
Table 1: Estimated values for the various errors in the components of the
self-force as described in the text. I show the truncation and discretization
errors for the mildly eccentric orbit ($p=7.2$, $e=0.5$). The truncation error
is calculated using a plot similar to the one shown in Fig. 5. The
discretization error is estimated using a plot similar to that in Fig. 4 for
the $\ell=2$, $m=2$ mode.

##  VIII Sample results

In this section I describe some results obtained from my numerical
calculation.

###  VIII.1 Mildly eccentric orbit

I choose a particle on an eccentric orbit with $p=7.2$, $e=0.5$ which starts
at $r=pM/(1-e^{2})$, halfway between periastron and apastron. The field is
evolved for $600\,M$ with a uniform resolution of 512 grid points per $M$,
both in the $t$ and $r^{*}$ directions, for all values of $\ell$. Multipole
coefficients for $1\leq\ell\leq 15$ are calculated and used to reconstruct the
regularized self-force $F_{\alpha}$ along the geodesic. Figure 7 shows the
result of the calculation.

Figure 6: Trajectory of a particle with $p=7.2$, $e=0.5$. The cross-hair
indicates the point where the data for Fig. 5 was extracted. Figure 7:
Regularized dimensionless self-force $\frac{M^{2}}{q^{2}}F_{t}$,
$\frac{M^{2}}{q^{2}}F_{r}$ and $\frac{M}{q^{2}}F_{\phi}$ on a particle on an
eccentric orbit with $p=7.2$, $e=0.5$.

For the choice of parameters used to calculate the force shown in Fig. 7, the
error bars corresponding to the truncation error Eq. (49) (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, I see that the self-force is most
important when the particle is closest to the black hole (ie. for
$200\,M\lesssim t\lesssim 400\,M$). The self-force acting on the particle is
very small once the particle has moved away to $r\approx 15\,M$.

###  VIII.2 Zoom-whirl orbit

Particles on highly eccentric orbits are of most interest as sources of
gravitational radiation. For nearly parabolic orbits with $e\lesssim 1$ and
$p\gtrsim 6+2e$, a particle revolves around the black hole a number of times,
moving on a nearly circular 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 spacetime,
emitting copious amounts of radiation. Figures 8 and 9 show the trajectory of
a particle and the force on such an orbit with $p=7.8001$, $e=0.9$ calculated
using a uniform step size of $h=1/256\,$ throughout the range $1\leq\ell\leq
15$.

Figure 8: 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.
10 and 11 was extracted. Figure 9: Self-force acting on a particle. Shown is
the dimensionless self-force $\frac{M^{2}}{q^{2}}F_{t}$,
$\frac{M^{2}}{q^{2}}F_{r}$ and $\frac{M}{q^{2}}F_{\phi}$ on a zoom-whirl orbit
with $p=7.8001$, $e=0.9$. No error bars showing an estimate error are shown,
as the errors shown are to small to show up on the graph. Notice that the
self-force is essentially zero during the zoom phase $900\,M\lesssim t\lesssim
1200\,M$ and reaches a constant value very quickly after the particle enters
into the whirl phase.

Even more so than for the mildly eccentric orbit discussed in Sec. VIII.1, 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. The force
graph is very similar to that obtained for the scalar self-force in paper II,
however the overshooting behaviour at the onset and near the end of the whirl
phase is not as pronounced.

Since the rates of change in energy $E$ and angular momentum $J$ of the
trajectory are directly related to the self-force

$\dot{E}=-a_{t}\text{,}\qquad\dot{J}=a_{\phi}\text{,}$ (51)

it is easy to see that the self-force shown in Fig. 9 confirms the expectation
that the self-force decreases both the energy and angular momentum of the
particle while radiation is emitted.

In Fig. 10 and Fig. 11 I show plots of $F^{\ell}_{(0)}$ constructed from
$F^{\ell}_{(\mu)(\nu)}$ after the removal of the $A_{(\mu)(\nu)}$,
$B_{(\mu)(\nu)}$, and $D_{(\mu)(\nu)}$ terms.

Figure 10: Multipole coefficients of
$\frac{M^{2}}{q}\operatorname{Re}F^{R}_{(0)}$ for a particle on a zoom-whirl
orbit ($p=7.8001$, $e=0.9$). The coefficients are extracted at $t=525\,M$ when
the particle is deep within the whirl phase. Here $\dot{r}\approx 0$ and the
behaviour of $F^{R}_{(\mu),\ell}$ is very close to that for a circular orbit,
requiring very little regularization. Red triangles are used for the
unregularized multipole coefficients $F_{(0),\ell}$, squares, diamonds and
disks are used for the partly regularized coefficients after the removal of
the $A_{(0)}$, $B_{(0)}$ and $D_{(0)}$ terms respectively. Figure 11:
Multipole coefficients of $\frac{M^{2}}{q}\operatorname{Re}F^{R}_{(0)}$ for a
particle on a zoom-whirl orbit ($p=7.8001$, $e=0.9$). The coefficients are
extracted at $t=1100\,M$ when the particle is far away from the black hole. As
$\dot{r}$ is non-zero, all components of the self-force require regularization
and I see that the dependence of the multipole coefficients on $\ell$ is as
predicted by Eq. 5. After the removal of the regularization parameters
$A_{(\mu)(\nu)}$, $B_{(\mu)(\nu)}$, and $D_{(\mu)(\nu)}$ the remainder is
proportional to $\ell^{0}$, $\ell^{-2}$ and $\ell^{-4}$ respectively.

##  IX Effects of the conservative self-force

In this section only, I will use the subscript “0” to denote quantities
evaluated on the unperturbed geodesic, and no subscript to denote quantities
evaluated on the perturbed world-line.

I follow the literature (see e.g. Pound and Poisson (2008)) and define the
dissipative part to be the half retarded minus half advanced force and the
conservative part to be the half retarded plus half advanced force

$\displaystyle
F^{\text{diss}}_{\alpha}\equiv\frac{1}{2}\left(F^{\text{ret}}_{\alpha}-F^{\text{adv}}_{\alpha}\right)\text{,}$
(52) $\displaystyle
F^{\text{cons}}_{\alpha}\equiv\frac{1}{2}\left(F^{\text{ret}}_{\alpha}+F^{\text{adv}}_{\alpha}\right)\text{.}$
(53)

The conservative force is the time reversal invariant part of the self-force.
It does not affect the radiated energy or angular momentum fluxes $\dot{E}$
and $\dot{J}$; it shifts the values of $E$ and $J$ away from their geodesic
values, affecting the orbital motion and the phase of the emitted waves.

To obtain expressions for $E$ and $J$ under the influence of the self-force, I
employ the procedure described in Diaz-Rivera et al. (2004). I begin by
writing down the normalization condition for the four velocity

$-1=u^{\alpha}u_{\alpha}=-\frac{E^{2}}{f}+\frac{J^{2}}{r^{2}}\text{,}$ (54)

as well as the $r$-component of the geodesic equation

$\frac{F^{r}}{m}=\ddot{r}-\frac{M}{(r-2M)r}\dot{r}^{2}-\frac{(r-2M)J^{2}}{r^{4}}+\frac{ME^{2}}{(r-2M)r}\text{,}$
(55)

where $F^{r}=qF^{r}_{\ \mu}u^{\mu}$ is the radial component of the self-force.
Solving Eq. (54) and (55) I find

$\displaystyle E^{2}=E_{0}^{2}-\frac{(r-2M)r}{r-3M}\frac{F^{r}}{m}\text{,}$
(56) $\displaystyle J^{2}=J_{0}^{2}-\frac{r^{4}}{r-3M}\frac{F^{r}}{m}\text{,}$
(57)

where

$\displaystyle E_{0}^{2}$
$\displaystyle=\dot{r}^{2}+\frac{(r-2M)r\ddot{r}}{r-3M}+\frac{(r-2M)^{2}}{(r-3M)r}\text{,}$
(58) $\displaystyle J_{0}^{2}$
$\displaystyle=\frac{r^{4}\ddot{r}}{r-3M}+\frac{Mr^{2}}{r-3M}\text{.}$ (59)

I stress that $E_{0}$ and $J_{0}$ are not the geodesic _values_ for energy and
angular momentum. They are of the correct form but are evaluated using the
_accelerated_ values for $r$, $\dot{r}$ and $\ddot{r}$ (instead of the
geodesic values $r_{0}$, $\dot{r}_{0}$, etc.).

For small perturbing force of order $\varepsilon$ I expand Eqs. (56) and (57)
in terms of the perturbation strength and find

$\displaystyle E=E_{0}+\Delta E\approx
E_{0}-\varepsilon\frac{(r-2M)r}{2(r-3M)E_{0}}\frac{F^{r}}{m}+O(\varepsilon^{2})\text{,}$
(60) $\displaystyle J=J_{0}+\Delta J\approx
J_{0}-\varepsilon\frac{r^{4}}{2(r-3M)J_{0}}\frac{F^{r}}{m}+O(\varepsilon^{2})\text{,}$
(61)

where $F^{r}$ is evaluated with the help of the unperturbed four velocity
$u_{0}^{\alpha}=[E_{0}/f,\dot{r}_{0},0,J_{0}/r_{0}^{2}]$. The fractional
changes $\Delta E/E_{0}$ and $\Delta J/J_{0}$ are given by

$\displaystyle\Delta
E/E_{0}=-\varepsilon\frac{(r-2M)r}{2(r-3M)E_{0}^{2}}\frac{F^{r}}{m}+O(\varepsilon^{2})\text{,}$
(62) $\displaystyle\Delta
J/J_{0}=-\varepsilon\frac{r^{4}}{2(r-3M)J_{0}^{2}}\frac{F^{r}}{m}+O(\varepsilon^{2})\text{.}$
(63)

Once the perturbations in $E$ and $J$ are known, I calculate the change in the
angular frequency

$\Omega\equiv\frac{\mathrm{d}\varphi}{\mathrm{d}t}=\frac{r-2M}{r^{3}}\frac{J}{E}\text{.}$
(64)

For small perturbing force I expand in powers of the perturbation strength

$\Omega=\frac{r_{0}-2M}{r_{0}^{3}}\frac{J_{0}}{E_{0}}\biggl{[}1-\varepsilon\biggl{(}\frac{r^{4}}{2(r-3M)J_{0}^{2}}\\\
-\frac{(r-2M)r}{2(r-3M)E_{0}^{2}}\biggr{)}\frac{F^{r}}{m}\biggr{]}+O(\varepsilon^{2})\text{.}$
(65)

The relative change $\Delta\Omega/\Omega_{0}$ is given by

$\Delta\Omega/\Omega_{0}=-\varepsilon\biggl{(}\frac{r^{4}}{2(r-3M)J_{0}^{2}}-\frac{(r-2M)r}{2(r-3M)E_{0}^{2}}\biggr{)}\frac{F^{r}}{m}\\\
+O(\varepsilon^{2})\text{.}$ (66)

###  IX.1 Circular orbits

The effect of the conservative self-force is most clearly observed for
circular orbits, where the unperturbed angular frequency $\Omega$ as well as
the shift due to the perturbation are constant in time.

For a particle in circular motion the self-force is constant in time and it
turns out that the radial component is entirely conservative whereas the $t$
and $\phi$ components are entirely dissipative. For circular orbits, the
unperturbed values of $E$ and $J$ are given by

$\displaystyle E_{0}$
$\displaystyle=\frac{r_{0}-2M}{\sqrt{r_{0}(r_{0}-3M)}}\text{,}$ (67)
$\displaystyle J_{0}$ $\displaystyle=r_{0}\sqrt{\frac{M}{r_{0}-3M}}\text{,}$
(68)

and substituting these into Eq. (65) I find

$\Omega=\sqrt{\frac{M}{r_{0}^{3}}}-\frac{(r_{0}-3M)}{2mM}\sqrt{\frac{M}{r_{0}}}F_{r}+O(\varepsilon^{2})\text{,}$
(69)

where the first term is just the angular frequency for an unperturbed geodesic
at radius $r_{0}$. The fractional change $\Delta\Omega/\Omega_{0}$ is then

$\frac{\Delta\Omega}{\Omega_{0}}=-\frac{(r_{0}-3M)r_{0}}{2mM}F_{r}+O(\varepsilon^{2})\text{.}$
(70)

Similarly the fractional changes in $E$ and $J$: $\Delta E/E_{0}$ and $\Delta
J/J_{0}$ are given by

$\displaystyle\Delta E/E_{0}$ $\displaystyle=-\frac{r_{0}}{2m}F_{r}\text{,}$
(71) $\displaystyle\Delta J/J_{0}$
$\displaystyle=-\frac{(r_{0}-2M)r_{0}}{2mM}F_{r}\text{.}$ (72)

Figure 12: Fractional change $\Delta\Omega/\Omega_{0}$ induced by the presence
of the conservative self-force. The effect of the self-force is to move the
radius of the orbit outward, decreasing its angular frequency.

Figure 12 shows the fractional change in $\Omega_{0}$, $E$ and $J$ as a
function of the orbit’s radius $r_{0}$.

###  IX.2 Eccentric orbits

For eccentric orbits the self-force is no longer constant in time and I have
to numerically calculate both the retarded and the advanced self-force in
order to construct the conservative self-force. I find the advanced force by
running the simulation backwards in time. That is I start the evolution on the
very last time slice and evolve backwards in time until I reach the slice
corresponding to $t=0$. I reverse the boundary condition at the event horizon
to be be outgoing radiation only $(\partial_{t}+\partial_{r^{*}})\psi=0$ and
adjust the outer boundary so as to simulate only the backwards domain of
dependence of the initial slice. I do not change the trajectory of the
particle. I do not change the regularization parameters, since they depend
only on the local behaviour of the field and are insensitive to the boundary
conditions far away.

####  IX.2.1 Conservative force on zoom-whirl orbits

I calculate the conservative self-force on a zoom-whirl orbit with $p=7.8001$,
$e=0.9$. Figs. 13 and 14

Figure 13: $r$ component of the dimensionless self-force acting on a particle
on a zoom-whirl orbit ($p=7.8001$, $e=0.9$) around a Schwarzschild black hole.
Shown are the retarded (solid, red), advanced (dashed, green), conservative
(dotted, blue) and dissipative (finely dotted, pink) force acting on the
particle. Figure 14: $\varphi$ component of the dimensionless self-force
acting on a particle on a zoom-whirl orbit ($p=7.8001$, $e=0.9$) around a
Schwarzschild black hole. Shown are the retarded (solid, red), advanced
(dashed, green), conservative (dotted, blue) and dissipative (finely dotted,
pink) force acting on the particle.

display the breakdown of the self-force into retarded and advanced, and
conservative and dissipative parts for a particle on a zoom-whirl orbit. In
both plots the force is very weak when the particle is in the zoom phase
$t\lesssim 400\,M$ or $t\gtrsim 800\,M$ and nearly constant while the particle
is in the whirl phase $400\,M\lesssim t\lesssim 800\,M$. Inspection of the
behaviour of the $r$ component reveals that it is almost exclusively
conservative, with only a tiny dissipative effect when the particle enters or
leaves the whirl phase. This result is consistent with the observation that
the particle moves on a nearly circular trajectory while in the whirl phase,
for which the radial component is precisely conservative. Similarly the $\phi$
component is almost entirely dissipative, with only a small conservative
contribution when the particle enters or leaves the whirl phase, its maximum
coinciding with that of $\ddot{r}$ (not shown on the graph).

I calculate the relative changes in $E$, $J$ and $\Omega$ under the influence
of the self-force using Eqs. (62), (63), (66). Fig. 15

Figure 15: Relative change in $\Omega$, $E$, $J$ for a particle on a zoom-
whirl orbit due to the conservative electromagnetic self-force.

displays the relative changes $\Delta E/E_{0}$, $\Delta J/J_{0}$ and
$\Delta\Omega/\Omega_{0}$ for a particle on a zoom whirl orbit $p=7.8001$,
$e=0.9$. The change in $E$, $J$ and $\Omega$ is strongest in the whirl phase
when $r\approx 4.1M$. It is consistent with the shift experienced by a
particle on a circular orbit at $4.1M$.

###  IX.3 Effects on the innermost stable orbit

In the gravitational case, considerable work has been done to identify gauge
invariant effects of the self-force Detweiler (2008a); Sago et al. (2008). The
electromagnetic self-force is not subject to the same ambiguity thus it can
help shed light on the gravitational case as well by providing a clear
distinction between kinetic and dynamic effects. In this section I calculate
the effect of the conservative self-foce on the location of the innermost
stable circular orbit around a Schwarzschild black hole. Such a calculation
was first performed for the scalar self-force by Diaz-Rivera et al. (2004),
where a highly accurate frequency domain numerical scheme was used. Recently
Barack and Sago (2009, 2011) have extended this calculation to gravity, using
their time domain code to perform the intergration of the wave equation. Since
the code presented in this paper is in the time domain as well, it is closest
in spirit to Barack and Sago (2009).

##  X Retardation of the self-force

For scalar perturbation in a weak gravitational field Poisson Pfenning and
Poisson (2002) showed that the self-force is delayed with respect to the
particle motion by the light travel time from the particle to the central body
and back to the particle again. In a spacetime where the central body is
compact the treatment of Pfenning and Poisson (2002) is no longer directly
applicable, but I still expect some retardation in the self-force when
compared to the particle’s motion. To study this effect, I calculate the self-
force on an eccentric orbit with $p=78$, $e=0.9$; ten times larger than the
zoom-whirl orbit discussed earlier. The large orbit was chosen so as to be
able to clearly see any possible retardation which might not be visible if the
particle’s orbit is deep within the strong field region close to the black
hole.

Figure 16: $r$ component of the dimensionless self-force acting on a particle
on an orbit with $p=78$, $e=0.9$. Shown are the retarded and advanced forces
as well as $\dot{r}$. The vertical line at $t\approx 2383\,M$ marks the time
of closest approach to the black hole. Figure 17: $\phi$ component of the
dimensionless self-force acting on a particle on an orbit with $p=78$,
$e=0.9$. Shown are the retarded and advanced forces as well as $\dot{r}$. The
vertical line at $t\approx 2383\,M$ marks the time of closest approach to the
black hole.

Figures 16 and 17 display plots of the $r$ and $\phi$ components of the self-
force acting on the particle close to periastron. Shown are the retarded and
advanced forces as well as the particle’s radial velocity $\dot{r}$. Without
considering retardation I expect the self-force to be strongest when the
particle is closest to the black hole, when $\dot{r}=0$, as evident in Fig.
13. Clearly for the $r$ component displayed in Fig. 16 the retarded and
advanced forces both peak at a time very close to the zero crossing of
$\dot{r}$, suggesting very little time delay in the $r$ component of the self-
force. In Fig. 17 on the other hand the retarded and advanced $\phi$-component
of the self-force peaks away from the time of closest approach
$t_{\text{min}}$. Inspection of the graph shows that the delay (advance)
between the time of closest approach and the peak in the retarded (advanced)
force is compatible with a delay of $\Delta t_{\text{min}}\approx
2(r_{\text{min}}-3.0\,M)\approx 74\,M$. Using a delay of $\Delta t\approx
2[r_{0}(t)-3.0\,M]$ and plotting $F^{\text{ret}}_{\varphi}(t+\Delta t)$ and
$-F^{\text{adv}}_{\varphi}(t-\Delta t)$ versus $t$ both curves visually lie on
top of each other and the maximum is located at $t_{\text{min}}$ as shown in
Fig. 18 below. This suggests that the self-force is in large parts due to
radiation that travels into the strong field region close to the black hole
and is scattered back to the particle. The time delay can then be loosely
interpreted as the time it takes the signal to travel to the light ring around
the black hole and back to the particle. This interpretation is loose for two
reasons: First $r^{*}$ and not $r$ is associated with the light travel time.
Using $r^{*}$, however, does not lead to a better overlap of the curves once a
suitable constant offset chosen. Second, for the zoom-whirl orbit shown in
Fig. 9 the (shallow) maximum in the self-force is offset by only $\Delta
t\approx 2[r_{0}(t)-1.0\,M]$ which leads to a reasonable overlap of the two
curves. Interestingly using $r^{*}$ instead of $r$ yields a worse overlap. For
very large orbits $p=780$, $e=0.9$ it is impossible to read off the small
constant offset to the dominant $2r_{0}(t)$ contribution.

##  XI Weak field limit

As a last application I use my code to compare the numerical self-force in the
weak field region to the self-force calculated using the weak field expression

$\bm{f}_{\text{self}}=\lambda_{c}\frac{q^{2}}{m}\frac{M}{r^{3}}\hat{\bm{r}}+\lambda_{rr}\frac{2}{3}\frac{q^{2}}{m}\frac{\mathrm{d}\bm{g}}{\mathrm{d}t}\text{,}\qquad\bm{g}=-\frac{M}{r^{2}}\hat{\bm{r}}\text{,}$
(73)

of Pound and Poisson (2008); DeWitt-Morette and Dewitt (1964). I calculate the
self-force for a particle on an eccentric orbits with $e=0.9$ and $p=78$ or
$p=780$. Fig. 18

Figure 18: $\phi$ component of the retarded (solid red line) and (dashed green
line) negative advanced self-forces acting on a particle with $p=78$, $e=0.9$.
The forces have been shifted by $\Delta t\approx 2[r_{0}(t)-3.0\,M]$. Also
shown is the self-force calculated using the weak field expression Eq. 73
(blue dotted line).

shows the retarded and (negative) advanced forces shifted by $\Delta t\approx
2[r_{0}(t)-3.0\,M]$ as well as the analytic force calculated using Eq. (73).
At this distance there are still some differences between the (shifted)
retarded field and the weak field expression. One reason for this lies in the
choice of a suitable $r$ coordinate to correspond to the $r$ coordinate in the
weak field expression. In this work I use the areal Schwarzschild $r$, but the
isotropic coordinate $\bar{r}=\frac{r-M}{2}+\frac{\sqrt{r\,(r-2*M)}}{2}$ or
even the tortoise $r^{*}$ could be used as well. Neither one yields a good
agreement between the two curves.

For $p=780$ using a shift of $\Delta t=2r_{0}(t)$ the agreement between
numerical data and analytic expression is excellent as is evident in Fig. 19.
At this distance $r$, $\bar{r}$ and $r^{*}$ are indistinguishable.

Figure 19: $\phi$ component of the retarded self-force acting on a particle on
an orbit with $p=780$, $e=0.9$ close to periastron. Shown are the numerical
(solid, red) and shifted analytical (dashed, green) forces. The agreement
between numerical and analytical calculation is excellent, the discrepancy for
$t\lesssim 7500\,M$ is due to initial data contamination.

##  XII Conclusions

I calculated the self-force acting a on an electromagnetic point charge in
orbit around a Schwarzschild black hole. To do so I calculated the
regularization parameters $A$,$B$, and $D$ in section B and implemented a
second order accurate numerical scheme in section II.

I find the behaviour of the electromagnetic self-force to be similar but not
identical to that of the scalar self-force. In both cases the self-force is
strongest when the particle is closest to the black hole. Further, during the
whirl phase of a zoom-whirl orbit with its nearly constant radius, the self-
force is very close to that of a particle in circular orbit at this radius. On
the other hand, the overshooting effect upon entering the whirl phase which
was observed in the scalar case is much weaker in the electromagnetic case.

I calculated the effects of the conservative self-force on circular orbits,
where it reduces the angular frequency and thus affects the phasing of the
observed waves. I find this effect to be much stronger in the electromagnetic
case than in the scalar case discussed in Diaz-Rivera et al. (2004). In
particular during the nearly circular whirl phase of a zoom-whirl orbit I find
that $\Omega$ decreases by $\approx 0.06\frac{q^{2}}{Mm}$. Due to the
smallness of the ratio $\frac{q^{2}}{Mm}$ this change is tiny for one orbit,
however since it accumulates over the inspiral, its effect on the total phase
shift during the full inspiral can be of order unity. This statement is not
directly transferable to the gravitational case since the radius $r_{0}$ of
the orbit is not a gauge invariant quantity. Therefore I cannot distinguish
between changes in $\Omega$ due to effects of the self-force and due to gauge
choices. To obtain a meaningful measure of the effect of the gravitational
self-force I need to compare two gauge invariant quantities, e.g. $\Omega$ and
the gauge invariant $u_{t}$ of Detweiler (2008b).

I investigated the retardation of the self-force with respect to the motion of
the particle. I found that the retardation is very weak for the $r$ component
of the force and strong in the $t$ and $\varphi$ components, which are linked
to radiated energy and angular momentum. In the later cases the retardation is
compatible with a delay of $\Delta t\approx 2(r_{0}(t)-R_{\text{delay}})$,
where $R_{\text{delay}}$ is a constant depending on the particle’s orbit.

###### Acknowledgements.

I thank Eric Poisson and Steve Detweiler for useful discussions and
suggestions. I gratefully acknowledge support by the Natural Sciences and
Engineering Council of Canada. This work was supported in part by NSF grants
PHY-0903973 and PHY-0904015. This work was made possible by the facilities of
the Shared Hierarchical Academic Research Computing Network
(SHARCNET:www.sharcnet.ca) as well as the e FoRCE cluster at Georgia Tech.

##  Appendix A Translation tables

I require coupling coefficients to translate between the tensor harmonic modes
of the Faraday tensor and the scalar harmonic modes of the tetrad components
of the Faraday tensor.

As a first step, I reconstruct the Faraday tensor modes from the numerical
variables. For the even mode auxiliary fields $\psi$, $\chi$ and $\xi$ this
reconstruction can be done algebraically while the odd sector requires a
numerical differentiation of the numerical variable $\tilde{v}$. The
reconstruction relations were already displayed in Eqs. (15) – (18), which
involves both the even and odd modes.

In terms of the vector potential the Faraday tensor modes are reconstructed
using the defining equation Eq. (7). In this case, the reconstruction of the
Faraday tensor reads

$\displaystyle F_{tr}$ $\displaystyle=\sum_{\ell,m}(A^{\ell m}_{r,t}-A^{\ell
m}_{t,r})\,Y^{\ell m}\text{,}$ (74a) $\displaystyle F_{tA}$
$\displaystyle=\sum_{\ell,m}[(v^{\ell m}_{,t}-A^{\ell m}_{t})\,Z^{\ell
m}_{A}+\tilde{v}^{\ell m}_{,t}\,X^{\ell m}_{A}]\text{,}$ (74b) $\displaystyle
F_{rA}$ $\displaystyle=\sum_{\ell,m}[(v^{\ell m}_{,r}-A^{\ell m}_{r})\,Z^{\ell
m}_{A}+\tilde{v}^{\ell m}_{,r}\,X^{\ell m}_{A}]\text{,}$ (74c) $\displaystyle
F_{\theta\varphi}$ $\displaystyle=\sum_{\ell,m}\tilde{v}_{\ell m}\,(X^{\ell
m}_{\varphi,\theta}-X^{\ell m}_{\theta,\varphi})\text{.}$ (74d)

Clearly both the expansion Eqs. (15) – (18) and the one in Eqs. (74a) – (74d)
are of the same form and it is only necessary to obtain one set of translation
coefficients to handle both the calculation using $\psi$, $\chi$ and $\xi$ in
the main text and the one using the vector potential that will be presented in
appendix C.

The tetrad components $F_{(\mu)(\nu)}$ are decomposed in terms of scalar
spherical harmonics

$F_{(\mu)(\nu)}=\sum_{\ell,m}F^{\ell m}_{(\mu)(\nu)}Y_{\ell m}\text{,}$ (75)

where each mode is given by

$F^{\ell m}_{(\mu)(\nu)}=\int F_{(\mu)(\nu)}\bar{Y}_{\ell
m}\,\mathrm{d}\Omega\text{.}$ (76)

To obtain expressions for the coupling coefficients I substitute
$F_{(\mu)(\nu)}=F_{\alpha\beta}e^{\alpha}_{\ (\mu)}e^{\beta}_{\ (\nu)}$ into
Eq. (76)

$\displaystyle F^{\ell m}_{(\mu)(\nu)}$
$\displaystyle=\int\mathrm{d}\Omega\,F_{(\mu)(\nu)}\bar{Y}^{\ell m}$
$\displaystyle=\int\mathrm{d}\Omega\,\left(A_{\beta,\alpha}-A_{\alpha,\beta}\right)e^{\alpha}_{\
(\mu)}e^{\beta}_{\ (\nu)}\bar{Y}^{\ell m}$
$\displaystyle=\int\sum_{\ell^{\prime},m^{\prime}}\Bigl{[}\left(A_{b,a}-A_{a,b}\right)e^{a}_{\
(\mu)}e^{b}_{\ (\nu)}\bar{Y}^{\ell m}$
$\displaystyle\quad+\left(A_{b,A}-A_{A,b}\right)e^{A}_{\ (\mu)}e^{b}_{\
(\nu)}\bar{Y}^{\ell m}$
$\displaystyle\quad+\left(A_{B,a}-A_{a,B}\right)e^{a}_{\ (\mu)}e^{B}_{\
(\nu)}\bar{Y}^{\ell m}$
$\displaystyle\quad+\left(A_{B,A}-A_{A,B}\right)e^{A}_{\ (\mu)}e^{B}_{\
(\nu)}\bar{Y}^{\ell m}\Bigr{]}\mathrm{d}\Omega$
$\displaystyle\equiv\sum_{\ell^{\prime},m^{\prime}}\left[C^{ab}_{(\mu)(\nu)}(\ell^{\prime}m^{\prime}|\ell
m)\left(A^{\ell^{\prime}m^{\prime}}_{b,a}-A^{\ell^{\prime}m^{\prime}}_{a,b}\right)\right.$
$\displaystyle\quad\left.+D^{a}_{(\mu)(\nu)}(\ell^{\prime}m^{\prime}|\ell
m)\left(\partial_{a}v^{\ell^{\prime}m^{\prime}}-A^{\ell^{\prime}m^{\prime}}_{a}\right)\right.$
$\displaystyle\quad\left.+E^{a}_{(\mu)(\nu)}(\ell^{\prime}m^{\prime}|\ell
m)\partial_{a}\tilde{v}^{\ell^{\prime}m^{\prime}}\right.$
$\displaystyle\quad\left.+E_{(\mu)(\nu)}(\ell^{\prime}m^{\prime}|\ell
m)\tilde{v}^{\ell^{\prime}m^{\prime}}\right]\text{,}$ (77)

which defines the coupling coefficients. It is often possible to express these
coupling coefficients in terms of linear combinations of the coupling
coefficients derived in paper I for the scalar field.

To simplify the notation of the coupling coefficients I use

$\displaystyle\gamma^{\ell m}$
$\displaystyle=\sqrt{\frac{(\ell+m)(\ell+m+1)}{(2\ell+1)(2\ell+3)}}\text{,}$
(78) $\displaystyle\epsilon^{\ell m}$
$\displaystyle=\sqrt{\frac{(\ell+m+1)(\ell-m+1)}{(2\ell+1)(2\ell+3)}}\text{,}$
(79)

as shorthands for recurring combinations of terms. With these the reusable
scalar coupling coefficients are written as

$\displaystyle C^{r}_{(+)}(\ell^{\prime}m^{\prime}|\ell m)$
$\displaystyle=-\gamma^{\ell-1,m}\sqrt{f}\delta_{\ell^{\prime}\ell-1}\delta_{m^{\prime}m-1}$
$\displaystyle\quad+\gamma^{\ell,-m+1}\sqrt{f}\delta_{\ell^{\prime}\ell+1}\delta_{m^{\prime}m-1}$
(80) $\displaystyle C_{(+)}(\ell^{\prime}m^{\prime}|\ell m)$
$\displaystyle=\gamma^{\ell-1,m}\frac{\ell-1}{r}\delta_{\ell^{\prime}\ell-1}\delta_{m^{\prime}m-1}$
$\displaystyle\quad+\gamma^{\ell,-m+1}\frac{\ell+2}{r}\delta_{\ell^{\prime}\ell+1}\delta_{m^{\prime}m-1}\text{,}$
(81) $\displaystyle C_{(-)}(\ell^{\prime}m^{\prime}|\ell m)$
$\displaystyle=-\gamma^{\ell-1,-m}\frac{\ell-1}{r}\delta_{\ell^{\prime}\ell-1}\delta_{m^{\prime}m+1}$
$\displaystyle\quad-\gamma^{\ell,m+1}\frac{\ell+2}{r}\delta_{\ell^{\prime}\ell+1}\delta_{m^{\prime}m+1}\text{.}$
(82)

Similarly it proves useful to define lower order coupling coefficients for the
odd sector, which is absent in the scalar case.

$\displaystyle E_{(+)}(\ell^{\prime}m^{\prime}|\ell
m)=-\frac{i}{r}\sqrt{(\ell-m+1)(\ell+m)}\delta_{\ell^{\prime}\ell}\delta_{m^{\prime}m-1}\text{,}$
(83) $\displaystyle E_{(-)}(\ell^{\prime}m^{\prime}|\ell
m)=-\frac{i}{r}\sqrt{(\ell+m+1)(\ell-m)}\delta_{\ell^{\prime}\ell}\delta_{m^{\prime}m+1}\text{.}$
(84)

In terms of the scalar coupling coefficients, the first coefficient for the
expansion of $F_{(0)(+)}$ is given by

$\displaystyle C^{tr}_{(0)(+)}(\ell^{\prime}m^{\prime}|\ell m)$
$\displaystyle=\int
Y^{\ell^{\prime}m^{\prime}}e^{t}_{(0)}e^{r}_{(+)}\bar{Y}^{\ell
m}\,\mathrm{d}\Omega$ $\displaystyle=\frac{1}{\sqrt{f}}\int
Y^{\ell^{\prime}m^{\prime}}e^{r}_{(+)}\bar{Y}^{\ell m}\,\mathrm{d}\Omega$ (85)
$\displaystyle=C^{r}_{(+)}(\ell^{\prime}m^{\prime}|\ell m)/\sqrt{f}\text{,}$

and all other combinations of $a$, $b$ and $(\mu)$, $(\nu)$ lead to vanishing
$C^{ab}_{(\mu)(\nu)}$. Similarly for the remaining non-vanishing coefficients
for the $F_{(0)(+)}$ component

$\displaystyle D^{t}_{(0)(+)}(\ell^{\prime}m^{\prime}|\ell m)$
$\displaystyle=C_{(+)}(\ell^{\prime}m^{\prime}|\ell m)/\sqrt{f}\text{,}$ (86)
$\displaystyle E^{t}_{(0)(+)}(\ell^{\prime}m^{\prime}|\ell m)$
$\displaystyle=E_{(+)}(\ell^{\prime}m^{\prime}|\ell m)/\sqrt{f}\text{.}$ (87)

The coupling coefficients for $F_{(+)(-)}$ contain both even and odd modes.
The first non-vanishing one is $D^{r}_{(+)(-)}(\ell^{\prime}m^{\prime}|\ell
m)$, which is given by

$\displaystyle D^{r}_{(+)(-)}(\ell^{\prime}m^{\prime}|\ell m)$
$\displaystyle=\sqrt{f}\gamma^{\ell,-m+1}C_{(-)}(\ell^{\prime}m^{\prime}|\ell+1,m-1)$
$\displaystyle\quad-\sqrt{f}\gamma^{\ell-1,m}C_{(-)}(\ell^{\prime}m^{\prime}|\ell-1,m-1)$
$\displaystyle\quad+\sqrt{f}\gamma^{\ell,m+1}C_{(+)}(\ell^{\prime}m^{\prime}|\ell+1,m+1)$
$\displaystyle-\sqrt{f}\gamma^{\ell-1,-m}C_{(+)}(\ell^{\prime}m^{\prime}|\ell-1,m+1)\text{,}$
(88)

while the coefficients coupling to odd modes are

$\displaystyle E^{r}_{(+)(-)}(\ell^{\prime}m^{\prime}|\ell m)$
$\displaystyle=\sqrt{f}\gamma^{\ell,-m+1}E_{(-)}(\ell^{\prime}m^{\prime}|\ell+1,m-1)$
$\displaystyle\quad-\sqrt{f}\gamma^{\ell-1,m}E_{(-)}(\ell^{\prime}m^{\prime}|\ell-1,m-1)$
$\displaystyle\quad+\sqrt{f}\gamma^{\ell,m+1}E_{(+)}(\ell^{\prime}m^{\prime}|\ell+1,m+1)$
$\displaystyle-\sqrt{f}\gamma^{\ell-1,-m}E_{(+)}(\ell^{\prime}m^{\prime}|\ell-1,m+1)\text{,}$
(89) $\displaystyle E_{(+)(-)}(\ell^{\prime}m^{\prime}|\ell m)$
$\displaystyle=-2i/r^{2}(\ell+1)(\ell+2)\epsilon^{\ell
m}\delta_{\ell^{\prime}\ell+1}\delta_{m^{\prime}m}$
$\displaystyle\quad-2i/r^{2}(\ell-1)\ell\epsilon^{\ell-1,m}\delta_{\ell^{\prime}\ell-1}\delta_{m^{\prime}m}\text{.}$
(90)

##  Appendix B Regularization parameters

I mimic the treatment in paper I and start from the covariant expression for
the singular vector potential tensor in Eq. (464) of Poisson (2004)

$\nabla_{\beta}A^{S}_{\alpha}(x)=-\frac{q}{2r^{2}}U_{\alpha\beta^{\prime}}u^{\beta^{\prime}}\nabla_{\beta}r-\frac{q}{2r_{\text{adv}}^{2}}U_{\alpha\beta^{\prime\prime}}u^{\beta^{\prime\prime}}\nabla_{\beta}r_{\text{adv}}\\\
\mbox{}+\frac{q}{2r}U_{\alpha\beta^{\prime};\beta}u^{\beta^{\prime}}+\frac{q}{2r}U_{\alpha\beta^{\prime};\gamma^{\prime}}u^{\beta^{\prime}}u^{\gamma^{\prime}}\nabla_{\beta}u+\frac{q}{2r_{\text{adv}}}U_{\alpha\beta^{\prime\prime};\beta}u^{\beta^{\prime\prime}}\\\
\mbox{}+\frac{q}{2r_{\text{adv}}}U_{\alpha\beta^{\prime\prime};\gamma^{\prime\prime}}u^{\beta^{\prime\prime}}u^{\gamma^{\prime\prime}}\nabla_{\beta}v+\frac{1}{2}qV_{\alpha\beta^{\prime}}u^{\beta^{\prime}}\nabla_{\beta}u\\\
\mbox{}-\frac{1}{2}qV_{\alpha\beta^{\prime\prime}}u^{\beta^{\prime\prime}}\nabla_{\beta}v-\frac{1}{2}q\int_{u}^{v}\nabla_{\beta}V_{\alpha\mu}\bm{(}x,z(\tau)\bm{)}u^{\beta}(\tau)\,\mathrm{d}\tau\text{.}$
(91)

I have introduced a large number of symbols. $x$ is the point where the field
is evaluated, $x^{\prime}$ and $x^{\prime\prime}$ are the retarded and
advanced points of $x$ on the world line $z(\tau)$. They are connected to $x$
with unique future-directed and past-directed null geodesics, respectively.
$u(x)$ and $v(x)$ are the retarded and advanced time functions such that
$x^{\prime}=z(\tau=u)$, $x^{\prime\prime}=z(\tau=v)$. $u^{\alpha^{\prime}}$
and $u^{\alpha^{\prime\prime}}$ are the four velocity at $x^{\prime}$ and
$x^{\prime\prime}$ respectively. Further I define Synge’s world function
$\sigma(x,\bar{x})$ which is numerically equal to half the squared geodesic
distance between two points $x$ and $\bar{x}$. Using its gradient
$\sigma_{\alpha}=\nabla_{\alpha}\sigma(x,\bar{x})$, I define
$r=u^{\alpha^{\prime}}\sigma_{\alpha^{\prime}}(x,x^{\prime})$ and
$r_{\text{adv}}=-u^{\alpha^{\prime\prime}}\sigma_{\alpha^{\prime\prime}}(x,x^{\prime\prime})$,
the affine parameter distances of $x$ away from the world line along its
future/past light cone. The potentials $U$ and $V$ appearing in Eq. (91) are
the direct and tail parts of the retarded Green function
$G_{\alpha\bar{\beta}}(x,\bar{x})$ associated with the wave operator.

From the definition of $r$, $r_{\text{adv}}$, $u$, and $v$ it follows that
(see Section 3.3.3 of Poisson (2004))

$\displaystyle\nabla_{\alpha}u=-\sigma_{\alpha}(x,x^{\prime})/r\text{,}$ (92)
$\displaystyle\nabla_{\alpha}v=\sigma_{\alpha}(x,x^{\prime\prime})/r_{\text{adv}}\text{,}$
(93)
$\displaystyle\nabla_{\alpha}r=-\sigma_{\alpha^{\prime}\beta^{\prime}}u^{\alpha^{\prime}}u^{\beta^{\prime}}\nabla_{\alpha}u+\sigma_{\alpha^{\prime}\alpha}u^{\alpha^{\prime}}\text{,}$
(94)
$\displaystyle\nabla_{\alpha}r_{\text{adv}}=-\sigma_{\alpha^{\prime\prime}\beta^{\prime\prime}}u^{\alpha^{\prime\prime}}u^{\beta^{\prime\prime}}\nabla_{\alpha}u-\sigma_{\alpha^{\prime\prime}\alpha}u^{\alpha^{\prime\prime}}\text{,}$
(95)

which are valid for geodesic motion.

The potentials $U_{\alpha\beta^{\prime}}$, $U_{\alpha\beta^{\prime\prime}}$
are determined by Eq. (322) of Poisson (2004)

$\displaystyle{U_{\alpha}}^{\beta^{\prime}}=g^{\beta^{\prime}}_{\phantom{\beta^{\prime}}\alpha}\Delta^{1/2}(x,x^{\prime})\text{,}$
(96)
$\displaystyle{U_{\alpha}}^{\beta^{\prime\prime}}=g^{\beta^{\prime\prime}}_{\phantom{\beta^{\prime\prime}}\alpha}\Delta^{1/2}(x,x^{\prime\prime})\text{,}$
(97)

where $g^{\bar{\mu}}_{\phantom{\bar{\mu}}\nu}$ is the parallel propagator from
$x^{\nu}$ to $\bar{x}^{\mu}$ and
$\Delta\equiv\det\bm{(}-g^{\alpha^{\prime}}_{\phantom{\alpha^{\prime}}\alpha}\sigma^{\alpha}_{;\beta^{\prime}}\bm{)}$
is the van Vleck determinant. Of the potentials $V_{\alpha\beta^{\prime}}$ and
$V_{\alpha\beta^{\prime\prime}}$ appearing in Eq. (91) I only need to know the
scaling behaviour following from Eq. (320) of Poisson (2004):

$\displaystyle V_{\alpha\beta^{\prime}}=O(\varepsilon^{2})\text{,}$ (98)
$\displaystyle V_{\alpha\beta^{\prime\prime}}=O(\varepsilon^{2})\text{,}$ (99)
$\displaystyle\nabla_{\beta}V_{\alpha\mu}=O(\varepsilon)\text{.}$ (100)

These expressions are valid in vacuum spacetimes where the Ricci tensor
vanishes.

Again mirroring the calculation in paper I, I introduce the arbitrary point
$\bar{x}\equiv z(\bar{\tau})$ on the world line and expand the quantities in
Eq. (91) in terms of a Taylor expansion around $\bar{x}$. I introduce the
convenient quantities

$\displaystyle\bar{r}\equiv\sigma_{\bar{\alpha}}(x,\bar{x})u^{\bar{\alpha}}\text{,}$
(101) $\displaystyle
s^{2}\equiv(g^{\bar{\alpha}\bar{\beta}}+u^{\bar{\alpha}}u^{\bar{\beta}})\sigma_{\bar{\alpha}}(x,\bar{x})\sigma_{\bar{\beta}}(x,\bar{x})\text{,}$
(102)

together with the time differences

$\Delta_{+}\equiv v-\bar{\tau}\text{,}\qquad\Delta_{-}\equiv u-\bar{\tau}$
(103)

from the advanced (retarded) point to the reference point $\bar{x}$.

I also use the expansion of the derivatives of the parallel propagator around
the point $\bar{x}$

$\displaystyle
g^{\bar{\alpha}}_{\phantom{\bar{\alpha}}\beta;\bar{\gamma}}=-g^{\bar{\beta}}_{\phantom{\bar{\beta}}\beta}\left(\frac{1}{2}R^{\bar{\alpha}}_{\phantom{\bar{\alpha}}\bar{\beta}\bar{\gamma}\bar{\delta}}\sigma^{\bar{\delta}}-\frac{1}{6}R^{\bar{\alpha}}_{\phantom{\bar{\alpha}}\bar{\beta}\bar{\gamma}\bar{\delta};\bar{\varepsilon}}\sigma^{\bar{\delta}}\sigma^{\bar{\varepsilon}}\right)+O(\varepsilon^{3})\text{,}$
(104)
$\displaystyle\begin{split}g^{\bar{\alpha}}_{\phantom{\bar{\alpha}}\beta;\gamma}&=-g^{\bar{\beta}}_{\phantom{\bar{\beta}}\beta}g^{\bar{\gamma}}_{\phantom{\bar{\gamma}}\gamma}\left(\frac{1}{2}R^{\bar{\alpha}}_{\phantom{\bar{\alpha}}\bar{\beta}\bar{\gamma}\bar{\delta}}\sigma^{\bar{\delta}}-\frac{1}{3}R^{\bar{\alpha}}_{\phantom{\bar{\alpha}}\bar{\beta}\bar{\gamma}\bar{\delta};\bar{\varepsilon}}\sigma^{\bar{\delta}}\sigma^{\bar{\varepsilon}}\right)\\\
&\quad+O(\varepsilon^{3})\text{,}\end{split}$ (105)

as well as an expansion for the second derivative of Synge’s world function

$\sigma_{\bar{\alpha}\bar{\beta}}=g_{\bar{\alpha}\bar{\beta}}-\frac{1}{3}R_{\bar{\alpha}\bar{\gamma}\bar{\beta}\bar{\delta}}\sigma^{\bar{\gamma}}\sigma^{\bar{\delta}}\\\
+\frac{1}{12}R_{\bar{\alpha}\bar{\gamma}\bar{\beta}\bar{\delta};\bar{\varepsilon}}\sigma^{\bar{\gamma}}\sigma^{\bar{\delta}}\sigma^{\bar{\varepsilon}}+O(\varepsilon^{4})\text{,}$
(106)

and the van Vleck determinant

$\Delta^{1/2}=1+O(\epsilon^{4})\text{,}$ (107)

which I calculate using the methods described in Sec. (2.4.2) of Poisson
(2004).

I make use of the fact that the bi-tensors

$\displaystyle U_{\alpha}(\tau)\equiv U_{\alpha\mu}u^{\mu}\text{,}$ (108)
$\displaystyle U_{\alpha\beta}(\tau)\equiv U_{\alpha\mu;\beta}u^{\mu}\text{,
and}$ (109) $\displaystyle\dot{U}_{\alpha}(\tau)\equiv
U_{\alpha\mu;\nu}u^{\mu}u^{\nu}$ (110)

appearing in Eq. (91) do not bear a free index on the world line, making them
scalars on the world line. With $w$ being either $u$ or $v$ and $\Delta\equiv
w-\bar{\tau}=\Delta_{\mp}$ I expand these as

$\displaystyle
U_{\alpha}(w)=U_{\alpha}+\dot{U}_{\alpha}\Delta+\frac{1}{2}\ddot{U}_{\alpha}\Delta^{2}+\frac{1}{6}U_{\alpha}^{(3)}\Delta^{3}+O(\varepsilon^{4})\text{,}$
(111) $\displaystyle
U_{\alpha\beta}(w)=U_{\alpha\beta}+\dot{U}_{\alpha\beta}\Delta+\frac{1}{2}\ddot{U}_{\alpha\beta}\Delta^{2}+O(\varepsilon^{3})\text{,
and}$ (112)
$\displaystyle\dot{U}_{\alpha}(w)=\dot{U}_{\alpha}+\ddot{U}_{\alpha}\Delta+\frac{1}{2}U_{\alpha}^{(3)}\Delta^{2}+O(\varepsilon^{3})\text{,}$
(113)

where it is understood that the coefficient functions are evaluated at
$\tau=\bar{\tau}$.

Repeatedly taking derivatives of Eq. (97) and contracting with $u^{\bar{\mu}}$
I find for the first set of coefficients

$\displaystyle
U_{\alpha}=g^{\bar{\alpha}}_{\phantom{\bar{\alpha}}\alpha}u_{\bar{\alpha}}+O(\varepsilon^{4})\text{,}$
(114)
$\displaystyle\dot{U}_{\alpha}=g^{\bar{\alpha}}_{\phantom{\bar{\alpha}}\alpha}\left(\frac{1}{2}R_{\bar{\alpha}uu\sigma}-\frac{1}{6}R_{\bar{\alpha}uu\sigma|\sigma}\right)+O(\varepsilon^{3})\text{,}$
(115)
$\displaystyle\ddot{U}_{\alpha}=\frac{1}{3}g^{\bar{\alpha}}_{\phantom{\bar{\alpha}}\alpha}R_{\bar{\alpha}uu\sigma|u}+O(\varepsilon^{2})\text{,}$
(116) $\displaystyle U^{(3)}_{\alpha}=0+O(\varepsilon)\text{,}$ (117)

where I have introduced the notation $R_{\bar{\alpha}uu\sigma}\equiv
R_{\bar{\alpha}\bar{\beta}\bar{\gamma}\bar{\delta}}u^{\bar{\beta}}u^{\bar{\gamma}}\sigma^{\bar{\delta}}$
and $R_{\bar{\alpha}uu\sigma|\sigma}\equiv
R_{\bar{\alpha}\bar{\beta}\bar{\gamma}\bar{\delta};\bar{\varepsilon}}u^{\bar{\beta}}u^{\bar{\gamma}}\sigma^{\bar{\delta}}\sigma^{\bar{\varepsilon}}$;
I will use this notation and its natural extension to higher derivatives and
different combinations of $u^{\bar{\mu}}$ and $\sigma^{\bar{\mu}}$ frequently
below.

Similarly I find for the second set

$\displaystyle
U_{\alpha\beta}=g^{\bar{\alpha}}_{\phantom{\bar{\alpha}}\alpha}g^{\bar{\beta}}_{\phantom{\bar{\beta}}\beta}\left(\frac{1}{2}R_{\bar{\alpha}u\bar{\beta}\sigma}-\frac{1}{3}R_{\bar{\alpha}u\bar{\beta}\sigma|\sigma}\right)+O(\varepsilon^{3})\text{,}$
(118)
$\displaystyle\begin{split}\dot{U}_{\alpha\beta}&=g^{\bar{\alpha}}_{\phantom{\bar{\alpha}}\alpha}g^{\bar{\beta}}_{\phantom{\bar{\beta}}\beta}\left(\frac{1}{2}R_{\bar{\alpha}u\bar{\beta}u}+\frac{1}{6}R_{\bar{\alpha}u\bar{\beta}\sigma|u}-\frac{1}{3}R_{\bar{\alpha}u\bar{\beta}u|\sigma}\right)\\\
&\quad+O(\varepsilon^{2})\text{,}\end{split}$ (119)
$\displaystyle\ddot{U}_{\alpha\beta}=\frac{1}{3}g^{\bar{\alpha}}_{\phantom{\bar{\alpha}}\alpha}g^{\bar{\beta}}_{\phantom{\bar{\beta}}\beta}R_{\bar{\alpha}u\bar{\beta}u|u}+O(\varepsilon)\text{.}$
(120)

Note that the third set does not involve new coefficients, but only those
already calculated for $U_{\alpha}$.

Finally I copy expressions for $\Delta_{\pm}$, $r$, $r_{\text{adv}}$, $u$, $v$
and their gradients from paper I

$\displaystyle\Delta_{\pm}=(\bar{r}\pm s)\mp\frac{(\bar{r}\pm
s)^{2}}{6s}R_{u\sigma u\sigma}\mp\frac{(\bar{r}\pm
s)^{2}}{24s}\left[(\bar{r}\pm s)R_{u\sigma u\sigma|u}-R_{u\sigma
u\sigma|\sigma}\right]+O(\varepsilon^{5})\text{,}$ (121) $\displaystyle
r=s-\frac{\bar{r}^{2}-s^{2}}{6s}R_{u\sigma
u\sigma}-\frac{\bar{r}-s}{24s}\left[(\bar{r}-s)(\bar{r}+2s)R_{u\sigma
u\sigma|u}-(\bar{r}+s)R_{u\sigma
u\sigma|\sigma}\right]+O(\varepsilon^{5})\text{,}$ (122) $\displaystyle
r_{\text{adv}}=s-\frac{\bar{r}^{2}-s^{2}}{6s}R_{u\sigma
u\sigma}-\frac{\bar{r}+s}{24s}\left[(\bar{r}+s)(\bar{r}-2s)R_{u\sigma
u\sigma|u}-(\bar{r}-s)R_{u\sigma
u\sigma|\sigma}\right]+O(\varepsilon^{5})\text{,}$ (123)
$\displaystyle\begin{split}\nabla_{\alpha}u&=\frac{1}{s}g^{\bar{\alpha}}_{\phantom{\bar{\alpha}}\alpha}\biggl{\\{}\left[\sigma_{\bar{\alpha}}+(\bar{r}-s)u_{\bar{\alpha}}\right]\\\
&\quad+\left[\frac{1}{6}(\bar{r}-s)R_{\bar{\alpha}\sigma
u\sigma}-\frac{1}{3}(\bar{r}-s)^{2}R_{\bar{\alpha}u\sigma
u}+\frac{\bar{r}^{2}-s^{2}}{6s^{2}}R_{u\sigma
u\sigma}\sigma_{\bar{\alpha}}+\frac{(\bar{r}-s)^{2}(\bar{r}+2s)}{6s^{2}}R_{\sigma
u\sigma u|u}u_{\bar{\alpha}}\right]\\\
&\quad+\biggl{[}-\frac{1}{12}(\bar{r}-s)R_{\bar{\alpha}u\sigma
u|\sigma}+\frac{1}{8}(\bar{r}-s)^{2}R_{\bar{\alpha}u\sigma
u|\sigma}+\frac{1}{24}(\bar{r}-s)^{2}R_{\bar{\alpha}\sigma
u\sigma|u}-\frac{1}{12}(\bar{r}-s)^{3}R_{\bar{\alpha}u\sigma u|u}\\\
&\quad+\frac{\bar{r}-s}{24s^{2}}\left((\bar{r}-s)(\bar{r}+2s)R_{u\sigma
u\sigma|u}-(\bar{r}+s)R_{u\sigma
u\sigma|\sigma}\right)\sigma_{\bar{\alpha}}\\\
&\quad+\frac{(\bar{r}-s)^{2}}{24s^{2}}\left((\bar{r}-s)(\bar{r}+3s)R_{u\sigma
u\sigma|\sigma}-(\bar{r}+2s)R_{u\sigma
u\sigma|\sigma}\right)u_{\bar{\alpha}}\biggr{]}+O(\varepsilon^{5})\biggr{\\}}\text{,}\end{split}$
(124)
$\displaystyle\begin{split}\nabla_{\alpha}v&=-\frac{1}{s}g^{\bar{\alpha}}_{\phantom{\bar{\alpha}}\alpha}\biggl{\\{}\left[\sigma_{\bar{\alpha}}+(\bar{r}+s)u_{\bar{\alpha}}\right]\\\
&\quad+\left[\frac{1}{6}(\bar{r}+s)R_{\bar{\alpha}\sigma
u\sigma}-\frac{1}{3}(\bar{r}+s)^{2}R_{\bar{\alpha}u\sigma
u}+\frac{\bar{r}^{2}-s^{2}}{6s^{2}}R_{u\sigma
u\sigma}\sigma_{\bar{\alpha}}+\frac{(\bar{r}+s)^{2}(\bar{r}-2s)}{6s^{2}}R_{\sigma
u\sigma u|u}u_{\bar{\alpha}}\right]\\\
&\quad+\biggl{[}-\frac{1}{12}(\bar{r}+s)R_{\bar{\alpha}u\sigma
u|\sigma}+\frac{1}{8}(\bar{r}+s)^{2}R_{\bar{\alpha}u\sigma
u|\sigma}+\frac{1}{24}(\bar{r}+s)^{2}R_{\bar{\alpha}\sigma
u\sigma|u}-\frac{1}{12}(\bar{r}+s)^{3}R_{\bar{\alpha}u\sigma u|u}\\\
&\quad+\frac{\bar{r}+s}{24s^{2}}\left((\bar{r}+s)(\bar{r}-2s)R_{u\sigma
u\sigma|u}-(\bar{r}-s)R_{u\sigma
u\sigma|\sigma}\right)\sigma_{\bar{\alpha}}\\\
&\quad+\frac{(\bar{r}+s)^{2}}{24s^{2}}\left((\bar{r}+s)(\bar{r}-3s)R_{u\sigma
u\sigma|\sigma}-(\bar{r}-2s)R_{u\sigma
u\sigma|\sigma}\right)u_{\bar{\alpha}}\biggr{]}+O(\varepsilon^{5})\biggr{\\}}\text{,}\end{split}$
(125)
$\displaystyle\begin{split}\nabla_{\alpha}r&=-\frac{1}{s}g^{\bar{\alpha}}_{\phantom{\bar{\alpha}}\alpha}\left\\{\left[\sigma_{\bar{\alpha}}+\bar{r}u_{\bar{\alpha}}\right]+\left[\frac{1}{6}\bar{r}R_{\bar{\alpha}\sigma
u\sigma}-\frac{1}{3}(\bar{r}^{2}-s^{2})R_{\bar{\alpha}u\sigma
u}+\frac{\bar{r}^{2}+s^{2}}{6s^{2}}R_{u\sigma
u\sigma}\sigma_{\bar{\alpha}}+\frac{\bar{r}(\bar{r}^{2}-s^{2})}{6s^{2}}R_{u\sigma
u\sigma}u_{\bar{\alpha}}\right]\right.\\\
&\quad\left.+\left[-\frac{1}{12}\bar{r}R_{\bar{\alpha}\sigma
u\sigma|\sigma}+\frac{1}{8}(\bar{r}^{2}-s^{2})R_{\bar{\alpha}u\sigma
u|\sigma}+\frac{1}{24}(\bar{r}^{2}-s^{2})R_{\bar{\alpha}\sigma
u\sigma|u}\right.\right.\\\
&\quad\left.\left.-\frac{1}{12}(\bar{r}-s)^{2}(\bar{r}+2s)R_{\bar{\alpha}u\sigma
u|u}+\frac{1}{24s^{2}}\left((\bar{r}-s)(\bar{r}^{2}+\bar{r}s+4s^{2})R_{u\sigma
u\sigma|u}-(\bar{r}^{2}+s^{2})R_{u\sigma
u\sigma|\sigma}\right)\sigma_{\bar{\alpha}}\right.\right.\\\
&\quad\left.\left.+\frac{\bar{r}-s}{24s^{2}}\left((\bar{r}-s)(\bar{r}^{2}+2\bar{r}s+3s^{2})R_{u\sigma
u\sigma|u}-\bar{r}(\bar{r}+s)R_{u\sigma
u\sigma|\sigma}\right)u_{\bar{\alpha}}\right]+O(\epsilon^{5})\right\\}\text{,}\end{split}$
(126)
$\displaystyle\begin{split}\nabla_{\alpha}r_{\text{adv}}&=-\frac{1}{s}g^{\bar{\alpha}}_{\phantom{\bar{\alpha}}\alpha}\left\\{\left[\sigma_{\bar{\alpha}}+\bar{r}u_{\bar{\alpha}}\right]+\left[\frac{1}{6}\bar{r}R_{\bar{\alpha}\sigma
u\sigma}-\frac{1}{3}(\bar{r}^{2}-s^{2})R_{\bar{\alpha}u\sigma
u}+\frac{\bar{r}^{2}+s^{2}}{6s^{2}}R_{u\sigma
u\sigma}\sigma_{\bar{\alpha}}+\frac{\bar{r}(\bar{r}^{2}-s^{2})}{6s^{2}}R_{u\sigma
u\sigma}u_{\bar{\alpha}}\right]\right.\\\
&\quad\left.+\left[-\frac{1}{12}\bar{r}R_{\bar{\alpha}\sigma
u\sigma|\sigma}+\frac{1}{8}(\bar{r}^{2}-s^{2})R_{\bar{\alpha}u\sigma
u|\sigma}+\frac{1}{24}(\bar{r}^{2}-s^{2})R_{\bar{\alpha}\sigma
u\sigma|u}\right.\right.\\\
&\quad\left.\left.-\frac{1}{12}(\bar{r}+s)^{2}(\bar{r}-2s)R_{\bar{\alpha}u\sigma
u|u}+\frac{1}{24s^{2}}\left((\bar{r}+s)(\bar{r}^{2}-\bar{r}s+4s^{2})R_{u\sigma
u\sigma|u}-(\bar{r}^{2}+s^{2})R_{u\sigma
u\sigma|\sigma}\right)\sigma_{\bar{\alpha}}\right.\right.\\\
&\quad\left.\left.+\frac{\bar{r}+s}{24s^{2}}\left((\bar{r}+s)(\bar{r}^{2}-2\bar{r}s+3s^{2})R_{u\sigma
u\sigma|u}-\bar{r}(\bar{r}-s)R_{u\sigma
u\sigma|\sigma}\right)u_{\bar{\alpha}}\right]+O(\epsilon^{5})\right\\}\text{.}\end{split}$
(127)

After substituting Eqs. (92) – (127) into Eq. (91) (all of them) and sorting
out the orders I find the final expression for the covariant expansion of
$A^{S}_{\alpha;\beta}$

$A^{S}_{\alpha;\beta}=qg^{\bar{\alpha}}_{\phantom{\bar{\alpha}}\alpha}g^{\bar{\beta}}_{\phantom{\bar{\beta}}\beta}\left\\{\left[\frac{1}{s^{3}}u_{\bar{\alpha}}\sigma_{\bar{\beta}}+\frac{\bar{r}}{s^{3}}u_{\bar{\alpha}}u_{\bar{\beta}}\right]+\left[\frac{\bar{r}}{6s^{3}}u_{\bar{\alpha}}R_{\bar{\beta}\sigma
u\sigma}+\left(\frac{\bar{r}}{2s^{3}}\sigma_{\bar{\beta}}+\frac{\bar{r}^{2}-s^{2}}{2s^{3}}u_{\bar{\beta}}\right)R_{\bar{\alpha}uu\sigma}\right.\right.\\\
\mbox{}\left.\left.+\frac{\bar{r}^{2}-s^{2}}{3s^{2}}u_{\bar{\alpha}}R_{\bar{\beta}uu\sigma}+\frac{1}{2s}R_{\bar{\alpha}u\bar{\beta}\sigma}+\frac{3\bar{r}^{2}-s^{2}}{6s^{5}}R_{u\sigma
u\sigma}u_{\bar{\alpha}}\sigma_{\bar{\beta}}+\frac{\bar{r}(\bar{r}^{2}-s^{2})}{2s^{5}}R_{u\sigma
u\sigma}u_{\bar{\alpha}}u_{\bar{\beta}}+\frac{\bar{r}}{2s}R_{\bar{\alpha}u\bar{\beta}u}\right]\right.\\\
\mbox{}\left.+\left[-\frac{\bar{r}}{12s^{3}}u_{\bar{\alpha}}R_{\bar{\beta}\sigma
u\sigma|\sigma}-\frac{\bar{r}^{2}-s^{2}}{24s^{3}}u_{\bar{\alpha}}R_{\bar{\beta}\sigma\sigma
u|u}-\left(\frac{\bar{r}}{6s^{3}}\sigma_{\bar{\beta}}+\frac{\bar{r}^{2}-s^{2}}{6s^{3}}u_{\bar{\beta}}\right)R_{\bar{\alpha}uu\sigma|\sigma}-\frac{\bar{r}^{2}-s^{2}}{8s^{3}}u_{\bar{\alpha}}R_{\bar{\beta}uu\sigma|\sigma}\right.\right.\\\
\mbox{}\left.\left.-\frac{1}{3s^{2}}R_{\bar{\alpha}u\bar{\beta}\sigma|\sigma}+\left(\frac{\bar{r}^{2}-s^{2}}{6s^{3}}\sigma_{\bar{\beta}}+\frac{\bar{r}(\bar{r}^{2}-3s^{2})}{6s^{3}}u_{\bar{\beta}}\right)R_{\bar{\alpha}uu\sigma|u}+\frac{\bar{r}(\bar{r}^{2}-3s^{2})}{12s^{3}}u_{\bar{\alpha}}R_{\bar{\beta}uu\sigma|u}-\frac{\bar{r}}{3s}R_{\bar{\alpha}u\bar{\beta}u|\sigma}\right.\right.\\\
\mbox{}\left.\left.+\frac{\bar{r}}{6s}R_{\bar{a}u\bar{b}\sigma|u}+\frac{\bar{r}^{2}+s^{2}}{6s^{3}}R_{\bar{a}u\bar{b}u|u}+\left(-\frac{3\bar{r}^{2}-s^{2}}{24s^{5}}R_{u\sigma
u\sigma|\sigma}+\frac{\bar{r}(\bar{r}^{2}-s^{2})}{8s^{5}}R_{u\sigma
u\sigma|u}\right)u_{\bar{\alpha}}\sigma_{\bar{\beta}}\right.\right.\\\
\mbox{}\left.\left.+\left(-\frac{\bar{r}(\bar{r}^{2}-s^{2})}{8s^{5}}R_{u\sigma
u\sigma|\sigma}+\frac{(\bar{r}^{2}-s^{2})^{2}}{8s^{5}}R_{u\sigma
u\sigma|u}\right)u_{\bar{\alpha}}u_{\bar{\beta}}\right]\right\\}+O(\varepsilon^{2})\text{,}$
(128)

where terms in square brackets are of the same power in $\varepsilon$.

I copy the results for the coordinate expansion of
$\sigma_{\bar{\alpha}}(x,\bar{x})$ and
$g^{\bar{\alpha}}_{\phantom{\bar{\alpha}}\beta}(x,\bar{x})$ from Eqs. (3.16) –
(3.19) and Eqs. (3.30) – (3.33) of paper I. I use

$\displaystyle\begin{split}-\sigma_{\bar{\alpha}}(x,\bar{x})&=g_{\alpha\beta}w^{\beta}+A_{\alpha\beta\gamma}w^{\beta}w^{\gamma}+A_{\alpha\beta\gamma\delta}w^{\beta}w^{\gamma}w^{\delta}\\\
&\mbox{}+A_{\alpha\beta\gamma\delta\varepsilon}w^{\beta}w^{\gamma}w^{\delta}w^{\varepsilon}+O(\varepsilon^{5})\text{,}\end{split}$
(129) $\displaystyle
A^{\alpha}_{\phantom{\alpha}\beta\gamma}\equiv\frac{1}{2}\Gamma^{\alpha}_{\phantom{\alpha}\beta\gamma}\text{,}$
(130) $\displaystyle
A^{\alpha}_{\phantom{\alpha}\beta\gamma\delta}\equiv\frac{1}{6}\left(\Gamma^{\alpha}_{\phantom{\alpha}\beta\gamma,\delta}+\Gamma^{\alpha}_{\phantom{\alpha}\beta\mu}\Gamma^{\mu}_{\phantom{\mu}\gamma\delta}\right)\text{,}$
(131)
$\displaystyle\begin{split}A^{\alpha}_{\phantom{\alpha}\beta\gamma\delta\epsilon}&\equiv\frac{1}{24}\left(\Gamma^{\alpha}_{\phantom{\alpha}\beta\gamma,\delta\epsilon}+\Gamma^{\alpha}_{\phantom{\alpha}\beta\gamma,\mu}\Gamma^{\mu}_{\phantom{\mu}\delta\epsilon}\right.\\\
&\mbox{}\left.+2\Gamma^{\alpha}_{\phantom{\alpha}\beta\mu}\Gamma^{\mu}_{\phantom{\mu}\gamma\delta,\epsilon}+\Gamma^{\alpha}_{\phantom{\alpha}\mu\nu}\Gamma^{\mu}_{\phantom{\mu}\beta\gamma}\Gamma^{\nu}_{\phantom{\nu}\delta\epsilon}\right)\text{,}\end{split}$
(132)

as well as

$\displaystyle\begin{split}g^{\bar{\mu}}_{\phantom{\bar{\mu}}\alpha}(x,\bar{x})&={\delta^{\mu}}_{\alpha}+B^{\mu}_{\phantom{\mu}\alpha\beta}w^{\beta}+B^{\mu}_{\phantom{\mu}\alpha\beta\gamma}w^{\beta}w^{\gamma}\\\
&\mbox{}+B^{\mu}_{\phantom{\mu}\alpha\beta\gamma\delta}w^{\beta}w^{\gamma}w^{\delta}+O(\varepsilon^{4})\text{,}\text{,}\end{split}$
(133) $\displaystyle
B^{\mu}_{\phantom{\mu}\alpha\beta}\equiv\Gamma^{\mu}_{\phantom{\mu}\alpha\beta}\text{,}$
(134) $\displaystyle
B^{\mu}_{\phantom{\mu}\alpha\beta\gamma}\equiv\frac{1}{2}\left(\Gamma^{\mu}_{\phantom{\mu}\alpha\beta,\gamma}+\Gamma^{\mu}_{\phantom{\mu}\beta\nu}\Gamma^{\nu}_{\phantom{\nu}\alpha\gamma}\right)\text{,}$
(135)
$\displaystyle\begin{split}B^{\mu}_{\phantom{\mu}\alpha\beta\gamma\delta}&\equiv\frac{1}{12}\left(2\Gamma^{\mu}_{\phantom{\mu}\alpha\beta,\gamma\delta}+2\Gamma^{\nu}_{\phantom{\nu}\alpha\beta}\Gamma^{\mu}_{\phantom{\mu}\nu\gamma,\delta}\right.\\\
&\mbox{}\left.-\Gamma^{\nu}_{\phantom{\nu}\beta\gamma}\Gamma^{\mu}_{\phantom{\mu}\alpha\nu,\delta}+4\Gamma^{\mu}_{\phantom{\mu}\beta\nu}\Gamma^{\nu}_{\phantom{\nu}\alpha\gamma,\delta}\right.\\\
&\mbox{}\left.+\Gamma^{\nu}_{\phantom{\nu}\beta\gamma}\Gamma^{\mu}_{\phantom{\mu}\alpha\delta,\nu}-\Gamma^{\mu}_{\phantom{\mu}\beta\nu}\Gamma^{\nu}_{\phantom{\nu}\alpha\lambda}\Gamma^{\lambda}_{\phantom{\lambda}\gamma\delta}\right.\\\
&\mbox{}\left.+\Gamma^{\mu}_{\phantom{\mu}\nu\lambda}\Gamma^{\nu}_{\phantom{\nu}\alpha\beta}\Gamma^{\lambda}_{\phantom{\lambda}\gamma\delta}+2\Gamma^{\mu}_{\phantom{\mu}\beta\nu}\Gamma^{\nu}_{\phantom{\nu}\gamma\lambda}\Gamma^{\lambda}_{\phantom{\lambda}\alpha\delta}\right)\text{,}\end{split}$
(136)

where $w^{\alpha}\equiv x^{\alpha}-x^{\bar{\alpha}}$ is the coordinate
distance between $x$ and $\bar{x}$. Together with Eq. (128) these equations
form an expansion in $w^{\alpha}$ of the singular part of the gradient of the
vector potential around a point $x$ near the world line of the particle. I
finally calculate the tetrad components of the singular Faraday tensor as

$F^{S}_{(\mu)(\nu)}=(A^{S}_{\beta,\alpha}-A^{S}_{\alpha,\beta})e^{\alpha}_{\
(\mu)}e^{\beta}_{\ (\nu)}\text{.}$ (137)

From this point on I proceed exactly as described in Section V of paper I
using Maple and GRTensorII to perform the calculations. I find, after an
extremely tedious calculation,

$\displaystyle
A_{(0)(+)}=\operatorname{sign}(\Delta)\Bigl{[}\frac{i\dot{r}_{0}J}{r_{0}\mathfrak{f}a^{2}}-\frac{1}{r_{0}^{2}}\Bigr{]}e^{i\varphi_{0}}\text{,}$
(138) $\displaystyle
B_{(0)(+)}=\biggl{\\{}\Bigl{[}-\frac{iE(J^{2}-r_{0}^{2})\dot{r}_{0}}{a^{3}\pi\mathfrak{f}J}+\frac{E(2-\mathfrak{f})}{\pi\mathfrak{f}ar_{0}}\Bigr{]}\mathcal{E}-\frac{ir_{0}^{2}E\dot{r}_{0}}{a^{3}J\mathfrak{f}\pi}\mathcal{K}\biggr{\\}}e^{i\varphi_{0}}\text{,}$
(139)
$\displaystyle\begin{split}D&{}_{(0)(+)}=\Bigg{\\{}\biggl{[}\frac{iEr_{0}^{2}(-14r_{0}^{2}J^{2}+J^{4}+r_{0}^{4})\dot{r}_{0}^{3}}{8\pi
J\mathfrak{f}a^{7}}-\frac{(-r_{0}\mathfrak{f}J^{2}+2r_{0}J^{2}+7r_{0}^{3}\mathfrak{f}-14r_{0}^{3})E\dot{r}_{0}^{2}}{8a^{5}\mathfrak{f}\pi}\\\
&\mbox{}+i\frac{\bigl{(}8MJ^{8}-14MJ^{6}r_{0}^{2}-3r_{0}^{5}J^{4}-80MJ^{4}r_{0}^{4}+4J^{4}r_{0}^{5}\mathfrak{f}-7r_{0}^{7}J^{2}-68Mr_{0}^{6}J^{2}+4r_{0}^{9}-26Mr_{0}^{8}-4r_{0}^{9}\mathfrak{f}\bigr{)}E\dot{r}_{0}}{8r_{0}^{5}a^{5}\mathfrak{f}J\pi}\\\
&\mbox{}-\frac{(8Mr_{0}\mathfrak{f}J^{6}-8r_{0}^{3}MJ^{4}+38J^{4}r_{0}^{3}\mathfrak{f}-2r_{0}^{6}J^{2}-16Mr_{0}^{5}J^{2}+3J^{2}r_{0}^{6}\mathfrak{f}+54J^{2}r_{0}^{5}\mathfrak{f}+20r_{0}^{7}\mathfrak{f}+5r_{0}^{8}\mathfrak{f}-6r_{0}^{8})E}{8r_{0}^{7}a^{3}\mathfrak{f}\pi}\biggr{]}\mathcal{E}\\\
&\mbox{}+\biggl{[}\frac{iEr_{0}^{4}(7J^{2}-r_{0}^{2})\dot{r}_{0}^{3}}{8\pi
J\mathfrak{f}a^{7}}-\frac{(2-\mathfrak{f})r_{0}^{3}E\dot{r}_{0}^{2}}{2a^{5}\mathfrak{f}\pi}\\\
&\mbox{}+\frac{\bigl{(}4Mr_{0}\mathfrak{f}J^{4}+20J^{2}r_{0}^{3}\mathfrak{f}+8Mr_{0}^{3}J^{2}+14r_{0}^{5}M\mathfrak{f}-2r_{0}^{6}+12Mr_{0}^{5}+r_{0}^{6}\mathfrak{f}\bigr{)}E}{8r_{0}^{5}a^{3}\mathfrak{f}\pi}\\\
&\mbox{}-\frac{i(2MJ^{6}-9Mr_{0}^{2}J^{4}-2J^{2}r_{0}^{5}\mathfrak{f}-20Mr_{0}^{4}J^{2}-2r_{0}^{7}\mathfrak{f}+2r_{0}^{7}-13Mr_{0}^{6})E\dot{r}_{0}}{4r_{0}^{3}a^{5}\mathfrak{f}J\pi}\biggr{]}\mathcal{K}\Bigg{\\}}e^{i\varphi_{0}}\text{,}\end{split}$
(140) $\displaystyle
A_{(+)(-)}=\operatorname{sign}(\Delta)\frac{2iEJ}{a^{2}r_{0}\mathfrak{f}}e^{i\varphi_{0}}\text{,}$
(141)
$\displaystyle\begin{split}B_{(+)(-)}&=-2i\Biggl{\\{}\biggl{[}-\frac{(r_{0}^{2}-J^{2})\dot{r}_{0}^{2}}{a^{3}\pi
J\mathfrak{f}}+\frac{-J^{2}r_{0}\mathfrak{f}+2r_{0}J^{2}+2r_{0}^{3}-2r_{0}^{3}\mathfrak{f}}{r_{0}^{3}aJ\pi}\biggr{]}\mathcal{E}+\biggl{[}\frac{r_{0}^{2}\dot{r}_{0}^{2}}{a^{3}\mathfrak{f}J\pi}-\frac{2(1-\mathfrak{f})}{aJ\pi}\biggr{]}\mathcal{K}\Biggr{\\}}e^{i\varphi_{0}}\text{,}\end{split}$
(142)
$\displaystyle\begin{split}D_{(+)(-)}&=-2i\Bigg{\\{}\biggl{[}-\frac{r_{0}^{2}(-14r_{0}^{2}J^{2}+J^{4}+r_{0}^{4})\dot{r}_{0}^{4}}{8\mathfrak{f}\pi
Ja^{7}}-\Bigl{(}4M\mathfrak{f}J^{8}-7Mr_{0}^{2}\mathfrak{f}J^{6}+2J^{4}r_{0}^{5}\mathfrak{f}+2J^{4}r_{0}^{4}M-J^{4}r_{0}^{5}\\\
&\mbox{}-43J^{4}r_{0}^{4}\mathfrak{f}M-7J^{2}r_{0}^{7}\mathfrak{f}-27J^{2}Mr_{0}^{6}\mathfrak{f}-11Mr_{0}^{8}\mathfrak{f}-r_{0}^{9}\mathfrak{f}+r_{0}^{9}-2r_{0}^{8}M\Bigr{)}r_{0}^{1/2}\dot{r}_{0}^{2}\Big{/}\Bigl{(}4r_{0}^{5}a^{5}J\mathfrak{f}\pi\Bigr{)}\\\
&\mbox{}-\Bigl{(}8M\mathfrak{f}J^{8}-8J^{6}Mr_{0}^{2}+30Mr_{0}^{2}\mathfrak{f}J^{6}-2J^{4}r_{0}^{5}+10J^{4}r_{0}^{4}\mathfrak{f}M-24J^{4}r_{0}^{4}M+3J^{4}r_{0}^{5}\mathfrak{f}\\\
&\mbox{}-28J^{2}Mr_{0}^{6}+J^{2}r_{0}^{7}\mathfrak{f}-28J^{2}Mr_{0}^{6}\mathfrak{f}-20Mr_{0}^{8}\mathfrak{f}-12r_{0}^{8}M\Bigr{)}\Big{/}\Bigl{(}8r_{0}^{7}a^{3}J\pi\Bigr{)}\biggr{]}\mathcal{E}\\\
&\mbox{}+\biggl{[}-\frac{r_{0}^{4}(7J^{2}-r_{0}^{2})\dot{r}_{0}^{4}}{8\mathfrak{f}\pi
Ja^{7}}\\\
&\mbox{}+\frac{4M\mathfrak{f}J^{6}-16J^{4}Mr_{0}^{2}+4r_{0}^{2}\mathfrak{f}MJ^{4}-18J^{2}r_{0}^{4}\mathfrak{f}M-28J^{2}r_{0}^{4}M-J^{2}\mathfrak{f}r_{0}^{5}-12r_{0}^{6}M-20r_{0}^{6}\mathfrak{f}M}{8r_{0}^{5}a^{3}J\pi}\\\
&\mbox{}+\Bigl{(}2M\mathfrak{f}J^{6}-9r_{0}^{2}\mathfrak{f}MJ^{4}+J^{2}r_{0}^{5}-2J^{2}r_{0}^{4}M-5J^{2}\mathfrak{f}r_{0}^{5}-14J^{2}r_{0}^{4}\mathfrak{f}M-2r_{0}^{6}M+r_{0}^{7}\\\
&\mbox{}-11r_{0}^{6}\mathfrak{f}M-\mathfrak{f}r_{0}^{7}\Bigr{)}\dot{r}_{0}^{2}\Big{/}\Bigl{(}4r_{0}^{5/2}a^{5}J\mathfrak{f}\pi\Bigr{)}\biggr{]}\mathcal{K}\Bigg{\\}}e^{i\varphi_{0}}\text{,}\end{split}$
(143)

where $\mathfrak{f}=\sqrt{\frac{r_{0}-2M}{r_{0}}}$, $a^{2}=r_{0}^{2}+J^{2}$.
Here, the rescaled elliptic integrals $\mathcal{E}$ and $\mathcal{K}$ are
defined by

$\mathcal{E}\equiv\frac{2}{\pi}\int_{0}^{\pi/2}(1-k\sin^{2}\psi)^{1/2}\,\mathrm{d}\psi=F\left(-{\frac{1}{2}},{\frac{1}{2}};1;k\right)\text{,}$
(144)

and

$\mathcal{K}\equiv\frac{2}{\pi}\int_{0}^{\pi/2}(1-k\sin^{2}\psi)^{-1/2}\,\mathrm{d}\psi=F\left({\frac{1}{2}},{\frac{1}{2}};1;k\right)\text{,}$
(145)

in which $F(a,b;c;x)$ are the hypergeometric functions and $k\equiv
J^{2}/(r_{0}^{2}+J^{2})$.

##  Appendix C Vector potential calculation

In this section I describe a variant of the numerical calculation discussed in
the main part of the paper that uses the vector potential instead of the
Faraday tensor. To this end I decompose the vector potential and the sources
in terms of vectorial spherical harmonics

$\displaystyle A_{a}(t,r^{*},\theta,\phi)$
$\displaystyle={\textstyle\frac{1}{r}}A^{\ell m}_{a}(t,r^{*})Y_{\ell
m}(\theta,\phi)\text{,}$ (146a) $\displaystyle j_{a}(t,r^{*},\theta,\phi)$
$\displaystyle=j^{\ell m}_{a}(t,r^{*})Y_{\ell m}(\theta,\phi)$ for
$a=t,r^{*}$, (146b) $\displaystyle A_{A}(t,r^{*},\theta,\phi)$
$\displaystyle=v_{\ell m}(t,r^{*})Z_{A}^{\ell m}(\theta,\phi)$
$\displaystyle\quad+\tilde{v}_{\ell m}(t,r^{*})X_{A}^{\ell
m}(\theta,\phi)\text{,}$ (146c) $\displaystyle j_{A}(t,r^{*},\theta,\phi)$
$\displaystyle=j^{\text{even}}_{\ell m}(t,r^{*})Z_{A}^{\ell m}(\theta,\phi)$
$\displaystyle\quad+j^{\text{odd}}_{\ell m}(t,r^{*})X_{A}^{\ell
m}(\theta,\phi)$ for $A=\theta,\phi$, (146d)

and substitute this into the Maxwell equations for the vector potential in the
Lorenz gauge $g^{\alpha\beta}A_{\alpha;\beta}=0$:

$g^{\mu\nu}A_{\alpha;\mu\nu}-{R^{\beta}}_{\alpha}A^{\beta}=-4\pi
j_{\alpha}\text{,}$ (147)

where $R_{\alpha\beta}$ is the spacetime’s Ricci tensor, which vanishes in
Schwarzschild spacetime. Substituting Eq. (146) into Eq. (147) I arrive at two
decoupled sets of equations for the even ($A^{\ell m}_{a}$, $v_{\ell m}$) and
odd ($\tilde{v}_{\ell m}$) modes

$\displaystyle-\frac{\partial^{2}A^{\ell
m}_{t}}{\partial{t}^{2}}+\frac{\partial^{2}A^{\ell
m}_{t}}{\partial{r^{*}}^{2}}+\frac{2M}{r^{2}}\left(\frac{\partial A^{\ell
m}_{r^{*}}}{\partial t}-\frac{\partial A^{\ell m}_{t}}{\partial
r^{*}}\right)-VA^{\ell m}_{t}=-4\pi rfj^{\ell m}_{t}\text{,}$ (148)
$\displaystyle\begin{split}-\frac{\partial^{2}A^{\ell
m}_{r^{*}}}{\partial{t}^{2}}+\frac{\partial^{2}A^{\ell
m}_{r^{*}}}{\partial{r^{*}}^{2}}+\frac{2M}{r^{2}}\left(\frac{\partial A^{\ell
m}_{t}}{\partial t}-\frac{\partial A^{\ell m}_{r^{*}}}{\partial
r^{*}}\right)-\left(V+2\frac{f^{2}}{r^{2}}\right)A^{\ell m}_{r^{*}}+fVv_{\ell
m}=-4\pi rfj^{\ell m}_{r^{*}}\text{,}\end{split}$ (149)
$\displaystyle-\frac{\partial^{2}v_{\ell
m}}{\partial{t}^{2}}+\frac{\partial^{2}v_{\ell
m}}{\partial{r^{*}}^{2}}-Vv_{\ell m}+2\frac{f}{r^{2}}A^{\ell m}_{r^{*}}=-4\pi
fj^{\text{even}}_{\ell m}\text{,}$ (150)
$\displaystyle-\frac{\partial^{2}\tilde{v}_{\ell
m}}{\partial{t}^{2}}+\frac{\partial^{2}\tilde{v}_{\ell
m}}{\partial{r^{*}}^{2}}-V\tilde{v}_{\ell m}=-4\pi fj^{\text{odd}}_{\ell
m}\text{,}$ (151)

where $V$ and $j^{\ell m}_{\alpha}$ is defined as in Eqs. (21b) and (10) in
the main text.

###  C.1 Numerical method

I discretize the set of reduced equations Eqs. (148) – (151) using Lousto’s
method as described in section II of the main text. Since the source terms on
the right hand side are less singular for the vector potential than they are
for the Faraday tensor, I do not have to distinguish between sourced and
vacuum cells in the integral over the potential terms.

Terms containing first derivatives $\frac{\partial\psi}{\partial t}$,
$\frac{\partial\psi}{\partial r^{*}}$, where now and in the remainder of the
appendix $\psi$ stands for any of $A^{\ell m}_{t}$, $A^{\ell m}_{r^{*}}$,
$v^{\ell m}$ or $\tilde{v}^{\ell m}$, were not treated in Lousto and Price
(1997), but, for generic vacuum cells, can be handled in a straightforward
manner

$\displaystyle\iint_{\text{\hbox
to0.0pt{cell\hss}}}\,\mathrm{d}u\,\mathrm{d}v\,V(r)\frac{\partial\psi}{\partial
t}=2h(\psi_{3}-\psi_{2})V_{0}+O(h^{4})\text{,}$ (152)
$\displaystyle\iint_{\text{\hbox
to0.0pt{cell\hss}}}\,\mathrm{d}u\,\mathrm{d}v\,V(r)\frac{\partial\psi}{\partial
r^{*}}=2h(\psi_{4}-\psi_{1})V_{0}+O(h^{4})\text{.}$ (153)

This fails for cells traversed by the particle, since the field is only
continuous across the world line but not differentiable. For these cells I
take recourse to Lousto’s original algorithm, which has to deal with a similar
issue, and use

$\displaystyle\iint_{\text{\hbox
to0.0pt{cell\hss}}}\,\mathrm{d}u\,\mathrm{d}v\,V(r)\frac{\partial\psi}{\partial
t}=V_{0}\,\sum_{i}A_{i}\partial_{t}\psi_{i}+O(h^{3})\text{,}$ (154)
$\displaystyle\iint_{\text{\hbox
to0.0pt{cell\hss}}}\,\mathrm{d}u\,\mathrm{d}v\,V(r)\frac{\partial\psi}{\partial
r^{*}}=V_{0}\,\sum_{i}A_{i}\partial_{r^{*}}\psi_{i}+O(h^{3})\text{,}$ (155)

where $A_{1}$,…,$A_{4}$ are the subareas indicated in Fig. 3 and
$\partial_{t}\psi_{1}$, …, $\partial_{t}\psi_{4}$, $\partial_{r^{*}}\psi_{1}$,
…, $\partial_{r^{*}}\psi_{4}$ are zeroth order accurate approximations to the
derivatives in the subareas. I calculate these using grid points outside of
the cell on the same side of the world line as the corresponding subarea, e.g.

$\partial_{r^{*}}\psi_{1}=\frac{\psi(t,r^{*}-h)-\psi(t,r^{*}-3h)}{2h}+O(h)\text{.}$
(156)

###  C.2 Gauge condition

In contrast to the scalar field, the electromagnetic vector potential has to
satisfy a gauge condition

$Z\equiv g^{\alpha\beta}A_{\alpha;\beta}=0\text{.}$ (157)

Analytically the gauge condition is preserved by the evolution equations, so
that it is sufficient to impose it on the initial data. Numerically, however,
small violations of the gauge condition due to the numerical approximation can
be amplified exponentially and come to dominate the numerical data. To handle
this situation I introduce a gauge damping scheme as described in Gundlach et
al. (2005); Barack and Lousto (2005). That is I add a term of the form

$\frac{4M}{r^{2}}Z=\frac{4M}{r^{2}}\left(-\frac{1}{r-2M}\frac{\partial
A_{t}}{\partial t}+\frac{1}{r-2M}\frac{\partial A_{r^{*}}}{\partial
r^{*}}\right.\\\
\mbox{}\left.+\frac{1}{r^{2}}A_{r^{*}}-\frac{\ell(\ell+1)}{r^{2}}v\right)$
(158)

to the $t$ components of the evolution equations Eqs. (147), which dampens out
violations of the gauge condition. This choice proved to be numerically stable
for the radiative ($\ell>0$) modes but unstable for the monopole ($\ell=0$)
mode.

###  C.3 Monopole mode

The monopole moment of an electromagnetic field is non-radiative. This makes
its behaviour sufficiently different from that of the radiative ($\ell>0$)
modes that the approach outlined earlier fails for $\ell=0$. In this case Eq.
(147) reduces to a set of coupled equations for $A_{a}^{0,0}$ only. Rather
than solving the system of equations directly for $A_{t}^{0,0}$ and
$A_{r^{*}}^{0,0}$ I use the analytical result for the $F_{tr}$ component of
the Faraday tensor derived in section II.3 in the main part of the paper. This
proves to be sufficient to reconstruct the combination
$A^{0,0}_{r,t}-A^{0,0}_{t,r}$ appearing in Eq. (77).

###  C.4 Initial values and boundary conditions

I handle the problem of initial data and boundary conditions the same way as
in the main text, that is I arbitrarily choose the fields to vanish on the
characteristic slices $u=u_{0}$ and $v=v_{0}$

$A_{\alpha}(u=u_{0})=A_{\alpha}(v=v_{0})=0\text{,}$ (159)

thereby adding a certain amount of spurious waves to the solution which show
up as an initial burst. Gauge violations in this initial data are damped out
along with those arising during the evolution.

I implement ingoing wave boundary conditions near the event horizon and choose
a numerical domain that covers the full domain of dependence of the initial
data near the outer boundary.

###  C.5 Extraction of the field data at the particle

In order to extract the value of the fields and their first derivatives at the
position of the particle, I use a variant of the extraction scheme described
in paper II. I introduce a piecewise polynomial

$p(x)=\begin{cases}c_{0}+c_{1}x+\frac{c_{3}}{2}x^{2}&\text{if $x<0$}\\\
c^{\prime}_{0}+c^{\prime}_{1}x+\frac{c^{\prime}_{3}}{2}x^{2}&\text{if
$x>0$}\end{cases}$ (160)

in $x\equiv r^{*}-r_{0}^{*}$ on the current slice. Its coefficients to the
left and right of the world line are linked by jump conditions
$c_{n}=c^{\prime}_{n}+\left[\partial^{n}_{r^{*}}\psi\right]$ listed in
Appendix D.2. Fitting this polynomial to the three grid points closest to the
particle, I extract approximations for $\psi(t_{0},r^{*}_{0})$ and
$\frac{\partial\psi(t_{0},r^{*}_{0})}{\partial r^{*}}$ which are just the
coefficients $c_{0}$, $c_{1}$ respectively. Once I have obtained these, I
proceed as in section V of the main part of the paper following Sago  to
obtain values for $\frac{\partial\psi(t_{0},r^{*}_{0})}{\partial t}$.

###  C.6 Results

Using the vector potential code described above I can reproduce the results
obtained from the Faraday tensor method discussed in the main paper. The
differences are small, typically of the order of $10^{-3}\%$ of the field
values as shown in Fig. 20.

Figure 20: Differences between $F_{tr}^{\ell m}$ calculated using the vector
potential and calculated using the Faraday tensor method for $\ell=2$, $m=2$
mode of field for the zoom-whirl orbit shown in Fig. 8. Displayed are the
difference and the actual field. The stepsizes were $h=1.041\bar{6}\times
10^{-2}\,M$ and $h=1/512\,M$ for the vector potential calculation and the
Faraday tensor calculation respectively.

I expect the Faraday tensor code to yield more accurate results since the
costly numerical differentiation that is necessary in the vector potential
calculation is absent. Nevertheless I can reproduce e.g. the correct decay
behaviour of the multipole coefficients for a zoom-whirl orbit as shown in
Fig. 21.

Figure 21: Multipole coefficients of
$\frac{M^{2}}{q}\operatorname{Re}F^{R}_{(0)}$ for a particle on a zoom-whirl
orbit ($p=7.8001$, $e=0.9$), calculated using a stepsize of $h=0.125M$ for the
$\ell=1$ modes and increasing the resolution linearly with $\ell$ for
$\ell>1$. The coefficients are extracted at $t=1100\,M$ when the particle is
deep within the zoom phase. Red triangles are used for the unregularized
multipole coefficients $F_{(0),\ell}$, squares, diamonds and disks are used
for the partly regularized coefficients after the removal of the $A_{(0)}$,
$B_{(0)}$ and $D_{(0)}$ terms respectively.

##  Appendix D Jump conditions

###  D.1 Faraday tensor calculation

Since the source term in Eqs. (21b) – (21f) contains a term proportional to
$\delta^{\prime}(r^{*}-r^{*}_{0})$, the field is discontinuous across the
world line of the particle. I use

$\left[\partial_{t}^{n}\partial_{r^{*}}^{m}\psi\right]=\lim_{\varepsilon\rightarrow
0^{+}}[\partial_{t}^{n}\partial_{r^{*}}^{m}\psi(t_{0},r_{0}^{*}+\varepsilon)\\\
-\partial_{t}^{n}\partial_{r^{*}}^{m}\psi(t_{0},r_{0}^{*}-\varepsilon)]$ (161)

to denote the jump in $\partial^{n}_{t}\partial^{m}_{r^{*}}\psi$ across the
world line. I only calculate jump conditions in the $r^{*}$ direction up to
$\left[\partial_{r^{*}}\psi\right]$, which I find by substituting the ansatz

$\displaystyle\psi$ $\displaystyle=\psi_{<}(t,r^{*})\theta(r^{*}_{0}-r^{*})$
$\displaystyle\quad+\psi_{>}(t,r^{*})\theta(r^{*}-r^{*}_{0})$ (162)

into Eqs. (21b) – (21f) and its $t$ and $r^{*}$ derivatives. Demanding in each
step that the singularity structure on the left hand side matches that of the
sources (and their derivatives) on the right hand side yields the jump
conditions

$\displaystyle\left[\psi\right]$
$\displaystyle=\frac{F_{\psi}}{f_{0}[(\partial_{t}r^{*}_{0})^{2}-1]}\text{,}$
(163) and $\displaystyle\left[\partial_{r^{*}}\psi\right]$
$\displaystyle=-\frac{G_{\psi}}{(\partial_{t}r^{*}_{0})^{2}-1}$
$\displaystyle\quad-\frac{\partial_{t}^{2}r^{*}_{0}\left[3\,(\partial_{t}r^{*}_{0})^{2}+1\right]F_{\psi}}{f_{0}\,[(\partial_{t}r^{*}_{0})^{2}-1]^{3}}$
$\displaystyle\quad+2\frac{\partial_{t}r^{*}_{0}\,\partial_{t}\left(F_{\psi}/f_{0}\right)}{[(\partial_{t}r^{*}_{0})^{2}-1]^{2}}\text{,}$
(164)

where $\psi$ stands for either one of $\psi$, $\chi$, or $\xi$.

###  D.2 Vector potential calculation

Since the source term in Eq. (147) is singular, the field is only continuous
across the world line of the particle, but not smooth. I use

$\left[\partial_{t}^{n}\partial_{r^{*}}^{m}\psi\right]=\lim_{\varepsilon\rightarrow
0^{+}}[\partial_{t}^{n}\partial_{r^{*}}^{m}\psi(t_{0},r_{0}^{*}+\varepsilon)\\\
-\partial_{t}^{n}\partial_{r^{*}}^{m}\psi(t_{0},r_{0}^{*}-\varepsilon)]$ (165)

to denote the jump in $\partial^{n}_{t}\partial^{m}_{r^{*}}\psi$ across the
world line. For my purposes I only need the jump conditions in the $r^{*}$
direction up to $\left[\partial^{2}_{r^{*}}\psi\right]$, which I find by
substituting the ansatz

$\displaystyle A^{\ell m}_{a}(t,r^{*})$ $\displaystyle=A^{\ell
m}_{a,<}(t,r^{*})\theta(r^{*}_{0}-r^{*})$ $\displaystyle\quad+A^{\ell
m}_{a,>}(t,r^{*})\theta(r^{*}-r^{*}_{0})\text{,}$ (166) $\displaystyle v^{\ell
m}(t,r^{*})$ $\displaystyle=v^{\ell m}_{<}(t,r^{*})\theta(r^{*}_{0}-r^{*})$
$\displaystyle\quad+v^{\ell m}_{>}(t,r^{*})\theta(r^{*}-r^{*}_{0})\text{,}$
(167) $\displaystyle\tilde{v}^{\ell m}(t,r^{*})$
$\displaystyle=\tilde{v}^{\ell m}_{<}(t,r^{*})\theta(r^{*}_{0}-r^{*})$
$\displaystyle\quad+\tilde{v}^{\ell m}_{>}(t,r^{*})\theta(r^{*}-r^{*}_{0})$
(168)

into Eqs. (148) – (151) and its $t$ and $r^{*}$ derivatives. Demanding in each
step that the singularity structure on the left hand side matches that of the
sources (and their derivatives) on the right hand side yields the jump
conditions

$\displaystyle\left[A^{\ell m}_{a}\right]=\left[w^{\ell m}\right]=0\text{,}$
(169) $\displaystyle\left[\partial_{r^{*}}A^{\ell
m}_{a}\right]=\frac{E^{2}}{E^{2}-\dot{r}_{0}^{2}}S_{a}\text{,}$ (170)
$\displaystyle\left[\partial_{r^{*}}w^{\ell
m}\right]=\frac{E^{2}}{E^{2}-\dot{r}_{0}^{2}}S_{\text{even/odd}}\text{,}$
(171) $\displaystyle\begin{split}\left[\partial^{2}_{r^{*}}A^{\ell
m}_{a}\right]=\left(\frac{2ME^{4}}{r_{0}^{2}(E^{2}-\dot{r}_{0}^{2})^{2}}-f_{0}\frac{(3\dot{r}_{0}^{2}+E^{2})E^{2}\ddot{r}_{0}}{(E^{2}-\dot{r}_{0}^{2})^{3}}\right)S_{a}\\\
\mbox{}+\frac{2ME^{3}\dot{r}_{0}}{r_{0}^{2}(E^{2}-\dot{r}_{0}^{2})^{2}}S_{b}-f_{0}\frac{2E^{2}\dot{r}_{0}}{(E^{2}-\dot{r}_{0}^{2})^{2}}\dot{S}_{a}\\\
\text{,}\end{split}$ (172)
$\displaystyle\begin{split}\left[\partial^{2}_{r^{*}}w^{\ell
m}\right]=-f_{0}\frac{(3\dot{r}_{0}^{2}+E^{2})E^{2}\ddot{r}_{0}}{(E^{2}-\dot{r}_{0}^{2})^{3}}S_{\text{even/odd}}\\\
\mbox{}-f_{0}\frac{2E^{2}\dot{r}_{0}}{(E^{2}-\dot{r}_{0}^{2})^{2}}\dot{S}_{\text{even/odd}}\text{,}\end{split}$
(173)

where $a,b\in\\{t,r^{*}\\},a\neq b$, $w\in\\{v,\tilde{v}\\}$.

## References

  * Haas (2007) R. Haas, Phys. Rev. D 75, 124011 (2007), eprint 0704.0797.
  * Haas and Poisson (2006) R. Haas and E. Poisson, Phys. Rev. D 74, 044009 (pages 29) (2006), eprint gr-qc/0605077, URL http://link.aps.org/abstract/PRD/v74/e044009.
  * Barack and Ori (2000) L. Barack and A. Ori, Phys. Rev. D 61, 061502 (2000), eprint gr-qc/9912010.
  * Vega and Detweiler (2008) I. Vega and S. Detweiler, Phys. Rev. D 77, 084008 (2008), eprint 0712.4405.
  * Barack et al. (2007) L. Barack, D. A. Golbourn, and N. Sago, Phys. Rev. D76, 124036 (2007), eprint 0709.4588.
  * Pound and Poisson (2008) A. Pound and E. Poisson, Phys. Rev. D 77, 044012 (2008), eprint 0708.3037.
  * Barack and Sago (2011) L. Barack and N. Sago, Phys. Rev. D83, 084023 (2011), eprint 1101.3331.
  * Warburton et al. (2011) N. Warburton, S. Akcay, L. Barack, J. R. Gair, and N. Sago (2011), * Temporary entry *, eprint 1111.6908.
  * Misner et al. (1973) C. W. Misner, K. S. Thorne, and J. A. Wheeler, _Gravitation_ (W. H. Freeman, San Francisco, 1973), ISBN 0716703343.
  * Lousto and Price (1997) C. O. Lousto and R. H. Price, Phys. Rev. D 56, 6439 (1997), eprint gr-qc/9705071.
  * Regge (1957) T. Regge, Phys. Rev. 108, 1063 (1957).
  * Martel and Poisson (2005) K. Martel and E. Poisson, Phys. Rev. D 71, 104003 (2005), eprint gr-qc/0502028.
  * Cunningham and Price (1979) C. T. Cunningham and R. H. Price, The Astrophysical Journal 230, 870 (1979).
  * (14) N. Sago, talk in 10th Capra meeting at UAH (2007).
  * Detweiler et al. (2003) S. Detweiler, E. Messaritaki, and B. F. Whiting, Phys. Rev. D 67, 104016 (2003), eprint gr-qc/0205079.
  * Diaz-Rivera et al. (2004) L. M. Diaz-Rivera, E. Messaritaki, B. F. Whiting, and S. Detweiler, Physical Review D (Particles, Fields, Gravitation, and Cosmology) 70, 124018 (pages 14) (2004), eprint gr-qc/0410011, URL http://link.aps.org/abstract/PRD/v70/e124018.
  * Detweiler (2008a) S. L. Detweiler, Phys. Rev. D77, 124026 (2008a), eprint 0804.3529.
  * Sago et al. (2008) N. Sago, L. Barack, and S. L. Detweiler, Phys. Rev. D78, 124024 (2008), eprint 0810.2530.
  * Barack and Sago (2009) L. Barack and N. Sago, Phys. Rev. Lett. 102, 191101 (2009), eprint 0902.0573.
  * Pfenning and Poisson (2002) M. J. Pfenning and E. Poisson, Phys. Rev. D 65, 084001 (2002), eprint gr-qc/0012057.
  * DeWitt-Morette and Dewitt (1964) C. DeWitt-Morette and B. S. Dewitt, _Relativite, groups et topologie = Relativity, groups and topology : lectures delivered at les Houches during the 1963 session of the Summer School of_ (Gordon and Breach, New York, 1964).
  * Detweiler (2008b) S. Detweiler, Phys. Rev. D 77, 124026 (2008b), eprint 0804.3529.
  * Poisson (2004) E. Poisson, Living Reviews in Relativity 7 (2004), URL http://www.livingreviews.org/lrr-2004-6.
  * Gundlach et al. (2005) C. Gundlach, J. M. Martin-Garcia, G. Calabrese, and I. Hinder, Class. Quant. Grav. 22, 3767 (2005), eprint gr-qc/0504114.
  * Barack and Lousto (2005) L. Barack and C. O. Lousto, Phys. Rev. D 72, 104026 (2005), eprint gr-qc/0510019.