Scalar self-force on eccentric geodesics in Schwarzschild spacetime: a time-domain computation

TL;DR

This paper presents a numerical time-domain method to compute the scalar self-force on particles in eccentric orbits around Schwarzschild black holes, using mode-sum regularization and high-order finite differences.

Contribution

It introduces a fourth-order convergent finite-difference scheme combined with mode-sum regularization for accurate self-force calculations on eccentric geodesics.

Findings

Validated the numerical method with various tests.

Provided results for mildly eccentric and zoom-whirl orbits.

Demonstrated the method's accuracy and applicability.

Abstract

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

Figures

Tables

Table 1. Table 1: Results for the self-force on a scalar particle with scalar charge q 𝑞 q on a circular orbit at r 0 = 6 ​ M subscript 𝑟 0 6 𝑀 r_{0}=6M . The first column lists the results as calculated in this work using time-domain numerical methods, while the second and third columns list the results as calculated in  Haas and Poisson ( 2006 ); Diaz-Rivera et al. ( 2004 ) using frequency-domain methods. For the t 𝑡 t and ϕ italic-ϕ \phi components the number of digits is limited by numerical roundoff error. For the r 𝑟 r component the number of digits is limited by the truncation error of the sum of multipole coefficients.

	This work:	Previous work:	Diaz-Rivera
	time-domain	frequency-domain Haas and Poisson (2006)	et. al. Diaz-Rivera et al. (2004)
$\frac{M^2}{q}{\Phi }_t^R$	$3.60339\times {10}^{-4}$	$3.60907254\times {10}^{-4}$
$\frac{M^2}{q}{\Phi }_r^R$	$1.6767\times {10}^{-4}$	$1.67730\times {10}^{-4}$	$1.6772834\times {10}^{-4}$
$\frac{M}{q}{\Phi }_{\varphi }^R$	$-5.30424\times {10}^{-3}$	$-5.30423170\times {10}^{-3}$

Table 2. Table 2: Estimated values for the various errors in the components of the self-force as described in the text. We show the truncation and discretization errors for a mildly eccentric orbit and the total error for a circular orbit. The truncation error is calculated using a plot similar to the one shown in Fig.  16 . The discretization error is estimated using a plot similar to that in Fig.  7 for the ℓ = 2 ℓ 2 \ell=2 , m = 2 𝑚 2 m=2 mode, and the total error is estimated as the difference between our values and those of  Haas and Poisson ( 2006 ) . We use p = 7.2 𝑝 7.2 p=7.2 , e = 0.5 𝑒 0.5 e=0.5 for the mildly eccentric orbit. Note that we use the tetrad component Φ ( + ) subscript Φ \Phi_{(+)} for the truncation error and the vector component ∂ r Φ subscript 𝑟 Φ \partial_{r}\Phi for the discretization error. Both are related by the translation table Eqs. ( 60 ) – ( 63 ), we expect corresponding errors to be comparable for Φ ( + ) subscript Φ \Phi_{(+)} and ∂ r Φ subscript 𝑟 Φ \partial_{r}\Phi .

error estimation	mildly eccentric orbit
truncation error ($\frac{M^2}{q}{\Phi }_{(+)}$)	$\approx 2\times {10}^{-3}\%$
discretization error ($\frac{M^2}{q}{\partial }_r{\Phi }_{\ell m}$)	$\approx {10}^{-5}\%$
comparison with reference values	circular orbit
$\frac{M^2}{q^2}{F}_t$	$0.2\%$
$\frac{M^2}{q^2}{F}_r$	$0.04\%$
$\frac{M}{q^2}{F}_{\varphi }$	$2\times {10}^{-4}\%$

Table 3. Table 3: Jump conditions for the derivatives purely in the u 𝑢 u or v 𝑣 v directions. r ˙ ˙ 𝑟 \dot{r} and r ¨ ¨ 𝑟 \ddot{r} are the particle’s radial velocity and acceleration, respectively. They are obtained from the equation of motion for the particle. ξ ¯ := E − r ˙ f assign ¯ 𝜉 𝐸 ˙ 𝑟 𝑓 \bar{\xi}:=\frac{E-\dot{r}}{f} and ξ := E + r ˙ f assign 𝜉 𝐸 ˙ 𝑟 𝑓 \xi:=\frac{E+\dot{r}}{f} were introduced for notational convenience. Quantities bearing a subscript p 𝑝 p are evaluated on the particle’s world line, while quantities bearing a subscript 0 0 are evaluated at the particle’s current position. Derivatives of V 𝑉 V with respect to either u 𝑢 u or v 𝑣 v are evaluated as ∂ u V = − 1 2 ​ f ​ ∂ r V subscript 𝑢 𝑉 1 2 𝑓 subscript 𝑟 𝑉 \partial_{u}V=-\frac{1}{2}f\partial_{r}V and ∂ v V = 1 2 ​ f ​ ∂ r V subscript 𝑣 𝑉 1 2 𝑓 subscript 𝑟 𝑉 \partial_{v}V=\frac{1}{2}f\partial_{r}V , respectively.

	$\left[\psi \right]=$	$0$
	$\left[{\partial }_u\psi \right]=$	$-\frac{1}{4}{\overline{\xi }}_0^{-1}\hat{S}\left({\tau }_0\right)\text{, }\left[{\partial }_v\psi \right]=\frac{1}{4}{\xi }_0^{-1}\hat{S}\left({\tau }_0\right)$
	$\left[{\partial }_u^2\psi \right]=$	$-\frac{1}{4}{\overline{\xi }}_0^{-1}\frac{d}{d\tau }{\left({\overline{\xi }}_p^{-1}\hat{S}\left(\tau \right)\right)}_{|\tau ={\tau }_0}$
	$\left[{\partial }_v^2\psi \right]=$	$\frac{1}{4}{\xi }_0^{-1}\frac{d}{d\tau }{\left({\xi }_p^{-1}\hat{S}\left(\tau \right)\right)}_{|\tau ={\tau }_0}$
	$\left[{\partial }_u^3\psi \right]=$	$\frac{1}{4}V{\xi }_0{\overline{\xi }}_0^{-1}\left[{\partial }_u\psi \right]-\frac{1}{4}{\overline{\xi }}_0^{-1}\frac{d}{d\tau }{\left[{\overline{\xi }}_p^{-1}\frac{d}{d\tau }\left({\overline{\xi }}_p^{-1}\hat{S}\left(\tau \right)\right)\right]}_{|\tau ={\tau }_0}$
	$\left[{\partial }_v^3\psi \right]=$	$\frac{1}{4}V{\overline{\xi }}_0{\xi }_0^{-1}\left[{\partial }_v\psi \right]+\frac{1}{4}{\xi }_0^{-1}\frac{d}{d\tau }{\left[{\xi }_p^{-1}\frac{d}{d\tau }\left({\xi }_p^{-1}\hat{S}\left(\tau \right)\right)\right]}_{|\tau ={\tau }_0}$
	$\left[{\partial }_u^4\psi \right]=$	$-\frac{1}{4}[-\frac{1}{2}{\overline{\xi }}_0^{-1}V\frac{{\ddot{r}}_0}{E}+\frac{1}{2}{\overline{\xi }}_0^{-1}V\frac{d}{d\tau }{\left(\frac{f_p}{E}{\xi }_p^2{\overline{\xi }}_p^{-1}\right)}_{|\tau ={\tau }_0}$
		$+3{\xi }_0{\overline{\xi }}_0^{-1}{\partial }_{u}V+{\xi }_0^2{\overline{\xi }}_0^{-2}{\partial }_{v}V\left]\right[{\partial }_u\psi ]+\frac{1}{2}\xi {}_0\overline{\xi }{}_0{}^{-1}V[{\partial }_u^2\psi ]$
		$-\frac{1}{4}{\overline{\xi }}_0^{-1}\frac{d}{d\tau }{\left({\overline{\xi }}_p^{-1}\frac{d}{d\tau }\left\{{\overline{\xi }}_p^{-1}\frac{d}{d\tau }\left[{\overline{\xi }}_p^{-1}\hat{S}\left(\tau \right)\right]\right\}\right)}_{|\tau ={\tau }_0}$
	$\left[{\partial }_v^4\psi \right]=$	$\frac{1}{4}[-\frac{1}{2}{\xi }_0^{-1}V\frac{{\ddot{r}}_0}{E}+\frac{1}{2}{\xi }_0^{-1}V\frac{d}{d\tau }{\left(\frac{f_p}{E}{\overline{\xi }}_p^2{\xi }_p^{-1}\right)}_{|\tau ={\tau }_0}$
		$+3{\overline{\xi }}_0{\xi }_0^{-1}{\partial }_{v}V+{\overline{\xi }}_0^2{\xi }_0^{-2}{\partial }_{u}V\left]\right[{\partial }_v\psi ]-\frac{1}{2}\overline{\xi }{}_0\xi {}_0{}^{-1}V[{\partial }_v^2\psi ]$
		$-\frac{1}{4}{\xi }_0^{-1}\frac{d}{d\tau }{\left({\xi }_p^{-1}\frac{d}{d\tau }\left\{{\xi }_p^{-1}\frac{d}{d\tau }\left[{\xi }_p^{-1}\hat{S}\left(\tau \right)\right]\right\}\right)}_{|\tau ={\tau }_0}$

Equations

$$\mathrm{d}{s}^2$$

$$g^{\alpha \beta }{\nabla }_{\alpha }{\nabla }_{\beta }\Phi (x)$$

$$\mu (x)$$

$$F_{\alpha }(\tau )$$

$$\Phi$$

$$F_{\alpha }$$

$${\Phi }_{(\mu )}(t,r,\theta ,\varphi )=\sum _{\ell ,m}{\Phi }_{(\mu )}^{\ell m}(t,r)Y_{\ell m}(\theta ,\varphi )\text{.}$$

$${\Phi }_{(\mu )}^S=$$

$${\Phi }_{(\mu ),\ell }^S$$

$$+\frac{D_{(\mu )}}{(\ell -\frac{1}{2})(\ell +\frac{3}{2})}+\cdots \text{,}$$

$${\Phi }_{(\mu )}^R$$

$$\Phi (x)$$

$$-{\partial }_t^2{\psi }_{\ell m}$$

$$-4\pi q\frac{f_0}{r_{0}E}{\bar{Y}}_{\ell m}(\pi /2,{\varphi }_0)\delta (r^{\ast }-r_0^{\ast })\text{,}$$

$$V_{\ell }=f\left(\frac{2M}{r^3}+\frac{\ell (\ell +1)}{r^2}\right)\text{.}$$

$$u^{\beta }{\nabla }_{\beta }{u}^{\alpha }$$

$${\Phi }_{(\mu )\ell }$$

$$\int {\int }_{\text{cell}}-4{\partial }_u{\partial }_v\psi \mathrm{d}u\mathrm{d}v=$$

$$\text{mbox}-\psi (t,r^{\ast }-h)-\psi (t,r^{\ast }+h)]\text{.}$$

$$G(t,r^{\ast })=-4\pi q\frac{f}{Er}{\bar{Y}}_{\ell m}(\pi /2,{\varphi }_0)\text{,}$$

$$\int {\int }_{\text{cell}}G\delta (r^{\ast }-r_0^{\ast }(t))\mathrm{d}u\mathrm{d}v=$$

$$\times {\bar{Y}}_{\ell m}(\pi /2,{\varphi }_0(t))\mathrm{d}t\text{,}$$

$${\int \int }_{\text{cell}}V\psi \mathrm{d}u\mathrm{d}v$$

$$g(t,r^{\ast })$$

$${\int \int }_{\text{cell}}g\mathrm{d}u\mathrm{d}v=$$

$$\text{mbox}{g}_{24}+g_{34}+g_{13})+16g_0]+O(h^6).$$

$$g_0=$$

$$\displaystyle+V_{10}\,\psi_{10}+V_{7}\,\psi_{7}-V_{9}\,\psi_{9}-V_{8}\,\psi_{8}\bigl{.}\bigr{]}+O(h^{4})\text{.}$$

$$g_{13}+g_{12}=$$

$$\text{mbox}\times \left[1-\frac{1}{2}{\left(\frac{h}{2}\right)}^{2}V(r_0^{\ast }-h/2)\right]+O(h^4)\text{,}$$

$$g_{24}+g_{34}=$$

$$\text{mbox}\times \left[1-\frac{1}{2}{\left(\frac{h}{2}\right)}^{2}V(r_0^{\ast }+h/2)\right]+O(h^4)$$

$$g_{13}+$$

$$+V(r_0^{\ast }-h/2){\psi }_1\left[1-\frac{1}{2}{\left(\frac{h}{2}\right)}^{2}V(r_0^{\ast }-h/2)\right]$$

$$+V(r_0^{\ast }+h/2){\psi }_4\left[1-\frac{1}{2}{\left(\frac{h}{2}\right)}^{2}V(r_0^{\ast }+h/2)\right]$$

$$\displaystyle+\frac{1}{2}\bigl{[}V(r_{0}^{*}-h/2)-2V(r_{0}^{*})$$

$$\displaystyle+V(r_{0}^{*}+h/2)\bigr{]}\,(\psi_{1}+\psi_{4})+O(h^{4})\text{.}$$

$${\psi }_3=$$

Full text

I Introduction

The inspiral and capture of solar-mass compact objects by supermassive black holes is one of the most promising and interesting sources of gravitational radiation to be detected by the future space-based gravitational-wave antenna LISA LIS . For these extreme mass-ratio inspirals, one can treat the compact object as a point mass and describe its influence on the spacetime perturbatively. Going beyond the test mass limit, its motion is no longer along a geodesic of the unperturbed spacetime of the central black hole; it is a geodesic of the perturbed spacetime created by the presence of the moving body. When viewed from the unperturbed spacetime, the small body is said to move under the influence of its gravitational self-force. The self-force induces radiative losses of energy and angular momentum, which will eventually drive the object into the black hole. To describe the motion of the body, including its inspiral toward the black hole, we seek to evaluate the self-force and calculate its effect on the motion. One way of doing this uses the mode-sum regularization procedure introduced by Barack and Ori Barack and Ori (2000). (For a comprehensive introduction of the problem, see the special issue of Classical and Quantum Gravity Lousto (2005a).)

In this paper, in an effort to build expertise to calculate the gravitational self-force, we retreat to the technically simpler problem of a point particle of mass $m$ endowed with a scalar charge $q$ orbiting a Schwarzschild black hole of mass $M$. Following up on a previous paper Haas and Poisson (2006), we implement the numerical part of the regularization procedure for generic orbits with a time-domain integration of the scalar-wave equation.

I.1 The problem

Our goal is to calculate the regularized self-force acting on a scalar point charge in orbit around a Schwarzschild black hole. In analogy with the gravitational case, where in a first-order (in $m/M$) perturbative calculation the particle moves on a geodesic of the background spacetime, we take the orbit of the particle to be a geodesic and calculate the self-force as a vector field on this geodesic. We start by writing the Schwarzschild metric using the tortoise coordinate $r^{*}=r+2M\ln\left(\frac{r}{2M}-1\right)$ as

[TABLE]

where $f=\left(1-\frac{2M}{r}\right)$, $\mathrm{d}\Omega^{2}=\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. Our task is to solve the scalar wave equation

[TABLE]

where $\nabla_{\alpha}$ is the covariant derivative compatible with the metric $g_{\alpha\beta}$, $\Phi(x)$ is the scalar field created by a scalar charge $q$ which moves along a world line $\gamma:\tau\mapsto z(\tau)$ parametrized by proper time $\tau$. The source term $\mu(x)$ appearing on the right-hand side is written in terms of a scalarized four-dimensional Dirac $\delta$-function $\delta_{4}(x,x^{\prime}):=\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})}$. Because of the singularity in the source term, the retarded solution to Eq. (2) is singular on the world line, and the naïve expression for the self-force,

[TABLE]

must be regularized. Following DeWitt and Brehme DeWitt and Brehme (1960), Mino, Sasaki, Tanaka Mino et al. (1997), Quinn and Wald Quinn and Wald (1997), Quinn Quinn (2000) carried out this regularization for the electromagnetic, scalar and gravitational radiation reaction. In later work, Detweiler and Whiting Detweiler and Whiting (2003) introduced a very useful decomposition of the retarded solution of Eq. (2) in terms of a singular part $\Phi^{S}$ and a regular remainder $\Phi^{R}$:

[TABLE]

$\Phi^{R}$ is regular and differentiable at the position of the particle, satisfies the homogeneous wave equation associated with Eq. (2), and is solely responsible for the self-force acting on the particle. $\Phi^{S}$, on the other hand, satisfies Eq. (2), is just as singular at the particle’s position as the retarded solution, and produces no force on the particle. Rearranging Eq. (5) and differentiating once, we can write the regularized self-force as

[TABLE]

In a previous paper Haas and Poisson (2006), we described our implementation of the regularization procedure to find a mode-sum representation of $\nabla_{\alpha}\Phi^{S}$ along a generic geodesic of the Schwarzschild spacetime. Schematically, we introduce a tetrad $e^{\alpha}_{(\mu)}$ and decompose the tetrad components $\Phi_{(\mu)}:=e^{\alpha}_{(\mu)}\nabla_{\alpha}\Phi$ of the field gradient in terms of ordinary scalar spherical harmonics $Y_{\ell m}$:

[TABLE]

Each mode $\Phi^{\ell m}_{(\mu)}(t,r)$ is finite at the position of the particle, but their sum diverges on the world line. In Haas and Poisson (2006), we derive analytic expressions for the mode-sum decomposition of $\Phi^{S}_{(\mu)}$,

[TABLE]

where the coefficients $A_{(\mu)}$, $B_{(\mu)}$, $C_{(\mu)}$, and $D_{(\mu)}$ are independent of $\ell$; they are listed in Appendix B for convenience.

As each mode of $\Phi$ is finite, it is straightforward to compute the modes of the retarded solution using numerical methods, and we will describe how this was done in Sec. IV. We use the numerical solutions in Eq. (6) to calculate the regularized self-force, regularizing mode-by-mode:

[TABLE]

where $\Phi_{(\mu),\ell}:=\sum_{m}\Phi^{\ell m}_{(\mu)}Y_{\ell m}$ (no summation over $\ell$ implied).

For numerical purposes it is convenient to define $\psi_{\ell m}$ by

[TABLE]

where $Y_{\ell m}$ are the usual scalar spherical harmonics. After substituting in Eq. (2), this yields a reduced wave equation for the multipole moments $\psi_{\ell m}$:

[TABLE]

where

[TABLE]

An overbar denotes complex conjugation, $E=-u_{t}$ is the particle’s conserved energy per unit mass, and $u^{\alpha}=\frac{\mathrm{d}z^{\alpha}}{\mathrm{d}\tau}$ is its four velocity. Quantities bearing a subscript “[math]” are evaluated at the particle’s position; they are functions of $\tau$ that are obtained by solving the geodesic equation

[TABLE]

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

Once we have numerically solved Eq. (12), we extract numerical estimates for $\psi_{\ell m}$, $\partial_{t}\psi_{\ell m}$ and $\partial_{r^{*}}\psi_{\ell m}$, which can then be used to find $\Phi_{\ell m}$, $\partial_{t}\Phi_{\ell m}$ and $\partial_{r}\Phi_{\ell m}$. These—together with the translation table displayed in Eqs. (1.23)–(1.26) of Haas and Poisson (2006), reproduced in Appendix A—allow us to find the tetrad components $\Phi_{(\mu)\ell m}$ with respect to the tetrad defined by Eqs. (1.18)–(1.21) of Haas and Poisson (2006). Eventually we regularize the multipole coefficients

[TABLE]

using Eq. (10); this involves the regularization parameters listed in Eqs. (1.30)–(1.45) of Haas and Poisson (2006), which are reproduced in Appendix B.

I.2 Organization of this paper

In Sec. II we introduce the main ideas behind the discretization scheme used in the numerical simulation. Sec. III describes the choices we make in order to handle the problems of specifying initial data and proper boundary conditions. The next section—Sec. IV—provides details on the concrete implementation of the ideas put forth in Secs. II and III. In Sec. V we describe the tests we performed in order to validate our implementation of the numerical method. Sec. VI finally presents sample results for a small number of representative simulations.

I.3 Future work

This work, which deals with a scalar charge moving in the Schwarzschild spacetime, is not intended to produce physically or astrophysically interesting results. Instead, its goal is to help us evaluate the merits of several strategies that could be used to tackle the more interesting (and difficult) problems of electromagnetism and gravity.

One future project we are currently exploring is to apply the formalism developed so far to the electromagnetic self-force acting on an electric charge. Beyond the technical complication of having to deal with a vector field instead of a single scalar quantity, we are also faced with the reality of having to impose a gauge (in our case: the Lorenz gauge) and to eliminate (or at least control) gauge violations in the numerical simulation. The first step, namely, the calculation of the regularization parameters $A_{(\mu)}$, $B_{(\mu)}$, $C_{(\mu)}$, and $D_{(\mu)}$ for the self-force, is currently underway. Also underway is the calculation of the regularization parameters for he gravitational self-force.

Another project is the implementation of a scheme to use the calculated self-force to update the orbital parameters of a particle on its inspiral toward the black hole. The standard proposed approach to this problem in the past has been to calculate the self-force on a set of geodesics which are momentarily tangent to the particle’s trajectory. The self-force calculated in this way is then used to update the orbital elements. This “after the fact” calculation of the motion requires one to build (in advance) a large database of self-force values for the anticipated set of orbital parameters that the particle’s trajectory will assume during its inspiral. Alternatively, and conceptually more simply, the self-force could be calculated self-consistently along the real, accelerated trajectory. Such an approach requires changes in the expressions of the regularization parameters, which so far have been derived only for geodesic orbits. We are currently investigating the merits of such an approach.

II Numerical method

In this section we describe the algorithm used to integrate the reduced wave equation [Eq. (12)] numerically. For the most part we use the fourth-order algorithm introduced by Lousto Lousto (2005b), with some modifications to suit our needs. We choose to implement a fourth-order convergent code because second-order convergence for the potential $\Phi$, while much easier to achieve, would guarantee only first-order convergence for $\nabla_{\alpha}\Phi$, the quantity in which we are ultimately interested. With a fourth-order convergent code we can expect to achieve third-order convergence for $\nabla_{\alpha}\Phi$, which is required for an accurate estimation of the self-force. Numerical experiments, however, show that in practice we do achieve fourth-order convergence for the derivatives of $\Phi$, a fortunate outcome that we exploit but cannot explain.

From now on, we will suppress the subscripts $\ell$ and $m$ on $V_{\ell}$ and $\psi_{\ell m}$ for convenience of notation. The wave equation consists of three parts: the wave-operator term $(\partial_{r^{*}}^{2}-\partial_{t}^{2})\psi$ and the potential term $V\psi$ on the left-hand side, and the source term on the right-hand side of the equation. Of these, the wave operator turns out to be easiest to handle, and the source term does not create a substantial difficulty. The term involving the potential $V$ turns out to be the most difficult one to handle.

Following Lousto we introduce a staggered grid with step sizes $\Delta t=\frac{1}{2}\Delta r^{*}\equiv h$, which follows the characteristic lines of the wave operator in Schwarzschild spacetime; see Fig. 1 for a sketch of a typical grid cell. The basic idea behind the method is to integrate the wave equation over a unit cell of the grid, which nicely deals with the Dirac-$\delta$ source term on the right-hand side. To this end, we introduce the Eddington-Finkelstein null coordinates $v=t+r^{*}$ and $u=t-r^{*}$ and use them as integration variables.

II.1 Differential operator

Rewriting the wave operator in terms of $u$ and $v$, we find $-\partial_{t}^{2}+\partial_{r^{*}}^{2}=-4\partial_{u}\partial_{v}$, which allows us to evaluate the integral involving the wave operator exactly. We find

[TABLE]

II.2 Source term

If we integrate over a cell traversed by the particle, then the source term on the right-hand side of the equation will have a non-zero contribution. Writing the source term as $G(t,r^{*})\delta\bm{(}r^{*}-r_{0}^{*}(t)\bm{)}$ with

[TABLE]

we find

[TABLE]

where $t_{1}$ and $t_{2}$ are the times at which the particle enters and leaves the cell, respectively. While we do not have an analytic expression for the trajectory of the particle (except when the particle follows a circular orbit), we can numerically integrate the first-order ordinary differential equations that govern the particle’s motion to a precision that is much higher than that of the partial differential equation governing $\psi$. In this sense we treat the integral over the source term as exact. To evaluate the integral we adopt a four-point Gauss-Legendre scheme, which has an error of order $h^{8}$.

II.3 Potential term

The most problematic term—from the point of view of implementing an approximation of sufficiently high order in $h$—turns out to be the term $V\psi$ in Eq. (12). Since this term does not contain a $\delta$-function, we have to approximate the double integral

[TABLE]

up to terms of order $h^{6}$ for a generic cell in order to achieve an overall $O(h^{4})$ convergence of the scheme.

Here we have to treat cells traversed by the particle (“sourced” cells) differently from the generic (“vacuum”) cells. While much of the algorithm can be transferred from the vacuum cells to the sourced cells, some modifications are required. We will describe each case separately in the following subsections.

II.3.1 Vacuum case

To implement Lousto’s algorithm to evolve the field across the vacuum cells, we use a double Simpson rule to compute the integral Eq. (19). We introduce the notation

[TABLE]

and label our points in the same manner (see Fig. 1) as in Lousto (2005b):

[TABLE]

Here, for example, $g_{1}$ is the value of $g$ at the grid point labeled $1$, and $g_{12}$ is the value of $g$ at the off-grid point labeled $12$, etc. Deviating from Lousto’s algorithm, we choose to calculate $g_{0}$ using an expression different from that derived in Lousto (2005b). Unlike Lousto’s approach, our expression exclusively involves points that are within the past light cone of the current cell. We find

[TABLE]

In order to evaluate the term in parentheses in Eq. (21), we again use a variant of the equations given in Lousto (2005b). Lousto’s equations (33) and (34),

[TABLE]

contain isolated occurrences of $\psi_{0}$, the value of the field at the central point. Since Eq. (22) only allows us to find $g_{0}=V_{0}\psi_{0}$, finding $\psi_{0}$ would involve a division by $V_{0}$, which will be numerically unstable very close to the event horizon where $V_{0}\approx 0$. Instead we choose to express the potential term appearing in the square brackets as a Taylor series around $r_{0}^{*}$. This allows us to eliminate the isolated occurrences of $\psi_{0}$, and we find

[TABLE]

Because of the $\left(\frac{h}{3}\right)^{2}$ factor in Eq. (21), this allows us to reach the required $O(h^{6})$ convergence for a generic vacuum cell. This—given that there is a number of order $N=1/h^{2}$ of such cells—yields the desired overall $O(h^{4})$ convergence of the full algorithm, at the end of the $N$ steps required to finish the simulation.

II.3.2 Sourced cells

For vacuum cells, the algorithm described above is the complete algorithm used to evolve the field forward in time. For cells traversed by the particle, however, we have to reconsider the assumptions used in deriving Eqs. (22) and (25). When deriving Eq. (25) we have employed the second-order evolution algorithm (see Lousto (2005b)), in which the single step equation

[TABLE]

is accurate only to $O(h^{3})$ for cells traversed by the particle. For these cells, therefore, the error term in Eq. (25) is $O(h^{3})$ instead of $O(h^{4})$. As there is a number of order $N^{\prime}=1/h$ of cells that are traversed by the particle in a simulation run, the overall error—after including the $\left(\frac{h}{3}\right)^{2}$ factor in Eq. (21)—is of order $h^{4}$. We can therefore afford this reduction of the convergence order in Eq. (25)

Equation (22), however, is accurate only to $O(h)$ for cells traversed by the particle. Again taking the $\left(\frac{h}{3}\right)^{2}$ factor into account, this renders the overall algorithm $O(h^{2})$.

Figure 2 shows the cells affected by the particle’s traversal and the reduced order of the single step equation for each cell. Cells whose convergence order is $O(h^{5})$ or higher do not need modifications, since there is only a number $N^{\prime}=1/h$ of such cells in the simulation. We are therefore concerned about cells neighboring the particle’s trajectory and those traversed by the particle.

Cells neighboring the particle

These cells are not traversed by the particle, but the particle might have traversed cells in their past light-cone, which are used in the calculation of $g_{0}$ in Eq. (22). For these cells, we use a one-dimensional Taylor expansion of $g(t,r^{*})$ within the current time-slice $t=t_{0}$,

[TABLE]

for the cell on the left-hand side, and

[TABLE]

for the cell on the right-hand side, where $(t_{0},r_{0}^{*})$ is the center of the cell traversed by the particle. Both of these are more accurate than is strictly necessary; we would need error terms of order $h^{3}$ to achieve the desired overall $O(h^{4})$ convergence of the algorithm. Keeping the extra terms, however, improves the numerical convergence slightly.

Cell traversed by the particle

We choose not to implement a fully explicit algorithm to handle cells traversed by the particle, because this would increase the complexity of the algorithm by a significant factor. Instead we use an iterative approach to evolve the field using the integrated wave equation

[TABLE]

In this equation the integral involving the source term can be evaluated to any desired accuracy at the beginning of the iteration, because the motion of the particle is determined by a simple system of ordinary differential equations, which are easily integrated with reliable numerical methods. It remains to evaluate the integral over the potential term, which we do iteratively. Schematically the method works as follows:

• •

Make an initial guess for $\psi_{3}$ using the second-order algorithm. This guess is correct up to terms of $O(h^{3})$.

• •

Match a second-order piecewise interpolation polynomial to the six points that make up the past light-cone of the future grid point, including the future point itself.

• •

Use this approximation for $\psi$ to numerically calculate

[TABLE]

using two-by-two point Gauss-Legendre rules for the six sub-parts indicated in Fig. 3.

• •

Update the future value of the field and repeat the process until the iteration has converged to a required degree of accuracy.

III Initial values and boundary conditions

As is typical for numerical simulations, we have to pay careful attention to specifying initial data and appropriate boundary conditions. These aspects of the numerical method are highly non-trivial problems in full numerical relativity, but they can be solved or circumvented with moderate effort in the present work.

III.1 Initial data

In this work we use a characteristic grid consisting of points lying on characteristic lines of the wave operator to evolve $\psi$ forward in time. As such, we need to specify characteristic initial data on the lines $u=u_{0}$ and $v=v_{0}$ shown in Fig. 4. We choose not to worry about specifying “correct” initial data, but instead arbitrarily choose $\psi$ to vanish on $u=u_{0}$ and $v=v_{0}$:

[TABLE]

This is equivalent to adding spurious initial waves in the form of a homogeneous solution of Eq. (12) to the correct solution. This produces an initial wave burst that moves away from the particle with the speed of light, and quickly leaves the numerical domain. Any remaining tails of the spurious initial data decay as $t^{-(2\ell+2)}$ as shown in Price (1972) and become negligible after a short time. We conclude that the influence of the initial-wave content on the self-force becomes negligible after a time of the order of the light-crossing time of the particle’s orbit.

III.2 Boundary conditions

On the analytical side we would like to impose ingoing boundary conditions at the event horizon $r^{*}\rightarrow-\infty$ and outgoing boundary conditions at spatial infinity $r^{*}\rightarrow\infty$, ie.

[TABLE]

Because of the finite resources available to a computer we can only simulate a finite region of the spacetime, and are faced with the reality of implementing boundary conditions at finite values of $r^{*}$. Two solutions to this problem present themselves:

choose the numerical domain to be the domain of dependence of the initial data surface. Since the effect of the boundary condition can only propagate forward in time with at most the speed of light, this effectively hides any influence of the boundary. This is what we choose to do in order to deal with the outer boundary condition. 2. 2.

implement boundary conditions sufficiently “far out” so that numerically there is no difference between imposing the boundary condition there or at infinity. Since the boundary conditions depend on the vanishing of the potential $V(r)$ appearing in the wave equation, this will happen once $1-2M/r\approx 0$. Near the horizon $r\approx 2M(1+\exp(r^{*}/2M))$, so this will happen—to numerical accuracy—for modestly large (negative) values of $r^{*}\approx-73\,M$. We choose to implement the ingoing waves condition $\partial_{u}\psi_{\ell m}=0$ there.

IV Implementation

Making more precise the ideas developed in the preceding sections, we implement the following numerical scheme.

IV.1 Particle motion

Following Darwin Darwin (1959) we introduce the dimensionless semi-latus rectum $p$ and the eccentricity $e$ such that for a bound orbit around a Schwarzschild black hole of mass $M$,

[TABLE]

are the radial positions of the periastron and apastron, respectively. Energy per unit mass and angular momentum per unit mass are then given by

[TABLE]

Together with these definitions it is useful to introduce an orbital parameter $\chi$ such that along the trajectory of the particle,

[TABLE]

where $\chi$ is single-valued along the orbit. We can then write down first-order differential equations for $\chi(t)$ and the azimuthal angle $\phi(t)$ of the particle,

[TABLE]

We 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. Intermediate values of the particle’s position are found using a Hermite interpolation of the nearest available calculated positions.

IV.2 Initial data

We do not specify initial data. The field is set to zero on the initial characteristic slices, $u=u_{0}$ and $v=v_{0}$.

IV.3 Boundary conditions

We adjust the outer boundary of the numerical domain at each time-step so that we cover the domain of dependence of the initial characteristic surfaces and the particle’s world line. The resulting numerical domain was already shown in Fig. 4.

Near the event horizon, at $r^{*}\approx-73\,M$, we implement an ingoing-wave boundary condition by imposing

[TABLE]

This allows us to drastically reduce the number of cells in the numerical domain, and consequently the running time of the simulation.

IV.4 Evolution in vacuum

Cells not traversed by the particle are evolved using Eqs. (16), (21) – (25). Explicitly written out, we use

[TABLE]

where $g_{0}$ is given by Eq. (22) and the sum $g_{12}+g_{24}+g_{34}+g_{13}$ is given by Eq. (25).

IV.5 Cells next to the particle

Vacuum cells close to the current position of the particle require a different approach to calculate $g_{0}$, since the cells in their past light cone could have been traversed by the particle. We use Eqs. (27) and (28) to find $g_{0}$ in this case. Other than this modification, the same algorithm as for generic vacuum cells is used.

IV.6 Cells traversed by the particle

We evolve cells traversed by the particle using the iterative algorithm described in Sec. II.3.2. Here

[TABLE]

where the initial guess for the iterative evolution of $\int\!\!\int_{\text{cell}}V\psi\,\mathrm{d}u\,\mathrm{d}v$ is obtained using the second order algorithm of Lousto and Price Lousto and Price (1997),

[TABLE]

Successive iterations use a four-point Gauss-Legendre rule to evaluate the integral of $V\psi$; this requires a second-order polynomial interpolation of the current field values as described in Appendix C.

IV.7 Extraction of the field data at the

particle

In order to extract the value of the field and its first derivatives at the position of the particle, we again use a polynomial interpolation at the points surrounding the particle’s position. Using a fourth-order polynomial, as described in Appendix C, we can estimate $\psi$, $\partial_{t}\psi_{t}$, and $\partial_{r^{*}}\psi$ at the position of the particle up to errors of order $h^{4}$. As was briefly mentioned in Sec. II, we would expect an error term of order $h^{3}$ for $\partial_{t}\psi_{t}$ and $\partial_{r^{*}}\psi$. The $O(h^{4})$ accuracy we actually achieve by using a fourth-order (instead of a third-order) piecewise polynomial shows up clearly in a regression plot such as Fig. 7.

IV.8 Regularization of the mode sum

We use the calculated multipole moments $\psi_{\ell m}$ to construct the multipole moments $\Phi_{\ell m}$, and first derivatives $\partial_{t}\Phi_{\ell m}$ and $\partial_{r}\Phi_{\ell m}$, of the scalar field. These, in turn, are used to calculate the tetrad components $\Phi_{(0)\ell m}$, $\Phi_{(+)\ell m}$, $\Phi_{(-)\ell m}$, and $\Phi_{(3)\ell m}$ of the field gradient according to Eqs. (1.23)–(1.26) of Haas and Poisson (2006), which are reproduced in Appendix A. These multipoles then give rise to the multipole coefficients of the retarded field,

[TABLE]

which are subjected to the regularization procedure described by Eq. (1.29) of Haas and Poisson (2006),

[TABLE]

using the regularization parameters $A_{(\mu)}$, $B_{(\mu)}$, $C_{(\mu)}$, and $D_{(\mu)}$ tabulated in Appendix B.

Finally we reconstruct the vector components of the field gradient using Eqs. (1.47)–(1.48) of Haas and Poisson (2006),

[TABLE]

and calculate the self-force

[TABLE]

We recall the discussion in Sec. I.1 concerning the definition of $\Phi^{R}$, its connection to the self-force acting on the particle, and its regularity at the particle’s position.

V Numerical tests

In this section we present the tests we have performed to validate our numerical evolution code. First, in order to check the fourth-order convergence rate of the code, we perform regression runs with increasing resolution for both a vacuum test case, where we seeded the evolution with a Gaussian wave packet, and a case where a particle is present. As a second test, we compute the regularized self-force for several different combinations of orbital elements $p$ and $e$ and check that the multipole coefficients decay with $\ell$ as expected. This provides a very sensitive check on the overall implementation of the numerical scheme, as well as the analytical calculations that lead to the regularization parameters. Finally, we calculate the self-force for a particle on a circular orbit and show that it agrees with the results presented in Haas and Poisson (2006); Diaz-Rivera et al. (2004).

V.1 Convergence tests: Vacuum

As a first test of the validity of our numerical code we estimate the convergence order by removing the particle and performing regression runs for several resolutions. We use a Gaussian wave packet as initial data,

[TABLE]

where $v_{p}=75\,M$ and $\sigma=10\,M$, $v_{0}=-u_{0}=6\,M+2\,M\ln 2$, and we extract the field values at $r^{*}=20\,M$. Several such runs were performed, with varying resolution of 2, 4, 8, 16, and 32 grid points per $M$. Figure 5 shows $\psi(2h)-\psi(h)$ rescaled by appropriate powers of $2$, so that in the case of fourth-order convergence the curves would lie on top of each other. As can be seen from the plots, they do, and the vacuum portion of the code is indeed fourth-order convergent.

V.2 Convergence tests: Particle

While the convergence test described in section V.1 clearly shows that the desired convergence is achieved for vacuum evolution, it does not test the parts of the code that are used in the integration of the inhomogeneous wave equation. To test these we perform a second set of regression runs, this time using a non-zero charge $q$. We extract the field at the position of the particle, thus also testing the implementation of the extraction algorithm described in section IV.7. For this test we 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. 6 the convergence is still of fourth order, but the two curves no longer lie precisely on top of each other at all times.

The region before $t\approx 100\,M$ is dominated by the initial wave burst and therefore does not scale as expected, yielding two very different curves. In the region $300\,M\lesssim t\lesssim 400\,M$ the two curves lie on top of each other, as expected for a fourth-order convergent algorithm. In the region between $t\approx 200\,M$ and $t\approx 300\,M$, however, the dashed curves have slightly smaller amplitudes than the solid one, indicating an order of convergence different from (but close to) four.

To explain this behavior we have to examine the terms that contribute significantly to the error in the simulation. The numerical error is almost completely dominated by that of the approximation of the potential term $\int\!\!\int_{\text{cell}}V\psi\,\mathrm{d}u\,\mathrm{d}v$ in the integrated wave equation. For vacuum cells the error in this approximation scales as $h^{6}$, where $h$ is the step size. For cells traversed by the particle, on the other hand, the approximation error depends also on the difference $t_{2}-t_{1}$ of the times at which the particle enters and leaves the cell. This difference is bounded by $h$ but does not necessarily scale as $h$. For example, if a particle enters a cell at its very left, then scaling $h$ by $\frac{1}{2}$ would not change $t_{2}-t_{1}$ at all, thus leading to a scaling behavior that differs from expectation.

To investigate this further we conducted test runs of the simulation for a particle on a circular orbit at $r=6\,M$. In order to observe the expected scaling behavior, we have to make sure that the particle passes through the tips of the cell it traverses. When this is the case, then $t_{2}-t_{1}\equiv h$ and a plot similar to the one shown in Fig. 6 shows the proper scaling behavior. As a further test we artificially reduced the convergence order of the vacuum algorithm to two by implementing the second-order algorithm described in Lousto (2005b). By keeping the algorithm that deals with sourced cells unchanged, we reduced the relative impact on the numerical error. This, too, allows us to recover the expected (second-order) convergence. Figures 8 and 9 illustrate the effects of the measures taken to control the convergence behavior.

V.3 High-$\ell$ behavior of the multipole coefficients

Inspection of Eq. (42) reveals that a plot of $\Phi_{(\mu)\ell}$ as a function of $\ell$ (for a selected value of $t$) should display a linear growth in $\ell$ for large $\ell$. Removing the $A_{(\mu)}$ term should produce a constant curve, removing the $B_{(\mu)}$ term (given that $C_{(\mu)}=0$) should produce a curve that decays as $\ell^{-2}$, and finally, removing the $D_{(\mu)}$ term should produce a curve that decays as $\ell^{-4}$. It is a powerful test of the numerical methods to check whether these expectations are borne out by the numerical data. Fig. 10 plots the remainders as obtained from our numerical simulation, demonstrating the expected behavior.

It displays, on a logarithmic scale, the absolute value of $\operatorname{Re}\Phi^{\mathrm{R}}_{(+)\ell}$, the real part of the $(+)$ component of the self-force. The orbit is eccentric ($p=7.2$, $e=0.5$), and all components of the self-force require regularization. The first curve (in 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)}$; this clearly approaches a constant for large values of $\ell$. The curve made up of diamonds shows the behavior after removal of $B_{(\mu)}$; because $C_{(\mu)}=0$, it decays as $\ell^{-2}$, a behavior that is confirmed by a fit to the $\ell\geq 5$ part of the curve. Finally, after removal of $D_{(\mu)}/[(\ell-\frac{1}{2})\,(\ell+\frac{3}{2})]$ the terms of the sum decrease in magnitude as $\ell^{-4}$ for large values of $\ell$, as derived in Detweiler et al. (2003). Each one of the last two curves would result in a converging sum, but the convergence is much faster after subtracting the $D_{(\mu)}$ terms. We thereby gain more than 2 orders of magnitude in the accuracy of the estimated sum.

Figure 10 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. (42) to have a vastly different behavior.

V.4 Self-force on a circular orbit

For the case of a circular orbit, the regularization parameters $A_{(0)}$, $B_{(0)}$, and $D_{(0)}$ all vanish identically, so that the $(0)$ (or alternatively the $t$) component of the self-force does not require regularization. Figure 11 thus shows only one curve, with the magnitude of the multipole coefficients decaying exponentially with increasing $\ell$.

As a final test, in Table 1 we compare our result for the self-force on a particle in a circular orbit at $r=6M$ to those obtained in Haas and Poisson (2006); Diaz-Rivera et al. (2004) using a frequency-domain code.

For a circular orbit, a calculation in the frequency domain is more efficient, and we expect the results of Haas and Poisson (2006); Diaz-Rivera et al. (2004) to be much more accurate than our own results. This fact is reflected in the number of regularization coefficients we can reliably extract from the numerical data, before being limited by the accuracy of the numerical method: the frequency-domain calculation found usable multipole coefficients up to $\ell=20$, whereas our data for $\Phi^{\mathrm{R}}_{(0)\ell}$ is dominated by noise by the time $\ell$ reaches $16$. Figure 11 shows this behavior.

V.5 Accuracy of the numerical method

Several figures of merit can be used to estimate the accuracy of numerical values for the self-force.

An estimate for the truncation error arising from cutting short the summation in Eq. (42) at some $\ell_{\text{max}}$ can be calculated by considering the behavior 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

[TABLE]

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

[TABLE]

where $\ell_{\text{max}}$ is the value at which we cut the summation short. For all but the special case of the $(0)$ component for a circular orbit, for which all regularization parameters vanish identically, we use this approach to calculate an estimate for the truncation error.

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, we can estimate this discretization error by extrapolating the differences of simulations using different step sizes down to $h=0$. This is what was done in the graphs shown in Sec. V.2.

We display results for mildly eccentric orbits. A high eccentricity causes $\partial_{r}\Phi$ (displayed in Fig. 7) to be plagued by high frequency noise produced by effects similar to those described in Sec. V.2. This makes it impossible to reliably estimate the discretization error for these orbits. We do not expect this to be very different from the errors for mildly eccentric orbits.

Finally we compare our final results for the self-force $F_{\alpha}$ to “reference values”. For circular orbits, frequency-domain calculations are much more accurate than our time-domain computations. We thus compare our results to the results obtained in Haas and Poisson (2006). Table 2 lists typical values for the various errors listed above.

VI Sample results

In this section we describe some results of our numerical calculation.

VI.1 Mildly eccentric orbit

We 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 $1000\,M$ with a resolution of 16 grid points per $M$, both in the $t$ and $r^{*}$ directions, for $\ell=0$. Higher values of $\ell$ (and thus $m$) require a corresponding increase in the number of grid points used to achieve the same fractional accuracy. Multipole coefficients for $0\leq\ell\leq 15$ are calculated and used to reconstruct the regularized self-force $F_{\alpha}$ along the geodesic. Figure 13 shows the result of the calculation.

For the choice of parameters used to calculate the force shown in Fig. 13, the error bars corresponding to the truncation error (which are already much larger than than the discretization error) would be of the order of the line thickness and have not been drawn.

Already for this small eccentricity, we see that the self-force is most important when the particle is closest to the black hole (ie. for $200\,M\lesssim t\lesssim 400\,M$ and $600\,M\lesssim t\lesssim 800\,M$); the self-force acting on the particle is very small once the particle has moved away to $r\approx 15\,M$.

VI.2 Zoom-whirl orbit

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 black hole, emitting copious amounts of radiation. Figures 14 and 15 show the trajectory of a particle and the force on such an orbit with $p=7.8001$, $e=0.9$.

Even more so than for the mildly eccentric orbit discussed in Sec. VI.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.

Defining energy $E$ per unit mass and angular momentum $L$ per unit mass in the usual way,

[TABLE]

and following eg. the treatment of Wald Wald (1984), Appendix C, it is easy to see that the rates of change $\dot{E}$ and $\dot{L}$ (per unit proper time) are directly related to components of the acceleration $a_{\alpha}$ (and therefore force) experienced by the particle via

[TABLE]

The self-force shown in Fig. 15 therefore confirms our naïve expectation that the self-force should decrease both the energy and angular momentum of the particle as radiation is emitted.

It is instructive to have a closer look at the force acting on the particle when it is within the zoom phase, and also when it is moving around the black hole on the nearly circular orbit of the whirl phase. In Fig. 16 and Fig. 17 we show plots of $\Phi_{(0)\ell}$ vs. $\ell$ after the removal of the $A_{(\mu)}$, $B_{(\mu)}$, and $D_{(\mu)}$ terms. While the particle is still zooming in toward the black hole, $\Phi_{(0)\ell}$ behaves exactly as for the mildly eccentric orbit described in Sec. VI.1 over the full range of $\ell$ plotted; ie. the magnitude of each term scales as $\ell^{0}$, $\ell^{-2}$ and $\ell^{-4}$, after removal of the $A_{(\mu)}$, $B_{(\mu)}$, and $D_{(\mu)}$ terms respectively.

Close to the black hole, on the other hand, the particle moves along a nearly circular trajectory. If the orbit were perfectly circular for all times, ie. $\dot{r}\equiv 0$, then the $(0)$ component would not require regularization at all, and the multipole coefficients would decay exponentially, resulting in a straight line on the semi-logarithmic plot shown in Fig. 17. As the real orbit is not precisely circular, curves eventually deviate from a straight line. Removal of the $A_{(\mu)}$ term is required almost immediately (beginning with $\ell\approx 3$), while the $D_{(\mu)}$ term starts to become important only after $\ell\approx 11$.

This shows that there is a smooth transition from the self-force on a circular orbit, which does not require regularization for the $t$ and $\phi$ components, to that of a generic orbit, for which all components of the self-force require regularization.

Acknowledgements.

We thank Eric Poisson and Eran Rosenthal for useful discussions and suggestions. This work was supported by the Natural Sciences and Engineering Council of Canada.

Appendix A Translation tables

We quote the results of Haas and Poisson (2006) for the translation table between the modes $\Phi_{\ell m}$ and the tetrad components $\Phi_{(\mu)\ell m}$ with respect to the pseudo-Cartesian basis

[TABLE]

and the complex combinations $e^{\alpha}_{(\pm)}:=e^{\alpha}_{(1)}\pm ie^{\alpha}_{(2)}$,

[TABLE]

With these, the spherical-harmonic modes $\Phi_{(\mu)\ell m}(t,r)$ are given in terms of $\Phi_{\ell m}(t,r)$ by

[TABLE]

Appendix B Regularization parameters

For completeness we list the regularization parameters as calculated in Haas and Poisson (2006). Quantities bearing a subscript “[math]” are evaluated at the particle’s position.

[TABLE]

where $f_{0}:=1-2M/r_{0}$ and $\operatorname{sign}(\Delta)$ is equal to $+1$ if $\Delta>0$ and to $-1$ if $\Delta<0$. We have, in addition, $A_{(-)}=\bar{A}_{(+)}$, $A_{(1)}=\operatorname{Re}[A_{(+)}]$, and $A_{(2)}=\operatorname{Im}[A_{(+)}]$.

We also use

[TABLE]

In addition, $B_{(-)}=\bar{B}_{(+)}$, $B_{(1)}=\operatorname{Re}[B_{(+)}]=B^{c}_{(+)}\cos\phi_{0}+B^{s}_{(+)}\sin\phi_{0}$, and $B_{(2)}=\operatorname{Im}[B_{(+)}]=B^{c}_{(+)}\sin\phi_{0}-B^{s}_{(+)}\cos\phi_{0}$.

Here, the rescaled elliptic integrals $\mathcal{E}$ and $\mathcal{K}$ are defined by

[TABLE]

and

[TABLE]

in which $k:=L^{2}/(r_{0}^{2}+L^{2})$.

We also use

[TABLE]

and

[TABLE]

And finally, $D_{(-)}=\bar{D}_{(+)}$, $D_{(1)}=\operatorname{Re}[D_{(+)}]=D^{c}_{(+)}\cos\phi_{0}+D^{s}_{(+)}\sin\phi_{0}$, and $D_{(2)}=\operatorname{Im}[D_{(+)}]=D^{c}_{(+)}\sin\phi_{0}-D^{s}_{(+)}\cos\phi_{0}$.

Appendix C Piecewise polynomials

In two places in the numerical simulation we introduce piecewise polynomials to approximate the scalar field $\psi_{\ell m}$ across the world line, where it is continuous but not differentiable. By a piecewise polynomial we mean a polynomial of the form

[TABLE]

where $u=t-r^{*}$, $v=t+r^{*}$ are characteristic coordinates, $r_{0}^{*}$ is the position of the particle at the time $t(u,v)$, and $N$ is the order of the polynomial, which for our purposes is $N=4$ or less. The two sets of coefficients $c_{nm}$ and $c^{\prime}_{nm}$ are not independent of each other, but are linked via jump conditions that can be derived from the wave equation [Eq. (12)]. To do so, we rewrite the wave equation in the characteristic coordinates $u$ and $v$ and reintroduce the integral over the world line on the right-hand side,

[TABLE]

where $\widehat{S}(\tau)=-8\pi q\frac{\bar{Y}_{\ell m}\bm{(}\pi/2,\phi_{p}(\tau)\bm{)}}{r_{p}(\tau)}$ is the source term and quantities bearing a subscript $p$ are evaluated on the world line at proper time $\tau$.

Here and in the following we use the notation

[TABLE]

to denote the jump in $\partial_{u}^{n}\partial_{v}^{m}\psi$ across the world line. First, we notice that the source term does not contain any derivatives of the Dirac $\delta$-function, causing the solution $\psi$ to be continuous. This means that the zeroth-order jump vanishes: $\left[\psi\right]=0$. Our task is then to find the remaining jump conditions at a point $(t_{0},r_{0}^{*})$ for $n,m\leq 4$. Alternatively, instead of crossing the world line along a line $t=t_{0}=\mathrm{const}$ we can also choose to cross along lines of $u=u_{0}=\mathrm{const}$ or $v=v_{0}=\mathrm{const}$, noting that for a line of constant $v$ the coordinate $u$ runs from $u_{0}+\epsilon$ to $u_{0}-\epsilon$ to cross from the left to the right of the world line. Figure 18 provides a clearer description of the paths taken.

In order to find the jump $\left[\partial_{u}\psi\right]$ we integrate the wave equation along the line $u=u_{0}$ from $v_{0}-\epsilon$ to $v_{0}+\epsilon$

[TABLE]

which, after involving $\int_{v_{0}-\epsilon}^{v_{0}+\epsilon}\delta(v-v_{p})\mathrm{d}v=\theta(v_{p}-v_{0}+\epsilon)\theta(v_{0}-v_{p}+\epsilon)$ and $\delta\bm{(}g(x)\bm{)}=\delta(x-x_{0})/\left\lvert{g^{\prime}(x_{0})}\right\rvert$, yields

[TABLE]

where the overdot denotes differentiation with respect to proper time $\tau$.

Similarly, after first taking a derivative of the wave equation with respect to $v$ and integrating from $u_{0}+\epsilon$ to $u_{0}-\epsilon$, we obtain

[TABLE]

We find

[TABLE]

Systematically repeating this procedure we find expressions for the jumps in all the derivatives that are purely in the $u$ or $v$ direction. Table 3 lists these results.

Jump conditions for derivatives involving both $u$ and $v$ are obtained directly from the wave equation [Eq. (81)]. We see that

[TABLE]

and taking an additional derivative with respect to $u$ on both sides reveals that

[TABLE]

Systematically repeating this procedure we can find jump conditions for each of the mixed derivatives by evaluating

[TABLE]

where $n,m\geq 0$ and derivatives of $V$ with respect to either $u$ or $v$ are evaluated as $\partial_{u}V=-\frac{1}{2}f\partial_{r}V$ and $\partial_{v}V=\frac{1}{2}f\partial_{r}V$, respectively.

The results of Table 3 and Eq. (89) allow us to express the coefficients of the left-hand polynomial in Eq. (80) in terms of the jump conditions and the coefficients of the right-hand side:

[TABLE]

For $N=4$ this leaves us with $25$ unknown coefficients $c_{nm}$ which can be uniquely determined by demanding that the polynomial match the value of the field on the 25 grid points surrounding the particle. When we are interested in integrating the polynomial, as in the case of the potential term in the fourth-order algorithm, we do not need all these terms. Instead, in order to calculate e.g. the integral ${\int\!\!\int}_{\text{cell}}V\psi\,\mathrm{d}u\,\mathrm{d}v$ up to terms of order $h^{5}$, as is needed to achieve overall $O(h^{4})$ convergence, it is sufficient to include only terms such that $n+m\leq 2$, thus reducing the number of unknown coefficients to $6$. In this case Eq. (80) becomes

[TABLE]

The six coefficients can then be determined by matching the polynomial to the field values at the six grid points which lie within the past light cone of the grid point whose field value we want to calculate.